Cover objectives

Penalty functions

Any hard cover is valid, but tighter covers better capture the scaling of A. The tightness of a cover 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 but finite and continuous 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

To override a solver's default penalty, pass a penalty as its first argument.

Functions and their objectives

For hard covers, symcover and cover are fast heuristics; their _min counterparts minimize the selected objective. cover_transversal minimizes the AbsLog{2} objective among the hard covers that are tight on a maximum-product transversal (see Transversals and conditioning). A soft cover has no constraint to satisfy, so there is no separate heuristic: soft_symcover and soft_cover always minimize the objective.

The method used to minimize the objective depends on the penalty:

• 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 a single linear least-squares solve for soft covers; for hard covers, see Objective-minimal covers.
• 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. The native iterative solvers cost roughly $O(mn)$ per iteration.

Objective-minimal covers

symcover_min and cover_min minimize the chosen penalty subject to the hard constraint. The built-in AbsLog{2} solver uses an augmented-Lagrangian iteration, with Newton steps between multiplier updates. For the matrix below, the heuristic cover is already close to minimal: the minimal cover lowers the objective by less than 1%.

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, and are available once those packages are loaded:

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.

Transversals and conditioning

Minimizing a penalty does not ensure that the cover-normalized matrix Ahat = A ./ (a .* b') is well conditioned. What matters for conditioning is whether the cover is tight on a suitable transversal, which plays the role for a general matrix that the diagonal plays for a symmetric positive-definite one.

For a symmetric positive-definite matrix A, a * a' covers A exactly when a[i] ≥ sqrt(A[i,i]) for every i. The diagonal entries require this, and it suffices because abs(A[i,j]) ≤ sqrt(A[i,i] * A[j,j]). The minimal symmetric cover is therefore a = sqrt.(diag(A)), which is tight on the whole diagonal and gives Ahat a unit diagonal; symcover and symcover_min both return it, to within roundoff. A general square matrix has no distinguished diagonal, but transversals provide an analog that matters for the asymmetric covers computed by cover, cover_min, and cover_transversal.

A transversal of an n×n matrix is a set of n entries with exactly one in each row and each column: A[1,σ[1]], A[2,σ[2]], …, A[n,σ[n]] for some permutation σ. Permuting the columns turns any transversal into the main diagonal, since A[:, σ][i, i] == A[i, σ[i]], so a transversal is a diagonal up to the ordering of the columns; some authors call it a diagonal or generalized diagonal. "Transversal" is the usual term in the literature on sparse direct solvers. It comes from combinatorics, where a transversal of a collection of sets picks a distinct element from each set; here each row contributes the set of columns of its nonzero entries. The products prod(A[i,σ[i]] for i in 1:n) are, up to sign, the terms of det(A), so a nonsingular matrix has at least one transversal of nonzero entries. A matrix with such a transversal is structurally nonsingular. A maximum-product transversal is one for which prod(abs(A[i,σ[i]]) for i in 1:n) is largest. By the inequality above, the main diagonal of a symmetric positive-definite matrix is a maximum-product transversal.

Transversals limit how well conditioned Ahat can be. If every transversal of Ahat contains an entry of magnitude at most t, then the smallest singular value of Ahat is at most (n+1)/2 * t, however well its rows and columns are otherwise balanced.

For any hard cover, prod(a) * prod(b) is at least the product of the magnitudes on any transversal of A, with equality exactly when abs(Ahat[i,σ[i]]) == 1 on every entry of that transversal. Since det(Ahat) = det(A) / (prod(a) * prod(b)), the hard covers that are tight on a maximum-product transversal are those that maximize abs(det(Ahat)). Such scalings are a standard preprocessing step for Gaussian elimination, where they favor the maximum-product transversal as the pivots (M. Olschowka and A. Neumaier, Linear Algebra Appl. 240:131–151, 1996; I. S. Duff and J. Koster, SIAM J. Matrix Anal. Appl. 22:973–996, 2001).

A cover that minimizes a penalty need not be tight on any transversal. The matrix of Conditioning and the error in A \ y has a single transversal of nonzero entries, A[1,2], A[2,3], and A[3,1], which includes the small entry A[1,2] = η. The AbsLog{2}-minimal cover does not rescale A at all, so Ahat keeps the condition number of A. In contrast, cover_transversal and the default start of cover are tight on the transversal:

julia> using MatrixCovers, LinearAlgebra

julia> η = 1e-8;

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

julia> diag(A[:, [2, 3, 1]])   # the transversal as a diagonal
3-element Vector{Float64}:
 1.0e-8
 1.0
 1.0

julia> scaledcond(a, b) = round(cond(A ./ (a .* b')); sigdigits=3);

julia> a, b = cover_min(A); scaledcond(a, b)
2.83e8

julia> a, b = cover_transversal(A); scaledcond(a, b)
4.05

julia> a, b = cover(A); scaledcond(a, b)
4.05

For a square, structurally nonsingular matrix, cover is tight on a transversal by default (start=:transversal). The :covariant and :geomean starts, and the _min solvers, do not guarantee this. Tightness on a transversal prevents only this one cause of ill-conditioning, and it comes at the cost of a somewhat larger objective than that of the minimal cover.