有人会问:非线性方程组不是用Newton循环方法就能解了吗?那为什么还需要JFNK方法呢?
我觉得,JFNK方法能解决Newton循环方法中最大的计算问题,让Newton循环方法不再痛苦。
Newton循环
Newton循环方法最大的痛苦,就是计算斜率。对于简单的方程式子还好,对于大型的非线性方程组,计算其斜率是非常痛苦的,这个过程往往需要手算,可以说是体力活,就算是码代码,你也要先手推出来公式,这个公式有时候就推不出来。
打个比方:
f(x) = x * log( x ² e x \frac{x²}{e^x} exx²)
让你求导x=100处的导数,然后写出代码,想想就很困难。
而JFNK方法让我们可以不用计算这个斜率,就能把Newton循环继续下去,最终得到x的值。
在Newton循环中:
x n + 1 \mathbf{x}_{n+1} xn+1 = x n \mathbf{x}_{n} xn - f ( x n ) f ′ ( x n ) \frac{f(\mathbf{x}_{n)}}{f'(\mathbf{x}_{n})} f′(xn)f(xn)
每一步计算都是在求
x
n
+
1
\mathbf{x}_{n+1}
xn+1,此时
x
n
\mathbf{x}_{n}
xn是已知。
我们稍微变形一下就可以得到
f
′
(
x
n
)
f'(\mathbf{x}_{n})
f′(xn) * (
x
n
\mathbf{x}_{n}
xn -
x
n
+
1
\mathbf{x}_{n+1}
xn+1 ) =
f
(
x
n
)
f(\mathbf{x}_{n)}
f(xn)
简洁一点:u = A − 1 A^{-1} A−1b
其中 A =
f
′
(
x
n
)
f'(\mathbf{x}_{n})
f′(xn) 就是我们一直在说的Jacobian matrix。理解Jacobian矩阵与行列式——gwave的回答,理论上可以求但是实际上比较复杂。
u = (
x
n
\mathbf{x}_{n}
xn -
x
n
+
1
\mathbf{x}_{n+1}
xn+1 )
b =
f
(
x
n
)
f(\mathbf{x}_{n)}
f(xn)
u是未知量,b和A都是已知量。但我们不想求A,更不想求
A
−
1
A^{-1}
A−1。
那么这边就可以用Krylov子空间迭代的方法Krylov子空间迭代
根据Krylov的方法,对于我们上面u =
A
−
1
A^{-1}
A−1b的问题,有:
A − 1 A^{-1} A−1b ≈ ∑ i = 0 m − 1 \displaystyle\sum_{i=0}^{m-1} i=0∑m−1 β i \mathbf{β}_i βi A i A^i Aib = β 0 \mathbf{β}_0 β0b + β 1 \mathbf{β}_1 β1Ab + β 2 \mathbf{β}_2 β2A²b…+ β m − 1 \mathbf{β}_{m-1} βm−1 A m − 1 A^{m-1} Am−1b
可以看出来,这样处理的话等号右边就没有单独出现的A,且每个A后面都跟着一个b。
这也就意味着我们不需要去计算A的值,只需要把Ab求出来就可以。(题外话,写到这就想到torch.mx)
对于A²b,我们也只需要先算出Ab,然后再将其结果当作新的b,再算一遍Ab就可以得到A²b,其他的以此类推。
那么会有人说,你连A都算不出来,Ab又怎么算?
神奇的地方来了,经过上式的推导,我们可以得出:
Ab = f ′ ( x n ) f'(\mathbf{x}_{n}) f′(xn)b = f ( x + ϵ b ) − f ( x ) ϵ \frac{f(x+\epsilon{b}) - f(x)}{\epsilon} ϵf(x+ϵb)−f(x)
其中 ϵ {\epsilon} ϵ表示一个微小的扰动量,完全可以自己定义,前提是 ϵ {\epsilon} ϵ足够小。
这样一近似,Ab就可以通过f来计算出来的,这样也就不用计算A,也就是Jacobian-free。