# 密度フィッティングと恒等分解(RI)

2 電子積分テンソル $(\mu\nu|\lambda\sigma)$ は $O(n_{\text{ao}}^4)$ にスケールします — ほとんどの
量子化学計算における単一最大のコストです。**密度フィッティング**(**恒等分解(RI)**とも呼ばれる)は、
すべての 4 添字積分を、はるかに小さな 3 添字因子にわたる縮約で近似し、$O(n_{\text{ao}}^4)$ の対象を
制御可能な小誤差で $O(n_{\text{fit}}\cdot n_{\text{ao}}^2)$ のものに変えます。本章はこの近似自体、
qc-rs があらゆる RI バックエンドで共有する「因子 $B$」に変換するために使う厳密な計量の白色化、そして
その単一因子が平均場 J/K とポスト HF 相関の両方にどう奉仕するかを導出します。
`.design/54two-electron-integral-spine.md` と `.design/68ri-jk-strategies.md` に基づき、`AGENTS.md`
の統一 `eri=` RI 戦略記述と一致させます。

## 恒等分解近似

厳密な AO 積密度 $\rho_{\mu\nu}(\mathbf r)=\chi_\mu(\mathbf r)\chi_\nu(\mathbf r)$ は一般にどんな有限
補助基底でも厳密に展開できませんが、フィッティング関数 $\{\varphi_P\}$ の線形結合で**近似**できます:

$$
\rho_{\mu\nu}(\mathbf r) \approx \tilde\rho_{\mu\nu}(\mathbf r) = \sum_P c^{\mu\nu}_P\,\varphi_P(\mathbf r).
$$

Dunlap の変分 RI は、フィッティング係数 $c^{\mu\nu}_P$ を、残差 $\rho_{\mu\nu}-\tilde\rho_{\mu\nu}$ の
実空間ノルムではなく**Coulomb 自己エネルギー**を最小化するように選びます — $J$ と $K$ の精度を決める
のは点ごとの密度誤差ではなく Coulomb エネルギーだからです。
$\iint\frac{(\rho-\tilde\rho)(\mathbf r)(\rho-\tilde\rho)(\mathbf r')}{|\mathbf r-\mathbf r'|}\,d\mathbf
r\,d\mathbf r'$ を $c^{\mu\nu}_P$ について最小化すると、停留条件
$\sum_Q(P|Q)c^{\mu\nu}_Q=(P|\mu\nu)$ が得られます、すなわち

$$
c^{\mu\nu}_P = \sum_Q (V^{-1})_{PQ}\,(Q|\mu\nu), \qquad V_{PQ} = (P|Q) = \iint
\frac{\varphi_P(\mathbf r)\varphi_Q(\mathbf r')}{|\mathbf r-\mathbf r'|}\,d\mathbf r\,d\mathbf r',
$$

$(Q|\mu\nu)$ は 1 つの補助関数と 1 つの AO 対の間の 3 中心 Coulomb 積分、$V$ は補助基底だけにわたる
2 中心 Coulomb 計量です。このフィットをどんな 4 添字積分にも代入すると、あらゆる RI バックエンドが
実際に使う近似が得られます:

$$
(\mu\nu|\lambda\sigma) \approx \sum_{PQ}(\mu\nu|P)\,(V^{-1})_{PQ}\,(Q|\lambda\sigma).
$$

これは(場当たり的な打ち切りではなく)真の変分最小化から来るので、RI 誤差はなめらかで、補助基底を
拡大することで系統的に改善可能であり、適切に対応するフィッティング基底(`cc-pvdz-jkfit`、
`def2-svp-jkfit` など、軌道基底ごとに表になっている)では経験的に微小です — これが、RI が単にメモリ
逼迫時に頼るだけの安い近似ではなく、qc-rs の*推奨既定*平均場経路である理由です。

## 白色化:計量を単一の 3 添字因子にまとめる

$(V^{-1})_{PQ}$ を各 J/K 縮約の中で直接扱うと、Fock 構築ごとに $n_{\text{fit}}\times n_{\text{fit}}$
演算子を再適用することになります。qc-rs は代わりに、3 中心積分を一度だけ**白色化**し、計量全体を
単一の 3 添字因子に折り込みます:

$$
B^P_{\mu\nu} = \sum_Q\bigl(V^{-1/2}\bigr)_{PQ}\,(Q|\mu\nu), \qquad
(\mu\nu|\lambda\sigma) \approx \sum_P B^P_{\mu\nu}\,B^P_{\lambda\sigma},
$$

これは $V^{-1}=V^{-1/2}V^{-1/2}$ を使って $\sum_{PQ}(\mu\nu|P)(V^{-1})_{PQ}(Q|\lambda\sigma)$ を書き
換えたものにほかなりません — 非対称な分解ではなく対称な平方根分解を選ぶことで、$B$ が計量全体を運び、
下流のすべての縮約が $P$ にわたる*単純な*内積になり、ループの中に計量が残らなくなります。qc-rs は
$V^{-1/2}$ を 2 通りの等価な方法で計算します:$V$ の固有分解(ゼロに近い固有値を `ri.metric_cutoff`
IOP 未満で落とす — 大きい/拡散した補助基底では計量が軽度にランク落ちしうる)、あるいは Cholesky 解
$V=LL^{\mathsf T}$、$B=L^{-1}\cdot(\cdot)$(固有分解なし、数値的により安定、その場で実行)— どちら
も同じ演算子について*同じ*縮約済み $(\mu\nu|\lambda\sigma)$ を与えるので、どちらを使うかは物理の
選択ではなく数値安定性の選択です。

検証済み例 — 水での RI 近似誤差、そして 2 つの白色化実装の厳密な一致(`ri-ram` = B を保存、
`ri-recomp` = サイクルごとに B を再計算、どちらも同じ数式から構築):

```python
import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

m4c = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").ints(eri="4c-incore").scf(ref="r").run()
m4c.scf.energy   # -76.02679364497408 -- 厳密な 4 中心参照

mri = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").ints(eri="ri-ram").scf(ref="r").run()
mri.scf.energy   # -76.02677274798874
mri.scf.energy - m4c.scf.energy   # 2.09e-05 Ha -- RI フィッティング誤差、自動導出された jkfit で

mri2 = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").ints(eri="ri-recomp").scf(ref="r").run()
mri2.scf.energy - mri.scf.energy   # -1.4e-14 -- ri-ram と ri-recomp は浮動小数点簡約の底まで一致
```

厳密な 4 中心参照に対する $2\times10^{-5}\ E_h$ のフィッティング誤差は、手動調整のない既定補助基底
としては典型的な RI-JK の値であり — 背後の軌道基底自体の精度をはるかに下回ります。明示的な `rijk=`
補助基底が与えられない場合、qc-rs は既定の JK フィッティング基底を未設定のままにするのではなく AO
基底から自動導出します(`default_jk_aux`)— 要求された厳密な方法を静かに RI へ「昇格」させる
「auto-RI」は存在しませんが、RI が明示的に要求されればいつでも妥当な補助基底の既定値があります。

## 白色化因子からの $J$ と $K$

$B$ が与えられれば、Coulomb と交換の両方は密度行列との単純な縮約になり、$B$ が常駐していれば
(RAM でもディスク/RMA からストリームされていても — ストレージバックエンドはポリシーであり、この
数式を変えません)完全にローカルです:

$$
J_{\mu\nu} = \sum_P B^P_{\mu\nu}\,q_P, \qquad q_P = \sum_{\lambda\sigma}B^P_{\lambda\sigma}D_{\lambda\sigma},
$$

$$
K_{\mu\nu} = \sum_P\sum_i B^P_{\mu i}\,B^P_{\nu i}, \qquad B^P_{\mu i}=\sum_\nu B^P_{\mu\nu}C_{\nu i}
\ \text{(占有 MO へ半変換済み)}.
$$

$J$ は小さな補助空間ベクトル $q_P$ だけを必要とします(密度に対する縮約 1 回、$B$ に対する縮約もう
1 回);$K$ はまず $B$ を占有 MO 添字へ半変換する必要があります([次章](post-hf-correlation.md)で
相関法が使う AO→MO 半変換とまったく類似)、その後外積の和です。どちらも完全な 4 添字テンソルには
決して触れません。これが、大規模での通常の 4 中心交換より RI が $K$ について著しく安い理由でもあり
ます — 半変換は外積の前に片方の AO 添字を $n_{\text{occ}}$ に縮めます、$n_{\text{ao}}^2$ を直接
$n_{\text{ao}}^2$ と縮約するのではなく。

## $B$ を得る 2 つの方法:再計算対保存

qc-rs の本番 RI 機構は、`eri=` の値で選ばれる**1 つの μ 分散機構に 2 つの $L$ ソースモード**
(`MuRiJk`/`MuRiJkDistributed`)です:

- **`ri-recomp`**は何も保存しません:*SCF サイクルごとに* libcint から 3 中心 $L=(Q|\mu\nu)$ を
  再計算し、常駐する白色化 $W=V^{-1/2}$ とその場で縮約します — 完全な $B$ テンソルを一切材料化
  しません、CPU(サイクルごとの再評価)とメモリ(小さな計量因子だけが常駐)をトレードします。
  `ri-recomp-disk` はさらに厳しいメモリ予算向けのアウトオブコア兄弟であり、サイクルごとの中間物を
  チェックポイントの `scratch` ディレクトリにスプールします。
- **`ri-ram`**は白色化因子 $\tilde B=W\cdot L$ を**一度だけ**その場で構築し、その後のすべての SCF
  サイクルで常駐させ続けます — サイクルごとの $J$ と $K$ の両方が*同じ*常駐 $\tilde B$ を読み、
  サイクルごとの 3 中心再計算はまったくありません。これは、完全な因子をメモリに保持するコストで
  サイクルあたりより高速です。

どちらも数学的に同一です — 上の検証済み例が直接示す通り浮動小数点簡約順序まで互いにビット一致します
(`ri-ram` と `ri-recomp` は $\sim10^{-14}$ だけ異なり、これは通常の浮動小数点ノイズの底であって
物理的な違いではありません)— 両者の選択は純粋なメモリ/CPU のトレードオフであり、精度のトレード
オフではありません。`ri-ram` は大規模での incore 並みの速度の既定であり、`ri-recomp` は最大の基底の
ためのメモリ節約の選択であり、SCF 自体が保存したか再計算したかに関わらず相関法が同じ 3 中心構成要素を
必要とするため `MoTransformSource` が直接再利用するものでもあります。

## 分散計量なしで $B$ を分散する

分散メモリ RI を一般に難しくする唯一の要因は、素朴な方式では計量 $V^{-1/2}$ を、3 中心テンソルを保持
するどんな分散方式にも*またがって*適用する必要がある — SCF サイクルごとの全対全通信または遠隔メモリ
行列ベクトル積 — ことです。qc-rs の設計は、**$B$ が SCF 不変**(積分だけから一度構築され、すべての
SCF 反復にわたって変更なく再利用される)であることを利用してこれを完全に回避し、再分散コストは構築時
に一度だけ払われます:

1. **計量は分散配列ではなく複製配列です。** $V$ とその因子($V^{-1}$/$V^{-1/2}$/Cholesky 解演算子)は
   $O(n_{\text{fit}}^2)$ — $O(n_{\text{fit}}\cdot n_{\text{ao}}^2)$ の 3 中心テンソルに比べて微小 —
   なので、各ランクは単に自分の完全なコピーを持ちます([`Workspace` RMA プール](../30-hpc/threads-and-blas.md)
   の `dup` 配列の種類:`get`/`put` はローカルメモリコピーで MPI ゼロ)。これほど小さいものを複製
   するのは無料です;これがそれ以降のすべてのステップをローカルにするものです。
2. **白色化は一度だけ、ローカルに、1 回限りの転置の間に適用されます。** 各ランクの生の 3 中心テンソル
   のスラブは、複製された計量因子に対するローカル GEMM で白色化されます — *反復*SCF ループのどこにも
   分散計量行列ベクトル積はなく、この 1 回限りの構築パス中にのみ(せいぜい)存在します。
3. **白色化因子は補助添字 $P$ で分散され、ノードごとに完全な AO 対を持ちます。** 各ランクは $P$ の
   連続した行ブロックを所有し、自分の行については*完全な* $\mu\nu$ 範囲を保持します。これがまさに
   $J$(所有行のローカル GEMM を $D$ に対して行い、その後小さな $O(n_{\text{ao}}^2)$ の `Allreduce`)
   と $K$(所有行の半変換、その後ローカル外積)の両方を各ランクで完全にローカルにするものです —
   SCF サイクルごとのランク間通信は Fock 行列の 1 回の小さな `Allreduce` だけであり、分散テンソル
   縮約は決してありません。

メモリは線形に集約します:各ランクは $B$ の自分の $P$ 行分だけを保持するので、1 ノードの RAM に
大きすぎる $B$ テンソルを持つ分子も、数ノードにまたがれば快適に収まります — `AGENTS.md` の qc-hpc/RMA
規則が記述する、主要な分散メモリターゲットである≤4ノードのメモリ集約体制です(大規模 HPC スケール
アウトではありません)。マルチランク $K$ 構築自体にはさらなる選択、`ri.ram.kdist` IOP キーがあります:
`"p-transpose"`(既定 — 1 回限りの `MPI_Alltoallv` が $\tilde B$ をフィット添字分散レイアウトへ転置
するので、以降のすべての SCF サイクルの $J$/$K$ はランクローカルにプラス小さな
$\operatorname{Allreduce}(n_{\text{ao}}^2)$ が 1 回)または `"mu-stream"`(μ 分散、サイクルごとに
半変換を Gram ストリームする代わり)— サイクルごとの通信対 1 回限りの転置のトレードオフであり、
精度のトレードオフではありません。

```{mermaid}
flowchart TD
    L["生の 3 中心 L = (Q|mu nu)<br/>(libcint, 補助殻ごと)"] --> WHITEN["白色化: B = V^-1/2 . L<br/>(複製計量、ローカル GEMM)"]
    WHITEN --> DIST["B を補助添字 P で分散<br/>(ランクごとに完全な mu-nu)"]
    DIST --> J["J: q_P = sum B^P.D (所有行)<br/>J = sum_P B^P q_P<br/>+ Allreduce(nao^2)"]
    DIST --> K["K: B を占有 MO へ半変換<br/>K = sum_P sum_i B^P_mu-i B^P_nu-i<br/>+ Allreduce(nao^2)"]
    DIST --> MO["相関: AO->MO 半変換<br/>B^P_pq (任意の MO 空間対)"]
```

## 同じ因子が相関にも奉仕する

[RI-MP2](post-hf-correlation.md) が必要とする AO→MO 半変換は、密度行列の代わりに MO 係数と縮約する
だけで、*同じ*白色化因子から構築されます:

$$
B^P_{pq} = \sum_{\mu\nu}C^L_{\mu p}\,B^P_{\mu\nu}\,C^R_{\nu q}
$$

任意の左/右 MO 空間対について(RI-MP2 では占有×仮想ですが、同じ機構が将来のどんな相関法にも占有×占有
や仮想×仮想として奉仕します)。qc-rs の設計が意図的に利用する 1 つの機微:白色化と変換の*順序*は
入れ替え可能で、「先に白色化してから変換」(上記の J/K 形式)と「先に変換してから白色化」は*まったく
同じ*最終 $B^P_{pq}$ を与えます。白色化は補助添字だけへの線形演算であり、軌道添字への AO→MO 縮約と
可換だからです。相関は**変換してから白色化**を使います:まず $n_{\text{ao}}^2\to n_{\text{left}}\cdot
n_{\text{right}}$ に縮める(占有/仮想空間が $n_{\text{ao}}$ からかけ離れていればずっと小さな対象)、
*それから*計量を適用します — 完全な AO 対テンソルを先に白色化するより $5$〜$20$ 倍少ない白色化
FLOP であり、対応してより小さな分散簡約です。これが、[RI-MP2 の章](post-hf-correlation.md) の
$B^P_{ia}$ が別途導出される対象ではない理由です — それは本節の $B$ を、もう一度変換したものです。

:::{exercise}
:label: ex-ri-theory-ja

1. 検証済み例は、`ri-ram` と `ri-recomp` が $\sim10^{-14}$ で一致する一方、両者とも厳密な 4 中心
   から $\sim2\times10^{-5}$ 異なることを示しています。これらが完全に異なる種類の「誤差」— 1 つは
   真の物理近似、もう 1 つは浮動小数点ノイズ — である理由と、どちらか一方だけがより大きな補助基底で
   縮小する理由を一文で説明しなさい。
2. 分子や AO 基底がどれほど大きくなっても計量 $V$ を常に `dup`(完全複製、MPI ゼロ)配列として
   扱えるのに対して、3 中心テンソル $B$ にはそれができないのはなぜですか。
3. 相関法は AO→MO 変換の*後*に計量を適用し(「変換してから白色化」)、平均場 J/K はその*前*に適用
   します(「白色化してから変換」)。どちらも同じ演算の数学的に厳密な並べ替えです。相関が平均場とは
   逆の順序を好む実務上の(数学上ではない)理由は何ですか。
:::

:::{solution} ex-ri-theory-ja
:class: dropdown

1. `ri-ram` 対 `ri-recomp` の違い($\sim10^{-14}$)は、同じ厳密な和を異なる順序で(一度保存して
   再利用 対 サイクルごとに再計算)累積したことによる純粋な浮動小数点丸めであり — 機械イプシロンで
   有界であり、どちらの方法も同一の数学量を計算しているのでどんな基底変更でも縮小しません。
   `ri-ram`/`ri-recomp` 対 4 中心の違い($\sim2\times10^{-5}$)は、恒等分解を有限補助基底に打ち切る
   ことによる真の RI 近似誤差です — これは実在する物理的/数学的近似であり、補助基底を完全性へ向けて
   拡大すれば(系統的に、最終的には数値ノイズまで)*縮小します*。
2. 計量 $V=(P|Q)$ は $O(n_{\text{fit}}^2)$ にスケールします — *補助*基底サイズだけの 2 乗 — のに
   対して、3 中心テンソル $B$ は $O(n_{\text{fit}}\cdot n_{\text{ao}}^2)$ にスケールし、(典型的には
   はるかに大きく速く成長する)AO 基底サイズの 2 乗でもあります。現実的なサイズのどんな分子でも、
   計量は各ランクで安価に複製できるほど小さいままである一方、$B$ はまさに、そのサイズがそもそも
   ノードをまたぐメモリ集約を強いる対象です — それを分散する代わりに複製することは、分散する目的
   全体を台無しにします。
3. AO→MO 変換自体は、安価な計量適用よりやり直すのがはるかに高価であり、テンソルを
   $n_{\text{ao}}^2$ から $n_{\text{left}}\cdot n_{\text{right}}$ へ*縮小*するからです(小さな基底を
   過ぎれば占有×仮想は AO×AO よりはるかに小さい)。より小さな、既に変換済みの対象に計量を適用する
   ことは、より大きな AO 対の対象に適用するより厳密に安い GEMM であり、白色化中のどんな分散簡約も
   より小さなテンソルで動作することを意味します。対照的に平均場 J/K は SCF サイクルごとに*AO*基底の
   結果を必要とします(Fock 行列は AO の対象です)ので、計量適用を先送りする縮小変換が存在せず —
   AO レベルで前もって一度白色化することが、そこでは既に最も安い選択肢です。
:::

次の 2 つの章はこの因子の上に直接築かれます:[解析的微分](analytic-derivatives.md) と
[線形応答理論](linear-response-cphf.md) はどちらも、この $B$ と RI 整合な微分積分を必要とし、
[RI-MP2](post-hf-correlation.md) は上で示した通り、それをもう一度縮約したものにほかなりません。