# Optimizing function to evaluate elementary symmetric polynomials

**URL:** <https://discourse.julialang.org/t/optimizing-function-to-evaluate-elementary-symmetric-polynomials/118321>\
**Category:** Numerics\
**Created:** [August 17, 2024, 9:36pm UTC](https://discourse.julialang.org/t/optimizing-function-to-evaluate-elementary-symmetric-polynomials/118321 "2024-08-17T21:36:44Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Gattu\_Mytraya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gattu_mytraya/32/45429_2.png) [@Gattu\_Mytraya](https://discourse.julialang.org/u/Gattu_Mytraya)\
**Post date:** [August 17, 2024, 9:36pm UTC](https://discourse.julialang.org/t/optimizing-function-to-evaluate-elementary-symmetric-polynomials/118321/1 "2024-08-17T21:36:44Z")

</div>

Hi,  
I am working on optimizing the following function:

```julia
@inbounds function get_symmetric_polynomials!(dest::Matrix{Complex{MultiFloat{Float64, 2}}}, roots::Matrix{Complex{MultiFloat{Float64, 2}}}, b::Int64)

    if b == 0
        dest[1, :] .= one(Complex{MultiFloat{Float64, 2}})
        return
    elseif b == 1
        dest[1, :] .= one(Complex{MultiFloat{Float64, 2}})
        dest[2, :] .= sum.(eachcol(roots))
        return
    end

    dest[1, :] .= one(Complex{MultiFloat{Float64, 2}})
    dest[2:end, :] .= zero(Complex{MultiFloat{Float64, 2}})

    for i in axes(roots, 1)
        for j in min(i, b):-1:1
            dest[j+1, :] .+= roots[i, :] .* dest[j, :]
        end
    end

    return
end

```

This function calculates elementary symmetric polynomials up to a given order b from a vector of complex numbers (roots). Specifically, given a vector X = {x\_{1}, x\_{2}, ..., x\_{N}} and an integer b, it computes e\_{a} = \sum\_{i\_{1} \< i\_{2} \< \dots \< i\_{a+1}=1}^{N}x\_{i\_{1}}x\_{i\_{2}}\dots x\_{i\_{a}} for 0 \leq a \leq b.

Any suggestions on how to optimize this function further would be greatly appreciated.

I was first using the complex numbers from DoubleFloats.jl (using ComplexF64 leads to loss of precision for the sort of roots I am interested in, but in my testing, the precision offered by DoubleFloats is sufficient for the range of `b` I am interested in, which is for `b` up to 25), but then using MultiFloats.jl gives a nice speed up.

When I profile my code, I see that about half of the time is spent in `multidimensional.jl` `_getindex`. I am not sure what that is.

Any suggestions on how to optimize this function further would be greatly appreciated.

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [August 19, 2024, 8:36pm UTC](https://discourse.julialang.org/t/optimizing-function-to-evaluate-elementary-symmetric-polynomials/118321/2 "2024-08-19T20:36:56Z")

</div>

Two quick observations:

> [@Gattu\_Mytraya](#):
>
> ```julia
> for i in axes(roots, 1)
> for j in min(i, b):-1:1
> dest[j+1, :] .+= roots[i, :] .* dest[j, :]
> end
> end
> 
> ```

This allocates for `root[i,:]` and `dest[j,:]` because slices copy in Julia. The former is also constant for the inner loop. So try instead

```julia
    @views for i in axes(roots, 1)
        tmp = roots[i, :]
        for j in min(i, b):-1:1
            dest[j+1, :] .+= temp .* dest[j, :]
        end
    end

```

You order of indices should be suboptimal since Julia `Array`s are column-major, which means you should iterate the left-most index in the inner most loop. In your code you have written a function that takes `Matrix` inputs but conceptually it is a operation on vectors. You kinda baked in a form of vectorization. I would consider this an antipattern in Julia. And note that this antipattern is responsible for the issue with the index order in this case.  
So I would recommend you to either switch the indices on your matrices, i.e. transpose them effectively, or rewrite this function to take in Vectors instead and perform a broadcast over the different instances at a higher level in the code.  
If you come from Matlab or similar: Note that in Julia you don’t need to vectorize computations for performances. A for loop has the same speed as broadcasting.

Also another comment on style: You don’t need to type annotate the function arguments. You can write this function completely generically, which makes it a lot easier to test difference number types 🙂 you are also almost there already since you use `zero` and `one`. Just drop the type annotations from the signature and use `eltype(dest)` where you hardcoded the types in the body.

---

<div class="post-metadata">

**Author:** ![Gattu\_Mytraya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gattu_mytraya/32/45429_2.png) [@Gattu\_Mytraya](https://discourse.julialang.org/u/Gattu_Mytraya)\
**Post date:** [August 19, 2024, 9:40pm UTC](https://discourse.julialang.org/t/optimizing-function-to-evaluate-elementary-symmetric-polynomials/118321/3 "2024-08-19T21:40:50Z")

</div>

Thanks for the help! I originally had a function that operated on vectors instead of matrices, and after reverting to that approach, I found it to be faster (I broadcast over the columns of `dest` and `roots`). I’ve also modified the function to work with matrices:

```julia
@inbounds function get_symmetric_polynomials!(dest::Matrix{Complex{MultiFloat{Float64, 2}}}, roots::Matrix{Complex{MultiFloat{Float64, 2}}}, b::Int64)

    dest[:, 1] .= one(Complex{MultiFloat{Float64, 2}})
    
    if b == 0
        return
    elseif b == 1
        dest[:, 2] .= sum(roots; dims=2)
        return
    end

    dest[:, 2:end] .= zero(Complex{MultiFloat{Float64, 2}})

    @views for i in axes(roots, 2)
        tmp = roots[:, i]
        for j in min(i, b):-1:1
            dest[:, j+1] .+= tmp .* dest[:, j]
        end
    end

    return
end

```

This matrix-based function performs just as well as column-wise broadcasting. I use type annotations mainly for debugging, but I will keep it in mind.