Matrix multiplication with Numpy

Question

Assume that I have an affinity matrix A and a diagonal matrix D. How can I compute the Laplacian matrix in Python with nympy?

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

Currently, I use L = D**(-1/2) * A * D**(1/2). Is this a right way?

Thank you.

Answer 1

Please note that it is recommended to use numpy's array instead of matrix: see this paragraph in the user guide. The confusion in some of the responses is an example of what can go wrong... In particular, D**0.5 and the products are elementwise if applied to numpy arrays, which would give you a wrong answer. For example:

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]]

In your case, the matrix is diagonal, and so the square root of the matrix is just another diagonal matrix with the square root of the diagonal elements. Using numpy arrays, the equation becomes

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 allows you to exponentiate a diagonal "matrix" with positive elements and a positive exponent directly:

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

The result is what you expect in this case because NumPy actually applies the exponentiation to each element of the NumPy array individually.

However, it indeed does not allow you to exponentiate any NumPy matrix directly:

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

produces the TypeError that you have observed (the exception says that the exponent must be an integer–even for matrices that can be diagonalized with positive coefficients).

So, as long as your matrix D is diagonal and your exponent is positive, you should be able to directly use your formula.

Answer 3

Well, the only problem I see is that if you are using Python 2.6.x (without from __future__ import division), then 1/2 will be interpreted as 0 because it will be considered integer division. You can get around this by using D**(-.5) * A * D**.5 instead. You can also force float division with 1./2 instead of 1/2.

Other than that, it looks correct to me.

Edit:

I was trying to exponentiate a numpy array, not a matrix before, which works with D**.5. You can exponentiate a matrix element-wise using numpy.power. So you would just use

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