Ask Your Question

jdbates's profile - activity

2013-06-14 06:08:44 +0100 received badge  Editor (source)
2013-06-14 06:07:28 +0100 asked a question Tracking eigenvectors of a 1-parameter family of matrices

My problem is this: I'm attempting a spectral decomposition of a random process via a (truncated) Karhunen-Loeve transform, but my covariance matrix is actually a 1-parameter family of matrices, and I need a way to estimate / visualize how my random process depends on this parameter. To do this, I need a way to track the eigenvectors produced by numpy.linalg.eigh().

To give you an idea of my issue, here's a sample toy problem: Suppose I have a set of points {xs}, and a random process R with covariance C(x,y) = 1/(1+a*(x-y)^2) which depends on the parameter a. For a random sample of grid points in the range [0,1] and a given choice of a (say, a=1), I can populate my covariance matrix and implement the Karhunen-Loeve transform using:

num_x = 1000
xs = numpy.array([numpy.random.uniform() for i in range(num_x)])
a0 = 1
def cov(x, y, a=a0):
    return 1/(1+a*(x-y)^2)
cov_m = numpy.array(map(lambda y: map(lambda x: cov(x,y), xs), xs))
w,v = numpy.linalg.eigh(cov_m)
z = numpy.random.standard_normal(num_x)
R = numpy.dot(v, z * w^0.5)

This will give me realizations of R with values defined at each of my random grid points xs. What I need to be able to do, however, is - for a specific realization (which means a specific choice of my grid xs and my random coefficients z) - track how R varies with respect to the parameter a in my covariance function.

This would be easy to do if I could compute the covariance matrix symbolically and specify a after the fact. However, for large matrices this is not a plausible option. The alternative method is to find a way to keep track of each eigenvector returned by numpy.linalg.eigh(). Unfortunately, numpy appears to be reordering them so as to always list the smallest eigenvalue first; this means that, when I vary a, the eigenvectors get reordered unpredictably, and the dot product numpy.dot(v,z*w^0.5) is no longer assigning the same coefficient to the same eigenvector.

Is there a way around this?