Clarabel Insufficient Progress on small problems with PSD quadratic cost

[GGooeytoe]

Clarabel v0.9.0 reports InsufficientProgress when asked to find the distance between two (very nearly touching, long and thin) polyhedra, expressed as intersections of halfspaces, provided I use a quadratic cost to do so. If instead I use a SOC constraint to constrain a variable to be larger than the distance between the two points, it works fine. SCS v3.2.7 solves the quadratic cost version of the problem without apparent issues, so it seems like my problem formulation isn't totally ridiculous.

I've attached an image of the sets in question:

Code block to demonstrate the behavior:

import numpy as np
from scipy import sparse
import clarabel,scs

#sets are DEFINED by their halfspace representations, Ax<=b
#for convenience we provide their extreme points, obtained via the cdd library
#first set
alpha_A=np.array([[-4.11921381e-01, -6.64332521e-02, -9.08794476e-01],
       [-9.72744681e-01,  2.31878814e-01, -3.70742082e-05],
       [ 8.83203670e-01,  1.99760944e-01, -4.24319270e-01],
       [ 4.14416781e-01, -4.31691291e-02,  9.09062791e-01],
       [-7.81970762e-01,  1.50805745e-01,  6.04796953e-01],
       [-9.04451722e-01, -1.80402478e-01,  3.86551458e-01],
       [ 9.61311600e-01, -2.65706218e-01, -7.26650799e-02]])
alpha_b=np.array([-4220.90632955,  3567.14395948, -1128.27381502,  4158.41293366,
        5288.6319985 ,  4166.52825645,   201.20700731])

alpha_vertices=np.array([[  655.5086259 ,   595.64346631,  4303.85294329],
                        [-1325.02754866,  9825.81483979,  4526.82403285],
                        [-1379.0563613 ,  9599.16749269,  4567.88121847],
                        [-1310.7493929 ,  9885.72176826,  4584.7464179 ],
                        [ -801.24532049,  9327.83522528,  5382.61714336],
                        [-2294.25022068,   531.76863331,  5645.5338004 ]])

#second set
beta_A=np.array([[-9.60430356e-01, -2.78520254e-01, -2.68755228e-08],
       [ 8.83247067e-01,  1.99761726e-01, -4.24228559e-01],
       [ 4.11920623e-01,  6.64330805e-02,  9.08794832e-01],
       [-4.34250337e-01, -8.46788414e-01,  3.07206816e-01],
       [-3.33761819e-01,  2.40414084e-02,  9.42350815e-01],
       [ 4.14416485e-01, -4.31691300e-02,  9.09062926e-01],
       [ 6.81451653e-01, -1.11696004e-01, -7.23289463e-01],
       [-9.72758816e-01,  2.31819511e-01,  3.07910366e-06],
       [-7.95249048e-02, -1.76063549e-03, -9.96831325e-01],
       [ 3.47099628e-01, -9.51590092e-03, -9.37779983e-01],
       [ 9.61310917e-01, -2.65708536e-01, -7.26656384e-02]])
beta_b=np.array([ -745.03161722, -1127.86121616,  4220.90835374,   655.97763879,
        4995.66937437,  4158.40577418, -2562.9235196 ,  3566.85188352,
       -3028.80886202, -2874.55064316,   201.20818082])

beta_vertices=np.array([[  603.29804383,   594.59139481,  4327.59603455],
       [  655.48605197,   595.72175725,  4303.85865231],
       [-1325.04585702,  9826.12677399,  4526.81073321],
       [-1662.8002849 ,  8408.86307811,  3156.23905064],
       [-1662.80182282,  8408.85662467,  3156.23918473],
       [  648.87050235,   437.44242992,  4215.55409187],
       [-1662.78983853,  8408.81529884,  3156.23830164],
       [  650.12959194,   433.1006676 ,  4248.08693587],
       [  671.63495527,   498.92327686,  4291.90009849],
       [  657.48927883,   446.66058654,  4295.8668939 ],
       [  643.37208262,   456.40279923,  4302.76516409],
       [-1662.79950876,  8408.84851554,  4497.8258165 ],
       [-1379.14302434,  9599.12421171,  4567.92476365],
       [-1218.74484163,  6877.60049365,  4694.16655204]])
print(
"""
first, attempt to find the distance between the two sets by constraining points to lie in each set and penalizing the squared distance:

min (1/2)*[x.T, y.T][2*I,-2*I][x]  s.t.
x,y                 [-2*I,2*I][y]

        x in alpha
        y in beta
"""
)
dim=3
zero_cost_vector=np.zeros(2*dim)
sq_distance_cost=sparse.diags([2*np.ones(2*dim),-2*np.ones(dim),-2*np.ones(dim)],[0,dim,-dim],format="csc")
first_constraint_matrix=sparse.block_array([[alpha_A,None],[None,beta_A]],format="csc")
first_constraint_rhs=np.concatenate([alpha_b,beta_b])
first_cones=[clarabel.NonnegativeConeT(len(first_constraint_rhs))]

settings=clarabel.DefaultSettings()
settings.verbose=True
solver1=clarabel.DefaultSolver(sq_distance_cost,zero_cost_vector,first_constraint_matrix,first_constraint_rhs,first_cones,settings)
solution1=solver1.solve()

first_cones_scs={"l":len(first_constraint_rhs)}
scs_solver1=scs.SCS({"P":sq_distance_cost,"A":first_constraint_matrix,"b":first_constraint_rhs,"c":zero_cost_vector},first_cones_scs,verbose=True)
scs_solution1=scs_solver1.solve()
print(
"""
second, attempt to find the distance by constraining points to lie in each set and doing an epigraph trick to penalize the actual distance:

min     f                           s.t.
x,y,f

        x in alpha
        y in beta
        [f,x.T-y.T] in SOC(dim+1)
"""
)
f_cost_vector=np.zeros(2*dim+1)
f_cost_vector[-1]=1
zero_quadratic_cost=sparse.csc_array((2*dim+1,2*dim+1))
second_constraint_matrix=sparse.block_array([[alpha_A,None,None],
                                             [None,beta_A,None],
                                             [None,None,[[-1]]],
                                             [sparse.eye(dim),-sparse.eye(dim),None]],format="csc")
second_constraint_rhs=np.concatenate([alpha_b,beta_b,np.zeros(dim+1)])
second_cones=[clarabel.NonnegativeConeT(len(first_constraint_rhs)),clarabel.SecondOrderConeT(dim+1)]
solver2=clarabel.DefaultSolver(zero_quadratic_cost,f_cost_vector,second_constraint_matrix,second_constraint_rhs,second_cones,settings)
solution2=solver2.solve()
print(
"""
third, show that Clarabel is perfectly happy to find a point that is IN both sets
min     0                           s.t.
x,y,f

        x in alpha
        y in beta
"""
)
solver3=clarabel.DefaultSolver(sparse.csc_array((2*dim,2*dim)),zero_cost_vector,first_constraint_matrix,first_constraint_rhs,first_cones,settings)
solution3=solver3.solve()
print(
"""
fourth, show that Clarabel (but not SCS) also fails if asked to minimize distance between two points both in alpha:

min (1/2)*[x.T, y.T][2*I,-2*I][x]  s.t.
x,y                 [-2*I,2*I][y]

        x,y in alpha
"""
)
fourth_constraint_matrix=sparse.block_array([[alpha_A,None],[None,alpha_A]],format="csc")
fourth_constraint_rhs=np.concatenate([alpha_b,alpha_b])
fourth_cones=[clarabel.NonnegativeConeT(len(fourth_constraint_rhs))]

solver4=clarabel.DefaultSolver(sq_distance_cost,zero_cost_vector,fourth_constraint_matrix,fourth_constraint_rhs,fourth_cones,settings)
solution4=solver4.solve()

fourth_cones_scs={"l":len(fourth_constraint_rhs)}
scs_solver4=scs.SCS({"P":sq_distance_cost,"A":fourth_constraint_matrix,"b":fourth_constraint_rhs,"c":zero_cost_vector},fourth_cones_scs,verbose=True)
scs_solution4=scs_solver4.solve()

print(
"""
fifth, show that Clarabel (but not SCS) also fails if asked to minimize distance between two points both in beta:

min (1/2)*[x.T, y.T][2*I,-2*I][x]  s.t.
x,y                 [-2*I,2*I][y]

        x,y in beta
"""
)
fifth_constraint_matrix=sparse.block_array([[beta_A,None],[None,beta_A]],format="csc")
fifth_constraint_rhs=np.concatenate([beta_b,beta_b])
fifth_cones=[clarabel.NonnegativeConeT(len(fifth_constraint_rhs))]

solver5=clarabel.DefaultSolver(sq_distance_cost,zero_cost_vector,fifth_constraint_matrix,fifth_constraint_rhs,fifth_cones,settings)
solution5=solver5.solve()

fifth_cones_scs={"l":len(fifth_constraint_rhs)}
scs_solver5=scs.SCS({"P":sq_distance_cost,"A":fifth_constraint_matrix,"b":fifth_constraint_rhs,"c":zero_cost_vector},fifth_cones_scs,verbose=True)
scs_solution5=scs_solver5.solve()

[goulart-paul]

I think the issue here is related to the formatting of python sparse matrices, which are not guaranteed to be in "canonical" format (i.e. no duplicates, data within each column appearing in row order).

If I add this

sq_distance_cost.sort_indices()
sq_distance_cost.sum_duplicates()

then all the problems solve. Note that I'm not sure that the solutions are actually correct, because I didn't check whether the A matrices are also in non-standard format.

I will have a look and see where the most appropriate place is for a fix, since it's probably asking for trouble to expect users of the python interface to realise that they need to manually force a reshuffle of the matrix internals.

When using the rust interface directly, we have a public check_format utility that can be run on the CscMatrix type before passing to the solver, but we don't use it internally. We should at least provide user facing sorting and deduplication tools as above in rust if that isn't available already.

[goulart-paul]

I have fixed this issue permanently (I hope!) in #140, but a temporary fix is to manually sort indices and deduplicate as suggested above.

[GGooeytoe]

I pulled #140 and built and now Clarabel solves all the problems in the example script I posted. I confirm that this fixed the issue!