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

[ionone]

thank you for your thorough answer

What i'm trying to do is to implement a matlab algorithm in C++

starting from there, I just copied the matlab functions in C++. It's a DSP algorithm. And with all the inverse functions i used, i have altered outputs.

for exemple, if in the input i have the matrix : R = [1/1024;1/1024;1], i expect, after several calculations to have something also growing slowly like : 0.021 0.025 0.030, etc

With octave's inverse function everyting is accurate, but as soon as i use another inverse function (and i've tried a lot) the numbers breaks : for example : 0.021 0.025 -6e-12. And i checked every other part of the algorithm and i finally i spotted the inv() function so i'm sure it comes from there.

If i use other inv() function that might be close i've got an output so high it mutes the track i have music-to-be-treated on.

I don't know exactly yet why there is this inverse() function in the algorithm because i'm tring to figure how it works. And as it said in the comment of this algorithm, taking the full inverse is a <aste of CPU. It's just that i don't know how to do otherwise yet.

Pardon my ignorance

Bonne soirée :)

Jeff

Le 10/07/2013 19:42, Pascal Dupuis-3 [via Octave] a écrit :

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 <[hidden email]>:

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

---

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