[Math] Comparing a series expansion to polynomial regression

polynomialsregressiontaylor expansion

So I don't have a great background in mathematics but I have a quick and hopefully simple question for you guys.

I'm a graduate student and I'm doing some polynomial regression on some thermodynamic data to create a model and I'm kind of struggling a bit with understanding some of the differences between a series expansion and a regression. Maybe it's a subtlety or maybe I'm just ignorant.

But anyway. Let's say I take a Maclaurin series expansion to third order:

$$f(x) = f(0) + f^\prime(0) x + \frac{1}{2!}f^{\prime\prime}(0) x^2 + \frac{1}{3!}f^{\prime\prime\prime}(0)x^3 + O(x^4)$$

But this is just a polynomial:

$$f(x) = a + bx + cx^2 + dx^3$$

where the coefficients are the derivatives of $x$ and their associated coefficients.

In the series expansion, the coefficients are defined to be the $n^{th}$ derivatives computed at $0$. But when I compare this to the polynomial expansion, can I really make the identity that each of the coefficients is proportional to the derivative computed at $0$? Or are the values of the coefficients more globally defined since they have been obtained through regression?

The reason I ask is because it has important implications in the language I will use to describe the values of the coefficients in the polynomial fit.

Any help or insight would be greatly appreciated.

Edit.

I have a model where I can compute the internal energy of a system $U$ as a function of the number of particles $N$ and then some extensive property of the system $\xi$, so $U \equiv U(N,\xi)$.

I am ultimately interested in calculating the chemical potential $\mu$, which is the change in the internal energy as particles are added to they system, this is expressed as taking a derivative with respect to $N$:

$$ \mu(N,\xi) = \left(\frac{\partial U}{\partial N}\right)_{\xi}$$

So since I am able to calculate $U(N,\xi)$, I can compute $\mu(N,\xi)$ by fitting a 2D surface (polynomial) to the data and then computing its derivative with respect to $N$.

$$ U(N,\xi) = U(0,0) + U_{ ξ } \xi + U_N N + \frac{1}{2!}U_{ξξ} ξ ^2 + \frac{1}{2!}U_{NN} N^2 + U_{N ξ }N ξ + \dots$$

This can be rewritten as
$$ U(N,\xi) = U(0,0) + \left(\frac{\partial U(0,0)}{\partial \xi}\right)_{N} \xi + \mu(0,0) N + \frac{1}{2}\left(\frac{\partial^2 U(0,0)}{\partial ξ^2}\right)_N ξ ^2 + \frac{1}{2}\left(\frac{\partial \mu(0,0)}{\partial N}\right)_{\xi} N^2 + \left(\frac{\partial \mu(0,0)}{\partial \xi}\right) N ξ + \dots$$

So the derivatives in this model have physical meaning, since they represent thermodynamic variables and derivatives of thermodynamic variables.

So in this case,

$$ \mu(N, ξ ) = \mu(0,0) + \left(\frac{\partial \mu(0,0)}{\partial N}\right)_{\xi} N + \left(\frac{\partial \mu(0,0)}{\partial \xi}\right) ξ + \dots $$

The main question above was if I perform a polynomial regression on my $U(N, ξ )$ data, can I reasonably identify each of the polynomial's coefficients with the $\partial_i U$ values which are supposed be computed at $N = ξ = 0$ in the Taylor expansion, or will these coefficients correspond to something else since they are obtained by regression.

I suppose the fitted coefficients should truly correspond to the $N = ξ = 0$ case.

Best Answer

Also very late answer, but this thread did help me and might help others by adding the following explanation:

A Taylor expansion produces a best fit polynomial in the vicinity of some point $x_0$. In this regard, the fundamental difference between the polynomial obtained through regression and the one obtained through a Taylor expansion is that a regression fits the function not just at $x_0$ but over an etire domain, albeit at the expense of enough terms.

Hope this helps clarify.

Related Solutions

Multivariate Quadratic Regression, Surface Fitting and the Hessian

The basis for two dimensions would be

$$ p(x) = \left[\begin{array}{ccccc} x_1 & x_2 & x_1^2 & x_1x_2 & x_2^2 \end{array}\right]^T $$

$f(x) = p(x)^T a$

where $a = [a_0, a_1, ...]$. This would be akin to calculating

$$ f(x) = a_0 + a_1 x_1 + a_2 x_2 + a_3 x_1 x_2 + a_4 x_2^2 $$

Then with given data points

$$ \left[\begin{array}{ccccc} (x_1^1) & (x_2^1) & (x_1^1)^2 & (x_1^1)(x_2^1) & (x_2^1)^2 \\ \vdots & & & & \vdots \\ (x_1^n) & (x_2^n) & (x_1^n)^2 & (x_1^n)(x_2^n) & (x_2^n)^2 \\ \end{array}\right] \mathbf{a} = \left[\begin{array}{c} y^1 \\ \vdots \\ y^n \end{array}\right] $$

and solve for $a$.

But I need the

$$ f(x) = x^T A x $$

form, so I can easily get $\nabla f (x) = 2 A x$ and the Hessian $\nabla^2 f(x) = 2 A$.

It turns out if we start with 3 dimensions, and use $x_3 = 1$ this gives the necessary form in two dimensions,

$$ x^T A x = \left[\begin{array}{ccc} x_1 & x_2 & x_3 \end{array}\right] \left[\begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array}\right] \left[\begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array}\right] $$

$$ = \left[\begin{array}{c} x_1 a_{11} + x_2 a_{21} + x_3 a_{31} \\ x_1 a_{12} + x_2 a_{22} + x_3 a_{32} \\ x_1 a_{13} + x_2 a_{23} + x_3 a_{33} \end{array}\right]^T \left[\begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array}\right] $$

$$ = x_1 x_1 a_{11} + x_1 x_2 a_{21} + x_1 x_3 a_{31} + $$ $$ x_1 x_2 a_{12} + x_2 x_2 a_{22} + x_3 x_2 a_{32} + $$ $$ x_1 x_3 a_{13} + x_2 x_2 a_{23} + x_3 x_3 a_{33} $$

Using $x_3 = 1$

$$ x^T A x = \left[\begin{array}{ccc} x_1 & x_2 & 1 \end{array}\right] \left[\begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array}\right] \left[\begin{array}{c} x_1 \\ x_2 \\ 1 \end{array}\right] $$

$$ = a_{11} x_1^2 + a_{21} x_1 x_2 + a_{31}x_1 + $$ $$ a_{12}x_1 x_2 + a_{22}x_2^2 + a_{32} x_2+ $$ $$ a_{13} x_1 + a_{23} x_2 x_3 + a_{33} $$

from scipy.interpolate import Rbf
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
import autograd.numpy as anp
import autograd

def random_ball(num_points, dimension, radius=1):
    from numpy import random, linalg
    random_directions = random.normal(size=(dimension,num_points))
    random_directions /= linalg.norm(random_directions, axis=0)
    random_radii = random.random(num_points) ** (1/dimension)
    return radius * (random_directions * random_radii).T

np.random.seed(0)
N = 20

def rosenbrock(x):
    return (1 + x[0])**2 + 100*(x[1] - x[0]**2)**2

def Rosenbrock(x,y):
    return (1 + x)**2 + 100*(y - x**2)**2

def get_fvals_in_region(xcurr, f, radius):    
    b = random_ball(N, 2, radius)
    pts = xcurr+b
    vals = [f(p) for p in pts]
    return xcurr+b, np.array(vals)

x0 = [1.5,0]
xs,vs = get_fvals_in_region(x0, rosenbrock, 0.5)

res = []
for i in range(vs.shape[0]):
    res.append((xs[i,0],xs[i,1],vs[i]))
res = np.array(res).reshape(vs.shape[0], 3)

x = np.linspace(-2,2,250)
y = np.linspace(-1,3,250)
X, Y = np.meshgrid(x, y)
Z = Rosenbrock(X, Y)

import itertools
import numpy.linalg as lin

def quad_interpolate(xi, yi):
    xi = np.hstack((xi, np.ones((1,len(xi))).T  ))
    #print (xi)
    D = xi.shape[1]
    print (D)
    X_train = []
    for row in xi:
        X_train.append([row[i]*row[j] for i,j in itertools.product(range(D),range(D)) ])
    X_train = np.array(X_train)
    print (X_train.shape)
    print (yi.shape)
    coef,_,_,_ = lin.lstsq(X_train, yi)
    return coef

xi = res[:,[0,1]]
yi = res[:,[2]]
coef = quad_interpolate(xi,yi)

print (coefs)

x = np.linspace(-2,2,250)
y = np.linspace(-1,3,250)
X, Y = np.meshgrid(x, y)
Z = Rosenbrock(X, Y)

fig = plt.figure(figsize = (8,4))
ax = fig.gca(projection='3d')
ax.plot3D(res[:,0],res[:,1],res[:,2],'r.')
ax.plot_surface(X,Y,Z,rstride = 5, cstride = 5, cmap = 'jet', alpha = .4, edgecolor = 'none' )

def q_interp(x1,x2):
    x = np.array([[x1,x2,1]])
    A = coef.reshape(3,3)
    res = np.dot(np.dot(x,A),x.T)
    return np.float(res)

Zi = np.array([q_interp(xx,yy) for xx,yy in zip(X.flatten(),Y.flatten())])
Zi = Zi.reshape(X.shape)
ax.plot_wireframe(X,Y,Zi)

coefs = coef.reshape(3,3)

g = (2 * np.dot(coefs[:2,:2],np.array(x0).reshape(2,1)))

gnorm = g / np.sum(g)

ax.set_zlim(0,2500)

ax.quiver(x0[0], x0[1], 0, -gnorm[0], -gnorm[1], 0, color='red')

hess = 2*coefs[:2,:2]
print (hess)
newton_dir = -np.dot(lin.inv(hess),g)
print (newton_dir)

d = newton_dir
print (d)

ax.quiver(x0[0], x0[1], 0, d[0], d[1], 0, color='green')

ax.plot3D([x0[0]], [x0[1]], [0.0], 'b.')

ax.view_init(21, -133)

Edit.

Best Answer

Related Solutions

Multivariate Quadratic Regression, Surface Fitting and the Hessian

Related Question