From 5f3256f1441174b9fdfcc72c8ec39401e3ada130 Mon Sep 17 00:00:00 2001 From: dance858 Date: Wed, 13 Nov 2024 08:40:05 -0800 Subject: [PATCH] fixed build systems, added timing, added flags --- CMakeLists.txt | 4 + Makefile | 7 +- include/cones.h | 14 +++- include/scs.h | 6 ++ include/util_spectral_cones.h | 2 +- src/cones.c | 29 +++++-- src/scs.c | 6 +- .../logdeterminant/log_cone_IPM.c | 6 +- .../logdeterminant/log_cone_Newton.c | 77 +++++++++++-------- .../logdeterminant/log_cone_wrapper.c | 17 ++-- .../logdeterminant/logdet_cone.c | 27 +++++-- src/spectral_cones/nuclear/ell1_cone.c | 12 ++- src/spectral_cones/nuclear/nuclear_cone.c | 9 +++ .../sum-largest/sum_largest_cone.c | 12 ++- .../sum-largest/sum_largest_eval_cone.c | 17 +++- src/spectral_cones/util_spectral_cones.c | 2 +- test/spectral_cones_problems/exp_design.h | 4 +- .../graph_partitioning.h | 2 +- test/spectral_cones_problems/robust_pca.h | 2 +- .../several_logdet_cones.h | 4 +- .../several_nuc_cone.h | 2 +- .../several_sum_largest.h | 2 +- 22 files changed, 187 insertions(+), 76 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 669492b4..3164331c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -163,6 +163,10 @@ if(USE_OPENMP) set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fopenmp") endif() +if(SPECTRAL_TIMING_FLAG) + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DSPECTRAL_TIMING_FLAG") +endif() + message(STATUS "COMPILER_OPTS = ${COMPILER_OPTS}") message(STATUS "CMAKE_C_FLAGS = ${CMAKE_C_FLAGS}") diff --git a/Makefile b/Makefile index f7796604..fab4ab9a 100644 --- a/Makefile +++ b/Makefile @@ -1,7 +1,12 @@ # 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 +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 SCS_O = src/scs.o SCS_INDIR_O = src/scs_indir.o diff --git a/include/cones.h b/include/cones.h index cebaeab2..ce4afbe1 100644 --- a/include/cones.h +++ b/include/cones.h @@ -5,6 +5,13 @@ 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" @@ -14,6 +21,7 @@ extern "C" { + /* private data to help cone projection step */ struct SCS_CONE_WORK { /* @@ -31,12 +39,16 @@ struct SCS_CONE_WORK { scs_float box_t_warm_start; /* if the projection onto the logarithmic cone should be warmstarted*/ - bool log_cone_warmstart; + 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; diff --git a/include/scs.h b/include/scs.h index eb0460c7..72d46060 100644 --- a/include/scs.h +++ b/include/scs.h @@ -218,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; /* diff --git a/include/util_spectral_cones.h b/include/util_spectral_cones.h index eb347f14..c0bfff52 100644 --- a/include/util_spectral_cones.h +++ b/include/util_spectral_cones.h @@ -16,7 +16,7 @@ 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 diff --git a/src/cones.c b/src/cones.c index 73940b4c..cfc1c1c3 100644 --- a/src/cones.c +++ b/src/cones.c @@ -48,7 +48,8 @@ scs_float SCS(proj_pd_exp_cone)(scs_float *v0, scs_int primal); // Forward declare spectral matrix cone projections scs_int SCS(proj_logdet_cone)(scs_float *tvX, const scs_int n, ScsConeWork *c, - Newton_stats *stats, scs_int offset); + Newton_stats *stats, scs_int offset, + bool *warmstart); scs_int SCS(proj_nuclear_cone)(scs_float *tX, size_t m, size_t n, ScsConeWork *c); void SCS(proj_ell_one)(scs_float *tx, size_t n, ScsConeWork *c); scs_int SCS(proj_sum_largest_evals)(scs_float *tX, scs_int n, scs_int k, @@ -448,6 +449,9 @@ void SCS(finish_cone)(ScsConeWork *c) { if (c->work_sum_of_largest) { scs_free(c->work_sum_of_largest); } + if (c->log_cone_warmstarts) { + scs_free(c->log_cone_warmstarts); + } #endif if (c->work_ell1) { scs_free(c->work_ell1); @@ -556,6 +560,9 @@ static scs_int set_up_sd_cone_work_space(ScsConeWork *c, const ScsCone *k) { #define _STR(tok) _STR_EXPAND(tok) scs_printf("BLAS(func) = '%s'\n", _STR(BLAS(func))); #endif + + c->log_cone_warmstarts = (bool *)scs_calloc(k->dsize, sizeof(bool)); + // ---------------------------------------------------------------------- // compute max dimension needed for eigenvalue decomposition workspace // (all cones appearing in the next section of code require eigenvalue @@ -776,10 +783,6 @@ static scs_int proj_semi_definite_cone(scs_float *X, const scs_int n, return 0; } -#ifdef DEBUG - scs_printf("number of positive eigenvalues: %d \n", n - first_idx); -#endif - /* Z is matrix of eigenvectors with positive eigenvalues */ memcpy(Z, &Xs[first_idx * n], sizeof(scs_float) * n * (n - first_idx)); @@ -1005,8 +1008,9 @@ static scs_int proj_cone(scs_float *x, const ScsCone *k, ScsConeWork *c, scs_int normalize, scs_float *r_y) { scs_int i, status; scs_int count = 0; - scs_int offset_logdet_cone = 0; /* used for warmstarting log-cone projections */ + scs_int offset_log_cone = 0; /* used for warmstarting log-cone projections */ scs_float *r_box = SCS_NULL; + SPECTRAL_TIMING(SCS(timer) spec_mat_proj_timer;) if (k->z) { /* doesn't use r_y */ /* project onto primal zero / dual free cone */ @@ -1043,7 +1047,9 @@ static scs_int proj_cone(scs_float *x, const ScsCone *k, ScsConeWork *c, if (k->ssize && k->s) { /* doesn't use r_y */ /* project onto PSD cones */ for (i = 0; i < k->ssize; ++i) { + SPECTRAL_TIMING(SCS(tic)(&spec_mat_proj_timer);) status = proj_semi_definite_cone(&(x[count]), k->s[i], c); + SPECTRAL_TIMING(c->tot_time_mat_cone_proj += SCS(tocq)(&spec_mat_proj_timer);) if (status < 0) { return status; } @@ -1093,9 +1099,12 @@ static scs_int proj_cone(scs_float *x, const ScsCone *k, ScsConeWork *c, /* project onto logdet cones */ if (k->dsize && k->d) { for (i = 0; i < k->dsize; ++i) { + SPECTRAL_TIMING(SCS(tic)(&spec_mat_proj_timer);) status = SCS(proj_logdet_cone)(&(x[count]), k->d[i], c, &(c->newton_stats), - offset_logdet_cone); - offset_logdet_cone += k->d[i] + 2; + offset_log_cone, + c->log_cone_warmstarts + i); + SPECTRAL_TIMING(c->tot_time_mat_cone_proj += SCS(tocq)(&spec_mat_proj_timer);) + offset_log_cone += k->d[i] + 2; if (status < 0) { return status; } @@ -1105,7 +1114,9 @@ static scs_int proj_cone(scs_float *x, const ScsCone *k, ScsConeWork *c, /* project onto nuclear norm cones */ if (k->nucsize && k->nuc_m && k->nuc_n) { for (i = 0; i < k->nucsize; ++i) { + SPECTRAL_TIMING(SCS(tic)(&spec_mat_proj_timer);) status = SCS(proj_nuclear_cone)(&(x[count]), k->nuc_m[i], k->nuc_n[i], c); + SPECTRAL_TIMING(c->tot_time_mat_cone_proj += SCS(tocq)(&spec_mat_proj_timer);) if (status < 0) { return status; } @@ -1122,7 +1133,9 @@ static scs_int proj_cone(scs_float *x, const ScsCone *k, ScsConeWork *c, /* project onto sum-of-largest eigenvalues cone */ if (k->sl_size && k->sl_n && k->sl_k) { for (i = 0; i < k->sl_size; ++i) { + SPECTRAL_TIMING(SCS(tic)(&spec_mat_proj_timer);) status = SCS(proj_sum_largest_evals)(&(x[count]), k->sl_n[i], k->sl_k[i], c); + SPECTRAL_TIMING(c->tot_time_mat_cone_proj += SCS(tocq)(&spec_mat_proj_timer);) if (status < 0) { return status; } diff --git a/src/scs.c b/src/scs.c index 34375737..65592c20 100644 --- a/src/scs.c +++ b/src/scs.c @@ -550,6 +550,10 @@ static void finalize(ScsWork *w, ScsSolution *sol, ScsInfo *info, info->rejected_accel_steps = w->rejected_accel_steps; info->accepted_accel_steps = w->accepted_accel_steps; info->comp_slack = ABS(sty); +#ifdef SPECTRAL_TIMING_FLAG + info->ave_time_matrix_cone_proj = w->cone_work->tot_time_mat_cone_proj / iter; + info->ave_time_vector_cone_proj = w->cone_work->tot_time_vec_cone_proj / iter; +#endif if (info->comp_slack > 1e-5 * MAX(nm_s, nm_y)) { scs_printf("WARNING - large complementary slackness residual: %f\n", info->comp_slack); @@ -1050,7 +1054,7 @@ scs_int scs_solve(ScsWork *w, ScsSolution *sol, ScsInfo *info, ScsSettings *stgs = w->stgs; /* set warm start */ stgs->warm_start = warm_start; - w->cone_work->log_cone_warmstart = false; + /* initialize ctrl-c support */ scs_start_interrupt_listener(); diff --git a/src/spectral_cones/logdeterminant/log_cone_IPM.c b/src/spectral_cones/logdeterminant/log_cone_IPM.c index 723cea8b..ce353028 100644 --- a/src/spectral_cones/logdeterminant/log_cone_IPM.c +++ b/src/spectral_cones/logdeterminant/log_cone_IPM.c @@ -444,7 +444,7 @@ scs_int log_cone_IPM(scs_float t0, scs_float v0, scs_float *x0, scs_float *u1, scs_float pres0, dres0; scs_float phi0 = 0.0, dphi0 = 0.0, step_size0 = 0.0; -#ifdef DEBUG +#ifdef SPECTRAL_DEBUG printf("%-3s%-15s%-15s%-15s%-10s%-10s\n", "", "gap", "pres", "dres", "sig", "step"); #endif @@ -498,7 +498,7 @@ scs_int log_cone_IPM(scs_float t0, scs_float v0, scs_float *x0, scs_float *u1, if (dres < FEASTOL_IPM && pres < FEASTOL_IPM && (gap < ABSTOL_IPM || relgap <= RELTOL_IPM)) { -#ifdef DEBUG +#ifdef SPECTRAL_DEBUG printf("optimal solution found: \n"); printf("gap / pres / dres: %.7e, %.7e, %.7e \n", gap, pres, dres); #endif @@ -710,7 +710,7 @@ scs_int log_cone_IPM(scs_float t0, scs_float v0, scs_float *x0, scs_float *u1, r = rnew; memcpy(z, z_new, 3 * sizeof(*z)); memcpy(s, s_new, 3 * sizeof(*s)); -#ifdef DEBUG +#ifdef SPECTRAL_DEBUG printf("%ld: %.7e, %.7e, %.7e, %f, %.3f \n", iter, gap, pres, dres, sigma, step_size); #endif diff --git a/src/spectral_cones/logdeterminant/log_cone_Newton.c b/src/spectral_cones/logdeterminant/log_cone_Newton.c index fa5ae214..a7531922 100644 --- a/src/spectral_cones/logdeterminant/log_cone_Newton.c +++ b/src/spectral_cones/logdeterminant/log_cone_Newton.c @@ -18,11 +18,12 @@ */ -#define LINESEARCH_RELATIVE_TOL 1e-15 +#define LINESEARCH_RELATIVE_TOL 1e-14 #define MIN_INIT_LOG_CONE 1 #define MIN_DENOMINATOR 1e-14 #define MIN_X 1e-17 #define MIN_FLOAT MIN_X / 2 +#define MIN_V 1e-14 #define MAX_ITER_NEWTON 100 #define ALPHA_NEWTON 0.01 @@ -30,18 +31,21 @@ #define TOL_NEWTON 1e-12 #define MAX_GRAD_STEPS 5 -#define MIN_DENOMINATOR_FAILURE -5 +#define TERMINATE_DUE_TO_ZEROS -5 #define MAX_GRAD_STEPS_REACHED -6 - -static scs_float objVal(const scs_float *u, scs_float t0, scs_float v0, - const scs_float *x0, scs_int n) +// the CALLER of this function must make sure that the arguments belong to the +// domain +static scs_float obj_val(const scs_float *u, scs_float t0, scs_float v0, + const scs_float *x0, scs_int n) { - const scs_float *v = u; + scs_float v = u[0]; const scs_float *x = u + 1; - scs_float sx = -((*v) * sum_log(x, n) - (int)(n) * (*v) * log(*v)); - scs_float obj = 0.5 * (sx - t0) * (sx - t0) + 0.5 * (*v - v0) * (*v - v0); + assert(v > 0 && min_vec(x, n) > 0); + + scs_float sx = -(v * sum_log(x, n) - n * v * log(v)); + scs_float obj = 0.5 * (sx - t0) * (sx - t0) + 0.5 * (v - v0) * (v - v0); for (scs_int i = 0; i < n; ++i) { obj += 0.5 * (x[i] - x0[i]) * (x[i] - x0[i]); @@ -51,7 +55,7 @@ static scs_float objVal(const scs_float *u, scs_float t0, scs_float v0, scs_int log_cone_Newton(scs_float t0, scs_float v0, const scs_float *x0, scs_float *u, scs_int n, scs_float *workspace, - Newton_stats *stats, bool warm_start) + Newton_stats *stats, bool *warm_start) { scs_float *t = u; scs_float *v = u + 1; @@ -91,6 +95,7 @@ scs_int log_cone_Newton(scs_float t0, scs_float v0, const scs_float *x0, *v = v0; memcpy(x, x0, sizeof(*x0) * n); stats->iter = IN_CONE; + *warm_start = true; return 0; } @@ -108,6 +113,8 @@ scs_int log_cone_Newton(scs_float t0, scs_float v0, const scs_float *x0, { memset(u, 0, (n + 2) * sizeof(*x0)); stats->iter = IN_NEGATIVE_DUAL_CONE; + // if 0 is the solution we should not use it to warmstart the next iteration + *warm_start = false; return 0; } } @@ -119,16 +126,17 @@ scs_int log_cone_Newton(scs_float t0, scs_float v0, const scs_float *x0, *v = 0; non_neg_proj(x0, x, n); stats->iter = ANALYTICAL_SOL; + *warm_start = true; return 0; } // ---------------------------------------------------------------------- - // if 'initialize' is true we initialize in the point + // if 'warm_start' is false we initialize in the point // (v, x) = (max(v0, MIN_INIT_LOG_CONE), max(x0, MIN_INIT_LOG_CONE)), // otherwise it is assumed that 'proj' has already been // initialized / warmstarted. // ---------------------------------------------------------------------- - if (!warm_start) + if (!(*warm_start)) { *v = (v0 > MIN_INIT_LOG_CONE) ? v0 : MIN_INIT_LOG_CONE; @@ -138,6 +146,8 @@ scs_int log_cone_Newton(scs_float t0, scs_float v0, const scs_float *x0, } } + scs_float obj_old = obj_val(u + 1, t0, v0, x0, n); + // ----------------------------- // parse workspace // ----------------------------- @@ -153,6 +163,14 @@ scs_int log_cone_Newton(scs_float t0, scs_float v0, const scs_float *x0, scs_float newton_decrement = 0; for (iter = 1; iter <= MAX_ITER_NEWTON; ++iter) { + // A small value on v indicates that Newton's method converges to the + // origin. In this case we should abort and apply an IPM. + if (*v < MIN_V) + { + stats->iter = TERMINATE_DUE_TO_ZEROS; + return -1; + } + // ------------------------------------------------------------------- // To avoid pathological cases where some components of // x approach 0 we use a minimum threshold. @@ -165,23 +183,24 @@ scs_int log_cone_Newton(scs_float t0, scs_float v0, const scs_float *x0, // ---------------------------------------------------------------- // compute gradient and Hessian // ---------------------------------------------------------------- - assert(*v > MIN_FLOAT && minVector(x, n) > MIN_FLOAT); + assert(*v > MIN_FLOAT && min_vec(x, n) > MIN_FLOAT); scs_float temp0 = -sum_log(x, n) + n * log(*v); scs_float a = (*v) * temp0 - t0; scs_float c = temp0 + n; grad[0] = a * c + (*v) - v0; - scs_float vInv = 1 / (*v); - d[0] = 1 + a * (-a * (vInv * vInv) + n * vInv - 2 * c * vInv); - w[0] = -(a + (*v) * c) * vInv; + scs_float v_inv = 1 / (*v); + d[0] = 1 + a * (-a * (v_inv * v_inv) + n * v_inv - 2 * c * v_inv); + w[0] = -(a + (*v) * c) * v_inv; scs_float av = a * (*v); for (i = 1; i < n + 1; ++i) { - const scs_float xInv = 1 / x[i - 1]; - grad[i] = -av * xInv + x[i - 1] - x0[i - 1]; - d[i] = 1 + av * (xInv * xInv); - w[i] = (*v) * xInv; + assert(x[i-1] > 0); + scs_float x_inv = 1 / x[i - 1]; + grad[i] = -av * x_inv + x[i - 1] - x0[i - 1]; + d[i] = 1 + av * (x_inv * x_inv); + w[i] = (*v) * x_inv; } // ---------------------------------------------------------------------- @@ -206,7 +225,7 @@ scs_int log_cone_Newton(scs_float t0, scs_float v0, const scs_float *x0, if (fabs(denominator) < MIN_DENOMINATOR) { - stats->iter = MIN_DENOMINATOR_FAILURE; + stats->iter = TERMINATE_DUE_TO_ZEROS; return -1; } @@ -254,38 +273,36 @@ scs_int log_cone_Newton(scs_float t0, scs_float v0, const scs_float *x0, { if (du[i] < 0) { - scs_float maxStep = -0.99 * u[i + 1] / du[i]; - if (step_size > maxStep) - { - step_size = maxStep; - } + scs_float max_step = -0.99 * u[i + 1] / du[i]; + step_size = MIN(step_size, max_step); } } // ------------------------------------------------- // backtracking line search. First two lines do // u_new = u + t * du; - // TODO: we don't have to recompute obj_old here! // ------------------------------------------------- memcpy(u_new + 1, u + 1, n_plus_one * sizeof(*u)); SCS(add_scaled_array)(u_new + 1, du, n_plus_one, step_size); - scs_float obj_old = objVal(u + 1, t0, v0, x0, n); - scs_float new_obj = objVal(u_new + 1, t0, v0, x0, n); + scs_float new_obj = obj_val(u_new + 1, t0, v0, x0, n); while ((1 - LINESEARCH_RELATIVE_TOL) * new_obj > obj_old + ALPHA_NEWTON * step_size * dirDer) { step_size *= BETA_NEWTON; memcpy(u_new + 1, u + 1, n_plus_one * sizeof(*u)); SCS(add_scaled_array)(u_new + 1, du, n_plus_one, step_size); - new_obj = objVal(u_new + 1, t0, v0, x0, n); + new_obj = obj_val(u_new + 1, t0, v0, x0, n); } + obj_old = new_obj; memcpy(u + 1, u_new + 1, n_plus_one * sizeof(*u)); } - assert(minVector(x, n) > MIN_FLOAT && *v > MIN_FLOAT); + assert(min_vec(x, n) > MIN_FLOAT && *v > MIN_FLOAT); *t = -(*v) * (sum_log(x, n) - n*log(*v)); + + *warm_start = true; stats->iter = iter; return 0; } \ No newline at end of file diff --git a/src/spectral_cones/logdeterminant/log_cone_wrapper.c b/src/spectral_cones/logdeterminant/log_cone_wrapper.c index 073360c6..3b7a5018 100644 --- a/src/spectral_cones/logdeterminant/log_cone_wrapper.c +++ b/src/spectral_cones/logdeterminant/log_cone_wrapper.c @@ -31,7 +31,7 @@ // forward declare from log_cone_Newton.c scs_int log_cone_Newton(scs_float t0, scs_float v0, const scs_float *x0, scs_float *proj, scs_int n, scs_float *workspace, - Newton_stats *stats, bool warm_start); + Newton_stats *stats, bool *warm_start); // forward declare form log_cone_IPM.c scs_int log_cone_IPM(scs_float t0, scs_float v0, scs_float *x0, scs_float *u1, @@ -45,7 +45,7 @@ static void check_opt_cond_log_cone(const scs_float *tvx, scs_float t0, scs_int log_cone_proj_wrapper(scs_float t0, scs_float v0, scs_float *x0, scs_float *proj, scs_int n, scs_float *workspace, - Newton_stats *stats, bool warm_start) + Newton_stats *stats, bool *warm_start) { scs_int status; // ----------------------------------------------------------------------- @@ -85,6 +85,8 @@ scs_int log_cone_proj_wrapper(scs_float t0, scs_float v0, scs_float *x0, status = log_cone_IPM(t0, v0, x0, proj, n, workspace, stats, 0); check_opt_cond_log_cone(proj, t0, v0, x0, n, stats->residuals, workspace); + *warm_start = true; // next iteration Newton should be warmstarted + if (status == 0 && stats->residuals[0] < DUAL_FEAS_TOL && stats->residuals[1] < PRI_FEAS_TOL && fabs(stats->residuals[2]) < COMP_TOL) { @@ -94,7 +96,12 @@ scs_int log_cone_proj_wrapper(scs_float t0, scs_float v0, scs_float *x0, // ------------------------------------------------------------------------ // Solve problem with primal-dual IPM without Mehrotra's correction - // etc. + // etc. (In all experiments in the paper by Cederberg and Boyd, + // the IPM above solves the problem correctly so the code below never + // runs. However, during development I ran into a (possibly pathological + // case) where the first IPM fails. If this happens we run below another + // version. This version of the IPM solves the problem that the first version + // failed on.) // ------------------------------------------------------------------------ status = log_cone_IPM(t0, v0, x0, proj, n, workspace, stats, 1); check_opt_cond_log_cone(proj, t0, v0, x0, n, stats->residuals, workspace); @@ -106,8 +113,8 @@ scs_int log_cone_proj_wrapper(scs_float t0, scs_float v0, scs_float *x0, return 0; } -#ifdef DEBUG - scs_printf("FAILURE: logarithmic cone projection") +#ifdef SPECTRAL_DEBUG + scs_printf("FAILURE: logarithmic cone projection"); scs_printf("dual_res / pri_res / comp / iter: %.3e, %.3e, %.3e, %d\n", stats->residuals[0], stats->residuals[1], stats->residuals[2], stats->iter); diff --git a/src/spectral_cones/logdeterminant/logdet_cone.c b/src/spectral_cones/logdeterminant/logdet_cone.c index 732347bd..0da2077d 100644 --- a/src/spectral_cones/logdeterminant/logdet_cone.c +++ b/src/spectral_cones/logdeterminant/logdet_cone.c @@ -18,13 +18,24 @@ * Last modified: 25 August 2024. */ +void BLAS(syev)(const char *jobz, const char *uplo, blas_int *n, scs_float *a, + blas_int *lda, scs_float *w, scs_float *work, blas_int *lwork, + blas_int *info); +blas_int BLAS(syrk)(const char *uplo, const char *trans, const blas_int *n, + const blas_int *k, const scs_float *alpha, + const scs_float *a, const blas_int *lda, + const scs_float *beta, scs_float *c, const blas_int *ldc); +void BLAS(scal)(const blas_int *n, const scs_float *sa, scs_float *sx, + const blas_int *incx); + // forward declare from log_cone_wrapper.c scs_int log_cone_proj_wrapper(scs_float t0, scs_float v0, const scs_float *x0, scs_float *proj, scs_int n, scs_float *workspace, - Newton_stats *stats, bool warm_start); + Newton_stats *stats, bool *warm_start); scs_int SCS(proj_logdet_cone)(scs_float *tvX, scs_int n, ScsConeWork *c, - Newton_stats *stats, scs_int offset) + Newton_stats *stats, scs_int offset, + bool *warmstart) { // tvX = [t, v, X], where X represents the lower triangular part of a matrix // stored in a compact form and off-diagonal elements have been scaled by @@ -91,15 +102,17 @@ scs_int SCS(proj_logdet_cone)(scs_float *tvX, scs_int n, ScsConeWork *c, // (or equivalently, in c->saved_log_projs + offset) // ---------------------------------------------------------------------- scs_float *current_log_proj = c->saved_log_projs + offset; + SPECTRAL_TIMING(SCS(timer) _timer; SCS(tic)(&_timer);) scs_int status = log_cone_proj_wrapper(sqrt2 * tvX[0], sqrt2 * tvX[1], e, current_log_proj, n, c->work_logdet, - stats, c->log_cone_warmstart); - + stats, warmstart); + SPECTRAL_TIMING(c->tot_time_vec_cone_proj += SCS(tocq)(&_timer);) + if (status < 0) { return status; } - // return immediateely if the origin is the solution + // return immediately if the origin is the solution else if (status == IN_NEGATIVE_DUAL_CONE) { memset(tvX, 0, sizeof(scs_float) * ((n * (n + 1)) / 2 + 2)); @@ -112,7 +125,7 @@ scs_int SCS(proj_logdet_cone)(scs_float *tvX, scs_int n, ScsConeWork *c, scs_float *evals_proj = current_log_proj + 2; for (i = 0; i < n; ++i) { - assert(evals_proj[i] > 0); + assert(evals_proj[i] >= 0); sq_eig_pos = SQRTF(evals_proj[i]); BLAS(scal) (&nb, &sq_eig_pos, &Xs[i * n], &one_int); @@ -131,7 +144,5 @@ scs_int SCS(proj_logdet_cone)(scs_float *tvX, scs_int n, ScsConeWork *c, tvX[0] = sqrt2_inv * current_log_proj[0]; tvX[1] = sqrt2_inv * current_log_proj[1]; - // after the first iteration we can warm start the spectral projection - c->log_cone_warmstart = true; return 0; } diff --git a/src/spectral_cones/nuclear/ell1_cone.c b/src/spectral_cones/nuclear/ell1_cone.c index 635e1e5b..6e8712de 100644 --- a/src/spectral_cones/nuclear/ell1_cone.c +++ b/src/spectral_cones/nuclear/ell1_cone.c @@ -19,6 +19,13 @@ * Last modified: 25 August 2024. */ +void BLAS(axpy)(blas_int *n, const scs_float *a, const scs_float *x, + blas_int *incx, scs_float *y, blas_int *incy); + +scs_float BLAS(dot)(const blas_int *n, const scs_float *x, const blas_int *incx, + const scs_float *y, const blas_int *incy); + +#ifdef SPECTRAL_DEBUG static void compute_cone_residuals_ell1(const scs_float *tx, scs_float t0, const scs_float *x0, scs_int n, scs_float residuals[3]) @@ -77,14 +84,13 @@ static void compute_cone_residuals_ell1(const scs_float *tx, scs_float t0, residuals[2] = complementarity; free(dualx); } +#endif // Asssumes that all components of x0 are positive and // x0[0] >= x0[1] >= ... x0[n-1]. scs_int ell1_cone_proj_sorted(scs_float t0, const scs_float *x0, scs_float *proj, scs_int n) { - assert(isSortedAndNonNeg(x0, n)); - if (-t0 >= x0[0]) { memset(proj, 0, (n + 1) * sizeof(*x0)); @@ -145,7 +151,7 @@ scs_int ell1_cone_proj_sorted(scs_float t0, const scs_float *x0, scs_float *proj } memset(proj + 1 + k, 0, (n - k) * sizeof(*x0)); -#ifdef DEBUG +#ifdef SPECTRAL_DEBUG //------------------------------------------------------------------------- // Check residuals - not needed in production //------------------------------------------------------------------------- diff --git a/src/spectral_cones/nuclear/nuclear_cone.c b/src/spectral_cones/nuclear/nuclear_cone.c index 9e6c899e..7638500c 100644 --- a/src/spectral_cones/nuclear/nuclear_cone.c +++ b/src/spectral_cones/nuclear/nuclear_cone.c @@ -16,6 +16,13 @@ * * Last modified: 25 August 2024. */ +void BLAS(gemm)(const char *transa, const char *transb, blas_int *m, + blas_int *n, blas_int *k, scs_float *alpha, scs_float *a, + blas_int *lda, scs_float *b, blas_int *ldb, scs_float *beta, + scs_float *c, blas_int *ldc); + +void BLAS(scal)(const blas_int *n, const scs_float *sa, scs_float *sx, + const blas_int *incx); void BLAS(gesvd)(const char *jobu, const char *jobvt, const blas_int *m, const blas_int *n, scs_float *a, const blas_int *lda, @@ -58,7 +65,9 @@ scs_int SCS(proj_nuclear_cone)(scs_float *tX, size_t m, size_t n, ScsConeWork *c // ------------------------------------------------------------------------- // Project onto spectral *vector* cone // ------------------------------------------------------------------------- + SPECTRAL_TIMING(SCS(timer) _timer; SCS(tic)(&_timer);) scs_int status = ell1_cone_proj_sorted(tX[0], s, tX, n); + SPECTRAL_TIMING(c->tot_time_vec_cone_proj += SCS(tocq)(&_timer);) if (status < 0) { diff --git a/src/spectral_cones/sum-largest/sum_largest_cone.c b/src/spectral_cones/sum-largest/sum_largest_cone.c index e3e4a988..f63da1b6 100644 --- a/src/spectral_cones/sum-largest/sum_largest_cone.c +++ b/src/spectral_cones/sum-largest/sum_largest_cone.c @@ -22,9 +22,11 @@ * Last modified: 25 August 2024. */ +#ifdef SPECTRAL_DEBUG static void compute_cone_residuals(scs_float t, const scs_float *x, scs_float t0, const scs_float *x0, scs_float residuals[3], scs_int n, scs_int k); +#endif scs_int assert_sorted(scs_float *x, scs_int n) { @@ -43,7 +45,7 @@ scs_int proj_sum_largest_cone_sorted(scs_float *t, scs_float *x, scs_int n, scs_int k) { -#ifdef DEBUG +#ifdef SPECTRAL_DEBUG scs_float t00 = *t; scs_float *x0 = scs_malloc(n * sizeof(*x0)); memcpy(x0, x, n * sizeof(*x0)); @@ -134,7 +136,7 @@ scs_int proj_sum_largest_cone_sorted(scs_float *t, scs_float *x, scs_int n, x[i] = a_t; } -#ifdef DEBUG +#ifdef SPECTRAL_DEBUG scs_float residuals[3]; compute_cone_residuals(*t, x, t00, x0, residuals, n, k); scs_free(x0); @@ -158,6 +160,8 @@ scs_int cmp_desc(const void *a, const void *b) return 0; } + +#ifdef SPECTRAL_DEBUG // this function is not used in production so fine to allocate memory static scs_float sum_largest_val(const scs_float *x, scs_int n, scs_int k) { @@ -181,7 +185,6 @@ static void compute_cone_residuals(scs_float t, const scs_float *x, scs_float t0 scs_float lmbda_t = t - t0; scs_float *lmbda_x = scs_malloc(n * sizeof(*x)); scs_float sum_lmbda_x = 0.0; - scs_int int_one = 1; for (scs_int i = 0; i < n; ++i) { lmbda_x[i] = x[i] - x0[i]; @@ -212,4 +215,5 @@ static void compute_cone_residuals(scs_float t, const scs_float *x, scs_float t0 residuals[2] = comp; scs_free(lmbda_x); -} \ No newline at end of file +} +#endif \ No newline at end of file diff --git a/src/spectral_cones/sum-largest/sum_largest_eval_cone.c b/src/spectral_cones/sum-largest/sum_largest_eval_cone.c index 25128afa..7c6f97de 100644 --- a/src/spectral_cones/sum-largest/sum_largest_eval_cone.c +++ b/src/spectral_cones/sum-largest/sum_largest_eval_cone.c @@ -17,6 +17,18 @@ * Last modified: 25 August 2024. */ +void BLAS(syev)(const char *jobz, const char *uplo, blas_int *n, scs_float *a, + blas_int *lda, scs_float *w, scs_float *work, blas_int *lwork, + blas_int *info); + +void BLAS(scal)(const blas_int *n, const scs_float *sa, scs_float *sx, + const blas_int *incx); + +void BLAS(gemm)(const char *transa, const char *transb, blas_int *m, + blas_int *n, blas_int *k, scs_float *alpha, scs_float *a, + blas_int *lda, scs_float *b, blas_int *ldb, scs_float *beta, + scs_float *c, blas_int *ldc); + // forward declaration scs_int proj_sum_largest_cone_sorted(scs_float *t, scs_float *x, scs_int n, scs_int k); @@ -32,9 +44,6 @@ void flip(scs_float *x, int n) { scs_int SCS(proj_sum_largest_evals)(scs_float *tX, scs_int n, scs_int k, ScsConeWork *c) { - SCS(timer) - _timer; - // tvX = [t, X], where X represents the lower triangular part of a matrix // stored in a compact form and off-diagonal elements have been scaled by // sqrt(2) @@ -85,9 +94,11 @@ scs_int SCS(proj_sum_largest_evals)(scs_float *tX, scs_int n, scs_int k, // tvX[0] by sqrt(2). // ---------------------------------------------------------------------- tX[0] *= sqrt2; + SPECTRAL_TIMING(SCS(timer) _timer; SCS(tic)(&_timer);) flip(e, n); scs_int status = proj_sum_largest_cone_sorted(&tX[0], e, n, k); flip(e, n); + SPECTRAL_TIMING(c->tot_time_vec_cone_proj += SCS(tocq)(&_timer);) if (status < 0) { diff --git a/src/spectral_cones/util_spectral_cones.c b/src/spectral_cones/util_spectral_cones.c index 0b59737f..1a91ee4d 100644 --- a/src/spectral_cones/util_spectral_cones.c +++ b/src/spectral_cones/util_spectral_cones.c @@ -41,7 +41,7 @@ void print_vector(const scs_float *x, size_t n) printf("\n"); } -scs_float min_vec(const scs_float *vec, size_t n) +scs_float min_vec(const scs_float *vec, scs_int n) { scs_float minVal = vec[0]; diff --git a/test/spectral_cones_problems/exp_design.h b/test/spectral_cones_problems/exp_design.h index 44929ec3..73bc9eb7 100644 --- a/test/spectral_cones_problems/exp_design.h +++ b/test/spectral_cones_problems/exp_design.h @@ -7,7 +7,6 @@ #include "scs_matrix.h" #include "util.h" -double sqrt2 = sqrt(2); // for SpectralSCS static const char *exp_design(void) @@ -19,7 +18,7 @@ static const char *exp_design(void) ScsInfo info = {0}; scs_int exitflag; scs_float perr, derr; - scs_int success, read_status; + scs_int success; const char *fail; /* data */ @@ -109,6 +108,7 @@ static const char *exp_design(void) stgs->eps_rel = 1e-7; stgs->eps_infeas = 1e-9; + stgs->log_csv_filename="test_exp_design.csv"; exitflag = scs(d, k, stgs, sol, &info); perr = SCS(dot)(d->c, sol->x, d->n) - opt; diff --git a/test/spectral_cones_problems/graph_partitioning.h b/test/spectral_cones_problems/graph_partitioning.h index e2c7b235..a6f9fc6f 100644 --- a/test/spectral_cones_problems/graph_partitioning.h +++ b/test/spectral_cones_problems/graph_partitioning.h @@ -17,7 +17,7 @@ static const char *graph_partitioning(void) ScsInfo info = {0}; scs_int exitflag; scs_float perr, derr; - scs_int success, read_status; + scs_int success; const char *fail; /* data */ diff --git a/test/spectral_cones_problems/robust_pca.h b/test/spectral_cones_problems/robust_pca.h index 2d3ae50e..e526d404 100644 --- a/test/spectral_cones_problems/robust_pca.h +++ b/test/spectral_cones_problems/robust_pca.h @@ -17,7 +17,7 @@ static const char *robust_pca(void) ScsInfo info = {0}; scs_int exitflag; scs_float perr, derr; - scs_int success, read_status; + scs_int success; const char *fail; /* data */ diff --git a/test/spectral_cones_problems/several_logdet_cones.h b/test/spectral_cones_problems/several_logdet_cones.h index 15d69df2..c294fbaf 100644 --- a/test/spectral_cones_problems/several_logdet_cones.h +++ b/test/spectral_cones_problems/several_logdet_cones.h @@ -17,7 +17,7 @@ static const char *several_logdet_cones(void) ScsInfo info = {0}; scs_int exitflag; scs_float perr, derr; - scs_int success, read_status; + scs_int success; const char *fail; /* data */ @@ -184,6 +184,8 @@ static const char *several_logdet_cones(void) stgs->eps_rel = 1e-6; stgs->eps_infeas = 1e-9; + + stgs->log_csv_filename="several_logdet_cone.csv"; exitflag = scs(d, k, stgs, sol, &info); perr = SCS(dot)(d->c, sol->x, d->n) - opt; diff --git a/test/spectral_cones_problems/several_nuc_cone.h b/test/spectral_cones_problems/several_nuc_cone.h index 565e5dc8..9fc56032 100644 --- a/test/spectral_cones_problems/several_nuc_cone.h +++ b/test/spectral_cones_problems/several_nuc_cone.h @@ -17,7 +17,7 @@ static const char *several_nuc_cone(void) ScsInfo info = {0}; scs_int exitflag; scs_float perr, derr; - scs_int success, read_status; + scs_int success; const char *fail; /* data */ diff --git a/test/spectral_cones_problems/several_sum_largest.h b/test/spectral_cones_problems/several_sum_largest.h index 26e9de59..d9adbdc4 100644 --- a/test/spectral_cones_problems/several_sum_largest.h +++ b/test/spectral_cones_problems/several_sum_largest.h @@ -17,7 +17,7 @@ static const char *several_sum_largest(void) ScsInfo info = {0}; scs_int exitflag; scs_float perr, derr; - scs_int success, read_status; + scs_int success; const char *fail; /* data */