C++:Eigen conservativeResize 太贵了？

Question

我有一些我事先不知道其维度的特征矩阵，我只有一个上限。我有一个循环，其中我逐列填充这些矩阵(我使用上限初始化它们)直到满足停止条件(假设在 j 次迭代之后)。

我现在的问题是:在循环之后，我需要那些矩阵进行矩阵乘法(显然只使用前 j 列)。直接的解决方案是使用 Eigen 的 conservativeResize 并继续执行矩阵乘法。由于矩阵往往非常大(超过 100000 维)并且(据我所知，虽然不确定)Eigen 的 conservativeResize 为调整大小的矩阵重新分配内存并执行一个深拷贝，这个解决方案非常昂贵。

我正在考虑编写自己的自定义矩阵乘法函数，它使用旧的(大)矩阵，采用指定要使用的列数的参数。我担心 Eigen 的矩阵乘法优化得太多了，最终这个解决方案比仅使用保守的大小调整和标准 Eigen 乘法要慢......

我应该硬着头皮使用 conservativeResize 还是有人有更好的主意？顺便说一句:我们正在讨论的矩阵在循环/调整大小后用于 3 次乘法和 1 次转置

提前致谢!

编辑:这是代码的相关部分(其中 X 是 MatrixXd，y 是 VectorXd，numComponents 是 PLS1 应该使用的潜在变量的数量)。问题是:一开始，numComponents 将始终是 X (X.cols()) 中的维数，但停止标准应该检查输出 vector 中解释方差的相对改进(即，我有尚未实现)。如果相对改进太小，算法应该停止(因为我们对前 j 个分量很满意)然后计算回归系数。为此，我需要 conservativeResize:

using namespace Eigen;
MatrixXd W,P,T,B;
VectorXd c,xMean;
double xMean;

W.resize(X.cols(),numComponents);
P.resize(X.cols(),numComponents);
T.resize(X.rows(),numComponents);
c.resize(numComponents);
xMean.resize(X.cols());
xMean.setZero();
yMean=0;
VectorXd yCopy=y;
//perform PLS1
for(size_t j=0; j< numComponents; ++j){
    VectorXd tmp=X.transpose()*y;
    W.col(j)=(tmp)/tmp.norm();
    T.col(j)=X*W.col(j);
    double divisorTmp=T.col(j).transpose()*T.col(j);
    c(j)=(T.col(j).transpose()*y);
    c(j)/=divisorTmp;
    P.col(j)=X.transpose()*T.col(j)/divisorTmp;
    X=X-T.col(j)*P.col(j).transpose();
    y=y-T.col(j)*c(j);
    if(/*STOPPINGCRITERION(TODO)*/ && j<numComponents-1){
        numComponents=j+1;
        W.conservativeResize(X.cols(),numComponents);
        P.conservativeResize(X.cols(),numComponents);
        T.conservativeResize(X.rows(),numComponents);
        c.conservativeResize(numComponents);
    }
}
//store regression matrix
MatrixXd tmp=P.transpose()*W;
B=W*tmp.inverse()*c;
yCopy=yCopy-T*c;
mse=(yCopy.transpose()*yCopy);
mse/=y.size();//Mean Square Error

Answer 1

我想你可以分配一次大矩阵，然后乘法使用 block 创建它的部分的 View ，其中将包含有意义的数据。然后你可以重复使用一个大矩阵。这将使您免于重新分配。

下面的例子充分说明了这一点。

./eigen_block_multiply.cpp:

#include <Eigen/Dense>
#include <iostream>
using namespace std;
using namespace Eigen;
int main()
{
  Matrix<float, 2, 3> small;
  small << 1,2,3,
           4,5,6;

  Matrix<float, 4, 4> big = Matrix<float, 4, 4>::Constant(0.6);
  cout << "Big matrix:\n";
  cout << big << endl;
  cout << "Block of big matrix:\n";
  cout << big.block(0,0,3,2) << endl;
  cout << "Small matrix:\n";
  cout << small << endl;
  cout << "Product:\n";
  cout << small * big.block(0,0,3,2) << endl;

  Matrix<float, 3, 3> small2;
  small2 << 1,2,3,
            4,5,6,
            7,8,9;
  big = Matrix<float, 4, 4>::Constant(6.66);
  cout << "Product2:\n";
  cout << small * big.block(0,0,3,3) << endl;
}

输出:

Big matrix:
0.6 0.6 0.6 0.6
0.6 0.6 0.6 0.6
0.6 0.6 0.6 0.6
0.6 0.6 0.6 0.6
Block of big matrix:
0.6 0.6
0.6 0.6
0.6 0.6
Small matrix:
1 2 3
4 5 6
Product:
3.6 3.6
  9   9
Product2:
39.96 39.96 39.96
 99.9  99.9  99.9