不使用for循环构建对角矩阵

Question

我正在尝试在Python中构建以下矩阵，但不想使用for循环：

A 
 [[ 0.1  0.2  0.   0.   0. ]
 [ 1.   2.   3.   0.   0. ]
 [ 0.   1.   2.   3.   0. ]
 [ 0.   0.   1.   2.   3. ]
 [ 0.   0.   0.   4.   5. ]]

我试过在NumPy中使用fill_diagonal方法（见下面的矩阵B），但它生成的矩阵和矩阵A不一样：

B 
 [[ 1.   0.2  0.   0.   0. ]
 [ 0.   2.   0.   0.   0. ]
 [ 0.   0.   3.   0.   0. ]
 [ 0.   0.   0.   1.   0. ]
 [ 0.   0.   0.   4.   5. ]]

这是我用来构建这些矩阵的Python代码：

import numpy as np
import scipy.linalg as sp   # maybe use scipy to build diagonal matrix?

#---- build diagonal square array using "for" loop
m = 5

A = np.zeros((m, m))

A[0, 0] = 0.1
A[0, 1] = 0.2

for i in range(1, m-1):
    A[i, i-1] = 1   # m-1
    A[i, i] = 2     # m
    A[i, i+1] = 3   # m+1

A[m-1, m-2] = 4
A[m-1, m-1] = 5

print('A \n', A)

#---- build diagonal square array without loop
B = np.zeros((m, m))

B[0, 0] = 0.1
B[0, 1] = 0.2

np.fill_diagonal(B, [1, 2, 3])

B[m-1, m-2] = 4
B[m-1, m-1] = 5

print('B \n', B)

那么，有没有办法在不使用for循环的情况下，构建出像矩阵A那样的对角矩阵呢？

Answer 1

注意到数字0的数量总是3（或者说是一个常数，当你想要像这样创建一个对角矩阵时）：

In [10]:    
import numpy as np
A1=[0.1, 0.2]
A2=[1,2,3]
A3=[4,5]
SPC=[0,0,0] #=or use np.zeros #spacing zeros
np.hstack((A1,SPC,A2,SPC,A2,SPC,A2,SPC,A3)).reshape(5,5)
Out[10]:
array([[ 0.1,  0.2,  0. ,  0. ,  0. ],
       [ 1. ,  2. ,  3. ,  0. ,  0. ],
       [ 0. ,  1. ,  2. ,  3. ,  0. ],
       [ 0. ,  0. ,  1. ,  2. ,  3. ],
       [ 0. ,  0. ,  0. ,  4. ,  5. ]])

---

In [11]:    
import itertools #A more general way of doing it
np.hstack(list(itertools.chain(*[(item, SPC) for item in [A1, A2, A2, A2, A3]]))[:-1]).reshape(5,5)
Out[11]:
array([[ 0.1,  0.2,  0. ,  0. ,  0. ],
       [ 1. ,  2. ,  3. ,  0. ,  0. ],
       [ 0. ,  1. ,  2. ,  3. ,  0. ],
       [ 0. ,  0. ,  1. ,  2. ,  3. ],
       [ 0. ,  0. ,  0. ,  4. ,  5. ]])

Answer 2

在 scipy.sparse 里有一些函数可以用来处理这个问题，比如：

from scipy.sparse import diags

C = diags([1,2,3], [-1,0,1], shape=(5,5), dtype=float)

C = C.toarray()

C[0, 0] = 0.1
C[0, 1] = 0.2
C[-1, -2] = 4
C[-1, -1] = 5

对角矩阵通常是非常稀疏的，所以你也可以把它保持为稀疏矩阵。根据具体的应用，这样做甚至可能会带来很大的效率提升。

---

稀疏矩阵能带来的效率提升很大程度上取决于矩阵的大小。对于一个5x5的数组，你可能觉得没必要去麻烦。但是对于更大的矩阵，使用稀疏矩阵创建数组的速度可能会快很多，下面这个例子展示了一个单位矩阵的情况：

%timeit np.eye(3000)
# 100 loops, best of 3: 3.12 ms per loop

%timeit sparse.eye(3000)
# 10000 loops, best of 3: 79.5 µs per loop

不过，稀疏矩阵数据类型的真正优势在于当你需要对稀疏数组进行数学运算时：

%timeit np.eye(3000).dot(np.eye(3000))
# 1 loops, best of 3: 2.8 s per loop

%timeit sparse.eye(3000).dot(sparse.eye(3000))
# 1000 loops, best of 3: 1.11 ms per loop

或者当你需要处理一些非常大但又稀疏的数组时：

np.eye(1E6)
# ValueError: array is too big.

sparse.eye(1E6)
# <1000000x1000000 sparse matrix of type '<type 'numpy.float64'>'
# with 1000000 stored elements (1 diagonals) in DIAgonal format>