python - 使用 Numpy 进行矩阵乘法

Question

假设我有一个亲和矩阵 A 和一个对角矩阵 D。如何使用 nympy 在 Python 中计算拉普拉斯矩阵？

L = D^(-1/2) AD^(1/2)

目前，我使用 L = D**(-1/2) * A * D**(1/2)。这是正确的方法吗？

谢谢你。

Answer 1

请注意，建议使用 numpyarray代替matrix：请参阅用户指南中的此段落。一些响应中的混淆是可能出错的一个例子......特别是，如果应用于 numpy 数组，D**0.5 和产品是元素级的，这会给你一个错误的答案。例如：

import numpy as np
from numpy import dot, diag
D = diag([1., 2., 3.])
print D**(-0.5)
[[ 1.                 Inf         Inf]
 [        Inf  0.70710678         Inf]
 [        Inf         Inf  0.57735027]]

在您的情况下，矩阵是对角线，因此矩阵的平方根只是另一个对角线元素的平方根的对角线矩阵。使用 numpy 数组，方程变为

D = np.array([1., 2., 3.]) # note that we define D just by its diagonal elements
A = np.cov(np.random.randn(3,100)) # a random symmetric positive definite matrix
L = dot(diag(D**(-0.5)), dot(A, diag(D**0.5)))

Answer 2

Numpy 允许您直接对具有正元素和正指数的对角线“矩阵”求幂：

m = diag(range(1, 11))
print m**0.5

结果是您在这种情况下所期望的，因为 NumPy 实际上将幂运算单独应用于 NumPy 数组的每个元素。

但是，它确实不允许您直接对任何 NumPy 矩阵求幂：

m = matrix([[1, 1], [1, 2]])
print m**0.5

产生您观察到的 TypeError（例外情况是指数必须是整数——即使对于可以用正系数对角化的矩阵也是如此）。

因此，只要您的矩阵 D 是对角线并且您的指数是正数，您就应该能够直接使用您的公式。

Answer 3

好吧，我看到的唯一问题是，如果您使用的是 Python 2.6.x（不带from __future__ import division），那么 1/2 将被解释为 0，因为它将被视为整数除法。您可以改用 D**(-.5) * A * D**.5 来解决这个问题。您还可以使用 1./2 而不是 1/2 强制浮点除法。

除此之外，它对我来说看起来是正确的。

编辑：

我试图对一个 numpy 数组求幂，而不是之前的矩阵，它适用于D**.5. 您可以使用 numpy.power 对矩阵元素进行取幂。所以你只需使用

from numpy import power
power(D, -.5) * A * power(D, .5)

Answer 4

numpy 是否具有矩阵的平方根函数？然后你可以做 sqrt(D) 而不是 (D**(1/2))

也许公式真的应该写

L = (D**(-1/2)) * A * (D**(1/2)) 

根据之前的评论，这个公式应该适用于 D 是对角矩阵的情况（我现在没有机会证明它）。