时间戳

首页 > 解决方案 > 为什么 GNU 科学库不允许列多于行的矩阵进行奇异值分解？

c - 为什么 GNU 科学库不允许列多于行的矩阵进行奇异值分解？

问题描述

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

>> a_matrix = [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]

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 的缺点还是可以解决？

标签： cmatlablinear-algebragsl

解决方案

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, S, "%g");

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

  return 0;
}

推荐阅读