Вычислить (Aᵀ × A) * (Bᵀ × B) матрицу, используя gsl, Blas и Lapack - PullRequest
Новая
       102

Вычислить (Aᵀ × A) * (Bᵀ × B) матрицу, используя gsl, Blas и Lapack

0 голосов
/ 09 сентября 2018

Пусть у меня есть 2 симметричные матрицы: 

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

Могу ли я вычислить матрицу (A T × A) * (B T × B), используя gsl, Blas и Lapack?

Я использую 

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.

1 Ответ

0 голосов
/ 14 сентября 2018

Симметризируйте BTB, и проблема будет решена.

Как вы заметили, верхние треугольные части симметричных матриц вычисляются как dsyrk(). Затем dsymm() применяется. Согласно определению dsymm() выполняется следующая операция, поскольку используется флаг CblasLeft: 

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

где альфа и бета - скаляры, A - симметричная матрица, а B и C это m на n матриц.

Действительно, матрица B является общей матрицей , необязательно симметричной. В результате ATA умножается на верхнюю треугольную часть BTB , поскольку нижняя треугольная часть BTB не вычисляется.

Симметризируйте BTB, и проблема будет решена. Для этого для циклов это простое решение, см. Преобразование симметричной матрицы между упакованным и полным хранилищем?

...