Bordered System Solvers

Note

BorderedLinearSolve is a sub-library of LinearSolve.jl. Using these solvers requires adding it: using BorderedLinearSolve.

Solves for $v$ and $σ$ in

\[\begin{pmatrix} J & b \\ c^{T} & d \end{pmatrix} \begin{pmatrix} v \\ σ \end{pmatrix} = \begin{pmatrix} f \\ g \end{pmatrix}\]

where J is n × n and the border b, c is n × m with m small. Systems of this shape come up in continuation and bifurcation methods, where J is the Jacobian of a large system and the border carries a few extra constraints.

BorderingBLS solves by block elimination against J itself, so its inner solves share one LinearSolve cache, which carries over between steps where J changes only slightly. DirectBLS assembles the whole matrix and solves it, which is fastest for a small dense J.

using BorderedLinearSolve, LinearSolve, LinearAlgebra
n = 100
J = rand(n, n) + n * I
b, c, f = rand(n), rand(n), rand(n)
prob = BorderedLinearProblem(J, b, c, 1.0, f, 2.0)
sol = solve(prob, BorderingBLS(LUFactorization()))
[J b; c' 1.0] * [sol.v; sol.σ] ≈ [f; 2.0]

true

BorderedLinearSolve.BorderedLinearSolve — Module

BorderedLinearSolve

Solvers for bordered linear systems

┌      ┐┌   ┐   ┌   ┐
│ J  b ││ v │ = │ f │
│ c' d ││ σ │   │ g │
└      ┘└   ┘   └   ┘

where J is n x n and the border b, c is n x m with m small (usually 1). Systems of this shape come up in continuation and bifurcation work, where J is the Jacobian of a large system and the border carries a handful of extra constraints.

The point of solving them this way rather than assembling the whole (n+m) x (n+m) matrix is that every inner solve uses the sameJ, so one LinearSolve cache serves them all, and it carries over between outer steps where J changes only slightly. That is what makes factorization or Krylov subspace reuse available here.

See BorderedLinearProblem, BorderingBLS and DirectBLS.

source

BorderedLinearSolve.BorderedLinearProblem — Type

BorderedLinearProblem(J, b, c, d, f, g)

The bordered system

┌      ┐┌   ┐   ┌   ┐
│ J  b ││ v │ = │ f │
│ c' d ││ σ │   │ g │
└      ┘└   ┘   └   ┘

f is the state, and may be any array. For a 2D problem it is natural to keep the state as a matrix, and nothing here requires it to be flattened into a vector.

b and c live in the same space as the state, so for a width-1 border they have the same shape as f, and d and g are scalars.

For a border of width m > 1, pass b and c as a tuple (or vector) of m states, d as an m x m matrix and g as a length-m vector. The width is taken from that structure and never inferred from an array's shape, since a state can itself be a matrix and the two cannot be told apart otherwise.

J is only ever applied and solved against, so it can be a matrix or a matrix-free operator; no assumption is made about its shape.

source

CommonSolve.solve — Method

solve(prob::BorderedLinearProblem, alg)

Solve the bordered system in one shot. Use init plus solve! instead when the same J is solved against repeatedly, which is the case this package is for.

source

BorderedLinearSolve.BorderingBLS — Type

BorderingBLS(; inner = nothing, k = 0, tol = 1e-12)

Block elimination. Solves J x₁ = f and J Δᵢ = bᵢ, then

solve (d - c'Δ) σ = (g - c'x₁)
v = x₁ - Δ σ

Every inner solve uses the same J, so they share one LinearSolve cache; inner is the LinearSolve algorithm for them (nothing picks the default).

k is the number of correction sweeps in the sense of Govaerts' BEC+k: after the first pass the true residual of the bordered system is formed and used to solve for a correction, which recovers accuracy when J is close to singular. That happens exactly where continuation is interesting, near a fold, so it is worth having, but each sweep costs another pass and the default of 0 skips it. A sweep is skipped early once the residual is below tol.

Only dot, scaling and addition are used on the state, so the state can be any array, not just a vector.

Reference

Govaerts, W. "Stable Solvers and Block Elimination for Bordered Systems." SIAM Journal on Matrix Analysis and Applications 12, no. 3 (1991): 469-83.

source

BorderedLinearSolve.DirectBLS — Type

DirectBLS(; inner = nothing)

Assemble the full bordered matrix and hand it to LinearSolve.

This gives up the cache reuse that BorderingBLS exists for, and for a sparse J the assembled matrix has a dense row and column that the factorization has to carry. It is here because it is the reference the bordering result should be checked against, and because for a small J it is simply faster than repeated inner solves.

Assembling requires a J that is an actual matrix and a state that can be flattened, so unlike BorderingBLS this is not available for a matrix-free J.

source