KKT systems are linear system with a special KKT structure. MadNLP uses a special structure AbstractKKTSystem to represent internally the KKT system. The AbstractKKTSystem fulfills two goals:
- Store the values of the Hessian of the Lagrangian and of the Jacobian.
- Assemble the corresponding KKT matrix $K$.
We recall that at each iteration the interior-point algorithm aims at solving the following relaxed KKT equations ($\mu$ playing the role of a homotopy parameter) with a Newton method:
\[F_\mu(x, s, y, v, w) = 0\]
with
\[F_\mu(x, s, y, v, w) =
+\left\{
+\begin{aligned}
+ & \nabla f(x) + A^\top y + \nu + w & (F_1) \\
+ & - y - w & (F_2) \\
+ & g(x) - s & (F_3) \\
+ & X v - \mu e_n & (F_4) \\
+ & S w - \mu e_m & (F_5)
+\end{aligned}
+\right.\]
The Newton step associated to the KKT equations writes
\[\overline{
+\begin{pmatrix}
+ W & 0 & A^\top & - I & 0 \\
+ 0 & 0 & -I & 0 & -I \\
+ A & -I & 0& 0 & 0 \\
+ V & 0 & 0 & X & 0 \\
+ 0 & W & 0 & 0 & S
+\end{pmatrix}}^{K_{1}}
+\begin{pmatrix}
+ \Delta x \\
+ \Delta s \\
+ \Delta y \\
+ \Delta v \\
+ \Delta w
+\end{pmatrix}
+= -
+\begin{pmatrix}
+ F_1 \\ F_2 \\ F_3 \\ F_4 \\ F_5
+\end{pmatrix}\]
The matrix $K_1$ is unsymmetric, but we can obtain an equivalent symmetric system by eliminating the two last rows:
\[\overline{
+\begin{pmatrix}
+ W + \Sigma_x & 0 & A^\top \\
+ 0 & \Sigma_s & -I \\
+ A & -I & 0
+\end{pmatrix}}^{K_{2}}
+\begin{pmatrix}
+ \Delta x \\
+ \Delta s \\
+ \Delta y
+\end{pmatrix}
+= -
+\begin{pmatrix}
+ F_1 + X^{-1} F_4 \\ F_2 + S^{-1} F_5 \\ F_3
+\end{pmatrix}\]
with $\Sigma_x = X^{-1} v$ and $\Sigma_s = S^{-1} w$. The matrix $K_2$, symmetric, has a structure more favorable for a direct linear solver.
In MadNLP, the matrix $K_1$ is encoded as an AbstractUnreducedKKTSystem and the matrix $K_2$ is encoded as an AbstractReducedKKTSystem.
We note that both $K_1$ and $K_2$ depend on the Hessian of the Lagrangian $W$, the Jacobian $A$ and the diagonal matrices $\Sigma_x = X^{1}V$ and $\Sigma_s = S^{-1}W$. Hence, we have to update the KKT system at each iteration of the interior-point algorithm.
By default, MadNLP stores the KKT system as a SparseKKTSystem. The KKT system takes as input a SparseCallback wrapping a given NLPModel nlp. We instantiate the callback cb with the function create_callback:
cb = MadNLP.create_callback(
+ MadNLP.SparseCallback,
+ nlp,
+)
+ind_cons = MadNLP.get_index_constraints(nlp)
(ind_eq = Int64[], ind_ineq = [1, 2], ind_fixed = Int64[], ind_lb = [3, 4], ind_ub = [1], ind_llb = Int64[], ind_uub = [1])
The size of the KKT system depends directly on the problem's characteristics (number of variables, number of of equality and inequality constraints). A SparseKKTSystem stores the Hessian and the Jacobian in sparse (COO) format. The KKT matrix can be factorized using either a dense or a sparse linear solvers. Here we use the solver provided in Lapack:
linear_solver = LapackCPUSolver
LapackCPUSolver
We can instantiate a SparseKKTSystem using the function create_kkt_system:
kkt = MadNLP.create_kkt_system(
+ MadNLP.SparseKKTSystem,
+ cb,
+ ind_cons,
+ linear_solver,
+)
MadNLP.SparseKKTSystem{Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int32}, MadNLP.ExactHessian{Float64, Vector{Float64}}, LapackCPUSolver{Float64, SparseArrays.SparseMatrixCSC{Float64, Int32}}, Vector{Int64}, Vector{Int32}}([0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0], MadNLP.ExactHessian{Float64, Vector{Float64}}(), [0.0, 6.94061061948217e-310, 6.9405756793178e-310, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0], [0.0, 0.0], [0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], 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], [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0], 6, 6), [1, 5, 8, 10, 1, 2, 5, 3, 6, 4, 7, 9, 11, 12, 13], [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0], sparse(Int32[1, 2, 2], [1, 1, 2], [1.0, 2.0, 3.0], 4, 4), [1, 2, 3], [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0], sparse(Int32[1, 2, 1, 2, 1, 2], [1, 1, 2, 2, 3, 4], [1.0, 2.0, 3.0, 4.0, 5.0, 6.0], 2, 4), [1, 3, 2, 4, 5, 6], 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], [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0], 6, 6), [6.94061779244065e-310 6.94061775409207e-310 … 6.9406177952853e-310 6.94061777997597e-310; 6.9406177270829e-310 6.94061772745877e-310 … 6.940617795294e-310 6.94061777632384e-310; … ; 6.94061776319197e-310 6.9406177267813e-310 … 6.9406177757369e-310 6.94061779528094e-310; 6.9406177788483e-310 6.9406177540972e-310 … 6.94061777795307e-310 6.9406177540877e-310], [0.0], -1, Base.RefValue{Int64}(0), Dict{Symbol, Any}(), MadNLP.LapackOptions(MadNLP.BUNCHKAUFMAN), MadNLP.MadNLPLogger(MadNLP.INFO, MadNLP.INFO, nothing)), [1, 2], [3, 4], [1])
Once the KKT system built, one has to initialize it to use it inside the interior-point algorithm:
MadNLP.initialize!(kkt);
3-element Vector{Float64}:
+ 0.0
+ 0.0
+ 0.0
The user can query the KKT matrix inside kkt, simply as
kkt_matrix = MadNLP.get_kkt(kkt)
6×6 SparseArrays.SparseMatrixCSC{Float64, Int32} with 13 stored entries:
+ 1.0 ⋅ ⋅ ⋅ ⋅ ⋅
+ 2.0 5.0 ⋅ ⋅ ⋅ ⋅
+ ⋅ ⋅ 8.0 ⋅ ⋅ ⋅
+ ⋅ ⋅ ⋅ 10.0 ⋅ ⋅
+ 3.0 6.0 9.0 ⋅ 12.0 ⋅
+ 4.0 7.0 ⋅ 11.0 ⋅ 13.0
This returns a reference to the KKT matrix stores internally inside kkt. Each time the matrix is assembled inside kkt, kkt_matrix is updated automatically.
We suppose now we want to refresh the values stored in the KKT system.
The Hessian part of the KKT system can be queried as
hess_values = MadNLP.get_hessian(kkt)
3-element Vector{Float64}:
+ 0.0
+ 0.0
+ 0.0
For a SparseKKTSystem, hess_values is a Vector{Float64} storing the nonzero values of the Hessian. Then, one can update the vector hess_values by using NLPModels.jl:
n = NLPModels.get_nvar(nlp)
+m = NLPModels.get_ncon(nlp)
+x = NLPModels.get_x0(nlp) # primal variables
+l = zeros(m) # dual variables
+
+NLPModels.hess_coord!(nlp, x, l, hess_values)
3-element Vector{Float64}:
+ 2.0
+ 0.0
+ 200.0
Eventually, a post-processing step is applied to refresh all the values internally:
MadNLP.compress_hessian!(kkt)
By default, the function compress_hessian! does nothing. But it can be required for very specific use-case, for instance building internally a Schur complement matrix.
We proceed exaclty the same way to update the values in the Jacobian. One queries the Jacobian values in the KKT system as
jac_values = MadNLP.get_jacobian(kkt)
4-element Vector{Float64}:
+ 0.0
+ 0.0
+ 0.0
+ 0.0
We can refresh the values with NLPModels as
NLPModels.jac_coord!(nlp, x, jac_values)
4-element Vector{Float64}:
+ 0.0
+ 0.0
+ 1.0
+ 0.0
And then applies a post-processing step as
MadNLP.compress_jacobian!(kkt)
Once the Hessian and the Jacobian updated, the algorithm can apply primal and dual regularization terms on the diagonal of the KKT system, to improve the numerical behavior in the linear solver. This operation is implemented inside the regularize_diagonal! function:
pr_value = 1.0
+du_value = 0.0
+
+MadNLP.regularize_diagonal!(kkt, pr_value, du_value)
2-element Vector{Float64}:
+ -0.0
+ -0.0
Once the values updated, one can assemble the resulting KKT matrix. This translates to
MadNLP.build_kkt!(kkt)
By doing so, the values stored inside kkt will be transferred to the KKT matrix kkt_matrix (as returned by the function get_kkt):
kkt_matrix
6×6 SparseArrays.SparseMatrixCSC{Float64, Int32} with 13 stored entries:
+ 4.0 ⋅ ⋅ ⋅ ⋅ ⋅
+ 0.0 202.0 ⋅ ⋅ ⋅ ⋅
+ ⋅ ⋅ 2.0 ⋅ ⋅ ⋅
+ ⋅ ⋅ ⋅ 2.0 ⋅ ⋅
+ 0.0 0.0 -1.0 ⋅ 0.0 ⋅
+ 1.0 0.0 ⋅ -1.0 ⋅ 0.0
Internally, a SparseKKTSystem stores the KKT system in a sparse COO format. When build_kkt! is called, the sparse COO matrix is transferred to SparseMatrixCSC if the linear solver is sparse, or alternatively to a Matrix if the linear solver is dense.
The KKT system stores only the lower-triangular part of the KKT system, as it is symmetric.
Now the KKT system is assembled in a matrix $K$ (here stored in kkt_matrix), we want to solve a linear system $K x = b$, for instance to evaluate the next descent direction. To do so, we use the linear solver stored internally inside kkt (here an instance of LapackCPUSolver).
We start by factorizing the KKT matrix $K$:
MadNLP.factorize!(kkt.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], [4.0, 0.0, 0.0, 1.0, 202.0, 0.0, 0.0, 2.0, -1.0, 2.0, -1.0, 0.0, 0.0], 6, 6), [4.0 0.0 … 0.0 0.0; 0.0 202.0 … 0.0 0.0; … ; 0.0 0.0 … -0.5 0.0; 0.25 0.0 … -0.0 -0.75], [384.0, 8.487983164e-314, 6.9529018960701e-310, 1.5e-323, 2.24726905e-316, 2.4e-322, 2.2472556e-316, 2.4e-322, 2.24724218e-316, 2.4e-322 … 2.2472888e-316, 8.487983164e-314, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 384, Base.RefValue{Int64}(0), Dict{Symbol, Any}(:ipiv => [1, 2, 3, 4, 5, 6]), MadNLP.LapackOptions(MadNLP.BUNCHKAUFMAN), MadNLP.MadNLPLogger(MadNLP.INFO, MadNLP.INFO, nothing))
By default, MadNLP uses a LBL factorization to decompose the symmetric indefinite KKT matrix.
Once the KKT matrix has been factorized, we can compute the solution of the linear system with a backsolve. The function takes as input a AbstractKKTVector, an object used to do algebraic manipulation with a AbstractKKTSystem. We start by instantiating two UnreducedKKTVector (encoding respectively the right-hand-side and the solution):
b = MadNLP.UnreducedKKTVector(kkt)
+fill!(MadNLP.full(b), 1.0)
+x = copy(b)
MadNLP.UnreducedKKTVector{Float64, Vector{Float64}, Vector{Int64}}([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0], [1.0, 1.0], [1.0], [1.0, 1.0], [1.0, 1.0], [1.0])
The right-hand-side encodes a vector of 1:
MadNLP.full(b)
9-element Vector{Float64}:
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
We solve the system $K x = b$ using the solve! function:
MadNLP.solve!(kkt, x)
+MadNLP.full(x)
9-element Vector{Float64}:
+ 0.3333333333333333
+ 0.0049504950495049506
+ -1.0
+ -0.6666666666666666
+ -2.0
+ -1.3333333333333333
+ -1.0
+ -1.0
+ 1.0
We verify that the solution is correct by multiplying it on the left with the KKT system, using mul!:
mul!(b, kkt, x) # overwrite b!
+MadNLP.full(b)
9-element Vector{Float64}:
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
We recover a vector filled with 1, which was the initial right-hand-side.
