Efficient Evaluation of a Matrix of Function Products

General Usage Performance
1

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]
    end
    return Mf
end

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?
    end
    return y
end

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!

2

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.

3

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

4

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

5

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]?

6

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.

7

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

2 Likes
8

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

9

@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.