eigenvalue decomposition issue (newbie question)

Jun 27, 2011 at 6:40 PM
Edited Jun 27, 2011 at 6:54 PM

Hi

I'm trying to decompose a complex 3x3 (Hermitian) matrix into eigen-values and eigen-vectors and I'm doing it this way:

 

var eigVR = Css.Evd();

 

VS2010 compiles fine but when running the application it gets stuck at this line (halts). I put a break point after this line which is never reached!

Any hints on what I have done wrong?

ADDED:

Also, when I do this:

 

 UserEvd eigVR = new UserEvd(Css);

The same happens!

 

 

 

thx

Jun 28, 2011 at 5:35 AM

Hi,

Could you post the matrix you are trying to decompose? 

Thanks,
Marcus 

Jun 28, 2011 at 12:42 PM

Hi Marcus,

Thanks for your support!

The code looks like this:

DenseMatrix[] mArrayOld;

mArrayOld = new DenseMatrix[256];
for (int k = 0; k < 256; k++)
{
    mArrayOld[k] = new DenseMatrix(3,3);            
}
// Breakpoint before the decomposition shows the content of the first DenseMatrix in the array of dense matrices (index = [0] )
base {MathNet.Numerics.LinearAlgebra.Complex32.Matrix} = {(0.9275488, 0),(0, 0),(0, 0)(-0.9275488, 0),(0.9275488, 0),(0, 0)(0, 0),(-0.9275488, 0),(0.9275488, 0)}�

UserEvd eigVR = new UserEvd(mArrayOld[0]); // halts the program

This works though:

 

DenseMatrix M = new DenseMatrix(3);

UserEvd eigVR = new UserEvd(M);

thx

Jun 28, 2011 at 1:46 PM

Hi Marcus,

I'm getting further into this issue. 

Please wait looking into this till I report back...I think I will have some valuable information for you!

Thanks

Jun 28, 2011 at 1:51 PM

Hi,

I took a quick look.

I think there is a bug in the Complex/Complex32 version of our EVD code. At one point we are doing:

s = Math.Abs(matrixH[n, n - 1].Real) + Math.Abs(matrixH[n - 1, n - 2].Real);

but n = 1.

Not sure why we are getting there with the matrix you provided. If I use the real versions of the EVD code, then there is no problem.

>Please wait looking into this till I report back..

OK. I'll wait until you get back to me.

Thanks,
Marcus 

Jun 28, 2011 at 2:11 PM
Edited Jun 28, 2011 at 2:12 PM

Hi Marcus,

I'm buidling a lot of complex 3x3 covariance matrices every second and I have found a few issues:

When I build the covariance matrix having only unique entries in the lower triangle of the matrix the inversion step seems efficient when done like this:

 

mArrayOld = (DenseMatrix)mArrayOld.Inverse();

 

but the Evd then halts. I know the inversion of the partially populated covariance matrix doesn't make much sense ...but you probably want to know this.

If I build the full covariance matrix the inversion step above is very slow but the evd then works fine (and seems efficient)

I tested two use-cases - please read the comments in the entry of each if-state (i.e. partial vs. full covariance matrix)

 

                if(false)
                {
                    // // partial Covariance Matrix - inversion fast but decomp halts the program!
                    // Compute covarians matrix
                    Complex32 tau1 = new Complex32(0.1f, 0f);
                    Complex32 tau2 = new Complex32(0.9f, 0f);
                    for (int k = 0; k < N / 2 + 1; k += 2)
                    {
                        m1 = new Complex32(Y[0, k], Y[0, k + 1]);
                        m2 = new Complex32(Y[1, k], Y[1, k + 1]);
                        m3 = new Complex32(Y[2, k], Y[2, k + 1]);
                        mArray[k / 2][0, 0] = m1.Multiply(new Complex32(m1.Real, -m1.Imaginary));   // Diagonal term
                        mArray[k / 2][1, 1] = m2.Multiply(new Complex32(m2.Real, -m2.Imaginary));   // Diagonal term
                        mArray[k / 2][2, 2] = m3.Multiply(new Complex32(m3.Real, -m3.Imaginary));   // Diagonal term
                        mArray[k / 2][1, 0] = m1.Multiply(new Complex32(m2.Real, -m2.Imaginary));
                        mArray[k / 2][2, 0] = m1.Multiply(new Complex32(m3.Real, -m3.Imaginary));
                        mArray[k / 2][2, 1] = m2.Multiply(new Complex32(m3.Real, -m3.Imaginary));

                        // Smooth it (tau = 50ms)
                        mArrayOld[k / 2] = mArray[k / 2] * tau1 + mArrayOld[k / 2] * tau2;

                        // Matrix inversion
                        mArrayOld[k / 2] = (DenseMatrix)mArrayOld[k / 2].Inverse();

                        eigVR = new UserEvd(mArray[k/2]);
                    }
                }
                else
                {
                    // full Covariance Matrix - inversion step slow but decomp seems fast
                    // Compute covarians matrix
                    Complex32 tau1 = new Complex32(0.1f, 0f);
                    Complex32 tau2 = new Complex32(0.9f, 0f);
                    for (int k = 0; k < N / 2 + 1; k += 2)
                    {
                        m1 = new Complex32(Y[0, k], Y[0, k + 1]);
                        m2 = new Complex32(Y[1, k], Y[1, k + 1]);
                        m3 = new Complex32(Y[2, k], Y[2, k + 1]);
                        mArray[k / 2][0, 0] = m1.Multiply(new Complex32(m1.Real, -m1.Imaginary));   // Diagonal term
                        mArray[k / 2][1, 1] = m2.Multiply(new Complex32(m2.Real, -m2.Imaginary));   // Diagonal term
                        mArray[k / 2][2, 2] = m3.Multiply(new Complex32(m3.Real, -m3.Imaginary));   // Diagonal term
                        mArray[k / 2][1, 0] = mArray[k / 2][0, 1] = m1.Multiply(new Complex32(m2.Real, -m2.Imaginary));
                        mArray[k / 2][2, 0] = mArray[k / 2][0, 2] = m1.Multiply(new Complex32(m3.Real, -m3.Imaginary));
                        mArray[k / 2][2, 1] = mArray[k / 2][1, 2] = m2.Multiply(new Complex32(m3.Real, -m3.Imaginary));

                        // Smooth it (tau = 50ms)
                        mArrayOld[k / 2] = mArray[k / 2] * tau1 + mArrayOld[k / 2] * tau2;

                        // Matrix inversion
                        mArrayOld[k / 2] = (DenseMatrix)mArrayOld[k / 2].Inverse();

                        eigVR = new UserEvd(mArray[k / 2]);
                    }
                }

 

thx

 

 

 

 

Jun 28, 2011 at 2:50 PM
This discussion has been copied to a work item. Click here to go to the work item and continue the discussion.