研究のメモとか
小ネタをメモします.
目次
eigenで特異値分解をする話 (2025/11/13)
逆行列を近似で求める話 (2025/11/14)
acpicaのメモ(2026/01/11)
計算化学徒のための数学物理(2026/01/23)
- 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追記)友人とこの近似について話していて, この近似はどのような行列で成立するのか, もう少し直観的に理解できないかと聞かれた.
一つ言えるのは, この近似は最終的にある行列の逆行列が(特定の条件の下で)対角成分の逆数を並べたものに近似できるといっているので,
元の行列の非対角成分の絶対値が対角成分に比べて十分に小さい(無視できるレベル)ということである. また仮に条件が成り立っていても,
いわばまさにテイラー展開の初項のみを取り出しているので, あまり非対角成分が大きい行列だと精度が急激に悪化することが推測される.
- acpicaのメモ(2026/01/11)
- acpicaのビルドで詰まったので書き残し.
- 環境
- Cコンパイラ:GCC14.2.0
- bison:3.8
- flex:2.6.4
- acpica:2025年12月のバージョン
- どうやらオブジェクトの重複宣言がされているらしい.
- 以下の手順を順に実行することで(無理やり)makeを通せる.
- まず/acpica/source/aslcompiler.hにあるAslCompilerlexの宣言をコメントアウト.
- この状態でmake. おそらく怒られが入って止まる.
- 続いて/acpica/generate/unix/iasl/obj/aslcompiler.y.hにある
#if !defined AslCompilererror && !defined YYERROR_IS_DECLARED
void AslCompilererror (const char *msg);
#endif
の部分を丸まるコメントアウト. 再びmake(ここのmakeはなくてもいいかも).
- また怒られて止まる. どうやらbythonのライブラリヘッダとacpicaのソース内のヘッダで型が食い違ってるものがありますよというエラー.
- aslcompiler.hに戻り, プロトタイプ宣言AslCompilererrorの返り値の型をvoidに変更. また同ディレクトリのaslerror.cのAslCompilererrorの返り値の型もvoidに変更.
- 再度make. しばらく順調にコンパイルが進むがまた怒られる. 先ほどintからvoidに変えたせいで返り値を返してる部分がおかしくなってますよというまっとうな指摘.
- 先ほどintからvoidに変えた部分をもう一回intに戻してmake. これでコンパイルが走りきるはずである.
- 計算化学徒のための数学物理(2026/01/23)
- 計算化学という我らが過疎分野に進む人は化学以外の内容をどう勉強すればいいんだろうか.
- 何らかの手違いで計算化学に足を踏み入れそうな, あるいは踏み入れてしまったB3,4を対象に考える.
- 数学に関しては行列・微積分, 物理は力学・電磁気学を教養で修めているものと仮定.
- まず行列の復習.「行列の指数関数」「特異値分解」「擬逆行列」みたいな重要だがあんま教養の授業で触れられないトピックに触れる.
行列の指数関数については『一歩進んだ物理数学 －レクチャーのその先へ－』(裳華房), その他は『これからの線形代数―３重対角化、特異値分解、一般逆行列』(紀伊国屋書店)が詳しい.
- 続いて初歩的な物理数学・統計学に触れる. ベクトル解析とかフーリエ・ラプラス変換とかね. 物理数学は『基礎解析学』(裳華房), 統計学は『現代数理統計学の基礎』(共立出版)で事足りる.
- 並行して物理もできるといい. まず解析力学に手を出す. 『解析力学 ―基礎の基礎から発展的なトピックまで―』(共立出版)がいいのではないでしょうか. ポがB3のころはまだ出てなくてなんか別の本を読んだ気がするけど忘れた.
- さらに並行してプログラミング. まずはLinuxの使い方を学びましょう. VMwareにubuntuでも突っ込んで遊んでみましょう. あと何かしら言語が1つ読めたほうがこの先参考書を読む際に有利です.
おすすめはC++かPython. Fortranは古いのでパス. これ以外の言語が計算化学で出てくることはあんまないと思います.
C++入門としてはAtCoderのAPG4bがあげられますが, 計算化学でプログラムを読み書きするなら避けては通れないポインタ, クラス,
stdに関する話題がちょっとしょぼいので, できれば本を買うといいと思います. Pythonは良く知りません.
- 物理の続き. 解析力学が終わったら量子力学か統計力学に手を出しましょう. どっちが重要かは分野によります. 第一原理計算とかやるなら量子力学, MD計算とかやるなら統計力学が重要だと思います.
- 多分このへんまで終わればあとは何とかなると思います.
- ここからは追加ステージ. 暇だったらまず集合論・位相空間論に手を出します. ハウスドルフ空間とかそういうやつです. これ自体は何の役にも立ちません.
- 続いて群論をやります. おすすめは『代数学1 群論入門』(日本評論社)などです. 群で分子構造を整理する, みたいなことがわかってきます.
- さらにホモロジー計算とかにも手を出してみましょう. 分子構造を代数トポロジー的な計算で処理するみたいな研究があったりします.
- また実解析に手を出してモンテカルロシミュレーションとかで遊ぶ道もあります. まずルベーグ積分と関数解析. 私は昔『新版 ルベーグ積分と関数解析』(朝倉書店)を読みました. ここにきてようやく無限次元のベクトルが扱えます.
- 確率論にも手を出します. 『例と演習で学ぶ 確率論』(講談社)などを読んで演習するとなんとなく公理的確率論の実感が湧いてきます.
- すると実解析が絡むようなモンテカルロ法や機械学習の本が読めるようになってきます. 良かったですね.
- ただこんなオプショナルな話題まで足を延ばすくらいなら, 『分子動力学法と原子間ポテンシャル』(森北出版)や『第一原理計算の基礎と応用 ―計算物質科学への誘い―』(共立出版)とか読んだほうが191919419倍ためになると思います.
もどる