小清水 (@curekoshimizu) です。

久しぶりの投稿になります。

長期にわたり転職活動をしており、かなり投稿に時間が空いてしまいました。

今回の記事は

$\displaystyle \frac{1}{a}$ に関する内容です！

逆数の近似から精度を高めたい！！！

$\displaystyle \frac{1}{a}$ の近似値 ${x}_{1}$ が与えられたときに

精度を高めたいということありませんでしょうか？

このモチベーションは色んな場面で登場します。

それはなぜか？

これは、ハードウェアが逆数の近似を行う命令を持っていることが多く、それを利用したいためです。

例えばこちらのブログでは

qiita.com

AVX-512 の vrcp28pd 命令を使って $\displaystyle \frac{1}{a}$ の近似値を得た際のお話が書かれています。

AVX-512 サポートしていれば vrcp28pd命令・vrcp14命令などの逆数近似命令をサポートしています。

AVX-512 をサポートしているプロセッサはかなり限られていますが、 AVX をサポートしているだけでも rcpps という逆数近似命令を持っています。

また、もう今は死んでいると言っても過言ではない IA-64 にも frcpa 命令 (Floating-Point Reciprocal Approximation) という逆数近似命令があります。

ちなみに IA-64 は Hardware として普通の除算命令をもっておらず (IEEE-754準拠した除算の意)、この近似逆数命令からうまく精度を高めて IEEE-754 に適合した除算に変える必要がありました。

その点では IA-64 ユーザーからすれば、モチベーションである「逆数の近似から精度を高めたい！」ということは普通に行われていました。

ここでは代表的な逆数近似命令を4つ挙げてみました。あくまでも近似ではありますが、きちんと誤差保証もついています。

その対応表は例えば次です：

要求	命令	最大相対誤差
IA-64	frcpa	$2^{-8.886}$
AVX	rcpss	$1.5 \times 2^{-12}$
AVX-512	vrcp14pd	$2^{-14}$
AVX-512	vrcp28pd	$2^{-28}$

こうした話は CPU だけでなく、 NVIDIA GPU では fdivdef などで検索してみていただければと思います。

参考文献：

vrcp28pd・vrcp14pd・rcpss の精度について：https://www.intel.com/content/dam/www/public/us/en/documents/manuals/64-ia-32-architectures-software-developer-instruction-set-reference-manual-325383.pdf
frcpa の精度について： https://www.intel.co.jp/content/dam/www/public/ijkk/jp/ja/documents/developer/adag_j.pdf

これらはハードウェアの制約からきまった精度であり、

IEEE754-2008 が規定するような精度は不要なものの、もう少し精度を高めたい。そしてそれに保証をつけたい！

こういう時に役立つ記事を目指します。

精度を高める前に重要な事実

逆数の近似において

$e = 1 - a {x}_{1}$

という式は重要な意味をもちます。

$e = 1 - a {x}_{1}$ は $\displaystyle \frac{1}{a}$ と ${x}_{1}$ の(符号付き)相対誤差を表す。

(証明) $\displaystyle \frac{1}{a}$ と ${x}_{1}$ の(符号付き)相対誤差を記し、そのまま計算すると、

(符号つき)相対誤差 $\displaystyle = \frac{ \frac{1}{a}- {x}_{1}}{ \frac{1}{a} } = 1 - a {x}_{1} = e$

となることから証明完了。

□

また、この $e$ は FMA (Fused-Multiply-ADD) $(a \times b + c 型)$ の1命令で計算できる点も注目したい事実です。

FMA をご存知ない方はこちらの記事をご覧ください：

math-koshimizu.hatenablog.jp

精度を高めるアルゴリズム紹介

精度を高めるアルゴリズムはいくつか知られております。そのうちの3つを紹介します。

精度向上アルゴリズムその1 $e = 1 - a {x}_{1}$ $x_{2} = {x}_{1} + {x}_{1} e$ FMA 2命令で得られる $x_{2}$ は相対誤差 $e^2$ となる

(証明)

(符号つき)相対誤差 $\displaystyle = \frac{ \frac{1}{a}- {x}_{2}}{ \frac{1}{a} } = 1 - a ( {x}_{1} + {x}_{1} e ) = e - a{x}_{1}e = e^2$

□

精度向上アルゴリズムその2 $e = 1 - a {x}_{1}$ $x_{2} = {x}_{1} + {x}_{1} e$ $x_{3} = {x}_{1} + {x}_{2} e$ FMA 3命令で得られる $x_{3}$ は(符号つき)相対誤差 $e^3$ となる

(証明) 同様.

□

精度向上アルゴリズムその3 $e = 1 - a {x}_{1}$ $x_{2} = {x}_{1} + {x}_{1} e$ , $\quad E = e^{2}$ $x_{4} = {x}_{2} + {x}_{2} E$ FMA 3命令で得られる $x_{3}$ は相対誤差 $e^4$ となる。
ただし 2命令目は並列に実行すること。

(証明) 同様.

□

上の証明は証明になっているか？

一般的な書物には、上の内容だけが書かれているのですが、

証明としてこれでいいのでしょうか？

というのは、FMA であっても丸め誤差が生じるからです。

実数世界であれば上の証明で正しいのですが、

2・3命令であっても丸め誤差が生じる演算が数回発生した結果、

得られた結果は ${e}^{2}, {e}^{3}, {e}^{4}$ などの結果と比べて近い結果になるのでしょうか。

本ブログではここについて踏み込みたいと思っています。

丸めの影響も含めた証明

ここで丸め関数として RN (x の最近接偶数方向丸め) のみを取り上げることにします。

他の RU・RD・UZ についても同様に証明できますので割愛します。

また、丸めについては下の記事を参考ください：

math-koshimizu.hatenablog.jp

□

目的：「精度向上アルゴリズムその1」を RN も含めて誤差評価を行う： $e' = \mathrm{RN}(1 - a {x}_{1})$ $x'_{2} = \mathrm{RN}({x}_{1} + {x}_{1} e')$

これについて 2つの方針で考えていきたいと思います。

この評価は本などに書かれていなかったため、計算間違いをしている可能性もあり、間違っている場合にはご指摘いただけると助かります。

前者が大変、後者が簡単な式変形になります。

アルゴリズムその1 – 丸め誤差評価気合編

気合で三角不等式でなんとか結果をだした評価式がこちらになります。

わかりやすいように、次のように記号を定義します：

$e = (1 - a {x}_{1})$
$e' = \mathrm{RN}(e)$
$y = ({x}_{1} + {x}_{1} e')$
$y' = \mathrm{RN}(y) (= x'_{2})$

(方針) ： $y'$ の相対誤差を出す前に先に $y$ の評価を試みる。

$y$ と $\displaystyle \frac{1}{a}$ の相対誤差

$= | 1 - a y |$

$= | 1 - a ({x}_{1} + {x}_{1} e') |$

$= | e - a {x}_{1} e' |$

$\leq | e - e' | + | e' - a {x}_{1} e'|$

$=| e - e' | + | e' e |$

$= | e - e' | + | e (e' - e) + {e}^{2} |$

$\leq |e| \left| \frac{e - e'}{e} \right| + {|e|}^{2} \left| \frac{e - e'}{e} \right| + {|e|}^{2}$

丸め誤差 RN の性質より ( 以前の記事Proposition2より)

$\leq |e| {2}^{-p} + {|e|}^{2} {2}^{-p} + {|e|}^{2}$

と良い評価式を得る。

最後に

$y'$ と $\displaystyle \frac{1}{a}$ の相対誤差

$= | 1 - a y' |$

$\leq | 1 - ay | + | ay - a y' |$

$= | 1 - ay | + | a | | y - y' |$

$= | 1 - ay | + |a y| \left| \frac{y - y'}{y} \right|$

丸め誤差 RN の性質より ( 以前の記事Proposition2より)

$= | 1 - ay | + |a y | {2}^{-p}$

$= | 1 - ay | + |a ({x}_{1} + {x}_{1} e') | {2}^{-p}$

$= | 1 - ay | + | (1-e)(1+e') | {2}^{-p}$

$\leq | 1 - ay | + \left( | (1-e)(1+e') - (1-e)(1+e) | + |(1-e)(1+e)| \right) {2}^{-p}$

$= | 1 - ay | + \left( |1-e||e'-e| + |(1-{e}^{2}| \right) {2}^{-p}$

$\leq | 1 - ay | + \left( (1+|e|) |e| \left| \frac{e - e'}{e} \right| + 1 \right) {2}^{-p}$

丸め誤差 RN の性質より ( 以前の記事Proposition2より)

$\leq | 1 - ay | + \left( |e|{2}^{-p} + {|e|}^{2} {2}^{-p} + 1 \right) {2}^{-p}$

この結果に先ほど計算した $y$ の評価式を使って

$\leq |e| {2}^{-p} + {|e|}^{2} {2}^{-p} + {|e|}^{2} + \left( |e|{2}^{-p} + {|e|}^{2} {2}^{-p} + 1 \right) {2}^{-p}$

$= {|e|}^{2} (1 + {2}^{-p} + {2}^{-2p}) + |e| ( {2}^{-p} + {2}^{-2p} ) + {2}^{-p}$

という結果を得る。

この方針で得られた

丸め誤差を含めたアルゴリズムその1 の評価式は

${|e|}^{2} (1 + {2}^{-p} + {2}^{-2p}) + |e| ( {2}^{-p} + {2}^{-2p} ) + {2}^{-p}$

と ${|e|}^{2}$ に比べて少し大きくなっています。

2命令しかないのにこんなに評価が長くなってしまいました…。

同じような話がアルゴリズムその2, その3 にも適用できるとは思うのですが…。正直やりたくないです…。

そこで次のような式変形を思いついてみました！

アルゴリズムその1 – 丸め誤差評価数式処理編

ここで重要な性質である次の性質を利用してみたいと思います。

Proposition 1：相対誤差に関する性質
2進p桁精度環境とする。このとき、
(丸め可能な)実数 $x$ に対してある $\epsilon_x$ で $|\epsilon_x| \leq {2}^{-p}$ なるが存在し $\mathrm{RN}(x) = (1 + \epsilon_x)x$ と表せる.

(証明)

丸め誤差 RN の性質より ( 以前の記事Proposition2より)

$\displaystyle \left| \frac{ \mathrm{RN}(x) - x }{ x } \right| \leq {2}^{-p}$

ここで $x$ > 0 とすると

$\displaystyle -2^{-p} x \leq \mathrm{RN}(x) - x \leq 2^{-p} x$

なので

$\displaystyle (1 -2^{-p}) x \leq \mathrm{RN}(x) \leq (1 + 2^{-p}) x$

となる。 $x$ < 0 の場合も同様にすると

$\displaystyle (1 +2^{-p}) x \leq \mathrm{RN}(x) \leq (1 - 2^{-p}) x$

となる。すなわち $\epsilon_x$ で $|\epsilon_x| \leq {2}^{-p}$ とできる。

□

それでは評価に移ります。

再びわかりやすいように、変数を次のように定義します：

$e = (1 - a {x}_{1})$
$e' = \mathrm{RN}(e)$
$y = ({x}_{1} + {x}_{1} e')$
$y' = \mathrm{RN}(y) (= x'_{2})$

ここで上の 相対誤差に関する性質 を使い、

$\epsilon_y, \epsilon_e$ で $|\epsilon_y| \leq {2}^{-p}, |\epsilon_e| \leq {2}^{-p}$ となるもので

$y' = (1 + \epsilon_y) y, e' = (1+\epsilon_e)e$ とできる。

$y'$ と $\displaystyle \frac{1}{a}$ の相対誤差

$= | 1 - a y' |$

$= \left| 1 - a (1 + \epsilon_y) y \right|$

$= \left| 1 - a (1 + \epsilon_y) ({x}_{1} + {x}_{1} e')\right|$

$= \left| 1 - a (1 + \epsilon_y) \left({x}_{1} + {x}_{1} (1 + \epsilon_e) e \right) \right|$

これをすべて展開して $e = (1 - a {x}_{1})$ を使って $a {x}_{a}$ を削除すると

$= \left| {e}^{2} (1 + \epsilon_e + \epsilon_y + \epsilon_e \epsilon_y) + e (- \epsilon_e -\epsilon_e \epsilon_y ) - \epsilon_y \right|$

ここで $|\epsilon_y| \leq {2}^{-p}, |\epsilon_e| \leq {2}^{-p}$ を利用して

$\leq {|e|}^{2} (1 + {2}^{-p} + {2}^{-p} + {2}^{-2p} ) + |e| ( {2}^{-p} + {2}^{-2p} ) + {2}^{-p}$

$= {|e|}^{2} (1 + {2}^{-p+1} + {2}^{-2p} ) + |e| ( {2}^{-p} + {2}^{-2p} ) + {2}^{-p}$

となる。

よってこの方針で得られた

丸め誤差を含めたアルゴリズムその1 の評価式は

${|e|}^{2} (1 + {2}^{-p+1} + {2}^{-2p} ) + |e| ( {2}^{-p} + {2}^{-2p} ) + {2}^{-p}$

となりました。先ほどの評価より、式変形は簡単になりましたが、少しだけ結果が悪くなりました。

アルゴリズムその1 丸め誤差を含む評価まとめ

目的：「精度向上アルゴリズムその1」を RN も含めて誤差評価を行う： $e' = \mathrm{RN}(1 - a {x}_{1})$ $x'_{2} = \mathrm{RN}({x}_{1} + {x}_{1} e')$

ここで得られる $x'_{2}$ は $\frac{1}{a}$ に対する相対誤差は

${|e|}^{2} (1 + {2}^{-p} + {2}^{-2p}) + |e| ( {2}^{-p} + {2}^{-2p} ) + {2}^{-p}$
${|e|}^{2} (1 + {2}^{-p+1} + {2}^{-2p} ) + |e| ( {2}^{-p} + {2}^{-2p} ) + {2}^{-p}$