ENH : add get_precision method with matrix inverse lemma to FactorAnalysis + use precision in score

Author: agramfort

sklearn/decomposition/factor_analysis.py

@@ -267,6 +267,33 @@ def get_covariance(self):
         cov.flat[::len(cov) + 1] += self.noise_variance_  # modify diag inplace
         return cov
 
+    def get_precision(self):
+        """Compute data precision matrix with the FactorAnalysis model.
+
+        Returns
+        -------
+        precision : array, shape=(n_features, n_features)
+            Estimated precision of data.
+        """
+        n_features = self.components_.shape[1]
+
+        # handle corner cases first
+        if self.n_components == 0:
+            return np.diag(1. / self.noise_variance_)
+        if self.n_components == n_features:
+            return linalg.inv(self.get_covariance())
+
+        # Get precision using matrix inversion lemma
+        components_ = self.components_
+        precision = np.dot(components_ / self.noise_variance_, components_.T)
+        precision.flat[::len(precision) + 1] += 1.
+        precision = np.dot(components_.T,
+                           np.dot(linalg.inv(precision), components_))
+        precision /=  self.noise_variance_[:, np.newaxis]
+        precision /=  -self.noise_variance_[np.newaxis, :]
+        precision.flat[::len(precision) + 1] += 1. / self.noise_variance_
+        return precision
+
     def score(self, X, y=None):
         """Compute score of X under FactorAnalysis model.
 
@@ -281,10 +308,10 @@ def score(self, X, y=None):
             log-likelihood of each row of X under the current model
         """
         Xr = X - self.mean_
-        cov = self.get_covariance()
+        precision = self.get_precision()
         n_features = X.shape[1]
         log_like = np.zeros(X.shape[0])
-        self.precision_ = linalg.inv(cov)
-        log_like = -.5 * (Xr * (fast_dot(Xr, self.precision_))).sum(axis=1)
-        log_like -= .5 * (fast_logdet(cov) + n_features * log(2. * np.pi))
+        self.precision_ = precision  # should not store it I guess...
+        log_like = -.5 * (Xr * (np.dot(Xr, precision))).sum(axis=1)
+        log_like -= .5 * (-fast_logdet(precision) + n_features * log(2. * np.pi))
         return log_like

sklearn/decomposition/tests/test_factor_analysis.py

@@ -11,8 +11,8 @@
 from sklearn.utils.testing import assert_less
 from sklearn.utils.testing import assert_raises
 from sklearn.utils.testing import assert_almost_equal
+from sklearn.utils.testing import assert_array_almost_equal
 from sklearn.utils import ConvergenceWarning
-
 from sklearn.decomposition import FactorAnalysis
 
 
@@ -62,6 +62,7 @@ def test_factor_analysis():
                             noise_variance_init=np.ones(n_features))
         assert_raises(ValueError, fa.fit, X[:, :2])
 
+
     f = lambda x, y: np.abs(getattr(x, y))  # sign will not be equal
     fa1, fa2 = fas
     for attr in ['loglike_', 'components_', 'noise_variance_']:
@@ -76,3 +77,26 @@ def test_factor_analysis():
         warnings.simplefilter('always', DeprecationWarning)
         FactorAnalysis(verbose=1)
         assert_true(w[-1].category == DeprecationWarning)
+
+    fa2 = FactorAnalysis(n_components=n_components,
+                         noise_variance_init=np.ones(n_features))
+    assert_raises(ValueError, fa2.fit, X[:, :2])
+
+    # Test get_covariance and get_precision with n_components < n_features
+    cov = fa.get_covariance()
+    precision = fa.get_precision()
+    assert_array_almost_equal(np.dot(cov, precision), np.eye(X.shape[1]), 12)
+
+    # Test get_covariance and get_precision with n_components == n_features
+    fa.n_components = n_features
+    fa.fit(X)
+    cov = fa.get_covariance()
+    precision = fa.get_precision()
+    assert_array_almost_equal(np.dot(cov, precision), np.eye(X.shape[1]), 12)
+
+    # Test get_covariance and get_precision with n_components == 0
+    fa.n_components = 0
+    fa.fit(X)
+    cov = fa.get_covariance()
+    precision = fa.get_precision()
+    assert_array_almost_equal(np.dot(cov, precision), np.eye(X.shape[1]), 12)