Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 98.48%. Comparing base (b9f325f) to head (df1bb37).
⚠️ Report is 20 commits behind head on main.

@@            Coverage Diff             @@
##             main     #275      +/-   ##
==========================================
+ Coverage   98.45%   98.48%   +0.03%     
==========================================
  Files          28       28              
  Lines        4659     4764     +105     
==========================================
+ Hits         4587     4692     +105     
  Misses         72       72              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

Benchmark Results (Julia v1)

Hi @franckgaga,
@amontoison alerted me to this PR, and I have a couple of questions/suggestions:

• If the bulk of the time is spent in actual autodiff, you could try GreedyColoringAlgorithm(; decompression=:substitution, postprocessing=true) to see if it reduces the number of colors (aka the number of autodiff passes, up to ForwardDiff chunking).
• If the bulk of the time is spent in decompression, it might be because your function fill_lowertriangle! is not optimized for SparseMatrixCSC (individual indexing is very slow on such matrices). I suggest you keep track of the issue I opened about this: Decompress into a single triangle gdalle/SparseMatrixColorings.jl#280

Many thanks for the suggestions @gdalle. I tried to differentiate actual autodiff to decompression using @profview, but I was a bit lost in the call stack. Here it is in attached files, in case you can decode it better than me, see "thread 1 (interactive)": profile.html

If the bulk of the time is spent in actual autodiff, you could try GreedyColoringAlgorithm(; decompression=:substitution, postprocessing=true) to see if it reduces the number of colors (aka the number of autodiff passes, up to ForwardDiff chunking).

I benchmarked my case study with theses two options and it essentially changes nothing on the results, compared to GreedyColoringAlgorithm().

If the bulk of the time is spent in decompression, it might be because your function fill_lowertriangle! is not optimized for SparseMatrixCSC (individual indexing is very slow on such matrices). I suggest you keep track of the issue I opened about this: gdalle/SparseMatrixColorings.jl#280

Alright thanks, I'm following it. The JuMP.@operator docstring says that I must fill in-place only the non-zero lower triangular entries only. Do you have any suggestions for a more efficient implementation of fill_lowertriangle! function?

I tried to differentiate actual autodiff to decompression using @profview, but I was a bit lost in the call stack. Here it is in attached files, in case you can decode it better than me, see "thread 1 (interactive)": profile.html

This doesn't look too bad, it's all blue and you're spending most of your time in the actual differentiation (gradients and hessians).

Do you have any suggestions for a more efficient implementation of fill_lowertriangle! function?

Here's the better implementation, which can be a thousand times faster (although it doesn't seem to impact your profiling). That's similar to what I'd like to do in SMC but I'm still figuring out API design for only filling one triangle.

using LinearAlgebra, SparseArrays, StableRNGs

function copyto_lowertriangle_naive!(T::SparseMatrixCSC, A::SparseMatrixCSC)
    for j in axes(A, 2)
        for i in axes(A, 1)
            if i >= j
                T[i, j] = A[i, j]
            end
        end
    end
    return
end

function copyto_lowertriangle!(T::SparseMatrixCSC, A::SparseMatrixCSC)
    @assert size(T) == size(A)
    kT = 0
    rvA, rvT = rowvals(A), rowvals(T)
    nzA, nzT = nonzeros(A), nonzeros(T)
    for j in axes(A, 2)
        for kA in nzrange(A, j)
            i = rvA[kA]
            if i >= j
                kT += 1
                @assert i == rvT[kT]
                nzT[kT] = nzA[kA]
            end
        end
    end
    return T
end

Benchmark and test:

julia> n, p = 1_000, 0.005;

julia> A = sparse(Symmetric(sprand(n, n, p)));

julia> T = similar(sparse(LowerTriangular(A)));

julia> copyto_lowertriangle!(T, A)

julia> T == LowerTriangular(A)
true

julia> using Chairmarks

julia> @b (T, A) copyto_lowertriangle_naive!(_[1], _[2])
7.881 ms

julia> @b (T, A) copyto_lowertriangle!(_[1], _[2])
12.958 μs

Impressive gains!

I just understood that MathOptInterface will call the in-place ∇²f function of @operator with the first argument as H::MathOptInterface.Nonlinear.ReverseAD._UnsafeLowerTriangularMatrixView, instead of a Matrix or SparseMatrixCSC. Its docstring is:

  _UnsafeLowerTriangularMatrixView(x, N)

  Lightweight unsafe view that converts a vector x into the lower-triangular component of a symmetric N-by-N matrix.

  Motivation
  ==========

  _UnsafeLowerTriangularMatrixView is needed as an allocation-free equivalent of view. Other alternatives, like reshape(view(x, 1:N^2), N, N) or a struct like

  struct _SafeView{T}
      x::Vector{T}
      len::Int
  end

  will allocate so that x can be tracked by Julia's GC. _UnsafeLowerTriangularMatrixView relies on the fact that the use-cases of _UnsafeLowerTriangularMatrixView only temporarily wrap a
  long-lived vector like d.jac_storage so that we don't have to worry about the GC removing d.jac_storage while _UnsafeLowerTriangularMatrixView exists. This lets us use a Ptr{T} and
  create a struct that is isbitstype and therefore does not allocate.

  Unsafe behavior
  ===============

  _UnsafeLowerTriangularMatrixView is unsafe because it assumes that the vector x remains valid during the usage of _UnsafeLowerTriangularMatrixView.

  ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

  _UnsafeLowerTriangularMatrixView(x::Vector{Float64}, N::Int)

  Create a new _UnsafeLowerTriangularMatrixView from x, zero the elements in x, and resize x if needed to ensure it has a length of at least N * (N + 1) / 2.

  Unsafe behavior
  ===============

  In addition to the usafe behavior of _UnsafeLowerTriangularMatrixView, this constructor is additionally unsafe because it may resize x. Only call it if you are sure that the usage of
  _UnsafeLowerTriangularMatrixView(x, N) is short-lived, and that there are no other views to x while the returned value is within scope.

That's probably why individual indexing of H does not seems to be slow here. This tutorial provides some additional details.

That being said, I just noticed that I filled all the entries of the lower-triangular part, and the docstring of @operator says that I "must fill in the non-zero lower-triangular entries only". I'm opening a PR to correct this.

Come to think of it, the destination matrix is not a SparseMatrixCSC, but the source A matrix will still be a SparseMatrixCSC, and I'm indexing to read the non-zero entries.

I need to do some benchmarks to verify if this part is slow...

Okay I did some more benchmarks. The indexing of A::SparseMatrixCSC for reading the lower-triangular non-zero entries is not the bottleneck here. As a matter of fact, I don't think there is much more room for improvements on my side with DI.jl computations and interfacing it with VectorNonlinearOracle/@operator API.

I'm not sure why it did not appear on my first profile above, but I ran multiple profiling and now, systematically, a fair amount of time is spent in MOI matrix coloring algorithms. See under thread "1 (interactive) and under initialize on the new profile here : profile.html. Is it expected @odow ? It feels weird that MOI spends a lot of time in its own coloring codebase when I'm delegating Hessian matrix coloring to DifferentiationInterface.jl.

edit: maybe also @amontoison can answer my question.

Maybe MOI does the coloring regardless of whether it actually has to compute a hessian?

Are there other scalar nonlinear constraints in the model? Our AD probably needs to kick in to compute the Hessian for those parts (or the objective, etc)

What's in the model:

• 1 @operator with a gradient and Hessian functions for the objective.
• 1 VectorNonlinearOracle with a Jacobian and a Hessian of the Lagrangian functions for nonlinear equality constraints
• 40 AffExpr in MOI.LessThan{Float64} for the linear equality constraints

I provide the Hessian/Hessian of the Langrangian function for all the nonlinear ingredients. I know that MOI still need to sum the 2 Hessians in the end, but AD should not be require for this, no ?

It is also worth mentioning that the Hessian sparsity pattern is different in the @operator and in the VectorNonlinearOracle, would it be the cause ?

1 @operator with a gradient and Hessian functions for the objective.

The AD in MOI runs for this. We could probably look at improving the performance of initialize if it's a bottleneck? It hasn't typically been an issue because most time is spent in the iterations.

The AD in MOI runs for this. We could probably look at improving the performance of initialize if it's a bottleneck? It hasn't typically been an issue because most time is spent in the iterations.

Oh that's why! Well here the performances with Uno are already quite good, about 4 ms per @optimize! call. But about 1/3 of the time is spent in initialize (in lime), 1/4 of the time in DI.jl and my user-defined functions (in red) and the rest in Uno iterrations (in grey):

So the MOI AD is a significant part of the profile, it could be probably improved, yes.

What I still don't get why AD is still needed if I provide user-defined functions for the gradient and Hessian of my nonlinear operator.

Another workaround on my side would be the replace the @operator with another NonlinearVectorOracle by using the epigraph reformulation min t; t >= f(x). The code would loose on readability IMO.

What I still don't get why AD is still needed if I provide user-defined functions for the gradient and Hessian of my nonlinear operator.

AD runs if there are any scalar nonlinear components, even if the result of the AD is just "query the operator, fill in the Hessian". We still need an algorithm for that. There's no special handling for exploiting the case where the constraint or objective body is a single user-defined operator.

Why is initialise being called each solve? Do you use parameters?

No parameters.

There are 35 solves in the profile above. I'm not sure if it's being called at each optimize! calls, or only at the first call and it takes that much of time over the 35 optimize! calls. I will verify right now.

Okay, funnily enough, initialize is not called at all if I do only 1 solve (1 optimize! call), but with 2 or more solves initialize is called at the 2nd solve, 3rd solve and so on. I call JuMP.set_normalized_rhs for the linear constraints before each solve, would it be the cause ?

MWE?

Here's what I'm benchmarking/profiling on Julia v1.12.1 and the latest stable releases of the packages:

using ModelPredictiveControl, JuMP
using UnoSolver

N = 35 # number of solves/optimize! calls

function f!(ẋ, x, u, _ , p)
    g, L, K, m = p          # [m/s²], [m], [kg/s], [kg]
    θ, ω = x[1], x[2]       # [rad], [rad/s]
    τ  = u[1]               # [Nm]
    ẋ[1] = ω
    ẋ[2] = -g/L*sin(θ) - K/m*ω + τ/m/L^2
    return nothing
end
h!(y, x, _ , _ ) = (y[1] = 180/π*x[1]; nothing) # [°]
p = [9.8, 0.4, 1.2, 0.3]
nu, nx, ny, Ts = 1, 2, 1, 0.1
model = NonLinModel(f!, h!, Ts, nu, nx, ny; p)
p_plant = copy(p); p_plant[3] = p[3]*1.25
plant = NonLinModel(f!, h!, Ts, nu, nx, ny; p=p_plant)
Hp, Hc, Mwt, Nwt = 20, 2, [0.5], [2.5]
α=0.01; σQ=[0.1, 1.0]; σR=[5.0]; nint_u=[1]; σQint_u=[0.1]
σQint_ym = zeros(0)
umin, umax = [-1.5], [+1.5]
nc = Hp+1
transcription = MultipleShooting()
optim = Model(()->UnoSolver.Optimizer(preset="filtersqp")); 
oracle = true
hessian = true
nmpc = NonLinMPC(model; 
    Hp, Hc, Mwt, Nwt, Cwt=Inf, transcription, oracle, hessian, optim,
    α, σQ, σR, nint_u, σQint_u, σQint_ym
)
nmpc = setconstraint!(nmpc; umin, umax)
unset_time_limit_sec(nmpc.optim)
sim!(nmpc, N, [180.0]; plant=plant, x_0=[0, 0], x̂_0=[0, 0, 0])
@time sim!(nmpc, N, [180.0]; plant=plant, x_0=[0, 0], x̂_0=[0, 0, 0])
@profview sim!(nmpc, N, [180.0]; plant=plant, x_0=[0, 0], x̂_0=[0, 0, 0])

Erratum: I made some improvements in DI.jl part of MPC.jl and I forget to register it.

The MWE is ran on MPC.jl v1.13.1. I'm registering it right now, but you can dev the package in the meantime.

There are a number of things going on here:

1. Modifying the RHS of a scalar constraint invalidates the inner model of Uno, so we
2. call initialise multiple times which
3. loops over every scalar constraint to compute the Hessian sparsity, which
4. is dense because we don't know the structure of the operator

We could fix this entirely in Uno (@amontoison) by supporting modifying the RHS.

Or we could improve MOI somewhat by trying to detect and exploit dense Hessians during the colouring.

I can support the modification of the RHS in UnoSolver.jl, but I think it would be more relevant to improve MOI instead.
I discussed this with @odow a long time ago, but I'm not able to find the related discussion or issue.
If the Jacobian or Hessian is dense and has to be specified or can be easily detected, it's not relevant to perform sparsity detection and coloring.
We can directly perform AD, which can by default run in vector mode in this case.

However, I am quite reluctant to change things in the nonlinear API of MOI.

It works, and all the coloring and sparsity-pattern detection were done back when MathProgBase.jl was still around.
I will give a talk about coloring at JuMP-dev, maybe we will get ideas about what we can easily improve.

Your take. I think it would make sense to support modifying the RHS of affine constraints in UnoSolver.jl. JuMP is designed to specialize on the problem structure, including the type of the constraints.

Oh but the RHS of affine / quadratic contraints is independent of the sparsity pattern, right?
We just update the lcon / ucon vectors.
If it is that, yes I should be quite easy to support in UnoSolver.jl.

Yes. I'm updating the RHS of affine inequality constraints with new constants. Or, in other words, updating the $b$ vector in $Ax \le b$ before calling optimize!.

About this, I think it's already supported in Ipopt.jl ? Because I never noticed this long initialize portion when I profiled Ipopt in the past.

The MOI interface in UnoSolver.jl is very similar to the one in Ipopt.jl.
Can you quickly check with Ipopt.jl.