OSQP vs simple SVD - Accuracy and Speed

[ipa-danb]

Hi All,

I'm having a bit of trouble, getting a good, reliable result from osqp, that is comparable to SVD.

Im using the python bindings to solve certain test cases with only equality constraints with OSQP and on the other hand with SVD from scipy linalg.

I create test cases with this function

def create_test_case(n: int = 100, ran=(0.0, 1.0)) -> (np.array, np.array):
    if isinstance(n, int):
        n1, n2 = n, n
    elif len(n) == 2:
        n1 = n[0]
        n2 = n[1]

    A = ran[1] * (np.random.random((n1, n2)) - ran[0])
    b = ran[1] * (np.random.random((n1, 1)) - ran[0])
    return A, b

then i use this function to solve this random problem with SVD

def solve_svd(A, b):

    U, s, Vh = svd(A, full_matrices=False)
    c = np.dot(U.T, b)
    w = np.dot(np.diag(1 / s), c)
    x = np.dot(Vh.conj().T, w)

    return x

Ill create a list of all permutations of a config file

vss = {
    "h": [20,25], # height of matrix
    "w": [15,20], # length of matrix
    "slack_no": [0.0,0.5,1.0], # how big is the added sparse matrix as a fraction of the original one
    "ran_max": [ 1.0, 100.0], # Scale of Problem to check for numerical issues
    "nu_scale":[1], # Scale on the sparse part of the extended matrix, 
    "p_scale":[1,100], # Scale of the optimization objective defined by P
    "eps_abs": [1e-6,1e-16,1e-20], # absolute termination value
    "eps_rel": [1e-6,1e-16,1e-20], # relative termination value
    "delta": [1e-6], # delta
    "rho": [0.1], # rho
    "cohort":[0,1], # how many iteration per parameter set
    "iterations": [1,10,100] # how many iterations are solved after setup (to check for influence of warmstart)
}

In the end, i test all those parameters and solve an extended problem with it

The extended problem is to test the influence of sparsity a

svd_err_list = []
osqp_err_list = []
conf_list = []
osqp_time_list = []
svd_time_list = []


for vs in track(list(dict_configs(vss)), description="Running...."):
    A, b = create_test_case((vs['h'],vs['w']),ran=(0,vs['ran_max']))
    A = A
    b = b
    slack_no = int(vs['h']*vs['slack_no'])
    A_ext = np.hstack((A, vs['nu_scale']*np.eye(vs['h'])[:,:slack_no]*A.max()))
    
    As = sparse.csc_matrix(A_ext)
    prob = osqp.OSQP()
    P = sparse.csc_matrix(vs['p_scale']*np.eye(As.shape[1]))
    q = np.zeros(As.shape[1])
    l = b
    u = b
    
    prob.setup(P, q, As, l, u, eps_abs=vs['eps_abs'] , eps_rel=vs['eps_rel'], verbose=False, max_iter=20000, scaling=10, delta=vs['delta'],rho=vs['rho'])
    osqp_t1_start = perf_counter() 
    for i in range(0,vs['iterations']):
        res = prob.solve()
    osqp_t1_time = (perf_counter() - osqp_t1_start)/vs['iterations']

    if res.info.status == "primal infeasible" or None in res.x:
        continue
        
    svd_t1_start = perf_counter() 
    x_svd = solve_svd(A_ext, b)
    svd_t1_time = perf_counter() - svd_t1_start

    osqp_err_list.append(np.linalg.norm(A_ext.dot(res.x.flatten()) - b.flatten()))
    svd_err_list.append(np.linalg.norm(A_ext.dot(x_svd) - b))
    conf_list.append(vs.copy())
    osqp_time_list.append(osqp_t1_time)
    svd_time_list.append(svd_t1_time)

    

I calculate the error of either SVD or OSQP solutions and the time to solve.

As you can see below, i can (almost only) reach parity with the SVD approach in terms of accuracy, but in terms of solve time, im definitly a magnitude off.

Advantage of SVD vs OSQP

(value is osqp time to solve - svd time to solve) there are values where osqp has an advantage, but only very minimal
(x and y axis are accuracy of the solution)

On the other hand solve time is a magnitude higher when using tighter tolerances, but ...

Accuracy is at most at the same level as the SVD solution.

I tinkered with it a bit, to see what the issue is, but honestly not what i expected.

My expectation would have been that OSQP has some advantages over SVD. In the sense of: high accuracy requirement --> slower solve time and low accuracy requirement --> faster solve time.

But from the comparison right now i would assume, that for this kind of problem, where only equality constraints are concerned, SVD is faster and more accurate than OSQP?
Or maybe i have some misconception and something is wrong in the setup?

Thanks in advance for your insight. If it would help, i can also provide the jupyter notebook for my setup.

[imciner2]

The matrices you are using are actually very tiny - it looks like they are only on the order of 20~30 in each dimension, so I wouldn't expect the SVD to really have a significant runtime with them. I think you might start to see the difference between OSQP and the SVD that you are expecting if you increase your matrix size.

Also, I would expect the SVD to usually be more accurate than OSQP when solving this problem - this is the nature of OSQP as a first-order solver.

[ipa-danb]

Also, I would expect the SVD to usually be more accurate than OSQP when solving this problem - this is the nature of OSQP as a first-order solver.

Yes, i agree. Its just that my assumption was that it should also come with a speed benefit, even if a tiny one. Especially for repeatedly solving the same (or similar) problem. I hope im doing it correctly, as in just repeatedly calling prob.solve() and OSQP automatically does warmstarting?

The matrices you are using are actually very tiny - it looks like they are only on the order of 20~30 in each dimension, so I wouldn't expect the SVD to really have a significant runtime with them. I think you might start to see the difference between OSQP and the SVD that you are expecting if you increase your matrix size.

Yes, the matrices are about 20-40 in each dimension. However, the same structured problem needs to be solved at a high frequency and sometimes also includes inequality constraints.
I'll try to increase the dimension of my testcases and hopefully see a speedup. Thank you!

[ipa-danb]

So i increased the dimensionality to 200 x 200 for the normal problem and 200 x 400 for the extended problem. I see still similar results.

Ill try increasing the dimensionality to 2000 x 2000 now

(i reduced the number of other parameters, so that it would run a bit quicker)

[goulart-paul]

If you have equality constraints only then OSQP is not needed. You can just solve the KKT problems for your equality constrained problem directly. This is actually what OSQP should end up doing anyway assuming you have the "polishing" option enabled, but it's going to waste a lot of time before getting there.

See, for example, §10.1.1 here. I think that is what you are describing -- apologies if I have misunderstood.

@imciner2 do you think it would be worth pre-checking whether a problem is entirely equality constrained before solving? In that case one could presumably just proceed directly to something very similar to the polishing step, bypassing the solver entirely. This could also be done in a slightly fancier way by checking that all bound are either 1) equalities, or 2) +/- inf bounds on both sides.

More complicated suggestion : in CVXOPT / ECOS / Clarabel, the initial iterate is chosen in a similar way, which means that problems that turn out to be purely equality constrained terminate in one step. A very similar initialization procedure is probably possible in OSQP, which might give better starting points and, hence, lower solve times overall. I think most of the code for doing that already lives in the solver.

[ipa-danb]

If you have equality constraints only then OSQP is not needed. You can just solve the KKT problems for your equality constrained problem directly. This is actually what OSQP should end up doing anyway assuming you have the "polishing" option enabled, but it's going to waste a lot of time before getting there.

Yes, that was my original solution. I now wanted to also be able to handle inequality constraints, if needed.

[imciner2]

@imciner2 do you think it would be worth pre-checking whether a problem is entirely equality constrained before solving? In that case one could presumably just proceed directly to something very similar to the polishing step, bypassing the solver entirely. This could also be done in a slightly fancier way by checking that all bound are either 1) equalities, or 2) +/- inf bounds on both sides.

More complicated suggestion : in CVXOPT / ECOS / Clarabel, the initial iterate is chosen in a similar way, which means that problems that turn out to be purely equality constrained terminate in one step. A very similar initialization procedure is probably possible in OSQP, which might give better starting points and, hence, lower solve times overall. I think most of the code for doing that already lives in the solver.

This is certainly work doing. I've opened #682 to track this so I don't forget about it.

[goulart-paul]

Another possible issue for you here is that your randomly generated A matrix is dense, but OSQP represents everything internally in sparse matrix format. SVD on a dense matrix is probably a lot faster than factorising some KKT version of the same data represented in sparse format.

Suggestion : compare your SVD run time to the time it takes to solve as an equality constrained QP (i.e. one shot KKT system solve as in the link above), but forcing the matrix to be sparse format. That might account for a lot of the difference.

Note also that, even in the best case, OSQP with polishing enabled will factor something twice. The first is the initial KKT factorisation, and the second will be from the polishing step. For an equality constrained QP this will turn out to be the same system, but I don't think the solver explicitly checks for that case.