# Generate faa di bruno matrix

**URL:** <https://discourse.julialang.org/t/generate-faa-di-bruno-matrix/38059>\
**Category:** Performance\
**Tags:** differentiation, function, algebra, symbolic\
**Created:** [April 23, 2020, 2:13am UTC](https://discourse.julialang.org/t/generate-faa-di-bruno-matrix/38059 "2020-04-23T02:13:21Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [April 23, 2020, 2:13am UTC](https://discourse.julialang.org/t/generate-faa-di-bruno-matrix/38059/1 "2020-04-23T02:13:21Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![tfehring](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tfehring/32/14271_2.png) [@tfehring](https://discourse.julialang.org/u/tfehring)\
**Post date:** [April 23, 2020, 3:23am UTC](https://discourse.julialang.org/t/generate-faa-di-bruno-matrix/38059/2 "2020-04-23T03:23:14Z")

</div>

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](https://github.com/JuliaCollections/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 `x`s across those matrices of polynomials.

---

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [April 23, 2020, 6:07am UTC](https://discourse.julialang.org/t/generate-faa-di-bruno-matrix/38059/3 "2020-04-23T06:07:44Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [April 23, 2020, 2:31pm UTC](https://discourse.julialang.org/t/generate-faa-di-bruno-matrix/38059/4 "2020-04-23T14:31:06Z")

</div>

there always is the `@memoze` macro [https://github.com/JuliaCollections/Memoize.jl](https://github.com/JuliaCollections/Memoize.jl)

---

<div class="post-metadata">

**Author:** ![tfehring](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tfehring/32/14271_2.png) [@tfehring](https://discourse.julialang.org/u/tfehring)\
**Post date:** [April 23, 2020, 5:30pm UTC](https://discourse.julialang.org/t/generate-faa-di-bruno-matrix/38059/5 "2020-04-23T17:30:40Z")

</div>

You could so something like the following:

```julia
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](https://docs.julialang.org/en/v1/manual/faq/index.html#What-does-the-...-operator-do?-1) 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.

---

<div class="post-metadata">

**Author:** ![tfehring](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tfehring/32/14271_2.png) [@tfehring](https://discourse.julialang.org/u/tfehring)\
**Post date:** [April 23, 2020, 5:59pm UTC](https://discourse.julialang.org/t/generate-faa-di-bruno-matrix/38059/6 "2020-04-23T17:59:24Z")

</div>

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

```julia
map(f -> f(x), A)

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 23, 2020, 6:04pm UTC](https://discourse.julialang.org/t/generate-faa-di-bruno-matrix/38059/7 "2020-04-23T18:04:30Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![tfehring](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tfehring/32/14271_2.png) [@tfehring](https://discourse.julialang.org/u/tfehring)\
**Post date:** [April 23, 2020, 11:28pm UTC](https://discourse.julialang.org/t/generate-faa-di-bruno-matrix/38059/8 "2020-04-23T23:28:33Z")

</div>

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&nbsp;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.

---

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [April 24, 2020, 12:39am UTC](https://discourse.julialang.org/t/generate-faa-di-bruno-matrix/38059/9 "2020-04-24T00:39:40Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![tfehring](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tfehring/32/14271_2.png) [@tfehring](https://discourse.julialang.org/u/tfehring)\
**Post date:** [April 24, 2020, 2:25am UTC](https://discourse.julialang.org/t/generate-faa-di-bruno-matrix/38059/10 "2020-04-24T02:25:31Z")

</div>

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)`:

```julia
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.