.ipynb_checkpoints/original_paper_implementation-checkpoint.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "51d6cfbf-a61a-4fd9-a73d-492c0b43dd3b",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Legacy Code from https://arxiv.org/pdf/1806.07572.pdf section 4.1 page 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5cd8e4cc-5537-4b66-9b37-f20a1e6a1def",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8e7f4c76-71d3-44ba-81e8-98fa4e16807a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "1d input x_i and x_j, alias to x and x' in the paper\n",
    "for 2d input we should use np.meshgrid or np.mgrid\n",
    "\"\"\"\n",
    "\n",
    "num_input = 100\n",
    "\n",
    "x_i = np.linspace(0.0, 1.0, num=num_input)\n",
    "x_j = np.linspace(0.0, 1.0, num=num_input)\n",
    "# x = np.random.uniform(size=num_input)\n",
    "\n",
    "# take outer product of x_i and x_j to get 100 x 100 matrix\n",
    "# in the paper sigma was defined entry-wise, but we want to compute the matrix at once\n",
    "n_0 = 1\n",
    "beta = 0.01\n",
    "\n",
    "def calc_sigma_1(x_i, x_j):\n",
    "    return np.outer(x_i, x_j) / n_0 + beta**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "84deb069-acd6-449a-a812-f4892abf6ad0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "f(x_i) and f(x_j) follow N(0, cov) where cov is a contatenation of 4 smaller matrices\n",
    "\"\"\"\n",
    "cov_ii = calc_sigma_1(x_i, x_i)\n",
    "cov_ij = calc_sigma_1(x_i, x_j)\n",
    "cov_ji = calc_sigma_1(x_j, x_i)\n",
    "cov_jj = calc_sigma_1(x_j, x_j)\n",
    "\n",
    "# concatenate the 4 sigmas\n",
    "v1 = np.vstack((cov_ii, cov_ji))\n",
    "v2 = np.vstack((cov_ij, cov_jj))\n",
    "\n",
    "cov = np.hstack((v1, v2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "aaac0ad8-2001-4da4-b9a0-441254a62fab",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\"\"\"\n",
    "activation functions and their derivatives\n",
    "\"\"\"\n",
    "\n",
    "def relu_kth(x, k):\n",
    "    return np.maximum(x**k, 0)\n",
    "\n",
    "def relu(x):\n",
    "    return np.maximum(x, 0)\n",
    "\n",
    "def d_relu(x):\n",
    "    return 1.0 * (x > 0)\n",
    "\n",
    "def d_relu_kth(x, k):\n",
    "    return k * x**(k-1) * (x > 0)\n",
    "\n",
    "# def is_pos_def(x):\n",
    "#     return np.all(np.linalg.eigvals(x) > 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b870448c-89d1-4f6a-ad08-ec852c91acb8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "gaol: approximate expectation using a finite sum and take the mean\n",
    "now sample Y_i, Y_j from this distribution N(0, cov)\n",
    "sample and compute outer product first, then take mean\n",
    "\"\"\"\n",
    "def calc_sigma_2(activation):\n",
    "    mean = [0] * 2 * num_input\n",
    "    num_samples = 1000\n",
    "    sum_N = np.zeros((num_input, num_input))\n",
    "    for i in range(num_samples):\n",
    "        sample = np.random.multivariate_normal(mean, cov, size=1)[0]\n",
    "        f_xi = sample[:num_input]\n",
    "        prod = np.outer(activation(f_xi), activation(f_xi))\n",
    "        sum_N += np.outer(activation(f_xi), activation(f_xi))\n",
    "\n",
    "    sigma_2 = (sum_N / num_samples) + beta**2\n",
    "    \n",
    "    return sigma_2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c97228ae-ae3c-41d2-b1e9-000907dd767c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1800beac5b0>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAWb0lEQVR4nO3df7Cc1X3f8fd3f0gICYEEFyokGUmuaiO3CcaKosaZjifEraRmIqfTZGAmRvV4RmEKDclkpiXNH03/6VC3sWumFCrbiiHJmHgSJ2gymlAPtaFpioNwXIr4KQQYWQIuyBYgYd1f3/6xz7137969965+Lrrn/ZrZ2d3nOWf3HH7s557nPOd5IjORJJWp1u8GSJL6xxCQpIIZApJUMENAkgpmCEhSwRr9bsCpuOKKK3LNmjX9boYkXVCeeOKJNzNzoNu+CyoE1qxZw759+/rdDEm6oETEKzPt83CQJBXMEJCkghkCklQwQ0CSCmYISFLBDAFJKpghIEkFKyIEHn7mde759ov9boYkve8UEQLffm6QXY8aApLUqYgQaNSDkVFvniNJnYoIgWa9xvDYWL+bIUnvO4WEQDDsSECSpikiBBq1GqNjifdTlqSpigiBZj0AHA1IUociQqBRb3VzxHkBSZqijBCoORKQpG6KCIHm+Ehg1JGAJLUrIgQa1ZzAyJgjAUlqV0QINGutbg6NOBKQpHZlhEDDkYAkdVNECDRqzglIUjdFhIDrBCSpu55CICK2RMRzEXEgIu7osj8i4q5q/5MRcX21fXVEfCsinomI/RFxe1ud5RHxzYh4oXpedva6NdXESMB1ApI0xZwhEBF14G5gK7ABuCkiNnQU2wqsrx47gXuq7SPAb2XmtcBm4Na2uncAD2fmeuDh6v050XAkIEld9TIS2AQcyMyDmTkEPABs7yizHbg/Wx4DLouIFZl5JDO/C5CZ7wDPACvb6txXvb4P+NSZdWVmrhOQpO56CYGVwKtt7w8x+UPec5mIWAN8FPhOtemqzDwCUD1f2e3LI2JnROyLiH2Dg4M9NHe68RXDnh0kSVP1EgLRZVvnr+msZSJiCfCnwG9k5tu9Nw8yc1dmbszMjQMDA6dSdUKzUa0TcCQgSVP0EgKHgNVt71cBh3stExFNWgHwR5n5jbYyr0fEiqrMCuCNU2t675oTp4g6EpCkdr2EwOPA+ohYGxELgBuBPR1l9gA3V2cJbQaOZeaRiAjgK8Azmfn5LnV2VK93AA+edi/mMHHZCEcCkjRFY64CmTkSEbcBDwF1YHdm7o+IW6r99wJ7gW3AAeAE8Jmq+seBTwP/LyK+V237t5m5F7gT+HpEfBb4PvDLZ61XHSbWCTgnIElTzBkCANWP9t6Obfe2vU7g1i71/oru8wVk5lvADafS2NPlimFJ6q6IFcOTh4McCUhSuyJCYHydwLArhiVpiiJCYOLOYl5KWpKmKCIExtcJuFhMkqYqIwSqiWGvHSRJUxURAq4TkKTuygiBmusEJKmbIkIgImjUwpGAJHUoIgSgdUjIiWFJmqqYEGjWagw7EpCkKYoJgUY9DAFJ6lBMCDTrNS8bIUkdigoB1wlI0lTFhEBrYtjDQZLUrpwQqIWHgySpQzEh0Doc5EhAktoVEwKuE5Ck6coJAdcJSNI0xYTAAg8HSdI0xYRAo+7EsCR1KigEal5FVJI6FBMCTa8iKknTFBMCHg6SpOkKCoEaw64YlqQpigmBpiuGJWmaYkKg4SmikjRNMSHgVUQlabqCQsCriEpSp2JCoFHzpjKS1KmYEGh6e0lJmqaYEPAqopI0XTkhUKsxOpaMGQSSNKGYEGjWA8AFY5LUpqAQaHXVyWFJmlRMCDQMAUmappgQ8HCQJE3XUwhExJaIeC4iDkTEHV32R0TcVe1/MiKub9u3OyLeiIinOur8bkT8ICK+Vz22nXl3ZtaoORKQpE5zhkBE1IG7ga3ABuCmiNjQUWwrsL567ATuadv3VWDLDB//hcy8rnrsPcW2n5LG+EjAtQKSNKGXkcAm4EBmHszMIeABYHtHme3A/dnyGHBZRKwAyMxHgaNns9GnY/xwkGsFJGlSLyGwEni17f2hatuplunmturw0e6IWNatQETsjIh9EbFvcHCwh4/sbvJwkCMBSRrXSwhEl22df073UqbTPcAHgeuAI8DvdSuUmbsyc2NmbhwYGJjjI2c2PhIYMgQkaUIvIXAIWN32fhVw+DTKTJGZr2fmaGaOAV+iddjpnHGdgCRN10sIPA6sj4i1EbEAuBHY01FmD3BzdZbQZuBYZh6Z7UPH5wwqvwQ8NVPZs2FinYCniErShMZcBTJzJCJuAx4C6sDuzNwfEbdU++8F9gLbgAPACeAz4/Uj4mvAJ4ArIuIQ8O8y8yvA5yLiOlqHjV4Gfu3sdWu6Zm387CBHApI0bs4QAKhO39zbse3ettcJ3DpD3Ztm2P7p3pt55lwxLEnTFbNiuOGKYUmappgQaLpiWJKmKSYEXDEsSdMVEwJNQ0CSpikoBDwcJEmdigkB1wlI0nTFhIDrBCRpumJCYHKdgCMBSRpXUAh4KWlJ6lRMCIyvE/BwkCRNKiYEXCcgSdOVEwLVxLBzApI0qZgQiAia9WDYOQFJmlBMCEDrFpOOBCRpUlkhUA8nhiWpTVEh0KzXXDEsSW2KCoFGLbx2kCS1KSoEmvUaQ84JSNKEokKgUXckIEntigoB5wQkaaqiQqBR8+wgSWpXVAg0664TkKR2RYVAox5eRVSS2hQVAs1azQvISVKbokLAs4MkaarCQsCRgCS1KyoEFnjtIEmaoqgQaNRcJyBJ7coKAecEJGmKokKgWa8x7EhAkiYUFQJeRVSSpiorBOo1J4YlqU1RIdCsh6eISlKbokLAewxL0lRFhUCzEQx77SBJmlBWCDgSkKQpegqBiNgSEc9FxIGIuKPL/oiIu6r9T0bE9W37dkfEGxHxVEed5RHxzYh4oXpedubdmV2jHowljDkakCSghxCIiDpwN7AV2ADcFBEbOoptBdZXj53APW37vgps6fLRdwAPZ+Z64OHq/TnVrLe661oBSWrpZSSwCTiQmQczcwh4ANjeUWY7cH+2PAZcFhErADLzUeBol8/dDtxXvb4P+NRptP+UNGoB4FoBSar0EgIrgVfb3h+qtp1qmU5XZeYRgOr5ym6FImJnROyLiH2Dg4M9NHdmjWokYAhIUksvIRBdtnX+ivZS5rRk5q7M3JiZGwcGBs7os5r1VjOHnByWJKC3EDgErG57vwo4fBplOr0+fsioen6jh7ackUatGgk4JyBJQG8h8DiwPiLWRsQC4EZgT0eZPcDN1VlCm4Fj44d6ZrEH2FG93gE8eArtPi3jIwEPB0lSy5whkJkjwG3AQ8AzwNczc39E3BIRt1TF9gIHgQPAl4B/OV4/Ir4G/B/gQxFxKCI+W+26E/hkRLwAfLJ6f05NnB3k4SBJAqDRS6HM3Evrh759271trxO4dYa6N82w/S3ghp5behY0xkcCrhOQJKCwFcPjcwKOBCSppagQcE5AkqYqKgQazglI0hRFhUCzWjHsjWUkqaWsEGi4TkCS2hUVAl47SJKmKioEXCcgSVMVFQKuE5CkqcoKAdcJSNIURYWA6wQkaaqiQsB1ApI0VVEhMLFOwDkBSQJKC4GJO4s5EpAkKCwEGs4JSNIURYXAxDoBVwxLElBYCLhiWJKmKioE6hMh4EhAkqCwEIgImvVgyJGAJAGFhQC0Vg07EpCkluJCoFkPrx0kSZUCQ6DmimFJqhQXAo16eHaQJFXKC4FazXUCklQpLgSajgQkaUJxIdCo17zHsCRVyguBWjA04khAkqDAEGg6EpCkCQWGgHMCkjSuuBBouE5AkiYUFwKuGJakScWFgNcOkqRJxYVAsx4MOycgSUCBIdCoOScgSePKCwHnBCRpQnEh4FVEJWlSgSHgOgFJGtdTCETEloh4LiIORMQdXfZHRNxV7X8yIq6fq25E/G5E/CAivlc9tp2dLs3OawdJ0qQ5QyAi6sDdwFZgA3BTRGzoKLYVWF89dgL39Fj3C5l5XfXYe6ad6UWz5tlBkjSul5HAJuBAZh7MzCHgAWB7R5ntwP3Z8hhwWUSs6LHuedWou05Aksb1EgIrgVfb3h+qtvVSZq66t1WHj3ZHxLKeW30GGvVg2LODJAnoLQSiy7bOX9GZysxW9x7gg8B1wBHg97p+ecTOiNgXEfsGBwd7aO7smq4TkKQJvYTAIWB12/tVwOEey8xYNzNfz8zRzBwDvkTr0NE0mbkrMzdm5saBgYEemju7Rj3IhFFHA5LUUwg8DqyPiLURsQC4EdjTUWYPcHN1ltBm4FhmHpmtbjVnMO6XgKfOsC89uXRRE4C33j15Pr5Okt7XGnMVyMyRiLgNeAioA7szc39E3FLtvxfYC2wDDgAngM/MVrf66M9FxHW0Dg+9DPzaWezXjK5dsRSAp4+8zZVLLzofXylJ71tzhgBAdfrm3o5t97a9TuDWXutW2z99Si09S8ZDYP/ht/nEh67sRxMk6X2juBXDly5qsmrZIp4+8na/myJJfVdcCAB85OqlPHPYEJCkIkNgw4pLeemt4xw/OdLvpkhSX5UZAlcvJROefe2dfjdFkvqq2BAAePrwsT63RJL6q8gQuPrSi7js4qaTw5KKV2QIRAQbVizlaSeHJRWuyBAA2LBiKc++9o5XFJVUtGJD4CMrl3JyZIyDbx7vd1MkqW+KDYENKy4F8JCQpKIVGwLrBhazoFFzclhS0YoNgWa9xoeuusSRgKSiFRsC0Lp8xNNH3qZ1/TtJKk/RIbDh6qUcPT7E4WM/7ndTJKkvig6Bj13Tuq3x4y8d7XNLJKk/ig6Ba//OUi5d1OSxg2/1uymS1BdFh0CtFmxau9wQkFSsokMAYPO6y3n5rRMcOfZev5siSeedIbBuOQDfOei8gKTyFB8CzgtIKlnxIVCrBT+1xnkBSWUqPgSgdUjIeQFJJTIEaE0Og/MCkspjCADXrljK0osaHhKSVBxDAKjXgk1rLzcEJBXHEKiMzwsc/pHzApLKYQhUfu7DVxIBf/jYK/1uiiSdN4ZAZd3AEn7hJ67mq3/9Mm+9e7LfzZGk88IQaHP7Dev58fAou/7XwX43RZLOC0Ogzd+9cgm/+JNXc/9fv8KbjgYkFcAQ6PDrN6zn5Mgo//2RF/vdFEk65wyBDusGlvCpj67kDx57hUM/PNHv5kjSOWUIdHH7DeupRbD9v/5vHn1+sN/NkaRzxhDo4prLF/PgrR/n8iUL2PH7f8N/euhZRkbH+t0sSTrrDIEZrL/qEh689Wf5lY+t5u5vvcjPf/4R/uxvDzE6lv1umiSdNYbALBYtqPMf//lP8OWbN3JRs85v/vH/5R9/4RG+8lcv8eLgu2QaCJIubHEh/ZBt3Lgx9+3b15fvHhtL/nL/a9z18As8+9o7AKxatoif+eDlfOyaZXzsmmWsu2IJtVr0pX2SNJOIeCIzN3bd10sIRMQW4ItAHfhyZt7ZsT+q/duAE8C/yMzvzlY3IpYDfwysAV4GfiUzfzhbO/oZAu1ePXqCR54f5JHnB3n85aP86MQwAEsWNvh7Vy3hwyuWsvbyxVx6cZNLFzW5ZGGDRr1GvRYsatZZvXwRl1zU7HMvJJXijEIgIurA88AngUPA48BNmfl0W5ltwL+iFQI/DXwxM396troR8TngaGbeGRF3AMsy89/M1pb3Swi0y0xeevM4+175Ift/cIxnX3uHZ197h2PvDc9a7/LFC1i1/GIuWdhg0YI6ixfUuXhhg8UL6ixe2GBJ9Vi8sMGCRo16BPV6sLBeY2GzxsJGnQWNGs16jQWNGo1aEAABtYjqAVE9T2yrQaNWm9gnaf6bLQQaPdTfBBzIzIPVhz0AbAeebiuzHbg/W4nyWERcFhEraP2VP1Pd7cAnqvr3Ad8GZg2B96OIYN3AEtYNLIGNq4FWMLxzcoRjJ4Y59t4w754cYXQsGRlLjp8c4ftHT/DKW8c59MP3OH5yhDffPcmJoVFODI1w/OQo7w2Pnpe216IKhBoEk6ExJRpi8ikiiJh8PRk0rfrt+1r/bFrhU6/FtNAxfqRT8x/+2T/gp9YsP+uf20sIrARebXt/iNZf+3OVWTlH3asy8whAZh6JiCu7fXlE7AR2AnzgAx/oobn9FxEsvajJ0ouarD6N+qNjyfGhEd798QjHT44wNDo2ESLDI2MMjY7x4+ExhkbGGBodZWhkjJGxJBOSVghlwljmxNlMrdeT20bGkrGxZLR6P1lnsh1J6834YDFzfEvrc8aytW1srFV2/PtpLzOWjGZrTqXzcyX1blGzfk4+t5cQ6PZHW+f/xTOV6aXurDJzF7ALWoeDTqXuhapemwwRSTqXejlF9BBM+YN2FXC4xzKz1X29OmRE9fxG782WJJ0NvYTA48D6iFgbEQuAG4E9HWX2ADdHy2bgWHWoZ7a6e4Ad1esdwINn2BdJ0ima83BQZo5ExG3AQ7RO89ydmfsj4pZq/73AXlpnBh2gdYroZ2arW330ncDXI+KzwPeBXz6rPZMkzcnFYpI0z812iqiXjZCkghkCklQwQ0CSCmYISFLBLqiJ4YgYBF45zepXAG+exeZcKOx3eUrtu/2e2TWZOdBtxwUVAmciIvbNNDs+n9nv8pTad/t9ejwcJEkFMwQkqWAlhcCufjegT+x3eUrtu/0+DcXMCUiSpitpJCBJ6mAISFLBigiBiNgSEc9FxIHqfsbzUkSsjohvRcQzEbE/Im6vti+PiG9GxAvV87J+t/Vsi4h6RPxtRPxF9X7e9xmgupXrn0TEs9W/939YQt8j4jer/8afioivRcRF87HfEbE7It6IiKfats3Yz4j47ep37rmI+Ce9fMe8D4HqZvd3A1uBDcBNEbGhv606Z0aA38rMa4HNwK1VX+8AHs7M9cDD1fv55nbgmbb3JfQZ4IvAX2bmh4GfpPXPYF73PSJWAr8ObMzMv0/rMvU3Mj/7/VVgS8e2rv2s/l+/EfhIVee/Vb9/s5r3IQBsorrZfWYOAeM3u593MvNIZn63ev0OrR+ElbT6e19V7D7gU31p4DkSEauAfwp8uW3zvO4zQEQsBf4R8BWAzBzKzB9RQN9p3QtlUUQ0gItp3bFw3vU7Mx8FjnZsnqmf24EHMvNkZr5E6/4um+b6jhJCoNvN7lf2qS3nTUSsAT4KfAe4qrrTG9XzlX1s2rnwX4B/DYy1bZvvfQZYBwwCv18dCvtyRCxmnvc9M38A/GdaN6M6QutOhv+Ded7vNjP187R+60oIgTO+2f2FJiKWAH8K/EZmvt3v9pxLEfELwBuZ+US/29IHDeB64J7M/ChwnPlxCGRW1THw7cBa4GpgcUT8an9b9b5wWr91JYTAbDe7n3ciokkrAP4oM79RbX49IlZU+1cAb/SrfefAx4FfjIiXaR3q+7mI+EPmd5/HHQIOZeZ3qvd/QisU5nvffx54KTMHM3MY+AbwM8z/fo+bqZ+n9VtXQgjMdrP7eSUigtbx4Wcy8/Ntu/YAO6rXO4AHz3fbzpXM/O3MXJWZa2j9u/2fmfmrzOM+j8vM14BXI+JD1aYbgKeZ/33/PrA5Ii6u/pu/gdb813zv97iZ+rkHuDEiFkbEWmA98DdzflpmzvsHsA14HngR+J1+t+cc9vNnaQ3/ngS+Vz22AZfTOovghep5eb/beo76/wngL6rXpfT5OmBf9e/8z4FlJfQd+PfAs8BTwB8AC+djv4Gv0Zr3GKb1l/5nZ+sn8DvV79xzwNZevsPLRkhSwUo4HCRJmoEhIEkFMwQkqWCGgCQVzBCQpIIZApJUMENAkgr2/wHSMfQQRlGwuwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"\n",
    "Compute Theta 1 and 2 using Sigma\n",
    "\"\"\"\n",
    "\n",
    "THETA_1 = calc_sigma_1(x_i, x_j)\n",
    "\n",
    "sigma_2 = calc_sigma_2(relu)\n",
    "sigma_2_prime = calc_sigma_2(d_relu)\n",
    "\n",
    "THETA_2 = np.multiply(THETA_1, sigma_2_prime) + sigma_2\n",
    "\n",
    "eigenvalues = np.linalg.eigvals(THETA_2)\n",
    "plt.plot(eigenvalues[1:])"
   ]
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}