ONLY DO WHAT ONLY YOU CAN DO

こけたら立ちなはれ 立ったら歩きなはれ

Go 言語で固有値を求める (ヤコビ法)

f:id:fornext1119:20150711180736p:plain
対称行列 f:id:fornext1119:20150805175030p:plain の非対角成分のうち絶対値が最大の成分を f:id:fornext1119:20150805175330p:plain とするとき,
f:id:fornext1119:20150805175341p:plain
その他の対角成分を f:id:fornext1119:20150805175723p:plain, その他の非対角成分を f:id:fornext1119:20150805175731p:plain とした行列 f:id:fornext1119:20150805175119p:plain を考える.
f:id:fornext1119:20150805175405p:plain

f:id:fornext1119:20150805175030p:plainf:id:fornext1119:20150805175119p:plain と, 転置行列 f:id:fornext1119:20150805175211p:plain を左右からかけてできた行列
f:id:fornext1119:20150805175710p:plain
固有値固有ベクトルは元の行列 f:id:fornext1119:20150805175030p:plain から変化しない. これを相似変換という.

f:id:fornext1119:20150807123424p:plain となるように f:id:fornext1119:20150805175429p:plain の値を設定して
処理を反復すると, 固有値が対角成分に並んだ対角行列に収束する.
f:id:fornext1119:20150806113934p:plain
また、f:id:fornext1119:20150807123856p:plain の列ベクトルが固有ベクトルになる.

package main

import "fmt"
import "math"

const N = 4

// ヤコビ法で固有値を求める
func main() {
    var a [N][N]float64 = [N][N]float64{{5.0, 4.0, 1.0, 1.0},
                                        {4.0, 5.0, 1.0, 1.0},
                                        {1.0, 1.0, 4.0, 2.0},
                                        {1.0, 1.0, 2.0, 4.0}}
    var v [N][N]float64 = [N][N]float64{{1.0 ,0.0 ,0.0 ,0.0},
                                        {0.0 ,1.0 ,0.0 ,0.0},
                                        {0.0 ,0.0 ,1.0 ,0.0},
                                        {0.0 ,0.0 ,0.0 ,1.0}}

    // ヤコビ法
    jacobi(&a, &v)

    fmt.Println("\neigenvalue")
    disp_eigenvalue(a)

    fmt.Println("\neigenvector")
    disp_eigenvector(v)
}

// ヤコビ法
func jacobi(a *[N][N]float64, v *[N][N]float64) {
    for k := 1; k < 100; k++ {
        // 最大値を探す
        var p int = 0
        var q int = 0
        var max_val float64 = 0.0
        for i := 0; i < N; i++ {
            for j := i + 1; j < N; j++ {
                if (max_val < math.Abs(a[i][j])) {
                    max_val = math.Abs(a[i][j])
                    p = i
                    q = j
                }
            }
        }

        // θ を求める
        var t float64 = 0.0
        if (math.Abs(a[p][p] - a[q][q]) < 0.00000000001) {
            // a_{pp} = a_{qq} のとき、回転角θをπ/4にする
            t = math.Pi / 4.0
            if (a[p][p] < 0) {
                t = -t
            }
        } else {
            // a_{pp} ≠ a_{qq} のとき
            t = math.Atan(2.0 * a[p][q] / (a[p][p] - a[q][q])) / 2.0
        }

        // θ を使って 行列 P を作成し、A = P^t × A × P
        var c float64 = math.Cos(t)
        var s float64 = math.Sin(t)
        // P^t × A
        var t1 float64 = 0.0
        var t2 float64 = 0.0
        for i := 0; i < N; i++ {
            t1      =  a[p][i] * c + a[q][i] * s
            t2      = -a[p][i] * s + a[q][i] * c
            a[p][i] = t1
            a[q][i] = t2
            // 固有ベクトル
            t1      =  v[p][i] * c + v[q][i] * s
            t2      = -v[p][i] * s + v[q][i] * c
            v[p][i] = t1
            v[q][i] = t2
        }
        // A × P
        for i := 0; i < N; i++ {
            t1      =  a[i][p] * c + a[i][q] * s
            t2      = -a[i][p] * s + a[i][q] * c
            a[i][p] = t1
            a[i][q] = t2
        }

        // 対角要素を表示
        fmt.Printf("%3d\t", k)
        disp_eigenvalue(*a)

        // 収束判定
        if (max_val < 0.00000000001) {
            break
        }
    }
}

// 対角要素を表示
func disp_eigenvalue(a[N][N]float64) {
    for i := 0; i < N; i++ {
        fmt.Printf("%14.10f\t", a[i][i])
    }
    fmt.Println("")
}
// 固有ベクトルを表示
func disp_eigenvector(matrix[N][N]float64) {
    for i := 0; i < N; i++ {
        var row = []float64 {0.0 ,0.0 ,0.0 ,0.0}
        for j := 0; j < N; j++ {
            row[j] = matrix[i][j]
        } 
        normarize(row)
        disp_vector(row)
    }
}
// 1次元配列を表示
func disp_vector(row[]float64) {
    for _, col := range row {
        fmt.Printf("%14.10f\t", col)
    }
    fmt.Println("")
}
// 正規化 (ベクトルの長さを1にする)
func normarize(x[]float64) {
    var s float64 = 0.0

    for i := 0; i < N; i++ {
        s += x[i] * x[i]
    }
    s = math.Sqrt(s)
    
    for i := 0; i < N; i++ {
        x[i] /= s
    }
}
z:\go>8g GO1105.go

z:\go>8l -o GO1105.exe GO1105.8

z:\go>GO1105
  1       9.0000000000    1.0000000000    4.0000000000    4.0000000000
  2       9.0000000000    1.0000000000    6.0000000000    2.0000000000
  3      10.0000000000    1.0000000000    5.0000000000    2.0000000000
  4      10.0000000000    1.0000000000    5.0000000000    2.0000000000

eigenvalue
 10.0000000000    1.0000000000    5.0000000000    2.0000000000

eigenvector
  0.6324555320    0.6324555320    0.3162277660    0.3162277660
 -0.7071067812    0.7071067812   -0.0000000000   -0.0000000000
 -0.3162277660   -0.3162277660    0.6324555320    0.6324555320
  0.0000000000    0.0000000000   -0.7071067812    0.7071067812
参考文献