LinearMaps: Transpose function not working the same as forward function

General Usage
1

I am running into an error while trying to create a LinearMap where if I put the same function in both the forward and transpose operator arguments, the forward map works but the transpose does not. My eventual goal is to create a LinearMap that does convolution (and the transpose does convolution with a flipped version of the filter), but I’ve simplified the example code to be a little simpler while still getting the error.

Here’s my MWE:

using LinearMaps
using ImageFiltering 

N = 10
x = randn(N)

Ck = LinearMap(x -> similar(x,map(n->0:n-1,size(x))),
                           x -> similar(x,map(n->0:n-1,size(x))), N, N)
# if I change the code to map(n->1:n,size(x)), then both lines below work fine 

Ck*x # this line works fine 
Ck' * x # this line gives the error below 
# note: I also get the error if doing Ck' * (Ck*x)

The error is:
BoundsError: attempt to access 10-element Array{Float64,1} at index [Base.IdentityUnitRange(0:9)]
throw_boundserror(::Array{Float64,1}, ::Tuple{Base.IdentityUnitRange{UnitRange{Int64}}}) at abstractarray.jl:538
checkbounds at abstractarray.jl:503 [inlined]
copyto!(::Array{Float64,1}, ::OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}}) at multidimensional.jl:902
At_mul_B!(::Array{Float64,1}, ::LinearMaps.FunctionMap{Float64,getfield(Main, Symbol(“##351#355”)),getfield(Main, Symbol(“##353#357”))}, ::Array{Float64,1}) at functionmap.jl:65
A_mul_B! at transpose.jl:33 [inlined]
mul!(::Array{Float64,1}, ::LinearMaps.TransposeMap{Float64,LinearMaps.FunctionMap{Float64,getfield(Main, Symbol(“##351#355”)),getfield(Main, Symbol(“##353#357”))}}, ::Array{Float64,1}, ::Float64, ::Float64) at LinearMaps.jl:26
mul!(::Array{Float64,1}, ::LinearMaps.TransposeMap{Float64,LinearMaps.FunctionMap{Float64,getfield(Main, Symbol(“##351#355”)),getfield(Main, Symbol(“##353#357”))}}, ::Array{Float64,1}) at LinearMaps.jl:24
*(::LinearMaps.TransposeMap{Float64,LinearMaps.FunctionMap{Float64,getfield(Main, Symbol(“##351#355”)),getfield(Main, Symbol(“##353#357”))}}, ::Array{Float64,1}) at LinearMaps.jl:22
top-level scope at bilevel_1d.jl:55

Any ideas on what is going on or even how to start debugging this would be much appreciated!

2

Your MWE doesn’t work for me in a fresh REPL with Julia 1.3.1

3

Apologies - I left out the ImageFiltering package. I’ve edited the code above and it works for me on a clean Julia 1.2 REPL.

I didn’t realize the ImageFiltering was needed for the particular syntax of similar to work. I am having trouble coming up with another example of this problem with LinearMaps without the OffsetArrays (which is loaded with ImageFiltering) coming into the picture, which makes me think the interaction between the two packages is at the heart of the problem (…?)

4

ok, works now.
I didn’t check for any logic but only changed 0:n-1 to 1:n

Ck = LinearMap(x -> similar(x,map(n->1:n,size(x))),
                           x -> similar(x,map(n->1:n,size(x))), N, N)

Julia indices are based 1 not 0.
Does this make sense in your case?

5

It doesn’t make sense in my case, though I noted that works as well. I’d prefer to use the OffsetArrays (which allow for the 0 or other index-based indexing) because of the way it interacts with the convolution function. There are definitely other ways to implement what I’d like to do, and I will probably have to go with one of those, but LinearMaps are simply more elegant and I would like to understand why this doesn’t work so I can use them as often as possible.

Also, I could understand there being some issues between the OffsetArray indexing and LinearMaps (there’s no documentation I can find on whether LinearMaps supports OffsetArrays or not), but I’m surprised there seems to be a difference between the forward and transpose functionality.

6

Ok, lets wait for the experts to step in…

7 
help?> LinearMap
...
LinearMap{T=Float64}(f, [fc,], M::Int, N::Int = M; kwargs...)
...
 Optionally, a corresponding function fc can be specified that
  implements the transpose/adjoint of f.
...

So maybe it should be:

Ck = LinearMap(x -> similar(x,map(n->0:n-1,size(x))),
                           x -> adjoint(similar(x,map(n->0:n-1,size(x)))) , N, N)
8

I have no idea what your function is doing? In particular, this does not seem to be a linear function. What kind of arguments do you apply this to?

9

@oheil - Unfortunately, adding in the adjoint operation won’t work (we may have imaginary data and we don’t want the complex conjugate). Interesting that that works though and I appreciate all the ideas for how to try thinking about it!

@juthohaegeman - Ultimately I am trying to use LinearMap to do a convolution, where the adjoint operation is equivalent to convolving with a flipped version of the filter. These are both linear operations, though the way we go about getting there is a bit complicated…

I know there are multiple ways to do convolution. The one we already have in other parts of the code calls imfilter!(z,xpad,(h,),NoPad(),Algorithm.FIR()) to do z = x \conv h. Because of the size of our signals, it makes sense to use the FIR option, pre-pad our input (xpad), and pre-allocate the output (z). This means that the output z has to have the correct indexes that we want the convolution function to return, which in turn leads us to the similar command I included in the original code sample. The final LinearMap code should look something like

ztemp = x -> similar(x,map(n->1:n,size(x))) # would eventually just preallocate these
_conv! = (z,xpad,h) -> imfilter!(z,xpad,(h,),NoPad(),Algorithm.FIR())
conv! = (z,x,h) -> _conv!(z,padarray(x,Pad(:circular,ntuple(_->0,ndims(x)),
        ntuple(_->size(h,1),ndims(x)) )),h)
C = (h,N) -> LinearMap(x -> conv!(ztemp(x),x,h),
                       x -> conv!(ztemp(x),x,h[end:-1:1]),  N, N)

(I did not include all the using commands here as I think debugging with the original MWE would be a lot easier, but I can add them if that would be helpful)

Based on how complicated this is becoming, I should probably think about reindexing the h variable so that z and x always have the same indexes. But, it seems simpler to integrate this with our existing variable set-up by just defining the adjoint separately:

Ck = (h,N) -> LinearMap(x -> conv!(ztemp(x),x,h), N, N)
CkT = (h,N) -> LinearMap(x -> conv!(ztemp(x),x,h[end:-1:1]), N, N)

The fact that this works (aka, both Ck*x and CkT*x execute) while C'*x gives the original error is really at the root of my question. I don’t see why the forward and adjoint operator slots for LinearMap should have different behaviors.

10

The short answer is that LinearMaps.jl is indeed not ready yet to work with generalized indices, i.e. with axes instead of size. That doesn’t sound like a change that is too difficult. I’ll take a look soon, but feel free to beat me to it. Also notifying @dkarrasch who has been the active maintainer lately !

1 Like
11

In a first quick step, I think we should remove the asymmetry showing up in the original example, i.e., where Ck*x works, but Ck'x doesn’t. Implementing the axes comparison sounds like a project for a day or two, to make sure it really works in CompositeMaps etc.

12

The reason is that we have a path that allows to simply apply the map to a vector x, but when we use the adjoint or the transpose of the map, then we first allocate a vector y, and then mul! into it. That’s when the indices of y and those of the output of your map don’t match up and the error is thrown. I think I have a fix for that, which should be available by tomorrow.

13

Thanks so much! Great package by the way!

1 Like
14

For reference: https://github.com/Jutho/LinearMaps.jl/pull/83

This solves the Ck*x vs Ck'x issue. A transition to axes will require some more work.