# Nearest correlation

This example illustrates the sensitivity analysis of the nearest correlation problem studied in [H02].

Higham, Nicholas J. Computing the nearest correlation matrix—a problem from finance. IMA journal of Numerical Analysis 22.3 (2002): 329-343.

using DiffOpt, JuMP, SCS, LinearAlgebra
solver = SCS.Optimizer

function proj(A, dH = Diagonal(ones(size(A, 1))), H = ones(size(A)))
n = LinearAlgebra.checksquare(A)
model = Model(() -> DiffOpt.diff_optimizer(solver))
@variable(model, X[1:n, 1:n] in PSDCone())
@constraint(model, [i in 1:n], X[i, i] == 1)
@objective(model, Min, sum((H .* (X - A)) .^ 2))
MOI.set(
model,
DiffOpt.ForwardObjectiveFunction(),
sum((dH .* (X - A)) .^ 2),
)
optimize!(model)
DiffOpt.forward_differentiate!(model)
dX = MOI.get.(model, DiffOpt.ForwardVariablePrimal(), X)
return value.(X), dX
end
proj (generic function with 3 methods)

Example from [H02, p. 334-335]:

A = LinearAlgebra.Tridiagonal(ones(2), ones(3), ones(2))
3×3 LinearAlgebra.Tridiagonal{Float64, Vector{Float64}}:
1.0  1.0   ⋅
1.0  1.0  1.0
⋅   1.0  1.0

The projection is computed as follows:

X, dX = proj(A)
------------------------------------------------------------------
SCS v3.2.4 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 6, constraints m: 9
cones: 	  z: primal zero / dual free vars: 3
s: psd vars: 6, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
max_iters: 100000, normalize: 1, rho_x: 1.00e-06
acceleration_lookback: 10, acceleration_interval: 10
compiled with openmp parallelization enabled
lin-sys:  sparse-direct-amd-qdldl
nnz(A): 9, nnz(P): 6
------------------------------------------------------------------
iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
0| 9.98e-01  3.35e+00  2.27e+01 -1.27e+01  1.00e-01  1.14e-04
75| 5.98e-05  2.10e-06  1.25e-04 -6.72e+00  1.00e-01  3.60e-04
------------------------------------------------------------------
status:  solved
timings: total: 3.61e-04s = setup: 4.59e-05s + solve: 3.15e-04s
lin-sys: 2.36e-05s, cones: 2.02e-04s, accel: 3.92e-06s
------------------------------------------------------------------
objective = -6.721507
------------------------------------------------------------------

The projection of A is:

X
3×3 Matrix{Float64}:
1.0       0.760748  0.157264
0.760748  1.0       0.760748
0.157264  0.760748  1.0

The derivative of the projection with respect to a uniform increase of the weights of the diagonal entries is:

dX
3×3 Matrix{Float64}:
6.73573e-8  -0.320654     0.388481
-0.320654     1.03972e-7  -0.320654
0.388481    -0.320654     7.60728e-8

Example from [H02, Section 4, p. 340]:

A = LinearAlgebra.Tridiagonal(-ones(3), 2ones(4), -ones(3))
4×4 LinearAlgebra.Tridiagonal{Float64, Vector{Float64}}:
2.0  -1.0    ⋅     ⋅
-1.0   2.0  -1.0    ⋅
⋅   -1.0   2.0  -1.0
⋅     ⋅   -1.0   2.0

The projection is computed as follows:

X, dX = proj(A)
------------------------------------------------------------------
SCS v3.2.4 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 10, constraints m: 14
cones: 	  z: primal zero / dual free vars: 4
s: psd vars: 10, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
max_iters: 100000, normalize: 1, rho_x: 1.00e-06
acceleration_lookback: 10, acceleration_interval: 10
compiled with openmp parallelization enabled
lin-sys:  sparse-direct-amd-qdldl
nnz(A): 14, nnz(P): 10
------------------------------------------------------------------
iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
0| 1.02e+00  3.57e+00  3.24e+01 -3.42e+01  1.00e-01  8.84e-05
75| 2.52e-05  7.03e-07  7.14e-05 -1.74e+01  1.00e-01  4.12e-04
------------------------------------------------------------------
status:  solved
timings: total: 4.13e-04s = setup: 4.12e-05s + solve: 3.72e-04s
lin-sys: 3.09e-05s, cones: 2.58e-04s, accel: 4.49e-06s
------------------------------------------------------------------
objective = -17.447239
------------------------------------------------------------------

The projection of A is:

X
4×4 Matrix{Float64}:
1.0       -0.808425   0.191576   0.10677
-0.808425   1.0       -0.656259   0.191576
0.191576  -0.656259   1.0       -0.808425
0.10677    0.191576  -0.808425   1.0

The derivative of the projection with respect to a uniform increase of the weights of the diagonal entries is:

dX
4×4 Matrix{Float64}:
0.101794    0.228658   -0.0880952  -0.0490976
0.228658    0.0565823   0.158685   -0.0880952
-0.0880952   0.158685    0.0565823   0.228658
-0.0490976  -0.0880952   0.228658    0.101794