はじめに
以前、sinc関数を使うことで四元数と回転ベクトルの相互変換における特異点を無くせるという旨の記事を書きました。
そこでは数学的な特徴にのみ着目して式展開を示したのですが、sinc関数は多くの場合プログラミング言語の標準ライブラリに存在しないため、実際に計算する際には誰かが作成したライブラリを信じて使用するか自分で実装する必要があります。
折角なら自分で作るかーということで、この記事ではsinc関数を使うための一つの方法として、テレスコーピング法で作った多項式を元にして最良近似式を作る方法についてまとめます。
本記事で使用したプログラムは以下のリンクより閲覧できます。
sinc関数と特異点
sinc関数は次式で定義される関数です。
$[ \mathrm{sinc}(x) := \frac{\mathrm{sin}(x)}{x} $]
上式を見ると明らかに$( x=0 $)に特異点を持つように見えるのですが、実はこの特異点は除去することができて以下のような曲線を描きます。
なぜ特異点を無くせるかという点については、sin関数のマクローリン展開が
$[ \begin{equation*} \sin (x) = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \cdots \end{equation*} $]ですので、sinc関数の定義通りこれを$( x $)で除すると
$[ \begin{equation*} \mathrm{sinc} (x) := \frac{\sin (x)}{x} = 1 - \frac{x^2}{3!} + \frac{x^4}{5!} - \frac{x^6}{7!} + \cdots \end{equation*} $]となることから、確かに$( x = 0 $)における特異点を除去できることがわかります（こういうのを可除特異点というらしいです）。
夢を叶えて最良近似式
さて、sinc関数を使うことで$( \mathrm{sin}(x) / x $)をそのまま計算する場合に生じる特異点を回避できることがわかりました。
プログラム上でsinc関数を計算しようと思ったときに、一番手っ取り早いのは先に示したマクローリン展開の式をそのまま実装する方法だと思います。 しかし、マクローリン展開では$( x=0 $)から離れるにつれて誤差が大きくなっていくため、計算範囲を広く取るためには近似式の次数をかなり高く取る必要があります。
次数が上がると計算量も増えますから、工学的視点からは近似式としての精度要求を満たす範囲でできるだけ低い次数の近似式を作りたいところです。 この問題をいい感じに解決してくれるのが、本記事のテーマとなっている最良近似式です。
簡単に説明すると、近似式の計算区間を$( a \leq x \leq b $)としたとき、被近似関数$( f(x) $)に対して
$[ \begin{equation*} \min \left \{ \max_{a \leq x \leq b} | g(x) - f(x) | \right \} \end{equation*} $]となるように近似式$( g(x) $)を定めたものを最良近似式といいます[1]。 また、$( g(x) $)が$( n $)次の最良近似式であるときは計算区間内に$( n+2 $)個の偏差点があり、偏差点の符号は交互に正負の値を取るものとします。
このような定義のもと最良近似式を求めると、誤差曲線（$( g(x) - f(x) $)）は例えば以下のようになります。
図中の青点で示した箇所が偏差点です。 マクローリン展開による近似式と異なり、最良近似式では計算区間内に誤差を分布させることで最大誤差を小さく抑えることができます。
チェビシェフ多項式
最良近似式を作る際には、予め最良近似式に近い特徴を持つ近似式を作っておいて、その近似式の係数を改良することで最良近似式に近づけていきます。
最良近似式の元となる式はどう求めても良いのですが、一般にはチェビシェフ多項式の性質を利用することが多いようです。
$( n $)次のチェビシェフ多項式$( T_n(x) $)は
$[ \begin{equation*} T_n(x) = \cos (n \cos^{-1} (x)) \end{equation*} $]と表され、$( g(x) $)を近似式としたとき
$[ \begin{equation*} g(x) = \sum_{k=0}^n a_k x_k = \sum_{k=0}^n b_k T_k (x) \end{equation*} $]が成り立つようにすることができます。イメージとしてはフーリエ級数展開が近いかもしれません。 チェビシェフ多項式ではコサインの組み合わせで関数を近似しているわけですね。
また、コサインの加法定理を用いるとチェビシェフ多項式は$( x $)のベキを用いて表すことができ、$( n=4 $)までを列挙すると
$[ \begin{align*} T_0(x) &= x^0 \\ T_1(x) &= x^1 \\ T_2(x) &= 2 x^2 - x^0 \\ T_3(x) &= 4 x^3 - 3 x \\ T_4(x) &= 8 x^4 - 8 x^2 + x^0 \end{align*} $]となります。この関係は以下の漸化式で表すことができます。
$[ \begin{equation*} T_{n+1}(x) = 2 x T_n(x) - T_{n-1}(x) \quad (n = 1, 2, 3, ...) \end{equation*} $]一方で$( x $)のベキもチェビシェフ多項式で表現することができ、以下の式で表されます。
$[ \begin{equation*} x^n = \frac{1}{2^{n-1}} \sum_{k=0}^{\lfloor n/2 \rfloor} \binom{n}{k} T_{n - 2k}(x) \end{equation*} $]ただし、$( \lfloor n/2 \rfloor $)は$( n/2 $)の整数部分（床関数）、$( \binom{n}{k} $)は二項係数であり、$( n $)が偶数のときには最後の項$( T_0(x) $)の係数を$( 2 $)で割る必要があります。
$( n=4 $)までを列挙すると以下のようになります。
$[ \begin{align*} x^0 &= T_0(x) \\ x^1 &= T_1(x) \\ x^2 &= \frac{1}{2} T_2(x) + \frac{1}{2} T_0(x) \\ x^3 &= \frac{1}{4} T_3(x) + \frac{3}{4} T_1(x) \\ x^4 &= \frac{1}{8} T_4(x) + \frac{1}{2} T_2(x) + \frac{3}{8} T_0(x) \end{align*} $]チェビシェフ多項式がこういった性質を持つことはわかりましたが、どうやってチェビシェフ多項式の係数を求めるんだという一番重要な問題が残っています。 これをいい感じに解決してくれる方法の一つが、次の項で述べるテレスコーピング法です。
テレスコーピング法
テレスコーピング法とは
前述したように$( x $)のベキとチェビシェフ多項式は互いに変換することができます。
そこで、テレスコーピング法ではマクローリン展開などで求めた冪級数を一度チェビシェフ多項式に変換し、高次の項をいくつかカットしてから再び冪級数に戻します。 こうすることでチェビシェフ多項式の特徴を持った冪級数を簡単に作ることができます。
式変形を見たほうがわかりやすいと思うので、偶関数を表す4次の冪級数をテレスコーピング法で2次に縮小する場合の例を示します。 まずは前述した関係式を用いて冪級数からチェビシェフ多項式に変換し、
とします。次に$( T_4(x) $)の項を切り落としてもう一度冪級数の形に書き直すと、
$[ \begin{align*} g(x)_{N=2} &= \left(a_0 + \frac{a_2}{2} + \frac{3 a_4}{8} \right) T_0(x) + \left(\frac{a_2}{2} + \frac{4 a_4}{8} \right) T_2(x) \\ &= \left(a_0 + \frac{a_2}{2} + \frac{3 a_4}{8} \right) x^0 + \left(\frac{a_2}{2} + \frac{4 a_4}{8} \right) (2x^2 - x^0) \\ &= \left(a_0 - \frac{a_4}{8} \right) x^0 + \left(a_2 + a_4 \right) 2x^2 \end{align*} $]となり、チェビシェフ多項式の性質を持った冪級数が得られました。これがテレスコーピング法です。
冪級数からチェビシェフ多項式への変換
前述したテレスコーピング法は、要するに冪級数の係数を一度チェビシェフ多項式の係数に変換し、高次の項の係数を消してから再度冪級数の係数に変換し直しているわけです。
つまり係数変換さえできれば頑張って式変形する必要はないのですが、残念ながら冪級数の係数をチェビシェフ多項式の係数に（またはその逆に）直接変換する式は無いようです。
直接変換できる式がないのは残念ですが、式変形には規則性があるのでちょっとしたプログラムで自動変換できます。 変換のために、$( T_n(x) $)と$( x^n $)の関係から係数を抜き出して以下のような表を作ります。
| $( T_5(x) $) | $( T_4(x) $) | $( T_3(x) $) | $( T_2(x) $) | $( T_1(x) $) | $( T_0(x) $) | |
|---|---|---|---|---|---|---|
| $( x^0 = $) | $( 1 $) | |||||
| $( x^1 = $) | $( 1 $) | |||||
| $( x^2 = $) | $( \frac{1}{2} $) | $( \frac{1}{2} $) | ||||
| $( x^3 = $) | $( \frac{1}{4} $) | $( \frac{3}{4} $) | ||||
| $( x^4 = $) | $( \frac{1}{8} $) | $( \frac{4}{8} $) | $( \frac{3}{8} $) | |||
| $( x^5 = $) | $( \frac{1}{16} $) | $( \frac{5}{16} $) | $( \frac{10}{16} $) |
先に示したテレスコーピング法の計算例を見てもらえばわかると思いますが、上記の表の行ごとにチェビシェフ多項式の係数$( a $)を掛けて、 同じ列の要素を全て足せば冪級数の係数$( b $)が求まります。
これをプログラムに実装すると以下のようになります（Rustで書いてあります）。
// Nは多項式の最高次数（0次から数えるので，係数の数はN+1個）
/// 階乗（n!）
fn factorial(n: i32) -> f64 {
assert!(n >= 0);
let mut y = 1.0;
for i in 1..=n {
y *= i as f64;
}
y
}
/// 二項係数（nCk）
fn binomial(n: i32, k: i32) -> f64 {
assert!(n >= k);
factorial(n) / ( factorial(k) * factorial(n - k) )
}
/// x^nの多項式係数列aからチェビシェフ多項式Tn(x)の係数列bを求める．
///
/// a: [a0, a1, a2, ..., aN]
fn convert_exp2chebyshev(a: [f64; N+1]) -> [f64; N+1] {
// 多項式の係数をまとめた配列
// x^0 = 1*T_0(x)
// x^1 = 1*T_1(x)
// x^2 = 0.5*T_0(x) 0.5*T_2(x)
// x^3 = 0.75*T_1(x) 0.25*T_3(x)
// ...
let mut coefs = [[0.0; N+1]; N+1];
coefs[0][0] = 1.0;
coefs[1][1] = 1.0;
for n in 2..(N+1) {
let tmp = 2.0f64.powi(n as i32 - 1);
for k in 0..=(n / 2) {
// x^nとT_n(x)の関係式を実装する
let index = n - 2*k;
coefs[n][index] = binomial(n as i32, k as i32) / tmp;
if index == 0 {
coefs[n][index] *= 0.5;
}
}
}
// coefsの行ごとにaをかけて同じ列の要素を全て足せばbが求まる
let mut b = [0.0; N+1];
for i in 0..N {
for j in i..N {
b[i] += a[j] * coefs[j][i];
}
}
b // [b0, b1, b2, ..., bN]
}
チェビシェフ多項式から冪級数への変換
先ほどと同様に係数を抜き出して表にすると以下のようになります。
| $( x^5 $) | $( x^4 $) | $( x^3 $) | $( x^2 $) | $( x^1 $) | $( x^0 $) | |
|---|---|---|---|---|---|---|
| $( T_0(x) = $) | $( 1 $) | |||||
| $( T_1(x) = $) | $( 1 $) | |||||
| $( T_2(x) = $) | $( 2 $) | $( -1 $) | ||||
| $( T_3(x) = $) | $( 4 $) | $( -3 $) | ||||
| $( T_4(x) = $) | $( 8 $) | $( -8 $) | $( 1 $) | |||
| $( T_5(x) = $) | $( 16 $) | $( -20 $) | $( 5 $) |
xのベキからチェビシェフ多項式への変換の際と同様に、上記の表の行ごとにチェビシェフ多項式の係数$( b $)を掛けて、 同じ列の要素を全て足せば冪級数の係数$( a $)が求まります。
試しに3次までの式で計算してみると
$[ \begin{align*} g(x) &= b_0 T_0(x) + b_1 T_1(x) + b_2 T_2(x) + b_3 T_3(x) \\ &= b_0 x^0 + b_1 x^1 + b_2 (2 x^2 - x^0) + b_3 (4 x^3 - 3x) \\ &= (b_0 - b_2) x^0 + (b_1 - 3 b_3) x^1 + 2 b_2 x^2 + 4 b_3 x^3 \end{align*} $]となり、確かに表を用いて変換できることがわかります。
この変換をプログラムに実装すると以下のようになります（こちらもRustで書いてあります）。
// Nは多項式の最高次数（0次から数えるので，係数の数はN+1個）
/// チェビシェフ多項式Tn(x)の係数列bからx^nの多項式係数列aを求める
///
/// b: [b0, b1, b2, ..., bN]
fn convert_chebyshev2exp(b: [f64; N+1]) -> [f64; N+1] {
// 多項式の係数をまとめた配列
// T0(x) = 1
// T1(x) = 1*x^1
// T2(x) = -1 2*x^2
// T3(x) = -3*x^1 4*x^3
// ...
let mut coefs = [[0.0; N+1]; N+1];
coefs[0][0] = 1.0;
coefs[1][1] = 1.0;
for i in 2..(N+1) {
for j in 0..i {
// T_n+1(x) = 2*x*T_n(x) - T_n-1(x)
coefs[i][j+1] = 2.0 * coefs[i-1][j];
coefs[i][j] -= coefs[i-2][j];
}
}
// coefsの行ごとにbをかけて同じ列の要素を全て足せばsが求まる
let mut a = [0.0; N+1];
for i in 0..(N+1) {
for j in i..(N+1) {
a[i] += b[j] * coefs[j][i];
}
}
a // [a0, a1, a2, ..., aN]
}
テレスコーピング法の計算例
試しに20次までのマクローリン展開で近似した多項式をテレスコーピング法で12次の多項式に縮小したところ、次の係数が得られました（偶関数なので奇数次の項は省略しています）。
| 記号 | 値 |
|---|---|
| $( s_0 $) | 0.9999999999999999 |
| $( s_2 $) | -0.16666666666665766 |
| $( s_4 $) | 0.008333333333188874 |
| $( s_6 $) | -0.00019841269754547076 |
| $( s_8 $) | 2.7557294426482064e-6 |
| $( s_{10} $) | -2.5048466629256012e-8 |
| $( s_{12} $) | 1.579349121783874e-10 |
$( s_n $)は$( n $)次の項の係数です。この係数を用いた近似式をグラフにすると以下のようになります。
上段が関数の出力、下段が真値との誤差を表しています。 また、凡例『maclaurin』が14次のマクローリン展開の結果、凡例『chebyshev』がテレスコーピング法で求めた12次の多項式の結果を表します。
ここで、真値をどのように計算するかですが、Rustの組み込み関数として提供されるsin関数は十分誤差が小さいものとして、大部分の値は$( \sin (x) / x $)で計算し、0付近のみ20次のマクローリン展開で計算することにしました。
話を先程のグラフに戻すと、マクローリン展開で求めた近似式はテレスコーピング法で求めた近似式よりも最高次数が2次高い（偶関数なので実質的には1次違うようなものですが）にも関わらず、計算区間両端における誤差はテレスコーピング法の誤差よりも随分大きいことがわかります。
テレスコーピング法で求めた近似式は計算区間全体に誤差を分布させることで最大誤差を小さく抑えているのですね。
というか、テレスコーピング方で求めた式は既に十分誤差が小さいし最良近似式っぽい形をしているような…。 テレスコーピング法で求めた近似式の誤差のみを拡大して表示すると以下のようになります。
誤差曲線の両端が偏差点になっていないのと、偏差点によって最大誤差が1e-16だったり2e-16だったりするので厳密には最良近似式ではありません。 しかし、倍精度で表現できる限界のところで値がふらついているようなのでこれ以上係数を改良してもあまり意味は無いでしょう。
すべての場合にここまで上手く行くとは限りませんが、テレスコーピング方でこれだけ特性の良い近似式を作ることができるのですね。
近似式の改良
計算の制約と下準備
前述したように、テレスコーピング法は上手く使えば最良近似式に近い近似式を作ることができます。 先の例の場合にはあれを最終的な近似式としても実用上問題ないと思いますが、せっかくなので係数を改良して最良近似式を作ってみましょう。
ただし、係数を改良して最良近似式とする場合には近似式の要求精度よりも多い桁数で計算する必要があるため、ここでは近似式の要求精度を1e-6程度まで落として最良近似式を計算します。
まずはテレスコーピング法で最良近似式の元となる近似式を作ります。 8次までのマクローリン展開をテレスコーピング法で4次の多項式に縮小したところ、以下の係数が得られました。
| 記号 | 値 |
|---|---|
| $( s_0 $) | 0.9999937996031747 |
| $( s_2 $) | -0.16655505952380953 |
| $( s_4 $) | 0.008035714285714285 |
先ほどと同様に$( s_n $)は$( n $)次の項の係数です。この係数を用いた近似式をグラフにすると以下のようになります。
凡例『maclaurin』が6次のマクローリン展開の結果、凡例『chebyshev』がテレスコーピング法で求めた4次の多項式の結果を表します。
上図では見やすさのために6次のマクローリン展開の結果と比較していますが、テレスコーピング法の結果と同じ4次で計算すると、マクローリン展開で求めた近似式の最大誤差はテレスコーピング法で求めた近似式の２倍近くになります。
さて、テレスコーピング法で求めた近似式の偏差点の数を数えてみると、全部で7個あることがわかります。あれ…？最良近似式の偏差点は$( n+2 $)個になるはずだったような…？
ここは私もよくわかっていないのですが、どうにも近似式の次数と偏差点の個数が合わないのですよね。 最高次数に係数が0の項を含めて5次の近似式として扱えば最良近似式の計算に持ち込むことができ、後ほど示すようにこれで問題なく計算できるものの、本当にそれで良いのでしょうか。
偶関数だからたまたま最高次数の項が0になってるだけなんですかね。
計算法
元となる近似式が得られたので、係数を改良するための計算法を見ていきましょう。
まず、被近似関数$( f(x) $)の近似式$( g(x) $)が任意に定めることのできる$( n+1 $)個の係数$( c_0 $)〜$( c_n $)を用いて
$[ \begin{align*} g(x) = c_0 x^0 + c_1 x^1 + c_2 x^2 + \cdots + c_n x^n \end{align*} $]と表されるものとします。また、誤差関数$( E(x) $)を
$[ \begin{align*} E(x) = g(x) - f(x) \end{align*} $]とします。
導出の詳細は文献[1]を見ていただくとして、全部で$( n+2 $)個あるうちの$( j $)番目の偏差点を$( x_j (j=0, 1, …, n+1) $)とすると、 以下の連立一次方程式を解くことで最良近似式に近づけるための補正量$( \Delta c_i $)を求めることができます。これを元の係数$( c_i $)に足して新しい係数$( c_i (i = 0, 1, \cdots, n) $)とします。
偏微分のところは一見複雑に見えますが、誤差関数の定義より
$[ \begin{align*} \frac{\partial E(x_j)}{\partial c_i} = x_j^i \quad (i = 0, 1, \cdots, n) \end{align*} $]となるので簡単に計算できます。
ここで偏差点をどう求めるかですが、初回の計算ステップのみチェビシェフ多項式の性質を利用して以下の式より求めます。
$[ \begin{align*} x_j = \cos \frac{ j \pi }{ n + 1 } \quad (j=0, 1, ..., n+1) \end{align*} $]2回目の計算ステップからは多項式の係数を改良したことで偏差点の位置が変わってくるので、次のオイラー法で新しい偏差点を探します（右辺の式で左辺の値を更新しています）。
$[ \begin{align*} x_j \leftarrow x_j - \frac{E'(x_j)}{E''(x_j)} \end{align*} $]文献[1]の実装例では4次精度の微分[3]を用いていたので、ここでもそれに倣うことにすると
となります。また、計算点は一つの山を8分割したうちの偏差点の周り5点を使うことにして$( h = 0.25 / n $)とします（間隔が広すぎなければ適当でいいと思います）。
結局、計算の流れとしては以下のようになります。
- 初期の偏差点$( x_j $)を計算する。
- 連立方程式を解いて$( \Delta c_i $)を求め、係数$( c_i $)を更新する。
- オイラー法で偏差点$( x_j $)を更新する。
- 偏差点の揃い方を見て収束判定。揃っていなかったら2〜4を繰り返す。
偏差点の揃い方は次のようにして調べます。
$( j $)番目の偏差点における誤差を$( \varepsilon_j $)としたとき、全ての偏差点における$( | \varepsilon_j | $)の最大値を$( \varepsilon_\mathrm{max} $)として
$[ \begin{align*} \varepsilon_\mathrm{min} = \min_{j = 0, \cdots, n+1} \frac{ | \varepsilon_j | }{ \varepsilon_\mathrm{max} } \end{align*} $]を求めます。係数を改良して最良近似式に近づくと$( \varepsilon_\mathrm{min} $)は$( 1 $)に収束していくので、各計算ステップの最後で$( \varepsilon_\mathrm{min} $)を調べて適当なところで計算を終えます。
計算結果
テレスコーピング法で求めた5次の近似式（計算のために最高次数に0の項を含めています）の係数を上述した方法で改良し、最良近似式を作ってみました。
反復計算は$( \varepsilon_\mathrm{min} > 0.999999 $)となって時点で打切ることにして計算したところ、以下の係数が得られました。 比較のため左の列に改良前の値、右の列に改良後の値を表示しています。
| 記号 | 改良前の値 | 改良後の値 |
|---|---|---|
| $( s_0 $) | 0.9999937996031747 | 0.999993983687626 |
| $( s_2 $) | -0.16655505952380953 | -0.16655775414525806 |
| $( s_4 $) | 0.008035714285714285 | 0.008040771577902516 |
これをグラフにして確認すると以下のようになります。
上段はsinc関数の計算結果、下段は誤差を表します。凡例『origin』が改良する前の係数を使用した近似式、『fined』が係数を改良した近似式の計算結果を示します。また、『deviation point』は偏差点です。
係数を改良する前の近似式では偏差点における誤差が不揃いですが、係数を改良した後の近似式では偏差点における誤差が均一になり最良近似式となっています。 係数は小数点以下4桁目辺りまで同じですが、誤差曲線で見るとちゃんと改良されていることがわかりますね。
おわりに
この記事では、テレスコーピング法で作った多項式を元にしてsinc関数の最良近似式を作る方法についてまとめました。
マクローリン展開で求めた近似式も、そこからテレスコーピング法で求めた近似式も、更に係数を改良して求めた最良近似式も、全て式の形は同じなのに係数が少し違うだけで誤差曲線が全く異なってくるというのは面白いですね。
今回は$( x=0 $)における特異点を回避することが目的だったので別に最良近似式じゃなくても良かったわけですが、マクローリン展開よりも最大誤差が小さく計算量も少ない近似式を作る方法があるということを知っておくと何処かで役に立つと思います。
計算機性能の向上により計算コストをそこまで気にしなくなってきている昨今においても、マイコンなどの限られた計算資源を活用する場合にはこういった知識が解決の糸口になるかもしれません。