Re: [eigen] TriangularView::solve() interface

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

OK, here goes.

In a solve() method, the error measurement is really a real number:

(rhs - matrix * returned_solution).norm()  /  rhs.norm()

The problem of our API where solve() returns a bool, is that this real
number is mapped to a bool; and if the bool is false, the computation
is stopped and the result is left undefined. Even if we agree on a
good threshold, and even if we do that check rigorously (which is much
more costly than what we do at the moment) there remains a fundamental
flaw here: if the precision is just a little coarser than the
threshold, the user is probably still interested in the result, rather
than having no result at all.

The argument in favor of exiting the computation and leaving the
result undefined was the fear of overflow/nan/other bad things, but
actually if you look at the implementation, that risk has nothing to
do with the existence of solutions (there are things to improve in
this respect but it's unrelated).

So my proposal is to just drop all the bool return values in solve and
not check at all for the existence of solutions.

If the user wants to check existence, he can compute

(rhs - matrix * returned_solution).norm()  /  rhs.norm()

and compare to his own favorite threshold.

I'll propose a patch shortly...

Then combined with Jitse's email and my previous email this could give
a complete uniformization of the solve() methods:

result = matrix.solve(rhs);

where we use ReturnByValue to make sure that the RVO happens.

ok?
Benoit

2009/8/25 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
> As I'm going to explain ASAP in a mail to this list, I believe that we
> can drop altogether the "bool success" in all solve() methods.
>
> There remains the choice of an API between
> API 1: a.solve(b,&result)
> and
> API 2: result = a.solve(b)
>
> The main argument in favor of API 1 is the one Keir mentions.
>
> But what about the wonderful little ReturnByValue class that Gael
> wrote? I suggest we investigate that, at the moment I don't see why it
> wouldn't work here! It would allow to have API 2 with the guarantee of
> it unfolding to the same as API 1.
>
> Maybe Gael can comment.
> Meanwhile, is there consensus that, if it works, then we want it?
> Again, this is assuming that the bool return value are removed, i'm
> going to write to the list to advocate that.
>
> Benoit
>
>
> 2009/8/25 Keir Mierle <mierle@xxxxxxxxx>:
>> I had some thoughts about this, so before I forget:
>>
>> On Mon, Aug 24, 2009 at 7:36 AM, Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>
>> wrote:
>>>
>>> Hello,
>>>
>>> While rewriting part of the tutorial, I noticed that the new
>>> TriangularView::solve() function returns the solution vector, while the
>>> other solve functions pass the solution via one of the parameters. Compare
>>>
>>>   x = A.triangularView<UpperTriangular>().solve(b);
>>
>> With this interface, are we guaranteed that NVRO will work? In other words,
>> are we sure this won't force memory copies and reallocations? With the
>> pointer interface, I am confident that eigen will do the right thing if my
>> matrix is the right size. With the return value, I have less control. In
>> particular, I believe you are only guaranteed that the following:
>> Matrix<...> x = A.lu().solve(b)
>> will directly construct x instead of via a temporary. Consider:
>> Matrix<double, HUGE, HUGE> x = A.lu().solve(b);  // Ok, NVRO works here
>> // ...
>> // ...
>> // form new b
>> x = A.lu().solve(b)  // Oops! x is already constructed, so here we must
>> create a temporary and copy.
>> Many numerical problems involve solving systems with the same dimensions
>> repeatedly; it would be a shame if there was no ability to avoid the price
>> of large temporaries.
>>>
>>> with
>>>
>>>   A.partialLu().solve(b, &x);
>>
>> An alternative, which I am not necessarily advocating because of the reasons
>> above, is the following:
>> bool success;
>> x = A.partialLU().solve(b, &success);
>> This has the benefit that success can be an optional parameter, and the API
>> becomes uniform.
>> Yet another thought:
>> If we make the solver interfaces like the following:
>> bool maybe_success = A.lu().solve(b, &x);
>> and changed the meaning of the return value to the following:
>> If the return value is true, then the solve *may* have succeeded. If the
>> return value is false, then the system *definitely* did not solve. This way
>> we can use the same API everywhere, and ignore the return value in cases
>> where it's known that the algorithm always returns true. Each solve()
>> implementation will declare what behaviour it has regarding success, with
>> the default behaviour to always return true.
>> Thoughts?
>> Keir
>>
>>> Was this done on purpose?
>>>
>>> The docs on the eigen website [1] were magically updated, but some of the
>>> latex was not parsed correctly. Am I doing anything wrong? Everything works
>>> fine on my own computer.
>>>
>>> [1]
>>> http://eigen.tuxfamily.org/dox-devel/TutorialAdvancedLinearAlgebra.html
>>>
>>> Something to perhaps add to the to-do list: a routine to estimate the
>>> condition number of a matrix based on the LU decomposition.
>>>
>>> Cheers,
>>> Jitse
>>>
>>>
>>
>>
>