Computing moments of discrete probability distribution

Question

I am wondering whether or not there is a computationally efficient way to compute the first $N$ moments $$m_k=\sum_{n=1}^{N}p_nx_n^k,\;\;\;\;k=1,...,N$$ of a probability mass function with mass $p_1,...,p_N$ on $N$ points $x_1\leq x_2\leq...\leq x_N$, $x_i\in\mathbb{R}$. By computationally efficient, I mean $O(N\log N)$ complexity or similar. Computing the moments directly gives $O(N^2)$ complexity. Algorithms which return the moments to within a tolerance $\epsilon$, similar to fast multipole methods, would be of interest.

I am particularly interested in the case $-1\leq x_1\leq...\leq x_n\leq 1$. If the $|x_i|$ are uniformly bounded away from 1, then only the first $O(\log N)$ of the $m_k$ are significantly different from 0, since $|x_i^k|\leq \epsilon/N$ for $k\geq C\log N$. Then the interesting moments can be computed in $O(N\log N)$ time. I am wondering what can be said in the general case where $|x_i|$ may be close to 1.

Answer 1

In the either the finite or infinite field case this problem can be solved in $O(M(n)\log(n))$ base field operations, where $M(n)$ is the time it takes to multiply polynomials in degree $n$ (so using FFT, this is quasi-linear).

This is in Chapter 3 of Pan's Structured Matrices and Polynomials.

You can adapt this technique here. Since the book is not publicly available I have some screenshots below:

Answer 2

Here's a bit more direct way than with Vandermonde transpose. Consider the generating function

$$ f(x) = \sum\limits_{k=0}^\infty m_k x^k $$

It rewrites as

$$ f(x) = \sum\limits_{k=0}^\infty \sum\limits_{n=1}^N p_n (x_n x)^k = \sum\limits_{n=1}^N p_n \sum_{k=0}^\infty (x_n x)^k = \sum\limits_{n=1}^N \frac{p_n}{1-x_n x}. $$

This sum can be computed as $\frac{P(x)}{Q(x)}$ in $O(M(N) \log N)$ with divide and conquer algorithm.

Then one can compute the first $N$ coefficients of $\frac{P(x)}{Q(x)}$ as

$$ m_0 + m_1 x + \dots + m_N x^N \equiv P(x) Q(x)^{-1} \pmod{x^{N+1}}, $$

where $Q(x)^{-1}$ is the inverse polynomial modulo $x^{N+1}$.