Search code examples
pythoncombinatoricspolynomials

Is there a simple and efficient way to evaluate Elementary Symmetric Polynomials in Python?

I'm currently working on a project that involves evaluating Elementary Symmetric Polynomials (or "ESPs") using Python. So essentially I need to create a function that :

  • takes a non-empty list x (of numeric values) as an input
  • outputs a list containing all the corresponding ESP evaluations, where each ESP is in len(x) variables (also, by convention, we could say that the "0-th" ESP is always equal to 1)

For example, if my input list is x = [x1, x2, x3], then I want the function to output the list [1, x1 + x2 + x3, x1 * x2 + x1 * x3 + x2 * x3, x1 * x2 * x3].

I initially wrote a function that iterated over all the k-tuples of x (for all integers k in the range [1, len(x)]) using the "itertools.combinations" method, but it obviously didn't scale up well, as the time complexity of my function was exponential : it was in O(len(x) * 2^len(x)). So I was wondering if there was any algorithm that could do the same job, but in polynomial time (say, with a quadratic or a cubic time complexity).

I then stumbled upon Newton's identities (linking ESPs to power sums), and it does the job in a cubic time complexity ! For my purposes this time complexity is good enough, but is there a more efficient way of evaluating ESPs ? After some research, I didn't see any Stack Overflow posts with explicit code to solve my problem (although this Math Stack Exchange post does mention some algorithms). I also went through the following resources :

  • The pySymmPol package (and the corresponding article) : it looked promising, but in fact it doesn't support polynomial evaluation ...
  • This article, describing a fast algorithm to evaluate ESPs (with floating-point inputs) : unfortunately, it doesn't seem like the authors provide any code, and their algorithms seem very tedious to implement (also I think they said their code was written in C)

Additionally, here is the improved function that I implemented to solve my initial problem, using Newton's identities :

def evaluate_all_ESPs(x: list) -> list:
    """
    This function evaluates all the Elementary Symmetric Polynomials (or
    "ESPs") in `len(x)` variables on the given input list of numeric values,
    and returns the corresponding `len(x) + 1` evaluations in a list

    If we label these evaluated ESPs as `e[k]` (where `0 <= k <= len(x)`), we have :
        • `e[0] = 1` (by convention)
        • For all integers `1 <= k <= len(x)`, `e[k]` will be equal to the sum of all
          the products of k-tuples of `x` :
              - `e[1] = x[0] + x[1] + ... + x[-1]` (sum of all the elements of `x`)
              - `e[2] = x[0] * x[1] + x[0] * x[2] + ... + x[-2] * x[-1]` (sum of all
                the "double-products" of elements of `x`)
              - ...
              - `e[len(x)] = x[0] * x[1] * ... * x[-1]` (product of all the elements of `x`)

    Note : This function uses Newton's identities to make the computations *much*
           faster than if we had manually iterated over all the k-tuples of `x`
           (for each integer `k` in the range `[1, len(x)]`) ! Also, this algorithm's
           complexity is in `O(len(x)**3)`
    """
    # doing some quick input checks
    if not(isinstance(x, list)):
        raise TypeError(f"\n\nExpected a list, but got a '{type(x).__name__}' instead : '{x}'\n")
    if len(x) == 0:
        raise ValueError("\n\nThe input cannot be the empty list !\n")

    # initialization
    nb_elements = len(x)
    evaluations = [1] + [0] * nb_elements   # corresponds to the list labelled as `e` in the docstring
    powers_of_x = [1] * nb_elements         # list that'll contain all the `x[k]**i` values
    sums_of_powers_of_x = [0] * nb_elements # list that'll contain all the `sum(x[k]**i over k)` values

    for i in range(nb_elements):
        # updating `powers_of_x` and `sums_of_powers_of_x`
        powers_of_x = [x_k * previous_power_of_x_k for (x_k, previous_power_of_x_k) in zip(x, powers_of_x)]
        sums_of_powers_of_x[i] = sum(powers_of_x)

        # applying Newton's identity for the current evaluation
        # (SOURCE : https://en.wikipedia.org/wiki/Newton%27s_identities)
        current_evaluation = 0
        alternating_sign = 1
        for j in range(i, -1, -1):
            current_evaluation += alternating_sign * evaluations[j] * sums_of_powers_of_x[i - j]
            alternating_sign *= -1
        if isinstance(current_evaluation, int):
            current_evaluation //= i + 1
        else:
            current_evaluation /= i + 1

        # updating `evaluations`
        evaluations[i + 1] = current_evaluation

    return evaluations

Correct me if I'm wrong, but I haven't seen any explicit Python code on the internet that implements Newton's identities for the evaluation of ESPs. So hopefully this post will save some time for other people in my situation !

Solution

  • Surely you can do it just from the coefficients of the polynomial constructed from those roots?

    e.g. x1=1, x2=2, x3=3

    The polynomial is (x-1)(x-2)(x-3) = x3-6x2+11x-6

    The coefficients are the ESPs with alternating signs: 1, -(x1+x2+x3), +(x2.x3+x3.x1+x1.x2), -x1.x2.x3 etc.

    Similarly for larger lists.

    Code: (Credit @no comment for the alternating-sign correction)

    import numpy as np
from numpy.polynomial import Polynomial

def getESP( a ):
   c = np.array( Polynomial.fromroots(a).coef )
   esp = c[-1::-1]           # reverse array so that x^N is first
   esp[1::2] *= -1           # correct for alternating signs
   return esp

c = getESP( [ 1, 2, 3 ] )
print( c )

    Output:

    [ 1.  6. 11.  6.]