osqp / osqpth

The differentiable OSQP solver layer for PyTorch

The Operator Splitting QP Solver

Visit our GitHub Discussions page for any questions related to the solver!

The documentation is available at osqp.org

The OSQP (Operator Splitting Quadratic Program) solver is a numerical optimization package for solving problems in the form

minimize        0.5 x' P x + q' x

subject to      l <= A x <= u

where x in R^n is the optimization variable. The objective function is defined by a positive semidefinite matrix P in S^n_+ and vector q in R^n. The linear constraints are defined by matrix A in R^{m x n} and vectors l and u so that l_i in R U {-inf} and u_i in R U {+inf} for all i in 1,...,m.

Citing OSQP

If you are using OSQP for your work, we encourage you to

• Cite the related papers,
• Put a star on this repository.

We are looking forward to hearing your success stories with OSQP! Please share them with us.

Bug reports and support

Please report any issues via the Github issue tracker. All types of issues are welcome including bug reports, documentation typos, feature requests and so on.

Numerical benchmarks

Numerical benchmarks against other solvers are available here.

osqpth's Issues

Make first release and upload to pypi

It would be nice to release the first (minor ?) version and upload it to pypi so that people can try it out. What do you think @bamos ?

Making the CSC shape/indices optional

The current interface requires that the user passes in the CSC shape/indices for A and P explicitly in the constructor and data passed into the forward pass is expected to be the values. Sometimes getting these values and reshaping the data from PyTorch can be cumbersome and it would be more convenient to also provide an interface that is closer to qpth's interface that takes A and P as dense inputs. We have a few options here:

1. Make the OSQP interface handle both cases. It may be slightly messy to have the code for the CSC and dense cases in the same class, but one potential advantage of this is that we could make the backward pass in the dense case slightly more efficient by computing dA and dP with a PyTorch batched dense outer product instead of indexing into all of the sparse elements. Also for another micro-optimization we could do here for the forward pass, we could convert the data to a numpy format and pass that directly into the csc_matrix constructor without needing to create a dense index map

2. Add an OSQP_Dense interface that infers the shape/indices and internally calls into OSQP (or a renamed version like OSQP_CSC). This would be slightly nicer to maintain and could be more reasonable so users understand what interface they're using. The backward pass may not be as efficient if we don't specifically override the dP and dA computations to be batched/dense

What are your thoughts/what do you want to move forward with here? 1) seems nicer so we could optimize for some of the dense operations but I could also understand 2) for simplicity for now.

---

I've currently hacked in 2) with:

class OSQP_Dense(Module):
    def __init__(self,
                 eps_rel=1e-05,
                 eps_abs=1e-05,
                 verbose=False,
                 max_iter=10000):
        super().__init__()
        self.eps_abs = eps_abs
        self.eps_rel = eps_rel
        self.verbose = verbose
        self.max_iter = max_iter

    def forward(self, P, q, A, l, u):
        # TODO: Better batch checks
        if A.ndimension() == 2:
            m, n = A.shape
        elif A.ndimension() == 3:
            n_batch, m, n = A.shape
        else:
            raise RuntimeError("A.ndimension() unexpected")

        P_idx = list(map(lambda I: np.array(I, dtype=np.int32),
                        zip(*product(range(n), range(n)))))
        A_idx = list(map(lambda I: np.array(I, dtype=np.int32),
                        zip(*product(range(m), range(n)))))

        x = OSQP(P_idx, P.shape, A_idx, A.shape)(P.view(-1), q, A.view(-1), l, u)
        return x

and am running some quick tests with:

def main():
    torch.manual_seed(0)
    n, m = 10, 5
    L = torch.randn(n, n, requires_grad=True)
    P = L.t().mm(L)
    q = torch.randn(n, requires_grad=True)
    A = torch.randn(m, n, requires_grad=True)
    l = -10*torch.ones(m, requires_grad=True)
    u = 10*torch.ones(m, requires_grad=True)

    P_idx = list(map(lambda I: np.array(I, dtype=np.int32),
                    zip(*product(range(n), range(n)))))
    A_idx = list(map(lambda I: np.array(I, dtype=np.int32),
                    zip(*product(range(m), range(n)))))

    # P, q, A, l, u = [x.cuda() for x in [P, q, A, l, u]]
    x = OSQP(P_idx, P.shape, A_idx, A.shape)(P.view(-1), q, A.view(-1), l, u)
    print(x)
    print(grad(x.sum(), L, retain_graph=True))

    x = OSQP_Dense()(P, q, A, l, u)
    print(x)
    print(grad(x.sum(), L))

Is this project still on-going?

Is this project still on-going or not in maintenance anymore?, What is your road map on future development? @bamos @sbarratt

Quick note on calling into OSQP again for the backward pass

Just leaving this here in case it's useful in the future (probably won't be)

I just tried replacing the sparse solve in the bw pass with another OSQP solve to see if re-using the factorizations helps and it seems ~twice as slow on some tests that I'm looking at. The relevant part of the bw pass is:

m = self.solvers[i]
l = np.zeros(self.m)
u = np.zeros(self.m)
l[ind_inactive] = -np.inf
u[ind_inactive] = np.inf
m.update(q=-dl_dx.squeeze(), l=l, u=u)
result = m.solve()
r_x = result.x
r_y = result.y
r_y[ind_inactive] = 0.
r_yl = r_y[ind_low]
r_yu = r_y[ind_upp]

And line-profiling this gives:

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   175        10    1017499.0 101749.9     67.9              m.update(q=-dl_dx.squeeze(), l=l, u=u)
   176        10     386027.0  38602.7     25.8              result = m.solve()

Debug active vs full KKT backward pass