Limited-memory BFGS

Sometimes, the second-order derivatives are just too expensive to evaluate. In that case, it is often a good idea to approximate the Hessian matrix. The BFGS algorithm uses the first-order derivatives (gradient and tranposed-Jacobian vector product) to approximate the Hessian of the Lagrangian. LBFGS is a variant of BFGS that computes a low-rank approximation of the Hessian matrix from the past iterates. That way, LBFGS does not have to store a large dense matrix in memory, rendering the algorithm appropriate in the large-scale regime.

MadNLP implements several quasi-Newton approximation for the Hessian matrix. In this page, we focus on the limited-memory version of the BFGS algorithm, commonly known as LBFGS. We refer to this article for a detailed description of the method.

How to set up LBFGS in MadNLP?

We look at the elec optimization problem from the COPS benchmark:

function elec_model(np)
    Random.seed!(1)
    # Set the starting point to a quasi-uniform distribution of electrons on a unit sphere
    theta = (2pi) .* rand(np)
    phi = pi .* rand(np)

    core = ExaModels.ExaCore(Float64)
    x = ExaModels.variable(core, 1:np; start = [cos(theta[i])*sin(phi[i]) for i=1:np])
    y = ExaModels.variable(core, 1:np; start = [sin(theta[i])*sin(phi[i]) for i=1:np])
    z = ExaModels.variable(core, 1:np; start = [cos(phi[i]) for i=1:np])
    # Coulomb potential
    itr = [(i,j) for i in 1:np-1 for j in i+1:np]
    ExaModels.objective(core, 1.0 / sqrt((x[i] - x[j])^2 + (y[i] - y[j])^2 + (z[i] - z[j])^2) for (i,j) in itr)
    # Unit-ball
    ExaModels.constraint(core, x[i]^2 + y[i]^2 + z[i]^2 - 1 for i=1:np)

    return ExaModels.ExaModel(core)
end

elec_model (generic function with 1 method)

The problem computes the positions of the electrons in an atom. The potential of each electron depends on the positions of all the other electrons: all the variables in the problem are coupled, resulting in a dense Hessian matrix. Hence, the problem is good candidate for a quasi-Newton algorithm.

We start by solving the problem with the default options in MadNLP, using a dense linear solver from LAPACK:

nh = 10
nlp = elec_model(nh)
results_hess = madnlp(
    nlp;
    linear_solver=LapackCPUSolver,
)
nothing

This is MadNLP version v0.8.6, running with Lapack-CPU (BUNCHKAUFMAN)

Number of nonzeros in constraint Jacobian............:       30
Number of nonzeros in Lagrangian Hessian.............:      975

Total number of variables............................:       30
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:       10
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  4.8810842e+01 2.22e-16 1.71e+01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  4.8681505e+01 3.65e-06 1.69e+01  -1.0 1.69e-03   4.0 1.00e+00 1.00e+00h  1
   2  4.8316232e+01 3.04e-05 1.62e+01  -1.0 4.87e-03   3.5 1.00e+00 1.00e+00h  1
   3  4.7384770e+01 2.21e-04 1.47e+01  -1.0 1.31e-02   3.0 1.00e+00 1.00e+00h  1
   4  4.5463791e+01 1.23e-03 1.18e+01  -1.0 3.06e-02   2.6 1.00e+00 1.00e+00h  1
   5  4.2506911e+01 4.51e-03 7.99e+00  -1.0 5.67e-02   2.1 1.00e+00 1.00e+00h  1
   6  3.7327650e+01 3.00e-01 4.13e+00  -1.0 4.12e-01   1.6 1.00e+00 1.00e+00h  1
   7  3.4989135e+01 4.08e-02 2.26e+00  -1.0 1.67e-01   1.1 1.00e+00 1.00e+00h  1
   8  3.3060619e+01 5.97e-02 8.75e-01  -1.0 1.85e-01   0.7 1.00e+00 1.00e+00h  1
   9  2.2921574e+01 7.16e+00 6.98e+00  -1.7 2.14e+00   0.2 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  2.6930954e+01 1.70e+00 1.43e+00  -1.7 1.15e+00    -  1.00e+00 1.00e+00h  1
  11  3.1634663e+01 2.92e-01 9.84e-01  -1.7 4.65e-01   0.6 1.00e+00 1.00e+00h  1
  12  3.1839102e+01 2.77e-01 5.39e-01  -1.7 4.30e-01   0.1 1.00e+00 1.00e+00h  1
  13  3.2320199e+01 7.38e-02 2.82e-01  -1.7 2.54e-01    -  1.00e+00 1.00e+00h  1
  14  3.2674733e+01 9.98e-03 5.99e-02  -1.7 9.29e-02  -0.3 1.00e+00 1.00e+00h  1
  15  3.2684885e+01 5.24e-03 1.60e-02  -2.5 5.78e-02  -0.8 1.00e+00 1.00e+00h  1
  16  3.2560375e+01 1.69e-02 6.92e-02  -3.8 1.11e-01    -  1.00e+00 1.00e+00h  1
  17  3.2706370e+01 1.35e-03 4.25e-03  -3.8 2.48e-02    -  1.00e+00 1.00e+00h  1
  18  3.2707565e+01 9.97e-04 3.60e-03  -3.8 2.80e-02    -  1.00e+00 1.00e+00h  1
  19  3.2716885e+01 6.42e-06 2.14e-05  -3.8 2.43e-03    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  3.2716949e+01 1.19e-08 4.33e-08  -5.7 1.01e-04    -  1.00e+00 1.00e+00h  1
  21  3.2716949e+01 4.44e-16 2.22e-15  -8.6 1.81e-08    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 21

                                   (scaled)                 (unscaled)
Objective...............:   3.2716949460147582e+01    3.2716949460147582e+01
Dual infeasibility......:   2.2204460492503131e-15    2.2204460492503131e-15
Constraint violation....:   4.4408920985006262e-16    4.4408920985006262e-16
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   2.2204460492503131e-15    2.2204460492503131e-15

Number of objective function evaluations             = 22
Number of objective gradient evaluations             = 22
Number of constraint evaluations                     = 22
Number of constraint Jacobian evaluations            = 22
Number of Lagrangian Hessian evaluations             = 21
Total wall-clock secs in solver (w/o fun. eval./lin. alg.)  =  4.322
Total wall-clock secs in linear solver                      =  0.001
Total wall-clock secs in NLP function evaluations           =  0.000
Total wall-clock secs                                       =  4.322

EXIT: Optimal Solution Found (tol = 1.0e-08).

We observe that MadNLP converges in 21 iterations.

To replace the second-order derivatives by an LBFGS approximation, you should pass the option hessian_approximation=CompactLBFGS to MadNLP.

results_qn = madnlp(
    nlp;
    linear_solver=LapackCPUSolver,
    hessian_approximation=MadNLP.CompactLBFGS,
)
nothing

This is MadNLP version v0.8.6, running with Lapack-CPU (BUNCHKAUFMAN)

Number of nonzeros in constraint Jacobian............:       30
Number of nonzeros in Lagrangian Hessian.............:      975

Total number of variables............................:       30
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:       10
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  4.8810842e+01 2.22e-16 1.71e+01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  2.2956342e+00 1.50e+03 2.04e+02  -1.0 2.74e+02    -  1.00e+00 1.25e-01f  4
   2  4.4967927e+00 3.76e+02 5.48e+01  -1.0 1.71e+01    -  1.00e+00 1.00e+00h  1
   3  8.7093834e+00 9.38e+01 1.37e+01  -1.0 8.37e+00    -  1.00e+00 1.00e+00h  1
   4  1.4770500e+01 2.32e+01 3.20e+00  -1.0 4.18e+00    -  1.00e+00 1.00e+00h  1
   5  1.8176499e+01 9.94e+00 1.50e+00  -1.0 3.09e+00    -  1.00e+00 1.00e+00h  1
   6  2.5109575e+01 2.32e+00 1.92e+00  -1.0 1.28e+00    -  1.00e+00 1.00e+00h  1
   7  2.8512602e+01 2.52e+00 2.52e+00  -1.0 1.14e+00    -  1.00e+00 1.00e+00h  1
   8  1.2760325e+01 3.54e+01 1.13e+01  -1.0 4.99e+00    -  1.00e+00 1.00e+00h  1
   9  1.5623867e+01 1.58e+01 2.83e+00  -1.0 3.43e+00    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  2.1982619e+01 5.68e+00 1.87e+00  -1.0 2.04e+00    -  1.00e+00 1.00e+00h  1
  11  2.2550708e+01 4.14e+00 1.15e+00  -1.0 2.61e+00    -  1.00e+00 5.00e-01h  2
  12  2.6617499e+01 1.76e+00 1.64e+00  -1.0 1.19e+00    -  1.00e+00 1.00e+00h  1
  13  2.8301250e+01 1.78e+00 2.36e+00  -1.0 9.16e-01    -  1.00e+00 1.00e+00h  1
  14  2.9440900e+01 8.43e-01 1.73e+00  -1.0 8.69e-01    -  1.00e+00 1.00e+00h  1
  15  3.0789962e+01 4.26e-01 1.73e+00  -1.0 6.18e-01    -  1.00e+00 1.00e+00h  1
  16  3.2172088e+01 1.77e-01 1.13e+00  -1.0 3.85e-01    -  1.00e+00 1.00e+00h  1
  17  3.2605615e+01 3.60e-02 3.88e-01  -1.0 1.55e-01    -  1.00e+00 1.00e+00h  1
  18  3.2790290e+01 9.45e-03 2.48e-01  -1.7 7.30e-02    -  1.00e+00 1.00e+00h  1
  19  3.2777796e+01 6.22e-03 1.87e-01  -1.7 7.70e-02    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  3.2204263e+01 1.44e-01 3.75e-01  -2.5 3.76e-01    -  1.00e+00 1.00e+00h  1
  21  3.2584660e+01 3.57e-02 2.45e-01  -2.5 1.74e-01    -  1.00e+00 1.00e+00h  1
  22  3.2700371e+01 5.30e-03 9.80e-02  -2.5 6.16e-02    -  1.00e+00 1.00e+00h  1
  23  3.2715038e+01 8.46e-04 4.00e-02  -2.5 2.51e-02    -  1.00e+00 1.00e+00h  1
  24  3.2717788e+01 4.35e-04 3.50e-02  -2.5 2.00e-02    -  1.00e+00 1.00e+00h  1
  25  3.2720180e+01 1.32e-04 2.08e-02  -2.5 9.20e-03    -  1.00e+00 1.00e+00h  1
  26  3.2720685e+01 3.56e-05 4.95e-03  -3.8 4.31e-03    -  1.00e+00 1.00e+00h  1
  27  3.2720835e+01 3.10e-06 2.66e-03  -3.8 1.57e-03    -  1.00e+00 1.00e+00h  1
  28  3.2720457e+01 4.49e-05 5.46e-03  -3.8 4.73e-03    -  1.00e+00 1.00e+00h  1
  29  3.2714888e+01 6.03e-04 2.59e-02  -3.8 2.23e-02    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  3.1599309e+01 1.26e-01 5.67e-01  -3.8 3.34e-01    -  1.00e+00 1.00e+00h  1
  31  3.1274890e+01 1.40e-01 3.39e-01  -3.8 3.65e-01    -  1.00e+00 1.00e+00h  1
  32  3.2666066e+01 1.48e-02 5.14e-01  -3.8 9.31e-02    -  1.00e+00 1.00e+00h  1
  33  3.2638320e+01 1.65e-02 1.35e-01  -3.8 9.20e-02    -  1.00e+00 1.00e+00h  1
  34  3.2704002e+01 5.69e-03 1.34e-01  -3.8 5.75e-02    -  1.00e+00 1.00e+00h  1
  35  3.2714347e+01 1.20e-03 1.56e-02  -3.8 3.01e-02    -  1.00e+00 1.00e+00h  1
  36  3.2719760e+01 7.06e-05 6.75e-03  -3.8 7.95e-03    -  1.00e+00 1.00e+00h  1
  37  3.2713492e+01 1.02e-03 3.18e-02  -3.8 2.40e-02    -  1.00e+00 1.00e+00h  1
  38  3.2700419e+01 2.46e-03 6.71e-02  -3.8 3.78e-02    -  1.00e+00 1.00e+00h  1
  39  3.2704367e+01 1.58e-03 4.67e-02  -3.8 4.79e-02    -  1.00e+00 5.00e-01h  2
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  40  3.2710343e+01 9.34e-04 2.87e-02  -3.8 2.91e-02    -  1.00e+00 1.00e+00h  1
  41  3.2708148e+01 1.46e-03 3.30e-02  -3.8 3.07e-02    -  1.00e+00 1.00e+00h  1
  42  3.2661678e+01 5.95e-03 5.50e-02  -3.8 6.67e-02    -  1.00e+00 1.00e+00h  1
  43  3.2678493e+01 5.20e-03 2.38e-02  -3.8 5.86e-02    -  1.00e+00 1.00e+00h  1
  44  3.2696745e+01 2.70e-03 3.89e-02  -3.8 3.33e-02    -  1.00e+00 5.00e-01h  2
  45  3.2715793e+01 1.98e-04 6.86e-03  -3.8 1.19e-02    -  1.00e+00 1.00e+00h  1
  46  3.2716827e+01 1.19e-05 2.56e-03  -3.8 3.30e-03    -  1.00e+00 1.00e+00h  1
  47  3.2711738e+01 6.09e-04 1.98e-02  -3.8 1.87e-02    -  1.00e+00 1.00e+00h  1
  48  3.2716512e+01 1.02e-04 1.95e-02  -3.8 7.33e-03    -  1.00e+00 1.00e+00h  1
  49  3.2715191e+01 3.29e-04 4.43e-03  -3.8 1.38e-02    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  50  3.2716812e+01 3.78e-05 1.33e-02  -3.8 5.06e-03    -  1.00e+00 1.00e+00h  1
  51  3.2716764e+01 3.58e-05 9.62e-04  -3.8 4.77e-03    -  1.00e+00 1.00e+00h  1
  52  3.2716948e+01 7.95e-07 1.63e-03  -5.7 7.15e-04    -  1.00e+00 1.00e+00h  1
  53  3.2716950e+01 1.15e-07 3.10e-04  -5.7 3.28e-04    -  1.00e+00 1.00e+00h  1
  54  3.2716951e+01 9.76e-09 2.13e-04  -5.7 9.30e-05    -  1.00e+00 1.00e+00h  1
  55  3.2716936e+01 1.60e-06 7.83e-04  -5.7 9.59e-04    -  1.00e+00 1.00e+00h  1
  56  3.2716869e+01 1.23e-05 5.09e-03  -5.7 2.77e-03    -  1.00e+00 1.00e+00h  1
  57  3.2716912e+01 6.32e-06 1.35e-03  -5.7 1.80e-03    -  1.00e+00 1.00e+00h  1
  58  3.2716927e+01 3.76e-06 4.73e-04  -5.7 1.38e-03    -  1.00e+00 5.00e-01h  2
  59  3.2716949e+01 1.53e-07 5.53e-04  -5.7 2.78e-04    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  60  3.2716949e+01 3.55e-08 6.97e-05  -5.7 1.80e-04    -  1.00e+00 1.00e+00h  1
  61  3.2716949e+01 1.64e-09 7.05e-05  -5.7 3.54e-05    -  1.00e+00 1.00e+00h  1
  62  3.2716949e+01 7.95e-10 4.53e-05  -5.7 2.43e-05    -  1.00e+00 1.00e+00h  1
  63  3.2716949e+01 5.08e-08 2.33e-04  -5.7 1.96e-04    -  1.00e+00 1.00e+00h  1
  64  3.2716948e+01 1.89e-07 5.26e-04  -5.7 3.13e-04    -  1.00e+00 1.00e+00h  1
  65  3.2716949e+01 1.58e-07 1.74e-04  -5.7 3.21e-04    -  1.00e+00 1.00e+00h  1
  66  3.2716949e+01 8.10e-08 2.11e-04  -5.7 1.92e-04    -  1.00e+00 5.00e-01h  2
  67  3.2716949e+01 8.24e-09 1.56e-05  -5.7 7.33e-05    -  1.00e+00 1.00e+00h  1
  68  3.2716949e+01 2.11e-10 2.76e-05  -8.6 1.15e-05    -  1.00e+00 1.00e+00h  1
  69  3.2716949e+01 8.34e-11 3.36e-06  -8.6 7.01e-06    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  70  3.2716949e+01 2.09e-12 2.47e-06  -8.6 1.04e-06    -  1.00e+00 1.00e+00h  1
  71  3.2716949e+01 2.47e-10 1.04e-05  -8.6 1.37e-05    -  1.00e+00 1.00e+00h  1
  72  3.2716949e+01 7.76e-10 2.84e-05  -8.6 2.10e-05    -  1.00e+00 1.00e+00h  1
  73  3.2716949e+01 3.25e-10 2.58e-05  -8.6 1.44e-05    -  1.00e+00 1.00e+00h  1
  74  3.2716949e+01 3.48e-10 9.13e-06  -8.6 1.43e-05    -  1.00e+00 1.00e+00h  1
  75  3.2716949e+01 4.80e-11 1.09e-05  -8.6 6.67e-06    -  1.00e+00 1.00e+00h  1
  76  3.2716949e+01 1.21e-11 1.31e-06  -8.6 3.13e-06    -  1.00e+00 1.00e+00h  1
  77  3.2716949e+01 5.24e-13 4.86e-07  -8.6 5.30e-07    -  1.00e+00 1.00e+00h  1
  78  3.2716949e+01 1.01e-13 3.66e-07  -8.6 2.57e-07    -  1.00e+00 1.00e+00h  1
  79  3.2716949e+01 1.06e-13 4.60e-07  -8.6 2.46e-07    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  80  3.2716949e+01 2.91e-14 2.56e-07  -8.6 1.49e-07    -  1.00e+00 1.00e+00h  1
  81  3.2716949e+01 1.98e-14 1.24e-07  -8.6 1.02e-07    -  1.00e+00 1.00e+00h  1
  82  3.2716949e+01 2.33e-13 2.21e-07  -8.6 3.75e-07    -  1.00e+00 1.00e+00h  1
  83  3.2716949e+01 2.96e-12 1.12e-06  -8.6 1.26e-06    -  1.00e+00 1.00e+00h  1
  84  3.2716949e+01 1.54e-12 1.84e-06  -8.6 8.76e-07    -  1.00e+00 1.00e+00h  1
  85  3.2716949e+01 2.68e-12 6.57e-07  -8.6 1.42e-06    -  1.00e+00 1.00e+00h  1
  86  3.2716949e+01 9.15e-13 1.45e-06  -8.6 8.49e-07    -  1.00e+00 1.00e+00h  1
  87  3.2716949e+01 3.18e-13 2.79e-07  -8.6 4.07e-07    -  1.00e+00 1.00e+00h  1
  88  3.2716949e+01 8.08e-14 6.82e-08  -8.6 2.11e-07    -  1.00e+00 1.00e+00h  1
  89  3.2716949e+01 4.31e-14 4.28e-07  -8.6 1.76e-07    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  90  3.2716949e+01 3.26e-14 2.05e-08  -8.6 1.58e-07    -  1.00e+00 1.00e+00h  1
  91  3.2716949e+01 4.44e-16 3.23e-08  -9.0 1.67e-08    -  1.00e+00 1.00e+00h  1
  92  3.2716949e+01 2.22e-16 5.88e-08  -9.0 2.20e-08    -  1.00e+00 1.00e+00H  1
  93  3.2716949e+01 4.44e-16 9.14e-09  -9.0 1.43e-08    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 93

                                   (scaled)                 (unscaled)
Objective...............:   3.2716949460147596e+01    3.2716949460147596e+01
Dual infeasibility......:   9.1355145115379344e-09    9.1355145115379344e-09
Constraint violation....:   4.4408920985006262e-16    4.4408920985006262e-16
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   9.1355145115379344e-09    9.1355145115379344e-09

Number of objective function evaluations             = 106
Number of objective gradient evaluations             = 94
Number of constraint evaluations                     = 106
Number of constraint Jacobian evaluations            = 94
Number of Lagrangian Hessian evaluations             = 0
Total wall-clock secs in solver (w/o fun. eval./lin. alg.)  =  6.293
Total wall-clock secs in linear solver                      =  0.006
Total wall-clock secs in NLP function evaluations           =  0.000
Total wall-clock secs                                       =  6.298

EXIT: Optimal Solution Found (tol = 1.0e-08).

We observe that MadNLP converges in 93 iterations. As expected, the number of Hessian evaluations is 0.

How to tune the options in LBGS?

You can tune the LBFGS options by using the option quasi_newton_options. The option takes as input a QuasiNewtonOptions structure, with the following attributes (we put on the right their default values):

init_strategy::BFGSInitStrategy = SCALAR1
max_history::Int = 6
init_value::Float64 = 1.0
sigma_min::Float64 = 1e-8
sigma_max::Float64 = 1e+8

The most important parameter is max_history, which encodes the number of vectors used in the low-rank approximation. For instance, we can increase the history to use the 20 past iterates using:

qn_options = MadNLP.QuasiNewtonOptions(; max_history=20)
results_qn = madnlp(
    nlp;
    linear_solver=LapackCPUSolver,
    hessian_approximation=MadNLP.CompactLBFGS,
    quasi_newton_options=qn_options,
)
nothing

This is MadNLP version v0.8.6, running with Lapack-CPU (BUNCHKAUFMAN)

Number of nonzeros in constraint Jacobian............:       30
Number of nonzeros in Lagrangian Hessian.............:      975

Total number of variables............................:       30
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:       10
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  4.8810842e+01 2.22e-16 1.71e+01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  2.2956342e+00 1.50e+03 2.04e+02  -1.0 2.74e+02    -  1.00e+00 1.25e-01f  4
   2  4.4967927e+00 3.76e+02 5.48e+01  -1.0 1.71e+01    -  1.00e+00 1.00e+00h  1
   3  8.7093834e+00 9.38e+01 1.37e+01  -1.0 8.37e+00    -  1.00e+00 1.00e+00h  1
   4  1.4770500e+01 2.32e+01 3.20e+00  -1.0 4.18e+00    -  1.00e+00 1.00e+00h  1
   5  1.8176499e+01 9.94e+00 1.50e+00  -1.0 3.09e+00    -  1.00e+00 1.00e+00h  1
   6  2.5109575e+01 2.32e+00 1.92e+00  -1.0 1.28e+00    -  1.00e+00 1.00e+00h  1
   7  2.8512602e+01 2.52e+00 2.52e+00  -1.0 1.14e+00    -  1.00e+00 1.00e+00h  1
   8  1.2759872e+01 3.54e+01 1.13e+01  -1.0 4.99e+00    -  1.00e+00 1.00e+00h  1
   9  1.5624022e+01 1.58e+01 2.83e+00  -1.0 3.43e+00    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  2.1982310e+01 5.68e+00 1.86e+00  -1.0 2.04e+00    -  1.00e+00 1.00e+00h  1
  11  2.2542843e+01 4.22e+00 1.14e+00  -1.0 2.56e+00    -  1.00e+00 5.00e-01h  2
  12  2.7255011e+01 1.47e+00 2.96e+00  -1.0 1.07e+00    -  1.00e+00 1.00e+00h  1
  13  2.9232894e+01 1.25e+00 2.66e+00  -1.0 8.60e-01    -  1.00e+00 1.00e+00h  1
  14  2.9261067e+01 9.81e-01 3.03e+00  -1.0 1.07e+00    -  1.00e+00 5.00e-01h  2
  15  2.7142323e+01 2.39e+00 3.27e+00  -1.0 1.42e+00    -  1.00e+00 1.00e+00h  1
  16  2.7776481e+01 1.35e+00 1.98e+00  -1.0 1.02e+00    -  1.00e+00 1.00e+00h  1
  17  3.0166747e+01 7.86e-01 1.80e+00  -1.0 6.79e-01    -  1.00e+00 1.00e+00h  1
  18  3.0609464e+01 3.09e-01 8.42e-01  -1.0 4.90e-01    -  1.00e+00 1.00e+00h  1
  19  3.2450081e+01 6.00e-02 5.04e-01  -1.7 2.27e-01    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  3.2636130e+01 3.88e-02 3.94e-01  -1.7 1.39e-01    -  1.00e+00 1.00e+00h  1
  21  3.2706987e+01 9.35e-03 1.66e-01  -1.7 6.88e-02    -  1.00e+00 1.00e+00h  1
  22  3.2699721e+01 6.51e-03 9.41e-02  -2.5 6.69e-02    -  1.00e+00 1.00e+00h  1
  23  3.2682575e+01 5.11e-03 1.14e-01  -2.5 6.73e-02    -  1.00e+00 1.00e+00h  1
  24  3.2701207e+01 4.15e-03 1.48e-01  -2.5 4.17e-02    -  1.00e+00 1.00e+00h  1
  25  3.2712954e+01 1.75e-03 4.88e-02  -2.5 3.79e-02    -  1.00e+00 1.00e+00h  1
  26  3.2719381e+01 2.29e-04 2.32e-02  -2.5 1.40e-02    -  1.00e+00 1.00e+00h  1
  27  3.2720051e+01 3.26e-05 6.65e-03  -3.8 4.78e-03    -  1.00e+00 1.00e+00h  1
  28  3.2720061e+01 2.64e-05 7.96e-03  -3.8 4.67e-03    -  1.00e+00 1.00e+00h  1
  29  3.2715805e+01 6.76e-04 2.85e-02  -3.8 2.41e-02    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  3.2702239e+01 2.73e-03 5.82e-02  -3.8 4.09e-02    -  1.00e+00 1.00e+00h  1
  31  3.2697960e+01 2.28e-03 7.47e-02  -3.8 3.88e-02    -  1.00e+00 1.00e+00h  1
  32  3.2707411e+01 2.46e-03 1.24e-01  -3.8 3.57e-02    -  1.00e+00 1.00e+00h  1
  33  3.2709995e+01 9.33e-04 1.60e-02  -3.8 2.89e-02    -  1.00e+00 1.00e+00h  1
  34  3.2712746e+01 6.45e-04 2.11e-02  -3.8 2.08e-02    -  1.00e+00 1.00e+00h  1
  35  3.2704283e+01 3.24e-03 5.49e-02  -3.8 4.93e-02    -  1.00e+00 1.00e+00h  1
  36  3.2709923e+01 1.48e-03 2.22e-02  -3.8 3.35e-02    -  1.00e+00 1.00e+00h  1
  37  3.2715749e+01 4.96e-04 5.36e-02  -3.8 1.79e-02    -  1.00e+00 1.00e+00h  1
  38  3.2714941e+01 4.33e-04 2.54e-03  -3.8 1.74e-02    -  1.00e+00 1.00e+00h  1
  39  3.2716951e+01 2.19e-06 1.15e-03  -3.8 1.21e-03    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  40  3.2716838e+01 3.01e-05 4.88e-03  -5.7 4.19e-03    -  1.00e+00 1.00e+00h  1
  41  3.2716805e+01 2.19e-05 3.16e-03  -5.7 4.44e-03    -  1.00e+00 1.00e+00h  1
  42  3.2716579e+01 5.15e-05 1.04e-02  -5.7 6.11e-03    -  1.00e+00 1.00e+00h  1
  43  3.2716654e+01 3.23e-05 1.67e-03  -5.7 5.14e-03    -  1.00e+00 1.00e+00h  1
  44  3.2716894e+01 1.17e-05 3.18e-03  -5.7 2.93e-03    -  1.00e+00 1.00e+00h  1
  45  3.2716900e+01 7.34e-06 1.19e-03  -5.7 2.14e-03    -  1.00e+00 1.00e+00h  1
  46  3.2716935e+01 2.79e-06 2.69e-03  -5.7 1.54e-03    -  1.00e+00 1.00e+00h  1
  47  3.2716941e+01 1.70e-06 1.41e-04  -5.7 1.14e-03    -  1.00e+00 1.00e+00h  1
  48  3.2716949e+01 4.56e-08 2.20e-04  -5.7 1.72e-04    -  1.00e+00 1.00e+00h  1
  49  3.2716949e+01 6.08e-08 1.21e-04  -5.7 2.04e-04    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  50  3.2716949e+01 1.17e-07 4.18e-04  -5.7 2.59e-04    -  1.00e+00 1.00e+00h  1
  51  3.2716949e+01 8.97e-08 1.03e-04  -5.7 2.36e-04    -  1.00e+00 1.00e+00h  1
  52  3.2716949e+01 3.50e-08 1.97e-04  -5.7 1.50e-04    -  1.00e+00 1.00e+00h  1
  53  3.2716949e+01 1.82e-08 4.78e-05  -5.7 1.07e-04    -  1.00e+00 1.00e+00h  1
  54  3.2716949e+01 2.88e-09 6.40e-05  -5.7 4.17e-05    -  1.00e+00 1.00e+00h  1
  55  3.2716949e+01 9.33e-10 8.16e-06  -5.7 2.42e-05    -  1.00e+00 1.00e+00h  1
  56  3.2716949e+01 5.70e-11 1.28e-05  -8.6 6.71e-06    -  1.00e+00 1.00e+00h  1
  57  3.2716949e+01 2.92e-11 5.40e-07  -8.6 4.75e-06    -  1.00e+00 1.00e+00h  1
  58  3.2716949e+01 8.10e-14 3.94e-07  -8.6 2.74e-07    -  1.00e+00 1.00e+00h  1
  59  3.2716949e+01 1.82e-14 1.96e-08  -8.6 1.28e-07    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  60  3.2716949e+01 2.22e-16 7.68e-09  -9.0 9.04e-09    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 60

                                   (scaled)                 (unscaled)
Objective...............:   3.2716949460147582e+01    3.2716949460147582e+01
Dual infeasibility......:   7.6805297677395856e-09    7.6805297677395856e-09
Constraint violation....:   2.2204460492503131e-16    2.2204460492503131e-16
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   7.6805297677395856e-09    7.6805297677395856e-09

Number of objective function evaluations             = 69
Number of objective gradient evaluations             = 61
Number of constraint evaluations                     = 69
Number of constraint Jacobian evaluations            = 61
Number of Lagrangian Hessian evaluations             = 0
Total wall-clock secs in solver (w/o fun. eval./lin. alg.)  =  0.005
Total wall-clock secs in linear solver                      =  0.009
Total wall-clock secs in NLP function evaluations           =  0.000
Total wall-clock secs                                       =  0.014

EXIT: Optimal Solution Found (tol = 1.0e-08).

We observe that the total number of iterations has been reduced from 93 to 60.

How does LBFGS is implemented in MadNLP?

MadNLP implements the compact LBFGS algorithm described in this article. At each iteration, the Hessian $W_k$ is approximated by a low rank positive definite matrix $B_k$, defined as

\[B_k = \xi_k I + U_k V_k^\top \]

with $\xi > 0$ a scaling factor, $U_k$ and $V_k$ two $n \times 2p$ matrices. The number $p$ denotes the number of vectors used when computing the limited memory updates (the parameter max_history in MadNLP): the larger, the more accurate is the low-rank approximation.

Replacing the Hessian of the Lagrangian $W_k$ by the low-rank matrix $B_k$, the KKT system solved in MadNLP rewrites as

\[\begin{bmatrix} \xi I + U V^\top + \Sigma_x & 0 & A^\top \\ 0 & \Sigma_s & -I \\ A & -I & 0 \end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta s \\ \Delta y \end{bmatrix} = \begin{bmatrix} r_1 \\ r_2 \\ r_3 \end{bmatrix} \]

The KKT system has a low-rank structure we can exploit using the Sherman-Morrison formula. The method is detailed e.g. in Section 3.8 of that paper.

Info

As MadNLP is designed to solve generic constrained optimization problems, it does not approximate the inverse of the Hessian matrix, as done in other implementations of the LBFGS algorithm specialized on solving nonlinear problems with bound constraints. If your problem has no generic nonlinear constraints, we recommend for optimal performance using LBFGSB or the LBFGS implemented in JSOSolvers.jl.