oxfordcontrol / cosmo.jl Goto Github PK

Just an issue to keep track of the JuMP / MOI integration:

• I am working on it on the mg/MOptInterface branch. The wrapper code is in src/MOIWrapper.jland unit tests in test/UnitTests/jump.jl. I also added the PsdConeTriangle type in src/convexset.hlthat only considers the upper triangle of the decision matrix in question.
• The unit tests pass for simple LPs, QPs, SOCPs and SPDs (with PSDSquare constraint)

To Do:

• support MOI.PositiveSemidefiniteConeTriangle constraint

• make solver pass the JuMP unit tests

Follow up on the slack discussion!

Small suggestion...

I believe the code is mature enough to be structured in a module so that it can be cloned and used as any other package. It would be just a matter of:

• Create the basic structure using PkgDev.generate("QOCS", "ASL")
• Change the repository name to QOCS.jl

No need to distribute it to METADATA.jl yet. The package creation would just make your life easier and add automatically the right scripts for the unittests and coverage.

We want to test CG on problems that are quite big, as this it only makes sense to use CG in those cases.

We would also want to apply to a challenging, ill conditioned dataset that will allow us to compare e.g. preconditioners. #30 is perhaps suitable.

Might actually improve convergence since symmetry is implicitly enforced! (not true since LHS is always symmetric by definition)

Is currently implemented to be a matrix of size nxn if x=vec(X) \in R^n

This is issue is to track the testing of a block-lanczos projection into the SemiDefinite Cone.

ToDo List:

• Code Block-Lanczos for QOCS that works on a simple example
• Check its performance on a medium-sized dataset
• Make it efficient
• Try its performance on a bigger dataset

The tag name "v0.4" is not of the appropriate SemVer form (vX.Y.Z).
cc: @migarstka

Should the array K.s in the cone contain the size of the matrix n or the number of elements of s (n^2)? Investigate what causes less overhead

I believe that JLD is deprecated and is not compatible with Julia 1.0. @migarstka should I create a PR for switching to JLD2?

• Show solver solves 4-5 problems from SDPLib.
• Compare the results to the results obtained by MOSEK.
• Check if infeasible problems are detected. (see #5 )

Some things to check:

• difference in optimal value
• difference in hinf of optimal decision variable
• solve time
• steps to convergence
• parameter changes
• comparison of OSSDP and SCS against MOSEK (if both OSSDP and SDP are off more than 1e-3 from MOSEK solution, which one is worse?)

@goulart-paul pointed me to COSMO. I was wondering if there any plans to combine this package with your osqp package (which I use and love a lot) or to release Python bindings? I am currently working on a trust region optimisation problem and so can't use osqp.

• understand fixed point acceleration

• Settings
• Results

Add ability to warm start problems

Todo:

• Separately test all logical parts of the code
• Plot the change of values of x,s,z,mu,lambda over the iterations
• Check problem formulation, transformation of the matrices to vectors [correct]
• Run several tests with different combinations of alpha, sigma and rho and determine which one leads to the lowest deviation from the true solution. [done]
• Check if you get residual convergence r^k -> 0, objective convergence, dual variable convergence
• Check a larger variety of example problems.
• Try Version without auxiliary variable z, leads to 0 in left side of linear system.
• Compare to Yang's code

Further Debugging suggestions:

• Solve a 1 or 2 variable QP or a 2 variable LP turned into an SDP
• Check if solutions fulfill the KKT eq

Observations:

• \mu is not equal to \nu as it should be by definition of the algorithm

• Objective and residuals converge

• Solution is not pos semidef (has eigenvalue at -6)

• Objective value

  • is smaller than true solution
  • is close to true solution if sigma is set high

• Constraint Ax=z not exactly satisfied (primal residual about 0.54, but not moving)

• Constraint x=s not satisfied (only if sigma is set very high, but then z=b constraint less satisfied)

• Constraint z=b satisfied

• Constraint s in pos semidef cone satisfied

• Dual variable \mu=[0,0] (relating to z=b constraint)

• Cost with s: q's=25.6, q'x = -38.2, true value = 13.9 (from SCS)

Follow up on the discussion in the one_big_module pull request.

I think we should use Julia 1.0.1's mul! functions for left- and right multiplication with a diagonal matrix.
We probably have to add mul! method definitions for the identity matrix (i.e. no scaling case)
mul!(L::IdentityMatrix, M::AbstractMatrix).

Is there a native Julia function for left-right-multiplication? If not we should keep the lrmul!() definitions.

The problem

I am trying to benchmark some things in SDPLib and it seems that QOCS is producing incorrect results. For example, in problem truss1 QOCS gives dual infeasibility while SCS (called via JuMP 0.18.3) returns the correct result: primal and dual feasible of optimal objective ~= -9.

I will try to understand why this is happening. @migarstka do you have any thoughts?

How to reproduce

• unzip and save truss1.jld.zip in test/Benchmarks/sdplib/truss1.jld. You might have to create the folder sdplib.
• run sampleProblem.jl (as seen on branch Lanczos).

Alternatively you can create the whole SDPLib dataset by running sdpLibImport.jl. This takes some time.

Info

• QOCS: 0b3a6ba (devel).
• Julia: 1.0
• JuMP: 0.18.3
• SCS: 0.4.0

https://kul-forbes.github.io/scs/page_benchmark_results.html

The tag name "0.4.0" is not of the appropriate SemVer form (vX.Y.Z).
cc: @migarstka

In order to support the box constraint with infeasibility detection and scaling, we have to implement some of the functions in src/convexset.jl

• in_dual(): Right now this seems only to check whether a vector is outside the box
• in_recc(): Evaluates always to true
• scale!(): Multiplies the box boundaries with a SplitVector. What was the rational behind that @goulart-paul

Once the code works, decompose the code into different files and modules in a logical way

• Add code files to a folder /src and all tests to /test
• Move all functionality in the main loop into separate functions and possibly also separate files
• Double check and improve naming of variables functions etc, especially the !-convention

In order to implement the psd projection of a matrix in packed storage format, i.e. upper triangular part of the matrix stacked in a vector, we should make use of the corresponding LAPACK function. Otherwise, we will have to "unpack" the matrix every time we want to project, check infeasibility, etc.

http://www.icl.utk.edu/~mgates3/docs/lapack.html
gives a good overview on the available LAPACK functions for eigenvalue computation. It seems like the relevant functions are spev, spevx and spevd.

Unfortunately, the julia LAPACK package doesn't provide wrappers for any of these functions. However, the wrapper for e.g. syevr doesn't seem to be too complicated:

LAPACK.jl Documentation - Line 4949

So I think we should write our own wrapper for one of the above functions.

Added in d4789f8

Wrong implementation of relaxation seems to be the major bug causing the convergence problems!

• Double-check that the relaxation of x and s is correctly implemented. Current implementation for x and s:
• Relaxation of s shouldnt be a problem concerning pos sem-definiteness, since if A, B pos semdef => (A+B) pos semdef

      xt from solving linear system
      xNew = α*xt + (1-α)*xPrev
      st from projection onto pos def cone
      sNew = α*st + (1-α)*sPrev

• Symmetric matrix equilibration
• Ruiz equilibration
• Scale termination criteria
• Scale Infeasibility checks
• Note: One has to scale the variables of the psd and soc cones with one scalar value, otherwise one destroys symmetry

Next steps:

• Investigate why my rho and sigma lead to different convergence behavior

• Choose rho intelligently in the beginning [different since I only have equality constraints]

• Check that scaling of OSQP and OSSDP standard scaling is the same for several LP, QPs [seems to be the case]

• Check if scaling helps with SOCP and SDP problems

• Implement scaling similar to CDCS and compare matrices.

• Run SCS, OSSDP and OSQP for several problems (LP, QP, SOCPs SDPs) with and without scaling and with Ruiz and SCS scaling strategies.

Bug 1 [Solved]:

• Scaling by e.g. 0.5 leads to the correct solution although it needs way more iterations

• Scaling by the Ruiz equilibration leads to a different solution for some reason, Ax=b fulfilled, but x=s in SDP cone not fulfilled at (what the solver considers a) "optimal" solution

• Scaling back the solution seems to move solution vector x and s out of the psd-cone. Why?
  The solution was to scale the psd-cone projection as well

Bug2 [solved]:

Scaling doesn't seem to improve convergence by much. Possible error sources:

• Wrong test problem (already pretty well scaled, but why does SCS need way fewer iterations)
• Wrong implementation of Ruiz equilibration algorithm (or the additional cost normalization)
• Wrong equilibration Ansatz
• Wrong / different residual formulation or scaling or residuals wrong
  The current scaling matrix destroyed the symmetry property of X

Solved in 5e95364

I am getting unexpectedly large memory allocations when running some problems of the Maros-Mescazos dataset.

For example, when solving PRIMALC2 (2076 non-zeros) @time reports:

Running QOCS:  7.148138 seconds (8.61 M allocations: 6.332 GiB, 7.95% gc time)

Or, when running QSCTAP2 (11417 non-zeros) I get:

Running QOCS:  1226.066138 seconds (93.82 M allocations: 1.013 TiB, 61.93% gc time)

Note that the memory allocation for LDL is very small (hundreds of KiB for PRIMALC2)

You can reproduce the results by running this script. Just change the for loop in order to run only the particular problem, e.g.:

for file in ["PRIMALC2"]

The tag name "v0.3" is not of the appropriate SemVer form (vX.Y.Z).
cc: @migarstka

This is to prepare Maros-Meszaros dataset to be used in #29

Look at different SEDUMI to find a good way to let the user describe his SDP problem, e.g.

• number of free, cone variables, etc.

"The rows of the data matrix A correspond to the cones in K. The rows of A must be in the order of the cones given above, i.e., first come the rows that correspond to the zero/free cones, then those that correspond to the positive orthants, then SOCs, etc. "
Implemented via merging cones branch into devel: e89ae5f

Find the correct residuals for the SDP case to have a working termination criterion.
Currently, the following is used:

Problem:

min 1/2x'Px+q'x
s.t. Ax=z, z=b, x=s and s in S+

r_prim = Ax - z
r_dual = Px+q+A'mu+\lambda

Argument for feasibility of other primal variables: z = b, x = s is ensured by projections and therefore not necessarily considered in the primal residual

Strictly following the Boyd ADMM paper the dual residual should be:

r_dual = sigma * (s_k+1 - s_k)

• remove unnecessary and outdated tests
• update relevant tests to new solver format

It would be nice to adopt uniform coding guidelines (e.g. this or something else).

I believe that only storing the upper triangular elements of PSD variables is the way to go.

This would result in a smaller linear system. Furthermore, converting a vector that stores the upper triangular elements of a symmetric matrix to a symmetric matrix is faster than performing (A + A')/2, thus I don't see any reason not to go for it. Indeed running:

using LinearAlgebra
using BenchmarkTools

function vector_to_symmatrix!(A::Array{Float64, 2}, a::Array{Float64, 1})
    # Indexing as seen on https://github.com/JuliaLang/julia/blob/870f99568aab0847962e453abb1c4b29d6906fe7/stdlib/LinearAlgebra/src/triangular.jl#L250-L252
    n = size(A, 1)
    @inbounds @simd for j in 1:n
        @inbounds @simd for i in 1:j
            idx = Int(j*(j-1)/2) + i
            A[i,j] = a[idx]
        end
    end
    # The matrix would passed to the eigensolver as Symmetric(A)
end

function make_symmetric_loop!(A::Array{Float64, 2})
    n = size(A, 1)
    @inbounds @simd for j in 1:n
        @inbounds @simd for i in 1:j
            A[j,i] = (A[j,i] + A[i,j])/2
        end
    end
end

n = 1000;

# Convert a vector to symmetric matrix
a = randn(Int(n*(n+1)/2));
A = zeros(n, n)
@btime vector_to_symmatrix!(A, a)

# Try two different versions of (A + A')/2
A = randn(n, n)
@btime @. A = (A + A')/2
@btime make_symmetric_loop!(A)

gives:

 578.334 μs (1 allocation: 32 bytes)
  4.528 ms (7 allocations: 7.63 MiB)
  1.393 ms (0 allocations: 0 bytes)

• Implement an infeasibility detection routine.
• Figure out a way to do Goran's check without the need for another eigenvalue decomposition
• Only check infeasibility every now and then (also check the conditions in order of computation cost, starting with the cheapest)

Where can I obtain the .MAT files required for the Maros and Mecrazos tests? Using these did not work.