📂統計的分析

STLSQとは何か？

アルゴリズム

STLSQ

フルランクの行列 $X \in \mathbb{R}^{m \times p}$ と $Y \in \mathbb{R}^{m \times n}$ が与えられたとき、割と スパース な $B \in \mathbb{R}^{p \times n}$ を見つけるスパース回帰問題を考えよう。STLSQ^{sequentially Thresholded Least-Squares Algorithm} は、「SINDy」Algorithm を提案したブルントンが使っていた方法で、名前の通り、 最小二乗法 を繰り返し利用し、各成分の 絶対値 が与えられたしきい値 $\lambda$ より下がったら $0$ として扱い、その列(独立変数) を排除する仕組みになっている。

擬似コード

共通して、括弧で囲まれた上付きのテキスト $\cdot ^{(k)}$ は $k$ 番目を、 $[p] = \left\{ 1 , \cdots , p \right\} \subset \mathbb{N}$ は $1$ から $p$ までの 自然数 の 集合 を示す。 $X^{\dagger}$ は行列 $X$ の 擬似逆行列 であり、最小二乗法 によって得られる。インデックス集合 $S$ に対して、 $\mathbf{x}_{S}$ はベクトル $\mathbf{x}$ の該当インデックスだけを残したベクトル、 $X_{S}$ は行列 $X$ の該当列だけを残した行列である。

説明

特に STLSQ が 動的系 $\dot{\mathbf{x}} = f \left( \mathbf{x} \right)$ の微分係数から作られた行列 $Y = \dot{X}$ とライブラリ行列 $X = \Theta (X)$ に対して行われる場合、我々はそれを「SINDy」と呼ぶ¹。一見するとSTLSQ自体がSINDyのために出たと言えるが、アルゴリズムの観点からはSTLSQ自体が動的系と全く関係なくそのまま作動できるため、「SINDy」という名前の代わりに言及されることもままある。

アルゴリズムの擬似コードを読めば分かるように、STLSQは一種の グリーディアルゴリズム だ。アルゴリズムの収束性が証明されたこと²とは別に、回帰分析の統計的側面から言えば正直、意味がないが、実装が簡単で収束速度が非常に速く、何はともあれ スパース回帰 を行うという点では疑いようのない有用な方法と言える。

コード

Matlab

以下は、ブルントン教授が自分の研究で使っている実際の Matlab コードだ³。

%% compute Sparse regression: sequential least squares
Xi = Theta\dXdt; % initial guess: Least-squares

% lambda is our sparsification knob.
for k=1:10
    smallinds = (abs(Xi)<lambda); % find small coefficients
    Xi(smallinds)=0; % and threshold
    for ind = 1:n % n is state dimension
        biginds = ˜smallinds(:,ind);
        % Regress dynamics onto remaining terms to find sparse Xi
        Xi(biginds,ind) = Theta(:,biginds)\dXdt(:,ind);
    end
end

Python

ブルントン教授の著書⁴ にある Python コードだ。

def sparsifyDynamics(Theta,dXdt,lamb,n):
    # Initial guess: Least-squares
    Xi = np.linalg.lstsq(Theta,dXdt,rcond=None)[0]

    for k in range(10):
        smallinds = np.abs(Xi) < lamb # Find small coeffs.
        Xi[smallinds]=0 # and threshold
        for ind in range(n): # n is state dimension
            biginds = smallinds[:,ind] == 0
            # Regress onto remaining terms to find sparse Xi
            Xi[biginds,ind] = np.linalg.lstsq(Theta[:,
                biginds],dXdt[:,ind],rcond=None)[0]
        return Xi

Julia

私が使っている Julia コードだ。

using SparseArrays

function STLSQ(Θ, dXdt, λ = 10^(-6))
    if Θ isa AbstractSparseMatrix
        Θ = Matrix(Θ)
    end
    Ξ = Θ \ dXdt
    dim = size(Ξ, 2)
    __🚫 = 0
    
    while true
        print(".")
        🚫 = abs.(Ξ) .< λ
        Ξ[🚫] .= 0
        for j in 1:dim
            i_ = .!🚫[:, j]
            Ξ[i_, j] = Θ[:,i_] \ dXdt[:,j]
        end
        if __🚫 == 🚫 println("Stopped!"); break end # Early stopping
        __🚫 = deepcopy(🚫)
    end
    println("MSE: ", sum(abs2, Θ * Ξ - dXdt) / length(dXdt))

    return Ξ
end