英文:
Efficient iteration over kron products
问题
我想要获得一个更紧凑和高效的算法来完成我正在做的任务。我的任务是多次将矩阵与单位矩阵进行 Kronecker 乘积。
现在我正在做的虽然有效,但无法扩展。
import scipy.sparse as sp
from scipy.sparse import csr_matrix
A=[[0., 0.],
[1., 0.]]
Id=sp.identity(2)
A_1=csr_matrix(sp.kron(A,sp.kron(Id,sp.kron(Id,Id))))
A_2=csr_matrix(sp.kron(Id,sp.kron(A,sp.kron(Id,Id))))
A_3=csr_matrix(sp.kron(Id,sp.kron(Id,sp.kron(A,Id))))
A_4=csr_matrix(sp.kron(Id,sp.kron(Id,sp.kron(Id,A))))
注意,对于每个A_i,矩阵A向右移动了一个位置。虽然这个方法有效,但是如果我需要对 A_20 进行操作,就需要写下非常长的一串 sp.kron。
这就是为什么我想要像这样的东西:
n_q=20
A_list=[]
for i in range(i,n_q):
A_i=csr_matrix(sp.kron(Id,...sp.kron(A...,sp.kron(Id,Id))))
A_list.append(A_i)
有没有人知道是否可能实现类似这样的功能?
英文:
I would like to get a more compact, and efficient, algorithm for what I am doing. My task is to take the kron product of a matrix with the identity multiple times.
Right now I am doing something that, while it works, I cannot scale it up.
import scipy.sparse as sp
from scipy.sparse import csr_matrix
A=[[0., 0.],
[1., 0.]]
Id=sp.identity(2)
A_1=csr_matrix(sp.kron(A,sp.kron(Id,sp.kron(Id,Id))))
A_2=csr_matrix(sp.kron(Id,sp.kron(A,sp.kron(Id,Id))))
A_3=csr_matrix(sp.kron(Id,sp.kron(Id,sp.kron(A,Id))))
A_4=csr_matrix(sp.kron(Id,sp.kron(Id,sp.kron(Id,A))))
Notice how for each A_i, the A matrix is moved one space to the right. While this works, at some point if I need to do it for the A_20, it would require writing down an extremely long sequence of sp.kron.
That is why I would like to get something like:
n_q=20
A_list=[]
for i in range(i,n_q):
A_i=csr_matrix(sp.kron(Id,...sp.kron(A...,sp.kron(Id,Id))))
A_list.append(A_i)
Does anybody know if something like this is possible?
答案1
得分: 1
这是一个试验性的想法。
定义一个函数来链式调用一些克朗矩阵的操作,这些矩阵存储在一个列表中:
def foo(alist):
M = alist[0]
for a in alist[1:]:
M = sparse.kron(M, a)
return M
你有两个矩阵:
A = sparse.csr_matrix([[0., 0.],
[1., 0.]])
Id = sparse.identity(2)
创建一个列表,例如,在这个示例中,有3个单位矩阵,然后是1个A矩阵,然后再有2个单位矩阵:
alist = [Id]*3 + [A] + [Id]*2
将它传递给 foo
函数:
M = foo(alist)
然后,只需遍历可能的Id/A混合列表,以不同数量的前后矩阵。
在NumPy中,提高“迭代”效率的主要方法是将其从Python移动到编译代码中,也就是使用内置的整个数组方法,或者使用编译器如numba
。
Scipy的稀疏矩阵是在NumPy的基础上工作的,大多数矩阵属性都存储为数组。有限的编译代码(Cython?)主要用于矩阵乘法(关键的稀疏计算)和不同格式之间的转换。
你有两个迭代级别,不同的 A_i
序列以及每个序列内的 kron
链接。
我对于kron
的工作原理了解不够,不能立刻提出改进建议。我猜想你可以预先计算一些链式序列,例如连续2个和4个I_d
。因此,上面的示例实际上可以写成:
sparse.kron(Id_3, sparse.kron(A, Id_2))
试试这个:
链式连接2个和3个单位矩阵:
Id2 = sparse.kron(Id, Id)
Id3 = sparse.kron(Id, Id2)
然后比较它与完全迭代的结果:
np.allclose(Id3.A, foo(alist).A)
这就是我在最初建议检查kron序列模式时所考虑的内容。
实际上,直接使用eye
矩阵也可以:
M2 = sparse.kron(sparse.eye(8), sparse.kron(A, sparse.eye(4)))
再次比较它与foo函数的结果:
np.allclose(M2.A, foo(alist).A)
这就是我在最初建议检查kron序列模式时所考虑的内容。
英文:
Here's an trial idea.
Define a function to chain some kron calls from a list of matrices:
In [50]: def foo(alist):
...: M = alist[0]
...: for a in alist[1:]:
...: M = sparse.kron(M,a)
...: return M
...:
Your two matrices:
In [51]: A=sparse.csr_matrix([[0., 0.],
...: [1., 0.]])
...: Id=sparse.identity(2)
Make a list, in this example, 3 id, and 2 after A:
In [52]: alist = [Id]*3+[A]+[Id]*2
In [53]: alist
Out[53]:
[<2x2 sparse matrix of type '<class 'numpy.float64'>'
with 2 stored elements (1 diagonals) in DIAgonal format>,
<2x2 sparse matrix of type '<class 'numpy.float64'>'
with 2 stored elements (1 diagonals) in DIAgonal format>,
<2x2 sparse matrix of type '<class 'numpy.float64'>'
with 2 stored elements (1 diagonals) in DIAgonal format>,
<2x2 sparse matrix of type '<class 'numpy.float64'>'
with 1 stored elements in Compressed Sparse Row format>,
<2x2 sparse matrix of type '<class 'numpy.float64'>'
with 2 stored elements (1 diagonals) in DIAgonal format>,
<2x2 sparse matrix of type '<class 'numpy.float64'>'
with 2 stored elements (1 diagonals) in DIAgonal format>]
And pass it to foo
:
In [54]: M = foo(alist)
In [55]: M
Out[55]:
<64x64 sparse matrix of type '<class 'numpy.float64'>'
with 512 stored elements (blocksize = 2x2) in Block Sparse Row format>
Then it's just a matter of iterating over the possible lists of Id/A mixes, with different numbers of before and after.
It saves typing if not speed.
In numpy the main way to make 'iteration' more efficient is to move it from python to compiled code, that is, by using the builtin whole-array methods. Or by using compilers like numba
.
Scipy sparse works on top of numpy, with most of the matrix attributes stored as arrays. There's limited amount of compiled code (cython?), mainly for matrix multiplication (a key sparse calculation), and conversion between formats.
You have to iteration levels, the different A_i
sequences, and kron
chaining within each.
I don't know enough of how kron
works to offhand suggest improvements on that. I suppose you could precompute some chained sequences, say for 2 and 4 I_d
in a row. So in effect the above example could be written as
sparse.kron(Id_3, sparse.kron(A, Id_2))
In fact let's try that
Chain 2 and 3 Id matrices:
In [68]: Id2 = sparse.kron(Id,Id); Id3 = sparse.kron(Id,Id2)
In [69]: Id3
Out[69]:
<8x8 sparse matrix of type '<class 'numpy.float64'>'
with 32 stored elements (blocksize = 4x4) in Block Sparse Row format>
In [70]: Id3.A
Out[70]:
array([[1., 0., 0., 0., 0., 0., 0., 0.],
[0., 1., 0., 0., 0., 0., 0., 0.],
[0., 0., 1., 0., 0., 0., 0., 0.],
[0., 0., 0., 1., 0., 0., 0., 0.],
[0., 0., 0., 0., 1., 0., 0., 0.],
[0., 0., 0., 0., 0., 1., 0., 0.],
[0., 0., 0., 0., 0., 0., 1., 0.],
[0., 0., 0., 0., 0., 0., 0., 1.]])
Actually we don't need to use kron to get that. It's just a eye
matrix for a desired size.
Now kron
3 matrices:
In [71]: M1 = sparse.kron(Id3,sparse.kron(A,Id2))
In [72]: M1
Out[72]:
<64x64 sparse matrix of type '<class 'numpy.float64'>'
with 512 stored elements in COOrdinate format>
and compare that the result of the full blown iteration:
In [73]: np.allclose(M1.A, foo(alist).A)
Out[73]: True
That's what I had in mind when first suggesting examining the kron sequence for patterns.
In fact using eye
directly:
In [92]: M2 = sparse.kron(sparse.eye(8),sparse.kron(A,sparse.eye(4)))
In [93]: M2
Out[93]:
<64x64 sparse matrix of type '<class 'numpy.float64'>'
with 32 stored elements in COOrdinate format>
In [94]: np.allclose(M2.A, foo(alist).A)
Out[94]: True
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论