Operator algebra design

General Usage
1

Although I can easily construct an operator algebra in Julia using closures, I wonder if this is the best/most “Julian” way of doing it. As an example, take the following code

import Base: +, -, *, /

abstract type AbstractOperator end

#Some concrete operator type:
struct MyOperator{T} <: AbstractOperator
    x::T
end
MyOperator() = MyOperator(rand())

(f::MyOperator)(x) = f.x + x

#Operations on operators:
struct OperatorOperation{O, F<:AbstractOperator, G<:AbstractOperator} <: AbstractOperator
    op::O
    f::F
    g::G
end

(h::OperatorOperation)(x) = h.op(h.f(x), h.g(x))

for op ∈ (:+, :-, :*, :/)
    @eval Base.$op(f::AbstractOperator, g::AbstractOperator) = OperatorOperation($op, f, g)
end

After which I can do

op = MyOperator()                 #Make an operator
op(1.0)                           #Evaluate for some input.

op = MyOperator() + MyOperator(1) #Make a slightly more involved operator
op(1.0)                           #Evaluate for some input.

op = let                          #An even more involved operator
    op = MyOperator()
    for i = 1:100
        ⋆ = rand((+,-,*,/))       #Randomly select binary operations, just to make the case non-trivial.
        op = op ⋆ MyOperator()
    end
    op
end
op(1.0)                           #Evaluate for some input.

All of the above works just fine, but the (massive) type signature* of the last example made me wonder if this approach is the recommended way to go in Julia… Or is there some better pattern to use**?

*Sure, I can make my own pretty printing of AbstractOperator, but that’s beside the point…
**Sure, the example above could (and should) be optimized by e.g. combining the sum of two MyOperators into one, but that’s also somewhat beside the point…

2

This could be an idiomatic approach, provided your application is type stable, ie operator types are predictable from inputs that construct them, are stored in homogeneous containers, etc.

The large type signature is not problematic per se, only if you get into some combinatorial explosion in your application.

These questions are difficult to answer without some context though.

3

Thanks. Good to know!

(In my application, I will only construct a few operators during a “modeling stage”, after which they will be used in some heavier calculations. I was a bit worried that the huge type signature might be a sign of a bad design… :sweat_smile:)

4

If you need to create just a few of these operators and evaluate them repetitively, this might be ok. but please note that you need to compile the composition every time:

julia> myops = [MyOperator() for i in 1:100];

julia> ops = rand((+,-,*,/), 100);

julia> @time op = let res = MyOperator(), ops = ops, X = myops
           for (op, x) in zip(ops, X)
              res = op(res, x)
           end
           res
       end;
 0.770064 seconds (514.76 k allocations: 29.287 MiB, 1.12% gc time)
julia> @time op = let res = MyOperator(), ops = ops, X = myops
           for (op, x) in zip(ops, X)
              res = op(res, x)
           end
           res
       end;
  0.000255 seconds (1.30 k allocations: 67.563 KiB)
julia> ops = rand((+,-,*,/), 100);

julia> myops = [MyOperator() for i in 1:100];

julia> @time op = let res = MyOperator(), ops = ops, X = myops
           for (op, x) in zip(ops, X)
              res = op(res, x)
           end
           res
       end;
  0.780709 seconds (514.70 k allocations: 29.269 MiB, 1.13% gc time)
5

I played with that idea in the context of static operator algebra:

https://github.com/mschauer/StaticLinearMaps.jl/blob/master/src/StaticLinearMaps.jl

1 Like