Re: [eigen] solve API

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

2010/6/28 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
> On Mon, Jun 28, 2010 at 5:15 PM, FMDSPAM <fmdspam@xxxxxxxxx> wrote:
>>  Just my 2 cent:
>>
>> I know, that expression do not have to be evaluated but nevertheless
>> I would try to avoid the word "inverse()", because of lession #1: "Never use
>> matrix inversion".
>
> This rule exists only because of the poor implementation/API of the
> other libs ;)
>
> I'm kidding, I see your point, but I don't see it as a strong argument
> against inverse.
>
> In books you never use "solve" in your mathematical formula, but
> always write, e.g., A^1 B. Wouldn't it be convenient to be able to
> write code which closely follow what you found in a textbook, or to
> write code as you write formulas on your paper sheets without worrying
> about performance?

In math, one indeed says A^-1 B when A is invertible. Or at least when
A is injective (so it makes sense as the unique element of A^-1({B})
). But I've never seen A^-1 B used when A isn't injective. In that
case, one says in plain text "let X be a solution of AX=B". There is
no way to write that as a formula.

So I see your point in the invertible case, but I feel that there is a
problem in the non-invertible case. Since invertibility is a numerical
issue, I agree not to let it impact the API. But at least, non-square
matrices should be banned.

How about this proposal:
 - call it solver() in the general case
 - add inverse() as another method that checks that the matrix is
square, and returns solver().

Of course, in all decompositions that are restricted to square
matrices, no need for solver(), so this amounts to just your proposal:
inverse().

Benoit

> Actually, I'm pretty sure that many newcomers to
> linear algebra do not know about these details, and simply write what
> they read in their favorite book:  A.inverse() * B (I've seen that
> several times). So if A.inverse() would return a A.lu().inverse()
> proxy, then no more worry.
>
>>
>> What about
>>
>> x = b * A.adjoint().lu().invertor()     (like benoit's solver())
>>
>
> As I said to Manoj, this is not equivalent to have the adjoint applied
> on the inverse proxy (lu.inverse().adjoint()) because I want to be
> able to use the LU dec to solve both adjoint and non adjoint problems.
>
> So it should be:
>
> LU<> lu(A);
>
> x = lu.invertor() * b;
> y = lu.invertor().adjoint() * b;
>
> but "invertor().adjoint()" is a bit weird. Nonetheless, I prefer
> invertor()  to solve() or solver().
>
>
>>
>> But cleaner for me to read is, what Manoj suggesting:
>>>
>>>     x = solve(A,b); // left
>>>     x = solve(b,A); // right
>>>     x = solve(b, A.triangularView<Lower>()); // Gael's example
>>
>> It's OK to read and understand w/o scarring about "inverse()" too.
>
> it is really ambiguous because
>
> x = solve(b,A);
>
> can mean x = b^-1 A, and anyway this proposal works at a different
> level that the one we are considering.
>
> cheers,
> gael
>
>>
>> --Frank
>>
>>
>>
>
>
>