search

giaf / hpipm

High-performance interior-point-method QP and QCQP solvers
Other
553 stars 130 forks source link

all matrices in C interface are in form of column major? #128

Closed rushahead closed 1 year ago

rushahead commented 1 year ago

Hello. I am using HPIPM now and just wanna to make sure if all matrices must be in column major? and matrix H seems not.

For example. I solve a problem with OSQP.

import osqp as op import numpy as np from scipy import sparse P = sparse.csc_matrix([[0,0,0,0,0,0], [0,0,0,0,0,0], [0,0,20,0,0,0], [0,0,0,20,0,0], [0,0,0,0,40,0], [0,0,0,0,-40,40]]) q = np.array([0,0,-200,-200,0,0]) A = sparse.csc_matrix([[1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0], [0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 1], [0, 0, 0, 0, -1, 1], [0, 0,-1, 1, -0.5, -0.5], [-1,1,-1, 0, -0.333333,-0.1666667], [1, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 1, 0]]) l = np.array([-1, -1, -1, -1,-5,-5,-6,0,0,0,1,0.5]) u = np.array([200,200,45,45,5, 5, 6, 0,0,0,1,0.5]) prob = op.OSQP() prob.setup(P, q, A, l, u, alpha=1.0) res = prob.solve() print(res.x)'

results are : [2.36370683e-08 1.49054196e+00 9.99999979e-01 2.22162614e+00 4.99999850e-01 1.94325252e+00] but if I use HPIPM to solve it can't get same resolution.

int nv = 6; int ne = 0; int nb = 0; int ng = 12; int ns = 0; int nsb = 0; int nsg = 0; double H[] = {0,0,0,0,0,0, 0,0,0,0,0,0, 0,0,20,0,0,0, 0,0,0,20,0,0, 0,0,0,0,40, -40, 0,0,0,0, 0,40}; double g[] = {0,0,-200,-200,0,0}; double C[] = {1,0,0,0,0,0,0, 0, -1, 1,0,0, 0,1,0,0,0,0,0, 0, 1, 0,0,0, 0,0,1,0,0,0,0,-1, -1, 0,1,0, 0,0,0,1,0,0,0, 1, 0, 0,0,0, 0,0,0,0,1,0,-1,-0.5,-0.333333,0,0,1, 0,0,0,0,0,1,1, -0.5,-0.166667,0,0,0}; double d_lg[] = {-1,-1,-1,-1,-5,-5,-6,0,0,0,1,0.5}; double d_ug[] = {200,200,45,45,5,5,6,0,0,0, 1,0.5};

results are : [0.00000 1.56482 1.00000 2.44444 0.50000 2.38889] By the way, I don't use H with column-major, then I could get same result.

double H[] = {0,0,0,0,0,0, 0,0,0,0,0,0, 0,0,20,0,0,0, 0,0,0,20,0,0, 0,0,0,0,40, 0, 0,0,0,0, -40,40};

It's confusing to me. Could anyone help me with this ? Thanks in advance.

giaf commented 1 year ago

Yes all matrices in the C interface are expected to be in column-major, unless another format is explicitly used in the routine naming, as e.g. here https://github.com/giaf/hpipm/blob/master/ocp_qp/x_ocp_qp.c#L610 where row-major is expected. Furthermore, in case of symmetric matrices (as e.g. the Hessian matrix), only the lower triangular part of the symmetric matrix is accessed, and the upper triangular part is disregarded.

So in your example above, in the HPIPM case, when you store the H matrix in column-major, the Hessian used internally looks like

double Hc[] = {0,0,0,0,0,0,
0,0,0,0,0,0,
0,0,20,0,0,0,
0,0,0,20,0,0,
0,0,0,0,40, -40,
0,0,0,0, -40,40};

while when you store the matrix in row-major (i.e. you pass the transposed of your matrix H above), the Hessian used internally looks like

double Hr[] = {0,0,0,0,0,0,
0,0,0,0,0,0,
0,0,20,0,0,0,
0,0,0,20,0,0,
0,0,0,0,40, 0,
0,0,0,0, 0,40};

i.e. without the off-diagonal terms. You can double check that the solutions returned by HPIPM and that you report above are indeed correct and consistent with this Hessian matrices by using another tool as e.g. qp in octave.

I am not familiar with the OSQP interface but there may be similar assumptions on e.g. only the upper or lower triangular part of the symmetric Hessian matrix H being accessed.