在科学和工程领域，理解和预测动态系统的行为至关重要。然而，很多时候我们并不完全知道描述这些系统的精确数学模型。近年来，一种名为SINDy（Sparse Identification of Nonlinear Dynamics，非线性动力学的稀疏识别）的系统辨识方法受到了广泛关注，它能够从数据中自动发现动态系统的数学模型。
SINDy的基本原理
SINDy方法的核心思想是利用稀疏回归技术从测量数据中识别出动态系统的非线性微分方程。简单来说，它假设系统的动态可以由一个函数库中的少数几个函数的线性组合来描述。通过优化算法，SINDy能够选择出对系统动态影响最大的几个函数，并确定它们的系数，从而构建出一个简洁而精确的数学模型。
一个简单的例子
假设我们有一个简单的一维动态系统，其动态方程为：
为了模拟这个系统，我们生成一些模拟数据，并使用SINDy方法来恢复这个方程。
步骤1：生成模拟数据
首先，我们使用北太天元来生成模拟数据。假设初始条件为(x(0) = 1)，我们使用欧拉方法或更高级的数值方法来求解这个微分方程，并收集一段时间内的数据。
% 模拟数据生成 dt = 0.1; % 时间步长 t = 0:dt:2; % 时间向量 t = t'; x = zeros(size(t)); % 初始化状态向量 x(1) = 0.6; % 初始条件 for i = 1:length(t)-1 x(i+1) = x(i) + dt*(-2*x(i) + x(i).^3 ); % 欧拉方法更新状态 end dxdt = diff(x) ./diff(t); % 计算速度
步骤2：构建函数库
接下来，我们构建一个函数库，其中包含可能描述系统动态的候选函数。在这个例子中，我们可以选择多项式函数作为候选函数，例如(x), (x^2), (x^3), ...等。
% 构建函数库 x = x(1:end-1); Theta = [x, x.^2, x.^3]; % 候选函数库
步骤3：稀疏回归
最后，我们使用稀疏回归方法来确定哪些函数对系统动态有显著影响。在这个例子中，我们可以使用LASSO（Least Absolute Shrinkage and Selection Operator）回归，它是一种能够产生稀疏解的线性回归方法。
针对这个具体的例子，LASSO回归的优化问题可以表示为：
使用迭代软阈值算法(ISTA)求解LASSO问题，其数学原理可以详细阐述如下：
LASSO问题的定义
LASSO（Least Absolute Shrinkage and Selection Operator）是一种回归分析方法，它通过构造一个惩罚函数来得到一个较为精炼的模型。LASSO的基本思想是在回归系数的绝对值之和小于一个常数的约束条件下，使残差平方和最小化，从而能够产生某些严格等于0的回归系数，得到可以解释的模型。其数学形式通常表示为最小化以下目标函数：
其中，A 是设计矩阵，b 是观测向量，\beta 是待求的回归系数向量，λ 是正则化参数，控制惩罚项的影响。
2. 迭代软阈值算法(ISTA)的原理
迭代软阈值算法(ISTA)是用于求解LASSO问题的一种有效方法。它通过迭代的方式逐步逼近最优解。算法的基本步骤如下：
初始化
* 设定初始估计值 ，通常可以设为0向量或随机向量。
* 选择合适的步长 t_k 和正则化参数 λ 。
迭代过程
对于每一次迭代 ，执行以下步骤：
1. 梯度下降：首先，对当前估计值 进行梯度下降，以减小残差平方和的部分。这可以通过计算
并按一定步长
更新 β 来实现。
2. 软阈值处理：然后，对梯度下降后的结果进行软阈值处理，以实现L1正则化的效果。软阈值函数的形式为：
其中，\(\text{sgn}(x)\) 是符号函数，返回 \(x\) 的符号。这个软阈值函数会将绝对值小于 λ 的系数缩减到0，而对大于λ 的系数进行一定程度的缩减。
3. 更新估计值：将软阈值处理后的结果作为下一次迭代的起始点，即
。
终止条件
* 当 小于某个预设的阈值时，或者达到最大迭代次数时，算法停止。
3. 算法收敛性
在一定的条件下，ISTA算法可以保证收敛到LASSO问题的最优解。收敛速度取决于步长 的选择和问题的具体性质。通常，为了加速收敛，可以采用一些优化策略，如Nesterov加速或采用更复杂的步长选择策略。这里简单实现一个线性搜索确定步长方法
function [alpha] = linesearch(Theta, y, beta, grad, lambda) % 线性搜索寻找最佳步长 % 初始化参数 alpha = 1; % 初始步长 c = 1e-4; % Armijo条件中的常数 rho = 0.5; % 步长缩减因子 max_ls_iter = 50; % 线性搜索的最大迭代次数 % 计算梯度下降的初始方向 d = -grad; for ls_iter = 1:max_ls_iter % 尝试新的步长 beta_new = beta + alpha * d; % 应用软阈值操作 beta_new = sign(beta_new) .* max(0, abs(beta_new) - lambda); % 计算目标函数的值 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 <= f_old + c * alpha * grad' * d break; % 如果满足，则退出线性搜索 end % 如果不满足，则缩减步长并继续搜索 alpha = rho * alpha; end % 如果达到最大迭代次数仍未找到满足条件的步长，则可能需要调整参数或检查算法实现 if ls_iter == max_ls_iter warning('线性搜索达到最大迭代次数，未找到满足条件的步长。'); end end
北太天元中暂时没有内置的LASSO函数，可以用m函数实现一个简单的版本。
function [beta] = simpleLassoSolver(Theta, y, lambda, maxIter, tol) % simpleLassoSolver: 使用迭代软阈值算法(ISTA)求解LASSO问题 % 输入: % Theta - 设计矩阵 (n x p) % y - 响应向量 (n x 1) % lambda - LASSO正则化参数 % maxIter - 最大迭代次数 % tol - 收敛容忍度 % 输出: % beta - 回归系数 (p x 1) % 初始化 [n, p] = size(Theta); beta = zeros(p, 1); % 初始化解 x_old = zeros(p, 1); % 上一步的迭代解 for iter = 1:maxIter % 梯度下降步骤 grad = Theta' * (Theta * beta - y) ; % 使用线性搜索确定最佳步长 alpha = linesearch(Theta, y, beta, grad, lambda); % 更新解 x_new = beta - alpha * grad ; % 应用软阈值操作 beta = sign(x_new) .* max(0, abs(x_new) - lambda); % 检查收敛性 fprintf(['iter = ', num2str(iter), ' beta = ', num2str(beta(1)) ', ' num2str(beta(2)), ', ', num2str(beta(3)), ' norm = ' num2str( norm( Theta*beta - y) ), newline]); if norm( beta -x_old ) < tol && norm ( Theta*beta - y) < tol break; end x_old = beta; end end
% 稀疏回归 % 注意：在实际应用中，你可能需要对数据进行标准化处理 % LASSO回归 lambda = 1e-9; maxIter = 100000;tol = 1e-9; [coef] = simpleLassoSolver(Theta, dxdt, lambda, maxIter, tol) norm( Theta*coef - dxdt)
结果分析
通过检查回归系数coef，我们可以发现哪些函数对系统动态有显著影响。在这个例子中，我们希望恢复出原始的微分方程(\dot{x} = -2x + x^3)，因此我们应该期望在回归系数中看到-2和1对应于(x)和(x^3)的系数，而其他函数的系数应该接近于0。
结论
SINDy是一种强大的方法，能够从数据中自动发现动态系统的数学模型。通过结合稀疏回归技术和适当的函数库，SINDy可以准确地识别出描述系统动态的关键函数和它们的系数。这种方法在科学和工程领域具有广泛的应用前景，特别是在对复杂系统进行建模和分析时。
北太天元代码汇总
% 模拟数据生成 dt = 0.1; % 时间步长 t = 0:dt:2; % 时间向量 t = t'; x = zeros(size(t)); % 初始化状态向量 x(1) = 0.6; % 初始条件 for i = 1:length(t)-1 x(i+1) = x(i) + dt*(-2*x(i) + x(i).^3 ); % 欧拉方法更新状态 end dxdt = diff(x) ./diff(t); % 计算速度% 构建函数库 x = x(1:end-1); Theta = [x, x.^2, x.^3]; % 候选函数库% 稀疏回归 % 注意：在实际应用中，你可能需要对数据进行标准化处理 % 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 - 设计矩阵 (n x p) % y - 响应向量 (n x 1) % lambda - LASSO正则化参数 % maxIter - 最大迭代次数 % tol - 收敛容忍度 % 输出: % beta - 回归系数 (p x 1) % 初始化 [n, p] = size(Theta); beta = zeros(p, 1); % 初始化解 x_old = zeros(p, 1); % 上一步的迭代解 for iter = 1:maxIter % 梯度下降步骤 grad = Theta' * (Theta * beta - y) ; % 使用线性搜索确定最佳步长 alpha = linesearch(Theta, y, beta, grad, lambda); % 更新解 x_new = beta - alpha * grad ; % 应用软阈值操作 beta = sign(x_new) .* max(0, abs(x_new) - lambda); % 检查收敛性 fprintf(['iter = ', num2str(iter), ' beta = ', num2str(beta(1)) ', ' num2str(beta(2)), ', ', num2str(beta(3)), ' norm = ' num2str( norm( Theta*beta - y) ), newline]); if norm( beta -x_old ) < tol && norm ( Theta*beta - y) < tol break; end x_old = beta; end endfunction [alpha] = linesearch(Theta, y, beta, grad, lambda) % 线性搜索寻找最佳步长 % 初始化参数 alpha = 1; % 初始步长 c = 1e-4; % Armijo条件中的常数 rho = 0.5; % 步长缩减因子 max_ls_iter = 50; % 线性搜索的最大迭代次数 % 计算梯度下降的初始方向 d = -grad; for ls_iter = 1:max_ls_iter % 尝试新的步长 beta_new = beta + alpha * d; % 应用软阈值操作 beta_new = sign(beta_new) .* max(0, abs(beta_new) - lambda); % 计算目标函数的值 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 <= f_old + c * alpha * grad' * d break; % 如果满足，则退出线性搜索 end % 如果不满足，则缩减步长并继续搜索 alpha = rho * alpha; end % 如果达到最大迭代次数仍未找到满足条件的步长，则可能需要调整参数或检查算法实现 if ls_iter == max_ls_iter warning('线性搜索达到最大迭代次数，未找到满足条件的步长。'); end end
公众号