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