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

Just reaching out to give an update on the progress. The tolerance plan
(with max() as the decision) is implemented in the PR (still in draft
state). I am fixing a few final bugs in GMRES but no blockers. I would
really appreciate you give it a go and engage in the PR. I know about these
solvers from the past but I don't use them in production since quite a long
time hence I am wondering if they are usable.

GMRES has a mind of its own due to many legacy behaviors piled up and I
tried to bring along as much as I can but it is a hot mess and I might have
introduced some bugs though existing tests pass.

Anyways, all feedback is most welcome.


On Wed, May 3, 2023 at 3:40 AM Charles R Harris <[email protected]>
wrote:

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

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