参考: 三上直樹, (2017), 「知っ得! 軽量sin/cos計算アルゴリズム 第3回 軽量ミニマックス近似式の求め方」, 『Interface(インターフェース)』(CQ出版) 2017年 09 月号, pp.109-111
Interface (インターフェース) 誌に、“軽量” なミニマックス近似多項式の求め方が紹介されているので Python で再現してみる。参考記事と同じく sin(pi * x/2) を -1 ≦ x ≦ 1 の範囲で近似してみる。
Python に置き換えた計算手順は下のとおりである。
- ミニマックス近似多項式に近い多項式を、チェビシェフ多項式 sympy.chebyshevt_poly() を利用して求める。
- scipy.signal.argrelmax() を使って真値とチェビシェフ近似との差の極値を見つける。
- いくつかあるその極値の大きさが全部同じになるまで (すなわちミニマックスになるまで) 多項式の各次数の係数を補正する。
多項式の各次数の補正値については行列を作って聯立方程式を解いて求めるが、正方行列にならないので numpy.linalg.solve() は使わずに擬似逆行列 numpy.linalg.pinv() を利用して解く。
あくまで、全体の精度は低いが (すなわち計算量が少なくて済むが) 区間内での誤差がミニマックスになっている近似式を求めることが目的であるので、この方法で高次まで近似してもあまり意味がない。
`ソースコード (.py) を見る/隠す`
MinMax.py3
# 参考: 三上直樹, (2017), 「知っ得! 軽量sin/cos計算アルゴリズム 第3回 軽量ミニマックス近似式の求め方」, 『Interface(インターフェース)』(CQ出版) 2017年 09 月号, pp.109-111
from matplotlib import pyplot as plt
from scipy.signal import argrelmax
import sympy as sym
import numpy as np
import scipy as sci
class MinMax:
def __init__(self, func, FromTo, N, tol):
self.a, self.b = FromTo
self.func = func
self.N = N
self.tol = tol
self.chebyPoly = self.__cheby(self.func, self.N)
self.c0List = self.chebyPoly.copy()
self.x0List, self.EPS0List = self.__epsilon(self.func, self.c0List)
# チェビシェフ補間による近似多項式を求める函数
def __cheby(self, func, N):
x = sym.Symbol("x")
k = np.arange(N+1)
s = [2 / (N + 1), 2 * k + 1, sci.pi / (2 * (N + 1)), (1 / (N + 1))]
p = s[1] * s[2]
T = [sym.chebyshevt_poly(n, x) for n in k] # 第 1 種チェビシェフ多項式を求める。
r = np.cos(p) # 数列 r を求める。
fVals = func(r) # 数列 r を使って係数 c[0] を求める。
c = np.zeros(N+1)
c[0] = s[3] * np.sum(fVals)
c[1:N+1] = [s[0] * np.sum(fVals * np.cos(p * n)) # 係数 c[1] ～ c[N] を求める。
for n in np.arange(1, N+1)]
g = np.sum([c[n] * T[n] for n in k]) # 係数 c とチェビシェフ多項式 T とで近似式 g(x) を組み立てる。
return np.array([round(g.coeff(x, n), 14) for n in k]) # 0 次 ～ N 次の係数を取り出す。
# 真値と近似式との差の極値を見つける函数
def __epsilon(self, f, coeffs):
dim = len(coeffs)
def g(x) : return sum([coeffs[n] * x**n for n in np.arange(dim)])
def EPS(f, g): return (f - g)
xList = np.linspace(self.a, self.b, 2**16+1) # x 軸を細かく分割して、
EPSList = EPS(f(xList), g(xList)) # 各 x の誤差 EPS を計算して、
peaksIndex = argrelmax(np.abs(EPSList)) # 誤差 EPS が極値となるところの x のインデックスを求めて、
x0List = np.concatenate(([self.a], xList[peaksIndex], [self.b])) # 区間の両端も含めて誤差 EPS が極値となるところの x の値をリスト化する。
EPSa = EPSList[0] # 区間の両端も含めて誤差 EPS の極値をリスト化する。
EPSb = EPSList[-1]
EPS0List = np.concatenate(([EPSa], EPSList[peaksIndex], [EPSb]))
return x0List, EPS0List
# ミニマックス状態にあるかどうかを判定する函数
def isMinMax(self):
vals = np.abs(self.EPS0List)
Max = np.max(vals)
Min = np.min(vals)
return ((Max - Min) < self.tol)
# 多項式を組み立てる函数
def buildSymPoly(self, coeffs):
x = sym.Symbol("x")
dim = len(coeffs)
coeffs = np.round(coeffs, 12)
return sum([coeffs[n] * x**n for n in np.arange(dim)])
# 近似式の係数を補正する函数
def updateCoeffs(self, coeffs):
self.x0List, self.EPS0List = self.__epsilon(self.func, coeffs)
dimX = len(self.x0List)
dimC = len(coeffs)
# 聯立方程式を解くための行列を組み立てるが、
# 正方行列にならないので numpy.linalg.solve() は使えない。
# 擬似逆行列を計算して解く。
A = [np.concatenate(([self.x0List[k]**n for n in range(dimC)],
[np.sign(self.EPS0List[k])]))
for k in np.arange(dimX)]
Ainv = np.linalg.pinv(A) # A X = B の A の擬似逆行列 Ainv を求めて、
B = (self.EPS0List.T).reshape(dimX, 1) # A X = B の B を多行 1 列 の行列に変換して、
sol = Ainv @ B # X = Ainv B を計算して、
sol = np.round(sol, 14) # 本来 0 になるはずなのに計算誤差として残っている項を丸めて排除して、
sol = (sol[0:-1]).flatten() # 多項式の係数の補正値だけを抜き出して 1 行の行列に変換して、
self.c0List = coeffs + sol # 多項式の係数を補正する。
# 計算結果を表示する函数
def showErr(self, chebyPoly, mmPoly):
dim = len(mmPoly)
x_numbers = np.linspace(self.a, self.b, 2**16+1)
# 多項式を組み立てる。
def poly(x, coeffs): return sum([coeffs[n] * x**n for n in np.arange(dim)])
cheby = poly(x_numbers, chebyPoly) # チェビシェフ補間による近似多項式
MinMax = poly(x_numbers, mmPoly) # 係数を補正した多項式
OriginFunc = self.func(x_numbers)
EPSCheby = OriginFunc - cheby # チェビシェフ近似との差
EPSMinMax = OriginFunc - MinMax # 係数を補正した多項式との差
self.maxErrCheby = np.max(np.abs(EPSCheby))
self.maxErrMinMax = np.max(np.abs(EPSMinMax))
plt.grid(color="black")
plt.plot(x_numbers, EPSCheby)
plt.plot(x_numbers, EPSMinMax)
plt.legend(["εCheby", "εMinMax"])
plt.show()
########
# test #
########
if __name__ == "__main__":
np.set_printoptions(precision=14, linewidth=200)
def func(x): return np.sin(np.pi * x/2) # sin(pi * x/2) の近似多項式を求める。
a, b = -1, 1 # 区間。区間は u = (2 * x - (b + a))/(b - a) の式で変換する。
N = 5 # N 次まで近似する。
tol = 1e-14 # 全部の極値の大きさの差が tol 未満になるまで (すなわちミニマックスになるまで) 補正する。
numOfCor = 20 # ただし補正は numOfCor 回で打ち切る。
mm = MinMax(func, (a, b), N, tol)
for i in range(numOfCor+1):
isMinMax = mm.isMinMax()
print("#%s isMinMax? %s" % (i, isMinMax))
print("0 ～ N 次の係数: %s" % (list(mm.c0List)))
print("誤差の極値の位置 x: %s" % (list(mm.x0List)))
print("誤差の極値の大きさ ε: %s" % (list(mm.EPS0List)))
if isMinMax: break
mm.updateCoeffs(mm.c0List)
print("")
print("")
mm.showErr(mm.chebyPoly, mm.c0List)
print("チェビシェフ近似多項式: %s" % (mm.buildSymPoly(mm.chebyPoly)))
print("ミニマックス近似多項式: %s" % (mm.buildSymPoly(mm.c0List)))
print("")
print("真値とチェビシェフ近似との差の最大値: %0.9f" % (mm.maxErrCheby))
print("真値とミニマックス近似との差の最大値: %0.9f" % (mm.maxErrMinMax))
#0 isMinMax? False
0 ～ N 次の係数: [0.0, 1.5706573558985499, -0.0, -0.64345777331469001, 0.0, 0.072934648358349993]
誤差の極値の位置 x: [-1.0, -0.873260498046875, -0.535888671875, -0.1424560546875, 0.1424560546875, 0.535888671875, 0.873260498046875, 1.0]
誤差の極値の大きさ ε: [0.00013423094220987863, -0.00011772771359030987, 7.140718466347451e-05, -1.294246151678502e-05, 1.294246151678502e-05, -7.140718466347451e-05, 0.00011772771359030987, -0.00013423094220987863]
#1 isMinMax? False
0 ～ N 次の係数: [0.0, 1.57031993870501, 0.0, -0.64204567341571006, 0.0, 0.0717827295452]
誤差の極値の位置 x: [-1.0, -0.873260498046875, -0.535888671875, -0.1424560546875, 0.1424560546875, 0.535888671875, 0.873260498046875, 1.0]
誤差の極値の大きさ ε: [0.00013423094220987863, -0.00011772771359030987, 7.140718466347451e-05, -1.294246151678502e-05, 1.294246151678502e-05, -7.140718466347451e-05, 0.00011772771359030987, -0.00013423094220987863]
#2 isMinMax? False
0 ～ N 次の係数: [0.0, 1.57032045080772, 0.0, -0.64211499896079005, 0.0, 0.071862167810559999]
誤差の極値の位置 x: [-1.0, -0.905487060546875, -0.625152587890625, -0.218994140625, 0.218994140625, 0.625152587890625, 0.905487060546875, 1.0]
誤差の極値の大きさ ε: [5.6994834499946023e-05, -6.4969990688279466e-05, 7.6652020768008811e-05, -6.7043717586368068e-05, 6.7043717586368068e-05, -7.6652020768008811e-05, 6.4969990688279466e-05, -5.6994834499946023e-05]
#3 isMinMax? False
0 ～ N 次の係数: [0.0, 1.57032001920971, 0.0, -0.64211316694128007, 0.0, 0.071860854119300002]
誤差の極値の位置 x: [-1.0, -0.900115966796875, -0.6214599609375, -0.221435546875, 0.221435546875, 0.6214599609375, 0.900115966796875, 1.0]
誤差の極値の大きさ ε: [6.7619657490025631e-05, -6.7877736576016368e-05, 6.7656670622362469e-05, -6.7630008892760607e-05, 6.7630008892760607e-05, -6.7656670622362469e-05, 6.7877736576016368e-05, -6.7619657490025631e-05]
#4 isMinMax? False
0 ～ N 次の係数: [0.0, 1.5703200191564399, 0.0, -0.64211316698839005, 0.0, 0.071860854234140001]
誤差の極値の位置 x: [-1.0, -0.900115966796875, -0.62158203125, -0.221466064453125, 0.221466064453125, 0.62158203125, 0.900115966796875, 1.0]
誤差の極値の大きさ ε: [6.7706387729904449e-05, -6.7706387732124895e-05, 6.7706435955217081e-05, -6.77063899374164e-05, 6.77063899374164e-05, -6.7706435955217081e-05, 6.7706387732124895e-05, -6.7706387729904449e-05]
#5 isMinMax? True
0 ～ N 次の係数: [0.0, 1.5703200191564399, 0.0, -0.64211316698839005, 0.0, 0.071860854234140001]
誤差の極値の位置 x: [-1.0, -0.900115966796875, -0.62158203125, -0.221466064453125, 0.221466064453125, 0.62158203125, 0.900115966796875, 1.0]
誤差の極値の大きさ ε: [6.7706402189893211e-05, -6.7706402182232672e-05, 6.7706402185341297e-05, -6.7706402185452319e-05, 6.7706402185452319e-05, -6.7706402185341297e-05, 6.7706402182232672e-05, -6.7706402189893211e-05]
チェビシェフ近似多項式: 0.072934648358*x**5 - 0.643457773315*x**3 + 1.570657355899*x
ミニマックス近似多項式: 0.071860854234*x**5 - 0.642113166988*x**3 + 1.570320019156*x
真値とチェビシェフ近似との差の最大値: 0.000134231
真値とミニマックス近似との差の最大値: 0.000067706
求まった多項式を実装するときは c1 x + c3 x^3 + c5 x^5 = x (c1 + x (x (c3 + c5 x^2))) のように x でどんどん括って計算量を減らす。
水色が真値とチェビシェフ近似との差。
オレンジが真値とミニマックス近似との差。
チェビシェフ近似は中心付近は精度が高いが区間両端に行くほど誤差が大きくなる。
ミニマックス近似は区間全体での誤差が一定の範囲に収まっている。
最初のチェビシェフ近似だけで十分な気もする。
関連: http://ti-nspire.hatenablog.com/archive/category/%E3%83%9F%E3%83%8B%E3%83%9E%E3%83%83%E3%82%AF%E3%82%B9%E5%A4%9A%E9%A0%85%E5%BC%8F%E8%BF%91%E4%BC%BC