This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
Bug 614 - Matrix power segfaults
Summary: Matrix power segfaults
Status: RESOLVED FIXED
Alias: None
Product: Eigen
Classification: Unclassified
Component: Unsupported modules (show other bugs)
Version: 3.2
Hardware: x86 - 64-bit Linux
: Normal Unknown
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on:
Blocks: 3.2
  Show dependency treegraph
 
Reported: 2013-06-13 20:21 UTC by Dale Lukas Peterson
Modified: 2019-12-04 12:23 UTC (History)
5 users (show)



Attachments
Matrix power segfaults for trivial 2x2 matrix. (606 bytes, text/x-c++src)
2013-06-13 20:21 UTC, Dale Lukas Peterson
no flags Details
Proposed patch, including regression test (1.34 KB, patch)
2013-06-16 00:04 UTC, Jitse Niesen
chtz: review+
Details | Diff
Proposed patch, including regression test (2.27 KB, patch)
2013-06-21 06:52 UTC, Chen-Pang He
no flags Details | Diff
Proposed patch, including regression test (3.39 KB, patch)
2013-06-22 07:13 UTC, Chen-Pang He
jdh8: review? (jitseniesen)
Details | Diff

Description Dale Lukas Peterson 2013-06-13 20:21:03 UTC
Created attachment 347 [details]
Matrix power segfaults for trivial 2x2 matrix.

The attached file demonstrates pow() segfaulting on trivial matrix powers.  I am using g++ (Gentoo 4.7.3 p1.0, pie-0.5.5) 4.7.3.  Not all powers segfault, for example A.pow(1.0) does not.  When I run the debugger (after compiling with -ggdb3), and step through until it segfaults, I get the following backtrace:

~/tmp/c++/eigen$ gdb ./matrix_power
GNU gdb (Gentoo 7.6 p1) 7.6
Copyright (C) 2013 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-pc-linux-gnu".
For bug reporting instructions, please see:
<http://bugs.gentoo.org/>...
Reading symbols from /home/hazelnusse/tmp/c++/eigen/matrix_power...done.
(gdb) 
(gdb) list
1	// Compiled with:
2	//
3	// $ g++ -Wall matrix_power.cc -O2 -I/home/hazelnusse/usr/include/eigen3 -o matrix_power
4	//
5	#include <iostream>
6	#include <Eigen/Dense>
7	#include <unsupported/Eigen/MatrixFunctions>
8	int main()
9	{
10	  Eigen::MatrixXcd A(2,2);
(gdb) b main
Breakpoint 1 at 0x401cd0: file matrix_power.cc, line 9.
(gdb) run
Starting program: /home/hazelnusse/tmp/c++/eigen/matrix_power 
warning: Could not load shared library symbols for linux-vdso.so.1.
Do you need "set solib-search-path" or "set sysroot"?

Breakpoint 1, main () at matrix_power.cc:9
9	{
(gdb) n
10	  Eigen::MatrixXcd A(2,2);
(gdb) 
9	{
(gdb) 
10	  Eigen::MatrixXcd A(2,2);
(gdb) 
11	  A << 1, 1, 4, -2;
(gdb) 
20	  std::cout << A.pow(3.1) << "\n\n";
(gdb) print A
$1 = Eigen::Matrix<std::complex<double>,2,2,ColMajor> (data ptr: 0x632010) = {[0,0] = {_M_value = 1 + 0 * I}, [1,0] = {_M_value = 4 + 0 * I}, [0,1] = {_M_value = 1 + 0 * I}, [1,1] = {_M_value = -2 + 0 * I}}
(gdb) n

Program received signal SIGSEGV, Segmentation fault.
0x00000000004068b7 in Eigen::MatrixPowerTriangularAtomic<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::compute2x2 (this=0x7fffffffd5d0, 
    res=Eigen::Matrix<std::complex<double>,0,0,ColMajor> (data ptr: 0x0), p=0.10000000000000009) at /home/hazelnusse/usr/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h:330
330	  res.coeffRef(0,0) = pow(m_A.coeff(0,0), p);
(gdb) bt
#0  0x00000000004068b7 in Eigen::MatrixPowerTriangularAtomic<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::compute2x2 (this=0x7fffffffd5d0, 
    res=Eigen::Matrix<std::complex<double>,0,0,ColMajor> (data ptr: 0x0), p=0.10000000000000009) at /home/hazelnusse/usr/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h:330
#1  0x00000000004165db in Eigen::MatrixPowerTriangularAtomic<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::compute (this=this@entry=0x7fffffffd5d0, 
    res=Eigen::Matrix<std::complex<double>,0,0,ColMajor> (data ptr: 0x0), p=0.10000000000000009) at /home/hazelnusse/usr/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h:304
#2  0x0000000000417432 in Eigen::MatrixPower<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::computeFracPower<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> > (this=this@entry=0x7fffffffd6e0, 
    res=Eigen::Matrix<std::complex<double>,2,2,ColMajor> (data ptr: 0x6320c0) = {...}, p=<optimized out>) at /home/hazelnusse/usr/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h:459
#3  0x00000000004187ed in Eigen::MatrixPower<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::compute (this=this@entry=0x7fffffffd6e0, 
    res=Eigen::Matrix<std::complex<double>,2,2,ColMajor> (data ptr: 0x6320c0) = {...}, p=1.0717734625362931) at /home/hazelnusse/usr/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h:335
#4  0x0000000000418919 in Eigen::MatrixPowerReturnValue<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::evalTo<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> > (
    this=this@entry=0x7fffffffd8f0, res=Eigen::Matrix<std::complex<double>,2,2,ColMajor> (data ptr: 0x6320c0) = {...})
    at /home/hazelnusse/usr/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h:511
#5  0x00000000004189a6 in evalTo<Eigen::Matrix<std::complex<double>, -1, -1> > (dst=Eigen::Matrix<std::complex<double>,2,2,ColMajor> (data ptr: 0x6320c0) = {...}, this=0x7fffffffd8f0)
    at /home/hazelnusse/usr/include/eigen3/Eigen/src/Core/ReturnByValue.h:61
#6  Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1>::Matrix<Eigen::MatrixPowerReturnValue<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> > > (this=0x7fffffffd830, other=...)
    at /home/hazelnusse/usr/include/eigen3/Eigen/src/Core/Matrix.h:296
#7  0x0000000000418aaa in eval (this=0x7fffffffd8f0) at /home/hazelnusse/usr/include/eigen3/Eigen/src/Core/DenseBase.h:357
#8  Eigen::operator<< <Eigen::ReturnByValue<Eigen::MatrixPowerReturnValue<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> > > > (s=..., m=...)
    at /home/hazelnusse/usr/include/eigen3/Eigen/src/Core/IO.h:244
#9  0x0000000000401de9 in main () at matrix_power.cc:20
(gdb) 

I'm not sure how to proceed next...
Comment 1 Gael Guennebaud 2013-06-13 22:23:13 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.
Comment 2 Jitse Niesen 2013-06-16 00:04:39 UTC
Created attachment 348 [details]
Proposed patch, including regression test
Comment 3 Christoph Hertzberg 2013-06-17 11:30:28 UTC
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.
Comment 4 Chen-Pang He 2013-06-19 06:53:24 UTC
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.
Comment 5 Christoph Hertzberg 2013-06-19 08:09:03 UTC
(In reply to comment #4)
>     res.resize(m_A.rows(), m_A.cols());
This is equivalent to:
    res.resizeLike(m_A);
Comment 6 Chen-Pang He 2013-06-21 06:52:26 UTC
Created attachment 352 [details]
Proposed patch, including regression test

This also fixes the infinite loop in computeBig (see mailing list).
Comment 7 Gael Guennebaud 2013-06-21 08:39:22 UTC
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?
Comment 8 Chen-Pang He 2013-06-21 12:51:50 UTC
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.
Comment 9 Chen-Pang He 2013-06-22 07:13:20 UTC
Created attachment 355 [details]
Proposed patch, including regression test

Sorry that compute2x2 produced wrong result for adjacent equal eigenvalues.
Comment 10 Chen-Pang He 2013-06-22 07:16:01 UTC
Comment on attachment 355 [details]
Proposed patch, including regression test

Missing semicolon at eigen_assert(m_A(i,i) != RealScalar(0))
Comment 11 Jitse Niesen 2013-06-24 15:12:13 UTC
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).
Comment 12 Nobody 2019-12-04 12:23:47 UTC
-- 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.

Note You need to log in before you can comment on or make changes to this bug.