Reported by marcus on 20 Apr 43350065 14:46 UTC
I've tested the Kernel PCA code and got some unexpected results.
So I looked into the code (kernel_pca_impl.hpp). I suppose it's not the right way to compute the covariance matrix (arma::mat transData = ccov(data);) and then compute the Kernel Matrix. Maybe I got something wrong. But consider the data points vec(x) and vec(y) in the input space I = R^n. When we map the data non linearly into a feature space F by
I = R^n
Phi: R^n -> F, x -> X
there is no covariance on which to perform eigendecomposition explicitly as we would in linear PCA.
I've written some code which provides the expected results. Maybe someone can explain if I made a wrong assumption.
using namespace mlpack;
using namespace mlpack::kpca;
using namespace mlpack::kernel;
int main(int argc, char** argv)
arma::mat data, transformedData, eigvec;
// Load the data.
data::Load("circle.txt", data, true);
// Using the Gaussian Kernel to construct the Kernel Matrix.
arma::mat kernelMat = GetKernelMatrix(kernel, trans(data));
// For PCA the data has to be centered, even if the data is centered.
// It is nit guarantee the data when mapped is also centered. Since
// we actually never work in the feature space we cannot center the data.
// Since centered data is required to perform an effective principal
// component analysis we perform a pseudo center method using the Kernel Matrix.
arma::mat oneMat = arma::ones<arma::mat>(kernelMat.n_rows, kernelMat.n_cols);
arma::mat kernelMatCenter = kernelMat - oneMat*kernelMat - kernelMat*oneMat + oneMat*kernelMat*oneMat;
// Compute eigenvectors and the corresponding eigenvalues.
arma::eig_sym(eigVal, eigvec, kernelMatCenter);
// The eigenvectors and the corresponding eigenvalues are
// already sorted but in the wrong order.
// Since descend is required, we reverse the eigenvectors
// and the corresponding eigenvalues.
// To avoid temporary matrices we use swap.
int n_eigVal = eigVal.n_elem;
for(int i = 0; i < floor(n_eigVal / 2.0); i++)
eigVal.swap_rows(i, (n_eigVal - 1) - i);
eigvec = arma::fliplr(eigvec);
// Dimension of output data.
size_t dim = 2;
// Projecting the data in lower dimensions.
transformedData = eigvec.submat(0, 0, eigvec.n_rows-1, dim-1).t() * kernelMatCenter.t();
std::cout << transformedData.t() << std::endl;
I've attached the data set and the plots of the results.
Thanks and regards,
Commented by rcurtin on 4 Jan 43354293 18:37 UTC
You are right, something seems very odd in the first line of KernelPCA::Apply():
arma::mat transData = ccov(data);
That's not how you transpose data...
Intuitively it seems like your implementation is correct, but I need to take a look at the test in kernel_pca_test.cpp and see if it is also flawed. If it is, I will devise a new test, and I suspect that your implementation will pass the test and the existing one will not -- in which case, I'll put your implementation in place.
Thanks for reporting this.
Commented by rcurtin on 23 Dec 43359018 22:35 UTC
r15050 should fix this. This also helps #204 towards completion. There is a much more robust test in place now, which basically generates the circle dataset (like you attached), then uses Kernel PCA with the Gaussian Kernel to map it to one dimension and ensures that the three "classes" are linearly separable.
I took your implementation and adapted it into kernel_pca_impl.hpp, so unless you ask me to do otherwise I've added your name to the authors list in that file and core.hpp and on the website.
Thanks again for pointing this out.
Commented by marcus on 24 Jul 43361472 23:16 UTC
Thank you. It was great to help.