線形代数-固有値(LAEigen)

file:/Users/bob/Github/TeamNishitani/jupyter_num_calc/laeigenvectors
https://github.com/daddygongon/jupyter_num_calc/tree/master/notebooks_python
cc by Shigeto R. Nishitani 2017

固有値

$A$を対称正方行列,$x$をベクトルとしたときに, $$ Ax = \lambda x $$ の解,$\lambda$を固有値,$x$を固有ベクトルという.$x$がゼロベクトルではない意味のある解は特性方程式det$(A-\lambda E)=0$が成り立つときにのみ得られる.

この特性方程式をpythonで特には,sympyってのを使わないといけない.一応示しておくと,

$$ \begin{align} A\, &= \, \left( \begin {array}{cc} 3&2/3\\2/3&2\end {array} \right) \\ \det( A\, - t E) &= \, \det \left( \begin {array}{cc} 3-t&2/3\\2/3& 2-t \end {array} \right) \\ & = t^2 -5t + 6 - 4/9 \\ & = t^2 -5t + 50/9 \\ t (\lambda) &= 10/3,\,5/3 \end{align} $$
In [1]:
from sympy import *
t = symbols('t')
A = Matrix([[3, 2/3],[2/3,2]]
)
xx = diag(t,t)
pprint(xx)
eq1 = simplify((A-xx).det())
pprint(eq1)
solve(eq1,t)
⎡t  0⎤
⎢    ⎥
⎣0  t⎦
 2                         
t  - 5⋅t + 5.55555555555556
Out[1]:
[1.66666666666667, 3.33333333333333]

固有値を求めるコマンドEigenvectorsを適用すると,固有値と固有ベクトルが求まる.ここで,固有ベクトルは行列の列(Column)ベクトルに入っている.

In [2]:
import pprint
import scipy.linalg   # SciPy Linear Algebra Library

A = scipy.array([[3, 2/3],[2/3,2]])
l, V = scipy.linalg.eig(A)
pprint.pprint(l)
pprint.pprint(V)
pprint.pprint(V[:,0]) # columnの取り出し方
array([ 3.33333333+0.j,  1.66666667+0.j])
array([[ 0.89442719, -0.4472136 ],
       [ 0.4472136 ,  0.89442719]])
array([ 0.89442719,  0.4472136 ])

固有ベクトルは規格化されている.

Googleのページランク

多くの良質なページからリンクされているページはやはり良質なページである

Googleのpage rankは上のような非常に単純な仮定から成り立っている.ページランクを実際に求めよう.つぎのようなリンクが張られたページを考える.

link structure

計算手順は以下の通り

  1. リンクを再現する隣接行列を作る.ページに番号をつけて,その間が結ばれているi-j要素を1,そうでない要素を0とする.
  2. 隣接行列を転置する
  3. 列ベクトルの総和が1となるように規格化する.
  4. こうして得られた推移確率行列の最大固有値に属する固有ベクトルを求め,適当に規格化する.

課題

  1. 上記手順を参考にして,pythonでページランクを求めよ.
  2. このような問題ではすべての固有値・固有ベクトルを求める必要はなく,最大の固有値を示す固有ベクトルを求めるだけでよい.初期ベクトルを適当に決めて,何度も推移確率行列を掛ける反復法でページランクを求めよ.
    隣接行列
$$ {\it A1}\, := \, \left[ \begin {array}{c|c|c|c|c|c|c|c} &1&2&3&4&5&6&7\\ 1&0&1&1&1&1&0&1\\ 2&1&0&0&0&0&0&0\\ 3& & & & & & & \\ 4& & & & & & & \\ 5& & & & & & & \\ 6& & & & & & & \\ 7& & & & & & & \end {array} \right] $$
転置行列
$$ {Transpose}({\it A1})\, := \, \left[ \begin {array}{c|c|c|c|c|c|c} \, \, &\, \, &\, \, &\, \, &\, \, &\, \, &\, \, \\ & & & & & & \\ & & & & & & \\ & & & & & & \\ & & & & & & \\ & & & & & & \\ & & & & & & \end {array} \right] $$
規格化
$$ \left[ \begin {array}{c|c|c|c|c|c|c} \, \, &\, \, &\, \, &\, \, &\, \, &\, \, &\, \, \\ & & & & & & \\ & & & & & & \\ & & & & & & \\ & & & & & & \\ & & & & & & \\ & & & & & & \end {array} \right] $$
遷移
$$ \left( \begin {array}{ccccccc} 0 &1 &1/2 &0 &1/4 &1/2 &0 \\ 1/5 &0 &1/2 &1/3 &0 &0 &0 \\ 1/5 &0 &0 &1/3 &1/4 &0 &0 \\ 1/5 &0 &0 &0 &1/4 &0 &0 \\ 1/5 &0 &0 &1/3 &0 &1/2 &1 \\ 0 &0 &0 &0 &1/4 &0 &0 \\ 1/5 &0 &0 &0 &0 &0 &0 \end {array} \right) \left( \begin {array}{c} 1/7\\ 1/7\\ 1/7\\ 1/7\\ 1/7\\ 1/7\\ 1/7 \end {array} \right) \, = \, \left( \begin {array}{ccccccc} \, \, &\, \, &\, \, &\, \, &\, \, &\, \, &\, \, \\ & & & & & & \\ & & & & & & \\ & & & & & & \\ & & & & & & \\ & & & & & & \\ & & & & & & \end {array} \right) \, = \, \left( \begin {array}{c} 0.32\\ 0.15\\ 0.11\\ 0.06\\ 0.29\\ 0.04\\ 0.03 \end {array} \right) $$
In [27]:
from pprint import pprint
from numpy import array, zeros, diagflat, dot, transpose
from scipy.linalg import eig

A = array([[0,1,1,1,1,0,1],
           [1,0,0,0,0,0,0],
           [1,1,0,0,0,0,0],
           [0,1,1,0,1,0,0],
           [1,0,1,1,0,1,0],
           [1,0,0,0,1,0,0],
           [0,0,0,0,1,0,0]])

diag = []
for i in range(0,7):
    tmp = 0.0
    for j in range(0,7):
        tmp += A[i,j]
    diag.append(1.0/tmp)

D = diagflat(diag)
tA = dot(transpose(A),D)

x = array([1,0,0,0,0,0,0])
pprint(dot(tA,dot(tA,x)))

l, V = eig(tA)
v0 = V[:,0]
pprint(v0)
array([ 0.350,  0.167,  0.117,  0.050,  0.267,  0.050,  0.000])
array([ 0.69945653+0.j,  0.38286042+0.j,  0.32395882+0.j,  0.24296911+0.j,
        0.41231122+0.j,  0.10307780+0.j,  0.13989131+0.j])
In [32]:
import numpy as np
from pprint import pprint

np.set_printoptions(formatter={'float': '{: 0.3f}'.format}) 
#pprint(A)
pprint(D)
#pprint(transpose(A))
pprint(dot(transpose(A),D))
x = array([1,0,0,0,0,0,0])
#pprint(x)
#pprint(dot(tA,dot(tA,dot(tA,x))))
l, V = eig(tA)
v0 = V[:,0]
#pprint(v0)
array([[ 0.200,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000],
       [ 0.000,  1.000,  0.000,  0.000,  0.000,  0.000,  0.000],
       [ 0.000,  0.000,  0.500,  0.000,  0.000,  0.000,  0.000],
       [ 0.000,  0.000,  0.000,  0.333,  0.000,  0.000,  0.000],
       [ 0.000,  0.000,  0.000,  0.000,  0.250,  0.000,  0.000],
       [ 0.000,  0.000,  0.000,  0.000,  0.000,  0.500,  0.000],
       [ 0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  1.000]])
array([[ 0.000,  1.000,  0.500,  0.000,  0.250,  0.500,  0.000],
       [ 0.200,  0.000,  0.500,  0.333,  0.000,  0.000,  0.000],
       [ 0.200,  0.000,  0.000,  0.333,  0.250,  0.000,  0.000],
       [ 0.200,  0.000,  0.000,  0.000,  0.250,  0.000,  0.000],
       [ 0.200,  0.000,  0.000,  0.333,  0.000,  0.500,  1.000],
       [ 0.000,  0.000,  0.000,  0.000,  0.250,  0.000,  0.000],
       [ 0.200,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000]])

累乗(べき乗)法により最大固有値が求まる原理

累乗(べき乗)法は,最大固有値とその固有ベクトルを効率的に見つける算法である.すこし,固有値について復習しておく.正方行列$A$に対して, $$ A x = \lambda x $$ の解$\lambda$を固有値,$x$を固有ベクトルという.$\lambda$は, $$ \det( A - \lambda E) =0 $$ として求まる永年方程式の解である.

では,なぜ適当な初期ベクトル$x_0$から始めて,反復 $$ x_{k+1} = A x_k $$ を繰り返すと,$A$の絶対値最大の固有値に属する固有ベクトルに近づいていくのかを見ておこう.

すべての固有値がお互いに異なる場合を考える.今,行列の固有値を絶対値の大きなもの順に並べて,$|\lambda_1|>|\lambda_2|>\cdots>|\lambda_n|$とし,対応する長さを1に規格化した固有ベクトルを$x_1, x_2, \ldots, x_n$とする.初期ベクトルは固有ベクトルの線形結合で表わせて, $$ X_0 = c_1x_1+c_2x_2+\cdots+c_nx_n $$ となるとする.これに行列$A$を$N$回掛けると, $$ A^N X_0 = c_1 \lambda_1^N x_1+ c_2 \lambda_2^N x_2+\cdots+ c_n \lambda_n^N x_n $$ となる.これを変形すると, $$ A^NX_0 = X_{N} = c_1 \lambda_1^N \left\{ x_1+ \frac{c_2}{c_1}\left(\frac{\lambda_2}{\lambda_1}\right)^N x_2+\cdots+ \frac{c_n}{c_1}\left(\frac{\lambda_n}{\lambda_1}\right)^N x_n \right\} $$ となる.$|\lambda_1|>|\lambda_i|(i\ge2)$だから括弧の中は$x_1$だけが生き残る.

こうして最大固有値に属する固有ベクトルが,反復計算を繰り返すだけで求められる.

Jacobi回転による固有値の求め方

固有値を求める手法として,永年方程式を解くというやり方は回りくどすぎる.少し古めかしいが非対角要素を0にする回転行列を反復的に作用させるJacobi(ヤコビ)法を紹介する.現在認められている最適の方策は,ハウスホルダー(Householder)変換で行列を単純な三重対角化行列に変形してから,反復法で解を追い込んでいくやり方である.Jacobi法は,Householder法ほど万能ではないが,10次程度までの行列には今でも役に立つ.

Mapleでみる回転行列

行列の軸回転の復習をする.対称行列$B$に回転行列$U$を作用すると $$ B.U = \left( \begin{array}{cc} {a_{1\,1}} & {a_{1\,2}}\\ {a_{2\,1}(={a_{1\,2}})} & {a_{2\,2}} \end{array} \right) \left( \begin{array}{cc} \cos(\theta) & -\sin(\theta)\\ \sin(\theta) & \cos(\theta) \end{array} \right) $$

となる.回転行列を4x4の行列に $$ U^t B U $$ と作用させたときの各要素の様子を以下に示した.

maple
> restart:
> n:=4:
> with(LinearAlgebra):
> B:=Matrix(n,n,shape=symmetric,symbol=a);
$$ B := \left[{ \begin{array}{cccc} {a_{1, \,1}} & {a_{1, \,2}} & {a_{1, \,3}} & {a_{1, \,4}} \\ {a_{1, \,2}} & {a_{2, \,2}} & {a_{2, \,3}} & {a_{2, \,4}} \\ {a_{1, \,3}} & {a_{2, \,3}} & {a_{3, \,3}} & {a_{3, \,4}} \\ {a_{1, \,4}} & {a_{2, \,4}} & {a_{3, \,4}} & {a_{4, \,4}} \end{array}} \right] $$
maple
> U:=Matrix(n,n,[[c,-s,0,0],[s,c,0,0],[0,0,1,0],[0,0,0,1]]);

#U:=Matrix(n,n,[[c,-s],[s,c]]);
$$ U := \left[ {\begin{array}{ccrr} c & - s & 0 & 0 \\ s & c & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array}} \right] $$
maple
>TT:=Transpose(U).B.U;
$$ \mathit{TT} := \\ \notag {\begin{array}{c} \left[ \right. (c\,{a_{1, \,1}} + s\,{a_{1, \,2}})\,c + (c\,{a_{1, \,2}} + s \,{a_{2, \,2}})\,s\,, \, - (c\,{a_{1, \,1}} + s\,{a_{1, \,2}})\,s + (c\,{a_{1, \,2}} + s\,{a_{2, \,2}})\,c\,, \\ c\,{a_{1, \,3}} + s\,{a_{2, \,3}}\,, \,c\,{a_{1, \,4}} + s\,{a_{2 , \,4}} \left. \right] \\ \left[ \right. ( - s\,{a_{1, \,1}} + c\,{a_{1, \,2}})\,c + ( - s\,{a_{1, \,2 }} + c\,{a_{2, \,2}})\,s\,, \, - ( - s\,{a_{1, \,1}} + c\,{a_{1, \,2}})\,s + ( - s\,{a_{1, \,2}} + c\,{a_{2, \,2}})\,c\,, \\ - s\,{a_{1, \,3}} + c\,{a_{2, \,3}}\,, \, - s\,{a_{1, \,4}} + c \,{a_{2, \,4}} \left. \right] \\ \left[ c\,{a_{1, \,3}} + s\,{a_{2, \,3}}\,, \, - s\,{a_{1, \,3}} + c\,{a_{2, \,3}}\,, \,{a_{3, \,3}}\,, \,{a_{3, \,4}} \right] \\ \left[ c\,{a_{1, \,4}} + s\,{a_{2, \,4}}\,, \, - s\,{a_{1, \,4}} + c\,{a_{2, \,4}}\,, \,{a_{3, \,4}}\,, \,{a_{4, \,4}} \right] \end{array}} $$
maple
>expand(TT[1,1]);
expand(TT[2,2]);
expand(TT[1,2]);
expand(TT[2,1]);
$$ c^{2}\,{a_{1, \,1}} + 2\,c\,s\,{a_{1, \,2}} + s^{2}\,{a_{2, \,2}} $$$$ s^{2}\,{a_{1, \,1}} - 2\,c\,s\,{a_{1, \,2}} + c^{2}\,{a_{2, \,2}} $$$$ - s\,c\,{a_{1, \,1}} - s^{2}\,{a_{1, \,2}} + c^{2}\,{a_{1, \,2}} + c\,s\,{a_{2, \,2}} $$$$ - s\,c\,{a_{1, \,1}} - s^{2}\,{a_{1, \,2}} + c^{2}\,{a_{1, \,2}} + c\,s\,{a_{2, \,2}} $$

この非対角要素を0にする$\theta$は以下のように求まる.

        

このとき注目している$i,j=1,2$以外の要素も変化する.

maple
>expand(TT[3,1]);
expand(TT[3,2]);
$$ c\,{a_{1, \,3}} + s\,{a_{2, \,3}} $$$$ - s\,{a_{1, \,3}} + c\,{a_{2, \,3}} $$

これによって一旦0になった要素も値を持つが,なんども繰り返すことによって,徐々に0へ近づいていく.

C 言語での実装

lapack_codes/c_versions/Jacobi2.c にコードを入れている.

$ gcc Jacobi2.c
$ a.out < input.txt

で固有値を求めていく過程が表示される.

数値計算ライブラリーについて

  • この辺りに記述しているC, Python, Maple, ruby言語のコードは,numerical_calc/lapack_codesにおいている

一般の数値計算ライブラリーについては,時間の関係で講義ではその能力を紹介するにとどめる.昔の演習で詳しく取り上げていたので,研究や今後のために必要と思うときは,テキストを取りにおいで.

行列の計算は,数値計算の中でも特に利用する機会が多く,また,律速ルーチンとなる可能性が高い.そこで,古くから行列計算の高速ルーチンが開発されてきた.なかでもBLASとLAPACKはフリーながら非常に高速である.

前回に示した,逆行列を求める単純なLU分解法をC言語でコーディングしたものと,LAPACKのルーチンを比べた場合,1000次元の行列で計測すると

maple
>  1000 [dim]     2.5200 [sec] #BOB
>  1000 [dim]     0.4700 [sec] #LAPACK

となった.2006年に初めてこの計算に用いたPCはMacBook(2GHz, Intel Core Duo)であるが, この計算での0.47秒は1.4GFLOPに相当する. 07年のMacBook(2GHz, Intel Core 2 Duo)ではさらに速くなって

maple
bob% gcc -O3 bob.c -o bob
bob% ./bob
1000
 1000 [dim]     1.7543 [sec] #BOB
bob% gcc -O3 lapack.c -llapack -lblas -o lapack
bob% ./lapack
1000
 1000 [dim]     0.1893 [sec] #LAPACK

で,3.5GFLOPSが出ている.今(2016年)は,MacBookAir(2.2GHz, Intel Core i7)で... top500.orgが毎年2回High Performance Computerのランクを発表している. 今は,Top1は100PFlopsであるが, 初回の1994年6月の500位は0.4GFlopsで,今のlaptopがはるかに凌いでいる. まさにlaptopスパコンの時代なんですよ.

ライブラリーは世界中の計算機屋さんがよってたかって検証しているので,バグがほとんど無く,また,高速である. 初学者はライブラリーを使うべきである.ただし,下のサンプルプログラムの行列生成の違いのように,ブラックボックス化すると思わぬ間違い(ここではFortranとCでの行列の並び順の違いが原因)をしでかすことがあるので,プログラムに組み込む前に必ず小さい次元(サンプルコード)で検証しておくこと\footnote{少し前(2002年ごろ)GotoBLASが開発されて,性能が10%ほども上がった}.

添付のコードはちょっと長いが時間があればフォローせよ.コンパイルは,OSXでは

maple
> gcc -O3 -UPRINT lapack.c -llapack -lblas

とすればできる.linuxではLAPACK, BLASがインストールされていれば,

maple
> #include <vecLib/vecLib.h>

をコメントアウトして,

maple
> gcc -O3 -DPRINT lapack.c -L/usr/local/lib64 -llapack -lblas -lg2c

などとすればコンパイルできるはず.

python謹製 lapack利用逆行列

In [11]:
import numpy as np
import scipy.linalg as la
import time

sizes = [1000, 2000, 4000]
for n in sizes:
    A = np.random.random((n, n))
    b = np.random.random((n))
    start = time.time()
    la.lapack.dgesv(A, b)
    print('%s [dim], %10.5f [sec] # python' % (n,time.time()-start))
1000 [dim],    0.04338 [sec] # python
2000 [dim],    0.13708 [sec] # python
4000 [dim],    0.95414 [sec] # python
size elapsed_time[sec]
100 0.004878
1000 0.085786
10000 46.031718

計算速度競争

Mapleと,C, pythonで競争させてみました.Cはそのままです.Mapleは

with(LinearAlgebra):
data:=[1000,2000,4000];
for n in data do
  A:=RandomMatrix(n,n,generator = 0..1.0):
  b:=RandomVector(n,generator = 0..1.0):
  st:=time();
  LUDecomposition(<A|b>):
  print(time()-st);
end:

です.pythonのcodeは,lapackにあるdgesvを指定して 呼び出すようにしています.scipy.laのsolveでは何を使っているかはよくわからないので.

結果は,次の通りです.

size Maple[sec] C[sec] python[sec]
1000 0.339 0.0589 0.04101
2000 1.384 0.3826 0.12682
4000 45.482 3.0111 0.81373

pythonの圧勝ですね.絶対値はあまり意味がありません.MacBook Air(13-inch, Early 2015/2.2GHz Core i7, ElCapitan 10.11.6)ですが,softのversionにもよります.

ただし,どんなrandomを生成しているかで結果は大きく変わるので,ちょっと怪しいです. たとえば,MapleのRandomMatrixでgeneratorを指定しないと size=1000でも9.773[sec]ととんでもないぐらい時間がかかります.

MapleではさらにRandomMatrixの生成にも時間がかかります. MapleではNAGのライブラリを使っているので...

課題

  1. 次の行列$A$ の固有値をpythonで求めよ. また,対角化行列$P$ を求めて,対角化せよ. $$ A\, = \, \left[ \begin {array}{ccc} 0&1&-2\\ -3&7&-3 \\ 3&-5&5 \end {array} \right] $$
In [5]:
import numpy as np
A = np.array([[0,1,-2],[-3,7,-3],[3,-5,5]])
l, P = np.linalg.eig( A )
print(l)
print(P)
[ 1.  9.  2.]
[[ -9.04534034e-01   2.29415734e-01  -7.07106781e-01]
 [ -3.01511345e-01   6.88247202e-01   1.11022302e-16]
 [  3.01511345e-01  -6.88247202e-01   7.07106781e-01]]
In [6]:
from pprint import pprint
dA = np.dot(np.linalg.inv(P),np.dot(A,P))
np.set_printoptions(precision=3, suppress=True)
pprint(format(dA))
'[[ 1.  0.  0.]\n [-0.  9. -0.]\n [ 0.  0.  2.]]'

行列のpprint設定イロイロ

In [7]:
from pprint import pprint
dA = np.dot(np.linalg.inv(P),np.dot(A,P))
pprint(format(dA))
'[[ 1.  0.  0.]\n [-0.  9. -0.]\n [ 0.  0.  2.]]'
In [8]:
np.set_printoptions(precision=3)
pprint(dA)
array([[ 1.,  0.,  0.],
       [-0.,  9., -0.],
       [ 0.,  0.,  2.]])
In [9]:
np.set_printoptions(precision=3, suppress=True)
pprint(dA)
array([[ 1.,  0.,  0.],
       [-0.,  9., -0.],
       [ 0.,  0.,  2.]])
In [10]:
np.set_printoptions(formatter={'float': '{: 0.3f}'.format}) 
pprint(dA)
array([[ 1.000,  0.000,  0.000],
       [-0.000,  9.000, -0.000],
       [ 0.000,  0.000,  2.000]])
In [ ]: