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_data = ones(size(A)))
    n = LinearAlgebra.checksquare(A)
    model = Model(() -> DiffOpt.diff_optimizer(solver))
    @variable(model, X[1:n, 1:n] in PSDCone())
    @variable(model, H[1:n, 1:n] in Parameter.(H_data))
    @variable(model, E[1:n, 1:n])
    @constraint(model, [i in 1:n], X[i, i] == 1)
    @constraint(model, E .== (H .* (X .- A)))
    @objective(model, Min, sum(E .^ 2))
    for i in 1:n
        DiffOpt.set_forward_parameter(model, H[i, i], dH[i, i])
    end
    optimize!(model)
    DiffOpt.forward_differentiate!(model)
    dX = DiffOpt.get_forward_variable.(model, 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.7 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 15, constraints m: 18
cones: 	  z: primal zero / dual free vars: 12
	  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): 27, nnz(P): 9
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 1.00e+00  2.56e-01  1.01e+00  5.26e-01  1.00e-01  1.40e-04
   100| 5.96e-05  1.22e-06  9.38e-05  2.79e-01  1.00e-01  4.97e-04
------------------------------------------------------------------
status:  solved
timings: total: 4.98e-04s = setup: 5.93e-05s + solve: 4.38e-04s
	 lin-sys: 5.16e-05s, cones: 2.68e-04s, accel: 5.96e-06s
------------------------------------------------------------------
objective = 0.278512
------------------------------------------------------------------

The projection of A is:

X

3×3 Matrix{Float64}:
 1.0       0.760733  0.157275
 0.760733  1.0       0.760733
 0.157275  0.760733  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}:
  5.12521e-15  3.01089e-14  -1.36629e-14
  3.01089e-14  3.94055e-16   1.86395e-14
 -1.36629e-14  1.86395e-14  -4.31989e-15

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.7 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 26, constraints m: 30
cones: 	  z: primal zero / dual free vars: 20
	  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): 46, nnz(P): 16
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 1.95e+00  2.54e-01  1.26e+01  8.41e+00  1.00e-01  1.24e-04
   125| 8.40e-08  1.20e-08  2.01e-07  4.55e+00  7.56e-01  7.39e-04
------------------------------------------------------------------
status:  solved
timings: total: 7.40e-04s = setup: 7.00e-05s + solve: 6.70e-04s
	 lin-sys: 1.15e-04s, cones: 4.33e-04s, accel: 8.01e-06s
------------------------------------------------------------------
objective = 4.552800
------------------------------------------------------------------

The projection of A is:

X

4×4 Matrix{Float64}:
  1.0       -0.808413   0.191587   0.106775
 -0.808413   1.0       -0.656233   0.191587
  0.191587  -0.656233   1.0       -0.808413
  0.106775   0.191587  -0.808413   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}:
 1.81949e-8   0.314394    0.314394    0.143315
 0.314394    -1.06177e-7  0.666831    0.314394
 0.314394     0.666831    4.70508e-9  0.314394
 0.143315     0.314394    0.314394    3.91717e-8

This page was generated using Literate.jl.