From f4e84f20b94239df9e34b15c023446d206de319e Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 22 Aug 2026 17:32:34 -0500 Subject: [PATCH 1/7] Make SparseArrays a hard dependency The sparse support traversal and the sparse `linsolve` defaults live in `src/sparse_support.jl`, and the Unitful extension carries the sparse methods that resolve its ambiguity with them. Assisted-by: Claude Fable 5 Assisted-by: Claude Opus 5 --- Project.toml | 7 +--- ext/MatrixCoversSparseArraysUnitfulExt.jl | 35 ------------------ ext/MatrixCoversUnitfulExt.jl | 21 ++++++++++- src/MatrixCovers.jl | 2 + src/minimal_covers.jl | 16 ++++---- .../sparse_support.jl | 37 ++++++++----------- test/runtests.jl | 4 +- test/unitful.jl | 2 +- 8 files changed, 49 insertions(+), 75 deletions(-) delete mode 100644 ext/MatrixCoversSparseArraysUnitfulExt.jl rename ext/MatrixCoversSparseArraysExt.jl => src/sparse_support.jl (72%) diff --git a/Project.toml b/Project.toml index a528950..be9d392 100644 --- a/Project.toml +++ b/Project.toml @@ -7,19 +7,17 @@ authors = ["Tim Holy and contributors"] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [weakdeps] HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [extensions] MatrixCoversIpoptExt = ["Ipopt", "JuMP"] MatrixCoversJuMPExt = ["HiGHS", "JuMP"] -MatrixCoversSparseArraysExt = "SparseArrays" -MatrixCoversSparseArraysUnitfulExt = ["SparseArrays", "Unitful"] MatrixCoversUnitfulExt = "Unitful" [compat] @@ -46,11 +44,10 @@ HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [targets] -test = ["Aqua", "ExplicitImports", "HiGHS", "Ipopt", "JuMP", "OffsetArrays", "SparseArrays", "StableRNGs", "Statistics", "Test", "Unitful"] +test = ["Aqua", "ExplicitImports", "HiGHS", "Ipopt", "JuMP", "OffsetArrays", "StableRNGs", "Statistics", "Test", "Unitful"] diff --git a/ext/MatrixCoversSparseArraysUnitfulExt.jl b/ext/MatrixCoversSparseArraysUnitfulExt.jl deleted file mode 100644 index cda665e..0000000 --- a/ext/MatrixCoversSparseArraysUnitfulExt.jl +++ /dev/null @@ -1,35 +0,0 @@ -module MatrixCoversSparseArraysUnitfulExt - -# `MatrixCoversSparseArraysExt` types the matrix slot of its refiners, where `MatrixCoversUnitfulExt` types -# the element: neither is more specific for a sparse matrix of quantities, so the two -# are ambiguous there. These methods resolve that pair. They are the only overlap -- -# every other sparse method leaves its matrix slot untyped. - -using LinearAlgebra: LinearAlgebra, Hermitian, Symmetric -using MatrixCovers -using MatrixCovers: AbsLog -using SparseArrays: SparseArrays, SparseMatrixCSC -using Unitful: Quantity - -const MC = MatrixCovers - -# Sparse storage synthesizes structural zeros with `zero(eltype)`, so the element -# type is concrete and every entry carries the same unit. -const QSparse = SparseMatrixCSC{<:Quantity} -const QSparseSym = Union{QSparse, - Symmetric{<:Quantity,<:SparseMatrixCSC}, - Hermitian{<:Quantity,<:SparseMatrixCSC}} - -# `Unitful` triggers `MatrixCoversUnitfulExt` as well as this extension, so it is loaded whenever -# these methods can be called. Reaching for it here rather than at load time leaves -# the two extensions' load order free. -unitfulext() = Base.get_extension(MatrixCovers, :MatrixCoversUnitfulExt)::Module - -MC.symcover_min!(ϕ::AbsLog{2}, a::AbstractVector{<:Quantity}, A::QSparseSym; kwargs...) = - unitfulext().symstart!(MC.symcover_min!, a, A, ϕ; kwargs...) - -MC.cover_min!(ϕ::AbsLog{2}, a::AbstractVector{<:Quantity}, b::AbstractVector{<:Quantity}, - A::QSparse; kwargs...) = - unitfulext().asymstart!(MC.cover_min!, a, b, A, ϕ; kwargs...) - -end # module MatrixCoversSparseArraysUnitfulExt diff --git a/ext/MatrixCoversUnitfulExt.jl b/ext/MatrixCoversUnitfulExt.jl index 88e8877..f63402a 100644 --- a/ext/MatrixCoversUnitfulExt.jl +++ b/ext/MatrixCoversUnitfulExt.jl @@ -1,8 +1,9 @@ module MatrixCoversUnitfulExt -using LinearAlgebra: LinearAlgebra +using LinearAlgebra: LinearAlgebra, Hermitian, Symmetric using MatrixCovers using MatrixCovers: AbsLog, AbsLinear +using SparseArrays: SparseMatrixCSC using Unitful: Unitful, FreeUnits, Quantity, Unit, unit, ustrip const MC = MatrixCovers @@ -319,6 +320,24 @@ MC.soft_symcover_min!(a::QVector, A::QMatrix; kwargs...) = symstart!(MC.soft_sym MC.soft_cover_min(A::QMatrix; kwargs...) = asym(MC.soft_cover_min, A; kwargs...) MC.soft_cover_min!(a::QVector, b::QVector, A::QMatrix; kwargs...) = asymstart!(MC.soft_cover_min!, a, b, A; kwargs...) +# MatrixCovers types the matrix slot of its sparse refiners, where the methods above +# type the element: neither is more specific for a sparse matrix of quantities, so the +# two are ambiguous there. These methods resolve that pair. They are the only overlap -- +# every other sparse method leaves its matrix slot untyped. +# +# Sparse storage synthesizes structural zeros with `zero(eltype)`, so the element type +# is concrete and every entry carries the same unit. +const QSparse = SparseMatrixCSC{<:Quantity} +const QSparseSym = Union{QSparse, + Symmetric{<:Quantity,<:SparseMatrixCSC}, + Hermitian{<:Quantity,<:SparseMatrixCSC}} + +MC.symcover_min!(ϕ::AbsLog{2}, a::QVector, A::QSparseSym; kwargs...) = + symstart!(MC.symcover_min!, a, A, ϕ; kwargs...) + +MC.cover_min!(ϕ::AbsLog{2}, a::QVector, b::QVector, A::QSparse; kwargs...) = + asymstart!(MC.cover_min!, a, b, A, ϕ; kwargs...) + for P in PENALTIES @eval begin MC.soft_symcover(ϕ::$P, A::QMatrix; kwargs...) = sym(MC.soft_symcover, A, ϕ; kwargs...) diff --git a/src/MatrixCovers.jl b/src/MatrixCovers.jl index 7faa988..24fa38b 100644 --- a/src/MatrixCovers.jl +++ b/src/MatrixCovers.jl @@ -4,6 +4,7 @@ using LinearAlgebra: LinearAlgebra, Adjoint, Bidiagonal, Diagonal, Hermitian, SymTridiagonal, Symmetric, Transpose, Tridiagonal, dot, norm using PrecompileTools: PrecompileTools, @compile_workload using Random: Random, AbstractRNG, MersenneTwister +using SparseArrays: SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals export AbsLog, AbsLinear export cover_objective, iscover @@ -27,6 +28,7 @@ include("gram_covers.jl") # symmetric covers of A'*W*A from an asymmetric cov include("initializers.jl") # the start menu; consumed by both solver families below include("soft_covers.jl") include("minimal_covers.jl") +include("sparse_support.jl") # sparse traversal and the sparse solver defaults # True only when a MethodError's argument types are consistent with the calling diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 50cf6fa..069a7fc 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -22,9 +22,8 @@ Supported ϕ values: normal equations; `:lsqr` uses matrix-free LSQR (per-iteration cost O(nnz), intended for large sparse supports)). `linsolve` defaults to `:auto` for dense `A`; the `SparseMatrixCSC`/`Symmetric`/`Hermitian` sparse methods - (from the SparseArrays extension) default to `:lsqr` instead, since a dense - factorization of the reweighted normal equations is the wrong solve when - `nnz ≪ n²`. + default to `:lsqr` instead, since a dense factorization of the reweighted + normal equations is the wrong solve when `nnz ≪ n²`. - `AbsLog{1}()`: requires JuMP and HiGHS. - `AbsLinear{1}()`, `AbsLinear{2}()`: requires JuMP and Ipopt. These objectives are nonconvex. Each strategy in `strategies` is refined, and the best local @@ -61,9 +60,9 @@ Supported ϕ values: linear solve: `:auto`/`:dense` use a dense factorization of the reweighted normal equations; `:lsqr` uses matrix-free LSQR (per-iteration cost O(nnz), intended for large sparse supports)). `linsolve` defaults to `:auto` for - dense `A`; the `SparseMatrixCSC` sparse method (from the SparseArrays - extension) defaults to `:lsqr` instead, since a dense factorization of the - reweighted normal equations is the wrong solve when `nnz ≪ n²`. + dense `A`; the `SparseMatrixCSC` sparse method defaults to `:lsqr` instead, + since a dense factorization of the reweighted normal equations is the wrong + solve when `nnz ≪ n²`. - `AbsLog{1}()`: requires JuMP and HiGHS. - `AbsLinear{1}()`, `AbsLinear{2}()`: requires JuMP and Ipopt. These objectives are nonconvex. Each strategy in `strategies` is refined, and the best local @@ -370,9 +369,8 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), boost::Bool=true, fname=:symcover_min) linsolve in (:auto, :dense, :lsqr) || throw(ArgumentError("linsolve must be :auto, :dense, or :lsqr; got :$linsolve")) - # The shared entry to the native solve, reached from the sym `*_min` methods in - # this package and in the SparseArrays extension, so the precondition is checked - # once here rather than at each of them. + # The shared entry to the native solve, reached from every sym `*_min` method, + # so the precondition is checked once here rather than at each of them. require_abs_symmetric(A, fname) ax = axes(A, 1) axes(A, 2) == ax || throw(ArgumentError("symcover_min requires a square matrix")) diff --git a/ext/MatrixCoversSparseArraysExt.jl b/src/sparse_support.jl similarity index 72% rename from ext/MatrixCoversSparseArraysExt.jl rename to src/sparse_support.jl index 30c1452..bef8c67 100644 --- a/ext/MatrixCoversSparseArraysExt.jl +++ b/src/sparse_support.jl @@ -1,15 +1,10 @@ -module MatrixCoversSparseArraysExt - -using LinearAlgebra: LinearAlgebra, Hermitian, Symmetric -using SparseArrays: SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals -using MatrixCovers -using MatrixCovers: AbsLog, _symcover_min_abslog2, _cover_min_abslog2 +# Sparse-storage support traversal and the sparse defaults of the native solvers. # ============================================================ # Support traversal # ============================================================ -function MatrixCovers.foreach_support(f, A::SparseMatrixCSC) +function foreach_support(f, A::SparseMatrixCSC) rv, nzs = rowvals(A), nonzeros(A) for j in axes(A, 2) for k in nzrange(A, j) @@ -20,7 +15,7 @@ function MatrixCovers.foreach_support(f, A::SparseMatrixCSC) return nothing end -function MatrixCovers.foreach_support_sym(f, A::SparseMatrixCSC) +function foreach_support_sym(f, A::SparseMatrixCSC) ax = axes(A, 1) axes(A, 2) == ax || throw(DimensionMismatch("foreach_support_sym requires a square matrix, got axes $(axes(A))")) rv, nzs = rowvals(A), nonzeros(A) @@ -39,7 +34,7 @@ end # the stored (i, j) with i >= j is reported as (j, i). Complex `Hermitian` is # admitted alongside the real case because only `abs` of a stored value is ever # read, and `abs(A[i,j]) == abs(conj(A[j,i]))`. -function MatrixCovers.foreach_support_sym(f, +function foreach_support_sym(f, S::Union{Symmetric{<:Any,<:SparseMatrixCSC},Hermitian{<:Any,<:SparseMatrixCSC}}) P = parent(S) ax = axes(P, 1) @@ -69,9 +64,9 @@ end # orientations and the diagonal once; the magnitudes agree in both, including for a # complex `Hermitian`. Without this the wrappers fall back to the generic # `AbstractMatrix` method and its full-grid `getindex` scan. -function MatrixCovers.foreach_support(f, +function foreach_support(f, S::Union{Symmetric{<:Any,<:SparseMatrixCSC},Hermitian{<:Any,<:SparseMatrixCSC}}) - MatrixCovers.foreach_support_sym(S) do i, j, v + foreach_support_sym(S) do i, j, v f(i, j, v) i == j || f(j, i, v) end @@ -88,43 +83,41 @@ end # is the intended path when nnz ≪ n²; pass `linsolve=:auto`/`:dense` to force the # dense factorization. Only AbsLog{2} is native; other penalties dispatch to the # JuMP extension. -function MatrixCovers.symcover_min(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) +function symcover_min(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) a, _ = _symcover_min_abslog2(A; linsolve, kwargs...) return a end -function MatrixCovers.cover_min(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) +function cover_min(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) a, b, _ = _cover_min_abslog2(A; linsolve, kwargs...) return a, b end -function MatrixCovers.symcover_min(ϕ::AbsLog{2}, S::Symmetric{<:Any, <:SparseMatrixCSC}; linsolve::Symbol=:lsqr, kwargs...) +function symcover_min(ϕ::AbsLog{2}, S::Symmetric{<:Any, <:SparseMatrixCSC}; linsolve::Symbol=:lsqr, kwargs...) a, _ = _symcover_min_abslog2(S; linsolve, kwargs...) return a end -function MatrixCovers.symcover_min(ϕ::AbsLog{2}, H::Hermitian{<:Any, <:SparseMatrixCSC}; linsolve::Symbol=:lsqr, kwargs...) +function symcover_min(ϕ::AbsLog{2}, H::Hermitian{<:Any, <:SparseMatrixCSC}; linsolve::Symbol=:lsqr, kwargs...) a, _ = _symcover_min_abslog2(H; linsolve, kwargs...) return a end # The refiners take the same sparse `linsolve` default as the solvers above. -function MatrixCovers.symcover_min!(ϕ::AbsLog{2}, a::AbstractVector, +function symcover_min!(ϕ::AbsLog{2}, a::AbstractVector, S::Union{SparseMatrixCSC,Symmetric{<:Any,<:SparseMatrixCSC},Hermitian{<:Any,<:SparseMatrixCSC}}; linsolve::Symbol=:lsqr, kwargs...) - MatrixCovers._prepare_symcover_start!(a, S) + _prepare_symcover_start!(a, S) anew, _ = _symcover_min_abslog2(S; start=a, linsolve, kwargs...) a .= anew return a end -function MatrixCovers.cover_min!(ϕ::AbsLog{2}, a::AbstractVector, b::AbstractVector, - A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) - MatrixCovers._prepare_cover_start!(a, b, A) +function cover_min!(ϕ::AbsLog{2}, a::AbstractVector, b::AbstractVector, + A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) + _prepare_cover_start!(a, b, A) anew, bnew, _ = _cover_min_abslog2(A; start=(a, b), linsolve, kwargs...) a .= anew b .= bnew return a, b end - -end # module MatrixCoversSparseArraysExt diff --git a/test/runtests.jl b/test/runtests.jl index 00412e0..81efbad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,7 +2,7 @@ using MatrixCovers using MatrixCovers: foreach_support, foreach_support_sym, unconstrained_min!, tighten_cover! using MatrixCovers: ncomponents, support_components using JuMP, HiGHS, Ipopt # triggers MatrixCoversJuMPExt and MatrixCoversIpoptExt extensions -using SparseArrays # triggers MatrixCoversSparseArraysExt extension +using SparseArrays using Unitful # triggers MatrixCoversUnitfulExt extension using LinearAlgebra using OffsetArrays @@ -68,7 +68,7 @@ include("helpers.jl") # isbalanced, covaries, PENALTIES # Aqua checks the package alone; the extensions need their own sweep. @testset "method ambiguities" begin @test isempty(detect_ambiguities(MatrixCovers; recursive=true)) - for extname in (:MatrixCoversSparseArraysExt, :MatrixCoversJuMPExt, :MatrixCoversIpoptExt, :MatrixCoversUnitfulExt, :MatrixCoversSparseArraysUnitfulExt) + for extname in (:MatrixCoversJuMPExt, :MatrixCoversIpoptExt, :MatrixCoversUnitfulExt) ext = Base.get_extension(MatrixCovers, extname) @test ext !== nothing @test isempty(detect_ambiguities(MatrixCovers, ext; recursive=true)) diff --git a/test/unitful.jl b/test/unitful.jl index 38d160e..f5de585 100644 --- a/test/unitful.jl +++ b/test/unitful.jl @@ -100,7 +100,7 @@ @test all(unit.(v) == fill(u"m^-1", 3) for v in cover(S)) @test unit.(symcover(Symmetric(S))) == fill(u"m^-1", 3) - # The refiners are where MatrixCoversSparseArraysExt and MatrixCoversUnitfulExt overlap. + # The refiners are where the sparse methods and MatrixCoversUnitfulExt overlap. a = initialize_symcover(S) @test unit.(symcover_min!(AbsLog{2}(), a, S)) == fill(u"m^-1", 3) @test unit.(symcover_min!(AbsLog{2}(), a, Symmetric(S))) == fill(u"m^-1", 3) From 9f07263a0f73023756be4b3aeda79ce6eacf6d64 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 22 Aug 2026 17:42:37 -0500 Subject: [PATCH 2/7] Add a sparse Woodbury inner solve to the MMC MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Each AbsLog{2} Newton step solves reweighted normal equations whose matrix, on a near-complete support, is the complete-support matrix corrected by the zero set and by the currently violated entries. That split is `C + U*Uᵀ` with `C` sparse and positive definite and `U` of rank one (symmetric) or two (asymmetric), so a sparse Cholesky of `C` plus a Woodbury update replaces the dense factorization. `linsolve=:woodbury` selects it, `:auto` takes it wherever its requirements hold, and `stats.linsolve` reports which path ran. Assisted-by: Claude Fable 5 Assisted-by: Claude Opus 5 --- src/MatrixCovers.jl | 5 +- src/minimal_covers.jl | 327 +++++++++++++++++++++++++++++++++++++---- test/minimal_covers.jl | 103 ++++++++++++- test/soft_covers.jl | 14 ++ 4 files changed, 414 insertions(+), 35 deletions(-) diff --git a/src/MatrixCovers.jl b/src/MatrixCovers.jl index 24fa38b..2e24aa7 100644 --- a/src/MatrixCovers.jl +++ b/src/MatrixCovers.jl @@ -1,10 +1,11 @@ module MatrixCovers using LinearAlgebra: LinearAlgebra, Adjoint, Bidiagonal, Diagonal, Hermitian, - SymTridiagonal, Symmetric, Transpose, Tridiagonal, dot, norm + SymTridiagonal, Symmetric, Transpose, Tridiagonal, cholesky, + dot, norm using PrecompileTools: PrecompileTools, @compile_workload using Random: Random, AbstractRNG, MersenneTwister -using SparseArrays: SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals +using SparseArrays: SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals, sparse export AbsLog, AbsLinear export cover_objective, iscover diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 069a7fc..cdf088c 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -18,12 +18,19 @@ Supported ϕ values: - `AbsLog{2}()`: solved natively (no external solver). Accepts keyword arguments `κs` (the penalty-continuation schedule, default `(1e2, 1e4, 1e6, 1e8)`), `maxiter` (Newton steps per stage, default `40`), and `linsolve` (the inner - linear solve: `:auto`/`:dense` use a dense factorization of the reweighted - normal equations; `:lsqr` uses matrix-free LSQR (per-iteration cost O(nnz), - intended for large sparse supports)). `linsolve` defaults to `:auto` for - dense `A`; the `SparseMatrixCSC`/`Symmetric`/`Hermitian` sparse methods - default to `:lsqr` instead, since a dense factorization of the reweighted - normal equations is the wrong solve when `nnz ≪ n²`. + linear solve). `:dense` factorizes the reweighted normal equations densely, + at O(n³) per Newton step. `:woodbury` solves the same equations as a sparse + correction of the complete-support ones: the matrix is a sparse symmetric + positive-definite matrix plus `e*eᵀ`, so a sparse Cholesky and a + Sherman–Morrison update replace the dense factorization. It requires + `Float64` arithmetic and a support missing at most `n ÷ 4` entries in any + row, and raises an `ArgumentError` otherwise. `:lsqr` uses matrix-free LSQR + (per-iteration cost O(nnz), intended for large sparse supports). `:auto` + selects `:woodbury` where its requirements hold and `:dense` elsewhere. + `linsolve` defaults to `:auto` for dense `A`; the + `SparseMatrixCSC`/`Symmetric`/`Hermitian` sparse methods default to `:lsqr` + instead, since neither factorization of the reweighted normal equations is + the right solve when `nnz ≪ n²`. - `AbsLog{1}()`: requires JuMP and HiGHS. - `AbsLinear{1}()`, `AbsLinear{2}()`: requires JuMP and Ipopt. These objectives are nonconvex. Each strategy in `strategies` is refined, and the best local @@ -57,12 +64,19 @@ Supported ϕ values: - `AbsLog{2}()`: solved natively (no external solver). Accepts keyword arguments `κs` (the penalty-continuation schedule, default `(1e2, 1e4, 1e6, 1e8)`), `maxiter` (Newton steps per stage, default `40`), and `linsolve` (the inner - linear solve: `:auto`/`:dense` use a dense factorization of the reweighted - normal equations; `:lsqr` uses matrix-free LSQR (per-iteration cost O(nnz), - intended for large sparse supports)). `linsolve` defaults to `:auto` for - dense `A`; the `SparseMatrixCSC` sparse method defaults to `:lsqr` instead, - since a dense factorization of the reweighted normal equations is the wrong - solve when `nnz ≪ n²`. + linear solve). `:dense` factorizes the reweighted normal equations densely, + at O((m+n)³) per Newton step. `:woodbury` solves the same equations as a + sparse correction of the complete-support ones: the matrix is a sparse + symmetric positive-definite matrix plus a rank-two term, so a sparse Cholesky + and a Woodbury update replace the dense factorization. It requires `Float64` + arithmetic and a support missing at most `min(m, n) ÷ 4` entries in any row + or column, and raises an `ArgumentError` otherwise. `:lsqr` uses matrix-free + LSQR (per-iteration cost O(nnz), intended for large sparse supports). + `:auto` selects `:woodbury` where its requirements hold and `:dense` + elsewhere. `linsolve` defaults to `:auto` for dense `A`; the + `SparseMatrixCSC` sparse method defaults to `:lsqr` instead, since neither + factorization of the reweighted normal equations is the right solve when + `nnz ≪ n²`. - `AbsLog{1}()`: requires JuMP and HiGHS. - `AbsLinear{1}()`, `AbsLinear{2}()`: requires JuMP and Ipopt. These objectives are nonconvex. Each strategy in `strategies` is refined, and the best local @@ -285,12 +299,34 @@ function _prepare_cover_start!(a::AbstractVector, b::AbstractVector, A::Abstract end -# Inner linear solve for the AbsLog{2} MMC Newton steps. `:auto` (the default) -# forms and factorizes the reweighted normal equations densely, which is fastest -# for dense supports: an LAPACK Cholesky beats the matrix-free path because each -# LSQR iteration costs O(nnz) = O(n²) there. `:lsqr` forces the matrix-free path, -# whose per-iteration cost is O(nnz); it is the intended solve for large sparse -# supports (where nnz ≪ n²) and is used by the structured/sparse methods. +# Inner linear solve for the AbsLog{2} MMC Newton steps. `:dense` forms and +# factorizes the reweighted normal equations densely, at O(n³) per step. +# `:woodbury` splits the same matrix as `C + U Uᵀ`, where `C` is sparse (its +# off-diagonal pattern is the zero set of `A` together with the currently violated +# entries) and symmetric positive definite, and `U` has one column (symmetric) or +# two (asymmetric); a sparse Cholesky of `C` plus a Woodbury update then costs far +# less than the dense factorization whenever `A` is close to fully supported. +# `:auto` takes `:woodbury` where it applies and `:dense` otherwise. `:lsqr` forces +# the matrix-free path, whose per-iteration cost is O(nnz); it is the intended +# solve for large sparse supports (where nnz ≪ n²) and is used by the +# structured/sparse methods. +# +# `C` is positive definite because the complete-support matrix contributes `n` (or +# `m`) to each diagonal while the zero set subtracts a signless Laplacian `L_Z` +# with λmax(L_Z) ≤ 2·maxdeg(Z); requiring at most a quarter of a row to be zero +# keeps the difference bounded below by half the diagonal. CHOLMOD is the sparse +# factorization behind it, and it is reliable only in `Float64`, so that is the +# only working type the path accepts. + +# Append one COO triplet of the sparse Woodbury matrix `C`, tracking its diagonal +# in `diagacc` so the ridge can be sized without a second pass over `C`. +function _push_coo!(Ci, Cj, Cv, diagacc, p, q, v) + push!(Ci, p) + push!(Cj, q) + push!(Cv, v) + p == q && (diagacc[p] += v) + return nothing +end # Matrix-free LSQR (Paige & Saunders) for the weighted least-squares problem # `min ‖M x - b‖` underlying the reweighted normal equations `MᵀM x = Mᵀb`. @@ -360,15 +396,15 @@ end # Worker for `symcover_min(::AbsLog{2})`. Returns `(a, stats)` where `stats` is a # NamedTuple `(; nsolves, lsqriters, linsolve)` recording the number of inner linear # solves, the total LSQR iterations (0 on the dense path), and which path ran. -# `linsolve` is `:auto`/`:dense` (dense factorization) or `:lsqr` -# (matrix-free, for sparse supports). `start`, when given, is a positive cover of `A` +# `linsolve` reports the path that ran: `:dense`, `:woodbury`, or `:lsqr`. +# `start`, when given, is a positive cover of `A` # indexed like `axes(A, 1)` and supplies the first iterate in place of the cold # unweighted solve; the objective is convex, so it changes the path but not the result. function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), maxiter::Int=40, linsolve::Symbol=:auto, start=nothing, boost::Bool=true, fname=:symcover_min) - linsolve in (:auto, :dense, :lsqr) || - throw(ArgumentError("linsolve must be :auto, :dense, or :lsqr; got :$linsolve")) + linsolve in (:auto, :dense, :lsqr, :woodbury) || + throw(ArgumentError("linsolve must be :auto, :dense, :lsqr, or :woodbury; got :$linsolve")) # The shared entry to the native solve, reached from every sym `*_min` method, # so the precondition is checked once here rather than at each of them. require_abs_symmetric(A, fname) @@ -389,14 +425,46 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), edges = Tuple{Int,Int}[] cvals = T[] hassupp = falses(n) + maxzero = 0 # largest number of zeros in any row of `A` for (ip, i) in enumerate(ax) - for s in _slots(G, i) + slots = _slots(G, i) + for s in slots push!(edges, (ip, G.idx[s] - first(ax) + 1)) push!(cvals, log(G.val[s])) end - hassupp[ip] = !isempty(_slots(G, i)) + hassupp[ip] = !isempty(slots) + maxzero = max(maxzero, n - length(slots)) end ne = length(edges) + # The Woodbury path splits the normal equations around the complete-support + # matrix `n·I + e·eᵀ`, so its cost is set by the zero set `Z` rather than by `n`, + # and `n·I − L_Z` is positive definite only while `Z` stays thin. + use_woodbury = false + if !use_lsqr && linsolve !== :dense + ok = T === Float64 && maxzero <= n ÷ 4 + if linsolve === :woodbury && !ok + T === Float64 || + throw(ArgumentError("linsolve=:woodbury requires Float64 arithmetic, but `A` works in $T; use :dense or :lsqr")) + throw(ArgumentError("linsolve=:woodbury requires every row of `A` to have at most n ÷ 4 = $(n ÷ 4) zeros; got $maxzero")) + end + use_woodbury = ok + end + # Zero set of `A` in the same convention as `edges`: both orientations of an + # off-diagonal pair, the diagonal once. It is the off-diagonal pattern of the + # sparse `C` the Woodbury path factorizes. + zedges = Tuple{Int,Int}[] + if use_woodbury + mark = falses(n) + for (ip, i) in enumerate(ax) + for s in _slots(G, i) + mark[G.idx[s] - first(ax) + 1] = true + end + for jp in 1:n + mark[jp] || push!(zedges, (ip, jp)) + end + fill!(mark, false) + end + end fκ = function (α, κ) v = zero(T) for (e, (ip, jp)) in enumerate(edges) @@ -410,13 +478,23 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # equations are the signless Laplacian system `B α = f`. The dense path forms and # factorizes `B` (a support-free variable gets an identity row; a minimal # scale-relative ridge lifts the bipartite gauge null space, e.g. the `[0 1; 1 0]` - # support graph whose signless Laplacian is singular). The LSQR path applies `√W R` + # support graph whose signless Laplacian is singular). The Woodbury path solves the + # same regularized system exactly, splitting `B` as `C + e·eᵀ` around the + # complete-support matrix `n·I + e·eᵀ` and correcting `C` for the zero set and the + # violated entries. The LSQR path applies `√W R` # and its transpose matrix-free and warm-starts from the incoming iterate; it # solves the least-squares form directly, so its accuracy tracks the conditioning # of `√W R` (≈ √κ) rather than that of `B` (≈ κ). ws = zeros(T, ne) # √weight per support entry, frozen during one solve cv = zeros(T, ne) # √weight · log|A_ij| (LSQR right-hand side) f = zeros(T, n) + # COO triplets of `C`, refilled each Woodbury solve; `diagacc` accumulates its + # diagonal as they are appended. + Ci = Int[] + Cj = Int[] + Cv = T[] + diagacc = zeros(T, n) + rhs = zeros(T, n, 2) nsolves = Ref(0) nlsqr = Ref(0) solve_weighted = function (α, κ) @@ -447,6 +525,60 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), sol, it = _lsqr(Amul!, Atmul!, cv, α) nlsqr[] += it return sol + elseif use_woodbury + # `B = C + e·eᵀ` with `C = n·I − L_Z + (κ−1)·L_V`: the complete-support + # matrix, corrected by the zero set `Z` and by the currently violated + # entries `V`. `sparse` sums the duplicate triplets. + fill!(f, zero(T)) + fill!(diagacc, zero(T)) + empty!(Ci) + empty!(Cj) + empty!(Cv) + for p in 1:n + _push_coo!(Ci, Cj, Cv, diagacc, p, p, T(n)) + end + for (p, q) in zedges + _push_coo!(Ci, Cj, Cv, diagacc, p, p, -oneunit(T)) + _push_coo!(Ci, Cj, Cv, diagacc, p, q, -oneunit(T)) + end + for (e, (ip, jp)) in enumerate(edges) + c = cvals[e] + w = κ === nothing ? oneunit(T) : ((α[ip] + α[jp] - c) < 0 ? T(κ) : oneunit(T)) + f[ip] += w * c + if w != oneunit(T) + _push_coo!(Ci, Cj, Cv, diagacc, ip, ip, w - oneunit(T)) + _push_coo!(Ci, Cj, Cv, diagacc, ip, jp, w - oneunit(T)) + end + end + # Same ridge as the dense path, so both solve the same regularized system: + # `e·eᵀ` puts 1 on every diagonal of `B`, and every variable has support + # here, so no identity row arises. + dmax = zero(T) + for p in 1:n + dmax = max(dmax, diagacc[p] + oneunit(T)) + end + ridge = (dmax > 0 ? dmax : oneunit(T)) * eps(T) + for p in 1:n + push!(Ci, p) + push!(Cj, p) + push!(Cv, ridge) + end + F = cholesky(Symmetric(sparse(Ci, Cj, Cv, n, n))) + for p in 1:n + rhs[p, 1] = f[p] + rhs[p, 2] = oneunit(T) + end + # Sherman–Morrison: with y = C\f and u = C\e, (C + e·eᵀ)\f is + # y − u·(eᵀy)/(1 + eᵀu). + YU = F \ rhs + sy = zero(T) + su = zero(T) + for p in 1:n + sy += YU[p, 1] + su += YU[p, 2] + end + r = sy / (oneunit(T) + su) + return [YU[p, 1] - r * YU[p, 2] for p in 1:n] else fill!(f, zero(T)) B = zeros(T, n, n) @@ -501,7 +633,8 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), for (ip, i) in enumerate(ax) a[i] = hassupp[ip] ? exp(α[ip] + γ) : zero(T) end - return a, (; nsolves=nsolves[], lsqriters=nlsqr[], linsolve=(use_lsqr ? :lsqr : :dense)) + return a, (; nsolves=nsolves[], lsqriters=nlsqr[], + linsolve=(use_lsqr ? :lsqr : use_woodbury ? :woodbury : :dense)) end # Worker for `cover_min(::AbsLog{2})`. Returns `(a, b, stats)` with `stats` a @@ -511,8 +644,8 @@ end function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), maxiter::Int=40, linsolve::Symbol=:auto, start=nothing, boost::Bool=true) - linsolve in (:auto, :dense, :lsqr) || - throw(ArgumentError("linsolve must be :auto, :dense, or :lsqr; got :$linsolve")) + linsolve in (:auto, :dense, :lsqr, :woodbury) || + throw(ArgumentError("linsolve must be :auto, :dense, :lsqr, or :woodbury; got :$linsolve")) axr = axes(A, 1) axc = axes(A, 2) # The problem only ever depends on abs.(A), a real quantity, so the working type @@ -544,6 +677,43 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), ne = length(edges) hasrow = nzrow .> 0 hascol = nzcol .> 0 + # The Woodbury path splits the normal equations around the complete-support + # matrix `D + u_r·u_rᵀ + u_c·u_cᵀ`, so its cost is set by the zero set `Z` rather + # than by `N`, and `D − L_Z` is positive definite only while `Z` stays thin. The + # bound is taken against `min(m, n)`, the smaller of the two diagonal blocks. + maxzero = 0 + for ip in 1:m + maxzero = max(maxzero, n - nzrow[ip]) + end + for jp in 1:n + maxzero = max(maxzero, m - nzcol[jp]) + end + zbound = min(m, n) ÷ 4 + use_woodbury = false + if !use_lsqr && linsolve !== :dense + ok = T === Float64 && maxzero <= zbound + if linsolve === :woodbury && !ok + T === Float64 || + throw(ArgumentError("linsolve=:woodbury requires Float64 arithmetic, but `A` works in $T; use :dense or :lsqr")) + throw(ArgumentError("linsolve=:woodbury requires every row and column of `A` to have at most min(m, n) ÷ 4 = $zbound zeros; got $maxzero")) + end + use_woodbury = ok + end + # Zero set of `A` as (row position, column position) pairs missing from the + # support: the off-diagonal pattern of the sparse `C` the Woodbury path factorizes. + zedges = Tuple{Int,Int}[] + if use_woodbury + mark = falses(n) + for (ip, i) in enumerate(axr) + for s in _slots(G, i) + mark[G.idx[s] - first(axc) + 1] = true + end + for jp in 1:n + mark[jp] || push!(zedges, (ip, jp)) + end + fill!(mark, false) + end + end # Gauge vector: ±1 on supported variables, 0 on support-free ones (which carry # no constraint and are decoupled with an identity row in `solve_weighted`). v0 = zeros(T, N) @@ -563,10 +733,13 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end # Each Newton step solves the reweighted least-squares problem for the stacked # scales x = (α; β), residuals z_ij = α_i + β_j - log|A_ij|. Row and column scales - # share the global (e; −e) gauge; both paths pin it. The dense path adds the rank-1 + # share the global (e; −e) gauge; every path pins it. The dense path adds the rank-1 # term v0*v0ᵀ to the normal equations `B x = f` and factorizes (support-free # variables get an identity row; a support with more than one connected component - # carries additional per-component gauges, lifted by the ridge below). The LSQR + # carries additional per-component gauges, lifted by the ridge below). The Woodbury + # path solves the same regularized system exactly: `B + v0·v0ᵀ = C + U·Uᵀ` with `U` + # the row and column indicators, so a sparse Cholesky of `C` and a rank-two update + # replace the dense factorization. The LSQR # path appends one gauge row `v0ᵀ x = 0` to the least-squares system so `√W R` has # full column rank, applies it matrix-free, and warm-starts from the incoming # iterate. After the solve a closed-form shift, applied within each component, @@ -574,6 +747,13 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), f = zeros(T, N) ws = zeros(T, ne) # √weight per support entry (LSQR path) cv = zeros(T, ne + 1) # √weight · log|A_ij|, with a trailing 0 gauge target + # COO triplets of `C`, refilled each Woodbury solve; `diagacc` accumulates its + # diagonal as they are appended. + Ci = Int[] + Cj = Int[] + Cv = T[] + diagacc = zeros(T, N) + rhs = zeros(T, N, 3) nsolves = Ref(0) nlsqr = Ref(0) solve_weighted = function (x, κ) @@ -607,6 +787,90 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), sol, it = _lsqr(Amul!, Atmul!, cv, x) nlsqr[] += it return sol + elseif use_woodbury + # `B + v0·v0ᵀ = C + U·Uᵀ` with `C = D − L_Z + (κ−1)·L_V`, + # `D = diag(n·1_m, m·1_n)` and `U = [u_r u_c]` the row and column + # indicators: the complete-support matrix, corrected by the zero set `Z` + # and by the currently violated entries `V`. `sparse` sums the duplicate + # triplets. + fill!(f, zero(T)) + fill!(diagacc, zero(T)) + empty!(Ci) + empty!(Cj) + empty!(Cv) + for ip in 1:m + _push_coo!(Ci, Cj, Cv, diagacc, ip, ip, T(n)) + end + for jp in 1:n + _push_coo!(Ci, Cj, Cv, diagacc, m + jp, m + jp, T(m)) + end + for (ip, jp) in zedges + q = m + jp + _push_coo!(Ci, Cj, Cv, diagacc, ip, ip, -oneunit(T)) + _push_coo!(Ci, Cj, Cv, diagacc, q, q, -oneunit(T)) + _push_coo!(Ci, Cj, Cv, diagacc, ip, q, -oneunit(T)) + _push_coo!(Ci, Cj, Cv, diagacc, q, ip, -oneunit(T)) + end + for (e, (p, q)) in enumerate(edges) + c = cvals[e] + w = κ === nothing ? oneunit(T) : ((x[p] + x[q] - c) < 0 ? T(κ) : oneunit(T)) + f[p] += w * c + f[q] += w * c + if w != oneunit(T) + dw = w - oneunit(T) + _push_coo!(Ci, Cj, Cv, diagacc, p, p, dw) + _push_coo!(Ci, Cj, Cv, diagacc, q, q, dw) + _push_coo!(Ci, Cj, Cv, diagacc, p, q, dw) + _push_coo!(Ci, Cj, Cv, diagacc, q, p, dw) + end + end + # Same ridge as the dense path, so both solve the same regularized system: + # `U·Uᵀ` puts 1 on every diagonal of `B + v0·v0ᵀ`, and every variable has + # support here, so no identity row arises. + dmax = zero(T) + for p in 1:N + dmax = max(dmax, diagacc[p] + oneunit(T)) + end + ridge = (dmax > 0 ? dmax : oneunit(T)) * eps(T) + for p in 1:N + push!(Ci, p) + push!(Cj, p) + push!(Cv, ridge) + end + F = cholesky(Symmetric(sparse(Ci, Cj, Cv, N, N))) + fill!(rhs, zero(T)) + for p in 1:N + rhs[p, 1] = f[p] + end + for ip in 1:m + rhs[ip, 2] = oneunit(T) + end + for jp in 1:n + rhs[m+jp, 3] = oneunit(T) + end + # Woodbury with a 2x2 capacitance: with y = C\f and Y = C\U, + # (C + U·Uᵀ)\f is y − Y·((I₂ + UᵀY)\(Uᵀy)). + YU = F \ rhs + ty = zero(T) + tz = zero(T) + k11 = zero(T) + k12 = zero(T) + k21 = zero(T) + k22 = zero(T) + for ip in 1:m + ty += YU[ip, 1] + k11 += YU[ip, 2] + k12 += YU[ip, 3] + end + for jp in 1:n + q = m + jp + tz += YU[q, 1] + k21 += YU[q, 2] + k22 += YU[q, 3] + end + K = [oneunit(T)+k11 k12; k21 oneunit(T)+k22] + g = K \ T[ty, tz] + return [YU[p, 1] - g[1] * YU[p, 2] - g[2] * YU[p, 3] for p in 1:N] else fill!(f, zero(T)) B = v0 * v0' @@ -717,7 +981,8 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), for (jp, j) in enumerate(axc) b[j] = hascol[jp] ? exp(x[m+jp] - s[colcomp[jp]]) : zero(T) end - return a, b, (; nsolves=nsolves[], lsqriters=nlsqr[], linsolve=(use_lsqr ? :lsqr : :dense)) + return a, b, (; nsolves=nsolves[], lsqriters=nlsqr[], + linsolve=(use_lsqr ? :lsqr : use_woodbury ? :woodbury : :dense)) end # Workers for the soft (unconstrained) AbsLog{2} covers. The soft objective diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index 4e53aa2..ef5661f 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -104,8 +104,8 @@ end @testset "MMC native AbsLog{2} matrix-free LSQR path" begin # Invalid solver selection is rejected. - @test_throws "linsolve must be :auto, :dense, or :lsqr" symcover_min(AbsLog{2}(), [2.0 1.0; 1.0 3.0]; linsolve=:qr) - @test_throws "linsolve must be :auto, :dense, or :lsqr" cover_min(AbsLog{2}(), [2.0 1.0; 1.0 3.0]; linsolve=:qr) + @test_throws "linsolve must be :auto, :dense, :lsqr, or :woodbury" symcover_min(AbsLog{2}(), [2.0 1.0; 1.0 3.0]; linsolve=:qr) + @test_throws "linsolve must be :auto, :dense, :lsqr, or :woodbury" cover_min(AbsLog{2}(), [2.0 1.0; 1.0 3.0]; linsolve=:qr) # The matrix-free LSQR path reproduces the dense path and the HiGHS reference # across the committed symmetric library, and returns a feasible cover. @@ -143,6 +143,105 @@ end @test a[2] * b[1] ≈ 1.0 end +# The Woodbury path solves the same regularized normal equations as the dense path, +# through a sparse Cholesky of `C` and a low-rank update, so the two must agree to +# roundoff. The tolerance is loose relative to `eps`: a converged cover pins the +# objective far more tightly than its own entries, so the two solves separate at +# roughly the square root of the working precision. +@testset "MMC native AbsLog{2} Woodbury path" begin + rng = StableRNG(9) + lognormal(m, n) = exp.(randn(rng, m, n)) + symlognormal(n) = (X = lognormal(n, n); (X .+ X') ./ 2) + + @testset "symmetric, n = $n" for n in (6, 30, 120) + A = symlognormal(n) + variants = ["all nonzero" => A, + "zero diagonal" => A - Diagonal(A)] + # A sparse symmetric zero set, thin enough to stay inside the n ÷ 4 guard. + Z = symlognormal(n) + mask = rand(rng, n, n) .< 0.5 / max(n ÷ 8, 1) + mask = mask .| mask' + Z[mask] .= 0.0 + maximum(count(iszero, Z; dims=2)) <= n ÷ 4 && push!(variants, "sparse zeros" => Z) + for (name, M) in variants + @testset "$name" begin + ad, sd = MatrixCovers._symcover_min_abslog2(M; linsolve=:dense) + aw, sw = MatrixCovers._symcover_min_abslog2(M; linsolve=:woodbury) + aa, sa = MatrixCovers._symcover_min_abslog2(M) + @test sd.linsolve === :dense + @test sw.linsolve === :woodbury + @test sa.linsolve === :woodbury + @test aw ≈ ad rtol=1e-7 + @test aa == aw + @test iscover(aw, M; atol=1e-8) + end + end + end + + @testset "asymmetric, ($m, $n)" for (m, n) in ((8, 6), (30, 22), (120, 90)) + A = lognormal(m, n) + # The gauge direction (e; −e) leaves every product a[i]*b[j] fixed, so only the + # products are required to agree between the two solves. + k = min(m, n) ÷ 4 # the largest zero band the guard admits + for (name, M) in ("all nonzero" => A, "one zero band" => (B = copy(A); B[1:k, 1] .= 0.0; B)) + @testset "$name" begin + ad, bd, sd = MatrixCovers._cover_min_abslog2(M; linsolve=:dense) + aw, bw, sw = MatrixCovers._cover_min_abslog2(M; linsolve=:woodbury) + _, _, sa = MatrixCovers._cover_min_abslog2(M) + @test sd.linsolve === :dense + @test sw.linsolve === :woodbury + @test sa.linsolve === :woodbury + @test aw .* bw' ≈ ad .* bd' rtol=1e-7 + @test iscover(aw, bw, M; atol=1e-7) + end + end + end + + # Only abs.(A) is read, so a complex Hermitian takes the same path and lands on + # the same cover as its magnitude matrix. + n = 8 + M = randn(rng, ComplexF64, n, n) + H = Hermitian(M + M') + aw, sw = MatrixCovers._symcover_min_abslog2(H; linsolve=:woodbury) + @test sw.linsolve === :woodbury + @test aw ≈ MatrixCovers._symcover_min_abslog2(abs.(Matrix(H)); linsolve=:dense)[1] rtol=1e-7 + + # Offset axes and views index the support through `axes(A)`, not `1:n`, on this + # path as on the others. + A = symlognormal(12) + aref = symcover_min(AbsLog{2}(), A; linsolve=:woodbury) + Ao = OffsetArray(A, -3, -3) + ao = symcover_min(AbsLog{2}(), Ao; linsolve=:woodbury) + @test axes(ao, 1) == axes(Ao, 1) + @test collect(ao) ≈ aref rtol=1e-10 + Av = view(symlognormal(16), 3:14, 3:14) + @test symcover_min(AbsLog{2}(), Matrix(Av); linsolve=:woodbury) ≈ + symcover_min(AbsLog{2}(), Av; linsolve=:woodbury) rtol=1e-10 + Ag = lognormal(14, 10) + Agv = view(Ag, 2:13, 2:9) + av, bv = cover_min(AbsLog{2}(), Agv; linsolve=:woodbury) + am, bm = cover_min(AbsLog{2}(), Matrix(Agv); linsolve=:woodbury) + @test av .* bv' ≈ am .* bm' rtol=1e-10 + + # The zero set must stay inside the guard, and the arithmetic must be Float64. + holey = symlognormal(12) + holey[1, 1:5] .= 0.0 + holey[1:5, 1] .= 0.0 + @test_throws "at most n ÷ 4 = 3 zeros; got 5" symcover_min(AbsLog{2}(), holey; linsolve=:woodbury) + @test MatrixCovers._symcover_min_abslog2(holey)[2].linsolve === :dense + gholey = lognormal(12, 12) + gholey[1, 1:5] .= 0.0 + @test_throws "at most min(m, n) ÷ 4 = 3 zeros; got 5" cover_min(AbsLog{2}(), gholey; linsolve=:woodbury) + @test MatrixCovers._cover_min_abslog2(gholey)[3].linsolve === :dense + + A32 = Float32.(symlognormal(8)) + @test_throws "requires Float64 arithmetic" symcover_min(AbsLog{2}(), A32; linsolve=:woodbury) + @test MatrixCovers._symcover_min_abslog2(A32)[2].linsolve === :dense + G32 = Float32.(lognormal(8, 6)) + @test_throws "requires Float64 arithmetic" cover_min(AbsLog{2}(), G32; linsolve=:woodbury) + @test MatrixCovers._cover_min_abslog2(G32)[3].linsolve === :dense +end + @testset "MMC disconnected-support gauge" begin # A support graph that splits into k connected components carries k independent # (e; −e) gauges. The asymmetric dense normal equations pin only the global one diff --git a/test/soft_covers.jl b/test/soft_covers.jl index 6130888..bf4839b 100644 --- a/test/soft_covers.jl +++ b/test/soft_covers.jl @@ -439,6 +439,20 @@ end # The `:lsqr` and dense paths solve the same problem. @test soft_symcover_min(AbsLog{2}(), sym_zeros; linsolve=:lsqr) ≈ soft_symcover_min(AbsLog{2}(), sym_zeros; linsolve=:dense) rtol=1e-6 + + # The soft cover is the hard worker's cold solve, so it reaches the Woodbury path + # too: one sparse Cholesky in place of the dense factorization, same answer. + rng = StableRNG(77) + X = exp.(randn(rng, 40, 40)) + Adense = (X .+ X') ./ 2 + @test soft_symcover_min(AbsLog{2}(), Adense; linsolve=:woodbury) ≈ + soft_symcover_min(AbsLog{2}(), Adense; linsolve=:dense) rtol=1e-8 + @test MatrixCovers._soft_symcover_min_abslog2(Adense)[2].linsolve === :woodbury + Y = exp.(randn(rng, 40, 30)) + aw, bw = soft_cover_min(AbsLog{2}(), Y; linsolve=:woodbury) + ad, bd = soft_cover_min(AbsLog{2}(), Y; linsolve=:dense) + @test aw .* bw' ≈ ad .* bd' rtol=1e-8 + @test MatrixCovers._soft_cover_min_abslog2(Y)[3].linsolve === :woodbury end @testset "soft_symcover!/soft_cover! refiners" begin From 7895e4f5a0aac3fd70a95a704215410fc01acb65 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 22 Aug 2026 17:46:03 -0500 Subject: [PATCH 3/7] End a penalty stage on a sign-stable full step The dense and Woodbury Newton steps minimize the quadratic model of the current stage exactly. When such a step is taken whole and the violated set at the new iterate is the one the step froze its weights on, the gradient of the convex stage objective there is the model's, which vanishes: the iterate is the stage minimizer and the stage ends without a confirmation solve. The inexact `:lsqr` steps keep the decrease test as their sole criterion. Assisted-by: Claude Fable 5 Assisted-by: Claude Opus 5 --- src/minimal_covers.jl | 46 ++++++++++++++++++++++++++++++++++++------ test/minimal_covers.jl | 29 ++++++++++++++++++++++++++ 2 files changed, 69 insertions(+), 6 deletions(-) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index cdf088c..14293e9 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -318,6 +318,16 @@ end # factorization behind it, and it is reliable only in `Float64`, so that is the # only working type the path accepts. +# True when the residuals of `x` are violated on exactly the entries `pat` marks. +# `edges` and `cvals` are the support list and its `log|A_ij|`, as gathered by the +# workers below. +function _violated_matches(pat, edges, cvals, x) + for (e, (p, q)) in enumerate(edges) + ((x[p] + x[q] - cvals[e]) < zero(eltype(cvals))) == pat[e] || return false + end + return true +end + # Append one COO triplet of the sparse Woodbury matrix `C`, tracking its diagonal # in `diagacc` so the ridge can be sized without a second pass over `C`. function _push_coo!(Ci, Cj, Cv, diagacc, p, q, v) @@ -488,6 +498,9 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), ws = zeros(T, ne) # √weight per support entry, frozen during one solve cv = zeros(T, ne) # √weight · log|A_ij| (LSQR right-hand side) f = zeros(T, n) + # Entries the frozen weights of the current solve treat as violated. A full Newton + # step that leaves this pattern intact has landed on the stage's minimizer. + vpat = falses(ne) # COO triplets of `C`, refilled each Woodbury solve; `diagacc` accumulates its # diagonal as they are appended. Ci = Int[] @@ -543,9 +556,11 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end for (e, (ip, jp)) in enumerate(edges) c = cvals[e] - w = κ === nothing ? oneunit(T) : ((α[ip] + α[jp] - c) < 0 ? T(κ) : oneunit(T)) + viol = κ !== nothing && (α[ip] + α[jp] - c) < 0 + vpat[e] = viol + w = viol ? T(κ) : oneunit(T) f[ip] += w * c - if w != oneunit(T) + if viol _push_coo!(Ci, Cj, Cv, diagacc, ip, ip, w - oneunit(T)) _push_coo!(Ci, Cj, Cv, diagacc, ip, jp, w - oneunit(T)) end @@ -584,7 +599,9 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), B = zeros(T, n, n) for (e, (ip, jp)) in enumerate(edges) c = cvals[e] - w = κ === nothing ? oneunit(T) : ((α[ip] + α[jp] - c) < 0 ? T(κ) : oneunit(T)) + viol = κ !== nothing && (α[ip] + α[jp] - c) < 0 + vpat[e] = viol + w = viol ? T(κ) : oneunit(T) f[ip] += w * c B[ip, ip] += w B[ip, jp] += w @@ -615,6 +632,11 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), fnew = fκ(α .+ t .* (αnew .- α), κ) end α = α .+ t .* (αnew .- α) + # `f_κ` is convex and the dense and Woodbury steps solve its quadratic model + # exactly, so a whole step that leaves the violated set unchanged has reached + # the stage's minimizer: the gradient there is the model's, which is zero. + # The `:lsqr` solves are inexact and carry no such guarantee. + !use_lsqr && isone(t) && _violated_matches(vpat, edges, cvals, α) && break fcur - fnew <= 5000 * eps(T) * max(fcur, one(T)) && break fcur = fnew end @@ -747,6 +769,9 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), f = zeros(T, N) ws = zeros(T, ne) # √weight per support entry (LSQR path) cv = zeros(T, ne + 1) # √weight · log|A_ij|, with a trailing 0 gauge target + # Entries the frozen weights of the current solve treat as violated. A full Newton + # step that leaves this pattern intact has landed on the stage's minimizer. + vpat = falses(ne) # COO triplets of `C`, refilled each Woodbury solve; `diagacc` accumulates its # diagonal as they are appended. Ci = Int[] @@ -813,10 +838,12 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end for (e, (p, q)) in enumerate(edges) c = cvals[e] - w = κ === nothing ? oneunit(T) : ((x[p] + x[q] - c) < 0 ? T(κ) : oneunit(T)) + viol = κ !== nothing && (x[p] + x[q] - c) < 0 + vpat[e] = viol + w = viol ? T(κ) : oneunit(T) f[p] += w * c f[q] += w * c - if w != oneunit(T) + if viol dw = w - oneunit(T) _push_coo!(Ci, Cj, Cv, diagacc, p, p, dw) _push_coo!(Ci, Cj, Cv, diagacc, q, q, dw) @@ -876,7 +903,9 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), B = v0 * v0' for (e, (p, q)) in enumerate(edges) c = cvals[e] - w = κ === nothing ? oneunit(T) : ((x[p] + x[q] - c) < 0 ? T(κ) : oneunit(T)) + viol = κ !== nothing && (x[p] + x[q] - c) < 0 + vpat[e] = viol + w = viol ? T(κ) : oneunit(T) f[p] += w * c f[q] += w * c B[p, p] += w @@ -931,6 +960,11 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), fnew = fκ(x .+ t .* (xnew .- x), κ) end x = x .+ t .* (xnew .- x) + # `f_κ` is convex and the dense and Woodbury steps solve its quadratic model + # exactly, so a whole step that leaves the violated set unchanged has reached + # the stage's minimizer: the gradient there is the model's, which is zero. + # The `:lsqr` solves are inexact and carry no such guarantee. + !use_lsqr && isone(t) && _violated_matches(vpat, edges, cvals, x) && break fcur - fnew <= 5000 * eps(T) * max(fcur, one(T)) && break fcur = fnew end diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index ef5661f..49fa0e4 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -242,6 +242,35 @@ end @test MatrixCovers._cover_min_abslog2(G32)[3].linsolve === :dense end +# A Newton step is exact on the dense and Woodbury paths, so a whole step that leaves +# the violated set unchanged has already reached the minimizer of the current +# penalty stage, and the stage ends without a confirmation solve. The `:lsqr` steps +# are inexact and keep the decrease test as their sole criterion, which is what makes +# their solve counts the reference here. +@testset "MMC exact paths stop on a sign-stable Newton step" begin + rng = StableRNG(31) + A = (X = exp.(randn(rng, 60, 60)); (X .+ X') ./ 2) + ad, sd = MatrixCovers._symcover_min_abslog2(A; linsolve=:dense) + aw, sw = MatrixCovers._symcover_min_abslog2(A; linsolve=:woodbury) + al, sl = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr) + @test ad ≈ al rtol=1e-6 + @test aw ≈ al rtol=1e-6 + # One solve per κ stage is saved; `κs` has four stages by default. + @test sd.nsolves == sw.nsolves + @test sd.nsolves <= sl.nsolves - length((1e2, 1e4, 1e6, 1e8)) + @test sd.nsolves <= 26 + + G = exp.(randn(rng, 60, 45)) + gd, hd, td = MatrixCovers._cover_min_abslog2(G; linsolve=:dense) + gw, hw, tw = MatrixCovers._cover_min_abslog2(G; linsolve=:woodbury) + gl, hl, tl = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr) + @test gd .* hd' ≈ gl .* hl' rtol=1e-6 + @test gw .* hw' ≈ gl .* hl' rtol=1e-6 + @test td.nsolves == tw.nsolves + @test td.nsolves <= tl.nsolves - length((1e2, 1e4, 1e6, 1e8)) + @test td.nsolves <= 24 +end + @testset "MMC disconnected-support gauge" begin # A support graph that splits into k connected components carries k independent # (e; −e) gauges. The asymmetric dense normal equations pin only the global one From bcf844185d12590b32d7e996f4bbcdf567fb3d3d Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 22 Aug 2026 18:05:31 -0500 Subject: [PATCH 4/7] Solve well-conditioned Woodbury stages by CG MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `C + U*Uᵀ` applies at O(n + |Z| + |V|) without being formed, and Gershgorin bounds its condition number by `1 + (κ−1)*2*maxdeg(V)/n`. Where that estimate stays under 1000, Jacobi-preconditioned conjugate gradients reach the same answer to rounding in a few hundred such applications; the sparse Cholesky serves the ill-conditioned stages and any run that exhausts the iteration cap. `stats.cgiters` counts the iterations. The objective sweep now reports whether the violated set still matches the one the step froze its weights on, so the line search and the stage's stopping test share one pass over the support. Assisted-by: Claude Fable 5 Assisted-by: Claude Opus 5 --- src/minimal_covers.jl | 374 ++++++++++++++++++++++++++++++++--------- test/minimal_covers.jl | 18 ++ 2 files changed, 312 insertions(+), 80 deletions(-) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 14293e9..20a4995 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -306,6 +306,14 @@ end # entries) and symmetric positive definite, and `U` has one column (symmetric) or # two (asymmetric); a sparse Cholesky of `C` plus a Woodbury update then costs far # less than the dense factorization whenever `A` is close to fully supported. +# It has two exact sub-paths. `C + U Uᵀ` can also be applied without being formed, +# at O(n + |Z| + |V|) per application, and Gershgorin bounds its condition number by +# `1 + (κ−1)·2·maxdeg(V)/n` against the complete-support diagonal; while that +# estimate stays under 1000, Jacobi-preconditioned conjugate gradients converge to +# rounding in a few hundred such applications, which beats a factorization whose +# fill on the near-random violated pattern of the early stages approaches dense. +# Above it, the factorization runs. Either way the answer is exact to rounding, +# which is what the sign-stability stopping test in the continuation loop requires. # `:auto` takes `:woodbury` where it applies and `:dense` otherwise. `:lsqr` forces # the matrix-free path, whose per-iteration cost is O(nnz); it is the intended # solve for large sparse supports (where nnz ≪ n²) and is used by the @@ -318,24 +326,90 @@ end # factorization behind it, and it is reliable only in `Float64`, so that is the # only working type the path accepts. -# True when the residuals of `x` are violated on exactly the entries `pat` marks. -# `edges` and `cvals` are the support list and its `log|A_ij|`, as gathered by the -# workers below. -function _violated_matches(pat, edges, cvals, x) - for (e, (p, q)) in enumerate(edges) - ((x[p] + x[q] - cvals[e]) < zero(eltype(cvals))) == pat[e] || return false +# `B*x` for the symmetric Woodbury system +# `B = dbase·I + e·eᵀ − L_Z + dκ·L_V` (`dbase = n + ridge`), applied without forming +# a matrix: `L_X·x` accumulates `x[p] + x[q]` into `y[p]` over the directed edges of +# `X`, so one application costs O(n + |Z| + |V|) rather than O(nnz). +function _symwoodbury_mul!(y, x, dbase, zedges, vedges, dκ) + s = zero(eltype(y)) + for p in eachindex(x) + s += x[p] end - return true + @. y = dbase * x + s + for (p, q) in zedges + y[p] -= x[p] + x[q] + end + for (p, q) in vedges + y[p] += dκ * (x[p] + x[q]) + end + return y end -# Append one COO triplet of the sparse Woodbury matrix `C`, tracking its diagonal -# in `diagacc` so the ridge can be sized without a second pass over `C`. -function _push_coo!(Ci, Cj, Cv, diagacc, p, q, v) - push!(Ci, p) - push!(Cj, q) - push!(Cv, v) - p == q && (diagacc[p] += v) - return nothing +# `B*x` for the asymmetric Woodbury system +# `B = D + ridge·I + u_r·u_rᵀ + u_c·u_cᵀ − L_Z + dκ·L_V` on the stacked positions +# `1:m` (rows) and `m+1:m+n` (columns), with `drow`/`dcol` the two diagonal blocks of +# `D + ridge·I`. Every edge of `L_X` adds `x[p] + x[q]` at both of its ends. +function _woodbury_mul!(y, x, m, drow, dcol, zedges, vedges, dκ) + sr = zero(eltype(y)) + for p in 1:m + sr += x[p] + end + sc = zero(eltype(y)) + for p in (m+1):length(x) + sc += x[p] + end + for p in 1:m + y[p] = drow * x[p] + sr + end + for p in (m+1):length(x) + y[p] = dcol * x[p] + sc + end + for (p, q) in zedges + t = x[p] + x[q] + y[p] -= t + y[q] -= t + end + for (p, q) in vedges + t = dκ * (x[p] + x[q]) + y[p] += t + y[q] += t + end + return y +end + +# Jacobi-preconditioned conjugate gradients for the symmetric positive-definite +# Woodbury system `B x = f`, with `Bmul!(y, x)` applying `B` and `dg` holding its +# diagonal. `x` carries the warm start in and the iterate out; `r`, `z`, `d`, `Ad` +# are work vectors of the same length. Returns `(iters, converged)`. +# +# The Newton step has to be exact to rounding for the stage's sign-stability +# stopping test to mean what it says, so `tol` sits at the level of `eps` and a run +# that exhausts `maxiter` reports failure instead of a partial answer; the caller +# then falls back to the factorization, which is exact. +function _pcg!(Bmul!, x, dg, f, r, z, d, Ad, maxiter::Int, tol) + Bmul!(r, x) + @. r = f - r + nrm = norm(r) + @. z = r / dg + copyto!(d, z) + rz = dot(r, z) + iters = 0 + for k in 1:maxiter + nrm <= tol && break + iters = k + Bmul!(Ad, d) + dAd = dot(d, Ad) + dAd > 0 || break + a = rz / dAd + @. x += a * d + @. r -= a * Ad + nrm = norm(r) + @. z = r / dg + rznew = dot(r, z) + @. d = z + (rznew / rz) * d + rz = rznew + end + return iters, nrm <= tol end # Matrix-free LSQR (Paige & Saunders) for the weighted least-squares problem @@ -404,8 +478,9 @@ function _lsqr(Amul!, Atmul!, b::AbstractVector{T}, x0::AbstractVector{T}; end # Worker for `symcover_min(::AbsLog{2})`. Returns `(a, stats)` where `stats` is a -# NamedTuple `(; nsolves, lsqriters, linsolve)` recording the number of inner linear -# solves, the total LSQR iterations (0 on the dense path), and which path ran. +# NamedTuple `(; nsolves, lsqriters, cgiters, linsolve)` recording the number of inner +# linear solves, the total LSQR and conjugate-gradient iterations (0 on paths that run +# neither), and which path ran. # `linsolve` reports the path that ran: `:dense`, `:woodbury`, or `:lsqr`. # `start`, when given, is a positive cover of `A` # indexed like `axes(A, 1)` and supplies the first iterate in place of the cold @@ -463,6 +538,10 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # off-diagonal pair, the diagonal once. It is the off-diagonal pattern of the # sparse `C` the Woodbury path factorizes. zedges = Tuple{Int,Int}[] + # Diagonal of `C` before any entry is violated: `n` from the complete-support + # matrix, less what `L_Z` puts there — one per zero entry of the row, and one more + # for a zero on the diagonal, which `L_Z` counts twice. + czero = fill(T(n), use_woodbury ? n : 0) if use_woodbury mark = falses(n) for (ip, i) in enumerate(ax) @@ -470,7 +549,10 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), mark[G.idx[s] - first(ax) + 1] = true end for jp in 1:n - mark[jp] || push!(zedges, (ip, jp)) + mark[jp] && continue + push!(zedges, (ip, jp)) + czero[ip] -= oneunit(T) + ip == jp && (czero[ip] -= oneunit(T)) end fill!(mark, false) end @@ -483,6 +565,20 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end return v end + # The objective and the violated set at `α` from one sweep: the line search needs + # the value and the stage's stopping test needs to know whether the set still + # matches `pat`, and both read the same residuals. + fκpat = function (α, κ, pat) + v = zero(T) + same = true + for (e, (ip, jp)) in enumerate(edges) + z = α[ip] + α[jp] - cvals[e] + viol = z < 0 + v += (viol ? T(κ) : oneunit(T)) * z^2 + same &= viol == pat[e] + end + return v, same + end # Each Newton step freezes the weights at the current α and solves the reweighted # least-squares problem `min ‖√W (Rα - c)‖`, `(Rα)_e = α_i + α_j`, whose normal # equations are the signless Laplacian system `B α = f`. The dense path forms and @@ -501,15 +597,22 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # Entries the frozen weights of the current solve treat as violated. A full Newton # step that leaves this pattern intact has landed on the stage's minimizer. vpat = falses(ne) - # COO triplets of `C`, refilled each Woodbury solve; `diagacc` accumulates its - # diagonal as they are appended. + vedges = Tuple{Int,Int}[] # the violated entries of the current solve + degV = zeros(Int, n) # violated entries per row + dg = zeros(T, n) # diagonal of `B`, for the ridge and the CG preconditioner + # COO triplets of `C`, refilled whenever a Woodbury solve is factorized. Ci = Int[] Cj = Int[] Cv = T[] - diagacc = zeros(T, n) rhs = zeros(T, n, 2) + cgx = zeros(T, n) + cgr = zeros(T, n) + cgz = zeros(T, n) + cgd = zeros(T, n) + cgAd = zeros(T, n) nsolves = Ref(0) nlsqr = Ref(0) + ncg = Ref(0) solve_weighted = function (α, κ) nsolves[] += 1 if use_lsqr @@ -541,42 +644,78 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), elseif use_woodbury # `B = C + e·eᵀ` with `C = n·I − L_Z + (κ−1)·L_V`: the complete-support # matrix, corrected by the zero set `Z` and by the currently violated - # entries `V`. `sparse` sums the duplicate triplets. + # entries `V`. One O(nnz) sweep collects the right-hand side, the violated + # set, and the diagonal of `B`; everything after it is O(n + |Z| + |V|). + dκ = κ === nothing ? zero(T) : T(κ) - oneunit(T) fill!(f, zero(T)) - fill!(diagacc, zero(T)) - empty!(Ci) - empty!(Cj) - empty!(Cv) - for p in 1:n - _push_coo!(Ci, Cj, Cv, diagacc, p, p, T(n)) - end - for (p, q) in zedges - _push_coo!(Ci, Cj, Cv, diagacc, p, p, -oneunit(T)) - _push_coo!(Ci, Cj, Cv, diagacc, p, q, -oneunit(T)) - end + copyto!(dg, czero) + fill!(degV, 0) + empty!(vedges) for (e, (ip, jp)) in enumerate(edges) c = cvals[e] viol = κ !== nothing && (α[ip] + α[jp] - c) < 0 vpat[e] = viol - w = viol ? T(κ) : oneunit(T) - f[ip] += w * c + f[ip] += (viol ? T(κ) : oneunit(T)) * c if viol - _push_coo!(Ci, Cj, Cv, diagacc, ip, ip, w - oneunit(T)) - _push_coo!(Ci, Cj, Cv, diagacc, ip, jp, w - oneunit(T)) + push!(vedges, (ip, jp)) + degV[ip] += 1 + dg[ip] += dκ + ip == jp && (dg[ip] += dκ) end end # Same ridge as the dense path, so both solve the same regularized system: # `e·eᵀ` puts 1 on every diagonal of `B`, and every variable has support # here, so no identity row arises. dmax = zero(T) + maxdegV = 0 for p in 1:n - dmax = max(dmax, diagacc[p] + oneunit(T)) + dmax = max(dmax, dg[p] + oneunit(T)) + maxdegV = max(maxdegV, degV[p]) end ridge = (dmax > 0 ? dmax : oneunit(T)) * eps(T) + for p in 1:n + dg[p] += oneunit(T) + ridge + end + # Gershgorin on `(κ−1)·L_V` against a diagonal of at least `n` bounds the + # condition number of `B`. While that bound is small the structured matvec + # plus conjugate gradients reaches the same answer in a few hundred O(n + + # |Z| + |V|) iterations, which is far cheaper than a factorization whose + # fill, on the near-random violated pattern of the early stages, is close + # to dense. + κest = oneunit(T) + dκ * 2 * maxdegV / n + if κest <= 1000 + copyto!(cgx, α) + dbase = T(n) + ridge + Bmul! = (y, x) -> _symwoodbury_mul!(y, x, dbase, zedges, vedges, dκ) + it, ok = _pcg!(Bmul!, cgx, dg, f, cgr, cgz, cgd, cgAd, + 50 + 20 * ceil(Int, sqrt(κest)), 100 * eps(T) * norm(f)) + ncg[] += it + ok && return copy(cgx) + end + # `sparse` sums the duplicate triplets; the ridge rides on `C`'s diagonal. + empty!(Ci) + empty!(Cj) + empty!(Cv) for p in 1:n push!(Ci, p) push!(Cj, p) - push!(Cv, ridge) + push!(Cv, T(n) + ridge) + end + for (p, q) in zedges + push!(Ci, p) + push!(Cj, p) + push!(Cv, -oneunit(T)) + push!(Ci, p) + push!(Cj, q) + push!(Cv, -oneunit(T)) + end + for (p, q) in vedges + push!(Ci, p) + push!(Cj, p) + push!(Cv, dκ) + push!(Ci, p) + push!(Cj, q) + push!(Cv, dκ) end F = cholesky(Symmetric(sparse(Ci, Cj, Cv, n, n))) for p in 1:n @@ -626,17 +765,20 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), for _ in 1:maxiter αnew = solve_weighted(α, κ) t = one(T) - fnew = fκ(αnew, κ) + αt = αnew + fnew, stable = fκpat(αt, κ, vpat) while fnew > fcur && t > 500_000 * eps(T) t /= 2 - fnew = fκ(α .+ t .* (αnew .- α), κ) + αt = α .+ t .* (αnew .- α) + fnew = fκ(αt, κ) + stable = false end - α = α .+ t .* (αnew .- α) + α = αt # `f_κ` is convex and the dense and Woodbury steps solve its quadratic model # exactly, so a whole step that leaves the violated set unchanged has reached # the stage's minimizer: the gradient there is the model's, which is zero. # The `:lsqr` solves are inexact and carry no such guarantee. - !use_lsqr && isone(t) && _violated_matches(vpat, edges, cvals, α) && break + !use_lsqr && stable && break fcur - fnew <= 5000 * eps(T) * max(fcur, one(T)) && break fcur = fnew end @@ -655,12 +797,12 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), for (ip, i) in enumerate(ax) a[i] = hassupp[ip] ? exp(α[ip] + γ) : zero(T) end - return a, (; nsolves=nsolves[], lsqriters=nlsqr[], + return a, (; nsolves=nsolves[], lsqriters=nlsqr[], cgiters=ncg[], linsolve=(use_lsqr ? :lsqr : use_woodbury ? :woodbury : :dense)) end # Worker for `cover_min(::AbsLog{2})`. Returns `(a, b, stats)` with `stats` a -# NamedTuple `(; nsolves, lsqriters, linsolve)` (see `_symcover_min_abslog2`). +# NamedTuple `(; nsolves, lsqriters, cgiters, linsolve)` (see `_symcover_min_abslog2`). # `start`, when given, is a positive cover `(a, b)` indexed like the rows and columns # of `A`, supplying the first iterate in place of the cold unweighted solve. function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), @@ -721,17 +863,30 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end use_woodbury = ok end - # Zero set of `A` as (row position, column position) pairs missing from the - # support: the off-diagonal pattern of the sparse `C` the Woodbury path factorizes. + # Zero set of `A` as stacked-position pairs missing from the support: the + # off-diagonal pattern of the sparse `C` the Woodbury path factorizes. `czero` is + # the diagonal of `C` before any entry is violated — the complete-support value + # `n` on rows and `m` on columns, less one per zero entry at each of its ends. zedges = Tuple{Int,Int}[] + czero = zeros(T, use_woodbury ? N : 0) if use_woodbury + for ip in 1:m + czero[ip] = T(n) + end + for jp in 1:n + czero[m+jp] = T(m) + end mark = falses(n) for (ip, i) in enumerate(axr) for s in _slots(G, i) mark[G.idx[s] - first(axc) + 1] = true end for jp in 1:n - mark[jp] || push!(zedges, (ip, jp)) + mark[jp] && continue + q = m + jp + push!(zedges, (ip, q)) + czero[ip] -= oneunit(T) + czero[q] -= oneunit(T) end fill!(mark, false) end @@ -753,6 +908,20 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end return v end + # The objective and the violated set at `x` from one sweep: the line search needs + # the value and the stage's stopping test needs to know whether the set still + # matches `pat`, and both read the same residuals. + fκpat = function (x, κ, pat) + v = zero(T) + same = true + for (e, (p, q)) in enumerate(edges) + z = x[p] + x[q] - cvals[e] + viol = z < 0 + v += (viol ? T(κ) : oneunit(T)) * z^2 + same &= viol == pat[e] + end + return v, same + end # Each Newton step solves the reweighted least-squares problem for the stacked # scales x = (α; β), residuals z_ij = α_i + β_j - log|A_ij|. Row and column scales # share the global (e; −e) gauge; every path pins it. The dense path adds the rank-1 @@ -772,15 +941,22 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # Entries the frozen weights of the current solve treat as violated. A full Newton # step that leaves this pattern intact has landed on the stage's minimizer. vpat = falses(ne) - # COO triplets of `C`, refilled each Woodbury solve; `diagacc` accumulates its - # diagonal as they are appended. + vedges = Tuple{Int,Int}[] # the violated entries of the current solve + degV = zeros(Int, N) # violated entries per row and per column + dg = zeros(T, N) # diagonal of `B`, for the ridge and the CG preconditioner + # COO triplets of `C`, refilled whenever a Woodbury solve is factorized. Ci = Int[] Cj = Int[] Cv = T[] - diagacc = zeros(T, N) rhs = zeros(T, N, 3) + cgx = zeros(T, N) + cgr = zeros(T, N) + cgz = zeros(T, N) + cgd = zeros(T, N) + cgAd = zeros(T, N) nsolves = Ref(0) nlsqr = Ref(0) + ncg = Ref(0) solve_weighted = function (x, κ) nsolves[] += 1 if use_lsqr @@ -818,24 +994,11 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # indicators: the complete-support matrix, corrected by the zero set `Z` # and by the currently violated entries `V`. `sparse` sums the duplicate # triplets. + dκ = κ === nothing ? zero(T) : T(κ) - oneunit(T) fill!(f, zero(T)) - fill!(diagacc, zero(T)) - empty!(Ci) - empty!(Cj) - empty!(Cv) - for ip in 1:m - _push_coo!(Ci, Cj, Cv, diagacc, ip, ip, T(n)) - end - for jp in 1:n - _push_coo!(Ci, Cj, Cv, diagacc, m + jp, m + jp, T(m)) - end - for (ip, jp) in zedges - q = m + jp - _push_coo!(Ci, Cj, Cv, diagacc, ip, ip, -oneunit(T)) - _push_coo!(Ci, Cj, Cv, diagacc, q, q, -oneunit(T)) - _push_coo!(Ci, Cj, Cv, diagacc, ip, q, -oneunit(T)) - _push_coo!(Ci, Cj, Cv, diagacc, q, ip, -oneunit(T)) - end + copyto!(dg, czero) + fill!(degV, 0) + empty!(vedges) for (e, (p, q)) in enumerate(edges) c = cvals[e] viol = κ !== nothing && (x[p] + x[q] - c) < 0 @@ -844,25 +1007,73 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), f[p] += w * c f[q] += w * c if viol - dw = w - oneunit(T) - _push_coo!(Ci, Cj, Cv, diagacc, p, p, dw) - _push_coo!(Ci, Cj, Cv, diagacc, q, q, dw) - _push_coo!(Ci, Cj, Cv, diagacc, p, q, dw) - _push_coo!(Ci, Cj, Cv, diagacc, q, p, dw) + push!(vedges, (p, q)) + degV[p] += 1 + degV[q] += 1 + dg[p] += dκ + dg[q] += dκ end end # Same ridge as the dense path, so both solve the same regularized system: # `U·Uᵀ` puts 1 on every diagonal of `B + v0·v0ᵀ`, and every variable has # support here, so no identity row arises. dmax = zero(T) + maxdegV = 0 for p in 1:N - dmax = max(dmax, diagacc[p] + oneunit(T)) + dmax = max(dmax, dg[p] + oneunit(T)) + maxdegV = max(maxdegV, degV[p]) end ridge = (dmax > 0 ? dmax : oneunit(T)) * eps(T) for p in 1:N + dg[p] += oneunit(T) + ridge + end + # Gershgorin on `(κ−1)·L_V` against the smaller diagonal block bounds the + # condition number of `B`. While that bound is small the structured matvec + # plus conjugate gradients reaches the same answer in a few hundred + # O(N + |Z| + |V|) iterations, which is far cheaper than a factorization + # whose fill on the near-random violated pattern of the early stages + # approaches dense. + κest = oneunit(T) + dκ * 2 * maxdegV / min(m, n) + if κest <= 1000 + copyto!(cgx, x) + drow = T(n) + ridge + dcol = T(m) + ridge + Bmul! = (yy, xx) -> _woodbury_mul!(yy, xx, m, drow, dcol, zedges, vedges, dκ) + it, ok = _pcg!(Bmul!, cgx, dg, f, cgr, cgz, cgd, cgAd, + 50 + 20 * ceil(Int, sqrt(κest)), 100 * eps(T) * norm(f)) + ncg[] += it + ok && return copy(cgx) + end + # `sparse` sums the duplicate triplets; the ridge rides on `C`'s diagonal. + empty!(Ci) + empty!(Cj) + empty!(Cv) + for p in 1:N + push!(Ci, p) + push!(Cj, p) + push!(Cv, czero[p] + ridge) + end + for (p, q) in zedges + push!(Ci, p) + push!(Cj, q) + push!(Cv, -oneunit(T)) + push!(Ci, q) + push!(Cj, p) + push!(Cv, -oneunit(T)) + end + for (p, q) in vedges + push!(Ci, p) + push!(Cj, p) + push!(Cv, dκ) + push!(Ci, q) + push!(Cj, q) + push!(Cv, dκ) push!(Ci, p) + push!(Cj, q) + push!(Cv, dκ) + push!(Ci, q) push!(Cj, p) - push!(Cv, ridge) + push!(Cv, dκ) end F = cholesky(Symmetric(sparse(Ci, Cj, Cv, N, N))) fill!(rhs, zero(T)) @@ -954,17 +1165,20 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), for _ in 1:maxiter xnew = solve_weighted(x, κ) t = one(T) - fnew = fκ(xnew, κ) + xt = xnew + fnew, stable = fκpat(xt, κ, vpat) while fnew > fcur && t > 500_000 * eps(T) t /= 2 - fnew = fκ(x .+ t .* (xnew .- x), κ) + xt = x .+ t .* (xnew .- x) + fnew = fκ(xt, κ) + stable = false end - x = x .+ t .* (xnew .- x) + x = xt # `f_κ` is convex and the dense and Woodbury steps solve its quadratic model # exactly, so a whole step that leaves the violated set unchanged has reached # the stage's minimizer: the gradient there is the model's, which is zero. # The `:lsqr` solves are inexact and carry no such guarantee. - !use_lsqr && isone(t) && _violated_matches(vpat, edges, cvals, x) && break + !use_lsqr && stable && break fcur - fnew <= 5000 * eps(T) * max(fcur, one(T)) && break fcur = fnew end @@ -1015,7 +1229,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), for (jp, j) in enumerate(axc) b[j] = hascol[jp] ? exp(x[m+jp] - s[colcomp[jp]]) : zero(T) end - return a, b, (; nsolves=nsolves[], lsqriters=nlsqr[], + return a, b, (; nsolves=nsolves[], lsqriters=nlsqr[], cgiters=ncg[], linsolve=(use_lsqr ? :lsqr : use_woodbury ? :woodbury : :dense)) end diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index 49fa0e4..4f57150 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -174,6 +174,11 @@ end @test aw ≈ ad rtol=1e-7 @test aa == aw @test iscover(aw, M; atol=1e-8) + # Both Woodbury sub-paths run within a continuation: the early stages + # are well enough conditioned for conjugate gradients, the late ones + # are not, and both are exact. + @test sw.cgiters > 0 + @test sd.cgiters == 0 end end end @@ -193,6 +198,8 @@ end @test sa.linsolve === :woodbury @test aw .* bw' ≈ ad .* bd' rtol=1e-7 @test iscover(aw, bw, M; atol=1e-7) + @test sw.cgiters > 0 + @test sd.cgiters == 0 end end end @@ -223,6 +230,17 @@ end am, bm = cover_min(AbsLog{2}(), Matrix(Agv); linsolve=:woodbury) @test av .* bv' ≈ am .* bm' rtol=1e-10 + # A single well-conditioned stage exercises the conjugate-gradient sub-path alone. + A1 = symlognormal(24) + c1, s1 = MatrixCovers._symcover_min_abslog2(A1; κs=(1e2,), linsolve=:woodbury) + @test s1.cgiters > 0 + @test c1 ≈ MatrixCovers._symcover_min_abslog2(A1; κs=(1e2,), linsolve=:dense)[1] rtol=1e-8 + G1 = lognormal(24, 18) + p1, q1, t1 = MatrixCovers._cover_min_abslog2(G1; κs=(1e2,), linsolve=:woodbury) + pd1, qd1, _ = MatrixCovers._cover_min_abslog2(G1; κs=(1e2,), linsolve=:dense) + @test t1.cgiters > 0 + @test p1 .* q1' ≈ pd1 .* qd1' rtol=1e-8 + # The zero set must stay inside the guard, and the arithmetic must be Float64. holey = symlognormal(12) holey[1, 1:5] .= 0.0 From 5681591fc817ea120df7a0f9bc6e08578e6bb9b2 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 22 Aug 2026 19:08:17 -0500 Subject: [PATCH 5/7] Right-precondition the LSQR inner solve MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The penalty strength enters the least-squares matrix only through the rows it reweights, so `M = diag(RᵀR) + (κ−1)*Σ_{e∈V} rₑ*rₑᵀ` carries all of it. Every generalized eigenvalue of `(RᵀWR, M)` is a mediant of eigenvalues of `(RᵀR, diag(RᵀR))` and so lies in their range, which holds the iteration count fixed as the continuation raises κ. LSQR runs on `√W*R*K⁻ᵀ`, `K` the permuted sparse Cholesky factor of `M`, applied through CHOLMOD's factor components; the violated rows join `M` only once diagonal scaling alone would leave the system ill conditioned. CHOLMOD is Float64-only, so other working types keep the plain matrix-free iteration. Assisted-by: Claude Fable 5 Assisted-by: Claude Opus 5 --- src/MatrixCovers.jl | 2 +- src/minimal_covers.jl | 236 +++++++++++++++++++++++++++++++++++++++-- test/minimal_covers.jl | 26 +++++ 3 files changed, 256 insertions(+), 8 deletions(-) diff --git a/src/MatrixCovers.jl b/src/MatrixCovers.jl index 2e24aa7..056a566 100644 --- a/src/MatrixCovers.jl +++ b/src/MatrixCovers.jl @@ -2,7 +2,7 @@ module MatrixCovers using LinearAlgebra: LinearAlgebra, Adjoint, Bidiagonal, Diagonal, Hermitian, SymTridiagonal, Symmetric, Transpose, Tridiagonal, cholesky, - dot, norm + dot, mul!, norm using PrecompileTools: PrecompileTools, @compile_workload using Random: Random, AbstractRNG, MersenneTwister using SparseArrays: SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals, sparse diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 20a4995..fc018e7 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -25,7 +25,10 @@ Supported ϕ values: Sherman–Morrison update replace the dense factorization. It requires `Float64` arithmetic and a support missing at most `n ÷ 4` entries in any row, and raises an `ArgumentError` otherwise. `:lsqr` uses matrix-free LSQR - (per-iteration cost O(nnz), intended for large sparse supports). `:auto` + (per-iteration cost O(nnz), intended for large sparse supports), right + preconditioned in `Float64` by a sparse Cholesky of the diagonal of the + unweighted normal matrix plus the rows the penalty currently weights, which + keeps its iteration count from growing with the penalty strength. `:auto` selects `:woodbury` where its requirements hold and `:dense` elsewhere. `linsolve` defaults to `:auto` for dense `A`; the `SparseMatrixCSC`/`Symmetric`/`Hermitian` sparse methods default to `:lsqr` @@ -71,7 +74,10 @@ Supported ϕ values: and a Woodbury update replace the dense factorization. It requires `Float64` arithmetic and a support missing at most `min(m, n) ÷ 4` entries in any row or column, and raises an `ArgumentError` otherwise. `:lsqr` uses matrix-free - LSQR (per-iteration cost O(nnz), intended for large sparse supports). + LSQR (per-iteration cost O(nnz), intended for large sparse supports), right + preconditioned in `Float64` by a sparse Cholesky of the diagonal of the + unweighted normal matrix plus the rows the penalty currently weights, which + keeps its iteration count from growing with the penalty strength. `:auto` selects `:woodbury` where its requirements hold and `:dense` elsewhere. `linsolve` defaults to `:auto` for dense `A`; the `SparseMatrixCSC` sparse method defaults to `:lsqr` instead, since neither @@ -317,7 +323,8 @@ end # `:auto` takes `:woodbury` where it applies and `:dense` otherwise. `:lsqr` forces # the matrix-free path, whose per-iteration cost is O(nnz); it is the intended # solve for large sparse supports (where nnz ≪ n²) and is used by the -# structured/sparse methods. +# structured/sparse methods. It is right preconditioned, which is what keeps its +# iteration count from growing as the continuation raises κ. # # `C` is positive definite because the complete-support matrix contributes `n` (or # `m`) to each diagonal while the zero set subtracts a signless Laplacian `L_Z` @@ -377,6 +384,19 @@ function _woodbury_mul!(y, x, m, drow, dcol, zedges, vedges, dκ) return y end +# Right preconditioner for the LSQR path. `M = diag(RᵀR) + (κ−1)·Σ_{e∈V} rₑ·rₑᵀ` +# holds the diagonal of the unweighted normal matrix together with the exact +# contribution of the rows LSQR weights by κ. All of the κ-dependence of `RᵀWR` sits +# in those rows, and every generalized eigenvalue of `(RᵀWR, M)` is a mediant of +# eigenvalues of `(RᵀR, diag(RᵀR))` and so lies in their range, which is what keeps +# the preconditioned iteration count from growing with κ. +# +# `M = K·Kᵀ` with `K` the permuted Cholesky factor `PtL` of `M`. LSQR then runs on +# `√W·R·K⁻ᵀ` in the variable `y = Kᵀ·x`, applying `K⁻¹` and `K⁻ᵀ` through the CHOLMOD +# factor components `F.PtL` and `F.UP`. `Kᵀ` itself is never needed as a product: the +# warm start is `Kᵀ·x₀ = K⁻¹·(M·x₀)`, which the same two objects and one sparse +# matrix-vector product supply. + # Jacobi-preconditioned conjugate gradients for the symmetric positive-definite # Woodbury system `B x = f`, with `Bmul!(y, x)` applying `B` and `dg` holding its # diagonal. `x` carries the warm start in and the iterate out; `r`, `z`, `d`, `Ad` @@ -420,6 +440,10 @@ end # strength κ) rather than that of `MᵀM` (≈ κ); at κ = 1e8 the squared conditioning # breaks CG while LSQR stays accurate. # +# The callers pass a right-preconditioned operator, so `x` here is the preconditioned +# variable and `M` is `√W·R·K⁻ᵀ`; the preconditioner is described above `_lsqr`'s +# callers in each worker. +# # The penalty least-squares problem is inconsistent (its optimal residual is # nonzero), so the stopping test is on the normal-equations residual # ‖Mᵀ(b - Mx)‖ ≤ atol · ‖M‖ · ‖b - Mx‖, both estimated from the bidiagonalization @@ -501,6 +525,9 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), T = float(real(eltype(A))) n = length(ax) use_lsqr = linsolve === :lsqr + # CHOLMOD, which factors the LSQR preconditioner, is reliable only in Float64; + # other working types run the plain matrix-free iteration. + use_precond = use_lsqr && T === Float64 # Support entries, one per residual z_ij = α_i + α_j - log|A_ij|, with `cvals` # holding log|A_ij| alongside. The gather reports each off-diagonal pair in both # orientations and the diagonal once, which is the full-grid weighting the @@ -600,6 +627,30 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), vedges = Tuple{Int,Int}[] # the violated entries of the current solve degV = zeros(Int, n) # violated entries per row dg = zeros(T, n) # diagonal of `B`, for the ridge and the CG preconditioner + # Diagonal of the unweighted normal matrix `RᵀR`, the base of the LSQR + # preconditioner: each directed support entry puts 1 at each of its ends, and a + # diagonal entry, whose row of `R` is `2·e_p`, puts 4. A support-free variable is + # given 1 so the preconditioner stays positive definite. + dpart = zeros(T, use_lsqr ? n : 0) + if use_lsqr + for (ip, jp) in edges + if ip == jp + dpart[ip] += 4 * oneunit(T) + else + dpart[ip] += oneunit(T) + dpart[jp] += oneunit(T) + end + end + for p in 1:n + dpart[p] > 0 || (dpart[p] = oneunit(T)) + end + end + mdiag = zeros(T, use_lsqr ? n : 0) # the violated rows' diagonal, per unit of κ−1 + Mi = Int[] # COO triplets of the preconditioner + Mj = Int[] + Mv = T[] + px = zeros(T, use_lsqr ? n : 0) # scale vector recovered from the LSQR variable + pg = zeros(T, use_lsqr ? n : 0) # `Rᵀ√W y` before the preconditioner is applied # COO triplets of `C`, refilled whenever a Woodbury solve is factorized. Ci = Int[] Cj = Int[] @@ -616,12 +667,90 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), solve_weighted = function (α, κ) nsolves[] += 1 if use_lsqr + dκ = κ === nothing ? zero(T) : T(κ) - oneunit(T) + empty!(vedges) + fill!(mdiag, zero(T)) for (e, (ip, jp)) in enumerate(edges) c = cvals[e] - w = κ === nothing ? oneunit(T) : ((α[ip] + α[jp] - c) < 0 ? T(κ) : oneunit(T)) - sw = sqrt(w) + viol = κ !== nothing && (α[ip] + α[jp] - c) < 0 + vpat[e] = viol + sw = sqrt(viol ? T(κ) : oneunit(T)) ws[e] = sw cv[e] = sw * c + if viol && use_precond + push!(vedges, (ip, jp)) + if ip == jp + mdiag[ip] += 4 * oneunit(T) + else + mdiag[ip] += oneunit(T) + mdiag[jp] += oneunit(T) + end + end + end + if use_precond + # Diagonal scaling alone leaves a conditioning that grows with κ once + # the violated rows dominate a variable's diagonal; past that point + # they enter the preconditioner in full, and its Cholesky pays for + # itself in the iterations it removes. + κest = oneunit(T) + for p in 1:n + κest = max(κest, oneunit(T) + dκ * 2 * mdiag[p] / dpart[p]) + end + empty!(Mi) + empty!(Mj) + empty!(Mv) + for p in 1:n + push!(Mi, p) + push!(Mj, p) + push!(Mv, dpart[p]) + end + if κest > 1000 + for (p, q) in vedges + if p == q + push!(Mi, p) + push!(Mj, p) + push!(Mv, 4 * dκ) + else + push!(Mi, p) + push!(Mj, p) + push!(Mv, dκ) + push!(Mi, q) + push!(Mj, q) + push!(Mv, dκ) + push!(Mi, p) + push!(Mj, q) + push!(Mv, dκ) + push!(Mi, q) + push!(Mj, p) + push!(Mv, dκ) + end + end + end + Msp = sparse(Mi, Mj, Mv, n, n) + MF = cholesky(Symmetric(Msp)) + Kc = MF.PtL + Uc = MF.UP + Pmul! = function (y, yv) + xv = Uc \ yv + for (e, (ip, jp)) in enumerate(edges) + y[e] = ws[e] * (xv[ip] + xv[jp]) + end + return y + end + Ptmul! = function (z, y) + fill!(pg, zero(T)) + for (e, (ip, jp)) in enumerate(edges) + t = ws[e] * y[e] + pg[ip] += t + pg[jp] += t + end + copyto!(z, Kc \ pg) + return z + end + mul!(px, Msp, α) + soly, it = _lsqr(Pmul!, Ptmul!, cv, Kc \ px) + nlsqr[] += it + return Uc \ soly end Amul! = function (y, x) for (e, (ip, jp)) in enumerate(edges) @@ -820,6 +949,9 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), n = length(axc) N = m + n use_lsqr = linsolve === :lsqr + # CHOLMOD, which factors the LSQR preconditioner, is reliable only in Float64; + # other working types run the plain matrix-free iteration. + use_precond = use_lsqr && T === Float64 # Support entries as edges linking a row position ip to a column position m+jp, # with `cvals` holding log|A_ij| alongside. Internal positions 1:m index rows, # m+1:m+n index columns, and results are scattered back through axr/axc so A's @@ -944,6 +1076,26 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), vedges = Tuple{Int,Int}[] # the violated entries of the current solve degV = zeros(Int, N) # violated entries per row and per column dg = zeros(T, N) # diagonal of `B`, for the ridge and the CG preconditioner + # Diagonal of the unweighted normal matrix of the gauge-augmented system, + # `RᵀR + v0·v0ᵀ`: the support degree at each position, plus the gauge row's 1. A + # support-free variable takes that 1 alone, which keeps the preconditioner + # positive definite. + dpart = zeros(T, use_lsqr ? N : 0) + if use_lsqr + for (p, q) in edges + dpart[p] += oneunit(T) + dpart[q] += oneunit(T) + end + for p in 1:N + dpart[p] += oneunit(T) + end + end + mdiag = zeros(T, use_lsqr ? N : 0) # the violated rows' diagonal, per unit of κ−1 + Mi = Int[] # COO triplets of the preconditioner + Mj = Int[] + Mv = T[] + px = zeros(T, use_lsqr ? N : 0) # scale vector recovered from the LSQR variable + pg = zeros(T, use_lsqr ? N : 0) # `Rᵀ√W y` before the preconditioner is applied # COO triplets of `C`, refilled whenever a Woodbury solve is factorized. Ci = Int[] Cj = Int[] @@ -960,14 +1112,84 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), solve_weighted = function (x, κ) nsolves[] += 1 if use_lsqr + dκ = κ === nothing ? zero(T) : T(κ) - oneunit(T) + empty!(vedges) + fill!(mdiag, zero(T)) for (e, (p, q)) in enumerate(edges) c = cvals[e] - w = κ === nothing ? oneunit(T) : ((x[p] + x[q] - c) < 0 ? T(κ) : oneunit(T)) - sw = sqrt(w) + viol = κ !== nothing && (x[p] + x[q] - c) < 0 + vpat[e] = viol + sw = sqrt(viol ? T(κ) : oneunit(T)) ws[e] = sw cv[e] = sw * c + if viol && use_precond + push!(vedges, (p, q)) + mdiag[p] += oneunit(T) + mdiag[q] += oneunit(T) + end end g = ne + 1 # index of the appended gauge row + if use_precond + # Diagonal scaling alone leaves a conditioning that grows with κ once + # the violated rows dominate a variable's diagonal; past that point + # they enter the preconditioner in full, and its Cholesky pays for + # itself in the iterations it removes. + κest = oneunit(T) + for p in 1:N + κest = max(κest, oneunit(T) + dκ * 2 * mdiag[p] / dpart[p]) + end + empty!(Mi) + empty!(Mj) + empty!(Mv) + for p in 1:N + push!(Mi, p) + push!(Mj, p) + push!(Mv, dpart[p]) + end + if κest > 1000 + for (p, q) in vedges + push!(Mi, p) + push!(Mj, p) + push!(Mv, dκ) + push!(Mi, q) + push!(Mj, q) + push!(Mv, dκ) + push!(Mi, p) + push!(Mj, q) + push!(Mv, dκ) + push!(Mi, q) + push!(Mj, p) + push!(Mv, dκ) + end + end + Msp = sparse(Mi, Mj, Mv, N, N) + MF = cholesky(Symmetric(Msp)) + Kc = MF.PtL + Uc = MF.UP + Pmul! = function (y, yv) + xv = Uc \ yv + for (e, (p, q)) in enumerate(edges) + y[e] = ws[e] * (xv[p] + xv[q]) + end + y[g] = dot(v0, xv) + return y + end + Ptmul! = function (z, y) + fill!(pg, zero(T)) + for (e, (p, q)) in enumerate(edges) + t = ws[e] * y[e] + pg[p] += t + pg[q] += t + end + @. pg += v0 * y[g] + copyto!(z, Kc \ pg) + return z + end + mul!(px, Msp, x) + soly, it = _lsqr(Pmul!, Ptmul!, cv, Kc \ px) + nlsqr[] += it + return Uc \ soly + end Amul! = function (y, xx) for (e, (p, q)) in enumerate(edges) y[e] = ws[e] * (xx[p] + xx[q]) diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index 4f57150..a397ecd 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -289,6 +289,32 @@ end @test td.nsolves <= 24 end +# The LSQR preconditioner absorbs the rows the continuation weights by κ, so the +# generalized spectrum it iterates on is the unweighted one and the iteration count +# stops growing as κ rises. CHOLMOD factors it, so it applies only in Float64; +# narrower and wider types run the plain matrix-free iteration. +@testset "MMC :lsqr iteration count is bounded across the continuation" begin + rng = StableRNG(5) + A = (X = exp.(randn(rng, 120, 120)); (X .+ X') ./ 2) + ad, _ = MatrixCovers._symcover_min_abslog2(A; linsolve=:dense) + al, sl = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr) + @test al ≈ ad rtol=1e-6 + @test sl.lsqriters <= 60 * sl.nsolves + + G = exp.(randn(rng, 120, 90)) + gd, hd, _ = MatrixCovers._cover_min_abslog2(G; linsolve=:dense) + gl, hl, tl = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr) + @test gl .* hl' ≈ gd .* hd' rtol=1e-6 + @test tl.lsqriters <= 60 * tl.nsolves + + # A working type CHOLMOD cannot factor keeps the plain matrix-free iteration. + A32 = Float32.([4.0 1.0 0.5; 1.0 3.0 1.0; 0.5 1.0 2.5]) + a32 = symcover_min(AbsLog{2}(), A32; linsolve=:lsqr) + @test a32 isa Vector{Float32} + @test a32 ≈ symcover_min(AbsLog{2}(), A32; linsolve=:dense) rtol=1e-5 + @test iscover(a32, A32; rtol=1e-5) +end + @testset "MMC disconnected-support gauge" begin # A support graph that splits into k connected components carries k independent # (e; −e) gauges. The asymmetric dense normal equations pin only the global one From a27703b6b4347010ecffffda21d3bd6ac8495e39 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 22 Aug 2026 19:35:34 -0500 Subject: [PATCH 6/7] Tighten the Woodbury and LSQR paths MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The conjugate-gradient residual recurrence drifts from `f − B x`, so convergence is confirmed against a freshly computed residual and a disagreement restarts the iteration; only a confirmed result is reported exact, which is what the sign-stability stopping test needs. A support thin enough per row can still carry a quadratic number of zeros, which would make the sparse correction dense work, so the Woodbury path also requires the total zero count to be O(n). The LSQR preconditioner is a plain scaling while its condition-number estimate stays low, forming and factorizing nothing. `stats.cholsolves` counts the Woodbury solves that reached the factorization, the two condition-number gates are named constants, and the Woodbury workspace is allocated only on that path. Assisted-by: Claude Fable 5 Assisted-by: Claude Opus 5 --- src/minimal_covers.jl | 398 +++++++++++++++++++++++++++-------------- src/sparse_support.jl | 7 +- test/minimal_covers.jl | 56 ++++-- 3 files changed, 313 insertions(+), 148 deletions(-) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index fc018e7..a29217e 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -21,14 +21,18 @@ Supported ϕ values: linear solve). `:dense` factorizes the reweighted normal equations densely, at O(n³) per Newton step. `:woodbury` solves the same equations as a sparse correction of the complete-support ones: the matrix is a sparse symmetric - positive-definite matrix plus `e*eᵀ`, so a sparse Cholesky and a - Sherman–Morrison update replace the dense factorization. It requires - `Float64` arithmetic and a support missing at most `n ÷ 4` entries in any - row, and raises an `ArgumentError` otherwise. `:lsqr` uses matrix-free LSQR - (per-iteration cost O(nnz), intended for large sparse supports), right - preconditioned in `Float64` by a sparse Cholesky of the diagonal of the - unweighted normal matrix plus the rows the penalty currently weights, which - keeps its iteration count from growing with the penalty strength. `:auto` + positive-definite matrix `C` plus `e*eᵀ`. Well-conditioned penalty stages + apply that sum without forming it and solve by Jacobi-preconditioned + conjugate gradients; the rest take a sparse Cholesky of `C` and a + Sherman–Morrison update. Both are exact to rounding. `:woodbury` requires + `Float64` arithmetic, a support missing at most `n ÷ 4` entries in any row, + and at most `4n` zero entries in total, and raises an `ArgumentError` + otherwise. `:lsqr` uses matrix-free LSQR (per-iteration cost O(nnz), + intended for large sparse supports), right preconditioned in `Float64` by + the diagonal of the unweighted normal matrix, joined by the rows the penalty + currently weights — through a sparse Cholesky — once diagonal scaling alone + would leave the system ill conditioned; this keeps its iteration count from + growing with the penalty strength. `:auto` selects `:woodbury` where its requirements hold and `:dense` elsewhere. `linsolve` defaults to `:auto` for dense `A`; the `SparseMatrixCSC`/`Symmetric`/`Hermitian` sparse methods default to `:lsqr` @@ -70,13 +74,17 @@ Supported ϕ values: linear solve). `:dense` factorizes the reweighted normal equations densely, at O((m+n)³) per Newton step. `:woodbury` solves the same equations as a sparse correction of the complete-support ones: the matrix is a sparse - symmetric positive-definite matrix plus a rank-two term, so a sparse Cholesky - and a Woodbury update replace the dense factorization. It requires `Float64` - arithmetic and a support missing at most `min(m, n) ÷ 4` entries in any row - or column, and raises an `ArgumentError` otherwise. `:lsqr` uses matrix-free + symmetric positive-definite matrix `C` plus a rank-two term. Well-conditioned + penalty stages apply that sum without forming it and solve by + Jacobi-preconditioned conjugate gradients; the rest take a sparse Cholesky of + `C` and a Woodbury update. Both are exact to rounding. `:woodbury` requires + `Float64` arithmetic, a support missing at most `min(m, n) ÷ 4` entries in + any row or column, and at most `4·max(m, n)` zero entries in total, and + raises an `ArgumentError` otherwise. `:lsqr` uses matrix-free LSQR (per-iteration cost O(nnz), intended for large sparse supports), right - preconditioned in `Float64` by a sparse Cholesky of the diagonal of the - unweighted normal matrix plus the rows the penalty currently weights, which + preconditioned in `Float64` by the diagonal of the unweighted normal matrix, + joined by the rows the penalty currently weights — through a sparse Cholesky + — once diagonal scaling alone would leave the system ill conditioned; this keeps its iteration count from growing with the penalty strength. `:auto` selects `:woodbury` where its requirements hold and `:dense` elsewhere. `linsolve` defaults to `:auto` for dense `A`; the @@ -305,33 +313,65 @@ function _prepare_cover_start!(a::AbstractVector, b::AbstractVector, A::Abstract end -# Inner linear solve for the AbsLog{2} MMC Newton steps. `:dense` forms and -# factorizes the reweighted normal equations densely, at O(n³) per step. +# Inner linear solve for the AbsLog{2} MMC Newton steps. +# +# `:dense` forms and factorizes the reweighted normal equations densely, at O(n³) +# per step. +# # `:woodbury` splits the same matrix as `C + U Uᵀ`, where `C` is sparse (its -# off-diagonal pattern is the zero set of `A` together with the currently violated -# entries) and symmetric positive definite, and `U` has one column (symmetric) or -# two (asymmetric); a sparse Cholesky of `C` plus a Woodbury update then costs far -# less than the dense factorization whenever `A` is close to fully supported. -# It has two exact sub-paths. `C + U Uᵀ` can also be applied without being formed, -# at O(n + |Z| + |V|) per application, and Gershgorin bounds its condition number by -# `1 + (κ−1)·2·maxdeg(V)/n` against the complete-support diagonal; while that -# estimate stays under 1000, Jacobi-preconditioned conjugate gradients converge to -# rounding in a few hundred such applications, which beats a factorization whose -# fill on the near-random violated pattern of the early stages approaches dense. -# Above it, the factorization runs. Either way the answer is exact to rounding, +# off-diagonal pattern is the zero set `Z` of `A` together with the currently +# violated entries `V`) and symmetric positive definite, and `U` has one column +# (symmetric) or two (asymmetric). It has two sub-paths, both exact to rounding, # which is what the sign-stability stopping test in the continuation loop requires. -# `:auto` takes `:woodbury` where it applies and `:dense` otherwise. `:lsqr` forces -# the matrix-free path, whose per-iteration cost is O(nnz); it is the intended -# solve for large sparse supports (where nnz ≪ n²) and is used by the -# structured/sparse methods. It is right preconditioned, which is what keeps its -# iteration count from growing as the continuation raises κ. +# A sparse Cholesky of `C` plus a Sherman–Morrison (symmetric) or Woodbury +# (asymmetric) update costs far less than the dense factorization whenever `A` is +# close to fully supported. `C + U Uᵀ` can alternatively be applied without being +# formed, at O(n + |Z| + |V|) per application; Gershgorin on `(κ−1)·L_V` against the +# complete-support diagonal gives `1 + (κ−1)·2·maxdeg(V)/n` as an estimate of its +# condition number (the sharp bound is a small multiple of that), and while the +# estimate stays under `WOODBURY_CG_KAPPA`, Jacobi-preconditioned conjugate +# gradients converge to rounding in a few hundred such applications — cheaper than a +# factorization whose fill, on the near-random violated pattern of the early stages, +# approaches dense. Above it the factorization runs, as it does for any CG run that +# exhausts its iteration cap. # # `C` is positive definite because the complete-support matrix contributes `n` (or -# `m`) to each diagonal while the zero set subtracts a signless Laplacian `L_Z` -# with λmax(L_Z) ≤ 2·maxdeg(Z); requiring at most a quarter of a row to be zero -# keeps the difference bounded below by half the diagonal. CHOLMOD is the sparse -# factorization behind it, and it is reliable only in `Float64`, so that is the -# only working type the path accepts. +# `m`) to each diagonal while the zero set subtracts a signless Laplacian `L_Z` with +# λmax(L_Z) ≤ 2·maxdeg(Z); requiring at most a quarter of a row to be zero keeps the +# difference bounded below by half the diagonal. A second requirement is about cost +# rather than definiteness: `Z` enters every matvec and every factorization, so the +# path is taken only while the total number of zeros is O(n). CHOLMOD is the sparse +# factorization behind it, and it is reliable only in `Float64`, so that is the only +# working type the path accepts. +# +# `:auto` takes `:woodbury` where it applies and `:dense` otherwise. +# +# `:lsqr` forces the matrix-free path, whose per-iteration cost is O(nnz); it is the +# intended solve for large sparse supports (where nnz ≪ n²) and is used by the +# structured/sparse methods. In `Float64` it is right preconditioned, which is what +# keeps its iteration count from growing as the continuation raises κ. The +# preconditioner is `M = diag(RᵀR) + (κ−1)·Σ_{e∈V} rₑ·rₑᵀ`: the diagonal of the +# unweighted normal matrix, together with the exact contribution of the rows LSQR +# weights by κ. All of the κ-dependence of `RᵀWR` sits in those rows, and every +# generalized eigenvalue of `(RᵀWR, M)` is a mediant of eigenvalues of +# `(RᵀR, diag(RᵀR))` and so lies in their range. `M = K·Kᵀ` and LSQR runs on +# `√W·R·K⁻ᵀ` in the variable `y = Kᵀ·x`. While the same condition-number estimate, +# taken against the unweighted diagonal, stays under `LSQR_PRECOND_KAPPA` the +# violated rows are left out and `K` is the diagonal `sqrt.(diag(RᵀR))`, applied +# without forming anything; above it they are included and `K` is the permuted sparse +# Cholesky factor of `M`, applied through the CHOLMOD factor components `F.PtL` and +# `F.UP`. `Kᵀ` is never needed as a product: the warm start `Kᵀ·x₀` is `K⁻¹·(M·x₀)`, +# which those same components and one sparse matrix-vector product supply. + +# Condition-number estimate above which a Woodbury solve is factorized rather than +# iterated: past it conjugate gradients need more applications than the sparse +# Cholesky costs. +const WOODBURY_CG_KAPPA = 1000 + +# Condition-number estimate above which the LSQR preconditioner takes in the rows the +# penalty currently weights; below it diagonal scaling alone leaves the system well +# enough conditioned, and no factorization is formed. +const LSQR_PRECOND_KAPPA = 1000 # `B*x` for the symmetric Woodbury system # `B = dbase·I + e·eᵀ − L_Z + dκ·L_V` (`dbase = n + ridge`), applied without forming @@ -384,19 +424,6 @@ function _woodbury_mul!(y, x, m, drow, dcol, zedges, vedges, dκ) return y end -# Right preconditioner for the LSQR path. `M = diag(RᵀR) + (κ−1)·Σ_{e∈V} rₑ·rₑᵀ` -# holds the diagonal of the unweighted normal matrix together with the exact -# contribution of the rows LSQR weights by κ. All of the κ-dependence of `RᵀWR` sits -# in those rows, and every generalized eigenvalue of `(RᵀWR, M)` is a mediant of -# eigenvalues of `(RᵀR, diag(RᵀR))` and so lies in their range, which is what keeps -# the preconditioned iteration count from growing with κ. -# -# `M = K·Kᵀ` with `K` the permuted Cholesky factor `PtL` of `M`. LSQR then runs on -# `√W·R·K⁻ᵀ` in the variable `y = Kᵀ·x`, applying `K⁻¹` and `K⁻ᵀ` through the CHOLMOD -# factor components `F.PtL` and `F.UP`. `Kᵀ` itself is never needed as a product: the -# warm start is `Kᵀ·x₀ = K⁻¹·(M·x₀)`, which the same two objects and one sparse -# matrix-vector product supply. - # Jacobi-preconditioned conjugate gradients for the symmetric positive-definite # Woodbury system `B x = f`, with `Bmul!(y, x)` applying `B` and `dg` holding its # diagonal. `x` carries the warm start in and the iterate out; `r`, `z`, `d`, `Ad` @@ -406,30 +433,48 @@ end # stopping test to mean what it says, so `tol` sits at the level of `eps` and a run # that exhausts `maxiter` reports failure instead of a partial answer; the caller # then falls back to the factorization, which is exact. +# +# `r` is carried by a recurrence that drifts from `f − B x`, so success is never +# declared on it: a claim of convergence is confirmed against a freshly computed +# residual, and a disagreement restarts the iteration there. The confirming +# application counts against `maxiter` like any other. function _pcg!(Bmul!, x, dg, f, r, z, d, Ad, maxiter::Int, tol) + iters = 0 Bmul!(r, x) @. r = f - r nrm = norm(r) + fresh = true # `r` holds `f − B x`, not the recurrence's estimate of it @. z = r / dg copyto!(d, z) rz = dot(r, z) - iters = 0 - for k in 1:maxiter - nrm <= tol && break - iters = k + while iters < maxiter + if nrm <= tol + fresh && return iters, true + iters += 1 + Bmul!(Ad, x) + @. r = f - Ad + nrm = norm(r) + fresh = true + @. z = r / dg + copyto!(d, z) + rz = dot(r, z) + continue + end + iters += 1 Bmul!(Ad, d) dAd = dot(d, Ad) dAd > 0 || break a = rz / dAd @. x += a * d @. r -= a * Ad + fresh = false nrm = norm(r) @. z = r / dg rznew = dot(r, z) @. d = z + (rznew / rz) * d rz = rznew end - return iters, nrm <= tol + return iters, fresh && nrm <= tol end # Matrix-free LSQR (Paige & Saunders) for the weighted least-squares problem @@ -440,9 +485,9 @@ end # strength κ) rather than that of `MᵀM` (≈ κ); at κ = 1e8 the squared conditioning # breaks CG while LSQR stays accurate. # -# The callers pass a right-preconditioned operator, so `x` here is the preconditioned -# variable and `M` is `√W·R·K⁻ᵀ`; the preconditioner is described above `_lsqr`'s -# callers in each worker. +# A `Float64` caller passes a right-preconditioned operator, so `x` is then the +# preconditioned variable and `M` is `√W·R·K⁻ᵀ`; see the description of `:lsqr` in the +# inner-solve overview above. Other working types pass `√W·R` itself. # # The penalty least-squares problem is inconsistent (its optimal residual is # nonzero), so the stopping test is on the normal-equations residual @@ -502,9 +547,10 @@ function _lsqr(Amul!, Atmul!, b::AbstractVector{T}, x0::AbstractVector{T}; end # Worker for `symcover_min(::AbsLog{2})`. Returns `(a, stats)` where `stats` is a -# NamedTuple `(; nsolves, lsqriters, cgiters, linsolve)` recording the number of inner -# linear solves, the total LSQR and conjugate-gradient iterations (0 on paths that run -# neither), and which path ran. +# NamedTuple `(; nsolves, lsqriters, cgiters, cholsolves, linsolve)` recording the +# number of inner linear solves, the total LSQR and conjugate-gradient iterations (0 on +# paths that run neither), how many Woodbury solves fell to the sparse factorization, +# and which path ran. # `linsolve` reports the path that ran: `:dense`, `:woodbury`, or `:lsqr`. # `start`, when given, is a positive cover of `A` # indexed like `axes(A, 1)` and supplies the first iterate in place of the cold @@ -548,16 +594,25 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), maxzero = max(maxzero, n - length(slots)) end ne = length(edges) - # The Woodbury path splits the normal equations around the complete-support - # matrix `n·I + e·eᵀ`, so its cost is set by the zero set `Z` rather than by `n`, - # and `n·I − L_Z` is positive definite only while `Z` stays thin. + # The Woodbury path splits the normal equations around the complete-support matrix + # `n·I + e·eᵀ`, so its cost is set by the zero set `Z` rather than by `n`. Two + # separate conditions gate it. Per row: `n·I − L_Z` is positive definite only + # while no row carries more than `n ÷ 4` zeros. In total: `Z` is materialized and + # then traversed by every matvec and every factorization, so the path is worth + # taking only while `|Z|` stays O(n) — a support that is merely thin per row can + # still carry Θ(n²) zeros, and the split would then be dense work under a name + # that promises otherwise. + nzero = n * n - ne + zbudget = 4 * n use_woodbury = false if !use_lsqr && linsolve !== :dense - ok = T === Float64 && maxzero <= n ÷ 4 + ok = T === Float64 && maxzero <= n ÷ 4 && nzero <= zbudget if linsolve === :woodbury && !ok T === Float64 || throw(ArgumentError("linsolve=:woodbury requires Float64 arithmetic, but `A` works in $T; use :dense or :lsqr")) - throw(ArgumentError("linsolve=:woodbury requires every row of `A` to have at most n ÷ 4 = $(n ÷ 4) zeros; got $maxzero")) + maxzero <= n ÷ 4 || + throw(ArgumentError("linsolve=:woodbury requires every row of `A` to have at most n ÷ 4 = $(n ÷ 4) zeros; got $maxzero")) + throw(ArgumentError("linsolve=:woodbury requires `A` to have at most 4n = $zbudget zeros in total; got $nzero")) end use_woodbury = ok end @@ -624,9 +679,9 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # Entries the frozen weights of the current solve treat as violated. A full Newton # step that leaves this pattern intact has landed on the stage's minimizer. vpat = falses(ne) - vedges = Tuple{Int,Int}[] # the violated entries of the current solve - degV = zeros(Int, n) # violated entries per row - dg = zeros(T, n) # diagonal of `B`, for the ridge and the CG preconditioner + vedges = Tuple{Int,Int}[] # the violated entries of the current solve + degV = zeros(Int, use_woodbury ? n : 0) # violated entries per row + dg = zeros(T, use_woodbury ? n : 0) # diagonal of `B`, for the ridge and the CG preconditioner # Diagonal of the unweighted normal matrix `RᵀR`, the base of the LSQR # preconditioner: each directed support entry puts 1 at each of its ends, and a # diagonal entry, whose row of `R` is `2·e_p`, puts 4. A support-free variable is @@ -651,19 +706,22 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), Mv = T[] px = zeros(T, use_lsqr ? n : 0) # scale vector recovered from the LSQR variable pg = zeros(T, use_lsqr ? n : 0) # `Rᵀ√W y` before the preconditioner is applied + # `K` of the diagonal preconditioner, which is κ-independent and so built once. + psqrt = use_lsqr ? sqrt.(dpart) : T[] # COO triplets of `C`, refilled whenever a Woodbury solve is factorized. Ci = Int[] Cj = Int[] Cv = T[] - rhs = zeros(T, n, 2) - cgx = zeros(T, n) - cgr = zeros(T, n) - cgz = zeros(T, n) - cgd = zeros(T, n) - cgAd = zeros(T, n) + rhs = zeros(T, use_woodbury ? n : 0, 2) + cgx = zeros(T, use_woodbury ? n : 0) + cgr = zeros(T, use_woodbury ? n : 0) + cgz = zeros(T, use_woodbury ? n : 0) + cgd = zeros(T, use_woodbury ? n : 0) + cgAd = zeros(T, use_woodbury ? n : 0) nsolves = Ref(0) nlsqr = Ref(0) ncg = Ref(0) + nchol = Ref(0) solve_weighted = function (α, κ) nsolves[] += 1 if use_lsqr @@ -696,6 +754,30 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), for p in 1:n κest = max(κest, oneunit(T) + dκ * 2 * mdiag[p] / dpart[p]) end + if κest <= LSQR_PRECOND_KAPPA + # `K` is diagonal here, so it is applied by a scaling and nothing + # is assembled or factorized. + Dmul! = function (y, yv) + @. px = yv / psqrt + for (e, (ip, jp)) in enumerate(edges) + y[e] = ws[e] * (px[ip] + px[jp]) + end + return y + end + Dtmul! = function (z, y) + fill!(pg, zero(T)) + for (e, (ip, jp)) in enumerate(edges) + t = ws[e] * y[e] + pg[ip] += t + pg[jp] += t + end + @. z = pg / psqrt + return z + end + soly, it = _lsqr(Dmul!, Dtmul!, cv, psqrt .* α) + nlsqr[] += it + return soly ./ psqrt + end empty!(Mi) empty!(Mj) empty!(Mv) @@ -704,32 +786,33 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), push!(Mj, p) push!(Mv, dpart[p]) end - if κest > 1000 - for (p, q) in vedges - if p == q - push!(Mi, p) - push!(Mj, p) - push!(Mv, 4 * dκ) - else - push!(Mi, p) - push!(Mj, p) - push!(Mv, dκ) - push!(Mi, q) - push!(Mj, q) - push!(Mv, dκ) - push!(Mi, p) - push!(Mj, q) - push!(Mv, dκ) - push!(Mi, q) - push!(Mj, p) - push!(Mv, dκ) - end + for (p, q) in vedges + if p == q + push!(Mi, p) + push!(Mj, p) + push!(Mv, 4 * dκ) + else + push!(Mi, p) + push!(Mj, p) + push!(Mv, dκ) + push!(Mi, q) + push!(Mj, q) + push!(Mv, dκ) + push!(Mi, p) + push!(Mj, q) + push!(Mv, dκ) + push!(Mi, q) + push!(Mj, p) + push!(Mv, dκ) end end Msp = sparse(Mi, Mj, Mv, n, n) MF = cholesky(Symmetric(Msp)) Kc = MF.PtL Uc = MF.UP + # CHOLMOD exposes no in-place solve for a factor component, so each + # application returns a fresh vector; the transpose product copies it + # into the buffer LSQR hands over, which is the only copy avoidable here. Pmul! = function (y, yv) xv = Uc \ yv for (e, (ip, jp)) in enumerate(edges) @@ -750,7 +833,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), mul!(px, Msp, α) soly, it = _lsqr(Pmul!, Ptmul!, cv, Kc \ px) nlsqr[] += it - return Uc \ soly + return (Uc \ soly)::Vector{T} end Amul! = function (y, x) for (e, (ip, jp)) in enumerate(edges) @@ -805,14 +888,14 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), for p in 1:n dg[p] += oneunit(T) + ridge end - # Gershgorin on `(κ−1)·L_V` against a diagonal of at least `n` bounds the - # condition number of `B`. While that bound is small the structured matvec - # plus conjugate gradients reaches the same answer in a few hundred O(n + - # |Z| + |V|) iterations, which is far cheaper than a factorization whose - # fill, on the near-random violated pattern of the early stages, is close - # to dense. + # Gershgorin on `(κ−1)·L_V` against a diagonal of at least `n` estimates + # the condition number of `B`. While that estimate is small the structured + # matvec plus conjugate gradients reaches the same answer in a few hundred + # O(n + |Z| + |V|) iterations, which is far cheaper than a factorization + # whose fill, on the near-random violated pattern of the early stages, is + # close to dense. κest = oneunit(T) + dκ * 2 * maxdegV / n - if κest <= 1000 + if κest <= WOODBURY_CG_KAPPA copyto!(cgx, α) dbase = T(n) + ridge Bmul! = (y, x) -> _symwoodbury_mul!(y, x, dbase, zedges, vedges, dκ) @@ -822,6 +905,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), ok && return copy(cgx) end # `sparse` sums the duplicate triplets; the ridge rides on `C`'s diagonal. + nchol[] += 1 empty!(Ci) empty!(Cj) empty!(Cv) @@ -927,11 +1011,13 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), a[i] = hassupp[ip] ? exp(α[ip] + γ) : zero(T) end return a, (; nsolves=nsolves[], lsqriters=nlsqr[], cgiters=ncg[], + cholsolves=nchol[], linsolve=(use_lsqr ? :lsqr : use_woodbury ? :woodbury : :dense)) end # Worker for `cover_min(::AbsLog{2})`. Returns `(a, b, stats)` with `stats` a -# NamedTuple `(; nsolves, lsqriters, cgiters, linsolve)` (see `_symcover_min_abslog2`). +# NamedTuple `(; nsolves, lsqriters, cgiters, cholsolves, linsolve)` (see +# `_symcover_min_abslog2`). # `start`, when given, is a positive cover `(a, b)` indexed like the rows and columns # of `A`, supplying the first iterate in place of the cold unweighted solve. function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), @@ -985,13 +1071,23 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), maxzero = max(maxzero, m - nzcol[jp]) end zbound = min(m, n) ÷ 4 + # Two separate conditions gate the path. Per row and column: `D − L_Z` is positive + # definite only while neither carries more than `min(m, n) ÷ 4` zeros. In total: + # `Z` is materialized and then traversed by every matvec and every factorization, + # so the path is worth taking only while `|Z|` stays O(m + n) — a support that is + # merely thin per row can still carry Θ(m·n) zeros, and the split would then be + # dense work under a name that promises otherwise. + nzero = m * n - ne + zbudget = 4 * max(m, n) use_woodbury = false if !use_lsqr && linsolve !== :dense - ok = T === Float64 && maxzero <= zbound + ok = T === Float64 && maxzero <= zbound && nzero <= zbudget if linsolve === :woodbury && !ok T === Float64 || throw(ArgumentError("linsolve=:woodbury requires Float64 arithmetic, but `A` works in $T; use :dense or :lsqr")) - throw(ArgumentError("linsolve=:woodbury requires every row and column of `A` to have at most min(m, n) ÷ 4 = $zbound zeros; got $maxzero")) + maxzero <= zbound || + throw(ArgumentError("linsolve=:woodbury requires every row and column of `A` to have at most min(m, n) ÷ 4 = $zbound zeros; got $maxzero")) + throw(ArgumentError("linsolve=:woodbury requires `A` to have at most 4·max(m, n) = $zbudget zeros in total; got $nzero")) end use_woodbury = ok end @@ -1073,9 +1169,9 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # Entries the frozen weights of the current solve treat as violated. A full Newton # step that leaves this pattern intact has landed on the stage's minimizer. vpat = falses(ne) - vedges = Tuple{Int,Int}[] # the violated entries of the current solve - degV = zeros(Int, N) # violated entries per row and per column - dg = zeros(T, N) # diagonal of `B`, for the ridge and the CG preconditioner + vedges = Tuple{Int,Int}[] # the violated entries of the current solve + degV = zeros(Int, use_woodbury ? N : 0) # violated entries per row and per column + dg = zeros(T, use_woodbury ? N : 0) # diagonal of `B`, for the ridge and the CG preconditioner # Diagonal of the unweighted normal matrix of the gauge-augmented system, # `RᵀR + v0·v0ᵀ`: the support degree at each position, plus the gauge row's 1. A # support-free variable takes that 1 alone, which keeps the preconditioner @@ -1096,19 +1192,22 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), Mv = T[] px = zeros(T, use_lsqr ? N : 0) # scale vector recovered from the LSQR variable pg = zeros(T, use_lsqr ? N : 0) # `Rᵀ√W y` before the preconditioner is applied + # `K` of the diagonal preconditioner, which is κ-independent and so built once. + psqrt = use_lsqr ? sqrt.(dpart) : T[] # COO triplets of `C`, refilled whenever a Woodbury solve is factorized. Ci = Int[] Cj = Int[] Cv = T[] - rhs = zeros(T, N, 3) - cgx = zeros(T, N) - cgr = zeros(T, N) - cgz = zeros(T, N) - cgd = zeros(T, N) - cgAd = zeros(T, N) + rhs = zeros(T, use_woodbury ? N : 0, 3) + cgx = zeros(T, use_woodbury ? N : 0) + cgr = zeros(T, use_woodbury ? N : 0) + cgz = zeros(T, use_woodbury ? N : 0) + cgd = zeros(T, use_woodbury ? N : 0) + cgAd = zeros(T, use_woodbury ? N : 0) nsolves = Ref(0) nlsqr = Ref(0) ncg = Ref(0) + nchol = Ref(0) solve_weighted = function (x, κ) nsolves[] += 1 if use_lsqr @@ -1138,6 +1237,32 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), for p in 1:N κest = max(κest, oneunit(T) + dκ * 2 * mdiag[p] / dpart[p]) end + if κest <= LSQR_PRECOND_KAPPA + # `K` is diagonal here, so it is applied by a scaling and nothing + # is assembled or factorized. + Dmul! = function (y, yv) + @. px = yv / psqrt + for (e, (p, q)) in enumerate(edges) + y[e] = ws[e] * (px[p] + px[q]) + end + y[g] = dot(v0, px) + return y + end + Dtmul! = function (z, y) + fill!(pg, zero(T)) + for (e, (p, q)) in enumerate(edges) + t = ws[e] * y[e] + pg[p] += t + pg[q] += t + end + @. pg += v0 * y[g] + @. z = pg / psqrt + return z + end + soly, it = _lsqr(Dmul!, Dtmul!, cv, psqrt .* x) + nlsqr[] += it + return soly ./ psqrt + end empty!(Mi) empty!(Mj) empty!(Mv) @@ -1146,26 +1271,27 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), push!(Mj, p) push!(Mv, dpart[p]) end - if κest > 1000 - for (p, q) in vedges - push!(Mi, p) - push!(Mj, p) - push!(Mv, dκ) - push!(Mi, q) - push!(Mj, q) - push!(Mv, dκ) - push!(Mi, p) - push!(Mj, q) - push!(Mv, dκ) - push!(Mi, q) - push!(Mj, p) - push!(Mv, dκ) - end + for (p, q) in vedges + push!(Mi, p) + push!(Mj, p) + push!(Mv, dκ) + push!(Mi, q) + push!(Mj, q) + push!(Mv, dκ) + push!(Mi, p) + push!(Mj, q) + push!(Mv, dκ) + push!(Mi, q) + push!(Mj, p) + push!(Mv, dκ) end Msp = sparse(Mi, Mj, Mv, N, N) MF = cholesky(Symmetric(Msp)) Kc = MF.PtL Uc = MF.UP + # CHOLMOD exposes no in-place solve for a factor component, so each + # application returns a fresh vector; the transpose product copies it + # into the buffer LSQR hands over, which is the only copy avoidable here. Pmul! = function (y, yv) xv = Uc \ yv for (e, (p, q)) in enumerate(edges) @@ -1188,7 +1314,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), mul!(px, Msp, x) soly, it = _lsqr(Pmul!, Ptmul!, cv, Kc \ px) nlsqr[] += it - return Uc \ soly + return (Uc \ soly)::Vector{T} end Amul! = function (y, xx) for (e, (p, q)) in enumerate(edges) @@ -1249,14 +1375,14 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), for p in 1:N dg[p] += oneunit(T) + ridge end - # Gershgorin on `(κ−1)·L_V` against the smaller diagonal block bounds the - # condition number of `B`. While that bound is small the structured matvec - # plus conjugate gradients reaches the same answer in a few hundred + # Gershgorin on `(κ−1)·L_V` against the smaller diagonal block estimates + # the condition number of `B`. While that estimate is small the structured + # matvec plus conjugate gradients reaches the same answer in a few hundred # O(N + |Z| + |V|) iterations, which is far cheaper than a factorization # whose fill on the near-random violated pattern of the early stages # approaches dense. κest = oneunit(T) + dκ * 2 * maxdegV / min(m, n) - if κest <= 1000 + if κest <= WOODBURY_CG_KAPPA copyto!(cgx, x) drow = T(n) + ridge dcol = T(m) + ridge @@ -1267,6 +1393,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), ok && return copy(cgx) end # `sparse` sums the duplicate triplets; the ridge rides on `C`'s diagonal. + nchol[] += 1 empty!(Ci) empty!(Cj) empty!(Cv) @@ -1452,6 +1579,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), b[j] = hascol[jp] ? exp(x[m+jp] - s[colcomp[jp]]) : zero(T) end return a, b, (; nsolves=nsolves[], lsqriters=nlsqr[], cgiters=ncg[], + cholsolves=nchol[], linsolve=(use_lsqr ? :lsqr : use_woodbury ? :woodbury : :dense)) end diff --git a/src/sparse_support.jl b/src/sparse_support.jl index bef8c67..5807b69 100644 --- a/src/sparse_support.jl +++ b/src/sparse_support.jl @@ -80,9 +80,10 @@ end # Native AbsLog{2} MMC solvers on sparse supports default to the matrix-free LSQR # inner solve, whose per-iteration cost is O(nnz) and whose accuracy tracks the # conditioning of √W·R (≈ √κ) rather than that of the normal equations (≈ κ). This -# is the intended path when nnz ≪ n²; pass `linsolve=:auto`/`:dense` to force the -# dense factorization. Only AbsLog{2} is native; other penalties dispatch to the -# JuMP extension. +# is the intended path when nnz ≪ n². Pass `linsolve=:dense` to force the dense +# factorization, or `linsolve=:auto` to let the solver choose between the Woodbury +# split and the dense factorization. Only AbsLog{2} is native; other penalties +# dispatch to the JuMP extension. function symcover_min(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) a, _ = _symcover_min_abslog2(A; linsolve, kwargs...) return a diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index a397ecd..7f9eb99 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -155,14 +155,18 @@ end @testset "symmetric, n = $n" for n in (6, 30, 120) A = symlognormal(n) - variants = ["all nonzero" => A, - "zero diagonal" => A - Diagonal(A)] - # A sparse symmetric zero set, thin enough to stay inside the n ÷ 4 guard. + # A symmetric zero set placed to sit inside both guards at every size tested: + # pairing consecutive indices gives exactly one zero per row, against a + # per-row allowance of `n ÷ 4` (which is 1 already at n = 6) and a total + # allowance of `4n`. Z = symlognormal(n) - mask = rand(rng, n, n) .< 0.5 / max(n ÷ 8, 1) - mask = mask .| mask' - Z[mask] .= 0.0 - maximum(count(iszero, Z; dims=2)) <= n ÷ 4 && push!(variants, "sparse zeros" => Z) + for k in 1:(n ÷ 2) + Z[2k-1, 2k] = 0.0 + Z[2k, 2k-1] = 0.0 + end + variants = ["all nonzero" => A, + "zero diagonal" => A - Diagonal(A), + "paired zeros" => Z] for (name, M) in variants @testset "$name" begin ad, sd = MatrixCovers._symcover_min_abslog2(M; linsolve=:dense) @@ -178,7 +182,9 @@ end # are well enough conditioned for conjugate gradients, the late ones # are not, and both are exact. @test sw.cgiters > 0 + @test sw.cholsolves > 0 @test sd.cgiters == 0 + @test sd.cholsolves == 0 end end end @@ -199,7 +205,9 @@ end @test aw .* bw' ≈ ad .* bd' rtol=1e-7 @test iscover(aw, bw, M; atol=1e-7) @test sw.cgiters > 0 + @test sw.cholsolves > 0 @test sd.cgiters == 0 + @test sd.cholsolves == 0 end end end @@ -230,18 +238,22 @@ end am, bm = cover_min(AbsLog{2}(), Matrix(Agv); linsolve=:woodbury) @test av .* bv' ≈ am .* bm' rtol=1e-10 - # A single well-conditioned stage exercises the conjugate-gradient sub-path alone. + # A single stage at κ = 1e2 stays inside the conjugate-gradient regime throughout, + # so the factorization is never reached. A1 = symlognormal(24) c1, s1 = MatrixCovers._symcover_min_abslog2(A1; κs=(1e2,), linsolve=:woodbury) @test s1.cgiters > 0 + @test s1.cholsolves == 0 @test c1 ≈ MatrixCovers._symcover_min_abslog2(A1; κs=(1e2,), linsolve=:dense)[1] rtol=1e-8 G1 = lognormal(24, 18) p1, q1, t1 = MatrixCovers._cover_min_abslog2(G1; κs=(1e2,), linsolve=:woodbury) pd1, qd1, _ = MatrixCovers._cover_min_abslog2(G1; κs=(1e2,), linsolve=:dense) @test t1.cgiters > 0 + @test t1.cholsolves == 0 @test p1 .* q1' ≈ pd1 .* qd1' rtol=1e-8 - # The zero set must stay inside the guard, and the arithmetic must be Float64. + # No row may carry more than a quarter zeros — that is what keeps `C` positive + # definite — and the arithmetic must be Float64. holey = symlognormal(12) holey[1, 1:5] .= 0.0 holey[1:5, 1] .= 0.0 @@ -252,6 +264,23 @@ end @test_throws "at most min(m, n) ÷ 4 = 3 zeros; got 5" cover_min(AbsLog{2}(), gholey; linsolve=:woodbury) @test MatrixCovers._cover_min_abslog2(gholey)[3].linsolve === :dense + # A support thin enough per row can still carry a quadratic number of zeros, which + # would make the "sparse" correction dense work; the total budget rejects it. + wide = symlognormal(40) + for i in 1:40, j in 1:40 + (i != j && (i + j) % 4 == 0) && (wide[i, j] = 0.0) + end + @test maximum(count(iszero, wide; dims=2)) <= 40 ÷ 4 + @test count(iszero, wide) > 4 * 40 + @test_throws "at most 4n = 160 zeros in total" symcover_min(AbsLog{2}(), wide; linsolve=:woodbury) + @test MatrixCovers._symcover_min_abslog2(wide)[2].linsolve === :dense + gwide = lognormal(40, 40) + for i in 1:40, j in 1:40 + (i + j) % 4 == 0 && (gwide[i, j] = 0.0) + end + @test_throws "at most 4·max(m, n) = 160 zeros in total" cover_min(AbsLog{2}(), gwide; linsolve=:woodbury) + @test MatrixCovers._cover_min_abslog2(gwide)[3].linsolve === :dense + A32 = Float32.(symlognormal(8)) @test_throws "requires Float64 arithmetic" symcover_min(AbsLog{2}(), A32; linsolve=:woodbury) @test MatrixCovers._symcover_min_abslog2(A32)[2].linsolve === :dense @@ -273,7 +302,10 @@ end al, sl = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr) @test ad ≈ al rtol=1e-6 @test aw ≈ al rtol=1e-6 - # One solve per κ stage is saved; `κs` has four stages by default. + # One solve per κ stage is saved; `κs` has four stages by default. The absolute + # bound guards against a regression in the count itself: this matrix takes 24 + # solves on the exact paths against 28 on `:lsqr`, so 26 leaves two solves of + # headroom while still failing if the early stop stops firing. @test sd.nsolves == sw.nsolves @test sd.nsolves <= sl.nsolves - length((1e2, 1e4, 1e6, 1e8)) @test sd.nsolves <= 26 @@ -286,6 +318,7 @@ end @test gw .* hw' ≈ gl .* hl' rtol=1e-6 @test td.nsolves == tw.nsolves @test td.nsolves <= tl.nsolves - length((1e2, 1e4, 1e6, 1e8)) + # 22 solves measured here against 26 on `:lsqr`; 24 leaves two of headroom. @test td.nsolves <= 24 end @@ -299,12 +332,15 @@ end ad, _ = MatrixCovers._symcover_min_abslog2(A; linsolve=:dense) al, sl = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr) @test al ≈ ad rtol=1e-6 + # 29.4 iterations per solve measured here; the bound doubles that, and an + # unpreconditioned run would sit in the hundreds by the last κ stage. @test sl.lsqriters <= 60 * sl.nsolves G = exp.(randn(rng, 120, 90)) gd, hd, _ = MatrixCovers._cover_min_abslog2(G; linsolve=:dense) gl, hl, tl = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr) @test gl .* hl' ≈ gd .* hd' rtol=1e-6 + # 41.4 iterations per solve measured here, against the same bound. @test tl.lsqriters <= 60 * tl.nsolves # A working type CHOLMOD cannot factor keeps the plain matrix-free iteration. From c45171894c2f69365ad7cce579f5eb6c1ec6f6dd Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 23 Aug 2026 03:18:01 -0500 Subject: [PATCH 7/7] Share one Woodbury solve between both workers `C` is assembled sparsely on every Woodbury solve, in full rather than in one triangle, so the conjugate-gradient sub-path applies it with `mul!` and adds the low-rank term directly instead of walking the edge lists. The symmetric Sherman-Morrison and the asymmetric 2x2 Woodbury update become one `_woodbury_solve!` over a low-rank block of k columns, taking `[f U]` through the factor in a single multi-right-hand-side solve. Assisted-by: Claude Fable 5 Assisted-by: Claude Opus 5 --- src/minimal_covers.jl | 252 ++++++++++++++++++------------------------ 1 file changed, 109 insertions(+), 143 deletions(-) diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index a29217e..e9af960 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -321,19 +321,19 @@ end # `:woodbury` splits the same matrix as `C + U Uᵀ`, where `C` is sparse (its # off-diagonal pattern is the zero set `Z` of `A` together with the currently # violated entries `V`) and symmetric positive definite, and `U` has one column -# (symmetric) or two (asymmetric). It has two sub-paths, both exact to rounding, -# which is what the sign-stability stopping test in the continuation loop requires. -# A sparse Cholesky of `C` plus a Sherman–Morrison (symmetric) or Woodbury -# (asymmetric) update costs far less than the dense factorization whenever `A` is -# close to fully supported. `C + U Uᵀ` can alternatively be applied without being -# formed, at O(n + |Z| + |V|) per application; Gershgorin on `(κ−1)·L_V` against the -# complete-support diagonal gives `1 + (κ−1)·2·maxdeg(V)/n` as an estimate of its -# condition number (the sharp bound is a small multiple of that), and while the -# estimate stays under `WOODBURY_CG_KAPPA`, Jacobi-preconditioned conjugate -# gradients converge to rounding in a few hundred such applications — cheaper than a -# factorization whose fill, on the near-random violated pattern of the early stages, -# approaches dense. Above it the factorization runs, as it does for any CG run that -# exhausts its iteration cap. +# (symmetric) or two (asymmetric). `C` is assembled sparsely on every such solve, and +# two sub-paths then take it, both exact to rounding, which is what the +# sign-stability stopping test in the continuation loop requires. A sparse Cholesky +# of `C` plus a Woodbury update — Sherman–Morrison, in the one-column symmetric case +# — costs far less than the dense factorization whenever `A` is close to fully +# supported. Alternatively `C + U Uᵀ` is applied as `C·x` plus the low-rank term, at +# O(nnz(C)) per application; Gershgorin on `(κ−1)·L_V` against the complete-support +# diagonal gives `1 + (κ−1)·2·maxdeg(V)/n` as an estimate of its condition number +# (the sharp bound is a small multiple of that), and while the estimate stays under +# `WOODBURY_CG_KAPPA`, Jacobi-preconditioned conjugate gradients converge to rounding +# in a few hundred such applications — cheaper than a factorization whose fill, on the +# near-random violated pattern of the early stages, approaches dense. Above it the +# factorization runs, as it does for any CG run that exhausts its iteration cap. # # `C` is positive definite because the complete-support matrix contributes `n` (or # `m`) to each diagonal while the zero set subtracts a signless Laplacian `L_Z` with @@ -373,55 +373,26 @@ const WOODBURY_CG_KAPPA = 1000 # enough conditioned, and no factorization is formed. const LSQR_PRECOND_KAPPA = 1000 -# `B*x` for the symmetric Woodbury system -# `B = dbase·I + e·eᵀ − L_Z + dκ·L_V` (`dbase = n + ridge`), applied without forming -# a matrix: `L_X·x` accumulates `x[p] + x[q]` into `y[p]` over the directed edges of -# `X`, so one application costs O(n + |Z| + |V|) rather than O(nnz). -function _symwoodbury_mul!(y, x, dbase, zedges, vedges, dκ) - s = zero(eltype(y)) - for p in eachindex(x) - s += x[p] +# `(C + U·Uᵀ) x = f` solved from a factorization `F` of the sparse `C`, by the +# Woodbury identity `x = y − Y·((I + Uᵀ·Y) \ (Uᵀ·y))` with `y = C\f` and `Y = C\U`. +# One multi-right-hand-side solve of `[f U]` supplies both, and the capacitance is +# `k×k` for `U` of `k` columns: `k = 1` for the symmetric gauge `e`, where this is +# Sherman–Morrison, and `k = 2` for the asymmetric row and column indicators. `rhs` +# is the `size(U, 1)×(k+1)` buffer the block right-hand side is staged in. +function _woodbury_solve!(x, F, U, f, rhs) + k = size(U, 2) + copyto!(view(rhs, :, 1), f) + copyto!(view(rhs, :, 2:k+1), U) + sol = F \ rhs + y = view(sol, :, 1) + Y = view(sol, :, 2:k+1) + K = U' * Y + for i in axes(K, 1) + K[i, i] += oneunit(eltype(K)) end - @. y = dbase * x + s - for (p, q) in zedges - y[p] -= x[p] + x[q] - end - for (p, q) in vedges - y[p] += dκ * (x[p] + x[q]) - end - return y -end - -# `B*x` for the asymmetric Woodbury system -# `B = D + ridge·I + u_r·u_rᵀ + u_c·u_cᵀ − L_Z + dκ·L_V` on the stacked positions -# `1:m` (rows) and `m+1:m+n` (columns), with `drow`/`dcol` the two diagonal blocks of -# `D + ridge·I`. Every edge of `L_X` adds `x[p] + x[q]` at both of its ends. -function _woodbury_mul!(y, x, m, drow, dcol, zedges, vedges, dκ) - sr = zero(eltype(y)) - for p in 1:m - sr += x[p] - end - sc = zero(eltype(y)) - for p in (m+1):length(x) - sc += x[p] - end - for p in 1:m - y[p] = drow * x[p] + sr - end - for p in (m+1):length(x) - y[p] = dcol * x[p] + sc - end - for (p, q) in zedges - t = x[p] + x[q] - y[p] -= t - y[q] -= t - end - for (p, q) in vedges - t = dκ * (x[p] + x[q]) - y[p] += t - y[q] += t - end - return y + g = K \ (U' * y) + copyto!(x, y) + return mul!(x, Y, g, -1, 1) end # Jacobi-preconditioned conjugate gradients for the symmetric positive-definite @@ -713,6 +684,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), Cj = Int[] Cv = T[] rhs = zeros(T, use_woodbury ? n : 0, 2) + Umat = ones(T, use_woodbury ? n : 0, 1) # the gauge `e`, as the low-rank block cgx = zeros(T, use_woodbury ? n : 0) cgr = zeros(T, use_woodbury ? n : 0) cgz = zeros(T, use_woodbury ? n : 0) @@ -857,7 +829,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # `B = C + e·eᵀ` with `C = n·I − L_Z + (κ−1)·L_V`: the complete-support # matrix, corrected by the zero set `Z` and by the currently violated # entries `V`. One O(nnz) sweep collects the right-hand side, the violated - # set, and the diagonal of `B`; everything after it is O(n + |Z| + |V|). + # set, and the diagonal of `B`; everything after it is O(|Z| + |V|). dκ = κ === nothing ? zero(T) : T(κ) - oneunit(T) fill!(f, zero(T)) copyto!(dg, czero) @@ -888,24 +860,10 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), for p in 1:n dg[p] += oneunit(T) + ridge end - # Gershgorin on `(κ−1)·L_V` against a diagonal of at least `n` estimates - # the condition number of `B`. While that estimate is small the structured - # matvec plus conjugate gradients reaches the same answer in a few hundred - # O(n + |Z| + |V|) iterations, which is far cheaper than a factorization - # whose fill, on the near-random violated pattern of the early stages, is - # close to dense. - κest = oneunit(T) + dκ * 2 * maxdegV / n - if κest <= WOODBURY_CG_KAPPA - copyto!(cgx, α) - dbase = T(n) + ridge - Bmul! = (y, x) -> _symwoodbury_mul!(y, x, dbase, zedges, vedges, dκ) - it, ok = _pcg!(Bmul!, cgx, dg, f, cgr, cgz, cgd, cgAd, - 50 + 20 * ceil(Int, sqrt(κest)), 100 * eps(T) * norm(f)) - ncg[] += it - ok && return copy(cgx) - end - # `sparse` sums the duplicate triplets; the ridge rides on `C`'s diagonal. - nchol[] += 1 + # `zedges` and `vedges` carry both orientations of every pair, so these + # triplets store `C` in full rather than in one triangle: the same matrix + # then serves the matvec below and the factorization after it. `sparse` + # sums the duplicates, and the ridge rides on the diagonal. empty!(Ci) empty!(Cj) empty!(Cv) @@ -930,22 +888,29 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), push!(Cj, q) push!(Cv, dκ) end - F = cholesky(Symmetric(sparse(Ci, Cj, Cv, n, n))) - for p in 1:n - rhs[p, 1] = f[p] - rhs[p, 2] = oneunit(T) - end - # Sherman–Morrison: with y = C\f and u = C\e, (C + e·eᵀ)\f is - # y − u·(eᵀy)/(1 + eᵀu). - YU = F \ rhs - sy = zero(T) - su = zero(T) - for p in 1:n - sy += YU[p, 1] - su += YU[p, 2] + C = sparse(Ci, Cj, Cv, n, n) + # Gershgorin on `(κ−1)·L_V` against a diagonal of at least `n` estimates + # the condition number of `B`. While that estimate is small, conjugate + # gradients on `C·x + e·(eᵀx)` reach the same answer in a few hundred + # O(nnz(C)) applications, which is far cheaper than a factorization whose + # fill, on the near-random violated pattern of the early stages, is close + # to dense. + κest = oneunit(T) + dκ * 2 * maxdegV / n + if κest <= WOODBURY_CG_KAPPA + copyto!(cgx, α) + Bmul! = function (y, x) + mul!(y, C, x) + s = sum(x) + y .+= s + return y + end + it, ok = _pcg!(Bmul!, cgx, dg, f, cgr, cgz, cgd, cgAd, + 50 + 20 * ceil(Int, sqrt(κest)), 100 * eps(T) * norm(f)) + ncg[] += it + ok && return copy(cgx) end - r = sy / (oneunit(T) + su) - return [YU[p, 1] - r * YU[p, 2] for p in 1:n] + nchol[] += 1 + return _woodbury_solve!(zeros(T, n), cholesky(Symmetric(C)), Umat, f, rhs) else fill!(f, zero(T)) B = zeros(T, n, n) @@ -1199,6 +1164,16 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), Cj = Int[] Cv = T[] rhs = zeros(T, use_woodbury ? N : 0, 3) + # The row and column indicators, as the low-rank block. + Umat = zeros(T, use_woodbury ? N : 0, 2) + if use_woodbury + for ip in 1:m + Umat[ip, 1] = oneunit(T) + end + for jp in 1:n + Umat[m+jp, 2] = oneunit(T) + end + end cgx = zeros(T, use_woodbury ? N : 0) cgr = zeros(T, use_woodbury ? N : 0) cgz = zeros(T, use_woodbury ? N : 0) @@ -1340,8 +1315,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # `B + v0·v0ᵀ = C + U·Uᵀ` with `C = D − L_Z + (κ−1)·L_V`, # `D = diag(n·1_m, m·1_n)` and `U = [u_r u_c]` the row and column # indicators: the complete-support matrix, corrected by the zero set `Z` - # and by the currently violated entries `V`. `sparse` sums the duplicate - # triplets. + # and by the currently violated entries `V`. dκ = κ === nothing ? zero(T) : T(κ) - oneunit(T) fill!(f, zero(T)) copyto!(dg, czero) @@ -1381,19 +1355,10 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # O(N + |Z| + |V|) iterations, which is far cheaper than a factorization # whose fill on the near-random violated pattern of the early stages # approaches dense. - κest = oneunit(T) + dκ * 2 * maxdegV / min(m, n) - if κest <= WOODBURY_CG_KAPPA - copyto!(cgx, x) - drow = T(n) + ridge - dcol = T(m) + ridge - Bmul! = (yy, xx) -> _woodbury_mul!(yy, xx, m, drow, dcol, zedges, vedges, dκ) - it, ok = _pcg!(Bmul!, cgx, dg, f, cgr, cgz, cgd, cgAd, - 50 + 20 * ceil(Int, sqrt(κest)), 100 * eps(T) * norm(f)) - ncg[] += it - ok && return copy(cgx) - end - # `sparse` sums the duplicate triplets; the ridge rides on `C`'s diagonal. - nchol[] += 1 + # `zedges` and `vedges` carry both ends of every pair, so these triplets + # store `C` in full rather than in one triangle: the same matrix then + # serves the matvec below and the factorization after it. `sparse` sums + # the duplicates, and the ridge rides on the diagonal. empty!(Ci) empty!(Cj) empty!(Cv) @@ -1424,40 +1389,41 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), push!(Cj, p) push!(Cv, dκ) end - F = cholesky(Symmetric(sparse(Ci, Cj, Cv, N, N))) - fill!(rhs, zero(T)) - for p in 1:N - rhs[p, 1] = f[p] - end - for ip in 1:m - rhs[ip, 2] = oneunit(T) - end - for jp in 1:n - rhs[m+jp, 3] = oneunit(T) - end - # Woodbury with a 2x2 capacitance: with y = C\f and Y = C\U, - # (C + U·Uᵀ)\f is y − Y·((I₂ + UᵀY)\(Uᵀy)). - YU = F \ rhs - ty = zero(T) - tz = zero(T) - k11 = zero(T) - k12 = zero(T) - k21 = zero(T) - k22 = zero(T) - for ip in 1:m - ty += YU[ip, 1] - k11 += YU[ip, 2] - k12 += YU[ip, 3] - end - for jp in 1:n - q = m + jp - tz += YU[q, 1] - k21 += YU[q, 2] - k22 += YU[q, 3] + C = sparse(Ci, Cj, Cv, N, N) + # Gershgorin on `(κ−1)·L_V` against the smaller diagonal block estimates + # the condition number of `B`. While that estimate is small, conjugate + # gradients on `C·x + u_r·(u_rᵀx) + u_c·(u_cᵀx)` reach the same answer in a + # few hundred O(nnz(C)) applications, which is far cheaper than a + # factorization whose fill on the near-random violated pattern of the early + # stages approaches dense. + κest = oneunit(T) + dκ * 2 * maxdegV / min(m, n) + if κest <= WOODBURY_CG_KAPPA + copyto!(cgx, x) + Bmul! = function (yy, xx) + mul!(yy, C, xx) + sr = zero(T) + for p in 1:m + sr += xx[p] + end + sc = zero(T) + for p in (m+1):N + sc += xx[p] + end + for p in 1:m + yy[p] += sr + end + for p in (m+1):N + yy[p] += sc + end + return yy + end + it, ok = _pcg!(Bmul!, cgx, dg, f, cgr, cgz, cgd, cgAd, + 50 + 20 * ceil(Int, sqrt(κest)), 100 * eps(T) * norm(f)) + ncg[] += it + ok && return copy(cgx) end - K = [oneunit(T)+k11 k12; k21 oneunit(T)+k22] - g = K \ T[ty, tz] - return [YU[p, 1] - g[1] * YU[p, 2] - g[2] * YU[p, 3] for p in 1:N] + nchol[] += 1 + return _woodbury_solve!(zeros(T, N), cholesky(Symmetric(C)), Umat, f, rhs) else fill!(f, zero(T)) B = v0 * v0'