Racket(Scheme) によるKrasnosel'skii-Mann iterations の数値実験

Go back

契約を用いたRacket(Scheme)プログラミングと，それらを使用した数値実験手法の入門的解説を行います． ここでは例として，非拡大写像$T$ の不動点$u\in{}F(T) \left(=\left\{x:Tx=x\right\}\right)$ を探索するプログラムを作成します． アルゴリズムとしては，次に示すKrasnosel'skii-Mann iterations を実装することにより，不動点近似を行います．

実装には，Racket¹⁾ を使用します． また，行列表現に用いるデータ構造，およびそれらに対する各種演算には，Math Library[(N. Toronto, J. A. Søgaard: Math Library. Racket Documentation: (2013).)] を使用します． 読者は，R⁵RS仕様による基礎的なSchemeプログラミングを習得していることを想定します．

今回は，より数学表現に則した形でプログラムを実装したいと思います． そのために必要となるいくつかの文法，特にScheme 標準からRacket により独自拡張された文法に関して解説します．

契約とは，定義する値に対してある一定の制約(保証)を与えるものです． 例えば，ある変数$x$が実数値であることを保証することや，ある関数$f$が$n$次元ベクトル$\mathbb{R}^n$から$m$次元ベクトル$\mathbb{R}^m$への写像であることを保証することができます．

Racketにおける契約の表現形態は種々存在します． 例えば，述語は契約として使用することができます． 述語とは，任意の1つの値を受け取り，真偽値を返す関数です． 例えば，関数number?は述語の1種で，与えられた引数が数値であれば真，そうでなければ偽を返します． この関数number?を変数$x$に対する契約として用いると，変数$x$が数値であることを保証することができます．

契約付き変数定義は，契約としてその変数が取りうる値の制約を与えて，変数を定義します． 以下に，契約付き変数定義の記述形式を載せます．

(define/contract <変数名> <契約> <初期値式>)

<変数名>には，この契約付き変数定義において定義する変数の名前を指定します． <契約>には，<変数名>に指定した名前の変数に対して，その変数が取りうる値の制約を表現した契約を指定します． <初期値式>には，<変数名>に指定した名前の変数に対して，初期値を与える式を指定します．

例として，契約付き変数定義を用いて，初期値$3$を取る変数$x\in\mathbb{N}$を定義してみましょう． <変数名>がx，<初期値式>が3であることは自明なので，ここに指定すべき<契約>を考えます． ここで，変数$x$は自然数です． Racketにおいては，自然数かどうかを判定する述語exact-positive-integer?が定義されています． (ここでは「自然数」を，「正確な正整数」と考えます．$\mathbb{N}:=\left\{x\in\mathbb{Z}:x>0\right\}$．) 述語は契約として使用することができるので，ここでは次の記述，

(define/contract x exact-positive-integer? 3)

により，初期値$3$を取る変数$x\in\mathbb{N}$を定義することができます．

契約付き変数定義により変数定義を行うと， これにより定義した変数に対して契約に反する操作を行った際に，エラーとして検出することができます． 例えば，先の例における変数$x$に自然数$5$を代入することは，

(set! x 5)

により行うことができます． しかし，自然数ではない値$0$を代入するような次のコード，

(set! x 0)

を実行すると，エラーとなります．

契約付き関数定義は，契約として定義域および値域を与えて(引数と戻り値に対して契約を与えて)，関数を定義します． 以下に，契約付き関数定義の記述形式を載せます．

(define/contract (<関数名> <引数>...)
  (-> <引数契約>... <値域契約>)
  <関数実装>)

この形式は，名前が<関数名>である関数を定義します． これにより定義される関数は，<引数>…にて定義した引数を受け取り，それらを用いて実装される式<関数実装>によって，その戻り値が計算されます． 契約は，受け取る各引数，およびそれらによって計算される戻り値に対して，制約をかけることができます． 各引数に対する契約は，定義した引数の個数分だけ順に，<引数契約>にて指定します． 戻り値に対する契約は，一連の契約を表すリストの最後の項<値域契約>にて指定します．

例として，実数$x,y\in\mathbb{R}$に対する，$\mathbb{R}$における通常の距離$d:\mathbb{R}^2\to\mathbb{R}_+$， $$ d(x,y):=\left|x-y\right| $$ を，契約付き関数定義により定義してみましょう．

(define/contract (d x y)
  (-> real? real? (not/c negative?))
  (abs (- x y)))

real?は実数かを判定する述語であり，negative?は負の実数かを判定する述語です． ここで，関数$d$は，2つの実数$x$，$y$を受け取り，非負実数を返す関数です． したがって，否定の契約を求める関数not/cを用いて，非負実数を保証する契約(not/c negative?)が構築できます． なお，関数absは実数の絶対値を計算する関数であり，関数-は減算を計算する関数です．

契約付き関数定義により定義された関数も，契約付き変数定義同様に， これにより定義した関数に対して契約に反する操作を行った際，または契約に反する結果が得られた際に，エラーとして検出することができます． 具体的には，関数呼び出し時に与えられた引数に関して，契約に反する操作を検出することができます． また，計算結果としての関数の戻り値に関して，契約に反する結果を検出することができます． 例えば，先に定義した関数$d$について，その引数に複素数を指定して呼び出すような式，

(d 3+2i 5)

を実行すると，エラーとなります．

Racketで行列およびその演算を表現する方法はいくつか考えられますが，ここでは線形代数関数群であるmath/matrixライブラリを用いましょう． math/matrixライブラリを利用するためには，次の式を評価すればよいです．

(require math/matrix)

今回は，行列のなかでも，特にn行1列行列，いわゆる列ベクトルが使用できればよいです． 一般の行列を生成する命令以外に，math/matrixライブラリでは列ベクトルを生成する命令col-matrixが定義されています． 例えば，列ベクトル$\left(1, 2, 3\right)^\mathrm{T}$を生成するためには，次の式を記述します．

(col-matrix (1 2 3))

命令col-matrixは，(col-matrix (<値>…))形式で記述し，<値>…の部分に生成したい列ベクトルの各要素を順に記述していきます．

行列に対しては，いくつかの演算が関数として定義されています． 代表的なものとしては，和matrix+，差matrix-，積matrix*が定義されています． これらは，2つの行列を引数に取り，その(線形代数学で定義される一般的な)計算結果となる行列を返します． また，2つの行列に対する内積は関数matrix-dotにより求めることができ，行列のノルムは関数matrix-normにより求めることができます．

例えば，2つの$n$次元列ベクトル$x,y\in\mathbb{R}^n$のノルムによる距離$d(x,y):=\|x-y\|$は，

(define/contract (d x y)
  (-> col-matrix? col-matrix? (not/c negative?))
  (matrix-norm (matrix- x y)))

により定義することができます． ここでは，列ベクトルかどうかを判定する述語としてはcol-matrix?を使用することができるので，これを用いて引数に契約を与えました． (ちなみに，一般の行列かどうかを判定する述語としてはmatrix?を使用することができます。)

(書きかけ.)

#lang racket
 
(require math/matrix racket/stream plot)
 
(define/contract (km-seq t a-seq x0)
  (-> (-> col-matrix? col-matrix?) stream? col-matrix? stream?)
 
  (define/contract (iter a x)
    (-> (and/c (>=/c 0) (</c 1)) col-matrix? col-matrix?)
    #:freevar t (-> col-matrix? col-matrix?)
    (matrix+ (matrix-scale x a)
             (matrix-scale (t x) (- 1 a))))
 
  (let lp ((a-seq a-seq) (x x0))
    (stream-cons x (lp (stream-rest a-seq) (iter (stream-first a-seq) x)))))
 
 
(define/contract alpha stream? (sequence->stream (in-cycle '(0.5))))
 
(define/contract (rotate90 x)
  (-> col-matrix? col-matrix?)
  (matrix* (matrix ((0 -1) (1 0))) x))
 
 
(let ((x-seq (km-seq rotate90 alpha (col-matrix (100 50)))))
  (plot
    (lines
      (for/list ((x (in-stream x-seq))
                 (i (in-range 30)))
        (vector i (matrix-norm (matrix- (rotate90 x) x)))))
    #:x-label "Number of iterations" #:y-label "||Tx-x||"))