使用gsl，Blas和Lapack计算（Aᵀ×A）*（Bᵀ×B）矩阵

Question

让我有2个对称矩阵：

A = {{1,2}, {2,3}}
B = {{2,3},{3,4}}

我可以使用gsl，Blas和Lapack计算矩阵（AT×A）*（BT×B）吗？

我正在使用

gsl_blas_dsyrk(CblasUpper, CblasTrans, 1.0, A, 0.0, ATA);
gsl_blas_dsyrk(CblasUpper, CblasTrans, 1.0, B, 0.0, BTB);
gsl_blas_dsymm(CblasLeft, CblasUpper, 1.0, ATA, BTB, 0.0, ATABTB); // It doesn't work

它返回：

(Aᵀ·A) = ATA = {{5, 8}, {0, 13}} -- ok, gsl_blas_dsyrk returns symmetric matrix as upper triangular matrix.
(Bᵀ·B) = BTB = {{13, 8}, {0, 25}} -- ok.
(Aᵀ·A)·(Bᵀ·B) = ATABTB = {{65, 290}, {104, 469}} -- it's wrong.

Answer 1

Symmetrize BTB和问题将得到解决。

正如您所注意到的，对称矩阵的上三角形部分由dsyrk()计算。然后应用dsymm()。根据dsymm()的定义，自使用标志CblasLeft以来执行以下操作：

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

其中alpha和beta是标量，A是对称矩阵，B和C是m×n矩阵。

实际上，B矩阵是一般矩阵，不一定是对称矩阵。结果，ATA乘以BTB的上三角部分，因为不计算BTB的下三角部分。

Symmetrize BTB和问题将得到解决。要做到这一点，for循环是一个简单的解决方案，请参阅Convert symmetric matrix between packed and full storage?