Add Arbitrary Precision support for PSD constrains (SDPs)

[Mark-L-Stone]

Per https://oxfordcontrol.github.io/COSMO.jl/stable/literate/build/arbitrary_precision/ , "Since we call the LAPACK function syevr for eigenvalue decompositions, we currently only support solving problems with PSD constraints in Float32 and Float64."

I would like a COSMO version which can support SDP cone (PSD constraints) in arbitrary precision, in order to successfully solve extremely ill-conditioned SDP problems. Hopefully, there would not be a weak chain in the link preventing SDPs from being solved to arbitrary accuracy, given enough computing time and memory.

[migarstka]

The main obstacle is that we need an eigenvalue decomposition method for BigFloat types.
It seems like https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl could provide this.

[lrnv]

An old topic on julia base indeed links to GenericLinearAlgebra, see there: JuliaLang/LinearAlgebra.jl#61

What exactly do we need here ? Looks like Cholesky, LU and QR are currently already implemented, maybe this is enough ?

[migarstka]

I am fairly sure that we just need to add a projection method for PsdCone and PsdConeTriangle for precision types besides Float32, Float64.

COSMO.jl/src/convexset.jl

Lines 186 to 204 in bedbed0

	function _project!(X::AbstractMatrix, ws::PsdBlasWorkspace{T}) where{T}
	#computes the upper triangular part of the projection of X onto the PSD cone
	#allocate additional workspace arrays if the ws
	#work and iwork have not yet been sized
	if ws.lwork == -1
	_syevr!(X,ws)
	ws.lwork = BLAS.BlasInt(real(ws.work[1]))
	resize!(ws.work, ws.lwork)
	ws.liwork = ws.iwork[1]
	resize!(ws.iwork, ws.liwork)
	end
	# below LAPACK function does the following: w,Z = eigen!(Symmetric(X))
	_syevr!(X, ws)
	# compute upper triangle of: X .= Z*Diagonal(max.(w, 0.0))*Z'
	rank_k_update!(X, ws)
	end

@lrnv Would you be interested in implementing this? I am currently in the process of finishing my thesis, but might be able to look at this in June.

[lrnv]

@migarstka with good guidance, yes I can try for sure.

[migarstka]

I would just make up a small SDP with BigFloat precision data and follow the error messages.

You can take the following:

using COSMO, LinearAlgebra

# Small example SDP

T = Float64
# T = BigFloat

# Problem data
C =  T[1. 2; 2. 2]
A = T[1. 5.; 5 2]
b = T(4);


# define the cost function
P = zeros(T, 4, 4)
q = vec(C)
# define the constraints
# A x = b
cs1 = COSMO.Constraint(vec(A)', -b, COSMO.ZeroSet)
# X in PSD cone
cs2 = COSMO.Constraint(Matrix{T}(1.0I, 4, 4), zeros(T, 4), COSMO.PsdCone)
constraints = [cs1; cs2]

# assemble and solve the model
model = COSMO.Model{T}();
settings = COSMO.Settings{T}(decompose = false, verbose = true, accelerator = EmptyAccelerator)
assemble!(model, P, q, constraints, settings = settings)
result = COSMO.optimize!(model);

X_sol = reshape(result.x, 2, 2)
obj_value = result.obj_val


# solve the same problem but with PsdConeTriangle (scale off-diagonals by √2)
# The solution is stored in x = [X11; sqrt(2) X12; X22]
At = T[1. sqrt(2)*5 2.]
Ct =  T[1.;2*sqrt(2); 2]
Pt = zeros(T, 3, 3)
qt = vec(Ct)
cs1t = COSMO.Constraint(At, -b, COSMO.ZeroSet)
cs2t = COSMO.Constraint(Matrix{T}(1.0I, 3, 3), zeros(T, 3), COSMO.PsdConeTriangle)
constraints2 = [cs1t; cs2t]
model2 = COSMO.Model{T}();
assemble!(model2, Pt, qt, constraints2, settings = settings)
result2 = COSMO.optimize!(model2);

# recover the solution in matrix form by unscaling the off-diagonals by 1 / sqrt(2)
X_sol2 = zeros(2, 2)
COSMO.populate_upper_triangle!(X_sol2, result2.x, 1 / sqrt(2))
X_sol2 = Symmetric(X_sol2)

This script solves an SDP once using the PsdCone constraint type and once with the (more compact and faster) PsdConeTriangle constraint type.
Changing T = BigFloat gives the following error:

ERROR: MethodError: no method matching _syevr!(::Base.ReshapedArray{BigFloat,2,SubArray{BigFloat,1,Array{BigFloat,1},Tuple{UnitRange{Int64}},true},Tuple{}}, ::COSMO.PsdBlasWorkspace{BigFloat})
Stacktrace:
 [1] _project!(::Base.ReshapedArray{BigFloat,2,SubArray{BigFloat,1,Array{BigFloat,1},Tuple{UnitRange{Int64}},true},Tuple{}}, ::COSMO.PsdBlasWorkspace{BigFloat}) at /Users/Micha/Dropbox/Research/OSSDP/Code/src/convexset.jl:0
 [2] project!(::SubArray{BigFloat,1,Array{BigFloat,1},Tuple{UnitRange{Int64}},true}, ::COSMO.PsdCone{BigFloat}) at /Users/Micha/Dropbox/Research/OSSDP/Code/src/convexset.jl:275
 [3] project!(::COSMO.SplitVector{BigFloat}, ::COSMO.CompositeConvexSet{BigFloat}) at /Users/Micha/Dropbox/Research/OSSDP/Code/src/convexset.jl:816
 [4] macro expansion at ./timing.jl:233 [inlined]
 [5] admm_z!(::COSMO.SplitVector{BigFloat}, ::Array{BigFloat,1}, ::COSMO.CompositeConvexSet{BigFloat}, ::Int64) at /Users/Micha/Dropbox/Research/OSSDP/Code/src/solver.jl:15
 [6] optimize!(::COSMO.Workspace{BigFloat}) at /Users/Micha/Dropbox/Research/OSSDP/Code/src/solver.jl:153
 [7] top-level scope at REPL[128]:1

As predicted the error happens because _sveyr! can't handle BigFloat.
To get this error message, you have to comment-out the type-check here:

COSMO.jl/src/interface.jl

Line 391 in a138a5f

type_checks(convex_set::Union{PsdCone{BigFloat}, PsdConeTriangle{BigFloat}}) = throw(ArgumentError("COSMO currently does not support the combination of PSD constraints and BigFloat."))