MKL庫奇異值分解(LAPACKE_dgesvd)

• 2022 年 4 月 24 日
• 筆記
• 高性能計算MKL

對任意一個\(m\times n\)的實矩陣，總可以按照SVD算法對其進行分解。即：

其中\(U、V\)分別為\(m\times m、n\times n\)的方陣，由\(A\)的左奇異向量和右奇異向量組成，且\(U\)與\(V\)均為正交陣。\(\Sigma\)為\(m\times n\)的對角矩陣，對角線上的元素為矩陣\(A\)的奇異值。

在MKL庫中求解奇異值和奇異向量的函數為LAPACKE_dgesvd。

1 參數詳解

lapack_int LAPACKE_dgesvd(
                            matrix_layout,	// (input)行優先(LAPACK_ROW_MAJOR)或列優先(LAPACK_COL_MAJOR)
                            jobu,		// (input)計算矩陣U的全部或部分並返回。
                                                /*"A":返回U的所有M列到U,
                                                  "S":返回U的前min(m,n)列到U,
                                                  "O":返回U的前min(m,n)列到A矩陣(覆蓋),
                                                  "N":不計算矩陣U*/
                            jobvt,		// (input)計算矩陣VT的全部或部分並返回；選項列表與jobu相同；
                            m, 			// (input)A矩陣的行，m>=0
                            n,  		// (input)A矩陣的列，n>=0
                            a, 			// (input/output)A矩陣
                            lda,		// (input)A矩陣的第一維大小
                            s,			// (output)A矩陣的奇異值，並按照從大到小的順序排列
                            u,			// (output) 矩陣U元素的一維數組
                            ldu, 		// (input) U矩陣的第一維大小
                            vt,			// (output) 矩陣VT元素的一維數組
                            ldvt,		// (input) VT矩陣的第一維大小
                            superb,		// (output)工作空間
                        )

2 定義待處理矩陣

#include <stdlib.h>
#include <stdio.h>
#include "mkl_lapacke.h"

#define min(a,b) ((a)>(b)?(b):(a))

// 矩陣維度參數
#define M 6
#define N 5
#define LDA N
#define LDU M
#define LDVT N

// 聲明需要的參數
MKL_INT m = M, n = N, lda = LDA, ldu = LDU, ldvt = LDVT, info;
double superb[min(M,N)-1];				

double s[N], u[LDU*M], vt[LDVT*N];		//聲明奇異值與奇異向量
double a[LDA*M] = {				//定義待分解的A矩陣
    8.79,  9.93,  9.83, 5.45,  3.16,
    6.11,  6.91,  5.04, -0.27,  7.98,
    -9.15, -7.93,  4.86, 4.85,  3.01,
    9.57,  1.64,  8.83, 0.74,  5.80,
    -3.49,  4.02,  9.80, 10.00,  4.27,
    9.84,  0.15, -8.99, -6.02, -5.31
};

3 執行SVD分解

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

結果如圖：

完整代碼

#include <stdlib.h>
#include <stdio.h>
#include "mkl_lapacke.h"

#define min(a,b) ((a)>(b)?(b):(a))

// 展示奇異向量
extern void print_matrix(const char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda);

#define M 6
#define N 5
#define LDA N
#define LDU M
#define LDVT N


int main() {
	//聲明、定義輸入
    MKL_INT m = M, n = N, lda = LDA, ldu = LDU, ldvt = LDVT, info;
    double superb[min(M, N) - 1];

    double s[N], u[LDU * M], vt[LDVT * N];
    double a[LDA * M] = {
        8.79,  9.93,  9.83, 5.45,  3.16,
        6.11,  6.91,  5.04, -0.27,  7.98,
       -9.15, -7.93,  4.86, 4.85,  3.01,
        9.57,  1.64,  8.83, 0.74,  5.80,
       -3.49,  4.02,  9.80, 10.00,  4.27,
        9.84,  0.15, -8.99, -6.02, -5.31
    };

    printf("LAPACKE_dgesvd (row-major, high-level) Example Program Results\n");
    //計算SVD
    info = LAPACKE_dgesvd(LAPACK_ROW_MAJOR, 'A', 'A', m, n, a, lda,
        s, u, ldu, vt, ldvt, superb);

    if (info > 0) {
        printf("The algorithm computing SVD failed to converge.\n");
        exit(1);
    }
    //奇異值
    print_matrix("Singular values", 1, n, s, 1);
    //左奇異向量
    print_matrix("Left singular vectors (stored columnwise)", m, n, u, ldu);
    //右奇異向量
    print_matrix("Right singular vectors (stored rowwise)", n, n, vt, ldvt);
    exit(0);
} 

void print_matrix(const char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda) {
    MKL_INT i, j;
    printf("\n %s\n", desc);
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) printf(" %6.2f", a[i * lda + j]);
        printf("\n");
    }
}

補充：SVD分解求逆

由之前的介紹，對於任意的實數矩陣\(A\)，可以進行SVD分解：

其中，\(U\)、\(V^T\)為正交矩陣，\(\Sigma\)為對角矩陣。若\(A\)矩陣可逆，易得

即當使用LAPACKE_dgesvd，將矩陣\(A\)分解出三部分後，再經過簡單的轉置、對角陣求逆，最後通過LAPACKE_dgemm完成各矩陣相乘即可得到\(A\)的逆矩陣。