Проблема с использованием cblas_dgemm для вычисления псевдообратного - PullRequest
Новая
       46

Проблема с использованием cblas_dgemm для вычисления псевдообратного

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

Я пытаюсь вычислить псевдообратную матрицу, сохраненную в макете LAPACK_ROW_MAJOR, используя Intel MKL.

A_5x4 =
     1     2     3     4
     5     6     7     8
     9    10    11    12
    13    14    15    16
    17    18    19    20

Я использую функцию gesvd для вычисления компактной формы SVD: 

info = LAPACKE_dgesvd(LAPACK_ROW_MAJOR, 'S', 'S', m, n, A, lda, s, u, ldu, vt, ldvt, superb);

, где m=5, n=4, lda=4, ldu=5ldvt=4.Я могу успешно использовать функцию MKL для получения SVD матрицы, A = U*S*VT: 

u_5x4 = 
     0.0965          0.7686          0.6323          0.0034
     0.2455          0.4896         -0.6208          0.0412
     0.3945          0.2107         -0.3285         -0.4681
     0.5435         -0.0683         -0.0097          0.7989
     0.6924         -0.3472          0.3267         -0.3754

s_4x1 = 
    53.520222
     2.363426
     0.000000
     0.000000

vt_4x4 = 
     0.4430          0.4799          0.5167          0.5536
    -0.7097         -0.2640          0.1816          0.6273
     0.0912         -0.5242          0.7747         -0.3417
     0.5401         -0.6521         -0.3160          0.4280

Поскольку s имеет только два ненулевых элемента, мне нужно рассмотреть первые два столбца u и два столбца v (не vt), а также инверсия элементов s

v_4x2_needed_for_pinv = 
     0.4430    0.4799
    -0.7097   -0.2640
     0.0912   -0.5242
     0.5401   -0.6521

u_2x5_needed_for_pinv = 
     0.0965   0.2455   0.3945   0.5435   0.6924
     0.7686   0.4896   0.2107  -0.0683  -0.3472

Я могу выполнить умножение матриц с for-loop без проблем ивычислить псевдообратное значение A. Однако мне очень интересно использовать dscal и cblas_dgemm главным образом потому, что фактическая матрица, обратная величина которой будет вычислена, очень велика.

Мне удалось успешно выяснитьиспользуя dscal и умножьте первые два столбца V на обратное к S: 

MKL_INT  k = ((m) < (n) ? (m) : (n));
// Computing VT = vt*(s^-1)
MKL_INT incx = 1;
MKL_INT r = 0;
for (int i = 0; i < k; i++)
{
    double ss;
    if (s[i] > 1.0e-9)
    {
        ss = 1.0 / s[i];
        r++;
    }
    else
        ss = s[i];
    dscal(&n, &ss, &vt[i*n], &incx); // this replaces vt with new values.
    }

Моя задача - выполнить умножение матриц v_4x2_needed_for_pinv на u_2x5_needed_for_pinv, которые являются подмножеством u и vt массивы, вычисленные LAPACKE_dgesvd.Может ли кто-нибудь помочь мне разобраться, как использовать cblas_dgemm?Я был бы признателен.

Я попробовал следующее, вход в функцию имеет для меня смысл, но он не работает 

// inv(A) = VT^T * U^T = V * U^T
double* inva = (double*)malloc(n*m * sizeof(double));
double alpha = 1.0, beta = 0.0;
MKL_INT ld_inva = n;
cblas_dgemm(CblasRowMajor, CblasTrans, CblasTrans, n, m, r, alpha, vt, n, u, m, beta, inva, ld_inva);

, где r=2, потому что s имеет дватолько ненулевые элементы (53.520222 и 2.363426).

1 Ответ

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

Поскольку последние три единичных значения равны нулю, мы можем сказать, что SVD выдает:

  • u(5,2) с ldu=4
  • vt(2,4) с ldvt=4
  • invA(5,4)

Инверсия вычисляется как invA = vt ^ T * invS * u ^ T и следующий цикл может быть преобразован в invA = (invS * vt) ^ T* u ^ T 

MKL_INT ma = mu = 5;
MKL_INT na = nvt = 4;
MKL_INT nu = mvt = ms = 2;

MKL_INT lda = ldu = ldvt = 4;

// vt = (invS * vt)
for(MKL_INT i=0; i<ms; i++){
  cblas_dscal (nvt, s[i], vt+(i*ldvt), 1);
}

// invA = vt^T * u^T
cblas_dgemm (CblasRowMajor, CblasTrans, CblasTrans, ma, na, nu, 1.0, vt, ldvt, u, ldu, 0.0, invA, lda);
...