当前位置: 首页 > news >正文

Python稀疏矩阵最小二乘法

文章目录

    • 最小二乘法
    • 返回值
    • 测试

最小二乘法

scipy.sparse.linalg实现了两种稀疏矩阵最小二乘法lsqrlsmr,前者是经典算法,后者来自斯坦福优化实验室,据称可以比lsqr更快收敛。

这两个函数可以求解Ax=bAx=bAx=b,或arg min⁡x∥Ax−b∥2\argmin_x\Vert Ax-b\Vert^2argminxAxb2,或arg min⁡x∥Ax−b∥2+d2∥x−x0∥2\argmin_x\Vert Ax-b\Vert^2+d^2\Vert x-x_0\Vert^2argminxAxb2+d2xx02,其中AAA必须是方阵或三角阵,可以有任意秩。

通过设置容忍度at,bta_t, b_tat,bt,可以控制算法精度,记r=b−Axr=b-Axr=bAx为残差向量,如果Ax=bAx=bAx=b是相容的,lsqr在∥r∥⩽at∗∥A∥⋅∥x∥+bt∥b∥\Vert r\Vert\leqslant a_t*\Vert A\Vert\cdot\Vert x\Vert + b_t\Vert b\VertratAx+btb时终止;否则将在∥ATr∥⩽at∥A∥⋅∥r∥\Vert A^T r\Vert\leqslant a_t\Vert A\Vert \cdot\Vert r\VertATratAr
如果两个容忍度都是10−610^{-6}106,最终的∥r∥\Vert r\Vertr将有6位精度。

lsmr的参数如下

lsmr(A, b, damp=0.0, atol=1e-06, btol=1e-06, conlim=100000000.0, maxiter=None, show=False, x0=None)

参数解释:

  • A 可谓稀疏矩阵、数组以及线性算子
  • b 为数组
  • damp 阻尼系数,默认为0
  • atol, btol 截止容忍度,是lsqr迭代的停止条件,即at,bta_t, b_tat,bt
  • conlim 另一个截止条件,对于最小二乘问题,conlim应该小于10810^8108,如果Ax=bAx=bAx=b是相容的,则conlim最大可以设到101210^{12}1012
  • iter_limint 迭代次数
  • show 如果为True,则打印运算过程
  • calc_var 是否估计(A.T@A + damp**2*I)^{-1}的对角线
  • x0 阻尼系数相关

lsqrlsmr相比,没有maxiter参数,但多了iter_lim, calc_va参数。

上述参数中,damp为阻尼系数,当其不为0时,记作δ\deltaδ,待解决的最小二乘问题变为

[AδI]x=[bδx0]\begin{bmatrix}A\\\delta I\end{bmatrix} x=\begin{bmatrix}b\\\delta x_0 \end{bmatrix} [AδI]x=[bδx0]

返回值

lsmr的返回值依次为:

  • xAx=bAx=bAx=b中的xxx
  • istop 程序结束运行的原因
  • itn 迭代次数
  • normr ∥b−Ax∥\Vert b-Ax\VertbAx
  • normar ∥AT(b−Ax)∥\Vert A^T(b-Ax)\VertAT(bAx)
  • norma ∥A∥\Vert A\VertA
  • conda A的条件数
  • normx ∥x∥\Vert x\Vertx

lsqr的返回值为

  • xAx=bAx=bAx=b中的xxx
  • istop 程序结束运行的原因
  • itn 迭代次数
  • r1norm ∥b−Ax∥\Vert b-Ax\VertbAx
  • r2norm ∥b−Ax∥2+δ2∥x−x0∥2\sqrt{\Vert b-Ax\Vert^2+\delta^2\Vert x-x_0\Vert^2}bAx2+δ2xx02
  • anorm 估计的Frobenius范数Aˉ\bar AAˉ
  • acond Aˉ\bar AAˉ的条件数
  • arnorm ∥ATr−δ2(x−x0)∥\Vert A^Tr-\delta^2(x-x_0)\VertATrδ2(xx0)
  • xnorm ∥x∥\Vert x\Vertx
  • var (ATA)−1(A^TA)^{-1}(ATA)1

二者的返回值较多,而且除了前四个之外,剩下的意义不同,调用时且须注意。

测试

下面对这两种算法进行验证,第一步就得先有一个稀疏矩阵

import numpy as np
from scipy.sparse import csr_arraynp.random.seed(42)  # 设置随机数状态
mat = np.random.rand(500,500)
mat[mat<0.9] = 0
csr = csr_array(mat)

然后用这个稀疏矩阵乘以一个xxx,得到bbb

xs = np.arange(500)
b = mat @ xs

接下来对这两个最小二乘函数进行测试

from scipy.sparse.linalg import lsmr, lsqr
import matplotlib.pyplot as plt
mx = lsmr(csr, b)[0]
qx = lsqr(csr, b)[0]
plt.plot(xs, lw=0.5)
plt.plot(mx, lw=0, marker='*', label="lsmr")
plt.plot(qx, lw=0, marker='.', label="lsqr")
plt.legend()
plt.show()

为了对比清晰,对图像进行放大,可以说二者不分胜负

在这里插入图片描述

接下来比较二者的效率,500×500500\times500500×500这个尺寸显然已经不合适了,用2000×20002000\times20002000×2000

from timeit import timeitnp.random.seed(42)  # 设置随机数状态
mat = np.random.rand(500,500)
mat[mat<0.9] = 0
csr = csr_array(mat)
timeit(lambda : lsmr(csr, b), number=10)
timeit(lambda : lsqr(csr, b), number=10)

测试结果如下

>>> timeit(lambda : lsqr(csr, b), number=10)
0.5240591000001587
>>> timeit(lambda : lsmr(csr, b), number=10)
0.6156221000019286

看来lsmr并没有更快,看来斯坦福也不靠谱(滑稽)。

相关文章:

Python稀疏矩阵最小二乘法

文章目录最小二乘法返回值测试最小二乘法 scipy.sparse.linalg实现了两种稀疏矩阵最小二乘法lsqr和lsmr&#xff0c;前者是经典算法&#xff0c;后者来自斯坦福优化实验室&#xff0c;据称可以比lsqr更快收敛。 这两个函数可以求解AxbAxbAxb&#xff0c;或arg min⁡x∥Ax−b…...

mac本前端Homebrew下载,操作

1、打开电脑终端 2、下载Homebrew&#xff0c;在终端中输入 /usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"开始下载Homebrew&#xff0c;因为这个地址是国外网站&#xff0c;下载失败的话&#xff0c;输入…...

Linux系统之查看进程监听端口方法

Linux系统之查看进程监听端口方法一、端口监听介绍二、使用netstat命令1.netstat命令介绍2.netstat帮助3.安装netstat工具4.列出所有监听 tcp 端口5.显示TCP端口的统计信息6.查看某个服务监听端口三、使用ss命令1.ss命令介绍2.ss命令帮助3.查看某个服务监听端口四、使用lsof命令…...

使用命令别名一键启动arthas

1. 使用命令别名启动arthas 确保单板上有jdk和arthas jdk目录&#xff1a;/home/xinliushijian/arthas/jdk arthas目录&#xff1b;/home/xinliushijian/arthas su xinliushijian编写脚本messi.sh cd /home/xinliushijian/arthas vi messi.sh 内容如下&#xff1a; #!/bin/ba…...

python+pytest接口自动化(2)-HTTP协议基础

HTTP协议简介HTTP 即 HyperText Transfer Protocol&#xff08;超文本传输协议&#xff09;&#xff0c;是互联网上应用最为广泛的一种网络协议。所有的 WWW 文件都必须遵守这个标准。设计 HTTP 最初的目的是为了提供一种发布和接收 HTML 页面的方法。HTTP 协议在 OSI 模型中属…...

操作系统权限提升(十五)之绕过UAC提权-基于白名单DLL劫持绕过UAC提权

系列文章 操作系统权限提升(十二)之绕过UAC提权-Windows UAC概述 操作系统权限提升(十三)之绕过UAC提权-MSF和CS绕过UAC提权 操作系统权限提升(十四)之绕过UAC提权-基于白名单AutoElevate绕过UAC提权 注&#xff1a;阅读本编文章前&#xff0c;请先阅读系列文章&#xff0c;以…...

非常好看的html网页个人简历

一. 前言 文末获取gitee链接 在前几天逛b站的时候&#xff0c;发现了个比较实用的东西-----个人简介网页版&#xff0c;相当于网页版的个人简历&#xff0c;相较于PDF形式的&#xff0c;网页版所能呈现内容更加丰富&#xff0c;而且更加美观&#xff0c;在BOOS上被HR小姐姐要…...

轻量级网络模型ShuffleNet V2

在学习ShuffleNet V2内容前需要简单了解卷积神经网络和MobileNet,以及Shuffnet V1的相关内容&#xff0c;大家可以出门左转&#xff0c;去看我之前的几篇博客MobileNet发展脉络&#xff08;V1-V2-V3&#xff09;&#xff0c;轻量级网络模型ShuffleNet V1&#x1f197;&#xff…...

分享美容美发会员管理系统功能的特点_美容美发会员管理系统怎么做

人们越来越关心美发&#xff0c;美发行业发展迅速&#xff0c;小程序可以连接在线场景&#xff0c;许多美发院也开发了会员卡管理系统。那么一个实用的美发会员管理系统怎么制作呢&#xff1f;它有什么功能&#xff1f;我们一起来看看~&#xff08;干货满满&#xff0c;耐心看完…...

Oracle-05-DCL篇

🏆一、简介 Oracle的DCL代表数据库控制语言,用于管理数据库对象的访问和安全性。DCL的两个主要命令是GRANT和REVOKE。 GRANT命令用于授予用户或角色对数据库对象的访问权限,例如表、视图或存储过程。GRANT命令的语法如下: GRANT privilege_name [, privilege_name]... …...

tess4j简单使用入门

tess4j下载 下载地址: https://sourceforge.net/projects/tess4j/ 不要直接下载,点击files,然后下载最新版 下载解压后放到指定的目录即可,这里放到d:\jar目录下 tess4j根目录: d:\jar\tess4j tess4j使用 把test4j项目目录中dist和lib目录下的所有jar包导入到需要的项目中…...

WebGPU学习(4)---使用 UniformBuffer

接下来让我们使用 UniformBuffer。UniformBuffer 是一个只读内存区域&#xff0c;可以在着色器上访问。 这次&#xff0c;我们将传递给着色器的矩阵存储在 UniformBuffer 中。演示示例 1.在顶点着色器中的 UniformBuffer 这次我们在顶点着色器里定义一个名为Uniforms的新结构体…...

Http客户端Feign-远程调用

Feign的使用步骤 引入依赖添加EnableFeignClients注解编写FeignClient接口使用FeignClient中定义的方法代替RestTemplate Feign的日志配置 1.方式一是配置文件&#xff0c;feign.client.config.xxx.loggerLevel 如果xxx是default则代表全局如果xxx是服务名称,例如userservi…...

RK3568镜像的拆包和打包

文章目录 前言一、window上分包和打包分包打包二、Linux上分包和打包分包打包总结前言 本文记录在win10上利用瑞芯微提供的工具进行分包和打包,同样也有Linux教程 提示:以下是本篇文章正文内容,下面案例可供参考 一、window上分包和打包 分包 window下一般直接利用工具即…...

《设计模式》适配器模式

《设计模式》适配器模式 适配器&#xff08;Adapter&#xff09;是一种结构型设计模式&#xff0c;它允许我们将一个类的接口转换成另一个类的接口&#xff0c;从而使得原本由于接口不兼容而无法合作的类能够一起工作。适配器模式通常用于以下情况&#xff1a; 在已有的类中添…...

linux 随笔 5-服务管理

0. 装到虚拟机与物理机&#xff0c;感觉各有各的不方便 Linux下systemctl命令和service、chkconfig命令的区别 1. service 根据/etc/init.d目录下的配置&#xff0c;做服务相关的&#xff1a; 启动停止重新启动关闭系统服务 2. chkconfig 用于维护 /etc/rc[0-6].d 的命令…...

【java基础】枚举类(enum)

文章目录基本介绍快速使用字段、方法、构造器枚举类方法toString方法valueOf方法values方法ordinal方法基本介绍 在java中有一种特殊的类型就是枚举类&#xff0c;对于一个有限的有固定值的集合&#xff0c;我们就可以考虑使用枚举类来进行表示&#xff0c;例如服装的大小为 小…...

Linux2

&#xff08;1&#xff09;root用户的主目录&#xff1a; &#xff08;3&#xff09;查看 &#xff08;4&#xff09;远程登陆系统&#xff1a;CentOS7上使用ifconfig查看IP&#xff0c;使用putty远程登陆 &#xff08;5&#xff09;查询目前用户登录情况&#xff1a;who命令…...

C语言基础应用(二)数据的转换与输入输出

学习了C语言的基本数据类型后&#xff0c;我们可能会想这些数据如何进行运算&#xff0c;是否可以让不同类型的数据直接进行运算呢&#xff1f; 一、数据类型转换 1.1 int类型与float类型之间的转换 int i 5; // j值为2.000000 因为左右操作数均为整型float j i/2; // …...

C# 用NPOI读取EXCEL

1. 复制DLL文件 ICSharpCode.SharpZipLib.dll NPOI.dll NPOI.OOXML.dll NPOI.OpenXml4Net.dll NPOI.OpenXmlFormats.dll 2. 在工程中添加引用 3. using System.IO; using NPOI.HSSF.UserModel; using NPOI.XSSF.UserModel; using NPOI.SS.UserModel; using NPOI.OpenXml4Ne…...

《高性能MySQL》——MySQL基准测试(笔记)

文章目录二、MySQL基准测试2.1 为什么需要基准测试2.2 基准测试的策略2.2.1 测试何种指标2.3 基准测试方法2.3.1设计和规划基准测试2.3.2 基准测试应该运行多长时间2.3.3 获取系统性能和状态2.3.4 获得准确的测试结果2.3.5 运行基准测试并分析结果2.3.6 绘图的重要性2.4 基准测…...

微服务注册到Nacos后如何读取外网IP

背景 微服务部署后&#xff0c;各服务是需要相互间调用的&#xff0c;其中服务A在去调用服务B的时候发现无法调用成功。其中服务注册和发现中心以及配置中心使用的是Nacos。Nacos客户端在注册服务时会从机器网卡中选择其中一个IP来注册&#xff0c;当我们要部署的机器存在多个…...

【华为OD机试模拟题】用 C++ 实现 - 匿名信(2023.Q1)

最近更新的博客 【华为OD机试模拟题】用 C++ 实现 - 货币单位换算(2023.Q1) 【华为OD机试模拟题】用 C++ 实现 - 选座位(2023.Q1) 【华为OD机试模拟题】用 C++ 实现 - 停车场最大距离(2023.Q1) 【华为OD机试模拟题】用 C++ 实现 - 重组字符串(2023.Q1) 【华为OD机试模…...

nginx配置https域名

如果已安装nginx若没有配置ssl模块需要添加ssl模块。注意&#xff1a;默认是不会安装ssl模块的查看nginx是否有安装ssl模块./nginx -V若执行结果显示下图中的-with-http_ssl_module则说明nginx已经安装了ssl模块若没有安装则重新打包编译./configure --with-http_ssl_module编译…...

c语言tips-大端小端存储介绍和使用union判断大小端

1. 大小端介绍 大端&#xff08;Big Endian&#xff09;和小端&#xff08;Little Endian&#xff09;是两种CPU或者计算机系统存储数据的方式。 在大端系统中&#xff0c;数据的高位字节&#xff08;MSB&#xff09;存储在内存地址的低位&#xff0c;低位字节&#xff08;LSB…...

DevOps落地与转型:提升研发效能的方法与实践

❤️作者主页&#xff1a;小虚竹 ❤️作者简介&#xff1a;大家好,我是小虚竹。Java领域优质创作者&#x1f3c6;&#xff0c;CSDN博客专家&#x1f3c6;&#xff0c;华为云享专家&#x1f3c6;&#xff0c;掘金年度人气作者&#x1f3c6;&#xff0c;阿里云专家博主&#x1f3…...

【测试】loadrunner安装

努力经营当下&#xff0c;直至未来明朗&#xff01; 文章目录备注一、下载安装包二、安装loadrunner三、修改浏览器配置今天搬砖不努力&#xff0c;明天地位不稳定&#xff01; 备注 电脑最好有IE浏览器&#xff0c;但是没有也没事儿。&#xff08;注意&#xff1a;IE浏览器不…...

物联网的新应用--触摸物联网

摘要&#xff1a;本文介绍一下触摸物联网的新进展--电子皮肤的物联网应用。还以为物联网的作用领域单单是从现场采集数据或者传输命令到执行设备吗&#xff1f;不&#xff0c;物联网的应用范围远比控制一盏灯大的多。据网上报道&#xff0c;香港城市大学&#xff08;城大&#…...

用 Python 画如此漂亮的插图 ,So easy

人生苦短&#xff0c;快学Python&#xff01; 今天我们进行一次实战案例分享&#xff0c;以全球预期寿命与人均 GPD数据为例&#xff0c;写一篇 Python 中漂亮散点图的快速指南。除了正常的数据清洗/处理、还会进行简单的统计分析&#xff0c;实现数据处理-统计分析-可视化一条…...

vue-知识点总结

历史 2015年 10月27 1.0.0 Evangelion 新福音战士2016年 10月1日 2.0.0 Ghost in the Shell 攻壳机动队2019年 2月 2.62020年 9月18日 3.0.0 One Piece(海贼王)2021年 8月10日 3.2.0 .sync where 父向子传递props, 需要双向绑定的时候(子组件想更新这个值) how <one-comp…...