subsm – lbfgsb

private subroutine subsm(n, m, Nsub, Ind, l, u, Nbd, x, d, Xp, Ws, Wy, Theta, Xx, Gg, Col, Head, Iword, Wv, Wn, Iprint, Info)

This routine contains the major changes in the updated version. The changes are described in the accompanying paper

Given xcp, l, u, r, an index set that specifies the active set at xcp, and an l-BFGS matrix B (in terms of WY, WS, SY, WT, head, col, and theta), this subroutine computes an approximate solution of the subspace problem

(P) min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp)

    subject to l<=x<=u
              x_i=xcp_i for all i in A(xcp)

along the subspace unconstrained Newton direction

 d = -(Z'BZ)^(-1) r.

The formula for the Newton direction, given the L-BFGS matrix and the Sherman-Morrison formula, is

 d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r.

where K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] [L_a -R_z theta*S'AA'S ]

Note that this procedure for computing d differs from that described in [1]. One can show that the matrix K is equal to the matrix M^[-1]N in that paper.

References

R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, "A limited memory algorithm for bound constrained optimization", SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.

Credits

NEOS, November 1994. (Latest revision June 1996.) Optimization Technology Center. Argonne National Laboratory and Northwestern University. Written by Ciyou Zhu in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
Jose Luis Morales, Jorge Nocedal "Remark On Algorithm 788: L-BFGS-B: Fortran Subroutines for Large-Scale Bound Constrained Optimization". Decemmber 27, 2010.
J.L. Morales & J, Nocedal, January 17, 2011

Arguments

Type	Intent		Name
integer,	intent(in)	::	n	the dimension of the problem
integer,	intent(in)	::	m	the maximum number of variable metric corrections used to define the limited memory matrix.
integer,	intent(in)	::	Nsub	the number of free variables
integer,	intent(in)	::	Ind(Nsub)	specifies the coordinate indices of free variables.
real(kind=wp),	intent(in)	::	l(n)	the lower bound of x.
real(kind=wp),	intent(in)	::	u(n)	the upper bound of x
integer,	intent(in)	::	Nbd(n)	represents the type of bounds imposed on the variables, and must be specified as follows: `nbd(i)=0` if `x(i)` is unbounded, `nbd(i)=1` if `x(i)` has only a lower bound, `nbd(i)=2` if `x(i)` has both lower and upper bounds, and `nbd(i)=3` if `x(i)` has only an upper bound.
real(kind=wp),	intent(inout)	::	x(n)	On entry x specifies the Cauchy point xcp. On exit x(i) is the minimizer of Q over the subspace of free variables.
real(kind=wp),	intent(inout)	::	d(n)	On entry d is the reduced gradient of Q at xcp. On exit d is the Newton direction of Q.
real(kind=wp)		::	Xp(n)	used to safeguard the projected Newton direction
real(kind=wp),	intent(in)	::	Ws(n,m)	variable that stores the information defining the limited memory BFGS matrix: ws(n,m) stores S, a set of s-vectors
real(kind=wp),	intent(in)	::	Wy(n,m)	variable that stores the information defining the limited memory BFGS matrix:
real(kind=wp),	intent(in)	::	Theta	variable that stores the information defining the limited memory BFGS matrix: theta is the scaling factor specifying B_0 = theta I
real(kind=wp),	intent(in)	::	Xx(n)	the current iterate
real(kind=wp),	intent(in)	::	Gg(n)	the gradient at the current iterate
integer,	intent(in)	::	Col	variable that stores the information defining the limited memory BFGS matrix: col is the number of variable metric corrections stored
integer,	intent(in)	::	Head	variable that stores the information defining the limited memory BFGS matrix: head is the location of the 1st s- (or y-) vector in S (or Y)
integer,	intent(out)	::	Iword	On exit iword specifies the status of the subspace solution: `iword = 0` if the solution is in the box `iword = 1` if some bound is encountered
real(kind=wp)		::	Wv(2*m)	a real(wp) working array of dimension 2m
real(kind=wp),	intent(in)	::	Wn(2m,2m)	The upper triangle of wn stores the LEL^T factorization of the indefinite matrix `K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] [L_a -R_z theta*S'AA'S ] where E = [-I 0] [ 0 I]`
integer,	intent(in)	::	Iprint	If iprint >= 99, some minimum debug printing is done.
integer,	intent(out)	::	Info	On exit: `info = 0` for normal return, `info /= 0` for abnormal return when the matrix `K` is ill-conditioned.

Calls

Help

Called by

Help

Source Code

      subroutine subsm(n,m,Nsub,Ind,l,u,Nbd,x,d,Xp,Ws,Wy,Theta,Xx,Gg, &
                       Col,Head,Iword,Wv,Wn,Iprint,Info)

      implicit none

      integer,intent(in) :: n !! the dimension of the problem
      integer,intent(in) :: m !! the maximum number of variable metric corrections
                              !! used to define the limited memory matrix.
      integer,intent(in) :: Nsub !! the number of free variables
      real(wp),intent(in) :: Ws(n,m) !! variable that stores the information defining the limited memory BFGS matrix:
                                     !! ws(n,m) stores S, a set of s-vectors
      real(wp),intent(in) :: Wy(n,m) !! variable that stores the information defining the limited memory BFGS matrix:
                                     ! !wy(n,m) stores Y, a set of y-vectors
      real(wp),intent(in) :: Theta !! variable that stores the information defining the limited memory BFGS matrix:
                                   !! theta is the scaling factor specifying B_0 = theta I
      integer,intent(in) :: Col !! variable that stores the information defining the limited memory BFGS matrix:
                                !! col is the number of variable metric corrections stored
      integer,intent(in) :: Head !! variable that stores the information defining the limited memory BFGS matrix:
                                 !! head is the location of the 1st s- (or y-) vector in S (or Y)
      integer,intent(out) :: Iword !! On exit iword specifies the status of the subspace solution:
                                   !!
                                   !!  * `iword = 0` if the solution is in the box
                                   !!  * `iword = 1` if some bound is encountered
      integer,intent(in) :: Iprint !! If iprint >= 99, some minimum debug printing is done.
      integer,intent(out) :: Info !! On exit:
                                  !!
                                  !!  * `info = 0`  for normal return,
                                  !!  * `info /= 0` for abnormal return
                                  !!    when the matrix `K` is ill-conditioned.
      integer,intent(in) :: Ind(Nsub) !! specifies the coordinate indices of free variables.
      integer,intent(in) :: Nbd(n) !! represents the type of bounds imposed on the
                                   !! variables, and must be specified as follows:
                                   !!
                                   !!  * `nbd(i)=0` if `x(i)` is unbounded,
                                   !!  * `nbd(i)=1` if `x(i)` has only a lower bound,
                                   !!  * `nbd(i)=2` if `x(i)` has both lower and upper bounds, and
                                   !!  * `nbd(i)=3` if `x(i)` has only an upper bound.
      real(wp),intent(in) :: l(n) !! the lower bound of x.
      real(wp),intent(in) :: u(n) !! the upper bound of x
      real(wp),intent(inout) :: x(n) !! On entry x specifies the Cauchy point xcp.
                                     !! On exit x(i) is the minimizer of Q over the subspace of
                                     !! free variables.
      real(wp),intent(inout) :: d(n) !! On entry d is the reduced gradient of Q at xcp.
                                     !! On exit d is the Newton direction of Q.
      real(wp) :: Xp(n) !! used to safeguard the projected Newton direction
      real(wp),intent(in) :: Xx(n) !! the current iterate
      real(wp),intent(in) :: Gg(n) !! the gradient at the current iterate
      real(wp) :: Wv(2*m) !! a real(wp) working array of dimension 2m
      real(wp),intent(in) :: Wn(2*m,2*m) !! The upper triangle of wn stores the LEL^T factorization
                                         !! of the indefinite matrix
                                         !!```
                                         !!  K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
                                         !!      [L_a -R_z           theta*S'AA'S ]
                                         !!  where E = [-I  0]
                                         !!            [ 0  I]
                                         !!```

      integer :: pointr , m2 , col2 , ibd , jy , js , i , j , k
      real(wp) :: alpha , xk , dk , temp1 , temp2
      real(wp) :: dd_p

      if ( Nsub<=0 ) return
      if ( Iprint>=99 ) write (output_unit,'(/,A,/)') '----------------SUBSM entered-----------------'

      ! Compute wv = W'Zd.

      pointr = Head
      do i = 1 , Col
         temp1 = zero
         temp2 = zero
         do j = 1 , Nsub
            k = Ind(j)
            temp1 = temp1 + Wy(k,pointr)*d(j)
            temp2 = temp2 + Ws(k,pointr)*d(j)
         enddo
         Wv(i) = temp1
         Wv(Col+i) = Theta*temp2
         pointr = mod(pointr,m) + 1
      enddo

      ! Compute wv:=K^(-1)wv.

      m2 = 2*m
      col2 = 2*Col
      call dtrsl(Wn,m2,col2,Wv,11,Info)
      if ( Info/=0 ) return
      do i = 1 , Col
         Wv(i) = -Wv(i)
      enddo
      call dtrsl(Wn,m2,col2,Wv,01,Info)
      if ( Info/=0 ) return

      ! Compute d = (1/theta)d + (1/theta**2)Z'W wv.

      pointr = Head
      do jy = 1 , Col
         js = Col + jy
         do i = 1 , Nsub
            k = Ind(i)
            d(i) = d(i) + Wy(k,pointr)*Wv(jy)/Theta + Ws(k,pointr)*Wv(js)
         enddo
         pointr = mod(pointr,m) + 1
      enddo

      call dscal(Nsub,one/Theta,d,1)

      !-----------------------------------------------------------------
      ! Let us try the projection, d is the Newton direction

      Iword = 0

      call dcopy(n,x,1,Xp,1)

      do i = 1 , Nsub
         k = Ind(i)
         dk = d(i)
         xk = x(k)
         if ( Nbd(k)/=0 ) then

            if ( Nbd(k)==1 ) then          ! lower bounds only
               x(k) = max(l(k),xk+dk)
               if ( x(k)==l(k) ) Iword = 1
            else

               if ( Nbd(k)==2 ) then       ! upper and lower bounds
                  xk = max(l(k),xk+dk)
                  x(k) = min(u(k),xk)
                  if ( x(k)==l(k) .or. x(k)==u(k) ) Iword = 1
               else

                  if ( Nbd(k)==3 ) then    ! upper bounds only
                     x(k) = min(u(k),xk+dk)
                     if ( x(k)==u(k) ) Iword = 1
                  endif
               endif
            endif

         else                                ! free variables
            x(k) = xk + dk
         endif
      enddo

      main : block

         if ( Iword==0 ) exit main

         ! check sign of the directional derivative

         dd_p = zero
         do i = 1 , n
            dd_p = dd_p + (x(i)-Xx(i))*Gg(i)
         enddo
         if ( dd_p<=zero ) exit main

         call dcopy(n,Xp,1,x,1)
         if ( Iprint>=0 ) write (output_unit,'(A)') ' Positive dir derivative in projection '
         if ( Iprint>=0 ) write (output_unit,'(A)') ' Using the backtracking step '

         !-----------------------------------------------------------------

         alpha = one
         temp1 = alpha
         ibd = 0
         do i = 1 , Nsub
            k = Ind(i)
            dk = d(i)
            if ( Nbd(k)/=0 ) then
               if ( dk<zero .and. Nbd(k)<=2 ) then
                  temp2 = l(k) - x(k)
                  if ( temp2>=zero ) then
                     temp1 = zero
                  else if ( dk*alpha<temp2 ) then
                     temp1 = temp2/dk
                  endif
               else if ( dk>zero .and. Nbd(k)>=2 ) then
                  temp2 = u(k) - x(k)
                  if ( temp2<=zero ) then
                     temp1 = zero
                  else if ( dk*alpha>temp2 ) then
                     temp1 = temp2/dk
                  endif
               endif
               if ( temp1<alpha ) then
                  alpha = temp1
                  ibd = i
               endif
            endif
         enddo

         if ( alpha<one ) then
            dk = d(ibd)
            k = Ind(ibd)
            if ( dk>zero ) then
               x(k) = u(k)
               d(ibd) = zero
            else if ( dk<zero ) then
               x(k) = l(k)
               d(ibd) = zero
            endif
         endif
         do i = 1 , Nsub
            k = Ind(i)
            x(k) = x(k) + alpha*d(i)
         enddo

      end block main

      if ( Iprint>=99 ) write (output_unit,'(/,A,/)') '----------------exit SUBSM --------------------'

      end subroutine subsm

subsm Subroutine

Contents

Source Code

private subroutine subsm(n, m, Nsub, Ind, l, u, Nbd, x, d, Xp, Ws, Wy, Theta, Xx, Gg, Col, Head, Iword, Wv, Wn, Iprint, Info)

References

Credits

Arguments

Calls

Called by

Source Code