Re: [eigen] Matrix power

[ Thread Index | Date Index | More Archives ]

Sorry for my late reply.

Powering rank deficient matrices whose size >= 3x3 will reproduce this bug.

To compute A^p, we Schur-decompose A to U T U* and compute U T^p U*.  T is upper
triangular with eigenvalues of A on the diagonal.

For triangular matrix T whose size >= 3x3, We use continued fraction of (I - T)
to compute powers of T, which requires norm of (I - T) to be sufficiently small.
If norm(I-T) is too large, we sqrt (I - T) and then square the result.

If A has a zero eigenvalue, there is a zero on diag of T.  Therefore, there is an 1
on diag of (I - T), enlarging (I - T).norm(), making Eigen sqrt'ing (I - T).
However, no matter how many times (I - T) is sqrt'ed, the 1 stands still on its
diag, leading to an infinite loop.

I will finish my final exam at 8:50 tomorrow UTC.  I will work on this since

On Tuesday, 2013-06-18 20:08:32,Dale Lukas Peterson wrote:
> A simpler test case that works is:
> Eigen::MatrixXd A(2, 2);
> A <<  0, 0, 0, 1;

Powers of 2x2 triangular matrices are directly computed. :)

> Yet a just barely more complicated case that does not work (i.e., it hangs) is:
>   Eigen::MatrixXd A(3, 3);
>   A <<  0, 0, 0,
>         0, 0, 0,
>         0, 0, 1;
> Maybe this will help with tracking down the bug.
> Luke

Thanks for your discovery!

jdh8@xxxxxxxxxxxxxx |

Mail converted by MHonArc 2.6.19+