[Python] 最小二乗法によるデータの区分的多項式近似
目次
理論
多項式の最小二乗法
データ\(\{(x_n, y_n)\}_{n=1}^N\)をP次多項式\(f(x)=c^{(P)}x^P + c^{(P-1)}x^{P-1} + \cdots + c^{(1)}x + c^{(0)}\)で近似したい。
関数によるデータの近似の例
最小二乗法では、多項式の値とデータの値との差の2乗の和の最小化によって多項式近似を行う。
\begin{align} \mathrm{minimize} \quad \mathrm{Error} &= \frac{1}{2} \sum_{n=1}^N \left( y_n - f(x_n) \right)^2 \\ \end{align}
各係数\(c^{(p)}\)でErrorを偏微分した値が0になるという条件から解いていくことで、\(c^{(p)}\)を求めることができる。
\begin{align} 0 &= \frac{\partial \mathrm{Error}}{\partial c^{(p)}} \\ &= \frac{1}{2} \sum_{n=1}^N \frac{\partial}{\partial c^{(p)}} \left( y_n - f(x_n) \right)^2 \\ &= - \sum_{n=1}^N \left( y_n - f(x_n) \right) \frac{\partial f(x_n)}{\partial c^{(p)}} \\ &= - \sum_{n=1}^N \left( y_n - f(x_n) \right) \frac{\partial}{\partial c^{(p)}} \left( c^{(P)}x_n^P + c^{(P-1)}x_n^{P-1} + \cdots + c^{(1)}x_n + c^{(0)} \right) \\ &= - \sum_{n=1}^N \left( y_n - f(x_n) \right) x_n^p \\ &= - \sum_{n=1}^N y_nx_n^p + \sum_{n=1}^N x_n^pf(x_n) \\ &= - \sum_{n=1}^N y_nx_n^p + \sum_{n=1}^N x_n^p \sum_{q=0}^P c^{(q)}x_n^q \\ &= - \sum_{n=1}^N y_nx_n^p + \sum_{q=0}^P \left( \sum_{n=1}^N x_n^{p+q} \right) c^{(q)} \\ \end{align}
これをまとめると、
\begin{align} \begin{pmatrix} \sum_{n=1}^N y_nx_n^0 \\ \sum_{n=1}^N y_nx_n^1 \\ \vdots \\ \sum_{n=1}^N y_nx_n^p \\ \vdots \\ \sum_{n=1}^N y_nx_n^P \\ \end{pmatrix} = \begin{pmatrix} \sum_{n=1}^N x_n^{0} & \sum_{n=1}^N x_n^{1} & \cdots & \sum_{n=1}^N x_n^{q} & \cdots & \sum_{n=1}^N x_n^{P} \\ \sum_{n=1}^N x_n^{1} & \sum_{n=1}^N x_n^{2} & \cdots & \sum_{n=1}^N x_n^{1+q} & \cdots & \sum_{n=1}^N x_n^{1+P} \\ \vdots & \vdots & \ddots & \vdots & & \vdots \\ \sum_{n=1}^N x_n^{p} & \sum_{n=1}^N x_n^{p+1} & \cdots & \sum_{n=1}^N x_n^{p+q} & \cdots & \sum_{n=1}^N x_n^{p+P} \\ \vdots & \vdots & & \vdots & \ddots & \vdots \\ \sum_{n=1}^N x_n^{P} & \sum_{n=1}^N x_n^{P+1} & \cdots & \sum_{n=1}^N x_n^{P+q} & \cdots & \sum_{n=1}^N x_n^{2P} \\ \end{pmatrix} \cdot \begin{pmatrix} c^{(0)} \\ c^{(1)} \\ \vdots \\ c^{(q)} \\ \vdots \\ c^{(P)} \\ \end{pmatrix} \end{align}
となるので、逆行列を使って解くことで多項式の係数を求めることができる。
区分的多項式の最小二乗法
次の関数のように、複数の区間ごとに多項式で定義される関数を区分的多項式という。
特に、区間の節点で関数が滑らかに(=微分可能に)接続されている区分多項式はsplineと呼ぶこともある。
\begin{align} f(x) = \begin{cases} f_1(x) = c_1^{(P)}x^P + \cdots + c_1^{(1)}x + c_1^{(0)} \quad & (a_0 \leq x < a_1) \\ f_2(x) = c_2^{(P)}x^P + \cdots + c_2^{(1)}x + c_2^{(0)} \quad & (a_1 \leq x < a_2) \\ \quad \vdots \\ f_K(x) = c_k^{(P)}x^P + \cdots + c_k^{(1)}x + c_k^{(0)} \quad & (a_{k-1} \leq x < a_k) \\ \quad \vdots \\ f_K(x) = c_K^{(P)}x^P + \cdots + c_K^{(1)}x + c_K^{(0)} \quad & (a_{K-1} \leq x < a_K) \\ \end{cases} \end{align}
区分的多項式についても通常の多項式と同様の考え方で、データを近似する係数を求めることができる。
区間の個数を\(K\)とする。\(k\)番目の区間を\(I_k = [a_{k-1}, a_k]\)と表し、区間\(I_k\)内のデータを\(\{(x_{k,n}, y_{k,n})\}_{n=1}^{N_k}\)と表す。
各区間の多項式を\(f_k(x)=c_k^{(P)}x^P + c_k^{(P-1)}x^{P-1} + \cdots + c_k^{(1)}x + c_k^{(0)}\)とすると、
\begin{align} \begin{pmatrix} \sum_{n=1}^{N_1} y_{1,n}x_{1,n}^0 \\ \vdots \\ \sum_{n=1}^{N_1} y_{1,n}x_{1,n}^P \\ \vdots \\ \sum_{n=1}^{N_k} y_{k,n}x_{k,n}^0 \\ \vdots \\ \sum_{n=1}^{N_k} y_{k,n}x_{k,n}^P \\ \vdots \\ \sum_{n=1}^{N_K} y_{K,n}x_{K,n}^0 \\ \vdots \\ \sum_{n=1}^{N_K} y_{K,n}x_{K,n}^P \\ \end{pmatrix} = \begin{pmatrix} \sum_{n=1}^{N_1} x_{1,n}^{0} & \cdots & \sum_{n=1}^{{N_1}_1} x_{1,n}^{P} \\ \vdots & \ddots & \vdots \\ \sum_{n=1}^{N_1} x_{1,n}^{P} & \cdots & \sum_{n=1}^{N_1} x_{1,n}^{2P} \\ &&& \ddots \\ &&&& \sum_{n=1}^{N_k} x_{k,n}^{0} & \cdots & \sum_{n=1}^{N_k} x_{k,n}^{P} \\ &&&& \vdots & \ddots & \vdots \\ &&&& \sum_{n=1}^{N_k} x_{k,n}^{P} & \cdots & \sum_{n=1}^{N_k} x_{k,n}^{2P} \\ &&&&&&& \ddots \\ &&&&&&&& \sum_{n=1}^{N_K} x_{K,n}^{0} & \cdots & \sum_{n=1}^{N_K} x_{K,n}^{P} \\ &&&&&&&& \vdots & \ddots & \vdots \\ &&&&&&&& \sum_{n=1}^{N_K} x_{K,n}^{P} & \cdots & \sum_{n=1}^{N_K} x_{K,n}^{2P} \\ \end{pmatrix} \cdot \begin{pmatrix} c_1^{(0)} \\ \vdots \\ c_1^{(P)} \\ \vdots \\ c_k^{(0)} \\ \vdots \\ c_k^{(P)} \\ \vdots \\ c_K^{(0)} \\ \vdots \\ c_K^{(P)} \\ \end{pmatrix} \end{align}
となって係数を求めることができる。
条件付きの区分的多項式の最小二乗法
- 境界での関数の具体的な値
- 区間の節点で関数が連続
- 区間の節点で関数が滑らか
などの境界条件が課される場合を考える。
境界での関数の具体的な値
\(f(a_0) = b_0\)あるいは\(f(a_K) = b_K\)といった条件を課す場合、
\begin{align} b_0 &= f(a_0) = f_1(a_0) = c_1^{(P)}a_0^P + \cdots + c_1^{(1)}a_0 + c_1^{(0)} \\ b_K &= f(a_K) = f_K(a_K) = c_K^{(P)}a_K^P + \cdots + c_K^{(1)}a_K + c_K^{(0)} \\ \end{align}
とならなければならない。
Lagrangeの未定乗数法の考え方で、最小化するErrorに
\begin{align} \lambda_0 \left( c_1^{(P)}a_0^P + \cdots + c_1^{(1)}a_0 + c_1^{(0)} - b_0 \right) \\ \end{align}
あるいは
\begin{align} \lambda_K \left( c_K^{(P)}a_K^P + \cdots + c_K^{(1)}a_K + c_K^{(0)} - b_K \right) \\ \end{align}
の項を追加すれば良い。
\(\frac{\partial\mathrm{Error}}{\partial c_1^{(p)}}\)には\(a_0^pc_1^{(p)}\)が追加され、\(\frac{\partial\mathrm{Error}}{\partial c_K^{(p)}}\)には\(a_K^pc_K^{(p)}\)が追加される。
また、Errorを\(\lambda_0\)あるいは\(\lambda_K\)で偏微分することで、
\begin{align} b_0 &= c_1^{(P)}a_0^P + \cdots + c_1^{(1)}a_0 + c_1^{(0)} \\ b_K &= c_K^{(P)}a_K^P + \cdots + c_K^{(1)}a_K + c_K^{(0)} \\ \end{align}
が方程式に追加される。
区間の節点で関数が連続
区間の節点で関数が連続であるという条件を課す場合、
\begin{align} 0 &= f_k(a_k) - f_{k+1}(a_k) \\ &= c_k^{(P)}a_k^P + \cdots + c_k^{(1)}a_k + c_k^{(0)} - c_{k+1}^{(P)}a_k^P - \cdots - c_{k+1}^{(1)}a_k - c_{k+1}^{(0)} \\ \end{align}
とならなければならない。(\(k=1, 2, \cdots K-1\))
これも同様にLagrangeの未定乗数法を適用して、
\begin{align} \mu_k \left( c_k^{(P)}a_k^P + \cdots + c_k^{(1)}a_k + c_k^{(0)} - c_{k+1}^{(P)}a_k^P - \cdots - c_{k+1}^{(1)}a_k - c_{k+1}^{(0)} \right) \\ \end{align}
の項を追加する。
その結果、\(\frac{\partial\mathrm{Error}}{\partial c_k^{(p)}}\)には\(a_k^p\mu_k - a_{k-1}^p\mu_{k-1}\)が追加される。
区間の節点で関数が滑らか
区間の節点で関数が滑らかという条件を課す場合、
\begin{align} 0 &= f_k'(a_k) - f_{k+1}'(a_k) \\ &= Pc_k^{(P)}a_k^{P-1} + \cdots + c_k^{(1)} - Pc_{k+1}^{(P)}a_k^{P-1} - \cdots - c_{k+1}^{(1)} \\ \end{align}
とならなければならない。(\(k=1, 2, \cdots K-1\))
これも同様にLagrangeの未定乗数法を適用して、
\begin{align} \nu_k \left( Pc_k^{(P)}a_k^{P-1} + \cdots + c_k^{(1)} - Pc_{k+1}^{(P)}a_k^{P-1} - \cdots - c_{k+1}^{(1)} \right) \\ \end{align}
の項を追加する。
その結果、\(\frac{\partial\mathrm{Error}}{\partial c_k^{(p)}}\)には\(pa_k^{p-1}\nu_k - pa_{k-1}^{p-1}\nu_{k-1}\)が追加される。
全て適用する場合
これらの条件を全て適用する場合、連立方程式は次のようになる。
\begin{align} \begin{pmatrix} \sum_{n=1}^{N_1} y_{1,n}x_{1,n}^0 \\ \vdots \\ \sum_{n=1}^{N_1} y_{1,n}x_{1,n}^P \\ \vdots \\ \sum_{n=1}^{N_k} y_{k,n}x_{k,n}^0 \\ \vdots \\ \sum_{n=1}^{N_k} y_{k,n}x_{k,n}^P \\ \vdots \\ \sum_{n=1}^{N_K} y_{K,n}x_{K,n}^0 \\ \vdots \\ \sum_{n=1}^{N_K} y_{K,n}x_{K,n}^P \\ 0 \\ \vdots \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} = \begin{pmatrix} \sum_{n=1}^{N_1} x_{1,n}^{0} & \cdots & \sum_{n=1}^{{N_1}_1} x_{1,n}^{P} &&&&&&&&& 1 &&&&&& 0 &&&&&& 1 \\ \vdots & \ddots & \vdots &&&&&&&&& \vdots &&&&&& \vdots &&&&&& \vdots \\ \sum_{n=1}^{N_1} x_{1,n}^{P} & \cdots & \sum_{n=1}^{N_1} x_{1,n}^{2P} &&&&&&&&& a_1^P &&&&&& Pa_1^{P-1} &&&&&& a_0^P \\ &&& \ddots \\ &&&& \sum_{n=1}^{N_k} x_{k,n}^{0} & \cdots & \sum_{n=1}^{N_k} x_{k,n}^{P} &&&&&&& -1 & 1 &&&&& 0 & 0 \\ &&&& \vdots & \ddots & \vdots &&&&&&& \vdots & \vdots &&&&& \vdots & \vdots \\ &&&& \sum_{n=1}^{N_k} x_{k,n}^{P} & \cdots & \sum_{n=1}^{N_k} x_{k,n}^{2P} &&&&&&& -a_{k-1}^P & a_k^P &&&&& -Pa_{k-1}^{P-1} & Pa_k^{P-1} \\ &&&&&&& \ddots \\ &&&&&&&& \sum_{n=1}^{N_K} x_{K,n}^{0} & \cdots & \sum_{n=1}^{N_K} x_{K,n}^{P} &&&&&& -1 &&&&&& 0 && 1 \\ &&&&&&&& \vdots & \ddots & \vdots &&&&&& \vdots &&&&&& \vdots && \vdots \\ &&&&&&&& \sum_{n=1}^{N_K} x_{K,n}^{P} & \cdots & \sum_{n=1}^{N_K} x_{K,n}^{2P} &&&&&& -a_{K-1}^P &&&&&& -Pa_{K-1}^{P-1} && a_K^P \\ 1 & \cdots & a_1^P \\ \\ &&&& -1 & \cdots & -a_{k-1}^P \\ &&&& 1 & \cdots & a_k^P \\ \\ &&&&&&&& -1 & \cdots & -a_{K-1}^P \\ 0 & \cdots & Pa_1^{P-1} \\ \\ &&&& 0 & \cdots & -Pa_{k-1}^{P-1} \\ &&&& 0 & \cdots & Pa_k^{P-1} \\ \\ &&&&&&&& 0 & \cdots & -Pa_{K-1}^{P-1} \\ 1 & \cdots & a_0^P \\ &&&&&&&& 1 & \cdots & a_K^P \\ \end{pmatrix} \cdot \begin{pmatrix} c_1^{(0)} \\ \vdots \\ c_1^{(P)} \\ \vdots \\ c_k^{(0)} \\ \vdots \\ c_k^{(P)} \\ \vdots \\ c_K^{(0)} \\ \vdots \\ c_K^{(P)} \\ \mu_1 \\ \vdots \\ \mu_{k-1} \\ \mu_k \\ \vdots \\ \mu_{K-1} \\ \nu_2 \\ \vdots \\ \nu_{k-1} \\ \nu_k \\ \vdots \\ \nu_{K-1} \\ \lambda_0 \\ \lambda_K \\ \end{pmatrix} \end{align}
逆行列によって解くことで、条件を満たした上でErrorを最小化する区分的多項式の係数を求めることができる。
ソースコード
言語はPython。必要なパッケージはNumPy。
Data
各区間に含まれるデータを扱うクラス。
class Data(list):
def sum(self, f):
s = 0
for (x, y) in self:
s += f(x, y)
return s
@staticmethod
def data_split(data_list, interval):
K = len(interval) - 1
d = [Data() for _ in range(K)]
for data in data_list:
x = data[0]
appended = False
for k in range(K):
if interval[k] <= x <= interval[k+1]:
d[k].append(data)
appended = True
if not appended:
print(f'data {data} is excluded.')
return d
make_mv
係数計算のための行列とベクトルを作る関数。
- data - (x, y)のリストを指定する。
- interval - 区間の端を指定する。例えばinterval=[0, 0.4, 0.7, 1]とすると、[0, 0.4], [0.4, 0.7], [0.7, 1]の3つの区間で多項式近似を行う。
- order - 多項式の次数。
- smooth_condition - Trueのとき区間の端で滑らかに繋がる(微分が一致する)。
- boundary_condition - 区間全体の左右の端における値を指定する。不要な場合はNoneを与える。片側だけ指定する場合は[y, None]のようにする。
def make_mv(data, interval, order, smooth_condition=True, boundary_condition=[0, 0]):
interval = sorted(interval)
K = len(interval) - 1
data = Data.data_split(data, interval)
size = K*(order+1) + K - 1
if smooth_condition:
size += K - 1
if boundary_condition is not None:
size += sum(1 for c in boundary_condition if c is not None)
m = np.zeros((size, size))
v = np.zeros(size)
for k in range(K):
d = data[k]
sumxp = [d.sum(lambda x, y: x**p) for p in range(2*order+1)]
i1 = K*(order+1) + k
i2 = i1 + K - 1
for p1 in range(order+1):
offset = k*(order+1)
for p2 in range(order+1):
m[offset + p1, offset + p2] = sumxp[p1 + p2]
v[offset + p1] = d.sum(lambda x, y: y*x**p1)
if k <= K-2:
ap = interval[k+1]**p1
m[i1, offset + p1 ] = ap
m[i1, offset + p1 + K - 1] = -ap
m[offset + p1, i1] = ap
m[offset + p1 + K - 1, i1] = -ap
if smooth_condition:
pap = p1*interval[k+1]**(p1-1) if p1 != 0 else 0
m[i2, offset + p1 ] = pap
m[i2, offset + p1 + K - 1] = -pap
m[offset + p1, i2] = pap
m[offset + p1 + K - 1, i2] = -pap
if boundary_condition is not None:
v1, v2 = boundary_condition
if v1 is not None:
i = -2 if v2 is not None else -1
v[i] = v1
for p1 in range(order+1):
ap = interval[0]**p1
m[i, p1] = ap
m[p1, i] = ap
if v2 is not None:
j = (K-1)*(order+1)
v[-1] = v2
for p1 in range(order+1):
ap = interval[-1]**p1
m[-1, j + p1] = ap
m[j + p1, -1] = ap
return m, v
fit
データを近似する多項式の係数を求める。
引数はmake_mvと同様。
def fit(data, interval, order, smooth_condition=True, boundary_condition=[0, 0]):
K = len(interval) - 1
m, v = make_mv(data, interval, order, smooth_condition, boundary_condition)
coeffs = np.linalg.inv(m)@v
return coeffs[:K*(order+1)].reshape([K, -1])
validate
近似の結果を評価。
実際のデータと近似関数とのMSE, L∞距離を測定する。
def validate(data, interval, coeffs):
interval = sorted(interval)
K = len(interval) - 1
def f(x, coeff):
y = 0
for i, c in enumerate(coeff):
y += c*x**i
return y
max_error = 0
sum_error = 0
n = 0
for x, y in data:
for k in range(K):
if interval[k] < x < interval[k+1] or (k == K-1 and x == interval[-1]):
error = (y - f(x, coeffs[k]))**2
if max_error < error:
max_error = error
sum_error += error
n += 1
return sum_error/n, np.sqrt(max_error)
ソースコード全体
import numpy as np
class Data(list):
def sum(self, f):
s = 0
for (x, y) in self:
s += f(x, y)
return s
@staticmethod
def data_split(data_list, interval):
K = len(interval) - 1
d = [Data() for _ in range(K)]
for data in data_list:
x = data[0]
appended = False
for k in range(K):
if interval[k] <= x <= interval[k+1]:
d[k].append(data)
appended = True
if not appended:
print(f'data {data} is excluded.')
return d
def make_mv(data, interval, order, smooth_condition=True, boundary_condition=[0, 0]):
interval = sorted(interval)
K = len(interval) - 1
data = Data.data_split(data, interval)
size = K*(order+1) + K - 1
if smooth_condition:
size += K - 1
if boundary_condition is not None:
size += sum(1 for c in boundary_condition if c is not None)
m = np.zeros((size, size))
v = np.zeros(size)
for k in range(K):
d = data[k]
sumxp = [d.sum(lambda x, y: x**p) for p in range(2*order+1)]
i1 = K*(order+1) + k
i2 = i1 + K - 1
for p1 in range(order+1):
offset = k*(order+1)
for p2 in range(order+1):
m[offset + p1, offset + p2] = sumxp[p1 + p2]
v[offset + p1] = d.sum(lambda x, y: y*x**p1)
if k <= K-2:
ap = interval[k+1]**p1
m[i1, offset + p1 ] = ap
m[i1, offset + p1 + K - 1] = -ap
m[offset + p1, i1] = ap
m[offset + p1 + K - 1, i1] = -ap
if smooth_condition:
pap = p1*interval[k+1]**(p1-1) if p1 != 0 else 0
m[i2, offset + p1 ] = pap
m[i2, offset + p1 + K - 1] = -pap
m[offset + p1, i2] = pap
m[offset + p1 + K - 1, i2] = -pap
if boundary_condition is not None:
v1, v2 = boundary_condition
if v1 is not None:
i = -2 if v2 is not None else -1
v[i] = v1
for p1 in range(order+1):
ap = interval[0]**p1
m[i, p1] = ap
m[p1, i] = ap
if v2 is not None:
j = (K-1)*(order+1)
v[-1] = v2
for p1 in range(order+1):
ap = interval[-1]**p1
m[-1, j + p1] = ap
m[j + p1, -1] = ap
return m, v
def fit(data, interval, order, smooth_condition=True, boundary_condition=[0, 0]):
K = len(interval) - 1
m, v = make_mv(data, interval, order, smooth_condition, boundary_condition)
coeffs = np.linalg.inv(m)@v
return coeffs[:K*(order+1)].reshape([K, -1])
def validate(data, interval, coeffs):
interval = sorted(interval)
K = len(interval) - 1
def f(x, coeff):
y = 0
for i, c in enumerate(coeff):
y += c*x**i
return y
max_error = 0
sum_error = 0
n = 0
for x, y in data:
for k in range(K):
if interval[k] < x < interval[k+1] or (k == K-1 and x == interval[-1]):
error = (y - f(x, coeffs[k]))**2
if max_error < error:
max_error = error
sum_error += error
n += 1
return sum_error/n, np.sqrt(max_error)
if __name__ == '__main__':
interval = [0, 0.35, 0.62, 1] # example
data = np.loadtxt('input.txt')
coeffs = fit(data, interval, 5,
smooth_condition=False,
boundary_condition=[0, 0])
mse, max_error = validate(data, interval, coeffs)
print(f'mse = {mse}, max error = {max_error}')
np.savetxt('output.txt', coeffs)