从计算的Angular 来看，有没有什么方法可以简化$ADA^T$，其中$D$是一个对角矩阵？

7gyucuyw 于 2023-05-04 发布在 Angular

关注(0)|答案(2)|浏览(133)

设$A$是一个$m \times n$矩阵，不一定对称，$D$是一个对角$n \times n$矩阵。在我的问题中，$m $〈$n$，但是$m$和$n$都有可能非常大，所以$ADA^T$的计算可能需要几秒钟，因为维度非常高。我正在编写的程序将对$A$和$D$的不同变体执行多次这种计算，因此我想找到一种方法来加速它，因为这可以在脚本过程中保存几分钟。
目前，我正在做tcrossprod(A %*% D, A)。但是我觉得应该有一种更快的方法来计算它，因为$D$是一个对角矩阵。从代数或计算的Angular 来看，有没有人在这里看到任何好的捷径？

来源：https://stackoverflow.com/questions/76149339/is-there-any-way-to-simplify-adat-from-a-computational-standpoint-where-d

2条答案

按热度按时间

iszxjhcz1#

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

对于密集`A`，隐式类`matrix`：

如果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)

如果您知道A和D是有限的，不包含IEEE特殊值Inf或NaN（包括R的NA_real_），则可以使用外部BLAS进行实验。如果R检测到任一矩阵因子是非有限的，则R将不使用外部BLAS。

对于形式类`dgCMatrix`的稀疏`A`：

如果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是对角的，你应该只需要一个矩阵相乘。

赞(0）回复(0）举报 2023-05-04

lxkprmvk2#

以下是几种可能实现的基准测试示例，其中D的对角线项可能为负（因此as.complex用于f4和f5）

set.seed(1)
m <- 100
n <- 500
A <- matrix(runif(m * n), m, n)
D <- diag(rnorm(n))

microbenchmark(
  f0 = A %*% D %*% t(A),
  f1 = A %*% (diag(D) * t(A)),
  f2 = tcrossprod(A %*% D, A),
  f3 = tcrossprod(A, tcrossprod(A, D)),
  f4 = Re(tcrossprod(A %*% `dim<-`(sqrt(as.complex(D)), dim(D)))),
  f5 = Re(tcrossprod(A * sqrt(as.complex(diag(D)))[col(A)])),
  check = "equal"
)

等着瞧吧

Unit: milliseconds
 expr     min       lq      mean   median       uq      max neval
   f0 14.2376 15.10425 15.411022 15.30060 15.48965  18.6084   100
   f1  2.5351  2.69190  2.837047  2.73960  2.80545   6.5885   100
   f2 14.1424 14.80185 15.114680 14.95035 15.10380  20.3685   100
   f3 13.5691 14.38270 15.200161 14.82995 15.45375  23.1657   100
   f4 32.7007 34.29920 37.068124 34.99865 36.99130 118.8519   100
   f5  5.4284  5.61460  5.989452  5.83605  5.92440  10.4077   100

赞(0）回复(0）举报 2023-05-04

我来回答

从计算的Angular 来看，有没有什么方法可以简化$ADA^T$，其中$D$是一个对角矩阵？

2条答案

对于密集`A`，隐式类`matrix`：

对于形式类`dgCMatrix`的稀疏`A`：

相关问题

热门标签

最新问答

从计算的Angular 来看，有没有什么方法可以简化$ADA^T$，其中$D$是一个对角矩阵？

2条答案

对于密集A，隐式类matrix：

对于形式类dgCMatrix的稀疏A：

相关问题

热门标签

最新问答

对于密集`A`，隐式类`matrix`：

对于形式类`dgCMatrix`的稀疏`A`：