search

JuliaLinearAlgebra / LinearMaps.jl

A Julia package for defining and working with linear maps, also known as linear transformations or linear operators acting on vectors. The only requirement for a LinearMap is that it can act on a vector (by multiplication) efficiently.
Other
303 stars 42 forks source link

Type issue with `kron` #185

Closed JeffFessler closed 2 years ago

JeffFessler commented 2 years ago

It seems that kron may not be inheriting the eltype of component linear maps:

using FFTW: fft
using LinearMaps
M,N = 64, 32 # image size
T = ComplexF64
Ax = LinearMap(T, x -> fft(x), M, M)
Ay = LinearMap(T, x -> fft(x), N, N)
Ak = kron(Ay, Ax) # has wrong eltype (non-complex)
Ak * rand(T, M*N) # fails with MethodError, perhaps due to type issue?

I was working on a suggestion for #183 when I hit this issue.

dkarrasch commented 2 years ago

The eltype parameter needs to be in curly brackets, otherwise it will interpret T as the "forward" map, and fft as the adjoint. 😄 When no eltype is given, the default Float64 is picked.

using FFTW: fft
using LinearMaps
M,N = 64, 32 # image size
T = ComplexF64
Ax = LinearMap{T}(fft, M, M)
Ay = LinearMap{T}(fft, N, N)
Ak = kron(Ay, Ax) # has wrong eltype (non-complex)
Ak * rand(T, M*N) # fails with MethodError, perhaps due to type issue?
JeffFessler commented 2 years ago

Oh gosh, I am so used to seeing T as the first (optional) type argument in the Julia ecosystem, e.g., ones(T, ...) and rand(T, ...), that I forgot the {T} syntax for LM. Doh!

But I am surprised that A = LinearMap(Float64, cumsum, 8, 8) even works and returns a FunctionMap. I'm proposing #186 so that next time anyone makes this mistake they will get an error at construct time, instead of a useless FunctionMap.

I realize that a type actually is a kind of function, cf, Float32(3), even though Float32 isa Function is false. Hard to imagine a use case for a type as the function in a FunctionMap though...