Understanding why `MovingHorizonEstimator` Hessians are always fully dense

[franckgaga]

Hello @gdalle !

I'm having the same problem as #193 but for the MovingHorizonEstimator (MHE) instead of the NonLinMPC, that is, a fully-dense sparsity pattern for the Hessian of the objective function, when I'm 99.9% sure that many zeros are in fact global zeros. Note that I applied the equivalent improvements of #202 but for the MHE (on PR #216), so the problem seems to be elsewhere. I put too much time on trying to figure out why, so I'm asking for your help. What debugging method did you use to find out the two results mentioned here ? To be clearer, I want to handle structural sparsity better in the MHE, but I cannot find where the problem is in update_prediction! and obj_nonlinprog methods.

By using ModelPredictiveControl.jl v2.4.1 (registered in a few minutes, else you can dev it), here's a MRE:

using ModelPredictiveControl, ControlSystemsBase

Ts = 4.0
A =  [  0.800737  0.0       0.0  0.0
        0.0       0.606531  0.0  0.0
        0.0       0.0       0.8  0.0
        0.0       0.0       0.0  0.6    ]
Bu = [  0.378599  0.378599
        -0.291167  0.291167
        0.0       0.0
        0.0       0.0                   ]
Bd = [  0; 0; 0.5; 0.5;;                ]
C =  [  1.0  0.0  0.684   0.0
        0.0  1.0  0.0    -0.4736        ]
Dd = [  0.19; -0.148;;                  ]
Du = zeros(2,2)
model = LinModel(ss(A,[Bu Bd],C,[Du Dd],Ts),Ts,i_d=[3])
model = setop!(model, uop=[10,10], yop=[50,30], dop=[5])

using DifferentiationInterface, SparseConnectivityTracer, SparseMatrixColorings
import ForwardDiff
hessian = AutoSparse(AutoForwardDiff(), sparsity_detector=TracerSparsityDetector(), coloring_algorithm=GreedyColoringAlgorithm())
# hessian = AutoSparse(AutoForwardDiff(), sparsity_detector=TracerLocalSparsityDetector(), coloring_algorithm=GreedyColoringAlgorithm())

function gc!(LHS, X̂e, V̂e, Ŵe, Ue, Yem, De, P̄, x̄, p, ε)
    N = length(X̂e)÷6
    LHS .= 0
    for i in eachindex(LHS)
        if i ≤ N
            LHS[i] = X̂e[6*(i-1)+1] - 0 - ε
        end
    end
    return nothing
end

mhe = MovingHorizonEstimator(model; He=10, hessian, gc!, nc=11)
using JuMP; unset_time_limit_sec(mhe.optim)
res = sim!(mhe, 15, x̂_0=zeros(6), d_step=[-2.0])
info = getinfo(mhe)
display(info[:∇²J])

giving with TracerSparsityDetector:

66×66 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4356 stored entries:
⎡⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎤
⎢⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎥
⎢⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎥
⎢⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎥
⎢⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎥
⎢⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎥
⎢⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎥
⎢⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎥
⎢⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎥
⎢⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎥
⎢⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎥
⎢⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎥
⎢⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎥
⎣⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⎦

and with TracerLocalSparsityDetector:

66×66 SparseArrays.SparseMatrixCSC{Float64, Int64} with 2178 stored entries:
⎡⣕⣽⣪⣪⣗⣕⣕⣯⣪⣺⣕⣕⣽⣪⣪⣗⣕⣕⣯⣪⣺⣕⣕⣽⣪⣪⣗⣕⎤
⎢⡪⣺⢕⢕⡯⡪⡪⣗⢕⢽⡪⡪⣺⢕⢕⡯⡪⡪⣗⢕⢽⡪⡪⣺⢕⢕⡯⡪⎥
⎢⢝⢽⡫⡫⣟⢝⢝⡯⡫⣻⢝⢝⢽⡫⡫⣟⢝⢝⡯⡫⣻⢝⢝⢽⡫⡫⣟⢝⎥
⎢⡵⣽⢮⢮⡷⡵⡵⣯⢮⢾⡵⡵⣽⢮⢮⡷⡵⡵⣯⢮⢾⡵⡵⣽⢮⢮⡷⡵⎥
⎢⣪⣺⣕⣕⣯⣪⣪⣗⣕⣽⣪⣪⣺⣕⣕⣯⣪⣪⣗⣕⣽⣪⣪⣺⣕⣕⣯⣪⎥
⎢⢕⢽⡪⡪⣗⢕⢕⡯⡪⣺⢕⢕⢽⡪⡪⣗⢕⢕⡯⡪⣺⢕⢕⢽⡪⡪⣗⢕⎥
⎢⡳⣻⢞⢞⡷⡳⡳⣟⢞⢾⡳⡳⣻⢞⢞⡷⡳⡳⣟⢞⢾⡳⡳⣻⢞⢞⡷⡳⎥
⎢⢮⢾⡵⡵⣯⢮⢮⡷⡵⣽⢮⢮⢾⡵⡵⣯⢮⢮⡷⡵⣽⢮⢮⢾⡵⡵⣯⢮⎥
⎢⢕⢽⡪⡪⣗⢕⢕⡯⡪⣺⢕⢕⢽⡪⡪⣗⢕⢕⡯⡪⣺⢕⢕⢽⡪⡪⣗⢕⎥
⎢⡫⣻⢝⢝⡯⡫⡫⣟⢝⢽⡫⡫⣻⢝⢝⡯⡫⡫⣟⢝⢽⡫⡫⣻⢝⢝⡯⡫⎥
⎢⢞⢾⡳⡳⣟⢞⢞⡷⡳⣻⢞⢞⢾⡳⡳⣟⢞⢞⡷⡳⣻⢞⢞⢾⡳⡳⣟⢞⎥
⎢⣕⣽⣪⣪⣗⣕⣕⣯⣪⣺⣕⣕⣽⣪⣪⣗⣕⣕⣯⣪⣺⣕⣕⣽⣪⣪⣗⣕⎥
⎢⡪⣺⢕⢕⡯⡪⡪⣗⢕⢽⡪⡪⣺⢕⢕⡯⡪⡪⣗⢕⢽⡪⡪⣺⢕⢕⡯⡪⎥
⎣⢝⢽⡫⡫⣟⢝⢝⡯⡫⣻⢝⢝⢽⡫⡫⣟⢝⢝⡯⡫⣻⢝⢝⢽⡫⡫⣟⢝⎦

Many thanks for your time !

[gdalle]

If this is due to the many zeros in your input data, SCT won't be able to do much (as you recall, we eventually reverted the changes inspired by #193). Have you tried storing A as a Diagonal for instance, instead of an artificially dense matrix?

[franckgaga]

Yes, the PR #216 introduced the changes to preserve Diagonal matrices for the weights. So the objective function is a sum of three dot(x, A, x), and the A argument is a sparse Diagonal matrix for two of them.

Now that I mention this, the first one, dot(x̄, invP̄, x̄), is with a dense matrix, it may be the cause. I will investigate this avenue. Edit: btw, invP̄ is really a dense (hermitian) matrix here, I cannot store it as a Diagonal.

Edit 2: I'm referring to the three dot calls here:

ModelPredictiveControl.jl/src/estimator/mhe/execute.jl

Line 796 in 5628046

return dot(x̄, invP̄, x̄) + dot(Ŵ, invQ̂_Nk, Ŵ) + dot(V̂, invR̂_Nk, V̂) + Jε

[franckgaga]

Oh sorry, the MRE example above is with a linear model, so it won't call this method of obj_nonlinprog. My bad.

Here's a better MRE with a NonLinModel, to be coherent with my explanation above and maybe simplify the debugging (it leads to the same results as above):

using ModelPredictiveControl, ControlSystemsBase

Ts = 4.0
A =  [  0.800737  0.0       0.0  0.0
        0.0       0.606531  0.0  0.0
        0.0       0.0       0.8  0.0
        0.0       0.0       0.0  0.6    ]
Bu = [  0.378599  0.378599
        -0.291167  0.291167
        0.0       0.0
        0.0       0.0                   ]
Bd = [  0; 0; 0.5; 0.5;;                ]
C =  [  1.0  0.0  0.684   0.0
        0.0  1.0  0.0    -0.4736        ]
Dd = [  0.19; -0.148;;                  ]
Du = zeros(2,2)
f(x,u,d,p) = p.A*x + p.Bu*u + p.Bd*d
h(x,d,p)   = p.C*x + p.Dd*d
model = NonLinModel(f, h, Ts, 2, 4, 2, 1, solver=nothing, p=(;A,Bu,Bd,C,Dd))
model = setop!(model, uop=[10,10], yop=[50,30], dop=[5])

using DifferentiationInterface, SparseConnectivityTracer, SparseMatrixColorings
import ForwardDiff
hessian = AutoSparse(AutoForwardDiff(), sparsity_detector=TracerSparsityDetector(), coloring_algorithm=GreedyColoringAlgorithm())
# hessian = AutoSparse(AutoForwardDiff(), sparsity_detector=TracerLocalSparsityDetector(), coloring_algorithm=GreedyColoringAlgorithm())

function gc!(LHS, X̂e, V̂e, Ŵe, Ue, Yem, De, P̄, x̄, p, ε)
    N = length(X̂e)÷6
    LHS .= 0
    for i in eachindex(LHS)
        if i ≤ N
            LHS[i] = X̂e[6*(i-1)+1] - 0 - ε
        end
    end
    return nothing
end

mhe = MovingHorizonEstimator(model; He=10, hessian, gc!, nc=11)
using JuMP; unset_time_limit_sec(mhe.optim)
res = sim!(mhe, 15, x̂_0=zeros(6), d_step=[-2.0])
info = getinfo(mhe)
display(info[:∇²J])

[franckgaga]

Now that I mention this, the first one, dot(x̄, invP̄, x̄), is with a dense matrix, it may be the cause. I will investigate this avenue. Edit: btw, invP̄ is really a dense (hermitian) matrix here, I cannot store it as a Diagonal.

So I just replaced:

return dot(x̄, invP̄, x̄) + dot(Ŵ, invQ̂_Nk, Ŵ) + dot(V̂, invR̂_Nk, V̂) + Jε

with:

return dot(Ŵ, invQ̂_Nk, Ŵ) + dot(V̂, invR̂_Nk, V̂)

and same issue. So this is not caused by the dense invP̄ matrix, nor the Jε term. And adding println(typeof(invR̂_Nk)) and println(typeof(invQ̂_Nk)) before the return really shows that the two matrices are sparse LinearAlgebra.Hermitian{Float64, LinearAlgebra.Diagonal{Float64, Vector{Float64}}}. Note sure what is the cause...