Linear algebra for multidimensional polynomial fitting

In [ ]:
import numpy as np
import itertools
# A small data set
#T_in = np.array([20,24,35,40])+273.15
#x_in = np.array([47,55,70,78,82])/100.0
#rho_in = np.array([[1047,1033,1020,1000,990],[997,983,970,950,940],[947,933,920,900,895],[987,883,870,850,845]])
#
# large data set
T_in   = np.array([   -45 ,    -40 ,    -35 ,    -30 ,    -25 ,    -20 ,    -15 ,    -10])+273.15 # Kelvin
x_in   = np.array([     5 ,     10 ,     15 ,     20 ,     25 ,     30 ,     35 ])/100.0 # mass fraction
        
rho_in = np.array([
          [1064.0,    1054.6,    1045.3,    1036.3,    1027.4,    1018.6,    1010.0],
          [1061.3,    1052.1,    1043.1,    1034.3,    1025.6,    1017.0,    1008.6],
          [1057.6,    1048.8,    1040.1,    1031.5,    1023.1,    1014.8,    1006.7],
          [1053.1,    1044.6,    1036.2,    1028.0,    1019.9,    1012.0,    1004.1],
          [1047.5,    1039.4,    1031.5,    1023.7,    1016.0,    1008.4,    1000.9],
          [1040.7,    1033.2,    1025.7,    1018.4,    1011.2,    1004.0,     997.0],
          [1032.3,    1025.3,    1018.5,    1011.7,    1005.1,     998.5,     992.0],
          [1021.5,    1015.3,    1009.2,    1003.1,     997.1,     991.2,     985.4]]) # kg/m3

Minimizing the squared error $\epsilon(\mathbf{c}) = \sqrt{\sum (\mathbf{z} - \mathbf{A} \cdot \mathbf{c})^2 }$ can be achieved by solving the system of orthogonal equations given by $\mathbf{A}^\text{T}\mathbf{A} \cdot \mathbf{c} =\mathbf{A}^\text{T}\mathbf{z}$. Using Python tools, we can leave this to the software and we only have to construct the Van der Monde matrix $\mathbf{A}$ of the independent variable and equate it with the result vector of the dependent variable.

In [ ]:
def getCoeffs1d(x,z,order):
    if (len(x)<order+1): 
        raise ValueError("You have only {0} elements and try to fit {1} coefficients, please reduce the order.".format(len(x),order+1))
    A = np.vander(x,order+1)[:,::-1]
    #Anew = np.dot(A.T,A)
    #znew = np.dot(A.T,z)
    #coeffs = np.linalg.solve(Anew, znew)
    coeffs = np.linalg.lstsq(A, z)[0]
    return coeffs

xorder = 2
x = T_in[0:4]#T_in#T_in[0:4]#
z = rho_in[0:4,0]#rho_in#rho_in[0:4,0]#
coeffs = getCoeffs1d(x,z,xorder)
print coeffs
zf = np.polynomial.polynomial.polyval(x,coeffs)
print z[0]
print zf[0]
print z[-1]
print zf[-1]

We can extend the whole procedure to 2D, given that we have a solution matrix $\mathbf{Z}$ instead of a vector. Since this potentially involves a large number of coefficients, we disregard the higher order terms in order to avoid overfitting. This is done by discarding terms with a sum of exponents higher than the largest single exponent. The pair of exponents for elements $\mathbf{x}$ and $\mathbf{y}$, $i$ and $j$, has to satisfy $ i+j \leq \max(k,l)$ with $k$ and $l$ being the highest exponents in $x$ and $y$ direction respectively. The matrix of exponent pairs $\mathbf{E}$ for $k<l$ is therefore defined as $$ \mathbf{E}_{i,j} = \begin{pmatrix} (0,0) & (0,1) & \cdots & \cdots & \cdots & (0,l) \\ (1,0) & \ddots & & & (1,l-1) & (0,0) \\ \vdots & & & \udots & \udots & \vdots \\ (k,0) & \cdots & (k,l-k) & (0,0) & \cdots & (0,0) \end{pmatrix} \text{ yielding, for example, } \begin{pmatrix} (0,0) & (0,1) & (0,2) & (0,3) & (0,4) \\ (1,0) & (1,1) & (1,2) & (1,3) & (0,0) \\ (2,0) & (2,1) & (2,2) & (0,0) & (0,0) \end{pmatrix} \text{ for $k=2$ and $l=4$. } $$ Following the matrix notation would result in a 4-dimensional functional matrix from the Cartesian product of all elements of $\mathbf{E}$ and the input vectors $\mathbf{x}$ and $\mathbf{y}$. However, having many linear algebra solvers available in 2D, it is more practical to manually reduce the dimensionality of $\mathbf{A}$ to two. Each row of $\mathbf{A}$ corresponds to an unique pair elements of the input vectors $\mathbf{x}$ and $\mathbf{y}$. Evaluating the resulting expression $x_n^i y_n^j$ for all non-zero entries of $\mathbf{E}$ fills the columns of $\mathbf{A}$ with values. Every row of $\mathbf{A}$ can therefore be mapped to an entry in the solution matrix $\mathbf{Z}$. A data set consisting of 10 entries in $x$ and 20 in $y$ requires $\mathbf{Z}$ to have 200 elements and thus $\mathbf{A}$ to have 200 rows. The example given above, $k=2$ and $l=4$, leads to 10 columns in $\mathbf{A}$. After minimizing $\epsilon(\mathbf{c})$, information from $\mathbf{E}$ can be used to convert the coefficient vector $\mathbf{c}$ to a matrix to be used with two-dimensional polynomials.

In [ ]:
def getCoeffs2dmatrix(x_in,y_in,z_in,x_order,y_order):
    
    x_order += 1
    y_order += 1
    
    #To avoid overfitting, we only use the upper left triangle of the coefficient matrix
    x_exp = range(x_order)
    y_exp = range(y_order)
    limit = max(x_order,y_order)
    
    xy_exp = []
    
    # Construct the upper left triangle of coefficients    
    for i in x_exp:
        for j in y_exp:
            if(i+j<limit): xy_exp.append((i,j))
                
    x_num = len(x_in)
    y_num = len(y_in)
                
    cols = len(xy_exp)
    eqns = x_num * y_num
    #if (eqns<cols):
    #    raise ValueError("You have only {0} equations and try to fit {1} coefficients, please reduce the order.".format(eqns,cols))   
    if (x_num<x_order):
        raise ValueError("You have only {0} x-entries and try to fit {1} x-coefficients, please reduce the x_order.".format(x_num,x_order))
    if (y_num<y_order):
        raise ValueError("You have only {0} y-entries and try to fit {1} y-coefficients, please reduce the y_order.".format(y_num,y_order))
    
    #Create functional matrix
    A = np.zeros((x_num,y_num,x_order,y_order))
    for i in range(x_num):
        for j in range(y_num):
            for (xk,yk) in xy_exp:
                A[i][j][xk][yk] = x[i]**xk * y[j]**yk
                
    raise NotImplementedError("No 4-dimensional solver implemented")
    
    
def getCoeffs2d(x_in,y_in,z_in,x_order,y_order):
    
    x_order += 1
    y_order += 1
    
    #To avoid overfitting, we only use the upper left triangle of the coefficient matrix
    x_exp = range(x_order)
    y_exp = range(y_order)
    limit = max(x_order,y_order)
    
    xy_exp = []
    
    # Construct the upper left triangle of coefficients    
    for i in x_exp:
        for j in y_exp:
            if(i+j<limit): xy_exp.append((i,j))
    
    # Construct input pairs   
    xx, yy = np.meshgrid(x_in,y_in,indexing='ij')
    xx = np.array(xx.flat)
    yy = np.array(yy.flat)
    zz = np.array(z_in.flat)
    
    x_num = len(x_in)
    y_num = len(y_in)
                
    cols = len(xy_exp)
    eqns = x_num * y_num
    #if (eqns<cols):
    #    raise ValueError("You have only {0} equations and try to fit {1} coefficients, please reduce the order.".format(eqns,cols))   
    if (x_num<x_order):
        raise ValueError("You have only {0} x-entries and try to fit {1} x-coefficients, please reduce the x_order.".format(x_num,x_order))
    if (y_num<y_order):
        raise ValueError("You have only {0} y-entries and try to fit {1} y-coefficients, please reduce the y_order.".format(y_num,y_order))
    
    # Build the functional matrix
    A = np.zeros((eqns,cols))
    for i in range(eqns): # row loop
        for j, (xj,yj) in enumerate(xy_exp): # makes columns
            A[i][j] = xx[i]**xj * yy[i]**yj
             
    coeffs = np.linalg.lstsq(A, zz)[0]
    
    #Rearrange coefficients to a matrix shape
    C = np.zeros((x_order,y_order))
    for i, (xi,yi) in enumerate(xy_exp): # makes columns
        C[xi][yi] = coeffs[i]
        
    return C

xorder = 3
yorder = 3

x = T_in
y = x_in
z = rho_in

coeffs = getCoeffs2d(x,y,z,xorder,yorder)
print coeffs

print z[0][0]
print np.polynomial.polynomial.polyval2d(x[0],y[0],coeffs)

print z[-1][-1]
print np.polynomial.polynomial.polyval2d(x[-1],y[-1],coeffs)
In [ ]:
%pylab inline
import matplotlib
matplotlib.rcParams['savefig.dpi'] = 2 * matplotlib.rcParams['savefig.dpi']
from mpl_toolkits.mplot3d import Axes3D

# Construct input pairs   
xx, yy = np.meshgrid(x,y,indexing='ij')
xx = np.array(xx.flat)
yy = np.array(yy.flat)
zz = np.array(z.flat)

zf = np.polynomial.polynomial.polyval2d(xx,yy,coeffs)

fig  = plt.figure()
ax   = plt.axes(projection='3d')
ax.scatter(xx, yy, zz, color='blue', label='Original data')
ax.scatter(xx, yy, zf, color='red',  label='Fitted data')
#ax.plot(xx[0], yy[0], zs=zz[0], label='Fitted line', zdir='z')
ax.legend()