This notebook presents the basic types of finite differences and interpolation methods that can be considered on a cartesian grid. The tools are presented in two dimensions, but apply in arbitrary dimension. They also apply in the context of:
Summary of volume Algorithmic tools, this series of notebooks.
Main summary of the Adaptive Grid Discretizations book of notebooks, including the other volumes.
Acknowledgement. Some of the experiments presented in these notebooks are part of ongoing research with Ludovic Métivier and Da Chen.
Copyright Jean-Marie Mirebeau, Centre Borelli, ENS Paris-Saclay, CNRS, University Paris-Saclay
import sys; sys.path.insert(0,"..") # Allow imports from parent directory
#from Miscellaneous import TocTools; print(TocTools.displayTOC('FiniteDifferences','Algo'))
from agd import FiniteDifferences as fd
from agd import AutomaticDifferentiation as ad
from agd import Interpolation
from agd import LinearParallel as lp
norm_infinity = ad.Optimization.norm_infinity
import numpy as np
from matplotlib import pyplot as plt
In order to test the finite difference and interpolation methods, we need some polynomial functions.
def u1(X): return X[0]+2*X[1]
def u2(X): return X[0]**2+2*(2*X[0]*X[1])+3*X[1]**2
def u3(X): return X[0]**3+X[0]*X[1]**2
def u123(X): return ad.array( (u1(X),u2(X),u3(X)) )
We also need to choose a direction, with integer coordinates, for the finite differences.
e = (1,2)
Let us also define a domain, here a square.
aX,h = np.linspace(-1,1,retstep=True)
X=np.array(np.meshgrid(aX,aX,indexing='ij'))
shape = X.shape[1:]
def interior(shape,k):
"""Boolean array excluding k boundary layers"""
interior = np.full(shape,True)
for i in range(len(shape)):
interior[*(slice(None),)*i,:k]=False
interior[*(slice(None),)*i,-k:]=False
return interior
def close(u,v,k,ndim=2,**kwargs):
"""
Wether u and v are close in the domain minus k boundary layers.
- **kwargs : passed to np.allclose
"""
dom = interior(u.shape[-ndim:],k)
# print(norm_infinity(u[...,dom]-v[...,dom]),kwargs)
return np.allclose(u[...,dom],v[...,dom],**kwargs)
interior((5,5),1)
array([[False, False, False, False, False], [False, True, True, True, False], [False, True, True, True, False], [False, True, True, True, False], [False, False, False, False, False]])
The following variables are used for validation, by comparison with automatic differentiation.
X_ad = ad.Dense.identity(constant=X,shape_free=(2,))
X_ad2 = ad.Dense2.identity(constant=X,shape_free=(2,))
du1 = u1(X_ad).gradient()
du2 = u2(X_ad).gradient()
ddu2 = u2(X_ad2).hessian()
du3 = u3(X_ad).gradient()
ddu3 = u3(X_ad2).hessian()
_e = fd.as_field(e,shape)
du1_e = lp.dot_VV(du1,_e)
du2_e = lp.dot_VV(du2,_e)
ddu2_e = lp.dot_VAV(_e,ddu2,_e)
du3_e = lp.dot_VV(du3,_e)
ddu3_e = lp.dot_VAV(_e,ddu3,_e)
Du1_e = fd.DiffUpwind(u1(X),e,h)
assert close(Du1_e,du1_e,2)
DDu2_e = fd.Diff2(u2(X),e,h)
assert close(DDu2_e,ddu2_e,2)
Centered finite differences are not degenerate elliptic, but they can be used within degenerate elliptic schemes if they are suitably compensated by second order finite differences. $$ \frac{u(x+h e)-u(x-h e)} {2 h} = <\nabla u(x),e> + O(h^2). $$
Du2_e = fd.DiffCentered(u2(X),e,h)
assert close(Du2_e,du2_e,2)
High order finite differences are not degenerate elliptic. They may be used within filtered schemes, which combine a stable degenerate elliptic scheme with a higher order scheme. By construction, they are exact on higher order polynomials.
Du2_e = fd.DiffUpwind(u2(X),e,h,order=2)
assert close(Du2_e,du2_e,4)
Du3_e = fd.DiffUpwind(u3(X),e,h,order=3)
assert close(Du3_e,du3_e,6)
DDu3_e = fd.Diff2(u3(X),e,h,order=4)
assert close(DDu3_e,ddu3_e,4)
Du3_e = fd.DiffCentered(u3(X),e,h,order=4)
assert close(Du3_e,du3_e,6)
The following finite differences can be used to estimate numerically the derivatives of a function, but they are rarely adequate for building numerical schemes. We denote by $e_i$ the $i$-th element of the canonical basis.
Du2 = fd.DiffGradient(u2(X),gridScale=h)
assert close(Du2,du2,2)
Du2.shape
(2, 50, 50)
for all $0\leq i < d$, and $$ \frac{u(x+h e_i+h e_j)+u(x-he_i-h e_j)-u(x+h e_i -h e_j) - u(x-h e_i+he_j)}{4h^2} = \frac {\partial^2 u} {\partial x_i \partial x_j} + O(h^2), $$ for all distinct $i,j$.
DDu2 = fd.DiffHessian(u2(X),gridScale=h)
assert close(DDu2,ddu2,2)
DDu2.shape
(2, 2, 50, 50)
DDu3 = fd.DiffHessian(u3(X),gridScale=h,order=4)
assert close(DDu3,ddu3,4)
The agd library contains a partial reimplementation of the scipy.ndimage.map_coordinates
function, which allows AD types, both for the coordinates and the input values. This is another avenue for numerically differentiating a function defined on a grid.
Boundary conditions.
We only support reflected (default) and periodic boundary conditions (grid-mirror
, grid-wrap
).
As a result, the spline interpolation methods does not exactly reproduce polynomials. In the tests below, we exclude some boundary layers, and choose a sufficiently high tolerance, to account for this non-exactness.
Exact reproduction of low degree polynomials could be achieved with the not a knot
boundary conditions, but these are substantially more difficult to implement, and numerically coslty.
Let us define a finer grid.
aX_ = np.linspace(-1,1,80)
X_=np.array(np.meshgrid(aX_,aX_,indexing='ij'))
shape_ = X_[0].shape
Interp = Interpolation.UniformGridInterpolation
(Piecewise) Linear splines are continuous, and reproduce linear functions. They are second order consistent $$ U_1^h(x) = u(x)+O(h^2). $$
U1 = Interp(X,u1(X),order=1)
assert np.allclose(U1(X_),u1(X_))
The spline can be differentiated, and yields the a first order consistent estimation of the gradient, except possibly at boundary points. $$ \nabla U_1^h(x) = \nabla u(x) + O(h) $$
dU1 = U1(X_ad).gradient()
assert close(dU1,du1,1)
Quadratic splines are continuously differentiable, and reproduce quadratic functions. We use a not-a-knot boundary condition : in one dimension, the second derivative is continuous accross the second node from the left. $$ U_2^h(x) = u(x)+O(h^3). $$
U2 = Interp(X,u2(X),order=2)
assert close(U2(X_),u2(X_),k=10,atol=1e-6)
The spline can be differentiated, one time or two times, and yields a second order consistent estimate of the gradient, and a first order consistent estimate of the hessian, except possibly at boundary points. $$ \nabla U_2^h(x) = \nabla u(x) + O(h^2), $$ $$ \nabla^2 U_2^h(x) = \nabla^2 u(x) + O(h). $$
dU2 = U2(X_ad).gradient()
assert close(dU2,du2,k=10,atol=1e-6)
ddU2 = U2(X_ad2).hessian()
assert close(ddU2,ddu2,k=10,atol=1e-6)
Cubic splines are twice continuously differentiable, and reproduce cubic functions. We use a not-a-knot boundary condition : in one dimension, the third derivative is continuous accross the second node from the left, and likewise from the right. $$ U_3^h(x) = u(x)+O(h^4). $$
U3 = Interp(X,u3(X),order=3)
assert close(U3(X_),u3(X_),k=10,atol=1e-6)
The spline can be differentiated two times, and yields consistent estimates of the gradient and hessian, except possibly at boundary points. $$ \nabla U_3^h(x) = \nabla u(x) + O(h^3), $$ $$ \nabla^2 U_3^h(x) = \nabla^2 u(x) + O(h^2). $$
dU3 = U3(X_ad).gradient()
assert close(dU3,du3,k=10,atol=1e-6)
ddU3 = U3(X_ad2).hessian()
assert close(ddU3,ddu3,k=12,atol=1e-6)
U123 = Interp(X,u123(X),order=3)
assert close(U123(X_),u123(X_),k=12,atol=1e-6)
As announced above, our spline interpolation function (partially) reproduces scipy.ndimage.map_coordinate
, and extends it to AD types. We make this transparent by providing a map_coordinates function.
agd_map_coordinates = Interpolation.map_coordinates
We also provide a small wrapper over ndimage.map_coordinates, which allows AD types for the input values (not the coordinates), as well as the interpolation of tensor data.
ndimage_map_coordinates = Interpolation.ndimage_map_coordinates
The implementations works both on the cpu and gpu, select by commenting/uncommenting the cell below. In the latter case, scipy.ndimage.map_coordinates
is replaced with cupyx.scipy.ndimage.map_coordinates
.
caster = lambda x:x # CPU
#caster = ad.cupy_generic.cupy_set # GPU
np.random.seed(42)
for i,(order,dom_shape,x_shape,c_shape,size_ad,mode) in enumerate([
# One dimensional tests
(1,(5,),(8,),(2,3),1,'reflect'),
(2,(6,),(4,),(3,),1,'reflect'),
(3,(4,),(2,),tuple(),2,'grid-mirror'),
(4,(10,),(3,),(1,),-1,'grid-wrap'),
(5,(10,),(3,2),(1,),-1,'grid-wrap'),
# Two dimensional tests
(1,(4,5),(3,),tuple(),-1,'grid-wrap'),
(2,(4,6),(3,4),(5,),2,'grid-wrap'),
(3,(4,6),(3,4),(5,),2,'reflect'),
(5,(4,3),(2,1),(2,1),1,'reflect'),
# Three dimensional tests
(1,(4,5,6),(2,5),(2,3),1,'reflect'),
(3,(6,2,3),(1,4),tuple(),2,'grid-mirror'),
]):
print(f"--- {order=}, {dom_shape=}, {x_shape=}, {c_shape=}, {size_ad=} ---")
vals = np.random.rand(*c_shape,*dom_shape)
if size_ad>0: vals=ad.Dense.denseAD(vals,np.random.rand(*vals.shape,size_ad))
pos = np.random.rand(len(dom_shape),*x_shape)*3*np.max(dom_shape)
vals,pos = map(caster,(vals,pos))
impl = agd_map_coordinates(vals,pos,order=order,mode=mode)
scpy = ndimage_map_coordinates(vals,pos,order=order,mode=mode) #scipy.ndimage
print("Absolute error : ",norm_infinity(impl-scpy))
assert np.allclose(impl,scpy,atol=1e-2)
# Note : the accuracy of scipy.ndimage uses float32 internally, hence its rather low accuracy
pos = np.array(np.meshgrid(*[np.arange(s,dtype=np.float64) for s in dom_shape], indexing='ij'))
pos = caster(pos)
impl = agd_map_coordinates(vals,pos,order=order,mode=mode)
scpy = ndimage_map_coordinates(vals,pos,order=order,mode=mode)
print("Reproduction error. impl:",norm_infinity(ad.remove_ad(impl-vals))," scpy:",norm_infinity(ad.remove_ad(scpy-vals)))
if ad.is_ad(impl) and impl.size_ad>0:
print("Reproduction max error. impl:",norm_infinity((impl-vals).coef)," scpy:",norm_infinity((scpy-vals).coef))
--- order=1, dom_shape=(5,), x_shape=(8,), c_shape=(2, 3), size_ad=1 --- Absolute error : denseAD(0.0,[0.]) Reproduction error. impl: 0.0 scpy: 0.0 Reproduction max error. impl: 0.0 scpy: 0.0 --- order=2, dom_shape=(6,), x_shape=(4,), c_shape=(3,), size_ad=1 --- Absolute error : denseAD(7.083111874806036e-12,[4.10549372e-12]) Reproduction error. impl: 2.220446049250313e-16 scpy: 1.1591035076197187e-10 Reproduction max error. impl: 2.220446049250313e-16 scpy: 6.718436917907411e-11 --- order=3, dom_shape=(4,), x_shape=(2,), c_shape=(), size_ad=2 --- Absolute error : denseAD(1.4526723157715082e-09,[-2.47407976e-07 -4.96368829e-07]) Reproduction error. impl: 1.3877787807814457e-17 scpy: 2.307997604145129e-08 Reproduction max error. impl: 2.220446049250313e-16 scpy: 7.886280012692204e-06 --- order=4, dom_shape=(10,), x_shape=(3,), c_shape=(1,), size_ad=-1 --- Absolute error : 5.551115123125783e-17 Reproduction error. impl: 3.3306690738754696e-16 scpy: 2.220446049250313e-16 --- order=5, dom_shape=(10,), x_shape=(3, 2), c_shape=(1,), size_ad=-1 --- Absolute error : 6.661338147750939e-16 Reproduction error. impl: 2.220446049250313e-16 scpy: 4.440892098500626e-16 --- order=1, dom_shape=(4, 5), x_shape=(3,), c_shape=(), size_ad=-1 --- Absolute error : 1.1102230246251565e-16 Reproduction error. impl: 0.0 scpy: 0.0 --- order=2, dom_shape=(4, 6), x_shape=(3, 4), c_shape=(5,), size_ad=2 --- Absolute error : denseAD(4.440892098500626e-16,[4.4408921e-16 1.2490009e-16]) Reproduction error. impl: 3.3306690738754696e-16 scpy: 2.220446049250313e-16 Reproduction max error. impl: 3.3306690738754696e-16 scpy: 3.3306690738754696e-16 --- order=3, dom_shape=(4, 6), x_shape=(3, 4), c_shape=(5,), size_ad=2 --- Absolute error : denseAD(8.41014820128494e-06,[3.36075072e-06 8.63728999e-06]) Reproduction error. impl: 5.551115123125783e-16 scpy: 8.262379805445974e-06 Reproduction max error. impl: 4.440892098500626e-16 scpy: 9.46352712638543e-06 --- order=5, dom_shape=(4, 3), x_shape=(2, 1), c_shape=(2, 1), size_ad=1 --- Absolute error : denseAD(0.0017610646499264804,[0.00317475]) Reproduction error. impl: 3.3306690738754696e-16 scpy: 0.0031852717644261785 Reproduction max error. impl: 3.3306690738754696e-16 scpy: 0.0042962930820822365 --- order=1, dom_shape=(4, 5, 6), x_shape=(2, 5), c_shape=(2, 3), size_ad=1 --- Absolute error : denseAD(2.220446049250313e-16,[1.11022302e-16]) Reproduction error. impl: 0.0 scpy: 0.0 Reproduction max error. impl: 0.0 scpy: 0.0 --- order=3, dom_shape=(6, 2, 3), x_shape=(1, 4), c_shape=(), size_ad=2 --- Absolute error : denseAD(0.0006403965303843351,[0.00030922 0.00044751]) Reproduction error. impl: 5.551115123125783e-16 scpy: 0.000487351300054617 Reproduction max error. impl: 5.551115123125783e-16 scpy: 0.0004580943068471788
The automatic differentiation classes store a Taylor expansion of the approximated function, which can be evaluated directly. The tangent, adjoint, and hessian operators may also be extracted.
x = np.array([0.2,0.5])
H = X-fd.as_field(x,shape)
x_ad = ad.Dense.identity(constant=x)
x_ad2 = ad.Dense2.identity(constant=x)
plt.title("Absolute difference between U3 and its first order taylor expansion")
plt.contourf(*X, np.abs(U3(x_ad).as_func(H) - U3(X)), levels=20);
plt.axis('equal'); plt.scatter(*x,color='red'); plt.colorbar();
plt.title("Absolute difference between U3 and its second order taylor expansion")
plt.contourf(*X, np.abs(U3(x_ad2).as_func(H) - U3(X)), levels=20);
plt.axis('equal'); plt.scatter(*x,color='red'); plt.colorbar();
When differentiating a high dimensional function, for instance in the context of a PDE discretization scheme, sparse automatic differentiation becomes for reasons of memory and computation cost. We check here the consistency of the related Taylor expansions.
x_sp = ad.Sparse.identity(constant=x)
x_sp2 = ad.Sparse2.identity(constant=x)
assert np.allclose(U3(x_sp ).as_func(H), U3(x_ad ).as_func(H))
assert np.allclose(U3(x_sp2).as_func(H), U3(x_ad2).as_func(H))
The sparse classes also provide tangent, adjoint, and hessian linear operators, stored as opaque sparse matrices.
tangent_op = U3(x_sp ).tangent_operator()
hessian_op = U3(x_sp2).hessian_operator()
H_ = H.reshape(2,-1) # Depth must be at most two
assert np.allclose(U3(x_sp ).as_func(H_), U3(x) + tangent_op*H_)
assert np.allclose(U3(x_sp2).as_func(H_), U3(x) + tangent_op*H_ + 0.5*lp.dot_VV(H_,hessian_op*H_))
Finally, we check the adjoint operator.
adjoint_op = U3(x_sp).adjoint_operator()
np.random.seed(42)
R_ = np.random.rand(1,H_.shape[1])
assert np.allclose(lp.dot_VV(R_,tangent_op*H_), lp.dot_VV(H_,adjoint_op*R_))