Peturn Normally to move Uniformly.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How should we peturb the individual components of a $d$ dimensional vector in order to ensure that all directions are treated as equally likely. \n",
    "\n",
    "In Random Search algorithms, we need to sample a few surrounding points representatively. Ideally we would like to sample on the unit ball (pick $v$) around the parameter vector $\\theta$ and take a step of size $\\alpha$ in that direction. This would be like a random walk.\n",
    "\n",
    "$$\\theta_{k+1} = \\theta_k + \\alpha v$$\n",
    "\n",
    "In simple random search used in RL, we approximate the gradient by picking $N$ random directions and taking a weighted average of them based on how much reward each of the directions give us.\n",
    "\n",
    "$$\\theta_{k+1} = \\theta_k + \\alpha \\frac{\\sum_{i=1}^N (R(\\theta_k + \\alpha v_i) - R(\\theta_k - \\alpha v_i))v_i}{N}$$\n",
    "\n",
    "Now the question is how to sample for $v$. In the lower dimensions we can sample from the unit circle using polar coordinates, it should be possible "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# returns a random d dimensional vector, a direction to peturb in \n",
    "def direction(d,t):\n",
    "    # if type == uniform\n",
    "    if(t == 'u'):\n",
    "        return np.random.uniform(-2/np.sqrt(d), 2/np.sqrt(d), d)\n",
    "    elif(t == 'n'):\n",
    "        return np.random.normal(0, 1/np.sqrt(d), d)\n",
    "    elif(t == 's'):\n",
    "        # a point on the N-Sphere r = 1 so it is ommited\n",
    "        angles = np.random.uniform(0, np.pi, d-2)\n",
    "        angleLast = np.random.uniform(0, 2*np.pi,1)[0]\n",
    "        x = np.zeros(d)\n",
    "        x[0] = np.cos(angles[0])\n",
    "        for i in range(1,d-1):\n",
    "            temp = 1\n",
    "            for j in range(i):\n",
    "                temp = temp * np.sin(angles[j])\n",
    "            if(i == d-2):\n",
    "                x[i] = temp * np.cos(angleLast)\n",
    "            else:\n",
    "                x[i] = temp*np.cos(angles[i])\n",
    "        x[d-1] = x[d-2]*np.tan(angleLast)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#N = 10000\n",
    "N = 100 # number of directions sampled AND the number of dimensions.\n",
    "d = N\n",
    "hN = []\n",
    "normal = direction(d,'n').reshape(d,1)\n",
    "\n",
    "for i in range(N-1):\n",
    "    hN.append(np.linalg.norm(direction(d,'n')))\n",
    "    normal = np.concatenate((normal,direction(d,'n').reshape(d,1)), axis = 1)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.hist(hN)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "hU = []\n",
    "\n",
    "uniform = direction(d,'u').reshape(d,1)\n",
    "\n",
    "for i in range(N-1):\n",
    "    hU.append(np.linalg.norm(direction(d,'u')))\n",
    "    uniform = np.concatenate((uniform,direction(d,'u').reshape(d,1)), axis = 1)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAARh0lEQVR4nO3df4xlZX3H8fengDZVWhd3QATW1ZZQ0QiS6aqlNaBlBSRiG9NCjG4tZtVoookm3baJNJomNI2aKFbcykZsFE2rKAk/N2iDRkAGwo9FVJCsYV3Crq4FjaZ29ds/5mxyHe6dvXPP3Z2Zp+9XcnPPeZ7n3Od5cmY/c/bMOeemqpAkteu3lnsAkqRDy6CXpMYZ9JLUOINekhpn0EtS445c7gEMs3bt2lq/fv1yD0OSVo277rrrR1U1M6xuRQb9+vXrmZubW+5hSNKqkeQHo+o8dSNJjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIad9CgT3JSkq8leTDJA0ne3ZUfk2R7koe69zUjtt/UtXkoyaZpT0CStLhxjuj3A++tqhcCLwfemeRUYAtwS1WdDNzSrf+GJMcAlwIvAzYAl476hSBJOjQOGvRV9VhV3d0t/xR4EDgBuBC4qmt2FfD6IZu/BtheVfuq6ifAduDcaQxckjSeJd0Zm2Q98FLgDuC4qnoM5n8ZJDl2yCYnAI8OrO/qyoZ99mZgM8C6deuWMizpsFm/5bpl6XfnZa9dln7VhrH/GJvkmcAXgfdU1ZPjbjakbOhXWlXV1qqararZmZmhj2uQJE1grKBPchTzIf/ZqvpSV/x4kuO7+uOBPUM23QWcNLB+IrB78uFKkpZqnKtuAlwJPFhVHx6ouhY4cBXNJuArQza/CdiYZE33R9iNXZkk6TAZ54j+TOBNwKuS3NO9zgcuA85J8hBwTrdOktkknwKoqn3AB4E7u9cHujJJ0mFy0D/GVtU3GH6uHeDVQ9rPAW8dWN8GbJt0gJKkfrwzVpIaZ9BLUuMMeklqnEEvSY1bkd8ZKx3Mct2hulyWc77elbv6eUQvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ17qDPukmyDbgA2FNVL+7KvgCc0jV5FvDfVXX6kG13Aj8FfgXsr6rZKY1bkjSmcR5q9mngcuAzBwqq6q8OLCf5EPDEItufXVU/mnSAkqR+xvkqwVuTrB9W131x+F8Cr5rusCRJ09L3HP2fAo9X1UMj6gu4OcldSTb37EuSNIG+z6O/GLh6kfozq2p3kmOB7Um+U1W3DmvY/SLYDLBu3bqew5IkHTDxEX2SI4G/AL4wqk1V7e7e9wDXABsWabu1qmaranZmZmbSYUmSFuhz6ubPgO9U1a5hlUmekeToA8vARmBHj/4kSRM4aNAnuRq4DTglya4kl3RVF7HgtE2S5ya5vls9DvhGknuBbwHXVdWN0xu6JGkc41x1c/GI8r8eUrYbOL9bfgQ4ref4JEk9eWesJDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNW6crxLclmRPkh0DZf+Y5IdJ7ule54/Y9twk303ycJIt0xy4JGk84xzRfxo4d0j5R6rq9O51/cLKJEcAHwfOA04FLk5yap/BSpKW7qBBX1W3Avsm+OwNwMNV9UhV/RL4PHDhBJ8jSeqhzzn6dyW5rzu1s2ZI/QnAowPru7qyoZJsTjKXZG7v3r09hiVJGjRp0H8C+H3gdOAx4END2mRIWY36wKraWlWzVTU7MzMz4bAkSQtNFPRV9XhV/aqqfg38G/OnaRbaBZw0sH4isHuS/iRJk5so6JMcP7D658COIc3uBE5O8vwkTwMuAq6dpD9J0uSOPFiDJFcDZwFrk+wCLgXOSnI686didgJv69o+F/hUVZ1fVfuTvAu4CTgC2FZVDxySWUiSRjpo0FfVxUOKrxzRdjdw/sD69cBTLr2UJB0+3hkrSY0z6CWpcQa9JDXOoJekxhn0ktS4g151I42yfst1yz0ESWPwiF6SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWrcQYM+ybYke5LsGCj7lyTfSXJfkmuSPGvEtjuT3J/kniRz0xy4JGk84xzRfxo4d0HZduDFVfUS4HvA3y2y/dlVdXpVzU42RElSHwcN+qq6Fdi3oOzmqtrfrd4OnHgIxiZJmoJpnKP/G+CGEXUF3JzkriSbF/uQJJuTzCWZ27t37xSGJUmCnkGf5B+A/cBnRzQ5s6rOAM4D3pnklaM+q6q2VtVsVc3OzMz0GZYkacDEQZ9kE3AB8MaqqmFtqmp3974HuAbYMGl/kqTJTBT0Sc4F/hZ4XVX9fESbZyQ5+sAysBHYMaytJOnQGefyyquB24BTkuxKcglwOXA0sL27dPKKru1zk1zfbXoc8I0k9wLfAq6rqhsPySwkSSMd9Dtjq+riIcVXjmi7Gzi/W34EOK3X6CRJvXlnrCQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDVurKBPsi3JniQ7BsqOSbI9yUPd+5oR227q2jzUfaG4JOkwGveI/tPAuQvKtgC3VNXJwC3d+m9IcgxwKfAyYANw6ahfCJKkQ2OsoK+qW4F9C4ovBK7qlq8CXj9k09cA26tqX1X9BNjOU39hSJIOoT7n6I+rqscAuvdjh7Q5AXh0YH1XV/YUSTYnmUsyt3fv3h7DkiQNOtR/jM2QshrWsKq2VtVsVc3OzMwc4mFJ0v8ffYL+8STHA3Tve4a02QWcNLB+IrC7R5+SpCXqE/TXAgeuotkEfGVIm5uAjUnWdH+E3diVSZIOk3Evr7wauA04JcmuJJcAlwHnJHkIOKdbJ8lskk8BVNU+4IPAnd3rA12ZJOkwOXKcRlV18YiqVw9pOwe8dWB9G7BtotFJknrzzlhJapxBL0mNM+glqXEGvSQ1zqCXpMaNddWNVrb1W65b7iFIWsE8opekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMZ5Z6ykRS3Xndc7L3vtsvTbIo/oJalxEwd9klOS3DPwejLJexa0OSvJEwNt3t9/yJKkpZj41E1VfRc4HSDJEcAPgWuGNP16VV0waT+SpH6mderm1cD3q+oHU/o8SdKUTCvoLwKuHlH3iiT3JrkhyYtGfUCSzUnmkszt3bt3SsOSJPUO+iRPA14H/MeQ6ruB51XVacDHgC+P+pyq2lpVs1U1OzMz03dYkqTONI7ozwPurqrHF1ZU1ZNV9bNu+XrgqCRrp9CnJGlM0wj6ixlx2ibJc5KkW97Q9ffjKfQpSRpTrxumkvwOcA7wtoGytwNU1RXAG4B3JNkP/AK4qKqqT5+SpKXpFfRV9XPg2QvKrhhYvhy4vE8fkqR+vDNWkhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGtc76JPsTHJ/knuSzA2pT5KPJnk4yX1JzujbpyRpfL2+SnDA2VX1oxF15wEnd6+XAZ/o3iVJh8HhOHVzIfCZmnc78Kwkxx+GfiVJTOeIvoCbkxTwyarauqD+BODRgfVdXdljg42SbAY2A6xbt24Kw5K0mq3fct2y9b3zstcuW9+HwjSO6M+sqjOYP0XzziSvXFCfIdvUUwqqtlbVbFXNzszMTGFYkiSYQtBX1e7ufQ9wDbBhQZNdwEkD6ycCu/v2K0kaT6+gT/KMJEcfWAY2AjsWNLsWeHN39c3LgSeq6jEkSYdF33P0xwHXJDnwWZ+rqhuTvB2gqq4ArgfOBx4Gfg68pWefkqQl6BX0VfUIcNqQ8isGlgt4Z59+JEmT885YSWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJatzEQZ/kpCRfS/JgkgeSvHtIm7OSPJHknu71/n7DlSQtVZ+vEtwPvLeq7u6+IPyuJNur6tsL2n29qi7o0Y8kqYeJj+ir6rGqurtb/inwIHDCtAYmSZqOqZyjT7IeeClwx5DqVyS5N8kNSV60yGdsTjKXZG7v3r3TGJYkiSkEfZJnAl8E3lNVTy6ovht4XlWdBnwM+PKoz6mqrVU1W1WzMzMzfYclSer0CvokRzEf8p+tqi8trK+qJ6vqZ93y9cBRSdb26VOStDR9rroJcCXwYFV9eESb53TtSLKh6+/Hk/YpSVq6PlfdnAm8Cbg/yT1d2d8D6wCq6grgDcA7kuwHfgFcVFXVo09J0hJlJebu7Oxszc3NLfcwlmT9luuWewiSVrmdl7124m2T3FVVs8PqvDNWkhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIa1+cRCCuSd6hK0m/yiF6SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMb1/XLwc5N8N8nDSbYMqX96ki909XckWd+nP0nS0vX5cvAjgI8D5wGnAhcnOXVBs0uAn1TVHwAfAf550v4kSZPpc0S/AXi4qh6pql8CnwcuXNDmQuCqbvk/gVcnSY8+JUlL1OcRCCcAjw6s7wJeNqpNVe1P8gTwbOBHCz8syWZgc7f6syTfXaTvtcM+oxGtzq3VeUG7c2t1XrBC55Z+5zyeN6qiT9APOzKvCdrMF1ZtBbaO1XEyN+rbzle7VufW6ryg3bm1Oi9oe27D9Dl1sws4aWD9RGD3qDZJjgR+D9jXo09J0hL1Cfo7gZOTPD/J04CLgGsXtLkW2NQtvwH4alUNPaKXJB0aE5+66c65vwu4CTgC2FZVDyT5ADBXVdcCVwL/nuRh5o/kL5rGoBnzFM8q1ercWp0XtDu3VucFbc/tKeIBtiS1zTtjJalxBr0kNW5FBX2SbUn2JNkxoj5JPto9UuG+JGcM1G1K8lD32jRs++XUc26/SnJP91r4B+9lNca8/jDJbUn+J8n7FtQt+giN5dZzbjuT3N/ts7nDM+LxjDGvN3Y/g/cl+WaS0wbqVvs+W2xuK3af9VZVK+YFvBI4A9gxov584Abmr89/OXBHV34M8Ej3vqZbXrPc85nG3Lq6ny33+HvM61jgj4B/At43UH4E8H3gBcDTgHuBU5d7PtOYW1e3E1i73HOYcF5/fODfD/OPODnw76yFfTZ0bit9n/V9ragj+qq6lcWvs78Q+EzNux14VpLjgdcA26tqX1X9BNgOnHvoRzy+HnNb0Q42r6raU1V3Av+7oGqcR2gsqx5zW9HGmNc3u39HALczf48MtLHPRs2taSsq6Mcw7LELJyxSvposNoffTjKX5PYkrz/8QzskWthniyng5iR3dY/3WK0uYf5/mtDePhucG7Szz56izyMQlsOoRyqM/aiFFWyxOayrqt1JXgB8Ncn9VfX9wzi2Q6GFfbaYM7t9diywPcl3uqPNVSPJ2cyH4Z8cKBrSbFXusyFzgwb22Sir7Yh+1GMXxnkcw0o3cg5VdeD9EeC/gJce7sEdAi3ss5EG9tke4BrmT3usGkleAnwKuLCqftwVN7HPRsxt1e+zxay2oL8WeHN3hcrLgSeq6jHm787dmGRNkjXAxq5sNRk6t25OTwdIshY4E/j2cg50SsZ5hMaqlOQZSY4+sMz8z+PQq0BWoiTrgC8Bb6qq7w1Urfp9Nmpuq32fHcyKOnWT5GrgLGBtkl3ApcBRAFV1BXA981enPAz8HHhLV7cvyQeZ/0EE+EBVraiHp006N+CFwCeT/Jr5X8yXVdWKCfqDzSvJc4A54HeBXyd5D/NXajyZIY/QWI45jDLp3Jh/BO41mf/qhSOBz1XVjYd/BsON8bP4fuYfJ/6v3Rz2V9VsjXjsyTJMYaRJ5wYcxwreZ335CARJatxqO3UjSVoig16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ17v8AaQ0RyTRVuaoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.hist(hU)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.220446049250313e-16\n",
      "2.220446049250313e-16\n",
      "2.220446049250313e-16\n",
      "2.220446049250313e-16\n",
      "2.220446049250313e-16\n"
     ]
    }
   ],
   "source": [
    "N = 1000\n",
    "hS = []\n",
    "\n",
    "spherical = direction(d,'s').reshape(d,1)\n",
    "\n",
    "for i in range(N-1):\n",
    "    hS.append(np.linalg.norm(direction(d,'s')))\n",
    "    spherical = np.concatenate((spherical,direction(d,'s').reshape(d,1)), axis = 1)\n",
    "    \n",
    "    \n",
    "for i in hS:\n",
    "    # All vectors are close enough to 1 in length. Even for small d\n",
    "    if((i-1)>10**-16):\n",
    "        print(i-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So from the histograms above we can see all these methods give us points on the unit sphere. (Uniform gives us almost) But are they all uncorrelated? Let us see,\n",
    "\n",
    "See how as $N$ increases the matrix tends to $I$ showing that they are indeed drawn i.i.d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.92668584, -0.15721723,  0.06571708, ...,  0.01862396,\n",
       "         0.10557481,  0.00388946],\n",
       "       [-0.15721723,  1.24468988,  0.1562552 , ...,  0.02279731,\n",
       "         0.05295795, -0.10228298],\n",
       "       [ 0.06571708,  0.1562552 ,  1.17864114, ...,  0.17327958,\n",
       "        -0.07316159, -0.04910568],\n",
       "       ...,\n",
       "       [ 0.01862396,  0.02279731,  0.17327958, ...,  0.91062053,\n",
       "         0.09141464, -0.11415045],\n",
       "       [ 0.10557481,  0.05295795, -0.07316159, ...,  0.09141464,\n",
       "         1.20160252,  0.01323173],\n",
       "       [ 0.00388946, -0.10228298, -0.04910568, ..., -0.11415045,\n",
       "         0.01323173,  1.11486184]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.matmul(normal.T, normal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.2516747 ,  0.10667293,  0.18255579, ...,  0.01486249,\n",
       "         0.1055461 ,  0.19918725],\n",
       "       [ 0.10667293,  1.35930531,  0.07551817, ..., -0.08116834,\n",
       "        -0.05050703, -0.10254823],\n",
       "       [ 0.18255579,  0.07551817,  1.3725409 , ..., -0.01580158,\n",
       "        -0.28052221, -0.12417085],\n",
       "       ...,\n",
       "       [ 0.01486249, -0.08116834, -0.01580158, ...,  1.79472098,\n",
       "        -0.2101487 ,  0.01505737],\n",
       "       [ 0.1055461 , -0.05050703, -0.28052221, ..., -0.2101487 ,\n",
       "         1.32117951,  0.03939992],\n",
       "       [ 0.19918725, -0.10254823, -0.12417085, ...,  0.01505737,\n",
       "         0.03939992,  1.39414934]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.matmul(uniform.T, uniform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        , -0.71569501, -0.86257711, ...,  0.3974408 ,\n",
       "         0.5421008 , -0.65220909],\n",
       "       [-0.71569501,  1.        ,  0.72344861, ..., -0.20731027,\n",
       "        -0.77324699,  0.8634745 ],\n",
       "       [-0.86257711,  0.72344861,  1.        , ..., -0.5127261 ,\n",
       "        -0.32940257,  0.75837026],\n",
       "       ...,\n",
       "       [ 0.3974408 , -0.20731027, -0.5127261 , ...,  1.        ,\n",
       "        -0.27848992, -0.55480492],\n",
       "       [ 0.5421008 , -0.77324699, -0.32940257, ..., -0.27848992,\n",
       "         1.        , -0.41714964],\n",
       "       [-0.65220909,  0.8634745 ,  0.75837026, ..., -0.55480492,\n",
       "        -0.41714964,  1.        ]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.matmul(spherical.T, spherical)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets see what happens to this vector if we rotate it, lets generate a arbitary orthogonal matrix and apply it and see if it still lies on sphere....\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}