svd.cpp

// ////////////////////////////////////////////////////////////////
// 特異値分解を用いて擬似逆行列を求める
// ////////////////////////////////////////////////////////////////
#include <iostream>

#include <tnt_array2d.h>
#include <tnt_array2d_utils.h>
#include <jama_svd.h>

template<class T>
TNT::Array2D<T>
transpose(const TNT::Array2D<T>& X)
{
  int m = X.dim1();
  int n = X.dim2();
  TNT::Array2D<T> Y(n, m);
  int i, j;
  for (i = 0; i < m; ++i) {
    for (j = 0; j < n; ++j) {
      Y[j][i] = X[i][j];
    }
  }
  return Y;
}

template<class T>
TNT::Array2D<T>
sinv(const TNT::Array2D<T>& X)
{
  int m = X.dim1();
  int n = X.dim2();
  TNT::Array2D<T> Y(m, n);
  int i, j;
  for (i = 0; i < m; ++i) {
    for (j = 0; j < n; ++j) {
      Y[i][j] = (X[i][j] == 0.0) ? 0.0 : 1.0/X[i][j];
    }
  }
  return Y;
}

int
main(int, char *[])
{
  TNT::Array2D<double> A(5, 4);
  A[0][0] = -2.0; A[0][1] =  1.0; A[0][2] =  2.0; A[0][3] =  4.0;
  A[1][0] = -3.0; A[1][1] =  1.0; A[1][2] =  2.0; A[1][3] =  2.0;
  A[2][0] =  1.0; A[2][1] = -2.0; A[2][2] = -3.0; A[2][3] = -2.0;
  A[3][0] = -4.0; A[3][1] =  1.0; A[3][2] =  3.0; A[3][3] =  2.0;
  A[4][0] = -5.0; A[4][1] =  0.0; A[4][2] = -5.0; A[4][3] =  4.0;
  std::cout << "A: " << std::endl;
  std::cout << A << std::endl;

  JAMA::SVD<double> svd(A);
  TNT::Array2D<double> U;
  svd.getU(U);
  TNT::Array2D<double> S;
  svd.getS(S);
  TNT::Array2D<double> V;
  svd.getV(V);
  TNT::Array2D<double> pinv = TNT::matmult(V, TNT::matmult(sinv(S), transpose(U)));
  std::cout << "pseudo-inverse matrix: " << std::endl;
  std::cout << pinv << std::endl;

  std::cout << "pinv(A) * A: " << std::endl;
  std::cout << TNT::matmult(pinv, A) << std::endl;

  return 0;
}