r - 从计算的角度来看，有没有办法简化 $ADA^T$，其中 $D$ 是对角矩阵？

Question

设 $A$ 为 $m\times n$ 矩阵，不一定对称，$D$ 为对角 $n\times n$ 矩阵。在我的问题中，$m < n$ 但 $m$ 和 $n$ 都有可能非常大，因此当维度变得非常高时，$ADA^T$ 的计算可能需要几秒钟。我正在编写的程序将针对 $A$ 和 $D$ 的不同变体多次执行此计算，因此我想找到一种加快计算速度的方法，因为这可以在脚本过程中节省几分钟。

目前，我只是在执行tcrossprod(A %*% D, A)。但我觉得应该有一种更快的方法来计算它，因为 $D$ 是对角矩阵。从代数或计算角度来看，有人在这里看到任何好的捷径吗？

Answer 1

只是为了枚举底层 C 代码所采用的不同路径tcrossprod ...

对于密集A，隐式类矩阵:

如果D是非负数，那么您可以获得BLAS例程DSYRK，它利用对称性，具有:

tcrossprod(A * rep(sqrt(diag(D, names = FALSE)), each = nrow(A)))

否则，您将陷入 BLAS 例程DGEMM:

tcrossprod(A * rep(diag(D, names = FALSE), each = nrow(A)), A)

您可以尝试使用 external BLAS如果您知道 A 和 D 是有限的，不包含 IEEE 特殊值 Inf 或 NaN(包括 R 的 NA_real_)。如果 R 检测到任一矩阵因子非有限，则不会使用外部 BLAS。

对于稀疏A，正式类dgCMatrix:

如果D是非负数，那么您可以获得CHOLMOD例程cholmod_aat，它利用对称性，具有:

A@x <- A@x * rep.int(sqrt(diag(D, names = FALSE)), A@p[-1L] - A@p[-ncol(A)])
tcrossprod(A)

否则，您将陷入 CHOLMOD 例程 cholmod_transpose 和 cholmod_ssmult:

AD <- A
AD@x <- A@x * rep.int(diag(D, names = FALSE), A@p[-1L] - A@p[-ncol(A)])
tcrossprod(AD, A)

无论如何，由于 D 是对角线，因此您应该只需要一次矩阵乘法。