# Zero pivot error in LDLt factorization for large matrices

I encountered a problem with the LDLt factorization in the finite-difference solution of a 2D Poisson equation: u_xx + u_yy = -1 on 0 < x < 1 and 0 < y < 1 with u = 0 Dirichlet boundary conditions. The factorization works fine for small problems (e.g. nx = ny = 64), but I get told my matrix has one or more zero pivots for bigger problems (e.g. nx = ny = 128). If I instead use an LU factorization, the solution works fine also for bigger problems. Any idea what could be going on?

Here’s the code:

``````# grid points
nx = 128
ny = 128

# grid spacing
dx = 1/(nx+1)
dy = 1/(ny+1)

# setup matrix
idx(j,i) = (i-1)*ny+j
B = spzeros(nx*ny,nx*ny)
for i = 1:nx
for j = 1:ny
# x-derivative
if i < nx
B[idx(j,i),idx(j,i+1)] += 1/dx^2
end
B[idx(j,i),idx(j,i)] -= 1/dx^2
B[idx(j,i),idx(j,i)] -= 1/dx^2
if i > 1
B[idx(j,i),idx(j,i-1)] += 1/dx^2
end
# y-derivative
if j < ny
B[idx(j,i),idx(j+1,i)] += 1/dy^2
end
B[idx(j,i),idx(j,i)] -= 1/dy^2
B[idx(j,i),idx(j,i)] -= 1/dy^2
if j > 1
B[idx(j,i),idx(j-1,i)] += 1/dy^2
end
end
end

# solve
u = -reshape(B\ones(ny*nx), (ny,nx))
``````

And here is the error I’m getting:

``````ERROR: LoadError: ArgumentError: matrix has one or more zero pivots
in #ldltfact!#10(::Float64, ::Function, ::Base.SparseArrays.CHOLMOD.Factor{Float64}, ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./sparse/cholmod.jl:1350
in (::Base.LinAlg.#kw##ldltfact!)(::Array{Any,1}, ::Base.LinAlg.#ldltfact!, ::Base.SparseArrays.CHOLMOD.Factor{Float64}, ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./<missing>:0
in #ldltfact#12(::Float64, ::Array{Int64,1}, ::Function, ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./sparse/cholmod.jl:1386
in #ldltfact#13(::Array{Any,1}, ::Function, ::Hermitian{Float64,SparseMatrixCSC{Float64,Int64}}) at ./sparse/cholmod.jl:1426
in \(::SparseMatrixCSC{Float64,Int64}, ::Array{Float64,1}) at ./sparse/linalg.jl:871
``````

The LDLt factorization we use is from CHOLMOD and, in contrast to UMFPACK which the sparse LU is based on, it doesn’t do any pivoting during the numerical factorization. Consequently, the `lufact` can factorize matrices for which `ldltfact` fails. A simple example is

``````julia> A = sparse([0 1; 1 1])
2×2 sparse matrix with 3 Int64 stored entries:
[2, 1]  =  1
[1, 2]  =  1
[2, 2]  =  1

julia> lufact(A)
UMFPACK LU Factorization of a (2, 2) sparse matrix
Ptr{Void} @0x00007f80c2240d80

julia> ldltfact(A)
ERROR: ArgumentError: matrix has one or more zero pivots
Stacktrace:
[1] #ldltfact!#12(::Float64, ::Function, ::Base.SparseArrays.CHOLMOD.Factor{Float64}, ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./sparse/cholmod.jl:1430
[2] (::Base.LinAlg.#kw##ldltfact!)(::Array{Any,1}, ::Base.LinAlg.#ldltfact!, ::Base.SparseArrays.CHOLMOD.Factor{Float64}, ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./<missing>:0
[3] #ldltfact#14(::Float64, ::Array{Int64,1}, ::Function, ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./sparse/cholmod.jl:1473
[4] #ldltfact#15(::Array{Any,1}, ::Function, ::SparseMatrixCSC{Int64,Int64}) at ./sparse/cholmod.jl:1516
[5] ldltfact(::SparseMatrixCSC{Int64,Int64}) at ./sparse/cholmod.jl:1516
``````

It is also worth noting that pivoting is not needed for well-conditioned positive definite systems.

The OP’s matrix is negative definite and well-conditioned.

So if you solve `(-B)\b` above, `ldltfact` should work with no trouble.

Perhaps the documentation for left-divide needs tweaking? Currently it says

When A is sparse, a similar polyalgorithm is used. For indefinite matrices,
the LDLt factorization does not use pivoting during the numerical
factorization and therefore the procedure can fail even for invertible
matrices.

It is actually a bug in Julia. I’d forgotten this one https://github.com/JuliaLang/julia/pull/19045. The fix is merged and should be included in next bug release of Julia. With that, the example works just fine. I’ve just tried it on master.

2 Likes