0

所以我正在尝试计算NxN矩阵的 SVD。奇怪的是，对于2x2矩阵的所有情况，来自 lapack 和 scipy 的 SVD 都匹配，但是当我选择 a3x3或4x4矩阵时它们会有所不同。

// LAPACK (C) Case: 2x2

#include <stdlib.h>
#include <stdio.h>
#include <Accelerate/Accelerate.h>

/* Complex datatype */
struct _dcomplex { double re, im; };
typedef struct _dcomplex dcomplex;

/* ZGESVD prototype */
extern void zgesvd( char* jobu, char* jobvt, int* m, int* n, dcomplex* a,
                   int* lda, double* s, dcomplex* u, int* ldu, dcomplex* vt, int* ldvt,
                   dcomplex* work, int* lwork, double* rwork, int* info );
/* Auxiliary routines prototypes */
extern void print_matrix( char* desc, int m, int n, dcomplex* a, int lda );
extern void print_rmatrix( char* desc, int m, int n, double* a, int lda );

/* Parameters */
#define M 2
#define N 2
#define LDA M
#define LDU M
#define LDVT N

/* Main program */
int main() {
    /* Locals */
    int m = M, n = N, lda = LDA, ldu = LDU, ldvt = LDVT, info, lwork;
    dcomplex wkopt;
    dcomplex* work;
    /* Local arrays */
    /* rwork dimension should be at least max( 1, 5*min(m,n) ) */
    double s[M], rwork[5*M];
    dcomplex u[LDU*M], vt[LDVT*N];
//    dcomplex a[LDA*N] = {
//        {0, 0}, {0, 0}, {1,  0},
//        {-0.36599657,  -0.27449743}, {-0.27449743,  0.36599657}, {0.76249285, 0},
//        {-0.36599657, 0.27449743}, {-0.27449743, -0.36599657}, {0.76249285, 0},
//    };
    dcomplex a[LDA*N] = {
        {0.70710678, 0}, {0, -0.70710678},
        {0.70710678,  0}, {0,  0.70710678},
    };
    /* Executable statements */
    printf( " ZGESVD Example Program Results\n" );
    /* Query and allocate the optimal workspace */
    lwork = -1;
    zgesvd( "All", "All", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &wkopt, &lwork,
           rwork, &info );
    lwork = (int)wkopt.re;
    work = (dcomplex*)malloc( lwork*sizeof(dcomplex) );
    /* Compute SVD */
    zgesvd( "All", "All", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork,
           rwork, &info );
    /* Check for convergence */
    if( info > 0 ) {
        printf( "The algorithm computing SVD failed to converge.\n" );
        exit( 1 );
    }
    /* Print singular values */
    print_rmatrix( "Singular values", 1, m, s, 1 );
    /* Print left singular vectors */
    print_matrix( "Left singular vectors (stored columnwise)", m, m, u, ldu );
    /* Print right singular vectors */
    print_matrix( "Right singular vectors (stored rowwise)", m, n, vt, ldvt );
    /* Free workspace */
    free( (void*)work );
    exit( 0 );
} /* End of ZGESVD Example */

/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, int m, int n, dcomplex* a, int lda ) {
    int i, j;
    printf( "\n %s\n", desc );
    for( i = 0; i < m; i++ ) {
        for( j = 0; j < n; j++ )
        printf( " (%6.2f,%6.2f)", a[i+j*lda].re, a[i+j*lda].im );
        printf( "\n" );
    }
}

/* Auxiliary routine: printing a real matrix */
void print_rmatrix( char* desc, int m, int n, double* a, int lda ) {
    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+j*lda] );
        printf( "\n" );
    }
}

产量

<!-- language: lang-bash -->
ZGESVD Example Program Results

 Singular values
   1.00   1.00

 Left singular vectors (stored columnwise)
 ( -0.71,  0.00) ( -0.71,  0.00)
 ( -0.00,  0.71) (  0.00, -0.71)

 Right singular vectors (stored rowwise)
 ( -1.00, -0.00) ( -0.00, -0.00)
 ( -0.00, -0.00) ( -1.00, -0.00)

和scipy.linalg.SVD产量

<!-- language: lang-bash -->
Singular values: [1. 1.]
Left singular vectors: [[-0.70710678+0.j         -0.70710678+0.j        ]
                        [ 0.        +0.70710678j  0.        -0.70710678j]]
Right singular vectors: [[-1.+0.j -0.+0.j]
                         [-0.+0.j -1.+0.j]]

到目前为止，一切都很好。现在当我尝试输入一个3x3或NxN矩阵时，结果就像

// LAPACK C Case: 3x3

#include <stdlib.h>
#include <stdio.h>
#include <Accelerate/Accelerate.h>

/* Complex datatype */
struct _dcomplex { double re, im; };
typedef struct _dcomplex dcomplex;

/* ZGESVD prototype */
extern void zgesvd( char* jobu, char* jobvt, int* m, int* n, dcomplex* a,
                   int* lda, double* s, dcomplex* u, int* ldu, dcomplex* vt, int* ldvt,
                   dcomplex* work, int* lwork, double* rwork, int* info );
/* Auxiliary routines prototypes */
extern void print_matrix( char* desc, int m, int n, dcomplex* a, int lda );
extern void print_rmatrix( char* desc, int m, int n, double* a, int lda );

/* Parameters */
#define M 3
#define N 3
#define LDA M
#define LDU M
#define LDVT N

/* Main program */
int main() {
    /* Locals */
    int m = M, n = N, lda = LDA, ldu = LDU, ldvt = LDVT, info, lwork;
    dcomplex wkopt;
    dcomplex* work;
    /* Local arrays */
    /* rwork dimension should be at least max( 1, 5*min(m,n) ) */
    double s[M], rwork[5*M];
    dcomplex u[LDU*M], vt[LDVT*N];
    dcomplex a[LDA*N] = {
        {0, 0}, {0, 0}, {1,  0},
        {-0.36599657,  -0.27449743}, {-0.27449743,  0.36599657}, {0.76249285, 0},
        {-0.36599657, 0.27449743}, {-0.27449743, -0.36599657}, {0.76249285, 0},
    };
//    dcomplex a[LDA*N] = {
//        {0.70710678, 0}, {0, -0.70710678},
//        {0.70710678,  0}, {0,  0.70710678},
//    };
    /* Executable statements */
    printf( " ZGESVD Example Program Results\n" );
    /* Query and allocate the optimal workspace */
    lwork = -1;
    zgesvd( "All", "All", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &wkopt, &lwork,
           rwork, &info );
    lwork = (int)wkopt.re;
    work = (dcomplex*)malloc( lwork*sizeof(dcomplex) );
    /* Compute SVD */
    zgesvd( "All", "All", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork,
           rwork, &info );
    /* Check for convergence */
    if( info > 0 ) {
        printf( "The algorithm computing SVD failed to converge.\n" );
        exit( 1 );
    }
    /* Print singular values */
    print_rmatrix( "Singular values", 1, m, s, 1 );
    /* Print left singular vectors */
    print_matrix( "Left singular vectors (stored columnwise)", m, m, u, ldu );
    /* Print right singular vectors */
    print_matrix( "Right singular vectors (stored rowwise)", m, n, vt, ldvt );
    /* Free workspace */
    free( (void*)work );
    exit( 0 );
} /* End of ZGESVD Example */

/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, int m, int n, dcomplex* a, int lda ) {
    int i, j;
    printf( "\n %s\n", desc );
    for( i = 0; i < m; i++ ) {
        for( j = 0; j < n; j++ )
        printf( " (%6.2f,%6.2f)", a[i+j*lda].re, a[i+j*lda].im );
        printf( "\n" );
    }
}

/* Auxiliary routine: printing a real matrix */
void print_rmatrix( char* desc, int m, int n, double* a, int lda ) {
    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+j*lda] );
        printf( "\n" );
    }
}

产量

<!-- language: lang-bash -->
 ZGESVD Example Program Results

 Singular values
   1.55   0.65   0.42

 Left singular vectors (stored columnwise)
 (  0.26,  0.00) (  0.49, -0.34) ( -0.75,  0.00)
 (  0.20, -0.00) ( -0.66,  0.46) ( -0.57,  0.00)
 ( -0.94, -0.00) (  0.00,  0.00) ( -0.33,  0.00)

 Right singular vectors (stored rowwise)
 ( -0.61, -0.00) ( -0.56, -0.00) ( -0.56, -0.00)
 (  0.00,  0.00) (  0.41, -0.58) ( -0.41,  0.58)
 ( -0.79, -0.00) (  0.43, -0.00) (  0.43, -0.00)

<!-- language: lang-bash -->
# Python
Singular values: [1.55161905 0.64699664 0.41698163]
Left singular vectors: [[ 0.26480555-9.68622857e-18j  0.57973136-1.54633603e-01j
  -0.75490266+2.76133169e-17j]
 [ 0.19860416+1.15286199e-17j -0.77297515+2.06178138e-01j
  -0.56617699-3.28655711e-17j]
 [-0.94362832+0.00000000e+00j  0.        +0.00000000e+00j
  -0.33100694+0.00000000e+00j]]
Right singular vectors: [[-0.60815722+0.j         -0.5613131 +0.j         -0.5613131 +0.j        ]
 [ 0.        +0.j          0.18223745-0.68321996j -0.18223745+0.68321996j]
 [-0.7938166 +0.j          0.43003209+0.j          0.43003209+0.j        ]]

现在，我知道 Scipy 计算的结果是完美的，因为我使用 SVD 的目的是完美的并给出完美的结果，而我的目标是生成像 scipy 这样的结果。现在我知道 Scipy 也使用 LAPACK 的驱动程序，但为什么会有不同呢？我在哪里搞砸了。

4

1 回答 1

0

您应该测试的条件不是输出矩阵相同。

u @ np.diag(s) @ vh = A带有,的条件u, s, vh = np.linalg.svd(A)是您应该关注的地方。

还要注意奇异值的顺序，如果你有两组排序时相同的奇异值 s1 和 s2，你可以构造一个矩阵 P，P @ np.diag(s1) @ P.T = s2如果幸运的话，你可以使用u1 @ P.inv() = u2and P.T.inv() @ vh1 = vh2since (u1 @ P.inv()) @ P @ np.diag(s1) @ P.T @ (P.T.inv() @ vh1) = u1 @ np.diag(1) @ vh1。

于 2021-04-02T17:40:55.930 回答