KKT systems

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:

1. Store the values of the Hessian of the Lagrangian and of the Jacobian.
2. Assemble the corresponding KKT matrix $K$.

A brief look at the math

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.

Assembling a KKT system, step by step

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])

Initializing a KKT system

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}}(), [6.94690576689115e-310, 6.94690576689194e-310, 6.94690576689273e-310, 6.9469057668935e-310], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0], [5.0e-324, 0.0], [6.9468751406827e-310], [5.0e-324, 0.0], [6.9467776208939e-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], 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), [1.0e-323 0.0 … 0.0 0.0; 5.0e-324 5.0e-324 … 0.0 0.0; … ; 0.0 5.0e-324 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [6.94690796608406e-310], -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.

Updating a KKT system

We suppose now we want to refresh the values stored in the KKT system.

Updating the values of the Hessian

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)

Updating the values of the Jacobian

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)

Updating the values of the diagonal matrices

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

Assembling the KKT matrix

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.

Solution of the KKT system

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, 1.897e-321, 3.06482477e-316, 0.0, 2.045e-321, 1.5174104e-316, 5.0e-324, 2.8e-322, 3.06468722e-316, 0.0  …  4.0e-323, 1.7740988276e-314, 1.18419314e-316, 8.74e-322, 2.97984825e-316, 1.44151043e-316, 6.94691429196146e-310, 4.43e-321, 1.71235564e-316, 6.37e-322], 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.