GSL gsl_linalg_SV_decomp 返回 U 和 V 反转

Question

我有以下 C 代码：

void calculate_svd_example() {
  int m = 4;
  int n = 5;
  double M[4][5] = {{1.0, 0.0, 0.0, 0.0, 2.0},
                     {0.0, 0.0, 3.0, 0.0, 0.0},
                     {0.0, 0.0, 0.0, 0.0, 0.0},
                     {0.0, 2.0, 0.0, 0.0, 0.0}};

  gsl_matrix *mat = gsl_matrix_alloc(m, n);
  for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
          double x = M[i][j];
          gsl_matrix_set(mat, i, j, x);
      }
  }
  printf("M = \n");
  pretty_print(mat);
  run_svd(mat);
}

---

#include <stdio.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <time.h>

#include "../include/run_svd.h"

/*
  gsl_matrix_printf prints a matrix as a column vector.  This function
  prints a matrix in block form.
*/
void pretty_print(const gsl_matrix * M)
{
  // Get the dimension of the matrix.
  int rows = M->size1;
  int cols = M->size2;
  // Now print out the data in a square format.
  for(int i = 0; i < rows; i++){
    for(int j = 0; j < cols; j++){
      printf("%f ", gsl_matrix_get(M, i, j));
    }
    printf("\n");
  }
  printf("\n");
}

void pretty_print_vector(const gsl_vector * M)
{
  int cols = M->size;
    for(int j = 0; j < cols; j++){
      printf("%f ", gsl_vector_get(M, j));
    }
  printf("\n");
}

int run_svd(const gsl_matrix * a) {
  gsl_matrix *aa;
  int m = a->size1;
  int n = a->size2;
  if (m >= n) {
    aa = gsl_matrix_alloc(m, n);
    gsl_matrix_memcpy(aa, a);
  } else {
    aa = gsl_matrix_alloc(n, m);
    // Need to transpose the input
    gsl_matrix_transpose_memcpy(aa, a);
  }

  m = aa->size2;
  gsl_matrix * V = gsl_matrix_alloc(m, m);
  gsl_vector * S = gsl_vector_alloc(m);
  gsl_vector * work = gsl_vector_alloc(m);

  /**
   * On output the matrix A is replaced by U. The diagonal elements of 
   * the singular value matrix S are stored in the vector S. The 
   * singular values are non-negative and form a non-increasing sequence 
   * from S_1 to S_N. The matrix V contains the elements of V in 
   * untransposed form. To form the product U S V^T it is necessary to 
   * take the transpose of V. A workspace of length N is required in 
   * work.
   */
  gsl_linalg_SV_decomp(aa, V, S, work);
  printf("U:\n");
  pretty_print(aa);
  printf("S:\n");
  pretty_print_vector(S);
  printf("V:\n");
  pretty_print(V);

  gsl_matrix_free(V);
  gsl_vector_free(S);
  gsl_vector_free(work);

  return 0;
}

---

它给了我以下结果：

U:
-0.000000 -0.447214  0.000000 -0.000000 
 0.000000 -0.000000 -1.000000 -0.000000 
-1.000000 -0.000000  0.000000  0.000000 
-0.000000 -0.000000  0.000000  1.000000 
-0.000000 -0.894427  0.000000  0.000000 

S:
3.000000 2.236068 2.000000 0.000000 

V:
-0.000000 -1.000000 -0.000000 0.000000 
-1.000000 -0.000000 -0.000000 0.000000 
-0.000000 -0.000000 -0.000000 1.000000 
-0.000000 -0.000000 -1.000000 0.000000 

U 和 V 不是颠倒过来了吗？这是 GSL 代码、GSL documentation 的问题，还是我做错了什么？

Answer 1

您没有将完整代码发布为 MWE，因此我在下面做了一个新的实现。不幸的是，目前 GSL 中的 SVD 需要 M >= N，所以我计算矩阵转置的 SVD，然后从中输出正确的 U 和 V 因子。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>

void pretty_print(const gsl_matrix * M)
{
  // Get the dimension of the matrix.
  int rows = M->size1;
  int cols = M->size2;
  // Now print out the data in a square format.
  for(int i = 0; i < rows; i++){
    for(int j = 0; j < cols; j++){
      printf("%f ", gsl_matrix_get(M, i, j));
    }
    printf("\n");
  }
  printf("\n");
}

void pretty_print_vector(const gsl_vector * M)
{
  int cols = M->size;
    for(int j = 0; j < cols; j++){
      printf("%f ", gsl_vector_get(M, j));
    }
  printf("\n");
}

int
main()
{
  const size_t M = 4;
  const size_t N = 5;
  double A_data[] = {
    1.0, 0.0, 0.0, 0.0, 2.0,
    0.0, 0.0, 3.0, 0.0, 0.0,
    0.0, 0.0, 0.0, 0.0, 0.0,
    0.0, 2.0, 0.0, 0.0, 0.0 };
  gsl_matrix_view A = gsl_matrix_view_array(A_data, 4, 5);
  gsl_matrix * B = gsl_matrix_alloc(N, M);
  gsl_matrix * V = gsl_matrix_alloc(M, M);
  gsl_vector * S = gsl_vector_alloc(M);
  gsl_vector * work = gsl_vector_alloc(M);

  gsl_matrix_transpose_memcpy(B, &A.matrix);

  gsl_linalg_SV_decomp(B, V, S, work);

  printf("U:\n");
  pretty_print(V);

  printf("S:\n");
  pretty_print_vector(S);

  printf("V:\n");
  pretty_print(B);

  gsl_matrix_free(B);
  gsl_matrix_free(V);
  gsl_vector_free(S);
  gsl_vector_free(work);

  return 0;
}

这是输出：

$ ./svdtest
U:
-0.000000 -1.000000 -0.000000 0.000000 
-1.000000 -0.000000 -0.000000 0.000000 
-0.000000 -0.000000 -0.000000 1.000000 
-0.000000 -0.000000 -1.000000 0.000000 

S:
3.000000 2.236068 2.000000 0.000000 
V:
-0.000000 -0.447214 0.000000 -0.000000 
0.000000 -0.000000 -1.000000 -0.000000 
-1.000000 -0.000000 0.000000 0.000000 
-0.000000 -0.000000 0.000000 1.000000 
-0.000000 -0.894427 0.000000 0.000000