利用梯度下降法实现线性回归的算法及matlab实现
1.线性回归算法概述
线性回归属于监督学习，因此方法和监督学习应该是一样的，先给定一个训练集，根据这个训练集学习出一个线性函数，然后测试这个函数训练的好不好（即此函数是否足够拟合训练集数据），挑选出最好的函数（cost function最小）即可；
注意：
(1)因为是线性回归，所以学习到的函数为线性函数，即直线函数；
(2)线性回归可分为单变量线性回归和多变量线性回归；对于单变量线性回归而言，只有一个输入变量x；
(1).单变量线性回归
我们能够给出单变量线性回归的模型：
我们常称x为feature，h(x)为hypothesis；上述模型中的θ0和θ1在代码中分别用theta0和theta1表示。
从上面“方法”中，我们肯定有一个疑问，怎么样能够看出线性函数拟合的好不好呢？我们需要使用到Cost Function（代价函数），代价函数越小，说明线性回归地越好（和训练集拟合地越好），当然最小就是0，即完全拟合。cost Function的内部构造如下面公式所述：
其中：
表示向量x中的第i个元素；
表示向量y中的第i个元素；
表示已知的假设函数；
m为训练集的数量；
虽然给定一个函数，我们能够根据cost function知道这个函数拟合的好不好，但是毕竟函数有这么多，总不可能一个一个试吧？因此我们引出了梯度下降：能够找出cost function函数的最小值；
梯度下降原理：将函数比作一座山，我们站在某个山坡上，往四周看，从哪个方向向下走一小步，能够下降的最快；当然解决问题的方法有很多，梯度下降只是其中一个，还有一种方法叫Normal Equation；
方法：
(1)先确定向下一步的步伐大小，我们称为Learning rate (alpha)；
(2)任意给定一个初始值：(用theta0和theta1表示)；
(3)确定一个向下的方向，并向下走预先规定的步伐，并更新；
(4)当下降的高度小于某个定义的值，则停止下降；
：
算法
特点：
(1)初始点不同，获得的最小值也不同，因此梯度下降求得的只是局部最小值；
(2)越接近最小值时，下降速度越慢；
梯度下降能够求出一个函数的最小值；
线性回归需要使得cost function的最小；
因此我们能够对cost function运用梯度下降，即将梯度下降和线性回归进行整合，如下图所示：
上式中右边的公式推导过程如下：
?
?
??
?
??
?
从上面的推导中可以看出，要想满足梯度下降的条件，则?项后面必须乘以对应的输入信号。
梯度下降是通过不停的迭代，而我们比较关注迭代的次数，因为这关系到梯度下降的执行速度，为了减少迭代次数，因此引入了Feature Scaling。
(2).Feature Scaling
此种方法应用于梯度下降，为了加快梯度下降的执行速度；
思想：将各个feature的值标准化，使得取值范围大致都在-1<=x<=1之间；
(3).常用的方法是Mean Normalization(均值归一化处理)
或者：
[X-mean(X)]/std(X);
(4).收获汇总
学习速率的大小对于系统是否收敛有决定性的影响。如果学习速率太大，那么可能导致系统震荡发撒；如果学习速率太小，那么可能导致系统收敛速度变慢。
(a)根据给定数据架设预测函数h(x)
(b)计算代价函数J
(c)计算各参数偏导
(d)更新参数
(e)重复2~4直到代价函数跟新的步长小于设定值或者是重复次数达到预设值。
为了在相同学习速率的前提下加快相同收敛，可采用训练数据归一化的方法来对样本数据进行预处理。采用均值均值归一化处理后，缩小了原来数据的变化幅度，从而可极大地提高学习速率，从而提高了梯度下降的执行速度。在第四部分的matlab代码中，对输入值向量进行归一化处理之后，将学习速率从0.01提高到1.9,从而将获得系统收敛的epoch次数由10000次减少到300次。
2.Matlab实现单变量梯度下降的线性回归
(1).Gradient descend code 1
clear all
clc
% training sample data;
p0=3;
p1=7;
x=1:3;
y=p0+p1*x;
num_sample=size(y,2);
% gradient descending process
% initial values of parameters
theta0=1;
theta1=3;
%learning rate
alpha=0.02;
% if alpha is too large, the final error will be much large.
% if alpha is too small, the convergence will be slow
epoch=500;
for k=1:epoch
v_k=k
h_theta_x=theta0+theta1*x; % hypothesis function
Jcost(k)=((h_theta_x(1)-y(1))^2+(h_theta_x(2)-y(2))^2+(h_theta_x(3)-y(3))^2)/num_sample
theta0=theta0-alpha*((h_theta_x(1)-y(1))+(h_theta_x(2)-y(2))+(h_theta_x(3)-y(3)))/num_sample;
theta1=theta1-alpha*((h_theta_x(1)-y(1))*x(1)+(h_theta_x(2)-y(2))*x(2)+(h_theta_x(3)-y(3))*x(3))/num_sample; end
plot(Jcost)
(2).Gradient descend code 2
clear all
clc
% training sample data;
p0=26;
p1=73;
x=1:3;
y=p0+p1*x;
num_sample=size(y,2);
% gradient descending process
% initial values of parameters
theta0=1;
theta1=3;
%learning rate
alpha=0.08;
% if alpha is too large, the final error will be much large.
% if alpha is too small, the convergence will be slow
epoch=500;
for k=1:epoch
v_k=k
h_theta_x=theta0+theta1*x; % hypothesis function
Jcost(k)=((h_theta_x(1)-y(1))^2+(h_theta_x(2)-y(2))^2+(h_theta_x(3)-y(3))^2)/num_sample;
theta0=theta0-alpha*((h_theta_x(1)-y(1))+(h_theta_x(2)-y(2))+(h_theta_x(3)-y(3)))/num_sample;
theta1=theta1-alpha*((h_theta_x(1)-y(1))*x(1)+(h_theta_x(2)-y(2))*x(2)+(h_theta_x(3)-y(3))*x(3))/num_sample;
% disp('*********comp 1**************');
r1=((h_theta_x(1)-y(1))+(h_theta_x(2)-y(2))+(h_theta_x(3)-y(3)));
r2=sum(h_theta_x-y);
% disp('*********comp 2**************');
r3=((h_theta_x(1)-y(1))^2+(h_theta_x(2)-y(2))^2+(h_theta_x(3)-y(3))^2);
r4=sum((h_theta_x-y).^2);
% disp('*********comp 3**************');
r5=((h_theta_x(1)-y(1))*x(1)+(h_theta_x(2)-y(2))*x(2)+(h_theta_x(3)-y(3))*x(3));
r6=sum((h_theta_x-y).*x);
if((r1~=r2)||(r3~=r4)||(r5~=r6))
disp('***wrong result******')
end
end
plot(Jcost)
3.线性回归与单神经元的对应关系
单变量线性回归的模型：
与单神经元模型具有对应的关系，可实现相同的功能。其中θ0相当于单神经元的偏置值bias值，x相当于单神经元的单个输入，θ
相当于单神经元的权值w。如果?不包含θ0，那么则可能无法准确完成线性回归的功能。因此，对于单神经元结构而言，偏置值1
bias是必不可少的。
x
y
图1. 单神经元单变量输入与线性回归的对应模型
对于双变量线性回归模型： ?
而言，可实现双输入单神经元的功能。其偏置值对应θ0。
x1
y
x2
图2. 单神经元双变量输入与线性回归的对应模型
4.Matlab实现多变量梯度下降的线性回归
(1).未进行归一化处理之前的matlab代码
两个输入变量x1和x2.
clear all
clc
% training sample data;
p0=6;
p1=7;
p2=2;
% the first group of inputs
x1=[7 9 12 5 4];
x2=[1 8 21 3 5];
% the second group of inputs
x1=[7 9 12 5 4 3];
x2=[1 8 21 3 5 26];
y=p0+p1*x1+p2*x2;
num_sample=size(y,2);
% gradient descending process
% initial values of parameters
theta0=9;
theta1=3;
theta2=9;
% theta0=19;theta1=23;theta2=91;
% theta0=0;theta1=0;theta2=0;
%learning rate
alpha=0.01; % good for the system with the first group of inputs
alpha=0.005; % good for the system the second group of inputs
% alpha=0.02; % bad for the system. the system will be unstable
% if alpha is too large, the final error will be much large.
% if alpha is too small, the convergence will be slow
epoch=10000;
for k=1:epoch
v_k=k
h_theta_x=theta0+theta1*x1+theta2*x2; % hypothesis function
Jcost(k)=((h_theta_x(1)-y(1))^2+(h_theta_x(2)-y(2))^2+(h_theta_x(3)-y(3))^2+(h_theta_x(4)-y(4))^2+(h_theta_x(5)-y(5) )^2)/num_sample;
theta0=theta0-alpha*((h_theta_x(1)-y(1)) +(h_theta_x(2)-y(2)) +(h_theta_x(3)-y(3))
+(h_theta_x(4)-y(4)) +(h_theta_x(5)-y(5)))/num_sample;
theta1=theta1-alpha*((h_theta_x(1)-y(1))*x1(1)+(h_theta_x(2)-y(2))*x1(2)+(h_theta_x(3)-y(3))*x1(3)+(h_theta_x(4)-y(4 ))*x1(4)+(h_theta_x(5)-y(5))*x1(5))/num_sample;
theta2=theta2-alpha*((h_theta_x(1)-y(1))*x2(1)+(h_theta_x(2)-y(2))*x2(2)+(h_theta_x(3)-y(3))*x2(3)+(h_theta_x(4)-y(4 ))*x2(4)+(h_theta_x(5)-y(5))*x2(5))/num_sample;
% % disp('*********comp 1**************');
% r1=((h_theta_x(1)-y(1))+(h_theta_x(2)-y(2))+(h_theta_x(3)-y(3)));
% r2=sum(h_theta_x-y);
% % disp('*********comp 2**************');
% r3=((h_theta_x(1)-y(1))^2+(h_theta_x(2)-y(2))^2+(h_theta_x(3)-y(3))^2);
% r4=sum((h_theta_x-y).^2);
% % disp('*********comp 3**************');
% r5=((h_theta_x(1)-y(1))*x(1)+(h_theta_x(2)-y(2))*x(2)+(h_theta_x(3)-y(3))*x(3));
% r6=sum((h_theta_x-y).*x);
% if((r1~=r2)||(r3~=r4)||(r5~=r6))
% disp('***wrong result******')
% end
end
yt=theta0+theta1*x1+theta2*x2
plot(Jcost)
上述matlab代码中，当学习速率alpha小于等于0.01时，系统能收敛；但是，如果学习速率alpha越小，系统收敛的速度越慢。当学习速率alpha大于等于0.02时，系统发散，无法收敛。为了使得系统尽快收敛下一步采用归一化的方法来加速系统收敛。
(2).对表达式进行简化后的matlab代码
clear all
clc
% training sample data;
p0=6;
p1=7;
p2=2;
x1=[7 9 12 5 4];
x2=[1 8 21 3 5];
y=p0+p1*x1+p2*x2;
num_sample=size(y,2);
% gradient descending process
% initial values of parameters
theta0=9;
theta1=3;
theta2=9;
% theta0=19;theta1=23;theta2=91;
% theta0=0;theta1=0;theta2=0;
%learning rate
alpha=0.005; % good for the system
alpha=0.01; % good for the system
% alpha=0.02; % bad for the system. the system will be unstable
% if alpha is too large, the final error will be much large.
% if alpha is too small, the convergence will be slow
epoch=10000;
for k=1:epoch
v_k=k
h_theta_x=theta0+theta1*x1+theta2*x2; % hypothesis function
Jcost(k)=sum((h_theta_x-y).^2)/num_sample;
r0=sum(h_theta_x-y);
theta0=theta0-alpha*r0/num_sample;
r1=sum((h_theta_x-y).*x1);
theta1=theta1-alpha*r1/num_sample;
r2=sum((h_theta_x-y).*x2);
theta2=theta2-alpha*r2/num_sample;
%
Jcost(k)=((h_theta_x(1)-y(1))^2+(h_theta_x(2)-y(2))^2+(h_theta_x(3)-y(3))^2+(h_theta_x(4)-y(4))^2+(h_theta_x(5)-y(5)
)^2)/num_sample;
%
theta1=theta1-alpha*((h_theta_x(1)-y(1))*x1(1)+(h_theta_x(2)-y(2))*x1(2)+(h_theta_x(3)-y(3))*x1(3)+(h_theta_x(4)-y(4 ))*x1(4)+(h_theta_x(5)-y(5))*x1(5))/num_sample;
%
theta2=theta2-alpha*((h_theta_x(1)-y(1))*x2(1)+(h_theta_x(2)-y(2))*x2(2)+(h_theta_x(3)-y(3))*x2(3)+(h_theta_x(4)-y(4 ))*x2(4)+(h_theta_x(5)-y(5))*x2(5))/num_sample;
% % disp('*********comp 1**************');
% r1=((h_theta_x(1)-y(1))+(h_theta_x(2)-y(2))+(h_theta_x(3)-y(3)));
% r2=sum(h_theta_x-y);
% % disp('*********comp 2**************');
% r3=((h_theta_x(1)-y(1))^2+(h_theta_x(2)-y(2))^2+(h_theta_x(3)-y(3))^2);
% r4=sum((h_theta_x-y).^2);
% % disp('*********comp 3**************');
% r5=((h_theta_x(1)-y(1))*x(1)+(h_theta_x(2)-y(2))*x(2)+(h_theta_x(3)-y(3))*x(3));
% r6=sum((h_theta_x-y).*x);
% if((r1~=r2)||(r3~=r4)||(r5~=r6))
% disp('***wrong result******')
% end
end
yt=theta0+theta1*x1+theta2*x2
plot(Jcost)
(3).进行归一化处理之后的matlab代码
clear all
clc
% training sample data;
p0=6;
p1=7;
p2=2;
x1=[7 9 12 5 4];
x2=[1 8 21 3 5];
x1_mean=mean(x1)
x1_max=max(x1)
x1_min=min(x1)
x2_mean=mean(x2)
x2_max=max(x2)
x2_min=min(x2)
x1=(x1-x1_mean)/(x1_max-x1_min) x2=(x2-x2_mean)/(x2_max-x2_min)
% x1=[7 9 12 5 4 3];
% x2=[1 8 21 3 5 26];
y=p0+p1*x1+p2*x2;
num_sample=size(y,2);
% gradient descending process
% initial values of parameters theta0=9;
theta1=3;
theta2=9;
% theta0=19;theta1=23;theta2=91;
% theta0=0;theta1=0;theta2=0;
%learning rate
alpha=1.9; % good for the system
% if alpha is too large, the final error will be much large.
% if alpha is too small, the convergence will be slow
epoch=300;
for k=1:epoch
v_k=k
h_theta_x=theta0+theta1*x1+theta2*x2; % hypothesis function Jcost(k)=sum((h_theta_x-y).^2)/num_sample;
r0=sum(h_theta_x-y);
theta0=theta0-alpha*r0/num_sample;
r1=sum((h_theta_x-y).*x1);
theta1=theta1-alpha*r1/num_sample;
r2=sum((h_theta_x-y).*x2);
theta2=theta2-alpha*r2/num_sample;
end
yt=theta0+theta1*x1+theta2*x2
plot(Jcost)
(4).进行归一化处理之后的matlab代码2
三个个输入变量x1、x2和x3.
clear all
clc
% training sample data;
p0=6;
p1=7;
p2=2;
p3=9;
x1=[7 9 12 5 4];
x2=[1 8 21 3 5];
x3=[3 2 11 4 8];
x1_mean=mean(x1)
x1_max=max(x1)
x1_min=min(x1)
x1=(x1-x1_mean)/(x1_max-x1_min)
x2_mean=mean(x2)
x2_max=max(x2)
x2_min=min(x2)
x2=(x2-x2_mean)/(x2_max-x2_min)
x3_mean=mean(x3)
x3_max=max(x3)
x3_min=min(x3)
x3=(x3-x3_mean)/(x3_max-x3_min)
y=p0+p1*x1+p2*x2+p3*x3;
num_sample=size(y,2);
% gradient descending process
% initial values of parameters theta0=9;
theta1=3;
theta2=9;
% theta0=19;theta1=23;theta2=91;
% theta0=0;theta1=0;theta2=0;
%learning rate
alpha=1.8; % good for the system
% alpha=0.01; % good for the system
% alpha=0.02; % bad for the system. the system will be unstable
% if alpha is too large, the final error will be much large.
% if alpha is too small, the convergence will be slow
epoch=2260;
for k=1:epoch
v_k=k
h_theta_x=theta0+theta1*x1+theta2*x2+theta3*x3; % hypothesis function Jcost(k)=sum((h_theta_x-y).^2)/num_sample;
r0=sum(h_theta_x-y);
theta0=theta0-alpha*r0/num_sample;
r1=sum((h_theta_x-y).*x1);
theta1=theta1-alpha*r1/num_sample;
r2=sum((h_theta_x-y).*x2);
theta2=theta2-alpha*r2/num_sample;
r3=sum((h_theta_x-y).*x3);
theta3=theta3-alpha*r3/num_sample;
end
yt=theta0+theta1*x1+theta2*x2+theta3*x3
plot(Jcost)
对输入值向量进行归一化处理之后，将学习速率从0.01提高到1.9,从而将获得系统收敛的epoch次数由10000次减少到300次。
图论算法及其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算法每对顶点之间的最短路径
相关主题
文本预览