ミニマックス近似多項式

MinMax.py

# 参考: 三上直樹, (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))