UP | HOME

Vandermonde matrices

Table of Contents

Solving polynomials

pages 185-186 introduces the Vandermonde matix which provides us with a way to solve for polynomials that fit a set of given points.

Polynomials are equations of the form:

\begin{equation} y=a_{1}+a_{2}x+a_{3}x^{2}+a_{4}x^{3}+... \end{equation}

It's a typicaly non-linear function of the form y=f(x) that takes inputs x and returns outputs y. Given a polynomial (ie. given its an factors) and given a value for y there can be multiple solutions for x. At first blush this doesn't seem related to the problems we've been tackling with matrices.

The Vandermonde matrix does a clever inversion of our typical problem. While we typically are looking for x given y here because we are ultimately interested in finding a polynomial that fits what we do instead is set up the problem to solve for the a_{n} polynomial factors. Instead of taking y=f(x) and solving for x we are taking a corresponding y=f(a) and solving for the polynomial factors an. We will see that the x's are baked into the f().

To see how this is built up we simply take the polynomial function and a series of points/measurements [ (x_1,y_1) (x_2,y_2) (x_3,y_3) ... ] and write out the equations.

\begin{equation} y_1=a_{1}+a_{2}x_1+a_{3}x_{1}^{2}+a_{4}x_{1}^{3}+...\\ y_2=a_{1}+a_{2}x_2+a_{3}x_{2}^{2}+a_{4}x_{2}^{3}+...\\ y_3=a_{1}+a_{2}x_3+a_{3}x_{3}^{2}+a_{4}x_{3}^{3}+...\\ ... \end{equation}

No we we want to pull out the a's to make this look like an Ax=b type of equation. While the x's in the equations look nonlinear they're actually known values coordinates of our x,y points so they're fixed scalar values and they not something we are solving for. Writen out in matrix form we get

\begin{equation} \begin{bmatrix} 1 & x_1 & x_{1}^2 & x_{1}^3 ..\\ 1 & x_2 & x_{2}^2 & x_{2}^3 ..\\ 1 & x_3 & x_{3}^2 & x_{3}^3 ..\\ ...\\ \end{bmatrix} \begin{bmatrix} a_1\\ a_2\\ a_3\\ a_4\\ ...\\ \end{bmatrix} = \begin{bmatrix} y_1\\ y_2\\ y_3\\ ...\\ \end{bmatrix} \end{equation}

Again, all the values of the x-matrix are known. (the x-matrix is the actual Vandermonde matrix - V)

Vandermonde matrices

(ns morelinear.vandermonde
  (:require [morelinear.leastsquares :as leastsquares]
            [clojure.core.matrix :refer :all]
            [clojure.core.matrix.linear :as linear]))
(set-current-implementation :vectorz)

So given a series of x's we can build up this matrix directly. 

(defn polynomial-matrix
  "Take a vector of input X's and build a vandermonde matrix.
  The DEGREE shouldn't be larger than the number of values given.
  If non provided it will be chosen to be the maximal value
  and your resulting matrix will be square"
  ([polynomial-input-values
    degree]
   (mapv (fn [x]
           (mapv (fn [power]
                   (pow x power))
                 (range degree)))
         polynomial-input-values))

  ;; Default square matrix case
  ([polynomial-input-values]
   (polynomial-matrix polynomial-input-values
                      (ecount polynomial-input-values))))

The number of polynomials terms (ie. the degree n in the last xn) is variable.

Full rank case

If we have as many polynomial terms as we have points then we are left with a square x-matrix which we can solve using Gaussian elimination (excluding singular/degenerate cases). In other words we can use the LU decomposition to solve for all the polynomial factors an in this equivalent system Va=y - where the V is the matrix of exponents of x.

Once we solve for the an's we can then stick them back into non-linear polynomial equation we started looking at in the beginning and get a solution.

\begin{equation} y=a_{1}+a_{2}x+a_{3}x^{2}+a_{4}x^{3}+... \end{equation}

It will not only hold true for all our (x,y) points but also for all other values of x we want to test - so we will have in the end fit a polynomial curve through all our points.

Least Squares

However with our least squares method we now have an alternate solution. We can choose a degree that is lower than the number of points/values we have and then solve for approximate solutions. This is often the preferrable solution. In general when we fit a polynomial we have a decent idea of the degree of the function that drives whatever we are measuring. We also know that the input values are not perfect, so if we overfit the points with a high-degree polynomial the results will also be very noisy. So typically we preselect the degree of the polynomial and then try to give the function as many data points as we can to get the best approximation of the underlying function as we can.

Note: That this doesn't always hold. If you have some non-linear function that you are approximating with a Taylor series, then cut off the number of polynomials serves no purpose

As before the Least Squares solutions come in a few flavors and it also works on the square case (even when the matrix is singular) - so we can skip writing a function for the square case and just have one polynomial solver for all square and skinny matrices 

(defn split-into-x-y
  [points]
  (let [vector-of-xs-and-ys (apply mapv vector points)]
    {:x (mutable (first vector-of-xs-and-ys))
     :y (mutable (second vector-of-xs-and-ys))}))


(defmulti polynomial-factors
  "Solve for a polynomial equation of DEGREE that fits the given POINTS
  and return the list of polynomial factors.
  the underlying least squares METHOD used can be:
  :lu
  :lu-jumbo
  :householder"
  (fn [degree
       points
       method]
    method))

(defmethod polynomial-factors :lu
  [degree
   points
   method]
  (let [xy (split-into-x-y points)]
    (leastsquares/lu-direct (polynomial-matrix (:x xy)
                                               degree)
                            (:y xy))))

(defmethod polynomial-factors :lu-jumbo
  [degree
   points
   method]
  (let [xy (split-into-x-y points)]
    (leastsquares/lu-jumbo (polynomial-matrix (:x xy)
                                              degree)
                           (:y xy))))

(defmethod polynomial-factors :householder
  [degree
   points
   method]
  (let [xy (split-into-x-y points)
        vandermonde-matrix (mutable (polynomial-matrix (:x xy)
                                                       degree))]
    (leastsquares/householder-qr vandermonde-matrix
                                 (:y xy))))

(polynomial-factors 3
                    [[138.0 165.0]
                     [275.0 197.0]
                     [277.0 262.0]
                     [186.0 277.0]
                     [154.0 240.0]
                     [236.0 141.0]
                     [322.0 123.0]
                     [358.0 122.0]
                     [182.0 289.0]
                     [118.0 294.0]
                     [108.0 292.0]
                     [54.0  259.0]]
                    :householder)
;;(251.1770462843548 0.3111353706549809 -0.001904254283169026)

Now that we've got the polynomial factors we want to just wrap this up and be able to get the polynomial function itself 

(defn polynomial-function
  [degree points method]
  (let [factors (polynomial-factors degree
                                    points
                                    method)]
    (fn [x] (reduce + (map-indexed (fn [index
                                        factor]
                                     (* factor
                                        (pow x index)))
                                   factors)))))


(let [our-function (polynomial-function 3
                                        [[138.0 165.0]
                                         [275.0 197.0]
                                         [277.0 262.0]
                                         [186.0 277.0]
                                         [154.0 240.0]
                                         [236.0 141.0]
                                         [322.0 123.0]
                                         [358.0 122.0]
                                         [182.0 289.0]
                                         [118.0 294.0]
                                         [108.0 292.0]
                                         [54.0  259.0]]
                                        :householder)]
  (println (our-function 11))
  (println (our-function 50))
  (println (our-function 70))
  (println (our-function 150))
  (println (our-function 250))
  (println (our-function 350)))
;;(251.1770462843548 0.3111353706549809 -0.001904254283169026)