Repository: pykale/pykale
Issue: [Bug] Complex number in MPCA eigenvalues/eigenvectors

# 🐛 Complex number in MPCA eigenvalues/eigenvectors  

<!-- A clear and concise description of what the bug is. -->


## To reproduce

In some cases, the eigen decompostion in MPCA can return complex eigenvalues and eignevectors.

```
>>> eig_values
Out[1]: 
array([ 2.03946400e+004+0.00000000e+00j,  8.53028325e+003+0.00000000e+00j,
        6.77832566e+003+0.00000000e+00j,  5.56850912e+003+0.00000000e+00j,
        4.14053054e+003+0.00000000e+00j,  3.85073231e+003+0.00000000e+00j,
        2.99306387e+003+0.00000000e+00j,  2.51583008e+003+0.00000000e+00j,
        2.07546552e+003+0.00000000e+00j,  1.72131535e+003+0.00000000e+00j,
        1.50710446e+003+0.00000000e+00j,  1.24574081e+003+0.00000000e+00j,
        1.06374939e+003+0.00000000e+00j,  8.82491750e+002+0.00000000e+00j,
        7.77560753e+002+0.00000000e+00j,  1.01041037e+002+0.00000000e+00j,
        6.41473082e+002+0.00000000e+00j,  5.66840523e+002+0.00000000e+00j,
        4.61197948e+002+0.00000000e+00j,  3.78630061e+002+0.00000000e+00j,
        2.13165912e+002+0.00000000e+00j,  2.55552782e+002+0.00000000e+00j,
        6.06324756e-013+0.00000000e+00j,  9.97100169e-015+0.00000000e+00j,
       -3.31739883e-016+0.00000000e+00j, -4.13037170e-022+0.00000000e+00j,
        7.86664595e-028+8.39380061e-27j,  7.86664595e-028-8.39380061e-27j,
        3.67760555e-027+4.43588239e-27j,  3.67760555e-027-4.43588239e-27j,
        5.10824776e-027+0.00000000e+00j, -4.05596311e-027+2.34425407e-27j,
       -4.05596311e-027-2.34425407e-27j, -1.58477788e-027+3.16041573e-27j,
       -1.58477788e-027-3.16041573e-27j, -3.58807233e-027+0.00000000e+00j,
        2.27118464e-027+1.39384460e-27j,  2.27118464e-027-1.39384460e-27j,
        9.98017692e-028+6.54044843e-28j,  9.98017692e-028-6.54044843e-28j,
       -1.27894122e-027+7.45212435e-28j, -1.27894122e-027-7.45212435e-28j,
       -8.19860299e-028+0.00000000e+00j,  4.44355796e-029+0.00000000e+00j,
       -4.33625284e-045+1.20085059e-44j, -4.33625284e-045-1.20085059e-44j,
       -1.96760884e-062+0.00000000e+00j, -3.83573351e-060+0.00000000e+00j,
        1.01157335e-076+0.00000000e+00j,  1.83716727e-091+0.00000000e+00j,
       -9.60847167e-092+0.00000000e+00j, -2.32451013e-092+0.00000000e+00j,
        7.71193314e-106+0.00000000e+00j, -1.86756379e-107+0.00000000e+00j,
       -4.61915882e-122+0.00000000e+00j, -1.13689708e-121+0.00000000e+00j,
        1.84293425e-138+0.00000000e+00j, -1.44155321e-137+0.00000000e+00j,
       -1.29151358e-152+0.00000000e+00j,  1.94857963e-155+0.00000000e+00j,
        4.47153100e-183+0.00000000e+00j,  1.06315569e-168+0.00000000e+00j,
        0.00000000e+000+0.00000000e+00j, -1.54240939e-199+0.00000000e+00j])
>>> eig_values.dtype
Out[2]: dtype('complex128')
```

## Expected Behaviour

In MPCA, we decompose the sum of covariance matrices for each mode. Because covaricance matrix is positive definite real symmetric, and therefore the sum of covariance matrices is a positive definite real symmetric matrix. The eigenvalues should be real and non-negative. This may be an python numerical error or the a problem of numpy `eig` function (ref: https://stackoverflow.com/questions/8765310/scipy-linalg-eig-return-complex-eigenvalues-for-covariance-matrix).

**So my question is, can we simply take the real part and ignore the complex part? Alternatively, can we use SVD to obatin eigenvalues and eigenvectors. I have tried the SVD implemented in Scipy for the same case, and no complex numbers returned.**

SVD singular values for the same input
```
array([2.03946400e+04, 8.53028325e+03, 6.77832566e+03, 5.56850912e+03,
       4.14053054e+03, 3.85073231e+03, 2.99306387e+03, 2.51583008e+03,
       2.07546552e+03, 1.72131535e+03, 1.50710446e+03, 1.24574081e+03,
       1.06374939e+03, 8.82491750e+02, 7.77560753e+02, 6.41473082e+02,
       5.66840523e+02, 4.61197948e+02, 3.78630061e+02, 2.55552782e+02,
       2.13165912e+02, 1.01041037e+02, 9.60507019e-12, 4.23097888e-12,
       3.72814677e-12, 2.92081445e-12, 2.24098733e-12, 1.39518126e-12,
       1.34449963e-12, 1.34449963e-12, 1.34449963e-12, 1.34449963e-12,
       1.34449963e-12, 1.34449963e-12, 1.34449963e-12, 1.34449963e-12,
       1.34449963e-12, 1.34449963e-12, 1.34449963e-12, 1.34449963e-12,
       1.34449963e-12, 1.34449963e-12, 1.34449963e-12, 1.34449963e-12,
       1.34449963e-12, 1.34449963e-12, 1.34449963e-12, 1.34449963e-12,
       1.34449963e-12, 1.34449963e-12, 1.34449963e-12, 1.34449963e-12,
       1.34449963e-12, 1.34449963e-12, 1.34449963e-12, 1.34449963e-12,
       1.34449963e-12, 1.34449963e-12, 1.34449963e-12, 1.34449963e-12,
       1.34449963e-12, 8.90409683e-13, 4.51888166e-13, 4.00569300e-13])
```
<!-- A clear and concise description of what you expected to happen. -->

## Environment
Windows 10 OS
```
[pip3] numpy==1.19.2
[pip3] numpydoc==1.1.0
[pip3] pytorch-lightning==0.8.5
[pip3] torch==1.7.1
[pip3] torchaudio==0.7.2
[pip3] torchvision==0.8.2
[conda] blas                      1.0                         mkl
[conda] cudatoolkit               10.1.243             h74a9793_0
[conda] mkl                       2020.2                      256
[conda] mkl-service               2.3.0            py37h196d8e1_0
[conda] mkl_fft                   1.2.1            py37h46781fe_0
[conda] mkl_random                1.1.1            py37h47e9c7a_0
[conda] numpy                     1.19.2           py37hadc3359_0
[conda] numpy-base                1.19.2           py37ha3acd2a_0
[conda] numpydoc                  1.1.0              pyhd3eb1b0_1
[conda] pytorch                   1.7.1           py3.7_cuda101_cudnn7_0    pytorch
[conda] pytorch-lightning         0.8.5                    pypi_0    pypi
[conda] torchaudio                0.7.2                      py37    pytorch
[conda] torchvision               0.8.2                py37_cu101    pytorch
```


Changed files:
setup.py, kale/embed/mpca.py, kale/embed/__init__.py

diff --git a/kale/embed/__init__.py b/kale/embed/__init__.py
index 662210f..929e5c6 100644
--- a/kale/embed/__init__.py
+++ b/kale/embed/__init__.py
@@ -10,9 +10,7 @@
 #     SmallCNNFeature,
 # )
 # from .linformer import LinearMultiheadAttention, LinearTransformerEncoderLayer
-# from ._mpca import MPCA
 # from .mpca import MPCA
-# from ._mpca_torch import MPCA
 # from .positional_encoding import PositionalEncoding
 # from .video_i3d import i3d, InceptionI3d, InceptionModule, MaxPool3dSamePadding, Unit3D
 # from .video_res3d import (
diff --git a/kale/embed/mpca.py b/kale/embed/mpca.py
index 75de815..4a27e24 100644
--- a/kale/embed/mpca.py
+++ b/kale/embed/mpca.py
@@ -14,6 +14,7 @@ import logging
 import warnings
 
 import numpy as np
+from scipy import linalg
 from sklearn.base import BaseEstimator, TransformerMixin
 
 # import tensorly as tl
@@ -74,7 +75,7 @@ class MPCA(BaseEstimator, TransformerMixin):
         n_components (int): Number of components to keep. Applies only when return_vector=True. Defaults to None.
 
     Attributes:
-        proj_mats (list of array-like): A list of transposed projection matrices, shapes (P_1, I_1), ...,
+        proj_mats (list of arrays): A list of transposed projection matrices, shapes (P_1, I_1), ...,
             (P_N, I_N), where P_1, ..., P_N are output tensor shape for each sample.
         idx_order (array-like): The ordering index of projected (and vectorised) features in decreasing variance.
         mean_ (array-like): Per-feature empirical mean, estimated from the training set, shape (I_1, I_2, ..., I_N).
@@ -82,7 +83,7 @@ class MPCA(BaseEstimator, TransformerMixin):
         shape_out (tuple): Output tensor shapes, i.e. (P_1, P_2, ..., P_N).
     Examples:
         >>> import numpy as np
-        >>> from kale.embed import MPCA
+        >>> from kale.embed.mpca import MPCA
         >>> x = np.random.random((40, 20, 25, 20))
         >>> x.shape
         (40, 20, 25, 20)
@@ -114,10 +115,10 @@ class MPCA(BaseEstimator, TransformerMixin):
         self.n_components = n_components
 
     def fit(self, x, y=None):
-        """Fit the model with input training data X.
+        """Fit the model with input training data x.
 
         Args
-            X (array-like tensor): Input data, shape (n_samples, I_1, I_2, ..., I_N), where n_samples is the number of
+            x (array-like tensor): Input data, shape (n_samples, I_1, I_2, ..., I_N), where n_samples is the number of
                 samples, I_1, I_2, ..., I_N are the dimensions of corresponding mode (1, 2, ..., N), respectively.
             y (None): Ignored variable.
 
@@ -131,7 +132,6 @@ class MPCA(BaseEstimator, TransformerMixin):
         """Solve MPCA"""
 
         shape_ = x.shape  # shape of input data
-        n_samples = shape_[0]  # number of samples
         n_dims = x.ndim
 
         self.shape_in = shape_[1:]
@@ -139,58 +139,50 @@ class MPCA(BaseEstimator, TransformerMixin):
         x = x - self.mean_
 
         # init
-        phi = dict()
         shape_out = ()
         proj_mats = []
 
+        # get the output tensor shape based on the cumulative distribution of eigen values for each mode
         for i in range(1, n_dims):
-            for j in range(n_samples):
-                # unfold the j_th tensor along the i_th mode
-                x_ji = unfold(x[j, :], mode=i - 1)
-                if i not in phi:
-                    phi[i] = 0
-                phi[i] = phi[i] + np.dot(x_ji, x_ji.T)
-
-        # get the output tensor shape based on the cumulative distribution of eigenvalues for each mode
-        for i in range(1, n_dims):
-            eig_values, eig_vectors = np.linalg.eig(phi[i])
+            mode_data_mat = unfold(x, mode=i)
+            singular_vec_left, singular_val, singular_vec_right = linalg.svd(mode_data_mat, full_matrices=False)
+            eig_values = np.square(singular_val)
             idx_sorted = (-1 * eig_values).argsort()
             cum = eig_values[idx_sorted]
-            var_tot = np.sum(cum)
+            tot_var = np.sum(cum)
 
             for j in range(1, cum.shape[0] + 1):
-                if np.sum(cum[:j]) / var_tot > self.var_ratio:
+                if np.sum(cum[:j]) / tot_var > self.var_ratio:
                     shape_out += (j,)
                     break
-            proj_mats.append(eig_vectors[:, idx_sorted][:, : shape_out[i - 1]].T)
+            proj_mats.append(singular_vec_left[:, idx_sorted][:, : shape_out[i - 1]].T)
 
         # set n_components to the maximum n_features if it is None
         if self.n_components is None:
             self.n_components = int(np.prod(shape_out))
 
         for i_iter in range(self.max_iter):
-            phi = dict()
             for i in range(1, n_dims):  # ith mode
-                if i not in phi:
-                    phi[i] = 0
-                for j in range(n_samples):
-                    xj = multi_mode_dot(
-                        x[j, :],  # jth tensor/sample
-                        [proj_mats[m] for m in range(n_dims - 1) if m != i - 1],
-                        modes=[m for m in range(n_dims - 1) if m != i - 1],
-                    )
-                    xj_unfold = unfold(xj, i - 1)
-                    phi[i] = np.dot(xj_unfold, xj_unfold.T) + phi[i]
-
-                eig_values, eig_vectors = np.linalg.eig(phi[i])
+                x_projected = multi_mode_dot(
+                    x,
+                    [proj_mats[m] for m in range(n_dims - 1) if m != i - 1],
+                    modes=[m for m in range(1, n_dims) if m != i]
+                )
+                mode_data_mat = unfold(x_projected, i)
+
+                singular_vec_left, singular_val, singular_vec_right = linalg.svd(mode_data_mat, full_matrices=False)
+                eig_values = np.square(singular_val)
                 idx_sorted = (-1 * eig_values).argsort()
-                proj_mats[i - 1] = (eig_vectors[:, idx_sorted][:, : shape_out[i - 1]]).T
+                proj_mats[i - 1] = (singular_vec_left[:, idx_sorted][:, : shape_out[i - 1]]).T
 
-        x_projected = multi_mode_dot(x, proj_mats,
-                                     modes=[m for m in range(1, n_dims)])
+        x_projected = multi_mode_dot(
+            x,
+            proj_mats,
+            modes=[m for m in range(1, n_dims)]
+        )
         x_proj_unfold = unfold(x_projected, mode=0)  # unfold the tensor projection to shape (n_samples, n_features)
         # x_proj_cov = np.diag(np.dot(x_proj_unfold.T, x_proj_unfold))  # covariance of unfolded features
-        x_proj_cov = np.sum(np.multiply(x_proj_unfold, x_proj_unfold), axis=1)  # memory saving computing covariance
+        x_proj_cov = np.sum(np.multiply(x_proj_unfold.T, x_proj_unfold.T), axis=1)  # memory saving computing covariance
         idx_order = (-1 * x_proj_cov).argsort()
 
         self.proj_mats = proj_mats
@@ -201,7 +193,7 @@ class MPCA(BaseEstimator, TransformerMixin):
         return self
 
     def transform(self, x):
-        """Perform dimension reduction on X
+        """Perform dimension reduction on x
 
         Args:
             x (array-like tensor): Data to perform dimension reduction, shape (n_samples, I_1, I_2, ..., I_N).
@@ -213,6 +205,9 @@ class MPCA(BaseEstimator, TransformerMixin):
                 (n_samples, P_1 * P_2 * ... * P_N) if self.n_components is None, and shape (n_samples, n_components)
                 if self.n_component is a valid integer.
         """
+        # reshape x to shape (1, I_1, I_2, ..., I_N) if x in shape (I_1, I_2, ..., I_N), i.e. n_samples = 1
+        if x.ndim == self.n_dims - 1:
+            x = x.reshape((1, ) + x.shape)
         _check_tensor_dim_shape(x, self.n_dims, self.shape_in)
         x = x - self.mean_
 
@@ -247,7 +242,7 @@ class MPCA(BaseEstimator, TransformerMixin):
             array-like tensor:
                 Reconstructed tensor in original shape, shape (n_samples, I_1, I_2, ..., I_N)
         """
-        # reshape X to tensor in shape (n_samples, self.shape_out) if X is vectorised
+        # reshape x to tensor in shape (n_samples, self.shape_out) if x has been unfolded
         if x.ndim <= 2:
             if x.ndim == 1:
                 # reshape x to a 2D matrix (1, n_components) if x in shape (n_components,)
diff --git a/setup.py b/setup.py
index 415f224..1eb72d9 100644
--- a/setup.py
+++ b/setup.py
@@ -13,11 +13,12 @@ from setuptools import find_packages, setup
 requirements = [
     'numpy>=1.18.0',
     'pytorch-lightning',
-    'scikit-image',
-    'scikit-learn',
-    'tensorly',
+    'scikit-image>=0.16.2',
+    'scikit-learn>=0.23.2',
+    'tensorly>=0.5.1',
     'torch>=1.7.0',
     'torchvision',  # >=0.8.1
+    'scipy>=1.5.4',
 ]
 
 # Additional dependencies for examples/tutorials