BLASFEO API explanation for Riccati recursion in HPIPM

[syangav]

Hello,

I was trying to understand the HPIPM code for two versions of Riccat recursion: square_root_alg = 0 for classical version (Algorithm 1 in paper, code starting from here) and square_root_alg = 1 for the square root version (Algorithm 3 in the same paper, code starting from here ).

For the square root version, I could understand the main pipeline per stage: trmm -> syrk -> potrf = chol, because the code matches clearly with the paper. However, I am still confused by the intermediate line of gead. Can you explain what it's for?

For the classical version, the code does not fully match with the paper, but I tried to figure out the pipeline per stage: gemm for line 3 of Algorithm 1 -> syrk for line 4 & 5, but with Q, R, S added -> potrf = chol for line 6 -> syrk for line 8. However, I am confused by the intermediate line of gead. I am also not so sure why the potrf in line returns the L in line 7 of Algorithm 3.

I think the primary reason of my confusion is that I don't fully get the BLASFEO API. For instance, in void blasfeo_dgead(int m, int n, double alpha, struct blasfeo_dmat *sA, int ai, int aj, struct blasfeo_dmat *sB, int bi, int bj), what are the roles of m, n, ai, aj, bi, bj? I have noticed that you post some related answers here, but unfortunately the BLAS API you pointed to is now 404 Not Found.

I also had the feeling that the algorithms on the paper is for easier reading while the code in HPIPM is of course for faster computation thus appears more compactly. Hence some information is lost. It would be very helpful if you may fill in the gap.

Thanks for your time in advance!

[giaf]

Hi,

regarding the line https://github.com/giaf/hpipm/blob/f7c4f502172b8ea5279c8f4d34afe882407526c3/ocp_qp/x_ocp_qp_kkt.c#L81 , the algorithms named FACT_SOLVE perform both the KKT factorization (by means of backward Riccati recursion) and the backward and forward solutions. For efficiency purposes, it is possible to write the algorithm such that the backward solution is merged with the backward factorization, by having specialized routines that operate on matrices of size n+1 x n, where the vector quantities of the backward solution are appended as the last row of the matrix. The specific line calling gead performs an operation specific to the backward solution.

About the classical Riccati recursion implementation, again I wrote specialized linear algebra routines, like a Cholesky factorization for non-squared (tall) matrices, that is equivalent to perform Cholesky of the top squared block and the solution with the Cholesky factor of the bottom block: this is what is happening here https://github.com/giaf/hpipm/blob/f7c4f502172b8ea5279c8f4d34afe882407526c3/ocp_qp/x_ocp_qp_kkt.c#L158 in the L computation at line 7 of the algorithm above. The advantage of the merged calls with the BLASFEO specialized routines is improved performance for small matrices, while with the standard BLAS and LAPACK routines this would require two distinct calls to dpotrf and dtrsm, that is less efficient.

The routine dgead just performs the sum (add) of two matrices of size m x n. The ai aj bi bj are just the row and column offsets of the submatrices that are added, like in all other BLASFEO routines.