はじめに
このページには，2019年10月24日の分離工学の課題の解説などが書いてあります。授業で扱った内容の解説は行っていません。
なお，非線形最小二乗法についての記述は，かなり簡略化してあるものの，求められているであろうレベルを遥かに上回ったものになっています。 理解できなくても嘆く必要はありません(というよりも端折りすぎていて理解しようにもできない; 興味がある人は自分で調べてください)。 線形の最小二乗法くらいは理解しておきましょう。
最終的な成果物だけが欲しい方はExcelファイルをダウンロードできます。
最小二乗法について
最小二乗法は，パラメータ$\v{b}=\rb{b_1, b_2, \cdots, b_m}^\mathrm{T}$を持つ函数$f_\v{b}(x)$について，既知の独立変数の組$X=\wb{x_1, x_2, \cdots, x_n}$の$f$による写像$f(X)$が， 既知のデータの組$Y=\wb{y_1, y_2, \cdots, y_n}$と最もよく合うように最適なパラメータ$\v{b}$を求める方法である。 ここで，「既知のデータの組と最もよく合う」とは，残差平方和 \[ \varepsilon = \sum_i (y_i - f_\v{b}(x_i))^2 \]が最小となることを言う。各残差の2乗をとっているのは，正負の異なる残差が打ち消しあうのを防ぐためである。
残差平方和$\varepsilon$の定義のうち，$x_i$と$y_i$は既知の値であるので，変更できる値は函数$f$のパラメータ$\v{b}$のみである。 したがって，$\varepsilon$は$\v{b}$の函数として$\varepsilon(\v{b})$と書くことができる。このとき，問題は「函数$\varepsilon$が最小値をとる$\v{b}$を求める」となる。
線形最小二乗法
一般に，函数$\varepsilon(\v{b})$が最小値を取る点は，連立方程式(正規方程式) \begin{align*} \frac{\partial\varepsilon}{\partial b_1} &= 0 \\ \frac{\partial\varepsilon}{\partial b_2} &= 0 \\ &\vdots \\ \frac{\partial\varepsilon}{\partial b_m} &= 0 \end{align*} を解くことで求められる。
例えば，モデル函数が$y=ax+b$であるとき，$\displaystyle\varepsilon=\sum_i(ax_i+b-y_i)^2$であるので， \begin{align*} \frac{\partial\varepsilon}{\partial a} &= \sum_i2(ax_i+b-y_i) \times x_i = 0 \\ \frac{\partial\varepsilon}{\partial b} &= \sum_i2(ax_i+b-y_i) \times 1 = 0 \end{align*} を解くことで最適な$a, b$を求めることができる。これを整理すると \[ \begin{pmatrix} \sum x_i^2 & \sum x_i \\ \sum x_i & \sum \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} \sum x_iy_i \\ \sum y_i \end{pmatrix} \]となるので，求める$a, b$は \begin{align*} \begin{pmatrix} a \\ b \end{pmatrix} &= \begin{pmatrix} \sum x_i^2 & \sum x_i \\ \sum x_i & n \end{pmatrix}^{-1}\begin{pmatrix} \sum x_iy_i \\ \sum y_i \end{pmatrix} \\ &= \frac1{n\sum x_i^2-\rb{\sum x_i}^2}\begin{pmatrix} n & -\sum x_i \\ -\sum x_i & \sum x_i^2 \end{pmatrix}\begin{pmatrix} \sum x_iy_i \\ \sum y_i \end{pmatrix} \\ &= \frac1{n\sum x_i^2-\rb{\sum x_i}^2}\begin{pmatrix} n\sum x_iy_i - \sum x_i\sum y_i \\ \sum x_i^2\sum y_i - \sum x_i\sum x_iy_i \end{pmatrix} \end{align*} となる。
非線形最小二乗法
モデル函数が簡単な場合は正規方程式を解析的に解くことが可能であるが，一般には正規方程式を解析的に解くことは困難である。 正規方程式による解法が使用できない場合には，数値計算によりパラメータを求めることになる。
数値計算により函数の最小値を求める方法としては様々な手法が提案されているが，ここではそのうちの一部の概要を紹介する。
最急降下法
最急降下法は，初期パラメータ$\v{b}_0$を \[ \v{b}_{i+1} = \v{b}_i - t_i\nabla\varepsilon(\v{b}_i) \] という反復計算により更新することで$\v{b}$を得るアルゴリズムである。ここで，$t_i$は1回の反復で更新する数値の重みを決定するパラメータである。
$\nabla\varepsilon$は$\varepsilon$の勾配ベクトルであり，勾配ベクトルの逆向きに進むときに函数値が最も大きく減少することを利用したアルゴリズムである。
Gauss–Newton法
Gauss–Newton法は，初期パラメータ$\v{b}_0$を \[ \v{b}_i = \v{b}_i - J(\v{b_i})^{-1}\varepsilon(\v{b}_i) \] という反復計算により更新することで$\v{b}$を得るアルゴリズムである。ここで，$J$はヤコビアンである。(逆ベクトル函数のTaylor展開から導出可能であるが詳細は省略。)
Gauss–Newton法は求める解の近傍でのみ適用可能であるため，適切な初期推測値を与える必要がある。
Levenberg–Marquardt法
最急降下法ではほとんどの場合に初期パラメータによらず目的の解に収束するのに対し，Gauss–Newton法では初期推測値が解から離れている場合に安定性が悪い。 一方，収束性の点から見ると，最急降下法では1次の収束性しか持たないのに対し，Gauss–Newton法は解の近傍で2次の収束性を持つ。 ここで，収束が遅いが安定である最急降下法と，収束が速いが(ただし解の付近でという制約がある)安定性が悪いGauss–Newton法を組み合わせるという方法が考えられる。 すなわち，始めは最急降下法により反復を行い，解の近傍でGauss–Newton法に切り替えるというアルゴリズムである。
これを実現するのがLevenberg–Marquardt法(LMA)であり，反復計算は \[ \v{b}_{i+1} = \v{b}_i - (J_i^\mathrm{T}+\lambda_iI)^{-1}J_i^\mathrm{T}\varepsilon(\v{b}_i) \] となる。ここで，$\lambda_i\geq0$はダンピング項であり，初期値は大きく取り，解に近づくにつれ小さくしていく。
Langmuirの吸着等温式
LMAにより，与えられたデータを基に函数の \[ N_A = \frac{a^\prime N_Sp}{1+a^\prime p} \] パラメータを求めると，$N_S=5.076228235424373$，$a^\prime=0.0009678519065705363$となり，自由度調整済み決定係数は$R_f^2=.999872$となる。このパラメータを基に$N_A$をプロットしたものは
Langmuirプロット
Langmuirの吸着等温式はそのままでは解析的にフィッティングを行うことができないが，両辺の逆数をとることで線形に変形することができる。すなわち，$N_A$についての式では， \[ \frac1{N_A} = \frac1{a^\prime N_S}\cdot\frac1p + \frac1{N_S} \] とすることで，独立変数$1/p$，従属変数$1/N_A$の1次函数とすることができる。(逆数をとるという発想は，Michaelis–Menten式におけるLineweaver–Burkプロットと同じ。) この式をプロットすると
LMAによる回帰よりも決定係数が小さくなった原因としては
- 逆数をとる操作などを含むため丸め誤差の影響を受けた
- $p=0$の点を無視したため(ゼロ除算となるため)，参照しているデータが少ない
などが考えられる。($R^2>.999$であるため近似の精度としては悪くはない; 相対的な話として)