Infeasibility for basic QP

[jamespinkerton]

Hi, I tried using the library and am getting infeasibility, so I got things down to a small working example. I can't understand why this solution is yielding infeasibility. Any help would be greatly appreciated. Here's the code I have:

import numpy as np
import osqp  # type: ignore
import piqp  # type: ignore
from scipy import sparse  # type: ignore

P = np.array(
    [
        [1.0, 0.0, 0.0, 0.0],
        [0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 1.0, 0.0],
        [0.0, 0.0, 0.0, 1.0],
    ]
)
c = np.array([1.0, 1.0, 1.0, 1.0])

A = np.array([[0.0, 0.0, 0.0, 0.0]])
b = np.array([0.0])

G = np.array(
    [
        [1.0, 0.0, 0.0, 0.0],
        [1.0, 0.0, -1.0, 0.0],
        [-1.0, 0.0, -1.0, 0.0],
        [-1.0, 0.0, 0.0, 0.0],
        [-1.0, 0.0, 1.0, 0.0],
        [1.0, 0.0, 1.0, 0.0],
    ]
)
h = np.array([1.0, 1.0, 1.0, 1.0, np.inf, np.inf])

problem = piqp.DenseSolver()
problem.settings.verbose = True
problem.settings.compute_timings = True
problem.setup(
    P=P,
    c=c,
    A=A,
    b=b,
    G=G,
    h=h,
)
problem.solve()
x = problem.result.x

easy_answer = np.array([0.0, 0.0, 0.0, 0.0])
assert A @ easy_answer == b
assert (G @ easy_answer <= h).all()


def loss(answer: np.ndarray) -> float:
    return 0.5 * answer.T @ P @ answer + c.T @ answer


print("easy loss", loss(easy_answer))
print("solved loss", loss(x))

print("x", x)

problem = osqp.OSQP()
lower = np.array([-np.inf, -np.inf, -np.inf, -np.inf, -np.inf, -np.inf])
problem.setup(sparse.csc_matrix(P), c, sparse.csc_matrix(G), lower, h)
r = problem.solve()
print("x", r.x)

[jamespinkerton]

If I replace the np.inf in h with something reasonable then I get the same solutions in osqp and piqp. This makes me think there's a bug on the piqp end.

[RSchwan]

Unfortunately, this is a limitation of interior-point methods in general. The inequality constraints are reformulated using slack variables, i.e., the constraints become $Gx − h + s = 0$. If $h$ contains infinity (or very big values), this leads to numerical issues during initialization and solving. On the contrary, OSQP projects into the inequality constraints, which is just a max/min operator. Hence, the infinity values are no issue numerically.

To deal with this issue, normally one would pre-filter the problems, i.e., remove the infinity/inactive inequality constraints. This is currently not happening internally in PIQP for general inequality constraints, only for box constraints. The reason for this is that I designed PIQP to be used for solving the same kind of problem over and over (control applications). PIQP allocates memory only ones in the setup step. If the active inequalities changed (e.g. one bound is changed from inf to something finite), the sparsity pattern would change, which would result in allocating new memory (which I try to avoid).

I have been thinking about this for a while, since it causes problems for people who don't know this (like you) and hence is not really user-friendly. I could add the preprocessing step internally, which can be disabled for people who really know what they are doing to get maximum performance and no subsequent allocations. This should fix the problem for the average user. What do you think?

In the meantime, try to avoid using inf constraints in the general inequalities...

Also, instead of passing a $A$ matrix filled with zeros, you can pass a zero matrix, removing one constraint:

A = np.zeros((0, 4))
b = np.zeros(0)

[jamespinkerton]

OK that makes sense! In that case I have a follow-up question. When I “update” the problem for solving something over and over again, which vectors can I update without allocating new memory? There’s c, h, b, x_lb and x_ub. You’re making me thing I can’t update h without new memory allocation?

[RSchwan]

You can update everything, but you should not change the sparsity pattern, i.e., only the values of sparse matrices are updated but not the structure. This means you can also update P, A and G as long as the sparsity pattern doesn't change. For example, if you want to update A, h, and x_ub you can do this using problem.update(P=None, c=None, A=A_new, b=None, G=None, h=h_new, x_lb=None, x_ub=x_ub_new). In the current version of PIQP there are no internal memory allocations for this operation, i.e., the memory is reused (also for the factorization).

But if I have to filter out inequality constraints based on the values in h, I have to remove/add rows from G internally which would change the sparsity pattern of the KKT matrix, requiring allocating new memory for the KKT matrix and its factorization.