将转置矩阵或原始矩阵保留在同一变量中而不使用内存

Question

我想编写以下函数：

Eigen::MatrixXd func(const Eigen::MatrixXd &matrix1, bool condition1,
                     const Eigen::MatrixXd &matrix2, bool condition2,
                     const Eigen::MatrixXd &matrix3, bool condition3) {
    const Eigen::MatrixXd &m1 = (condition1) ? matrix1 : matrix1.transpose();
    const Eigen::MatrixXd &m2 = (condition2) ? matrix2 : matrix2.transpose();
    const Eigen::MatrixXd &m3 = (condition3) ? matrix3 : matrix3.transpose();
    return m1 + m2 + m3;
}

不幸的是，它有时表现出缓慢的性能，10000x10000 矩阵示例：

0 transpose M1+M2+M3:    2.641621 seconds
1 transpose M1+M2+M3T:   4.142240 seconds
1 transpose M1+M2T+M3:   4.335527 seconds
2 transpose M1+M2T+M3T:  5.276609 seconds
1 transpose M1T+M2+M3:   3.844055 seconds
2 transpose M1T+M2+M3T:  6.080098 seconds
2 transpose M1T+M2T+M3:  5.448039 seconds
3 transpose M1T+M2T+M3T: 6.644677 seconds

问题在于

const Eigen::MatrixXd &m = (condition) ? matrix : matrix.transpose();

行，它创建了一个新的转置矩阵，导致额外的计算开销和内存使用。

使用

Eigen::Transpose<Eigen::MatrixXd> m = matrix.transpose()

可以避免不必要的内存分配，但它不能用于存储转置矩阵和原始矩阵。

我的问题是：如何优化上面的函数以避免矩阵复制而不诉诸大量 if-else 语句？

谢谢你。

编辑：我向 Eigen 开发人员 Discord 提出了同样的问题，他们回答说，在这种情况下使用 if-else 分支是最好的方法。

Answer 1

问题是您正在强制转换为

const Eigen::MatrixXd &

。 Transpose 返回

Eigen::Transpose< MatrixType >

代理对象。避免进行转换应该可以防止创建显式副本。

现在将三元运算符合并到带有加法的表达式中可能不会有帮助，因为三元运算符必须选择公共类型。这将强制转换。

因此，为了避免隐式转换，我会尝试使用这种代码：

void addMaybeTransposed(Eigen::MatrixXd &dest, const Eigen::MatrixXd &src, bool transpose)
{
    if (transpose) {
        dest += src.transpose();
    } else {
        dest += src;
    }
}

Eigen::MatrixXd func(const Eigen::MatrixXd &matrix1, bool condition1,
                     const Eigen::MatrixXd &matrix2, bool condition2,
                     const Eigen::MatrixXd &matrix3, bool condition3) {
    Eigen::MatrixXd result{ matrix1.rows(), matrix1.cols() };
    addMaybeTransposed(result, matrix1, condition1);
    addMaybeTransposed(result, matrix2, condition2);
    addMaybeTransposed(result, matrix3, condition3);
    return result;
}

Answer 2

您可以通过模板函数的递归来完成它。通用 lambda 使这变得更容易一些。

Eigen::MatrixXd func(
    const Eigen::MatrixXd &matrix1, bool condition1,
    const Eigen::MatrixXd &matrix2, bool condition2,
    const Eigen::MatrixXd &matrix3, bool condition3)
{
    auto add3 = [&](const auto& left) -> Eigen::MatrixXd {
        return condition3 ?
            (left + matrix3.transpose()).eval() :
            (left + matrix3).eval();
    };
    auto add2 = [&](const auto& left) -> Eigen::MatrixXd {
        return condition2 ?
            add3(left + matrix2.transpose()) :
            add3(left + matrix2);
    };
    return condition1 ? add2(matrix1.transpose()) : add2(matrix1);
}

请注意，这将导致严重的代码膨胀，因为每种可能的组合都将被实例化，但它将导致正确优化的代码路径。

理论上

std::visit

应该能够处理这个问题，但我无法得到合理的编译解决方案。