Spectral cones

.gitignore

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

CMakeLists.txt

@@ -163,6 +163,14 @@ if(USE_OPENMP)
   set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fopenmp")
 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}")
 
@@ -176,18 +184,36 @@ set(${PROJECT_NAME}_PUBLIC_HDR include/scs_types.h include/scs.h)
 
 # Common source files
 set(${PROJECT_NAME}_SRC
-    src/aa.c
-    src/cones.c
-    src/ctrlc.c
-    src/exp_cone.c
-    src/linalg.c
-    src/normalize.c
-    src/rw.c
-    src/scs.c
-    src/scs_version.c
-    src/util.c
-    ${LINSYS}/csparse.c
-    ${LINSYS}/scs_matrix.c)
+  src/aa.c
+  src/cones.c
+  src/ctrlc.c
+  src/exp_cone.c
+  src/linalg.c
+  src/normalize.c
+  src/rw.c
+  src/scs.c
+  src/scs_version.c
+  src/util.c
+  ${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_LAPACK)
+  list(APPEND ${PROJECT_NAME}_SRC ${SPECTRAL_SRC})
+endif()
+
 
 # Common header files
 set(${PROJECT_NAME}_HDR
@@ -207,6 +233,15 @@ set(${PROJECT_NAME}_HDR
     ${LINSYS}/csparse.h
     ${LINSYS}/scs_matrix.h)
 
+# Spectral header files
+set(SPECTRAL_HDR
+  include/util_spectral_cones.h)
+
+if(USE_LAPACK)
+  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)
 

Makefile

@@ -1,12 +1,27 @@
 # 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
+ifneq ($(USE_LAPACK), 0)
+    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 \
+             	  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
+else
+	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
+endif
+
 SCS_O = src/scs.o
 SCS_INDIR_O = src/scs_indir.o
 
 SRC_FILES = $(wildcard src/*.c)
-INC_FILES = $(wildcard include/*.h)
+
+ifneq ($(USE_LAPACK), 0)
+	INC_FILES = $(wildcard include/*.h)
+else
+	INC_FILES = $(filter-out include/util_spectral_cones.h, $(wildcard include/*.h))
+endif
 
 AMD_SOURCE = $(wildcard $(EXTSRC)/amd/*.c)
 LDL_SOURCE = $(EXTSRC)/qdldl/qdldl.c

include/cones.h

@@ -5,11 +5,22 @@
 extern "C" {
 #endif
 
+// macro for time measurements of SpectralSCS
+#ifdef SPECTRAL_TIMING_FLAG
+    #define SPECTRAL_TIMING(action) action
+#else
+    #define SPECTRAL_TIMING(action)
+#endif
+
 #include "glbopts.h"
 #include "scs.h"
 #include "scs_blas.h"
 #include "scs_work.h"
 #include <string.h>
+#include "util_spectral_cones.h" // for newton_stats
+
+
+
 
 /* private data to help cone projection step */
 struct SCS_CONE_WORK {
@@ -26,10 +37,40 @@ struct SCS_CONE_WORK {
   scs_int m;            /* total length of cone */
   /* box cone quantities */
   scs_float box_t_warm_start;
+
+  /* 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;)
+
 #ifdef USE_LAPACK
   /* workspace for eigenvector decompositions: */
   scs_float *Xs, *Z, *e, *work;
   blas_int lwork;
+
+  /* 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
 };
 

include/scs.h

@@ -139,6 +139,22 @@ typedef struct {
   scs_float *p;
   /** Number of (primal and dual) power cone triples. */
   scs_int psize;
+  /** 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;
 } ScsCone;
 
 /** Contains primal-dual solution arrays or a certificate of infeasibility.
@@ -202,6 +218,12 @@ 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/scs_types.h

@@ -1,3 +1,7 @@
+/*
+ * GENERATED BY CMAKE
+ */
+
 /*
  * Pulic header including definitions of primitive types used in SCS.
  * Make sure this file and `scs.h` are somewhere appropriate and then use
@@ -11,24 +15,8 @@
 extern "C" {
 #endif
 
-#ifdef DLONG
-/*#ifdef _WIN64
-#include <stdint.h>
-typedef int64_t scs_int;
-#else
-typedef long scs_int;
-#endif
-*/
-typedef long long scs_int;
-#else
-typedef int scs_int;
-#endif
-
-#ifndef SFLOAT
-typedef double scs_float;
-#else
-typedef float scs_float;
-#endif
+typedef int    scs_int;
+typedef double  scs_float;
 
 #ifdef __cplusplus
 }

include/util_spectral_cones.h

@@ -0,0 +1,47 @@
+#ifndef UTILSPECTRALCONES_H
+#define UTILSPECTRALCONES_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include <stddef.h>
+#include <stdbool.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include "scs_types.h"
+#include "scs_blas.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