I want to invert a 10X10 antisymmetric matrix in Python around 10,000 - 20,000 times. Is there a faster way to do it other than to use the built-in inverse function in Python?
Thanks.
$\begingroup$
$\endgroup$
5
-
1$\begingroup$ Why do you want to want to perform an involutive operation more than once? $\endgroup$– Nikolaj-KJul 13, 2012 at 13:34
-
$\begingroup$ I am confused why the fact that $(A^{-1})^{-1}=A$ for any matrix $A$ doesn't make this a lot easier... $\endgroup$– Benjamin HorowitzJul 13, 2012 at 13:49
-
2$\begingroup$ Presumably, the poster is inverting 10000 matrices once, not one matrix 10000 times. $\endgroup$– user1504Jul 13, 2012 at 14:00
-
$\begingroup$ This should be on scicomp.SE $\endgroup$– Colin KJul 13, 2012 at 14:05
-
$\begingroup$ Are you inverting the same matrix repeatedly or inverting multiple $10\times10$ matrices? $\endgroup$– PaulJul 13, 2012 at 18:26
Add a comment
|
3 Answers
$\begingroup$
$\endgroup$
I picked this trick up from Jack Poulson when he answered this related question on antisymmetric (or skew-symmetric) matrix exponentials.
An antisymmetric (more commonly called skew-symmetric matrix) $A$ is one in which $A^{T} = -A$. Since the matrix wasn't called skew-Hermitian, I'm assuming that the matrix is real.
Conveniently, $(iA)^{H} = -iA^{H}$, where $H$ denotes the Hermitian transpose, so you could compute $iA$, and invert it using the LAPACK routine ZHESV (or CHESV; unless it is also positive definite, in which case you could use ZPOSV or CPOSV). At this point, you have $(iA)^{-1} = -iA^{-1} = B$. It follows that $iB = A^{-1}$.
Unfortunately, NumPy and SciPy don't implement those functions (you'd have to call them from another language, like Fortran, C, C++, Java, etc.; there could be other libraries that provide a Python interface to LAPACK, but I don't know any that implement all of it). Based on the module scipy.linalg, your best option is probably to call scipy.linalg.solve (and ignore any symmetry), or if $iA$ is Hermitian positive definite, call scipy.linalg.solveh_banded after rearranging your data appropriately.
Based on looking at the source, at a high level, it doesn't matter whether you call scipy.linalg.solve or scipy.linalg.inv; both are ultimately LAPACK calls. If any conceivable speed difference matters, you may as well test the two, but before you do so, you're probably better off making sure you use a high-quality BLAS implementation (ATLAS, Intel MKL (if appropriate), GotoBLAS) and making sure you build LAPACK, NumPy, and SciPy accordingly. Also, if there's any task parallelism with these matrices, you could potentially exploit that as well.
All of this assumes that you want to invert 10,000 - 20,000 different matrices. I presume that if you wanted to invert the same matrix that many times, you know to just calculate the LU decomposition once, and use it to solve a linear system with 10,000 - 20,000 different right hand sides (each with the same coefficient matrix), in which case, the appropriate functions are marked in scipy.linalg. I can add further details if need be.
answered Jul 14, 2012 at 10:43
$\begingroup$
$\endgroup$
1
Yes, this is the wrong place to post the question, but can't resist answering.
Python doesn't have a built-in matrix inverse. Numpy does. Numpy's algorithm is written in a low-level language, and written by matrix-inversion experts, so it's about as fast as possible. Unless you are a matrix-inversion expert yourself, you cannot write one that is faster.
-
$\begingroup$ Yep, I meant I was using numpy's inverse function. Thanks for the help though! $\endgroup$– tut_einsteinJul 13, 2012 at 15:00
$\begingroup$
$\endgroup$
You should really be asking in the Maths or IT Stack Exchanges, however if you want a layman's view:
Python's matrix inversion is likely to be written in a high level language and compiled so it will be a lot faster than interpreted Python code
as far as I know there aren't any cunning short-cut algorithms specifically for anti-symmetric matrices
So I would use the Python implementation rather than writing your own (in Python).
