More intelligent initial point selection

Author: PTNobel

src/scs.c

@@ -174,21 +174,46 @@ static inline scs_int _is_nan(scs_float x) {
 }
 
 /* given x,y,s warm start, set v = [x; s / R + y; 1]
- * check for nans and set to zero if present
+ * check for nans:
+ * 	x NaN -> w->g
+ * 	y NaN -> -w->g
+ * 	s NaN -> b - Ax (primal feasibility)
+ * g is set by update_work_cache
  */
 static void warm_start_vars(ScsWork *w, ScsSolution *sol) {
   scs_int n = w->d->n, m = w->d->m, i;
   scs_float *v = w->v;
+  scs_float *g = w->g;
+  scs_int s_has_nan = 0;
   /* normalize the warm-start */
   if (w->stgs->normalize) {
     SCS(normalize_sol)(w->scal, sol);
   }
   for (i = 0; i < n; ++i) {
-    v[i] = _is_nan(sol->x[i]) ? 0. : sol->x[i];
+    v[i] = _is_nan(sol->x[i]) ? g[i] : sol->x[i];
   }
+  /* check if s has any NaN entries */
   for (i = 0; i < m; ++i) {
-    v[i + n] = sol->y[i] + sol->s[i] / w->diag_r[i + n];
-    v[i + n] = _is_nan(v[i + n]) ? 0. : v[i + n];
+    if (_is_nan(sol->s[i])) {
+      s_has_nan = 1;
+      break;
+    }
+  }
+  if (s_has_nan) {
+    /* compute Ax into v[n:n+m] as scratch space */
+    memset(v + n, 0, m * sizeof(scs_float));
+    SCS(accum_by_a)(w->d->A, v, v + n); /* v[n:n+m] = A * x */
+    /* for NaN entries of s, use b - Ax (primal feasibility) */
+    for (i = 0; i < m; ++i) {
+      scs_float si = _is_nan(sol->s[i]) ? (w->d->b[i] - v[i + n]) : sol->s[i];
+      scs_float yi = _is_nan(sol->y[i]) ? -g[i + n] : sol->y[i];
+      v[i + n] = yi + si / w->diag_r[i + n];
+    }
+  } else {
+    for (i = 0; i < m; ++i) {
+      scs_float yi = _is_nan(sol->y[i]) ? -g[i + n] : sol->y[i];
+      v[i + n] = yi + sol->s[i] / w->diag_r[i + n];
+    }
   }
   v[n + m] = 1.0; /* tau = 1 */
   /* un-normalize so sol unchanged */
@@ -339,9 +364,23 @@ static void populate_residual_struct(ScsWork *w, scs_int iter) {
 }
 
 static void cold_start_vars(ScsWork *w) {
-  scs_int l = w->d->n + w->d->m + 1;
-  memset(w->v, 0, l * sizeof(scs_float));
-  w->v[l - 1] = 1.;
+  scs_int n = w->d->n, m = w->d->m, i;
+  scs_float *v = w->v;
+  scs_float *g = w->g;
+
+  for (i = 0; i < n; ++i) {
+    v[i] = g[i];
+  }
+
+  memset(v + n, 0, m * sizeof(scs_float));
+  SCS(accum_by_a)(w->d->A, v, v + n); /* v[n:n+m] = A * x */
+  /* for NaN entries of s, use b - Ax (primal feasibility) */
+  for (i = 0; i < m; ++i) {
+    scs_float si = (w->d->b[i] - v[i + n]);
+    scs_float yi = -g[i + n];
+    v[i + n] = yi + si / w->diag_r[i + n];
+  }
+  v[n + m] = 1.;
 }
 
 /* utility function that computes x'Ry */
@@ -944,16 +983,17 @@ static void reset_tracking(ScsWork *w) {
 static scs_int update_work(ScsWork *w, ScsSolution *sol) {
   reset_tracking(w);
 
+  /* h = [c;b] */
+  memcpy(w->h, w->d->c, w->d->n * sizeof(scs_float));
+  memcpy(&(w->h[w->d->n]), w->d->b, w->d->m * sizeof(scs_float));
+  update_work_cache(w);
+
   if (w->stgs->warm_start) {
     warm_start_vars(w, sol);
   } else {
     cold_start_vars(w);
   }
 
-  /* h = [c;b] */
-  memcpy(w->h, w->d->c, w->d->n * sizeof(scs_float));
-  memcpy(&(w->h[w->d->n]), w->d->b, w->d->m * sizeof(scs_float));
-  update_work_cache(w);
   return 0;
 }
 

test/problems/partial_warm_start.h

@@ -0,0 +1,150 @@
+#include "glbopts.h"
+#include "linalg.h"
+#include "minunit.h"
+#include "problem_utils.h"
+#include "scs.h"
+#include "scs_matrix.h"
+#include "util.h"
+
+static const char *partial_warm_start(void) {
+  ScsCone *k = (ScsCone *)scs_calloc(1, sizeof(ScsCone));
+  ScsData *d = (ScsData *)scs_calloc(1, sizeof(ScsData));
+  ScsSettings *stgs = (ScsSettings *)scs_calloc(1, sizeof(ScsSettings));
+  ScsSolution *sol = (ScsSolution *)scs_calloc(1, sizeof(ScsSolution));
+  ScsInfo info = {0};
+  scs_int exitflag;
+  scs_int cold_iters;
+  scs_int i;
+  scs_float perr, derr;
+  scs_int success;
+  const char *fail;
+
+  /* data (same as hs21_tiny_qp) */
+  scs_float Ax[] = {-10., -1., 1., -1.};
+  scs_int Ai[] = {1, 2, 1, 3};
+  scs_int Ap[] = {0, 2, 4};
+
+  scs_float Px[] = {0.02, 2.};
+  scs_int Pi[] = {0, 1};
+  scs_int Pp[] = {0, 1, 2};
+
+  scs_float b[] = {1., 0., 0., 0.};
+  scs_float c[] = {0., 0.};
+
+  scs_int m = 4;
+  scs_int n = 2;
+
+  scs_float bl[] = {10.0, 2.0, -50.0};
+  scs_float bu[] = {1e+20, 50.0, 50.0};
+  scs_int bsize = 4;
+
+  scs_float opt = 0.04000000000000625;
+  /* end data */
+
+  d->m = m;
+  d->n = n;
+  d->b = b;
+  d->c = c;
+
+  d->A = (ScsMatrix *)scs_calloc(1, sizeof(ScsMatrix));
+  d->P = (ScsMatrix *)scs_calloc(1, sizeof(ScsMatrix));
+
+  d->A->m = m;
+  d->A->n = n;
+
+  d->A->x = Ax;
+  d->A->i = Ai;
+  d->A->p = Ap;
+
+  d->P->m = n;
+  d->P->n = n;
+
+  d->P->x = Px;
+  d->P->i = Pi;
+  d->P->p = Pp;
+
+  k->bsize = bsize;
+  k->bl = bl;
+  k->bu = bu;
+
+  scs_set_default_settings(stgs);
+  stgs->eps_abs = 1e-9;
+  stgs->eps_rel = 1e-9;
+  stgs->eps_infeas = 0.;
+  stgs->acceleration_lookback = 0; /* disable acceleration for consistent iters */
+
+  /* Step 1: Cold solve to get reference solution and iteration count */
+  exitflag = scs(d, k, stgs, sol, &info);
+
+  perr = info.pobj - opt;
+  derr = info.dobj - opt;
+
+  success = ABS(perr) < 1e-3 && ABS(derr) < 1e-3 && exitflag == SCS_SOLVED;
+  mu_assert("partial_warm_start: cold solve failed", success);
+  fail = verify_solution_correct(d, k, stgs, &info, sol, exitflag);
+  if (fail) {
+    SCS(free_sol)(sol);
+    scs_free(d->A);
+    scs_free(d->P);
+    scs_free(k);
+    scs_free(stgs);
+    scs_free(d);
+    return fail;
+  }
+
+  cold_iters = info.iter;
+  scs_printf("partial_warm_start: cold solve took %li iters\n",
+             (long)cold_iters);
+
+  /* Step 2: Partial warm start - keep x, set s and y to NaN */
+  for (i = 0; i < m; ++i) {
+    sol->s[i] = NAN;
+    sol->y[i] = NAN;
+  }
+
+  stgs->warm_start = 1;
+  exitflag = scs(d, k, stgs, sol, &info);
+
+  perr = info.pobj - opt;
+  derr = info.dobj - opt;
+
+  success = ABS(perr) < 1e-3 && ABS(derr) < 1e-3 && exitflag == SCS_SOLVED;
+  mu_assert("partial_warm_start: partial warm start (x only) failed to solve",
+            success);
+
+  scs_printf("partial_warm_start: partial warm start (x only) took %li iters\n",
+             (long)info.iter);
+  mu_assert(
+      "partial_warm_start: partial warm start should take fewer iters than "
+      "cold start",
+      info.iter < cold_iters);
+
+  /* Step 3: All-NaN warm start - should still converge */
+  for (i = 0; i < n; ++i) {
+    sol->x[i] = NAN;
+  }
+  for (i = 0; i < m; ++i) {
+    sol->s[i] = NAN;
+    sol->y[i] = NAN;
+  }
+
+  stgs->warm_start = 1;
+  exitflag = scs(d, k, stgs, sol, &info);
+
+  perr = info.pobj - opt;
+  derr = info.dobj - opt;
+
+  success = ABS(perr) < 1e-3 && ABS(derr) < 1e-3 && exitflag == SCS_SOLVED;
+  mu_assert("partial_warm_start: all-NaN warm start failed to solve", success);
+
+  scs_printf("partial_warm_start: all-NaN warm start took %li iters\n",
+             (long)info.iter);
+
+  SCS(free_sol)(sol);
+  scs_free(d->A);
+  scs_free(d->P);
+  scs_free(k);
+  scs_free(stgs);
+  scs_free(d);
+  return fail;
+}

test/run_tests.c

@@ -13,6 +13,7 @@
 #include "problems/small_lp.h"
 #include "problems/small_qp.h"
 #include "problems/test_exp_cone.h"
+#include "problems/partial_warm_start.h"
 #include "problems/unbounded_tiny_qp.h"
 
 int tests_run = 0;
@@ -86,6 +87,7 @@ static const char *all_tests(void) {
   mu_run_test(complex_PSD);
   mu_run_test(sd_and_complex_sd);
   mu_run_test(hs21_tiny_qp);
+  mu_run_test(partial_warm_start);
   mu_run_test(hs21_tiny_qp_rw);
   mu_run_test(qafiro_tiny_qp);
   mu_run_test(infeasible_tiny_qp);