In the article Relative Information Loss in the PCA, the authors make, at some point (in the introductory section), the following statement:
In case the orthogonal matrix is not known a priori, but has to be estimated from a set of input data vectors collected in the matrix $\underline{X}$, the PCA becomes a nonlinear operation:
$$\underline{Y} = \underline{w}(\underline{X})\underline{X}$$
Here, $\underline{w}$ is a matrix-valued function which computes the orthogonal matrix required for rotating the data (e.g., using the QR algorithm).
This statement contrasts with most statements about PCA, which is regarded as a linear transformation.
I designed a toy experiment to check the linearity (additivity property): $f(a + b) = f(a) + f(b)$.
import numpy
from sklearn.decomposition import PCA
if __name__ == '__main__':
numpy.random.seed(42)
m = 100
d = 3
X = numpy.random.normal(size = (m, d))
# Center data
X -= numpy.mean(X, axis = 0)
pca = PCA(n_components = d)
pca.fit(X)
Y = pca.transform(X)
# Check linearity, pca(a + b) = pca(a) + pca(b)
for i in range(0, m):
for j in range(0, m):
d = pca.transform([X[i] + X[j]]) - (Y[i] + Y[j])
assert numpy.allclose(d, numpy.array([0.0, 0.0, 0.0]))
The expression $f(a + b) - (f(a) + f(b))$, where $f = \mathrm{PCA}$, seems to be the zero vector, thus I assume the transform (PCA) is linear.
What am I missing then, that PCA is considered non-linear when the orthogonal matrix (the matrix of the principal components) is estimated from $X$ (see the quote above)?