Fix GP problem with scipy versions > 1.4.x

Author: krasserm

gaussian-processes/gaussian_processes.ipynb

@@ -112,17 +112,16 @@
     "import numpy as np\n",
     "\n",
     "def kernel(X1, X2, l=1.0, sigma_f=1.0):\n",
-    "    '''\n",
-    "    Isotropic squared exponential kernel. Computes \n",
-    "    a covariance matrix from points in X1 and X2.\n",
+    "    \"\"\"\n",
+    "    Isotropic squared exponential kernel.\n",
     "    \n",
     "    Args:\n",
     "        X1: Array of m points (m x d).\n",
     "        X2: Array of n points (n x d).\n",
     "\n",
     "    Returns:\n",
-    "        Covariance matrix (m x n).\n",
-    "    '''\n",
+    "        (m x n) matrix.\n",
+    "    \"\"\"\n",
     "    sqdist = np.sum(X1**2, 1).reshape(-1, 1) + np.sum(X2**2, 1) - 2 * np.dot(X1, X2.T)\n",
     "    return sigma_f**2 * np.exp(-0.5 / l**2 * sqdist)"
    ]
@@ -200,7 +199,7 @@
     "from numpy.linalg import inv\n",
     "\n",
     "def posterior_predictive(X_s, X_train, Y_train, l=1.0, sigma_f=1.0, sigma_y=1e-8):\n",
-    "    '''\n",
+    "    \"\"\"\n",
     "    Computes the suffifient statistics of the GP posterior predictive distribution \n",
     "    from m training data X_train and Y_train and n new inputs X_s.\n",
     "    \n",
@@ -214,7 +213,7 @@
     "    \n",
     "    Returns:\n",
     "        Posterior mean vector (n x d) and covariance matrix (n x n).\n",
-    "    '''\n",
+    "    \"\"\"\n",
     "    K = kernel(X_train, X_train, l, sigma_f) + sigma_y**2 * np.eye(len(X_train))\n",
     "    K_s = kernel(X_train, X_s, l, sigma_f)\n",
     "    K_ss = kernel(X_s, X_s, l, sigma_f) + 1e-8 * np.eye(len(X_s))\n",
@@ -398,42 +397,49 @@
     "from scipy.optimize import minimize\n",
     "\n",
     "def nll_fn(X_train, Y_train, noise, naive=True):\n",
-    "    '''\n",
+    "    \"\"\"\n",
     "    Returns a function that computes the negative log marginal\n",
-    "    likelihood for training data X_train and Y_train and given \n",
+    "    likelihood for training data X_train and Y_train and given\n",
     "    noise level.\n",
-    "    \n",
+    "\n",
     "    Args:\n",
     "        X_train: training locations (m x d).\n",
     "        Y_train: training targets (m x 1).\n",
     "        noise: known noise level of Y_train.\n",
-    "        naive: if True use a naive implementation of Eq. (7), if \n",
-    "               False use a numerically more stable implementation. \n",
-    "        \n",
+    "        naive: if True use a naive implementation of Eq. (7), if\n",
+    "               False use a numerically more stable implementation.\n",
+    "\n",
     "    Returns:\n",
     "        Minimization objective.\n",
-    "    '''\n",
+    "    \"\"\"\n",
+    "    \n",
+    "    Y_train = Y_train.ravel()\n",
+    "    \n",
     "    def nll_naive(theta):\n",
     "        # Naive implementation of Eq. (7). Works well for the examples \n",
     "        # in this article but is numerically less stable compared to \n",
     "        # the implementation in nll_stable below.\n",
     "        K = kernel(X_train, X_train, l=theta[0], sigma_f=theta[1]) + \\\n",
     "            noise**2 * np.eye(len(X_train))\n",
     "        return 0.5 * np.log(det(K)) + \\\n",
-    "               0.5 * Y_train.T.dot(inv(K).dot(Y_train)) + \\\n",
+    "               0.5 * Y_train.dot(inv(K).dot(Y_train)) + \\\n",
     "               0.5 * len(X_train) * np.log(2*np.pi)\n",
-    "\n",
+    "        \n",
     "    def nll_stable(theta):\n",
     "        # Numerically more stable implementation of Eq. (7) as described\n",
     "        # in http://www.gaussianprocess.org/gpml/chapters/RW2.pdf, Section\n",
     "        # 2.2, Algorithm 2.1.\n",
+    "        \n",
+    "        def ls(a, b):\n",
+    "            return lstsq(a, b, rcond=-1)[0]\n",
+    "        \n",
     "        K = kernel(X_train, X_train, l=theta[0], sigma_f=theta[1]) + \\\n",
     "            noise**2 * np.eye(len(X_train))\n",
     "        L = cholesky(K)\n",
     "        return np.sum(np.log(np.diagonal(L))) + \\\n",
-    "               0.5 * Y_train.T.dot(lstsq(L.T, lstsq(L, Y_train)[0])[0]) + \\\n",
+    "               0.5 * Y_train.dot(ls(L.T, ls(L, Y_train))) + \\\n",
     "               0.5 * len(X_train) * np.log(2*np.pi)\n",
-    "    \n",
+    "\n",
     "    if naive:\n",
     "        return nll_naive\n",
     "    else:\n",
@@ -717,7 +723,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.9"
+   "version": "3.7.9"
   }
  },
  "nbformat": 4,

gaussian-processes/requirements.txt

@@ -0,0 +1,5 @@
+matplotlib==3.3.1
+numpy==1.19.2
+scipy==1.5.2
+scikit-learn==0.23.2
+GPy==1.9.8