Support pinv in Woodbury and SymWoodbury (#50)

When building a `Woodbury` or `SymWoodbury` matrix, we precompute

    (C⁻¹ + V * (A\U))⁻¹

However, if the final `W` is rank-deficient, this inverse does not
exist. To handle this case, we now support an optional keyword argument
`use_pinv::Bool=false` to use the pseudo-inverse instead.

Also updates cache actions and adds dependabot support.

Author: timholy

.github/dependabot.yml

@@ -0,0 +1,7 @@
+# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "monthly"

.github/workflows/ci.yml

@@ -15,19 +15,19 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
+          - 'min'
           - '1'
         os:
           - ubuntu-latest
         arch:
           - x64
     steps:
       - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
+      - uses: actions/cache@v4
         env:
           cache-name: cache-artifacts
         with:

Project.toml

@@ -1,6 +1,6 @@
 name = "WoodburyMatrices"
 uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
-version = "1.0.0"
+version = "1.1.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -12,7 +12,7 @@ Random = "1"
 LinearAlgebra = "1"
 SparseArrays = "1"
 Test = "1"
-julia = "1.6"
+julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"

README.md

@@ -41,3 +41,9 @@ If passed the keyword argument `allocatetmp=true`, (Sym)Woodbury allocates inter
 This eliminates memory allocation for these common operations.
 
 However, using the same `W` across multiple threads can lead to race conditions. Hence, this optimization is opt-in and should only be used if you know it is safe.
+
+## Rank-deficiency
+
+While `A` is required to be invertible, rank-deficient `W` are supported if you supply the `use_pinv=true` keyword argument to the constructor,
+so that we compute `pinv(inv(C) + V*(A\U))` instead of the standard inverse.
+Then `x = W \ b` will effectively use the pseudo-inverse of `W`.

src/WoodburyMatrices.jl

@@ -12,6 +12,9 @@ abstract type AbstractWoodbury{T} <: Factorization{T} end
 safeinv(A) = inv(A)
 safeinv(A::SparseMatrixCSC) = safeinv(Matrix(A))
 
+safepinv(A) = pinv(A)
+safepinv(A::SparseMatrixCSC) = safepinv(Matrix(A))
+
 include("woodbury.jl")
 include("symwoodbury.jl")
 include("sparsefactors.jl")

src/symwoodbury.jl

@@ -13,7 +13,7 @@ struct SymWoodbury{T,AType,BType,DType,DpType} <: AbstractWoodbury{T}
 end
 
 """
-    W = SymWoodbury(A, B, D; allocatetmp::Bool=false)
+    W = SymWoodbury(A, B, D; allocatetmp::Bool=false, use_pinv::Bool=false)
 
 Represent a matrix of the form `W = A + BDBᵀ`, where `A` and `D` are symmetric.
 Equations `Wx = b` will be solved using the
@@ -26,9 +26,9 @@ semidefinite matrix).
 This checks that `A` is symmetric; you can elide the check by passing a `Symmetric` matrix
 or factorization.
 
-See also [Woodbury](@ref), where `allocatetmp` is explained.
+See also [Woodbury](@ref), where `allocatetmp` and `use_pinv` are explained.
 """
-function SymWoodbury(A, B::AbstractVecOrMat, D; allocatetmp::Bool=false)
+function SymWoodbury(A, B::AbstractVecOrMat, D; allocatetmp::Bool=false, use_pinv::Bool=false)
     @noinline throwdmm(B, D, A) = throw(DimensionMismatch("Sizes of B ($(size(B))) and/or D ($(size(D))) are inconsistent with A ($(size(A)))"))
 
     n = size(A, 1)
@@ -40,7 +40,8 @@ function SymWoodbury(A, B::AbstractVecOrMat, D; allocatetmp::Bool=false)
         issymmetric(A) || throw(ArgumentError("A must be symmetric"))
     end
     issymmetric(D) || throw(ArgumentError("D must be symmetric"))
-    Dp = safeinv(safeinv(D) .+ B'*(A\B))
+    Dpinv = safeinv(D) .+ B'*(A\B)
+    Dp = use_pinv ? safepinv(Dpinv) : safeinv(Dpinv)
     # temporary space for allocation-free solver (vector RHS only)
     T = typeof(float(zero(eltype(A)) * zero(eltype(B)) * zero(eltype(D))))
     if allocatetmp

src/woodbury.jl

@@ -14,7 +14,7 @@ struct Woodbury{T,AType,UType,VType,CType,CpType} <: AbstractWoodbury{T}
 end
 
 """
-    W = Woodbury(A, U, C, V; allocatetmp::Bool=false)
+    W = Woodbury(A, U, C, V; allocatetmp::Bool=false, use_pinv::Bool=false)
 
 Represent a matrix `W = A + UCV`.
 Equations `Wx = b` will be solved using the
@@ -23,18 +23,20 @@ Equations `Wx = b` will be solved using the
 If your main goal is to solve equations, it's often advantageous to supply
 `A` as a factorization (e.g., `Woodbury(lu(A), U, C, V)`).
 
+If `W` is rank-deficient or nearly so, setting `use_pinv` to true will use
+the pseudoinverse in the Woodbury formula to improve numerical stability.
+
 If `allocatetmp` is true, temporary storage used for intermediate steps in
 multiplication and division will be allocated.
 
 !!! warning
     If you'll use the same `W` in multiple threads, you should use `allocatetmp=false`
     or risk data races.
 
-
 See also [SymWoodbury](@ref).
 
 """
-function Woodbury(A, U::AbstractMatrix, C, V::AbstractMatrix; allocatetmp::Bool=false)
+function Woodbury(A, U::AbstractMatrix, C, V::AbstractMatrix; allocatetmp::Bool=false, use_pinv::Bool=false)
     @noinline throwdmm1(U, V, A) = throw(DimensionMismatch("Sizes of U ($(size(U))) and/or V ($(size(V))) are inconsistent with A ($(size(A)))"))
     @noinline throwdmm2(k) = throw(DimensionMismatch("C should be $(k)x$(k)"))
 
@@ -48,7 +50,8 @@ function Woodbury(A, U::AbstractMatrix, C, V::AbstractMatrix; allocatetmp::Bool=
             throwdmm2(k)
         end
     end
-    Cp = safeinv(safeinv(C) .+ V*(A\U))
+    Cpinv = safeinv(C) .+ V*(A\U)
+    Cp = use_pinv ? safepinv(Cpinv) : safeinv(Cpinv)
     # temporary space for allocation-free solver (vector RHS only)
     T = typeof(float(zero(eltype(A)) * zero(eltype(U)) * zero(eltype(C)) * zero(eltype(V))))
     if allocatetmp

test/symwoodbury.jl

@@ -26,7 +26,7 @@ for elty in (Float32, Float64, ComplexF32, ComplexF64, Int)
     ε = eps(abs2(float(one(elty))))
     A = Diagonal(a)
 
-    for W in (SymWoodbury(A, B, D), SymWoodbury(A, B, D; allocatetmp=true), SymWoodbury(A, B[:,1][:], 2.))
+    for W in (SymWoodbury(A, B, D), SymWoodbury(A, B, D; allocatetmp=true), SymWoodbury(A, B[:,1][:], 2.), SymWoodbury(A, B, D; use_pinv=true))
 
         @test issymmetric(W)
         F = Matrix(W)
@@ -131,48 +131,48 @@ end
 A = Diagonal((rand(n)))
 B = sprandn(n,2,1.)
 D = sprandn(2,2,1.); D = (D + D')/2
-W = SymWoodbury(A, B, D)
-v = randn(n)
-vdiag = Diagonal(v)
-V = randn(n,1)
-
-@test size(W) == (n,n)
-@test size(W,1) == n
-@test size(W,2) == n
-
-@test inv(W)*v ≈ inv(Matrix(W))*v
-@test (2*W)*v ≈ 2*(W*v)
-@test (W'W)*v ≈ Matrix(W)*(Matrix(W)*v)
-@test (W*W)*v ≈ Matrix(W)*(Matrix(W)*v)
-@test (W*W')*v ≈ Matrix(W)*(Matrix(W)*v)
-
-@test inv(W)*vdiag ≈ W\vdiag
-@test W\vdiag ≈ W\Matrix(vdiag)
-@test inv(W)*vdiag ≈ inv(Matrix(W))*vdiag
-
-@test inv(W)*V ≈ inv(Matrix(W))*V
-@test (2*W)*V ≈ 2*(W*V)
-@test (W'W)*V ≈ Matrix(W)*(Matrix(W)*V)
-@test (W*W)*V ≈ Matrix(W)*(Matrix(W)*V)
-@test (W*W')*V ≈ Matrix(W)*(Matrix(W)*V)
-
-@test Matrix(2*W) ≈ 2*Matrix(W)
-@test Matrix(W*2) ≈ 2*Matrix(W)
-R = Symmetric(rand(size(W)...))
-# Not sure when this got fixed, but check that failures are for a "known" reason
-canadd = try (W + R; true;) catch err; (@test err isa MethodError && err.f == ldiv!; false;) end
-if canadd
-    @test Matrix(W + R) ≈ Matrix(W) + R
-    @test Matrix(R + W) ≈ Matrix(W) + R
-else
-    @test_broken Matrix(W + R) ≈ Matrix(W) + R
-    @test_broken Matrix(R + W) ≈ Matrix(W) + R
-end
-Wm = SymWoodbury(A, Matrix(B), D)
-@test Matrix(Wm + R) ≈ Matrix(Wm) + R
-@test Matrix(R + Wm) ≈ Matrix(Wm) + R
+for W in (SymWoodbury(A, B, D), SymWoodbury(A, B, D; use_pinv=true))
+    v = randn(n)
+    vdiag = Diagonal(v)
+    V = randn(n,1)
+    @test size(W) == (n,n)
+    @test size(W,1) == n
+    @test size(W,2) == n
+
+    @test inv(W)*v ≈ inv(Matrix(W))*v
+    @test (2*W)*v ≈ 2*(W*v)
+    @test (W'W)*v ≈ Matrix(W)*(Matrix(W)*v)
+    @test (W*W)*v ≈ Matrix(W)*(Matrix(W)*v)
+    @test (W*W')*v ≈ Matrix(W)*(Matrix(W)*v)
+
+    @test inv(W)*vdiag ≈ W\vdiag
+    @test W\vdiag ≈ W\Matrix(vdiag)
+    @test inv(W)*vdiag ≈ inv(Matrix(W))*vdiag
+
+    @test inv(W)*V ≈ inv(Matrix(W))*V
+    @test (2*W)*V ≈ 2*(W*V)
+    @test (W'W)*V ≈ Matrix(W)*(Matrix(W)*V)
+    @test (W*W)*V ≈ Matrix(W)*(Matrix(W)*V)
+    @test (W*W')*V ≈ Matrix(W)*(Matrix(W)*V)
+
+    @test Matrix(2*W) ≈ 2*Matrix(W)
+    @test Matrix(W*2) ≈ 2*Matrix(W)
+    R = Symmetric(rand(size(W)...))
+    # Not sure when this got fixed, but check that failures are for a "known" reason
+    canadd = try (W + R; true;) catch err; (@test err isa MethodError && err.f == ldiv!; false;) end
+    if canadd
+        @test Matrix(W + R) ≈ Matrix(W) + R
+        @test Matrix(R + W) ≈ Matrix(W) + R
+    else
+        @test_broken Matrix(W + R) ≈ Matrix(W) + R
+        @test_broken Matrix(R + W) ≈ Matrix(W) + R
+    end
+    Wm = SymWoodbury(A, Matrix(B), D)
+    @test Matrix(Wm + R) ≈ Matrix(Wm) + R
+    @test Matrix(R + Wm) ≈ Matrix(Wm) + R
 
-@test transpose(W)*v ≈ transpose(Matrix(W))*v
+    @test transpose(W)*v ≈ transpose(Matrix(W))*v
+end
 
 # Factorization for A
 A = SymTridiagonal(rand(5).+2, rand(4))
@@ -205,6 +205,7 @@ W1 = SymWoodbury(A, B, D)
 @test_throws ArgumentError SymWoodbury(rand(5,5),rand(5,2),rand(2,2))
 
 # Display
+W = SymWoodbury(A, B, D)
 iob = IOBuffer()
 show(iob, MIME("text/plain"), W)
 str = String(take!(iob))
@@ -243,4 +244,23 @@ end
     @test all(logabsdet(W) .≈ logabsdet(Matrix(W)))
 end
 
+@testset "pinv" begin
+    # Build a matrix W that is equivalent to
+    #   A = [Diagonal(c) ones(length(c))]
+    #   A'A
+    # This is rank-deficient, so inv should fail but pinv should work
+    c = [-1.0, -2.0, -3.0, -4.0]
+    d = c.^2
+    push!(d, length(c))
+    U = [c zeros(length(c)); 0 1]
+    @test_throws SingularException SymWoodbury(Diagonal(d), U, [0 1; 1 0])
+    W = SymWoodbury(Diagonal(d), U, [0 1; 1 0]; use_pinv=true)
+    b = randn(length(c)+1)
+    # Project out the component along the singular direction
+    E = eigen(Matrix(W); sortby=+)
+    bproj = b - (E.vectors[:,begin]'*b)*E.vectors[:,begin]
+    x = W \ bproj
+    @test norm(W*x - bproj) ≤ 1e-10
+end
+
 end # @testset "SymWoodbury"

test/woodbury.jl

@@ -202,4 +202,24 @@ W = Woodbury([randpsd(50) for _ in 1:4]...)
     @test all(logabsdet(W) .≈ logabsdet(Matrix(W)))
 end
 
+@testset "pinv" begin
+    # Build a matrix W that is equivalent to
+    #   A = [Diagonal(c) ones(length(c))]
+    #   A'A
+    # This is rank-deficient, so inv should fail but pinv should work
+    c = [-1.0, -2.0, -3.0, -4.0]
+    d = c.^2
+    push!(d, length(c))
+    U = [c zeros(length(c)); 0 1]
+    @test_throws SingularException Woodbury(Diagonal(d), U, [0 1; 1 0], U')
+    W = Woodbury(Diagonal(d), U, [0 1; 1 0], U'; use_pinv=true)
+    b = randn(length(c)+1)
+    # Project out the component along the singular direction
+    E = eigen(Matrix(W); sortby=+)
+    bproj = b - (E.vectors[:,begin]'*b)*E.vectors[:,begin]
+    x = W \ bproj
+    @test norm(W*x - bproj) ≤ 1e-10
+end
+
+
 end  # @testset "Woodbury"