当A是下三角矩阵时，实现ATA的更好方法

Question

在特征库中实现A^T*A可以写：

X.template triangularView<Lower>().setZero(); 
X.template selfadjointView<Lower>().rankUpdate(A.transpose());

如果A是下三角矩阵，有没有更好(更有效(的方法来写它？我尝试了以下操作，但它给出了编译错误：

X.template selfadjointView<Lower>().rankUpdate(A.template triangularView<Lower>().transpose());

它给出错误：

error: no matching member function for call to 'rankUpdate'

Answer 1

Eigen通常只包括那些在BLAS中具有相应函数的矩阵例程。对于这种特殊情况，没有BLAS例程。另一个问题是，除非定义-DEIGEN_USE_BLAS并链接到执行此操作的BLAS库(例如OpenBLAS(，否则Eigen的rankUpdate和三角矩阵乘法的实现是不并行的。

我对各种实现进行了一些测试。

rankUpdate

void trtr(Eigen::MatrixXf& out, const Eigen::MatrixXf& in)
{
out = Eigen::MatrixXf::Zero(in.rows(), in.cols());
out.selfadjointView<Eigen::Lower>().rankUpdate(in.transpose());
out.triangularView<Eigen::StrictlyUpper>() = out.transpose();
}

仅与OpenBLAS并行。然后它是相当快，否则有点慢。

三角矩阵乘积

void trtr(Eigen::MatrixXf& out, const Eigen::MatrixXf& in)
{ out.noalias() = in.transpose() * in.triangularView<Eigen::Lower>(); }

仅与OpenBLAS并行。即便如此，它还是有些缓慢。

通用矩阵产品

void trtr(Eigen::MatrixXf& out, const Eigen::MatrixXf& in)
{ out.noalias() = in.transpose() * in; }

始终并行。持续快速。

自定义分块乘法

我尝试实现我自己的版本，将矩阵拆分为块，以删除多余的操作。对于大型矩阵来说，这要快一些，尤其是在使用OpenBLAS时。

它对小矩阵的调整不好。

void trtr_recursive
(Eigen::Ref<Eigen::MatrixXf> out,
const Eigen::Ref<const Eigen::MatrixXf>& in) noexcept
{
Eigen::Index size = out.rows();
if(size <= 16) {
/*
* Two notes:
* 1. The cutoff point is a possible tuning parameter
* 2. This is the only place where we assume that the upper triangle
*    is initialized to zero. We can remove this assumption by making
*    a copy of the input into a fixed size matrix. Use Eigen::Matrix
*    with the MaxRows and MaxCols template argument for this
*/
out.selfadjointView<Eigen::Lower>().rankUpdate(in.transpose());
return;
}
Eigen::Index full = (size + 1) >> 1;
Eigen::Index part = size >> 1;
const auto& fullin = in.bottomLeftCorner(full, full);
const auto& bottomright = in.bottomRightCorner(part, part);
const auto& bottomleft = in.bottomLeftCorner(part, full);
out.topLeftCorner(full, full).selfadjointView<Eigen::Lower>()
.rankUpdate(fullin.transpose());
out.bottomLeftCorner(part, full).noalias() +=
bottomright.transpose().triangularView<Eigen::Upper>() *
bottomleft;
trtr_recursive(out.topLeftCorner(part, part),
in.topLeftCorner(part, part));
trtr_recursive(out.bottomRightCorner(part, part),
bottomright);
}
void trtr(Eigen::MatrixXf& out, const Eigen::MatrixXf& in)
{
out = Eigen::MatrixXf::Zero(in.rows(), in.cols());
trtr4_recursive(out, in);
out.triangularView<Eigen::StrictlyUpper>() = out.transpose();
}