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$.
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.
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], [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.90899195008305e-310 6.9089919496791e-310 … 5.7293886388e-313 2.652494739364e-312; 6.90899195008305e-310 6.9089919496791e-310 … 5.09278989955e-313 2.88591427633e-312; … ; 6.9089919496791e-310 1.27319747489e-312 … 2.33419537061e-312 2.33419537077e-312; 6.9089919496791e-310 1.506617011976e-312 … 2.334195370605e-312 2.334195370605e-312], [2.0e-323], -1, Base.RefValue{Int64}(0), Dict{Symbol, Any}(), 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], [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], [384.0, 1.060997899e-314, 6.4e-323, 3.81288324e-316, 5.4e-323, 2.9834573e-316, 2.9834577e-316, 8.487983164e-314, 2.98358576e-316, 2.4e-322 … 0.0, 2.673714696616e-312, 2.98360394e-316, 4.2439915824e-314, 2.98360473e-316, 1.69759663287e-313, 5.4e-323, 3.70967886e-316, 2.1612977e-316, 2.9836075e-316], 384, Base.RefValue{Int64}(0), Dict{Symbol, Any}(:ipiv => [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.42071023896062926
0.0015444497363156934
-0.6499530324454892
0.27851501957973374
-0.20598500432852984
-0.11994588588809041
The values of b
being modified inplace to store the solution $x$ of the linear system $Kx =b$.