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

On Tue, May 2, 2023 at 1:06 PM Robert Kern <[email protected]> wrote:

> On Mon, May 1, 2023 at 9:04 PM Charles R Harris <[email protected]>
> wrote:
>
>>
>> On Mon, May 1, 2023 at 5:21 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?
>>>
>>
>> IIRC, I first saw it in *Numerical Recipes* (1986).
>>
>
> To be fair, I *think* it ultimately derives from a principled combination
> of sources of numerical error. If your inputs have absolute roundoff error
> (i.e. rounded to the nearest 0.0001) and the floating point computation
> contributes `rtol` amount of relative roundoff error during the computation
> (usually something like `n_flops * eps`), then the total error should be
> about the sum of the two.
>
> In general, I still lean towards `max()`. The principled application of
> the sum depends on knowing the details of the computation and the data
> sources to properly scale each of the tolerances. I'm _pretty_ okay with
> understanding the details of these things, but in reality, I'll just stick
> with the defaults and tweak up or down when it doesn't work. And in that
> scenario, `max()` is more interpretable to me.
>
>
I agree with this. In practice, the choice doesn't matter much. The sum
version has the advantage of being smoother, but I suspect it was just an
easy way to get something useful. Whether one or the other is easier to
reason about depends on how much time you want to spend reasoning about it
:) I will note that a lot of the early numerical stuff was either in Algol
68 (died) or Fortran, C wasn't a suitable language, not least because it
only had doubles at a time that single precision floats were about 4x
faster, at least on the popular VAX series, and the VAX file system was
record based (Fortran) rather than byte stream based. Of course, there were
a lot of different float versions before the IEEE standard came out.

Chuck

_______________________________________________
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]