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..056a566 100644 --- a/src/MatrixCovers.jl +++ b/src/MatrixCovers.jl @@ -1,9 +1,11 @@ module MatrixCovers using LinearAlgebra: LinearAlgebra, Adjoint, Bidiagonal, Diagonal, Hermitian, - SymTridiagonal, Symmetric, Transpose, Tridiagonal, dot, norm + SymTridiagonal, Symmetric, Transpose, Tridiagonal, cholesky, + dot, mul!, norm using PrecompileTools: PrecompileTools, @compile_workload using Random: Random, AbstractRNG, MersenneTwister +using SparseArrays: SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals, sparse export AbsLog, AbsLinear export cover_objective, iscover @@ -27,6 +29,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..e9af960 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -18,13 +18,26 @@ 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 - (from the SparseArrays extension) 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 `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` + 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 @@ -58,12 +71,26 @@ 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 (from the SparseArrays - extension) 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 `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 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` 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 @@ -286,12 +313,140 @@ 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 `Z` of `A` together with the currently +# violated entries `V`) and symmetric positive definite, and `U` has one column +# (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 +# λ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 + +# `(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 + g = K \ (U' * y) + copyto!(x, y) + return mul!(x, Y, g, -1, 1) +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. +# +# `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) + 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, fresh && nrm <= tol +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`. @@ -301,6 +456,10 @@ end # strength κ) rather than that of `MᵀM` (≈ κ); at κ = 1e8 the squared conditioning # breaks CG while LSQR stays accurate. # +# 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 # ‖Mᵀ(b - Mx)‖ ≤ atol · ‖M‖ · ‖b - Mx‖, both estimated from the bidiagonalization @@ -359,20 +518,21 @@ 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. -# `linsolve` is `:auto`/`:dense` (dense factorization) or `:lsqr` -# (matrix-free, for sparse supports). `start`, when given, is a positive cover of `A` +# 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 # 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")) - # 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. + 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) ax = axes(A, 1) axes(A, 2) == ax || throw(ArgumentError("symcover_min requires a square matrix")) @@ -382,6 +542,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 @@ -391,14 +554,62 @@ 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`. 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 && nzero <= zbudget + if linsolve === :woodbury && !ok + T === Float64 || + throw(ArgumentError("linsolve=:woodbury requires Float64 arithmetic, but `A` works in $T; use :dense or :lsqr")) + 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 + # 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}[] + # 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) + for s in _slots(G, i) + mark[G.idx[s] - first(ax) + 1] = true + end + for jp in 1:n + mark[jp] && continue + push!(zedges, (ip, jp)) + czero[ip] -= oneunit(T) + ip == jp && (czero[ip] -= oneunit(T)) + end + fill!(mark, false) + end + end fκ = function (α, κ) v = zero(T) for (e, (ip, jp)) in enumerate(edges) @@ -407,29 +618,194 @@ 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 # 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) + # 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, 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 + # 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 + # `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, 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) + 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 + 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 + 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) + for p in 1:n + push!(Mi, p) + push!(Mj, p) + push!(Mv, dpart[p]) + 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) + 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)::Vector{T} end Amul! = function (y, x) for (e, (ip, jp)) in enumerate(edges) @@ -449,12 +825,100 @@ 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`. One O(nnz) sweep collects the right-hand side, the violated + # 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) + fill!(degV, 0) + empty!(vedges) + for (e, (ip, jp)) in enumerate(edges) + c = cvals[e] + viol = κ !== nothing && (α[ip] + α[jp] - c) < 0 + vpat[e] = viol + f[ip] += (viol ? T(κ) : oneunit(T)) * c + if viol + 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, 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 + # `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) + for p in 1:n + push!(Ci, p) + push!(Cj, p) + 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 + 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 + nchol[] += 1 + return _woodbury_solve!(zeros(T, n), cholesky(Symmetric(C)), Umat, f, rhs) else fill!(f, zero(T)) 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 @@ -479,12 +943,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 && stable && break fcur - fnew <= 5000 * eps(T) * max(fcur, one(T)) && break fcur = fnew end @@ -503,18 +975,21 @@ 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[], 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, 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), 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 @@ -525,6 +1000,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 @@ -546,6 +1024,66 @@ 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 + # 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 && nzero <= zbudget + if linsolve === :woodbury && !ok + T === Float64 || + throw(ArgumentError("linsolve=:woodbury requires Float64 arithmetic, but `A` works in $T; use :dense or :lsqr")) + 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 + # 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] && continue + q = m + jp + push!(zedges, (ip, q)) + czero[ip] -= oneunit(T) + czero[q] -= oneunit(T) + 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,12 +1101,29 @@ 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; 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, @@ -576,19 +1131,166 @@ 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) + 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 + # 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 + # `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, 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) + 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 + 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 + 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) + for p in 1:N + push!(Mi, p) + push!(Mj, p) + push!(Mv, dpart[p]) + 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) + 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)::Vector{T} + end Amul! = function (y, xx) for (e, (p, q)) in enumerate(edges) y[e] = ws[e] * (xx[p] + xx[q]) @@ -609,12 +1311,127 @@ 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`. + dκ = κ === nothing ? zero(T) : T(κ) - oneunit(T) + fill!(f, zero(T)) + 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 + vpat[e] = viol + w = viol ? T(κ) : oneunit(T) + f[p] += w * c + f[q] += w * c + if viol + 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, 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 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. + # `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) + 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, dκ) + end + 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 + nchol[] += 1 + return _woodbury_solve!(zeros(T, N), cholesky(Symmetric(C)), Umat, f, rhs) else fill!(f, zero(T)) 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 @@ -663,12 +1480,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 && stable && break fcur - fnew <= 5000 * eps(T) * max(fcur, one(T)) && break fcur = fnew end @@ -719,7 +1544,9 @@ 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[], cgiters=ncg[], + cholsolves=nchol[], + linsolve=(use_lsqr ? :lsqr : use_woodbury ? :woodbury : :dense)) end # Workers for the soft (unconstrained) AbsLog{2} covers. The soft objective diff --git a/ext/MatrixCoversSparseArraysExt.jl b/src/sparse_support.jl similarity index 69% rename from ext/MatrixCoversSparseArraysExt.jl rename to src/sparse_support.jl index 30c1452..5807b69 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 @@ -85,46 +80,45 @@ 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. -function MatrixCovers.symcover_min(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) +# 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 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/minimal_covers.jl b/test/minimal_covers.jl index 4e53aa2..7f9eb99 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,214 @@ 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) + # 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) + 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) + 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) + # 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 sw.cholsolves > 0 + @test sd.cgiters == 0 + @test sd.cholsolves == 0 + 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) + @test sw.cgiters > 0 + @test sw.cholsolves > 0 + @test sd.cgiters == 0 + @test sd.cholsolves == 0 + 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 + + # 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 + + # 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 + @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 + + # 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 + 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 + +# 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. 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 + + 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)) + # 22 solves measured here against 26 on `:lsqr`; 24 leaves two of headroom. + @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 + # 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. + 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 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/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 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)