補間 # 今週の話と来週の話は、テキストにないので色々と解説する。 == 前回は、突然の休講だった。 補講を行う予定である。 来週は、出張なので、講義ができない ( しかし、補講は一度しかできないので.. ) 今週は、講義を二回行う 来週は、演習を二回行う == 関数近似 関数の計算 ( i.e. sin(x) ) は、基本演算 ( 四則 ) の 50 倍から 100 倍 の計算が必要 幾つかの関数値を計算し、その間の関数値はまともに計算せずに補間で補う 関数の誤差 word 長 ( 32 bit ) の最後 bit 辺に誤差が入ってもしょうがない 全ての場合で最後の bit だけで済むようにするのは大変 最後から二つめの bit 位までは有りうる 全ての区間で誤差の入る範囲を縮めたい => 最大最小近似 : 最大誤差を最小にする近似 # 商用のライブラリは、大体、最後の 2 bit 位に収めている ## もちろん、出来の悪いものもある 関数近似では、ほぼ、多項式での近似を行い、誤差を収めている 関数によっては、色々工夫が必要 i.e. 平方根などは、多項式では近似せず、ニュートン法で値を求める。 # 平方根の場合は、正規化 ( 1 < x < 4 ) して、値を値を求める。 ## x = y * 4^n => \sqrt{x} = \sqrt{y} * 2^n ## where 1 < y < 4 ニュートン法では、初期値が重要 y の先頭の数 bit を利用して、初期値の table を look up する。 # ニュートン法の繰り返し回数は 4 回程度で収束する ## メモリが多ければ、大きな 初期値 table を用意することにより、高速化できる。 # プログラマの腕により、ステップ数は変る ## 平方根の場合は、x から y * 4^n を計算する部分 # 意外に、面倒なのが、三角関数 三角関数は周期関数 ( \sin{T} = \sin{T + 2\pi} ) 区間 [0,2\pi] の間に正規化を行えばよさそうだが.. # そうは問屋が卸さない \sin{n} = \sin{ 10000*\2pi + x } ( x \in {0,2\pi} とすると、n の上位 bit は、意味がなくなり、本質的な x の bit が殆ど意味がない ( 入力の時点で、n [x] の精度が落ちている )。 # これは、周期関数一般の話 # 関数の実装側としてはどうしようもない # 関数を利用する人間が気をつけないといけない。 # 三角関数の正規化は更に狭めることができる 対称性を利用すれば、[0,\frac{\pi}{2}] ですむ 加法定理を利用すれば、更に半分とか三分の一にできる 区間を狭めるにはコストがかかるがそうすれば近似 # 関数近似の研究者は少い # 戸田 先生、二宮市三 先生、山下真一郎先生、浜田穂積先生 # 関数近似の計算は一度おこなえば、済む # 多数の人間が行う必要がない # 必要だが少数な人数しか不要な分野は、人材育成が難しい 誤差の振舞い 区間の端で大聴く、真中で小さいのが普通 区間の真中で誤差が最大になる場合は、その区間を含まなない区間まで縮めてしまうことが多い ( そのために、速度にしわ寄せがくるので、そこで苦労する )。 == Sin と Cos は殆ど中身が同じ 名前が違っているだけで、同じ所に入る # sin と cos を対にして利用することが多いので、sin x # を計算するときに、同時に cos x も計算してキャッシュ # すると、sin の後の cos の計算を高速化できる \tan{x} は \frac{\sin{x}}{\cos{x}} で対処 \exp{x} や \log{x} なども、同様に、正規化を行って、近似 == 課題 四則と条件判断だけで、\sin{x} を [ -\frac{\pi}{6}, \frac{\pi}{6} ] の 範囲で計算しなさい どう解くか ? # 最大最初誤差まで考えると大変だが、多少誤差がはいってもきにしないのでれば.. テーラー展開を使えばよい。 f(x+h) = f(x) + h f'(x) + \frac{1}{2!} h^2 f''(x) + .. + \frac{1}{n!}h^nf^{(n)}(x) + .. ここで、f = sin, x = 0 として、f(h) を計算すると.. \sin{h} = \sin{0} + h \sin'{0} + \frac{1}{2!} h^2 \sin''{0} + .. = 0 + h + 0 - ..\frac{1}{6} h^3 + .. h が 原点の近くであれば、h^3 は、ほとんど 0 と考えてよい。 => h^3 まで計算すれば十分 # 最大最小誤差にするには.. 条件に併せて、係数を調整する # どのような係数でも条件を満さなければ、更に、項数をふやす 多項式は、和差積だけで計算できるが、商が利用できるのであれば、更に工夫ができる 連分数による近似 多項式/多項式の形の計算 # 工夫の仕方に関しては、浜田先生の本がある # 「初等関数」 (教育出版(シリーズ物)) 一松 信 先生 # 昔は数値表現によって、同じ関数に対する様々な工夫があったが、 # 現在は、IEEE 方式が普通なので、関数の種類や、ハード化がテーマ ## ハードでの実装 ## コーデック ( HP の電卓で採用され、一次流行った )。 ## CORDIC (STL) ## 昔、Intel の CPU (ペンテイアム) の割り算で、特定のパターンだけが間違えた ## 割り算は複雑なので、関数近似と同じアプローチをした ## テーブルデータに誤りがあったため、このような不思議なことがおきた == [課題] \sin{x} の関数近似 ( テーラ展開による ) を計算する h^5 の項位まで \sin{h} = h - \frac{1}{6}h^3 + \frac{1}{120}h^5 区間 [0, \frac{\pi}{6}] で、近似計算と、ライブラリによる計算の比較 ( 近似 値の計算 ) を行う。 [工夫] オーナー法を適用すれば、少し計算量 ( かけ算回数 ) を減らせる \sin{h} = h ( 1 - h^2 (\frac{1}{6} + h^2 \frac{1}{120} ) ) [ポイント] 1. こんなに簡単に関数近似ができる 2. 最大最小近似でないので、0 付近は誤差がないが、0 から離れると誤差が増える # どう係数を弄れば、最大最小近似になるだろうか ??? テーラー展開による関数近似 原点は、(誤差が小さいという意味で..) 特異点 問題を上手に(理論的に..、数値的な変換を行っても、変換の時点で誤差が)変換して、原点付近にもっていってから計算すると茣蓙が少くてすむことが多い ==== # 今週は、二回講義なので、その二回目 数値微分 関数 f から f' を計算したい f が測定データなので、論理的に微分ができない 画像データでエッジの判定などで微分 ( 変化量が多い所 ) が計算したい # 画像処理 : 須藤先生 数学的には、微分可能だが、計算が大変 ( 例えば..非線形のために ) 関数が、プログラム ( による関数近似 ) として与えられている場合 # 最大最小の値が知りたい ( 最適化問題の解決 ) # 微係数が 0 の所 ## ヤコビアン、ハミルトニアンの計算 f' を近似するには ? (二点間の..) 差分を利用する ( 一次近似 [直線近似] ) f'(x) 〜 \frac{f(x+h)-f(x)}{h} # 誤差を減らすために、次数をふやすことを考えることができる # この形は、どこかでみたことが.. f'(x) = \frac{f(x+h)-f(x)}{h} を変形すると.. h f'(x) = f(x+h)-f(x) f(x+h) = f(x) + h f'(x) これはテーラ展開 !! 精度をあげるには.. ( 福井先生のテーマの一つ ) 近似の次数を上る 微分の階数を上る # 問題領域による。x が長さならば dx/dt は速度、d^2x/dt^2 は可束度になる。このような問題を解く場合は、階数をあげる テーラ展開を利用して微分の階数をあげた近似関数を考える 以下の項を無視すると.. f(x+h) = f(x) + h f'(x) + \frac{1}{2!} h^2 f''(x) + | + \frac{1}{3!}h^3f^{(3)}(x) + .. f'(x) = \frac{f(x+h)-f(x)}{h} - \frac{h}{2!}f''(x) これは、一次近似に、( 打ち切り ) 誤差項 ( \frac{h}{2!}f''(x) ) を加えた結果 数学的には、h を小くすれば、誤差を小くできる 数値計算的には h を小くすると f(x+h) と f(x) が近付くので、 f(x+h) - f(x) の計算により、桁落ち ( 有限桁の誤差 : 一般には、「丸め誤差」とよばれているが、正しいとは思えない。ここでは、有限桁での計算が原因になって生じる誤差なので、「有限桁の誤差」と呼ぶ」) が生じるので、h が小くなると誤差が大きくなる。 # 更に、h で割っているので、誤差の拡大が生じる 一概に、h を小くればよいということにならない。 # 本質的に計算が有限桁で行われているための現象 # 無限桁ならば、h が小さい方がよいという意味で単純 h と誤差の関係は、v 字のグラフになる # v の先は、まるまっているだけでなく、誤差のキャンセルが生れたりするので、「がたがた」する 誤差を最小にする h を考えることができる ( 福井先生のお仕事 ) 微分の階数 m, 近似の次数を n とすると h = \frac{m}{m+n} x bit 数 の場合に、誤差が最小になる。 # もちろん、関数によって、多少の影響を受けるが、一般的にはそんなもの m = n = 1, bit 数 24 の時には、 h = 24 x 1/2 = 12 bit 半分の長さ位の h が最適 h を小くすると良いということにならない !! # 関数の対称性によっても、色々と影響を受ける 次数を上るには.. f(x+h) = f(x) + h f'(x) + \frac{1}{2!} h^2 f''(x) + \frac{1}{3!}h^3f^{(3)}(x) + .. f(x-h) = f(x) - h f'(x) + \frac{1}{2!} h^2 f''(x) - \frac{1}{3!}h^3f^{(3)}(x) + .. の両辺を引くと.. 2 hf'(x) = f(x+h) - f(x-h) - \frac{h^3}{3!}f'''(x) f'(x) = \frac{f(x+h) - f(x-h)}{2h} - \frac{h^2}{6}f'''(x) となる。 公式 f'(x) = \frac{f(x+h) - f(x-h)}{2h} は、中間の地点の微分の 近似を両側の点の差分で近似している 中間の値 ( f(x) ) が利用されていない 計算が少くてすむ 誤差項目は h^2 のオーダ ( 二次の近似 ) 精度が良い # 対称性のある公式は、そうでない公式に比べて一次分、精度が良いことが多い 今度は両辺を足す f''(x) = \frac{f(x+h) - 2h f(x) + f(x-h)}{h^2} - \frac{h^2}{12}f^(4)(x) これはもニ階の微分の公式が、二次オーダで求めることができる 対称性質があると、近似の次数が大きくなる # 対称性があるということは良い # カメラのレンズも対称なもの ( ガウスレンズ ) の方が性能がよい == 微分方程式の離散化 ( の問題点は、少くても二つ ) 差分化の影響 : 解空間そのものが変ってしまう可能性がある 所々しか見ていない サンプルポイントの取りかたにより非どい結果になることがある # sin で、0 になる点を取ってしまうと... ## 差分化の結果は、元の微分方程式の解を含むがそうでない解 ( 幻影解 ) も含まれる ( 判定条件が必要になる )。 有限桁の影響 : 誤差が入ってしまう可能性がある == [課題] 微分の近似を行う 近似公式 f'(x) = \frac{f(x+h)-f(x)}{h} を利用し、 f(x) = \exp{x} の微分値を近似計算する。 x としては、1, \frac{1}{2}, \frac{1}{4}, .. \frac{1}{2^20} の結果を、真値 ( ライブラリ関数の結果 ) との比較 ( 誤差 ) を視る # 計算機や、コンパイラーによって誤差の出方の挙動が異るなる可能性がある ## 数値実験をする場合は、bit の長さを自分でコントロールする ( double で計算して、単精度になるように、下位ビットを 0 clear して、確実に 24 bit で計算する )。 # 打ち切り誤差の影響は明確。しかし、計算 == 「来年度は、科目名も変るし (?) 教科書が変るので、(単位を)落さないように(そうしないと、教科書を新しく買う羽目に..)」 == 数値計算で問題を解くときに、微分で解くより積分でといた方がよい 微分は、変化を拡大するので、不安定になりやすい 積分は、変化を均すので、安定に計算できる フーリエ変換 周期関数なので、特定の窓 ( 周期で.. ) 観測することになる 本来、区間の両側は同じでなければならないのに、そうでないと誤差が大変 三角波を加えることにより、両端の値を同じにして計算 ( 結果が安定する ) フーリエ変換が終ったあとで、三角波の影響を取り除けばよい 数値計算では、連続性が非常に重要 微分は、連続性を壊す 積分は、連続性を作る == 偏微分方程式を解く場合は、h (メッシュサイズ) を小くする場合は、誤差のことを考えて ( 倍長計算にする等の ) 工夫が必要となる。