Misc comments, enforce diagonal, etc (#165)

* Re-use indexed column
* Comment that Prec_r is always diagonal
* Fail if Prec_eps is not diagonal - in that case math does not check out
* Only log posterior precision log det if logging.INFO

Author: tommyod

graphite_maps/enif.py

@@ -169,6 +169,7 @@ def transport(
 
         # Work out residuals and associate unexplained variance
         residuals = self.response_residual(U, Y)
+
         # Due to observation error
         eps = self.generate_observation_noise(
             n,
@@ -178,9 +179,9 @@ def transport(
 
         # Update in canonical parametrization
         canonical_updated = self.update_canonical(
-            canonical,
-            residual_noisy,
-            d,
+            canonical=canonical,
+            residual_noisy=residual_noisy,
+            d=d,
         )
 
         # Bring realizations back
@@ -253,6 +254,11 @@ def Prec_residual_noisy(self) -> sparray:
         if self.unexplained_variance is None:
             raise ValueError("`unexplained_variance` is not set.")
 
+        # The equation below is only valid if Prec_eps is diagonal
+        row_idx, col_idx, _ = sp.sparse.find(self.Prec_eps)
+        if np.any(row_idx != col_idx):
+            raise ValueError("Precision matrix 'Prec_eps' must be diagonal")
+
         eps_variances = 1.0 / self.Prec_eps.diagonal()
         residual_noisy_var = self.unexplained_variance + eps_variances
         Prec_r = diags_array(1.0 / residual_noisy_var, offsets=0, format="csc")
@@ -327,7 +333,7 @@ def update_canonical(
             logdet_value = 2.0 * np.sum(np.log(chol_LLT.L().diagonal()))
             log.info("Prior precision log-determinant: %.3f", logdet_value)
 
-        Prec_r = self.Prec_residual_noisy()
+        Prec_r = self.Prec_residual_noisy()  # This is a diagonal matrix
 
         # posterior canonical params
         # this is equation (46), but transposed to update each row (realizations)
@@ -338,9 +344,10 @@ def update_canonical(
         # posterior precision, equation (47)
         self.Prec_u = self.Prec_u + self.H.T @ Prec_r @ self.H  # Eqn (47)
 
-        chol_LLT = cholesky(self.Prec_u, ordering_method="metis")
-        logdet_value = 2.0 * np.sum(np.log(chol_LLT.L().diagonal()))
-        log.info("Posterior precision log-determinant: %.3f", logdet_value)
+        if log.isEnabledFor(logging.INFO):
+            chol_LLT = cholesky(self.Prec_u, ordering_method="metis")
+            logdet_value = 2.0 * np.sum(np.log(chol_LLT.L().diagonal()))
+            log.info("Posterior precision log-determinant: %.3f", logdet_value)
 
         return updated_canonical
 

graphite_maps/linear_regression.py

@@ -193,7 +193,7 @@ def boost_linear_regression(
     """
     n_samples, n_features = X.shape
     coefficients = np.zeros(n_features)
-    residuals = y.copy()
+    residuals = y.copy()  # residuals = y - X @ coef = y - X @ 0 = y
     residuals_loo = y.copy()
 
     # A stricter criterion is the loo-adjustment: mse(residuals_loo)-mse
@@ -220,15 +220,12 @@ def boost_linear_regression(
         # Inlined influence: psi / M
         # where influence = -residuals * x and
         # M = -mean(x^2) = -1 since x standardized
-        influence = (residuals - beta_estimate * X[:, best_feature]) * X[
-            :, best_feature
-        ]
+        X_best = X[:, best_feature]
+        influence = (residuals - beta_estimate * X_best) * X_best
         beta_estimate_loo = beta_estimate - influence / n_samples
 
         # residuals_full = residuals - beta_estimate * X[:, best_feature]
-        residuals_full_loo = (
-            residuals_loo - learning_rate * beta_estimate_loo * X[:, best_feature]
-        )
+        residuals_full_loo = residuals_loo - learning_rate * beta_estimate_loo * X_best
 
         if mse(residuals_loo) < mse(residuals_full_loo):
             break
@@ -246,7 +243,7 @@ def boost_linear_regression(
             break
         else:
             # Update
-            residuals -= coef_change * X[:, best_feature]
+            residuals -= coef_change * X_best
             coefficients[best_feature] += coef_change
 
             # loo update