Vandermonde Matrices

Alexandre-Théophile Vandermonde was an 18th century French mathematician. As was common at the time, he contributed to a wide range of mathematical topics, including algebra and combinatorics. He is best known for a special kind of matrix that was named after him: the Vandermonde matrix.

The Vandermonde matrix is a matrix with the following form:

\[V(x_1, x_2, \ldots, x_m) = \begin{bmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^{n-1} \\ 1 & x_2 & x_2^2 & \cdots & x_2^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_m & x_m^2 & \cdots & x_m^{n-1} \end{bmatrix}\]

where \(x_1, x_2, \ldots, x_m\) are real numbers. The Vandermonde matrix is square if \(m = n\), and it is invertible if the \(x_i\) are distinct.

Vandermonde matrices are used in many areas of mathematics, including polynomial interpolation, numerical analysis, and signal processing. They are also used in cryptography and coding theory.

Interpolation

Vandermonde matrices appear in the context of polynomial interpolation. Given a set of distinct points \((x_1, y_1), (x_2, y_2), \ldots, (x_m, y_m)\), there is a unique polynomial of degree at most \(m-1\) that passes through all of these points. If the coefficients of the polynomial are \(c_0, c_1, \ldots, c_{m-1}\), then the polynomial can be written as

\[p(x) = c_0 + c_1 x + c_2 x^2 + \cdots + c_{m-1} x^{m-1}.\]

Each point \((x_i, y_i)\) gives a linear equation in the coefficients \(c_0, c_1, \ldots, c_{m-1}\):

\[y_i = c_0 + c_1 x_i + c_2 x_i^2 + \cdots + c_{m-1} x_i^{m-1}.\]

The coefficients of the polynomial can be found by solving a system of linear equations, where the matrix of the system is the Vandermonde matrix defined by the abscissas of the interpolation points.

Numerical Derivatives

A simple way to approximate the derivative of a function is to use a finite difference approximation. Given a function \(f(x)\), the derivative at a point \(x\) can be approximated by

\[f'(x) \approx \frac{f(x+h) - f(x)}{h}\]

for a small value of \(h\), the step size. This is known as the forward difference formula. In this formula, the function is evaluated at two points: \(x\) and \(x+h\). We can say that these points are at multiples of \(h\) or offsets of one and zero from the point \(x\). Looking at the numerator only, the coefficient of \(f(x)\) is \(-1\) and the coefficient of \(f(x+h)\) is \(1\).

We can use other combinations of function values at different points (with different offsets) to approximate the derivative. For example, the central difference formula is

\[f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}.\]

Here, the offsets are \(+1\) and \(-1\), and the coefficients are \(1/2\) and \(-1/2\). We can generalize this idea to approximate higher-order derivatives as well. We can approximate the \(k\)th derivative using \(n\) function evaluations as follows:

\[f^{(k)}(x) \approx \frac{\sum_i^n c_i f(x+g_ih)}{h^k}.\]

This leaves us with one question: Given a set of offsets, \(g_i\), how do we compute the coefficients \(c_i\) so that we have a good approximation?

Since \(h\) is very small, we can approximate \(f(x+h)\) by a Taylor series expansion around \(x\), and plug this into our equation:

\[f^{(k)}(x) \approx \frac{\sum_{i=1}^n c_i \sum_{j=0}^\infty f^{(j)}(x)\frac{(g_ih)^j}{j!}}{h^k}.\]

for some coefficients \(c_i\) and offsets \(g_i\). After some rearranging, we find that

\[f^{(k)}(x) \approx \sum_{j=0}^\infty f^{(j)}(x) \left( \frac{h^{j-k}}{j!} \sum_{i=1}^n c_i g_i^j \right).\]

We now compare the coefficients of derivatives of \(f(x)\) on both sides. This tells us that we’ll get a good approximation when the coefficient of \(f^{(j)}(x)\) for \(k=j\) is equal to \(1\) and zero for \(k \neq j\). Since we have only \(n\) coefficients, we can only write an equation for the first \(n\) derivatives. This gives us a system of \(n\) equations in \(n\) unknowns:

\[\sum_{i=1}^n c_i g_i^j = y_j\]

for \(j = 1, 2, \ldots, n-1\) where \(y_j\) is \(k!\) if \(k=j\) and zero otherwise. This is a Vandermonde-like system defined by the offsets. Unlike interpolation, we’re multiplying by the transpose of the Vandermonde matrix. Solving this system of equations gives us the coefficients \(c_i\) that are optimal for a given set of offsets \(g_i\).

Solving Vandermonde systems

Vandermonde matrices have a simple structure. This suggests that we might be able to exploit this structure to solve systems more efficiently. And indeed we can.

Whereas the general problem of solving a system of linear equations is \(O(n^3)\), the Vandermonde system can be solved in \(O(n^2)\) time. This is true for both the direct system we found in the interpolation example and the transposed system we found in the numerical derivatives example. The algorithms were first described by Åke Björck and Victor Pereyra in 1970.

Other applications

Vandermonde matrices and Vandermonde-like matrices appear in many other applications.

In signal processing, the matrix that defines a Discrete Fourier Transform is a Vandermonde matrix.
In cryptography, the matrix that defines the encryption and decryption functions in the RSA cryptosystem is a Vandermonde matrix.
In coding theory, they are used in error-correcting codes, particularly in the construction of Reed-Solomon codes.
In control theory, Vandermonde matrices, and specifically the Vandermonde decomposition are used in system identification.
In finance, Vandermonde matrices are used in option pricing.

The VandermondeMatrix<T> class

New in version 9.0.2 of Numerics.NET is the VandermondeMatrix<T> class. This class represents a Vandermonde-like matrix defined by the second column (the values with exponent 1), and optionally the maximum exponent.

The class provides methods to solve the system of equations for one or more right-hand sides using the Björck-Pereyra algorithms. Solving a Vandermonde system is as simple as:

var values = Vector.Create<double>(1, 2, 3, 4);
var vandermonde = new VandermondeMatrix<double>(values);
var b = Vector.Create<double>(1, 2, 3, 5);
var x = vandermonde.Solve(b);

It can also compute the determinant efficiently. These are the most common applications of Vandermonde matrices. For other operations, like arithmetic, it is more efficient to convert the matrix to a dense matrix.

Conclusion

The addition of the VandermondeMatrix<T> class to Numerics.NET gives you access to an easy and efficient way to solve systems of equations of a special form that appear in countless applications.