核主元分析 KPCA及matlab代码
f核主元分析 (Kernel Principal Component Analysis , KPCA) – MATLAB中文论坛 (ilovematlab.cn)发现一篇很好的讲解KPCA代码的文章KPCA简介KPCA、PCA应用:降维特征提取去噪故障检测假设训练数据由矩阵𝑋𝑛×𝑚表示,𝑚代表训练数据的大小,𝑛代表训练数据中特征向量的维度,KPCA具体实现过程如下:Step1:利用核函数
前言
发现一篇很好的分享KPCA代码的文章
转载链接:
f核主元分析 (Kernel Principal Component Analysis , KPCA) – MATLAB中文论坛 (ilovematlab.cn)
原博主处有打包matlab代码
对于简单的数据降维,可参考:
用matlab实现kpca(核主成分分析法)_Bbei的博客-CSDN博客_matlabkpca
kpca思路和pca很像,可以了解 MATLAB实例:PCA(主成成分分析)详解_FDA_sq的博客-CSDN博客_matlab pca
KPCA简介
KPCA、PCA应用:
- 降维
- 特征提取
- 去噪
- 故障检测
假设训练数据由矩阵𝑋𝑛×𝑚表示,𝑚代表训练数据的大小,𝑛代表训练数据中特征向量的维度,KPCA具体实现过程如下:
Step1:利用核函数计算核矩阵𝑲=[𝐾𝑖𝑗]𝑚×𝑚
Step2:核矩阵𝑲的中心化
Step3:对中心化核矩阵𝑲𝒄𝒆𝒏𝒕𝒆𝒓进行特征值分解
Step4:特征向量𝛼的标准化
Step5:提取训练数据的主成分
Step6:提取测试数据的主成分
KPCA的性能取决于核函数、核参数的选择。常用的核函数有线性核函数(Linear)、径向基(Radial Basis Function, RBF)和积分径向基(The Integrated Radial Basis Function, IRBF)核函数。当使用线性核函数时,KPCA退化为PCA。
主要功能
- 易于使用的 API
- 支持基于 KPCA 的数据降维、特征提取、数据重构
- 支持基于 KPCA 的故障检测和故障诊断
- 支持多种核函数 (linear, gaussian, polynomial, sigmoid, laplacian)
- 支持基于主元贡献率或给定数字的降维维度/主元个数选取
注意
- 仅支持基于高斯核函数的故障诊断
- 核函数的参数对KPCA 模型的性能影响很大
- 此代码仅供参考
使用说明
1. KPCA的建模过程(故障检测):
(1)获取训练数据(工业过程数据需要进行标准化处理)
(2)计算核矩阵
(3)核矩阵中心化
(4)特征值分解
(5)特征向量的标准化处理
(6)主元个数的选取
(7)计算非线性主成分(即降维结果或者特征提取结果)
(8)SPE和T2统计量的控制限计算
function model = kpca_train(X,options)
% DESCRIPTION
% Kernel principal component analysis (KPCA)
%
% mappedX = kpca_train(X,options)
%
% INPUT
% X Training samples (N*d)
% N: number of samples
% d: number of features
% options Parameters setting
%
% OUTPUT
% model KPCA model
%
%
% Created on 9th November, 2018, by Kepeng Qiu.
% number of training samples
L = size(X,1);
% Compute the kernel matrix
K = computeKM(X,X,options.sigma);
% Centralize the kernel matrix
unit = ones(L,L)/L;
K_c = K-unit*K-K*unit+unit*K*unit;
% Solve the eigenvalue problem
[V_s,D] = eigs(K_c/L);
lambda = diag(D);
% Normalize the eigenvalue
% V_s = V ./ sqrt(L*lambda)';
% Compute the numbers of principal component
% Extract the nonlinear component
if options.type == 1 % fault detection
dims = find(cumsum(lambda/sum(lambda)) >= 0.85,1, 'first');
else
dims = options.dims;
end
mappedX = K_c* V_s(:,1:dims) ;
% Store the results
model.mappedX = mappedX ;
model.V_s = V_s;
model.lambda = lambda;
model.K_c = K_c;
model.L = L;
model.dims = dims;
model.X = X;
model.K = K;
model.unit = unit;
model.sigma = options.sigma;
% Compute the threshold
model.beta = options.beta;% corresponding probabilities
[SPE_limit,T2_limit] = comtupeLimit(model);
model.SPE_limit = SPE_limit;
model.T2_limit = T2_limit;
end
2. KPCA的测试过程:
(1)获取测试数据(工业过程数据需要利用训练数据的均值和标准差进行标准化处理)
(2)计算核矩阵
(3)核矩阵中心化
(4)计算非线性主成分(即降维结果或者特征提取结果)
(5)SPE和T2统计量的计算
function [SPE,T2,mappedY] = kpca_test(model,Y)
% DESCRIPTION
% Compute the T2 statistic, SPE statistic,and the nonlinear component of Y
%
% [SPE,T2,mappedY] = kpca_test(model,Y)
%
% INPUT
% model KPCA model
% Y test data
%
% OUTPUT
% SPE the SPE statistic
% T2 the T2 statistic
% mappedY the nonlinear component of Y
%
% Created on 9th November, 2018, by Kepeng Qiu.
% Compute Hotelling's T2 statistic
% T2 = diag(model.mappedX/diag(model.lambda(1:model.dims))*model.mappedX');
% the number of test samples
L = size(Y,1);
% Compute the kernel matrix
Kt = computeKM(Y,model.X,model.sigma );
% Centralize the kernel matrix
unit = ones(L,model.L)/model.L;
Kt_c = Kt-unit*model.K-Kt*model.unit+unit*model.K*model.unit;
% Extract the nonlinear component
mappedY = Kt_c*model.V_s(:,1:model.dims);
% Compute Hotelling's T2 statistic
T2 = diag(mappedY/diag(model.lambda(1:model.dims))*mappedY');
% Compute the squared prediction error (SPE)
SPE = sum((Kt_c*model.V_s).^2,2)-sum(mappedY.^2 ,2);
end
3. demo1: 降维、特征提取
(1) 源代码
% Demo1: dimensionality reduction or feature extraction
% ---------------------------------------------------------------------%
clc
clear all
close all
addpath(genpath(pwd))
% 4 circles
load circledata
%
X = circledata;
for i = 1:4
scatter(X(1+250*(i-1):250*i,1),X(1+250*(i-1):250*i,2))
hold on
end
% Parameters setting
options.sigma = 5; % kernel width
options.dims = 2; % output dimension
options.type = 0; % 0:dimensionality reduction or feature extraction
% 1:fault detection
options.beta = 0.9; % corresponding probabilities (for ault detection)
options.cpc = 0.85; % Principal contribution rate (for ault detection)
% Train KPCA model
model = kpca_train(X,options);
figure
for i = 1:4
scatter(model.mappedX(1+250*(i-1):250*i,1), ...
model.mappedX(1+250*(i-1):250*i,2))
hold on
end
(2)结果 (分别为原图和特征提取后的图)
4. demo2: 故障检测(需要调节核宽度、主元贡献率和置信度等参数来提高故障检测效果)
(1)源代码
% Demo2: Fault detection
% X: training samples
% Y: test samples
% Improve the performance of fault detection by adjusting parameters
% 1. options.sigma = 16; % kernel width
% 2. options.beta % corresponding probabilities
% 3. options.cpc ; % principal contribution rate
% ---------------------------------------------------------------------%
clc
clear all
close all
addpath(genpath(pwd))
%
X = rand(200,10);
Y = rand(100,10);
Y(20:40,:) = rand(21,10)+3;
Y(60:80,:) = rand(21,10)*3;
% Normalization (if necessary)
% mu = mean(X);
% st = std(X);
% X = zscore(X);
% Y = bsxfun(@rdivide,bsxfun(@minus,Y,mu),st);
% Parameters setting
options.sigma = 16; % kernel width
options.dims = 2; % output dimension
options.type = 1; % 0:dimensionality reduction or feature extraction
% 1:fault detection
options.beta = 0.9; % corresponding probabilities (for ault detection)
options.cpc = 0.85; % principal contribution rate (for ault detection)
% Train KPCA model
model = kpca_train(X,options);
% Test a new sample Y (vector of matrix)
[SPE,T2,mappedY] = kpca_test(model,Y);
% Plot the result
plotResult(model.SPE_limit,SPE);
plotResult(model.T2_limit,T2);
(2)结果(分别是SPE统计量和T2统计量的结果图)
更新:
01 核函数
类 Kernel 用于计算核函数矩阵:
%{
type -
linear : k(x,y) = x'*y
polynomial : k(x,y) = (γ*x'*y+c)^d
gaussian : k(x,y) = exp(-γ*||x-y||^2)
sigmoid : k(x,y) = tanh(γ*x'*y+c)
laplacian : k(x,y) = exp(-γ*||x-y||)
degree - d
offset - c
gamma - γ
%}
kernel = Kernel('type', 'gaussian', 'gamma', value);
kernel = Kernel('type', 'polynomial', 'degree', value);
kernel = Kernel('type', 'linear');
kernel = Kernel('type', 'sigmoid', 'gamma', value);
kernel = Kernel('type', 'laplacian', 'gamma', value);
例如,计算 X 和 Y 的高斯核函数矩阵:
X = rand(5, 2);
Y = rand(3, 2);
kernel = Kernel('type', 'gaussian', 'gamma', 2);
kernelMatrix = kernel.computeMatrix(X, Y);
>> kernelMatrix
kernelMatrix =
0.5684 0.5607 0.4007
0.4651 0.8383 0.5091
0.8392 0.7116 0.9834
0.4731 0.8816 0.8052
0.5034 0.9807 0.7274
02. 数据降维
以 helix 数据为例:
clc
clear all
close all
addpath(genpath(pwd))
load('.\data\helix.mat', 'data')
kernel = Kernel('type', 'gaussian', 'gamma', 2);
parameter = struct('numComponents', 2, ...
'kernelFunc', kernel);
% build a KPCA object
kpca = KernelPCA(parameter);
% train KPCA model
kpca.train(data);
% mapping data
mappingData = kpca.score;
% Visualization
kplot = KernelPCAVisualization();
% visulize the mapping data
kplot.score(kpca)
训练结果为:
*** KPCA model training finished ***
running time = 0.2798 seconds
kernel function = gaussian
number of samples = 1000
number of features = 3
number of components = 2
number of T2 alarm = 135
number of SPE alarm = 0
accuracy of T2 = 86.5000%
accuracy of SPE = 100.0000%
香蕉型数据的降维结果:
03. 数据重构
同心圆数据的重构结果:
clc
clear all
close all
addpath(genpath(pwd))
load('.\data\circle.mat', 'data')
kernel = Kernel('type', 'gaussian', 'gamma', 0.2);
parameter = struct('numComponents', 2, ...
'kernelFunc', kernel);
% build a KPCA object
kpca = KernelPCA(parameter);
% train KPCA model
kpca.train(data);
% reconstructed data
reconstructedData = kpca.newData;
% Visualization
kplot = KernelPCAVisualization();
kplot.reconstruction(kpca)
04. 主元个数的确定
提供了两种方式来确定主元个数:主元贡献率和给定的数量
第一种:主元贡献率 (取值 0~1 之间)
通过键值对 'numComponents' - value 来设置。比如,将主元贡献率设置为 0.75:
clc
clear all
close all
addpath(genpath(pwd))
load('.\data\TE.mat', 'trainData')
kernel = Kernel('type', 'gaussian', 'gamma', 1/128^2);
parameter = struct('numComponents', 0.75, ...
'kernelFunc', kernel);
% build a KPCA object
kpca = KernelPCA(parameter);
% train KPCA model
kpca.train(trainData);
% Visualization
kplot = KernelPCAVisualization();
kplot.cumContribution(kpca)
从上图可以看出,当主元个数为 21 的时候,主元贡献率为 75.2656%, 超过了给定的阈值(0.75)。
第二种:给定的数目
同样通过键值对 'numComponents' - value 来设置。比如,将主元个数设置为 24:
clc
clear all
close all
addpath(genpath(pwd))
load('.\data\TE.mat', 'trainData')
kernel = Kernel('type', 'gaussian', 'gamma', 1/128^2);
parameter = struct('numComponents', 24, ...
'kernelFunc', kernel);
% build a KPCA object
kpca = KernelPCA(parameter);
% train KPCA model
kpca.train(trainData);
% Visualization
kplot = KernelPCAVisualization();
kplot.cumContribution(kpca)
从上图可以看出,当主元个数为 24 的时候,主元贡献率为 80.2539%。
05. 故障检测
数据来源于 TE 过程:
clc
clear all
close all
addpath(genpath(pwd))
load('.\data\TE.mat', 'trainData', 'testData')
kernel = Kernel('type', 'gaussian', 'gamma', 1/128^2);
parameter = struct('numComponents', 0.65, ...
'kernelFunc', kernel);
% build a KPCA object
kpca = KernelPCA(parameter);
% train KPCA model
kpca.train(trainData);
% test KPCA model
results = kpca.test(testData);
% Visualization
kplot = KernelPCAVisualization();
kplot.cumContribution(kpca)
kplot.trainResults(kpca)
kplot.testResults(kpca, results)
训练的结果为
*** KPCA model training finished ***
running time = 0.0986 seconds
kernel function = gaussian
number of samples = 500
number of features = 52
number of components = 16
number of T2 alarm = 16
number of SPE alarm = 17
accuracy of T2 = 96.8000%
accuracy of SPE = 96.6000%
测试的结果为
*** KPCA model test finished ***
running time = 0.0312 seconds
number of test data = 960
number of T2 alarm = 799
number of SPE alarm = 851
06. 故障诊断
故障诊断通过键值对 'diganosis' - value 来进行设置。
比如对某个采样点(假设为500)出现的故障进行诊断,则: 'diagnosis', [500, 500]
比如对某个时间段(假设为300-500)出现的故障进行诊断,则: 'diagnosis', [300, 500]
具体代码为
clc
clear all
close all
addpath(genpath(pwd))
load('.\data\TE.mat', 'trainData', 'testData')
kernel = Kernel('type', 'gaussian', 'gamma', 1/128^2);
parameter = struct('numComponents', 0.65, ...
'kernelFunc', kernel,...
'diagnosis', [300, 500]);
% build a KPCA object
kpca = KernelPCA(parameter);
% train KPCA model
kpca.train(trainData);
% test KPCA model
results = kpca.test(testData);
% Visualization
kplot = KernelPCAVisualization();
kplot.cumContribution(kpca)
kplot.trainResults(kpca)
kplot.testResults(kpca, results)
kplot.diagnosis(results)
故障诊断结果为
*** Fault diagnosis ***
Fault diagnosis start...
Fault diagnosis finished.
running time = 18.2738 seconds
start point = 300
ending point = 500
fault variables (T2) = 44 1 4
fault variables (SPE) = 1 44 18
可以看出,在300-500时间段,给出的故障变量参考为: 44,1(T2); 1,44(SPE)。
![Logo](https://devpress.csdnimg.cn/79de2bf0b7994defa4242ef90d5513fa.jpg)
开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!
更多推荐
所有评论(0)