Zero pivot error in LDLt factorization for large matrices

Specific Domains Numerics
1

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
 in include_from_node1(::String) at ./loading.jl:488
2

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
3

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.

4

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