实现特征库伪逆函数崩溃的MEX文件

Question

我试图在Matlab MEX文件中实现一个特征库伪逆函数.它编译成功，但当我运行它时会崩溃。

我正试着跟随利用特征库实现伪逆函数的常见问题。

FAQ建议将其作为一种方法添加到JacobiSVD类中，但是由于您不能在C++中这样做，所以我要将它添加到子类中。它编译成功，但随后在没有错误消息的情况下崩溃。如果我用.pinv调用注释掉行，那么它成功地输出了"hi“而不会崩溃，因此问题就出现在这里。要运行，我只需编译它(作为test.cpp)，然后在命令行输入test。我正在使用Matlab R2019a下的MacOS 10.14.5和特征3.3.7。在我的完整代码中，我还得到了许多关于pinv代码的奇怪错误消息，但在我排除故障之前，我需要这个简单的测试用例来工作。这一切都在我对C++的理解范围内。任何帮助都很感激。

#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;
}

期望的结果是输出随机矩阵的伪逆和"hi“。相反，它在没有错误消息的情况下崩溃。

Answer 1

在构造Eigen::JacobiSVD对象时，无法请求计算矩阵U和V。默认情况下，不会计算这些值。显然，如果不计算这些矩阵，则访问它们将导致分段冲突。

见构造函数的文档。第二个输入参数必须指定ComputeFullU | ComputeFullV或ComputeThinU | ComputeThinV。当计算伪逆时，薄矩阵更可取，因为不需要其余的矩阵。

我不会仅仅为了添加一个方法而从JacobiSVD类派生。相反，我只需要编写一个免费的函数。这两种方法都比较容易，并且允许您只使用特征API的文档部分。

我编写了以下MEX文件，该文件按预期工作(使用我已经为此计算编写的代码)。它也是这样做的，但以一种略为不同的方式避免了编写显式循环。不确定这种书写方式是非常清楚的，但它有效。

// 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";
}