#はじめに

Exvel VBAで最小自乗スプライン法を使ってフィッティングする関数です。
他の言語では既存のコードがあったりライブラリがあったりするらしいのですが、
Excel VBAでは見当たらなかったので書いてみました。
基本的な流れは下のリンク先を参考にしています。
ここでは細かい計算は載せていないので、気になる方は参照ください。

スタンフォード大学の学生が学ぶ制約付き最小二乗法とその応用
 数値流体力学大全第4章スプライン関数

#最小自乗スプライン法とは

そもそもこの呼び方が正式名称なのかも分かりませんが、とりあえずここではそう呼びます。
例えば、次のようなデータがあるとします。
sin関数に乱数を含めて適当に生成したデータです。

このデータをフィッティングしたい。しかし、多項式関数ではうまく行かない。
そういうときに最小自乗スプライン法が役に立ちます。
1つの関数で表現できないなら、複数の関数をつなげてやればいいのです。

このように、広い領域で曲線をうまくフィッティングできています。
この例では元データを5つの領域に分けて、各領域を3次関数でフィッティングしています。
つまり、ある領域iを表す関数は

S_i(\alpha) = \sum _{n=0}^3 a_{i,n}\alpha^n\quad(0 < \alpha < 1)

となっています。
境界には、関数が連続になるような境界条件が入っています。
各領域の幅は別に同じでなくてもいいですが、ここでは簡単のため同じとします。

#最小自乗スプライン法のコード

詳しいことは一旦後回しにして、まずはコードを載せておきます。
今回は、実験データと条件を入力すると各領域の係数が返ってくるFunctionの形にしています。
入力は、

ni(領域番号)：各領域番号のデータ数
nv：領域の数
min：最小値
max：最大値
x(領域番号，データ番号)：xのデータ
y(領域番号，データ番号)：yのデータ

出力は、

spl(領域番号，0~3)：各領域の係数．

となっています。

Function spl(ni As Variant, nv As Long, min As Double, max As Double, x As Variant, y As Variant) As Double()
    
    '入力データ
    'ni                         各領域内にデータがいくつあるか
    'nv                         領域を何分割するか
    'min                        フィッティングする領域の最小値
    'max                        フィッティングする領域の最大値
    'x(領域番号，データ番号)    各領域内にあるデータのx座標
    'y(領域番号，データ番号)    各領域内にあるデータのy座標
    
    '関数の定義
    Dim i As Long, j As Long, k As Long     'ループ計算用整数
    Dim l As Long, m As Long                'ループ計算用整数
    Dim temp1 As Double, temp2 As Double    '計算で一時的に値を入れる用の変数
    Dim dx As Double                        '領域の幅
    Dim Ai() As Double                      '最小自乗条件（Ai・ui=bi）のAi
    Dim bi() As Double                      '最小自乗条件（Ai・ui=bi）のbi
    Dim Ci(2, 8) As Double                  '境界条件の行列
    Dim A() As Double                       '各領域の最小自乗の行列を含む行列
    Dim AA() As Double                      'KKT行列の左上部分
    Dim b() As Double                       '各領域の最小自乗の右辺を含む右辺
    Dim C() As Double                       '各領域の境界条件を含む行列
    Dim KKT() As Double                     'KKT行列
    Dim KKT2() As Double                    '並び替え後のKKT行列
    Dim u() As Double                       'KKT条件の解
    Dim v() As Double                       'KKT条件の右辺
    Dim v2() As Double                      '並び替え後のKKT条件の右辺
    Dim nk As Long                          'KKT行列のサイズ
    Dim a_KKT() As Double                   'KKT条件の解
    
    '計算用定数の計算
    dx = (max - min) / nv                   '各領域の幅（今回は全領域を等分している）
    nk = 4 * nv + 2 * (nv - 1)              'KKT行列のサイズ（最小自乗部分+境界条件部分）
    
    '配列の再定義
    ReDim Ai(nv, 4, 4), bi(nv, 4)
    ReDim A(4 * nv, 4 * nv), AA(4 * nv, 4 * nv), b(4 * nv), C(2 * (nv - 1), 4 * nv)
    ReDim KKT(nk, nk), KKT2(nk, nk), u(nk), v(nk), v2(nk)
    ReDim a_KKT(nv, 4)
    
    '各領域での最小自乗行列の生成
    For j = 0 To nv - 1
        For k = 0 To 3
            For l = 0 To 3
                Ai(j, k, l) = 0
                For i = 0 To ni(j) - 1
                    Ai(j, k, l) = Ai(j, k, l) + x(j, i) ^ (k + l)
                Next i
            Next l
            bi(j, k) = 0
            For i = 0 To ni(j) - 1
                bi(j, k) = bi(j, k) + y(j, i) * x(j, i) ^ (k)
            Next i
        Next k
    Next j
    
    '境界条件
    Ci(0, 0) = 1
    Ci(0, 1) = 1
    Ci(0, 2) = 1
    Ci(0, 3) = 1
    Ci(0, 4) = -1
    Ci(0, 5) = 0
    Ci(0, 6) = 0
    Ci(0, 7) = 0
    Ci(1, 0) = 0
    Ci(1, 1) = 1
    Ci(1, 2) = 2
    Ci(1, 3) = 3
    Ci(1, 4) = 0
    Ci(1, 5) = -1
    Ci(1, 6) = 0
    Ci(1, 7) = 0
    
    '各領域の最小自乗条件を含む行列とKKT行列の左上の計算
    For k = 0 To 4 * nv - 1
        For l = 0 To 4 * nv - 1
            A(k, l) = 0
            AA(k, l) = 0
        Next l
    Next k
    For j = 0 To nv - 1
        For k = 0 To 3
            For l = 0 To 3
                A(4 * j + k, 4 * j + l) = Ai(j, k, l)
                For m = 0 To 3
                    AA(4 * j + k, 4 * j + l) = AA(4 * j + k, 4 * j + l) + 2 * Ai(j, k, m) * Ai(j, m, l)
                Next m
            Next l
            b(4 * j + k) = bi(j, k)
        Next k
    Next j
    
    '各領域の境界条件を含む行列の計算
    For k = 0 To 2 * (nv - 1) - 1
        For l = 0 To 4 * nv - 1
            C(k, l) = 0
        Next l
    Next k
    For j = 0 To nv - 1 - 1
        For k = 0 To 1
            For l = 0 To 7
                C(2 * j + k, 4 * j + l) = Ci(k, l)
            Next l
        Next k
    Next j
    
    'KKT行列の計算
    For k = 0 To nk - 1
        For l = 0 To nk - 1
            KKT(k, l) = 0
        Next l
    Next k
    For k = 0 To 4 * nv - 1
        For l = 0 To 4 * nv - 1
            KKT(k, l) = AA(k, l)
        Next l
    Next k
    For k = 0 To 2 * (nv - 1) - 1
        For l = 0 To 4 * nv - 1
            KKT(4 * nv + k, l) = C(k, l)
            KKT(l, 4 * nv + k) = C(k, l)
        Next l
    Next k
    
    'KKT条件の右辺の計算
    For k = 0 To nk - 1
        u(k) = 0
        v(k) = 0
    Next k
    For k = 0 To 4 * nv - 1
        For l = 0 To 4 * nv - 1
            v(k) = v(k) + 2 * A(k, l) * b(l)
        Next l
    Next k
    
    '行列の対角項が非0になるよう並び替え
    For k = 0 To nk - 1
        For l = 0 To nk - 1
            KKT2(k, l) = KKT(k, l)
        Next l
        v2(k) = v(k)
    Next k
    For j = 0 To nv - 2
        For k = 0 To nk - 1
            KKT2(4 * j, k) = KKT(4 * nv + 2 * j, k)
            KKT2(4 * j + 1, k) = KKT(4 * nv + 2 * j + 1, k)
            KKT2(4 * nv + 2 * j, k) = KKT(4 * j, k)
            KKT2(4 * nv + 2 * j + 1, k) = KKT(4 * j + 1, k)
        Next k
        v2(4 * j) = v(4 * nv + 2 * j)
        v2(4 * j + 1) = v(4 * nv + 2 * j + 1)
        v2(4 * nv + 2 * j) = v(4 * j)
        v2(4 * nv + 2 * j + 1) = v(4 * j + 1)
    Next j
    
    'ガウスの消去法でKKT条件を解く
    For k = 0 To nk - 1
        For l = k + 1 To nk - 1
            temp1 = KKT2(l, k) / KKT2(k, k)
            For m = k To nk - 1
                KKT2(l, m) = KKT2(l, m) - temp1 * KKT2(k, m)
            Next m
            v2(l) = v2(l) - temp1 * v2(k)
        Next l
    Next k
    u(nk - 1) = v2(nk - 1) / KKT2(nk - 1, nk - 1)
    For k = nk - 2 To 0 Step -1
        For l = k + 1 To nk - 1
            v2(k) = v2(k) - KKT2(k, l) * u(l)
        Next l
        u(k) = v2(k) / KKT2(k, k)
    Next k
    
    '解を整理
    m = 0
    For j = 0 To nv - 1
        For k = 0 To 3
            m = m + 1
            a_KKT(j, k) = u(m - 1)
        Next
    Next j
    
    '値を返す
    spl = a_KKT()

End Function

#コードの中身について

ここから、コードの中身をポイントを絞って説明していきます。
数学的な話をするとそれだけで長くなってしまうので、できるだけ簡潔に。

##解くべき式

このコードは何の式を解いているの？の部分です。
普通の3次の最小二乗法であれば、

\begin{pmatrix}
\sum 1 & \sum x_i & \sum x_i^2 & \sum x_i^3 \\
\sum x_i & \sum x_i^2 & \sum x_i^3 & \sum x_i^4 \\
\sum x_i^2 & \sum x_i^3 & \sum x_i^4 & \sum x_i^5 \\
\sum x_i^3 & \sum x_i^4 & \sum x_i^5 & \sum x_i^6 
\end{pmatrix}

\begin{pmatrix}
a_{i,0} \\
a_{i,1} \\
a_{i,2} \\
a_{i,3} 
\end{pmatrix}

=

\begin{pmatrix}
\sum y_i \\
\sum y_ix_i \\
\sum y_ix_i^2 \\
\sum y_ix_i^3 
\end{pmatrix}

\Leftrightarrow

A_i\vec{a_i}=\vec{b_i}

の行列式を解くことでフィッティングが可能です。
Aの逆行列を作って、左からかければ解けます。
複数の領域（iとj）でフィッティングするなら、これを並べてやればいいのです。

\begin{pmatrix}
A_i & 0 \\
0 & A_j
\end{pmatrix}

\begin{pmatrix}
\vec{a_i} \\
\vec{a_j} 
\end{pmatrix}

=

\begin{pmatrix}
\vec{b_i} \\
\vec{b_j} 
\end{pmatrix}


\Leftrightarrow

A\vec{a}=\vec{b}

バラせばそれぞれは普通の最小二乗法の形になっています。
しかし、このままでは領域の境界でうまくつながってくれません。
そこで、Lagrangeの未定乗数法を用いて、ある条件を満たしながら最小二乗法を行います。
細かい計算はすっ飛ばして、最終形はこちらです。

\begin{pmatrix}
2AA^T & C^T \\
C & 0
\end{pmatrix}

\begin{pmatrix}
\vec{a} \\
0
\end{pmatrix}

=

\begin{pmatrix}
2A\vec{b} \\
0
\end{pmatrix}

これをKKT条件、左辺の行列をKKT行列と呼ぶそうです。
行列中のCが境界条件です。領域i,jについて

\begin{pmatrix}
1 & 1 & 1 & 1 & -1 & 0 & 0 & 0 \\
0 & 1 & 2 & 3 & 0 & -1 & 0 & 0
\end{pmatrix}

\begin{pmatrix}
\vec{a_i} \\
\vec{a_j}
\end{pmatrix}

=

\begin{pmatrix}
0 \\
0
\end{pmatrix}

という境界条件を課しています。
1行目が0次の、2行目が1次の連続境界条件です。
各領域の分割幅を同じにしているので、少し簡単な形になっています。
コードではこの式を解いて各領域の関数の係数を計算しています。
KKT条件の式を作るまでがコードの前半です。
（'KKT条件の右辺の計算　まで）

##行列を解く

ガウスの消去法を使って行列を解こうと思ったのですが、1つ問題があります。
KKT行列の右下の部分が0であり、つまり対角項に0を含んでいます。
このままでは解けませんが、自動のピボット選択を入れるのも面倒です。
そこで、境界条件行列Cの転置の部分を使って、
非0の成分が対角項に来るように並べ替えます。
掃き出し法とかでよくやるアレです。
そのままのKKT行列が

\begin{pmatrix}
k_{i,00} & k_{i,01} & k_{i,02} & k_{i,03} & 0 & 0 & 0 & 0 & 1 & 0 \\
k_{i,10} & k_{i,11} & k_{i,12} & k_{i,13} & 0 & 0 & 0 & 0 & 1 & 1 \\
k_{i,20} & k_{i,21} & k_{i,22} & k_{i,23} & 0 & 0 & 0 & 0 & 1 & 2 \\
k_{i,30} & k_{i,31} & k_{i,32} & k_{i,33} & 0 & 0 & 0 & 0 & 1 & 3 \\
0 & 0 & 0 & 0 & k_{j,00} & k_{j,01} & k_{j,02} & k_{j,03} & -1 & 0 \\
0 & 0 & 0 & 0 & k_{j,10} & k_{j,11} & k_{j,12} & k_{j,13} & 0 & -1 \\
0 & 0 & 0 & 0 & k_{j,20} & k_{j,21} & k_{j,22} & k_{j,23} & 0 & 0 \\
0 & 0 & 0 & 0 & k_{j,30} & k_{j,31} & k_{j,32} & k_{j,33} & 0 & 0 \\
1 & 1 & 1 & 1 & -1 & 0 & 0 & 0 & 0 & 0 \\
0 & 1 & 2 & 3 & 0 & -1 & 0 & 0 & 0 & 0
\end{pmatrix}

とすれば、入れ替えた後は、

\begin{pmatrix}
1 & 1 & 1 & 1 & -1 & 0 & 0 & 0 & 0 & 0 \\
0 & 1 & 2 & 3 & 0 & -1 & 0 & 0 & 0 & 0 \\
k_{i,20} & k_{i,21} & k_{i,22} & k_{i,23} & 0 & 0 & 0 & 0 & 1 & 2 \\
k_{i,30} & k_{i,31} & k_{i,32} & k_{i,33} & 0 & 0 & 0 & 0 & 1 & 3 \\
0 & 0 & 0 & 0 & k_{j,00} & k_{j,01} & k_{j,02} & k_{j,03} & -1 & 0 \\
0 & 0 & 0 & 0 & k_{j,10} & k_{j,11} & k_{j,12} & k_{j,13} & 0 & -1 \\
0 & 0 & 0 & 0 & k_{j,20} & k_{j,21} & k_{j,22} & k_{j,23} & 0 & 0 \\
0 & 0 & 0 & 0 & k_{j,30} & k_{j,31} & k_{j,32} & k_{j,33} & 0 & 0 \\
k_{i,00} & k_{i,01} & k_{i,02} & k_{i,03} & 0 & 0 & 0 & 0 & 1 & 0 \\
k_{i,10} & k_{i,11} & k_{i,12} & k_{i,13} & 0 & 0 & 0 & 0 & 1 & 1
\end{pmatrix}

となっています。
行列が大きくなっても、これを繰り返せば対角項がすべて非0になります。
ここまで来たら、後は解くだけですね。

#使用例

最後に、簡単に使用例を紹介して終わりにします。
冒頭で出した例のデータを生成し、スプライン関数でフィッティングするマクロです。

Sub スプライン近似マクロ使用例()
    
    '関数の定義
    Dim i As Long, j As Long, k As Long
    Dim nd As Long, nv As Long, ni() As Long
    Dim min As Double, max As Double
    Dim x() As Double, y() As Double, a() As Double, dx As Double
    Dim temp1 As Double, temp2 As Double, temp3 As Double
    Dim cnt As Long
    
    'データ数,最小値,最大値,分割数,幅の入力
    nd = 201
    min = 0
    max = 20
    nv = 5
    dx = (max - min) / nv
    
    '配列の再定義
    ReDim ni(nv)
    ReDim x(nv, nd), y(nv, nd)
    
    '初期値入力
    For j = 0 To nv - 1
        ni(j) = 0
    Next j
    
    'データの生成
    For i = 0 To nd - 1
        temp3 = i / 10
        For j = 0 To nv - 1
            temp1 = min + dx * j
            temp2 = min + dx * (j + 1)
            If temp3 >= temp1 And temp3 < temp2 Then
                ni(j) = ni(j) + 1
                x(j, ni(j) - 1) = (temp3 - temp1) / dx
                y(j, ni(j) - 1) = Sin(temp3) + temp3 + (Rnd - 0.5)
            End If
            If j = nv - 1 And temp3 = temp2 Then
                ni(j) = ni(j) + 1
                x(j, ni(j) - 1) = (temp3 - temp1) / dx
                y(j, ni(j) - 1) = Sin(temp3) + temp3 + (Rnd - 0.5)
            End If
        Next j
    Next i
    
    'スプライン関数による係数の計算
    a = spl(ni(), nv, min, max, x(), y())
    
    'データの出力
    cnt = 0
    For j = 0 To nv - 1
        For i = 1 To ni(j)
            cnt = cnt + 1
            temp1 = 0
            For k = 0 To 3
                temp1 = temp1 + a(j, k) * x(j, i - 1) ^ k
            Next k
            Cells(cnt + 1, 1) = (cnt - 1) / 10
            Cells(cnt + 1, 2) = y(j, i - 1)
            Cells(cnt + 1, 3) = temp1
        Next i
    Next j
    
    'ラベルの出力
    Cells(1, 1) = "x"
    Cells(1, 2) = "y"
    Cells(1, 3) = "y_spl"
    
End Sub

追記(20220612)

pythonで書いた場合のメモ


#data:2×nのリスト、(x,y)
#bc:境界値、入力例（[0.0, 10.0, 15.0, 20.0]）

import numpy as np

def LSM_spl_a(data, bc):
    #ndarrayに変換
    dt = np.array(data)
    #初期設定
    bcs = len(bc) - 1
    dx = np.zeros(bcs)
    tmp_sum = np.zeros(7)
    Ai = np.zeros((bcs, 4, 4))
    A = np.zeros((4 * bcs, 4 * bcs))
    di = np.zeros((bcs, 4))
    d = np.zeros(4 * bcs + 2 * (bcs - 1))
    Ci = np.zeros((bcs - 1, 2, 8))
    C = np.zeros((2 * (bcs - 1), 4 * bcs))
    C_tpl = np.array([[1, 1, 1, 1, -1, 0, 0, 0],
                      [0, 0, 0, 0, 0, 0, 0, 0]])
    KKT = np.zeros((4 * bcs + 2 * (bcs - 1), 4 * bcs + 2 * (bcs - 1)))
    #最小二乗行列計算
    for i in range(bcs):
        dx[i] = bc[i + 1] - bc[i] 
        tmp_dt = dt[((dt[:, 0] <= bc[i + 1]) & (dt[:, 0] >= bc[i]))]
        dt = dt[~((dt[:, 0] <= bc[i + 1]) & (dt[:, 0] >= bc[i]))]
        tmp_dt[:, 0] = [(row[0] - bc[i]) / dx[i] for row in tmp_dt]
        for j in range(7):
            tmp_sum[j] = np.sum(tmp_dt[:,0] ** j)
        for r in range(4):
            for c in range(4):
                Ai[i, r, c] = tmp_sum[r + c]
        for j in range(4):
            di[i, j] = np.sum((tmp_dt[:,0] ** j) * tmp_dt[:,1])
        A[i * 4 : (i * 4) + 4, i * 4 : (i * 4) + 4] = Ai[i, :, :]
        d[i * 4 : (i * 4) + 4] = 2 * np.dot(Ai[i, :, :], di[i, :])
    #境界条件行列計算
    for i in range(bcs - 1):
        Ci[i, :, :] = C_tpl
        Ci[i, 1, 1] = 1 / dx[i]
        Ci[i, 1, 2] = 2 / dx[i]
        Ci[i, 1, 3] = 3 / dx[i]
        Ci[i, 1, 5] = -1 / dx[i + 1]
        C[i * 2 : (i * 2) + 2, i * 4 : (i * 4) + 8] = Ci[i, :, :]
    #KKT行列計算
    KKT[0 : 4 * bcs, 0 : 4 * bcs] = 2 * np.dot(A, A)
    KKT[4 * bcs : 4 * bcs + 2 * (bcs - 1), 0 : 4 * bcs] = C
    KKT[0 : 4 * bcs, 4 * bcs : 4 * bcs + 2 * (bcs - 1)] = C.T
    KKT_inv = np.linalg.inv(KKT)
    a = np.dot(KKT_inv, d)
    a = np.delete(a, slice(4 * bcs, 4 * bcs + 2 * (bcs - 1)), 0)
    a = a.reshape(bcs, 4)
    return a

def LSM_spl_d(data, a, bc):
    #初期設定
    bcs = len(bc) - 1
    dx = np.zeros(bcs)
    for i in range(bcs):
        dx[i] = bc[i + 1] - bc[i] 
    cal = np.zeros(len(data))
    #スプライン計算
    for k in range(len(data)):
        for i in range(bcs):
            if bc[i] <= data[k] <= bc[i + 1]:
                x = (data[k] - bc[i]) / dx[i]
                for j in range(4):
                    cal[k] += a[i, j] * (x ** j)
                break
    return cal

def LSM_spl_x(x, a, bc):
    #初期設定
    bcs = len(bc) - 1
    dx = np.zeros(bcs)
    for i in range(bcs):
        dx[i] = bc[i + 1] - bc[i] 
    #スプライン計算
    cal = 0.0
    for i in range(bcs):
        if bc[i] <= x <= bc[i + 1]:
            x = (x - bc[i]) / dx[i]
            for j in range(4):
                cal += a[i, j] * (x ** j)
            return cal
    return cal