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*: Sat, 17 Jul 2010 18:32:26 -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=UvjrZiERqx95Z7zrXx8yRo7EKO3iaqxoAutAdVeXn08=; b=JxYzKh6WdFFspcI/ekU7kt9baCnIGn9WuNNq1CgyxN2kBvotaC5jI41jyXwNrdadlN oLnYfCNxpc9CCaR4PSopx4IFRbI7UsxfXJeQwLw179d4hctu+HvSP876VREEOc0KORSJ cyixBcni9c3iY59XAJWn0EXwyadihQsbpsC+A=*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=vUTn81yUFec46otNDcBRXkraQ8Z3JuJHpYi4ydv5vKRulNDQ+X+QwkV5hkLelFTXdm C24TC1BAEeJAGbSOk7/7ldZZ01m6OY+GsaWbh9jxhdmtXVLb9hcMeu29ieRdnG8+xNJP F5o0hfEPcpoSfVOkthFdW/DmK7N/1UdjSf9XA=

2010/7/17 Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>: > On Sat, 17 Jul 2010, Benoit Jacob wrote: > >> Another reason why we decided not to expose condition numbers in Eigen >> is that for the main purpose they're used for, namely checking if a >> result is reliable, there is a better approach which is: check the >> result itself. For example, if you want to check how accurate your >> matrix inverse is, just compute matrix*inverse and see how close it is >> to the identity matrix. Nothing beats that! When it comes to more >> general solving with potentially non full rank matrices, this is even >> better, because the condition number of the lhs matrix alone doesn't >> tell all you need to know (it also depends on your particular rhs), so >> the approach we're recommending in Eigen, to compute lhs*solution and >> compare with rhs, is the only way to know for sure how good your >> solution is. > > I'm not sure I agree. The residual (lhs * solution - rhs) is certainly > useful, but so is the condition number (norm of A multiplied by norm of > A^{-1}, where I deliberately don't specify the norm used). If the residual > is zero but the condition number is huge then you probably should not trust > the solution, 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. > because even though Eigen may have found the exact solution to > the problem you gave it, it is likely that the lhs and rhs that you gave to > Eigen are not quite exact, In the sense of being not quite what the user had in mind? We can't catch PEBKAC mistakes :-) 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 just want to make it clear that trying to check this, is a much more advanced use case than just trying to check if the solution is correct. The residual check does guarantee that the solution is correct for the given lhs and rhs. You're trying to do something more involved than that: when the user knows that a certain condition number shouldn't exceed a certain value, you want to enable checking that. 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... Since the use case for that is that the user knows that a *certain* condition number shouldn't exceed a certain threshold, and since a mismatched notion of condition number would amount to changing that threshold, I think it's best to let the user do that himself. > This implies (if correct) that it is useful to have Eigen compute estimates > for the condition number. The operative word here is 'estimate': I think the > order of magnitude is what is usually important, not the exact value. Unfortunately, the different notions of condition numbers can give really different orders of magnitude, at least for large matrices. Let's check L1 and Linf norms since that's what LAPACK is doing. I don't know if they mean entrywise or induced norms so let's give examples in both cases. First let's consider entrywise norms: For example, consider a NxN matrix all of whose coefficients are equal to 1. Its Linf entrywise norm is 1. Its L1 entrywise norm is N^2. So if you have N=1000, the ratio between the two is 1e+6 ! Now let's consider induced norms. (http://en.wikipedia.org/wiki/Matrix_norm) Consider a NxN matrix whose first column is filled with 1's and the rest with 0's. Its L1 induced norm is N Its Linf induced norm is 1 So now if your have N=1000, the ratio between the two is still 1000. Not as bad, but still bad! In conclusion, inequivalent notions of condition numbers are really inequivalent, even if you just need an order of magnitude. Benoit > > Cheers, > Jitse > > > >

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

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

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0** - Next by Date:
**Re: [eigen] Behavior of sum() and prod() for empty vectors** - 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/ |