Implementing a convolution operator

[Jilhg]

Given a (periodic) function $f$ I would like to implement its convolution operator $C_f$, i.e. an operator that I can later apply on a function to give $C_f g = \int dy f(x-y)g(y)$.

After doing some search through the code, I don't think a convolution operator has been implemented in ApproxFun, though as a complete newbie, tbh, I am quite lost in the code. I could imagine a convolution was implemented as least implicitly, because I guess it is used for multiplying two functions over any Fourier space.

I managed to implement and explicit convolution for two functions on a Laurent() space. This is convenient, as the convolution basically just amounts to multiplying the coefficients of both functions:

#assume that f and g are Laurent over the same domain
function convolution_approxfun(f,g)
    x0=first(ApproxFunBase.domain(f))
    x1=last(ApproxFunBase.domain(f))
    L=x1-x0
    #fourier index of the n'th coefficient
    fourier_index = [if mod(a,2)==1 div(a-1,2) else -div(a,2) end for a=1:ncoefficients(f)]
    return Fun(Laurent(ApproxFunBase.domain(f)),f.coefficients.*g.coefficients.*L.*exp.(-2*pi*1im/L.*fourier_index*x0))
end

It works, but it is of course not very flexible code (Note the extra phase has to be added in case the interval starts at $x_0 \neq 0$).

The reason I need to implement the convolution operator is that I would like to solve some linear equation that contains a convolution operator. So the above explicit implementation is not enough for me, I need to implement an operator that I can add and multiply to other operators.

From what I understand about the functionality of ApproxFun.jl I think a good implementation would look something like this:

1. Using a conversion operator to convert a function on any periodic domain to the canonical domain with the canonical space (I think its Fourier()) on it.
2. Implementing a convolution operator on this canonical space, which is just multiplying the coefficients
3. Using a conversion operator to convert the function back to the original domain/space

I have the feeling all of these steps could be very easy to implement, but I am lacking the detailed understanding of how ApproxFun.jl works in the background. I would be great if someone could help me with that: How do I implement a custom operator in general? How do I take care of the conversion properly?

PS: In case others agree I would be great to add convolution operators to the general functionality of ApproxFun in the future. At least for me convolution operators appear frequently in scientific computing.

[MikaelSlevinsky]

This code is quite old, but in case it's helpful to you https://github.com/JuliaApproximation/SingularIntegralEquations.jl/blob/master/src/Operators/Convolution.jl

[Jilhg]

Thank you, that's indeed what I was looking for.

I tried to use the package, but already the ConvolutionHelmholtz.jl example did not work. Is the SingularIntegralEquations repository discontinued?

I managed to include the Convolution.jl file into my own project. However, I had to change spaces slightly, for instance rangespace(C::ConcreteConvolution{Laurent{D},BT}) where {D,BT} = space(C.G) to rangespace(C::ConcreteConvolution{Laurent{D,R},BT}) where {D,R,BT} = space(C.G). Apparently at least Laurent{PeriodicSegment{Float64},Float64} is not identified by Laurent{D} where D.

Btw fyi, after testing the code I realized the following: On a domain $(0,b)$ the Convolution operator of a Laurent function works, but if I represent the function using Fourier I got a wrong result! (I checked the code and somehow it only uses the odd coefficients of the function, i.e. only the cosine part is involved). Also, if use a domain $(a,b)$ the result for a Laurent function is shifted depending on $a$ (I think this is because the code misses the extra phase).

Anyway, I think I will be able to implement what I need using Convolution.jl as a template. So thank you very much.

Just for curiosity: Is there a way of mapping a function on a domain to a canonical domain? What I have in mind is to implement the convolution operator on the canonical domain and then use this to generalize to other domains.

[MikaelSlevinsky]

Hi @Jilhg, sorry for the slow reply! It isn't being actively maintained, but it wouldn't take much to get the code to work again. That code was quite experimental, so it's indeed possible that not all mapped domains are covered and all Fourier/Laurent variants.

I think the ApproxFun functions tocanonical and fromcanonical functions implement this mapping for you.

[Jilhg]

Hi @MikaelSlevinsky, no problem, you have been a great help already. Thank you again! I checked the two functions tocanonical and fromcanonical and from my understanding they return the identity function on the domain/space (maybe I am using them incorrectly).
But anyway, I don't need them anymore because I simply implemented the operator on any PeriodicSegment by hand (which is probably also much more efficient).

Btw, in case you decide to revive SingularIntegralEquations.jl, I would be happy to contribute my code: I defined the convolution operator for PeriodicSegment of any real interval for both Laurent and Fourier spaces; however, only for the usual line integral, not for the singular integral (not sure what that even does).

In case you don't want to revive SingularIntegralEquations.jl, but you or the community thinks having convolution operators implemented in ApproxFun would be a useful addition to the project in general, I could also try to make a pull request to implement a convolution operator in one of the maintained repositories (I guess ApproxFunFourier.jl would be natural candidate here). I have the code and it’s short and it works well, so I am quite confident that I could do that. However, just to warn you, it would be my first time contributing to such a big project. If you prefer not to include convolution operators or have doubts for any reason, I fully respect that. I just want to offer :)

What do you think?

[knightshrub]

Hi @Jilhg, did you post your code somewhere?

I'm currently in the need of a convolution operator myself and have tried to combine the code from @MikaelSlevinsky and your comments into some working code, unfortunately the type signatures are throwing me for a loop.

I've managed to implement a getindex(C::ConcreteConvolution{Fourier{D,R},T}, k::Integer, j::Integer) where {D <: PeriodicSegment, R, T} method with the correct behavior but unfortunately when I try to apply it to a function defined in domainspace(C) I get the following error

ConcreteConvolution : Fourier(【-1.0,1.0❫) → Fourier(【-1.0,1.0❫)
Error showing value of type ConcreteConvolution{ApproxFunBase.SumSpace{Tuple{CosSpace{PeriodicSegment{Float64}, Float64}, SinSpace{PeriodicSegment{Float64}, Float64}}, PeriodicSegment{Float64}, Float64}, Int64}:
ERROR: Override getindex(::ConcreteConvolution{ApproxFunBase.SumSpace{Tuple{CosSpace{PeriodicSegment{Float64}, Float64}, SinSpace{PeriodicSegment{Float64}, Float64}}, PeriodicSegment{Float64}, Float64}, Int64}, ::Integer, ::Integer)
Stacktrace:
  [1] defaultgetindex(B::ConcreteConvolution{ApproxFunBase.SumSpace{…}, Int64}, k::Int64, j::Int64)
    @ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:361
  [2] getindex(B::ConcreteConvolution{ApproxFunBase.SumSpace{…}, Int64}, k::Int64, j::Int64)
    @ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:350
  [3] getindex(V::ApproxFunBase.SubOperator{Int64, ConcreteConvolution{…}, Tuple{…}, Tuple{…}, Tuple{…}}, k::Int64, j::Int64)
    @ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/SubOperator.jl:313
  [4] default_BandedMatrix(S::ApproxFunBase.SubOperator{Int64, ConcreteConvolution{…}, Tuple{…}, Tuple{…}, Tuple{…}})
    @ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:970
  [5] BandedMatrix
    @ ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:878 [inlined]
  [6] AbstractArray
    @ ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:958 [inlined]
  [7] defaultgetindex(::ConcreteConvolution{…}, ::Val{…}, ::UnitRange{…}, ::UnitRange{…})
    @ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:397
  [8] defaultgetindex
    @ ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:378 [inlined]
  [9] getindex(B::ConcreteConvolution{ApproxFunBase.SumSpace{…}, Int64}, k::UnitRange{Int64}, j::UnitRange{Int64})
    @ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/Operators/Operator.jl:350
 [10] show(io::IOContext{Base.TTY}, mimetype::MIME{Symbol("text/plain")}, B::Operator; header::Bool)
    @ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/show.jl:99
 [11] show(io::IOContext{Base.TTY}, mimetype::MIME{Symbol("text/plain")}, B::Operator)
    @ ApproxFunBase ~/.julia/packages/ApproxFunBase/L5iEf/src/show.jl:89

I'm assuming this might be because I don't quite understand why Fourier gets expanded into ApproxFunBase.SumSpace{Tuple{CosSpace{PeriodicSegment{Float64}, Float64}, SinSpace{PeriodicSegment{Float64}, Float64}}, PeriodicSegment{Float64}, Float64} under the hood.

[Jilhg]

I have used the following preliminary code so far. I checked that it works for both Laurent and Fourier spaces. I am planning on cleaning it up a bit and creating a PR in the next couple of weeks. Any feedback is welcome :)

export Convolution

abstract type Convolution{S,T} <: Operator{T} end

struct ConcreteConvolution{S<:Space,T} <: Convolution{S,T}
    G::Fun{S,T}
end

Convolution(G::Fun{S,T}) where {S,T} = ConcreteConvolution(G)

ApproxFunBase.domain(C::ConcreteConvolution) = ApproxFunBase.domain(C.G)
ApproxFunBase.domainspace(C::ConcreteConvolution) = ApproxFunBase.space(C.G)
ApproxFunBase.rangespace(C::ConcreteConvolution) = ApproxFunBase.UnsetSpace()
ApproxFunBase.bandwidths(C::ConcreteConvolution) = error("No range space attached to Convolution")
ApproxFunBase.getindex(C::ConcreteConvolution,k::Integer,j::Integer) = error("No range space attached to Convolution")

ApproxFunBase.rangespace(C::ConcreteConvolution{Laurent{D,R}}) where {D,R} = ApproxFunBase.space(C.G)
ApproxFunBase.rangespace(C::ConcreteConvolution{Fourier{D,R}}) where {D,R} = ApproxFunBase.space(C.G)

ApproxFunBase.bandwidths(C::ConcreteConvolution{Laurent{D,R}}) where {D,R} = 0,0
ApproxFunBase.bandwidths(C::ConcreteConvolution{Fourier{D,R}}) where {D,R} = 1,1

function ApproxFunBase.getindex(C::ConcreteConvolution{Laurent{D,R},T},k::Integer,j::Integer) where {D<:PeriodicSegment,R,T}
    fourier_index::Integer = if isodd(k) div(k-1,2) else -div(k,2) end
    if k == j && k ≤ ncoefficients(C.G)
        return (exp(-2*pi*1im/arclength(domain(C.G))*fourier_index*first(domain(C.G)))*arclength(domain(C.G))*C.G.coefficients[k])::T
    else
        return zero(T)
    end
end

function ApproxFunBase.getindex(C::ConcreteConvolution{Fourier{D,R},T},k::Integer,j::Integer) where {D<:PeriodicSegment,R,T}
    fourier_index::Integer = if isodd(k) div(k-1,2) else div(k,2) end
    if k == 1 && k ≤ ncoefficients(C.G)
        if j==k
            return (arclength(domain(C.G))*C.G.coefficients[1])::T
        else
            return zero(T)
        end
    elseif 2*fourier_index ≤ ncoefficients(C.G)
        Gs::T = if 2*fourier_index ≤ ncoefficients(C.G) C.G.coefficients[2*fourier_index] else 0.0::T end # sine coefficient
        Gc::T = if 2*fourier_index+1 ≤ ncoefficients(C.G) C.G.coefficients[2*fourier_index+1] else 0.0::T end # cosine coefficient
        phase::T = 2*pi/arclength(domain(C.G))*fourier_index*first(domain(C.G))
        if iseven(k) && j==k
            return (arclength(domain(C.G))*(Gc*cos(phase)-Gs*sin(phase))/2)::T
        elseif iseven(k) && j==k+1
            return (arclength(domain(C.G))*(Gc*sin(phase)+Gs*cos(phase))/2)::T
        elseif isodd(k) && j==k
            return (arclength(domain(C.G))*(Gc*cos(phase)-Gs*sin(phase))/2)::T
        elseif isodd(k) && j==k-1
            return (arclength(domain(C.G))*(-Gc*sin(phase)-Gs*cos(phase))/2)::T
        else
            return zero(T)
        end
    else
        return zero(T)
    end
end

Just copy this code into any file MyConvolution.jl and then put include("MyConvolution.jl") in your code. The Convolution operator is then constructed for instance as

T = Convolution(Fun(sin,Fourier()))

which you can apply to a function normally.

[knightshrub]

Hi @Jilhg, thank you for the swift reply!

comparing to your code, it turns out I forgot to import ApproxFunBase.getindex in my module!
Julia didn't find the implementation, prompting me to override getindex (which I thought I did). I'll try out your code too!