EigenFace算法回顾及Matlab代码基于PCA的EigenFace算法发表自1987年，是第一种可行的人脸辨识算法。虽然已有20余年历史，但仍是人脸辨识算法研究中的经典，新算法都要与之作比较。
EigenFace是2D辨识算法，但为了进行3D表情辨识的研究，有必要对这一经典算法进行回顾，SIGGRAPH13的文献Online Modeling For Realtime Facial Animation实现表情3D重构的基础是SIGGRAPH99中A Morphable Model for the Synthesis of 3D Faces 提出的Morphable Facial Model，而建立这一模型的基础思想仍是PCA，与Eigenface有着天然联系。
学习EigenFace应该是研究生时代的事儿了，旧编重拾、开卷有益，并写了Matlab 代码附录于后。网上许多实例代码只实现了辨识，略去了一个重要环节:通过分解-重构，将一幅输入人脸照片表示为EigenFace基底的组合，这对于表情辨识及3D人脸模型分析都是很有用的(Online Modeling For Realtime Facial Animation中就利用了此思想)，附录的Matlab代码做了这一步。
训练集合包含20幅图片
生成20-1=19个特征脸
最近欧式距离法得出的前三位匹配
利用特征脸空间进行人脸重构
https://www.doczj.com/doc/8818330785.html,/~sis26/Eigenface Tutorial.htm中给出了EigenFace算法Matlab 代码，含有重构过程与比较，但其代码中有一个错误，训练集合照片未减去平均脸，就计算协方差矩阵了
具体算法文献教材上都有，捡要点写几句：
(1)将训练集合中的每幅图像拉伸为列向量，并减去所有图像的均值(称之为平均脸)，形成N*M矩阵A,其中N为单幅图像像素数，M为图像数目(训练集合容量)。
(2)求协方差矩阵AA'的特征向量，作为正交基底张成人脸空间，好是好但运算量过大，转而求替代矩阵(surrogate)A'A的特征向量，减少计算量
(3)矩阵A'A的秩等于M-1，这是由于减去平均脸所致，故有M-1个非零特征值(正)，去除属于0的特征向量，将M-1个属于非零特征值得特征向量(记住须作左乘A的修正)作为EigenFace基底(特征脸)，张成人脸空间。
(4)EigenFace基底由M-1个相互正交的向量构成，它们是协方差矩阵AA'的前M-1个最显著的特征向量方向，能量主要集中在这些向量方向上，但要记住虽正交但不完备，故存在重构误差。
(5)为了用EigenFace基底对人脸照片进行正确的分解-重构，需要对所得的基底向量进行规一化修正，因为A'A的特征向量左乘A之后，虽成为AA'的特征向量，但模不为1，需除以自身的模，修正为标准正交向量集合，才能进行投影分解-重构。
(6)训练集合及测试集合中的人脸照片都能利用EigenFace基底较好地实现分解-重构，但训练集合之外的人脸，重构误差变大
(7)EigenFace缺点：1拍摄时光照环境对识别效果(EigenFace基底)影响大
2训练集合扩容时，需重构EigenFace基底
(8) 为了有效显示EigenFace基底图像(特征脸)需要用imagesc函数
Eigenface算法识别人脸的步骤：
This section gives step-by-step instructions along with photos and formulas on how to recognize faces and implemented into Matlab. All the necessary files to complete this tutorial would be provided.
Steps
1． The first step is to obtain a set S with M face images. In our example M = 25 as shown at the beginning of the tutorial. Each image is transformed into a vector of size N and placed into the set.
2． After you have obtained your set, you will obtain the mean image Ψ
3． Then you will find the difference Φ between the input image and the mean image
4． Next we seek a set of M orthonormal vectors, u n , which best describes the distribution of the data. The k th vector, u k , is chosen such that
求特征值：
is a maximum, subject to
Note: u k and λk are the eigenvectors and eigenvalues of the covariance matrix C
5． We obtain the covariance matrix C in the following manner
Where the matrix
. The matrix C, however, is N 2
by N 2, and determining the N 2 eigenvectors and eigenvalues is an intractable task for typical image sizes. We need a computationally feasible method to find these eigenvectors.
We can solve for the N2 dimensional eigenvectors in this case by first solving for the eigenvectors of an M by M matrix- e.g., solving a 16*16 matrix rather than a 16,384*16,384 matrix and then taking appropriate linear combinations of the face images Φi ,Consider the eigenvectors v i of A T A such that A T Av i =μi v i Premultiplying both sides by A, we have
AA T Av i =μi Av i
From which we see that Av i (特征向量)are the eigenvectors of C=AA T .
{}
M S ΓΓΓΓ=, ......... , , , 321∑
=Γ=
ψM
n n
M 1
1ψ
-Γ=Φi i ()
∑=Φ=
M n n
T k
k u
M
1
2
1λotherwise if
01
k l u u lk k T l =???==δT
M n T
n n AA M C 11
=ΦΦ=
∑=}
{ , ....... , , , 321n A ΦΦΦΦ=
Following this analysis, we construct the M by M matrix L =A T A , where
(Φ是1*16,384维的向量), and find the M eigenvectors, v i (M 维的向
量)of L . These vectors determine linear combinations of the M training set face
images to form the eigenfaces u l (特征脸向量). Once we have found the eigenvectors, v l , u l
l =1,2,…M
6． These are the eigenfaces of our set of original images
Recognition Procedure
1. A new face is transformed into its eigenface components. First we compare our input image with our mean image and multiply their difference with each eigenvector of the L matrix. Each value would represent a weight and would be saved on a vector ?.
新脸向量在每个特征脸向量上的投影：
For k =1,2,…,M ’.
2. We now determine which face class provides the best description for the input
image. This is done by minimizing the Euclidean distance
3. The input face is consider to belong to a class if εk is bellow an established
threshold θε. Then the face image is considered to be a known face. If the difference is above the given threshold, but bellow a second threshold, the image can be
determined as a unknown face. If the input image is above these two thresholds, the image is determined NOT to be a face.
4. If the image is found to be an unknown face, you could decide whether or not you want to add the image to your training set for future recognitions. You would have to repeat steps 1 trough 7 to incorporate this new face image. The source code face recognition using Matlab is provided below:
Matlab 代码：
所用训练集合含20幅jpg 图片，存于名为TrainDataBase 的文件夹中，测试集合含10幅jpg 图片，存于名为TestDataBase 的文件夹中，命名均为1.jpg 2.jpg 3.jpg........ %---------------------------------------------------------------------- %人脸辨识：Face Recognition %EigenFace 算法
n
T m mn L ΦΦ=()ψ-Γ=T
k k u ω[]M T
ωωω , ....... , ,21=Ω2
k
k Ω-Ω=ε
%EigenFace空间的生成分解与合成
% 人脸辨识
%-----------------------------------------------------------------------
function meigenface()
%------------------------------------------------------
%step1:读取训练集合中的人脸照片，并通过列向量拉伸生成
% 人脸训练集合矩阵
%------------------------------------------------------
clear all
clc
close all
imset=[];
M=20;%训练集合中人脸照片的数目
figure(1);
for i=1:M
filename=['.\TrainDataBase\' num2str(i) '.jpg'];
im=imread(filename,'JPG');
im=im(:,:,2);
% im=im(:,:,3);
[imheight imwidth]=size(im);
%将图像矩阵按列拉伸为列向量
imv=reshape(im,imheight*imwidth,1); %将im中数据转化为imheight*imwidth高的列向量imset=[imset imv]; % 将每个图像的像素占据一列位置，共20列向量
if i==4
title('训练集合人脸照片','fontsize',18)
end
subplot(5,4,i);
imshow(im);
end
disp('============Output image is completed!===================')
% imset=double(imset);
imset=double(imset);
[w,h]=size(imset)
%------------------------------------------------------
%step2:生成训练集合的平均脸
%------------------------------------------------------
averageface=mean(imset,2);
% imset每行的平均值，由于imset为单列向量，每行包含一个元素，因此此操作为将每行的所有值取平均值
figure(2);
imshow(uint8(reshape(averageface,imheight,imwidth)));
title('平均脸');
%------------------------------------------------------
%step3:将人脸训练集合修正为与平均人脸的差:以便计算
% 协方差矩阵Covariance matrix
%------------------------------------------------------
imset_d=imset;
figure(3);
for i=1:20
imset_d(:,i)=imset(:,i)- averageface;
im=reshape(imset_d(:,i),imheight,imwidth);
if i==4
title('与平均脸的差','fontsize',18)
end
% subplot(4,5,i);
subplot(5,4,i);
im=imagesc(im); colormap('gray');
% imshow(im);
end
%
%------------------------------------------------------
%step4:生成EigenFace特征脸空间
%------------------------------------------------------
%互相关矩阵为A*A',为减少计算量，计算代理矩阵A'*A的特征值特征向量
A=double(imset_d); % include the data of twenty images. % matrix of 36000 hight and 20 width
[Ah,Aw]=size(A)
L=A'*A; % matrix of 20 hight and 20 width
[Lh,Lw]=size(L)
[v d]=eig(L) % v为L的特征向量，d为L的特征值
dL=L*v-v*d
%人脸照片修正为与平均脸的差，此时矩阵A'*A的秩为M-1,特征值中含有一个0
%需要去除属于0的特征向量,Eigenface基底由M-1组正交向量组成
[v d]= EliminateZero(v,d) %去除属于0的特征向量，同时将特征向量按特征值大小降序排列
%生成人脸空间矩阵:eigenface
eigenface=[];
for i=1:M-1
mv=A*v(:,i);%左乘矩阵A，从而修正为协方差矩阵AA'的特征向量
%对Eigenface向量进行规一化修正，使每个Eigenface向量模为1，便于进行分解与合成%是否进行此归一化修正对识别没有影响，但对利用Eigenface作为基底对人脸照片进行
%分解与重构却是必要的
mv=mv/norm(mv); %!!!!!!!!!!
eigenface=[eigenface mv]; % 特征人脸空间
[eh,ew]=size(eigenface)
end
%显示归一化之后的人脸空间eigenface
figure(4);
for i=1:M-1
im=eigenface(:,i);
im=reshape(im,imheight,imwidth);
if i==4
title('特征脸空间：EigenFace','fontsize',18)
end
subplot(5,4,i);
im=imagesc(im);colormap('gray');
%im=histeq(im,255);imshow(im);
end
%------------------------------------------------------
%step5: 计算训练集合中的照片在EigenFace基底上的投影坐标
% 以便进行人脸辨识
%注:如果用Eigenface基底直接对训练集合中的照片进行分解-重构
%则重构误差趋于0(1e-10)，如果对测试集合中的人脸(同一人，但表情姿态略有出入) %进行重构，则可观察到重构误差，这是由于EigenFace基底虽正交但不完备
%------------------------------------------------------
TrainCoeff=[];
for i=1:M
%计算训练集合中第i幅图像在EigenFace基底上的投影
disp('The %d th face image')
i=i
c=double(imset_d(:,i)')*eigenface;
TrainCoeff=[TrainCoeff c']
end
%------------------------------------------------------
%step6:人脸辨识，采用最近欧式距离法
%计算输入图像在EigenFace基底上的投影坐标
%与训练集合中各图像的投影坐标进行比较，选取欧式距离最小者作为匹配项
%------------------------------------------------------
%读入测试集合人脸照片
N=10;%测试集合人脸照片数目
figure(5);
testset=[];
for i=1:N
filename=['.\TestDataBase\' num2str(i) '.jpg'];
im=imread(filename,'JPG');
im=im(:,:,2);
[imheight imwidth]=size(im);
%将图像矩阵按列拉伸为列向量
imv=reshape(im,imheight*imwidth,1);
testset=[testset imv];
%
subplot(3,4,i);
imshow(im);title(int2str(i));
end
text(200,15,'请选择一幅照片(1-10)','fontsize',18);
testset=double(testset);
n= input('请选择一幅照片(1-10): ');
%选取训练集合中第n幅人脸照片作为输入进行分解
InputImage=testset(:,n)-averageface;
%计算输入照片在EigenFace基底上的投影坐标向量
coeff=InputImage'*eigenface;
%计算与训练集合中各图像的投影坐标向量的欧式距离
cdist=[];
for i=1:M
dist=(norm(coeff'-TrainCoeff(:,i)))^2;
cdist=[cdist dist];
end
[d index]=sort(cdist)
inum=[index(1),index(2),index(3)];
figure(6);
subplot(2,3,2);
InputImage=reshape(InputImage+averageface,imheight,imwidth);
imshow(uint8(InputImage));title('输入照片');
subplot(2,3,4);
im1=imset(:,index(1));
im1=reshape(im1,imheight,imwidth);
imshow(uint8(im1));title('匹配照片1');
subplot(2,3,5);
im2=imset(:,index(2));
im2=reshape(im2,imheight,imwidth);
imshow(uint8(im2));title('匹配照片2');
subplot(2,3,6);
im2=imset(:,index(3));
im2=reshape(im2,imheight,imwidth);
imshow(uint8(im2));title('匹配照片3');
pause;close;
%------------------------------------------------------
%step6:利用EigenFace基底对测试集合中的人脸照片进行
% 分解-合成(重构)
%注:如果用Eigenface基底直接对训练集合中的照片进行分解-重构
%则重构误差趋于0(1e-10)，如果对测试集合中的人脸(同一人，但表情姿态略有出入)
%进行重构，则可观察到重构误差，这是由于EigenFace基底虽正交但不完备%------------------------------------------------------
;%选取训练集合中第n幅人脸照片作为输入进行分解
InputImage=testset(:,n)-averageface;
%计算输入照片在EigenFace基底上的投影坐标
coeff=InputImage'*eigenface;
%在EigenFace基底上进行图像重构
ReconstructedImage=eigenface(:,1:M-1)*coeff';
%计算重构误差
diff=ReconstructedImage-InputImage;
norm(diff),
InputImage=reshape(InputImage+averageface,imheight,imwidth); ReconstructedImage=reshape(ReconstructedImage+averageface,imheight,imwidth); %figure(5);
%imagesc(ReconstructedImage); colormap('gray');
figure(7);
subplot(2,2,1);
im1=imset(:,n*2-1);
im1=reshape(im1,imheight,imwidth);
imshow(uint8(im1));title('训练集合中照片1');
subplot(2,2,2);
im2=imset(:,n*2);
im2=reshape(im2,imheight,imwidth);
imshow(uint8(im2));title('训练集合中照片2');
subplot(2,2,3);
InputImage=reshape(InputImage,imheight,imwidth);
imshow(uint8(InputImage));title('测试照片');
subplot(2,2,4);
ReconstructedImage=reshape(ReconstructedImage,imheight,imwidth);
imshow(uint8(ReconstructedImage));title('测试照片的重构');
pause;
return;
%----------------------------------------------------------------------------
%去除属于0的特征向量，同时将特征向量按特征值大小降序排列
%----------------------------------------------------------------------------
function [v d]= EliminateZero(v,d)
%取特征值
% dd=diag(d);
dd=abs(diag(d));
%找出等于0的特征值
ind=find(dd<1e-06);
% <1e-06);
%
% <1e-06);
% <1e-06);
%
%去除=0的特征值
dd(ind)=[];
%去除等于0的特征向量
v(:,ind)=[];
%特征值降序排序
[ddd ind]=sort(dd,'descend');
d=diag(ddd);
%重排特征向量
v=v(:,ind);
return;
/////////////////////////////////////////////////////Result Pictures///////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1.Example of reshape
tm=[]
ma=[1 9 7; 2 4 8;17 78 1; 22 4 79]
[imheight imwidth]=size(ma)
mac=reshape(ma,imheight*imwidth,1)
tm=[tm mac]
result:
mac =
1
2
17
22
9
4
78
4
7
8
1
79
2.Mean: compute mean value of matrix
mean(X,DIM) takes the mean along the dimension DIM of X.
Example: If X = [1 2 3; 3 3 6; 4 6 8; 4 7 7];
then mean(X,1) is [3.0000 4.5000 6.0000] and
mean(X,2) is [2.0000 4.0000 6.0000 6.0000].'
3.Imagesc
load clown
figure
subplot(1,2,1)
imshow(uint8(X))
subplot(1,2,2)
imagesc(X)
%colormap(gray)
axis image
title('Default CLim (= [1 81])')
4.Eig: Eigenvalues and eigenvectors.
E = eig(X) is a vector containing the eigenvalues of a square matrix X.
[V,D] = eig(X) produces a diagonal matrix D of eigenvalues and a
full matrix V whose columns are the corresponding eigenvectors so
that X*V = V*D.
B=[1 3 7; 31 8 4;11 89 3]
[VB,DB] = eig(B)
B*VB - VB*DB
[VN,DN] = eig(B,'nobalance')
B*VN - VN*DN
5.Norm
norm Matrix or vector norm.
norm(X,2) returns the 2-norm of X.
norm(X) is the same as norm(X,2).
图论算法及其MATLAB 程序代码 求赋权图G =(V ,E ,F )中任意两点间的最短路的Warshall-Floyd 算法： 设A =(a ij )n ×n 为赋权图G =(V ,E ,F )的矩阵,当v i v j ∈E 时a ij =F (v i v j ),否则取a ii =0,a ij =+∞(i ≠j ),d ij 表示从v i 到v j 点的距离,r ij 表示从v i 到v j 点的最短路中一个点的编号. ①赋初值.对所有i ,j ,d ij =a ij ,r ij =j .k =1.转向② ②更新d ij ,r ij .对所有i ,j ,若d ik +d k j ＜d ij ,则令d ij =d ik +d k j ,r ij =k ,转向③. ③终止判断.若d ii ＜0,则存在一条含有顶点v i 的负回路,终止;或者k =n 终止;否则令k =k +1,转向②. 最短路线可由r ij 得到. 例1求图6-4中任意两点间的最短路. 解：用Warshall-Floyd 算法,MATLAB 程序代码如下： n=8;A=[0281Inf Inf Inf Inf 206Inf 1Inf Inf Inf 8607512Inf 1Inf 70Inf Inf 9Inf Inf 15Inf 03Inf 8 Inf Inf 1Inf 3046 Inf Inf 29Inf 403 Inf Inf Inf Inf 8630];%MATLAB 中,Inf 表示∞ D=A;%赋初值 for (i=1:n)for (j=1:n)R(i,j)=j;end ;end %赋路径初值 for (k=1:n)for (i=1:n)for (j=1:n)if (D(i,k)+D(k,j) Floyd算法Matlab程序第一种: %floyd.m %采用floyd算法计算图a中每对顶点最短路 %d是矩离矩阵 %r是路由矩阵 function ,d,r,=floyd(a) n=size(a,1); d=a; for i=1:n for j=1:n r(i,j)=j; end end r for k=1:n for i=1:n for j=1:n if d(i,k)+d(k,j) end k d r end 第二种: %Floyd算法 %解决最短路径问题，是用来调用的函数头文件 %[D,path]=floyd(a) %输入参数a是求图的带权邻接矩阵，D(i,j)表示i到j的最短距 离,path(i,j)i,j之间最短路径上顶点i的后继点 %[D,path,min1,path1]=floyd(a,i,j) %输入参数a是所求图的带权邻接矩阵，i,j起点终点,min1表示i与j最短距离，path1为最短路径function [D,path,min1,path1]=floyd(a,start,terminal) D=a;n=size(D,1);path=zeros(n,n); for i=1:n for j=1:n if D(i,j)~=inf path(i,j)=j; end end end for k=1:n for i=1:n for j=1:n if D(i,k)+D(k,j) 实验四：Floyd 算法 一、实验目的 利用MATLAB 实现Floyd 算法，可对输入的邻接距离矩阵计算图中任意两点间的最短距离矩阵和路由矩阵，且能查询任意两点间的最短距离和路由。 二、实验原理 Floyd 算法适用于求解网络中的任意两点间的最短路径：通过图的权值矩阵求出任意两点间的最短距离矩阵和路由矩阵。优点是容易理解，可以算出任意两个节点之间最短距离的算法，且程序容易实现，缺点是复杂度达到，不适合计算大量数据。 Floyd 算法可描述如下： 给定图G 及其边(i , j )的权w i, j (1≤i≤n ,1≤j≤n) F0：初始化距离矩阵W(0)和路由矩阵R(0)。其中： F1：已求得W(k-1)和R(k-1)，依据下面的迭代求W(k)和R(k) F2：若k≤n，重复F1；若k>n，终止。 三、实验内容 1、用MATLAB 仿真工具实现Floyd 算法：给定图G 及其边(i , j )的权 w i , j (1≤i≤n ,1≤j≤n) ，求出其各个端点之间的最小距离以及路由。 （1）尽可能用M 函数分别实现算法的关键部分，用M 脚本来进行算法结 果验证； （2）分别用以下两个初始距离矩阵表示的图进行算法验证： 分别求出W(7)和R(7)。 2、根据最短路由矩阵查询任意两点间的最短距离和路由 （1）最短距离可以从最短距离矩阵的ω(i,j)中直接得出； （2）相应的路由则可以通过在路由矩阵中查找得出。由于该程序中使用的是前向矩阵，因此在查找的过程中，路由矩阵中r(i,j)对应的值为Vi 到Vj 路由上的下一个端点，这样再代入r(r(i,j),j)，可得到下下个端点，由此不断循环下去，即可找到最终的路由。 （3）对图1，分别以端点对V4 和V6, V3 和V4 为例，求其最短距离和路由；对图2，分别以端点对V1 和V7，V3 和V5，V1 和V6 为例，求其最短距离和路由。 3、输入一邻接权值矩阵，求解最短距离和路由矩阵，及某些点间的最短路径。 四、采用的语言 MatLab 源代码： 【func1.m】 function [w r] = func1(w) n=length(w); x = w; r = zeros(n,1);%路由矩阵的初始化 for i=1:1:n for j=1:1:n if x(i,j)==inf r(i,j)=0; else r(i,j)=j; end, end end; %迭代求出k次w值 for k=1:n a=w; s = w; for i=1:n 精心整理 图论算法matlab实现 求最小费用最大流算法的 MATLAB 程序代码如下： n=5;C=[0 15 16 0 0 0 0 0 13 14 for while for for(i=1:n)for(j=1:n)if(C(i,j)>0&f(i,j)==0)a(i,j)=b(i,j); elseif(C(i,j)>0&f(i,j)==C(i,j))a(j,i)=-b(i,j); elseif(C(i,j)>0)a(i,j)=b(i,j);a(j,i)=-b(i,j);end;end;end for(i=2:n)p(i)=Inf;s(i)=i;end %用Ford 算法求最短路, 赋初值 for(k=1:n)pd=1; %求有向赋权图中vs 到vt 的最短路 for(i=2:n)for(j=1:n)if(p(i)>p(j)+a(j,i))p(i)=p(j)+a(j,i);s(i)=j;pd=0;end ;end;end if(pd)break;end;end %求最短路的Ford 算法结束 if(p(n)==Inf)break;end %不存在vs 到vt 的最短路, 算法终止. 注意在求最小费用最大流时构造有 while if elseif if if pd=0; 值 t=n; if elseif if(s(t)==1)break;end %当t 的标号为vs 时, 终止调整过程 t=s(t);end if(pd)break;end%如果最大流量达到预定的流量值 wf=0; for(j=1:n)wf=wf+f(1,j);end;end %计算最大流量 zwf=0;for(i=1:n)for(j=1:n)zwf=zwf+b(i,j)*f(i,j);end;end%计算最小费用 一、实验目的 利用MATLAB实现Floyd算法，可对输入的邻接距离矩阵计算图中任意两点间的最短距离矩阵和路由矩阵，且能查询任意两点间的最短距离和路由。 二、实验原理 Floyd 算法适用于求解网络中的任意两点间的最短路径：通过图的权值矩阵求出任意两点间的最短距离矩阵和路由矩阵。优点是容易理解，可以算出任意两个 节点之间最短距离的算法，且程序容易实现，缺点是复杂度达到，不适合计算大量数据。 Floyd 算法可描述如下： 给定图G及其边(i , j ) 的权w, j (1 < i < n ,1 图论算法及其MATLAB程序代码 求赋权图G = (V, E , F )中任意两点间的最短路的Warshall-Floyd算法： 设A = (a ij )n×n为赋权图G = (V, E , F )的矩阵, 当v i v j∈E时a ij= F (v i v j), 否则取a ii=0, a ij = +∞(i≠j ), d ij表示从v i到v j点的距离, r ij表示从v i到v j点的最短路中一个点的编号. ①赋初值. 对所有i, j, d ij = a ij, r ij = j. k = 1. 转向② ②更新d ij, r ij . 对所有i, j, 若d ik + d k j＜d ij, 则令d ij = d ik + d k j, r ij = k, 转向③. ③终止判断. 若d ii＜0, 则存在一条含有顶点v i的负回路, 终止; 或者k = n终止; 否则令k = k + 1, 转向②. 最短路线可由r ij得到. 例1求图6-4中任意两点间的最短路. 图6-4 解：用Warshall-Floyd算法, MA TLAB程序代码如下： n=8;A=[0 2 8 1 Inf Inf Inf Inf 2 0 6 Inf 1 Inf Inf Inf 8 6 0 7 5 1 2 Inf 1 Inf 7 0 Inf Inf 9 Inf Inf 1 5 Inf 0 3 Inf 8 Inf Inf 1 Inf 3 0 4 6 Inf Inf 2 9 Inf 4 0 3 Inf Inf Inf Inf 8 6 3 0]; % MATLAB中, Inf表示∞ D=A; %赋初值 for(i=1:n)for(j=1:n)R(i,j)=j;end;end%赋路径初值 for(k=1:n)for(i=1:n)for(j=1:n)if(D(i,k)+D(k,j) Dijkstra 算法Matlab 实现。 %求一个点到其他各点的最短路径 function [min,path]=dijkstra(w,start,terminal) %W 是邻接矩阵 %start 是起始点 %terminal 是终止点 %min 是最短路径长度 %path 是最短路径 n=size(w,1); label(start)=0; f(start)=start; for i=1:n if i~=start label(i)=inf; end end s(1)=start; u=start; while length(s) 基于matlab的floyd算法 matlab计算最短路径 function [d,path]=floyd(a,sp,ep) % floyd - 最短路问题 % % Syntax: [d,path]=floyd(a,sp,ep) % % Inputs: % a - 距离矩阵是指i到j之间的距离，可以是有向的 % sp - 起点的标号 % ep - 终点的标号 % % Outputs: % d - 最短路的距离 % path - 最短路的路径 % a =[ 0 50 inf; 50 0 15 ; Inf 15 0 ];% a(i,j)，从节点i到j之间的距离 % [d,path]=floyd(a,2,5) sp=3; ep=1; n=size(a,1); D=a; path=zeros(n,n); for i=1:n for j=1:n if D(i,j)~=inf path(i,j)=j; %j是i的后续点 end end end for k=1:n for i=1:n for j=1:n if D(i,j)>D(i,k)+D(k,j) D(i,j)=D(i,k)+D(k,j); path(i,j)=path(i,k); end end end end p=[sp]; mp=sp; for k=1:n if mp~=ep d=path(mp,ep); p=[p,d]; mp=d; end end d=D(sp,ep) path=p 试计算下图的最短路径， 1.起点C点，终点A点。 2.起点A点，终点G点。 3.起点D点，终点F点。 试计算下图的最短路径， 1.起点F点，终点A点。 2. 起点E点，终点C点。 超全的图论程序 关注微信公众号“超级数学建模”，教你做有料、有趣的数模人 程序一：可达矩阵算法 function P=dgraf(A) n=size(A,1); P=A; for i=2:n P=P+A^i; end P(P~=0)=1; P; 程序二：关联矩阵和邻接矩阵互换算法 function W=incandadf(F,f) if f==0 m=sum(sum(F))/2; n=size(F,1); W=zeros(n,m); k=1; for i=1:n for j=i:n if F(i,j)~=0 W(i,k)=1; W(j,k)=1; k=k+1; end end end elseif f==1 m=size(F,2); n=size(F,1); W=zeros(n,n); for i=1:m a=find(F(:,i)~=0); W(a(1),a(2))=1; W(a(2),a(1))=1; end else fprint('Please imput the right value of f'); end W; 程序三：有向图关联矩阵和邻接矩阵互换算法 function W=mattransf(F,f) if f==0 m=sum(sum(F)); n=size(F,1); W=zeros(n,m); k=1; for i=1:n for j=i:n if F(i,j)~=0 W(i,k)=1; W(j,k)=-1; k=k+1; end end end elseif f==1 m=size(F,2); n=size(F,1); W=zeros(n,n); for i=1:m a=find(F(:,i)~=0); if F(a(1),i)==1 W(a(1),a(2))=1; else W(a(2),a(1))=1; end end else fprint('Please imput the right value of f'); end W; 中国数学建模-数学工具-Floyd最短路算法的MATLAB程序 wh-ee 重登录隐身用户控制面板搜索风格论坛状态论坛展区社区服务社区休闲网站首页退出 >> Matlab,Mathematica,maple,几何画板,word,sas,spss...使用方法技巧我的收件箱(0) 中国数学建模→数学建模→数学工具→Floyd最短路算法的MATLAB程序 您是本帖的第90 个阅读者 * 贴子主题：Floyd最短路算法的MATLAB程序 hanlong 等级：新手上路 文章：28 积分：125 门派：☆nudter☆ 注册：2004-5-20 鲜花(1) 鸡蛋(0) 楼主 Floyd最短路算法的MATLAB程序 %floyd.m %采用floyd算法计算图a中每对顶点最短路 %d是矩离矩阵 %r是路由矩阵 function [d,r]=floyd(a) n=size(a,1); d=a; for i=1:n for j=1:n r(i,j)=j; end end r for k=1:n for i=1:n for j=1:n if d(i,k)+d(k,j) r(i,j)=r(i,k) end end end k d r end 2004-5-24 1:04:35 wanggaoyang 等级：新手上路 文章：9 积分：106 门派：☆nudter☆ 注册：2004-5-24 第2 楼 顶 2004-5-28 23:06:16 feifei7 头衔：蓝魂行者 等级：新手上路 文章：36 积分：258 门派：桃花岛 注册：2004-4-27 第3 楼 ^ #include for (i=0;i 层次分析法： %层次分析法的matlab程序 disp('请输入判断矩阵A(n阶)'); A=input('A=');%判断矩阵 [n,n]=size(A);%n阶 x=ones(n,100); y=ones(n,100); m=zeros(1,100); m(1)=max(x(:,1)); y(:,1)=x(:,1); x(:,2)=A*y(:,1); m(2)=max(x(:,2)); y(:,2)=x(:,2)/m(2); p=0.0001;i=2;k=abs(m(2)-m(1)); while k>p i=i+1; x(:,i)=A*y(:,i-1); m(i)=max(x(:,i)); y(:,i)=x(:,i)/m(i); k=abs(m(i)-m(i-1)); end a=sum(y(:,i)); w=y(:,i)/a; t=m(i); disp('权向量');disp(w); disp('最大特征值');disp(t); %以下是一致性检验 CI=(t-n)/(n-1);RI=[0 0 0.52 0.89 1.12 1.26 1.36 1.41 1.46 1.49 1.52 1.54 1.56 1.58 1.59];%RI根据n 来选 CR=CI/RI(n);%判断是否通过一致性检验 if CR<0.10 disp('此矩阵的一致性可以接受!'); disp('CI=');disp(CI); disp('CR=');disp(CR); else disp('此矩阵的一致性不可以接受!'); end 主成分分析： function [lambda,T,fai]=MSA2(A)%lambda特征根,T特征向量,fai贡献率%求标准化后的协差矩阵,再求特征根和特征向量 %标准化处理 [p,n]=size(A);%p行n列,列为几种 for j=1:n mju(j)=mean(A(:,j)); sigma(j)=sqrt(cov(A(:,j))); end for i=1:p for j=1:n Y(i,j)=(A(i,j)-mju(j))/sigma(j); end end sigmaY=cov(Y); %求X标准化的协差矩阵的特征根和特征向量 [T,lambda]=eig(sigmaY); disp('特征根(由小到大):'); disp(lambda); disp('特征向量:'); disp(T); %方差贡献率; Xsum=sum(sum(lambda,2),1); for i=1:n fai(i)=lambda(i,i)/Xsum; end disp('方差贡献率:'); disp(fai); u=T(:,n); B=[]; h=length(A(:,1)); for k=1:n m1=mean(A(:,k)); t=(A(:,k)-m1).^2; m2=sqrt(sum(t))/(h-1); B=[B,(A(:,k)-m1)./m2]; end y=B*u; x1=1:1:length(y); plot(x1,y); xlabel('时间/小时') ylabel('综合指标') Floyd算法 算法能实现求一个图中任意两点间最断距离，并给出路径。 设图G的顶点数为n，D=(d ij)n×n，为其距离矩阵，Floyd算法步骤如下： Step 1：输入距离矩阵D； Step 2：k=1， Step 3：i=1； Step 4：d ij=min(d ij,d ik+d kj),j=1,2,…,n; Step 5：i++；如果i≤n，转Step 4； Step 6：k++；如果k≤n，转Step 3；否则转Step 7； Step 7：程序结束，输出结果。 用Floyd算法求图中各顶点之间的最短路的Matlab程序。 M=99999; a=[0,20,14,M,M,M,M M,0,M,15,12,M,M M,M,0,10,M,13,M M,M,M,0,8,M,9 M,M,M,M,0,8,10 M,M,M,M,M,0,12 M,M,M,M,M,M,0]; length=floyd(a) r代表任意城市之间经过的路径,这个结果比起Mathematica来说不好理解，这个也是我自己对比上面结果的理解。 第一行代表1经过1，2，3，4，5，6，7的路径； 第一个数是1，代表1直接到1最短 （1，1） 第二个数是2，代表1直接到2最短 （1，2） 第三个数是3，代表1直接到3最短 （1，3） 第四个数是3，代表1经过3再到4最短 （1，3，4） 第五个数是2，代表1经过2再到5最短 （1，2，5） 第六个数是3，代表1经过3再到6最短 （1，3，6） 第七个数是3，代表1经过3再到7最短 （1，3，7）？？和前面不一样 出现问题了，其实错误在于3不是直接到7最短造成的（我们取的3到7距离为99999） 而前面3个红色的路径是不会出现这种情况的（可以看下面的表）。 可以看第三行 3是经过4再到7最短！ 所以（1，3，4，7） 任意城市之间经过的路径我们只需要求红色括号里面的路径就可以了（用上面的分析可以得到答案），其他的已经是最短路了。 { {1,1} {1,2} {1,3} {} {} {} {} {2,1} {2,2} {2,3} {2,4} {2,5} {2，5，6} {2，5，7} {3,1} {3,2} {3,3} {3,4} {3，4，5} {3,6} {3，4，7} {4,1} {4,2} {4,3} {4,4} {4,5} {4，5，6} {4,7} {5,1} {5,2} {5,3} {5,4} {5,5} {5,6} {5,7} {6,1} {6,2} {6,3} {6,4} {6,5} {6,6} {6,7} {7,1} {7,2} {7,3} {7,4} {7,5} {7,6} {7,7} } 如果熟练后也可以很快得到结果！ Floyd最短路算法的MATLAB程序 2006-08-17 20:14 %floyd.m %采用floyd算法计算图a中每对顶点最短路 %d是矩离矩阵 %r是路由矩阵 function [d,r]=floyd(a) n=size(a,1); d=a; for i=1:n for j=1:n r(i,j)=j; end end r for k=1:n for i=1:n for j=1:n if d(i,k)+d(k,j) { double min=INF; int u=v0; for(int j=0;j function [DeD,aver_DeD]=Degree_Distribution(A) %% 求网络图中各节点的度及度的分布曲线 %% 求解算法：求解每个节点的度，再按发生频率即为概率，求P(k) %A————————网络图的邻接矩阵 %DeD————————网络图各节点的度分布 %aver_DeD———————网络图的平均度 N=size(A,2); DeD=zeros(1,N); for i=1:N % DeD(i)=length(find((A(i,:)==1))); DeD(i)=sum(A(i,:)); end aver_DeD=mean(DeD); if sum(DeD)==0 disp('该网络图只是由一些孤立点组成'); return; else figure; bar([1:N],DeD); xlabel('节点编号n'); ylabel('各节点的度数K'); title('网络图中各节点的度的大小分布图'); end figure; M=max(DeD); for i=1:M+1; %网络图中节点的度数最大为M,但要同时考虑到度为0的节点的存在性 N_DeD(i)=length(find(DeD==i-1)); % DeD=[2 2 2 2 2 2] end P_DeD=zeros(1,M+1); P_DeD(:)=N_DeD(:)./sum(N_DeD); bar([0:M],P_DeD,'r'); xlabel('节点的度 K'); ylabel('节点度为K的概率 P(K)'); title('网络图中节点度的概率分布图'); function [C,aver_C]=Clustering_Coefficient(A) %% 求网络图中各节点的聚类系数及整个网络的聚类系数 %% 求解算法：求解每个节点的聚类系数，找某节点的所有邻居，这些邻居节点构成一个子图 %% 从A中抽出该子图的邻接矩阵，计算子图的边数，再根据聚类系数的定义，即可算出该节点的聚类系数 %A————————网络图的邻接矩阵 %C————————网络图各节点的聚类系数 %aver———————整个网络图的聚类系数 N=size(A,2); C=zeros(1,N); for i=1:N aa=find(A(i,:)==1); %寻找子图的邻居节点 if isempty(aa) disp(['节点',int2str(i),'为孤立节点，其聚类系数赋值为0']); C(i)=0; else m=length(aa); if m==1 disp(['节点',int2str(i),'只有一个邻居节点，其聚类系数赋值为0']); C(i)=0; else B=A(aa,aa) % 抽取子图的邻接矩阵 C(i)=length(find(B==1))/(m*(m-1)); end end end aver_C=mean(C) 每对顶点之间的最短路径 计算赋权图中各对顶点之间最短路径，显然可以调用Dijkstra 算法。具体方法是：每次以不同的顶点作为起点，用Dijkstra 算法求出从该起点到其余顶点的最短路径，反复执行n 次这样的操作，就可得到从每一个顶点到其它顶点的最短路径。这种算法的时间复杂度为)(3n O 。第二种解决这一问题的方法是由Floyd R W 提出的算法，称之为Floyd 算法。 假设图G 权的邻接矩阵为0A ， ????? ???????=nn n n n n a a a a a a a a a A 2122221 112110 来存放各边长度，其中： 0=ii a n i ,,2,1 =； ∞=ij a j i ,之间没有边， 在程序中以各边都不可能达到的充分大的数代替； ij ij w a = ij w 是j i ,之间边的长度，n j i ,,2,1, =。 对于无向图，0A 是对称矩阵，ji ij a a =。 Floyd 算法的基本思想是：递推产生一个矩阵序列n k A A A A ,,,,,10 ，其中),(j i A k 表示从顶点i v 到顶点j v 的路径上所 经过的顶点序号不大于k 的最短路径长度。 计算时用迭代公式： )),(),(),,(m in(),(111j k A k i A j i A j i A k k k k ---+= k 是迭代次数，n k j i ,,2,1,, =。 最后，当n k =时，n A 即是各顶点之间的最短通路值。 例10某公司在六个城市621,,,c c c 中有分公司，从i c 到j c 的直接航程票价记在下述矩阵的),(j i 位置上。（∞表示无直接航路），请帮助该公司设计一张城市1c 到其它城市间的票价最便 宜的路线图。 ???????? ??????????∞∞∞∞∞∞05525251055010202525100102040 2010015252015050102540500 矩阵path 用来存放每对顶点之间最短路径上所经过的顶点的序号。Floyd 算法的Matlab 程序如下： clear; clc; M=10000; a(1,:)=[0,50,M,40,25,10]; a(2,:)=[zeros(1,2),15,20,M,25]; a(3,:)=[zeros(1,3),10,20,M]; a(4,:)=[zeros(1,4),10,25]; a(5,:)=[zeros(1,5),55]; a(6,:)=zeros(1,6); b=a+a';path=zeros(length(b)); for k=1:6 for i=1:6Floyd算法Matlab程序
Floyd算法_计算最短距离矩阵和路由矩阵_查询最短距离和路由_matlab实验报告
matlab图论程序算法大全
Floyd算法_计算最短距离矩阵和路由矩阵_查询最短距离和路由_matlab实验报告
图论算法及其MATLAB程序代码
Dijkstra、Floyd算法Matlab_Lingo实现
基于matlab的floyd算法+matlab计算最短路径
超全图论matlab程序
Floyd最短路算法
floyd算法C++实现
一些简单的算法MATLAB代码
floyd算法 matlab程序
Floyd最短路算法的MATLAB程序
复杂网络主要拓扑参数的matlab实现
Floyd算法每对顶点之间的最短路径
相关主题
文本预览