abstol/reltol for MINRES

Author: ranocha

src/minres.jl

@@ -32,16 +32,17 @@ mutable struct MINRESIterable{matT, solT, vecT <: DenseVector, smallVecT <: Dens
     # Bookkeeping
     mv_products::Int
     maxiter::Int
-    tolerance::realT
+    tol::realT
     resnorm::realT
 end
 
 function minres_iterable!(x, A, b;
-    initially_zero::Bool = false,
-    skew_hermitian::Bool = false,
-    tol = sqrt(eps(real(eltype(b)))),
-    maxiter = size(A, 2)
-)
+                         initially_zero::Bool = false,
+                         skew_hermitian::Bool = false,
+                         abstol::Real = zero(real(eltype(b))),
+                         reltol::Real = sqrt(eps(real(eltype(b)))),
+                         tol = nothing, # TODO: Deprecations introduced in v0.8
+                         maxiter = size(A, 2))
     T = eltype(x)
     HessenbergT = skew_hermitian ? T : real(T)
 
@@ -53,6 +54,12 @@ function minres_iterable!(x, A, b;
     w_curr = similar(x)
     w_next = similar(x)
 
+    # TODO: Deprecations introduced in v0.8
+    if tol !== nothing
+        Base.depwarn("The keyword argument `tol` is deprecated, use `reltol` instead.", :minres_iterable!)
+        reltol = tol
+    end
+
     mv_products = 0
 
     # For nonzero x's, we must do an MV for the initial residual vec
@@ -64,7 +71,7 @@ function minres_iterable!(x, A, b;
     end
 
     resnorm = norm(v_curr)
-    reltol = resnorm * tol
+    tolerance = max(reltol * resnorm, abstol)
 
     # Last active column of the Hessenberg matrix
     # and last two entries of the right-hand side
@@ -84,18 +91,21 @@ function minres_iterable!(x, A, b;
         w_prev, w_curr, w_next,
         H, rhs,
         c_prev, s_prev, c_curr, s_curr,
-        mv_products, maxiter, reltol, resnorm
+        mv_products, maxiter, tolerance, resnorm
     )
 end
 
-converged(m::MINRESIterable) = m.resnorm ≤ m.tolerance
+converged(m::MINRESIterable) = m.resnorm ≤ m.tol
 
 start(::MINRESIterable) = 1
 
 done(m::MINRESIterable, iteration::Int) = iteration > m.maxiter || converged(m)
 
 function iterate(m::MINRESIterable, iteration::Int=start(m))
-    if done(m, iteration) return nothing end
+    # Check for termination first
+    if done(m, iteration)
+        return nothing
+    end
 
     # v_next = A * v_curr - H[2] * v_prev
     mul!(m.v_next, m.A, m.v_curr)
@@ -172,7 +182,13 @@ Solve Ax = b for (skew-)Hermitian matrices A using MINRES.
   matrix-vector product can be saved when computing the initial
   residual vector;
 - `skew_hermitian::Bool = false`: if `true` assumes that `A` is skew-symmetric or skew-Hermitian;
-- `tol`: tolerance for stopping condition `|r_k| / |r_0| ≤ tol`. Note that the residual is computed only approximately;
+- `abstol::Real = zero(real(eltype(b)))`,
+  `reltol::Real = sqrt(eps(real(eltype(b))))`: absolute and relative
+  tolerance for the stopping condition
+  `|r_k| / |r_0| ≤ max(reltol * resnorm, abstol)`, where `r_k = A * x_k - b`
+  is the residual in the `k`th iteration
+  !!! note
+      The residual is computed only approximately.
 - `maxiter::Int = size(A, 2)`: maximum number of iterations;
 - `Pl`: left preconditioner;
 - `Pr`: right preconditioner;
@@ -191,20 +207,27 @@ Solve Ax = b for (skew-)Hermitian matrices A using MINRES.
 - `history`: convergence history.
 """
 function minres!(x, A, b;
-    skew_hermitian::Bool = false,
-    verbose::Bool = false,
-    log::Bool = false,
-    tol = sqrt(eps(real(eltype(b)))),
-    maxiter::Int = size(A, 2),
-    initially_zero::Bool = false
-)
+                 skew_hermitian::Bool = false,
+                 verbose::Bool = false,
+                 log::Bool = false,
+                 abstol::Real = zero(real(eltype(b))),
+                 reltol::Real = sqrt(eps(real(eltype(b)))),
+                 tol = nothing, # TODO: Deprecations introduced in v0.8
+                 maxiter::Int = size(A, 2),
+                 initially_zero::Bool = false)
     history = ConvergenceHistory(partial = !log)
     history[:tol] = tol
     log && reserve!(history, :resnorm, maxiter)
 
+    # TODO: Deprecations introduced in v0.8
+    if tol !== nothing
+        Base.depwarn("The keyword argument `tol` is deprecated, use `reltol` instead.", :minres!)
+        reltol = tol
+    end
+
     iterable = minres_iterable!(x, A, b;
         skew_hermitian = skew_hermitian,
-        tol = tol,
+        abstol = abstol, reltol = reltol,
         maxiter = maxiter,
         initially_zero = initially_zero
     )

test/minres.jl

@@ -29,25 +29,25 @@ n = 15
 
 @testset "Hermitian Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
     A, x, b = hermitian_problem(T, n)
-    tol = sqrt(eps(real(T)))
+    reltol = sqrt(eps(real(T)))
     x0 = rand(T, n)
 
-    x1, hist1 = minres(A, b, maxiter = 10n, tol = tol, log = true)
-    x2, hist2 = minres!(x0, A, b, maxiter = 10n, tol = tol, log = true)
+    x1, hist1 = minres(A, b, maxiter = 10n, reltol = reltol, log = true)
+    x2, hist2 = minres!(x0, A, b, maxiter = 10n, reltol = reltol, log = true)
 
     @test isa(hist1, ConvergenceHistory)
-    @test norm(b - A * x1) / norm(b) ≤ tol
+    @test norm(b - A * x1) / norm(b) ≤ reltol
     @test hist1.isconverged
-    @test norm(b - A * x2) / norm(b) ≤ tol
+    @test norm(b - A * x2) / norm(b) ≤ reltol
     @test x2 == x0
 end
 
 @testset "Skew-Hermitian Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
     A, x, b = skew_hermitian_problem(T, n)
-    tol = sqrt(eps(real(T)))
-    x_approx, hist = minres(A, b, skew_hermitian = true, maxiter = 10n, tol = tol, log = true)
+    reltol = sqrt(eps(real(T)))
+    x_approx, hist = minres(A, b, skew_hermitian = true, maxiter = 10n, reltol = reltol, log = true)
 
-    @test norm(b - A * x_approx) / norm(b) ≤ tol
+    @test norm(b - A * x_approx) / norm(b) ≤ reltol
     @test hist.isconverged
 end
 
@@ -59,12 +59,38 @@ end
 
     x = ones(T, n)
     b = A * x
-    tol = sqrt(eps(real(T)))
+    reltol = sqrt(eps(real(T)))
 
-    x_approx, hist = minres(A, b, maxiter = 10n, tol = tol, log = true)
+    x_approx, hist = minres(A, b, maxiter = 10n, reltol = reltol, log = true)
 
-    @test norm(b - A * x_approx) / norm(b) ≤ tol
+    @test norm(b - A * x_approx) / norm(b) ≤ reltol
     @test hist.isconverged
 end
 
+@testset "Termination criterion" begin
+    for T in (Float32, Float64, ComplexF32, ComplexF64)
+        A = T[ 2 -1  0
+              -1  2 -1
+               0 -1  2]
+        n = size(A, 2)
+        b = ones(T, n)
+        x0 = A \ b
+        perturbation = T[(-1)^i for i in 1:n]
+
+        # If the initial residual is small and a small relative tolerance is used,
+        # many iterations are necessary
+        x = x0 + sqrt(eps(real(T))) * perturbation
+        initial_residual = norm(A * x - b)
+        x, ch = minres!(x, A, b, log=true)
+        @test 2 ≤ niters(ch) ≤ n
+
+        # If the initial residual is small and a large absolute tolerance is used,
+        # no iterations are necessary
+        x = x0 + 10*sqrt(eps(real(T))) * perturbation
+        initial_residual = norm(A * x - b)
+        x, ch = minres!(x, A, b, abstol=2*initial_residual, reltol=zero(real(T)), log=true)
+        @test niters(ch) == 0
+    end
+end
+
 end