数値計算

Microsoft Cloud Numerics の数学関数をテストする

Stuart Brorson
Alan Edelman
Ben Moskowitz

 

数学の計算を実行する必要があるとします。たとえば、角度 34 度の正弦値を求める必要があるとします。どうするでしょうか。おそらく、電卓やコンピューターなどのスマート デバイスを使用するでしょう。デバイスに「sin(34)」と入力すれば、多くの場合、小数点以下 16 桁の精度で答えが得られます。では、その答えが正しいかどうか、どうすればわかるでしょう。

電子機器を使って計算結果を得ることに慣れているあまりに、だれもその答えが正しいかどうか疑わなくなっています。ほとんどだれもが、当然のように機械が正解を教えてくれると思っています。しかし、ソフトウェアの品質管理技術者という、正確さを当たり前とは考えない人々もいます。彼らの仕事は必ず正解を得られるようにすることです。今回は、新しい Microsoft Cloud Numerics 数学ライブラリの数学関数をテストする方法について説明します。

ほとんどの科学計算や工学計算には、浮動小数点演算を使用します。IEEE は 1985 年に浮動小数点演算の基本動作を標準化しました。この標準の大きな特徴は、すべての数をコンピューターで表せるわけではないことを認めているところです。実数のセットは無限で数えられませんが、IEEE の浮動小数点数のセットは有限で数えられます。これは、浮動小数点数の数値表現が、固定数のビットを使用しているためです。IEEE の浮動小数点数は構造上有理数で、π などよく使用される数値が正確に表されず、そうした数値を使用する sin(x) などの演算も正確に表されないことを意味します。

さらに、整数とは違い、IEEE の浮動小数点数の間隔は局所的には均一ですが、数値自体が大きくなると対数的に増加します。この対数的側面を図 1 に示します。この図では、実数を表す直線上で IEEE の浮動小数点数の均一な固まりの場所を図式的に示しています。図では、(実数を表す直線上に埋め込まれている) 浮動小数点数の有効な固まりを垂直線で示しています。有効な浮動小数点数の間隔は x が大きくなるにつれて対数的に増加していることに注意します。2 つの隣接する浮動小数点数の間隔は、多くの場合最小精度単位 (ULP: Unit of Least Precision) または最終桁単位 （ULP: Unit in Last Place) と呼ばれ、多くの一般的な数学プログラムでは eps(x) という組み込み関数として提供されます。このような間隔の違いから、0 に近い小さな数値を 1 などの比較的大きな数値に加算しても結果に影響しません。

図 1 図式化した浮動小数点数のグリッド

IEEE の標準は、浮動小数点数の形式の体系化以外に、最も基本的な数学演算の戻り値の正確性についての厳密なガイドを提供します。たとえば、IEEE の標準は 4 つの算術演算 (+、-、*、および /) の戻り値は、"最善の概数" でなければならないと指定しています。つまり、返される答えは "数学的に正しい" 結果に "最も近い" 値になる必要があります。この要件には最初はいくらか抵抗があったものの、現在ではハードウェアで最善の概数が返されるようになっており、非常に一般的になっています。さらに興味深いことに、IEEE の浮動小数点数の仕様は、平方根関数が最善の概数を返すことも求めています。これが興味深いのは、次の 2 つの疑問を生じるためです。1 つは「IEEE の浮動小数点数表記を使用して無理関数をどの程度の精度で計算できるのか」、もう 1 つは「関数実装の精度をどのようにテストすべきなのか」という疑問です。

関数で精度の高い計算を行う

どの要素が計算の精度を決めるのでしょうか。実際には、アルゴリズムのあらゆる段階で数値の丸めが起こり得るため、不正確さが生じることは予想すべきです。こうした誤差は、相殺されることもありますが、相乗効果をもたらすこともあります。実践では、アルゴリズムの設計者が、数値計算の本質的な不正確さを考慮して設計することを学ばなければなりません。

計算誤差を引き起こす要因は、次の 2 つのカテゴリに分けられます。

1. 計算の実行に使用するアルゴリズムに関連する要因。特に、アルゴリズム自体の安定性と、数値の丸めによる誤差に対するアルゴリズムの感度を意味します。脆弱で不安定なアルゴリズムは精度の低い結果を返す傾向があり、堅牢で安定したアルゴリズムはより精度の高い結果を返します。
2. 関数自体とその入力の定義域に内在する性質。計算で有限 (無限の反対として) のビット数を使用すると、すべての関数実装で実現できる精度には理論上の限界があります。これは、計算における数値の丸めが誤差の原因となり、その出力が次の計算の入力に使用されると、この誤差が増幅するか減衰するかは関数自体の動作によって決まるためです。たとえば、特異点、ゼロ、またはすばやい振幅の近辺での関数値の計算は、同じ入力定義域でゆっくり変化する関数の計算と比べて不正確になります。こうした本質的に実現可能な精度は、関数がどの程度 "適切に条件付けられているか" によって決まります。

したがって、関数実装の精度をテストすることには、アルゴリズムが理論上可能な限り精度の高い結果を返すことを検証する作業が関係します。理論上の可能性の限界は、それぞれの入力における関数の条件付けによって決まります。別の言い方をすると、浮動小数点数を使用して結果を計算すると、誤差の 2 つの要因が持ち込まれる可能性があります。適切で、安定したアルゴリズムを使用すれば回避できる誤差 (要因 1) と、入力値に応じた関数の動作に関連するため回避が困難な誤差 (要因 2) の 2 つです。

数値分析での条件付けの概念は、いわゆる "条件数" によって定量化され、入力の変異に対する関数の感度を測定します。検討中の関数の詳細な性質に応じて (入力数、スカラー対行列など)、条件数には多数の複雑な式があります。最も単純なケースは、1/x、x2、sin(x)、および高校の代数で習ったその他の関数など、変数 1 つの微分可能関数です。このように単純なケースでは、関数 f の条件数は次の式で計算します。

式 1

これは、後で式 1 として参照します。通例どおり、f'(x) は、関数 f(x) の導関数を意味します。大きな条件数 (Kf >> 1) は、相対的な摂動に対する高い感度を示し、小さな条件数 (Kf <= 1) は、関数が相対的な摂動に対して感度が低いことを意味します。関数が特異点を持つかどうかは、直感的にわかります。つまり、1/x が x = 0 に近づくなど、関数が発散すると、導関数はこの影響を受けて大きな条件数を返します。

たとえば、関数 f(x) = 1/(1-x) を考えてみます。この関数は、x = 1 で明らかに発散します (つまり、無限になります)。この導関数は、f'(x) = 1/(1-x)2 です。これを式 1 に組み込み、文字を消去すると、この関数の条件数が求められます。

条件数 Kf(x) はほぼ x = 1 のときに無限になります。つまり、この関数を含む計算は、x が 1 に近づいたときに数値の丸めによる誤差などの摂動を感知する関数に関係します。これは、f(x) が x → 1 のときに発散し、この動作が予想されるためです。また、x が 0 に近づくと、条件数 Kf(x) は 0 に向かいます。これは、この関数を使用した計算が x → 0 のときに摂動に対する感度が低くなることを意味します。x が 1 よりもはるかに小さいとき、f(x) は定数値 1 になる傾向があるため、直感的に理解できます。

関数の精度をテストする

では、数学関数の実用的なテストはどのようになるでしょう。x が入力で y が返される出力値となる、スカラー関数 y = f(x) の実装をテストするとします。x を浮動小数点数の入力値とすると、テストには次の 3 つの項目が必要です。

1. テストする関数で計算された値 (ycomputed = fcomputed (x))。
2. "数学的に正しい" 結果。正確で、数学的に正しい答えがあらかじめ求められているとします。この答えを最も近い浮動小数値に丸めます (それにより、コンピューターで表せるようになります)。この値を ytrue = ftrue (x) とします。
3. テストの許容誤差。関数によっては、この許容誤差が 0 の場合があります。これは、ycomputed が数学的に正しい ytrue とまったく同じ値であることを意味します。たとえば、浮動小数点の仕様では、関数 sqrt() は正確な答えに最も近い浮動小数点数値 (つまり、最善の概数) を返すことを求めています。ただし、一般的な場合は、すべての関数の戻り値が正確な値に対する最善の概数になると期待することはできません。そのため、テストの許容誤差には、f(x) の戻り値で誤差がどれだけ許されるかについての情報を含んでいる必要があります。許容誤差は、関数 f(x) の詳細と入力の正確な値である x に依存することがあります。

これらの項目を使用して、関数実装の精度は、次の場合に受け入れ可能と見なされます (つまり、テストに合格します)。

| ycomputed - ytrue | < tol(f; x),

ここでは、許容誤差は入力値 x と関数 f 自体の動作の両方に依存します。

言葉で説明すると、この式は、返される値と "正しい" 値の違いがテストの許容誤差の範囲内であれば、関数はテストに合格することを表しています。| ycomputed - ytrue | は、計算関数の "絶対" 誤差です。

このような定式化において、関数のテストは関数の入力定義域内に収まる多数の入力値を作成することで実行され、テスト対象の関数でこれらの値から計算を行い、関数出力の ycomputed を最善の概数 (数学的に正しい) 値の ytrue と比較します。入力値は、関数の有効な入力定義域のすべての関連する部分をカバーするように選択します。

ここで、テスト担当者は、関数の "正しい" 結果とは何で、妥当な許容誤差とは何なのか、という 2 つの疑問を持ちます。

最初の問いに答えるには、"無限精度" の数学パッケージを使用して、テストに使用する入力と出力の組み合わせを作成するのが最も簡単です。このようなパッケージは、値を計算するのに 32 ビットまたは 64 ビットの浮動小数点数を使用しません。代わりに、計算速度は遅くなりますが、数値表現か、もっと良いのはシンボル表現を使用して、任意の多くの桁数で計算する機能を提供します。無限精度の数学計算を実装する商用の数学パッケージがいくつかあります。また、多くの無限精度の数学ライブラリを一般的な言語に組み込むことができます。近代的な速度で、無限精度の数学計算を使用する機能は最新の発明で、この種のテストを使いやすくしています。これらのリソースはすべて、関数の浮動小数点実装をテストするのに役立つ、いわゆる "黄金値" の組み合わせを作成するのに十分です。

これらをまとめ、テストの許容誤差を取得する

黄金値を作成したら、「妥当なテストの許容誤差とは何か」というもう 1 つの疑問に答える必要があります。正しいテストの許容誤差を取得することは、テストの重要な構成要素です。許容誤差が非現実的なほど小さいと、最善のアルゴリズムで計算しても、テストには合格しません。一方、テストの許容誤差が大きすぎると、アルゴリズムに問題があっても、テストに合格してしまいます。

関数テストの許容誤差を判断する鍵となるのが、前述の条件数です。条件数の式は、次のように再配置できます。

式 2

この式は、後で式 2 として参照します。∆x をある浮動小数点数と次の浮動小数点数の間隔とすると、この式はある浮動小数点数 x からその隣の浮動小数点数 x + ∆x に移動する際に f(x) の出力がどれだけ変化するかを示します。従来、コンピューターの数値フィールドは ∆x に関数 eps(x) を使用し、隣接する 2 つの浮動小数点数の間隔を決めます。グリッドの間隔は一定ではないため (図 1 参照)、eps(x) は x の関数になります (既に述べたように、ある浮動小数点数と最も近い浮動小数点数の間隔は ULP、つまり最下位ビットで表される大きさにも関連しています)。

次に、テストの出力誤差 (ycomputed - ytrue) がこの変化の数倍よりも小さくなることを求めます。つまり、C を定数として、次のように求めます。

直感的に、この式から次のような考えが浮かびます。x が浮動小数点グリッドを 1 つ進む場合、出力の変化は式 2 以下になるべきです。精度関数 f(x) は、条件数から派生する量だけ出力を変化させ、それ以上は変化させません。そのため、f(x) の計算の間に発生する出力の摂動は、浮動小数点グリッドを 1 つ進んだときに発生する変化よりも小さくなります。定数 C は "仮因数" です。明らかに、C が増加するにつれて、許容誤差が大きくなります。そのため、この定数は理論上関数の最小許容誤差と解釈できます。実際のテストでは、許容誤差と解釈する C は 1 から 10 までの整数値を取り、ULP で表されます。

最後に、条件数の定義を使用して、テストの許容誤差を表すように式を再配置できます。

| ycomputed - ytrue | ≤ C | f'(x) | eps(x) = tolerance

これは、目的のテストで、あらゆる入力 x の関数 f(x) を正確に実装することで満たす必要がある式です。言い換えると、これはスカラー数学関数を検証するのに使用する式です。この式は、すべての有効な入力値 x に対して成り立つ必要があることに注意します。

答える必要がある最後の疑問は、「テストで使用する C の正しい値は何か」です。C は許容 ULP 誤差の数として解釈されるため、C は次数 1 の整数です。最初の C の値 (通常は 10) から開始してテストを実行します。テストに合格したら、C を 1 減らして、再度テストを実行します。テストに不合格になるまで、この手順を繰り返し C の値を減らします。次に、最後にテストに合格する C の値を選択します (ただし、このループには必ず手動で行う部分を加え、妥当性を保証します)。このプロセスの目的は、可能な限り多くテストを行い (関数がテストに合格できるようにしつつ)、関数ができるだけ厳密にテストされたことを保証することです。適切に動作する関数は、多くの場合 1 から 3 までの C 値を要求しますが、関数が複雑になるとより大きな C 値でテストに合格することを求められる場合もあります。合格するのに C > 10 を求める関数実装は非常に少なく、C > 10 を求める関数は多くの場合最適な実装ではないことがわかります。

ソフトウェアのテストの観点では、関数テスト用の C が決まった (テストに合格した) 後に、次のテストで不合格になる場合は、関数実装の何かが、おそらく悪い方へ変更されたことがわかります (回帰が発生したなど)。当然、アルゴリズムの設計にはトレードオフがあります (速度と精度のトレードオフなど)。つまり、回帰が小さければ、C に余裕を与えることに価値がある場合があります。同様に、アルゴリズムを改良した後、C をさらに狭められないか検討することにも価値があります。いずれにしても、ソフトウェアのテストの 1 つの仕事は、ソフトウェアの残りの部分を開発する際に、テストの許容誤差を管理することになります。

黄金値のテストとその先

ここまでで説明したテスト手法では、テスト対象の関数の入力定義域全体に分散し、事前に計算した入力と出力の組み合わせを使用しています。この種のテストを "黄金値のテスト" と呼びます。無限精度のテストは、非常に精度が高いことがわかっている入力と出力の黄金の組み合わせを提供することを意味します。組み合わせはファイルに格納し、実行時にテスト プログラムに読み取ります。関数の精度をテストするブラックボックス手法として、この方法は非常に適切に機能します。ただし、ほかにも関数の動作を検証する重要な方法を提供する関数テストもあります。たとえば、次のようなテストがあります。

• 特殊値のテスト: 多くの関数は、cos(0) => 1、sin(π) => 0、gamma(1/2) = √π など、特定の入力に対して、既知の理論上正確な値を返します。特殊値のテストでは、既知の固定入力をテスト対象の関数に渡し、関数の戻り値を既知の結果と比較します。当然、π などの無理数は浮動小数点グリッドにはそのまま存在しないため、それらの無理数に最も近い浮動小数点の近似値をテストで使用する必要があり、ゼロ以外の正弦値とテストの許容誤差 (以前に計算しました) を、計算結果の検証に使用します。
• 恒等性のテスト: 多くの関数には、既知の定義域のすべての入力に成り立つ恒等式があります。たとえば、すべての入力 x に対して sin(x)^2 + cos(x)^2 == 1 が成り立ちます。また、arcsin(x) == -i log(i*x + sqrt(1 – x2)) も成り立ちます。恒等性のテストと特殊値のテストの違いは、恒等性のテストはすべての入力に対して成り立ちますが、特殊な値テストは 1 つの特定の入力値にしか成り立たないことです。恒等性のテストは関数どうしの関係を検証します。

恒等性のテストには注意が必要です。1 つには、一部の恒等式は不適切に条件付けされていることです。たとえば、最終結果のサイズは適切に調整されるとしても、中間結果は非常に大きくなります。この場合、計算の中間手順で発生する小さな許容誤差により、誤ってテストが不合格になる場合があります。また、多くの関数 (たとえば、逆三角関数) は、複素平面で複数の値を持ちます。むやみに恒等性のテストを使用すると、恒等式の左辺と右辺が複素平面の異なるリーマン シート上に正しい値を返したとき、誤ってテストが不合格になる場合があります。テスト担当者は、この両方の問題を回避するために、注意深く恒等性のテストを作成する必要があります。

• 逆関数のテスト: このテストでは、逆関数の計算とその結果がテストの入力と同じかどうかを検証します。たとえば、正の数 x に、log(exp(x)) = x が成り立つかどうかをテストできます。これは、恒等性のテストの特殊な例といえます。ただし、非常に多くの関数に逆関数があり、数学的な一貫性のためには逆関数が本来の入力を返すことが求められるため、関数とその逆関数が一貫した方法で動作することを検証する、別のテスト カテゴリを使用します。恒等性のテストに当てはまるのと同じ注意 (条件付けや異なるリーマン シート上に値を返すことなど) が逆関数テストにも当てはまります。
• 級数総和テスト: ほぼすべての超越関数が、いくつかの入力定義域をまたがって級数の拡大を許容します。級数総和テストでは、テスト対象の関数の戻り値と、テストで明示的に総和した級数の値を比較します。級数総和を恒等性のテストの一種と見なすことができますが、適切な級数総和テストを作成する際の考慮事項は、複数の関数で似ています (たとえば、安定性と収束の基準)。そのため、この種のテストは恒等性のテストとは概念的に異なると考えます。
• 連続分数展開: いくつかの関数は連続分数展開を使用して既知の定義域で評価できることがあります。そうした展開が存在する場合、級数総和よりも非常にすばやく収束する場合があります。そのため、特定の関数には、級数総和テストの代わりに連続分数展開を使用することができます。
• 構成値テスト (別名モデルベース テスト): 構成値テストでは、テスト自体においてシンプルで検証可能な正しい関数の実装を作成します。この種のテストは、必ずしもスカラー数学関数には当てはまりません。ただし、優れた数学パッケージは分析関数のセットには制限されません。たとえば、床関数、天井関数、MOD 関数、ベクトル用の連結演算子などを検討します。これらのそれぞれの関数は、基本的なコンピューター言語の構造体を使用して関数を実装するコードを作成することでテストできます。たとえば、床関数、天井関数、およびさまざまな数学のチョップ関数は、(検討中の関数の正確な動作を模倣するいくつかの簡単な数学を実行して) 浮動小数点値を整数に変換し、結果を浮動小数に戻すキャスト演算子を使用してモデル化できます。この種のテストのメリットは、モデル実装が多くの場合簡単で、検査によって検証可能だということです。
• エラー テスト: このテストは、プログラマが明らかにすばやく問題を把握できるわかりやすいエラーを返すことで、無効な入力が適切に処理されるようにします。たとえば、実数で作業する際、sqrt(x) は x < 0 の場合、役に立つエラー (NaN) を返します。一部の入力エラーは、適切に処理しないとコンピューターをクラッシュさせることがわかっており、これが起こらないようにするのはテスト担当者の役割です。

この種のテストのメリットは、数学自体だけを使用して関数の正確さをクロスチェックできることです。これらの手法は、数学の自己一貫性をテストの基盤として使用するため、"適切" と考えられており、サードパーティ製ソフトウェアを使用した数学計算に依存しません。

方法を検証する

条件数の許容誤差のメリットを活かす方法の 1 つは、入力定義域に特異点を持つスカラー数学関数を検証するのに使用した許容誤差を調べる方法です。たとえば、よく知られた三角関数の正接と逆正接は、実数値の線に沿った無限数の点における特異点です。固定されると、定数の許容誤差が精度の検証に使用されることになり、これらの関数が入力値が特異点に近いかどうかのテストで不合格になるか、テストの許容誤差が必然的に非常に大きくなりテストが効果的でなくなります。条件数の許容誤差により、関数の入力定義域全体が少数の ULP に収まるようにスカラー関数を検証できます。図 2 に示すのは、Microsoft Cloud Numerics 製品で配布された関数の認定の中で見つけた、合格する ULP 値の例です。つまり、対応する関数がテストに合格する ULP 値です。

図 2 合格する ULP 値

関数	入力定義域	合格する ULP	コメント
Atan2	実数	1	2 入力、4 象限逆正接
Cot	実数	3
Coth	複素数	3
Cuberoot	実数	2
Ellipe	0 ≤ x ≤ 1 となる実数 x	4	楕円積分 E
Ellipk	0 ≤ x ≤ 1 となる実数 x	6	楕円積分 K
Expm1	実数	3	Exp(x) - 1
Gamma	実数	3
Log	実数	3
Log	複素数	3
Log1p	実数	1	Log(1+x)
Log1p	複素数	1	Log(1+x)
Psi	実数	6	Psi 関数 — ガンマ関数の導関数
Sqrt	実数	0	IEEE の仕様では、"最善の概数" を求めます。
Tan	実数	1
Tanh	実数	4
Tanh	複素数	4

図 2 の例から導き出せる結論は、少数の ULP の範囲内で正確なアルゴリズムを使用して特別な関数が実装されることです。図 2 で示される最悪の場合は、誤差が 6 ULP 以下で、これは数値の最低数位の小数点以下 log10(2^6) = 1.8 桁と一致します。倍精度浮動小数では、これは 1.3323e-015 の相対誤差に相当します。

最先端のテスト

ここまでで、厳密なテスト手法を Microsoft Cloud Numerics 数学パッケージで使用されるアルゴリズムの検証にうまく適用できました。関数の条件数から導き出したテストの許容誤差を使用することで、誤差が関数の戻り値の最後の数桁でのみ発生している場合でも、不正確な関数の実装を特定して修正できるようになります。テストの許容誤差を固定すると、故意の誤り (関数の特異点の近接など) を引き起こさない厳密なテストを提供できません。テストの大部分は、無限精度の数学パッケージを使用して事前に計算され、テスト時に読み取るファイルに格納された黄金値の入力値と出力値を使用します。黄金値の使用に加え、さまざまな数学の恒等式を使用して関数の実装を照合し、関数の動作を検証します。

Microsoft Cloud Numerics のユーザーは、たくさんの数学と統計のライブラリ関数を使用でき、関数は最新のソフトウェア テスト プラクティスを使用して徹底的に検証されていることを確信できると信じています。

詳細については、David Goldberg の浮動小数点数学の古典的な説明である「What Every Computer Scientist Should Know About Floating-Point Arithmetic」(すべてのコンピューター科学者が浮動小数点演算について知っておくべきこと) を一読することをお勧めします (bit.ly/vBhP9m (英語) などの Web サイトで入手可能)。

Stuart Brorson は、マサチューセッツ州ケンブリッジのマイクロソフト NERD センターの SDET です。彼がマイクロソフトに参加したのは、ユーザーが並列スーパーコンピューターで数学計算および行列計算を実行できるソフトウェアを開発した、対話型のスーパーコンピューティングをマイクロソフトが取得したときです。趣味として、Brorson は電子機器のハッキングと、ボストン界隈でのアイリッシュ フィドルの演奏を楽しんでいます。

Ben Moskowitz は、マイクロソフトで数学ベースのプロジェクトのテストを専門とする SDET です。彼は、Cloud Numerics と、2011 年に Visual Studio Lighthouse Award を受賞した最適化パッケージとプラットフォームである Solver Foundation のリリースに寄与しました。プライベートでは、妻と共に過ごし、やんちゃ坊主たちの世話をしています。

Alan Edelman は、数学の教授でマサチューセッツ州ケンブリッジにある MIT のコンピューター サイエンスおよび AI 研究室のメンバーです。Edelman 教授は、数値分析、並列コンピューティング、およびランダム行列理論に関する研究で数多くの賞を受賞しています。彼は、対話型のスーパーコンピューティングの創始者で、Julia プロジェクトを立ち上げた 1 人です。Velvel Kahan から多くを学び、浮動小数点の数学の問題に関する数値解析コンソーシアムに積極的に参加しました。

この記事のレビューに協力してくれた技術スタッフの Paul Ringseth に心より感謝いたします。