问题描述
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;
}