Scale the variables with variable bounds

benchmarks/linear_programming/run_mps_files.sh

@@ -360,7 +360,7 @@ else
         mapfile -t mps_files < <(find "$MPS_DIR" -type f \( -name "*.mps" -o -name "*.MPS" -o -name "*.SIF" \) | sort)
     else
         # Gather .mps/.MPS and .SIF files in the directory
-        mapfile -t mps_files < <(ls "$MPS_DIR"/*.mps "$MPS_DIR"/*.MPS "$MPS_DIR"/*.SIF 2>/dev/null)
+        mapfile -t mps_files < <(ls "$MPS_DIR"/*.mps "$MPS_DIR"/*.MPS "$MPS_DIR"/*.SIF "$MPS_DIR"/*.mps.gz "$MPS_DIR"/*.MPS.gz "$MPS_DIR"/*.SIF.gz 2>/dev/null)
     fi
 
     echo "Found ${#mps_files[@]} .mps and .SIF files in $MPS_DIR"

cpp/src/dual_simplex/scaling.cpp

@@ -115,6 +115,78 @@ i_t scaling(const lp_problem_t<i_t, f_t>& unscaled,
         q_ratio);
     }
 
+    // -----------------------------------------------------------------------
+    // Bound-magnitude column pre-scaling.
+    // -----------------------------------------------------------------------
+    // Ruiz equilibration only balances the *coefficient* magnitudes of A and Q.
+    // On problems whose variable BOUNDS span many orders of magnitude (e.g. a
+    // network-flow QP with |bound| ranging from ~1e4 to ~1e11), the constraint
+    // matrix can already be perfectly balanced (all +/-1) while the variables
+    // themselves live at wildly different scales. The interior-point diagonal
+    // D = z/x then spans those same many orders of magnitude, wrecking the
+    // conditioning of the KKT factorization and stalling convergence.
+    //
+    // We fix this by first scaling each column so that the variable it
+    // represents becomes O(1): c0[j] = (geometric mean of the finite bound
+    // magnitudes of x_j). Substituting x_j = c0[j] * x'_j leaves the feasible
+    // region shape unchanged but brings every variable to a common scale, which
+    // Ruiz then finishes off on the coefficient side. Columns with no finite,
+    // nonzero bound are left at scale 1.
+    {
+      std::vector<f_t> c0(n, 1.0);
+      const i_t cone_start0 =
+        unscaled.second_order_cone_dims.empty() ? n : unscaled.cone_var_start;
+      f_t geo_sum = 0.0;
+      i_t geo_count = 0;
+      for (i_t j = 0; j < cone_start0; ++j) {
+        f_t lo = std::abs(scaled.lower[j]);
+        f_t hi = std::abs(scaled.upper[j]);
+        f_t mag;
+        bool lo_fin = scaled.lower[j] > -1e20 && lo > 0;
+        bool hi_fin = scaled.upper[j] < 1e20 && hi > 0;
+        if (lo_fin && hi_fin) {
+          mag = std::sqrt(lo * hi);
+        } else if (lo_fin) {
+          mag = lo;
+        } else if (hi_fin) {
+          mag = hi;
+        } else {
+          continue;  // free / one-sided-zero: leave at scale 1
+        }
+        c0[j] = mag;
+        geo_sum += std::log(mag);
+        geo_count++;
+      }
+      if (geo_count > 0) {
+        // Normalize so the average column scale is 1, keeping the overall
+        // problem magnitude centered rather than uniformly shrinking it.
+        const f_t geo_mean = std::exp(geo_sum / static_cast<f_t>(geo_count));
+        for (i_t j = 0; j < cone_start0; ++j) {
+          c0[j] /= geo_mean;
+        }
+        // Apply x_j = c0[j] * x'_j : A(:,j) *= c0[j], obj[j] *= c0[j],
+        // bounds /= c0[j], Q(i,j) *= c0[i]*c0[j], accumulate into col_scale.
+        for (i_t j = 0; j < n; ++j) {
+          if (c0[j] == 1.0) continue;
+          for (i_t p = scaled.A.col_start[j]; p < scaled.A.col_start[j + 1]; ++p) {
+            scaled.A.x[p] *= c0[j];
+          }
+          scaled.objective[j] *= c0[j];
+          if (scaled.lower[j] > -1e20) scaled.lower[j] /= c0[j];
+          if (scaled.upper[j] < 1e20) scaled.upper[j] /= c0[j];
+          col_scale[j] *= c0[j];
+        }
+        if (scaled.Q.n > 0) {
+          for (i_t row = 0; row < scaled.Q.m; ++row) {
+            for (i_t p = scaled.Q.row_start[row]; p < scaled.Q.row_start[row + 1]; ++p) {
+              i_t col = scaled.Q.j[p];
+              scaled.Q.x[p] *= c0[row] * c0[col];
+            }
+          }
+        }
+      }
+    }
+
     // Apply Ruiz equilibration
     csr_matrix_t<i_t, f_t> Arow(0, 0, 0);
     scaled.A.to_compressed_row(Arow);