Mon intention est d'empiler, en tmp , les matrices générées comme orth (randn (lenghtDim, dimsSubsp (iK))) ' (en notation Matlab) Je simule nresample fois cette procédure et chaque fois que je calcule la norme au carré de tmp et la sauvegarde dans normVal. J'ai essayé différentes manières de générer ces matrices aléatoires mais leur valeur de norme est toujours la même (avec Matlab je n'ai pas ce comportement)!

Pourriez-vous m'aider à comprendre ce comportement étrange et à le résoudre? Merci

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

--- ÉDITER ---

Ce que j'essaie d'écrire en C ++ est le code Matlab suivant

nb = length(dimsSubsp);
tmp = zeros(sum(dimsSubsp), lengthDim);

normVal = zeros(1, nresample);

    for i = 1:nresample
        for ib = 1:nb
            irow_start = sum(dimsSubsp(1 : (ib - 1))) + 1;
            irow_end = sum(dimsSubsp(1 : ib));
            tmp(irow_start : irow_end, :) =  orth(randn(lengthDim,dimsSubsp(ib)))';
        end
        normVal(i) = norm(M, 2)^2;
    end

Pour obtenir le orth(), en C ++ je calcule le svd puis je prends le matrixU.

0
Virginie 3 sept. 2020 à 16:33

2 réponses

Meilleure réponse

ÉDITER:

Haha, je sais ce que c'est maintenant. Vous effectuez un SVD et regardez U, qui est toujours (quelle que soit l'entrée) une matrice unitaire. Une matrice unitaire a la propriété U^T * U == I, ce qui implique que la norme (et la norme au carré) de chacune de ses colonnes sont exactement 1. Par conséquent, la norme au carré de la matrice va être égale au nombre de colonnes (le minimum de lignes ou de colonnes dans le cas de votre "thin" U), quel que soit le générateur de nombres aléatoires que vous utilisez.

Informations non pertinentes ci-dessous:

Au lieu de std::default_random_engine, essayez std::mt19937. Je ne sais pas s'il est important pour vous d'utiliser nresample comme graine, mais vous pouvez essayer d'utiliser quelque chose comme std::chrono::high_resolution_clock::now().time_since_epoch().count(). Sinon, je pense que votre approche est similaire à la mienne.

#include <chrono>
#include <random> 
#include <Eigen/eigen>

MatrixX<double> random_matrix(int rows, int cols, double min, double max, unsigned seed)
{
   std::mt19937 generator(seed);
   std::uniform_real_distribution<double> distribution(min,max);
   MatrixX<double> result(rows,cols);
   for(double& val : result.reshaped())
      val = distribution(generator);
   return result;
}

MatrixX<double> rando = random_matrix(20, 34, -30.1, 122.3, std::chrono::high_resolution_clock::now().time_since_epoch().count());
1
Charlie S 4 sept. 2020 à 14:35

Selon la documentation default_random_engine, son constructeur par défaut crée moteur congruentiel linéaire qui fait un simple incrément et modulo à partir d'une graine.

Par conséquent, le moteur aléatoire que vous utilisez est déterministe.

C'est pourquoi vous obtenez les mêmes matrices et la même norme. Matlab est probablement indéterministe.

0
Łukasz Ślusarczyk 3 sept. 2020 à 13:45