BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
$ X* e1 u' P0 S    !!!输入函数信息，输出函数的稳定点及迭代次数；
6 a4 ^7 ]' p3 J# e' D, P8 {% j    !!!iter整型变量，存放迭代次数；
2 b3 Q6 C" y+ W, L- L* h# N4 j    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
    !!!dir实型变量，存放搜索方向；
    program main
    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
) V  n- n% k" Z" T( ^$ _4 O    real::x0,tol
    integer::n ,iter,i,j
, O0 V9 Q* a# q% V  F0 @% p    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
" C! P) h/ O  A3 P  o0 g$ c    print*,'请输入hessin矩阵'
    read*,hessin
    print*,'请输入矩阵b'
; x& v4 i, A4 p! T0 U    read*,b
8 M3 c! b* K5 E0 }5 s6 D$ E+ P    iter=0
2 g/ K7 c2 ]# f$ W tol=0.00001</P>
4 J% X2 F9 ~! }+ k* o) [. k<> do i=1,n
2 B* m# i3 n, P" N6 H    do j=1,n
       if (i==j)then
$ R% u5 R8 K$ B       B1(i,j)=1
, t9 C  V7 E  [2 R& T    else
& e9 c, [) V' i! Q0 A# A: Y       B1(i,j)=0
/ B, u, Z. Y* m* O    endif
    enddo
0 i1 I6 U5 H) s" P# x6 k, G% r1 _ enddo   
    gradt=matmul(hessin,x)+b
100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
        !print*,'极小值点为：',x
     !print*,'迭代次数:',iter
     goto 101
/ v4 X) Y6 f6 f    endif
5 w& j/ P! a& u+ @ call gaussj(B1,n,(-1)*gradt)
dir=gradt
    x0=golden(x,dir,hessin,b)
    x1=x+x0*dir
gradt1=matmul(hessin,x1)+b
; |3 y( y% r3 a( [6 j; B) k s=x1-x
5 X+ R6 V% p1 _0 F& S! D y=gradt1-gradt
, P+ a+ U4 t4 G7 I call vectorm(gradt,G)
G1=G
call vectorm(y,G)
) V4 @$ J- a3 G  { B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
x=x1
gradt=gradt1
$ r1 K7 G# [  _% y    iter=iter+1
. B0 B6 o& a" _* X  if(iter&gt;10*n)then
  z# l! P$ Z2 K& h9 m8 |    print*,"out"
    goto 101
% E7 {1 [7 p$ m' ~, h endif
6 |' W" Q9 ]' s; _& ^9 d; p8 _    print*,"第",iter,"次运行结果为",x
+ j6 J( U. B9 E" z$ t print*,"方向为",dir  
  t) t8 [  U, f2 C% w* x: \) s    goto 100
    contains</P>
<>    !!!子程序，返回函数值   
    function f(x,A,b) result(f_result)
& B* J- I" p& G! S) ]    real,dimension(,intent(in)::x,b
" T0 s0 e$ A7 H) a1 T0 }4 S3 q    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)
5 d  u; ]$ O+ U/ l" Z0 _. A real,dimension(,intent(in)::p
8 A- W( D# L* C9 y0 ` real,dimension(:,,intent(out)::G
n=size(p)
8 {5 c& u3 a9 B# H do i=1,n
6 i' q+ y  N: Q1 W& `    !do j=1,n
) N* S; W- X! k* y, p0 R       G(i,=p(i)*p
" Z! M; m0 X3 {$ Z, O0 q    !enddo
enddo
end subroutine

* h  D4 u/ J$ s% J, s) F& R, ?    !!!精确线搜索0.618法子程序 ，返回步长；
    function golden(x,d,A,b) result(golden_n)
& |5 X& L/ Y( p. `3 O# O$ d    real::golden_n
    real::x0
    real,dimension(,intent(in)::x,d
! m+ E% }6 A0 \/ {0 m    real,dimension(,intent(in)::b
    real,dimension(:,,intent(in)::A
3 C& L& _) ^4 J% H2 ~& X7 I    real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
    parameter(r=0.618)
5 m) u8 D  P- H. q$ p9 e: h  o    tol=0.0001
, Z1 y9 j- q( _  x% V* F7 R    dx=0.1
    x0=1
    x1=x0+dx
    f0=f(x+x0*d,A,b)
& C& a9 y/ M4 w    f1=f(x+x1*d,A,b)
  z3 i/ c/ y9 A6 S. U4 O    if(f0&lt;f1)then
/ Z9 i9 w% U# C% w! D& k4       dx=dx+dx
        x2=x0-dx
8 M* y5 \/ F# ?6 R; }- }5 |, P        f2=f(x+x2*d,A,b)
- D4 O: ]$ F  e) F5 y4 Z: P        if(f2&lt;f0)then
           x1=x0
        x0=x2
" z9 s9 ~6 Q% ?        f1=f0
        f0=f2
* ^5 B4 |' |7 Y. |4 F# Z/ T        goto 4
/ o* J+ g) m  W" n        else
           a1=x2
0 Z2 K9 ?' T, `) @6 v        b1=x1
2 N3 G" I7 o. p1 \        endif
6 n  F( z8 u& t3 Y7 n3 w. \5 t    else
2       dx=dx+dx
( M) V0 t/ a) \+ {0 h  Q* W        x2=x1+dx
2 e9 Q) u4 p! M9 _# ]' I        f2=f(x+x2*d,A,b)
        if(f2&gt;=f1)then
! ^% w( m1 N+ y- Y# v5 \- z           b1=x2
        a1=x0
- a7 x. j: M8 {* P) [        else
           x0=x1
        x1=x2
( H0 G% L! [  R        f0=f1
        f1=f2
        goto 2
( q2 }8 R5 L' m8 q& k        endif
7 N2 W" F; K$ K9 l& `    endif
    x1=a1+(1-r)*(b1-a1)
    x2=a1+r*(b1-a1)
- A  l# I. b- }7 c    f1=f(x+x1*d,A,b)
    f2=f(x+x2*d,A,b)
" M" v) x9 d, j# t6 E* d8 X. Y3   if(abs(b1-a1)&lt;=tol)then
        x0=(a1+b1)/2
; Q% \$ Q( I0 [) X; q3 i% w8 C% L8 v    else
2 p+ U; E+ t4 ?& J+ x& T& B        if(f1&gt;f2)then
3 g7 ^4 X8 h2 t        a1=x1
+ |& A+ @! O( W1 Q" @# q3 G        x1=x2
        f1=f2
        x2=a1+r*(b1-a1)
        f2=f(x+x2*d,A,b)
: e7 ?1 l8 J# q1 |        goto 3
3 ?6 _) J6 d8 c6 q  I) V     else
        b1=x2
        x2=x1
+ A  l5 j$ \# H- j        f2=f1
        x1=a1+(1-r)*(b1-a1)
4 _7 x( F" @0 o        f1=f(x+x1*d,A,b)
        goto 3
     endif
    endif
" [, e$ e. B% y; s" s+ Z    golden_n=x0
  l, t3 `* u, V* C/ m6 {    end  function golden</P>
+ ]# A4 W' t$ Z8 e& C5 x4 U, y<>
2 F6 I8 r. x) e    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
/ O0 Y& V' ^5 b" r! S6 A% y% G    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)
5 d( H) E. n7 I% P: N! [    real big,dum,pivinv  
2 Q# t3 B( i) y    do j=1,n
       ipiv(j)=0
% }: l2 J2 l% V6 T- ~& u    enddo
- b( N% `0 o. E, p' q    do i=1,n
+ W0 W* a/ {( {% B! g% k       big=0.
       do j=1,n
       if(ipiv(j)/=1)then
          do k=1,n
          if(ipiv(k)==0)then
8 F9 ^5 g, y  [; Z1 t          if(abs(a(j,k))&gt;=big)then
           big=abs(a(j,k))
           irow=j
           icol=k
       endif
4 _- E% i7 e. d, p3 o" P       else if(ipiv(k)&gt;1)then
          pause'singular matrix in gaussj'
       endif
9 z0 T$ a- a, h& J) ?0 ?2 |9 E       enddo
& y- {+ _; C! B/ p+ _: g    endif
    enddo
9 n* u' i" l& J3 d# E    ipiv(icol)=ipiv(icol)+1
1 h9 }9 K$ _, J: c' s4 l    if(irow/=icol)then
       do l=1,n
; C2 S* w  L( b! G/ f7 Q) X          dum=a(irow,l)
7 k6 E5 S3 q- N7 i/ z( V; B6 e9 n4 h       a(irow,l)=a(icol,l)
5 Y9 D" q% j% J$ {- _6 B0 S       a(icol,l)=dum
       enddo
       dum=b(irow)
       b(irow)=b(icol)
       b(icol)=dum
5 A! q2 |2 A: _0 \; e    endif
: ^; a! t2 Y' l7 a4 t" x* u    indxr(i)=irow
! \6 G& ]5 u' Y3 S    indxc(i)=icol
5 d7 ]2 d$ W0 F* C    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
3 l5 R( e/ w. b    pivinv=1./a(icol,icol)
    a(icol,icol)=1.
8 R0 i& a3 J* j* I* [% @0 W1 E& Y    do l=1,n
        a(icol,l)=a(icol,l)*pivinv
    enddo
    b(icol)=b(icol)*pivinv
    do ll=1,n
" D  g; z+ N# U, w6 F+ B       if(ll/=icol)then
          dum=a(ll,icol)
& h! B" W* E% V+ ~       a(ll,icol)=0
       do l=1,n
          a(ll,l)=a(ll,l)-a(icol,l)*dum
       enddo
. V5 r8 ]+ V& Z       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))
; y* |+ B3 h/ K" H' l# p- o       a(k,indxr(l))=a(k,indxc(l))
: _8 }$ s; j# \       a(k,indxc(l))=dum
, B* c. }4 l# S' W3 X6 [7 }       enddo
" z" @, j2 T6 s4 j$ `0 Z    endif
) F& F! s& u  x2 o    enddo
2 {2 M2 m1 N' y* |4 K    end subroutine gaussj
101 end
; Q2 f) q4 D/ S/ \</P>
/ e& D4 v8 |# m( i& G1 A4 E<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

楼主，怎么这么多笑脸，晕了，有没有matlab的