查看: 5109|回复: 1

BFGS算法

[复制链接]

字体大小: 正常放大

ilikenba

1万主题	49 听众	2万积分

TA的每日心情

	奋斗 2024-6-23 05:14

签到天数: 1043 天

[LV.10]以坛为家III

群组: 万里江山

群组: sas讨论小组

群组: 长盛证券理财有限公司

群组: C 语言讨论组

群组: Matlab讨论组

电梯直达

1^#

发表于 2004-4-30 10:51 |只看该作者 |倒序浏览

|招呼Ta 关注Ta

<

> !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
!!!输入函数信息，输出函数的稳定点及迭代次数；
!!!iter整型变量，存放迭代次数；
!!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
!!!dir实型变量，存放搜索方向；
program main
real,dimension(

,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
real,dimension(:,

,allocatable::hessin ,B1 ,G,G1
real::x0,tol
integer::n ,iter,i,j
print*,'请输入变量的维数'
read*,n
allocate(x(n),gradt(n),dir(n),b(n),s(n),y(n),gradt1(n),x1(n))
allocate(hessin(n,n),B1(n,n),G(n,n),G1(n,n))
print*,'请输入初始向量x'
read*,x
print*,'请输入hessin矩阵'
read*,hessin
print*,'请输入矩阵b'
read*,b
iter=0
tol=0.00001
<

> do i=1,n
do j=1,n
 if (i==j)then
 B1(i,j)=1
else
 B1(i,j)=0
endif
enddo
enddo
gradt=matmul(hessin,x)+b
100 if(sqrt(dot_product(gradt,gradt))<tol)then
 !print*,'极小值点为：',x
 !print*,'迭代次数:',iter
 goto 101
endif
call gaussj(B1,n,(-1)*gradt)
dir=gradt
x0=golden(x,dir,hessin,b)
x1=x+x0*dir
gradt1=matmul(hessin,x1)+b
s=x1-x
y=gradt1-gradt
call vectorm(gradt,G)
G1=G
call vectorm(y,G)
B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
x=x1
gradt=gradt1
iter=iter+1
 if(iter>10*n)then
print*,"out"
goto 101
endif
print*,"第",iter,"次运行结果为",x
print*,"方向为",dir
goto 100
contains
<

> !!!子程序，返回函数值
function f(x,A,b) result(f_result)
real,dimension(

,intent(in)::x,b
real,dimension(:,

,intent(in)::A
real::f_result
f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
end function f
!!!子程序，矩阵与向量相乘
subroutine vectorm(p,G)
real,dimension(

,intent(in)::p
real,dimension(:,

,intent(out)::G
n=size(p)
do i=1,n
!do j=1,n
G(i,

=p(i)*p
!enddo
enddo
end subroutine

!!!精确线搜索0.618法子程序，返回步长；
function golden(x,d,A,b) result(golden_n)
real::golden_n
real::x0
real,dimension(

,intent(in)::x,d
real,dimension(

,intent(in)::b
real,dimension(:,

,intent(in)::A
real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
parameter(r=0.618)
tol=0.0001
dx=0.1
x0=1
x1=x0+dx
f0=f(x+x0*d,A,b)
f1=f(x+x1*d,A,b)
if(f0<f1)then
4 dx=dx+dx
 x2=x0-dx
 f2=f(x+x2*d,A,b)
 if(f2<f0)then
 x1=x0
 x0=x2
 f1=f0
 f0=f2
 goto 4
 else
 a1=x2
 b1=x1
 endif
else
2 dx=dx+dx
 x2=x1+dx
 f2=f(x+x2*d,A,b)
 if(f2>=f1)then
 b1=x2
 a1=x0
 else
 x0=x1
 x1=x2
 f0=f1
 f1=f2
 goto 2
 endif
endif
x1=a1+(1-r)*(b1-a1)
x2=a1+r*(b1-a1)
f1=f(x+x1*d,A,b)
f2=f(x+x2*d,A,b)
3 if(abs(b1-a1)<=tol)then
 x0=(a1+b1)/2
else
 if(f1>f2)then
 a1=x1
 x1=x2
 f1=f2
 x2=a1+r*(b1-a1)
 f2=f(x+x2*d,A,b)
 goto 3
 else
 b1=x2
 x2=x1
 f2=f1
 x1=a1+(1-r)*(b1-a1)
 f1=f(x+x1*d,A,b)
 goto 3
 endif
endif
golden_n=x0
end function golden
<

>
!!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
subroutine gaussj(a,n,b)
integer n,nmax
real a(n,n),b(n)
parameter(nmax=50)
integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
real big,dum,pivinv
do j=1,n
 ipiv(j)=0
enddo
do i=1,n
 big=0.
 do j=1,n
 if(ipiv(j)/=1)then
 do k=1,n
 if(ipiv(k)==0)then
 if(abs(a(j,k))>=big)then
 big=abs(a(j,k))
 irow=j
 icol=k
 endif
 else if(ipiv(k)>1)then
 pause'singular matrix in gaussj'
 endif
 enddo
endif
enddo
ipiv(icol)=ipiv(icol)+1
if(irow/=icol)then
 do l=1,n
 dum=a(irow,l)
 a(irow,l)=a(icol,l)
 a(icol,l)=dum
 enddo
 dum=b(irow)
 b(irow)=b(icol)
 b(icol)=dum
endif
indxr(i)=irow
indxc(i)=icol
if(a(icol,icol)==0.)pause'singular matrix in gaussj'
pivinv=1./a(icol,icol)
a(icol,icol)=1.
do l=1,n
 a(icol,l)=a(icol,l)*pivinv
enddo
b(icol)=b(icol)*pivinv
do ll=1,n
 if(ll/=icol)then
 dum=a(ll,icol)
 a(ll,icol)=0
 do l=1,n
 a(ll,l)=a(ll,l)-a(icol,l)*dum
 enddo
 b(ll)=b(ll)-b(icol)*dum
 endif
enddo
enddo
do l=n,1,-1
 if(indxr(l)/=indxc(l))then
 do k=1,n
 dum=a(k,indxr(l))
 a(k,indxr(l))=a(k,indxc(l))
 a(k,indxc(l))=dum
 enddo
endif
enddo
end subroutine gaussj
101 end

<

>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！

zan