logger changed and option to separate appending data and rank-n update

fvgp/fvgp.py

@@ -151,12 +151,10 @@ class provides all the methods described for the GP (:py:class:`fvgp.GP`) class.
     calc_inv : bool, optional
         If True, the algorithm calculates and stores the inverse of the covariance
         matrix after each training or update of the dataset or hyperparameters,
-        which makes computing the posterior covariance faster (5-10 times).
-        For larger problems (>2000 data points), the use of inversion should be avoided due
-        to computational instability and costs. The default is
-        False. Note, the training will not the
+        which makes computing the posterior covariance faster (3-10 times).
+        The default is False. Note, the training will not use the
         inverse for stability reasons. Storing the inverse is
-        a good option when the dataset is not too large and the posterior covariance is heavily used.
+        a good option when the posterior covariance is heavily used.
     ram_economy : bool, optional
         Only of interest if the gradient and/or Hessian of the log marginal likelihood is/are used for the training.
         If True, components of the derivative of the log marginal likelihood are

fvgp/gp.py

@@ -142,12 +142,10 @@ class GP:
     calc_inv : bool, optional
         If True, the algorithm calculates and stores the inverse of the covariance
         matrix after each training or update of the dataset or hyperparameters,
-        which makes computing the posterior covariance faster (5-10 times).
-        For larger problems (>2000 data points), the use of inversion should be avoided due
-        to computational instability and costs. The default is
-        False. Note, the training will not the
+        which makes computing the posterior covariance faster (3-10 times).
+        The default is False. Note, the training will not use the
         inverse for stability reasons. Storing the inverse is
-        a good option when the dataset is not too large and the posterior covariance is heavily used.
+        a good option when the posterior covariance is heavily used.
     ram_economy : bool, optional
         Only of interest if the gradient and/or Hessian of the log marginal likelihood is/are used for the training.
         If True, components of the derivative of the log marginal likelihood are
@@ -326,7 +324,8 @@ def update_gp_data(
         x_new,
         y_new,
         noise_variances_new=None,
-        append=True
+        append=True,
+        gp_rank_n_update=None
     ):
         """
         This function updates the data in the gp object instance.
@@ -351,8 +350,13 @@ def update_gp_data(
         append : bool, optional
             Indication whether to append to or overwrite the existing dataset. Default=True.
             In the default case, data will be appended.
+        gp_rank_n_update : bool, optional
+            Indicates whether the GP marginal should be rank-n updated or recomputed. The default
+            is `gp_rank_n_update=append`, meaning if data is only appended, the rang_n_update will
+            be performed.
         """
         old_x_data = self.data.x_data.copy()
+        if gp_rank_n_update is None: gp_rank_n_update = append
         # update data
         self.data.update(x_new, y_new, noise_variances_new, append=append)
 
@@ -365,7 +369,7 @@ def update_gp_data(
                                self.prior.hyperparameters)
 
         # update marginal density
-        self.marginal_density.update_data(append)
+        self.marginal_density.update_data(gp_rank_n_update)
         ##########################################
         self.x_data = self.data.x_data
         self.y_data = self.data.y_data
@@ -840,7 +844,7 @@ def posterior_covariance(self, x_pred, x_out=None, variance_only=False, add_nois
         variance_only : bool, optional
             If True the computation of the posterior covariance matrix is avoided which can save compute time.
             In that case the return will only provide the variance at the input points.
-            Default = False.
+            Default = False. This is only relevant if `calc_inv` at initialization is True.
         add_noise : bool, optional
             If True the noise variances will be added to the posterior variances. Default = False.
 

fvgp/gp_lin_alg.py

@@ -93,7 +93,7 @@ def calculate_random_logdet(KV, info, compute_device):
 
 
 def calculate_sparse_conj_grad(KV, vec, info=False):
-    logger.debug("calculate_sparse_conj_grad")
+    logger.info("calculate_sparse_conj_grad")
     st = time.time()
     if info: logger.info("CG solve in progress ...")
     if np.ndim(vec) == 1: vec = vec.reshape(len(vec), 1)
@@ -108,13 +108,14 @@ def calculate_sparse_conj_grad(KV, vec, info=False):
     if info: logger.info("CG compute time: {} seconds, exit status {} (0:=successful)", time.time() - st, exit_code)
     return res
 
+
 def update_sparse_conj_grad(KV, vec, x0, info=False):
     assert np.ndim(vec) == 1
     assert np.ndim(KV) == 2
     assert np.ndim(x0) == 1
-    logger.debug("update_sparse_conj_grad")
+    logger.info("update_sparse_conj_grad")
     st = time.time()
-    if len(x0)<KV.shape[0]: x0=np.append(x0,np.zeros(KV.shape[0]-len(x0)))
+    if len(x0) < KV.shape[0]: x0 = np.append(x0, np.zeros(KV.shape[0] - len(x0)))
     if info: logger.info("CG solve in progress ...")
     vec = vec.reshape(len(vec), 1)
     res, exit_code = cg(KV.tocsc(), vec[:, 0], rtol=1e-8, x0=x0)
@@ -128,7 +129,7 @@ def update_sparse_conj_grad(KV, vec, x0, info=False):
 
 
 def calculate_sparse_minres(KV, vec, info=False):
-    logger.debug("calculate_sparse_minres")
+    logger.info("calculate_sparse_minres")
     st = time.time()
     if info: logger.info("MINRES solve in progress ...")
     if np.ndim(vec) == 1: vec = vec.reshape(len(vec), 1)
@@ -145,9 +146,9 @@ def update_sparse_minres(KV, vec, x0, info=False):
     assert np.ndim(vec) == 1
     assert np.ndim(KV) == 2
     assert np.ndim(x0) == 1
-    logger.debug("update_sparse_minres")
+    logger.info("update_sparse_minres")
     st = time.time()
-    if len(x0)<KV.shape[0]: x0=np.append(x0,np.zeros(KV.shape[0]-len(x0)))
+    if len(x0) < KV.shape[0]: x0 = np.append(x0, np.zeros(KV.shape[0] - len(x0)))
     if info: logger.info("MINRES solve in progress ...")
     vec = vec.reshape(len(vec), 1)
     res, exit_code = minres(KV.tocsc(), vec[:, 0], rtol=1e-8, x0=x0)
@@ -195,7 +196,7 @@ def cholesky_update_rank_n(L, b, c):
 
 
 def calculate_logdet(A, compute_device='cpu'):
-    logger.debug("calculate_logdet")
+    logger.info("calculate_logdet")
     if compute_device == "cpu":
         s, logdet = np.linalg.slogdet(A)
         return logdet
@@ -218,7 +219,7 @@ def calculate_logdet(A, compute_device='cpu'):
 
 
 def update_logdet(old_logdet, old_inv, new_matrix, compute_device="cpu"):
-    logger.debug("update_logdet")
+    logger.info("update_logdet")
     size = len(old_inv)
     KV = new_matrix
     kk = KV[size:, size:]
@@ -228,7 +229,7 @@ def update_logdet(old_logdet, old_inv, new_matrix, compute_device="cpu"):
 
 
 def calculate_inv(A, compute_device='cpu'):
-    logger.debug("calculate_inv")
+    logger.info("calculate_inv")
     if compute_device == "cpu":
         return np.linalg.inv(A)
     elif compute_device == "gpu":
@@ -241,7 +242,7 @@ def calculate_inv(A, compute_device='cpu'):
 
 
 def update_inv(old_inv, new_matrix, compute_device="cpu"):
-    logger.debug("update_inv")
+    logger.info("update_inv")
     size = len(old_inv)
     KV = new_matrix
     kk = KV[size:, size:]
@@ -254,7 +255,7 @@ def update_inv(old_inv, new_matrix, compute_device="cpu"):
 
 
 def solve(A, b, compute_device='cpu'):
-    logger.debug("solve")
+    logger.info("solve")
     if np.ndim(b) == 1: b = np.expand_dims(b, axis=1)
     if compute_device == "cpu":
         try:
@@ -307,13 +308,13 @@ def solve(A, b, compute_device='cpu'):
 
 ##################################################################################
 def is_sparse(A):
-    logger.debug("is_sparse")
+    logger.info("is_sparse")
     if float(np.count_nonzero(A)) / float(len(A) ** 2) < 0.01:
         return True
     else:
         return False
 
 
 def how_sparse_is(A):
-    logger.debug("how_sparse_is")
+    logger.info("how_sparse_is")
     return float(np.count_nonzero(A)) / float(len(A) ** 2)

fvgp/gp_marginal_density.py

@@ -56,7 +56,7 @@ def _update_KVinvY(self, K, V, m):
         y_mean = self.y_data - m
         KV = K + V
         self.KVlinalg.update_KV(KV)
-        KVinvY = self.KVlinalg.solve(y_mean, x0 = self.KVinvY)
+        KVinvY = self.KVlinalg.solve(y_mean, x0=self.KVinvY)
         return KVinvY.reshape(len(y_mean))
 
     def _set_KVinvY(self, K, V, m, mode):

fvgp/gp_posterior.py

@@ -97,7 +97,8 @@ def posterior_covariance(self, x_pred, x_out=None, variance_only=False, add_nois
         if self.marginal_density_obj.KVlinalg.KVinv is not None:
             if variance_only:
                 S = None
-                v = np.diag(kk) - np.einsum('ij,jk,ki->i', k.T, self.marginal_density_obj.KVlinalg.KVinv, k)
+                v = np.diag(kk) - np.einsum('ij,jk,ki->i', k.T,
+                                            self.marginal_density_obj.KVlinalg.KVinv, k, optimize=True)
             else:
                 S = kk - (k.T @ self.marginal_density_obj.KVlinalg.KVinv @ k)
                 v = np.array(np.diag(S))