z軸周りに回転された2つの姿勢のクォータニオンのSlerpと角度のLerpの関係性について

問題

z軸まわりの回転角 \theta_0 と \theta_1 によって得られる2つの姿勢を考える。この2つの姿勢を内分比 tで球面線形補間(Slerp)によって内分した姿勢は、角度を線形補間(Lerp)した姿勢と一致するか？

考察

2つの姿勢をクォータニオンで表すとそれぞれ次のようになる。

\left(\cos\frac{\theta_0}{2}, 0, 0, \sin \frac{\theta_0}{2}\right),\quad \left(\cos\frac{\theta_1}{2}, 0, 0, \sin \frac{\theta_1}{2}\right)

この2つの姿勢の間を t（ 0 \le t \le 1）の割合で内分する姿勢を考える。まず、Slerpで補間した姿勢の w 成分は次のように表される。

\frac{\sin (1-t)\frac{\theta_1-\theta_0}{2}}{\sin\frac{\theta_1-\theta_0}{2}} \cos\frac{\theta_0}{2} + \frac{\sin t \frac{\theta_1 - \theta_0}{2}}{\sin\frac{\theta_1-\theta_0}{2}}\cos \frac{\theta_1}{2}

次に、ヨー角を線形補間した姿勢の w 成分は以下のように表される。

\cos\frac{(1-t)\theta_0+t\theta_1}{2}

次のように定義する。

\begin{aligned} A &\equiv \left[ \sin (1-t)\frac{\theta_1-\theta_0}{2} \right] \cos\frac{\theta_0}{2} + \sin \left[ t \frac{\theta_1 - \theta_0}{2}\right] \cos \frac{\theta_1}{2} \\ B &\equiv \sin\frac{\theta_1 - \theta_0}{2} \cos \frac{(1-t)\theta_0 + t\theta_1}{2} \end{aligned}

問題は A = B が成り立つかどうかに帰着される。

簡単化のために、以下の定義を行う

\begin{aligned} c_i &\equiv \cos \frac{\theta_i}{2}, \quad s_i \equiv \sin \frac{\theta_i}{2}, \\ \hat{c}_i &\equiv \cos \frac{t\theta_i}{2}, \quad \hat{s}_i \equiv \sin \frac{t\theta_i}{2} \end{aligned}

この定義のもとで A と B を展開すると、次のようになる。

\begin{aligned} A &= (c_0^2 s_1 - c_0 c_1 s_0) \hat{c}_0 \hat{c}_1 + (c_1 s_0^2 - c_0 s_0 s_1) \hat{c}_0 \hat{s}_1 \\ &\quad + (c_0 s_0 s_1 - c_1 s_0^2) \hat{c}_1 \hat{s}_0 + (c_0^2 s_1 - c_0 c_1 s_0) \hat{s}_0 \hat{s}_1 \\ B &= (c_0^2 s_1 - c_0 c_1 s_0) \hat{c}_0 \hat{c}_1 + (c_1 s_0^2 - c_0 s_0 s_1) \hat{c}_0 \hat{s}_1 \\ &\quad + (c_0 s_0 s_1 - c_1 s_0^2) \hat{c}_1 \hat{s}_0 + (c_0^2 s_1 - c_0 c_1 s_0) \hat{s}_0 \hat{s}_1 \end{aligned}

よって A = B が成り立つので、クォータにオンの w成分が一致することが示された。

次に、z成分に関しても同様の確認を行う。

\begin{aligned} C &\equiv \left[ \sin (1-t)\frac{\theta_1-\theta_0}{2}\right] \sin\frac{\theta_0}{2} + \sin \left[ t \frac{\theta_1 - \theta_0}{2}\right] \sin \frac{\theta_1}{2} \\ D &\equiv \sin\frac{\theta_1 - \theta_0}{2} \sin \frac{(1-t)\theta_0 + t\theta_1}{2} \end{aligned}

この C および D を展開すると次のようになる。

\begin{aligned} C &= (c_0 s_0 s_1 - c_1 s_0^2) \hat{c}_0 \hat{c}_1 + (c_0 s_0 s_1 - c_1 s_0^2) \hat{s}_0 \hat{s}_1 \\ &\quad + (c_0^2 s_1 - c_0 c_1 s_0) \hat{c}_0 \hat{s}_1 + (c_0 c_1 s_0 - c_0^2 s_1) \hat{c}_1 \hat{s}_0 \\ D &= (c_0 s_0 s_1 - c_1 s_0^2) \hat{c}_0 \hat{c}_1 + (c_0 s_0 s_1 - c_1 s_0^2) \hat{s}_0 \hat{s}_1 \\ &\quad + (c_0^2 s_1 - c_0 c_1 s_0)\hat{c}_0 \hat{s}_1 + (c_0 c_1 s_0 - c_0^2 s_1) \hat{c}_1 \hat{s}_0 \end{aligned}

よって C = D も成り立ち、クォータニオンのz成分も一致することが示された。

結論

z軸まわりの回転においては、クォータニオンのSlerpによって得られる姿勢と、ヨー角を線形補間（Lerp）した姿勢は一致する。すなわち、この特定の軸回りの回転において、SlerpとLerpは同一の補間結果を与える。

SegmentTreeはモノイドを扱えるデータ構造という話について

TL;DR

SegmentTreeはモノイドが扱えるデータ構造だという理解をしていたが、モノイドではなく半群を扱えるデータ構造だと理解した方が良いのではないかと考えている。

前提知識

半群とモノイドについて

まず半群とモノイドについて整理する。

半群

結合法則を満たす群構造のことを半群と呼ぶ。結合法則とは、任意の群の要素 a、 b、 cについて

\begin{eqnarray} (a \circ b) \circ c = a \circ (b \circ c) \end{eqnarray}

が成立することを言う。

モノイド

結合法則を満たしかつ単位元が存在する群構造のことをモノイドと呼ぶ。ある群の要素 eとは単位元であるとは、任意の群の要素 aについて

\begin{eqnarray} a \circ e = e \circ a = a \end{eqnarray}

が成立することを言う。

すなわち、モノイドと半群の違いは単位元が存在するかどうかである。

SegmentTreeについて

AC libraryによると

モノイド 、つまり
- 結合律
- 単位元の存在
を満たす代数構造に対し使用できるデータ構造です。
長さNのSの配列に対し、
- 要素の1点変更
- 区間の要素の総積の取得
をO(logN) で行うことが出来ます。

と説明されている。コンストラクタは2種類あり、 a. 群の集合、要素間の演算、単位元、SegmentTreeサイズ b. 群の集合、要素間の演算、単位元、初期配列

の2つの与え方がある。

背景

ABC254-F問題をやっていた時に、正の整数列のある区間の最大公約数を高速で求めるためにSegmentTreeが使えないか検討していた。そのとき、SegmentTreeを使う条件として結合則の成立と単位元の存在が要求されるが、結合則の成立は問題なさそうだが、単位元については少し手が止まった。数学的には、あらゆる正の整数を掛け合わせたものすごい大きい整数をが単位元になりそうだが、実装する上ではoverflowが起きてしまうのでそのような単位元は実装上は使えない。そこで、0をそのものすごく大きい数の代用に使おうと考えていたが、ここでふと「あれ、単位元の定義がなぜ必要なのか」と思い始めた。0を単位元として定義するにしても、正の整数列の区間取得を行う限りにおいては0が登場することはなく、なぜ「0を単位元にします」という宣言を実装しないといけないのか疑問に感じ、SegmentTreeのデータ構造の一般性を失わないために、「SegmentTreeは半群に対して使えるデータ構造」という理解をした方がよいのではないかと考えた。

SegmentTreeに単位元を与える必要があるか

SegmentTreeの実装において単位元は次の2ケースで使われるようである。

  1. Treeノードの初期値

  2. 区間取得を行うときの区間値の初期設定値

1.について、ノードをどの値で初期化することに関してはユーザーの責務とする考え方も自然に感じる。単位元で初期化することにより部分的にノードの値が更新されたときにも区間取得クエリが正しく動くため便利という考え方もあるかもしれないが、部分的にノードの値を更新するような操作を許さず、初期化されていないノードを含む区間取得クエリに対しては返す値を保証しない（未定義動作）ということでも良いのではないか。

2.について、これは完全に実装をシンプルにするためだけの都合なので、それを理由にユーザーに単位元を要求してほしくない。

以上から、一般的に「SegmentTreeは半群を扱えるデータ構造」という理解をしたい。ただし、あるライブラリが実装の都合で単位元を要求しその結果として「このSegmentTreeはモノイドを扱えるデータ構造」というように使い方を限定すること自体は問題ないと思っている。

C++クラスのpublic、protected、privateの挙動確認

派生クラスから基底クラスの要素にアクセスできるか

基底クラスの指定 \ 派生クラスの継承の仕方 public protected private
public ◯ ◯ ◯
protected ◯ ◯ ◯
private × × ×

クラス外から基底クラスの要素にアクセスできるか

Baseクラスの指定 \ 派生クラスの継承の仕方 public protected private
public ◯ × ×
protected × × ×
private × × ×

派生クラスで基底クラスの要素がどう扱われるか

Baseクラスの指定 \ 派生クラスの継承の仕方 public protected private
public public protected private
protected protected protected private
private no access no access no access

まとめ

クラス外からはpublicのみアクセス可能。派生クラスからは（規定クラスをどのような指定子で継承したかによらず）publicとprotectedがアクセス可能。
派生クラスで規定クラスの要素がどう扱われるかは、「基底クラスの指定子」と「継承指定子」の厳しい方となる。

きたまさ法

概要

数列 {a}について、隣接 k項間の漸化式および初期条件（第 k項までの数列の値）が与えられた時に a_{n}を高速に求める方法。言い換えると、 a_{n}を a_{0}, a_{2}, \dots a_{k-1}で表現する方法。

内容

前提の確認

任意の mについて、隣接 k項間の漸化式が

\begin{eqnarray} a_{m+k} = \sum_{i=0}^{k-1}c_{i}a_{m+i} \end{eqnarray}

で与えられているとする。このとき、

\begin{eqnarray} a_{n} = \sum_{i=0}^{k-1} C(n, i)a_{i} \end{eqnarray}

となる係数 C(n, i)を、 C(k, i) (i=0, 1, \dots k-1)で表すことが最終的なゴールである。ここで、任意の0以上の整数 mについて

\begin{eqnarray} a_{n+m} = \sum_{i=0}^{k-1} C(n, i)a_{i+m} \ \ \ \ \ (1) \end{eqnarray}

となることを言及しておく。

(A) n \rightarrow n+1の関係

さて、 a_{n+1}について考えていくと

\begin{eqnarray} a_{n+1} &=& \sum_{i=0}^{k-1}C(n+1, i)a_{i} \\ &=& \sum_{i=0}^{k-1}C(n, i)a_{i+1} \ \ \ \ \ (1)より \\ &=& \sum_{i=0}^{k-2}C(n, i)a_{i+1} + C(n, k-1)a_{k} \\ &=& \sum_{i=1}^{k-1}C(n, i-1)a_{i} + C(n, k-1)\sum_{i=0}^{k-1}C(k, i)a_{i} \\ &=& C(n, k-1)C(k, 0)a_{0} + \sum_{i=0}^{k-1}\left[ C(N, i-1)+ C(N, k-1)C(k, i) \right]a_{i} \\ \end{eqnarray}

となるから、

\begin{eqnarray} C(n+1, 0) &=& C(n, k-1)C(k, 0) \\ C(n+1, i) &=& C(n, k-1)C(k, i) \ \ \ (i>0) \end{eqnarray}

という関係式を導くことができた。

(B) n \rightarrow 2nの関係

次に a_{2n}について考えてみると

\begin{eqnarray} a_{2n} &=& \sum_{j=0}^{k-1} C(2n, j)a_{j} \\ &=& \sum_{j=0}^{k-1} C(n, j) \sum_{i=0}^{k-1}C(n+j, i)a_{i} \\ &=& \sum_{i=0}^{k-1} \left[ \sum_{j=0}^{k-1} C(n, j) C(n+j, i) \right] a_{i} \\ \end{eqnarray}

となるから

\begin{eqnarray} C(2n, i) = \sum_{j=0}^{k-1}C(n, j)C(n+j, i) \end{eqnarray}

という関係式を導くことができた。

解法

(A)と(B)の漸化式を使い、 C(2n, i) (i=0. 1, k-1)を C(k, j) (j=0, 1, k-1)で表すために、以下のようなステップで \mathcal{O})(\log n)の計算量で求めることができる。ここで、基本的には現在の値が奇数の時は(A)、偶数の場合は(B)を使って遷移させることになるが、値が 2kを下回った場合は常に(A)を使って遷移させることに注意。

f:id:salpik:20220414225242p:plain
k=100、 n=819の場合の遷移の例

計算量

(A)の1回あたりの計算量は \mathcal{O}(k)であり、(B)の1回あたりの計算量は \mathcal{O}(k^{2})である。 nから kまで減らしていくステップで、(A)は \mathcal{O}(k + \log n)回、(B)は \mathcal{O}(\log n)回なので、合わせて \mathcal{O}(k(k+\log n) + k^{2}\log n) = \mathcal{O}(k^{2}\log n)である。

参考

smijake3.hatenablog.com

画像特徴量

特徴量

SIFT

Scale Invariant Feature Transformの略。2004年発表。Difference of Gaussian (DoG)の極値を探すことでエッジ点を検出。回転、拡大縮小、限定アフィン変換に対して頑強だが、計算コストがかかる。

SURF

Speed Up Robust Featureの略。2008年発表。回転、拡大縮小に対して頑強だが、アフィン変換に対して弱い。SIFTと比べると高速。

KAZE

2012年発表。回転、拡大縮小、限定アフィン変換に対して頑強。

AKAZE

Accelerated-KAZE。2013年発表。回転、拡大縮小、限定アフィン変換に対して頑強。

BRIEF

Binary Robust Independent Elementary Featuresの略。 特徴点の記述に使われ、検出には使えない。

ORB

Oriented FAST and Rotated BRIEF の略。2013年発表。回転、拡大縮小、限定アフィン変換に対して頑強。

BRISK

Binary Robust Invariant Scalable Keypointsの略。2011年発表。

FAST

Features from Accelerated Segment Testの略。

Harris

コーナー検出。1988年

http://www.bmva.org/bmvc/1988/avc-88-023.pdf

参考

ieeexplore.ieee.org

Harrisコーナー検出。

u, vを適当な変位とする。画像上のある点 Pがコーナーであるかを判定するために、点 P周りのみ0でない値を持つ窓関数 w_{P}(x, y)を用いて

\begin{eqnarray} E_{P} &\equiv& \sum_{x, y} w_{P}(x, y) \left[ I(x+u, y+v) - I(x, y) \right]^{2} \\ \end{eqnarray}

を使い、この値が大きい点 Pでコーナーと判定する。

f:id:salpik:20220403224050p:plain
窓関数のイメージ

E_{P}を変形すると

\begin{eqnarray} E_{P} &=& \sum_{x, y} w_{P}(x, y) \left[ \frac{\partial I}{\partial x} u + \frac{\partial I}{\partial y} v \right]^{2} \\ &=& (u, v) \sum_{x, y} w_{P}(x, y) \left( \begin{array}{cc} (\frac{\partial I}{\partial x})^{2} & \frac{\partial I}{\partial x} \frac{\partial I}{\partial y} \\ \frac{\partial I}{\partial x} \frac{\partial I}{\partial y} & (\frac{\partial I}{\partial y})^{2} \\ \end{array} \right) \left( \begin{array}{c} u \\ v \end{array} \right) \\ &\equiv& (u, v) M \left( \begin{array}{c} u \\ v \end{array} \right) \end{eqnarray}

となる。ここで行列 Mを M = P^{-1}\Sigma Pと対角化したとすると

\begin{eqnarray} E_{P} &=& (u, v) P^{-1}\Sigma P \left( \begin{array}{c} u \\ v \end{array} \right) \\ &=& \left[ P \left( \begin{array}{c} u \\ v \end{array} \right) \right]^{-1} \Sigma \left[ P \left( \begin{array}{c} u \\ v \end{array} \right) \right] \end{eqnarray}

と変形できる。これは、 \left( \begin{array}{c} u^{'} \ v^{'} \end{array} \right) \equiv P\left( \begin{array}{c} u \ v \end{array} \right)という新しい直行軸 u^{'}, v^{'}を取ったときに、それぞれの方向にどれだけ変化するかを表しており、その変化係数がそれぞれの固有値である。。新しい直行軸は変化の主成分を表しているため、以下のことが言える。

  • 第1、第2固有値のどちらも大きい：頂点である
  • 第1固有値のみ大きく、第2固有値は小さい：辺である
  • 第1、第2固有値のどちらも小さい：辺でも頂点でもない

具体例

コーナーの場合

窓関数領域において、 x方向と y方向がそれぞれ独立に大きな微分 \frac{\partial I}{\partial x}、 \frac{\partial I}{\partial y}を取るとする。そうすると (\frac{\partial I}{\partial x})^{2}や (\frac{\partial I}{\partial y})^{2}は大きな値をとるが \frac{\partial I}{\partial x}\frac{\partial I}{\partial y}はある程度互いに独立になり0に近い値になるので、行列 Mは対角行列に近くなり、固有ベクトルの方向も x軸と y軸に平行になる。

f:id:salpik:20220403224921p:plain
行列 Mの計算例

辺の場合

例えばx軸に平行な辺の場合、 |\frac{\partial I}{\partial y}| >> |\frac{\partial I}{\partial x}|である。 (\frac{\partial I}{\partial y})^{2}は (\frac{\partial I}{\partial x})^{2}と比べると大きくなり、 \frac{\partial I}{\partial x} \frac{\partial I}{\partial y}は互いにキャンセルして0に近くなる。 (\frac{\partial I}{\partial x})^{2}も小さい値であるが常に2乗で足していっているので \frac{\partial I}{\partial x} \frac{\partial I}{\partial y}の項よりは大きくなると思うと、やはり行列 Mは対角行列となり、第1固有値は |\frac{\partial I}{\partial y}|で固有ベクトルが[tex: y軸方向、第2固有値は |\frac{\partial I}{\partial x}|で固有ベクトルが[tex: x軸方向となる。

f:id:salpik:20220403225318p:plain
固有ベクトルの方向とその大きさのイメージ

参考

labs.eecs.tottori-u.ac.jp

SIFT

画像ピラミッド

画像ピラミッドとは、重み関数 w(m, n)を定義し、この重みに基づいて圧縮した画像をピラミッドのように並べたものである。元の画像のサイズを N\times Mとし、 G_{0}(x, y)\ (0 \le x \lt N, 0 \le y \lt M)とすると、下から i\ (i=0, 1, 2, \dots)層目のピラミッドの画像は \frac{N}{2^{i}} \times \frac{M}{2^{i}}であり、

\begin{eqnarray} G_{i+1}(x, y) = \sum_{m, n} w(m, n) G_{i}(2x+m, 2y+n) \end{eqnarray}

と表される。重み関数としてはガウシアンなどが使われる。またこのように画像サイズを圧縮することをOctaveと呼ぶ

f:id:salpik:20220327221720p:plain
画像ピラミッドのイメージ

参考

medium.com

labs.eecs.tottori-u.ac.jp

画像フィルターと微分の関係

画像の微分はフィルターをかけることと同値である。

1階微分

例えば I(x, y)のピクセルにおいて xと y方向の微分を考える。

\begin{eqnarray} I(x+1, y) &\sim& I(x, y) + I_{x}(x, y) \\ I(x-1, y) &\sim& I(x, y) - I_{x}(x, y) \\ \end{eqnarray}

となるから、

\begin{eqnarray} I_{x}(x, y) = I(x+1, y) - I(x-1, y) \end{eqnarray}

となるから、これは (-1\ 0\ 1)というフィルターを適用したのと同等である。

2階微分

\begin{eqnarray} I(x+1, y) &\sim& I(x, y) + I_{x}(x, y) + \frac{1}{2}I_{xx}(x, y) \\ I(x-1, y) &\sim& I(x, y) - I_{x}(x, y) + \frac{1}{2}I_{xx}(x, y) \\ \end{eqnarray}

となるから、

\begin{eqnarray} I_{xx}(x, y) = I(x-1, y) -2 I(x, y) + I(x+1, y) \end{eqnarray}

となり、 (1\ -2\ 1)というフィルターを適用したのと同等である。2次元に拡張すれば、2階微分の和であるラプラシアンは

\begin{eqnarray} \Delta \equiv I_{xx} + I_{yy} = \left( \begin{array}{ccc} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0 \\ \end{array} \right) \end{eqnarray}

というフィルターとなる。

エッジとラプラシアンの関係

エッジ位置というのは2階微分が0でありかつ1階微分の絶対値が大きいところとして定義できる。したがってラプラシアン \Delta \equiv I_{xx}+I_{yy}のフィルタをかけこの値の絶対値が小さいところをエッジ点とすればよい。

f:id:salpik:20220327224519p:plain
エッジ位置と1階微分、2階微分の関係（点線がエッジ位置）

参考

https://www.cse.psu.edu/~rtc12/CSE486/lecture11.pdf

LoGとDoGの関係

まず、ガウシアンフィルターのラプラシアンは、ガウシアンのラプラシアンのフィルターになる。すなわち、 G(x, y; \sigma)をガウシアンフィルタ

\begin{eqnarray} G(x, y; \sigma) \equiv \frac{1}{2\pi \sigma^{2}} \exp(-\frac{x^{2}+y^{2}}{2\sigma^{2}}) \end{eqnarray}

とし、ピクセルの値を I(x, y)とすると、

\begin{eqnarray} \Delta [ G(x, y; \sigma) \otimes I(x, y) ] = [ \Delta G(x, y; \sigma) ] \otimes I(x, y) \end{eqnarray}

となる。

証明

\begin{eqnarray} \Delta [ G(x, y; \sigma) \otimes I(x, y) ] &=& (\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}) \sum_{X, Y} G(x-X, y-Y)I(X, Y) \\ &=& \sum_{X, Y} [ (\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}})G(x-X, y-Y) ] I(X, Y) \\ &=& [ \Delta G(x, y, \sigma) ] \otimes I(x, y) \end{eqnarray}

1次元の場合は、

\begin{eqnarray} G(x; \sigma) = \frac{1}{2\pi \sigma^{2}} \exp(-\frac{x^{2}}{2\sigma^{2}}) \end{eqnarray}

となるので、

\begin{eqnarray} \frac{\partial G(x; \sigma)}{\partial x} &=& \frac{1}{2\pi \sigma^{2}} (-\frac{x}{\sigma^{2}})\exp(-\frac{x^{2}}{2\sigma^{2}}) \\ \frac{\partial^{2} G(x; \sigma)}{\partial x^{2}} &=& \frac{1}{2\pi \sigma^{2}} \left(-\frac{1}{\sigma^{2}} \right) \left(1-\frac{x^{2}}{\sigma^{2}}\right) \exp(-\frac{x^{2}}{2\sigma^{2}}) \end{eqnarray}

f:id:salpik:20220327231349p:plain
ガウシアンの1階微分、2階微分

一方、異なる標準偏差のガウシアンを引き算すると、ガウシアンの2階微分の形に近くなる。

f:id:salpik:20220327232930p:plain
ガウシアンの引き算
これが、LoG(Laplacian of Gaussian filter)をDoG(Difference of Gaussian filter)で近似するということである。すなわち、 \sigmaの異なるガウシアンフィルターを使って同じ画像サイズの画像をいくつか生成し、その差分をとることでラプラシアンフィルターの近似を得て、その絶対値が小さいところとしてエッジを検出する。また、その時に、1)同じ \sigmaのフィルターをかけた画像内での周辺8ピクセルおよび2)前後の \sigmaのフィルターをかけた画像内での周辺および自分の位置9ピクセルを使って、異なる \sigmaのフィルターをかけた画像においてもエッジとなっているかを確認したのち正式な特徴点として登録される。

f:id:salpik:20220327235213p:plain
赤い対象ピクセルが特徴点として登録するか判断する時にチェックされる周辺ピクセル

参考

https://www.cse.psu.edu/~rtc12/CSE486/lecture11.pdf

Descriptorの計算

特徴点の周辺を 4 \times 4=16領域に分割し、各領域ごとに 45^\circごとに合計8方向の特徴ベクトルを用意する。したがってSIFTの1つの特徴点は合計128次元のベクトルである。

マルチスレッド

マルチスレッドの例

C++

#include <thread>

void thread1() {
  printf("thread1\n");
}

void thread2() {
  printf("thread2\n");
}

int main() {
  std::thread th1(thread1);
  std::thread th2(thread2);
  // joinを呼ぶと、そのスレッドが終了するまで現在のコンテクストが中断する
  th1.join();
  th2.join();
  return 0;
}

Python

import threading

def thread1():
    print("thread1")

def thread2():
    print("thread2")

t1 = threading.Thread(target=thread1)
t2 = threading.Thread(target=thread2)
print("start")
t1.start()
t2.start()
t1.join()
t2.join()
print("end")

グローバル変数を複数スレッドから書き込む失敗例

#include <thread>

int cnt = 0;

void thread1() {
  printf("thread1 start cnt = %d\n", cnt);
  for (int i=0; i<100000; i++) {
    cnt++;
  }
  printf("thread1 end cnt = %d\n", cnt);
}

void thread2() {
  printf("thread2 start cnt = %d\n", cnt);
  for (int i=0; i<100000; i++) {
    cnt++;
  }
  printf("thread2 end cnt = %d\n", cnt);
}

int main() {
  std::thread th1(thread1);
  std::thread th2(thread2);
  // joinを呼ぶと、そのスレッドが終了するまで現在のコンテクストが中断する
  th1.join();
  th2.join();
  return 0;
}

出力例

thread1 start cnt = 0
thread2 start cnt = 422
thread2 end cnt = 77428
thread2 end cnt = 115025

排他処理を入れてグローバル変数を更新する例

C++

#include <thread>
#include <mutex>

int cnt = 0;
std::mutex mtx;

void add_count() {
  std::lock_guard<std::mutex> lock(mtx);
  cnt++;
}

void thread1() {
  printf("thread1 start cnt = %d\n", cnt);
  for (int i=0; i<100000; i++) {
    add_count();
  }
  printf("thread1 end cnt = %d\n", cnt);
}

void thread2() {
  printf("thread2 start cnt = %d\n", cnt);
  for (int i=0; i<100000; i++) {
    add_count();
  }
  printf("thread2 end cnt = %d\n", cnt);
}

int main() {
  std::thread th1(thread1);
  std::thread th2(thread2);
  // joinを呼ぶと、そのスレッドが終了するまで現在のコンテクストが中断する
  th1.join();
  th2.join();
  return 0;
}

出力例

thread1 start cnt = 0
thread2 start cnt = 0
thread2 end cnt = 195262
thread1 end cnt = 200000

最終的にグローバル変数が意図通りの値になっている。

Python

import threading
import time

lock = threading.Lock()
value = 0

def thread(n):
    lock.acquire()
    global value
    value = value + n
    value = value - n
    lock.release()

print(f"start value={value}")
ts = []
N = 4
for i in range(N):
    t = threading.Thread(target=thread, args=(i,))
    t.start()
    ts.append(t)

for i in range(N):
    ts[i].join()

print(f"end value={value}")

セマフォ

import threading
import time

# Semaphoreというオブジェクトもある                                                                                          
# Semaphore: acquireせずにreleaseして残セマフォ数を初期値以上に増やすことができる                                            
# BoundedSemaphore: qcquireしないとreleaseできない                                                                           
s = threading.BoundedSemaphore(2)

def thread(n):
    s.acquire()
    time.sleep(5)
    print(f"thread: {n}")
    s.release()

# こちらの獲得方法も可能                                                                                                     
def thread2(n):
    with s:
        time.sleep(5)
        print(f'thread2: {n}')

print("start")
for i in range(10):
    t = threading.Thread(target=thread, args=(i,))
    t.start()
print("end")

とやると、文字が2つずつ画面に出力される。

実行されるタイミングの確認

std::threadをコールした段階でマルチスレッドが立ち上がりスレッド処理が始まる。joinをコールした段階で開始するわけではないことに注意。

#include <thread>
#include <mutex>
#include <unistd.h>

int cnt = 0;
std::mutex mtx;

void add_count() {
  std::lock_guard<std::mutex> lock(mtx);
  cnt++;
}

void thread1() {
  printf("thread1 start cnt = %d\n", cnt);
  for (int i=0; i<100000; i++) {
    add_count();
  }
  printf("thread1 end cnt = %d\n", cnt);
}

void thread2() {
  printf("thread2 start cnt = %d\n", cnt);
  for (int i=0; i<100000; i++) {
    add_count();
  }
  printf("thread2 end cnt = %d\n", cnt);
}

int main() {
  std::thread th1(thread1);
  printf("cnt before sleep = %d\n", cnt);
  usleep(1000000); // 1 sec
  printf("cnt after sleep = %d\n", cnt);
  std::thread th2(thread2);
  // joinを呼ぶと、そのスレッドが終了するまで現在のコンテクストが中断する
  th1.join();
  th2.join();
  return 0;
}

出力例

cnt before sleep = 0
thread1 start cnt = 0
thread1 end cnt = 100000
cnt after sleep = 100000
thread2 start cnt = 100000
thread2 end cnt = 200000

usleep直後の段階ではまだth1.join()は呼んでいないが、thread1の関数の処理が実行完了していることがわかる。

マルチスレッドで実行する関数の返り値を設定する

PythonではThreadを継承したクラスを作りrunメソッドをオーバーロードすればよい

import threading
import time

def myfunc(n):
    ans = 0
    for i in range(n):
        ans += i
    return ans

class MyThread(threading.Thread):
    def __init__(self, func, args):
    # 親クラスのコンストラクタを呼ぶことを忘れない                                                     
        super().__init__()
        self.func = func
        self.args = args

    def run(self):
        self.result = self.func(*self.args)

    def get(self):
        return self.result

t1 = MyThread(myfunc, args=(10,))
t2 = MyThread(myfunc, args=(100,))
t1.start()
t2.start()
t1.join()
t2.join()

print(f"result1={t1.get()}")
print(f"result2={t2.get()}")

参考

Pythonマルチスレッドの戻り値を取得する2つの方法 - JPDEBUG.COM

デッドロック

Pythonでのデッドロックの例

import threading
import time

l1 = threading.Lock()
l2 = threading.Lock()

def func1():
    l1.acquire()
    time.sleep(2) # 時間のかかる処理を模擬                                                                                   
    g1 = 10
    l2.acquire()
    g2 = 10
    l2.release()
    l1.release()

def func2():
    l2.acquire()
    time.sleep(2) # 時間のかかる処理を模擬                                                                                   
    g2 = 10
    l1.acquire()
    g1 = 10
    l1.release()
    l2.release()

t1 = threading.Thread(target=func1)
t2 = threading.Thread(target=func2)
t1.start()
t2.start()
t1.join()
t2.join()

デッドロック対策

  • プロセスが利用するリソースをまとめて1つのロックで管理するようにし、実行開始前にリソース全体のロックを取得する。これによりリソースの利用率は低下する。
  • リソースのロックをかける順番を揃える。

参考

デッドロックと回避策

コメント

  • 複数スレッドから同じ変数の値を書き出す時だけでなく、読み込む時にむ排他処理を入れた方が良い。読み込み中に他スレッドによって値が更新されると値がundefinedになることがある。

参考

Pythonのthreadingとmultiprocessingを完全理解 - Qiita

高速フーリエ変換

問題概要

A_{i}\ \ (i=1, 2, \dots, N)、 B_{i}\ \ (i=1, 2, \dots, N)が与えられたとき、

\begin{eqnarray} \sum_{i, j, i+j=k} A_{i}B_{j} \end{eqnarray}

を高速に求めたい。普通に計算すると \mathcal{O}(N^{2})の計算量となるが、これを \mathcal{O}(N\log N)に高速化することを目指す。

問題の詳細は以下を参照 atcoder.jp

フーリエ変換

ある N-1次多項式 F(x) = f(0) + f(1)x + f(2)x^{2} + \dots f(N-1)x^{N-1}があったとして、この多項式の N個の xに対する値 F(x_{0}), F(x_{2}), \dots F(x_{N-1})がわかれば係数が決定できるはずである。適当に x_{0}, x_{1}, x_{N-1}を選んでしまうと F(x_{0}), F(x_{2}), \dots F(x_{N-1})と f(0), f(1), \dots f(N-1)の間の関係式は単純にならないが、1の N乗根を選ぶと直行性により綺麗な関係になるということである。

\xi_{N} \equiv e^{i\frac{2\pi}{N}}と定義すると

\begin{eqnarray} f(x) &=& \frac{1}{N}\sum_{k=0}^{N-1}\hat{f}(\xi_{N}^{-k}) x^{k} \\ \hat{f}(t) &=& \sum_{k=0}^{N-1}f(\xi_{N}^{k}) t^{k} \\ \end{eqnarray}

が成立する。

証明

\begin{eqnarray} f(x) = \frac{1}{N}\sum_{k=0}^{N-1}\hat{f}(\xi_{N}^{-k}) x^{k} \end{eqnarray}

より、

\begin{eqnarray} f(\xi_{N}^{m}) &=& \frac{1}{N}\sum_{k=0}^{N-1} \sum_{k^{'}=0}^{N-1} f(\xi_{N}^{k^{'}}) (\xi_{N}^{-k})^{k^{'}} (\xi_{N}^{m})^{k} \\ &=& \frac{1}{N}\sum_{k=0}^{N-1} \sum_{k^{'}=0}^{N-1} f(\xi_{N}^{k^{'}}) \xi_{N}^{k(m-k^{'})} \\ &=& \frac{1}{N} \sum_{k^{'}=0}^{N-1} f(\xi_{N}^{k^{'}}) \sum_{k=0}^{N-1} \xi_{N}^{k(m-k^{'})} \\ &=& \frac{1}{N} \sum_{k^{'}=0}^{N-1} f(\xi_{N}^{k^{'}}) N\delta (m-k^{'}) \\ &=& f(\xi_{N}^{m}) \end{eqnarray}

よって、 f(\xi_{N}^{m}) \ \ (m=0, 2, \dots , N-1)の異なる Nこの点で関数 fの値が一致することが示されたので、この変換の妥当性が示された。

関数の値の効率の良い求め方。

f(x)は N-1次式なので、 m=0, 1, \dots, N-1の Nこの値について

\begin{eqnarray} f(\xi_{N}^{k}) = f[k] = \hat{f}_{0} + \hat{f}_{1}\xi_{N}^{k} + \hat{f}_{2}(\xi_{N}^{k})^{2} + \dots + \hat{f}_{N-1}(\xi_{N}^{k})^{N-1} \end{eqnarray}

を普通に計算しようとすると \mathcal{O}(N^{2})の計算量となる。これを効率よく計算するため部分問題に分割する。 fの項を偶数部分 f^{(0)}と奇数部分 f^{(1)}に分割すると、それぞれ N/2の次数になるので、

\begin{eqnarray} f^{(0)}[k] &=& \hat{f}_{0} + \hat{f}_{2}\xi_{N/2}^{k} + \hat{f}_{4}(\xi_{N/2}^{k})^{2} + \hat{f}_{6}(\xi_{N/2}^{k})^{3} + \dots \\ &=& \hat{f}_{0} + \hat{f}_{2}(\xi_{N}^{k})^{2} + \hat{f}_{4}(\xi_{N}^{k})^{4} + \hat{f}_{6}(\xi_{N}^{k})^{6} + \dots \\ f^{(1)}[k] &=& \hat{f}_{1} + \hat{f}_{3}\xi_{N/2}^{k} + \hat{f}_{5}(\xi_{N/2}^{k})^{2} + \hat{f}_{7}(\xi_{N/2}^{k})^{3} + \dots \\ &=& \hat{f}_{1} + \hat{f}_{3}(\xi_{N}^{k})^{2} + \hat{f}_{5}(\xi_{N}^{k})^{4} + \hat{f}_{7}(\xi_{N}^{k})^{6} + \dots \\ \end{eqnarray}

したがって、元の定義と見比べれば

\begin{eqnarray} f[k] = f^{(0)}[k] + \xi_{N}^{k} f^{(1)}[k] \end{eqnarray}

と表されることがわかる。

計算量の評価

k=0, 1, \dots, N-1について、 f[k ]を求めるのにかかる計算量を T(N)とする。上で説明した解法により、サイズが半分になった問題 T(N/2)を2つ解き、それを k=0, 1, \dots, N-1について足し合わせればよいので

\begin{eqnarray} T(N) \sim 2T(N/2) + \mathcal{O}(N) \end{eqnarray}

という漸化式が成立する。ここから

\begin{eqnarray} T(N) &\sim& 2T(N/2) + \mathcal{O}(N) \\ &\sim& 2(2T(N/4)+ \mathcal{O}(N/2)) + \mathcal{O}(N) \\ &\sim& 4T(N/4) + 2 \mathcal{O}(N) \\ &\sim& 4(2T(N/8)+ \mathcal{O}(N/4)) + 2 \mathcal{O}(N) \\ &\sim& 8T(N/8) + 3 \mathcal{O}(N) \\ \vdots \\ &\sim& NT(1) + \log N \mathcal{O}(N) \\ &\sim& \mathcal{O}(N\log N) \end{eqnarray}

参考

https://qiita.com/ageprocpp/items/0d63d4ed80de4a35fe79