Inverting a Vandermonde matrix

Is there a quick way to invert the $n\times n$ Vandermonde matrix with columns $(1, x, x^2, ..., x^{n-1})$ where $x$ takes values $0,1,...,n-1$ (in ascending order from left to right)?

Perhaps by row operations $(A|I)\to (I|A^{-1})$ where $A$ is the Vandermonde matrix, and $I$ the identity matrix? I don't see how to do it though... or maybe there is another way?

Thank you.

-
See here. – Ragib Zaman Nov 19 '11 at 15:42
The "quick" ($O(n^2)$) method is to use divided differences (i.e., the Björck-Pereyra algorithm.) – J. M. Nov 19 '11 at 15:45
Somewhat related, the determinant of such a matrix, for fixed $n$, is $G(n+1)$, where $G$ the Barnes G-function – Sasha Nov 19 '11 at 15:53
Thanks, everyone! – harry Nov 19 '11 at 17:36

(Might as well...)

The Björck-Pereyra idea is based on the fact that the usual divided-difference algorithms for polynomial interpolation can be recast as matrix decompositions of the Vandermonde matrix. I'll suggest having a look at the original paper, as well as Golub and Van Loan's treatment for details; for convenience, I'll only give here the final results.

First, we set up notation. Let

$$\mathbf V(x_0,\dots,x_n)=\begin{pmatrix}1&1&\cdots&1\\x_0&x_1&\cdots&x_n\\\vdots&\vdots&\ddots&\vdots\\x_0^n&x_1^n&\cdots&x_n^n\end{pmatrix}$$

and introduce the diagonal matrix

$$\mathbf D_k=\mathrm{diag}\left(\underbrace{1,\dots,1}_{k+1},\frac1{x_{k+1}-x_0},\dots,\frac1{x_n-x_{n-k-1}}\right)$$

and the lower bidiagonal matrix

$$\mathbf L_k(h)=\left(\begin{array}{ccc|cccc} 1&&&&&&\\&\ddots&&&&&\\&&1&&&&\\\hline&&&1&&&\\&&&h&\ddots&&\\&&&&\ddots&\ddots&\\&&&&&h&1\end{array}\right)$$

where the upper left block is a $k\times k$ identity matrix.

Then,

$$\mathbf V^{-1}=\left(\mathbf L_0(1)^\top\cdot\mathbf D_0\cdots \mathbf L_{n-1}(1)^\top\cdot\mathbf D_{n-1}\right)\left(\mathbf L_{n-1}(x_{n-1})\cdots\mathbf L_1(x_1)\cdot\mathbf L_0(x_0)\right)$$

It takes less effort than usual to multiply diagonal and bidiagonal matrices together, giving Björck-Pereyra an edge over vanilla Gaussian elimination.

Here's some sundry Mathematica code:

lmat[n_Integer, k_Integer, h_] := SparseArray[{Band[{1, 1}] -> 1,
Band[{2, 1}] -> PadLeft[ConstantArray[-h, {n - k}], n]},
{n + 1, n + 1}]

dmat[k_Integer, vec_?VectorQ] :=
DiagonalMatrix[Join[ConstantArray[1, k + 1],
1/(Drop[vec, k + 1] - Drop[vec, -k - 1])]]

vaninv[vec_?VectorQ] :=
Module[{n = Length[vec] - 1},
Apply[Dot,
Table[Transpose[lmat[n, k, 1]].dmat[k, vec], {k, 0, n - 1}]] .
Apply[Dot,
Reverse[MapIndexed[lmat[n, First[#2] - 1, #1] &, Most[vec]]]]]


Here's a test:

vec = {a, b, c, d, e, f};

vaninv[vec].LinearAlgebraVandermondeMatrix[vec] ==
LinearAlgebraVandermondeMatrix[vec].vaninv[vec] ==
IdentityMatrix[Length[vec]]
// FullSimplify


should yield True.

-