Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.

share|cite|improve this question
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.


$$\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}, 
    Table[Transpose[lmat[n, k, 1]].dmat[k, vec], {k, 0, n - 1}]] .
    Reverse[MapIndexed[lmat[n, First[#2] - 1, #1] &, Most[vec]]]]]

Here's a test:

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

vaninv[vec].LinearAlgebra`VandermondeMatrix[vec] ==
LinearAlgebra`VandermondeMatrix[vec].vaninv[vec] ==
// FullSimplify

should yield True.

share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.