Hello all,

I am having a problem with eigen_unsymm.h. The problem is that it keeps returning incorrect eigenvectors but correct eigenvalues.

I am using the NAG C library documented test for checking.

http://http://www.nag.co.uk/numeric/CL/nagdoc_cl08/pdf/F02/f02agc.pdf

Here is my test 4x4 matrix.

1.5 0.1 4.5 −1.5
−22.5 3.5 12.5 −2.5
−2.5 0.3 4.5 −2.5
−2.5 0.1 4.5 2.5

Here is my simple program.

int main()
{

string filename;
ifstream input;
int n=0;

cout << "Enter the name of the matrix input file " << endl;
cin >> filename;
cout << endl;
input.open(filename.c_str());
if(input.fail())
cout << "Error opening matrix input file\n";
cout << "Enter the size of the matrix " << endl;
cin >> n;
cout << endl;

MatDoub a(n,n);
for(int i=0; i<a.nrows(); i++){
for(int j=0; j<a.ncols(); j++){
input >> a[i][j];
}
}

// print out the input matrix and vector
cout << fixed << setprecision(3);
cout << "The matrix read in is of size " << n << " x " << n <<" and contains the following elements" << endl << endl;
// print out the matrix that was read into the program
for(int i=0;i<n;i++) {
for(int j=0;j<n;j++) {
cout << setw(14) << a[i][j];
}
cout << endl;
}
cout << endl;

Unsymmeig h(a);
MatDoub res(n,n);
VecComplex comp(n);

comp = h.wri;
res = h.zz;

for(int i=0; i<n; i++){
cout << comp[i] << endl;
}
for(int i=0; i<n; i++){
for(int j=0; j<n; j++){
cout << res[i][j] << endl;
}
cout << endl;
}

double pause;
cin >> pause;

return 0;
}



Here are the returned eigenvalues. These are correct.

(4.000,0.000)
(3.000,4.000)
(3.000,-4.000)
(2.000,0.000)

Here are the returned eigenvectors. These are not correct.

0.103
-0.791
0.119
-0.404

-3.080
-2.852
-2.807
-6.434

-0.034
-0.570
-0.561
-0.037

-0.479
-0.791
0.119
-0.404

Here are the correct eigenvectors from the NAG library.

Eigenvectors
( 0.113, -0.151) ( 0.113, 0.151) ( -0.033, 0.000) ( 0.063, 0.000)
( 0.945, 0.000) ( 0.945, 0.000) ( 0.988, 0.000) ( 0.996, 0.000)
( 0.189, 0.000) ( 0.189, 0.000) ( 0.011, 0.000) ( 0.006, 0.000)
( 0.113, -0.151) ( 0.113, 0.151) ( 0.154, 0.000) ( 0.063, 0.000)

In reading the NR3 book, it mentions that some compiliers have trouble with Complex math. I am not sure if that is my problem or not.

I am using a Windows XP 32 bit machine with Microsoft Visual Studio 2008.

Any help would be appreciated.

Thanks
Matt

Hi Matt,

You're printing out the eigenvectors incorrectly. The line
cout << res[i][j] << endl;
should be
cout << res[j][i] << endl;

The eigenvectors returned by the NR code are unnormalized. To make them look like the NAG values you will have to divide each vector by its norm (and perhaps change the overall sign). For the eigenvectors belonging to the complex eigenvalues, you use two consecutive vectors to get the real and imaginary parts, as described in the text.

Thank you for your help.

Matt

Matt,

I can confirm that your returned eigenvalues are correct but I'm not sure of the notation in which you show them. MathCAD output of eigenvals(Q):

3 + 4i
3 - 4i
4
2

The inversion I get might be the "normalized" look which kibitzer refers to, as accounted for in MathCAD. Actually, I tried inverting the output of your cout per his suggestion and didn't get your "correct" values either, so ...

MathCAD eigenvecs(Q):

0.113 - 0.151i 0.113 + 0.151i -0.033 0.063
0.945 0.945 0.988 0.996
0.189 0.189 0.011 5.692e10^-3
0.113 - 0.151i 0.113 + 0.151i 0.154 0.063

Hope this helps.

(compiled your code from above on a Windows Server 2008 64-bit OS with Visual Studio 2010 and NR3(latest))
MathCAD 6.0e (whoopie !!!) on XP 32-bit ... good fer what ails ya anyway ...

Yes, even with the output changed the results are not what the NAG library example shows.

Thanks for checking it with MathCad. I will have to dust off my MathCad program and give it a try.

Thanks
Matt

Hi Matt,

The NR3 code is correct and gives the same results as NAG. After you fix your code as I described, you also have to normalize the eigenvectors so that their magnitude is unity. (Recall that the overall scale and sign of the eigenvectors is arbitrary.)

results are not what the NAG library example shows.

Matt: The eigenvectors are not unique.

For a given eigenvalue (real or complex) we know that multiplying an eigenvector by a constant gives us something that is still an eigenvector.

However, and additionally...

For complex eigenvalues of real matrices there is another degree of freedom in obtaining eigenvectors, so, for a given eigenvalue, there are other eigenvectors not necessarily related by constant multipliers.

For a simple example, completely worked out, see this Complex Eigenvalues (http://www.google.com/url?sa=t&source=web&cd=1&ved=0CBMQFjAA&url=http%3A%2F%2Fwww.math.vt.edu%2Fpeople%2Fafkham is%2Fclass_home%2Fnotes%2FF08W11b.pdf&ei=NfVCTYymJJP6sAPc5cy6Cg&usg=AFQjCNHOGnNG2W9Ob7_hZWOhAV4wIfea0w) paper. Look at the results that are shown little less than halfway down the third page.

Now, programs like the ones that you are comparing give numbers as outputs. That's what you get, and that's all that you get. (See Footnote.)

Then the question should not be whether the results from the NR program are "the same as" results from another program, but whether the results from the NR program are valid.

Bottom line: If you normalize the eigenvectors from the NR program, the ones corresponding to real eigenvalues should will match the normalized eigenvalues from other correct implementations. (Well...the signs may be different, but to me that's pretty much "the same as" in eigenspeak.)

However...

The normalized eigenvectors of complex eigenvalues may not be "the same as" the other ones, but they will still be valid.

I'll show some results in the next post if you are still interested. (It's kind of lengthy.)

LAAEFRIF (Left As An Exercise For Really Interested Folks): Add stuff to your NR main program to extract, normalize, and verify the results without writing a file to be used by a separate program.

Regards,

Dave

Footnote:
"The purpose of computing is insight, not numbers."
---Richard W. Hamming

I'll show some results in the next post if you are still interested.


I started with your NR program and printed results to four significant decimal digits:

//
// (Your stuff to read in the matrix)
//
// print out the input matrix and vector
cout << "The matrix read in is of size " << n
<< " x " << n << " and contains the following elements"
<< endl << endl;
// print out the matrix that was read into the program
cout << scientific << showpos << setprecision(3);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
cout << " "<< a[i][j];
}
cout << endl;
}
cout << endl;

Unsymmeig h(a);

VecComplex lambda(n); // hold the eigenvalues
MatDoub x(n, n); // The output with eigenvectors

lambda = h.wri;
x = h.zz;

cout << "Eigenvalues:" << endl;
for (int i = 0; i < n; i++) {
cout << " " << lambda[i];
if (lambda[i].imag() == 0) {
cout << " (real)";
}
else {
cout << " (complex)";
}
cout << endl;
}

cout << "Eigenvectors are obtained from the columns:" << endl;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
cout << " " << x[i][j];
}
cout << endl;
}

return 0;
}



Output:

The matrix read in is of size 4 x 4 and contains the following elements

+1.500e+00 +1.000e-01 +4.500e+00 -1.500e+00
-2.250e+01 +3.500e+00 +1.250e+01 -2.500e+00
-2.500e+00 +3.000e-01 +4.500e+00 -2.500e+00
-2.500e+00 +1.000e-01 +4.500e+00 +2.500e+00

Eigenvalues:
(+4.000e+00,+0.000e+00) (real)
(+3.000e+00,+4.000e+00) (complex)
(+3.000e+00,-4.000e+00) (complex)
(+2.000e+00,+0.000e+00) (real)
Eigenvectors are obtained from the columns:
+1.027e-01 -7.913e-01 +1.194e-01 -4.044e-01
-3.080e+00 -2.852e+00 -2.807e+00 -6.434e+00
-3.422e-02 -5.703e-01 -5.614e-01 -3.677e-02
-4.790e-01 -7.913e-01 +1.194e-01 -4.044e-01


Then I put the eigenvalues and eigenvectors obtained from this output into a separate C++ program that normalized the eigenvectors and printed results of a*x and lambda*x for each eigenvalue/vector. This uses standard C++ complex arithmetic for all arithmetic operations.

#include <iostream>
#include <iomanip>
#include <complex>

using namespace std;

typedef complex <double> dcomplex;

inline dcomplex dc(double rr, double ii)
{
return dcomplex(rr, ii);
}

double normal(dcomplex x[4])
{
double n = 0.0;
for (int i = 0; i < 4; i++) {
n += norm(x[i]);
}
return n;
}

int main()
{
// The real matrix formulated with complex elements
dcomplex a[4][4] = {
{dc( 1.5, 0),dc(0.1, 0),dc( 4.5,0), dc(-1.5, 0)},
{dc(-22.5, 0),dc(3.5, 0),dc(12.5,0), dc(-2.5, 0)},
{dc( -2.5, 0),dc(0.3, 0),dc( 4.5,0), dc(-2.5, 0)},
{dc( -2.5, 0),dc(0.1, 0),dc( 4.5,0), dc( 2.5, 0)}
};

// The eigenvalues
dcomplex lambda[4] = {
dc(4.000, 0.000),
dc(3.000, 4.000),
dc(3.000,-4.000),
dc(2.000, 0.000),
};

// Columns are eigenvectors
// These are obtained from the results of the NR program
dcomplex x[4][4] = {
{dc(+1.027e-01, 0), dc(-7.913e-01, +1.194e-01), dc(-7.913e-01, -1.194e-01), dc(-4.044e-01, 0)},
{dc(-3.080e+00, 0), dc(-2.852e+00, -2.807e+00), dc(-2.852e+00, +2.807e+00), dc(-6.434e+00, 0)},
{dc(-3.422e-02, 0), dc(-5.703e-01, -5.614e-01), dc(-5.703e-01, +5.614e-01), dc(-3.677e-02, 0)},
{dc(-4.790e-01, 0), dc(-7.913e-01, +1.194e-01), dc(-7.913e-01, -1.194e-01), dc(-4.044e-01, 0)}
};

dcomplex lambdax[4][4];
dcomplex ax[4][4];

cout << showpos << scientific << setprecision(3);
cout << "a: " << endl;
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
cout << a[i][j] << " ";
}
cout << endl;
}
cout << endl;
for (int i = 0; i < 4; i++) {
cout << noshowpos << "lambda[" << i
<< showpos << "] = " << lambda[i] << endl;
}
cout << endl << endl;

cout <<"Un-normalized eigenvectors are the columns:" << endl;
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
cout << x[i][j] << " ";
}
cout << endl;
}
cout << endl;
for (int j = 0; j < 4; j++) {
double sum = 0.0;
for (int i = 0; i < 4; i++) {
sum += norm(x[i][j]);
}
cout << noshowpos << "norm of x[" << j << "] = "
<< showpos << sum << endl;
for (int i = 0; i < 4; i++) {
x[i][j] /= sqrt(sum);
}
}
cout << endl;
cout <<"Normalized eigenvectors are the columns):" << endl;
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
cout << x[i][j] << " ";
}
cout << endl;
}
cout << endl;
for (int j = 0; j < 4; j++) {
double sum = 0.0;
for (int i = 0; i < 4; i++) {
sum += norm(x[i][j]);
}
cout << noshowpos << "norm of x[" << j << "] = "
<< showpos << sum << endl;
}
cout << endl;



// Multiply each j column of x by lambda[j];
for (int j = 0; j < 4; j++) {
for (int i = 0; i < 4; i++) {
lambdax[i][j] = lambda[j] * x[i][j];
}
}
cout << "Each column of x multiplied by the corresponding lambda:" << endl;
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
cout << lambdax[i][j] << " ";
}
cout << endl;
}
cout << endl;

for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
ax[i][j] = 0.0;
for (int k = 0; k < 4; k++) {
ax[i][j] += a[i][k] * x[k][j];
}
}
}

cout << "a * x: " << endl;
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
cout << ax[i][j] << " ";
}
cout << endl;
}

return 0;
}


Output:

a:
(+1.500e+00,+0.000e+00) (+1.000e-01,+0.000e+00) (+4.500e+00,+0.000e+00) (-1.500e+00,+0.000e+00)
(-2.250e+01,+0.000e+00) (+3.500e+00,+0.000e+00) (+1.250e+01,+0.000e+00) (-2.500e+00,+0.000e+00)
(-2.500e+00,+0.000e+00) (+3.000e-01,+0.000e+00) (+4.500e+00,+0.000e+00) (-2.500e+00,+0.000e+00)
(-2.500e+00,+0.000e+00) (+1.000e-01,+0.000e+00) (+4.500e+00,+0.000e+00) (+2.500e+00,+0.000e+00)

lambda[0] = (+4.000e+00,+0.000e+00)
lambda[1] = (+3.000e+00,+4.000e+00)
lambda[2] = (+3.000e+00,-4.000e+00)
lambda[3] = (+2.000e+00,+0.000e+00)


Un-normalized eigenvectors are the columns:
(+1.027e-01,+0.000e+00) (-7.913e-01,+1.194e-01) (-7.913e-01,-1.194e-01) (-4.044e-01,+0.000e+00)
(-3.080e+00,+0.000e+00) (-2.852e+00,-2.807e+00) (-2.852e+00,+2.807e+00) (-6.434e+00,+0.000e+00)
(-3.422e-02,+0.000e+00) (-5.703e-01,-5.614e-01) (-5.703e-01,+5.614e-01) (-3.677e-02,+0.000e+00)
(-4.790e-01,+0.000e+00) (-7.913e-01,+1.194e-01) (-7.913e-01,-1.194e-01) (-4.044e-01,+0.000e+00)

norm of x[0] = +9.728e+00
norm of x[1] = +1.793e+01
norm of x[2] = +1.793e+01
norm of x[3] = +4.172e+01

Normalized eigenvectors are the columns):
(+3.293e-02,+0.000e+00) (-1.869e-01,+2.819e-02) (-1.869e-01,-2.819e-02) (-6.261e-02,+0.000e+00)
(-9.875e-01,+0.000e+00) (-6.735e-01,-6.628e-01) (-6.735e-01,+6.628e-01) (-9.961e-01,+0.000e+00)
(-1.097e-02,+0.000e+00) (-1.347e-01,-1.326e-01) (-1.347e-01,+1.326e-01) (-5.692e-03,+0.000e+00)
(-1.536e-01,+0.000e+00) (-1.869e-01,+2.819e-02) (-1.869e-01,-2.819e-02) (-6.261e-02,+0.000e+00)

norm of x[0] = +1.000e+00
norm of x[1] = +1.000e+00
norm of x[2] = +1.000e+00
norm of x[3] = +1.000e+00

Each column of x multiplied by the corresponding lambda:
(+1.317e-01,+0.000e+00) (-6.733e-01,-6.628e-01) (-6.733e-01,+6.628e-01) (-1.252e-01,+0.000e+00)
(-3.950e+00,+0.000e+00) (+6.309e-01,-4.682e+00) (+6.309e-01,+4.682e+00) (-1.992e+00,+0.000e+00)
(-4.389e-02,+0.000e+00) (+1.263e-01,-9.364e-01) (+1.263e-01,+9.364e-01) (-1.138e-02,+0.000e+00)
(-6.143e-01,+0.000e+00) (-6.733e-01,-6.628e-01) (-6.733e-01,+6.628e-01) (-1.252e-01,+0.000e+00)

a * x:
(+1.316e-01,+0.000e+00) (-6.733e-01,-6.628e-01) (-6.733e-01,+6.628e-01) (-1.252e-01,+0.000e+00)
(-3.950e+00,+0.000e+00) (+6.309e-01,-4.682e+00) (+6.309e-01,+4.682e+00) (-1.992e+00,+0.000e+00)
(-4.400e-02,+0.000e+00) (+1.262e-01,-9.364e-01) (+1.262e-01,+9.364e-01) (-1.140e-02,+0.000e+00)
(-6.144e-01,+0.000e+00) (-6.733e-01,-6.628e-01) (-6.733e-01,+6.628e-01) (-1.252e-01,+0.000e+00)


Note that slight numeric differences in a few places are due to the fact that the eigenvector values from the NR program have been rounded to four significant decimal digits.

Bottom line:
The normalized eigenvectors for the real eigenvalues match the ones of your example (except for the signs). The complex ones are different, but they are still valid.



Regards,

Dave

Footnote:
"The purpose of computing is insight, not numbers."
---Richard W. Hamming

Thank you very much for all the work you put into this. I will get this up and running on my machine today.

Thanks
Matt

The eigenvectors *are* determined up to a constant factor even when complex (at least in cases like this one, where there are n distinct eigenvalues). It's just that the constant factor can be complex. For instance, here's one of Matt's complex eigenvectors after normalization by Dave's code:

(-1.869e-01,+2.819e-02)
(-6.735e-01,-6.628e-01)
(-1.347e-01,-1.326e-01)
(-1.869e-01,+2.819e-02)

and here's one of the complex eigenvectors from the NAG document:

( 0.113, -0.151)
( 0.945, 0.000)
( 0.189, 0.000)
( 0.113, -0.151)

and here are the ratios between corresponding elements of those two vectors, which you can see are all equal apart from rounding errors:

(-0.713,-0.704)
(-0.713,-0.701)
(-0.713,-0.702)
(-0.713,-0.704)

These eigenvectors *are* the same apart from a (complex) constant factor.

...
These eigenvectors *are* the same apart from a (complex) constant factor.
Thank you very much!

Regards,

Dave