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.


(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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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