diff --git a/README.md b/README.md index c57edc4..af6dad9 100644 --- a/README.md +++ b/README.md @@ -16,12 +16,10 @@ for all `i` and `j`, this simple construct finds applications that range from of data to the design of well-behaved numerical algorithms (thanks, e.g., to [bounds on $`\hat A`$'s eigenvalues](https://en.wikipedia.org/wiki/Gershgorin_circle_theorem)). -Fast O(mn) heuristics (`symcover`, `cover`) are provided for everyday use, along -with *soft* covers (`soft_symcover`, `soft_cover`) that penalize under-coverage -rather than forbid it. Objective-minimal hard covers (`symcover_min`, -`cover_min`) minimize a penalty subject to the coverage constraint: the default -squared-log-excess penalty is solved natively, with no external solver, while -the other penalties are available when JuMP and HiGHS (or Ipopt) are loaded. +The package provides O(mn) heuristics (`symcover`, `cover`), *soft* covers that +penalize violations (`soft_symcover`, `soft_cover`), and objective-minimal hard +covers (`symcover_min`, `cover_min`). The default squared-log penalty uses a +built-in solver; other penalties use JuMP with HiGHS or Ipopt. ## Example @@ -44,7 +42,7 @@ julia> iscover(a, A) true ``` -Covers are scale-covariant: rescaling the matrix rescales the cover the same way. +Covers are scale-covariant: ```julia julia> D = Diagonal([10.0, 0.5]); @@ -53,9 +51,8 @@ julia> symcover(D * A * D) ≈ D * a true ``` -Non-symmetric matrices get separate row and column scales from `cover`, and -`cover_min`/`symcover_min` trade the fast heuristic for a cover that minimizes a -penalty subject to the same constraint: +For nonsymmetric matrices, `cover` returns separate row and column scales. +`cover_min` minimizes a penalty subject to the coverage constraint: ```julia julia> M = [1.0 2.0 3.0; 6.0 5.0 4.0]; @@ -65,11 +62,11 @@ julia> a, b = cover(M); julia> iscover(a, b, M) true -julia> aq, bq = cover_min(AbsLog{2}(), M); # minimal, solved natively +julia> aq, bq = cover_min(AbsLog{2}(), M); julia> cover_objective(AbsLog{2}(), aq, bq, M) <= cover_objective(AbsLog{2}(), a, b, M) true ``` -See the [documentation](https://HolyLab.github.io/MatrixCovers.jl/dev/) -for motivation, examples, and a full API reference. +See the [documentation](https://HolyLab.github.io/MatrixCovers.jl/dev/) for the +algorithm guide and API reference. diff --git a/docs/src/index.md b/docs/src/index.md index 07ae7d3..46bb182 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -4,31 +4,25 @@ CurrentModule = MatrixCovers # MatrixCovers -This package computes *covers* of matrices. Given a matrix `A`, a cover (more -specifically, a *hard cover*) is a matrix `C` that can be defined as -`C = a * b'`, where `a` and `b` are non-negative vectors. `C` must satisfy +Given a matrix `A`, a *hard cover* is `C = a * b'`, where `a` and `b` are +nonnegative vectors satisfying ```math C_{ij} \;\geq\; |A_{ij}| \quad \text{for all } i, j. ``` -For a symmetric matrix the cover is symmetric (`b = a`), so a single vector -suffices: `a[i] * a[j] >= abs(A[i, j])`. - -A *minimal cover* chooses `C` "as tight as possible" in bounding `A`, by -criteria that will be described below. - -The package also supports *soft covers*, which penalize uncovered entries instead -of requiring every inequality to hold. +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 a natural **scale-covariant** "summary" of a matrix. If you +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`. Scalar summaries like `norm(A)` or -`maximum(abs, A)` do not have this property and therefore implicitly encode an -arbitrary choice of units. +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 @@ -54,12 +48,8 @@ julia> round.(typeof.(a), a; digits=6) 0.001 N^-1 ``` -`A[i,j]` has units `1/(u[i]*u[j])`, as in a Hessian whose parameters have units -`u[i]`. The cover has units `1/u[i]` and identifies scales of 1 mm, 1 m/s, and -1 kN. Had we expressed `A` in those units, we would have gotten the equivalent cover. - -Normalizing by the cover cancels the units along with the magnitudes, leaving a -matrix that is all-ones, dimensionless, and scale-invariant: +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: ```jldoctest coverunits julia> round.(A ./ (a .* a'); digits=6) @@ -73,9 +63,7 @@ An entry is 1 only where the cover bound is tight, and this is not guaranteed fo 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])`, and a matrix that fails this is rejected with a -`DimensionMismatch`. Without it, the terms in a row of `A*x` have incommensurate -units and cannot be added, so `A*x` is undefined for every `x`. +`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 @@ -91,17 +79,12 @@ A **penalty function** `ϕ` combines those ratios into a scalar objective \sum_{i,j} \phi\!\left(\frac{|A_{ij}|}{a_i\, b_j}\right), ``` -which the solvers minimize. Two penalty families are provided: +Two penalty families are provided: -- [`AbsLog`](@ref)`{p}` — `ϕ(r) = |log r|^p` (and `ϕ(0) = 0`). Convex in log - space, which makes them a favorable (and therefore default) penalty for *hard* - covers, where `r ≤ 1` and `|log r|` is the log-excess of a constraint. - `AbsLog{1}` sums the log-excesses (L1), `AbsLog{2}` sums their squares (L2). - Their principal disadvantage is the divergence and discontinuity at `r = 0`. -- [`AbsLinear`](@ref)`{p}` — `ϕ(r) = |1 - r|^p`. Non-convex, but unlike - `AbsLog` these are finite and continuous at `r = 0` (`ϕ(0) = 1`), so zero entries of `A` - contribute a bounded penalty. This is the penalty used by default for the - *soft* covers, where `r > 1` (an uncovered entry) is allowed but penalized. +- [`AbsLog`](@ref)`{p}`: `ϕ(r) = |log r|^p`, with `ϕ(0) = 0`. It is convex in + log space and is the default for hard covers. +- [`AbsLinear`](@ref)`{p}`: `ϕ(r) = |1-r|^p`. It is nonconvex, finite at zero, + and is the default for soft covers. [`cover_objective`](@ref) evaluates either penalty for a given cover: @@ -122,9 +105,7 @@ julia> cover_objective(AbsLog{2}(), a, A) # sum of squared log-excesses (L2) 3.843624111345611 ``` -Both objectives are zero if and only if every constraint is exactly tight. - -You can override the default penalty by supplying it as an argument to the solvers. +Pass a penalty as the first solver argument to override the default. ## Choosing a cover algorithm @@ -139,9 +120,8 @@ You can override the default penalty by supplying it as an argument to the solve | [`soft_symcover_min`](@ref) | yes | soft (penalized) | `AbsLog{2}`, `AbsLinear` | native for `AbsLog{2}`; else JuMP | | [`soft_cover_min`](@ref) | no | soft (penalized) | `AbsLog{2}`, `AbsLinear` | native for `AbsLog{2}`; else JuMP | -For hard covers, [`symcover`](@ref) trades optimality for speed while -[`symcover_min`](@ref) minimizes the selected objective. The soft solvers differ -in their convergence guarantees: +For hard covers, [`symcover`](@ref) and [`cover`](@ref) are fast heuristics; +their `_min` counterparts minimize the selected objective. For soft covers: - [`soft_symcover`](@ref) and [`soft_cover`](@ref) use native coordinate descent and multistart. For nonconvex or nonsmooth penalties they may stop at a fixed @@ -150,56 +130,48 @@ in their convergence guarantees: `AbsLog{2}` is native; `AbsLinear` requires JuMP and Ipopt; `AbsLog{1}` is not implemented. -Under `AbsLog{2}` the objective is convex with a unique minimizer, so both soft -tiers return the same result. Under `AbsLog{1}`, only the native coordinate-descent -solvers are available. - -The heuristic solvers run in ``O(mn)`` time for an ``m\times n`` matrix and -often come within a few percent of the minimum objective. Native iterative -solvers are also roughly ``O(mn)`` per iteration. JuMP-based methods are slower -and primarily useful when the native methods do not support the selected -penalty. +Under `AbsLog{2}`, the objective is convex, so both soft tiers reach the same minimum. +The heuristics cost ``O(mn)``; native iterative solvers cost roughly ``O(mn)`` +per iteration. ### Covariance of the heuristics -[`symcover`](@ref) and [`cover`](@ref), the two heuristic solvers, *are not universally covariant*. Both -are covariant when every row and column of -`A` has the same pattern of nonzeros, notably for any dense `A` lacking zero entries. But on an *irregular* sparse support they are only -approximately covariant. A symmetric three-node path is enough to show it: +The heuristic solvers are exactly covariant when every row and column has the +same nonzero pattern, including dense matrices without zeros. On irregular +sparse support they may be only approximately covariant: ```jldoctest julia> using MatrixCovers, LinearAlgebra -julia> A = [1.0 1 0; 1 1 1; 0 1 1]; # rows 1 and 3 supported on 2 columns, row 2 on all 3 +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'; # does scaling commute with cover-computation? +julia> P1 = (d .* a1) * (d .* a1)'; P2 = a2 * a2'; -julia> round.(extrema(P2 ./ P1); digits=3) # not for the heuristic solver +julia> round.(extrema(P2 ./ P1); digits=3) (1.0, 1.077) ``` -Both covers are valid, but the difference matters when covariance is required. Use -[`symcover_min`](@ref) or [`cover_min`](@ref), whose minimizer is scale-covariant +Use [`symcover_min`](@ref) or [`cover_min`](@ref) when exact covariance is +required. ### Objective-minimal covers -[`symcover_min`](@ref) and [`cover_min`](@ref) return a cover that minimizes the -chosen penalty subject to the hard constraint. For the default `AbsLog{2}` -penalty they are solved natively (no external solver) by penalty-continuation -with a damped semismooth Newton iteration: +[`symcover_min`](@ref) and [`cover_min`](@ref) minimize the chosen penalty +subject to the hard constraint. The built-in `AbsLog{2}` solver uses penalty +continuation with a damped semismooth Newton iteration: ```jldoctest qmin; filter = r"(\d+\.\d{6})\d+" => s"\1" julia> using MatrixCovers julia> A = [1 2 3; 6 5 4]; -julia> a, b = cover(A); # fast heuristic +julia> a, b = cover(A); -julia> aq, bq = cover_min(AbsLog{2}(), A); # AbsLog{2}-minimal, native +julia> aq, bq = cover_min(AbsLog{2}(), A); julia> a * b' 2×3 Matrix{Float64}: @@ -218,17 +190,15 @@ julia> round(cover_objective(AbsLog{2}(), aq, bq, A); digits=6) 1.141281 ``` -The native solver typically has relative objective error of a few -``\times 10^{-7}``, growing slowly with problem size. The other penalties — `AbsLog{1}` -(a linear program) and the non-convex `AbsLinear` variants — are solved through -[JuMP](https://jump.dev/) and are loaded on demand as a package extension: +`AbsLog{1}` and `AbsLinear` use [JuMP](https://jump.dev/) with HiGHS and Ipopt, +respectively: ```jldoctest jumpmin -julia> using MatrixCovers, JuMP, HiGHS # HiGHS for the AbsLog penalties +julia> using MatrixCovers, JuMP, HiGHS -julia> S = [4 1 0; 1 1 5; 0 5 2]; # symmetric +julia> S = [4 1 0; 1 1 5; 0 5 2]; -julia> round.(symcover_min(AbsLog{1}(), S); digits=6) # L1-minimal symmetric hard cover +julia> round.(symcover_min(AbsLog{1}(), S); digits=6) 3-element Vector{Float64}: 2.0 1.0 @@ -236,94 +206,70 @@ julia> round.(symcover_min(AbsLog{1}(), S); digits=6) # L1-minimal symmetric h julia> A = [1 2 3; 6 5 4]; -julia> a, b = cover_min(AbsLog{1}(), A); # L1-minimal general hard cover +julia> a, b = cover_min(AbsLog{1}(), A); -julia> round.(a * b'; digits=6) # tight on four of the six entries +julia> round.(a * b'; digits=6) 2×3 Matrix{Float64}: 2.4 2.0 3.0 6.0 5.0 7.5 ``` -The solver returns values good to roughly solver tolerance, so these examples -round before displaying. - [`soft_symcover_min`](@ref) and [`soft_cover_min`](@ref) solve `AbsLog{2}` natively and use JuMP with Ipopt for `AbsLinear`. They do not accept -`AbsLog{1}`; the soft `AbsLog{1}` covers are available through [`soft_symcover`](@ref) and -[`soft_cover`](@ref), which are native. +`AbsLog{1}`; use [`soft_symcover`](@ref) or [`soft_cover`](@ref) instead. ### Uniqueness -The `AbsLog{2}()` penalty generally has a unique minimum, with one exception: -row/column scaling `a → γ*a`, `b → b/γ` does not affect `C` and is thus -invisible to the objective function. For non-symmetric (i.e., not `symcover`) problems, -the scaling of each is pinned by the balance convention +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`, respectively. The gauge freedom, and hence this -convention, acts independently on each connected component of the bipartite -support graph of `A` (rows and columns as vertices, stored nonzeros as edges), -so the sums are taken within each component separately. This convention is not -scale-invariant but has no impact on the cover itself. +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. -The `AbsLog{1}()` minimum can be a face of the feasible polytope containing -different covers with the same objective. The implementation returns the member -that also minimizes the `AbsLog{2}` objective. - -`AbsLinear` penalties typically have isolated minima, so are not as degenerate -as `AbsLog{1}()`, but these minima occur in separate basins. There is no -guarantee of global optimality. +`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. 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 the starting -point. The interface has three layers: - -- **Initializers** construct starting points. [`initialize_symcover`](@ref) and - [`initialize_cover`](@ref) take a `strategy` — `:geomean`, `:leaveout`, `:diagfeasible`, - or `:hardcover` — and return that point. Each is a property of `A` alone; no objective - is involved, so an initializer takes no penalty. A second keyword, `feasible`, - controls how the point is made to cover `A`: `:inflate` (the default) applies one - common factor, `:boost` raises only the rows touching a violated entry, and `:none` - leaves it unchanged. Hard-cover solvers use one of the first two; soft-cover - solvers use `:none`. - - The two feasible routes reach different points on the boundary and can enter - different basins. The heuristic [`cover`](@ref) computes the geometric mean, - boosts it, and then tightens it. -- **Refiners** improve a starting point in place, and are the `!`-suffixed forms of the +For objectives with multiple minima, the result can depend on its starting +point. The interface separates initialization, refinement, and multistart +selection: + +- **Initializers** [`initialize_symcover`](@ref) and [`initialize_cover`](@ref) + build a named `strategy`. Their `feasible` keyword selects uniform inflation, + selective boosting, or no feasibility step. +- **Refiners** are the `!`-suffixed forms of the solvers: [`symcover_min!`](@ref), [`cover_min!`](@ref), [`soft_symcover!`](@ref), [`soft_cover!`](@ref), [`soft_symcover_min!`](@ref), and [`soft_cover_min!`](@ref) - validate a caller-provided start, then optimize from it. Hard refiners require - a cover; soft refiners accept starts built with `feasible=:none`. -- **Solvers** bundle the two. [`symcover_min`](@ref), [`cover_min`](@ref), + optimize a supplied point. Hard refiners require a cover; soft refiners do not. +- **Solvers** [`symcover_min`](@ref), [`cover_min`](@ref), [`soft_symcover`](@ref), [`soft_cover`](@ref), [`soft_symcover_min`](@ref), and - [`soft_cover_min`](@ref) refine a set of starts (the `strategies` keyword or the - multistart list) and return the best result by [`cover_objective`](@ref). + [`soft_cover_min`](@ref) refine several starts and return the best objective. -In general, a plain form chooses among several starts, while a `!` form refines -the supplied start. [`symcover!`](@ref) and [`cover!`](@ref) are initializers and -overwrite their vector arguments. - -For finer control, you can run these manually: +Plain forms choose their starts; `!` forms refine the supplied start, except +[`symcover!`](@ref) and [`cover!`](@ref), which are in-place heuristics. ```jldoctest manualstart -julia> using MatrixCovers, JuMP, Ipopt # Ipopt for the AbsLinear penalties +julia> using MatrixCovers, JuMP, Ipopt julia> S = [4 1 0; 1 1 5; 0 5 2]; -julia> round.(symcover_min(AbsLinear{2}(), S); digits=6) # use all default starts +julia> round.(symcover_min(AbsLinear{2}(), S); digits=6) 3-element Vector{Float64}: 2.0 1.0 5.0 -julia> round.(symcover_min(AbsLinear{2}(), S; strategies=(:geomean,)); digits=6) # use one strategy +julia> round.(symcover_min(AbsLinear{2}(), S; strategies=(:geomean,)); digits=6) 3-element Vector{Float64}: 2.0 1.0 5.0 -julia> a0 = initialize_symcover(S; strategy=:geomean); # construct a start explicitly +julia> a0 = initialize_symcover(S; strategy=:geomean); julia> symcover_min!(AbsLinear{2}(), a0, S); @@ -334,66 +280,43 @@ julia> round.(a0; digits=6) 5.0 ``` -The same menu supplies the starting points of the [`soft_symcover`](@ref) and -[`soft_cover`](@ref) multistarts, adding (by default) a few randomized -perturbations of a base point up to a user-controllable number of `starts`. - -For the convex `AbsLog` penalties the start cannot change the result, and the refiners -accept one only so that the two families share an interface. +For convex `AbsLog` penalties, the start does not change the result. ## Consuming one factor alone: gauges and Gram covers -For asymmetric covers, only the products `a[i]*b[j]` are determined by the -problem; the split into the pair is fixed by the balance convention described -under [Uniqueness](@ref). That convention makes the split *deterministic*, but -it is still a convention, and it is **not covariant** under one-sided -rescaling: if `a*b'` covers `A`, then `a*(D*b)'` covers `A*D` — but the -balanced representative of the rescaled problem is `(γ*a, D*b/γ)` for a -per-component constant `γ ≠ 1` that depends on `D`. +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`. -This matters when covers are composed. For example, the -[Levenberg-Marquardt algorithm](https://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm) -uses products `J'*J` of a Jacobian `J`. If `a*b'` covers `J`, then -`(a'*a) * b * b'` covers `J'*J`, but its tightness and invariance depends on the balance -convention for `a` and `b`. - -To do better, this package provides the Gram cover `s = ` [`gramcover`](@ref)`(a, b, J[, W])`, -a symmetric cover of `J'*W*J` built from the asymmetric cover of `J`. -Built this way, `s` covaries with right-scaling of `J`. +[`gramcover`](@ref)`(a, b, J[, W])` constructs a symmetric cover of `J'*W*J` +that covaries with right-scaling of `J`. ```jldoctest gauge julia> using MatrixCovers, LinearAlgebra julia> J = [1.0 2; 3 4; 5 6]; -julia> D = Diagonal([100.0, 1.0]); # reparametrize the second frame +julia> D = Diagonal([100.0, 1.0]); -julia> a1, b1 = cover(J); a2, b2 = cover(J * D); # `cover` is covariant because J has no zeros; `cover_min` is safer +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 # bare-factor consumers see this constant +julia> first(r) ≈ 1 false julia> s1 = gramcover(a1, b1, J); s2 = gramcover(a2, b2, J * D); -julia> s2 ≈ D.diag .* s1 # the Gram cover co-varies exactly +julia> s2 ≈ D.diag .* s1 true ``` ## Worked example: roundoff in `A \ b` -A cover provides scales for each variable with which to measure a solution independently -of the parameterization. - -Solving `x = A \ b` is *contravariant*: rescaling `A → D*A*D` and `b → D*b` sends -`x → x ./ d`, while the cover is covariant, `a → d .* a`. The products `x .* a` are -therefore unchanged, and `∑ᵢ |xᵢ * aᵢ|` is a measure of the solution's size that is the -same in every frame. - -That quantity can be estimated from the magnitudes of `A` and `b` alone, without -forming `x` at all: +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`: ```jldoctest roundoff julia> using MatrixCovers, LinearAlgebra @@ -413,8 +336,7 @@ julia> mag = sum(abs(bi / ai) for (bi, ai) in zip(b, a)) 4.5 ``` -The cover reports natural scales of 1000 and 2, and `mag` estimates the size of the -solution measured against them — here within a factor of 1.5 of the truth: +Here `mag` is within a factor of 1.5 of the scaled solution norm: ```jldoctest roundoff julia> x = A \ b; @@ -423,8 +345,7 @@ julia> sum(abs.(x .* a)) 3.0 ``` -Both numbers are scale-invariant, so the estimate is unchanged by any diagonal -rescaling of the problem: +The estimate is unchanged by diagonal rescaling: ```jldoctest roundoff julia> d = [0.05, 3.0]; @@ -437,8 +358,7 @@ julia> sum(abs(bi / ai) for (bi, ai) in zip(bd, ad)) 4.5 ``` -This makes `eps(mag)` a scale-invariant estimate of the roundoff floor of the sum. -For a well-conditioned `A`, the error of the `Float64` solve meets that floor: +For well-conditioned `A`, `eps(mag)` estimates the roundoff floor: ```jldoctest roundoff julia> xbig = big.(A) \ big.(b); @@ -447,9 +367,7 @@ julia> abs(sum(abs.(x .* a)) - sum(abs.(Float64.(xbig) .* a))) <= 2 * eps(mag) true ``` -The estimate is built from magnitudes only, so it knows nothing about the conditioning -of `A` or about cancellation during the solve. When `A` is ill-conditioned the true -error sits far above the floor: +For ill-conditioned `A`, the error can be much larger: ```jldoctest roundoff julia> Aill = [1.0 -0.9999; -0.9999 1.0]; @@ -470,8 +388,7 @@ julia> err > 1e6 * eps(magill) true ``` -Folding in the condition number of the *normalized* matrix `A ./ (a .* a')` — itself -scale-invariant, since normalizing cancels the frame — restores a usable bound: +The condition number of the normalized matrix provides a corresponding bound: ```jldoctest roundoff julia> κ = cond(Aill ./ (aill .* aill')); diff --git a/ext/MatrixCoversIpoptExt.jl b/ext/MatrixCoversIpoptExt.jl index bf47724..5894a5f 100644 --- a/ext/MatrixCoversIpoptExt.jl +++ b/ext/MatrixCoversIpoptExt.jl @@ -6,27 +6,15 @@ using MatrixCovers using MatrixCovers: AbsLinear using MatrixCovers: _edge_list, _sym_edge_list, _degrees -# The models are built over 1-based positions 1:n; `pr`/`pc` map each position to the -# corresponding axis index of `A`, and results are scattered back onto vectors -# whose axes match `A`'s so offset axes are honored. A row/column of `A` with no -# nonzero entry appears in no constraint or objective term; its scale is set to -# exactly 0, matching the native solvers. -# -# `A` is read through the support hook and gathered into a flat edge list in position -# space — `ei`/`ej` the endpoints, `elog` the log-magnitude — so a model costs O(nnz) -# to build rather than O(length(A)). The symmetric list is the full-grid reading; the -# hard-cover models below sum over its `ei <= ej` half instead, per the objective each -# one is defined by. +# Models use 1-based positions and scatter results back to `A`'s axes. Support is +# gathered as an O(nnz) edge list; unsupported scales are zero. -# Ipopt returns a local minimum selected by the start. These kernels therefore -# implement the mutating refiners; the main package supplies multistart drivers. +# Ipopt kernels refine one start; the main package supplies multistart drivers. check_solved(model, fname) = MatrixCovers.check_solved(JuMP.termination_status(model), "Ipopt", fname) -# `set_silent` alone still lets Ipopt print its startup banner, once per session, from -# its C++ core; `sb` ("suppress banner") is the option that covers it. Every model here -# is solved for a caller who asked for a cover, not for solver output. +# Suppress both solver output and Ipopt's startup banner. function _ipopt_model() model = JuMP.Model(Ipopt.Optimizer) JuMP.set_silent(model) @@ -34,10 +22,7 @@ function _ipopt_model() return model end -# The `i ≤ j` half of a symmetric gather, paired with the multiplicity `w` each entry -# stands for in the full grid: an off-diagonal pair is two entries of `A`, a diagonal -# entry one. Carrying the weight is equivalent to summing over both orientations and -# costs half the terms — and, for the AbsLinear{1} models, half the auxiliary variables. +# Symmetric triangle with full-grid multiplicities. function _triangle(fi, fj, flog) keep = [e for e in eachindex(fi) if fi[e] <= fj[e]] return fi[keep], fj[keep], flog[keep], [fi[e] == fj[e] ? 1 : 2 for e in keep] @@ -112,16 +97,8 @@ end # ============================================================ # Hard cover: cover_min!(::AbsLinear{p}, a, b, A) -# The bipartite analog of symcover_min!: row scales α = log a, column scales -# β = log b, residuals over every stored (i, j) rather than over i ≤ j. The product -# a[i]*b[j] is invariant under (α, β) → (α + s, β - s), so — unlike the symmetric -# problem, which has no such freedom — the model is degenerate along that direction -# until the balance constraint ∑ nzaᵢ αᵢ = ∑ nzbⱼ βⱼ pins it, exactly as -# cover_min(::AbsLog{1}) does. That constraint pins only the global gauge direction; -# a support with more than one connected component carries one such gauge per -# component (see MatrixCovers._support_components), so each kernel below finishes -# with a post-solve balance shift (`_balance_cover!`, then `inflate_feasible!` to -# restore exact coverage) that pins the rest. +# Asymmetric hard-cover model in row and column log scales. The model pins the +# global gauge; post-processing balances components and restores feasibility. # ============================================================ function MatrixCovers.cover_min!(::AbsLinear{2}, a::AbstractVector, b::AbstractVector, A) @@ -196,9 +173,7 @@ end # ============================================================ # Soft cover: soft_symcover_min!(::AbsLinear{p}, a, A) -# Same objective, no coverage constraints — so the start need not cover `A`, and the -# raw geometric mean (the exact soft AbsLog{2} optimum) is a natural one. The multistart -# driver over these kernels is native; see soft_symcover_min. +# Same objective without coverage constraints; starts need not cover `A`. # ============================================================ function MatrixCovers.soft_symcover_min!(::AbsLinear{2}, a::AbstractVector, A) @@ -255,14 +230,8 @@ end # ============================================================ # Soft cover: soft_cover_min!(::AbsLinear{p}, a, b, A) -# The bipartite analog of soft_symcover_min!, and the unconstrained analog of cover_min!: -# no coverage constraints, but the same row/column gauge, pinned in the model by the same -# balance constraint ∑ nzaᵢ αᵢ = ∑ nzbⱼ βⱼ. As in cover_min!, that constraint pins only the -# global gauge direction, so each kernel below finishes with a post-solve `_balance_cover!` -# that pins the rest (one per connected component of the support); unlike the hard-cover -# kernels, no `inflate_feasible!` follows, since the soft objective imposes no coverage -# constraint for it to restore. A zero entry of `A` contributes ϕ(0) = 1 whatever the -# scales, so the count of zeros enters the objective as a constant, matching cover_objective. +# Asymmetric soft-cover model. Post-processing balances component gauges; zero +# entries contribute the constant `ϕ(0) = 1`. # ============================================================ function MatrixCovers.soft_cover_min!(::AbsLinear{2}, a::AbstractVector, b::AbstractVector, A) diff --git a/ext/MatrixCoversJuMPExt.jl b/ext/MatrixCoversJuMPExt.jl index 05cffa3..20a666e 100644 --- a/ext/MatrixCoversJuMPExt.jl +++ b/ext/MatrixCoversJuMPExt.jl @@ -10,19 +10,10 @@ using LinearAlgebra: dot check_solved(model, fname) = MatrixCovers.check_solved(JuMP.termination_status(model), "HiGHS", fname) -# The models are built over 1-based positions 1:n; `pr`/`pc` map each position to -# the corresponding axis index of `A`, and results are scattered back onto vectors -# whose axes match `A`'s so offset axes are honored. A row/column of `A` with no -# nonzero entry carries no constraint or objective term; its scale is set to exactly -# 0, matching the native solvers. -# -# `A` is read through the support hook and gathered into a flat edge list in position -# space — `ei`/`ej` the endpoints, `elog` the log-magnitude — so a model costs O(nnz) -# to build rather than O(length(A)). The symmetric list is the full-grid reading, whose -# `ei <= ej` half is the constraint set. +# Models use 1-based positions and scatter results back to `A`'s axes. Support is +# gathered as an O(nnz) edge list; unsupported scales are zero. -# Exact reference for the native `symcover_min(::AbsLog{2})`: same QP, solved by -# HiGHS. Not exported; used by the test suite to cross-check the native solver. +# HiGHS reference for tests of native symmetric `AbsLog{2}`. function MatrixCovers.symcover_min_jump(::AbsLog{2}, A) axr = axes(A, 1) axes(A, 2) == axr || throw(ArgumentError("symcover_min_jump requires a square matrix")) @@ -56,16 +47,10 @@ function MatrixCovers.symcover_min!(::AbsLog{1}, a::AbstractVector, A) return a end -# Relative slack allowed on the AbsLog{1} optimum while the AbsLog{2} objective is -# minimized over it. The incumbent attains the bound exactly, so the face is never empty; -# the slack only has to absorb the rounding of re-evaluating the objective row, and it -# bounds how far the reported AbsLog{1} objective can drift above its true optimum. +# Slack for re-evaluating the `AbsLog{1}` optimum during tie-breaking. const LEX_L1_SLACK = 1e-9 -# Select a unique point on the optimal AbsLog{1} face by minimizing AbsLog{2} -# over it. Both objectives depend only on scale-invariant residuals. -# -# `residuals` uses the same support weighting as `cover_objective`. +# Break `AbsLog{1}` ties with `AbsLog{2}` using full-grid support weights. function _minimize_l2_over_l1_face!(model, lin, residuals, fname) isempty(residuals) && return nothing linopt = JuMP.value(lin) @@ -77,11 +62,8 @@ function _minimize_l2_over_l1_face!(model, lin, residuals, fname) return nothing end -# The AbsLog{1} hard cover is an LP: the coverage constraint forces every residual -# α[i]+α[j]-log|A[i,j]| to be nonnegative, so |·| drops away and the objective is -# linear in α. Its optimum is a face, not a point, so a second stage picks the canonical -# member of that face. `start`, when given, is a cover of `A` supplying the initial point; -# it is a hint to the solver, and the canonical selection keeps it out of the result. +# Symmetric `AbsLog{1}` LP. A second stage selects the canonical point on the +# optimal face; `start` is only a solver hint. function _symcover_min_abslog1(A, start) axr = axes(A, 1) axes(A, 2) == axr || throw(ArgumentError("symcover_min requires a square matrix")) @@ -136,10 +118,7 @@ function MatrixCovers.cover_min_jump(::AbsLog{2}, A) @constraint(model, α[ei[e]] + β[ej[e]] - elog[e] >= 0) end nza, nzb = _degrees(ei, m), _degrees(ej, n) - # Pins only the global gauge direction; a disconnected support carries one (e; -e) - # gauge per component (see MatrixCovers._support_components), so the remaining - # directions are left to whichever vertex HiGHS returns. The post-solve balance - # shift below fixes all of them, matching the native solver. + # The post-solve balance handles component gauges not pinned here. @constraint(model, sum(nza[i] * α[i] for i in 1:m) == sum(nzb[j] * β[j] for j in 1:n)) JuMP.optimize!(model) check_solved(model, "cover_min_jump") @@ -165,11 +144,8 @@ function MatrixCovers.cover_min!(::AbsLog{1}, a::AbstractVector, b::AbstractVect return a, b end -# Asymmetric counterpart of `_symcover_min_abslog1`, on the bipartite support: the -# same LP over row scales α and column scales β. The balance constraint below pins -# the global row/column gauge; a support with more than one connected component -# carries additional per-component gauges that the post-solve balance shift pins, -# so the split between `a` and `b` is deterministic. +# Asymmetric `AbsLog{1}` LP. Balance globally in the model and per component +# after solving. function _cover_min_abslog1(A, start) axr = axes(A, 1) axc = axes(A, 2) @@ -199,11 +175,7 @@ function _cover_min_abslog1(A, start) @constraint(model, α[ei[e]] + β[ej[e]] - elog[e] >= 0) end nza, nzb = rowcount, colcount - # Gauge pin: the products a[i]*b[j] are unchanged by a -> c*a, b -> b/c, so without this - # the split between `a` and `b` would be arbitrary. It is orthogonal to the AbsLog{1} - # degeneracy the second stage resolves, and stays in force there. It pins only the - # global gauge direction, one of possibly several (one per connected component of the - # support); the post-solve balance shift below pins the rest. + # Pin the global row/column gauge; post-processing handles components. @constraint(model, sum(nza[i] * α[i] for i in 1:m) == sum(nzb[j] * β[j] for j in 1:n)) JuMP.optimize!(model) check_solved(model, "cover_min") diff --git a/ext/MatrixCoversUnitfulExt.jl b/ext/MatrixCoversUnitfulExt.jl index f63402a..accc42b 100644 --- a/ext/MatrixCoversUnitfulExt.jl +++ b/ext/MatrixCoversUnitfulExt.jl @@ -12,10 +12,7 @@ const MC = MatrixCovers # halves the diagonal's exponents. const UnitExps = Dict{FreeUnits,Rational{Int}} -# A cover objective is dimensionless, so it accumulates in the quantity's own -# numeric type. This is stated for `Quantity{T}` rather than a concrete quantity -# type because a matrix whose entries carry different units has exactly that -# abstract element type. +# Cover objectives are dimensionless and accumulate in the quantity's numeric type. MC.scalar_type(::Type{<:Quantity{T}}) where {T} = MC.scalar_type(T) const QMatrix = AbstractMatrix{<:Quantity} const QVector = AbstractVector{<:Quantity} @@ -24,22 +21,13 @@ const QVector = AbstractVector{<:Quantity} # Unit algebra # ============================================================ -# Unitful exposes no public accessor for the atomic units of a `FreeUnits`, and no -# supported way to iterate them, so the code below reads its representation -# directly: `FreeUnits{N,D,A}` carries `N` as a tuple of `Unit{U,D}`, each with -# fields `tens::Int` and `power::Rational{Int}`. That layout is what the -# `Unitful.Unit` and `Unitful.FreeUnits` docstrings specify. The decomposition -# cannot be replaced by unit arithmetic: `gauge` takes a per-atom median over -# rational exponents, which `*`, `/`, and `^` cannot express. +# Decompose `FreeUnits{N,D,A}` through its documented type parameters. Gauge +# selection needs per-atom rational exponents, which unit arithmetic cannot expose. -# An atomic unit at the first power. A prefix belongs to the atom -- `mm` and `m` -# are distinct -- so a coordinate named in `mm` keeps `mm` in its cover. +# Atomic unit at the first power; prefixes remain distinct atoms. atomic(x::Unit{N,D}) where {N,D} = FreeUnits{(Unit{N,D}(x.tens, 1//1),), D, nothing}() -# The third parameter of `FreeUnits` is the affine offset, `nothing` for an -# ordinary unit. An affine unit measures from a shifted origin, so `a[i]*b[j]` -# does not scale it and there is no cover to find; Unitful likewise refuses to -# multiply affine units. +# Reject affine units because multiplicative covers require an absolute zero. function exps(u::FreeUnits{N,D,A}) where {N,D,A} A === nothing || throw(ArgumentError(""" affine units are not supported: `$u` measures from a shifted origin, so no \ @@ -54,17 +42,14 @@ function exps(u::FreeUnits{N,D,A}) where {N,D,A} return d end -# `ContextUnits` and `FixedUnits` carry a conversion context this code does not -# read, so they are refused by name rather than reaching `atomic` as a -# `MethodError` on an unexported internal. +# Context-dependent units are unsupported. exps(u::Unitful.Units) = throw(ArgumentError( "unsupported unit type $(nameof(typeof(u))) for `$u`: MatrixCovers reads `FreeUnits`. " * "Convert with `uconvert(FreeUnits(u), x)`.")) exps(q::Quantity) = exps(unit(q)) -# Zero exponents are pruned throughout so that `==` on a `UnitExps` compares -# units rather than representations. +# Remove zero exponents so equality compares units, not representations. function combine(f, d1::UnitExps, d2::UnitExps) d = UnitExps() for k in union(keys(d1), keys(d2)) @@ -83,26 +68,20 @@ freeunits(d::UnitExps) = isempty(d) ? Unitful.NoUnits : prod(k^v for (k, v) in d # Rank-1 unit factorization # ============================================================ -# A cover needs `unit(A[i,j]) == unit(a[i])*unit(b[j])`: the unit exponents form a -# rank-1 additive matrix, and any violation is witnessed by a 2x2 minor. The message -# quotes that minor, so it names only entries the caller wrote. +# Unit exponents must form a rank-1 additive matrix. A failing 2×2 minor +# identifies a mismatch. function throw_nofactor(lhs, rhs, lhsname, rhsname) throw(DimensionMismatch(""" units of `A` do not factor: $lhsname = $lhs, but $rhsname = $rhs. A cover requires `unit(A[i,j]) == unit(a[i])*unit(b[j])`, which forces these two \ - products to agree. Any matrix that models the physical world satisfies this: \ - without it the terms of a row of `A*x` do not share units, so `A*x` is undefined \ - for every `x`.""")) + products to agree. Without this factorization, the terms in a row of `A*x` \ + can have incompatible units.""")) end -# A concrete element type names one unit for every entry, structural zeros -# included, so there is nothing to verify and nothing to read: `A` contributes only -# `unit(eltype(A))`. This is the only shape a sparse `A` can take, since sparse -# storage synthesizes its structural zeros with `zero(eltype(A))`. +# A concrete element type gives every entry, including structural zeros, one unit. uniform_unit(A::QMatrix) = isconcretetype(eltype(A)) ? exps(unit(eltype(A))) : nothing -# `unit(a[i])` and `unit(b[j])` up to the gauge `a -> a*c`, `b -> b/c`, taken -# relative to the first row and column. +# Factor row and column units relative to the first row and column. function factor_units(A::QMatrix) ax1, ax2 = axes(A) uas = similar(Array{UnitExps}, ax1) @@ -138,8 +117,7 @@ function factor_units(A::QMatrix) return ua, ub end -# `a[i]*a[i] == A[i,i]` pins `unit(a[i])` outright: the symmetric gauge `a -> a*c` -# would scale every product by `c^2`, so only `c = 1` preserves them. +# Diagonal entries fix symmetric scale units directly. function factor_units_sym(A::QMatrix) ax = axes(A, 1) uas = similar(Array{UnitExps}, ax) @@ -163,23 +141,12 @@ function factor_units_sym(A::QMatrix) return ua end -# The gauge `a -> a*c`, `b -> b/c` leaves every product `a[i]*b[j]` unchanged, so -# the factorization fixes the units only up to `c`. Pin it by minimizing the total -# atomic-unit powers carried by the two scale vectors, +# Choose the unit gauge by minimizing total atomic-unit powers: # # minimize_c ∑_i ‖exps(ua[i]*c)‖₁ + ∑_j ‖exps(ub[j]/c)‖₁, # -# which separates over atoms into subproblems ∑_i |t + αᵢ| + ∑_j |t - βⱼ|, each -# minimized on the median interval of {-αᵢ} ∪ {βⱼ}. -# -# That interval is a single point only when the atom's exponents pin it; otherwise -# the objective is flat across it and the choice within it is the whole content of -# the convention. Take its midpoint, which makes `cover` reproduce `symcover` on -# symmetric input. There, `unit(A[i,j])` has exponents `dᵢ + dⱼ`, so `αᵢ = dᵢ - d₀` -# and `βⱼ = d₀ + dⱼ` relative to the reference row `i0`, and the points -# `{d₀ - dᵢ} ∪ {d₀ + dⱼ}` are distributed symmetrically about `d₀`. The midpoint is -# therefore `d₀` exactly, which is the shift that returns `unit(a[i])` with -# exponents `dᵢ` -- what `a[i]*a[i] == A[i,i]` demands. +# Each atom reduces to a median interval; its midpoint makes `cover` agree with +# `symcover` on symmetric input. function gauge(uas, ubs) atoms = Set{FreeUnits}() for d in uas @@ -211,17 +178,11 @@ end # Strip and reattach # ============================================================ -# `A` is stripped in the units the caller wrote, not in a canonical system. The -# cover itself is scale-invariant, but the balance convention that splits `a` from -# `b` is not, so the strip scale selects the parametrization; the caller's units -# are their statement of the scale they want it pinned to. Stripping requires the -# units to factor as written -- otherwise entries on incommensurate scales -# (`1.0mm^-2` and `1.0m^-2` both strip to `1.0`) would be covered as if comparable. +# Strip values in their written units so the balance convention follows the +# caller's parametrization. Units must factor before stripping. strip_matrix(A::QMatrix) = ustrip.(A) -# The `*_min!` family reads `a` as a start, so its units must be the cover's. Any -# dimensionally equivalent spelling is accepted and converted; `ustrip` raises on -# a start that is not. +# Refiners accept dimensionally equivalent start units and convert them. strip_start(a::QVector, ua) = map(ustrip, ua, a) reattach(a, ua) = a .* ua @@ -238,9 +199,7 @@ function asym(f, A::QMatrix, ϕ...; kwargs...) return reattach(a, ua), reattach(b, ub) end -# `a` is overwritten, so neither its values nor its units are read: the scratch it -# is stripped into takes its element type from `A`, matching what the unitless -# methods allocate. An `a` of undefined references is a valid destination. +# Allocating scratch from `A` permits uninitialized destination vectors. function sym!(f, a::QVector, A::QMatrix, ϕ...; kwargs...) ua = factor_units_sym(A) An = strip_matrix(A) @@ -275,12 +234,7 @@ function asymstart!(f, a::QVector, b::QVector, A::QMatrix, ϕ...; kwargs...) return a, b end -# Every penalty slot below mirrors MatrixCovers's own method table: where it -# accepts any `AbstractCoverPenalty` these do too, and where it dispatches on -# concrete penalties these enumerate the same ones. Each method is then strictly more specific than the -# one it shadows -- including those in the JuMP and Ipopt extensions, which leave the -# matrix slot untyped -- so no ambiguity arises. A penalty the package does not -# support raises a `MethodError` here exactly as it does on a unitless matrix. +# Mirror the core penalty dispatch while specializing on unitful matrices. const PENALTIES = (:(AbsLog{1}), :(AbsLog{2}), :(AbsLinear{1}), :(AbsLinear{2})) # Heuristic covers and initializers: `ϕ` is checked but not consulted. @@ -294,7 +248,7 @@ MC.cover(ϕ::MC.AbstractCoverPenalty, A::QMatrix; kwargs...) = asym(MC.cover, A, MC.cover!(a::QVector, b::QVector, A::QMatrix; kwargs...) = asym!(MC.cover!, a, b, A; kwargs...) MC.cover!(ϕ::MC.AbstractCoverPenalty, a::QVector, b::QVector, A::QMatrix; kwargs...) = asym!(MC.cover!, a, b, A, ϕ; kwargs...) -# `cover`/`cover!` dispatch on `Adjoint`/`Transpose` upstream without an eltype +# Core `cover`/`cover!` methods dispatch on `Adjoint`/`Transpose` without an eltype # bound, so a wrapped `Quantity` matrix needs these to stay unambiguous. for W in (:(LinearAlgebra.Adjoint{<:Quantity}), :(LinearAlgebra.Transpose{<:Quantity})) @eval begin @@ -308,7 +262,7 @@ MC.initialize_symcover!(a::QVector, A::QMatrix; kwargs...) = sym!(MC.initialize_ MC.initialize_cover(A::QMatrix; kwargs...) = asym(MC.initialize_cover, A; kwargs...) MC.initialize_cover!(a::QVector, b::QVector, A::QMatrix; kwargs...) = asym!(MC.initialize_cover!, a, b, A; kwargs...) -# Soft covers and the `*_min` family: `ϕ` is dispatched on upstream. +# Soft covers and the `*_min` family preserve core penalty dispatch. MC.soft_symcover(A::QMatrix; kwargs...) = sym(MC.soft_symcover, A; kwargs...) MC.soft_cover(A::QMatrix; kwargs...) = asym(MC.soft_cover, A; kwargs...) MC.symcover_min(A::QMatrix; kwargs...) = sym(MC.symcover_min, A; kwargs...) @@ -320,13 +274,7 @@ MC.soft_symcover_min!(a::QVector, A::QMatrix; kwargs...) = symstart!(MC.soft_sym MC.soft_cover_min(A::QMatrix; kwargs...) = asym(MC.soft_cover_min, A; kwargs...) MC.soft_cover_min!(a::QVector, b::QVector, A::QMatrix; kwargs...) = asymstart!(MC.soft_cover_min!, a, b, A; kwargs...) -# MatrixCovers types the matrix slot of its sparse refiners, where the methods above -# type the element: neither is more specific for a sparse matrix of quantities, so the -# two are ambiguous there. These methods resolve that pair. They are the only overlap -- -# every other sparse method leaves its matrix slot untyped. -# -# Sparse storage synthesizes structural zeros with `zero(eltype)`, so the element type -# is concrete and every entry carries the same unit. +# Resolve the overlap between sparse refiner and unitful matrix methods. const QSparse = SparseMatrixCSC{<:Quantity} const QSparseSym = Union{QSparse, Symmetric{<:Quantity,<:SparseMatrixCSC}, diff --git a/src/gram_covers.jl b/src/gram_covers.jl index ac33887..0bca308 100644 --- a/src/gram_covers.jl +++ b/src/gram_covers.jl @@ -16,83 +16,43 @@ const GRAMCOVER_DEGENERATE = :error s = gramcover(a, b, sc::SupportComponents, w::AbstractVector) s = gramcover(a, b, sc::SupportComponents, W::AbstractMatrix; degenerate=:$(GRAMCOVER_DEGENERATE)) -Given an asymmetric cover `a[i]*b[j] >= abs(A[i,j])`, return a symmetric cover -`s` of a weighted Gram matrix without forming it: `s[j]*s[k] >= abs(G[j,k])`, -where `G` is `A'*A`, `A'*Diagonal(w)*A`, or `A'*W*A`. Only `abs.(W)` enters the -bound, so `W` need not be symmetric or positive semidefinite. Passing -`W::Diagonal` is equivalent to passing `W.diag`. +Given a cover `(a, b)` of `A`, return a symmetric cover of `A'*A`, +`A'*Diagonal(w)*A`, or `A'*W*A` without forming the product. Only `abs.(W)` +enters the bound. Passing `Diagonal(w)` is equivalent to passing `w`. `(a, b)` must cover `A`; use [`iscover`](@ref)`(a, b, A)` to check it. -The methods accepting [`SupportComponents`](@ref) reuse a previous -[`support_components`](@ref)`(A)` computation. For these methods, `sc.rowax` and -`sc.colax` replace `axes(A, 1)` and `axes(A, 2)` in the axis requirements. +Pass [`SupportComponents`](@ref) to reuse a previous +[`support_components`](@ref)`(A)` computation. -Equivalent componentwise rescalings of `(a, b)` produce the same `s`. Some -coupling patterns make this impossible; the general form then throws an -`ArgumentError`. Set `degenerate=:uniform` to return a gauge-dependent cover. +The result is invariant under componentwise rescaling of `(a, b)`. If `W` makes +this impossible, the matrix form throws an `ArgumentError`; use +`degenerate=:uniform` to allow a gauge-dependent result. # Extended help -For `G[j,k] = Σ_{i,i'} A[i,j]*W[i,i']*A[i',k]`, the triangle inequality against -`a[i]*b[j] >= abs(A[i,j])` gives -`abs(G[j,k]) <= (Σ_{i,i'} a[i]*abs(W[i,i'])*a[i']) * b[j]*b[k]`. -Partitioning the rows and columns of `A` into the connected components of its -bipartite support graph, columns in different components share no supported -row, so for `W` diagonal the sum needed is exactly the one over the rows of -`j`'s own component: +For diagonal `W`, each support component has the scale s[j] = sqrt(Σ_{i ∈ rows(comp(j))} abs(w[i])*a[i]^2) * b[j] -(the unweighted form is this with `w[i] = 1`). A nonzero off-diagonal -`W[i,i']` can couple rows from two different components; components joined by a -chain of such couplings form a group. Within a group, writing +The unweighted form uses `w[i] = 1`. Off-diagonal entries of `W` may join +components. For each joined group, define `M[p,q] = Σ_{i ∈ rows(p), i' ∈ rows(q)} a[i]*abs(W[i,i'])*a[i']` for the block -sum over components `p` and `q`, `abs(G[j,k]) <= M[p,q]*b[j]*b[k]` for `j ∈ p` -and `k ∈ q`. When `abs.(W)` is not symmetric neither is `G`, and since -`s[j]*s[k]` is a single number bounding both `abs(G[j,k])` and `abs(G[k,j])`, the -block sums enter only through their symmetrization -`Ms[p,q] = max(M[p,q], M[q,p])`. What remains is to divide each bound between its -two components: any `σ` with `σ[p]*σ[q] >= Ms[p,q]` for every `p`, `q` yields +sum, and symmetrize it as `Ms[p,q] = max(M[p,q], M[q,p])`. A cover `σ` of `Ms` +yields s[j] = σ[p]*b[j], j ∈ p -This is a symmetric cover of the `k×k` matrix `Ms`, where `k` is the number of -components in the group. The implementation computes -[`symcover_min`](@ref)`(AbsLog{2}(), Ms)`. For a component that no `W` entry -couples to another, this reduces to the diagonal-`W` formula. Entries of `G` -across distinct groups vanish, and unsupported columns get `s[j] = 0`. - -Rescaling `a -> γ*a`, `b -> b/γ` within each support component leaves `s` -unchanged. Under this rescaling, `Ms[p,q] -> γ[p]*γ[q]*Ms[p,q]` and the minimal -cover changes as `σ[p] -> γ[p]*σ[p]`; therefore `σ[p]*b[j]` is invariant. This -makes `s` suitable as an absolute scale, such as in a Levenberg-Marquardt term -`λ*Diagonal(s.^2)`. - -The invariant cover does not exist when a connected part of `Ms`'s support graph -has an edge, no loop, and is bipartite. Scaling one color class by `t` and the -other by `1/t` leaves `Ms` unchanged but changes the individual `σ` values. -`gramcover` throws an `ArgumentError` in this case. With -`degenerate=:uniform`, it instead uses -`σ[p] = sqrt(Σ_{p,q} Ms[p,q])`, which depends on the gauge of `(a, b)`. A loop -or odd cycle removes this degeneracy. - -For uncoupled components, `s[j] <= norm(a)*b[j]`, with a strict inequality when -another component carries weight. There is no uniform comparison for coupled -components because changing the gauge redistributes tightness among them. +The implementation computes `σ` with +[`symcover_min`](@ref)`(AbsLog{2}(), Ms)`. Unsupported columns receive zero. -When a positive-semidefinite `W` is available only as an operator — `W[i,i]` -readable, `W[i,i']` for `i != i'` not — `abs(W[i,i']) <= sqrt(W[i,i]*W[i',i'])` -yields the looser diagonal-only bound -`s[j] = (Σ_{i ∈ rows(comp(j))} sqrt(W[i,i])*a[i]) * b[j]`, computable by hand -from `diag(W)`. The methods here always compute the tighter entrywise form -above, which requires `W`'s entries. +Under componentwise rescaling, `Ms[p,q]` and `σ[p]` transform so that +`σ[p]*b[j]` remains unchanged. -The minimal-cover solve has size `k`, the number of support components coupled -by `W`. Constructing `Ms` also costs `O(k^2)` and usually dominates the solve. - -Roundoff margins on the block sums and a final feasibility check preserve the -cover in floating-point arithmetic. +An invariant cover does not exist when a nontrivial connected component of +`Ms` is loopless and bipartite. With `degenerate=:uniform`, the fallback is +`σ[p] = sqrt(Σ_{p,q} Ms[p,q])`, which depends on the gauge of `(a, b)`. A loop +or odd cycle removes this degeneracy. See also: [`gramcover!`](@ref), [`symcover`](@ref), [`cover`](@ref), [`iscover`](@ref). @@ -161,9 +121,8 @@ end Mutating counterpart of [`gramcover`](@ref): writes the symmetric cover of the (weighted) Gram matrix into `s` and returns it, rather than allocating a new vector. `eachindex(s)` must match `axes(A, 2)` — `sc.colax` for the -[`SupportComponents`](@ref) forms — in addition to the axis requirements -[`gramcover`](@ref) places on `a`, `b`, and `w`/`W`, and shares its -`degenerate` keyword. +[`SupportComponents`](@ref) forms. The `degenerate` keyword is shared with +[`gramcover`](@ref). See also: [`gramcover`](@ref). """ @@ -222,11 +181,7 @@ function gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, sc: throw(DimensionMismatch("`W` couples support rows, so it must be square on the row axis: axes(W) must be $(string((sc.rowax, sc.rowax))), got $(string(axes(W)))")) ncomp = ncomponents(sc) - # Union-find over the component ids of `sc`: merge two of them whenever a - # nonzero `W[i,i']` couples a supported row of one to a supported row of the - # other, a coupling `A`'s own support graph does not carry but the product - # `A'*W*A` does. Rows with no support (component id 0) contribute nothing to - # `A'*W*A` regardless of `W`, so they are skipped. + # Merge support components coupled by nonzero entries of `W`. parent = collect(1:ncomp) function find(p) while parent[p] != p @@ -250,8 +205,7 @@ function gramcover!(s::AbstractVector, a::AbstractVector, b::AbstractVector, sc: end end - # With every component still its own group, each needs only its own block sum, - # and the group bookkeeping below would allocate per component to say so. + # Avoid group bookkeeping when all components remain independent. if !merged m = zeros(typeof(_gc_term(a, W)), ncomp) n = zeros(Int, ncomp) @@ -344,11 +298,9 @@ _gc_eltype(a, b, args...) = typeof(sqrt(_gc_term(a, args...)) * zero(eltype(b))) # Convert uncoupled component sums to scales. # -# Naive summation of `n` nonnegative terms computes `fl(Σ) >= Σ/(1+γ)` with -# `γ = n*ulp/(1-n*ulp)`, `ulp = eps(scalarT)/2`; inflating `sqrt(m[c])` by -# `1 + (n[c]+3)*eps(scalarT)` absorbs that shortfall together with the roundoff of -# the `sqrt` itself and of the multiply by `b[j]`, so the cover holds despite -# `A'*W*A` never being formed to check it. +# Inflate so the cover holds in floating point without forming `A'*W*A`: naive +# summation of `n` nonnegative terms can fall short of the true sum by a factor +# `1 - n*eps`, and the `sqrt` and the multiply by `b[j]` each add a rounding. function _write_gramcover!(s::AbstractVector, b::AbstractVector, colcomp::Vector{Int}, oc, m::AbstractVector, n::AbstractVector{Int}) scalarT = scalar_type(eltype(m)) sq = [sqrt(m[c]) * (1 + (n[c] + 3) * eps(scalarT)) for c in eachindex(m)] @@ -442,10 +394,7 @@ function _loopless_bipartite(Ms::AbstractMatrix) return false end -# Write `s[j] = sq[c]*b[j]` for `j` in component `c` (`colcomp[j - oc]`), and -# `s[j] = 0` for a column with no support. `sq` carries different units than the -# sums it came from (e.g. sums in `s` give `sq` in `s^(1/2)`), so its element type -# is inferred from the computation rather than taken from those sums. +# Write component scales to `s`, inferring the post-square-root element type. function _write_gramcover_sq!(s::AbstractVector{T}, b::AbstractVector, colcomp::Vector{Int}, oc, sq::AbstractVector) where T for j in eachindex(s) c = colcomp[j-oc] diff --git a/src/heuristic_covers.jl b/src/heuristic_covers.jl index b0c8a83..74b74ac 100644 --- a/src/heuristic_covers.jl +++ b/src/heuristic_covers.jl @@ -78,11 +78,8 @@ Given a matrix `A`, return vectors `a` and `b` such that and column geometric means, covers the most-violated entries first, then applies `maxiter` tightening iterations. -Only the products `a[i] * b[j]` are determined by the problem: `a -> c*a`, `b -> b/c` -leaves every one of them unchanged. The split is fixed by the balance convention -`∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j]` (`nzaᵢ`, `nzbⱼ` = nonzero counts of row `i`, -column `j`), imposed within each connected component of the support (the gauge acts -independently on each), as it is throughout the package; see [`cover_min`](@ref). +The factors use the per-component balance convention described by +[`cover_min`](@ref). `ϕ` is accepted for API compatibility but is currently ignored. For a cover that provably minimizes a given `ϕ`, use [`cover_min`](@ref). @@ -161,35 +158,9 @@ end # ============================================================ # Internal helpers # ============================================================ -# Shift `(a, b)` along the gauge `a -> c*a`, `b -> b/c`, which leaves every product -# `a[i]*b[j]` — and hence every objective and every coverage constraint — untouched, -# onto the balance convention `∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j]` that every asymmetric -# cover in the package reports its result in. Summing `log a[i]` over the support counts -# row `i` exactly `nzaᵢ` times, which is those weighted sums. Nothing else pins the gauge: -# the objective cannot see it, so without a convention the split between `a` and `b` would -# be an artifact of whichever pass last touched them. -# -# The gauge acts independently on each connected component of the bipartite -# support graph (`_support_components`), so the convention is imposed per -# component: within every component, the row-side and column-side weighted log -# sums agree to within the rounding described below. This makes the split a -# well-defined function of the support and the products — block-diagonal -# assembly commutes with balancing — rather than pinning only the global scalar -# and leaving the per-component splits to whichever solver internals ran last. -# Rows and columns with empty support belong to no component and are left -# untouched. -# -# The shift is rounded to a whole power of two before it is applied. Scaling `a` -# by `2^k` and `b` by `2^-k` is exact in binary floating point, so every product -# `a[i]*b[j]` is preserved bit for bit and no cover is perturbed into -# infeasibility by the act of pinning its gauge. The balance is therefore met to -# within a factor of `√2` rather than exactly — a bound on the residual -# imbalance, traded for exactness of the quantity that carries the meaning. -# -# Two points differing only by the gauge land on the same point here, so a refiner given -# either start cannot tell them apart. The uniform inflation to feasibility raises every -# supported scale of `a` and `b` alike, and within each component -# `∑ nzaᵢ = ∑ nzbⱼ = nnz`, so it preserves the balance it finds. +# Apply the row/column balance convention independently to each support +# component. Rounding the shift to a power of two preserves cover products +# exactly, at the cost of balancing only within a factor of `sqrt(2)`. function _balance_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix) T = float(promote_type(eltype(a), eltype(b))) rowcomp, colcomp, ncomp = _support_components(A) @@ -205,8 +176,7 @@ function _balance_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix Lβ[c] += log2(T(b[j])) nnz[c] += 1 end - # `Lα`, `Lβ` are log2 sums, so the shift is already an exponent: rounding it - # to an integer is what makes the rescaling below exact. + # An integer base-2 exponent makes the rescaling exact. gamma = [exp2(round((Lβ[c] - Lα[c]) / (2 * nnz[c]))) for c in 1:ncomp] for i in eachindex(a) c = rowcomp[i-or] @@ -222,10 +192,10 @@ function _balance_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix end -# Compute the analytical minimizer of the unconstrained AbsLog{2} symmetric objective +# Analytical minimizer of the unconstrained `AbsLog{2}` symmetric objective # ∑_{i,j: A[i,j]≠0} (log(a[i]*a[j]) - log|A[i,j]|)² -# Fills `a` in-place and returns nza[i] = number of nonzero entries in row i. -# For efficiency, uses a Sherman-Morrison approximation for the pattern of nonzeros. (It's exact when there are no zeros.) +# Returns row support counts. The Sherman-Morrison approximation is exact on +# complete support. function unconstrained_min!(::AbsLog{2}, a::AbstractVector{T}, A::AbstractMatrix) where T ax = eachindex(a) axes(A) == (ax, ax) || throw(DimensionMismatch("`unconstrained_min!(ϕ, a, A)` requires a square matrix with matching axes to `a` (got axes(A)=$(string(axes(A))), axes(a)=$(string(axes(a)))")) @@ -299,14 +269,8 @@ function init_feasible_diag!(a::AbstractVector{T}, A::AbstractMatrix) where T return boost_feasible_seq!(a, A) end -# Shrink x by exp(lr/2) (lr = log of the cover-to-entry ratio for the tightest -# entry touching x). The direct quotient degenerates to exact zero when -# exp(lr/2) overflows or the division underflows; recomputing in log space -# recovers any representable result, and the floatmin clamp handles genuine -# underflow: x is supported (nonzero going in), and an exact zero would make -# every entry through it permanently uncoverable, whereas floatmin keeps it -# representable and, being the smallest normal positive magnitude, changes the -# resulting cover products negligibly. +# Shrink `x` by `exp(lr/2)` in log space, clamping supported scales at +# `floatmin` on underflow. function _tighten_shrink(x, lr) T = float(promote_type(typeof(x), typeof(lr))) y = T(x) / exp(T(lr) / 2) @@ -393,26 +357,12 @@ function tighten_cover!(a::AbstractVector, b::AbstractVector, A::Transpose; kwar return a, b end -# Feasibility boost by approximate greedy max-deficit. `entries` holds the -# support entries; `deficit(e)` returns the log-deficit z = log|A[i,j]| minus -# the log of the current cover product, > 0 iff violated; `apply!(e, z)` grows -# the scales so the entry becomes exactly covered. Deficits only shrink as -# scales grow, so entries only move to lower buckets: total work is -# O(#entries + moves), moves per entry bounded by the bucket count. Bucket -# edges are anchored at 0 in log-deficit (a scale-invariant quantity), so -# processing order — hence the result — is covariant under diagonal rescaling -# of A, up to within-bucket ties. Working in log-deficit keeps every quantity -# finite for finite nonzero entries and positive scales, however extreme the -# dynamic range. +# Approximate greedy max-deficit boost. Deficits only decrease, so entries move +# to lower buckets. Log-deficit buckets preserve covariance except for ties. const BOOST_BUCKET_WIDTH = log(2) / 4 # quality indistinguishable from exact greedy; only bucket count grows as w shrinks -# Buckets are a flat singly-linked bucket queue (as in bucket-queue Dijkstra), -# not one growable Vector{Int} per bucket: `head[b]` is the top-of-stack index -# into `entries`, `nxt[k]` the next index below it in whatever bucket k -# currently occupies. Push/pop are O(1) pointer updates with no reallocation, -# and `deficit(entries[k])` is cheap enough (a couple of array reads and a -# subtraction) that recomputing it on each of the two passes below costs less -# than caching it in a separate array would. +# Flat linked bucket queue: `head[b]` is the first entry and `nxt[k]` links the +# rest. Deficits are recomputed instead of cached. function bucket_boost!(deficit::F, apply!::G, entries, ::Type{T}) where {F,G,T} n = length(entries) zmax = zero(T) @@ -546,21 +496,9 @@ function boost_feasible!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix return a, b end -# Feasibility by sequential nearest-neighbor propagation with deferral, -# processing off-diagonal pairs in order of increasing offset -# j = 1, …, n-1 (each nonzero A[k, k+j] requires a[k]*a[k+j] ≥ |A[k, k+j]|). -# `a[k] == 0` means "not yet resolved" rather than "unsupported": entries with -# both endpoints still zero are deferred until a later offset (or another -# deferred entry) supplies a scale for one of them, and any left unresolved -# after every offset is processed are equal-split as a[k]=a[l]=√|A[k,l]|. -# -# When both a[k] and a[l] are already nonzero but a[k]*a[l] < |A[k,l]|, both -# are scaled by the square root of the ratio √(|A[k,l]|/(a[k]*a[l])). Equal -# scaling is of course ad-hoc; while it might be better to do something tuned -# to a particular penalty function, that would risk making the algorithm -# O(n^3) (we'd likely need to revisit earlier offsets), and earlier decisions -# might be reversed by later ones anyway. For something intended as an -# initialization, a heuristic guaranteed to be O(n^2) seems reasonable. +# Sequential nearest-neighbor feasibility propagation in increasing diagonal +# offset. Zero scales are unresolved; deferred pairs are revisited, then split +# equally if neither endpoint acquires a scale. The method costs O(n²). function boost_feasible_seq!(a::AbstractVector{T}, A::AbstractMatrix) where T ax = eachindex(a) axes(A) == (ax, ax) || throw(DimensionMismatch("`boost_feasible_seq!(a, A)` requires a square matrix with matching axes to `a` (got axes(A)=$(string(axes(A))), axes(a)=$(string(axes(a))))")) @@ -638,18 +576,8 @@ function boost_feasible_seq!(a::AbstractVector{T}, A::AbstractMatrix) where T return a end -# Feasibility by the smallest uniform inflation: multiply every scale by the same -# factor until `a[i]*a[j] >= |A[i,j]|` for every entry visited by -# `foreach_support_sym`. Requires a start with strictly positive scale on every -# supported row (the geometric-mean init from `unconstrained_min!` guarantees this). -# -# Unlike `boost_feasible!`, this preserves the shape of the starting point. The -# two methods can reach different basins of the nonconvex AbsLinear objective. -# The shift depends -# on `A` only through the log-deficits at the starting point, which are invariant -# under a diagonal rescaling `D*A*D`, so the result is scale-covariant. Growing the -# log-scales directly (rather than multiplying by `exp(t)`) stays finite even when -# `exp(t)` alone would overflow. +# Apply the smallest uniform inflation that covers `A`. This preserves the +# starting point's shape and works in log space to avoid overflow. function inflate_feasible!(a::AbstractVector{T}, A::AbstractMatrix) where T la = map(log, a) tref = Ref(zero(T)) diff --git a/src/initializers.jl b/src/initializers.jl index 471cfa3..59cf9d8 100644 --- a/src/initializers.jl +++ b/src/initializers.jl @@ -1,9 +1,6 @@ # Starting covers shared by the soft-cover multistarts and `*_min` solvers. -# Starting covers the non-convex AbsLinear solvers refine, in the order they are tried. -# `:leaveout` and `:diagfeasible` have no asymmetric formulation, so the two menus differ. -# The hard-cover drivers take these starts as covers (`feasible=:inflate`), the soft-cover -# driver takes them raw (`feasible=:none`) — the soft objective constrains nothing. +# Default starts for nonconvex `AbsLinear` solvers. const SYMCOVER_MIN_STRATEGIES = (:hardcover, :geomean, :leaveout) const COVER_MIN_STRATEGIES = (:hardcover, :geomean) @@ -14,41 +11,26 @@ const COVER_MIN_STRATEGIES = (:hardcover, :geomean) """ a = initialize_symcover(A; strategy=:hardcover, feasible=:inflate, kwargs...) -Build a starting point for the symmetric cover of `A`, as consumed by -[`symcover_min`](@ref) and by the [`soft_symcover`](@ref) multistart. - -The strategies depend only on `A`, so this function takes no penalty. +Build a symmetric starting point for [`symcover_min`](@ref) or +[`soft_symcover`](@ref). Strategies depend only on `A`: `strategy` names the point: -- `:geomean` — the geometric mean of each row's nonzero entries. It is not - generally a cover. +- `:geomean` — geometric means of the nonzero entries in each row. - `:leaveout` — the geometric mean recomputed with the most-underweighted - support entry dropped, which lands in the basin that treats that entry as - effectively zero. Raises an `ArgumentError` when no entry can be dropped - (empty support, or dropping it would empty a row). Not a cover. + 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 tightened hard cover of [`symcover`](@ref), which is - `:geomean` boosted to feasibility and then tightened. Forwards `maxiter` to the - tightening pass. `feasible` has no effect on it. +- `:hardcover` — the result of [`symcover`](@ref). It forwards `maxiter` and + ignores `feasible`. -`feasible` names how the point is brought up to covering `A` — that is, to -`a[i]*a[j] >= abs(A[i,j])`, up to the roundoff of the log-domain arithmetic: +`feasible` controls whether and how the point is made into a cover: -- `:inflate` (the default) multiplies every scale by the smallest common factor - that achieves coverage. -- `:boost` raises only the rows that touch a violated entry, so it changes the - shape of the point. This is the route [`symcover`](@ref) itself takes. +- `:inflate` multiplies every scale by the smallest common factor that covers `A`. +- `:boost` raises scales that touch violated entries. - `:none` returns the strategy's point without a coverage guarantee. -The two feasible routes reach different points on the boundary and can enter -different basins of the nonconvex `AbsLinear` objectives. - -Under every setting the result is strictly positive on every row that carries -support and exactly zero on every row that carries none. - -An unrecognized `strategy` or `feasible` raises an `ArgumentError`. +Supported rows receive positive scales; unsupported rows receive zero. See also: [`initialize_symcover!`](@ref), [`initialize_cover`](@ref), [`symcover`](@ref), [`symcover_min`](@ref). """ @@ -83,22 +65,15 @@ end """ a, b = initialize_cover(A; strategy=:hardcover, feasible=:inflate, kwargs...) -Build a starting point for the cover of `A`, as consumed by [`cover_min`](@ref) -and by the [`soft_cover`](@ref) multistart. This is the asymmetric analog of -[`initialize_symcover`](@ref), and takes the same `feasible` keyword, under -which the result covers `A` as `a[i]*b[j] >= abs(A[i,j])`. +Build an asymmetric starting point for [`cover_min`](@ref) or +[`soft_cover`](@ref). The `feasible` keyword matches +[`initialize_symcover`](@ref). -Two of the strategies carry over: `:hardcover` (the tightened hard cover of -[`cover`](@ref), forwarding `maxiter`) and `:geomean` (the AbsLog{2} -unconstrained minimum). `:leaveout` and `:diagfeasible` have no asymmetric -formulation and raise an `ArgumentError`, as does any unrecognized `strategy` or -`feasible`. +Supported strategies are `:hardcover` (the result of [`cover`](@ref), forwarding +`maxiter`) and `:geomean` (the unconstrained `AbsLog{2}` minimum). -Under every `feasible` setting the result is strictly positive on every -supported row and column and exactly zero on the unsupported ones, and the split -between `a` and `b` is fixed by the balance convention -`∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j]`, imposed within each connected component of the -support, that every asymmetric cover in the package uses (see [`cover_min`](@ref)). +Supported rows and columns receive positive scales; unsupported ones receive +zero. The factors use the balance convention of [`cover_min`](@ref). See also: [`initialize_cover!`](@ref), [`initialize_symcover`](@ref), [`cover`](@ref), [`cover_min`](@ref). """ @@ -134,9 +109,7 @@ function initialize_cover!(a::AbstractVector, b::AbstractVector, A::AbstractMatr throw(ArgumentError("unknown strategy :$strategy; expected one of :hardcover, :geomean")) end _make_feasible!(feasible, a, b, A) - # `:boost` raises rows and columns independently and so moves the gauge; pin it, as every - # asymmetric cover in the package does. This is invisible to the refiners, which read a - # start only up to the gauge, but it means a start can be compared against a cover. + # Restore the package's balance convention after a selective boost. return _balance_cover!(a, b, A) end @@ -144,11 +117,7 @@ end # Internal helpers # ============================================================ -# Build the named symmetric start in `a` and return `true`, or return `false` — leaving `a` -# unspecified — when `A` admits no such start. Only `:leaveout` can decline, and only for -# want of a support entry it can drop. The two callers want opposite things there: -# `initialize_symcover!` raises the `ArgumentError`, since the caller named one strategy and -# did not get it, while a multistart forfeits the slot and refines the rest of its menu. +# Build a named start. `:leaveout` returns `false` when no entry can be removed. function _initialize_symcover!(a::AbstractVector, A::AbstractMatrix, strategy::Symbol, feasible::Symbol; kwargs...) if strategy === :hardcover @@ -169,14 +138,8 @@ function _initialize_symcover!(a::AbstractVector, A::AbstractMatrix, strategy::S return true end -# Raise a starting point onto the coverage boundary by the named route, or leave it -# where it is. The two routes land at different points — `inflate_feasible!` scales -# every entry by one common factor, `boost_feasible!` raises only the rows touching a -# violated entry — so which one is used is part of what names a start, not an -# implementation detail of reaching feasibility. -# Binding the argument count in `Vararg{Any,N}` lets Julia specialize the -# splatted calls below. Otherwise, juliac's trim verifier treats them as -# dynamic calls. +# Apply the selected feasibility step. The `Vararg` length keeps calls +# statically specialized for `juliac`. function _make_feasible!(feasible::Symbol, scales::Vararg{Any,N}) where N if feasible === :inflate inflate_feasible!(scales...) @@ -188,47 +151,23 @@ function _make_feasible!(feasible::Symbol, scales::Vararg{Any,N}) where N return nothing end -# Strategies with no tunables of their own must reject stray keywords rather than -# discard them: a forwarded `maxiter` that silently does nothing would misreport -# which start was built. +# Reject keywords unused by a strategy. function _reject_kwargs(strategy::Symbol, kwargs) isempty(kwargs) && return nothing throw(ArgumentError("strategy=:$strategy accepts no further keyword arguments, got $(string(join(keys(kwargs), ", ")))")) end -# Leave-one-out geometric mean. The geometric mean weights every nonzero entry equally, so -# an entry with |A[i,j]| far below the rest (in the scale-invariant sense of its log-residual -# z[i,j] = log|A[i,j]| - α[i] - α[j] at the unweighted minimum) skews the start into a worse -# basin than the exact-zero limit. Here the entry with the most negative residual is dropped -# from the support and the geometric mean recomputed, giving a start in the basin that treats -# that entry as effectively zero; the AbsLinear objective is finite at r = 0, so refinement -# then varies continuously as the entry vanishes. -# -# Scale-covariance: the residuals z are scale-invariant, so selecting the entry by argmin z -# is covariant, as is the reduced-support geometric mean. Residual ties are broken by -# ascending raw |A[i,j]| — NOT scale-invariant, but exact ties are precisely where -# covariance is unachievable: whenever A is scale-equivalent to a row/column permutation of -# itself (true of EVERY symmetric 2×2 with nonzero off-diagonal, via t² = A[2,2]/A[1,1]), -# the competing basins have exactly equal objectives, so no deterministic algorithm can be -# simultaneously scale-covariant, permutation-equivariant, and continuous there. The raw -# magnitude is the only continuity-relevant information left, and using it only on ties -# confines the covariance exception to that degenerate class. (Weighting all entries by raw -# |A[i,j]|² instead would carry per-entry physical units — incommensurate sums — and break -# covariance on an open set of matrices.) -# -# Returns `true` and fills `a` with the leave-one-out start, or returns `false` (leaving `a` -# unspecified) when no entry can be dropped: empty support, or dropping the selected entry -# would empty some row's support. +# Recompute the geometric mean after dropping the most negative log-residual. +# Residual ties use raw magnitude; this is the strategy's only covariance +# exception. Return `false` if no entry can be removed without emptying a row. function _leaveout_logmean_init!(a::AbstractVector{T}, A::AbstractMatrix) where T ax = eachindex(a) axes(A) == (ax, ax) || throw(DimensionMismatch("`_leaveout_logmean_init!(a, A)` requires a square matrix with matching axes to `a` (got axes(A)=$(string(axes(A))), axes(a)=$(string(axes(a))))")) nza = unconstrained_min!(AbsLog{2}(), a, A) sum(nza) == 0 && return false - # One gather serves all three passes below: the two pair scans read the `j >= i` - # half of it, the Gauss-Seidel sweeps read whole rows. + # Pair scans use the upper triangle; Gauss-Seidel uses complete rows. S = _sym_support(A, T) - # Most negative residual over the support, with a roundoff-tolerant tie set: exact ties - # (e.g. z[1,1] == z[2,2] for every 2×2) must not be ordered by floating-point noise. + # Use a roundoff-tolerant set for tied residuals. zmin = T(Inf) for i in ax, s in _slots(S, i) j = S.idx[s] @@ -247,16 +186,11 @@ function _leaveout_logmean_init!(a::AbstractVector{T}, A::AbstractMatrix) where ibest, jbest, Abest = i, j, Aij end end - # Dropping entry (i,j) removes one support count from row i and (if off-diagonal) row j. + # Account for both endpoints of an off-diagonal entry. nza[ibest] > 1 || return false ibest == jbest || nza[jbest] > 1 || return false - # Minimize the unconstrained AbsLog{2} objective over the reduced support by - # Gauss-Seidel on its normal equations, starting from the full-support solution - # already in `a`. The closed-form geometric-mean formula used by - # `unconstrained_min!` is exactly scale-covariant only for rank-1 support - # patterns, which the reduced support never is; a Gauss-Seidel update, by - # contrast, is exactly covariant from any covariant iterate, for any sweep - # count, so basin selection downstream cannot depend on the frame. + # Minimize the reduced-support `AbsLog{2}` objective by Gauss-Seidel. Starting + # from a covariant point preserves covariance at every sweep. α = similar(a) for i in ax α[i] = iszero(nza[i]) ? zero(T) : log(a[i]) diff --git a/src/iscover.jl b/src/iscover.jl index db871b0..52c48b7 100644 --- a/src/iscover.jl +++ b/src/iscover.jl @@ -1,5 +1,4 @@ -# The cover predicate. Kept beside the traversal it is built on: checking coverage is -# exactly a walk over the support, since an entry that is zero constrains nothing. +# Cover predicates. """ iscover(a, b, A; rtol=0, atol=0) @@ -13,12 +12,10 @@ requires `A` to be square. a[i]*b[j] >= abs(A[i,j])*(1 - rtol) - atol -The default for both tolerances is zero (no slack, test that the cover condition -holds); note that `atol != 0` breaks scale-invariance. +Both tolerances default to zero. A nonzero `atol` breaks scale invariance. `a` and `b` must be nonnegative; a negative scale raises an `ArgumentError`. -Zero is allowed, and is what every solver here returns for a row or column with -no support. +Zero is allowed for unsupported rows and columns. `eachindex(a)` must match `axes(A, 1)` and `eachindex(b)` must match `axes(A, 2)`. @@ -46,8 +43,7 @@ function iscover(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; rtol=0 throw(DimensionMismatch("indices of `b` must match column-indexing of `A`, got eachindex(b)=$(string(eachindex(b))), axes(A, 2)=$(string(axes(A, 2)))")) _require_nonneg(a, "a") _require_nonneg(b, "b") - # Zero entries of `A` are skipped by `foreach_support`, and need no check: they demand - # `a[i]*b[j] >= 0`, which nonnegative scales satisfy outright. + # Zero entries impose no constraint and are skipped by the traversal. covered = Ref(true) foreach_support(A) do i, j, v covered[] &= _iscovered(a[i] * b[j], v, rtol, atol) @@ -65,9 +61,8 @@ _iscovered(p, v, rtol, atol) = iszero(atol) ? p >= v * (1 - rtol) : p >= v * (1 function _require_nonneg(x::AbstractVector, name::String) for i in eachindex(x) - # `zero(x[i])`, not `zero(eltype(x))`: a dimensional scale carries its units in the - # value, and `zero(Quantity{Float64})` is undefined. Also rejects NaN, which fails - # every comparison. + # Ask the value for zero because dimensional abstract element types may + # not define `zero(eltype(x))`. This also rejects NaN. x[i] >= zero(x[i]) || throw(ArgumentError("iscover requires nonnegative scales, got $name[$(string(i))] = $(string(x[i]))")) end diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index 4da4c3c..c687b6f 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -1,6 +1,5 @@ -# Objective-minimal hard covers. The default AbsLog{2} penalty is solved natively -# here; the other penalties are provided by the MatrixCoversJuMPExt and MatrixCoversIpoptExt extensions, -# whose entry points are declared as stubs below. +# Objective-minimal hard covers. `AbsLog{2}` is native; extensions provide the +# other penalties. # ============================================================ # Public interface @@ -11,48 +10,32 @@ a = symcover_min(A; kwargs...) Return the ϕ-minimal symmetric hard cover of `A`: the vector `a` minimizing -`∑_{i,j} ϕ(|A[i,j]|/(a[i]*a[j]))` subject to `a[i]*a[j] >= |A[i,j]|` for every nonzero -entry of `A`. The no-ϕ form defaults to `AbsLog{2}()`, matching [`symcover`](@ref). +`∑ ϕ(|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}()`: solved natively (no external solver). Accepts keyword arguments - `κs` (the penalty-continuation schedule, default `(1e2, 1e4, 1e6, 1e8)`), - `maxiter` (Newton steps per stage, default `40`), and `linsolve` (the inner - linear solve). `:dense` factorizes the reweighted normal equations densely, - at O(n³) per Newton step. `:woodbury` solves the same equations as a sparse - correction of the complete-support ones: the matrix is a sparse symmetric - positive-definite matrix `C` plus `e*eᵀ`. Well-conditioned penalty stages - apply that sum without forming it and solve by Jacobi-preconditioned - conjugate gradients; the rest take a sparse Cholesky of `C` and a - Sherman–Morrison update. Both are exact to rounding. `:woodbury` requires - `Float64` arithmetic, a support missing at most `n ÷ 4` entries in any row, - and at most `4n` zero entries in total, and raises an `ArgumentError` - otherwise. `:lsqr` uses matrix-free LSQR (per-iteration cost O(nnz), - intended for large sparse supports), right preconditioned in `Float64` by - the diagonal of the unweighted normal matrix, joined by the rows the penalty - currently weights — through a sparse Cholesky — once diagonal scaling alone - would leave the system ill conditioned; this keeps its iteration count from - growing with the penalty strength. `:auto` - selects `:woodbury` where its requirements hold and `:dense` elsewhere. - `linsolve` defaults to `:auto` for dense `A`; the - `SparseMatrixCSC`/`Symmetric`/`Hermitian` sparse methods default to `:lsqr` - instead, since neither factorization of the reweighted normal equations is - the right solve when `nnz ≪ n²`. +- `AbsLog{2}()`: native. - `AbsLog{1}()`: requires JuMP and HiGHS. -- `AbsLinear{1}()`, `AbsLinear{2}()`: requires JuMP and Ipopt. These objectives are - nonconvex. Each strategy in `strategies` is refined, and the best local - minimum is returned. A strategy that cannot produce a start is skipped. +- `AbsLinear{1}()` and `AbsLinear{2}()`: require JuMP and Ipopt and return the + best local minimum found from `strategies`. -The `AbsLog` penalties are convex in the log-scales. `AbsLog{2}` has a unique -minimizer. When the `AbsLog{1}` optimum is a face, the method returns the member -that minimizes the `AbsLog{2}` objective. +`AbsLog` is convex in the log scales. If `AbsLog{1}` has multiple minima, the +method selects the one with the smallest `AbsLog{2}` objective. -The native `AbsLog{2}` solver runs its penalty continuation in `Float64` when `A` -works in a narrower type, whose resolution the continuation's tolerances outrun, and -returns the cover in the element type `A` calls for. +# Extended help -!!! note - Even the native solver is more expensive than the [`symcover`](@ref) heuristic. +The native solver accepts `κs` (penalty-continuation schedule), `maxiter` +(Newton steps per stage), 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 and `:dense` otherwise. + +The native solver computes in `Float64` for narrower input types, then converts +the result to the required element type. See also: [`cover_min`](@ref), [`symcover`](@ref), [`symcover_min!`](@ref). """ @@ -64,52 +47,27 @@ symcover_min(A::AbstractMatrix; kwargs...) = symcover_min(AbsLog{2}(), A; kwargs a, b = cover_min(A) Return the ϕ-minimal asymmetric hard cover of `A`: vectors `a`, `b` minimizing -`∑_{i,j} ϕ(|A[i,j]|/(a[i]*b[j]))` subject to `a[i]*b[j] >= |A[i,j]|` for every nonzero -entry of `A`. The split between `a` and `b` is set within each support component -by the balance convention +`∑ ϕ(|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 no-ϕ form defaults to `AbsLog{2}()`, matching -[`cover`](@ref). +column `j`). Supported ϕ values: -- `AbsLog{2}()`: solved natively (no external solver). Accepts keyword arguments - `κs` (the penalty-continuation schedule, default `(1e2, 1e4, 1e6, 1e8)`), - `maxiter` (Newton steps per stage, default `40`), and `linsolve` (the inner - linear solve). `:dense` factorizes the reweighted normal equations densely, - at O((m+n)³) per Newton step. `:woodbury` solves the same equations as a - sparse correction of the complete-support ones: the matrix is a sparse - symmetric positive-definite matrix `C` plus a rank-two term. Well-conditioned - penalty stages apply that sum without forming it and solve by - Jacobi-preconditioned conjugate gradients; the rest take a sparse Cholesky of - `C` and a Woodbury update. Both are exact to rounding. `:woodbury` requires - `Float64` arithmetic, a support missing at most `min(m, n) ÷ 4` entries in - any row or column, and at most `4·max(m, n)` zero entries in total, and - raises an `ArgumentError` otherwise. `:lsqr` uses matrix-free - LSQR (per-iteration cost O(nnz), intended for large sparse supports), right - preconditioned in `Float64` by the diagonal of the unweighted normal matrix, - joined by the rows the penalty currently weights — through a sparse Cholesky - — once diagonal scaling alone would leave the system ill conditioned; this - keeps its iteration count from growing with the penalty strength. - `:auto` selects `:woodbury` where its requirements hold and `:dense` - elsewhere. `linsolve` defaults to `:auto` for dense `A`; the - `SparseMatrixCSC` sparse method defaults to `:lsqr` instead, since neither - factorization of the reweighted normal equations is the right solve when - `nnz ≪ n²`. +- `AbsLog{2}()`: native. - `AbsLog{1}()`: requires JuMP and HiGHS. -- `AbsLinear{1}()`, `AbsLinear{2}()`: requires JuMP and Ipopt. These objectives are - nonconvex. Each strategy in `strategies` is refined, and the best local - minimum is returned. +- `AbsLinear{1}()` and `AbsLinear{2}()`: require JuMP and Ipopt and return the + best local minimum found from `strategies`. -The `AbsLog` penalties are convex in the log-scales. `AbsLog{2}` has a unique -minimizer. When the `AbsLog{1}` optimum is a face, the method returns the member -that minimizes the `AbsLog{2}` objective. +`AbsLog` is convex in the log scales. If `AbsLog{1}` has multiple minima, the +method selects the one with the smallest `AbsLog{2}` objective. -The native `AbsLog{2}` solver runs its penalty continuation in `Float64` when `A` -works in a narrower type, whose resolution the continuation's tolerances outrun, and -returns the cover in the element type `A` calls for. +# Extended help -!!! note - Even the native solver is more expensive than the [`cover`](@ref) heuristic. +The native solver accepts the same `κs`, `maxiter`, and `linsolve` keywords as +[`symcover_min`](@ref). 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`. See also: [`symcover_min`](@ref), [`cover`](@ref), [`cover_min!`](@ref). """ @@ -120,17 +78,12 @@ cover_min(A::AbstractMatrix; kwargs...) = cover_min(AbsLog{2}(), A; kwargs...) a = symcover_min!(ϕ, a, A; kwargs...) a = symcover_min!(a, A; kwargs...) -Refine the starting cover `a` into the ϕ-minimal symmetric hard cover of `A`, in -place. This is the second half of the initialize/refine pair: `a` must already be -a starting point, as produced by [`initialize_symcover`](@ref) (or by -[`symcover`](@ref)). The no-ϕ form defaults to `AbsLog{2}()`, matching -[`symcover_min`](@ref), whose keyword arguments and supported ϕ values these -methods share. +Refine the symmetric hard cover `a` in place. The no-ϕ form uses `AbsLog{2}()`; +supported penalties and keywords match [`symcover_min`](@ref). -`a` must be strictly positive on every row of `A` that carries support, and must -cover `A` — `a[i]*a[j] >= abs(A[i,j])` — to within the roundoff of the log-domain -arithmetic; otherwise an `ArgumentError` is raised. Scales on rows carrying no -support are inert: whatever they hold on input, they are zero on output. +`a` must cover `A` and be positive on supported rows. Unsupported scales are +ignored on input and set to zero. Use [`initialize_symcover`](@ref) or +[`symcover`](@ref) 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 @@ -146,17 +99,12 @@ symcover_min!(a::AbstractVector, A::AbstractMatrix; kwargs...) = a, b = cover_min!(ϕ, a, b, A; kwargs...) a, b = cover_min!(a, b, A; kwargs...) -Refine the starting cover `(a, b)` into the ϕ-minimal asymmetric hard cover of -`A`, in place. This is the asymmetric counterpart of [`symcover_min!`](@ref), and -carries the same contract on the start: strict positivity on every supported row -and column, coverage of `A` to within roundoff, and inert scales on the -unsupported rows and columns. The no-ϕ form defaults to `AbsLog{2}()`, matching -[`cover_min`](@ref), whose keyword arguments and supported ϕ values these methods -share. +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 and keywords match +[`cover_min`](@ref). -The product `a[i]*b[j]` is unchanged by `a -> c*a`, `b -> b/c`, so the start is -read only up to that gauge: `(a, b)` and `(2a, b/2)` give the same result, and -the result itself is pinned to the balance convention of [`cover_min`](@ref). +Equivalent starts `(c*a, b/c)` give the same balanced result. See also: [`initialize_cover`](@ref), [`cover_min`](@ref), [`symcover_min!`](@ref). """ @@ -164,40 +112,27 @@ function cover_min! end cover_min!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; kwargs...) = cover_min!(AbsLog{2}(), a, b, A; kwargs...) -# Symmetric AbsLog{2} hard cover via a one-sided quadratic penalty on the -# log-residuals z_ij = α_i + α_j - log|A_ij| (α = log a): +# Symmetric `AbsLog{2}` hard cover using a one-sided quadratic penalty on +# `z_ij = α_i + α_j - log|A_ij|`: # # f_κ(α) = ∑_{ij ∈ support} w(z_ij) z_ij², w = 1 for z ≥ 0, κ for z < 0. # -# As κ → ∞ the minimizer approaches the constrained (hard-cover) optimum. Each κ -# stage runs a damped semismooth Newton iteration: freeze the weights at the -# current α, solve the reweighted normal equations `B α = f` (an SDD system with -# the sparsity of the nonzero-pattern graph), and take a backtracking line -# search toward that point (which ensures convergence). A final uniform shift -# makes the cover exactly feasible. +# Each `κ` stage freezes the weights, solves the normal equations, and uses a +# backtracking line search. A final uniform shift restores feasibility. function symcover_min(::AbsLog{2}, A::AbstractMatrix; kwargs...) a, _ = _symcover_min_abslog2(A; kwargs...) return a end -# Asymmetric AbsLog{2} hard cover via the same one-sided quadratic penalty as -# `symcover_min`, on stacked log-scales x = (α; β) (α = log a over rows, β = log b -# over columns) with residuals z_ij = α_i + β_j - log|A_ij|. The row and column -# scales share a gauge freedom (α_i, β_j) → (α_i + s, β_j - s), one dimension per -# connected component of the support, that leaves every residual unchanged; during -# the solve the global one is fixed by adding v0*v0ᵀ, v0 = [ones(m); -ones(n)], to -# the normal equations, and afterwards the result is shifted, per component, to the -# balance convention ∑ nzaᵢ αᵢ = ∑ nzbⱼ βⱼ (nzaᵢ, nzbⱼ = nonzero counts of row i, -# column j, summed within the component) so it is deterministic — see -# `_cover_min_abslog2` for how the remaining per-component gauges are lifted and shifted. +# Asymmetric counterpart on stacked log scales `(α; β)`. The solve pins the +# global gauge; `_cover_min_abslog2` regularizes and balances the remaining +# component gauges. function cover_min(::AbsLog{2}, A::AbstractMatrix; kwargs...) a, b, _ = _cover_min_abslog2(A; kwargs...) return a, b end -# The AbsLog{2} objective is convex in the log-scales, so the continuation converges -# to the same cover from any start; the start is honored (it replaces the cold -# unweighted solve as the first iterate) but is not observable in the result. +# The convex objective has the same result from every valid start. function symcover_min!(::AbsLog{2}, a::AbstractVector, A::AbstractMatrix; kwargs...) _prepare_symcover_start!(a, A) anew, _ = _symcover_min_abslog2(A; start=a, kwargs...) @@ -213,10 +148,7 @@ function cover_min!(::AbsLog{2}, a::AbstractVector, b::AbstractVector, A::Abstra return a, b end -# Multistart drivers for nonconvex AbsLinear objectives. The `*_min!` kernels live -# in MatrixCoversIpoptExt; the main package owns the starts and selection. -# -# Use the same roundoff-tolerant selection rule as the soft-cover multistarts. +# Multistart drivers for the `AbsLinear` kernels in MatrixCoversIpoptExt. function symcover_min(ϕ::AbsLinear, A::AbstractMatrix; strategies=SYMCOVER_MIN_STRATEGIES) ax = axes(A, 1) axes(A, 2) == ax || throw(ArgumentError("symcover_min requires a square matrix")) @@ -224,8 +156,7 @@ function symcover_min(ϕ::AbsLinear, A::AbstractMatrix; strategies=SYMCOVER_MIN_ throw(ArgumentError("symcover_min: `strategies` must name at least one starting cover")) T = float(real(eltype(A))) starts = [similar(Array{T}, ax) for _ in strategies] - # A strategy for which `A` admits no start forfeits its slot; only a menu that yields - # no start at all leaves nothing to refine. + # Skip strategies that cannot produce a start for `A`. built = [_initialize_symcover!(a, A, strategy, :inflate) for (a, strategy) in zip(starts, strategies)] covers = [symcover_min!(ϕ, a, A) for (a, ok) in zip(starts, built) if ok] isempty(covers) && @@ -247,20 +178,13 @@ end # Internal helpers # ============================================================ -# Log-domain slack allowed of a start supplied to the `*_min!` refiners, in units of -# `eps(T)` scaled by the magnitudes entering the residual. The heuristics reach the -# coverage boundary through log-domain updates and so land on it only to within their -# own roundoff — a fresh `symcover` violates `a[i]*a[j] >= abs(A[i,j])` by a fraction -# of one such unit — and an exact test would reject them. This bound accepts that -# while still rejecting a start that misses coverage by any margin a solver would see. +# Roundoff allowance when validating log-domain heuristic starts. const START_FEASIBILITY_ULPS = 64 _start_slack(lv::T, li::T, lj::T) where {T} = START_FEASIBILITY_ULPS * eps(T) * max(oneunit(T), abs(lv), abs(li), abs(lj)) -# Shared prologue of the `symcover_min!` kernels: check that the caller's start is a -# cover of `A`, discard the inert scales on unsupported rows, and move the start onto -# the coverage boundary exactly, so every kernel begins from a feasible point. +# Validate and normalize a symmetric hard-cover start. function _prepare_symcover_start!(a::AbstractVector, A::AbstractMatrix, fname=:symcover_min!) ax = axes(A, 1) axes(A, 2) == ax || throw(ArgumentError("$fname requires a square matrix")) @@ -287,10 +211,7 @@ function _prepare_symcover_start!(a::AbstractVector, A::AbstractMatrix, fname=:s return inflate_feasible!(a, A) end -# Shared prologue of the `cover_min!` kernels; the asymmetric counterpart of -# `_prepare_symcover_start!`. The start is additionally pinned to the balance -# convention (imposed within each connected component of the support), so the -# refiners read it only up to the per-component row/column gauge. +# Validate, normalize, and balance an asymmetric hard-cover start. function _prepare_cover_start!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix) axes(A, 1) == eachindex(a) || throw(DimensionMismatch("indices of `a` must match row-indexing of `A`, got eachindex(a)=$(string(eachindex(a))), axes(A, 1)=$(string(axes(A, 1)))")) axes(A, 2) == eachindex(b) || throw(DimensionMismatch("indices of `b` must match column-indexing of `A`, got eachindex(b)=$(string(eachindex(b))), axes(A, 2)=$(string(axes(A, 2)))")) @@ -321,72 +242,27 @@ function _prepare_cover_start!(a::AbstractVector, b::AbstractVector, A::Abstract end -# Inner linear solve for the AbsLog{2} MMC Newton steps. -# -# `:dense` forms and factorizes the reweighted normal equations densely, at O(n³) -# per step. -# -# `:woodbury` splits the same matrix as `C + U Uᵀ`, where `C` is sparse (its -# off-diagonal pattern is the zero set `Z` of `A` together with the currently -# violated entries `V`) and symmetric positive definite, and `U` has one column -# (symmetric) or two (asymmetric). `C` is assembled sparsely on every such solve, and -# two sub-paths then take it, both exact to rounding, which is what the -# sign-stability stopping test in the continuation loop requires. A sparse Cholesky -# of `C` plus a Woodbury update — Sherman–Morrison, in the one-column symmetric case -# — costs far less than the dense factorization whenever `A` is close to fully -# supported. Alternatively `C + U Uᵀ` is applied as `C·x` plus the low-rank term, at -# O(nnz(C)) per application; Gershgorin on `(κ−1)·L_V` against the complete-support -# diagonal gives `1 + (κ−1)·2·maxdeg(V)/n` as an estimate of its condition number -# (the sharp bound is a small multiple of that), and while the estimate stays under -# `WOODBURY_CG_KAPPA`, Jacobi-preconditioned conjugate gradients converge to rounding -# in a few hundred such applications — cheaper than a factorization whose fill, on the -# near-random violated pattern of the early stages, approaches dense. Above it the -# factorization runs, as it does for any CG run that exhausts its iteration cap. +# Inner solves for `AbsLog{2}` Newton steps: # -# `C` is positive definite because the complete-support matrix contributes `n` (or -# `m`) to each diagonal while the zero set subtracts a signless Laplacian `L_Z` with -# λmax(L_Z) ≤ 2·maxdeg(Z); requiring at most a quarter of a row to be zero keeps the -# difference bounded below by half the diagonal. A second requirement is about cost -# rather than definiteness: `Z` enters every matvec and every factorization, so the -# path is taken only while the total number of zeros is O(n). CHOLMOD is the sparse -# factorization behind it, and it is reliable only in `Float64`, so that is the only -# working type the path accepts. -# -# `:auto` takes `:woodbury` where it applies and `:dense` otherwise. -# -# `:lsqr` forces the matrix-free path, whose per-iteration cost is O(nnz); it is the -# intended solve for large sparse supports (where nnz ≪ n²) and is used by the -# structured/sparse methods. In `Float64` it is right preconditioned, which is what -# keeps its iteration count from growing as the continuation raises κ. The -# preconditioner is `M = diag(RᵀR) + (κ−1)·Σ_{e∈V} rₑ·rₑᵀ`: the diagonal of the -# unweighted normal matrix, together with the exact contribution of the rows LSQR -# weights by κ. All of the κ-dependence of `RᵀWR` sits in those rows, and every -# generalized eigenvalue of `(RᵀWR, M)` is a mediant of eigenvalues of -# `(RᵀR, diag(RᵀR))` and so lies in their range. `M = K·Kᵀ` and LSQR runs on -# `√W·R·K⁻ᵀ` in the variable `y = Kᵀ·x`. While the same condition-number estimate, -# taken against the unweighted diagonal, stays under `LSQR_PRECOND_KAPPA` the -# violated rows are left out and `K` is the diagonal `sqrt.(diag(RᵀR))`, applied -# without forming anything; above it they are included and `K` is the permuted sparse -# Cholesky factor of `M`, applied through the CHOLMOD factor components `F.PtL` and -# `F.UP`. `Kᵀ` is never needed as a product: the warm start `Kᵀ·x₀` is `K⁻¹·(M·x₀)`, -# which those same components and one sparse matrix-vector product supply. - -# Condition-number estimate above which a Woodbury solve is factorized rather than -# iterated: past it conjugate gradients need more applications than the sparse -# Cholesky costs. +# - `:dense` factorizes the normal equations. +# - `:woodbury` represents them as sparse `C + U*U'`, using conjugate gradients +# while the condition estimate is small and sparse Cholesky otherwise. It is +# restricted to nearly dense `Float64` problems. +# - `:lsqr` applies the weighted residual operator `M` matrix-free. In `Float64`, its +# right preconditioner includes violated rows once diagonal scaling is inadequate. +# It is not interchangeable with CG on the normal equations: LSQR's accuracy +# tracks the condition number of `M` (≈ √κ), CG's that of `MᵀM` (≈ κ), and at +# κ = 1e8 the latter exhausts double precision. +# - `:auto` selects `:woodbury` when supported and `:dense` otherwise. + +# Condition estimate above which Woodbury uses sparse Cholesky instead of CG. const WOODBURY_CG_KAPPA = 1000 -# Condition-number estimate above which the LSQR preconditioner takes in the rows the -# penalty currently weights; below it diagonal scaling alone leaves the system well -# enough conditioned, and no factorization is formed. +# Condition estimate above which LSQR includes the weighted rows in its preconditioner. const LSQR_PRECOND_KAPPA = 1000 -# `(C + U·Uᵀ) x = f` solved from a factorization `F` of the sparse `C`, by the -# Woodbury identity `x = y − Y·((I + Uᵀ·Y) \ (Uᵀ·y))` with `y = C\f` and `Y = C\U`. -# One multi-right-hand-side solve of `[f U]` supplies both, and the capacitance is -# `k×k` for `U` of `k` columns: `k = 1` for the symmetric gauge `e`, where this is -# Sherman–Morrison, and `k = 2` for the asymmetric row and column indicators. `rhs` -# is the `size(U, 1)×(k+1)` buffer the block right-hand side is staged in. +# Solve `(C + U*U')x = f` from a sparse factorization of `C` using the +# Woodbury identity. `rhs` stores the combined `[f U]` solve. function _woodbury_solve!(x, F, U, f, rhs) k = size(U, 2) copyto!(view(rhs, :, 1), f) @@ -403,20 +279,8 @@ function _woodbury_solve!(x, F, U, f, rhs) return mul!(x, Y, g, -1, 1) end -# Jacobi-preconditioned conjugate gradients for the symmetric positive-definite -# Woodbury system `B x = f`, with `Bmul!(y, x)` applying `B` and `dg` holding its -# diagonal. `x` carries the warm start in and the iterate out; `r`, `z`, `d`, `Ad` -# are work vectors of the same length. Returns `(iters, converged)`. -# -# The Newton step has to be exact to rounding for the stage's sign-stability -# stopping test to mean what it says, so `tol` sits at the level of `eps` and a run -# that exhausts `maxiter` reports failure instead of a partial answer; the caller -# then falls back to the factorization, which is exact. -# -# `r` is carried by a recurrence that drifts from `f − B x`, so success is never -# declared on it: a claim of convergence is confirmed against a freshly computed -# residual, and a disagreement restarts the iteration there. The confirming -# application counts against `maxiter` like any other. +# Jacobi-preconditioned CG for `B*x = f`. Convergence is confirmed with a fresh +# residual; failure lets the caller fall back to factorization. function _pcg!(Bmul!, x, dg, f, r, z, d, Ad, maxiter::Int, tol) iters = 0 Bmul!(r, x) @@ -456,23 +320,9 @@ function _pcg!(Bmul!, x, dg, f, r, z, d, Ad, maxiter::Int, tol) return iters, fresh && nrm <= tol end -# Matrix-free LSQR (Paige & Saunders) for the weighted least-squares problem -# `min ‖M x - b‖` underlying the reweighted normal equations `MᵀM x = Mᵀb`. -# `Amul!(y, x)` overwrites `y` with `M*x`; `Atmul!(z, y)` overwrites `z` with -# `Mᵀ*y`. Warm-started from `x0`. LSQR is used in preference to CG on the normal -# equations because it works with the condition number of `M` (≈ √κ at penalty -# strength κ) rather than that of `MᵀM` (≈ κ); at κ = 1e8 the squared conditioning -# breaks CG while LSQR stays accurate. -# -# A `Float64` caller passes a right-preconditioned operator, so `x` is then the -# preconditioned variable and `M` is `√W·R·K⁻ᵀ`; see the description of `:lsqr` in the -# inner-solve overview above. Other working types pass `√W·R` itself. -# -# The penalty least-squares problem is inconsistent (its optimal residual is -# nonzero), so the stopping test is on the normal-equations residual -# ‖Mᵀ(b - Mx)‖ ≤ atol · ‖M‖ · ‖b - Mx‖, both estimated from the bidiagonalization -# scalars (‖Mᵀr‖ = ϕbar·α·|c|, ‖r‖ = ϕbar, ‖M‖ from the Frobenius norm of the -# bidiagonal). Returns `(x, iters)`. +# Matrix-free LSQR for `min norm(M*x-b)`, warm-started from `x0`. The stopping +# test uses the normal-equations residual estimated from the bidiagonalization. +# Returns `(x, iters)`. function _lsqr(Amul!, Atmul!, b::AbstractVector{T}, x0::AbstractVector{T}; atol=5000 * eps(T), maxiter::Int=2 * (length(b) + length(x0)) + 100) where {T} x = copy(x0) @@ -525,38 +375,21 @@ function _lsqr(Amul!, Atmul!, b::AbstractVector{T}, x0::AbstractVector{T}; return x, iters end -# The two layouts the AbsLog{2} continuation reads its support through. Both present -# the same abstract object: a set of support entries, each pairing two unknowns `(p, q)` -# with a value `c = log|A_ij|` and contributing a residual `z = x[p] + x[q] - c`. -# -# Every support entry is represented once: in the symmetric case an off-diagonal entry -# and its mirror are the single pair `(p, q)` with `p < q`, the diagonal is `(p, p)`, -# and the `symmetric` flag of `SupportSystem` records that an off-diagonal entry stands -# for two residuals. Its multiplicity is therefore `mult = (symmetric && p != q) ? 2 : 1`, -# and every weighted quantity below carries it. - -# One stored element per support entry. Suits a support that is sparse relative to the -# grid, which is what the LSQR and dense paths face. +# Support layouts for residuals `x[p] + x[q] - log|A[i,j]|`. Symmetric +# off-diagonal entries are stored once with multiplicity two. + +# One element per support entry for LSQR and dense solves. struct EdgeList{T} edges::Vector{Tuple{Int,Int}} # support entries as pairs of unknowns cvals::Vector{T} # log|A_ij| per stored entry end -# The values on a dense grid, with `-Inf` (the image of a zero entry under `log`) -# marking a position outside the support. The sweeps then run as contiguous column -# passes that vectorize, at the cost of visiting the whole grid; the Woodbury path, -# which is taken only when the zero set is thin, is where that trade pays. -# -# Symmetric: `C` is `n×n` and only the upper triangle (`i ≤ j`) is read, entry `(i, j)` -# coupling unknowns `i` and `j`. Asymmetric: `C` is `m×n`, entry `(i, j)` coupling -# unknowns `i` and `m + j`. +# Dense grid for Woodbury sweeps; `-Inf` marks entries outside the support. struct Grid{T} C::Matrix{T} end -# The support of `A` as the linear system the AbsLog{2} continuation solves. The -# unknowns are stacked log-scales `x[1:N]` — the row scales alone for a symmetric -# problem, the row scales followed by the column scales for an asymmetric one. +# Linear system over stacked log scales. struct SupportSystem{T,S} N::Int supp::S # support layout: `EdgeList` or `Grid` @@ -572,8 +405,7 @@ struct SupportSystem{T,S} end SupportSystem{T}(N, supp::S, args...) where {T,S} = SupportSystem{T,S}(N, supp, args...) -# The support sweeps, one method per layout. `symmetric` carries the multiplicity -# convention described above; on a `Grid` it also selects the layout's geometry. +# Support sweeps, specialized by layout. # Objective `f_κ(x) = Σ_e mult_e·w_e·z_e²`, `w_e = κ` where `z_e < 0` and 1 elsewhere. function _fκ(x, κ, supp::EdgeList{T}, symmetric::Bool) where {T} @@ -626,9 +458,7 @@ function _fκ(x, κ, supp::Grid{T}, symmetric::Bool) where {T} return v end -# The objective and the violated set at `x` from one sweep: the line search needs the -# value and the stage's stopping test needs to know whether the set still matches -# `pat`, and both read the same residuals. +# Compute the objective and violated set in one sweep. function _fκpat(x, κ, pat, supp::EdgeList{T}, symmetric::Bool) where {T} edges, cvals = supp.edges, supp.cvals v = zero(T) @@ -694,13 +524,11 @@ function _fκpat(x, κ, pat, supp::Grid{T}, symmetric::Bool) where {T} return v, ndiff == 0 end -# Storage for the violated set of one solve, in the layout's own shape. `Matrix{Bool}` -# rather than `BitMatrix` so that the column views the `Grid` sweeps take vectorize. +# Violated-set storage. `Matrix{Bool}` permits vectorized column views. _violation_pattern(supp::EdgeList) = falses(length(supp.edges)) _violation_pattern(supp::Grid) = fill(false, size(supp.C)) -# Shift to exact feasibility: the smallest γ ≥ 0 with `x[p] + x[q] + 2γ ≥ c` on the -# whole support. +# Smallest uniform log-scale shift that restores feasibility. function _boost_shift(x, supp::EdgeList{T}, symmetric::Bool) where {T} edges, cvals = supp.edges, supp.cvals γ = zero(T) @@ -740,14 +568,10 @@ function _boost_shift(x, supp::Grid{T}, symmetric::Bool) where {T} return γ end -# One Woodbury assembly sweep: the right-hand side `f`, the violated set `vedges` with -# its per-unknown count `degV`, and the diagonal `dg` of `B` on top of its -# already-initialized zero-set correction. `dκ = κ - 1`, or 0 for the cold unweighted -# solve, which `κ === nothing` marks and in which no entry counts as violated. -# -# The vectorized pass over a column writes `vpat` and accumulates `f`; a scalar scan of -# the same column then collects the violated entries. `w·c` is formed with `ifelse` on -# `isfinite(c)`: `0 * -Inf` is NaN. +# Assemble the Woodbury right-hand side, violated set, degrees, and diagonal. +# `κ === nothing` denotes the unweighted solve. Off-support entries of `C` are +# `-Inf`, so products with them go through `ifelse(isfinite(c), ...)`: `0 * -Inf` +# is NaN. function _assemble_woodbury!(f, dg, degV, vedges, vpat, x, κ, supp::Grid{T}, symmetric::Bool, dκ) where {T} C = supp.C @@ -827,22 +651,9 @@ function _assemble_woodbury!(f, dg, degV, vedges, vpat, x, κ, supp::Grid{T}, return f end -# The AbsLog{2} penalty continuation on `sys`: a sequence of stages of increasing κ, -# each a sequence of reweighted Newton steps with a backtracking line search, starting -# from `x0` (or from the cold unweighted solve when `x0 === nothing`). Returns the -# stacked log-scales and the `stats` NamedTuple the callers pass on. With `boost`, the -# result is shifted uniformly to exact feasibility `x[p] + x[q] ≥ c` on the support. -# -# The objective counts each support entry with its multiplicity, `f_κ(x) = -# Σ_e mult_e·w_e·(x[p] + x[q] - c_e)²`, so a symmetric problem is weighted on the full -# grid rather than on one triangle. The normal equations assembled below are that -# system scaled by ½ — uniformly, so they have the same solution: a support entry puts -# `w` on `B[p,p]` and `B[p,q]` and `w·c` on `f[p]`, and the same again transposed when -# `q != p`. A symmetric diagonal entry, whose row of `R` is `2·e_p`, thereby collects -# `2w` on `B[p,p]` and `w·c` on `f[p]`. -# -# The layout of `sys.supp` selects the linear-solve path: a `Grid` is built exactly -# when the Woodbury path applies, and the LSQR and dense paths run on an `EdgeList`. +# `AbsLog{2}` penalty continuation. Each stage freezes residual weights, solves +# the weighted least-squares problem, and backtracks. `boost=true` applies a final +# feasibility shift. The support layout selects the inner solver. function _abslog2_continuation(sys::SupportSystem{T}, x0; κs, maxiter::Int, linsolve::Symbol, boost::Bool) where {T} N = sys.N @@ -852,47 +663,29 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; use_woodbury = supp isa Grid ne = supp isa EdgeList ? length(supp.edges) : 0 use_lsqr = linsolve === :lsqr - # CHOLMOD, which factors the LSQR preconditioner, is reliable only in Float64; - # other working types run the plain matrix-free iteration. + # CHOLMOD preconditioning is limited to Float64. use_precond = use_lsqr && T === Float64 - # Number of residuals a stored entry stands for. + # Residual multiplicity of each stored entry. symmetric = sys.symmetric mult = (p, q) -> (symmetric && p != q) ? 2 : 1 - # Diagonal of `C` before any entry is violated: the complete-support value, less - # one per zero entry at each of its ends (twice over for a symmetric zero on the - # diagonal, whose signless Laplacian row counts it at both). + # Base diagonal of `C`, corrected for the zero set. czero = copy(sys.dfull) for (p, q) in sys.zedges czero[p] -= oneunit(T) czero[q] -= oneunit(T) end - # Each Newton step freezes the weights at the current `x` and solves the reweighted - # least-squares problem `min ‖√W (R x - c)‖`, `(R x)_e = x[p] + x[q]`, whose normal - # equations are the signless Laplacian system `B x = f`. The dense path forms - # `B + v0·v0ᵀ` and factorizes it (a support-free variable gets an identity row; a - # minimal scale-relative ridge lifts what the gauge term leaves singular — the - # bipartite null space of a symmetric support such as `[0 1; 1 0]`, and the extra - # gauge each connected component beyond the first carries in the asymmetric case). - # The Woodbury path solves the same regularized system exactly, splitting - # `B + v0·v0ᵀ` as `C + U·Uᵀ` around the complete-support matrix and correcting `C` - # for the zero set and the violated entries. The LSQR path applies `√W R` and its - # transpose matrix-free, with the gauge as an appended row, and warm-starts from - # the incoming iterate; it solves the least-squares form directly, so its accuracy - # tracks the conditioning of `√W R` (≈ √κ) rather than that of `B` (≈ κ). + # Each Newton step solves `min norm(sqrt(W)*(R*x-c))`. Dense and Woodbury + # paths solve regularized normal equations; LSQR applies `sqrt(W)*R` + # matrix-free with a gauge row. f = zeros(T, N) ws = zeros(T, ne) # √weight per stored entry, frozen during one solve cv = zeros(T, ne + 1) # √weight · log|A_ij|, with a trailing 0 gauge target - # Entries the frozen weights of the current solve treat as violated. A full Newton - # step that leaves this pattern intact has landed on the stage's minimizer. + # Violated entries under the current frozen weights. vpat = _violation_pattern(supp) vedges = Tuple{Int,Int}[] # the violated entries of the current solve degV = zeros(Int, use_woodbury ? N : 0) # violated entries per unknown dg = zeros(T, use_woodbury ? N : 0) # diagonal of `B`, for the ridge and the CG preconditioner - # Diagonal of the unweighted normal matrix of the gauge-augmented system, - # `RᵀR + v0·v0ᵀ`, the base of the LSQR preconditioner: each stored entry puts its - # multiplicity at each of its ends, and an entry whose row of `R` is `2·e_p` puts 4. - # A variable with neither support nor gauge is given 1 so the preconditioner stays - # positive definite. + # Diagonal of the unweighted, gauge-augmented normal matrix. dpart = zeros(T, use_lsqr ? N : 0) if supp isa EdgeList && use_lsqr for (p, q) in supp.edges @@ -937,20 +730,15 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; solve_weighted = function (x, κ) nsolves[] += 1 if supp isa Grid - # `B + v0·v0ᵀ = C + U·Uᵀ` with `C = D − L_Z + (κ−1)·L_V`: the - # complete-support diagonal `D`, corrected by the zero set `Z` and by the - # currently violated entries `V`. One sweep over the grid collects the - # right-hand side, the violated set, and the diagonal of `B`; everything - # after it is O(|Z| + |V|). + # `B + v0*v0' = C + U*U'`, with `C` corrected for zero and + # violated entries. dκ = κ === nothing ? zero(T) : T(κ) - oneunit(T) fill!(f, zero(T)) copyto!(dg, czero) fill!(degV, 0) empty!(vedges) _assemble_woodbury!(f, dg, degV, vedges, vpat, x, κ, supp, symmetric, dκ) - # Same ridge as the dense path, so both solve the same regularized system: - # `U·Uᵀ` puts 1 on every diagonal of `B + v0·v0ᵀ`, and every variable has - # support here, so no identity row arises. + # Match the ridge used by the dense path. dmax = zero(T) maxdegV = 0 for p in 1:N @@ -961,12 +749,7 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; for p in 1:N dg[p] += oneunit(T) + ridge end - # These triplets store `C` in full rather than in one triangle: the same - # matrix then serves the matvec below and the factorization after it. - # `sparse` sums the duplicates, and the ridge rides on the diagonal. The - # zero set's diagonal contribution is already in `czero`, so only its - # off-diagonal entries are pushed here; the violated entries are assembled - # afresh on every solve, both diagonal and off-diagonal. + # Store full `C` for both matrix-vector products and factorization. empty!(Ci) empty!(Cj) empty!(Cv) @@ -1001,19 +784,13 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; end end C = sparse(Ci, Cj, Cv, N, N) - # Gershgorin on `(κ−1)·L_V` against the smallest complete-support diagonal - # estimates the condition number of `B`. While that estimate is small, - # conjugate gradients on `C·x + U·(Uᵀx)` reach the same answer in a few - # hundred O(nnz(C)) applications, which is far cheaper than a factorization - # whose fill on the near-random violated pattern of the early stages - # approaches dense. + # Use CG while the Gershgorin condition estimate remains small. κest = oneunit(T) + dκ * 2 * maxdegV / dmin if κest <= WOODBURY_CG_KAPPA copyto!(cgx, x) Bmul! = function (yy, xx) mul!(yy, C, xx) - # The columns of `U` are indicator vectors, so the low-rank term is - # a block sum broadcast back over the same block. + # Indicator columns make the low-rank term a block sum. for k in axes(U, 2) s = zero(T) for p in 1:N @@ -1060,17 +837,13 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; end g = ne + 1 # index of the appended gauge row if use_precond - # Diagonal scaling alone leaves a conditioning that grows with κ once - # the violated rows dominate a variable's diagonal; past that point - # they enter the preconditioner in full, and its Cholesky pays for - # itself in the iterations it removes. + # Include violated rows once diagonal scaling is inadequate. κest = oneunit(T) for p in 1:N κest = max(κest, oneunit(T) + dκ * 2 * mdiag[p] / dpart[p]) end if κest <= LSQR_PRECOND_KAPPA - # `K` is diagonal here, so it is applied by a scaling and nothing - # is assembled or factorized. + # Diagonal `K` needs only elementwise scaling. Dmul! = function (y, yv) @. px = yv / psqrt for (e, (p, q)) in enumerate(edges) @@ -1127,9 +900,7 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; MF = cholesky(Symmetric(Msp)) Kc = MF.PtL Uc = MF.UP - # CHOLMOD exposes no in-place solve for a factor component, so each - # application returns a fresh vector; the transpose product copies it - # into the buffer LSQR hands over, which is the only copy avoidable here. + # CHOLMOD factor-component solves allocate their result. Pmul! = function (y, yv) xv = Uc \ yv for (e, (p, q)) in enumerate(edges) @@ -1193,14 +964,8 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; B[q, p] += w end end - # A minimal scale-relative ridge on the supported diagonals, sized by the - # largest of them, lifts the gauge directions `v0·v0ᵀ` does not pin: the - # bipartite null space of a symmetric support, and the independent gauge - # each connected component of an asymmetric support beyond the first - # carries. The right-hand side is orthogonal to every gauge null vector, so - # the ridge leaves the recovered scales essentially unperturbed, and the - # gauge it fixes is unobservable — no product a_i·b_j spans two components. - # Support-free variables get an identity row. + # A small scale-relative ridge lifts unpinned gauge directions. + # Support-free variables receive an identity row. dmax = zero(T) for p in 1:N dmax = max(dmax, B[p, p]) @@ -1227,18 +992,13 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; stable = false end x = xt - # `f_κ` is convex and the dense and Woodbury steps solve its quadratic model - # exactly, so a whole step that leaves the violated set unchanged has reached - # the stage's minimizer: the gradient there is the model's, which is zero. - # The `:lsqr` solves are inexact and carry no such guarantee. + # For exact inner solves, an unchanged violation pattern ends the stage. !use_lsqr && stable && break fcur - fnew <= 5000 * eps(T) * max(fcur, one(T)) && break fcur = fnew end end - # Uniform boost to exact feasibility: x[p] + x[q] ≥ log|A_ij| on the support. - # `boost=false` leaves the iterate untouched, for the soft objective, which - # imposes no coverage constraint and whose optimum the boost would move off. + # Hard covers receive a final uniform feasibility shift. if boost γ = _boost_shift(x, supp, symmetric) for p in 1:N @@ -1250,40 +1010,20 @@ function _abslog2_continuation(sys::SupportSystem{T}, x0; linsolve=(use_lsqr ? :lsqr : use_woodbury ? :woodbury : :dense)) end -# Worker for `symcover_min(::AbsLog{2})`. Returns `(a, stats)` where `stats` is a -# NamedTuple `(; nsolves, lsqriters, cgiters, cholsolves, linsolve)` recording the -# number of inner linear solves, the total LSQR and conjugate-gradient iterations (0 on -# paths that run neither), how many Woodbury solves fell to the sparse factorization, -# and which path ran. -# `linsolve` reports the path that ran: `:dense`, `:woodbury`, or `:lsqr`. -# A working type narrower than `Float64` is solved in `Float64` and the cover converted -# back, since the continuation's tolerances assume double precision. -# `start`, when given, is a positive cover of `A` -# indexed like `axes(A, 1)` and supplies the first iterate in place of the cold -# unweighted solve; the objective is convex, so it changes the path but not the result. +# Worker for `symcover_min(::AbsLog{2})`, returning `(a, stats)`. A supplied +# `start` replaces the cold initial solve. Narrow types compute in `Float64`. function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), maxiter::Int=40, linsolve::Symbol=:auto, start=nothing, boost::Bool=true, fname=:symcover_min) linsolve in (:auto, :dense, :lsqr, :woodbury) || throw(ArgumentError("linsolve must be :auto, :dense, :lsqr, or :woodbury; got :$linsolve")) - # The shared entry to the native solve, reached from every sym `*_min` method, - # so the precondition is checked once here rather than at each of them. + # Shared symmetry check for native symmetric minimal covers. require_abs_symmetric(A, fname) ax = axes(A, 1) axes(A, 2) == ax || throw(ArgumentError("symcover_min requires a square matrix")) - # The problem only ever depends on abs.(A), a real quantity, so the working type - # stays real even for complex A (e.g. a complex Hermitian) — Complex has no total - # order, and the reweighted Newton solve below compares residuals with `<`/`min`. + # The problem depends only on `abs.(A)`, so the working type is real. T = float(real(eltype(A))) - # The continuation's tolerances are multiples of `eps(T)` — the decrease test at - # `5000*eps(T)`, the line-search floor at `500_000*eps(T)` — and the objective - # `f_κ` itself must resolve differences of that order at κ up to 1e8. Both assume - # double precision: at `eps(Float32)` the stages past the first carry no - # resolvable descent, and the continuation halts far from the constrained optimum. - # A narrower working type therefore runs the whole solve in `Float64` and the - # cover is returned in the caller's type. `convert` keeps the wrapper — the - # `Symmetric`, `Hermitian`, sparse and structured storage all have their own - # support traversals — and widens a complex eltype to `ComplexF64`. + # Continuation tolerances require at least Float64 resolution. if eps(T) > eps(Float64) a64, stats = _symcover_min_abslog2(convert(AbstractMatrix{promote_type(eltype(A), Float64)}, A); κs, maxiter, linsolve, start, boost, fname) @@ -1291,9 +1031,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end n = length(ax) use_lsqr = linsolve === :lsqr - # The support is measured before it is laid out, so that only the layout the chosen - # path needs is built. The Newton solve runs on 1-based positions 1:n and is - # scattered back onto `a` through `ax` so `A`'s own axes are honored. + # Build only the support layout needed by the chosen solver. G = _sym_support(A, T) hassupp = falses(n) nsupp = 0 # support entries of `A`, counted in both orientations @@ -1304,14 +1042,8 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), nsupp += ns maxzero = max(maxzero, n - ns) end - # The Woodbury path splits the normal equations around the complete-support matrix - # `n·I + e·eᵀ`, so its cost is set by the zero set `Z` rather than by `n`. Two - # separate conditions gate it. Per row: `n·I − L_Z` is positive definite only - # while no row carries more than `n ÷ 4` zeros. In total: `Z` is materialized and - # then traversed by every matvec and every factorization, so the path is worth - # taking only while `|Z|` stays O(n) — a support that is merely thin per row can - # still carry Θ(n²) zeros, and the split would then be dense work under a name - # that promises otherwise. + # `n*I - L_Z` is positive definite only while no row carries more than + # `n ÷ 4` zeros; the total budget bounds cost. nzero = n * n - nsupp zbudget = 4 * n use_woodbury = false @@ -1326,10 +1058,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end use_woodbury = ok end - # The layout the chosen path reads: the Woodbury sweeps run on the grid, the LSQR - # and dense sweeps on the list. Both hold one element per unordered pair `{i, j}` — - # the residuals `z_ij = α_i + α_j - log|A_ij|` of a pair and its mirror are the - # same, and `SupportSystem` weights the stored entry for both. + # Woodbury uses a grid; dense and LSQR use an edge list. supp = if use_woodbury C = fill(T(-Inf), n, n) for (ip, i) in enumerate(ax) @@ -1353,8 +1082,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end EdgeList{T}(edges, cvals) end - # Zero set of `A`, one entry per unordered pair. It is the off-diagonal pattern of - # the sparse `C` the Woodbury path factorizes. + # Zero set defining the off-diagonal pattern of sparse `C`. zedges = Tuple{Int,Int}[] if use_woodbury mark = falses(n) @@ -1368,9 +1096,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), fill!(mark, false) end end - # Nothing here pins a gauge: `v0` is zero, so the dense path adds no rank-one term - # and the LSQR gauge row is inert. The ridge lifts what singularity remains — the - # null space of a bipartite support graph such as `[0 1; 1 0]`. + # The ridge handles singular symmetric support graphs. sys = SupportSystem{T}(n, supp, true, hassupp, use_woodbury ? fill(T(n), n) : T[], zedges, ones(T, use_woodbury ? n : 0, use_woodbury ? 1 : 0), @@ -1386,11 +1112,7 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), return a, stats end -# Worker for `cover_min(::AbsLog{2})`. Returns `(a, b, stats)` with `stats` a -# NamedTuple `(; nsolves, lsqriters, cgiters, cholsolves, linsolve)` (see -# `_symcover_min_abslog2`, whose promotion of narrow working types this shares). -# `start`, when given, is a positive cover `(a, b)` indexed like the rows and columns -# of `A`, supplying the first iterate in place of the cold unweighted solve. +# Worker for `cover_min(::AbsLog{2})`, returning `(a, b, stats)`. function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), maxiter::Int=40, linsolve::Symbol=:auto, start=nothing, boost::Bool=true) @@ -1398,19 +1120,9 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), throw(ArgumentError("linsolve must be :auto, :dense, :lsqr, or :woodbury; got :$linsolve")) axr = axes(A, 1) axc = axes(A, 2) - # The problem only ever depends on abs.(A), a real quantity, so the working type - # stays real even for complex A (e.g. a complex Hermitian) — Complex has no total - # order, and the reweighted Newton solve below compares residuals with `<`/`min`. + # The problem depends only on `abs.(A)`, so the working type is real. T = float(real(eltype(A))) - # The continuation's tolerances are multiples of `eps(T)` — the decrease test at - # `5000*eps(T)`, the line-search floor at `500_000*eps(T)` — and the objective - # `f_κ` itself must resolve differences of that order at κ up to 1e8. Both assume - # double precision: at `eps(Float32)` the stages past the first carry no - # resolvable descent, and the continuation halts far from the constrained optimum. - # A narrower working type therefore runs the whole solve in `Float64` and the - # cover is returned in the caller's type. `convert` keeps the wrapper — the - # `Symmetric`, `Hermitian`, sparse and structured storage all have their own - # support traversals — and widens a complex eltype to `ComplexF64`. + # Continuation tolerances require at least Float64 resolution. if eps(T) > eps(Float64) a64, b64, stats = _cover_min_abslog2(convert(AbstractMatrix{promote_type(eltype(A), Float64)}, A); κs, maxiter, linsolve, start, boost) @@ -1420,10 +1132,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), n = length(axc) N = m + n use_lsqr = linsolve === :lsqr - # Each support entry links a row position ip to a column position m+jp. Internal - # positions 1:m index rows, m+1:m+n index columns, and results are scattered back - # through axr/axc so A's axes are honored. The support is measured before it is - # laid out, so that only the layout the chosen path needs is built. + # Stack row positions before column positions; scatter results back to `A`'s axes. G = _row_support(A, T) nzrow = zeros(Int, m) # support entries per row, for the balance convention nzcol = zeros(Int, n) # ditto per column @@ -1438,10 +1147,8 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end hasrow = nzrow .> 0 hascol = nzcol .> 0 - # The Woodbury path splits the normal equations around the complete-support - # matrix `D + u_r·u_rᵀ + u_c·u_cᵀ`, so its cost is set by the zero set `Z` rather - # than by `N`, and `D − L_Z` is positive definite only while `Z` stays thin. The - # bound is taken against `min(m, n)`, the smaller of the two diagonal blocks. + # `min(m,n)*I - L_Z` is positive definite only while no row or column + # carries more than `min(m,n) ÷ 4` zeros; the total budget bounds cost. maxzero = 0 for ip in 1:m maxzero = max(maxzero, n - nzrow[ip]) @@ -1450,12 +1157,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), maxzero = max(maxzero, m - nzcol[jp]) end zbound = min(m, n) ÷ 4 - # Two separate conditions gate the path. Per row and column: `D − L_Z` is positive - # definite only while neither carries more than `min(m, n) ÷ 4` zeros. In total: - # `Z` is materialized and then traversed by every matvec and every factorization, - # so the path is worth taking only while `|Z|` stays O(m + n) — a support that is - # merely thin per row can still carry Θ(m·n) zeros, and the split would then be - # dense work under a name that promises otherwise. + # Enforce both per-axis and total zero-count limits. nzero = m * n - ne zbudget = 4 * max(m, n) use_woodbury = false @@ -1470,8 +1172,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end use_woodbury = ok end - # The layout the chosen path reads: the Woodbury sweeps run on the grid, the LSQR - # and dense sweeps on the list. Both hold `log|A_ij|` for each support entry. + # Woodbury uses a grid; dense and LSQR use an edge list. supp = if use_woodbury C = fill(T(-Inf), m, n) for (ip, i) in enumerate(axr) @@ -1491,9 +1192,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), end EdgeList{T}(edges, cvals) end - # Zero set of `A` as stacked-position pairs missing from the support: the - # off-diagonal pattern of the sparse `C` the Woodbury path factorizes. The - # complete-support diagonal it corrects is `n` on rows and `m` on columns. + # Zero set defining the off-diagonal pattern of sparse `C`. zedges = Tuple{Int,Int}[] dfull = zeros(T, use_woodbury ? N : 0) Umat = zeros(T, use_woodbury ? N : 0, use_woodbury ? 2 : 0) @@ -1517,11 +1216,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), fill!(mark, false) end end - # Row and column scales share the global (e; −e) gauge, which every path pins - # through `v0`: ±1 on supported variables, 0 on support-free ones (which carry no - # constraint and are decoupled with an identity row instead). After the solve a - # closed-form shift, applied within each component, moves the result to the balance - # convention, so the pinned gauge is not observable. + # Pin the global row/column gauge on supported variables. v0 = zeros(T, N) for ip in 1:m hasrow[ip] && (v0[ip] = one(T)) @@ -1544,12 +1239,7 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), s0 end x, stats = _abslog2_continuation(sys, x0; κs, maxiter, linsolve, boost) - # Shift along the (e; -e) gauges to the balance convention ∑ nzaᵢ αᵢ = ∑ nzbⱼ βⱼ, - # imposed within each connected component of the support: the gauge acts - # independently on each component, so a single global shift would leave the - # per-component splits wherever the ridge (or LSQR's gauge row) put them. It - # applies whether or not the iterate was boosted: the gauge is a convention, not a - # constraint, and every cover this package returns satisfies it. + # Apply the balance convention independently to each support component. rowcomp, colcomp, ncomp = _support_components(A) Lα = zeros(T, ncomp) Lβ = zeros(T, ncomp) @@ -1581,37 +1271,21 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), return a, b, stats end -# Workers for the soft (unconstrained) AbsLog{2} covers. The soft objective -# `∑_{i,j∈S} (log a_i + log a_j - log|A_ij|)²` is the hard workers' reweighted -# least-squares problem with every weight held at 1, which is the cold solve they -# already take as their first iterate: `κs=()` runs no penalty continuation, and -# `boost=false` keeps the unconstrained minimizer where it is. It is convex, so one -# linear solve settles it — no iteration and no multistart, unlike the non-convex -# `AbsLinear` soft covers. -# -# Both paths inherit the hard workers' handling of a singular signless Laplacian (the -# `[0 1; 1 0]` support graph among them) and of support-free rows and columns. +# Soft `AbsLog{2}` covers are the unweighted initial solve with no feasibility +# shift. Convexity makes continuation and multistart unnecessary. _soft_symcover_min_abslog2(A::AbstractMatrix; kwargs...) = _symcover_min_abslog2(A; κs=(), boost=false, fname=:soft_symcover_min, kwargs...) _soft_cover_min_abslog2(A::AbstractMatrix; kwargs...) = _cover_min_abslog2(A; κs=(), boost=false, kwargs...) -# Internal exact reference implemented by the MatrixCoversJuMPExt extension; used only to -# cross-check the native `symcover_min(::AbsLog{2})` in the test suite. +# JuMP reference used to test the native symmetric solver. function symcover_min_jump end -# Internal exact reference implemented by the MatrixCoversJuMPExt extension; used only to -# cross-check the native `cover_min(::AbsLog{2})` in the test suite. +# JuMP reference used to test the native asymmetric solver. function cover_min_jump end -# A solve that stops for any reason other than a solved one leaves the model holding -# a point that does not solve the problem posed -- the base of an unbounded ray, or -# whatever the solver last had. Handing that back would be a minimal cover in name -# only, so it is an error. `ALMOST_*` statuses are rejected along with the rest: -# they report a tolerance the caller did not ask for. -# -# Takes the status rather than the model so that both solver extensions can share it -# without the main package depending on JuMP. +# Reject all non-solved statuses, including `ALMOST_*`. Taking the status keeps +# this helper independent of JuMP. function check_solved(status, solver, fname) Symbol(status) in (:OPTIMAL, :LOCALLY_SOLVED) || error("$fname: $solver terminated with status $status") diff --git a/src/penalties.jl b/src/penalties.jl index c30d6b4..4488d81 100644 --- a/src/penalties.jl +++ b/src/penalties.jl @@ -5,14 +5,11 @@ """ AbstractCoverPenalty <: Function -Supertype of the penalty functions `ϕ` that score a cover, and the type of the -first argument of most of this package's API. The built-in subtypes are -[`AbsLog`](@ref) and [`AbsLinear`](@ref). +Supertype of cover penalties. Built-in subtypes are [`AbsLog`](@ref) and +[`AbsLinear`](@ref). -A penalty is a function of the single ratio `r = |A[i,j]| / (a[i]*b[j])`, and -[`cover_objective`](@ref) sums it over the entries of `A`. Because `ϕ` sees only -that ratio, and every diagonal rescaling of `A` leaves it fixed, any objective -built from a penalty is automatically scale-invariant. +[`cover_objective`](@ref) applies the penalty to +`r = |A[i,j]|/(a[i]*b[j])` and sums over `A`. # Extending @@ -20,9 +17,8 @@ A subtype must be callable on a nonnegative real: (::MyPenalty)(r::Real) -`r` ranges over `[0, Inf]`. Both endpoints occur and neither may error: `r = 0` -whenever `A[i,j]` is zero, and `cover_objective` passes `typemax` for an entry -left uncovered by a zero scale. Penalties are conventionally singleton structs. +The method must accept `r = 0` and `r = typemax(...)`. Penalties are usually +singleton structs. [`cover_objective`](@ref) works for any subtype, but solvers support only specific built-in penalties: `AbsLog{2}` natively and `AbsLinear` through JuMP. @@ -38,13 +34,8 @@ Penalty type for φ(r) = |log(r)|^p if r > 0 0 if r = 0 -The discontinuity at r=0 prevents zero entries in A from sending the objective -value to infinity. - -This leads to convex optimization problems in log space. `AbsLog{1}` typically -has a flat minimum-basin in which members of an entire family of solutions are -equally good. `AbsLog{2}`, except in degenerate cases like `[0 1; 1 0]`, has a -unique minimum. +The `r=0` convention keeps zero entries finite. The objective is convex in log +space; `AbsLog{1}` may have multiple minima. See also: [`AbsLinear`](@ref). """ @@ -72,15 +63,8 @@ struct AbsLinear{p} <: AbstractCoverPenalty end """ MatrixCovers.scalar_type(T) -The plain floating-point type underlying the element type `T`, with any units -removed. [`cover_objective`](@ref) sums the ratios `|A[i,j]| / (a[i]*b[j])`, -which are dimensionless because a cover requires -`unit(A[i,j]) == unit(a[i])*unit(b[j])`, so the score is an ordinary number -whatever the operands carry. - -This cannot be expressed as `float(real(T))`: a matrix whose entries carry -different units has an abstract `eltype`, for which `real` and `oneunit` are -undefined. A unit-carrying element type therefore needs its own method. +Return the unitless floating-point type underlying `T`. Unit-carrying element +types should specialize this method. """ scalar_type(::Type{T}) where {T<:Number} = float(real(T)) @@ -99,16 +83,11 @@ end cover_objective(ϕ, a, b, A) cover_objective(ϕ, a, A) -Compute the cover objective `∑_{i,j} ϕ(|A[i,j]| / (a[i] * b[j]))` for the given -penalty function `ϕ`. The two-argument form is for symmetric matrices where the cover -is `a*a'`. +Compute `∑ ϕ(|A[i,j]|/(a[i]*b[j]))`. The shorter form uses the symmetric cover +`a*a'`. -The sum runs over the full grid in both forms, so in the symmetric form each -off-diagonal pair contributes twice and each diagonal entry once. This weighting -is what the `sym` solvers minimize, so the score reported here is the quantity -they optimized; code that reads a symmetric matrix through -[`foreach_support_sym`](@ref), which reports each pair once, must apply the -factor of 2 itself to match. +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}`: zero entries contribute 0 (φ(0) = 0 by convention). @@ -116,8 +95,7 @@ Zero entries of `A` are handled according to `ϕ`: `eachindex(a)` must match `axes(A, 1)` and `eachindex(b)` must match `axes(A, 2)`. -`A` is read through [`foreach_support`](@ref), so the cost is proportional to the -support rather than to `length(A)` for a storage type that specializes it. +`A` is read through [`foreach_support`](@ref). See also: - Penalty types (options for `ϕ`): [`AbsLog`](@ref), [`AbsLinear`](@ref). @@ -143,8 +121,7 @@ function cover_objective(ϕ, a, b, A) # zero entry over a zero scale, which constrains nothing. They therefore share # one penalty value, and only their count is needed: zero for `AbsLog`, but a # nonzero constant for `AbsLinear`, which is continuous at `r = 0`. - # The guard is not just an optimization: a penalty that is infinite at zero - # would otherwise turn a fully-dense `A` into `0 * Inf`, i.e. `NaN`. + # The guard prevents `0 * Inf` for penalties that are infinite at zero. nzero = length(a) * length(b) - nsupport[] return iszero(nzero) ? s[] : s[] + nzero * T(ϕ(zero(T))) end diff --git a/src/soft_covers.jl b/src/soft_covers.jl index 05899e3..754fffa 100644 --- a/src/soft_covers.jl +++ b/src/soft_covers.jl @@ -1,6 +1,4 @@ -# Soft covers: `symcover`/`cover` variants that penalize under-coverage instead -# of forbidding it, plus the scale-covariant multistart and coordinate-descent -# machinery they use. +# Soft covers and their multistart/coordinate-descent implementations. # ============================================================ # Public interface @@ -10,11 +8,8 @@ a = soft_symcover(ϕ, A; maxiter=32, starts=5, σ=2.0, rng=MersenneTwister(0)) a = soft_symcover(A; maxiter=32, starts=5, σ=2.0, rng=MersenneTwister(0)) -Given a square matrix `A` assumed to be symmetric, return a vector `a` approximately -minimizing the soft-cover objective `∑_{i,j} ϕ(|A[i,j]| / (a[i]*a[j]))`. - -Unlike [`symcover`](@ref), there is no hard coverage constraint: `a[i]*a[j]` may be -less than `|A[i,j]|`, with violations penalized by `ϕ`. +Approximately minimize `∑ ϕ(|A[i,j]|/(a[i]*a[j]))` for symmetric `A`, without a +hard coverage constraint. Supported penalties are: @@ -25,10 +20,9 @@ Supported penalties are: - `AbsLinear{1}()`: weighted-median descent initialized from the `AbsLinear{2}` result. -For `AbsLinear`, `starts` controls the number of deterministic and perturbed -starting points. Perturbations have the form `a .* exp.(σ .* ξ)`, where `ξ` is -drawn from `rng`. The default RNG is reset for each call. Pass an explicit `rng` -to control reproducibility. `sigma` is an ASCII alias for `σ`. +For `AbsLinear`, `starts` controls the number of starting points and `σ` the +spread of log-normal perturbations. Pass `rng` for reproducibility. `sigma` is +an alias for `σ`. See also: [`symcover`](@ref), [`cover_objective`](@ref), [`soft_symcover_min`](@ref). @@ -58,7 +52,7 @@ soft_symcover(::AbsLog{2}, A::AbstractMatrix; kwargs...) = soft_symcover_min(Abs function soft_symcover(::AbsLog{1}, A::AbstractMatrix; maxiter::Int=20) require_abs_symmetric(A, :soft_symcover) - a = soft_symcover_min(AbsLog{2}(), A) # convex AbsLog{2} minimum: a good start + a = soft_symcover_min(AbsLog{2}(), A) _abslog1_iter!(a, A, maxiter) return a end @@ -79,8 +73,7 @@ function soft_symcover(::AbsLinear{1}, A::AbstractMatrix; maxiter::Int=20, kwarg ax = axes(A, 1) axes(A, 2) == ax || throw(ArgumentError("soft_symcover requires a square matrix")) require_abs_symmetric(A, :soft_symcover) - # Initialize from the AbsLinear{2} soft cover (a good starting point for AbsLinear{1}); - # the multistart is spent there, then the AbsLinear{1} weighted-median descent refines. + # Refine the `AbsLinear{2}` multistart result by weighted-median descent. a = soft_symcover(AbsLinear{2}(), A; maxiter=5, kwargs...) _abslinear1_iter!(a, A, maxiter) return a @@ -90,19 +83,14 @@ end a = soft_symcover!(ϕ, a, A; maxiter=...) a = soft_symcover!(a, A; maxiter=...) -Refine the starting point `a` into a symmetric soft cover of `A`, in place, and return it. -The no-ϕ form defaults to `AbsLinear{2}()`, matching [`soft_symcover`](@ref), whose -supported ϕ values these methods share. - -Unlike [`soft_symcover`](@ref), this method refines one caller-provided start. -Build it with [`initialize_symcover`](@ref) and `feasible=:none`. +Refine one symmetric soft-cover start in place. The no-ϕ form uses +`AbsLinear{2}()`. Build a start with [`initialize_symcover`](@ref) and +`feasible=:none`. -`a` must be finite and strictly positive on every row of `A` that carries support; scales -on rows carrying no support are inert, and are zero on output. Unlike [`symcover_min!`](@ref), -`a` need *not* cover `A` — the soft objective imposes no coverage constraint. +`a` must be finite and positive on supported rows. It need not cover `A`, and +unsupported scales are set to zero. -`maxiter` bounds the descent sweeps. Under `AbsLog{2}`, the unique minimizer is -independent of the start. +`maxiter` bounds the descent sweeps. See also: [`soft_symcover`](@ref), [`soft_symcover_min!`](@ref), [`initialize_symcover`](@ref), [`soft_cover!`](@ref). """ @@ -138,12 +126,8 @@ end a, b = soft_cover(ϕ, A; maxiter=200, starts=4, σ=2.0, rng=MersenneTwister(0)) a, b = soft_cover(A; maxiter=200, starts=4, σ=2.0, rng=MersenneTwister(0)) -Given a matrix `A`, return vectors `a` and `b` approximately minimizing the soft-cover -objective `∑_{i,j} ϕ(|A[i,j]| / (a[i]*b[j]))`. This is the asymmetric analog of -[`soft_symcover`](@ref). - -Unlike [`cover`](@ref), there is no hard coverage constraint: `a[i]*b[j]` may be less than -`|A[i,j]|`, with violations penalized by `ϕ`. +Approximately minimize `∑ ϕ(|A[i,j]|/(a[i]*b[j]))` without a hard coverage +constraint. This is the asymmetric form of [`soft_symcover`](@ref). Supported penalties are: @@ -154,15 +138,12 @@ Supported penalties are: - `AbsLinear{1}()`: alternating weighted-median updates initialized from the `AbsLinear{2}` result. -Rows or columns of `A` that are entirely zero receive scale `0`. As with [`cover`](@ref), -only the products `a[i] * b[j]` are determined by the problem; the split is fixed by the -balance convention `∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j]`, imposed within each connected -component of the support (the gauge acts independently on each). +Unsupported rows and columns receive zero scale. The factors use the balance +convention of [`cover_min`](@ref). -For `AbsLinear`, `starts` controls the number of deterministic and perturbed -starting points. Perturbations have spread `σ` and use `rng`. The default RNG is -reset for each call; pass one explicitly to control reproducibility. `sigma` is -an ASCII alias for `σ`. +For `AbsLinear`, `starts` controls the number of starting points and `σ` the +spread of perturbations. Pass `rng` for reproducibility. `sigma` is an alias for +`σ`. See also: [`cover`](@ref), [`soft_symcover`](@ref), [`cover_objective`](@ref). @@ -186,7 +167,7 @@ soft_cover(A::AbstractMatrix; kwargs...) = soft_cover(AbsLinear{2}(), A; kwargs. soft_cover(::AbsLog{2}, A::AbstractMatrix; kwargs...) = soft_cover_min(AbsLog{2}(), A; kwargs...) function soft_cover(::AbsLog{1}, A::AbstractMatrix; maxiter::Int=20) - a, b = soft_cover_min(AbsLog{2}(), A) # convex AbsLog{2} minimum: a good start + a, b = soft_cover_min(AbsLog{2}(), A) _abslog1_iter_asym!(a, b, A, maxiter) return _balance_cover!(a, b, A) end @@ -201,8 +182,7 @@ function soft_cover(ϕ::AbsLinear{2}, A::AbstractMatrix; maxiter::Int=200, start end function soft_cover(ϕ::AbsLinear{1}, A::AbstractMatrix; maxiter::Int=100, kwargs...) - # Spend the multistart on the AbsLinear{2} cover (a good basin selector), then refine with - # the AbsLinear{1} weighted-median descent. + # Refine the `AbsLinear{2}` multistart result by weighted-median descent. a, b = soft_cover(AbsLinear{2}(), A; maxiter=5, kwargs...) _abslinear1_iter_asym!(a, b, A, maxiter) return _balance_cover!(a, b, A) @@ -256,9 +236,8 @@ end a = soft_symcover_min(ϕ, A) a = soft_symcover_min(A) -Return the ϕ-minimal symmetric soft cover of `A`: minimizes `∑_{i,j} ϕ(|A[i,j]|/(a[i]*a[j]))` -with no coverage constraints. The no-ϕ form defaults to `AbsLinear{2}()`, matching -[`soft_symcover`](@ref). +Return a local minimum of `∑ ϕ(|A[i,j]|/(a[i]*a[j]))` without coverage +constraints. The no-ϕ form uses `AbsLinear{2}()`. Supported ϕ values and required extensions: - `AbsLog{2}()`: solved natively as linear least squares; `linsolve` has the same @@ -279,10 +258,7 @@ function soft_symcover_min(::AbsLog{2}, A::AbstractMatrix; kwargs...) return a end -# Multistart driver for the non-convex soft AbsLinear covers, the unconstrained counterpart -# of `symcover_min(::AbsLinear)`. The kernels (`soft_symcover_min!`) live in MatrixCoversIpoptExt; the -# menu and the selection are native. Starts are taken raw: a soft cover is under no -# obligation to cover `A`. +# Multistart driver for the symmetric Ipopt kernels. Starts need not cover `A`. function soft_symcover_min(ϕ::AbsLinear, A::AbstractMatrix; strategies=SYMCOVER_MIN_STRATEGIES) ax = axes(A, 1) axes(A, 2) == ax || throw(ArgumentError("soft_symcover_min requires a square matrix")) @@ -307,8 +283,7 @@ The no-ϕ form uses `AbsLinear{2}()`. `a` must be positive on supported rows; unsupported scales are zeroed. It need not cover `A`. Use `feasible=:none` with [`initialize_symcover`](@ref). -For `AbsLinear`, the result can depend on the start. The `AbsLog{2}` result is -unique. +For `AbsLinear`, the result can depend on the start. See also: [`initialize_symcover`](@ref), [`soft_symcover_min`](@ref), [`symcover_min!`](@ref). """ @@ -348,25 +323,16 @@ end a, b = soft_cover_min(ϕ, A) a, b = soft_cover_min(A) -Return the ϕ-minimal asymmetric soft cover of `A`: minimizes -`∑_{i,j} ϕ(|A[i,j]|/(a[i]*b[j]))` with no coverage constraints. This is the asymmetric -analog of [`soft_symcover_min`](@ref). The no-ϕ form defaults to `AbsLinear{2}()`, -matching [`soft_cover`](@ref). - -The row/column scales are pinned to the balance convention -`∑ nzaᵢ log a[i] = ∑ nzbⱼ log b[j]` (`nzaᵢ`, `nzbⱼ` = nonzero counts of row `i`, column -`j`), imposed within each connected component of the support, as in [`cover_min`](@ref): -the objective depends on `a` and `b` only through the products `a[i]*b[j]`, so without a -convention the split between them would be arbitrary. +Return a local minimum of `∑ ϕ(|A[i,j]|/(a[i]*b[j]))` without coverage +constraints. The no-ϕ form uses `AbsLinear{2}()`. Factors use the balance +convention of [`cover_min`](@ref). Supported ϕ values and required extensions: -- `AbsLog{2}()`: solved natively; the minimizer is unique. +- `AbsLog{2}()`: solved natively. - `AbsLinear{1}()`, `AbsLinear{2}()`: require JuMP and Ipopt. Each strategy in `strategies` is refined, and the best local minimum is returned. - `AbsLog{1}()`: not implemented. -The starting points and objective are scale-covariant, as is the selected result. - See also: [`soft_cover_min!`](@ref), [`soft_symcover_min`](@ref), [`soft_cover`](@ref). """ function soft_cover_min end @@ -377,8 +343,7 @@ function soft_cover_min(::AbsLog{2}, A::AbstractMatrix; kwargs...) return a, b end -# Multistart driver for the non-convex asymmetric soft covers, the unconstrained -# counterpart of `cover_min(::AbsLinear)`. The kernels (`soft_cover_min!`) live in MatrixCoversIpoptExt. +# Multistart driver for the asymmetric Ipopt kernels. function soft_cover_min(ϕ::AbsLinear, A::AbstractMatrix; strategies=COVER_MIN_STRATEGIES) isempty(strategies) && throw(ArgumentError("soft_cover_min: `strategies` must name at least one starting cover")) @@ -450,13 +415,10 @@ end # Internal helpers # ============================================================ -# Default seed for the multistart perturbation RNG. Callers wanting reproducibility across -# Julia versions (whose default RNG streams are not stable) should pass their own `rng`. +# Default multistart seed. Pass an explicit RNG for cross-version reproducibility. const _MULTISTART_SEED = 0 -# Resolve a Unicode keyword and its ASCII alias to a single value, falling back to -# `default` when neither is given. Passing both raises an error unless they agree, -# so one can never silently override the other. +# Resolve Unicode and ASCII keyword aliases; reject conflicting values. function _resolve_alias(primary, alias, default, primary_name::Symbol, alias_name::Symbol) primary === nothing && return alias === nothing ? default : alias alias === nothing && return primary @@ -465,8 +427,7 @@ function _resolve_alias(primary, alias, default, primary_name::Symbol, alias_nam return primary end -# Relative improvement required to replace the incumbent. This prevents -# roundoff-equivalent candidates from changing the selection after rescaling. +# Margin that prevents roundoff-equivalent candidates from replacing the incumbent. _multistart_switchtol(::Type{T}) where {T} = 5_000_000 * eps(T) # Symmetric AbsLinear{2} starts in selection order, followed by log-normal @@ -515,15 +476,8 @@ function _multistart_select(objs) return besti end -# Shared driver for the AbsLinear{2} soft-cover multistarts: build the candidate list with -# `inits_builder`, refine each candidate in place with `iterate!`, score it with `objective`, -# and return the `_multistart_select` winner. `iterate!`/`objective` take the candidate itself -# (a bare vector for the symmetric cover, an `(a, b)` tuple for the asymmetric one) so the same -# driver serves both shapes. -# -# For cheap provenance auditing (which initialization earned the selection), pass `labels` -# and/or `objs` as empty vectors: they are filled in place with every candidate's label and -# final objective, in candidate order, so the winner is `labels[_multistart_select(objs)]`. +# Shared multistart driver. Optional `labels` and `objs` collect candidate data +# for tests. function _multistart_run(inits_builder::F, iterate!::G, objective::H, A::AbstractMatrix, iter::Int, starts::Int, σ::Real, rng; labels=nothing, objs=nothing) where {F,G,H} labs, inits = inits_builder(A, starts, σ, rng) @@ -536,12 +490,7 @@ function _multistart_run(inits_builder::F, iterate!::G, objective::H, A::Abstrac return inits[_multistart_select(E)] end -# Scale-covariant multistart for the symmetric AbsLinear{2} soft cover. Runs the single-start -# coordinate descent `_abslinear2_iter!` from the candidate list built by -# `_soft_symcover_abslinear2_inits` and returns the candidate `_multistart_select` picks. Every -# start co-varies with a diagonal rescaling `D*A*D` and the objective is scale-invariant, so the -# selection is scale-covariant; passing the same `rng` state across the two frames (as the -# default fresh-seeded RNG does) makes it reproducible. +# Scale-covariant symmetric `AbsLinear{2}` multistart. function _soft_symcover_abslinear2(A::AbstractMatrix, iter::Int, starts::Int, σ::Real, rng; labels=nothing, objs=nothing) return _multistart_run(_soft_symcover_abslinear2_inits, @@ -550,22 +499,12 @@ function _soft_symcover_abslinear2(A::AbstractMatrix, iter::Int, starts::Int, σ A, iter, starts, σ, rng; labels, objs) end -# Coordinate-descent iteration for AbsLinear{2} soft cover. -# Each coordinate a[k] is updated to the exact minimizer of +# Coordinate descent for a symmetric `AbsLinear{2}` soft cover. Each update +# minimizes # ½(1 - d/x²)² + ∑_{j≠k} (1 - c_j/x)² -# where d = |A[k,k]| and c_j = |A[k,j]|/a[j]. This is half the part of the -# `cover_objective` sum that depends on a[k]: that sum runs over the full grid, -# so each off-diagonal pair contributes twice and the diagonal once. -# Closed form when d=0 (x = s2/s1); Newton on a cubic otherwise. -# -# `iter` bounds the sweeps; the descent exits early once every coordinate's -# stationarity residual r_k = ∑_j (1 - ρ)ρ (ρ = |A[k,j]|/(a[k]a[j])), half the -# gradient of the objective above in log a[k], has magnitude below `tol` at the -# start of a sweep. The residual is available for free from the sums already formed -# for the update (r_k = s1/a[k] - s2/a[k]² + d/a[k]² - d²/a[k]⁴), it is the exact -# quantity optimality demands be zero, and it is scale-invariant (each ρ is), so -# covariant restarts of a rescaled problem exit on the same sweep and the -# multistart selection stays covariant. +# where `d = |A[k,k]|` and `c_j = |A[k,j]|/a[j]`. The `d=0` case is closed form; +# otherwise safeguarded Newton solves a cubic. A scale-invariant stationarity +# residual controls early exit. function _abslinear2_iter!(a::AbstractVector{T}, A::AbstractMatrix, iter::Int; tol::Real=50_000_000 * eps(T)) where T ax = eachindex(a) ax == axes(A, 1) || throw(DimensionMismatch("row indices of `A` must match `a`, got $(string(axes(A, 1))) vs $(string(ax))")) @@ -598,12 +537,8 @@ function _abslinear2_iter!(a::AbstractVector{T}, A::AbstractMatrix, iter::Int; t elseif iszero(d) x = s2 / s1 else - # The cubic g(x) = s1*x³ + (d - s2)*x² - d² has exactly one positive root - # (one Descartes sign change), and it is bracketed by √d and s2/s1: - # g(√d) = d*(s1*√d - s2) and g(s2/s1) = d*(s2²/s1² - d) have opposite signs. - # Safeguarded Newton with geometric bisection: the bracket endpoints can be - # separated by hundreds of orders of magnitude, so fallback steps must bisect - # in log space to converge in O(60) iterations. + # The cubic has one positive root bracketed by `sqrt(d)` and + # `s2/s1`; geometric bisection handles wide dynamic range. lo, hi = minmax(sqrt(d), s2 / s1) x = sqrt(lo * hi) while hi - lo > 2 * eps(hi) @@ -626,12 +561,8 @@ function _abslinear2_iter!(a::AbstractVector{T}, A::AbstractMatrix, iter::Int; t return a end -# Weighted median of the values in `c` using the values themselves as weights: the -# point `m` with ∑_{cᵢm} cᵢ. This minimizes ∑ᵢ |1 - cᵢ/x| over x > 0 — -# the gradient is (1/x²)(∑_{cᵢx} cᵢ), whose positive 1/x² factor leaves -# the root at that balance point. Sorts `c` in place and returns the lower weighted -# median, a deterministic, scale-covariant tie-break on the flat basin the AbsLinear{1} -# objective admits. +# Weighted median using each value as its weight. Sorts in place and returns the +# lower median as a deterministic, covariant tie-break. function _weighted_self_median!(c::AbstractVector{T}) where T sort!(c) half = sum(c) / 2 @@ -645,26 +576,12 @@ function _weighted_self_median!(c::AbstractVector{T}) where T return wm end -# AbsLinear{1} coordinate objective at candidate `x`: |1 - d/x²| + 2∑ᵢ |1 - cᵢ/x|, -# the part of the full-grid `cover_objective` sum that depends on a[k] (each -# off-diagonal pair appears twice there, the diagonal once). -# This is a top-level function because a closure in `_abslinear1_iter!` would -# capture and box the reassigned `d`, allocating in the inner loop and causing -# juliac's trim verifier to report a dynamic call. +# Symmetric `AbsLinear{1}` coordinate objective. A top-level function avoids a +# captured, boxed variable in the inner loop. _abslinear1_obj(x, d, c) = abs(1 - d/x^2) + 2 * sum(abs(1 - ci/x) for ci in c) -# Coordinate-descent iteration for AbsLinear{1} soft cover. -# Each coordinate a[k] is updated to reduce |1 - d/x²| + 2∑_{j≠k} |1 - c_j/x|, -# where d = |A[k,k]| and c_j = |A[k,j]|/a[j]. The off-diagonal sum is minimized -# by the weighted median of the c_j with weights c_j (the factor 2 does not move -# it). When d ≠ 0 the update takes the better of that median and sqrt(d) under -# `_abslinear1_obj`; this is not an exact coordinate minimization. -# -# `iter` bounds the sweeps; the descent exits early once the largest relative -# coordinate movement in a sweep drops to `tol`. The median update reaches an -# exact fixed point (identical ordering selects the same value), so movement -# falls to zero there. Relative movement is scale-invariant, so covariant -# restarts of a rescaled problem exit on the same sweep. +# Symmetric `AbsLinear{1}` coordinate descent. Each update chooses the better of +# the weighted median and `sqrt(d)`. Relative movement controls early exit. function _abslinear1_iter!(a::AbstractVector{T}, A::AbstractMatrix, iter::Int; tol::Real=5000 * eps(T)) where T ax = eachindex(a) ax == axes(A, 1) || throw(DimensionMismatch("row indices of `A` must match `a`, got $(string(axes(A, 1))) vs $(string(ax))")) @@ -710,19 +627,9 @@ function _abslinear1_iter!(a::AbstractVector{T}, A::AbstractMatrix, iter::Int; t return a end -# Coordinate-descent iteration for AbsLog{1} soft cover, working in log space (α = log a). -# Each coordinate α[k] is updated to minimize ∑_{j: A[k,j]≠0} |α[k] + α[j] - log|A[k,j]||. -# Holding the neighbours fixed, the minimizer over α[k] is the median of the points -# log|A[k,j]| - α[j], one per off-diagonal nonzero entry. The diagonal term is -# |2α[k] - log|A[k,k]|| = 2|α[k] - log|A[k,k]|/2|, i.e. the point log|A[k,k]|/2 with weight two, -# represented here by inserting it twice so a plain median carries the weighting. The AbsLog{1} -# minimum is a flat basin; the lower median is chosen for a deterministic, scale-covariant result. -# -# `iter` bounds the sweeps; the descent exits early once the largest relative -# coordinate movement in a sweep drops to `tol`. The median update reaches an -# exact fixed point, so movement falls to zero there. Relative movement is -# scale-invariant, so covariant restarts of a rescaled problem exit on the same -# sweep. +# Symmetric `AbsLog{1}` coordinate descent in log space. Each update is a +# weighted median; the diagonal contributes twice. Relative movement controls +# early exit. function _abslog1_iter!(a::AbstractVector{T}, A::AbstractMatrix, iter::Int; tol::Real=5000 * eps(T)) where T ax = eachindex(a) ax == axes(A, 1) || throw(DimensionMismatch("row indices of `A` must match `a`, got $(string(axes(A, 1))) vs $(string(ax))")) @@ -760,21 +667,9 @@ function _abslog1_iter!(a::AbstractVector{T}, A::AbstractMatrix, iter::Int; tol: return a end -# Alternating weighted-median descent for the asymmetric AbsLog{1} soft cover, working in -# log space (α = log a, β = log b). Updating α[i] with β fixed minimizes -# ∑_{j: A[i,j]≠0} |α[i] + β[j] - log|A[i,j]||, whose minimizer is the median of the points -# log|A[i,j]| - β[j], one per nonzero entry of row i; the β-update is dual. Row and column -# scales are distinct variables, so no term is self-coupled — the symmetric solver's -# double-weighted diagonal has no counterpart here — and each half-sweep is an exact block -# minimization. The AbsLog{1} minimum is a flat basin; the lower median is chosen for a -# deterministic, scale-covariant result. -# -# Each half-sweep minimizes exactly, but the objective's nonsmoothness couples α[i] with -# β[j], so a fixed point of the sweeps need not minimize it. See `soft_cover`. -# -# `iter` bounds the sweeps; the descent exits early once the largest relative coordinate -# movement in a sweep drops to `tol` (scale-invariant, so covariant restarts of a rescaled -# problem exit on the same sweep). Rows/columns with empty support keep scale 0. +# Asymmetric `AbsLog{1}` alternating weighted-median descent in log space. Each +# half-sweep minimizes its block, though a fixed point need not be a local +# minimum. Relative movement controls early exit. function _abslog1_iter_asym!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, iter::Int; tol=nothing) T = float(promote_type(eltype(a), eltype(b))) @@ -826,21 +721,14 @@ function _abslog1_iter_asym!(a::AbstractVector, b::AbstractVector, A::AbstractMa return a, b end -# Labeled `(a, b)` candidate starts for the asymmetric AbsLinear{2} multistart, in selection -# order. Deterministic starts: the boosted geometric mean (also the perturbation base) and the -# tightened hard cover, obtained by tightening a copy of the former so the shared passes run -# once. Remaining slots, up to `starts` total, are multiplicative log-normal perturbations -# `a_g .* exp.(σ .* ξ)`, `b_g .* exp.(σ .* η)` of that base, `ξ`/`η` drawn from `rng` (drawn -# for every index so the stream is frame-independent). +# Asymmetric `AbsLinear{2}` starts: boosted geometric mean, tightened cover, and +# log-normal perturbations. function _soft_cover_abslinear2_inits(A::AbstractMatrix, starts::Int, σ::Real, rng) T = float(real(eltype(A))) ag, bg = initialize_cover(A; strategy=:geomean, feasible=:boost) labels = ["boost"] inits = [(copy(ag), copy(bg))] - # `cover(A)` is this point tightened, then balanced and re-inflated. Tighten a copy - # (at `tighten_cover!`'s own default `maxiter`) rather than recomputing the shared - # geometric-mean and boost passes; the balance is gauge-only and the re-inflation - # only recovers roundoff, so this start differs from `cover(A)` negligibly. + # Reuse the geometric-mean and boost passes when constructing the hard start. length(inits) < starts && (push!(labels, "hardcover"); push!(inits, tighten_cover!(copy(ag), copy(bg), A))) k = 0 while length(inits) < starts @@ -858,52 +746,31 @@ function _soft_cover_abslinear2_inits(A::AbstractMatrix, starts::Int, σ::Real, return labels, inits end -# Scale-covariant multistart for the asymmetric AbsLinear{2} soft cover. Runs the single-start -# alternating least squares `_msmc_als!` from the candidate list built by -# `_soft_cover_abslinear2_inits` and returns the pair `_multistart_select` picks. Every start -# co-varies with an independent row/column rescaling `D_r*A*D_c` and the objective is -# scale-invariant, so the selection is scale-covariant; passing the same `rng` state across the -# two frames (as the default fresh-seeded RNG does) makes it reproducible. +# Scale-covariant asymmetric `AbsLinear{2}` multistart. function _soft_cover_abslinear2(A::AbstractMatrix, iter::Int, starts::Int, σ::Real, rng; labels=nothing, objs=nothing) - # Index the candidate pair in the body because juliac's trim verifier - # treats a tuple-destructuring lambda as a dynamic call. + # Explicit indexing avoids a dynamic tuple-destructuring call under `juliac`. a, b = _multistart_run(_soft_cover_abslinear2_inits, (ab, A, iter) -> _msmc_als!(ab[1], ab[2], A, iter), (ab, A) -> cover_objective(AbsLinear{2}(), ab[1], ab[2], A), A, iter, starts, σ, rng; labels, objs) - # The alternating half-sweeps rescale rows and columns independently, so they leave the - # gauge where it falls; pin it to the package's convention. The objective cannot see the - # gauge, so this changes no product a[i]*b[j] and no candidate's score. + # Balance the gauge after alternating row and column updates. return _balance_cover!(a, b, A) end -# Alternating least squares for the AbsLinear{2} soft cover in the inverse-scale variables -# u = 1 ./ a, v = 1 ./ b. With M = |A| restricted to its nonzero support, the objective -# E = ∑ (1 - M[i,j] u[i] v[j])² is biconvex; each half-sweep sets u[i] (resp. v[j]) to its -# exact minimizer. Rows/columns with empty support keep scale 0 and are held fixed. -# Refines `a`, `b` in place starting from their incoming values. -# -# Both half-sweeps accumulate over `i` with `j` held fixed, so the support is gathered -# by column once up front and each sweep walks that gather rather than the full grid. -# -# The post-sweep objective costs no extra pass over `A`: once the v-half-sweep has set -# v[j] = num[j]/den[j] from num[j] = ∑_i M[i,j] u[i] and den[j] = ∑_i (M[i,j] u[i])², -# column j contributes +# Alternating least squares for `AbsLinear{2}` in inverse scales `u=1/a`, +# `v=1/b`. Each half-sweep is exact. Column-grouped support is reused throughout. +# After updating `v[j] = num[j]/den[j]`, column `j` contributes # ∑_i (1 - M[i,j] u[i] v[j])² = nnz[j] - 2 v[j] num[j] + v[j]² den[j] # = nnz[j] - num[j]²/den[j] -# with nnz[j] the number of nonzeros in column j, which is that column's size in the -# gather. An empty column has den[j] = nnz[j] = 0 and contributes nothing. +# so the objective needs no extra matrix pass. function _msmc_als!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, iter::Int; tol=nothing) axr, axc = axes(A, 1), axes(A, 2) eachindex(a) == axr || throw(DimensionMismatch("row indices of `A` must match `a`, got $(string(axr)) vs $(string(eachindex(a)))")) eachindex(b) == axc || throw(DimensionMismatch("column indices of `A` must match `b`, got $(string(axc)) vs $(string(eachindex(b)))")) T = float(promote_type(eltype(a), eltype(b), real(eltype(A)))) - # The convergence test is on a relative movement, so its floor is set by the - # precision of `T`: a fixed Float64-scaled literal can never be reached in - # Float32 (every call would run to `iter`) and stops far short of what a wider - # type can resolve. + # Scale the convergence floor to the precision of `T`. rtol = tol === nothing ? 50 * eps(T) : T(tol) # Invert to inverse-scale variables; empty-support rows/columns (scale 0) stay at 0. u = map(x -> x > 0 ? inv(T(x)) : zero(T), a) @@ -952,9 +819,7 @@ function _msmc_als!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, ite return a, b end -# Soft-cover objective in inverse-scale variables: ∑_{i,j: A[i,j]≠0} (1 - |A[i,j]| u[i] v[j])². -# Takes the column-grouped support rather than the matrix, so the sweeps and the -# objective share one gather. +# Soft-cover objective in inverse scales, using column-grouped support. function _msmc_objective(C::GroupedSupport{T}, u::AbstractVector, v::AbstractVector) where T E = zero(T) for j in C.ax @@ -967,15 +832,8 @@ function _msmc_objective(C::GroupedSupport{T}, u::AbstractVector, v::AbstractVec return E end -# Alternating weighted-median descent for the asymmetric AbsLinear{1} soft cover. -# Updating a[i] with b fixed minimizes ∑_j |1 - |A[i,j]|/(a[i] b[j])| over a[i] > 0, whose -# minimizer is the weighted median of c_j = |A[i,j]|/b[j] weighted by the same c_j (see -# `_weighted_self_median!`); the b-update is dual. There is no self-coupled diagonal term (the -# (i,i) entry enters the row update through b[i] like any other column), so each full sweep is -# an exact block minimization and the objective decreases monotonically. Refines `a`, `b` in -# place from their incoming values; exits early once the largest relative coordinate movement -# in a sweep drops to `tol` (scale-invariant, so covariant restarts of a rescaled problem exit -# on the same sweep). Rows/columns with empty support keep scale 0. +# Asymmetric `AbsLinear{1}` alternating weighted-median descent. Each half-sweep +# minimizes its block; relative movement controls early exit. function _abslinear1_iter_asym!(a::AbstractVector, b::AbstractVector, A::AbstractMatrix, iter::Int; tol=nothing) T = float(promote_type(eltype(a), eltype(b))) diff --git a/src/sparse_support.jl b/src/sparse_support.jl index 5807b69..e863ef4 100644 --- a/src/sparse_support.jl +++ b/src/sparse_support.jl @@ -77,13 +77,8 @@ end # Native minimal-cover (MMC) solvers # ============================================================ -# Native AbsLog{2} MMC solvers on sparse supports default to the matrix-free LSQR -# inner solve, whose per-iteration cost is O(nnz) and whose accuracy tracks the -# conditioning of √W·R (≈ √κ) rather than that of the normal equations (≈ κ). This -# is the intended path when nnz ≪ n². Pass `linsolve=:dense` to force the dense -# factorization, or `linsolve=:auto` to let the solver choose between the Woodbury -# split and the dense factorization. Only AbsLog{2} is native; other penalties -# dispatch to the JuMP extension. +# Sparse `AbsLog{2}` solvers default to matrix-free LSQR. Use `linsolve=:dense` +# or `:auto` to request factorization-based paths. function symcover_min(ϕ::AbsLog{2}, A::SparseMatrixCSC; linsolve::Symbol=:lsqr, kwargs...) a, _ = _symcover_min_abslog2(A; linsolve, kwargs...) return a diff --git a/src/support.jl b/src/support.jl index 15189ab..2e6d796 100644 --- a/src/support.jl +++ b/src/support.jl @@ -1,24 +1,13 @@ -# Traversal of a matrix's stored support, shared by the cover heuristics -# (geometric-mean init, feasibility boost, tightening) so each is written once -# instead of once per storage type. -# -# Both are higher-order functions rather than iterators so that `f` is -# specialized and inlined into a tight loop at each call site; every index -# used is the matrix's own (`axes`, `eachindex`), so offset axes are honored. +# Traversal hooks for matrix support. Callbacks specialize at each call site, +# and indices follow the matrix axes. """ foreach_support(f, A) -Call `f(i, j, v)` once for every entry of `A` whose magnitude -`v = abs(A[i, j])` is nonzero, and return `nothing`. Entries that are zero are -skipped, so `f` never sees `v == 0`. The order is whatever suits `A`'s storage -and is not part of the contract; `i` and `j` are `A`'s own indices, so offset -axes are honored. +Call `f(i, j, abs(A[i,j]))` once per nonzero entry and return `nothing`. +Traversal order is unspecified; indices follow `axes(A)`. -This is the hook through which cover algorithms read a matrix. Specializing it -is what lets a storage type be covered in time proportional to its support -rather than to `length(A)` — the package's own `SparseMatrixCSC` methods, which -walk `nzrange` instead of the full grid, are the model. +Specialize this hook to support custom sparse storage in O(nnz) time. # Extending @@ -26,12 +15,8 @@ To support a new matrix type, define MatrixCovers.foreach_support(f, A::MyMatrix) -which must call `f(i, j, abs(A[i, j]))` exactly once for each `(i, j)` with -`abs(A[i, j]) != 0`, must not call `f` for any other entry (a stored zero is -still a zero), and must return `nothing`. Emitting an entry twice double-counts -it in the objective; omitting one silently drops a constraint, yielding a -"cover" that does not cover. Whatever `f` returns is ignored, so a traversal -runs to completion and must not be stopped early on the strength of it. +It must emit each nonzero entry exactly once, skip stored zeros, ignore callback +return values, and return `nothing`. See also: [`foreach_support_sym`](@ref). """ @@ -49,31 +34,15 @@ end foreach_support_sym(f, A) Symmetric counterpart of [`foreach_support`](@ref): call `f(i, j, v)` once per -unordered index pair rather than once per entry, and return `nothing`. Pairs are -reported in the canonical orientation `i <= j`, the diagonal included, with -`v = abs(A[i, j])`; zero pairs are skipped. `A` must be square, or a -`DimensionMismatch` is thrown. - -`abs.(A)` must also be **symmetric**, not merely square. That is what makes -reporting one member of each pair sufficient: a symmetric cover constrains -`a[i]*a[j]` by a single magnitude, so visiting `(j, i)` as well would only -duplicate it. Note the predicate is on the magnitudes, so a complex `Hermitian` -satisfies it — `|A[i,j]| == |conj(A[j,i])|`. +nonzero unordered pair in canonical order `i <= j`, including the diagonal. -This traversal does not check the precondition; the public `sym` entry points do, -before they call it (`MatrixCovers.require_abs_symmetric`). +`A` must be square and `abs.(A)` symmetric. Public symmetric solvers check this +precondition before calling the traversal. # Objective weighting -Because each pair is reported once, a caller accumulating a cover objective must -supply the multiplicity itself: `w = (i == j) ? 1 : 2`. That reproduces the -`∑_{i,j}` convention of [`cover_objective`](@ref), which runs over the full grid -and so counts each off-diagonal pair twice and each diagonal entry once. The -constraint set needs no such correction — `a[i]*a[j] >= |A[i,j]|` and its -transpose are the same constraint, so imposing it on the `i <= j` triangle alone -is equivalent to imposing it everywhere. Every solver in this package minimizes -the full-grid objective, so a cover's reported score and the quantity that was -minimized agree. +For full-grid objective weighting, use multiplicity 1 on the diagonal and 2 +off-diagonal. Constraints need no multiplicity. # Extending @@ -81,16 +50,9 @@ To support a new matrix type, define MatrixCovers.foreach_support_sym(f, A::MyMatrix) -which must call `f(i, j, v)` exactly once for each pair `i <= j` with -`v = abs(A[i, j]) != 0`, must not call `f` for zero pairs, and must return -`nothing`. Whatever `f` returns is ignored, so a traversal runs to completion -and must not be stopped early on the strength of it. -Reporting the same pair in both orientations double-counts it: the off-diagonal -weight of 2 is the caller's to apply, per *Objective weighting* above, so a pair -emitted twice is weighted 4. A -type whose storage is triangular (`Symmetric{<:Any,<:SparseMatrixCSC}` in this -package's own extension) must map stored `(i, j)` with `i > j` back to `(j, i)` -rather than emit it as found. +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`](@ref). """ @@ -117,16 +79,8 @@ const ASYMMETRY_ULPS = 8 """ MatrixCovers.require_abs_symmetric(A, fname) -Throw unless `abs.(A)` is symmetric to within roundoff, naming `fname` and the -first offending index pair. Return `nothing` otherwise. - -This is the precondition of [`foreach_support_sym`](@ref), enforced at the public -`sym` entry points rather than inside the traversal, which runs many times per -solve. An unchecked violation is not a visible failure: the cover returned would -be a cover of a symmetrization of `A`, plausible-looking and wrong. - -The predicate is on the magnitudes rather than on `A` itself, which is both what -the traversal reads and what admits a complex `Hermitian`. +Throw unless `abs.(A)` is symmetric to within roundoff. The error names `fname` +and the first offending pair. """ function require_abs_symmetric(A::AbstractMatrix, fname) ax = axes(A, 1) @@ -208,20 +162,9 @@ Connected components of a matrix's bipartite support graph, as returned by [`support_components`](@ref): one vertex per row and one per column, one edge per stored nonzero. -Component ids run `1:ncomponents(sc)`. A row or column of empty support belongs -to no component and reports `0`. Query an id with [`rowcomponent`](@ref) or -[`colcomponent`](@ref), which take the matrix's own indices, so offset axes need -no special case at the call site. - -The gauge orbit of an asymmetric cover has one dimension per component: the -rescaling `a -> γ*a`, `b -> b/γ` acts independently on each, because no product -`a[i]*b[j]` spans two components. Any convention pinning the split between `a` -and `b` must therefore be imposed per component; a single global constraint -leaves `ncomponents(sc) - 1` directions unpinned. - -Constructing this once and passing it to [`gramcover!`](@ref) lets a caller that -already knows the component structure — or that obtains it by some route other -than traversing a matrix — skip the traversal entirely. +Component ids run `1:ncomponents(sc)`; unsupported rows and columns report `0`. +Use [`rowcomponent`](@ref) and [`colcomponent`](@ref) with the matrix's own +indices. Pass this object to [`gramcover`](@ref) to reuse the traversal. """ struct SupportComponents{R<:AbstractUnitRange,C<:AbstractUnitRange} rowcomp::Vector{Int} @@ -234,9 +177,8 @@ end """ support_components(A) -> sc::SupportComponents -Connected components of the bipartite support graph of `A`, read through -[`foreach_support`](@ref) so a sparse storage type costs its support rather than -`length(A)`. +Return the connected components of `A`'s bipartite support graph. The matrix is +read through [`foreach_support`](@ref). See also: [`SupportComponents`](@ref). """ diff --git a/test/element_types.jl b/test/element_types.jl index 31b4f12..1ffb0ac 100644 --- a/test/element_types.jl +++ b/test/element_types.jl @@ -1,7 +1,4 @@ -# Element types other than Float64. The internal convergence tolerances are -# multiples of `eps(T)`, so both a narrower and a wider type must behave: a -# Float64-scaled literal is unreachable in Float32 (every descent would run to -# `maxiter`) and stops far short of what BigFloat can resolve. +# Element types other than Float64, including precision-scaled tolerances. @testset "element types" begin @@ -22,12 +19,7 @@ @test iscover(a, b, B; rtol=8eps(Float32)) end - # The AbsLog{2} penalty continuation resolves descent of order `eps(T)` at penalty - # strengths up to 1e8, which `Float32` cannot represent: carried out in `Float32` - # throughout, every stage past the first makes no progress and the cover lands tens - # of percent from the optimum. The solve therefore runs in `Float64` whenever the - # working type is narrower, so a narrow answer is the `Float64` answer rounded, and - # the element and container types still follow the input. + # Narrow inputs compute in Float64 and convert back to the requested type. @testset "narrow working types solve in Float64" begin rng = StableRNG(77) X = exp.(randn(rng, 60, 60)) @@ -94,9 +86,7 @@ end @testset "row and column scales may differ in element type" begin - # The public entry points take plain `AbstractVector`s, so a mismatch must - # not surface as a MethodError from an unexported internal. Each vector - # keeps its own element type; the arithmetic promotes. + # Public refiners accept mixed vector element types and promote arithmetic. A = [4.0 1.5 0.5; 1.5 1.0 2.0] a0, b0 = Float64[3.0, 3.0], Float32[3.0, 3.0, 3.0] @@ -119,15 +109,12 @@ end @testset "tolerances follow the element type" begin - # Each convergence threshold must sit above the type's resolution, or the - # test it guards can never fire. + # Convergence thresholds must exceed the type's resolution. for T in (Float32, Float64, BigFloat) @test MatrixCovers._multistart_switchtol(T) > eps(T) end - # And must sit below it for a wider type, so the extra precision is used. - # Running the ALS kernel from one start under the eltype-scaled tolerance - # and under the Float64-scaled literal it replaced separates the two. + # Wider types must use their extra precision. A = BigFloat[4 1.5 0.3; 1.5 1 0.7; 0.3 0.7 2.0] start = initialize_symcover(A; feasible=:none) u1, v1 = copy(start), copy(start) diff --git a/test/extensions.jl b/test/extensions.jl index 5743382..409b9f4 100644 --- a/test/extensions.jl +++ b/test/extensions.jl @@ -50,10 +50,7 @@ a_soft = soft_symcover_min(AbsLog{2}(), A_rank1) @test cover_objective(AbsLog{2}(), a_soft, A_rank1) < 1e-8 - # A solve that does not reach an optimum is an error, not a cover: the point - # such a model holds is the base of a ray rather than a minimizer. The guard is - # exercised directly on the status, since the symmetry precondition rejects the - # asymmetric input that is what left this LP unbounded. + # Non-optimal solver statuses are errors. @test MatrixCovers.check_solved(JuMP.OPTIMAL, "HiGHS", "symcover_min") === nothing @test MatrixCovers.check_solved(JuMP.LOCALLY_SOLVED, "Ipopt", "symcover_min!") === nothing @test_throws "terminated with status" MatrixCovers.check_solved(JuMP.DUAL_INFEASIBLE, "HiGHS", "symcover_min") @@ -117,9 +114,7 @@ end Aasym = [1.0 2.0 3.0; 4.0 5.0 6.0] for ϕ in PENALTIES - # Refining a start yields a hard cover no worse than the start itself, to - # within the tolerance the solvers converge to (on this matrix the :hardcover - # start is already all but optimal, so there is nothing else separating them). + # Refinement does not worsen the start beyond solver tolerance. a0 = initialize_symcover(A) a = symcover_min!(ϕ, copy(a0), A) @test iscover(a, A; rtol=1e-6) @@ -136,9 +131,7 @@ end @test ga ≈ ab && gb ≈ bb end - # AbsLog{1}'s optimum is a whole face of equally-scoring covers, but the solver pins the - # member minimizing the AbsLog{2} objective over it, so the start cannot be read off the - # result — the same point comes back from any of them. + # `AbsLog{2}` tie-breaking makes the `AbsLog{1}` result start-independent. a_cold = symcover_min(AbsLog{1}(), A) for strategy in (:hardcover, :geomean, :diagfeasible) @test symcover_min!(AbsLog{1}(), initialize_symcover(A; strategy), A) ≈ a_cold @@ -149,9 +142,7 @@ end @test ah ≈ ab_cold && bh ≈ bb_cold end - # The AbsLinear objectives are non-convex: on this matrix the :hardcover and - # :geomean starts descend into genuinely different local minima, which is what - # makes a menu of starts worth having. + # Different starts reach different `AbsLinear` local minima. Abasin = [81.892035218799 1.06622031288736 29.4700945830419 0.0181293142917846; 1.06622031288736 0.243512973596586 38.0236584552296 0.0279078887878805; 29.4700945830419 38.0236584552296 8.96405068596511 26.5775238859338; @@ -171,9 +162,7 @@ end end @testset "soft_symcover_min multistart and refiner (Ipopt)" begin - # A start for the soft cover need not cover `A`: the objective imposes no coverage - # constraint, and the raw geometric mean — the exact soft AbsLog{2} optimum — does not - # cover. The refiner must accept it where symcover_min! would reject it. + # Soft refiners accept non-covering starts. A = [4.0 2.0 1.0; 2.0 3.0 2.0; 1.0 2.0 5.0] a0 = initialize_symcover(A; strategy=:geomean, feasible=:none) @test !iscover(a0, A) @@ -197,11 +186,7 @@ end end end - # For a non-convex objective, scale covariance is a property of the *start*: a start that - # does not co-vary with `A` can reach a different basin once the frame is rescaled, and - # the objective is scale-invariant only within a basin. This matrix and rescaling separate - # the basins sharply enough to catch that — an `A`-independent start scores 4.00 in the - # rescaled frame where the co-varied answer scores 1.47 — so it pins the menu's covariance. + # Covariant starts must select corresponding nonconvex basins after rescaling. Abasins = [0.020358630644342735 0.53352144014074843 5.8899714528796077 0.23770314779348869 3.0721768720180109; 0.53352144014074843 3.5416395788642903 37.199280652497748 49.972622569225109 333.53567816710364; 5.8899714528796077 37.199280652497748 0.76014027958709862 2.4189759139690739 0.92571970793600067; @@ -263,10 +248,7 @@ end end @testset "AbsLog{1} canonical selection on the optimal face (HiGHS)" begin - # The AbsLog{1} optimum is a face of the feasible polytope, not a point: its members are - # different covers scoring the same objective. The solver returns the member that also - # minimizes the AbsLog{2} objective — L1-optimal still, but the tightest such cover - # rather than whichever vertex the LP happened to reach. + # Select the `AbsLog{2}`-minimal member of the `AbsLog{1}` optimal face. ϕ1, ϕ2 = AbsLog{1}(), AbsLog{2}() rng = StableRNG(17) for n in (3, 5, 8) @@ -293,9 +275,7 @@ end @test symcover_min(ϕ1, D * A * D) ≈ D * a end - # Asymmetric: the balance convention pins the gauge `a -> c*a, b -> b/c`, which leaves - # every product a[i]*b[j] untouched; the canonical selection resolves the L1 face, whose - # members have genuinely different products. The two rules are independent and both hold. + # Gauge balancing and `AbsLog{1}` face selection are independent. Aasym = [3.0 1.0 7.0; 2.0 5.0 1.0; 8.0 1.0 4.0] a, b = cover_min(ϕ1, Aasym) nza = vec(count(!iszero, Aasym, dims=2)) @@ -309,9 +289,7 @@ end end @testset "AbsLinear multistart drivers (Ipopt)" begin - # The matrix on which the starts genuinely disagree: :geomean reaches a better - # AbsLinear{2} minimum than :hardcover, so a driver that tried only the latter - # would report the worse of the two. + # `:geomean` reaches the better `AbsLinear{2}` basin here. Abasin = [81.892035218799 1.06622031288736 29.4700945830419 0.0181293142917846; 1.06622031288736 0.243512973596586 38.0236584552296 0.0279078887878805; 29.4700945830419 38.0236584552296 8.96405068596511 26.5775238859338; @@ -345,8 +323,7 @@ end cover_objective(ϕ, symcover_min(ϕ, Abasin), Abasin) * (1 + 1e-6) + 1e-8 end - # On Abasin the :geomean start wins, so restricting the menu to :hardcover is - # observable — the driver is refining the menu it is given, not a fixed start. + # Restricting `strategies` changes the selected basin. ϕ = AbsLinear{2}() @test symcover_min(ϕ, Abasin; strategies=(:hardcover,)) ≈ symcover_min!(ϕ, initialize_symcover(Abasin; strategy=:hardcover), Abasin) @@ -372,8 +349,7 @@ end end @testset "error hint gated on argument types" begin - # Wrong-argument-type MethodError: no extension load would fix this, - # so the hint must not fire. + # Do not hint when loading an extension cannot fix the argument type. A = [4.0 1.0; 1.0 4.0] e = try symcover_min(AbsLog{2}(), "not a matrix") @@ -384,9 +360,7 @@ end @test e isa MethodError @test !occursin("loading JuMP", sprint(showerror, e)) - # Genuine missing-extension MethodError: run in a fresh process with - # JuMP/HiGHS/Ipopt unloaded (this test file loads them itself, which - # would otherwise mask the failure), so the hint should fire. + # Test missing-extension hints in a fresh process. script = """ using MatrixCovers A = [4.0 1.0; 1.0 4.0] @@ -399,10 +373,7 @@ end out = read(`$(Base.julia_cmd()) --project=$(Base.active_project()) -e $script`, String) @test occursin("loading JuMP", out) - # The no-ϕ wrapper `soft_symcover_min(A)` exists in the base package, but the - # `AbsLinear{2}` method it forwards to lives in the MatrixCoversIpoptExt extension. The - # MethodError raised (and hinted on) is for the inner call, so the hint must - # still fire even though the outer, no-ϕ call is what the user wrote. + # The no-ϕ wrapper should preserve the inner missing-extension hint. script_noϕ = """ using MatrixCovers A = [4.0 1.0; 1.0 4.0] @@ -415,9 +386,7 @@ end out_noϕ = read(`$(Base.julia_cmd()) --project=$(Base.active_project()) -e $script_noϕ`, String) @test occursin("loading JuMP", out_noϕ) - # soft_cover_min's AbsLog{1} is genuinely unimplemented rather than gated behind an - # extension, so its hint must say so rather than claim a package load would help — while - # still advising the load for the AbsLinear penalties, which Ipopt does provide. + # Distinguish an unimplemented penalty from a missing extension. e3 = try soft_cover_min(AbsLog{1}(), A) nothing @@ -474,8 +443,7 @@ end @test !occursin("loading JuMP", sprint(showerror, e5)) end -# `HookOnlyMatrix` (defined in soft_covers.jl) throws from `getindex`, so a model -# builder that materialized `A` in position space could not get this far. +# `HookOnlyMatrix` verifies that model builders use support hooks. @testset "the solver entry points read through the support hook" begin entries = [(1, 1, 2.0), (1, 3, 1.5), (2, 2, 3.0), (2, 4, 0.5), (3, 4, 4.0), (4, 4, 1.0)] M = HookOnlyMatrix(entries, 4) @@ -519,12 +487,7 @@ end end end -# The `sym` AbsLinear solvers minimize the full-grid objective — each off-diagonal -# pair twice, each diagonal entry once — which is what `cover_objective` reports. -# This matrix discriminates between that and weighting each unordered pair once: -# the latter convention puts the optimum near [2.0, 3.177, 1.574] instead. Most -# matrices do not discriminate, because the binding constraints have zero residual -# at the optimum and a zero residual is weight-independent. +# This matrix distinguishes full-grid symmetric weighting from triangle weighting. @testset "cover_min balances a block-diagonal support (JuMP/HiGHS)" begin # Two connected components: the model constraint pins only the global gauge # direction, so the per-component balance must come from the post-solve diff --git a/test/gram_covers.jl b/test/gram_covers.jl index 823aeeb..fd1f061 100644 --- a/test/gram_covers.jl +++ b/test/gram_covers.jl @@ -46,16 +46,8 @@ end @testset "one component: the global bound is attained" begin - # With a single support component there is nothing to accumulate separately, so the - # per-component construction returns exactly the global pair bound `norm(a)*b`. The - # agreement is two-sided: an inequality alone would also pass for a construction that - # gave up slack it did not have to. - # - # It is not bitwise, and the deviation is one-sided by design. Each `sqrt` is inflated by - # `1 + (n+3)*eps` so the cover holds despite naive summation; `norm` accumulates its own - # `n` terms with no such margin. So the tolerance is `(2n+3)*eps` — the declared inflation - # plus the reference's own roundoff — and `all(s .<= norm(a) .* b)` is *false* by a few - # ulps. That inequality is exposition, not a test. + # A single component matches `norm(a)*b` within the documented roundoff + # inflation. for (seed, m, k) in ((3, 5, 4), (7, 4, 4), (11, 12, 3)) rng = StableRNG(seed) J = randn(rng, m, k) @@ -77,10 +69,7 @@ a, b = cover(J) s = gramcover(a, b, J) - # gramcover on the joint (a, b) restricted to a block equals gramcover on - # that block alone with the corresponding sub-vectors: the computation is - # a function of the connected component, and block-diagonal J has disjoint - # components per block. + # Each block-diagonal support component is independent. sB = gramcover(a[1:4], b[1:3], B) sC = gramcover(a[5:7], b[4:5], C) @test isapprox(s, vcat(sB, sC); rtol=1e-9) @@ -106,9 +95,7 @@ s2 = gramcover(a2, b2, J) @test isapprox(s, s2; rtol=1e-12) - # A `W` coupling the two components merges them, and the gauge then acts - # with a different `γ` on each half of every cross-block sum. The bound - # must still not depend on which gauge the solver happened to return. + # Coupled components remain invariant under independent input gauges. m = size(J, 1) for W in (Matrix(2.0I, m, m) + [i == 1 && ip == 5 for i in 1:m, ip in 1:m], fill(0.5, m, m) + Diagonal(1:m)) @@ -295,9 +282,7 @@ end @testset "two coupled components, both diagonal blocks nonzero" begin - # For a 2x2 `Ms` with positive diagonal, the minimal cover is - # `σ[p] = sqrt(Ms[p,p]) * sqrt(max(1, κ))`, `κ = Ms[1,2]/sqrt(Ms[1,1]*Ms[2,2])`: - # below `κ = 1` the diagonal constraints bind alone, above it the coupling does. + # Closed form for a 2×2 `Ms` with positive diagonal. rng = StableRNG(9) B = randn(rng, 4, 3); C = randn(rng, 3, 2) J = [B zeros(4, 2); zeros(3, 3) C] diff --git a/test/heuristic_covers.jl b/test/heuristic_covers.jl index 6f08e44..1a40f03 100644 --- a/test/heuristic_covers.jl +++ b/test/heuristic_covers.jl @@ -192,9 +192,7 @@ end dr, dc = exp.(randn(rng, n)), exp.(randn(rng, m)) @test covaries(A -> cover(AbsLog{2}(), A; maxiter=0), Ag, dr, dc; rtol=1e-10) - # Quality gate: median log-optimality-gap of the 3-iteration heuristic over - # this fixed corpus, with a generous 1.5x margin over the measured value; a - # tighter algorithm may lower it, a regression will trip it. + # Bound the median objective gap across the fixed corpus. qrng = StableRNG(20260708) gaps = Float64[] for _ in 1:15 @@ -226,11 +224,7 @@ end M = Matrix(T40) @test iscover(a, M; rtol=8eps()) - # Float32 dynamic range wide enough that linear-domain deficit ratios overflow. - # The boost's apply! step shifts log(a[i]) by h = z/2, where z ~ 120 for this - # matrix; exp(log(a[i]) + h) then carries forward the rounding error already - # present in h at that magnitude, so the achievable relative precision is set - # by eps(Float32) scaled by |h|, not by a fixed few-ulp bound. + # Float32 range where linear-domain deficit ratios overflow. A32 = fill(1f-35, 6, 6); A32[1, 2] = A32[2, 1] = 3f37 a32 = symcover(AbsLog{2}(), A32) @test all(isfinite, a32) diff --git a/test/initializers.jl b/test/initializers.jl index e88e625..470d2ab 100644 --- a/test/initializers.jl +++ b/test/initializers.jl @@ -89,10 +89,7 @@ end @testset "no penalty argument" begin - # Every start on the menu is a property of `A` alone, so none of them needs a ϕ and - # none is offered one. A regression check on that, not a promise it will never - # change: a ϕ-tuned start would arrive as a new method, with the fallback dropping ϕ - # and calling these — additive, so callers of the current forms are unaffected. + # Initializers depend on `A`, not on a penalty. @test_throws MethodError initialize_symcover(AbsLog{2}(), Asyms[1]) @test_throws MethodError initialize_cover(AbsLog{2}(), Aasyms[1]) end diff --git a/test/invariants.jl b/test/invariants.jl index e9788e6..cf3c0a2 100644 --- a/test/invariants.jl +++ b/test/invariants.jl @@ -1,9 +1,4 @@ -# Cross-notion invariants: conventions documented for every cover notion, -# checked uniformly across all of them. Each entry supplies the solver as -# `A -> a` (symmetric) or `A -> (a, b)` (general), whether it promises hard -# feasibility, and the tolerance its algorithm warrants. The themed files pin -# each notion's algorithm-specific precision; this file pins the shared -# conventions: +# Shared invariants across cover algorithms: # - repeated calls return identical results # - hard covers are feasible # - results co-vary with diagonal rescaling of A @@ -87,9 +82,7 @@ const GEN_NOTIONS = ( @test axes(ao, 1) == axes(Ao, 1) @test axes(bo, 1) == axes(Ao, 2) @test collect(ao) .* transpose(collect(bo)) ≈ a .* transpose(b) rtol=nt.rtol - # The gauge a -> c*a, b -> b/c is invisible to every objective and every coverage - # constraint, so nothing in the problem fixes the split between `a` and `b`. The - # balance convention does, and every asymmetric cover reports its result in it. + # Every asymmetric cover uses the balance convention. @test isbalanced(a, b, Agen) @test isbalanced(az, bz, Azgen) # The balance convention is imposed per connected component: on a support with @@ -105,11 +98,7 @@ const GEN_NOTIONS = ( @test isbalanced(initialize_cover(Ablk; strategy, feasible)..., Ablk) end - # The gauge factor is a whole power of two, so imposing the convention is exact - # in binary floating point: every product the coverage constraints see is - # preserved bit for bit. A cover cannot be perturbed into infeasibility by the - # act of pinning its gauge. Products across two components are not preserved, - # and are not constrained either — the support is exactly where both hold. + # Power-of-two gauge shifts preserve supported products exactly. @testset "balancing preserves on-support products exactly" begin rng = StableRNG(11) for A in (Agen, Ablk, Azgen) @@ -143,13 +132,7 @@ const GEN_NOTIONS = ( @test bbd == vcat(b1, b2) end - # The sym objective is summed over the full grid: each off-diagonal pair counts - # twice, each diagonal entry once. `cover_objective` is the reference, and every - # sym solver must minimize that same weighting even though its constraints live on - # the `i <= j` triangle. The references below impose the constraints on the full - # grid explicitly, so agreeing with them pins both halves of the convention: that - # the triangle is the equivalent constraint set, and that the objective is not - # halved along with it. + # Symmetric solvers use full-grid objective weights with triangle constraints. @testset "sym solvers minimize the full-grid objective" begin function ref_min(pow, A) n = size(A, 1) @@ -167,10 +150,7 @@ const GEN_NOTIONS = ( return JuMP.objective_value(model) end - # The two conventions share a minimizer whenever the diagonal residuals - # vanish at the optimum, which is the common case — so a discriminating - # matrix is committed rather than left to the random draws. Weighting this - # one's off-diagonal pairs once instead of twice moves the optimum. + # This matrix distinguishes full-grid from triangle weighting. Mdisc = [2.4 1.2 1.2; 1.2 1.4 2.4; 1.2 2.4 0.6] rng = StableRNG(20) diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index b71f141..a95b131 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -1,7 +1,6 @@ # The *_min family: native AbsLog{2} minimal-cover solvers and their edge cases. -# Committed 5x5 matrix libraries (`symmetric_matrices`, `general_matrices`); the guard -# permits re-inclusion of this file in an already-initialized session. +# Reuse the committed 5×5 matrix libraries when already loaded. if !isdefined(@__MODULE__, :symmetric_matrices) include("testmatrices.jl") end @@ -10,8 +9,7 @@ end # Non-square rejected. @test_throws "symcover_min requires a square matrix" symcover_min(AbsLog{2}(), [1.0 2.0; 3.0 4.0; 5.0 6.0]) - # Native solver matches the HiGHS reference in objective across the whole - # committed symmetric library, and returns a feasible cover. + # Match HiGHS across the symmetric corpus. for (_, A) in symmetric_matrices Af = Float64.(A) a = symcover_min(AbsLog{2}(), Af) @@ -48,9 +46,7 @@ end end @testset "cover_min native AbsLog{2}" begin - # Native solver returns a feasible cover across the whole committed general - # library, and matches the HiGHS reference in objective on a deterministic - # subsample (the full 4367-matrix JuMP cross-check is slow). + # Match HiGHS on a deterministic sample of the general corpus. idx_sub = Set(round.(Int, range(1, length(general_matrices), length=500))) for (k, (_, A)) in enumerate(general_matrices) Af = Float64.(A) @@ -143,11 +139,8 @@ end @test a[2] * b[1] ≈ 1.0 end -# The Woodbury path solves the same regularized normal equations as the dense path, -# through a sparse Cholesky of `C` and a low-rank update, so the two must agree to -# roundoff. The tolerance is loose relative to `eps`: a converged cover pins the -# objective far more tightly than its own entries, so the two solves separate at -# roughly the square root of the working precision. +# Dense and Woodbury solve the same regularized equations. Factor comparisons use +# `sqrt(eps)` because the objective is more tightly determined than the factors. @testset "MMC native AbsLog{2} Woodbury path" begin rng = StableRNG(9) lognormal(m, n) = exp.(randn(rng, m, n)) @@ -155,10 +148,7 @@ end @testset "symmetric, n = $n" for n in (6, 30, 120) A = symlognormal(n) - # A symmetric zero set placed to sit inside both guards at every size tested: - # pairing consecutive indices gives exactly one zero per row, against a - # per-row allowance of `n ÷ 4` (which is 1 already at n = 6) and a total - # allowance of `4n`. + # One zero per row satisfies both Woodbury sparsity guards. Z = symlognormal(n) for k in 1:(n ÷ 2) Z[2k-1, 2k] = 0.0 @@ -167,8 +157,7 @@ end variants = ["all nonzero" => A, "zero diagonal" => A - Diagonal(A), "paired zeros" => Z] - # Diagonal and off-diagonal zeros together put two per row, which needs a - # per-row allowance of at least two. + # Diagonal and off-diagonal zeros give two per row. n ÷ 4 >= 2 && push!(variants, "zero diagonal and paired zeros" => Z - Diagonal(Z)) for (name, M) in variants @testset "$name" begin @@ -181,9 +170,7 @@ end @test aw ≈ ad rtol=1e-7 @test aa == aw @test iscover(aw, M; atol=1e-8) - # Both Woodbury sub-paths run within a continuation: the early stages - # are well enough conditioned for conjugate gradients, the late ones - # are not, and both are exact. + # Exercise both CG and factorized Woodbury solves. @test sw.cgiters > 0 @test sw.cholsolves > 0 @test sd.cgiters == 0 @@ -308,11 +295,8 @@ end @test MatrixCovers._cover_min_abslog2(Gbig)[3].linsolve === :dense end -# A Newton step is exact on the dense and Woodbury paths, so a whole step that leaves -# the violated set unchanged has already reached the minimizer of the current -# penalty stage, and the stage ends without a confirmation solve. The `:lsqr` steps -# are inexact and keep the decrease test as their sole criterion, which is what makes -# their solve counts the reference here. +# Exact inner solves stop a stage when the violated set is unchanged; LSQR uses +# the decrease test. @testset "MMC exact paths stop on a sign-stable Newton step" begin rng = StableRNG(31) A = (X = exp.(randn(rng, 60, 60)); (X .+ X') ./ 2) @@ -321,10 +305,7 @@ end al, sl = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr) @test ad ≈ al rtol=1e-6 @test aw ≈ al rtol=1e-6 - # One solve per κ stage is saved; `κs` has four stages by default. The absolute - # bound guards against a regression in the count itself: this matrix takes 24 - # solves on the exact paths against 28 on `:lsqr`, so 26 leaves two solves of - # headroom while still failing if the early stop stops firing. + # Exact paths save one solve per continuation stage. @test sd.nsolves == sw.nsolves @test sd.nsolves <= sl.nsolves - length((1e2, 1e4, 1e6, 1e8)) @test sd.nsolves <= 26 @@ -337,29 +318,26 @@ end @test gw .* hw' ≈ gl .* hl' rtol=1e-6 @test td.nsolves == tw.nsolves @test td.nsolves <= tl.nsolves - length((1e2, 1e4, 1e6, 1e8)) - # 22 solves measured here against 26 on `:lsqr`; 24 leaves two of headroom. + # Leave a small margin in the solve-count bound. @test td.nsolves <= 24 end -# The LSQR preconditioner absorbs the rows the continuation weights by κ, so the -# generalized spectrum it iterates on is the unweighted one and the iteration count -# stops growing as κ rises. CHOLMOD factors it, so it applies only in Float64; -# narrower and wider types run the plain matrix-free iteration. +# The Float64 LSQR preconditioner includes κ-weighted rows; other types use the +# plain matrix-free iteration. @testset "MMC :lsqr iteration count is bounded across the continuation" begin rng = StableRNG(5) A = (X = exp.(randn(rng, 120, 120)); (X .+ X') ./ 2) ad, _ = MatrixCovers._symcover_min_abslog2(A; linsolve=:dense) al, sl = MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr) @test al ≈ ad rtol=1e-6 - # 29.4 iterations per solve measured here; the bound doubles that, and an - # unpreconditioned run would sit in the hundreds by the last κ stage. + # Bound the preconditioned iteration count with margin. @test sl.lsqriters <= 60 * sl.nsolves G = exp.(randn(rng, 120, 90)) gd, hd, _ = MatrixCovers._cover_min_abslog2(G; linsolve=:dense) gl, hl, tl = MatrixCovers._cover_min_abslog2(G; linsolve=:lsqr) @test gl .* hl' ≈ gd .* hd' rtol=1e-6 - # 41.4 iterations per solve measured here, against the same bound. + # Apply the same iteration bound to asymmetric problems. @test tl.lsqriters <= 60 * tl.nsolves # A working type CHOLMOD cannot factor keeps the plain matrix-free iteration. @@ -371,11 +349,7 @@ end end @testset "MMC disconnected-support gauge" begin - # A support graph that splits into k connected components carries k independent - # (e; −e) gauges. The asymmetric dense normal equations pin only the global one - # with v0*v0ᵀ; a minimal scale-relative ridge lifts the remaining k−1, so - # `cover_min` no longer hits a SingularException on block-disconnected supports. - # The dense (`:auto`) and matrix-free (`:lsqr`) paths must agree. + # Dense and LSQR must handle independent gauges on disconnected support. singletons(vals) = Matrix(sparse(1:length(vals), 1:length(vals), float.(vals))) # k singleton components block2(k) = cat(([2.0+i i; i 3.0+i] for i in 1:k)...; dims = (1, 2)) # k dense 2×2 components for M in (singletons([4.0, 9.0, 1.0]), singletons(1.0:6.0), block2(3), block2(6)) @@ -386,7 +360,7 @@ end @test cover_objective(AbsLog{2}(), ad, bd, M) ≈ cover_objective(AbsLog{2}(), al, bl, M) rtol = 1e-6 atol = 1e-8 end - # The canonical failure mode: `cover_min` on a Diagonal (n singleton components). + # Diagonal input has one support component per entry. D = Diagonal([4.0, 9.0, 1.0]) a, b = cover_min(AbsLog{2}(), D) @test cover_objective(AbsLog{2}(), a, b, Matrix(D)) ≈ 0.0 atol = 1e-10 @@ -563,10 +537,7 @@ end @test median(gen_ratios) < 1.02 end -# The `:lsqr` path is the intended solve when nnz ≪ n²; its working set must -# therefore be sized by the support, not by `length(A)`. Doubling `n` at fixed -# nnz-per-row doubles the support, so allocation must roughly double too — a -# working set carrying any n×n array would quadruple instead. +# LSQR allocation should scale with support rather than matrix area. @testset ":lsqr allocates in proportion to the support" begin function lsqr_alloc(n) rng = StableRNG(4) @@ -576,10 +547,6 @@ end return @allocated MatrixCovers._symcover_min_abslog2(A; linsolve=:lsqr) end small, large = lsqr_alloc(200), lsqr_alloc(800) - # Measured ratio is ≈6: the support grows 4x, and the Newton continuation takes - # more LSQR iterations at the larger size. Quadratic growth would put it near 16. - # The bound discriminates sharply despite the gap, because one n×n Float64 array - # at n = 800 is 4.9 MB against a whole measured working set of ≈3.9 MB — a single - # reintroduced one lands the ratio above 13. + # Allow iteration growth while rejecting an added dense workspace. @test large < 10 * small end diff --git a/test/runtests.jl b/test/runtests.jl index 81efbad..3b03953 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -34,13 +34,7 @@ include("helpers.jl") # isbalanced, covaries, PENALTIES Aqua.test_all(MatrixCovers) @testset "ExplicitImports" begin - # The public-ness checks consult `Base.ispublic` only on 1.11+; before that they - # fall back to `isexported` and flag every `public`-but-unexported binding, so - # they are meaningful only on 1.11+. The other five checks run on every version. - # - # These are this package's own internals, which its extensions legitimately - # extend and call: extension and package ship from one repo at one version, so - # there is no cross-package promise to break. + # Public-name checks require Julia 1.11. Extensions may use package internals. internals = (:_cover_min_abslog2, :_symcover_min_abslog2, :_prepare_cover_start!, :_prepare_symcover_start!, :_prepare_soft_cover_start!, :_prepare_soft_symcover_start!, @@ -49,12 +43,7 @@ include("helpers.jl") # isbalanced, covaries, PENALTIES :require_abs_symmetric, :_edge_list, :_sym_edge_list, :_degrees, :_balance_cover!, :inflate_feasible!) - # Non-public names owned by other packages, each with no public equivalent: - # `FreeUnits`/`Unit` are Unitful's unit representation and `Units` their - # abstract supertype, needed to reject the unit types this package cannot - # read; `Optimizer` is the solver handle JuMP's own documented - # `Model(HiGHS.Optimizer)` entry point names; and `register_error_hint` is - # Base-internal. + # External non-public names with no usable public equivalent. foreign = (:FreeUnits, :Unit, :Units, :Optimizer, :Experimental, :register_error_hint) test_explicit_imports( MatrixCovers; diff --git a/test/soft_covers.jl b/test/soft_covers.jl index bf4839b..5663114 100644 --- a/test/soft_covers.jl +++ b/test/soft_covers.jl @@ -1,7 +1,6 @@ # Unconstrained AbsLinear soft covers: descent, multistart, and start provenance. -# Committed 5x5 matrix libraries (`symmetric_matrices`, `general_matrices`); the guard -# permits re-inclusion of this file in an already-initialized session. +# Reuse the committed 5×5 matrix libraries when already loaded. if !isdefined(@__MODULE__, :symmetric_matrices) include("testmatrices.jl") end @@ -54,9 +53,7 @@ end end end - # Continuity as a near-zero entry vanishes: AbsLinear has no discontinuity at r=0, so - # the soft cover varies continuously as A[2,2] → 0. The leave-one-out start drops the - # most-outlying small entry, reaching the same basin as the exact-zero case. + # `AbsLinear` remains continuous as an entry vanishes. γ = 0.5 A_zero = [γ 1.0; 1.0 0.0] A_small = [γ 1.0; 1.0 1e-10] @@ -66,14 +63,7 @@ end @test a_small ≈ a_zero atol=1e-5 end - # Covariance must survive the regime where the leave-one-out start wins: the entry - # dropped is selected by the scale-invariant log-residuals, so which basin wins - # cannot depend on the frame. (Weighting entries by raw |A[i,j]|² fails here: the - # entries' physical units differ, so their sums are incommensurate and a rescaling - # can flip the winning basin.) For A_small the residuals of the two diagonal entries - # tie exactly (true of every symmetric 2×2), where the tie-break uses raw magnitude; - # the scaling below preserves the magnitude ordering, as covariance under - # order-flipping scalings is unachievable on that degenerate class. + # The leave-one-out basin selection is covariant away from its documented tie. for ϕ in (AbsLinear{1}(), AbsLinear{2}()) for (B, d) in ((A_small, [50.0, 0.02]), ([3.0 7.6e-10; 7.6e-10 80.0], [35.0, 3400.0])) @@ -89,17 +79,13 @@ end end @testset "soft_cover" begin - # Closed form on Aε = [1 ε; ε 1]: the uniform-product critical point has - # a*b' ≡ (1+ε²)/(1+ε) on every entry. A single geometric-mean start converges to - # it. For small ε this is only a local minimizer — a strongly asymmetric solution - # that covers three entries and sacrifices one off-diagonal has lower objective — - # so the default multistart may (correctly) return a different, better product. + # Closed-form uniform-basin critical point for `Aε = [1 ε; ε 1]`. for ε in (0.5, 0.1, 1e-3) Aε = [1.0 ε; ε 1.0] target = (1 + ε^2) / (1 + ε) a, b = soft_cover(Aε; starts=1, maxiter=200) @test all(≈(target; atol=1e-10), a * b') - # The multistart never does worse than this uniform-basin local minimizer. + # Multistart cannot be worse than the included single start. am, bm = soft_cover(Aε; maxiter=200) @test cover_objective(AbsLinear{2}(), am, bm, Aε) <= cover_objective(AbsLinear{2}(), a, b, Aε) + 1e-12 @@ -201,9 +187,7 @@ end @test az[1] == 0 @test all(isfinite, az) && all(isfinite, bz) - # On a symmetric matrix this does not reduce to `soft_symcover`: freeing `a` from `b` - # relaxes the problem, so the two descents are minimizing over different sets and their - # fixed points differ. Only exact coverability forces them to agree. + # Free row and column factors need not match the symmetric solution. S = [4.0 1.0 2.0; 1.0 9.0 3.0; 2.0 3.0 16.0] S1 = (v = [2.0, 0.5, 3.0]; v * v') a, b = soft_cover(AbsLog{1}(), S1) @@ -234,8 +218,7 @@ end @test_throws "specify only one" soft_cover(A; σ=1.5, sigma=2.0) @test_throws "specify only one" soft_symcover(As; σ=1.5, sigma=2.0) - # best-of-8 never exceeds the single-start objective on the committed libraries - # (the multistart's incumbent is the single start, replaced only on improvement). + # Best-of-eight cannot exceed the included single-start objective. for (_, M) in general_matrices Mf = float.(M) a1, b1 = soft_cover(Mf; starts=1) @@ -259,15 +242,11 @@ end d = [2.0, 0.5, 3.0] @test covaries(soft_symcover, As, d; rtol=1e-7) - # The objective is the sharp covariant: it depends on the cover only through the - # scale-invariant ratios |A[i,j]|/(a[i]*b[j]), so it matches across frames to roundoff - # even where the cover itself does not (see "converged cover covariance" below). + # Objectives match across rescaled frames to roundoff. @test covaries_objective(AbsLinear{2}(), soft_cover, Ac, dr, dc; rtol=1e-12) @test covaries_objective(AbsLinear{2}(), soft_symcover, As, d; rtol=1e-12) - # On a hard lognormal-σ=5 ensemble the multistart strictly lowers the objective on a - # substantial fraction of a fixed corpus. Both the corpus and the solver's internal - # perturbation draws use `StableRNG`, so the count is fixed across Julia versions. + # Stable RNGs make the fixed-corpus multistart comparison reproducible. rng = StableRNG(2024) imp_sym = 0; imp_gen = 0 for k in 1:40 @@ -280,37 +259,26 @@ end cover_objective(AbsLinear{2}(), g8a, g8b, G) < cover_objective(AbsLinear{2}(), g1a, g1b, G) - 1e-9 && (imp_gen += 1) end - # Gate set a modest margin below the measured counts on this corpus (29 and 24 of 40): - # multistart must beat the single start on a solid fraction of instances. The counts - # depend on `maxiter`: the better each start converges, the less room a rival start has - # to improve on it, so raising `maxiter` lowers them. + # Require improvement on a substantial fraction of the corpus. @test imp_sym >= 24 @test imp_gen >= 19 end @testset "feasible start and provenance" begin - # The multistart fills caller-supplied `labels`/`objs` in place; the winner is - # `labels[_multistart_select(objs)]`, using the same selection rule as the solver. - # The positional arguments mirror `soft_symcover`'s `maxiter`, `starts` and `σ` defaults, - # so the instrumented call reproduces the public entry point exactly (asserted below). + # Instrumented multistart exposes candidate labels and objectives. function provenance(A; rng=StableRNG(0)) labels = String[]; objs = Float64[] a = MatrixCovers._soft_symcover_abslinear2(A, 32, 5, 2.0, rng; labels, objs) return a, labels[MatrixCovers._multistart_select(objs)], labels, objs end - # The greedy feasible cover (`init_feasible_diag!`) is offered as a start only when `A` has a - # zero entry. On this matrix every geometric-mean-derived start lands in one basin while - # `feasible` reaches a distinctly better one, so it is the selected winner. + # The feasible start reaches the better basin for this matrix with zeros. Afe = Float64[0 11 18 0 12; 11 0 1 0 20; 18 1 0 3 18; 0 0 3 18 0; 12 20 18 0 17] a, winner, labels, objs = provenance(Afe) @test winner == "feasible" # The instrumented call returns exactly what the public entry point selects. @test a == soft_symcover(Afe; rng=StableRNG(0)) - # `feasible` wins by a genuine basin gap, not descent-tolerance noise: it is - # co-optimal in the best basin (a perturbed start may also reach that basin and - # tie it to within the descent tolerance), and that basin beats every start in a - # different basin by a wide margin. + # The winning basin is separated beyond descent tolerance. fi = findfirst(==("feasible"), labels) @test objs[fi] <= minimum(objs) * (1 + 1e-6) other_basin = minimum(o for o in objs if o > objs[fi] * (1 + 1e-6)) @@ -337,32 +305,21 @@ end end @testset "converged cover covariance" begin - # A cover driven to convergence pins the objective to `eps` but its own entries only to - # `sqrt(eps)`. The objective is stationary at the minimizer, so a displacement `δ` along a - # direction of low curvature changes it by only `O(δ²)`; two frames of the same problem, - # whose entries differ by roundoff, therefore settle `O(sqrt(eps))` apart in `a` and `b` - # while agreeing on the objective to `O(eps)`. Most matrices have no such soft direction - # and co-vary to roundoff; this one does. + # A low-curvature direction permits O(sqrt(eps)) factor differences while + # objectives agree to O(eps). rng = StableRNG(1) B = exp.(2 .* randn(rng, 40, 40)) .* randn(rng, 40, 40) A = (B + B') / 2 dr = exp.(randn(rng, 40)); dc = exp.(randn(rng, 40)) @test covaries_objective(AbsLinear{2}(), soft_cover, A, dr, dc; rtol=1e-12) - # With the row/column gauge pinned, the two frames converge to the same cover and not - # merely to the same objective: the agreement is roundoff, far inside the `sqrt(eps)` a - # low-curvature direction would otherwise allow. Leaving the gauge free costs six orders - # of magnitude here, which is what makes the balance convention worth enforcing rather - # than merely documenting. + # Gauge balancing makes the rescaled factors agree to roundoff. @test covaries(soft_cover, A, dr, dc; rtol=1e-9) end @testset "soft AbsLog{2} is the exact unconstrained minimum" begin - # The soft AbsLog{2} objective is a linear least-squares in log space, so an - # oracle needs no solver: `pinv(M) * z` settles it directly. `M` always carries - # the (e; −e) gauge null direction in the asymmetric case and can be singular in - # the symmetric one (bipartite support), so compare objectives — which the gauge - # cannot move — rather than the scale vectors. + # Direct log-space least squares supplies the oracle; compare gauge-invariant + # objectives because `M` may be singular. function exact_sym(A) n = size(A, 1) S = [(i, j) for i in 1:n, j in 1:n if !iszero(A[i, j])] @@ -413,8 +370,7 @@ end @test isbalanced(a, b, A) end - # The geometric mean is the minimum only on a full support. Were the solvers to - # fall back to it, the sparse cases above would silently regress. + # Sparse support distinguishes the exact minimum from the geometric mean. @test cover_objective(AbsLog{2}(), initialize_symcover(sym_dense; strategy=:geomean, feasible=:none), sym_dense) ≈ cover_objective(AbsLog{2}(), soft_symcover_min(AbsLog{2}(), sym_dense), sym_dense) rtol=1e-8 @test cover_objective(AbsLog{2}(), initialize_symcover(sym_zeros; strategy=:geomean, feasible=:none), sym_zeros) > @@ -473,10 +429,7 @@ end @test soft_cover!(b1, c1, Agen) == soft_cover!(AbsLinear{2}(), b2, c2, Agen) end - # The refiner descends from the start it is handed; the multistart owns a menu. - # Under the convex AbsLog{2} the start is honored but cannot be seen in the result, - # while a non-convex AbsLinear{2} objective with two basins reports whichever the - # start lies in. + # Convex refiners ignore the basin; nonconvex refiners do not. @testset "start-dependence" begin Abasins = [0.021778451276962405 1.5690256886348526 1.5690256886348526 0.20473123461805692] @@ -541,9 +494,7 @@ end end end -# A matrix readable only through the support hook: any full-grid scan hits the -# throwing `getindex`. Both traversals are defined so the sym and asym kernels -# can be driven from the same fixture. +# Matrix readable only through support hooks; `getindex` throws. struct HookOnlyMatrix{T} <: AbstractMatrix{T} entries::Vector{Tuple{Int,Int,T}} # `i <= j` only; the transpose is implied n::Int @@ -565,8 +516,7 @@ function MatrixCovers.foreach_support(f, M::HookOnlyMatrix) end return nothing end -# Storing one member of each pair makes `abs`-symmetry structural, so the -# precondition check is a no-op — as it must be, since it cannot index `M`. +# One stored orientation makes magnitude symmetry structural. MatrixCovers.require_abs_symmetric(::HookOnlyMatrix, fname) = nothing @testset "the soft-cover kernels read through the support hook" begin @@ -578,9 +528,7 @@ MatrixCovers.require_abs_symmetric(::HookOnlyMatrix, fname) = nothing end start() = [1.3, 0.7, 2.1, 0.9] - # Reaching a result at all proves the kernel never indexed `M`; matching the - # dense run proves the gathered support is the same set of entries carrying the - # same full-grid multiplicity. + # Match the dense result without indexing `M`. for kernel! in (MatrixCovers._abslog1_iter!, MatrixCovers._abslinear1_iter!, MatrixCovers._abslinear2_iter!) @test kernel!(start(), M, 50) ≈ kernel!(start(), dense, 50) rtol=1e-10 @@ -594,8 +542,7 @@ MatrixCovers.require_abs_symmetric(::HookOnlyMatrix, fname) = nothing end end -# The kernels above are only half the path: the public entry points reach them -# through the initializers, so a full-grid scan in either one would surface here. +# Exercise support-only traversal through public initialization paths. @testset "the cover entry points read through the support hook" begin entries = [(1, 1, 2.0), (1, 3, 1.5), (2, 2, 3.0), (2, 4, 0.5), (3, 4, 4.0), (4, 4, 1.0)] M = HookOnlyMatrix(entries, 4) diff --git a/test/storage_types.jl b/test/storage_types.jl index e5b4231..62de122 100644 --- a/test/storage_types.jl +++ b/test/storage_types.jl @@ -112,15 +112,8 @@ end end @testset "traversal-based kernels match dense reference" begin - # unconstrained_min! and tighten_cover! are order-insensitive folds over - # foreach_support(_sym) (sum/min accumulations), so structured/sparse - # storage must agree with the dense form up to floating-point summation - # order (rtol=1e-12). The full symcover pipeline includes the bucketed - # feasibility boost, whose within-bucket processing order follows the - # storage type's traversal order: storages that traverse the canonical - # triangle in the dense fallback's column-major order are compared - # elementwise, the rest on feasibility and objective value. cover's boost - # order is likewise storage-dependent and is not compared elementwise. + # Order-insensitive kernels match dense results directly. For bucketed boosts, + # traversal-order differences are compared by feasibility and objective. rng = StableRNG(11) n = 8 Adense = randn(rng, n, n); Adense = Adense + Adense' @@ -142,11 +135,7 @@ end @test symcover(AbsLog{2}(), D) ≈ symcover(AbsLog{2}(), Matrix(D)) rtol=1e-12 @test symcover(AbsLog{2}(), St) ≈ symcover(AbsLog{2}(), Tsym) rtol=1e-12 # same symmetric-valued matrix - # Traversal order differs from the dense fallback's (Ssp_L keys pairs by - # the smaller index; St/Tsym visit all diagonal entries before any - # off-diagonal, dense interleaves them), so the bucketed boost's - # within-bucket order can differ: compare on feasibility and objective - # value rather than elementwise. + # Different traversal orders may choose different within-bucket ties. objclose(a1, M1, a2, M2) = isapprox(cover_objective(AbsLog{2}(), a1, M1), cover_objective(AbsLog{2}(), a2, M2); rtol=1e-2, atol=1e-10) for A in (Ssp_L,) @@ -198,11 +187,7 @@ end end @testset "native solvers on sparse and structured inputs" begin - # The native AbsLog{2} MMC solvers (`symcover_min`/`cover_min`) and the AbsLinear - # soft covers must agree with the dense reference on `Matrix(A)` when handed a - # sparse-backed or structured input, and the hard MMC covers must stay feasible. - # On a `SparseMatrixCSC`/`Symmetric`/`Hermitian`-sparse the MMC solvers default to - # the matrix-free LSQR inner solve; structured inputs use the generic dense path. + # Sparse and structured solvers match dense references; hard covers remain feasible. symdenses = [[2.0 1.0 0.0; 1.0 3.0 2.0; 0.0 2.0 5.0], [4.0 0.0 1.0; 0.0 0.0 0.0; 1.0 0.0 2.0]] # second has a zero row/column @@ -266,10 +251,7 @@ end end end -# `Symmetric`/`Hermitian` over a sparse parent store one triangle, so the -# asymmetric traversal must reconstitute the other. Without its own method these -# fall back to the generic full-grid `getindex` scan, which is correct but defeats -# the sparse specialization the wrapper exists to enable. +# Asymmetric traversal of sparse symmetric wrappers must emit both orientations. @testset "asymmetric traversal of wrapped sparse storage" begin P = sparse([1, 2, 1, 3], [1, 2, 3, 3], [2.0, 3.0, 1.5, 4.0], 3, 3) for W in (Symmetric(P, :U), Symmetric(sparse(transpose(P)), :L), diff --git a/test/support.jl b/test/support.jl index 1374a88..af2a12b 100644 --- a/test/support.jl +++ b/test/support.jl @@ -154,12 +154,7 @@ end @test symcover(Bidiagonal([3.0, 2.0, 1.0], [0.0, 0.0], :U)) isa AbstractVector end -# The kernels that gather the support into per-group neighbor lists read those -# lists in place of the matrix, so the gather must reproduce the matrix exactly — -# including multiplicity. `_sym_support` in particular enters each off-diagonal -# pair in both orientations, which is what gives a kernel accumulating over its -# groups the full-grid weighting of `cover_objective` (each off-diagonal pair -# twice, each diagonal entry once) rather than a halved one. +# Grouped support must preserve entries and symmetric full-grid multiplicity. @testset "grouped support reproduces the matrix" begin function regroup(S, groups_are_rows::Bool, sz) R = zeros(Float64, sz) @@ -198,10 +193,7 @@ end @test sort([(SO.idx[s], SO.val[s]) for s in MatrixCovers._slots(SO, -1)]) == [(-1, 2.0), (0, 1.0)] end -# The gauge freedom of an asymmetric cover (`a -> γ*a`, `b -> b/γ`) acts independently -# on each connected component of the bipartite support graph, so the balance -# convention that pins it must be imposed per component; these tests check the -# labeling `_support_components` builds directly. +# Check connected-component labels used for per-component gauge balancing. @testset "_support_components" begin rng = StableRNG(11) @@ -256,9 +248,7 @@ end @test ocolcomp == colcomp end - # The public `SupportComponents` wraps `_support_components` with accessors that - # take the matrix's own indices, so `gramcover` (and any other caller) can hold - # the structure and query it without re-traversing. + # Public component accessors use the matrix's own indices. @testset "SupportComponents accessors" begin B = randn(rng, 3, 2) C = randn(rng, 2, 4) diff --git a/test/unitful.jl b/test/unitful.jl index f5de585..3cde2c9 100644 --- a/test/unitful.jl +++ b/test/unitful.jl @@ -31,9 +31,7 @@ end @testset "coordinates on mixed scales" begin - # Coordinate 1 in m, coordinate 2 in mm. Julia promotes this to a common - # unit whenever the entries share a dimension, so an eltype that admits - # heterogeneous units is what preserves the scales as written. + # An abstract quantity eltype preserves heterogeneous written units. H = Quantity[4.0u"m^-2" 1.0u"m^-1*mm^-1"; 1.0u"m^-1*mm^-1" 4.0u"mm^-2"] a = symcover(H) @@ -51,9 +49,7 @@ @test unit.(b) == [u"kg^-1", u"K^-1", u"s/m"] @test iscover(a, b, B; rtol=8eps()) - # The unit gauge `a -> a*c`, `b -> b/c` is pinned by minimizing the total - # atomic-unit powers the two vectors carry, so a factor shared by every - # entry lands on whichever side has fewer of them. + # Minimize total atomic-unit powers to select the unit gauge. aJ, bJ = cover(B .* u"J") @test unit.(aJ) == [u"J/m", u"J/s"] @test unit.(bJ) == unit.(b) @@ -63,8 +59,7 @@ end @testset "cover reproduces symcover on symmetric input" begin - # The gauge takes its smallest-magnitude optimizer, which is trivial here: - # a symmetric matrix pins `unit(a[i])` outright via `a[i]^2 == A[i,i]`. + # Symmetric diagonal entries pin each scale unit. a, b = cover(A) @test unit.(a) == unit.(b) == UA @test a == b @@ -72,9 +67,7 @@ end @testset "uniform units" begin - # A concrete element type names one unit for every entry. The gauge's - # median interval is not a single point here -- `∑|t| + ∑|t+2|` is flat - # across `t ∈ [-2, 0]` -- so only its midpoint reproduces symcover. + # The median-interval midpoint makes `cover` agree with `symcover`. Uni = [4.0 1.0; 1.0 4.0] .* u"m^-2" @test isconcretetype(eltype(Uni)) @test unit.(symcover(Uni)) == [u"m^-1", u"m^-1"] @@ -109,9 +102,7 @@ end @testset "balance convention holds in the caller's units" begin - # `A` is stripped as written rather than in a canonical system, so the - # (non-scale-invariant) balance that splits `a` from `b` is pinned to the - # scale the caller named. + # Written units determine the balance convention's parametrization. ru, cu = (u"m", u"s"), (u"kg", u"K", u"m/s") B = [1.0 / (r * c) for r in ru, c in cu] .* [1e3 1e-2 5.0; 2.0 1e4 1e-1] a, b = cover(B) @@ -124,12 +115,11 @@ @test_throws "units of `A` do not factor" symcover(E) @test_throws "unit(A[2,2])*unit(A[1,1]) = kg^2 m^2" symcover(E) @test_throws "unit(A[2,1])^2 = s^2" symcover(E) - @test_throws "`A*x` is undefined for every `x`" symcover(E) + @test_throws "terms in a row of `A*x` can have incompatible units" symcover(E) @test_throws "units of `A` do not factor" cover(E) @test_throws DimensionMismatch cover(E) - # Dimensionally consistent, but the off-diagonal is named at a scale the - # diagonal contradicts: every entry is 𝐋^-2, yet no `a` covers them. + # Dimensions agree, but written unit scales do not factor symmetrically. M = Quantity[1.0u"m^-2" 1.0u"m^-2"; 1.0u"m^-2" 1.0u"mm^-2"] @test all(==(dimension(u"m^-2")), dimension.(M)) @test_throws "units of `A` do not factor" symcover(M) @@ -223,8 +213,7 @@ @test initialize_cover!(a, b, A) == (a, b) @test unit.(a) == unit.(b) == UA - # The `*_min!` family reads `a` as a start, so `initialize_*` feeds them - # directly: both sides name the units the same way. + # Initializer output is valid refiner input with matching units. a = initialize_symcover(A) @test symcover_min!(AbsLog{2}(), a, A) === a @test unit.(a) == UA @@ -272,10 +261,7 @@ end @testset "gramcover carries units" begin - # A rectangular matrix with a single uniform unit: the row/column split of - # `unit(A[i,j]) = unit(a[i])*unit(b[j])` need not be even (more rows than - # columns here pushes the whole unit onto `b`), but `gramcover` must carry - # whatever split `cover` returns through correctly regardless. + # `gramcover` must handle an uneven row/column unit split. J = [4.0 1.0; 1.0 3.0; 2.0 0.5] .* u"N/m" a, b = cover(J) s = gramcover(a, b, J) @@ -288,8 +274,7 @@ w = [1.0, 2.0, 0.5] .* u"s" sw = gramcover(a, b, J, w) @test unit.(sw) == unit(a[1]) .* unit.(b) .* unit(sqrt(1.0u"s")) - # Two-step product: the three-term `J' * Diagonal(w) * J` routes through a - # mixed-unit intermediate that throws a DimensionError on Julia 1.10. + # Form the weighted Gram product without a mixed-unit intermediate. Gw = J' * (Diagonal(w) * J) @test all(ustrip.(sw * sw') .>= ustrip.(abs.(Gw))) end