PCA方法从原理到实现
深度神经网路已经在语音识别,图像识别等领域取得前所未有的成功。本人在多年之前也曾接触过神经网络。本系列文章主要记录自己对深度神经网络的一些学习心得。
第五篇,谈谈PCA模型。本来PCA模型与深度学习是没有任何联系的。通常我们只是用PCA来对机器学习的数据做预处理。
本来想详细记录一下PCA的原理,但网上已经有一篇不错的文章,链接如下:
http://hi.baidu.com/ifengzh/item/8851b6387aebefc4382ffa60
本文前面部分内容引用了这篇文章的内容。
一、简介
PCA(Principal Components Analysis)即主成分分析,是图像处理中经常用到的降维方法,大家知道,我们在处理有关数字图像处理方面的问题时,比如经常用的图像的查询问题,在一个几万或者几百万甚至更大的数据库中查询一幅相近的图像。这时,我们通常的方法是对图像库中的图片提取响应的特征,如颜色,纹理,sift,surf,vlad等等特征,然后将其保存,建立响应的数据索引,然后对要查询的图像提取相应的特征,与数据库中的图像特征对比,找出与之最近的图片。这里,如果我们为了提高查询的准确率,通常会提取一些较为复杂的特征,如sift,surf等,一幅图像有很多个这种特征点,每个特征点又有一个相应的描述该特征点的128维的向量,设想如果一幅图像有300个这种特征点,那么该幅图像就有300*vector(128维)个,如果我们数据库中有一百万张图片,这个存储量是相当大的,建立索引也很耗时,如果我们对每个向量进行PCA处理,将其降维为64维,是不是很节约存储空间啊?对于学习图像处理的人来说,都知道PCA是降维的,但是,很多人不知道具体的原理,为此,我写这篇文章,来详细阐述一下PCA及其具体计算过程:
二、PCA原理
1、原始数据:
为了方便,我们假定数据是二维的,借助网络上的一组数据,如下:
x=[2.5, 0.5, 2.2, 1.9, 3.1, 2.3, 2, 1,1.5, 1.1]T
y=[2.4, 0.7, 2.9, 2.2, 3.0, 2.7, 1.6, 1.1, 1.6, 0.9]T
2、计算协方差矩阵
什么是协方差矩阵?相信看这篇文章的人都学过数理统计,一些基本的常识都知道,但是,也许你很长时间不看了,都忘差不多了,为了方便大家更好的理解,这里先简单的回顾一下数理统计的相关知识,当然如果你知道协方差矩阵的求法你可以跳过这里。
(1)协方差矩阵:
首先我们给你一个含有n个样本的集合,依次给出数理统计中的一些相关概念:
均值:
标准差:
方差:
既然我们都有这么多描述数据之间关系的统计量,为什么我们还要用协方差呢?我们应该注意到,标准差和方差一般是用来描述一维数据的,但现实生活我们常常遇到含有多维数据的数据集,最简单的大家上学时免不了要统计多个学科的考试成绩。面对这样的数据集,我们当然可以按照每一维独立的计算其方差,但是通常我们还想了解这几科成绩之间的关系,这时,我们就要用协方差,协方差就是一种用来度量两个随机变量关系的统计量,其定义为:
从协方差的定义上我们也可以看出一些显而易见的性质,如:
需要注意的是,协方差也只能处理二维问题,那维数多了自然就需要计算多个协方差,比如n维的数据集就需要计算CN2【此乃组合数基本公式】个协方差,那自然而然的我们会想到使用矩阵来组织这些数据。给出协方差矩阵的定义:
这个定义还是很容易理解的,我们可以举一个简单的三维的例子,假设数据集有三个维度{x,y,z},则协方差矩阵为
可见,协方差矩阵是一个对称的矩阵,而且对角线是各个维度上的方差。
(2)协方差矩阵的求法:
协方差矩阵计算的是不同维度之间的协方差,而不是不同样本之间的。下面我们将在matlab中用一个例子进行详细说明:
首先,随机产生一个10*3维的整数矩阵作为样本集,10为样本的个数,3为样本的维数。
- MySample = fix(rand(10,3)*50)
根据公式,计算协方差需要计算均值,那是按行计算均值还是按列呢,我一开始就老是困扰这个问题。前面我们也特别强调了,协方差矩阵是计算不同维度间的协方差,要时刻牢记这一点。样本矩阵的每行是一个样本,每列为一个维度,所以我们要按列计算均值。为了描述方便,我们先将三个维度的数据分别赋值:
- dim1 = MySample(:,1);
- dim2 = MySample(:,2);
- dim3 = MySample(:,3);
- %计算dim1与dim2,dim1与dim3,dim2与dim3的协方差:
- sum( (dim1-mean(dim1)) .*(dim2-mean(dim2)) ) / ( size(MySample,1)-1 ) %得到 74.5333
- sum( (dim1-mean(dim1)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到 -10.0889
- sum( (dim2-mean(dim2)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到 -10***000
- %搞清楚了这个后面就容易多了,协方差矩阵的对角线就是各个维度上的方差,下面我们依次计算:
- std(dim1)^2 % 得到 108.3222
- std(dim2)^2 % 得到 260.6222
- std(dim3)^2 % 得到 94.1778
- %这样,我们就得到了计算协方差矩阵所需要的所有数据,调用Matlab自带的cov函数进行验证:
- cov(MySample)
可以看到跟我们计算的结果是一样的,说明我们的计算是正确的。但是通常我们不用这种方法,而是用下面简化的方法进行计算:
先让样本矩阵中心化,即每一维度减去该维度的均值,然后直接用新的到的样本矩阵乘上它的转置,然后除以(N-1)即可。其实这种方法也是由前面的公式通道而来,只不过理解起来不是很直观而已。大家可以自己写个小的矩阵看一下就明白了。其Matlab代码实现如下:
- X = MySample –repmat(mean(MySample),10,1); %中心化样本矩阵
- C = (X’*X)./(size(X,1)-1)
- %为方便对matlab不太明白的人,小小说明一下各个函数,同样,对matlab有一定基础的人直接跳过:
- %B = repmat(A,m,n ) %%将矩阵 A复制 m×n块,即把 A 作为 B的元素,B由 m×n个 A平铺而成。B的维数是 [size(A,1)*m, (size(A,2)*n]
- %B = mean(A)的说明:
- %如果你有这样一个矩阵:A = [1 2 3; 3 36; 4 6 8; 4 7 7];
- %用mean(A)(默认dim=1)就会求每一列的均值
- % ans =
- % 3.0000 4.5000 6.0000
- % 用mean(A,2)就会求每一行的均值
- % ans =
- % 2.0000
- % 4.0000
- % 6.0000
- % 6.0000
- size(A,n)%% 如果在size函数的输入参数中再添加一项n,并用1或2为n赋值,则 size将返回矩阵的行数或列数。其中r=size(A,1)该语句返回的是矩阵A的行数, %c=size(A,2)该语句返回的是矩阵A的列数
上面我们简单说了一下协方差矩阵及其求法,言归正传,我们用上面简化求法,求出样本的协方差矩阵为:
3、计算协方差矩阵的特征向量和特征值
因为协方差矩阵为方阵,我们可以计算它的特征向量和特征值,如下:
- [eigenvectors,eigenvalues] = eig(cov)
我们可以看到这些矢量都是单位矢量,也就是它们的长度为1,这对PCA来说是很重要的。
4、选择成分组成模式矢量
求出协方差矩阵的特征值及特征向量之后,按照特征值由大到小进行排列,这将给出成分的重要性级别。现在,如果你喜欢,可以忽略那些重要性很小的成分,当然这会丢失一些信息,但是如果对应的特征值很小,你不会丢失很多信息。如果你已经忽略了一些成分,那么最后的数据集将有更少的维数,精确地说,如果你的原始数据是n维的,你选择了前p个主要成分,那么你现在的数据将仅有p维。现在我们要做的是组成一个模式矢量,这只是几个矢量组成的矩阵的一个有意思的名字而已,它由你保持的所有特征矢量构成,每一个特征矢量是这个矩阵的一列。
对于我们的数据集,因为有两个特征矢量,因此我们有两个选择。我们可以用两个特征矢量组成模式矢量:
我们也可以忽略其中较小特征值的一个特征矢量,从而得到如下模式矢量:
5、得到降维后的数据
其中rowFeatureVector是由模式矢量作为列组成的矩阵的转置,因此它的行就是原来的模式矢量,而且对应最大特征值的特征矢量在该矩阵的最上一行。rowdataAdjust是每一维数据减去均值后,所组成矩阵的转置,即数据项目在每一列中,每一行是一维,对我们的样本来说即是,第一行为x维上数据,第二行为y维上的数据。FinalData是最后得到的数据,数据项目在它的列中,维数沿着行。
这将给我们什么结果呢?这将仅仅给出我们选择的数据。我们的原始数据有两个轴(x和y),所以我们的原始数据按这两个轴分布。我们可以按任何两个我们喜欢的轴表示我们的数据。如果这些轴是正交的,这种表达将是最有效的,这就是特征矢量总是正交的重要性。我们已经将我们的数据从原来的xy轴表达变换为现在的单个特征矢量表达。
说明:如果要恢复原始数据,只需逆过程计算即可,即:
到此为止,相信你已经掌握了PCA的原理了。
三 . PCA的应用
PCA及其改进算法主要应用的人脸识别领域,是人脸识别的经典算法之一。OpenCv2.4以后的版本实现了三种经典的人脸识别算法,其中就包括PCA。对openCv比较老的版本也可以调用PCA的算法去做,只是稍显复杂而已,网上有一篇博文如下:
http://www.cognotics.com/opencv/servo_2007_series/part_5/index.html
该代码运行在openCv2.1之前的版本当中,但是该代码有个重要的bug就是特征数K被设置为固定的值,而选择更小的值的时候,代码会crash。
PCA另外一个主要的用途是作为其他算法的预处理,术语叫做数据的白化。由于PCA具有压缩数据的作用,所以可以认为经过PCA处理过之后的数据是不相关的,但一般未必是独立的。实际可用的PCA算法一般不是以解析解的形式给出的,而是在线学习算法。有很多的原因决定了只能使用在线学习算法。在线学习算法主要有基于神经网络学习的算法和递归最小二乘法,相关的文献如下:
http://wenku.baidu.com/view/c91f31c058f5f61fb73666f8.html
要注意的是openCv的实现不是在线学习算法。
四. PCA 的实现
- #ifndef _FCE__PCA__H__
- #define _FCE__PCA__H__
- #define HIGH_PRECISON
- #ifdef HIGH_PRECISON
- #define real float
- #else
- #define real double
- #endif
- #ifdef _cplusplus
- {
- extern "C"
- #endif
- typedef struct _FCE_PCA{
- int count; //the number of sample
- int n; // the number of features
- real *covariance;
- real *mean;
- real *z;
- }FCE_PCA;
- FCE_PCA *fce_pca_init(int n);
- void fce_pca_push_add(FCE_PCA *pca, real *v);
- int fce_pca_solve_eig(FCE_PCA *pca, real *eigenvector, real *eigenvalue);
- void fce_pca_free(FCE_PCA *pca);
- #ifdef _cplusplus
- }
- #endif
- #endif
- #include "fce_pca.h"
- #define FCE_MIN(i,j) (((i) > (j)) ? (j) : (i))
- #define FCE_MAX(i,j) (((i) > (j)) ? (i) : (j))
- FCE_PCA *fce_pca_init(int n){
- FCE_PCA *pca;
- real zero = 0.0;
- if(n <= 1)
- return NULL;
- pca = (FCE_PCA* )malloc(sizeof(FCE_PCA));
- if (pca == NULL){
- return NULL;
- }
- pca->n = n;
- pca->z = (real* )malloc(sizeof(*pca->z) * n);
- if (pca->z == NULL){
- free(pca);
- return NULL;
- }
- memset(pca->z, zero, sizeof(*pca->z) * n);
- pca->count=0;
- pca->covariance= (real* )malloc(sizeof(real) * n * n);
- if (pca->covariance == NULL){
- free(pca->z);
- free(pca);
- return NULL;
- }
- memset(pca->covariance, zero, sizeof(real) * n * n);
- pca->mean = (real* )malloc(sizeof(real) * n);
- if (pca->mean == NULL){
- free(pca->covariance);
- free(pca->z);
- free(pca);
- return NULL;
- }
- memset(pca->mean, zero, sizeof(real) * n);
- return pca;
- }
- void fce_pca_free(FCE_PCA *pca){
- free(pca->covariance);
- free(pca->mean);
- free(pca->z);
- free(pca);
- }
- void fce_pca_push_add(FCE_PCA *pca, real *v){
- int i, j;
- const int n = pca->n;
- for(i = 0; i < n; i++){
- pca->mean[i] += v[i];
- for(j = i; j < n; j++)
- pca->covariance[j + i * n] += v[i]*v[j];
- }
- pca->count++;
- }
- int fce_pca_solve_eig(FCE_PCA *pca, real *eigenvector, real *eigenvalue){
- int i, j, pass;
- int k = 0;
- const int n = pca->n;
- real *z = pca->z;
- real zero = 0.0;
- memset(eigenvector, zero, sizeof(real)*n*n);
- for(j = 0; j < n; j++){
- pca->mean[j] /= pca->count;
- eigenvector[j + j * n] = 1.0;
- for(i = 0; i <= j; i++){
- pca->covariance[j + i * n] /= pca->count;
- pca->covariance[j + i * n] -= pca->mean[i] * pca->mean[j];
- pca->covariance[i + j * n] = pca->covariance[j + i * n];
- }
- eigenvalue[j] = pca->covariance[j + j*n];
- z[j] = 0;
- }
- for(pass=0; pass < 50; pass++){
- real sum = 0;
- for(i = 0; i < n; i++)
- for(j = i+1; j < n; j++)
- sum += fabs(pca->covariance[j + i * n]);
- if(sum == 0){
- for(i = 0; i < n; i++){
- real maxvalue = -1;
- for(j = i; j < n; j++){
- if(eigenvalue[j] > maxvalue){
- maxvalue = eigenvalue[j];
- k= j;
- }
- }
- eigenvalue[k] = eigenvalue[i];
- eigenvalue[i] = maxvalue;
- for(j = 0; j < n; j++){
- real tmp = eigenvector[k + j * n];
- eigenvector[k + j * n] = eigenvector[i + j * n];
- eigenvector[i + j * n] = tmp;
- }
- }
- return pass;
- }
- for(i = 0; i < n; i++){
- for(j = i + 1; j < n; j++){
- real covar = pca->covariance[j + i * n];
- real t,c,s,tau,theta, h;
- if(pass < 3 && fabs(covar) < sum / (5*n*n))
- continue;
- if(fabs(covar) <= 0.00000000001)
- continue;
- if(pass >=3 && fabs((eigenvalue[j]+z[j])/covar) > (1LL<<32) && fabs((eigenvalue[i]+z[i])/covar) > (1LL<<32)){
- pca->covariance[j + i * n]=0.0;
- continue;
- }
- h = (eigenvalue[j] + z[j]) - (eigenvalue[i] + z[i]);
- theta = 0.5 * h/covar;
- t = 1.0 /(fabs(theta) + sqrt(1.0 + theta * theta));
- if(theta < 0.0) t = -t;
- c = 1.0 /sqrt(1 + t * t);
- s = t * c;
- tau = s /(1.0 + c);
- z[i] -= t * covar;
- z[j] += t * covar;
- #define ROTATE(a,i,j,k,l) {\
- real g =a[j + i*n];\
- real h =a[l + k*n];\
- a[j + i*n] = g - s * (h + g * tau);\
- a[l + k*n] = h + s * (g - h * tau); }
- for(k = 0; k < n; k++) {
- if(k != i && k != j){
- ROTATE(pca->covariance,FCE_MIN(k,i),FCE_MAX(k,i),FCE_MIN(k,j),FCE_MAX(k,j))
- }
- ROTATE(eigenvector,k,i,k,j)
- }
- pca->covariance[j + i * n]=0.0;
- }
- }
- for (i = 0; i < n; i++) {
- eigenvalue[i] += z[i];
- z[i]=0.0;
- }
- }
- return 0;
- }
PCA方法从原理到实现相关推荐
- PCA的数学原理(通俗易懂)
PCA(Principal Component Analysis)是一种常用的数据分析方法.PCA通过线性变换将原始数据变换为一组各维度线性无关的表示,可用于提取数据的主要特征分量,常用于高维数据的降 ...
- PCA的数学原理(非常值得阅读)!!!!
PCA(Principal Component Analysis)是一种常用的数据分析方法.PCA通过线性变换将原始数据变换为一组各维度线性无关的表示,可用于提取数据的主要特征分量,常用于高维数据的降 ...
- 主成分分析(PCA)数学原理详解
文章目录 1 数据的向量表示及降维问题 2 向量的表示及基变换 2.1 内积与投影 2.2 基 2.3 基变换的矩阵表示 3 协方差矩阵及优化目标 3.1 方差 3.2 协方差 3.3 协方差矩阵 3 ...
- 一文读懂PCA分析 (原理、算法、解释和可视化)
生物信息学习的正确姿势 NGS系列文章包括NGS基础.高颜值在线绘图和分析.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流 ...
- PCA方法的分解过程_协方差矩阵+特征值分解+降维投影
文章目录 问题引入 PCA的实现 1. 数据集加载 2. 数据预处理 3. 主成分分解 3.1 协方差矩阵构造 3.2 特征值分解 3.3 方差解释 3.4 降维投影 3.5 降维结果 问题总结 附录 ...
- 模式识别PCA方法实现人脸识别-Python
模式识别PCA方法实现人脸识别-Python 一.学习目标 二.学习内容 三.学习时间 四.学习产出 五.正文部分 1.什么是PCA? 2.使用PCA方法进行人脸识别的步骤 1)前期尝试及准备 (1) ...
- 异常检测-PCA方法
文章目录 异常检测-PCA方法 线性回归-最小二乘拟合 目标函数的数学意义 正规方程解决最小二乘问题 PCA用于异常检测 PyOD实例 异常检测-PCA方法 记录DataWhale的异常检测的学习过程 ...
- 多线程下ArrayList类线程不安全的解决方法及原理
多线程下ArrayList类线程不安全的解决方法及原理 参考文章: (1)多线程下ArrayList类线程不安全的解决方法及原理 (2)https://www.cnblogs.com/fangting ...
- 可逆加密算法 php,php可逆加密的方法及原理
本篇文章主要介绍php可逆加密的方法及原理,感兴趣的朋友参考下,希望对大家有所帮助. PHP代码如下:<?php class encryptCalss { var $key=12; functi ...
- 【Python-ML】无监督线性降维PCA方法
# -*- coding: utf-8 -*- ''' Created on 2018年1月18日 @author: Jason.F @summary: 特征抽取-PCA方法,无监督.线性 ''' i ...
最新文章
- 开发日记-20190429 关键词 患病 NDK
- 什么是计算机网络?—Vecloud微云
- VTK修炼之道38:图像平滑_中值滤波器
- stm32中断优先级快速入门
- SQL Server 查询数据库中所有的表名及行数
- python generator输出_python 高级特性:Generator(生成器)
- html5 - history 历史管理
- What is AJAX?(转)(二)
- 这是我的卡,去买个包吧
- PacMan开发-碰撞检测实现
- ObjectDataSource 如何传递查询参数
- 蓝色简洁的企业cms网站权限后台管理模板——后台
- maya中英文对比_maya2017中英文对照表.doc
- 2毫秒c51汇编语言延时函数,单片机精确毫秒延时函数
- 【爬虫实战】国家企业公示网-crawler爬虫抓取数据
- 练习2-6 编写一个函数setbits(x, p ,n, y),该函数返回对x执行下列操作后的结果值: 将x中从第p位开始的n个(二进制)位设置为y中最右边n位的值,x的其余各位保持不变。
- 使用Java API访问HFDS
- 【Lintcode】1647. Path Search
- 贪心算法-大整数乘法/加法/减法
- 小米6 android p是什么,小米6什么配置参数?小米6标配有什么?
热门文章
- Objective-C 继承新的认识以及作用
- centos5.5安装csvn,以及问题处理
- 开发小技巧: 如何在jQuery中禁用或者启用滚动事件.scroll
- 局域网邮件服务器搭建地址薄更新,搭建局域网邮件服务器
- Linux执行U盘里内程序,Linux 最小系统挂载U盘(SD、TF卡)并执行程序
- hdu 2094 “产生冠军”——set容器的应用
- php网页制作头部和尾部,用phpcms如何将静态页面制作成企业网站,头部加尾部
- 图解think php,图解ThinkPHP5框架(三):配置类Config.php源码解读
- pytorch中模型结构图的可视化
- C程序设计--排序(冒泡、选择、插入)--插入