# PCA for complex-valued data

I'm quite shocked for encountering this error on PCA from sklearn

ValueError: Complex data not supported

After trying to fit complex-valued data. Is this just unimplemented thing? Should I just go ahead and do it 'manually' with SVD or is their a catch for complex-values?

Apparently this functionality is left out intentionally, see here. I'm afraid you have to use SVD, but that should be fairly straightforward:

def pca(X):
mean = X.mean(axis=0)
center = X - mean
_, stds, pcs = np.linalg.svd(center/np.sqrt(X.shape[0]))

return stds**2, pcs

• Thanks. Although there's no need for sorting, it is sorted already.
– Alex
Jun 12, 2020 at 3:38
• Oh yeah, that's true, thanks! I was thinking of np.linalg.eig. I'll fix the answer. Jun 12, 2020 at 10:42

While I'm late in the game, I felt that an deeper insight to this question could be provided for those who would like to use the above mentioned algorithms (matthiaw91 and Alex) on the data. Especially since sklearn could easily implement their algorithms but chose not to, specifying that it would require substantially more work.

As explained very neatly by amoeba here: relationship between SVD and PCA The relationship between SVD and PCA lies in the Covariance of the matrix $$X$$. For real multivariate cases, lets assume $$X$$ is from a normal distribution, i.e. $$X \sim N(\mu, \Gamma)$$ where $$\mu = E[X]$$ is the mean and $$\Gamma = E[(X-\mu)(X-\mu)^T]$$ is the covariance. This means that if we center $$X$$, $$X - \mu => X \sim N(0, \Gamma)$$, then $$X$$ can be completely described by the covariance $$\Gamma$$. This is one reason why we decompose the covariance matrix in PCA to gain the load vectors.

Lets view the complex cases. when we have a complex multivariate vector or matrix, we say that $$X$$ is complex normal (wiki) if $$X \sim CN(\mu, \Gamma, C)$$ where yet again $$\mu$$ is the mean, $$\Gamma = E[(X-\mu)(X-\mu)^H]$$ is the covariance matrix (note the hermitian transpose H instead of regular transpose T) and $$C = E[(X-\mu)(X-\mu)^T]$$ is the pseudo-covariance matrix. Yet again we center and get $$X \sim CN(0, \Gamma, C)$$. Here though, we get that we both need the covariance and the pseudo-covariance matrix to accurately describe the data. The svd, which relies on the covariance matrix may thus not be accurate if not also $$C = 0$$. When $$C = 0$$ we say that $$X$$ is circular-symmetric, which is further described in the wiki.

In the case of non-circular-symmetric data, i.e. $$C \not = 0$$ for instance widely linear transformation can be adapted to enhance the data analysis. A good example and implementation I found is Component Analysis of complex-valued data for machine learning and computer vision tasks by A. D. Papaioannou which is free to download. In this case we let

$$\underline X = [X, X^*]$$ as we can show that $$cov(\underline X) = E[\underline X \underline X^H] = [[\Gamma, C], [C^*, \Gamma^*]]$$ and thus accurately describes $$X$$. Though, I will be quick to say that by simply plugging $$\underline X$$ into the svd is not the way to go, (or at least I don't think so) but the widely linear estimation are needed (see the article). But I have not done the calculations of svd using the augmented form ($$\underline X$$), so maybe it actually is the same, check for yourself.

Complex numbers are more 'complex' to deal with, as one would expect.

My implementation exactly mimicks the original PCA so any existing code that deals with PCA would work seamlessly.

class ComplexPCA:
def __init__(self, n_components):
self.n_components = n_components
self.u = self.s = self.components_ = None
self.mean_ = None

@property
def explained_variance_ratio_(self):
return self.s

def fit(self, matrix, use_gpu=False):
self.mean_ = matrix.mean(axis=0)
if use_gpu:
import tensorflow as tf  # torch doesn't handle complex values.
tensor = tf.convert_to_tensor(matrix)
u, s, vh = tf.linalg.svd(tensor, full_matrices=False)  # full=False ==> num_pc = min(N, M)
# It would be faster if the SVD was truncated to only n_components instead of min(M, N)
else:
_, self.s, vh = np.linalg.svd(matrix, full_matrices=False)  # full=False ==> num_pc = min(N, M)
# It would be faster if the SVD was truncated to only n_components instead of min(M, N)
self.components_ = vh  # already conjugated.
# Leave those components as rows of matrix so that it is compatible with Sklearn PCA.

def transform(self, matrix):
data = matrix - self.mean_
result = data @ self.components_.T
return result

def inverse_transform(self, matrix):
result = matrix @ np.conj(self.components_)
return self.mean_ + result


• Caveat: GPU doesn't save you any time! SVD is not a paralleizable aglo.
– Alex
Jun 12, 2020 at 7:46