MEX-файл, реализующий сбой псевдообратной функции библиотеки Eigen

Question

Я пытаюсь реализовать псевдообратную функцию библиотеки Eigen в MEX-файле Matlab.Он успешно компилируется, но вылетает при его запуске.

Я пытаюсь следовать часто задаваемым вопросам о том, как реализовать псевдообратную функцию с использованием библиотеки Eigen .

.FAQ предлагает добавить его как метод к классу JacobiSVD, но так как вы не можете сделать это в C ++, я добавляю его в дочерний класс.Он успешно компилируется, но затем падает без сообщения об ошибке.Он успешно выдает «привет» без сбоев, если я закомментирую строку с вызовом .pinv, так что именно здесь возникает проблема.Для запуска я просто скомпилирую его (как test.cpp), а затем введите test в командной строке.Я использую Matlab R2019a под MacOS 10.14.5 и Eigen 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 только для добавления метода.Вместо этого я просто написал бы бесплатную функцию.Это и проще, и позволяет вам использовать только документированные части Eigen 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";
}