BLAS 矩阵乘以矩阵转置

Question

我必须以 A'A 或更一般的 A'DA 形式计算一些产品，其中 A 是一般的 mxn 矩阵，D 是对角 mxm 矩阵。他们都是满级；即rank(A)=min(m,n).

我知道这样的对称乘积可以节省大量时间：鉴于 A'A 是对称的，您只需计算乘积矩阵的下对角线部分或上对角线部分。这增加了要计算的 n(n+1)/2 个条目，这大约是大型矩阵的典型 n^2 的一半。

这是我想要利用的巨大节省，而且我知道我可以在 for 循环中实现矩阵-矩阵乘法。但是，到目前为止，我一直在使用 BLAS，它比我自己编写的任何 for 循环实现都要快得多，因为它优化了缓存和内存管理。

有没有办法使用 BLAS 高效地计算 A'A 甚至 A'DA？谢谢！

Answer 1

您正在寻找 BLAS 的 dsyrk 子程序。

如文档中所述：

SUBROUTINE dsyrk(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)

DSYRK performs one of the symmetric rank k operations

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

or

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

where alpha and beta are scalars, C is an n by n symmetric matrix and A is an n by k matrix in the first case and a k by n matrix in the second case.

在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 performs the symmetric rank 1 operation

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

where alpha is a real scalar, x is an n element vector and A is an n by n symmetric matrix.

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

Answer 2

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

Answer 3

@ztik 建议的 BLAS 中的 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 是真实的。

BLAS 矩阵乘以矩阵转置

BLAS matrix by matrix transpose multiply

matrix

linear-algebra

blas