Hi guys,

I have the following puzzle. My function computes a matrix sign, here is the full code

Code: Select all

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>

using namespace Eigen;
using std::cout; using std::endl;

template <typename T>
inline int sgn(T val)
{
    return (T() < val) - (val < T());
}


void sign(const Ref<const MatrixXd> A, Ref<MatrixXd> signA)
{
    RealSchur<MatrixXd> schur(A);

    MatrixXd T = schur.matrixT();
    //auto U = schur.matrixU();

    int n = static_cast<int>(A.rows());

    MatrixXd Q(n, n);

    for(int i = 0; i < n; ++i)
        Q(i,i) = sgn(T(i,i));

    for(int p = 1; p < n; ++p)
    {
        for (int i = 0; i < n-p; ++i)
        {
            int j = i + p;

            if ( Q(i,i) + Q(j,j) != 0 )
            {
                MatrixXd temp = -Q.block(i,i+1,1,j-i-1) * Q.block(i+1,j,j-i-1,1) / (Q(i,i) + Q(j,j));
                Q(i,j) = temp(0,0);
            }
            else
            {
                double d = T(j,j) - T(i,i);
                double q = T(i,j) * (Q(j,j) - Q(i,i));
                if (p > 1)
                {
                    MatrixXd temp = T.block(i,i+1,1,j-i-1) * Q.block(i+1,j,j-i-1,1) - Q.block(i,i+1,1,j-i-1) * T.block(i+1,j,j-i-1,1);
                    q += temp(0,0);
                }
                Q(i,j) = q/d;
            }
        }
    }
    signA = schur.matrixU() * Q * schur.matrixU().transpose();
}

int main()
{
    int n = 3;
    MatrixXd A(n,n), signA(n,n);

    A << 1.0/1.0, 1.0/1.0, 1.0/1.0,
         1.0/2.0, 1.0/3.0, 1.0/4.0,
         1.0/3.0, 1.0/4.0, 1.0/5.0;

    sign(A, signA);
    cout << signA << endl;
}


This code calculates the result incorrectly. However, if I simply uncomment the line "auto U = schur.matrixU();", the function produces a correct result. Could you please explain to me why this happens exactly and how I can trigger the correct behavior without adding redundant code lines?.. By the way, is there a better way to assign a 1x1 matrix to a double?

I get the following answer if the auto line is commented out:

0.25039 1.0942 1.01117
0.521153 -0.543945 0.421448
0.363004 0.31766 -0.706445

I think this is the correct answer. So I cannot reproduce.

oh wait, valgrind produces a lot of errors (accessing an uninitialized value). The issue seems to be that you never set the lower triangular part of Q to zero. Once I fixed that, the valgrind errors disappear.

For your information, bug 31 http://eigen.tuxfamily.org/bz/show_bug.cgi?id=31 asks for an implementation of the matrix sign function, so if you get it working I would be interested.

Thanks, you were right. Initializing Q to a zero matrix did the trick. This particular algorithm for calculating the matrix sign is slow but stable. The scaled Newton iteration that you mention in Bug 31 should be faster as long as the matrix is not pathological. It should be straightforward to code.

Does your algorithm also work if the matrix has complex eigenvalues?

Regarding your other question, you can also write "q += temp(0,0)" as "q += temp.value()" but whether that is clearer is open to discussion. Perhaps better is to write "T.block(i,i+1,1,j-i-1)" as "T.row(i).segment(i+1,j-i-1)" so that you get a vector and can compute the dot product which gives you a scalar - this also gives Eigen and the compiler opportunities to optimize.

Yeah, it will work for complex eigenvalues, if you use the ComplexSchur decomposiiton rather than RealSchur.