Static regularization

[feodorp]

Hello,

I have two questions:

1. The first one is related to the implementation of the 'regularize_kkt()' function. The diagonal is saved and updated with the regularized diagonal at these lines:

piqp/include/piqp/sparse/kkt.hpp

Line 260 in a9e2cd9

kkt_diag(col) = PKPt.valuePtr()[PKPt.outerIndexPtr()[col + 1] - 1];

piqp/include/piqp/sparse/kkt.hpp

Line 267 in a9e2cd9

PKPt.valuePtr()[PKPt.outerIndexPtr()[ordering.inv(col) + 1] - 1] += rho_reg;

However, you are saving the permuted diagonal and then updating the original one. In the second case, ordering.inv(col) is used to get the correct position.

Another issue appears in 'solver.hpp'. The function 'm_kkt.regularize_and_factorize()' is called inside a 'while' loop:

piqp/include/piqp/solver.hpp

Line 400 in a9e2cd9

while (!m_kkt.regularize_and_factorize(m_enable_iterative_refinement))

But 'rho' and 'delta' are not updated in the KKT class. So, if 'regularize_and_factorize()' returns 'false' and iterative refinement is enabled, the loop will continue until 'max_factor_retries' is reached and will then exit with the 'PIQP_NUMERIC' status. Am I missing something?

2. The paper doesn’t provide an explanation, but why is static regularization needed if 'rho' and 'delta' are already inserted into the diagonal of the KKT matrix? Can they become too small?

[RSchwan]

Hi,

1. kkt_diag in the code stores the already permuted diagonal. There, I don't care about the ordering since I just save it temporary, since I use it to restore it later. In the second case on Line 267, I update the regularization, which is dependent on the ordering, since some are positive and some are negative.
   I think you are right about the issue in solver.hpp. I will fix this in the future. I never encountered a numeric's error in the first iteration, especially since delta and rho are already relatively big, but this should be fixed nevertheless.

2. In theory, regularization is not needed with rho and delta already inserted. In fact, it is disabled by default iterative_refinement_always_enabled = false. But as rho and delta go to zero to get faster convergence, sometimes numerical issues can arise. This is not in the paper, since it's not needed/enabled by default and because I didn't have enough space to add additional explanations in the paper. In the end, they are all heuristics to catch some very rare specific numerical issues which might happen. In general, these are all issues due to finite precision arithmetic of floating point numbers.