问题描述
void calculate_svd_example() {
int m = 4;
int n = 5;
double M[4][5] = {{1.0,0.0,2.0},{0.0,3.0,0.0},2.0,0.0}};
gsl_matrix *mat = gsl_matrix_alloc(m,n);
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
double x = M[i][j];
gsl_matrix_set(mat,i,j,x);
}
}
printf("M = \n");
pretty_print(mat);
run_svd(mat);
}
#include <stdio.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <time.h>
#include "../include/run_svd.h"
/*
gsl_matrix_printf prints a matrix as a column vector. This function
prints a matrix in block form.
*/
void pretty_print(const gsl_matrix * M)
{
// Get the dimension of the matrix.
int rows = M->size1;
int cols = M->size2;
// Now print out the data in a square format.
for(int i = 0; i < rows; i++){
for(int j = 0; j < cols; j++){
printf("%f ",gsl_matrix_get(M,j));
}
printf("\n");
}
printf("\n");
}
void pretty_print_vector(const gsl_vector * M)
{
int cols = M->size;
for(int j = 0; j < cols; j++){
printf("%f ",gsl_vector_get(M,j));
}
printf("\n");
}
int run_svd(const gsl_matrix * a) {
gsl_matrix *aa;
int m = a->size1;
int n = a->size2;
if (m >= n) {
aa = gsl_matrix_alloc(m,n);
gsl_matrix_memcpy(aa,a);
} else {
aa = gsl_matrix_alloc(n,m);
// Need to transpose the input
gsl_matrix_transpose_memcpy(aa,a);
}
m = aa->size2;
gsl_matrix * V = gsl_matrix_alloc(m,m);
gsl_vector * S = gsl_vector_alloc(m);
gsl_vector * work = gsl_vector_alloc(m);
/**
* On output the matrix A is replaced by U. The diagonal elements of
* the singular value matrix S are stored in the vector S. The
* singular values are non-negative and form a non-increasing sequence
* from S_1 to S_N. The matrix V contains the elements of V in
* untransposed form. To form the product U S V^T it is necessary to
* take the transpose of V. A workspace of length N is required in
* work.
*/
gsl_linalg_SV_decomp(aa,V,S,work);
printf("U:\n");
pretty_print(aa);
printf("S:\n");
pretty_print_vector(S);
printf("V:\n");
pretty_print(V);
gsl_matrix_free(V);
gsl_vector_free(S);
gsl_vector_free(work);
return 0;
}
它给了我以下结果:
U:
-0.000000 -0.447214 0.000000 -0.000000
0.000000 -0.000000 -1.000000 -0.000000
-1.000000 -0.000000 0.000000 0.000000
-0.000000 -0.000000 0.000000 1.000000
-0.000000 -0.894427 0.000000 0.000000
S:
3.000000 2.236068 2.000000 0.000000
V:
-0.000000 -1.000000 -0.000000 0.000000
-1.000000 -0.000000 -0.000000 0.000000
-0.000000 -0.000000 -0.000000 1.000000
-0.000000 -0.000000 -1.000000 0.000000
这里的 U
和 V
不是颠倒了吗?这是 GSL 代码、GSL documentation 的问题还是我做错了什么?
解决方法
您没有将完整代码作为 MWE 发布,所以我在下面做了一个新的实现。不幸的是,目前 GSL 中的 SVD 需要 M >= N,所以我计算矩阵转置的 SVD,然后从中输出正确的 U 和 V 因子。
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
void pretty_print(const gsl_matrix * M)
{
// Get the dimension of the matrix.
int rows = M->size1;
int cols = M->size2;
// Now print out the data in a square format.
for(int i = 0; i < rows; i++){
for(int j = 0; j < cols; j++){
printf("%f ",gsl_matrix_get(M,i,j));
}
printf("\n");
}
printf("\n");
}
void pretty_print_vector(const gsl_vector * M)
{
int cols = M->size;
for(int j = 0; j < cols; j++){
printf("%f ",gsl_vector_get(M,j));
}
printf("\n");
}
int
main()
{
const size_t M = 4;
const size_t N = 5;
double A_data[] = {
1.0,0.0,2.0,3.0,0.0 };
gsl_matrix_view A = gsl_matrix_view_array(A_data,4,5);
gsl_matrix * B = gsl_matrix_alloc(N,M);
gsl_matrix * V = gsl_matrix_alloc(M,M);
gsl_vector * S = gsl_vector_alloc(M);
gsl_vector * work = gsl_vector_alloc(M);
gsl_matrix_transpose_memcpy(B,&A.matrix);
gsl_linalg_SV_decomp(B,V,S,work);
printf("U:\n");
pretty_print(V);
printf("S:\n");
pretty_print_vector(S);
printf("V:\n");
pretty_print(B);
gsl_matrix_free(B);
gsl_matrix_free(V);
gsl_vector_free(S);
gsl_vector_free(work);
return 0;
}
输出如下:
$ ./svdtest
U:
-0.000000 -1.000000 -0.000000 0.000000
-1.000000 -0.000000 -0.000000 0.000000
-0.000000 -0.000000 -0.000000 1.000000
-0.000000 -0.000000 -1.000000 0.000000
S:
3.000000 2.236068 2.000000 0.000000
V:
-0.000000 -0.447214 0.000000 -0.000000
0.000000 -0.000000 -1.000000 -0.000000
-1.000000 -0.000000 0.000000 0.000000
-0.000000 -0.000000 0.000000 1.000000
-0.000000 -0.894427 0.000000 0.000000