0

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?

1 Answer 1

1

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 

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.