非線性動(dòng)力學(xué)的稀疏識(shí)別(SINDy)簡(jiǎn)介

社區(qū)小助手 2024-08-16 17:06:51

在科學(xué)和工程領(lǐng)域，理解和預(yù)測(cè)動(dòng)態(tài)系統(tǒng)的行為至關(guān)重要。然而，很多時(shí)候我們并不完全知道描述這些系統(tǒng)的精確數(shù)學(xué)模型。近年來(lái)，一種名為SINDy（Sparse Identification of Nonlinear Dynamics，非線性動(dòng)力學(xué)的稀疏識(shí)別）的系統(tǒng)辨識(shí)方法受到了廣泛關(guān)注，它能夠從數(shù)據(jù)中自動(dòng)發(fā)現(xiàn)動(dòng)態(tài)系統(tǒng)的數(shù)學(xué)模型。

SINDy的基本原理

SINDy方法的核心思想是利用稀疏回歸技術(shù)從測(cè)量數(shù)據(jù)中識(shí)別出動(dòng)態(tài)系統(tǒng)的非線性微分方程。簡(jiǎn)單來(lái)說(shuō)，它假設(shè)系統(tǒng)的動(dòng)態(tài)可以由一個(gè)函數(shù)庫(kù)中的少數(shù)幾個(gè)函數(shù)的線性組合來(lái)描述。通過(guò)優(yōu)化算法，SINDy能夠選擇出對(duì)系統(tǒng)動(dòng)態(tài)影響最大的幾個(gè)函數(shù)，并確定它們的系數(shù)，從而構(gòu)建出一個(gè)簡(jiǎn)潔而精確的數(shù)學(xué)模型。

一個(gè)簡(jiǎn)單的例子

假設(shè)我們有一個(gè)簡(jiǎn)單的一維動(dòng)態(tài)系統(tǒng)，其動(dòng)態(tài)方程為：

$\dot{x} = -2x + x^3$

為了模擬這個(gè)系統(tǒng)，我們生成一些模擬數(shù)據(jù)，并使用SINDy方法來(lái)恢復(fù)這個(gè)方程。

步驟1：生成模擬數(shù)據(jù)

首先，我們使用北太天元來(lái)生成模擬數(shù)據(jù)。假設(shè)初始條件為(x(0) = 1)，我們使用歐拉方法或更高級(jí)的數(shù)值方法來(lái)求解這個(gè)微分方程，并收集一段時(shí)間內(nèi)的數(shù)據(jù)。

% 模擬數(shù)據(jù)生成  dt = 0.1; 
% 時(shí)間步長(zhǎng)  t = 0:dt:2; 
% 時(shí)間向量  t = t&#39;; x = zeros(size(t)); 
% 初始化狀態(tài)向量  x(1) = 0.6; 
% 初始條件  for i = 1:length(t)-1      x(i+1) = x(i) + dt*(-2*x(i)  +  x(i).^3 ); 
% 歐拉方法更新狀態(tài)  end  dxdt =  diff(x) ./diff(t); % 計(jì)算速度

步驟2：構(gòu)建函數(shù)庫(kù)

接下來(lái)，我們構(gòu)建一個(gè)函數(shù)庫(kù)，其中包含可能描述系統(tǒng)動(dòng)態(tài)的候選函數(shù)。在這個(gè)例子中，我們可以選擇多項(xiàng)式函數(shù)作為候選函數(shù)，例如(x), (x^2), (x^3), ...等。

% 構(gòu)建函數(shù)庫(kù)  x = x(1:end-1);
 Theta = [x,  x.^2, x.^3]; % 候選函數(shù)庫(kù)

步驟3：稀疏回歸

最后，我們使用稀疏回歸方法來(lái)確定哪些函數(shù)對(duì)系統(tǒng)動(dòng)態(tài)有顯著影響。在這個(gè)例子中，我們可以使用LASSO（Least Absolute Shrinkage and Selection Operator）回歸，它是一種能夠產(chǎn)生稀疏解的線性回歸方法。

針對(duì)這個(gè)具體的例子，LASSO回歸的優(yōu)化問(wèn)題可以表示為：

$\min_{\beta} \left\{ \frac{1}{2N} \sum_{i=1}^{N} (y_i - \beta_1 x_i - \beta_2 x_i^2 - \beta_3 x_i^3)^2 + \lambda \sum_{j=1}^{3} |\beta_j| \right\}$

使用迭代軟閾值算法(ISTA)求解LASSO問(wèn)題，其數(shù)學(xué)原理可以詳細(xì)闡述如下：

LASSO問(wèn)題的定義

LASSO（Least Absolute Shrinkage and Selection Operator）是一種回歸分析方法，它通過(guò)構(gòu)造一個(gè)懲罰函數(shù)來(lái)得到一個(gè)較為精煉的模型。LASSO的基本思想是在回歸系數(shù)的絕對(duì)值之和小于一個(gè)常數(shù)的約束條件下，使殘差平方和最小化，從而能夠產(chǎn)生某些嚴(yán)格等于0的回歸系數(shù)，得到可以解釋的模型。其數(shù)學(xué)形式通常表示為最小化以下目標(biāo)函數(shù)：

$\min_{\beta} \left( \frac{1}{2} \| A\beta - b \|_2^2 + \lambda \| \beta \|_1 \right)$

其中，A 是設(shè)計(jì)矩陣，b 是觀測(cè)向量，\beta 是待求的回歸系數(shù)向量，λ 是正則化參數(shù)，控制懲罰項(xiàng)的影響。

2. 迭代軟閾值算法(ISTA)的原理

迭代軟閾值算法(ISTA)是用于求解LASSO問(wèn)題的一種有效方法。它通過(guò)迭代的方式逐步逼近最優(yōu)解。算法的基本步驟如下：

初始化

* 設(shè)定初始估計(jì)值 $\beta^0$ ，通?？梢栽O(shè)為0向量或隨機(jī)向量。

* 選擇合適的步長(zhǎng) t_k 和正則化參數(shù) λ 。

迭代過(guò)程

對(duì)于每一次迭代 $k = 0, 1, 2, \ldots$ ，執(zhí)行以下步驟：

1. 梯度下降：首先，對(duì)當(dāng)前估計(jì)值 $\beta^k$ 進(jìn)行梯度下降，以減小殘差平方和的部分。這可以通過(guò)計(jì)算 $A^T A \beta^k - b$ 并按一定步長(zhǎng) $t_k$ 更新 β 來(lái)實(shí)現(xiàn)。

2. 軟閾值處理：然后，對(duì)梯度下降后的結(jié)果進(jìn)行軟閾值處理，以實(shí)現(xiàn)L1正則化的效果。軟閾值函數(shù)的形式為：

$\text{SoftThreshold}_{\lambda}(x) = \text{sgn}(x) \cdot (\max(|x| - \lambda, 0))$

其中，$\text{sgn}(x)$ 是符號(hào)函數(shù)，返回 $x$ 的符號(hào)。這個(gè)軟閾值函數(shù)會(huì)將絕對(duì)值小于 λ 的系數(shù)縮減到0，而對(duì)大于λ 的系數(shù)進(jìn)行一定程度的縮減。

3. 更新估計(jì)值：將軟閾值處理后的結(jié)果作為下一次迭代的起始點(diǎn)，即

$\beta^{k+1} = \text{SoftThreshold}_{\mu}(\text{梯度下降后的結(jié)果})$ 。

終止條件

* 當(dāng) $\| \beta^{k+1} - \beta^k \|_2$ 小于某個(gè)預(yù)設(shè)的閾值時(shí)，或者達(dá)到最大迭代次數(shù)時(shí)，算法停止。

3. 算法收斂性

在一定的條件下，ISTA算法可以保證收斂到LASSO問(wèn)題的最優(yōu)解。收斂速度取決于步長(zhǎng) $t_k$ 的選擇和問(wèn)題的具體性質(zhì)。通常，為了加速收斂，可以采用一些優(yōu)化策略，如Nesterov加速或采用更復(fù)雜的步長(zhǎng)選擇策略。這里簡(jiǎn)單實(shí)現(xiàn)一個(gè)線性搜索確定步長(zhǎng)方法

function [alpha] = linesearch(Theta, y, beta, grad, lambda)  % 線性搜索尋找最佳步長(zhǎng)        
% 初始化參數(shù)      alpha = 1; % 初始步長(zhǎng)      
c = 1e-4;  % Armijo條件中的常數(shù)      
rho = 0.5; % 步長(zhǎng)縮減因子      
max_ls_iter = 50; % 線性搜索的最大迭代次數(shù)            
% 計(jì)算梯度下降的初始方向      d = -grad;            
for ls_iter = 1:max_ls_iter          % 嘗試新的步長(zhǎng)          
beta_new = beta + alpha * d;                    % 應(yīng)用軟閾值操作          beta_new = sign(beta_new) .* max(0, abs(beta_new) - lambda);                    
% 計(jì)算目標(biāo)函數(shù)的值          f_new = 0.5 * norm(Theta * beta_new - y)^2 + lambda * norm(beta_new, 1);         
 f_old = 0.5 * norm(Theta * beta - y)^2 + lambda * norm(beta, 1);                    % 檢查Armijo條件是否滿足         
  if f_new &lt;= f_old + c * alpha * grad&#39; * d              break; % 如果滿足，則退出線性搜索          
  end                    
  % 如果不滿足，則縮減步長(zhǎng)并繼續(xù)搜索          alpha = rho * alpha;      
  end            
  % 如果達(dá)到最大迭代次數(shù)仍未找到滿足條件的步長(zhǎng)，則可能需要調(diào)整參數(shù)或檢查算法實(shí)現(xiàn)      
  if ls_iter == max_ls_iter          warning(&#39;線性搜索達(dá)到最大迭代次數(shù)，未找到滿足條件的步長(zhǎng)。&#39;);      
  end  end

北太天元中暫時(shí)沒(méi)有內(nèi)置的LASSO函數(shù)，可以用m函數(shù)實(shí)現(xiàn)一個(gè)簡(jiǎn)單的版本。

function [beta] = simpleLassoSolver(Theta, y, lambda, maxIter, tol)  % simpleLassoSolver: 
使用迭代軟閾值算法(ISTA)求解LASSO問(wèn)題  % 輸入:  
%   Theta - 設(shè)計(jì)矩陣 (n x p)  %   y - 響應(yīng)向量 (n x 1)  %   lambda - LASSO正則化參數(shù)  %   maxIter - 最大迭代次數(shù)  %   tol - 收斂容忍度  % 輸出:  
%   beta - 回歸系數(shù) (p x 1)    % 初始化  [n, p] = size(Theta);  beta = zeros(p, 1);  
% 初始化解  x_old = zeros(p, 1); % 上一步的迭代解    for iter = 1:maxIter      % 梯度下降步驟      grad = Theta&#39; * (Theta * beta - y) ;        
% 使用線性搜索確定最佳步長(zhǎng)      alpha = linesearch(Theta, y, beta, grad, lambda);      
% 更新解      x_new = beta - alpha * grad ;          
% 應(yīng)用軟閾值操作      beta = sign(x_new) .* max(0, abs(x_new) - lambda);          
% 檢查收斂性      fprintf([&#39;iter = &#39;, num2str(iter), &#39; 
beta = &#39;, num2str(beta(1)) &#39;, &#39; num2str(beta(2)),  &#39;, &#39;, 
num2str(beta(3)),  &#39;  norm = &#39; num2str( norm( Theta*beta - y) ), newline]);    
 if  norm( beta -x_old ) &lt; tol  &amp;&amp; norm ( Theta*beta - y) &lt; tol        break;      
end      x_old = beta; 
end    
end

% 稀疏回歸  % 注意：在實(shí)際應(yīng)用中，你可能需要對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理  % LASSO回歸  lambda = 1e-9; 
maxIter = 100000;tol = 1e-9; 
[coef] = simpleLassoSolver(Theta, dxdt, lambda, maxIter, tol) norm( Theta*coef - dxdt)

結(jié)果分析

通過(guò)檢查回歸系數(shù)coef，我們可以發(fā)現(xiàn)哪些函數(shù)對(duì)系統(tǒng)動(dòng)態(tài)有顯著影響。在這個(gè)例子中，我們希望恢復(fù)出原始的微分方程(\dot{x} = -2x + x^3)，因此我們應(yīng)該期望在回歸系數(shù)中看到-2和1對(duì)應(yīng)于(x)和(x^3)的系數(shù)，而其他函數(shù)的系數(shù)應(yīng)該接近于0。

結(jié)論

SINDy是一種強(qiáng)大的方法，能夠從數(shù)據(jù)中自動(dòng)發(fā)現(xiàn)動(dòng)態(tài)系統(tǒng)的數(shù)學(xué)模型。通過(guò)結(jié)合稀疏回歸技術(shù)和適當(dāng)?shù)暮瘮?shù)庫(kù)，SINDy可以準(zhǔn)確地識(shí)別出描述系統(tǒng)動(dòng)態(tài)的關(guān)鍵函數(shù)和它們的系數(shù)。這種方法在科學(xué)和工程領(lǐng)域具有廣泛的應(yīng)用前景，特別是在對(duì)復(fù)雜系統(tǒng)進(jìn)行建模和分析時(shí)。

北太天元代碼匯總

% 模擬數(shù)據(jù)生成  dt = 0.1; 
% 時(shí)間步長(zhǎng)  t = 0:dt:2; 
% 時(shí)間向量  t = t&#39;; x = zeros(size(t)); 
% 初始化狀態(tài)向量  x(1) = 0.6; 
% 初始條件  for i = 1:length(t)-1      x(i+1) = x(i) + dt*(-2*x(i)  +  x(i).^3 ); 
% 歐拉方法更新狀態(tài)  end  dxdt =  diff(x) ./diff(t); 
% 計(jì)算速度% 構(gòu)建函數(shù)庫(kù)  x = x(1:end-1); Theta = [x,  x.^2, x.^3]; 
% 候選函數(shù)庫(kù)% 稀疏回歸  % 注意：在實(shí)際應(yīng)用中，你可能需要對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理  
% LASSO回歸  lambda = 1e-9; maxIter = 100000;tol = 1e-9;
 [coef] = simpleLassoSolver(Theta, dxdt, lambda, maxIter, tol) norm( Theta*coef - dxdt)function [beta] = simpleLassoSolver(Theta, y, lambda, maxIter, tol)  % simpleLassoSolver: 
使用迭代軟閾值算法(ISTA)求解LASSO問(wèn)題  
% 輸入:  %   Theta - 設(shè)計(jì)矩陣 (n x p)  %   y - 響應(yīng)向量 (n x 1)  %   lambda - LASSO正則化參數(shù)  %   maxIter - 最大迭代次數(shù)  %   tol - 收斂容忍度  
% 輸出:  %   beta - 回歸系數(shù) (p x 1)    % 初始化  [n, p] = size(Theta);  beta = zeros(p, 1);  % 初始化解  x_old = zeros(p, 1); 
% 上一步的迭代解    for iter = 1:maxIter      
% 梯度下降步驟      grad = Theta&#39; * (Theta * beta - y) ;        
% 使用線性搜索確定最佳步長(zhǎng)      alpha = linesearch(Theta, y, beta, grad, lambda);      
% 更新解      x_new = beta - alpha * grad ;         
 % 應(yīng)用軟閾值操作      beta = sign(x_new) .* max(0, abs(x_new) - lambda);          
 % 檢查收斂性      fprintf([&#39;iter = &#39;, 
 num2str(iter), &#39; beta = &#39;, num2str(beta(1)) &#39;, 
 &#39; num2str(beta(2)),  &#39;, &#39;, num2str(beta(3)),  &#39;  
 norm = &#39; num2str( norm( Theta*beta - y) ), newline]);     
 if  norm( beta -x_old ) &lt; tol  &amp;&amp; norm ( Theta*beta - y) &lt; tol        break;      
 end      
 x_old = beta;  end    endfunction [alpha] = linesearch(Theta, y, beta, grad, lambda)  % 線性搜索尋找最佳步長(zhǎng)        % 初始化參數(shù)      alpha = 1; 
 % 初始步長(zhǎng)      c = 1e-4;  % Armijo條件中的常數(shù)      rho = 0.5; % 步長(zhǎng)縮減因子      max_ls_iter = 50; % 線性搜索的最大迭代次數(shù)            
 % 計(jì)算梯度下降的初始方向      d = -grad;            for ls_iter = 1:max_ls_iter          
 % 嘗試新的步長(zhǎng)          beta_new = beta + alpha * d;                   
  % 應(yīng)用軟閾值操作          beta_new = sign(beta_new) .* max(0, abs(beta_new) - lambda);                    
  % 計(jì)算目標(biāo)函數(shù)的值          f_new = 0.5 * norm(Theta * beta_new - y)^2 + lambda * norm(beta_new, 1);          f_old = 0.5 * norm(Theta * beta - y)^2 + lambda * norm(beta, 1);                    
  % 檢查Armijo條件是否滿足          if f_new &lt;= f_old + c * alpha * grad&#39; * d              break; 
  % 如果滿足，則退出線性搜索          
  end                   
% 如果不滿足，則縮減步長(zhǎng)并繼續(xù)搜索          
alpha = rho * alpha;      
end            
% 如果達(dá)到最大迭代次數(shù)仍未找到滿足條件的步長(zhǎng)，則可能需要調(diào)整參數(shù)或檢查算法實(shí)現(xiàn)      
if ls_iter == max_ls_iter          
warning(&#39;線性搜索達(dá)到最大迭代次數(shù)，未找到滿足條件的步長(zhǎng)。&#39;);      
end  
end

452 0 0 收藏回復(fù)

非線性動(dòng)力學(xué)的稀疏識(shí)別(SINDy)簡(jiǎn)介

回復(fù)

回復(fù)

推薦話題

熱門標(biāo)簽

社區(qū)牛人

產(chǎn)品介紹

生態(tài)合作

服務(wù)支持

關(guān)于我們