Clarabel 7.1 errors with insufficient progress on a seemingly simple SOCP

[RussTedrake]

trafficstars

While upgrading Drake to clarabel 7.1 (https://github.com/RobotLocomotion/drake/pull/21053), one of our (seemingly simple) unit tests started falling down with "insufficient progress". Mosek solves it fine.

I've carved out the following reproduction

import clarabel
import numpy as np
from scipy import sparse

q = [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1]
b = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 1, -0, -0, -0, 1, -0, -0, -0, 1, 1, -0, 1, -0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
P = sparse.csc_matrix((18, 18))

data = [-1, 1, 1, -1, 1, -1, -1, -1, -1, 1, -0.5, -0.5, 1, -0.7071067811865476, 1, -0.7071067811865476, 1, -1, 0.7071067811865475, -1.0536712127723509e-08, 1, -1, 0.7071067811865475, -1.0536712127723509e-08, -0.5, 0.5, -1, -1, 1, -1, 1, -1, -1, -1, -1, 1, -0.5, -0.5, 1, -0.7071067811865476, 1, -0.7071067811865476, 1, -1, 0.7071067811865475, -1.0536712127723509e-08, 1, -1, 0.7071067811865475, -1.0536712127723509e-08, -0.5, 0.5, -1, 1, -2, 1, -1, 1, 1, -1, 1, -0.5, -0.5, 1, 1, -0.7071067811865476, 1, 1, -0.7071067811865476, 1, 0.7071067811865475, -1.0536712127723509e-08, 1, 0.7071067811865475, -1.0536712127723509e-08, -0.5, 0.5]
rows = [2, 3, 12, 13, 20, 21, 22, 23, 30, 31, 34, 35, 0, 36, 1, 37, 2, 14, 36, 38, 3, 15, 37, 39, 34, 35, 6, 7, 12, 16, 24, 25, 26, 27, 30, 31, 40, 41, 4, 42, 5, 43, 6, 17, 42, 44, 7, 18, 43, 45, 40, 41, 8, 9, 10, 13, 19, 28, 29, 32, 33, 46, 47, 8, 14, 48, 9, 15, 49, 10, 48, 50, 11, 49, 51, 46, 47]
cols = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 13, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 16, 17, 17]
A = sparse.csc_matrix((data, (rows, cols)))
cones = [
  clarabel.ZeroConeT(20),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.ZeroConeT(1),
  clarabel.NonnegativeConeT(4),
  clarabel.SecondOrderConeT(6),
  clarabel.SecondOrderConeT(6),
  clarabel.SecondOrderConeT(6),
]

settings = clarabel.DefaultSettings()
solver = clarabel.DefaultSolver(P, q, A, b, cones, settings)
solver.solve()

Results in

-------------------------------------------------------------
           Clarabel.rs v0.7.1  -  Clever Acronym              

                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

problem:
  variables     = 18
  constraints   = 52
  nnz(P)        = 0
  nnz(A)        = 77
  cones (total) = 15
    :        Zero = 2,  numel = (20,1)
    : Nonnegative = 10,  numel = (1,1,1,1,...,4)
    : SecondOrder = 3,  numel = (6,6,6)

settings:
  linear algebra: direct / qdldl, precision: 64 bit
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,
  static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-5, max_scale = 1.0e5
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step      
---------------------------------------------------------------------------------------------
  0  +0.0000e+00  -5.7518e+00  5.75e+00  9.93e-01  6.89e+03  1.00e+00  1.33e+03   ------   
  1  +4.4363e-02  +7.9443e+01  7.94e+01  9.93e-01  3.89e-01  8.61e+01  2.96e+02  7.90e-01  
  2  +4.4016e-01  +4.4876e+02  4.48e+02  9.92e-01  4.14e-01  4.55e+02  2.93e+02  8.17e-02  
  3  +2.1750e+00  +1.5761e+03  7.24e+02  9.91e-01  1.30e+00  1.58e+03  2.79e+02  3.00e-01  
  4  +5.5621e+00  +2.2199e+03  3.98e+02  9.69e-01  1.81e-01  2.22e+03  4.48e+01  8.55e-01  
  5  +6.7417e+00  +1.0322e+02  1.43e+01  4.21e-01  4.90e-03  9.66e+01  1.22e+00  9.83e-01  
  6  +5.0782e+00  +6.1775e+00  2.16e-01  4.76e-02  5.32e-04  1.10e+00  1.31e-01  9.05e-01  
  7  +5.0782e+00  +6.1775e+00  2.16e-01  4.76e-02  5.32e-04  1.10e+00  1.31e-01  0.00e+00  
---------------------------------------------------------------------------------------------
Terminated with status = InsufficientProgress
solve time = 286.812µs

Admittedly there are some ugly small coefficients in there... this was generated through our higher-level GCS abstraction. But I thought it was best to capture the instance and report.

[goulart-paul]

In 0.7.1 we fixed a bug in the data equilibration that allowed the matrix scaling to exceed the limits set by the user. In the course of fixing that I also relaxed the bounds slightly (because hey, at least the bounds always work now....).

Probably that was a mistake. A temporary fix is to go back to the old limits:

settings.equilibrate_min_scaling = 1e-4;
settings.equilibrate_max_scaling = 1e+4;

We will make this reversion permanent in the next minor release assuming it doesn't break something elsewhere. I'll also add this case to our unit testing.

[goulart-paul]

Old limits have been restored in new v0.8.0.