目录
1. 数值分析
1.1 误差
- 理解绝对误差、相对误差、有效数字
1.2.线性方差组的直接求解
-
线性方程组的直接求解、求解矛盾（超定）方程组
-
对于方程（如下），如果无解（矛盾方程组），也就是 x = A − 1 x=A^{-1} x=A−1不成立，也就是 A A A不可逆
A x = b Ax=b Ax=b
我们一般会选择求解满足最小化方程的解 x x x，如下
F ( x ) = a r g m i n ( ∣ ∣ A x − b ∣ ∣ 2 ) F(x)=argmin(||Ax-b||^2) F(x)=argmin(∣∣Ax−b∣∣2)
通过求导等操作，得到 x x x的解如下
x = ( A T A ) − 1 A T b x=(A^TA)^{-1}A^Tb x=(ATA)−1ATb
在求解 ( A T A ) − 1 (A^TA)^{-1} (ATA)−1一般使用cholesky,也就是平方根法将 A T A A^TA ATA 分解为下三角和上三角的相乘 ( A T A = L L T A^TA=LL^T ATA=LLT ), 再展开求三角的逆。分解的前提是假设 A T A A^TA ATA是正定的。该方法相比LU分解更加高效，数值更加稳定。
1.3 线性方程组的迭代求解
- 对线性方程
组求解
A X = b AX=b AX=b - 雅可比迭代方法，Gaussian-Siderl方法（重点）
- Gaussian-Siderl求解过程
- 将系数矩阵构造满足 A = D − L − U A=D-L-U A=D−L−U。其中：对角 D D D、负下三角 L L L、负上三角 U U U
- 得到
迭代矩阵M = ( D − L ) − 1 U M=(D-L)^{-1}U M=(D−L)−1U、 F = ( D − L ) − 1 b F=(D-L)^{-1}b F=(D−L)−1b - 若没有给出
x
k
x_k
xk,则一般为
x
k
=
[
0
,
0...
]
T
x_k=[0,0...]^T
xk=[0,0...]T。然后进入迭代, 直到
最大迭代次数或者最大范数才停止。
x n + 1 = M x n + F x_{n+1} = M x_n+F xn+1=Mxn+F
1.5 插值法
- 拉格朗日多项式插值、牛顿插值、埃尔米特插值（了解）
1.7 数值微分和数值积分
- 数值微分：一阶、二阶；数值积分；
1.8 非线性方程及非线性方程组的求解（包含LM算法）
-
非线性方程，比如求解下面的非线性方程，实际上这也是一个残差函数，残差值r ( x ) = f ( x ) = f ( x ) − 0 r(x)=f(x)=f(x)-0 r(x)=f(x)=f(x)−0
f ( x ) = x 3 + 10 x − 20 = 0 f(x)=x^3+10x-20=0 f(x)=x3+10x−20=0 -
非线性方程：对分区间法、简单迭代法、牛顿法(重点)
-
牛顿法求解
非线性方程，求解过程 = 线性化 + 迭代的过程- 开始前会给出 n n n时刻满足方程的粗略解，即 x n = x 0 x_n=x_0 xn=x0。这里的 x n x_n xn是一个迭代变量， x 0 x_0 x0是一个超参数
- 根据非线性方程 f ( x ) = 0 f(x)=0 f(x)=0，对其求导，得到 x n x_n xn的导数 f ′ ( x n ) f'(x_n) f′(xn)值，意味着会随着迭代而更新
- (
这里只是为了方便理解)通过 f ( x ) = 0 f(x)=0 f(x)=0在 x n x_n xn的一阶线性展开，然后建立更精确的 x n + 1 x_{n+1} xn+1满足以下式
f ( x n + 1 ) = f ( x n ) + f ′ ( x n ) ( x n + 1 − x n ) = 0 f(x_{n+1}) = f(x_n)+f'(x_n)(x_{n+1}-x_n)=0 f(xn+1)=f(xn)+f′(xn)(xn+1−xn)=0 - 得到非线性方程求解迭代式, 然后得到
x
n
+
1
x_{n+1}
xn+1,并将其作为新的
x
n
x_n
xn，进行下一次迭代，直到
变化量小于某个阈值
x n + 1 = x n − f ( x n ) f ′ ( x n ) x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} xn+1=xn−f′(xn)f(xn)
-
（拓展内容）
非线性方程组的求解: 需要引入雅可比矩阵- 设有非线性方程组
F ( x 1 , x 2 , . . . , x n ) { f 1 ( x 1 , x 2 , . . . , x n ) = 0 f 2 ( x 1 , x 2 , . . . , x n ) = 0 . . . f n ( x 1 , x 2 , . . . , x n ) = 0 F(x_1, x_2,...,x_n)\left\{\begin{array}{cc} f_1(x_1, x_2,... ,x_n)=0 \\ f_2(x_1, x_2,... ,x_n)=0 \\ ... \\ f_n(x_1, x_2,... ,x_n)=0 \\ \end{array}\right. F(x1,x2,...,xn)⎩ ⎨ ⎧f1(x1,x2,...,xn)=0f2(x1,x2,...,xn)=0...fn(x1,x2,...,xn)=0
- 则它在
第k次迭代时的雅可比矩阵如下（意味着每一次迭代都需要更新），注意左边的x是一个向量
F ′ ( x k ) = [ ∂ f 1 ∂ x 1 ∂ f 1 ∂ x 2 . . . ∂ f 1 ∂ x n ∂ f 2 ∂ x 1 ∂ f 2 ∂ x 2 . . . ∂ f 2 ∂ x n . . . ∂ f n ∂ x 1 ∂ f n ∂ x 2 . . . ∂ f n ∂ x n ] F'(x^k)=\begin{bmatrix} \frac {\partial f_1}{\partial x_1} & \frac {\partial f_1}{\partial x_2} & ... & \frac {\partial f_1}{\partial x_n} \\ \frac {\partial f_2}{\partial x_1} & \frac {\partial f_2}{\partial x_2} & ... & \frac {\partial f_2}{\partial x_n} \\ &&... \\ \frac {\partial f_n}{\partial x_1} & \frac {\partial f_n}{\partial x_2} & ... & \frac {\partial f_n}{\partial x_n} \\ \end{bmatrix} F′(xk)= ∂x1∂f1∂x1∂f2∂x1∂fn∂x2∂f1∂x2∂f2∂x2∂fn............∂xn∂f1∂xn∂f2∂xn∂fn - 根据
非线性方程求解过程进行拓展得到非线性方程组的求解，其中的 x k x^{k} xk是一个向量。迭代的标准一般为函数、x中的最大增量等
x k + 1 = x k − ( F ′ ( x k ) ) − 1 F ( x k ) x^{k+1} = x^{k} - (F'(x^k))^{-1}F(x^k) xk+1=xk−(F′(xk))−1F(xk)
-
（拓展内容）
优化问题的求解:先进行梯度下降得到较精确的解，然后进行非线性方程组求解得到更好的解。本质是一个求极值的问题- 如果
已经建立非线性方程组 F ( x 1 , x 2 , . . . , x n ) F(x_1, x_2,...,x_n) F(x1,x2,...,xn) ，此时构造了模函数如下, 显然就是求解零点。其中 m m m是 f f f的个数
ϕ ( x 1 , x 2 , . . . , x n ) = ∑ i = 0 m [ f i ( x 1 , x 2 , . . . , x n ) ] 2 \phi(x_1, x_2,...,x_n)=\sum_{i=0}^m [ f_i(x_1, x_2,...,x_n) ]^2 ϕ(x1,x2,...,xn)=i=0∑m[fi(x1,x2,...,xn)]2 - 给出满足方程的粗略解，作为迭代的初始值，因此赋值给 x k = ( x 1 , x 2 , . . . , x n ) x^k=(x_1,x_2,...,x_n) xk=(x1,x2,...,xn)
- 计算当前的
x
k
x^k
xk 在
ϕ
(
x
1
,
x
2
,
.
.
.
,
x
n
)
\phi(x_1, x_2,...,x_n)
ϕ(x1,x2,...,xn) 的梯度值，因此会随着迭代过程而更新
G ( x k ) = ( ∂ ϕ ∂ x 1 ， ∂ ϕ ∂ x 2 ， . . . ， ∂ ϕ ∂ x n ) G(x^k)=( \frac {\partial \phi}{\partial x_1}， \frac {\partial \phi}{\partial x_2} ，...， \frac {\partial \phi}{\partial x_n}) G(xk)=(∂x1∂ϕ，∂x2∂ϕ，...，∂xn∂ϕ) - 得到迭代值
x
k
+
1
x^{k+1}
xk+1，其中
λ
\lambda
λ 是一个超参数（按理说最优值
x
k
x^k
xk 满足任意的
λ
\lambda
λ 都使得
∣
x
k
−
x
k
+
1
∣
|x^k - x^{k+1}|
∣xk−xk+1∣ 足够小）
x k + 1 = x k − λ G ( x k ) x^{k+1} = x^{k} - \lambda G(x^k) xk+1=xk−λG(xk) - 等待上面的迭代差不多了(称为最速下降法，收敛慢，因为它是线性收敛)。需要进一步
利用求解非线性方程组进行优化，此时只需要将 x k + 1 x^{k+1} xk+1 作为非线性方程组的初始解就行。
- 如果
-
（拓展内容） 关于SLAM方向的同学可以参考这篇文章：非线性优化方法小结-(最小二乘，梯度下降，高斯牛顿, LM)
-
损失函数：是所有残差函数平方的和，是一个单值。残差函数，是一个 m m m 维列向量，一般是值当前函数值与目标函数值的差值，默认差值期望为0.
F ( x ) = ∑ i = 1 m ( f 2 ( x 1 . . . x n ) ) = f T ( x ) f ( x ) ，其中 f ( x ) = [ f 1 ( x 1 . . . x n ) f 2 ( x 1 . . . x n ) . . . f m ( x 1 . . . x n ) ] ，且 J = [ δ f 1 ( x 1 . . . x n ) δ x δ f 2 ( x 1 . . . x n ) δ x . . . δ f m ( x 1 . . . x n ) δ x ] m × n F(x)=\sum_{i=1}^{m}( f^2(x_1...x_n))=f^T(x)f(x)，其中f(x)= \begin{bmatrix} f_1(x_1...x_n) \\ f_2(x_1...x_n) \\ ... \\ f_m(x_1...x_n) \end{bmatrix}， 且\ J= \begin{bmatrix} \frac{\delta f_1(x_1...x_n)}{\delta x}\\ \frac{\delta f_2(x_1...x_n)}{\delta x} \\ ... \\ \frac{\delta f_m(x_1...x_n)}{\delta x} \end{bmatrix}_{m\times n} F(x)=i=1∑m(f2(x1...xn))=fT(x)f(x)，其中f(x)= f1(x1...xn)f2(x1...xn)...fm(x1...xn) ，且 J= δxδf1(x1...xn)δxδf2(x1...xn)...δxδfm(x1...xn) m×n -
非线性的优化问题一般归结于最小二乘问题，其中 x = [ x 1 ， x 2 ， . . . ， x n ] x=[x_1，x_2，...，x_n] x=[x1，x2，...，xn] 是待优化的值
a r g m i n ⏟ x 1 2 ∑ ( ∣ ∣ f ( x 1 ， x 2 ， . . . ， x n ) ∣ ∣ 2 ) \underbrace{argmin}_{x} \frac{1}{2} \sum(||f(x_1，x_2，...，x_n)||^2) x argmin21∑(∣∣f(x1，x2，...，xn)∣∣2) -
f ( x ) f(x) f(x)是子残差函数 f i ( x ) f_i(x) fi(x)构成的向量，并且是一个非线性函数，它的一阶泰勒展开近似有 , x x x
在每一轮泰勒展开中一般是固定的
f ( x + △ x ) ≈ ι ( x + △ x ) = f ( x ) + J ⋅ △ x f(x+\triangle x) \approx \iota(x + \triangle x)=f(x) + J \cdot \triangle x \\ f(x+△x)≈ι(x+△x)=f(x)+J⋅△x
需要注意这里的 J J J 是残差函数 f f f 的雅可比矩阵，带入损失函数中可得
F ( x + △ x ) ≈ L ( x + △ x ) = 1 2 ι T ι = 1 2 f T f + △ x T ⋅ J T f + 1 2 △ x T ⋅ J T J ⋅ △ x = F ( x ) + △ x T ⋅ J T f + 1 2 △ x ⋅ J T J ⋅ △ x \begin{align} F(x+\triangle x) & \approx L(x + \triangle x) \\ &=\frac{1}{2}\iota^T \iota \\ &= \frac{1}{2}f^Tf + \triangle x^T \cdot J^Tf + \frac{1}{2}\triangle x^T \cdot J^TJ\cdot \triangle x \\ &=F(x) + \triangle x^T \cdot J^Tf + \frac{1}{2}\triangle x \cdot J^TJ \cdot \triangle x \end{align} F(x+△x)≈L(x+△x)=21ιTι=21fTf+△xT⋅JTf+21△xT⋅JTJ⋅△x=F(x)+△xT⋅JTf+21△x⋅JTJ⋅△x
此时可以得到 F ( x ) F(x) F(x)的雅可比矩阵 F ′ ( x 1 . . . x n ) = ( J T f ) T F'(x_1...x_n)=(J^Tf)^T F′(x1...xn)=(JTf)T，二阶导 F ′ ‘ ( x 1 . . . x n ) = J T J F'‘(x_1...x_n)=J^TJ F′‘(x1...xn)=JTJ -
使得展开后的一阶项、或者二阶项为0，可以求解出 △ x \triangle x △x的增量。 -
保留不同的高阶项进行优化， x k x_k xk 进行反梯度方向前进，衍生出了不同的算法
最速下降法：仅保留一阶，省去高阶后，则 F ( x ) F(x) F(x)在 x k x_k xk处的梯度为 J k J_k Jk，此时期望的 △ x = − J k T \triangle x= -J_k^T △x=−JkT，即 x k + 1 − x k = − J k T x_{k+1} - x_k= -J_k^T xk+1−xk=−JkT
牛顿法：仅保留一二阶，省去高阶后，为了得到极值，则 F ( x ) F(x) F(x)对 △ x \triangle x △x处的梯度应为0，此时新的关系满足 J k + H k △ x = 0 J_k+H_k \triangle x=0 Jk+Hk△x=0，即 x k + 1 − x k = − H k − 1 J k T x_{k+1} - x_k= -H_k^{-1} J_k^T xk+1−xk=−Hk−1JkT
高斯牛顿法：它构造了一个新的优化函数，就是上面的式子，为了取得新函数的极值，因此得到了新的关系满足 J k T J k △ x = − J k T f ( x k ) J_k^T J_k \triangle x=-J_k^T f(x_k) JkTJk△x=−JkTf(xk)，即 x k + 1 − x k = − ( J k T J k ) − 1 J k T f ( x k ) x_{k+1} - x_k=-(J_k^T J_k )^{-1}J_k^T f(x_k) xk+1−xk=−(JkTJk)−1JkTf(xk)。
Levenberg-Marquadt：在高斯牛顿法中，给 △ x \triangle x △x 添加一个信赖区域(Trust Region)，不能让它变化太大使得近似不准确。LM 算法如下，所以很多地方近似为 H ⋅ δ x = b H\cdot\delta x = b H⋅δx=b
( J T J + u I ) ⋅ △ x = − J T f ，且 u > 0 (J^TJ+uI)\cdot \triangle x=-J^Tf，且u>0 (JTJ+uI)⋅△x=−JTf，且u>0- 下面是
u
u
u的更新过程
- u u u 的作用： u > 0 u>0 u>0能够保证 J T J + u J J^TJ+uJ JTJ+uJ 是一个正定矩阵，迭代朝着下降方向进行
-
u
u
u的初始化：一个简单的
u
0
u_0
u0 的初始化就是
u 0 = τ ⋅ m a x ( [ J T J ] i i ) ，其中 τ ～ [ 10 − 8 ， 1 ] u_0=\tau \cdot max([J^TJ]_{ii})，其中 \tau ～[10^{-8}，1] u0=τ⋅max([JTJ]ii)，其中τ～[10−8，1]
然后通过计算 ρ \rho ρ，来控制 u u u的变化，如下
ρ = F ( x ) − F ( x + △ x l m ) L ( 0 ) − L ( △ x l m ) = F ( x ) − F ( x + △ x l m ) 1 2 △ x l m T ( u △ x l m − J T f ) \begin{align} \rho &=\frac{F(x)-F(x+\triangle x_{lm})}{L(0)-L(\triangle x_{lm})} \\ &=\frac{F(x)-F(x+\triangle x_{lm})}{ \frac{1}{2}\triangle x^T_{lm}(u\triangle x_{lm}-J^Tf ) } \end{align} ρ=L(0)−L(△xlm)F(x)−F(x+△xlm)=21△xlmT(u△xlm−JTf)F(x)−F(x+△xlm)
最后， u u u的更新如下
u = { u ∗ m a x ( 1 3 ， 1 − ( 2 ρ − 1 ) 3 ) ; v = 2 ， i f ρ > 0 u ∗ v ； v = 2 ∗ v ， e l s e u=\begin{cases} u*max(\frac{1}{3}，1-(2\rho-1)^3);v=2，& if \ \rho > 0 \\ u*v；v=2*v，&else \end{cases} u={u∗max(31，1−(2ρ−1)3);v=2，u∗v；v=2∗v，if ρ>0else
- 下面是
u
u
u的更新过程
-
下面是相关算法的总结
算法类型 算法介绍 迭代过程 最速下降法 负梯度方向，收敛速度慢 x k + 1 − x k = − J k T x_{k+1} - x_k= -J_k^T xk+1−xk=−JkT Newton 法 保留泰勒级数一阶和二阶项，二次收敛速度，但每步都计算Hessian矩阵，复杂 x k + 1 − x k = − H k − 1 J k T x_{k+1} - x_k= -H_k^{-1} J_k^T xk+1−xk=−Hk−1JkT 高斯牛顿法 目标函数的 J k T J k J_k^T J_k JkTJk 近似H矩阵，提高算法效率，但H矩阵不满秩则无法迭代 x k + 1 − x k = − ( J k T J k ) − 1 J k T f ( x k ) x_{k+1} - x_k=-(J_k^T J_k )^{-1}J_k^T f(x_k) xk+1−xk=−(JkTJk)−1JkTf(xk) Levenberg-Marquadt 信赖域算法，解决H矩阵不满秩或非正定 △ x = − ( J k T J k + u I ) − 1 J k T f ( x k ) \triangle x= -(J_k^T J_k + u I )^{-1} J_k^T f(x_k) △x=−(JkTJk+uI)−1JkTf(xk) 注意： f ( x k ) f(x_k) f(xk)是一个列向量，代表每个子方程的残差。
-
1.8 核函数的引入
-
为什么引入核函数
核函数（Kernel Function）在优化问题中，尤其是在存在噪声或异常测量的情况下，起到了增强优化稳健性的作用。通过对误差或残差进行修正，核函数能够减小异常值对优化过程的影响，使得优化结果更加可靠。具体来说，核函数有以下几种主要作用：-
增强优化的稳健性
核函数的一个关键作用是增强优化过程的稳健性。没有核函数时，较大的残差（如由异常值或错误测量引起）可能会主导优化过程，导致优化结果偏离真实值。核函数通过减少大残差的权重，避免它们对优化过程产生过大影响，从而提高优化的可靠性。 -
减少异常值的影响
在实际应用中，观测数据往往受到噪声的干扰，可能会出现极端值。没有核函数时，这些异常值可能会导致优化结果的失真。通过引入核函数，尤其是当残差较大时，核函数能够有效减小异常数据对优化结果的影响，使得优化过程更加可靠。 -
控制残差对优化过程的影响
核函数还能够控制残差对优化过程的影响。当残差较小（误差较小）时，核函数的影响较小，优化过程类似于标准的最小二乘方法；但当残差较大时，核函数会对残差进行缩放，使它们对目标函数的贡献减小，从而减少大误差对优化过程的负面影响。 -
改善数值稳定性
在某些情况下，优化过程中可能会出现数值不稳定的问题，尤其是在处理大尺度数据或噪声较大的情况下。使用核函数可以有效提高数值稳定性，避免由于大残差或异常值引起的梯度更新过大或计算不准确的问题。
-
-
核函数的常见类型
- 若定义残差函数如下：
r ( x ) = f ( x ) − 期望值 r(x)=f(x) - 期望值 r(x)=f(x)−期望值 - 则引入核函数后的残差为,
s
s
s
一般代表计算的残差值平方
ρ ( s ( x ) ) ρ( s(x) ) ρ(s(x))
- 若定义残差函数如下：
-
常见的核函数
- Huber 核函数，
s
s
s
一般代表计算的残差值平方, δ \delta δ是超参数，Huber 核函数通常用于最小化二次误差
ρ ( s ) = { 1 2 s 2 i f ∣ s ∣ ≤ δ δ ⋅ ( ∣ s ∣ − 1 2 δ ) i f ∣ s ∣ ≥ δ ρ( s ) = \left\{ \begin{matrix} \frac{1}{2} s^2 & if |s| \leq\delta\\ \delta \cdot(|s| - \frac{1}{2}\delta) & if |s| \geq\delta\\ \end{matrix} \right. ρ(s)={21s2δ⋅(∣s∣−21δ)if∣s∣≤δif∣s∣≥δ - Cauchy 核函数，
s
s
s
一般代表计算的残差值平方， σ \sigma σ是超参数，Cauchy 核函数是一种常见的稳健核函数
ρ ( s ) = l o g ( 1 + s 2 σ 2 ) ρ(s ) = log(1+\frac{s^2}{\sigma ^2}) ρ(s)=log(1+σ2s2)
- Huber 核函数，
s
s
s
-
引入核函数后的雅克比矩阵计算-
一般而言，我们将其展开到二阶，然后左右都乘一个 1 2 \frac{1}{2} 21，这里把 s k s_k sk当作一个整体看
1 2 ρ ( s k ) ≈ 1 2 ( 常数 值 ρ ( s k = 具体 ) + ∂ ρ ∂ s k ⋅ △ s k + 1 2 ∂ 2 ρ ∂ s k 2 ⋅ ( △ s k ) 2 ) \frac{1}{2}ρ( s_k) \approx \frac{1}{2}( 常数值_{ρ( s_k=具体)}+ \frac {\partial ρ}{ \partial s_k} \cdot \triangle s_k + \frac{1}{2}\frac {\partial ^2 ρ}{ \partial s_k^2} \cdot (\triangle s_k)^2 ) 21ρ(sk)≈21(常数值ρ(sk=具体)+∂sk∂ρ⋅△sk+21∂sk2∂2ρ⋅(△sk)2) -
下面推导具有核函数的优化算法
- 核函数
ρ
(
⋅
)
ρ(\cdot)
ρ(⋅)作用于残差结果上
r
(
x
)
r(x)
r(x)上，因此最小二乘的形式如下
m i n ⏟ x 1 2 ∑ k ρ ( ∣ ∣ r k ( x ) ∣ ∣ 2 ) \underbrace{min}_{x} \frac{1}{2} \sum_k ρ(||r_k(x)||^2) x min21k∑ρ(∣∣rk(x)∣∣2) - 首先定义第k个残差的二次项的增量表达 s k ( x ) = ∣ ∣ r k ( x ) ∣ ∣ 2 s_k(x) =||r_k(x)||^2 sk(x)=∣∣rk(x)∣∣2,由此可以得到这个二次项的增量与x的关系如下， J J J是当前第k个残差向量 r k r_k rk 关于参数 x x x 的雅克比矩阵
△ s k = ∣ ∣ r k ( x + △ x ) ∣ ∣ 2 − ∣ ∣ r k ( x ) ∣ ∣ 2 ≈ ∣ ∣ r k + J k △ x ∣ ∣ 2 − ∣ ∣ r ( x ) ∣ ∣ 2 = 2 r k T ( J k △ x ) + ( J k △ x ) T J k △ x \begin{matrix} \triangle s_k&=& ||r_k(x+\triangle x)||^2 - ||r_k(x)||^2 \\ \\ &\approx&||r_k+J_k\triangle x||^2-||r(x)||^2 \\ \\ &=& 2r_k^T(J_k\triangle x) + (J_k\triangle x)^TJ_k\triangle x \end{matrix} △sk=≈=∣∣rk(x+△x)∣∣2−∣∣rk(x)∣∣2∣∣rk+Jk△x∣∣2−∣∣r(x)∣∣22rkT(Jk△x)+(Jk△x)TJk△x
- 将
△
s
k
\triangle s_k
△sk带入函数
1
2
ρ
(
s
k
)
\frac{1}{2}ρ(s_k)
21ρ(sk)中，然后展开里面的平方项可以得到如下近似
1 2 ρ ( ∣ ∣ r k ( x ) ∣ ∣ 2 ) ≈ ∂ ρ ∂ r ⋅ r k T ⋅ ( J k △ x ) + 1 2 ( J k △ x ) T ( ∂ ρ ∂ r I + 2 ∂ 2 ρ ∂ r 2 r k r k T ) ( J k △ x ) + 常数项 \frac{1}{2} ρ(||r_k(x)||^2) \approx \frac {\partial ρ}{ \partial r} \cdot r_k^T \cdot (J_k\triangle x)+\frac{1}{2}(J_k\triangle x)^T(\frac {\partial ρ}{ \partial r}I+2\frac {\partial^2 ρ}{ \partial r^2}r_kr_k^T)(J_k\triangle x)+常数项 21ρ(∣∣rk(x)∣∣2)≈∂r∂ρ⋅rkT⋅(Jk△x)+21(Jk△x)T(∂r∂ρI+2∂r2∂2ρrkrkT)(Jk△x)+常数项 - 因为
最小二乘的展开形式是上面式子的求和形式，因此我们先对上面式子求和。
然后对上式子左右两边同时对 △ x \triangle x △x求导，最终可以得到该最小二乘的
∑ k [ J k T ( ∂ ρ ∂ r I + 2 ∂ 2 ρ ∂ r 2 r k r k T ) ( J k △ x ) ] = − ∑ k ( ∂ ρ ∂ r J k T r k ) \sum_k[J_k^T(\frac {\partial ρ}{ \partial r}I+2\frac {\partial^2 ρ}{ \partial r^2}r_kr_k^T)(J_k\triangle x)]=-\sum_k(\frac {\partial ρ}{ \partial r} J_k^T r_k) k∑[JkT(∂r∂ρI+2∂r2∂2ρrkrkT)(Jk△x)]=−k∑(∂r∂ρJkTrk) - 更一般的定义
W
=
(
∂
ρ
∂
r
I
+
2
∂
2
ρ
∂
r
2
r
k
r
k
T
)
W=(\frac {\partial ρ}{ \partial r}I+2\frac {\partial^2 ρ}{ \partial r^2}r_kr_k^T)
W=(∂r∂ρI+2∂r2∂2ρrkrkT)
是一个权重因子，得到如下的简洁表达式
∑ k J k T W J k △ x = − ∑ k ( ∂ ρ ∂ r J k T r k ) \sum_kJ_k^TWJ_k\triangle x=-\sum_k(\frac {\partial ρ}{ \partial r} J_k^T r_k) k∑JkTWJk△x=−k∑(∂r∂ρJkTrk) - 最终引入核函数后的LM的增量迭代表达式为：
△ x = − ( J k T W J k + u I ) − 1 ⋅ ( ∂ ρ ∂ r J k T r k ) \triangle x=-(J_k^TWJ_k + uI)^{-1}\cdot (\frac {\partial ρ}{ \partial r}J_k^T r_k) △x=−(JkTWJk+uI)−1⋅(∂r∂ρJkTrk)
- 核函数
ρ
(
⋅
)
ρ(\cdot)
ρ(⋅)作用于残差结果上
r
(
x
)
r(x)
r(x)上，因此最小二乘的形式如下
-
1.9 常微分方程的数值解法
- 欧拉法、改进欧拉、龙格库塔
2. 常微分方程数值解法-Python实现
2.1 一阶微分方程
-
一阶微分方程简介，下面便是一般求解的微分方程
-
四阶龙格库塔方法
一阶微分方程解法代码class Runge_Kutta: def __init__(self) -> None: pass # 原函数的导函数 def f_xy(self, x, y): value = x - y return value # 由当前点和步进值，计算下一个点 def run(self, start=[0, 0], step=0.01 ): # cal x0 y0 x0, y0 = start K1 = self.f_xy(x0, y0) K2 = self.f_xy(x0 + step/2, y0 + (step/2)*K1 ) K3 = self.f_xy(x0 + step/2, y0 + (step/2)*K2 ) K4 = self.f_xy(x0 + step, y0 + step*K3) Next_y = y0 + (step / 6) * (K1 + 2*K2 + 2*K3 + K4) Next_x = x0 + step return round(Next_x, 6), round(Next_y, 6)用法：
# 实例化 RK = Runge_Kutta() # 传入起始值， 步进值， 获得下一个 xy 点 x, y = RK.run(start=(0, 0), step=0.2) # 可以使用此语句进行迭代 x, y = RK.run(start=(x, y), step=0.2) # 0 --- : 0.6 0.148817 # 1 --- : 0.8 0.249335 # 2 --- : 1.0 0.367886 # 3 --- : 1.2 0.501201 # 4 --- : 1.4 0.646603 # 5 --- : 1.6 0.801902 # 6 --- : 1.8 0.965304 # 7 --- : 2.0 1.13534 # 8 --- : 2.2 1.310807 # 9 --- : 2.4 1.490721
2.2. 广义微分方程
- 相关代码
class Equ: def __init__(self) -> None: pass # 原函数的导函数 def f_xy(self, x, y): value = x - y return value # 欧拉法 def euler(self, start=[0, 0], step=0.01): # cal x0 y0 x0, y0 = start Next_y = y0 + step * self.f_xy(x0, y0) Next_x = x0 + step return round(Next_x, 6), round(Next_y, 6) # 改进欧拉法 def euler_improve(self, start=[0, 0], step=0.01): # cal x0 y0 x0, y0 = start yp = y0 + step * self.f_xy(x0, y0) yq = y0 + step * self.f_xy(x0 + step, yq) Next_y = 0.5 * (yq + yp) Next_x = x0 + step return round(Next_x, 6), round(Next_y, 6) # 四阶龙哥库达算法 def RK(self, start=[0, 0], step=0.01 ): # cal x0 y0 x0, y0 = start K1 = self.f_xy(x0, y0) K2 = self.f_xy(x0 + step/2, y0 + (step/2)*K1 ) K3 = self.f_xy(x0 + step/2, y0 + (step/2)*K2 ) K4 = self.f_xy(x0 + step, y0 + step*K3) Next_y = y0 + (step / 6) * (K1 + 2*K2 + 2*K3 + K4) Next_x = x0 + step return round(Next_x, 6), round(Next_y, 6)
2.3. 高阶微分方程
- 简介
引进新的变量将高阶微分方程归结为一阶微分方程组的初值问题来求解。
一般而言，高阶微分方程可以改写为状态方程。
例如：
y ( x ) ′ ′ + 3 y ( x ) ′ + 2 y ( x ) = 4 x y(x)''+3y(x)'+2y(x)=4x y(x)′′+3y(x)′+2y(x)=4x
令
{ x 1 = y x 2 = x 1 ′ \begin{cases} x_1=y\\ x_2=x_1' \end{cases} {x1=yx2=x1′
得到新的微分方程
{ x 1 ′ = x 2 x 2 ′ = − 2 x 1 − 3 x 2 + 4 x \begin{cases} x_1'=x_2 \\ x_2'=-2x_1-3x_2+4x \end{cases} {x1′=x2x2′=−2x1−3x2+4x
故求解 y y y 的过程就是求解 x 1 x_1 x1 、求解 x 1 x_1 x1 就是求解 x 2 x_2 x2
相当于求解2个一阶微分方程
3. 奇异值分解(SVD)
3.1 简介
- 参考
- 简介：特征值分解仅适用于提取方阵特征，但在实际应用中，大部分数据对应的矩阵都不是方阵；
奇异值分解是将任意形状矩阵A用更小、更简单的 3个子矩阵乘积 U ξ V T U\xi V^T UξVT 表示，如下
A = U ξ V T A=U\xi V^T A=UξVT - SVD分解后子矩阵的特点
-
ξ
\xi
ξ 中的元素会按大到小排好顺序，因此代表该维度重要性越高(值越大)。它是
A
A
T
AA^T
AAT的
特征值的根号结果。 -
V
V
V 的列向量表示了
A
A
T
AA^T
AAT矩阵
的标准正交基(特征向量归一化后的结果，顺序需要根据 ξ \xi ξ 对应)，在使用过程中一般会进行使用转置结果 V T V^T VT； - U U U 的列向量表示变换A后的新标准正交基，一般需要等待 ξ \xi ξ、 V V V计算完成后再计算， U = A V ξ − 1 U=AV\xi ^{-1} U=AVξ−1, 若维度不够直接补充正交基就行；
-
ξ
\xi
ξ 中的元素会按大到小排好顺序，因此代表该维度重要性越高(值越大)。它是
A
A
T
AA^T
AAT的
- (补充知识)
- 对于任何一个矩阵
A
A
A，它的
A
A
T
AA^T
AAT 一定是一个
对称矩阵。然而对称矩阵的特征根一定是实数，且一定可以相似对角化和正交对角化，即 A A T = P D P − 1 AA^T=PDP^{-1} AAT=PDP−1， D D D 是 A A T AA^T AAT 的特征值， P P P 是 A A T AA^T AAT 的正交基
- 对于任何一个矩阵
A
A
A，它的
A
A
T
AA^T
AAT 一定是一个
3.1 应用1-图像压缩
- 代码
import numpy as np import numpy.linalg as la import matplotlib.pyplot as plt from sklearn import datasets from skimage import io # -------------------------定义一些函数----------------------------------- def getImgAsMatFromFile(filename): img = io.imread(filename, as_gray=True) return np.mat(img) def plotImg(imgMat): plt.imshow(imgMat, cmap=plt.cm.gray) plt.show() def recoverBySVD(imgMat, k): # 奇异值分解，这里的结果V已经经历了转置 U, s, V = la.svd(imgMat) # 选取重要的异常值 Uk = U[:, 0:k] Sk = np.diag(s[0:k]) Vk = V[0:k, :] # 相关信息的输出 print("当前U形状:", Uk.shape, " S形状:", Sk.shape, " V形状:", Vk.shape," 占据空间:", Uk.shape[0]*Uk.shape[1]+Sk.shape[0]*Sk.shape[1]+ Vk.shape[0]*Vk.shape[1]) print("原图占空间：",imgMat.shape,"=",imgMat.shape[0]*imgMat.shape[1]) print("压缩率:", (1-(Uk.shape[0]*Uk.shape[1]+Sk.shape[0]*Sk.shape[1]+ Vk.shape[0]*Vk.shape[1])/(imgMat.shape[0]*imgMat.shape[1]))*100 , "%") # 使用矩阵信息恢复图像 imgMat_new = Uk @ Sk @ Vk return imgMat_new # ------------------------- 进行操作 ----------------------------------- A = getImgAsMatFromFile('/home/txz/Downloads/python_demo/data/traffic.png') plotImg(A) A_new = recoverBySVD(A, 50) plotImg(A_new) A_new = recoverBySVD(A, 2) plotImg(A_new) - 结果可视化
4. LM的非线性方程求解和数据拟合（Python实现）
4.1 非线性方程的求解
- 求解这个方程
0 = x 2 + y − 5 0 = 4 ∗ x + y 2 − 9 0 = ( z − 3 ) 2 0 =x^2+y-5 \\ 0=4*x+y^2-9 \\ 0=(z-3)^2 0=x2+y−50=4∗x+y2−90=(z−3)2 - Python代码
-
使用到的所有库
import numpy as np -
实现的类
class LM_Opt: def __init__(self) -> None: self.u = 1.0 self.v = 2.0 # 论文中 u0 缩放因子， τ 属于 [1e-8, 1] self.scale_for_u0 = 1 # 获取此时变量的关于残差函数的残差值，是一个向量 def get_residual_vector_value(self, var): # 获取输入， xyz也就是需要求解的未知数 x, y, z = var # 残差函数 f1 = x*x + y - 5 f2 = 4*x + y*y - 9 f3 = (z-3)**2 # 构建残差向量 return np.array([f1, f2, f3], dtype=np.float32) # 获取残差向量函数的雅可比矩阵，我这里使用了数值微分进行计算 def get_residual_func_Jocabi(self, var, step=0.001): x, y, z = var # 这里使用数值求导的方法进行雅可比求导，当然也可以使用 vec1 = (self.get_residual_vector_value( [x + step, y , z ] ) - self.get_residual_vector_value( [x, y, z ] )) / step + 0.00001 vec2 = (self.get_residual_vector_value( [x , y + step, z ] ) - self.get_residual_vector_value( [x, y, z ] )) / step + 0.00001 vec3 = (self.get_residual_vector_value( [x , y , z + step ] ) - self.get_residual_vector_value( [x, y, z ] )) / step + 0.00001 return np.array([vec1, vec2, vec3], dtype=np.float32).transpose() # 进行求解，设置初始估计，最大迭代次数 def run(self, start_value, max_step = 1000): x, y, z = start_value start_flag = True for iter_indx in range(max_step): # 算法的初始化 if start_flag : cur_xyz = [x, y, z] Jacobi = self.get_residual_func_Jocabi(var = cur_xyz) # 初始化 u self.u = self.scale_for_u0 * max( *((Jacobi.transpose() @ Jacobi).diagonal()), 1) start_flag = False # 每个残差函数 构成的 残差向量 residual_vec = self.get_residual_vector_value(var = cur_xyz) # 残差向量函数 关于变量 的雅可比矩阵 resudual_Jacobi = self.get_residual_func_Jocabi(var = cur_xyz) # LM 中的 Hessian 矩阵，也就是残差关于变量的二阶导 hessian_mat = resudual_Jacobi.transpose() @ resudual_Jacobi + self.u * np.eye( resudual_Jacobi.shape[1] ) # 计算每一个残差的平方和 F_func = (residual_vec @ residual_vec) # / len(residual_vec) # X 期望增量 delta_x = - np.linalg.inv(hessian_mat) @ resudual_Jacobi.transpose() @ residual_vec # 更新 x， 并计算 新x 对应的 F值 cur_xyz = cur_xyz + delta_x residual_vec_x_update = self.get_residual_vector_value( cur_xyz ) F_func_x_update = residual_vec_x_update @ residual_vec_x_update # 计算 p(读作 rou)， 然后由此更新下一轮的 u p = (F_func - F_func_x_update) / (0.5 * delta_x.transpose() @ (self.u * delta_x - resudual_Jacobi.transpose()@residual_vec) ) if p > 0: self.u = self.u * max([1/3, (1 - pow((2*p-1),3) )]) self.v = 2 else: self.u = self.u * self.v self.v = 2 * self.v print(iter_indx, "当前残差绝对值:{: <25},\t x 最大变化量绝对值：{: < 25}, \t 得到最优x结果:{} ".format(sum( abs(residual_vec)), max( abs(delta_x) ), cur_xyz) ) # 结束标志 if abs(sum(abs(residual_vec))) < 1e-8 or max( abs(delta_x) ) < 1e-8: print("".center(150, "=")) print("优化结束！ 当前残差:{}, 得到结果:{} ".format(sum(abs(residual_vec)), cur_xyz) ) break -
进行运行
opt = LM_Opt() opt.run(start_value=[1.9, 0.5, 2], max_step=1000) -
运行结果
4.2 数据拟合求解参数
- 假如我有一堆
x
、
y
x、y
x、y 的数据，它是关于函数
y
=
e
x
p
0.1
x
+
5
y=exp^{0.1x}+5
y=exp0.1x+5的数据，采集的结果
y
y
y 有白噪声。现在需要根据采集的
x
、
y
x、y
x、y数据求解下面方程的未知数。其中
a
、
b
a、b
a、b 是未知数。
y = e x p a x + 5 ∗ b y=exp^{ax}+5*b y=expax+5∗b - Pyhton实现代码，基于LM
- 导入所有使用到的库
import numpy as np - 定义一个类，包含生成模拟数据
class LM_Opt: def __init__(self) -> None: self.u = 1.0 self.v = 2.0 # 论文中 u0 缩放因子， τ 属于 [1e-8, 1] self.scale_for_u0 = 1 # 假如这是采集的数据 def generate_data(self, data_num = 500): # 这是真实的函数, 需要拟合的数据就是 exp(a*x) + 5*b func = lambda x:np.exp( 0.1*x ) + 5*1 # 0到50，产生多少个数据 x = np.linspace(0, 50, data_num) y = func(x) # 生成标准正态分布的数据，期望为0 np.random.seed(1) noise = np.random.random(data_num) - 0.5 y = y+noise*10 return [x, y] # 获取此时变量的关于残差函数的残差值，是一个向量 def get_residual_vector_value(self, var, fine_xy_data): # 获取输入， xyz也就是需要求解的未知数 a, b = var data_for_x, data_for_y = fine_xy_data[0], fine_xy_data[1] # 残差函数, 这里面包含了待拟合的函数 f1 = np.exp( a * data_for_x) + 5*b - data_for_y # 构建残差向量 return f1 # 获取残差向量函数的雅可比矩阵，我这里使用了数值微分进行计算 def get_residual_func_Jocabi(self, var, step=0.001, fine_xy_data=None): a, b = var # 这里使用数值求导的方法进行雅可比求导，当然也可以使用 vec1 = (self.get_residual_vector_value([a + step, b ], fine_xy_data=fine_xy_data) - self.get_residual_vector_value([a, b] , fine_xy_data = fine_xy_data)) / step + 0.00001 vec2 = (self.get_residual_vector_value([a, b+step], fine_xy_data=fine_xy_data) - self.get_residual_vector_value([a, b] , fine_xy_data = fine_xy_data)) / step + 0.00001 return np.array([vec1, vec2], dtype=np.float32).transpose() # 进行拟合 def run(self, init_value, max_step = 1000): # 生成数据 fine_xy_data = self.generate_data(500) print("采集得到的数据XY, 共计：", len(fine_xy_data[0])) # 初始解 a, b = init_value print("参数拟合的初始值：", init_value) start_flag = True for iter_indx in range(max_step): # 算法的初始化 if start_flag : cur_ab = [a, b] Jacobi = self.get_residual_func_Jocabi(var = cur_ab, fine_xy_data=fine_xy_data) # 初始化 u self.u = self.scale_for_u0 * max( *((Jacobi.transpose() @ Jacobi).diagonal()), 1) start_flag = False # 每个残差函数 构成的 残差向量 residual_vec = self.get_residual_vector_value(var = cur_ab, fine_xy_data=fine_xy_data) # 残差向量函数 关于变量 的雅可比矩阵 resudual_Jacobi = self.get_residual_func_Jocabi(var = cur_ab, fine_xy_data=fine_xy_data) # LM 中的 Hessian 矩阵，也就是残差关于变量的二阶导 hessian_mat = resudual_Jacobi.transpose() @ resudual_Jacobi + self.u * np.eye( resudual_Jacobi.shape[1] ) # 计算每一个残差的平方和 F_func = (residual_vec @ residual_vec) # / len(residual_vec) # X 期望增量 delta_x = - np.linalg.inv(hessian_mat) @ resudual_Jacobi.transpose() @ residual_vec # 更新 x， 并计算 新x 对应的 F值 cur_ab = cur_ab + delta_x residual_vec_x_update = self.get_residual_vector_value( cur_ab, fine_xy_data= fine_xy_data) F_func_x_update = residual_vec_x_update @ residual_vec_x_update # 计算 p(读作 rou)， 然后由此更新下一轮的 u p = (F_func - F_func_x_update) / (0.5 * delta_x.transpose() @ (self.u * delta_x - resudual_Jacobi.transpose()@residual_vec) ) if p > 0: self.u = self.u * max([1/3, (1 - pow((2*p-1),3) )]) self.v = 2 else: self.u = self.u * self.v self.v = 2 * self.v print(iter_indx, "当前残差绝对值:{: <25},\t x 最大变化量绝对值：{: < 25}, \t 得到最优x结果:{} ".format(sum( abs(residual_vec)), max( abs(delta_x) ), cur_ab) ) # 结束标志 if abs(sum(abs(residual_vec))) < 1e-8 or max( abs(delta_x) ) < 1e-8: print("".center(150, "=")) print("优化结束！ 当前残差:{}, 得到结果:{} ".format(sum(abs(residual_vec)), cur_ab) ) break - 调用并运行
opt = LM_Opt() opt.run(init_value=[0.5, 2], max_step=1000) - 最终结果
扫一扫
4196
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
