📂통계적분석

STLSQ이란?

알고리즘

STLSQ

풀 랭크인 행렬 $X \in \mathbb{R}^{m \times p}$ 와 $Y \in \mathbb{R}^{m \times n}$ 에 대해 행렬방정식 $$ Y = X B $$ 이 주어졌을 때, 가능한 한 스파스한 $B \in \mathbb{R}^{p \times n}$ 을 찾는 스파스 회귀 문제를 생각해보자. STLSQ^{sequentially Thresholded Least-Squares Algorithm}은 신디^sINDy 알고리즘을 제안한 브런턴^brunton의 논문에서 사용된 메소드으로써, 이름 그대로 순차적으로 최소제곱법을 반복하며 각 성분의 절대값이 주어진 쓰레숄딩 $\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)$ 에 대해 수행될 때, 우리는 그것을 신디라 부른다¹. 어찌보면 STLSQ 자체가 신디를 위해서 나왔다고 할 수 있어서 조금 아이러니하지만, 알고리즘의 측면에서 STLSQ는 동역학계와 전혀 상관 없이 그 자체로 작동할 수 있어 ‘신디’라는 이름 대신 언급되는 경우도 제법 있다.

알고리즘의 의사코드를 읽어보면 알 수 있듯, STLSQ는 일종의 그리디 알고리즘이다. 알고리즘의 수렴성이 증명된 것²과는 별개로 회귀분석의 통계적 측면에서는 까놓고 말해 아무런 의미가 없지만, 구현하기 쉽고 수렴속도가 매우 빠르면서 어쨌거나 스파스 회귀를 수행한다는 점에서는 의심의 여지 없이 쓸만한 방법이라 할 수 있다.

코드

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

파이썬

브런턴 교수의 저서⁴에 있는 파이썬 코드다.

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

내가 사용하는 줄리아 코드다.

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