Why does GNU Scientific Library not allow matrices with more columns than rows for singular value decomposition?

Question

MATLAB allows a matrix with more columns than rows when computing the singular value decomposition.

>> 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

But the GNU Scientific Library (GSL) does not. It gives this error:

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

Is this a shortcoming of the GSL or can it be worked around?

Answer 1

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;
}