非线性最小二乘问题的数值方法 |
您所在的位置：网站首页 › picard迭代法求解非线性方程组例子 › 非线性最小二乘问题的数值方法 |
非线性最小二乘问题的数值方法
|
Title: 非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (实例篇 V) 姊妹博文 非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (I) 非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (II) 非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (III) 非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (IV) ↑ \uparrow ↑ 理论部分 ↓ \downarrow ↓ 实例部分 非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (实例篇 V) ⟵ \longleftarrow ⟵ 本篇 文章目录 0.前言1. 最优问题实例2. 符号计算3. 高斯-牛顿法计算4. 结论 0.前言本篇博文作为对前述 “非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法” 的实践扩展, 理论部分参见前述博文, 此处不再重复. 1. 最优问题实例m i n i m i z e g ( x ) = 1 2 ∥ r ( x ) ∥ 2 2 = 1 2 ∑ i = 1 3 r i ( x ) 2 (I-1) {\rm minimize}\quad {g}(\mathbf{x}) = \frac{1}{2}\|\mathbf{r}(\mathbf{x})\|_2^2 = \frac{1}{2}\sum_{i=1}^{3} r_i(\mathbf{x})^2 \tag{I-1} minimizeg(x)=21∥r(x)∥22=21i=1∑3ri(x)2(I-1) 其中 x = [ x 1 , x 2 ] T \mathbf{x} = \begin{bmatrix} x_1, x_2 \end{bmatrix}^{\small\rm T} x=[x1,x2]T r ( x ) = [ r 1 ( x ) , r 2 ( x ) , r 3 ( x ) ] T \mathbf{r}(\mathbf{x}) = \begin{bmatrix} r_1(\mathbf{x}), \, r_2(\mathbf{x}) ,\,r_3(\mathbf{x}) \end{bmatrix}^{\small\rm T} r(x)=[r1(x),r2(x),r3(x)]T r 1 ( x ) = sin x 1 − 0.4 r_1(\mathbf{x}) = \sin x_1 -0.4 r1(x)=sinx1−0.4 r 2 ( x ) = cos x 2 + 0.8 r_2(\mathbf{x}) = \cos x_2 + 0.8 r2(x)=cosx2+0.8 r 3 ( x ) = x 1 2 + x 2 2 − 1 r_3(\mathbf{x}) = \sqrt{x_1^2 +x_2^2} -1 r3(x)=x12+x22 −1 2. 符号计算Maxima 符号运算代码: (%i10) r_1(x_1,x_2) := sin(x_1)-0.4; r_2(x_1,x_2) :=cos(x_2)+0.8; r_3(x_1,x_2) := sqrt((x_1)^2 +(x_2)^2)-1; r: matrix([r_1(x_1,x_2)], [r_2(x_1,x_2)], [r_3(x_1,x_2)]); dr: matrix([diff(r_1(x_1,x_2), x_1), diff(r_1(x_1,x_2), x_2)], [diff(r_2(x_1,x_2), x_1), diff(r_2(x_1,x_2), x_2)], [diff(r_3(x_1,x_2), x_1), diff(r_3(x_1,x_2), x_2)]); sH: transpose(dr) . dr; dg: transpose(dr) . r; g : (1/2) * (r_1(x_1,x_2)^2+r_2(x_1,x_2)^2+r_3(x_1,x_2)^2); dg: ratexpand(matrix([diff(g,x_1)], [diff(g,x_2)])); plot3d(g, [x_1,-4,4], [x_2,-4,4]); (%o1) r_1(x_1,x_2):=sin(x_1)-0.4 (%o2) r_2(x_1,x_2):=cos(x_2)+0.8 (%o3) r_3(x_1,x_2):=sqrt(x_1^2+x_2^2)-1 (r) matrix( [sin(x_1)-0.4], [cos(x_2)+0.8], [sqrt(x_2^2+x_1^2)-1] ) (dr) matrix( [cos(x_1), 0], [0, -sin(x_2)], [x_1/sqrt(x_2^2+x_1^2), x_2/sqrt(x_2^2+x_1^2)] ) (sH) matrix( [x_1^2/(x_2^2+x_1^2)+cos(x_1)^2, (x_1*x_2)/(x_2^2+x_1^2)], [(x_1*x_2)/(x_2^2+x_1^2), sin(x_2)^2+x_2^2/(x_2^2+x_1^2)] ) (dg) matrix( [(x_1*(sqrt(x_2^2+x_1^2)-1))/sqrt(x_2^2+x_1^2)+cos(x_1)*(sin(x_1)-0.4)], [(x_2*(sqrt(x_2^2+x_1^2)-1))/sqrt(x_2^2+x_1^2)-(cos(x_2)+0.8)*sin(x_2)] ) (g) ((cos(x_2)+0.8)^2+(sqrt(x_2^2+x_1^2)-1)^2+(sin(x_1)-0.4)^2)/2 rat: replaced -0.4 by -2/5 = -0.4 rat: replaced 0.8 by 4/5 = 0.8 (dg) matrix( [-x_1/sqrt(x_2^2+x_1^2)+cos(x_1)*sin(x_1)-(2*cos(x_1))/5+x_1], [-cos(x_2)*sin(x_2)-(4*sin(x_2))/5-x_2/sqrt(x_2^2+x_1^2)+x_2] ) (%o10) ["/tmp/maxout347947.gnuplot_pipes"]
Fig. 1 Maxima 公式推导输出
Fig. 2 Maxima 中代价函数绘制
3. 高斯-牛顿法计算
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np
from numpy.linalg import inv, det
from math import cos
from math import sin
from math import sqrt
from math import pow
# multiplication of two matrixs
def multiply_matrix(A, B):
if A.shape[1] == B.shape[0]:
C = np.zeros((A.shape[0], B.shape[1]), dtype = float)
[rows, cols] = C.shape
for row in range(rows):
for col in range(cols):
for elt in range(len(B)):
C[row, col] += A[row, elt] * B[elt, col]
return C
else:
return "Cannot multiply A and B. Please check whether the dimensions of the inputs are compatible."
# g(x) = (1/2) ||r(x)||_2^2
def g(x_1, x_2):
return ( pow(sin(x_1)-0.4, 2)+ pow(cos(x_2)+0.8, 2) + pow(sqrt(pow(x_2,2)+pow(x_1,2))-1, 2) ) /2
# r(x) = [r_1, r_2, r_3]^{T}
def r(x_1, x_2):
return np.array([[sin(x_1)-0.4],
[cos(x_2)+0.8],
[sqrt(pow(x_1,2)+pow(x_2,2))-1]])
# \partial r(x) / \partial x
def dr(x_1, x_2):
if sqrt(pow(x_2,2)+pow(x_1,2)) == 0: ## 人为设置
return np.array([[cos(x_1), 0],
[0, -sin(x_2)],
[0, 0]])
else:
return np.array([[cos(x_1), 0],
[0, -sin(x_2)],
[x_1/sqrt(pow(x_2,2)+pow(x_1,2)), x_2/sqrt(pow(x_2,2)+pow(x_1,2))]])
# Simplified Hessian matrix in Gauss-Newton method, refer to eq. (IV-1-2)
def sH(x_1, x_2):
return multiply_matrix(np.transpose(dr(x_1, x_2)), dr(x_1, x_2))
# \nabla g(x_1, x_2), refer to eq. (III-1-2)
def dg(x_1, x_2):
return multiply_matrix(np.transpose(dr(x_1, x_2)), r(x_1, x_2))
def gauss_newton_method(x_1, x_2, epsilon, max_iter):
iter = 0
array_x_1 = []
array_x_2 = []
array_x_3 = []
new_x = np.matrix([[0],[0]])
while iter |
【本文地址】
今日新闻 |
推荐新闻 |