Constructing a sum of many functions (as a new function)

Hello all. I have question about function construction. Specifically, suppose I have constructed a Vector containing a functions. That is, my vector would look like v =[f1,f2,f3,f4,...]. I now want to produce a new function, which is the sum of the functions in that vector. So f_{tot} = f1 + f2 + f3 + .... As a toy example, see the following.

fun_vec = Vector(undef,5)
fun_vec[1] = (t->0.0)
for i = 2:5
     fun_vec[i] = (t -> t + rand())
end

This code block produces a vector of functions. If you extract each one, one at a time, they evaluate. My issue is that I need to create a function which is the sum of these functions, and I need it to be a function, not an evaluation of that function.

My end goal is that I’m going to have N functions where I need to evaluate the sum of those functions at M points. I could in principle evaluate each function on the vector of M points f_i(pts) and then sum the resulting vectors \sum f_i(pts). However, for reasons beyond this discussion, it is more convenient to first evaluate the sum function f_{tot}, then evaluate that on the vector of points \left( \sum f_i \right) (pts). I’ve tried a number of different things but none work. For example, I’ve tried defining a seed function h(t) = 0 and then iterating the step h(t) = h(t) + f_i(t). This doesn’t work (with a number of different syntax attempts). I’ve tried other things, but I’m missing something.

Any ideas for this simple task?

If I understand what you want, some options:


julia> funcs = [sin, cos]
2-element Vector{Function}:
 sin (generic function with 14 methods)
 cos (generic function with 14 methods)

julia> mapreduce(f -> f(1.0), +, funcs)
1.3817732906760363

julia> function sum_funcs(funcs,x)
           s = zero(eltype(x))
           for f in funcs
               s += f(x)
           end
           return s
       end
sum_funcs (generic function with 1 method)

julia> sum_funcs(funcs, 1.0)
1.3817732906760363
1 Like

Note that this is an abstractly typed container. A tuple of functions may be more efficient, depending on how you construct it.

1 Like

Thanks for the response. I’ve found similar solutions to this, but it is not quite what I’m looking for. Both of those examples essentially construct the function and evaluate it in one pass. What I would like to do is create a function, that I can later pass arguments to. So I would like to be able to construct f_{tot} = \sum f_i(t) as a function itself. Then be able to evaluate that at a later date. I can always fall back on this construct + evaluate strategy, but I would prefer to construct and then evaluate. The side benefit of this approach is that the function is guaranteed to be compiled only once even though it will be evaluated many times later, which may (or may not) come with some performance. I of course realize that for many evaluations, I could do some version of

mapreduce(f -> f.(pts), +, funcs)

from the example above (not sure if that exact syntax is right). But I am going to reuse this sum of functions in different places, thus the desire to have the function itself.

You then probably need a generated function or macro there (but not because of the semantics, but to avoid looping over the function array at all). Semantically you can write:

juila> funcs = (sin, cos, exp)

julia> sum_of_funcs(x,funcs=funcs) = mapreduce(f -> f(x), +, funcs)
sum_of_funcs (generic function with 2 methods)

julia> sum_of_funcs.([1,2])
2-element Vector{Float64}:
 4.100055119135082
 7.882206689209189
1 Like

If you prefer to use Arrays then consider using FunctionWrappers.jl - this will ensure type stability.

1 Like

What you need is a higher-order function: a function that returns a function.

function make_fs_sum(fs)
    function fs_sum(x)
        mapreduce(f -> f(x), +, fs)
    end
end
julia> fs = (sin, cos, exp);

julia> fs_sum = make_fs_sum(fs)
fs_sum (generic function with 1 method)

julia> fs_sum(1)
4.100055119135082

julia> fs_sum.([1, 2])
2-element Vector{Float64}:
 4.100055119135082
 7.882206689209189

Alternatively, if you’re only doing this once, you could use a let block to create a closure:

fs_sum = let
    fs = (sin, cos, exp)
    x -> mapreduce(f -> f(x), +, fs)
end

…although the performance is probably about the same as @lmiq’s sum_of_funcs function with the default argument.

1 Like

Yeah, with a few more than 3 functions that starts to do dynamic dispatch. If the need is performance critical I guess the only solution is a macro there.

One alternative is to define your own function type, and define + on that type:

julia> import Base:+

julia> struct Fun{F <: Function}
           f::F
       end

julia> (f::Fun)(x) = f.f(x)

julia> +(f1::Fun, f2::Fun) = Fun(x -> let f1 = f1, f2 = f2
                                          f1(x) + f2(x)
                                      end)

julia> f1 = Fun(sin); f2 = Fun(cos);

julia> f3 = f1 + f2

julia> f3(1)
1.3817732906760363

This approach requires more work at the start, but is ultimately more flexible (you can define any operation on your function type), and (in my experience) it is quite fast, since the compiler knows all types. A quick benchmark:

julia> using BenchmarkTools

julia> x = randn(100_000);

julia> y = zeros(100_000);

julia> @btime $y .= cos.($x) .+ sin.($x);
  1.391 ms (0 allocations: 0 bytes)

julia> @btime $y .= $f3.($x);
  1.318 ms (0 allocations: 0 bytes)
6 Likes

This right here shows the ease of Julia: just doing math and treating functions as points in a vector space by defining a + op. All in a couple of lines : )

This approach is cool! I have a question here: For defining the sum of two Funs, why not just

 +(f1::Fun, f2::Fun) = Fun(x -> f1(x)+f2(x))

It seems this produces the same results. Is there a benefit to use a let block here?

Thank you!

First and foremost: this is called type piracy, as you to neither own the function type nor the +method (both are from base)
It is a good idea to not fall into the habit of using type piracy since it can lead to very complicated bugs in other packages (which can also affected by your custom methods)
Edit: what you can do, however is to define your own plus method. For instance you can use a different Unicode symbol which you define as the plus in function space

According to Style Guide · The Julia Language, type piracy is about “extending or redefining methods in Base or other packages on types that you have not defined.”

Here we have a new type Fun. And according to the Style Guide, wrapping things in another type is a solution, so I think the above code should be fine.

1 Like

Oops sorry you’re right, i thought the dispatch was on the function type

2 Likes

You are right; there is no benefit in this case. I tend to add let blocks automatically when defining functions like this, because closures can be slow in some cases; see Performance Tips · The Julia Language.

Thank you. Indeed the performance of closures is an issue. I should also have the habit of using let.