カルーネン・レーベ展開の理論と導出 — 確率過程を最適な基底で表す

センサから流れてくる波形を1000点サンプリングすると、1本の観測が1000次元のベクトルになります。これを100本集めれば1000×100の数値の山です。ところが実際に眺めてみると、どの波形も「ゆっくりした上下」と「わずかな細かい揺れ」の重ね合わせでしかなく、本質的な自由度は10個もないように見える。ここで自然な問いが生まれます。この確率的な波形の集まりを、たった $N$ 本の関数の重ね合わせで表すとしたら、どんな関数を選ぶのが一番得なのでしょうか。

フーリエ級数を使えばよい、と考えるかもしれません。実際それでも表せます。しかしフーリエ基底は「その波形がどんな統計的性質を持つか」を一切見ずに決め打ちした基底です。データ自身の相関構造に合わせて基底をあつらえれば、同じ $N$ 本でもっと誤差を減らせるはずです。この「あつらえた最適な基底」を与えるのがカルーネン・レーベ展開（Karhunen–Loève expansion, KL展開）です。

KL展開は、主成分分析（PCA）を有限次元ベクトルから無限次元の関数へ持ち上げたものだと考えるとしっくりきます。実際、応用範囲もPCAと同じくらい広いです。たとえば、

• 不確実性定量化（UQ）と確率有限要素法 — 材料定数や境界条件が空間的にランダムに揺らぐとき、そのランダム場をKL展開して有限個の確率変数に落とし、モンテカルロやスペクトル法の計算量を劇的に減らします。
• 気象・海洋データのEOF解析 — 海面水温や気圧場の時空間データを共分散の固有関数（経験直交関数）で分解し、エルニーニョのような支配的な変動モードを取り出します。
• 信号圧縮とKLT — 画像・音声の圧縮で使われるKLT（Karhunen–Loève変換）は、離散版のKL展開そのものです。DCTはこれを高速に近似する道具として生まれました。
• 関数データ解析（fPCA） — 成長曲線や日内の電力需要カーブといった「1サンプル＝1本の曲線」のデータを、少数の主成分関数で要約します。

この記事では、「係数が無相関になるような直交展開を作りたい」という素朴な要求から出発し、それが第二種フレドホルム積分方程式という固有値問題に必然的に行き着くことを示します。そのうえで、打ち切り誤差がちょうど残りの固有値の和になること、そしてあらゆる正規直交系の中でこの基底が誤差最小であることをラグランジュ乗数法で証明します。最後に、ブラウン運動という具体例で固有関数と固有値を手計算で解き切り、その解析解が数値計算と一致することをPythonで確かめます。

本記事の内容

• 直交展開の係数が無相関になる条件から積分方程式が出てくる筋道
• マーサーの定理とKL展開の定理（証明つき）
• 打ち切り誤差 $\sum_{k>N}\lambda_k$ の導出と、その最適性の証明
• ブラウン運動のKL展開の解析解 $\sqrt{2}\sin((k-\tfrac12)\pi t)$ の完全な導出
• Pythonによる数値固有分解・再構成誤差の測定・他の基底との比較

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• ブラウン運動（ウィーナー過程）の理論を解説
• 固有値分解の理論と応用を完全解説する
• 共分散行列とは？定義・性質・固有値分解からPython実装まで完全ガイド

特にPCAを分散最大化と最小二乗誤差の両面から理解しておくと、この記事の「最適性」の議論がそのまま連続版の焼き直しであることに気づけます。まだの方は主成分分析（PCA）の理論と導出も先に眺めておくとよいでしょう。

カルーネン・レーベ展開とは — PCAを関数に持ち上げる

この図がKL展開の全体像です。左の①は「1本1本は違うが統計的な癖を共有するランダムな波形の束」で、これが出発点のデータにあたります。そこから2次統計量だけを抜き出したのが中央の②の共分散関数 $R(s,t)$ で、色が明るいほど「その2時刻が一緒に大きく動く」ことを表します。この $R$ の固有値問題を解いて得られるのが右の③の固有関数 $\varphi_k$ であり、$k$ が増えるほど振動が細かく、対応する固有値（そのモードの分散）は $0.4053 \to 0.0450 \to 0.0162 \to 0.0083$ と急速に小さくなります。以降の節では、なぜこの3ステップが「最適」なのかを一つずつ詰めていきます。

ギターの弦を弾くところを思い浮かべてください。弦の形は無限個の点の集まりですが、実際の振動は基本振動・2倍振動・3倍振動……という決まった「モード」の重ね合わせで書けます。しかも、強く出るのは低次のモードで、高次はどんどん小さくなる。だから最初の数モードだけ取れば、弦の形はほぼ再現できます。

KL展開がやりたいのは、これと同じことをランダムな波形に対して行うことです。違うのは、モードを決めるのが弦の物理（波動方程式）ではなく、データの統計的な相関構造（共分散関数） だという点です。「どの時刻とどの時刻が一緒に動きやすいか」という情報だけからモードを決めるので、対象がブラウン運動でも株価でも風速場でも、同じ手続きが使えます。

もう一つのイメージはPCAです。PCAでは、$d$ 次元のデータ点の雲に対して共分散行列 $\bm{\Sigma}$ の固有ベクトルを求め、分散の大きい方向から順に軸を並べ替えました。KL展開では、データ点が「点」ではなく「関数」になり、共分散行列 $\bm{\Sigma}$ が「共分散関数 $R(s,t)$」になり、固有ベクトルが「固有関数 $\varphi_k(t)$」になります。添字 $i=1,\dots,d$ が連続変数 $t\in[0,T]$ に置き換わっただけで、骨格は完全に同じです。

結論を先に書いてしまうと、KL展開とは次の表現のことです。

$$ \begin{equation} X(t) = m(t) + \sum_{k=1}^{\infty} \sqrt{\lambda_k}\,\xi_k\, \varphi_k(t) \end{equation} $$

ここで $m(t)=E[X(t)]$ は平均関数、$\varphi_k$ は正規直交な固有関数、$\lambda_k$ はそれに対応する固有値、$\xi_k$ は平均0・分散1・互いに無相関な確率変数です。ランダムさは全部 $\xi_k$ に押し込められ、時間依存性は全部 $\varphi_k(t)$ が引き受けます。この「ランダムさと時間依存性の分離」がKL展開の最大の売りで、たとえばランダム場を有限個の確率変数でパラメータ化したいUQの応用では、まさにこの形が欲しかったわけです。

では、なぜ固有関数でなければならないのでしょうか。次の節から順に、必要な道具を用意したうえで、その必然性を導きます。

準備：2次確率過程と共分散関数

まず、扱う対象をきちんと決めましょう。区間 $[0,T]$ 上で定義された確率過程 $\{X(t)\}_{t\in[0,T]}$ を考えます。全ての $t$ について2次モーメントが有限、すなわち

$$ E[X(t)^2] < \infty \quad (\forall t \in [0,T]) $$

を満たすものを2次確率過程と呼びます。この仮定があると、$X(t)$ たちを二乗可積分な確率変数のヒルベルト空間 $L^2(\Omega)$ の元として扱え、「内積＝共分散」という幾何が使えるようになります。

議論を見やすくするため、以降は平均を引いて中心化してあるとします。つまり $E[X(t)] = 0$ とします。平均が0でない場合は、最後に $m(t)$ を足し戻すだけなので一般性は失われません。中心化した過程の共分散関数は

$$ \begin{equation} R(s,t) = E[X(s)X(t)], \qquad s,t \in [0,T] \end{equation} $$

で定義されます。この $R$ は共分散行列の連続版なので、共分散行列と同じ2つの性質を持ちます。

性質1（対称性）：$R(s,t) = E[X(s)X(t)] = E[X(t)X(s)] = R(t,s)$。積の可換性から明らかです。

性質2（非負定値性）：任意の実数 $c_1,\dots,c_n$ と時刻 $t_1,\dots,t_n$ に対して

$$ \sum_{i=1}^n\sum_{j=1}^n c_i c_j R(t_i,t_j) = E\left[\left(\sum_{i=1}^n c_i X(t_i)\right)^2\right] \ge 0 $$

が成り立ちます。左辺が二乗の期待値に化けるのは、$R(t_i,t_j)=E[X(t_i)X(t_j)]$ を代入して和と期待値を交換し、二重和を「和の二乗」にまとめ直したからです。二乗の期待値は負になりえないので、非負定値性が言えます。これは共分散行列が半正定値であることの、そのままの連続版です。

これから主役になるブラウン運動を例に、共分散関数の姿を先に見ておきます。左は標本路30本で、原点から出発して $\pm\sqrt{t}$ の破線に沿って扇形に広がっていく様子がわかります。中央のヒートマップが共分散関数 $R(s,t)=\min(s,t)$ で、対角線 $s=t$ に沿った値がそのまま各時刻の分散 $R(t,t)=t$ になり、右上（遅い時刻どうし）ほど明るくなっています。右は $s$ を固定した断面で、$ts$ になると $s$ の値で頭打ちになる折れ線です。この折れ点が、後でブラウン運動の積分方程式を解くときの場合分けの正体になります。

さらに $R$ が $[0,T]^2$ 上で連続だと仮定します。この連続性は後で使うマーサーの定理に必要です。実用上出てくる過程の多くはこれを満たします（ブラウン運動の $R(s,t)=\min(s,t)$ も連続です）。

なぜ $R$ だけあれば十分なのかという点も押さえておきましょう。KL展開が使うのは $X$ の2次統計量だけです。ですから、分布の細かい形（歪度・尖度）は結果に影響しません。逆に言えば、$X$ がガウス過程なら2次統計量が分布を完全に決めるので、KL展開は「情報を一切捨てていない」表現になります。ガウスでない場合でも展開自体は成立しますが、後で見るように係数の「無相関」が「独立」まで強くはならない、という違いが出ます。

道具が揃いました。次は、展開の形を仮定するところから始めて、基底がどう決まるかを追いかけます。

展開の形を仮定する — 何を要求すれば基底が決まるか

いきなり固有関数を持ち出すのではなく、「こういう展開が欲しい」という要求から出発します。$[0,T]$ 上の正規直交系 $\{\psi_k\}_{k=1}^\infty$、つまり

$$ \langle \psi_j, \psi_k \rangle = \int_0^T \psi_j(t)\psi_k(t)\,dt = \delta_{jk} $$

を満たす関数列を持ってきて、$X$ をこれで展開したとします。

$$ X(t) = \sum_{k=1}^{\infty} c_k \psi_k(t) $$

係数 $c_k$ は正規直交性から一意に決まります。両辺に $\psi_j(t)$ を掛けて $t$ で積分すると、右辺は $\sum_k c_k \delta_{jk} = c_j$ となるので、

$$ \begin{equation} c_j = \int_0^T X(t)\psi_j(t)\,dt = \langle X, \psi_j\rangle \end{equation} $$

です。$X$ がランダムなので、$c_j$ もまたランダムな量（確率変数）になります。ここまではフーリエ級数と何も変わりません。

直交展開の中身を1本の標本で追いかけた図です。①の黒い曲線が展開したい関数 $X(t)$、②はそれを基底ごとの成分 $c_k\varphi_k(t)$ に分けたもので、係数 $c_k=\int X\varphi_k\,dt$ は「$X$ を $\varphi_k$ 方向へ射影した長さ」にあたります。この例では $c_2=+0.289$ が最も大きく、$c_1=+0.022$ は小さい——つまり「どの基底が効くか」は関数ごとに違います。③は成分を1本ずつ足し上げた部分和で、$N$ を増やすほど黒い曲線に近づいていくのが見えます。問題は「この基底 $\{\psi_k\}$ をどう選べば、少ない $N$ で最も速く近づくか」です。

さて、ここでどんな要求を追加すれば基底が絞り込まれるかを考えます。KL展開が課す要求はただ一つ、係数どうしが無相関であることです。

$$ \begin{equation} E[c_j c_k] = \lambda_j \delta_{jk} \end{equation} $$

なぜこれが欲しいのでしょうか。無相関だと、各モードの寄与が「独立した情報の断片」として扱えるからです。全体のエネルギーがモードごとにきれいに足し算で分解され、あるモードを捨てたときの損失が他のモードと干渉しません。逆に係数が相関していると、モード1とモード2が同じ情報を重複して持っていることになり、「$N$ 本で表す」という圧縮の目的から見て無駄が生じます。PCAで主成分どうしが無相関になるように軸を選んだのと、まったく同じ発想です。

「無相関にせよ」という要求がどれだけ強いかを、ブラウン運動の係数 $(c_1,c_2)$ の散布図で見てみます。左はKL基底（固有関数）を使った場合で、点群は座標軸に沿った横長の楕円になり、相関係数は $\rho=-0.006$ とほぼゼロです。右は同じ2次元空間を張る基底を35度だけ回転させた場合で、まったく同じ空間を張っているにもかかわらず点群は斜めに傾き、$\rho=-0.770$ と強い相関が出ています。張る空間が同じでも係数の相関は基底の取り方で変わる——だからこそ「無相関」という要求が、回転の自由度を潰して基底を一意に決められるのです。

では、この要求から何が導かれるか計算しましょう。$c_j$ の定義式を代入して期待値を取ります。

$$ E[c_j c_k] = E\left[\int_0^T X(s)\psi_j(s)\,ds \int_0^T X(t)\psi_k(t)\,dt\right] $$

積分どうしの積を二重積分にまとめ、期待値と積分の順序を交換すると（2次モーメントの有限性と $R$ の連続性からフビニの定理が使えます）、

$$ E[c_j c_k] = \int_0^T\!\!\int_0^T E[X(s)X(t)]\,\psi_j(s)\psi_k(t)\,ds\,dt = \int_0^T\!\!\int_0^T R(s,t)\psi_j(s)\psi_k(t)\,ds\,dt $$

を得ます。ここで内側の $s$ 積分を先に実行した量に名前を付けましょう。

$$ \begin{equation} (\mathcal{R}\psi_j)(t) := \int_0^T R(s,t)\psi_j(s)\,ds \end{equation} $$

この $\mathcal{R}$ は、関数を入れると関数が返ってくる作用素で、共分散作用素と呼ばれます。共分散行列がベクトルに掛かるのとまったく同じ役割です。この記号を使うと、先ほどの式は

$$ E[c_j c_k] = \int_0^T (\mathcal{R}\psi_j)(t)\,\psi_k(t)\,dt = \langle \mathcal{R}\psi_j, \psi_k\rangle $$

と簡潔に書けます。要求 $E[c_jc_k]=\lambda_j\delta_{jk}$ は、したがって

$$ \langle \mathcal{R}\psi_j, \psi_k\rangle = \lambda_j\delta_{jk} = \langle \lambda_j \psi_j, \psi_k \rangle \quad (\forall k) $$

と同値です。移項すると、全ての $k$ について

$$ \langle \mathcal{R}\psi_j – \lambda_j \psi_j,\ \psi_k \rangle = 0 $$

が成り立つ、ということになります。ここで $\{\psi_k\}$ が $L^2([0,T])$ の完全正規直交系であれば、全ての基底ベクトルと直交する関数はゼロ関数しかありません。よって

$$ \begin{equation} \mathcal{R}\psi_j = \lambda_j \psi_j, \qquad \text{すなわち} \qquad \int_0^T R(s,t)\psi_j(s)\,ds = \lambda_j \psi_j(t) \end{equation} $$

が結論されます。これが第二種フレドホルム積分方程式であり、$\psi_j$ は共分散作用素 $\mathcal{R}$ の固有関数、$\lambda_j$ はその固有値です。

大事なのは論理の向きです。私たちは固有関数を仮定したのではありません。「係数が無相関になってほしい」という要求だけを置いたところ、基底は積分方程式の解、すなわち共分散作用素の固有関数でなければならないという結論が強制されたのです。KL展開の基底に選択の余地はありません。

もっとも、この積分方程式が本当に解を持つのか、解が完全系をなすのか、はまだ確かめていません。次の節でその保証を与えます。

共分散作用素とマーサーの定理

積分方程式に解があるかどうかは、作用素 $\mathcal{R}$ がどんな性質を持つかで決まります。$\mathcal{R}$ には3つの良い性質があります。

（1）対称（自己共役）である：$R(s,t)=R(t,s)$ なので、

$$ \langle \mathcal{R}f, g\rangle = \int\!\!\int R(s,t)f(s)g(t)\,ds\,dt = \langle f, \mathcal{R}g\rangle $$

が成り立ちます。二重積分の中で $s$ と $t$ の役割を入れ替えても、核が対称だから値が変わらない、というだけの話です。

（2）非負である：任意の $f\in L^2$ に対して

$$ \langle \mathcal{R}f, f\rangle = \int\!\!\int R(s,t)f(s)f(t)\,ds\,dt = E\left[\left(\int_0^T X(t)f(t)\,dt\right)^2\right] \ge 0 $$

です。$R(s,t)=E[X(s)X(t)]$ を代入して期待値を外に出すと、二重積分が $\left(\int Xf\right)^2$ の期待値に化けます。二乗の期待値なので非負、というのは前節の非負定値性の連続版そのものです。

（3）コンパクト（ヒルベルト・シュミット）である：$R$ が連続で $[0,T]^2$ が有界閉集合なので $\int\!\int R(s,t)^2\,ds\,dt<\infty$ が成り立ち、$\mathcal{R}$ はヒルベルト・シュミット作用素、したがってコンパクトです。

この3つが揃うと、コンパクト自己共役作用素のスペクトル定理が使えます。結論はこうです。$\mathcal{R}$ は可算個の非負固有値

$$ \lambda_1 \ge \lambda_2 \ge \lambda_3 \ge \cdots \ge 0, \qquad \lambda_k \to 0 $$

を持ち、対応する固有関数 $\{\varphi_k\}$ は $L^2([0,T])$ の正規直交系をなします（同じ固有値に属する固有関数はグラム・シュミットで直交化できます）。固有値が非負なのは性質（2）から、$\lambda_k\to0$ はコンパクト性から来ます。

さらに、核 $R$ が連続かつ非負定値であるとき、次のマーサーの定理が成り立ちます。

$$ \begin{equation} R(s,t) = \sum_{k=1}^{\infty} \lambda_k \varphi_k(s)\varphi_k(t) \end{equation} $$

しかもこの級数は $[0,T]^2$ 上で絶対一様収束します。これは対称行列のスペクトル分解 $\bm{\Sigma}=\sum_k \lambda_k \bm{u}_k\bm{u}_k^\top$ の連続版です。行列版で $\bm{u}_k\bm{u}_k^\top$ という「ランク1行列」を足し上げたのと同じように、ここでは $\varphi_k(s)\varphi_k(t)$ という「ランク1の核」を足し上げています。

マーサーの定理を、部分和を実際に足し上げて確かめた図です。$K=1$ 項だけでは $\varphi_1(s)\varphi_1(t)$ という滑らかな山しか作れず、真の核（右下）が持つ対角線 $s=t$ の「折れ目」がまったく再現できていません（最大誤差0.1882）。項数を増やすと折れ目が徐々に立ち上がり、$K=5$ で0.0392、$K=20$ で0.0089、$K=100$ で0.0013と誤差が単調に減っていきます。誤差が「どこか1点だけ大きい」のではなく画面全体で一様に小さくなっていく点が、マーサー級数が絶対一様収束することの見た目上の意味です。

マーサーの定理から、直ちに便利な系が2つ出ます。1つ目は、$s=t$ とおいて

$$ R(t,t) = \sum_{k=1}^\infty \lambda_k \varphi_k(t)^2 $$

これは「時刻 $t$ における分散が、各モードの寄与の和に分解される」ことを意味します。2つ目は、これを $t$ で積分して $\int\varphi_k^2=1$ を使うと

$$ \begin{equation} \int_0^T R(t,t)\,dt = \sum_{k=1}^\infty \lambda_k \end{equation} $$

となることです。左辺は過程の全エネルギー（時間積分した分散、$E\int X^2dt$）ですから、固有値の総和は過程の総分散に等しいという、行列のトレースとまったく同じ関係が成り立ちます。この関係は後で打ち切り誤差を計算するときに主役になります。

「固有値の総和＝総分散」を2つの見方で描いたものです。左は固有値を大きい順に足していった累積和で、$N=1$ で早くも0.405に達し、そのまま総分散 $\int_0^1 t\,dt=0.5$（赤い破線）へ下から漸近していきます。右は同じ等式を積分する前の姿、つまり式 $R(t,t)=\sum_k\lambda_k\varphi_k(t)^2$ で、各時刻の分散が色分けされたモードの積み重ねに分解される様子です。最も濃い第1モードの層が全体の大半を占め、上に乗る第2モード以降の層は薄く、それらを足し上げると赤い破線 $R(t,t)=t$ にぴったり近づいていくことが読み取れます。

道具立てが完成しました。ここからが本題で、いよいよKL展開の定理を述べて証明します。

カルーネン・レーベ展開の定理

改めて定理の形にまとめます。

定理（Karhunen–Loève） $\{X(t)\}_{t\in[0,T]}$ を平均0・連続な共分散関数 $R$ を持つ2次確率過程とし、$(\lambda_k,\varphi_k)$ を $\mathcal{R}$ の固有値・固有関数（$\lambda_1\ge\lambda_2\ge\cdots\ge0$、$\{\varphi_k\}$ は正規直交）とする。このとき $$X(t) = \sum_{k=1}^{\infty}\sqrt{\lambda_k}\,\xi_k\,\varphi_k(t)$$ が $[0,T]$ 上一様に平均二乗収束する。ここで $$\xi_k = \frac{1}{\sqrt{\lambda_k}}\int_0^T X(t)\varphi_k(t)\,dt$$ であり、$E[\xi_k]=0$、$E[\xi_j\xi_k]=\delta_{jk}$ を満たす。

まず係数の統計的性質を確かめます。$c_k=\langle X,\varphi_k\rangle$ とおくと、期待値は

$$ E[c_k] = \int_0^T E[X(t)]\varphi_k(t)\,dt = 0 $$

です（中心化を仮定したので $E[X(t)]=0$）。共分散は前節の計算そのままで、

$$ E[c_jc_k] = \langle \mathcal{R}\varphi_j, \varphi_k\rangle = \langle \lambda_j\varphi_j,\varphi_k\rangle = \lambda_j\delta_{jk} $$

となります。2つ目の等号で固有方程式 $\mathcal{R}\varphi_j=\lambda_j\varphi_j$ を使い、3つ目で正規直交性を使いました。したがって $c_k$ の分散は $\lambda_k$、$j\ne k$ なら無相関。$\xi_k=c_k/\sqrt{\lambda_k}$ と正規化すれば分散1になります。固有値 $\lambda_k$ は、そのモードが持つ分散そのものだとわかります。PCAで固有値が主成分の分散だったのと同じです。

次に収束を示します。$N$ 項で打ち切った近似を

$$ X_N(t) = \sum_{k=1}^{N}\sqrt{\lambda_k}\,\xi_k\varphi_k(t) = \sum_{k=1}^N c_k \varphi_k(t) $$

とおき、各時刻での平均二乗誤差 $E[(X(t)-X_N(t))^2]$ を計算します。二乗を展開すると3つの項が出ます。

$$ E[(X-X_N)^2] = E[X(t)^2] – 2E[X(t)X_N(t)] + E[X_N(t)^2] $$

第1項は定義から $R(t,t)$ です。第2項の交差項を計算するために、$E[X(t)c_k]$ を求めましょう。$c_k$ の定義を代入して期待値と積分を交換すると

$$ E[X(t)c_k] = \int_0^T E[X(t)X(s)]\varphi_k(s)\,ds = \int_0^T R(s,t)\varphi_k(s)\,ds = \lambda_k\varphi_k(t) $$

となります。最後の等号は固有方程式そのものです。したがって

$$ E[X(t)X_N(t)] = \sum_{k=1}^N E[X(t)c_k]\varphi_k(t) = \sum_{k=1}^N \lambda_k\varphi_k(t)^2 $$

です。第3項は $E[c_jc_k]=\lambda_k\delta_{jk}$ を使って

$$ E[X_N(t)^2] = \sum_{j=1}^N\sum_{k=1}^N E[c_jc_k]\varphi_j(t)\varphi_k(t) = \sum_{k=1}^N \lambda_k\varphi_k(t)^2 $$

となり、第2項とぴったり同じ値になります。ここが気持ちのいいところで、交差項と二乗項が同じ大きさなので、まとめると

$$ \begin{equation} E[(X(t)-X_N(t))^2] = R(t,t) – \sum_{k=1}^{N}\lambda_k\varphi_k(t)^2 \end{equation} $$

という簡潔な式が残ります。マーサーの定理より右辺は $\sum_{k>N}\lambda_k\varphi_k(t)^2$ に等しく、しかもマーサー級数は一様収束するので、$N\to\infty$ で $t$ について一様に0へ向かいます。これで定理が示せました。

なお、$X$ がガウス過程の場合は話がもっと良くなります。$c_k=\langle X,\varphi_k\rangle$ はガウス確率変数の（積分という）線形結合なのでガウス分布に従い、ガウスの世界では無相関と独立が一致します。したがって $\xi_k \sim \mathcal{N}(0,1)$ が独立同分布になり、KL展開は「独立な標準正規乱数から確率過程を組み立てるレシピ」として使えるようになります。ブラウン運動のシミュレーションにKL展開が使えるのはこのためです。

さて、式(9)を $t$ で積分すると、区間全体で見た誤差が得られます。次はそれを見ていきましょう。

打ち切り誤差は $\sum_{k>N}\lambda_k$

$N$ 項で打ち切ったときの誤差を、区間全体で測ります。自然な尺度は $L^2$ ノルムの期待値です。

$$ \varepsilon_N^2 := E\left[\int_0^T \big(X(t)-X_N(t)\big)^2 dt\right] = E\big[\|X-X_N\|^2\big] $$

期待値と積分を交換して式(9)を代入すると

$$ \varepsilon_N^2 = \int_0^T\left(R(t,t) – \sum_{k=1}^N \lambda_k\varphi_k(t)^2\right)dt $$

です。第1項は式(8)より $\sum_{k=1}^\infty\lambda_k$、第2項は $\int\varphi_k^2\,dt=1$ より $\sum_{k=1}^N\lambda_k$ になります。したがって

$$ \begin{equation} \varepsilon_N^2 = \sum_{k=1}^{\infty}\lambda_k – \sum_{k=1}^{N}\lambda_k = \sum_{k>N}\lambda_k \end{equation} $$

打ち切り誤差は、捨てた固有値の和にちょうど等しい。 これがKL展開の最も使いやすい性質です。誤差を評価するのに実際の波形を一切見る必要がなく、固有値のリストを見て「上から何個まで足せば全体の99%に届くか」を数えるだけで済みます。PCAの寄与率・累積寄与率とまったく同じ道具立てです。

打ち切り誤差の正体を絵にすると左のようになります。棒グラフが固有値 $\lambda_k$（＝そのモードの分散）で、赤い破線より左の $N=5$ 本を採用すると全分散の95.96%を回収でき、右側に残った灰色の棒の面積の合計 $\sum_{k>5}\lambda_k=0.02020$ がそのまま誤差になります。右は $N$ を横軸にとった誤差の減り方を両対数で見たもので、直線に乗っており、傾きが $-1$、すなわち $\varepsilon_N^2\approx 1/(\pi^2N)$ です。指数関数的に落ちるのではなく $N$ に反比例してしか減らないので、ブラウン運動は「最初の数本は劇的に効くが、そこから先が苦しい」タイプの過程だとわかります。

累積寄与率も同じように定義できます。

$$ \rho_N = \frac{\sum_{k=1}^N \lambda_k}{\sum_{k=1}^\infty \lambda_k} $$

分母は式(8)より $\int_0^T R(t,t)dt$、つまり総分散です。$\rho_N=0.99$ なら「$N$ 項で全分散の99%を説明できた」ことになります。

固有値の減衰の速さが、そのまま圧縮の効きやすさになります。共分散関数が滑らかな過程（たとえばガウスカーネル $R(s,t)=\exp(-(s-t)^2/2\ell^2)$）では固有値が指数的に減るので、数項で十分な精度が出ます。逆に、ブラウン運動のように標本路がギザギザで滑らかさが低い過程では、後で見るように固有値は $\lambda_k\sim 1/(k^2\pi^2)$ 程度の多項式オーダーでしか減らず、同じ精度を出すのにずっと多くの項が要ります。「過程が滑らかなら少数モードで表せる」という直感が、固有値の減衰率として定量化されるわけです。

ここまでは「固有関数を使ったときの誤差」を計算しただけです。まだ「それが最良である」ことは示していません。次節でそれを証明します。

最適性の証明（ラグランジュ乗数法）

いよいよ本丸です。示したいのは次の主張です。

主張（最適性） 任意の正規直交系 $\{\psi_1,\dots,\psi_N\}$ に対し、その張る空間へ射影したときの平均二乗誤差は $$E\Big[\big\|X – \textstyle\sum_{k=1}^N\langle X,\psi_k\rangle\psi_k\big\|^2\Big] \ \ge\ \sum_{k>N}\lambda_k$$ を満たし、等号は $\{\psi_k\}$ が上位 $N$ 個の固有関数と同じ空間を張るときに成り立つ。

まず、目的関数を扱いやすい形に書き直します。$\hat{X}=\sum_{k=1}^N\langle X,\psi_k\rangle\psi_k$ とおくと、これは $X$ の部分空間への直交射影です。射影の性質（ピタゴラスの定理）から

$$ \|X-\hat{X}\|^2 = \|X\|^2 – \|\hat{X}\|^2 = \|X\|^2 – \sum_{k=1}^N \langle X,\psi_k\rangle^2 $$

が成り立ちます。$\|\hat X\|^2=\sum_k\langle X,\psi_k\rangle^2$ となるのは $\{\psi_k\}$ が正規直交だからです。両辺の期待値を取ると、$E\|X\|^2=\int R(t,t)dt=\sum_k\lambda_k$（$N$ に依存しない定数）で、

$$ E\big[\langle X,\psi_k\rangle^2\big] = \langle \mathcal{R}\psi_k,\psi_k\rangle $$

は前に計算した通りです。よって

$$ \begin{equation} J(\psi_1,\dots,\psi_N) := E\|X-\hat X\|^2 = \sum_{k=1}^\infty \lambda_k – \sum_{k=1}^{N}\langle \mathcal{R}\psi_k,\psi_k\rangle \end{equation} $$

となります。第1項は定数なので、誤差 $J$ を最小化することは、$\sum_{k=1}^N\langle\mathcal{R}\psi_k,\psi_k\rangle$ を最大化することと同値です。PCAで「誤差最小化＝分散最大化」だったのと完全に同じ構図が現れました。$\langle\mathcal{R}\psi_k,\psi_k\rangle$ は係数 $\langle X,\psi_k\rangle$ の分散なので、これは文字通り「射影後の分散の総和を最大にせよ」という問題です。

そこで、制約付き最大化問題

$$ \max_{\{\psi_k\}} \ \sum_{k=1}^N \langle \mathcal{R}\psi_k,\psi_k\rangle \quad \text{subject to} \quad \langle\psi_j,\psi_k\rangle = \delta_{jk} $$

をラグランジュ乗数法で解きます。制約は $N(N+1)/2$ 本（対称なので）あるので、乗数も対称行列 $\mu_{jk}=\mu_{kj}$ の形で用意します。ラグランジアンは

$$ L = \sum_{k=1}^N \langle \mathcal{R}\psi_k,\psi_k\rangle – \sum_{j=1}^N\sum_{k=1}^N \mu_{jk}\big(\langle\psi_j,\psi_k\rangle – \delta_{jk}\big) $$

です。これを関数 $\psi_l$ について変分（汎関数微分）します。$\psi_l\to\psi_l+\epsilon\eta$ と摂動させて $\epsilon$ の1次の項を拾いましょう。第1項からは、$\mathcal{R}$ が自己共役であることを使って

$$ \langle \mathcal{R}(\psi_l+\epsilon\eta),\psi_l+\epsilon\eta\rangle = \langle\mathcal{R}\psi_l,\psi_l\rangle + 2\epsilon\langle \mathcal{R}\psi_l,\eta\rangle + O(\epsilon^2) $$

が出ます（交差項が2つとも $\langle\mathcal{R}\psi_l,\eta\rangle$ にまとまるのは自己共役性のおかげです）。第2項からは、$l$ を含む制約を集めて同様に展開すると $-2\epsilon\sum_j \mu_{jl}\langle\psi_j,\eta\rangle$ が出ます（$\mu$ の対称性を使いました）。合わせて、任意の摂動 $\eta$ について1次の項が消える条件は

$$ \left\langle \mathcal{R}\psi_l – \sum_{j=1}^N \mu_{jl}\psi_j,\ \eta \right\rangle = 0 \quad (\forall \eta) $$

すなわち

$$ \begin{equation} \mathcal{R}\psi_l = \sum_{j=1}^N \mu_{jl}\,\psi_j \end{equation} $$

です。この式が意味するのは、「$\mathcal{R}$ で写しても $\{\psi_1,\dots,\psi_N\}$ の張る空間 $V$ から出ない」、つまり $V$ が $\mathcal{R}$ の不変部分空間であるということです。

ここで一手間かけます。$\mu=(\mu_{jl})$ は実対称行列なので、直交行列 $\bm{Q}$ で対角化できます：$\bm{Q}^\top\mu\bm{Q}=\mathrm{diag}(\lambda’_1,\dots,\lambda’_N)$。新しい基底 $\tilde\psi_k=\sum_j Q_{jk}\psi_j$ を取ると、これは正規直交性を保ったまま（直交変換だから）同じ空間 $V$ を張り、しかも式(12)は

$$ \mathcal{R}\tilde\psi_k = \lambda’_k \tilde\psi_k $$

という対角形になります。つまり停留点では、基底を空間 $V$ の中で回転させるだけで、必ず $\mathcal{R}$ の固有関数に取り直せるのです。目的関数 $\sum_k\langle\mathcal{R}\psi_k,\psi_k\rangle=\mathrm{tr}(\mu)$ は回転で不変なので、この取り直しで値は変わりません。

こうして停留点では $\{\tilde\psi_k\}$ が固有関数であり、そのとき目的関数の値は

$$ \sum_{k=1}^N \langle \mathcal{R}\tilde\psi_k,\tilde\psi_k\rangle = \sum_{k=1}^N \lambda’_k $$

つまり「選んだ $N$ 個の固有値の和」になります。これを最大にするには、当然大きい方から $N$ 個を選べばよい。すなわち $\lambda_1,\dots,\lambda_N$ を選び、$V$ は上位 $N$ 個の固有関数が張る空間になります。このとき式(11)より誤差は

$$ J = \sum_{k=1}^\infty\lambda_k – \sum_{k=1}^N\lambda_k = \sum_{k>N}\lambda_k $$

となり、主張が示されました。

この証明の読みどころを整理しておきます。ラグランジュ乗数法が教えてくれたのは「固有関数そのもの」ではなく、「最適な部分空間は $\mathcal{R}$ の不変部分空間である」という構造でした。個々の $\psi_k$ が固有関数である必要はなく（回転の自由度がある）、張る空間が正しければ誤差は同じです。実際、上位 $N$ 個の固有関数を任意に直交回転させても、誤差 $\sum_{k>N}\lambda_k$ は変わりません。「係数を無相関にしたい」という追加要求を課すと、その回転の自由度が潰れて固有関数そのものに固定される、という関係になっています。

ラグランジュ乗数法は停留点を見つける手法なので、厳密には「それが最大値である」ことの保証が要ります。次節でその穴を埋める、乗数法によらない直接評価を与えます。

最適性のもう一つの証明（不等式による直接評価）

停留点の議論を経由せず、いきなり不等式で押さえる方法を紹介します。行列に対するKy Fanの定理の連続版です。

任意の正規直交系 $\{\psi_1,\dots,\psi_N\}$ を取り、それぞれを固有関数系 $\{\varphi_j\}$ で展開します（$\{\varphi_j\}$ は完全系とします）。

$$ \psi_k = \sum_{j=1}^{\infty} a_{jk}\,\varphi_j, \qquad a_{jk} = \langle \psi_k, \varphi_j\rangle $$

$\psi_k$ が単位ノルムなので、パーセバルの等式から $\sum_j a_{jk}^2 = 1$ です。また

$$ \langle \mathcal{R}\psi_k,\psi_k\rangle = \Big\langle \sum_j a_{jk}\lambda_j\varphi_j,\ \sum_i a_{ik}\varphi_i \Big\rangle = \sum_{j=1}^\infty \lambda_j a_{jk}^2 $$

となります（$\mathcal{R}\varphi_j=\lambda_j\varphi_j$ を代入し、正規直交性で二重和が対角成分だけ残ります）。したがって最大化したい量は

$$ S := \sum_{k=1}^N \langle\mathcal{R}\psi_k,\psi_k\rangle = \sum_{j=1}^{\infty}\lambda_j\, d_j, \qquad d_j := \sum_{k=1}^{N} a_{jk}^2 $$

と書けます。この重み $d_j$ には2つの制約があります。

制約A：$\sum_j d_j = N$。 各 $k$ について $\sum_j a_{jk}^2=1$ なので、$k$ について足すと $N$ になります。

制約B：$0 \le d_j \le 1$。 $d_j=\sum_{k=1}^N \langle\varphi_j,\psi_k\rangle^2$ は、単位ベクトル $\varphi_j$ を部分空間 $V=\mathrm{span}\{\psi_k\}$ に射影した長さの二乗です。射影は長さを伸ばさないので $d_j\le\|\varphi_j\|^2=1$。これはベッセルの不等式そのものです。

これで問題は、「$0\le d_j\le1$、$\sum_j d_j=N$ という条件のもとで $\sum_j \lambda_j d_j$ を最大化せよ」という単純な線形計画に化けました。$\lambda_1\ge\lambda_2\ge\cdots$ と並んでいるので、限られた重みの総量 $N$ は大きい $\lambda_j$ に優先的に、上限1まで割り当てるのが最善です。

この線形計画がどれほど単純かを、$N=3$ の場合の3通りの割り当てで比べたのが上の図です。灰色の棒が固有値 $\lambda_j$、色付きの棒が実際に回収できる分 $\lambda_j d_j$ を表します。左の「上位3個に $d=1$ を集中させる」割り当てが $S=0.46653$ で最大、中央のように6個へ薄くばらまくと $S=0.28657$ まで落ち、右のように大きい2個を外すと $S=0.02949$ と壊滅的です。重みの総量が $N$ に固定されている以上、大きい $\lambda_j$ から順に上限まで詰めるしかない——これが最適性の本質で、式で書けば

$$ S \le \lambda_1 + \lambda_2 + \cdots + \lambda_N $$

であり、等号は $d_1=\cdots=d_N=1$、$d_j=0\ (j>N)$ のとき、すなわち $V$ が上位 $N$ 個の固有関数の張る空間と一致するときに限ります。式(11)に戻せば

$$ E\|X-\hat X\|^2 = \sum_j \lambda_j – S \ \ge\ \sum_{k>N}\lambda_k $$

が得られ、最適性が不等式として直接示されました。ラグランジュ乗数法の議論と違い、こちらは停留点かどうかを問わず全ての正規直交系を一網打尽にしています。2つの証明を並べて見ると、乗数法は「なぜ固有関数という形になるのか」を、直接評価は「なぜそれが最良なのか」を、それぞれ説明していることがわかります。

理論はここまでです。次は、この連続の話が離散のPCAとどう繋がるかを見ておきましょう。

PCAとの対応 — 離散化するとただの固有値分解

実務では、確率過程を連続関数として扱うことはまずありません。$[0,T]$ を $n$ 分割して $t_i=(i-\tfrac12)\Delta t$（$\Delta t=T/n$）でサンプリングし、ベクトル $\bm{X}=(X(t_1),\dots,X(t_n))^\top$ として扱います。このとき積分方程式はどうなるでしょうか。

積分を中点則で近似すると

$$ \int_0^T R(s,t_i)\varphi(s)\,ds \approx \sum_{j=1}^n R(t_j,t_i)\varphi(t_j)\,\Delta t = \lambda\varphi(t_i) $$

です。ここで共分散行列 $\bm{R}$ を $R_{ij}=R(t_i,t_j)$ で定義し、ベクトル $\bm{\varphi}=(\varphi(t_1),\dots,\varphi(t_n))^\top$ を置くと、これは

$$ \begin{equation} (\Delta t)\,\bm{R}\,\bm{\varphi} = \lambda\,\bm{\varphi} \end{equation} $$

というただの行列固有値問題になります。共分散行列に $\Delta t$ を掛けた行列を固有分解するだけです。実務でKL展開を「やる」というのは、ほとんどの場合この行列固有値問題を解くことを意味します。

有限次元と無限次元を並べると、対応関係が一目でわかります。左は $\mathbb{R}^2$ の点群に対するPCAで、共分散「行列」の固有ベクトル $u_1,u_2$ が点群の広がりの主軸として引かれ、矢印の長さは $\sqrt{\lambda_k}$ に比例しています。右は同じことを曲線の束に対して行ったもので、灰色の細線1本1本がデータ曲線 $X^{(m)}(t)$、太い曲線が共分散「関数」の固有関数 $\varphi_1,\varphi_2$（振幅は $\sqrt{\lambda_k}$ でスケール）です。第1固有関数が曲線束の最も広がっている方向、つまり「終端に向かって単調に伸びる成分」を捉えていることが見て取れます。矢印が曲線に変わっただけで、やっていることは完全に同じです。

対応関係を表で整理しておきましょう。

正規化の対応には注意が必要です。連続版の $\int \varphi_k^2 dt = 1$ を離散化すると $\sum_i \varphi(t_i)^2\Delta t = 1$、つまり $\|\bm{u}_k\|=1$ となる固有ベクトルに対して $\varphi(t_i) = u_{k,i}/\sqrt{\Delta t}$ という関係になります。この $1/\sqrt{\Delta t}$ のスケーリングを忘れると、固有関数の振幅が合わずに首をかしげることになります。あとで実装するときに実際に確認します。

もう一つ実務上の注意です。ここでは $R$ が理論的にわかっている前提で書きましたが、現実には $M$ 本の観測サンプルから標本共分散行列 $\hat{\bm{R}} = \frac{1}{M}\sum_{m}\bm{x}^{(m)}\bm{x}^{(m)\top}$ を作って固有分解します。これはまさにPCAです。つまり「経験的KL展開＝PCA」であり、気象学のEOF解析も、関数データ解析のfPCAも、名前が違うだけで中身は同じ手続きです。

対応関係が見えたところで、具体例に進みましょう。KL展開の解析解が手で求まる、最も有名な例を扱います。

具体例：ブラウン運動のKL展開

ブラウン運動（ウィーナー過程）$\{W(t)\}_{t\in[0,1]}$ を考えます。定義から $E[W(t)]=0$ で、共分散関数は

$$ R(s,t) = E[W(s)W(t)] = \min(s,t) $$

でした。$\min$ という素朴な関数なのに、固有関数がきれいな三角関数になり、しかも手計算で完全に解けます。KL展開の教科書的な例として必ず登場するのはこのためです。

解くべきは積分方程式

$$ \begin{equation} \int_0^1 \min(s,t)\,\varphi(s)\,ds = \lambda\,\varphi(t) \end{equation} $$

です。左辺の $\min$ が扱いにくいので、$s$ の積分区間を $t$ で切って場合分けします。$st$ なら $\min(s,t)=t$ なので、

$$ \begin{equation} \int_0^t s\,\varphi(s)\,ds + t\int_t^1 \varphi(s)\,ds = \lambda\varphi(t) \end{equation} $$

と書けます。これで $\min$ が消え、普通の積分の和になりました。

左が、いま行った場合分けの図解です。$t=0.6$ を固定して核 $\min(s,t)$ を $s$ の関数として描くと、$st$ では高さ $t$ の定数（橙の領域）という折れ線になります。この折れ点で切ることで、積分が $\int_0^t s\varphi(s)ds$ と $t\int_t^1\varphi(s)ds$ の2項に分かれるわけです。右は、これから導く2つの境界条件が固有関数に課す形を先取りしたもので、どのモードも $t=0$ で必ず0を通り（赤い丸）、$t=1$ では接線が水平（点線）になっています。この $\varphi'(1)=0$ が $\cos\omega=0$ を強制し、固有値が $\omega_k=(k-\tfrac12)\pi$ という半端な位置に飛び飛びで並ぶ理由になります。

両辺を $t$ で微分します。 ライプニッツの積分則を使います。第1項は上端が $t$ なので $t\varphi(t)$。第2項は積の微分で、$t$ の微分から $\int_t^1\varphi(s)ds$ が、積分の下端が $t$ であることから $t\cdot(-\varphi(t))$ が出ます。まとめると

$$ \lambda\varphi'(t) = t\varphi(t) + \int_t^1\varphi(s)\,ds – t\varphi(t) = \int_t^1 \varphi(s)\,ds $$

となり、$t\varphi(t)$ が打ち消し合ってきれいな形が残ります。

$$ \begin{equation} \lambda\varphi'(t) = \int_t^1\varphi(s)\,ds \end{equation} $$

もう一度 $t$ で微分します。 右辺は下端が $t$ の積分なので、微分すると $-\varphi(t)$ です。

$$ \begin{equation} \lambda\varphi”(t) = -\varphi(t), \qquad \text{すなわち}\qquad \varphi”(t) + \frac{1}{\lambda}\varphi(t) = 0 \end{equation} $$

積分方程式が、単振動の常微分方程式に化けました。 これがブラウン運動のKL展開が解ける理由です。$\min(s,t)$ という核は、実は微分作用素のグリーン関数になっているのです。

境界条件を拾いましょう。境界条件は元の積分方程式(15)から読み取ります。

$t=0$ を代入すると、左辺の第1項は $\int_0^0=0$、第2項は $t=0$ が掛かるので0。よって

$$ \lambda\varphi(0) = 0 \ \Longrightarrow\ \varphi(0) = 0 \qquad (\lambda\neq0 \text{ のとき}) $$

$t=1$ を式(16)に代入すると、右辺は $\int_1^1=0$ なので

$$ \lambda\varphi'(1) = 0 \ \Longrightarrow\ \varphi'(1) = 0 $$

（$\lambda=0$ の場合は式(15)から $\int_0^t s\varphi ds + t\int_t^1\varphi ds=0$ が全ての $t$ で必要となり、2回微分して $\varphi\equiv0$。固有関数にならないので除外します。また $\lambda<0$ だと $\varphi''=|\lambda|^{-1}\varphi$ の指数解となり、$\varphi(0)=0,\varphi'(1)=0$ を同時に満たすのは $\varphi\equiv0$ のみ。共分散作用素が非負であることとも整合します。）

$\lambda>0$ として $\omega=1/\sqrt{\lambda}$ とおくと、式(17)の一般解は

$$ \varphi(t) = A\sin(\omega t) + B\cos(\omega t) $$

です。$\varphi(0)=0$ から $B=0$。残るのは $\varphi(t)=A\sin(\omega t)$ です。もう一つの条件 $\varphi'(1)=A\omega\cos(\omega)=0$ から、$A\neq0,\omega\neq0$ なので

$$ \cos\omega = 0 \ \Longrightarrow\ \omega = \left(k-\tfrac{1}{2}\right)\pi, \quad k=1,2,3,\dots $$

が固有値条件になります。$\omega=1/\sqrt{\lambda}$ に戻すと、固有値は

$$ \begin{equation} \lambda_k = \frac{1}{\omega_k^2} = \frac{1}{\left(k-\frac{1}{2}\right)^2\pi^2}, \qquad k=1,2,3,\dots \end{equation} $$

最後に振幅 $A$ を正規化条件から決めます。

$$ \int_0^1 A^2\sin^2(\omega_k t)\,dt = A^2\int_0^1 \frac{1-\cos(2\omega_k t)}{2}dt = \frac{A^2}{2} – \frac{A^2\sin(2\omega_k)}{4\omega_k} $$

ここで $\omega_k=(k-\tfrac12)\pi$ なので $2\omega_k=(2k-1)\pi$ となり $\sin(2\omega_k)=0$。よって $A^2/2=1$、$A=\sqrt2$ が得られます。固有関数は

$$ \begin{equation} \varphi_k(t) = \sqrt{2}\,\sin\!\left(\left(k-\tfrac{1}{2}\right)\pi t\right) \end{equation} $$

まとめると、ブラウン運動のKL展開は

$$ \begin{equation} W(t) = \sum_{k=1}^{\infty} \frac{\sqrt{2}}{\left(k-\frac{1}{2}\right)\pi}\,\xi_k \sin\!\left(\left(k-\tfrac{1}{2}\right)\pi t\right), \qquad \xi_k \overset{\text{iid}}{\sim}\mathcal{N}(0,1) \end{equation} $$

となります（ブラウン運動はガウス過程なので、$\xi_k$ は無相関にとどまらず独立です）。

検算をしておきましょう。 式(8)より固有値の総和は総分散に一致するはずです。左辺は

$$ \int_0^1 R(t,t)\,dt = \int_0^1 t\,dt = \frac{1}{2} $$

右辺は

$$ \sum_{k=1}^\infty \frac{1}{(k-\frac12)^2\pi^2} = \frac{1}{\pi^2}\sum_{k=1}^\infty\frac{1}{(k-\frac12)^2} = \frac{1}{\pi^2}\cdot\frac{\pi^2}{2} = \frac{1}{2} $$

ぴったり一致します（奇数の逆二乗和 $\sum_{m\ \text{odd}}1/m^2=\pi^2/8$ を4倍しました）。理論が正しく機能している確かな証拠です。

固有関数の形にも意味があります。 $\varphi_1(t)=\sqrt2\sin(\pi t/2)$ は $t=0$ で0から始まり $t=1$ で最大となる、単調増加のなだらかな曲線です。ブラウン運動が「原点から出発してだんだん広がっていく」という性質を、これ1本でかなり捉えています。実際 $\lambda_1/\sum\lambda_k = 0.405285/0.5 = 81.06\%$ で、たった1モードでブラウン運動の全分散の8割を説明できます。境界条件 $\varphi(0)=0$ は「$W(0)=0$ で不確かさゼロ」、$\varphi'(1)=0$ は「終端では制約がなく、そこで各モードが極値をとる」という物理的な意味に対応しています。

固有値の減衰は $\lambda_k \approx 1/(k^2\pi^2)$ で、たった2次のオーダーです。累積寄与率を並べると 81.06%、90.06%、93.31%、94.96%、95.96% と、最初こそ勢いよく伸びますが、そこからの伸びが鈍い。95%は $N=5$ で越えるのに、99%に届くには $N=21$、99.9%となると $N=203$ も要ります。ブラウン運動の標本路が至る所微分不可能なほどギザギザであることが、この遅い減衰に表れています。

具体例2：ブラウン橋のKL展開

もう一つ、ほとんど同じ計算で解ける例を挙げます。ブラウン橋 $B(t)=W(t)-tW(1)$（両端を0に固定したブラウン運動）の共分散関数は

$$ R(s,t) = \min(s,t) – st $$

です。積分方程式に代入して $\min$ を分解すると

$$ \int_0^t s\varphi(s)ds + t\int_t^1\varphi(s)ds – t\int_0^1 s\varphi(s)ds = \lambda\varphi(t) $$

となります。増えたのは最後の項ですが、これは $t$ の1次式なので2回微分すると消えます。したがって微分方程式は先ほどとまったく同じ $\lambda\varphi”=-\varphi$ です。

変わるのは境界条件です。$t=0$ を代入すると全項が0なので $\varphi(0)=0$。$t=1$ を代入すると、第2項は $\int_1^1=0$、第1項と第3項が $\int_0^1 s\varphi ds$ で打ち消し合って0。よって $\varphi(1)=0$ です。両端がディリクレ条件になるので

$$ \varphi_k(t) = \sqrt2\sin(k\pi t), \qquad \lambda_k = \frac{1}{k^2\pi^2}, \qquad k=1,2,\dots $$

が得られます。これは通常のフーリエ正弦級数そのものです。検算すると $\int_0^1(t-t^2)dt = 1/2-1/3=1/6$ で、$\sum_k 1/(k\pi)^2 = (1/\pi^2)(\pi^2/6)=1/6$。やはり一致します。

2つの過程を3つの側面から並べた図です。左の標本路では、青のブラウン運動が終端で自由に散らばるのに対し、赤のブラウン橋は $t=1$ に向かって一点に絞られていきます。中央の固有関数を見ると、実線（ブラウン運動）は $t=1$ で山や谷の頂点に達して傾きが0になり、破線（ブラウン橋）は $t=1$ でちょうど0に戻る——これが $\varphi'(1)=0$ と $\varphi(1)=0$ という境界条件の違いそのものです。右の固有値は、ブラウン橋のほうが同じ $k$ で常に小さく、15項までの和も0.493と0.160（極限は $1/2$ と $1/6$）と3倍の差がありますが、どちらも $k^{-2}$ の多項式減衰であることは共通しています。

ブラウン運動とブラウン橋を比べると、KL展開の本質が見えます。 微分方程式は同じで、違いは境界条件だけ。「終端を固定するかしないか」という確率過程の性質の差が、そのまま固有関数の境界条件の差として現れています。共分散関数が過程の性質を丸ごと持っており、固有関数はそれを忠実に読み取っているのです。

ここまで手で解いてきました。次は、この解析解が本当に正しいかをPythonで確かめます。

Pythonでの実装

やることは4つです。①解析解と離散固有分解を突き合わせる、②固有関数と固有値スペクトルを可視化する、③実際のブラウン運動サンプルを再構成して誤差が $\sum_{k>N}\lambda_k$ になるか測る、④KL基底が他の基底より本当に優れているかを比べる。

まずは共通の設定と、①の照合からです。

import numpy as np
import matplotlib
import matplotlib.pyplot as plt

# 日本語フォント設定
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

# 区間 [0,1] を n 分割（中点則）
n = 400
dt = 1.0 / n
t_grid = (np.arange(1, n + 1) - 0.5) / n

# ブラウン運動の共分散行列 R_ij = min(t_i, t_j)
R = np.minimum.outer(t_grid, t_grid)

# 離散化した固有値問題 (Δt) R u = λ u を解く
w, V = np.linalg.eigh(R * dt)
order = np.argsort(w)[::-1]          # 大きい順に並べ替え
lam_num = w[order]
U = V[:, order]

# 解析解 λ_k = 1/((k-1/2)π)^2
K = 8
lam_ana = np.array([1.0 / ((k - 0.5) * np.pi) ** 2 for k in range(1, K + 1)])

print("k :   数値解      解析解     相対誤差")
for k in range(K):
    rel = abs(lam_num[k] - lam_ana[k]) / lam_ana[k]
    print(f"{k+1} : {lam_num[k]:.6f}  {lam_ana[k]:.6f}  {rel:.2e}")
print(f"\n固有値の総和（数値）: {lam_num.sum():.6f}  / 理論値: 0.5")

出力は次のようになります。

k :   数値解      解析解     相対誤差
1 : 0.405285  0.405285  1.29e-06
2 : 0.045032  0.045032  1.16e-05
3 : 0.016212  0.016211  3.21e-05
4 : 0.008272  0.008271  6.30e-05
5 : 0.005004  0.005004  1.04e-04
6 : 0.003350  0.003349  1.56e-04
7 : 0.002399  0.002398  2.17e-04
8 : 0.001802  0.001801  2.89e-04

固有値の総和（数値）: 0.500000  / 理論値: 0.5

$400\times400$ の行列を固有分解しただけで、手計算した $1/((k-\frac12)\pi)^2$ が小数点以下6桁まで再現されました。相対誤差が高次のモードで大きくなるのは、振動が速い固有関数ほど有限グリッドでの離散化誤差を受けやすいからです。固有値の総和も0.500000となり、式(8)の「固有値の総和＝総分散 $\int_0^1 t\,dt = 1/2$」が数値的にも確認できます。

次に固有関数そのものを比べます。離散版の固有ベクトルは $\|\bm{u}_k\|=1$ に正規化されているので、連続版の $\int\varphi^2dt=1$ に合わせるには $1/\sqrt{\Delta t}$ を掛ける必要があります。前節で触れたスケーリングです。

# 固有ベクトル → 固有関数へのスケール変換（∫φ²dt = 1 に合わせる）
phi_num = U / np.sqrt(dt)

fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))

ax = axes[0]
for k in range(1, 5):
    v = phi_num[:, k - 1]
    ana = np.sqrt(2) * np.sin((k - 0.5) * np.pi * t_grid)
    if np.dot(v, ana) < 0:            # 固有ベクトルの符号は不定なので揃える
        v = -v
    ax.plot(t_grid, ana, lw=3, alpha=0.35, label=f"解析解 k={k}")
    ax.plot(t_grid, v, lw=1.2, ls="--", color="k")
ax.set_xlabel("時刻 t"); ax.set_ylabel("固有関数 $\\varphi_k(t)$")
ax.set_title("ブラウン運動の固有関数（太線=解析解, 破線=数値解）")
ax.legend(fontsize=9); ax.grid(alpha=0.3)

ax = axes[1]
ks = np.arange(1, 21)
ax.semilogy(ks, lam_num[:20], "o", label="数値解")
ax.semilogy(ks, 1.0 / ((ks - 0.5) * np.pi) ** 2, "-", label="解析解 $1/((k-1/2)\\pi)^2$")
ax.set_xlabel("モード番号 k"); ax.set_ylabel("固有値 $\\lambda_k$（対数軸）")
ax.set_title("固有値スペクトルの減衰")
ax.legend(); ax.grid(alpha=0.3, which="both")

plt.tight_layout(); plt.show()

左のグラフでは、数値固有ベクトル（破線）が解析解（太い薄線）の上に完全に重なります。実際、両者の最大絶対誤差は $10^{-14}$ 程度で、機械精度の範囲です。$k$ が増えるごとに振動が1回ずつ増え、どのモードも $t=0$ で0（$\varphi(0)=0$）、$t=1$ で傾きが0（$\varphi'(1)=0$）になっている点にも注目してください。手計算で課した境界条件が、行列固有分解でも自動的に満たされています。

右のグラフは固有値の減衰です。対数軸で見ると直線的に落ちておらず、緩やかに曲がっています。これは指数減衰ではなく $\lambda_k\propto k^{-2}$ の多項式減衰であることを示しています。ブラウン運動の標本路が滑らかでないことが、この遅い減衰として現れているわけです。

続いて、累積寄与率を見てみましょう。何モードでどれだけ説明できるかの指標です。

# 解析解の固有値を 200 項ぶん用意し、累積寄与率を計算（分母は総分散 = 0.5）
lam_full = np.array([1.0 / ((k - 0.5) * np.pi) ** 2 for k in range(1, 201)])
cum = np.cumsum(lam_full) / 0.5

plt.figure(figsize=(7, 4.5))
plt.plot(np.arange(1, 51), cum[:50] * 100, "o-", ms=4)
plt.axhline(95, color="r", ls="--", lw=1, label="95%")
plt.axhline(99, color="g", ls="--", lw=1, label="99%")
plt.xlabel("採用モード数 N"); plt.ylabel("累積寄与率 [%]")
plt.title("ブラウン運動：累積寄与率（分散の何%を説明できるか）")
plt.legend(); plt.grid(alpha=0.3); plt.tight_layout(); plt.show()

for N in [1, 2, 3, 5, 10]:
    print(f"N={N:2d}: 累積寄与率 {cum[N-1]*100:.2f}%  打ち切り誤差 {0.5-cum[N-1]*0.5:.5f}")

出力は次の通りです。

N= 1: 累積寄与率 81.06%  打ち切り誤差 0.09472
N= 2: 累積寄与率 90.06%  打ち切り誤差 0.04968
N= 3: 累積寄与率 93.31%  打ち切り誤差 0.03347
N= 5: 累積寄与率 95.96%  打ち切り誤差 0.02020
N=10: 累積寄与率 97.98%  打ち切り誤差 0.01012

曲線の形が「最初だけ急で、あとは寝る」という $k^{-2}$ 減衰の特徴をよく表しています。$N=1$ でいきなり81.06%まで跳ね上がり、$N=5$ で95%線（赤）を越えますが、そこから99%線（緑）に届くまでには $N=21$ まで待たされます。グラフを右端まで伸ばしても $N=50$ で99.59%にしかならず、最後の0.4%を削るコストが非常に高いことがわかります。

第1モードだけで81%、2モードで90%を説明できるのは驚くほど効率的です。一方で、そこから先の伸びは鈍く、10モード使っても98%止まり。「最初の数モードは劇的に効くが、残りを削り切るのは苦しい」 という、$k^{-2}$ 減衰の典型的な振る舞いです。ここで打ち切り誤差の値（$N=10$ で 0.01012）を覚えておいてください。後でモンテカルロ実験の測定値と突き合わせます。

では実際にブラウン運動のサンプルパスを生成し、KL展開で再構成してみます。

rng = np.random.default_rng(0)

n2 = 1024
dt2 = 1.0 / n2
tt = np.arange(1, n2 + 1) / n2

# ブラウン運動のサンプルパスを生成（増分の累積和）
M = 20000
dW = rng.normal(0.0, np.sqrt(dt2), size=(M, n2))
W = np.cumsum(dW, axis=1)

def kl_basis(K):
    """ブラウン運動のKL基底 φ_k(t) = √2 sin((k-1/2)πt) を K 本返す"""
    return np.array([np.sqrt(2) * np.sin((k - 0.5) * np.pi * tt) for k in range(1, K + 1)])

# 1本のパスを N 項で再構成して重ね描き
plt.figure(figsize=(9, 5))
plt.plot(tt, W[0], color="k", lw=1.6, label="真のパス")
for K in [1, 2, 5, 20]:
    B = kl_basis(K)
    coef = W[0] @ B.T * dt2          # ξ_k√λ_k = <W, φ_k>
    plt.plot(tt, coef @ B, lw=1.4, label=f"KL再構成 N={K}")
plt.xlabel("時刻 t"); plt.ylabel("W(t)")
plt.title("ブラウン運動のKL展開による再構成（項数を増やす）")
plt.legend(); plt.grid(alpha=0.3); plt.tight_layout(); plt.show()

$N=1$ の再構成は、真のパスの「全体としての行き先」だけをなぞる、なめらかな1本の弧です。$N=2$ で大きなうねりが1つ加わり、$N=5$ でだいたいの形が追えるようになり、$N=20$ ではかなり細かい起伏まで再現されます。それでも真のパスのギザギザは完全には再現されません。ブラウン運動は至る所微分不可能なので、有限個の滑らかな正弦波では原理的に追いつけないのです。ここに固有値の遅い減衰の意味が視覚的に現れています。

次に、打ち切り誤差の理論値 $\sum_{k>N}\lambda_k$ が本当に当たるかをモンテカルロで検証します。

def mse_of_basis(B, W, dt):
    """基底 B（行が基底関数）へ射影したときの平均二乗誤差 E||W - Ŵ||² を測る"""
    G = B @ B.T * dt                      # グラム行列（直交なら単位行列）
    coef = W @ B.T * dt
    rec = np.linalg.solve(G, coef.T).T @ B
    return np.mean(np.sum((W - rec) ** 2, axis=1) * dt)

print(" N   モンテカルロ   理論値 Σ_{k>N}λ_k")
for K in [1, 2, 4, 8, 16, 32]:
    mc = mse_of_basis(kl_basis(K), W, dt2)
    th = 0.5 - sum(1.0 / ((k - 0.5) * np.pi) ** 2 for k in range(1, K + 1))
    print(f"{K:2d}   {mc:.5f}       {th:.5f}")

結果は次の通りです。

 N   モンテカルロ   理論値 Σ_{k>N}λ_k
 1   0.09417       0.09472
 2   0.04955       0.04968
 4   0.02521       0.02520
 8   0.01267       0.01265
16   0.00635       0.00633
32   0.00317       0.00317

左の両対数プロットでは、モンテカルロ実測（丸）が理論曲線（実線）の上に完全に乗っており、目視ではずれが判別できません。右はその相対差を取り出した棒グラフで、$N=1$ の0.58%が最大、$N=4$ では0.04%まで小さくなり、どの $N$ でも0.6%未満に収まっています。

理論式 $\varepsilon_N^2=\sum_{k>N}\lambda_k$ が小数点以下4桁レベルで的中しています。 わずかなずれは有限サンプル（$M=20000$）と有限グリッド（$n=1024$）による誤差で、どの $N$ でも相対差は0.6%未満に収まっています。注目したいのは、誤差が $N$ を倍にするたびにほぼ半分になっている点です。これは $\sum_{k>N}\lambda_k \approx 1/(\pi^2 N) \approx 0.1013/N$ という漸近形と一致します（$N=32$ なら $0.1013/32=0.00317$）。「精度を1桁上げるには項数を10倍にせよ」という、圧縮としては厳しい条件になっているわけです。

最後に、最適性の主張を実験で確かめます。KL基底と、KL基底ではない2つの正規直交系（フーリエ基底と区分定数基底）を、同じ項数で比べます。

def fourier_basis(K):
    """周期1の三角関数系 {1, √2cos(2πkt), √2sin(2πkt), ...} を K 本"""
    B = [np.ones_like(tt)]
    k = 1
    while len(B) < K:
        B.append(np.sqrt(2) * np.cos(2 * np.pi * k * tt))
        if len(B) < K:
            B.append(np.sqrt(2) * np.sin(2 * np.pi * k * tt))
        k += 1
    return np.array(B[:K])

def block_basis(K):
    """区間を K 等分した区分定数の正規直交基底（Haar型）"""
    B = np.zeros((K, n2))
    edges = np.linspace(0, n2, K + 1).astype(int)
    for i in range(K):
        B[i, edges[i]:edges[i+1]] = 1.0 / np.sqrt((edges[i+1] - edges[i]) * dt2)
    return B

print(" N      KL基底   フーリエ   区分定数")
res = {}
for K in [1, 2, 4, 8, 16, 32]:
    row = [mse_of_basis(f(K), W, dt2) for f in (kl_basis, fourier_basis, block_basis)]
    res[K] = row
    print(f"{K:2d}   {row[0]:.5f}   {row[1]:.5f}   {row[2]:.5f}")

Ks = list(res.keys())
plt.figure(figsize=(7.5, 5))
for i, name in enumerate(["KL基底（最適）", "フーリエ基底", "区分定数基底"]):
    plt.loglog(Ks, [res[K][i] for K in Ks], "o-", label=name)
plt.xlabel("採用する基底の本数 N"); plt.ylabel("平均二乗誤差 $E\\|X-\\hat{X}\\|^2$")
plt.title("基底の選び方による再構成誤差の比較（ブラウン運動）")
plt.legend(); plt.grid(alpha=0.3, which="both"); plt.tight_layout(); plt.show()

出力は次のようになります。

 N      KL基底   フーリエ   区分定数
 1   0.09417   0.16700   0.16700
 2   0.04955   0.14180   0.08300
 4   0.02521   0.05899   0.04150
 8   0.01267   0.02726   0.02080
16   0.00635   0.01314   0.01042
32   0.00317   0.00646   0.00520

左の両対数プロットで、赤のKL基底の線が常に最も下にあります。3本の傾きはどれもほぼ $-1$ で、$N^{-1}$ という減衰の速さ自体は共通なのに、切片（同じ $N$ での誤差の大きさ）が違うわけです。右はKL基底の誤差を1としたときの倍率で、フーリエは1.8〜2.9倍、区分定数は1.6〜1.8倍。どの $N$ でも1を下回る基底が存在しないという事実こそ、最適性の定理を実験で確かめたことになります。

全ての $N$ でKL基底が最小の誤差を与えています。 これが前節で証明した最適性の実験的な裏付けです。特に $N=32$ を見ると、KLの0.00317に対して区分定数は0.00520、フーリエは0.00646。区分定数基底の誤差はちょうど $1/(6N)$（$N=32$ なら0.00521）で、KLの $1/(\pi^2N)$ と比べると比が $6/\pi^2 = 0.608$。同じ精度を出すのに、KL基底なら約6割の本数で済む計算になります。フーリエ基底が振るわないのは、周期1の三角関数がどれも $t=0$ と $t=1$ で同じ値を取るのに対し、ブラウン運動は $W(0)=0$ で固定・$W(1)$ は自由という非対称な性質を持つためです。KL基底の $\sin((k-\frac12)\pi t)$ は、まさにこの非対称性を境界条件として取り込んでいます。

念のため、理論の要である「係数の無相関性」も直接確認しておきましょう。

K = 5
B = kl_basis(K)
coef = W @ B.T * dt2                       # c_k = <W, φ_k>
lam5 = np.array([1.0 / ((k - 0.5) * np.pi) ** 2 for k in range(1, K + 1)])

print("係数の分散（実測）:", np.round(coef.var(axis=0), 6))
print("固有値 λ_k（理論）:", np.round(lam5, 6))
print("\n正規化係数 ξ_k の相関行列:")
print(np.round(np.corrcoef((coef / np.sqrt(lam5)).T), 3))

係数の分散（実測）: [0.405309 0.044667 0.016124 0.008234 0.005032]
固有値 λ_k（理論）: [0.405285 0.045032 0.016211 0.008271 0.005004]

正規化係数 ξ_k の相関行列:
[[ 1.    -0.016  0.001 -0.012  0.017]
 [-0.016  1.    -0.007 -0.006 -0.006]
 [ 0.001 -0.007  1.     0.012 -0.002]
 [-0.012 -0.006  0.012  1.    -0.009]
 [ 0.017 -0.006 -0.002 -0.009  1.   ]]

左の対数軸の棒グラフでは、実測の分散（青）と理論固有値（赤）が5つのモードすべてでほぼ同じ高さに揃っています。右の相関行列は対角だけが濃い赤（$\rho=1$）で、非対角は色がほとんど付かない白に近い——つまり異なるモードの係数が互いにまったく相関していないことが視覚的にわかります。

係数 $c_k$ の分散が固有値 $\lambda_k$ と3桁レベルで一致し、$\xi_k=c_k/\sqrt{\lambda_k}$ の相関行列は対角成分が1、非対角成分がどれも $|\rho|\le0.017$ とほぼゼロになりました。$M=20000$ サンプルでの相関係数の標準誤差が $1/\sqrt{M}\approx0.007$ なので、この程度のばらつきは統計的な揺らぎの範囲に収まっています。「係数が無相関になるように基底を選んだら固有関数が出てきた」という出発点の要求が、実データ上でもきちんと満たされていることが確認できました。

数値実験がすべて理論と一致しました。最後に、実務でKL展開を使うときに気をつけたい点をまとめておきます。

KL展開の限界と実務上の注意

理論的に美しく最適でも、使う場面には向き不向きがあります。

（1）$R$ を知らないと始まらない。 KL展開は共分散関数を完全に知っている前提の理論です。実務では標本共分散で代用しますが、サンプル数 $M$ がグリッド点数 $n$ より少ないと標本共分散行列のランクが $M$ 未満になり、推定された固有関数は高次モードでほとんどノイズになります。$M \ll n$ の状況では、共分散関数に滑らかさを仮定した平滑化や、パラメトリックなカーネルへの当てはめが必要です。

（2）「最適」なのは二乗誤差についてだけ。 最適性の証明で最小化したのは平均二乗誤差でした。しかし用途によっては、二乗誤差が小さいことと「役に立つ」ことは違います。分類が目的なら、分散が小さくてもクラスを分ける方向のほうが重要になりえます（PCAとLDAの関係と同じ構図です）。KL展開は「分散を残す」ことに特化した圧縮であって、「情報を残す」ことを保証するわけではありません。

（3）基底がデータ依存で、計算コストがかかる。 KL基底は共分散ごとに違うので、圧縮側と復元側の両方が基底を共有する必要があります。$n\times n$ の固有分解は $O(n^3)$ で、$n$ が大きいと現実的でありません。画像圧縮でKLTではなくDCTが使われるのは、DCTが「多くの自然画像の共分散（1次マルコフ的な相関）に対してKLTの良い近似になり、しかも高速変換が使える」からです。理論的最適解より、少し劣るが安い近似が勝つ典型例です。

（4）非定常でも使えるのが強み。 一方で、KL展開はフーリエ解析と違って定常性を仮定しません。ブラウン運動の $R(s,t)=\min(s,t)$ は $s-t$ だけの関数ではない、完全な非定常過程です。それでもKL展開は成立し、しかも最適でした。定常過程しか扱えないスペクトル解析に対して、これはKL展開の大きな利点です。実際、成長曲線や過渡応答のような「時刻によって性質が変わる」データにこそKL展開は向いています。

（5）打ち切り誤差の見積もりは固有値だけで済む。 これは限界ではなく利点ですが、実務上とても効きます。データを再構成して誤差を測らなくても、固有値スペクトルを見れば必要な項数が決まります。UQで「入力ランダム場を何次元に落とすか」を決めるとき、この性質が設計を大きく楽にします。

まとめ

本記事では、カルーネン・レーベ展開の理論と導出を扱いました。

• 出発点は要求だった — 「直交展開の係数を無相関にしたい」という要求だけを課すと、基底は共分散作用素の固有関数でなければならないという結論が強制され、第二種フレドホルム積分方程式 $\int_0^T R(s,t)\varphi(s)ds=\lambda\varphi(t)$ が現れる
• マーサーの定理が土台 — 共分散作用素は対称・非負・コンパクトなので固有関数が正規直交完全系をなし、$R(s,t)=\sum_k\lambda_k\varphi_k(s)\varphi_k(t)$ が一様収束する。系として「固有値の総和＝総分散」が成り立つ
• 打ち切り誤差は残りの固有値の和 — $N$ 項打ち切りの平均二乗誤差はちょうど $\sum_{k>N}\lambda_k$ になる。誤差評価にデータは要らず、固有値のリストだけで済む
• 最適性は2通りに証明できる — ラグランジュ乗数法は「最適部分空間が共分散作用素の不変部分空間である」ことを、重み $d_j$ による直接評価は「$\sum_{k=1}^N\lambda_k$ を超えられない」ことを示す
• ブラウン運動は手で解ける — $\min(s,t)$ の積分方程式は2回微分すると $\lambda\varphi”=-\varphi$ に化け、境界条件 $\varphi(0)=0,\varphi'(1)=0$ から $\varphi_k=\sqrt2\sin((k-\frac12)\pi t)$、$\lambda_k=1/((k-\frac12)\pi)^2$ が得られる。第1モードだけで全分散の81%
• 数値実験で全て裏付けた — 離散固有分解は解析解を6桁再現し、モンテカルロで測った再構成誤差は理論値 $\sum_{k>N}\lambda_k$ と4桁一致、KL基底はフーリエ・区分定数基底より常に小さい誤差を与えた

KL展開は、「PCAは有限次元の話だ」という思い込みを解いてくれる定理です。固有値分解という一つの道具が、ベクトルから関数へ、そしてランダム場へと、そのまま持ち上がっていく。この構造が見えると、EOF解析もfPCAもKLTも確率有限要素法も、すべて同じ定理の別名だとわかります。

次のステップとして、以下の記事も参考にしてください。

• ブラウン運動（ウィーナー過程）の理論を解説 — この記事で共分散 $\min(s,t)$ を導いた過程そのものの定義と性質
• ブラウン運動の性質 — 連続性・非微分可能性・自己相似性 — 固有値が $k^{-2}$ でしか減らない理由の背景にある、標本路の粗さ
• 主成分分析（PCA）の理論と導出 — 本記事の離散版。分散最大化と誤差最小化の等価性を確認する
• 固有値分解の理論と応用を完全解説する — スペクトル分解の基礎。マーサーの定理はこの連続版
• 共分散行列とは？定義・性質・固有値分解からPython実装まで完全ガイド — 共分散作用素の有限次元版としての共分散行列

参考文献

• M. Loève, Probability Theory II, 4th ed., Springer, 1978（KL展開の原典的扱い）
• R. G. Ghanem, P. D. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer, 1991（UQへの応用の定番）
• I. Karatzas, S. E. Shreve, Brownian Motion and Stochastic Calculus, 2nd ed., Springer, 1991
• J. Ramsay, B. W. Silverman, Functional Data Analysis, 2nd ed., Springer, 2005（fPCAとしてのKL展開）