Bugs and performance issues related to PCov-ized sample selection methods (#118)

* Possible bug in the pcovr orthogonalization concept.

Loooots of printf debug statements. Also changed regression to lstsq which makes the y orthogonalization 10x faster

* Force y to be 2D in pcov selection routines

Author: ceriottm

skcosmo/_selection.py

@@ -8,8 +8,8 @@
 
 import numpy as np
 import scipy
-from scipy.linalg import eig
-from scipy.sparse.linalg import eigs as speig
+from scipy.linalg import eigh
+from scipy.sparse.linalg import eigsh
 from sklearn.base import (
     BaseEstimator,
     MetaEstimatorMixin,
@@ -142,6 +142,11 @@ def fit(self, X, y=None, warm_start=False):
                 force_all_finite=not tags.get("allow_nan", True),
                 multi_output=True,
             )
+            if len(y.shape) == 1:
+                # force y to have multi_output 2D format even when it's 1D, since
+                # many functions, most notably PCov routines, assume an array storage
+                # format, most notably to compute (y @ y.T)
+                y = y.reshape((len(y), 1))
         else:
             X = check_array(
                 X,
@@ -659,6 +664,8 @@ def _init_greedy_search(self, X, y, n_to_select):
         features and computes their initial importance score.
         """
 
+        self.X_ref_ = X
+        self.y_ref_ = y
         self.X_current_ = X.copy()
         if y is not None:
             self.y_current_ = y.copy()
@@ -760,15 +767,16 @@ def _compute_pi(self, X, y=None):
             )
 
         if self.k < pcovr_distance.shape[0] - 1:
-            v, U = speig(pcovr_distance, k=self.k, tol=1e-12)
+            v, U = eigsh(pcovr_distance, k=self.k, tol=1e-12)
         else:
-            v, U = eig(pcovr_distance)
+            v, U = eigh(pcovr_distance)
         U = U[:, np.flip(np.argsort(v))]
         pi = (np.real(U)[:, : self.k] ** 2.0).sum(axis=1)
 
         return pi
 
     def _orthogonalize(self, last_selected):
+
         if self._axis == 1:
             self.X_current_ = X_orthogonalizer(
                 x1=self.X_current_, c=last_selected, tol=self.tolerance
@@ -777,18 +785,17 @@ def _orthogonalize(self, last_selected):
             self.X_current_ = X_orthogonalizer(
                 x1=self.X_current_.T, c=last_selected, tol=self.tolerance
             ).T
-
         if self.y_current_ is not None:
             if self._axis == 1:
                 self.y_current_ = Y_feature_orthogonalizer(
                     self.y_current_, X=self.X_selected_, tol=self.tolerance
                 )
             else:
                 self.y_current_ = Y_sample_orthogonalizer(
-                    self.y_current_,
-                    self.X_current_,
-                    y_ref=self.y_selected_,
-                    X_ref=self.X_selected_,
+                    self.y_ref_,
+                    self.X_ref_,
+                    y_ref=self.y_selected_[: self.n_selected_],
+                    X_ref=self.X_selected_[: self.n_selected_],
                     tol=self.tolerance,
                 )
 

skcosmo/utils/_orthogonalizers.py

@@ -135,9 +135,7 @@ def Y_sample_orthogonalizer(y, X, y_ref, X_ref, tol=1e-12, copy=True):
 
     """
 
-    y_frag = (
-        X @ (np.linalg.pinv(X_ref.T @ X_ref, rcond=tol) @ X_ref.T) @ y_ref
-    ).reshape(y.shape)
+    y_frag = (X @ (np.linalg.lstsq(X_ref, y_ref, rcond=tol)[0])).reshape(y.shape)
 
     if copy:
         return y.copy() - y_frag

tests/test_sample_pcov_cur.py

@@ -11,7 +11,7 @@
 class TestPCovCUR(unittest.TestCase):
     def setUp(self):
         self.X, self.y = load_boston(return_X_y=True)
-        self.idx = [488, 283, 183, 380, 41, 438, 368, 374, 123, 353]
+        self.idx = [492, 450, 183, 199, 380, 228, 399, 126, 412, 368]
 
     def test_known(self):
         """
@@ -50,7 +50,7 @@ def test_non_it(self):
         """
         This test checks that the model can be run non-iteratively
         """
-        self.idx = [488, 492, 491, 374, 398, 373, 386, 400, 383, 382]
+        self.idx = [492, 488, 491, 489, 374, 373, 386, 398, 383, 382]
         selector = PCovCUR(n_to_select=10, iterative=False)
         selector.fit(self.X, self.y)
 

tests/test_sample_pcov_fps.py

@@ -8,7 +8,7 @@
 class TestPCovFPS(unittest.TestCase):
     def setUp(self):
         self.X, self.y = load_boston(return_X_y=True)
-        self.idx = [39, 410, 492, 102, 54, 413, 34, 346, 126, 134, 433, 380]
+        self.idx = [39, 410, 492, 102, 54, 413, 34, 126, 346, 134, 433, 380]
 
     def test_restart(self):
         """