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:

https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/blob/40d9e9db7e0c6ee2047d121125005b217023baf9/src/simple.jl#L53 https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/blob/40d9e9db7e0c6ee2047d121125005b217023baf9/src/simple.jl#L119-L120

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:

https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/blob/40d9e9db7e0c6ee2047d121125005b217023baf9/src/lobpcg.jl#L751

---

(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.)