C++ 实现PCA 算法

ViperEk
371 阅读3分钟

C++ 实现PCA算法

1.公式

用二维数据举例

白数据：举个例子，我们有一组二维平面上的点 ,用(xi,yi)(x_i,y_i) 表示，xi,yix_i, y_i 的分布符合正态分布，则这一组数据 W=(x1x2...xny1y2...yn)W=\begin{pmatrix}x_1&x_2&...&x_n\\y_1&y_2&...&y_n \end{pmatrix} 则被称为白数据

对于其他各方向分量均值为0的来说数据DD来说，通过拉伸S=(a00b)S = \begin{pmatrix}a&0\\0&b\end{pmatrix}和旋转R=(cosθ−sinθsinθcosθ)R = \begin{pmatrix}cos\theta&-sin\theta\\sin\theta&cos\theta\end{pmatrix}这种数据可以被分解为D=RSWD=RSW

对于多个维度一般的可以得到公式D=RSWD=RSW

协方差 :cov(x,y)=1n−1Σi=1n(x−x^)(y−y^)cov(x,y)=\frac{1}{n-1}\Sigma_{i=1}^{n}(x-\hat{x})(y-\hat{y})

问题有来了为什么是n−1n-1 ？这和无偏估计有关，目前我个人暂且给不出证明

协方差矩阵：C=(cov(x,x)cov(x,y)cov(y,x)cov(y,y))C =\begin{pmatrix}cov(x,x)&cov(x,y)\\cov(y,x)&cov(y,y)\end{pmatrix} ，如果数据的平均值为0的话可以推理得到

C=1n−1(Σi=1nxi2Σi=1nxiyiΣi=1nxiyiΣi=1nyi2)=1n−1DDTC =\frac{1}{n-1}\begin{pmatrix}\Sigma_{i=1}^{n}x_i^{2}&\Sigma_{i=1}^{n}x_iy_i\\\Sigma_{i=1}^{n}x_iy_i&\Sigma_{i=1}^{n}y_i^{2}\end{pmatrix} = \frac{1}{n-1}DD^{T} 协方差矩阵的特征值就是对应轴向的方差，特征向量就是坐标轴方向，证明如下：

D=RSWD=RSW

cov(D)=1n−1DDT=1n−1RSWWTSTRTcov(D)=\frac{1}{n-1}DD^{T}=\frac{1}{n-1}RSWW^TS^TR^T

白数据协方差矩阵为单位矩阵1n−1WWT=E\frac{1}{n-1}WW^T = E

cov(D)=RS(1n−1WWT)STRTcov(D)=RS(\frac{1}{n-1}WW^T)S^TR^T

cov(D)=RSSTRTcov(D)=RSS^TR^T

S=(a00b)S = \begin{pmatrix}a&0\\0&b\end{pmatrix} ,令L=SST=(a200b2)L = SS^T = \begin{pmatrix}a^2&0\\0&b^2\end{pmatrix}

cov(D)=RLRTcov(D)=RLR^T

又因为旋转矩阵是正交矩阵，对于正交矩阵来说RRT=ERR^T = E 则可推 RT=R−1R^T = R^{-1}

cov(D)=RLR−1cov(D)=RLR^{-1} 令C=cov(D)C = cov(D)

实际上由特征向量公式AX=XΣAX=X\Sigma ，A=XΣX−1A = X\Sigma X^{-1} 可得到，CC的特征向量的集合就是旋转矩阵，CC 的特征值的集合就是方差

协方差矩阵是对称矩阵，自然而然的我们想到了Jacobi求特征值和特征向量

我们应该实现的求协方差矩阵特征向量的方法 jacobi.h

参考之前写的Jacobi方法：如何求对称矩阵的特征值和特征向量

2.代码实现

//mat.h 偷懒了懒得自己写了用嵌套vector
#pragma once
#include<vector>
#include<iostream>
typedef std::vector<std::vector<double>> mat_d;

void print(mat_d mat) {
    for (int i = 0; i < mat.size(); i++) {
        for (int j = 0; j < mat[0].size(); j++) {
            printf("%+5.3f\t", mat[i][j]);
        }
        printf("\n");
    }
    printf("\n");
}

mat_d get_mat_d(int i, int j) {
    return std::vector<std::vector<double>>
        (i, std::vector<double>(j, 0));
}


mat_d matmul(mat_d a, mat_d b) {
    if (a[0].size() != b.size()) {
        printf("Dim-matmul\n");
        exit(-1);
    };
    int m = a.size();
    int n = a[0].size();
    int p = b[0].size();
    mat_d res = get_mat_d(m, p);

    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            for (int k = 0; k < p; k++) {
                res[i][k] += a[i][j] * b[j][k];
            }
        }
    }
    return res;
}

mat_d transpose(mat_d a) {
    mat_d res = get_mat_d(a[0].size(), a.size());
    for (int i = 0; i < a.size(); i++) {
        for (int j = 0; j < a[0].size(); j++) {
            res[j][i] = a[i][j];
        }
    }
    return res;
}

mat_d mul(mat_d a, double t) {
    mat_d r = a;
    for (int i = 0; i < a.size(); i++) {
        for (int j = 0; j < a[0].size(); j++) {
            r[i][j] = r[i][j] * t;
        }
    }
    return r;
}


mat_d cov(mat_d a) {
    mat_d b = transpose(a);
    mat_d r = matmul(b,a);
    if (a.size() < 2) {
        printf("Dim\n");
        exit(-1);
    }
    r = mul(r, 1.0 / (a.size() - 1));
    return r;
}

PCA 算法主体

#pragma once
#include<algorithm>
#include"jacobi.h"
#include"mat.h"

/*
    PCA降维实现
    batch:      输入行特征行数目
    feature：   输入特征列数目
    Data:       输入特征的列方向平均后的数据
    TargetData: 返回的目标数据
    Cov:        协方差矩阵
    idx:        存方差最大的方向
*/

class PCA
{
private:
    using mat_d =
        std::vector<std::vector<double>>;
    struct Vec2d {
        int x;
        int y;
    };
    int batch;
    int feature;
    mat_d Data;
    mat_d TargetData;
    mat_d Cov;
    mat_d EigenVal;
    mat_d EigenVec;
    std::vector<Vec2d> idx;
private:
    mat_d mean(mat_d);
    void FindMaxFeatue();
public:
    PCA(mat_d&, int);
    mat_d getPCA();
    ~PCA() {}
};

PCA::PCA(mat_d& input, int target_dim) {
    //输入特征都是列向量 （批次，特征） 需要降低到的特征维度
    Data = input;
    Data = mean(Data);
    Cov = cov(Data);
    Jacobi jac(Cov, 1e-10, 100);
    EigenVal = jac.EigenVal();
    EigenVec = jac.EigenVec();

    Data = matmul(Data,EigenVec);
    FindMaxFeatue();

    TargetData = get_mat_d(Data.size(),0);
    for (int j = 0; j < target_dim; j++) {
        int cur_col = idx[j].x;
        for (int i = 0; i < Data.size(); i++) {
            TargetData[i].push_back(Data[i][cur_col]);
        }
    }

}

mat_d PCA::mean(mat_d data) {
    //求列方向上的平均值
    mat_d res = data;
    batch = data.size();
    feature = data[0].size();
    std::vector<double>  mean(feature);
    for (int i = 0; i < batch; i++) {
        for (int j = 0; j < feature; j++) {
            mean[j] += data[i][j];
        }
    }
    for (int i = 0; i < feature; i++) {
        mean[i] = mean[i] / batch;
    }

    for (int i = 0; i < batch; i++) {
        for (int j = 0; j < feature; j++) {
            res[i][j] = data[i][j] - mean[j];
        }
    }
    return res;
}
void PCA::FindMaxFeatue() {
    for (int i = 0; i < EigenVal.size(); i++) {
        Vec2d temp;
        temp.x = i;
        temp.y = EigenVal[i][i];
        idx.push_back(temp);
    }
    std::sort(idx.begin(), idx.end(),
        [](Vec2d a, Vec2d b) {
            return a.y > b.y;
        });
}

mat_d PCA::getPCA() {
    return TargetData;
}
avatar
文章
阅读
粉丝