Can anyone throw light on the following behaviour?
Consider the following program to diagonalize a block-
diagonal matrix.
#include <oxstd.h>
#include <oxfloat.h>
main()
{
decl eps = 0;
// decl eps = DBL_EPSILON;
decl mat1 = < 1,.5; .5, 1 >;
decl mat2 = < 2,.2; .2, 2 >;
decl mat3 = constant(eps,2,2);
decl mat = (mat1~mat3)|(mat3~mat2);
decl eval, evec;
eigensym(mat, &eval, &evec);
println(mat, evec*diag(eval)*evec');
}
This produces the following output, in both Ox3.4
and Ox 4.0.
1.0000 0.50000 0.00000 0.00000
0.50000 1.0000 0.00000 0.00000
0.00000 0.00000 2.0000 0.20000
0.00000 0.00000 0.20000 2.0000
1.0000 0.50000 0.00000 0.00000
0.50000 1.0000 0.00000 0.00000
0.00000 0.00000 2.0000 -0.20000
0.00000 0.00000 -0.20000 2.0000
Note the sign flip in the lower block!
Now comment out the first line of main(), and uncomment
the second line. This gives (as expected)
1.0000 0.50000 2.2204e-016 2.2204e-016
0.50000 1.0000 2.2204e-016 2.2204e-016
2.2204e-016 2.2204e-016 2.0000 0.20000
2.2204e-016 2.2204e-016 0.20000 2.0000
1.0000 0.50000 2.2204e-016 2.2204e-016
0.50000 1.0000 2.2204e-016 2.2204e-016
2.2204e-016 2.2204e-016 2.0000 0.20000
2.2204e-016 2.2204e-016 0.20000 2.0000
Replace mat1 and mat2 with the positive definite matrices
of your choice, you will observe the same effect.
The eigen() function appears to work correctly
in both situations, but replace evec' by invert(evec) to
recover the matrix.
Any comments/suggestions?
Regards
James
---------------------------------------------------------
James Davidson
Professor of Econometrics
University Of Exeter, School of Business and Economics
Exeter EX4 4PU
Web: http://www.ex.ac.uk/sobe/
Email: [log in to unmask]
Personal page: http://people.ex.ac.uk/jehd201/
---------------------------------------------------------
|