RcppArmadillo中协方差矩阵快速秩2更新的优化问询

核心优化：利用对称性减少冗余计算 #

更新项 B %*% H %*% t(B) 是对称矩阵（因H对称，转置后与原矩阵相等），而Sigma本身也是对称矩阵。我们可以利用这一特性，避免计算整个矩阵的所有元素，只计算三角部分再复制到另一侧，直接砍掉一半计算量。

1. 用Armadillo内置对称矩阵类型自动优化 #

Armadillo提供symmat类型，显式声明矩阵为对称后，底层会自动跳过冗余计算：

// 声明Sigma为对称矩阵
arma::symmat<double> Sigma = ...; 
// H是对称矩阵，同样声明为symmat以触发优化
arma::symmat<double> H = ...;
arma::mat B = ...;

// 直接执行原地更新，Armadillo会自动利用对称性优化乘法
Sigma -= B * H * B.t();

这种方式无需手动写循环，依赖库的底层优化即可提升效率。

2. 手动计算三角部分（精细控制场景） #

如果需要更底层的控制，可以手动计算更新矩阵的上三角（或下三角），再复制到对称位置：

arma::mat Sigma = ...; // 假设Sigma本身是对称的
arma::mat B = ...;
arma::mat H = ...; // H为对称矩阵

// 先计算中间项 C = B %*% H，避免重复计算
arma::mat C = B * H;

// 只计算上三角并同步下三角
for (int i = 0; i < Sigma.n_rows; ++i) {
    for (int j = i; j < Sigma.n_cols; ++j) {
        Sigma(i,j) -= arma::dot(C.row(i), B.row(j));
        if (i != j) {
            Sigma(j,i) = Sigma(i,j); // 利用对称性复制，无需重复计算
        }
    }
}

这里用dot计算行向量内积，比生成完整中间矩阵后取元素更节省内存和计算时间。

其他RcppArmadillo加速技巧 #

• 复用内存，减少分配开销：迭代过程中提前声明好中间矩阵（比如作为类成员或全局变量），每次迭代仅重新赋值，避免频繁的内存分配与释放。
• 启用多线程线性代数库：确保R链接到OpenBLAS或MKL这类多线程BLAS/LAPACK库，Armadillo会自动调用其多线程优化。可通过R中sessionInfo()查看当前使用的线性代数库。
• 优先使用原地操作符：用+=、-=这类原地修改运算符代替Sigma = Sigma - ...，避免生成临时矩阵带来的额外开销。
• 避免冗余数据复制：如果B是从Sigma中提取的列，直接用Armadillo的列视图（Sigma.cols(col_indices)），无需复制整列数据。

内容的提问来源于stack exchange，提问作者Søren Højsgaard