非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (实例篇 V)

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)=21r(x)22=21i=13ri(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)=sinx10.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.4rat: 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"]
maxima_output
Fig. 1 Maxima 公式推导输出
funtion_img
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.shapefor 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 Celse: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 = 0array_x_1 = []array_x_2 = []array_x_3 = []new_x = np.matrix([[0],[0]])while iter < max_iter:array_x_1.append(x_1)array_x_2.append(x_2)array_x_3.append(g(x_1, x_2))# print("iter: %d"%iter)# print("x_1: %f"%x_1)# print("x_2: %f"%x_2)# print("g: %f"%g(x_1, x_2))# print("sH:")# print(sH(x_1, x_2))# print("inv(sH):")# print(inv(sH(x_1, x_2)))sH_i = sH(x_1, x_2)if det(sH_i) == 0:print("The simplified Hessian matrix in Guass-Newton Method is singular.")breakelse:inv_sH_i = inv(sH_i)dg_i = dg(x_1, x_2)new_x = np.matrix([[x_1], [x_2]]) - multiply_matrix(inv_sH_i, dg_i)update_distance = sqrt(pow(x_1-new_x[0,0], 2) + pow(x_2-new_x[1,0], 2))# print("update_distance: %f"%update_distance)if update_distance < epsilon:  breakiter += 1x_1 = new_x[0,0]x_2 = new_x[1,0]return array_x_1, array_x_2, array_x_3def result_plot(trajectory):fig = plt.figure()ax3 = plt.axes(projection='3d')xx = np.arange(-5,5,0.1)yy = np.arange(-4,4,0.1)X, Y = np.meshgrid(xx, yy)Z = np.zeros((X.shape[0], Y.shape[1]), dtype = float)for i in range(X.shape[0]):for j in range(Y.shape[1]):Z[i,j] = g(X[0,j], Y[i,0])ax3.plot_surface(X, Y, Z, rstride = 1, cstride = 1, cmap='rainbow', alpha=0.25)ax3.contour(X, Y, Z, offset=-1, cmap = 'rainbow')ax3.plot(trajectory[0], trajectory[1], trajectory[2], "r--")offset_data = -1*np.ones(len(trajectory[0]))ax3.plot(trajectory[0], trajectory[1], offset_data,'k--')ax3.set_title('Guass-Newton Method (Initial point [%.1f, %.1f])' %(trajectory[0][0], trajectory[1][0]))file_name_prefix = "./gauss_newton"file_extension = ".png"file_name = f"{file_name_prefix}_{trajectory[0][0]}_{trajectory[1][0]}{file_extension}"print(file_name)plt.draw()plt.savefig(file_name)if __name__ == "__main__":test_data = np.array([[4.9, 3.9], [-2.9, 1.9], [0.1, -0.1], [-0.1, 0.1], [0,-3.8],[1,2.5], [0,0]])for inital_data in test_data:print("\nInitial point:")print(inital_data)x_1 = inital_data[0]x_2 = inital_data[1]epsilon = 1e-10max_iter = 1000trajectory = gauss_newton_method(x_1, x_2, epsilon, max_iter)result_plot(trajectory)

得到的结果:

测试显示测试显示
gauss_newton_4.9_3.9gauss_newton_-2.9_1.9
gauss_newton_0.1_-0.1gauss_newton_-0.1_0.1
gauss_newton_0.0_-3.8gauss_newton_1.0_2.5

4. 结论

围绕一个最小二乘问题实例, 利用 Maxima 进行了公式推导, 利用 Python 进行了数值计算和结果显示.

(如有问题请指出, 谢谢! )

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/229461.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

Docker 从入门到实践:Docker介绍

前言 在当今的软件开发和部署领域&#xff0c;Docker已经成为了一个不可或缺的工具。Docker以其轻量级、可移植性和标准化等特点&#xff0c;使得应用程序的部署和管理变得前所未有的简单。无论您是一名开发者、系统管理员&#xff0c;还是IT架构师&#xff0c;理解并掌握Dock…

CSS 纵向底部往上动画

<template><div class"container" mouseenter"startAnimation" mouseleave"stopAnimation"><!-- 旋方块 --><div class"box" :class"{ scale-up-ver-bottom: isAnimating }"><!-- 元素内容 --&g…

c++_STL容器总结

STL容器总结 1.STL的基本概念1.2STL的六大组件 2.string类2.1string的基本概念2.2string容器常用操作 3.vector容器3.1vector容器基本概述 4.deque容器4.1deque容器的基本概念4.2deque容器的实现原理4.3deque常用API 5. stack容器5.2stack常用API 6.queue容器6.1 queue 容器基本…

electron进程通信之预加载脚本和渲染进程对主进程通信

主进程和预加载脚本通信 主进程 mian,js 和预加载脚本preload.js,在主进程中创建预加载脚本, const createWindow () > {// Create the browser window.const mainWindow new BrowserWindow({width: 300,height: 300,// 指定预加载脚本webPreferences: {preload: path.j…

dns主从搭建测试

一、DNS的介绍 1、DNS&#xff1a;Domain Name System&#xff0c;域名系统。将主机名解析为IP地址的过程&#xff0c;完成从域名到主机识别ip地址之间的转换&#xff0c;如&#xff1a;www.baidu.com, 其中 www为主机名&#xff0c;baidu.com为域名。 2、DNS无论是走TCP,还是走…

EBDP:解锁大数据的奥秘✨

大数据时代已经来临&#xff0c;你是否也想掌握这门“显学”&#xff1f;&#x1f31f; EBDP&#xff0c;这个让众多专业人士趋之若鹜的认证&#xff0c;究竟有何魅力&#xff1f;今天就带你一探究竟&#xff01; &#x1f31f;EBDP&#xff1a;大数据的“敲门砖”&#x1faa…

算法28:力扣64题,最小路径和------------样本模型

题目&#xff1a; 给定一个二维数组matrix&#xff0c;一个人必须从左上角出发&#xff0c;最后到达右下角 。沿途只可以向下或者向右走&#xff0c;沿途的数字都累加就是距离累加和 * 返回累加和最小值 思路&#xff1a; 1. 既然是给定二维数组matrix&#xff0c;那么二维数…

element el-table实现可进行横向拖拽滚动

【问题】表格横向太长&#xff0c;表格横向滚动条位于最底部&#xff0c;需将页面滚动至最底部才可左右拖动表格&#xff0c;用户体验感不好 【需求】基于elment的el-table组件生成的表格&#xff0c;使其可以横向拖拽滚动 【实现】灵感来源于这篇文章【Vue】表格可拖拽滚动&am…

C++摸版(初阶)----函数模版与类模版

本专栏内容为&#xff1a;C学习专栏&#xff0c;分为初阶和进阶两部分。 通过本专栏的深入学习&#xff0c;你可以了解并掌握C。 &#x1f493;博主csdn个人主页&#xff1a;小小unicorn ⏩专栏分类&#xff1a;C &#x1f69a;代码仓库&#xff1a;小小unicorn的代码仓库&…

大数据StarRocks(二) StarRocks集群部署

一、生产机器资源评估 1.梳理数据量&#xff0c;包括每天增量数据接入和全量数据接入 2.数据存储时间长度&#xff08;1个月/3个月/半年/1年/三年等&#xff09; 3.报表的SQL查询数量&#xff0c;SQL查询占用资源的统计&#xff0c;需要提前做好压测 4.压测可以采用官网提供的…

C++Qt6 多种排序算法的比较 数据结构课程设计 | JorbanS

一、 问题描述 在计算机科学与数学中&#xff0c;一个排序算法&#xff08;英语&#xff1a;Sorting algorithm&#xff09;是一种能将一串资料依照特定排序方式排列的算法。最常用到的排序方式是数值顺序以及字典顺序。有效的排序算法在一些算法&#xff08;例如搜索算法与合…

STM32F407-14.3.10-表73具有有断路功能的互补通道OCx和OCxN的输出控制位-01x00

如上表所示&#xff0c;MOE0&#xff0c;OSSI1&#xff0c;CCxE0&#xff0c;CCxNE0时&#xff0c;OCx与OCxN的输出状态取决于GPIO端口上下拉状态。 ---------------------------------------------------------------------------------------------------------------------…

web前端开发网页制作html/css结课作业

效果图展示&#xff1a; 注意事项&#xff1a; 引用JQuery文件地址和图片地址要更换一下。 百度网盘链接&#xff1a; http://链接&#xff1a;https://pan.baidu.com/s/1wYkmLr7csjBwQY6GmlYm4Q?pwd4332 提取码&#xff1a;4332 html界面展示&#xff1a; main.css代码部…

JavaScript:作用域变量回收

JavaScript&#xff1a;作用域&变量回收 局部作用域函数作用域块作用域 全局作用域作用域链变量在浏览器模型中的位置浏览器模型全局变量的产生情况直接赋值全局对象与var全局对象的区别 垃圾回收机制引用计数法标记清除法 闭包变量提升&函数提升 作用域规定了变量能够…

Zookeeper之Java客户端实战

ZooKeeper应用的开发主要通过Java客户端API去连接和操作ZooKeeper集群。可供选择的Java客户端API有&#xff1a; ZooKeeper官方的Java客户端API。第三方的Java客户端API&#xff0c;比如Curator。 接下来我们将逐一学习一下这两个java客户端是如何操作zookeeper的。 1. ZooKe…

[DAU-FI Net开源 | Dual Attention UNet+特征融合+Sobel和Canny等算子解决语义分割痛点]

文章目录 概要I Introduction小结 概要 提出的架构&#xff0c;双注意力U-Net与特征融合&#xff08;DAU-FI Net&#xff09;&#xff0c;解决了语义分割中的挑战&#xff0c;特别是在多类不平衡数据集上&#xff0c;这些数据集具有有限的样本。DAU-FI Net 整合了多尺度空间-通…

C#使用 OpenHardwareMonitor获取CPU或显卡温度、使用率、时钟频率相关方式

C# 去获取电脑相关的基础信息&#xff0c;还是需要借助 外部的库&#xff0c;我这边尝试了自己去实现它 网上有一些信息&#xff0c;但不太完整&#xff0c;都比较零碎&#xff0c;这边尽量将代码完整的去展示出来 OpenHardwareMonitor获取CPU的温度和频率需要管理员权限 在没…

opencv003图像裁剪(应用NumPy矩阵的切片)

这一部分相对于马上要学习的二值化是要更更更简单一些的&#xff0c;只需三行&#xff0c;便能在opencv上裁剪图像啦&#xff08;顺便云吸猫&#xff0c;太可爱了&#xff01;&#xff09; 出处见水印&#xff01; 1、复习图像的显示 前几天期末考试&#xff0c;太久没有看…

Unity中Shader的Reversed-Z(DirectX平台)

文章目录 前言一、在对裁剪坐标归一化设置NDC时&#xff0c;DirectX平台Z的特殊二、在图形计算器中&#xff0c;看一下Z值反转前后变化1、在图形计算器创建两个变量 n 和 f 分别 控制近裁剪面 和 远裁剪面2、带入公式得到齐次裁剪空间下Z值3、进行透视除法4、用 1 - Z 得出Z值反…

是否需要跟上鸿蒙(OpenHarmony)开发岗位热潮?

前言 自打华为2019年发布鸿蒙操作系统以来&#xff0c;网上各种声音百家争鸣。尤其是2023年发布会公布的鸿蒙4.0宣称不再支持Android&#xff0c;更激烈的讨论随之而来。 本文没有宏大的叙事&#xff0c;只有基于现实的考量。 通过本文&#xff0c;你将了解到&#xff1a; Har…