Reference

AbsLinear{p}

Penalty type for φ(r) = |1 - r|^p. Unlike AbsLog, this penalty is continuous at r = 0, so each zero entry of A contributes φ(0) = 1 to the objective.

The resulting optimization problems are not convex and may have multiple local minima.

AbsLog{p}

Penalty type for

φ(r) = |log(r)|^p  if r > 0
       0           if r = 0

Defining φ(0) = 0 makes zero entries of A, for which log(r) = -∞, contribute nothing to the objective. The objective is convex in the logarithms of the scales; AbsLog{1} may have multiple minima.

See also: AbsLinear.

AbstractCoverPenalty <: Function

Supertype of cover penalties. Built-in subtypes are AbsLog, AbsLinear, and PowerMean.

cover_objective applies the penalty to r = |A[i,j]|/(a[i]*b[j]) and sums over A.

Extending

A subtype must be callable on a nonnegative real:

(::MyPenalty)(r::Real)

The method must accept r = 0 and r = typemax(...). Penalties are usually singleton structs.

cover_objective works for any subtype, but solvers support only specific built-in penalties:

• hard covers (symcover_min, cover_min): AbsLog{2} natively, AbsLog{1} and AbsLinear through JuMP;
• soft covers (soft_symcover, soft_cover): PowerMean and AbsLog{2} natively, AbsLog{1} and AbsLinear through JuMP.

Passing a custom subtype to a solver raises a MethodError.

AbstractCoverStats

Supertype of the statistics that the optimizing cover solvers return as their last output, e.g. a, stats = symcover_min(A) or a, b, stats = cover_min(A). Each solver family has its own subtype, documented with its fields.

A subtype implements those of converged, iterations, residual, and tolerance that apply to its solver. A solver without a stopping criterion implements none of them, so calling one throws a MethodError.

ActiveSetStats{T}

Result of the active-set method that finishes an AbsLog{2} minimal-cover solve whose multiplier iteration stopped before converging.

Fields

• certified::Bool: the KKT conditions hold to the tolerances below, so the result is the minimizer to rounding error. When false, the solver returns the multiplier iterate instead.
• nsteps::Int: active-set steps, each a graph-Laplacian solve.
• primal::T, ptol::T: the largest constraint violation of the result, in the log domain, and its tolerance.
• dual::T, dtol::T: the magnitude of the most negative multiplier of an active constraint, and its tolerance.

primal and dual are NaN when certified is false.

AugmentedLagrangianStats{T} <: AbstractCoverStats

Statistics of the native AbsLog{2} minimal-cover solver, returned by symcover_min, cover_min, cover_transversal, and their mutating forms.

The multiplier iteration stops when its KKT residual reaches tol. If it stops earlier, the active-set method described by polish finishes the solve. The accessors describe the returned cover: converged is true when either stage met its criterion, and residual and tolerance are those of the stage that produced the result (polish.primal and polish.ptol after a certified finish, otherwise last(kkt) and tol). iterations counts multiplier updates.

Fields

• converged::Bool, tol::T: as above.
• nouter::Int: multiplier updates.
• kkt::Vector{T}, κs::Vector{T}, exits::Vector{Symbol}, drops::Vector{T}: per update, the KKT residual after the update, the penalty weight, how the inner minimization ended, and its last relative decrease of the objective.
• polish::Union{Nothing,ActiveSetStats{T}}: the active-set finish, or nothing if it did not run.
• linsolve::Symbol: the inner solver (:dense, :woodbury, or :lsqr).
• The remaining fields count the inner linear-algebra work and describe the preconditioner; their names and meanings may change between releases: nsolves, lsqriters, lsqrtrace, cgiters, cholsolves, precond, nrefactor, nforest, fill_entries, factor_flops, nzeroed.

ExternalSolverStats <: AbstractCoverStats

Statistics of the covers computed by an external optimizer through JuMP: HiGHS for AbsLog{1} (symcover_min, cover_min, soft_symcover, soft_cover, and their mutating forms), and Ipopt for one local refinement under AbsLinear (symcover_min!, cover_min!, soft_symcover!, and soft_cover!).

A solve may consist of several stages, each an optimization of one model. Ipopt runs one. HiGHS runs two unless A has no nonzero entries: a linear program for the AbsLog{1} objective, then a quadratic program that selects the point of least AbsLog{2} objective among its minimizers. The fields hold one entry per stage, in order.

A stage that does not end with a solved status throws a SolverFailure, so every status recorded here is a solved one and converged returns true. The accessors iterations, residual, and tolerance are not implemented.

Fields

• solver::Symbol: :HiGHS or :Ipopt.
• status::Vector{Symbol}: the termination status of each stage, as the name of the MathOptInterface.TerminationStatusCode (:OPTIMAL or :LOCALLY_SOLVED).
• objective::Vector{Float64}: the optimal value of each stage's model. The models work in log scales and may include slack variables or constant terms, so this need not equal cover_objective of the returned cover.
• niters::Vector{Int}: the iterations of each stage as the solver reports them through MathOptInterface: for HiGHS, its simplex and interior-point iterations (its quadratic-program iterations are not reported, so the second stage counts 0); for Ipopt, its interior-point iterations.
• solvetime::Vector{Float64}: the solver's run time of each stage, in seconds.

niters and solvetime may change meaning between releases.

LeastSquaresStats <: AbstractCoverStats

Statistics of the soft AbsLog{2} covers (soft_symcover and soft_cover with AbsLog{2}()), which are computed by one linear least-squares solve. There is no stopping criterion, so none of the accessors of AbstractCoverStats apply.

Fields

• linsolve::Symbol: the solver (:dense, :woodbury, or :lsqr).
• The remaining fields count the linear-algebra work and describe the preconditioner, as for AugmentedLagrangianStats; their names and meanings may change between releases: lsqriters, cgiters, cholsolves, precond, fill_entries, factor_flops.

MultistartStats{S,T} <: AbstractCoverStats

Statistics of a cover selected from several local refinements: the AbsLinear forms of symcover_min and cover_min, which refine one start per entry of strategies, and of soft_symcover and soft_cover, which refine deterministic and random starts. The cover with the least objective is returned.

A failed refinement of a random start (one that throws SolverFailure) is left out of the selection and recorded in failed; any other failure throws. The accessors converged, iterations, residual, and tolerance return those of the selected start, stats[selected], where implemented for S.

Fields

• labels::Vector{Symbol}: the starts whose refinement succeeded, in the order they ran: the strategy names, or for the soft covers names such as :geomean and :rand1.
• stats::Vector{S}: the statistics of each of these refinements.
• objectives::Vector{T}: the cover_objective of each refined cover.
• selected::Int: the index (into the vectors above) of the returned cover.
• failed::Vector{Symbol}: the labels of the random starts whose refinement failed.

PowerMean{p}

Penalty type for

φ(r) = (r^p - 1 - p*log(r)) / p  if r > 0
       0                         if r = 0

with p > 0. As for AbsLog, defining φ(0) = 0 makes zero entries of A contribute nothing.

PowerMean{2} is the default penalty for soft covers. For entries modeled as independent normal variables with standard deviations a[i]*b[j], it gives the maximum-likelihood scales.

The objective is strictly convex in the log scales, so the soft cover it defines has unique products a[i]*b[j] on the support of A, and a symmetric A has a symmetric minimizer. At the minimum, the p-power mean of the ratios r = |A[i,j]|/(a[i]*b[j]) over the nonzero entries of each row and each column equals one: ∑_j r[i,j]^p = n_i and ∑_i r[i,j]^p = m_j, where n_i and m_j count the nonzeros of row i and column j. Consequently no ratio exceeds min(n_i, m_j)^(1/p).

See also: soft_cover, soft_symcover.

PowerMeanStats{T} <: AbstractCoverStats

Statistics of the PowerMean soft covers (soft_symcover, soft_cover, and their mutating forms with a PowerMean penalty or no penalty).

The imbalance of a nonzero row (or column) is |mean(r^p) - 1|, the mean taken over the ratios r of its nonzero entries. The solve stops when the largest imbalance is at most tol. Power-mean updates run first; when they are slow, Newton steps finish the solve. converged, residual, and tolerance return the fields of the same names (residual returns imbalance); iterations returns nsweeps + nnewton.

Fields

• converged::Bool: imbalance <= tol.
• tol::T: the stopping tolerance.
• imbalance::T: the largest imbalance of the returned cover.
• nsweeps::Int: power-mean updates. For the asymmetric solver each counts an update of all rows and then all columns (a sweep); for the symmetric solver, one simultaneous update of all scales.
• newton::Bool: whether Newton steps ran.
• nnewton::Int: Newton steps.
• The remaining fields count the work of the Newton stage; their names and meanings may change between releases: npasses (passes over the support), nfactor (factorizations), and linsolve (the solver for the Newton equations: :dense, :cholesky, or :cg, and :none when Newton did not run).

SolverFailure(msg)

Thrown by the AbsLog{1} and AbsLinear forms of symcover_min, cover_min, soft_symcover, soft_cover, and their mutating forms when the external solver (HiGHS or Ipopt) stops without reaching a usable local minimum: it reports a status other than solved, or it returns scales that diverge on a plateau of the objective. A different start may succeed. The AbsLinear multistarts of soft_cover and soft_symcover leave out perturbed starts that throw this and record them in MultistartStats.

The native solvers do not throw this: one that stops before meeting its stopping criterion warns and returns a result whose statistics report converged(stats) == false.

SupportComponents

Connected components of a matrix's bipartite support graph, as returned by support_components: one vertex per row and one per column, one edge per stored nonzero.

Component ids run 1:ncomponents(sc); unsupported rows and columns report 0. Use rowcomponent and colcomponent with the matrix's own indices.

colcomponent(sc::SupportComponents, j) -> Int

Component id of column j, or 0 if that column has empty support. j is the matrix's own column index.

converged(stats::AbstractCoverStats) -> Bool

Report whether the returned cover meets the solver's stopping criterion.

a, b = cover!(ϕ, a, b, A; maxiter, cgiter, start=:auto)
a, b = cover!(a, b, A; maxiter, cgiter, start=:auto)

Mutating counterpart of cover: writes the hard cover into a and b and returns them, rather than allocating new vectors. eachindex(a) must match axes(A, 1) and eachindex(b) must match axes(A, 2). ϕ has the same meaning as in cover, and is likewise ignored by the current heuristic covers.

See also: cover.

a, b = cover(A, a::AbstractVector)

Return the tightest cover (a, b) of A with the row scale a held fixed: b[j] = max_i abs(A[i,j]) / a[i], the column maxima of Diagonal(1 ./ a) * A. The returned a is the argument itself. Columns with no stored nonzero get b[j] = 0.

a must be positive on every row that has a stored nonzero and nonnegative elsewhere; eachindex(a) must equal axes(A, 1). For each supported column, some i attains a[i]*b[j] ≈ abs(A[i,j]) (up to one rounding), and a[i]*b[j] >= abs(A[i,j]) holds exactly in floating point.

See also: cover, gramcover, iscover.

Examples

julia> A = [1.0 2.0; 4.0 1.0];

julia> a, b = cover(A, [1.0, 2.0]);

julia> b
2-element Vector{Float64}:
 2.0
 2.0

julia> iscover(a, b, A)
true

a, b = cover(ϕ, A; maxiter, cgiter, start=:auto)
a, b = cover(A; maxiter, cgiter, start=:auto)

Given a matrix A, return vectors a and b such that a[i] * b[j] >= abs(A[i, j]) for all i, j. This is a fast heuristic; for the cover that minimizes a penalty ϕ, use cover_min.

The heuristic refines a starting point with up to cgiter conjugate-gradient iterations (cgiter=0 disables refinement), raises the scales until every entry is covered, and then applies maxiter improvement passes.

start selects the kind of cover and the point the refinement begins from:

• :transversal returns a cover that is tight on every maximum-product transversal of the square matrix A (see the extended help), a heuristic counterpart of cover_transversal. A must be square and structurally nonsingular. Defaults: cgiter=20, maxiter=5.
• :covariant ensures the result is scale-covariant (typically at the cost of permutation-equivariance). Defaults: cgiter=16, maxiter=3.
• :geomean makes the result permutation-equivariant up to floating-point round-off (typically at the cost of scale-covariance). Defaults as for :covariant.
• :auto (the default) uses :transversal when A is square and structurally nonsingular, and :covariant otherwise.

The factors use the per-component balance convention described by cover_min. ϕ is accepted for API compatibility but is ignored.

See also: cover!, cover_transversal, cover_min, symcover.

Examples

julia> A = [1 2 3; 6 5 4];

julia> a, b = cover(A)
([1.2605791487949616, 3.475905976751536], [1.7261686708847903, 1.6168981055127307, 2.3798584982725366])

julia> a * b'
2×3 Matrix{Float64}:
 2.17597  2.03823  3.0
 6.0      5.62019  8.27216

Extended help

For square A, every hard cover has prod(a) * prod(b) >= π*, where π* is the largest product ∏ᵢ |A[i,σ(i)]| over permutations σ. Equality holds exactly when the cover is tight on every entry of every transversal attaining π*, and these covers maximize |det(A ./ (a .* b'))|. A cover that leaves such a transversal loose can make the scaled matrix far more ill-conditioned than necessary. The :transversal start prevents this; the :covariant start does not.

The :transversal start first finds a maximum-product transversal by shortest augmenting paths. Tightness on the maximum-product transversals fixes every row and column scale up to one unknown per group of rows and columns joined by their entries. These unknowns are then fitted to the remaining entries by a layer-average start on the graph of groups, cgiter conjugate-gradient iterations on the least-squares fit, a feasibility boost, and maxiter projected Gauss–Seidel sweeps that keep every entry covered. The transversal search typically dominates the cost.

With start=:covariant, initialization and refinement approximate the least-squares fit of log(a[i]) + log(b[j]) to log(abs(A[i, j])) over nonzero entries; coverage is then restored by boosting and maxiter tightening passes. The initialization alone gives the exact least-squares fit when each connected component of the support is a tree or is complete bipartite, so complete support needs no refinement. Each refinement iteration makes one pass over the support. The :covariant and :geomean starts coincide on complete support, where the geometric mean already solves the least-squares problem.

With start=:transversal or :covariant the cover products are scale-covariant on the support: for positive diagonal D1, D2, covering D1 * A * D2 multiplies each supported product by D1[i, i] * D2[j, j]. This holds in exact arithmetic for every cgiter and maxiter. In floating point, the products are covariant to within the roundoff accumulated by the conjugate-gradient iterations, typically near machine precision. cgiter caps the iterations, which stop earlier once the residual is at roundoff.

With start=:geomean the cover is permutation-equivariant up to floating-point round-off: covering A[p, q] for permutations p, q gives a[p] and b[q].

a, b, stats = cover_min(ϕ, A; kwargs...)
a, b, stats = cover_min(A; kwargs...)

Return the ϕ-minimal asymmetric hard cover of A: vectors a, b minimizing ∑ ϕ(|A[i,j]|/(a[i]*b[j])) subject to a[i]*b[j] >= |A[i,j]|. The default penalty is AbsLog{2}().

The objective and the constraints depend only on the products a[i]*b[j], which do not change when a is multiplied by c and b divided by c on one connected component of the support of A. The returned factors are fixed by the balance convention ∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j] within each component, where nzaᵢ and nzbⱼ are the nonzero counts of row i and column j.

The factors of a nonsymmetric band matrix can vary geometrically along the band. For the tridiagonal Toeplitz matrix with subdiagonal s and superdiagonal u, the AbsLog{2} minimizer has a[i+1]/a[i] ≈ sqrt(s/u), so for n = 1000 and s/u = 1/2 the entries of a span about 150 decades. The products a[i]*b[j] are unaffected, but quantities that combine factors from distant rows or columns, such as norm(a), inherit the spread. cover and soft_cover behave the same way. Factors outside the floating-point range throw an error; converting A to BigFloat extends the range.

Supported ϕ values, stats, and the failure behavior are as for symcover_min:

• AbsLog{2}(): native.
• AbsLog{1}(): requires JuMP and HiGHS.
• AbsLinear{1}() and AbsLinear{2}(): require JuMP and Ipopt and return the best local minimum found from strategies, a collection of initialize_cover strategy names (default (:hardcover, :covariant)).

The AbsLog objectives are convex in the logarithms of the scales, so their local minima are global. When the AbsLog{1} minimum is not unique, the method returns the minimizer with the smallest AbsLog{2} objective.

Extended help

The native solver accepts the same κ, maxouter, maxiter, tol, fillbudget, flopbudget, and linsolve keywords as symcover_min. For :woodbury, an m × n matrix may omit at most min(m,n) ÷ 4 entries per row or column and 4 * max(m,n) entries in total. :dense costs O((m+n)³) per Newton step; sparse matrices default to :lsqr. The accuracy and scale covariance of the result are as described for symcover_min. For A ./ (d * e'), the products a[i]*b[j] approximate those of cover_min(A) divided by d[i]*e[j].

When the multiplier iteration stops before its KKT residual reaches the tolerance, an active-set method (the dual method of Goldfarb and Idnani) continues from its result. The entries that are tight at the solution form a graph. Each step minimizes the objective with the entries of a spanning forest of this graph held tight, which requires solving a sparse graph-Laplacian system with one unknown per tree. The method ends when the KKT conditions hold to rounding error, which makes the result the exact minimizer; stats.polish records the outcome. If the method fails, the solver returns the multiplier iterate and warns.

See also: symcover_min, cover, cover_min!.

a, b, stats = cover_min!(ϕ, a, b, A; kwargs...)
a, b, stats = cover_min!(a, b, A; kwargs...)

Refine the asymmetric hard cover (a, b) in place. The start must cover A and be positive on supported rows and columns; unsupported scales are zeroed. The no-ϕ form uses AbsLog{2}(); supported penalties, keywords, stats, and failure behavior match cover_min.

Equivalent starts (c*a, b/c) give the same balanced result.

See also: initialize_cover, cover_min, symcover_min!.

cover_objective(ϕ, a, b, A)
cover_objective(ϕ, a, A)

Compute ∑ ϕ(|A[i,j]|/(a[i]*b[j])). The shorter form uses the symmetric cover a*a'.

Both forms sum over every entry of A, so in the symmetric form each off-diagonal pair contributes twice and each diagonal entry once.

Zero entries of A are handled according to ϕ:

• AbsLog{p}, PowerMean{p}: zero entries contribute 0 (φ(0) = 0 by convention).
• AbsLinear{p}: zero entries contribute 1 (φ(0) = |1-0|^p = 1).

eachindex(a) must match axes(A, 1) and eachindex(b) must match axes(A, 2).

The nonzero entries of A are enumerated by foreach_support.

See also:

• Penalty types (options for ϕ): AbsLog, AbsLinear, PowerMean.
• Solvers that minimize it: symcover_min, cover_min, soft_symcover, soft_cover.

a, b, stats = cover_transversal!(a, b, A; kwargs...)

Mutating counterpart of cover_transversal: writes the cover into a and b and returns them with the solver's statistics. The initial contents of a and b are ignored. eachindex(a) must match axes(A, 1) and eachindex(b) must match axes(A, 2).

a, b, stats = cover_transversal(A; kwargs...)

Return the minimal hard cover of the square matrix A among those that are tight on its maximum-product transversals.

A transversal of A is a set of entries A[i,σ(i)], one in each row and each column, for a permutation σ. Let

π* = max over permutations σ of ∏ᵢ |A[i,σ(i)]|.

Every hard cover has prod(a) * prod(b) >= π*, with equality exactly when it is tight (a[i]*b[j] == |A[i,j]|) on every entry of every transversal whose product is π*. Among hard covers, these maximize |det(A ./ (a .* b'))|, whereas a cover that leaves such a transversal loose can make the scaled matrix far more ill-conditioned than necessary. cover_transversal returns the hard cover with prod(a) * prod(b) == π* that minimizes the AbsLog{2} objective ∑ log(a[i]*b[j]/|A[i,j]|)^2 over the support of A.

The products a[i]*b[j] are unique on the support, they do not depend on which maximum-product transversal the solver finds, and they are scale-covariant: for D₁ * A * D₂ with positive diagonal D₁, D₂, they are multiplied by D₁[i,i] * D₂[j,j]. When the minimal cover cover_min is already tight on a maximum-product transversal, the two coincide. The factors use the balance convention of cover_min.

A must be square and structurally nonsingular, i.e., have at least one transversal of nonzero entries; otherwise an ArgumentError is thrown.

stats is an AugmentedLagrangianStats; as for cover_min, the solver warns when MatrixCovers.converged(stats) is false.

Extended help

A maximum-product transversal σ is found by shortest augmenting paths on the costs log(max_k |A[k,j]|) - log|A[i,j]|. Tightness on σ determines b from a, b[σ[i]] = |A[i,σ[i]]| / a[i], so the remaining problem is posed on the row scales alone, with one constraint per entry off the transversal. The augmented-Lagrangian solver of cover_min solves it, starting from the dual variables of the transversal, which are a hard cover tight on σ.

The keywords κ, maxouter, maxiter, tol, fillbudget, flopbudget, and linsolve are as for cover_min. linsolve=:woodbury solves the row-scale system as a sparse matrix (conjugate gradients, or sparse Cholesky when poorly conditioned). For n×n A, it requires at most n ÷ 4 zeros in every row and column, at most 4n zeros in total, and Float64 arithmetic (narrower types promote to Float64), and throws an ArgumentError otherwise. :auto chooses :lsqr when the stored support fills at most a quarter of the grid, :woodbury when its requirements hold, and :dense otherwise. Sparse matrices default to :lsqr. As for cover_min, an active-set method finishes a solve whose multiplier iteration stops early.

See also: cover_transversal!, cover_min.

foreach_support(f, A)

Call f(i, j, abs(A[i,j])) once for each nonzero entry of A, in unspecified order, and return nothing. The indices i, j are those of axes(A).

Extending

To visit only the stored entries of a new sparse matrix type, in O(nnz) time, define

MatrixCovers.foreach_support(f, A::MyMatrix)

The method must visit each nonzero entry exactly once, skip stored zeros, ignore the return values of f, and return nothing.

See also: foreach_support_sym.

foreach_support_sym(f, A)

Symmetric counterpart of foreach_support: call f(i, j, abs(A[i,j])) once for each nonzero pair {i, j}, with i <= j, including diagonal entries.

A must be square and abs.(A) symmetric. The traversal checks only that A is square; the package's symmetric solvers check symmetry before calling it.

Objective weighting

Because each off-diagonal pair is visited once, a sum over all entries of A must weight off-diagonal pairs by 2 and diagonal entries by 1. A coverage constraint a[i]*a[j] >= |A[i,j]| is the same for (i, j) and (j, i), so constraints need no weighting.

Extending

To support a new matrix type, define

MatrixCovers.foreach_support_sym(f, A::MyMatrix)

The method must visit each nonzero pair once with i <= j, skip zero pairs, ignore the return values of f, and return nothing. For storage that holds only the lower triangle, an entry stored at (i, j) with i > j must be passed as (j, i).

See also: foreach_support.

s = gramcover!(s, A)
s = gramcover!(s, A, w::AbstractVector)
s = gramcover!(s, A, W::Diagonal)
s = gramcover!(s, A, W::Cholesky)
s = gramcover!(s, A, W::AbstractMatrix; v=<symmetric cover of abs.(W)>)

Mutating counterpart of gramcover: writes the symmetric cover of the (weighted) Gram matrix into s and returns it. eachindex(s) must equal axes(A, 2).

See also: gramcover.

s = gramcover(A)
s = gramcover(A, w::AbstractVector)
s = gramcover(A, W::Diagonal)
s = gramcover(A, W::Cholesky)
s = gramcover(A, W::AbstractMatrix; v=<symmetric cover of abs.(W)>)

Return a vector s whose outer product s*s' covers the Gram matrix G = A'*W*A: s[j]*s[k] >= abs(G[j,k]) for all j, k. G is never formed. gramcover(A) covers A'*A, and passing w is equivalent to passing Diagonal(w). Columns of A with no stored nonzero get s[j] = 0.

One might expect to start from a cover a, b = cover(A), but none is needed: W supplies the row scale, which is held fixed while s is computed from the columns of A.

Every method is scale-covariant: replacing A by D1*A*D2 and W by inv(D1)*W*inv(D1) (for positive diagonal D1, D2) multiplies s by diag(D2).

• Diagonal weight (w, Diagonal(w), or none): s[j] = sqrt(Σ_i abs(w[i]) * abs(A[i,j])^2). The weights may have any sign. When w .>= 0 this is sqrt(G[j,j]), which is the minimal symmetric cover of G. The cost is one pass over the stored nonzeros of A, plus a second pass when the squares in some column underflow or overflow.
• Cholesky weight: for W = L*L', s[j] = norm(L'*A[:,j]) = sqrt(G[j,j]), the minimal symmetric cover of the positive semidefinite G. Passing the factorization is the caller's assertion that W is positive semidefinite. The cost is that of the dense product L'*A, O(m²n) for A of size m×n.
• General weight W: s[j] = Σ_i v[i] * abs(A[i,j]), where v .>= 0 satisfies abs(W[i,k]) <= v[i]*v[k] for all i, k. W may be indefinite or asymmetric. This bound is looser than the other two: for diagonal W with w .>= 0 it exceeds sqrt(G[j,j]) by up to a factor sqrt(n_j), where n_j is the number of nonzeros in column j. Prefer the Diagonal or Cholesky method whenever W has that structure. The default v is symcover(max.(abs.(W), abs.(transpose(W)))); a supplied v is checked for these conditions. The cost is that of computing v plus one pass over A.

W must be square with both axes equal to axes(A, 1), and s has axes (axes(A, 2),). The bounds hold in floating-point arithmetic: s is inflated by a relative O(n_j * eps) to absorb the rounding of the sums. For the Cholesky method the bound is relative to the computed factor. With a diagonal weight, s[j] has this accuracy whenever the largest sqrt(abs(w[i])) * abs(A[i,j]) in column j is at least floatmin, even if its square is not representable; an s[j] beyond floatmax is Inf.

See also: gramcover!, symcover, cover, iscover.

Examples

julia> using LinearAlgebra

julia> A = [4.0 1.0 0.0; 1.0 3.0 2.0; 0.0 2.0 5.0; 1.0 0.0 1.0];

julia> w = [1.0, -2.0, 0.5, 3.0];   # mixed signs

julia> s = gramcover(A, w);

julia> all(s * s' .>= abs.(A' * Diagonal(w) * A))
true

julia> W = [2.0 1.0 0.0 0.0; 1.0 2.0 1.0 0.0; 0.0 1.0 2.0 1.0; 0.0 0.0 1.0 2.0];

julia> s = gramcover(A, cholesky(W));

julia> all(s * s' .>= abs.(A' * W * A))
true

a, b = initialize_cover!(a, b, A; strategy=:hardcover, feasible=:inflate, kwargs...)

Mutating counterpart of initialize_cover: writes the starting cover into a and b and returns them, rather than allocating new vectors. eachindex(a) must match axes(A, 1) and eachindex(b) must match axes(A, 2).

See also: initialize_cover.

a, b = initialize_cover(A; strategy=:hardcover, feasible=:inflate, kwargs...)

Build an asymmetric starting point for cover_min or soft_cover. The feasible keyword matches initialize_symcover.

strategy names the point:

• :geomean — geometric means of the nonzero entries in each row and column.
• :covariant — the scale-covariant start of cover, without conjugate-gradient refinement. It agrees with :geomean on complete support.
• :hardcover — the result of cover, which is already a cover, so feasible is ignored. Additional keywords (maxiter, cgiter, start) are passed to cover; the other strategies accept none.

Supported rows and columns receive positive scales; unsupported ones receive zero. The factors use the balance convention of cover_min.

See also: initialize_cover!, initialize_symcover, cover, cover_min.

a = initialize_symcover!(a, A; strategy=:hardcover, feasible=:inflate, kwargs...)

Mutating counterpart of initialize_symcover: writes the starting cover into a and returns it, rather than allocating a new vector. eachindex(a) must match axes(A, 1) (and A must be square).

See also: initialize_symcover.

a = initialize_symcover(A; strategy=:hardcover, feasible=:inflate, kwargs...)

Build a symmetric starting point for symcover_min or soft_symcover. The point depends only on A, not on a penalty.

strategy selects the point:

• :geomean — the unconstrained AbsLog{2} start of symcover: geometric means of the diagonally normalized entries in each row.
• :leaveout — the geometric mean recomputed with one support entry omitted: the entry that the :geomean scales overestimate the most, i.e., the one with the smallest ratio abs(A[i, j]) / (a[i] * a[j]). It throws an ArgumentError if removing that entry empties a row.
• :diagfeasible — a cover grown from the diagonal by nearest-neighbor propagation.
• :hardcover — the result of symcover, which is already a cover, so feasible is ignored. Additional keywords (maxiter) are passed to symcover; the other strategies accept none.

feasible controls whether and how the point is made into a cover:

• :inflate multiplies every scale by the smallest common factor that covers A.
• :boost splits each violated entry's shortfall (in log magnitude) between the scales it touches, and raises each scale by the largest share it is asked for.
• :none returns the strategy's point without a coverage guarantee.

Supported rows receive positive scales; unsupported rows receive zero.

See also: initialize_symcover!, initialize_cover, symcover, symcover_min.

iscover(a, b, A; rtol=0, atol=0)
iscover(a, A; rtol=0, atol=0)

Test whether a and b cover A, that is, whether a[i]*b[j] >= abs(A[i,j]) for every entry. The two-argument form tests the symmetric cover a*a', and requires A to be square.

rtol and atol allow for small violations, testing

a[i]*b[j] >= abs(A[i,j])*(1 - rtol) - atol

Both tolerances default to zero. A nonzero atol breaks scale invariance.

The results of cover, symcover, and the native AbsLog{2} minimal covers pass this test at zero tolerance. Initializers and extension solvers may require a nonzero rtol. Soft covers are not hard covers and in general fail it.

a and b must be nonnegative; a negative scale raises an ArgumentError. Zero is allowed for unsupported rows and columns. eachindex(a) must match axes(A, 1) and eachindex(b) must match axes(A, 2).

See also: cover_objective, cover, symcover.

Examples

julia> A = [1.0 2.0; 3.0 4.0];

julia> a, b = cover(A);

julia> iscover(a, b, A)
true

julia> iscover([1.0, 1.0], [1.0, 1.0], A)   # a*b' = ones, which does not reach A[2,2]
false

iterations(stats::AbstractCoverStats) -> Int

Return the number of iterations of the solver's main loop; each stats type states which loop that is.

level_profile((a, b), (aref, bref), A)
level_profile(a, aref, A)

Compare the product a*b' with the reference product aref*bref' on the nonzero entries of A, and return the named tuple (; level, profile). With

ℓ[i,j] = log((a[i]*b[j]) / (aref[i]*bref[j]))

level is the mean of ℓ over the nonzero entries of A, and profile is its standard deviation, normalized by the number of nonzero entries. The three-argument form compares the symmetric products a*a' and aref*aref'.

A positive level means that a*b' is larger than the reference on average. Multiplying either product by a constant c shifts the level by ±log(c) and leaves the profile unchanged. Bounds that rest on the cover inequality depend on both, whereas cover-normalized quantities such as condition numbers depend only on the profile. Both are unchanged when A and all four scales are rescaled together, as under a change of units.

The products need not be covers, but each must be positive and finite on every nonzero entry of A, and A must have at least one nonzero entry. The scales must be nonnegative, and their indices must match those of A, as for iscover.

See also: max_coverage_ratio, iscover.

Examples

julia> A = [1.0 2.0; 3.0 4.0];

julia> a, b = cover(A);

julia> level_profile((2a, b), (a, b), A)   # a pure level
(level = 0.6931471805599453, profile = 0.0)

julia> level_profile(([2.0, 0.5], [1.0, 1.0]), ([1.0, 1.0], [1.0, 1.0]), [1.0 0.0; 0.0 1.0])   # a pure profile
(level = 0.0, profile = 0.6931471805599453)

max_coverage_ratio(a, b, A)
max_coverage_ratio(a, A)

Return the largest coverage ratio abs(A[i,j])/(a[i]*b[j]) over the nonzero entries of A. The two-argument form uses the symmetric cover a*a', and requires A to be square.

a and b cover A exactly when the result r is at most one, so this is the quantitative form of iscover. Scaling the product a*b' by r makes it a cover that is tight at the worst entry, whether r is above or below one; in the symmetric form that scaling is sqrt(r)*a.

Zero entries of A impose no constraint and are skipped. A nonzero entry over a zero product gives Inf, and a matrix with no nonzero entries gives zero. The ratio is unchanged when A, a, and b are rescaled together, as under a change of units.

a and b must be nonnegative, and their indices must match those of A, as for iscover.

See also: iscover, cover_objective.

Examples

julia> A = [1.0 2.0; 3.0 4.0];

julia> max_coverage_ratio([1.0, 1.0], [1.0, 2.0], A)   # A[2,1] = 3 against a[2]*b[1] = 1
3.0

julia> max_coverage_ratio([2.0, 2.0], A)
1.0

ncomponents(sc::SupportComponents) -> Int

Number of connected components, so component ids run 1:ncomponents(sc).

residual(stats::AbstractCoverStats)

Return the residual of the returned cover that the stopping criterion compares with tolerance(stats).

rowcomponent(sc::SupportComponents, i) -> Int

Component id of row i, or 0 if that row has empty support. i is the matrix's own row index.

MatrixCovers.scalar_type(T)

Return the unitless floating-point type underlying T. Unit-carrying element types should specialize this method.

a, b, stats = soft_cover!(ϕ, a, b, A; kwargs...)
a, b, stats = soft_cover!(a, b, A; kwargs...)

Refine the starting point (a, b) into a soft cover of A in place. Scales must be finite and positive on supported rows and columns; unsupported scales are zeroed. The start need not cover A. Build one with initialize_cover and feasible=:none. The no-ϕ form uses PowerMean{2}().

For the convex penalties (PowerMean, AbsLog) the result is the minimizer returned by soft_cover, with the same keywords, and does not depend on the start. For AbsLinear it is the local minimum reached from the start, and there are no keywords. The result uses the balance convention of cover_min, so equivalent rescalings (c*a, b/c) give the same result.

See also: soft_cover, initialize_cover, soft_symcover!.

a, b, stats = soft_cover(ϕ, A; kwargs...)
a, b, stats = soft_cover(A; kwargs...)

Return the scales a, b that minimize ∑ ϕ(|A[i,j]|/(a[i]*b[j])). This is the asymmetric counterpart of soft_symcover: it penalizes violations of a[i]*b[j] >= |A[i,j]| rather than forbidding them. The default penalty is PowerMean{2}(). stats and the failure behavior are as for soft_symcover.

Supported penalties and their keywords are:

• PowerMean{p}() (default): the convex minimum. The solver alternates power-mean updates of the row and column scales, with adaptive overrelaxation, and switches to Newton steps with backtracking on the objective when these sweeps are slow. The sweeps amount to Sinkhorn scaling of abs.(A).^p to row sums n_i and column sums m_j, the nonzero counts. The iteration stops when, in every nonzero row and column, the mean of r^p over the ratios r of the nonzero entries is within tol of one. The computation uses at least Float64 precision; tol defaults to 4096*eps of the working type, and a warning is issued if the iteration ends without reaching it. maxiter (default 10_000) bounds the number of sweeps, each an update of all rows and then all columns. The sweeps hand over to Newton when their rate of progress predicts more than 100 further sweeps, and after at most min(maxiter, 200) sweeps. newton, maxnewton, and linsolve have the meanings given in the extended help of soft_symcover.
• AbsLog{2}(): the convex minimum, computed by one linear solve; linsolve has the same meaning as in cover_min.
• AbsLog{1}(): the convex minimum, computed as a linear program; requires JuMP and HiGHS. Among the minimizers, the one with the least AbsLog{2} objective is returned.
• AbsLinear{1}(), AbsLinear{2}(): requires JuMP and Ipopt. The objective is not convex, so several starts are each refined to a local minimum, and the result with the least objective is returned. starts (default 4) is the number of starts: the boosted covariant start of cover and the hard cover tightened from it, then log-normal perturbations of the first with spread σ (default 2.0; sigma is an alias), drawn from rng (default MersenneTwister(0)). A perturbed start whose refinement throws a SolverFailure is left out of the selection; a failing deterministic start throws.

When abs.(A) is exactly symmetric, the PowerMean minimizer has b == a under the balance convention. In that case soft_cover computes a by the algorithm of soft_symcover (with the same keywords; maxiter then counts its updates) and returns (a, copy(a), stats). The in-place soft_cover! always uses the alternating iteration.

Unsupported rows and columns receive zero scale. The factors use the balance convention of cover_min.

See also: cover, soft_symcover, soft_cover!, cover_objective, PowerMean.

Examples

julia> A = [1 2 3; 6 5 4];

julia> a, b, stats = soft_cover(A);

julia> a * b'
2×3 Matrix{Float64}:
 1.73729  1.93617  2.37048
 4.64478  5.17648  6.33766

a, stats = soft_symcover!(ϕ, a, A; kwargs...)
a, stats = soft_symcover!(a, A; kwargs...)

Refine one symmetric soft-cover start in place. The no-ϕ form uses PowerMean{2}(). Build a start with initialize_symcover and feasible=:none.

a must be finite and positive on supported rows. It need not cover A, and unsupported scales are set to zero.

For the convex penalties (PowerMean, AbsLog) the result is the minimizer returned by soft_symcover, with the same keywords, and does not depend on the start. For AbsLinear it is the local minimum reached from the start, and there are no keywords.

See also: soft_symcover, initialize_symcover, soft_cover!.

a, stats = soft_symcover(ϕ, A; kwargs...)
a, stats = soft_symcover(A; kwargs...)

Return the scales a that minimize ∑ ϕ(|A[i,j]|/(a[i]*a[j])) for symmetric A. Unlike a hard cover, this soft cover penalizes violations of a[i]*a[j] >= |A[i,j]| rather than forbidding them. The default penalty is PowerMean{2}().

Supported penalties and their keywords are:

• PowerMean{p}() (default): the convex minimum. The solver applies damped simultaneous power-mean updates a[k] ← sqrt(a[k] * M_p(|A[k,j]|/a[j])), where M_p is the p-power mean over the nonzeros of row k, and switches to Newton steps with backtracking on the objective when the updates are slow. Every update decreases the objective. The iteration stops when, in every nonzero row, the mean of r^p over the ratios r = |A[k,j]|/(a[k]*a[j]) of the nonzero entries is within tol of one. The computation uses at least Float64 precision; tol defaults to 4096*eps of the working type, and a warning is issued if the iteration ends without reaching it. The keywords maxiter, newton, maxnewton, and linsolve are described in the extended help.
• AbsLog{2}(): the convex minimum, computed by one linear solve; linsolve has the same meaning as in symcover_min.
• AbsLog{1}(): the convex minimum, computed as a linear program; requires JuMP and HiGHS. The minimizer need not be unique, and among the minimizers the one with the least AbsLog{2} objective is returned.
• AbsLinear{1}(), AbsLinear{2}(): requires JuMP and Ipopt. The objective is not convex, so several starts are each refined to a local minimum, and the result with the least objective is returned. starts (default 5) is the number of starts: a few deterministic starts, then log-normal perturbations of the first deterministic start with spread σ (default 2.0; sigma is an alias), drawn from rng (default MersenneTwister(0)). A perturbed start whose refinement throws a SolverFailure is left out of the selection; a failing deterministic start throws.

stats is an AbstractCoverStats. The JuMP-backed solves (AbsLog{1}, AbsLinear) throw a SolverFailure when the external solver fails.

Unsupported rows receive zero scale. When a connected component of the support is bipartite with no diagonal entry, the products on the support do not determine the scales; PowerMean and AbsLog then apply the balance convention of cover_min between the two sides of the bipartition.

See also: symcover, cover_objective, soft_symcover!, PowerMean.

Examples

julia> A = [4 -1; -1 1];

julia> a, stats = soft_symcover(A);

julia> MatrixCovers.converged(stats)
true

julia> R = abs.(A) ./ (a .* a');

julia> round.(sum(R .^ 2; dims=2); digits=8)   # 2 nonzeros per row: each mean square ratio is one
2×1 Matrix{Float64}:
 2.0
 2.0

julia> round.(first(soft_symcover([0 1; 1 0])); digits=4)
2-element Vector{Float64}:
 1.0
 1.0

Extended help

PowerMean keywords

• maxiter (default 10_000): the largest number of power-mean updates.
• newton (default true): whether Newton steps finish the solve. The updates hand over to Newton when the rate of decrease of the imbalance over recent updates predicts more than 200 further updates, and after at most min(maxiter, 400) updates. With newton=false the updates continue until convergence or maxiter.
• maxnewton (default 100): the largest number of Newton steps. Each step solves the Newton equations and backtracks on the objective.
• linsolve (default :auto): the solver for the Newton equations. :dense factorizes the dense Hessian; :cholesky factorizes the sparse Hessian with CHOLMOD (Float64 only); :cg uses conjugate gradients with a diagonal preconditioner, at one pass over the support per iteration, for problems too large to factor. :auto chooses :dense for support that fills at least a quarter of the grid and otherwise :cholesky, each when its predicted flop count is at most 8e3 per stored entry (and, for :cholesky, its predicted storage at most 2^30 bytes), and :cg when neither fits. Poorly conditioned problems can need many conjugate-gradient iterations.

support_components(A) -> sc::SupportComponents

Return the connected components of A's bipartite support graph. The matrix is read through foreach_support.

See also: SupportComponents.

a = symcover!(ϕ, a, A; maxiter=3)
a = symcover!(a, A; maxiter=3)

Mutating counterpart of symcover: writes the symmetric hard cover into a and returns it, rather than allocating a new vector. eachindex(a) must match axes(A, 1) (and A must be square). ϕ has the same meaning as in symcover, and is likewise ignored by the current heuristic covers.

See also: symcover.

a = symcover(ϕ, A; maxiter=3)
a = symcover(A; maxiter=3)

Return a vector a whose outer product a*a' is a hard cover of the square, symmetric matrix A: a[i] * a[j] >= abs(A[i, j]) for all i, j. This is a fast heuristic; for the cover that minimizes a penalty ϕ, use symcover_min.

The heuristic starts from the per-row geometric means of the diagonally normalized entries abs(A[i, j]) / sqrt(abs(A[i, i] * A[j, j])), rescaled by sqrt(abs(A[i, i])). It then raises the scales until every entry is covered and applies maxiter tightening iterations.

ϕ is accepted for API compatibility but is ignored.

See also: symcover!, symcover_min, soft_symcover, cover.

Examples

julia> A = [4 1; 1 4];

julia> a = symcover(A)
2-element Vector{Float64}:
 2.0
 2.0

julia> a * a'   # covers |A|: a[i]*a[j] >= abs(A[i, j])
2×2 Matrix{Float64}:
 4.0  4.0
 4.0  4.0

Extended help

The result is scale-covariant up to floating-point round-off (covering D * A * D for positive diagonal D gives diag(D) .* a) whenever every connected component of the support has a nonzero diagonal entry; rows with a zero diagonal take their reference from neighbors that have one. Covariance is not exact even for power-of-two scalings, because the method works with logarithms of the entries.

The result is also permutation-equivariant up to floating-point round-off: covering A[p, p] for a permutation p gives a[p].

a, stats = symcover_min(ϕ, A; kwargs...)
a, stats = symcover_min(A; kwargs...)

Return the ϕ-minimal symmetric hard cover of A: the vector a minimizing ∑ ϕ(|A[i,j]|/(a[i]*a[j])) subject to a[i]*a[j] >= |A[i,j]|. The default penalty is AbsLog{2}().

Supported ϕ values:

• AbsLog{2}(): native.
• AbsLog{1}(): requires JuMP and HiGHS.
• AbsLinear{1}() and AbsLinear{2}(): require JuMP and Ipopt and return the best local minimum found from strategies, a collection of initialize_symcover strategy names (default (:hardcover, :geomean, :leaveout); a strategy that cannot build a start is skipped).

The AbsLog objectives are convex in the logarithms of the scales, so their local minima are global. When the AbsLog{1} minimum is not unique, the method returns the minimizer with the smallest AbsLog{2} objective.

stats is an AbstractCoverStats; MatrixCovers.converged(stats) reports whether a meets the solver's stopping criterion. The native solver warns when it does not, and returns a cover of A that may not minimize the objective. JuMP-backed solves throw a SolverFailure instead.

Extended help

The native solver accepts κ (initial augmented-Lagrangian penalty, default 1e2), maxouter (multiplier updates, default 32; 0 returns the unconstrained fit, shifted to cover A, with converged(stats) == false and no warning), maxiter (Newton steps per update), tol (the tolerance on the KKT residual, default 1000*eps of the working type), fillbudget and flopbudget (see below), and linsolve:

• :dense factorizes dense normal equations at O(n³) per Newton step.
• :woodbury handles nearly dense Float64 support as a sparse correction. It requires at most n ÷ 4 missing entries per row and 4n in total.
• :lsqr is matrix-free with O(nnz) work per iteration and is the sparse-matrix default.
• :auto chooses :woodbury when supported, :lsqr when the stored support fills at most a quarter of the grid, and :dense otherwise.

The solver increases κ when the KKT residual contracts slowly, and when a positive multiplier remains on an entry with slack.

The budgets choose the preconditioner of :lsqr. For Float64, it is a sparse Cholesky factorization when the factor's predicted storage is at most fillbudget bytes (default 2^30) and its predicted flop count is at most flopbudget (default 8e3) times the number of stored entries in the support of A, counting each symmetric off-diagonal pair twice, as nnz does for a fully stored matrix. Otherwise, and for every other element type, it is a factorization of a maximum-weight spanning forest of the couplings among entries whose constraints are active, together with the diagonal contributions of all entries; this factorization has no fill and needs no budget. fillbudget=Inf and flopbudget=Inf together always select the Cholesky preconditioner; fillbudget=0 or flopbudget=0 always select the spanning-forest one. The same budgets apply to the active-set method that finishes an unconverged solve (see cover_min): graph-Laplacian systems whose sparse Cholesky factor would exceed them are solved by conjugate gradients instead.

If the solver warns that the result may not minimize the objective, follow the advice in the warning: increase maxouter when the update limit was reached, and κ when the residual stopped contracting before the limit. A caller who accepts a larger residual can pass a larger tol, or compare MatrixCovers.residual(stats) with its own threshold.

When converged(stats) is true, the result typically matches the exact minimizer to a relative accuracy of about 1e-8, and occasionally only 1e-7: the stopping test measures feasibility and complementary slackness but not stationarity. Scale covariance holds to the same accuracy. For a diagonal rescaling A ./ (d * d'), the result approximates symcover_min(A) ./ d but does not equal it bitwise, even when every d[i] is a power of two.

The native solver computes in Float64 for narrower input types, then converts the result to the required element type.

See also: cover_min, symcover, symcover_min!.

a, stats = symcover_min!(ϕ, a, A; kwargs...)
a, stats = symcover_min!(a, A; kwargs...)

Refine the symmetric hard cover a in place. The no-ϕ form uses AbsLog{2}(); supported penalties, keywords, stats, and failure behavior match symcover_min.

a must cover A and be positive on supported rows. Unsupported scales are ignored on input and set to zero. Use initialize_symcover or symcover to construct a start.

The AbsLog result is independent of the start. For AbsLog{1}, ties are broken by the AbsLog{2} objective. Local minima under AbsLinear can depend on the start.

See also: initialize_symcover, symcover_min, cover_min!.

tolerance(stats::AbstractCoverStats)

Return the tolerance that residual(stats) is compared with.