1 #include "cones.h"
2 #include "linalg.h"
3 #include "scs.h"
4 #include "scs_blas.h" /* contains BLAS(X) macros and type info */
5 #include "util.h"
6
7 #define CONE_TOL (1e-9)
8 #define CONE_THRESH (1e-8)
9 #define EXP_CONE_MAX_ITERS (100)
10 #define BOX_CONE_MAX_ITERS (25)
11 #define POW_CONE_MAX_ITERS (20)
12
13 /* In the box cone projection we penalize the `t` term additionally by this
14 * factor. This encourages the `t` term to stay close to the incoming `t` term,
15 * which should provide better convergence since typically the `t` term does
16 * not appear in the linear system other than `t = 1`. Setting to 1 is
17 * the vanilla projection.
18 */
19 #define BOX_T_SCALE (1.)
20
21 /* Box cone limits (+ or -) taken to be INF */
22 #define MAX_BOX_VAL (1e15)
23
24 #ifdef USE_LAPACK
25
26 #ifdef __cplusplus
27 extern "C" {
28 #endif
29
30 void BLAS(syev)(const char *jobz, const char *uplo, blas_int *n, scs_float *a,
31 blas_int *lda, scs_float *w, scs_float *work, blas_int *lwork,
32 blas_int *info);
33 blas_int BLAS(syrk)(const char *uplo, const char *trans, const blas_int *n,
34 const blas_int *k, const scs_float *alpha,
35 const scs_float *a, const blas_int *lda,
36 const scs_float *beta, scs_float *c, const blas_int *ldc);
37 void BLAS(scal)(const blas_int *n, const scs_float *sa, scs_float *sx,
38 const blas_int *incx);
39
40 #ifdef __cplusplus
41 }
42 #endif
43
44 #endif
45
46 /* set the vector of rho y terms, based on scale and cones */
SCS(set_rho_y_vec)47 void SCS(set_rho_y_vec)(const ScsCone *k, scs_float scale, scs_float *rho_y_vec,
48 scs_int m) {
49 scs_int i, count = 0;
50 /* f cone */
51 for (i = 0; i < k->z; ++i) {
52 /* set rho_y small for z, similar to rho_x term, since z corresponds to
53 * dual free cone, this effectively decreases penalty on those entries
54 * and lets them be determined almost entirely by the linear system solve
55 */
56 rho_y_vec[i] = 1.0 / (1000. * scale);
57 }
58 count += k->z;
59 /* others */
60 for (i = count; i < m; ++i) {
61 rho_y_vec[i] = 1.0 / scale;
62 }
63
64 /* Note, if updating this to use different scales for other cones (e.g. box)
65 * then you must be careful to also include the effect of the rho_y_vec
66 * in the cone projection operator.
67 */
68
69 /* Increase rho_y_vec for the t term in the box cone */
70 if (k->bsize) {
71 rho_y_vec[k->z + k->l] *= BOX_T_SCALE;
72 }
73 }
74
get_sd_cone_size(scs_int s)75 static inline scs_int get_sd_cone_size(scs_int s) {
76 return (s * (s + 1)) / 2;
77 }
78
79 /*
80 * boundaries will contain array of indices of rows of A corresponding to
81 * cone boundaries, boundaries[0] is starting index for cones of size strictly
82 * larger than 1, boundaries malloc-ed here so should be freed.
83 */
SCS(set_cone_boundaries)84 scs_int SCS(set_cone_boundaries)(const ScsCone *k, scs_int **cone_boundaries) {
85 scs_int i, s_cone_sz, count = 0;
86 scs_int cone_boundaries_len =
87 1 + k->qsize + k->ssize + k->ed + k->ep + k->psize;
88 scs_int *b = (scs_int *)scs_calloc(cone_boundaries_len, sizeof(scs_int));
89 /* cones that can be scaled independently */
90 b[count] = k->z + k->l + k->bsize;
91 count += 1; /* started at 0 now move to first entry */
92 for (i = 0; i < k->qsize; ++i) {
93 b[count + i] = k->q[i];
94 }
95 count += k->qsize;
96 for (i = 0; i < k->ssize; ++i) {
97 s_cone_sz = get_sd_cone_size(k->s[i]);
98 b[count + i] = s_cone_sz;
99 }
100 count += k->ssize; /* add ssize here not ssize * (ssize + 1) / 2 */
101 /* exp cones */
102 for (i = 0; i < k->ep + k->ed; ++i) {
103 b[count + i] = 3;
104 }
105 count += k->ep + k->ed;
106 /* power cones */
107 for (i = 0; i < k->psize; ++i) {
108 b[count + i] = 3;
109 }
110 count += k->psize;
111 /* other cones */
112 *cone_boundaries = b;
113 return cone_boundaries_len;
114 }
115
get_full_cone_dims(const ScsCone * k)116 static scs_int get_full_cone_dims(const ScsCone *k) {
117 scs_int i, c = k->z + k->l + k->bsize;
118 if (k->qsize) {
119 for (i = 0; i < k->qsize; ++i) {
120 c += k->q[i];
121 }
122 }
123 if (k->ssize) {
124 for (i = 0; i < k->ssize; ++i) {
125 c += get_sd_cone_size(k->s[i]);
126 }
127 }
128 if (k->ed) {
129 c += 3 * k->ed;
130 }
131 if (k->ep) {
132 c += 3 * k->ep;
133 }
134 if (k->psize) {
135 c += 3 * k->psize;
136 }
137 return c;
138 }
139
SCS(validate_cones)140 scs_int SCS(validate_cones)(const ScsData *d, const ScsCone *k) {
141 scs_int i;
142 if (get_full_cone_dims(k) != d->m) {
143 scs_printf("cone dimensions %li not equal to num rows in A = m = %li\n",
144 (long)get_full_cone_dims(k), (long)d->m);
145 return -1;
146 }
147 if (k->z && k->z < 0) {
148 scs_printf("free cone dimension error\n");
149 return -1;
150 }
151 if (k->l && k->l < 0) {
152 scs_printf("lp cone dimension error\n");
153 return -1;
154 }
155 if (k->bsize) {
156 if (k->bsize < 0) {
157 scs_printf("box cone dimension error\n");
158 return -1;
159 }
160 for (i = 0; i < k->bsize - 1; ++i) {
161 if (k->bl[i] > k->bu[i]) {
162 scs_printf("infeasible: box lower bound larger than upper bound\n");
163 return -1;
164 }
165 }
166 }
167 if (k->qsize && k->q) {
168 if (k->qsize < 0) {
169 scs_printf("soc cone dimension error\n");
170 return -1;
171 }
172 for (i = 0; i < k->qsize; ++i) {
173 if (k->q[i] < 0) {
174 scs_printf("soc cone dimension error\n");
175 return -1;
176 }
177 }
178 }
179 if (k->ssize && k->s) {
180 if (k->ssize < 0) {
181 scs_printf("sd cone dimension error\n");
182 return -1;
183 }
184 for (i = 0; i < k->ssize; ++i) {
185 if (k->s[i] < 0) {
186 scs_printf("sd cone dimension error\n");
187 return -1;
188 }
189 }
190 }
191 if (k->ed && k->ed < 0) {
192 scs_printf("ep cone dimension error\n");
193 return -1;
194 }
195 if (k->ep && k->ep < 0) {
196 scs_printf("ed cone dimension error\n");
197 return -1;
198 }
199 if (k->psize && k->p) {
200 if (k->psize < 0) {
201 scs_printf("power cone dimension error\n");
202 return -1;
203 }
204 for (i = 0; i < k->psize; ++i) {
205 if (k->p[i] < -1 || k->p[i] > 1) {
206 scs_printf("power cone error, values must be in [-1,1]\n");
207 return -1;
208 }
209 }
210 }
211 return 0;
212 }
213
SCS(finish_cone)214 void SCS(finish_cone)(ScsConeWork *c) {
215 #ifdef USE_LAPACK
216 if (c->Xs) {
217 scs_free(c->Xs);
218 }
219 if (c->Z) {
220 scs_free(c->Z);
221 }
222 if (c->e) {
223 scs_free(c->e);
224 }
225 if (c->work) {
226 scs_free(c->work);
227 }
228 #endif
229 if (c->s) {
230 scs_free(c->s);
231 }
232 if (c->bu) {
233 scs_free(c->bu);
234 }
235 if (c->bl) {
236 scs_free(c->bl);
237 }
238 if (c) {
239 scs_free(c);
240 }
241 }
242
SCS(get_cone_header)243 char *SCS(get_cone_header)(const ScsCone *k) {
244 char *tmp = (char *)scs_malloc(sizeof(char) * 512);
245 scs_int i, soc_vars, sd_vars;
246 sprintf(tmp, "cones: ");
247 if (k->z) {
248 sprintf(tmp + strlen(tmp), "\t z: primal zero / dual free vars: %li\n",
249 (long)k->z);
250 }
251 if (k->l) {
252 sprintf(tmp + strlen(tmp), "\t l: linear vars: %li\n", (long)k->l);
253 }
254 if (k->bsize) {
255 sprintf(tmp + strlen(tmp), "\t b: box cone vars: %li\n", (long)(k->bsize));
256 }
257 soc_vars = 0;
258 if (k->qsize && k->q) {
259 for (i = 0; i < k->qsize; i++) {
260 soc_vars += k->q[i];
261 }
262 sprintf(tmp + strlen(tmp), "\t q: soc vars: %li, qsize: %li\n",
263 (long)soc_vars, (long)k->qsize);
264 }
265 sd_vars = 0;
266 if (k->ssize && k->s) {
267 for (i = 0; i < k->ssize; i++) {
268 sd_vars += get_sd_cone_size(k->s[i]);
269 }
270 sprintf(tmp + strlen(tmp), "\t s: psd vars: %li, ssize: %li\n",
271 (long)sd_vars, (long)k->ssize);
272 }
273 if (k->ep || k->ed) {
274 sprintf(tmp + strlen(tmp), "\t e: exp vars: %li, dual exp vars: %li\n",
275 (long)(3 * k->ep), (long)(3 * k->ed));
276 }
277 if (k->psize && k->p) {
278 sprintf(tmp + strlen(tmp), "\t p: primal + dual power vars: %li\n",
279 (long)(3 * k->psize));
280 }
281 return tmp;
282 }
283
exp_newton_one_d(scs_float rho,scs_float y_hat,scs_float z_hat,scs_float w)284 static scs_float exp_newton_one_d(scs_float rho, scs_float y_hat,
285 scs_float z_hat, scs_float w) {
286 scs_float t_prev, t = MAX(w - z_hat, MAX(-z_hat, 1e-9));
287 scs_float f = 1., fp = 1.;
288 scs_int i;
289 for (i = 0; i < EXP_CONE_MAX_ITERS; ++i) {
290 t_prev = t;
291 f = t * (t + z_hat) / rho / rho - y_hat / rho + log(t / rho) + 1;
292 fp = (2 * t + z_hat) / rho / rho + 1 / t;
293
294 t = t - f / fp;
295
296 if (t <= -z_hat) {
297 t = -z_hat;
298 break;
299 } else if (t <= 0) {
300 t = 0;
301 break;
302 } else if (ABS(t - t_prev) < CONE_TOL) {
303 break;
304 } else if (SQRTF(f * f / fp) < CONE_TOL) {
305 break;
306 }
307 }
308 if (i == EXP_CONE_MAX_ITERS) {
309 scs_printf("warning: exp cone newton step hit maximum %i iters\n", (int)i);
310 scs_printf("rho=%1.5e; y_hat=%1.5e; z_hat=%1.5e; w=%1.5e; f=%1.5e, "
311 "fp=%1.5e, t=%1.5e, t_prev= %1.5e\n",
312 rho, y_hat, z_hat, w, f, fp, t, t_prev);
313 }
314 return t + z_hat;
315 }
316
exp_solve_for_x_with_rho(const scs_float * v,scs_float * x,scs_float rho,scs_float w)317 static void exp_solve_for_x_with_rho(const scs_float *v, scs_float *x,
318 scs_float rho, scs_float w) {
319 x[2] = exp_newton_one_d(rho, v[1], v[2], w);
320 x[1] = (x[2] - v[2]) * x[2] / rho;
321 x[0] = v[0] - rho;
322 }
323
exp_calc_grad(const scs_float * v,scs_float * x,scs_float rho,scs_float w)324 static scs_float exp_calc_grad(const scs_float *v, scs_float *x, scs_float rho,
325 scs_float w) {
326 exp_solve_for_x_with_rho(v, x, rho, w);
327 if (x[1] <= 1e-12) {
328 return x[0];
329 }
330 return x[0] + x[1] * log(x[1] / x[2]);
331 }
332
exp_get_rho_ub(const scs_float * v,scs_float * x,scs_float * ub,scs_float * lb)333 static void exp_get_rho_ub(const scs_float *v, scs_float *x, scs_float *ub,
334 scs_float *lb) {
335 *lb = 0;
336 *ub = 0.125;
337 while (exp_calc_grad(v, x, *ub, v[1]) > 0) {
338 *lb = *ub;
339 (*ub) *= 2;
340 }
341 }
342
343 /* project onto the exponential cone, v has dimension *exactly* 3 */
proj_exp_cone(scs_float * v)344 static scs_int proj_exp_cone(scs_float *v) {
345 scs_int i;
346 scs_float ub, lb, rho, g, x[3];
347 scs_float r = v[0], s = v[1], t = v[2];
348
349 /* v in cl(Kexp) */
350 if ((s * exp(r / s) - t <= CONE_THRESH && s > 0) ||
351 (r <= 0 && s == 0 && t >= 0)) {
352 return 0;
353 }
354
355 /* -v in Kexp^* */
356 if ((r > 0 && r * exp(s / r) + exp(1) * t <= CONE_THRESH) ||
357 (r == 0 && s <= 0 && t <= 0)) {
358 memset(v, 0, 3 * sizeof(scs_float));
359 return 0;
360 }
361
362 /* special case with analytical solution */
363 if (r < 0 && s < 0) {
364 v[1] = 0.0;
365 v[2] = MAX(v[2], 0);
366 return 0;
367 }
368
369 /* iterative procedure to find projection, bisects on dual variable: */
370 exp_get_rho_ub(v, x, &ub, &lb); /* get starting upper and lower bounds */
371 for (i = 0; i < EXP_CONE_MAX_ITERS; ++i) {
372 rho = (ub + lb) / 2; /* halfway between upper and lower bounds */
373 g = exp_calc_grad(v, x, rho, x[1]); /* calculates gradient wrt dual var */
374 if (g > 0) {
375 lb = rho;
376 } else {
377 ub = rho;
378 }
379 if (ub - lb < CONE_TOL) {
380 break;
381 }
382 }
383 #if VERBOSITY > 10
384 scs_printf("exponential cone proj iters %i\n", (int)i);
385 #endif
386 if (i == EXP_CONE_MAX_ITERS) {
387 scs_printf("warning: exp cone outer step hit maximum %i iters\n", (int)i);
388 scs_printf("r=%1.5e; s=%1.5e; t=%1.5e\n", r, s, t);
389 }
390 v[0] = x[0];
391 v[1] = x[1];
392 v[2] = x[2];
393 return 0;
394 }
395
set_up_sd_cone_work_space(ScsConeWork * c,const ScsCone * k)396 static scs_int set_up_sd_cone_work_space(ScsConeWork *c, const ScsCone *k) {
397 scs_int i;
398 #ifdef USE_LAPACK
399 blas_int n_max = 0;
400 blas_int neg_one = -1;
401 blas_int info = 0;
402 scs_float wkopt = 0.0;
403 #if VERBOSITY > 0
404 #define _STR_EXPAND(tok) #tok
405 #define _STR(tok) _STR_EXPAND(tok)
406 scs_printf("BLAS(func) = '%s'\n", _STR(BLAS(func)));
407 #endif
408 /* eigenvector decomp workspace */
409 for (i = 0; i < k->ssize; ++i) {
410 if (k->s[i] > n_max) {
411 n_max = (blas_int)k->s[i];
412 }
413 }
414 c->Xs = (scs_float *)scs_calloc(n_max * n_max, sizeof(scs_float));
415 c->Z = (scs_float *)scs_calloc(n_max * n_max, sizeof(scs_float));
416 c->e = (scs_float *)scs_calloc(n_max, sizeof(scs_float));
417
418 /* workspace query */
419 BLAS(syev)
420 ("Vectors", "Lower", &n_max, c->Xs, &n_max, SCS_NULL, &wkopt, &neg_one,
421 &info);
422
423 if (info != 0) {
424 scs_printf("FATAL: syev failure, info = %li\n", (long)info);
425 return -1;
426 }
427 c->lwork = (blas_int)(wkopt + 1); /* +1 for int casting safety */
428 c->work = (scs_float *)scs_calloc(c->lwork, sizeof(scs_float));
429
430 if (!c->Xs || !c->Z || !c->e || !c->work) {
431 return -1;
432 }
433 return 0;
434 #else
435 for (i = 0; i < k->ssize; i++) {
436 if (k->s[i] > 1) {
437 scs_printf(
438 "FATAL: Cannot solve SDPs without linked blas+lapack libraries\n");
439 scs_printf(
440 "Install blas+lapack and re-compile SCS with blas+lapack library "
441 "locations\n");
442 return -1;
443 }
444 }
445 return 0;
446 #endif
447 }
448
449 /* size of X is get_sd_cone_size(n) */
proj_semi_definite_cone(scs_float * X,const scs_int n,ScsConeWork * c)450 static scs_int proj_semi_definite_cone(scs_float *X, const scs_int n,
451 ScsConeWork *c) {
452 /* project onto the positive semi-definite cone */
453 #ifdef USE_LAPACK
454 scs_int i, first_idx;
455 blas_int nb = (blas_int)n;
456 blas_int ncols_z;
457 blas_int nb_plus_one = (blas_int)(n + 1);
458 blas_int one_int = 1;
459 scs_float zero = 0., one = 1.;
460 scs_float sqrt2 = SQRTF(2.0);
461 scs_float sqrt2_inv = 1.0 / sqrt2;
462 scs_float *Xs = c->Xs;
463 scs_float *Z = c->Z;
464 scs_float *e = c->e;
465 scs_float *work = c->work;
466 blas_int lwork = c->lwork;
467 blas_int info = 0;
468 scs_float sq_eig_pos;
469
470 #endif
471
472 if (n == 0) {
473 return 0;
474 }
475 if (n == 1) {
476 X[0] = MAX(X[0], 0.);
477 return 0;
478 }
479
480 #ifdef USE_LAPACK
481
482 /* copy lower triangular matrix into full matrix */
483 for (i = 0; i < n; ++i) {
484 memcpy(&(Xs[i * (n + 1)]), &(X[i * n - ((i - 1) * i) / 2]),
485 (n - i) * sizeof(scs_float));
486 }
487 /*
488 rescale so projection works, and matrix norm preserved
489 see http://www.seas.ucla.edu/~vandenbe/publications/mlbook.pdf pg 3
490 */
491 /* scale diags by sqrt(2) */
492 BLAS(scal)(&nb, &sqrt2, Xs, &nb_plus_one); /* not n_squared */
493
494 /* Solve eigenproblem, reuse workspaces */
495 BLAS(syev)("Vectors", "Lower", &nb, Xs, &nb, e, work, &lwork, &info);
496 if (info != 0) {
497 scs_printf("WARN: LAPACK syev error, info = %i\n", (int)info);
498 if (info < 0) {
499 return info;
500 }
501 }
502
503 first_idx = -1;
504 /* e is eigvals in ascending order, find first entry > 0 */
505 for (i = 0; i < n; ++i) {
506 if (e[i] > 0) {
507 first_idx = i;
508 break;
509 }
510 }
511
512 if (first_idx == -1) {
513 /* there are no positive eigenvalues, set X to 0 and return */
514 memset(X, 0, sizeof(scs_float) * get_sd_cone_size(n));
515 return 0;
516 }
517
518 /* Z is matrix of eigenvectors with positive eigenvalues */
519 memcpy(Z, &Xs[first_idx * n], sizeof(scs_float) * n * (n - first_idx));
520
521 /* scale Z by sqrt(eig) */
522 for (i = first_idx; i < n; ++i) {
523 sq_eig_pos = SQRTF(e[i]);
524 BLAS(scal)(&nb, &sq_eig_pos, &Z[(i - first_idx) * n], &one_int);
525 }
526
527 /* Xs = Z Z' = V E V' */
528 ncols_z = (blas_int)(n - first_idx);
529 BLAS(syrk)("Lower", "NoTrans", &nb, &ncols_z, &one, Z, &nb, &zero, Xs, &nb);
530
531 /* undo rescaling: scale diags by 1/sqrt(2) */
532 BLAS(scal)(&nb, &sqrt2_inv, Xs, &nb_plus_one); /* not n_squared */
533
534 /* extract just lower triangular matrix */
535 for (i = 0; i < n; ++i) {
536 memcpy(&(X[i * n - ((i - 1) * i) / 2]), &(Xs[i * (n + 1)]),
537 (n - i) * sizeof(scs_float));
538 }
539 return 0;
540
541 #else
542 scs_printf("FAILURE: solving SDP but no blas/lapack libraries were found!\n");
543 scs_printf("SCS will return nonsense!\n");
544 SCS(scale_array)(X, NAN, n);
545 return -1;
546 #endif
547 }
548
pow_calc_x(scs_float r,scs_float xh,scs_float rh,scs_float a)549 static scs_float pow_calc_x(scs_float r, scs_float xh, scs_float rh,
550 scs_float a) {
551 scs_float x = 0.5 * (xh + SQRTF(xh * xh + 4 * a * (rh - r) * r));
552 return MAX(x, 1e-12);
553 }
554
pow_calcdxdr(scs_float x,scs_float xh,scs_float rh,scs_float r,scs_float a)555 static scs_float pow_calcdxdr(scs_float x, scs_float xh, scs_float rh,
556 scs_float r, scs_float a) {
557 return a * (rh - 2 * r) / (2 * x - xh);
558 }
559
pow_calc_f(scs_float x,scs_float y,scs_float r,scs_float a)560 static scs_float pow_calc_f(scs_float x, scs_float y, scs_float r,
561 scs_float a) {
562 return POWF(x, a) * POWF(y, (1 - a)) - r;
563 }
564
pow_calc_fp(scs_float x,scs_float y,scs_float dxdr,scs_float dydr,scs_float a)565 static scs_float pow_calc_fp(scs_float x, scs_float y, scs_float dxdr,
566 scs_float dydr, scs_float a) {
567 return POWF(x, a) * POWF(y, (1 - a)) * (a * dxdr / x + (1 - a) * dydr / y) -
568 1;
569 }
570
571 /*
572 * Routine to scale the limits of the box cone by the scaling diagonal mat D > 0
573 *
574 * want (t, s) \in K <==> (t', s') \in K'
575 *
576 * (t', s') = (d0 * t, D s) (overloading D to mean D[1:])
577 * (up to scalar scaling factor which we can ignore due to conic prooperty)
578 *
579 * K = { (t, s) | t * l <= s <= t * u, t >= 0 } =>
580 * { (t, s) | d0 * t * D l / d0 <= D s <= d0 * t D u / d0, t >= 0 } =>
581 * { (t', s') | t' * l' <= s' <= t' u', t >= 0 } = K'
582 * where l' = D l / d0, u' = D u / d0.
583 */
normalize_box_cone(ScsConeWork * c,scs_float * D,scs_int bsize)584 static void normalize_box_cone(ScsConeWork *c, scs_float *D, scs_int bsize) {
585 scs_int j;
586 for (j = 0; j < bsize - 1; j++) {
587 if (c->bu[j] >= MAX_BOX_VAL) {
588 c->bu[j] = INFINITY;
589 } else {
590 c->bu[j] = D ? D[j + 1] * c->bu[j] / D[0] : c->bu[j];
591 }
592 if (c->bl[j] <= -MAX_BOX_VAL) {
593 c->bl[j] = -INFINITY;
594 } else {
595 c->bl[j] = D ? D[j + 1] * c->bl[j] / D[0] : c->bl[j];
596 }
597 }
598 }
599
600 /* project onto { (t, s) | t * l <= s <= t * u, t >= 0 }, Newton's method on t
601 tx = [t; s], total length = bsize
602 uses Moreau since \Pi_K*(tx) = \Pi_K(-tx) + tx
603 */
proj_box_cone(scs_float * tx,const scs_float * bl,const scs_float * bu,scs_int bsize,scs_float t_warm_start)604 static scs_float proj_box_cone(scs_float *tx, const scs_float *bl,
605 const scs_float *bu, scs_int bsize,
606 scs_float t_warm_start) {
607 scs_float *x, gt, ht, t_prev, t = t_warm_start;
608 scs_int iter, j;
609
610 if (bsize == 1) { /* special case */
611 tx[0] = MAX(tx[0], 0.0);
612 return tx[0];
613 }
614
615 x = &(tx[1]);
616
617 /* should only require about 5 or so iterations, 1 or 2 if warm-started */
618 for (iter = 0; iter < BOX_CONE_MAX_ITERS; iter++) {
619 t_prev = t;
620 /* incorporate the additional BOX_T_SCALE factor into the projection */
621 gt = BOX_T_SCALE * (t - tx[0]); /* gradient */
622 ht = BOX_T_SCALE; /* hessian */
623 for (j = 0; j < bsize - 1; j++) {
624 if (x[j] > t * bu[j]) {
625 gt += (t * bu[j] - x[j]) * bu[j]; /* gradient */
626 ht += bu[j] * bu[j]; /* hessian */
627 } else if (x[j] < t * bl[j]) {
628 gt += (t * bl[j] - x[j]) * bl[j]; /* gradient */
629 ht += bl[j] * bl[j]; /* hessian */
630 }
631 }
632 t = MAX(t - gt / MAX(ht, 1e-8), 0.); /* newton step */
633 #if VERBOSITY > 3
634 scs_printf("iter %i, t_new %1.3e, t_prev %1.3e, gt %1.3e, ht %1.3e\n", iter,
635 t, t_prev, gt, ht);
636 scs_printf("ABS(gt / (ht + 1e-6)) %.4e, ABS(t - t_prev) %.4e\n",
637 ABS(gt / (ht + 1e-6)), ABS(t - t_prev));
638 #endif
639 /* TODO: sometimes this check can fail (ie, declare convergence before it
640 * should) if ht is very large, which can happen with some pathological
641 * problems.
642 */
643 if (ABS(gt / MAX(ht, 1e-6)) < 1e-12 * MAX(t, 1.) ||
644 ABS(t - t_prev) < 1e-11 * MAX(t, 1.)) {
645 break;
646 }
647 }
648 if (iter == BOX_CONE_MAX_ITERS) {
649 scs_printf("warning: box cone proj hit maximum %i iters\n", (int)iter);
650 }
651 for (j = 0; j < bsize - 1; j++) {
652 if (x[j] > t * bu[j]) {
653 x[j] = t * bu[j];
654 } else if (x[j] < t * bl[j]) {
655 x[j] = t * bl[j];
656 }
657 /* x[j] unchanged otherwise */
658 }
659 tx[0] = t;
660 #if VERBOSITY > 3
661 scs_printf("box cone iters %i\n", (int)iter + 1);
662 #endif
663 return t;
664 }
665
666 /* project onto SOC of size q*/
proj_soc(scs_float * x,scs_int q)667 static void proj_soc(scs_float *x, scs_int q) {
668 if (q == 0) {
669 return;
670 }
671 if (q == 1) {
672 x[0] = MAX(x[0], 0.);
673 return;
674 }
675 scs_float v1 = x[0];
676 scs_float s = SCS(norm_2)(&(x[1]), q - 1);
677 scs_float alpha = (s + v1) / 2.0;
678
679 if (s <= v1) {
680 return;
681 } else if (s <= -v1) {
682 memset(&(x[0]), 0, q * sizeof(scs_float));
683 } else {
684 x[0] = alpha;
685 SCS(scale_array)(&(x[1]), alpha / s, q - 1);
686 }
687 }
688
proj_power_cone(scs_float * v,scs_float a)689 static void proj_power_cone(scs_float *v, scs_float a) {
690 scs_float xh = v[0], yh = v[1], rh = ABS(v[2]);
691 scs_float x = 0.0, y = 0.0, r;
692 scs_int i;
693 /* v in K_a */
694 if (xh >= 0 && yh >= 0 &&
695 CONE_THRESH + POWF(xh, a) * POWF(yh, (1 - a)) >= rh) {
696 return;
697 }
698
699 /* -v in K_a^* */
700 if (xh <= 0 && yh <= 0 &&
701 CONE_THRESH + POWF(-xh, a) * POWF(-yh, 1 - a) >=
702 rh * POWF(a, a) * POWF(1 - a, 1 - a)) {
703 v[0] = v[1] = v[2] = 0;
704 return;
705 }
706
707 r = rh / 2;
708 for (i = 0; i < POW_CONE_MAX_ITERS; ++i) {
709 scs_float f, fp, dxdr, dydr;
710 x = pow_calc_x(r, xh, rh, a);
711 y = pow_calc_x(r, yh, rh, 1 - a);
712
713 f = pow_calc_f(x, y, r, a);
714 if (ABS(f) < CONE_TOL) {
715 break;
716 }
717
718 dxdr = pow_calcdxdr(x, xh, rh, r, a);
719 dydr = pow_calcdxdr(y, yh, rh, r, (1 - a));
720 fp = pow_calc_fp(x, y, dxdr, dydr, a);
721
722 r = MAX(r - f / fp, 0);
723 r = MIN(r, rh);
724 }
725 v[0] = x;
726 v[1] = y;
727 v[2] = (v[2] < 0) ? -(r) : (r);
728 }
729
730 /* project onto the primal K cone in the paper */
proj_cone(scs_float * x,const ScsCone * k,ScsConeWork * c,scs_int normalize)731 static scs_int proj_cone(scs_float *x, const ScsCone *k, ScsConeWork *c,
732 scs_int normalize) {
733 scs_int i, status;
734 scs_int count = 0;
735
736 if (k->z) {
737 /* project onto primal zero / dual free cone */
738 memset(x, 0, k->z * sizeof(scs_float));
739 count += k->z;
740 }
741
742 if (k->l) {
743 /* project onto positive orthant */
744 for (i = count; i < count + k->l; ++i) {
745 x[i] = MAX(x[i], 0.0);
746 }
747 count += k->l;
748 }
749
750 if (k->bsize) {
751 /* project onto box cone */
752 if (normalize) {
753 c->box_t_warm_start = proj_box_cone(&(x[count]), c->bl, c->bu, k->bsize,
754 c->box_t_warm_start);
755 } else {
756 c->box_t_warm_start = proj_box_cone(&(x[count]), k->bl, k->bu, k->bsize,
757 c->box_t_warm_start);
758 }
759 count += k->bsize; /* since b = (t,s), len(s) = bsize - 1 */
760 }
761
762 if (k->qsize && k->q) {
763 /* project onto second-order cones */
764 for (i = 0; i < k->qsize; ++i) {
765 proj_soc(&(x[count]), k->q[i]);
766 count += k->q[i];
767 }
768 }
769
770 if (k->ssize && k->s) {
771 /* project onto PSD cones */
772 for (i = 0; i < k->ssize; ++i) {
773 status = proj_semi_definite_cone(&(x[count]), k->s[i], c);
774 if (status < 0) {
775 return status;
776 }
777 count += get_sd_cone_size(k->s[i]);
778 }
779 }
780
781 if (k->ep) {
782 /*
783 * exponential cone is not self dual, if s \in K
784 * then y \in K^* and so if K is the primal cone
785 * here we project onto K^*, via Moreau
786 * \Pi_C^*(y) = y + \Pi_C(-y)
787 */
788 #ifdef _OPENMP
789 #pragma omp parallel for
790 #endif
791 for (i = 0; i < k->ep; ++i) {
792 proj_exp_cone(&(x[count + 3 * i]));
793 }
794 count += 3 * k->ep;
795 }
796
797 if (k->ed) { /* dual exponential cone */
798 /*
799 * exponential cone is not self dual, if s \in K
800 * then y \in K^* and so if K is the primal cone
801 * here we project onto K^*, via Moreau
802 * \Pi_C^*(y) = y + \Pi_C(-y)
803 */
804 scs_int idx;
805 scs_float r, s, t;
806 SCS(scale_array)(&(x[count]), -1, 3 * k->ed); /* x = -x; */
807 #ifdef _OPENMP
808 #pragma omp parallel for private(r, s, t, idx)
809 #endif
810 for (i = 0; i < k->ed; ++i) {
811 idx = count + 3 * i;
812 r = x[idx];
813 s = x[idx + 1];
814 t = x[idx + 2];
815
816 proj_exp_cone(&(x[idx]));
817
818 x[idx] -= r;
819 x[idx + 1] -= s;
820 x[idx + 2] -= t;
821 }
822 count += 3 * k->ed;
823 }
824
825 if (k->psize && k->p) {
826 scs_float v[3];
827 scs_int idx;
828 /* don't use openmp for power cone
829 ifdef _OPENMP
830 pragma omp parallel for private(v, idx)
831 endif
832 */
833 for (i = 0; i < k->psize; ++i) {
834 idx = count + 3 * i;
835 if (k->p[i] >= 0) {
836 /* primal power cone */
837 proj_power_cone(&(x[idx]), k->p[i]);
838 } else {
839 /* dual power cone, using Moreau */
840 v[0] = -x[idx];
841 v[1] = -x[idx + 1];
842 v[2] = -x[idx + 2];
843
844 proj_power_cone(v, -k->p[i]);
845
846 x[idx] += v[0];
847 x[idx + 1] += v[1];
848 x[idx + 2] += v[2];
849 }
850 }
851 count += 3 * k->psize;
852 }
853 /* project onto OTHER cones */
854 return 0;
855 }
856
SCS(init_cone)857 ScsConeWork *SCS(init_cone)(const ScsCone *k, const ScsScaling *scal,
858 scs_int cone_len) {
859 ScsConeWork *c = (ScsConeWork *)scs_calloc(1, sizeof(ScsConeWork));
860 c->cone_len = cone_len;
861 c->s = (scs_float *)scs_calloc(cone_len, sizeof(scs_float));
862 if (k->bsize && k->bu && k->bl) {
863 c->box_t_warm_start = 1.;
864 if (scal) {
865 c->bu = (scs_float *)scs_calloc(k->bsize - 1, sizeof(scs_float));
866 c->bl = (scs_float *)scs_calloc(k->bsize - 1, sizeof(scs_float));
867 memcpy(c->bu, k->bu, (k->bsize - 1) * sizeof(scs_float));
868 memcpy(c->bl, k->bl, (k->bsize - 1) * sizeof(scs_float));
869 /* also does some sanitizing */
870 normalize_box_cone(c, scal ? &(scal->D[k->z + k->l]) : SCS_NULL,
871 k->bsize);
872 }
873 }
874 if (k->ssize && k->s) {
875 if (set_up_sd_cone_work_space(c, k) < 0) {
876 SCS(finish_cone)(c);
877 return SCS_NULL;
878 }
879 }
880 return c;
881 }
882
883 /* outward facing cone projection routine
884 performs projection in-place
885 if normalize > 0 then will use normalized (equilibrated) cones if applicable.
886 */
SCS(proj_dual_cone)887 scs_int SCS(proj_dual_cone)(scs_float *x, const ScsCone *k, ScsConeWork *c,
888 scs_int normalize) {
889 scs_int status;
890 /* copy x, s = x */
891 memcpy(c->s, x, c->cone_len * sizeof(scs_float));
892 /* negate x -> -x */
893 SCS(scale_array)(x, -1., c->cone_len);
894 /* project -x onto cone, x -> Pi_K(-x) */
895 status = proj_cone(x, k, c, normalize);
896 /* return Pi_K*(x) = s + Pi_K(-x) */
897 SCS(add_scaled_array)(x, c->s, c->cone_len, 1.);
898 return status;
899 }
900