Thanks Greg! These recommendations worked great. Dominique, thanks for pointing out your operator infrastructure! Your lazy evaluation protocol is exactly what I need...
Madeleine On Thursday, April 21, 2016 at 6:36:57 AM UTC-7, Dominique Orban wrote: > > I have a whole operator infrastructure here: > https://github.com/JuliaOptimizers/LinearOperators.jl (develop branch) > but I've never tried SVD. It should work with normest & co. as all that's > required is matvecs. > > On Wednesday, April 20, 2016 at 9:17:32 PM UTC-4, Madeleine Udell wrote: >> >> Hi, >> >> I'm trying to define my own custom operator type that will allow me to >> implement my own * and '* operations for use inside eigenvalue or singular >> value routines like eigs and svds. But I'm having trouble making this work. >> >> Here's a simple example reimplementing matrix multiplication, with >> questions sprinkled as comments in the code. >> >> ----------- >> >> import Base: *, Ac_mul_b, size >> # Ac_mul_b doesn't seem to live in base; where does it live? >> >> type MyOperator >> A::Array{Float64} >> end >> >> o = MyOperator(randn(5,4)) >> >> *(o::MyOperator, v::AbstractVector) = (println("using custom * method"); >> o.A*v) >> o*rand(4) # this works >> >> Ac_mul_b(o::MyOperator, v::AbstractVector) = (println("using custom '* >> method"); o.A'*b) >> o'*rand(5) # this doesn't work; instead, it calls ctranspose(o)*v, which >> is not what I want >> >> size(o::MyOperator) = size(o.A) >> size(o::MyOperator, i::Int) = size(o)[i] >> >> svds(o, nsv=1) >> # this doesn't work; svds requires a `zero` method for the parametrized >> type of my operator. If I know the type will always be Float64, what's the >> easiest way to tell this to svds? >> >> ----------- >> >> Questions, collected: >> 1. Ac_mul_b doesn't seem to live in base; where does it live? Or should I >> be extending a different function? This causes o'*v to call >> ctranspose(o)*v, which is not what I want. (The default ctranspose for new >> types seems to be the identity.) >> 2. svds requires a `zero` method for the parametrized type of my >> operator. If I know the type will always be Float64, what's the easiest way >> to tell this to svds? Here's the eigs/svds >> <https://github.com/JuliaLang/julia/blob/master/base/linalg/arnoldi.jl> >> code, but I haven't been able to track down enough information on >> AbstractMatrices and parametrized types to make this work. >> 3. Any other methods I should implement for my operator? >> >> Thanks! >> Madeleine >> >
