contents

8.3

テキストp.262の表8.2より、共分散行列Vを作成する。 ただし、分散は、不偏分散を用いる。

x1 = (11.91, 18.37, 3.64, 24.37, 30.42, -1.45, 20.11, 9.28, 17.63, 15.71)
x2 = (29.59, 15.25, 3.53, 17.67, 12.74, -2.56, 25.46, 6.92, 9.73, 25.09)
他のxkも同様に設定する。

E(x1) = (11.91 + 18.37 + ... + 15.71) / 10
= 149.99 / 10
= 15.00
E(x2) = (29.59 + 15.25 + ... + 25.09)
= 143.42 / 10
= 14.34
他のE(xk)も同様に求める。

cov(x1, x1) = var(x1)
= E[(x1 - E(x1))2]
= [ (11.91 - 15)2 + (18.37 - 15)2 + ... + (15.71 - 15)2] / (10 - 1)
= [9.54 + 11.36 + ... + 0.50] / 9
= 812.34 / 9
= 90.26
cov(x1, x2) = E[(x1 - E(x1))(x2 - E(x2))]
= E[ (11.91 - 15.00)(29.59 - 14.34) + (18.37 - 15.00)(15.25 - 14.34) + ... + (15.71 - 15.00)(25.09 - 14.34)]
=[-47.12 + 3.07 + ... + 7.63 ] / (10 - 1)
= 458.01 / 9
= 50.89
他のcov(xi, xj)も同様に求める。 

    [ 90.26  50.89  79.02 40.18 ]
V = [ 50.89 107.23 105.38 30.99 ]
    [ 79.02 105.38 162.20 56.54 ]
    [ 40.18  30.99  56.54 68.25 ]

続いて、共分散行列Vの固有値と固有ベクトルを累乗法(power method)で求める。

出発ベクトルx0を次とし、xk+1 = Vxkを反復する。 

x0
=[ 1 ]
 [ 0 ]
 [ 0 ]
 [ 0 ]

x1 = Vx0
= [ 90.26 ] = 260.35[ 0.347 ]
  [ 50.89 ]         [ 0.195 ]
  [ 79.02 ]         [ 0.303 ]
  [ 40.18 ]         [ 0.154 ]

x1
 = [ 0.347 ]と置き直す。以下同様。
   [ 0.195 ]
   [ 0.303 ]
   [ 0.154 ]

x2 = Vx1
= [  71.42 ] = 300.42[ 0.238 ]
  [  75.37 ]         [ 0.251 ]
  [ 105.95 ]         [ 0.353 ]
  [  47.68 ]         [ 0.159 ]

x3 = Vx2
= [  68.47 ] = 309.05[ 0.222 ]
  [  81.08 ]         [ 0.262 ]
  [ 111.40 ]         [ 0.360 ]
  [  48.10 ]         [ 0.156 ]

x4 = Vx3
= [  68.09 ] = 310.76[ 0.219 ]
  [  82.22 ]         [ 0.265 ]
  [ 112.42 ]         [ 0.362 ]
  [  48.04 ]         [ 0.155 ]

x5 = Vx4
= [  68.04  ] = 310.09[ 0.219 ]
  [  82.43  ]         [ 0.265 ]
  [ 112.61  ]         [ 0.362 ]
  [  48.01  ]         [ 0.154 ]

x6 = Vx5
= [  68.03  ] = 311.14[ 0.219 ]
  [  82.47  ]         [ 0.265 ]
  [ 112.64  ]         [ 0.362 ]
  [  48.00  ]         [ 0.154 ]

x7 = Vx6
= [  68.03  ] = 311.16[ 0.219 ]
  [  82.48  ]         [ 0.265 ]
  [ 112.65  ]         [ 0.362 ]
  [  48.00  ]         [ 0.154 ]

x8 = Vx7
= [  68.03  ]
  [  82.48  ]
  [ 112.65  ]
  [  48.00  ]

収束したので、 次のレイレイ(Rayleigh)商により固有値λ（近似値）を算定する。

λ = xkTxk+1 / |xk|2 = xkTAxk / xkTxk

上式に、次を代入する。 

x7
= [ 0.219 ]
  [ 0.265 ]
  [ 0.362 ]
  [ 0.154 ]

Vx7 
= [  68.03 ]
  [  82.48 ]
  [ 112.65 ]
  [  48.00 ]

x7TVx7 = 84.927
x7Tx7 = 0.27293

∴ λ = 84.9266 / 0.272932 = 311.16

以上により、
固有値は、311.16 である。
固有ベクトルは、(0.219, 0.265, 0.362, 0.154) である。

累乗法

累乗法は、次を用いた。

累乗法 出発ベクトルx0を適当に選び

xk+1 = Axk

を反復してx1, x2, x3, ... を作ると、 k → ∞ のとき xkは行列Aの絶対値最大の固有値に対応する固有ベクトルに収束する。

十分に反復を行った後、レイレイ(Rayleigh)商

λ = xkTxk+1 / |xk|2 = xkTAxk / xkTxk

により固有値（近似値）算定する。

出典 戸川隼人、情報処理入門コース7、数値計算、岩波書店、1991、pp.124-125.

Python

added on 2019-09-22 8.3.py

Julia

add on 2023-04-30 8.3.jl

JavaScript

2024-11-02 8.3.js

Fortran

2024-11-02 8.3.f90

history

2017-01-15 create.
2019-09-22 add a Python code.
2021-02-17 change XHTML to html5, change shift-jis to utf-8, add viewport.
2023-04-30 move the Python code to link, add a Julia code.
2024-11-02 add a JavaScript code and a Fortran code.