| Summary: | Matrix power segfaults | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Product: | Eigen | Reporter: | Dale Lukas Peterson <hazelnusse> | ||||||||||
| Component: | Unsupported modules | Assignee: | Nobody <eigen.nobody> | ||||||||||
| Status: | RESOLVED FIXED | ||||||||||||
| Severity: | Unknown | CC: | chtz, gael.guennebaud, hazelnusse, jdh8, jitseniesen | ||||||||||
| Priority: | Normal | ||||||||||||
| Version: | 3.2 | ||||||||||||
| Hardware: | x86 - 64-bit | ||||||||||||
| OS: | Linux | ||||||||||||
| Whiteboard: | |||||||||||||
| Bug Depends on: | |||||||||||||
| Bug Blocks: | 387 | ||||||||||||
| Attachments: |
|
||||||||||||
|
Description
Dale Lukas Peterson
2013-06-13 20:21:03 UTC
With your example given in bug 465 I observed it works fine if you use fixed size matrices, like Matrix2d instead of MatrixXd(2,2), and with MatrixXd it also works fine with integer powers. Created attachment 348 [details]
Proposed patch, including regression test
Comment on attachment 348 [details] Proposed patch, including regression test Review of attachment 348 [details]: ----------------------------------------------------------------- The patch solves the issue, however I think it should be specified better which method resizes its output parameter and which does not. E.g., computeBig() called in the default case automatically resizes, but case 1 does not (actually, case 1 would fail, but apparently is never called, because that case is handled elsewhere). I would accept the patch as a fix, but still leave the entire MatrixPower-module in a "needs review/cleanup"-status. The function compute2x2 is frequently used in computeBig, so calling resize in compute2x2 is inefficient. I suggest calling
res.resize(m_A.rows(), m_A.cols());
before switch. Sorry for my overlooking that res was not allocated.
(In reply to comment #4) > res.resize(m_A.rows(), m_A.cols()); This is equivalent to: res.resizeLike(m_A); Created attachment 352 [details]
Proposed patch, including regression test
This also fixes the infinite loop in computeBig (see mailing list).
just one question: why limiting the iteration count to 64? I mean why 64? Is there any good justification for this choice? Is it reasonable for any matrix orders? Cannot we easily detect the non convergence? Hmm... maybe just test if there is 0 on T's diag. T(i,i) == 0 seems okay because nonzero (even 1e-307) entries converge to 1 if continually sqrt'ed. I shouldn't have worried about "==" too much. The hard limit 64 is given in http://www.mathworks.com/matlabcentral/fileexchange/30527-arbitrary-real-power-of-a-matrix-by-schur-pade-algorithm/content/powerm_pade.m but not in the article. The MATLAB code is for double precision matrices. Created attachment 355 [details]
Proposed patch, including regression test
Sorry that compute2x2 produced wrong result for adjacent equal eigenvalues.
Comment on attachment 355 [details]
Proposed patch, including regression test
Missing semicolon at eigen_assert(m_A(i,i) != RealScalar(0))
Patch looks fine, I pushed it to the repository. changeset: 5380:01867724d43d tag: tip parent: 5378:1e6d2684d0b5 user: Chen-Pang He <jdh8@ms63.hinet.net> date: Mon Jun 24 13:58:51 2013 +0100 summary: Fix segfault and bug with equal eivals in matrix power (bug 614). Regarding your FIXME, the documentation says that the matrix power is defined via the matrix logarithm, and the logarithm is not defined for singular matrices, so one might argue that the matrix power does not need to accept singular matrices. Nevertheless, this would of course be good to implement. If you want, you can open a new bug/enhancement request for this. I also want to repeat what Christoph said: (In reply to comment #3) > The patch solves the issue, however I think it should be specified better which > method resizes its output parameter and which does not. E.g., computeBig() > called in the default case automatically resizes, but case 1 does not > (actually, case 1 would fail, but apparently is never called, because that case > is handled elsewhere). > I would accept the patch as a fix, but still leave the entire > MatrixPower-module in a "needs review/cleanup"-status. Thus I intend to make a cleaning pass through the whole MatrixFunctions module (not just the matrix power). -- GitLab Migration Automatic Message -- This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity. You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/614. |