On Fri, Apr 25, 2014 at 11:54 AM, Ondřej Čertík <[email protected]> wrote: > Hi Peter, > > On Fri, Apr 25, 2014 at 11:35 AM, Peter Brady <[email protected]> wrote: >> So I start with two sparse matrices, L, and R each with data on just a few >> bands (ie 3 to 5) >> >> My goal is compute the largest and smallest eigenvalues of the matrix A >> given by: >> >> A = -c*(L^-1*R)+d*(L^-1*R)**2 where c and d are constants >> >> In my code this is written as: >> >> L = SparseMatrix(...) >> R = SparseMatrix(...) >> >> B = L.inv()*R >> >> A = np.array(-c*B+d*B**2).astype('double') >> >> I can then use scipy/ARPACK to get the values I want. If I convert L,R or B >> to numpy arrays before computing A, I get crappy eigenvalues so this has to >> be done symbolically. My problem is that while computing B is manageable >> for the matrices I'm interested (from 20x20 to 160x160), computing A takes >> about 5 minutes and eats up a 15-30% of my memory so I need to run this in >> serial. In contrast, if I convert B to a numpy array, it takes < 1s to >> compute A (although it is the wrong A, so it's essentially worthless). > > Can you post code that does the above? You can use gist: > https://gist.github.com/ > or IPython notebook or whatever you prefer. > > Essentialy doing L.inv() means you can't use Lapack, because the > condition number > will be bad. But I wonder if there is a way to rewrite the problem > without doing L.inv(), > possibly using some generalized eigenvalue problem or something, so > that you can still > use Lapack.
There is! The idea is that for arpack, you only need to know A*x, not the A directly. So here is how to calculate A*x without explicitly calculating L^-1: from numpy import dot from numpy.linalg import inv, solve from numpy.random import random n = 5 R = random((n, n)) L = random((n, n)) c, d = 1.5, 2.5 S = dot(inv(L), R) A = -c*S + d*dot(S, S) x = random(n) print "A*x calculation using L^-1:" print dot(A, x) print "A*x calculation without using L^-1:" Sx = solve(L, dot(R, x)) S2x = solve(L, dot(R, Sx)) print -c*Sx + d*S2x This prints: A*x calculation using L^-1: [-0.21153976 2.50748822 2.03177497 4.24144355 -4.15743541] A*x calculation without using L^-1: [-0.21153976 2.50748822 2.03177497 4.24144355 -4.15743541] You can see that both approaches exactly agree. So you can continue using numerics for this, which is always better, if you can. Ondrej > > It's always good to do solid numerics, it's faster, you can use > Fortran and so on. But a separate issue is how to do this > symbolically, and SymPy should have the best algorithms. I need to > play with your code > to get better opinion how to speed this up. > > > >> Is there some way to speed this up and/or reduce the memory footprint? >> Ideally, I would like to run hundreds (maybe thousands) of different cases. >> I'm fine with installing the necessary libraries on my machine (linux). > > Perfect. This might be another great application for CSymPy. We are > implementing the matrix > support this summer as part of GSoC. > > Ondrej -- You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To post to this group, send email to [email protected]. Visit this group at http://groups.google.com/group/sympy. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CADDwiVCnLwS2Ds50aK0fgbYN6V2mtPda5uvqdoRN7NBM2yipHQ%40mail.gmail.com. For more options, visit https://groups.google.com/d/optout.
