Because of the popularity of this area of research, we want to clarify the objectives and contributions of this paper from the outset. Our contributions are the following:

This paper is organized as follows. In Section 2, we give a brief outline of related work combining tensors and deep learning. In Section 3, we give the background notation on tensor-tensor products. In Section 4, we formally introduce tensor neural networks (t-NNs) and tubal loss functions. In Section 5, we extend t-NNs to stable architectures and outline a Hamiltonian-inspired architecture. In Section 6, we provide numerical support for using t-NNs over comparable traditional fully-connected neural networks. In Section 7, we discuss future work including implementations for higher-order data and new t-NN designs.

Let A∈ℝm1×m2×n be a real-valued, third-order tensor. Fixing the third-dimension, frontal slices A(k)∈ℝm1×m2 are matrices for k = 1, …, n. Fixing the second-dimension, lateral slices A→j∈ℝm1×1×n are matrices oriented along the third dimension for j = 1, …, m₂. Fixing the first and second dimensions, tubes aij∈ℝ1×1×n are vectors oriented along the third dimension for i = 1, …, m₁ and j = 1, …, m₂. We depict these partitions in Figure 1. While this paper focuses on real-valued, third-order tensors (three indexes), we note all of the presented concepts generalize to higher-order and complex-valued tensors.

We interpret tensors as t-linear operators (Kilmer and Martin, 2011; Kernfeld et al., 2015). Through our operator lens, it is possible to define analogous matrix algebraic properties for tensors, such as orthogonality and rank. Thus, this framework has been described as matrix-mimetic. We describe the fundamental tools to understand how tensors operate for this paper, and refer the reader to Kilmer et al. (2013, 2021) and Kernfeld et al. (2015) for details about the underlying algebra.

We define a product to apply matrices along the third dimension of a tensor (i.e., along the tubes).

Definition 3.1 (mode-3 product). Given A∈ℝm1×m2×n and M ∈ ℝ^{ℓ × n}, the mode-3 product, denoted A^≡A×3M, outputs an m₁×m₂×ℓ tensor with entries

for i₁ = 1, …, m₁, i₂ = 1, …, m₂, and k = 1, …, ℓ.

The mode-3 product can be generalized along any mode; see Kolda and Bader (2009) for details.

Next, we define the facewise product to multiply the frontal slices of two third-order tensors in parallel.

Definition 3.2 (facewise product). Given A∈ℝm1×ℓ×n and B∈ℝℓ×m2×n, the facewise product, denoted C≡A▵B, returns an m₁×m₂×n tensor where

for k = 1, …, n.

Combining Definition 3.1 and Definition 3.2, we define our tensor operation, the ⋆_M-product, as follows:

Definition 3.3 (⋆_M-product). Given A∈ℝm1×ℓ×n, B∈ℝℓ×m2×n, and an invertible n×n matrix M, the ⋆_M-product outputs an m₁×m₂×n tensor of the following form:

where X^≡X×3M.

We say that A and B live in the spatial domain and A^ and B^ live in the transform domain. We perform the facewise product in the transform domain, then return to the spatial domain by applying M⁻¹ along the tubes. If M is the identity matrix, the ⋆_M-product is exactly facewise product. If M were the discrete Fourier transformation matrix (DFT), we obtain the t-product (Kilmer and Martin, 2011). In this case, the frontal slices of A^ correspond to different frequencies in the Fourier domain and are therefore decoupled.

We can interpret the ⋆_M-product as a block-structured matrix product via

where ⊗ is the Kronecker product (Petersen and Pedersen, 2012). The block matrix structure, struct(A), depends on the choice of transformation, M. We consider two familiar examples: the facewise product (M = I_n) and the t-product (M = F_n, the DFT matrix):

While we never explicitly form Equation (2), the block structure will be helpful for subsequent analysis.

The ⋆_M-product yields a well-defined algebraic structure. Specifically, suppose we have tubes a, b∈ℝ^{1 × 1 × n}. Then,

where vec:ℝ1×1×n→ℝn×1 turns a tube into a column vector. The ⋆_M-product of tubes (3) is equivalent to post-multiplying by the structured matrix, R[b]. Note that R[·] implicitly depends on the choice of M, but we omit explicitly writing this dependence for notational simplicity. The tubes form a matrix subalgebra which dictates the algebraic structure imposed on the high-dimensional space (Kernfeld et al., 2015). As a result, the ⋆_M-product is matrix-mimetic and yields several familiar concepts.

Definition 3.4 (⋆_M-identity tube). The identity tube e ∈ ℝ^{1 × 1 × n} under the ⋆_M-product is

where 1 is the 1 × 1 × n tube containing all ones.

This gives rise to the notion of an identity tensor.

Definition 3.5 (⋆_M-identity tensor). A tensor I∈ℝm×m×n is the identity tensor if I(i,i,:)=e for i = 1, …, m where e is the ⋆_M-identity tube.

Note that if the size of third dimension is equal to one (i.e., n = 1), then Definition 3.5 collapses into the identity matrix. This is a hallmark of our ⋆_M-framework and matrix-mimeticity; the product and definitions reduce to the equivalent matrix definitions when the third dimension is removed.

Definition 3.6 (⋆_M-transpose). Given A∈ℝm1×m2×n, its transpose B≡A⊤∈ℝm2×m1×n is

for k = 1, …, n.

Note that if our transformation M is complex-valued, the transpose operator in the transform domain performs the conjugate transpose. However, because we are working with real-valued tensors in the spatial domain, the transpose will be real-valued as well.

When designing a neural network, we seek to balance a simple parameter space with an expressive feature space. However, traditionally fully-connected layers like (5) use parameters in a highly inefficient manner. We propose new tensor fully-connected layers for a more efficient parametrization while still creating a rich feature space. Specifically, we consider the following tensor forward propagation scheme:

for j = 1, …, d. Suppose our input features Y→0=Y→∈Y is of size m₀×1 × n. Here, the weight tensor Wj is of size m_j+1×m_j×n and the bias B→j is of size m_j+1×1 × n. The forward propagation through Equation (6) results in a tensor neural network F_tnn(·, θ) where θ ≡ (Wj,B→j)j=1d.

Through the illustration in Figure 2, we depict the number of weight parameters required to preserve the size of our feature space using either a dense matrix (Figure 2A) or a dense tensor under the ⋆_M-product (Figure 2B). Using the ⋆_M-algebra, we can reduce the number of weight parameters by a factor of n while maintaining the same number of features (i.e., maintaining a rich feature space). Beyond the parametric advantages, the multilinear ⋆_M-product incorporates the structure of the data into the features. This enables our t-NNs to extract more meaningful features, and hence improve the richness of our feature space.

Training a (tensor) neural network is posed as a stochastic optimization problem given by

where L:ℝmd+1×1×n×C→ℝ is the loss function that measures the misfit between the network prediction and the true target. The expectation is taken over all input-target pairs (Y→,c)∈D. The additional function R:Θ→ℝ regularizes the weights to promote desirable properties (e.g., smoothness), weighted by a regularization parameter λ > 0.

The loss function is chosen based on the given task. For regression, we often use mean squared error, and for classification, which is the focus of this paper, we often use cross entropy. Cross entropy loss, related to the Kullback–Leibler (KL) divergence (Kullback and Leibler, 1951), measures the distance between two probability distributions. In practice, we first transform the network outputs into a set of probabilities using exponential normalization. Specifically, we use the softmax function h:ℝp→Δp, defined entrywise as

where Δ^p is the p-dimensional unit simplex.

To preserve the algebraic integrity of t-NNs, we introduce a tubal variant of the softmax function. Drawing inspiration from Lund (2020), a tubal function, we start by defining tubal functions generally.

Definition 4.1 (tubal function). Given b ∈ ℝ^{1 × 1 × n}, a tubal function f:ℝ1×1×n→ℝ1×1×n acts on the action of b under the ⋆_M-product; that is,

In practice, a tubal function is applied pointwise in the transform domain.

We note that the pointwise function in Definition 4.1 is equivalent to applying a matrix function to the eigenvalues of the matrix R[b]. We provide a visualization of the effects of tubal functions compared to applying an entry-wise function in Example 4.2.

Example 4.2 (Visualizations of tubal functions). Consider the following RGB image B∈ℝ150×169×3 of a Tufts Community Apeal elephant (TCA, 2023), where 3 is the number of color channels (Figure 3A). We rescale each entry of B between 0 to 1 and consider the following transformation matrix:

To illustrate the effect of using tubal functions, we compare applying various functions as pointwise functions (i.e., independent of M) and as tubal functions using Equation (9) in Figure 3B.

Tubal functions are able to capture shared patterns among the tubes that entrywise operators ignore. Each choice of tubal function highlights different features of the image, such as the body with f(x) = max(x, 0) and the heart with f(x) = tanh(x). In comparison, the entrywise counterparts do not bring new insights or structure. A striking difference occurs for f(x) = sign(x). The entrywise operator turns all pixels white because the RGB entries are all nonnegative. In contrast, the tubal equivalent is applied in the transform domain, and hence reveals meaningful features of the image.

Note that the results depend on the chosen tubal function and transform M. We deliberately chose M to emphasize green and blue channels and diminish the effect of the red channel, enabling easy-to-interpret distinctions between tubal and entrywise functions.

We can now define a tubal softmax function based on the traditional softmax function in (8).

Definition 4.3 (tubal softmax). Consider a lateral slice X→∈ℝp×1×n as a p×1 vector of 1 × 1 × n tubes. Then, the tubal softmax function h:ℝp×1×n→Δp×1×n performs the following mapping:

for i = 1, …, p where h(X→)∈ℝp×1×n and the exponential functions are applied as tubal functions. Here, Δ^{p×1 × n} is the tubal-equivalent of p-dimensional unit simplex.

We interpret h(X→) as a vector of “tubal probabilities” in that the tubes sum to the identity tube

Through Definition 4.3, we demonstrate the parallels between tubal functions and traditional functions; the similarities are a direct consequence of a matrix-mimetic tensor framework.

The last step to define the tubal cross entropy function that converts the output of tubal function to a scalar. Recall, traditional cross entropy Lce:ℝp×{1,…,p}→ℝ+ for one sample is given by

where x is the output of the network, c is the corresponding target class, and h is the softmax function. We generalize (11) to a tubal variant as follows:

Definition 4.4 (tubal cross entropy (t-cross entropy)). The tubal cross entropy function Ltce:Δp×1×n×{1,…,p}→ℝ+ is given by

where h is the tubal softmax function, log is applied as a tubal function, c is the index corresponding to the target class, and ||·||_q is a vector norm.

The intuition behind Definition 4.4 is the following. If we have good features X→, then [h(X→)]c,1,:≈e, the identity tube, and the remaining tubes will be closer to 0. In the transform domain, [h(X→)×M]c,1,:≈1 and the remaining entries will be close to zero. When we apply the log pointwise in the transform domain, the tube log[h(X→)×M]c,1,:≈0 and the remaining entries will be large negative numbers. As a result, L_tce is smallest when the [h(X→)]c,1,:≈e, as desired.

In practice, if q-norm corresponds to a finite integer, we can instead use ||log[h(X→)]c,1,:||qq for t-cross entropy for easier derivative computations. For numerical benefits when training, we consider normalized versions based on the number of tubal entries, e.g., multiply by 1/n. These suggested modifications should not change performance in theory, but could change preferred training hyperparameters.

The workhorse of neural network training is backward propagation (Rumelhart et al., 1986; Bengio et al., 1994; Shalev-Shwartz et al., 2017; Nielsen, 2018), a method to calculate the gradient of the objective function (7) with respect to the weights. With gradient information, one can apply standard stochastic gradient optimization techniques to train.

In the ⋆_M-framework, for an orthogonal transformation M, the backpropagation formulas are analogous to the matrix case. For example, the derivatives of the ⋆_M-product are

where ∂X indicates a direction or perturbation of the same size as X. For full details on the derivation, we refer the reader to Newman (2019).

The simplicity of the back-propagation formulas (12) is one of the hallmarks of our choice of leveraging a matrix-mimetic tensor framework. Other tensor-based neural network designs (Wang et al., 2023) often require complicated indexing and non-traditional notation which, in addition to being cumbersome, can preclude extending more sophisticated neural network architectures to higher dimensions. The ⋆_M-framework yields derivative formulations that are easy to interpret, implement, and analyze.

Consider the residual neural network (He et al., 2016) with tensor operations, given by

where Y→j∈ℝm×1×n, Wj∈ℝm×m×n, and B→j∈ℝm×1×n. With the addition of the step size h, we can interpret Equation (13) as a forward Euler discretization of the continuous ordinary differential equation (ODE)

for all t ∈ [0, T] where T is the final time corresponding to the depth of the discretized network.

The stability of non-autonomous ODEs like (14) depends on the eigenvalues of the Jacobian with repsect to the features. To perform analogous analysis for tensor operators, it is useful to consider the equivalent block matrix version of the ⋆_M-product in Equation (1) where W⋆MY→ ≡ struct(W)unfold(Y→). It follows that we can matricize (14) via

The Jacobian of the matricized system (15) with respect to unfold(Y→(t)) is

where x(t)=struct(W(t))unfold(Y→(t))+unfold(B→(t)) and J(t) ∈ ℝ^mn×mn. In most cases, the activation function σ is non-decreasing, and thus the entries in diag(σ′(x(t))) are nonnegative. As a result, the stability and well-posedness of the ODE relies on the eigenvalues of struct(W(t)). As described in Haber and Ruthotto (2017), the learning problem is well-posed if

where λ_i(A) is the i-th eigenvalue of A. This criterion implies that the imaginary part of the eigenvalues drive the dynamics, promoting rotational movement of features. This produces stable forward propagation that avoids features diverging and prevents inputs from distinct classes from converging to indistinguishable points; the latter would lead to ill-posed back propagation and hence an ill-posed learning problem.

We can be more concrete about the eigenvalues because in Equation (1), struct(W(t)) is block-diagonalized in the transform domain. Thus, we can equivalently write Equation (16) as follows:

for i = 1, …, m and k = 1, …, n. In short, if the eigenvalues of each frontal slice in the transform domain have a real part close to zero, the t-NN learning problem will be well-posed, save one more requirement.

A subtle, yet important requirement for stability is that the weights change gradually over time (i.e., layers). This ensures that small perturbations of the weights yield small perturbations of the features. To promote this desired behavior, we imposed a smoothing regularizer in (discrete) time via

While theoretically useful, it is impractical to evaluate the eigenvalues of the weight tensors to satisfy the well-posedness condition (17) as we train a network. Instead, motivated by the presentation in Haber and Ruthotto (2017), we implement a forward propagation that inherently satisfies (17) independent of the weights. This forward propagation scheme is inspired by Hamiltonian dynamics, which we briefly describe and refer to Ascher (2010) and Brooks et al. (2011) for further details. We define a Hamiltonian as follows:

Definition 5.1 (Hamiltonian). Let y(t)∈ℝmy×1, z(t)∈ℝmz×1, and t ∈ [0, T]. A Hamiltonian H:ℝmy×1×ℝmz×1×[0,T]→ℝ is a system governed by the following dynamics:

Intuitively, the Hamiltonian H describes the total energy of the system with y as the position and z as the momentum or velocity. We can separate total energy into potential energy U and kinetic energy T; that is, H(y, z, t) = U(y)+T(z). This separability ensures the Hamiltonian dynamics conserve energy; i.e.,

In terms of neural networks, energy conservation (19) ensures that network features are preserved during forward propagation, thereby avoiding the issue of exploding/vanishing gradients and enabling the use of deeper networks. Additionally, Hamiltonians are symplectic or volume-preserving in the sense that the dynamics are divergence-free. For neural networks, this ensures the distance between features does not change significantly, avoiding converging and diverging behaviors. We also note that Hamiltonians are time-reversible. This ensures that if we have well-posed dynamics during forward propagation, we will have similar dynamics for backward propagation.

To preserve the benefits of Hamiltonians in the discretized setting, we symmetrize the Hamiltonian (Definition 5.1) and use a leapfrog integration method (Skeel, 1993; Ascher, 2010; Haber and Ruthotto, 2017). For t-NNs, we write the new system in terms of the ⋆_M-product (with slight abuse of notation)

where Y→(0)=Y→0 and Z→(0)=0. Here, Y→(t) indicates the data features and Z→(t) is an auxilary variable not related to the data directly. We add the bias tube, b(t), to each element. The equivalent matricized version of Equation (20) is

The (matricized) system (21) is inherently stable, independent of the weight tensors W(t), because of the block antisymmetric structure. The eigenvalues of antisymmetric matrices are purely imaginary, which exactly satisfy the stability condition in Equation (17).

We discretize Equation (20) using the leapfrog method, a symplectic integration technique, defined as

for j = 0, …, d−1. We demonstrate the benefits of stable forward propagation (22) in Example 5.2.

Example 5.2 (Trajectories of stable t-NNs). We construct a dataset in ℝ³ randomly drawn from a multivariate normal distribution with mean of 0 and a covariance matrix of 3I₃. The points are divided into three classes based on distance to the origin with yellow points inside a sphere with radius r = 3.5, green points inside a sphere with radius R = 5.5, and purple points outside both spheres. We train with 1200 data points and store the data as 1 × 1 × 3 tubes.

We forward propagate using one of two integrators with weights and biases as 1 × 1 × 3 tubes:

We create a network with d = 32 layers and use a step size of h = 2. We use the discrete cosine transform for M. We train for 50 epochs using Adam (Kingma and Ba, 2015) with a batch size of 10 and a learning rate of γ = 10⁻². To create smoother dynamics, we regularize the weights using Equation (18) with regularization parameter λ = 10⁻⁴/h. We illustrate the dynamics of trained networks in Figure 4.

The dynamics in Figure 4 show the topological benefits of leapfrog integration. The data points exhibit smoother, rotational dynamics to reach the desired linearly-separable final configuration. In comparison, forward Euler propagation significantly changes topology during forward propagation. Such topological changes may yield ill-posed learning problems and poor network generalization.

We implement both tensor and matrix frameworks using PyTorch (RRID:SCR_018536) (Paszke et al., 2017). All of the code to reproduce the experiments is provided in https://github.com/elizabethnewman/tnn. All experiments were run on an Exxact server with four RTX A6000 GPUs, each with 48 GB of RAM. Only one GPU was used to generate each result. The results we report are for the networks that yielded the best accuracy on the validation data.

We train all models using the stochastic gradient method Adam (Kingma and Ba, 2015), which uses a default learning rate of 10⁻³. In most cases, we select hyperparameters and weight initialization based on the default settings in PyTorch. We indicate the few exceptions to this at the start of the corresponding sections. We also pair the stochastic optimizers with a learning rate scheduler that decreases the learning rate by a factor of γ every M steps. In all cases, we used the default γ = 0.9 and M = 100. This is a common practice to ensure convergence of stochastic optimizer in idealized settings (Bottou et al., 2018). We utilize a weight decay parameter in some experiments to reduce overfitting.

For the tensor networks in all experiments, we use the discrete cosine transform for M, which is a close, real-valued version of the discrete Fourier transform (DFT). The DFT matrix corresponds to the t-product, which has been shown to be effective for natural image applications (Kilmer and Martin, 2011; Hao et al., 2013; Newman et al., 2018).

The MNIST dataset (LeCun et al., 2010) is composed of 28 × 28 grayscale images of handwritten digits. We train on 50, 000 images and reserve 10, 000 for validation. We report test accuracy on 10, 000 images not used for training nor validation. For NNs, we vectorize images and store as columns, resulting in a matrix of size 28²×b where b is the number of images. For t-NNs, we store the images as lateral slices, resulting in a tensor of size 28 × b×28.

In this experiment, we train an autoencoder to efficiently represent the high-dimensional MNIST data in a low-dimensional subspace. Autoencoders can be thought of as nonlinear, parameterized extensions of the (truncated) singular value decomposition. Our goal is to solve the (unsrpervised) learning problem

where fenc:Y→Z is the encoder and fdec:Z→Y is the decoder. Here, Z is the latent space that is smaller than the data space; in terms of dimension, we say dim(Z)<dim(Y). Note that the mean squared error (MSE) tubal loss is the same as the MSE loss function when using an orthogonal transformation matrix, as we have done in our experiments. We describe the NN and t-NN autoencoder architectures used in Figure 5 and report results in Figure 6.

We observe that the t-NN outperforms the NN autoencoders with similar numbers of weights with an order of magnitude smaller training and validation loss and test error as well as qualitative improvements of the approximations. The neural network autoencoder with the same feature space dimensions, NN(560,280), performs best in terms of the loss and error metrics, but requires over 20 times more network weights than the t-NN autoencoder. The t-NN layers are able capture spatial correlations more effectively using multilinear operations, resulting in quality approximations with significantly fewer weights.

We use the same MNIST dataset as for the autoencoder example. We train for 20 epochs using Adam with a batch size of 32 and a learning rate of 10⁻². We add Tikhonov regularization (weight decay) with a regularization parameter of λ = 10⁻⁴. We use the PyTorch defaults for the other optimizer hyperparameters. For the t-cross entropy loss, we use a squared ℓ₂-norm and normalize by the number of entries in the tube.

We compare four different two-layer neural network architectures, described in Table 1. We use either cross entropy loss or t-cross entropy loss, depending on the architecture.

We report the convergence and accuracy results in Figure 7 and Table 2, respectively.

The t-NN architecture with cross entropy loss outperforms all networks in terms of test accuracy and accuracy per class. The second-best performing network is the t-NN with t-loss. These results are evidence that the features learned from the tensor linear layer (layer 1) are better than those learned by a dense matrix layer. We further note that matrix network NN with square weights has the same final layer shape as the t-NN with cross entropy loss; the only difference between the networks is the first layer. We depict the learned features of each network that preserves the size of the images in Figure 8. The t-NN features from the first layer contain more structure than the NN features with square weights. This reflects the ⋆_M-operation, which first acts along the tubes (rows of the images in Figure 8). We also observe that the features of NN are more extreme, close to +1 and −1, the limits of the range of the activation function σ(x) = tanh(x). In comparison, the features extracted from the t-NN with t-loss offer more variety of entries, but still hits the extreme values often. This demonstrates that t-NNs still produce rich feature spaces and are able to achieve these features with about 20 times fewer weights.

We note that the t-NN network with t-cross entropy loss performs well, and we gain insight into the beneifts of t-losses from the accuracy per class in Figure 9. We observe that when we use tubal losses, we require high values for many frontal slices. This creates a more rigorous classification requirement and, as we will see in subsequent experiments, can yield networks that generalize better. Additionally, the distribution of values in the tubal softmax function is reflective of the predicted class. For the second image (true = 3), the two most likely predicted classes were, in order, 3 and 8. Qualitatively, the particularly handwritten 3 has similarities to the digit 8, and the tubal softmax captures this similarity. For the cases that were incorrectly predicted the digit, the handwritten image had structure emblematic of the predicted class and second most likely predicted class matched the true label.

The CIFAR10 dataset (Krizhevsky and Hinton, 2009) is composed of 32 × 32 × 3 RGB natural images belonging to ten classes. We train on 40, 000 images and reserve 10, 000 for validation. We report test accuracy on 10, 000 images not used for training nor validation. For NNs, we vectorize images and store as columns, resulting in a matrix of size (3·32²) × b where b is the number of images. For t-NNs, we store the images as lateral slices and stack the color channels vertically, resulting in a tensor of size (3·32) × b×32. We train for 500 epochs using Adam with a batch size of 32. We use a learning rate of 10⁻³ that decays after every 100 epochs by a factor of 0.9. We use the PyTorch defaults for the other optimizer hyperparameters. For the t-cross entropy loss, we use a squared ℓ₂-norm and normalize by the number of entries in the tube.

We compare the performance of the Hamiltonian networks with dense matrix operators and dense tensor operators for various numbers of layers (d = 4, 8). In conjunction with the Hamiltonian network, we use the smoothing regularizer (20) with a regularization parameter of λ = 10⁻². We describe the network architectures and number of parameters in Table 3. The NN architectures require more than 30 times the number of weights than the t-NN architectures. We compare the convergence and accuracy results in Figure 10 and Table 4, respectively.

There are several key takeaways from the numerical results. First, the depth of the network did not significantly change performance in this experiment. We state this observation cautiously. We observe this behavior for a certain set of fixed hyperparameters (e.g., step size h, learning rate, regularization parameter λ,...). The interaction of the hyperparameters and performance is complex and a complete ablation study is outside of the scope of this paper. A second takeaway is that the t-NN trained with the tubal loss generalizes better the NN networks and better than t-NNs with cross entropy loss (not shown for simplicity; see Figure 10 for details). This behavior is especially apparent when looking at the test loss in Table 4. The t-NN with four Hamiltonian layers and t-loss performs well overall, obtaining almost 5% better overall test accuracy. In comparison, the matrix NN quickly overfits the training data and thus does not generalize as well. In terms of the test accuracy per class, the t-NN architectures achieve the best performance in all but two classes.

To look into the performance further, we examine the extracted features of Hamiltonian NNs and t-NNs in Figure 11. We see that the NN and t-NN features share similarities. Both features gradually remove the structure of the original image at similar rates. The t-NN architecture achieves this pattern with significantly fewer network weights (over 30 times fewer). The noisy artifacts differ between the two architectures. In particular, we see that the t-NN layers produce blockier artifacts because of the structured ⋆_M-operation.

In the last numerical study in Figure 12, we explore how quickly we can train Hamiltonian t-NNs compared to NNs. In addition to an order of magnitude fewer weights, training t-NNs takes less time than training NNs. As we increase the depth of the networks, we see that each NN epoch takes approximately 1.75 times longer to complete. This performance could potentially be further improved if we optimized the ⋆_M-product, e.g., using fast transforms instead of matrix multiplication.

The CIFAR100 dataset (Krizhevsky and Hinton, 2009) is composed of 32 × 32 × 3 RGB natural images belonging to 100 classes. We train on 40, 000 images and reserve 10, 000 for validation. We report test accuracy on 10, 000 images not used for training nor validation. We use the same setup and training parameters as the CIFAR10 experiment (Section 6.4). For all experiments, we use Hamiltonian networks with a depth of d = 16, a step size of h = 0.25, and a regularization parameter λ = 10⁻². We report the accuracy results in Figure 13.

Observing the results in Figure 13, the conclusions are the same as in the CIFAR-10 experiment. Specifically, training with the t-NN and t-loss produces the best test accuracy, about a 4% improvement from the comparable NN network. This demonstrates that the benefits of dense tensor operations over dense matrix operations can be realized for more challenging classification problems and motivates further development of these tools to improve state-of-the-art convolutional neural networks and other architectures.