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

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] TriangularView::solve() interface*From*: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>*Date*: Tue, 25 Aug 2009 22:37:21 -0400*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :date:message-id:subject:from:to:content-type :content-transfer-encoding; bh=EgOIgrCi9PukOMaXOU2Y+xK+jqRhpLJOIYcb7bwCJEE=; b=mqABNRutYpF8soRizrzfHstbIN29F3ijgq/9YzdRfshZA2iU96TH0fmqpu8x02TEUq hbYZGYOHAM0hxz6OZWLxCLfnoeoRQpnoK92VM0nc1Rmd2CEOep685xVbw3j0xLWMRd9x TiX7QOO+7ztO8ISFLaSaOx/O3ScwXEu6cXoUY=*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:content-transfer-encoding; b=Bi0OEg2r9lEGeX77PMPkjqeF4CRiOTS9TyDr+HO+TPwcg2skq88SH52OfQv+jXIEiG trBY0xi1gwJeUDw3XRU+GSr0+sZ1Q25aHSBB5hldeOXVPgFeyF03h7800mxC/EoEaR7C EBWNjFweKtxeNd+C4dlZ1Opyt0pLe7HEMDujw=

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

**References**:**[eigen] TriangularView::solve() interface***From:*Jitse Niesen

**Re: [eigen] TriangularView::solve() interface***From:*Keir Mierle

**Re: [eigen] TriangularView::solve() interface***From:*Benoit Jacob

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] TriangularView::solve() interface** - Next by Date:
**Re: [eigen] TriangularView::solve() interface** - Previous by thread:
**Re: [eigen] TriangularView::solve() interface** - Next by thread:
**Re: [eigen] TriangularView::solve() interface**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |