[iOS 10] AccelerateフレームワークのQuadratureを用いて数値積分を行う

データ分析

#iPhone

#数学

#統計

#iOS

yad

2016.09.14

この記事は公開されてから1年以上経過しています。情報が古い可能性がありますので、ご注意ください。

はじめに

iOS 10では数値計算を主に行うフレームワークであるAccelerateフレームワークに数値積分のためのQuadrature APIが追加されました。本記事はその紹介になります。

数値積分について

数値積分は何らかの定積分をコンピュータ等を用いて数値的に求める手法のことをいいます。ここで言う「数値的に求める」というフレーズはよく「解析的に求める」というフレーズと対比されますが、私は次のように捉えています。

解析的に求める: 積分を式変形を用いて積分が現れない形まで展開し、展開の結果に値を代入するだけで値が求められる状態にすること
数値的に求める: コンピュータ等の計算機を用いて計算結果を近似的に求めること

積分の対象になる関数についてはその不定積分が容易に求められないことも多く、また不定積分が存在しないような被積分関数についても特定の積分区間では定積分を求めることができるものの、それ以外の積分区間では定積分が解析的に求められないことがあります。

例えばガウシアン

[latex] p(x) = \frac{1}{\sqrt{2 \pi}}e^{-\frac{x^2}{2}} [/latex]

の定積分は\([-\infty, \infty]\)の積分区間では1となりますが、積分区間が有限区間\([a, b]\)になった場合は解析的には求められず、数値的に求めることが必要になります。

数値積分の手法としてはいろいろなアルゴリズムが提案されていますが、代表的なものとしては次のようになります。

ニュートン・コーツの公式
ガウス求積法
モンテカルロ法

iOS 10で導入されたQuadratureはこのうちガウス求積法をより発展させたガウス＝クロンロッド求積法が用いられています。使い方の概要をガウシアンの定積分を例にとって確認していきます。

Swiftによるコード例

環境情報

Xcode 8 GM

import Foundation
import Accelerate

struct FunctionArrays {

    private static func writeWithFunction(n: Int, xs: UnsafePointer<Double>, ys: UnsafeMutablePointer<Double>, f: (Double) -> Double) {
        for i in 0..<n {
            let xi = xs[i]
            ys[i] = f(xi)
        }
    }

    static let gaussian: quadrature_function_array = { arg, n, x, y in
        writeWithFunction(n: n, xs: x, ys: y) { xi in exp( -xi * xi / 2.0 ) / sqrt(2 * Double.pi) }
    }
}


struct IntegralSample {

    static func runSample() {

        var gaussian = quadrature_integrate_function(fun: FunctionArrays.gaussian, fun_arg: nil)

        var options = quadrature_integrate_options()
        options.integrator = QUADRATURE_INTEGRATE_QAGS
        options.abs_tolerance = 1.0e-8
        options.max_intervals = 20

        let gaussianResult = quadrature_integrate(&gaussian, 0, 1.0, &options, nil, nil, 0, nil)

        print("gaussian result : \(gaussianResult)")
    }
}

順に上からコードを見ていきます。

l2

Accelerateフレームワークをインポートします。

l4

名前空間的にStructを用いるため、FunctionArraysを定義し、その下にstaticな関数を宣言していきます。

l6

この関数は積分区間の区切り点の個数nと、区切り点の入力xsをもとに区切り点における被積分関数の値を出力として可変ポインタysを書き換える処理を記述しています。被積分関数は最後の引数fとして表わせ、この関数のクライアントが自由に被積分関数を渡せるようになっています。

l13

quadrature_function_arrayは実際には次のようなタイプエイリアスです。

@convention(c) (UnsafeMutableRawPointer!, Int, UnsafePointer<Double>, UnsafeMutablePointer<Double>) -> Swift.Void

@convention(c) はCの関数ポインタをSwiftから扱えるようにしたものです。この型に対してはクロージャリテラルをSwiftからシームレスにあてがうことが可能です。ただ、@convention(c)のつけられた関数オブジェクトに対してはObjective-CのBlocksに互換性のあるconvention(blocks)やSwiftのデフォルトの関数オブジェクトとはことなり、外部変数をキャプチャする機能がないです。このため、何らかの参照や外部変数をキャプチャする処理は書けません。書くと次のようなエラーが出ます。

A C function pointer cannot be formed from a closure that captures context

l14

この行では実際の被積分関数としてガウシアンを指定して可変ポインタを書き換える処理をプライベート関数に委譲しています。外部から使う際のインターフェイス的に関数オブジェクトを渡してquadrature_function_arrayを生成する仕組みを生成したかったですが、そうすると、quadrature_function_arrayを生成する部分で先ほどのキャプチャのエラーが出てしまうので致し方無いところです。

l21

先ほど定義したロジックで生成されたquadrature_function_arrayをもとに数値積分を行い、結果を集中力する処理を行う関数です。

l23

ガウシアンに対応するquadrature_integrate_functionのstructを生成します。

l25-28

数値積分の計算に対する設定を表すquadrature_integrate_optionsを定義します。設定しているプロパティの内容は次の通りです。

#integrator : 数値積分のアルゴリズムを表します
#abs_tolerance : 絶対許容誤差
#max_intervals : 数値積分の際に積分区間を小区間に区切る個数の最大値です

上記コードのIntegralSample.runSample()をコードのどこかで実行すると下記の出力をログに得ます

l30

quadrature_integrate関数で実際に数値積分を行います。第一引数に被積分関数に対応するquadrature_function_arrayをポインタで渡し、第二、第三引数に積分区間の範囲を指定します。第四引数には積分の際のquadrature_integrate_optionsをポインタで渡します。

それ以降の引数については数値積分に失敗した際のエラーステータスをポインタ経由で取得したり、積分の際に用いる計算空間のサイズの指定などを行えます。

この関数呼び出しの結果は下記の数式に対応しており

[latex] \frac{1}{\sqrt{2 \pi}} \int^1_0 e^{-\frac{x^2}{2}} dx [/latex]

上記コードの計算結果はこちらになります。