I have the following C code:
void calculate_svd_example() {
int m = 4;
int n = 5;
double M[4][5] = {{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}};
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, 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 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;
}
It is giving me the following results:
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
Aren't U and V reversed here? Is this a problem with the GSL code, the GSL documentation, or something I'm doing wrong?
You didn't post your full code as a MWE, so I made a new implementation below. Unfortunately, currently the SVD in GSL requires M >= N, so I compute the SVD of the transpose of the matrix, and then output the correct U and V factors from that.
#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, 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 };
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;
}
Here is the output:
$ ./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