tcb@breakable

# Measuring dependencies between variables of a dynamical system using fuzzy
affiliations

Niklas Wulkow  Department of Mathematics and Computer Science, Freie
Universität Berlin, Germany  Zuse Institute Berlin, Germany

###### Contents

  1. 1 Introduction
  2. 2 Scalable Probabilistic Approximation (SPA)
    1. 2.1 Transformation of coordinates
    2. 2.2 Estimation of dynamics
    3. 2.3 Model estimation between two processes
  3. 3 Dependency Measures
    1. 3.1 Schatten-1 norm
    2. 3.2 Average row variance
    3. 3.3 Minimizers and maximizers of the dependency measures
    4. 3.4 Differences to related practices
  4. 4 Numerical examples
    1. 4.1 Coupled two-species logistic map
    2. 4.2 Continuous movement by stochastic diffusions
    3. 4.3 Multidimensional autoregressive processes
  5. 5 Real-World Example 1
  6. 6 Real-world example: basketball player movement
    1. 6.1 Applying the dependency analysis to the basketball data
  7. 7 Conclusion
  8. A Appendix

##  1 Introduction

  * •

Problem formulation

    * –

Detect driving variables in a dynamical system from time series data.

    * –

Where is this important?

    * –

There exist: Granger causality [granger, kirchgaessner]: If $X$ forces $Y$,
then values of $X$ should help predict $Y$. Used in economics since 1960s.
Usually linear AR models [stokes]. Can be extended, of course. But especially
for stochastic processes or chaotic systems prediction error are not always
meaningful.

    * –

Convergent Cross Mapping [sugihara]. But it requires dynamics to stay on an
attracting manifold and suffers in the presence of noise. Plus, if two
variables of a dynamical system are too strongly connected, it fails to
indicate a uni-directional influence [sugiharaCCM2015].

    * –

Furthermore,… lots of literature

  * •

Our contribution

    * –

Almost-non-parametric method (in contrast to Granger causality. Build a simple
linear model to predict one variable from the other with direct
interpretability of system matrix.

    * –

We define two measures that quantify the strength of influence.

    * –

We use the fuzzy clustering method SPA [spa].

    * –

We introduce the two measures and prove that they are in line with a desired
interpretation of the system matrix. We then apply it to multiple examples,
including player movements from a Basketball game.

  * •

In Section…

##  2 Scalable Probabilistic Approximation (SPA)

Scalable Probabilistic Approximation (SPA) [spa] is a versatile method that
transforms points from a $D$-dimensional state space into an alternative
representation using new $K$-dimensional coordinates. If $K$ is at most than
$D$, SPA serves as a dimension reduction method by representing points as
closely as possible. If $K>D$, SPA can be seen as a fuzzy clustering method
which assign points to landmark points in $D$-dimensional space. For points
coming from a dynamical system $X_{t}=F(X_{t-1})\in\mathbb{R}^{D}$, SPA then
finds a linear model which propagates states in these $K$-dimensional
coordinates. The first step, the transformation to other coordinates, will be
called SPA I while the reconstruction of the dynamics in these coordinates
will be called SPA II. We assume dynamical systems of the form

###  2.1 Transformation of coordinates

The mathematical formulation of SPA I is: let
$\textbf{X}=[X_{1},\dots,X_{T}]\in\mathbb{R}^{D\times T}$. Then solve

$\begin{split}[\Sigma,\Gamma]&=\operatorname*{arg\,min}\|\textbf{X}-\Sigma\Gamma\|_{F}\\\
\text{subject to }&\Sigma=[\sigma_{1},\dots,\sigma_{K}]\in\mathbb{R}^{D\times
K},\quad\Gamma=[\gamma_{1},\dots,\gamma_{T}]\in\mathbb{R}^{K\times T},\\\
&(\gamma_{t})_{i}\geq 0,\quad\sum_{i=1}^{K}(\gamma_{t})_{i}=1.\end{split}$
(SPA I)

Is was discussed in [wulkowmSPA] that for $K\leq D$, the representation of
points in this way is the orthogonal projection onto a convex polytope with
vertices given by the columns of $\Sigma$. The coordinates $\gamma$ then
specify the position of this projection with respect to the vertices of the
polytope and are called barycentric coordinates (BCs). A high entry in such a
coordinate then signals closeness of the projected point to the vertex.

######  Remark 2.1.

This representation of points has already been introduced in PCCA+ [weber].

For $K>D$, however, the interpretation is different. In [spa], the
interpretation of a probabilistic clustering is introduced. According to the
authors, the entries of a $K$-dimensional coordinate of a point denote the
probabilities to be inside a certain box around a landmark point, given by a
column in $\Sigma$. One can generalize this interpretation to affiliations to
these landmark points, again in the sense of closeness. A BC $\gamma$ then
denotes the distribution of affiliations to each landmark point. Figure 1
shows the representation of a point in $\mathbb{R}^{2}$ with respect to four
landmark points.

Figure 1: Representation of a point $X$ by barycentric coordinates
$\gamma_{1},\dots,\gamma_{4}$ with respect to the vertices of a polytope
$\sigma_{1},\dots,\sigma_{4}$.

For $K>D+1$, the representation with barycentric coordinates is not unique.
Therefore, we define the representation of a point $X$ analogously to
[wulkowmSPA] in the following way:

$\begin{split}\rho_{\Sigma}(X,\gamma)&:=\operatorname*{arg\,min}\limits_{\gamma^{*}}\|\gamma-\gamma^{*}\|_{2}\\\
\text{ s.t.
}\gamma^{*}&=\operatorname*{arg\,min}\limits_{\gamma^{\prime}}\|X-\Sigma\gamma^{\prime}\|_{2}\text{
with }\gamma^{\prime}_{\bullet},\gamma^{*}_{\bullet}\geq 0\text{ and
}\|\gamma^{\prime}\|_{1},\|\gamma^{*}\|_{1}=1.\end{split}$ (2.1)

$\gamma$ should be selected among all barycentric coordinates that represent
$X$ without loss so that it is closest to the reference coordinate
$\gamma^{\prime}$. In the case of a dynamical system $X_{t}=F(X_{t-1})$ we
select $\gamma_{t-1}$ as reference coordinate for the point $X_{t}$. This
gives

$\gamma_{t}=\rho_{\Sigma}(X_{t},\gamma_{t-1})=\rho_{\Sigma}(F(X_{t-1}),\gamma_{t-1})=\rho_{\Sigma}(F(\Sigma\gamma_{t-1}),\gamma_{t-1})=:\textbf{v}(\gamma_{t-1}).$
(2.2)

With this, $\gamma_{t}$ solely depends on $\gamma_{t-1}$. We therefore have
formulated a time-discrete dynamical system in the barycentric coordinates. By
choosing the reference coordinate as $\gamma_{t-1}$ we assert that the steps
taken are as small as possible.

######  Remark 2.2.

Note that the solution for SPA I is never unique [spa] if $K>1$. Landmark
points can be placed differently. Especially if $K>D+1$, there are vastly
different sets of landmark points which all allow to solve SPA I exactly. In
order to solve SPA I we iteratively minimize the objective function by
separately solving for $\Sigma$ and $\Gamma$. For this, we randomly draw
initial values. Therefore, for $K>D+1$, an exact solution can be reached by
placing $D+1$ points so that all data points lie inside their convex hull (all
convex combinations of them) while the other landmark points can be chosen
arbitrarily. As a consequence, the placement of landmark points depends
strongly on the randomly chosen initial values of the optimization process.

###  2.2 Estimation of dynamics

For the estimation of the dynamics, we use SPA II which estimates a linear
approximation of v that by construction maps BCs to BCs.

$\begin{split}\Lambda&=\operatorname*{arg\,min}_{\Lambda^{*}\in\mathbb{R}^{K\times
K}}\|[\gamma_{2}|\cdots|\gamma_{T}]-\Lambda^{*}[\gamma_{1}|\cdots|\gamma_{T-1}]\|_{F},\\\
\text{subject to }&\Lambda\geq 0\text{ and
}\sum_{k=1}^{K}\Lambda_{k,\bullet}=1.\end{split}$ (SPA II)

$\Lambda$ is a column-stochastic matrix as typically used in Markov State
Models [sarich, husic] for stochastic processes. With $\Lambda$, we can
construct dynamics on the level of the BCs with the possibility to transform
back to the original state space by multiplication with $\Sigma$, since

$X_{t}=\Sigma\gamma_{t}=\Sigma\textbf{v}(\gamma_{t-1})\approx\Sigma\Lambda\gamma_{t-1}.$
(2.3)

Note that such a linear model is generally not suited to reconstruct the long-
term behaviour of nonlinear dynamics. In [wulkowmSPA] an extension to SPA for
this purpose is presented.

######  Remark 2.3.

In [spa] is shown a way to combine SPA I and SPA II in a single SPA I problem.
The landmark points are then selected so that the training error of the SPA II
problem can be set to $0$.

###  2.3 Model estimation between two processes

Given two processes, $X_{t}=F(X_{t-1})\in\mathbb{R}^{D}$ and
$Y_{t}=G(Y_{t-1})\in\mathbb{R}^{E}$, we can perform SPA I for both finding
landmark points as the columns of $\Sigma^{X}\in\mathbb{R}^{D\times
K_{X}},\Sigma^{Y}\in\mathbb{R}^{E\times K_{Y}}$ and BCs $\gamma^{X}_{t}$ and
$\gamma^{Y}_{t}$ for $t=1,\dots,T$. Let us find a column-stochastic matrix
$\Lambda_{XY}$ that fulfils

$\begin{split}\Lambda_{XY}=\operatorname*{arg\,min}_{\Lambda^{*}\in\mathbb{R}^{K_{Y}\times
K_{X}}}\|[\gamma^{Y}_{2}|\cdots|\gamma^{Y}_{T}]-\Lambda^{*}[\gamma^{X}_{1}|\cdots|\gamma^{X}_{T-1}]\|_{F}.\end{split}$
(2.4)

Rather than approximating a dynamical system, $\Lambda_{XY}$ represents a
model from $X_{t-1}$ to $Y_{t}$ on the level of the BCs, i.e., tries to
predict subsequent values of $Y$ using only $X$.

Now, let us assume that $X$ in fact has direct influence on $Y$, meaning that
there exists a function $\textbf{w}(Y_{t-1},X_{t-1})=Y_{t}$. Then similarly as
when constructing the dynamical system in the BCs previously in Eq. (2.2), we
can observe

$\gamma^{Y}_{t}=\rho_{\Sigma^{Y}}(Y_{t},\gamma^{Y}_{t-1})=\rho_{\Sigma^{Y}}(\textbf{w}(Y_{t-1},X_{t-1}),\gamma^{Y}_{t-1})=\rho_{\Sigma^{Y}}(\textbf{w}(\Sigma^{X}\gamma^{X}_{t-1},\Sigma^{Y}\gamma^{Y}_{t-1}),\gamma^{Y}_{t-1})=:\textbf{u}(\gamma^{X}_{t-1},\gamma^{Y}_{t-1}).$
(2.5)

$\gamma^{Y}_{t}$ therefore directly depends on $\gamma^{X}_{t-1}$ and
$\gamma^{Y}_{t-1}$ while $\Lambda_{XY}$ attempts to predict $\gamma_{t}^{Y}$
using only $\gamma_{t-1}^{X}$. One can therefore interpret $\Lambda_{XY}$ as
an approximation to the conditional expectation of $\gamma^{Y}_{t}$ given
$\gamma^{X}_{t-1}$, assuming that $\gamma_{t-1}^{Y}$ is distributed by a
function $\mu_{Y}$, i.e.,

$\Lambda_{XY}\gamma\approx\mathbb{E}_{\mu_{Y}}[\gamma^{Y}_{t}|\gamma^{X}_{t-1}=\gamma].$
(2.6)

Assuming Eq. (2.6) is a close approximation, we can assert

$\mathbb{E}_{\mu_{Y}}[\gamma^{Y}_{t}|\gamma^{X}_{t-1}]\approx\Lambda_{XY}\gamma^{X}_{t-1}\Leftrightarrow\mathbb{E}_{\mu_{Y}}[\gamma^{Y}_{t}|\gamma^{X}_{t-1}]\approx\sum_{j=1}^{K_{X}}(\Lambda_{XY})_{|j}(\gamma^{X}_{t-1})_{j}$
(2.7)

where $(\Lambda_{XY})_{|j}$ is the $j$th column of $\Lambda_{XY}$. A
prediction for $\gamma^{Y}_{t}$ is therefore constructed using a weighted
average of the columns of $\Lambda_{XY}$. The weights are the entries of
$\gamma^{X}_{t-1}$.

######  Remark 2.4.

Note that the same argumentation starting in Eq. (2.4) holds if we choose a
time shift of length $\tau>0$ and consider information of $X_{t-\tau}$ about
$Y_{t}$. If $\tau=1$ we will simply write $\Lambda_{XY}$ but generally write
$\Lambda_{XY}^{(\tau)}$.

##  3 Dependency Measures

In the following, we will use our interpretation of $\Lambda_{XY}$ introduced
above to define two methods that quantify the strength of dependence of
$Y_{t+1}$ on $X_{t}$. The intuition can be illustrated as follows: if a
variable $X$ carries zero information on the future state of another variable,
$Y$, this means, that weights $\gamma^{X}_{t-1}$ in Eq. (2.7) should have no
influence on the composition of $\gamma^{Y}_{t}$ through the weighted average
of columns of $\Lambda_{XY}$. In turn, this means that the columns of
$\Lambda_{XY}$ should be equal to each other. Since in this case
$\Lambda_{XY}\gamma^{X}\equiv(\Lambda_{XY})_{|1}=\dots(\Lambda_{XY})_{|K_{X}}$,
in SPA II we then essentially strive to minimize

$\lambda=\operatorname*{arg\,min}_{\lambda^{*}\in\mathbb{R}^{K_{Y}}}\|[\gamma^{Y}_{2}|\cdots|\gamma^{Y}_{T}]-[\lambda^{*},\dots,\lambda^{*}]\|_{F}$
(3.1)

From classical theory on statistics [spokoiny], the solution of this problem
is given by the mean of the time series $\gamma^{Y}_{2},\dots,\gamma^{Y}_{T}$.
This makes sense noting that if $\gamma^{Y}_{t}$ is independent of
$\gamma^{X}_{t-1}$, this means that

$\mathbb{E}_{\mu_{Y}}[\gamma^{Y}_{t+}|\gamma^{X}_{t-1}]=\mathbb{E}_{\mu_{Y}}[\gamma^{Y}_{t}].$
(3.2)

Therefore, each column of $\Lambda_{XY}$ should be an approximation of
$\mathbb{E}_{\mu_{Y}}[\gamma^{Y}_{t}]$.

In the opposite case, if $X$ has strong influence on $Y$, the columns of
$\Lambda_{XY}$ should deviate from each other strongly so that $\Lambda_{XY}$
has full rank. This is the main idea behind the shortly presented methods.

Clearly, even if $X$ provided no information on $Y$, with limited data it
seems unlikely that all columns of $\Lambda_{XY}$ are exactly equal to each
other so that generally it will be of full rank. We will therefore quantify
differences between its columns to measure the strength of dependence of one
variable on the other.

Note that we do not attempt to quantify causation. Rather, we aim to quantify
how much information one variable carries for the prediction of another
variable which one can view as an indication for causal influence.

###  3.1 Schatten-1 norm

For the first measure, we consider the Singular Value Decomposition (SVD)
[golub] of a matrix $\Lambda$, given by
$\Lambda=USV^{T}\in\mathbb{R}^{K_{Y}\times K_{X}}$.
$S\in\mathbb{R}^{K_{Y}\times K_{X}}$ is a matrix which is only non-zero in the
entries $(i,i)$ for $i=1,\dots,min(K_{X},K_{Y})$ which are given by
$s_{1},\dots,s_{min(K_{X},K_{Y})})\geq 0$. $U\in\mathbb{R}^{K_{Y}\times
K_{Y}}$ and $V\in\mathbb{R}^{K_{X}\times K_{X}}$ fulfil some orthogonality
properties and consist of columns $u_{i},v_{i}$. We can thus write $\Lambda$
as

$\Lambda=\sum\limits_{i=1}^{r}u_{i}v_{i}^{T}s_{i}$

A classic linear algebra result asserts that $rank(\Lambda)=\\#\\{s_{i}>0\\}$.
As a consequence, if some of the $s_{i}$ are close to $0$, then this means
that only a small perturbation is sufficient to make $\Lambda$ a matrix of
lower rank. Therefore, we use the sum of singular values as a continuous
measure of the rank and thus of the difference in the rows of $\Lambda$.

######  Definition 3.1 (Schatten-1 norm).

Let the SVD of a matrix $\Lambda\in\mathbb{R}^{K_{Y}\times K_{X}}$ be given by
$\Lambda=USV^{T}$ with singular values $s_{1},\dots,s_{min(K_{X},K_{Y})}$.
Then the Schatten-1 norm of $\Lambda$ is defined as

$\|\Lambda\|_{1}:=\sum\limits_{i=1}^{min(K_{X},K_{Y})}s_{i}.$ (3.3)

###  3.2 Average row variance

As our second dependency measure, we directly quantify the difference of
columns of a matrix $\Lambda$ using the mean statistical variance per row. We
therefore consider every row and compute the variance between its entries,
thereby comparing the columns with respect to this particular row. We then
take the mean of these variances across all rows.

######  Definition 3.2 (Average row variance).

For a matrix $\Lambda\in\mathbb{R}^{K_{Y}\times K_{X}}$, let
$\bar{\Lambda}_{i|}$ denote the mean of the $i$th row of $\Lambda$. Let

$\nu_{i}:=\frac{1}{K_{X}}\sum\limits_{j=1}^{K_{X}}(\Lambda_{ij}-\bar{\Lambda}_{i|})^{2}$

be the variance of the entries of the $i$th row. Then we define as the average
row variance

$\nu(\Lambda):=\frac{1}{K_{Y}}\sum\limits_{i=1}^{K_{Y}}\nu_{i}.$ (3.4)

###  3.3 Minimizers and maximizers of the dependency measures

About $\|\cdot\|_{1}$ and $\nu$, we can proof properties that validate why
they represent sensible measures for the strengths of dependency between two
processes. For this we need the following definition.

######  Definition 3.3 (Permutation matrix).

As a permutation matrix we define a matrix $A\in\\{0,1\\}^{n\times n}$ such
that every row and column contains exactly one $1$.

Then we obtain the following results, whose proofs can be found in the
Appendix. For the Schatten-1-norm, we find for the maximizers and minimisers:
NW: Rot bedeutet ”noch nicht bewiesen”, grün bedeutet ”bewiesen” in Appendix

######  Proposition 3.1 (Maximal Schatten-1 norm, $K_{Y}\geq K_{X}$).

Let $\Lambda\in\mathbb{R}^{K_{Y}\times K_{X}}$ with $K_{Y}\geq K_{X}$. Then
the Schatten-1 norm of $\Lambda$ obtains the maximal value $K_{X}$ if and only
if deletion of $K_{Y}-K_{X}$ rows of $\Lambda$ yields a $K_{X}\times K_{X}$
permutation matrix.

######  Proposition 3.2 (Maximal Schatten-1 norm, $K_{Y}<K_{X}$).

Let $\Lambda\in\mathbb{R}^{K_{Y}\times K_{X}}$ with $K_{Y}<K_{X}$. Then the
Schatten-1 norm of $\Lambda$ is maximized if and only if $\Lambda$ contains an
$K_{Y}\times K_{Y}$ permutation matrix and the matrix of the remaining
$K_{X}-K_{Y}$ columns can be extended by $K_{Y}$ columns to a permutation
matrix.

######  Proposition 3.3 (Minimal Schatten-1 norm).

The Schatten-1 norm of a column stochastic $K_{X}\times K_{Y}$-matrix $A$ is
minimal if and only if $A_{ij}\equiv\frac{1}{n}$ and its minimal value is
equal to $1$.

For the average row variance, we can derive the following results:

######  Proposition 3.4 (Maximal average row variance, $K_{Y}\geq K_{X}$).

Let $\Lambda\in\mathbb{R}^{K_{Y}\times K_{X}}$ with $K_{Y}\geq K_{X}$. Then
the average row variance of a column-stochastic $K_{Y}\times K_{X}$-matrix
$\Lambda$ is maximal if and only if deletion of $K_{Y}-K_{X}$ rows of
$\Lambda$ yields a $K_{X}\times K_{X}$ permutation matrix.

######  Proposition 3.5 (Maximal average row variance, $K_{Y}<K_{X}$).

Let $\Lambda\in\mathbb{R}^{K_{Y}\times K_{X}}$ with $K_{Y}<K_{X}$. Then the
average row variance of a column-stochastic $K_{Y}\times K_{X}$-matrix
$\Lambda$ is maximal if and only if $\Lambda$ contains an $K_{Y}\times K_{Y}$
permutation matrix and the matrix of the remaining $K_{X}-K_{Y}$ columns can
be extended by $K_{Y}$ columns to a permutation matrix.

######  Proposition 3.6 (Minimal average row variance).

The average row variance of a column-stochastic $K_{Y}\times K_{X}$-matrix
$\Lambda$ obtains the minimal value $0$ if and only if all columns are equal
to each other.

In summary, the maximizing and minimizing matrices of $\|\cdot\|$ and $\nu$
are identical and are of the following forms:

$\text{Maximal: }\begin{pmatrix}0&1&0\\\ 0&0&1\\\ 1&0&0\\\
0&0&0\end{pmatrix},\begin{pmatrix}0&1&0&1\\\ 0&0&1&0\\\ 1&0&0&0\\\
\end{pmatrix},\quad\text{Minimal:}\begin{pmatrix}\frac{1}{n}&\dots&\frac{1}{n}\\\
\vdots&&\vdots\\\ \frac{1}{n}&\dots&\frac{1}{n}\end{pmatrix}.$ (3.5)

These results show that the two dependency measures $\|\cdot\|_{1}$ and $\nu$
fulfil important intuitions: they are minimal, when information about
$X_{t-1}$ gives us no information about $Y_{t}$ because in this case all
columns of $\Lambda$ should be identical. Maximal dependence is detected if
the information about $X_{t-1}$ yields maximally strong information about
$Y_{t}$. This happens when $\Lambda$ is, respectively can be reduced or
extended to, a permutation matrix. This is also illustrated in Figure 2.

Figure 2: Illustration of the intuition behind the dependency measures. If the
distribution of $Y$ is independent of $X$, the matrix $\Lambda_{XY}$ will have
similar columns (left). If the distribution of $Y$ is strongly dependent of
$X$, the columns of $Y$ will be very different from each other.

We will store the calculated values for $\|\cdot\|_{1}$ and $\nu$ in tables,
resp. matrices of the form

$M_{\|\cdot\|_{1}}=\begin{pmatrix}\|\Lambda_{XX}\|_{1}&\|\Lambda_{XY}\|_{1}\\\
\|\Lambda_{YX}\|_{1}&\|\Lambda_{YY}\|_{1}\end{pmatrix},M_{v}=\begin{pmatrix}\nu(\Lambda_{XX})&\nu(\Lambda_{XY})\\\
\nu(\Lambda_{YX})&\nu(\Lambda_{YY})\end{pmatrix}.$ (3.6)

Then for each of these matrices, the property $M-M^{T}$ should be interesting
for us, because they matrices contain the differences between dependency
measures, stating how strongly $X$ depends on $Y$ compared to $Y$ depending on
$X$. We therefore define

$\delta(M)_{ij}=\frac{M_{ij}-M_{ji}}{max(M_{ij},M_{ji})}$ (3.7)

as the relative difference between dependencies in both directions.

###  3.4 Differences to related practices

#### Simple linear correlations

The presented measures might seem strongly related to the computation of the
linear correlation

$C(X,Y)=\frac{1}{T-1}\sum_{t=1}^{T-1}(X_{t}-\bar{X})(Y_{t+1}-\bar{Y})^{T}$
(3.8)

(for $\tau=1$) where $\bar{X},\bar{Y}$ are the component-wise averages.
However, $C$ can only detect global linear patterns between $X$ and $Y$. In
contrast, we transform points into a higher-dimensional space by expressing
them by barycentric coordinates with $K>D$. While we still determine a linear
operator between the variables, given by $\Lambda_{XY}^{(\tau)}$, itself a
local approximation of a potentially nonlinear function, denoted earlier by w.
Furthermore, upon perturbations to the function w, $\Lambda_{XY}^{(\tau)}$
should react in a nonlinear way by construction of the SPA II problem. The
dependency measures then are nonlinear functions on $\Lambda_{XY}^{(\tau)}$.
In the examples that will follow, using linear correlations could generally
not uncover unidirectional influences while our measures were able to do so.

#### Granger causality

A prominent method to measure the influence of one variable on another is by
employing the Granger causality framework [granger, kirchgaessner]. It
functions by determining two models of the form

$\begin{split}Y_{t+1}&=f(X_{t},\dots,X_{t-p+1},Y_{t},\dots,Y_{t-q+1})\\\
Y_{t+1}&=g(Y_{t},\dots,Y_{t-q+1})\end{split}$ (3.9)

from training data and using them to compute subsequent values of $Y_{t}$ on
testing data which was not used for training. The prediction errors of $f$ and
$g$ are then compared. If $f$, which uses information of $X$, gives a
significantly better prediction error, then it is deduced that $X$ influences
$Y$.

Typical model forms for $f$ and $g$ are linear autoregressive models
[brockwell] which are described in more detail in the next section. It is
pointed out in [sugihara] that using past terms of $X$ and $Y$ can constraint
the interpretability of the result, since if $Y$ forces $X$ information about
$Y$ is stored in past terms of $X$ due to the delay-embedding theorem of
Takens [takens] (please see [sugihara] including its supplement for details).
Then if $Y$ can be predicted from past terms of $X$ it actually is a sign that
$Y$ forces $X$, not vice versa. In [sugihara] examples are shown where the
Granger causality method fails to detect influences between variables. This
makes the interpretation of the Granger results more difficult.

This effect does not occur when dispensing of the past terms and instead
fitting models $Y_{t+1}=f(X_{t},Y_{t})$ and $Y_{t+1}=g(Y_{t})$. However, in
systems which are highly stochastic or chaotic, meaning that from similar
initial conditions diverging trajectories emerge, even an accurate model can
be prone to give weak prediction errors. In such cases the prediction error
often times has limited meaning.

Furthermore, even if $X$ influences $Y$, one has to select a suitable model
family for $f$ and $g$ so that this actually shows. The selection of the model
family can be a challenging task of its own.

Nevertheless, Granger causality can be a strong tool for the detection of
influences, e.g., as shown for a Gene Regulatory Network in [sincerities].

#### Discretization by boxes instead of landmark points

Earlier the similarity between the model constructed solving SPA II and Markov
State Models (MSMs) was mentioned. In MSMs, one discretizes the state space
into boxes and statistically estimates the transition probabilities of the
state of a dynamical system between the boxes, typically by considering the
relative frequencies of transitions. One then obtains a column-stochastic
transition matrix that contains these relative frequencies. In the same
manner, one could compute this matrix for the frequencies that a variable $Y$
is in a certain box at time $t+\tau$ given that a variable $X$ is in a certain
box at time $t$ and apply the dependency measures to this transition matrix to
assess how meaningful the information about a variable $X$ is for the future
value of $Y$. However, MSMs only describe the dynamics with a precision that
depends on the fineness of their discretization. However, fixing the edge
length of each box, the number of boxes increases exponentially with the
dimension of points. One then requires very long time series for a robust
estimation of the transition probabilities. The advantage in the SPA I
representation of points is that we can derive a column-stochastic matrix but
can maintain a lossless representation of points with $K\sim D$ while $K>D$ is
the only prerequisite.

##  4 Numerical examples

We now apply the two measures to different examples of dynamical systems with
varying complexity. In order to assess their efficacy, we explicitly install
unidirectional dependencies in the formulation of the dynamics and
investigate, if these are detected.

###  4.1 Coupled two-species logistic map

The first example is a coupled two-species logistic difference system as used
in [sugihara]. It describes oscillating behaviour of two quantities while each
is forced by the other, although with different strengths. The behaviour of
each variable is mainly determined by the logistic map [may] which for one
variable reads

$X_{t+1}=rX_{t}(1-X_{t}).$ (4.1)

There is extensive literature on the dependence of the behaviour of the system
on the parameter $r$, e.g., [tsuchuya, nagashima, may]. For $r<3$, the system
will converge to the value $\frac{r-1}{r}$ while for $r>4$ the system diverges
for most initial values. The most interesting behaviour can be observed for
for $r\in(3,4)$. Then the system converges to a set of points which it
oscillates between and whose number increases with $r$ approaching $4$ (see
the bifurcation diagram in Figure 3, right).

Figure 3: Left: Trajectory of the coupled two-species logistic system (4.2).
$X_{t}$ in blue, $Y_{t}$ in orange. Right: Bifurcation diagram of the one-
species logistic map (4.1) (taken from the Wikipedia entry on the logistic map
[lmapwiki]).

We define a coupled two-species logistic map, taken from [sugihara], as

$\begin{split}X_{t+1}&=3.8X_{t}(1-X_{t})-0.02Y_{t}X_{t}\\\
Y_{t+1}&=3.5Y_{t}(1-Y_{t})-0.1X_{t}Y_{t}\end{split}$ (4.2)

with $X_{0}=Y_{0}=0.8$. The trajectory is depicted in Figure 3 ,left.

#### Dependency analysis

Through the prefactor $0.1$ of $X_{t}$ in the equation for $Y_{t+1}$ compared
to the lower prefactor $0.02$ of $Y_{t}$ in the equation for $X_{t+1}$, $X$
has a direct influence on $Y$ that is higher than the direct influence that
flows in the other direction. We will uncover this using SPA and the
dependency measures on the resulting transition matrices.
For the analysis, we use $K_{X}=K_{Y}=10$ landmark points for both processes
on realisations of length 1800 to produce representations of points in
barycentric coordinates in $\mathbb{R}^{10}$. We then perform the analysis
explained in Section 3 on the first 200,400,600… up to all 1800 data points
with a time shift of one time step. We see in Figure 4 that the two dependency
measures detect the different strengths of influences.

Figure 4: Dependency measures $\|\cdot\|_{1}$ and $\nu$ for the logistic map
system Eq. (4.2) for increasing time. The dependency measures of $Y$ on $X$ is
coloured in blue and the measures of $X$ on $Y$ in red. We see that both
measures uncover the higher influence of $X$ to $Y$ that is anchored in the
system.

For the trajectory of length $1800$, we obtain the result

$M_{\|\cdot\|_{1}}:\begin{tabular}[]{|l|l|l|l|l|l|}\hline\cr From $\downarrow$
to $\rightarrow$&$X_{t+1}$&$Y_{t+1}$\\\ \hline\cr$X_{t}$&5.76&2.92\\\
\hline\cr$Y_{t}$&2.45&6.09\\\
\hline\cr\end{tabular}\Rightarrow\delta(M_{\|\cdot\|_{1}})=\begin{pmatrix}\hbox{\pagecolor[rgb]{1.00,1.00,0.5}0}&\hbox{\pagecolor[rgb]{0.00,1.00,0.5}0.16}\\\
\
\hbox{\pagecolor[rgb]{1.00,0.00,0.5}-0.16}&\hbox{\pagecolor[rgb]{1.00,1.00,0.5}0}\\\
\end{pmatrix}\\\ $ $M_{\nu}:\begin{tabular}[]{|l|l|l|l|l|l|}\hline\cr From
$\downarrow$ to $\rightarrow$&$X_{t+1}$&$Y_{t+1}$\\\
\hline\cr$X_{t}$&0.0521&0.0131\\\ \hline\cr$Y_{t}$&0.0053&0.0542\\\
\hline\cr\end{tabular}\Rightarrow\delta(M_{\nu})\begin{pmatrix}\hbox{\pagecolor[rgb]{1.00,1.00,0.5}0}&\hbox{\pagecolor[rgb]{-0.00,1.00,0.5}0.6}\\\
\
\hbox{\pagecolor[rgb]{1.00,-0.00,0.5}-0.6}&\hbox{\pagecolor[rgb]{1.00,1.00,0.5}0}\\\
\end{pmatrix}$

While the Schatten-1 norm produces as relative difference of $0.16$, the
average row variance emphasises the stronger dependence of $Y$ on $X$ by a
relative difference in dependencies of $0.6$.

###  4.2 Continuous movement by stochastic diffusions

This model describes a continuous evolution of processes $A,B,C$ along
solutions of a stochastic differential equation (SDE) [oeksendal] where $C$
acts autonomously and $A$ and $B$ hierarchically depend on each other, in
short:

The SDE is given by

$\begin{split}dC_{t}&=G(C_{t})dt+\sigma_{C}dW_{Ct}\\\ dB_{t}&=(\alpha
H(B_{t})-10(B_{t}-C_{t}))dt+\sigma_{B}dW_{Bt}\\\ dA_{t}&=(\alpha
H(A_{t})-5(A_{t}-B_{t})-5(A_{t}-C_{t}))dt+\sigma_{A}dW_{At}\\\ \end{split}$
(4.3)

where $G(x)=\begin{pmatrix}0\\\ -10(x_{2}^{3}-x_{2})\end{pmatrix}$ and
$H(x)=\begin{pmatrix}-(x_{1}^{3}-x_{1})\\\ -1\end{pmatrix}$.

The choice of the function $G$ ensures a metastable behaviour of $C$ between
regions around the values $-1$ and $1$ in the $x_{2}$, coordinate. In $x_{1}$
direction, the movement of $C$ is determined by noise with variance
$\sigma_{C}$. The movements of the other processes are partly governed by the
function $H$, which gives a metastable behaviour in the $x_{1}$-coordinates,
and by a difference function between themselves and $C$, respectively (Figure
5 for $\sigma_{A}=\sigma_{B}=0.01$). The movement of $A$ depends equally on
its difference to $C$ as on its difference to $B$. Since diffusion processes
are attracted to low values of their governing potential, this urges $B$ to
move into the direction of $C$. Furthermore, it urges $A$ into the directions
of both $C$ and $B$. The parameter $\alpha$ governs how autonomous the
processes $B$ and $C$ are.

Figure 5: Realisation of the system described by Eq. (4.3). $A$ in blue, $B$
orange, $C$ in yellow. We see the metastable behaviour in the $y$ coordinate
in $C$ that the other processes emulate. Parameters: $\alpha=5$,
$\sigma_{C}=diag((0.01,0.05))$, $\sigma_{B}=\sigma_{A}=0.01$.

We create realisations of the processes with the Euler-Maruyama scheme
[kloeden] with a time step of size $\Delta t=0.1$ for $1000$ time steps. The
parameter we use are $\alpha=5,\sigma_{C}=diag((0.01,0.05))$ For the noise in
the evolution of $A$ and $B$, we use $\sigma_{B}=\sigma_{A}=0.01,0.05$ and
$0.2$.

#### Dependency analysis

For SPA I, we use $K_{X}=K_{Y}=10$. Contrary to the previous example, we
compute the dependencies of the time-differences $\Delta A,\Delta B,\Delta C$
between time steps NW: Give reason. Instead of $\Lambda_{XY}^{(\Delta t)}$ we
therefore compute $\Lambda_{\Delta X\Delta Y}^{(\Delta t)}$ for $A,B,C=X,Y$.

For different values for the noise variances $\sigma_{A},\sigma_{B}$, we
create 50 realisations of the SDE Eq. (4.3) with the same parameters and
initial conditions and perform the analysis each time, computing a SPA I
solution anew, so that these solutions are generally different from each
other, to test the robustness of the method with regard to the SPA I solution.

The results of the dependency analysis then well reflect the hierarchical
dependencies between the three processes. An exemplary result is given in Eq.
(4.4) and Eq. (4.5). The statistics of the overall analysis are given in Table
1. It shows that in the vast majority of the realisations the more influential
variable was measured as the stronger influence. The average row variance
$\nu$ gives more lob-sided results with the minimal relative difference at
most $0.19$ for large noise but generally around $0.4$. For the Schatten-1
norm, the relative differences around mostly around $0.2$. Note that the
results seem quite independent of the strength of noise which loosens the
strict dependence of $A$ and $B$ on $C$.

$M_{\|\cdot\|_{1}}:\begin{tabular}[]{|l|l|l|l|l|l|}\hline\cr From $\downarrow$
to $\rightarrow$&$A_{t+\Delta t}$&$B_{t+\Delta t}$&$C_{t+\Delta t}$\\\
\hline\cr$A_{t}$&5.51&3.86&3.46\\\ \hline\cr$B_{t}$&5.29&4.11&2.94\\\
\hline\cr$C_{t}$&4.51&7.72&4.26\\\
\hline\cr\end{tabular}\Rightarrow\delta(M_{\|\cdot\|_{1}})=\begin{pmatrix}\hbox{\pagecolor[rgb]{1.00,1.00,1}0}&\hbox{\pagecolor[rgb]{1.00,0.71,0}-0.27}&\hbox{\pagecolor[rgb]{1.00,0.77,0}-0.23}\\\
\
\hbox{\pagecolor[rgb]{0.77,1.00,0}0.27}&\hbox{\pagecolor[rgb]{1.00,1.00,1}0}&\hbox{\pagecolor[rgb]{1.00,0,0.5}-0.61}\\\
\
\hbox{\pagecolor[rgb]{0.77,1.00,0}0.23}&\hbox{\pagecolor[rgb]{0,1.00,0}0.61}&\hbox{\pagecolor[rgb]{1.00,1.00,1}0}\end{pmatrix}$
(4.4) $M_{\nu}:\begin{tabular}[]{|l|l|l|l|l|l|}\hline\cr From $\downarrow$ to
$\rightarrow$&$A_{t+\Delta t}$&$B_{t+\Delta t}$&$C_{t+\Delta t}$\\\
\hline\cr$A_{t}$&0.039&0.021&0.018\\\ \hline\cr$B_{t}$&0.046&0.024&0.016\\\
\hline\cr$C_{t}$&0.039&0.065&0.026\\\
\hline\cr\end{tabular}\Rightarrow\delta(M_{\nu})=\begin{pmatrix}\hbox{\pagecolor[rgb]{1.00,1.00,1}0}&\hbox{\pagecolor[rgb]{1.00,0.49,0}-0.54}&\hbox{\pagecolor[rgb]{1.00,0.57,0}-0.53}\\\
\
\hbox{\pagecolor[rgb]{0.49,1.00,0}0.54}&\hbox{\pagecolor[rgb]{1.00,1.00,1}0}&\hbox{\pagecolor[rgb]{1.00,0.00,0}-0.75}\\\
\
\hbox{\pagecolor[rgb]{0.57,1.00,0}0.53}&\hbox{\pagecolor[rgb]{0.00,1.00,0}0.75}&\hbox{\pagecolor[rgb]{1.00,1.00,1}0}\
\end{pmatrix}$ (4.5)

$\sigma_{A},\sigma_{B}$ | Variables | Average $\delta(M_{\|\cdot\|_{1}})$ | Incorrect | Average $\delta(M_{\nu})$ | Incorrect
---|---|---|---|---|---
$0.01$ | $A,B$ | -0.21 | 0.12 | -0.36 | 0.2
$A,C$ | -0.22 | 0.08 | -0.56 | 0.08
$B,C$ | -0.46 | 0 | -0.83 | 0.00
$0.2$ | $A,B$ | -0.2 | 0.06 | -0.4 | 0.1
$A,C$ | -0.17 | 0.16 | -0.51 | 0.14
$B,C$ | -0.36 | 0 | -0.76 | 0
$0.5$ | $A,B$ | -0.2 | 0 | -0.4 | 0.02
$A,C$ | -0.12 | 0.14 | -0.43 | 0.12
$B,C$ | -0.30 | 0.02 | -0.71 | 0.04
$1$ | $A,B$ | -0.25 | 0.02 | -0.65 | 0
$A,C$ | -0.04 | 0.28 | -0.19 | 0.30
$B,C$ | -0.22 | 0.08 | -0.62 | 0.02
Table 1: Results of the dependency analysis for the diffusion processes in
Eq.(4.3). The third and fifth columns denote the average relative difference
between the first variable (in the first row $A$) and the second variable (in
the first row $B$) in the Schatten-1 norm and the average row variance. The
fourth and sixth columns denote the relative number of occurrences when one
variable was falsely identified as more influential. The average relative
differences are always negative, meaning that the first variable (with lower
influence) was correctly measures as less influential than the second
variable.

###  4.3 Multidimensional autoregressive processes

In order to demonstrate that dependency can be detected for processes whose
evolution depends not only on present but on past terms, too, we simulate
realisations of multidimensional linear autoregressive processes (AR)
[brockwell] in which some variables are coupled with others. An
$n$-dimensional linear AR($q$) process is a dynamical system of the form

$X_{t+1}=\sum\limits_{i=0}^{q-1}\phi_{i}X_{t-1}+\varepsilon_{t+1}$ (4.6)

where the $\phi_{i}\in\mathbb{R}^{n\times n}$ and $\varepsilon_{t+1}$ is a
stochastic term which we will set to be normally distributed with mean $0$ and
(positive semi-definite) covariance matrix $C\in\mathbb{R}^{n\times n}$.
We now consider AR($q$) processes of the form

$\begin{pmatrix}X_{t+1}\\\
Y_{t+1}\end{pmatrix}=\sum\limits_{i=0}^{q-1}\begin{pmatrix}\phi_{i}^{XX}&\phi_{i}^{YX}\\\
0&\phi_{i}^{YY}\end{pmatrix}\begin{pmatrix}X_{t-i}\\\
Y_{t-i}\end{pmatrix}+\varepsilon_{t+1}.$ (4.7)

Specifically, we let $X$ and $Y$ be variables in $\mathbb{R}^{4}$. Thus, the
$\phi_{i}^{XX},\phi_{i}^{YX}$ and $\phi_{i}^{YY}$ are matrices in
$\mathbb{R}^{4\times 4}$. Through the structure of the coefficient matrices we
impose that $X$ is influenced by $Y$ but not vice versa. We let $q=3$ and
construct the $\phi_{i}$ by drawing normally distributed values with mean $0$
and variance $0.1$. $\varepsilon_{t+1}$ is normally distributed with mean $0$
and variance $0.01$. We then create a realisation of length $T=1000$ for such
a process. We obtain the results shown in Table 2. This procedure is executed
50 times and the dependency analysis done on each realisation.

#### Dependency analysis

We choose $K_{X}=K_{Y}=10$ and compute $\Lambda_{XY}^{(\tau)}$ and
$\Lambda_{YX}^{(\tau)}$ for $\tau=1,3,10,50$ to investigate how the dependence
between the processes evolves with increasing time shift.

We can see in Table 2 that the stronger influence of $Y$ on $X$ is recovered
for $\tau=1$ and $\tau=3$. For $\tau=1$ the relative differences are stronger
which is in line with the fact that the AR coefficients $\phi_{3}$ were
selected to be smaller in magnitude than for $\phi_{1}$ so that there should
be stronger influence of $Y_{t}$ on $X_{t+1}$ than on $X_{t+3}$. Moreover, not
only the relative differences are smaller for $\tau=3$ but also the average
absolute number of the measures, again correctly indicating a smaller cross-
influence with bigger time shift. For $\tau\geq 10$, only negligible
differences can be seen. This is consistent with the construction of the
processes which include direct influence up to $\tau=3$.

$\tau$ |  | Average $\|\cdot\|_{1}$ | Average $\delta(M_{\|\cdot\|_{1}})$ | Average $\nu\cdot 10^{-2}$ | Average $\delta(M_{\nu})$
---|---|---|---|---|---
1 | $X_{t}\rightarrow Y_{t+1}$ | 1.94 | -0.19 | 1.67 | -0.58
$Y_{t}\rightarrow X_{t+1}$ | 2.43 | 0.19 | 4.40 | 0.58
3 | $X_{t}\rightarrow Y_{t+3}$ | 1.88 | -0.15 | 1.43 | -0.49
$Y_{t}\rightarrow X_{t+3}$ | 2.22 | 0.15 | 3.20 | 0.49
10 | $X_{t}\rightarrow Y_{t+10}$ | 1.82 | -0.02 | 1.27 | -0.06
$Y_{t}\rightarrow X_{t+10}$ | 1.86 | 0.02 | 1.38 | 0.06
50 | $X_{t}\rightarrow Y_{t+50}$ | 1.86 | 0.01 | 1.39 | 0.02
$Y_{t}\rightarrow X_{t+50}$ | 1.86 | -0.01 | 1.36 | -0.02
Table 2: Results for dependency measures of $X$ on $Y$ and vice versa for 50
realisations of processes of the form Eq. (4.7).

##  5 Real-World Example 1

NW: Vielleicht wäre ein reales Beispiel schön, das aus der Naturwissenschaft
kommt und bei der man schon über Einflüsse Bescheid weiß, die man mit dieser
Methode rekonstruieren kann.

##  6 Real-world example: basketball player movement

We will now apply the dependency measures to the movement of basketball
players during a game and quantify influences between players in the same
manner as in the previous examples. For this, we use player tracking data from
SportsVU from a game of the 2015/16 NBA season between the Dallas Mavericks
and the Cleveland Cavaliers, played on 12th January in Dallas. The data give
the $x$\- and $y$ coordinates of each player on the court in 25 time frames
per second for most of the 48 minutes of play. The ball data sometimes seem
out of sync with the positions of the players and is not always available,
therefore we were forced to omit it. The positions are measured in the unit
feet (ft). The basketball court has a rectangular shape with a width ($x$
axis) of 94 and a length ($y$ axis) of 50 ft. We install a coordinate system
whose origin is at the middle point of both axes. If a player is sitting on
the bench and not actively participating in the game, we assign to him the
coordinate $(-48.5,-27.5)$, which is slightly outside of the court.

In basketball, each team has five players on court. Typically, all five
players attack or defend simultaneously so that all ten players are in one
half of the court for several, around 20, seconds before moving into the other
half. Each of the two teams in a game has five players on court at all times.

Figure 6 shows the distribution of positions of Cleveland player LeBron James
depending on whether he is attacking (in the opponent’s half) or defending (in
his own team’s half) during the first twenty minutes of play and his position
on court over time. We can see that during attack James can mostly be found
around the 3-point-line and occasionally closer to the basket, including often
at the edge of the small rectangular area around the basket. On defense, he is
typically positioned slightly to the left and right of the basket.

Figure 6: Left: distribution of positions of LeBron James during the first
half of the game between the Cleveland Cavaliers and the Dallas Mavericks on
the 12th January 2016 depending on whether he is in his own team’s (defending)
or the opponents half (attacking). Right: $x$ and $y$ coordinates over time,
measured in ft.

###  6.1 Applying the dependency analysis to the basketball data

We now perform the dependency analysis on player coordinate data during the
first half of the game. We consider only the ten players in the starting
lineups of the teams. For the representation of each two-dimensional position
of a player, instead of solving SPA I, we choose landmark points in advance
and solely compute the barycentric coordinates of player coordinates with
respect to them. For clarity of both visualization and numerical computations,
we consider only the absolute value of the $x$-coordinates, meaning that we
reflect coordinates along the half-court line. For this reason, we can
consider only the right half of the court and use the following landmark
points,

$\Sigma=\begin{bmatrix}48.5&48.5&0&0&20&20&40\\\
-27.5&27.5&-27.5&27.5&15&-15&0\end{bmatrix}$ (6.1)

so that $K_{X}=K_{Y}=7$. The landmark points are depicted in Figure 7.

Figure 7: Landmark points chosen with respect to which the barycentric
coordinates of the player coordinates are computed.

To measure dependencies between each two players, we use those points in time
during the game, when both players were on court and compute
$\Lambda^{(\tau)}_{XY}$ for each pair of players $X,Y$. We choose
$\tau=1\text{sec}$, so that we investigate the amount of information of a
player $X$’s position for the position of a player $Y$ one second later.

Note that, as mentioned, basketball games are frequently interrupted for
various reason such as fouls or time-outs. We take this into account by
defining an event as the part of play between interruptions and denote the
number of events considered by $L$. We then construct the training data in the
form of multiple short time series, i.e., storing coordinates from the $k$th
events
$\Gamma^{X}_{k}=[\gamma^{X}_{k,1},\dots,\gamma^{X}_{k,T_{k}-\tau}],\Gamma^{Y}_{k}=[\gamma^{Y}_{k,1+\tau},\dots,\gamma^{Y}_{k,T_{k}}]$
and minimizing
$\|[\Gamma^{Y}_{1},\dots,\Gamma^{Y}_{L}]-\Lambda[\Gamma^{X}_{1},\dots,\Gamma^{X}_{L}]$.

Furthermore, we distinguish between which team is attacking since decisions of
players should be strongly influenced by whether they are attacking or
defending. We therefore define three different scenarios:

$\begin{split}\text{Dallas attacking }&\Leftrightarrow\text{ at least 9
players in Cleveland's half.}\\\ \text{Cleveland attacking
}&\Leftrightarrow\text{ at least 9 players in Dallas' half.}\\\
\text{Transition}&\Leftrightarrow\text{ otherwise.}\end{split}$

We omit the analysis on the transition phase.

##  7 Conclusion

##  Appendix A Appendix

######  Lemma A.1.

The maximal Schatten-1 norm of a column-stochastic $n\times m$ matrix is $m$.

###### Proof.

Since $\|\cdot\|_{1}$ is a norm, the triangle inequality holds and yields

$\|A+B\|_{1}\leq\|A\|_{1}+\|B\|_{1}.$

Thus, for $A=\sum\limits_{i}A_{i}$ it holds

$\|A\|_{1}\leq\sum\limits_{i}\|A_{i}\|_{1}.$

Note that $A$ can be written as
$A=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\tilde{A}^{ij}$ where

$(\tilde{A}^{ij})_{kl}\left\\{\begin{array}[]{ll}A_{ij}&(k,l)=(i,j)\\\
0&\,\textrm{else}\\\ \end{array}\right.$

A matrix with only one non-zero entry $a$ has only one non-zero singular value
that is equal $a$. This means that $\|\tilde{A}^{ij}\|_{1}=A_{ij}$. Thus,

$\|A\|_{1}\leq\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\|\tilde{A}^{ij}\|=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\tilde{A}^{ij}=m,$

because the entries in every column sum up to $1$ since $A$ is a column-
stochastic matrix. ∎

Proof of Lemma 3.1:

######  Lemma A.2.

Let $A\in\mathbb{R}^{n\times m}$ with $n\geq m$. Then the Schatten-1 norm of
$A$ obtains the maximal value $m$ if deletion of $n-m$ rows of $A$ yields an
$m\times m$ permutation matrix.

###### Proof.

Let $A$ be of the form described in the Lemma. Then by deletion of $n-m$ rows
we can derive a permutation matrix $P$. For those matrices, it holds that

$PP^{T}=Id.$

Thus, the singular values of $P$, which are the square roots of the
eigenvalues of $PP^{T}$ are given by the square roots of the eigenvalues of
the $m\times m$ identity matrix. These are, of course, given by
$\sigma_{1}=\dots=\sigma_{m}=1$. Thus, $\|P\|_{1}=m$.
All deleted rows must be identical to the zero-vector of length $m$, since the
column sums of $A$ have to be equal to $1$ and the column sums of $P$ are
already equal to $1$. Therefore, the singular values of $P$ are equal to the
singular values of $A$ and their sum is equal to $m$ because the sum of
singular values of a matrix cannot shrink by adding zero rows. This is because

$AA^{T}=\begin{pmatrix}P\\\
0\end{pmatrix}\begin{pmatrix}P^{T}&0\end{pmatrix}=\begin{pmatrix}PP^{T}&0\\\
0&0\end{pmatrix},$

whose eigenvalues are the eigenvalues of $PP^{T}$, i.e. the singular values of
$P$, and additional zeros. Thus, the sum of singular values does not change.
Since $m$ is the maximal value for $\|A\|_{1}$ by Lemma A.1, it holds that
$\|A\|_{1}=\|P\|_{1}=m$. ∎

Proof of Lemma 3.2:

######  Lemma A.3.

Let $A\in\mathbb{R}^{n\times m}$ with $n<m$. Then the Schatten-1 norm of $A$
is maximized if and only if $A$ contains an $n\times n$ permutation matrix and
the matrix of the remaining $m-n$ columns can be extended by $n$ columns to a
permutation matrix.

###### Proof.

Fehlt noch. ∎

Proof of Lemma 3.3:

######  Lemma A.4.

The Schatten-1 norm of a column stochastic $n\times m$-matrix $A$ is minimal
if and only if $A_{ij}\equiv\frac{1}{n}$ and in this case is equal to
$\sqrt{\frac{m}{n}}$.

###### Proof.

If all entries of $A$ are given by $\frac{1}{n}$, then

$(AA^{T})_{ij}=\frac{m}{n^{2}}\quad\text{for all }i,j=1,\dots,n.$ (A.1)

Then $AA^{T}$ has exactly one non-zero eigenvalue, since it is a rank-$1$
matrix. This is equal to $\frac{m}{n}$ (corresponding to the eigenvector
$(1,\dots,1)^{T}$), since

$AA^{T}\begin{pmatrix}1\\\ \vdots\\\
1\end{pmatrix}=\frac{m}{n^{2}}\begin{pmatrix}1&\dots&1\\\ \vdots&&\vdots\\\
1&\dots&1\end{pmatrix}\begin{pmatrix}1\\\ \vdots\\\
1\end{pmatrix}=\frac{m}{n^{2}}\begin{pmatrix}n\\\ \vdots\\\
n\end{pmatrix}=\frac{m}{n}\begin{pmatrix}1\\\ \vdots\\\ 1\end{pmatrix}.$

The square root of this eigenvalue, which is the only positive singular value,
is then $\sqrt{\frac{m}{n}}$. This yields

$\|A\|_{1}=\sqrt{\frac{m}{n}}$

The reason there cannot be a column-stochastic matrix $B$ with
$\|B\|_{1}<\sqrt{\frac{m}{n}}$ is the following: It holds that

$\|B\|_{1}\geq\|B\|_{2}\quad\text{where
}\|B\|_{2}:=\sqrt{\sum\limits_{i=1}^{min(n,m)}\sigma_{i}^{2}}\text{ is the
}\textit{Schatten-2-norm}$

For this, it holds

$\|B\|_{2}=\|B\|_{F}^{2}:=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}B_{ij}^{2},$
(A.2)

which is the squared Frobenius norm of $B$. Note that, if
$B_{ij}=\frac{1}{n}$, then

$\|B\|_{F}^{2}=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\frac{1}{n^{2}}=\frac{nm}{n^{2}}=\frac{m}{n}\Rightarrow\|B\|_{2}=\sqrt{\frac{m}{n}}.$
(A.3)

Assume that there exists a $B_{ij}<\frac{1}{n}$. Then since
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}B_{ij}=1$, there also exists
$B_{kj}>\frac{1}{n}$. Immediately, this increases the sum over the squared
entries of $B$. Therefore, $\|B\|_{F}$ is minimal only for the choice given
above. As a consequence, $\|B\|_{2}$ is minimal in this case and thus
$\sqrt{\frac{m}{n}}$ is the minimal value for $\|B\|_{1}$.
∎

Proof of Lemma 3.6:

######  Lemma A.5.

The average row variance of a column-stochastic $n\times m$-matrix $A$ obtains
the minimal value $0$ if and only if all columns are equal to each other.

###### Proof.

Trivially, if all columns in $A$ are identical, then the average row variance
of $A$ is $0$. Since the variance is a non-negative property from its
definition, this is the minimal value. If at least two values differ in a
column, then the average row variance immediately becomes positive. ∎

######  Corollary A.1.

For two processes $X$ and $Y$, let
$\Lambda^{(\tau)}_{ij}=\mathbb{P}(Y_{t+\tau}=S^{Y}_{i}|X_{t}=S^{X}_{j})$. Then
$\|\Lambda^{(\tau)}\|_{1}$ and $v(\Lambda^{(\tau)})$ are minimal if and only
if
$\mathbb{P}(Y_{t+\tau}=S^{Y}_{i}|X_{t}=S^{X}_{j})=\mathbb{P}(Y_{t+\tau}=S^{Y}_{i})$.

###### Proof.

If
$\mathbb{P}(Y_{t+\tau}=S^{Y}_{i}|X_{t}=S^{X}_{j})=\mathbb{P}(Y_{t+\tau}=S^{Y}_{i})$,
then all rows of $\Lambda^{(\tau)}$ are equal to each other. The result then
follows from Lemmas 3.3 and 3.6. ∎

Proof of Lemma 3.4:

######  Lemma A.6.

Let $A\in\mathbb{R}^{n\times m}$ with $n\geq m$. Then the average row variance
of a column-stochastic $n\times m$-matrix $A$ is maximal if (and only if) $A$
contains an $n\times n$ permutation matrix and the matrix of the remaining
$m-n$ columns can be extended by $n$ columns to a permutation matrix.

###### Proof.

The variance of a row that contains only values between $0$ and $1$ is
maximized if exactly one value is $1$ and all other entries are $0$. ∎


#  What’s for Lunch?
A systematic ordering of foods in the Soup-Salad-Sandwich phase space

Madelyn Leembruggen  Department of Physics, Harvard University, Cambridge, MA
02138, USA Caroline Martin  Harvard John A. Paulson School of Engineering and
Applied Sciences, Harvard University, Cambridge, MA 02138, USA

(April Fools’ Day 2022)

###### Abstract

The statistical physics of phase transitions has been hugely successful at
describing numerous natural, physical, and technological phenomena, and now a
rigorous examination of the phase space of the culinary regime is likewise
ripe for the picking. Despite great demand for the resolution of many
scientific debates over the taxonomy of food, past attempts have failed to
account for complex phase behavior and co-existence, and have thus left the
public hungry for a more substantial theory. By applying the principles of
statistical physics and thermodynamics, we here map out the complete phase
space of all culinary dishes and find three distinct phase regimes: Soup,
Salad, and Sandwich. We consider the effect of different state variables on
these phase boundaries, as well as regions of co-existence and triple points.
With this complete 3-dimensional phase diagram of all foods, we can
conclusively answer many bitter debates, including the imperishable question
“is a hotdog a sandwich?” The answer: yes.

##  1 Introduction

Phase transitions are ubiquitous across the physical world. By understanding
the underlying physics of phases, we have been able to make connections
between such disparate behavior as the ordering of spins in a lattice [1], the
melting of sea ice [2], and the firing of neurons in a brain [3]. The similar
behavior of such different systems which share a single universality class
enables not only physical comparison, but also increased predictive power for
such asymptotic phenomena as the critical exponents.

In recent years, the tools developed in the study of the phase behavior of
physical systems have been zestfully applied to food science, though
unfortunately with great controversy and dissent. Namely, the as-yet
unresolved scientific debate of “Is a hotdog a sandwich?” has periodically
sparked a resurgence of interest in this field. However, all efforts to this
point have floundered and fallen short of a complete examination of the full
phase behavior of foods, neglecting such important concepts as phase
transitions, co-existence, and critical points within the phase space. Without
this statistical physical framework, there have been grave misunderstandings
and misclassifications. In fact, many fierce disputes over the phase of a
particular dish can be solved simply by consider regions of phase co-
existence.

Here, we fully map out the phase space of the soup-salad-sandwich (Triple S)
transitions along the axes that we have identified as essential for studying
these three phases: effective temperature, effective pressure, and
carbohydrate enclosure. We then identify order parameters, including effective
density of the matrix and carbohydrate ordering.

Some might wonder, “why address this question at all?” In answer, we point to
the historic scientific progress and advancements made in the fields of
thermodynamics and statistical physics. Innumerable interesting properties of
systems have been explored and examined through the lens of phase transitions
and critical phenomena. The same progress can also be made in the vast menu of
culinary phenomena. By identifying distinct phases, we reveal common
properties for each phase, and although some of these characteristics might be
abstract, practical conclusions can also be drawn. For example, knowing the
phase of a dish may allow you to adequately prepare for your meal by
predicting the types of utensils which will be necessary, or the potential for
the dish as a finger food. More excitingly, studying the transitions between
phases may lead to bizarre and novel dishes which are incapable of existing
anywhere except at the phase boundaries.

Finally, we are motivated to carry out this analysis to achieve a more
universal understanding of dishes, one that rises above the arbitrary
limitations of human language. Often, different cultures will categorize
fundamentally similar foods in vastly disparate ways. While some types of
cuisine might not formally have “sandwiches”, they certainly have foods that
could, in a broader sense, be classified as such. Thus, although we frame our
problem in terms of American-centric definitions of “soup”, “salad”, and
“sandwich”, we aim to abstract away our colloquial connotations and provide
objective and quantitative definitions useful for characterizing these phases
and their critical behavior.

Figure 1: The full phase space of all foods can be broken down into three
distinct regions of phase space: soup, salad, and sandwich. These phase
regions depend on such physical parameters as effective temperature, effective
pressure, and carbohydrate enclosure.

###  1.1 Background

In many ways, the viral popularity of The Great Hotdog-Sandwich Debate
revolutionized the field of culinary taxonomy. Because it has never been fully
settled, interest in this question periodically re-emerges in the zeitgeist,
as demonstrated in Figure 2 (b). Previous attempts to formalize an answer to
this question have made some progress despite various flaws in logic and
scientific understanding. One important intellectual jump was the realization
that food exists in more categories than just sandwich alone. Broader
categories were introduced, such as the Soup-Salad-Sandwich (Triple S)
trichotomy [4] and the Cube Rule [5], which included such categories as Toast,
Sandwich, and Taco. Others have introduced the concept of foods existing in
regions of many dimensions, included proposed dimensions such as Soupiness,
Arrangement Entropy, and Ingredient Entropy [6]. While we acknowledge the
intellectual progress made in each of these attempts, we feel that none have
satisfyingly resolved the debate; of them, however, the authors (along with
many other citizen scientists participating in social media discourse) are
most compelled by the attempt to classify all foods as either a soup, salad,
or sandwich.

Despite our ostensibly innate human desire to sort and classify information,
past efforts to apply this Triple S trichotomy to the culinary space have
continued to lose steam. This has resulted in widespread skepticism in the
validity of the approach at all. In a Twitter poll conducted over a 24 hour
period with 406 responses, only $21.4\%$ of our respondents agreed that all
foods can be classified as a soup, salad, or sandwich (Fig. 2). Most of the
remaining $70.6\%$ of respondents rather jarringly disagreed. For example,
Twitter user @I_Am_Stan, who revealed they strongly disagreed, stated “Please
vote in this nonsense [poll]. There is only one correct choice.” Amongst even
the more mild dissenters, distaste was immediate, as reported by Twitter user
@BeckyDouglas: “I have clicked ‘disagree’ but I want you to know that this has
sparked a heated debate in the household this evening.”

Figure 2: Attempts of culinary taxonomy and phase categorization have
generated widespread public interest, though with limited success. (a) A
Twitter poll reveals broad distrust in the Triple S trichotomy, likely due to
misguided past attempts at applying these descriptors without a full
understanding of the Triple S phase behavior. (b) Many of the previous
research attempts can be traced back to an effort to resolve The Great Hotdog-
Sandwich Debate: “Is a hotdog a sandwich?”. This issue periodically recaptures
the zeitgeist, indicating strong public appetite for a conclusive answer.

Our public opinion poll (Figure 2 (a)) reveals a strong distrust of the Triple
S trichotomy. We draw the conclusion that, as a trichotomy, soup-salad-
sandwich is emphatically unfavorable; attempts to sort foods into three
distinct, non-overlapping bins are controversial and fail time and again.
However, these authors believe that all hope is not lost. We assert that
public discourse has been focusing on the wrong question. The question should
not be “can all foods be sorted”, but rather “can all foods be classified”.
Although subtle, the distinction of classification unfurls an entire phase
space for foods where soup, salad, and sandwich are phases separated by first
or second order transitions, with numerous regions of co-existence and
criticality. It is in these transitional regions that the most wonderful (and
sometimes the most appalling) dishes reside.

Other attempts to resolve this question of food classification have included
mapping all dishes onto a 3-dimensional space where each dish has its own
values for the soupiness, sandwichness, and saladness bounded between 0 and 1
[7]. While interesting, this effort misunderstands the role here of Soup,
Salad, and Sandwich. They are not state variables, they are distinct phases.
Because this previous work demonstrates a fundamental misunderstanding of
phase behavior, the result is a vast oversimplification of the Triple S space.

###  1.2 Methods

In this paper we take a different approach, and fully map the Triple S space
and phase transitions along three axes identified as essential for capturing
the appropriate phenomena of the space: effective temperature ($\tau$),
effective pressure ($\rho$), and carbohydrate enclosure ($\Theta$). By varying
these state variables, we can move into different regions of the phase space.
Importantly, in this work we only classify complete dishes– a dish defined
colloquially as something one would actually eat and prepared as such, not
just a random individual ingredient. This is because a single ingredient of a
dish111Twitter user @gooseus replied to our poll (Fig. 2) with the question,
“What about a single berry or steak alone? Is that a salad of one item? a
sandwich of organic tissues? or sandwich, holding all but one component?”,
which we assert are single ingredients that cannot hold a phase on their own.
should be considered as analogous to an individual spin. Individual spins do
not have a phase; rather, you must consider them as a collective group in
order to characterize their bulk behavior. Therefore you cannot qualify the
phase of single spin systems (i.e. single ingredient dishes). The individual
component must be coupled to other components in order for the system (dish)
as the whole to have a phase.

##  2 Soup-Salad phase diagram

First, we consider a projection of the phase space which involves only the
soup and salad phases. This is equivalent to a projection along the $\Theta=0$
plane, with effective temperature and pressure ($\tau$ and $\rho$) as the
remaining axes. We find that in this limit the soup-salad phase behavior maps
almost directly to the phase diagram for water.

The order parameter for the soup-salad transition is the density of the
embedding matrix. More explicitly, soups are a suspension of ingredients in a
liquid matrix, whereas salads are a suspension of ingredients in a gaseous
matrix. Thus, the soup-salad phase transition looks quite similar to the
liquid-gas transition for water, with the primary difference that the soup-
salad transition is a line of second order phase transitions for all $\tau$,
(i.e. the critical temperature is $\tau_{c}=0$).

Intuition for this distinction of soup from salad can be gained by considering
an argument about pore size and penetration depth of the embedding matrix. If
the bulkiest components of the dish form a scaffold, then the ability of the
remaining ingredients to efficiently permeate this scaffold determines the
phase of the dish. That is, if the constituents are completely incapable of
penetrating the bulky scaffold, the dish must be a salad because the embedding
matrix is considered to be air. On the other hand, complete penetration of the
scaffold with another of the dish’s ingredients must be a soup, because the
embedding matrix is liquid in nature and thus able to permeate the system
without flow restrictions.

Importantly, we identify the region of the soup-salad phase that is analogous
to solid water as a glassy soup. In the case of a glassy soup, the embedding
matrix is not air and yet the dish, at its serving temperature, does not flow
like a (liquid) soup. Glassy soups are separated from the soup phase by a line
of first order phase transitions as the embedding matrix melts from a glass to
a liquid. However, it should be noted that soups and glassy soups share an
overwhelming number of properties and should not be treated as separate
phases.

Figure 3: The Soup-Salad projection of the phase space (setting carbohydrate
enclosure $\Theta=0$) maps almost directly to the phase diagram for water,
where the density of the dish’s embedding matrix determines the phase of the
dish. A liquid matrix indicates a soup (a), here shown as a bowl of cereal. A
soup at low effective temperatures (frozen or solid embedding matrix) is
classified as a glassy soup (b), here shown as vanilla ice cream. If the
components of a dish are embedded in a gaseous matrix, they are in the salad
phase. The transition between the soup and salad phases is continuous, as
demonstrated in (c) where an identical set of ingredients is presented as a
chunky salad, a salsa, then finally as a gazpacho, smoothly transitioning from
one phase to the other as the components of the dish are chopped more finely.
The triple point of this phase diagram is embodied in an ice cream sundae (d),
a dish in which all three phases coexist. Finally a salad is defined as a
collection of components which are embedded in air, here shown as a bowl of
mixed nuts (e).

###  2.1 Case Studies

In addition to colloquially agreed upon soups, some other examples of soup
include smoothies, cold cereal, and hot oatmeal. Many types of salads are
already in the common vernacular, such as leafy green salad, fruit salad, and
tuna salad. Our results do indicate that many unexpected dishes are in fact
salads, including trail mix, a handful of M&M’s, a plate of chicken
wings222Thank you to Twitter user @supernova_mads for raising this question
and to Twitter user @biguyreacts who suggested “chicken skin as the crust so
chicken wings are a sandwich” but was unfortunately wrong, as will be
explained in more detail in Section 3., or even shredded cheese eaten straight
from the bag while you stand in front of your open fridge at 2:30 AM. Glassy
soups are actually quite prevalent; a few examples include ice cream,
chocolate, and hummus333Thank you to Twitter user @QuantumYakar for asking
about the phase of hummus.. Note that some dishes colloquially referred to as
salads (such as tuna, chicken, egg, and potato “salad”) often exist close to
the salad-glassy soup boundary and will experience a second order phase
transition to a glassy soup as a function of the amount of mayonnaise added to
the salad (note that mayonnaise on its own should be classified as a glassy
soup). Similarly, a bowl of cubed Jello sits near the salad-glassy soup
transition, as the Jello itself is a glassy soup, but the cubes in a bowl
constitute a salad (as the embedding matrix is air). On the other hand, a
Jello “salad” (which is comprised of fruits and vegetables embedded in a Jello
matrix) is actually not a salad at all, but rather a glassy soup (though as
the Jello matrix cools from a liquid into a gel, the Jello salad will cross
the soup-glassy soup boundary).

We now consider the second order phase transition which occurs across the
soup-salad boundary. On one side of the transition, we find such dishes as a
chunky salad with tomatoes, onions, and peppers with a dressing of lime juice
and olive oil. If the same ingredients are diced more finely, the same dish
might be called salsa, and salsas span a wide range from chunky to soupy. By
blending the same ingredients you can make gazpacho, which falls solidly in
the soup phase (Fig 3 (c)).

Finally, as an example of the triple point of the soup-salad phase space, we
present the ice cream sundae. It is built from both a soup (hot fudge) and
glassy soup (ice cream) with salad components (sprinkles or candy bits) and
arranged as a salad overall. In this narrow co-existence region, we find this
nonpareil of food, the pinnacle of culinary ingenuity, a shining example of
the edible richness discovered through the insights of phase behavior.

##  3 The Sandwich axis

We now introduce a third axis and state variable: the carbohydrate enclosure
($\Theta$) of the system. $\Theta$ is bounded from $[0,1]$, with $0$
indicating that the carbohydrate is completely enclosed by the other elements
of the dish (henceforth referred to as “elements”), and $1$ indicating the
state where all elements are contained entirely within the carbohydrate. Our
carbohydrate enclosure parameter is similar in spirit to the Cube Rule, which
was introduced and popularized in 2018 [5]. The Cube Rule differentiates
between sandwiches and not-sandwiches based on where the carb is oriented with
respect to the other elements. While we draw a great amount of inspiration
from the Cube Rule, we assert their definition of “sandwich” is in fact too
narrow, and most of the food categories they identify actually exist in the
sandwich phase. The following analysis will elucidate the breadth of the
sandwich phase, and provide compelling arguments for carbohydrate enclosure as
the proper parameter for this third axis.

Figure 4: (a) We introduce the state variable of the carbohydrate enclosure
$\Theta$, inspired by past work which considered the effect of carb enclosure
in food taxonomy: the Cube Rule [5]. We reject the differentiation of carb-
containing foods into these six categories and feel it is an unnecessary
complication of the phase space. (b) Instead, we introduce the parameter of
the the carbohydrate enclosure ($\Theta$). Here, as an introduction to the
role of $\Theta$ as a state variable, we consider the transition of the same
ingredients from salad (poke bowl) to sandwich (nigiri, sushi roll, and
onigiri) simply by varying the enclosure of carbohydrate (rice).

Interestingly, the sandwich region is not bounded by a constant carbohydrate
enclosure. Instead, the location of this phase transition depends also on the
effective pressure and temperature of the system. This is because the efficacy
of a carbohydrate at enclosing the elements depends on the effective
temperature and pressure of the elements. If the elements are in a soup phase
(high effective pressure and temperature), a carbohydrate enclosure close to
$1$ will be necessary for that dish to be in sandwich phase. In the more solid
region of the TP plane (generally lower effective pressure and temperature),
on the other hand, a dish can be in the sandwich phase with a relatively lower
value of $\Theta$, such as an open-face sandwich.

We capture these phenomena by taking slices at high, moderate, and low
effective pressures and examining the sandwich transition as a function of
effective temperature and carbohydrate enclosure only (See Fig. 5). We then
connect these projections to compile the full, 3-dimensional diagram of the
phase space (Fig. 1). In each of the slices, the critical $\Theta_{c}$ line
represents transitional sandwich foods. In certain regimes, the phase of the
dish might be ambiguous and difficult to classify analytically. Thus, in these
regions it is necessary that we to turn to an experimental probe of the
system. The design of the experiment is simple; namely, the research question
is: Can you pick it up by the carb and eat it? If so, the dish is a sandwich.

(a) The soup-sandwich space at high effective pressure

(b) The soup-salad-sandwich space at moderate effective pressure

(c) The salad-sandwich space at low effective pressure

Figure 5: 2-dimensional slices of the full 3-dimensional phase space as a
function of effective temperature and carbohydrate enclosure $\Theta$ for
high, moderate, and low effective pressure. These regions of phase space
reveal a rich dependence on these state variables, including a triple point of
soup-salad-sandwich co-existence.

###  3.1 Case Studies

At high effective pressure (Figure 5(a)), we find ourselves in a region of
phase space with a transition between soup and sandwich phases that depends
nonlinearly on the carbohydrate enclosure. At high effective temperature and
low carb enclosure, we have a soups, including dishes such as chicken noodle
soup; when this fluid at the same effective temperature and pressure is
enclosed completely by a carbohydrate, it becomes a sandwich, namely a soup
dumpling. At lower effective temperature, the necessary carbohydrate enclosure
for the experimental probe of sandwichness (i.e. can you pick it up and eat
it?) becomes much less due to the solid (or glassy) state of the internal
components. Here, an ice cream sandwich or mochi are both solidly within the
sandwich phase even at varying carbohydrate enclosure.

At moderate effective pressures (Figure 5(b)), there exists a triple point of
phase co-existence between all three phases in the Triple S space. This triple
point includes dishes such as an ice cream sundae in a waffle cone- a glassy
soup topped with soup and salad components and contained within a
carbohydrate. At lower carbohydrate enclosures, we find soups (such as sopa de
conchas) that transition to a salad dish (shells with pesto) as the liquid
penetrates less of the scaffold and the ingredient components are therefore
embedded in a less dense matrix. Finally, at very high carbohydrate
enclosures, the liquid sauce is entirely enclosed within the carb shell
(stuffed shells), allowing a courageous eater to potentially pick this
sandwich up and eat it.

At low effective pressure (Figure 5(c)), the behavior simplifies to a salad-
sandwich transition with linear dependence on carbohydrate enclosure. For an
example of a phase transition in this region, consider the constant effective
temperature dishes of a salad with croutons (salad), tostada (sandwich), and
taco (sandwich). This transition occurs simply due to the increasing carb
enclosure, until the internal components can be fully supported by the
enclosing carbohydrate.

###  3.2 The Effects of Carbohydrate Component Ordering

Upon introduction of carbohydrate enclosure, questions will naturally arise
regarding continuous encasing carbs versus discrete carb elements.
Particularly, one might wonder about the effects of the discrete carbs’ order
and symmetry on the phase of the system. It could be that there exists a yet
unexplored, fourth axis of the Triple S space. This fourth axis could describe
the ordering of the discrete carbs, mapping to the liquid crystal universality
class of phase transitions. For example, the smectic or nematic transition
could indicate a transition of the dish to a sandwich phase. A potentially
compelling case study that may justify further analysis is that of spaghetti,
which is solidly in the soup region when served with a lot of sauce, but is in
the salad region when served only with mountains of shakey parmesan cheese.
Cooked spaghetti could never fall in the sandwich phase because the carb
fragments are firmly disordered or, at best, nematic (due to their lack of
positional or orientational order). An additional case study may be lasagna,
which sits at the soup-salad-sandwich triple point, but only because of the
smectic ordering of the carb components. Baked ziti, however, has all the same
components but sits on the soup-salad transition line since the noodles are
disordered. While these cases all point to potentially rich areas of further
study, a full exploration is beyond the scope of this work.

##  4 Conclusion

It is tempting to get bogged down in the microscopic status of a dish’s
components when attempting to categorize the dish. However, at a certain
point, we must borrow the philosophy of renormalization, and account for the
fluctuations from the base, mean phase by zooming out and examining their
overall, cumulative effect. At the end of the day, we must put details aside
and honestly ask ourselves: “Is this a soup, a salad, or a sandwich?” Most of
the time, the answer you blurt out with the least amount of thought is the
most true. Concurrently we must allow for, and frankly encourage, a
willingness to straddle the divides between phases and permit simultaneous
truths. This is the delicacy of phase behavior: our ability to categorize the
intrinsically similar while also relishing the messiness of phase co-
existence, triple points, and criticality.

Finally, to remove any of the reader’s lingering uncertainty: yes, a hotdog is
a sandwich.

##  5 Acknowledgements

ML and CM would both like to thank Dr. Dinah Melk, Mr. Bigg Snoup Melk, Karrot
Cayke, Tayto Conqueso, S.N. Vinny, and E. Moureg for their thoughts and
feedback on this project (and for tolerating innumerable heated conversations
on the topic). This research was made possible by Good Soup™, Massaged Kale
Salad, and Hot Off the Panini Presses, Inc. The only competing interest to
declare is our actual research desperately begging us for the time we
dedicated to this paper instead. We apologize to our advisors, and promise to
them that we never worked on this paper during times that we would have been
otherwise productive.

## References

  * [1] John Cardy.  Scaling and Renormalization in Statistical Physics.  Cambridge University Press, April 1996.  Google-Books-ID: Wt804S9FjyAC.
  * [2] K. M. Golden, H. Eicken, A. L. Heaton, J. Miner, D. J. Pringle, and J. Zhu.  Thermal evolution of permeability and microstructure in sea ice.  Geophysical Research Letters, 34(16), 2007.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2007GL030447.
  * [3] Gašper Tkačik, Thierry Mora, Olivier Marre, Dario Amodei, Stephanie E. Palmer, Michael J. Berry, and William Bialek.  Thermodynamics and signatures of criticality in a network of neurons.  Proceedings of the National Academy of Sciences, 112(37):11508–11513, September 2015.  Publisher: Proceedings of the National Academy of Sciences.
  * [4] Soup Salad Sandwich (@LunchTaxonomy) / Twitter.  https://twitter.com/LunchTaxonomy.
  * [5] The Cube Rule.  https://cuberule.com/.
  * [6] Salad Theory.  https://saladtheory.github.io/.
  * [7] Soup-Salad-Sandwich Space.  http://sandwichspace.xyz.