Run narrow-precision MMC solves in Float64

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

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

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