DIIS
direct inversion in the iterative subspaceの略でDIIS。$x$を変数とする関数/汎関数 $f(x)$/$f[x]$ を反復的に最適化する際に、近似解列 $x_i, i=1,\dots,N$ の重み付き線形結合で、$N+1$ 反復の近似解 $x_N$ を求める手法。このとき、解 $\bar{x}$ と $N+1$ 反復の近似解 $x_{N+1}$ のerror $e_{N+1} = \bar{x}-x_{N+1}$ を $x$ が満たすべき制限の下で最小化することに特徴がある。
元はPulayによる提案(P. Pulay, Chem. Phys. Lett. 73 (1980) 393)で、主にKohn-Sham方程式、Hartree-Fock方程式の自己無撞着解を求める為に使われている。前者では密度の最適化で、後者は密度行列の最適化に使われている。また、平面波や実空間グリッドのような基底で一体ハミルトニアンを評価した際はその対角化も疎行列を念頭に置いた反復的に求めるため、この求解にも使われている。ここではKohn-Sham方程式に対する密度の最適化について考える。
多分、物理や化学で現れる第一原理計算に限らず、一般の最適化問題にも適用できると思うが、この分野以外での適用例は検索してもあまり見当たらない。
密度行列の最適化、対角化へのDIISの適用に対するコメント
DIISをキーワードについて色々探していた際に、様々な文献で一見異なる問題を解いているように思えて混乱したので、少しコメント。
とにかく共通していることは、上で述べたように $N+1$ 反復での近侍解を $i\leq N$ の近似解列の重み付き線形結合で近似して、 $N+$ 反復の近似解列に対するerrorを最小化することである。errorを評価するには解を何かしらの方法で評価する(もちろんこの解はこの手続きで求めたいと思っているので、何かしらの方法で代替物を用意する)必要があり、このerrorを評価する際に流儀が二つある。
一つの方法は解を持っている情報の中から適当に推測する。下に書いたように密度汎関数理論のKohn-Sham法では、入力した密度に対してKohn-Sham方程式を解くことで、出力密度を得る。この出力密度をとりあえずの解候補としてerrorが評価できる。
もう一つの方法は以下では述べないので、少し詳しく説明する。この方法ではNewton-Raphson法で解を更新することを念頭においている：$x_{N+1}=x_N - \left(H_N\right)^{-1}g_N$。$H_N$, $g_N$ は$x_N$で評価したへシアンと勾配。$x_{N+1}$ が解に近いとみなすと、errorは $- \left(H_N\right)^{-1}g_N$ になる。へシアンを求めるのは通常困難なので、このへシアンを現実的な計算コストでうまく近似する。
この二つの方法はDIISの特徴(近似解列の重み付き線形結合で評価したerrorを拘束条件の下最適化)を備えているが、errorの評価方法が実際上全く異なっている。
準備：拘束条件付き最適化問題: Lagrangeの未定乗数法
$A$を$n$次の実対称行列とするときの、$\sum_{ij}^n c_i A_{ij} c_j$の値を$\sum_i^n c_i = 1$の下で最小化する$c_i$を求める。(拘束条件が$\sum_i \left|c_i\right|^2=1$であれば、行列$A$の最低固有値問題を求めることに対応するが、本問題とは設定が異なる。)
未定乗数$\lambda$を導入して
$$
L(c_i) = \frac{1}{2}\sum_{ij}^n c_i A_{ij} c_j - \lambda \left(1 - \sum_i c_i \right)
$$
以下の関数の停留条件を解く:
$$
\frac{\partial L}{\partial c_i} = 0.
$$
この条件は$$ B_{ij} = A_{ij} (1 \leq i,j \leq n),\ -1 (i=n+1, 1 \leq j \leq n), \ -1 (1 \leq i \leq n, j=n+1), \ 0 (i=j=n+1) $$$$ d_{i} = c_{i} (1 \leq i \leq n),\ \lambda (i=n+1),\quad b_{i} = 0 (1 \leq i \leq n), \ -1 (i=n+1)$$として未定乗数も含めた変数に対する線型方程式$$\sum_j^{n+1} B_{ij} d_{j} = b_{i}$$を考える。こ)の線型方程式を解いて、$c_i$をえる。
$A$が対称行列であれば、$B$も対称行列である。($n+1$行目の定義のしようによっては対称行列にならないようにもなるので、注意)
$B$が対称行列で、右辺ベクトルが実数なので、解$c_i, \lambda$は実数になる。
DFTでの実装
(書きかけ)
SCF収束解を与える密度を$\bar{n}$として、ある密度$n$のerrorを$e[\bar{n},n]=\bar{n}-n$とする。反復計算$i=1,2,\dots,N$での入力(出力密度は後で出てくる)密度$n_i^{\rm in}$における収束解とのerrorを$e_i$とする：$e_i = e[\bar{n},n_i^{\rm in}]$。
$N+1$反復で入力として使う密度$n_{N+1}^{\rm in}$をこれまでの入力密度のの重み付した和で近似する事を考える：
$$
n_{N+1}^{\rm in} = \sum_i^N c_i n_i^{\rm in}, \quad \sum_i^N c_i = 1.
$$
ここで導入された係数$c_i$は$N+1$反復でのerror $\left\langle \left. e[\bar{n},n_{N+1}^{\rm in}] \right| e[\bar{n},n_{N+1}^{\rm in}]\right\rangle$が最小になるように決める。このerrorを係数と反復$i$での入力密度でのerrorで陽に書くと、
$$
\left\langle \left. e[\bar{n},n_{N+1}^{\rm in}] \right| e[\bar{n},n_{N+1}^{\rm in}]\right\rangle = \sum_{ij} c_i^* A_{ij} c_j ,\qquad A_{ij} = \left\langle \left. e_i \right| e_j\right\rangle.
$$
この値を拘束条件付きで最小化すればよいので、前節の方法で解けることが分かった。
SCF収束解を与える密度を$\bar{n}$どう求めたら良いか という問題を一つ残している。
一つの解決策は$\bar{n}$として$n_i^{\rm in}$で構成されるハミルトニアンを対角化して得られる出力密度$n_i^{\rm out}$を使うことである。十分にSCFが収束してくれば入力と出力は近づくことが期待できる。