矩阵自乘N次的指数运算？

Question

我正在使用FOR循环实现矩阵的乘方：

import numpy as np
fl=2
cl=2
fl2=fl
cl2=cl
M = random.random((fl,cl))
M2 = M
Result = np.zeros((fl,cl))
Temp = np.zeros((fl,cl))
itera = 2
print('Matriz A:\n',M)
print('Matriz AxA:\n',M2)

for i in range (0,itera):
    for a in range(0,fl):
        for b in range (0,cl):  
            Result[a,b]+=M[a,b]*M[a,b]
            temp[a,b]=Result[a,b]
            Res[a,k]=M[a,b]
print('Potencia:\n',temp)
print('Matriz:\n', Result)

错误在于在 Result[a,b]+=M[a,b]*M[a,b] 中没有正确执行乘法，当我将其保存在临时矩阵中与原始矩阵相乘时，在 for i in range (0,itera): 中它不会进行下一次跳跃。

我知道我可以使用函数np.matmul来执行此操作，但我尝试使用FOR循环来完成。 示例

- Pedro

---

但是在这里，您使用了 numpy 数组进行矩阵乘法？Numpy 支持矩阵乘法，甚至支持矩阵幂。 - Willem Van Onsem

关于为什么这可能不起作用：您从未重置“Res”，也没有在乘法后将“Res”重新分配回“M”。 - Willem Van Onsem

Answer 1

你想要使用 np.linalg.matrix_power 来进行翻译。

如果你在使用numpy，不要使用for循环，使用向量化操作。

arr = np.arange(16).reshape((4,4))

np.linalg.matrix_power(arr, 3)

array([[ 1680,  1940,  2200,  2460],
       [ 4880,  5620,  6360,  7100],
       [ 8080,  9300, 10520, 11740],
       [11280, 12980, 14680, 16380]])

这与显式乘法相同：

arr @ arr @ arr

---

>>> np.array_equal(arr @ arr @ arr, np.linalg.matrix_power(arr, 3))
True

---

既然您提出了要求

如果您真的想要一种使用循环的天真解决方案，我们可以很容易地组合这些部分。首先，我们需要一种实际上可以将矩阵相乘的方法。虽然有打败n^3复杂度的选项，但本答案不会这样做。这里是一个基本的矩阵乘法函数：

def matmultiply(a, b):
  res = np.zeros(a.shape)
  size = a.shape[0]

  for i in range(size):
    for j in range(size):
      for k in range(size):
        res[i][j] += a[i][k] * b[k][j]

  return res

现在您需要一个指数函数。该函数接受一个矩阵和一个幂，将矩阵提高到该幂。

def loopy_matrix_power(a, n):
  res = np.identity(a.shape[0])
  while n > 0:
    if n % 2 == 0:
      a = matmultiply(a, a)
      n /= 2
    else:
      res = matmultiply(res, a)
      n -= 1

  return res

实际应用场景：

loopy_matrix_power(arr, 3)

array([[ 1680.,  1940.,  2200.,  2460.],
       [ 4880.,  5620.,  6360.,  7100.],
       [ 8080.,  9300., 10520., 11740.],
       [11280., 12980., 14680., 16380.]])

- user3483203

---

我知道 np.linalg.matrix_power ... 但我正在尝试使用 For 循环完成所有操作... - Pedro

@Pedro 为什么？这正是你不想用 numpy 做的相反之处。 - roganjosh

@Pedro - 如果你想手动完成，你可以使用平方取幂算法。这是一个更快的算法。 - Mr. Llama

Answer 2

这里有一些问题：

1. 在矩阵乘法完成后，你没有重置result矩阵，因此你会不断地添加更多的值；以及
2. 你从未将result赋值回m以执行下一个乘法代数。

朴素幂实现

我认为最好将矩阵乘法“封装”在一个单独的函数中，例如：

def matmul(a1, a2):
    m, ka = a1.shape
    kb, n = a2.shape
    if ka != kb:
        raise ValueError()
    res = np.zeros((m, n))
    for i in range(m):
        for j in range(n):
            d = 0.0
            for k in range(ka):
                d += a1[i,k] * a2[k,j]
            res[i, j] = d
    return res

然后，我们可以使用以下公式计算这个矩阵的幂：

m2 = m
for i in range(topow-1):
    m = matmul(m, m2)

请注意，我们不能仅使用m作为矩阵。因为如果我们写m = matmul(m, m)，那么m现在是m2。但这意味着如果我们再次执行乘法，我们得到的是m4而不是m3。

然后这将产生预期的结果：

>>> cross = np.array([[1,0,1],[0,1,0], [1,0,1]])
>>> matmul(cross, cross)
array([[2., 0., 2.],
       [0., 1., 0.],
       [2., 0., 2.]])
>>> matmul(cross, matmul(cross, cross))
array([[4., 0., 4.],
       [0., 1., 0.],
       [4., 0., 4.]])
>>> matmul(cross, matmul(cross, matmul(cross, cross)))
array([[8., 0., 8.],
       [0., 1., 0.],
       [8., 0., 8.]])

对数幂乘法

上述方法可以在O(n)（线性时间）内计算出Mⁿ，但我们可以做得更好，可以在对数时间内计算出此矩阵：我们首先检查幂是否为1，如果是，则直接返回矩阵；如果不是，则检查幂是否为偶数，如果是，则将矩阵乘以自身，并计算该矩阵的幂，但将幂除以二，因此M^{2 n}=(M×M)ⁿ。如果幂是奇数，我们也差不多这样做，只不过要用原始值M将其乘以：M^{2 n + 1}=M×(M×M)ⁿ。像这样：

def matpow(m, p):
    if p <= 0:
        raise ValueError()
    if p == 1:
        return m
    elif p % 2 == 0:  # even
        return matpow(matmul(m, m), p // 2)
    else:             # odd
        return matmul(m, matpow(matmul(m, m), p // 2))

上述内容可以更加优雅地编写，但我将其留给你练习 :).

然而需要注意的是，使用numpy数组进行标量计算通常比使用矩阵乘法（和其他函数）效率低。这些被优化过，并且是非解释性的，通常比Python等效方法表现出色。因此，我真的建议你使用它们。numpy函数也经过测试，使得其中出现bug的可能性较小。

- Willem Van Onsem