GPy vs Tensorflow.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import gpflow\n",
    "from GPy.util import linalg\n",
    "from gpflow.ci_utils import ci_niter\n",
    "from gpflow.utilities import print_summary\n",
    "import tensorflow as tf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=1000\n",
    "X = np.sort(0.2*np.random.rand(N))[:, None]\n",
    "def true_f(X_input):\n",
    "    return (4.5 * np.cos(2 * np.pi * X_input + 1.5*np.pi) - \\\n",
    "           3 * np.sin(4.3 * np.pi * X_input + 0.3 * np.pi) + \\\n",
    "           5 * np.cos(7 * np.pi * X_input + 2.4 * np.pi))\n",
    "    \n",
    "\n",
    "Y = np.random.normal(true_f(X))\n",
    "data = (X,Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 3  # Number of inducing locations\n",
    "\n",
    "kernel = gpflow.kernels.SquaredExponential()\n",
    "Z = np.linspace(0,max(X),M)  # Initialize inducing locations to the first M inputs in the dataset\n",
    "m = gpflow.models.SVGP(kernel, gpflow.likelihoods.Gaussian(), Z, num_data=N)\n",
    "gpflow.set_trainable(m.inducing_variable, False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_adam(model, iterations):\n",
    "    \"\"\"\n",
    "    Utility function running the Adam optimizer\n",
    "\n",
    "    :param model: GPflow model\n",
    "    :param interations: number of iterations\n",
    "    \"\"\"\n",
    "    # Create an Adam Optimizer action\n",
    "    logf = []\n",
    "    training_loss = model.training_loss_closure(data)\n",
    "    optimizer = tf.optimizers.Adam()\n",
    "\n",
    "    @tf.function\n",
    "    def optimization_step():\n",
    "        optimizer.minimize(training_loss, model.trainable_variables)\n",
    "\n",
    "    for step in range(iterations):\n",
    "        optimization_step()\n",
    "        if step % 10 == 0:\n",
    "            elbo = -training_loss().numpy()\n",
    "            logf.append(elbo)\n",
    "    return logf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "maxiter = ci_niter(10000)\n",
    "run_adam(m, maxiter);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead>\n",
       "<tr><th>name                    </th><th>class    </th><th>transform       </th><th>prior  </th><th>trainable  </th><th>shape    </th><th>dtype  </th><th>value                   </th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "<tr><td>SVGP.kernel.variance    </td><td>Parameter</td><td>Softplus        </td><td>       </td><td>True       </td><td>()       </td><td>float64</td><td>5.215057865828718       </td></tr>\n",
       "<tr><td>SVGP.kernel.lengthscales</td><td>Parameter</td><td>Softplus        </td><td>       </td><td>True       </td><td>()       </td><td>float64</td><td>0.09564411898268795     </td></tr>\n",
       "<tr><td>SVGP.likelihood.variance</td><td>Parameter</td><td>Softplus + Shift</td><td>       </td><td>True       </td><td>()       </td><td>float64</td><td>1.0486152193402796      </td></tr>\n",
       "<tr><td>SVGP.inducing_variable.Z</td><td>Parameter</td><td>                </td><td>       </td><td>False      </td><td>(3, 1)   </td><td>float64</td><td>[[0.        ]\n",
       " [0.09976028]\n",
       " [0.19952057]]                         </td></tr>\n",
       "<tr><td>SVGP.q_mu               </td><td>Parameter</td><td>                </td><td>       </td><td>True       </td><td>(3, 1)   </td><td>float64</td><td>[[-0.54408742]\n",
       " [-1.93647905]\n",
       " [ 7.29369837]]                         </td></tr>\n",
       "<tr><td>SVGP.q_sqrt             </td><td>Parameter</td><td>FillTriangular  </td><td>       </td><td>True       </td><td>(1, 3, 3)</td><td>float64</td><td>[[[0.03728874, 0., 0....</td></tr>\n",
       "</tbody>\n",
       "</table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print_summary(m, fmt=\"notebook\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "Znew = np.linspace(0,2*max(X),2*M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(3, 3), dtype=float64, numpy=\n",
       "array([[2.28365012, 0.        , 0.        ],\n",
       "       [1.32553053, 1.85957702, 0.        ],\n",
       "       [0.25922104, 1.44303922, 1.75085697]])>"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chol = tf.linalg.cholesky(kernel(Z))\n",
    "chol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(6, 3), dtype=float64, numpy=\n",
       "array([[5.21505787e+00, 3.02704796e+00, 5.91970153e-01],\n",
       "       [3.68184618e+00, 5.10281165e+00, 2.38272420e+00],\n",
       "       [1.29564244e+00, 4.28758549e+00, 4.78036171e+00],\n",
       "       [2.27257230e-01, 1.79568049e+00, 4.78036171e+00],\n",
       "       [1.98684171e-02, 3.74850723e-01, 2.38272420e+00],\n",
       "       [8.65808503e-04, 3.90032376e-02, 5.91970153e-01]])>"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "conditional = kernel(Znew,Z)\n",
    "conditional"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.00000000e+00,  1.71201884e-01, -1.37323609e-01,\n",
       "         1.44960851e-01,  1.73465418e-01,  5.54698263e-02],\n",
       "       [ 2.42658650e-16,  9.42839390e-01,  5.44054288e-01,\n",
       "        -3.95614327e-01, -4.26159590e-01, -1.31130352e-01],\n",
       "       [-1.12430621e-16, -1.09805542e-01,  6.16440862e-01,\n",
       "         1.12982308e+00,  6.84564478e-01,  1.83329037e-01]])"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "linalg.dpotrs(np.asfortranarray(chol.numpy()), conditional.numpy().T)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(3, 6), dtype=float64, numpy=\n",
       "array([[ 1.00000000e+00,  1.71201884e-01, -1.37323609e-01,\n",
       "         1.44960851e-01,  1.73465418e-01,  5.54698263e-02],\n",
       "       [ 1.87554634e-16,  9.42839390e-01,  5.44054288e-01,\n",
       "        -3.95614327e-01, -4.26159590e-01, -1.31130352e-01],\n",
       "       [-7.62005196e-17, -1.09805542e-01,  6.16440862e-01,\n",
       "         1.12982308e+00,  6.84564478e-01,  1.83329037e-01]])>"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf.linalg.cholesky_solve(chol.numpy(), conditional.numpy().T)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:root] *",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}