Linear solvers

We suppose that the KKT system has been assembled previously into a given AbstractKKTSystem. Then, it remains to compute the Newton step by solving the KKT system for a given right-hand-side (given as a AbstractKKTVector). That's exactly the role of the linear solver.

If we do not assume any structure, the KKT system writes in generic form

\[K x = b\]

with $K$ the KKT matrix and $b$ the current right-hand-side. MadNLP provides a suite of specialized linear solvers to solve the linear system.

Inertia detection

If the matrix $K$ has negative eigenvalues, we have no guarantee that the solution of the KKT system is a descent direction with regards to the original nonlinear problem. That's the reason why most of the linear solvers compute the inertia of the linear system when factorizing the matrix $K$. The inertia counts the number of positive, negative and zero eigenvalues in the matrix. If the inertia does not meet a given criteria, then the matrix $K$ is regularized by adding a multiple of the identity to it: $K_r = K + \alpha I$.

Note

We recall that the inertia of a matrix $K$ is given as a triplet $(n,m,p)$, with $n$ the number of positive eigenvalues, $m$ the number of negative eigenvalues and $p$ the number of zero eigenvalues.

Factorization algorithm

In nonlinear programming, it is common to employ a LBL factorization to decompose the symmetric indefinite matrix $K$, as this algorithm returns the inertia of the matrix directly as a result of the factorization.

Note

When MadNLP runs in inertia-free mode, the algorithm does not require to compute the inertia when factorizing the matrix $K$. In that case, MadNLP can use a classical LU or QR factorization to solve the linear system $Kx = b$.

Solving a KKT system with MadNLP

We suppose available a AbstractKKTSystem kkt, properly assembled following the procedure presented previously. We can query the assembled matrix $K$ as

K = MadNLP.get_kkt(kkt)

6×6 SparseArrays.SparseMatrixCSC{Float64, Int32} with 13 stored entries:
 2.0     ⋅     ⋅     ⋅    ⋅    ⋅ 
 0.0  200.0    ⋅     ⋅    ⋅    ⋅ 
  ⋅      ⋅    0.0    ⋅    ⋅    ⋅ 
  ⋅      ⋅     ⋅    0.0   ⋅    ⋅ 
 0.0    0.0  -1.0    ⋅   0.0   ⋅ 
 1.0    0.0    ⋅   -1.0   ⋅   0.0

Then, if we want to pass the KKT matrix K to Lapack, this translates to

linear_solver = LapackCPUSolver(K)

LapackCPUSolver{Float64, SparseArrays.SparseMatrixCSC{Float64, Int32}}(sparse(Int32[1, 2, 5, 6, 2, 5, 6, 3, 5, 4, 6, 5, 6], Int32[1, 1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 6], [2.0, 0.0, 0.0, 1.0, 200.0, 0.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0], 6, 6), [6.32624201605996e-310 6.32624438756083e-310 … 6.32624453767853e-310 6.32624517025647e-310; 6.32624637047315e-310 6.3262443875632e-310 … 6.3262445376801e-310 6.3262451702802e-310; … ; 6.32589128938435e-310 6.32624637047315e-310 … 6.3262458591057e-310 6.3262451703039e-310; 6.32624646270334e-310 6.32624438756795e-310 … 6.32624453768327e-310 6.3262451703197e-310], 6, Float64[], Float64[], Float64[], [384.0, 1.22280453e-315, 4.243991582e-314, 8.4879831644e-314, NaN, 1.222804687e-315, 0.0, 2.121995791e-314, 1.32434862e-315, 1.222804845e-315  …  0.0, 1.22281931e-315, 4.2439915824e-314, 6.32625150696557e-310, NaN, 1.22281947e-315, 0.0, NaN, 1.324347593e-315, 1.222819628e-315], 384, [128044627922608], -1, Base.RefValue{Int64}(0), [128044616693160, 4, 563204, 4, 4, 128044627582896], MadNLP.LapackOptions(MadNLP.BUNCHKAUFMAN), MadNLP.MadNLPLogger(MadNLP.INFO, MadNLP.INFO, nothing))

The instance linear_solver does not copy the matrix $K$ and instead keep a reference to it.

linear_solver.A === K

true

That way every time we re-assemble the matrix $K$ in kkt, the values are directly updated inside linear_solver.

To compute the factorization inside linear_solver, one simply as to call:

MadNLP.factorize!(linear_solver)

LapackCPUSolver{Float64, SparseArrays.SparseMatrixCSC{Float64, Int32}}(sparse(Int32[1, 2, 5, 6, 2, 5, 6, 3, 5, 4, 6, 5, 6], Int32[1, 1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 6], [2.0, 0.0, 0.0, 1.0, 200.0, 0.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0], 6, 6), [2.0 0.0 … 0.0 0.0; 0.0 200.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.5 0.0 … -1.0 -0.5], 6, Float64[], Float64[], Float64[], [384.0, 1.22280453e-315, 4.243991582e-314, 8.4879831644e-314, NaN, 1.222804687e-315, 0.0, 2.121995791e-314, 1.32434862e-315, 1.222804845e-315  …  0.0, 1.22281931e-315, 4.2439915824e-314, 6.32625150696557e-310, NaN, 1.22281947e-315, 0.0, NaN, 1.324347593e-315, 1.222819628e-315], 384, [128044627922608], -1, Base.RefValue{Int64}(0), [1, 2, -5, -5, -6, -6], MadNLP.LapackOptions(MadNLP.BUNCHKAUFMAN), MadNLP.MadNLPLogger(MadNLP.INFO, MadNLP.INFO, nothing))

Once the factorization computed, computing the backsolve for a right-hand-side b amounts to

nk = size(kkt, 1)
b = rand(nk)
MadNLP.solve!(linear_solver, b)

6-element Vector{Float64}:
  0.8318843119784843
  0.004541507586279584
 -0.830130488834118
  0.40606281675462175
 -0.7767075774951446
 -0.6690029233317701

The values of b being modified inplace to store the solution $x$ of the linear system $Kx =b$.