Spurious scaling of the default tolerance in powm

[danielwe]

From the docs for powm!:

tol::Real = eps(real(eltype(B))) * size(B, 2) ^ 3: stopping tolerance for the residual norm

that is, the tolerance scales as the cube of the dimension of the vector space. This is not a documentation typo, it's what the actual code does:

IterativeSolvers.jl/src/simple.jl

Line 53 in 40d9e9d

function powm_iterable!(A, x; tol = eps(real(eltype(A))) * size(A, 2) ^ 3, maxiter = size(A, 1))

IterativeSolvers.jl/src/simple.jl

Lines 119 to 120 in 40d9e9d

	function powm!(B, x;
	tol = eps(real(eltype(B))) * size(B, 2) ^ 3,

As usual, the residual norm is ‖Ax - λx‖₂, and x is normalized at every iteration (see the iterate function starting on line 28 in the same file).

This seems wrong. With normalized x, there's no reason the tolerance should scale with the dimension in any way, is there? If there is a logic behind this I'd love to learn.

For comparison, the other eigensolver in this package, lobpcg, does no such scaling, but sets the default tolerance to eps^0.3:

IterativeSolvers.jl/src/lobpcg.jl

Line 751 in 40d9e9d

default_tolerance(::Type{T}) where {T<:Number} = eps(real(T))^(real(T)(3)/10)

---

(I'd argue that the better way to assess convergence of an eigenvalue-eigenvector pair is through a scale-invariant tolerance, the way Arpack.jl and ArnoldiMethod.jl do, by using the criterion ‖Ax - λx‖₂ < tol * |λ| (modulo some abstol for when λ ~ 0). That way you effectively impose a relative tolerance on the eigenvalue and an absolute tolerance on the direction of the eigenvector, and the convergence of an eigenpair is invariant to rescaling and inversion of A. This would be a nice improvement for both powm and lobpcg, but is orthogonal to the spurious scaling issue with powm.)