Hi, First, excuse me if I am over-optimizing, but I am curious if there exist a way to apply a function to an ndarray over a given dimension. In case I don't make myself clear, I have an array of shape( n,2,2) where each row represents a 2 by 2 covariance matrix, and I want to perform the eigenvalue decomposition of each row. Right now I do it with list comprehensions: import numpy as np import scipy as sp C = np.arange(10*2*2).reshape(10,2,2) ed = [sp.linalg.eig(r) for r in C[:]] Is there a better way, along the lines of vectorize, of doing this? Cheers, Jorge