使用Numpy进行矩阵乘法

Question

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

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

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

谢谢。

- mrcuongnv

Answer 1

请注意，建议使用numpy的array而不是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)))

- pberkes

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是对角线的且指数为正，您就可以直接使用公式。

- Eric O Lebigot

---

谢谢。D实际上是一个对角矩阵。我现在可以做到。 - mrcuongnv

diag(range(1, 11)) 返回的是一个数组而不是矩阵，这就是原因。 - ZDL-so

好的观点。这正是pberkes在他的回答中详细阐述的内容。我更新了这个答案以便让它更清晰明了。 - Eric O. Lebigot

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)

- Justin Peel

---

嗨Justin，感谢您的回复。但是，我收到了以下错误：“TypeError：指数必须是整数”。 - mrcuongnv

你正在使用哪个版本的Numpy？这些是普通矩阵还是稀疏矩阵。对于普通的Python矩阵并将它们提升至幂方，我没有任何问题。 - Justin Peel

算了，我太蠢了，使用的是np.array而不是np.matrix。 - Justin Peel

谢谢。我需要的是真正的矩阵运算，而不是点积或类似的东西（数组乘法）。 - mrcuongnv

对于对角矩阵不是也一样吗？ - Justin Peel

Answer 4

NumPy是否有矩阵平方根函数？那么您可以使用sqrt(D)代替(D **（1/2）)

也许公式应该真正写成

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

根据之前的评论，如果D是一个对角矩阵，这个公式应该有效（我现在没有证明的机会）。

- Tony Veijalainen