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 eigenvalues and eigenvectors 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



Hi,
Could you post the matrix you are trying to decompose?
Thanks,
Marcus



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



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



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 usecases  please read the comments in the entry of each ifstate (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



This discussion has been copied to a work item. Click
here to go to the work item and continue the discussion.

