complex eigenvalues
Peter Dobcsanyi
p.dobcsanyi at designtheory.org
Sat May 29 21:07:23 BST 2004
Dear All,
As we are crunching out more and more designs it happens some times that
the numerical library used to compute the eigenvalues of the information
matrix returns complex numbers. Here is an example:
complex eigenvalue found ...
id: v7-b6-k3-109
blocks:
[[0, 1, 2], [0, 1, 2], [0, 1, 3], [2, 4, 6], [3, 4, 5], [3, 4, 5]]
incidence_matrix:
[[1 1 1 0 0 0]
[1 1 1 0 0 0]
[1 1 0 1 0 0]
[0 0 1 0 1 1]
[0 0 0 1 1 1]
[0 0 0 0 1 1]
[0 0 0 1 0 0]]
concurrence_matrix:
[[3 3 2 1 0 0 0]
[3 3 2 1 0 0 0]
[2 2 3 0 1 0 1]
[1 1 0 3 2 2 0]
[0 0 1 2 3 2 1]
[0 0 0 2 2 2 0]
[0 0 1 0 1 0 1]]
eigenvalues:
[ 2.22044605e-16 +0.00000000e+00j 6.66666667e-01 -1.85103892e-16j
6.66666667e-01 +1.85103892e-16j 2.33333333e+00 +0.00000000e+00j
2.33333333e+00 +0.00000000e+00j 3.00000000e+00 +0.00000000e+00j
3.00000000e+00 +0.00000000e+00j]
... skipped.
Of course, the imaginary parts are very small. Nevertheless, the
appearance of complex numbers causes problems in the computation of
statistical properties since the algorithms involved expect real
numbers.
To avoid these problems while maintaining as much numerical precision as
possible for the further computation, I have modified the block_design
Python module. From now on, all eigenvalues of the information matrix
returned by the numerical library are turned into nonnegative reals by
taking their absolute values. While this is not the perfect solution, we
certainly can do this in a "numerical sense" since, in theory, all
eigenvalues are nonnegative real numbers.
Using this method the eigenvalues for the above example become:
[2.2204460492503131e-16, 0.66666666666666674,
0.66666666666666674, 2.3333333333333335,
2.3333333333333335, 3.0000000000000004,
3.0000000000000027]
The perfect solution would be the exact computation, however writing
such code would take some time and money! ... you know what I mean :-)
Please comment on how you find this numerical solution for the time
being.
-- ,
Peter Dobcsanyi
More information about the Developers
mailing list