Hessian approximation

[edljk]

Is it possible to use DCISolver.jl if only analytical cost gradient and Jacobian of the constraints are available in the model?
Is there a way to use for instance BFGS / L-BFGS approximation?
Thanks for the nice work!

[abelsiqueira]

Hi @edljk, it is possible, but it is not immediately available.

It is possible to wrap a model in a way that the Hessian is substituted by a L-BFGS approximation, but it is only available as matrix-vector product, because we don't form the matrix in the L-BFGS case. A user recently asked about this for the Percival solver, in which case it is possible, and here is an example: https://gist.github.com/abelsiqueira/c4e4c590fb7f0d755f756e21e11e40ad
I plan to write a tutorial for that once I'm back from vacation.

The issue is that DCISolver.jl uses factorizations, so the strategy above doesn't work.

What is possible is to create a similar NLPModel wrapper that provides a dense BFGS - or maybe just a simple diagonal matrix - when the Hessian is asked for.
This can be implemented the same way that the LBFGSModel is, but overriding different hess methods. This wouldn't be cheap, since you will be building dense matrices, so maybe using just a diagonal instead of the BFGS approximation might be desirable.

So I see three possibilities:

• Use Percival.jl instead, with the link that I mentioned
• Implement a SimpleQPTransform model wrapper that returns the diagonal matrix instead of the Hessian.
• Implement a DenseBFGSModel model wrapper that returns the BFGS approximation of the matrix instead of the Hessian.

Let me know if this is what you are looking for, and maybe your use case, please.

[edljk]

Hi @abelsiqueira, thank you so much for the very precise answer.
Percival.jl would definitely be an option for me: I am able to provide the gradient of the cost function and the Jacobian of the equality constraints of my problem.
I tried to adapt your example without using ADNLPModels but run into troubles.. I guess it could be related to the constraint.
If you have any suggestion? Btw I do not understand what exactly does the call push!(nlp, xt, gt)

using NLPModels, ADNLPModels, Percival, NLPModelsModifiers, ManualNLPModels

A = [1.0 1.0]
b = [1.0]

f(x) = (x[1] - 1)^2 + 4 * (x[2] - x[1]^2)^2
g!(gx, x) = begin
  gx[1] = 2 * (x[1] - 1) - 16 * x[1] * (x[2] - x[1]^2)
  gx[2] = 8 * (x[2] - x[1]^2)
  gx
end

nlp = NLPModel(
  [-1.2; 1.0],
  f,
  grad = g!,
  cons = (x -> A * x .- b, [0.0], [0.0])
  )

# Wrap into LBFGS Model
qnlp = LBFGSModel(nlp)

xt = zeros(2)
gt = zeros(2)

cb = (nlp, solver, stats) -> begin
        @show stats.iter solver.x
        # Handle LBFGS update via callback
        if stats.iter > 1
                xt .-= solver.x
                xt .*= -1  # = xₖ₊₁ - xₖ
                gt .-= solver.gx
                gt .*= -1  # = ∇f(xₖ₊₁) - ∇f(xₖ)
                push!(nlp, xt, gt)
        end
        xt .= solver.x
        gt .= solver.gx
end
stats = percival(qnlp, callback=cb, verbose=1)

@show hess_op(qnlp, stats.solution) |> Matrix

[abelsiqueira]

Hi @edljk, the changes that you need are:

  cons = ((cx, x) -> cx .= A * x .- b, [0.0], [0.0]),
  jprod = (jv, x, v) -> jv .= A * v,
  jtprod = (jtv, x, v) -> jtv .= A' * v,

The cons first argument is a function (cx, x) -> ..., like the grad function.
The jprod and jtprod functions are what Percival uses instead of the full Jacobian: Jacobian- and Jacobian transposed-vector products.

Most things inside the callback are related to the L-BFGS update. The update of the L-BFGS approximation uses the vectors $s_k = x_{k+1} - x_k$ and $y_k = \nabla f(x_{k+1}) - \nabla f(x_k)$.
The callback function uses a vector xt as both the temporary storage for $x_k$, and for the temporary storage for $s_k$. Similarly for gt.

At the moment that the function push! is called, the values of xt and gt are $s_k$ and $y_k$. The push! function tells the L-BFGS model to use these two vectors for the update.

Hope this works now.

[edljk]

It works perfectly.
Thanks a lot for your help and all your work!