I would like to solve a system of linear equations given in a vector form as
The matrix Q is symmetric and positive definite. The matrix A has full row rank. The matrix is called KKT matrix and the full equation is just an optimality condition for a linear equality-constrained quadratic optimization. I would like to solve the problem for sparse data. I am rather confused as for which tool/solver to use (and perhaps even how to use it). I do solve the problem successfuly but I wanted some more insight and overview anyway.
I will build the KKT matrix K and the right hand side vector k first:
using BlockArrays
using LinearAlgebra
using SparseArrays
n = 2
m = 1
A = rand(n,n)
B = rand(n,m)
x₀ = [1.0, 3.0]
N = 100
s = [1.0, 2.0]
q = [1.0, 2.0]
r = [1.0]
S = diagm(0=>s)
Q = diagm(0=>q)
R = diagm(0=>r)
Qbar = BlockArray(spzeros(N*n,N*n),repeat([n],N),repeat([n],N))
for i=1:(N-1)
Qbar[Block(i,i)] = Q
end
Qbar[Block(N,N)] = S
Rbar = BlockArray(spzeros(N*m,N*m),repeat([m],N),repeat([m],N))
for i=1:N
Rbar[Block(i,i)] = R
end
Qtilde = blockdiag(sparse(Qbar),sparse(Rbar))
Bbar = BlockArray(spzeros(N*n,N*m),repeat([n],N),repeat([m],N))
for i=1:N
Bbar[Block(i,i)] = B
end
Abar = BlockArray(sparse(-1.0*I,n*N,n*N),repeat([n],N),repeat([n],N))
for i=2:N
Abar[Block(i,(i-1))] = A
end
Atilde = sparse([Abar Bbar])
A0bar = spzeros(n*N,n)
A0bar[1:n,1:n] = A
btilde = A0bar*sparse(x₀)
K = [Qtilde Atilde'; Atilde spzeros(size(Atilde,1),size(Atilde,1))] # Sparse KKT matrix.
k = [spzeros(size(Qtilde,1)); -btilde] # Right hand side of the KKT system
This is how K looks like
julia> K
500×500 SparseMatrixCSC{Float64, Int64} with 1892 stored entries:
⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠈⠙⠲⢤⣀⠀⠀⠀⠀⠀⠀⠀⠈
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠈⠙⠲⢤⣀⠀⠀⠀⠀
⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠈⠙⠲⢤⣀
⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⢳⠀
Now, the provided code restricts the upper left block (in the code named as) Qtilde of the KKT matrix to a diagonal matrix. I know that this could be exploited because its Cholesky factorization is for free, which suggests some elimination scheme. But let’s pretend that this left upper block in the KKT matrix is a general symmetric positive definite matrix (but sparse) and aim at a general procedure.
The immediate solution in Julia is to invoke the backslash:
julia> xtildeλ = K\Array(k)
500-element Vector{Float64}:
0.9038685869909274
1.9104601797113756
⋮
1.7348358419940635e-14
4.2401702681320544e-14
Upon checking the norm of the residual, the solution seems to be accurate:
julia> norm(K*xtildeλ-k)
1.5176942235117913e-14
But I was wondering if I the solver invoked by the backslash takes advantage of the fact that K is symmetric (although it is indefinite). Therefore I tried to do LDLᵀ factorization:
julia> LDLT = ldlt(K,check=true)
ERROR: ZeroPivotException: factorization encountered one or more zero pivots. Consider switching to a pivoted LU factorization.
Stacktrace:
[1] #ldlt!#10
@ ~/julia-1.7.2/share/julia/stdlib/v1.7/SuiteSparse/src/cholmod.jl:1308 [inlined]
[2] ldlt(A::SuiteSparse.CHOLMOD.Sparse{Float64}; shift::Float64, check::Bool, perm::Nothing)
@ SuiteSparse.CHOLMOD ~/julia-1.7.2/share/julia/stdlib/v1.7/SuiteSparse/src/cholmod.jl:1347
[3] #ldlt#13
@ ~/julia-1.7.2/share/julia/stdlib/v1.7/SuiteSparse/src/cholmod.jl:1396 [inlined]
[4] top-level scope
@ REPL[27]:1
Can anyone help me understand what is going on?
Having failed with LDLᵀ factorization, I resorted to LU factorization:
julia> LU = lu(K,check=true)
SuiteSparse.UMFPACK.UmfpackLU{Float64, Int64}
L factor:
500×500 SparseMatrixCSC{Float64, Int64} with 1590 stored entries:
⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀
⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄
U factor:
500×500 SparseMatrixCSC{Float64, Int64} with 1692 stored entries:
⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈
⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦
julia> xtildeλ_lu = LU\Array(k)
500-element Vector{Float64}:
0.9038685869909274
1.9104601797113756
⋮
1.7348358419940635e-14
4.2401702681320544e-14
The solution seems accurate as well and after timing it is seems even a bit faster compared to the solver based on backslash (of course the factorization is included in the measurement of time). But I was wondering if there is any other solver that takes the symmetry into consideration.
What else is out there? I have also tried MINRES from IterativeSolvers.jl:
julia> using IterativeSolvers
julia> xtildeλ_minres = minres(K,Array(k))
500-element Vector{Float64}:
0.9038685869833563
1.9104601797222807
0.5444506437441838
1.425402384616977
0.42137309546601803
0.9619176467853587
0.278065185300917
⋮
0.0
0.0
0.0
0.0
0.0
0.0
0.0
julia> norm(K*xtildeλ_minres-k)
5.0481687737604944e-8
I also tried LDLFactorizations.jl package:
julia> LDLT = ldl(K)
LDLFactorizations.LDLFactorization{Float64, Int64, Int64, Int64}(true, false, false, 500, [4, 5, 4, 5, 6, 7, 9, 9, 10, 11 … 497, 493, 494, 495, 496, 499, 499, 499, 500, -1], [1, 1, 2, 3, 2, 2, 2, 2, 3, 2 … 0, 0, 2, 2, 3, 2, 2, 2, 1, 0], [4, 5, 5, 7, 7, 9, 10, 10, 12, 12 … 491, 492, 495, 496, 500, 500, 500, 500, 500, 500], [199, 200, 300, 499, 500, 198, 197, 299, 498, 497 … 4, 302, 201, 301, 2, 1, 3, 202, 304, 303], [496, 495, 497, 491, 487, 486, 482, 481, 477, 476 … 25, 24, 20, 19, 15, 14, 10, 9, 4, 5], [1, 2, 3, 5, 8, 10, 12, 14, 16, 19 … 1079, 1081, 1083, 1085, 1088, 1090, 1092, 1094, 1095, 1095], Int64[], Int64[], [4, 5, 4, 5, 5, 6, 7, 6, 7, 7 … 85482544, -5353651815063748608, 112, 110808136, 17179869185, 85042200, 85042600, 110808160, 11403440, 0], [-1.0, -0.5, 0.11475043095954407, 0.04739161352364707, 0.0053675302547291385, -0.8362922137587595, -0.35878159855262715, -1.157090992397963, -1.9477015745917226, 0.4246762957067049 … 1.93e-322, 4.0589489e-316, 5.09783336e-316, 4.31171524e-316, 4.2439915824e-314, 5.09783336e-316, 0.0, 4.5365092e-316, 2.5463949504e-313, 1.14e-322], [1.0, 2.0, 1.0, -1.013167661405401, -0.5022167752859968, 3.3809916604166665, 2.425837675665596, 1.0, -0.37236261123419273, -0.3481253955555499 … 41.30042883848547, 0.0, 3.93298033e-316, 2.4e-322, NaN, 2.016e-320, 2.0316e-320, 9.54e-322, 3.4844387e-316, 4.3504235e-316], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 2.0316e-320, 1.03e-321, 3.4772269e-316, 3.4872655e-316, 0.0, 3.4872655e-316, 0.0, 5.0e-324], [489, 489, 6, 0, 0, 78199376, 103186912, 0, 0, 139855653765121 … 0, 13814488, 4112, 224, 102002560, 83321744, 76587200, 488, 489, 490], 0.0, 0.0, 0.0, 500)
julia> xtildeλ_ldl = LDLT\Array(k)
ERROR: LDLFactorizations.SQDException("LDL' factorization was not computed or failed")
Stacktrace:
[1] \(LDL::LDLFactorizations.LDLFactorization{Float64, Int64, Int64, Int64}, b::Vector{Float64})
@ LDLFactorizations ~/.julia/packages/LDLFactorizations/UlVgR/src/LDLFactorizations.jl:579
[2] top-level scope
@ REPL[43]:1
Obviously failing here.
Any other solver or perhaps even tip on more appropriate usage of the ones that I have mentioned?