GitHub

MatrixCovers

Given a matrix A, a hard cover is C = a * b', where a and b are nonnegative vectors satisfying

\[C_{ij} \;\geq\; |A_{ij}| \quad \text{for all } i, j.\]

For symmetric A, a single vector suffices (b = a). A minimal cover minimizes a chosen penalty, while a soft cover penalizes violations instead of enforcing every inequality.

Why covers?

Covers provide the "natural scales" of a matrix. If you rescale rows by a positive diagonal factor D_r and columns by D_c, the optimal cover transforms as a → D_r * a, b → D_c * b, so the product a * b' transforms identically to A. Moreover, Ahat = A ./ (a * b') is scale-invariant. Scalar metrics like norm(A) or maximum(abs, A) implicitly encode an arbitrary choice of units, but applying them to Ahat rather than A fixes this deficiency.

While most users will employ matrices that store pure numbers, we'll start with an example of a 3×3 matrix whose rows and columns correspond to physical variables with different units — position in meters, velocity in m/s, force in Newtons. Loading Unitful lets the matrix express those units directly:

julia> using MatrixCovers, Unitful

julia> L, V, F = u"m", u"m/s", u"N";  # length, velocity, force

julia> A = [1e6/L^2   1e3/(L*V)  1.0/(L*F)
            1e3/(L*V) 1.0/V^2    1e-3/(V*F)
            1.0/(L*F) 1e-3/(V*F) 1e-6/F^2];

julia> a = symcover(A);

julia> round.(typeof.(a), a; digits=6)
3-element Vector{Quantity{Float64}}:
 1000.0 m^-1
    1.0 s m^-1
    0.001 N^-1

Here A[i,j] has units 1/(u[i]*u[j]), as in a Hessian. Its cover identifies scales of 1 mm, 1 m/s, and 1 kN. Normalization removes both units and the magnitudes affected by choice of units:

julia> round.(A ./ (a .* a'); digits=6)
3×3 Matrix{Float64}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

An entry is 1 only where the cover bound is tight, and this is not guaranteed for all matrices. For example, given diagonal A, the normalized matrix is also diagonal.

A cover exists only when the units of A factor as unit(A[i,j]) == unit(a[i])*unit(b[j]). But this is not an onerous requirement, as it is the same one that lets expressions like A*x be well-defined. If a matrix can be used in matrix-vector multiplication, it has a cover.

Penalty functions

A cover is valid as long as every constraint is satisfied, but tighter covers better capture the scaling of A. Cover quality is measured through the ratios r[i, j] = |A[i,j]| / (a[i] * b[j]): a hard cover has every 0 ≤ r[i, j] ≤ 1, and r[i, j] == 1 means the constraint is exactly tight.

A penalty function ϕ combines those ratios into a scalar objective

\[\sum_{i,j} \phi\!\left(\frac{|A_{ij}|}{a_i\, b_j}\right),\]

Three penalty families are provided:

  • AbsLog{p}: ϕ(r) = |log r|^p, with ϕ(0) = 0. It is convex in log space and is the default for hard covers.
  • PowerMean{p}: ϕ(r) = (r^p - 1 - p log r)/p, with ϕ(0) = 0. It is strictly convex in log space and is the default for soft covers. Its minimizer makes the mean of r^p over the nonzero entries of every row and column equal to one.
  • AbsLinear{p}: ϕ(r) = |1-r|^p. It is nonconvex and finite at zero, so an entry that shrinks to zero drops out of the soft cover continuously.

cover_objective evaluates any penalty for a given cover:

julia> using MatrixCovers

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

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

julia> cover_objective(AbsLog{1}(), a, A)   # sum of log-excesses (L1)
2.772588722239781

julia> cover_objective(AbsLog{2}(), a, A)   # sum of squared log-excesses (L2)
3.843624111345611

Pass a penalty as the first solver argument to override the default.

Choosing a cover algorithm

FunctionSymmetricConstraintDefault (or alternative) objectiveRequires
symcoveryeshard (r ≤ 1)heuristic—
covernohard (r ≤ 1)heuristic—
symcover_minyeshard (r ≤ 1)AbsLog{2} (or AbsLog{1}, AbsLinear)native for AbsLog{2}; else JuMP
cover_minnohard (r ≤ 1)AbsLog{2} (or AbsLog{1}, AbsLinear)native for AbsLog{2}; else JuMP
cover_transversalnohard, tight on a maximum-product transversalAbsLog{2}—
soft_symcoveryessoft (penalized)PowerMean{2} (or PowerMean{p}, AbsLog, AbsLinear)native for PowerMean, AbsLog{2}; else JuMP
soft_covernosoft (penalized)PowerMean{2} (or PowerMean{p}, AbsLog, AbsLinear)native for PowerMean, AbsLog{2}; else JuMP

For hard covers, symcover and cover are fast heuristics; their _min counterparts minimize the selected objective. cover_transversal minimizes the AbsLog{2} objective of a square, structurally nonsingular matrix among the hard covers that are tight on a transversal of maximum product, the covers that maximize |det| of the scaled matrix. For such matrices cover returns a heuristic cover with the same tightness by default (start=:transversal); a cover that leaves a maximum-product transversal loose can make the scaled matrix far more ill-conditioned than necessary. A soft cover has no constraint to satisfy, so there is no separate heuristic tier: soft_symcover and soft_cover always minimize the objective.

  • PowerMean (native) uses power-mean scaling iterations (damped simultaneous updates for the symmetric problem, overrelaxed alternating updates for the asymmetric one), followed by damped Newton steps when those iterations converge slowly.
  • AbsLog{2} (native) is one linear least-squares solve.
  • 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.

The heuristics cost $O(mn)$, except that the transversal start of cover also solves an assignment problem by shortest augmenting paths; native iterative solvers cost roughly $O(mn)$ per iteration.

Covariance of the heuristics

symcover is exactly covariant whenever every connected component of the support contains a nonzero diagonal entry:

julia> using MatrixCovers, LinearAlgebra

julia> A = [1.0 1 0; 1 1 1; 0 1 1];

julia> d = [1.0, 6.0, 0.5]; D = Diagonal(d);

julia> a1 = symcover(A); a2 = symcover(D * A * D);

julia> P1 = (d .* a1) * (d .* a1)'; P2 = a2 * a2';

julia> round.(extrema(P2 ./ P1); digits=3)
(1.0, 1.0)

cover is covariant on any support: positive diagonal rescaling of A rescales the cover products identically wherever A is nonzero. This holds for both the default transversal start (square, structurally nonsingular A) and the :covariant start, and for any number of conjugate-gradient refinement iterations (cgiter):

julia> e = [2.0, 0.3, 5.0]; E = Diagonal(e);

julia> function covariance_spread(cgiter; start=:auto)
           r1, c1 = cover(A; cgiter, start)
           r2, c2 = cover(D * A * E; cgiter, start)
           Q1 = (d .* r1) * (e .* c1)'
           Q2 = r2 * c2'
           return round.(extrema(Q2 ./ Q1); digits=3)
       end;

julia> covariance_spread(0)
(1.0, 1.0)

julia> covariance_spread(4)
(1.0, 1.0)

The start=:geomean alternative begins from the row and column geometric means. Unless the support is rank one, that point depends on the units A is written in, so the cover becomes scale-covariant only as the refinement converges:

julia> covariance_spread(0; start=:geomean)
(0.846, 1.236)

julia> covariance_spread(4; start=:geomean)
(1.0, 1.0)

For symmetric support components without a nonzero diagonal entry, use symcover_min when covariance is required.

Objective-minimal covers

symcover_min and cover_min minimize the chosen penalty subject to the hard constraint. The built-in AbsLog{2} solver uses penalty continuation with a damped semismooth Newton iteration:

julia> using MatrixCovers

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

julia> a, b = cover(A);

julia> aq, bq = cover_min(AbsLog{2}(), A);

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

julia> aq * bq'
2×3 Matrix{Float64}:
 2.21042  2.0      3.0
 6.0      5.42884  8.14325

julia> round(cover_objective(AbsLog{2}(), a, b, A); digits=6)
1.146449

julia> round(cover_objective(AbsLog{2}(), aq, bq, A); digits=6)
1.141281

AbsLog{1} and AbsLinear use JuMP with HiGHS and Ipopt, respectively:

julia> using MatrixCovers, JuMP, HiGHS

julia> S = [4 1 0; 1 1 5; 0 5 2];

julia> a, stats = symcover_min(AbsLog{1}(), S);

julia> round.(a; digits=6)
3-element Vector{Float64}:
 2.0
 1.0
 5.0

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

julia> a, b = cover_min(AbsLog{1}(), A);

julia> round.(a * b'; digits=6)
2×3 Matrix{Float64}:
 2.4  2.0  3.0
 6.0  5.0  7.5

soft_symcover and soft_cover solve PowerMean and AbsLog{2} natively, and use JuMP with HiGHS for AbsLog{1} and with Ipopt for AbsLinear.

Uniqueness

For asymmetric covers, a → γ*a, b → b/γ leaves a*b' unchanged. The package chooses a unique representative using ∑ n_i log a[i] = ∑ m_j log b[j], where n_i, m_j are the nonzero counts of row i and column j. The convention is applied separately to each connected component of the bipartite support graph. It affects the factors but not their products.

AbsLog{2} has a unique minimizer except when the support pattern leaves a scaling freedom, as in [0 1; 1 0], where every a with a[1]*a[2] = 1 is optimal. PowerMean has unique products a[i]*b[j] on the support, and a symmetric matrix has a symmetric minimizer; where a symmetric support leaves a scaling freedom, the two sides of the bipartition are balanced as above. If AbsLog{1}() has multiple minima, the implementation chooses the one with the smallest AbsLog{2} objective. AbsLinear may have several local minima.

Starting points: initialize and refine

For objectives with multiple minima, the result can depend on its starting point. The interface separates initialization, refinement, and multistart selection:

Plain forms choose their starts; ! forms refine the supplied start, except symcover!, cover!, and cover_transversal!, which write their result in place.

julia> using MatrixCovers, JuMP, Ipopt

julia> S = [4 1 0; 1 1 5; 0 5 2];

julia> a, stats = symcover_min(AbsLinear{2}(), S);

julia> round.(a; digits=6)
3-element Vector{Float64}:
 2.0
 1.0
 5.0

julia> a, stats = symcover_min(AbsLinear{2}(), S; strategies=(:geomean,));

julia> round.(a; digits=6)
3-element Vector{Float64}:
 2.0
 1.0
 5.0

julia> a0 = initialize_symcover(S; strategy=:geomean);

julia> symcover_min!(AbsLinear{2}(), a0, S);

julia> round.(a0; digits=6)
3-element Vector{Float64}:
 2.0
 1.0
 5.0

For the convex AbsLog and PowerMean penalties, the start does not change the result.

Convergence and solver statistics

The optimizing solvers return statistics as their last output, reporting whether the solver met its stopping criterion (the heuristics, which run a fixed number of steps, return none):

julia> using MatrixCovers

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

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

julia> MatrixCovers.converged(stats)
true

julia> MatrixCovers.residual(stats) <= MatrixCovers.tolerance(stats)
true

a, b = cover_min(A) discards the statistics, but a = symcover_min(A) binds the pair; write a, stats = symcover_min(A).

Each solver family returns its own subtype of AbstractCoverStats. The accessors converged, iterations, residual, and tolerance are public but not exported, and each subtype implements those that apply to its solver.

A native solver that stops before meeting its criterion warns and returns its result; a hard cover returned this way still covers A. Test MatrixCovers.converged(stats) to reject such results. The JuMP-backed solvers throw a SolverFailure instead.

Consuming one factor alone: gauges and Gram covers

The balanced factors of an asymmetric cover are deterministic but not individually covariant under one-sided scaling. This matters when consuming one factor, for example when covering J'*J from a cover of J.

gramcover(J[, W]) constructs a symmetric cover of J'*W*J directly from J and W, without a cover of J, and covaries with right-scaling of J. cover(J, a) gives the tightest cover of J for a fixed row scale a.

julia> using MatrixCovers, LinearAlgebra

julia> J = [1.0 2; 3 4; 5 6];

julia> D = Diagonal([100.0, 1.0]);

julia> a1, b1 = cover(J); a2, b2 = cover(J * D);

julia> r = b2 ./ (D.diag .* b1); all(x -> x ≈ first(r), r)
true

julia> first(r) ≈ 1
false

julia> s1 = gramcover(J); s2 = gramcover(J * D);

julia> s2 ≈ D.diag .* s1
true

Worked example: roundoff in A \ b

For x = A \ b, diagonal rescaling sends x → x ./ d and a symmetric cover a → d .* a. Thus sum(abs.(x .* a)) is invariant. The quantity can be estimated without solving for x:

julia> using MatrixCovers, LinearAlgebra

julia> A = [1e6 1e3; 1e3 4.0];

julia> b = [1.5e3, 6.0];

julia> a = symcover(A);

julia> round.(a; digits=6)
2-element Vector{Float64}:
 1000.0
    2.0

julia> mag = sum(abs(bi / ai) for (bi, ai) in zip(b, a));

julia> round(mag; digits=6)
4.5

Here mag is within a factor of 1.5 of the scaled solution norm:

julia> x = A \ b;

julia> round(sum(abs.(x .* a)); digits=6)
3.0

The estimate is unchanged by diagonal rescaling:

julia> d = [0.05, 3.0];

julia> Ad, bd = d .* A .* d', d .* b;

julia> ad = symcover(Ad);

julia> round(sum(abs(bi / ai) for (bi, ai) in zip(bd, ad)); digits=6)
4.5

For well-conditioned A, eps(mag) estimates the roundoff floor:

julia> xbig = big.(A) \ big.(b);

julia> abs(sum(abs.(x .* a)) - sum(abs.(Float64.(xbig) .* a))) <= 2 * eps(mag)
true

For ill-conditioned A, the error can be much larger:

julia> Aill = [1.0 -0.9999; -0.9999 1.0];

julia> bill = [0.75, 7.0];

julia> aill = symcover(Aill);

julia> magill = sum(abs(bi / ai) for (bi, ai) in zip(bill, aill));

julia> xill = Aill \ bill;

julia> xbigill = big.(Aill) \ big.(bill);

julia> err = abs(sum(abs.(xill .* aill)) - sum(abs.(Float64.(xbigill) .* aill)));

julia> err > 1e6 * eps(magill)
true

The condition number of the normalized matrix provides a corresponding bound:

julia> κ = cond(Aill ./ (aill .* aill'));

julia> err <= 1e3 * eps(κ * magill)
true

Index of available tools

Reference documentation

MatrixCovers.AbsLinear — Type
AbsLinear{p}

Penalty type for φ(r) = |1 - r|^p. Unlike AbsLog, this penalty is continuous at r = 0 (φ(0) = 1), so zero entries in A naturally contribute a constant penalty.

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

source
MatrixCovers.AbsLog — Type
AbsLog{p}

Penalty type for

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

The r=0 convention keeps zero entries finite. The objective is convex in log space; AbsLog{1} may have multiple minima.

See also: AbsLinear.

source
MatrixCovers.AbstractCoverPenalty — Type
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:

Passing a custom subtype to a solver raises a MethodError.

source
MatrixCovers.AbstractCoverStats — Type
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.

source
MatrixCovers.ActiveSetStats — Type
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.

source
MatrixCovers.AugmentedLagrangianStats — Type
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.
source
MatrixCovers.ExternalSolverStats — Type
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.

source
MatrixCovers.LeastSquaresStats — Type
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.
source
MatrixCovers.MultistartStats — Type
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.
source
MatrixCovers.PowerMean — Type
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, the r=0 convention keeps zero entries finite, so structural zeros of A contribute nothing.

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).

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.

See also: soft_cover, soft_symcover.

source
MatrixCovers.PowerMeanStats — Type
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).
source
MatrixCovers.SolverFailure — Type
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.

source
MatrixCovers.colcomponent — Method
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.

source
MatrixCovers.converged — Function
converged(stats::AbstractCoverStats) -> Bool

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

source
MatrixCovers.cover! — Method
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.

source
MatrixCovers.cover — Method
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
source
MatrixCovers.cover — Method
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. The heuristic uses up to cgiter conjugate-gradient iterations to refine its starting point, then ensures coverage and applies maxiter improvement passes. cgiter=0 disables refinement.

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, 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 currently ignored. For a cover that provably minimizes a given ϕ, use cover_min.

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

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, and no worse than about 1e-9 on ill-conditioned supports. cgiter caps the iterations, which stop earlier once the residual is at roundoff.

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, and the :covariant start has no mechanism to prevent this. The :transversal start finds a maximum-product transversal by shortest augmenting paths, fixes every row and column scale up to one unknown per group of rows and columns joined by entries on maximum-product transversals, and then fits the remaining entries: 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 :covariant start is exact on trees and complete bipartite components. Each refinement iteration makes one pass over the support; complete support needs no refinement.

The :covariant and :geomean starts coincide on complete support, where the geometric mean already solves the least-squares problem.

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].

source
MatrixCovers.cover_min — Function
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}(). Within each support component, the factors use the balance convention ∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j] (nzaᵢ, nzbⱼ = nonzero counts of row i, 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.

AbsLog is convex in the log scales. If AbsLog{1} has multiple minima, the method selects the one 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, and each step solves the objective with a spanning forest of them held tight, a sparse graph Laplacian 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 it. When it fails, the solver returns the multiplier iterate and warns.

See also: symcover_min, cover, cover_min!.

source
MatrixCovers.cover_min! — Function
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!.

source
MatrixCovers.cover_objective — Method
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 use full-grid weighting: symmetric off-diagonal pairs contribute twice and diagonal entries 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).

A is read through foreach_support.

See also:

source
MatrixCovers.cover_transversal! — Method
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).

source
MatrixCovers.cover_transversal — Method
a, b, stats = cover_transversal(A; kwargs...)

Return the transversal-tight minimal cover of the square matrix A: the hard cover minimizing ∑ log(a[i]*b[j]/|A[i,j]|)^2 over the support of A, among the hard covers with prod(a) * prod(b) == π*, where

π* = 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'))|.

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.

source
MatrixCovers.foreach_support — Method
foreach_support(f, A)

Call f(i, j, abs(A[i,j])) once per nonzero entry and return nothing. Traversal order is unspecified; indices follow axes(A).

Specialize this hook to support custom sparse storage in O(nnz) time.

Extending

To support a new matrix type, define

MatrixCovers.foreach_support(f, A::MyMatrix)

It must emit each nonzero entry exactly once, skip stored zeros, ignore callback return values, and return nothing.

See also: foreach_support_sym.

source
MatrixCovers.foreach_support_sym — Method
foreach_support_sym(f, A)

Symmetric counterpart of foreach_support: call f(i, j, v) once per nonzero unordered pair in canonical order i <= j, including the diagonal.

A must be square and abs.(A) symmetric. Public symmetric solvers check this precondition before calling the traversal.

Objective weighting

For full-grid objective weighting, use multiplicity 1 on the diagonal and 2 off-diagonal. Constraints need no multiplicity.

Extending

To support a new matrix type, define

MatrixCovers.foreach_support_sym(f, A::MyMatrix)

It must emit each nonzero pair once in canonical order, skip zero pairs, ignore callback return values, and return nothing. Triangular storage must map lower entries back to (j, i).

See also: foreach_support.

source
MatrixCovers.gramcover! — Method
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.

source
MatrixCovers.gramcover — Method
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.

No cover of A is required. 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. One pass over the stored nonzeros of A.
  • 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. 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. 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.

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
source
MatrixCovers.initialize_cover! — Method
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.

source
MatrixCovers.initialize_cover — Method
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. It forwards maxiter, cgiter, and start, and ignores feasible.

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.

source
MatrixCovers.initialize_symcover — Method
a = initialize_symcover(A; strategy=:hardcover, feasible=:inflate, kwargs...)

Build a symmetric starting point for symcover_min or soft_symcover. Strategies depend only on A:

strategy names 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 the most-underweighted support entry omitted. It fails if removing that entry empties a row.
  • :diagfeasible — a cover grown from the diagonal by nearest-neighbor propagation.
  • :hardcover — the result of symcover. It forwards maxiter and ignores feasible.

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.

source
MatrixCovers.iscover — Method
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.

cover, symcover, and native AbsLog{2} minimal covers certify their results at zero tolerance. Initializers and extension solvers may require a nonzero rtol.

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; rtol=8eps())
true

julia> iscover([1.0, 1.0], [1.0, 1.0], A)   # a*b' = ones, which does not reach A[2,2]
false
source
MatrixCovers.iterations — Function
iterations(stats::AbstractCoverStats) -> Int

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

source
MatrixCovers.ncomponents — Method
ncomponents(sc::SupportComponents) -> Int

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

source
MatrixCovers.rowcomponent — Method
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.

source
MatrixCovers.scalar_type — Method
MatrixCovers.scalar_type(T)

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

source
MatrixCovers.soft_cover! — Function
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!.

source
MatrixCovers.soft_cover — Method
a, b, stats = soft_cover(ϕ, A; kwargs...)
a, b, stats = soft_cover(A; kwargs...)

Minimize ∑ ϕ(|A[i,j]|/(a[i]*b[j])) without a hard coverage constraint. This is the asymmetric form of soft_symcover. 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, computed by alternating power-mean updates of the row and column scales (Sinkhorn scaling of abs.(A).^p to row sums n_i and column sums m_j, the nonzero counts), with adaptive overrelaxation, followed when the sweeps are slow by Newton steps with backtracking on the objective. 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. It computes in at least Float64, and tol defaults to 4096*eps of that type; a warning reports an iteration that ends without reaching it. maxiter (default 10_000) bounds the sweeps, each an update of all rows and then all columns. The sweeps hand over to Newton when their rate 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. When abs.(A) is exactly symmetric, the minimizer has b == a under the balance convention, and soft_cover instead 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.
  • 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}(): the objective is not convex. Several starts are each refined to a local minimum, and the one with the least objective is returned; requires JuMP and Ipopt. 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.

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
source
MatrixCovers.soft_symcover! — Function
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!.

source
MatrixCovers.soft_symcover — Method
a, stats = soft_symcover(ϕ, A; kwargs...)
a, stats = soft_symcover(A; kwargs...)

Minimize ∑ ϕ(|A[i,j]|/(a[i]*a[j])) for symmetric A, without a hard coverage constraint. The default penalty is PowerMean{2}().

Supported penalties and their keywords are:

  • PowerMean{p}() (default): the convex minimum, computed by 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, followed when they are slow by Newton steps with backtracking on the objective. 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. It computes in at least Float64, and tol defaults to 4096*eps of that type, and a warning reports an iteration that 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}(): the objective is not convex. Several starts are each refined to a local minimum, and the one with the least objective is returned; requires JuMP and Ipopt. 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)   # each row's 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.
source
MatrixCovers.symcover! — Method
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.

source
MatrixCovers.symcover — Method
a = symcover(ϕ, A; maxiter=3)
a = symcover(A; maxiter=3)

Given a square matrix A assumed to be symmetric, return a vector a representing a symmetric hard cover of A: a[i] * a[j] >= abs(A[i, j]) for all i, j.

The method initializes from 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])), raises the scales until every entry is covered, then applies maxiter tightening iterations.

ϕ is accepted for API compatibility but is currently ignored. For a cover that provably minimizes a given ϕ, use symcover_min.

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 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].

source
MatrixCovers.symcover_min — Function
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.

AbsLog is convex in the log scales. If AbsLog{1} has multiple minima, the method selects the one 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.

For Float64, :lsqr uses a Cholesky preconditioner when its predicted storage does not exceed fillbudget bytes (default 2^30) and its predicted numeric- factorization flop count does not exceed flopbudget (default 8e3) times the number of stored entries in A's support, where a symmetric off-diagonal entry counts twice, as nnz counts it for a fully stored matrix. Otherwise, and for every other element type, it preconditions with a maximum-weight spanning forest of the couplings among entries whose constraints are active, together with the diagonal contributions of all entries; that factorization has no fill and needs no budget. fillbudget=Inf and flopbudget=Inf together always select the full Cholesky preconditioner; fillbudget=0 or flopbudget=0 always select the diagonal one. The same budgets apply to the active-set method that finishes an unconverged solve: 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!.

source
MatrixCovers.symcover_min! — Function
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!.

source