MPD78
01-26-2011, 07:40 AM
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
davekw7x
01-28-2011, 11:23 AM
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
davekw7x
01-28-2011, 11:31 AM
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
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.