当前位置: 首页 > 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…...

12.找到字符串中所有字母异位词

&#x1f9e0; 题目解析 题目描述&#xff1a; 给定两个字符串 s 和 p&#xff0c;找出 s 中所有 p 的字母异位词的起始索引。 返回的答案以数组形式表示。 字母异位词定义&#xff1a; 若两个字符串包含的字符种类和出现次数完全相同&#xff0c;顺序无所谓&#xff0c;则互为…...

(转)什么是DockerCompose?它有什么作用?

一、什么是DockerCompose? DockerCompose可以基于Compose文件帮我们快速的部署分布式应用&#xff0c;而无需手动一个个创建和运行容器。 Compose文件是一个文本文件&#xff0c;通过指令定义集群中的每个容器如何运行。 DockerCompose就是把DockerFile转换成指令去运行。 …...

docker 部署发现spring.profiles.active 问题

报错&#xff1a; org.springframework.boot.context.config.InvalidConfigDataPropertyException: Property spring.profiles.active imported from location class path resource [application-test.yml] is invalid in a profile specific resource [origin: class path re…...

管理学院权限管理系统开发总结

文章目录 &#x1f393; 管理学院权限管理系统开发总结 - 现代化Web应用实践之路&#x1f4dd; 项目概述&#x1f3d7;️ 技术架构设计后端技术栈前端技术栈 &#x1f4a1; 核心功能特性1. 用户管理模块2. 权限管理系统3. 统计报表功能4. 用户体验优化 &#x1f5c4;️ 数据库设…...

Python基于历史模拟方法实现投资组合风险管理的VaR与ES模型项目实战

说明&#xff1a;这是一个机器学习实战项目&#xff08;附带数据代码文档&#xff09;&#xff0c;如需数据代码文档可以直接到文章最后关注获取。 1.项目背景 在金融市场日益复杂和波动加剧的背景下&#xff0c;风险管理成为金融机构和个人投资者关注的核心议题之一。VaR&…...

【SSH疑难排查】轻松解决新版OpenSSH连接旧服务器的“no matching...“系列算法协商失败问题

【SSH疑难排查】轻松解决新版OpenSSH连接旧服务器的"no matching..."系列算法协商失败问题 摘要&#xff1a; 近期&#xff0c;在使用较新版本的OpenSSH客户端连接老旧SSH服务器时&#xff0c;会遇到 "no matching key exchange method found"​, "n…...

C#学习第29天:表达式树(Expression Trees)

目录 什么是表达式树&#xff1f; 核心概念 1.表达式树的构建 2. 表达式树与Lambda表达式 3.解析和访问表达式树 4.动态条件查询 表达式树的优势 1.动态构建查询 2.LINQ 提供程序支持&#xff1a; 3.性能优化 4.元数据处理 5.代码转换和重写 适用场景 代码复杂性…...

[ACTF2020 新生赛]Include 1(php://filter伪协议)

题目 做法 启动靶机&#xff0c;点进去 点进去 查看URL&#xff0c;有 ?fileflag.php说明存在文件包含&#xff0c;原理是php://filter 协议 当它与包含函数结合时&#xff0c;php://filter流会被当作php文件执行。 用php://filter加编码&#xff0c;能让PHP把文件内容…...

Visual Studio Code 扩展

Visual Studio Code 扩展 change-case 大小写转换EmmyLua for VSCode 调试插件Bookmarks 书签 change-case 大小写转换 https://marketplace.visualstudio.com/items?itemNamewmaurer.change-case 选中单词后&#xff0c;命令 changeCase.commands 可预览转换效果 EmmyLua…...

Qt的学习(一)

1.什么是Qt Qt特指用来进行桌面应用开发&#xff08;电脑上写的程序&#xff09;涉及到的一套技术Qt无法开发网页前端&#xff0c;也不能开发移动应用。 客户端开发的重要任务&#xff1a;编写和用户交互的界面。一般来说和用户交互的界面&#xff0c;有两种典型风格&…...