使用Numpy进行矩阵乘法

Question

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

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

现在，我是这样做的：L = D**(-1/2) * A * D**(1/2)。这样做对吗？

谢谢。

Answer 1

我看到的唯一问题是，如果你在用 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 2

请注意，建议使用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)))

Answer 3

Numpy 让你可以直接对一个包含正数的对角“矩阵”进行指数运算，也就是把每个元素都进行幂运算：

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

在这种情况下，结果是你预期的，因为 NumPy 实际上是对 NumPy 数组中的每个元素单独进行幂运算。

不过，它确实不允许你直接对任何 NumPy 矩阵进行幂运算：

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

这会产生你看到的类型错误（错误信息说明指数必须是一个整数——即使是那些可以用正系数对角化的矩阵）。

所以，只要你的矩阵 D 是对角矩阵，并且你的指数是正数，你就可以直接使用你的公式。