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