[julia-users] Re: a'*b and svds for custom operators
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
>>
>> 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
>>
>
[julia-users] Re: a'*b and svds for custom operators
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
>
> 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
>
[julia-users] Re: a'*b and svds for custom operators
> > > 3. Any other methods I should implement for my operator? > > http://docs.julialang.org/en/release-0.4/manual/interfaces/#abstract-arrays
[julia-users] Re: a'*b and svds for custom operators
On Thursday, April 21, 2016 at 11:17:32 AM UTC+10, 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?
>
> It's Ac_mul_B (capital B)
> type MyOperator
> A::Array{Float64}
> end
>
type can subtype from AbstractMatrix
can also parameterise to use concrete Array
type MyOperator{T} <: AbstractMatrix{T}
A::Matrix{T}
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)
>
Ac_mul_B(o::MyOperator, v::AbstractVector) = (println("using custom '*
method"); Ac_mul_B(o.A, v))
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?
>
you get zero() for free when MyOperator is a subtype of AbstractMatrix
>
> ---
>
> 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
>
> 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
>
import Base: *, Ac_mul_B, size, show
type MyOperator{T} <: AbstractMatrix{T}
A::Matrix{T}
end
*(o::MyOperator, v::AbstractVector) = (println("using custom * method");
o.A*v)
Ac_mul_B(o::MyOperator, v::AbstractVector) = (println("using custom '*
method"); Ac_mul_B(o.A, v))
size(o::MyOperator) = size(o.A)
o = MyOperator(randn(5,4))
o * rand(4)
o' * rand(5)
svds(o, nsv=1)
