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

news2024/12/23 14:57:53

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​=21​i=1∑3​ri​(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"]
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.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 < 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.")
            break
        else:
            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:  
            break
        iter += 1
        x_1 = new_x[0,0]
        x_2 = new_x[1,0]
        
    return array_x_1, array_x_2, array_x_3
      

def 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-10
        max_iter = 1000
        trajectory = 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.coloradmin.cn/o/1352606.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

RPA如何入门？

RPA如何入门？

许多人可能认为RPA等人工智能软件的入门非常困难&#xff0c;不知道该如何开始。作为一名刚接触RPA软件一年多的新手&#xff0c;回答这个问题也是对自己学习RPA历程的回顾。同时&#xff0c;这也是分享给那些打算学习RPA初学者的经验&#xff0c;希望可以帮助大家更好地、更迅…
阅读更多...
SANSAN新鲜事｜物联网平台就只能做数据监控吗？

SANSAN新鲜事｜物联网平台就只能做数据监控吗？

​物联网&#xff08;Internet of Things&#xff0c;简称IoT&#xff09;是指通过网络将各种物理设备、智能终端、传感器等连接起来&#xff0c;实现信息的交换和通信的技术。物联网平台&#xff08;IoT Platform&#xff09;是指为物联网应用提供基础设施、服务和工具的软件平…
阅读更多...
【springboot配置文件加载源码分析】

【springboot配置文件加载源码分析】

在Spring Boot的源码中&#xff0c;配置文件的加载是在应用程序启动的早期阶段进行的。具体来说&#xff0c;配置文件加载的主要步骤发生在SpringApplication类的run()方法中的prepareEnvironment方法中&#xff0c;真正读取我们的配置文件还是PropertySourceLoader。 本篇博客…
阅读更多...
RFID智能生产制造全周期管理系统

RFID智能生产制造全周期管理系统

一、MES系统简述 RFID/条码技术基于的MES制造执行系统可以加强ERP/MRP计划的执行功能。该系统包括MES与ERP的对接、MES报表与看板、MES物料管理、MES设备与工具管理、MES品质管理和MES生产过程管理WIP等功能&#xff0c;通过将MES系统与ERP计划和车间作业现场控制系统联系起来…
阅读更多...
【面试】面向对象编程的三大概念（实例辅助记忆）

【面试】面向对象编程的三大概念（实例辅助记忆）

【面试】面向对象编程的三大概念&#xff08;实例辅助记忆&#xff09; 虑面向对象编程的三大特性&#xff0c;它们是&#xff1a; 封装&#xff08;Encapsulation&#xff09;&#xff1a; 将对象的状态和行为封装在一起&#xff0c;对外部隐藏对象的内部实现细节。这样可以防…
阅读更多...
第二部分组件化编程:vue学习（53-60）

第二部分组件化编程:vue学习（53-60）

文章目录 53.对组件的理解54 非单文件组件55 组件的几个注意事项56 组件的嵌套57 vuecomponent构造函数58 vue实例与组件实例59 一个重要的内置关系60 单文件组件 53.对组件的理解 左侧2个页面&#xff0c;如果要复用js和css的话&#xff0c;引用的路线十分混乱。使用js模块化&…
阅读更多...
熔断、隔离、重试、降级、超时、限流，高可用架构流量治理核心策略全掌握

熔断、隔离、重试、降级、超时、限流，高可用架构流量治理核心策略全掌握

可用性的定义 在探讨高可用架构之前&#xff0c;让我们以 O2 系统为例&#xff0c;解释一下何谓可用性。O2 是腾讯内部的一个广告投放系统&#xff0c;专注于提升投放效率、分析广告效果&#xff0c;拥有自动化广告投放、AIGC 自动化素材生产等多种功能。 其整体架构概览如下&…
阅读更多...
abaqus复合材料 19个实例

abaqus复合材料 19个实例

实例操作: 1.复合材料层结构的三种常用建模方法、静力分析中强度准则和损伤判据的引入、数据输入与输出 2.层合结构的热-力耦合分析 3.基于虚裂纹闭合技术(VCCT)的分层扩展模拟 4.基于cohesive单元的分层/界面损伤扩展模拟 5.基于XFEM方法的裂纹扩展模拟 6.复合材料加筋板…
阅读更多...
在Go语言中实现HTTP请求的缓存

在Go语言中实现HTTP请求的缓存

大家好&#xff0c;我是你们可爱的编程小助手&#xff0c;今天我们要一起探讨如何使用Go语言实现HTTP请求的缓存。听起来是不是很酷&#xff1f;让我们开始吧&#xff01; 首先&#xff0c;我们要明白什么是缓存。简单来说&#xff0c;缓存就是将数据存储在内存中&#xff0c;…
阅读更多...
中国供应链，出海大时代

中国供应链，出海大时代

尽量优化、打通跨境电商每一个参与方的物流、商流、信息流、资金流是电商供应链出海的解题之法。这个过程中数智化便是打通这些节点的钥匙。 作者|斗斗 编辑|皮爷 出品|产业家 “速卖通加持&#xff0c;阿里国际零售商业收入同比上涨73%”“拼多多发布Q3财报同比增长94%…
阅读更多...
STM32MP157/linux驱动学习记录

STM32MP157/linux驱动学习记录

1. uboot烧录 2.linux安装nfs服务 sudo apt-get install nfs-kernel-server rpcbind安装nfs服务 在用户根目录下创建一个名为“linux”的文件夹&#xff0c;以后所有的东西都放到这个“linux”文件夹里面&#xff0c;在“linux”文件夹里面新建一个名为“nfs”的文件夹&#…
阅读更多...
文件属性信息

文件属性信息

文件的属性信息 Linux是一个基于文件的操作系统&#xff0c;因此作为文件本身也就有很多属性&#xff0c;如果想要查看某一个文件的属性有两种方式&#xff1a;命令和函数。虽然有两种方式但是它们对应的名字是相同的&#xff0c;叫做stat。另外使用file命令也可以查看文件的一…
阅读更多...
从fuzz视角看CTF堆题--qwb2023_chatting

从fuzz视角看CTF堆题--qwb2023_chatting

前言 这个题目是一个c的堆题&#xff0c;而我自己对于c的一些内存分配不太了解&#xff0c;同时也不太会c的逆向&#xff0c;硬看是没有办法了&#xff0c;所以就想能不能通过fuzz的角度去进行利用 fuzz 大概思路 函数选择 可以看到有add delete switch read listuser mes…
阅读更多...
大创项目推荐 深度学习动物识别 - 卷积神经网络 机器视觉 图像识别

大创项目推荐 深度学习动物识别 - 卷积神经网络 机器视觉 图像识别

文章目录 0 前言1 背景2 算法原理2.1 动物识别方法概况2.2 常用的网络模型2.2.1 B-CNN2.2.2 SSD 3 SSD动物目标检测流程4 实现效果5 部分相关代码5.1 数据预处理5.2 构建卷积神经网络5.3 tensorflow计算图可视化5.4 网络模型训练5.5 对猫狗图像进行2分类 6 最后 0 前言 &#…
阅读更多...
小红书如何高效引流？

小红书如何高效引流？

近年来&#xff0c;公域流量价格不断上涨&#xff0c;私域流量的优势逐渐凸显。企业正花费大量资源和成本来获取新流量&#xff0c;但与其如此&#xff0c;不如将精力放在留存和复购上&#xff0c;从而实现业绩的新增长。其中关键在于如何有效地将公域流量转化为私域流量。 然而…
阅读更多...
html5实现好看的个人博客模板源码

html5实现好看的个人博客模板源码

文章目录 1.设计来源1.1 主界面1.2 认识我界面1.3 我的文章界面1.4 我的模板界面1.5 文章内容界面 2.结构和源码2.1 目录结构2.2 源代码 源码下载 作者&#xff1a;xcLeigh 文章地址&#xff1a;https://blog.csdn.net/weixin_43151418/article/details/135368653 html5实现好看…
阅读更多...
ExecutorCompletionService详解

ExecutorCompletionService详解

本文已收录至Github&#xff0c;推荐阅读 &#x1f449; Java随想录 微信公众号&#xff1a;Java随想录 文章目录 摘要ExecutorCompletionService适用场景ExecutorCompletionService使用ExecutorCompletionService原理解析注意事项总结 摘要 ExecutorCompletionService 是Jav…
阅读更多...
LeetCode做题总结 15. 三数之和（未完）

LeetCode做题总结 15. 三数之和（未完）

不会做&#xff0c;参考了代码随想录和力扣官方题解&#xff0c;对此题进行整理。 代码思路 思想&#xff1a;利用双指针法&#xff0c;对数组从小到大排序。先固定一个数&#xff0c;找到其他两个。 &#xff08;1&#xff09;首先对数组从小到大排序。 &#xff08;2&…
阅读更多...
【网络技术】【Kali Linux】Wireshark嗅探（四）域名系统（DNS）

【网络技术】【Kali Linux】Wireshark嗅探（四）域名系统（DNS）

一、实验目的 本次实验使用wireshark流量分析工具进行网络嗅探&#xff0c;旨在了解域名系统&#xff08;DNS&#xff09;的工作原理。 二、域名系统概述 简单来说&#xff0c;域名系统&#xff08;Domain Name System, DNS&#xff09;将域名&#xff08;可以理解为“网址”…
阅读更多...
为什么要用扫码出入库？

为什么要用扫码出入库？

一、什么是扫码出入库管理系统 传统的仓库管理模式存在很多问题&#xff0c;如&#xff1a;货物积压、过期、丢失等。这些问题不仅影响了企业的正常运营&#xff0c;还给企业带来了经济损失。为了解决这些问题&#xff0c;扫码出入库管理系统应运而生。该系统采用先进的二维码…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023