概要

今回は、Fortranの有名なライブラリLAPACKを用いて、特異値分解する方法を紹介したいと思います。
数学的には以下のようなこのことを行っています。

A = U \Sigma V^{t}

Aがm×n行列であれば、Uはm×min(m,n)行列,シグマはmin(m,n)×min(m,n)行列（対角成分に特異値、それ以外は0の行列,対角成分は特異値の値の大きいものから並んでいる。）、Vtはmin(m,n)×n。

特異値分解の具体的な説明は他のページにお任せするとして、使い方だけを説明します。LAPACKはFORTRAN77で書かれており、とても使いづらいものです。そのため、Fortran95で書かれ使いやすくしたLAPACK95を用いることにします。
これらのダウンロードとインストール方法は以下のページを参照しました。
(https://qiita.com/AnchorBlues/items/69c1744de818b5e045ab)

方法1：テストデータ作成

まず、以下のようなデータファイルを用意し、作業ディレクトリに保存する。

test_data_lapack.dat

2.27  -1.54   1.15  -1.94
0.28  -1.67   0.94  -0.78
0.48  -3.09   0.99  -0.21
1.07   1.22   0.79   0.63
-2.35   2.93  -1.45   2.30
0.62  -7.39   1.03  -2.57

方法2：プログラムファイル作成

さらにプログラムファイルを用意する。

test_lapack95.f90

program lapack_test
  use f95_lapack !ここでLAPACK95を用いる宣言

  implicit none
  integer :: i, m, n, info, small_mn
  real(8) :: A(1:6,1:4)
  real(8),allocatable :: S(:), U(:,:), Vt(:,:)
  character :: filename*120

  m = ubound(A,1)
  n = ubound(A,2)
  small_mn = min(m,n)!small_mnは整数m,nの小さい方の整数値を取る。
  allocate(S(1:small_mn), U(1:m,1:small_mn),Vt(1:small_mn,1:n))!small_mnを決めてからallocateする。

  filename = "./test_data_lapack.dat"
  open(20,file=filename, status = 'old')
  do i = 1, 6
     read(20,*) A(i,1:4)!用意したデータの読み込み
  enddo
  close(20)

  call LA_GESVD(A(1:m,1:n), S(1:small_mn), U=U(1:m,1:small_mn), & Vt = Vt(1:small_mn,1:n), JOB ="N", INFO = info )!ここで特異値分解を行っている。

  if(info < 0) then !もし特異値分解が失敗していれば、Failと表示される。
     write(*,*) info
     stop "Fail"
  endif

!以下結果の書き出し
  write(*,*) "sigma="
  write(*,*) sigma(1:small_mn)
  write(*,*) "U="
  write(*,*) U(1:m,1:small_mn)
  write(*,*) "Vt="
  write(*,*) Vt(1:small_mn,n)
  stop
end program lapack_test

実行方法

gfortran test_lapack95.f90 -I/usr/local/include -llapack95 -llapack -lblas
./a.out

解説1:変数の説明

Aという配列(行列)が今回特異値分解したいと考えている行列。（この行列をバラバラにする）subroutineではintent(inout)この計算ではinと思えば良い。
Sは特異値を集めた配列である。大きい値の特異値から順に並べられている。通常ではシグマという行列の対角成分だけをとってきたものである。subroutineではintent(out)
Uは左特異ベクトル。subroutineではintent(out)
Vtは右特異ベクトル。subroutineではintent(out)
JOB="N"はAに値を返すかということを決めるものであるが、特に気にしなくてよい。
INFO=infoはsubroutine LA_GESVDの計算が成功した場合には0それ以外には失敗の原因により、負の整数が与えられる。subroutineではintent(out)

解説2：重要な点

small_mnをしっかりと計算することである。これを計算しallocate。もしくは自分の頭でどちらが小さいかを予め計算してからU,Vt等の定義をすることである。これらの配列のサイズをうまく指定しておかないとおかしなことになる。

解説3：subroutine LA_GESVD

LA_GESVDというのがLAPACK95内に含まれているサブルーチンである。目的は特異値分解であるので、右特異ベクトルVt、左特異ベクトルUを求める際にはキーワードを用いてオプションで指定することで計算できるようになっている。ただし右特異ベクトルは転置された状態になって格納されていることに注意(Vtのtは転置のtの意味)キーワードの部分は省略でき、もし特異値のみを得たい場合はU=,Vt=といった＝の入った部分は入力する必要はない。
またこのサブルーチンを使用する際には

このサブルーチンではプログラムの冒頭において
use f95_lapack
という宣言が必要である。

参考文献

元のデータ
 LAPACK95

最後に

ご質問等あればお気軽にどうぞ。
返信が遅い可能性もあります。
間違った記述がある場合もあります。
他のサブルーチンについては詳しくないです。
以上、ご了承ください。

LAPACK95で特異値分解する

概要