[ Thread Index |
| More lists.tuxfamily.org/eigen Archives
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: Re: [eigen] solve API
- From: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
- Date: Mon, 28 Jun 2010 08:35:18 -0400
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:received:in-reply-to :references:date:message-id:subject:from:to:content-type; bh=MFtoemdbE6uMZcf6r/wg5DAVsTgy4Uy3uxvVzmRBkOU=; b=DIMLr5+EPEemfSWOC63Jby5lAzOJewA+d0/o6Elayv1fPiyLWhV9Ie1FU05k0ZWeDR O2X35QEk85a6fthk3qYsJc15x2rFJuXuW3VJRFK5CqnuuWTOZse23uSRCpDWdUIhBfCa IiEcuntsXdqhHeK4NB25oP/0sqsYG/hhz8fPY=
- Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type; b=d/PIJUuF0AL3krTR2X2Zp3pBn66GgtbUf2g8TUgXusAgX3txSw3ZxXn37MfTDmxTHZ qV7PIRXaYic5ZrWZ+gQNTPDs6Y2kadwkBjYh0ksxgNcPQ6V/JZ6n6engBce6Bcxb1l7d 7WlY0008J+gIxXjF5b/qAxqta4FCuUAjsHiqY=
2010/6/28 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
> I was thinking that we did not discussed the solved API for a while,
> at so it would be perfect time to discuss it again ;)
OK to add such new API, but let's keep the existing solve() API
- as you mention below, computer programmers will be expecting such
an imperative API;
- we're about to release beta1.
> Recall that one goal is to be able to write code as you write math. So
> here is the proposal:
> x = m.lu().inverse() * b;
> x = b * m.lu().inverse();
> x = m.lu().adjoint().inverse() * b;
> x = b * m.lu().adjoint().inverse();
I would much prefer:
x = m.lu().inverse().adjoint() * b;
x = b * m.lu().inverse().adjoint();
Other than that, I completely agree with your proposal (as long as the
existing solve() API remains available) although I would like to
discuss the name of this method (see below).
Your proposal is Eigen-ish because it is object-oriented. The API only
returns an object, doesn't modify existing objects. I think the most
object-oriented API is one where all methods are const-qualified
except for operator=.
> .inverse() and .adjoint() would return simple proxy with operator *
> overloaded to return a "solve" proxy that we already have (so
> implementation-wise that's not much).
> I know that computer scientists used to other linear algebra packages
> might be shocked by such an API and will look forever for a function
> having "solve" in its name, but this is pure Eigen style IMO and it is
> consistent with the rest of the API (diagonal, permutation,
> transposition already use this API).
Yes, keeping the existing imperative solve() API is very important too.
The two can coexist, the only potential conflict that I see is in LU
classes that already have a inverse() method, we have to make sure
(say for beta2) that this doesn't break existing behavior (it's ok to
break some corner cases until 3.0 is released).
> It remains the problem of convenient "inplace" solving. Currently, when doing:
What do you think about:
as we already do in Householder?
Finally, let's discuss another issue: for decompositions that allow
non-invertible matrices, such as FullPivLU, calling this method
inverse() is a bit problematic. So here's a proposal. Since the
imperative API is:
how about calling this method solve() so that the new API is really
very consistent with the old one:
m.lu().solve() * x;
Now you might object that since this method just returns an object, it
should be named not as a verb, but as a noun. So how about solver() ?
m.lu().solver() * x;
Either way (calling it solve() or solver()) we also solve the problem
of the potential conflict with existing inverse() methods (yes I know
these inverse() methods already return proxies, but those
ReturnByValue proxies offer much more API than we have to here, since
they are real expressions).
> v = m.lu().solve(v);
> we can detect that we have the same object for the "x" and the "b" via
> one runtime pointer comparison, and then skip the matrix copy of b to
> x before the inplace solving path. So the same mechanism would still
> b = b * m.lu().adjoint().inverse();
> => inplace solving, no temporary, no matrix copy (except for the
> decomposition itself but that's another story).
> The pb, is when "b" is a complex expression:
> b.block(i+s, j+block_size-1, rows-r, std::min(blocksize,cols))
> =m.lu().inverse() * b.block(i+s, j+block_size-1, rows-r,
> You don't want to write it twice!!
> Here a solveInPlace() function could make sense:
> m.lu().solveInPlace(b.block(i+s, j+block_size-1, rows-r,
> better ? not so sure:
> 1) what about "ontheright" and transposed variants ?
> 2) and even worse, this does not work because block(...) returns a
> temporary object which has to be passed by const reference, so you
> have to name it anyway:
> Block<...> b_block(n,i+s, j+block_size-1, rows-r, std::min(blocksize,cols));
> and then repeating twice b_block is fine.
> As some of you probably guessed, I did not slept well this night....