to salt 1 and the reference-electrode potential from equation (15.6), one
finds a concise MacInnes equation,

$\displaystyle\frac{\vec{i}}{\kappa^{0}}=$
$\displaystyle-\nabla\Phi^{\ominus}+\frac{s_{\textrm{e}_{1}}\nabla\mu_{\textrm{e}_{1}}^{0}}{n_{\textrm{e}^{-}}Fz_{\textrm{e}^{-}}}+\left(\frac{s_{\textrm{e}_{2}}}{n_{\textrm{e}^{-}}Fz_{\textrm{e}^{-}}}-\frac{\kappa^{0}\xi}{F\|\textbf{z}\|}\right)\nabla\mu_{\textrm{e}_{2}}^{0},$
(16.65)

in agreement with Pollard and Newman’s equation 32 [30].

###  16.3 The cosolvent electrolyte

For a final illustration of our method of developing flux laws for
concentrated multicomponent electrolytes, we consider a four-species liquid
solution comprising two uncharged species, which we call a _cosolvent_ , and a
single salt. Electrolytic solutions of this type are used in lithium-ion
batteries, which typically rely on blends of linear and cyclic carbonates to
dissolve simple lithium salts [26]. We let subscripts $+$ and $-$ describe the
cation and anion, respectively, a subscript e, the salt they associate to
form, and let subscripts $0$ and o stand for the two solvents. As was the case
in the binary-electrolyte example, the stoichiometric coefficients for the
salt satisfy equation (16.1), where $\nu$ is the total salt stoichiometry.

A salt–charge basis for the cosolvent electrolyte is embodied by the
transformation matrix

$\textbf{Z}=\ \begin{bmatrix}1&0&0&0\\\ \ 0&1&0&0\\\ 0&0&-\frac{\nu
z_{-}}{z_{+}-z_{-}}&\frac{\nu z_{+}}{z_{+}-z_{-}}\\\
0&0&\frac{z_{+}}{\|\textbf{z}\|}&\frac{z_{-}}{\|\textbf{z}\|}\end{bmatrix}.$
(16.66)

The species concentrations in the neutral state can then be written as

$\left[\begin{array}[]{c}c_{0}^{0}\\\ c_{\textrm{o}}^{0}\\\ c_{+}^{0}\\\
c_{-}^{0}\end{array}\right]=\textbf{Z}\left[\begin{array}[]{c}c_{0}\\\
c_{\textrm{o}}\\\ c_{\textrm{e}}\\\
0\end{array}\right]=\left[\begin{array}[]{c}c_{0}\\\ c_{\textrm{o}}\\\
-\frac{z_{-}\nu}{z_{+}-z_{-}}\cdot c_{\textrm{e}}\\\
\frac{z_{+}\nu}{z_{+}-z_{-}}\cdot c_{\textrm{e}}\end{array}\right],$ (16.67)

parametric functions of the three component concentrations $c_{0}$,
$c_{\textrm{o}}$, and $c_{\textrm{e}}$.

The nonideal parts of the electroneutral component chemical potentials are
parametrised by a symmetric $2\times 2$ block matrix

$\bm{\Lambda}=\left[\begin{array}[]{cc}\Lambda_{00}&\Lambda_{0\text{o}}\\\
\Lambda_{0\text{o}}&\Lambda_{\text{o}\text{o}},\end{array}\right],$ (16.68)

whose three independent entries are functions of only two of the three
component fractions $y_{0}$, $y_{\textrm{o}}$, and $y_{\textrm{e}}$. Following
equation (8.3), the gradients of the electroneutral component chemical
potentials are expanded as

$\displaystyle\nabla\left[\begin{array}[]{c}\mu_{0}^{0}\\\
\mu_{\textrm{o}}^{0}\\\ \mu_{\text{e}}^{0}\end{array}\right]=\,$
$\displaystyle RT\nabla\ln\textbf{y}_{\nu}$ (16.72)
$\displaystyle+RT(\textbf{I}_{3}-\bm{\nu}\textbf{y}^{\top}_{\nu})\left[\begin{array}[]{cc}\bm{\Lambda}&\textbf{o}\\\
\textbf{o}^{\top}&0\end{array}\right]\nabla\textbf{y}_{\nu},$ (16.75)

where in this case $\bm{\nu}=[1,1,\nu]^{\top}$ and
$\textbf{y}_{\nu}=[y_{0},y_{\textrm{o}},y_{\textrm{e}}]^{\top}$.

To derive flux explicit laws, we again use the ansatz (14.1) to write

$\frac{c_{\textrm{T}}^{0}\textbf{M}_{\nu}}{RT}=\begin{bmatrix}\tfrac{c_{\text{o}}}{c_{0}\mathfrak{D}_{0\text{o}}}+\tfrac{\nu
c_{\text{e}}}{c_{0}\mathfrak{D}_{0\textrm{e}}}&-\tfrac{1}{\mathfrak{D}_{0\text{o}}}&-\tfrac{\nu}{\mathfrak{D}_{0\textrm{e}}}\\\
-\tfrac{1}{\mathfrak{D}_{0\text{o}}}&\tfrac{c_{0}}{c_{\text{o}}\mathfrak{D}_{0\text{o}}}+\tfrac{\nu
c_{\text{e}}}{c_{\textrm{o}}\mathfrak{D}_{\textrm{oe}}}&-\tfrac{\nu}{\mathfrak{D}_{\textrm{oe}}}\\\
-\tfrac{\nu}{\mathfrak{D}_{0\textrm{e}}}&-\tfrac{\nu}{\mathfrak{D}_{\textrm{oe}}}&\tfrac{\nu
c_{0}}{c_{\text{e}}\mathfrak{D}_{0\text{e}}}+\tfrac{\nu
c_{\text{o}}}{c_{\text{e}}\mathfrak{D}_{\textrm{oe}}}\end{bmatrix}.$ (16.76)

These entries straightforwardly relate to the Stefan–Maxwell coefficients
through

$\displaystyle\mathfrak{D}_{0\text{o}}=\mathscr{D}_{0\text{o}},$ (16.77)
$\displaystyle\mathfrak{D}_{0\textrm{e}}=\frac{(z_{+}-z_{-})\mathscr{D}_{0+}\mathscr{D}_{0-}}{\mathscr{D}_{0+}z_{+}-\mathscr{D}_{0-}z_{-}},$
(16.78)
$\displaystyle\mathfrak{D}_{\textrm{oe}}=\frac{(z_{+}-z_{-})\mathscr{D}_{\text{o}+}\mathscr{D}_{\text{o}-}}{\mathscr{D}_{\text{o}+}z_{+}-\mathscr{D}_{\text{o}-}z_{-}}.$
(16.79)

Observe that the second two definitions take forms similar to equation
(16.22), despite the addition of a solvent.

Taking the reference velocity to be that of the first solvent, i.e.,
$\bm{\psi}^{\prime}=\textbf{i}_{0}/c_{0}$, the diffusivity matrix and
migration coefficients assume the forms

$\overline{\textbf{L}}^{0}_{\nu}=\begin{bmatrix}0&0&0\\\
0&\mathscr{L}_{\textrm{o}}&\mathscr{L}_{\times}\\\
0&\mathscr{L}_{\times}&\mathscr{L}_{\textrm{e}}\end{bmatrix},\>\>\>\>\>\>\bm{\xi}^{0}=\left[\begin{array}[]{c}0\\\
\xi_{\textrm{o}}\\\ \xi_{\textrm{e}}\end{array}\right].$ (16.80)

The terms appearing in the diffusion matrix are understood as follows.
$\mathscr{L}_{\textrm{o}}$ denotes the Onsager diffusivity of the second
solvent relative to the velocity of the first, $\mathscr{L}_{\times}$ the
cross-diffusivity between the second solvent and the salt in a frame moving
with the first solvent, and $\mathscr{L}_{\textrm{e}}$ the diffusivity of the
salt relative to the first solvent’s velocity. Each of these parameters is
independent, and the $2\times 2$ sub-block of Onsager diffusivities is
positive definite. The set of migration coefficients depends on two
independent parameters $\xi_{\textrm{o}}$ and $\xi_{\textrm{e}}$. A
distinction from the previous cases is that a neutral species—the second
solvent—may carry current. Thus transference numbers as typically understood
do not suffice to write complete flux-explicit transport laws for this system.

By relation (14.2), the Onsager diffusivities that make up
$\overline{\textbf{L}}^{0}_{\nu}$ can be understood in terms of component
diffusivities, as

$\displaystyle\mathscr{L}_{\textrm{o}}=\frac{c_{\text{o}}\mathfrak{D}_{0\text{o}}(\mathfrak{D}_{\text{oe}}c_{0}+\mathfrak{D}_{0\text{e}}c_{\text{o}})}{c_{0}(\mathfrak{D}_{\text{oe}}c_{0}+\mathfrak{D}_{0\text{o}}\nu
c_{\text{e}}+\mathfrak{D}_{0\text{e}}c_{\text{o}})}$ (16.81)
$\displaystyle\mathscr{L}_{\textrm{e}}=\frac{c_{\text{e}}\mathfrak{D}_{0\text{e}}(\mathfrak{D}_{\text{oe}}c_{0}+\mathfrak{D}_{0\text{o}}\nu
c_{\text{e}})}{c_{0}(\mathfrak{D}_{\text{oe}}c_{0}+\mathfrak{D}_{0\text{o}}\nu
c_{\text{e}}+\mathfrak{D}_{0\text{e}}c_{\text{o}})}$ (16.82)
$\displaystyle\mathscr{L}_{\times}=\frac{\mathfrak{D}_{0\textrm{o}}\mathfrak{D}_{0\textrm{e}}c_{\text{o}}c_{\text{e}}}{c_{0}(\mathfrak{D}_{\textrm{oe}}c_{0}+\nu\mathfrak{D}_{0\textrm{o}}c_{\textrm{e}}+\mathfrak{D}_{0\textrm{e}}c_{\textrm{o}})}.$
(16.83)

Each of these expressions agrees within a prefactor with the thermodynamic
oxygen diffusivity, the thermodynamic electrolyte diffusivity, and the cross
diffusivity proposed to describe lithium/air-battery electrolytes by Monroe
[29].

The benefits of relying on two migration coefficients and directly writing
flux laws in the salt–charge basis are made clearest by this example. Equation
(12.13) reads for the two-solvent case as

$\displaystyle\vec{J}_{0}^{\,0}=\vec{0},$ (16.84)
$\displaystyle\vec{J}_{\textrm{o}}^{\,0}=-\mathscr{L}_{\textrm{o}}\frac{c_{\text{T}}^{0}\nabla\mu_{\textrm{o}}}{RT}-\mathscr{L}_{\times}\frac{c_{\text{T}}^{0}\nabla\mu_{\textrm{e}}}{RT}+\xi_{\textrm{o}}\frac{\vec{i}}{F\|z\|}$
(16.85)
$\displaystyle\vec{J}_{\textrm{e}}^{\,0}=-\mathscr{L}_{\times}\frac{c_{\text{T}}^{0}\nabla\mu_{\textrm{o}}}{RT}-\mathscr{L}_{\textrm{e}}\frac{c_{\text{T}}^{0}\nabla\mu_{\textrm{e}}}{RT}+\xi_{\textrm{e}}\frac{\vec{i}}{F\|z\|}.$
(16.86)

Equations (16.84)-(16.86) are equivalent to equations 14-16 in the paper by
Monroe [29]. These laws combine with the MacInnes equation for a reference
electrode, written through equations (12.17) and (15.6) as

$\vec{i}=-\kappa^{0}\nabla\Phi^{\ominus}+\frac{\kappa^{0}\textbf{s}_{\nu}^{\top}\nabla\bm{\mu}_{\nu}^{0}}{n_{\textrm{e}^{-}}Fz_{\textrm{e}^{-}}}-\frac{\kappa^{0}{\bm{\xi}^{0}}^{\top}\nabla\bm{\mu}_{\nu}^{0}}{F\left\|\textbf{z}\right\|},$
(16.87)

to offer a complete description of isothermal, isobaric mass transport in the
cosolvent electrolyte.

One could proceed as in the previous examples to elucidate how the migration
coefficients, diffusivities, and conductivity relate to the Stefan–Maxwell
coefficients. This would show, for example, that the migration coefficient
$\xi_{k}$ is expected to be proportional to component concentration
$c_{\nu,k}$, as observed by Monroe [29]. We refrain from completing this
exercise here because of the meteoric rise in complexity of the intermediate
equations. On account of the matrix inversion involved when moving between
force- and flux-explicit forms, writing out entries within both the
thermodynamic- and transport-property matrices becomes increasingly untenable
as the number of species increases. Nevertheless, this process was implemented
and found to produce relationships among parameters consistent with those
reported in detail elsewhere [29].

Matrix forms that connect flux- and force-explicit representations of
transport laws are compact and generally useful. The operations described in
prior sections are readily implemented with symbolic-manipulation software for
systems containing large numbers of species. Perhaps most important is that
the structures of the coefficient matrices within the electroneutral transport
laws are very general. As shown in these examples, structural knowledge can be
leveraged to create consistent sets of thermodynamic and transport
constitutive laws a priori in either the Stefan–Maxwell or Fick
representation, for electrolytes involving any number of components.

##  17 Numerical implementation: Hull cell

Recently, the authors presented finite-element methods to simulate
Onsager–Stefan–Maxwell transport of multiple uncharged species [8, 18]. The
efficacy of these approaches generally relies on the spectral structure of the
coefficient matrix M. A significant advantage of the electroneutral
Onsager–Stefan–Maxwell equations developed here is that they can be solved
easily within this numerical scheme. Indeed, it was discussed in sections 5
and 8.4 how Sylvester’s law of inertia guarantees that the spectral structure
of electrolyte-transport laws remains essentially intact, in that the
electroneutral transport matrix over the salt–charge basis,
$\textbf{M}_{Z}^{0}$, is symmetric positive semidefinite with one null
eigenvalue corresponding to the column $\textbf{c}^{0}_{Z}$.

In the prior papers [8] and [18], the numerical challenge posed by singular
transport matrices was resolved with a matrix augmentation. In the present
context, for a given reference velocity in the component space,
$\bm{\psi}^{\prime}_{Z}$, we solve the augmented electroneutral
Onsager–Stefan–Maxwell equations as

$-\vec{\nabla}\mu_{Z}=\textbf{M}_{Z}^{0,\psi^{\prime}}\left(\gamma\right)\vec{\textbf{n}}_{Z},$
(17.1)

where
$\textbf{M}_{Z}^{0,\psi^{\prime}_{Z}}\left(\gamma\right)=\textbf{M}_{Z}^{0}+\gamma\bm{\psi}^{\prime}_{Z}\bm{\psi}^{\prime\top}_{Z}$
proves to be a symmetric positive definite matrix whenever $\gamma\neq 0$ (see
reference [8]). These equations are coupled to the component balances (5.3)
and electroneutral charge balance (8.28).

We demonstrate how the framework created above can facilitate numerical
analysis by solving for the steady-state current distribution in a symmetric
electrode Hull cell containing a nonideal binary electrolytic solution, from
which cations are plated and stripped under potentiostatic conditions,
assuming that the electrode reactions have very low overpotential (fast
kinetics). To show how thermodynamic consistency and system nonideality can be
brought into simulations in a meaningful way, we take the electrolyte to be
lithium hexafluorophosphate (LiPF6) in ethyl–methyl carbonate (EMC)—a
formulation used in lithium-ion batteries—whose properties were measured as
functions of composition by Wang et al. [32, 33]. This electrolyte has total
stoichiometry $\nu=2$ and the species charges are
$\textbf{z}=\left[z_{0},z_{+},z_{-}\right]^{\top}=\left[0,1,-1\right]^{\top}$;
the transformation matrix to the salt-charge basis we employ is

$\textbf{Z}=\begin{bmatrix}1&0&0\\\ 0&1&1\\\
0&\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}\end{bmatrix},$ (17.2)

in which the first two rows correspond to the simple association equilibria

$\text{EMC}\rightleftharpoons\text{EMC}\>\>\>\>\>\>\text{and}\>\>\>\>\>\>\text{Li}^{+}+\text{PF}_{6}^{-}\rightleftharpoons\text{LiPF}_{6}.$
(17.3)

Wang’s correlations for the composition dependences of ionic conductivity
$\kappa$, Fickian diffusivity $D$, and Darken factor $\chi$ are expressed in
terms of the salt fraction $y_{\textrm{e}}$ and solvent fraction $y_{0}=1-\nu
y_{\textrm{e}}$ in Table 1.

Properties
---
Term |  Function |  Units
$\kappa$ |  $\left(48.93y_{\text{e}}^{3/2}-284.8y_{\text{e}}^{5/2}817.8y_{\text{e}}^{4}\right)^{2}$ |  $\text{S m}^{-1}$
$D$ |  $\left(4.998-29.96y_{\text{e}}+53.78y_{\text{e}}^{2}\right)\times 10^{-10}$ |  $\text{m}^{2}\text{s}^{-1}$
$\chi$ |  $1-18.38y_{\text{e}}^{\frac{1}{2}}+155.3y_{\text{e}}-450.6y_{\text{e}}^{\frac{3}{2}}+1506y_{\text{e}}^{\frac{5}{2}}$ |  $-$
$t_{+}^{0}$ |  $0.4107-1.487y_{\text{e}}+2.547y_{\text{e}}^{2}$ |  $-$
$\rho$ |  $1007.1+10^{5}\times\big{(}0.0180y_{\text{e}}-0.1946y_{\text{e}}^{2}+1.960y_{\text{e}}^{3}-7.008y_{\text{e}}^{4}+8.004y_{\text{e}}^{5}\big{)}$ |  $\text{g}\>\text{L}^{-1}$
Table 1: Phenomenological property data necessary to specify the transport
matrix M and thermodynamic factor $\bm{\Lambda}$ for solutions of lithium
hexafluorophosphate in ethyl–methyl carbonate at room temperature, derived
from data reported by Wang et al. [32, 33].

Note that $D$ and $\chi$ determine the properties $\Lambda_{00}$ and
$\mathfrak{D}$ defined in section 16.1 through

$\displaystyle\mathfrak{D}=\frac{D}{\chi},$ (17.4)
$\displaystyle\Lambda_{00}=\frac{\chi-1}{\nu y_{\textrm{e}}y_{0}}.$ (17.5)

The migration coefficient $\xi$ is determined by the cation transference
number relative to the solvent velocity $t_{+}^{0}$ through equation (16.26),
which reduces to

$\xi=\frac{1-2t_{+}^{0}}{\sqrt{2}}$ (17.6)

for this electrolyte system. Finally, it should be noted that the composition
dependence of $\rho$ appearing in Table 1 was not reported by Wang et al., but
instead derives from a least-squares fit of the densitometric data reported in
the supplementary information of those papers as a function of mass fraction
[32, 33]. For this system, $\rho\left(y_{\textrm{e}}\right)$ acts as the
volumetric equation of state. When using density to compute properties
relevant to concentration from $\rho$, the molar mass of EMC was taken to be
$\overline{m}_{0}=104.105$ $\text{g}\,\text{mol}^{-1}$, and that of LiPF6,
$\overline{m}_{\text{e}}=151.905$ $\text{g}\,\text{mol}^{-1}$. After inferring
the component diffusivity and migration coefficients by inserting the
correlations from table 1 into equations (17.4) and (17.6), equations
(16.32)-(16.34) were used locally to assemble the matrix $\textbf{M}_{Z}$.
Simulations were performed assuming an equilibrium electrolyte composition of
$\left\langle y_{\text{e}}\right\rangle=0.15$ (ca.
$1.9~{}\textrm{mol}\,\,\textrm{L}^{-1}$).

To model potentiostatic electrodes with fast heterogeneous reaction kinetics,
the voltages at the positive electrode (p) and negative electrode (n) were
respectively assumed to be $\text{V}^{\text{p}}=10~{}\textrm{mV}$ and
$\text{V}^{\text{n}}=0$. At each electrode surface, we applied the Robin
boundary condition

$-\left.\left(\vec{i}\cdot\vec{n}\right)\right|_{k}=\frac{Fz_{+}i_{0}}{RT}\left(\text{V}^{k}-\left.\Phi_{Z}^{0}\right|_{k}\right);$
(17.7)

here $\vec{n}$ is an outward unit normal vector, anodic overpotentials are
positive, and $k\in\left\\{\textrm{p},\textrm{n}\right\\}$. In simulations the
exchange-current density $i_{0}$ was chosen to be very large
($1~{}\textrm{A}\,\textrm{cm}^{-2}$) so as to ensure that the solution-phase
electrode potentials $\Phi_{\textrm{Z}}^{0}$ matched $V^{\textrm{a}}$ or
$V^{\textrm{c}}$ to an appropriate degree of precision ($4$ digits). To relate
the electrode currents to interfacial fluxes, the flux of lithium ions was
taken to be

$\left.\left(\vec{N}_{\text{Li}^{+}}\cdot\vec{n}\right)\right|_{k}=\frac{1}{Fz_{+}}\left.\left(\vec{i}\cdot\vec{n}\right)\right|_{k}$
(17.8)

at electrode surface $k\in\left\\{\textrm{p},\textrm{n}\right\\}$; all other
boundary conditions were set as homogenous Neumann conditions.

The stripping and plating of lithium at the boundary as expressed in (17.8)
implies a transfer of mass across the solution. It would therefore be
inconsistent to take the reference velocity as the mass-average velocity and
then assume this to be zero. Unfortunately, it is also not possible to solve
for the distribution of mass-average velocity in two dimensions, because
closure of problems in systems with more than one spatial dimension requires
consideration of momentum continuity [34]. Instead we follow Newman and choose
the reference velocity to be that of the solvent, and thence assume that
velocity to equal zero uniformly. Practically, this amounts to a formal
neglect of the Faradaic convection phenomenon.

A steady–state simulation of the potentiostatic Hull cell was performed with
Firedrake software [35], using the MUMPS direct linear solver [36, 37] via
PETSc [38, 39]. The mesh of the geometry was constructed using the Gmsh
software [40]. Each linear system had $183,831$ degrees of freedom. The finite
elements were chosen as $\text{CG}_{2}$ for component mole fractions/voltage
and $\text{DG}_{1}$ for fluxes/current. The computed current distribution is
plotted in Figure 1. Observe that the current concentrates at the upper corner
because of the obtuse angle with an insulating surface, and is dispersed at
the lower corner, where the angle is acute, as might be expected from
electrostatic analysis of the wedge geometry in cylindrical coordinates [1].

Figure 1: Steady-state current density $|\vec{i}|$ within a Hull cell with
reversible electrode reactions subjected to a constant voltage drop. The
streamlines denote trajectories of current flow from the positive electrode
(left) to the negative electrode (right).

This example illustrates the potential of the structurally electroneutral
Onsager–Stefan–Maxwell framework to model concentrated electrolytes with full
composition dependence of all independent thermodynamic and transport
properties, on complex geometries. Similarly we might consider modelling
molten salts, or cosolvent electrolytes, or transient dynamics of any of
these. Subject to availability of the necessary property data, the scope of
applications can be extended much further.

## Conclusion

We have detailed how local electroneutrality affects Onsager–Stefan–Maxwell
mass transport within electrolytic solutions. Choosing a salt–charge basis
according to the process we outline in section 2 allows one to establish
constitutive equations for component chemical-potential gradients in terms of
a symmetric thermodynamic-factor matrix. Crucially, the use of a salt–charge
basis also obviates the need to impose local electroneutrality as an algebraic
constraint on the electrolytic mass-transport problem, allowing it to be woven
into the parameter space of the differential governing system. This leads to
the electroneutral Onsager–Stefan–Maxwell equations (8.34), which phrase $n-1$
component chemical potential gradients and the salt–charge potential gradient
explicitly in terms of $n-1$ uncharged component fluxes and current density.
Salt–charge potential relates in a straightforward way to the electric
potential measured by a reference electrode of a given kind.

A process was established to invert the electroneutral Onsager–Stefan–Maxwell
laws into the flux-explicit form familiar from concentrated solution theory,
allowing the identification of standard transport properties such as
conductivity, component diffusivity, and species transference numbers. We
introduced migration coefficients, a useful alternative that remedies the
issue that transference numbers as commonly understood cannot describe all the
electro-osmotic interactions in systems containing multiple uncharged species.
Several examples were given to substantiate how the salt–charge basis concept
enables the direct formulation of flux-explicit constitutive laws for mass and
charge transport in locally electroneutral concentrated electrolytic
solutions.

Establishing a salt–charge basis allows one to formulate the transport
equations for locally electroneutral electrolytic solutions in a way that
retains the spectral structure afforded by the Onsager–Stefan–Maxwell laws for
nonelectrolytes. The utility of this perspective was underlined by using a
code designed for nonelectrolytes to simulate electroneutral transport within
a nonideal binary electrolytic solution with composition-dependent properties
in a Hull cell.

## Acknowledgements

This work was supported by the Engineering and Physical Sciences Research
Council Centre for Doctoral Training in Partial Differential Equations:
Analysis and Applications (EP/L015811/1), the Engineering and Physical
Sciences Research Council (grants EP/R029423/1 and EP/W026163/1); the
Clarendon fund scholarship; and the Faraday Institution Multiscale Modelling
project (subaward FIRG025 under grant EP/P003532/1).

## References

  * Newman and Thomas-Alyea [2004] J. Newman, K. E. Thomas-Alyea, Electrochemical Systems, John Wiley & Sons, Hoboken, New Jersey, 2004.
  * Newman et al. [1965] J. Newman, D. Bennion, C. W. Tobias,  Mass transfer in concentrated binary electrolytes,  Berichte der Bunsengesellschaft für physikalische Chemie 69 (1965) 608–612.
  * Doyle et al. [1993] M. Doyle, T. F. Fuller, J. Newman,  Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell,  Journal of The Electrochemical Society 140 (1993) 1526–1533.
  * Doyle et al. [1996] M. Doyle, J. Newman, A. S. Gozdz, C. N. Schmutz, J. Tarascon,  Comparison of Modeling Predictions with Experimental Data from Plastic Lithium Ion Cells,  Journal of the Electrochemical Society 143 (1996) 1890–1903.
  * Hirschfelder et al. [1954] J. O. Hirschfelder, C. F. Curtiss, R. B. Bird, The Molecular Theory of Gases and Liquids, John Wiley & Sons, New York, 1954.
  * Monroe et al. [2015] C. W. Monroe, D. R. Wheeler, J. Newman,  Nonequilibrium linear response theory: application to Onsager–Stefan–Maxwell diffusion,  Industrial & Engineering Chemistry Research 54 (2015) 4460–4467.
  * Onsager [1945] L. Onsager,  Theories and problems of liquid diffusion,  Annals of the New York Academy of Sciences 46 (1945) 241–265.
  * Van-Brunt et al. [2022] A. J. Van-Brunt, P. E. Farrell, C. W. Monroe,  Augmented saddle point formulation of the steady-state Stefan–Maxwell diffusion problem,  IMA Journal of Numerical Analysis 42 (2022) 3272–3305.
  * Debye and Hückel [1954] P. J. W. Debye, E. Hückel,  On the theory of electrolytes. I. Freezing point depression and related phenomena, _trans._ from _Physikalische Zeitschrift_ 24 (1923) 185-206,  in: The Collected Papers of Peter J. W. Debye, Interscience Publishers, Inc., London, 1954, pp. 217–263.
  * Newman and Chapman [1973] J. Newman, T. W. Chapman,  Restricted diffusion in binary solutions,  AIChE Journal 19 (1973) 343–348.
  * Pollard and Newman [1979] R. Pollard, J. Newman,  Transport equations for a mixture of two binary molten salts in a porous electrode,  Journal of The Electrochemical Society 126 (1979) 1713–1717.
  * Crothers et al. [2020a] A. R. Crothers, R. M. Darling, A. Kusoglu, C. J. Radke, A. Z. Weber,  Theory of multicomponent phenomena in cation-exchange membranes: Part II. transport model and validation,  Journal of The Electrochemical Society 167 (2020a) 013548\.
  * Crothers et al. [2020b] A. R. Crothers, R. M. Darling, D. I. Kushner, M. L. Perry, A. Z. Weber,  Theory of multicomponent phenomena in cation-exchange membranes: Part III. transport in vanadium redox-flow-battery separators,  Journal of The Electrochemical Society 167 (2020b) 013549\.
  * Guggenheim [1928] E. A. Guggenheim,  The conceptions of electrical potential difference between two phases and the individual activities of ions,  The Journal of Physical Chemistry 33 (1928) 842–849.
  * Monroe and Delacourt [2013] C. W. Monroe, C. Delacourt,  Continuum transport laws for locally non-neutral electrolytes,  Electrochimica Acta 114 (2013) 649–657.
  * Goyal and Monroe [2021] P. Goyal, C. W. Monroe,  Thermodynamic factors for locally non-neutral, concentrated electrolytic fluids,  Electrochimica Acta 371 (2021) 137638.
  * Guggenheim [1985] E. A. Guggenheim, Thermodynamics–An advanced treatment for chemists and physicists, North-Holland, Amsterdam, 1985\.
  * Van-Brunt et al. [2022] A. J. Van-Brunt, P. E. Farrell, C. W. Monroe,  Consolidated theory of fluid thermodiffusion,  AIChE Journal 68 (2022) e17599.
  * Monroe and Newman [2006] C. W. Monroe, J. Newman,  Onsager reciprocal relations for Stefan–Maxwell diffusion,  Industrial & Engineering Chemistry Research 45 (2006) 5361–5367.
  * Monroe and Newman [2009] C. W. Monroe, J. Newman,  Onsager’s shortcut to proper forces and fluxes,  Chemical Engineering Science 64 (2009) 4804–4809.
  * MacInnes [1961] D. A. MacInnes, The Principles of Electrochemistry, Dover books on advanced mathematics, Dover, 1961.
  * Miller [1959] D. G. Miller,  Ternary isothermal diffusion and the validity of the Onsager reciprocity relations,  The Journal of Physical Chemistry 63 (1959) 570–578.
  * Fong et al. [2020a] K. D. Fong, J. Self, B. D. McCloskey, K. A. Persson,  Ion correlations and their impact on transport in polymer-based electrolytes,  Macromolecules 54 (2020a) 2575–2591.
  * Fong et al. [2020b] K. D. Fong, H. K. Bergstrom, B. D. McCloskey, K. K. Mandadapu,  Transport phenomena in electrolyte solutions: Nonequilibrium thermodynamics and statistical mechanics,  AIChE Journal 66 (2020b) e17091.
  * Bard and Faulkner [2000] A. J. Bard, L. R. Faulkner, Electrochemical Methods: Fundamentals and Applications, 2nd Edition, John Wiley & Sons, Incorporated, 2000\.
  * Xu [2004] K. Xu,  Nonaqueous liquid electrolytes for lithium-based rechargeable batteries,  Chemical Reviews 104 (2004) 4303–4418.
  * Valøen and Reimers [2005] L. O. Valøen, J. N. Reimers,  Transport properties of lipf6-based li-ion battery electrolytes,  Journal of The Electrochemical Society 152 (2005) A882–A891.
  * Weber [2008] A. Z. Weber,  Gas-crossover and membrane-pinhole effects in polymer-electrolyte fuel cells,  Journal of the Electrochemical Society 155 (2008) B521–B531.
  * Monroe [2017] C. W. Monroe,  Does oxygen transport affect the cell voltages of metal/air batteries?,  Journal of The Electrochemical Society 164 (2017) E3547–E3551.
  * Pollard and Newman [1979] R. Pollard, J. Newman,  Transport equations for a mixture of two binary molten salts in a porous electrode,  Journal of The Electrochemical Society 126 (1979) 1713–1717.
  * Doyle and Newman [1997] M. Doyle, J. Newman,  Analysis of capacity rate data for lithium batteries using simplified models of the discharge process,  Journal of Applied Electrochemistry 27 (1997) 846–856.
  * Wang et al. [2020] A. A. Wang, T. Hou, M. Karanjavala, C. W. Monroe,  Shifting-reference concentration cells to refine composition-dependent transport characterization of binary lithium-ion electrolytes,  Electrochimica Acta 358 (2020) 136688.
  * Wang et al. [2021] A. A. Wang, A. B. Gunnarsdóttir, J. Fawdon, M. Pasta, C. P. Grey, C. W. Monroe,  Potentiometric mri of a superconcentrated lithium electrolyte: Testing the irreversible thermodynamics approach,  ACS Energy Letters 6 (2021) 3086–3095.
  * Liu and Monroe [2014] J. Liu, C. W. Monroe,  Solute-volume effects in electrolyte transport,  Electrochimica Acta 135 (2014) 447–460.
  * Rathgeber et al. [2016] F. Rathgeber, D. A. Ham, L. Mitchell, M. Lange, F. Luporini, A. T. T. McRae, G. Bercea, G. R. Markall, P. H. J. Kelly,  Firedrake: Automating the finite element method by composing abstractions,  ACM Transactions on Mathematical Software (TOMS) 43 (2016) 1–27.
  * Amestoy et al. [2001] P. R. Amestoy, I. S. Duff, L. J.-Y., J. Koster,  A fully asynchronous multifrontal solver using distributed dynamic scheduling,  SIAM Journal on Matrix Analysis and Applications 23 (2001) 15–41.
  * Amestoy et al. [2006] P. R. Amestoy, A. Guermouche, J.-Y. L’Excellent, S. Pralet,  Hybrid scheduling for the parallel solution of linear systems,  Parallel Computing 32 (2006) 136–156.
  * Balay et al. [2019] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Karpeyev, D. Kaushik, M. G. Knepley, D. May, L. C. McInnes, R. T. Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, PETSc Users Manual, Technical Report ANL-95/11 - Revision 3.11, 2019\.
  * Balay et al. [1997] S. Balay, W. D. Gropp, L. C. McInnes, B. F. Smith,  Efficient management of parallelism in object oriented numerical software libraries,  in: E. Arge, A. M. Bruaset, H. P. Langtangen (Eds.), Modern Software Tools in Scientific Computing, Birkhäuser Press, 1997, pp. 163–202.
  * Geuzaine and Remacle [2009] C. Geuzaine, J. F. Remacle,  Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities,  International Journal for Numerical Methods in Engineering 79 (2009) 1309–1331.


2022

1]Department of Computer Science & Engineering, South Dakota Mines, Rapid
City, SD

2]Department of Mathematics, South Dakota Mines, Rapid City, SD

3]Naval Surface Warfare Center, Dahlgren Division, Dahlgren, VA

# Forecasting Multilinear Data via Transform-Based Tensor Autoregression

Jackson Cates  Randy C. Hoover  Kyle Caudle  Cagri Ozdemir  Karen Braman
David Marchette  [  [  [

###### Abstract

In the era of big data, there is an increasing demand for new methods for
analyzing and forecasting 2-dimensional data. The current research aims to
accomplish these goals through the combination of time-series modeling and
multilinear algebraic systems. We expand previous autoregressive techniques to
forecast multilinear data, aptly named the $\mathcal{L}$-Transform Tensor
autoregressive ($\mathcal{L}$-TAR for short). Tensor decompositions and
multilinear tensor products have allowed for this approach to be a feasible
method of forecasting. We achieve statistical independence between the columns
of the observations through invertible discrete linear transforms, enabling a
divide and conquer approach. We present an experimental validation of the
proposed methods on datasets containing image collections, video sequences,
sea surface temperature measurements, stock prices, and networks.

The current research was supported in part by the Department of the Navy,
Naval Engineering Education Consortium under Grant No. (N00174-19-1-0014) and
the National Science Foundation under Grant No. (2007367). Marchette was
funded by the NSWC Naval Innovative Science and Engineering (NISE) program.
Any opinions, findings, and conclusions or recommendations expressed in this
material are those of the authors and do not necessarily reflect the views of
the Naval Engineering Education Consortium or the National Science Foundation.
Portions of this work were presented in part at the 2021 IEEE/ACM
International Conference on Machine Learning and Applications art:cates2021 .

##  1 Introduction

Forecasting is known to be among the most challenging and problematic problems
within machine learning. It involves extrapolation — prediction of the future
from only past data de200625 . There are numerous methods that exist to meet
the challenges of forecasting. Some of the more classical forecasting
techniques include Box-Jenkins Autoregressive Integrated Moving Average
(ARIMA) book:Box  modeling and exponential smoothing Boxexp ; Brown1 ; Brown2
. Other novel methods have been used that provide forecasts based on
historical pattern matching. This method of forecasting, referred to as “Flow
Field” forecasting, bases forecasts on the previous slopes and positions in
the data record, similar to a slope field solution of a differential equation
Frey ; Caudle1 ; Caudle2 ; Caudle3 ; Caudle4 . More recently, recurrent neural
networks have also provided very competitive forecasts Haykin ; Hill , with
the cost of interpretablity.

In the context of the current work, we will focus on the autoregressive (AR)
portion of ARIMA book:Box  modeling and demonstrate how such methods can be
extended to model multilinear observations. Because our aim is to forecast in
multiple steps in the future, we excluded the moving average (MA) portion from
the ARIMA process as it generally contains unobservable error terms. In this
model, future values are forecasted using a linear (or multilinear)
combination of previous time series observations. The number of previous
values (also known as lags) that are used to forecast the present value is
known as the “order” of the model. For example, given the univariate AR model
of order $p$,

$y_{t}=\beta+\alpha_{1}y_{t-1}+\alpha_{2}y_{t-2}+...+\alpha_{p}y_{t-p}+\epsilon_{t},$
(1)

our goal is to estimate the model parameters
$\theta=\\{\alpha_{1},\alpha_{2},\dots,\alpha_{p},\beta\\}$, from the prior
observations $y_{j}\in\mathbb{R}$ ($j=1,\dots,n$), where in general $n>>p$.

This can be extended to a multivariate times series, where the observations
are represented as a vector. For example, given the multivariate Vector
Autoregressive (VAR) Asteriou  model of order $p$,

$\mathbf{y}_{t}=\mathbf{c}+A_{1}\mathbf{y}_{t-1}+...+A_{p}\mathbf{y}_{t-p}+\mbox{\boldmath{$\epsilon$}}_{t},$
(2)

where the collection of observations are $\mathbf{y}_{j}\in\mathbb{R}^{k}$
($j=1,2,\dots,n$) are vectors. Similar to the univariate case, the goal is to
estimate the model parameters
$\mbox{\boldmath{$\theta$}}=\\{\boldmath{A}_{1},\boldmath{A}_{2},\dots,\boldmath{A}_{p},\mathbf{c}\\}$,
$\boldmath{A}_{i}\in\mathbb{R}^{k\times k},\mathbf{c}\in\mathbb{R}^{k}$.
Again, it is generally assumed that $n>>p$.

In order to forecast 2-dimensional observations, their representation will
need to be viewed as a lateral slice of a tensor (e.g. ${\cal
Y}_{t}\in\mathbb{R}^{\ell\times 1\times m}$) 111Note: It’s customary in the
literature to represent tensors with upper-case calligraphic letters. instead
of a vector. Tensor in this context is a multi-dimensional array, often
referred to as $n$-mode or $n$-way array as defined in Section 2. Dynamic
networks, video sequencing, correlated image sets, and distributed sensing are
specific examples where tensor-based forecasting are of upmost importance. In
NIPS2013_5117  the authors develop a method to forecast higher-order tensors
based on the Tucker decomposition and $n$-mode products (referred to as
multilinear dynamical systems (MLDS)) art:Tucker66 ; art:Lathauwer00 . The
MLDS approach (based on dynamical systems theory and system identification
methods) was extended in art:weijun2018  by transitioning from the Tucker
decomposition to a recently defined tensor product based on discrete
transforms and mod-$n$ convolution, referred to as the ${\cal L}$-transform
KilMP08 ; KilM09 ; Bra10 ; art:Hoover11 ; art:Hoover18 ; ozdemir20212dtpca ;
ozdemir2021fast  (the details of which are outlined in Section 2). While both
methods outlined in NIPS2013_5117  and art:weijun2018  (MLDS and ${\cal
L}$-MLDS respectively) show promise, they are both based on multilinear
dynamical systems modeling as opposed to an autoregressive model as defined
above. In other words, they attempt to find a single state-transition tensor
to obtain their forecast.

The contributions of the current work are twofold: 1) We extend the results in
art:weijun2018  by transitioning from a traditional dynamical systems model to
an autoregressive model. Building on art:Kilmer13 ; art:Hao13 ; art:Hoover11 ;
art:Hoover18 , and the ${\cal L}$-transform outlined in art:weijun2018  we
show that we can extend a VAR model by estimating the model parameters
$\Theta=\\{{\cal A}_{1},{\cal A}_{2},\dots,{\cal A}_{p},{\cal C}\\}$ of the
tensor autoregressive model (${\cal L}$-TAR)

${\cal Y}_{t}={\cal C}+{\cal A}_{1}\bullet{\cal Y}_{t-1}+\cdots+{\cal
A}_{p}\bullet{\cal Y}_{t-p}+{\cal E}_{t},$

where $\bullet$ denotes the ${\cal L}$-product outlined in definition 3,
${\cal A}_{i}\in\mathbb{R}^{\ell\times\ell\times m}$ is the parameter tensor
for lag $i$, and ${\cal C}\in\mathbb{R}^{\ell\times 1\times m}$ is a tensor of
centers. 2) We extend these results by adding the capability for modeling
seasonal and non-stationary tensor data by adding a differencing step that
results in an extension to the classical autoregressive integrated (ARI) model
and seasonal autoregressive (SAR) model in a tensor framework. We refer the
integration step (differencing) to this model as an ${\cal L}$-TARI model and
we refer to the seasonal differencing to this model as an ${\cal L}$-STAR.
Experimental results on benchmark datasets are presented to compare the
proposed approach against both the traditional MLDS models and a long-short
term memory artificial neural network (LSTM) in NIPS2013_5117 ; art:weijun2018
. The results suggest that in most multilinear forecasting problems, the
current approach outperforms previous methods in their ability to execute both
long-term and short-term forecasts.

The remainder of the paper is organized as follows. First we provide some
mathematical background for the tensor linear algebra in Section 2 . Next, we
provide some preliminary information regarding the ${\cal L}$-transform Tensor
AutoRegressive ($\mathcal{L}$-TAR) method and outline the different variants
of $\mathcal{L}$-TAR ($\mathcal{L}$-TARI, $\mathcal{L}$-STAR, and
$\mathcal{L}$-STARI) in Section 3 . Experimental results of our proposed
method are shown first with synthetic data and then with 4 standard benchmark
data sets in Section 4. In Section 5, we provide some interpretive remarks and
outline some directions of future work.

##  2 Mathematical Preliminaries

In order to keep this paper self contained, we will outline some of the
mathematical foundations of the tensor decompositions presented in KilMP08 ;
KilM09 ; Bra10 ; art:Hao13 ; art:Kilmer13 ; art:Hoover11 ; art:Hoover18 ;
liu2017fourthorder .

###  2.1 Mathematical Preliminaries

In the context of the current work, the term _tensor_ refers to a multi-
dimensional array of numbers, sometimes called an _n-way_ or _n-mode_ array.
For example, we say ${\cal A}$ is a third-order tensor if ${\cal
A}\in\mathbb{R}^{\ell\times m\times n}$ where _order_ is the number of ways or
modes of the tensor. Thus, matrices are second-order tensors and vectors are
first-order tensors.

First, we will introduce some basic notation and review the basic definitions
from KilMP08 ; KilM09 ; Bra10 ; art:Hao13 ; art:Kilmer13 . It will be
convenient to have an indexing on our tensor by breaking the tensor ${\cal
A}\in\mathbb{R}^{\ell\times m\times n}$ up into various slices and tubal
elements. We will denote ${\cal A}_{(i)}$ as the $i^{\text{th}}$ lateral slice
whereas will denoted ${\cal A}^{(j)}$ as the $j^{\text{th}}$ frontal slice .
In terms of Python slicing, this means ${\cal A}_{(i)}\equiv{\cal A}[:,:,i]$
while ${\cal A}^{(j)}\equiv{\cal A}[j,:,:]$. We will denote the
$i,k^{\text{th}}$ frontal tube in ${\cal A}$ as ${\bf a}_{ik}$; i.e., ${\bf
a}_{ik}={\cal A}[:,i,k]$, and we will denote the $i,k^{\text{th}}$ vertical
tube in ${\cal A}$ as $\mathbf{a}^{ik}$; i.e., $\mathbf{a}^{ik}={\cal
A}[i,:,k]$. Indeed, these tubes will play a role similar to scalars in
$\mathbb{R}$ so they will have special meaning for us in the present work.
Thus, we make the following definition:

######  Definition 1.

An element ${\bf e}\in\mathbb{R}^{1\times 1\times n}$ is called a tubal-scalar
of length $n$.

The t-product for multiplying tensors, developed by Kilmer et al. KilMP08 ;
KilM09 ; Bra10 ; art:Hao13 ; art:Kilmer13 , performs a product on two third-
order tensors which produces a third-order tensor. The resulting complex
arithmetic associated with the t-product was built around the discrete Fourier
transform (DFT) and an algebra of circulants. However, because of the complex
arithmetic, this becomes computationally prohibitive for large datasets.
Therefore, the research community found that two variations on the original
formulation that utilize either the discrete cosine transform (DCT), or the
discrete wavelet transform (DWT) KERNFELD2015545  gave alternative solutions.
We define the following operators by combining the notation outlined in
KERNFELD2015545  with the prior work outlined in KilMP08 ; KilM09 ; Bra10 ;
art:Hao13 ; art:Kilmer13 ; art:Hoover11 ; art:Hoover18 :

We anchor the MatVec command to the frontal slices of the tensor such that
$\mbox{\tt MatVec}({\cal A})$ takes an $\ell\times m\times n$ tensor and
returns a block $\ell n\times m$ matrix

$\mbox{\tt MatVec}({\cal A})=\left[\begin{array}[]{c}{\cal A}^{(1)}\\\ {\cal
A}^{(2)}\\\ \vdots\\\ {\cal A}^{(n)}\end{array}\right].$

We anchor the MatView command to the frontal slices of the tensor such that
$\mbox{\tt MatView}({\cal B})$ takes an $\ell\times m\times n$ tensor and
returns a block diagonal $\ell n\times mn$ matrix

$\mbox{\tt MatView}({\cal B})=\left[\begin{array}[]{cccc}{\cal
B}^{(1)}&\textbf{0}&\cdots&\textbf{0}\\\ \textbf{0}&{\cal
B}^{(2)}&\cdots&\textbf{0}\\\ \vdots&\vdots&\ddots&\vdots\\\
\textbf{0}&\textbf{0}&\cdots&{\cal B}^{(n)}\\\ \end{array}\right],$

where the 0’s in the previous matrix represent an $\ell\times n$ zero matrix.

The operation that takes both $\mbox{\tt MatVec}({\cal A})$ and/or $\mbox{\tt
MatView}({\cal B})$ back to tensor form ($\ell\times m\times n$) is the fold
command:

$\mbox{\tt fold}(\mbox{\tt MatVec}({\cal A}))={\cal A}\text{ and/or }\mbox{\tt
fold}(\mbox{\tt MatView}({\cal B}))={\cal B}.$

Finally, we anchor the Collect command to the collection of either mode-1 or
mode-2 tensors (i.e., vectors and matrices respectively) such that $\mbox{\tt
Collect}(\\{\boldmath{{\cal A}}_{1},\cdots,\boldmath{{\cal
A}}_{i},\cdots,\boldmath{{\cal A}}_{m}\\})$, $i=1,2,\dots,m$ and
$\boldmath{{\cal A}}_{i}\in\mathbb{R}^{k\times k}$ returns a tensor ${\cal
A}\in\mathbb{R}^{k\times k\times m}$ with the $\boldmath{{\cal A}}_{i}$ as its
frontal slices with increasing $i$ from front to back.

The above operators enable a generalized tensor product to be defined via any
invertible discrete transform ${\cal
L}:\mathbb{C}^{n}\rightarrow\mathbb{C}^{n}$. As such, we have the following
definition:

######  Definition 2.

The ${\cal L}$-transform of the tensor ${\cal A}$, given by

$\tilde{{\cal A}}={\cal L}({\cal A})\in\mathbb{C}^{\ell\times m\times n},$

is computed by applying the discrete transform of your choice along the tubes
${\bf a}_{ik}$ of ${\cal A}$. 222Note: the current work focuses on the DWT and
DCT. However, the DFT framework also applies here.

Using this formulation, given two third order tensors ${\cal
A}\in\mathbb{C}^{\ell\times m\times n}$ and ${\cal B}\in\mathbb{C}^{m\times
p\times n}$, we define the ${\cal L}$-product between 2 tensors as follows.

######  Definition 3.

The ${\cal L}$-product between ${\cal A}$ and ${\cal B}$ can be defined via
traditional convolution as

${\cal C}={\cal A}\bullet_{\cal L}{\cal B}={\cal L}^{-1}(\mbox{\tt
fold}(\mbox{\tt MatView}(\tilde{{\cal A}})\cdot\mbox{\tt MatVec}(\tilde{{\cal
B}}))),$

where we denote $\bullet_{\cal L}$ as the ${\cal L}$-product (henceforth we
will drop the subscript ${\cal L}$ in the ${\cal L}$-product). $\cdot$ is
computed via classical matrix multiplication, and the resulting tensor ${\cal
C}={\cal A}\bullet{\cal B}\in\mathbb{C}^{\ell\times m\times n}$.

##  3 Transform-Based Tensor Autoregression

In this section. we will discuss the details of building the proposed
extensions to the classical AR, ARI, and SARI models using the ${{\cal
L}\text{-transform}}$ and ${{\cal L}\textbf{-product}}$. Namely, we will show
how we can divide and conquer by recasting the multilinear time-series problem
into a subset of linear VAR problems by using the ${{\cal
L}\text{-transform}}$.

###  3.1 Model Overview

Our overarching goal is to construct the $p^{\text{th}}$ order tensor
autoregressive model (referred to as a ${\cal L}$-TAR($p$)) given by

${\cal Y}_{t}={\cal C}+{\cal A}_{1}\bullet{\cal Y}_{t-1}+\cdots+{\cal
A}_{p}\bullet{\cal Y}_{t-p}+{\cal E}_{t},$ (3)

by estimating the model parameters $\Theta=\\{{\cal A}_{1},{\cal
A}_{2},\dots,{\cal A}_{p},{\cal C}\\}$ from a collection of multilinear
observations ${\cal Y}_{j}\in\mathbb{R}^{\ell\times 1\times m}$,
$j=1,2,\dots,n$ with $n>>p$. $\bullet$ denotes the ${{\cal
L}\textbf{-product}}$ outlined in definition 3, ${\cal
A}_{i}\in\mathbb{R}^{\ell\times\ell\times m}$ is the model coefficient tensor
for lag $i$, ${\cal C}\in\mathbb{R}^{\ell\times 1\times m}$ is a tensor of
centers, and ${\cal E}_{t}$ represents the model errors. We assume the model
errors have zero mean, with constant variance, and are uncorrelated (i.e.,
E$\\{{\cal E}\\}=0$, E$\\{{\cal L}({\cal E},{\cal E}^{T})\\}=\Psi$, and
E$\\{{\cal E}_{i},{\cal E}_{j}\\}=0$ for $i\neq j$). A graphical illustration
of the ${\cal L}$-TAR($p$) model is shown in Fig. 1.

Figure 1: Graphical illustration of the proposed ${\cal L}$-TAR($p$) model.

###  3.2 Parameter Estimation

It is assumed the multilinear observations ${\cal Y}_{j}$ are correlated in
the sampling domain. However, the vertical tubes $\mathbf{y}^{ik}$ are
uncorrelated in the transform domain, therefore we proceed by computing
$\tilde{{\cal Y}}_{j}={\cal L}({\cal Y}_{j})$ for each $j=1,2,\dots,n$. As
such, we receive a collection of $m$ vector observations
$\mathbf{y}_{j}^{k}\in\mathbb{C}^{\ell\times 1\times 1}$ for $j=1,2,\dots,n$
and $k=1,2,\dots,m$ from the transformed multilinear observation $\tilde{{\cal
Y}}_{j}$. In other words, we sampled from $m$ different VAR processes in the
transform domain where each VAR process has a collection of $n$ multivariate
observations $\mathbf{y}_{j}^{k}$. As such, by applying the techniques of a
standard VAR process (least squares regression, maximum likelihood, or
expectation maximization), we estimate $m$ different VAR model parameters
$\mbox{\boldmath{$\theta$}}_{k}=\\{\tilde{{\cal A}}_{1}^{k},\tilde{{\cal
A}}_{2}^{k},\dots,\tilde{{\cal A}}_{p}^{k},\mathbf{c}^{i}\\}$, $k=1,2,\dots,m$
as outlined in Eqn. (2). This enables us to reconstruct the parameter tensors
$\\{\tilde{{\cal A}}_{1},\tilde{{\cal A}}_{2},\dots,\tilde{{\cal
A}}_{p},\tilde{{\cal C}}\\}$ by applying the Collect$(\cdot)$ command to each
of the parameter matrices/vectors in $\mbox{\boldmath{$\theta$}}_{k}$ for each
$k$. Finally, the inverse of the ${\cal L}$-transform is applied resulting in
the ${\cal L}$-TAR($p$) model parameters

${\cal L}^{-1}\\{\tilde{{\cal A}}_{1},\tilde{{\cal A}}_{2},\dots,\tilde{{\cal
A}}_{p},\tilde{{\cal C}}\\}\rightarrow\Theta=\\{{{\cal A}}_{1},{{\cal
A}}_{2},\dots,{{\cal A}}_{p},{{\cal C}}\\}.$

The entire process for constructing the ${\cal L}$-TAR($p$) model is
illustrated graphically in Fig. 2. The process is also shown in a algorithmic
fashion in Algorithm 1.

Figure 2: Graphical illustration of the overall process for computing the
${\cal L}$-TAR($p$) model. The top row illustrates how the original
observations ${\cal Y}_{j}$ are encoded using the ${\cal L}$-transform to
construct the individual $p$ observations $\mathbf{y}_{j}^{i}$ for estimating
the $m$ VAR models in the transform domain. The middle row illustrates how the
$m$ VAR model parameters $\mbox{\boldmath{$\theta$}}_{k}$ are estimated and
collected back to tensor form using the Collect$(\cdot)$ command. The bottom
row illustrates how the inverse ${\cal L}$-transform is applied to compute the
${\cal L}$-TAR($p$) model parameters $\Theta=\\{{\cal A}_{1},{\cal
A}_{2},\dots,{\cal A}_{p},{\cal C}\\}$. Algorithm 1 ${\cal L}$-TAR training
algorithm

1:procedure ${\cal L}$-TAR(${\cal Y}_{j}$, $p$)

2: Let $\tilde{{\cal Y}}_{j}$ be a new collection of tensors $\tilde{{\cal
Y}}_{j}\in\mathbb{C}^{\ell\times 1\times m}$ for $j=1,2,\dots,n$

3: for $j=1$ to $n$ do

4: $\tilde{{\cal Y}}_{j}={\cal L}({\cal Y}_{j})$

5: end for

6: Let $\mbox{\boldmath{$\theta$}}_{k}$ be a new collection of VAR model
parameters for $k=1,2,\dots,m$

7: for $k=1$ to $m$ do

8: $\mbox{\boldmath{$\theta$}}_{k}$ = VAR($\mathbf{y}_{j}^{i}$, $p$)

9: end for

10: $\\{\tilde{{\cal A}}_{1},\tilde{{\cal A}}_{2},\dots,\tilde{{\cal
A}}_{p},\tilde{{\cal C}}\\}=$ Collect$(\cdot)$

11: $\Theta={\cal L}^{-1}\\{\tilde{{\cal A}}_{1},\tilde{{\cal
A}}_{2},\dots,\tilde{{\cal A}}_{p},\tilde{{\cal C}}\\}$

12: return $\Theta$

13:end procedure

###  3.3 Complexity of Training

We will consider the time complexity of each step in the ${\cal L}$-TAR($p$)
training process for the multilinear observation ${\cal
Y}_{j}\in\mathbb{R}^{\ell\times 1\times m}$, $j=1,2,\dots,n$. When we compute
$\tilde{{\cal Y}}_{j}={\cal L}({\cal Y}_{j})$ for each $j=1,2,\dots,n$, the
time complexity of that observation is $O(n\ell m\log m)$. Performing the
transform of a $2^{\text{nd}}$-order tensor has time complexity of $O(\ell
m\log m)$ and that computation is performed for every multilinear observation
($n$ times).

To consider the complexity of training $m$ VAR models, we consider the time
complexity of each VAR model independently. We consider a VAR model described
in Eqn. (2) that is trained with ordinary least squares for vector
observations $\mathbf{y}_{j}\in\mathbb{C}^{\ell\times 1\times 1}$ for
$j=1,2,\dots,n$ from the transformed multilinear observation $\tilde{{\cal
Y}_{j}}$. How the $k^{\text{th}}$ VAR model is trained is structuring the
transformed observations into

$Y=XA+\epsilon,$ (4)

where we denote the block matrices
$Y\in\mathbb{R}^{(n-p)\times\ell},X\in\mathbb{R}^{(n-p)\times(\ell
p+1)},A\in\mathbb{R}^{(\ell p+1)\times\ell}$ as

$Y=\begin{pmatrix}y^{T}\\\ y_{p+1}^{T}\\\ \vdots\\\ y_{n}^{T}\end{pmatrix},$
$X=\begin{pmatrix}1&y_{p-1}^{T}&\cdots&y_{1}^{T}\\\
1&y_{p}^{T}&\cdots&y_{2}^{T}\\\ \vdots&\vdots&\vdots&\vdots\\\
1&y_{n-1}^{T}&\cdots&y_{n-p}^{T}\\\ \end{pmatrix},$
$A=\begin{pmatrix}\mathbf{c}^{T}\\\ A_{1}\\\ \vdots\\\ A_{p}\\\
\end{pmatrix}.$

Ordinary least squares is then used to compute and estimate of $A$ via

$\hat{A}=(X^{T}X)^{-1}X^{T}Y.$

The time complexity of computing $X^{T}X$ is $O((\ell p+1)(n-p)(\ell p+1))$,
since the complexity of multiplying a $p\times q$ matrix with a $q\times r$
matrix is done iteratively in $O(pqr)$. With the same reasoning, the time
complexity of computing $X^{T}Y$ is $O((\ell p+1)(n-p)\ell)$, the time
complexity of computing $(X^{T}X)^{-1}$ using LU factorization is $O((\ell
p+1)^{3})$, and finally to compute the final product $(X^{T}X)^{-1}X^{T}Y$ is
$O((\ell p+1)^{2}\ell)$ this gives an overall complexity, after reducing, of
$O(\ell^{2}p^{2}n+\ell^{3}p^{3})$ for training a single VAR model. Therefore,
the time complexity will be $O(m(\ell^{2}p^{2}n+\ell^{3}p^{3}))$ for training
$m$ VAR models.

Applying the Collect$(\cdot)$ command is simply restructuring the data, so the
time complexity of that computation is $O(\ell^{2}mp)$. And finally,
performing the inverse transformation ${\cal L}^{-1}\\{\tilde{{\cal
A}}_{1},\tilde{{\cal A}}_{2},\dots,\tilde{{\cal A}}_{p},\tilde{{\cal C}}\\}$
has a time complexity of $O(p\ell^{2}m\log m)$ since we are performing the
inverse on $p$ third order tensors.

Note that in general, $p<<n$, therefore, $p$ is insignificant in terms of
computational cost and can be omitted. Combining all computations together, we
arrive at the final complexity being

$O(n\ell m\log m+m\ell^{2}n+m\ell^{3}+\ell^{2}m\log m).$

Note that this complexity can be reduced even further, because in general, $n$
will be much larger than $\ell,m$. This will simply reduce to $O(n)$ because
$\ell,m$ will generally be less then $n$ in most cases. In the current work,
we include the specifics of $\ell,m$ in the complexity, however we “expect” to
see linear complexity in $n$.

###  3.4 Illustrative Example of the $\mathcal{L}$-TAR($p$) model

As an initial evaluation of the proposed approach, we construct a ground truth
$\mathcal{L}$-TAR($1$) model using,

${\cal Y}_{t}={\cal A}_{1}\bullet{\cal Y}_{t-1}+{\cal C}+{\cal E}_{t}.$ (5)

The parameters $\Theta=\\{{\cal A}_{1},{\cal C}\\}$ were arbitrarily selected
as

$\mathcal{A}_{1}^{(1)}=\mathcal{A}_{1}^{(3)}=\begin{pmatrix}-0.2&0&0\\\
0&-0.2&0\\\ 0&0&-0.2\end{pmatrix},$
$\mathcal{A}_{1}^{(2)}=\begin{pmatrix}0.2&0&0\\\ 0&0.2&0\\\
0&0&0.2\end{pmatrix},$

and

$\mathcal{C}=\begin{pmatrix}0.1&0.1&0.1\\\ 0.1&0.1&0.1\\\
0.1&0.1&0.1\end{pmatrix}.$

We generated $n=2000$ observations, i.e., ${\cal Y}_{j}$, $j=1,2,\dots,n$ with
${\cal Y}_{0}$ randomly initialized. ${\cal E}_{t}$ is white noise generated
under a uniform distribution between -1 and 1. Using the proposed
$\mathcal{L}$-TAR($1$) model outlined in Eqn. (5), our goal was to estimate
the model parameters $\hat{\Theta}=\\{\hat{{\cal A}}_{1},\hat{{\cal C}}\\}$
from the observations and compare to the ground truth parameters outlined
above. The resulting estimates are

$\hat{\mathcal{A}}_{1}^{(1)}=\begin{pmatrix}-0.185&-0.004&-0.000\\\
-0.006&-0.189&-0.014\\\ -0.000&-0.004&-0.188\end{pmatrix},$
$\hat{\mathcal{A}}_{1}^{(2)}=\begin{pmatrix}0.203&0.004&-0.005\\\
-0.007&0.197&-0.006\\\ -0.000&-0.011&0.193\end{pmatrix},$
$\hat{\mathcal{A}}_{1}^{(3)}=\begin{pmatrix}-0.203&-0.001&-0.005\\\
-0.007&-0.209&-0.001\\\ -0.001&0.007&-0.198\end{pmatrix},$

and

$\hat{\mathcal{C}}=\begin{pmatrix}0.090&0.100&0.106\\\ 0.101&0.094&0.110\\\
0.111&0.092&0.102\end{pmatrix}.$

While not exact, due to the addition of noise terms, the above example
illustrates the effectiveness and accuracy of the proposed approach. As can be
seen, the estimates of the resulting model parameters
$\hat{\Theta}=\\{\hat{{\cal A}}_{1},\hat{{\cal C}}\\}\approx{\Theta}=\\{{{\cal
A}}_{1},{{\cal C}}\\}$.

###  3.5 Considering both seasonality and non-stationarity

For the formulation of $\mathcal{L}$-TAR($p$), the multilinear observations
${\cal Y}_{j}$ must fulfill two conditions: 1) the observations ${\cal Y}_{j}$
are stationary and 2) there is no seasonal trend within the observations.
There are many real-world applications (video sequences for example are non-
stationary) where these assumptions are either invalid or violated. We can
overcome these restrictions by extending the traditional VAR($p$) techniques
to account for seasonality, non-stationarity, or both. We will illustrate how
such extensions can be applied to a multilinear framework in the following
subsections.

####  3.5.1 Non-stationary derivation ($\mathcal{L}$-TARI)

We enforce stationarity within a time series utilizing an integration step. A
time series is stationary if the observations ${\cal Y}_{j}$ have constant
mean and variance, i.e., E$\\{{\cal Y}\\}={\cal M}$ and E$\\{({\cal Y}-{\cal
M})^{2}\\}=\Psi$, where ${\cal M}\in\mathbb{R}^{\ell\times 1\times m}$ is the
mean tensor. Similar to how stationarity is enforced in the VAR process,
enforcing stationarity within the tensor time-series can be performed by
applying $d\in\mathbb{Z}^{+}$ lagged differences to our observations, ${\cal
Y}_{j}$. The resulting multilinear model is referred to as a
$\mathcal{L}$-TARI($p$, $d$) model, where $d$ is the order of differencing
(i.e., the amount of times that Eqn. (6) is applied to the observations ${\cal
Y}_{j}$). We apply the lagged difference $d$ times as

$\mathcal{Y}_{j}^{\prime}={\cal Y}_{j}-{\cal Y}_{j-1}.$ (6)

Then using the differenced observations ${\cal Y}_{j}^{\prime}$, the ${\cal
L}-$TAR model is constructed and the forecast is performed for
$w\in\mathbb{Z}^{+}$ steps. This results in a multilinear response $\hat{{\cal
Y}}_{k}^{\prime}$ for $k=n+1,n+2,\dots,n+w$. The differencing must be removed
from the response ${\cal Y}_{j}^{\prime}$ by inverting Eqn. (6) as,

$\begin{array}[]{ccl}{\cal Y}_{k}&=&\mathcal{Y}_{k}^{\prime}+{\cal Y}_{k-1}\\\
{\cal Y}_{k}&=&\mathcal{Y}_{k}^{\prime}+\mathcal{Y}_{k-1}^{\prime}+{\cal
Y}_{k-2}\\\ &\vdots&\\\ {\cal
Y}_{k}&=&\mathcal{Y}_{k}^{\prime}+\mathcal{Y}_{k-1}^{\prime}+\cdots+\mathcal{Y}_{n}.\end{array}$

####  3.5.2 Seasonality derivation ($\mathcal{L}$-STAR)

We can also enforce no seasonal trend within a time-series by utilizing an
integration step Hyndman . This is done by applying a seasonal difference to
our observations ${\cal Y}_{j}$. The resulting multilinear model is referred
to as ${\cal L}$-STAR($p$, $s$), where we consider $1<s<n$ as the period of
the seasonal trend. We apply the seasonal difference to our observations
${\cal Y}_{j}$ as

$\mathcal{Y}_{j}^{\prime}=\mathcal{Y}_{j}-\mathcal{Y}_{j-s}.$ (7)

As before, the model is constructed using the differenced observations ${\cal
Y}_{j}^{\prime}$, and the forecast is performed for $w$ steps to obtain the
multilinear response $\hat{\cal Y}_{k}^{\prime}$. The difference must be
removed to recover the response $\hat{\cal Y}_{k}$ by inverting Eqn. (7) as,

$\begin{array}[]{ccl}\mathcal{Y}_{k}&=&\mathcal{Y}_{k}^{\prime}+\mathcal{Y}_{k-s}\\\
\mathcal{Y}_{k}&=&\mathcal{Y}_{k}^{\prime}+\mathcal{Y}_{k-s}^{\prime}+\mathcal{Y}_{k-2s}\\\
&\vdots&\\\
\mathcal{Y}_{k}&=&\mathcal{Y}_{k}^{\prime}+\mathcal{Y}_{k-s}^{\prime}+\cdots+{\cal
Y}_{n-s+k}.\end{array}$

####  3.5.3 Combining both non-stationarity and seasonality
($\mathcal{L}$-STARI)

A combination of ${\cal L}$-TARI and ${\cal L}$-STAR can be done when
presented with non-stationary observations after applying a seasonal
difference, or vice-versa. This results in a multilinear model, referred as
${\cal L}$-STARI($p$, $d$, $s$), where we consider both the order of
difference $d$ and the period of seasonality $s$. We apply a number of
differences in the observations ${\cal Y}_{j}$ to form ${\cal
Y}_{j}^{\prime}$, construct the model based on the differenced observations
${\cal Y}_{j}^{\prime}$, then once the forecasted response $\hat{\cal
Y}_{k}^{\prime}$ is made, remove the differencing, as a similar fashion as the
two above methods. The order of which difference to apply is up to the user,
where we can apply the the lagged difference from Eqn. (6) then seasonal
difference first from Eqn. (7), or vice-versa.

##  4 Experimental Results

###  4.1 Qualitative Analysis

To validate the effectiveness of the proposed approach, we ran qualitative
evaluations on two data sets. First is the MNIST 333The MNIST dataset analysed
during the current study is available in the MNIST repository,
http://yann.lecun.com/exdb/mnist/ dataset where the goal is to forecast a
sequence of ordered image data (handwritten digits in this case) and the
second is a dynamic time-varying synthetic weighted graph 444The synthetic
graph dataset generated and analysed during this study are included in this
published article.. We will revisit both datasets in the quantitative section.
Both evaluations are presented in the following subsections.

####  4.1.1 MNIST

To evaluate the effectiveness of the proposed approach when forecasting image
data, we use $n=2000$ samples from the MNIST dataset web:MNIST18  which
contains 60,000 samples of handwritten digits ranging from 0 - 9. We sequence
the images in a repeating pattern from 0-9 throughout all 2000 samples, i.e.,
different observations were selected and ordered 0 - 9 in a repeating fashion,
a random sample of this sequence is shown in the top of Fig. 3. The collection
of this sequence can be represented as multilinear observations ${\cal
Y}_{j}\in\mathbb{R}^{28\times 1\times 28}$ for $j=1,2,\dots,n$. By
construction, we notice a few things about this particular dataset: 1) the
data is non-stationary due to the different representations of individual
digits and 2) the data is seasonal due to our particular sequencing (this was
intentional to illustrate seasonality within the data). The goal is to
estimate the model parameters $\Theta=\\{{\cal A}_{1},{\cal A}_{2},\dots,{\cal
A}_{p},{\cal C}\\}$ for the ${\cal L}$-STARI($p$, $d$, $s$) model from the
collection of multilinear observations ${\cal Y}_{j}$. When estimating the
model parameters, we set $s=10$ due to the number of different digits in the
sequence, $d=1$ for the lagged differences, and $p=10$ is found empirically.
Once the model parameters $\Theta$ are estimated from the multilinear
observations, the ${\cal L}$-STARI($p$, $d$, $s$) model is used to forecast
the next 10 images in the sequence. As can be seen in the bottom row of Fig.
3, the resulting 10 step forecast is qualitatively quite good.

Figure 3: MNIST qualitative evaluation of the proposed ${\cal L}$-STARI($p$,
$d$, $s$) model for $p=10$, $s=10$, and $d=1$. Top: a sampling 10 images from
$n=2000$ observations of the MNIST dataset ordered from 0-9 (repeating). Each
image was transformed into a multilinear observation ${\cal Y}_{j}$ and used
to estimate the model parameters of the ${\cal L}$-STARI(10, 1, 10) model.
Bottom: illustration of a 10 step forecast, i.e., ${\cal Y}_{t}$,
$t=1,2,\dots,10$ for the proposed model.

####  4.1.2 Synthetic Graph

To evaluate the effectiveness of the proposed approach when forecasting graph-
states and community separation in a dynamic graph, we use $n=2000$ samples
from a user generated synthetic weighted graph. The goal is to generated a
weighted, undirected graph with deterministic edges and time-varying community
separation. This graph contains 20 nodes, resulting in multilinear
observations ${\cal Y}_{j}\in\mathbb{R}^{20\times 1\times 20}$ using the
adjacency matrix representation of the graph. Deterministic edge weights are
generated by altering the edge weights between $[0,1]$ in a sinusoidal
fashion. We then apply a shift to each edge. The goal is to simulate community
separation by creating a repeating sinusoidal pattern of the graph starting
with one large community (20-nodes), separating into two smaller communities
(10=nodes each), and combining back into the original large community.
Mathematically, the collection of graph observations (i.e., graph adjacency
matrices $Y_{\ell,m}$) are constructed using

$Y_{\ell,m}=\begin{cases}0&\text{if }\ell=m\\\
\frac{1+\text{sin}(xp_{1}+S_{\ell,m})}{2}+\epsilon&Y_{\ell,m}\in\text{block
diagonal}\\\
\frac{1+\text{cos}(xp_{2})}{2}*\frac{1+\text{sin}(xp_{1}+S_{\ell,m})}{2}+\epsilon&Y_{\ell,m}\notin\text{block
diagonal}\end{cases},$

where $p_{1}$ is the period of the edges and $p_{2}$ is the period of the
community separation where in general, $p_{2}>>p_{1}$. $\epsilon$ is random
white noise with $\epsilon\sim N(0,\sigma=0.02)$. Graphically, a subset of
adjacency matrices for selected time-instances are illustrated in the top row
of Fig. 4 with the corresponding graphs illustrated in the top row of Fig. 5.

Figure 4: Adjacency matrix representation of the graph qualitative evaluation
using ${\cal L}$-STAR($p$, $s$) for $p=40$ and $s=200$. generation for
evaluation of the proposed models. $S_{\ell,m}$ is the shift designated for
that edge. Top: a sampling of 5 adjacency matrices from $n=2000$ observations
of the synthetic graph. Each image was transformed into a multilinear
observation ${\cal Y}_{j}$ and used to estimate the model parameters of the
${\cal L}$-STARI(40, 200) model. Bottom: sampling of the multi-step forecast
for the adjacency matrix from the ${\cal L}$-STARI(40, 200) model. Figure 5:
Graph representation of the graph qualitative evaluation using ${\cal
L}$-STAR($p$, $s$) for $p=40$ and $s=200$. generation for evaluation of the
proposed models. Top: a sampling of 5 graphs from $n=2000$ observations of the
synthetic graph where the edge weight is represented by the opacity. Each
image was transformed into a multilinear observation ${\cal Y}_{j}$ and used
to estimate the model parameters of the ${\cal L}$-STARI(40, 200) model.
Bottom: sampling of the multi-step forecast for the graph from the ${\cal
L}$-STARI(40, 200) model.

The goal is to estimate model parameters $\Theta=\\{{\cal A}_{1},{\cal
A}_{2},\dots,{\cal A}_{p},{\cal C}\\}$ for the ${\cal L}$-STAR($p$, $s$) from
the observations $Y_{\ell,m}$ where each $Y_{\ell,m}$ is treated as a
$\ell\times 1\times m$ tensor. When estimating the model parameters, as
expected, the period that gives the best fit is the period for the community
separation, with $s=200$ and $p=40$ lags. Once the model parameters $\Theta$
are estimated, we forecast the next set of graph states, the results of which
are illustrated in the bottom row of Fig. 4 (adjacency matrix) & 5 (graph
state).

(a) SST

(b) Video

(c) NASDAQ-100

(d) Graph

Figure 6: Results of single-step forecasting for the quantitative evaluation.
Error is the absolute error, i.e., $\|{\cal Y}_{t}-\hat{{\cal Y}}_{t}\|_{F}$.

(a) SST

(b) Video

(c) NASDAQ-100

(d) Graph

Figure 7: Results of multi-step forecasting for the quantitative evaluation.
Error is the absolute error, i.e., $\|{\cal Y}_{t}-\hat{{\cal Y}}_{t}\|_{F}$.
Notice that ${\cal L}$-STAR is not trained for the NASDAQ-100 dataset due to
there being no seasonality.

###  4.2 Quantitative Evaluation and Experimental Results

As a quantitative evaluation, we compare our proposed approach to current
state-of-the-art in multilinear time-series methods. Namely, the proposed
approach is compared against the $\mathcal{L}$-MLDS model proposed in
art:weijun2018  and a convolutional Long Short-Term Memory Neural Network
model. In an effort to compare and contrast both methods, we use a subset of
the same datasets proposed in the $\mathcal{L}$-MLDS model in art:weijun2018
555The SST, NASDAQ-100, and Video datasets generated during and analysed
during the current study are available in the L-MLDS-for-Tensor-Time-Series
repository, https://github.com/XiaoYangLiu-FinRL/L-MLDS-for-Tensor-Time-Series
as well as the synthetic graph generated in section 4.1. The information
pertaining to each dataset is outlined in Table 1 (additional details on the
datasets can be found in art:weijun2018 ), and a tabulated list of all models
used in our evaluation are outlined in Table 2, with the details of each
provided in the following subsections At the end, we will revist the MNIST
dataset outlined in section 4.1.1 for a separate quantitative evaluation.

Table 1: Datasets Used in the Qualitative Evaluation of the Proposed Multilinear Forecasting Approach. Dataset | Notes
---|---
SST | A $5\times 6$ grid of sea-surface temperatures. The first
| 1800 hours are used for training and the
| last 200 hours are used for testing. art:weijun2018
Video | A $10\times 10$ video of the ocean. The first
| 1000 hours are used for training and the
| last 171 hours are used for testing. art:weijun2018
NASDAQ-100 | Opening, closing, high, and low for 50
| randomly-chosen NASDAQ-100 companies ($50\times 4$).
| The first 2000 days are used for training and the
| last 186 days are used for testing. art:weijun2018
Synthetic Graph | Graph synthetically created with deterministic edge
| and communities as described in section 4.1
| ($20\times 20$). The first 1800 time
| slices are used for training and the
| last 200 time slices are used for testing.
Table 2: Different Variations of the Models Used in our Quantitative Evaluation Model | Notes
---|---
$\mathcal{L}$-TAR |  ${\cal L}$-transform computed using the DWT and DCT
$\mathcal{L}$-TARI |  $\mathcal{L}$-TAR for non-stationary data
$\mathcal{L}$-STAR |  $\mathcal{L}$-TAR for seasonal data
LSTM | A Long Short-Term Memory Neural Network
$\mathcal{L}$-MLDS | Outlined in art:weijun2018  using DWT and DCT

Two different evaluations are performed on both real and synthetic datasets:
1) single-step forecasting, where we estimate the multilinear response
$\hat{\cal Y}_{t}$ for $t=p+1,p+2,\dots,w$ using ground truth observations
${\cal Y}_{j}$ for $j=1,2,\dots,p$. The assumption here is that we’re only
interested in forecasting the next time-step using observed historical data
and 2) multi-step forecasting, where we estimate the multilinear response
$\hat{{\cal Y}}_{t}$ for $t=p+1,p+2,\dots,w$ using estimated observations
$\hat{{\cal Y}}_{j}$ for $j=1,2,\dots,p$. In some situations we are interested
in longer term forecasting (e.g., weather prediction). However, in general,
the single-step solution will be much more accurate because the forecast is
using the true observations of ${\cal Y}_{j}$ as opposed to our forecasted
estimates. We report the result of both of these evaluations in Fig. 6 and 7.
These figures illustrates the absolute error in the forecast, i.e., $\|{\cal
Y}_{t}-\hat{{\cal Y}}_{t}\|_{F}$.

Because both $\mathcal{L}$-TARI and $\mathcal{L}$-STARI models are more suited
for multi-step forecasting, we evaluate the original $\mathcal{L}$-MLDS and
$\mathcal{L}$-TAR for the single-step forecasting evaluation. To illustrate
the model’s ability to make long-term predictions, we present multi-step
forecasting evaluations for all models presented in Table 2. In art:weijun2018
$\mathcal{L}$-MLDS only evaluates single-step forecasting, therefore, we
modify their proposed method to make it more suitable for a multi-step
forecasting evaluation. The results of the evaluations for the datasets
outlined in Table 1. The details of each experiment will be outlined in the
following subsections.

####  4.2.1 SST

The SST dataset is a $5\times 6$ grid of sea-surface temperatures, where the
observations were recorded every hour art:weijun2018 . Each observation can be
represented as a multilinear observation ${\cal Y}_{j}\in\mathbb{R}^{5\times
1\times 6}$ for $n=2000$. The first 1800 hours are used to construct the
proposed models and the last 200 hours are used for evaluation. The model
configurations for this dataset can be seen in Table 3. Fig. 6(a) illustrates
the single-step forecasts and the Fig. 7(a) illustrates multi-step
forecasting. As seen in the figure, for single-step forecasting, ${\cal
L}$-TAR outperforms all other methods and is nearly identical to the ground
truth data. For multi-step forecasting, ${\cal L}$-TAR is comparable with the
${\cal L}$-STAR model. Both methods however, outperform the other multilinear
forecasting methods.

Table 3: Model configuration for SST dataset Model | Single-step | Multi-Step
---|---|---
$\mathcal{L}$-TAR | $p=5$ | $p=19$
$\mathcal{L}$-TARI | NA |  $p=19$ & $d=1$
$\mathcal{L}$-STAR | NA |  $p=3$ & $s=24$
LSTM | NA | 2 LSTM layers & relu activation

####  4.2.2 Video

The video dataset is a $10\times 10$ gray-scale video of the ocean, where the
observations were recorded every frame art:weijun2018 . Each observation can
be represented as a multilinear observation ${\cal
Y}_{j}\in\mathbb{R}^{10\times 1\times 10}$ for $n=1171$. The first 1000 frames
are used to construct the models and the last 171 frames are used for testing.
The model configurations for this dataset can be seen in Table 4. Fig. 6(b)
illustrates the single-step forecasts and Fig. 7(b) illustrates multi-step
forecasting. For single-step forecasting, ${\cal L}$-TAR and ${\cal L}$-MLDS
have the same performance. For multi-step forecasting, ${\cal L}$-TARI
performs the best until around the $25^{\text{th}}$ frame, then the LSTM
performs the best afterwards.

Table 4: Model configuration for Video dataset Model | Single-step | Multi-Step
---|---|---
$\mathcal{L}$-TAR | $p=10$ | $p=13$
$\mathcal{L}$-TARI | NA |  $p=9$ & $d=1$
$\mathcal{L}$-STAR | NA |  $p=9$ & $s=10$
LSTM | NA | 2 LSTM layers & sigmoid activation

####  4.2.3 NASDAQ-100

The NASDAQ-100 dataset contains the opening, closing, high, and low stock
price of the day for 50 random NASDAQ-100 companies, resulting in a $50\times
4$ grid art:weijun2018 . Each observation can be represented as a multilinear
observation ${\cal Y}_{j}\in\mathbb{R}^{50\times 1\times 4}$ for $n=2186$.
${\cal L}$-STAR was not trained since there was no seasonality. The model
configurations for this dataset can be seen in Table 5. Fig. 6(c) illustrates
the single-step forecasts and Fig. 7(c) illustrates multi-step forecasting.
For single-step forecasting, ${\cal L}$-TAR and ${\cal L}$-MLDS show equal
performance. For multi-step forecasting, ${\cal L}$-LTARI outperforms the
other methods until the $100^{\text{th}}$ day where the LSTM begins to
outperform all methods.

Table 5: Model configuration for NASDAQ-100 dataset Model | Single-step | Multi-Step
---|---|---
$\mathcal{L}$-TAR | $p=10$ | $p=5$
$\mathcal{L}$-TARI | NA |  $p=16$ & $d=1$
LSTM | NA | 2 LSTM layers & relu activation

####  4.2.4 Synthetic Graph

The synthetic graph we used is the same generated in the quantitative
evaluation section 4.1, which results in a $20\times 20$ adjacency matrix.
Each observation can be represented as a multilinear observation ${\cal
Y}_{j}\in\mathbb{R}^{20\times 1\times 20}$ for $n=2000$. The model
configurations for this dataset can be seen in Table 6. Fig. 6(d) illustrates
the single-step forecasts and Fig. 7(d) illustrates multi-step forecasting.
For single-step forecasting, ${\cal L}$-MLDS performs slightly better
throughout. For multi-step forecasting, ${\cal L}$-STAR, ${\cal L}$-TAR and
${\cal L}$-TARI are comparable throughout.

Table 6: Model configuration for Synthetic Graph dataset Model | Single-step | Multi-Step
---|---|---
$\mathcal{L}$-TAR | $p=40$ | $p=200$
$\mathcal{L}$-TARI | NA |  $p=200$ & $d=1$
$\mathcal{L}$-STAR | NA |  $p=40$ & $s=200$
LSTM | NA | 2 LSTM layers & sigmoid activation

####  4.2.5 MNIST

As a quantitative evaluation of the MNIST dataset, outlined in section 4.1.1,
we compare the distance to the forecasted digit with the (correct) ground
truth digit and (incorrect) every other digit. We compute the distance via
absolute error, which is normalized by the amount of pixels, i.e.,
$\frac{\|{\cal Y}_{t}-\hat{{\cal Y}}_{t}\|_{F}}{28^{2}}$. We use the proposed
${\cal L}$-STARI($p$, $d$, $s$) model for $p=10$, $s=10$, and $d=1$. We also
trained with $n=2000$ observations and tested with 200 observations. The
result of this evaluation can be seen in Fig. 8. We can see that, similar to
the qualitative results presented in Fig. 3, quantitatively, the forecasts are
very close.

Figure 8: Qualitative evaluation of the MNIST dataset outlined in section
4.1.1. We use the proposed ${\cal L}$-STARI($p$, $d$, $s$) model for $p=10$,
$s=10$, and $d=1$. Trained with $n=2000$ observations and tested with 200
observations. We compare the result quantitatively using the absolute error,
which is normalized with the amount of pixels, i.e., $\frac{\|{\cal
Y}_{t}-\hat{{\cal Y}}_{t}\|_{F}}{28^{2}}$. We compare the distance to the
forecast with the true digit and the distance to the forecast with every other
digit.

###  4.3 Speed Evaluation

####  4.3.1 Execution Time between Models

As an evaluation of speed, we compared the execution time for all multi-step
experiments performed in the previous section. Each model is trained 20 times
and the average execution time is recorded in Table 7. We can see that for
this implementation, our proposed model has significant speedup. All code was
implemented in Python.

Table 7: Execution time (in seconds) in multi-step quantitative analysis Model | SST | Video | NASDAQ-100 | Synthetic Graph
---|---|---|---|---
$\mathcal{L}$-TAR | 0.357 | 0.979 | 0.776 | 382.297
$\mathcal{L}$-TARI | 0.285 | 0.349 | 5.089 | 342.313
$\mathcal{L}$-STAR | 0.106 | 0.362 | NA | 80.97
LSTM | 97.033 | 95.415 | 267.246 | 927.801
$\mathcal{L}$-MLDS | 36.665 | 42.375 | 449.817 | 698.258

####  4.3.2 Execution Time for Parallelization

As another evaluation for speed, we considered how much speedup could be
achieved if multiprocessing were performed via distributed computing. Because
estimating the model parameters can be divided up into multiple VAR sub-
problems, we estimate these parameters by computing each sub-problem in
parallel. In order to compare the speedup ${\cal L}$-TAR achieved with multi-
processing, we used the generated synthetic graph dataset outlined earlier in
section 4.2.4. We used this dataset because we were able to scale the size of
the multilinear observations ${\cal Y}_{j}$ by simply selecting the number of
nodes $n$ in the graph. The number of nodes for each graph were incremented by
5, with the exception of nodes 45 to 48 as the maximum number of CPU cores on
our system was 48. Ultimately, this increased the number of VAR models
trained. In a similar fashion to the previous time test, each model was
trained 100 times. The speedup ($\text{Speedup}=\frac{\text{Sequential
time}}{\text{Parallel time}}$) for each test was recorded and is displayed in
Fig. 9. Inspection of Fig. 9 shows that speedup occurs almost linearly with
the number of nodes.

Figure 9: Swarm plot of speedup for multi-processing results where
$\text{speedup}=\frac{\text{Sequential time}}{\text{Parallel time}}$. 100
trials were ran for each number of nodes.

Referring to the time complexity (outlined earlier in section 3.3), we see
that the experimental results confirm our time complexity calculations. First,
we need to consider the time complexity achieved by training the VAR models in
parallel. For, we simply drop the $m$ term for our VAR training portion, so
the complexity of training a VAR model in parallel is now:

$O(n\ell m\log(m)+\ell^{2}n+\ell^{3}+\ell^{2}m\log(m)).$

Since we are using the adjacency matrix of a graph, $n$ is a constant and
$\ell=m$. Therefore, the complexity for this experiment is $O(m^{4})$ for
single-processing and $O(m^{3}\log(m))$ for multiprocessing. Thus, the overall
speedup is,

$\text{Speedup}=\frac{m^{4}}{m^{3}\log(m)}=\frac{m}{log(m)}.$

Fig. 9 verifies this result.

##  5 Conclusions and Future Directions

From both the qualitative results and the quantitative results presented in
our experiments, ${\cal L}$-TAR($p$), ${\cal L}$-TARI($p$, $d$), and ${\cal
L}$-STAR($p$, $s$) have been shown to be excellent methods for forecasting a
multilinear time series. In our experiments, our methods provided extremely
competitive forecasts and in most situations they outperformed the current
state of the art. Furthermore, our methods were shown to require less training
time than the other forecasting methods.

Future work includes applying extensions to ${\cal L}$-TAR in a similar
fashion to its autoregressive predecessors, such as applying moving averages
(${\cal L}$-TARMA, ${\cal L}$-STARMA, ${\cal L}$-TARIMA, ${\cal L}$-STARIMA)
and considering non-linearity with exogenous observations (${\cal L}$-NTARX).
Also, in the section 4, $p$, $d$, and $s$ was picked via trial and error.
Future work will also include creating similar tensor versions of auto-
correlation factor (ACF) and partial auto-correlation factor (PACF) plots to
have a more precise method of estimating these parameters.

## References

  *   * (1) Cates, J., Hoover, R.C., Caudle, K.: Transform-based tensor auto regression for multilinear time series forecasting. In: 2021 IEEE International Conference on Machine Learning and Applications (ICMLA) (2021)
  * (2) De Gooijer, J.G., Hyndman, R.J.: 25 years of time series forecasting. International journal of forecasting 22(3), 443–473 (2006)
  * (3) Box, P.J., GM, R.: Time Series Analysis: Forecasting & Control, (2008)
  * (4) Box., G.E.P.: Understanding exponential smoothin-a simple way to forecast sales and inventory. Quality Engineering 4(3), 561–566 (1991)
  * (5) Brown, R.G.: Statistical Forecasting for Inventory Control. McGraw-Hill, New York (1959)
  * (6) Brown, R.G.: Smoothing, Forecasting, and Prediction. Pretice Hall, Englewood Cliffs, NJ (1963)
  * (7) Frey, M., Caudle, K.: Introducing flow field forecasting. In: 2011 10th International Conference on Machine Learning and Applications and Workshops, vol. 1, pp. 395–400 (2011). IEEE
  * (8) Caudle, K.A.: Flowfield: Forecasts Future Values of a Univariate Time Series. (2014). R package version 1.0. https://CRAN.R-project.org/package=flowfield
  * (9) Caudle, K.A., Fleming, P.S.: Closest history flow field forecasting for ieee csci-iscs. In: 2016 International Conference on Computational Science and Computational Intelligence (CSCI), pp. 1202–1207 (2016). IEEE
  * (10) Caudle, K., Fleming, P., Pyeatt, L., Hoover, R.C.: Flow field forecasting with many predictors. In: Proceedings of the 2019 3rd International Conference on Compute and Data Analysis, pp. 109–115 (2019). ACM
  * (11) Caudle, K., Fleming, P., Hoover, R.C.: A review of flow field forecasting: A high-dimensional forecasting procedure. Wiley Interdisciplinary Reviews: Computational Statistics 13 (2021)
  * (12) Haykin, S.: Neural Networks and Learning Machines (3rd Ed.). Pearson, New York (2009)
  * (13) Hill, T., Marquez, L., O’Connor, M., Remus, W.: Artificial neural network models for forecasting and decision making. International Journal of Forecasting 10, 5–15 (1994)
  * (14) Asteriou D., S.G. & Hall: Applied Econometrics: Macmillan International Higher Education. SG Hall, London (2015)
  * (15) Rogers, M., Li, L., Russell, S.J.: Multilinear dynamical systems for tensor time series. In: in Neural Information Processing Systems (NIPS), pp. 2634–2642 (2013)
  * (16) Tucker, L.R.: Some mathematical notes on three-mode factor analysis. Psychometrika 31(3), 279–311 (1966)
  * (17) Lathauwer, L.D., Moor, B.D., Vandewalle, J.: A multilinear singular value decomposition. SIAM J. Matrix Anal. Appl. 21(4), 1253–1278 (2000)
  * (18) Lu, W., Liu, X.-Y., Wu, Q., Sun, Y., Elwalid, A.: Transform-Based Multilinear Dynamical System for Tensor Time Series Analysis. In: Neural Information Processing (NIPS) Workshop on Spatiotemporal Data, (2018)
  * (19) Kilmer, M.E., Martin, C.D., Perrone, L.: A third-order generalization of the matrix SVD as a product of third-order tensors. Technical Report TR-2008-4, Tufts University, Department of Computer Science (October 2008)
  * (20) Kilmer, M.E., Moravitz Martin, C.D.: Factorization strategies for third-order tensors. Linear Algebra and Its Applications (Special Issue in Honer of G.W.Stewart’s 75th birthday) (2009)
  * (21) Braman, K.: Third-order tensors as linear operators on a space of matrices. Linear Algebra and its Applications 433(7), 1241–1253 (2010)
  * (22) Hoover, R.C., Braman, K.S., Hao, N.: Pose estimation from a single image using tensor decomposition and an algebra of circulants. In: Int. Conf. on Intel. Robots and Sys. (2011)
  * (23) Hoover, R.C., Caudle, K., Braman, K.: In: 2018 17th IEEE International Conference on Machine Learning and Applications (ICMLA), pp. 578–584 (2018)
  * (24) Ozdemir, C., Hoover, R.C., Caudle, K.: 2DTPCA: A new framework for multilinear principal component analysis. In: 2021 IEEE International Conference on Image Processing (ICIP), pp. 344–348. IEEE, ??? (2021)
  * (25) Ozdemir, C., Hoover, R.C., Caudle, K.: Fast tensor singular value decomposition using the low-resolution features of tensors. In: 2021 20th IEEE International Conference on Machine Learning and Applications (ICMLA), pp. 527–533 (2021). IEEE
  * (26) Kilmer, M.E., Braman, K.S., Hao, N., Hoover, R.C.: Third order tensors as operators on matrices: A theoretical and computational framework with applications in imaging. SIAM Journal on Matrix Analysis and Applications (SIMAX) 34(1), 148–172 (2013)
  * (27) Hao, N., Kilmer, M.E., Braman, K.S., Hoover, R.C.: New tensor decompositions with applications in facial recognition. SIAM Journal on Imaging Science (SIIMS) 6(1), 437–463 (2013)
  * (28) Liu, X.-Y., Wang, X.: Fourth-order Tensors with Multidimensional Discrete Transforms (2017)
  * (29) Kernfeld, E., Kilmer, M., Aeron, S.: Tensor–tensor products with invertible linear transforms. Linear Algebra and its Applications 485, 545–570 (2015)
  * (30) Hyndman R. J., G. & Athanasopoulous: Forecasting: Principles and Practice. OTexts, Melbourne (2018)
  * (31) MNIST Database of Handwritten Digits. http://yann.lecun.com/exdb/mnist/