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*: Keir Mierle <mierle@xxxxxxxxx>*Date*: Tue, 25 Aug 2009 13:10:34 -0700*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; bh=9TZGQ/inm8jVFgycOgA4YiboUwqW300N4R7BIFTM4HQ=; b=OjkM6e6hSg1dyXE5zU2HBna9tbJqluLUvXVl27VMU7ZURI/aGGK3MYyi3FTT+ums8b 86+zyC27AvHRhGSpQvoDeWYX6glDqio6qohfrPQingFR9O57evay7+bpt6qEALOOmPeZ yjVgv8sVm17qfLu/BptfFLt3SsmDMxvsm2pfA=*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=XV5+hiHpsbIr75ZcrLuj4s2q3r6eiOjkqp00erSTgg6alCDAWYgcPczP5DQh1VFa7G EwLMitngliooyJXhTK6JMUA3HANzAoeVIwVxyAFDmEyQu4fQ70sTi6U5F3n9VZxp0sdT NWC5lPvoovA5958OjnC4MPfRPUVvwiHTl0fTg=

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

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

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

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] eigen short support: discussion thread** - Next by Date:
**Re: [eigen] eigen short support: discussion thread** - 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/ |