Gnimuc / LBFGSB.jl

Julia wrapper for L-BFGS-B Nonlinear Optimization Code
MIT License
38 stars 6 forks source link

Use a modified version of L-BFGS-B which is already used in scipy.optimize? #11

Open rht opened 3 years ago

rht commented 3 years ago

Thank you for wrapping the library. Such a lifesaver. Someone in discourse.julialang.org mentioned that the lbfgsb.f used in scipy.optimize is actually a 2011 modification of the original method, see: https://github.com/JuliaNLSolvers/Optim.jl/issues/927.

There is https://github.com/JuliaBinaryWrappers/L_BFGS_B_jll.jl but it is not clear to me how you generated the .so binary.

rht commented 3 years ago

You basically just need to use this: https://github.com/scipy/scipy/blob/master/scipy/optimize/lbfgsb_src/lbfgsb.f and recompile.

Gnimuc commented 3 years ago

There is https://github.com/JuliaBinaryWrappers/L_BFGS_B_jll.jl but it is not clear to me how you generated the .so binary.

The BB script: https://github.com/JuliaPackaging/Yggdrasil/tree/master/L/L_BFGS_B

You basically just need to use this: https://github.com/scipy/scipy/blob/master/scipy/optimize/lbfgsb_src/lbfgsb.f and recompile.

I guess we can apply a patch over the original version, but another thing I'm concerning is the license.

rht commented 3 years ago

You can view the licensing requirement of the 2011 version at https://github.com/scipy/scipy/blob/master/scipy/optimize/lbfgsb_src/README.

I additionally have noticed 2 changes to the API:

Gnimuc commented 3 years ago

That's the license of the original fortran code, I'm not sure whether it applies to the modified version or not. Anyway, feel free to open a PR and ignore the license issue at the moment.

rht commented 3 years ago

That 2011 version already includes the modification. I have inspected the source code of http://www.ece.northwestern.edu/~nocedal/lbfgsb.html. It has the required c-jlm-jn comments.

c===========   L-BFGS-B (version 3.0.  April 25, 2011  ===================
c
c     This is a modified version of L-BFGS-B. Minor changes in the updated 
c     code appear preceded by a line comment as follows 
c  
c     c-jlm-jn 
c

OK, I will make a PR.

Gnimuc commented 3 years ago

Feel free to ping me if you have any other questions :)

rht commented 3 years ago

Ah, I see. https://github.com/JuliaPackaging/Yggdrasil/blob/master/L/L_BFGS_B/build_tarballs.jl#L10 already uses the 3.0 version. I have also inspected the comment, also contains c-jlm-jn. But scipy.optimize's, lbfgsb.f is different since it has the extra maxls argument.

the first argument of setulb, n, is gone

This argument in scipy.optimize's lbfgsb.f is still there. It's just that scipy auto-computes it as len(x).

Gnimuc commented 3 years ago

You could download both the code from the official 3.0 tarball and the one from scipy, and then do a git diff to see the changes.

rht commented 3 years ago

OK, basically the main differences are

Seems like a useful change.

$ diff original.f lbfgsb.f 
1,5d0
< c                                                                                      
< c  L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License”        
< c  or “3-clause license”)                                                              
< c  Please read attached file License.txt                                               
< c                                        
48c43
<      +                 task, iprint, csave, lsave, isave, dsave)
---
>      +                 task, iprint, csave, lsave, isave, dsave, maxls)
52c47
<       integer          n, m, iprint, 
---
>       integer          n, m, iprint, maxls,
143,144c138
< c       When iprint > 0, the file iterate.dat will be created to
< c                        summarize the iteration.
---
> c       
279c273
<      +  csave,lsave,isave(22),dsave)
---
>      +  csave,lsave,isave(22),dsave, maxls)
290c284
<      +                  iprint, csave, lsave, isave, dsave)
---
>      +                  iprint, csave, lsave, isave, dsave, maxls)
295c289,290
<      +                 iwhere(n), indx2(n), isave(23)
---
>      +                 iwhere(n), indx2(n), isave(23),
>      +                 maxls
429,430c424
< c       When iprint > 0, the file iterate.dat will be created to
< c                        summarize the iteration.
---
> c       
486c480
<       integer          i,k,nintol,itfile,iback,nskip,
---
>       integer          i,k,nintol,iback,nskip,
494a489,490
>       double precision dlamch
>       external         dlamch
499c495
<          epsmch = epsilon(one)
---
>          epsmch = 2 * dlamch('e')
547,553c543
<          info = 0
< 
<          itfile = 8
<          if (iprint .ge. 1) then
< c                                open a summary file 'iterate.dat'
<             open (8, file = 'iterate.dat', status = 'unknown')
<          endif            
---
>          info = 0          
559c549
<             call prn3lb(n,x,f,task,iprint,info,itfile,
---
>             call prn3lb(n,x,f,task,iprint,info,
566c556
<          call prn1lb(n,m,l,u,x,iprint,itfile,epsmch)
---
>          call prn1lb(n,m,l,u,x,iprint,epsmch)
583d572
<          itfile = isave(3)
649d637
<          write (itfile,1003) iter,nfgv,sbgnrm,f
684c672
<          if(iprint .ge. 1) write (6, 1005)
---
>          if(iprint .ge. 1) write (0, 1005)
732c720
<          if(iprint .ge. 1) write (6, 1006)
---
>          if(iprint .ge. 1) write (0, 1006)
758c746
<          if(iprint .ge. 1) write (6, 1005)
---
>          if(iprint .ge. 1) write (0, 1005)
789,790c777,778
<      +            boxed,cnstnd,csave,isave(22),dsave(17))
<       if (info .ne. 0 .or. iback .ge. 20) then
---
>      +            boxed,cnstnd,csave,isave(22),dsave(17), iprint)
>       if (info .ne. 0 .or. iback .ge. maxls) then
809c797
<             if(iprint .ge. 1) write (6, 1008)
---
>             if(iprint .ge. 1) write (0, 1008)
837c825
<          call prn2lb(n,x,f,g,iprint,itfile,iter,nfgv,nact,
---
>          call prn2lb(n,x,f,g,iprint,iter,nfgv,nact,
907c895
<          if(iprint .ge. 1) write (6, 1007)
---
>          if(iprint .ge. 1) write (0, 1007)
930c918
<       call prn3lb(n,x,f,task,iprint,info,itfile,
---
>       call prn3lb(n,x,f,task,iprint,info,
943,944c931
<       isave(1)  = nintol 
<       isave(3)  = itfile 
---
>       isave(1)  = nintol
1090c1077
<          if (prjctd) write (6,*)
---
>          if (prjctd) write (0,*)
1093c1080
<      +      write (6,*) 'This problem is unconstrained.'
---
>      +      write (0,*) 'This problem is unconstrained.'
1364,1365c1351
< c       When iprint > 0, the file iterate.dat will be created to
< c                        summarize the iteration.
---
> c       
2456c2442
<      +                  isave, dsave)
---
>      +                  isave, dsave, iprint)
2461c2447
<      +                 nbd(n), isave(2)
---
>      +                 nbd(n), isave(2), iprint
2472a2459,2462
> c     Be mindful that the dcsrch subroutine being called is a copy in
> c       this file (lbfgsb.f) and NOT in the Minpack2 copy distributed
> c       by scipy.
> c      
2553c2543,2545
<             write(6,*)' ascent direction in projection gd = ', gd
---
>             if (iprint .ge. 0) then
>                 write(0,*)' ascent direction in projection gd = ', gd
>             endif
2569a2562
> c        take step and prevent rounding error beyond bound
2571a2565,2566
>                if (nbd(i).eq.1.or.nbd(i).eq.2) x(i) = max(x(i), l(i))
>                if (nbd(i).eq.2.or.nbd(i).eq.3) x(i) = min(x(i), u(i))
2671c2666
<       subroutine prn1lb(n, m, l, u, x, iprint, itfile, epsmch)
---
>       subroutine prn1lb(n, m, l, u, x, iprint, epsmch)
2673c2668
<       integer n, m, iprint, itfile
---
>       integer n, m, iprint
2703,2705d2697
<             write (itfile,2001) epsmch
<             write (itfile,*)'N = ',n,'    M = ',m
<             write (itfile,9001)
2742c2734
<       subroutine prn2lb(n, x, f, g, iprint, itfile, iter, nfgv, nact, 
---
>       subroutine prn2lb(n, x, f, g, iprint, iter, nfgv, nact, 
2746c2738
<       integer          n, iprint, itfile, iter, nfgv, nact, nseg,
---
>       integer          n, iprint, iter, nfgv, nact, nseg,
2796,2797c2788
<       if (iprint .ge. 1) write (itfile,3001)
<      +          iter,nfgv,nseg,nact,word,iback,stp,xstep,sbgnrm,f
---
>       
2810c2801
<       subroutine prn3lb(n, x, f, task, iprint, info, itfile, 
---
>       subroutine prn3lb(n, x, f, task, iprint, info, 
2817c2808
<       integer          n, iprint, info, itfile, iter, nfgv, nintol,
---
>       integer          n, iprint, info, iter, nfgv, nintol,
2860,2865c2851,2856
<             if (info .eq. -1) write (6,9011)
<             if (info .eq. -2) write (6,9012)
<             if (info .eq. -3) write (6,9013)
<             if (info .eq. -4) write (6,9014)
<             if (info .eq. -5) write (6,9015)
<             if (info .eq. -6) write (6,*)' Input nbd(',k,') is invalid.'
---
>             if (info .eq. -1) write (0,9011)
>             if (info .eq. -2) write (0,9012)
>             if (info .eq. -3) write (0,9013)
>             if (info .eq. -4) write (0,9014)
>             if (info .eq. -5) write (0,9015)
>             if (info .eq. -6) write (0,*)' Input nbd(',k,') is invalid.'
2868,2888c2859,2860
<             if (info .eq. -8) write (6,9018)
<             if (info .eq. -9) write (6,9019)
<          endif
<          if (iprint .ge. 1) write (6,3007) cachyt,sbtime,lnscht
<          write (6,3008) time
<          if (iprint .ge. 1) then
<             if (info .eq. -4 .or. info .eq. -9) then
<                write (itfile,3002)
<      +             iter,nfgv,nseg,nact,word,iback,stp,xstep
<             endif
<             write (itfile,3009) task
<             if (info .ne. 0) then
<                if (info .eq. -1) write (itfile,9011)
<                if (info .eq. -2) write (itfile,9012)
<                if (info .eq. -3) write (itfile,9013)
<                if (info .eq. -4) write (itfile,9014)
<                if (info .eq. -5) write (itfile,9015)
<                if (info .eq. -8) write (itfile,9018)
<                if (info .eq. -9) write (itfile,9019)
<             endif
<             write (itfile,3008) time
---
>             if (info .eq. -8) write (0,9018)
>             if (info .eq. -9) write (0,9019)
2931,2932c2903,2905
<      +' Line search cannot locate an adequate point after 20 function',/
<      +,'  and gradient evaluations.  Previous x, f and g restored.',/,
---
>      +' Line search cannot locate an adequate point after MAXLS',/
>      +,'  function and gradient evaluations.',/
>      +,'  Previous x, f and g restored.',/,
2974a2948,2952
>         if (gi.ne.gi) then
> c          NaN value in gradient: propagate it
>            sbgnrm = gi
>            return
>         endif
3147,3148c3125
< c       When iprint > 0, the file iterate.dat will be created to
< c                        summarize the iteration.
---
> c       
3282,3283c3259,3262
<          write(6,*) ' Positive dir derivative in projection '
<          write(6,*) ' Using the backtracking step '
---
>          if (iprint .ge. 0) then
>             write(6,*) ' Positive dir derivative in projection '
>             write(6,*) ' Using the backtracking step '
>          endif
Gnimuc commented 3 years ago

Looks like this is the feature mentioned in #2

rht commented 3 years ago

OK, this issue is a duplicate then.

Gnimuc commented 3 years ago

I think we should close #2 and keep this issue open.