2018-12-15

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

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

記事を書く時間が十分に取れ無さそう（そもそも現時点で締切を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） であるといいます

反射律: $\forall x \in X,$ \begin{align} x \succeq x \end{align}
推移律: $\forall x, y, z \in X,$ \begin{align} x \succeq y \land y \succeq z \implies x \succeq z \end{align}
歪対称律: $\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） であるといいます

完全律: $\forall x, y \in X,$ \begin{align} x \succeq y \lor y \succeq x \end{align}
推移律: $\forall x, y, z \in X,$ \begin{align} x \succeq y \land y \succeq z \implies x \succeq z \end{align}
歪対称律: $\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$ を任意に取る.

反射律: $x \succeq_1 y \implies x \succeq_1 y$ が成り立つので, $(\succeq_1) \subseteq (\succeq_1)$
推移律: $(\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)$
歪対称律: $(\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} で定義する. この時

$(\succeq_U)$ は半順序
$(\succeq_U) \in \mathcal{P}$
$(\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, 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の補題は選択公理と同値なので, 選択公理フリークの方は上記の命題が選択公理と同値なのか/真に弱い命題なのか気になるのだと思いますが, 詳しくないのでわかりません（すみません）
整列可能定理（任意の集合は整列順序付け可能である）は選択公理と同値で, それよりも弱い「任意の集合は全順序で順序付け可能である」という命題は選択公理よりも真に弱いそうですが, 後者の命題よりもちょっと強い今回の命題はどうなんでしょうか

2018-12-04

多項式の掛け算を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_{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を行う操作は,

$\mathbf{x}$ を奇数列と偶数列に分割する操作（ $O(k)$ ）
$k/2$ 次元ベクトルに対しFFTを行う計算2つ（計算量は後述）
$D F_{k/2} \mathbf{x}^{\text{odd}}$ の計算（ $D$ が対角行列なので $O(k)$ で可能）
ベクトルの和/差を計算する操作（ $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があるのはどうしてなんでしょう

2017-03-13

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人で入るアジア人はおおよそ中国人だということなんだろうか.

2017-03-13

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だそうだ.

2017-03-13

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