diff --git a/src/minimal_covers.jl b/src/minimal_covers.jl index e9af960..5bb3408 100644 --- a/src/minimal_covers.jl +++ b/src/minimal_covers.jl @@ -47,6 +47,10 @@ 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. +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. + !!! note Even the native solver is more expensive than the [`symcover`](@ref) heuristic. @@ -100,6 +104,10 @@ 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. +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. + !!! note Even the native solver is more expensive than the [`cover`](@ref) heuristic. @@ -523,6 +531,8 @@ end # 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. @@ -540,6 +550,20 @@ function _symcover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # 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`. 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`. + if eps(T) > eps(Float64) + a64, stats = _symcover_min_abslog2(convert(AbstractMatrix{promote_type(eltype(A), Float64)}, A); + κs, maxiter, linsolve, start, boost, fname) + return T.(a64), stats + end n = length(ax) use_lsqr = linsolve === :lsqr # CHOLMOD, which factors the LSQR preconditioner, is reliable only in Float64; @@ -982,7 +1006,7 @@ end # Worker for `cover_min(::AbsLog{2})`. Returns `(a, b, stats)` with `stats` a # NamedTuple `(; nsolves, lsqriters, cgiters, cholsolves, linsolve)` (see -# `_symcover_min_abslog2`). +# `_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. function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), @@ -996,6 +1020,20 @@ function _cover_min_abslog2(A::AbstractMatrix; κs=(1e2, 1e4, 1e6, 1e8), # 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`. 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`. + if eps(T) > eps(Float64) + a64, b64, stats = _cover_min_abslog2(convert(AbstractMatrix{promote_type(eltype(A), Float64)}, A); + κs, maxiter, linsolve, start, boost) + return T.(a64), T.(b64), stats + end m = length(axr) n = length(axc) N = m + n diff --git a/test/element_types.jl b/test/element_types.jl index c507850..31b4f12 100644 --- a/test/element_types.jl +++ b/test/element_types.jl @@ -22,6 +22,67 @@ @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. + @testset "narrow working types solve in Float64" begin + rng = StableRNG(77) + X = exp.(randn(rng, 60, 60)) + Asym = (X .+ X') ./ 2 + Agen = exp.(randn(rng, 60, 45)) + aref = symcover_min(AbsLog{2}(), Asym) + sref = soft_symcover_min(AbsLog{2}(), Asym) + gref, href = cover_min(AbsLog{2}(), Agen) + A32 = Float32.(Asym) + + # Every storage the solvers specialize on reaches the same cover. + @testset "$name" for (name, M) in ("Matrix" => A32, + "Symmetric" => Symmetric(A32), + "SparseMatrixCSC" => sparse(A32), + "Symmetric{SparseMatrixCSC}" => Symmetric(sparse(triu(A32))), + "Hermitian{ComplexF32}" => Hermitian(ComplexF32.(A32))) + a = symcover_min(AbsLog{2}(), M) + @test a isa Vector{Float32} + @test a ≈ aref rtol=1e-5 + end + + # Offset axes survive the promotion and the conversion back. + Ao = OffsetArray(A32, -1, -1) + ao = symcover_min(AbsLog{2}(), Ao) + @test ao isa OffsetVector{Float32} + @test axes(ao, 1) == axes(Ao, 1) + @test collect(ao) ≈ aref rtol=1e-5 + + # `Float16` pins the returned cover only to its own precision, which is what + # the looser tolerance measures; the solve behind it is the same `Float64` one. + a16 = symcover_min(AbsLog{2}(), Float16.(Asym)) + @test a16 isa Vector{Float16} + @test a16 ≈ aref rtol=2e-3 + + # The soft cover is the same worker with no continuation and no boost. + s32 = soft_symcover_min(AbsLog{2}(), A32) + @test s32 isa Vector{Float32} + @test s32 ≈ sref rtol=1e-5 + + # Asymmetric: the product is the gauge-invariant object to compare. + G32 = Float32.(Agen) + g32, h32 = cover_min(AbsLog{2}(), G32) + @test g32 isa Vector{Float32} && h32 isa Vector{Float32} + @test g32 .* h32' ≈ gref .* href' rtol=1e-5 + gs, hs = cover_min(AbsLog{2}(), sparse(G32)) + @test gs .* hs' ≈ gref .* href' rtol=1e-5 + + # A type at least as wide as `Float64` is solved in itself. + Abig = BigFloat.(Asym[1:8, 1:8]) + abig = symcover_min(AbsLog{2}(), Abig) + @test abig isa Vector{BigFloat} + @test Float64.(abig) ≈ symcover_min(AbsLog{2}(), Float64.(Abig)) rtol=1e-6 + @test abig != BigFloat.(Float32.(abig)) + end + @testset "BigFloat flows through the family" begin A = BigFloat[4 1.5; 1.5 1] for a in (symcover(A), soft_symcover(A), symcover_min(AbsLog{2}(), A), diff --git a/test/minimal_covers.jl b/test/minimal_covers.jl index 7f9eb99..eb239cf 100644 --- a/test/minimal_covers.jl +++ b/test/minimal_covers.jl @@ -281,12 +281,20 @@ end @test_throws "at most 4·max(m, n) = 160 zeros in total" cover_min(AbsLog{2}(), gwide; linsolve=:woodbury) @test MatrixCovers._cover_min_abslog2(gwide)[3].linsolve === :dense + # A working type narrower than Float64 is promoted, so it reaches the Woodbury + # path; a wider one keeps its precision and is refused by the CHOLMOD-backed solve. A32 = Float32.(symlognormal(8)) - @test_throws "requires Float64 arithmetic" symcover_min(AbsLog{2}(), A32; linsolve=:woodbury) - @test MatrixCovers._symcover_min_abslog2(A32)[2].linsolve === :dense + a32, s32 = MatrixCovers._symcover_min_abslog2(A32) + @test s32.linsolve === :woodbury && eltype(a32) === Float32 + @test symcover_min(AbsLog{2}(), A32; linsolve=:woodbury) ≈ symcover_min(AbsLog{2}(), A32; linsolve=:dense) rtol=1e-6 + Abig = BigFloat.(symlognormal(8)) + @test_throws "requires Float64 arithmetic" symcover_min(AbsLog{2}(), Abig; linsolve=:woodbury) + @test MatrixCovers._symcover_min_abslog2(Abig)[2].linsolve === :dense G32 = Float32.(lognormal(8, 6)) - @test_throws "requires Float64 arithmetic" cover_min(AbsLog{2}(), G32; linsolve=:woodbury) - @test MatrixCovers._cover_min_abslog2(G32)[3].linsolve === :dense + @test MatrixCovers._cover_min_abslog2(G32)[3].linsolve === :woodbury + Gbig = BigFloat.(lognormal(8, 6)) + @test_throws "requires Float64 arithmetic" cover_min(AbsLog{2}(), Gbig; linsolve=:woodbury) + @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