Add NIPALS vector difference norm calculation (#10)

* Add norm calculation for NIPALS score vector convergence

When using the NIPALS algorithm the difference of subsequent score
vectors is now calculated by using the second order norm of the
difference vector. i.e. the frobenius norm as default. Optionally, other
norms can be used as defined in the function's docstring.

* Update new test data to match precise decimals of NIPALS change

* Add scikit-learn array copy behaviour to MBPLS methods

Author: lvermue

mbpls/mbpls.py

@@ -51,6 +51,26 @@ class MBPLS(TransformerMixin, RegressorMixin, MultiOutputMixin, BaseEstimator):
         max_tol : non-negative float (default 1e-14)
             Maximum tolerance allowed when using the iterative NIPALS algorithm
 
+        nipals_convergence_norm : {non-zero int, inf, -inf, 'fro', 'nuc'} (default 2)
+            Order of the norm that is used to calculate the difference of the superscore vectors between subsequent
+            iterations of the NIPALS algorithm. Following orders are available:
+
+            =====  ============================  ==========================
+            ord    norm for matrices             norm for vectors
+            =====  ============================  ==========================
+            None   Frobenius norm                2-norm
+            'fro'  Frobenius norm                --
+            'nuc'  nuclear norm                  --
+            inf    max(sum(abs(x), axis=1))      max(abs(x))
+            -inf   min(sum(abs(x), axis=1))      min(abs(x))
+            0      --                            sum(x != 0)
+            1      max(sum(abs(x), axis=0))      as below
+            -1     min(sum(abs(x), axis=0))      as below
+            2      2-norm (largest sing. value)  as below
+            -2     smallest singular value       as below
+            other  --                            sum(abs(x)**ord)**(1./ord)
+            =====  ============================  ==========================
+
         calc_all : bool (default True)
             Calculate all internal attributes for the used method. Some methods do not need to calculate all attributes,
             i.e. scores, weights etc., to obtain the regression coefficients used for prediction. Setting this parameter
@@ -62,7 +82,11 @@ class MBPLS(TransformerMixin, RegressorMixin, MultiOutputMixin, BaseEstimator):
             allowed.
             Without setting this parameter to 'True', sparse data will not be accepted.
 
-        
+        copy : bool (default True)
+            Whether the deflation should be done on a copy. Not using a copy might alter the input data and have
+            unforeseeable consequences.
+
+
         Model attributes after fitting
         ------------------------------
 
@@ -97,7 +121,6 @@ class MBPLS(TransformerMixin, RegressorMixin, MultiOutputMixin, BaseEstimator):
         
         explained_var_y_ : list, explained variance in :math:`Y` :math:`[k]`
         
-        
 
         Notes
         -----
@@ -216,15 +239,17 @@ class MBPLS(TransformerMixin, RegressorMixin, MultiOutputMixin, BaseEstimator):
         https://github.com/DTUComputeStatisticsAndDataAnalysis/MBPLS/tree/master/examples
     """
 
-    def __init__(self, n_components=2, full_svd=False, method='NIPALS', standardize=True, max_tol=1e-14, calc_all=True,
-                 sparse_data=False):
+    def __init__(self, n_components=2, full_svd=False, method='NIPALS', standardize=True, max_tol=1e-14,
+                 nipals_convergence_norm=2, calc_all=True, sparse_data=False, copy=True):
         self.n_components = n_components
         self.full_svd = full_svd
         self.method = method
         self.standardize = standardize
         self.max_tol = max_tol
+        self.nipals_convergence_norm = nipals_convergence_norm
         self.calc_all = calc_all
         self.sparse_data = sparse_data
+        self.copy = copy
 
     def check_sparsity_level(self, data):
         total_rows, total_columns = data.shape
@@ -264,7 +289,7 @@ def fit(self, X, Y):
                 self.method = 'NIPALS'
 
         global U_, T_, R_
-        Y = check_array(Y, dtype=np.float64, ensure_2d=False, force_all_finite=not self.sparse_data)
+        Y = check_array(Y, dtype=np.float64, ensure_2d=False, force_all_finite=not self.sparse_data, copy=self.copy)
         if self.sparse_data is True:
             self.sparse_Y_info_ = {}
             self.sparse_Y_info_['Y'] = self.check_sparsity_level(Y)
@@ -279,14 +304,15 @@ def fit(self, X, Y):
                     self.x_scalers_.append(StandardScaler(with_mean=True, with_std=True))
                     # Check dimensions
                     check_consistent_length(X[block], Y)
-                    X[block] = check_array(X[block], dtype=np.float64, copy=True, force_all_finite=not self.sparse_data)
+                    X[block] = check_array(X[block], dtype=np.float64, copy=self.copy,
+                                           force_all_finite=not self.sparse_data)
                     if self.sparse_data is True:
                         self.sparse_X_info_[block] = self.check_sparsity_level(X[block])
                     X[block] = self.x_scalers_[block].fit_transform(X[block])
             else:
                 self.x_scalers_.append(StandardScaler(with_mean=True, with_std=True))
                 # Check dimensions
-                X = check_array(X, dtype=np.float64, copy=True, force_all_finite=not self.sparse_data)
+                X = check_array(X, dtype=np.float64, copy=self.copy, force_all_finite=not self.sparse_data)
                 if self.sparse_data is True:
                     self.sparse_X_info_ = {}
                     self.sparse_X_info_[0] = self.check_sparsity_level(X)
@@ -302,12 +328,13 @@ def fit(self, X, Y):
                 for block in range(len(X)):
                     # Check dimensions
                     check_consistent_length(X[block], Y)
-                    X[block] = check_array(X[block], dtype=np.float64, copy=True, force_all_finite=not self.sparse_data)
+                    X[block] = check_array(X[block], dtype=np.float64, copy=self.copy,
+                                           force_all_finite=not self.sparse_data)
                     if self.sparse_data is True:
                         self.sparse_X_info_[block] = self.check_sparsity_level(X[block])
             else:
                 # Check dimensions
-                X = check_array(X, dtype=np.float64, copy=True, force_all_finite=not self.sparse_data)
+                X = check_array(X, dtype=np.float64, copy=self.copy, force_all_finite=not self.sparse_data)
                 if self.sparse_data is True:
                     self.sparse_X_info_ = {}
                     self.sparse_X_info_[0] = self.check_sparsity_level(X)
@@ -852,7 +879,7 @@ def fit(self, X, Y):
                     if run == 1:
                         pass
                     else:
-                        diff_t = np.sum(superscores_old - superscores)
+                        diff_t = np.linalg.norm((superscores_old - superscores), ord=self.nipals_convergence_norm)
                     superscores_old = np.copy(superscores)
                     # 6. Regress superscores agains Y_calc
                     if self.sparse_data:
@@ -1017,27 +1044,35 @@ def fit(self, X, Y):
         else:
             raise NameError('Method you called is unknown')
 
-    def transform(self, X, Y=None, return_block_scores=False):
+    def transform(self, X, Y=None, return_block_scores=False, copy=True):
         """ Obtain scores based on the fitted model
 
+
          Parameters
         ----------
         X : list
             of arrays containing all xblocks x1, x2, ..., xn. Rows are observations, columns are features/variables
+
         (optional) Y : array
             1-dim or 2-dim array of reference values
+
         return_block_scores: bool (default False)
             Returning block scores T_ when transforming the data
 
+        copy : bool (default True)
+            Whether to perform in-place transformation. Not using a copy might alter the input data and have
+            unforeseeable consequences.
+
+
         Returns
         ----------
         Super_scores : np.array
 
         Block_scores : list
-        List of np.arrays containing the block scores
+            List of np.arrays containing the block scores
 
         Y_scores : np.array (optional)
-        Y-scores, if y was given
+            Y-scores, if y was given
         """
         check_is_fitted(self, 'beta_')
 
@@ -1049,13 +1084,13 @@ def transform(self, X, Y=None, return_block_scores=False):
             if isinstance(X, list) and not isinstance(X[0], list):
                 for block in range(len(X)):
                     # Check dimensions
-                    X[block] = check_array(X[block], dtype=np.float64, force_all_finite=not self.sparse_data)
+                    X[block] = check_array(X[block], dtype=np.float64, force_all_finite=not self.sparse_data, copy=copy)
                     if self.sparse_data:
                         sparse_X_info_[block] = self.check_sparsity_level(X[block])
                     X[block] = self.x_scalers_[block].transform(X[block])
             else:
                 # Check dimensions
-                X = check_array(X, dtype=np.float64, force_all_finite=not self.sparse_data)
+                X = check_array(X, dtype=np.float64, force_all_finite=not self.sparse_data, copy=copy)
                 if self.sparse_data:
                     sparse_X_info_[0] = self.check_sparsity_level(X)
                 X = [self.x_scalers_[0].transform(X)]
@@ -1075,7 +1110,7 @@ def transform(self, X, Y=None, return_block_scores=False):
                 Ts_ = X_comp.dot(self.R_)
 
             if Y is not None:
-                Y = check_array(Y, dtype=np.float64, ensure_2d=False, force_all_finite=not self.sparse_data)
+                Y = check_array(Y, dtype=np.float64, ensure_2d=False, force_all_finite=not self.sparse_data, copy=copy)
                 if self.sparse_data:
                     sparse_Y_info_['Y'] = self.check_sparsity_level(Y)
                 if Y.ndim == 1:
@@ -1181,12 +1216,12 @@ def transform(self, X, Y=None, return_block_scores=False):
             if isinstance(X, list) and not isinstance(X[0], list):
                 for block in range(len(X)):
                     # Check dimensions
-                    X[block] = check_array(X[block], dtype=np.float64, force_all_finite=not self.sparse_data)
+                    X[block] = check_array(X[block], dtype=np.float64, force_all_finite=not self.sparse_data, copy=copy)
                     if self.sparse_data:
                         sparse_X_info_[block] = self.check_sparsity_level(X[block])
             else:
                 # Check dimensions
-                X = check_array(X, dtype=np.float64, force_all_finite=not self.sparse_data)
+                X = check_array(X, dtype=np.float64, force_all_finite=not self.sparse_data, copy=copy)
                 if self.sparse_data:
                     sparse_X_info_[0] = self.check_sparsity_level(X)
                 X = [X]
@@ -1206,7 +1241,7 @@ def transform(self, X, Y=None, return_block_scores=False):
                 Ts_ = X_comp.dot(self.R_)
 
             if Y is not None:
-                Y = check_array(Y, dtype=np.float64, ensure_2d=False, force_all_finite=not self.sparse_data)
+                Y = check_array(Y, dtype=np.float64, ensure_2d=False, force_all_finite=not self.sparse_data, copy=copy)
                 if self.sparse_data:
                     sparse_Y_info_['Y'] = self.check_sparsity_level(Y)
                 if Y.ndim == 1:
@@ -1290,19 +1325,25 @@ def transform(self, X, Y=None, return_block_scores=False):
                 else:
                     return Ts_
 
-    def predict(self, X):
+    def predict(self, X, copy=True):
         """Predict y based on the fitted model
 
+
         Parameters
         ----------
+
         X : list
             of all xblocks x1, x2, ..., xn. Rows are observations, columns are features/variables
 
+        copy : bool (default True)
+            Whether to perform in-place transformation. Not using a copy might alter the input data and have
+            unforeseeable consequences.
+
+
         Returns
         ----------
         y_hat : np.array
-        Predictions made based on trained model and supplied X
-
+            Predictions made based on trained model and supplied X
         """
         check_is_fitted(self, 'beta_')
 
@@ -1313,10 +1354,10 @@ def predict(self, X):
             if isinstance(X, list) and not isinstance(X[0], list):
                 for block in range(len(X)):
                     # Check dimensions
-                    X[block] = check_array(X[block], dtype=np.float64, force_all_finite=not self.sparse_data)
+                    X[block] = check_array(X[block], dtype=np.float64, force_all_finite=not self.sparse_data, copy=copy)
                     X[block] = self.x_scalers_[block].transform(X[block])
             else:
-                X = check_array(X, dtype=np.float64, force_all_finite=not self.sparse_data)
+                X = check_array(X, dtype=np.float64, force_all_finite=not self.sparse_data, copy=copy)
                 X = [self.x_scalers_[0].transform(X)]
 
 
@@ -1336,9 +1377,9 @@ def predict(self, X):
             if isinstance(X, list) and not isinstance(X[0], list):
                 for block in range(len(X)):
                     # Check dimensions
-                    X[block] = check_array(X[block], dtype=np.float64, force_all_finite=not self.sparse_data)
+                    X[block] = check_array(X[block], dtype=np.float64, force_all_finite=not self.sparse_data, copy=copy)
             else:
-                X = check_array(X, dtype=np.float64, force_all_finite=not self.sparse_data)
+                X = check_array(X, dtype=np.float64, force_all_finite=not self.sparse_data, copy=copy)
                 X = [X]
             X = np.hstack(X)
             if self.sparse_data:

mbpls/tests/test_data/A_NIPALS.csv

@@ -1,2 +1,2 @@
-9.460050930379139134e-01,3.864865617219123001e-02
-5.399490696208620460e-02,9.613513438278087353e-01
+9.274646237881400967e-01,5.716118652731461136e-02
+7.253537621185997264e-02,9.428388134726853886e-01

mbpls/tests/test_data/P1_NIPALS.csv

@@ -1,25 +1,25 @@
-5.808427406751311750e+00,6.520579921256950406e-01
-5.926596442662182440e+00,3.414058821124497101e-01
-6.086453323360574430e+00,4.960972689404420932e-01
-5.673416006491299513e+00,2.196348688272564686e-02
-6.059110719312798210e+00,5.575722505091080805e-01
-5.650942259585248095e+00,8.035470872887177096e-01
-6.193024415711314568e+00,3.761258217688800976e-01
-5.718477717250506132e+00,1.180765043320053109e-01
-5.974501227518040025e+00,1.105125183151397239e-01
-3.851918076221314990e+00,3.173299523953762025e-01
-6.038798238605442847e+00,1.152914973439507662e+00
-5.346208049272602736e+00,1.285607400057692695e+00
-6.140535096631946743e+00,6.281127826818750925e-01
-6.094854212815943306e+00,6.500160775864066709e-01
-5.111949514516579640e+00,1.041565457794748850e+00
-5.862289305226957126e+00,5.993825610514952329e-02
-5.608125399568422154e+00,-1.001429840936495008e+00
-1.626524804393373547e+00,9.296548505843652555e-01
-5.993136960784692491e+00,8.656563868700076769e-01
--6.865050824961729248e-01,-1.347598947582476225e-01
-6.157950800160661764e+00,7.322040498902356864e-01
-1.147230076524884401e-01,8.922509625715264736e-01
-5.348783867045700191e+00,7.416084505076506739e-01
-5.884084660201166272e+00,1.894637312204782709e-01
-5.689278078360489488e+00,1.176970173160049482e+00
+5.778141171141287558e+00,8.807390506152754650e-01
+5.908470722787161122e+00,5.750364384512615112e-01
+6.062100478068037646e+00,7.358715443699839209e-01
+5.668136154400002624e+00,2.457974398229465507e-01
+6.032348050525556893e+00,7.961299169925415420e-01
+5.614807645731139552e+00,1.025857546053548397e+00
+6.173329201680968659e+00,6.202103429515453126e-01
+5.709381844946635809e+00,3.437633844666128646e-01
+5.965492877132286154e+00,3.462350148801632832e-01
+3.836469051155537269e+00,4.693360014527241386e-01
+5.988570606931146401e+00,1.390411268493505270e+00
+5.291293080231936230e+00,1.495354912435294281e+00
+6.110934167362700897e+00,8.698676960424125415e-01
+6.064406950376159067e+00,8.899853132326137617e-01
+5.066888370390126184e+00,1.242888876859258307e+00
+5.855372046908557948e+00,2.912183806269581510e-01
+5.643261210523061067e+00,-7.792454675832711786e-01
+1.588547991938236148e+00,9.930027998805039946e-01
+5.954338133724771609e+00,1.101669613466050412e+00
+-6.806038670955768533e-01,-1.622982367317185592e-01
+6.124247472541748394e+00,9.747324336857725591e-01
+7.941412319847895862e-02,8.956274678154049207e-01
+5.315367390122011315e+00,9.521689148940255532e-01
+5.871993582189810823e+00,4.215876038627815459e-01
+5.638359679673937919e+00,1.400667425441476599e+00