BLAS矩阵对矩阵转置乘法

Question

我需要计算形如 A'A 或更一般的 A'DA 的一些乘积，其中 A 是一个一般的大小为 mxn 的矩阵，而 D 则是一个对角线大小为 mxm 的矩阵。它们都是满秩的；即 rank(A)=min(m,n)。

我知道在这样的对称乘积中可以节省大量时间：由于 A'A 是对称的，你只需要计算乘积矩阵的下方--或上方--对角线部分。这添加了要计算的项数为 n(n+1)/2 ，这大致相当于大型矩阵的典型 n^2 的一半。

这是我想利用的一个巨大的优势，而我知道我可以在循环中实现矩阵-矩阵乘积。然而，到目前为止我一直在使用 BLAS，因为它比我自己编写的任何循环实现方式都快得多，因为它优化了高速缓存和内存管理。

有没有一种方法可以使用 BLAS 高效地计算 A'A 或甚至 A'DA ？

谢谢！

- enanone

Answer 1

您正在寻找与BLAS的dsyrk子例程相关的信息。

根据文档描述：

SUBROUTINE dsyrk（UPLO，TRANS，N，K，ALPHA，A，LDA，BETA，C，LDC）

DSYRK执行对称秩k运算之一

C := alpha*A*A**T + beta*C,

或

C := alpha*A**T*A + beta*C，

其中alpha和beta是标量，C是一个n乘n的对称矩阵，A在第一种情况下是一个n乘k的矩阵，在第二种情况下是一个k乘n的矩阵。

如果以上三角形式存储A'A，则：

CALL dsyrk( 'U' , 'T' ,  N , M ,  1.0  , A , M , 0.0 , C , N )

对于A'DA，在BLAS中没有直接的等效物。但是您可以在for循环中使用dsyr。

SUBROUTINE dsyr(UPLO,N,ALPHA,X,INCX,A,LDA)

DSYR执行对称秩1操作

A := alpha*x*x**T + A,

其中alpha是实数标量，x是n元素向量，A是n乘以n对称矩阵。

do i = 1, M
    call dsyr('U',N,D(i,i),A(1,i),M,C,N)
end do

- ztik

Answer 2

BLAS中由@ztik建议使用的dsyrk例程是用于A'A。对于A'DA，一种可能性是使用dsyr2k例程，它可以执行对称秩2k运算：

C := alpha*A**T*B + alpha*B**T*A + beta*C.

设置alpha = 0.5, beta = 0.0，并让B = DA。请注意，这种方式假设您的对角矩阵D是实数。

- Qimen Xu

---

由于D是对角矩阵，您也可以对其进行平方根，并将其预乘到A上，然后在新矩阵上使用dsyrk。 - John

Answer 3

对于A'A，可以使用SYRK。对于A'DA，您可以在一侧使用SYMM，例如V = A'D，然后使用英特尔MKL的GEMMT来计算W = VA。GEMMT与GEMM类似，但它利用了结果矩阵是对称的这一事实，因此只需要完成大约一半的工作。

- Luke Mazur