LinearDiscriminantAnalysis factor loadings

docs/sources/CHANGELOG.md

@@ -15,13 +15,18 @@ The CHANGELOG for the current development version is available at
 ##### New Features
 
 - New `max_len` parameter for the frequent itemset generation via the `apriori` function to allow for early stopping. ([#270](https://github.com/rasbt/mlxtend/pull/270))
+- Added a `loadings_` attribute to `LinearDiscriminantAnalysis` to compute the factor loadings of the features on the components (discrimnants). ([#277](https://github.com/rasbt/mlxtend/pull/277))
+- Added a `tol` parameter to `LinearDiscriminantAnalysis` to threshold small eigenvalues, which are due to floating point imprecision, to zero for numerical stability. ([#277](https://github.com/rasbt/mlxtend/pull/277))
 
 ##### Changes
 
 - All feature index tuples in `SequentialFeatureSelector` or now in sorted order [#262](https://github.com/rasbt/mlxtend/pull/262)
 - The `SequentialFeatureSelector` now runs the continuation of the floating inclusion/exclusion as described in Novovicova & Kittler (1994). 
 Note that this didn't cause any difference in performance on any of the test scenarios but could lead to better performance in certain edge cases. 
 [#262](https://github.com/rasbt/mlxtend/pull/262)
+- Added "SVD solver" option to the `LinearDiscriminantAnalysis`. ([#277](https://github.com/rasbt/mlxtend/pull/277))
+- `LinearDiscriminantAnalysis` now uses NumPy's `cov` function instead of computing the scatter matrices manually to improve numerical stability. ([#277](https://github.com/rasbt/mlxtend/pull/277))
+
 
 
 ### Version 0.9.0 (2017-10-21)

docs/sources/user_guide/feature_extraction/PrincipalComponentAnalysis.ipynb

@@ -317,7 +317,9 @@
   {
    "cell_type": "code",
    "execution_count": 3,
-   "metadata": {},
+   "metadata": {
+    "collapsed": true
+   },
    "outputs": [],
    "source": [
     "from mlxtend.data import iris_data\n",

mlxtend/feature_extraction/linear_discriminant_analysis.py

@@ -19,6 +19,19 @@ class LinearDiscriminantAnalysis(_BaseModel):
     n_discriminants : int (default: None)
         The number of discrimants for transformation.
         Keeps the original dimensions of the dataset if `None`.
+        Note that the number of meaningful discriminants is
+        is max. n_classes - 1. In other words,
+        in LDA, the number of linear discriminants is at
+        most c-1, where c is the number of class labels,
+        since the in-between scatter matrix SB is
+        the sum of c matrices with rank 1 or less.
+        We can indeed see that we only have two nonzero eigenvalues
+    solver : str (default: 'eigen')
+        Method for performing the matrix decomposition.
+        {'eigen', 'svd'}
+    tol : float (default: 1-e8)
+        Tolerance value for thresholding small eigenvalues, which
+        are due to floating point imprecision, to zero.
 
     Attributes
     ----------
@@ -30,10 +43,17 @@ class LinearDiscriminantAnalysis(_BaseModel):
        Eigenvectors in sorted order.
 
     """
-    def __init__(self, n_discriminants=None):
+    def __init__(self, n_discriminants=None, solver='eigen', tol=1e-8):
+        valid_solver = {'eigen', 'svd'}
+        if solver not in valid_solver:
+            raise AttributeError('Must be in %s. Found %s'
+                                 % (valid_solver, solver))
+        self.solver = solver
+
         if n_discriminants is not None and n_discriminants < 1:
             raise AttributeError('n_discriminants must be > 1 or None')
         self.n_discriminants = n_discriminants
+        self.tol = tol
 
     def fit(self, X, y, n_classes=None):
         """ Fit the LDA model with X.
@@ -63,6 +83,7 @@ def fit(self, X, y, n_classes=None):
 
     def _fit(self, X, y, n_classes=None):
 
+        n_samples = X.shape[0]
         if self.n_discriminants is None or self.n_discriminants > X.shape[1]:
             n_discriminants = X.shape[1]
         else:
@@ -82,11 +103,17 @@ def _fit(self, X, y, n_classes=None):
         between_scatter = self._between_scatter(X=X,
                                                 y=y,
                                                 mean_vectors=mean_vecs)
-        self.e_vals_, self.e_vecs_ = self._eigendecom(
-            within_scatter=within_scatter, between_scatter=between_scatter)
+        self.e_vals_, self.e_vecs_ = self._decomposition(
+            within_scatter=within_scatter,
+            between_scatter=between_scatter,
+            n_samples=n_samples)
+
+        self.e_vals_ = self.e_vals_.copy()
+        self.e_vals_[abs(self.e_vals_) < self.tol] = 0.0
         self.w_ = self._projection_matrix(eig_vals=self.e_vals_,
                                           eig_vecs=self.e_vecs_,
                                           n_discriminants=n_discriminants)
+        self.loadings_ = self._loadings()
         return self
 
     def transform(self, X):
@@ -118,11 +145,13 @@ def _mean_vectors(self, X, y, n_classes):
     def _within_scatter(self, X, y, n_classes, mean_vectors):
         S_W = np.zeros((X.shape[1], X.shape[1]))
         for cl, mv in zip(range(n_classes), mean_vectors):
-            class_sc_mat = np.zeros((X.shape[1], X.shape[1]))
-            for row in X[y == cl]:
-                row, mv = row.reshape(X.shape[1], 1), mv.reshape(X.shape[1], 1)
-                class_sc_mat += (row - mv).dot((row - mv).T)
-            S_W += class_sc_mat
+            class_sc_mat = np.cov((X[y == cl] - mv).T)
+            # class_sc_mat = np.zeros((X.shape[1], X.shape[1]))
+            # for row in X[y == cl]:
+            #    row, mv = row.reshape(X.shape[1], 1),
+            #                          mv.reshape(X.shape[1], 1)
+            #    class_sc_mat += (row - mv).dot((row - mv).T)
+            S_W += y[y == cl].shape[0] * class_sc_mat
         return S_W
 
     def _between_scatter(self, X, y, mean_vectors):
@@ -136,9 +165,15 @@ def _between_scatter(self, X, y, mean_vectors):
                 (mean_vec - overall_mean).T)
         return S_B
 
-    def _eigendecom(self, within_scatter, between_scatter):
-        e_vals, e_vecs = np.linalg.eig(np.linalg.inv(within_scatter).dot(
-            between_scatter))
+    def _decomposition(self, within_scatter, between_scatter, n_samples):
+        combined_scatter = np.linalg.inv(within_scatter).dot(
+            between_scatter)
+        if self.solver == 'eigen':
+            e_vals, e_vecs = np.linalg.eig(combined_scatter)
+        elif self.solver == 'svd':
+            u, s, v = np.linalg.svd(combined_scatter)
+            e_vecs, e_vals = u, s
+            e_vals = e_vals ** 2 / n_samples
         sort_idx = np.argsort(e_vals)[::-1]
         e_vals, e_vecs = e_vals[sort_idx], e_vecs[sort_idx]
         return e_vals, e_vecs
@@ -147,3 +182,9 @@ def _projection_matrix(self, eig_vals, eig_vecs, n_discriminants):
         matrix_w = np.vstack([eig_vecs[:, i] for
                               i in range(n_discriminants)]).T
         return matrix_w
+
+    def _loadings(self):
+        """Compute factor loadings"""
+
+        return (self.e_vecs_ *
+                np.sqrt(np.abs(self.e_vals_)))

mlxtend/feature_extraction/tests/test_linear_discriminant_analysis.py

@@ -18,7 +18,7 @@
 def test_default_components():
     lda = LDA()
     lda.fit(X, y)
-    res = lda.fit(X).transform(X)
+    res = lda.fit(X, y).transform(X)
     assert res.shape[1] == 4
 
 
@@ -30,18 +30,28 @@ def test_default_2components():
 
 
 @raises(AttributeError)
-def test_default_components():
+def test_default_components_0():
     lda = LDA(n_discriminants=0)
     lda.fit(X, y)
-    res = lda.fit(X).transform(X)
 
 
-def test_evals():
+def test_evals_eigen():
     lda = LDA(n_discriminants=2)
-    res = lda.fit(X, y).transform(X)
+    lda.fit(X, y).transform(X)
     np.set_printoptions(suppress=True)
     print('%s' % lda.e_vals_)
-    assert_almost_equal(lda.e_vals_, [20.90, 0.14, 0.0, 0.0], decimal=2)
+    assert_almost_equal(lda.e_vals_, [20.49, 0.14, 0.0, 0.0], decimal=2)
+
+
+def test_evecs_eigen_vs_svd():
+
+    lda = LDA(n_discriminants=2)
+    lda.fit(X, y).transform(X)
+    eigen_vecs = lda.e_vecs_
+    lda = LDA(n_discriminants=2, solver='svd')
+    lda.fit(X, y).transform(X)
+    assert_almost_equal(lda.e_vecs_[:, 0],
+                        eigen_vecs[:, 0], decimal=2)
 
 
 @raises(ValueError)
@@ -54,4 +64,16 @@ def test_fail_array_fit():
 def test_fail_array_transform():
     lda = LDA()
     lda.fit(X, y)
-    exp = lda.transform(X[1])
+    lda.transform(X[1])
+
+
+def test_loadings():
+
+    expect = np.abs(np.array([[0.7, 0., 0., 0.],
+                              [0.7, 0.1, 0., 0.],
+                              [3.9, 0.3, 0., 0.],
+                              [2.1, 0.2, 0., 0.]]))
+
+    lda = LDA(n_discriminants=2)
+    lda.fit(X, y)
+    assert_almost_equal(np.abs(lda.loadings_), expect, decimal=1)