線形代数ライブラリEigenの遅延評価について

EigenはC++で書かれた線形代数ライブラリである。特徴の一つとして、式の実装にtemplate機能を活用している点があり、これにより 遅延評価 (lazy evaluation) を実現している。遅延評価を適切に使うことで、計算コストを削減できることがある。

Eigenの遅延評価については公式ドキュメントの次のページが詳しいため、本記事では厳密性・網羅性よりも概要の理解しやすさを重視する。

なお、以下のコード例では

using namespace Eigen;

が暗黙的に実行されているとする。

Eigenにおける遅延評価

演算の返り値の型に注目する

行列を表す MatrixXd 型の変数 a があり、その転置を表現するために

MatrixXd a = ...;

... = a.transpose();

と記述したとする。ただし ... は適当な実装が入るとする。ナイーブに考えれば、メンバ関数呼び出し a.transpose() の戻り値は、列数と行数は入れ替わるものの、「行列」であると期待される。実際、例えば

MatrixXd b = a.transpose();

と別の行列 b に代入しても期待通りに動作する。

しかし、厳密にその型について考えると、少し見方が変わる。ここで呼ばれている aメンバ関数 transpose の定義はこのページ(最終アクセス: 2020/02/01) にある。その記述によると、返り値の型は a の型である MatrixXd ではなく、 Transpose<MatrixXd> という別の型であると分かる。これは、演算結果として得られる行列そのものではなく、転置という式を表す型である。

つまり、上記のように別の MatrixXd 型の行列オブジェクトに転置の演算結果を代入するという例は、

MatrixXd b = static_cast<MatrixXd>(a.transpose());

のような、式から行列への暗黙的なキャストを実行していたのだと解釈できる。

遅延評価の観点で重要なのは、MatrixXdメンバ関数 transpose の呼び出し自体には計算コストがかからないという点である。実際に転置の操作が実行される(式が評価される)のは、転置の結果が必要になったときである。上の例では static_cast 内で Transpose<MatrixXd> 型を MatrixXd 型に変換する際に初めて転置に相当する操作が実行されている(式が評価されている)はずである。

遅延評価の利点

遅延評価によって、余計な計算コストを削減できる場合がある。例えば以下のような行列の転置と行列の積を考える。

MatrixXd a = ...;
MatrixXd b = ...;

... = a.transpose() * b;

この例では遅延評価の効果により、いわゆる転置の操作(行列の行と列を入れ替えた行列オブジェクトの生成)が直ちに実行されるわけではない。その変わり、 Transpose<MatrixXd> 型のメンバ関数 operator* が呼ばれ、 Transpose<MatrixXd> 型と MatrixXd 型の積を直接計算するための操作が実行される(もっといえば、実際には Product<Transpose<MatrixXd>, MatrixXd> という型の式が生成され、必要があるまで積の操作は計算されない)ということになる。

もしこれを

MatrixXd c = a.transpose();

... = c * b;

のように、転置の操作(行列の行と列を入れ替えた行列オブジェクトの生成)をまず実行し、そのあとで MatrixXd 型同士の積を計算するようにしたとすると、転置の操作がある分だけ計算コストが増えてしまうと考えられる。

遅延評価の挙動を意識して適切にコーディングすることで、計算コストを削減することができる。逆に、上記のような無用な型のキャスト(無用なタイミングでの式の評価)を行ってしまうと、遅延評価の利点を得られず計算コストが増えてしまう場合がある。

遅延評価を活用するコーディング

型推論で無用なキャストを避ける

C++11から導入された型推論のための auto キーワードを使うことで、容易に無用な型のキャスト(無用なタイミングでの式の評価)を避けることができる場合がある。

例えば

MatrixXd a = ...;
MatrixXd b = ...;

auto c = a.transpose() * b;

とすれば、変数 c の型はコンパイラによって Product<Transpose<MatrixXd>, MatrixXd> のように推論されると思われる。注意点として、変数 c はあたかも行列のように振舞うが、実際には評価されるのを待っている式であるという点は常に意識しなくてはならない。

実装例と簡単なパフォーマンス比較

例として、以下のような関数を考える。なおこの関数自体はナンセンスで、あくまで計算の様子を明らかにするための例にすぎない。

Func1

明示的に const MatrixXd 型を指定した中間オブジェクトを3つ生成してしまっており、遅延評価の恩恵を受けられていない。

double Func1(const MatrixXd& a)
{
    const MatrixXd b = a.transpose();
    const MatrixXd c = MatrixXd::Ones(a.cols(), a.rows());
    const MatrixXd d = b + c;

    return d.maxCoeff();
}
Func2

最後の和の演算を auto で式として受けることで、中間オブジェクトの生成を抑制しようとしている。

double Func2(const MatrixXd& a)
{
    const MatrixXd b = a.transpose();
    const MatrixXd c = MatrixXd::Ones(a.cols(), a.rows());
    const auto     d = b + c;

    return d.maxCoeff();
}
Func3

さらに cauto で宣言することで、定数で埋められた行列であるということを表す型( MatrixXd とは異なる)のまま扱おうとしている。

double Func3(const MatrixXd& a)
{
    const MatrixXd b = a.transpose();
    const auto     c = MatrixXd::Ones(a.cols(), a.rows());
    const auto     d = b + c;

    return d.maxCoeff();
}
Func4

本記事の冒頭の例のように、転置を行った結果を評価し MatrixXd 型の中間オブジェクトとして明示的に生成するのではなく、式としてそのまま受け取ろうとしている。

double Func4(const MatrixXd& a)
{
    const auto b = a.transpose();
    const auto c = MatrixXd::Ones(a.cols(), a.rows());
    const auto d = b + c;

    return d.maxCoeff();
}
Func5

最後に、ここでは変数の宣言をやめて一行で演算を書こうとしている。これによって自動的に適切な型が推論されて遅延評価の利点が得られることが期待できる。

double Func5(const MatrixXd& a)
{
    return (a.transpose() + MatrixXd::Ones(a.cols(), a.rows())).maxCoeff();
}
計測結果

以上の5つの異なる実装を実行した際の計算時間の例を下記に示す。ただし値は500回計算を行った際の平均を表し、エラーバーはその際の標準偏差を表す。

これらの例では auto を使って無用な型のキャストを避けるほど実行性能が向上している。なお Func4 (全てに auto を使用)と Func5 (途中変数を定義せずに一行で記述)は、この場合は実質的に同じ処理が走っていると思われる。

一行で記述できるときはそれで良いが、一行で書くのが憚られるほど長い式や、同じ部分式が複数回含まれるような式を扱う際には、式が評価されるタイミングを意識しつつ、うまく auto を活用して実装すると良いかもしれない。

おまけ:遅延評価を活用した関数の定義

上の例では関数の引数を const MatrixXd& としているため、関数呼び出し側で暗黙的なキャストが発生してしまい、遅延評価の恩恵が失われてしまうことがある。

下記のように MatrixBase クラスを活用して引数の宣言をすることで、この問題を解決できることがある。

template <typename Derived>
double Func6(const MatrixBase<Derived>& a)
{
    const auto b = a.transpose();
    const auto c = MatrixBase<Derived>::Ones(a.cols(), a.rows());
    const auto d = b + c;

    return d.maxCoeff();
}

このことは公式ドキュメントの次のページが詳しい。

最後に

C++についてもEigenについても大した知識なく本記事を書いているので、致命的な勘違いが含まれている可能性がある。何か気付いた方はご指摘いただけると大変ありがたい。

比較データと確率モデル

n 個のアイテムからなる集合 \mathcal{P} = \{ 1, \ldots, n \} を考える。本記事では、これらのアイテム間の比較に基づく観測データ(例えば二者間の勝敗など)に対して、適当な確率モデルを導入することで、その事象を説明することについて考える。また、その確率モデルの潜在変数を、観測データを元に推定することについて考える。

前提:Luce's Choice Axiom

この記事では特に、 Luce's choice axiom (ルースの選択公理) [Luce08] に基づく確率モデルのみを考える(公理の定義はこの記事では触れない)。

この公理のもとでは、各アイテムに紐づけられた潜在的なパラメタ  v_{1}, \ldots, v_{n} > 0 (これらの値が 0 である場合もあるが、この記事では簡単のため正の値として考える。)が存在し、比較はこれらの値に基づいて行われる。

なお、このパラメタは論文や文脈により異なる呼ばれ方をする。例えば "skill rating" [Caron+12] 、 "goodness value" [Koyama+17] 、 "quality scale" [Tsukida+11] などがある。

本題:比較データの形式と主な確率モデル

比較データにはいくつかの形式が考えられる。それぞれについて、それを説明するための確率モデルが考案されている。

なお、この節の内容は http://choix.lum.li/en/latest/data.html に簡潔にまとまっている。

Pairwise Comparison(一対比較)

Pairwise comparisonは日本語では一対比較などと呼ばれる。Paired comparisonとも表記される(例えば [神嶌09; Handley01] など)。これは2つのアイテムが比較された結果どちらかが選択されたという事象を表すデータ形式である。

アイテム i がアイテム j と比較された結果 i が選ばれるという事象は

\displaystyle i \succ j

などと表現される。

Pairwise comparisonのデータを扱う確率モデルとしては Bradley-Terryモデル がある [Brochu+07; Hunter04; Tsukida+11] 。このモデルは、上記の事象が起こる確率を

\displaystyle\text{Pr}(i \succ j) = \frac{v_{i}}{v_{i} + v_{j}}

とするものである。

指数関数やロジスティック関数を用いた表現

Bradley-Terryモデルは、その潜在変数を

v_{i} = \exp(\beta_{i}) \:\:\:\: (i = 1, \ldots, n)

として

\displaystyle\text{Pr}(i \succ j) = \frac{\exp(\beta_{i})}{\exp(\beta_{i}) + \exp(\beta_{j})}

と表現することも多い。この表現を用いてさらに式変形をすると、 \sigma を標準ロジスティック関数として

\displaystyle\text{Pr}(i \succ j) = \frac{1}{\exp(- (\beta_{i} - \beta_{j})) + 1} = \sigma(\beta_{i} - \beta_{j})

と表現することもできる。

en.wikipedia.org

その他のモデル

Pairwise comparisonを扱う確率モデルとしては、他にThurstone-Mostellerモデル [Chu+05; Handley01; Tsukida+11] が存在し、Bradley-Terryモデルと同様の文脈で議論されることがある [Handley01] 。Thurstone-Mostellerモデルは "Thurstone's Law of Comparative Judgement, Case V" を指している [Handley01] 。Bradley-Terryモデルは、Thurstone-Mostellerモデルが仮定しているGaussian分布をGumbel分布に置き換えたものと解釈することが可能である [Tsukida+11] 。

Generalized Bradley-Terryモデルというより表現力の高い確率モデルも存在する [Caron+12; Hunter04] 。

Top-1 List

Top-1 listはあまり共通認識のある呼称ではないかもしれないが、ここでは http://choix.lum.li/ に倣い、このように呼ぶことにする。これは、2つ以上のアイテムが同時に比較された結果ある1つのアイテムが選択されたという事象を表すデータ形式である。

アイテム群 \{ i, j, \ldots, k \} からアイテム i が選ばれるという事象は

\displaystyle i \succ \{ j, \ldots, k \}

などと表現される [Koyama+17] 。

Top-1 listのデータを扱う確率モデルとしては Bradley-Terry-Luceモデル がある [Koyama+17; Luce08; Tsukida+11] 。このモデルでは、上記の事象が起こる確率を

\displaystyle\text{Pr}(i \succ \{ j, \ldots, k \}) = \frac{v_{i}}{v_{i} + v_{j} + \cdots + v_{k}}

とする。

なお、比較対象がちょうど2つのときはBradley-Terryモデルと実質的に同じになる。

Partial Ranking

Partial rankingは、2つ以上のアイテムについて、それらの順序が与えられたという事象を表すデータ形式である。

アイテム群  \{ i, j, \ldots, k \} に対し、 i \rightarrow j \rightarrow \cdots \rightarrow k という順序が与えられるという事象は

\displaystyle i \succ j \succ \cdots \succ k

などと表現される [神嶌09] 。

Partial rankingを扱う確率モデルとしては Plackett-Luceモデル がある [Guiver+09; 神嶌09] 。このモデルでは、上記の事象が起こる確率を

\displaystyle\text{Pr}(i \succ j \succ \cdots \succ k) = \frac{v_{i}}{v_{i} + v_{j} + \cdots + v_{k}} \cdot \frac{v_{j}}{v_{j} + \cdots + v_{k}} \cdot \cdots \cdot \frac{v_{k}}{v_{k}}

とする。これは、順位が早いものから逐次的に選択し取り除くという事象をBradley-Terry-Luceモデルを用いて表したものと解釈することができる。

潜在的なパラメタの推定

最尤推定 (maximum likelihood estimation)により、比較データから確率モデルの潜在的なパラメタを推定することができる [Tsukida+11] 。

また、潜在的なパラメタにある事前分布を仮定することで 最大事後確率推定 (maximum a posteriori estimation) により潜在的なパラメタを推定することもできる [Tsukida+11; Caron+12] 。

choix http://choix.lum.li/ というPythonライブラリは、これらの実装を提供している。

また、各アイテム間に類似度を定義できる場合には、ガウス過程などの分布を仮定して最大事後確率推定を行うこともある [Chu+05; Koyama+17] 。

余談:主観的な好み(preference)と比較データ

比較データを扱う確率モデルは preference learning という分野でも活用される(例えば [Chu+05; Brochu+07] など)。これは、人間(評価者)の 主観的な好み(preference) をデータから学習することに関連する分野である。

好みに関するデータを集める上で、あるアイテムに対する好みの評価値をそのまま評価者に聞く方法は望ましくなく、比較に基づくデータを集めるような方法が良いと考えられている [Tsukida+11; Koyama+17; Brochu+10] 。あるアイテムに対する好みの評価値を直接聞く方法が望ましくない理由として、

  • (評価者が複数人いる場合)評価者によってスケールが異なる
  • (複数のアイテムの評価を行う場合)スケールが時間経過と共に変化する(この効果はdriftと呼ばれる)
  • (複数のアイテムの評価を行う場合)早い段階で評価を行ったアイテムによって全体のスケールが影響を受ける(この効果はanchoringと呼ばれる)

などにより、一貫したデータを取得できないことが考えられるためである [Brochu+10] 。

文献

以下に示す文献は本記事を書く上で実際に参考にしたもの(私がもともと良く知っていた文献や、キーワードで検索してたまたま発見した文献など)であって、記事の内容を説明する上で適切な文献を厳選したわけではない。

注意とお願い

私は確率・統計・機械学習等の専門家ではありません。本記事は、私が学びながらとったメモを整形したものに過ぎず、根本的な誤りが含まれるかもしれません。また、些細なことでも構いませんので、誤りの訂正・表現に違和感のある箇所の指摘・改善の提案等ありましたらコメント頂けると幸いです。

共分散行列 (Covariance Matrix) の行列式 (determinant) の対数 (log) を計算する

f:id:yuki-koyama:20190815150220p:plain

下記ライブラリのガウス過程回帰 (Gaussian Process Regression) の実装の際に悩んだ点をメモしておく。

github.com

いつ行列式の対数 (log-determinant) を計算する必要があるか

例えば、多変量正規分布 (multivariate normal distribution) が関連する統計的推定において、対数尤度を計算する際などに、共分散行列の行列式の対数を計算する必要が生じる。

行列式 (determinant) を計算する際の問題点

特にサイズが大きな行列の行列式は大変小さな値をとることがあり、計算機の浮動小数点数の表現力がシビアに効いてくることも多い。これの対数をとると、計算が不安定になったり、infが問題となったりする。

正定値対称行列 (symmetric positive-definite matrix) の行列式の対数 (log-determinant) を直接計算する方法

共分散行列など、計算の対象とする行列が正定値対称行列である場合、行列式を明示的に計算せずに直接行列式の対数を計算することができる。正定値対称行列  \mathbf{K} をコレスキー分解によって

\displaystyle
\mathbf{K} = \mathbf{L}\mathbf{L}^{\textsf{T}}

と分解すると、  \mathbf{K} 行列式の対数は

\displaystyle
\log ( \det ( \mathbf{K} ) ) = 2 \cdot \text{sum} ( \log ( \text{diag2vec} ( \mathbf{L} ) ) )

によって計算することができる。ただし  \text{diag2vec} は対角要素を取り出してベクトルにする変換、右辺の  \log は要素毎の対数、  \text{sum} はベクトルの各要素の和を表す。(この式が成立することの証明は こちら を参照。)

C++でEigenを用いる場合の実装

ナイーブな方法(このようにするべきではない):

// #include <Eigen/LU>
// #include <cmath>

// const Eigen::MatrixXd K = ...
const double log_det_K = std::log(K.determinant());

工夫した方法(こちらの方が好ましいと思われる):

// #include <Eigen/Cholesky>

// const Eigen::MatrixXd K = ...
const double log_det_K = 2.0 * Eigen::LLT<Eigen::MatrixXd>(K).matrixL().toDenseMatrix().diagonal().array().log().sum();

参考

既に実行中のジョブを走らせたまま端末を閉じる際の手順メモ

計算サーバ等で既に実行し始めてしまった時間のかかるプロセスについて、走らせ続けたまま端末を閉じたいことがある。

実行し始める前のプロセスについては、 nohup を使った方法がある。

ここでは、通常のコマンドで実行し始めてしまったプロセスについて考える。

disown を使う

時間がかかるタスクを通常のコマンドで実行し始めてしまったとする:

$ ./heavy_task.sh

対象のジョブを一旦停止する:

^Z

対象のジョブをバックグラウンド実行する:

$ bg

ジョブ番号を確認する:

$ jobs

disown コマンドによってバックグラウンド実行中のジョブをジョブテーブルから除外する:

$ disown %1

(これは対象のジョブ番号が1だったときの入力例)

以上の操作によって、端末を閉じてもプロセスは残ったままになる。

参考: