Re: tol, atol and rtol behavior of scipy.sparse.linalg solvers

On Mon, May 1, 2023 at 8:17 PM Robert Kern <[email protected]> wrote:

> On Mon, May 1, 2023 at 7:22 PM Stefan van der Walt <[email protected]>
> wrote:
>
>> On Mon, May 1, 2023, at 11:49, Ilhan Polat wrote:
>>
>> If atol is actually set to something, then we return max(atol,
>> tol*norm(b)). Why max is used I don't understand.
>>
>> Typically, atol, rtol pair is used as tol = atol + rtol*<some metric>.
>> This is more or less what everybody expects from this pair and their naming
>> is chosen to reflect this. But I'm a bit lost in all the issues I could
>> read.
>>
>>
>> The first formulation — max(atol, rtol*norm(b)) — is what I'd expect.
>> I.e., you use atol when answers are small, or rtol for larger answers.  I'm
>> a bit surprised to see them summed; is this standard?
>>
>
> I concur. `max()` is the appropriate way to combine the two that I'd
> expect to see in almost all circumstances. `+` does have a history in our
> ecosystem as we do use `+` in `np.isclose()`, though, but I think that is a
> bit of a special case of comparing two numbers rather than the size of one
> number, and I've never been really convinced by the arguments for it.
> `math.isclose()` works like `max()` (though it is implemented with explicit
> `||` operations).
>
> Julia's iterative solvers use `max()`.
>
>
> https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/blob/master/src/cg.jl#L141
>
> MATLAB's numerical integration routines use `max()`.
>
>   https://www.mathworks.com/help/matlab/ref/integral.html#btc_m8o-4
>

That said, `+` is not at all uncommon outside of our own ecosystem either.
Julia's Krylov solvers use it:

  https://docs.juliahub.com/Krylov/0fcC3/0.5.5/solvers/

Here's some numerical analysis course that uses it:

  http://web.mit.edu/10.001/Web/Tips/Converge.htm

I've never really seen any *comparison* between the two outside of some
half-remembered arguments over `np.isclose()` early in development. I think
`max()` is the most justified and interpretable. It's just a way to code
the ORing of the two conditions. If either the absolute tolerance is met OR
the relative tolerance is met, then stop. With `+`, it's a weird mix of
both that's close to the `max()` result in most sane circumstances, but
diverges in some cases. When the user chooses really bad choices for the
two tolerances, neither will work particularly well, but `max()` is easier
to reason about why it doesn't work well in that case.

I strongly suspect `+` gets its lineage either from C codes or pre-F77
codes where there was no builtin `max()` to use.

--
Robert Kern

_______________________________________________
SciPy-Dev mailing list -- [email protected]
To unsubscribe send an email to [email protected]
https://mail.python.org/mailman3/lists/scipy-dev.python.org/
Member address: [email protected]