最优化

论文信息：Liu, L., Yang, B., Zhang, Y., Xu, Y., Peng, Z., & Wang, F. (2023). Solving Electromagnetic Inverse Problem Using Adaptive Gradient Descent Algorithm. IEEE Transactions on Geoscience and Remote Sensing, 61, 5902415. https://doi.org/10.1109/TGRS.2023.3239106 大地电磁（MT）等频率域电磁（EM）反演属于典型的非线性、非唯一地球物理反问题，工程上多采用线性化迭代求解。其中非线性共轭梯度（NLCG）算法因实现简单、内存需求低而广受欢迎，但它每一步都要做费力费时的线搜索来确定步长——这正是本文想要"动刀"的地方。本文作者之一张壹参与的研究把深度学习里早已成熟的自适应梯度下降（AGD）算法引入 EM 反演，用简单的代数累积代替线搜索，在保持结果可比的同时，把计算成本砍掉三分之一以上。下面围绕这篇 IEEE TGRS 论文，梳理其动机、方法要点与主要结论。 为什么 NLCG 反演"贵"在线搜索 在 NLCG 的每次迭代里，模型梯度决定了搜索方向，而线搜索负责确定步长。线搜索每次通常要多次求解相应的正演问题；即便采用高效的伴随方程法，梯度估计本身也大致要花一次正演两倍的计算量，而线搜索至少还要再额外付出一到两次正演。粗略估算下来，线搜索约占 NLCG 总成本的三分之一。对于需要反复迭代、每个样本都要正演的大型 EM 反演，这一块开销相当可观。 一个自然的想法是：能不能不做线搜索？深度学习社区在过去十年里给出了答案——自适应梯度下降（AGD）。它的核心是把历次迭代累积的梯度与模型更新用简单的代数运算来估计当前迭代的模型参数，从而彻底躲开耗时的线搜索。从 Adagrad 到 Adadelta、RMSProp，再到广泛使用的 Adam 及其变体（Adamax、NAdam、AMSGrad），这些算法在深度学习里已被验证得非常成熟。尽管在地震全波形反演中已有研究指出 AGD 相对 NLCG / L-BFGS 在时间与分辨率上的优势，但不同研究也存在分歧，且此前尚无将 AGD 用于 EM 数据反演的系统工作。这正是本研究要填补的空白。 方法：用累积梯度替代线搜索的 AGD 反演框架 AGD 家族共享一个朴素而漂亮的思想：不再费力去"试"一个标量步长，而是把每次迭代的梯度（及其平方）累积起来，用它们自适应地调节每个模型参数的学习率，从而把步长从标量变成向量（逐参数不同）。以最简洁的 Adagrad 为例，模型更新形如 ...

View on Github 线性共轭梯度算法库（C++ Library of Linear Conjugate Gradient，LIBLCG） 张壹（zhangyiss@icloud.com） 浙江大学地球科学学院·地球物理研究所 简介 liblcg是一个简单的C++线性共轭梯度算法库，包含了一般形式的共轭梯度与预优共轭梯度算法。可用于求解如下形式的线性方程组： Ax = B （可用共轭梯度求解）或者 PAx = PB（可用预优共轭梯度求解） 其中，A是一个N阶的对称矩阵、x为N*1的待求解的模型向量，B为N*1需拟合的目标向量。P则为预优算法中的预优矩阵，一般是一个N阶的对角阵。共轭梯度法广泛应用于无约束的线性最优化问题，拥有优良收敛与计算效率。 安装 GCC（MacOS） mkdir build cd build cmake .. make make install Windows 请自行拷贝代码新建项目并编译…嘿嘿😁！ 使用说明 自定义Ax计算函数 通常我们在使用共轭梯度法求解线性方程组Ax=B时A的维度可能会很大，直接储存A将消耗大量的内存空间，因此一般并不直接计算并储存A而是在需要的时候计算Ax的乘积。因此用户在使用liblcg时需要定义Ax的计算函数，同时提供初始解x与共轭梯度的B项（即拟合的对象）。如果使用预优方法还需要提供预优矩阵P项。Ax计算函数的形式必须满足算法库定义的一般形式： typedef void (*lcg_axfunc_ptr)(void* instance, const lcg_float* x, lcg_float* prod_Ax, const int n_size); 函数接受4个参数，分别为： void *instance 传入的实例对象； const lcg_float *input_array Ax计算中的x数组的指针； lcg_float *output_array Ax的乘积； const int n_size 矩阵的大小。 自定义进程监控函数 用户可以以下面的模版创建函数来显示共轭梯度迭代中的参数，并可以在适当的情况下停止迭代的进程。 typedef int (*lcg_progress_ptr)(void* instance, const lcg_float* m, const lcg_float converge, const lcg_para* param, const int n_size, const int k); 函数接收6个参数，分别为： ...

非线性共轭梯度法-DY法 算法可用于求解形如$min(f(x))$的非线性最优化问题，解的迭代形式为$x_{k+1}=x_k+a_kd_k$。其中$a_k$为迭代步长，$d_k$为迭代方向。（不同的$d_k$计算方法即为不同的非线性共轭梯度算法，它们对于$a_k$计算的要求也各有区别。） 要求目标函数$f(x)$连续可微且梯度$g(x)$已知。 DY法具有全局收敛性，但有一定几率陷入局部极小值。 DY法具有二次终止性，迭代终止于零梯度值位置或者$f(x) \leq threshold$处。 算法描述 step 1 -> 给出解的初值$f(x_0)$并计算其梯度$g(x_0)$，如果$\left\lVert g(x_0) \right\rVert = 0$则算法终止，说明$x_0$是$f(x)$的全局或局部极小解。否则设置$d_0 = -g_0$作为初始迭代方向（即解的最速下降方向）。 step 2 -> 执行满足如下所示的标准Wolfe条件的非精确线性搜索，计算得到$a_k$。 $f(x_k+a_k d_k)-f(x_k) \leq c_1 a_k {g_k}^T d_k$， ${g(x_k+a_k d_k)}^T d_K > c_2 {g_k}^T d_k$ 其中$0< c_1 < c_2 < 1$，一般取$c_1=0.1$，$c_2 \in (0.6,0.8)$。对于重磁非线形反演一般取$c_1=1e-4$，$c_2=0.9$。 step 3 -> 更新$x_{k+1}=x_k+a_kd_k $，若$\left\lVert g(x_{k+1}) \right\rVert = 0 $则算法终止，说明$x_{k+1}$是$f(x)$的全局或局部极小解。否则继续step 4。 step 4 -> 更新$d_{k+1}=-g_{k+1}+\beta_k d_k$，其中$\beta_k = \left\lVert g_{k+1} \right\rVert / {d_k}^T (g_{k+1}-g_k)$。$k:=k+1$返回step 2。 C++代码 func.h ...