为什么 GNU 科学图书馆不允许列数多于行数的矩阵进行奇异值分解?
Why does GNU Scientific Library not allow matrices with more columns than rows for singular value decomposition?
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 的缺点还是可以解决?
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;
}
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 的缺点还是可以解决?
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;
}