高效迭代Kronecker积

huangapple go评论62阅读模式
英文:

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

huangapple
  • 本文由 发表于 2023年6月6日 00:13:42
  • 转载请务必保留本文链接:https://go.coder-hub.com/76408250.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定