A Perspective on Numpy

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
This entry was posted in Enthought Tool Suite, NumPy, Open Source, Python, SciPy on by .
avatar

About Corran Webster

Corran is Enthought's European consulting director, and has been with Enthought as a scientific software developer and Python instructor since 2008. Corran has a B.S. from the University of New South Wales and a Ph.D. in pure mathematics from UCLA and has held positions at the University of Nevada, Las Vegas and Texas A&M University, working primarily in operator algebras and functional analysis. Prior to joining Enthought, he was Chief Scientist at Compudigm International working on machine learning with self-organizing maps and large data visualization methods. Corran wrote his first Python code back in the days of Python 1.3 and has used it as his primary programming language for over 20 years.

One thought on “A Perspective on Numpy

Leave a Reply

Your email address will not be published. Required fields are marked *

Please leave these two fields as-is:

Protected by Invisible Defender. Showed 403 to 103,441 bad guys.