Haskellのおおよその数の比較

おそらく誰もが、精度が限られている計算では、数学的に等価な2つの式が互いに等しくない可能性があることを知っています。 たとえば、Haskellで計算するときの次の明白な数学的な同等性は、予想外に偽であることが判明しました。

ghci> 3 * sqrt(24 ^ 2 + 16 ^ 2) == sqrt(72 ^ 2 + 48 ^ 2) False
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

この違反の理由は、この等式の式がおよそのみ計算されることです。

 ghci> 3 * sqrt(24 ^ 2 + 16 ^ 2) 86.53323061113574 ghci> sqrt(72 ^ 2 + 48 ^ 2) 86.53323061113575 ghci> sqrt(72 ^ 2 + 48 ^ 2) - 3 * sqrt(24 ^ 2 + 16 ^ 2) 1.4210854715202004e-14
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

ここの違いは最後の（14番目！）10進数の場所にありますが、これは比較を偽にするのに十分です。



この問題はよく知られていますが、プログラマーはほとんど注意を払いません。 第一に、この種の比較は数値的手法の狭い領域でのみ生じると考えられ、第二に、平等の違反は非常にまれにしか起こらないと考えられています。 結局のところ、両方とも完全に真実ではありません。 整数座標を使用してベクトルの長さを計算するための関数を実装する必要があるときに、特定のケースが発生しました。 同時に、単体テストでは、 QuickCheckパッケージのツールが使用されます。これにより、ベクターの長さのスケーリングの不変条件に違反するケースがかなり迅速に見つかりました。 これは、テスト中に違反が発見された唯一の不変条件とはほど遠いことに注意してください。



疑問が生じます：限られた精度での計算の結果として得られた2つの数値のおおよその等価性の検証を記述する最も簡単な方法は何ですか？ Haskellでこの問題を解決するには、別の比較演算子（たとえば〜=）を定義するだけで十分です。これは、通常の等価演算子と同じ方法で使用されます。 このような演算子の実装を検討することを提案します。この演算子は、かなり単純なCircaモジュールとして配置できます。







2つの数値のおおよその比較で最初に頭に浮かぶのは、比較した数値の差の絶対値を計算し、それが特定のしきい値を超えているかどうかを確認することです。

 ghci> abs(sqrt(72 ^ 2 + 48 ^ 2) - 3 * sqrt(24 ^ 2 + 16 ^ 2)) < 1e-12 True
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

もちろん、このようなソリューションは完全に機能します。 しかし、彼には2つの重大な欠陥があります。 まず第一に、このレコードを読んでいるとき、2つの数値の等価性（近似的であっても）をチェックしていることはまったく明らかではありません。 さらに、予測の不正確さを示すために、「マジック」番号1e-12を使用する必要がありました。 もちろん、これらの不便とのそのような比較の少数で、人は条件に達することができました。 しかし、不変式の数が10単位で測定される場合、それらを記述するためのより単純で明確な方法を取得したいと思います。 近似比較の2桁演算子を通常の等号と同じ方法で使用するコードは、たとえば次のようになります。

 sqrt(72 ^ 2 + 48 ^ 2) ~= 3 * sqrt(24 ^ 2 + 16 ^ 2)
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

残念ながら、標準のHaskell構文はそのようなステートメントを提供しません。 ただし、言語自体は、標準の構文の一部であるかのように、新しい演算子を独自に導入する素晴らしい機会を提供します。



最初に、概算の比較レコードを削減するcircaEq関数を定義します

 circaEq :: (Fractional t, Ord t) => t -> t -> t -> Bool circaEq txy = abs (x - y) < t
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

これで、例は少し短くなりましたが、それほど明確ではありません。

 ghci> circaEq 1e-12 (sqrt(72 ^ 2 + 48 ^ 2)) (3 * sqrt(24 ^ 2 + 16 ^ 2)) True
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

まだ多くの不必要な情報があります。引数を分離するには、括弧で囲む必要があり、最も重要なことは、以前と同様に、比較の正確さを示す必要があります。 余分なパラメーターを取り除くために、 Data.Fixedモジュールで使用されたトリックを使用して、固定された精度で数値を表します。 一連の二重関数を定義します。それぞれの関数は、事前に決定されたエラーと数値を比較します。 実際には、最も一般的なエラー値に対してこれらの関数のうち7つだけを決定すれば十分であることがわかります。

 picoEq :: (Fractional t, Ord t) => t -> t -> Bool picoEq = circaEq 1e-12 infix 4 `picoEq` nanoEq :: (Fractional t, Ord t) => t -> t -> Bool nanoEq = circaEq 1e-9 infix 4 `nanoEq` microEq :: (Fractional t, Ord t) => t -> t -> Bool microEq = circaEq 1e-6 infix 4 `microEq` milliEq :: (Fractional t, Ord t) => t -> t -> Bool milliEq = circaEq 1e-3 infix 4 `milliEq` centiEq :: (Fractional t, Ord t) => t -> t -> Bool centiEq = circaEq 1e-2 infix 4 `centiEq` deciEq :: (Fractional t, Ord t) => t -> t -> Bool deciEq = circaEq 1e-1 infix 4 `deciEq` uniEq :: (Fractional t, Ord t) => t -> t -> Bool uniEq = circaEq 1 infix 4 `uniEq`
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

これらの関数はいずれも、2桁の比較演算子として使用できます。次に例を示します。

 ghci> sqrt(72 ^ 2 + 48 ^ 2) `picoEq` 3 * sqrt(24 ^ 2 + 16 ^ 2) True
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

あとは、少し砂糖を加えるだけです：

 (~=) :: (Fractional t, Ord t) => t -> t -> Bool (~=) = picoEq infix 4 ~=
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

必要なものを取得します。

 ghci> sqrt(72 ^ 2 + 48 ^ 2) ~= 3 * sqrt(24 ^ 2 + 16 ^ 2) True
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

別の重要な質問に答えます：近似比較のために、なぜいくつかの異なる関数が必要なのですか？ ピコ精度との比較は十分ではありませんか？ 実のところ、十分ではありません。 同じQuickCheckパッケージを使用して、適切な反例が見つかります。

 ghci> sqrt(5588 ^ 2 + 8184 ^ 2) ~= 44 * sqrt(127 ^ 2 + 186 ^ 2) False ghci> sqrt(5588 ^ 2 + 8184 ^ 2) `nanoEq` 44 * sqrt(127 ^ 2 + 186 ^ 2) True
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

明らかに、必要な精度のレベルは、使用する必要がある数値の規模に大きく依存します。 そのため、Circaモジュールは近似比較の演算子だけでなく、同義語になり得る一連の関数もエクスポートします。 アプリケーションがピコ精度の使用を受け入れない場合、必要な関数をインポートし、それに応じて近似比較演算子を定義できます。 誰かが近似比較演算子の異なる表記を好む場合、同じことができます。