Switch off eigensolver per default

doc/gp_custom.ipynb

@@ -27,9 +27,9 @@
     "        for j in np.arange(len(xb)):\n",
     "            K[i, j] = kern_sqexp(xa[i], xb[j], hyp[0])\n",
     "\n",
-    "noise_train = 0.01\n",
+    "noise_train = 0.0\n",
     "\n",
-    "ntrain = 20\n",
+    "ntrain = 30\n",
     "xtrain = halton(1, ntrain)\n",
     "ftrain = f(xtrain)\n",
     "np.random.seed(0)\n",
@@ -50,7 +50,7 @@
     "Ky = K + hyp[-1]*np.eye(ntrain)\n",
     "Kyinv = invert(Ky, 4, 1e-6)       # using gp_functions.invert\n",
     "\n",
-    "ntest = 30\n",
+    "ntest = 20\n",
     "xtest = np.linspace(0, 1, ntest)\n",
     "ftest = f(xtest)\n",
     "\n",
@@ -59,8 +59,7 @@
     "build_K(xtrain, xtest, hyp, Ks)\n",
     "build_K(xtest, xtest, hyp, Kss)\n",
     "\n",
-    "fmean = Ks.T.dot(Kyinv.dot(ytrain)) # predictive mean\n",
-    " "
+    "fmean = Ks.T.dot(Kyinv.dot(ytrain)) # predictive mean"
    ]
   },
   {
@@ -85,7 +84,7 @@
     "# Negative log likelihood over length scale\n",
     "ls = np.linspace(1e-3, 3, 50)\n",
     "nlls = np.array(\n",
-    "    [nll([l, 0.00694534], xtrain, ytrain, 0) for l in ls]\n",
+    "    [nll([l, 1e-3], xtrain, ytrain) for l in ls]\n",
     "    ).flatten()"
    ]
   },
@@ -110,11 +109,28 @@
    "source": [
     "from scipy.optimize import minimize\n",
     "\n",
-    "def nll_transform(log10hyp):\n",
+    "# Prior to cut out range\n",
+    "def cutoff(x, xmin, xmax, slope=1e3):\n",
+    "    if x < xmin:\n",
+    "        return slope*(x - xmin)**2\n",
+    "    if x > xmax:\n",
+    "        return slope*(x - xmax)**2\n",
+    "    \n",
+    "    return 0.0\n",
+    "\n",
+    "def nlprior(log10hyp):\n",
+    "    return cutoff(log10hyp[0], -2, 1) + cutoff(log10hyp[-1], -8, 0)\n",
+    "\n",
+    "x = np.linspace(-10, 1, 100)\n",
+    "plt.figure()\n",
+    "plt.plot(x, [cutoff(xi, -6, 0) for xi in x])\n",
+    "plt.show()\n",
+    "\n",
+    "def nlp_transform(log10hyp):\n",
     "    hyp = 10**log10hyp\n",
-    "    return nll(hyp, xtrain, ytrain, 0)\n",
+    "    return nll(hyp, xtrain, ytrain) + nlprior(log10hyp)\n",
     "\n",
-    "res = minimize(nll_transform, np.array([0, -6]), method='BFGS')"
+    "res = minimize(nlp_transform, np.array([-1, -6]), method='BFGS')"
    ]
   },
   {
@@ -141,7 +157,7 @@
     "[Ll, Ls2] = np.meshgrid(log10l, log10s2)\n",
     "\n",
     "nlls = np.array(\n",
-    "    [nll([10**ll, 10**ls2], xtrain, ytrain, 0) for ls2 in log10s2 for ll in log10l]\n",
+    "    [nlp_transform(np.array([ll, ls2])) for ls2 in log10s2 for ll in log10l]\n",
     "    ).reshape([ns2, nl])\n",
     "\n",
     "# Do some cut for visualization\n",
@@ -165,16 +181,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Trying out priors to cut values\n",
-    "\n",
-    "def sigmoid(x):\n",
-    "    return 1.0/(1.0 + np.exp(-x))\n",
-    "\n",
-    "def prior(hyp):\n",
-    "    return sigmoid(hyp[0]-6)*sigmoid(hyp[-1]-6)\n",
-    "\n",
-    "x = np.logspace(-10, -5, 100)\n",
-    "plt.semilogx(x, np.log(sigmoid(1e9*x - 10)))"
+    "nlls"
    ]
   },
   {

profit/sur/backend/gp_functions.py

@@ -29,15 +29,21 @@ def invert_cholesky(L):
         L.T, solve_triangular(L, np.eye(L.shape[0]), lower=True, check_finite=False),
         lower=False, check_finite=False)
 
-def invert(K, neig=8, tol=1e-10):
+def invert(K, neig=0, tol=1e-10):
     """Inverts a positive-definite matrix A using either an eigendecomposition or
        a Cholesky decomposition, depending on the rapidness of decay of eigenvalues"""
-    if (neig <= 0):
-        return invert_cholesky(np.linalg.cholesky(K))
+    if (neig <= 0 or neig > 0.05*K.shape[0]):
+        try:
+            return invert_cholesky(np.linalg.cholesky(K))
+        except:
+            print('Warning! Fallback to eig solver!')
     w, Q = eigsh(K, neig, tol=tol)
     while np.abs(w[0]-tol) > tol:
         if neig > 0.05*K.shape[0]:  # TODO: get more stringent criterion
-            return invert_cholesky(np.linalg.cholesky(K))
+            try:
+                return invert_cholesky(np.linalg.cholesky(K))
+            except:
+                print('Warning! Fallback to eig solver!')
         neig = 2*neig
         w, Q = eigsh(K, neig, tol=tol)
     return Q.dot(np.diag(1.0/w).dot(Q.T))
@@ -56,18 +62,20 @@ def nll_chol(hyp, x, y, build_K=build_K):
     return ret.item()
 
 
-def nll(hyp, x, y, neig=8, build_K=build_K):
+def nll(hyp, x, y, neig=0, build_K=build_K):
     K = np.empty((len(x), len(x)))
     build_K(x, x, np.abs(hyp[:-1]), K)
     Ky = K + np.abs(hyp[-1])*np.diag(np.ones(len(x)))
-    if (neig <= 0):
-        return nll_chol(hyp, x, y, build_K)
+    if (neig <= 0 or neig > 0.05*len(x)):
+        try:
+            return nll_chol(hyp, x, y, build_K)
+        except:
+            print('Warning! Fallback to eig solver!')
     w, Q = eigsh(Ky, neig, tol=max(1e-6*np.abs(hyp[-1]), 1e-15))
     while np.abs(w[0]-hyp[-1])/hyp[-1] > 1e-6 and neig < len(x):
         if neig > 0.05*len(x):  # TODO: get more stringent criterion
             try:
                 return nll_chol(hyp, x, y, build_K)
-
             except:
                 print('Warning! Fallback to eig solver!')
         neig =  2*neig
@@ -79,7 +87,7 @@ def nll(hyp, x, y, neig=8, build_K=build_K):
     return ret.item()
 
 
-def predict_f(hyp, x, y, xtest, neig=8):
+def predict_f(hyp, x, y, xtest, neig=0):
     Ktest = sklearn.metrics.pairwise.rbf_kernel(xtest, x, 0.5/hyp[0]**2)
     Ktest2 = sklearn.metrics.pairwise.rbf_kernel(xtest, xtest, 0.5/hyp[0]**2)
     K = sklearn.metrics.pairwise.rbf_kernel(x, x, 0.5/hyp[0]**2)