C++ SVDによる低ランク近似

O(n³)で入力行列Pのランクをrにするアルゴリズムです．eigenをつかっています．

#include <iostream>
#include <Eigen/Dense>
#include <vector>
#include <cmath>
#include <time.h>
#include <cstdio>
#include <fstream>
#include "utils.h"

using namespace std;
using namespace Eigen;

MatrixXd SVD(MatrixXd& P, int r){
    int const N = P.rows();
    int const rank = r;

    //JacobiSVD<MatrixXd> svd(P, ComputeFullU | ComputeFullV);
    BDCSVD<MatrixXd> svd(P, ComputeFullU | ComputeFullV);
    MatrixXd SVD_P(N,N);
    MatrixXd U(N,N), Ur(N,rank-1);
    MatrixXd V(N,N), Vr(N,rank-1);
    MatrixXd S(N,N), Sr(rank-1, rank-1);
    VectorXd sigma(N);

    U = svd.matrixU();
    V = svd.matrixV();
    sigma = svd.singularValues();
    S = sigma.array().matrix().asDiagonal();

    Ur = U(all, seq(0,rank-1));
    Vr = V(all, seq(0,rank-1));
    Sr = S(seq(0,rank-1),seq(0,rank-1));

    SVD_P = Ur * Sr * Vr.transpose();

    return SVD_P;
}