1

64 ビット マシン上の Matlab で double として表される一連の無理数Nのサイズの厳密な下限は何ですか? たとえば、さまざまなランダムな pi のチャンクをエンコードする ~10^12 の double を乗算した後、どのような精度を期待できますか?

4

3 に答える 3

3

タイト バウンドを要求した場合、@EricPostpischil の応答は、すべての操作が IEEE 754 倍精度で実行される場合の絶対エラー バウンドです。

自信を持って言えば、誤差の統計的分布と理解しています。[-e/2,e/2] の誤差の分布が均一であると仮定すると、数学スタック交換で M 回の操作を行った後の理論上の誤差の分布を求めることができます...タイトな境界は、どういうわけか非常に保守的だと思います。

いくつかの Smalltalk コードを使用して、これらの統計の実験的推定を説明しましょう (大きな整数/小数演算を行う言語であればどれでも実行できます)。

nOp := 500.
relativeErrorBound := ((1 + (Float epsilon asFraction / 2)) raisedTo: nOp * 2 - 1) - 1.0.
nToss := 1000.
stats := (1 to: nToss)
    collect: [:void |
        | fractions exactProduct floatProduct relativeError  |
        fractions := (1 to: nOp) collect: [:e | 10000 atRandom / 3137].
        exactProduct := fractions inject: 1 into: [:prod :element | prod * element].
        floatProduct := fractions inject: 1.0 into: [:prod :element | prod * element].
        relativeError := (floatProduct asFraction - exactProduct) / exactProduct.
        relativeError].
s1 := stats detectSum: [:each | each].
s2 := stats detectSum: [:each | each squared].
maxEncounteredError  := (stats detectMax: [:each | each abs]) abs asFloat.
estimatedMean := (s1 /nToss) asFloat.
estimatedStd := (s2 / (nToss-1) - (s1/nToss) squared) sqrt.

nOp=20 double の乗算で次の結果が得られます。

relativeErrorBound -> 4.440892098500626e-15
maxEncounteredError -> 1.250926201710214e-15
estimatedMean -> -1.0984634797115124e-18
estimatedStd -> 2.9607828266493842e-16

nOp=100 の場合:

relativeErrorBound -> 2.220446049250313e-14
maxEncounteredError -> 2.1454964094158273e-15
estimatedMean -> -1.8768492273800676e-17
estimatedStd -> 6.529482793500846e-16

nOp=500 の場合:

relativeErrorBound -> 1.1102230246251565e-13
maxEncounteredError -> 4.550696454362764e-15
estimatedMean -> 9.51007740905571e-17
estimatedStd -> 1.4766176010100097e-15

標準偏差の増加は、誤差範囲の増加よりもはるかに遅いことがわかります。

更新:最初の近似(1+e)^m = 1+m*e+O((m*e)^2)で、分布は m*e が十分に小さい限り、[-e,e] で均一な m のほぼ合計であり、この合計は分散の正規分布 (ガウス) に非常に近いm*(2e)^2/12です。Matlabstd(sum(rand(100,5000)))で近くにあることを確認できます。sqrt(100/12)

についてはまだ真であると考えることができます。m=2*10^12-1つまり、およそm=2^41ですm*e=2^-12。この場合、グローバル誤差は準正規分布であり、グローバル誤差の標準偏差はsigma=(2^-52*sqrt(2^41/12))またはおおよそsigma=10^-10

P(abs(error)>k*sigma) を計算するには、 http://en.wikipedia.org/wiki/Normal_distributionを参照してください。

  • ケースの 68% (1 シグマ) では、10 桁以上の精度が得られます。
  • erfc(10/sqrt(2))6*10^22 のうち約 1 桁の精度を持つ確率が 9 桁未満になる確率が得られるので、4 桁の精度しか持たない確率を計算させます (倍精度で評価することはできません。アンダーフロー) !!!

私の実験的な標準偏差は理論的なものより少し小さかった (2e-15 9e-16 4e-16 for 20 100 & 500 double) が、これは入力エラー i/3137 i=1..10000 の偏った分布によるものに違いない...

これは、結果が入力の誤差の分布によって支配されることを思い出させる良い方法です。M_PI*num/den のような浮動小数点演算の結果である場合、誤差は e を超える可能性があります。

また、Eric が言ったように、* だけを使用するのは非常に理想的なケースであり、+ を混在させると事態がより早く悪化する可能性があります。

最後の注意: 最大誤差範囲に達する入力のリストを作成し、すべての要素を (1+e) に設定して 1.0 に丸め、理論上の最大誤差範囲を取得しますが、入力分布はかなり偏っています。 ! HEM WRONGすべての乗算は正確であるため、(1+e)^(2n-1) ではなく (1+e)^n しか得られないため、誤差は約半分にすぎません...

更新 2: 逆の問題

逆数が必要なので、一定レベルの信頼度 10^-c で k 桁の精度が得られるようなシーケンスの長さ n はいくらですか?

上記の近似では が必要なためk>=8、についてのみ回答します。(m*e) << 1

を取りましょうc=7。10^-7 の信頼度で k 桁を取得すると、5.3*sigma < 10^-k.
sigma = 2*e*sqrt((2*n-1)/12)それはn=0.5+1.5*(sigma/e)^2となりe=2^-53ます。
したがって、n ~ 3*2^105*sigma^2、sigma^2 < 10^-2k/5.3^2 として、次のように記述できます。n < 3*2^105*10^-(2*k)/5.3^2

AN k=9 桁未満の確率は、長さ n=4.3e12 の場合は 10^-7 未満、10 桁の場合は n=4.3e10 付近です。

15 桁で n=4 の数に到達しますが、ここでの正規分布の仮説は非常に大雑把であり、特に 5 シグマでの分布の裾が成り立たないため、注意して使用してください (Berry-Esseen の定理は、そのような分布が正規からどれだけ離れているかを制限しますhttp ://en.wikipedia.org/wiki/Berry-Esseen_theorem )

于 2012-08-21T09:54:11.537 に答える
2

すべての入力値、中間値、および最終値がアンダーフローまたはオーバーフローしないと仮定すると、前述のM演算の相対誤差は最大で (1+2 -53 ) M -1 です。

実数a0を倍精度に変換することを検討してください。結果はある数a0 •(1+ e ) で、ここで -2 -53 ≤ e ≤ 2 -53です (倍精度への変換は常に表現可能な最も近い値を生成する必要があり、倍精度値の量子は 2 -53であるため)最も近い値は常にクォンタムの半分以内です)。さらに分析するために、 eの最悪の場合の値2 -53を検討します。

1 つの (以前に変換された) 値を別の値で乗算すると、数学的に正確な結果はa0 •(1+ e ) • a1 •(1+ e ) になります。計算結果には別の丸め誤差があるため、計算結果はa0 •(1+ e ) • a1 •(1+ e ) • (1+ e ) = a0 • a1 • (1+ e ) 3となります。明らかに、これは (1+ e ) 3の相対誤差です。エラーが単純に (1+ e ) Mとして累積されることがわかります。これらの操作の場合: 各操作は、前のすべての誤差項を 1+ eで乗算します。

N個の入力がある場合、N回の変換とN -1 回の乗算が行われるため、最悪のエラーは (1+ e ) 2 N - 1になります。

この誤差の等式は、 N ≤1の場合にのみ達成されます。それ以外の場合、エラーはこの境界より小さくなければなりません。

この単純なエラー バウンドは、同種操作を伴うこのような単純な問題でのみ発生する可能性があることに注意してください。加算、減算、乗算、およびその他の演算が混在する典型的な浮動小数点演算では、境界を単純に計算することは通常不可能です。

N =10 12 ( M =2•10 12 -1) の場合、上限は 2 -53 の 2.000222062•10 12単位未満であり、.0002220693 未満です。したがって、計算結果は 10 進数 4 桁以下で十分です。(ただし、オーバーフローとアンダーフローを避ける必要があることに注意してください。)

(上記の計算の厳密さに注意してください: Maple を使用して、二項式 (1+2 -53 ) 2•10 12 -1の 1000 項を正確に計算し (最初の 1 項を削除して)、証明可能な値を追加しました。残りのすべての項の合計よりも大きい. 次に、Maple にその正確な結果を 1000 桁の 1000 桁に評価させたところ、上記で報告した範囲よりも小さかった.)

于 2012-08-21T03:28:43.083 に答える
1

64ビット浮動小数点数の場合、標準IEEE 754を想定すると、52+1ビットの仮数があります。

つまり、相対精度は1.0000...0から1.0000...1の間であり、小数点以下の2進数は52です（1.000 ... 0は、2進数で格納されているものと考えることができます。 mantissa AKA仮数）。

誤差は、52の1/2の2で割った値（解像度の半分）です。最悪の場合であるため、可能な限り1.0に近い相対精度を選択することに注意してください（そうでない場合は、1.111..11と1.111..01の間で、より正確になります）。

10進数では、doubleの最悪の場合の相対精度は1.11E-16です。

この精度でNdoubleを乗算すると、新しい相対精度（中間の丸めによる追加のエラーがないと仮定）は次のようになります。

1 - (1 - 1.11E-16)^N

したがって、pi（または任意のdouble 10 ^ 12）を乗算すると、エラーの上限は次のようになります。

1.1102e-004

つまり、約4〜5桁の自信を持つことができます。

CPUが中間結果の拡張精度浮動小数点数をサポートしている場合は、中間丸め誤差を無視できます。

拡張精度FPU（浮動小数点ユニット）が使用されていない場合、中間ステップでの丸めにより、追加のエラーが発生します（乗算によるものと同じ）。つまり、厳密な下限は次のように計算されます。

1 -
((1 - 1.11E-16) * (1 - 1.11E-16) * (1 - 1.11E-16)
                * (1 - 1.11E-16) * (1 - 1.11E-16) % for multiplication, then rounding

               ... (another N-4 lines here) ...

                * (1 - 1.11E-16) * (1 - 1.11E-16))

= 1-(1-1.11E-16)^(N*2-1)

Nが大きすぎると、実行に時間がかかりすぎます。発生する可能性のあるエラー（中間丸めあり）は2.2204e-012であり、中間丸めなし1-(1 - 1.11E-16)^N=1.1102e-012と比較して2倍になります。

概算すると、中間の丸めによってエラーが2倍になると言えます。

piを10^12倍に乗算し、拡張精度FPUがなかった場合。これは、続行する前に（FPU結果が蓄積されないようにコンパイラが命令を並べ替えていないことを確認してください）、メモリに中間ステップを書き込む（そしておそらく何か他のことをする）ためである可能性があります。エラーは次のとおりです。

2.22e-004

小数点以下の桁数に自信があるからといって、正確に小数点以下の桁数になることはありません。

たとえば、答えが次の場合：

1.999999999999、エラーが1E-5の場合、実際の回答は2.000001234になる可能性があります。

この場合、小数点以下1桁も間違っていました。しかし、それは本当にあなたがどれほど幸運であるかに依存します（答えがこのような境界にあるかどうか）。

このソリューションは、double（答えを含む）がすべて正規化されていることを前提としています。非正規化された結果の場合、明らかに、非正規化される2進数の桁数は、その数桁の精度を低下させます。

于 2012-08-21T02:07:00.620 に答える