jeudi 3 septembre 2020

Something strange when generating random matrices ( Eigen Library )

My intention is to stack, in tmp, matrices generated as orth(randn(lenghtDim, dimsSubsp(iK)))' (in Matlab notation) I simulate nresample times this procedure and every time I compute the squared norm of tmp and save it in normVal. I tried different ways to generate these random matrices but their norm value is always the same (with Matlab I don't have this behaviour)!

Could you help me to understand this strange behaviour and to solve it? Thank you

Eigen::VectorXd EmpDistrLargSV(const std::size_t lenghtDim, const std::vector<std::size_t>& dimsSubsp, int nresample ){
 Eigen::VectorXd normVal;
 Eigen::MatrixXd tmp(std::accumulate(dimsSubsp.cbegin(),dimsSubsp.cend(),0),lenghtDim);
 normVal.resize(nresample);
 std::normal_distribution<double> distribution(0,1);
 std::default_random_engine engine (nresample );
    for (int i = 0; i <nresample ; ++i) {
        for (int iK = 0; iK < dimsSubsp.size(); ++iK) {
            std::size_t row_start=std::accumulate(dimsSubsp.begin(),dimsSubsp.begin()+iK,0);
            Eigen::MatrixXd myRandMat = Eigen::MatrixXd::NullaryExpr(lenghtDim,dimsSubsp[iK],[&](){return distribution(engine);});
            Eigen::JacobiSVD<Eigen::MatrixXd> svd(myRandMat, Eigen::ComputeThinU );
            tmp.block(row_start,0,dimsSubsp[iK],lenghtDim)=svd.matrixU().transpose();

        }
        normVal(i)=tmp.squaredNorm();
    }
    return normVal;
}



Aucun commentaire:

Enregistrer un commentaire