QP: cvxopt vs lava

[elrond91]

I compared the performance of the QP solver provided by Lava against the optimization solver cvxopt (https://cvxopt.org/), and I got very different values; I have tried modifying iterations numbers with alpha and beta values, but still is the same behavior. Am I missing something?

Packages:
lava-nc 0.2.0
cvxopt 1.2.7

Code:

import cvxopt
import numpy as np
from lava.lib.optimization.problems.problems import QP
from lava.lib.optimization.solvers.qp.solver import QPSolver

if __name__ == "__main__":
    """
    Lava:
    """
    Q = np.array([[100, 0, 0], [0, 15, 0], [0, 0, 5]])
    p = np.array([[1, 2, 1]]).T
    A = -np.array([[1, 2, 2], [2, 100, 3]])
    k = -np.array([[-50, 50]]).T

    alpha, beta = 0.001, 1
    alpha_d, beta_g = 10000, 10000
    iterations = 400
    problem = QP(Q, p, A, k)
    solver = QPSolver(
        alpha=alpha,
        beta=beta,
        alpha_decay_schedule=alpha_d,
        beta_growth_schedule=beta_g,
    )
    print('Loihi: ', solver.solve(problem, iterations=iterations))

    """
    cvxopt:
    """
    Q = np.array([[100, 0, 0], [0, 15, 0], [0, 0, 5]]).astype('float')
    p = np.array([1, 2, 1]).astype('float')
    A = -np.array([[1, 2, 2], [2, 100, 3]]).astype('float')
    k = -np.array([-50, 50]).astype('float')

    sol = cvxopt.solvers.qp(cvxopt.matrix(Q), cvxopt.matrix(p), cvxopt.matrix(A), cvxopt.matrix(k))
    print('cvxopt: ', sol['x'])

Output:

Loihi:  [[-1.00000000e-02]
 [ 8.51815265e+00]
 [ 1.44819123e+01]]

cvxopt:  [-8.09e-03]
[ 5.04e-01]
[-1.43e-01]

[ashishrao7]

Hi @elrond91,

Thanks for initiating a discussion. You do not seem to be preconditioning the problem before using it in the Loihi Solver. That could be one of the reasons why the solutions don't match. The QP solver for Lava is a solver that uses first-order dynamics and hence ill-conditioned systems will not work well without some form of preconditioning. cvxopt would internally take care of these issues.

[elrond91]

Thanks for the response, I reduce the search space and also added the preconditioner:

"""
Lava:
"""
Q = 2 * np.array([[2, .5], [.5, 1]])
p = np.array([[1.0, 1.0]]).T
A = np.array([[-1.0, 0.0], [0.0, -1.0]])
k = np.array([[-1.0, -1.0]]).T

pre_mat_Q = np.diag(1 / np.linalg.norm(Q, axis=1))
Q_pre = pre_mat_Q @ Q @ pre_mat_Q
p_pre = pre_mat_Q @ p

print(np.linalg.cond(Q), np.linalg.cond(Q_pre))

pre_mat_A = np.diag(1 / np.linalg.norm(A, axis=1))
A_pre = pre_mat_A @ A @ pre_mat_Q
k_pre = pre_mat_A @ k

alpha, beta = 0.1, 1
alpha_d, beta_g = 10000, 10000
iterations = 5000
problem = QP(Q_pre, p_pre, pre_mat_A, k_pre)
solver = QPSolver(
    alpha=alpha,
    beta=beta,
    alpha_decay_schedule=alpha_d,
    beta_growth_schedule=beta_g,
)
sol = solver.solve(problem, iterations=iterations)
print('Lava: ', pre_mat_Q @ sol)
"""
cvxopt:
"""
Q = 2 * np.array([[2, .5], [.5, 1]]).astype('float')
p = np.array([1, 1]).astype('float')
A = np.array([[-1.0, 0.0], [0.0, -1.0]]).astype('float')
k = np.array([-1.0, -1.0]).astype('float')

sol = cvxopt.solvers.qp(cvxopt.matrix(Q), cvxopt.matrix(p), cvxopt.matrix(A), cvxopt.matrix(k))
print('cvxopt: ', sol['x'])

Output
`Lava: [[-0.2401519 ]
[-0.44462554]]

cvxopt: [ 1.00e+00]
[ 1.00e+00]

`

Any other suggestion that I could try?

Thanks!

[ashishrao7]

You seem to be not decaying your alpha/ increasing your beta at all (This is because alpha_d = 10000 and beta_g = 10000 while iterations = 5000). This means that alpha_d/beta_g is halved/doubled every 10000 timesteps but your simulation runs only 5000 steps. Thus it could be getting stuck at a sub-optimal solution. I now notice that you also did the same thing in your earlier example. Perhaps try reducing that to a smaller number (say 250 if you are using 5000 iterations).

Also preconditioners can be different based on the type of problems you intend to solve. For example, we envision solving sparse block diagonal QP when it comes to robotic problems and therefore we use a preconditioner that works well with these kinds of problems. Please check if the condition number of your problem reduces after you precondition it with your preconditioner. If your problem doesn't have a very large condition number to begin with, preconditioning may not be necessary.

[LVAjay]

Hi,
I tried using QPSolver in another program having lower alpha_d and beta_g values as suggested above. But I am getting an Error as explained below:
While running the QPSolver with alpha =10, and alpha_d value less than solver iterations: I hit an exception in the file runtime.py. Further analyzing the exception, printing the traceback for the exception in the file runtime.py- ran into the below statement.

TypeError: ufunc 'right_shift' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''.

I am afraid that the exception is due to right shift operation performed on alpha which is float64 dtype.

[ashishrao7]

Hi @LVAjay,

Thank you for pointing that out. This will be fixed soon in a future release sometime this month. Could you please file an issue for this since that would be the appropriate place for this? In the meanwhile, you can work around this problem by replacing the right-shift and left-shift of the constants by halving (alpha=alpha/2) or doubling (beta=beta*2) respectively.

[LVAjay]

Hi @ashishrao7
Thanks for your response, I will file this issue. Also, I have encountered another issue after updating work around for alpha and beta. Which I am afraid to see an attribute error for self.decay_counter.shape. The detailed error is as shown below:

self.decay_counter = np.zeros(self.decay_counter.shape)
AttributeError: 'int' object has no attribute 'shape'.

May I know commenting this code part for execution, will have an effect on results? IF Yes, What can be done to have an work around for this. Thanks in advance.

[ashishrao7]

Hi @LVAjay,

Thanks for rightly pointing out these issues. These issues have been fixed in my version of the fork. However, they've not been incorporated into the main lava-optimization repo yet. This will happen soon this month. In the meanwhile, try replacing it with something like this,

self.alpha = self.alpha/2
self.decay_counter = 0

Let's continue having this discussion in an issue thread that you file :)

[LVAjay]

Hi @ashishrao7
Yeah sure:). Thanks for your quick response.