区間代数と無限小と無限大

区間代数で無限大と無限小を統一的に扱う

このドキュメントのコードのテストに必要な拡張とライブラリは以下の通りです。

1. numeric-prelude
2. QuickCheck

{-# LANGUAGE TypeOperators        #-}
{-# LANGUAGE TypeFamilies         #-}
{-# LANGUAGE DataKinds            #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE RebindableSyntax     #-}
{-# LANGUAGE NoImplicitPrelude    #-}
import           Test.QuickCheck
import           NumericPrelude
import qualified Algebra.Field      as Field
import qualified Algebra.Ring       as Ring
import qualified Algebra.Additive   as Additive
import qualified Algebra.ZeroTestable as ZeroTestable
import qualified Number.Ratio       as Ratio
import qualified MathObj.Polynomial as Polynomial
import qualified MathObj.LaurentPolynomial as Laurent
import           Data.List

何をしたいか

どんな言語でも区間の計算で必ず現れる面倒に、閉区間と開区間の区別があります。ノートでは簡単に(a,b),(a,b],[a,∞)などと書けても、これほど柔軟な文法をサポートする実用プログラミング言語はなさそうです。[^1]

現状、区間代数を扱うHaskellパッケージはいくつかありますが、ほぼ全て次のような直和型を使います。

data Interval a = Open   a a
                | Closed a a
                | OpenClosed a a
                | ClosedOpen a a

でもこれだとcaseがあちこちに必要になって、コードも重複しがちだしめんどくさい、そもそも無限区間の扱い考えたらパターンがまた膨れあがって、憂鬱だ・・・というのが本記事のモチベーションです。

パフォーマンスを検討するために、GHC-8.2.1で導入された-XUnboxedSums拡張も調べてみました。

無限小と無限大

"無限小"(infinitesimal)dというシンボルを用意して,

open   = (0   , 1  )  -- 開区間
closed = (0-d , 1+d)  -- 閉区間
half   = (-1/d, 1  )  -- 無限区間 (-∞, 1)

のように書けたら嬉しくないでしょうか？これならつまらないパターンマッチを省けますが、いくつかの疑問も浮かびます。

1. ポリモーフィックに定義できるのか？
   どんな型なら無限小dを加えてももともとの型の意味を変えず拡張できるのか？

2. 無限小や無限大を四則演算することの意味は？
   例えばd+dや1/d^2は何を意味するのか(意味すべきか)？

3. もとの型に全順序がある場合、
   dで拡張された型にも、全順序は拡張できるのか？

モデル

公理を満たす構造をモデルといいます。任意の全順序集合は、全順序の公理を満たすモデルです。HaskellではOrdクラスとして順序比較のインターフェースが用意されています。 QuickCheckで公理を書きくだすと：

prop_reflexibity :: Ord a
                 => a -> Bool
prop_reflexibity a = a <= a

prop_antisymmetry :: Ord a
                   => a -> a -> Property
prop_antisymmetry a b = a <= b && b <= a
                   ==> a == b

prop_transitivity :: (Ord a)
                     => a -> a -> a -> Property
prop_transitivity a b c = a <= b && b <= c
                          ==> a <= c
-- totallityはOrdインスタンスで保証される

のようになります。これらの公理を満たすOrd a => aが全順序のモデルです。

∞を加えてみる

公理が満たされることを強調するため、以下OrdインスタンスではなくOrdモデルと呼びましょう。任意のOrdモデルに、「最大値」「最小値」を加えて、全順序を拡張することはできるでしょうか。

data BoundedOrd a = Min | X a | Max deriving (Eq, Ord)

instance Bounded  (BoundedOrd a) where
  minBound = Min
  maxBound = Max

これで任意の全順序の拡張ができました。この拡張は何度でも繰り返すことができます[^fn3]。

-- Here Nat is promoted to a kind by -XDataKinds
data Nat = Z | S Nat
type family   Power (f :: * -> *) (n :: Nat) (a :: *) :: *
type instance Power f Z a     = a
type instance Power f (S n) a = f (Power f n a)

-- Min < X Min < X (X Min) < .. < X (..(X a) ..) < .. < X Max < Max
type BoundedOrd_n n a = Power BoundedOrd n a

このようなモデルの拡張は、もっと複雑な公理系でも可能でしょうか？区間代数で使いたい(Ord a, Num a) => aをOrdやNumの公理を壊さないよう無限大を加えて拡張するには、どうしたらよいでしょうか。

順序体

区間代数では、四則演算と順序の関係が大事です。順序体OrdFieldクラスを定義してみましょう。
numeric-preludeのAlgebra.Field.Cクラスを使うと便利です。

class (Ord a, ZeroTestable.C a, Field.C a)
    => OrdField a where
newtype  OF a = OF a
   deriving (Eq, Ord,
             Field.C,
             Ring.C,
             Additive.C,
             ZeroTestable.C)
instance (Ord a, Field.C a, ZeroTestable.C a)
        => OrdField (OF a) where

prop_OrderedAddition :: OrdField a
                      => a -> a -> a -> Property
prop_OrderedAddition a b c =
   a <= b     ==>   a + c <= b + c

prop_LawOfSign :: OrdField a
                  => a -> a -> a -> Property
prop_LawOfSign a b c =
  0 <= a && 0 <= b    ==>   0 <= a * b

OrdFieldのモデルaに、無限小d(0 < x ==> 0 < d < x)を加えたモデルは多項式環 MathObj.Polynomial aの商体に、適当な順序を入れたものになっているはずです。適当なライブラリが見つからなかったので、便宜的にクラスQuotField aと書きます。

どのような順序か、調べてみましょう。(0/=1な)OrdFieldは整数(1 < 2=1+1
< ..)を含みます(つまり標数0)。上の公理から導きだせる、幾つかの性質を見てみましょう。詳しくはこちらを参照してください。

-- | 命題1
prop1 :: OrdField a => a -> a -> a -> Property
prop1 a b c =  (a < b) && (c > 0) ==> a*c < b*c

0 < d < 1から帰納法ですぐに次が言えます。

 0 < .. < d^n < .. d^3 < d^2 < d < .. < 1/n < ...  <  1/2  <  1

-- | 命題2
prop2 :: OrdField a => a -> Property
prop2 a =  a /= 0   ==>   0 < a^2

1/dは必ず無限大でしょうか？

-- | 命題3
-- 0 < 1/(a*b) に注意
prop3 :: OrdField a => a -> a -> Property
prop3 a b =
   sorted [0, a, b] ==> sorted [0, 1/b, 1/a]

なので、

0 < 1 < 2 .. < n < .. < 1/d < 1/d^2 < 1/d^3 < .. < 1/d^n ...

も言えます。

定理1

OrdFieldのモデルa、a係数の多項式環Polynomial aとその商体QuotField aに対し、次がなりたつ。

• aの順序を以下のようにQuotField aに拡張し、QuotField a をOrdFieldのモデルにできる。

instance OrdField a => Ord (Polynomial.T a) where
  x <= y  =  lc (y-x) > 0
    where
      -- the leading coefficient of f
      lc f = last (Polynomial.coeffs f)

newtype QuotField a = QF (Ratio.T (Polynomial.T a))
    deriving(Eq, Additive.C, ZeroTestable.C, Ring.C, Field.C)

instance OrdField a => Ord (QuotField a) where
  QF x <= QF y = let
     (a, b) = (numerator x, denominator x)
     (c, d) = (numerator y, denominator y)
     in  a * d <= b * c
instance OrdField a => OrdField (QuotField a)

• このような拡張はただ一つしかない。

これで任意のOrdFieldのモデルaに、無限小dを加えてOrdFieldの拡張モデルを作ることができました。aにシンボルdを加えて多項式環を作り、その商体に上のように順序を拡張すればできあがりです。

ローラン環

Rationalなどの順序体でなく、Integerなどの順序整域OrdIntegralDomainを考えるとどうなるでしょうか。numeric-preludeにはAlgebra.IntegralDomain.Cという型クラスが用意されています。順序整域をI、その商体をQと書くと、次の環拡大の図式を書けます。

  I
  /\
Q  I[d,1/d]
 \  /
 Q(d)

ここから上の定理1と同様に、I[d,1/d]の順序はただ一つ定まることが言えます(分母を払って比較するものしかない)。 I[d,1/d]をローラン環といいます。ローラン環の計算は商体の計算(ユークリッドの互助法またはグレブナー基底などを使う)よりも格段に易しく高速なので、区間の計算には望ましいでしょう。

-- | 0 < .. < d^n < .. d^3 < d^2 < d < .. < 1/n < ...  <  1/2  <  1
prop_laurenOrderInfinitesimal n =
   n > 0   ==>   sorted [0,  d^(n+1), d^n, 1/n]
-- | 0 < 1 < 2 .. < n < .. < 1/d < 1/d^2 < 1/d^3 < .. < 1/d^n ...
prop_laurenOrderInfinitude    n =
   n > 0   ==>   sorted [n, 1/d, 1/d^n, 1/d^(n+1)]

sorted :: Ord a => [a] -> Bool
sorted xs = xs == sort xs

dを一次の無限小、d^nをn次の無限小と呼びましょう。 二次の無限小d^2は、どんな整数nに対しても,n * d^2 < dが成りたちます。n次の無限小は、(n-1)次の無限小より「無限に小さい」というわけです。同様に、n次の無限大は(n-1)次の無限大より、「無限に大きい」と言えます。

HaskellのRationalやIntegerで考えると、Laurent RationalやLaurent Integerはそれぞれ順序体、順序整域になり、これはいわゆる超実数(hyperreal)の一部分になっています。

これで晴れて、無限小と無限大の正しい扱いができました。

例えば、任意の区間の補集合の、「無限に大きい」上下二つの部分区間を定義してみます。

type HyperRational = Laurent.T Rational
type HyperInterval = (HyperRational, HyperRational)

complement :: Integer
           -> HyperInterval
           -> (HyperInterval, HyperInterval)
complement n (a,b) = let
  ω = 1/d
  in ( (a - ω^n , a - d^n)
     , (b + d^n  , b + ω^n) )

これで、どんなx::HyperRationalも適当なcomplement nでカバーできますし、有限部分にしか興味がないならcomplement 1で間に合います。

• Ex.1
  (Ord a => Maybe a)の公理を書きだして、無限小と無限大を加えて拡張してください。

• Ex.2
  順序リスト(Ord a => [a])の公理を書きだして、無限小と無限大を加えて拡張してください。

• Ex.3
  Haskellの任意の順序モデルに対して、無限小と無限大を加えた順序モデルは存在するでしょうか。(単にOrdというだけでなく、他のクラスと順序の関係を保存すること)

パフォーマンスの考察

numeric-preludeによるローラン環の定義では、ボックスの使用は避けられず、Int64やFloat32に比べてパフォーマンスの低下が想定されます。

GHC 8.2.1から-XUnboxedSumsという拡張がサポートされました。この拡張を使って、パフォーマンスを向上させられるか、考えてみます

GHC 8.2.2のインストール

11月に、GHC8.2.2がリリースされ、各ディストリビューション向けにバイナリが用意されています。

Debian 9

GHC-8.2.2のDebian用ビルドはここにあります。

Ubuntu Zesty

$ # Should print the messages showing where the binaries are to be installed
$ sudo add-apt-repository ppa:hvr/ghc
$ sudo apt-get update
$ sudo apt-get install ghc-8.2.2
$ PATH=/opt/ghc/bin:$PATH

バイナリは/opt/ghc/bin以下にインストールされます。
インストール時のメッセージに注意してください。

UnboxedSums拡張

ghc-8.2.1でサポートされたUnboxedSums拡張を使うと、直和型の不要なポインタ使用を省けます。文法はUnboxedTupleの類似で、

 data T a1 .. ak =  (# t_1 | t_2 | ... | t_N #)

のように書けます。t_iは任意の型で、ボックス型でもアンボックス型でもどちらでも大丈夫です。C/C++/機械語プログラマであれば、"タグ付構造体"と言えばすぐにパフォーマンスのイメージが掴めるでしょう。タグは32-bitです。ユニオンではなく構造体を使うのは、ガーベッジコレクタとの連携のためと思われます。不要な重複フィールドを省略するアルゴリズムを実装してある、とドキュメントにはありますが、未確認です。

最小例を示します。

{-# LANGUAGE UnboxedSums #-}

data T a = (# Int | a #)

x0 :: T String
x0 = (# 0 | #) -- 最初の選択肢

x1 :: T String
x1 = (# | "abc" #) -- 二番目の選択肢

f :: T a -> Either Int String
f x = case x of
  (# i |   #) -> Left  i
  (#   | s #) -> Right s

main = do
    print $ f x0
    print $ f x1

ドキュメントはこちら。

この拡張を使い、ローラン環HyperInteger を最適化してみましょう。

type HyperInteger = (#
{- |  Default value for 'small' enough finite values. 
      For simplicity we will restrict ourselves to 32-bit
      arithmetics (except the most significant bit), so  
      there will be no loss of information in arithmetic
      operations. 
-}
   Int64# 
{- | Fallback for infinite or infinitesimal values -}

   | Laurent Integer 
   #) 

この定義によると、HyperIntegerのメモリ上の表現は

// This is a C code
struct HyperInteger {
  uint32_t tag;
  int64_t  first;
  box*     second;
};

のようになり、amd64/x86_64アーキテクチャでは(アラインメントによって)64*3ビットの構造体で表されます。メモリのサイズとしては最適ではないものの、IntelのCPUでは分岐予測が効きやすく、パフォーマンスが期待できそうです。

欠点としては、unboxed型の使用によってパラメトリックポリモーフィズムが失われることがあげられます。AssociatedTypes/TypeFamiliesなどの型関数を使い、アドホックに必要なケースを実装する必要があります。

もう一つの注意として、UnboxedSumはunionではなくstructとして表されるので、使用頻度の低い選択肢を含めてしまうと、キャッシュを浪費してしまうことです。選択肢は可能な限り少数にすることが望ましいでしょう。

整数、実数を部分として含む構造は無数にあります。そのような構造の計算では、大半の時間をInt#やFloat#として計算、オーバーフロー/アンダーフローやエラー時のみ他のボックスデータとして扱う、というのは自然で適用範囲の広い最適化といえます。将来的にはもっと使いやすい文法でGHCに実装されることを期待したいと思います。同様にInteger(巡回群の多項式)やRationalの高速化も期待できそうです。

結論

コンピュータの計算では、メモリは常に有限で、実数計算の精度は限られます。浮動点小数の計算では、NaNやエプシロン、無限大などの特別の値と性質がアドホックに定められ、足し算の結合法則すら保証されません。

ここでは、Haskellでは割合簡単に、またパフォーマンスへの犠牲を少なく、「正しい」無限小・無限大の扱いができることを示唆してみました。

ローラン環は大変高速で、実数の近似にも役立ちます。極限をとって形式的ローラン級数を考えれば体にもなります。Rationalほか商体の代替案として、プログラミング上もっと注目されてもいいのではないでしょうか。

[fn1]: 残念ながらDoubleやIntは、体や整域ではないので、任意の精度の計算が求められる場合は、 そのままは使えません。

補足：Haskellと型レベルの文法解析

Haskellでは、型クラスを使って型レベルのプッシュダウンオートマトンを実装するテクニックがあって、openPar openBracket a b closePar closeBracketなどの文脈自由なシンボルの並びを処理することはできますが、残念ながら記号は変数として扱えません(エラーメッセージも理解不能なくらい複雑になりがち)。