Improve the performance of cover_transversal

docs/src/objectives.md

@@ -87,9 +87,9 @@ The method used to minimize the objective depends on the penalty:
   for hard covers, see [Objective-minimal covers](@ref).
 - `PowerMean{p}` with a target ratio, for [`cover_transversal`](@ref) (native),
   uses an active-set Newton method that keeps every iterate a hard cover
-  tight on the transversal. With `maxbalance=true`, max-plus policy iteration
-  max-balances each block first, and the Newton method acts on one scale per
-  block.
+  tight on the transversal. With `maxbalance=true`, a parametric shortest-path
+  algorithm max-balances each block first, and the Newton method acts on one
+  scale per block.
 - `AbsLog{1}` is a linear program solved with JuMP and HiGHS.
 - `AbsLinear` is not convex; several starts are each refined to a local minimum
   with JuMP and Ipopt, and the best is returned.

src/active_set.jl

@@ -386,23 +386,50 @@ end
 function _laplacian_factor(L::SparseMatrixCSC{T,Int}; fillbudget::Real=LSQR_FILL_BUDGET,
                            flopbudget::Real=LSQR_FLOP_BUDGET, cgmaxiter::Int=10 * size(L, 1) + 100) where {T}
     size(L, 1) == 0 && return _laplacian_factor(Matrix(L))
+    return _laplacian_cg(L; cgmaxiter)
+end
+# Conjugate-gradient solver for `L`, built without a symbolic analysis, for a
+# pattern presumed to exceed the budgets. The fallback factorization is analyzed
+# only if CG fails (see `_cg_fallback`).
+function _laplacian_cg(L::SparseMatrixCSC{T,Int}; cgmaxiter::Int=10 * size(L, 1) + 100) where {T}
     S = L + triu(L, 1)'
     return _LaplacianCG{T}(L, S, Vector(diag(S)), nothing, Ref{Any}(nothing), cgmaxiter, Ref(zero(T)))
 end
+# The factorization that replaces conjugate gradients once they fail. In
+# `Float64` it reuses the symbolic analysis in `C.F` when there is one.
+function _cg_fallback(C::_LaplacianCG{Float64})
+    F = C.F
+    if F === nothing
+        F = SparseCholesky()
+        analyze!(F, C.L)
+    end
+    return factorize!(F, C.L)
+end
+_cg_fallback(C::_LaplacianCG) = _laplacian_factor(Matrix(C.L))
 _laplacian_factor(L::Matrix) = LinearAlgebra.cholesky!(Symmetric(L, :U))
 
-_laplacian_solve!(x::Vector{Float64}, F::SparseCholesky, b::Vector{Float64}) =
+# Solves `L x = b`. The keywords apply only to a `_LaplacianCG`: `rtol` is the
+# relative residual at which conjugate gradients stop, and `refine` enables the
+# iterative refinement that follows, which reduces the forward error to the
+# resolution of the arithmetic. `_solve_error` then estimates the remaining
+# error; it is `NaN` when refinement is skipped.
+_laplacian_solve!(x::Vector{Float64}, F::SparseCholesky, b::Vector{Float64}; rtol=nothing, refine=true) =
     isempty(b) ? x : solve!(x, F, CHOLMOD_A, b)
-_laplacian_solve!(x, F::Cholesky, b) = ldiv!(x, F, b)
+_laplacian_solve!(x, F::Cholesky, b; rtol=nothing, refine=true) = ldiv!(x, F, b)
 
-function _laplacian_solve!(x::Vector{T}, C::_LaplacianCG{T}, b::Vector{T}) where {T}
+function _laplacian_solve!(x::Vector{T}, C::_LaplacianCG{T}, b::Vector{T}; rtol::Real=sqrt(eps(T)),
+                           refine::Bool=true) where {T}
     if C.fallback[] === nothing
         n = length(b)
         fill!(x, zero(T))
         rdd, r, z, d, Ad, e = ntuple(_ -> zeros(T, n), 6)
         Smul! = (y, v) -> mul!(y, C.S, v)
-        η = sqrt(eps(T))
+        η = T(rtol)
         _, ok = _pcg!(Smul!, x, C.dg, b, r, z, d, Ad, C.maxiter, η * norm(b))
+        if ok && !refine
+            C.err[] = T(NaN)
+            return x
+        end
         # Iterative refinement with the residual in compensated arithmetic.
         # CG leaves a forward error of its normwise backward error times the
         # condition number, and the multipliers, routed sums of the solution
@@ -426,7 +453,7 @@ function _laplacian_solve!(x::Vector{T}, C::_LaplacianCG{T}, b::Vector{T}) where
             δ <= δprev / 100 || break
             δprev = δ
         end
-        C.fallback[] = C.F === nothing ? _laplacian_factor(Matrix(C.L)) : factorize!(C.F, C.L)
+        C.fallback[] = _cg_fallback(C)
     end
     C.err[] = zero(T)
     return _laplacian_solve!(x, C.fallback[], b)

src/balanced.jl

@@ -15,6 +15,11 @@ const BALANCED_NOISE = 64
 # Newton systems on at most this many unknowns, or at least an eighth full, are
 # factored densely.
 const BALANCED_DENSE_MAX = 128
+# Relative residual of a conjugate-gradient solve for a Newton direction. Every
+# iterate of conjugate gradients from zero is a descent direction, and the
+# reduced gradient is recomputed at each step, so the direction need not be
+# accurate.
+const BALANCED_CG_RTOL = 1e-4
 # Ridge retries for a Newton matrix that roundoff makes indefinite; each retry
 # scales the ridge by `BALANCED_RIDGE_GROWTH`. The ridge is relative to each
 # diagonal entry, so that it does not swamp the curvature of a group whose
@@ -93,12 +98,17 @@ descent. Where the groups are stationary, a maximum-flow problem on each group's
 active entries either proves the KKT conditions or yields a set of rows whose
 scales can decrease, releasing the active entries that enter it.
 
-With `maxbalance=true`, each block is max-balanced by Algorithm 3.10 of Hook et
-al.: max-plus policy iteration finds a cycle of largest geometric-mean ratio in
-the block, a diagonal similarity brings every ratio of the block down to that
-mean, and the cycle's rows, whose relative scales are then fixed, are merged
-into one; this repeats until one row remains. The active-set Newton method then
-runs with one scale per block.
+With `maxbalance=true`, each block is max-balanced as in Algorithm 3.10 of Hook
+et al.: a cycle of largest geometric-mean ratio in the block is found, a
+diagonal similarity brings every ratio of the block down to that mean, and the
+cycle's rows, whose relative scales are then fixed, are merged into one; this
+repeats until one row remains. The cycles are found by the parametric
+shortest-path algorithm of Young, Tarjan & Orlin (Networks 21:205–221, 1991).
+For a block with `n` rows and `m` entries, its cost is at most
+`O(nm + n² log n)`. The cost is roughly linear in `m` when the block's
+shortest-path tree stays shallow, as it does for most sparse matrices, and
+about `O(nm)` for a banded block, whose tree is a path. The active-set Newton
+method then runs with one scale per block.
 
 `tol` (default `1e-9`) is the stopping tolerance: the groups are stationary when
 the derivative of the objective with respect to each group's scale is at most
@@ -110,7 +120,8 @@ nonnegative but within the roundoff of the groups' derivatives. `maxiter`
 the active set.
 `fillbudget` and `flopbudget` are as for [`cover_min`](@ref): a Newton system
 whose sparse Cholesky factor would exceed them is solved by conjugate
-gradients. The Newton systems have condition numbers of order `target^(-p)`
+gradients to a loose tolerance, since a descent direction need not be
+accurate. The Newton systems have condition numbers of order `target^(-p)`
 times that of the graph, and of order `r^(-p)` where cycles of entries hold
 ratios `r` far below `target`, so a very small `target^p` or a large `p`
 degrades the linear solves and the convergence.
@@ -663,7 +674,7 @@ function _balanced_direction!(δv, cache, d, gid, uidx, nu, gG, P, Q, p::T, tp::
             L = sparse([I; 1:nu], [J; 1:nu], [V; dg .+ ridge .* max.(dg, dfloor)], nu, nu)
             F = _balanced_factor!(cache, L, fillbudget, flopbudget)
             if F !== nothing
-                _laplacian_solve!(x, F, rhs)
+                _laplacian_solve!(x, F, rhs; rtol=BALANCED_CG_RTOL, refine=false)
                 kind = F isa SparseCholesky ? :cholesky : :cg
                 break
             end
@@ -682,17 +693,27 @@ end
 # made indefinite.
 _balanced_ridge(ridge::T) where {T} = iszero(ridge) ? eps(T) : ridge * T(BALANCED_RIDGE_GROWTH)
 
-# Factor the pinned Laplacian `L` (upper triangle), reusing the symbolic analysis
-# in `cache` when the pattern is unchanged. Returns `nothing` when `L` is not
-# numerically positive definite.
+# Factor the pinned Laplacian `L` (upper triangle), or return a
+# conjugate-gradient solver for it when its factor would exceed the budgets.
+# Returns `nothing` when `L` is not numerically positive definite. `cache` holds
+# one of:
+# - the analyzed factorization and its pattern, reused when the pattern of `L`
+#   is unchanged;
+# - the order and number of nonzeros of the last pattern whose factor exceeded
+#   the budgets. A pattern more than half that size in both is presumed to
+#   exceed them too and gets conjugate gradients without a symbolic analysis,
+#   which would cost more than the solve.
 function _balanced_factor!(cache, L::SparseMatrixCSC{Float64,Int}, fillbudget, flopbudget)
     c = cache[]
     try
-        if c !== nothing && c[2] == SparseArrays.getcolptr(L) && c[3] == rowvals(L)
-            return factorize!(c[1], L)
+        if c isa Tuple{SparseCholesky,Vector{Int},Vector{Int}}
+            c[2] == SparseArrays.getcolptr(L) && c[3] == rowvals(L) && return factorize!(c[1], L)
+        elseif c isa Tuple{Int,Int}
+            2 * size(L, 1) > c[1] && 2 * nnz(L) > c[2] && return _laplacian_cg(L)
         end
         F = _laplacian_factor(L; fillbudget, flopbudget)
-        cache[] = F isa SparseCholesky ? (F, copy(SparseArrays.getcolptr(L)), copy(rowvals(L))) : nothing
+        cache[] = F isa SparseCholesky ? (F, copy(SparseArrays.getcolptr(L)), copy(rowvals(L))) :
+                  (size(L, 1), nnz(L))
         return F
     catch err
         err isa LinearAlgebra.PosDefException || rethrow()