Eikonal

A modal spectral method for the Eikonal equation applied to continuous shortest path problems on simplicial tesselations.

The primarily goal in developing this library is for use in solving continuous shortest path problems, with application in navigation, infrastructure design, and GR. This README is more of a stream of conciousness blogging of progress, obstacles, solutions, and new ideas, all as they relate to this project.

About

The implementation offers a backend library for computing near-optimal Gaussian-like quadratures on the simplex in dimensions 1,2, and 3 (see ngjquad.cpp for driver code), built on top of NLOPT, BLAS and LAPACK, my own work, theoretical results due to Koornwinder regarding a family of orthogonal polynomials on the triangle, and numerical results due to Vioreanu and Rokhlin on constructing near-optimal Gaussian-like quadratures over convex regions (VR-quadrature).

I term the VR-quadrature approach "Discretize-Orthogonalize-Optimize" (DOO), while my approach is more "Orthogonalize-Discretize-Optimize" (ODO). This nomenclature draws analogy from similar jargon in variational/infinite-dimensional optimzation - Discretize-then-Optimize (DtO) or Optimize-then-Discretize (OtD) - as both concern a choice of an analytical or numerical first step.

The generatable quadrature rules are used to represent to high accuracy and with high efficiency functions appearing in Eikonal models under the Koornwinder polynomial basis family.

Recent results by Townsend, Oliver and Vasil enable the construction of sparse differential operators on the standard triangle, interoperating within the Koornwinder polynomial family parameter space, and leading to low-storage, banded discrete differential operators. The library provides an implementation of these discrete PDOs. They can be composed together to form solvers for a variety of linear PDEs of arbitrary order, even with variable coefficients (coming soon (i hope)). Non-linearity is handled by passing an appropriately defined problem to nlopt in terms of these operators acting on minimization variables (solution coefficients in a Koornwinder expansion).

At the time of writing, the following capabilities exist:

• representing and manipulating functions in a modal sense under Koornwinder polynomial expansions (global approximation) as well the so-called Jacobi polynomials on the Tetrahedron.

• differentiation up to any order via promotion and ladder operators acting only on modes of a function on the triangle

• evaluation at a point in real space of a function (dim=1,2,3) represented in the modal basis (coefficients). This could be accelerated with Clenshaw's algorithm, generalized to 2D and 3D

• Implementation of analytical entries to the Jacobi matrices appearing in the 2D-recurrence relation for the orthogonal Koornwinder polynomials, generalized for any choice of parameters. The 3D ones remain to be derived.

• Quadrature discovery with optimization parameters from nlopt exposed for such use. In particular, I developed an interlaced scheme where first a constrained non-linear optimization problem is solved, then, if needed, a Newton relaxation, followed by another non-linear problem with a different objective, and a last on-demand Newton relaxation.

  • The gist is: With z=(x,y,w) in R^(3N), minimize over z norm(P(x,y).w - e1)^2, where e1 is the first basis vector in R^M with N < M, P(x,y) is the NxM Vandermonde matrix, (x,y) is constrained to lie within the triangle, and sum(w) is constrained to 1, followed by a Newton relaxation to find roots of P(z).w - e1 = 0. Then the argmin is passed to another constrained non-linear problem, wherein we try to minimize cond(P(x,y)) (in the 2-norm) with w fixed, and constrain norm(P(x,y).w-e1) to be at most the objective minimum from the first optimization problem. This seems to find Gaussian-like quadrature rules at the limit proposed by Vioreanu and Rokhlin (in terms of how much N < M), but with the added benefit of O(1) condition number for the interpolation matrix. This latter feature is important for modal discretizations of PDEs, as when we impose boundary conditions nodally, the interpolation operator on those nodes must be well conditioned if we hope for the assembled systems to be invertible.

For example, the image below depicts the output when attempting to generate an N = 253 (n max = 21) node quadrature rule which can exactly integerate 630 polynomials (m max = 34). In this case, we use the Subplex algorithm, strict tolerances of 1e-14 on the relative change in iterates of the minimization varialbes, and 6 threads.

The red nodes are our new found abscissa, while the black ones are the initial points. We see the new nodes are well conditioned, and we correctly integrate the test function sin(x^2+y^2) over the standard right triangle with vertices (0,0), (0,1), (1,0) (plug it into Wolfram alpha for reference). See triaintq.f for comparison at https://github.com/zgimbutas/triasymq. We exceed the largest rule generated there, and ours is better conditioned by 2 two orders of magnitude.

There remains much to be done for this project, on both paper and metal!

TODO Ideas for Quadrature Generation

We seeks roots to the function F(x,y,w) = P(x,y).w - I = 0, where P is the transpose of the vandermonde matrix P(x,y)^T, and (x,y) are a set of points in the standard right triangle. There is some theoretical limit on the solvability of the equation in terms of the number of points (#columns of P) and the number of integrals for which we want exact quadrature (#rows of P or length of I). We know it is not solvable at the optimal Gaussian quadrature limit, as the Jacobi matrices for Koornwinder polynomials on the simplex are not a commuting family. The limit is explored numerically for more general convex regions by Rokhlin and Vioreanu, and an empirical cutoff formula is provided in their work. Assuming we are under this limit in terms of the system size, and further assuming a root exists (ignoring conditioning of P), then ||F||^2 in the 2-norm should also admit a zero at some point (x,y,w). Rather than trying to minimize ||F||^2 using non-linear optimization techniques like SQP, or any method where quadratic approximations to sub-problems are used, we might be as well successful building a quadratic approximation of the objective G=||F||^2 and using the Newton direction -(\nabla^2 G(x,y,w))^{-1}\nabla G(x,y,w) to find the root. The only "approximation" to be made here in the model is for the gradient and hessian entries concerning variable w. We can hopefully use simple centered finite difference schemes there. Variables x,y enter only in P, and we should know by now how to take derivatives of polynomials. Ill-conditioning in the argmin can be handled with our current approach (nlopt call on min cond(P)).

• the first step is seeing if incorporating analytical gradient information even helps convergence in the nlopt routines currently used. After all, if it isn't broken, don't fix it, and we are already able to generate high order quadrature with a large increase in efficiency relative to the past (days->hours on the same OTS equipment). I will need to check if the ones currently in use even need derivatives (many are derivative-free algorithms), and whether there are others that require derivative info.
• if nlopt responds positively to analytical gradient info, then we can just plug and play the Hessian into my current Newton method implementation, discarding the backtracking line search part. Let's see what happens!

Newton's method has fast quadratic convergence, but there is difficulty of step length control for gradient descent which leads to constraint violations for the quadrature rule. If we can get Newton descent right, it may open up the possibility for on-line p-adaptivity in PDE discretizations using the Koornwinder bases (which also reminds me of the possiblity of (a,b,c) adaptivity - interoperating between basis families parameterized by (a,b,c)).

TODO Ideas for generalizing to d>2

The only challenge with all of this has been in deriving matrix entries for the 3-term matrix recurrence relations defining spaces of polynomials orthogonal to all spaces of lower order (for d > 1). We need these entries to assemble Jacobi matrices, and we need those to initialize any descent method for finding Gaussian-like quadrature/cubature. By definition, Jn_i.P = x_i.P + [0\\A_{n-1}_i \bb{P}_n]. In 2D, if {x_1,x_2}_j is a set of nodes admitting a well conditioned interpolation matrix on the triangle, then we should be able to solve to mach_eps for the entries in Jn_1,Jn_2. All we would need to do so (and, in fact, already have) are the entries to A_{n-1}_{1,2}. In higher dimensions, this reduces the courageous pencil/paper work by quite a bit, even if we use symbolic algebra systems. Then, we can consider JEVD techniques used in SP problems like BSS to generalize the initialization method for d>2 (complexification doesn't work for 3D, which is used in 2D for the eigenvalue matching problem. Maybe we can use quaternions and projection to dim under for 3D? Is there even an eigenvalue solver that supports quaternions??).

Computing approximate joint eigenvalue decompositions

Consider orthogonal polynomials up to total degree $n-1$ in dimension $d$ defined in a convex region $\Omega \subset \mathbb{R}^d$ with boundary $\partial \Omega$. There are $\text{dim}\Pi_{n-1}^d=N$ of them. Suppose that the symmetric, real-valued Jacobi matrices $J^k \in \mathbb{R}^{N\times N}$, $k=1,\cdots d$ corresponding to a 3-term recurrence relation on subspaces of homogoneous polynomials with total degree $m=0,\cdots n-1$ do NOT form a commuting family. Then, there exists an approximate joint eigenvalue decomposition of the matrices $J_k$ such that the $d$-tuple of matched eigenvalues $(\lambda_m^1,\cdots \lambda_m^d)\in (\Omega \setminus \partial\Omega)$ (not proven, but seems to hold numerically! We need it to hold only finite manly times, so that's ok.).

That is, we claim that there exists a single unitary matrix $U$ such that the matrices $$D^k = U J^k U^T, \quad 1 \leq k \leq d$$ are approximately diagonal. Here, we define approximate diagonality in the sense that $D^k$ is the best representation of $J^k$, entry-wise, in the basis $U$. From a subspace fitting approach, we can cast the task of finding $U$ and $D^k$ into a quadratic optimization problem with quadratic equality constraints:

Solve

$$ \min_{U,F^1,\cdots F^d} \sum_{k=1}^d ||J^k - U^T F^k U||_F^2 $$

constrained by

$$U^TU=UU^T=I,$$

where $I \in \mathbb{R}^{N \times N}$ is the identity matrix.

The idea is to throw this into an optimization routine on vectorizations of the minimization variables, so there are $N^2(d+1)$ of them. The main question has to do with initialization:

• Ensure the constraint is satisfied to begin with by setting $U=I$.
  • Try initializing $F^k$ with a perturbation of $J^k$ drawn randomly from a distribution to be determined.
  • Try initializing $F^k = (J^k)_{ii}$ the diagonal elements of $J^k$.
• Since we seek an 'average eigenstructure' for $(J^k)$, initialize the QR decomposition of one of the $J^k$.

A better idea would be to use an old, but working implementation (Jacobi angles for simultaneous diagonalization), and port it into the library with threading/simd additions. In fact, this has already been implemented with simd optimizations in the library, and tested successfully (see include/jevd.h and gtest code for testing). There are some challenges with threading the algorithm, mainly to do with an appropriate packing of the input matrices into a combined matrix so as to minimize cache incoherencies, and maximize the number of loops the compiler deems vectorizable. There's also some more thinking that needs to be done in terms of restructuring the algorithm to avoid race conditions on matrix updates, but not require any critical barriers on execution per thread. For now, I'm happy with a serial but vectorized code that can handle $d\geq 2$.

UPDATE

UPDATE: It works! JEVD can be used to initialize nodes for Gaussian-like quadrature optimization on the Tetrahedron. See the image below, which I generated for n=5. Implementation of the definitions for the Jacobi polynomials on the tetrahedron remains, along with the structural constants, Jacobi matrix block defintions, etc (UPDATE: all that is done too). However, I have written procedures in Mathematica to symbolically generate most of this stuff. I need to tweak that code a bit so that it gives me nicer formulas, which I can then just port into C (I rly don't like algebra).

The jacobi matrices, the approximate joint eigenvalues of which appear as the nodes above, have interesting expanding block band sparsity patterns ($J_x,J_y,J_z$, resp., depicted below):

Since I've been waiting for 2 days on the Wolfram Kernel (using 90% of my computer's memory) to evaluate Jacobi matrix entries on the tetrahedron for n = 15, I had to do some thinking..

I realized there exists a straightforward way of numerically computing the entries to the Jacobi matrices for the tetrahedron. However, it relies on getting the near Gaussian-like quadrature optimization working correctly. The idea is to bootstrap:

• Generate Jacobi matrices for a low order, like n = 4.
• Generate the initial nodes for adapting the quadrature rule via JEVD
• Push this adapatation as far as we can, resulting in a higher order rule using the same number of nodes.
• Use this new rule to approximate the required weighted inner products for entries to the Jacobi matrix, but for a higher order, like n = 5 (if we started with n = 4).
• Repeat this process until we get SUPER high order!!

The justification here is that we will certainly be able to compute the required integrals for the next order, to within some relatively acceptable (but still unacceptable) epsilon. Now, if this epsilon error is small enough, the JEVD process should still produce nodes that live in the tetrahedron, and would serve well for initializing the quadrature optimization routine. Conditioning is not a concern here. We only want to integrate accurately, so we can eventually have accurate interpolation and function expansion. Think of the epsilon in terms of the epsilon-commutativity of the Jacobi matrices. We never really find exact joint eigenvalues, as they don't exist in this case! We only find the joint eigenvalues of $J_{ij}^k + \epsilon_{ij}^k$ for $k=1,2,...,d$. So now, we are approximating the approximate joint eigenvalules, but iterating towards the nearly-epsilon-commuting family of Jacobi matrices.

Remark: I'm assuming that this epsilon commuting family exists. That is, I assume the Jacobi matrices live on some kind of manifold of non-commuting matrices, but are just far enough away from a point on a manifold of commuting matrices. So, I assume that we need only push these matrices a small amount to be on the manifold of commuting matrices, and their projections on that manifold yield eigenvalues that live in the tetrahedron. The vaguery of language here is intentional. The point is that we should still be able to push the matrices to the point with eigenvalues which are in a basin of convergence for the underdetermined linear system / quadratic root finding problem we seek to solve at the end of the day.

UPDATE

(UPDATE:) EHHHH!!! - we can't bootstrap because the inner products we need to compute involve order 2n+1 total degree polynomials. BUT, we can kinda-sorta bootstrap by incorporating the below process, when said process starts lacking in efficacy.

(FIX THE WRONG) We can use mapped quadrature rules to achieve the end goal of computing these inner products. I have defined a deformation mapping between the unit cube (in R3+) and the standard-right tetrahedron. By making use of the fact that in $\mathbb{R}^{3+}$, both regions can be expressed in terms of norms, and exploiting the equivalence of norms in finite dimensions up to a constant (which can depend on x, but is uniformly bounded), we have a simple mapping : C2T : $[0,1]^3 \to T$ which takes $x \in [0,1]^3$ and maps it to $y \in T$ ( there's a stack exchange post on this matter). The trouble here is in computing the deformation gradient of this map. It's Frankenstein's monster's little brother analytically, so I just use centered second order finite differences to compute each component of the deformation gradient. The integral mapping then appears as:

$ \int_T f(y) dy= \int_{[0,1]^3} f(C2T[x]) det[\nabla_x C2T]dx$

where the det term is the Jacobian determinant, or change of variable form. The utility here is that we have recast an integral over a domain for which we do not have a sufficiently accurate quadrature rule into an integral over a domain for which we can use a tensor product of simple 1D quadrature rules. In my case, I just generate the (0,0) Jacobi polynomials (i.e. Legendre basis), and use the Golub-Welsch algorithm to compute an optimal Gaussian quadrature rule of order $2n-1$ with only $n$ nodes. This can be made even fancier with Smolyak's quadrature approach for the case of symmetric weight functions (which lead to symmetric nodes, and nesting of nodes between certain orders). Such multi-level schemes might be the key if inefficiency prevails for higher order, but let me try and generate n=20 total order on the tet first. NOTE: the Pochhammer symbol and hypergeometric 2F1 series functions which appear in the definition of the Jacobi polynomials perform unstably for high order - (factorials and what not) - so we really can't go beyond a source basis order of 20-30. It becomes a real problem when we have to evaluate polynomials at such high order. (recall that the polynomials or d>1 are still definied in terms of the polynomials for d=1).

UPDATE

(UPDATE:) This worked swimmingly, and I can now generate high order quadrature on the tetrahedron relatively quickly! It turns out that the error in numerically approximating the inner products is within the tolerable error required for JEVD routines on the approximate Jacobi matrices to converge. There are research cookies here in terms of proving the relationship between the integration error, the off-diagonality minimization error, and interiority of matched eigenvalues to the simplex. I'll leave it to a numerical analyst to prove this stuff, though the relations are readily observable in numerical experiments.

I've moved on to making the library routines callable from Python, and playing around with solving Eikonal problems on the triangle. I migrated from plain-old make to CMake as the build system, which has made compilation a lot faster, and I guess I was linking things incorrectly, so the code runs faster too! I can add support for the Intel icc compilers too, but don't have the free version yet, so haven't bothered.

Docker Containerization

It is likely most simple to install and use the libarary from within a Docker container. For user convenience, a Dockerfile is provided which can generate a Docker image with all required dependencies and the Eikonal library installed. Note, the Docker Engine must be installed on your system (See https://docs.docker.com/engine/), along with git.

Administrative privileges may be required to run Docker commands. On Linux systems, you must add your user to the Docker group by executing:

sudo groupadd docker
sudo usermod -aG docker <user> # replace <user> with your username
newgrp docker # or log out and log back in

Since you need root privilege to execute the above commands, you could instead forgo their execution, and prefix docker commands with sudo.

Once these dependencies are met, executing (on Linux or Unix systems) the following commands will build the Docker base image for the project and test the installation within a running instance of the image (a Docker container):

git clone git@github.com:snatesh/Eikonal.git
cd Eikonal
docker build -t ngj_tri_opt:latest .
docker run -it ngj_tri_opt:latest bash
cd Eikonal
cd build && ctest --verbose

Above, the docker run command with bash post-fixed will instantiate a bash shell running in the container, within which you can read/write/execute files as you please. The last two commands above change directory into the /Eikonal folder in the container context, compile the gtest code (located in /Eikonal/testing/gtest.cpp), and run the tests therein.

Dependencies

• g++ compiler (tested on V13.2.0)
• build-essentials (particularly GNU make utility)
• google testing framework (available via dpkg and apt as libgtest-dev)
• dh-autoreconf,autoreconf,autotools for SNOPT installation (currently not used)
• NLOPT (open source non linear optimization library, https://github.com/stevengj/nlopt.git)
• cblas and lapack - Installation is easiest via package manager as:

sudo apt install libopenblas-openmp-dev liblapacke-dev

which ensures the C wrapper to lapack (in headers lapacke.h) is installed in a sane location.

• multipledispatch library in Python.

• The c++ compiler should have support for the OpenMP shared memory parallelization library, and the library must exist on your system. That is, omp.h and libomp.so must exist somewhere in the filesystem, the compiler must understand OpenMP directives, and the linker should be able to find and link to libomp given the -fopenmp flag. Most modern compilers will ship with the header and library files, as well as support/implementation of OpenMP directives. In case the files don't make it, you can use (on Linux with dpkg)

sudo apt install libomp-dev

Thread settings for OpenMP

An example threading config file cpuconfig.sh is included to show some of the environment variables that OpenMP exposes for users to set from the shell. In general, using the number of physical cores on the system improves the performance of most algorithsm that Eikonal uses (over hyperthreading), while the binding of threads spawned by OpenMP to those physical core IDs is something with which you should experiment on your system by setting the OMP_PROC_BIND and OMP_PLACES environement variables. Proc-binding essentially disables hyperthreading when you set the number of threads to less than or equal to the number of cores per socket.

I find that enabling threads significantly reduces the convergence time for quadrature search when n,m are large enough (>~10), so using and playing around with OpenMP settings is well worth the effort if you need high order quadrature.

# number of threads for OpenMP
num_threads=6

# let the shell use the maximum amount of stack memory
ulimit -s unlimited
# set mem for thread stack
export OMP_STACKSIZE=256m

#thread pinning settings
export OMP_PLACES="{0}:${num_threads}:1"
export OMP_PROC_BIND=true
export OMP_DISPLAY_ENV=true
export OMP_NUM_THREADS=${num_threads}

Note, the OMP_NUM_THREADS variable is set here, though this will change. It is not advisable to set such an environmental variable if you link to other programs which also use OpenMP. They may have their own tested/working heuristics for setting the number of threads. and fixing it in the shell context can mess up the performance of their threaded functions when called from within the same context. The num_threads variable set above will eventually be passed to the program at runtime (set by calling omp_set_num_threads), while OMP_NUM_THREADS will be unset/empty.

Python Bindings

The wrapper libraries are compiled by default, and the Python bindings are available in the python folder. You have to set the relevant environment variables so the interpreter can link to the libraries declared in the bindings files. From the Eikonal directory, run

export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:$PWD/lib
export PYTHONPATH=${PYTHONPATH}:/usr/local/lib/python3.12/site-packages

ensuring to replace the location of your Python installation accordingly. If using Docker, this is handled on image creation.

NLOPT Installation

The open-source nonlinear optimization libary which we use is nlopt, by our favorite FFTW co-creator Steven Johnson! For our purposes, it is used to generate near-optimal Gaussian-like quadrature on the triangle, and is to be investigated for use in solving non-linear PDEs/PDE-constrained optimization problems.

To install NLOPT, simply execute

git clone https://github.com/stevengj/nlopt.git
cd nlopt
mkdir build && cd build
cmake ..
make
sudo make isntall 

which (on Linux systems) will by default install header and shared library files in /usr/local/include and /usr/local/lib.

Threading in nlopt

Following the discussion of threading with OpenMP above, it seems that NLOPT is threaded using the lower level pthreads library, but responds to the OMP_NUM_THREADS variable setting (i.e. if num_threads=6, CPU utilization will not exceed 600%). I have yet to test whether this response is in favor of performance, in terms of time or memory use. The faster we can go, the farther we can push the order of generated quadratures, so this is worth looking into at some point.