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.3.1 - 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: 5
	  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  3.59e-01  1.39e+00  7.37e-01  1.00e-01  1.97e-04
    75| 7.16e-11  6.49e-11  3.13e-10  2.79e-01  1.00e-01  7.74e-04
------------------------------------------------------------------
status:  solved
timings: total: 7.76e-04s = setup: 8.87e-05s + solve: 6.87e-04s
	 lin-sys: 3.91e-05s, cones: 4.22e-04s, accel: 1.04e-04s
------------------------------------------------------------------
objective = 0.278563
------------------------------------------------------------------

The projection of A is:

X

3×3 Matrix{Float64}:
 1.0       0.76069  0.157298
 0.76069   1.0      0.76069
 0.157298  0.76069  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}:
 -4.4511e-19   -7.53425e-19  -4.36151e-21
 -7.53425e-19  -7.27453e-21   7.73997e-19
 -4.36151e-21   7.73997e-19   4.2487e-19

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.3.1 - 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: 5
	  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.94e+00  3.56e-01  1.20e+01  8.68e+00  1.00e-01  1.59e-04
    75| 2.55e-11  3.48e-12  5.15e-12  4.55e+00  1.00e-01  9.00e-04
------------------------------------------------------------------
status:  solved
timings: total: 9.02e-04s = setup: 8.97e-05s + solve: 8.12e-04s
	 lin-sys: 5.79e-05s, cones: 5.46e-04s, accel: 1.16e-04s
------------------------------------------------------------------
objective = 4.552800
------------------------------------------------------------------

The projection of A is:

X

4×4 Matrix{Float64}:
  1.0       -0.808412   0.191588   0.106775
 -0.808412   1.0       -0.656233   0.191588
  0.191588  -0.656233   1.0       -0.808412
  0.106775   0.191588  -0.808412   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}:
 -6.80725e-8  0.316747     0.316747    0.157095
  0.316747    5.33931e-8   0.630911    0.316747
  0.316747    0.630911    -5.60168e-8  0.316747
  0.157095    0.316747     0.316747    6.4224e-8

This page was generated using Literate.jl.