MEX-file implementing Eigen library pseudo-inverse function crashes

Question

I am trying to implement an Eigen library pseudo-inverse function in a Matlab MEX-file. It compiles successfully but crashes when I run it.

I am trying to follow the FAQ on how to implement a pseudo-inverse function using the Eigen library.

The FAQ suggests adding it as a method to the JacobiSVD class but since you can't do that in C++ I'm adding it to a child class. It compiles successfully but then crashes without an error message. It successfully outputs "hi" without crashing if I comment out the line with the .pinv call so that's where the problem is arising. To run, I just compile it (as test.cpp) and then type test at the command line. I am using Matlab R2019a under MacOS 10.14.5 and Eigen 3.3.7. In my full code I also get lots of weird error messages regarding the pinv code but before I can troubleshoot I need this simple test case to work. This is all at the far limits of my understanding of C++. Any help appreciated.

#include "mex.h"
#include <Eigen/Dense>
#include <Eigen/SVD>
#include <cstdlib>
#include <cmath>
#include <iostream>

using namespace Eigen;
using namespace std;

//https://stackoverflow.com/questions/18804402/add-a-method-to-existing-c-class-in-other-file
class JacobiSVDext : public JacobiSVD<MatrixXf> {
    typedef SVDBase<JacobiSVD<MatrixXf>> Base;
    public:
    using JacobiSVD::JacobiSVD; //inherit constructors //https://stackoverflow.com/questions/347358/inheriting-constructors
    MatrixXf pinv() //http://eigen.tuxfamily.org/index.php?title=FAQ
    {
        eigen_assert(m_isInitialized && "SVD is not initialized.");
        double  pinvtoler=1.e-6; // choose your tolerance wisely!
        JacobiSVDext::SingularValuesType singularValues_inv=m_singularValues;
        for ( long i=0; i<m_workMatrix.cols(); ++i) {
            if ( m_singularValues(i) > pinvtoler )
                singularValues_inv(i)=1.0/m_singularValues(i);
            else singularValues_inv(i)=0;
        }
        return m_matrixV*singularValues_inv.asDiagonal()*m_matrixU.transpose();
    };
};

/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
             int nrhs, const mxArray *prhs[])
{
    MatrixXf X = MatrixXf::Random(5, 5);
    JacobiSVDext svd(X);
    MatrixXf Y=svd.pinv();
    cout << Y << endl;
    cout << "hi" << endl;
}

The expected result is to output the pseudo-inverse of the random matrix and also "hi". Instead it crashes without an error message.

Answer 1

When constructing the Eigen::JacobiSVD object, you fail to request that matrices U and V should be computed. By default, these are not computed. Obviously, accessing these matrices if they are not computed will cause a segmentation violation.

See the documentation to the constructor. A second input argument must specify either ComputeFullU | ComputeFullV, or ComputeThinU | ComputeThinV. The thin ones are preferable when computing the pseudo-inverse, as the rest of the matrices are not needed.

---

I would not derive from the JacobiSVD class just to add a method. Instead, I would simply write a free function. This is both easier, and allows you to use only the documented portions of the Eigen API.

I wrote the following MEX-file, which works as intended (using code I already had for this computation). It does the same, but in a slightly different way that avoids writing an explicit loop. Not sure this way of writing it is very clear, but it works.

// Compile with:
//    mex -v test.cpp -I/usr/local/include/eigen3

#include "mex.h"
#include <Eigen/Dense>
#include <Eigen/SVD>
#include <cstdlib>
#include <cmath>
#include <iostream>

Eigen::MatrixXf PseudoInverse(Eigen::MatrixXf matrix) {
   Eigen::JacobiSVD< Eigen::MatrixXf > svd( matrix, Eigen::ComputeThinU | Eigen::ComputeThinV );
   float tolerance = 1.0e-6f * float(std::max(matrix.rows(), matrix.cols())) * svd.singularValues().array().abs()(0);
   return svd.matrixV()
         * (svd.singularValues().array().abs() > tolerance).select(svd.singularValues().array().inverse(), 0).matrix().asDiagonal()
         * svd.matrixU().adjoint();
}

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    Eigen::MatrixXf X = Eigen::MatrixXf::Random(5, 5);
    Eigen::MatrixXf Y = PseudoInverse(X);
    std::cout << Y << '\n';
    std::cout << "hi\n";
}