On Fri, Apr 25, 2014 at 3:17 PM, Peter Brady <[email protected]> wrote:
> 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
It's given to you by arpack. I've updated the script for you:
from numpy import dot
from numpy.linalg import inv, solve
from numpy.random import random
from scipy.sparse.linalg import LinearOperator, eigs
n = 5
R = random((n, n))
L = random((n, n))
def matvec(x):
Sx = solve(L, dot(R, x))
S2x = solve(L, dot(R, Sx))
return -c*Sx + d*S2x
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:"
print matvec(x)
print "Largest real part n-2 eigenvalues of A using L^-1:"
print eigs(A, k=n-2, which="LR")[0]
print "Largest real part n-2 eigenvalues of A without using L^-1:"
A = LinearOperator((n, n), matvec=matvec)
print eigs(A, k=n-2, which="LR")[0]
This prints:
A*x calculation using L^-1:
[-33.82564316 -8.22498638 7.44165407 36.72209849 -1.70563463]
A*x calculation without using L^-1:
[-33.82564316 -8.22498638 7.44165407 36.72209849 -1.70563463]
Largest real part n-2 eigenvalues of A using L^-1:
[ 473.88681348+0.j 1.82296261+0.j 1.06790730+0.j]
Largest real part n-2 eigenvalues of A without using L^-1:
[ 473.88681348+0.j 1.82296261+0.j 1.06790730+0.j]
As you can see, the last line is computed completely without the
actual knowledge of A, only using the "solve" mechanism.
So this should fix the problem for you.
Can you please try this on your matrices and report back?
Ondrej
>
> 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]>
>> 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/714de0b2-3915-4455-a72e-4c06b8367f1e%40googlegroups.com.
>
> For more options, visit https://groups.google.com/d/optout.
--
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/CADDwiVBLTUjJ3z%2B-dgKOo3dMgW7SVcTEm_VqKwNc3hRJpeQVqQ%40mail.gmail.com.
For more options, visit https://groups.google.com/d/optout.
