差分

このページの2つのバージョン間の差分を表示します。

--- kaz:km-iter [2014/01/10 00:14] – [行列表現とその演算] kaz
+++ kaz:km-iter [2014/05/16 22:26] – [参考文献] kaz
@@ 行 1: / 行 1: @@
+====== 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 を実装することにより，不動点近似を行います．
+<WRAP box>
+//**__Krasnosel'skii-Mann iterations[(高橋 渉: **[[http://www.ybook.co.jp/noco.htm|非線形・凸解析学入門]]**. 横浜図書: 126--129 (2005).)]__**//
+$H:\text{Hilbert Space}$, $C\subset{}H: \text{nonempty}, \text{closed}, \text{convex}$, $T:C\to{}C,\text{nonexpansive}$, $F(T)\not=\emptyset$.
+$\left\{\alpha_i\right\}\subset[0,1)$, $\sum_{i=1}^\infty\alpha_i\left(1-\alpha_i\right)=\infty$.
+$x_1\in{}C$.
+$$ x_{i+1}=\alpha_ix_i+\left(1-\alpha_i\right)Tx_i\quad\left(i=1,2,\cdots\right) $$
+</WRAP>
+実装には，Racket((PLT Design Inc.: **[[http://racket-lang.org/|The Racket Language]]**.)) を使用します．
+また，行列表現に用いるデータ構造，およびそれらに対する各種演算には，Math Library[(N. Toronto, J. A. Søgaard: **[[http://docs.racket-lang.org/math/|Math Library]]**. Racket Documentation: (2013).)] を使用します．
+読者は，R<sup>5</sup>RS仕様による基礎的なSchemeプログラミングを習得していることを想定します．
+===== Racket(Scheme) の文法 =====
+今回は，より数学表現に則した形でプログラムを実装したいと思います．
+そのために必要となるいくつかの文法，特にScheme 標準からRacket により独自拡張された文法に関して解説します．
+==== 定義と契約 ====
+契約とは，定義する値に対してある一定の制約(保証)を与えるものです．
+例えば，ある変数$x$が**実数値であることを保証**することや，ある関数$f$が**$n$次元ベクトル$\mathbb{R}^n$から$m$次元ベクトル$\mathbb{R}^m$への写像であることを保証**することができます．
+Racketにおける契約の表現形態は種々存在します．
+例えば，述語は契約として使用することができます．
+述語とは，任意の1つの値を受け取り，真偽値を返す関数です．
+例えば，関数''number?''は述語の1種で，与えられた引数が数値であれば真，そうでなければ偽を返します．
+この関数''number?''を変数$x$に対する契約として用いると，変数$x$が数値であることを保証することができます．
+**契約付き変数定義**は，契約としてその変数が取りうる値の制約を与えて，変数を定義します．
+以下に，契約付き変数定義の記述形式を載せます．
+<code scheme>
+(define/contract <変数名> <契約> <初期値式>)
+</code>
+''<変数名>''には，この契約付き変数定義において定義する変数の名前を指定します．
+''<契約>''には，''<変数名>''に指定した名前の変数に対して，その変数が取りうる値の制約を表現した契約を指定します．
+''<初期値式>''には，''<変数名>''に指定した名前の変数に対して，初期値を与える式を指定します．
+例として，契約付き変数定義を用いて，初期値$3$を取る変数$x\in\mathbb{N}$を定義してみましょう．
+''<変数名>''が''x''，''<初期値式>''が''3''であることは自明なので，ここに指定すべき''<契約>''を考えます．
+ここで，変数$x$は自然数です．
+Racketにおいては，自然数かどうかを判定する述語''exact-positive-integer?''が定義されています．
+(ここでは「自然数」を，「正確な正整数」と考えます．$\mathbb{N}:=\left\{x\in\mathbb{Z}:x>0\right\}$．)
+述語は契約として使用することができるので，ここでは次の記述，
+<code scheme>
+(define/contract x exact-positive-integer? 3)
+</code>
+により，初期値$3$を取る変数$x\in\mathbb{N}$を定義することができます．
+契約付き変数定義により変数定義を行うと，
+これにより定義した変数に対して契約に反する操作を行った際に，**エラーとして検出する**ことができます．
+例えば，先の例における変数$x$に自然数$5$を代入することは，
+<code scheme>
+(set! x 5)
+</code>
+により行うことができます．
+しかし，自然数ではない値$0$を代入するような次のコード，
+<code scheme>
+(set! x 0)
+</code>
+を実行すると，エラーとなります．
+**契約付き関数定義**は，契約として定義域および値域を与えて(引数と戻り値に対して契約を与えて)，関数を定義します．
+以下に，契約付き関数定義の記述形式を載せます．
+<code scheme>
+(define/contract (<関数名> <引数>...)
+  (-> <引数契約>... <値域契約>)
+  <関数実装>)
+</code>
+この形式は，名前が''<関数名>''である関数を定義します．
+これにより定義される関数は，''<引数>...''にて定義した引数を受け取り，それらを用いて実装される式''<関数実装>''によって，その戻り値が計算されます．
+契約は，受け取る各引数，およびそれらによって計算される戻り値に対して，制約をかけることができます．
+各引数に対する契約は，定義した引数の個数分だけ順に，''<引数契約>''にて指定します．
+戻り値に対する契約は，一連の契約を表すリストの最後の項''<値域契約>''にて指定します．
+例として，実数$x,y\in\mathbb{R}$に対する，$\mathbb{R}$における通常の距離$d:\mathbb{R}^2\to\mathbb{R}_+$，
+$$ d(x,y):=\left|x-y\right| $$
+を，契約付き関数定義により定義してみましょう．
+<code scheme>
+(define/contract (d x y)
+  (-> real? real? (not/c negative?))
+  (abs (- x y)))
+</code>
+''real?''は実数かを判定する述語であり，''negative?''は負の実数かを判定する述語です．
+ここで，関数$d$は，2つの実数$x$，$y$を受け取り，非負実数を返す関数です．
+したがって，否定の契約を求める関数''not/c''を用いて，非負実数を保証する契約''(not/c negative?)''が構築できます．
+なお，関数''abs''は実数の絶対値を計算する関数であり，関数''-''は減算を計算する関数です．
+契約付き関数定義により定義された関数も，契約付き変数定義同様に，
+これにより定義した関数に対して契約に反する操作を行った際，または契約に反する結果が得られた際に，**エラーとして検出する**ことができます．
+具体的には，関数呼び出し時に与えられた引数に関して，契約に反する操作を検出することができます．
+また，計算結果としての関数の戻り値に関して，契約に反する結果を検出することができます．
+例えば，先に定義した関数$d$について，その引数に複素数を指定して呼び出すような式，
+<code scheme>
+(d 3+2i 5)
+</code>
+を実行すると，エラーとなります．
+==== 行列表現とその演算 ====
+Racketで行列およびその演算を表現する方法はいくつか考えられますが，ここでは線形代数関数群である''math/matrix''ライブラリを用いましょう．
+''math/matrix''ライブラリを利用するためには，次の式を評価すればよいです．
+<code scheme>
+(require math/matrix)
+</code>
+今回は，行列のなかでも，特にn行1列行列，いわゆる列ベクトルが使用できればよいです．
+一般の行列を生成する命令以外に，''math/matrix''ライブラリでは列ベクトルを生成する命令''col-matrix''が定義されています．
+例えば，列ベクトル$\left(1, 2, 3\right)^\mathrm{T}$を生成するためには，次の式を記述します．
+<code scheme>
+(col-matrix (1 2 3))
+</code>
+命令''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\|$は，
+<code scheme>
+(define/contract (d x y)
+  (-> col-matrix? col-matrix? (not/c negative?))
+  (matrix-norm (matrix- x y)))
+</code>
+により定義することができます．
+ここでは，列ベクトルかどうかを判定する述語としては''col-matrix?''を使用することができるので，これを用いて引数に契約を与えました．
+(ちなみに，一般の行列かどうかを判定する述語としては''matrix?''を使用することができます。)
+==== 遅延評価と無限列 ====
+===== Racket(Scheme) による実装 =====
+(書きかけ.)
+<code scheme>
+#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||"))
+</code>