Re: Octave/C++ matrix Inv() comparison

[Pascal Dupuis]

Hello,
Ed Meyer already gave you a thorough explanation. With the condition
number of your matrix (the 7e-13 number you computed), the
computations of the smallest values of the inverse occur roughly with
2 significant digits.
 The differences between algorithms reflect their noise sensitivities.
It is NORMAL that different approaches give results which don't
compare.

To take an example from real life: you have to advance 10 m north. But
you decide to take another approach
- walk 10 km east, in a direction making an angle of 89.9425 degrees with North;
- walk 10 km west.
In theory, the distance in the north direction is
cos(89.9425/180*pi)*10000 = 10 m. In practice, what are your chances
of arriving where you wish to ?

The case is similar. With infinite accuracy, your matrix is
invertible. With double-precision arithmetic, your matrix is nearly
singular. Instead of using one eigenvector, you are using two
eigenvectors which are nearly aligned, and go forward in the direction
resulting from their arithmetic substraction.

In the first place, why do you need to invert the matrix ?

If you are trying to solve an overdetermined problem Ax = b, textbooks
says that the "best" (i.e. minimising the square of the sum of errors)
solution is to compute x=(A'*A)^-1 (A'*b). This is true in textbooks
and usefull to study the structure and properties of the solution.
>From a numerical point of view, you are SQUARING the condition number.
So a perfectly valid problem may become invalid, due to rounding
issues.

To get around: use the qr transform: A=Q*R, with Q orthonormal (Q' =
Q^-1) and R upper triangular.  The previous solution is equivalent to
solving R X = Q'b; as R is upper triangular you can use
backsubstitution. With this approach, there is no squaring of the
condition number.

If A is symmetric, it is even simpler: use the Cholesky factorisation:
A= R' R where R is upper triangular. In this case, you have to solve
two times by back-substitution a linear system of equation. No
squaring of the condition number.

If your matrix is REALLY badly conditioned on 64 bits arithmetic, the
only cure is to do it in 128 bits or 256 bits arithmetic.

Regards

Pascal

2013/7/10 ionone <address@hidden>:
> I computed the ratio of S as you told me and i end up with  this number :
> 1.9902e+000/2.6573e-014 = 7.49e13 so it's bigger than 1e8.
> The thing is that i have good datas with Inv(), i mean they make sense. For
> example with datas growing from 1/1024 to 1 in steps of 1/1024 i have an e
> (in the Octave function i posted before) coherent : from example, the
> numbers follow  pattern : 0.024 0.0234 0.0231 0.0224, etc which i don't have
> using Lapack dgetri/dsytri
>
> Jeff
>
> Le 09/07/2013 15:46, CdeMills [via Octave] a écrit :
>
> Hello,
>
> before starting to invest your time into comparing algorithms and toolboxes
> ... first look at your data. Perform a singular value decomposition, look at
> the ratio of the biggest singular value (top left of the S matrix) to the
> smallest one (bottom right). A ratio greater than 1/sqrt(eps), around 1e8,
> is calling for rounding issues. If this is the case, try understanding why
> this ratio is so low: it is the symptom that some of your columns are
> strongly linearly correlated. Common cures: compute pseudo-inverse (pinv),
> i.e. neglecting "directions" where noise sensitivity is great; or
> minimal-norm solutions.
>
> Regards
>
> Pascal
>
> ________________________________
> If you reply to this email, your message will be added to the discussion
> below:
> http://octave.1599824.n4.nabble.com/Octave-C-matrix-Inv-comparison-tp4655291p4655433.html
> To unsubscribe from Octave/C++ matrix Inv() comparison, click here.
> NAML
>
>
>
> ________________________________
> View this message in context: Re: Octave/C++ matrix Inv() comparison
> Sent from the Octave - Maintainers mailing list archive at Nabble.com.