In a former life, before joining Enthought, I was a researcher in operator algebras, a field of mathematics which can be briefly summed up as the study of certain sets of “infinite-dimensional matrices.” My first reaction when seeing numpy (actually, numerical python, this was a while back) was “A matrix module! In Python! Show me more!”

But then I looked a little deeper, and my second reaction was “Why on earth would you not have the product of two arrays be the matrix product?” And then: “And why is function application element-wise?”

From the point of view of an operator algebraist everything interesting comes out of non-commutative multiplication, and clearly the “right” way to apply a function to a matrix is

*f*(*A*) = *U* *f*(*D*) *U* *

where *A* = *U* *D* *U* * is a unitary diagonalization of *A*, and *f* is applied just to entries on the diagonal of *D*. The fact that this only works for unitarily diagonalizable matrices was a minor concern for me at the time… after all I was a pure mathematician!

Since then, particularly in a more recent life as a scientist at a data visualization company, I’ve come to appreciate the way that numpy slices and dices arrays of numbers. And numpy now meets me half-way with a first-class matrix type, so I can at least write A*B and have it mean what I want it to.

But for those who still think that numpy made the wrong choice for applying functions to arrays…

from numpy import * def mfunc(func): def new_func(x): x = mat(x) if not allclose(x*x.H, x.H*x): # matrix is not normal-ish raise ValueError, "matrix must be normal" # matrix is normal-ish e, U = linalg.eig(x) return U*mat(diag(func(e)))*U.H return new_func

Which you can then use as:

mexp = mfunc(exp)

for any ufunc, or as a decorator on your own functions:

@mfunc def f(x): return cos(x) + 1.0j*sin(x)

When I showed this to Travis Oliphant, he pointed out that scipy already has a function linalg.funm which does something similar to the above, but not as a decorator. So if you have scipy, you can instead use:

import scipy.linalg def mfunc(func): def new_func(x): return scipy.linalg.funm(x, func) return new_func

OndrejGreat, thanks, I was just asking exactly this question recently on the numpy list.