ecos crash on a very sparse problem

[uricohen]

The code below reproduces a non-deterministic bug which causes ecos to crash with the following message, on both Windows (python 3.10) and Gentoo Linux (python 3.8). Despite using the same inputs (by setting numpy's seed), the crash occurs after a varying number of trials. The code uses the correct parameters to trigger a crash 100% of the time (on my machine).

RuntimeError: Internal problem occurred in ECOS while setting up the problem.
Please send a bug report with data to Alexander Domahidi.
Email: domahidi@control.ee.ethz.ch

The code:

import numpy as np
import time
import math
import scipy
import qpsolvers

np.random.seed(1)

def relu_inv(x): d=np.empty_like(x); d[:]=x[:]; d[x[:]<=0]=np.nan; return d
def relu(x): return (x>0)*x
def invH(x): return math.sqrt(2)*scipy.special.erfcinv(2*x)

solver = "ecos"
def optimize_fixed_points(X, H, threshold=0):
    (N, P) = X.shape
    assert np.all(H.shape == (N, P))
    
    # Solve quadratic problem for each row
    W = np.empty((N, N)); W[:] = np.nan
    Q = np.eye(N-1)
    found = 0
    for ii in range(N):
        I = np.concatenate((np.array(range(ii),dtype=int),np.array(range(ii+1,N),dtype=int)))
        I = I.reshape(I.shape[0],1)
        J = np.nonzero(np.isfinite(H[ii,:])*1)[0]
        coJ = np.nonzero(np.isnan(H[ii,:])*1)[0]
        assert len(J)+len(coJ) == P
        c = np.zeros(I.shape[0])
        lb = None
        A = X[I,J].T
        b = H[ii,J]
        G = X[I,coJ].T
        z = np.ones(G.shape[0])*threshold
        w = qpsolvers.solve_qp(Q, c, G, z, A, b, lb, solver=solver, abstol=1e-8, reltol=1e-6)
        if w is not None:
            assert np.all(np.isfinite(w))
            assert np.linalg.norm(A@w-b)<1e-6
            assert np.all(G@w<z+1e-6)
            W[ii,ii] = 0
            W[ii,I.T] = w
            found += 1
    return (W, found)

# Problem parameters
N = 128; f = 1/64; K = np.round(f*N); P = 1
print("N=%d K=%d P=%d" % (N,K,P))

for i in range(100):
    tictocT = time.time()
    # Create patterns
    X = relu(np.random.randn(N, P) - invH(f))
    H = relu_inv(X)
    (W, nfound) = optimize_fixed_points(X, H)
    d = X - relu(W@X)
    assert np.linalg.norm(d[np.isfinite(d[:])]) < 1e-8, "Fixed-points violation: %1.1e"%np.linalg.norm(d[np.isfinite(d[:])])
    assert nfound<N or np.all(np.isfinite(W))
    print("Done #%d (found %d/%d, took %1.1fsec)."%(i, nfound, N, time.time()-tictocT))

[uricohen]

The problem seems to occur mainly when A=[1, N], ie when there is a single equality constraint.
However, I could not use this knowledge to find a working work-around. I tried changing A to be [2, N] (by duplication or by adding zero constraint) or absorbing it into G (i.e., adding two appropriate inequality constraints), but the crash still happens.