Article

Finding an Interpolating Polynomial Using the Vandermonde Method

Use LAPACK to solve a linear system and find an interpolating polynomial to construct new points between a series of known data points.

Overview

This article demonstrates how you can generate a continuous curve that passes through a small set of points by computing an interpolating polynomial. A polynomial is the sum of a series of terms constructed from variables, coefficients, and exponents (for example, 6x³ + 7x², where 6 and 7 and the coefficients, and x is the variable); and an interpolating polynomial fills in the gaps between the supplied variables and coeefficients.

For any number of data points, there is a unique interpolating polynomial of order (that is, the largest exponent) which is the number of data points minus one. However, for large numbers of data points, this solution can become numerically unstable.

The image below shows five known points, as white dots, and the values generated by evaluating the found interpolating polynomial, as a red line:

Diagram showing a series of dots, that indicate known data points, joined by a continuous curve, that indicates interpolated values.

The code in this article determines the polynomial coefficients using a Vandermonde matrix based on the x-components of the known points. The coefficients are the solution to Ax=b, where A is the Vandermonde matrix and b is a vector of the y-components of the known points. You'll use LAPACK to solve Ax=b. LAPACK is an acronym for Linear Algebra Package and is a standard software library for numerical linear algebra.

Generate Known Data

Create an array containing five two-element vectors that describe the known data points between which the code will interpolate.

In a real-world app, you will most likely acquire data points from an external source such a meteorological or financial data source. For this example, specify x-components that are evenly distributed between 0 and 1023, and generate random y-components:

import simd

let points: [simd_double2] = [
    simd_double2(0,    Double.random(in: -10...10)),
    simd_double2(256,  Double.random(in: -10...10)),
    simd_double2(512,  Double.random(in: -10...10)),
    simd_double2(768,  Double.random(in: -10...10)),
    simd_double2(1023, Double.random(in: -10...10)),
    ]

Create a Vandermonde Matrix

Constuct a Vandermonde matrix where the rows are defined by the elements in a source vector that are successively raised to each integer power up to the the source vector's element count, minus one. For example, in the case of a five-element source vector, x, the Vandermonde matrix is of the form:

Equation that shows the form of a 5 x 5 Vandermonde matrix. Each row corresponds to an element in the source vector, x, raised to the power of 0, 1, 2, 3, and 4.

The Vandermonde matrix used in this article derives from the x-components of the points you're interpolating. For example, given the following points:

Equation that shows a 5 row matrix containing the values 0, 10, 1, -2, 2, 3, 3, -10, 4, 16.

The resulting Vandermonde matrix contains the following values:

Equation that shows the derived Vandermonde matrix.

The following code constructs a Vandermonde matrix from the points array:

import Accelerate

let exponents = (0 ..< points.count).map {
    return Double($0)
}

let vandermonde: [[Double]] = points.map { point in
    let bases = [Double](repeating: point.x,
                         count: points.count)
    return vForce.pow(bases: bases,
                      exponents: exponents)
}

Calculate Coefficients

The coefficients for the polynomial are the solution to Ax=b, where A is the Vandermonde matrix and b is the y-components of the known points. For example, using the matrix created in Create a Vandermonde Matrix, the coefficients are the x in the following:

Equation that shows Ax=b.

Create the function solveLinearSystem(a:a_rowCount:a_columnCount:b:b_count:) to encapsulate the LAPACK routines that solve Ax=b. Note that LAPACK overwrites b with the solution vector, x:

let coefficients: [Double] = {
    var a = vandermonde.flatMap{ $0 }
    var b = points.map{ $0.y }
    
    do {
        try ViewController.solveLinearSystem(a: &a,
                                             a_rowCount: points.count,
                                             a_columnCount: points.count,
                                             b: &b,
                                             b_count: points.count)
    } catch {
        fatalError("Unable to solve linear system.")
    }
    
    vDSP.reverse(&b)
    
    return b
}()

On return, coefficients contains the polynomial coefficients.

Use LAPACK to Solve Linear System

Use the LAPACK dgels routine to perform the solve. The dgels name derives from double-precision, general-matrix, least-squares.

  1. By default, LAPACK expects matrices in column-major format. Specify transpose to support the row-major Vandermonde matrix.

  2. Define constants used by LAPACK routines.

  3. Pass -1 to the LWORK parameter of dgels_ to calculate the optimal size for the workspace array.

  4. Create the workspace array based on the workspace query result.

  5. Perform the solve by passing the workspace array size to the LWORK parameter of dgels_.

static func solveLinearSystem(a: inout [Double],
                              a_rowCount: Int, a_columnCount: Int,
                              b: inout [Double],
                              b_count: Int) throws {
    
    var info = Int32(0)
    
    // 1: Specify transpose.
    var trans = Int8("T".utf8.first!)
    
    // 2: Define constants.
    var m = __CLPK_integer(a_rowCount)
    var n = __CLPK_integer(a_columnCount)
    var lda = __CLPK_integer(a_rowCount)
    var nrhs = __CLPK_integer(1) // assumes `b` is a column matrix
    var ldb = __CLPK_integer(b_count)
    
    // 3: Workspace query.
    var workDimension = Double(0)
    var minusOne = Int32(-1)
    
    dgels_(&trans, &m, &n,
           &nrhs,
           &a, &lda,
           &b, &ldb,
           &workDimension, &minusOne,
           &info)
    
    if info != 0 {
        throw LAPACKError.internalError
    }
    
    // 4: Create workspace.
    var lwork = Int32(workDimension)
    var workspace = [Double](repeating: 0,
                             count: Int(workDimension))
    
    // 5: Solve linear system.
    dgels_(&trans, &m, &n,
           &nrhs,
           &a, &lda,
           &b, &ldb,
           &workspace, &lwork,
           &info)
    
    if info < 0 {
        throw LAPACKError.parameterHasIllegalValue(parameterIndex: abs(Int(info)))
    } else if info > 0 {
        throw LAPACKError.diagonalElementOfTriangularFactorIsZero(index: Int(info))
    }
}

public enum LAPACKError: Swift.Error {
    case internalError
    case parameterHasIllegalValue(parameterIndex: Int)
    case diagonalElementOfTriangularFactorIsZero(index: Int)
}

Evaluate the Polynomial

The vDSP evaluatePolynomial(usingCoefficients:withVariables:) function evaluates a polynomial. For example, the following code evaluates a simple polynomial that consists of three variables and three coefficients:

let coefficients: [Float] = [5, 6, 7] 
let variables: [Float] = [1, 2, 3]  

var c = [Float](repeating: .nan,
                count: variables.count)

let result = vDSP.evaluatePolynomial(usingCoefficients: coefficients,
                                     withVariables: variables)

On return, result contains [18.0, 39.0, 70.0] by performing the following:

Equation that shows the polynomial evaluation for the coefficients five, six, and seven, and the variables one, two, and three.

Note that the number of elements returned by the polynomial evaluation is the same as the number of elements in the variables array.

To create an interpolation result that contains 1024 elements, use ramp(in:count:) to create the variables:

let ramp = vDSP.ramp(withInitialValue: 0,
                     increment: Double(1),
                     count: 1024)

let polynomialResult = vDSP.evaluatePolynomial(usingCoefficients: coefficients,
                                               withVariables: ramp)

On return, polynomialResult contains 1024 elements with the indices corresponding the the x-components, and the values corresponding to the interpolated y-components:

Diagram showing a series of dots, that indicate known data points, joined by a continuous curve, that indicates interpolated values.

See Also

Linear Algebra

BLAS

Apple’s implementation of the Basic Linear Algebra Subprograms (BLAS).