gslで特異値分解してみた

特異値分解を使うことになったからとりあえずライブラリで特異値分解を実行してみた．

gslには特異値分解を行う

 int gsl_linalg_SV_decomp (gsl_matrix * A, gsl_matrix * V, gsl_vector * S, gsl_vector * work)

という関数が用意されている．

この関数は

と分解する関数で，引数で与えられるAは計算で得られるUで置き換えられるようにできているみたい

とりあえず特異値分解するソースコードはこんな感じ

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_linalg.h>

void my_gsl_matrix_print(const gsl_matrix *matrix, int n, int m);//行列を行列の形で出力

int main (void)
{
  //int n ;
  double a_data[4][3] = { {1.00, 1.00, 0.00},
              {1.00, 1.00, 1.00},
              {0.00, 1.00, 1.00},
              {0.00, 1.00, 1.00}} ;
  int n=4,m=3;
  int i,j;
  double vec[n*m];
  for(i=0;i<n;i++){
    for(j=0;j<m;j++){
      vec[i*m+j]=a_data[i][j];
    }
  }
  gsl_matrix_view A = gsl_matrix_view_array (vec, n, m);
  printf("\n") ;
  gsl_matrix *V = gsl_matrix_alloc (m,m);
  gsl_vector *S = gsl_vector_alloc (m);
  gsl_vector *work = gsl_vector_alloc (m);
  gsl_matrix *S_matrix = gsl_matrix_alloc(m,m);
  gsl_matrix *temp = gsl_matrix_alloc(n,m);
  gsl_matrix *ans = gsl_matrix_alloc(n,m);
  printf("A is\n");
  my_gsl_matrix_print(&A.matrix,n,m);
  
  //  int gsl_linalg_SV_decomp (gsl_matrix * A, gsl_matrix * V, gsl_vector * S, gsl_vector * work)
  gsl_linalg_SV_decomp (&A.matrix, V, S, work) ;//特異値分解
  
  for(i=0;i<m;i++){
    gsl_matrix_set(S_matrix,i,i,gsl_vector_get(S,i));//ベクトルSを行列S_matrixに
  }
  
  printf("U is\n");
  my_gsl_matrix_print(&A.matrix,n,m);
  printf("S is\n");
  my_gsl_matrix_print(S_matrix,m,m);
  printf("V is\n") ;
  my_gsl_matrix_print(V,m,m);
  printf("\n") ;
  
  gsl_matrix_transpose(V);//Vを転置
  
  /*分解したものをかけて確認*/
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, &A.matrix, S_matrix, 0.0, temp);
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, temp,V, 0.0, ans);
  my_gsl_matrix_print(ans,n,m);
  return 0;
}

void my_gsl_matrix_print(const gsl_matrix *matrix, int n, int m){
  int i,j;
  for(i=0;i<n;i++){
    for(j=0;j<m;j++){
      printf("%15g",gsl_matrix_get(matrix,i,j));
    }
    printf("\n");
  }
}

実行結果はこんな感じ

 A is
          1              1              0
          1              1              1
          0              1              1
          0              1              1 
 U is
  -0.403846       0.725356      -0.557464
  -0.615815       0.235084       0.752002
  -0.478373       -0.45749      -0.248723
  -0.478373       -0.45749      -0.248723
 S is
    2.72375              0              0
          0        1.17761              0
          0              0       0.440905
 V is
  -0.374359       0.815583       0.441225
  -0.725619      0.0386051      -0.687013
   -0.57735       -0.57735        0.57735

          1              1   -8.88178e-16
          1              1              1
-5.41234e-16              1              1
-4.44089e-16              1              1

まあ，とりあえず特異値分解できてるっぽい