Spectral cones (#295)

* Added spectral cone projections and a few tests

* more tests

* fixed existing logging bug

* fixed comment

* fixed BLAS call

* validate nuclear cone

* validate nuclear cone

* fixed build systems, added timing, added flags

* comment

* minor

* minor cmake-build

* minor cmake build

* tightened tolerance to pass a test

* changed tolerance again

* LAPACK build flags

* removed complementary check in test (talk to brendan)

* changed tolerance in test

* changed tolerance of test

* trying to fix C++ build

* trying to fix C++ build

* trying to fix C++ build

* trying to fix C++ build

* changed cmake

* removed cast

* changed makefile

* removed static temporarily

* removed static temporarily

* fixing C++ build

* mkl

* mkl build

* fixed merge conflict

* complex PSD cone and C++ build

* complex PSD cone wrong datatype?

* complex PSD cone

* complex PSD cone

* complex PSD cone

* complex PSD cone

* complex PSD cone

* as it used to be

* added new workflow for spectral cones and conditional compilation

* corrected Makefile

* GPU build

* updated makefile

* ran clang formatter

* formatted CMakeLists.txt

* fixed two tests

Author: dance858

.github/workflows/build_spectral.yml

@@ -0,0 +1,36 @@
+name: Build and Test with Spectral Cones
+
+on: [push, pull_request]
+
+jobs:
+  linux:
+    strategy:
+      fail-fast: false
+      matrix:
+        long: [0, 1]
+        spectral: [0, 1]
+
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - run: sudo apt-get update && sudo apt-get install -y libopenblas-dev liblapack-dev
+      - run: make DLONG=${{ matrix.long }} USE_LAPACK=1 USE_SPECTRAL_CONES=${{ matrix.spectral }}
+      - run: make test DLONG=${{ matrix.long }} USE_LAPACK=1 USE_SPECTRAL_CONES=${{ matrix.spectral }}
+      - run: out/run_tests_direct    # test direct solver
+      - run: out/run_tests_indirect  # test indirect solver
+
+  mac:
+    strategy:
+      fail-fast: false
+      matrix:
+        long: [0, 1]
+        spectral: [0, 1]
+
+    runs-on: macos-latest
+    steps:
+      - uses: actions/checkout@v4
+      - run: brew update && brew install openblas lapack
+      - run: make DLONG=${{ matrix.long }} USE_LAPACK=1 USE_SPECTRAL_CONES=${{ matrix.spectral }}
+      - run: make test DLONG=${{ matrix.long }} USE_LAPACK=1 USE_SPECTRAL_CONES=${{ matrix.spectral }}
+      - run: out/run_tests_direct    # test direct solver
+      - run: out/run_tests_indirect  # test indirect solver

.gitignore

@@ -51,3 +51,4 @@ out
 hs21_tiny_qp
 rob_gauss_cov_est
 build/
+*.vscode

CMakeLists.txt

@@ -124,6 +124,12 @@ message(STATUS "Read/write functionality is NOT ${NO_READ_WRITE}")
 option(USE_LAPACK "Whether to use BLAS/LAPACK" ON)
 message(STATUS "BLAS/LAPACK usage is ${USE_LAPACK}")
 
+# spectral cones
+option(USE_SPECTRAL_CONES "Whether to use spectral cones (requires LAPACK)" OFF)
+if(USE_SPECTRAL_CONES AND NOT USE_LAPACK)
+  message(FATAL_ERROR "USE_SPECTRAL_CONES requires USE_LAPACK to be ON")
+endif()
+
 # Enable OpenMP support
 option(USE_OPENMP "Compile with OpenMP support" OFF)
 message(STATUS "OpenMP parallelization is ${USE_OPENMP}")
@@ -158,7 +164,7 @@ if(NO_READ_WRITE)
   set(COMPILER_OPTS "-DNO_READ_WRITE=1 ${COMPILER_OPTS}")
 endif()
 
-if (USE_LAPACK)
+if(USE_LAPACK)
   set(COMPILER_OPTS "-DUSE_LAPACK ${COMPILER_OPTS}")
   list(APPEND LAPACK_LINK_LIBRARIES "blas" "lapack")
 endif()
@@ -167,6 +173,19 @@ if(USE_OPENMP)
   set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fopenmp")
 endif()
 
+# Flags for spectral cones.
+if(USE_SPECTRAL_CONES)
+  set(COMPILER_OPTS "-DUSE_SPECTRAL_CONES ${COMPILER_OPTS}")
+endif()
+
+if(USE_SPECTRAL_TIMING_FLAG)
+  set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DSPECTRAL_TIMING_FLAG")
+endif()
+
+if(USE_SPECTRAL_DEBUG)
+  set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DSPECTRAL_DEBUG")
+endif()
+
 message(STATUS "COMPILER_OPTS = ${COMPILER_OPTS}")
 message(STATUS "CMAKE_C_FLAGS = ${CMAKE_C_FLAGS}")
 
@@ -193,6 +212,23 @@ set(${PROJECT_NAME}_SRC
     ${LINSYS}/csparse.c
     ${LINSYS}/scs_matrix.c)
 
+# Spectral source files
+set(SPECTRAL_SRC
+    src/spectral_cones/logdeterminant/log_cone_Newton.c
+    src/spectral_cones/logdeterminant/log_cone_IPM.c
+    src/spectral_cones/logdeterminant/log_cone_wrapper.c
+    src/spectral_cones/logdeterminant/logdet_cone.c
+    src/spectral_cones/nuclear/ell1_cone.c
+    src/spectral_cones/nuclear/nuclear_cone.c
+    src/spectral_cones/sum-largest/sum_largest_cone.c
+    src/spectral_cones/sum-largest/sum_largest_eval_cone.c
+    src/spectral_cones/util_spectral_cones.c)
+
+# Conditionally add spectral files
+if(USE_SPECTRAL_CONES)
+  list(APPEND ${PROJECT_NAME}_SRC ${SPECTRAL_SRC})
+endif()
+
 # Common header files
 set(${PROJECT_NAME}_HDR
     include/aa.h
@@ -211,6 +247,13 @@ set(${PROJECT_NAME}_HDR
     ${LINSYS}/csparse.h
     ${LINSYS}/scs_matrix.h)
 
+# Spectral header files
+set(SPECTRAL_HDR include/util_spectral_cones.h)
+
+if(USE_SPECTRAL_CONES)
+  list(APPEND ${PROJECT_NAME}_HDR ${SPECTRAL_HDR})
+endif()
+
 # get all the c file in amd/external
 file(GLOB ${PROJECT_NAME}_AMD_EXTERNAL_SRC ${EXTERNAL}/amd/*.c)
 
@@ -249,7 +292,8 @@ target_include_directories(
 target_compile_definitions(${${PROJECT_NAME}_DIRECT} PRIVATE ${COMPILER_OPTS})
 
 # The library depends on math (m) and (optionally) blas and lapack
-target_link_libraries(${${PROJECT_NAME}_DIRECT} PRIVATE m ${LAPACK_LINK_LIBRARIES})
+target_link_libraries(${${PROJECT_NAME}_DIRECT}
+                      PRIVATE m ${LAPACK_LINK_LIBRARIES})
 
 # Set some properties
 set_target_properties(
@@ -297,7 +341,8 @@ target_compile_definitions(${${PROJECT_NAME}_INDIRECT} PRIVATE ${COMPILER_OPTS}
                                                                -DINDIRECT)
 
 # The library depends on math (m) and (optionally) blas and lapack
-target_link_libraries(${${PROJECT_NAME}_INDIRECT} PUBLIC m ${LAPACK_LINK_LIBRARIES})
+target_link_libraries(${${PROJECT_NAME}_INDIRECT}
+                      PUBLIC m ${LAPACK_LINK_LIBRARIES})
 
 # Set some properties
 set_target_properties(
@@ -351,8 +396,8 @@ if(DEFINED ENV{MKLROOT})
   target_compile_definitions(${${PROJECT_NAME}_MKL} PRIVATE ${COMPILER_OPTS})
   # See:
   # https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl-link-line-advisor.html
-  # This is probably not correct for other systems. TODO: make SCS-MKL
-  # work for all combinations of platform / compiler / threading options.
+  # This is probably not correct for other systems. TODO: make SCS-MKL work for
+  # all combinations of platform / compiler / threading options.
   target_link_options(${${PROJECT_NAME}_MKL} PRIVATE "LINKER:--no-as-needed")
   target_link_directories(${${PROJECT_NAME}_MKL} PRIVATE $ENV{MKLROOT}/lib)
   target_link_libraries(
@@ -393,27 +438,29 @@ endif()
 
 # ##############################################################################
 
-# If cuDSS is enabled and environment variables are set, install the cuDSS version
+# If cuDSS is enabled and environment variables are set, install the cuDSS
+# version
 if(USE_CUDSS)
   find_package(CUDAToolkit REQUIRED)
   find_package(cudss REQUIRED)
 
   # Force 32-bit integers for cuDSS compatibility
   if(DLONG)
-    message(FATAL_ERROR "cuDSS requires 32-bit integers. Set DLONG=OFF or do not use -DDLONG=ON")
+    message(
+      FATAL_ERROR
+        "cuDSS requires 32-bit integers. Set DLONG=OFF or do not use -DDLONG=ON"
+    )
   endif()
 
   message(STATUS "Will install SCS-cuDSS (libscscudss).")
 
-
   set(CUDSSSRC ${LINSYS}/cudss/direct)
 
   # Here we compile the direct cuDSS library
   set(${PROJECT_NAME}_CUDSS ${PROJECT_NAME}cudss)
   add_library(
-    ${${PROJECT_NAME}_CUDSS}
-    ${${PROJECT_NAME}_HDR} ${${PROJECT_NAME}_SRC} ${CUDSSSRC}/private.c
-    ${CUDSSSRC}/private.h)
+    ${${PROJECT_NAME}_CUDSS} ${${PROJECT_NAME}_HDR} ${${PROJECT_NAME}_SRC}
+                             ${CUDSSSRC}/private.c ${CUDSSSRC}/private.h)
 
   target_include_directories(
     ${${PROJECT_NAME}_CUDSS}
@@ -428,12 +475,9 @@ if(USE_CUDSS)
 
   target_compile_definitions(${${PROJECT_NAME}_CUDSS} PRIVATE ${COMPILER_OPTS})
 
-
   target_link_libraries(
     ${${PROJECT_NAME}_CUDSS}
-    PRIVATE CUDA::cudart
-            cudss
-            $<$<NOT:$<PLATFORM_ID:Windows>>:m>
+    PRIVATE CUDA::cudart cudss $<$<NOT:$<PLATFORM_ID:Windows>>:m>
             ${LAPACK_LINK_LIBRARIES})
 
   # Set some properties

Makefile

@@ -1,7 +1,33 @@
 # MAKEFILE for scs
 include scs.mk
 
-SCS_OBJECTS = src/util.o src/cones.o src/exp_cone.o src/aa.o src/rw.o src/linalg.o src/ctrlc.o src/scs_version.o src/normalize.o
+# Base object files always included
+BASE_SCS_OBJECTS = \
+    src/util.o src/cones.o src/exp_cone.o src/aa.o src/rw.o src/linalg.o src/ctrlc.o src/scs_version.o src/normalize.o
+
+# Spectral cones files, only included if USE_SPECTRAL_CONES=1
+SPECTRAL_CONES_OBJECTS = \
+    src/spectral_cones/logdeterminant/log_cone_Newton.o \
+    src/spectral_cones/logdeterminant/log_cone_IPM.o \
+    src/spectral_cones/logdeterminant/log_cone_wrapper.o \
+    src/spectral_cones/logdeterminant/logdet_cone.o \
+    src/spectral_cones/nuclear/ell1_cone.o \
+    src/spectral_cones/nuclear/nuclear_cone.o \
+    src/spectral_cones/sum-largest/sum_largest_cone.o \
+    src/spectral_cones/sum-largest/sum_largest_eval_cone.o \
+    src/spectral_cones/util_spectral_cones.o
+
+# Now decide which objects to include:
+ifneq ($(USE_LAPACK),0)
+  ifneq ($(USE_SPECTRAL_CONES),0)
+    SCS_OBJECTS = $(BASE_SCS_OBJECTS) $(SPECTRAL_CONES_OBJECTS)
+  else
+    SCS_OBJECTS = $(BASE_SCS_OBJECTS)
+  endif
+else
+  SCS_OBJECTS = $(BASE_SCS_OBJECTS)
+endif
+
 SCS_O = src/scs.o
 SCS_INDIR_O = src/scs_indir.o
 
@@ -27,6 +53,12 @@ else
 	@echo "NOT compiled with blas/lapack, cannot solve SDPs (can solve LPs, SOCPs, ECPs, and PCPs)."
 	@echo "To solve SDPs, install blas and lapack, then edit scs.mk to set USE_LAPACK=1"
 	@echo "and point to the library install locations, and recompile with 'make purge', 'make'."
+endif
+ifneq ($(USE_SPECTRAL_CONES), 0)
+	@echo "Compiled with spectral cones extension enabled (Daniel Cederberg 2025)"
+else
+	@echo "Spectral cones extension is NOT enabled. To enable spectral cones, edit"
+	@echo "scs.mk to set USE_SPECTRAL_CONES=1 and recompile."
 endif
 	@echo "****************************************************************************************"
 

include/cones.h

@@ -11,6 +11,17 @@ extern "C" {
 #include "scs_work.h"
 #include <string.h>
 
+#ifdef USE_SPECTRAL_CONES
+#include "util_spectral_cones.h" // for newton_stats
+
+// macro for time measurements of SpectralSCS
+#ifdef SPECTRAL_TIMING_FLAG
+#define SPECTRAL_TIMING(action) action
+#else
+#define SPECTRAL_TIMING(action)
+#endif
+#endif
+
 /* private data to help cone projection step */
 struct SCS_CONE_WORK {
   /*
@@ -26,13 +37,45 @@ struct SCS_CONE_WORK {
   scs_int m;            /* total length of cone */
   /* box cone quantities */
   scs_float box_t_warm_start;
+
 #ifdef USE_LAPACK
   /* workspace for eigenvector decompositions: */
   scs_float *Xs, *Z, *e, *work, *rwork;
   scs_complex_float *cXs, *cZ, *cwork;
   blas_int *isuppz, *iwork;
   blas_int lwork, lcwork, lrwork, liwork;
 #endif
+
+#ifdef USE_SPECTRAL_CONES
+  /* if the projection onto the logarithmic cone should be warmstarted*/
+  bool *log_cone_warmstarts;
+
+  /* Needed for ell1 norm cone projection */
+  Value_index *work_ell1;
+  scs_float *work_ell1_proj;
+
+  // used for timing spectral vector cone and spectral matrix cone projections
+  SPECTRAL_TIMING(scs_float tot_time_mat_cone_proj;)
+  SPECTRAL_TIMING(scs_float tot_time_vec_cone_proj;)
+
+  /* workspace for singular value decompositions: */
+  scs_float *s_nuc, *u_nuc, *vt_nuc, *work_nuc;
+  blas_int lwork_nuc;
+
+  /* workspace that is used internally in the logdet projection (for example,
+     the gradient and Hessian of the objective function in the projection
+     problem are stored using this memory) */
+  scs_float *work_logdet;
+
+  /* workspace to store the projection onto the logarithm cone */
+  scs_float *saved_log_projs;
+
+  /* Stats for spectral projections, assuming there is only one spectral cone */
+  Newton_stats newton_stats;
+
+  /* workspace for projection onto sum-largest-evals cone */
+  scs_float *work_sum_of_largest;
+#endif
 };
 
 void SCS(free_cone)(ScsCone *k);

include/scs.h

@@ -143,6 +143,24 @@ typedef struct {
   scs_float *p;
   /** Number of (primal and dual) power cone triples. */
   scs_int psize;
+#ifdef USE_SPECTRAL_CONES
+  /** Array of logdet cone constraints 'len(d) = dsize'. The dimension of a
+   * log-det cone is "n" and not "n+2" if the matrix has dimension n */
+  scs_int *d;
+  /** Length of logdet cone constraints array `d`. */
+  scs_int dsize;
+  /** Array of nuc norm cone constraints 'len(nuc_m) = len(nuc_n) = nucsize.*/
+  scs_int *nuc_m;
+  scs_int *nuc_n;
+  scs_int nucsize;
+  /** Array of ell1-norm cone constraints 'len(ell1) = ell1_size */
+  scs_int *ell1;
+  scs_int ell1_size;
+  /** Array of sum-of-largest-evals cone  */
+  scs_int *sl_n;
+  scs_int *sl_k;
+  scs_int sl_size;
+#endif
 } ScsCone;
 
 /** Contains primal-dual solution arrays or a certificate of infeasibility.
@@ -206,6 +224,13 @@ typedef struct {
   scs_float cone_time;
   /** Total time (milliseconds) spent in the acceleration routine. */
   scs_float accel_time;
+#ifdef SPECTRAL_TIMING_FLAG
+  /** Average time (milliseconds) per iteration matrix cone projection */
+  scs_float ave_time_matrix_cone_proj;
+  /** Average time (milliseconds) per iteration for spectral vector cone
+   * projection */
+  scs_float ave_time_vector_cone_proj;
+#endif
 } ScsInfo;
 
 /*

include/util_spectral_cones.h

@@ -0,0 +1,45 @@
+#ifndef UTILSPECTRALCONES_H
+#define UTILSPECTRALCONES_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include "scs_blas.h"
+#include "scs_types.h"
+#include <assert.h>
+#include <math.h>
+#include <stdbool.h>
+#include <stddef.h>
+#include <stdio.h>
+
+#define IN_CONE -1
+#define IN_NEGATIVE_DUAL_CONE -2
+#define ANALYTICAL_SOL -3
+
+bool is_pos(const scs_float *x, scs_int n);
+bool is_negative(const scs_float *x, scs_int n);
+void non_neg_proj(const scs_float *src, scs_float *dst, scs_int n);
+scs_float sum_log(const scs_float *x, scs_int n);
+scs_float min_vec(const scs_float *vec, scs_int n);
+
+// used for sorting in ell1-norm cone and sum of largest cone.
+typedef struct {
+  scs_float value;
+  int index;
+} Value_index;
+
+typedef struct {
+  int iter;
+
+  // if plain Newton computed the projection or if an IPM was used
+  int newton_success;
+
+  // dual_res, pri_res, complementarity for the projection problem
+  scs_float residuals[3];
+} Newton_stats;
+
+#ifdef __cplusplus
+}
+#endif
+#endif