Remezのアルゴリズム

Remezのアルゴリズム（英: Remez algorithm）、Remez法またはRemezの交換アルゴリズム（英: Remez exchange algorithm）とはEvgeny Yakovlevich Remez（英語版）によって1934年に発表された、関数に簡単な近似関数を見つけるために使用される反復アルゴリズムであり、具体的には、関数をチェビシェフ空間で一様ノルム L_∞を最適化することで求める。^[1]

チェビシェフ空間の典型的な例は、 区間 ${\displaystyle C=[a,b]}$上の実連続関数空間における次数${\displaystyle n}$の チェビシェフ多項式の部分空間であり、与えられた部分空間内の最良近似の多項式は、多項式と関数の間の最大絶対差を最小にするものと定義される。 この場合、解の形式は等振動定理（英語版）によって正確性が保証される。

手順

まず近似する関数fに対し、近似区間中の${\displaystyle n+2}$のサンプル点の集合${\displaystyle X=\{x_{1},x_{2},...,x_{n+2}\}}$を考える。${\displaystyle X}$は通常、区間中に線形に射影されたチェビシェフ多項式の極値を用いる。

1. 未知数${\displaystyle b_{0},b_{1}...b_{n}}$および${\displaystyle E}$について以下の線形連立方程式を解く

   ${\displaystyle b_{0}+b_{1}x_{i}+...+b_{n}x_{i}^{n}+(-1)^{i}E=f(x_{i})}$ （${\displaystyle i=1,2,...n+2}$ ）

2. ${\displaystyle b_{i}}$を多項式${\displaystyle P_{n}(x)=\sum b_{i}x^{i}}$を形成する係数とする。
3. 絶対誤差 ${\displaystyle |P_{n}(x)-f(x)|}$が極大となる点${\displaystyle m}$の集合${\displaystyle M}$を見つける。
4. 誤差がすべての${\displaystyle m}$について大きさが等しく、符号が交互である場合、 ${\displaystyle P_{n}}$は、最大誤差最小の近似多項式となる。 そうでない場合は、 ${\displaystyle X}$を${\displaystyle M}$に置き換え、上記の手順を繰り返す。
5. 上記の手順を繰り返す。

結果は、最良近似多項式またはミニマックス近似（英語版）と呼ばれる。

Remezアルゴリズムの実装における技術のレビューはW. Fraserによってなされた。 ^[2]

初期値の選択

多項式補間の理論における役割のため、近似の初期値としてチェビシェフノードが一般的に使用される。 ラグランジュ補間${\displaystyle L_{n}(f)}$を用いた関数${\displaystyle f}$の最適化問題の初期化では、この近似の初期値が次の範囲にあることが示される。

${\displaystyle \lVert f-L_{n}(f)\rVert _{\infty }\leq (1+\lVert L_{n}\rVert _{\infty })\inf _{p\in P_{n}}\lVert f-p\rVert }$

ここで、ノード（${\displaystyle t_{1},\ldots ,t_{n+1}}$）のラグランジュ補間演算子${\displaystyle L_{n}}$のノルム（ルベーグ定数（英語版））は

${\displaystyle \lVert L_{n}\rVert _{\infty }={\overline {\Lambda }}_{n}(T)=\max _{-1\leq x\leq 1}\lambda _{n}(T;x)}$

と表される。

${\displaystyle T}$はチェビシェフ多項式の零点であり、ルベーグ関数は

${\displaystyle \lambda _{n}(T;x)=\sum _{j=1}^{n+1}\left|l_{j}(x)\right|,\quad l_{j}(x)=\prod _{\stackrel {i=1}{i\neq j}}^{n+1}{\frac {(x-t_{i})}{(t_{j}-t_{i})}}}$

となる。

Theodore A. Kilgore、 ^[3] Carl de Boor、およびAllan Pinkus ^[4]は、各${\displaystyle L_{n}}$に一意の${\displaystyle t_{i}}$が存在することを証明したが、（通常の）多項式については明示的に知られていない。 同様に、 ${\displaystyle {\underline {\Lambda }}_{n}(T)=\min _{-1\leq x\leq 1}\lambda _{n}(T;x)}$ 、およびノードの選択の最適性は${\displaystyle {\overline {\Lambda }}_{n}-{\underline {\Lambda }}_{n}\geq 0}$のように表現できる。

準最適ではあるが分析的に明示的な選択肢を提供するチェビシェフノードの場合、漸近的挙動は以下の様になることが知られている^[5]

${\displaystyle {\overline {\Lambda }}_{n}(T)={\frac {2}{\pi }}\log(n+1)+{\frac {2}{\pi }}\left(\gamma +\log {\frac {8}{\pi }}\right)+\alpha _{n+1}}$

（ γはオイラー・マスケローニ定数 ）

ここで、

${\displaystyle 0<\alpha _{n}<{\frac {\pi }{72n^{2}}}}$for ${\displaystyle n\geq 1,}$

であり、および上限^[6]

${\displaystyle {\overline {\Lambda }}_{n}(T)\leq {\frac {2}{\pi }}\log(n+1)+1}$

があたえられる。

Lev Brutman ^[7]は${\displaystyle n\geq 3}$ かつ${\displaystyle {\hat {T}}}$が拡張されたチェビシェフ多項式の零点であるとき以下を示した。

${\displaystyle {\overline {\Lambda }}_{n}({\hat {T}})-{\underline {\Lambda }}_{n}({\hat {T}})<{\overline {\Lambda }}_{3}-{\frac {1}{6}}\cot {\frac {\pi }{8}}+{\frac {\pi }{64}}{\frac {1}{\sin ^{2}(3\pi /16)}}-{\frac {2}{\pi }}(\gamma -\log \pi )\approx 0.201}$

また、Rüdiger Günttner ^[8]は、 ${\displaystyle n\geq 40}$について

${\displaystyle {\overline {\Lambda }}_{n}({\hat {T}})-{\underline {\Lambda }}_{n}({\hat {T}})<0.0196}$

を示した。

詳細な議論

この節では、上記の手順の詳細について説明する。${\displaystyle i}$ は${\displaystyle 0\leq i\leq n+1}$である。

ステップ1

与えられた${\displaystyle x_{0},x_{1},...x_{n+1}}$に対して 、以下の線形連立方程式を未知数${\displaystyle b_{0},b_{1}...b_{n}}$および${\displaystyle E}$について解く

${\displaystyle b_{0}+b_{1}x_{i}+...+b_{n}x_{i}^{n}+(-1)^{i}E=f(x_{i})}$

${\displaystyle (-1)^{i}E}$の項が存在しているためノード${\displaystyle x_{0},...,x_{n+1}}$ は整列（単調増加もしくは単調減少）している必要があることがわかる。 このとき、この線形方程式には一意の解が存在する。また、標準的なライブラリのソルバーで線形方程式を解く場合は${\displaystyle O(n^{3})}$の計算量となるが、この処理は${\displaystyle O(n^{2})}$の計算量で得られることがわかっている。

以下に簡単な証明を次に示す。

${\displaystyle f(x)}$に対して${\displaystyle n}$次補間関数${\displaystyle p_{1}(x)}$を、${\displaystyle (-1)^{i}}$について${\displaystyle n}$次補間関数${\displaystyle p_{2}(x)}$を最初の${\displaystyle n+1}$ノード（${\displaystyle x_{0},...,x_{n}}$）について計算する。

${\displaystyle p_{1}(x_{i})=f(x_{i}),p_{2}(x_{i})=(-1)^{i},i=0,...,n.}$

この計算のため、各補間関数の計算において${\displaystyle 0,...,n}$階の差商を使用したニュートン補間を用いるとそれぞれ${\displaystyle O(n^{2})}$の計算量となる。

多項式${\displaystyle p_{2}(x)}$の${\displaystyle x_{i-1}}$と${\displaystyle x_{i}}$の間に${\displaystyle i}$番目の零点がある（ ${\displaystyle i=1,\ldots ,n}$）。したがって、${\displaystyle x_{n}}$と${\displaystyle x_{n+1}}$の間に零点はなく、${\displaystyle p_{2}(x_{n})}$と${\displaystyle p_{2}(x_{n+1})}$は${\displaystyle (-1)^{n}}$ と同符号である。

線形結合${\displaystyle p(x):=p_{1}(x)-p_{2}(x)\!\cdot \!E}$は、${\displaystyle n}$次の多項式であり、

${\displaystyle p(x_{i})=p_{1}(x_{i})-p_{2}(x_{i})\!\cdot \!E\ =\ f(x_{i})-(-1)^{i}E,\ \ \ \ i=0,\ldots ,n.}$

と変形できる。

これは任意の${\displaystyle E}$に対して、${\displaystyle i=0,\ldots ,n}$において式${\displaystyle b_{0}+b_{1}x_{i}+...+b_{n}x_{i}^{n}+(-1)^{i}E}$と同じである。${\displaystyle i=n+1}$とのときは

${\displaystyle p(x_{n+1})\ =\ p_{1}(x_{n+1})-p_{2}(x_{n+1})\!\cdot \!E\ =\ f(x_{n+1})-(-1)^{n+1}E}$

となり、${\displaystyle E}$について以下の式を得ることができる。

${\displaystyle E\ :=\ {\frac {p_{1}(x_{n+1})-f(x_{n+1})}{p_{2}(x_{n+1})+(-1)^{n}}}.}$

上述の通り、分母の2項は同じ符号を持つため、${\displaystyle E}$及び ${\displaystyle p(x)\equiv b_{0}+b_{1}x+\ldots +b_{n}x^{n}}$常に一意に定まる。

与えられた${\displaystyle n+2}$の整列したノードでの誤差は、以下の通り正負の交互の順になる。

${\displaystyle p(x_{i})-f(x_{i})\ =\ -(-1)^{i}E,\ \ i=0,...,n\!+\!1.}$

de La Vallée Poussinの定理によれば、この条件下では誤差が${\displaystyle E}$より小さい${\displaystyle n}$次の多項式は存在しない。もしそのような多項式${\displaystyle {\tilde {p}}(x)}$が存在した場合、${\displaystyle p(x)-{\tilde {p}}(x)=(p(x)-f(x))-({\tilde {p}}(x)-f(x))}$はノード${\displaystyle x_{i}}$において正負交互のままであり${\displaystyle n+1}$個の零点を有することになるが、次数${\displaystyle n}$の多項式ではありえない。 したがって、この${\displaystyle E}$は、次数${\displaystyle n}$の多項式におけるで達成できる誤差の下限となる。

ステップ2

${\displaystyle p(x)}$を${\displaystyle b_{0}+b_{1}x+...+b_{n}x^{n}}$とする。

ステップ3

次のように入力ノード${\displaystyle x_{0},...,x_{n+1}}$とそのエラー${\displaystyle \pm E}$を更新する。

${\displaystyle p(x)-f(x)}$について、零点によって区切られる正負の領域にそれぞれついて、正の領域で${\displaystyle x_{i}}$を極大値${\displaystyle {\bar {x}}_{i}}$に置き換え、負の領域では${\displaystyle x_{i}}$極小値${\displaystyle {\bar {x}}_{i}}$に置き換える。 （期待する${\displaystyle {\bar {x}}_{0}}$ Aで${\displaystyle {\bar {x}}_{i}}$近く${\displaystyle x_{i}}$ 、そして${\displaystyle {\bar {x}}_{n+1}}$ Bで 。）このステップでは高い精度は必要なく2つの2次近似を使用した直線探索で十分である。 （ ^[9]）

ここで${\displaystyle z_{i}:=p({\bar {x}}_{i})-f({\bar {x}}_{i})}$ とする。 各${\displaystyle z_{i}}$に対して絶対値${\displaystyle |z_{i}|}$はE以上となる。 de La Vallée Poussinの定理とその証明は、 ${\displaystyle z_{0},...,z_{n+1}}$についての${\displaystyle \min\{|z_{i}|\}\geq E}$についても適用可能で、次数${\displaystyle n}$の最良多項式近似の最大誤差の下限と考えることができる。

また、 ${\displaystyle \max\{|z_{i}|\}}$は最大誤差の上限となる。

ステップ4

${\displaystyle \min \,\{|z_{i}|\}}$と${\displaystyle \max \,\{|z_{i}|\}}$を最良多項式近似の最大誤差の下限と上限として、十分な精度、すなわち ${\displaystyle \max\{|z_{i}|\}-\min\{|z_{i}|\}}$十分に小さいか、減少しなくなったら繰り返しを終了する。 これらの上下限は反復処理の収束状況を示していると考えられる。

派生

以下のような派生が有る。

• 複数のサンプル点を、同時に最大絶対差の近くの場所に置き換える。
• すべてのサンプル点を、すべての位置、交互の符号、最大差となるよう単一の反復で置き換える。 ^[10]
• 浮動小数点演算を使用するコンピューターで関数を計算するために近似を使用する場合は、 相対誤差を使用して近似と関数の差を定義することがある。
• 修正されたRemez交換アルゴリズムには、無誤差点制約が含まれることがある。 ^[10]

関連項目

• 近似法

参考文献