smile.projection.KPCA Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of hivemall-all
The newest version!
/*******************************************************************************
 * Copyright (c) 2010 Haifeng Li
 *   
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *  
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *******************************************************************************/
package smile.projection;

import smile.math.Math;
import smile.math.kernel.MercerKernel;
import smile.math.matrix.EigenValueDecomposition;

/**
 * Kernel principal component analysis. Kernel PCA is an extension of
 * principal component analysis (PCA) using techniques of kernel methods.
 * Using a kernel, the originally linear operations of PCA are done in a
 * reproducing kernel Hilbert space with a non-linear mapping.
 * 
 * In practice, a large data set leads to a large Kernel/Gram matrix K, and
 * storing K may become a problem. One way to deal with this is to perform
 * clustering on your large dataset, and populate the kernel with the means
 * of those clusters. Since even this method may yield a relatively large K,
 * it is common to compute only the top P eigenvalues and eigenvectors of K.
 * 

 * Kernel PCA with an isotropic kernel function is closely related to metric MDS.
 * Carrying out metric MDS on the kernel matrix K produces an equivalent configuration
 * of points as the distance (2(1 - K(x_i, x_j)))^1/2
 * computed in feature space.
 * 

 * Kernel PCA also has close connections with Isomap, LLE, and Laplacian eigenmaps.
 *
 * 
References
 * 
 * Bernhard Scholkopf, Alexander Smola, and Klaus-Robert Muller. Nonlinear Component Analysis as a Kernel Eigenvalue Problem. Neural Computation, 1998.
 * 
 *
 * @see smile.math.kernel.MercerKernel
 * @see PCA
 * @see smile.manifold.IsoMap
 * @see smile.manifold.LLE
 * @see smile.manifold.LaplacianEigenmap
 * @see smile.mds.SammonMapping
 *
 * @author Haifeng Li
 */
public class KPCA implements Projection {

    /**
     * The dimension of feature space.
     */
    private int p;
    /**
     * Learning data.
     */
    private T[] data;
    /**
     * Mercer kernel.
     */
    private MercerKernel kernel;
    /**
     * The row mean of kernel matrix.
     */
    private double[] mean;
    /**
     * The mean of kernel matrix.
     */
    private double mu;
    /**
     * Eigenvalues of kernel principal components.
     */
    private double[] latent;
    /**
     * Projection matrix.
     */
    private double[][] projection;
    /**
     * The coordinates of projected training data.
     */
    private double[][] coordinates;

    /**
     * Constructor. Learn kernel principal component analysis.
     * @param data learning data.
     * @param kernel Mercer kernel to compute kernel matrix.
     * @param threshold only principal components with eigenvalues larger than
     * the given threshold will be kept.
     */
    public KPCA(T[] data, MercerKernel kernel, double threshold) {
        this(data, kernel, data.length, threshold);
    }

    /**
     * Constructor. Learn kernel principal component analysis.
     * @param data learning data.
     * @param kernel Mercer kernel to compute kernel matrix.
     * @param k choose upto k principal components (larger than 0.0001) used for projection.
     */
    public KPCA(T[] data, MercerKernel kernel, int k) {
        this(data, kernel, k, 0.0001);
    }

    /**
     * Constructor. Constructor. Learn kernel principal component analysis.
     * @param data learning data.
     * @param kernel Mercer kernel to compute kernel matrix.
     * @param k choose top k principal components used for projection.
     * @param threshold only principal components with eigenvalues larger than
     * the given threshold will be kept.
     */
    public KPCA(T[] data, MercerKernel kernel, int k, double threshold) {
        if (threshold < 0) {
            throw new IllegalArgumentException("Invalid threshold = " + threshold);
        }

        if (k < 1 || k > data.length) {
            throw new IllegalArgumentException("Invalid dimension of feature space: " + k);
        }

        this.data = data;
        this.kernel = kernel;
        int n = data.length;

        double[][] K = new double[n][n];
        for (int i = 0; i < n; i++) {
            for (int j = 0; j <= i; j++) {
                K[i][j] = kernel.k(data[i], data[j]);
                K[j][i] = K[i][j];
            }
        }

        mean = Math.rowMean(K);
        mu = Math.mean(mean);

        for (int i = 0; i < n; i++) {
            for (int j = 0; j <= i; j++) {
                K[i][j] = K[i][j] - mean[i] - mean[j] + mu;
                K[j][i] = K[i][j];
            }
        }

        EigenValueDecomposition eigen = Math.eigen(K, k);

        p = 0;
        for (int i = 0; i < k; i++) {
            double e = (eigen.getEigenValues()[i] /= n);
            if (e > threshold) {
                p++;
            } else {
                break;
            }
        }

        latent = new double[p];
        projection = new double[p][n];
        for (int j = 0; j < p; j++) {
            latent[j] = eigen.getEigenValues()[j];
            double s = Math.sqrt(latent[j]);
            for (int i = 0; i < n; i++) {
                projection[j][i] = eigen.getEigenVectors()[i][j] / s;
            }
        }

        coordinates = new double[n][p];
        for (int i = 0; i < n; i++) {
            Math.ax(projection, K[i], coordinates[i]);
        }
    }

    /**
     * Returns the eigenvalues of kernel principal components, ordered from largest to smallest.
     */
    public double[] getVariances() {
        return latent;
    }

    /**
     * Returns the projection matrix. The dimension reduced data can be obtained
     * by y = W * K(x, ·).
     */
    public double[][] getProjection() {
        return projection;
    }

    /**
     * Returns the nonlinear principal component scores, i.e., the representation
     * of learning data in the nonlinear principal component space. Rows
     * correspond to observations, columns to components.
     */
    public double[][] getCoordinates() {
        return coordinates;
    }

    @Override
    public double[] project(T x) {
        int n = data.length;

        double[] y = new double[n];
        for (int i = 0; i < n; i++) {
            y[i] = kernel.k(x, data[i]);
        }

        double my = Math.mean(y);
        for (int i = 0; i < n; i++) {
            y[i] = y[i] - my - mean[i] + mu;
        }

        double[] z = new double[p];
        Math.ax(projection, y, z);
        return z;
    }

    @Override
    public double[][] project(T[] x) {
        int m = x.length;
        int n = data.length;
        double[][] y = new double[m][n];
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                y[i][j] = kernel.k(x[i], data[j]);
            }

            double my = Math.mean(y[i]);
            for (int j = 0; j < n; j++) {
                y[i][j] = y[i][j] - my - mean[j] + mu;
            }

        }

        double[][] z = new double[x.length][p];
        for (int i = 0; i < y.length; i++) {
            Math.ax(projection, y[i], z[i]);
        }
        return z;
    }
}