BFGS算法

ilikenba

<>    !!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点；
: M9 p' K: E5 t& c    !!!输入函数信息，输出函数的稳定点及迭代次数；
! m$ j  V0 B! e4 t    !!!iter整型变量，存放迭代次数；
    !!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度;
    !!!dir实型变量，存放搜索方向；
    program main
3 h; J2 o; z( T    real,dimension(,allocatable::x,gradt,dir,b ,s,y,gradt1,x1
: `7 f; v' U0 n6 w7 Y    real,dimension(:,,allocatable::hessin ,B1 ,G,G1
8 C- z4 ]2 H% ]+ ?    real::x0,tol
    integer::n ,iter,i,j
    print*,'请输入变量的维数'
    read*,n
0 A' P+ k: G; b    allocate(x(n),gradt(n),dir(n),b(n),s(n),y(n),gradt1(n),x1(n))
! O. }) _% t+ i# q. @    allocate(hessin(n,n),B1(n,n),G(n,n),G1(n,n))
: C: g2 |+ B2 n: B7 A1 @# F    print*,'请输入初始向量x'
    read*,x
    print*,'请输入hessin矩阵'
    read*,hessin
    print*,'请输入矩阵b'
9 q& P- l' i) \& j4 J    read*,b
- j3 f. G  ~. X6 e: F; D- f    iter=0
0 l" k: a/ f) V; i% o" u tol=0.00001</P>
<> do i=1,n
    do j=1,n
9 S9 |: O6 \5 S7 M0 Z       if (i==j)then
1 g) L2 U% P0 b8 p8 b; D       B1(i,j)=1
+ D, u2 t9 ?5 p: A    else
1 W7 H5 }8 b1 F0 C, C' E7 g       B1(i,j)=0
( R1 n9 }$ F1 `( B6 E2 q7 X) X    endif
) ^& f  |* z8 @7 Q! k0 v/ b* F. B    enddo
5 I! c+ G* y0 V- W8 A enddo   
& P8 K2 B& G% C" N! j( h    gradt=matmul(hessin,x)+b
! k, w: a2 C3 j. x- o100 if(sqrt(dot_product(gradt,gradt))&lt;tol)then
4 F; i& N5 U2 O. u' o* z: e1 @( s        !print*,'极小值点为：',x
' d5 u, ~  l; P     !print*,'迭代次数:',iter
8 q  N% l" \/ Q     goto 101
5 c7 i- B2 [3 a& I# D. w( n5 Y    endif
call gaussj(B1,n,(-1)*gradt)
dir=gradt
' V2 ^$ W' e6 j& q    x0=golden(x,dir,hessin,b)
( s; s+ W& }% `: c- o% k) X    x1=x+x0*dir
gradt1=matmul(hessin,x1)+b
: h: g+ [5 p6 @5 \! {. Z s=x1-x
y=gradt1-gradt
call vectorm(gradt,G)
* F5 J9 t$ P9 N# U' s G1=G
call vectorm(y,G)
9 ]6 `/ V/ G3 ~& m9 t# p B1=B1+1/dot_product(gradt,dir)*G1+1/(x0*dot_product(y,dir))*G
5 t8 L* r% S2 j7 t( n* d x=x1
, G; I+ {# Q8 D5 b8 R% s gradt=gradt1
    iter=iter+1
- k+ c( U8 Z! m% o: G( C$ E  if(iter&gt;10*n)then
) [  S- D' M1 M! S+ ?    print*,"out"
/ L# }( B2 H& i3 C5 T4 D    goto 101
; q3 a9 y' G1 H4 F5 ^ endif
    print*,"第",iter,"次运行结果为",x
) t5 e6 _6 ?* t9 r. B* K! Z5 u/ O print*,"方向为",dir  
/ U6 i3 Y  ?) r- `  F    goto 100
    contains</P>
<>    !!!子程序，返回函数值   
    function f(x,A,b) result(f_result)
    real,dimension(,intent(in)::x,b
; V  C# }0 K6 }$ m$ l4 r6 @) p    real,dimension(:,,intent(in)::A
" C9 G/ k: s; [& @    real::f_result
    f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x)
2 b- y4 a8 h5 ^# B    end function f
7 f" k1 O3 D% H: _. d !!!子程序，矩阵与向量相乘
subroutine vectorm(p,G)
' t5 K! O5 H( ?+ ? real,dimension(,intent(in)::p
real,dimension(:,,intent(out)::G
* A/ o' p: s& t8 n4 ^& } n=size(p)
do i=1,n
2 U1 l" ^; l' O& y! K: X" ?    !do j=1,n
       G(i,=p(i)*p
( S" D: b) M2 ?- e    !enddo
enddo
5 c+ x1 [5 O$ g end subroutine

    !!!精确线搜索0.618法子程序 ，返回步长；
+ R# ~) F+ u7 X0 `6 V    function golden(x,d,A,b) result(golden_n)
* ]" l6 w2 n' a. q8 `+ H    real::golden_n
/ s. F$ W5 y# ~3 i9 D    real::x0
    real,dimension(,intent(in)::x,d
; D6 E4 ^& ~4 e% V$ D1 x* N! M; }    real,dimension(,intent(in)::b
    real,dimension(:,,intent(in)::A
    real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx
* A# b- a0 s; i# c& j    parameter(r=0.618)
4 l0 X# G2 x+ w% S/ K$ N. F    tol=0.0001
    dx=0.1
% c$ I& E* {  U3 u% ~8 y; x    x0=1
    x1=x0+dx
    f0=f(x+x0*d,A,b)
    f1=f(x+x1*d,A,b)
0 \, L7 }+ C& ~8 I    if(f0&lt;f1)then
4       dx=dx+dx
        x2=x0-dx
% b, u4 ^5 M0 _7 I2 D, p7 ~/ T2 o        f2=f(x+x2*d,A,b)
9 {' U+ E+ x8 R        if(f2&lt;f0)then
           x1=x0
7 B: h, z2 x" s$ R; [+ `: {        x0=x2
, B7 L4 L! M' s0 B        f1=f0
        f0=f2
        goto 4
6 a( {; v! b6 n+ `; s        else
- V& d2 U0 ?  k/ g           a1=x2
        b1=x1
- ^; c* U$ G" @. W( a4 {        endif
    else
4 j: m8 ~! T- u+ F$ V2       dx=dx+dx
        x2=x1+dx
        f2=f(x+x2*d,A,b)
0 U0 _: F6 ]& a8 N        if(f2&gt;=f1)then
4 U: M8 f' H; O7 V           b1=x2
        a1=x0
        else
% ^# u2 {) Q! v# i% K           x0=x1
' K/ u% m' M9 g        x1=x2
. ]4 m' ]% @% Y2 X: Y/ o2 e6 _5 r* y        f0=f1
# G  z3 {, @' K7 ?5 }  W        f1=f2
        goto 2
        endif
% C' I2 l5 J. v3 U# c9 u0 f. ]    endif
4 [6 l# H3 z  _    x1=a1+(1-r)*(b1-a1)
% U+ c! d" L& y; c4 e7 z    x2=a1+r*(b1-a1)
    f1=f(x+x1*d,A,b)
7 q. ~4 F2 [/ n" x    f2=f(x+x2*d,A,b)
3   if(abs(b1-a1)&lt;=tol)then
        x0=(a1+b1)/2
2 k' \: q& R  p( m# |    else
        if(f1&gt;f2)then
        a1=x1
8 d% H  q6 K( ~# J        x1=x2
9 [5 q& `, a) C/ V        f1=f2
        x2=a1+r*(b1-a1)
        f2=f(x+x2*d,A,b)
2 N7 A7 P/ |; N1 h( }        goto 3
     else
        b1=x2
/ p& v6 }, Y% R' w" v        x2=x1
        f2=f1
        x1=a1+(1-r)*(b1-a1)
$ J9 z# b/ b/ [4 Z        f1=f(x+x1*d,A,b)
2 y5 Q3 i) e1 B: W7 f8 j/ _        goto 3
     endif
, w" y6 `; e; `6 E4 _/ i  }( g    endif
    golden_n=x0
0 y7 W/ @: N$ j$ l1 z" k    end  function golden</P>
<>
    !!!A为二维数，返回值为的A逆;b为一维数组，返回值为方程Ax=b的解
; `# T. N7 U& `! C# h' D) y2 u    subroutine gaussj(a,n,b)
5 B( b9 [  u, \, \5 O) O    integer n,nmax
    real a(n,n),b(n)
    parameter(nmax=50)
# v# Y& y5 P( C. {    integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
, }' G2 \) ]7 @4 G% z9 `    real big,dum,pivinv  
    do j=1,n
       ipiv(j)=0
    enddo
    do i=1,n
0 \5 b. I4 v+ ~6 v; L, W* s8 }3 c       big=0.
       do j=1,n
       if(ipiv(j)/=1)then
8 U! c8 H, _. S7 V1 |2 p: d          do k=1,n
% B. D$ f; m  w) a5 @9 b          if(ipiv(k)==0)then
          if(abs(a(j,k))&gt;=big)then
$ Z- `) E# @" _/ V$ t           big=abs(a(j,k))
& d4 {) U1 x" {% q7 \  U           irow=j
           icol=k
       endif
9 N- B0 T9 w9 \6 `) s6 z       else if(ipiv(k)&gt;1)then
          pause'singular matrix in gaussj'
       endif
       enddo
    endif
* w7 n0 r: n# R- n    enddo
    ipiv(icol)=ipiv(icol)+1
) o9 Z8 l' _1 V1 b0 N    if(irow/=icol)then
8 J# S/ o/ G' q( a; q/ V       do l=1,n
8 n2 y; d& k$ M          dum=a(irow,l)
       a(irow,l)=a(icol,l)
' r* T, P3 u3 _3 n5 H/ j       a(icol,l)=dum
       enddo
7 d* h, C  Z9 f; y/ a2 j       dum=b(irow)
       b(irow)=b(icol)
       b(icol)=dum
2 Q* w- T) r+ q- O# E- a/ a# z    endif
    indxr(i)=irow
    indxc(i)=icol
    if(a(icol,icol)==0.)pause'singular matrix in gaussj'
- U) [7 L6 Y0 Y( C$ X6 O+ Y7 s    pivinv=1./a(icol,icol)
5 ^+ u. K2 r5 O& x: ]5 C7 K/ [    a(icol,icol)=1.
. _! U$ ?( J( a! c    do l=1,n
        a(icol,l)=a(icol,l)*pivinv
    enddo
6 V( S* ]- ?. t6 t: r    b(icol)=b(icol)*pivinv
    do ll=1,n
1 K* ]: q; r" v! Z1 z% ~       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
, @/ i& l9 a5 G3 r$ l  l% W    enddo
, c: e/ o. V, @    do l=n,1,-1
% S; D# B7 D* |       if(indxr(l)/=indxc(l))then
: f8 I# s8 P) ^( N5 [. a       do k=1,n
          dum=a(k,indxr(l))
       a(k,indxr(l))=a(k,indxc(l))
       a(k,indxc(l))=dum
  w' ]# R0 K5 O, N2 u3 k       enddo
    endif
) n( ]6 S2 }# n+ _7 P) u    enddo
9 S: B  s. q" o! E    end subroutine gaussj
+ I! v8 @# t. P0 ]101 end
* U* W5 V* m4 J) f</P>
<>本程序用Fortran 90编写，在Virual Fortran 5上编译通过，本程序由沙沙提供！</P>

hxjean

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