myumori diary

Live like a cat

半順序を全順序に拡張するやつ

この記事は「みゅーもり Advent Calendar 2018」の11日目分となります:

adventar.org

記事を書く時間が十分に取れ無さそう(そもそも現時点で締切を2日超過している)なので, 計画を変更して簡単に書ける内容にしました.

(当初予定していた「banditとcoarse correlated equilibriumの関係性」の記事は(みゅーもりのやる気があれば)後日改めて書こうと思います……)

Jupyter Notebook: http://nbviewer.jupyter.org/github/myuuuuun/blog_contents/blob/master/notebooks/topological_sort.ipynb

定義

 X を任意の集合とし,  \succeq X 上で定義された二項関係とします.

このとき  \succeq が以下の3つを満たせば,  \succeq半順序(partial order) であるいいます

  1. 反射律:  \forall x \in X, \begin{align} x \succeq x \end{align}
  2. 推移律:  \forall x, y, z \in X, \begin{align} x \succeq y \land y \succeq z \implies x \succeq z \end{align}
  3. 歪対称律:  \forall x, y \in X, \begin{align} x \succeq y \land y \succeq x \implies x = y \end{align}

また  \succeq が以下の3つを満たせば,  \succeq全順序, または線形順序(total order / linear order) であるいいます

  1. 完全律:  \forall x, y \in X, \begin{align} x \succeq y \lor y \succeq x \end{align}
  2. 推移律:  \forall x, y, z \in X, \begin{align} x \succeq y \land y \succeq z \implies x \succeq z \end{align}
  3. 歪対称律:  \forall x, y \in X, \begin{align} x \succeq y \land y \succeq x \implies x = y \end{align}

 \succeq X 上の半順序である時,  (X, \succeq) 半順序集合(partially ordered set / poset) といいます. 同様に  \succeq X 上の全順序である時,  (X, \succeq) 全順序集合(totally ordered set) といいます.


半順序と全順序の違いは1.のみです. 完全律が成り立てば反射律も自動的に成り立つので, 全順序集合であれば必ず半順序集合になります.

イメージ

順序関係を有向グラフで表現してみます. たとえば  X = \{ a, b, c, d, e \} として,  X上の二項関係を \begin{align} \forall x, y \in X, x \succeq y \Longleftrightarrow \ \text{$x = y$, または$x$から$y$へ向かうpathがある} \end{align} と定義します.

f:id:myumori:20181214005655p:plain
順序の例

  • 例1は半順序ですが,  b \succeq c c \succeq b のいずれも成り立たないため, 全順序ではありません
  • 例2は全順序(したがって半順序)です. 全順序である時, 順序を崩さずに集合の各元を一列に並べることができます(例では a \to b \to c \to d \to e
  • 例3は b \succeq c \land c \succeq b ですが  b \neq c なので, 歪対称律を満たさず, 半順序ではありません
  • 例4は, 仮に推移律を満たすとすると  c \succeq d \land d \succeq e より  c \succeq e となります. しかし  e \succeq c でもあるので, 歪対称律が成り立ちません. したがって推移律と歪対称律が同時に成り立たないので, これは半順序ではありません

問題設定

半順序集合  (X, \succeq) が与えられた時, この半順序  \succeq を拡張して全順序  \succeq^{L} を構成することはできるでしょうか?

より厳密に書けば \begin{align} \forall x, y \in X, x \succeq y \implies x \succeq^{L} y \end{align} を満たす  X 上の全順序  \succeq^{L} は存在するでしょうか?

 X が有限の場合

 X の要素が有限個の場合は, 以下で具体的に説明するように, 半順序集合を非巡回有向グラフ(Directed Acyclic Graph, DAG)に置き換えてトポロジカルソートを実行することで, 条件を満たす全順序を具体的に構成できます.

Step 1: 半順序集合からDAGを構成する

まず半順序集合  (X, \succeq) を所与として, 次のような有向グラフ (V, E) を構成します( V は頂点の集合,  E が有向辺の集合です )

  • 頂点は  X の各元( V = X
  • 頂点  x \in V から 頂点  y \in V への有向辺  (x, y) は, \begin{align} (x, y) \in E \iff x \neq y \land x \succeq y \end{align}

このようにして作られた有向グラフ (V, E) には閉路(サイクル)はありません.

(∵)仮にサイクル  x_1 \to x_2 \to \ldots \to x_n \to x_1 があったとする. グラフの構成から \begin{align} x_1 \succeq x_2, \hspace{0.5em} x_2 \succeq x_3, \hspace{0.5em} \ldots, \hspace{0.5em} x_{n-1} \succeq x_n, \hspace{0.5em} x_n \succeq x_1 \end{align} が成り立つ.  \succeq は半順序なので, 推移律から  x_1 \succeq x_n が成り立つ. しかし  x_n \succeq x_1 でもあるので,  x_1 \neq x_n よりこれは歪対称律を破り,  \succeq が半順序であることに矛盾する.


 (V, E) に閉路が無いので, この有向グラフはDAGになっています. よって半順序集合からDAGを構成できました.

Step2: DAGにトポロジカルソートを適用する

トポロジカルソートとは, 以下のルールを満たすようDAGの各点を並び替えることを言います

どの点も, その点から有向パスが伸びている任意の点より前に並んでいる

下図の例でトポロジカルソートを実行すると,

  •  a, b, c, d, e
  •  a, c, b, d, e
  •  e, a, c, b, d
  •  a, e, c, b, d

などが解になります.

f:id:myumori:20181215054149j:plain
DAGの例

解はユニークではありませんが, 有向グラフに閉路が存在しなければ(DAGであれば)解は必ず1つ以上存在します. このことは次のアルゴリズム停止条件からわかります.

アルゴリズム

深さ優先探索を応用してトポロジカルソートを実現するアルゴリズムを紹介します. 以下の擬似コードwikipediaからのコピペです.

L ← トポロジカルソートされた結果の入る空の連結リスト

for each ノード n do
    if n に印が付いていない then
        visit(n)

function visit(Node n)
    if n に「一時的」の印が付いている then
        閉路があり DAG でないので中断
    
    else if n に印が付いていない then
        n に「一時的」の印を付ける
        for each n の出力辺 e とその先のノード m do
            visit(m)
        n に「恒久的」の印を付ける
        n を L の先頭に追加

補足: 深さ優先探索をする過程で「一時的」の印が付いたノードに辿り着いたということは, 探索をしてきたパスの中に閉路が存在するということです. したがってDAGであれば, visit() の最初のifの中身が実行されることはありません.

擬似コードPythonで書くと以下のようになります

class DiGraph(object):
    false_flag, temp_flag, perm_flag = 0, 1, 2
    
    def __init__(self, nodes, edges):
        self.nodes = nodes
        self.edges = edges
        self.node_size = len(nodes)

    def topological_sort(self):
        self.flags = {node:self.false_flag for node in self.nodes}
        self.sorted_nodes = []
        
        for node in self.nodes:
            if self.flags[node] == self.false_flag:
                self.visit(node)
        
        return self.sorted_nodes[::-1]
            
    def visit(self, node):
        if self.flags[node] == self.temp_flag:
            raise ValueError("DiGraph has a cycle")
            
        elif self.flags[node] == self.false_flag:
            self.flags[node] = self.temp_flag
            for node_to in self.edges[node]:
                self.visit(node_to)
                
            self.flags[node] = self.perm_flag
            self.sorted_nodes.append(node)

上の例を解いてみると

nodes = ["a", "b", "c", "d", "e"]
edges = {
    "a": ["b", "c"], 
    "b": ["d"], 
    "c": ["d"], 
    "d": [], 
    "e": [], 
}

digraph = DiGraph(nodes, edges)
digraph.topological_sort()
>> ['e', 'a', 'c', 'b', 'd']

となりました.

停止条件と計算量

上述のようにDAGであればアルゴリズムが途中でエラーを出力することは無いので,

  • 全ての辺を探索し( O(|E|)
  • 全ての頂点を sorted_nodes (L) に追加した( O(|V|)

時点でアルゴリズムは終了します. したがって計算量は  O(|V| + |E|) となります.

Step 3: トポロジカルソートの結果を全順序とする

トポロジカルソートの結果,  X の各元の順番が  x_1, x_2, \ldots, x_n と定まったら, あとは  X 上の二項関係  \succeq_L を \begin{align} x_1 \succeq_L x_2 \succeq_L \ldots \succeq_L x_n \end{align} で定義すれば, これは  \succeq を拡張した全順序になっています. なぜならば, 任意の  x, y \in X s.t.  x \succeq y に対し,

  • DAGの構成から,  (x, y) \in E
  • トポロジカルソートの定義から,  x y よりも前に並んでいる
  •  \succeq_L の構成から,  x \succeq_L y

となっているためです.

 X が無限の場合

 X の要素が無限個の場合はトポロジカルソートのアルゴリズムの停止が保証できないため, 上記のような構成的な証明は使えません.

代わりに有名なZorn補題を用いて証明します.

Zorn補題
 Sを半順序集合とする.  Sの任意の全順序部分集合が Sの中に上界を持つ時,  Sは少なくとも1つの極大元を持つ

いま X上で定義された2つの半順序  (\succeq_1) (\succeq_2) の包含関係  \subseteq を \begin{align} (\succeq_1) \subseteq (\succeq_2) \iff [ \forall x, y \in X, \ x \succeq_1 y \implies x \succeq_2 y ] \end{align} と定義します(つまり  \succeq_2 \succeq_1 を拡張したものだということです).

また与えられた半順序集合  (X, \succeq) に対し, 「  \succeq を拡張した X上の半順序」全体の集合を \begin{align} \mathcal{P} = \left\{ \succeq^{\prime} \mid \text{$\succeq^{\prime}$は$X$上の半順序} \land (\succeq) \subseteq (\succeq^{\prime}) \right\} \end{align} で定義します.

このとき,  (\mathcal{P}, \subseteq) は半順序集合になっています.

(∵)半順序集合が満たすべき性質1-3を順にチェックすれば良い.  (\succeq_1), (\succeq_2), (\succeq_3) \in \mathcal{P} x, y \in X を任意に取る.

  1. 反射律:  x \succeq_1 y \implies x \succeq_1 y が成り立つので,  (\succeq_1) \subseteq (\succeq_1)
  2. 推移律:  (\succeq_1) \subseteq (\succeq_2), \ (\succeq_2) \subseteq (\succeq_3) が成り立つとすると, \begin{align} x \succeq_1 y \implies x \succeq_2 y \implies x \succeq_3 y \end{align} が成り立つ. よって  (\succeq_1) \subseteq (\succeq_3)
  3. 歪対称律:  (\succeq_1) \subseteq (\succeq_2) \ \land \ (\succeq_2) \subseteq (\succeq_1) が成り立つとすると, \begin{align} x \succeq_1 y \iff x \succeq_2 y \end{align} が成り立つ. よって  (\succeq_1) = (\succeq_2)

また半順序集合 (\mathcal{P}, \subseteq)Zorn補題の仮定を満たします.

(∵) (\mathcal{P}, \subseteq) の(空でない)全順序部分集合  (\mathcal{D}, \subseteq) を任意に取る.

いま新たな X上の二項関係  \succeq_U を, \begin{align} \forall x, y \in X, \ x \succeq_U y \iff x \succeq_D y \hspace{0.5em} \text{for some $(\succeq_D) \in \mathcal{D}$} \end{align} で定義する. この時

  1.  (\succeq_U) は半順序
  2.  (\succeq_U) \in \mathcal{P}
  3.  (\succeq_U) (\mathcal{D}, \subseteq) の上界

が成り立つ. まず1.に関しては, 任意の  x, y, z \in X に対して

  • 反射律:  \forall (\succeq_D) \in \mathcal{D}, x \succeq_D x なので,  x \succeq_U x
  • 推移律:  x \succeq_U y かつ  y \succeq_U z だと仮定する. この時  (\succeq_{D1}), (\succeq_{D2}) \in \mathcal{D} が存在して \begin{align} x \succeq_{D1} y, \ y \succeq_{D2} z \end{align} が成り立つ. いま  (\mathcal{D}, \subseteq) は全順序集合だから, 完全律より  (\succeq_{D1}) \subseteq (\succeq_{D2}) または  (\succeq_{D2}) \subseteq (\succeq_{D1}) の少なくとも一方が成り立つ. 一般性を失わず  (\succeq_{D1}) \subseteq (\succeq_{D2}) が成り立つとすると, \begin{align} x \succeq_{D1} y \implies x \succeq_{D2} y \end{align} が成り立つ.  (\succeq_{D2}) が半順序であることから,  y \succeq_{D2} z と合わせて推移律より \begin{align} x \succeq_{D2} z \end{align} が成り立つ. よって x \succeq_{U} z
  • 歪対称律:  x \succeq_U y かつ  y \succeq_U x だと仮定する. この時  (\succeq_{D1}), (\succeq_{D2}) \in \mathcal{D} が存在して \begin{align} x \succeq_{D1} y, \ y \succeq_{D2} x \end{align} が成り立つ. よって推移律と全く同様の議論から,  i = 1 or  2に対して \begin{align} x \succeq_{Di} y \land y \succeq_{Di} x \end{align} が成り立つので,  \succeq_{Di} の歪対称律から  x = y となる

2.に関しては, すでに  (\succeq_U) が半順序であることがわかったので,  (\succeq) \subseteq (\succeq_U) が成り立つことだけを調べればよい. これは  \forall x, y \in X, \begin{align} x \succeq y \implies \forall (\succeq_{D}) \in \mathcal{D}, x \succeq_{D} y \implies x \succeq_{U} y \end{align} より成り立つ.

3.に関しては,  (\succeq_{U}) の構成から任意の  (\succeq_{D}) \in \mathcal{D} に対して, \begin{align} \forall x, y \in X, \ x \succeq_{D} y \implies x \succeq_{U} y \end{align} となり  (\succeq_D) \subseteq (\succeq_U) が成り立つので,  (\succeq_{U}) は確かに  \mathcal{D} の上界となる.


したがってZorn補題が適用でき, 半順序集合 (\mathcal{P}, \subseteq) は極大元を持ちます. 極大元の1つを  (\succeq_{M}) とすると, これは  (\succeq) を拡張した X上の全順序になっています.

(∵) (\succeq_{M}) \in \mathcal{P} なので,  (\succeq) を拡張した X上の半順序であることは明らか. したがって完全律だけを確かめれば良い.

任意の  x, y \in X をとる. 背理法の仮定として,  x \succeq_{M} y y \succeq_{M} x のいずれも成り立たないとする. この時  X 上の新しい二項関係  \succeq_{M^{\prime}} を \begin{align} \forall v, w \in X, \ v \succeq_{M^{\prime}} w \iff \begin{cases} v \succeq_{M} w, \ \text{ or } \\ v \succeq_{M} x \land y \succeq_{M} w \end{cases}\tag{1} \end{align} で定義する(つまり (\succeq_{M}) x \succeq_{M^{\prime}} y を加え, 推移律を満たすように拡張したもの).

この時  \succeq_{M^{\prime}} は半順序になっている. なぜならば任意の  u, v, w \in X に対し,

  • 反射律:  \succeq_{M} の反射律から  u \succeq_{M} u が成り立つので,  u \succeq_{M^{\prime}} u もOK
  • 推移律:  u \succeq_{M^{\prime}} v および  v \succeq_{M^{\prime}} w が成り立つと仮定する.
    •  (u, v) = (x, y) の場合:  x \succeq_{M} x かつ  y \succeq_{M} w なので, (1)より  x \succeq_{M^{\prime}} w
    •  (v, w) = (x, y) の場合:  u \succeq_{M} x かつ  y \succeq_{M} y なので, (1)より  u \succeq_{M^{\prime}} y
    • その他の場合:  (\succeq_{M}) が推移律を満たすことからOK
  • 歪対称律:  u \succeq_{M^{\prime}} v および  v \succeq_{M^{\prime}} u が成り立つと仮定する. 背理法の仮定から  (u, v) \neq (x, y), (y, x) なので,  (\succeq_{M}) が歪対称律を満たすことからOK

加えて  (\succeq) \subseteq (\succeq_{M}) \subseteq (\succeq_{M^\prime}) なので,  (\succeq_{M^\prime}) \in \mathcal{P}. これは  (\succeq_{M}) \mathcal{P} の極大元であることに反する.


したがって完全律が成り立ち,  \succeq_{M} \succeq を拡張した X上の全順序になっています. よって任意の半順序集合  (X, \succeq) に対して, それを拡張した全順序を構成できることがわかりました.

所感

  • トポロジカルソートのやり方は, 上で紹介したもの以外にもあります
  • Zorn補題選択公理と同値なので, 選択公理フリークの方は上記の命題が選択公理と同値なのか/真に弱い命題なのか気になるのだと思いますが, 詳しくないのでわかりません(すみません)
  • 整列可能定理(任意の集合は整列順序付け可能である)は選択公理と同値で, それよりも弱い「任意の集合は全順序で順序付け可能である」という命題は選択公理よりも真に弱いそうですが, 後者の命題よりもちょっと強い今回の命題はどうなんでしょうか

多項式の掛け算をFFTでやるやつを実装

この記事は「みゅーもり Advent Calendar 2018」の4日目分となります:

adventar.org

こんなにふざけたタイトルのAdvent Calendarが25日分すべて埋まるとは正直思っていませんでした. 登録してくださったみなさま, 本当にありがとうございます.

さて今日の内容は表題の通り, FFTを使って多項式同士の積の計算を高速化しようという話です.

専攻によっては学部の講義で扱う有名事実だそうですが, 無知なわたしはつい先日まで知りませんでしたし, これほど単純な問題に対して計算量を改善するテクニックがあるのかと驚いた話でもあったので, まとめておこうと思います.

Jupyter Notebook: http://nbviewer.jupyter.org/github/myuuuuun/blog_contents/blob/master/notebooks/polynomial_multiplication.ipynb

追記

  • (2018/12/5)添字のミスを修正しました. @stmfj さんありがとうございます!

問題設定

 f, g: \mathbb{R} \to \mathbb{R} n-1多項式

\begin{align} f(x) = a_0 + a_1 x + a_2 x^{2} + \ldots + a_{n - 1} x^{n - 1}, \\ g(x) = b_0 + b_1 x + b_2 x^{2} + \ldots + b_{n - 1} x^{n - 1} \end{align}

とする(  a_0, \ldots, a_{n - 1}, b_0, \ldots, b_{n - 1} \in \mathbb{R}*1. このとき, 2つの多項式の積  h = f \cdot g 2n-2多項式 \begin{align} h(x) = c_0 + c_1 x + c_2 x^{2} + \ldots + c_{2n-2} x^{2n-2} \end{align} となる( c_0, \ldots, c_{2n-2} \in \mathbb{R}). この  c_0, c_1, \ldots, c_{2n-2} を求め, 多項式  hを計算したい.

 f(x) = x^{2} + 2x + 3,  g(x) = 3x + 4 のとき,  h = f \cdot gは \begin{align} h(x) = 3x^{3} + 10 x^{2} + 17x + 12 \end{align} となる

素直な解法

単純に考えれば,  h k次の係数 c_kは \begin{align} c_k = \sum_{i=0}^{k} a_i b_{k-i} \end{align} で与えられるので, これを k=0, 1, \ldots, 2n-2について計算をすればよいでしょう. Pythonで実装すると

import numpy as np


def get_poly_product(f_coefs, g_coefs):
    """
    Parameters:
        f_coefs: array_like. fの係数のリスト. [0次の係数, 1次の係数, ..., 最高次の係数]
        g_coefs: array_like. gの係数のリスト. [0次の係数, 1次の係数, ..., 最高次の係数]
    Returns:
        h_coefs: ndarray. h=fgの係数のリスト. [0次の係数, 1次の係数, ..., 最高次の係数]
    """
    f_size, g_size = len(f_coefs), len(g_coefs)
    h_coefs = np.zeros([f_size+g_size-1], dtype=float)
    
    for i in range(f_size):
        for j in range(g_size):
            h_coefs[i+j] += f_coefs[i] * g_coefs[j]
            
    return h_coefs

のようになります. 上の例題を解いてみると

get_poly_product([3, 2, 1], [4, 3])

>> array([12., 17., 10.,  3.])

となり, 手で解いた結果と一致しました.

計算量

コードから明らかなように, この方法の計算量は  O(n^{2}) になります. さて, これよりも少ない計算量で  h の係数を求めることはできるでしょうか?

実は以下で説明するように, FFT高速フーリエ変換)を用いることによって計算量のオーダーを  O(n \log n) に改善することができます.

多項式補間を用いた解法

いま

相異なる k+1個の点  (x_0, y_0), (x_1, y_1), \ldots, (x_k, y_k) を通る k多項式は一意に存在する

という多項式の性質を用い, 次のような手順で hを計算することを考えます:

  • Step 1: 相異なる 2n-1個の点  x_0, x_1, \ldots, x_{2n-2} を任意にとり,  f(x_k), g(x_k) k = 0, 1, \ldots, 2n-2)を計算する
  • Step 2:  h(x_k) = f(x_k) \times g(x_k) k = 0, 1, \ldots, 2n-2)を計算する
  • Step 3:  h(x_0), h(x_1), \ldots, h(x_{2n-2}) を通る一意な 2n-2多項式が求める hなので, 点の座標から  c_0, c_1, \ldots, c_{2n-2} を逆算する

ここでVandermonde行列 \begin{align} V = \left(\begin{array}{ccccc} 1 & x_0 & x_0^{2} & \ldots & x_0^{2n-2} \\ 1 & x_1 & x_1^{2} & \ldots & x_1^{2n-2} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{2n-2} & x_{2n-2}^{2} & \ldots & x_{2n-2}^{2n-2} \end{array}\right) \end{align} ( 2n-1次正方行列)を導入します. また a_n = a_{n+1} = \ldots = a_{2n - 2} = b_n = b_{n+1} = \ldots = b_{2n - 2} = 0とし, \begin{align} \mathbf{a} = \left(\begin{array}{c} a_0 \\ a_1 \\ \vdots \\ a_{2n-2} \end{array}\right), \ \mathbf{b} = \left(\begin{array}{c} b_0 \\ b_1 \\ \vdots \\ b_{2n-2} \end{array}\right) \end{align} とします.

するとstep 1の計算は \begin{align} \left(\begin{array}{c} f(x_0) \\ f(x_1) \\ \vdots \\ f(x_{2n-2}) \end{array}\right) = V \ \mathbf{a}, \hspace{1em} \left(\begin{array}{c} g(x_0) \\ g(x_1) \\ \vdots \\ g(x_{2n-2}) \end{array}\right) = V \ \mathbf{b} \end{align} のような行列の掛け算に帰着できます.

また (2n - 1) \times (2n - 1) Vandermonde行列の行列式が \begin{align} \det(V) = \prod_{0 \leq s < t\leq 2n - 2} (x_t - x_s) \end{align} であるという有名事実から,  x_kが相異なれば  \det(V) \neq 0 となり,  V逆行列を持ちます. よってstep 3の計算も \begin{align} \left(\begin{array}{c} c_0 \\ c_1 \\ \vdots \\ c_{2n-2} \end{array}\right) = V^{-1} \ \left(\begin{array}{c} h(x_0) \\ h(x_1) \\ \vdots \\ h(x_{2n-2}) \end{array}\right) \end{align} という連立一次方程式の問題に帰着できます.

以上をまとめると, \begin{align} \left(\begin{array}{c} c_0 \\ c_1 \\ \vdots \\ c_{2n-2} \end{array}\right) = V^{-1} \ \left( V \ \mathbf{a} \circ V \ \mathbf{b} \right) \end{align} となります. ただし  \circ はベクトルの要素積です.

計算量

上の解法の計算量を各step毎に考えてみます.


前準備

まずstep 1, step 3の両方で使うVandermonde行列 Vを先に求めておきます. 各 x_kに対し,  x_k^{2}, x_k^{3}, \ldots, x_k^{2n-2}を求めるには 2n-3回の掛け算が必要になるので, 合計の計算回数は \begin{align} (2n - 1) \times (2n - 3) \end{align} 回程度になり, オーダーは O(n^{2})になります

Step 1

 f, gそれぞれに対して行列とベクトルの掛け算を行います.  Vの各行に対してそれぞれ 2n - 1スカラーの掛け算を行うので, 合計の計算回数は \begin{align} 2 \times (2n - 1) \times (2n - 1) \end{align} 回程度になり, オーダーは O(n^{2})になります

Step 2

 k = 0, 1, \ldots, 2n-2 に対して  h(x_k) = f(x_k) \times g(x_k) を求めればよいので, 計算回数は 2n-1回, オーダーは O(n)です

Step 3

 k変数の連立一次方程式をLU分解で解く場合の計算量は O(\frac{2}{3} k^{3})程度なので,  k = 2n-1である今回も計算量のオーダーは O(n^{3})になります


以上をまとめると, 多項式補間を用いた解法の計算量は \begin{align} O(n^{2}) + O(n^{2}) + O(n) + O(n^{3}) = O(n^{3}) \end{align} 程度となり, 素直な解法よりも遅くなってしまいました*2.

しかし x_0, \ldots, x_{2n-2}任意にとってよいという性質を活かすことで, 各stepでの計算量を大幅にへらすことができます.

多項式補間 + DFT(離散フーリエ変換)を使った解法

以降  f, g, h の定義域と値域を \mathbb{C}に拡張して考えます.

前準備の改良

まず標本点として  x_k = (\omega_{2n-1})^{k} k = 0, 1, \ldots, 2n-2)をとってきます. ただし \begin{align} \omega_t = e^{2 \pi i / t} \end{align} です( i虚数単位). このとき \begin{align} \omega_{t}^{k} = \exp \left(\frac{2 \pi k}{t} i \right) = \cos \left(\frac{2 \pi k}{t} \right) + i \sin \left(\frac{2 \pi k}{t} \right) \end{align} となります.

たとえば t=5のケースを複素平面上にプロットすると, 以下のようになります.

f:id:myumori:20181201202545j:plain
 \omega_{5}^{0}, \omega_{5}^{1}, \ldots, \omega_{5}^{4}複素平面上にプロット

いま \omega_tのべき乗は t周期になっている, すなわち \begin{align} \forall \lambda \in \mathbb{Z}, \ \omega_t^{\lambda t + k} = \omega_t^{k} \end{align} が成り立つので, 上記のように定めた  x_0, \ldots, x_{2n-2} に対するVandermonde行列は, 次のように簡単に求めることができます:

 \omega \equiv \omega_{2n-1} とすると, \begin{align} V = \left(\begin{array}{ccccc} 1 & \omega^{0 \times 1} & \omega^{0 \times 2} & \ldots & \omega^{0 \times (2n - 2)} \\ 1 & \omega^{1 \times 1} & \omega^{1 \times 2} & \ldots & \omega^{1 \times (2n - 2)} \\ 1 & \omega^{2 \times 1} & \omega^{2 \times 2} & \ldots & \omega^{2 \times (2n - 2)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{(2n - 2) \times 1} & \omega^{(2n - 2) \times 2} & \ldots & \omega^{(2n - 2) \times (2n - 2)} \end{array}\right) \\ = \left(\begin{array}{ccccc} 1 & \omega^{0} & \omega^{0} & \ldots & \omega^{0} \\ 1 & \omega^{1} & \omega^{2} & \ldots & \omega^{2n - 2} \\ 1 & \omega^{2} & \omega^{4} & \ldots & \omega^{2n - 3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{2n - 2} & \omega^{2n - 3} & \ldots & \omega^{0} \end{array}\right) \end{align} すなわち  V (j, k)-要素( j, k \in \{ 0, 1, \ldots, 2n-2 \})が \begin{align} V_{jk} = \omega^{k^{\prime}} \hspace{1em} \text{s.t.} \hspace{1em} k^{\prime} & \in \{ 0, 1, \ldots, 2n-2 \} \tag{1} \\ k^{\prime} & \equiv j \times k \text{ (mod $2n-1$)} \end{align} となります.

なおこのVandermonde行列  V をDFT行列といい,  V \mathbf{x} \in \mathbb{C}^{2n-1}に左から掛けることを離散フーリエ変換,  V^{-1} \mathbf{x}に左から掛けることを離散フーリエ逆変換といいます(詳細は後述).

 Vの要素に出てくる値は  \omega^{0}, \omega^{1}, \ldots, \omega^{2n - 2} だけなので, 各 x_kに対して陽にその2乗, 3乗, ...を計算する必要はなくなりました. しかしmodの計算をする必要があるため, 計算量自体は  O(n^{2}) のままです.

Step 3の改良

実は U \equiv \frac{1}{\sqrt{2n - 1}} V はユニタリ行列になっています. すなわち \begin{align} U^{*} U = U U^{*} = I \end{align} が成り立ちます( I単位行列で,  X^{*}は行列 Xの共役転置を表します).


(証明)

まず \omega^{k}複素共役 \bar{\omega^{k}}は, \begin{align} \bar{\omega^{k}} = \omega^{(2n - 1) - k} = \omega^{ - k} \end{align} となる(実数軸に関して線対称な位置への移動).

ここで行列 X j行目を切り出した行ベクトルを X(j, \cdot),  k列目を切り出した列ベクトルを X(\cdot, k)とする( j, k \in \{ 0, 1, \ldots, 2n-2 \}). このとき \begin{align} U^{*}(j, \cdot) \ U(\cdot, k) & = \frac{1}{2n - 1} \left(\begin{array}{ccccc} (\bar{\omega^{j}})^{0} & (\bar{\omega^{j}})^{1} &(\bar{\omega^{j}})^{2} & \ldots & (\bar{\omega^{j}})^{2n - 2} \end{array}\right) \left(\begin{array}{c} \omega^{0 \times k} \\ \omega^{1 \times k} \\ \omega^{2 \times k} \\ \vdots \\ \omega^{(2n - 2) \times k} \end{array}\right) \\ & = \frac{1}{2n - 1} \left(\begin{array}{ccccc} \omega^{- j \times 0} & \omega^{- j \times 1} & \omega^{- j \times 2} & \ldots & \omega^{- j \times (2n - 2)} \end{array}\right) \left(\begin{array}{c} \omega^{0 \times k} \\ \omega^{1 \times k} \\ \omega^{2 \times k} \\ \vdots \\ \omega^{(2n - 2) \times k} \end{array}\right) \\ & = \frac{1}{2n - 1} \sum_{t=0}^{2n - 2} \omega^{t \times (k - j)} \\ & = \begin{cases} \cfrac{1}{2n - 1} \cfrac{1 - (\omega^{k - j})^{2n - 1} }{1 - \omega^{k - j}} & \text{if $j \neq k$}\\ \cfrac{1}{2n - 1} (\sum_{t=0}^{2n - 2} 1) & \text{if $j = k$} \end{cases} \\ & = \begin{cases} 0 & \text{if $j \neq k$}\\ 1 & \text{if $j = k$} \end{cases} \end{align}  j \neq kのケースの最後の変形は,  (k - j) \times (2n - 1) \equiv 0 (mod  2n - 1)より, \begin{align} (\omega^{k - j})^{2n - 1} = \omega^{0} = 1 \end{align} となることを用いた.

 U(j, \cdot) \ U^{*}(\cdot, k) = Iに関しても, 全く同様にして示すことができる.


したがって \begin{align} V^{-1} = \frac{1}{\sqrt{2n - 1}} U^{-1} = \frac{1}{\sqrt{2n - 1}} U^{ * } = \frac{1}{2n - 1} V^{ * } \end{align} となり,  V逆行列を計算する処理は,  Vの共役転置を取る操作に置き換えられます. これを用いるとstep 3の計算量は

  •  V^{*}の計算に O((2n - 1)^{2}) = O(n^{2})
  •  V^{*} とベクトル  \left( h(x_0), \ldots, h(x_{2n - 2}) \right)^{\text{T}} の掛け算に O(n^{2})

かかるので, 合わせて  O(n^{2}) に改善します.

したがって前処理〜step 3までの全体の計算量のオーダーは \begin{align} O(n^{2}) + O(n^{2}) + O(n) + O(n^{2}) = O(n^{2}) \end{align} となります. ただこれでもまだ「単純な方法」と同じオーダーになってしまう(かつこちらの方が複雑な処理をする分時間がかかる)ので, さらに改良して O(n \log n)のオーダーを目指します.

余談: DFTについて

離散フーリエ変換の意味について少し補足をしておきます.  \mathbf{x} \in \mathbb{C}^{k}とし,  k \times kのDFT行列を F_kとします.  \mathbf{x}を離散フーリエ変換した結果が \mathbf{y} \in \mathbb{C}^{k}であるとすると, 上の議論から \begin{align} \mathbf{y} = F_k \mathbf{x} & \Longleftrightarrow F_k^{- 1} \mathbf{y} = I \mathbf{x} \\ & \Longleftrightarrow \cfrac{1}{k} F_k^{*} \mathbf{y} = I \mathbf{x} \\ & \Longleftrightarrow \cfrac{1}{k} \ \left(\begin{array}{ccccc} 1 & 1 & 1 & \ldots & 1 \\ \omega^{- 0} & \omega^{- 1} & \omega^{- 2} & \ldots & \omega^{- (k - 1)} \\ \omega^{- 0} & \omega^{- 2} & \omega^{- 4} & \ldots & \omega^{- (k - 1) \times 2} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \omega^{- 0} & \omega^{- (k - 1)} & \omega^{- (k - 1) \times 2} & \ldots & \omega^{- (k - 1) \times (k - 1)} \end{array}\right) \mathbf{y} = I \mathbf{x} \\ & \Longleftrightarrow y_0 \phi_0 + \ldots + y_{k - 1} \phi_{k - 1} = x_0 e_0 + \ldots x_{k - 1} e_{k - 1} \end{align} が成り立ちます. ただし \begin{align} \phi_j = \cfrac{1}{k} \ \left(\begin{array}{c} \omega^{- j \times 0} \\ \omega^{- j \times 1} \\ \omega^{- j \times 2} \\ \vdots \\ \omega^{- j \times (k - 1)} \end{array}\right) \end{align} であり,  e_j j番目の要素が1, それ以外が0のベクトルです.

したがって離散フーリエ変換は, 標準基底での座標 \mathbf{x}を,  \{ \phi_0, \phi_1, \ldots, \phi_{k - 1} \}を基底としたときの座標に, 座標変換する操作になります.

なお  \phi_s^{*} \phi_t = \mathbb{1} \{ s = t \} であることから,  \{ \phi_0, \phi_1, \ldots, \phi_{k - 1} \} \mathbb{C}^{k}の正規直交基底になります.

多項式補間 + FFTを使った解法

分割統治法を使って離散フーリエ変換を高速化しようというのが, FFT高速フーリエ変換)の考え方です.

いま  kは2のべき乗 であると仮定し,  \mathbf{x} \in \mathbb{C}^{k}フーリエ変換した結果である  \mathbf{y} = F_k \mathbf{x} を求めることを考えます( kが2のべき乗でない場合は,  kよりも大きい2のべき乗の中で一番小さいものを kと再定義し, それに合わせて \mathbf{x}に0を付け足せばよいです).

まず事実として, 2の倍数であるような tに対しては \begin{align} \omega_t^{2} = \exp \left( 2 \times \frac{2 \pi}{t} i \right) = \exp \left( \frac{2 \pi}{t / 2} i \right) = \omega_{t / 2} \tag{2} \end{align} が成り立ちます.

次に行列 F_kの偶数番号列だけを抜き出した行列を F_k^{\text{even}}, 奇数番号列だけを抜き出した行列を F_k^{\text{odd}} とします. すなわち \begin{align} F_k^{\text{even}} & = \left(\begin{array}{ccccc} 1 & \omega^{0 \times 2} & \omega^{0 \times 4} & \ldots & \omega^{0 \times (k - 2)} \\ 1 & \omega^{1 \times 2} & \omega^{1 \times 4} & \ldots & \omega^{1 \times (k - 2)} \\ 1 & \omega^{2 \times 2} & \omega^{2 \times 4} & \ldots & \omega^{2 \times (k - 2)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{(k - 1) \times 2} & \omega^{(k - 1) \times 4} & \ldots & \omega^{(k - 1) \times (k - 2)} \end{array}\right) \\ F_k^{\text{odd}} & = \left(\begin{array}{cccc} \omega^{0 \times 1} & \omega^{0 \times 3} & \ldots & \omega^{0 \times (k - 1)} \\ \omega^{1 \times 1} & \omega^{1 \times 3} & \ldots & \omega^{1 \times (k - 1)} \\ \omega^{2 \times 1} & \omega^{2 \times 3} & \ldots & \omega^{2 \times (k - 1)} \\ \vdots & \vdots & \ddots & \vdots \\ \omega^{(k - 1) \times 1} & \omega^{(k - 1) \times 3} & \ldots & \omega^{(k - 1) \times (k - 1)} \end{array}\right) \end{align} とし,  F_k^{\prime} = [F_k^{\text{even}}, \ F_k^{\text{odd}}] と定義します.

同様に \mathbf{x}の偶数番号行だけ抜き出したベクトルを \mathbf{x}^{\text{even}}, 奇数番号行だけ抜き出したベクトルをそれぞれ \mathbf{x}^{\text{odd}}とし, \begin{align} \mathbf{x}^{\prime} = \left(\begin{array}{c} \mathbf{x}^{\text{even}} \\ \mathbf{x}^{\text{odd}} \end{array}\right) \end{align} とします.

このとき定義より \begin{align} \mathbf{y} = F_k^{\prime} \mathbf{x}^{\prime} \end{align} が成り立ちます.  \mathbf{x}から \mathbf{x}^{\prime}を計算するコストは O(k)なので,  F_k^{\prime} \mathbf{x}^{\prime} の演算さえ O(k \log k)で実現できれば, 全体の計算量も O(k \log k)にできます.

いま F_k^{\text{even}}は \begin{align} F_k^{\text{even}} = \left(\begin{array}{ccc} 1 & \ldots & \omega_k^{0 \times (k - 2)} \\ \vdots & \ddots & \vdots \\ 1 & \ldots & \omega_k^{(k/2 - 1) \times (k - 2)} \\ \hline 1 & \ldots & \omega_k^{k/2 \times (k - 2)} \\ \vdots & \ddots & \vdots \\ 1 & \ldots & \omega_k^{(k - 1) \times (k - 2)} \end{array}\right) & = \left(\begin{array}{ccc} 1 & \ldots & \omega_{k/2}^{0 \times (k/2 - 1)} \\ \vdots & \ddots & \vdots \\ 1 & \ldots & \omega_{k/2}^{(k/2 - 1) \times (k/2 - 1)} \\ \hline 1 & \ldots & \omega_{k/2}^{k/2 \times (k/2 - 1)} \\ \vdots & \ddots & \vdots \\ 1 & \ldots & \omega_{k/2}^{(k - 1) \times (k/2 - 1)} \end{array}\right) \hspace{0.5em} \text{($\because$ (2))} \\ & = \left(\begin{array}{ccc} 1 & \ldots & \omega_{k/2}^{0 \times (k/2 - 1)} \\ \vdots & \ddots & \vdots \\ 1 & \ldots & \omega_{k/2}^{(k/2 - 1) \times (k/2 - 1)} \\ \hline 1 & \ldots & \omega_{k/2}^{0 \times (k/2 - 1)} \\ \vdots & \ddots & \vdots \\ 1 & \ldots & \omega_{k/2}^{(k/2 - 1) \times (k/2 - 1)} \end{array}\right) \hspace{0.5em} \text{($\because$ (1))} \\ & = \left(\begin{array}{c} F^{\prime}_{k/2} \\ F^{\prime}_{k/2} \end{array}\right) \end{align} と変形できます. 一方 D = \text{diag}(1, \omega_k, \omega_k^{2}, \ldots, \omega_k^{k/2 - 1})とすると,  F_k^{\text{odd}}も \begin{align} F_k^{\text{odd}} & = \left(\begin{array}{ccc} \omega_k^{0 \times 1} & \ldots & \omega_k^{0 \times (k - 1)} \\ \vdots & \ddots & \vdots \\ \omega_k^{(k/2 - 1) \times 1} & \ldots & \omega_k^{(k/2 - 1) \times (k - 1)} \\ \hline \omega_k^{k/2 \times 1} & \ldots & \omega_k^{k/2 \times (k - 1)} \\ \vdots & \ddots & \vdots \\ \omega_k^{(k - 1) \times 1} & \ldots & \omega_k^{(k - 1) \times (k - 1)} \end{array}\right) \\ & = \left(\begin{array}{ccc} \omega_k^{0 \times 1} & \ldots & \omega_k^{0 \times (k - 1)} \\ \vdots & \ddots & \vdots \\ \omega_k^{(k/2 - 1) \times 1} & \ldots & \omega_k^{(k/2 - 1) \times (k - 1)} \\ \hline \omega_k^{(k - k/2) \times 1} & \ldots & \omega_k^{(k - k/2) \times (k - 1)} \\ \vdots & \ddots & \vdots \\ \omega_k^{(k - (k - 1) ) \times 1} & \ldots & \omega_k^{(k - (k - 1) ) \times (k - 1)} \end{array}\right) \\ & = \left(\begin{array}{ccc} \omega_k^{0 \times 1} & \ldots & \omega_k^{0 \times (k - 1)} \\ \vdots & \ddots & \vdots \\ \omega_k^{(k/2 - 1) \times 1} & \ldots & \omega_k^{(k/2 - 1) \times (k - 1)} \\ \hline - \omega_k^{0 \times 1} & \ldots & - \omega_k^{0 \times (k - 1)} \\ \vdots & \ddots & \vdots \\ - \omega_k^{(k/2 - 1) \times 1} & \ldots & - \omega_k^{(k/2 - 1) \times (k - 1)} \end{array}\right) \\ & = \left(\begin{array}{c} D F^{\prime}_{k/2} \\ \hline - D F^{\prime}_{k/2} \end{array}\right) \end{align} と変形できます. 3番目の等式は, 複素単位円上で原点対称な位置への移動を行っていることから来ています.

まとめると \begin{align} F_k = [F_k^{\text{even}}, F_k^{\text{odd}}] = \left(\begin{array}{cc} F_{k/2} & D F_{k/2} \\ F_{k/2} & - D F_{k/2} \end{array}\right) \end{align} となるので, \begin{align} \left(\begin{array}{c} y_{0} \\ \vdots \\ y_{k/2 - 1} \end{array} \right) & = F_{k/2} \mathbf{x}^{\text{even}} + D F_{k/2} \mathbf{x}^{\text{odd}} \\ \left(\begin{array}{c} y_{k/2} \\ \vdots \\ y_{k - 1} \end{array} \right) & = F_{k/2} \mathbf{x}^{\text{even}} - D F_{k/2} \mathbf{x}^{\text{odd}} \end{align} が成り立ちます. よって k次元ベクトルにFFTを行う操作は,

  1.  \mathbf{x}を奇数列と偶数列に分割する操作( O(k)
  2.  k/2次元ベクトルに対しFFTを行う計算2つ(計算量は後述)
  3.  D  F_{k/2} \mathbf{x}^{\text{odd}} の計算( Dが対角行列なので O(k)で可能)
  4. ベクトルの和/差を計算する操作( O(k)

に分割できます

FFTの計算量

 \mathbf{x} k次元の時にFFTを行う計算量を T(k)とします. すると上の分割計算の手順より, 十分大きな kに対して \begin{align} T(k) \leq 2 T(k/2) + c k \end{align} という関係が成り立ちます( c > 0). また1点に対するFFTは元の値そのものなので,  T(1) = O(1)です. この時, \begin{align} T(k) = O(k \log k) \end{align} が成り立ちます(帰納法で証明できます).

高速フーリエ逆変換に関しても同様の分割統治法で計算でき, 計算量のオーダーは O(k \log k)です.

アルゴリズム全体の計算量

長々と説明してきましたが, 結局のところ hを求める手順としては

  • Step 1:  \mathbf{a} および  \mathbf{b}高速フーリエ変換する( O(n \log n)
  • Step 2: 変換したベクトル  V\mathbf{a}, V\mathbf{b} 同士の要素積  V\mathbf{a} \circ V\mathbf{b} を計算する( O(n)
  • Step 3:  V\mathbf{a} \circ V\mathbf{b} を高速フーリエ逆変換する( O(n \log n)

となります. したがってトータルの計算量は O(n \log n)となり, 単純な方法の場合よりも少ない計算量を実現することができました.

実装

これをPythonで実装すると, 以下のようになります.

def get_poly_product_fft(f_coefs, g_coefs):
    """
    Parameters:
        f_coefs: array_like. fの係数のリスト. [0次の係数, 1次の係数, ..., 最高次の係数]
        g_coefs: array_like. gの係数のリスト. [0次の係数, 1次の係数, ..., 最高次の係数]
    Returns:
        h_coefs: ndarray. h=fgの係数のリスト. [0次の係数, 1次の係数, ..., 最高次の係数]
    """
    h_size = len(f_coefs) + len(g_coefs) - 1
    h_fft_coefs = np.fft.fft(f_coefs, n=h_size) * np.fft.fft(g_coefs, n=h_size)
    h_coefs = np.fft.ifft(h_fft_coefs)
    
    return np.real(h_coefs)

例題を解いてみると,

get_poly_product_fft([3, 2, 1], [4, 3])
>> array([12., 17., 10.,  3.])

となり, 手で解いた結果と一致しました.

速度比較: 単純な方法 vs 多項式補間 + FFTを使った解法

私は「多項式補間 + FFT」で hを求める方法を初めて聞いたとき, 「いくら計算量のオーダーが小さくなるといっても, こんなに複雑なことをして計算が速くなるのか?」と疑問に思いました(はじめて聞いた人は普通そう思うでしょう).

というわけで簡単に実験をしてみました *3 *4:

f:id:myumori:20181203211534p:plain
speed comparison

結論: FFTしゅごい.

*1: f, gの次数が違う場合は, 大きい方の次数を n-1とし, 小さい次数の多項式の係数を適当に 0にすれば問題ありません

*2:またVandermonde行列は条件数が大きい行列として有名で, step 3の線形方程式をそのまま解くのは得策ではありません

*3: 本当はFFTも自分で実装するつもりでしたが, 面倒くさくなりました

*4:  2^{12}付近にkinkがあるのはどうしてなんでしょう

Mon, Mar 6, 2017 to Fri, Mar 10, 2017

平日は9:00-18:00まで作業.

なのだが, 実際のところ全員が揃うのは9時15分くらい. 帰りは17時45分くらいには帰る. そんな感じ.

RAに関して:

週の前半でレクチャーサイトの移行作業はおおよそ型がついた. 水曜日以降はrstファイルからコードブロックだけを抽出し, Jupyter Notebookに変換するパーサー(Sphinxのextension)を書いていた.

作業時間中ひたすらSphinxやらdocutilsのドキュメントを読んでいた. docutilsはドキュメントがきちんと揃っていないようなので, コードを読みつつ.

このextensionを拡張し, テキストブロックや数式などもパースできるようにした上で, QuantEconとは別口でオープンソースとして公開しようということになった. 大したものではないが, Private repoでの作業とは違い, 成果物として他の人に見せられるのでよかった.

キャンベラ情報:

  • キャンベラでは, おおよその店は5時には閉まってしまう. 例えばたいていのカフェの営業時間は7:00〜17:00という感じ. Canberra Centreというショッピングモールの中の店も, 22:00まで営業しているスーパー以外は17:00に閉まってしまう. ちょっとカルチャーショックだった.

そんな中, 中心部にいくつかあるコンビニは24時間やっている.

f:id:myumori:20170313220628j:plain
f:id:myumori:20170313220644j:plain

それなりに必要なものは揃っていて, コーヒーやドーナツも売っている. キャンベラにきて初めてどんどんドーナツした. してない.

  • キャンベラは中心部でも基本的に人が少ない. 日本の地方都市並みで, 例えば倉敷くらいのイメージ. For leaseの表示が出ている空きテナントがたくさんある.

土地がたっぷりあるのに, 高いビルを建てるのはなぜなんだろう. しかもどれもデザインが良かったりする.

f:id:myumori:20170313221345j:plain
f:id:myumori:20170313221401j:plain

  • 試しにフードコートにあったMee's sushiという寿司屋のパック寿司を食べてみた.

f:id:myumori:20170313221919j:plain

A$11ほど. サーモンは脂がのっていておいしかった. 残りのネタはいまいち. タコはパサパサ, いなり寿司は味が薄い. 食べられない人が多いからか, 山葵は入っていなかった(必要な場合は山葵の小袋をもらって付けて食べる).

  • 夕食を自炊しないときは, 主にアジア料理を食べに行っている. 安いので.

f:id:myumori:20170313222347j:plain
f:id:myumori:20170313222354j:plain

おいしい. どこの店もその国の人がやっているようだ.

中華料理屋に入った際, 別の店員が来る度に中国語で話しかけられた. 中華料理屋に1人で入るアジア人はおおよそ中国人だということなんだろうか.

Sun, Mar 5, 2017

12時起床.

教授 + イングランド銀行?から来たフランス人 + わたしたち2人 で野生のカンガルーを見に行った. ちょっと郊外に行くだけでカンガルーは見られるらしい. 一方でコアラは生息数が極めて少なく, オーストラリア人でも野生のコアラを見たことのある人は少ないようだ.

f:id:myumori:20170313173013j:plain

中心部から車で15分ほど. こんなちょっとした丘でもカンガルーがみられた.

f:id:myumori:20170313212847j:plain
f:id:myumori:20170313212907j:plain
f:id:myumori:20170313212909j:plain

カンガルーはそれほど臆病な動物では無いようで, 30mくらいの距離までは近づけた. しかし18-135mmのキットレンズだとこの大きさが限界.

家に返って改めて写真をLightroomでみてみると, ピントがぜんぜん合っていなくて悲しかった. 修業が必要.

その後, 教授の家へ招かれホームパーティーをした.

f:id:myumori:20170313213451j:plain

このビール, 今まで飲んだビールの中で一番フルーティーな香りがした.

帰りにApartmentまで車で送ってもらった時に驚いたのだが, OZは飲酒運転が許容されているらしい. およそ缶ビール1杯程度まではOKだそうだ.

Mon, Feb 27, 2017 to Sat, Mar 4, 2017

平日は9:00〜18:00まで労働. 大した内容がないのでまとめて書く.

2/27(月):

9時半に教授の研究室に行き, QuantEconの現状と今後の作業方針を聞いた.

quantecon.org

QuantEconのレクチャーサイトはSphinxというドキュメントビルダーを用いて作られている. SphinxはreST(reStructuredText)というマークアップ言語でドキュメントを書き, それをコンパイルしてHTMLやLaTeX等複数の形式に出力することができるツールらしい.

現在はPythonとJuliaのレクチャーで別々のドキュメントを用意しており, それゆえ2つのレクチャーサイトで内容が一部異なっている(Python版にはあるがJulia版にはない, というトピックがある). これを1つに統一し, コード等の言語によって異なる部分だけを別々のソースから読み込むよう改装することが, 主な作業目標だということだ.

その後実質監督となるポスドクの人と会い, 作業場所に移動.

f:id:myumori:20170313025302j:plain

新しくQuantEconのラボを作るとかで, 教室が2つ割り当てられていた. 新品のPCも10台弱?搬入されるということ. お金ありますね.

2/28(火)〜3/4(土):

当初ポスドクの人+RAの人1人+わたしたち2人の4人で作業をしていたのだが, 毎日RAやらプリドク?やらの人が増えていき, 計8人ほどで作業をすることに. マネーイズパワー.

わたしの作業内容は, レクチャーサイトの中に数式がinvalidだったりコードが間違っていたりする箇所がないかをチェックするプログラムを書く, という感じ. 特に大変なこともなければ勉強になることもなかった.

気を遣ってくれているのか, 昼食は毎日ラボのメンバーみんなでどこかに食べに行っている. 近くのKorean Restaurant or ANU内のフードコートが定番. 韓国料理でお米を食べているので, 特に日本食が恋しくはならなかった.

キャンベラの物価は(おそらくシドニーメルボルンもだが)東京と比べると1.3倍くらいだろうか. ミネラルウォーター, ペットボトル入りジュースなどは高い(500mlでA$3とか). 一方でチーズ, 果物などは日本より(はるかに)安い気がする. OZビーフを始めとした肉製品は, 残念ながら日本とそれほど変わらない値段.

f:id:myumori:20170313112353j:plain

このチェダーチーズなんて, 1個A$6ほど. そのため朝食には毎日こんなものを食べている.

f:id:myumori:20170313112510j:plain

オーストラリアは人件費が高いので(最低賃金がA$17/hとか?), レストランはかなり高い. 前述の韓国料理などアジア料理の店は安めで, ランチでA$12くらい. ハンバーガーも安め.

f:id:myumori:20170313111247j:plain

いつだかの夕食に食べに行った, Burger Heroという店のハンバーガー. 飲み物付きでA$15. $1=¥85なので, 1300円ほど. 東京でこのハンバーガーをこの値段では食べられないんじゃないだろうか.

他のレストラン, 特にディナーは高い. 夕食だと1人A$30以上はすると思う.

f:id:myumori:20170313113350j:plain

インドカレー. ナン + カレーでA$25くらい. ラッシーとチキンティッカを合わせてA$44だった(ただし, 1人分にはちょっと多すぎる量だった).

Sun, Feb 26, 2017

昼頃起床.

この日は特に予定が無かったので, National Gallery of Australia に行った.

Apartmentから4kmほどあったのでバスを使ったのだが, 仕組みが非常に分かりづらい. 系統がかなりたくさんある, 1本毎の時間もあいている, 平日と土日で運行する路線が違うなど, かなり初見殺しな印象.

1時間くらい検索して結局このサイトに現在地と目的地のバス停(付近の地名などではなく, 正確にバス停名を入れる必要がある)を打ち込み, google mapに飛ぶというソリューションに辿り着いた.

1日券を購入して乗車. 調べた限り$9.6だと書いてあったのだが, 運転手が$4.8だと言ってきたので, その通りの額を払った. 後からわかったのだが, これは学生料金のようだ. 他の地域でも種々の割引が大学生に適用されるようで, オーストラリアは学生に優しい.

15分ほどで, バーリー・グリフィン湖を挟んで対岸にある国会議事堂や最高裁判所などがある地域へ.

f:id:myumori:20170309224232j:plain
f:id:myumori:20170309224242j:plain

入場料は無料だった. この他の国立美術館や博物館も無料だそうで, これが文化政策の違いかと思った(が, 日本の美術館はいつも混んでいるので, 需要を抑えるためにもっとお金をとってもいいかもしれない).

それほど大きな美術館ではないが, オーストラリアの美術品に加え, 現代美術やアジア/アメリカ/ヨーロッパのセクションが設けられていた.

4時間ほどで, アジア, オーストラリアと現代美術を見て回った.

f:id:myumori:20170310070601j:plain
f:id:myumori:20170310070630j:plain
f:id:myumori:20170310070625j:plain

作品を撮影できるのも嬉しい.

f:id:myumori:20170310070620j:plain

この写真をみてびっくりした. こんなに歪んでるんだね……

この日も適当に何か料理をして食べました.

Sat, Feb 25, 2017

適当に寝たり起きたりした後にシドニーに到着. 9時間以上のフライトはかなり疲れる. 初めて機内にスリッパを持ち込んだが, とても快適だった.

f:id:myumori:20170308223229j:plain

f:id:myumori:20170308223248j:plain

機上から雲の様子を見て予想していたが, シドニーは雨. そのためか半袖では寒いほどの気温だった. 真夏のオーストラリアはどこへ.

入国審査場に行くと, このような機械が設置されていた.

f:id:myumori:20170309203703j:plain

ICチップのついたパスポートを持つ人の入国審査は自動化されているらしい. パスポートを機械に入れ犯罪歴の有無など3つ程の質問に答えるとチケットが発券され, あとはゲートで写真を撮影するだけで審査官に会う必要はない. おかげで5分もかからず入国できた.

乗り継ぎ便の手荷物検査場でPCを機内に置き忘れたことに気付いて慌てるも, サービスカウンターですぐに対応してもらい事なきを得た. 飛行機が1時間弱遅れていたことも幸いした.

f:id:myumori:20170309204534j:plain

プロペラ機でシドニーからキャンベラへ. 直線距離で350kmほどで, フライト時間は1時間弱.

f:id:myumori:20170309204651j:plain

キャンベラ周辺はこのような低木の林が広がっていた. 年間を通じて乾燥しているらしい.

空港でANUの教授と待ち合わせた後, 車でキャンベラ市街地を案内してもらった. 首都とはおもえないような小さな町.

その後宿泊するapartmentへ.

f:id:myumori:20170309205304j:plain


大学構内にこのような宿泊施設がいくつもある. すごい.

この日は適当に料理を作って寝た.