Problem with TDecompSVD

ROOT
1

Hi there,

I am using svd to invert matrices. For the majority of events it appears to work ok but every so often I get the following message which causes my simulations to crash:

[quote]Error in : no convergence after 31 steps
Fatal in TMatrixD::: IsValid() violated at line 354 of `matrix/src/TMatrixDBase.cxx’
aborting
Generating stack trace…
0x0033df84 in TMatrixDBase::ResizeTo(int, int) + 0x62 from /usr/local/root/lib/libMatrix.so
0x00329c58 in TDecompBase::Invert() + 0x78 from /usr/local/root/lib/libMatrix.so
0x08053bc9 in Hit::ReconstructVertex(Neutrino, Array*, Visuals*, int, int) + 0x875 from ./NeutrinoPulse
0x080521b0 in Hit::RegisterHits(Spectrum*, Array*, Visuals*, int, int, int) + 0x166c from ./NeutrinoPulse
0x0804c831 in main + 0x609 from ./NeutrinoPulse
0x05020770 in __libc_start_main + 0xf0 from /lib/tls/libc.so.6
0x0804c199 in TFile::TFile[in-charge](char const*, char const*, char const*, int) + 0x41 from ./NeutrinoPulse
Aborted[/quote]

Is there a way I can test the matrix before I invert it to prevent this from happening? I have played around with the IsValid(), tolerance and the condition but haven’t had any success.

Cheers

Jon

I run root 4.00/04 on Linux Fedora Core 1 with g++ 3.3.2

2

Hi Jon,

Step through the decomposition by hand so that you can catch
these cases where the diagonalization fails :

I assume you did something like

TDecompSVD svd(a);
TMatrixD b = svd.Invert();

change that to :

TDecompSVD svd(a);
Bool_t ok = svd.Decompose();
TMatrixD b;
if (ok)
b = svd.Invert();
else {
cout << "SVD failed, condition: " << svd.Condition() <<endl;
a.Print();
}

It is a curious that SVD fails, would be interesting to see what the print out gives in these conditions !

Eddy

3

Many thanks, this works a treat!

Here is an exapmple of a failed matrix:

[code] | 0 | 1 | 2 |

0 | 0 1000 0
1 | 0 2000 0
2 | 0 3000 0
3 | 0 4000 0
4 | 0 5000 0
5 | 0 6000 0
6 | 0 7000 0
7 | 0 8000 0[/code]

The columns represent the difference in x, y and z coordinate between 1 reference sensor and the other sensors (in a 3-dimensional lattice-like array) that are hit during an event. By measuring the difference in arrival times of the signals we can reconstruct the vertex location. A matrix like the one above indicates all the sensors that are hit lie in the same plane along one axis direction- which would not allow for a 3 dimensional vertex reconstruction.

The Decompose() method thus proves to be good method for cutting these events!

Regards

Jon

4

I’m facing the same problem. I have to SVD decompose a lot of 3x3 matrices. Strangely, some fail. This happens for matrices with always the same configuration : all elements equal to 0 and 2 different of 0 on the second column. Exemple:

root [0] TMatrixD aa(3,3);
root [1] aa(1,1)=1;
root [2] aa(0,1)=1;
root [3] aa.Print()
      |        0  |        1  |        2  |
------------------------------------------------------------------
   0 |          0           1           0
   1 |          0           1           0
   2 |          0           0           0

root [4] TDecompSVD svd(aa);
root [5] svd.Decompose();
Error in <TDecompSVD::Diagonalize>: no convergence after 31 steps

To my knowledge, it is always possible to SVDecompose an arbitrary matrix. For instance, it this case, the solution can easily be computed “by hand” :

U =     Sqrt(1/2)   -Sqrt(1/2)   0
        Sqrt(1/2)   Sqrt(1/2)    0
        0           0            1

V =     0           1          0
        -1          0          0
        0           0          1

S = Sqrt(2)      0         0

That’s really strange that Root cannot find a solution for this kind of matrices. Does someone know any way to bypass the problem and find a solution?

Thanks a lot

PS : I’m also looking for a SVDecomposition of imaginary martix.
PS : I’m using Root 5.06/00

5

Hi Goffinet,

You seem to be facing the same problem as Jon above, the
diagonalization is not converging in case of a particular singular
matrix .

Most SVD algorithms on the market are variations on the
work by Golub and Kahn . So currently I am looking into
alternatives for my implementation of the diagonalization part .
Fairly tough problem , so it is going to take me some time .

Concerning your inquiry about decomposing complex matrices .
We are releasing in January a templated version of the matrix
package . However, things like matrix norms have not been
adjusted (yet ?) for complex numbers .
If you need this feature to solve complex linear systems, you
can also the current implementation in a slightly less
efficient way:

suppose : (A + i B) (x + i y) = w + i z

this can be cast in the following form :

        ( A  -B )  ( x )  = ( w )
        ( B   A )  ( y )       ( z )

of course now your n x n matrix is 2n x 2n .

Eddy

6

Hi Goffinet,

You seem to be facing the same problem as Jon above, the
diagonalization is not converging in case of a particular singular
matrix .

Most SVD algorithms on the market are variations on the
work by Golub and Kahn . So currently I am looking into
alternatives for my implementation of the diagonalization part .
Fairly tough problem , so it is going to take me some time .

Concerning your inquiry about decomposing complex matrices .
We are releasing in January a templated version of the matrix
package . However, things like matrix norms have not been
adjusted (yet ?) for complex numbers .
If you need this feature to solve complex linear systems, you
can also the current implementation in a slightly less
efficient way:

suppose : (A + i B) (x + i y) = w + i z

this can be cast in the following form :

        ( A  -B )  ( x )  = ( w )
        ( B   A )  ( y )      ( z )

of course now your n x n matrix is 2n x 2n .

Eddy

7

Hi,

In CVS a patch has been installed to TDecompSVD.cxx so that the
diagonalization process converges in case of a singular matrix with
lots of zeros .

{
  TMatrixD m(3,3);
  m(0,1) = m(1,1) = 1;
  TDecompSVD svd(m);
  svd.Decompose();
  TMatrixD U = svd.GetU();
  TMatrixD V = svd.GetV();
  TVectorD S = svd.Get();
  U.Print();
  V.Print();
  S.Print();
}

gives

3x3 matrix  is as follows

     |        0  |        1  |        2  |
------------------------------------------------------------------
   0 |     0.7071      0.7071     0 
   1 |     0.7071     -0.7071     0 
   2 |          0           0           1 

3x3 matrix  is as follows

     |        0  |               1  |        2  |
------------------------------------------------------------------
   0 | -1.009e-16              1               -0 
   1 |          1       8.865e-17              -0 
   2 |         -0                  -0           1 

Vector (3)  is as follows

     |        1  |
------------------
   0 |1.41421 
   1 |2.23711e-17 
   2 |0

Eddy