# JuMP: Squared Frobenius norm objective
**URL:** https://discourse.julialang.org/t/jump-squared-frobenius-norm-objective/114400
**Category:** Optimization (Mathematical)
**Tags:** question, jump
**Created:** [May 17, 2024, 3:07pm UTC](https://discourse.julialang.org/t/jump-squared-frobenius-norm-objective/114400 "2024-05-17T15:07:14Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@smallpondtom](https://discourse.julialang.org/u/smallpondtom)
#### Post date: [May 17, 2024, 3:07pm UTC](https://discourse.julialang.org/t/jump-squared-frobenius-norm-objective/114400/1 "2024-05-17T15:07:14Z")
Hi. I can’t seem to find a definitive answer so I’m posting here for help. What is the standard and latest way in JuMP to construct an objective function for a Frobenius norm?
Say I have a (fictional) minimization of
\min\_{P\in\mathbb{R}^{n\times n}} \| \mathbf{PAX} - \mathbf{QX} \|\_F^2
where \mathbf{X}\in\mathbb{R}^{n\times K}, K \> n, is some data matrix. And, the matrices \mathbf{A}\in\mathbb{R}^{n\times n}, \mathbf{Q}=\mathbf{Q}^\top\in\mathbb{R}^{n\times n} are given. Also assume \mathbf{P} = \mathbf{P}^\top.
Would the code to solve this be like the following?
```julia
n = 5
K = 100
X = rand(n, K)
A = rand(5,5)
Q = rand(5,5)
Q *= Q'
model = Model(Ipopt.Optimizer)
@variable(model, P[1:n, 1:n], Symmetric)
@objective(model, sum((P*A*X - Q*X).^2))
optimize!(model)
```
Or do I have to use the `MOI.SecondOrderCone()` as in the example below (for a vector case)?
```julia
@constraint(model, [t; x] in SecondOrderCone())
```
---
### Author:  [@jd-foster](https://discourse.julialang.org/u/jd-foster)
#### Post date: [May 17, 2024, 4:08pm UTC](https://discourse.julialang.org/t/jump-squared-frobenius-norm-objective/114400/2 "2024-05-17T16:08:46Z")
I think what you’ve done is reasonable.
Below is how I would write your example, plus a `SecondOrderCone` version as you suggested:
```julia
using JuMP
import Ipopt
import LinearAlgebra
import SCS
n = 5
K = 100
## Generate data:
X = rand(n, K)
A = rand(5,5)
Q = LinearAlgebra.Symmetric(rand(5,5))
## Nonlinear version:
solver = Ipopt
model = Model(solver.Optimizer)
@variable(model, P[1:n, 1:n], Symmetric)
@objective(model, Min, sum((P*A*X - Q*X).^2));
optimize!(model)
solution_summary(model)
z_nl = objective_value(model)
P_nl = value.(P)
## SoC version:
solver = SCS
model = Model(solver.Optimizer)
@variable(model, P[1:n, 1:n], Symmetric)
@variable(model, t)
@constraint(model, [t; vec(P*A*X - Q*X)] in SecondOrderCone());
@objective(model, Min, t)
optimize!(model)
solution_summary(model)
z_soc = objective_value(model)^2
P_soc = value.(P)
## Compare
isapprox(z_soc, z_nl, rtol=1e-3)
LinearAlgebra.norm(P_soc - P_nl) <= 1e-3
```
While you’ve probably just given a simplified problem, this type of problem has an “exact” solution using the [Moore–Penrose inverse](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse):
(see Ben-Israel, A., & Greville, T. N. (2003). _Generalized inverses: theory and applications_ (Vol. 15). Springer)
Also the [Ref 636](https://www.cambridge.org/core/journals/mathematical-proceedings-of-the-cambridge-philosophical-society/article/on-best-approximate-solutions-of-linear-matrix-equations/FDB8762BB13CB93E9E70080B2CFD986A)
```julia
## "Exact" version
P_pi = Q*X*LinearAlgebra.pinv(A*X)
z_pi = LinearAlgebra.norm(P_pi*A*X - Q*X, 2)
# Compare to:
LinearAlgebra.norm(P_soc*A*X - Q*X, 2)
LinearAlgebra.norm(P_nl *A*X - Q*X, 2)
```
Hopefully you observe how much smaller the residual is in the exact pseudoinverse version.
---
### Author:  [@smallpondtom](https://discourse.julialang.org/u/smallpondtom)
#### Post date: [May 18, 2024, 11:40pm UTC](https://discourse.julialang.org/t/jump-squared-frobenius-norm-objective/114400/3 "2024-05-18T23:40:24Z")
@jd-foster Thanks for the detailed explanation. And, yes I agree for this simplified example an exact solution does exist. Thanks for that clarification as well.
Do you know if using the `SecondOrderCone()` would make the optimization more or less efficient (in the case of using a solver like SCS)? I’m curious about what the most efficient implementation would be like.
---
### Author:  [@jd-foster](https://discourse.julialang.org/u/jd-foster)
#### Post date: [May 20, 2024, 1:39am UTC](https://discourse.julialang.org/t/jump-squared-frobenius-norm-objective/114400/4 "2024-05-20T01:39:58Z")
In the case of SCS, I would guess that the second one would be most efficient since the solver targets cones in the algorithm design. Using a general nonlinear solver like Ipopt, the first form is best I believe.
The way to really know would be to benchmark over increasing test problem sizes while tracking the accuracy against the known solution.
@odow Thoughts?
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [May 20, 2024, 2:36am UTC](https://discourse.julialang.org/t/jump-squared-frobenius-norm-objective/114400/5 "2024-05-20T02:36:59Z")
You probably want `RotatedSecondOrderCone` instead of `SecondOrderCone`. See [Tips and Tricks · JuMP](https://jump.dev/JuMP.jl/stable/tutorials/conic/tips_and_tricks/#Rotated-Second-Order-Cone).
“Most efficient” is probably problem-dependent. You should just try a few different formulations. Ipopt.jl does not support the conic formulation.
```Julia
using JuMP
import Ipopt
import LinearAlgebra
import SCS
function test_ipopt(X, A, Q)
n = size(X, 1)
model = Model(Ipopt.Optimizer)
set_silent(model)
@variable(model, P[1:n, 1:n], Symmetric)
@objective(model, Min, sum((P*A*X - Q*X).^2));
optimize!(model)
@assert is_solved_and_feasible(model)
return objective_value(model), value.(P)
end
function test_ipopt2(X, A, Q)
n, K = size(X)
model = Model(Ipopt.Optimizer)
set_silent(model)
@variable(model, P[1:n, 1:n], Symmetric)
@variable(model, residuals[1:n, 1:K])
@constraint(model, residuals .== P*A*X - Q*X)
@objective(model, Min, sum(residuals.^2));
optimize!(model)
@assert is_solved_and_feasible(model)
return objective_value(model), value.(P)
end
function test_scs(X, A, Q)
n = size(X, 1)
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, P[1:n, 1:n], Symmetric)
@variable(model, t)
@constraint(model, [t; 0.5; vec(P*A*X - Q*X)] in RotatedSecondOrderCone());
@objective(model, Min, t)
optimize!(model)
return objective_value(model), value.(P)
end
n, K = 5, 100
X = rand(n, K);
A = rand(5, 5);
Q = LinearAlgebra.Symmetric(rand(5, 5))
@time test_ipopt(X, A, Q)
@time test_ipopt2(X, A, Q)
@time test_scs(X, A, Q)
```
---
### Author:  [@smallpondtom](https://discourse.julialang.org/u/smallpondtom)
#### Post date: [May 20, 2024, 1:32pm UTC](https://discourse.julialang.org/t/jump-squared-frobenius-norm-objective/114400/6 "2024-05-20T13:32:32Z")
@odow For `Ipopt`, since I’m dealing with high-dimensionalities, constraining the `residual` as in the second example you gave me took much less time to set up the optimization problem compared to the first one. I’m guessing this is because `Ipopt` works well with constraints.
Thanks for the tips!