Ruiz clamping

benchmarks/benchmark_runner.c

@@ -60,6 +60,12 @@ static void apply_setting(QOCOSettings* settings, const char* arg)
   else if (strcmp(key, "ruiz_iters") == 0) {
     settings->ruiz_iters = atoi(val);
   }
+  else if (strcmp(key, "ruiz_scaling_min") == 0) {
+    settings->ruiz_scaling_min = atof(val);
+  }
+  else if (strcmp(key, "ruiz_scaling_max") == 0) {
+    settings->ruiz_scaling_max = atof(val);
+  }
   else if (strcmp(key, "max_ir_iters") == 0) {
     settings->max_ir_iters = atoi(val);
   }

docs/qoco/api/settings.rst

@@ -14,6 +14,10 @@ The settings are defined in the :code:`include/structs.h` file.
 +--------------------------------+-----------------------+------------------------------------------------------------+---------------------+---------------+
 | :code:`ruiz_iters`             |  :code:`QOCOInt`      | Number of Ruiz equilibration iterations performed          | :math:`(0, \infty)` | 0             |
 +--------------------------------+-----------------------+------------------------------------------------------------+---------------------+---------------+
+| :code:`ruiz_scaling_min`       |  :code:`QOCOFloat`    | Minimum cumulative Ruiz scaling                            | :math:`(0, \infty)` | 1e-4          |
++--------------------------------+-----------------------+------------------------------------------------------------+---------------------+---------------+
+| :code:`ruiz_scaling_max`       |  :code:`QOCOFloat`    | Maximum cumulative Ruiz scaling                            | :math:`(0, \infty)` | 1e4           |
++--------------------------------+-----------------------+------------------------------------------------------------+---------------------+---------------+
 | :code:`max_ir_iters`           |  :code:`QOCOInt`      | Maximum number of iterative refinement iterations          | :math:`(0, \infty)` | 5             |
 +--------------------------------+-----------------------+------------------------------------------------------------+---------------------+---------------+
 | :code:`ir_tol`                 |  :code:`QOCOFloat`    | Iterative refinement stopping tolerance: stop when         | :math:`(0, \infty)` | 1e-6          |

include/equilibration.h

@@ -39,7 +39,8 @@
 
  */
 void ruiz_equilibration(QOCOProblemData* data, QOCOScaling* scaling,
-                        QOCOInt ruiz_iters);
+                        QOCOInt ruiz_iters, QOCOFloat ruiz_scaling_min,
+                        QOCOFloat ruiz_scaling_max);
 
 /**
  * @brief Undo variable transformation induced by ruiz equilibration.

include/structs.h

@@ -111,6 +111,12 @@ typedef struct {
   /** Number of Ruiz equilibration iterations. */
   QOCOInt ruiz_iters;
 
+  /** Minimum cumulative Ruiz scaling. */
+  QOCOFloat ruiz_scaling_min;
+
+  /** Maximum cumulative Ruiz scaling. */
+  QOCOFloat ruiz_scaling_max;
+
   /** Maximum number of iterative refinement iterations. */
   QOCOInt max_ir_iters;
 

src/equilibration.c

@@ -1,7 +1,37 @@
 #include "equilibration.h"
 
+static QOCOFloat clip_ruiz_scale(QOCOFloat scale, QOCOFloat cumulative,
+                                 QOCOFloat scaling_min, QOCOFloat scaling_max)
+{
+  if (!isfinite(scale) || scale <= 0.0) {
+    scale = 1.0;
+  }
+  if (!isfinite(cumulative) || cumulative <= 0.0) {
+    cumulative = 1.0;
+  }
+
+  QOCOFloat lower = scaling_min / cumulative;
+  QOCOFloat upper = scaling_max / cumulative;
+  if (scale < lower) {
+    return lower;
+  }
+  if (scale > upper) {
+    return upper;
+  }
+  return scale;
+}
+
+static QOCOFloat norm_to_ruiz_scale(QOCOFloat norm)
+{
+  if (!isfinite(norm) || norm <= 0.0) {
+    return 1.0;
+  }
+  return 1.0 / qoco_sqrt(norm);
+}
+
 void ruiz_equilibration(QOCOProblemData* data, QOCOScaling* scaling,
-                        QOCOInt ruiz_iters)
+                        QOCOInt ruiz_iters, QOCOFloat ruiz_scaling_min,
+                        QOCOFloat ruiz_scaling_max)
 {
   // This function runs on the CPU.
   set_cpu_mode(1);
@@ -48,7 +78,8 @@ void ruiz_equilibration(QOCOProblemData* data, QOCOScaling* scaling,
 
     // g = 1 / max(mean(Pinf), norm(c, "inf"));
     g = qoco_max(Pinf_mean, g);
-    g = safe_div(1.0, g);
+    g = (isfinite(g) && g > 0.0) ? 1.0 / g : 1.0;
+    g = clip_ruiz_scale(g, scaling->k, ruiz_scaling_min, ruiz_scaling_max);
     scaling->k *= g;
 
     // Compute column infinity norms of A and G
@@ -76,8 +107,10 @@ void ruiz_equilibration(QOCOProblemData* data, QOCOScaling* scaling,
 
     // d(i) = 1 / sqrt(max([Pinf(i), Atinf(i), Gtinf(i)]));
     for (QOCOInt j = 0; j < data->n; ++j) {
-      QOCOFloat temp = qoco_sqrt(get_element_vectorf(scaling->delta, j));
-      temp = safe_div(1.0, temp);
+      QOCOFloat temp =
+          norm_to_ruiz_scale(get_element_vectorf(scaling->delta, j));
+      temp = clip_ruiz_scale(temp, Druiz_data[j], ruiz_scaling_min,
+                             ruiz_scaling_max);
       set_element_vectorf(scaling->delta, j, temp);
     }
 
@@ -88,9 +121,10 @@ void ruiz_equilibration(QOCOProblemData* data, QOCOScaling* scaling,
       col_inf_norm_matrix(data->At, &delta_data[data->n]);
       // d(i) = 1 / sqrt(Ainf(i));
       for (QOCOInt k = 0; k < data->p; ++k) {
-        QOCOFloat temp =
-            qoco_sqrt(get_element_vectorf(scaling->delta, data->n + k));
-        temp = safe_div(1.0, temp);
+        QOCOFloat temp = norm_to_ruiz_scale(
+            get_element_vectorf(scaling->delta, data->n + k));
+        temp = clip_ruiz_scale(temp, Eruiz_data[k], ruiz_scaling_min,
+                               ruiz_scaling_max);
         set_element_vectorf(scaling->delta, data->n + k, temp);
       }
     }
@@ -102,9 +136,10 @@ void ruiz_equilibration(QOCOProblemData* data, QOCOScaling* scaling,
       col_inf_norm_matrix(data->Gt, &delta_data[data->n + data->p]);
       // d(i) = 1 / sqrt(Ginf(i));
       for (QOCOInt k = 0; k < data->m; ++k) {
-        QOCOFloat temp = qoco_sqrt(
+        QOCOFloat temp = norm_to_ruiz_scale(
             get_element_vectorf(scaling->delta, data->n + data->p + k));
-        temp = safe_div(1.0, temp);
+        temp = clip_ruiz_scale(temp, Fruiz_data[k], ruiz_scaling_min,
+                               ruiz_scaling_max);
         set_element_vectorf(scaling->delta, data->n + data->p + k, temp);
       }
     }

src/input_validation.c

@@ -18,9 +18,21 @@ QOCOInt qoco_validate_settings(const QOCOSettings* settings)
     return QOCO_SETTINGS_VALIDATION_ERROR;
   }
 
-  // ruiz_iters must be positive.
+  // ruiz_iters must be non-negative.
   if (settings->ruiz_iters < 0) {
-    printf("ruiz_iters must be positive.\n");
+    printf("ruiz_iters must be non-negative.\n");
+    return QOCO_SETTINGS_VALIDATION_ERROR;
+  }
+  if (settings->ruiz_scaling_min <= 0) {
+    printf("ruiz_scaling_min must be positive.\n");
+    return QOCO_SETTINGS_VALIDATION_ERROR;
+  }
+  if (settings->ruiz_scaling_max <= 0) {
+    printf("ruiz_scaling_max must be positive.\n");
+    return QOCO_SETTINGS_VALIDATION_ERROR;
+  }
+  if (settings->ruiz_scaling_min > settings->ruiz_scaling_max) {
+    printf("ruiz_scaling_min must be <= ruiz_scaling_max.\n");
     return QOCO_SETTINGS_VALIDATION_ERROR;
   }
 

src/qoco_api.c

@@ -93,7 +93,9 @@ QOCOInt qoco_setup(QOCOSolver* solver, QOCOInt n, QOCOInt m, QOCOInt p,
   // Compute scaling statistics before equilibration and regularization.
   compute_scaling_statistics(data);
 
-  ruiz_equilibration(data, work->scaling, solver->settings->ruiz_iters);
+  ruiz_equilibration(data, work->scaling, solver->settings->ruiz_iters,
+                     solver->settings->ruiz_scaling_min,
+                     solver->settings->ruiz_scaling_max);
 
   // Regularize P.
   set_cpu_mode(1);
@@ -243,11 +245,13 @@ void qoco_set_csc(QOCOCscMatrix* A, QOCOInt m, QOCOInt n, QOCOInt Annz,
 void set_default_settings(QOCOSettings* settings)
 {
   settings->max_iters = 200;
-  settings->ruiz_iters = 0;
+  settings->ruiz_iters = 1;
+  settings->ruiz_scaling_min = 1e-4;
+  settings->ruiz_scaling_max = 1e4;
   settings->max_ir_iters = 5;
   settings->ir_tol = 1e-6;
   settings->kkt_static_reg_P = 1e-13;
-  settings->kkt_static_reg_A = 1e-8;
+  settings->kkt_static_reg_A = 1e-13;
   settings->kkt_static_reg_G = 1e-13;
   settings->kkt_dynamic_reg = 1e-11;
   settings->abstol = 1e-7;
@@ -266,6 +270,8 @@ QOCOInt qoco_update_settings(QOCOSolver* solver,
 
   solver->settings->max_iters = new_settings->max_iters;
   solver->settings->ruiz_iters = new_settings->ruiz_iters;
+  solver->settings->ruiz_scaling_min = new_settings->ruiz_scaling_min;
+  solver->settings->ruiz_scaling_max = new_settings->ruiz_scaling_max;
   solver->settings->max_ir_iters = new_settings->max_ir_iters;
   solver->settings->ir_tol = new_settings->ir_tol;
   solver->settings->kkt_static_reg_P = new_settings->kkt_static_reg_P;
@@ -440,7 +446,9 @@ void qoco_update_matrix_data(QOCOSolver* solver, QOCOFloat* Pxnew,
   compute_scaling_statistics(data);
 
   // Equilibrate new matrix data.
-  ruiz_equilibration(data, scaling, solver->settings->ruiz_iters);
+  ruiz_equilibration(data, scaling, solver->settings->ruiz_iters,
+                     solver->settings->ruiz_scaling_min,
+                     solver->settings->ruiz_scaling_max);
 
   // Regularize P.
   unregularize(Pcsc, -solver->settings->kkt_static_reg_P);
@@ -627,4 +635,4 @@ QOCOInt qoco_cleanup(QOCOSolver* solver)
   qoco_free(solver);
 
   return 1;
-}
+}

src/qoco_utils.c

@@ -179,6 +179,7 @@ void print_header(QOCOSolver* solver)
   printf("|     max_iters: %-3d abstol: %3.2e reltol: %3.2e  |\n", settings->max_iters, (double)settings->abstol, (double)settings->reltol);
   printf("|     abstol_inacc: %3.2e reltol_inacc: %3.2e     |\n", (double)settings->abstol_inacc, (double)settings->reltol_inacc);
   printf("|     kkt_static_reg_P: %3.2e ruiz_iters: %-2d         |\n", (double)settings->kkt_static_reg_P, settings->ruiz_iters);
+  printf("|     ruiz_scale_min: %3.2e max: %3.2e            |\n", (double)settings->ruiz_scaling_min, (double)settings->ruiz_scaling_max);
   printf("|     kkt_static_reg_A: %3.2e max_ir_iters: %-2d       |\n", (double)settings->kkt_static_reg_A, settings->max_ir_iters);
   printf("|     kkt_static_reg_G: %3.2e ir_tol: %3.2e       |\n", (double)settings->kkt_static_reg_G, (double)settings->ir_tol);
   printf("|     kkt_dynamic_reg: %3.2e                         |\n", (double)settings->kkt_dynamic_reg);
@@ -474,6 +475,8 @@ QOCOSettings* copy_settings(QOCOSettings* settings)
   new_settings->reltol = settings->reltol;
   new_settings->reltol_inacc = settings->reltol_inacc;
   new_settings->ruiz_iters = settings->ruiz_iters;
+  new_settings->ruiz_scaling_min = settings->ruiz_scaling_min;
+  new_settings->ruiz_scaling_max = settings->ruiz_scaling_max;
   new_settings->verbose = settings->verbose;
 
   return new_settings;