目录
粒子群算法
简介
粒子群算法(Particle Swarm Optimization,简称PSO)是1995年Eberhart博士和Kennedy博士一起提出的[1]。粒子群算法是通过模拟鸟群捕食行为设计的一种群智能算法。区域内有大大小小不同的食物源，鸟群的任务是找到最大的食物源（全局最优解），鸟群的任务是找到这个食物源。鸟群在整个搜寻的过程中，通过相互传递各自位置的信息，让其他的鸟知道食物源的位置最终，整个鸟群都能聚集在食物源周围，即我们所说的找到了最优解，问题收敛。学者受自然界的启发开发了诸多类似智能算法，如蚁群算法、布谷鸟搜索算法、鱼群算法、捕猎算法等等。
算法策略
粒子群算法的目标是使所有粒子在多维超体（multi-dimensional hyper-volume）中找到最优解。首先给空间中的所有粒子分配初始随机位置和初始随机速度。然后根据每个粒子的速度、问题空间中已知的最优全局位置和粒子已知的最优位置依次推进每个粒子的位置。随着计算的推移，通过探索和利用搜索空间中已知的有利位置，粒子围绕一个或多个最优点聚集或聚合。该算法设计玄妙之处在于它保留了最优全局位置和粒子已知的最优位置两个信息。后续的实验发现，保留这两个信息对于较快收敛速度以及避免过早陷入局部最优解都具有较好的效果。这也奠定了后续粒子群算法改进方向的基础。
算法步骤
基础粒子群算法步骤较为简单。粒子群优化算法是由一组粒子在搜索空间中运动，受其自身的最佳过去位置pbest和整个群或近邻的最佳过去位置gbest的影响。每次迭代粒子i的第d维速度更新公式为:
每个粒子的速度和位置都以随机方式进行初始化。而后粒子就朝着全局最优和个体最优的方向靠近。
参数解释：
有关c1，c2的补充：
-
低的值允许粒子在被拉回前可以在目标区域外徘徊，高的值则导致粒子突然冲向或超越目标区域；
-
若c1 = 0，则粒子没有认知能力，在粒子的相互作用下，容易陷入局部极值点；
-
若c2 = 0，粒子间没有社会信息共享，算法变成一个多起点的随机搜索；
-
若c1 = c2 = 0，粒子将一直以当前速度飞行，直到到达边界。通常c1，c2在[0，4]之间，一般取c1 = c2 = 2
粒子群算法的流程图如下：
支持向量机
简介
支持向量机（support vector machines, SVM）是一种二分类模型，它的基本模型是定义在特征空间上的间隔最大的线性分类器，间隔最大使它有别于感知机；SVM还包括核技巧，这使它成为实质上的非线性分类器。SVM的的学习策略就是间隔最大化，可形式化为一个求解凸二次规划的问题，也等价于正则化的合页损失函数的最小化问题。SVM的的学习算法就是求解凸二次规划的最优化算法。
原理
SVM是一种在特征空间中寻找最大间隔超平面的判别式、线性、二分类器。我们之所以去寻找最大间隔的超平面就是为了追求模型的最大鲁棒性，对未知数据有最强的泛化能力。例如我们对样本施加一些噪声，间隔越大，就越不容易被错误分到另一类，如下图所示 ，a图为所有的数据点，b图为支持向量机最终的最大间隔超平面，可以看到它对噪声有着很强的鲁棒性。但是c图的分类效果的鲁棒性很差。
SVM的基本形式：
公式二是是约束条件，也就是“函数距离”，它确保所有正负样本都被正确分类。
由于(2)约束条件的存在，导致目标函数(1)的定义域可能奇形怪状，不一定是凸集，不能直接用导数为0求极值。所以我们需要用拉格朗日算子将原来问题转化成对偶问题后，约束条件被整合到下式(3)中，这个新目标函数不存在约束条件，可以直接用凸优化求解。
分别求拉式函数 关于 和 的偏导数，化为：
可以看到(4)的最后两项是两两样本内积<xi,xj>，这与后面提到的核函数形式一致，方便直接引入核函数的。
当然，拉式函数有解的充要条件是KKT条件：
核函数
有些问题比较付下，在低维空间线性不可分的样本，映射到高维空间后就线性可分了 ，如下图：
为此我们引入核函数，核函数并不直接对样本升维，它提供了一种计算高维空间两向量内积的方法，在低维空间通过计算核函数就可以得出了高维空间两向量内积。
该方法可以非常丝滑地嵌入公式(4)：
其中的 就是核函数，常见的核函数：
线性核：
线性核等价于原生SVM，直接计算内积；
选择线性核的场景：
-
数据线性可分
-
样本量特别大（此时再引入非线性核，计算量太大）
多项式核：
多项式核将数据映射到d维空间
选择多项式核的场景：
-
样本线性不可分时
-
特征较少时
高斯核：
选择高斯核的场景：
-
样本线性不可分时
-
特征较少时
粒子群优化算法的优化
粒子群算法的速度迭代公式中有 3 个重要参数：惯性权重 w 和学习因子 c1、c2。惯性权重决定着粒子先前的飞行速度对于当前的飞行速度的影响程度，惯性权重的选取对于平衡算法的全局搜索和局部搜索能力有重大意义。在迭代过程中，要考虑到算法的全局性和局部性，采取合适的惯性权重来进行搜索。构建改进公式的原因是想利用改进后的幂指函数算子，将其加入到惯性权重中，以总的迭代次数为基础，在 D 维空间中搜索时让每个粒子扩大搜索空间，以此增大种群多样性。算法在进化过程中，w 的值动态变化调整，算法在进化初期搜索能力较强，能在较大范围的解空间内搜索，解空间是所有可能的解构成的集合，并不断搜索到新的区域，然后到后期逐渐收敛到较好的区域并进行更精细的搜索，以加快收敛速度。表达式为：
上式中：g 为当前的迭代次数；a、b、d 为参数；a 为惯性权重的最大值值，b为惯性权重的最小值，d 为参数，取值为 2；gmax 代表迭代次数的最大值。
作为算法运行中的重要式子，c1 是粒子具有自我学习总结能力的体现，代表的是个体自我认知，学习因子 c2 是粒子具有向表现更好的粒子学习的本领体现，代表的是粒子的社会认知，需要设置恰当的 c1、c2 以便于粒子进行搜寻。初期要关注个体自我认识的能力，后期则应注重个体获取社会信息的能力。由此设置学习因子为：
式中：g 为当前迭代次数；gmax 为最大迭代次数；e1、e2、f1、f2 为参数，e1 取值 1，e2 取值 2，f1 取值 1.，f2 取值 2。
改进后新的粒子速度和位置迭代公式为：
与基本粒子群算法相比较，对基于惯性权重和学习因子动态调整的粒子群算法的速度公式作出调整，做到了随着迭代次数变化而动态变化，其全局搜索能力有效提高，粒子的收敛性也得到了加强。然而当遇到一些复杂的多峰函数曲线的寻优问题时，此改进的粒子群算法也存在收敛速度及精度达不到要求等问题。
基于粒子群优化算法的SVM的实现
本方法通过对人脸库人脸进行特征提取，同时通过改进的PSO算法对 SVM 参数进行优化后将所提取的特征用于分类识别。由于支持向量机的分类准确率在一定程度上会受到惩罚因子 C 与高斯核函数参数 y 的作用，因此把 C 与 y 作为粒子，通过改进的 PSO 算法对参数 C 与 y 进行优化，来增强支持向量机对于人脸识别的准确程度。
本程序设置粒子群体个体数目为30，迭代代数为300，设置速度更新的惯性权重w的最大值为0.9，最小值为0.2，使得pso早期搜索具有大的搜索范围，后期减少位移，加快收敛速度。所有的程序如下所示：
%///
%main.m 优化算法的主程序
%///
clc
clear all
[xm,fv,jilu,Pbest] = PSO_adaptation(@AdaptFunc,30,0.9,0.2,300,1000)
%function[xm,fv] = PSO_adaptation(fitness,N,wmax,wmin,M,D)
% wmax：惯性权重最大值
% wmin：惯性权重最小值
% M：最大迭代次数
% D：搜索空间维数
% N：初始化群体个体数目
%///
%AdaptFunc.m 适应度函数的程序
%///
function y = AdaptFunc(x);
x(1)=x(1)/1000;
train(x(1),x(2));
y=test();
end
%///
%PSO_adaptation.m 粒子群优化算法的程序
%///
function[xm,fv,jilu,Pbest] = PSO_adaptation(fitness,N,wmax,wmin,M,flag)
% c1,c2：学习因子
% wmax：惯性权重最大值
% wmin：惯性权重最小值
% M：最大迭代次数
% D：搜索空间维数
% N：初始化群体个体数目
popcmax=10^(6);
popcmin=10;
popgmax=10^(3);
popgmin=10^(-4);
k = 0.6; % k belongs to [0.1,1.0];
Vcmax = k*popcmax;
Vcmin = -Vcmax ;
Vgmax = k*popgmax;
Vgmin = -Vgmax ;
% 初始化种群的个体（可以在这里限定位置和速度的范围）
for i = 1:N
x(i,1) = (popcmax-popcmin)*rand+popcmin; % 随机初始化位置
x(i,2) = (popgmax-popgmin)*rand+popgmin; % 随机初始化位置
v(i,1) = Vcmax*rands(1); % 随即初始化速度
v(i,2) = Vgmax*rands(1); % 随即初始化速度
end
%先计算各个粒子的适应度，并初始化个体最优解pi和整体最优解pg %
%初始化pi %
display("初始参数计算")
for i = 1:N
p(i) = fitness(x(i,:));
y(i,:) = x(i,:);
end
[pg_fit,pg_indx]=max(p);
pg = x(pg_indx,:);
%主循环函数，进行迭代，直到达到精度的要求 %
for t = 1:M
display(['优化代数---第',num2str(t), '代']);
display('......................................................................................');
if t==1
fv=p;
end
w=wmax+(wmin-wmax)*t/M;
c1=2-t/M;
c2=1+t/M;
for i = 1:N %更新函数,其中v是速度向量，x为位置,i为迭代特征
display(['优化第',num2str(t), '代的第',num2str(i), '个粒子']);
v(i,:) = w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg-x(i,:));
if v(i,1) > Vcmax
v(i,1) = Vcmax;
end
if v(i,1) < Vcmin
v(i,1) = Vcmin;
end
if v(i,2) > Vgmax
v(i,2) = Vgmax;
end
if v(i,2) < Vgmin
v(i,2) = Vgmin;
end
x(i,:) = x(i,:)+v(i,:);
if x(i,1) > popcmax
x(i,1) = popcmax-rand*10;
end
if x(i,1) < popcmin
x(i,1) = popcmin+rand;
end
if x(i,2) > popgmax
x(i,2) = popgmax-rand*10;
end
if x(i,2) < popgmin
x(i,2) = popgmin+rand;
end
fv(i)=fitness(x(i,:));
if fv(i) > p(i)
p(i) = fv(i) ;
y(i,:) = x(i,:) ;
end
end
jilu(1:30,2*t-1:2*t)=x;
[pg_fit,pg_indx]=max(fv);
pg=y(pg_indx,:);
Pbest(t) =pg_fit;
end
%给出最后的计算结果 %
xm = pg';
xm(1)=xm(1)*flag;
fv = Pbest(M);
for t=1:M
plot(t,Pbest(t),'ro');
hold on
end
xlabel('进化次数') ;
ylabel('适应度值') ;
%///
%test.m 测试对于整个测试集的识别率，输出：accuracy --- 对于测试集合的识别率
%///
function y=test()
display('测试开始...');
nFacesPerPerson = 5;
nPerson = 40;
bTest = 1;
% 读入测试集合
[imgRow,imgCol,TestFace,testLabel] = ReadFaces(nFacesPerPerson, nPerson, bTest);
% 读入相关训练结果
load('Mat/PCA.mat');
load('Mat/scaling.mat');
load('Mat/trainData.mat');
load('Mat/multiSVMTrain.mat');
% PCA降维
[m n] = size(TestFace);
TestFace = (TestFace-repmat(meanVec, m, 1))*V; % 经过pca变换降维
TestFace = scaling(TestFace,1,A0,B0);
% 多类 SVM 分类
classes = multiSVMClassify(TestFace);
% 计算识别率
nError = sum(classes ~= testLabel);
accuracy = 1 - nError/length(testLabel);
y=accuracy*100;
rate=num2str(accuracy*100);
display(['对于测试集200个人脸样本的识别率为',rate, '%']);
display('..............................');
%///
%test.m
%整个训练过程，包括读入图像，PCA降维以及多类 SVM 训练，各个阶段的处理结果分别保存至文件：
% 将 PCA 变换矩阵 W 保存至 Mat/PCA.mat
% 将 scaling 的各维上、下界信息保存至 Mat/scaling.mat
% 将 PCA 降维并且 scaling 后的数据保存至 Mat/trainData.mat
% 将多类 SVM 的训练信息保存至 Mat/multiSVMTrain.mat
%///
global imgRow;
global imgCol;
display(' ');
display(' ');
display('训练开始...');
nPerson=40;
nFacesPerPerson = 5;
[imgRow,imgCol,FaceContainer,faceLabel]=ReadFaces(nFacesPerPerson,nPerson);
nFaces=size(FaceContainer,1);%样本（人脸）数目
[pcaFaces, W] = fastPCA(FaceContainer, 20);% 主成分分析PCA
% pcaFaces是200*20的矩阵, 每一行代表一张主成分脸(共40人，每人5张)，每个脸20个维特征
% W是分离变换矩阵, 10304*20 的矩阵
X = pcaFaces;
[X,A0,B0] = scaling(X);
save('Mat/scaling.mat', 'A0', 'B0');
% 保存 scaling 后的训练数据至 trainData.mat
TrainData = X;
trainLabel = faceLabel;
W = W;
save('Mat/trainData.mat', 'TrainData', 'trainLabel','W');
for iPerson = 1:nPerson
nSplPerClass(iPerson) = sum( (trainLabel == iPerson) );
end
multiSVMStruct = multiSVMTrain(TrainData, nSplPerClass, nPerson, C, gamma);
save('Mat/multiSVMTrain.mat', 'multiSVMStruct');
display('训练结束。');
display('..............................');
%///
%scaling.m 对人脸图像进行归一化处理函数
%///
function [SVFM, lowVec, upVec] = scaling(VecFeaMat, bTest, lRealBVec, uRealBVec)
if nargin < 2
bTest = 0;
end
% 缩放目标范围[-1, 1]
lTargB = -1;
uTargB = 1;
[m n] = size(VecFeaMat);
SVFM = zeros(m, n);
if bTest
if nargin < 4
error('To do scaling on testset, param lRealB and uRealB are needed.');
end
if nargout > 1
error('When do scaling on testset, only one output is supported.');
end
for iCol = 1:n
if uRealBVec(iCol) == lRealBVec(iCol)
SVFM(:, iCol) = uRealBVec(iCol);
SVFM(:, iCol) = 0;
else
SVFM(:, iCol) = lTargB + ( VecFeaMat(:, iCol) - lRealBVec(iCol) ) / ( uRealBVec(iCol)-lRealBVec(iCol) ) * (uTargB-lTargB); % 测试数据的scaling
end
end
else
upVec = zeros(1, n);
lowVec = zeros(1, n);
for iCol = 1:n
lowVec(iCol) = min( VecFeaMat(:, iCol) );
upVec(iCol) = max( VecFeaMat(:, iCol) );
if upVec(iCol) == lowVec(iCol)
SVFM(:, iCol) = upVec(iCol);
SVFM(:, iCol) = 0;
else
SVFM(:, iCol) = lTargB + ( VecFeaMat(:, iCol) - lowVec(iCol) ) / ( upVec(iCol)-lowVec(iCol) ) * (uTargB-lTargB); % 训练数据的scaling
end
end
end
%///
%ReadFaces.m 用于对人脸数据进行读取的函数
%///
function [imgRow,imgCol,FaceContainer,faceLabel]=ReadFaces(nFacesPerPerson, nPerson, bTest)
if nargin==0 %default value
nFacesPerPerson=5;%前5张用于训练
nPerson=40;%要读入的人数（每人共10张，前5张用于训练）
bTest = 0;
elseif nargin < 3
bTest = 0;
end
img=imread('Data/ORL/S1/1.pgm');%为计算尺寸先读入一张
[imgRow,imgCol]=size(img);
FaceContainer = zeros(nFacesPerPerson*nPerson, imgRow*imgCol);
faceLabel = zeros(nFacesPerPerson*nPerson, 1);
% 读入训练数据
for i=1:nPerson
i1=mod(i,10); % 个位
i0=char(i/10);
strPath='Data/ORL/S';
if( i0~=0 )
strPath=strcat(strPath,'0'+i0);
end
strPath=strcat(strPath,'0'+i1);
strPath=strcat(strPath,'/');
tempStrPath=strPath;
for j=1:nFacesPerPerson
strPath=tempStrPath;
if bTest == 0 % 读入训练数据
strPath = strcat(strPath, '0'+j);
else
strPath = strcat(strPath, num2str(5+j));
end
strPath=strcat(strPath,'.pgm');
img=imread(strPath);
%把读入的图像按列存储为行向量放入向量化人脸容器faceContainer的对应行中
FaceContainer((i-1)*nFacesPerPerson+j, :) = img(:)';
faceLabel((i-1)*nFacesPerPerson+j) = i;
end % j
end % i
% 保存人脸样本矩阵
save('Mat/FaceMat.mat', 'FaceContainer')
%///
%multiSVMTrain.m 用于对人脸数据进行读取的函数
% 采用1对1投票策略将 SVM 推广至多类问题的训练过程，将多类SVM训练结果保存至multiSVMStruct中
% 输入:--TrainData:每行是一个样本人脸
% --nClass:人数，即类别数
% --nSampPerClass:nClass*1维的向量，记录每类的样本数目，如 nSampPerClass(iClass)
% 给出了第iClass类的样本数目
% --C:错误代价系数，默认为 Inf
% --gamma:径向基核函数的参数 gamma，默认值为1
% 输出:--multiSVMStruct:一个包含多类SVM训练结果的结构体
%///
function multiSVMStruct = multiSVMTrain(TrainData, nSampPerClass, nClass, C, gamma)
% 默认参数
if nargin < 4
C = Inf;
gamma = 1;
elseif nargin < 5
gamma = 1;
end
%开始训练，需要计算每两类间的分类超平面，共(nClass-1)*nClass/2个
for ii=1:(nClass-1)
for jj=(ii+1):nClass
clear X;
clear Y;
startPosII = sum( nSampPerClass(1:ii-1) ) + 1;
endPosII = startPosII + nSampPerClass(ii) - 1;
X(1:nSampPerClass(ii), :) = TrainData(startPosII:endPosII, :);
startPosJJ = sum( nSampPerClass(1:jj-1) ) + 1;
endPosJJ = startPosJJ + nSampPerClass(jj) - 1;
X(nSampPerClass(ii)+1:nSampPerClass(ii)+nSampPerClass(jj), :) = TrainData(startPosJJ:endPosJJ, :);
% 设定两两分类时的类标签
Y = ones(nSampPerClass(ii) + nSampPerClass(jj), 1);
Y(nSampPerClass(ii)+1:nSampPerClass(ii)+nSampPerClass(jj)) = 0;
% 第ii个人和第jj个人两两分类时的分类器结构信息
CASVMStruct{ii}{jj}= svmtrain( X, Y, 'Kernel_Function', @(X,Y) kfun_rbf(X,Y,gamma), 'boxconstraint', C );
end
end
% 已学得的分类结果
multiSVMStruct.nClass = nClass;
multiSVMStruct.CASVMStruct = CASVMStruct;
% 保存参数
save('Mat/params.mat', 'C', 'gamma');
%///
%multiSVMClassify.m 用于对人脸数据进行读取的函数
% 采用1对1投票策略将 SVM 推广至多类问题的分类过程
% 输入:--TestFace:测试样本集。m*n 的2维矩阵，每行一个测试样本
% --multiSVMStruct:多类SVM的训练结果，由函数 multiSVMTrain 返回，默认是从Mat/multiSVMTrain.mat文件中读取
% 输出:--class: m*1 列向量，对应 TestFace 的类标签
% 读入训练结果
%///
function class = multiSVMClassify(TestFace, multiSVMStruct)
if nargin < 2
t = dir('Mat/multiSVMTrain.mat');
if length(t) == 0
error('没有找到训练结果文件，请在分类以前首先进行训练！');
end
load('Mat/multiSVMTrain.mat');
end
nClass = multiSVMStruct.nClass; % 读入类别数
CASVMStruct = multiSVMStruct.CASVMStruct; % 读入两两类之间的分类器信息
%%%%%%%%%%%%%%%%%%%%%%%% 投票策略解决多类问题 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m = size(TestFace, 1);
Voting = zeros(m, nClass); % m个测试样本，每个样本nPerson 个类别的投票箱
for iIndex = 1:nClass-1
for jIndex = iIndex+1:nClass
classes = svmclassify(CASVMStruct{iIndex}{jIndex}, TestFace);
% 投票
Voting(:, iIndex) = Voting(:, iIndex) + (classes == 1);
Voting(:, jIndex) = Voting(:, jIndex) + (classes == 0);
end % for jClass
end % for iClass
% final decision by voting result
[vecMaxVal, class] = max( Voting, [], 2 );
%display(sprintf('TestFace对应的类别是:%d',class));
%///
%kfun_rbf.m rbf核函数程序
%///
function K = kfun_rbf(U, V, gamma)
% rbf 核函数
[m1 n1] = size(U);
[m2 n2] = size(V);
K = zeros(m1, m2);
for ii = 1:m1
for jj = 1:m2
K(ii, jj) = exp( -gamma * norm(U(ii, :)-V(jj, :))^2 );
end
end
%///
%fastPCA.m 快速PCA程序
%///
function [pcaA V] = fastPCA( A, k )
[r c] = size(A);
% 样本均值
meanVec = mean(A);
% 计算协方差矩阵的转置 covMatT
Z = (A-repmat(meanVec, r, 1));
covMatT = Z * Z';
% 计算 covMatT 的前 k 个本征值和本征向量
[V D] = eigs(covMatT, k);
% 得到协方差矩阵 (covMatT)' 的本征向量
V = Z' * V;
% 本征向量归一化为单位本征向量
for i=1:k
V(:,i)=V(:,i)/norm(V(:,i));
end
% 线性变换（投影）降维至 k 维
pcaA = Z * V;
% 保存变换矩阵 V 和变换原点 meanVec
save('Mat/PCA.mat', 'V', 'meanVec');通过优化后的PSO优化算法得到的结果如下图所示：
适应度值为PSO每一代的参数训练出的SVM对人脸库的200张人脸进行识别的准确率，经过300代的迭代次数最后我们可以看到，收敛到83.5%的准确率。
第一代的各个粒子的的位置分布如下图所示：
从上图我们可以观察到，粒子的分布较分散，有较强的全局搜索能力。
经过300代后的粒子分布如下图所示：
从上图我们可以观察到，粒子的分布比较集中，而且还有少量的粒子具有一定的探索能力，防止算法陷入局部最优解。
861
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
