4

I'm trying to calculate an expression of the form K = P*C.T*S^-1 (implementation of a Kalman filter)

All the involved matrices are sparse, and I would of course like to avoid calculating the actual inverse.

I tried using 

import scipy.sparse.linalg as spln

self.K = self.P.dot(spln.spsolve(S.T, C).T)

The problem is that spsolve expects it's second argument to be a vector and not a matrix.

edit: Clarification, the problem could in Matlab be solved by K = P * (C / S), so what I'm looking for is a method similar to spsolve but which can accept a matrix as its second argument. This could of course be done by splitting C into a number of column vectors c1..cn and solving the problem for each of them and then reassembling them into a matrix, but I suspect doing that will be both cumbersome and inefficient.

edit2&3: The dimensions of the matrices will typically be around P~10⁶x10^6, S~100x100, C=100x10⁶. P diagonal and S symmetric and C will only have one element per row. It will be used for an implementation of a Kalman filter using sparse matrices, see

http://en.wikipedia.org/wiki/Kalman_filter#The_Kalman_filter

Saullo G. P. Castro
  • 56,802
  • 26
  • 179
  • 234
ajn
  • 225
  • 1
  • 11
  • Without computing the inverse, you can't compute `K`. What you *can* do without computing the inverse is computing `Kx` for some vector `x`, which would involve solving a linear system. – Sven Marnach Oct 14 '11 at 13:18
  • I disagree, in matlab the solution to my question would simply be: K = P*(S' \ C)' or equivantly K = P*(C / S) The fact that C is a matrix instead of a vector does not change the reasoning, you do along the lines of what you are saying by solving once for each column in C. My question is about the fact that spsolve restrict me to C being a vector whereas in Matlab it can also be done for matrices. Depending on the dimension of the matrices this can still be significantly more efficient then calculating the actual inverse. – ajn Oct 14 '11 at 14:01
  • Sorry, I implicitely assumed all matrices to be square. So why not simply iterate over the columns of C and solve for each column? Since you have to solve a linear system each time, the loop overhead wil be negligible. – Sven Marnach Oct 14 '11 at 14:13
  • Just saw your edit. No, I don't think iterating over the colmuns of C will be inefficient. Depending on the size of the linear system, consider using an iterative solver instead of `spsolve()`, though. – Sven Marnach Oct 14 '11 at 14:16
  • If there isn't any better method I will definitely try solving it by iterating over the columns, but I have a hunch it will not be so efficient. The dimensions of the matrices will typically be around P~10⁶x10^6, S~100x100, C=100x10⁶. P and S will be diagonal and C will only have one element per row. I will update my question with this information aswell. – ajn Oct 14 '11 at 14:32
  • With the given dimensions, this will of course be inefficient! But if S is diagonal, why not simply invert it? – Sven Marnach Oct 14 '11 at 14:46
  • Good point, unfortunately my statement was incorrect :( S is merely symmetric, not diagonal. However, S has rather small dimensions compared to the others so I'm thinking that it might be worthwhile to just invert it as a dense matrix. But I would still like an answer to the general question, someone pointed out to me that the underlying code for spsolve seems to handle the general problem and that it is just the scipy wrapper to unecessary restrict the parameters. – ajn Oct 14 '11 at 14:57
  • The general rule of thumb should be: If the number of coulmns of C is less than the dimension of S, iterate over the columns of C, otherwise invert S. And if S is only 100x100, don't worry about inverting it -- it will take a millisecond or so. – Sven Marnach Oct 14 '11 at 15:08

1 Answers1

2

As a workaround can do

import numpy as np
from scipy.sparse.linalg import splu

def spsolve2(a, b):
    a_lu = splu(a)
    out = np.empty((A.shape[1], b.shape[1]))
    for j in xrange(b.shape[1]):
        out[:,j] = a_lu.solve(b[:,j])
    return out

self.K = self.P.dot(spsolve2(S.T, C).T)

But yes, it's a bug that spsolve does not accept matrices.

However, as your S is not very large, you can as well use a dense inverse.

pv.
  • 33,875
  • 8
  • 55
  • 49