C++/Eigen で特異値分解(SVD)を使って一般化逆行列を計算する

一般化逆行列（擬似逆行列，最小二乗法に対応するやつ）をc++で計算するサンプル

特異値分解(SVD)が必要なので，実装は Eigen を使います

python/numpy の実装も用意しました．こちらです https://pyopyopyo.hatenablog.com/entry/2021/09/14/161955

#include <Eigen/Dense>
#include <iostream>
#include <random>
#include <assert.h>

int
main()
{

    //ランダムな 3x3 行列を生成

    std::random_device rd;
    std::mt19937 gen(rd());
    std::normal_distribution<float> dis(0, 10);
    //ラムダ式で乱数を生成，生成した値を行列の各要素に格納する
    Eigen::MatrixXf A = Eigen::MatrixXf::NullaryExpr(3,3,  
                              [&](){return dis(gen);});

    // std::cout << "A: " << A << std::endl;

    Eigen::JacobiSVD<Eigen::MatrixXf> svd(A,
					  Eigen::ComputeFullU
					  | Eigen::ComputeFullV);

    //std::cout << "U: " << svd.matrixU() << std::endl;
    //std::cout << "S: " << svd.singularValues()  << std::endl;
    //std::cout << "V: " << svd.matrixV() << std::endl;


    // 一般化逆行列を計算する
    //
    //
    //  svdで  行列A を U * S * V.T と分解しているから
    //
    //  一般化逆行列 A+  は V * S.T * U.T

    Eigen::VectorXf s = svd.singularValues();
    s = s.array().inverse();

    // 一般化逆行列の変数名は Ainv 
    Eigen::MatrixXf Ainv =
                svd.matrixV()
             * s.asDiagonal()  
             * svd.matrixU().transpose();

    // A と A+ をかけると単位行列になる
    Eigen::MatrixXf m= A * Ainv;
    std::cout << m << std::endl;

    // 実際は，計算誤差があるので「ほぼ」単位行列になる．
    // そこで mとIdentity(3,3)が「ほぼ」等しいかEigenの isApprox()　で確認
    // 結果は assert()でチェック

    assert( Eigen::MatrixXf::Identity(3,3).isApprox(m) );

    return 0;