Pythonで機械学習を学ぶための前知識として線形代数を勉強中。
理解を深めるためにも、アルゴリズムを自分で組んでみることにしました。
折角なのでGo言語でプログラミング。
今回やりたいこと
- n元連立1次方程式の解を求める。 n元はスペックの許す限り無限。
- 数学的ライブラリや組み込み関数の利用は禁止。ランダム、四則演算のみ可能。
- 今回は正則行列のみを対象。
- ロジックが簡単そうな掃き出し法（ガウス・ジョルダン法）を実装する。
- 行列積を求めるロジックを作り、2の逆行列と行列の積が単位行列になることを確認する
- パフォーマンスは2の次。やっていることが明確なロジックにすること
今後のこの方面（数学的プログラミング）の学習予定
- 正則でない行列の判定を行い「解なし」or「不定」の表現をできるようにする
- 余因子展開（クラメール）、LU分解、共役勾配法をプログラミング
- 値の小さい2次行列の固有値・固有ベクトルを求め、xyとXYの基底関係をplotで表現してみる
- 行列式の公式をプログラミングし、n元の行列式を出せるようにする（サラスは実装済）
- 統計学にも手を広げてみる
- 微積もやっとくか…
- Go言語でも3次元を表現できるplotが欲しいです安西先生
前提条件
数学的知識を発信する記事ではありません。
間違っている場合、ご指摘があればなるべく早く訂正しますのでご了承ください。
コーディング開始
fox/main.go
自作関数を呼んでるだけのロジック。
fox/main.go
package main
import "fox/matrix"
func main() {
simultaneousLinearEquations(6)
}
// n元連立1次方程式を解くための関数。
// 引数でn次元の数を指定すると、ランダムな整数値を持った正方行列を生成し解を求める。
func simultaneousLinearEquations(unknownsNum int) {
// 左辺をランダムに生成
lfsMatrix := matrix.CreateRamdomMatrix([]int{unknownsNum, unknownsNum})
// 左辺をランダムに生成
rfsMatrix := matrix.CreateRamdomMatrix([]int{unknownsNum, 1})
// ガウス・ジョルダンによる解を求める
ge := matrix.NewGaussJordanElimination(lfsMatrix, rfsMatrix)
ge.Solve()
printEquations(lfsMatrix, rfsMatrix)
printAnswer(lfsMatrix, ge.RfsMatrix)
printInverse(lfsMatrix, ge.InverseMatrix)
}
fox/matrix/matrix.go
行列を扱うための関数群。今回必要なもののみを用意。
多次元配列をコピーする関数が無くて衝撃を受けるなど。
ランダム行列の生成とか、DeepCopyのロジックとか、処理量がorderの乗数になっててかなりイケていない気がしたけど今回の本題ではないので愚直に実装。
あとGoって演算子の再定義できないんですね…これはかなり悲しい。
fox/matrix/matrix.go
package matrix
import (
"math/rand"
"time"
)
// Multiplication 行列積を求める関数。(関数名Dotは微妙に意味が違う気がした）
// もちろん左の列数と右の行数が一致しない行列を渡すとpanic
func Multiplication(leftMatrix [][]float64, rightMatrix [][]float64) [][]float64 {
matrix := make([][]float64, len(leftMatrix))
for i := range leftMatrix {
matrix[i] = make([]float64, len(rightMatrix[0]))
for n := range rightMatrix[0] {
matrix[i][n] = 0
for m := range rightMatrix {
matrix[i][n] = matrix[i][n] + leftMatrix[i][m]*rightMatrix[m][n]
}
}
}
return matrix
}
// DeepCopy 行列のディープコピーを行う関数
// ありそうでなかった…（ちゃんと探せばあるのかも）
func DeepCopy(matrix [][]float64) [][]float64 {
newMatrix := make([][]float64, len(matrix), cap(matrix))
for i, item := range matrix {
newMatrix[i] = make([]float64, len(item), cap(item))
copy(newMatrix[i], item)
}
return newMatrix
}
// CreateRamdomMatrix 指定された次元の正方行列をランダム生成する関数
// かなり愚直なロジック…。
func CreateRamdomMatrix(index []int) [][]float64 {
matrix := make([][]float64, index[0])
rand.Seed(time.Now().UnixNano())
// 生成が1ナノミリ秒以内に終わってしまい同一シードでの生成が行われてしまうためスリープ
time.Sleep(1)
for i := range matrix {
matrix[i] = make([]float64, index[1])
for n := range matrix[i] {
if rand.Intn(2) == 1 {
matrix[i][n] = float64(rand.Intn(100))
} else {
matrix[i][n] = float64(rand.Intn(100)) * -1
}
}
}
return matrix
}
// CreateUnitMatrix 指定された次元の単位行列を生成する関数
func CreateUnitMatrix(index int) [][]float64 {
unitMatrix := make([][]float64, index)
for i := range unitMatrix {
unitMatrix[i] = make([]float64, index)
unitMatrix[i][i] = 1
}
return unitMatrix
}
// ReplaceMatrixRow 行列の行を入れ替える関数
func ReplaceMatrixRow(matrix [][]float64, index1 int, index2 int) {
if index1 == index2 {
return
}
tmp := matrix[index1]
matrix[index1] = matrix[index2]
matrix[index2] = tmp
}
fox/matrix/gaussian.go
今回のメインとなるパッケージ。
掃き出し法自体が非常に単純な処理ではあるんですが、左辺や右辺を間違えていたり、添え字を間違っていたりと無駄にハマったのである程度処理単位でメソッドを分けました。
自分的には分かりやすいつもりです笑
テストケースもあるのでそのうちテストケース回の記事を書くかもしれません（たぶん書かない）
fox/matrix/gaussian.go
/*
gaussianパッケージは正方行列に対してガウス・ジョルダン法の適用を行うライブラリです。
行列式 |A|=0 となるような行列には対応しておりません。
[使い方]
lfsMatrix := [][]float64{{2, 5}, {1, 3}}
rfsMatrix := [][]float64{{3},{2}}
ge := gaussian.NewGaussJordanElimination(lfsMatrix, rfsMatrix)
result := ge.Solve()
*/
package matrix
// GaussJordanElimination ガウス・ジョルダン法（掃き出し法）を行う構造体
type GaussJordanElimination struct {
LfsMatrix [][]float64 // 左辺の行列
RfsMatrix [][]float64 // 右辺の行列
InverseMatrix [][]float64 // 逆行列
}
// NewGaussJordanElimination GaussJordanEliminationの構造体を生成するためのコンストラクタ。
func NewGaussJordanElimination(lfsMatrix [][]float64, rfsMatrix [][]float64) *GaussJordanElimination {
ge := new(GaussJordanElimination)
ge.InverseMatrix = CreateUnitMatrix(len(lfsMatrix))
ge.LfsMatrix = DeepCopy(lfsMatrix)
ge.RfsMatrix = DeepCopy(rfsMatrix)
return ge
}
// Solve ガウスジョルダン法を実行するメソッド
func (ge *GaussJordanElimination) Solve() [][]float64 {
for n := range ge.LfsMatrix {
ge.pivoting(n)
ge.diagonalComponents(n)
ge.forwardElimination(n)
}
return ge.RfsMatrix
}
// ピボット選択を行うメソッド
// 選択基準は「対角要素以降の行で最大の値を持つ行を選択」
func (ge *GaussJordanElimination) pivoting(n int) {
pivot := n
for i := n; i < len(ge.LfsMatrix); i++ {
if ge.LfsMatrix[i][n] != 0 && ge.LfsMatrix[pivot][n] < ge.LfsMatrix[i][n] {
pivot = i
}
}
ReplaceMatrixRow(ge.LfsMatrix, pivot, n)
ReplaceMatrixRow(ge.RfsMatrix, pivot, n)
ReplaceMatrixRow(ge.InverseMatrix, pivot, n)
}
// 対角成分の単位ベクトル化と、その行の除算処理を行う関数
func (ge *GaussJordanElimination) diagonalComponents(n int) {
denominator := ge.LfsMatrix[n][n]
for i := range ge.LfsMatrix[n] {
ge.LfsMatrix[n][i] = ge.LfsMatrix[n][i] / denominator
ge.InverseMatrix[n][i] = ge.InverseMatrix[n][i] / denominator
}
ge.RfsMatrix[n][0] = ge.RfsMatrix[n][0] / denominator
}
// 前進消去を行うメソッド。ガウスの消去法ではないので後進のメソッドはない。
func (ge *GaussJordanElimination) forwardElimination(n int) {
for rowIndex := range ge.LfsMatrix {
if rowIndex == n || ge.LfsMatrix[rowIndex][n] == 0 {
continue
}
sign := float64(-1)
magnification := ge.LfsMatrix[rowIndex][n]
for colIndex := range ge.LfsMatrix[rowIndex] {
if colIndex >= n {
ge.LfsMatrix[rowIndex][colIndex] += ge.LfsMatrix[n][colIndex] * magnification * sign
}
ge.InverseMatrix[rowIndex][colIndex] += ge.InverseMatrix[n][colIndex] * magnification * sign
}
ge.RfsMatrix[rowIndex][0] += ge.RfsMatrix[n][0] * magnification * sign
}
}
fox/print.go
結果の確認をするためのprint処理群です。
fox/print.go
package main
import (
"fmt"
"fox/matrix"
)
var eqVarNames = []string{"x", "y", "z", "a", "b", "c", "d", "e", "f", "g",
"h", "i", "j", "k", "m", "n", "o", "p", "q", "r", "s", "t", "u", "v", "w"}
func printEquations(lfs [][]float64, rfs [][]float64) {
fmt.Printf("■今回与えられる %d元連立一次方程式\n", len(lfs))
for i, eq := range lfs {
for n, elm := range eq {
fmt.Printf("%4.0f%s", elm, eqVarNames[n%len(eqVarNames)])
if n < len(eq)-1 {
fmt.Printf(" +")
}
}
fmt.Printf(" = %4.0f\n", rfs[i][0])
}
}
func printAnswer(lfs [][]float64, answer [][]float64) {
fmt.Printf("\n■ガウス・ジョルダン法（掃き出し法）を利用し解を求める\n")
for i, eq := range answer {
fmt.Printf(" %s =%12.9f\n", eqVarNames[i%len(eqVarNames)], eq[0])
}
fmt.Printf("\n■与式に代入し確認をしてみる\n")
for i := range lfs {
fmt.Printf(" ")
ans := float64(0)
for n := range lfs[i] {
fmt.Printf("(%3.0f * %.8f)", lfs[i][n], answer[n][0])
ans += lfs[i][n] * answer[n][0]
if n < len(lfs[i])-1 {
fmt.Printf(" + ")
}
}
fmt.Printf(" = %7.3f\n", ans)
}
}
func printInverse(lfs [][]float64, inv [][]float64) {
fmt.Printf("\n■ガウス・ジョルダン法により求めれた逆行列\n")
for i := range inv {
for n := range inv[i] {
fmt.Printf(" %12.9f", inv[i][n])
if n < len(inv[i])-1 {
fmt.Printf(" ")
}
}
fmt.Printf("\n")
}
fmt.Printf("\n■逆行列・行列の結果が単位行列になることを確認\n")
unit := matrix.Multiplication(inv, lfs)
for i := range unit {
for n := range unit[i] {
v := unit[i][n]
if v < 0.00000001 {
v = 0
} else if v > 0.9999999 && v < 1.00000001 {
v = 1
}
fmt.Printf(" %2.2f", v)
}
fmt.Printf("\n")
}
fmt.Printf(" ※最後の単位行列は丸め誤差が生じてしまうため、0と1との差が0.00000001以内のものをそれぞれ0、1と判断")
}
実行
今回もEclipse上でコーディング、実行を行いました。
（早くもemacsにしたくなっているところ…）
行列は任意の行列を与えることも可能ですが、面倒なので行列の次元数を指定すればその次元の正方行列を自動生成して毎回ランダムで与えられた方程式を解くようにしてあります。
とりあえず6元を指定して実行
コンソールのログ
■今回与えられる 6元連立一次方程式
-17x + -61y + 46z + -31a + 60b + 69c = 57
-73x + -79y + 89z + -29a + 84b + -89c = -93
-50x + 97y + -76z + 99a + -73b + 79c = 23
-17x + -70y + -75z + -51a + -55b + -58c = 44
-2x + 95y + -48z + 74a + 98b + 3c = 51
-75x + -24y + -4z + 43a + -65b + -10c = 43
■ガウス・ジョルダン法（掃き出し法）を利用し解を求める
x = 1.887769457
y =-3.022774844
z =-1.099671032
a = 3.289637186
b = 0.453206328
c = 0.435855883
■与式に代入し確認をしてみる
(-17 * 1.88776946) + (-61 * -3.02277484) + ( 46 * -1.09967103) + (-31 * 3.28963719) + ( 60 * 0.45320633) + ( 69 * 0.43585588) = 57.000
(-73 * 1.88776946) + (-79 * -3.02277484) + ( 89 * -1.09967103) + (-29 * 3.28963719) + ( 84 * 0.45320633) + (-89 * 0.43585588) = -93.000
(-50 * 1.88776946) + ( 97 * -3.02277484) + (-76 * -1.09967103) + ( 99 * 3.28963719) + (-73 * 0.45320633) + ( 79 * 0.43585588) = 23.000
(-17 * 1.88776946) + (-70 * -3.02277484) + (-75 * -1.09967103) + (-51 * 3.28963719) + (-55 * 0.45320633) + (-58 * 0.43585588) = 44.000
( -2 * 1.88776946) + ( 95 * -3.02277484) + (-48 * -1.09967103) + ( 74 * 3.28963719) + ( 98 * 0.45320633) + ( 3 * 0.43585588) = 51.000
(-75 * 1.88776946) + (-24 * -3.02277484) + ( -4 * -1.09967103) + ( 43 * 3.28963719) + (-65 * 0.45320633) + (-10 * 0.43585588) = 43.000
■ガウス・ジョルダン法により求めれた逆行列
-0.000108434 -0.015973602 -0.019576388 -0.004692837 0.007260210 0.016159913
-0.007959774 0.014968416 0.020698025 -0.000243264 -0.010079301 -0.026239805
-0.002876637 0.003037332 0.000157517 -0.008081124 -0.005545227 -0.000429715
0.002783137 -0.018545382 -0.024584853 -0.006964579 0.015173789 0.034983900
0.003985050 0.000706373 -0.001945663 0.001452416 0.006005228 -0.000783059
0.007132027 -0.001673686 0.004016588 -0.000375835 -0.001829849 -0.002531279
■逆行列・行列の結果が単位行列になることを確認
1.00 0.00 0.00 0.00 0.00 0.00
0.00 1.00 0.00 0.00 0.00 0.00
0.00 0.00 1.00 0.00 0.00 0.00
0.00 0.00 0.00 1.00 0.00 0.00
0.00 0.00 0.00 0.00 1.00 0.00
0.00 0.00 0.00 0.00 0.00 1.00
※最後の単位行列は丸め誤差が生じてしまうため、0と1との差が0.00000001以内のものをそれぞれ0、1と判断
問題なく結果がでました。
100元連立や1000元連立を指定しても（正則であれば）まともに解が返って来るので問題はなさそうですが、与式のバラつきが大きかったり少数以下が大きいとどうしても誤差が発生してしまいます。
まあでも、アルゴリズムを組んだことである程度の水準で理解はできたつもりなので、今回の目的は十分に果たせたかな。。
～ではまた次回～
次回は数学的記事ではなくGAEの話題の予定です。