非線性動力學的稀疏識別(SINDy)簡介

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

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

SINDy的基本原理

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

一個簡單的例子

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

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

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

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

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

% 模擬數(shù)據(jù)生成  dt = 0.1; 
% 時間步長  t = 0:dt:2; 
% 時間向量  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); % 計算速度

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

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

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

步驟3：稀疏回歸

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

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

$\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問題，其數(shù)學原理可以詳細闡述如下：

LASSO問題的定義

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

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

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

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

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

初始化

* 設(shè)定初始估計值 $\beta^0$ ，通常可以設(shè)為0向量或隨機向量。

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

迭代過程

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

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

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

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

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

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

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

終止條件

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

3. 算法收斂性

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

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

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

function [beta] = simpleLassoSolver(Theta, y, lambda, maxIter, tol)  % simpleLassoSolver: 
使用迭代軟閾值算法(ISTA)求解LASSO問題  % 輸入:  
%   Theta - 設(shè)計矩陣 (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) ;        
% 使用線性搜索確定最佳步長      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

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

結(jié)果分析

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

結(jié)論

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

北太天元代碼匯總

% 模擬數(shù)據(jù)生成  dt = 0.1; 
% 時間步長  t = 0:dt:2; 
% 時間向量  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); 
% 計算速度% 構(gòu)建函數(shù)庫  x = x(1:end-1); Theta = [x,  x.^2, x.^3]; 
% 候選函數(shù)庫% 稀疏回歸  % 注意：在實際應(yīng)用中，你可能需要對數(shù)據(jù)進行標準化處理  
% 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問題  
% 輸入:  %   Theta - 設(shè)計矩陣 (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) ;        
% 使用線性搜索確定最佳步長      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)  % 線性搜索尋找最佳步長        % 初始化參數(shù)      alpha = 1; 
 % 初始步長      c = 1e-4;  % Armijo條件中的常數(shù)      rho = 0.5; % 步長縮減因子      max_ls_iter = 50; % 線性搜索的最大迭代次數(shù)            
 % 計算梯度下降的初始方向      d = -grad;            for ls_iter = 1:max_ls_iter          
 % 嘗試新的步長          beta_new = beta + alpha * d;                   
  % 應(yīng)用軟閾值操作          beta_new = sign(beta_new) .* max(0, abs(beta_new) - lambda);                    
  % 計算目標函數(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                   
% 如果不滿足，則縮減步長并繼續(xù)搜索          
alpha = rho * alpha;      
end            
% 如果達到最大迭代次數(shù)仍未找到滿足條件的步長，則可能需要調(diào)整參數(shù)或檢查算法實現(xiàn)      
if ls_iter == max_ls_iter          
warning(&#39;線性搜索達到最大迭代次數(shù)，未找到滿足條件的步長。&#39;);      
end  
end

434 0 0 收藏回復(fù)

非線性動力學的稀疏識別(SINDy)簡介

回復(fù)

回復(fù)

推薦話題

熱門標簽

社區(qū)牛人

產(chǎn)品介紹

生態(tài)合作

服務(wù)支持

關(guān)于我們