python numpy逆_python-具有numpy的N * M * M张量的矢量化(部分)逆
所以我有一个张量的索引为a [n,i,j]的维度为(N,M,M)的张量,我想为N中的每个n求M * M方阵部分的值.
例如,假设我有
In [1]: a = np.arange(12)
a.shape = (3,2,2)
a
Out[1]: array([[[ 0, 1],
[ 2, 3]],
[[ 4, 5],
[ 6, 7]],
[[ 8, 9],
[10, 11]]])
然后for循环反转将如下所示:
In [2]: inv_a = np.zeros([3,2,2])
for m in xrange(0,3):
inv_a[m] = np.linalg.inv(a[m])
inv_a
Out[2]: array([[[-1.5, 0.5],
[ 1. , 0. ]],
[[-3.5, 2.5],
[ 3. , -2. ]],
[[-5.5, 4.5],
[ 5. , -4. ]]])
根据github上的this issue,这显然将在NumPy 2.0中实现.
我猜我需要按照github问题线程中的seberg的说明安装开发版本,但是现在还有另一种矢量化的方法吗?
解决方法:
更新:
在NumPy 1.8和更高版本中,numpy.linalg中的函数是通用通用函数.
这意味着您现在可以执行以下操作:
import numpy as np
a = np.random.rand(12, 3, 3)
np.linalg.inv(a)
这将反转每个3×3数组,并将结果作为12x3x3数组返回.
参见numpy 1.8 release notes.
原始答案:
由于N相对较小,我们如何手动计算所有矩阵的LU分解.
这确保了所涉及的for循环相对较短.
这是使用常规NumPy语法可以完成的方法:
import numpy as np
from numpy.random import rand
def pylu3d(A):
N = A.shape[1]
for j in xrange(N-1):
for i in xrange(j+1,N):
#change to L
A[:,i,j] /= A[:,j,j]
#change to U
A[:,i,j+1:] -= A[:,i,j:j+1] * A[:,j,j+1:]
def pylusolve(A, B):
N = A.shape[1]
for j in xrange(N-1):
for i in xrange(j+1,N):
B[:,i] -= A[:,i,j] * B[:,j]
for j in xrange(N-1,-1,-1):
B[:,j] /= A[:,j,j]
for i in xrange(j):
B[:,i] -= A[:,i,j] * B[:,j]
#usage
A = rand(1000000,3,3)
b = rand(3)
b = np.tile(b,(1000000,1))
pylu3d(A)
# A has been replaced with the LU decompositions
pylusolve(A, b)
# b has been replaced to the solutions of
# A[i] x = b[i] for each A[i] and b[i]
如我所写,pylu3d修改A来计算LU分解.
用LU分解替换每个NxN矩阵后,可以使用pylusolve求解代表矩阵系统右侧的MxN数组b.
它在适当位置修改b并进行适当的后替换以解决系统问题.
在撰写本文时,此实现不包括数据透视,因此它在数值上并不稳定,但在大多数情况下应能很好地工作.
根据数组在内存中的排列方式,使用Cython可能仍会更快.
这是两个执行相同功能的Cython函数,但它们首先沿M进行迭代.
它不是向量化的,但是相对较快.
from numpy cimport ndarray as ar
cimport cython
@cython.boundscheck(False)
@cython.wraparound(False)
def lu3d(ar[double,ndim=3] A):
cdef int n, i, j, k, N=A.shape[0], h=A.shape[1], w=A.shape[2]
for n in xrange(N):
for j in xrange(h-1):
for i in xrange(j+1,h):
#change to L
A[n,i,j] /= A[n,j,j]
#change to U
for k in xrange(j+1,w):
A[n,i,k] -= A[n,i,j] * A[n,j,k]
@cython.boundscheck(False)
@cython.wraparound(False)
def lusolve(ar[double,ndim=3] A, ar[double,ndim=2] b):
cdef int n, i, j, N=A.shape[0], h=A.shape[1]
for n in xrange(N):
for j in xrange(h-1):
for i in xrange(j+1,h):
b[n,i] -= A[n,i,j] * b[n,j]
for j in xrange(h-1,-1,-1):
b[n,j] /= A[n,j,j]
for i in xrange(j):
b[n,i] -= A[n,i,j] * b[n,j]
您也可以尝试使用Numba,尽管在这种情况下,我无法让它像Cython一样快地运行.
标签:scipy,matrix,vectorization,python,numpy
来源: https://codeday.me/bug/20191123/2063993.html
python numpy逆_python-具有numpy的N * M * M张量的矢量化(部分)逆相关推荐
- python正态分布随机数_Python使用numpy产生正态分布随机数的向量或矩阵操作示例...
本文实例讲述了Python使用numpy产生正态分布随机数的向量或矩阵操作.分享给大家供大家参考,具体如下: 简单来说,正态分布(Normal distribution)又名高斯分布(Gaussian ...
- python傅里叶变换库_python的numpy库和cv2库实现图像傅里叶变换
码字不易,如果对您有所帮助,记着点赞哦! 一. 图像傅里叶变换原理: 对二维图像进行傅里叶变换用如下公式进行: 图像长M,高N.F(u,v)表示频域图像,f(x,y)表示时域图像.u的范围为[0,M- ...
- python阈值计算_python – 在numpy中计算超过阈值的数组值的最快方法
使用cython可能是一个不错的选择. import numpy as np cimport numpy as np cimport cython from cython.parallel impor ...
- python numpy逆_Python使用numpy计算矩阵特征值、特征向量与逆矩阵
原标题:Python使用numpy计算矩阵特征值.特征向量与逆矩阵 Python扩展库numpy.linalg的eig()函数可以用来计算矩阵的特征值与特征向量,而numpy.linalg.inv() ...
- python自定义随机数_python:numpy.random模块生成随机数
简介 所谓生成随机数,即按照某种概率分布,从给定的区间内随机选取一个数.常用的分布有:均匀分布(uniform distribution),正态分布(normal distribution),泊松分布 ...
- python创建矩阵_python中Numpy的属性与创建矩阵
本篇文章给大家带来的内容是关于python中Numpy的属性与创建矩阵,有一定的参考价值,有需要的朋友可以参考一下,希望对你有所帮助. ndarray.ndim:维度 ndarray.shape:形状 ...
- python之numpy基础_Python之NumPy学习(基础篇)
NumPy(NumericalPython的缩写)是一个开源的Python科学计算库.使用NumPy,就可以很自然的使用数组和矩阵.NumPy包含很多实用的数学函数,涵盖线性代数运算.傅里叶变换和随机 ...
- python numpy数据类型_Python之numpy数组学习(一)
原标题:Python之numpy数组学习(一) 我回来了. 前言 前面已经安装并学习了Python中的科学计算库,今天主要学习下numpy数组. Numpy数组对象 Numpy中的多维数组称为ndar ...
- python条件替换_Python中Numpy条件替换操作一例
为了数据分析快捷方便,实际操作中,我们往往要对字符串标签进行0和1的转换操作,如性别:男和女.还有根据条件进行转换,比如:大于60的归为1,60以下的归为2. 以下是在Numpy中进行转换的例子: & ...
最新文章
- Java基础-内部类
- Couchbase 101:从Java应用程序创建视图(MapReduce)
- [BZOJ 4300]绝世好题
- 软件工程(可行性研究讲解)
- 信号与系统 matlab实验报告,信号与系统 MATLAB实验报告
- ******:突破空格的限制
- matlab模拟光栅,matlab对光栅的仿真代码
- [转载] Python(析构函数)
- [数据结构]P1.3 栈 Stack
- System Center 2016组件将发生什么变化?
- POJ 3617 Best Cow Line
- 《云计算核心技术剖析》迷你书连载一 – 首席的推荐和前言
- 基于MATLAB实现四阶龙格库塔法求解一、二阶微分方程实例
- 机房维护 网拷_利用网络还原系统(远志)快速维护机房
- matlab求统计量:均值/中位数/极值/方差和标准差
- epcs1s是epcs1系列的么_fpga的EPCS 配置的2种方法 FPGA下载程序的方法(EPCS)
- Scala之特质特质Trait
- 金士顿服务器内存条怎么看型号,Win10怎么查看内存条型号?
- MT4软件IOS版如何下载
- 移动指数加权平均笔记
热门文章
- [实变函数]5.1 Riemann 积分的局限性, Lebesgue 积分简介
- 处理音视频合并的简单方法
- rsync复制软件应用于实践
- linux命令deploy_Linux 命令每天必学(34)之du命令
- mysql varchar类型实例_Mysql实例MySQL数据类型varchar详解
- FreeSurfer和FSL的安装和使用教程
- 【转载】Ubuntu顶部的任务栏-标题栏-菜单栏-启动器消失不见7个解决办法
- Au 入门系列之四:降噪与修复
- 老大易主,汽车互联网行业格局变了
- Mac打包Android的apk,【ReactNative】Mac下分分钟打包 Android apk