I've created a gist with the text for my L and R matrices as well as a simple function to read them in and turn them into SparseMatrices at https://gist.github.com/pbrady/11303375
I can switch from B = L.inv()*R to B = solve(L,R) but most of the time is not spent computing B but rather A from B. (i.e in computing A=-c*B+d*B**2) @Ondrej, I don't think I know what 'x' is On Friday, April 25, 2014 2:30:24 PM UTC-6, Ondřej Čertík wrote: > > On Fri, Apr 25, 2014 at 11:54 AM, Ondřej Čertík > <[email protected]<javascript:>> > wrote: > > Hi Peter, > > > > On Fri, Apr 25, 2014 at 11:35 AM, Peter Brady > > <[email protected]<javascript:>> > 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/714de0b2-3915-4455-a72e-4c06b8367f1e%40googlegroups.com. For more options, visit https://groups.google.com/d/optout.
