Question

Let $(x_1 \ldots ,x_n) \in \mathbb{R}^n$ and $f_i = \Pi_{j=1, j \neq i }^n ( x_i - x_j )$

I'm trying to evaluate $(f_1, \ldots, f_n)$. A trivial algorithm runs in $\mathcal{O}(n^2)$ but given the very specific form of the problem, there's got to be something faster. Maybe I've overlooked something simple, maybe a fourier transform is in order... What are your thoughts?

Answer 1

The following may be of some help, if you haven't thought of it already. Let $V$ be the Vandermonde matrix with $(i,j)$th entry $x^{i-1}_{j}$, $i,j=1,\ldots,n$. Its inverse $W$ has $(i,1)$th entry $$ (-1)^{n} \frac{x_1 \ldots x_{i-1} x_{i+1} … x_n}{f_i}. $$ Hence, to find $f_1,\ldots,f_n$ we need to solve $V \alpha = e_1$, where $e_1=(1,0,\ldots,0)$. Thus, the question is, how fast can one solve a system of equations with a Vandermonde matrix? There's discussion on this topic in Nick Higham's book "Accuracy and stability of numerical algorithms" (see chpt 22). It appears that there are algorithms that are $O(n (log n)^2)$. However, it's noted that these may not be numerically stable or practical.

Answer 2

I heard there's a method called non-equispaced FFT (though I know very little about it), which should be able to compute the coefficients of the polynomial $p(x)=\prod(x-x_j)$ from the $x_j$ in quasi-linear time.

Using it, we can get a divide-and-conquer algorithm as follows:

1. Divide your points in two sets $A$ and $B$ of more or less equal magnitude
2. Compute with non-equispaced FFT $g_{A,i}=\prod_{j\in A} (x_i-x_k)$ for all $i\in B$ and $g_{B,i}$ for all $i\in A$.
3. solve the two smaller sub-problems for $A$ and $B$.
4. compute the $f_i$ by multiplying the results in points 2 and 3.

I am not sure about the practicality of this nFFT though.