核主成份分析的MatLab实现
来源:互联网 发布:数据库decode函数功能 编辑:程序博客网 时间:2024/04/27 20:30
最近准备往核方法方向发展,于是仔细研读了经典的核主成份分析方法。通过参考别人的源代码,自己实现了KPCA,基本了解了核方法的实质。个人体会亲手实现一下还是很有益的。matlab源代码如下:
function [eigvector, eigvalue,Y] = KPCA(X,r,opts)
%
% Kernel Principal Component Analysis
% [eigvector, eigvalue,Y] = KPCA(X,r,opts)
% Input:
% X: d*N data matrix;Each column vector of X is a sample vector.
% r: Dimensionality of reduced space (default: d)
% opts: Struct value in Matlab. The fields in options that can be set:
% KernelType - Choices are:
% 'Gaussian' - exp{-gamma(|x-y|^2)}
% 'Polynomial' - (x'*y)^d
% 'PolyPlus' - (x'*y+1)^d
% gamma - parameter for Gaussian kernel
% d - parameter for polynomial kernel
%
% Output:
% eigvector: N*r matrix;Each column is an embedding function, for a new
% data point (column vector) x, y = eigvector'*K(x,:)'
% will be the embedding result of x.
% K(x,:) = [K(x1,x),K(x2,x),...K(xN,x)]
% eigvalue: The sorted eigvalue of KPCA eigen-problem.
% Y : Data matrix after the nonlinear transform
if nargin<1
error('Not enough input arguments.')
end
[d,N]=size(X);
if nargin<2
r=d;
end
%% Ensure r is not bigger than d
if r>d
r=d;
end;
% Construct the Kernel matrix K
K =ConstructKernelMatrix(X,[],opts);
% Centering kernel matrix
One_N=ones(N)./N;
Kc = K - One_N*K - K*One_N + One_N*K*One_N;
clear One_N;
% Solve the eigenvalue problem N*lamda*alpha = K*alpha
if N>1000 && r
% using eigs to speed up!
opts.disp=0;
[eigvector, eigvalue] = eigs(Kc,r,'la',opts);
eigvalue = diag(eigvalue);
else
[eigvector, eigvalue] = eig(Kc);
eigvalue = diag(eigvalue);
[junk, index] = sort(-eigvalue);
eigvalue = eigvalue(index);
eigvector = eigvector(:,index);
end
if r < length(eigvalue)
eigvalue = eigvalue(1:r);
eigvector = eigvector(:, 1:r);
end
% Only reserve the eigenvector with nonzero eigenvalues
maxEigValue = max(abs(eigvalue));
eigIdx = find(abs(eigvalue)/maxEigValue < 1e-6);
eigvalue (eigIdx) = [];
eigvector (:,eigIdx) = [];
% Normalizing eigenvector
for i=1:length(eigvalue)
eigvector(:,i)=eigvector(:,i)/sqrt(eigvalue(i));
end;
if nargout >= 3
% Projecting the data in lower dimensions
Y = eigvector'*K;
end
function K=ConstructKernelMatrix(X,Y,opts)
%
% function K=ConstructKernelMatrix(X,Y,opts)
% Usage:
% opts.KernelType='Gaussian';
% K = ConstructKernelMatrix(X,[],opts)
% K = ConstructKernelMatrix(X,Y,opts)
%
% Input:
% X: d*N data matrix;Each column vector of X is a sample vector.
% Y: d*M data matrix;Each column vector of Y is a sample vector.
% opts: Struct value in Matlab. The fields in options that can be set:
% KernelType - Choices are:
% 'Gaussian' - exp{-gamma(|x-y|^2)}
% 'Polynomial' - (x'*y)^d
% 'PolyPlus' - (x'*y+1)^d
% gamma - parameter for Gaussian kernel
% d - parameter for polynomial kernel
% Output:
% K N*N or N*M matrix
if nargin<1
error('Not enough input arguments.')
end
if (~exist('opts','var'))
opts = [];
else
if ~isstruct(opts)
error('parameter error!');
end
end
N=size(X,2);
if isempty(Y)
K=zeros(N,N);
else
M=size(Y,2);
if size(X,1)~=size(Y,1)
error('Matrixes X and Y should have the same row dimensionality!');
end;
K=zeros(N,M);
end;
%=================================================
if ~isfield(opts,'KernelType')
opts.KernelType = 'Gaussian';
end
switch lower(opts.KernelType)
case {lower('Gaussian')} % exp{-gamma(|x-y|^2)}
if ~isfield(opts,'gamma')
opts.gamma = 0.5;
end
case {lower('Polynomial')} % (x'*y)^d
if ~isfield(opts,'d')
opts.d = 1;
end
case {lower('PolyPlus')} % (x'*y+1)^d
if ~isfield(opts,'d')
opts.d = 1;
end
otherwise
error('KernelType does not exist!');
end
switch lower(opts.KernelType)
case {lower('Gaussian')}
if isempty(Y)
for i=1:N
for j=i:N
dist = sum(((X(:,i) - X(:,j)).^2));
temp=exp(-opts.gamma*dist);
K(i,j)=temp;
if i~=j
K(j,i)=temp;
end;
end
end
else
for i=1:N
for j=1:M
dist = sum(((X(:,i) - Y(:,j)).^2));
K(i,j)=exp(-opts.gamma*dist);
end
end
end
case {lower('Polynomial')}
if isempty(Y)
for i=1:N
for j=i:N
temp=(X(:,i)'*X(:,j))^opts.d;
K(i,j)=temp;
if i~=j
K(j,i)=temp;
end;
end
end
else
for i=1:N
for j=1:M
K(i,j)=(X(:,i)'*Y(:,j))^opts.d;
end
end
end
case {lower('PolyPlus')}
if isempty(Y)
for i=1:N
for j=i:N
temp=(X(:,i)'*X(:,j)+1)^opts.d;
K(i,j)=temp;
if i~=j
K(j,i)=temp;
end;
end
end
else
for i=1:N
for j=1:M
K(i,j)=(X(:,i)'*Y(:,j)+1)^opts.d;
end
end
end
otherwise
error('KernelType does not exist!');
end
参考文献:
[1]B. Sch¨olkopf, A.J. Smola, and K.-R M¨uller, “Nonlinear component analysis as a kernel eigenvalue problem,”Neural Computation, vol. 10, pp. 1299–1319, 1998.
[2]KERNEL PRINCIPAL COMPONENT ANALYSIS FOR FEATURE REDUCTION IN HYPERSPECTRALE IMAGES ANALYSIS
[3]http://www.cs.uiuc.edu/homes/dengcai2/Data/data.html
KPCA 算法:
引用地址:http://blog.sciencenet.cn/blog-457187-342046.html
注意:
1.如果出现错误:
Undefined function or method 'ConstructKernelMatrix' for input arguments of type 'double'
可以把ConstructKernelMatrix函数另存为一个*.m文件试试。
2.若有其他错误,请回复,共同解决。
- 核主成份分析的MatLab实现
- 核主成份分析的MatLab实现(转)
- PCA:matlab主成份分析
- C++与matlab混合编程基于主成份分析算法的数值分析(一)
- C++与matlab混合编程基于主成份分析算法的数值分析(二)
- PCA-主成份分析实现流程
- 主成份分析
- PCA(主成份分析)
- R主成份分析
- 主成份分析
- 主成份分析图解
- PCA(主成份分析法)技术及其Python实现
- 主成份分析(PCA)——原理、实现步骤
- 主成份分析与因子分析
- 独立主成份分析(ICA)
- 主成份分析算法原理
- 主成份分析(PCA)详解
- [图像] PCA主成份分析
- WINDOWS 2008 R2无法安装显卡驱动问题解决
- ibatis 入门
- StringBuffer类
- C语言中的.c 和.h 文件区别
- 笔记2
- 核主成份分析的MatLab实现
- 川大校赛总结
- 职业生涯规划
- Linux目录结构
- UVa 10617 - Again Palindrome 字符串dp
- buttongroup控件的SelectionChangeFcn的用法
- Gridview用法大总结(牛年珍藏版)
- Deep learning:四(logistic regression练习)
- 汇编程序结构