Seto's Page

Switch to

Last update:

最小二乗法と格子定数精密化

最小二乗法の概要

　観測方程式とは、実験によって観測される量 (観測量) と未知のパラメータとの関係を表す条件方程式のことです。例えば未知のパラメータを$ a^1, a^2, \cdots, a^m$ とし、観測量を $Y$としたとき、$$
f(a^1, a^2, \cdots, a^m) = Y$$ のような関係式が観測方程式ということになります。なお、この等式が成り立つのは全く誤差のない完璧な観測ができたときだけです。もし観測方程式が線形的であれば、次のような行列のかけ算で表現することができます。$$
f(a^1, a^2, \cdots, a^m) =
\begin{pmatrix}X^1 & X^2 & \cdots & X^m \end{pmatrix}
\begin{pmatrix}a^1 \\ a^2 \\ \vdots \\ a^m \end{pmatrix}
$$ここで $ X^1, X^2, \cdots, X^m$ は観測方程式中の係数です。念のため、このページで出てくる上付き文字は乗数ではなく、パラメータの個数を表しています。

　もし観測量 $Y$ のばらつきが正規分布に従っており、かつ観測方程式が線形であれば、最小二乗法を適用することで、最も尤(もっと)もらしい未知の値を推定することができます。以下に例を挙げて説明します。
　$n$ 個の観測の結果、次のような観測方程式の係数($X$)および観測量 ($Y$)が得られたとします。$$
(X^1_1, X^2_1, \cdots, X^m_1, Y_1), (X^1_2, X^2_2, \cdots, X^m_2, Y_2), \cdots,(X^1_n, X^2_n, \cdots, X^m_n, Y_n)$$これらを使って、未知のパラメータ $ a^1, a^2, \cdots, a^m$ を推定するのが目標です。まず以下のように記号を定義しましょう。$$
\mathbf{a} =
\begin{pmatrix} a^1 \\ a^2 \\ \vdots \\ a^m \end{pmatrix}
, \,\,\,
\mathbf{X}=
\begin{pmatrix}
X^1_1 & X^2_1 & \ldots & X^m_1 \\
X^1_2 & X^2_2 & \ldots & X^m_2 \\
\vdots & \vdots & \ddots & \vdots \\
X^1_n & X^2_n & \ldots & X^m_n
\end{pmatrix}
, \,\,\,
\mathbf{Y}=\begin{pmatrix}Y_1 \\ Y_2 \\ \vdots \\ Y_n \end{pmatrix}
$$また、すべての観測が等しい精度であるとは限りませんので、重み行列$\mathbf{W}$ (重み [誤差の二乗の逆数] を対角成分に持つ対角行列)を以下のように定義します。
$$\mathbf{W}=
\begin{pmatrix}
w_1 & 0 & \ldots & 0 \\
0 & w_2 & \ldots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \ldots & w_n \\
\end{pmatrix}
$$各観測における、「観測方程式から得られる量」と「実際の観測量」との差は、以下のような列ベクトルで表現することができます。$$
\begin{pmatrix}
f_1(a_1,a_2,\cdots,a_m)-Y_1\\
f_2(a_1,a_2,\cdots,a_m)-Y_2 \\
\vdots \\
f_n(a_1,a_2,\cdots,a_m) -Y_n
\end{pmatrix}
=\begin{pmatrix}
X^1_1 & X^2_1 & \ldots & X^m_1 \\
X^1_2 & X^2_2 & \ldots & X^m_2 \\
\vdots & \vdots & \ddots & \vdots \\
X^1_n & X^2_n & \ldots & X^m_n
\end{pmatrix}
\begin{pmatrix} a^1 \\ a^2 \\ \vdots \\ a^m \end{pmatrix}- \mathbf{Y}
= \mathbf{X}\mathbf{a}-\mathbf{Y}
$$この列ベクトルの各成分をそれぞれ二乗して重みを掛け、さらにそれらの総和をとった量が「重み付き残差の二乗の和」ということになります。これを$\Delta$と書くことにして、行列演算で表現すると、$$
\Delta=(\mathbf{a}^{tr} \, \mathbf{X}^{tr} – \mathbf{Y}^{tr})\, \mathbf{W} \, (\mathbf{X}\, \mathbf{a} -\mathbf{Y})
$$となります¹。なお上付き$tr$は転置行列を意味しています。

　さて、$\Delta$を最小にするような $\mathbf{a}$ が何なのかを知りたいわけですが、ここで簡単な微分積分の定理を思い出しましょう。$\Delta$が最小であるということは、「 $\mathbf{a}$ の微小変化 $\delta\mathbf{a}$ に対して、$\Delta$の応答する量 $\delta \Delta$がゼロである」、ということに他なりません。これを数式で表現すると、$$
\delta \Delta = \delta \mathbf{a}^{tr} \, \mathbf{X}^{tr} \, \mathbf{W}\, (\mathbf{X} \, \mathbf{a} – \mathbf{Y}) + (\mathbf{a}^{tr} \, \mathbf{X}^{tr} – \mathbf{Y}^{tr}) \, \mathbf{W} \, \mathbf{X}\, \delta \mathbf{a}
= 2\delta\mathbf{a}^{tr} \, (\mathbf{X}^{tr} \, \mathbf{W} \, \mathbf{X} \, \mathbf{a} – \mathbf{X}^{tr} \, \mathbf{W} \, \mathbf{Y}) = 0$$となります。この式を変形すると、$$
\mathbf{a} = (\mathbf{X}^{tr} \, \mathbf{W} \, \mathbf{X})^{-1} \, \mathbf{X}^{tr}\, \mathbf{W}\, \mathbf{Y} $$という関係が得られます。これが未知パラメータ $\mathbf{a}$ を推定する方程式です。もし重み行列 $ \mathbf{W}$ の対角成分をすべて1とみなしてよければ (つまりすべて観測の誤差が等しければ)、$$
\mathbf{a} = (\mathbf{X}^{tr} \, \mathbf{X})^{-1} \, \mathbf{X}^{tr} \, \mathbf{Y} $$という式になります。

格子定数の最適化

　なんらかの回折実験をして、面指数と面間隔の値をたくさん得たとしましょう。このデータから格子定数を計算することを考えます。格子定数の計算は、結晶系によって観測方程式が変わります。たとえば、三斜晶系であれば、
$$
\frac{1}{d^2}=
\begin{pmatrix}h^2 & k^2 & l^2 & hk & kl & lh \end{pmatrix}
\begin{pmatrix}{a^*}^2 \\ {b^*}^2 \\ {c^*}^2 \\ b^*c^*\cos\alpha^* \\ c^*a^*\cos\beta^* \\ a^*b^*\cos\gamma^* \end{pmatrix}
$$のようになります。ここで、$a^*, b^*$ などは未知パラメータであるこの結晶の逆格子定数示しており、$h,k,l$ は観測方程式中の係数と考えることができます。観測方程式が分かれば、あとは上で説明したとおりの手順で逆格子定数を求めればいいのですが、少し悩ましいのは誤差行列をどのように取り扱うかということです。

　まず、簡単のため散乱角度 $\theta$ のみが誤差を含むと考えることにしましょう (実際には様々な要因が関与します)。$\theta$ の変化に対する$1/d^2$ の応答は、ブラッグ条件から$$1/d^2 = 4\sin^2\theta/\lambda^2$$であるので、$$\partial{(1/d^2)}/\partial{\theta} = 4\sin 2\theta/\lambda^2$$となります。つまり、$\theta$ が $\partial{\theta}$ 変化すると、 $1/d^2$ は $\partial{\theta}\times 4\sin 2\theta/\lambda^2$ だけ変化するということです。ということで、とりあえず観測値 $1/d^2$ に対する重み行列対角成分は、$\sin 2\theta$ の二乗の逆数である $1/\sin^{2}2\theta$ということにしておきましょう。すなわち、
$$\mathbf{W’}=
\begin{pmatrix}
1/\sin^{2}2\theta_1 & 0 & \ldots & 0 \\
0 & 1/\sin^{2}2\theta_2 & \ldots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \ldots & 1/\sin^{2}2\theta_n \\
\end{pmatrix}$$ とします。ただし、 $1/\sin^{2}2\theta$ は各 $1/d^2$ の分散の逆数そのものではなく、それを何倍かした量です。本当の誤差行列 $\mathbf{W}$は、 $$
\mathbf{W}= z \, \mathbf{W’}
$$のように、$\mathbf{W’}$に適当な実数 $z$ を掛けた量です。

　ところが、本当の誤差行列ではないこの $\mathbf{W’}$ を上述の式に当てはめ、$\mathbf{a} = (\mathbf{X}^{tr} \, \mathbf{W’} \, \mathbf{X})^{-1} \, \mathbf{X}^{tr}\, \mathbf{W’}\, \mathbf{Y}$ のように計算しても、問題なく $\mathbf{a}$ の最適値を求めることができます。それは、重み行列 $\mathbf{W’}$が2回作用することによって、うまい具合に $z$ の効果がキャンセルアウトされるからです。うまくできてますね。

　さてこれで一応パラメータ $ \mathbf{a}$は求まりました。ここで誤差について改めて考えてみましょう。誤差の伝播式を考えると、もし正しい $ \mathbf{W}$ を使っていれば、行列$ (\mathbf{X}^{tr} \mathbf{W} \mathbf{X})^{-1} $ の対角成分が各パラメータの分散になるはずです。しかしながら、上で述べたように $\mathbf{W’}$ は割合しか決めていないのでその絶対値が分かりません。ここで、少し目先を変えて、分散の定義を思い出してみましょう。すなわち$$
N – P = \sum^{N}_{i=0}\frac{\delta_i^2}{s_i}
$$ (ただし$N$ : データ個数, $P$ : パラメータ数, $\delta$ : データの偏差, $s$ : データの分散) という関係です。この式を用いて、最適化したパラメータから以下の量を計算すると、なんと各パラメータの分散の絶対値がわかってしまいます。$$
s_i = \frac{\displaystyle\sum^{N}_{i=0}\delta_i^2w_i}{w_i (N-P)}
$$この分散の平方根が誤差になります。なお、これらはあくまで $\mathbf{a}$ (すなわち逆格子定数)の誤差なので、格子定数の誤差に戻すためにはさらに誤差の伝播を解いていく必要がありますが、根気さえあればできるのでここでは省略します。

列ベクトルである$\mathbf{X}\, \mathbf{a} -\mathbf{Y}$を行ベクトルに変換するため、$(\mathbf{X}\, \mathbf{a} -\mathbf{Y})^{tr} = \mathbf{a}^{tr} \, \mathbf{X}^{tr} – \mathbf{Y}^{tr}$ という関係を使っています。 ↩︎