Go 言語で固有値を求める (QR分解)

対称行列 を 正規直交行列 と 右三角行列 の
積に分解し,

逆順に掛けたものを新たな として

処理を反復すると, 元の行列の固有値が対角成分に並んだ右三角行列に収束する.

package main

import "fmt"
import "math"

const N = 4

// QR分解で固有値を求める
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 q [N][N]float64 = [N][N]float64{{0.0 ,0.0 ,0.0 ,0.0},
                                        {0.0 ,0.0 ,0.0 ,0.0},
                                        {0.0 ,0.0 ,0.0 ,0.0},
                                        {0.0 ,0.0 ,0.0 ,0.0}}
    var r [N][N]float64 = [N][N]float64{{0.0 ,0.0 ,0.0 ,0.0},
                                        {0.0 ,0.0 ,0.0 ,0.0},
                                        {0.0 ,0.0 ,0.0 ,0.0},
                                        {0.0 ,0.0 ,0.0 ,0.0}}

    for k := 1; k < 200; k++ {
        // QR分解
        decomp(a, &q, &r)
        // 行列の積
        multiply(r, q, &a)
        // 対角要素を表示
        fmt.Printf("%3d\t", k)
        disp_eigenvalue(a)
        
        // 収束判定
        var e float64 = 0.0
        for i := 1; i < N; i++ {
            for j := 0; j < i; j++ {
                e += math.Abs(a[i][j])
            }
        }
        if (e < 0.00000000001) {
            break
        }
    }

    fmt.Println("\neigenvalue")
    disp_eigenvalue(a)
}
// QR分解
func decomp(a[N][N]float64, q *[N][N]float64, r *[N][N]float64) {
    var x = []float64 {0.0 ,0.0 ,0.0 ,0.0}

    for k := 0; k < N; k++ {
        for i := 0; i < N; i++ {
            x[i] = a[i][k]
        }

        for j := 0; j < k; j++ {
            var t float64 = 0.0
            for i := 0; i < N; i++ {
                t += a[i][k] * q[i][j]
            }
            r[j][k] = t
            r[k][j] = 0.0
            for i := 0; i < N; i++ {
                x[i] -= t * q[i][j]
            }
        }

        var s float64 = 0.0
        for i := 0; i < N; i++ {
            s += x[i] * x[i]
        }
        r[k][k] = math.Sqrt(s)
        for i := 0; i < N; i++ {
            q[i][k] = x[i] / r[k][k]
        }
    }
}
// 行列の積
func multiply(a[N][N]float64, b[N][N]float64, c *[N][N]float64) {
    for i := 0; i < N; i++ {
        for j := 0; j < N; j++ {
            var s float64 = 0.0
            for k := 0; k < N; k++ {
                s += a[i][k] * b[k][j]
            }
            c[i][j] = s
        }
    }
}
// 対角要素を表示
func disp_eigenvalue(a[N][N]float64) {
    for i := 0; i < N; i++ {
        fmt.Printf("%14.10f\t", a[i][i])
    }
    fmt.Println("")
}

z:\go>8g GO1104.go

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

z:\go>GO1104
  1       9.6046511628    1.1012311902    4.8997514499    2.3943661972
  2       9.9219788334    1.0010980897    5.0142271522    2.0626959248
  3       9.9805335651    1.0000111820    5.0095550207    2.0099002322
  4       9.9951218389    1.0000001123    5.0033019159    2.0015761328
  5       9.9987795937    1.0000000011    5.0009686042    2.0002518010
  6       9.9996948428    1.0000000000    5.0002648858    2.0000402713
  7       9.9999237072    1.0000000000    5.0000698501    2.0000064427
  8       9.9999809266    1.0000000000    5.0000180426    2.0000010308
  9       9.9999952316    1.0000000000    5.0000046034    2.0000001649
 10       9.9999988079    1.0000000000    5.0000011657    2.0000000264
(中略)
 85      10.0000000000    5.0000000000    2.0000000000    1.0000000000
 86      10.0000000000    5.0000000000    2.0000000000    1.0000000000
 87      10.0000000000    5.0000000000    2.0000000000    1.0000000000
 88      10.0000000000    5.0000000000    2.0000000000    1.0000000000
 89      10.0000000000    5.0000000000    2.0000000000    1.0000000000
 90      10.0000000000    5.0000000000    2.0000000000    1.0000000000
 91      10.0000000000    5.0000000000    2.0000000000    1.0000000000
 92      10.0000000000    5.0000000000    2.0000000000    1.0000000000
 93      10.0000000000    5.0000000000    2.0000000000    1.0000000000
 94      10.0000000000    5.0000000000    2.0000000000    1.0000000000

eigenvalue
 10.0000000000    5.0000000000    2.0000000000    1.0000000000

参考文献