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()
In [ ]: