Python实现两种稀疏矩阵的最小二乘法

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

最小二乘法

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

这两个函数可以求解Ax=b,或arg minx ∥Ax−b∥2,或arg minx ∥Ax−b∥2 +d2∥x−x0∥2,其中A必须是方阵或三角阵,可以有任意秩。

通过设置容忍度at ,bt,可以控制算法精度,记r=b-A为残差向量,如果Ax=b是相容的,lsqr在∥r∥⩽at∗∥A∥⋅∥x∥+bt∥b∥时终止;否则将在∥ATr∥⩽at∥A∥⋅∥r∥。

如果两个容忍度都是10−6 ,最终的∥r∥将有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 ,bt 。
  • conlim 另一个截止条件,对于最小二乘问题,conlim应该小于108,如果Ax=b是相容的,则conlim最大可以设到1012
  • iter_limint 迭代次数
  • show 如果为True,则打印运算过程
  • calc_var 是否估计(A.T@A + damp**2*I)^{-1}的对角线
  • x0 阻尼系数相关

lsqr和lsmr相比,没有maxiter参数,但多了iter_lim, calc_va参数。

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

返回值

lsmr的返回值依次为:

  • x 即Ax=b中的x
  • istop 程序结束运行的原因
  • itn 迭代次数
  • normr ∥b−Ax
  • normar ∥AT (b−Ax)∥
  • norma ∥A∥
  • conda A的条件数
  • normx ∥x∥

lsqr的返回值为

  1. x 即Ax=b中的x
  2. istop 程序结束运行的原因
  3. itn 迭代次数
  4. r1norm
  5. anorm 估计的Frobenius范数Aˉ
  6. acond Aˉ的条件数
  7. arnorm ∥ATr−δ2(x−x0)∥
  8. xnorm ∥x∥
  9. var (ATA)−1

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

测试

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

import numpy as np
from scipy.sparse import csr_array

np.random.seed(42)  # 设置随机数状态
mat = np.random.rand(500,500)
mat[mat<0.9] = 0
csr = csr_array(mat)

然后用这个稀疏矩阵乘以一个x,得到b

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 × 500 500\times500500×500这个尺寸显然已经不合适了,用2000×2000

from timeit import timeit

np.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实现两种稀疏矩阵的最小二乘法的详细内容,更多关于Python稀疏矩阵最小二乘法的资料请关注我们其它相关文章!

(0)

相关推荐

  • Python最小二乘法矩阵

    最小二乘法矩阵 #! /usr/bin/env python # -*- coding: utf-8 -*- import numpy as np def calc_left_k_mat(k): """ 获得左侧k矩阵 :param k: :return: """ k_mat = [] for i in range(k + 1): now_line = [] for j in range(k + 1): now_line.append(j + i

  • python实现稀疏矩阵示例代码

    工程实践中,多数情况下,大矩阵一般都为稀疏矩阵,所以如何处理稀疏矩阵在实际中就非常重要.本文以Python里中的实现为例,首先来探讨一下稀疏矩阵是如何存储表示的. 1.sparse模块初探 python中scipy模块中,有一个模块叫sparse模块,就是专门为了解决稀疏矩阵而生.本文的大部分内容,其实就是基于sparse模块而来的. 第一步自然就是导入sparse模块 >>> from scipy import sparse 然后help一把,先来看个大概 >>> h

  • 最小二乘法及其python实现详解

    最小二乘法Least Square Method,做为分类回归算法的基础,有着悠久的历史(由马里·勒让德于1806年提出).它通过最小化误差的平方和寻找数据的最佳函数匹配.利用最小二乘法可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小.最小二乘法还可用于曲线拟合.其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达. 那什么是最小二乘法呢?别着急,我们先从几个简单的概念说起. 假设我们现在有一系列的数据点 ,那么由我们给出的拟合函数h(x)得到的估计量就是

  • Python 如何解决稀疏矩阵运算

    用Python求解微分线性方程 因为之前用matlab也编写过,所以前不久试着用python写,感觉之间互通点也蛮多的,易理解. 题目:稀疏线性方程组的求解方法 简单的方程如: AX=b 其中 python有很多功能库,这些库对于编程很有帮助,可以在pycharm的Project Interpreter导入库,例如numpy.os.scipy等比较基础的库, 下面是用来求解的代码: import numpy as np from scipy import linalg import os #输入

  • Python实现曲线拟合的最小二乘法

    本文实例为大家分享了Python曲线拟合的最小二乘法,供大家参考,具体内容如下 模块导入 import numpy as np import gaosi as gs 代码 """ 本函数通过创建增广矩阵,并调用高斯列主元消去法模块进行求解. """ import numpy as np import gaosi as gs shape = int(input('请输入拟合函数的次数:')) x = np.array([0.6,1.3,1.64,1

  • Python实现两种稀疏矩阵的最小二乘法

    目录 最小二乘法 返回值 测试 最小二乘法 scipy.sparse.linalg实现了两种稀疏矩阵最小二乘法lsqr和lsmr,前者是经典算法,后者来自斯坦福优化实验室,据称可以比lsqr更快收敛. 这两个函数可以求解Ax=b,或arg minx ∥Ax−b∥2,或arg minx ∥Ax−b∥2 +d2∥x−x0∥2,其中A必须是方阵或三角阵,可以有任意秩. 通过设置容忍度at ,bt,可以控制算法精度,记r=b-Ax 为残差向量,如果Ax=b是相容的,lsqr在∥r∥⩽at∗∥A∥⋅∥x∥

  • python处理两种分隔符的数据集方法

    在做机器学习的时候,遇到这样一个数据集... 一共399行10列, 1-9列是用不定长度的空格分割, 第9-10列之间用'\t'分割, 前九列都是数值类型,其中第三列有若干个'?'填充的缺失值... 第十列是字符串类型,.. 部分数据截图: 之前我是用python强写的...很麻烦,代码如下: 至此,可以已平均值,填充缺失值... 今天再回顾此数据库;决定用pandas库来试试; 1,导包,用pandas.read_table导入数据集, 2,数据处理 最后输出如下: 以上这篇python处理两

  • Python之两种模式的生产者消费者模型详解

    第一种使用queue队列实现: #生产者消费者模型 其实服务器集群就是这个模型 # 这里介绍的是非yield方法实现过程 import threading,time import queue q = queue.Queue(maxsize=10) def Producer(anme): # for i in range(10): # q.put('骨头%s'%i) count = 1 while True: q.put('骨头%s'%count) print('生产了骨头',count) cou

  • 对python中两种列表元素去重函数性能的比较方法

    测试函数: 第一种:list的set函数 第二种:{}.fromkeys().keys() 测试代码: #!/usr/bin/python #-*- coding:utf-8 -*- import time import random l1 = [] leng = 10L for i in range(0,leng): temp = random.randint(1,10) l1.append(temp) print '测试列表长度为:',leng #first set last = time.

  • Python实现两种多分类混淆矩阵

    目录 1.什么是混淆矩阵 2.分类模型评价指标 3.两种多分类混淆矩阵 3.1直接打印出每一个类别的分类准确率. 3.2打印具体的分类结果的数值 4.总结 1.什么是混淆矩阵 深度学习中,混淆矩阵是ROC曲线绘制的基础,同时它也是衡量分类型模型准确度中最基本,最直观,计算最简单的方法.它可以直观地了解分类模型在每一类样本里面表现,常作为模型评估的一部分.它可以非常容易的表明多个类别是否有混淆(也就是一个class被预测成另一个class). 首先要明确几个概念: T或者F:该样本 是否被正确分类

  • python使用两种发邮件的方式smtp和outlook示例

    smtp是直接调用163邮箱的smtp服务器,需要在163邮箱中设置一下.outlook发送就是Python直接调用win32方式.调用程序outlook直接发送邮件. import win32com.client as win32 import xlrd outlook = win32.Dispatch('outlook.application') mail = outlook.CreateItem(0) receivers = ['Yutao.A.Wang@alcatel-sbell.com

  • 总结python实现父类调用两种方法的不同

    python中有两种方法可以调用父类的方法: super(Child, self).method(args) Parent.method(self, args) 我用其中的一种报了如下错误: 找不到 classobj.当我把调用改为 super(B, self).f(name) 就能正确运行,且结果正确. 分析错误 因为基类没有继承 object , 在python中,一个可以这样创建: class A: pass 也可以这样创建: class A(object): pass 这两者的区别就是:

  • Python中的is和==比较两个对象的两种方法

    Python中的is和==比较两个对象的两种方法 在Python中有两种方式比较两个对象是否相等,分别是is和==,两者之间是不同的 ==比较的是值(如同java中的equals方法) is比较的是引用(可以看作比较内存地址, 类似于java中的==) 对于: >>> n = 1 >>> n is 1 True >>> b = '1' >>> b is 1 False >>> n == b False 由于1和'1'

  • python机器学习MATLAB最小二乘法的两种解读

    目录 最小二乘法 代价函数与最小二乘法 向量到子空间的距离与最小二乘法 最小二乘法与多项式拟合 多项式拟合结果绘图: 最小二乘法与多元线性回归 多元线性回归结果绘图: 最小二乘法 大部分的最小二乘法公式推导,都是使用的 代价函数偏导 的方式来求得的,在这里首先展示如何通过代价函数求偏导的方式得到最小二乘公式,再展示李扬老师讲解的如何由向量到子空间的距离得来最小二乘法公式. 代价函数与最小二乘法 假设我们的拟合结果为: 则平方损失函数为: 平方损失函数的形式只有极小值,没有极大值,我们要使代价函数

  • python学习之第三方包安装方法(两种方法)

    这篇文章主要介绍了python学习之第三方包安装方法,最近在学习QQ空间.微博(爬虫)模拟登录,都涉及到了RSA算法.这样需要下一个RSA包(第三方包),在网上搜了好多资料,具体有以下两种方法: 第一种方法(不使用pip或者easy_install): Step1:在网上找到的需要的包,下载下来.eg. rsa-3.1.4.tar.gz Step2:解压缩该文件. Step3:命令行工具cd切换到所要安装的包的目录,找到setup.py文件,然后输入python setup.py install

随机推荐