Multilinear Extension Polynomial: Compute the coefficients of the expanded polynomial Let $\textbf{a}=\left[10, 32, 57, 81\right]$ and $x=\{0, 1, 2, 3\}$.
Then multilinear extension of vector $\textbf{a}$ is the polynomial, $f_\textbf{a}(x_2, x_1) = 10(1-x_2)(1-x_1) + 32(1-x_2)x_1 + 57x_2(1-x_1) + 81x_2x_1$, where
$x_2$ and $x_1$ is the second and first bit of $x$ respectively.
Thus,
$$ x = 0, f_\textbf{a}(0, 0) = \textbf{a}[0] = 10$$
$$ x = 1, f_\textbf{a}(0, 1) = \textbf{a}[1] = 32$$
$$ x = 2, f_\textbf{a}(1, 0) = \textbf{a}[2] = 57$$
$$ x = 3, f_\textbf{a}(1, 1) = \textbf{a}[3] = 81$$
Expanding the polynomial:
\begin{align}
    f_\textbf{a}(x_2, x_1) &amp;= 10(1-x_2)(1-x_1) + 32(1-x_2)x_1 + 57x_2(1-x_1) + 81x_2x_1 \\
     &amp;=  2 x_2 x_1 + 22 x_1 + 47 x_2 + 10
\end{align}
Is there an algorithm to directly compute the coefficients of the expanded polynomial $f_\textbf{a} = 2 x_2 x_1 + 22 x_1 + 47 x_2 + 10$ without naively expanding the polynomial from $10(1-x_2)(1-x_1) + 32(1-x_2)x_1 + 57x_2(1-x_1) + 81x_2x_1$?
This can be done in $O(n2^n)$, where $n$ is the number of variables (naive expansion would take $O(2^{2n})$ in the worst case).
This is a standard technique called &quot;sum over subsets/supersets/submasks/supermasks&quot;. A usual example is the ANF computation, which does it in $GF(2)$ so it uses the XOR operation. Here it is basically the same but since we work with integers, we have to care about the sign now.
The idea is to apply recursively the matrix
$$
M=\begin{pmatrix}
1 &amp; 0\\
-1 &amp; 1\\
\end{pmatrix}
$$
which corresponds to applying the matrix $M^{\otimes n}$ to the input.
To illustrate the idea a bit further, see what happens when we apply the transformation for one bit (variable), say $x_1$. We have
$$
f_a = 10(1-x_2)(1-x_1) + 32(1-x_2)x_1 + 57x_2(1-x_1) + 81x_2x_1.
$$
What $M$ tells us to do, is to consider coefficients for terms containing $(1-x_1)$, and subtract those from coefficients of terms which are the same except that $(1-x_1)$ is replaced by $x_1$. For the example above, it tells to subtract 10 from 32, and 57 from 81. It is easy to see that this action corresponds to expanding the term $(1-x_1)$. By repeating this procedure for all variables, we obtain the desired result.
Here is python code.
def ext(a):
    if len(a) == 1:
        return (a[0],)
    n = len(a)
    h = n // 2
    l = ext(a[:h])
    r = ext(a[h:])
    return l + tuple(vr - vl for vl, vr in zip(l, r))


a = [10, 32, 57, 81]
b = ext(a)
print(b)  # (10, 22, 47, 2)

To clarify the notation.

in the input array $a$, $a_i$ is the coefficient of the product of: $x_j$ if $j$-th bit of $i$ in the binary representation is 1 or $1-x_j$ otherwise;
in the output array $b$, $b_i$ is the coefficient of the product of $x_j$  if $j$-th bit of $i$ in the binary representation is 1 or $1$ otherwise.

Bits are counted starting from the most significant bit. E.g. for $n=3$, $i=3$ in the input corresponds to $(1-x_1)x_2x_3$, and $i=3$ in the output corresponds to $x_2x_3$.