生成随机矩阵时有些奇怪(特征库)
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());
我的意图是在 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());