cl-num-utils: 数値計算ユーティリティ

概要

このページでは数値計算ライブラリあるcl-num-utilsについて簡単に説明します。

cl-num-utilsは大きなライブラリで、数学的に専門性の高い関数なども定められています。ただ、多くの場合、バリバリの数値計算(線形代数・行列計算)はBLASやLAPACKなどの専用のライブラリを用います。また、最近ではCPUではなくGPUで計算を行うことも増えています。GPUで計算する場合はnVidiaのCUDAというライブラリが必須です。Common Lispでバリバリの数値計算をやる場合、BLAS/LAPACK/CUDAに対応したMGL-MATを使うのが最も高速だと思いますが、MGL-MATを導入すると依存関係でLLAというライブラリが入ります。これがBLAS/LAPACKを呼び出すのですが、このLLAの実装にcl-num-utilsが使われています。Common Lisp製の数値計算ライブラリはいくつかありますが、cl-num-utilsが最も応用性が高いと判断したため、紹介します。

cl-num-utilsは多くのパッケージに分割され、名前空間が別れています。インポートせずに使用する場合はニックネームのclnuでシンボルにアクセスでき、インポートする場合は細かいパッケージを指定してインポートすることもできます。ただし、meansumなどの衝突を起こしやすそうな名前のシンボルも定義されているため、それらをインポートする場合はshadowing-importを使用してください。(mean, variance, medianの3種はAlexandriaと競合します。)
(asdf:load-system :cl-num-utils)
; => T

このページでは基本的に上の状態でサンプルを示します。

紹介するコンテンツ
  1. ユーティリティ
  2. 算術関数
  3. 要素毎の計算
  4. インターバル
  5. 行列
  6. チェビシェフの多項式
  7. 求根アルゴリズム
  8. 数値積分(ロンバーグ法)

ユーティリティ

cl-num-utilsにはcl-num-utils.utilitiesパッケージの中にいくつかの便利なオペレータが定められています。ここでは、汎用性が高いと思われるものをいくつか紹介します。

Function: within?(範囲内の判定)

3つの値を引数に取り、2番目の値が「1番目の引数以上、3番目の引数未満」に入っているかどうかを判定します。
(clnu:within? 1 2 3)
; => T

(clnu:within? 1 1 3)
; => T

(clnu:within? 1 3 3)
; => NIL

Function: as-double-float(double-flaot型変換)

数値計算は'double-float型で行うことが圧倒的に多いので、as-double-float関数で型変換を行うことができます。
(map 'vector #'clnu:as-double-float #(1 2.0s0 1/3))
; => #(1.0d0 2.0d0 0.3333333333333333d0)

Macro: with-double-floats (double-float型の束縛)

as-double-float関数を用いて型変換を行ってから値を束縛するマクロです。
(clnu:with-double-floats ((x 1)
                          (y 2.0s0)
                          (z 1/3))
  (vector x y z))
; => #(1.0d0 2.0d0 0.3333333333333333d0)

Function: binary-search(二分探索)

ソート済の数列ベクトルから二分探索法によってマッチするインデックスを返します。
;; 乱数による数列をベクトルで生成したいので iterate パッケージを使う
(asdf:load-system :iterate)
; => T

;; 乱数の数列ベクトルを生成
(iter (repeat 10) (collect (random 100) :result-type 'vector))
; => #(18 13 81 0 20 19 40 74 52 68)

;; ソートしておく
;; (* は前の評価結果を束縛したスペシャル変数)
(sort * #'<)
; => #(0 13 18 19 20 40 52 68 74 81)

;; 二分探索をする
(clnu:binary-search * 40)
; => 5

;; 2つ前の評価結果(配列)の中で、
;; 1つ前の評価結果(5)の要素番号の要素を取得する
(aref ** *)
; => 40

算術関数

cl-num-utilsにはいくつかの有用な算術関数が定められています。ここではcl-num-utils.arithmeticパッケージのオペレータを紹介します。

Function: same-sign?(符号の等価性判定)

same-sign?関数は可変長引数の符号が同じかを判定します。
(clnu:same-sign? 1 -2 3)
; => NIL

(clnu:same-sign? 1 2 3)
; => T

Function: square(2乗)

square関数は引数を2乗します。
(map 'vector #'clnu:square '(1 2 3 4 5))
; => #(1 4 9 16 25)

Function: abs-diff(距離)

abs-diff関数は2引数間の距離を求めます。
(clnu:abs-diff 2 -3)
; => 5

Function: log10, log2(対数)

ANSI Common Lisp標準のlogは対数の底がe(ネイピア数)であり、「自然対数」を意味しますが、第2引数を指定すると底を変更できます。log10log2はそのエイリアスとしての関数です。
;; cl-num-utils版
(clnu:log10 10000)
; => 4

;; ANSI Common Lisp標準版
(log 10000 10)
; => 4

Function: 1c(1からの減算)

ANSI標準の1-関数は(- number 1)という意味ですが、cl-num-utilsの1c(- 1 number)という意味です。
(clnu:1c 0.3)
; => 0.7

Function: divides?(割り切れるかの判定)

関数名の通り、割り切れるかどうかを判定します。割り切れる場合は商を、割り切れない場合はnilを返します。
(clnu:divides? 10 3)
; => NIL

(clnu:divides? 9 3)
; => 3

Function: numseq(数列の生成)

数列を生成します。Alexandriaのiotaと似ていて、iotaの方がメジャーですが、numseqも直感的です。
;; from 1 to 9
(clnu:numseq 1 9)
; => #(1 2 3 4 5 6 7 8 9)

;; from 1 to 9 length 5 => by 2 (= (9 - 1) / 5)
(clnu:numseq 1 9 :length 5)
; => #(1 3 5 7 9)

;; from 1 to 9 by 2 => length 5 (= (9 - 1) / 2 + 1)
(clnu:numseq 1 9 :by 2)
; => #(1 3 5 7 9)

;; from 1 to 9 result-type 'list
(clnu:numseq 1 9 :type 'list)
; => (1 2 3 4 5 6 7 8 9)

以下がAlexandriaのiotaの場合ですが、返り値がリスト固定なので大量の数列を生成する場合は不利かもしれません。
(alexandria:iota 5 :start 1 :step 2)
; => (1 3 5 7 9)

Function: ivec(整数ベクトルの生成)

numseqと似ていますが、整数のベクトルを生成することに特化したものです。使い方はiotaによく似ていますが、キーワード引数ではなくオプショナル引数を使用します。
;; from 0 below 10
(clnu:ivec 10)
; => #(0 1 2 3 4 5 6 7 8 9)

;; from 1 below 10
(clnu:ivec 1 10)
; => #(1 2 3 4 5 6 7 8 9)

;; from 1 below 10 by 2
(clnu:ivec 1 10 2)
; => #(1 3 5 7 9)

Generic Function: sum, product(加算・乗算)

sumは加算の縮約を、productは乗算の縮約を行いますが、どちらも総称関数であり、リスト・ベクトルなどのシーケンスに加え、行列をはじめとする多次元配列にも対応しています。
(clnu:sum #(1 2 3 4))
; => 10

(clnu:product #2a((1 2)(3 4)))
; => 24

Function: cumulative-sum, cumulative-product(累積和・累積)

sumproductの計算過程が必要であれば、cumulative-sumcumulative-productを使います。
(clnu:cumulative-sum #(1 2 3 4))
; => #(1 3 6 10) ;
;    10

(clnu:cumulative-product #(1 2 3 4))
; => #(1 2 6 24) ;
;    24

Function: l2norm(L2ノルム)

原点0からの距離の2乗の和の平方根を求めます。
;; (sqrt (+ 4 1 0 1 4))に等しい
(clnu:l2norm #(-2 -1 0 1 2))
; => 3.1622777

;; (sqrt (/ (+ 4 1 0 1 2) 5))に等しい
;; (今回はどちらも平均が0にしている)
(clnu:sd #(-2 -1 0 1 2))
; => 1.5811388300841898d0

Function: normalize-probabilities(発生確率の正規化)

複数の事象が発生するイベントがあるとき、各事象の発生回数をもとに発生確率を計算する関数です。サンプルを示した方がわかりやすいと思うので、サイコロの例で試します。

まず、繰り返しを使うので、iterateパッケージをロードします。また、ANSI標準の乱数ではなく、せっかくなのでメルセンヌ・ツイスタ法を使いたいので、MT19937もロードします。
(asdf:load-system :iterate)
; => T
(use-package :iterate)
; => T

(asdf:load-system :mt19937)
; => T
(shadowing-import 'mt19937:random)
; => T

サイコロを10000回降って、どの目が何回出たかを記録してみます。
(iter (repeat 10000)
      (with result = (make-array 6 :initial-element 0))
      (for r = (random 6))
      (incf (svref result r))
      (finally (return result)))
; => #(1682 1613 1683 1663 1688 1671)

どれも似たような数字ですが、これを正規化します。
;; * は変数。手前の式の評価結果を使うことができる。
(clnu:normalize-probabilities *)
; => #(841/5000 1613/10000 1683/10000 1663/10000 211/1250 1671/10000)

分数になってちょっとよく分からなくなったので、小数に変換してみます。
(map 'vector (lambda (x) (coerce x 'float)) *)
; => #(0.1682 0.1613 0.1683 0.1663 0.1688 0.1671)

正規化すると、どの目も16.13%〜16.88%(幅0.75%pt)で発生したことが分かります。ちなみに、100万回に増やしたところ、16.6354%〜16.7261%(幅0.09%pt)でした。許容範囲が1%ptなら1万回、0.1%ptなら100万回という感じでしょうか。

要素毎の計算

cl-num-utilsには配列の要素毎の計算を行う便利な関数が複数定められています。パッケージはcl-num-utils.elementwiseです。

Function: e+, e-, e*, e/(四則演算)

要素毎の四則演算は以下の通りです。
(clnu:e+ #2a((10 20) (30 40)) #2a((1 2)(3 4)))
; => #2A((11 22) (33 44))

(clnu:e- #2a((10 20) (30 40)) #2a((1 2)(3 4)))
; => #2A((9 18) (27 36))

(clnu:e* #2a((10 20) (30 40)) #2a((1 2)(3 4)))
; => #2A((10 40) (90 160))

(clnu:e/ #2a((10 20) (30 40)) #2a((1 2)(3 4)))
; => #2A((10 10) (10 10))

cl-num-utilsにはe+系オペレータの裏に引数の数が固定されたe2+e1-などのオペレータも定められています。e+e-は引数の個数によってこれらのオペレータを使い分けるように設計されているので、一般的にはe+系オペレータを使った方が便利ですが、引数の個数が1または2で明確ならe1-e2+を直接使う方が若干早いはずです。

Function: eexpt, eexp, elog, esqrt(指数・対数)

これらの関数も要素毎の処理を行いますが、複数の配列に対して処理をするのではなく、1つの配列内の要素に対して処理を行います。
(clnu:eexpt #2a((1 2)(3 4)) 2)
; => #2A((1 4) (9 16))

(clnu:esqrt *)
; => #2A((1 2) (3 4))

(clnu:eexp #2a((1 2)(3 4)))
; => #2A((2.7182817 7.389056) (20.085537 54.59815))

(clnu:elog *)
; => #2A((0.99999994 2.0) (3.0 4.0))

(clnu:elog #2a((1 2)(4 8)(16 32)) 2)
; => #2A((0 1) (2 3) (4 5))

Function: emax, emin(最大値・最小値)

ANSI Common Lisp標準のmax関数とmin関数は可変長引数をとるオペレータですが、cl-num-utilsのemaxeminは配列・リストを引数に取ります。
(clnu:emax #2a((3 1)(4 2)))
; => 4
(clnu:emin #2a((3 1)(4 2)))
; => 1

(clnu:emax '(3 1 4 2))
; => 4
(clnu:emin '(3 1 4 2))
; => 1

インターバル

cl-num-utilsには「インターバル」というオブジェクトがCLOS (Common Lisp Object System)で実装されています。これは、実数における範囲を示すものです。特徴としては、無限大も指定できることです。

パッケージはcl-num-utils.intervalです。

Function: interval(インターバルの生成)

interval関数で範囲オブジェクト(インターバル)を生成できます。生成されるインターバルは無限大の存在によって4種類に分けられます。
(clnu:interval -1 1)
; => #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [-1,1]>

(clnu:interval -1 :plusinf)
; => #<CL-NUM-UTILS.INTERVAL:PLUSINF-INTERVAL [-1,∞)>

(clnu:interval :minusinf 1)
; => #<CL-NUM-UTILS.INTERVAL:MINUSINF-INTERVAL (-∞,1]>

(clnu:interval :minusinf :plusinf)
; => #<CL-NUM-UTILS.INTERVAL:REAL-LINE (-∞,∞)>

種類によってカッコの付き方も違うのは仕様です。

Function: left, right(インターバルのアクセッサ)

インターバルオブジェクトのスロットにアクセスするには、leftrightを使います。
(defparameter *ivl* (clnu:interval -1 1))
; => *IVL*

(clnu:left *ivl*)
; => -1
(clnu:right *ivl*)
; => 1

これらはアクセッサですが、読み取り専用であるためスロットの変更はできません。

Function: open-left?, open-right?(無限大かの判定)

インターバルに無限大が含まれるかを判断するには、open-left?open-right?を使用します。
(clnu:open-left? *ivl*)
; => NIL
(clnu:open-right? *ivl*)
; => NIL

(clnu:open-left? (clnu:interval :minusinf 1))
; => T

Function: interval-length(インターバルの長さ)

インターバルの左右の距離を求めるにはinterval-length関数を用います。
(clnu:interval-length *ivl*)
; => 2

Function: interval-midpoint(インターバルの中間点)

インターバルの中間点・中心を求めるにはinterval-midpoint関数を使います。
(clnu:interval-midpoint *ivl*)
; => 0

Function: in-interval?(インターバルの範囲内判定)

ある数がインターバルに含まれるかどうかを判定するにはin-interval?関数を使います。
(clnu:in-interval? *ivl* 0.5)
; => T

(clnu:in-interval? *ivl* 3)
; => NIL

Function: split-interval(インターバルの分割)

インターバルを分割して複数のインターバルにするにはsplit-interval関数を使います。どのように分割するかの指定にはrelativespacerを使用します。
;; 相対的な長さで半分 + 残り
;;   => 半分 + 半分
(clnu:split-interval *ivl* (list (clnu:relative 1/2) (clnu:spacer)))
; => #(#<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [-1,0]>
;      #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [0,1]>>

;; 絶対的な長さで1/4 + 残り
;;   => 1/4 + 7/4
(clnu:split-interval *ivl* (list 1/4 (clnu:spacer)))
; => #(#<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [-1,-3/4]>
;      #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [-3/4,1]>)

Function: shrink-interval (インターバルの縮小・拡大)

インターバルは縮小したり、拡大したりすることができます。
;; 負の数だと拡大する
(clnu:shrink-interval *ivl* -1)
; => #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [-2,2]>
(clnu:shrink-interval *ivl* -2)
; => #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [-3,3]>

;; 正の数だと縮小する
(clnu:shrink-interval *ivl* 1/2)
; => #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [-1/2,1/2]>
(clnu:shrink-interval *ivl* 1/3)
; => #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [-2/3,2/3]>

;; relative も使用できる
(clnu:shrink-interval *ivl* (clnu:relative 1/3))
; => #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [-1/3,1/3]>

;; 左右を分けて伸縮できる
;; (左は1/2縮小、右は2拡大)
(clnu:shrink-interval *ivl* 1/2 -2)
; => #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [-1/2,3]>

Function: grid-in(グリッドを取得する)

インターバルにおいて等間隔に位置する点を「グリッド」と呼ぶとすると、そのグリッドを取得するのがgrid-in関数です。
(clnu:grid-in *ivl* 2)
; => #(-1 1)

(clnu:grid-in *ivl* 3)
; => #(-1 0 1)

(clnu:grid-in *ivl* 9)
; => #(-1 -3/4 -1/2 -1/4 0 1/4 1/2 3/4 1)

第2引数はグリッドの数で、間隔の数ではありません。そのため、最小は2です(左端と右端の2箇所)。

Function: subintervals-in(インターバルの等分割)

グリッドと同様に等分割でインターバルを生成する場合はsubintervals-in関数を使用できます。
(clnu:subintervals-in *ivl* 2)
; => #(#<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [-1,0.0)>
;      #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [0.0,1]>)

(clnu:subintervals-in *ivl* 4)
; => #(#<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [-1,-0.5)>
;      #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [-0.5,0.0)>
;      #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [0.0,0.5)>
;      #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [0.5,1]>)

こちらは間隔の数を指定します。

Generic Function: shift-interval(インターバルのシフト)

関数名の通り、インターバルを左右にシフトします。
(clnu:shift-interval *ivl* 2)
; => #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [1,3]>

(clnu:shift-interval *ivl* -2)
; => #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [-3,-1]>

行列

行列はcl-num-utilsの中心的な機能であり、パッケージは3つに別れています。
  • cl-num-utils.matrix-shorthand: 行列やベクトルを手軽に生成するためにオペレータ
  • cl-num-utils.print-matrix: print-matrixを中心とする出力オペレータ
  • cl-num-utils.matrix: transposeなどの処理オペレータ
このうち、上の2つは標準ではclnuパッケージにexportされていません。ただ、パッケージ名を付けて使用するのは長くて冗長なので、使用する場合はシンボルをインポートします。ここでは、この2つのパッケージについてはインポート済である前提でサンプルを示します。
(use-package :cl-num-utils.matrix-shorthand)
; => T
(use-package :cl-num-utils.print-matrix)
; => T

このページの最初に述べている通り、行列の計算は他の線形代数ライブラリを使うことが多いので、行列計算に関するオペレータは含まれていないのが特徴です。あくまでも行列を扱いやすくする、というのが主眼です。

ここではオペレータを絞って紹介しますので、興味があればソースコードを確認してください。

Function: vec (型付きベクトル)

vecは関数で、型指定ありのベクトルを生成します。要素は型に合うように変換されます。
(vec 'double-float 1 2 3 4 5)
; => #(1.0d0 2.0d0 3.0d0 4.0d0 5.0d0)

Macro: mx(型付き行列)

mxはマクロで、型指定ありの行列を生成します。要素は型に合うように変換されます。
(mx 'double-float (1 2) (3 4))
; => #2A((1.0d0 2.0d0) (3.0d0 4.0d0))

要素はリストの可変長引数で指定します。

ちなみに、行列の内容がリストのリストで指定された場合は以下のようなマクロを別途定義しておき、こちらを通すと使えるようになります。
(defmacro mx* (element-type rows)
  `(cl-num-utils.matrix-shorthand:mx ,element-type ,@rows))
; => MX*

(mx* 'double-float ((1 2)(3 4)))
; => #2A((1.0d0 2.0d0) (3.0d0 4.0d0))

どのように指定するかは好みに応じて使い分けてください。

Function: diagonal-mx(対角行列)

対角行列とは、左上から右下の対角成分以外は全て0である行列です。cl-num-utilsでは対角行列専用の構造体が定義されています。
(diagonal-mx 'double-float 1 2 3 4 5)
; => #S(CL-NUM-UTILS.MATRIX:DIAGONAL-MATRIX
;       :ELEMENTS #(1.0d0 2.0d0 3.0d0 4.0d0 5.0d0))

(aops:as-array *)
; => #2A((1.0d0 0 0 0 0)
;        (0 2.0d0 0 0 0)
;        (0 0 3.0d0 0 0)
;        (0 0 0 4.0d0 0)
;        (0 0 0 0 5.0d0))

array-operationsas-array総称関数がdiagonal-matrix構造体用に再定義されているため、as-arrayを通せば普通の行列(ANSI Common Lisp標準の2次元配列)に変わります。

Function: print-matrix(行列の表示)

print-matrix関数はとても便利な出力オペレータです。行列とストリームを指定すると、綺麗に表示してくます。
(print-matrix #2a((1 2 3)(4 5 6)) t)
;   1 2 3
;   4 5 6
; => NIL

(print-matrix #2a((1 2 3)(40000 50000 60000)) t)
;       1     2     3
;   40000 50000 60000
; => NIL

標準では右寄せされますが、:aligned?キーワード引数をnilにすると右寄せされなくなります。
(print-matrix #2a((1 2 3)(40000 50000 60000)) t :aligned? nil)
;   1 2 3
;   40000 50000 60000
; => NIL

行列をCSVで出力したければ、:paddingキーワード引数を","に、:indentキーワード引数を""にします。
(print-matrix #2a((1 2 3)(40000 50000 60000)) t 
              :aligned? nil
              :padding ","
              :indent "")
; 1,2,3
; 40000,50000,60000
; => NIL

Generic Function: diagonal-vector(対角成分の取得)

行列から対角成分を取得します。
(clnu:diagonal-vector #2a((1 2)(3 4)))
; => #(1 4)

この総称関数は(setf diagonal-vector)総称関数も合わせて定義されているため、アクセッサとして機能します。
(let ((a #2a((1 2)(3 4))))
  (setf (clnu:diagonal-vector a) #(10 40))
  a)
; => #2A((10 2) (3 40))

Function: diagonal-matrix(対角行列)

cl-num-utils.matrix-shorthandパッケージに定められたdiagonal-mx関数は、この関数を適用します。
(clnu:diagonal-matrix #(1 2 3))
; => #S(CL-NUM-UTILS.MATRIX:DIAGONAL-MATRIX :ELEMENTS #(1 2 3))

(aops:as-array *)
; => #2A((1 0 0) (0 2 0) (0 0 3))

diagonal-mxと異なり、引数は予めベクトルにしておく必要があり、また、型の変換も行われません。

Function: lower-triangular-matrix(下三角行列)

行列から下三角行列を抜き出して生成します。返り値は独自の構造体です。
(clnu:lower-triangular-matrix #2a((1 2 3)(4 5 6)(7 8 9)))
; => #<CL-NUM-UTILS.MATRIX:LOWER-TRIANGULAR-MATRIX
;      element-type T
;      1 . .
;      4 5 .
;      7 8 9>

print-object総称関数も定義されているため、構造体ですが独自の表示形式で表示されます。

三角行列などの特殊な行列はcl-sliceslice総称関数も定義されているため、手軽に切り出すことができます。
(cl-slice:slice * t 1)
; => #(0 5 8)

Function: upper-triangular-matrix(上三角行列)

上三角行列を抜き出します。
(clnu:upper-triangular-matrix #2a((1 2 3)(4 5 6)(7 8 9)))
; => #<CL-NUM-UTILS.MATRIX:UPPER-TRIANGULAR-MATRIX
;      element-type T
;      1 2 3
;      . 5 6
;      . . 9>
print-object, aops:as-array, cl-slice:sliceなどの総称関数が合わせて定義されています。

Generic Function: transpose(転置行列)

行列の転置を行います。

ANSI Common Lisp標準の2次元配列のほか、ここで紹介してきた三角行列にも対応しています。
(defparameter *a* #2a((1 2 3)(4 5 6)(7 8 9)))
; => *A*

(clnu:transpose *a*)
; => #2A((1 4 7) (2 5 8) (3 6 9))

(clnu:transpose (clnu:lower-triangular-matrix *a*))
; => #<CL-NUM-UTILS.MATRIX:UPPER-TRIANGULAR-MATRIX
;      element-type T
;      1 4 7
;      . 5 8
;      . . 9>

(clnu:transpose (clnu:upper-triangular-matrix *a*))
; => #<CL-NUM-UTILS.MATRIX:LOWER-TRIANGULAR-MATRIX
;      element-type T
;      1 . .
;      2 5 .
;      3 6 9>

チェビシェフの多項式近似

cl-num-utilsには「チェビシェフの多項式」に関する関数も定められており、関数近似を行うことができます。チェビシェフの多項式についてはこちらのサイトなどを参考にしてみてください。

このページでは、cl-num-utilsのテストコードに掲載された式の一つである y = x / (x + 4) のチェビシェフ多項式近似を行うサンプルを示します。チェビシェフの多項式はcl-num-utils.chebyshevパッケージに収録されています。

Function: chebyshev-roots, chebyshev-root(チェビシェフ多項式の根)

チェビシェフ多項式の根、すなわち0に等しくなるときのxの値を求めます。

チェビシェフ多項式はこちらのサイトでも計算したり、グラフ化できますので、0になっていそうか確かめてみてください。
(clnu:chebyshev-roots 3)
; => #(-0.8660254037844387D0 -6.123233995736766D-17 0.8660254037844387D0)

(clnu:chebyshev-root 3 2)
; => 0.8660254037844387D0

Function: chebyshev-approximate(チェビシェフ多項式による近似)

ある関数の特定の区間において、チェビシェフ多項式を用いて近似関数を得ることができます。

近似化したい元の関数が y = x / (x + 4) だとします(cl-num-utilsのテストコードの例1)。この時、x = 1から x = 10までの値を求めておきます。
(defun real-fn1 (x)
  (float (/ x (+ 4 x))))
; => REAL-FN1

(defparameter *real-value1*
  (map 'vector #'real-fn1 (clnu:numseq 1 10)))
; => *REAL-VALUE1*

*read-value1*
; => #(0.2 0.33333334 0.42857143 0.5
;      0.5555556 0.6 0.6363636 0.6666667
;      0.6923077 0.71428573)

次に、元の関数から15次元のチェビシェフ多項式の近似関数を得て、その関数で近似値を取得してみます。)

(defparameter *app-fn1*
  (clnu:chebyshev-approximate #'real-fn1
                              (clnu:interval 1 10) 15))
; => *APP-FN1*

(defparameter *app-value1*
  (map 'vector *app-fn1* (clnu:numseq 1 10)))
; => *APP-VALUE1*

これだけでは近似できているかどうか分からないので、差分を取って比較してみます。
(defparameter *diff1*
  (map 'vector #'clnu:abs-diff *real-value1* *app-value1*))
; => *DIFF1*

(every #'(lambda (x) (< x 1e-5)) *diff1*)
; => T

1e-50.00001ですので、誤差は十分に小さいと言えます。

Function: chebyshev-regression(チェビシェフ多項式の係数)

Common Lispは関数をラムダ式で直接表現できるため、前節で得られる近似式は関数の実体そのものです。ただ、関数の実体ではなく、係数が欲しい場合もあります。

その場合はchebyshev-regression関数を使います。
(clnu:chebyshev-regression #'real-fn1 15)
; => #(-0.03279555898864451d0 0.26236447190915607d0 -0.03332465729595915d0
;      0.004232786458517312d0 -5.376343721793371d-4 6.828851891812111d-5
;      -8.673779165265424d-6 1.1017144053308373d-6 -1.3993607700009866d-7
;      1.777421220205966d-8 -2.2576206755905066d-9 2.867554582675069d-10
;      -3.64223170142471d-11 4.6253204969796496d-12 -5.776545908275921d-13)

係数の一覧を見ると、後半の係数はほとんど0に近くなっており、項の説明力にあまり寄与していないことが分かります。テストコードが15を採用していたので15にしたのですが、試して見ると10でも1e-5の範囲に誤差が収まりました。

求根アルゴリズム

cl-num-utilsには単一根を求めるアルゴリズムも実装されています。

求根アルゴリズムを使って、2の平方根((sqrt 2))を近似的に求めてみます。なお、求根アルゴリズムはcl-num-utils.rootfindingパッケージに収録されています。

Function: root-bisection(単一解の求根)

2の平方根は「ひとよ ひとよに ひとみごろ(1.41421356)」として覚えると思いますが、数値解析的に求める場合、関数を想定する必要があります。

まず、√2 を求めたい訳ですから、これを x とおきます( x = √2 )。次に、両辺を2乗します( x^2 = 2 )。最後に、根を求めたいので、右辺が 0 になるように移項します( x^2 - 2 = 0 )。こうして得られる左辺が、今回の関数です ( f(x) = x^2 - 2 )。

関数を想定すれば、あとは解析的に求める際の範囲を指定します。試しに0から10までを指定します。
(clnu:root-bisection #'(lambda (x) (- (expt x 2) 2))
                     (clnu:interval 0 10))
; => 1.4141845703125d0
;    -8.200109004974365d-5
;    T
;    #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [0,10]>

root-bisection関数は第1引数に関数を、第2引数に範囲を指定します。すると範囲内の値を適用して関数を評価し、値が0に近くなるかどうかを調べます。返り値は1つ目が求めた解、2つ目がその解を関数に適用した場合の計算値、3つ目がその計算値が許容誤差の範囲に収まっているかどうか、4つ目が解析に使用した範囲です。

標準では許容誤差の範囲が小数点3〜4桁程度しかないため、1.414までしか一致していません。もう少し許容誤差を狭めるには、*rootfinding-epsilon*スペシャル変数を変更します。
(setf clnu:*rootfinding-epsilon* 1.0d-8)
; => 1.0D-8

(clnu:root-bisection #'(lambda (x) (- (expt x 2) 2))
                     (clnu:interval 0 10))
; => 1.4142135623842478d0
;    3.154454475406965d-11
;    T
;    #<CL-NUM-UTILS.INTERVAL:FINITE-INTERVAL [0,10]>

これで1.4142135623まで一致しました。

数値積分

cl-num-utilsにはロンバーグ法による数値積分を行う関数も定義されています。

数値積分に関してはこちらのサイトが参考になります。

Function: romberg-quadrature(数値積分)

ある関数の特定の範囲を積分した値を返します。ロンバーグ法を使うのに適切かどうかはさておき、ここでは f(x) = x^2 + 2x + 1 という2次関数の1から2までの範囲の数値積分をしてみます。
(clnu:romberg-quadrature #'(lambda (x) (+ (expt x 2) (* 2 x) 1))
                         (clnu:interval 1 2))
; => 6.333333333333333d0
;    12

こちらの高精度計算サイトでも計算できますが、おそらく返り値の2つ目はステップ(繰り返しの回数)だと思います。(ソースコードのdocstringには明示されていませんが。)

0 件のコメント :

コメントを投稿