\section{Introduction} \label{sec:intro}
Natural phenomena that are susceptible to a mathematical physics treatment typically lead to initial and boundary value problems. The solution of partial differential equations that govern these problems has occupied the better part of the last century of research in applied mathematics, leading to numerical methods, and from there to computational science and high performance computing. With Fourier methods, time integration algorithms, discretization frameworks of finite elements, finite difference, Green's functions-based formulations, and hybrids of these techniques among other approaches, it is now possible to attain a degree of fidelity that is sufficient for many applications in mathematical physics. While the need to incorporate additional physics and gain a finer resolution of phenomena will never be met, there is also the recognition that the fidelity that is already achievable comes paired with a complexity that outstrips easy representation and analysis. Reduced-order models must therefore be derived on a principled basis with controlled approximations.

For the purpose of distinguishing the approach adopted here, we take account of areas including projection methods, coarse graining techniques, dimensionality reduction, sub-manifold extraction and numerous others that are arguably encompassed by the field of dynamical systems. Within this field, is the graph theoretic approach for representation, exploration and analysis of computed states of physical systems that we have proposed recently \cite{Banerjee2019}. This formalism is predicated on computed solutions that are typically high-dimensional in the sense that the primary unknown fields are represented by degrees of freedom numbering $\gtrsim O(10^8)$. Functionals on these discrete solutions are extracted as \emph{states}. Some examples are the lift, drag or pressure difference in computational fluid dynamics, measures of load, deformation and total energies in solid mechanics, phase volume fractions and free energies in materials physics, and cross sections in nuclear physics \cite{Kochunas2020}. These states undergo \emph{transitions} as some \emph{parameter} of the system is varied, such as its physical, possibly time-like, parameters, boundary or initial conditions.

\subsection{Quantities of interest and graph representations} \label{sec:intro_qoi}
The states, of which there is typically a small number $\sim O(10)$ relevant to a system, form a low-dimensional vector of quantities of interest. However, these state vectors can be of high fidelity limited only by the dimensionality of the full-field numerical solution. For completeness of the description, the state vector can be extended to include the parameters referred to above. We have shown that by defining the states to be vertices on a graph, and the transitions to be edges between the vertices, a near one-to-one correspondence is uncovered between the properties of the computed physical system and the elements of graph theory. Weighted edges between vertices can then be assigned, or found through graph theoretic principles, indicating relationships or transitions between states, and the magnitude of such correlations. Reversible linear systems lead to undirected graphs that are fully connected, whereas dissipative dynamical systems are represented by directed trees. Path-dependence can be thought of as being associated with an absence of cycles. Notions of centrality reveal insights to the relations between states and path traversal properties are induced on graphs by equilibrium or dynamical transitions on the physical system.

Given the above graph theoretic framework we seek to extend it to reduced order modelling based on computed states and transitions of physical systems. We take the view of the graph of states (vertices) and transitions (edges) as a discrete manifold describing the physical system. On this manifold we seek to define differential equations and functional representations for the evolution of certain components of the state vector. An illustrative example that we pursue in this communication is of an ordinary differential equation for the phase fractions of a material system, driven by itself and other components of the state vector such as the free energy and strains. Guided by our knowledge of first-order systems in materials physics, it is of interest to compute derivatives on these states, such as of the free energy, with respect to the phase fractions, yielding generalized chemical potentials, or with respect to the strains, yielding stresses. The motivation comes from the rates of change of phase fractions driven by chemical potentials and stresses in materials. An alternative to the differential equation is a direct functional representation, such as for the free energy with phase fractions and strains as arguments. In this case, we are guided by the universality of the Taylor series representation, to which we seek approximations. This too makes necessary the computation of derivatives of state vector components with respect to each other.

\subsection{Non-local calculus for reduced order modelling on graphs} \label{sec:intro_nonlocalcalculus}
We have adopted a discrete non-local calculus on weighted graphs to develop the reduced order models outlined above. This formalism, which has been applied previously to image analysis \cite{Gilboa2008,Desquesnes2013,Hein2007} begins with definition of a vector space on the graph, followed by the gradient operator, inner products and partial derivatives. Also important is the consistency and rigor reflected in a divergence operator, adjointness relations and a divergence theorem. This non-local theory bears a relation to classical differential calculus, that is determined by the graph weights. Through its self consistently defined inner products, the weights appear as kernels determining the support of integral operators on the graph manifold, and crucially controlling the behavior of partial derivatives relative to their classical definitions.

This forms a central consideration on theoretical and practical grounds. However, the machinery that translates the non-local calculus on graphs to reduced order models can be formally put in place independently of the correspondence between non-local and local treatments. With functions and non-local partial derivatives in hand, algebraic and non-local differential operators can be introduced and viewed as basis terms in the reduced order representation. Since these can be evaluated at the vertex data points, a stepwise regression approach can be invoked to choose the best combination of operators representing the evolution or dependence sought. Model parsimony can be dialed in by using thresholds on the stepwise regression.

\subsection{Physical systems of interest} \label{sec:intro_physicalsystems}
The system of interest for this work are multi-crystalline solids, where the dynamical quantities of interest are the volume fraction of each phase, which are assumed to undergo first-order dynamics. These dynamics are selected to be forced by a linear functional of a $\ird[3]$ order polynomial of the phase volume fractions themselves, the interfacial lengths between the phases, the number of particles, and the first derivatives of the free energy, with respect to the phase volumes, and the mechanical strains. In order to complete the model of the system to be able to forward solve the first order dynamics, a functional representation for the free energy is found as a Taylor series in the phase volume fraction and mechanical strains, up to $\ith[4]$ order. This example allows the general graph theoretic method to be shown, and analyzed for its effectiveness at developing functional representations, and investigating any numerical issues associated with these derivatives. Further examples of physical systems will be explored in future communications.

Computations on each of these examples relies upon a general graph theory library code that accepts general data with inputs and outputs, and determines appropriate models based on these stepwise regression procedures, as well as visualizes fits with plots of loss curves, fitted curves for various models, and trends in the relative magnitudes of fitted coefficients over the stepwise regression. These examples and library allow conjectures on future applications and examples of this graph theoretic approach to be explored.

The non-local calculus on finite weighted graphs is laid out in Section \ref{sec:graphtheory}, illustrative physical systems and discussions are presented in Section \ref{sec:physicalsystems}, and conclusions appear in Section \ref{sec:concl}. Appendices outline the implementation details and example calculations for these methods, including details of stepwise linear regression in \cref{app:implementation_stepwiseregression,app:implementation_linearregression}, implementation of the graph theory code library \cref{app:implementation_graphtheorylibrary}, and error analyses of the method in \cref{app:error}.

\section{Graph theoretic approach} \label{sec:graphtheory}
The major formalism used in this work relies on methods from graph theory \cite{Banerjee2019,West2001,Newman2010}. A graph, denoted by $G = (V,E)$, consists of a set of vertices, $i \in V$, with $\abs{V} = n$. The vertices are connected by a set of edges, $e \in E$, where each edge $e = (i,j)$ is a pair of vertices. A vertex is denoted as being the part of an edge by $i \in e$. The graph can also have global and local attributes on the vertices and edges, as we describe below. Of importance are edge weights $\weight[][]$, that are generally functions of the local vertex attributes, and admit a notion of distance between states on the graph as will be elaborated on in \cref{sec:graphtheory_weightconstraints}.

As outlined in \cref{sec:intro_qoi} of the Introduction, in our treatment, a state is a low-dimensional vector of quantities of dimension $p$, obtained after solving a high-fidelity computation of dimension $P$, where ideally $p \ll P$. These high fidelity computations are typically direct numerical simulations (DNS), for example the solution of partial differential equations (PDE) using a finite element-like approach. The DNS simulations can yield some functionals, such as volume averaged fields, or resulting scalar outputs. These, collected into a low-dimensional vector will be denoted by $\func[][]$. Separately, \emph{parameters} such as physical parameters, initial and boundary conditions of the DNS can be collected into another low-dimensional vector denoted by $\var[][] \in \reals[p]$. Each state is represented as a vertex $i \in V$ on the graph. It follows that the notation $\var[][i]$ is unambiguous in labelling a state by its vertex. Transitions between these states, including changes in system parameters, or steps of the numerical solver, allow for connections to be made between vertices on the graph. These transitions therefore can be thought of as inducing edges $e = (i,j)$ between states $i$ and $j$.

Stated formally now, the vertices, $i \in V$ represent states $\{\var[\varindex][i],\func[][](\var[\varindex][i])\}$, where, as introduced above, $\var[][i]$ are parameters (physical parameters, initial and boundary conditions) and $\func[][](\var[][i])$ are quantities of interest, obtained as functionals from the high-dimensional numerical solution of a system of PDEs. The graphs naturally inherit the mathematical structure of the PDEs \cite{West2001,Newman2010}, and even confer new machinery to it \cite{Banerjee2019}.

Several reoccurring graph structures arise in formal graph theory, including fully connected clique graphs, such as the subgraphs shown in \cref{fig:cliquegraph}, or minimally connected tree graphs, such as shown in \cref{fig:treegraph}. There are intricate correspondences between the physical processes being studied and the resulting graph structure. As identified by Banerjee \emph{et al.}~\cite{Banerjee2019}, directed graphs and less densely connected tree graphs correspond to processes with causality or irreversibility induced by energy dissipation, whereas undirected graphs correspond to reversible processes. Fully connected clique graphs generally represent linear processes in equilibrium. These correspondences, centrality measures and graph traversal will be revisited in future communications.

\begin{figure}[ht]
	\centering

	\includegraphics[width=0.6\textwidth,trim={0 0 0 6cm},clip]{figures/clique_graph.pdf}

	\caption{Example graph, with a central clique, and branching clique subgraphs. Vertex coloring may represent various local properties.}
	\label{fig:cliquegraph}
\end{figure}

\begin{figure}[ht]
	\centering

	\includegraphics[width=0.6\textwidth,trim={0 0 0 6cm},clip]{figures/tree_graph.pdf}

	\caption{Example directed tree graphs.}
	\label{fig:treegraph}
\end{figure}

\subsection{Non-local calculus on finite, weighted graphs} \label{sec:graphtheory_nonlocalcalculus}
As outlined in \cref{sec:intro_nonlocalcalculus} of the Introduction, we seek to develop reduced order models in the form of differential equations for the evolution of some component of the state vector on the graph, or as functional representations. In the latter case, we explore forms that are motivated by the Taylor Series of local, differential calculus. In both cases, we need to first define a calculus on the discrete manifold that is the graph. Specifically, we adopt a discrete, non-local calculus on finite, weighted graphs \cite{Hein2007,Gilboa2008,Elmoataz2008,Desquesnes2013,Lozes2014}.

Gilboa \emph{et al.}~ define a discrete calculus, consisting of non-local operators, based on differences between states $i,j \in V$, and an edge weight $\weight[][](\var[][i],\var[][j])$\cite{Gilboa2008}. The edge weights $\weight[][]$ can be defined in a number of ways, and $\weight[][](\var[][i],\var[][j]) = \weight[][](\var[][j] - \var[][i])$ is natural for our purposes. In the case of symmetric, undirected graphs, a radially symmetric weight, $\weight[][](\var[][i],\var[][j]) = \weight[][](\radius[][](\var[][i],\var[][j]))$, where $\radius[][](\var[][i],\var[][j]) = \abs{\var[][j] - \var[][i]}$ is also appropriate. As will be discussed, there are several constraints placed on the allowed forms of the weight functions, to ensure that the non-local calculus converges to true differential operators in certain limits.

Scalars, $\func[][](\var[][i])$ correspond to a function at a single vertex $i$. Vectors $\functwo[](\var[][i],\var[][j])$ are functions on vertex pairs $i,j$ or the edge $e = (i,j)$. Accordingly, we have a vector space $\mathnormal{L}$ such that there is a mapping $V\times V \mapsto \mathnormal{L}$ or $E\mapsto \mathnormal{L}$. Integration on the discrete manifold $G$ is a sum over vertices $i \in V$, and is relevant to operations such as inner products, which we define to be normalized by the size of the graph, $\abs{V} = n$. There should be no cause for confusion between the state vectors $\var[][i] \in \reals[p]$ represented as vertices $i \in V$, and the edge vectors $\functwo[](\var[][i],\var[][j]) \in \mathnormal{L}$. We note that for this treatment, there are no self edges and so vectors $\functwo[](\var[][i],\var[][j])$ at vertex $i$ are defined as being between vertices distinct from $i$, and $j \neq i$. These are distinct manifestations of vector spaces that are not copies of each other.

\subsubsection{Non-local calculus definitions} \label{sec:graphtheory_nonlocalcalculus_definitions}
The non-local gradient operator $\gradient[][{\func[][]}] = \gradient[{\weight[][]}][[{\func[][]}]](\var[][i],\var[][j]) $, with respect to the weight $\weight[][]$, can now be defined as the vector functional of scalars $\func[][]$ at state $\var[][i]$, and represents the vector of weighted differences with \textit{all other vertices} $j \in V$:
\begin{equation}
	 \gradient[{\weight[][]}][[{\func[][]}]](\var[][i],\var[][j]) \equiv [\func[][](\var[][j])-\func[][](\var[][i])]\sqrt{\weight[][](\var[][i],\var[][j])},\quad \forall ~ j \neq i \in V.
\end{equation}
The inner product between scalars $\func[][1]$ and $\func[][2]$ is defined as the sum over the vertices:
\begin{align}
	\langle \func[][1],\func[][2]\rangle \equiv&~ \frac{1}{n}\sum_{i\in V} \func[][1](\var[][i])\func[][2](\var[][i]),
\end{align}
and the contraction between vectors $\functwo[][1], \functwo[][2]$ at vertex $i$ is defined as the sum over disparate vertices:
\begin{align}
	[\dotproduct[{\functwo[][1]}][{\functwo[][2]}]](\var[][i]) \equiv&~ \frac{1}{n-1}\sum_{j \in V\setminus \{i\}} \functwo[][1](\var[][i],\var[][j])\functwo[][2](\var[][i],\var[][j]).
\end{align}
Inner products, and norms between pairs of vectors are likewise defined.
From these definitions, in the case of symmetric weights, there are parallels to standard differential calculus. There is an adjoint relation between the gradient and the divergence operators, as well as a divergence theorem, leading to a definition of the Laplacian being found to be self-adjoint.

\noindent For this work, the most important operators are the partial derivatives of $\func[][](\var[][i])$ with respect to $\var$, the $\ith[\varindex]$ component of the state $\var[][]$. Here, we denote differential calculus partial derivatives with $\partial$, and non-local partial derivatives with $\delta$. The argument of the function in $\func[][](\var[][i])$ also makes clear that the function's partial derivative is evaluated at vertex $i$. These non-local operations can be used to rigorously define the partial derivatives as the contraction between the function gradient and the specifically defined unit vector,
\begin{align}
	\difference[1]{\func[][](\var[][i])}{\var} ~=&~ [\dotproduct[{\gradient[{\weight[][]}][{[{\func[][]}](\var[][i])}]}][{\hat{\var[][]}^{\varindex}_{i}}]] \label{eq:diff{1}} \\
	 ~=&~~ \frac{1}{n-1}\sum_{j \in V\setminus \{i\}} [\func[][](\var[][j]) - \func[][](\var[][i])](\var[\varindex][j] - \var[\varindex][i])\weight[][](\var[][i],\var[][j]),
\end{align}
where the unit vector at state $i$ in the $\varindex$ direction is defined as a vector gradient between states $i$ and $j$:
\begin{equation}
	\hat{\var[][]}^{\varindex}(\var[][i],\var[][j]) \equiv \gradient[{\weight[][]}][[\var]](\var[][i],\var[][j]) = [\var[\varindex][j] - \var[\varindex][i]]\sqrt{\weight[][](\var[][i],\var[][j])}.
\end{equation}
We will proceed with the weights $\weight[][](\var[][i],\var[][j])$ being chosen such that the unit vectors are normalized:
\begin{equation}
 	[\dotproduct[{\hat{\var[][]}^{\varindex}}][{\hat{\var[][]}^{\varindexone}}]] = \kronecker[\varindex][\varindexone], \label{eq:unitnorm}
\end{equation}
where $\kronecker[\varindex][\varindexone]$ is the Kronecker delta. As will be discussed, this constraint together with the inner products also being normalized by the volume of the space, leads to other important behaviors of the convergence and consistency of these non-local derivatives to their counterparts in differential calculus.

Higher order derivatives also can be computed as approximations to $\nderivative[p]{\func[][](\var[][])}{\var[{\varindex[1]}][]}{\var[{\varindex[p]}][]}$ using an extension of \cref{eq:diff{1}} and can be defined recursively:
\begin{align}
\ndifference[q]{\func[][](\var[][i])}{\var[{\varindex[1]}]}{\var[{\varindex[q]}]} ~=&~
[\dotproduct[{\gradient[{\weight[][]}][{[
	{\dotproduct[{\gradient[{\weight[][]}][{[
	{\cdots[\dotproduct[{\gradient[{\weight[][]}][{[{\dotproduct[{\gradient[{\weight[][]}][{[\func[][](\var[][i])]}]}][{\hat{\var[][]}^{\varindex[1]}}] }]}]}][{\hat{\var[][]}^{\varindex[2]}}]]\cdots}]}]}][{\hat{\var[][]}^{\varindex[q-1]}}]}]}]}][{\hat{\var[][]}^{\varindex[q]}}]] \label{eq:diff_p}
\end{align}

It will now be shown, that under specific constraints on the weight functions, the above partial derivatives are, to leading order, the corresponding partial derivatives of differential calculus. However certain properties of differential calculus, such as the commutativity of mixed partial derivatives, do not automatically hold to all orders of approximation with all possible choices of edge weights. Fortunately, with weights that decay suitably, issues of commutativity and convergence can be guaranteed.

\subsection{Constraints on edge weights} \label{sec:graphtheory_weightconstraints}
When using the given non-local calculus definitions, there is an arbitrary normalization of the quantities that can be chosen. For numerical purposes, and to ensure the operator definitions are consistent and independent of the size of the graph, as discussed, all integrations over the graph are normalized by the size of the discrete space we are in $\vert V\vert = n$. In what follows repeated Greek indices $\varindex \in \{1,\dots p\}$ will follow Einstein summation convention.

For the non-local calculus, we will observe order-by-order how the partial derivatives converge, given derived constraints on the weight function $\weight[][](\var[][i],\var[][j]) = \weight[][](\var[][i],\var[][j])$ for the $p$ dimensional representation of the state vector $\var[][i]$. We begin by examining the normalization condition in \cref{eq:unitnorm} for this graph, with the normalization constant denoted as $\Weight[][][p,n]$
\begin{align}
	\sum_{j \in V\setminus \{i\}} (\var[\varindex][j]-\var[\varindex][i])(\var[{\varindexone}][j]-\var[{\varindexone}][i])\weight[][](\var[][i],\var[][j]) = \Weight[\varindex\varindexone][][p,n] = (n-1)\kronecker[\varindex][\varindexone] ~~~ \forall~ i. \label{eq:weightnorm}
\end{align}
Due to the quadratic form of the constraint, posed as a second order moment of this weight function, it is natural to assume that the weights are decaying, and take a form that decays at least as rapidly the inverse square of the distance between graph vertices:
\begin{align}
	\weight[][](\var[][i],\var[][j]) = O\left(\frac{1}{\abs{\var[][j]-\var[][i]}^2}\right).
\end{align}
For $p=1$, the form of the weight may be chosen as precisely the relation above, and for $\var[][] = \var[1][]$,
\begin{align}
	\weight[][](\var[][i],\var[][j]) = \frac{1}{\abs{\var[][j]-\var[][i]}^2}.
\end{align}
For $p>1$ dimensional graphs, the form of the weight is less obvious if the Kronecker delta representation is to be extended or generalized. The weight function must decay rapidly enough to minimize consistency error and ensure that off-diagonal terms on the left-hand side of \cref{eq:weightnorm} are zero; however, too rapid a decay may fail to satisfy the finite diagonal constraints.

In the case of radially symmetric metrics $\weight[][]({\radius[][ij]})$, where ${\radius[][ij]}^2 = \sum_{\varindex} \abs{\var[\varindex][j]-\var[\varindex][i]}^2$, the normalization constraint can be written as
\begin{align}
	\sum_{j \in V\setminus \{i\}} {\radius[][ij]}^2 \weight[][]({\radius[][ij]}) = \Weight[][][p,n] = (n-1)p ~~~ \forall~ i,
	\label{eq:discreteweightnormalization}
\end{align}
which follows by computing the trace in \cref{eq:weightnorm}.

In generalizing \cref{eq:weightnorm} to the case of continuous manifolds, we consider the case of graphs whose states $\var[][i]$ have support over a subset of $\reals[p]$, where the discrete sums approach integrals. Here we desire a form for the weights to guarantee that all off-diagonal terms, $\varindex \neq \varindexone$ are identically zero. In this setting, \cref{eq:discreteweightnormalization} leads to the form
\begin{align}
	\int \integrand[]{\radius} \radius[p+1][] \weight[][](\radius[][]) = \Weight[][][p,R] = p\frac{V_p}{S_p}\Radius^p = \Radius^p, \label{eq:weightrad}
\end{align}
where $V_p$ and $S_p$ are the volume and surface area of the unit ball in $p$ dimensions, $S_p = pV_p$, and $\Radius$ is the radius of the sphere that encloses the graph.

In this continuous case of a graph whose states $\var[][i]$ are dense in $\mathbb{R}^p$, the off-diagonal terms in \cref{eq:weightnorm} can be shown to cancel due to symmetry considerations, since $\weight[][](\radius)$ is an even symmetric function. Letting $\var[\varindex][] = rf^{\varindex}(\theta)$, where $f$ are angular functions, the spherical symmetry ensures that only diagonal terms are non-zero in the tensor components
\begin{align}
	\Weight[][\varindex\varindexone][p,R] &\equiv \int \integrand[p]{\varone[][]} (\varone[\varindex] - \var[\varindex][])(\varone[\varindexone][] - \var[{\varindex}])\weight[][](\var[][],\varone[][]).
\end{align}
For example in $p=2$ dimensions, where $\var[1][] = \radius\cos{\theta}, \var[2][] = \radius\sin{\theta}$, the angular integral trivially evaluates to zero.

However, in the case of graphs that are not dense in $\reals[p]$, with radius
\begin{align}
	\Radius[][] = \max_{ij} \frac{{\radius[][ij]}}{2}
\end{align}
in order for these off-diagonal terms to cancel in these finite integrals for an arbitrary graph structure, in addition to being symmetric functions, the weight functions must also decay rapidly so that the odd, antisymmetric terms are negligible.

\noindent It will now be assumed that for $p>1$, the weights decay faster than the above assumed inverse square relationship:
\begin{align}
	\weight[][](\radius) = o\left(\frac{1}{\radius^2}\right).
\end{align}
There still remains some choice in the exact form of the weights, and for this work, we will consider radially symmetric weights that decay like a Gaussian function $\sim e^{-r^2}$, or decay polynomially $\sim 1/r^2$. Based on these forms of the weights, there are some free parameters that can be calculated or chosen. We choose the weights to have the form $\weight[][](\radius) = C ~\tilde{\weight[][]}(\radius/\sigma)$, where $C(p,\Radius)$ is the weight scale, and $\sigma(p,\Radius)$ is the weight decay parameter, or an appropriate length scale in the graph manifold. Based on \cref{eq:discreteweightnormalization}, the weights must scale as the inverse of length-squared, and so we will discuss the weight properties in terms of the scaled weight constraints
\begin{align}
\widetilde{\Weight[][]}[p,z] = \frac{\Weight[][][p,\Radius]}{C\sigma^{p+2}}, \label{eq:scaled_weight_constraint}
\end{align}
where $z = {R}/{\sigma}$ is a scaled radius. The weight constraints therefore take the form of $p$-order moments of the weight function:
\begin{align}
	\int_0^z \integrand[]{\varone[][]} {\varone[][]}^{p+1} \tilde{\weight[][]}(\varone[][]) = \widetilde{\Weight[][]}[p,z]. \label{eq:constraint_scaled}
\end{align}

\subsubsection{Gaussian weights}\label{sec:graphtheory_weightconstraints_gaussian}
\noindent We first may consider Gaussian decaying weights
\begin{align}
	\weight[][](\radius) = Ce^{-\frac{1}{2}\frac{\radius^2}{\sigma^2}}, \label{eq:weight_gauss}
\end{align}
where in the case of graphs with a finite radius $\Radius$, the weight is a truncated Gaussian for $0\leq \radius\leq \Radius$. The scaled weight constraints are the $p+1$ order moments of the Gaussian distribution. These moments have a known recursion relationship \cite{Orjebin2014} for $p>0$:
\begin{align}
	\widetilde{\Weight[][]}[p,z] ~=&~ p\widetilde{\Weight[][]}[p-2,z] - z^pe^{-\frac{1}{2}z^2}, \label{eq:recurrence_gauss}
\end{align}
which also can be expressed in the analytic form
\begin{align}
	\widetilde{\Weight[][]}[p,z] ~=&~ p!!\left[\widetilde{\Weight[][]}[\nmodm{p}{2},z] - z^{p}e^{-\frac{1}{2}z^2}\sum_{i=0}^{\floor{p/2}-1} \frac{z^{-2i}}{(p-2i)!!}\right].
\end{align}

\subsubsection{Polynomial weights}\label{sec:graphtheory_weightconstraints_polynomial}
We also may consider polynomially decaying weights
\begin{align}
	\weight[][](\radius) = C\left(\frac{\sigma}{\radius}\right)^{2+\epsilon}, \label{eq:weight_poly}
\end{align}
where $0\geq\epsilon<p$ is the super-quadratic decay scaling. As previously discussed, the weights must decay faster than quadratically, however cannot decay so fast as for the integrals to vanish over discrete spherical volumes. Here the scaled weight constraints are tractable integrals for $p>0$:
\begin{align}
	\widetilde{\Weight[][]}[p,z] = \frac{1}{p-\epsilon} z^{p-\epsilon}.	\label{eq:constraint_poly}
\end{align}
Values for $\widetilde{\Weight[][]}[p,z]$ for radially symmetric Gaussian and polynomially decaying weights, with a finite radius, continuous graph are shown in \cref{tab:weights}. Here, $\Phi(x)$ is the standard normal cumulative distribution function.
\begin{table}[htb!]
	\centering
	\begin{tabular}{|c|c|c|} \hline\tableheight{20pt}
		$p$ & Gaussian & Polynomial \\ \hline\tableheight{20pt}
		$-2$ & $0$ & Undefined \\ \hline\tableheight{20pt}
		$-1$ & $\sqrt{2\pi}\left[\Phi(z) - \Phi(0)\right]$ & Undefined\\ \hline\tableheight{20pt}
		$0$ & $1 - e^{-\frac{1}{2}z^2}$ & Undefined \\ \hline\tableheight{20pt}
		$1$ & $\sqrt{2\pi}[\Phi(z) - \Phi(0)] - ze^{-\frac{1}{2}z^2}$ & $\frac{1}{1-\epsilon} z^{1-\epsilon}$ \\ \hline\tableheight{20pt}
		$2$ & $2 - (2+z^2)e^{-\frac{1}{2}z^2}$ & $\frac{1}{2-\epsilon} z^{2-\epsilon}$ \\ \hline

	\end{tabular}
	\caption{Weight constraint values for a radially symmetric weight $\widetilde{\Weight[][]}[p,z]$, given a finite radius $z = {\Radius}/{\sigma}$, continuous graph in $p$ dimensions.}
	\label{tab:weights}
\end{table}

\subsubsection{Free parameters in weight functions}\label{sec:graphtheory_weightconstraints_freeparameters}

The free parameters $C$ and $\sigma$, and other hyperparameters such as the scaling $\epsilon$, may now be constrained based on for example \cref{eq:constraint_scaled} and \cref{eq:constraint_poly}, and there is some freedom, depending on the graph of interest. Generally, the decay parameter of the weight $\sigma$ can be chosen, after fixing the scale
\begin{align}
 C ~=&~ \frac{z^p}{\sigma^2}\frac{1}{\widetilde{\Weight[][]}[p,z]}, \label{eq:weightscale}
\end{align}
where the scaling of the weight by $1/\sigma^2$ is evident.

\begin{figure}[hpt]
\centering
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=0.6\textwidth]{figures/Wtilde_p_sigma_gauss.pdf}

	\subcaption{Gaussian}
	\label{fig:Wtilde_p_sigma_gauss}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=0.6\textwidth]{figures/Wtilde_p_sigma_poly.pdf}

	\subcaption{Polynomial}
	\label{fig:Wtilde_p_sigma_poly}
\end{subfigure}
\caption{Scaled weight constraints for various $p$ values, as a function of weight decay parameter $\sigma$; $\epsilon = {p}/{2}$ for the polynomial weights.}
\label{fig:Wtilde_p_sigma}
\end{figure}

\begin{figure}[hpt]
\centering
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=0.6\textwidth]{figures/C_p_sigma_gauss.pdf}

	\subcaption{Gaussian}
	\label{fig:C_p_sigma_gauss}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=0.6\textwidth]{figures/C_p_sigma_poly.pdf}

	\subcaption{Polynomial}
	\label{fig:C_p_sigma_poly}
\end{subfigure}
\caption{Weight constraint scales for various $p$ values, as a function of weight decay parameter $\sigma$; $\epsilon = {p}/{2}$ for the polynomial weights.}
\label{fig:C_p_sigma}
\end{figure}

\begin{figure}[hpt]
\centering
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[trim={0 2 0 0},clip,width=0.8\textwidth]{figures/w_p_sigma_gauss.pdf}
	\vspace{-1.5cm}
	\subcaption{Gaussian}
	\label{fig:w_p_sigma_gauss}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[trim={0 -10pt 0 0},clip,width=0.8\textwidth]{figures/w_p_sigma_poly.pdf}
	\vspace{-1.5cm}
	\subcaption{Polynomial}
	\label{fig:w_p_sigma_poly}
\end{subfigure}
\vspace{-0.5cm}
\caption{Weight functions for various $p$ values, as a function of weight decay parameter $\sigma$ and radius $r$; $\epsilon = {p}/{2}$ for the polynomial weights.}
\label{fig:w_p_sigma}
\end{figure}
Using Gaussian and polynomially decaying weights, the calculated weight constraints $\widetilde{\Weight[][]}[p,\sigma]$, and the weight scalings $C[p,\sigma]$ are shown in \cref{fig:Wtilde_p_sigma} and \cref{fig:C_p_sigma}. For the polynomial weights, $\epsilon = {p}/{2}$. Here it can be seen that the Gaussian weight constraints are constant for small $\sigma > 0$, and then decrease for $\sigma > {R}/{z_{\textrm{inflection}}}$, where $z_{\textrm{inflection}}^2 = p+{3}/{2} + \sqrt{2p+{9}/{2}}$ is the point at which the Gaussian case of the integrand in \cref{eq:constraint_scaled} changes curvature and tends to zero. The Gaussian weight scalings also decrease polynomially in $\sigma$, as per \cref{eq:weightscale}. The polynomial weight constraints decay as ${1}/{\sigma^{p-\epsilon}}$, and the polynomial weight scalings decay as ${1}/{\sigma^{2+\epsilon}}$.

Importantly, for a fixed $\epsilon$, the polynomial weight scalings scale much slower with $\sigma$ and $p$ compared to the Gaussian weight scalings due to this polynomial form. For $\epsilon=0$ and a given $p$ value, how the polynomial scalings $C$ vary with $\sigma$ exactly cancel how the weight itself varies with $\sigma$. Therefore the polynomial weight function is constant with respect to $\sigma$:
\begin{align}
	\weight[][](r) = (p-\epsilon)\frac{R^{\epsilon}}{r^{2+\epsilon}}, \label{eq:poly_weight_full}
\end{align}
and suggests this weight form is more easily applied to graphs with different weight scales. In contrast, the Gaussian weights scale as ${1}/{\sigma^{2}} e^{{R^2}/{2\sigma^2}}$ with respect to $\sigma$ and increases exponentially for $\sigma \ll R$. Both of these calculations show the weight values also increasing with $p$, indicating that functions that depend on a large number of variables may still have very large weights, regardless of the $\sigma$ dependence.

We note that the Gaussian weights bear connections to the diffusion kernel, which furnishes the probability distribution between states on a graph in the context of discovering the data embedding in a continuous manifold. Currently the form of the weights also has a single length scale $\sigma$, and the manifold is assumed to be homogeneous along each state variable dimension. Future communications will investigate the effects of using each independent equation along each dimension in \cref{eq:discreteweightnormalization} to determine constraints on length scales in each direction $\sigma^{\varindex}$. For the case of discrete manifolds, certain approximations on these continuous weight constraints must be made. Please refer to \cref{app:implementation_gauss} of the Appendix for implementation details of these approximations.

\subsection{Consistency of non-local derivatives} \label{sec:graphtheory_consistency}
With this insight into the form of the weights, and their constraints, we now examine the consistency of the partial derivatives in this non-local calculus. The analysis will follow a similar approach to any numerical scheme, substituting Taylor series of differential calculus terms for various non-local terms, to determine the leading order behavior. Due to the weight constraint, we will see that depending on the decay rate of the weights, the leading order error is first order in the spacing of the vertices.

\noindent For the first order derivatives the weight constraint yields
\begin{align}
	\difference[1]{\func[][](\var[][i])}{\var}
	~=&~ \derivative[1]{\func[][](\var[][i])}{\var[\varindex][]} + O\left(\derivative[2]{\func[][](\var[][i])}{\var[\varindexone],\var[\varindextwo][]}(\var[\varindexone][j]-\var[\varindexone][i])(\var[\varindextwo][j]-\var[\varindextwo][i])(\var[\varindex][j]-\var[\varindex][i]) \weight[][](\var[][i],\var[][j])\right). \label{eq:non-local1}
\end{align}
For second order derivatives, each non-local first derivative can be expanded in the form found above in \cref{eq:non-local1}, yielding
\begin{align}
	\difference[2]{\func[][](\var[][i])}{\var,\var[\varindexone]}
	~=&~ \derivative[2]{\func[][](\var[][i])}{\var,\var[\varindexone]} + O\left(\derivative[2]{\func[][](\var[][i])}{\var[\varindextwo][],\var[\varindexthree]} (\var[\varindextwo][j]-\var[\varindextwo][k])(\var[\varindexthree][j]-\var[\varindexthree][k])(\var[\varindex][j]-\var[\varindex][k])(\var[\varindexone][k]-\var[\varindexone][i]) \weight[][](\var[][k],\var[][j])\weight[][](\var[][i],\var[][k])\right).\label{eq:non-local2}
\end{align}
Therefore, given this weight constraint, the first order non-local derivatives to leading order correspond with the corresponding local derivatives. Although the second order non-local derivatives do not contain first order terms, the second order terms contain the other second order derivatives along the other dimensions. It can be concluded that in order for the non-local calculus definitions for the partial derivatives to equal the local calculus definitions to leading order, a constraint on the higher order moments of the weight function for each order of derivative is necessary. These constraints are assumed to be of the form of products of all differences in state variables at all combinations of dimension and pairs of vertices, weighted by the edge weight between those vertices:
\begin{align}
\ndifference[q]{\func[][](\var[][i])}{\var[{\varindex[1]}]}{\var[{\varindex[q]}]}  = \nderivative[q]{\func[][](\var[][i])}{\var[{\varindex[1]}]}{\var[{\varindex[q]}]} + O((\var[][j]-\var[][i])) &\longleftrightarrow \\
\sum_{j_{1}\dots j_{q} \in V\setminus \{i\}} (\var[{\varindex[q]}][j_{q}]-\var[{\varindex[q]}][{j_{q-1}}]) \cdots (\var[{\varindex[1]}][j_{1}]-\var[{\varindex[1]}][{i}]) \weight[][](\var[][j_{q}],\var[][{j_{q-1}}]) \cdots \weight[][](\var[][{j_{1}}],\var[][{{i}}]) ~=&~ (n-1)^q\kronecker[{\varindex[q]}][{\varindex[q-1]}]\cdots \kronecker[{\varindex[1]}][{\varindex[i]}]. \nonumber
\end{align}

As alluded to in the previous example of decaying Gaussian and polynomial weights, an alternative to choosing weights based on these difficult to construct higher order moment constraints, is to choose weights that decay faster than any finite polynomial, and still satisfy the primary constraint in \cref{eq:weightnorm}. Given a fast enough decaying weight, this will allow several issues to be resolved. First, all off-diagonal terms in the primary constraint will decay rapidly. Second, all higher order non-local derivatives will equal their differential derivative counterparts to leading order due to the weights decaying faster than the polynomials in the Taylor series.

\subsubsection{Compromises with a non-local calculus approach}\label{sec:graphtheory_consistency_compromise}

However, due to this calculus involving sums over the entire graph, the non-local nature means the terms beyond the local calculus derivatives in the consistency relationships above in \cref{eq:non-local1,eq:non-local2} may be unreasonably large. This is particularly true if the weight scales are large as we have seen for the Gaussian weights. Therefore there may be no real similarity in the values and trends of the non-local derivatives compared to their differential counterparts. However, from studies of the total scale of the weights, due to the polynomially decaying weights having a closed form that is independent of the decay parameter, these weights are reasonable scales of order $O(1)$ for small $p$.

Due to these numerical issues, we choose to use the polynomially decaying weights for the current studies, with $\epsilon = {p}/{2}$. Although not necessarily similar in form, the non-local derivative values will be similar in magnitude to what would be expected for differential derivatives, and decay reasonably quickly with $\epsilon>0$. Future communications will study in depth the effect of using different forms and parameter choices for the weights, and how these large non-local error terms can be understood and controlled.

This extends the compromise regarding the weight functions to be made, between obeying the consistency and constraint equations, and having scales of derivatives that correspond with those from differential calculus in the model. As to be discussed when these methods are applied to physical systems, the model fitting procedure may be able to overcome the large scales of these derivatives, such as associating small or large linear coefficients to terms in the model that contain these large values. Please refer to \cref{app:implementation_linearregression} of the Appendix for details on the regression implementation.

This non-local calculus has the objective of being a general, accurate framework for any modelling procedure, and the exact form of the non-local error term will be developed. An important aspect of numerical method development is an analysis of the effect of data sparsity on model error. As can be seen in \cref{eq:non-local1}, the non-local derivatives have a persistent error term for all weight functions with non-local extent, and whether this is detrimental to subsequent reduced order modelling is an important question in determining the validity of this approach.

For example, a modified Taylor series, as discussed in \cref{app:implementation_modifiedtaylorseries} of the Appendix, may be constructed using a training dataset as a functional representation, where a different Taylor series model is developed at each possible base point. Please refer to \cref{app:error_mesh} for details concerning the mesh of data used. The point-wise error for evaluation points within a distance $h$ of each base point model can then be computed, and refining the spacing of the data, we are able to get a sense of how the method error scales with $h$. As discussed in \cref{app:error_pointwisederivative} of the Appendix, in addition to the expected effects of a $\ith[k]$ order Taylor series having error that scales as $h^{k+1}$, the non-local calculus derivatives, and any fit model parameters are also dependent on the spatial distribution of the training and testing data. Initial error estimates for this modified Taylor series approach can be referred to in \cref{app:error_totalerror}. Future communications will complete a full study of these errors, to be compared to standard interpolation and other functional representations.

With this understanding of the behaviors of the non-local calculus and partial derivatives, reduced order models can now be developed for systems of interest. The modelling will consist of finding a basis of operators for a given model, where the operators are defined in terms of this non-local calculus. In the following section, an example of a physical systems is presented, and several reduced order modelling approaches are examined.

\section{Physical systems of interest} \label{sec:physicalsystems}
\subsection{Microstructures in a gradient-regularized model of non-convex elasticity} \label{sec:physicalsystems_microstructures}
We apply the graph theoretic approach to studying reduced order models for the mechanochemical response of solids that form microstructures. This physical system is driven by a free energy density function that is parameterized by composition (chemical) and strain (elastic) variables. The underlying functional form $\energydensity[][{\textrm{hom}}]$ is non-convex in both these quantities and is regularized by composition and strain gradients. Microstructures develop as phases and symmetry-breaking structural variants arise corresponding to negative eigenvalues of the Hessian of the free energy density in the strain-composition space. These models are ill-posed, however; a condition arising from the absence of penalization on inter-phase and inter-variant interfaces and manifesting in PDE solutions that are pathologically discretization-dependent. Well-posedness and penalization of interfaces are restored by free energy functionals $\energydensity[][{\textrm{tot}}] = \energydensity[][{\textrm{hom}}] + \energydensity[][{\textrm{grad}}]$, dependent on the composition and strain gradients. When decoupled, the chemical component of this problem causes spinodal decomposition and Ostwald ripening, described by the Cahn-Hilliard equation \cite{Cahn1958}, and the non-linear elasticity problem by a non-convex strain energy density stabilized by incorporating strain gradient effects as described by many authors.

Our treatment follows Toupin \cite{Toupin1962} and has appeared as decoupled, gradient-regularization of non-convex, non-linear elasticity \cite{Rudraraju2014,Teichertetal2017}, as well as mechanochemical spinodal decomposition \cite{Rudraraju2016,Sagiyamaetal2016,SagiyamaGarikipati2017b}. The free energy density functional is $\energydensity[][]= \energydensity[][{\textrm{tot}}]$, with a dependence on the fields $\{\composition[][], \Estrainspatial, \gradient[]{\composition[][]},\gradient[]{\Estrainspatial[][]}\}$ where $\composition[][]$ is the composition and $\Estrainspatial[][]$ is the Green-Lagrange strain tensor. This dependence, however, is expressed in terms of symmetry-coding strain parameters that are introduced below.

At low temperatures, these systems consist of stable phases related to the local crystalline symmetry, specifically cubic, or tetragonal symmetry. The system considered is restricted to $d=2$ dimensions, and so the material will consist of phases with either square or rectangular structure. This pattern of regions of phases of a field are generally referred to as microstructures. This system is further modeled as a binary mixture, with scalar composition $0\leq \composition[][] \leq 1$ of one of the two constituents.

Additionally, the introduction of mechanical strain couplings to the composition allows for a degeneracy in the rectangular phase. Linear combinations of the strain components $\Estrainspatial$ yield relevant parameters $\estrainspatial$ relating to the original high symmetry of the square phase, namely $\estrainspatialindex[][1,2] \equiv {\Estrainspatialindex[][11] \pm \Estrainspatialindex[][22]}/{\sqrt{2}}, \estrainspatialindex[][6] \equiv \sqrt{2} \Estrainspatialindex[][12]$, that allow this degeneracy to be determined. Specifically, the sign of the $\estrainspatialindex[][2]$ field can determine the orientation of the rectangular phase, either in the horizontal "plus", or the vertical "minus" variant. The order parameters will thus be defined to be $\orderparameter[][] = (\composition[][],\estrainspatialindex[][2])$, and the microstructure can be described by its phase volume fractions $\volfrac[][\vardim]$, where $\vardim = \{\square,\hrectangle,\vrectangle\}$ corresponds to the phase of the material.

\subsubsection{Free energy description} \label{sec:physicalsystems_microstructures_freeenergy}
To describe these phase transitions, we take the description used by Zhang \& Garikipati\cite{Zhang2020}. The total free energy density description of the system $\energydensity[][]$ can be comprised of a homogeneous component, $\energydensity[][{\textrm{hom}}](\composition[][],\estrainspatial)$, that is not entirely concave up over the entire $\composition[][]-\estrainspatialindex[][2]$ domain, yielding potential wells with stable phases. An example two dimensional energy landscape from Rudraraju \emph{et al.}~ \cite{Rudraraju2016} is shown in \cref{fig:energyphase}.
\begin{figure}[ht]
	\centering

	\includegraphics[width=0.6\textwidth]{energy_curvature_mechanochem.pdf}

	\caption{Free energy density landscape over the composition and strain component $\estrainspatialindex[][2]$ order parameters, showing changes in curvature and resulting phases of square and rectangular lattice symmetries, from Rudraraju \emph{et al.}~ \cite{Rudraraju2016}.}
	\label{fig:energyphase}
\end{figure}

During the mechanical-chemical coupled dynamical processes, within the domain where the free energy is concave, such as around the saddle point $B$ in \cref{fig:energyphase}, phase separation will tend to occur, where a mixture of the stable phases, will yield a lower energetic state. This can be seen from a typical lever law analysis for a binary mixture, to obtain the phase volume fractions $\volfrac[][\vardim]$. This simple phase separation within the non-convex region of the free energy density is mechano-chemical spinodal decomposition, where phases tend to separate. Furthermore, the gradient regularization that represents the energy of interfaces also manifests in a number $\numberphase[][\vardim]$ of each variant, and total interfacial lengths $\interfacial[][\vardim]$.

With the composition $\composition[][]$ and symmetry-coding strain parameter $\estrainspatial$ and their gradients, a non-uniform free energy density component $\energydensity[][{\textrm{grad}}](\composition[][],\estrainspatial,\gradient[]{\composition[][]}, \gradient[]{\estrainspatial})$ allows for gradient regularization, typically through (coupled) quadratic terms in the order parameter gradients. These terms yield phase regions that are limited by interfaces being less energetically favorable, and so equilibrium states will tend to have less numbers of domains of phases, and furthermore, larger phase regions grow at the expense of smaller ones. This effect is known as Oswald ripening \cite{Rudraraju2016}.

This gradient contribution, in addition to the homogeneous contribution yields the total free energy functional of the order parameters over the material volume $\Omega_t$ at time $t$,
\begin{equation}
	\energy[][][\composition[][],\estrainspatial] = \int_{\Omega_t} \integrand[]{V} (\energydensity[][{\textrm{hom}}]+\energydensity[][{\textrm{grad}}]) + \int_{\partial \Omega_t} \integrand[]{S}\ \dotproduct[{\traction[][]}][{\udeform[][]}], \label{eq:Psi_tot}
\end{equation}
where $\traction[][]$ are mechanical traction terms, depending on which boundary conditions are specified and $\udeform[][]$ is the spatial displacement field of the material.

A variational formalism $\difference[1]{\energy[][]}{{\composition[][],\estrainspatial}}$ can then be applied to yield the relevant generalized chemical potentials. As discussed, for this particular work, the homogeneous free energy density is a smooth polynomial function of the order parameters,
\begin{equation}
	\energydensity[][{\textrm{hom}}](\composition[][],\estrainspatial) \equiv 16 \const[][{\composition[][]}] \composition[][]^4 - 32\const[][{\composition[][]}] \composition[][]^3 + 16\const[][{\composition[][]}] \composition[][]^2 + 2\frac{\const[][\estrainspatialindex]}{{\constone[][\estrainspatialindex]}^2}({\estrainspatialindex[][1]}^2 + {\estrainspatialindex[][6]}^2) + \frac{\const[][\estrainspatialindex]}{{\constone[][\estrainspatialindex]}^2}{\estrainspatialindex[][2]}^4 + 2\frac{\const[][\estrainspatialindex]}{{\constone[][\estrainspatialindex]}^2}(2c-1){\estrainspatialindex[][2]}^2,
	\label{eq:psi_hom}
\end{equation}
where $\{\const[][{\composition[][]}],\const[][\estrainspatialindex],\constone[][\estrainspatialindex]\}$ are constant coefficients.

The gradient contributions will take the form of strictly quadratic forms in the order parameter gradients,
\begin{equation}
	\energydensity[][{\textrm{grad}}](\gradient[]{\composition[][]},\gradient[]{\estrainspatial}) = \frac{1}{2} \gradient[]{\composition[][]} \cdot\vec{\consttwo[][]} \gradient[]{\composition[][]} + \frac{1}{2} \gradient[]{\estrainspatial} \cdot \boldsymbol{\constthree[][]} \gradient[]{\estrainspatial} + \gradient[]{\composition[][]}\cdot \boldsymbol{\constfour} \gradient[]{\estrainspatial}, \label{eq:psi_grad}
\end{equation}
and for this work, it is assumed that there are no gradient couplings between the composition and strain fields $\boldsymbol{\constfour} = 0$, and only constant isotropic tensors $\vec{\consttwo[][]} = \consttwo[][]\identity$ and $\boldsymbol{\constthree[][]} = \constthree[][]\identity$ for the gradient terms. Furthermore, the additive terms mean the total free energy density can be separated into the sum of independent mechanical and chemical terms, plus a coupling energy term:
\begin{equation}
	\energydensity[][]= \energydensity[][{\textrm{mech}}] + \energydensity[][{\textrm{chem}}] + \energydensity[][{\textrm{couple}}].
\end{equation}

When stationarity of the free energy is imposed in a variational framework, the Euler-Lagrange equations yield the weak form of the balance of linear momentum, and further variational arguments lead to the strong form. Using this variational approach, based on this form of the free energy, and given its dependence on the strain and composition fields, the governing PDE's are chosen to be coupled Cahn-Hilliard dynamics for the composition \cite{Cahn1958}, and non-linear gradient elasticity equations in equilibrium \cite{Toupin1962,Rudraraju2016}.

\subsubsection{Cahn-Hilliard dynamics} \label{sec:physicalsystems_microstructures_CahnHilliard}
Cahn-Hilliard dynamics are first order for the composition and take the form of a transport equation
\begin{align}
	\derivative[1]{\composition[][]}{t} + \divergence[][{\current[][]}] ~=&~ 0, \label{eq:pde_cahnhilliard}
\end{align}
where
\begin{align}
	\current[][] ~=&~ -\mobility[][]\gradient[][\chemicalpotential], \label{eq:cflux}
\end{align}
and $\mobility[][] = \mobility[][](\composition[][],\estrainspatial)$ is a mobility transport tensor, and is assumed to be isotropic $\mobility[][] = \mobilityscalar\Identity$. The variational treatment yields the chemical potential $\chemicalpotential = \difference[1]{\energy[][]}{\composition[][]}$, whose gradient yields the flux in a form that is guided by the thermodynamic dissipation inequality \cite{degrootmazur1984}, and leads to there being fourth-order mechano-chemically coupled equations.

\subsubsection{Toupin model of elasticity} \label{sec:physicalsystems_microstructures_Toupin}
The strain gradient elasticity can be modeled in equilibrium as equations for the first Piola-Kirchhoff $\Pstressspatial$ and higher order $\Bstressspatial$ stress tensors,
\begin{align}
	\Pstressspatial ~=&~~ \derivative[1]{\energydensity[][]}{\fgrad[][]} \label{eq:piola1}\\
	\Bstressspatial ~=&~~ \derivative[1]{\energydensity[][]}{\gradient[]{\fgrad[][]}} \label{eq:piola3}
\end{align}
which are derivatives of the strain energy and are conjugate to the deformation gradient $\fgrad[][]$. Here the Green-Lagrange strains are $\Estrainspatial = {1}/{2}\left(\fgrad[T][]\fgrad[][] - \Identity\right)$ and gradients $\gradient[]$ are derivatives with respect to the material indices. The elasticity governing equations, shown here in the strong form, are
\begin{align}
	\divergence[]{\Pstressspatial} - \divergence[]{(\divergence[]{\Bstressspatial})} = 0, \label{eq:pde_elasticity}
\end{align}
plus additional Dirichlet and higher order traction boundary conditions.

\subsubsection{Direct numerical simulations} \label{sec:physicalsystems_microstructures_DNS}
These governing PDEs in \cref{eq:pde_cahnhilliard,eq:pde_elasticity} may now be solved, yielding time series data, as per Zhang \& Garikipati \cite{Zhang2020}. The weak form is chosen here for implementation, motivated by numerical methods such as finite elements and the isogeometric approaches. \cite{Rudraraju2014,Rudraraju2016} A given system is described by the static mechanical Dirichlet boundary conditions for the displacement fields $\udeform[][]$, uniform initial $\estrainspatialindex[][2] = 0$ fields, corresponding to the material being initially in the high temperature square phase, random initial compositions $\composition[][] = 0.46 \pm 0.05$, and zero chemical flux boundaries. The DNS proceeds by time stepping with steps of $\Delta t = 10^{-5}~\textrm{s}$, and iteratively solving the PDE's at each time step, with tolerances placed on the residual of the solver. The data generated for this work is at times $50-1500 ~\mu\textrm{s}$, with example DNS observables shown in \cref{fig:DNSobservables}. In this work, we consider a single boundary and initial condition and one evolution of the microstructure and aim to model its quantities of interest. In future communications we will consider the effectiveness of the graph theoretic approach at producing general models for a family of microstructures with different boundary and initial conditions.

To develop the graph theoretic model for the given mechano-chemical system, the state vector is assembled from quantities of interest that are functionals of the high-dimensional solution to \cref{eq:pde_cahnhilliard,eq:pde_elasticity}. These quantities, denoted as barred variables unless otherwise specified, are effective volume averaged quantities $\Functional[][]$ over the material domain. They therefore contain no spatial dependence of their high dimensional counterparts $\functional[](\orderparameter[][](x))$. Chemical quantities such as the local composition and chemical potential are represented by the phase volume fractions, interfacial lengths, and numbers of phase domains
\begin{align}
\{\composition[][](x),~\chemicalpotential(x)\} ~\to&~ \{\volfrac[][\vardim],~\interfacial[][\vardim],~\numberphase[][\vardim] \}_{\vardim =\left\{\square,\hrectangle,\vrectangle,\dots\right\}}; \nonumber \\
\intertext{mechanical quantities such as the deformation gradients, Green-Langrange strains, and Piola-Kirchoff stresses are represented by their volume averaged quantities}
\{\fgrad[][](x),\Estrainspatial[][](x),~\Pstressspatial[][](x),\dots\} ~\to&~ \{\Fgrad[][],\Estrain[][],~\Pstress[][],\dots\}; \nonumber \\
\intertext{and free energy quantities are represented by their total values}
\{\energydensity[][](x),~\energydensity[][{\textrm{mech}}](x),~\energydensity[][{\textrm{chem}}](x),\dots\} ~\to&~ \{\energy[][],~\energydensity[][{\textrm{mech}}],~\energydensity[][{\textrm{chem}}],\dots \}.
\end{align}
The phase volume fractions $\volfrac$ are calculated by summing over volume elements that are within the associated phase, as per the phase diagram in \cref{fig:energyphase}, the interfacial lengths $\interfacial[][\vardim]$ are calculated by computing the total lengths of the interfaces between the square-rectangular-(plus,minus) fields, the number of domains $\numberphase[][\vardim]$ are counted through a gradient edge detection algorithm to identify the number of distinct domains of each field. The volume like terms of $\interfacial[][\vardim]$ and $\numberphase[][\vardim]$ are essential for distinguishing changes in topology of the microstructure, such as when field domains merge, separate, appear, or disappear. The specific forms of the free energy $\energy[][], \energydensity[][{\textrm{mech}}]$ are the integrated free energy densities over the material volume. For the strain quantities, they are chosen to be defined in terms of the volume averaged Green-Lagrange strains
\begin{align}
	\Estrain = \frac{1}{2}\left(\Fgrad[T][]\Fgrad[][] - \Identity\right),
\end{align}
where it is the deformation gradient $\fgrad[][]$ that has been explicitly volume averaged
\begin{align}
\Fgrad[][] = \int_{\Omega_t} \integrand[]{V} \fgrad[][].
\end{align}

Due to the high dimensional data being generated by a finite element DNS method, each of the spatially dependent quantities can be thought of as fluctuating about their volume averaged value, and could be represented by being expanded as a Taylor series
\begin{align}
\functional[](\orderparameter[][]) = \Functional[] + \epsilon + \dotproduct[{\derivative[1]{\functional[](\tilde{\orderparameter[][]})}{\orderparameter[][]}}][{(\orderparameter[][]-\tilde{\orderparameter[][]})}] + O((\orderparameter[][]-\tilde{\orderparameter[][]})^2), \label{eq:local_avg_quantities}
\end{align}
where $\tilde{\orderparameter[][]} = \argmin_{\orderparameter[][]} \abs{\functional[](\orderparameter[][])-\Functional[]}$, such that $\functional[](\tilde{\orderparameter[][]}) = \Functional[] + \epsilon$. This polynomial series means any dependencies of the spatially dependent quantities on the order parameters will become higher order polynomial dependencies in the volumed averaged quantities; $\Functional[]$ will be a higher order polynomial in $\Orderparameter[]$. For example, the free energy density contains up to fourth order polynomials of the order parameters, however the total free energy will be a function of much higher polynomials of the effective order parameter quantities. These effects will be discussed when determining forms of the reduced order models while being informed of the original dependencies that yield the training data.

The relevant $p$-dimensional state vector $\var[][] = \{\var\}$ associated with each vertex $i$ of the graph is chosen to be
\begin{equation}
	\var[][] = \{\energy[][],\energy[][{\textrm{mech}}] ,\Estrain,\volfrac[][\vardim],\interfacial[][\vardim],\numberphase[][\vardim]\}, \label{eq:observables}
\end{equation}
and vertices in this graph are associated with the time index. From physical arguments, the microstructure states at different times are all related, and therefore the graph is considered to be fully connected when considering edges and edge weights.

\begin{figure}[hpt]
\centering
\includegraphics[width=0.6\textwidth]{figures/observables_DNS_updated.pdf}
\caption{DNS calculated state vector components.}
\label{fig:DNSobservables}
\end{figure}

\subsection{Reduced order modelling} \label{sec:physicalsystems_reducedordermodelling}
To determine a model for the system, two different aspects of the model are considered, the time dependence of the phase volume fractions, when zero additional mechanical stress is applied to the system, and the free energy dependence on the strain and phase volume fractions. The forms of these reduced order models will follow the general principals of how the various fields are coupled from this variational framework. This choice of PDE's involves polynomial couplings between the composition and strain fields in the free energy density in \cref{eq:psi_hom}, however the constitutive relations for the chemical potential do not involve derivatives of the free energy density with respect to the strains or strain gradients, and the constitutive relations for the elasticity do not involve derivatives of the free energy density with respect to the compositions or composition gradients. These polynomial and derivative couplings indicate important quantities of interest are those dominated by either algebraic or differential forms of the chemical and mechanical effects.

\subsubsection{Selection of physically relevant basis} \label{sec:physicalsystems_reducedordermodelling_Basis}
Given the microstructure data is generated from explicit forms for the free energy density in \cref{eq:psi_hom,eq:psi_grad}, these exact forms give an indication for which algebraic and differential terms should be proposed for the model. The free energy contains no explicit spatial dependence, and contains monomials up to $\composition[][]^4$ and ${{\estrainspatialindex[][\vardim]}}^4$ order, plus coupled polynomial terms that are up to $\composition[][] {\estrainspatialindex[][2]}^2$ order, plus quadratic gradient terms. The chemical potential and associated current are of the form
\begin{align}
	\chemicalpotential ~=&~ \derivative[1]{\energydensity[][]}{\composition[][]} - \consttwo[][] \laplacian[][{\composition[][]}] \\
	\divergence[][{\current[][]}] ~=&~ -L\left[\uniderivative[3]{\energydensity[][]}{\composition[][]}\dotproduct[{\gradient[]{\composition[][]}}][{\gradient[]{\composition[][]}}] + 2\derivative[3]{\energydensity[][]}{\composition[][],\composition[][],{\estrainspatialindex[][\vardim]}}\dotproduct[{\gradient[]{\composition[][]}}][{\gradient[]{{\estrainspatialindex[][\vardim]}}}] + \derivative[3]{\energydensity[][]}{\composition[][],{\estrainspatialindex[][\vardim]},{\estrainspatialindex[][\vardimone]}}\dotproduct[{\gradient[]{{\estrainspatialindex[][\vardim]}}}][{\gradient[]{{\estrainspatialindex[][\vardimone]}}}] \right.\\
	&~~~~~~~~~~~ + \left.\uniderivative[2]{\energydensity[][]}{\composition[][]}\laplacian[][{\composition[][]}] +\derivative[2]{\energydensity[][]}{\composition[][],{\estrainspatialindex[][\vardim]}}\laplacian[][{\estrainspatialindex[][\vardim]}] - \consttwo[][] \laplacian[][]\laplacian[][{\composition[][]}]\right] \nonumber.
\end{align}
\noindent The equations are solved in the weak form over the material domain, and therefore any spatial dependencies, and particularly any spatial derivatives of the order parameters, are removed, yielding terms that are proportional to the parameter being differentiated, similar to \cref{eq:local_avg_quantities}. For example $\int_{\Omega_t}\integrand[]{V} \nlaplacian[][{\orderparameter[][]}][n] \sim O(\orderparameter[][])$. Therefore the Cahn-Hilliard first order dynamics for the composition should be of the order
\begin{align}
	\derivative[1]{\composition[][]}{t} &\sim O(\uniderivative[3]{\energydensity[][]}{\composition[][]} \composition[][]^2) + O(\derivative[3]{\energydensity[][]}{\composition[][],\composition[][],{\estrainspatialindex[][\vardim]}} \composition[][] {\estrainspatialindex[][\vardim]})
	+ O(\derivative[3]{\energydensity[][]}{\composition[][],{{\estrainspatialindex[][\vardim]}},{{\estrainspatialindex[][\vardimone]}}}~ {{\estrainspatialindex[][\vardim]}}{\estrainspatialindex[][\vardimone]}) \label{eq:cahnhilliard_approx}	\\
	&~ + O(\uniderivative[2]{\energydensity[][]}{\composition[][]} \composition[][]) + O(\derivative[2]{\energydensity[][]}{\composition[][],{\estrainspatialindex[][\vardim]}} {\estrainspatialindex[][\vardim]}) + O(\composition[][]) \nonumber.
\end{align}
For the elasticity equations, given the free energy and its derivatives, upon being integrated over the material volume are polynomials in the order parameters, it will be assumed that a functional representation of the free energy can be expressed as a finite polynomial series in the phase volume fractions and volume averaged strain parameters.

When performing this reduced order modelling of the volume averaged quantities, we will be informed by the various equations and relations used to produce this data in determining which operators should be present in the models. Both the specific terms present, the form of the higher dimensional representations, such as polynomial expansions, differential relationships, or general dependencies, as well as the number of terms required to form the model will inform the hypothesis of the model. Specific implementation details are discussed in \cref{app:implementation_stepwiseregression} of the Appendix, as well as in Wang \emph{et al.}~ \cite{Wang2019,Wang2020} and Kochunas \emph{et al.}~ \cite{Kochunas2020}.

Here, all partial derivatives are approximated using the non-local calculus definitions in \cref{eq:diff_p}, with the exception of the time derivatives, which are calculated using the standard backwards Euler finite difference method. It is also important to note when computing these operators and forming a basis, finite datasets mean numerically some terms may be not linearly independent. This will affect the span of this basis, and if there are more linearly dependent terms in the model, some models with different subsets of operators may be degenerate. This is an effect that will be discussed for the specific models below.

\subsubsection{Weighted residual loss functions} \label{sec:physicalsystems_reducedordermodelling_lossfuncs}
An important aspect of the fitting procedure is the choice of loss function to characterize the error between the known data and predictions. For general models, we choose the root mean square error with the $\normspace[][2]$ norm. The $\normspace[][2]$ norm loss generally minimizes fluctuations in the data and leads to smooth curves, whereas the $\normspace[][1]$ norm loss generally minimizes error at specific points, and the $\normspace[][\infty]$ norm bounds the maximum error of a single point. For highly oscillatory data, such as mechano-chemical systems which can change their local phase rapidly during their transient response, a weighted loss function with several norms is more appropriate. Here, we define the total loss function as a weighted sum of loss functions
\begin{align}
	\loss[][] = \sum_{\chi \in \{1,2,\infty,\cdots\}} \weight[][\chi]\cdot \loss[][\chi], \label{eq:lossfunc}
\end{align}
where $\chi$ denotes the type of norm $\normspace[][\chi]$ used in the loss function, with weighting $\weight[][\chi]$. Each loss function $\loss[][\chi]$ may have its own sub-domain to evaluate the loss. In general, we choose the $\loss[][2]$ domain to be over the entire dataset. For oscillatory problems we choose the $\loss[][1]$ loss to be evaluated at the local maxima and minima of the data, through a specified number of passes of a peak finding algorithm.

\subsubsection{Regression approaches} \label{sec:physicalsystems_reducedordermodelling_regression}
The general modelling approach will consist of proposing a basis of algebraic and differential operators. System inference methods, specifically linear stepwise regression approaches, allow the model to be made as parsimonious as possible, while remaining an accurate representation of the underlying physics.

Here, we choose to perform regression using Ordinary Least Squares (OLS), as well as Ridge regression (Ridge)\cite{Lv2009} and Lasso regression (Lasso)  \cite{Tibshirani1996}, schemes which constrain the norm of the coefficients. Regularization is added to the regression with the objective of improving the robustness of the method, particularly for low rank basis of operators that induce very high condition numbers in the linear problems to be solved. In the case of Ridge regression, it can be directly shown that regularization increases the rank and decreases the collinearity of the data. Please refer to \cref{app:implementation_linearregression_ridge} and \cref{app:implementation_linearregression_lasso} of the Appendix for details on these regression approaches and the residual equations being minimized.

The Ridge and Lasso regression approaches generally use a cross validation scheme to determine the optimal ratio $\lambda$ of the primary constraint and the additional norm constraints in the total residual loss function. At each iteration of the stepwise regression, we perform an independent cross validation to determine the optimal $\lambda$ hyper-parameter for each given reduced basis of operations with a proposed operator removed. This approach increases the computational complexity at least linearly with the size of the data, and linearly with the number of possible hyperparameters, however potentially ensures more optimal coefficients that are away from a local basin of the loss function. We observe that the Lasso regression induces sparsity in the regression coefficients due to using a $1$-norm constraint on the coefficients. This conflicts with the stepwise regression approach, and appears less suitable for these types of problems. We therefore will be comparing the OLS and Ridge approaches.

This effect of the independent cross validations yielding optimal coefficients for a given proposed model at a stepwise iteration, but not accounting for the overall stepwise procedure that is selecting which coefficients to remove from the model, will be seen to affect the overall trends of the model loss. There is potentially a competing interest between the penalization of the loss function with the coefficient norm constraints, and the desire to remove insignificant terms and produce the most effective model.

\textbf{Remark}: It is also important to acknowledge that the cross validation may be performed using a different optimality metric compared to the loss function for the statistical tests in the stepwise regression. The global minima of the two metrics do not necessarily coincide, and as will be discussed in \cref{sec:physicalsystems_vol}, this greatly affects the model, and the order of operators present after the overall stepwise regression is complete. Ultimately, we desire a robust, accurate, and efficient method, and the following example studies indicate the compromises to be made between these objectives.

\subsection{Phase volume fraction dynamics}\label{sec:physicalsystems_vol}
To model the phase volume fraction dynamics, we are guided by the Cahn-Hilliard equation in \cref{eq:pde_cahnhilliard}, and the order of terms present in \cref{eq:cahnhilliard_approx}. Here the composition evolves under first-order dynamics driven by a conservative flux that itself has a non-linear, and high-order gradient dependence on the composition and strains. Importantly, the flux is dependent on the chemical potential of $\composition[][]$ and stresses $\Estrain$ that arise as variational derivatives of the free energy with respect to the composition and strains, respectively in \cref{eq:cflux}. We therefore seek evolution equations for the volume phase volume fractions $\volfrac$ that are also first-order in time, and dependent on polynomials of products of the free energy derivatives up to $\difference[3]{\energy[][]}{\volfrac[][\vardim],\volfrac[][\vardimone],\Estrainindex[][{\vardim \vardimone}]}$, with phase dependent terms in $\{\volfrac[][\vardim],\interfacial[][\vardim], \numberphase[][\vardim],\Estrain\}$. We further restrict the possible dependencies based on the form of the free energy in $d$ dimensions, a function of a single chemical parameter and ${d(d+1)}/{2}$ mechanical parameters, which we choose to be $\Volfrac[][] = \volfrac[][\square]$, and $\Estrain$.

We also assume that any dependencies on the microstructure parameter terms $\interfacial[][\vardim],\numberphase[][\vardim]$ show up only in the terms that contain polynomials of the state variables, the last term in \cref{eq:cahnhilliard_approx}, and not in terms that contain derivatives. These polynomial terms without derivatives arise in the Cahn-Hilliard equation from the interfacial gradient terms in the free energy. The interfacial lengths and particle numbers, here defined to be $\Interfacial[][1] = \interfacial{}{\square},~\Numberphase[][1] = \numberphase{}{\square},~\Interfacial[][2] = {(\interfacial[][\hrectangle]+
\interfacial[][\vrectangle]-\interfacial{}{\square})}/{2},~\Numberphase[][2] = \numberphase[][\hrectangle]+\numberphase[][\vrectangle]$ will suffice to model the $\consttwo[][]$ and $\constthree[][]$ interfacial terms.

\subsubsection{Proposed model for phase volume fraction dynamics}\label{sec:physicalsystems_vol_proposedmodel}
The model for the phase volume fraction first order dynamics is assumed to be of the following form
\begin{equation}
	\derivative[1]{\Volfrac[][]}{t} = f(\nderivative[q]{\energy[][]}{\Volfrac[][]}{\Estrainindex[][{\vardim \vardimone}]},\Volfrac[][],\Estrainindex[][{\vardim \vardimone}],\Interfacial[][\vardim], \Numberphase[][\vardim]), \label{eq:model_vol}
\end{equation}
where $\nderivative[q]{\energy[][]}{\Volfrac[][]}{\Estrainindex[][{\vardim \vardimone}]}  = \{\derivative[2]{\energy[][]}{\Volfrac[][],\Estrainindex[][{\vardim \vardimone}]},\derivative[2]{\energy[][]}{\Volfrac[][],\Volfrac[][]},\derivative[3]{\energy[][]}{\Volfrac[][],\Estrainindex[][{\vardim \vardimone}],\Estrainindex[][\vardimtwo \vardimthree]},\derivative[3]{\energy[][]}{\Volfrac[][],\Volfrac[][],\Estrainindex[][{\vardim \vardimone}]},\uniderivative[3]{\energy[][]}{{\Volfrac[][]}}\}$ are the relevant derivatives for this model.

The function $f$ for this work is an algebraic, linear function of polynomials of the state vector components and derivatives of the free energy
\begin{equation}
	f = \mspace{-20mu}\sum_{{\substack{a = \left\{a_{\nderivative[q]{\energy[][]}{\Volfrac[][]}{\Estrainindex[][{\vardim \vardimone}]}}\cdots~ a_{\Numberphase[][\vardimone]}\right\}\\\abs{a}\leq A\\ 0 \leq a_{\nderivative[q]{\energy[][]}{\Volfrac[][]}{\Estrainindex[][{\vardim \vardimone}]}} \neq a_{\Interfacial[][\vardim]} \neq a_{\Numberphase[][\vardim]} \leq 1}}} \mspace{-40mu} \coef[a][]~ ({\nderivative[q]{\energy[][]}{\Volfrac[][]}{\Estrainindex[][{\vardim \vardimone}]})^{a_{\nderivative[q]{\energy[][]}{\Volfrac[][]}{\Estrainindex[][{\vardim \vardimone}]}}}({\Volfrac[][]})^{a_{\Volfrac[][]}} ({\Estrainindex[][{\vardim \vardimone}]})^{a_{\Estrainindex[][{\vardim \vardimone}]}}({\Estrainindex[][{\vardimtwo \vardimthree}]})^{a_{\Estrainindex[][{\vardimtwo \vardimthree}]}}({\Interfacial[][\vardim]})^{a_{\Interfacial[][\vardim]}} ({\Numberphase[][\vardimone]})^{a_{\Numberphase[][\vardimone]}}} \label{eq:model_volfunc}
\end{equation}
where the powers of the quantities are represented by components of the multicomponent index $a$, with the indices being denoted by the quantities themselves. Here, summation over repeated indices in a given term should not be assumed. Instead, the function is explicitly a sum over all combinations of products of powers of the quantities, including all combinations of components of the tensor quantities. The powers are constrained such that the derivative $a_{\nderivativeinline[q]{\energy[][]}{\Volfrac[][]}{\Estrainindex[][{\vardim \vardimone}]}}$, interfacial length $a_{\Interfacial[][\vardim]}$, and number $a_{\Numberphase[][\vardim]}$ terms are all separate from each other and never appear more than once in a term, as per \cref{eq:cahnhilliard_approx}. The polynomial terms are up to order $A=3$ powers, and there are a maximum of $308$ terms in the model.

\subsubsection{Numerical implementation of phase volume fraction dynamics}\label{sec:physicalsystems_vol_numerical}
Given this basis for the terms that drive the phase volume fraction dynamics, stepwise regression is performed using a mixture of the $\loss[][2]$ and $\loss[][1]$ loss functions, as per \cref{eq:lossfunc}. Various loss weights are chosen and a single pass of a peak finding algorithm is used to find the points at which the $\loss[][1]$ loss is computed. The resulting loss curves with a comparison of OLS and Ridge regression fitting are shown in \cref{fig:vol_loss_OLS_Ridge17}. From initial leave-one-out cross-validation studies, it is found that any ridge parameter $\lambda > 0$ improves the rank of the matrix, and for this problem, the fit losses increase with increased $\lambda$. Therefore, $\lambda = 10^{-17}$ is fixed to perform Ridge regression. The optimal direct solver for Ridge regression for this problem is found to be an SVD based solver, whereas OLS regression uses a standard least squares solver. Please refer to a description of the graph theory library in \cref{app:implementation_graphtheorylibrary} of the Appendix for details on the specific direct linear solvers.

Initial studies are also performed using an $R^2$ correlation coefficient for the cross validation metric, while using the weighted loss function for the stepwise regression statistical criteria. This choice of different metrics results in significantly higher errors, of the order $O(10^{-1})$ and the stepwise regression procedure yields a significantly more oscillatory loss curve, as opposed to the desired monotonically increasing curve. The optimal coefficients given the correlation coefficient or the weighted loss metric do not coincide and the cross validation and stepwise regression are competing to the detriment of the overall method. Therefore all cross validation studies will be performed using an identical loss function for all metrics.

Each choice of loss weighting yields similar trends in the loss curves increasing with increased model parsimony, offset by how much the $\loss[][1]$ contributes to the total loss. The high frequency oscillations, as per \cref{fig:DNSobservables} in the phase volume fractions suggests why the fits with solely $\loss[][1]$ have the lowest loss. When looking at the exact ordering of operators that the stepwise regression presents as most important to the model, the different loss weights do not result in identical orderings, even at low numbers of operators where the loss curves are more similar. However, the general trends of which operators are present at various stages of the stepwise regression, denoted by the various plateaus in the loss curves, are very similar across choices of loss weights. We therefore choose to focus on analyzing the fits where solely the $\loss[][1]$ loss is used.

\subsubsection{Reduced order model for phase volume fraction dynamics}\label{sec:physicalsystems_vol_finalmodel}
Looking at the plateaus and sharp increases in the weighted loss curves, particularly for the OLS model in \cref{fig:vol_loss_OLS} over the stepwise iterations, the first plateau extends from $308$ down to approximately past $187$ operators, which corresponds to the rank of the matrix of this basis of operators. In general, as the model becomes more parsimonious, the loss should increase monotonically, however in these cases of lower rank basis, there may be some decreases when certain operators are rejected from the model. The second plateau ends at approximately $47$ operators, which is also the point in the stepwise iterations where the model first has a term that contains a derivative. Practically all iterations between $47$ and $308$ have derivative terms, and this suggests that the purely algebraic polynomial terms in $\{\volfrac[][\vardim],\interfacial[][\vardim], \numberphase[][\vardim],\Estrain\}$ are over-fitting the model. The derivatives are not seen to play as important a role in the first order dynamics. However, the suggested importance of the algebraic polynomial terms may also be explained by higher order polynomials of the volume averaged quantities being required to represent dynamics controlled by spatially varying fields, as per \cref{eq:local_avg_quantities}.

With loss weights $\weight[][2] = 0,\weight[][1] = 1$, the first 10 terms from the stepwise regression for the first order dynamics of $\Volfrac[][]$ are found using OLS and Ridge regression to be
\input{./figures/model_vol_OLS_Ridge17.tex}
It can be seen that both models consist of solely algebraic terms. The OLS and Ridge approaches also generate models with similar terms, with the Ridge regression model having 9 terms compared to 8 terms for OLS regression that contain predominantly strain dependencies.

The first order dynamics, and fits for models with $20$, $50$, and $200$ terms is shown in \cref{fig:vol_fit_OLS_Ridge17}. For the practically full model with $200$ terms, the fits correspond to the data very well at early times, but have difficulty matching all of the small oscillations at late times. As the model becomes more parsimonious to $50$ terms where the sharp increase in the loss occurs, the general trend of the dynamics is captured, however this model misses capturing the full amplitude of the major peaks compared to the more complex models. Finally, with $20$ terms, which corresponds to the number of terms present in the original Cahn-Hilliard equation in $d=2$ dimensions, the large peaks, and small oscillations are not entirely captured, however the general trends are still shown. As discussed, the algebraic polynomial terms are potentially over-fitting, however manage to retain the general behavior of the phase volume fractions changing over time.

\subsubsection{Regression approaches for phase volume fraction dynamics}\label{sec:physicalsystems_vol_regression}
When comparing the performance of the OLS and Ridge regression approaches, both approaches are shown to yield similar trends in the loss curves in \cref{fig:vol_loss_OLS_Ridge17}. However the Ridge regression loss is shown to be greater, by less than an order of magnitude for more complex models than the OLS regression. The resulting fits for the Ridge regression in \cref{fig:vol_fit_OLS_Ridge17} are also shown to capture the general trends of fluctuations quite well, however do not fit the sharper oscillations at later time quite as well as the OLS approach. This increased Ridge loss is reasonable considering the additional constraints, and is balanced by improving the smoothness of the loss curves over the stepwise iterations. The smoother Ridge curves subsequently have less plateaus, and the losses for different loss weights are more distinct. This approach appears to be less affected by changes in the rank of the matrix. Both approaches also behave similarly at the most parsimonious models with less than approximately $47$ terms. These comparable results indicate that Ridge regression is the most appropriate approach for future studies, increasing the robustness of the modelling to more collinear data, and still capturing the essential physical trends. Future cross validation studies, as well as the use of additional constraints on the coefficients in conjunction with the stepwise regression will be useful to determine how to best select the most parsimonious model.

The importance of choosing a physics informed basis is shown clearly by this example. Other choices of basis, such as including other derivatives, or a different set of basis variables like the other phase volume fractions, were shown to have much worse regression losses, and were much less intuitive. Future applications of this graph theoretic method should ensure the forms of the original models in the higher dimensional data are considered when determining how to capture effective physical processes.

\begin{figure}[hpt]
\centering
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=\textwidth]{figures/plot_Loss_vol_OLS.pdf}

	\subcaption{OLS regression.}
	\label{fig:vol_loss_OLS}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=\textwidth]{figures/plot_Loss_vol_Ridge17.pdf}

	\subcaption{Ridge regression with $\lambda = 10^{-17}$.}
	\label{fig:vol_loss_Ridge17}
\end{subfigure}
\caption{Stepwise regression loss curves for the phase volume fractions first order dynamics. The iterations use a weighted proportion of $\loss[][2]$ and $\loss[][1]$ loss functions for the total loss, and shown are the curves with the $\loss[][2]$ weight.}
\label{fig:vol_loss_OLS_Ridge17}
\end{figure}

\begin{figure}[hpt]
\centering
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=\textwidth]{figures/plot_BestFit_vol_200_50_20_OLS.pdf}

	\subcaption{OLS regression.}
	\label{fig:vol_fit_OLS}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=\textwidth]{figures/plot_BestFit_vol_200_50_20_Ridge17.pdf}

	\subcaption{Ridge regression with $\lambda = 10^{-17}$.}
	\label{fig:vol_fit_Ridge17}
\end{subfigure}
\caption{First order dynamics of phase volume fractions for $20$, $50$, and $200$ terms in the forcing function, and use of the $\loss[][1]$ loss in stepwise regression. Backward Euler first time derivative data is shown with the dashed blue curve.}
\label{fig:vol_fit_OLS_Ridge17}
\end{figure}

\subsection{Free energy functional representation} \label{sec:physicalsystems_psi}
To model the free energy of the system, we are guided by the polynomial form of the free energy density in \cref{eq:psi_hom} and write
\begin{align}
	\energy[][] = \energy[][](\Volfrac[][],{\Estrain}) \label{eq:model_psi}
\end{align}
as a function of the single chemical and ${d(d+1)}/{2}$ mechanical variables in $d$ dimensions.

\subsubsection{Proposed model for free energy functional representation}\label{sec:physicalsystems_psi_proposedmodel}

The form of this model suggests our functional representation should be a Taylor series at time $j$, about a base value at the time $i$:
\begin{align}
	\energy[][](\Volfrac[][j],\Estrain[][j]) = \energy[][](\Volfrac[][i],\Estrain[][i])
	~&+~ \coef[{\Volfrac[][]}][] \derivative[1]{\energy[][](\Volfrac[][i],\Estrain[][i])}{\Volfrac[][]} \Delta {{}{\Volfrac[][]}}_{ij} ~+~ \coef[{\Estrainindex[][{\vardim \vardimone}]}][] \derivative[1]{\energy[][](\Volfrac[][i],\Estrain[][i])}{\Estrainindex[][{\vardim \vardimone}]} \Delta {{}\Estrainindex[][{\vardim \vardimone}]}_{ij} \label{eq:model_psitaylor} \\
	~&+~ \coef[{\Volfrac[][]}{\Volfrac[][]}][] \uniderivative[2]{\energy[][](\Volfrac[][i],\Estrain[][i])}{{\Volfrac[][]}} {\Delta {{}{\Volfrac[][]}}_{ij}}^2 \nonumber\\
	 ~&+~ \coef[{\Estrainindex[][{\vardim \vardimone}]}{\Estrainindex[][{\vardimtwo \vardimthree}]}][] \derivative[2]{\energy[][](\Volfrac[][i],\Estrain[][i])}{{\Estrainindex[][{\vardim \vardimone}]},{\Estrainindex[][{\vardimtwo \vardimthree}]}} \Delta {{}{\Estrainindex[][{\vardim \vardimone}]}}_{ij} \Delta {{}{\Estrainindex[][{\vardimtwo \vardimthree}]}}_{ij} \nonumber\\
	 ~&+~ \coef[{\Volfrac[][]}{\Estrainindex[][{\vardim \vardimone}]}][] \derivative[2]{\energy[][](\Volfrac[][i],\Estrain[][i])}{{\Volfrac[][]},{\Estrainindex[][{\vardim \vardimone}]}} \Delta {{}{\Volfrac[][]}}_{ij} \Delta {{}{\Estrainindex[][{\vardim \vardimone}]}}_{ij} \nonumber \\
	 ~&+~ O(\Delta {{}{\Volfrac[][]}}_{ij}^3) + O(\Delta {{}{\Estrainindex[][{\vardim \vardimone}]}}_{ij} ^3), \nonumber
\end{align}
where ${\Dvar[][ij]} = \var[][j] - \var[][i]$ represents the change in the $\volfrac[][\vardim]$ or $\Estrain$ state variable between times $j$ and $i$. Please refer \cref{app:implementation_modifiedtaylorseries} of the Appendix for details on how the Taylor series are formed. Similar to the model for the phase volume fraction evolution, we choose to restrict the dependencies to be $q=4$ variables of a single chemical and three mechanical order parameters in $d=2$ dimensions. We further restrict the Taylor series to be up to $k = \ith[4]$ order, and there are $p = {(q^{k+1}-1)}/{(q-1)} = 341$ possible terms in the model. Here, the expansion is about the initial time $i=0$, and all partial derivatives in the Taylor series are approximated using the non-local calculus definitions in \cref{eq:diff_p}.

\subsubsection{Numerical implementation for free energy functional representation}\label{sec:physicalsystems_psi_numerical}
From initial leave-one-out cross-validation studies, it is found that any ridge parameter $\lambda > 0$ improves the rank of the matrix, and for this problem, decreases the fit losses compared to the OLS method. The fit losses also increase with increased $\lambda$. Therefore, $\lambda = 10^{-17}$ is fixed to perform Ridge regression. Similarly to for the phase volume fractions, the optimal direct solver for Ridge regression for this problem is found to be an SVD based solver, and the OLS regression is performed with a standard least squares solver.

The Ridge regression loss curve in \cref{fig:psi_loss_OLS_Ridge17} is slightly smoother than the OLS regression, and unlike in the phase volume fractions problem, does not increase the loss with $\lambda > 0$. The regularization is shown to be able to increase the rank of the operators, and the robustness of the method, while still producing an accurate model. Here, the free energy is a monotonically decreasing function and is very smooth, so the loss function is the standard $\loss[][2]$ error.

\subsubsection{Reduced order model for free energy functional representation}\label{sec:physicalsystems_psi_finalmodel}
From the loss curves in \cref{fig:psi_loss_OLS_Ridge17}, the sharp increase in the loss occurs at around $70$ terms, which corresponds to the rank of this basis of Taylor series terms. The loss increasing sharply for models with less than this number of terms indicates the importances of these terms in this truncated Taylor series. The loss for the full model is of order $O(10^{-5})$, comparable to other methods \cite{Zhang2020}.

The first 10 terms from the stepwise regression for the $\ith[4]$ order Taylor series functional representation for $\energy[][]$ are found for OLS and Ridge regression to be:
\input{./figures/model_psi_OLS_Ridge17.tex}

\noindent Here, we have imposed for model consistency, that the coefficient associated with the base Taylor series term is $\gamma_0 = 1$. The most important higher order terms for both models are primarily second and third order derivatives, and both the OLS and Ridge models depend most strongly on the second derivative with respect to phase volume fraction. The OLS model has only one fourth order term, whereas the Ridge model has two fourth order terms, and a first order term. This possibly corresponds to the original free energy density $\energydensity[][]$ in \cref{eq:psi_hom} being second,third and fourth order powers of the order parameters. However care must be taken on the interpretation of the results of the stepwise regression as being directly induced by the forms of the models for the higher dimensional data. In this example, the stepwise regression indicates a direct correspondence between the form of the reduced order model for $\energy[][]$, and the form of the spatially dependent $\energydensity[][]$, however this correspondence should not always be assumed. Physics informed initial choices of basis will ensure the most appropriate reduced order model is found, and later interpretation is specific to each problem.

During the modelling process, it is observed that this Taylor series basis without any regularization is very singular, with rank $70 \ll p$, and condition number on the order of $O(10^{18})$. This can be attributed to the sparsity of the data, as well as potentially the form of the free energy, being a monotonically decreasing function, and derivatives with respect to various quantities being possible numerically linearly dependent. This singular nature of this basis means standard linear regression, computed through an SVD algorithm to solve the normal equations, as per \cref{app:implementation_linearregression} of the Appendix, sets certain singular values to zero, as per a given tolerance. The low rank of this matrix means many singular values are close to zero, and floating point issues related to the machine precision may occur.

Regarding the fits for the free energy, Taylor series with $10$, $30$, and $70$ terms are shown in \cref{fig:psi_fit_OLS_Ridge17}. Models beyond the rank of this basis fit the free energy very well and model a very smooth function. More parsimonious models capture the overall decrease of the free energy, however oscillations in the phase volume fractions and strains are not able to be as smoothed out. Free energies at earlier times also appear to be fit better by all models. When comparing the OLS and Ridge fits, the fits are almost identical for more complex models, however the OLS method has smaller oscillations for the most parsimonious models. This Taylor series approach is shown to be a logical and effective basis of terms to model this smooth free energy function.

\begin{figure}[hpt]
\centering
\includegraphics[width=0.5\textwidth]{figures/plot_Loss_psi_OLS_Ridge17.pdf}
\caption{Stepwise regression loss curve for the $\ith[4]$ order free energy Taylor series functional representation, using OLS and Ridge regression.}
\label{fig:psi_loss_OLS_Ridge17}
\end{figure}

\begin{figure}[hpt]
\centering
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=\textwidth]{figures/plot_BestFit_psi_70_30_10_OLS.pdf}

	\subcaption{OLS regression.}
	\label{fig:psi_fit_OLS}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=\textwidth]{figures/plot_BestFit_psi_70_30_10_Ridge17.pdf}

	\subcaption{Ridge regression with $\lambda = 10^{-17}$.}
	\label{fig:psi_fit_Ridge17}
\end{subfigure}
\caption{Fitted curves for the $\ith[4]$ order free energy Taylor series functional representation, with $10$, $30$, and $70$ terms. DNS free energy data is shown with the dashed blue curve, and is essentially coincident with the higher order models.}
\label{fig:psi_fit_OLS_Ridge17}
\end{figure}

\section{Conclusions}
\label{sec:concl}
In this work, we develop a graph theoretic approach for reduced order modelling of physical systems. After extending a non-local calculus as the foundation of this approach, we observe there are challenges and compromises when imposing consistency with differential calculus definitions of partial derivatives. Constraints placed on the graph's edge weights lead to numerical issues, but also interpretation issues if the non-local derivatives have error constants that are consistently large. A more in depth comparison of the behavior of non-local and differential derivatives is necessary in future works.

Applying this calculus to studying volume averaged quantities from high dimensional direct numerical simulations by representing these quantities on a graph shows the intuitive, general, and effective nature of this approach. With an example physical system of the phase evolution of a mechano-chemical multi-crystalline solid, we see that using a physics informed basis of operators can model both differential equations, and functional representations. Highly oscillatory or smooth monotonically evolving functions can be fit using linear regression and an adequate number of terms, with errors on the order of $O(10^{-4})$. As a comparison, recent studies of the same microstructure system using neural networks \cite{Zhang2020} achieve errors on the order of $O(10^{-6})$.

We also note that near the beginning of the stepwise regression when the models contain the full set of operators, there is sometimes, for example in \cref{fig:vol_loss_OLS_Ridge17}, decrease in the loss function as the model becomes initially more parsimonious. This suggests that although we are choosing an intuitive, physics informed basis, there may be some operators that worsen the model, and must be rejected empirically. A more rigorous study to more precisely determine which operators from the appropriate basis should be included in the initial model will be conducted in future communications.

\textbf{Remark}: There remains an open question of how accurate these graph theoretic methods can be. This question may be related to the non-local nature of the derivatives, and whether they are appropriate to model strictly locally behaving systems. Whether or not these methods can achieve the error tolerances of other direct numerical simulations, or methods such as neural networks, will show the true effectiveness of this approach. There is also the compromise of the robustness and consistency of this modelling approach versus its accuracy. As discussed in \cref{app:implementation_linearregression} of the Appendix, the improvement of the multicollinearity of the basis of operators from regularization approaches may come at the cost of decreased loss. The physical systems studied indicate a physically appropriate basis of operators may still be numerically collinear, and regularization aids in the robustness and consistency of the modelling, without significant loss in accuracy.

We note, however, that in any reduced-order model, the goal is a parsimony of representation, which is achieved here by $O(10)$ components in the state vector, at the cost of some loss of fidelity. This approach appears to require $O(10^2)$ terms in order to get accurate fits of the data, however these higher order models involve higher order derivatives that may suffer numerical errors due to their large values. This compromise of model accuracy versus model complexity is another open question to be studied. To justify the requirements of these higher order models, we note that since the original PDE itself in the Cahn-Hilliard equation has nearly 20 terms when fully expanded, this number could be regarded as a lower bound for the complexity of the reduced-order models. Then, with the significant saving in dimensionality of the state vector over that of a DNS, even $O(10^2)$ operators in the reduced-order model is not unreasonable. This is especially reasonable when this representation is required to reproduce the dynamics of the vastly higher dimensional DNS.

Given these findings, this graph theoretic approach is shown to be very general, and many applications to a wide range of physical system, and system scales are possible. There are also several interesting studies that can be conducted in future works regarding the methods themselves. First it will be important to develop edge weight constraints for finite, discrete graphs, as opposed to current use of constraints for finite continuous graphs. Rigorous approximations for these discrete sums over the volume of the graph will possibly constrain the numerical issues brought on by the consistency demands, and will give allow the local and non-local calculus to be better compared. The forms of the decaying edge weights can also be generalized to include several decay parameters that reflect independent behaviors of the state vector components. Finally, aspects of the stepwise regression, including the effects of different loss functions, and how the independence of operators in the basis plays a role in how and which operators are rejected from the model, will aid in the fitting process.

\pagestyle{bibliography}

\section{Introduction} \label{sec:intro}
Natural phenomena that are susceptible to a mathematical physics treatment typically lead to initial and boundary value problems. The solution of partial differential equations that govern these problems has occupied the better part of the last century of research in applied mathematics, leading to numerical methods, and from there to computational science and high performance computing. With Fourier methods, time integration algorithms, discretization frameworks of finite elements, finite difference, Green's functions-based formulations, and hybrids of these techniques among other approaches, it is now possible to attain a degree of fidelity that is sufficient for many applications in mathematical physics. While the need to incorporate additional physics and gain a finer resolution of phenomena will never be met, there is also the recognition that the fidelity that is already achievable comes paired with a complexity that outstrips easy representation and analysis. Reduced-order models must therefore be derived on a principled basis with controlled approximations.

For the purpose of distinguishing the approach adopted here, we take account of areas including projection methods, coarse graining techniques, dimensionality reduction, sub-manifold extraction and numerous others that are arguably encompassed by the field of dynamical systems. Within this field, is the graph theoretic approach for representation, exploration and analysis of computed states of physical systems that we have proposed recently \cite{Banerjee2019}. This formalism is predicated on computed solutions that are typically high-dimensional in the sense that the primary unknown fields are represented by degrees of freedom numbering $\gtrsim O(10^8)$. Functionals on these discrete solutions are extracted as \emph{states}. Some examples are the lift, drag or pressure difference in computational fluid dynamics, measures of load, deformation and total energies in solid mechanics, phase volume fractions and free energies in materials physics, and cross sections in nuclear physics \cite{Kochunas2020}. These states undergo \emph{transitions} as some \emph{parameter} of the system is varied, such as its physical, possibly time-like, parameters, boundary or initial conditions.

\subsection{Quantities of interest and graph representations} \label{sec:intro_qoi}
The states, of which there is typically a small number $\sim O(10)$ relevant to a system, form a low-dimensional vector of quantities of interest. However, these state vectors can be of high fidelity limited only by the dimensionality of the full-field numerical solution. For completeness of the description, the state vector can be extended to include the parameters referred to above. We have shown that by defining the states to be vertices on a graph, and the transitions to be edges between the vertices, a near one-to-one correspondence is uncovered between the properties of the computed physical system and the elements of graph theory. Weighted edges between vertices can then be assigned, or found through graph theoretic principles, indicating relationships or transitions between states, and the magnitude of such correlations. Reversible linear systems lead to undirected graphs that are fully connected, whereas dissipative dynamical systems are represented by directed trees. Path-dependence can be thought of as being associated with an absence of cycles. Notions of centrality reveal insights to the relations between states and path traversal properties are induced on graphs by equilibrium or dynamical transitions on the physical system.

Given the above graph theoretic framework we seek to extend it to reduced order modelling based on computed states and transitions of physical systems. We take the view of the graph of states (vertices) and transitions (edges) as a discrete manifold describing the physical system. On this manifold we seek to define differential equations and functional representations for the evolution of certain components of the state vector. An illustrative example that we pursue in this communication is of an ordinary differential equation for the phase fractions of a material system, driven by itself and other components of the state vector such as the free energy and strains. Guided by our knowledge of first-order systems in materials physics, it is of interest to compute derivatives on these states, such as of the free energy, with respect to the phase fractions, yielding generalized chemical potentials, or with respect to the strains, yielding stresses. The motivation comes from the rates of change of phase fractions driven by chemical potentials and stresses in materials. An alternative to the differential equation is a direct functional representation, such as for the free energy with phase fractions and strains as arguments. In this case, we are guided by the universality of the Taylor series representation, to which we seek approximations. This too makes necessary the computation of derivatives of state vector components with respect to each other.

\subsection{Non-local calculus for reduced order modelling on graphs} \label{sec:intro_nonlocalcalculus}
We have adopted a discrete non-local calculus on weighted graphs to develop the reduced order models outlined above. This formalism, which has been applied previously to image analysis \cite{Gilboa2008,Desquesnes2013,Hein2007} begins with definition of a vector space on the graph, followed by the gradient operator, inner products and partial derivatives. Also important is the consistency and rigor reflected in a divergence operator, adjointness relations and a divergence theorem. This non-local theory bears a relation to classical differential calculus, that is determined by the graph weights. Through its self consistently defined inner products, the weights appear as kernels determining the support of integral operators on the graph manifold, and crucially controlling the behavior of partial derivatives relative to their classical definitions.

This forms a central consideration on theoretical and practical grounds. However, the machinery that translates the non-local calculus on graphs to reduced order models can be formally put in place independently of the correspondence between non-local and local treatments. With functions and non-local partial derivatives in hand, algebraic and non-local differential operators can be introduced and viewed as basis terms in the reduced order representation. Since these can be evaluated at the vertex data points, a stepwise regression approach can be invoked to choose the best combination of operators representing the evolution or dependence sought. Model parsimony can be dialed in by using thresholds on the stepwise regression.

\subsection{Physical systems of interest} \label{sec:intro_physicalsystems}
The system of interest for this work are multi-crystalline solids, where the dynamical quantities of interest are the volume fraction of each phase, which are assumed to undergo first-order dynamics. These dynamics are selected to be forced by a linear functional of a $\ird[3]$ order polynomial of the phase volume fractions themselves, the interfacial lengths between the phases, the number of particles, and the first derivatives of the free energy, with respect to the phase volumes, and the mechanical strains. In order to complete the model of the system to be able to forward solve the first order dynamics, a functional representation for the free energy is found as a Taylor series in the phase volume fraction and mechanical strains, up to $\ith[4]$ order. This example allows the general graph theoretic method to be shown, and analyzed for its effectiveness at developing functional representations, and investigating any numerical issues associated with these derivatives. Further examples of physical systems will be explored in future communications.

Computations on each of these examples relies upon a general graph theory library code that accepts general data with inputs and outputs, and determines appropriate models based on these stepwise regression procedures, as well as visualizes fits with plots of loss curves, fitted curves for various models, and trends in the relative magnitudes of fitted coefficients over the stepwise regression. These examples and library allow conjectures on future applications and examples of this graph theoretic approach to be explored.

The non-local calculus on finite weighted graphs is laid out in Section \ref{sec:graphtheory}, illustrative physical systems and discussions are presented in Section \ref{sec:physicalsystems}, and conclusions appear in Section \ref{sec:concl}. Appendices outline the implementation details and example calculations for these methods, including details of stepwise linear regression in \cref{app:implementation_stepwiseregression,app:implementation_linearregression}, implementation of the graph theory code library \cref{app:implementation_graphtheorylibrary,app:implementation_gauss}, and error analyses of the method in \cref{app:error,app:error_p}.

\section{Graph theoretic approach} \label{sec:graphtheory}
The major formalism used in this work relies on methods from graph theory \cite{Banerjee2019,West2001,Newman2010}. A graph, denoted by $G = (V,E)$, consists of a set of vertices, $i \in V$, with $\abs{V} = n$. The vertices are connected by a set of edges, $e \in E$, where each edge $e = (i,j)$ is a pair of vertices. A vertex is denoted as being the part of an edge by $i \in e$. The graph can also have global and local attributes on the vertices and edges, as we describe below. Of importance are edge weights $w$, that are generally functions of the local vertex attributes, and admit a notion of distance between states on the graph as will be elaborated on in \cref{sec:graphtheory_weightconstraints}.

As outlined in \cref{sec:intro_qoi} of the Introduction, in our treatment, a state is a low-dimensional vector of quantities of dimension $p$, obtained after solving a high-fidelity computation of dimension $P$, where ideally $p \ll P$. These high fidelity computations are typically direct numerical simulations (DNS), for example the solution of partial differential equations (PDE) using a finite element-like approach. The DNS simulations can yield some functionals, such as volume averaged fields, or resulting scalar outputs. These, collected into a low-dimensional vector will be denoted by $u$. Separately, \emph{parameters} such as physical parameters, or initial and boundary conditions of the DNS can be collected into another low-dimensional vector denoted by $x \in \reals[p]$. Each state is represented as a vertex $i \in V$ on the graph. It follows that the notation $x_{i}$ is unambiguous in labelling a state by its vertex. Transitions between these states, including changes in system parameters, or steps of the numerical solver, allow for connections to be made between vertices on the graph. These transitions therefore can be thought of as inducing edges $e = (i,j)$ between states $i$ and $j$.

Stated formally now, the vertices, $i \in V$ represent states $\{x^{\mu}_{i},u(x_{i})\}$, where, as introduced above, $x_{i}$ are parameters (physical parameters, initial and boundary conditions) and $u(x_{i})$ are quantities of interest, obtained as functionals from the high-dimensional numerical solution of a system of PDEs. The graphs naturally inherit the mathematical structure of the PDEs \cite{West2001,Newman2010}, and even confer new machinery to it \cite{Banerjee2019}.

Several reoccurring graph structures arise in formal graph theory, including fully connected clique graphs, such as the subgraphs shown in \cref{fig:cliquegraph}, or minimally connected tree graphs, such as shown in \cref{fig:treegraph}. There are intricate correspondences between the physical processes being studied and the resulting graph structure. As identified by Banerjee \emph{et al.}~\cite{Banerjee2019}: directed graphs and less densely connected tree graphs correspond to processes with causality or irreversibility induced by energy dissipation, whereas undirected graphs correspond to reversible processes. Fully connected clique graphs generally represent linear processes in equilibrium. These correspondences, centrality measures and graph traversal will be revisited in future communications.

\begin{figure}[ht]
	\centering

	\includegraphics[width=0.6\textwidth,trim={0 0 0 6cm},clip]{figures/clique_graph.pdf}

	\caption{Example graph, with a central clique, and branching clique subgraphs. Vertex coloring may represent various local attributes.}
	\label{fig:cliquegraph}
\end{figure}

\begin{figure}[ht]
	\centering

	\includegraphics[width=0.6\textwidth,trim={0 0 0 6cm},clip]{figures/tree_graph.pdf}

	\caption{Example directed tree graphs.}
	\label{fig:treegraph}
\end{figure}

\subsection{Non-local calculus on finite, weighted graphs} \label{sec:graphtheory_nonlocalcalculus}
As outlined in \cref{sec:intro_nonlocalcalculus} of the Introduction, we seek to develop reduced order models in the form of differential equations for the evolution of some component of the state vector on the graph, or as functional representations. In the latter case, we explore forms that are motivated by the Taylor Series of local, differential calculus. In both cases, we need to first define a calculus on the discrete manifold that is the graph. Specifically, we adopt a discrete, non-local calculus on finite, weighted graphs \cite{Hein2007,Gilboa2008,Elmoataz2008,Desquesnes2013,Lozes2014}.

Gilboa \emph{et al.}~ define a discrete calculus, consisting of non-local operators, based on differences between states $i,j \in V$, and an edge weight $w(x_{i},x_{j})$\cite{Gilboa2008}. The edge weights $w$ can be defined in a number of ways, and $w(x_{i},x_{j}) = w(x_{j} - x_{i})$ is natural for our purposes. In the case of symmetric, undirected graphs, a radially symmetric weight, $w(x_{i},x_{j}) = w(r(x_{i},x_{j}))$, where $r(x_{i},x_{j}) = \abs{x_{j} - x_{i}}$ is also appropriate. As will be discussed, there are several constraints placed on the allowed forms of the weight functions, to ensure that the non-local calculus converges to true differential operators in certain limits.

Scalars, $u(x_{i})$ correspond to a function at a single vertex $i$. Vectors $v(x_{i},x_{j})$ are functions on vertex pairs $i,j$ or the edge $e = (i,j)$. Accordingly, we have a vector space $\mathnormal{L}$ such that there is a mapping $V\times V \mapsto \mathnormal{L}$ or $E\mapsto \mathnormal{L}$. Integration on the discrete manifold $G$ is a sum over vertices $i \in V$, and is relevant to operations such as inner products, which we define to be normalized by the size of the graph, $\abs{V} = n$. There should be no cause for confusion between the state vectors $x_{i} \in \reals[p]$ represented as vertices $i \in V$, and the edge vectors $v(x_{i},x_{j}) \in \mathnormal{L}$. These are distinct manifestations of vector spaces that are not copies of each other. We also note that for this treatment, there are no self edges and so vectors $v(x_{i},x_{j})$ at vertex $i$ are defined as being between vertices distinct from $i$, and $j \neq i$.

\subsubsection{Non-local calculus definitions} \label{sec:graphtheory_nonlocalcalculus_definitions}
The non-local gradient operator $\gradient[][{u}] = \gradient[{w}][[{u}]](x_{i},x_{j}) $, with respect to the weight $w$, can now be defined as the vector functional of scalars $u$ at state $x_{i}$, and represents the vector of weighted differences with \textit{all other vertices} $j \in V$:
\begin{equation}
	 \gradient[{w}][[{u}]](x_{i},x_{j}) \equiv [u(x_{j})-u(x_{i})]\sqrt{w(x_{i},x_{j})},\quad \forall ~ j \neq i \in V.
\end{equation}
The inner product between scalars $u_{\alpha}$ and $u_{\beta}$ is defined as the sum over the vertices:
\begin{align}
	\langle u_{\alpha},u_{\beta}\rangle \equiv&~ \frac{1}{n}\sum_{i\in V} u_{\alpha}(x_{i})u_{\beta}(x_{i}),
\end{align}
and the contraction between vectors $v_{\alpha}$ and $v_{\beta}$ at vertex $i$ is defined as the sum over disparate vertices:
\begin{align}
	[\dotproduct[{v_{\alpha}}][{v_{\beta}}]](x_{i}) \equiv&~ \frac{1}{n-1}\sum_{j \in V\setminus \{i\}} v_{\alpha}(x_{i},x_{j})v_{\beta}(x_{i},x_{j}).
\end{align}
Inner products, and norms between pairs of vectors are likewise defined.
From these definitions, in the case of symmetric weights, there are parallels to standard differential calculus. There is an adjoint relation between the gradient and the divergence operators, as well as a divergence theorem, leading to a definition of the Laplacian being found to be self-adjoint.

\noindent For this work, the most important operators are the partial derivatives of $u(x_{i})$ with respect to $x^{\mu}$, the $\ith[\mu]$ component of the state $x$. Here, we denote differential calculus partial derivatives with $\partial$, and non-local partial derivatives with $\delta$. The argument of the function in $u(x_{i})$ also makes clear that the function's partial derivative is evaluated at vertex $i$. These non-local operations can be used to rigorously define the partial derivatives as the contraction between the function gradient and the specifically defined unit vector,
\begin{align}
	\difference[1]{u(x_{i})}{x^{\mu}} ~=&~ [\dotproduct[{\gradient[{w}][{[{u}](x_{i})}]}][{\hat{x}^{\mu}_{i}}]] \label{eq:diff{1}} \\
	 ~=&~~ \frac{1}{n-1}\sum_{j \in V\setminus \{i\}} [u(x_{j}) - u(x_{i})](x^{\mu}_{j} - x^{\mu}_{i})w(x_{i},x_{j}),
\end{align}
where the unit vector at state $i$ in the $\mu$ direction is defined as a vector gradient between states $i$ and $j$:
\begin{equation}
	\hat{x}^{\mu}(x_{i},x_{j}) \equiv \gradient[{w}][[x^{\mu}]](x_{i},x_{j}) = [x^{\mu}_{j} - x^{\mu}_{i}]\sqrt{w(x_{i},x_{j})}.
\end{equation}
We will proceed with the weights $w(x_{i},x_{j})$ being chosen such that the unit vectors are normalized:
\begin{equation}
 	[\dotproduct[{\hat{x}^{\mu}}][{\hat{x}^{\nu}}]] = \delta^{\mu\nu}, \label{eq:unitnorm}
\end{equation}
where $\delta^{\mu\nu}$ is the Kronecker delta. As will be discussed, this constraint together with the inner products also being normalized by the volume of the space, leads to other important behaviors of the convergence and consistency of these non-local derivatives to their counterparts in differential calculus.

Higher order derivatives also can be computed as approximations to $\nderivative[q]{u(x)}{x^{\mu_{0}}}{x^{\mu_{q-1}}}$ using an extension of \cref{eq:diff{1}} and can be defined recursively:
\begin{align}
\ndifference[q]{u(x_{i})}{x^{\mu_{0}}}{x^{\mu_{q-1}}} ~=&~
[\dotproduct[{\gradient[{w}][{[
	{\dotproduct[{\gradient[{w}][{[
	{\cdots[\dotproduct[{\gradient[{w}][{[{\dotproduct[{\gradient[{w}][{[u(x_{i})]}]}][{\hat{x}^{\mu_{0}}}] }]}]}][{\hat{x}^{\mu_{1}}}]]\cdots}]}]}][{\hat{x}^{\mu_{q-2}}}]}]}]}][{\hat{x}^{\mu_{q-1}}}]] \label{eq:diff_p}
\end{align}

It will now be shown, that under specific constraints on the weight functions, the above partial derivatives are, to leading order, the corresponding partial derivatives of differential calculus. However certain properties of differential calculus, such as the commutativity of mixed partial derivatives, do not automatically hold to all orders of approximation with all possible choices of edge weights. Fortunately, with weights that decay suitably, issues of commutativity and convergence can be guaranteed.

\subsection{Constraints on edge weights} \label{sec:graphtheory_weightconstraints}
When using the given non-local calculus definitions, there is an arbitrary normalization of the quantities that can be chosen. For numerical purposes, and to ensure the operator definitions are consistent and independent of the size of the graph, as discussed, all integrations over the graph are normalized by the size of the discrete space we are in $\abs{V} = n$. In what follows repeated Greek indices $\mu \in \{ 0,\dots p-1\}$ will follow Einstein summation convention.

For the non-local calculus, we will observe order-by-order how the partial derivatives converge, given derived constraints on the weight function $w(x_{i},x_{j})$ for the $p$ dimensional representation of the state vector $x_{i}$. We begin by examining the normalization condition in \cref{eq:unitnorm} for this graph, with the normalization constant denoted as $W[p,n]$
\begin{align}
	\sum_{j \in V\setminus \{i\}} (x^{\mu}_{j}-x^{\mu}_{i})(x^{\nu}_{j}-x^{\nu}_{i})w(x_{i},x_{j}) = W^{\mu\nu}_{}[p,n] = (n-1)\delta^{\mu\nu} ~~~ \forall~ i. \label{eq:weightnorm}
\end{align}
Due to the quadratic form of the constraint, posed as a second order moment of this weight function, it is natural to assume that the weights are decaying, and take a form that decays at least as rapidly the inverse square of the distance between graph vertices:
\begin{align}
	w(x_{i},x_{j}) = O\left(\frac{1}{\abs{x_{j}-x_{i}}^2}\right).
\end{align}
For $p=1$, the form of the weight may be chosen as precisely the relation above, and for $x = x^{1}$,
\begin{align}
	w(x_{i},x_{j}) = \frac{1}{\abs{x_{j}-x_{i}}^2}.
\end{align}
For $p>1$ dimensional graphs, the form of the weight is less obvious if the Kronecker delta representation is to be extended or generalized. The weight function must decay rapidly enough to minimize consistency error and ensure that off-diagonal terms on the left-hand side of \cref{eq:weightnorm} are zero; however, too rapid a decay may fail to satisfy the finite diagonal constraints.

In the case of radially symmetric metrics $w({r_{ij}})$, where ${r_{ij}}^2 = \sum_{\mu} \abs{x^{\mu}_{j}-x^{\mu}_{i}}^2$, the normalization constraint can be written as
\begin{align}
	\sum_{j \in V\setminus \{i\}} {r_{ij}}^2 w({r_{ij}}) = W[p,n] = (n-1)p ~~~ \forall~ i,
	\label{eq:discreteweightnormalization}
\end{align}
which follows by computing the trace in \cref{eq:weightnorm}.

In generalizing \cref{eq:weightnorm} to the case of continuous manifolds, we consider the case of graphs whose states $x_{i}$ have support over a subset of $\reals[p]$, where the discrete sums approach integrals. Here we desire a form for the weights to guarantee that all off-diagonal terms, $\mu \neq \nu$ are identically zero. In this setting, \cref{eq:discreteweightnormalization} leads to the form
\begin{align}
	\int dr~ r^{p+1} w(r) = W[p,R] = p\frac{V_p}{S_p}R^p = R^p, \label{eq:weightrad}
\end{align}
where $V_p$ and $S_p$ are the volume and surface area of the unit ball in $p$ dimensions, $S_p = pV_p$, and $R$ is the radius of the sphere that encloses the graph.

In this continuous case of a graph whose states $x_{i}$ are dense in $\mathbb{R}^p$, the off-diagonal terms in \cref{eq:weightnorm} can be shown to cancel due to symmetry considerations, since $w(r)$ is an even symmetric function. Letting $x^{\mu} = rf^{\mu}(\theta)$, where $f$ are angular functions, the spherical symmetry ensures that only diagonal terms are non-zero in the tensor components
\begin{align}
	W^{\mu\nu}[p,R] &\equiv \int d^{p}y~ (y^{\mu} - x^{\mu})(y^{\nu} - x^{\nu})w(x,y).
\end{align}
For example in $p=2$ dimensions, where $x^{1} = r\cos{\theta}, x^{2} = r\sin{\theta}$, the angular integral trivially evaluates to zero.

However, in the case of graphs that are not dense in $\reals[p]$, with radius
\begin{align}
	R = \max_{ij} \frac{{r_{ij}}}{2}
\end{align}
in order for these off-diagonal terms to cancel in these finite integrals for an arbitrary graph structure, in addition to being symmetric functions, the weight functions must also decay rapidly so that the odd, antisymmetric terms are negligible.

\noindent It will now be assumed that for $p>1$, the weights decay faster than the above assumed inverse square relationship:
\begin{align}
	w(r) = o\left(\frac{1}{r^2}\right).
\end{align}
There still remains some choice in the exact form of the weights, and for this work, we will consider radially symmetric weights that decay like a Gaussian function $\sim e^{-r^2}$, or decay polynomially $\sim 1/r^2$. Based on these forms of the weights, there are some free parameters that can be calculated or chosen. We choose the weights to have the form $w(r) = C ~\tilde{w}(r/\sigma)$, where $C(p,R)$ is the weight scale, and $\sigma(p,R)$ is the weight decay parameter, or an appropriate length scale in the graph manifold. Based on \cref{eq:discreteweightnormalization}, the weights must scale as the inverse of length-squared, and so we will discuss the weight properties in terms of the scaled weight constraints
\begin{align}
\widetilde{W}[p,z] = \frac{W[p,R]}{C\sigma^{p+2}}, \label{eq:scaled_weight_constraint}
\end{align}
where $z = {R}/{\sigma}$ is a scaled radius. The weight constraints therefore take the form of $p$-order moments of the weight function:
\begin{align}
	\int_0^z dy~ {y}^{p+1} \tilde{w}(y) = \widetilde{W}[p,z]. \label{eq:constraint_scaled}
\end{align}

\subsubsection{Gaussian weights}\label{sec:graphtheory_weightconstraints_gaussian}
We first may consider Gaussian decaying weights
\begin{align}
	w(r) = Ce^{-\frac{1}{2}\frac{r^2}{\sigma^2}}, \label{eq:weight_gauss}
\end{align}
where in the case of graphs with a finite radius $R$, the weight is a truncated Gaussian for $0\leq r\leq R$. The scaled weight constraints are the $p+1$ order moments of the Gaussian distribution. These moments have a known recursion relationship \cite{Orjebin2014} for $p>0$:
\begin{align}
	\widetilde{W}[p,z] ~=&~ p\widetilde{W}[p-2,z] - z^pe^{-\frac{1}{2}z^2}, \label{eq:recurrence_gauss}
\end{align}
which also can be expressed in the analytic form
\begin{align}
	\widetilde{W}[p,z] ~=&~ p!!\left[\widetilde{W}[\nmodm{p}{2},z] - z^{p}e^{-\frac{1}{2}z^2}\sum_{i=0}^{\floor{p/2}-1} \frac{z^{-2i}}{(p-2i)!!}\right].
\end{align}

\subsubsection{Polynomial weights}\label{sec:graphtheory_weightconstraints_polynomial}
We also may consider polynomially decaying weights
\begin{align}
	w(r) = C\left(\frac{\sigma}{r}\right)^{2+\epsilon}, \label{eq:weight_poly}
\end{align}
where $0\leq\epsilon<p$ is the super-quadratic decay scaling. As previously discussed, the weights must decay faster than quadratically, however cannot decay so fast as for the integrals to vanish over discrete spherical volumes. Here the scaled weight constraints are tractable integrals for $p>0$:
\begin{align}
	\widetilde{W}[p,z] = \frac{1}{p-\epsilon} z^{p-\epsilon}.	\label{eq:constraint_poly}
\end{align}
Values for $\widetilde{W}[p,z]$ for radially symmetric Gaussian and polynomially decaying weights, with a finite radius, continuous graph are shown in \cref{tab:weights}. Here, $\Phi(x)$ is the standard normal cumulative distribution function.
\begin{table}[htb!]
	\centering
	\begin{tabular}{|c|c|c|} \hline\tableheight{20pt}
		$p$ & Gaussian & Polynomial \\ \hline\tableheight{20pt}
		$-2$ & $0$ & Undefined \\ \hline\tableheight{20pt}
		$-1$ & $\sqrt{2\pi}\left[\Phi(z) - \Phi(0)\right]$ & Undefined\\ \hline\tableheight{20pt}
		$0$ & $1 - e^{-\frac{1}{2}z^2}$ & Undefined \\ \hline\tableheight{20pt}
		$1$ & $\sqrt{2\pi}[\Phi(z) - \Phi(0)] - ze^{-\frac{1}{2}z^2}$ & $\frac{1}{1-\epsilon} z^{1-\epsilon}$ \\ \hline\tableheight{20pt}
		$2$ & $2 - (2+z^2)e^{-\frac{1}{2}z^2}$ & $\frac{1}{2-\epsilon} z^{2-\epsilon}$ \\ \hline

	\end{tabular}
	\caption{Weight constraint values for a radially symmetric weight $\widetilde{W}[p,z]$, given a finite radius $z = {R}/{\sigma}$, continuous graph in $p$ dimensions.}
	\label{tab:weights}
\end{table}

\subsubsection{Free parameters in weight functions}\label{sec:graphtheory_weightconstraints_freeparameters}

The free parameters $C$ and $\sigma$, and other hyperparameters such as the scaling $\epsilon$, may now be constrained based on for example \cref{eq:constraint_scaled} and \cref{eq:constraint_poly}, and there is some freedom, depending on the graph of interest. Generally, the decay parameter of the weight $\sigma$ can be chosen, after fixing the scale
\begin{align}
 C ~=&~ \frac{z^p}{\sigma^2}\frac{1}{\widetilde{W}[p,z]}, \label{eq:weightscale}
\end{align}
where the scaling of the weight by $1/\sigma^2$ is evident.

\begin{figure}[hpt]
\centering
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=0.6\textwidth]{figures/Wtilde_p_sigma_gauss.pdf}

	\subcaption{Gaussian}
	\label{fig:Wtilde_p_sigma_gauss}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=0.6\textwidth]{figures/Wtilde_p_sigma_poly.pdf}

	\subcaption{Polynomial}
	\label{fig:Wtilde_p_sigma_poly}
\end{subfigure}
\caption{Scaled weight constraints for various $p$ values, as a function of weight decay parameter $\sigma$; $\epsilon = {p}/{2}$ for the polynomial weights.}
\label{fig:Wtilde_p_sigma}
\end{figure}

\begin{figure}[hpt]
\centering
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=0.6\textwidth]{figures/C_p_sigma_gauss.pdf}

	\subcaption{Gaussian}
	\label{fig:C_p_sigma_gauss}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=0.6\textwidth]{figures/C_p_sigma_poly.pdf}

	\subcaption{Polynomial}
	\label{fig:C_p_sigma_poly}
\end{subfigure}
\caption{Weight scales for various $p$ values, as a function of weight decay parameter $\sigma$; $\epsilon = {p}/{2}$ for the polynomial weights.}
\label{fig:C_p_sigma}
\end{figure}

\begin{figure}[hpt]
\centering
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[trim={0 2 0 0},clip,width=0.8\textwidth]{figures/w_p_sigma_gauss.pdf}
	\vspace{-1.5cm}
	\subcaption{Gaussian}
	\label{fig:w_p_sigma_gauss}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[trim={0 -10pt 0 0},clip,width=0.8\textwidth]{figures/w_p_sigma_poly.pdf}
	\vspace{-1.5cm}
	\subcaption{Polynomial}
	\label{fig:w_p_sigma_poly}
\end{subfigure}
\vspace{-0.5cm}
\caption{Weight functions for various $p$ values, as a function of weight decay parameter $\sigma$ and radius $r$; $\epsilon = {p}/{2}$ for the polynomial weights.}
\label{fig:w_p_sigma}
\end{figure}
Using Gaussian and polynomially decaying weights, the calculated weight constraints $\widetilde{W}[p,\sigma]$, and the weight scalings $C[p,\sigma]$ are shown in \cref{fig:Wtilde_p_sigma} and \cref{fig:C_p_sigma}. For the polynomial weights, $\epsilon = {p}/{2}$. Here it can be seen that the Gaussian weight constraints are constant for small $\sigma > 0$, and then decrease for $\sigma > {R}/{z_{\textrm{inflection}}}$, where $z_{\textrm{inflection}}^2 = p+{3}/{2} + \sqrt{2p+{9}/{2}}$ is the point at which the Gaussian case of the integrand in \cref{eq:constraint_scaled} changes curvature and tends to zero. The Gaussian weight scalings also decrease polynomially in $\sigma$, as per \cref{eq:weightscale}. The polynomial weight constraints decay as ${1}/{\sigma^{p-\epsilon}}$, and the polynomial weight scalings decay as ${1}/{\sigma^{2+\epsilon}}$.

Importantly, for a fixed $\epsilon$, the polynomial weight scalings scale much slower with $\sigma$ and $p$ compared to the Gaussian weight scalings due to this polynomial form. For $\epsilon=0$ and a given $p$ value, how the polynomial scalings $C$ vary with $\sigma$ exactly cancel how the weight itself varies with $\sigma$. Therefore the polynomial weight function is constant with respect to $\sigma$:
\begin{align}
	w(r) = (p-\epsilon)\frac{R^{\epsilon}}{r^{2+\epsilon}}, \label{eq:poly_weight_full}
\end{align}
and suggests this weight form is more easily applied to graphs with different weight scales. In contrast, the Gaussian weights scale as ${1}/{\sigma^{2}} e^{{R^2}/{2\sigma^2}}$ with respect to $\sigma$ and increases exponentially for $\sigma \ll R$. Both of these calculations show the weight values also increasing with $p$, indicating that functions that depend on a large number of variables may still have very large weights, regardless of the $\sigma$ dependence.

We note that the Gaussian weights bear connections to the diffusion kernel, which furnishes the probability distribution between states on a graph in the context of discovering the data embedding in a continuous manifold. Currently the form of the weights also has a single length scale $\sigma$, and the manifold is assumed to be homogeneous along each state variable dimension. Future communications will investigate the effects of using each independent equation along each dimension in \cref{eq:discreteweightnormalization} to determine constraints on length scales in each direction $\sigma^{\mu}$. For the case of discrete manifolds, certain approximations on these continuous weight constraints must be made. Please refer to \cref{app:implementation_gauss} of the Appendix for implementation details of these approximations.

\subsection{Consistency of non-local derivatives} \label{sec:graphtheory_consistency}
With this insight into the form of the weights, and their constraints, we now examine the consistency of the partial derivatives in this non-local calculus. The analysis will follow a similar approach to any numerical scheme, substituting Taylor series of differential calculus terms for various non-local terms, to determine the leading order behavior. Due to the weight constraint, we will see that depending on the decay rate of the weights, the leading order error is potentially first order in the spacing of the vertices.

\noindent For the first order derivatives the weight constraint yields
\begin{align}
	\difference[1]{u(x_{i})}{x^{\mu}}
	~=&~ \derivative[1]{u(x_{i})}{x^{\mu}} + O\left(\derivative[2]{u(x_{i})}{x^{\nu},x^{\theta}}(x^{\nu}_{j}-x^{\nu}_{i})(x^{\theta}_{j}-x^{\theta}_{i})(x^{\mu}_{j}-x^{\mu}_{i}) w(x_{i},x_{j})\right). \label{eq:non-local1}
\end{align}
For second order derivatives, each non-local first derivative can be expanded in the form found above in \cref{eq:non-local1}, yielding
\begin{align}
	\difference[2]{u(x_{i})}{x^{\mu},x^{\nu}}
	~=&~ \derivative[2]{u(x_{i})}{x^{\mu},x^{\nu}} + O\left(\derivative[2]{u(x_{i})}{x^{\theta},x^{\epsilon}} (x^{\theta}_{j}-x^{\theta}_{k})(x^{\epsilon}_{j}-x^{\epsilon}_{k})(x^{\mu}_{j}-x^{\mu}_{k})(x^{\nu}_{k}-x^{\nu}_{i}) w(x_{k},x_{j})w(x_{i},x_{k})\right).\label{eq:non-local2}
\end{align}
Therefore, given this weight constraint, the first order non-local derivatives to leading order correspond with the corresponding local derivatives. Although the second order non-local derivatives do not contain first order terms, the second order terms contain the other second order derivatives along the other dimensions. It can be concluded that in order for the non-local calculus definitions for the partial derivatives to equal the local calculus definitions to leading order, a constraint on the higher order moments of the weight function for each order of derivative is necessary. These constraints are assumed to be of the form of products of all differences in state variables at all combinations of dimension and pairs of vertices, weighted by the edge weight between those vertices:
\begin{align}
\ndifference[q]{u(x_{i})}{x^{\mu_{0}}}{x^{\mu_{q-1}}}  = \nderivative[q]{u(x_{i})}{x^{\mu_{0}}}{x^{\mu_{q-1}}} + O(x_{j}-x_{i}) \longleftrightarrow & \\
\sum_{j_{0}\dots j_{q-1} \in V\setminus \{i\}} (x^{\mu_{q-1}}_{j_{q-1}}-x^{\mu_{q-1}}_{j_{q-2}}) \cdots (x^{\mu_{0}}_{j_{0}}-x^{\mu_{0}}_{i}) w(x_{j_{q-2}},x_{j_{q-1}}) \cdots w(x_{i},x_{j_{0}}) ~=&~  \nonumber \\
(n-1)^q \delta^{\mu_{q-1}\mu_{q-2}}\cdots \delta^{\mu_{0}\mu_{i}}. \nonumber
\end{align}

As alluded to in the previous example of decaying Gaussian and polynomial weights, an alternative to choosing weights based on these difficult to construct higher order moment constraints, is to choose weights that decay faster than any finite polynomial, and still satisfy the primary constraint in \cref{eq:weightnorm}. Given a fast enough decaying weight, this will allow several issues to be resolved. First, all off-diagonal terms in the primary constraint will decay rapidly. Second, all higher order non-local derivatives will equal their differential derivative counterparts to leading order due to the weights decaying faster than the polynomials in the Taylor series.

\subsubsection{Compromises with a non-local calculus approach}\label{sec:graphtheory_consistency_compromise}

However, due to this calculus involving sums over the entire graph, the non-local nature means the terms beyond the local calculus derivatives in the consistency relationships above in \cref{eq:non-local1,eq:non-local2} may be unreasonably large. This is particularly true if the weight scales are large as we have seen for the Gaussian weights. Therefore there may be no real similarity in the values and trends of the non-local derivatives compared to their differential counterparts. However, from studies of the total scale of the weights, due to the polynomially decaying weights having a closed form that is independent of the decay parameter, these weights are reasonable scales of order $O(1)$ for small $p$.

Due to these numerical issues, we choose to use the polynomially decaying weights for the current studies, with $\epsilon = {p}/{2}$. Although not necessarily similar in form, the non-local derivative values will be similar in magnitude to what would be expected for differential derivatives, and decay reasonably quickly with $\epsilon>0$. Future communications will study in depth the effect of using different forms and parameter choices for the weights, and how these large non-local error terms can be understood and controlled.

This extends the compromise regarding the weight functions to be made, between obeying the consistency and constraint equations, and having scales of derivatives that correspond with those from differential calculus in the model. As to be discussed when these methods are applied to physical systems, the model fitting procedure may be able to overcome the large scales of these derivatives, such as associating small or large linear coefficients to terms in the model that contain these large values. Please refer to \cref{app:implementation_linearregression} of the Appendix for details on the regression implementation.

This non-local calculus has the objective of being a general, accurate framework for any modelling procedure, and the exact form of the non-local error term will be developed. An important aspect of numerical method development is an analysis of the effect of data sparsity on model error. As can be seen in \cref{eq:non-local1}, the non-local derivatives have a persistent error term for all weight functions with non-local extent, and whether this is detrimental to subsequent reduced order modelling is an important question in determining the validity of this approach.

For example, a modified Taylor series, as discussed in \cref{app:implementation_modifiedtaylorseries} of the Appendix, may be constructed using a training dataset as a functional representation, where a different Taylor series model is developed at each possible base point. Please refer to \cref{app:error_mesh} for details concerning the mesh of data used. The point-wise error for evaluation points within a distance $h$ of each base point model can then be computed, and refining the spacing of the data, we are able to get a sense of how the method error scales with $h$. A complete error analysis for this modified Taylor series approach can be referred to in \cref{app:error} for $p=1$, and in \cref{app:error_p} for higher dimensions. As discussed in \cref{app:error_pointwisederivative,app:error_pointwisgamma,app:error_totalerror} of the Appendix, in addition to the expected effects of a $\ith[k]$ order Taylor series having error that scales as $h^{k+1}$, the non-local calculus derivatives, and any fit model parameters are also dependent on the spatial distribution of the training and testing data.

With this understanding of the behaviors of the non-local calculus and partial derivatives, reduced order models can now be developed for systems of interest. The modelling will consist of finding a basis of operators for a given model, where the operators are defined in terms of this non-local calculus. In the following section, an example of a physical systems is presented, and several reduced order modelling approaches are examined.

\section{Physical systems of interest} \label{sec:physicalsystems}
\subsection{Microstructures in a gradient-regularized model of non-convex elasticity} \label{sec:physicalsystems_microstructures}
We apply the graph theoretic approach to studying reduced order models for the mechanochemical response of solids that develop distinct regions of local phases over their evolution. This physical system is driven by a free energy density function that is parameterized by composition (chemical) and strain (elastic) variables. The underlying functional form of the free energy is non-convex in both these quantities and is regularized by composition and strain gradients. Microstructures develop as phases and symmetry-breaking structural variants arise corresponding to negative eigenvalues of the Hessian of the free energy density in the strain-composition space. These models are ill-posed, however; a condition arising from the absence of penalization on inter-phase and inter-variant interfaces and manifesting in PDE solutions that are pathologically discretization-dependent. Well-posedness and penalization of interfaces are restored by free energy functionals $\psi_{\textrm{tot}} = \psi_{\textrm{hom}} + \psi_{\textrm{grad}}$, dependent on the composition and strain gradients. When decoupled, the chemical component of this problem causes spinodal decomposition and Ostwald ripening, described by the Cahn-Hilliard equation \cite{Cahn1958}. The non-linear elasticity problem is described by a non-convex strain energy density stabilized by incorporating strain gradient effects as described by many authors.

Our treatment follows Toupin \cite{Toupin1962} and has appeared as decoupled, gradient-regularization of non-convex, non-linear elasticity \cite{Rudraraju2014,Teichertetal2017}, as well as mechanochemical spinodal decomposition \cite{Rudraraju2016,Sagiyamaetal2016,SagiyamaGarikipati2017b}. The free energy density functional is $\psi= \psi_{\textrm{tot}}$, with a dependence on the fields $\{c,\vec{E},\gradient[]{c},\gradient[]{\vec{E}}\}$ where $c$ is the composition and $\vec{E}$ is the Green-Lagrange strain tensor. This dependence, however, is expressed in terms of symmetry-coding strain parameters that are introduced below.

At low temperatures, these systems consist of stable phases related to the local crystalline symmetry, specifically cubic, or tetragonal symmetry. The system considered is restricted to $d=2$ spatial dimensions, and so the material will consist of phases with either square or rectangular structure. This pattern of regions of phases of a field are generally referred to as microstructures. This system is further modeled as a binary mixture, with scalar composition $0\leq c \leq 1$ of one of the two constituents.

Additionally, the introduction of mechanical strain couplings to the composition allows for a degeneracy in the rectangular phase. Linear combinations of the strain components $\vec{E}$ yield relevant parameters $\vec{e}$ relating to the original high symmetry of the square phase, namely $e_{1,2} \equiv ({E_{11} \pm E_{22}})/{\sqrt{2}}, e_{6} \equiv \sqrt{2} E_{12}$, that allow this degeneracy to be determined. Specifically, the sign of the $e_{2}$ field can determine the orientation of the rectangular phase, either in the horizontal "plus", or the vertical "minus" variant. The order parameters will thus be defined to be $\eta = (c,e_{2})$, and the microstructure can be described by its phase volume fractions $\varphi_{\alpha}$, where $\alpha = \{\square,\hrectangle,\vrectangle\}$ corresponds to the phase of the material.

\subsubsection{Free energy description} \label{sec:physicalsystems_microstructures_freeenergy}
To describe these phase transitions, we take the description used by Zhang \& Garikipati\cite{Zhang2020}. The total free energy density description of the system $\psi$ can be comprised of a homogeneous component, $\psi_{\textrm{hom}}(c,\vec{e})$, that is not entirely concave up over the entire $c-e_{2}$ domain, yielding potential wells with stable phases. An example two dimensional energy landscape from Rudraraju \emph{et al.}~ \cite{Rudraraju2016} is shown in \cref{fig:energyphase}.
\begin{figure}[ht]
	\centering

	\includegraphics[width=0.6\textwidth]{energy_curvature_mechanochem.pdf}

	\caption{Free energy density landscape over the composition $c$ and strain component $e_{2}$ order parameters, showing changes in curvature and resulting phases of square and rectangular lattice symmetries, from Rudraraju \emph{et al.}~ \cite{Rudraraju2016}. On the right, the effect of temperature and composition on the phase diagram of stable regions of purely square or rectangular phases, or a mixture of phases is shown.}
	\label{fig:energyphase}
\end{figure}

During the mechanical-chemical coupled dynamical processes, within the domain where the free energy is concave, such as around the saddle point $B$ in \cref{fig:energyphase}, phase separation will tend to occur, where a mixture of the stable phases, will yield a lower energetic state. This can be seen from a typical lever law analysis for a binary mixture, to obtain the phase volume fractions $\varphi_{\alpha}$. This simple phase separation within the non-convex region of the free energy density is mechano-chemical spinodal decomposition, where phases tend to separate. Furthermore, the gradient regularization that represents the energy of interfaces also manifests in a number $N_{\alpha}$ of each variant, and total interfacial lengths $l_{\alpha}$.

With the composition $c$ and symmetry-coding strain parameter $\vec{e}$ and their gradients, a non-uniform free energy density component $\psi_{\textrm{grad}}(c,\vec{e},\gradient[]{c}, \gradient[]{\vec{e}})$ allows for gradient regularization, typically through (coupled) quadratic terms in the order parameter gradients. These terms yield phase regions that are limited by interfaces being less energetically favorable, and so equilibrium states will tend to have less numbers of domains of phases, and furthermore, larger phase regions grow at the expense of smaller ones. This effect is known as Oswald ripening \cite{Rudraraju2016}.

This gradient contribution, in addition to the homogeneous contribution yields the total free energy functional of the order parameters over the material volume $\Omega_t$ at time $t$,
\begin{equation}
	\Psi[c,\vec{e}] = \int_{\Omega_t} dV~ (\psi_{\textrm{hom}}+\psi_{\textrm{grad}}) + \int_{\partial \Omega_t} dS~ \dotproduct[{\vec{T}}][{\vec{u}}], \label{eq:Psi_tot}
\end{equation}
where $\vec{T}$ are mechanical traction terms, depending on which boundary conditions are specified, and $\vec{u}$ is the spatial displacement field of the material.

A variational formalism $\difference[1]{\Psi}{{c,\vec{e}}}$ can then be applied to yield the relevant generalized chemical potentials. As discussed, for this particular work, the homogeneous free energy density is a smooth polynomial function of the order parameters,
\begin{equation}
	\psi_{\textrm{hom}}(c,\vec{e}) \equiv 16 \alpha_{c} c^4 - 32\alpha_{c} c^3 + \alpha_{c} c^2 + 2\frac{\alpha_{e}}{{\beta_{e}}^2}({e_{1}}^2 + {e_{6}}^2) + \frac{\alpha_{e}}{{\beta_{e}}^2}{e_{2}}^4 + 2\frac{\alpha_{e}}{{\beta_{e}}^2}(2c-1){e_{2}}^2,
	\label{eq:psi_hom}
\end{equation}
where $\{\alpha_{c},\alpha_{e},\beta_{e}\}$ are constant coefficients.

The gradient contributions will take the form of strictly quadratic forms in the order parameter gradients,
\begin{equation}
	\psi_{\textrm{grad}}(\gradient[]{c},\gradient[]{\vec{e}}) = \frac{1}{2} \gradient[]{c} \cdot\vec{\kappa} \gradient[]{c} + \frac{1}{2} \gradient[]{\vec{e}} \cdot \boldsymbol{\gamma} \gradient[]{\vec{e}} + \gradient[]{c}\cdot \boldsymbol{\theta} \gradient[]{\vec{e}}, \label{eq:psi_grad}
\end{equation}
and for this work, it is assumed that there are no gradient couplings between the composition and strain fields $\boldsymbol{\theta} = 0$, and only constant isotropic tensors $\vec{\kappa} = \kappa \textrm{I}$ and $\boldsymbol{\gamma} = \gamma \vec{I}$ for the gradient terms. Furthermore, the additive terms mean the total free energy density can be separated into the sum of independent mechanical and chemical terms, plus a coupling term:
\begin{equation}
	\psi= \psi_{\textrm{mech}} + \psi_{\textrm{chem}} + \psi_{\textrm{couple}}.
\end{equation}

When stationarity of the free energy is imposed in a variational framework, the Euler-Lagrange equations yield the weak form of the balance of linear momentum, and further variational arguments lead to the strong form. Using this variational approach, based on this form of the free energy, and given its dependence on the strain and composition fields, the governing PDE's are chosen to be coupled Cahn-Hilliard dynamics for the composition \cite{Cahn1958}, and non-linear gradient elasticity equations in equilibrium \cite{Toupin1962,Rudraraju2016}.

\subsubsection{Cahn-Hilliard dynamics} \label{sec:physicalsystems_microstructures_CahnHilliard}
Cahn-Hilliard dynamics are first order for the composition and take the form of a transport equation
\begin{align}
	\derivative[1]{c}{t} + \divergence[][{\vec{J}}] ~=&~ 0, \label{eq:pde_cahnhilliard}
\end{align}
where
\begin{align}
	\vec{J} ~=&~ -\vec{L}\gradient[][\mu], \label{eq:cflux}
\end{align}
and $\vec{L} = \vec{L}(c,\vec{e})$ is a mobility transport tensor, and is assumed to be isotropic $\vec{L} = \textrm{L}\vec{I}$. The variational treatment yields the chemical potential $\mu = \difference[1]{\Psi}{c}$, whose gradient yields the flux in a form that is guided by the thermodynamic dissipation inequality \cite{degrootmazur1984}, and leads to there being fourth-order mechano-chemically coupled equations.

\subsubsection{Toupin model of elasticity} \label{sec:physicalsystems_microstructures_Toupin}
The strain gradient elasticity can be modeled in equilibrium as equations for the first Piola-Kirchhoff $\vec{P}$ and higher order $\vec{B}$ stress tensors,
\begin{align}
	\vec{P} ~=&~~ \derivative[1]{\psi}{\vec{F}} \label{eq:piola1}\\
	\vec{B} ~=&~~ \derivative[1]{\psi}{\gradient[]{\vec{F}}} \label{eq:piola3}
\end{align}
which are derivatives of the strain energy and are conjugate to the deformation gradient $\vec{F}$. Here the Green-Lagrange strains are $\vec{E} = {1}/{2}\left(\vec{F}^{T}\vec{F} - \vec{I}\right)$ and gradients $\gradient[]$ are derivatives with respect to the material indices. The elasticity governing equations, shown here in the strong form, are
\begin{align}
	\divergence[]{\vec{P}} - \divergence[]{(\divergence[]{\vec{B}})} = 0, \label{eq:pde_elasticity}
\end{align}
plus additional Dirichlet and higher order traction boundary conditions.

\subsubsection{Direct numerical simulations} \label{sec:physicalsystems_microstructures_DNS}
These governing PDEs in \cref{eq:pde_cahnhilliard,eq:pde_elasticity} may now be solved, yielding time series data, as per Zhang \& Garikipati \cite{Zhang2020}. The weak form is chosen here for implementation, motivated by numerical methods such as finite elements and the isogeometric approaches. \cite{Rudraraju2014,Rudraraju2016} A given system is described by the static mechanical Dirichlet boundary conditions for the displacement fields $\vec{u}$, uniform initial $e_{2} = 0$ fields, corresponding to the material being initially in the high temperature square phase, random initial compositions $c = 0.46 \pm 0.05$, and zero chemical flux boundaries. The DNS proceeds by time stepping with steps of $\Delta t = 10^{-5}~\textrm{s}$, and iteratively solving the PDE's at each time step, with tolerances placed on the residual of the solver. The data generated for this work is at times $50-1500 ~\mu\textrm{s}$, with example DNS observables shown in \cref{fig:DNSobservables}. In this work, we consider a single boundary and initial condition and one evolution of the microstructure and aim to model its quantities of interest. In future communications we will consider the effectiveness of the graph theoretic approach at producing general models for a family of microstructures with different boundary and initial conditions.

To develop the graph theoretic model for the given mechano-chemical system, the state vector is assembled from quantities of interest that are functionals of the high-dimensional solution to \cref{eq:pde_cahnhilliard,eq:pde_elasticity}. These quantities, denoted as barred variables unless otherwise specified, are effective volume averaged quantities $\bar{f}$ over the material domain. They therefore contain no spatial dependence of their high dimensional counterparts $f(\eta(x))$. Chemical quantities such as the local composition and chemical potential are represented by the phase volume fractions, interfacial lengths, and numbers of phase domains
\begin{align}
\{c(x),~\mu(x)\} ~\to&~ \{\varphi_{\alpha},~l_{\alpha},~N_{\alpha} \}_{\alpha =\{\square,\hrectangle,\vrectangle,\dots~\}}, \nonumber \\
\intertext{mechanical quantities such as the deformation gradients, Green-Langrange strains, and Piola-Kirchoff stresses are represented by their volume averaged quantities}
\{\vec{F}(x),~\vec{E}(x),~\vec{P}(x),\dots\} ~\to&~ \{\vecbar{F},~\vecbar{E},~\vecbar{P},\dots\}, \nonumber \\
\intertext{and free energy quantities are represented by their total values}
\{\psi(x),~\psi_{\textrm{mech}}(x),~\psi_{\textrm{chem}}(x),\dots\} ~\to&~ \{\Psi,~\Psi_{\textrm{mech}},~\Psi_{\textrm{chem}},\dots \}.
\end{align}
The phase volume fractions $\varphi$ are calculated by summing over volume elements that are within the associated phase, as per the phase diagram in \cref{fig:energyphase}; the interfacial lengths $l_{\alpha}$ are calculated by computing the total lengths of the interfaces between the square-rectangular-(plus,minus) fields; the number of domains $N_{\alpha}$ are counted through a gradient edge detection algorithm to identify the number of distinct domains of each field. The volume like terms of $l_{\alpha}$ and $N_{\alpha}$ are essential for distinguishing changes in topology of the microstructure, such as when field domains merge, separate, appear, or disappear. The specific forms of the free energy $\Psi, \Psi_{\textrm{mech}}$ are the integrated free energy densities over the material volume. For the strain quantities, they are chosen to be defined in terms of the volume averaged Green-Lagrange strains
\begin{align}
	\vecbar{E} = \frac{1}{2}\left(\vecbar{F}^{T}\vecbar{F} - \vec{I}\right),
\end{align}
where it is the deformation gradient $\vec{F}$ that has been explicitly volume averaged
\begin{align}
\vecbar{F} = \int_{\Omega_t} dV~ \vec{F}.
\end{align}

Due to the high dimensional data being generated by a finite element DNS method, each of the spatially dependent quantities can be thought of as fluctuating about their volume averaged value, and could be represented by being expanded as a Taylor series
\begin{align}
f(\eta) = \bar{f} + \epsilon + \dotproduct[{\derivative[1]{f(\tilde{\eta})}{\eta}}][{(\eta-\tilde{\eta})}] + O((\eta-\tilde{\eta})^2) \label{eq:local_avg_quantities}
\end{align}
in the order parameters, where $\tilde{\eta} = \argmin_{\eta} \abs{f(\eta)-\bar{f}}$, such that $f(\tilde{\eta}) = \bar{f} + \epsilon$. This polynomial series means any dependencies of the spatially dependent quantities on the order parameters will become higher order polynomial dependencies in the volumed averaged quantities; $\bar{f}$ will be a higher order polynomial in $\bar{\eta}$. For example, the free energy density contains up to fourth order polynomials of the order parameters, however the total free energy will be a function of much higher polynomials of the effective order parameter quantities. These effects will be discussed when determining forms of the reduced order models while being informed of the original dependencies that yield the training data.

The relevant $p$-dimensional state vector $x = \{x^{\mu}\}$ associated with each vertex $i$ of the graph is chosen to be
\begin{equation}
	x = \{\Psi,\Psi_{\textrm{mech}} ,\vecbar{E},\varphi_{\alpha},l_{\alpha},N_{\alpha}\}, \label{eq:observables}
\end{equation}
and vertices in this graph are associated with the time index. From physical arguments, the microstructure states at different times are all related, and therefore the graph is considered to be fully connected when considering edges and edge weights.

\begin{figure}[hpt]
\centering
\includegraphics[width=0.6\textwidth]{figures/observables_DNS_updated.pdf}
\caption{DNS calculated state vector components.}
\label{fig:DNSobservables}
\end{figure}

\subsection{Reduced order modelling} \label{sec:physicalsystems_reducedordermodelling}
To determine a model for the system, two different aspects of the model are considered, the time dependence of the phase volume fractions, when zero additional mechanical stress is applied to the system, and the free energy dependence on the strain and phase volume fractions. The forms of these reduced order models will follow the general principals of how the various fields are coupled from this variational framework. This choice of PDE's involves polynomial couplings between the composition and strain fields in the free energy density in \cref{eq:psi_hom}, however the constitutive relations for the chemical potential do not involve explicit derivatives of the free energy density with respect to the strains or strain gradients, and the constitutive relations for the elasticity do not involve derivatives of the free energy density with respect to the compositions or composition gradients. These polynomial and derivative couplings indicate important quantities of interest are those dominated by either algebraic or differential forms of the chemical and mechanical effects.

\subsubsection{Selection of physically relevant basis} \label{sec:physicalsystems_reducedordermodelling_Basis}
Given the microstructure data is generated from explicit forms for the free energy density in \cref{eq:psi_hom,eq:psi_grad}, these exact forms give an indication for which algebraic and differential terms should be proposed for the model. The free energy contains no explicit spatial dependence, and contains monomials up to $c^4$ and ${{e_{\alpha}}}^4$ order, plus coupled polynomial terms that are up to $c {e_{2}}^2$ order, plus quadratic gradient terms. The chemical potential and associated current are of the form
\begin{align}
	\mu ~=&~ \derivative[1]{\psi}{c} - \kappa \laplacian[][{c}] \\
	\divergence[][{\vec{J}}] ~=&~ -L\left[\uniderivative[3]{\psi}{c}\dotproduct[{\gradient[]{c}}][{\gradient[]{c}}] + 2\derivative[3]{\psi}{c,c,{e_{\alpha}}}\dotproduct[{\gradient[]{c}}][{\gradient[]{{e_{\alpha}}}}] + \derivative[3]{\psi}{c,{e_{\alpha}},{e_{\beta}}}\dotproduct[{\gradient[]{{e_{\alpha}}}}][{\gradient[]{{e_{\beta}}}}] \right.\\
	&~~~~~~~~~~~ + \left.\uniderivative[2]{\psi}{c}\laplacian[][{c}] +\derivative[2]{\psi}{c,{e_{\alpha}}}\laplacian[][{e_{\alpha}}] - \kappa \laplacian[][]\laplacian[][{c}]\right] \nonumber.
\end{align}
\noindent The equations are solved in the weak form over the material domain, and therefore any spatial dependencies, and particularly any spatial derivatives of the order parameters, are removed, yielding terms that are proportional to the parameter being differentiated, similar to \cref{eq:local_avg_quantities}. For example $\int_{\Omega_t} dV~ \nlaplacian[][{\eta}][n] \sim O(\eta)$. Therefore the Cahn-Hilliard first order dynamics for the composition should be of the order
\begin{align}
	\derivative[1]{c}{t} &\sim O(\uniderivative[3]{\psi}{c} c^2) + O(\derivative[3]{\psi}{c,c,{e_{\alpha}}} c {e_{\alpha}})
	+ O(\derivative[3]{\psi}{c,{{e_{\alpha}}},{{e_{\beta}}}}~ {{e_{\alpha}}}{e_{\beta}}) \label{eq:cahnhilliard_approx}	\\
	&~ + O(\uniderivative[2]{\psi}{c} c) + O(\derivative[2]{\psi}{c,{e_{\alpha}}} {e_{\alpha}}) + O(c) \nonumber.
\end{align}
For the elasticity equations, given the free energy and its derivatives, upon being integrated over the material volume are polynomials in the order parameters, it will be assumed that a functional representation of the free energy can be expressed as a finite polynomial series in the phase volume fractions and volume averaged strain parameters.

When performing this reduced order modelling of the volume averaged quantities, we will be informed by the various equations and relations used to produce this data in determining which operators should be present in the models. Both the specific terms present, the form of the higher dimensional representations, such as polynomial expansions, differential relationships, or general dependencies, as well as the number of terms required to form the model will inform the hypothesis of the model. Specific implementation details are discussed in \cref{app:implementation_stepwiseregression} of the Appendix, as well as in Wang \emph{et al.}~ \cite{Wang2019,Wang2020} and Kochunas \emph{et al.}~ \cite{Kochunas2020}.

Here, all partial derivatives are approximated using the non-local calculus definitions in \cref{eq:diff_p}, with the exception of the time derivatives, which are calculated using the standard backwards Euler finite difference method. It is also important to note when computing these operators and forming a basis, finite datasets mean numerically some terms may be not linearly independent. This will affect the span of this basis, and if there are more linearly dependent terms in the model, some models with different subsets of operators may be degenerate. This is an effect that will be discussed for the specific models below.

\subsubsection{Weighted residual loss functions} \label{sec:physicalsystems_reducedordermodelling_lossfuncs}
An important aspect of the fitting procedure is the choice of loss function to characterize the error between the known data and predictions. For general models, we choose the root mean square error with the $L_{2}$ norm. The $L_{2}$ norm loss generally minimizes fluctuations in the data and leads to smooth curves, whereas the $L_{1}$ norm loss generally minimizes error at specific points, and the $L_{\infty}$ norm bounds the maximum error of a single point. For highly oscillatory data, such as mechano-chemical systems which can change their local phase rapidly during their transient response, a weighted loss function with several norms is more appropriate. Here, we define the total loss function as a weighted sum of loss functions
\begin{align}
	\ell = \sum_{\chi \in \{1,2,\infty,\cdots\}} w_{\chi}\cdot \ell_{\chi}, \label{eq:lossfunc}
\end{align}
where $\chi$ denotes the type of norm $L_{\chi}$ used in the loss function, with weighting $w_{\chi}$. Each loss function $\ell_{\chi}$ may have its own sub-domain to evaluate the loss. In general, we choose the $\ell_{2}$ domain to be over the entire dataset. For oscillatory problems we choose the $\ell_{1}$ loss to be evaluated at the local maxima and minima of the data, through a specified number of passes of a peak finding algorithm.

\subsubsection{Regression approaches} \label{sec:physicalsystems_reducedordermodelling_regression}
The general modelling approach will consist of proposing a basis of algebraic and differential operators. System inference methods, specifically linear stepwise regression approaches, allow the model to be made as parsimonious as possible, while remaining an accurate representation of the underlying physics.

Here, we choose to perform regression using Ordinary Least Squares (OLS), as well as Ridge regression (Ridge)\cite{Lv2009} and Lasso regression (Lasso) \cite{Tibshirani1996}, schemes which constrain the norm of the coefficients. Regularization is added to the regression with the objective of improving the robustness of the method, particularly for low rank basis of operators that induce very high condition numbers in the linear problems to be solved. In the case of Ridge regression, it can be directly shown that regularization increases the rank and decreases the collinearity of the data. Please refer to \cref{app:implementation_linearregression_ridge} and \cref{app:implementation_linearregression_lasso} of the Appendix for details on these regression approaches and the residual equations being minimized.

The Ridge and Lasso regression approaches generally use a cross validation scheme to determine the optimal ratio $\lambda$ of the primary constraint and the additional norm constraints in the total residual loss function. At each iteration of the stepwise regression, we perform an independent cross validation to determine the optimal $\lambda$ hyper-parameter for each given reduced basis of operations with a proposed operator removed. This approach increases the computational complexity at least linearly with the size of the data, and linearly with the number of possible hyperparameters, however potentially ensures more optimal coefficients that are away from a local basin of the loss function. We observe that the Lasso regression induces sparsity in the regression coefficients due to using a $1$-norm constraint on the coefficients. This conflicts with the stepwise regression approach, and appears less suitable for these types of problems. We therefore will be comparing the OLS and Ridge approaches.

This effect of the independent cross validations yielding optimal coefficients for a given proposed model at a stepwise iteration, but not accounting for the overall stepwise procedure that is selecting which coefficients to remove from the model, will be seen to affect the overall trends of the model loss. There is potentially a competing interest between the penalization of the loss function with the coefficient norm constraints, and the desire to remove insignificant terms and produce the most effective model.

\textbf{Remark}: It is also important to acknowledge that the cross validation may be performed using a different optimality metric compared to the loss function for the statistical tests in the stepwise regression. The global minima of the two metrics do not necessarily coincide, and as will be discussed in \cref{sec:physicalsystems_vol}, this greatly affects the model, and the order of operators present after the overall stepwise regression is complete. Ultimately, we desire a robust, accurate, and efficient method, and the following example studies indicate the compromises to be made between these objectives.

\subsection{Phase volume fraction dynamics}\label{sec:physicalsystems_vol}
To model the phase volume fraction dynamics, we are guided by the Cahn-Hilliard equation in \cref{eq:pde_cahnhilliard}, and the order of terms present in \cref{eq:cahnhilliard_approx}. Here the composition evolves under first-order dynamics driven by a conservative flux that itself has a non-linear, and high-order gradient dependence on the composition and strains. Importantly, the flux is dependent on the chemical potential of $c$ and stresses $\vecbar{E}$ that arise as variational derivatives of the free energy with respect to the composition and strains, respectively in \cref{eq:cflux}. We therefore seek evolution equations for the volume phase volume fractions $\varphi$ that are also first-order in time, and dependent on polynomials of products of the free energy derivatives up to $\difference[3]{\Psi}{\varphi_{\alpha},\varphi_{\beta},\bar{E}_{\xi \kappa}}$, with phase dependent terms in $\{\varphi_{\alpha},l_{\alpha}, N_{\alpha},\vecbar{E}\}$. We further restrict the possible dependencies based on the form of the free energy in $d$ spatial dimensions, a function of a single chemical parameter and ${d(d+1)}/{2}$ mechanical parameters, which we choose to be $\bar{\varphi} = \varphi_{\square}$, and $\vecbar{E}$.

We also assume that any dependencies on the microstructure parameter terms $l_{\alpha},N_{\alpha}$ show up only in the terms that contain polynomials of the state variables, the last term in \cref{eq:cahnhilliard_approx}, and not in terms that contain derivatives. These polynomial terms without derivatives arise in the Cahn-Hilliard equation from the interfacial gradient terms in the free energy. The interfacial lengths and particle numbers, here defined to be $\bar{l}_{1} = l_{\square},~\bar{N}_{1} = N_{\square},~\bar{l}_{2} = {(l_{\vrectangle}+
l_{\vrectangle}-l_{\square})}/{2},~\bar{N}_{2} = N_{\hrectangle}+N_{\vrectangle}$ will suffice to model the $\kappa$ and $\gamma$ interfacial terms.

\subsubsection{Proposed model for phase volume fraction dynamics}\label{sec:physicalsystems_vol_proposedmodel}
The model for the phase volume fraction first order dynamics is assumed to be of the following form
\begin{equation}
	\derivative[1]{\bar{\varphi}}{t} = f(\nderivative[q]{\Psi}{\bar{\varphi}}{\bar{E}_{\alpha \beta}},\bar{\varphi},\bar{E}_{\alpha \beta},\bar{l}_{\alpha}, \bar{N}_{\alpha}), \label{eq:model_vol}
\end{equation}
where $\nderivative[q]{\Psi}{\bar{\varphi}}{\bar{E}_{\alpha \beta}}  = \{\derivative[2]{\Psi}{\bar{\varphi},\bar{E}_{\alpha \beta}},\derivative[2]{\Psi}{\bar{\varphi},\bar{\varphi}},\derivative[3]{\Psi}{\bar{\varphi},\bar{E}_{\alpha \beta},\bar{E}_{\xi \kappa}},\derivative[3]{\Psi}{\bar{\varphi},\bar{\varphi},\bar{E}_{\alpha \beta}},\uniderivative[3]{\Psi}{{\bar{\varphi}}}\}$ are the relevant derivatives for this model.

The function $f$ for this work is an algebraic, linear function of polynomials of the state vector components and derivatives of the free energy
\begin{equation}
	f = \mspace{-20mu}\sum_{{\substack{a = \left\{a_{\nderivative[q]{\Psi}{\bar{\varphi}}{\bar{E}_{\alpha \beta}}}\cdots~ a_{\bar{N}_{\beta}}\right\}\\\abs{a}\leq A\\ 0 \leq a_{\nderivative[q]{\Psi}{\bar{\varphi}}{\bar{E}_{\alpha \beta}}} \neq a_{\bar{l}_{\alpha}} \neq a_{\bar{N}_{\alpha}} \leq 1}}} \mspace{-40mu} \gamma^{a}~ ({\nderivative[q]{\Psi}{\bar{\varphi}}{\bar{E}_{\alpha \beta}})^{a_{\nderivative[q]{\Psi}{\bar{\varphi}}{\bar{E}_{\alpha \beta}}}}({\bar{\varphi}})^{a_{\bar{\varphi}}} ({\bar{E}_{\alpha \beta}})^{a_{\bar{E}_{\alpha \beta}}}({\bar{E}_{\xi \kappa}})^{a_{\bar{E}_{\xi \kappa}}}({\bar{l}_{\alpha}})^{a_{\bar{l}_{\alpha}}} ({\bar{N}_{\beta}})^{a_{\bar{N}_{\beta}}}} \label{eq:model_volfunc}
\end{equation}
where the powers of the quantities are represented by components of the multicomponent index $a$, with the indices being denoted by the quantities themselves. Here, summation over repeated indices in a given term should not be assumed. Instead, the function is explicitly a sum over all combinations of products of powers of the quantities, including all combinations of components of the tensor quantities. The powers are constrained such that the derivative $a_{\nderivativeinline[q]{\Psi}{\bar{\varphi}}{\bar{E}_{\alpha \beta}}}$, interfacial length $a_{\bar{l}_{\alpha}}$, and number $a_{\bar{N}_{\alpha}}$ terms are all separate from each other and never appear more than once in a term, as per \cref{eq:cahnhilliard_approx}. The polynomial terms are up to order $A=3$ powers, and there are a maximum of $308$ terms in this model.

\subsubsection{Numerical implementation of phase volume fraction dynamics}\label{sec:physicalsystems_vol_numerical}
Given this basis for the terms that drive the phase volume fraction dynamics, stepwise regression is performed using a mixture of the $\ell_{2}$ and $\ell_{1}$ loss functions, as per \cref{eq:lossfunc}. Various loss weights are chosen and a single pass of a peak finding algorithm is used to find the points at which the $\ell_{1}$ loss is computed. The resulting loss curves with a comparison of OLS and Ridge regression fitting are shown in \cref{fig:vol_loss_OLS_Ridge17}. From initial leave-one-out cross-validation studies, it is found that any ridge parameter $\lambda > 0$ improves the rank of the matrix, and for this problem, the fit losses increase with increased $\lambda$. Therefore, $\lambda = 10^{-17}$ is fixed to perform Ridge regression. The optimal direct solver for Ridge regression for this problem is found to be an SVD based solver, whereas OLS regression uses a standard least squares solver. Please refer to a description of the graph theory library in \cref{app:implementation_graphtheorylibrary} of the Appendix for details on the specific direct linear solvers.

Initial studies are also performed using an $R^2$ correlation coefficient for the cross validation metric, while using the weighted loss function for the stepwise regression statistical criteria. This choice of different metrics results in significantly higher errors, of the order $O(10^{-1})$ and the stepwise regression procedure yields a significantly more oscillatory loss curve, as opposed to the desired monotonically increasing curve. The optimal coefficients given the correlation coefficient or the weighted loss metric do not coincide and the cross validation and stepwise regression are competing to the detriment of the overall method. Therefore all cross validation studies will be performed using an identical loss function for all metrics.

Each choice of loss weighting yields similar trends in the loss curves increasing with increased model parsimony, offset by how much the $\ell_{1}$ contributes to the total loss. The high frequency oscillations, as per \cref{fig:DNSobservables} in the phase volume fractions suggests why the fits with solely $\ell_{1}$ have the lowest loss. When looking at the exact ordering of operators that the stepwise regression presents as most important to the model, the different loss weights do not result in identical orderings, even at low numbers of operators where the loss curves are more similar. However, the general trends of which operators are present at various stages of the stepwise regression, denoted by the various plateaus in the loss curves, are very similar across choices of loss weights. We therefore choose to focus on analyzing the fits where solely the $\ell_{1}$ loss is used.

\subsubsection{Reduced order model for phase volume fraction dynamics}\label{sec:physicalsystems_vol_finalmodel}
Looking at the plateaus and sharp increases in the weighted loss curves, particularly for the OLS model in \cref{fig:vol_loss_OLS} over the stepwise iterations, the first plateau extends from $308$ down to approximately past $187$ operators, which corresponds to the rank of the matrix of this basis of operators. In general, as the model becomes more parsimonious, the loss should increase monotonically, however in these cases of lower rank basis, there may be some decreases when certain operators are rejected from the model. The second plateau ends at approximately $47$ operators, which is also the point in the stepwise iterations where the model first has a term that contains a derivative. Practically all iterations between $47$ and $308$ have derivative terms, and this suggests that the purely algebraic polynomial terms in $\{\varphi_{\alpha},l_{\alpha}, N_{\alpha},\vecbar{E}\}$ are over-fitting the model. The derivatives are not seen to play as important a role in the first order dynamics. However, the suggested importance of the algebraic polynomial terms may also be explained by higher order polynomials of the volume averaged quantities being required to represent dynamics controlled by spatially varying fields, as per \cref{eq:local_avg_quantities}.

With loss weights $w_2 = 0,~ w_1 = 1$, the first 10 terms from the stepwise regression for the first order dynamics of $\bar{\varphi}$ are found using OLS and Ridge regression to be
\input{./figures/model_vol_OLS_Ridge17.tex}
It can be seen that both models consist of solely algebraic terms. The OLS and Ridge approaches also generate models with similar terms, with the Ridge regression model having 9 terms compared to 8 terms for OLS regression that contain predominantly strain dependencies.

The first order dynamics, and fits for models with $20$, $50$, and $200$ terms is shown in \cref{fig:vol_fit_OLS_Ridge17}. For the practically full model with $200$ terms, the fits correspond to the data very well at early times, but have difficulty matching all of the small oscillations at late times. As the model becomes more parsimonious to $50$ terms where the sharp increase in the loss occurs, the general trend of the dynamics is captured, however this model misses capturing the full amplitude of the major peaks compared to the more complex models. Finally, with $20$ terms, which corresponds to the number of terms present in the original Cahn-Hilliard equation in $d=2$ dimensions, the large peaks, and small oscillations are not entirely captured, however the general trends are still shown. As discussed, the algebraic polynomial terms are potentially over-fitting, however manage to retain the general behavior of the phase volume fractions changing over time.

\subsubsection{Regression approaches for phase volume fraction dynamics}\label{sec:physicalsystems_vol_regression}
When comparing the performance of the OLS and Ridge regression approaches, both approaches are shown to yield similar trends in the loss curves in \cref{fig:vol_loss_OLS_Ridge17}. However the Ridge regression loss is shown to be greater, by less than an order of magnitude for more complex models than the OLS regression. The resulting fits for the Ridge regression in \cref{fig:vol_fit_OLS_Ridge17} are also shown to capture the general trends of fluctuations quite well, however do not fit the sharper oscillations at later time quite as well as the OLS approach. This increased Ridge loss is reasonable considering the additional constraints, and is balanced by improving the smoothness of the loss curves over the stepwise iterations. The smoother Ridge curves subsequently have less plateaus, and the losses for different loss weights are more distinct. This approach appears to be less affected by changes in the rank of the matrix. Both approaches also behave similarly at the most parsimonious models with less than approximately $47$ terms. These comparable results indicate that Ridge regression is the most appropriate approach for future studies, increasing the robustness of the modelling to more collinear data, and still capturing the essential physical trends. Future cross validation studies, as well as the use of additional constraints on the coefficients in conjunction with the stepwise regression will be useful to determine how to best select the most parsimonious model.

The importance of choosing a physics informed basis is shown clearly by this example. Other choices of basis, such as including other derivatives, or a different set of basis variables like the other phase volume fractions, were shown to have much worse regression losses, and were much less intuitive. Future applications of this graph theoretic method should ensure the forms of the original models in the higher dimensional data are considered when determining how to capture effective physical processes.

\begin{figure}[hpt]
\centering
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=\textwidth]{figures/plot_Loss_vol_OLS.pdf}

	\subcaption{OLS regression.}
	\label{fig:vol_loss_OLS}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=\textwidth]{figures/plot_Loss_vol_Ridge17.pdf}

	\subcaption{Ridge regression with $\lambda = 10^{-17}$.}
	\label{fig:vol_loss_Ridge17}
\end{subfigure}
\caption{Stepwise regression loss curves for the phase volume fractions first order dynamics, labeled by their weighting of the $\ell_{2}$ loss in the weighted residual loss function.}
\label{fig:vol_loss_OLS_Ridge17}
\end{figure}

\begin{figure}[hpt]
\centering
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=\textwidth]{figures/plot_BestFit_vol_200_50_20_OLS.pdf}

	\subcaption{OLS regression.}
	\label{fig:vol_fit_OLS}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=\textwidth]{figures/plot_BestFit_vol_200_50_20_Ridge17.pdf}

	\subcaption{Ridge regression with $\lambda = 10^{-17}$.}
	\label{fig:vol_fit_Ridge17}
\end{subfigure}
\caption{First order dynamics of phase volume fractions for $20$, $50$, and $200$ terms in the forcing function, and use of the $\ell_{1}$ loss in stepwise regression. Backward Euler first time derivative data is shown with the dashed blue curve.}
\label{fig:vol_fit_OLS_Ridge17}
\end{figure}

\subsection{Free energy functional representation} \label{sec:physicalsystems_psi}
To model the free energy of the system, we are guided by the polynomial form of the free energy density in \cref{eq:psi_hom} and write
\begin{align}
	\Psi = \Psi(\bar{\varphi},{\vecbar{E}}) \label{eq:model_psi}
\end{align}
as a function of the single chemical and ${d(d+1)}/{2}$ mechanical variables in $d$ spatial dimensions.

\subsubsection{Proposed model for free energy functional representation}\label{sec:physicalsystems_psi_proposedmodel}

The form of this model suggests our functional representation should be a modified Taylor series at time $j$, about a base value at the time $i$:
\begin{align}
	\Psi(\bar{\varphi}_{j},\vecbar{E}_{j}) = \Psi(\bar{\varphi}_{i},\vecbar{E}_{i})
	~&+~ \gamma^{\bar{\varphi}} \difference[1]{\Psi(\bar{\varphi}_{i},\vecbar{E}_{i})}{\bar{\varphi}} \Delta {{}{\bar{\varphi}}}_{ij} ~+~ \gamma^{\bar{E}_{\alpha \beta}} \difference[1]{\Psi(\bar{\varphi}_{i},\vecbar{E}_{i})}{\bar{E}_{\alpha \beta}} \Delta {{}\bar{E}_{\alpha \beta}}_{ij} \label{eq:model_psitaylor} \\
	~&+~ \gamma^{{\bar{\varphi} \bar{\varphi}}} \unidifference[2]{\Psi(\bar{\varphi}_{i},\vecbar{E}_{i})}{{\bar{\varphi}}} {\Delta {{}{\bar{\varphi}}}_{ij}}^2 \nonumber\\
	 ~&+~ \gamma^{{\bar{E}_{\alpha \beta}}{\bar{E}_{\xi \kappa}}} \difference[2]{\Psi(\bar{\varphi}_{i},\vecbar{E}_{i})}{{\bar{E}_{\alpha \beta}},{\bar{E}_{\xi \kappa}}} \Delta {{}{\bar{E}_{\alpha \beta}}}_{ij} \Delta {{}{\bar{E}_{\xi \kappa}}}_{ij} \nonumber\\
	 ~&+~ \gamma^{{\bar{\varphi}}{\bar{E}_{\alpha \beta}}} \difference[2]{\Psi(\bar{\varphi}_{i},\vecbar{E}_{i})}{{\bar{\varphi}},{\bar{E}_{\alpha \beta}}} \Delta {{}{\bar{\varphi}}}_{ij} \Delta {{}{\bar{E}_{\alpha \beta}}}_{ij} \nonumber \\
	 ~&+~ O(\Delta {{}{\bar{\varphi}}}_{ij}^3) + O(\Delta {{}{\bar{E}_{\alpha \beta}}}_{ij} ^3), \nonumber
\end{align}
where ${\Delta x_{ij}} = x_{j} - x_{i}$ represents the change in the $\varphi_{\alpha}$ or $\vecbar{E}$ state variable between times $j$ and $i$. Please refer \cref{app:implementation_modifiedtaylorseries} of the Appendix for details on how the Taylor series with fit linear coefficients $\gamma$ are formed. Similar to the model for the phase volume fraction evolution, we choose to restrict the dependencies to be $q=4$ variables of a single chemical and three mechanical order parameters in $d=2$ spatial dimensions. We further restrict the Taylor series to be up to $k = \ith[4]$ order, and there are $p = {(q^{k+1}-1)}/{(q-1)} = 341$ possible terms in this model. Here, the expansion is about the initial time $i=0$, and all partial derivatives in the Taylor series are approximated using the non-local calculus definitions in \cref{eq:diff_p}.

\subsubsection{Numerical implementation for free energy functional representation}\label{sec:physicalsystems_psi_numerical}
From initial leave-one-out cross-validation studies, it is found that any ridge parameter $\lambda > 0$ improves the rank of the matrix, and for this problem, decreases the fit losses compared to the OLS method. The fit losses also increase with increased $\lambda$. Therefore, $\lambda = 10^{-17}$ is fixed to perform Ridge regression. Similarly to for the phase volume fractions, the optimal direct solver for Ridge regression for this problem is found to be an SVD based solver, and the OLS regression is performed with a standard least squares solver.

The Ridge regression loss curve in \cref{fig:psi_loss_OLS_Ridge17} is slightly smoother than the OLS regression, and unlike in the phase volume fractions problem, does not increase the loss with $\lambda > 0$. The regularization is shown to be able to increase the rank of the operators, and the robustness of the method, while still producing an accurate model. Here, the free energy is a monotonically decreasing function and is very smooth, so the loss function is the standard $\ell_{2}$ error.

\subsubsection{Reduced order model for free energy functional representation}\label{sec:physicalsystems_psi_finalmodel}
From the loss curves in \cref{fig:psi_loss_OLS_Ridge17}, the sharp increase in the loss occurs at around $70$ terms, which corresponds to the rank of this basis of Taylor series terms. The loss increasing sharply for models with less than this number of terms indicates the importances of these terms in this truncated Taylor series. The loss for the full model is of order $O(10^{-5})$, comparable to other methods \cite{Zhang2020}.

The first 10 terms from the stepwise regression for the $\ith[4]$ order Taylor series functional representation for $\Psi$ are found for OLS and Ridge regression to be:
\input{./figures/model_psi_OLS_Ridge17.tex}

\noindent Here, we have imposed for model consistency, that the coefficient associated with the base Taylor series term is $\gamma_0 = 1$. The most important higher order terms for both models are primarily second and third order derivatives, and both the OLS and Ridge models depend most strongly on the second derivative with respect to phase volume fraction. The OLS model has only one fourth order term, whereas the Ridge model has two fourth order terms, and a first order term. This possibly corresponds to the original free energy density $\psi$ in \cref{eq:psi_hom} being second, third and fourth order powers of the order parameters. However care must be taken on the interpretation of the results of the stepwise regression as being directly induced by the forms of the models for the higher dimensional data. In this example, the stepwise regression indicates a direct correspondence between the form of the reduced order model for $\Psi$, and the form of the spatially dependent $\psi$, however this correspondence should not always be assumed. Physics informed initial choices of basis will ensure the most appropriate reduced order model is found, and later interpretation is specific to each problem.

During the modelling process, it is observed that this Taylor series basis without any regularization is very singular, with rank $70 \ll p$, and condition number on the order of $O(10^{18})$. This can be attributed to the sparsity of the data, as well as potentially the form of the free energy, being a monotonically decreasing function, and derivatives with respect to various quantities being possible numerically linearly dependent. This singular nature of this basis means standard linear regression, computed through an SVD algorithm to solve the normal equations, as per \cref{app:implementation_linearregression} of the Appendix, sets certain singular values to zero, as per a given tolerance. The low rank of this matrix means many singular values are close to zero, and floating point issues related to the machine precision may occur.

Regarding the fits for the free energy, Taylor series with $10$, $30$, and $70$ terms are shown in \cref{fig:psi_fit_OLS_Ridge17}. Models beyond the rank of this basis fit the free energy very well and model a very smooth function. More parsimonious models capture the overall decrease of the free energy, however oscillations in the phase volume fractions and strains are not able to be as smoothed out. Free energies at earlier times also appear to be fit better by all models. When comparing the OLS and Ridge fits, the fits are almost identical for more complex models, however the OLS method has smaller oscillations for the most parsimonious models. This Taylor series approach is shown to be a logical and effective basis of terms to model this smooth free energy function.

\begin{figure}[hpt]
\centering
\includegraphics[width=0.5\textwidth]{figures/plot_Loss_psi_OLS_Ridge17.pdf}
\caption{Stepwise regression loss curve for the $\ith[4]$ order free energy Taylor series functional representation, using OLS and Ridge regression.}
\label{fig:psi_loss_OLS_Ridge17}
\end{figure}

\begin{figure}[hpt]
\centering
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=\textwidth]{figures/plot_BestFit_psi_70_30_10_OLS.pdf}

	\subcaption{OLS regression.}
	\label{fig:psi_fit_OLS}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.49\textwidth}
	\centering

	\includegraphics[width=\textwidth]{figures/plot_BestFit_psi_70_30_10_Ridge17.pdf}

	\subcaption{Ridge regression with $\lambda = 10^{-17}$.}
	\label{fig:psi_fit_Ridge17}
\end{subfigure}
\caption{Fitted curves for the $\ith[4]$ order free energy Taylor series functional representation, with $10$, $30$, and $70$ terms. DNS free energy data is shown with the dashed blue curve, and is essentially coincident with the higher order models.}
\label{fig:psi_fit_OLS_Ridge17}
\end{figure}

\section{Conclusions}
\label{sec:concl}
In this work, we develop a graph theoretic approach for reduced order modelling of physical systems. After extending a non-local calculus as the foundation of this approach, we observe there are challenges and compromises when imposing consistency with differential calculus definitions of partial derivatives. Constraints placed on the graph's edge weights lead to numerical issues, but also interpretation issues if the non-local derivatives have error constants that are consistently large. A more in depth comparison of the behavior of non-local and differential derivatives is necessary in future works.

Applying this calculus to studying volume averaged quantities from high dimensional direct numerical simulations by representing these quantities on a graph shows the intuitive, general, and effective nature of this approach. With an example physical system of the phase evolution of a mechano-chemical multi-crystalline solid, we see that using a physics informed basis of operators can model both differential equations, and functional representations. Highly oscillatory or smooth monotonically evolving functions can be fit using linear regression and an adequate number of terms, with errors on the order of $O(10^{-4})$. As a comparison, recent studies of the same microstructure system using neural networks \cite{Zhang2020} achieve errors on the order of $O(10^{-6})$.

We also note that near the beginning of the stepwise regression when the models contain the full set of operators, there is sometimes, for example in \cref{fig:vol_loss_OLS_Ridge17}, decrease in the loss function as the model becomes initially more parsimonious. This suggests that although we are choosing an intuitive, physics informed basis, there may be some operators that worsen the model, and must be rejected empirically. A more rigorous study to more precisely determine which operators from the appropriate basis should be included in the initial model will be conducted in future communications.

\textbf{Remark}: There remains an open question of how accurate these graph theoretic methods can be. This question may be related to the non-local nature of the derivatives, and whether they are appropriate to model strictly locally behaving systems. Whether or not these methods can achieve the error tolerances of other direct numerical simulations, or methods such as neural networks, will show the true effectiveness of this approach. An error analysis in \cref{app:error,app:error_p} also indicates that the error in a one dimensional linear Taylor series model only scales linearly with the sparsity of the data, compared to quadratically for standard interpolation methods, and further only scales sub-linearly for higher dimensional systems. This sub-optimal error scaling is directly attributed to the non-local calculus definitions and the form of the edge weights, and suggests modifications to these definitions may improve the scaling.

There is also the compromise of the robustness and consistency of this modelling approach versus its accuracy. As discussed in \cref{app:implementation_linearregression} of the Appendix, the improvement of the multicollinearity of the basis of operators, and ensuring a consistent order of operators yielded by the stepwise regression, from regularization approaches may come at the cost of decreased loss. The physical systems studied indicate a physically appropriate basis of operators may still be numerically collinear, and regularization aids in the robustness and consistency of the modelling, without significant loss in accuracy.

We note, however, that in any reduced-order model, the goal is a parsimony of representation, which is achieved here by $O(10)$ components in the state vector, at the cost of some loss of fidelity. This approach appears to require $O(10^2)$ terms in order to get accurate fits of the data, however these higher order models involve higher order derivatives that may suffer numerical errors due to their large values. This compromise of model accuracy versus model complexity is another open question to be studied. To justify the requirements of these higher order models, we note that since the original PDE itself in the Cahn-Hilliard equation has nearly 20 terms when fully expanded, this number could be regarded as a lower bound for the complexity of the reduced-order models. Then, with the significant saving in dimensionality of the state vector over that of a DNS, even $O(10^2)$ operators in the reduced-order model is not unreasonable. This is especially reasonable when this representation is required to reproduce the dynamics of the vastly higher dimensional DNS.

Given these findings, this graph theoretic approach is shown to be very general, and many applications to a wide range of physical system, and system scales are possible. There are also several interesting studies that can be conducted in future works regarding the methods themselves. First it will be important to develop edge weight constraints for finite, discrete graphs, as opposed to current use of constraints for finite continuous graphs. Rigorous approximations for these discrete sums over the volume of the graph will possibly constrain the numerical issues brought on by the consistency demands, and will allow the local and non-local calculus to be better compared. The forms of the decaying edge weights can also be generalized to include several decay parameters that reflect independent behaviors of the state vector components. Finally, aspects of the stepwise regression, including the effects of different loss functions, and how the independence of operators in the basis plays a role in how and which operators are rejected from the model, will aid in the fitting process. \\

\section*{Acknowledgments}
\noindent MD and KG gratefully acknowledge the support of the National Science Foundation, United States grant \textbf{\#1729166}. "DMREF/GOALI: Integrated Framework for Design of Alloy-Oxide Structures".

\pagestyle{bibliography}