Efficient Evaluation of a Matrix of Function Products

Hi all,
I’m trying to compute every product of two functions of a set of different functions fv (e.g., a set of wavefunctions) at different x to obtain a matrix y for each x of shape NxN where N is the number of functions. Here, fv is just a vector of arbitrary functions.

The best I could come up with is the following:

using BenchmarkTools
import Base.*

*(f::Function, g::Function) = @. (x...) -> f(x...) * g(x...)

function get_function_matrix(fv::Vector{Function})
    N = length(fv)
    Mf = Array{Function, 2}(undef, N, N)
    for j in 1:N, i in 1:N
        Mf[i, j] = fv[i] * fv[j]
    return Mf

function evaluate_function_matrix(Mf::Array{Function, 2}, x)
    N = size(Mf, 1)
    y = zeros(length(x), N, N)
    for j in 1:N, i in 1:N
        y_ = @view y[:, i, j]
        map!(Mf[i, j], y_, x)    # <- Why does cause an allocation?
    return y

x = rand(100)
fv = [sin, cos, sinh, cosh, tan, tanh]
Mf = get_function_matrix(fv)
@btime y = evaluate_function_matrix(Mf, x)

My questions are:

1.) Why does the map! function cause an allocation and is there a way to optimize this? Am I wrong to expect that this step should not cause any allocations at all?
2.) If one replaces the map! function with y[i, j] = M[i, j](x), where x` is just a float, the allocations increase even more which I don’t understand.
3.) Is there a more Julian way to solve this problem?

Thank you very much in advance for your help!

It is because Functionis an abstract type and every function has their separate types. Allocation occurs because it causes boxing and unboxing if some container has element type as an abstract type.

So if my container has elements of concrete types no allocation occurs?

Yes, if they have the same concrete type. However, two distinct functions in julia have different concrete types.

If I am not mistaken, you are complicating your life unnecessarily.

Do not build the array of functions. Build a matrix of function values
V = [fv[j](x[i]) for i in 1:length(x), j in 1:N]

Wouldn’t your y[i,j,k] be V[i,j] * V[i,k]?

As others have suggested, you really should be doing this specific example a different way. However, in cases where one really does require the ability to do what you’re attempting in the original post, one can look into the FunctionWrappers.jl package.

In case you’re not aware, this is called type piracy, and is strongly discouraged.


Thanks! Is there any documentation or explanation how to use this package?

@DNF I was aware that it is far from good practice what I’m doing but I couldn’t yet find a different way to get the desired functionality.

@pitsianis For this simple example that is true. However, in practice I don’t know any of the entries of my array of functions explicitly. They are constructed on the go depending on my system as sums/products of some base functions. Of course, one could still evaluate all components at the desired points and sum/multiply them in the required way but I feel like this would be even more complicated whereas constructing arrays of functions is just really convenient.