The source files for all examples can be found in /examples.

Closest Correlation Matrix

We consider the problem of finding the closest correlation matrix $X$ to a given random matrix $C$. With closest correlation matrix we mean a positive semidefinite matrix with ones on the diagonal. The problem is given by:

\[\begin{array}{ll} \text{minimize} & \frac{1}{2}||X - C||_F^2\\ \text{subject to} & X_{ii} = 1, \quad i=1,\dots,n \\ & X \succeq 0. \end{array}\]

Notice that we use JuMP to model the problem. COSMO is chosen as the backend solver. And COSMO-specific settings are passed using the optimizer_with_attributes() function.

using COSMO, JuMP, LinearAlgebra, SparseArrays, Test, Random
rng = Random.MersenneTwister(12345);
# create a random test matrix C
n = 8;
C = -1 .+ rand(rng, n, n) .* 2;
c = vec(C);

Define problem in JuMP:

q = -vec(C);
r = 0.5 * vec(C)' * vec(C);
m = JuMP.Model(optimizer_with_attributes(COSMO.Optimizer, "verbose" => true));
@variable(m, X[1:n, 1:n], PSD);
x = vec(X);
@objective(m, Min, 0.5 * x' * x  + q' * x + r);
for i = 1:n
  @constraint(m, X[i, i] == 1.);
end

Solve the JuMP model with COSMO and query the solution X_sol:

status = JuMP.optimize!(m);
obj_val = JuMP.objective_value(m);
X_sol = JuMP.value.(X);
Minimize ScalarQuadraticFunction{Float64}:
 9.153897960749184 - 0.5868445948072374 X[1,1] - 0.6979328738566593 X[1,2] - 0.7086271918729499 X[1,3] - 0.34791857373710666 X[1,4] + 0.39881354761120935 X[1,5] - 1.5157405234804902 X[1,6] - 0.2779934009144056 X[1,7] - 0.724869793873145 X[1,8] + 0.0007824842282477817 X[2,2] + 1.1688231933811437 X[2,3] + 0.39138130091300827 X[2,4] + 1.1754088610406903 X[2,5] - 0.49033659501878457 X[2,6] - 0.2433110958503466 X[2,7] - 0.0888332573675128 X[2,8] - 0.5841006783457146 X[3,3] - 1.0859021758454164 X[3,4] - 0.21982762884318996 X[3,5] + 0.77349753680036 X[3,6] + 0.006240353381980146 X[3,7] + 0.35917684338625033 X[3,8] - 0.03464464750524243 X[4,4] + 0.7620066249958741 X[4,5] - 1.0971751275764996 X[4,6] - 0.4168944118970974 X[4,7] - 1.3903367522952355 X[4,8] - 0.3781801564883356 X[5,5] - 0.512564646198213 X[5,6] - 0.00030128280174812616 X[5,7] - 0.3510799342418949 X[5,8] + 0.6507354932966494 X[6,6] - 0.7986461394322726 X[6,7] + 0.30662204970370954 X[6,8] - 0.5281451186073247 X[7,7] + 1.3797582344661254 X[7,8] + 0.7932198369405121 X[8,8] + 0.5 X[1,1]² + 1.0 X[1,2]² + 1.0 X[1,3]² + 1.0 X[1,4]² + 1.0 X[1,5]² + 1.0 X[1,6]² + 1.0 X[1,7]² + 1.0 X[1,8]² + 0.5 X[2,2]² + 1.0 X[2,3]² + 1.0 X[2,4]² + 1.0 X[2,5]² + 1.0 X[2,6]² + 1.0 X[2,7]² + 1.0 X[2,8]² + 0.5 X[3,3]² + 1.0 X[3,4]² + 1.0 X[3,5]² + 1.0 X[3,6]² + 1.0 X[3,7]² + 1.0 X[3,8]² + 0.5 X[4,4]² + 1.0 X[4,5]² + 1.0 X[4,6]² + 1.0 X[4,7]² + 1.0 X[4,8]² + 0.5 X[5,5]² + 1.0 X[5,6]² + 1.0 X[5,7]² + 1.0 X[5,8]² + 0.5 X[6,6]² + 1.0 X[6,7]² + 1.0 X[6,8]² + 0.5 X[7,7]² + 1.0 X[7,8]² + 0.5 X[8,8]²

Subject to:

VectorAffineFunction{Float64}-in-Zeros
 ┌                 ┐
 │-1.0 + 1.0 X[1,1]│
 └                 ┘ ∈ Zeros(1)
 ┌                 ┐
 │-1.0 + 1.0 X[2,2]│
 └                 ┘ ∈ Zeros(1)
 ┌                 ┐
 │-1.0 + 1.0 X[3,3]│
 └                 ┘ ∈ Zeros(1)
 ┌                 ┐
 │-1.0 + 1.0 X[4,4]│
 └                 ┘ ∈ Zeros(1)
 ┌                 ┐
 │-1.0 + 1.0 X[5,5]│
 └                 ┘ ∈ Zeros(1)
 ┌                 ┐
 │-1.0 + 1.0 X[6,6]│
 └                 ┘ ∈ Zeros(1)
 ┌                 ┐
 │-1.0 + 1.0 X[7,7]│
 └                 ┘ ∈ Zeros(1)
 ┌                 ┐
 │-1.0 + 1.0 X[8,8]│
 └                 ┘ ∈ Zeros(1)

VectorAffineFunction{Float64}-in-Scaled{PositiveSemidefiniteConeTriangle}
 ┌                               ┐
 │0.0 + 1.0 X[1,1]               │
 │0.0 + 1.4142135623730951 X[1,2]│
 │0.0 + 1.0 X[2,2]               │
 │0.0 + 1.4142135623730951 X[1,3]│
 │0.0 + 1.4142135623730951 X[2,3]│
 │0.0 + 1.0 X[3,3]               │
 │0.0 + 1.4142135623730951 X[1,4]│
 │0.0 + 1.4142135623730951 X[2,4]│
 │0.0 + 1.4142135623730951 X[3,4]│
 │0.0 + 1.0 X[4,4]               │
 │0.0 + 1.4142135623730951 X[1,5]│
 │0.0 + 1.4142135623730951 X[2,5]│
 │0.0 + 1.4142135623730951 X[3,5]│
 │0.0 + 1.4142135623730951 X[4,5]│
 │0.0 + 1.0 X[5,5]               │
 │0.0 + 1.4142135623730951 X[1,6]│
 │0.0 + 1.4142135623730951 X[2,6]│
 │0.0 + 1.4142135623730951 X[3,6]│
 │0.0 + 1.4142135623730951 X[4,6]│
 │0.0 + 1.4142135623730951 X[5,6]│
 │0.0 + 1.0 X[6,6]               │
 │0.0 + 1.4142135623730951 X[1,7]│
 │0.0 + 1.4142135623730951 X[2,7]│
 │0.0 + 1.4142135623730951 X[3,7]│
 │0.0 + 1.4142135623730951 X[4,7]│
 │0.0 + 1.4142135623730951 X[5,7]│
 │0.0 + 1.4142135623730951 X[6,7]│
 │0.0 + 1.0 X[7,7]               │
 │0.0 + 1.4142135623730951 X[1,8]│
 │0.0 + 1.4142135623730951 X[2,8]│
 │0.0 + 1.4142135623730951 X[3,8]│
 │0.0 + 1.4142135623730951 X[4,8]│
 │0.0 + 1.4142135623730951 X[5,8]│
 │0.0 + 1.4142135623730951 X[6,8]│
 │0.0 + 1.4142135623730951 X[7,8]│
 │0.0 + 1.0 X[8,8]               │
 └                               ┘ ∈ Scaled{PositiveSemidefiniteConeTriangle}(PositiveSemidefiniteConeTriangle(8))

------------------------------------------------------------------
          COSMO v0.8.10 - A Quadratic Objective Conic Solver
                         Michael Garstka
                University of Oxford, 2017 - 2022
------------------------------------------------------------------

Problem:  x ∈ R^{36},
          constraints: A ∈ R^{44x36} (44 nnz),
          matrix size to factor: 80x80,
          Floating-point precision: Float64
Sets:     DensePsdConeTriangle of dim: 36 (8x8)
          ZeroSet of dim: 8
Settings: ϵ_abs = 1.0e-05, ϵ_rel = 1.0e-05,
          ϵ_prim_inf = 1.0e-04, ϵ_dual_inf = 1.0e-04,
          ρ = 0.1, σ = 1e-06, α = 1.6,
          max_iter = 5000,
          scaling iter = 10 (on),
          check termination every 25 iter,
          check infeasibility every 40 iter,
          KKT system solver: QDLDL
Acc:      Anderson Type2{QRDecomp},
          Memory size = 15, RestartedMemory,
          Safeguarded: true, tol: 2.0
Setup Time: 845.84ms

Iter:	Objective:	Primal Res:	Dual Res:	Rho:
1	 5.7298e+00	5.9277e-01	7.1084e-01	1.0000e-01
25	-3.9953e-01	4.6201e-06	2.1964e-06	1.0000e-01

------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 25
Optimal objective: -0.3995
Runtime: 4.278s (4277.85ms)

Double check result against known solution:

known_opt_val = 8.75427
known_solution =  [1.0         0.354428    0.238389    0.323015  -0.117854    0.620975    0.0969092   0.231657
  0.354428    1.0        -0.590998   -0.172179  -0.571414    0.232656    0.109805    0.0229169
  0.238389   -0.590998    1.0         0.379909   0.0205233  -0.237496    0.0431976  -0.0364916
  0.323015   -0.172179    0.379909    1.0       -0.250821    0.349762    0.13669     0.49255
 -0.117854   -0.571414    0.0205233  -0.250821   1.0         0.145922   -0.0418897   0.0604382
  0.620975    0.232656   -0.237496    0.349762   0.145922    1.0         0.457861    0.0215352
  0.0969092   0.109805    0.0431976   0.13669   -0.0418897   0.457861    1.0        -0.626215
  0.231657    0.0229169  -0.0364916   0.49255    0.0604382   0.0215352  -0.626215    1.0   ];

@test isapprox(obj_val, known_opt_val , atol=1e-3)
Test Passed
@test norm(X_sol - known_solution, Inf) < 1e-3
Test Passed

This page was generated using Literate.jl.