Merge pull request scikit-learn#2396 from dengemann/insert_fast_dot

ENH: put fast dot to PCA, ICA and FA

Author: ogrisel

doc/developers/performance.rst

@@ -98,7 +98,7 @@ replacement for `numpy.dot`. example::
 
   >>> X = np.random.random_sample([2, 10])
 
-  >>> (np.dot(X, X.T) == fast_dot(X, X.T)).all()
+  >>> np.allclose(np.dot(X, X.T), fast_dot(X, X.T))
   True
 
 However this requires data to be exactly 2 dimensional and of the same type,

sklearn/decomposition/factor_analysis.py

@@ -19,7 +19,7 @@
 from ..base import BaseEstimator, TransformerMixin
 from ..externals.six.moves import xrange
 from ..utils import array2d, check_arrays
-from ..utils.extmath import fast_logdet
+from ..utils.extmath import fast_logdet, fast_dot
 
 
 class FactorAnalysis(BaseEstimator, TransformerMixin):
@@ -193,8 +193,8 @@ def transform(self, X):
 
         Wpsi = self.components_ / self.noise_variance_
         cov_z = linalg.inv(Ih + np.dot(Wpsi, self.components_.T))
-        tmp = np.dot(X_transformed, Wpsi.T)
-        X_transformed = np.dot(tmp, cov_z)
+        tmp = fast_dot(X_transformed, Wpsi.T)
+        X_transformed = fast_dot(tmp, cov_z)
 
         return X_transformed
 

sklearn/decomposition/fastica_.py

@@ -16,6 +16,7 @@
 from ..externals import six
 from ..externals.six import moves
 from ..utils import array2d, as_float_array, check_random_state, deprecated
+from ..utils.extmath import fast_dot
 
 __all__ = ['fastica', 'FastICA']
 
@@ -66,7 +67,7 @@ def _ica_def(X, tol, g, fun_args, max_iter, w_init):
         w /= np.sqrt((w ** 2).sum())
 
         for _ in moves.xrange(max_iter):
-            gwtx, g_wtx = g(np.dot(w.T, X), fun_args)
+            gwtx, g_wtx = g(fast_dot(w.T, X), fun_args)
 
             w1 = (X * gwtx).mean(axis=1) - g_wtx.mean() * w
 
@@ -94,12 +95,12 @@ def _ica_par(X, tol, g, fun_args, max_iter, w_init):
     del w_init
     p_ = float(X.shape[1])
     for ii in moves.xrange(max_iter):
-        gwtx, g_wtx = g(np.dot(W, X), fun_args)
-        W1 = _sym_decorrelation(np.dot(gwtx, X.T) / p_
+        gwtx, g_wtx = g(fast_dot(W, X), fun_args)
+        W1 = _sym_decorrelation(fast_dot(gwtx, X.T) / p_
                                 - g_wtx[:, np.newaxis] * W)
         del gwtx, g_wtx
         # builtin max, abs are faster than numpy counter parts.
-        lim = max(abs(abs(np.diag(np.dot(W1, W.T))) - 1))
+        lim = max(abs(abs(np.diag(fast_dot(W1, W.T))) - 1))
         W = W1
         if lim < tol:
             break
@@ -322,7 +323,7 @@ def g(x, fun_args):
 
     if whiten:
         if compute_sources:
-            S = np.dot(np.dot(W, K), X).T
+            S = fast_dot(fast_dot(W, K), X).T
         else:
             S = None
         if return_X_mean:
@@ -331,7 +332,7 @@ def g(x, fun_args):
             return K, W, S
     else:
         if compute_sources:
-            S = np.dot(W, X).T
+            S = fast_dot(W, X).T
         else:
             S = None
         if return_X_mean:
@@ -508,7 +509,7 @@ def transform(self, X, y=None, copy=True):
         if self.whiten:
             X -= self.mean_
 
-        return np.dot(X, self.components_.T)
+        return fast_dot(X, self.components_.T)
 
     @deprecated('To be removed in 0.16. Use the `mixing_` attribute.')
     def get_mixing_matrix(self):
@@ -543,7 +544,7 @@ def inverse_transform(self, X, copy=True):
         """
         if copy:
             X = X.copy()
-        X = np.dot(X, self.mixing_.T)
+        X = fast_dot(X, self.mixing_.T)
         if self.whiten:
             X += self.mean_
 

sklearn/decomposition/pca.py

@@ -4,6 +4,8 @@
 # Author: Alexandre Gramfort <[email protected]>
 #         Olivier Grisel <[email protected]>
 #         Mathieu Blondel <[email protected]>
+#         Denis A. Engemann <[email protected]>
+#
 # License: BSD 3 clause
 
 from math import log, sqrt
@@ -16,9 +18,8 @@
 from ..base import BaseEstimator, TransformerMixin
 from ..utils import array2d, check_random_state, as_float_array
 from ..utils import atleast2d_or_csr
-from ..utils.extmath import fast_logdet
-from ..utils.extmath import safe_sparse_dot
-from ..utils.extmath import randomized_svd
+from ..utils.extmath import fast_logdet, safe_sparse_dot, randomized_svd, \
+                            fast_dot
 
 
 def _assess_dimension_(spectrum, rank, n_samples, n_features):
@@ -303,7 +304,7 @@ def transform(self, X):
         X = array2d(X)
         if self.mean_ is not None:
             X = X - self.mean_
-        X_transformed = np.dot(X, self.components_.T)
+        X_transformed = fast_dot(X, self.components_.T)
         return X_transformed
 
     def inverse_transform(self, X):
@@ -325,7 +326,7 @@ def inverse_transform(self, X):
         If whitening is enabled, inverse_transform does not compute the
         exact inverse operation as transform.
         """
-        return np.dot(X, self.components_) + self.mean_
+        return fast_dot(X, self.components_) + self.mean_
 
 
 class ProbabilisticPCA(PCA):