Nearest Global Section: Iterative Method

This example demonstrates the iterative conjugate-gradient (CG) backend for projecting a 0-cochain onto the nearest global section. The CG solver is stochastic in its convergence behaviour and may occasionally fail to find a global section for ill-conditioned sheaves or small tolerances. The example below wraps the call in a try/catch so that a docs build is never broken by a solver failure.

For production use, prefer nearest_global_section with method=:ldl (the direct LDLt backend), which is deterministic and typically much faster for small-to-medium sheaves.

using CellularSheaves
using LinearAlgebra

Building the Sheaf

We build the same identity-restriction-map sheaf on a 4-cycle that appears in the main Code Example, but this time with 3-dimensional stalks so the sheaf Laplacian is slightly larger.

n = 4
stalk_dim = 3

s = EuclideanSheaf{Float64}(Int[])
for _ in 1:n
    add_vertex_stalk!(s, stalk_dim)
end
for i in 1:n
    v1 = i
    v2 = mod1(i + 1, n)
    rm = Matrix{Float64}(I, stalk_dim, stalk_dim)
    add_sheaf_edge!(s, v1, v2, rm, rm)
end

println(s)

A network sheaf with 4 vertex stalks and 4 edge stalks.

Running the Iterative Solver

nearest_global_section with method=:cg uses the Krylov CG method from Krylov.jl. Because the edge-space operator $E = d d^\top$ is positive semi-definite (not strictly positive definite when the sheaf has non-trivial global sections), the iterative solver can sometimes stagnate. We therefore wrap the call and print any error rather than letting it propagate.

x0 = rand(sum(vertex_stalks(s)))

gs = try
    nearest_global_section(s, x0; method=:cg, verbose=true)
catch e
    println("Iterative solver did not converge: ", e)
    nothing
end

if gs !== nothing
    d = coboundary_map(s)
    println("‖d · gs‖ = ", norm(d * gs))
else
    println("No global section returned by the iterative solver.")
end

SimpleStats
 niter: 2
 solved: true
 inconsistent: false
 indefinite: false
 npcCount: 0
 residuals: []
 Aresiduals: []
 κ₂(A): []
 allocation timer: 0.59μs
 timer: 10.88μs
 status: solution good enough given atol and rtol

‖d · gs‖ = 2.220446049250313e-16

Comparison with the Direct Method

For reference, method=:ldl uses a deterministic direct backend and is typically faster on small-to-medium sheaves. Internally it may fall back to a pseudoinverse if the LDL factorization fails, and in rare edge cases either path can still throw, so we guard this call as well.

gs_ldl = try
    nearest_global_section(s, x0; method=:ldl)
catch e
    println("Direct LDL backend failed: ", e)
    nothing
end

if gs_ldl !== nothing
    d = coboundary_map(s)
    println("‖d · gs_ldl‖ = ", norm(d * gs_ldl))
else
    println("No global section returned by the direct solver.")
end

‖d · gs_ldl‖ = 3.8459253727671276e-16