生成随机矩阵时有些奇怪(特征库)

Something strange when generating random matrices ( Eigen Library )

我的意图是在 tmp 中堆叠生成的矩阵 orth(randn(lenghtDim, dimsSubsp(iK)))' (在 Matlab 符号中) 我模拟 nresample 次此过程,每次我计算 tmp 的平方范数并将其保存在 normVal 中。我尝试了不同的方法来生成这些随机矩阵,但它们的范数值始终相同(使用 Matlab 我没有这种行为)!

你能帮我理解这个奇怪的行为并解决它吗?谢谢

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

--- 编辑 ---

我试图用 C++ 编写的是以下 Matlab 代码

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

为了得到 orth() ,在 C++ 中我计算 svd 然后取矩阵 U。

根据 default_random_engine documentation its de fault constructor creates linear congruential engine 从一些种子开始进行简单的递增和取模。

因此您使用的随机引擎是确定性的。

这就是为什么您会得到相同的矩阵和相同的范数。 Matlab 可能是不确定的。

编辑:

哈哈,我知道是什么了。您正在执行 SVD 并查看 U,它始终(无论输入如何)都是酉矩阵。酉矩阵的 属性 即 U^T * U == I,这意味着其每一列的范数(和平方范数)恰好为 1。因此,矩阵的平方范数将等于列数(在您的“瘦”U 情况下,行数或列数的最小值),与您使用的随机数生成器无关。

以下无关信息:

尝试 std::mt19937 而不是 std::default_random_engine。我不确定使用 nresample 作为种子对你来说是否重要,但你可能想尝试使用 std::chrono::high_resolution_clock::now().time_since_epoch().count() 之类的东西。否则,我认为您的方法与我的相似。

#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());