2次フィッティングをPythonからGoへ移植
プログラム言語によって、数学ライブラリが充実したり
そうでなかったりするのは何故でしょうか？
実はプログラム言語仕様に「演算子のオーバーロード」が可能か否か
によって数値計算時におけるプログラムの読みやすさが大きく変わってきます。
「演算子のオーバーロード」に対応していない場合、関数コールを繰り返すことにより数式表現を行う必要があります。
演算子を用いた自然な表記には敵わず、オープンソースで活動されているライブラリ開発者のモチベーションも変わってくるでしょう。
AI案件が増えてきている昨今、「Pythonからある数式ライブラリを処理スピードに有利なプログラム言語への移植したい」などのプロジェクトも必然的に増えます。
※Python+GPUの組み合わせが必ず問題解決にはならない。
ある数値計算処理を「Python」から「Go」に移植するなどの業務案件の機会も増えるでしょう。
例えばPythonのnumpyライブラリに2次フィッティングを、実現している処理があります。
np.polyfit(x,y,2)
業務要件にて「Go」言語で同等の機能を実現したいとします。
ネット上の自作している人も見かけない場合、ゼロから自力で実装することになります。
Go言語にて2次フィッティングを実装してみます。
package main
import (
"fmt"
"math"
)
type Element struct {
x float64
y float64
}
func gauss(p_A [][]float64, dimension int) []float64 {
for i := 0; i < dimension; i++ {
l_m := 0.0
l_pivot := i
for l := i; l < dimension; l++ {
if l_m < math.Abs(p_A[l][i]) {
l_m = math.Abs(p_A[l][i])
l_pivot = l
}
}
if l_pivot != i {
l_b := 0.0
for j := 0; j < dimension+1; j++ {
l_b = p_A[i][j]
p_A[i][j] = p_A[l_pivot][j]
p_A[l_pivot][j] = l_b
}
}
}
for k := 0; k < dimension; k++ {
l_p := p_A[k][k]
p_A[k][k] = 1
for j := k + 1; j < dimension+1; j++ {
p_A[k][j] /= l_p
}
for i := k + 1; i < dimension; i++ {
l_q := p_A[i][k]
for j := k + 1; j < dimension+1; j++ {
p_A[i][j] -= l_q * p_A[k][j]
}
p_A[i][k] = 0
}
}
var l_Result []float64 = make([]float64, dimension)
for i := dimension - 1; 0 <= i; i-- {
l_Result[i] = p_A[i][dimension]
for j := dimension - 1; i < j; j-- {
l_Result[i] -= p_A[i][j] * l_Result[j]
}
}
return l_Result
}
func polyfit(elementList []Element, dimension int) []float64 {
n := dimension + 1
m := dimension + 2
l_a := make([][]float64, n)
for i := 0; i < n; i++ {
l_a[i] = make([]float64, m)
}
for i := 0; i < n; i++ {
for j := 0; j < n; j++ {
for k := 0; k < len(elementList); k++ {
l_a[i][j] += math.Pow(elementList[k].x, (float64)(i+j))
}
}
}
for i := 0; i < n; i++ {
for k := 0; k < len(elementList); k++ {
l_a[i][n] += math.Pow(elementList[k].x, (float64)(i)) * elementList[k].y
}
}
return gauss(l_a, n)
}
func main() {
x := [...]float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19}
y := [...]float64{-1, -3, -1, 9, 21, 30, 37, 39, 67, 65, 95, 123, 142, 173, 191, 216, 256, 292, 328, 358}
var elementList []Element = make([]Element, len(x))
for i := 0; i < len(x); i++ {
elementList[i].x = x[i]
elementList[i].y = y[i]
}
result := polyfit(elementList, 2)
a := result[2]
b := result[1]
c := result[0]
fmt.Printf("a = %v\n", a)
fmt.Printf("b = %v\n", b)
fmt.Printf("c = %v\n", c)
x_2 := [...]float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25}
for i := 0; i < len(x_2); i++ {
fmt.Printf("x = %v / y = %v\n", i, math.Pow(a * (float64)(i),2) + (b * (float64)(i)) + c )
}
}
■実行結果
a = 1.0281954887218046
b = -0.5033834586466206
c = -0.34999999999996767
x = 0 / y = -0.34999999999996767
x = 1 / y = 0.20380250438128233
x = 2 / y = 2.8719769348182735
x = 3 / y = 7.654523291311005
x = 4 / y = 14.55144157385948
x = 5 / y = 23.562731782463693
x = 6 / y = 34.68839391712365
x = 7 / y = 47.92842797783935
x = 8 / y = 63.28283396461079
x = 9 / y = 80.75161187743795
x = 10 / y = 100.33476171632088
x = 11 / y = 122.03228348125954
x = 12 / y = 145.84417717225392
x = 13 / y = 171.7704427893041
x = 14 / y = 199.81108033240997
x = 15 / y = 229.9660898015716
x = 16 / y = 262.235471196789
x = 17 / y = 296.61922451806214
x = 18 / y = 333.1173497653909
x = 19 / y = 371.72984693877544
x = 20 / y = 412.4567160382158
x = 21 / y = 455.297957063712
x = 22 / y = 500.2535700152637
x = 23 / y = 547.3235548928711
x = 24 / y = 596.5079116965345
x = 25 / y = 647.8066404262536
ふ～
Pythonでは1行で済んだコードも
ゼロから自作すると、それなりの規模なります。
元ソースの調査分析の方法としては、Pythonのnumpyはオープンソースであるため、コードを見ることができます。それを参考にしながら実装することも可能です。
それでは「ax2+bx+c」の結果を可視化してみます。
用途として、ある程度「未来」や「過去」を予測できたり、実データの欠損部を補完したり、滑らかしたり、いろいろ便利ですね。
ただ、線形近似あるため精度はそれなりです。
行列計算の場合、おそらくもっとプログラムが複雑になりますね。
次回はニューラルネットワークの移植などもチャレンジしてみます。
Discussion