核主成分分析KPCA+C代码

KPCA用非线性变换将输入数据空间映射到高维空间，使非线性问题转为线性问题，然后在高维空间中使用PCA方法提取主成分，在保持原数据信息量的基础上达到降维的目的。

常用的核函数有以下几种：

核函数化后的得到m*m的样本矩阵（m为样本个数）。用核函数将原始样本投射到高维空间，再用PCA进行降维。

实现步骤：

1. 将数据进行核函数化；

2. 对核矩阵样本进行归一化；归一化方法如下：

2. 之后用PCA进行降维，我上篇博客有：点击打开链接

实现代码：（代码是参考matlab代码和按自己的理解编写的，如有错误，欢迎指正）

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "vector"
using namespace std;


#define C 2   //常数
#define G 3   //gama值
#define S 1   //高斯核函数的方差
#define P 2   //多项式阶数
#define E 0.0000001
#define INF 99999 
#define dimNum 4     //维数
#define MAXITER 10000   //最大迭代次数


typedef vector<double> doubleVector;
typedef vector<doubleVector> dim2Vector;


vector<doubleVector> getInputSample(char* File);  //获取输入样本
void PAC(vector<doubleVector> nomalKernel, vector<doubleVector> kernelTrain);  //主成分分析法PAC
vector<doubleVector> normalizationSPSS(vector<doubleVector> inputTrain);  //采用z-score法标准数据
vector<doubleVector> calCovariation(vector<doubleVector> inputTrain);   //计算协方差
vector<dim2Vector> Jacobi(vector<doubleVector> Array);   //使用Jacobi计算协方差的特征值和特征矩阵
bool QueryArray(vector<doubleVector> Array);   //检查是否满足
vector<doubleVector> matTran(vector<doubleVector> Array);   //矩阵转置
vector<doubleVector> matMul(vector<doubleVector> mat1, vector<doubleVector> mat2);   //矩阵相乘


vector<doubleVector> kernelLinearFun(vector<doubleVector> inputTrain);  //使用线性核函数
vector<doubleVector> kernelPolyFun(vector<doubleVector> inputTrain);  //使用多项式核函数
vector<doubleVector> kernelGuassFun(vector<doubleVector> inputTrain);  //使用高斯径向基(RBF)核函数
vector<doubleVector> kernelTanhFun(vector<doubleVector> inputTrain);  //多层感知器(MLP)核函数
vector<doubleVector> normalizationKernel(vector<doubleVector> kernel);   //标准化核矩阵


double Input_Meam[dimNum] = {0};    //每一维的均值
double Input_Dev[dimNum] = {0};     //每一维的标准差


void main()
{
char *File = "input.txt";


vector<doubleVector> inputTrain;
vector<doubleVector> kernelTrain;
vector<doubleVector> nomalKernel;


inputTrain = getInputSample(File);
inputTrain = normalizationSPSS(inputTrain);


// kernelTrain = kernelLinearFun(inputTrain);    //使用线性核函数
// kernelTrain = kernelPolyFun(inputTrain);    //使用多项式核函数
// kernelTrain = kernelGuassFun(inputTrain);  //使用高斯径向基(RBF)核函数
kernelTrain = kernelTanhFun(inputTrain);  //使用多层感知器(MLP)核函数


nomalKernel = normalizationKernel(kernelTrain);   //标准化数据



PAC(nomalKernel,  kernelTrain);  //主成分分析法PAC


} 




//主成分分析法PAC
void PAC(vector<doubleVector> inputTrain,  vector<doubleVector> kernelTrain)
{
int i, j, m, n;
vector<doubleVector> input_Cov;  //协方差
vector<dim2Vector> jacobi;   //1为特征值，2为特征矩阵
double rate;  //贡献率
double rateSum1=0;
double rateSum2=0;
doubleVector tempVector;
vector<doubleVector> redTemp;  
vector<doubleVector> reduce_Dim_Mat;  //降维矩阵
vector<doubleVector> reduce_Dim_Sample;  //降维数据


input_Cov = calCovariation(inputTrain);   //计算协方差


jacobi = Jacobi(input_Cov);   //使用Jacobi计算协方差的特征值和特征矩阵


//计算贡献率
for(i=0; i<jacobi[0].size(); i++)
{
for(j=0; j<jacobi[0][i].size(); j++)
rateSum1 += jacobi[0][i][j];


for(j=0; j<jacobi[0][i].size(); j++)
{
rateSum2 += jacobi[0][i][j];
rate = rateSum2/rateSum1;
if(rate>=0.85)
break;
}


//获取将维矩阵并归一化（除以sqrt(特征值)）
for(m=0; m<=j; m++)
{
tempVector.clear();
for(n=0; n<jacobi[1][m].size(); n++)
tempVector.push_back(jacobi[1][n][m]/sqrtf(jacobi[0][0][m]));


reduce_Dim_Mat.push_back(tempVector);
}
}


reduce_Dim_Mat = matTran(reduce_Dim_Mat);


reduce_Dim_Sample = matMul(kernelTrain, reduce_Dim_Mat);  //计算降维结果


printf("协方差为:\n");
for(i=0; i<input_Cov.size(); i++)
{
for(j=0; j<input_Cov[i].size(); j++)
printf("%lf  ", input_Cov[i][j]);
printf("\n");
}


printf("\n特征值：\n");
for(i=0; i<jacobi[0].size(); i++)
{
for(j=0; j<jacobi[0][i].size(); j++)
printf("%lf  ", jacobi[0][i][j]);
printf("\n");
}

printf("\n特征向量：\n");
for(i=0; i<jacobi[1].size(); i++)
{
for(j=0; j<jacobi[1][i].size(); j++)
printf("%lf  ", jacobi[1][i][j]);
printf("\n");
  }




printf("\n降维矩阵：\n");
for(i=0; i<reduce_Dim_Mat.size(); i++)
{
for(j=0; j<reduce_Dim_Mat[i].size(); j++)
printf("%lf  ", reduce_Dim_Mat[i][j]);
printf("\n");
}


printf("\n降维结果：\n");
for(i=0; i<reduce_Dim_Sample.size(); i++)
{
for(j=0; j<reduce_Dim_Sample[i].size(); j++)
printf("%lf  ", reduce_Dim_Sample[i][j]);
printf("\n");
}

}




//计算协方差
vector<doubleVector> calCovariation(vector<doubleVector> inputTrain)
{
int i, j, k;
doubleVector tempDst(inputTrain[0].size(), 0);
vector<doubleVector> dst(inputTrain.size(), tempDst);

for(i=0 ; i<inputTrain.size(); i++)
Input_Meam[i] = 0;

//计算均值
for(i=0; i<inputTrain[0].size(); i++)
{
for(j=0; j<inputTrain.size(); j++)
Input_Meam[i] += inputTrain[j][i];

Input_Meam[i] = Input_Meam[i]/inputTrain.size();
}

//计算协方差
for(i=0; i<inputTrain[0].size(); i++)
for(j=0; j<inputTrain[0].size(); j++)
{
for(k=0; k<inputTrain.size(); k++)
dst[i][j] += (inputTrain[k][i]-Input_Meam[i])*(inputTrain[k][j]-Input_Meam[j]);

dst[i][j] = dst[i][j]/(inputTrain.size()-1);

}


return dst;
}






//使用Jacobi计算协方差的特征值和特征矩阵
vector<dim2Vector> Jacobi(vector<doubleVector> Array)
{
int i, j;
int count;
bool flag = false;
vector<dim2Vector> dst;
doubleVector tempArray(Array.size(), 0);
vector<doubleVector> charatMat(Array.size(), tempArray);   //特征向量
vector<doubleVector> sortArray;  //排序后的特征值
vector<doubleVector> dim2Jac;
vector<doubleVector> dim2JacT;
vector<dim2Vector> dim3Jac;
double maxArrayNum;
int laber_j, laber_i;

double theta;


//开始迭代
count = 0;
tempArray.clear();
tempArray.resize(Array.size(), 0);
while(count<MAXITER && !flag)
{
count++;
dim2Jac.clear();
dim2Jac.resize(Array.size(), tempArray);
maxArrayNum = 0;
laber_i = laber_j = 0;


//寻找非对角元中绝对值最大的A[i][j]
for(i=0; i<Array.size(); i++)
for(j=0; j<Array.size(); j++)
{
if(i==j)
continue;


if(maxArrayNum<fabs(Array[i][j]))
{
maxArrayNum = fabs(Array[i][j]);
laber_i = i;
laber_j = j;
}
}


theta = atanf(Array[laber_i][laber_j]*2/(Array[laber_i][laber_i]-Array[laber_j][laber_j]+E));


//构造雅克比矩阵
for(i=0; i<Array.size(); i++)
dim2Jac[i][i] = 1;


dim2Jac[laber_i][laber_i] = dim2Jac[laber_j][laber_j] = cosf(theta/2);
dim2Jac[laber_i][laber_j] = sinf(theta/2);
dim2Jac[laber_j][laber_i] = -sinf(theta/2);


dim2JacT = matTran(dim2Jac);  //矩阵转置
dim3Jac.push_back(dim2JacT);  //保存矩阵


Array = matMul(matMul(dim2Jac, Array), dim2JacT);


if(QueryArray(Array))
flag = true;

}


//初始化特征矩阵
for(i=0; i<Array.size(); i++)
charatMat[i][i] = 1;


//计算特征矩阵
for(i=0; i<dim3Jac.size(); i++)
charatMat = matMul(charatMat, dim3Jac[i]);


//排序
doubleVector sortA;
double tempNum;
for(i=0; i<Array.size(); i++)
sortA.push_back(Array[i][i]);


for(i=0; i<sortA.size(); i++)
{
maxArrayNum = sortA[i];
laber_j = i;


for(j=i; j<sortA.size(); j++)
if(maxArrayNum<sortA[j])
{
maxArrayNum = sortA[j];
laber_j = j;
}


tempNum = sortA[i];
sortA[i] = sortA[laber_j];
sortA[laber_j] = tempNum;


for(j=0; j<charatMat[laber_j].size(); j++)
tempArray[j] = charatMat[j][i];


for(j=0; j<charatMat[laber_j].size(); j++)
charatMat[j][i] = charatMat[j][laber_j];


for(j=0; j<charatMat[laber_j].size(); j++)
charatMat[j][laber_j] = tempArray[j]; 


}


sortArray.push_back(sortA);


dst.push_back(sortArray);
dst.push_back(charatMat);


return dst;
}




//检查是否满足
bool QueryArray(vector<doubleVector> Array)
{
int i, j;


for(i=0; i<Array.size(); i++)
for(j=0; j<Array.size(); j++)
{
if(i==j)
continue;


if(fabs(Array[i][j])>E)
return false;
}


return true;
}




//矩阵转置
vector<doubleVector> matTran(vector<doubleVector> Array)
{
int i, j;
doubleVector temp(Array.size(), 0);
vector<doubleVector> dst(Array[0].size(), temp);


for(i=0; i<Array.size(); i++)
for(j=0; j<Array[0].size(); j++)
dst[j][i] = Array[i][j];


return dst;
}


//矩阵相乘
vector<doubleVector> matMul(vector<doubleVector> mat1, vector<doubleVector> mat2)
{
  int i, j, k;
  doubleVector temp(mat2[0].size(), 0);
  vector<doubleVector> dst(mat1.size(), temp);


  for(i=0; i<mat1.size(); i++)
  for(j=0; j<mat2[0].size(); j++)
  for(k=0; k<mat2.size(); k++)
  dst[i][j] += mat1[i][k]*mat2[k][j];


return dst;
}




//采用z-score法标准数据
vector<doubleVector> normalizationSPSS(vector<doubleVector> inputTrain)
{
int i, j;
vector<doubleVector> dst;
doubleVector tempDst;

//初始化
for(i=0 ; i<dimNum; i++)
{
Input_Meam[i] = 0;
Input_Dev[i] = 0;
}

//计算均值
for(i=0; i<dimNum; i++)
{
for(j=0; j<inputTrain.size(); j++)
Input_Meam[i] += inputTrain[j][i];

Input_Meam[i] = Input_Meam[i]/inputTrain.size();
}

//计算标准差
for(i=0; i<dimNum; i++)
{
for(j=0; j<inputTrain.size(); j++)
Input_Dev[i] += (inputTrain[j][i]-Input_Meam[i])*(inputTrain[j][i]-Input_Meam[i]);

Input_Dev[i] = sqrtf(Input_Dev[i]/(inputTrain.size()-1));
}

//标准化
for(i=0; i<inputTrain.size(); i++)
{
tempDst.clear();
for(j=0; j<inputTrain[i].size(); j++)
tempDst.push_back((inputTrain[i][j]-Input_Meam[j])/Input_Dev[j]);

dst.push_back(tempDst);
}


return dst;
}




//核矩阵标准化
vector<doubleVector> normalizationKernel(vector<doubleVector> kernel)
{
int i, j;
doubleVector tempKernel(kernel[0].size(), 1);
vector<doubleVector> dstKernel(kernel.size(), tempKernel);
vector<doubleVector> onesMat(kernel.size(), tempKernel);
vector<doubleVector> K1, K2, K3;

K1 = matMul(onesMat, kernel);
K2 = matMul(kernel, onesMat);
K3 = matMul(matMul(onesMat, kernel), onesMat);

for(i=0; i<kernel.size(); i++)
for(j=0; j<kernel[i].size(); j++)
dstKernel[i][j] = kernel[i][j]-K1[i][j]-K2[i][j]+K3[i][j];

return dstKernel;

}




//获取输入样本
vector<doubleVector> getInputSample(char* File)
{
vector<doubleVector> dst;
doubleVector temp;
int i;
double num;


FILE *fp = fopen(File, "r");

if(fp == NULL)
{
printf("OPEN FILE ERROR!!\n");
exit(0);
}


//从文件读取样本
i=1;
temp.clear();
dst.clear();
while(fscanf(fp, "%lf", &num)!=EOF)
{
temp.push_back(num);
if(i%dimNum==0)
{
dst.push_back(temp);
temp.clear();
}
i++;
}


return dst;
}




//多层感知器(MLP)核函数
vector<doubleVector> kernelTanhFun(vector<doubleVector> inputTrain)
{
int i, j, k;
doubleVector tempKernel(inputTrain.size(), 0);
vector<doubleVector> dstKernel(inputTrain.size(), tempKernel);

for(i=0; i<inputTrain.size(); i++)
for(j=0; j<inputTrain.size(); j++)
{
for(k=0; k<inputTrain[i].size(); k++)
dstKernel[i][j] += inputTrain[i][k]*inputTrain[j][k];

dstKernel[i][j] = tanh(G*dstKernel[i][j]+C);
}

return dstKernel;
}




//使用高斯径向基(RBF)核函数
vector<doubleVector> kernelGuassFun(vector<doubleVector> inputTrain)
{
int i, j, k;
double sum;
doubleVector tempKernel(inputTrain.size(), 0);
vector<doubleVector> dstKernel(inputTrain.size(), tempKernel);

for(i=0; i<inputTrain.size(); i++)
for(j=0; j<inputTrain.size(); j++)
{
sum = 0;
//计算向量的2范数
for(k=0; k<inputTrain[i].size(); k++)
sum += (inputTrain[i][k]-inputTrain[j][k])*(inputTrain[i][k]-inputTrain[j][k]);

//高斯径向基(RBF)核函数计算公式
dstKernel[i][j] = exp(-0.5*(sqrtf(sum)/S)*(sqrtf(sum)/S));
}

return dstKernel;
}




//使用多项式核函数
vector<doubleVector> kernelPolyFun(vector<doubleVector> inputTrain)
{
int i, j, k, p;
double sum;
doubleVector tempKernel(inputTrain.size(), 0);
vector<doubleVector> dstKernel(inputTrain.size(), tempKernel);

for(i=0; i<inputTrain.size(); i++)
for(j=0; j<inputTrain.size(); j++)
{
for(k=0; k<inputTrain[i].size(); k++)
dstKernel[i][j] += inputTrain[i][k]*inputTrain[j][k];

sum = 1;
for(p=0; p<P; p++)
sum = sum*(dstKernel[i][j]+1);

dstKernel[i][j] = sum;
}

return dstKernel;
}


//使用线性核函数
vector<doubleVector> kernelLinearFun(vector<doubleVector> inputTrain)
{
int i, j, k;
doubleVector tempKernel(inputTrain.size(), 0);
vector<doubleVector> dstKernel(inputTrain.size(), tempKernel);

for(i=0; i<inputTrain.size(); i++)
for(j=0; j<inputTrain.size(); j++)
for(k=0; k<inputTrain[i].size(); k++)
dstKernel[i][j] += inputTrain[i][k]*inputTrain[j][k];

return dstKernel;
}