Mat-Op-Ex Toolkit

Matrix Operator Extraction Toolkit - Tools to obtain the matrix form of common MATLAB linear operations (e.g. trapz, interp) performed over a 1-D or 2-D grid as generated in ndgrid. This solution is based heavily on a wonderful bit of code presented by Bruno Luong, and expanded here to be more generally applied.

Description

When performing repeated integration or interpolation of large data-sets over fixed grids, it is often advantageous to access the bare-metal matrix multiplications behind linear functions like interp1, interp2, and trapz. Mat-Op-Ex Toolkit seeks to provide these matrices efficiently, and in a way accessible to the user to build on.

Mat-Op-Ex functions are designed to be analogous to their conventional counterparts, except that they receive no function data and output a matrix instead. To avoid dealing with high-order tensor operations, physical dimensions are always vectorized within these functions, such that the operator matrix is limited to order N=2. Further details are available in the source code files, and the user should make use of the test scripts available.

Function List

• Interpolation
  • M=interp1_matrix(X,Xq) - 1D interp on 1D mesh
  • M=interp2_matrix(X,Y,Xq,Yq) - 1-2D interp on 2D mesh
• Integration
  • M=trapz_matrix1(X) - 1D int. on 1D mesh
  • M=trapz_matrix2(X,Y,DIM) - 1D int. on a 2D mesh
• Auxiliary
  • mesh_type - determines if X/Y inputs are invalid, vectors, meshgrid matrices, or ndgrid matrices

Example Test Figures

(A) A demonstration of the speed test for Mat-Op-Ex using trapz_matrix2 compared to common alternatives. Mat-Op-Ex on a CPU is the fastest option on small grids, and Mat-Op-Ex on a consumer GPU is the fastest on large grids (~3.5x faster than interp2).

(B) An example of the stencils used in bi-cubic interpolation in interp2_matrix, and the associated deviations from MATLAB interp2 with the cubic method.

Current Status

Mat-Op-Ex is currently tested to support 1D and 2D interpolation (nearest, linear, and cubic) and integration (trapz, simpsons rule) on 1-D and 2-D based grids that are evenly spaced along each axis (i.e. generated with ndgrid). All functions test within machine precision of MATLAB functions, except for the cubic interpolation on the edges of the domain, where a difference in boundary conditions incur a typical error of <1%.

Because the matrix operators depend heavily on the grid makeup, Mat-Op-Ex functions are less generalized than conventional counterparts, and expansion to higher dimensions (such as with interpn), are not straightforward. Mat-Op-Ex is similar in concept FUNC2MAT, but implements the schemes natively and is far more efficient.

Motivation

This project arose from the need to calculate the stability matrix for a certain set of discretized ODE's where interpolation and spatial integration played a critical role. In addition to furnishing the full matrix operator, the pure matrix form enabled highly-efficient solutions, especially using gpuArrays (which to present knowledge, must recompute stencils for each interpolation and integration call, unlike the cpu-based griddedInterpolant, which can store grid data). Any ideas to improve or expanded on this technique would be greatly welcomed.