# Least squares fitting that accounts for covariance matrix of all input parameters

**URL:** <https://discourse.julialang.org/t/least-squares-fitting-that-accounts-for-covariance-matrix-of-all-input-parameters/64775>\
**Category:** Optimization (Mathematical)\
**Tags:** curve-fitting\
**Created:** [July 16, 2021, 9:21pm UTC](https://discourse.julialang.org/t/least-squares-fitting-that-accounts-for-covariance-matrix-of-all-input-parameters/64775 "2021-07-16T21:21:11Z")\
**Posts on this page:** 3\
**Page:** 2

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [July 29, 2021, 8:05pm UTC](https://discourse.julialang.org/t/least-squares-fitting-that-accounts-for-covariance-matrix-of-all-input-parameters/64775/21 "2021-07-29T20:05:38Z")

</div>

> [@cgeoga](#):
>
> I have tried to get Hessians out of JuMP objects and it’s pretty brutal. It might be easiest to just write the objective as a standalone function

@cgeoga is unfortunately correct.

There’s limited documentation and you’ll have to get pretty low-level:

> **[Nonlinear Modeling (Legacy) · JuMP](https://jump.dev/JuMP.jl/stable/manual/nlp/#Querying-derivatives-from-a-JuMP-model)**
>
> Documentation for JuMP.

---

<div class="post-metadata">

**Author:** ![nluetts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nluetts/32/24967_2.png) [@nluetts](https://discourse.julialang.org/u/nluetts)\
**Post date:** [July 29, 2021, 11:39pm UTC](https://discourse.julialang.org/t/least-squares-fitting-that-accounts-for-covariance-matrix-of-all-input-parameters/64775/22 "2021-07-29T23:39:25Z")

</div>

Thank you @cgeoga and @odow !

> [@cgeoga](#):
>
> It might be easiest to just write the objective as a standalone function and then just call `ForwardDiff.hessian(obj, mle)` with the computed `mle` .

I am not sure how to do this in a way that accounts for the constraints. The objective function itself does not depend on the unknown parameters, so its hessian would not be enough. Maybe taking the hessian of the function Φ in the paper mentioned above (page 3, objective function + Lagrangian multiplier term) would do the trick, but I am missing the mathematical background and are merely guessing 🤷‍♂️ Plus, I do not have the values of the Lagrangian multipliers.

> [@odow](#):
>
> There’s limited documentation and you’ll have to get pretty low-level:  
> [https://jump.dev/JuMP.jl/stable/manual/nlp/#Querying-derivatives-from-a-JuMP-model](https://jump.dev/JuMP.jl/stable/manual/nlp/#Querying-derivatives-from-a-JuMP-model)

Thank you for pointing me to the right place in the docs. Based on the given examples, I tried the following:

```julia
# get variable values in correct order
vars = [zs; ζs; [f, A, m1, m2, m3, m4]]
indx = [JuMP.index(v).value for v in vars]
vals = zeros(length(vars))
for (v, i) in zip(vars, indx)
    vals[i] = value(v) # looks good, order is just as I entered the variables in the model
end

# make evaluator object
evl = NLPEvaluator(model)
# request hessian calculation
MathOptInterface.initialize(evl, [:Hess])
# get index to place elements from vector returned by eval_hessian_lagrangian into a matrix
hess_struct = MathOptInterface.hessian_lagrangian_structure(evl)
# flat vector, to be filled with hessian elements
hess_flat = zeros(length(hess_struct))
# retrieve the hessian
MathOptInterface.eval_hessian_lagrangian(evl, hess_flat, vals, 1.0, ones(19))
# reshape from vector to matrix
hessian = zeros((length(vars), length(vars)))
for (i, (k, l)) in enumerate(hess_struct)
    hessian[k, l] = hess_flat[i]
end

```

It runs through w/o an error, but the numbers I get span many orders of magnitude, which seems wrong:

```julia
# print diagonal elements together with corresponding variables
for (v, h) in zip(vars, [hessian[i, i] for i in 1:52])
    println(v, " ", h)
end

```

```julia
Iₒ[1] 3.780718336483931e9
Iₒ[2] 3.780718336483932e9
Iₒ[3] 3.7807183364839325e9
Iₒ[4] 3.780718336483932e9
Iₒ[5] 3.780718336483932e9
Iₒ[6] 3.780718336483931e9
Iₒ[7] 3.7807183364839315e9
Iₒ[8] 3.7807183364839315e9
Iₒ[9] 3.780718336483932e9
Iₒ[10] 3.780718336483931e9
Iₒ[11] 3.7807183364839315e9
Iₒ[12] 3.7807183364839325e9
Iₒ[13] 3.7807183364839315e9
Iₒ[14] 3.780718336483931e9
Iₒ[15] 3.780718336483931e9
Iₒ[16] 3.780718336483932e9
Irₒ[1] 3.780718336483931e9
Irₒ[2] 3.780718336483932e9
msₒ 2.0e10
mrₒ 3.125e10
ρₒ 23.781212841854938
ρrₒ 3.9674667724657806
aₒ 163265.306122449
I[1] 8.580578122520386e-9
I[2] 8.580578122520386e-9
I[3] 8.580578122520386e-9
I[4] 8.580578122520386e-9
I[5] 8.580578122520386e-9
I[6] 8.580578122520386e-9
I[7] 8.580578122520386e-9
I[8] 8.580578122520386e-9
I[9] 8.580578122520386e-9
I[10] 8.580578122520386e-9
I[11] 8.580578122520386e-9
I[12] 8.580578122520386e-9
I[13] 8.580578122520386e-9
I[14] 8.580578122520386e-9
I[15] 8.580578122520386e-9
I[16] 8.580578122520386e-9
Ir[1] 8.580578122520386e-9
Ir[2] 8.580578122520386e-9
ms 2.0e10
mr 0.0
ρ 3.9568198152688134e-12
ρr 3.9674667724657806
a 0.0
f 0.0
A 0.0
m1 0.0
m2 0.0
m3 0.0
m4 0.0

```

The returned hessian cannot be inverted, I get a LAPACKException.

Maybe I did something obviously wrong, but I don’t see it and my missing background knowledge in this field is not helping … 🙂 I am also not sure what the weighting factors passed to [`eval_hessian_lagrangian`](https://jump.dev/MathOptInterface.jl/dev/reference/nonlinear/#MathOptInterface.eval_hessian_lagrangian) (sigma and mu) should be, so I set them to 1.  
I noticed that I can manipulate these values to “compress” the range that the values in the hessian span, but even if I do, I still cannot invert it (plus I am pretty sure I cannot just plug in whatever values I like here).

---

<div class="post-metadata">

**Author:** ![nluetts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nluetts/32/24967_2.png) [@nluetts](https://discourse.julialang.org/u/nluetts)\
**Post date:** [August 4, 2021, 10:24am UTC](https://discourse.julialang.org/t/least-squares-fitting-that-accounts-for-covariance-matrix-of-all-input-parameters/64775/23 "2021-08-04T10:24:39Z")

</div>

Since I did not make considerable progress on the JuMP and ModelingToolkit fronts, I post here my “homebrew” solution which does what I want:

[https://github.com/nluetts/ConstrainedLeastSquares.jl](https://github.com/nluetts/ConstrainedLeastSquares.jl)

It is basically the code fragment I posted above together with some convenience functionality to declare the input parameters and model more easily.  
I don’t suspect it to be very performant on larger problems, but for my personal use it is good enough. Perhaps it will be helpful for others as well.

[Previous page](https://discourse.julialang.org/t/least-squares-fitting-that-accounts-for-covariance-matrix-of-all-input-parameters/64775.md?page=1)