为什么 GNU Scientific Library 不允许列多于行的矩阵用于奇异值分解?

问题描述

MATLAB 在计算奇异值分解时允许矩阵的列数多于行数。

>> a_matrix = [1.0,0.0,2.0;
     0.0,3.0,0.0;
     0.0,2.0,0.0]

a_matrix =

     1     0     0     0     2
     0     0     3     0     0
     0     0     0     0     0
     0     2     0     0     0

>> [u,s,v] = svd(a_matrix)

u =

     0    -1     0     0
    -1     0     0     0
     0     0     0    -1
     0     0    -1     0


s =

    3.0000         0         0         0         0
         0    2.2361         0         0         0
         0         0    2.0000         0         0
         0         0         0         0         0


v =

         0   -0.4472         0         0   -0.8944
         0         0   -1.0000         0         0
   -1.0000         0         0         0         0
         0         0         0    1.0000         0
         0   -0.8944         0         0    0.4472

但 GNU 科学图书馆 (GSL) 没有。它给出了这个错误

gsl: svd.c:61: ERROR: svd of MxN matrix,M<N,is not implemented
Default GSL error handler invoked.

这是 GSL 的缺点还是可以解决

解决方法

int run_svd(const gsl_matrix * a) {
  // Need to transpose the input
  gsl_matrix *a_transpose = gsl_matrix_alloc(a->size2,a->size1);
  gsl_matrix_transpose_memcpy(a_transpose,a);
  printf("a_matrix'\n");
  pretty_print(a_transpose);

  int m = a->size1;
  gsl_matrix * V = gsl_matrix_alloc(m,m);
  gsl_vector * S = gsl_vector_alloc(m);
  gsl_vector * work = gsl_vector_alloc(m);

  gsl_linalg_SV_decomp(a_transpose,V,S,work);
  printf("U\n");
  pretty_print(V); printf("\n");
  printf("V\n");
  pretty_print(a_transpose); printf("\n");
  printf("S\n");
  gsl_vector_fprintf(stdout,"%g");

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

  return 0;
}