Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0 |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0*From*: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>*Date*: Mon, 19 Jul 2010 10:03: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 :content-transfer-encoding; bh=rkStwVFPLMXo4B7kxgwh8+5HGfZvAsc7OkzSBO1Oick=; b=wvNcChTWXmeoeFHI61NxUYGT2arVB0WParlvHDmULjmlDjSmkq6VkqLjxoKcJlPWIV yzaOxw/kxVDgaY+o3xL4V6rIxrrh1uJeXbT6gEdp1lSFFPpVFgnXZISintjjmn0J0/iQ NCRKOQ3ww16t7PshiiKo+Vr/WLwX9muZ2FkbY=*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=qSq3na4zkTh6YO8XU1zz/mpxSbngZOXHZk3dxSi/OLTiTJUKNwDWjz+rLgHPgolebZ MfoKbk1wSTSfXA5so1pSdELZsySqFH4+iEkwzMfObIeV5MKcdm/uAKX1REUa9pPxCHLW eaiVOb2NaHME0buF3wErTjUD05UirUodQMVVo=

2010/7/19 Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>: > On Sat, 17 Jul 2010, Benoit Jacob wrote: > >> If (lhs * solution - rhs) is small compared to rhs, then what more >> could you ask for? >> >> It might be that lhs is ill-conditioned but (perhaps just by chance) >> your solution is still good. > > Exactly. And in that case you want to know that the lhs is ill-conditioned, > because it means that rounding and other errors in the construction of the > lhs will have a big effect on the solution. OK, I understand. > > Perhaps a better reason is what Helmut said. One measure for the accuracy of > the solution x_computed, as computed by Eigen, is the distance to the true > solution, x_true. The relative error in the solution is bounded in terms of > the residual by > > | x_computed - x_true | / | x_true | <= \kappa(A) | r | / | b | > > where r = A x_computed - b is the residual and \kappa(A) is the condition > number, | A | * | A^{-1} |. Ah, that makes a lot of sense. I had understood part of that in my reply to Helmut but still couldn't see where the actual condition number came into play. Thanks! > >> ok, honestly I think I understand what you mean, so here is my more >> serious answer. I understand you mean that it is still useful to >> compute condition numbers and not just residuals, because in some >> situations you know that your condition number should never be higher >> than a certain value, so you want to check that to detect potential >> problems that could have happened earlier (leading to the formation of >> this ill-conditioned matrix). > > I hadn't actually thought of this. I'm not that such situations arise often > in practice. OK. > >> I am tempted to leave that for the user to do, because this is a >> really advanced use case and even in the context of a given >> decomposition, there are different inequivalent notions of condition >> numbers. For example in the case of LU, in this discussion we've >> already seen 3 inequivalent notions of condition number: >> - max/min absolute value of diagonal coeff of U >> - matrix.L1_norm() / inverse.L1_norm() (as proposed by LAPACK) >> - matrix.Linf_norm() / inverse.Linf_norm() (as proposed by LAPACK too) >> and this is only considering the condition numbers that can be >> computed directly from LU! Not even mentioning the "good" condition >> number max/min singval... > > For the second and third one, how is the user supposed to compute these? > LAPACK does not actually construct the inverse, because that costs O(n^3) > even if you have the LU decomposition. Indeed. > Instead, LAPACK does something > smarter which costs not much more than computing the residual and gives you > an estimate; the book of Golub and Van Loan has a description of the > algorithm. It is this *estimate* that I'd like to see in Eigen; it's a > nontrivial algorithm with (I believe) many use cases. I see. Sounds reasonable! Do I understand correctly that it would be OK to add as a new, separate method in existing decomposition classes, i.e. that it would not require changing any existing API? > I agree that it is not > very useful to have a one-line function which just computes A.someNorm() * > A.inverse().someNorm(), OK. What I didn't know was that there existed a nontrivial algorithm giving a good, quick estimate for that. >> Unfortunately, the different notions of condition numbers can give >> really different orders of magnitude, at least for large matrices. > > Yes, there are several things that can go wrong. One is the fact that the > norms are inequivalent, as you point out; it's tempting to conclude that the > error in the solution is small in the 2-norm while you compute the condition > number in the 1-norm. The second one is that the estimate for the condition > number may not be accurate. The third one is that the bound for the error in > the solution may not be tight. > > But this is how scientific computing usually works. You cannot prove that > your computation is correct, but you collect different pieces of data which > give you more and more confidence in the reliability of your computation. > One such piece is the residual, another one is the (estimate for the) > condition number. If they are low, it suggests you may be fine. If they are > high, it serves as a warning light. OK, I don't want to go against the current of what most users seem to be expecting. (Plus, all that sounds reasonable, even if personnally if I wanted something 100% reliable for a scientific computation, I would just use a fail-proof algorithm that is insensitive to the condition number, e.g. FullPivLU instead of PartialPivLU, and be done with it.) > Of course, this whole discussion is academic until somebody actually > implements the routine for estimating the condition number. Let me re-ask: do you confirm it's OK to add as new API without changing existing API? Then it can be done at any time in the future. Cheers, Benoit > > Cheers, > Jitse > > >

**Follow-Ups**:**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0***From:*Jitse Niesen

**References**:**[eigen] Problem inverting a Matrix4f with Eigen 2.0.0***From:*Robert Lupton the Good

**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0***From:*Christoph Hertzberg

**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0***From:*Gael Guennebaud

**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0***From:*Robert Lupton the Good

**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0***From:*Benoit Jacob

**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0***From:*Aron Ahmadia

**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0***From:*Benoit Jacob

**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0***From:*Benoit Jacob

**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0***From:*Aron Ahmadia

**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0***From:*Aron Ahmadia

**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0***From:*Benoit Jacob

**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0***From:*Jitse Niesen

**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0***From:*Benoit Jacob

**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0***From:*Jitse Niesen

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0** - Next by Date:
**Re: [eigen] Compiling Eigen2-based code in Eigen3** - Previous by thread:
**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0** - Next by thread:
**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0**

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