LoginSignup
8
1

More than 3 years have passed since last update.

@ken_sumi1019(すみ けん)

Goで作った数値計算の関数をまとめとく

Posted at

まだ学部二年生なんですけど、ちょっと背伸びをして三年生推奨の授業を何個か取りました。そこで出された課題に行列演算が必要になりました。先生方はRなりpyhtonなり行列演算の関数がすでに豊富にある言語を使用することをおそらく考えてたと思うんですが、ここはあえてGoで書いて行列演算の仕組みを勉強してしまおうとたくさん書いてたら「おっ！Qiitaの記事一つ書けるんじゃね？」と思ったので書きました。

ここからコードを書いていくよ

三角関数、べき乗、絶対値

行列計算に必要なスカラーの計算をmath.goに書きます。

aの平方根は

a - x^2 = 0

この関数を満たすxの値をニュートン法を使って求めます。更新式は

x_{i+1} := \frac{x_i + \frac{x_i}{2}} {2}

三角関数はマクローリン展開をします。
マクローリン展開の公式はこのようなものです。

f(x) = f(0) + f'(0)(x-0) + \frac{f''(0)(x-0)^2}{2!} + ... + \frac{f^{(n)}(0)(x-0)^n}{n!}

sinとcosの微分は

(sin\theta)' = cos\theta \\
(cos\theta)' = sin\theta

ここで、sinやcosに関して0の場合は値が計算できるので近似的に数値を算出できるということを使います。一般式は

sinx = \sum_{i=0}^{i=n} (-1)^{n-1}\frac{x^{2n-1}}{(2n-1)!} \\

cosx = \sum_{i=0}^{i=n} (-1)^n \frac{x^{2n}}{(2n)!}

これを20次の項まで計算してます。

Expは簡単に展開できるので一般式は省きます。30回ループさせました。

math.go
package math


//ルート計算
func Sqrt(n float64) float64{
    x := n / 2
    prev := -1.1
    for i:=0;i<100;i++{
        prev = x
        x = (x + n / x) / 2
        if float32(x) == float32(prev){
            return x
        }
    }
    return x
}

//べき乗
func Pow(n float64,m int) float64{
    if m == 1{
        return n
    }else if m == 0 {
        return 1
    }else if m % 2 == 0{
        k := Pow(n,m/2)
        return  k*k
    }else {
        return n*Pow(n,m-1)
    }

}

//三角関数
//円周率
var pi = 3.141592653589793
//sin
func Sin(x float64) float64 {
    //符号を調節するために後からかける
    sign := 1
    //xが負の時
    if x < 0{
        x *= -1
        sign *= -1
    }
    //xが2piより大きかった時
    if x > 2*pi{
        x = x - 2*pi*float64(int(x/(2*pi)))
    }
    if pi < x && x <= 2*pi {
        x -= pi
        sign *= -1
    }
    if pi/2 < x && x <= pi{
        x = pi - x
    }
    if x == 0{return 0
    } else if x == pi/2 {
        return float64(sign)
    }
    n := 0.0
    k := 1
    xc := x
    I := 1
    for i:=0;i<20;i++{
        n += (float64(k)*xc)/float64(I)
        k *= -1
        xc *= x*x
        I *= (i+1)*2 * ((i+1) * 2 + 1)
    }
    return n * float64(sign)
}
//cos
func Cos(x float64) float64 {
    //符号を調節するために後からかける
    sign := 1
    if x < 0{x *= -1}
    //xが2piより大きかった時
    if x > 2*pi{
        x = x - 2*pi*float64(int(x/(2*pi)))
    }
    //xがpi<x<2pi
    if pi < x && x <= 2*pi {
        x = 2*pi - x
    }
    //pi/2<x<pi
    if pi/2 < x && x <= pi{
        x = pi - x
        sign *= -1
    }
    if x == 0{return float64(sign)
    } else if x == pi/2 {
        return 0
    }
    n := 0.0
    k := 1
    xc := 1.0
    I := 1
    for i:=0;i<20;i++{
        n += (float64(k)*xc)/float64(I)
        k *= -1
        xc *= x*x
        I *= ((i+1)*2-1) * (i+1) * 2
    }
    return n * float64(sign)
}
//tan
func Tan(x float64) float64 {
    return Sin(x)/Cos(x)
}

//Exp
func Exp(x float64) float64 {
    ans := 1.0
    check := 1.0
    ex := x
    for i:=0;i<30;i++{
        ans += ex/check
        check *= float64(i+2)
        ex *= x
    }
    return ans
}

//絶対値
func Abs(x float64) float64 {
    if x < 0{
        return -x
    }
    return x
}

行列演算

行列の演算をします。matrix.goというファイルに記述します。

内積

行列の授業を受けたことがある方ならだれもが知ってると思われる、内積です。左からかけた行列に注目すると内積は各行ごとの計算が独立していることがわかると思います。それを使って並列化しています。これ以降の行列計算は全部行で分割して並列にしました。GPUとかでめちゃめちゃ並列させることができるマシンを使っている場合はさらに右の行列を列で分解することでさらに行列を分解することができます。まあGoでGPUは使えないようなのでやるとしたらCUDAですかね。

matrix.go
//どんだけ並列するかについての変数
var thread = 4

func Multi(a [][]float64,b [][]float64) [][]float64{
    if len(a[0]) != len(b){
        println("ふええ～おにいちゃーん、内積が計算できないよぉ～～")
        os.Exit(1)
    }
    ans := MakeMatrix(len(a),len(b[0]))
    t := int(Min(float64(thread),float64(len(a))))
    finish := make(chan bool)
    f := func(a [][]float64,b [][]float64,ans [][]float64,s int,g int) {
        for i:=s;i<g;i++{
            for j:=0;j<len(b[0]);j++{
                for k:=0;k<len(b);k++{
                    ans[i][j] += a[i][k] * b[k][j]
                }
            }
        }
        finish <- true
    }
    for i:=0;i<t;i++{
        g := (len(a)/t)*(i+1)
        if i == t -1{
            g = len(a)
        }
        go f(a,b,ans,(len(a)/t)*i,g)
    }
    for i:=0;i<t;i++{
        <-finish
    }
    return ans
}

転置

転置です。トランスです。行で分割して転置してます。

//行列の転置
func Trans(a [][]float64) [][]float64{
    ans := MakeMatrix(len(a[0]),len(a))
    t := int(Min(float64(thread),float64(len(a))))
    finish := make(chan bool)
    f := func(a [][]float64,ans [][]float64,s int,g int) {
        for i:=s;i<g;i++{
            for j:=0;j<len(a[0]);j++{
                ans[j][i] = a[i][j]
            }
        }
        finish <- true
    }
    for i:=0;i<t;i++{
        g := (len(a)/t)*(i+1)
        if i == t -1{
            g = len(a)
        }
        go f(a,ans,(len(a)/t)*i,g)
    }
    for i:=0;i<t;i++{
        <-finish
    }
    return ans
}

シュミットの直交化

直交化なるものをします。固有値を求めるのに使います。
グラムシュミットの直交化と具体例
このサイトに書いてあることを素直に実装してあります。もう少し早い方法があるのかもしれません。

//シュミットの直交化法
func Schmidt(vec [][]float64) [][]float64{
    u := make([][]float64,len(vec))
    //uの初期化
    for i:=0;i<len(u);i++{
        u[i] = make([]float64,len(vec[0]))
    }
    for j:=0;j<len(vec[0]);j++{
        //v := make([]float64,len(vec))
        for i:=0;i<len(vec);i++{
            u[i][j] = vec[i][j]
        }
        for i:=0;i<j;i++{
            //ないせき
            inpro := 0.0
            for k:=0;k<len(vec);k++{
                inpro += vec[k][j] * u[k][i]
            }
            for k:=0;k<len(vec);k++{
                u[k][j] -= inpro * u[k][i]
            }
        }
        //おおきさ
        vecSize := 0.0
        for i:=0;i<len(vec);i++{
            vecSize += u[i][j] * u[i][j]
        }
        vecSize /= math.Sqrt(vecSize)
        for i:=0;i<len(vec);i++{
            u[i][j] /= vecSize
        }
    }
    return u
}

固有値、固有ベクトルの算出

QR法

QR分解というものを使って固有値を算出します。正則行列Aは直交行列Qと上三角行列Rに分解ができるものがQR分解です。直交行列の算出に上で作ったシュミットの直交化を用います。たぶん検索するともうちょっとまともな実装が出てくると思います。

まず、Aを直交行列Qと上三角行列Rを用いて下のようにあらわせます。

A_0 = Q_0R_0

ここで、直交行列の性質を使って、

Q_0^{-1} = Q_0^T \\ 
R_0 = Q_0^TA_0  \\

のように変形できます。さらに、別の行列を 

A_1 = R_0Q_0 \\

と定義します。これを変形してやると

A_1 = Q_0^TA_0Q_0

というように変形ができます。式を見ていただければわかるんですが、A1とA0は相似な関係にあるので固有値が一致します。こいつをAが収束するまでやってあげると対角成分に固有値が出てきます。

//QR法（固有値を求めるやつ）
//収束判定を行っていないので適宜回数調節
func QrMethod(A [][]float64) []float64{
    for i:=0;i<20;i++{
        tyokkou := Schmidt(A)
        A = Multi(Multi(Trans(tyokkou),A),tyokkou)
    }
    ans := make([]float64,len(A))
    for i:=0;i<len(A);i++{
        ans[i] = A[i][i]
    }
    return ans
}

収束判定がいまいちわからなかったので固定で20回ループさせてます。まあまあな精度は出ます。

逆反復法

逆反復法で固有ベクトルを求めます。しかし、この方法は事前に固有値を求めないといけないので上で作ったQR法の関数にて固有ベクトルを求めてから実行します。理論は下のページを見てください。めっちゃ詳しく親切に書いてあります。（証明は読者にお任せしようとか書いてありません。）
その３ 固有値と固有ベクトルを数値解析で求める

//第一引数に固有ベクトルを求めたい行列、第二引数に固有値の配列
//たぶん動く
//固有ベクトルを求めるやつ
//逆反復法？とかいうやつ
func GetEigenvector(A [][]float64,value []float64) [][]float64{
    I := make([][]float64,len(A))
    ans := make([][]float64,len(value))
    for i:=0;i<len(A);i++{
        I[i] = make([]float64,len(A))
        ans[i] = make([]float64,len(A))
    }
    for k:=0;k<len(value);k++{
        for i:=0;i<len(A);i++{
            I[i][i] = value[k]
        }
        //fmt.Println("A",A)
        //fmt.Println("I",I)
        B := Gauss(Add(A,ConstMult(-1,I)))
        //fmt.Println(Add(I,ConstMult(-1,A)))
        //fmt.Println(B)
        x := make([][]float64,len(A))
        for i:=0;i<len(x);i++{
            x[i] = make([]float64,1)
            x[i][0] = 1.0
        }
        for i:=0;i<19;i++{
            x = Multi(B,x)
            //fmt.Println(x)
        }
        for i:=0;i<len(x);i++{
            ans[k][i] = x[i][0]
        }
    }
    return ans
}

ヤコビ法

ヤコビ法というものです。対称行列の固有値と固有ベクトルを求めることができます。遅いです。もっと効率的な実装にできるのかもしれません。あんま意味わかってないけど実装したらうまくいったので載せときます。

//ヤコビ法（対称行列の固有値固有ベクトルを求める）
func Jacobi(mat [][]float64) ([]float64,[][]float64) {
    if len(mat) != len(mat[0]){
        println("ふええ～おにいちゃーん、正方行列じゃないよぉ～～")
        os.Exit(1)
    }
    A := MakeMatrix(len(mat),len(mat))
    B := Eye(len(mat),len(mat))
    for i:=0;i<len(A);i++ {
        for j := 0; j < len(A); j++ {
            A[i][j] = mat[i][j]
        }
    }
    for {
        //Aの非対角成分の絶対値最大値と最大の位置を求める
        maxNum := math.Abs(A[0][1])
        max_I := 0
        max_J := 1
        for i:=0;i<len(A);i++{
            for j:=0;j<len(A);j++{
                if i == j{continue}
                if math.Abs(A[i][j]) > maxNum{
                    maxNum = math.Abs(A[i][j])
                    max_I = i
                    max_J = j
                }
            }
        }
        if maxNum < 0.000001{
            break
        }
        //回転角を求める
        theta := (1/(2*math.Tan((2*A[max_I][max_J])/(A[max_I][max_I] - A[max_J][max_J]))))
        P := Eye(len(mat),len(mat))
        P[max_I][max_I] = math.Cos(theta)
        P[max_J][max_J] = P[max_I][max_I]
        P[max_I][max_J] = math.Sin(theta)
        P[max_J][max_I] = -P[max_I][max_J]

        //
        A = Multi(Multi(Trans(P),A),P)
        B = Multi(B,P)
    }
    //固有ベクトルの配列
    vec := make([]float64,len(mat))
    for i:=0;i<len(mat);i++{
        vec[i] = A[i][i]
    }
    return vec,B
}

逆行列を求める

吐き出し法

吐き出し法です。吐き出します。またはガウスの消去法とも言います。名前がカッコいいので変数名は後者の呼び方に合わせてあります。ピポット選択は実装してないので対角成分に0がきたら詰みます。そのうち付け加えます。

func Gauss(m [][]float64) [][]float64{
    ans := Eye(len(m),len(m[0]))
    mat := MakeMatrix(len(m),len(m[0]))
    copy(mat,m)
    for i:=0;i<len(mat);i++{
        for j:=0;j<len(mat);j++{
            if i == j || mat[j][i] == 0{continue}
            n := mat[j][i] / mat[i][i]
            for k:=0;k<len(mat[0]);k++{
                ans[j][k] -= ans[i][k] * n
                mat[j][k] -= mat[i][k] * n
            }
        }
    }
    for i:=0;i<len(mat);i++{
        n := mat[i][i]
        for j:=0;j<len(mat);j++{
            ans[i][j] /= n
        }
    }
    return ans
}

LU分解（未実装）

実装したはいいけど動きを確認していない上に本当に計算量を削減した実装であるかどうかすらも怪しいので今後追加します。

行列式を求める

なんか逆行列と行列式、やることあんま変わらないから同じカテゴリーに入れてしまいました。行列式は行基本変形をしても結果は変わりません。ですので、行基本変形をして上三角行列を作ります。上三角行列は対角成分の積が行列式なので簡単に計算できます。これも、対角成分に0がきたら詰むのでその辺の対策をそのうちします。

//行列式
func Det(mat [][]float64) float64 {
    buf := MakeMatrix(len(mat),len(mat[0]))
    copy(buf,mat)
    for i:=0;i<len(buf);i++{
        for j:=i+1;j<len(buf);j++{
            n := buf[j][i] / buf[i][i]
            for k:=0;k<len(buf[0]);k++{
                buf[j][k] -= buf[i][k] * n
            }
        }
    }
    ans := 1.0
    for i:=0;i<len(buf);i++{
        ans *= buf[i][i]
    }
    return ans
}

最後に

まだまだ実装も理論もガバガバなので精進していきたいと思います。

8
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
8
1