I am trying to implement the funcion A(x1, x2, x3, ..xn)
Such that [A]_{m,j} = B_{m,j} , where B_{m,j} corresponds to the incomplete bell polynomial.
The thing is, that I don’t want to have to construct that matrix every time that I want to use it, especially because constructing the polynomials can be expensive.
What would be the Julia way of doing this?
Is the expensive part constructing the polynomials or evaluating them (or both)? Which of those are you doing many times? (Again, it could be both.)
One option that comes to mind is to memoize the calculations (or a subset of them) using Memoize.jl. If it’s quick to evaluate the polynomials once they’ve been created, or if the creation of those polynomials is the only part that you’re repeating multiple times for the same input, you could create a function A(m::Integer, j::Integer)::Matrix{Function} to create a matrix of polynomials given m and j; memoize that function, or just store the matrices of polynomials as local variables if you prefer; and then broadcast the tuples of xs across those matrices of polynomials.
I’m evaluating the same polynomials several times.
So I could technically calculate symbolically the polynomials and then save them.
The issue is that I don’t know how I would do this in Julia.
You could so something like the following:
using Base.Iterators
using Memoize
@memoize function A(nrow::Integer, ncol::Integer)::Matrix{Function}
map(args -> bell_polynomial(args...), product(1:nrow, 1:ncol))
end
function bell_polynomial(m::Integer, j::Integer)
function(x::Vector{<:Real})
# Polynomial expression in x, given m and j
end
end
Here bell_polynomial() is a function that returns a function, and A() returns an nrow by ncol matrix of those functions. Note that args -> bell_polynomial(args...) is there to splat the (n, k) pairs returned by product(1:nrow, 1:ncol) to the arguments of bell_polynomial(). There’s probably a cleaner way to do this, but this does work.
Depending on your usage, it might be more performant to memoize bell_polynomial() instead of A(). It probably isn’t helpful (and might be harmful) to memoize both. But it should be pretty easy to check.
Edit: I originally collect()ed the iterator returned by product(1:nrow, 1:ncol), but it turns out that’s not needed.
One more thing I’ll add to my previous comment - once you have a matrix A of functions, if you have a vector x at which you want to evaluate all of the polynomials in that matrix, you can do that with
map(f -> f(x), A)
Why not just precompute a matrix of the polynomial coefficients, which you can then use to evaluate the polynomial via Horner’s method or whatever? Messing around with map and @memoize seems overly complicated, and probably slower.
I’m not that familiar with Bell polynomials. so my solution might be overly general and overly complicated because of that. But my impression is that all of the coefficients need to be constructed independently for each pair (m, j) - meaning that a matrix of coefficients would really need to be a matrix of vectors of coefficients, or some other effectively 3-dimensional data structure. I’d consider a matrix of polynomials simpler and cleaner than a matrix of vectors of coefficients (or the equivalent 3d array or whatever), though that claim is subjective, and you’re probably right about performance.
If there’s a recurrence or something that lets one of those dimensions be collapsed - e.g., if it’s possible to construct a single matrix of coefficients such that the j th column contains the coefficients of B_{m,j} for any (small enough) m, then I agree that what you described is the better design. It’s not immediately obvious to me that that’s possible, but I could easily be missing something.
For this case I think it would be better to save the polynomials instead of the matrices.
If I already have calculated B_{n,m} and now I want to store this functions in an array. Should I just say something like:
A = Array{Function,2}
A[n,m] = B_{n,m}
PD: I’m calculating the polynomials using a symbolic library
Assigning to all the values in an array isn’t quite that straightforward, though it’s still a one-liner. The following code from my earlier comment returns an nrow by ncol 2D array of whatever type is returned by bell_polynomial(::Integer, ::Integer):
map(args -> bell_polynomial(args...), product(1:nrow, 1:ncol))
In my example, bell_polynomial(::Integer, ::Integer) returned a function, so the result of that code was an Array{Function,2}. Replacing bell_polynomial with B will result in a 2D array of whatever type is returned by B(::Integer, ::Integer) - probably either a function or a custom type defined in the symbolic library you’re using.