add complex psd cone

[araujoms]

Closes #290

In the end I didn't use Hypatia's indexing, but something that was close to what SCS was already doing in the real case.

You can test it via the following code:

#include "scs.h"    /* SCS API */
#include <stdio.h>  /* printf */
#include <stdlib.h> /* memory allocation */
#include <math.h>   // sqrt

//C = {{1, 2 + 3 I, 4 + 5 I}, {2 - 3 I, 6, 7 - 8 I}, {4 - 5 I, 7 + 8 I, 9}}, mineig -5.22892944

/* Set up and solve basic qp */
int main(int argc, char **argv) {
  /* Set up the problem data */
  /* A and P must be in compressed sparse column format */
  double A_x[12] = {1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0}; //nzval
  int A_i[12] = {0,  1,  2,  3, 4, 5, 0, 6, 7, 8, 0, 9}; //rowval
  int A_p[10] = {0, 2, 3, 4, 5, 6, 8, 9, 10, 12}; //colptr
  double b[10] = {1.0, 0., 0., 0., 0., 0., 0., 0., 0., 0.};
  double c[9] = {1, 2*sqrt(2), 3*sqrt(2), 4*sqrt(2), 5*sqrt(2), 6, 7*sqrt(2), -8*sqrt(2), 9};
  /* data shapes */
  int n = 9; /* number of variables */
  int m = 10; /* number of constraints */
  int cs[] = {3};
  int cssize = 1;
  /* Allocate SCS structs */
  
  ScsCone *k = (ScsCone *)calloc(1, sizeof(ScsCone));
  ScsData *d = (ScsData *)calloc(1, sizeof(ScsData));
  ScsSettings *stgs = (ScsSettings *)calloc(1, sizeof(ScsSettings));
  ScsSolution *sol = (ScsSolution *)calloc(1, sizeof(ScsSolution));
  ScsInfo *info = (ScsInfo *)calloc(1, sizeof(ScsInfo));

  /* Fill in data struct */
  d->m = m;
  d->n = n;
  d->b = b;
  d->c = c;
  d->A = &(ScsMatrix){A_x, A_i, A_p, m, n};
  d->P = SCS_NULL; //&(ScsMatrix){P_x, P_i, P_p, n, n};

  /* Cone */
  k->z = 1;
  k->cs = cs;
  k->cssize = cssize;

  /* Utility to set default settings */
  scs_set_default_settings(stgs);

  /* Modify tolerances */
  stgs->eps_abs = 1e-9;
  stgs->eps_rel = 1e-9;

  /* Initialize SCS workspace */
  ScsWork *scs_work = scs_init(d, k, stgs);

  /* Solve! */
  int exitflag = scs_solve(scs_work, sol, info, 0);

  /*
   * If we wanted to solve many related problems with different
   * b / c vectors we could update the SCS workspace as follows:
   *
   * int success = scs_update(scs_work, new_b, new_c)
   * int new_exitflag = scs_solve(scs_work, sol, info, 1);
   *
   */

  /* Free SCS workspace */
  scs_finish(scs_work);

  /* Verify that SCS solved the problem */
  printf("SCS solved successfully: %i\n", exitflag == SCS_SOLVED);

  /* Print some info about the solve */
  printf("SCS took %i iters, using the %s linear solver.\n", info->iter,
         info->lin_sys_solver);

  /* Print solution x */
  printf("Optimal solution vector x*:\n");
  for (int i = 0; i < n; ++i) {
    printf("x[%i] = %4f\n", i, sol->x[i]);
  }

  /* Print dual solution y */
  printf("Optimal dual vector y*:\n");
  for (int i = 0; i < m; ++i) {
    printf("y[%i] = %4f\n", i, sol->y[i]);
  }

  /* Free allocated memory */
  free(k);
  free(d);
  free(stgs);
  free(info);
  /* SCS allocates sol->x,y,s if NULL on entry, need to be freed */
  free(sol->x);
  free(sol->y);
  free(sol->s);
  free(sol);
  return 0; /* returning exitflag will set bash exit code to 1 */
}

I didn't do benchmarking because that's just painful in C. Once it's merged I can write the Julia interface and benchmark it.

[kalmarek]

@araujoms you can also benchmark this before merging by compiling and using SCS_jll locally. If you need any pointers on how to do this let me know

[araujoms]

The problem is that SCS.jl can't access the new cone, neither through JuMP nor through the low-level interface. That's the part that I still need to write.

[araujoms]

How do I use the binaries I compiled locally, though? I tried simply replacing the files libscsdir.so and libscsindir.so in .julia/artifacts and all I got were errors.

[kalmarek]

For sure you need to write the interface code in julia, we'll also need to rethink how we interact with the solver in our C-wrapper. @blegat can maybe help you with advertising the complex functionality through MOI apis ;)

To compile jll locally you can either override the library manually, or compile SCS_jll through BinaryBuilder and deploy locally. Both are described in details here:
https://docs.binarybuilder.org/stable/building/#Building-and-testing-JLL-packages-locally

I don't actually recommend the first option (due to compiler options and linking expected by SCS.jl). To run the latter you can:

1. direct Yggdrasil/S/SCS/build_tarballs.jl to use your modified git tree of scs
2. initialize environment in Yggdrasil/ with just BinaryBuilder added
3. inside Yggdrasil/S/SCS folder run

julia --project=@. build_tarballs.jl --debug --verbose x86_64-linux-gnu --deploy=local

4. ask SCS.jl to use deved SCS_jll (created by --deploy=local)

[araujoms]

Thanks! I did it manually in the end because life is too short to deal with BinaryBuilder. Your instructions were essential anyway so I could find out the compilations options that needed to be changed. (Btw since you have mastered the dark arts, perhaps you're interested in ressurrecting SDPAFamily? Someone needs to move it from BinaryProvider to BinaryBuilder).

In any case, I hacked together a bare Julia interface and did some benchmarking, via a simple SDP that computes the minimal eigenvalue of a random complex matrix. The blue line does it with the real psd cone, via the usual mapping C -> [real(C) imag(C); -imag(C) real(C)], and the orange line does it with the complex psd cone.

The speedup is roughly a factor of 5 on average.

[kalmarek]

Thanks! I did it manually in the end because life is too short to deal with BinaryBuilder. Your instructions were essential anyway so I could find out the compilations options that needed to be changed. (Btw since you have mastered the dark arts, perhaps you're interested in ressurrecting SDPAFamily? Someone needs to move it from BinaryProvider to BinaryBuilder).

By no means I'm interested in hacking on some old C/C++ code to get it compiling :D
Similarly to Eric, I'm no longer working in academia, I'm doing this in my spare time for fun.

In any case, I hacked together a bare Julia interface and did some benchmarking

looking forward to it!

The speedup is roughly a factor of 5 on average.

Looks great to me 🚀. What are the axes labels?

[araujoms]

Maybe there are some people who find C++ archaeology fun 🤷 Personally I find it about as fun as warts.

Vertical axis is the time in milliseconds, horizontal axis is the dimension of the matrix.

[kalmarek]

maybe some llms will do it some day for us :D
Can you also try something larger (2^{9,10,11}) to see how it scales there?

[araujoms]

My patience only lasted until 2^9. The ratio of the time taken with the real psd cone over the complex psd cone was [4.3, 2.3, 4.3, 6.3, 7.7] for d = [2^5, 2^6, 2^7, 2^8, 2^9]. The advantage should be constant, and the factor upper bounded by 8 (because we're diagonalizing a d x d complex matrix instead of 2d x 2d real matrix). If you want to know more numbers you'll have to run the code yourself. My C wrapper is here, and the code for benchmarking is below:

using SCS
using SparseArrays
using LinearAlgebra
import Random

function extract_lower_triangle(C::AbstractMatrix{T}) where {T}
    R = real(T)
    is_complex = T <: Complex
    d = size(C,1)
    vec_dim = is_complex ? d^2 : div(d*(d+1),2)
    c = Vector{R}(undef, vec_dim)
    if is_complex
        extract_lower_triangle_complex!(c, C, sqrt(R(2)))
    else
        extract_lower_triangle_real!(c, C, sqrt(R(2)))
    end
    return c
end

function extract_lower_triangle_real!(c, C, sqrt2)
    d = size(C,1)
    counter = 0
    for j = 1:d
        counter += 1
        c[counter] = C[j,j]
        for i = j+1:d
            counter += 1
            c[counter] = C[i,j]*sqrt2
        end
    end
end

function extract_lower_triangle_complex!(c, C, sqrt2)
    d = size(C,1)
    counter = 0
    for j = 1:d
        counter += 1
        c[counter] = real(C[j,j])
        for i = j+1:d
            counter += 1
            c[counter] = real(C[i,j])*sqrt2
            counter += 1
            c[counter] = imag(C[i,j])*sqrt2
        end
    end
end

function benchmark(::Type{T}, d) where {T}
    Random.seed!(1337)
    C = Hermitian(randn(complex(T), d, d))
    if T <: Real
        real_C = Hermitian([real(C) imag(C); -imag(C) real(C)])
        return mineig_raw(real_C)
    else
        return mineig_raw(C)
    end
end

function mineig_raw(C::AbstractMatrix{T}) where {T}
    is_complex = T <: Complex
    R = real(T)
    d = size(C, 1)
    c = extract_lower_triangle(C)
    vec_dim = is_complex ? d^2 : div(d*(d+1), 2)
    id_vec = sparse(extract_lower_triangle(T.(I(d))))
    A = [id_vec'; -sparse(I, vec_dim, vec_dim)]
    b = zeros(R, vec_dim + 1)
    b[1] = 1

    n = vec_dim
    m = vec_dim + 1

    z = 1
    bu = Float64[]
    bl = Float64[]
    l = 0
    ep = 0
    ed = 0
    q = Int[]
    if is_complex
        s = Int[]
        cs = [d]
    else
        s = [d]
        cs = Int[]
    end
    
    P = spzeros(n, n)
    solver = SCS.DirectSolver
    sol = scs_solve(solver, m, n, A, P, b, c, z, l, bu, bl, q, s, cs, ep, ed, Float64[], Float64[])
    display(minimum(eigvals(Hermitian(C))))
    display(sol.info.pobj)
    return sol.info.solve_time
end

run benchmark(Float64,d) to time a d x d complex matrix over the real psd cone, and benchmark(ComplexF64,d) to time the same matrix over the complex psd cone.

[kalmarek]

looks great!
I'm happy to see the very minimal changes to c_wrapper.jl ;)

[bodono]

Awesome, thanks for doing this! I will take a look.