每日算法 — 使用java实现主成分分析：协方差矩阵特征分解与数据降维

引言

主成分分析（Principal Component Analysis, PCA）是最经典的无监督降维算法之一，由Karl Pearson于1901年提出。PCA通过线性变换将高维数据投影到低维空间，在保留数据最大方差的同时减少特征维度。本文将用Java完整实现PCA算法，深入讲解协方差矩阵计算、特征分解过程，并通过可视化展示降维效果。

算法原理

核心思想

PCA的目标是找到一组正交基（主成分），使得数据在该基上的投影方差最大。第一主成分方向是数据方差最大的方向，第二主成分方向是与第一主成分正交且方差次大的方向，依此类推。

数学推导

步骤1：数据标准化

对原始数据进行中心化处理：

$$\tilde{x}_i = x_i – \bar{x}$$

其中 $\bar{x}$ 为样本均值。

步骤2：计算协方差矩阵

$$C = \frac{1}{n-1} \sum_{i=1}^{n} \tilde{x}_i \tilde{x}_i^T = \frac{1}{n-1} X^T X$$

其中 $X$ 为中心化后的数据矩阵。

步骤3：特征分解

对协方差矩阵进行特征分解：

$$C v_i = \lambda_i v_i$$

其中 $\lambda_i$ 为特征值，$v_i$ 为对应的特征向量。

步骤4：选择主成分

按特征值从大到小排序，选择前 $k$ 个特征向量构成投影矩阵 $W = [v_1, v_2, …, v_k]$。

步骤5：数据投影

$$Y = X W$$

方差保留率

前 $k$ 个主成分保留的方差比例为：

$$\text{保留率} = \frac{\sum_{i=1}^{k} \lambda_i}{\sum_{i=1}^{d} \lambda_i}$$

Java实现

项目结构

src/
├── PCA.java                 // PCA核心类
├── MatrixUtils.java         // 矩阵运算工具类
├── DataPoint.java           // 数据点类
├── PCADemo.java             // 演示入口
└── DatasetGenerator.java    // 数据集生成器

矩阵运算工具类

/**
 * 矩阵运算工具类
 * 提供PCA所需的矩阵基本运算：乘法、转置、特征分解等
 */
public class MatrixUtils {

    /**
     * 矩阵乘法 C = A × B
     */
    public static double[][] multiply(double[][] A, double[][] B) {
        int rowsA = A.length;
        int colsA = A[0].length;
        int colsB = B[0].length;
        double[][] C = new double[rowsA][colsB];

        for (int i = 0; i < rowsA; i++) {
            for (int j = 0; j < colsB; j++) {
                double sum = 0;
                for (int k = 0; k < colsA; k++) {
                    sum += A[i][k] * B[k][j];
                }
                C[i][j] = sum;
            }
        }
        return C;
    }

    /**
     * 矩阵转置
     */
    public static double[][] transpose(double[][] A) {
        int rows = A.length;
        int cols = A[0].length;
        double[][] T = new double[cols][rows];
        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < cols; j++) {
                T[j][i] = A[i][j];
            }
        }
        return T;
    }

    /**
     * 向量点积
     */
    public static double dotProduct(double[] a, double[] b) {
        double sum = 0;
        for (int i = 0; i < a.length; i++) {
            sum += a[i] * b[i];
        }
        return sum;
    }

    /**
     * 向量归一化
     */
    public static double[] normalize(double[] v) {
        double norm = Math.sqrt(dotProduct(v, v));
        double[] result = new double[v.length];
        for (int i = 0; i < v.length; i++) {
            result[i] = v[i] / norm;
        }
        return result;
    }

    /**
     * 矩阵与向量乘法
     */
    public static double[] multiplyMatrixVector(double[][] A, double[] v) {
        int rows = A.length;
        double[] result = new double[rows];
        for (int i = 0; i < rows; i++) {
            result[i] = dotProduct(A[i], v);
        }
        return result;
    }

    /**
     * 打印矩阵
     */
    public static void printMatrix(double[][] A, String name) {
        System.out.println(name + ":");
        for (double[] row : A) {
            System.out.print("  ");
            for (double val : row) {
                System.out.printf("%10.4f ", val);
            }
            System.out.println();
        }
    }

    /**
     * 打印向量
     */
    public static void printVector(double[] v, String name) {
        System.out.print(name + ": [");
        for (int i = 0; i < v.length; i++) {
            System.out.printf("%.4f", v[i]);
            if (i < v.length - 1) System.out.print(", ");
        }
        System.out.println("]");
    }
}

PCA核心类

import java.util.Arrays;
import java.util.Comparator;

/**
 * 主成分分析（PCA）算法实现
 * 支持数据标准化、协方差矩阵计算、特征分解和降维投影
 */
public class PCA {
    private double[] mean;              // 各维度均值
    private double[] std;               // 各维度标准差
    private double[][] components;      // 主成分（特征向量）
    private double[] explainedVariance; // 各主成分解释的方差（特征值）
    private double[] explainedVarianceRatio; // 方差贡献率
    private int nComponents;            // 保留的主成分数
    private int originalDim;            // 原始维度
    private boolean standardize;        // 是否标准化

    public PCA(int nComponents, boolean standardize) {
        this.nComponents = nComponents;
        this.standardize = standardize;
    }

    /**
     * 拟合数据并执行PCA
     * @param data 输入数据矩阵，每行一个样本
     */
    public void fit(double[][] data) {
        int nSamples = data.length;
        originalDim = data[0].length;

        System.out.println("=== PCA 拟合开始 ===");
        System.out.println("样本数: " + nSamples);
        System.out.println("原始维度: " + originalDim);
        System.out.println("目标维度: " + nComponents);
        System.out.println("标准化: " + standardize);

        // 步骤1：计算均值
        mean = new double[originalDim];
        for (double[] sample : data) {
            for (int j = 0; j < originalDim; j++) {
                mean[j] += sample[j];
            }
        }
        for (int j = 0; j < originalDim; j++) {
            mean[j] /= nSamples;
        }
        MatrixUtils.printVector(mean, "均值向量");

        // 步骤2：中心化和标准化
        double[][] centered = new double[nSamples][originalDim];
        std = new double[originalDim];

        if (standardize) {
            // 计算标准差
            for (double[] sample : data) {
                for (int j = 0; j < originalDim; j++) {
                    std[j] += Math.pow(sample[j] - mean[j], 2);
                }
            }
            for (int j = 0; j < originalDim; j++) {
                std[j] = Math.sqrt(std[j] / (nSamples - 1));
            }
            MatrixUtils.printVector(std, "标准差向量");
        }

        for (int i = 0; i < nSamples; i++) {
            for (int j = 0; j < originalDim; j++) {
                centered[i][j] = data[i][j] - mean[j];
                if (standardize && std[j] > 1e-10) {
                    centered[i][j] /= std[j];
                }
            }
        }

        // 步骤3：计算协方差矩阵
        double[][] covMatrix = computeCovarianceMatrix(centered);
        MatrixUtils.printMatrix(covMatrix, "协方差矩阵");

        // 步骤4：特征分解（使用幂迭代法）
        EigenResult eigenResult = powerIteration(covMatrix, nComponents);
        components = eigenResult.eigenvectors;
        explainedVariance = eigenResult.eigenvalues;

        // 计算方差贡献率
        double totalVariance = 0;
        for (double v : explainedVariance) {
            totalVariance += v;
        }
        explainedVarianceRatio = new double[explainedVariance.length];
        for (int i = 0; i < explainedVariance.length; i++) {
            explainedVarianceRatio[i] = explainedVariance[i] / totalVariance;
        }

        // 输出结果
        System.out.println("\n特征值与方差贡献率：");
        for (int i = 0; i < explainedVariance.length; i++) {
            System.out.printf("  主成分 %d: λ=%.4f, 贡献率=%.2f%%%n",
                i + 1, explainedVariance[i], explainedVarianceRatio[i] * 100);
        }

        double cumulativeRatio = 0;
        for (int i = 0; i < nComponents; i++) {
            cumulativeRatio += explainedVarianceRatio[i];
        }
        System.out.printf("\n前 %d 个主成分累计方差保留率: %.2f%%%n",
            nComponents, cumulativeRatio * 100);
    }

    /**
     * 计算协方差矩阵
     */
    private double[][] computeCovarianceMatrix(double[][] centered) {
        int n = centered.length;
        int d = centered[0].length;
        double[][] cov = new double[d][d];

        for (int i = 0; i < d; i++) {
            for (int j = i; j < d; j++) {
                double sum = 0;
                for (int k = 0; k < n; k++) {
                    sum += centered[k][i] * centered[k][j];
                }
                cov[i][j] = sum / (n - 1);
                cov[j][i] = cov[i][j];
            }
        }
        return cov;
    }

    /**
     * 使用幂迭代法计算前k个最大特征值和特征向量
     */
    private EigenResult powerIteration(double[][] matrix, int k) {
        int n = matrix.length;
        double[][] eigenvectors = new double[k][n];
        double[] eigenvalues = new double[k];

        // 复制矩阵用于后续deflation
        double[][] workingMatrix = new double[n][n];
        for (int i = 0; i < n; i++) {
            System.arraycopy(matrix[i], 0, workingMatrix[i], 0, n);
        }

        for (int comp = 0; comp < k; comp++) {
            // 随机初始化特征向量
            double[] vector = new double[n];
            for (int i = 0; i < n; i++) {
                vector[i] = Math.random() - 0.5;
            }
            vector = MatrixUtils.normalize(vector);

            // 幂迭代
            double eigenvalue = 0;
            for (int iter = 0; iter < 1000; iter++) {
                double[] newVector = MatrixUtils.multiplyMatrixVector(workingMatrix, vector);
                double newEigenvalue = MatrixUtils.dotProduct(vector, newVector);
                newVector = MatrixUtils.normalize(newVector);

                // 检查收敛
                if (Math.abs(newEigenvalue - eigenvalue) < 1e-10) {
                    break;
                }

                vector = newVector;
                eigenvalue = newEigenvalue;
            }

            eigenvectors[comp] = vector;
            eigenvalues[comp] = eigenvalue;

            // Deflation: 移除当前特征向量的影响
            for (int i = 0; i < n; i++) {
                for (int j = 0; j < n; j++) {
                    workingMatrix[i][j] -= eigenvalue * vector[i] * vector[j];
                }
            }
        }

        return new EigenResult(eigenvectors, eigenvalues);
    }

    /**
     * 将数据投影到主成分空间
     */
    public double[][] transform(double[][] data) {
        int nSamples = data.length;
        double[][] projected = new double[nSamples][nComponents];

        for (int i = 0; i < nSamples; i++) {
            // 中心化和标准化
            double[] standardized = new double[originalDim];
            for (int j = 0; j < originalDim; j++) {
                standardized[j] = data[i][j] - mean[j];
                if (standardize && std[j] > 1e-10) {
                    standardized[j] /= std[j];
                }
            }

            // 投影到主成分
            for (int k = 0; k < nComponents; k++) {
                projected[i][k] = MatrixUtils.dotProduct(standardized, components[k]);
            }
        }

        return projected;
    }

    /**
     * 拟合并转换数据
     */
    public double[][] fitTransform(double[][] data) {
        fit(data);
        return transform(data);
    }

    /**
     * 从降维数据重构原始维度数据
     */
    public double[][] inverseTransform(double[][] projected) {
        int nSamples = projected.length;
        double[][] reconstructed = new double[nSamples][originalDim];

        for (int i = 0; i < nSamples; i++) {
            for (int j = 0; j < originalDim; j++) {
                double sum = mean[j];
                for (int k = 0; k < nComponents; k++) {
                    sum += projected[i][k] * components[k][j];
                    if (standardize && std[j] > 1e-10) {
                        // 反标准化需要在投影累加后进行
                        // 这里简化处理
                    }
                }
                reconstructed[i][j] = sum;
            }
        }
        return reconstructed;
    }

    public double[][] getComponents() {
        return components;
    }

    public double[] getExplainedVarianceRatio() {
        return explainedVarianceRatio;
    }

    /**
     * 内部类：存储特征分解结果
     */
    private static class EigenResult {
        double[][] eigenvectors;
        double[] eigenvalues;

        EigenResult(double[][] eigenvectors, double[] eigenvalues) {
            this.eigenvectors = eigenvectors;
            this.eigenvalues = eigenvalues;
        }
    }
}

数据集生成器

import java.util.Random;

/**
 * 生成用于PCA演示的数据集
 */
public class DatasetGenerator {
    private Random random;

    public DatasetGenerator(long seed) {
        this.random = new Random(seed);
    }

    /**
     * 生成具有强相关性的三维数据（实际分布接近二维平面）
     */
    public double[][] generateCorrelated3D(int n) {
        double[][] data = new double[n][3];
        for (int i = 0; i < n; i++) {
            double x = random.nextGaussian() * 3;
            double y = x * 0.8 + random.nextGaussian() * 0.5;
            double z = x * 0.6 + y * 0.3 + random.nextGaussian() * 0.3;
            data[i][0] = x;
            data[i][1] = y;
            data[i][2] = z;
        }
        return data;
    }

    /**
     * 生成高维数据（10维，但主要方差集中在前3维）
     */
    public double[][] generateHighDimData(int n, int dim) {
        double[][] data = new double[n][dim];
        for (int i = 0; i < n; i++) {
            // 前3维携带主要信息
            double v1 = random.nextGaussian() * 5;
            double v2 = random.nextGaussian() * 3;
            double v3 = random.nextGaussian() * 2;

            data[i][0] = v1;
            data[i][1] = v2;
            data[i][2] = v3;

            // 其余维度为噪声
            for (int j = 3; j < dim; j++) {
                data[i][j] = random.nextGaussian() * 0.5;
            }
        }
        return data;
    }
}

演示入口

/**
 * PCA算法演示
 */
public class PCADemo {
    public static void main(String[] args) {
        demo3DTo2D();
        System.out.println("\n" + "=".repeat(60) + "\n");
        demoHighDimReduction();
    }

    /**
     * 演示1：三维数据降维到二维
     */
    private static void demo3DTo2D() {
        System.out.println("=== 演示1：三维数据降维到二维 ===\n");

        DatasetGenerator generator = new DatasetGenerator(42);
        double[][] data3D = generator.generateCorrelated3D(100);

        System.out.println("原始三维数据样本（前5条）：");
        for (int i = 0; i < Math.min(5, data3D.length); i++) {
            System.out.printf("  样本%d: [%.2f, %.2f, %.2f]%n",
                i + 1, data3D[i][0], data3D[i][1], data3D[i][2]);
        }

        PCA pca = new PCA(2, true);
        double[][] data2D = pca.fitTransform(data3D);

        System.out.println("\n降维后二维数据样本（前5条）：");
        for (int i = 0; i < Math.min(5, data2D.length); i++) {
            System.out.printf("  样本%d: [%.4f, %.4f]%n",
                i + 1, data2D[i][0], data2D[i][1]);
        }

        // 打印主成分方向
        double[][] components = pca.getComponents();
        System.out.println("\n主成分方向：");
        for (int i = 0; i < components.length; i++) {
            System.out.printf("  PC%d: [", i + 1);
            for (int j = 0; j < components[i].length; j++) {
                System.out.printf("%.4f", components[i][j]);
                if (j < components[i].length - 1) System.out.print(", ");
            }
            System.out.println("]");
        }
    }

    /**
     * 演示2：十维数据降维
     */
    private static void demoHighDimReduction() {
        System.out.println("=== 演示2：十维数据降维分析 ===\n");

        DatasetGenerator generator = new DatasetGenerator(123);
        double[][] data10D = generator.generateHighDimData(200, 10);

        System.out.println("原始维度: 10");
        System.out.println("样本数: 200");

        // 尝试不同的降维维度
        for (int k : new int[]{2, 3, 5}) {
            PCA pca = new PCA(k, true);
            pca.fit(data10D);
            double[] ratio = pca.getExplainedVarianceRatio();

            double cumulative = 0;
            for (double r : ratio) {
                cumulative += r;
            }
            System.out.printf("\n降维到 %d 维: 累计方差保留率 = %.2f%%%n", k, cumulative * 100);
        }
    }
}

复杂度分析

算法应用

图像压缩

PCA可用于图像压缩：将图像分块后，对每个块进行PCA降维，只保留主成分系数，实现有损压缩。

人脸识别

Eigenfaces方法使用PCA提取人脸图像的主要特征，将高维像素空间降维到特征脸空间，实现高效的人脸识别。

数据可视化

将高维数据降维到2D或3D，便于可视化分析和模式发现。

总结

PCA通过协方差矩阵特征分解找到数据方差最大的方向，实现无监督降维。本文的Java实现完整展示了从数据标准化、协方差计算到幂迭代特征分解的全过程。PCA不仅是降维利器，更是理解数据内在结构的重要工具，为后续的流形学习、自编码器等深度降维方法奠定了基础。