研究のメモとか
小ネタをメモします.
目次
eigenで特異値分解をする話 (2025/11/13)
逆行列を近似で求める話 (2025/11/14)
C++で統計処理をしたい話 (2025/12/06)
- eigenで特異値分解をする話 (2025/11/13)
- C++で数値計算をする人ならお世話になったことがない人はいないであろうeigenについて.
- eigenで特異値分解をするとき, よくやりそうな操作である割に意外と日本語の情報が少なかったのでメモしておこうと思う.
- そもそも特異値分解とは何か...(ﾈｯﾄﾘ). 教員によっては教養レベルの線形代数で触れなかったりする. 実際私が学部1年生の時の線形代数の講義では登場しなかった.
- 詳しい説明は数学のサイトに譲るとして, 簡単に言えば対角化の一般化である.
- 特に, 扱いやすい直交行列で対角化されるのは実対称行列のみ, これはいけない. 正方行列でなくても, 実対称行列でなくても, このような分解が存在してほしいところである.
- 要は正方行列とは限らない行列Aを, USV*(U, V*は直交行列, Sは対角行列と零行列の組み合わせ, このSの非零成分を特異値と呼称)という形に分解できるという話である.
- この分解は有用なので非常によく使われており, 当然eigenにも機能として存在するが, 少し使いにくい.
- 公式ドキュメントはこちらにあるので,
ドキュメントを読むのに慣れてるなら最初っからこれを見ろという話なのだが, 一応メモ程度に情報を残しておこうと思う.
- まずeigenでの固有値分解には二種類のメソッドが提供されており, 一つはJacobiSVD,
もう一つはBDCSVDである.
前者は後者に比べ正確性で勝るが, 速度で劣るようである. ドキュメントでは, サイズが16次よりも小さい場合は前者を推奨している. 今回はJacobiSVDについて説明する.
- まず利用するためにはヘッダをインクルードする必要があるが, eigenのおまじないとして有名なDenseの中に含まれているようである.
- さて, 先ほど挙げたドキュメントのサンプルコードを見ると, JacobiSVDの関数本体はテンプレートとして実装されていることがわかる.
テンプレートは二つの引数をとっていて, 一つは分解する行列の型(サンプルコードではMatrixXfであるが, たいていの場合MatrixXdだろう),
もう一つは計算オプションの指定である. "ComputeThinU | ComputeThinV"という変な形の指定がされているが, これは左右の直交行列も計算してねというおまじないである.
どうやらこのおまじないを入れないと特異値しか計算してくれないようである. また別途計算アルゴリズムも指定できるようであるが, これはデフォのやつでいいんじゃないかと思う.
- 計算結果の取り出しも少しわかりにくい. JacobiSVDは分解する行列を引数として特異値分解した結果をオブジェクトとして返すコンストラクタである(にわかだから用語おかしかったらごめんね)ので,
値の取り出しはオブジェクトのメンバ関数を利用することになる. singularValuesメンバ関数は特異値を取り出すメンバ関数であるが,
どうやらベクトルの形で取り出しているらしく, 行列Sに変換するには(多分何らかの関数があるのかもしれないが)少し手を加えてやる必要がある.
- 残りの二つはmatrixU,matrixVで取り出す. これらは普通に対応する行列の型で取り出される. Vのほうを転置してUSV*を計算すると元の行列に戻るはずである.
- 逆行列を近似で求める話 (2025/11/14)
- 前の記事と続けて書いているのだが, 書いてる間に日付が変わってしまった.
- そこそこサイズの大きい正則行列Aがあるとしよう. 逆行列を求めるとき, 厳密に求めると計算量はO(N^3)程度なので,
計算化学としては避けたい壁の一つである. 一般に我々の界隈には次のような格言がある. 「逆行列を求めようとしたら負け」.
- ただやっぱり逆行列の値が欲しいケースというのは, ある. 実際, 先日私はそのような状況に陥っていた. しかも相手はかなりサイズの大きい行列である.
- こういう時, 何らかの近似をぶっこむというのは生きていく上で欠かせない. 自然科学とはつまるところ近似であり, 我々は近似で飯を食っているのである.
今回私が欲しかったのは逆行列の対角成分の値のみであったので, 次のような手法で近似を行った.
- さて, 件の行列Aの対角成分のみを取り出した行列をB, それ以外の成分を取り出した行列をCとする. このとき, \[A^{-1}=(B+C)^{-1}\]
Bの対角成分がすべて非零であることを仮定すると, Bは正則である. このとき, 同サイズの単位行列をIとして, \[=(B(I+B^{-1}C))^{-1}\]
さらに\[=(I+B^{-1}C)^{-1}B^{-1}\]であり, これを\[=(I-(-B^{-1}C))^{-1}B^{-1}\]とする.
さて, ここで\[(I-(-B^{-1}C))^{-1}\]の部分を無限和として\[I+\sum_{n=1}^{\infty}(-1)^n(B^{-1}C)^n\]と書き換える.
(1-x)^{-1} = 1+x+x^2+...のノリである. 実はこの変形は「ノイマン級数」として知られており, 有界線形作用素においてはある条件の下で正当化される.
行列はフロべニウスノルムを取ればわかる通り, 有界線形作用素である.
その条件とやらは置いておき, ひとまず話を進める. このとき\[\left(I+\sum_{n=1}^{\infty}(-1)^n(B^{-1}C)^n\right)B^{-1}\]の
第二項以降を落とすと, \[A^{-1} \simeq B^{-1}\]と近似できる. B^{-1}は対角行列なので逆行列はただ各値の逆数をとるだけであり, 計算量はO(N)である.
- さて, 途中でごまかした条件というのは「作用素ノルムが1未満であること」である. 行列の積から誘導される作用素ノルムはスペクトルノルムとして知られ,
今回この近似が成り立つためには, \[(-B^{-1}C) = D\]として, 転置行列との積D*Dの最大固有値の正の平方根が1未満であることが成立の条件である.
- 今回, この近似はかなりうまくいき, 短い計算時間で実験値とよく一致する値を得ることができた. しかし残念ながら先の条件については, 実は未検証である.
いつかこの近似が成立しない行列を投げつけられるんじゃないかと, ひやひやする日々である.
- (20251116追記)友人とこの近似について話していて, この近似はどのような行列で成立するのか, もう少し直観的に理解できないかと聞かれた.
一つ言えるのは, この近似は最終的にある行列の逆行列が(特定の条件の下で)対角成分の逆数を並べたものに近似できるといっているので,
元の行列の非対角成分の絶対値が対角成分に比べて十分に小さい(無視できるレベル)ということである. また仮に条件が成り立っていても,
いわばまさにテイラー展開の初項のみを取り出しているので, あまり非対角成分が大きい行列だと精度が急激に悪化することが推測される.
- C++で統計処理をしたい話 (2025/12/06)
- いろいろあって, 研究で少し統計的な処理をしたいという話になった.
- ところが統計処理といえばPythonとRの天下, C++で統計処理なんて誰もやらない.
- しかし私は学部時代に授業でRでの統計処理を行い, 奇怪な文法で発狂した経験がある. 今更Rstudioをインストールするのも面倒.
なんとかC++で統計処理を行えないものか.
- 調べたところ, いくつかC++で統計処理を行うライブラリがあった. どうやらalglibというライブラリが一番機能的に充実していそうだったので,
これを採用.
- しかしalglibは日本語の情報が少ないので, インストールの方法だけでもひとまずここにメモしようと思う.
- cppソースとヘッダファイルの形式で配布されており, まずこれを落としてくる. そんで展開してワーキングディレクトリに全部おいて,
ヘッダファイルを都度インクルードしつつコンパイルは全部のcppソースをまとめてやってねというのが公式の案内となっている.
- でもこれだとディレクトリ内がゴチャゴチャしてなかなか厳しい. ここはライブラリにしてしまうのが得策だろう.
まず落としてきたファイルを展開した中身のcppファイルを-fPICオプションをつけてまとめて全部オブジェクトファイル(.o)にコンパイルする.
- ここからはちょっとOSによって異なる.
- Linuxの場合
- g++に-sharedオプションをつけて先ほどのオブジェクトファイルをすべてリンク. ファイル名は「libalg.so」とでもする.
- ヘッダはinclude先のディレクトリに置き, libalg.soは/usr/lib/x86.../に配置. この場所は環境や設定で変わるのでダメだったらコンパイラの設定見て.
- あとは計算したいプログラムを書き, コンパイル時にオプション-lalgを指定してコンパイル. undefined referenceエラーが出る場合は件のライブラリファイルの場所を見直したり,
ライブラリのパスが通っているかを確認する.
- Windowsの場合(MinGWがインストールされていることを想定)
- arコマンドでアーカイブ化. ファイル名は「libalg.a」
- ヘッダはinclude先のディレクトリに置き, libalg.aは(前略)\mingw64\x86_64-w64-mingw32\libに配置(コンパイル時に指定するだけなので, 場所にこだわりがあれば別の場所でも可).
- あとは計算したいプログラムを書き, コンパイル時にオプション-Iオプションで.aファイルが置いてあるディレクトリを指定し, -lalgを指定してコンパイル. undefined referenceエラーが出る場合はライブラリのパスを打ち間違えてないか確認する.
でもたぶん素直にPythonで書いたほうがいい.
もどる