Matrix multiplication with Numpy

Question

Hi all,

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

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)

Answer 2

Numpy allows you to exponentiate a diagonal matrix directly:

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

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

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

produces the TypeError that you have observed (exponent must be an integer).

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

Answer 3

Does numpy have square root function for matrixes? Then you could do sqrt(D) instead of (D**(1/2))

Maybe the formula should realy be written

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

Based on previous comment this formula should work in case of D being diagonal matrix (I have not chance to prove it now).

Answer 4

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