Striking The Right Balance: Three-Dimensional Ocean Sound Speed Field
  Reconstruction Using Tensor Neural Networks
Accurately reconstructing a three-dimensional ocean sound speed field (3D
SSF) is essential for various ocean acoustic applications, but the sparsity and
uncertainty of sound speed samples across a vast ocean region make it a
challenging task. To tackle this challenge, a large body of reconstruction
methods has been developed, including spline interpolation, matrix/tensor-based
completion, and deep neural networks-based reconstruction. However, a
principled analysis of their effectiveness in 3D SSF reconstruction is still
lacking. This paper performs a thorough analysis of the reconstruction error
and highlights the need for a balanced representation model that integrates
both expressiveness and conciseness. To meet this requirement, a 3D
SSF-tailored tensor deep neural network is proposed, which utilizes tensor
computations and deep neural network architectures to achieve remarkable 3D SSF
reconstruction. The proposed model not only includes the previous tensor-based
SSF representation model as a special case, but also has a natural ability to
reject noise. The numerical results using the South China Sea 3D SSF data
demonstrate that the proposed method outperforms state-of-the-art methods. The
code is available at https://github.com/OceanSTARLab/Tensor-Neural-Network.
Siyuan Li, Lei Cheng, Ting Zhang, Hangfang Zhao, Jianlong Li
Striking The Right Balance: Three-Dimensional Ocean Sound Speed Field Reconstruction Using Tensor Neural Networks

###### Abstract

Accurately reconstructing a three-dimensional ocean sound speed field (3D SSF) is essential for various ocean acoustic applications, but the sparsity and uncertainty of sound speed samples across a vast ocean region make it a challenging task. To tackle this challenge, a large body of reconstruction methods has been developed, including spline interpolation, matrix/tensor-based completion, and deep neural networks-based reconstruction. However, a principled analysis of their effectiveness in 3D SSF reconstruction is still lacking. This paper performs a thorough analysis of the reconstruction error and highlights the need for a balanced representation model that integrates both expressiveness and conciseness. To meet this requirement, a 3D SSF-tailored tensor deep neural network is proposed, which utilizes tensor computations and deep neural network architectures to achieve remarkable 3D SSF reconstruction. The proposed model not only includes the previous tensor-based SSF representation model as a special case, but also has a natural ability to reject noise. The numerical results using the South China Sea 3D SSF data demonstrate that the proposed method outperforms state-of-the-art methods. The code is available at [https://github.com/OceanSTARLab/Tensor-Neural-Network](https://github.com/OceanSTARLab/Tensor-Neural-Network).

[[https://doi.org](https://doi.org)(DOI number)]

Pages: 1-18

## I Introduction

Three-dimensional (3D) ocean sound speed fields (SSFs), which characterize sound propagation over 3D geographical regions [1], form the stepping stone towards the success of a myriad of ocean acoustic applications, including underwater detection, [2] localization, [3] and communications. [4] To realize sound speed awareness, ocean observing systems have been rapidly developing in recent years, leading to the proliferation of intelligent underwater floats [5] and vehicles. [6] Despite these advancements, sound speed samples are still sparsely scattered across a vast ocean region (see Fig. 1), making crafting an accurate and fine-grained 3D SSF a considerable challenge.

The crux of the matter is that the number of samples is much less than the unknowns, resulting in a highly under-determined inverse problem. [7] To make the problem tractable, the primary approach involves modeling the prior knowledge of the unknowns and incorporating that information into the reconstruction process. This has led to a sizable body of related works on reconstruction methods, [8; 9; 10; 11; 12] not limited to SSF reconstruction, that have a near-universal three-step approach: 1) proposing a representation model that incorporates prior information, 2) learning the model parameters from limited samples, and 3) reconstructing the underlying signal using the learned model.

For instance, the simple yet widely-adopted spline interpolation method assumes that the underlying signal is a linear combination of smooth Green functions centered at the samples, [13] which correspond to the linear regression model in machine learning (see detailed discussion in Sec. II of Ref. [8]). Using linear regression as a representation model, the spline interpolation method starts by learning the optimal weights for the Green functions from the samples, and then uses the learned model to reconstruct the unknowns. A step forward is the Gaussian process regression (GPR), [14] which models the underlying signal as a Gaussian random process with a kernel function describing spatial correlations. GPR first learns the hyper-parameters of the model from the samples and then reconstructs the unobserved values. In recent years, matrix/tensor decompositions and deep neural networks have been increasingly used as representation models for various data completion tasks, leading to state-of-the-art (SOTA) results. Our task is closely related to recent advances in 3D data inverse problems and field reconstruction, which includes video/image inpainting [15; 10], hyperspectral reconstruction [11; 16], and radio map cartography [17; 12].

The burgeoning literature on reconstruction methods invites the question: _Among numerous ways, which one
is the most effective for 3D ocean SSF reconstruction?_ Rather than exhaustively comparing various methods, this paper sets out to answer this question in a principled manner by delving into the analysis of reconstruction error. Specifically, inspired by the well-known bias-variance trade-off in machine learning,[18] the reconstruction error can be mainly decomposed into two parts: the _representation error_, which measures the model's capacity to fit the underlying signal (e.g., 3D SSF), and the _identification error_, which assesses whether the desired model can be uniquely learned from limited samples.

To clarify this notion, consider the case where an over-parameterized deep neural network (i.e., one with more parameters than the sound speed samples) is used as the representation model. Its high expressive power[19] makes it likely to precisely represent the underlying 3D SSF, resulting in negligible representation error. However, insufficient samples and numerous unknown model parameters can lead to difficulties in identifying the optimal parameter configuration across a vast solution space. Unless strong regularizations (e.g., early stopping[20]) are carefully employed, the overwhelming identification error may lead to poor reconstruction. On the other hand, choosing a simpler representation model with fewer parameters makes it easier to uniquely learn the model parameters from limited samples, but may result in a higher representation error.

The reconstruction error analysis above _highlights the need for a balanced representation model_ in 3D SSF reconstruction, with both high expressive power to reduce representation error and conciseness (in terms of both parameter numbers and mathematical operations) to minimize identification error. To meet this requirement, we propose a 3D SSF-tailored tensor deep neural network, drawing on the advantages of both recent _succinct tensor-based_ ocean SSF representation[21, 22] and _powerful deep neural network-based_ image representation.[23] The network, shown in Fig. 1, aims to provide a concise and accurate representation of SSFs.

The proposed 3D SSF-tailored tensor deep neural network leverages the strengths of tensor computations and deep neural network architectures to achieve remarkable 3D SSF reconstruction. The tensor computations, which are with a few parameters and serve as the model's backbone, effectively exploit the multi-dimensional correlations among sound speeds.[21, 22] On the other hand, the enrollment of deep neural network architectures enhances the model's expressive power, especially in capturing fine-grained variations of sound speeds. Due to its conciseness, a simple gradient descent algorithm[24] is sufficient to identify the parameters that yield a quite good reconstruction from limited samples. Relying on its high expressive power, the learned model successfully reconstructs 3D SSFs with intricate details. Two key features of the proposed method are highlighted in this paper. First, the model has a natural ability to reject noise. This property is backed by theoretical proof under the one-layer assumption, and experimental verification in multiple layer scenarios. Second, the proposed model includes the tensor Tucker decomposition model as a specific case, which is strongly linked to classical SSF basis functions such as empirical orthogonal functions (EOFs).[22]

Moreover, the proposed model can easily incorporate additional structural assumptions to boost the reconstruction performance further. As an example, we integrate total variation[25] regularizations in this paper to exploit the spatial smoothness of SSF. Finally, to corroborate the advantages of the proposed model, we conduct comprehensive numerical experiments using real-life 3D SSF data and compare the performance of the proposed approach with SOTA reconstruction methods. The encouraging results of the proposed algorithm not only affirm our theoretical analysis, but also demonstrate the effectiveness of viewing the reconstruction task from the reconstruction error analysis perspective.

_Notation_: Lower- and uppercase bold letters (e.g., \(\mathbf{x}\) and \(\mathbf{X}\)) are used to denote the vectors and matrices. Higher-order tensors (order three or higher) are denoted by upper-case bold calligraphic letters. For a tensor \(\boldsymbol{\mathcal{X}}\), \(\mathbf{X}_{(p)}\) stands for its model-\(p\) unfolding matrix. \(\boldsymbol{\mathcal{X}}\times_{p}\mathbf{A}\) denotes the \(p\)-mode product between tensor \(\boldsymbol{\mathcal{X}}\) and matrix \(\mathbf{A}\). The Kronecker product is denoted by \(\otimes\). The Hadamard product is denoted by \(*\). \(\|\cdot\|_{\mathrm{F}}\) stands for the Frobenius norm. \(\langle\cdot,\cdot\rangle\) denotes the tensor inner product. The superscript \(\mathrm{T}^{\mathrm{T}}\) stands for transposition. \(\mathbb{R}\) and \(\mathbb{C}\) are the field of real and complex numbers, respectively.

Figure 1: Illustration of the reconstruction system considered in this paper, with the proposed 3D SSF-tailored tensor deep neural network serving as the representation model.
## II Reconstruction and error analysis

In this section, we formulate the ocean 3D SSF reconstruction task as an optimization problem, and present a unified framework to analyze the reconstruction error. Based on this framework, two widely used reconstruction methods are analyzed.

### Reconstruction Problem Formulation

The reconstruction system considered in this paper is illustrated in Fig. 1. Assuming that the sound speed samples are contaminated by independently and identically distributed (i.i.d) Gaussian noise, the SSF reconstruction problem can be formulated as

\[\min_{\boldsymbol{\mathcal{X}}}\|\boldsymbol{\mathcal{Y}}-\boldsymbol{ \mathcal{O}}\ast\boldsymbol{\mathcal{X}}\|_{\mathrm{F}}^{2}, \tag{1}\]

where \(\boldsymbol{\mathcal{X}}\in\mathbb{R}^{I\times J\times K}\) is the ground-truth 3D SSF (to be reconstructed) and \(\boldsymbol{\mathcal{O}}\) is a binary indicator tensor with \(\boldsymbol{\mathcal{O}}_{i,j,k}=1\) if entry \((i,j,k)\) is observed. Observation tensor \(\boldsymbol{\mathcal{Y}}\in\mathbb{R}^{I\times J\times K}\) collects noisy and limited SSF samples, i.e., \(\boldsymbol{\mathcal{Y}}_{i,j,k}\) is the sampled sound speed value if \(\boldsymbol{\mathcal{O}}_{i,j,k}=1\), and otherwise equals to zero.

Denote \(N\) as the number of observed (non-zero) entries in \(\boldsymbol{\mathcal{Y}}\), and \(T=IJK\) as the total number of entries in \(\boldsymbol{\mathcal{X}}\). The SSF reconstruction problem in Eq. (1) is under-determined because \(N\) is usually much smaller than \(T\). To tackle this challenge, existing methods (e.g., Refs. [23, 26, 8, 13, 27]) usually take the following three steps.

#### Ii-A1 SSF Representation

Since ocean sound speeds exhibit strong spatial correlations, SSF can be effectively represented by one model with only _a few_ parameters, as demonstrated in prior literature. [22, 28] Mathematically, assume that ocean SSF is represented by:

\[\boldsymbol{\mathcal{X}}\approx\boldsymbol{\mathcal{D}}(\boldsymbol{\Theta}), \tag{2}\]

where \(\boldsymbol{\mathcal{D}}(\cdot)\) denotes an SSF representation model and \(\boldsymbol{\Theta}\) is the set of model parameters.

#### Ii-A2 Parameter Learning

Using the representation model \(\boldsymbol{\mathcal{D}}(\cdot)\), the SSF reconstruction problem defined in Eq. (1) can be recast as

\[\min_{\boldsymbol{\Theta}}\|\boldsymbol{\mathcal{Y}}-\boldsymbol{\mathcal{O}} \ast\boldsymbol{\mathcal{D}}(\boldsymbol{\Theta})\|_{\mathrm{F}}^{2}. \tag{3}\]

The model parameter estimate \(\hat{\boldsymbol{\Theta}}\) is then obtained by solving problem (3).

#### Ii-A3 SSF Reconstruction

Using the learned model parameter \(\hat{\boldsymbol{\Theta}}\), the SSF is reconstructed by \(\hat{\boldsymbol{\mathcal{X}}}=\boldsymbol{\mathcal{D}}(\hat{\boldsymbol{ \Theta}})\).

### Reconstruction Error Analysis

We proceed to conduct a theoretical analysis of the reconstruction error, which is defined as \(E=\|\hat{\boldsymbol{\mathcal{X}}}-\boldsymbol{\mathcal{X}}\|_{\mathrm{F}}^{2}\). Since the reconstruction error is a least-squares function, its associated error analysis, such as the bias-variance trade-off, has already been developed in the machine learning and signal processing literature. For example, refer to Page 149 of Ref. [18]. The existing results are primarily developed for supervised regression tasks, but here we re-interpreted them in the context of our unsupervised SSF reconstruction task.

Specifically, denote the set of parameters that can best represent the SSF data as \(\boldsymbol{\Theta}^{*}\), i.e., \(\boldsymbol{\Theta}^{*}=\min_{\boldsymbol{\Theta}}\|\boldsymbol{\mathcal{X}}- \boldsymbol{\mathcal{D}}(\boldsymbol{\Theta})\|_{\mathrm{F}}^{2}\). The reconstruction error can be shown (see Appendix B) to be decomposed as

\[E=E_{1}+E_{2}+\epsilon, \tag{4}\]

where \(E_{1}=\|\boldsymbol{\mathcal{X}}-\boldsymbol{\mathcal{D}}(\boldsymbol{\Theta}^ {*})\|_{\mathrm{F}}^{2}\) is the representation error; is the cross term defined as \(\epsilon=2\langle\boldsymbol{\mathcal{X}}-\boldsymbol{\mathcal{D}}(\boldsymbol{ \Theta}^{*}),\boldsymbol{\mathcal{D}}(\boldsymbol{\Theta}^{*})-\boldsymbol{ \mathcal{D}}(\hat{\boldsymbol{\Theta}})\rangle\).

From the definition of the three terms in Eq. (4), it is easy to show (see Appendix B) that the cross term \(\epsilon\) is equal to zero if \(E_{1}\) or \(E_{2}\) is equal to zero, i.e.,

\[\epsilon=0\Leftarrow\boldsymbol{\mathcal{X}}=\boldsymbol{\mathcal{D}}( \boldsymbol{\Theta}^{*})\text{ or }\boldsymbol{\mathcal{D}}(\boldsymbol{\Theta}^{*})=\boldsymbol{\mathcal{D}}( \hat{\boldsymbol{\Theta}}). \tag{5}\]

The error decomposition results show that minimizing the reconstruction error \(E\) requires simultaneously reducing the representation error \(E_{1}\) and the identification error \(E_{2}\). Like neural networks, models with high expressive powers usually have low representation errors but high identification errors. Conversely, concise models, like linear regression, have low identification errors but high representation errors due to their limited parameters and simple operations.

The result in (5) suggests two possible ways to minimize the reconstruction error. To null the cross term \(\epsilon\), one can choose a complicated model with universal approximation property [19] to zero the representation error; or a relatively simple model to induce zero identification error. In the first case, the reconstruction error is up to the identification error, thus necessitating the model complexity optimization. In the latter case, the effort should be paid to enhance the expressive power of the model while not increasing the model's complexity.

To draw more insights into the reconstruction error, two commonly used reconstruction methods are analyzed in the following subsection.

### Commonly Used Approaches

#### Ii-C1 Spline Interpolation

The biharmonic spline interpolation assumes that the SSF can be well represented by a linear combination of smooth Green functions centered at the samples, i.e., [13]

\[\boldsymbol{\mathcal{X}}_{\mathbf{i}}=[\boldsymbol{\mathcal{D}}(\mathbf{w})]_ {\mathbf{i}}=\sum_{n=1}^{N}w_{n}g(\mathbf{i},\mathbf{i}_{n}), \tag{6}\]

where \(\mathbf{i}_{n}=(i_{n},j_{n},k_{n})\in\mathbb{R}^{3}\) is the index of the sampled location; \(w_{n}\) is the weight associated with each sample; \(\boldsymbol{\Theta}=\mathbf{w}=[w_{1},\cdots,w_{N}]^{\mathrm{T}}\); and \(g(\cdot,\cdot)\) is the Green function. The representation error \(E_{1}\) is given by

\[E_{1}=\|\boldsymbol{\mathcal{X}}-\boldsymbol{\mathcal{D}}(\mathbf{w}^{*})\|_{ \mathrm{F}}^{2}, \tag{7}\]
where the optimal weights \(\mathbf{w}^{*}\) are determined by solving the following problem:

\[\min_{\mathbf{w}}\sum_{n=1}^{N}\left[\boldsymbol{\mathcal{Y}}_{\mathbf{l}_{n}}- \sum_{k=1}^{N}w_{k}g(\mathbf{l}_{n},\mathbf{l}_{k})\right]^{2}. \tag{8}\]

Note that the optimal weights can be determined uniquely from Eq. (8) since it is a well-posed problem, which leads to zero identification error, i.e., \(E_{2}=0\). However, the predefined Green function is not flexible and may not reflect the correlations among the sound speeds at different regions, resulting in a significant \(E_{1}\). According to the results in (4) and (5), since \(E_{2}=0\) and \(\epsilon=0\), the reconstruction error \(E=E_{1}\).

## 2 Deep Neural Networks

In recent years, deep neural networks (DNN) have become highly effective tools in various areas, including computer vision[29], natural language processing[30], and acoustic signal processing[31, 32, 33]. Here we assume that the SSF can be represented by a \(L\)-layer neural network, i.e.,

\[\boldsymbol{\mathcal{X}}=\boldsymbol{\mathcal{D}}(\boldsymbol{\Theta})=F_{ \boldsymbol{\theta}_{L}}(F_{\boldsymbol{\theta}_{L-1}}(\cdots F_{\boldsymbol{ \theta}_{1}}(\boldsymbol{\mathcal{G}}))), \tag{9}\]

where \(F_{\boldsymbol{\theta}_{1}}(\cdot)\) denotes the function of the \(l\)th layer with parameters (e.g., weights and biases) \(\boldsymbol{\theta}_{l}\); and \(\boldsymbol{\mathcal{G}}\) is the core tensor. The parameter set \(\boldsymbol{\Theta}\) includes \(\boldsymbol{\mathcal{G}}\) and \(\{\boldsymbol{\theta}_{l}\}_{l}\). For the functional forms of \(\{F_{\boldsymbol{\theta}_{l}}(\cdot)\}_{l}\), commonly used ones include the multi-layer perceptron (MLP)[34] and convolutional neural network (CNN)[35].

For the DNN-based SSF representation model, the identification error is given by

\[E_{2}=\|F_{\boldsymbol{\theta}_{L}^{*}}(F_{\boldsymbol{\theta}_{L-1}^{*}}( \cdots F_{\boldsymbol{\theta}_{1}^{*}}(\boldsymbol{\mathcal{G}}^{*})))-F_{ \hat{\boldsymbol{\theta}}_{L}}(F_{\hat{\boldsymbol{\theta}}_{L-1}}(\cdots F_{ \hat{\boldsymbol{\theta}}_{1}}(\hat{\boldsymbol{\mathcal{G}}})))\|_{\mathrm{F }}^{2}, \tag{10}\]

where \(\boldsymbol{\Theta}^{*}\) (includes \(\boldsymbol{\mathcal{G}}^{*}\) and \(\{\boldsymbol{\theta}_{l}^{*}\}_{l}\)) is the solution to the representation problem

\[\min_{\boldsymbol{\Theta}}\|\boldsymbol{\mathcal{X}}-F_{\boldsymbol{\theta}_{ L}}(F_{\boldsymbol{\theta}_{L-1}}(\cdots F_{\boldsymbol{\theta}_{1}}( \boldsymbol{\mathcal{G}})))\|_{\mathrm{F}}^{2}; \tag{11}\]

and \(\hat{\boldsymbol{\Theta}}\) (includes \(\hat{\boldsymbol{\mathcal{G}}}\) and \(\{\hat{\boldsymbol{\theta}}_{l}\}_{l}\)) is the solution to the reconstruction problem

\[\min_{\boldsymbol{\Theta}}\|\boldsymbol{\mathcal{Y}}-\boldsymbol{\mathcal{O}}* F_{\boldsymbol{\theta}_{L}}(F_{\boldsymbol{\theta}_{L-1}}(\cdots F_{\boldsymbol{ \theta}_{1}}(\boldsymbol{\mathcal{G}})))\|_{\mathrm{F}}^{2}. \tag{12}\]

Due to the universal approximation property[19] of DNN, the representation error \(E_{1}\) can be zero, i.e., \(E_{1}=0\), if a sufficient number of neurons are available and/or a highly sophisticated neural architecture is utilized. According to the results in (4) and (5), since \(E_{1}=0\) and \(\epsilon=0\), the reconstruction error \(E=E_{2}\), which, however, will be large since the solution space of problem (12) is vast (that is, there exists a large number of solutions that attain the same objective value of problem (12)).

## III Tensor Neural Network-Aided Reconstruction

The reconstruction error analysis introduced in the last section shows that a highly expressive representation model is necessitated to accurately capture the complex details of sound speeds' spatial distribution, thus minimizing the representation error. But this often leads to increased parameters and complicated mathematic operations, resulting in a more complicated model with increased identification error. Consequently, the key to approaching the optimal reconstruction is to enhance the model's expressive power while maintaining the model's conciseness (in terms of both parameter number and mathematical operations). Towards this goal, in this section, inspired by the recent success of succinct tensor-based ocean SSF representation and powerful deep neural network-based image representation, we _take the best from them_ and propose a 3D SSF-tailored tensor deep neural network. Then, an effective reconstruction algorithm is derived.

In the following subsections, we first briefly review the recent tensor-based SSF representation model (Sec. III.1), and then devise our proposed model that integrates the tensor model and deep neural network (Sec. III.2). Next, theoretical insights are drawn by analyzing the proposed model's noise rejection property and particular form (Sec. III.3). Finally, the reconstruction algorithm is derived (Sec. III.4).

### Tensor-based SSF Representation

This subsection briefly reviews the tensor-based SSF representation framework established in Ref. [22], which relies on the Tucker decomposition model. For readers unfamiliar with tensors, a brief review in the context of SSF representation was provided in Sec. III of Ref. [22]. The definitions of the tensor operations utilized in this paper are elucidated in Appendix A.

Mathematically, the tensor-based SSF representation model is:

\[\boldsymbol{\mathcal{X}}\approx\boldsymbol{\mathcal{S}}\times_{1}\mathbf{U}^{( 1)}\times_{2}\mathbf{U}^{(2)}\times_{3}\mathbf{U}^{(3)}, \tag{13}\]

where the columns of factor matrices \(\{\mathbf{U}^{(p)}\}_{p=1}^{3}\) can be interpreted as the spatial basis functions; and core tensor \(\boldsymbol{\mathcal{S}}\) contains the weighting coefficients. Symbol \(\times_{p}\) denotes the mode-\(p\) product; see the definition in Appendix A. Various constraints can be specified on factor matrices and the core tensor to incorporate additional knowledge, e.g., orthogonal constraints[22] and sparsity constraints[21].

Ref. [22] has shown the conciseness of the tensor-based representation model, which is mathematically succinct and has a small number of parameters, leading to a small identification error \(E_{2}\). On the other hand, the tensor-based representation model is also effective in characterizing the three-dimensional interactions among sound speeds. Notably, it includes the classical SSF basis functions, such as EOFs and Fourier basis functions, as special cases[22]. However, since the Tucker decomposition model is based on a multi-linear form, the representation model in Eq. (13) is difficult to capture the highly nonlinear variations of sound speeds (caused by small-scale/mesoscale ocean processes), resulting in a non
negligible representation error \(E_{1}\), as demonstrated by the numerical results presented in Sec. IV.2.

### Integrating Tensor and Neural Networks for SSF Representation

In this subsection, we propose an integrated model that can take both advantages of DNN and tensor-based representation models, i.e., expressiveness and conciseness, respectively. The proposed model is based on the tensor contraction layer (TCL),[36] as shown in Fig. 2.

Given an input tensor \(\mathbf{\mathcal{X}}_{l}\in\mathbb{R}^{R_{1}^{l}\times R_{2}^{l}\times R_{3}^{l}}\), the output of TCL is:

\[\mathbf{\mathcal{X}}_{l+1}=\varsigma(\mathbf{\mathcal{X}}_{l}\times_{1}\mathbf{W}_{l}^ {(1)}\times_{2}\mathbf{W}_{l}^{(2)}\times_{3}\mathbf{W}_{l}^{(3)}), \tag{14}\]

where the output \(\mathbf{\mathcal{X}}_{l+1}\in\mathbb{R}^{R_{1}^{l+1}\times R_{2}^{l+1}\times R_{3} ^{l+1}}\); each matrix \(\mathbf{W}_{l}^{(i)}\in\mathbb{R}^{R_{i}^{l+1}\times R_{i}^{l}}\); and \(\varsigma(\cdot)\) is the activation function. Note that the TCL is based on the Tucker decomposition, enabling it to effectively capture the multidimensional interactions of sound speeds even with a small number of parameters.

To further enhance the expressive power, we propose a tensor neural network (TNN) composed of multiple TCLs; see Fig. 3. Mathematically, a \(L\)-layer TNN can be defined as

\[\mathbf{\mathcal{X}}=\varsigma\bigg{(} \cdots\varsigma\left(\mathbf{\mathcal{G}}\times_{1}\mathbf{W}_{1}^{( 1)}\times_{2}\mathbf{W}_{1}^{(2)}\times_{3}\mathbf{W}_{1}^{(3)}\right)\cdots \tag{15}\] \[\times_{1}\mathbf{W}_{L}^{(1)}\times_{2}\mathbf{W}_{L}^{(2)} \times_{3}\mathbf{W}_{L}^{(3)}\bigg{)},\]

where \(\mathbf{\mathcal{G}}\in\mathbb{R}^{R_{1}\times R_{2}\times R_{3}}\) is the core tensor, \(\{\{\mathbf{W}_{l}^{(i)}\}_{i=1}^{3}\}_{l=1}^{L}\) are the factor matrices. The choice of activation function depends on the data. The commonly used ones are rectified linear unit (ReLU), sigmoid, tanh, etc.[37].

### Theoretical Insights

Before deriving the reconstruction algorithm, theoretical analyses were performed to reveal the insights of the proposed TNN-based representation model.

First, we show that the particular architecture of the proposed model enables it to reject noise effectively. Here the Gaussian noise is considered due to the least-squares loss function used in Eq. (1), which is a common choice in physical field reconstruction research[12, 38]. For simplicity, consider a _one-layer_ TNN equipped with a ReLU activation function. The TNN output can be written as:

\[\mathbf{\mathcal{X}}=\text{ReLU}(\mathbf{\mathcal{G}}\times_{1}\mathbf{W}_{1}^{(1)} \times_{2}\mathbf{W}_{1}^{(2)}\times_{3}\mathbf{W}_{1}^{(3)}). \tag{16}\]

Then we have the following proposition.

**Proposition 1**: Consider a _one-layer_ TNN with parameters \(\mathbf{\Theta}=\{\mathbf{\mathcal{G}},\mathbf{W}_{1}^{(1)},\mathbf{W}_{1}^{(2)}, \mathbf{W}_{1}^{(3)}\}\) and a noise tensor \(\mathbf{\mathcal{E}}\) with each element following a zero-mean i.i.d Gaussian distribution, i.e., \(\mathbf{\mathcal{E}}_{i,j,k}\sim\mathcal{N}(0,\sigma^{2})\). Then with probability at least \(1-e^{-R}-e^{-\frac{T}{16}}\),

\[\min_{\mathbf{\Theta}}\|\mathbf{\mathcal{X}}-\mathbf{\mathcal{E}}\|_{\rm F}^{2}\geq\|\bm {\mathcal{E}}\|_{\rm F}^{2}(1-\frac{10R}{T}), \tag{17}\]

where \(R=R_{1}R_{2}R_{3}\) and \(T=IJK\).

_Proof_: See Appendix B.

The term \(\|\mathbf{\mathcal{X}}-\mathbf{\mathcal{E}}\|_{\rm F}^{2}\) in Eq. (17) serves as a metric for assessing the model's capacity to fit the noise. A higher value of this term indicates that the model is less prone to fit the noise. **Proposition 1** provides a lower bound for \(\|\mathbf{\mathcal{X}}-\mathbf{\mathcal{E}}\|_{\rm F}^{2}\), which is a product of the noise power \(\|\mathbf{\mathcal{E}}\|_{\rm F}^{2}\) and a constant \((1-\frac{10R}{T})\). Since \(R\) is typically much smaller than \(T\) to promote model conciseness (i.e., \(R\ll T\)), the constant \((1-\frac{10R}{T})\) is approximately equal to \(1\). As a consequence, the lower bound is close to the noise power, which is significantly greater than zero. Therefore, **Proposition 1** points out that the one-layer TNN rejects fitting the noise with a high probability. Note that different from classical low-rank tensor decomposition works (Refs. [39, 40, 41] and references therein) that are primarily based on multi-linear forms, **Proposition 1** theoretically quantifies the noise-fitting ability of the proposed tensor model with the ReLU nonlinear activation function, which is novel and has a standalone value beyond SSF reconstruction tasks.

The key concept underlying **Proposition 1** is the intuition that a non-linear tensor decomposition model with intrinsically low-rank property, such as the one defined in Eq. (15), can effectively reject noise. While this proposition is specifically formulated for the _one

Figure 3: Illustration of the tensor neural network (TNN), which is composed of multiple TCLs. \(\mathbf{\mathcal{G}}\) is the core tensor and \(\mathbf{\mathcal{X}}\) is the output tensor.
layer_ TNN, the intuition behind noise rejection extends to TNN models with multiple layers. However, extending the results in **Proposition 1** to multiple-layer TNNs poses a daunting challenge. As a result, we experimentally validate its noise rejection property instead.

The multiple-layer TNN model used in this experiment has the same architecture as the one used in Sec. IV (see Fig. 7(c)), which has three layers with dimensions being (5, 5, 5), (10, 10, 10) and (20, 20, 20) respectively. In our experiment, we fit the proposed model to three distinct datasets: white Gaussian noise with standard deviation \(\sigma=0.5\), SSF data, and SSF data with added noise. In Fig. 4, we present the fitting error, quantified by the mean squared error (MSE), with respect to the iteration number of the parameter learning process. Fig. 4 demonstrates that the proposed model effectively fits the SSF data, as evidenced by the smallest MSE values across all iteration numbers. Conversely, the MSE values of the noise fitting are the highest, substantiating the model's noise rejection property. Furthermore, when the SSF data is contaminated by noise, the proposed model tends to fit the SSF data while fitting only a minor portion of the noise, resulting in the MSE curve (red curve) being much closer to that of fitting the noise-free SSF data (blue curve).

Subsequently, we establish a connection between the proposed TNN and the tensor Tucker model defined in Eq. (13). Specifically, if the activation function \(\varsigma(\cdot)\) is chosen to be the identity or linear function, following the definition of tensor Tucker decomposition,[39, 41] the \(L\)-layer TNN (defined in Eq. (15)) can be easily shown (see Appendix B) to be equivalent to the Tucker decomposition model (defined in Eq. (13)), that is,

\[\begin{split}\mathbf{\mathcal{X}}&=\varsigma\left( \cdots\varsigma\left(\mathbf{\mathcal{G}}\times_{1}\mathbf{W}_{1}^{(1)}\times_{2} \mathbf{W}_{1}^{(2)}\times_{3}\mathbf{W}_{1}^{(3)}\right)\cdots\right.\\ &\quad\times_{1}\mathbf{W}_{L}^{(1)}\times_{2}\mathbf{W}_{L}^{(2) }\times_{3}\mathbf{W}_{L}^{(3)}\right),\\ &=\mathbf{\tilde{\mathcal{G}}}\times_{1}\mathbf{\tilde{W}}^{(1)} \times_{2}\mathbf{\tilde{W}}^{(2)}\times_{3}\mathbf{\tilde{W}}^{(3)},\end{split} \tag{18}\]

where \(\varsigma(x)=ax\) with \(a\in\mathbb{R}\) being a constant, \(\mathbf{\tilde{\mathcal{G}}}=a^{L}\mathbf{\mathcal{G}}\) and \(\mathbf{\tilde{W}}^{(i)}=\mathbf{W}_{L}^{(i)}\mathbf{W}_{L-1}^{(i)}\cdots \mathbf{W}_{1}^{(i)},i=1,2,3\).

This result reveals that the tensor Tucker decomposition model, which is the cornerstone of the recently developed tensor-based SSF basis function learning framework that includes classical ones (such as EOFs and Fourier basis functions) as special cases, is _a particular instance_ of the proposed TNN model. In other words, _the proposed TNN model is the most general ocean SSF representation model to date_. The introduction of nonlinear activation functions significantly improves the tensor model's representation capability without compromising its concieness, allowing simultaneous decreases in both the representation error \(E_{1}\) and identification error \(E_{2}\). Specifically, the nonlinearity and hierarchical construction enable the effective capture of sound speed variations, leading to a small \(E_{1}\). Meanwhile, the Tucker decomposition-based backbone, or TCL, preserves the tensor model's conciseness, leading to a small \(E_{2}\).

After establishing a strong link between the proposed tensor neural network and the tensor Tucker decomposition model, one may contemplate extending the recoverability analysis[42] using Tucker models under either systematic sampling[43, 44] or random sampling[27]. However, due to the involvement of nonlinear activation functions, the analysis becomes a formidable challenge. Instead, in Sec. IV, we present a numerical analysis of the recoverability of the proposed model, while preserving its theoretical analysis as an intriguing avenue for future research.

### DNN-aided 3D SSF Reconstruction Algorithm

Based on the proposed TNN model, the 3D SSF reconstruction problem can be formulated as

\[\begin{split}\min_{\mathbf{\Theta}}&\quad\underbrace{ \left\|\mathbf{\mathcal{Y}}-\mathbf{\mathcal{O}}\ast\mathbf{\mathcal{X}}\right\|_{\rm F}^ {2}+\lambda R(\mathbf{\mathcal{X}})}_{\triangleq f(\mathbf{\Theta})},\\ \text{s.t.}&\quad\mathbf{\mathcal{X}}=\varsigma\left( \cdots\varsigma\left(\mathbf{\mathcal{G}}\times_{1}\mathbf{W}_{1}^{(1)}\times_{2} \mathbf{W}_{1}^{(2)}\times_{3}\mathbf{W}_{1}^{(3)}\right)\cdots\right.\\ &\quad\left.\times_{1}\mathbf{W}_{L}^{(1)}\times_{2}\mathbf{W}_{L}^ {(2)}\times_{3}\mathbf{W}_{L}^{(3)}\right),\end{split} \tag{20}\]

where \(R(\mathbf{\mathcal{X}})\) is a regularization term that allows the incorporation of side information for further performance enhancement, and \(\lambda\) is a hyper-parameter that balances the importance of the data fitting and regularization terms. Various types of structural information can be

Figure 4: The mean squared error (MSE) curves were obtained by fitting the proposed model to white Gaussian noise with standard deviation \(\sigma=0.5\), SSF data, and SSF data with added noise. The three-layer TNN model used in this experiment has the same architecture as the one used in Sec. IV (see Fig. 7(c)). The results demonstrate that the proposed model is inclined to fit the SSF data while effectively rejecting noise fitting.
incorporated by choosing different \(R(\mathbf{\mathcal{X}})\). If \(\lambda=0\), \(R(\mathbf{\mathcal{X}})\) is no longer effective. In the experimental results (Sec. IV.3), we demonstrate using the total variation (TV) regularizer to capture spatial correlations among sound speeds.

The model parameters \(\mathbf{\Theta}\) can be updated using the calculated gradients denoted as \(\nabla_{\mathbf{\Theta}}f(\mathbf{\Theta})\) via the gradient descent method [24]. This algorithm is widely used in deep learning and is straightforward. The gradients of the loss function with respect to the model parameters \(\mathbf{\Theta}\) can be obtained by utilizing the automatic differentiation technique in TensorFlow [45]. Appendix C provides some key steps of the derivation of the gradients. The update rule for the model parameters is given as follows:

\[\mathbf{\Theta}_{k+1}\leftarrow\mathbf{\Theta}_{k}-\eta_{k}\nabla_{\mathbf{\Theta}}f(\bm {\Theta}_{k}), \tag{21}\]

where \(\mathbf{\Theta}_{k}\) and \(\mathbf{\Theta}_{k+1}\) are the model parameters at iterations \(k\) and \(k+1\), respectively, and \(\eta_{k}\) is the learning rate. It's worth noting that the gradient descent method is a first-order method that can scale well to tackle big data [46]. Advanced acceleration schemes such as AdaGrad, RMSProp, and Adam can be used to further speed up the convergence [47]. We summarize the proposed algorithm in **Algorithm 1**. Note that the tensor \(\mathbf{\mathcal{G}}\) is a model parameter that is optimized during the reconstruction, while not the input of the reconstruction algorithm.

## IV Numerical results and discussions

In this section, numerical results using real-life ocean 3D SSF data are presented to showcase the encouraging reconstruction performance of the proposed tensor neural network-aided reconstruction method (labeled as TNN).

### Experimental Settings

**3D SSF Data**: The 3D South China Sea (SCS) data \(\mathbf{\mathcal{X}}\in\mathbb{R}^{20\times 20\times 20}\) on December 21, 2012, is considered in this paper and illustrated in Fig. 5. The data were derived from 3D conductivity, temperature, and depth (CTD) measurements using the hybrid coordinate ocean model (HYCOM). The spatial coverage of the dataset is 152 km \(\times\) 152 km \(\times\) 190 m, with a horizontal resolution of 8 km and a vertical resolution of 10 m.

**Sound Speed Measurements**: Following the configuration of the recent ocean observing systems (e.g., the Argo program [5] and Ocean Internet of Things [48]), the sound speed measurements are randomly collected over the 3D spatial region. The sampling ratio is defined as \(\rho=\frac{\sum_{i,j,k}\mathbf{\mathcal{O}}_{i,j,k}}{IJK}\). Examples of the sampled SSF under different sampling ratios in one single Monte-Carlo trial are shown in Fig. 6. In addition, the sampled data are assumed to be corrupted by i.i.d. Gaussian noise with zero mean and standard deviation \(\sigma\).

**Performance measure**: The reconstruction performances of different methods are assessed by the root mean squared error (RMSE), given by

\[\text{RMSE}=\sqrt{\frac{1}{T}\|\mathbf{\mathcal{\hat{X}}}-\mathbf{\mathcal{X}}\|_{ \text{F}}^{2}}, \tag{22}\]

where \(\mathbf{\mathcal{X}}\) is the ground-truth; \(\mathbf{\mathcal{\hat{X}}}\) is the reconstructed SSF; and \(T=IJK\) denotes the total number of SSF data entries. The RMSEs were averaged over three Monte-Carlo trials with different sampling patterns. All the experiments are conducted in a computer with a 2.2GHz 6-Core Intel i7 CPU.

**Training Optimizer**: The gradient computation was carried out using the automatic differentiation mechanism of TensorFlow. To optimize the proposed TNN model, we use the popular Adam optimizer with an initial learning rate of 0.005 and fix the number of iterations at 15,000.

### Validating Reconstruction Error Analysis

To validate the reconstruction error analysis results presented in Sec. II, numerical results are first presented in this subsection before comparing the proposed algorithm with SOTA methods. Specifically, we compare the reconstruction performance of the proposed tensor neural network-aided algorithm with the methods that use only deep learning models or tensor models. The results demonstrate the importance of striking the right balance between model conciseness and expressiveness for achieving superior reconstruction performance. In this subsection, no extra prior information is included in the TNN for a fair comparison, i.e., the regularization parameter \(\lambda\) is set to zero. And unless otherwise specified, \(\sigma\) is set to 0.1, indicating that the RMSEs of sound speed measurements/estimates mostly lie in the range of 0.1-0.3 m/s, a moderate range according to our real-life sea trials. The reconstruction performances under a wider range of noise powers are reported in the next subsection.

## 1 Conciseness Matters

First, the importance of conciseness in reconstruction is demonstrated by comparing the reconstruction performances of different neural network-based methods. The baseline networks used for this comparison were a
three-layer multi-layer perceptron (MLP) and a three-layer convolutional neural network (CNN). In addition, a three-layer tensor neural network (TNN) was used for the proposed method. The detailed network architectures are illustrated in Fig. 7. The reconstruction performances of these network-based methods under different sampling ratios are presented in Table 1. For visualization, Fig. 8 presents the recovered SSF of different methods under a sampling ratio of 0.3 in one Monte-Carlo trial.

The results presented in Table 1 and Fig. 8 indicate that the SSF reconstruction performance of the MLP is poor. This is due to the large number of parameters in the MLP (around 8 million), which make it difficult to learn the optimal model from limited and noisy measurements (e.g., 2400), resulting in a significant identification error \(E_{2}\). In contrast, the CNN, which has fewer parameters (1565), can leverage the spatial smoothness of the SSF and exhibit better reconstruction performance. However, the CNN does not exploit the 3D spatial correlations among sound speeds as well as the TNN. Specifically, the proposed TNN model, utilizing tensor computation, can effectively capture the 3D spatial correlations with fewer parameters (875), resulting in significantly better reconstruction performance.

Additionally, we compare the reconstruction performance of TNN and CNN with 3D convolutions. Fig. 9 shows the architecture of the three-layer CNN model with 3D convolutions, while the three-layer TNN model has the same neural network architectures as the one in Fig. 7(c). Here we denote the CNN models with \(C=10\) and \(C=20\) as CNN1 and CNN2 respectively, where \(C\) is the number of kernels in the hidden layer. Table 2 presents the averaged reconstruction RMSEs of different models under various sampling ratios. Our findings are as follows:

Figure 5: Illustration of the 3D SSF data.

Figure 6: Measurements of SSF under different sampling ratios in one single Monte-Carlo trial.
1) TNN achieves the best performance in all scenarios. While CNN models are capable of exploiting spatial correlations, they do not fully leverage the intrinsic low-rankness of the SSF. On the other hand, the TNN model, which nonlinearly concatenates several low-rank tensor models, better exploits the global coherence of the SSF. Notably, even with the minimum number of model parameters, the TNN model achieves the best reconstruction performance among the three models when the SSF is fully observed (\(\rho=1\)). This result indicates that the TNN model has a lower representation error because the identification error is negligible when \(\rho=1\), demonstrating the strong expressiveness of the TNN model for SSF representation.

2) CNN2 outperforms CNN1 in all scenarios. CNN2 has more parameters and higher expressive power, resulting in a lower representation error. Surprisingly, CNN2 also yields better reconstruction results under sparse sampling conditions (e.g., \(\rho=0.3\)). This could be attributed to the inductive bias of the network architecture, allowing it to produce reasonable outcomes even for highly ill-posed inverse problems, as evidenced by recent works in the deep learning literature (see Ref. [49]).

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline \hline \(\rho\) & \(N\) & **MLP**(\(\approx\) 8M) & **CNN**(1565) & **TNN**(875) \\ \hline
0.2 & 1600 & 5.67 & 0.63 & **0.50** \\ \hline
0.3 & 2400 & 5.37 & 0.59 & **0.34** \\ \hline
0.4 & 3200 & 4.94 & 0.56 & **0.28** \\ \hline \hline \end{tabular}
\end{table}
Table 1: RMSEs of different neural network-aided methods under different sampling ratio \(\rho\). The number of measurements \(N\) is presented in the second column. And the number of parameters for different models is in parentheses.

Figure 8: Visual effects of reconstructed SSF of different neural network-based methods at depths 0 m, 90 m, and 190 m under a sampling ratio of 0.3 in one single Monte-Carlo trial. The RMSEs are shown at the top.

Figure 7: Th network architecture of (a) MLP, (b) CNN, (c) TNN. The number of channels C is set to 10. The MLP is composed of three layers with dimensions of 125, 1000, and 8000, respectively. The convolutional layers have a kernel size of (2,2), and the Upsampling layers use the nearest-neighbor interpolation method. For all neural networks, the ReLU activation function is chosen for all layers except for the final layer, which uses the tanh activation function.
## 2 Expressiveness Matters

Then, experiments were conducted to demonstrate the effectiveness of enhancing the model's representation capability (expressiveness) in improving reconstruction performance. We compare the proposed model with two tucker models with different architectures. The first model, named Tucker1, has the same structure as the TNN model but without non-linearity (i.e., using linear activation functions). The second model, named Tucker2, is just a tensor Tucker decomposition model (defined in Eq. (13)) with the dimension of the core tensor setting to \((7,8,8)\).

The number of parameters in Tucker2 is \(20\times(7+8+8)+7\times 8\times 8=908\), while TNN and Tucker1 both have \((20\times 10+10\times 5)\times 3+5\times 5\times 5=875\) parameters. Therefore, the three models have a similar number of parameters. Table 3 presents the average RMSEs of the three models under different sampling ratios. Visualization of the reconstructed SSFs in one Monte-Carlo trial under \(\rho=0.3\) is shown in Fig. 10.

Based on the experimental results shown in Table 3 and Fig. 10, we can draw the following conclusions. First, when the number of measurements is relatively large (e.g., \(N=2400\)), Tucker1 has the highest RMSE (poorest reconstruction) due to its limited degree of freedom. According to Eq. (18), Tucker1 is mathematically equivalent to a Tucker decomposition model (defined in Eq. (13)) with a core tensor of dimension (5, 5, 5), which results in a larger representation error \(E_{1}\) compared to the other two models. Second, although Tucker2 has slightly more parameters than TNN, it still has a larger RMSE due to its limited representation capability imposed by the multi-linear form. This demonstrates the importance of introducing non-linear activation functions in the model. Finally, the proposed TNN model achieves the best performance when the number of measurements is not too small, and shows a similar RMSE to Tucker1 under a small sampling ratio (e.g., \(\rho=0.2\)), thanks to its high expressive power and concise structure.

## 3 Recoverability, Trade-offs between being "Deeper" or "Wider"

To evaluate the recoverability and draw more insights into the proposed model, we compare the reconstruction performance of TNN models with different configurations (see Table 4) assuming no observation noise. Fig. 11 shows the RMSEs under three different sampling ratios. More experimental results can be found in Appendix D.

It can be seen from Fig. 11 that the recovery/reconstruction RMSEs of the proposed TNN models (with different layer numbers/configurations) keep decreasing to a very small number (e.g., \(<0.1\)) as the sam

\begin{table}
\begin{tabular}{|c|c|c|c|} \hline \hline \(\rho\) & TNN(875) & CNN1(1106) & CNN2(3686) \\ \hline
0.3 & **0.35** & 0.57 & 0.46 \\ \hline
0.5 & **0.23** & 0.51 & 0.40 \\ \hline
1 & **0.21** & 0.48 & 0.37 \\ \hline \hline \end{tabular}
\end{table}
Table 2: The averaged RMSEs of TNN and 3D CNN models under different sampling ratios. The number of model parameters is in parentheses.

Figure 10: Visual effects of reconstructed SSF of different tensor-based models at depths 0 m, 90 m, and 190 m under a sampling ratio of 0.3 in one single Monte-Carlo trial. The RMSEs are shown at the top.

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline \hline \(\rho\) & \(N\) & **Tucker1**(875) & **Tucker2**(908) & **TNN**(875) \\ \hline
0.2 & 1600 & **0.47** & 1.20 & 0.50 \\ \hline
0.3 & 2400 & 0.43 & 0.70 & **0.34** \\ \hline
0.4 & 3200 & 0.42 & 0.34 & **0.28** \\ \hline \hline \end{tabular}
\end{table}
Table 3: RMSEs of different tensor-based methods under different sampling ratios. The number of measurements \(N\) is presented in the second column. And the number of parameters for different models is in parentheses.

Figure 9: The architecture of the 3D CNN model, where \(C\) is the number of kernels in the hidden layer.
pling ratio increases (i.e., with more samples). This suggests that the recoverability of the proposed model is highly likely to be guaranteed when samples are abundant and clean.

On the other hand, despite having the largest number of layers, TNN1 gives the highest error when the \(\rho=1\) (i.e., the SSF is fully observed). The reason is that TNN1 has the smallest number of model parameters among the three models, as shown in Table 4. Generally, making the TNN "deeper" can significantly reduce the number of model parameters at the cost of slightly increasing the representation error. On the other hand, making the TNN "wider" can effectively increase expressiveness while also increasing the identification error, as seen in TNN3 in Table 4 and Fig. 11. The main idea behind the TNN model is to hierarchically reduce the tensor dimensionality through non-linear approximation, allowing for effective parameter reduction without significantly hampering expressiveness.

The experimental results indicate that designing the TNN model involves trade-offs between being wider or deeper. Although it is challenging to determine the optimal sweet spot theoretically, we can offer some practical suggestions on how to strike a good balance in practice:

* It is advisable to keep the number of unknown model parameters less than the number of samples in order to potentially prevent ill-posed reconstruction problems.
* Consider increasing the dimensionality of the core tensor layer by layer when the number of layers is fixed. For example, the dimensionality of the \(l\)th layer can be set to twice that of the \((l-1)\)th layer. This technique has shown promise in reducing the number of parameters while maintaining expressiveness.
* Recent studies have suggested that a one-layer Tucker model can effectively capture the main pattern of sound speeds, and adding just two or three more layers may be sufficient to capture non-linear variations and reduce the number of parameters. Consequently, the layer number of TNN is recommended to be set to 3 or 4.

It should be noted that these configurations were tested under the experimental settings outlined in Sec. IV.1, and the optimal model configuration may vary depending on the specific application or scale.

### Enhancement by TV Regularization

In this subsection, experiments were conducted to demonstrate how the reconstruction performance of TNN can be further boosted by incorporating the structural information of the regularizer. For illustration purposes, the widely-used total variation (TV) [25] regularizer was used to enforce spatial smoothness among sound speeds. Specifically, the TV regularizer is defined as follows:

\[\begin{split} R(\mathbf{\mathcal{X}})&=\sum_{i,j,k}(| \mathbf{\mathcal{X}}(i+1,j,k)-\mathbf{\mathcal{X}}(i,j,k)|+|\mathbf{\mathcal{X}}(i,j+1,k) \\ &-\mathbf{\mathcal{X}}(i,j,k)|+|\mathbf{\mathcal{X}}(i,j,k+1)-\mathbf{ \mathcal{X}}(i,j,k)|).\end{split} \tag{23}\]

The reconstructed results of the proposed TNN model with TV regularizer, denoted as TNN-TV, are presented in Fig. 12. The hyper-parameter \(\lambda\) was manually tuned to optimize performance.

The results demonstrate that TNN-TV consistently outperforms TNN in SSF reconstruction, particularly under very sparse sampling (e.g., \(\rho=0.1\)). Note that when the sampling ratio \(\rho\geq 0.3\), the measurements can provide sufficient information for learning the model parameters, and hence TNN and TNN-TV have similar performances. In this case, \(\lambda\) can be set to be nearly zero, and TNN-TV gradually reduces to TNN (see Eq. (19)). Since TNN is a special case of TNN-TV, in the next subsection, we only consider TNN-TV and refer to it as TNN for simplicity.

### Comparisons with SOTA methods

Previous subsections have demonstrated the conciseness and expressiveness of the proposed model. In this

\begin{table}
\begin{tabular}{|c|c|c|c|} \hline \hline Name & L & Dimensionality & M \\ \hline TNN1 & 3 & (5, 5, 5), (10, 10, 10), (20, 20, 20) & 875 \\ \hline TNN2 & 2 & (10, 10, 10), (20, 20, 20) & 1600 \\ \hline TNN3 & 2 & (15, 15, 15), (20, 20, 20) & 4275 \\ \hline \hline \end{tabular}
\end{table}
Table 4: The configurations of TNN models. M is the number of parameters.

Figure 11: The averaged RMSEs of TNN models with different configurations under different sampling ratios.
subsection, we compare the reconstruction performance of the proposed algorithm with several SOTA SSF reconstruction methods.

**SOTA Methods**: The selected SOTA methods include: 1) low-rank matrix-based method represented by the recently proposed graph-guided Bayesian matrix completion (BMCG)[8]; 2) tensor-based methods represented by alternating least squares for the Tucker model (Tucker-ALS), low-rank tensor completion (LRTC)[50], and LRTC with total variation (LRTC-TV)[51]; and 3) statistical learning-based methods represented by Gaussian process regression (GPR)[14].

**Model Settings**: The kernel of GPR is the widely used radial basis function (RBF), with the hyper-parameters being learned via evidence optimization[18]. In BMCG, the graph Laplacian matrix is constructed using the same kernel[8]. The tensor rank surrogate function used in LRTC is the tensor nuclear norm, defined as \(\sum_{k=1}^{K}\|\widehat{\mathcal{X}}(:,:,k)\|_{*}\)[50], where \(\widehat{\mathcal{X}}\) is the Fourier transformation of \(\mathcal{X}\) along the third mode.

**Results and Discussions**: The reconstruction RMSEs of different methods under different sampling ratios and noise powers are shown in Table 5. Visual inspections of reconstructed SSFs and error surfaces under different scenarios are presented in Figs. 13-16. As can be seen, the proposed tensor neural network-aided reconstruction algorithm outperforms all other methods in terms of reconstruction accuracy, regardless of the sampling ratio and noise power. The reason is that the proposed model has a high expressive power with relatively few parameters, achieving an outstanding trade-off between \(E_{1}\) and \(E_{2}\), as discussed in Sec. III. Additionally, the TNN model exhibits strong noise rejection (as proved in **Proposition 1**), making it robust against measurement noise. In the following, we provide additional discussions to draw further insights from the results.
* The matrix-based method, BMCG, yields good reconstruction results for the horizontal SSF slices by explicitly exploring the spatial smoothness through graph models, as discussed in Sec. III of Ref. [8]. However, the vertical correlations of sound speeds are ignored in this method, leading to the loss of fine-grained details in the vertical dimension. This can be observed in the error surfaces of the vertical dimension shown in Fig. 17.
* Among tensor completion methods (i.e., Tucker-ALS, LRTC, and LRTC-TV), Tucker-ALS gives the best reconstruction performance. It outperforms the matrix-based method BMCG when the observation is moderate (e.g., \(\rho=0.3\)) since it can effectively exploit the multi-dimensional correlations. However, the multi-linear form has limited the expressive power of Tucker-ALS, resulting in a considerable \(E_{1}\), making the associated reconstruction performance inferior to the proposed one. Moreover, these tensor completion methods are sensitive to the sampling pattern and exhibit significantly degraded performance when the measurements are sparse.
* GPR is a non-parametric statistical model that paves another promising path for reconstructing SSFs. The expressiveness of the kernel determines the representation error \(E_{1}\), while the learning process of the kernel hyper-parameters determines the identification error \(E_{2}\). Although GPR has demonstrated great success in many areas in recent years[38, 52, 3], research on its applications in ocean SSF reconstruction is still in its infancy (including optimal kernel design and hyper-parameter learning). In our experiments, GPR, using the RBF kernel and evidence maximization-based hyper-parameter learning, achieved the second-best performance in almost all scenarios. However, to enhance its performance, more powerful kernels and advanced optimization techniques can be employed[53, 53], which require significant research efforts and are beyond the scope of this paper. Moreover, the high computational complexity of GPR limits its widespread application. The scalability of GPR is an active research area in machine learning, as evidenced by recent works[54, 55].

## V Conclusions and Future Directions

In this paper, under the unified framework of analyzing the reconstruction error, a tensor neural network-aided reconstruction algorithm is proposed. By leveraging the succinct form of tensor models and the high expressive power of deep neural networks, the proposed TNN model can concisely and accurately represent the 3D SSF, achieving an outstanding balance between the representation error \(E_{1}\) and identification error \(E_{2}\). A simple and efficient algorithm is developed based on the

Figure 12: The RMSEs of TNN-TV and TNN under different sampling ratios.
widely adopted gradient descent method. Experiments using real-life 3D SSF data demonstrate the effectiveness of the proposed reconstruction error analysis framework and showcase the encouraging performance of the proposed algorithm compared to other SOTA reconstruction methods. Our study demonstrated how the proposed model could effectively integrate structural assumptions, as exemplified by the use of the TV regularizer.

Future research could explore the incorporation of additional prior information, such as the physical knowledge of ocean processes that drive sound speed variations, to achieve even better reconstruction performance. It would also be valuable to extend the analysis in Ref. [12] to quantify the recoverability of the proposed TNN model and establish the theoretical connection between TNN and CNN with 3D convolutions. Finally, rethinking the reconstruction problem from a Bayesian learning perspective [56, 57] would hold significant research value. This approach allows for automatic control of model complexity to prevent both overfitting and underfitting, and enables the utilization of uncertainty information to guide future sensing planning.

\begin{table}
\begin{tabular}{|c|c|c c c c c|} \hline \hline \(\rho\) & \(\sigma\) & **BMCG** & **GPR** & **Tucker-ALS** & **LRTC** & **LRTC-TV** & **TNN** \\ \hline  & 0.1 & 1.33 & 0.68 & 1.32 & 1.47 & 1.23 & **0.47** \\
0.1 & 0.3 & 1.34 & 0.71 & 1.36 & 1.50 & 1.25 & **0.51** \\  & 0.5 & 1.36 & 0.77 & 1.45 & 1.56 & 1.30 & **0.58** \\ \hline  & 0.1 & 0.72 & 0.52 & 0.76 & 0.96 & 0.77 & **0.35** \\
0.2 & 0.3 & 0.74 & 0.56 & 0.80 & 1.02 & 0.82 & **0.39** \\  & 0.5 & 0.77 & 0.61 & 0.84 & 1.12 & 0.89 & **0.46** \\ \hline  & 0.1 & 0.56 & 0.44 & 0.42 & 0.64 & 0.53 & **0.30** \\
0.3 & 0.3 & 0.59 & 0.48 & 0.44 & 0.73 & 0.60 & **0.33** \\  & 0.5 & 0.64 & 0.53 & 0.47 & 0.85 & 0.69 & **0.40** \\ \hline  & 0.1 & 0.48 & 0.37 & 0.41 & 0.46 & 0.39 & **0.28** \\
0.4 & 0.3 & 0.52 & 0.44 & 0.42 & 0.57 & 0.48 & **0.31** \\  & 0.5 & 0.57 & 0.49 & 0.45 & 0.71 & 0.59 & **0.34** \\ \hline \hline \end{tabular}
\end{table}
Table 5: The average RMSEs over three Monte-Carlo trials of different algorithms under different sampling ratios and noise powers.

Figure 13: Visual effects of the reconstructed SSF of different SOTA reconstruction methods in one single Monte-Carlo trial with sampling ratio \(\rho=0.2\) and \(\sigma=0.1\). The RMSEs are shown above the top subfigures.
## VI Acknowledgement

This work was supported in part by the National Key R&D Program of China, in part by the National Natural Science Foundation of China under Grant 62001309, in part by the Fundamental Research Funds for the Central Universities, in part by the Zhejiang University Education Foundation Qizhen Scholar Foundation, and in part by Science and Technology on Sonar Laboratory under Grant 6142109KF212204.

## Appendix A Tensor operation definitions

**Definition 1** (mode-n product): _The mode-\(n\) product of a tensor \(\mathbf{\mathcal{A}}\in\mathbb{R}^{I_{1}\times\cdots\times I_{N}}\) and a matrix \(\mathbf{B}\in\mathbb{R}^{J\times I_{n}}\) is defined by_

\[\mathbf{\mathcal{C}}=\mathbf{\mathcal{A}}\times_{n}\mathbf{B}\in\mathbb{R}^{I_{1} \times\cdots\times I_{n-1}\times J\times I_{n+1}\times\cdots\times I_{N}}, \tag{11}\]

_whose entries are given by_

\[\mathbf{\mathcal{C}}(i_{1},\cdots,j,\cdots,i_{N})=\sum_{i_{n}=1}^{I_{n}}\mathbf{ \mathcal{A}}(i_{1},\cdots,i_{n},\cdots,i_{N})\mathbf{B}(j,i_{n}). \tag{12}\]

**Definition 2** (mode-\(n\) unfolding): _Given a tensor \(\mathbf{\mathcal{A}}\in\mathbb{R}^{I_{1}\times\cdots\times I_{N}}\), its mode-\(n\) unfolding gives a matrix \(\mathbf{A}_{(n)}\in\mathbb{R}^{I_{1}\times\cdots\times I_{N}}\)._

Figure 14: Visual effects of the reconstructed SSF of different SOTA reconstruction methods in one single Monte-Carlo trial with sampling ratio \(\rho=0.3\) and \(\sigma=0.1\). The RMSEs are shown above the top subfigures.

Figure 15: Visual effects of the reconstructed SSF of different SOTA reconstruction methods in one single Monte-Carlo trial with sampling ratio \(\rho=0.3\) and \(\sigma=0.5\). The RMSEs are shown above the top subfigures.
\(\mathbb{R}^{I_{n}\times\prod_{i=1,k\neq n}^{N}I_{k}}\). Each tensor element \(\mathbf{\mathcal{A}}(i_{1},\cdots,i_{n})\) is mapped to the matrix element \([\mathbf{A}_{(n)}](i_{n},j)\) where \(j=1+\sum_{k=1,k\neq n}^{N}(i_{k}-1)J_{k}\) with \(J_{k}=\prod_{m=1,m\neq n}^{k-1}I_{m}\).

**Definition 3** (Tucker decomposition): _For a tensor \(\mathbf{\mathcal{A}}\in\mathbb{C}^{I_{1}\times\cdots\times I_{N}}\), the Tucker decomposition is defined as_

\[\mathbf{\mathcal{A}}=\mathbf{\mathcal{G}}\times_{1}\mathbf{U}^{(1)}\times_{2}\mathbf{ U}^{(2)}\times_{3}\cdots\times_{N}\mathbf{U}^{(N)}, \tag{10}\]

_where each factor matrix \(\mathbf{U}^{(n)}\in\mathbb{C}^{I_{n}\times R_{n}},\forall n=1,\cdots,N\). The core tensor \(\mathbf{\mathcal{G}}\in\mathbb{C}^{R_{1}\times\cdots\times R_{N}}\). The tuple \((R_{1},\cdots,R_{N})\) is known as the multi-linear rank._

**Definition 4** (tensor inner product): _The inner product of tensor \(\mathbf{\mathcal{A}}\) and \(\mathbf{\mathcal{B}}\) with the same size is formulated as follows:_

\[c=\langle\mathbf{\mathcal{A}},\mathbf{\mathcal{B}}\rangle=\langle\mathrm{vec}(\mathbf{ \mathcal{A}}),\mathrm{vec}(\mathbf{\mathcal{B}})\rangle\in\mathbb{R}, \tag{11}\]

_where \(\mathrm{vec}(\mathbf{\mathcal{A}})\) is the vectorization of \(\mathbf{\mathcal{A}}\)._

## Appendix B Proof and derivation

### Derivation of Eq. (4)

According to its definition, the reconstruction error \(E\) can be written as

\[E =\|\mathbf{\mathcal{X}}-\hat{\mathbf{\mathcal{X}}}\|_{\mathrm{F}}^{2} \tag{12}\] \[=\|\mathbf{\mathcal{X}}-\mathbf{\mathcal{D}}(\hat{\mathbf{\Theta}})\|_{ \mathrm{F}}^{2}\] \[=\|\mathbf{\mathcal{X}}-\mathbf{\mathcal{D}}(\mathbf{\Theta}^{*})+\mathbf{ \mathcal{D}}(\mathbf{\Theta}^{*})-\mathbf{\mathcal{D}}(\hat{\mathbf{\Theta}})\|_{ \mathrm{F}}^{2}\] \[=\underbrace{\|\mathbf{\mathcal{X}}-\mathbf{\mathcal{D}}(\mathbf{\Theta}^{*}) \|_{\mathrm{F}}^{2}}_{\triangleq\ E_{1}}+\underbrace{\|\mathbf{\mathcal{D}}(\mathbf{ \Theta}^{*})-\mathbf{\mathcal{D}}(\hat{\mathbf{\Theta}})\|_{\mathrm{F}}^{2}}_{ \triangleq\ E_{2}}+\] \[\underbrace{2(\mathbf{\mathcal{X}}-\mathbf{\mathcal{D}}(\mathbf{\Theta}^{*}),\mathbf{\mathcal{D}}(\mathbf{\Theta}^{*})-\mathbf{\mathcal{D}}(\hat{\mathbf{\Theta}}))}_{ \triangleq\ \epsilon}.\]

Thus the reconstruction error can be decomposed as \(E=E_{1}+E_{2}+\epsilon\) and Eq. (4) holds.

### Derivation of Eq. (5)

If \(\mathbf{\mathcal{X}}=\mathbf{\mathcal{D}}(\mathbf{\Theta}^{*})\), we have \(E_{1}=\|\mathbf{\mathcal{X}}-\mathbf{\mathcal{D}}(\mathbf{\Theta}^{*})\|_{\mathrm{F}}^{2}=0\) and \(\epsilon=2\langle\mathbf{0},\mathbf{\mathcal{D}}(\mathbf{\Theta}^{*})-\mathbf{\mathcal{D}}( \hat{\mathbf{\Theta}})\rangle=0\). Similarly, if \(\mathbf{\mathcal{D}}(\mathbf{\Theta}^{*})=\mathbf{\mathcal{D}}(\hat{\mathbf{\Theta}})\), then \(E_{2}=\|\mathbf{\mathcal{D}}(\mathbf{\Theta}^{*})-\mathbf{\mathcal{D}}(\hat{\mathbf{\Theta}}) \|_{\mathrm{F}}^{2}=0\) and \(\epsilon=2\langle\mathbf{\mathcal{X}}-\mathbf{\mathcal{D}}(\mathbf{\Theta}^{*}),\mathbf{0} \rangle=0\).

Figure 16: Visual effects of the error surfaces in one single Monte-Carlo trial of different methods under \(\rho=0.2\) and \(\sigma=0.1\). The RMSEs are shown above the top subfigures.

Figure 17: The error surfaces of tensor-based method TNN and matrix-based method BMCG in the vertical dimension in one Monte-Carlo trial. The TNN can exploit multi-dimensional correlations thus having smoother error surfaces.
### Proof of Proposition 1

We start by rewriting the term \(\|\mathbf{\mathcal{X}}-\mathbf{\mathcal{E}}\|_{\mathrm{F}}^{2}\) in a convenient form. Specifically, since each element of \(\mathbf{\mathcal{E}}\) follows a zero-mean i.i.d Gaussian distribution, we have \(\|\mathbf{\mathcal{X}}-\mathbf{\mathcal{E}}\|_{\mathrm{F}}^{2}=\|\mathrm{vec}(\mathbf{ \mathcal{X}})-\zeta\|_{\mathrm{F}}^{2}\), where \(\mathrm{vec}(\mathbf{\mathcal{X}})\in\mathbb{R}^{T}\) is the vectorization of \(\mathbf{\mathcal{X}}\) and \(\zeta\sim\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I})\) where \(\mathbf{I}\) is a \(T\times T\) identity matrix. Consider a set of parameters \(\mathbf{\Theta}=\{\mathbf{\mathcal{G}},\mathbf{W}_{1}^{(1)},\mathbf{W}_{1}^{(2)}, \mathbf{W}_{1}^{(3)}\}\), denote \(\mathbf{\mathcal{O}}\) as the indicator tensor that contains one if \(\mathbf{\mathcal{X}}_{i,j,k}>0\) and zero otherwise. Then we have

\[\mathrm{vec}(\mathbf{\mathcal{X}}) =\mathrm{vec}\left(\mathrm{ReLU}(\mathbf{\mathcal{G}}\times_{1} \mathbf{W}_{1}^{(1)}\times_{2}\mathbf{W}_{1}^{(2)}\times_{3}\mathbf{W}_{1}^{( 3)})\right) \tag{14}\] \[=\mathrm{vec}\left(\mathbf{\mathcal{O}}*(\mathbf{\mathcal{G}}\times_{1} \mathbf{W}_{1}^{(1)}\times_{2}\mathbf{W}_{1}^{(2)}\times_{3}\mathbf{W}_{1}^{( 3)})\right)\] \[=\underbrace{\left[\mathrm{diag}(\mathrm{vec}(\mathbf{\mathcal{O}}))( \mathbf{W}_{1}^{(3)}\otimes\mathbf{W}_{1}^{(2)}\otimes\mathbf{W}_{1}^{(1)}) \right]}_{\triangleq\mathbf{B}\in\mathbb{R}^{T\times R}}\mathrm{vec}(\mathbf{ \mathcal{G}})\] \[=\mathbf{B}\mathrm{vec}(\mathbf{\mathcal{G}}).\]

Therefore, \(\mathrm{vec}(\mathbf{\mathcal{X}})\) is a weighted combination of \(R\)\(T\)-dimensional vectors and lies in a at-most-\(R\)-dimensional subspace \(S\) of \(\mathbb{R}^{T}\). It follows that

\[\min_{\mathbf{\Theta}}\|\mathbf{\mathcal{X}}-\mathbf{\mathcal{E}}\|_{\mathrm{F}}^{2}\geq\|P _{S^{c}\zeta}\|_{2}^{2}, \tag{15}\]

where \(P_{S^{c}\zeta}\) is the projection of \(\zeta\) outside the subspace \(S\).

Next, we make use of the following lemma to give a bound on the projection of the noise \(\zeta\) onto a subspace.

**Lemma 1**: Let \(S\subset\mathbb{R}^{T}\) be a subspace with dimension \(R\). Let \(\zeta\sim\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I})\) and \(\beta\geq 1\). Then,

\[P\left[\frac{\|P_{S^{c}\zeta}\|_{2}^{2}}{\|\zeta\|_{2}^{2}}\geq 1-\frac{10 \beta R}{T}\right]\geq 1-e^{-\beta R}-e^{-T/16}. \tag{16}\]

Proof.: From the Lemma 1 in Ref. [58], if \(X\sim\chi_{T}^{2}\), then

\[P[X-T\geq 2\sqrt{Tx}+2x] \leq e^{-x}, \tag{17}\] \[P[X\leq T-2\sqrt{Tx}] \leq e^{-x}. \tag{18}\]

We have \(\frac{\|P_{S^{c}\zeta}\|_{2}^{2}}{\|\zeta\|_{2}^{2}}=1-\frac{\|P_{S^{c}}\|_{2} ^{2}}{\|\zeta\|_{2}^{2}}\). Note that \(\|P_{S}\zeta\|_{2}\sim\chi_{R}^{2}\) and \(\|\zeta\|_{2}^{2}\sim\chi_{T}^{2}\). Applying the inequality (17) to bound \(\|P_{S}\zeta\|_{2}\) and inequality (18) to \(\|\zeta\|_{\mathrm{F}}^{2}\), a union bound gives that claim.

With Lemma 1 and inequality (15), we have

\[P\left[\frac{1}{\|\mathbf{\mathcal{E}}\|_{\mathrm{F}}^{2}}(\min_{\mathbf{\Theta}}\|\bm {\mathcal{X}}-\mathbf{\mathcal{E}}\|_{\mathrm{F}}^{2})\geq 1-\frac{10\beta R}{T} \right]\geq 1-e^{-\beta R}-e^{-T/16}, \tag{19}\]

for \(\beta\geq 1\). Let \(\beta=1\), then

\[P\left[\min_{\mathbf{\Theta}}\|\mathbf{\mathcal{X}}-\mathbf{\mathcal{E}}\|_{\mathrm{F}}^{2 }\geq\|\mathbf{\mathcal{E}}\|_{\mathrm{F}}^{2}(1-\frac{10R}{T})\right]\geq 1-e^{-R}-e^{-T/16}. \tag{20}\]

Thus Proposition 1 holds.

### Derivation of Eq. (18)

Before giving the derivation of Eq. (18), we present a few facts regarding \(n\)-mode matrix products. Specifically, for distinct modes in a series of multiplications, the order of the multiplication is irrelevant, i.e.,

\[\mathbf{\mathcal{X}}\times_{m}\mathbf{A}\times_{n}\mathbf{B}=\mathbf{\mathcal{X}} \times_{n}\mathbf{B}\times_{m}\mathbf{A}\ \ (m\neq n). \tag{21}\]

If the modes are the same, then

\[\mathbf{\mathcal{X}}\times_{m}\mathbf{A}\times_{m}\mathbf{B}=\mathbf{\mathcal{X}} \times_{m}(\mathbf{B}\mathbf{A}). \tag{22}\]

Given these facts, a \(L\)-layer TNN with a linear activation function can be written as

\[\mathbf{\mathcal{X}} =\varsigma(\cdots\varsigma\left(\mathbf{\mathcal{G}}\times_{1}\mathbf{W }_{1}^{(1)}\times_{2}\mathbf{W}_{1}^{(2)}\times_{3}\mathbf{W}_{1}^{(3)}\right)\cdots \tag{23}\] \[\quad\times_{1}\mathbf{W}_{L}^{(1)}\times_{2}\mathbf{W}_{L}^{(2)} \times_{3}\mathbf{W}_{L}^{(3)})\] \[=a(\cdots a\left(\mathbf{\mathcal{G}}\times_{1}\mathbf{W}_{1}^{(1)} \times_{2}\mathbf{W}_{1}^{(2)}\times_{3}\mathbf{W}_{1}^{(3)}\right)\cdots\] \[\quad\times_{1}\mathbf{W}_{L}^{(1)}\times_{2}\mathbf{W}_{L}^{(2)} \times_{3}\mathbf{W}_{L}^{(3)})\] \[=a^{L}\mathbf{\mathcal{G}}\times_{1}\mathbf{W}_{1}^{(1)}\cdots\times_{ 1}\mathbf{W}_{L}^{(1)}\times_{2}\mathbf{W}_{1}^{(2)}\cdots\] \[\quad\times_{2}\mathbf{W}_{L}^{(2)}\times_{3}\mathbf{W}_{1}^{(3)} \cdots\times_{3}\mathbf{W}_{1}^{(3)}\] \[=a^{L}\mathbf{\mathcal{G}}\times_{1}(\mathbf{W}_{L}^{(1)}\cdots\mathbf{ \mathcal{W}}_{1}^{(1)})\times_{2}(\mathbf{W}_{L}^{(2)}\cdots\mathbf{W}_{1}^{(2)})\] \[\quad\times_{3}(\mathbf{W}_{L}^{(3)}\cdots\mathbf{\mathcal{W}}_{1}^{(3)})\] \[=\mathbf{\mathcal{G}}\times_{1}\tilde{\mathbf{W}}^{(1)}\times_{2} \tilde{\mathbf{W}}^{(2)}\times_{3}\tilde{\mathbf{W}}^{(3)}.\]

## Appendix C Derivation of the gradient

The gradient with respect to the model parameters can be obtained via the backpropagation technique.[59] Since the gradient of the activation function depends on the specific choice, here we mainly introduce the gradient propagation in the TCL. By unfolding the tensor to the matrix, we can calculate the gradients with respect to the factor matrices, given by

\[\frac{\partial\mathbf{X}_{(k)}^{l+1}}{\partial\mathbf{W}_{l}^{(k)}}=\frac{ \partial\mathbf{W}_{l}^{(k)}\mathbf{X}_{(k)}^{l}(\mathbf{W}_{l}^{(-k)})^{ \mathrm{T}}}{\partial\mathbf{W}_{l}^{(k)}} \tag{24}\]

where \(\mathbf{W}_{l}^{(-k)}=\mathbf{W}_{l}^{(3)}\otimes\cdots\mathbf{W}_{l}^{(k+1)} \otimes\mathbf{W}_{l}^{(k-1)}\cdots\mathbf{W}_{l}^{(1)}\), \(k=1,2,3\); \(\mathbf{X}_{(k)}^{l}\) is the mode-\(k\) unfolding of \(\mathbf{\mathcal{X}}_{l}\) and \(\mathbf{X}_{(k)}^{l+1}\) is the mode-\(k\) unfolding of \(\mathbf{\mathcal{X}}_{l+1}\). Similarly, the gradient with respect to the input tensor can be obtained by vectorizing the tensor, given by

\[\frac{\partial\mathrm{vec}(\mathbf{\mathcal{X}}_{l+1})}{\partial\mathrm{vec}(\mathbf{ \mathcal{X}}_{l})}=\frac{\partial(\mathbf{W}_{l}^{(3)}\otimes\mathbf{W}_{l}^{(2)} \otimes\mathbf{W}_{l}^{(1)})\mathrm{vec}(\mathbf{\mathcal{X}}_{l})}{\partial \mathrm{vec}(\mathbf{\mathcal{X}}_{l})}. \tag{25}\]
Then the gradient of the loss function with respect to the parameters can be calculated by the chain rule.

## Appendix D Further experimental results

Table 6 shows the RMSEs of various TNN models at different sampling ratios. The model configurations used are the same as in Table 4.
