使用 Numpy 进行矩阵乘法

Question

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

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

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

谢谢。

Answer 1

请注意，建议使用 numpy 的 array 而不是 matrix：参见 用户指南中的此段落。一些响应中的混乱是可能出错的一个例子......特别是，D**0.5 和乘积如果应用于 numpy 数组，则元素级，这会给你一个错误的答案。例如：

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 div），那么 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 是对角矩阵的情况下工作（我没有机会证明它现在）。