FractalDB を作ってみよう(理論編)

はじめに

どうも. terasaki です. 機械学習の PoC を黙々やっています. テックブログを書く当番がやってきました. どうしよっかな? 何かこうかな? 行っている業務や扱っているデータの性質上, やってることをなかなか社外にオープンに出せないし 面白いこと書きたいよね? 面白いってなんだろう? 白い犬は尻尾も白いはずだから面白いよね! と考えてるうちに社内で書いた下書きがいっぱいできてしまったのでそろそろ真面目に書くとします.

今回のトピックは FractalDB と呼ばれる幾何学模様の人工データを生成するロジックを Julia で書いたというお話です. 結果として既存のコードよりも高速に生成できたよって話です.

何を作ったの?

IFS (反復関数系, Iterated Function System) による人工データを作っていました. 下記のようなフラクタル画像を生成するプログラムです.

こんな感じのやつ

近年,機械学習モデルの事前学習を行う際に自然画像の代わりに人工的に生成された画像を用いる試みを行う研究が行われています([1], [2], [3]). 面白そうだったので社内で行っている業務・実験でのダミーデータとして使えるのではと思って自分でコードを書いてました.

Julia ってなに?

Julia は LLVM に基づいた JIT コンパイル(Just-in-time compilation)で動作するプログラミング言語です. 関数の実行時に入力引数の型に応じてコンパイルされ機械語に変換され高速に動作することができます. また, Python のように高水準な記述もでき, 速度と生産性を両立することができます. 開発初期に公開された Julia の開発者のブログ Why We Created Julia では多くのプログラミング言語の良いところを取り込んだプログラミング言語を作りたかった様子がうかがえます.

IFS によるフラクタル画像構成法(大まかな方針)

さて, タイトルが「FractalDB を作ってみよう」 ということなのでコードも紹介すべきなのですが, 何を作るのかを書く必要があります. FractalDB は [1] で提案された比較的素直なアルゴリズムによって生成されるデータセットです. [2] で生成方法の改善版が提案され, 得られる画像はこちらの方が見栄えが良いのでこの方法を紹介することにします. 実際に動くコードは別の記事で紹介しようと思います.

Step1: IFS の構成

大まかな方針を説明します. はじめに正の整数 $N$ をランダムに選びます. $1\leq i \leq N$ に対してアフィン変換

\mathrm{Affine}(W_i, b_i): \mathbb{R}^2 \ni x \mapsto W_i x+b_i \in \mathbb{R}^2\quad i = 1,2,\dots,N

を用意します. ここで $W _ i$ は 2x2 の実行列, $b _ i$ は2次元ベクトルです.

次に $\sum _ {i=1} ^ N p _ i = 1$ を満たす正の実数 $p _ 1, \dots, p _ N$ を用意します. これは後述する点の生成アルゴリズムにおいて $i$ 番目のアフィン変換 $\mathrm{Affine}(W _ i, b _ i)$ の選択する確率値と解釈させます. 例えば $\mathrm{Affine}(W _ i, b _ i)$ を用意したもとでは $p _ i$ として $W _ i$ の行列式の絶対値 $| \det W _ i|$ に比例した値として定めることができます:

p_i = \frac{|\det W_i|}{\sum_{j=1}^N |\det W_j| }.

このようにしてカテゴリカル分布を $C = \mathrm{Cat}(p _ 1, \dots, p _ N)$ と置くことにします.

Step2: 点の生成

初期値 $x _ 0 \in \mathbb{R} ^ 2$ を用意します. 例えば $x _ 0 = (0, 0)$ とします. アフィン変換 $A _ 0$ を上述の $C$ の分布に従って選択します. つまり $\mathrm{Affine}(W _ i, b _ i)$ が選ばれる確率が $p _ i$ となるように選択します. $x _ 1\in\mathbb{R} ^ 2$ を次のようにして定義します:

x_1 = A_0(x_0)

再び $C$ の分布に従ってアフィン変換を選択します.それを $A _ 1$ とします. $x _ 2\in\mathbb{R} ^ 2$ を次のようにして定めます.

x_2 = A_1(x_1)

以下同様に $x _ {k+1} = A _ k(x _ k)$ のようにして定まる点列を $x _ 0, x _ 1, \dots$ を構成します. 例えば 100000 点分生成します.

Step3: 画像の描画

得られた点を特定のサイズ $(H, W)$ をもつ画像内に収まるようにします.例えば $H=W=368$ とします.

以上が大まかな作成方法です. 興味のある方は各自の得意なプログラミング言語で細部を補いながらトライしてみると良いでしょう.

Step 1 の具体的な方法

Step1 で述べたアフィン変換の生成の仕方で得られる描像が変化しますのでもう少し詳しく説明します.

成分を一様サンプルする方法

[1] の方法はアフィン変換 $\mathrm{Affine}(W, b)$ の $W$, $b$ の構成法として全ての成分を -1 から1 の間で値をランダムにサンプルする方法をとっています. この方法をとると点の生成のアルゴリズムを実施する際に座標値が発散するケースが出てきます. [1] に対応する公開実装では点が発散する時点で点の生成アルゴリズムを中断し画像の描画ステップに移行しています. こうなるとアルゴリズムが途中で中断してしまうので画像内にしめる図形の割合が小さいフラクタル画像が生成されてしまいます. この問題を解決するために画像内に占める図形の割合が一定以上のIFSを収集するパラメータサーチを実行しています. 公開されている実装は Python で書かれており,パラメータサーチを行うのに3時間ほどかかります.

アフィン変換の選び方によって点の生成アルゴリズム実行の際に座標が発散する理由を直感的に理解してみましょう. 話を単純にするために一次元の場合で考えてみます.$\lambda$ を正の実数 として $x _ {k+1} = \lambda x _ k$, $x _ 0 = 1$ で定まる漸化式を考えてみましょう. この場合 $x _ k$ は明示的に $x _ k = \lambda ^ k$ とかけます. $\lambda \gt 1$ ならば $x _ k$ の値は $k$ が大きくなるにつれて膨大に大きくなり値が発散します. 一方で $\lambda \lt 1$ の場合は有限の値に収まるので点の生成アルゴリズムを続けることができます.

このようにむやみに $\lambda$, 二次元の場合,行列の成分を選ぶのではなく適切な性質を持った行列を生成する必要があります.

行列の特異値・固有値に基づく方法

[2] の研究は上述の点の発散する問題を解決しています. 反復関数系の理論において系を構成するアフィン変換 $\mathrm{Affine}(W, b)$ の $W$ は縮小写像である必要があります. すなわち:

ある $s \lt 1$ が存在して任意の $x, y \in \mathbb{R} ^ 2$ に対して$||Wx - Wy || \leq s || x -y ||$ を満たす

を要請します. ここで $||x||$ ベクトルは $x\in \mathbb{R} ^ 2$ のノルムを表しています.

つまり下記で定める $W$ のスペクトルノルム ||W|| が 1 より小さいことを要請します.

||W|| \underset{\mathrm{def}}{=} \underset{z\neq 0}{\sup} \frac{||Wz||}{||z||} = \underset{||e||=1}{\sup} ||We||

線形代数の一般論よりスペクトルノルムは行列 $W$ の特異値の最大値に一致します.

[2] ではこの理論的な要請を満足するように $W$ を構成しています. 特異値分解では $W$ を $W = U \Sigma V ^ \top$ のように直交行列 $U, V$, 特異値を格納する対角行列 $\Sigma$ によって分解します. よって右辺の式に基づいて $W$ を構成すれば所望の行列を生成することができます. 直交行列は回転行列

R(\theta) = \begin{bmatrix}\cos \theta & -\sin \theta \\\\\sin \theta & \cos \theta\end{bmatrix}

または 行列式が -1 となる適当な直交行列 $D$ と $R(\theta)$ の積として得ることができます. 例えば $D$ として

D = \begin{bmatrix}1 & 0 \\\\0 & -1\end{bmatrix}

があります. $d _ 1, d _ 2$ をランダムに -1 または 1 と選択し, $0\leq \sigma _ 2 \leq \sigma _ 1 \lt 1$ となる $\sigma _ 1, \sigma _ 2$ を生成すれば

W \underset{\mathrm{def}}{=} R(\theta) \begin{bmatrix}\sigma_1 & 0 \\\\0 & \sigma_2\end{bmatrix} R(\phi) \begin{bmatrix}d_1 & 0 \\\\0 & d_2\end{bmatrix}

によって最大特異値が $\sigma _ 1$ となる $W$ を構成することができます.

特異値の制御

上記の説明では特異値として $0\leq \sigma _ 2 \leq \sigma _ 1 \lt 1$ を満たすようにしました. [2] の研究では良い形状のフラクタルを得るためにIFSを構成するアフィン変換の行列の特異値にも一定の制約を設けることを提案しています.

$\sigma$-factor のお話

[2] の説明に従って $\sigma$-factor という量を定義します.

整数 $N$ と $0\leq \sigma _ {i,2} \leq \sigma _ {i, 1} \leq 1\ (1\leq i\leq N)$ を満たす列

\Sigma = \{ \sigma_{1,1}, \sigma_{1,2}, \sigma_{2,1}, \sigma_{2,2},\dots, \sigma_{N,1}, \sigma_{N,2} \}

に対して $\Sigma$ の $\sigma$-factor という量 $\alpha$ を次で定義する:

\alpha = \alpha(N, \Sigma) \underset{\mathrm{def}}{=} \sum_{i=1}^N (\sigma_{i, 1} + 2\sigma_{i, 2}).

このとき $\sigma _ {i,1}\leq 1$, $\sigma _ {i, 2} \leq 1$ が成立するので

0\leq \alpha \leq 3N

が成立する.

[2] での実験的な調査では $N = 2,3,4$ の場合 $0.5(5 + N) \leq \alpha(N, \Sigma) \leq 0.5(6 + N)$ を満たすように

\Sigma = \{\, \sigma_{1,1}, \sigma_{1,2}, \sigma_{2,1}, \sigma_{2,2},\dots, \sigma_{N,1}, \sigma_{N,2} \mid 0\leq \sigma_{i,2} \leq \sigma_{i, 1} \leq 1, 1\leq i \leq N\,\}

を構成しそれらを $\sigma _ {i,1} \leq \sigma _ {i, 2}$ を特異値とする $W _ i$ を構成するとバラエティのあるフラクタル画像が得られると報告しています.

具体的な構成法

構成法を論文よりも丁寧に追っていきましょう. 数学的な説明をするので「です,ます」調を抜かします.

与えられた整数 $N$ と $0\leq \alpha \leq 3N$ を満たす $\alpha$ に対して $\sigma$-factor が $\alpha$ である列

\Sigma = \{\, \sigma_{1,1}, \sigma_{1,2}, \sigma_{2,1}, \sigma_{2,2},\dots, \sigma_{N,1}, \sigma_{N,2} \mid 0\leq \sigma_{i,2} \leq \sigma_{i, 1} \leq 1, 1\leq i \leq N\,\}

を構成することを考える. この制約問題を解くことを 「制約問題 $C(N, \alpha)$ 解く」と呼ぶことにする. この問題を解くために $\sigma _ {1,1}$ がとりうる範囲を考える. 簡単な不等式評価から

\begin{aligned}\alpha&= \sigma_{1, 1} + 2\sigma_{1, 2} + \sum_{i=2}^N(\sigma_{1, 1} + 2\sigma_{1, 2}) \\&\leq \sigma_{1, 1} + 2\sigma_{1, 2} + 3(N-1) \quad \because \sigma_{i,j} \leq 1 \\&\leq \sigma_{1, 1} + 2\sigma_{1, 1} + 3(N-1) \quad \because \sigma_{1,2}\leq \sigma_{1,1} \\&\leq 3\sigma_{1, 1} + 3(N-1)\end{aligned}

一方で

\begin{aligned}\alpha&= \sigma_{1, 1} + 2\sigma_{1, 2} + \sum_{i=2}^N(\sigma_{1, 1} + 2\sigma_{1, 2}) \geq \sigma_{1,1} \quad \because \sigma_{i,j} \geq 0\end{aligned}

が従う. よって $(\alpha - 3(N-1))/3 \leq \sigma _ {1,1} \leq \alpha$ を満たす. $0\leq \sigma _ {1,1}\leq 1$ という条件と合わせると

\max\left(0, \frac{1}{3}(\alpha - 3(N-1))\right) \leq \sigma_{1,1} \leq \min(1, \alpha)

と書くことができる. いま, 上記の条件を満たすように $\sigma _ {1,1}$ を決めたとする. 例えば上記の範囲から乱数を生成する方法がある. $\sigma _ {1,1}$ が固定されたもとで $\sigma _ {1,2}$ のとりうる範囲を調べたい.

簡単な不等式評価によって

\begin{aligned}\alpha&= \sigma_{1, 1} + 2\sigma_{1, 2} + \sum_{i=2}^N(\sigma_{1, 1} + 2\sigma_{1, 2}) \\&\leq \sigma_{1, 1} + 2\sigma_{1, 2} + 3(N-1)\end{aligned}

そして

\begin{aligned}\alpha&= \sigma_{1, 1} + 2\sigma_{1, 2} + \sum_{i=2}^N(\sigma_{1, 1} + 2\sigma_{1, 2}) \geq \sigma_{1,1} + 2\sigma_{1,2}\end{aligned}

が従う. よって

\frac{1}{2}(\alpha - 3(N-1) - \sigma_{1,1}) \leq \sigma_{1,2} \leq \frac{1}{2}(\alpha - \sigma_{1,1})

を得る. さらに $0\leq \sigma _ {1,2} \leq \sigma _ {1,1}$ という制約と合わせると次を得る:

\max\left(0, \frac{1}{2}(\alpha - 3(N-1) - \sigma_{1,1})\right) \leq \sigma_{1,2} \leq \min\left(\sigma_{1,1}, \frac{1}{2}(\alpha - \sigma_{1,1})\right)

いま, 上記の条件を満たすように $\sigma _ {2,2}$ を決めたとする. このとき $\sigma _ {i, j}\ (i\geq 2)$ らを決める方法は制約問題 $C(N-1, \alpha-\sigma _ {1,1} - 2\sigma _ {1,2})$ を解くことに帰着される. このとき帰納的に解いてくと $N=1$ に関する問題 $C(1, \alpha)$ を解くことになる. 次の不等式

\alpha = \sigma_{1,1} + 2 \sigma_{1, 2} \geq \sigma_{1,2} + 2 \sigma_{1, 2} = 3\sigma_{1, 2}

および

\alpha = \sigma_{1,1} + 2 \sigma_{1, 2} \leq 1 + 2 \sigma_{1, 2}
\max\left(0, \frac{1}{2}(\alpha - 1)\right) \leq \sigma_{1,2} \leq \min\left(1, \frac{\alpha}{3}\right)

がえられる. この条件を満たすように $\sigma _ {1,2}$ を定めた後 $\sigma _ {1,1} = \alpha - 2\sigma _ {1, 2}$ と定めれば良い.

生成結果

[2] の方法に従った方法で IFS を構成し [1] のプロトコルに合わせてグレー調の画像を生成しました.

いくつかの例.

特徴的な幾何学模様が得られている様子がわかります. 機械学習の話は抜きにしても面白いトピックだと感じました. また人工的に作られた画像ですので気軽にクラウド環境で画像をアップロードしたり業務でクラウド環境の機能を触る練習のデータとして使うといった運用もできます. 初めて使う技術などで取り合えず画像を入力させて遊んでみようといった時に使えるのでとても良いです. アイコンに使うというのもアリですね.

Julia で書いて良かったこと

まず速い

  • Julia は JIT コンパイル方式で動作するプログラミング言語のため高速に動作します. [1], [3] に対応する公開されている実装は Python で書かれているためとても時間がかかります. 例えば [1] におけるパラメータサーチや実際の画像生成も含めると数時間以上待たないといけません. Python のコードを動かしている間に Julia で実装ができたのでそちらを動かすことにしました. 上記の画像は Julia で書いたコードから得られたものになります. vCPU が 16を持つ GCP インスタンスで並列処理することで 15 分もかからずに生成することができました. かなり速いと思います.[2] に対応する実装は一部を Numba で高速化し洗練された実装になっていますが, メモリ上で画像を生成する部分と比較すると手元の実装の方が高速に動くことがわかりました. 今回は [1] で提案されていた FractalDB を改善手法 [2] によって生成しましたが, [3] で提案されている RCDB(Radial Contour DataBase) も高速に生成できます.

書きやすい

  • Julia は単に高速に動くだけでなく Ruby, Python のような動的な言語であり, 高水準な記述も可能です. また Unicode 文字列の入力をサポートしているので $\sigma _ {i, 1}$ という数式に対応する変数として σᵢ₁ を使えるようにアルゴリズムと実装の乖離を抑えることができます. 今回のトピックは数式ベースでデータを行うアルゴリズムなので Python を使わずに Julia を使うことは技術選択の観点から合理的だと思いました. 本来 Julia で書いたならそのコードも併せて公開するのが筋ですが, 多くの人に読んでもらおうと前提知識を埋めようとすると文章が膨れてしまったので別の記事にしてまとめようと思います.

Python との相互連携

  • Python の良いところはライブラリが豊富なところです. Julia は後発の言語であり新しい分野に手を出すと細かい部分で「Python にあった便利なアレ」に相当するライブラリがなくて詰まるという現象に出会います. が,PyCall.jl や PythonCall.jl を使うと pyimport("ホニャララ") とすることで大抵の Python ライブラリを呼び出すことができますので「ひとまずは Python のライブラリで凌いで後で自前で書き直す, Python 側の実装をリファレンス実装として用いてテストを書く」といった運用もできます. Julia らしく書けるように IJulia.jl, PyPlot.jl, Pandas.jl, PyPlotly.jl, SciPy.jl, ScikitLearn.jl, Seaborn.jl などメジャーなライブラリは有志によって整備がされています. 逆に Julia のライブラリも Python から呼び出す仕組みも最近では快適にできるようになったので Python を普段使っているメンバーに負い目を感じることなく速度が必要な場面で Julia を書くことができます.

まとめ

  • FractalDB というものを紹介しました.
  • アルゴリズムの解説を書きました.
  • お気持ちしか書いてませんが, Julia で書くメリットや Julia の近況を紹介しました. せっかくなので今後は Julia の話ができればと思います.

References

  • [1] Kataoka H, Okayasu K, Matsumoto A, Yamagata E, Yamada R, Inoue N, Nakamura A, Satoh Y. Pre-training without natural images. InProceedings of the Asian Conference on Computer Vision 2020.
  • [2]Anderson C, Farrell R. Improving Fractal Pre-training. InProceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision 2022 (pp. 1300-1309).
  • [3] Kataoka H, Hayamizu R, Yamada R, Nakashima K, Takashima S, Zhang X, Martinez-Noriega EJ, Inoue N, Yokota R. Replacing Labeled Real-Image Datasets With Auto-Generated Contours. InProceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition 2022 (pp. 21232-21241).