c++ - 使用 cblas_dgemm 计算伪逆的问题

Question

我正在尝试使用 Intel MKL 计算存储在 LAPACK_ROW_MAJOR 布局中的矩阵的伪逆。

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=5, ldvt=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.
    }

我的问题是用 u_2x5_needed_for_pinv 执行矩阵乘法 v_4x2_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)。

Answer 1

由于最后三个奇异值为零，我们可以说 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);