\section{Introduction}

{D}{eep} neural networks (DNNs) have achieved great success in solving complex tasks, such as image processing \cite{RN69, RN27}, language processing \cite{RN44, NIPS2017_3f5ee243}, object detection \cite{RN70, 7780460}, and medical diagnostics \cite{RN68, Rn81}. Empowered by the backpropagation (BP) algorithm, DNNs have become the mainstream approach in almost all cognitive applications. However, the standard BP suffers from the well-known backward locking problem that only permits the updates of a module after all dependent modules finish execution of both forward and backward passes \cite{10.5555/3305381.3305549}. This problem arises as errors are propagated backward from the top in a layer-by-layer fashion to update downstream module parameters. It restricts the network to perform training in a sequential manner. Intermediate tensors and operations necessary for module updates must be saved during the forward pass, causing high memory cost and frequent memory access \cite{RN18}. The memory constraints generally impose limitations in the training of the state-of-the-art DNNs on high-resolution inputs and large batch sizes. The strong inter-layer backward dependency prohibits training parallelization, slowing down the training process. This inefficient training process also holds back DNNs from being deployed in resource-constrained platforms such as edge devices.

The major obstacle imposed by backward locking is the heavy dependency between layers due to the error feedback signals. Various local training methods are introduced to tackle this difficulty by cutting off the feedback path \cite{10.5555/3305381.3305549, RN18, RN19, pmlr-v80-huang18b, RN41, pmlr-v119-belilovsky20a, wang2021revisiting, RN76, RN82}. These methods split the network into multiple modules, attach an auxiliary network to each module, and train each module separately and simultaneously or asynchronously. Since training happens locally, intermediate states can be saved in buffers temporally before parameter updates, eliminating the need for memory storage and access. Moreover, local training allows for the parallel execution of the forward and backward pass, significantly accelerating the training process.

Some works were reported to train networks in a layer-wise cascading fashion \cite{RN19, pmlr-v80-huang18b}. In other words, one layer is fully trained with an auxiliary classifier before moving on to train the next layer. For example, Marquez et al. proposed to train each layer using a multiple-layer fully-connected classifier in a cascading fashion to improve training efficiency \cite{RN19}. However, the performance of these approaches is usually limited to small datasets. Other works focus on designing locally supervised auxiliary networks to produce error feedback signals to train individual network modules synchronously \cite{10.5555/3305381.3305549,  pmlr-v80-huang18b, RN41, pmlr-v119-belilovsky20a, wang2021revisiting, RN76, RN82}. Jaderberg et al. proposed a decoupled neural interface method that synthesizes gradients locally to approximate the true gradients and eliminates the need to receive feedback signals \cite{10.5555/3305381.3305549}. Belilovsky et al. reported using a convolutional classifier to scale the network performance to ImageNet \cite{RN76}. Nokland et al. applied two cla  ssifiers equipped with different loss functions to train each layer \cite{RN41}. Wang et al. analyzed the defect of greedy local training methods from an information perspective and proposed a two-branch auxiliary network aiming to compute information propagation loss \cite{wang2021revisiting}. Significant improvement was demonstrated in various tasks. However, all these works are built on top of the same conventional greedy local training configuration that completely cuts off the backward path between local modules and prevents errors from flowing between them and hence the exchange of information.

This work aims to restore the inter-layer backward dependency and facilitate the information flow to improve network performance while preserving the computational advantages of local training methods. We introduce BackLink, a new local training algorithm that permits the errors generated in the current network module to flow back to its predecessor along a restricted propagation path. We perform extensive experiments to evaluate its effectiveness and analyze the computational costs in terms of memory consumption and computational time.
The main contributions of this work are summarized as follows.
\begin{itemize}
    \item Proposes a novel local training algorithm that introduces inter-module backward dependency to facilitate information flow and improve network performance with small computational overheads.
    \item Extensive experiments demonstrate that the proposed method outperforms the state-of-the-art methods \cite{pmlr-v119-belilovsky20a, wang2021revisiting} in classification performance. Increases in accuracy of up to 3.82\% and 4.05\% are observed on CIFAR10 for ResNet32 and ResNet110, respectively.
    \item Memory analysis and multi-GPU implementation reveal that the proposed method achieves largely reduced memory costs and significant training acceleration over BP. A 79\% reduction in memory cost and a 52\% decrease in runtime can be achieved in ResNet110.
\end{itemize}

\begin{figure*}[!t]
    \centering
    \includegraphics[width=\linewidth]{Figures/Fig_1.jpg}
    \caption{Training process for (a) the standard BP, (b) the greedy local training method and (c) the proposed local training method with the backward link. $l$ is the propagation length of the errors. $\delta_{l}^{n}$ and $\delta_{g}^{n}$ are the local and global error of the n-th hidden layer, respectively. $W_{c}^{n}$ and $W^{n+1}$ are the weights of the local classifier and the (n+1)-th hidden layer, respectively.}
    \label{fig:1}
\end{figure*}

The remainder of this work consists of the following sections. Section II describes the details of the proposed local training algorithm. Section III discusses the experimental results in various networks. Section IV analyzes the computational costs. Section V discusses performance trade-off and concludes this work.

\section{Method}
The backpropagation algorithm has been the standard training method in deep learning. The training process is depicted in Fig. \ref{fig:1} (a). Assume that a network consists of N layers. In the forward pass, the n-th hidden layer receives an input $x^{n}$ and computes its activations $y^{n}$  by
\begin{equation}
    y^{n}=f(z^{n})=f(W^{n}x^{n}+b^{n})
\end{equation}
where $W^{n}$ is the weight matrix of this layer, $b^{n}$ is the bias vector, $z^{n}$ is the synaptic input, and $f$ is the non-linear activation function, such as Rectified Linear Units (ReLUs). Instead of classifying outputs at the final layer, local training attaches a classifier to the end of each pre-defined local module of the network and predicts the local outputs. The errors generated locally are used to update the module parameters. We refer to the conventional local training method as the greedy local training (GLL) hereafter \cite{RN76}. In the GLL setting, as depicted in Fig. \ref{fig:1} (b), the error flow is prohibited from entering the preceding module, leading to independent training processes. In the backward pass, at the n-th hidden layer, the errors are computed by
\begin{equation}
    \delta^{n} = \frac{\partial L}{\partial z^{n}} = f^{'}(z^{n})\odot (W^{n+1})^{T} \delta^{n+1}
\end{equation}
where $L$ is the loss function, $\odot$ is the element-wise multiplication. The weight gradients are calculated by
\begin{equation}
    \Delta W^{n} = \delta^{n} \otimes x^{n}
\end{equation}
where $\otimes$ is the outer product.

In the proposed local training method, as illustrated in Fig. \ref{fig:1} (c), the errors are propagated backward from one module to its predecessor. Information carried by the error signals can thus be passed between modules. Instead of only focusing on improving the local classification performance, the current module updates its parameters accordingly to improve the performance of subsequent modules. However, to preserve the computational benefits of local training, the errors are restricted from traveling all the way down the module. We define the propagation length of the errors in a module as a hyperparameter, denoted as $l$. For example, $l=2$ is shown in Fig. \ref{fig:1} (c). In this particular case, the greedy local training method corresponds to $l=0$. Any layers located in the propagation range, e.g., the layers in green color in Fig. \ref{fig:1} (c), receive two error contributions we term as local and global contributions, respectively, when updating their errors. These layers can be unfolded into two identical parallel branches that share the same parameters and synchronize parameter updates resulting from different error contributions. If the n-th layer is the last layer of the module, the errors are derived as
\begin{align}
    \delta^{n} &= \delta_{l}^{n} + \delta_{g}^{n} \\
    &= f^{'}(z^{n})\odot [\alpha (W_{c}^{n})^{T} \delta_{c}^{n} + (1-\alpha)(W^{n+1})^{T} \delta^{n+1}]
\end{align}
where $\delta_{l}^{n}$ is the local error vector derived from the errors $\delta_{c}^{n}$ of the local classifier, $\delta_{g}^{n}$ is the global error vector derived from the errors $\delta^{n+1}$ of the next hidden layer, $W_{c}^{n}$ is the weights of the local classifier, and $\alpha$ is the weighting factor. For the layers located before, the errors are derived as
\begin{align}
    \delta^{i} &= \delta_{l}^{i} + \delta_{g}^{i} \\
    &= f^{'}(z^{i})\odot [(W^{i+1})^{T}(\alpha \delta_{l}^{i+1} + (1-\alpha)\delta_{g}^{i+1})]
\end{align}
The unfolding creates two separate network paths and decouples the two adjacent modules so that training can be performed independently.

\section{Classification experiments and results}
\subsection{Experiment setup}
The proposed training method was applied to train three types of widely-used CNNs, namely, AlexNet \cite{RN27}, VGG16 \cite{SimonyanZ14a}, and ResNet \cite{RN42}. The network structures are adapted from the original implementations. In all the networks, convolution operations are performed with a kernel size of 3×3, a stride of 1, and a padding of 1. In ResNets, the residual block consists of a stack of two convolution layers with an identity or down-sampling shortcut. The down-sampling happens at the first residual block where the number of channels doubles and is performed by a 1×1 convolution with a stride of 2. Two types of local classifiers are studied: a simple fully-connected (FC) layer and a convolutional network. The local convolutional network consists of one 3×3 convolutional layer and two FC layers with a hidden feature size of 128.  CIFAR10 and CIFAR100 datasets are used to evaluate the classification performance of the proposed training method \cite{RN84}.
For local training, the networks are split into multiple modules. The number of modules, denoted as K, is chosen as a multiple of 2, varying from 2 to 16 (8 in AlexNet). In ResNets, the residual block is treated as a basic unit (or layer) to preserve the residual structure. We divide the network evenly so that each local module consists of the same number of layers. If the number of basic layers is not divisible by K, we assign one more layer to lower modules. For example, in ResNet110, the number of basic layers is 55, and it can be split into $\{4\}\times7+\{3\}\times9$ for $K=16$ modules. The backward propagation length of the errors, $l$, is adjusted from 0 to 4. The weighting factor $\alpha$ is selected from $\{0,\ 0.25,\ 0.5,\ 0.75,\ 1\}$ and optimized in each network for each classification task.
All the experiments are conducted in the Pytorch framework. A cross-entropy loss function is used at each local classification layer. Stochastic gradient descent (SGD) is used as the optimizer. The momentum is set to 0.9, and the weight decay is set to 0.0001 for VGG16 and 0.0005 for ResNets, respectively. The batch size is set to 512 in all the experiments.  Training is run for 100 epochs in AlexNet, 150 epochs in VGG16, and 200 epochs in ResNets. Learning rates and their schedules are optimized in each network for each classification task. Dropout is only applied in fully-connected layers with the probability of 0.5 \cite{RN34}. ReLU is used as the non-linear function, before which batch normalization was applied.

\subsection{Classification results with local linear classifiers}
CIFAR10 is a collection of RGB frame images from different objects of 10 classes. It is divided into 50,000 training images and 10,000 testing images. We apply the four CNNs with different local training configurations to classify the dataset. Experiments on the impact of different propagation lengths and the number of modules are conducted. Learning rates are optimized in each network, which are 0.01 in AlexNet, 0.01 in VGG16, 0.5 in ResNet32, 0.3 in ResNet110, respectively.
The average classification errors in 5 trials are presented in Table \ref{tab:1}. The errors resulting from the standard BP are also reported. The number of local modules is changed from 2 to 16. The maximum propagation length is limited by the number of layers inside each local module. It can be observed that the classification performance of the conventional method degrades rapidly with the number of local modules. This is because in the conventional GLL method, without backward connection, local modules only focus on performing their classification task and ignore the demand from the next module, thus causing information loss while updates progress along with the network \cite{wang2021revisiting}. Our method provides a pathway for error signals to travel between modules, therefore retaining useful information to be used by the next modules. It outperforms the GLL method in all cases. The improvement increases with the propagation length, suggesting that more information is passed between modules. The difference becomes more and more prominent as the number of local modules doubles. For example, in ResNet32, with $K=16$ and $l=1$, our method surpasses the conventional method by 4.00\% in accuracy.

\begin{table}[t]
\centering
\caption{CLASSIFICATION ERRORS ON CIFAR10 AND CIFAR100 DATASETS WHEN LINEAR CLASSIFIERS ARE APPLIED. ERRORS ARE PERCENTAGES. K IS THE NUMBER OF LOCAL MODULES. GLL REFERS TO THE CONVENTIONAL GREEDY LOCAL TRAINING.}
\label{tab:1}
\begin{tabular}{|cccccc|}
\hline
\multicolumn{6}{|l|}{\textbf{CIFAR10}}                                                                                \\ \hline
\multicolumn{1}{|c|}{Network}                                                                                & \multicolumn{1}{c|}{Method}       & \multicolumn{1}{c|}{K=2}            & \multicolumn{1}{c|}{K=4}            & \multicolumn{1}{c|}{K=8}            & K=16           \\ \hline
\multicolumn{1}{|c|}{\multirow{3}{*}{\begin{tabular}[c]{@{}c@{}}AlexNet\\    \\ (BP: 11.44)\end{tabular}}}   & \multicolumn{1}{c|}{\textit{GLL}} & \multicolumn{1}{c|}{12.43}          & \multicolumn{1}{c|}{12.44}          & \multicolumn{1}{c|}{13.89}          &                \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=1}} & \multicolumn{1}{c|}{12.24}          & \multicolumn{1}{c|}{12.43}          & \multicolumn{1}{c|}{\textbf{12.45}} &                \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=2}} & \multicolumn{1}{c|}{\textbf{11.93}} & \multicolumn{1}{c|}{\textbf{11.96}} & \multicolumn{1}{c|}{}               &                \\ \hline
\multicolumn{1}{|c|}{\multirow{3}{*}{\begin{tabular}[c]{@{}c@{}}VGG16\\    \\ (BP: 7.41)\end{tabular}}}      & \multicolumn{1}{c|}{\textit{GLL}} & \multicolumn{1}{c|}{10.93}          & \multicolumn{1}{c|}{11.34}          & \multicolumn{1}{c|}{13.65}          & 18.02          \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=1}} & \multicolumn{1}{c|}{10.68}          & \multicolumn{1}{c|}{11.18}          & \multicolumn{1}{c|}{11.74}          & \textbf{12.72} \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=2}} & \multicolumn{1}{c|}{\textbf{10.54}} & \multicolumn{1}{c|}{\textbf{10.58}} & \multicolumn{1}{c|}{\textbf{10.82}} &                \\ \hline
\multicolumn{1}{|c|}{\multirow{5}{*}{\begin{tabular}[c]{@{}c@{}}ResNet32\\    \\ (BP: 7.51)\end{tabular}}}   & \multicolumn{1}{c|}{\textit{GLL}} & \multicolumn{1}{c|}{10.35}          & \multicolumn{1}{c|}{15.02}          & \multicolumn{1}{c|}{20.59}          & 24.21          \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=1}} & \multicolumn{1}{c|}{10.02}          & \multicolumn{1}{c|}{14.27}          & \multicolumn{1}{c|}{17.53}          & \textbf{20.21} \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=2}} & \multicolumn{1}{c|}{9.03}           & \multicolumn{1}{c|}{12.92}          & \multicolumn{1}{c|}{\textbf{16.23}} &                \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=3}} & \multicolumn{1}{c|}{9.22}           & \multicolumn{1}{c|}{10.65}          & \multicolumn{1}{c|}{}               &                \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=4}} & \multicolumn{1}{c|}{\textbf{9.18}}  & \multicolumn{1}{c|}{\textbf{10.63}} & \multicolumn{1}{c|}{}               &                \\ \hline
\multicolumn{1}{|c|}{\multirow{5}{*}{\begin{tabular}[c]{@{}c@{}}ResNet110\\    \\ (BP: 6.63)\end{tabular}}}  & \multicolumn{1}{c|}{\textit{GLL}} & \multicolumn{1}{c|}{8.55}           & \multicolumn{1}{c|}{13.19}          & \multicolumn{1}{c|}{14.97}          & 17.9           \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=1}} & \multicolumn{1}{c|}{8.34}           & \multicolumn{1}{c|}{12.6}           & \multicolumn{1}{c|}{14.92}          & 17.51          \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=2}} & \multicolumn{1}{c|}{7.99}           & \multicolumn{1}{c|}{12.56}          & \multicolumn{1}{c|}{14.64}          & 16.18          \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=3}} & \multicolumn{1}{c|}{\textbf{7.96}}  & \multicolumn{1}{c|}{\textbf{12.38}} & \multicolumn{1}{c|}{14.38}          & \textbf{15.62} \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=4}} & \multicolumn{1}{c|}{8.25}           & \multicolumn{1}{c|}{12.54}          & \multicolumn{1}{c|}{\textbf{13.62}} &                \\ \hline
\multicolumn{6}{|l|}{\textbf{CIFAR100}}                                                                                \\ \hline
\multicolumn{1}{|c|}{Network}                                                                                & \multicolumn{1}{c|}{Method}       & \multicolumn{1}{c|}{K=2}            & \multicolumn{1}{c|}{K=4}            & \multicolumn{1}{c|}{K=8}            & K=16           \\ \hline
\multicolumn{1}{|c|}{\multirow{5}{*}{\begin{tabular}[c]{@{}c@{}}ResNet32\\    \\ (BP: 29.68)\end{tabular}}}  & \multicolumn{1}{c|}{\textit{GLL}} & \multicolumn{1}{c|}{34.02}          & \multicolumn{1}{c|}{42.44}          & \multicolumn{1}{c|}{48.32}          & 51.75          \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=1}} & \multicolumn{1}{c|}{33.60}          & \multicolumn{1}{c|}{41.92}          & \multicolumn{1}{c|}{45.12}          & \textbf{46.76} \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=2}} & \multicolumn{1}{c|}{32.56}          & \multicolumn{1}{c|}{38.95}          & \multicolumn{1}{c|}{\textbf{42.59}} &                \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=3}} & \multicolumn{1}{c|}{32.64}          & \multicolumn{1}{c|}{36.61}          & \multicolumn{1}{c|}{}               &                \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=4}} & \multicolumn{1}{c|}{\textbf{32.47}} & \multicolumn{1}{c|}{\textbf{35.99}} & \multicolumn{1}{c|}{}               &                \\ \hline
\multicolumn{1}{|c|}{\multirow{5}{*}{\begin{tabular}[c]{@{}c@{}}ResNet110\\    \\ (BP: 28.33)\end{tabular}}} & \multicolumn{1}{c|}{\textit{GLL}} & \multicolumn{1}{c|}{31.37}          & \multicolumn{1}{c|}{37.22}          & \multicolumn{1}{c|}{41.73}          & 45.23          \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=1}} & \multicolumn{1}{c|}{30.85}          & \multicolumn{1}{c|}{37.5}           & \multicolumn{1}{c|}{41.49}          & 44.82          \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=2}} & \multicolumn{1}{c|}{30.95}          & \multicolumn{1}{c|}{36.94}          & \multicolumn{1}{c|}{41.34}          & 43.85          \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=3}} & \multicolumn{1}{c|}{30.54}          & \multicolumn{1}{c|}{36.37}          & \multicolumn{1}{c|}{41.05}          & \textbf{42.09} \\ \cline{2-6}
\multicolumn{1}{|c|}{}                                                                                & \multicolumn{1}{c|}{\textit{l=4}} & \multicolumn{1}{c|}{\textbf{30.28}} & \multicolumn{1}{c|}{\textbf{35.63}} & \multicolumn{1}{c|}{\textbf{38.88}} &                \\ \hline
\end{tabular}
\end{table}

CIFAR100 has the same dataset size as CIFAR10, but it has 100 classes. It is a more challenging classification task, a good benchmark for evaluating the proposed method. We apply ResNet32 and ResNet110 with different local training configurations to classify the dataset. The same experiments are performed to analyze the impact of our method. The average classification errors in 5 trials are presented in Table \ref{tab:1}. The same observations can be made as follows. Increasing the number of local modules degrades the classification performance. The degradation becomes more severe than in the case of CIFAR10. Our method reduces classification errors in all cases. Especially, the reduction becomes significant when the number of local modules is large. For example, in ResNet32, with $K=16$ and $l=1$, our method surpasses the GLL method by 4.99\% in accuracy. Therefore, our method is demonstrated to be consistent in improving classification performance of local training algorithms. Additionally, ResNet32 benefits more from the backward dependency than ResNet110, as the improvement is more evident. Particularly, when $K=4$, ResNet32 achieves better accuracy on CIFAR10 and similar accuracy on CIFAR100.

\subsection{Classification results with local convolutional classifiers}
Increasing the complexity of the local classifiers helps train the local modules to produce features containing useful information for the subsequent modules \cite{pmlr-v119-belilovsky20a, RN76}. We perform experiments in ResNet32 and ResNet110 with local convolutional classifiers consisting of one convolutional layer and two FC layers. The classification results on CIFAR10 are obtained in Table \ref{tab:2}. We include two state-of-the-art local training methods for comparisons. The decoupled greedy learning (DGL) was proposed to parallelize module computations in both forward and backward pass by means of local replay buffers and classifiers \cite{pmlr-v119-belilovsky20a}. Wang et al. analyzed the drawbacks of local training methods from information perspectives and proposed an information propagation (InfoPro) loss in the local optimization objective to preserve information in local modules \cite{wang2021revisiting}. The same local classifiers are adopted for fair comparisons. It is clear that in all the methods, the accuracy drops as the number of local modules increases. Overall, our method achieves the best classification performance. Compared with the InfoPro method, our method shows notably better accuracy in ResNet32 and slight increase in ResNet110. The InfoPro method relies on a complex local network architecture composed of a convolutional classifier and a convolutional decoder and introduces four additional hyperparameters for tuning. The decoder is built on a bilinear interpolation layer and two convolutional layers to reconstruct the input image locally. In contrast, our method has much lower complexity.
\begin{table}[t]
\centering
\caption{CLASSIFICATION ERRORS ON CIFAR10 DATASET WHEN CONVOLUTIONAL CLASSIFIERS ARE APPLIED. ERRORS ARE PERCENTAGES. K IS THE NUMBER OF LOCAL MODULES. GLL REFERS TO THE CONVENTIONAL GREEDY LOCAL TRAINING.}
\label{tab:2}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
Network                                                                               & Method                   & K=2           & K=4           & K=8            & K=16           \\ \hline
\multirow{7}{*}{\begin{tabular}[c]{@{}c@{}}ResNet32\\    \\ (BP: 6.94)\end{tabular}}  & \textit{GLL}             & 7.96          & 10.61         & 13.61          & 14.97          \\ \cline{2-6}
                                                                                & \textit{DGL{[}14{]}}     & 8.69          & 11.48         & 14.17          & 16.22          \\ \cline{2-6}
                                                                                & \textit{InfoPro{[}15{]}} & 8.13          & 8.64          & 11.4           & 14.23          \\ \cline{2-6}
                                                                                & \textit{l=1}             & 7.44          & 9.65          & 11.7           & \textbf{12.40} \\ \cline{2-6}
                                                                                & \textit{l=2}             & 7.17          & 8.66          & \textbf{10.60} &                \\ \cline{2-6}
                                                                                & \textit{l=3}             & 7.13          & 7.95          &                &                \\ \cline{2-6}
                                                                                & \textit{l=4}             & \textbf{6.97} & \textbf{7.55} &                &                \\ \hline
\multirow{7}{*}{\begin{tabular}[c]{@{}c@{}}ResNet110\\    \\ (BP: 6.27)\end{tabular}} & \textit{GLL}             & 7.45          & 10.07         & 11.87          & 12.99          \\ \cline{2-6}
                                                                                & \textit{DGL{[}14{]}}     & 7.70          & 10.50         & 12.46          & 13.80          \\ \cline{2-6}
                                                                                & \textit{InfoPro{[}15{]}} & 7.01          & 7.96          & 9.40           & 10.78          \\ \cline{2-6}
                                                                                & \textit{l=1}             & 6.47          & 8.65          & 9.95           & 10.77          \\ \cline{2-6}
                                                                                & \textit{l=2}             & 6.39          & 7.96          & 9.56           & 10.61          \\ \cline{2-6}
                                                                                & \textit{l=3}             & \textbf{6.36} & 7.92          & 9.29           & \textbf{9.75}  \\ \cline{2-6}
                                                                                & \textit{l=4}             & 6.38          & \textbf{7.79} & \textbf{9.25}  &                \\ \hline
\end{tabular}
\end{table}

\section{Computational costs}
Local training algorithms are the favored alternative approach to the standard BP because of their low memory cost and excellent model parallelization. In this section, we will analyze the computational cost of the proposed method and compare it with other methods.
\subsection{GPU memory cost}
The greedy local training splits networks into gradient-isolated modules, eliminating the necessity to save the whole computational graph and intermediate tensors and leading to reduced memory cost. While our method allows errors to flow between modules and hence imposes the backward dependency, we limit the propagation length within the local module, reduce the extent of dependency, and preserve the advantage of low memory usage.

\begin{figure}[b]
    \centering
    \includegraphics[width=\linewidth]{Figures_1/Fig_2.jpg}
    \caption{GPU memory cost of training (a) ResNet32 and (b) ResNet110 with local convolutional classifiers on a single GPU. The standard BP corresponds to 1 local module.}
    \label{fig:2}
\end{figure}

We measure the maximum GPU memory occupied by ResNets with local convolutional classifiers under different settings on a single Nvidia Titan RTX GPU. The measurement is done with the commonly-used command \textit{max\_memory\_allocated} in Pytorch \cite{wang2021revisiting, RN57}. The results of training ResNet32 and ResNet110 are shown in Fig.\ref{fig:2} (a) and (b), respectively. As expected, the memory cost decreases significantly with the number of local modules. Compared with the standard BP, local training can lead up to 69\% in ResNet32 and 81\% in ResNet110. Applying our method introduces a small overhead, which is 23\% at maximum in ResNet32 and 12\% in ResNet110. In most cases, the overhead is smaller than 5\%. Moreover, due to the complexity of the local network, the InfoPro method consumes much more GPU memory in ResNet32 with an overhead of 71\% when the number of modules is 16. However, the overhead becomes less in ResNet110 because the size of the local network is much smaller than the local module. Therefore, compared with the GLL method, our method causes a smaller overhead in GPU memory cost than the InfoPro method.
\subsection{Runtime cost on multiple GPUs}
\begin{figure}
    \centering
    \includegraphics[width=0.6\linewidth]{Figures/Fig_3.jpg}
    \caption{Multi-GPU implementation of the proposed local training method. $g_{l}$ and $g_{g}$ are the local gradients and global gradients, respectively. The green circles indicate the hidden layers with the backward link.}
    \label{fig:3}
\end{figure}

Local training algorithms have inherent model parallelization that enables the parallel execution of the forward pass of a module and the backward pass of its previous module. To leverage the parallelization, we can implement networks on multiple GPUs, each running a local module as a single process. Each GPU receives the inputs from the previous one except the very first one that gets the inputs from the dataset. The inputs are saved in shared memory. In this way, each GPU performs training in the local module independently, leading to reduced training runtime.
Similarly, our training method can be implemented on multiple GPUs to achieve training acceleration. The implementation is illustrated in Fig. \ref{fig:3}. The example shows an 8-layer network split into two local modules with an error propagation length of 2. Each GPU trains one module. We duplicate the last two layers (in green) of the first module into the second module to compute their gradients resulting from the global errors on the same GPU. The local gradients are computed on GPU0. Then the final gradients of these two layers are synchronized across GPUs by exchanging the gradients through shared memory. The gradients from the two layers on the same GPU can be concatenated and sent out to the shared memory at once to reduce the communication time. In this way, training on different GPUs is performed independently with an overhead of gradient communication time.
We utilize the multiprocessing package in Pytorch to create multiple subprocesses and necessary shared memories \cite{pmlr-v119-belilovsky20a, osti_10074580}. We split ResNet110 into K modules and measure the simulation runtime on K GPUs. The convolutional classifier is used in local modules. Runtime is measured on two different types of GPU models, namely Nvidia Tesla P100 and V100, as presented in Fig. \ref{fig:4} (a) and (b), respectively. The runtime is normalized over the standard BP. From Fig. \ref{fig:4}, we can tell that GPU models affect the local training speed. The GLL method significantly reduce the runtime by up to 53\% on P100 GPUs and 23\% on V100 GPUs compared with the standard BP, respectively. For the P100 model, the runtime is improved as the number of GPUs increases, while little or even no runtime improvement can be observed with more than two GPUs of the V100 models. The difference could be attributed to the GPU processing speed and the communication time between GPUs and the shared memory. More importantly, our method runs faster than the InfoPro method and causes a small runtime overhead, which are less than 5\% for the P100 model and 10\% for the V100 model in most of the cases, respectively.
\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{Figures_1/Fig_4.jpg}
    \caption{Runtime comparison among different local training methods in ResNet110 on multiple GPUs for two different GPU models: (a) Nvidia Tesla P100 and (b) Nvidia Tesla V100. The runtime is normalized over the standard BP.}
    \label{fig:4}
\end{figure}

\section{Discussion and conclusion}
The proposed BackLink local training algorithm is demonstrated to produce accuracy close to the standard BP, outperforming the state-of-the-art local training methods. It has shown a great advantage in lowering GPU memory cost and runtime by increasing the number of local modules. However, classification performance is also impacted. Thus, we present the performance trade-off among the accuracy drop, GPU memory and runtime, as shown in Fig. \ref{fig:5}. For the same propagation length $l$, increasing the number of local modules incurs more accuracy drop, but leads to higher memory and runtime reduction. Specifically, in the case of $l=2$, when the number of local modules doubles, 1.57\% drop in accuracy is caused, while around 23\% more reduction in memory and runtime are achieved. On the other hand, increasing the propagation length from 1 to 4 leads to improvement in accuracy, i.e. 0.87\%, with negligible overhead, since the memory reduction remains unchanged and the runtime is slightly affected (< 6\%). Accuracy improevement is more significant with a larger K, as shown in Table \ref{tab:1} and Table \ref{tab:1}.
\begin{figure}[b]
    \centering
    \includegraphics[width=0.6\linewidth]{Figures/Fig_5.jpg}
    \caption{Performance trade-off resulting from different number of local modules for the error propagation length changing from 1, 2, 3, to 4. Accuracy drop, GPU memory and runtime reduction are computed with the standard BP as the baseline. Results are obtained from ResNet110 on CIFAR10.}
    \label{fig:5}
\end{figure}
In this work, we have demonstrated a new local training algorithm that introduces a backward dependency between local modules by allowing errors to flow between them and guides information to travel along with the network from top to bottom. The proposed algorithm tackles the shortcoming of the conventional local training algorithm that loses information while features are propagated forward along with the network. The backward dependency ensures that the preceding modules preserve enough information in the output features to be utilized by the subsequent modules. The network can achieve better global classification performance instead of only focusing on improving individual local performance. Extensive experiments on different types of CNNs and benchmark datasets have demonstrated the effectiveness of the proposed algorithm and considerable improvement in classification performance. For example, in ResNet32 with 16 local modules, our method surpasses the conventional greedy local training method by 4.00\% and the recent InfoPro method by 1.83\% in accuracy on CIFAR10, respectively. In order to preserve the computational advantages of local training algorithms, we lighten the backward dependency between modules by restricting the propagation length of the errors within a module. As a result, computational advantages of local training in GPU memory and multi-GPU simulation runtime, are marginally impacted, while significant reductions are still achieved in comparison to the standard BP. Our method can lead up to a 79\% reduction in memory cost and 53\% in simulation runtime in ResNet110 compared to the standard BP. Therefore, we believe that our algorithm brings new possibilities for developing more efficient and biologically plausible deep learning algorithms.


\section{Introduction}\label{sec:introduction}

Secure distributed matrix multiplication (SDMM) has been studied as a way to compute a matrix product using the help of worker servers such that the computation is information-theoretically secure against colluding workers. SDMM was first studied by Chang and Tandon in \cite{chang2018capacity}. Their scheme was improved by D'Oliveira \emph{et al.}\ in \cite{d2020gasp, d2021degree, d2020notes} using GASP codes. Different codes have also been introduced in \cite{kakar2019capacity, aliasgari2020private, jia2019capacity, mital2020secure, kim2019private, yu2020entangled}. Furthermore, different modes of SDMM, such as private, batch or cooperative SDMM, have been studied in \cite{chang2019upload, jia2021cross, yu2019lagrange, yu2020entangled, chen2021gcsa, zhu2021secure, li2022efficient, li2022private}. The information-theoretic capacity of SDMM has been studied in \cite{chang2018capacity, jia2019capacity, kakar2019capacity, yang2019secure}, but overall capacity results are still scarce. In addition to considering SDMM over finite fields, SDMM has also been utilized over the analog domain (real/complex numbers) in \cite{makkonen2022analog}.

The workers in an SDMM scheme are thought of as untrustworthy-but-useful, which means that some of them might not work as intended. The main robustness has been against providing security against colluding workers, which share the information they receive and try to infer the contents of the original matrices. Tools from secret sharing have been used to guarantee information-theoretic security against such colluding workers. Additionally, robustness against so-called straggling workers has been considered. Stragglers are workers that respond slowly or not at all. Such workers cause an undesired straggler effect if the computation is limited by the slowest worker.

Byzantine workers are workers that return erroneous results either intentionally or as a result of a fault. Such errors can be difficult to detect directly without further analysis. To guarantee correctness of the matrix product, it is crucial to be able to detect the errors and correct them with minimal overhead in communication and computation. Tools from classical coding theory can be used to correct errors caused by the Byzantine workers and erasures caused by stragglers.

A coded computation scheme that accounts for stragglers and Byzantine workers has been presented in \cite{yu2019lagrange} using so-called Lagrange coded computation. This scheme considers stragglers as erasures and Byzantine workers as errors in some linear codes. This means that a straggling worker requires one additional worker and a Byzantine worker requires two additional workers.

As the main contribution, this paper introduces a general framework for linear SDMM schemes that can be used to construct many SDMM schemes in a unified way. We show a strong connection between star product codes and SDMM schemes and relate the associated codes to the proofs of the security of the schemes as well as to the recovery threshold and collusion tolerance. Previously, star product codes have been successfully utilized in private information retrieval (PIR) \cite{freij2017private}. Furthermore, we present a bounded-distance decoding strategy utilizing interleaved codes, which provides robustness against stragglers and Byzantine workers. Finally, we analyze the error-correcting capabilities of the proposed strategy under some natural assumptions about the error distributions.

\section{Preliminaries on Interleaved Codes}\label{sec:preliminaries_interleaved_codes}

We write $[n] = \{1, \dots, n\}$. We consider matrices and scalars over a finite field $\ensuremath{\mathbb{F}}_q$ with $q$ elements. Vectors in $\ensuremath{\mathbb{F}}_q^n$ are considered to be row vectors. If $G$ is a matrix, then $G^{\leq m}$ and $G^{> m}$ denote the first $m$ rows and the the rest of the rows. Furthermore, if $\mathcal{I}$ is a set of indices, then $G_\mathcal{I}$ is the submatrix of $G$ with the columns indexed by $\mathcal{I}$.

Interleaved codes have been used to correct burst errors in a stream of codewords in many applications. Burst errors are errors where multiple consecutive symbols are affected instead of single symbol errors distributed randomly.

\begin{definition}[Homogeneous interleaved codes]
Let $\mathcal{C}$ be a linear code over the field $\ensuremath{\mathbb{F}}_q$. Then the $\ell$-interleaved code of $\mathcal{C}$ is the code
\begin{align*}
    \mathcal{IC}^{(\ell)} = \left\{ \begin{pmatrix}c_1 \\ \vdots \\ c_\ell \end{pmatrix} \colon c_i \in \mathcal{C}~\forall i \in [\ell] \right\}.
\end{align*}
\end{definition}

The codewords in an interleaved code are matrices, where each row is a codeword in the code $\mathcal{C}$. Instead of the Hamming weight as the measure of the size of an error, the column weight is used. The column weight of an error matrix is defined to be the number of nonzero columns.

When many codewords need to be transmitted, they can be sent such that first the first symbol of each codeword is sent, then the second symbol of each codeword, and so on. If a burst error occurs, then multiple codewords are affected, but only a small number of symbols is affected in a codeword. This transforms the burst error into single symbol errors in the individual codewords, which means that regular error correction algorithms can be used to correct up to half the minimum distance of errors.

Even more efficient error correction algorithms can be performed for interleaved codes by considering collaborative decoding, where all of the codewords in the interleaved code are considered at the same time. This is advantageous, since the error locations in each of the codewords are the same. Collaborative decoding algorithms have been studied in \cite{krachkovsky2003reed, schmidt2009collaborative} and more recently in \cite{holzbaur2021success}. Collaborative decoding algorithms can achieve beyond half the minimum distance decoding by correcting the errors as a system of simultaneous equations.

\section{Examples of SDMM Schemes}\label{sec:examples_of_SDMM_schemes}

In this section we recall some examples of previous SDMM schemes adopting the presentation typically used in the literature. Later, we will show how these schemes arise as special cases from the general framework proposed in this paper.

The goal is to compute the matrix product of the matrices $A \in \ensuremath{\mathbb{F}}_q^{t \times s}$ and $B \in \ensuremath{\mathbb{F}}_q^{s \times r}$ using a total of $N$ workers, while protecting against any $X$ colluding workers. Furthermore, we denote by $S$ the number of stragglers and by $B$ the number of Byzantine workers. These should not be confused with any matrices. The \emph{recovery threshold} defined as the number of responses from workers that are required in order to decode the intended product.

The schemes are based on different matrix partitioning techniques. The most popular partitionings are the \emph{inner product partitioning} (IPP) and \emph{outer product partitioning} (OPP). In IPP the matrices are partitioned into $p$ pieces such that
\begin{align*}
    A = \begin{pmatrix} A_1 & \dots & A_p \end{pmatrix}, \quad B = \begin{pmatrix} B_1 \\ \vdots \\ B_p \end{pmatrix}.
\end{align*}
Then the product can be expressed as $AB = \sum_{j=1}^p A_j B_j$. In OPP the matrices are partitioned into $m$ and $n$ pieces, respectively, such that
\begin{align*}
    A = \begin{pmatrix} A_1 \\ \vdots \\ A_m \end{pmatrix}, \quad B = \begin{pmatrix} B_1 & \dots & B_n \end{pmatrix}.
\end{align*}
Then the product can be expressed as
\begin{align*}
    AB = \begin{pmatrix}
    A_1B_1 & \dots & A_1B_n \\
    \vdots & \ddots & \vdots \\
    A_mB_1 & \dots & A_mB_n
    \end{pmatrix}.
\end{align*}
It is also possible to consider combinations of these partitions.

\textit{Secure MatDot:} A well-known example of an SDMM scheme is the secure MatDot scheme introduced in \cite{aliasgari2020private}. The secure MatDot scheme uses the inner product partitioning to partition the matrices to $p$ pieces. Define the polynomials
\begin{align*}
    f(x) &= \sum_{j=1}^p A_j x^{j-1} + \sum_{k=1}^X R_k x^{p+k-1} \\
    g(x) &= \sum_{j'=1}^p B_{j'} x^{p-j'-1} + \sum_{k'=1}^X S_{k'}x^{p+k'-1},
\end{align*}
where $R_1, \dots, R_X$ and $S_1, \dots, S_X$ are matrices of appropriate size that are chosen uniformly at random from matrices over $\ensuremath{\mathbb{F}}_q$. Let $\alpha_1, \dots, \alpha_N \in \ensuremath{\mathbb{F}}_q^\times$ be distinct nonzero points. Evaluate the polynomials $f(x)$ and $g(x)$ at these points to get the encoded matrices
\begin{align*}
    \widetilde{A}_i = f(\alpha_i), \quad \widetilde{B}_i = g(\alpha_i).
\end{align*}
These encoded matrices can be sent to each worker node. The worker nodes compute the matrix products $\widetilde{C}_i = \widetilde{A}_i \widetilde{B}_i$ and return these  to the user. The user receives evaluations of the polynomial $h(x) = f(x)g(x)$ from each worker node. Using the definition of $f(x)$ and $g(x)$ we can write out the coefficients of $h(x)$ as
\begin{align*}
    h(x) &= f(x)g(x) \\
    &= \sum_{j=1}^p \sum_{j'=1}^p A_jB_{j'} x^{p+j-j'-1} + (\text{terms of degree $\geq p$}).
\end{align*}
The degree of $h(x)$ is at most $2p + 2X - 2$. Furthermore, the coefficient of the term $x^{p-1}$ is exactly the product $AB$, which we wish to recover. Using polynomial interpolation we can compute the required coefficient, given that we have at least $2p + 2X - 1$ evaluations. The recovery threshold for the secure MatDot code is $K = 2p + 2X - 1$.

\textit{GASP code:} Another SDMM scheme that has been widely researched is the GASP scheme, which was first presented in \cite{d2020gasp}. This scheme is also based on polynomial evaluation, but the choice of the polynomials and the evaluation points is more involved. Additionally, the matrices are partitioned according to the outer product partitioning. The following example will give an idea of the general construction described in \cite{d2020gasp, d2021degree}.

The matrices $A \in \ensuremath{\mathbb{F}}_q^{t \times s}$ and $B \in \ensuremath{\mathbb{F}}_q^{s \times r}$ are partitioned into $m = n = 3$ submatrices with the outer product partitioning. We wish to protect against $X = 2$ colluding workers. Let us define the polynomials
\begin{align*}
    f(x) &= A_1 + A_2x + A_3x^2 + R_1x^9 + R_2x^{12}, \\
    g(x) &= B_1 + B_2x^3 + B_3x^6 + S_1x^9 + S_2x^{10},
\end{align*}
where $R_1, R_2, S_1, S_2$ are matrices of appropriate size that are chosen uniformly at random from matrices over $\ensuremath{\mathbb{F}}_q$. The exponents are chosen carefully so that the total number of workers needed is as low as possible. Let $\alpha_1, \dots, \alpha_N \in \ensuremath{\mathbb{F}}_q^\times$ be distinct nonzero points. Evaluate the polynomials $f(x)$ and $g(x)$ at these points to get the encoded matrices
\begin{align*}
    \widetilde{A}_i = f(\alpha_i), \quad \widetilde{B}_i = g(\alpha_i).
\end{align*}
These encoded matrices can be sent to each worker node. The worker nodes compute the matrix products $\widetilde{C}_i = \widetilde{A}_i \widetilde{B}_i$ and send these  to the user. The user receives evaluations of the polynomial $h(x) = f(x)g(x)$ from each worker. Using the definition of $f(x)$ and $g(x)$ we can write out the coefficients of $h(x)$ as
\begin{align*}
    h(x) &= A_1B_1 + A_2B_1x + A_3B_1x^2 + A_1B_2x^3 + A_2B_2x^4 \\
    &+ A_2B_3x^5 + A_1B_3x^6 + A_2B_3x^7 + A_3B_3x^8 \\
    &+ (\text{terms of degree $\geq 9$}).
\end{align*}
We notice that the coefficients of the first $9$ terms are exactly the submatrices we wish to recover. We need 18 responses from the workers, since $h(x)$ has $18$ nonzero coefficients, provided that the corresponding linear equations are solvable. In this case, the recovery threshold is $K = 18$.

The general choice of the exponents in the polynomials $f(x)$ and $g(x)$ is explained in \cite{d2021degree}. A so-called degree table is used to analyze the recovery threshold of the scheme. Furthermore, the choice of the evaluation points is not as simple as with the secure MatDot code, but it has been shown that a suitable choice can be made in a large enough field.

\textit{SDMM based on FFT:} A novel SDMM scheme based on the Fourier transform over finite fields was introduced by Mital \emph{et al}.\ in \cite{mital2020secure}. This scheme is not able to protect against straggling workers. The matrices are partitioned into $p = N - 2X$ pieces with the inner product partitioning. Define the functions
\begin{align*}
    f(x) &= \sum_{j=1}^p A_j x^{j-1} + \sum_{k=1}^X R_k x^{p + k - 1} \\
    g(x) &= \sum_{j'=1}^p B_{j'} x^{-j'+1} + \sum_{k'=1}^X S_{k'} x^{-p-X-k'+1},
\end{align*}
where $R_1, \dots, R_X$ and $S_1, \dots, S_X$ are matrices of appropriate size that are chosen uniformly at random from matrices over $\ensuremath{\mathbb{F}}_q$. Let $\alpha_N \in \ensuremath{\mathbb{F}}_q^\times$ be a primitive $N$th root of unity. The functions $f(x)$ and $g(x)$ are evaluated at the points $1, \alpha_N, \alpha_N^2, \dots, \alpha_N^{N-1}$ and the results are sent to the workers such that worker $i \in [N]$ receives the encoded matrices
\begin{align*}
    \widetilde{A}_i = f(\alpha_N^{i-1}), \quad \widetilde{B}_i = f(\alpha_N^{i-1}).
\end{align*}
The workers compute the matrix products of the encoded matrices and return the results $\widetilde{C}_i = \widetilde{A}_i \widetilde{B}_i$. The user receives evaluations of the function
\begin{align*}
    h(x) = f(x)g(x) = \sum_{j=1}^p A_j B_j + (\text{non-constant terms}).
\end{align*}
The other terms have degree in $[-N + 1, N - 1]$, which means that the average of the responses equals the constant term, since $\sum_{i=1}^N (\alpha_N)^s = 0$ for $N \nmid s$. Hence, the product $AB$ can be computed as the average of all the responses. This means that no stragglers can be tolerated, since all of the responses are needed. Furthermore, the field has to be such that the appropriate $N$th root of unity exists.

\section{Linear SDMM}\label{sec:linear_SDMM}

Many SDMM schemes in the literature use concepts from coding theory and secret sharing, but are usually presented as concrete constructions based on polynomial interpolation. This makes it easy to argue that the schemes compute the desired matrix product, but the comparison of different schemes is difficult. A more general and abstract way can provide simpler comparisons between SDMM schemes and constructions that are not based on any particular SDMM scheme, while losing some detail about why each scheme works as they do. In this section we present a general linear SDMM framework that can be used to describe the earlier SDMM schemes in a compact way. Furthermore, we prove a general security result for linear SDMM schemes.

\subsection{General Linear SDMM Framework via Star Products}\label{sec:general_linear_SDMM_framework}

Linear SDMM schemes can be constructed in general with the following formula. The input matrices $A \in \ensuremath{\mathbb{F}}_q^{t \times s}$ and $B \in \ensuremath{\mathbb{F}}_q^{s \times r}$ are partitioned into block matrices of equal size  using some partitioning such that the blocks can be arranged as the vectors
\begin{equation*}
    (A_1, \dots, A_m),\quad (B_1, \dots, B_n)
\end{equation*}
of length $m$ and $n$, respectively. Matrices $R_1, \dots, R_X$ and $S_1, \dots, S_X$ are drawn uniformly at random such that the matrices $R_k$ and $S_{k'}$ have the same size as $A_j$ and $B_{j'}$, respectively. These random matrices are appended to the partition vectors to get the following vectors of matrices
\begin{equation*}
    (A_1, \dots, A_m, R_1, \dots, R_X), \quad (B_1, \dots, B_n, S_1, \dots, S_X)
\end{equation*}
of length $m + X$ and $n + X$, respectively. These vectors are encoded using $[N, m + X]$ and $[N, n + X]$ linear codes $\mathcal{C}_A$ and $\mathcal{C}_B$, respectively, where $N$ denotes the number of workers. Let $F$ and $G$ be generator matrices for $\mathcal{C}_A$ and $\mathcal{C}_B$, respectively. The associated codewords are then
\begin{align*}
    \widetilde{A} &= (\widetilde{A}_1, \dots, \widetilde{A}_N) = (A_1, \dots, A_m, R_1, \dots, R_X) F \\
    \widetilde{B} &= (\widetilde{B}_1, \dots, \widetilde{B}_N) = (B_1, \dots, B_n, S_1, \dots, S_X) G.
\end{align*}
Each worker is sent one component of both vectors, \emph{i.e.}, worker $i$ receives matrices $\widetilde{A}_i$ and $\widetilde{B}_i$. The worker then multiplies these matrices together and sends the result to the user. In coding-theoretic terms, this can be interpreted as a star product of the vectors $\widetilde{A}$ and $\widetilde{B}$. Hence, we may write
\begin{align*}
    \widetilde{C} = \widetilde{A} \star \widetilde{B} = (\widetilde{A}_1 \widetilde{B}_1, \dots, \widetilde{A}_N \widetilde{B}_N).
\end{align*}
Let $H$ be a suitable generator matrix for $\mathcal{C}_A \star \mathcal{C}_B$ and $\mathcal{I} \subseteq [N]$ the indices of the fastest responses. The result can be obtained by computing $(\widetilde{A} \star \widetilde{B})_\mathcal{I} H_\mathcal{I}^{-1}$, provided that $\mathcal{I}$ is an information set of $\mathcal{C}_A \star \mathcal{C}_B$. The generator matrix $H$ has to be chosen such that the desired product can be obtained as one of the components of the decoded vector. It is often convenient to assume that $\mathcal{C}_A \star \mathcal{C}_B$ is an MDS code, which means that any set of sufficient size is an information set.

\begin{remark}
The vectors $\widetilde{A}$, $\widetilde{B}$ and $\widetilde{C}$ are vectors of matrices that are obtained as the encoding of the partition vectors. Another way to interpret these is that $\widetilde{A}$, $\widetilde{B}$ and $\widetilde{C}$ are matrices of codewords in $\mathcal{C}_A$, $\mathcal{C}_B$ and $\mathcal{C}_A \star \mathcal{C}_B$, respectively. By abuse of notation, we denote
\begin{align*}
    \widetilde{A} \in \mathcal{C}_A,~\widetilde{B} \in \mathcal{C}_B,~\widetilde{C} \in \mathcal{C}_A \star \widetilde{C}_B.
\end{align*}
\end{remark}

\begin{remark}
This construction of linear SDMM schemes is quite abstract as it does not provide a general way of constructing new SDMM schemes from any linear codes. However, it provides a robust and general way to look at different SDMM schemes and prove general results. This general linear structure can be used to easily show the security of the schemes by checking properties of the generator matrices $F$ and $G$ as we show in Section~\ref{sec:security_linear_SDMM}.
\end{remark}

With the linear SDMM framework, an SDMM scheme can be described using the partition of the matrices and the three generator matrices $F$, $G$ and $H$. It is worth noting that $H$ is a generator matrix of $\mathcal{C}_A \star \mathcal{C}_B$ and not a parity-check matrix.

If $\mathcal{C}_A \star \mathcal{C}_B$ is an MDS code, only a certain number of fastest responses are required for decoding. This is possible since any set that is large enough is an information set. In other words, the dimension $K$ of $\mathcal{C}_A \star \mathcal{C}_B$ is the recovery threshold of the SDMM scheme, which means that any $K$ responses are enough to decode the answer. Recovery threshold has been denoted by $R_c$ in the literature, but we choose to denote it by $K$ as it conveniently matches our coding-theoretic approach.

As an immediate consequence of \cite[Theorem 7]{mirandola2015critical} we get the following lower bound for the recovery threshold for a linear SDMM scheme.

\begin{proposition}\label{prop:recovery_thresold}
If $\mathcal{C}_A$ and $\mathcal{C}_B$ are full-support MDS codes, then the linear SDMM scheme using the codes $\mathcal{C}_A$ and $\mathcal{C}_B$ has recovery threshold
\begin{align*}
    K \geq \min\{N, (m + n) + 2X - 1\}.
\end{align*}
\end{proposition}

As a corollary of Proposition~\ref{prop:recovery_thresold} we get an upper bound on the number of colluding workers.

\begin{corollary}
The number of colluding workers a linear SDMM scheme can protect against is $X < \frac{N}{2}$, when $\mathcal{C}_A$ and $\mathcal{C}_B$ are full-support MDS codes.
\end{corollary}

\subsection{Constructing SDMM Schemes Using Linear SDMM}\label{sec:constructing_linear_SDMM_schemes}

The examples of SDMM schemes presented in Section~\ref{sec:examples_of_SDMM_schemes} can be described using the linear SDMM framework by describing the partitioning of the matrices and the associated codes $\mathcal{C}_A$ and $\mathcal{C}_B$. Here we describe only the generator matrices of $\mathcal{C}_A$ and $\mathcal{C}_B$. The generator matrix of $\mathcal{C}_A \star \mathcal{C}_B$ is also of importance, but in the interest of space we will not describe it in detail.

\textit{Secure MatDot:} The secure MatDot scheme can be described using the linear SDMM framework by the following. The matrices $A \in \ensuremath{\mathbb{F}}_q^{t \times s}$ and $B \in \ensuremath{\mathbb{F}}_q^{s \times r}$ are partitioned into $p$ pieces with the inner product partitioning such that the order of the pieces of $B$ is reversed. The generator matrices $F$ and $G$ are defined as $(p + X) \times N$ Vandermonde matrices on the distinct evaluation points $\alpha_1, \dots, \alpha_N \in \ensuremath{\mathbb{F}}_q^\times$. These matrices generate Reed-Solomon codes of dimension $p + X$ and length $N$ on the locators $\alpha = (\alpha_1, \dots, \alpha_N)$. We denote this by $\mathcal{C}_A = \RS_{p + X}(\alpha)$ and $\mathcal{C}_B = \RS_{p + X}(\alpha)$. It is easy to see that this produces the same encoding as the general description of the secure MatDot scheme. It was shown in \cite{mirandola2015critical} that the resulting star product code is then $\mathcal{C}_A \star \mathcal{C}_B = \RS_{2p + 2X - 1}(\alpha)$, provided that $N \geq 2p + 2X - 1$.

\textit{GASP code:} The GASP code can similarly be described using the linear SDMM framework. The matrices $A \in \ensuremath{\mathbb{F}}_q^{t \times s}$ and $B \in \ensuremath{\mathbb{F}}_q^{s \times r}$ are partitioned into $m$ and $n$ pieces, respectively, using the outer product partitioning. The encoding is done using a polynomial evaluation code, where the terms of the polynomials are determined by optimizing the degree table. Let $\varphi$ and $\gamma$ be the vectors that contain the exponents of the encoding polynomials $f(x)$ and $g(x)$, respectively. The generator matrices are
\begin{equation*}
    F = \begin{pmatrix}
    \alpha_1^{\varphi_1} & \dots & \alpha_N^{\varphi_1} \\
    \vdots & \ddots & \vdots  \\
    \alpha_1^{\varphi_{m + X}} & \dots & \alpha_N^{\varphi_{m + X}}
    \end{pmatrix}, \quad G = \begin{pmatrix}
    \alpha_1^{\gamma_1} & \dots & \alpha_N^{\gamma_1} \\
    \vdots & \ddots & \vdots  \\
    \alpha_1^{\gamma_{n + X}} & \dots & \alpha_N^{\gamma_{n + X}}
    \end{pmatrix}.
\end{equation*}
Thus, $\mathcal{C}_A$ is an $[N, m + X]$ code and $\mathcal{C}_B$ is an $[N, n + X]$ code.

\textit{SDMM based on FFT:} The SDMM scheme based on FFT that was first presented in \cite{mital2020secure} uses the inner product partitioning to partition the matrices to $p = N - 2X$ pieces. The generator matrices can be expressed as
\begin{align*}
    F &= \begin{pmatrix}
    1 & 1 & \dots & 1 \\
    1 & \alpha_N & \dots & \alpha_N^{N - 1} \\
    1 & \alpha_N^2 & \dots & \alpha_N^{2(N - 1)} \\
    \vdots & & \ddots & \vdots \\
    1 & \alpha_N^{p + X - 1} & \dots & \alpha_N^{(p + X - 1)(N - 1)}
    \end{pmatrix}, \quad
    G = \begin{pmatrix}
    1 & 1 & \dots & 1 \\
    1 & \alpha_N^{-1} & \dots & \alpha_N^{-(N - 1)} \\
    \vdots & & \ddots & \vdots \\
    1 & \alpha_N^{-(p - 1)} & \dots & \alpha_N^{-(p - 1)(N - 1)} \\
    1 & \alpha_N^{-(p + X)} & \dots & \alpha_N^{-(p + X)(N - 1)} \\
    \vdots & & \ddots & \vdots \\
    1 & \alpha_N^{-(p + 2X - 1)} & \dots & \alpha_N^{-(p + 2X - 1)(N - 1)}
    \end{pmatrix}.
\end{align*}
These follow directly from the general description in Section~\ref{sec:examples_of_SDMM_schemes}.

\subsection{Security of Linear SDMM Schemes}\label{sec:security_linear_SDMM}

The following theorem will highlight the usefulness of the linear SDMM framework, since the security of the schemes can be proven by checking the properties of the generator matrices. A matrix has the MDS property if and only if it generates an MDS code, \emph{i.e.}, all of its maximal submatrices are invertible. This approach has been taken in \cite{d2020gasp} in reference to GASP codes. Here we apply the theorem in a more general setting.

\begin{theorem}\label{thm:linear_SDMM_security}
Let $F$ and $G$ be the generator matrices used in the encoding of some linear SDMM scheme. Then the scheme is secure against $X$-collusion if $F^{> m}$ and $G^{> n}$ have the MDS property.
\end{theorem}

\begin{proof}
Let $\mathcal{X} \subseteq [N]$, $|\mathcal{X}| = X$, be a set of $X$ colluding nodes. Writing the generator matrix $F$ as
\begin{equation*}
    F = \begin{pmatrix}
    F^{\leq m} \\
    F^{>m}
    \end{pmatrix}
\end{equation*}
allows us to write the shares the colluding nodes have about the encoded matrix $\widetilde{A}$ as
\begin{equation*}
    \widetilde{A}_\mathcal{X} = (A_1, \dots, A_m)F^{\leq m}_\mathcal{X} + (R_1, \dots, R_X)F^{> m}_\mathcal{X}.
\end{equation*}
If $F^{> m}$ has the MDS property, then any $X \times X$ submatrix of $F^{> m}$ is invertible. As $(R_1, \dots, R_X)$ is uniformly distributed, we get that $(R_1, \dots, R_X)F^{> m}_\mathcal{X}$ is also uniformly distributed. Therefore,
\begin{equation*}
    I(\mathbf{A}; \mathbf{\widetilde{A}}_\mathcal{X}) = 0.
\end{equation*}
The idea is that the confidential data of $A$ is hidden by adding uniformly random noise. A similar argument works for the matrix $B$. Finally, we get that
\begin{align*}
    I(\mathbf{A}, \mathbf{B}; \mathbf{\widetilde{A}}_\mathcal{X}, \mathbf{\widetilde{B}}_\mathcal{X}) &= I(\mathbf{A}, \mathbf{B}; \mathbf{\widetilde{A}}_\mathcal{X}) + I(\mathbf{A}, \mathbf{B}; \mathbf{\widetilde{B}}_\mathcal{X} \mid \mathbf{\widetilde{A}}_\mathcal{X}) \\
    &\leq I(\mathbf{A}; \mathbf{\widetilde{A}}_\mathcal{X}) + I(\mathbf{B}; \mathbf{\widetilde{B}}_\mathcal{X}) = 0.
\end{align*}
The inequality follows from $\mathbf{\widetilde{A}}_\mathcal{X}$ being conditionally independent of $\mathbf{B}$ given $\mathbf{A}$, and $\mathbf{\widetilde{B}}_\mathcal{X}$ being conditionally independent of $\mathbf{\widetilde{A}}_\mathcal{X}$ and $\mathbf{A}$ given $\mathbf{B}$. This proves the claim.
\end{proof}

The security of the SDMM schemes presented earlier follows as simple corollaries of the above theorem.

\section{Error Correction in SDMM}\label{sec:error_correction_SDMM}

Protecting against straggling workers has been a subject of research in many SDMM schemes. Another form of robustness is protection against so-called Byzantine workers, which return erroneous responses as a result of a fault or on purpose. This error can occur during the computation or during transmission. Robustness against Byzantine workers has been studied in the context of private information retrieval (PIR) and other distributed computation systems such as Lagrange coded computation in \cite{yu2019lagrange}.

The difference between straggling workers and Byzantine workers is that a straggling worker is simple to detect, while noticing erroneous responses from a Byzantine worker is not as straightforward. In coding-theoretic terms the straggling workers correspond to erasures in codes and Byzantine workers correspond to errors. It is well-known that erasures require one additional code symbol to fix with MDS codes, while errors typically require two additional code symbols to fix. The authors of \cite{yu2019lagrange} devised a coded computation scheme, where each additional straggler requires one additional response and each Byzantine worker requires two additional responses. This disparity between the costs can be fixed using interleaved codes by utilizing the structure of the error patterns.

\subsection{Interleaved Codes in SDMM}\label{sec:interleaved_codes_SDMM}

\begin{figure}[t]
    \centering
    \includegraphics[width=0.4\columnwidth]{drawing3x4.png}
    \caption{Diagram depicting the responses from the worker nodes. The Byzantine worker is depicted by the purple layer and the straggler by the blurred layer. Each response is a matrix, which is represented as a rectangular array in the figure. The codewords are the length $N$ vectors formed by stacking the responses and looking at the corresponding matrix entries. Hence, a Byzantine worker and stragglers can only affect their own position in the codewords.}
    \label{fig:response_diagram}
\end{figure}

The responses from the workers in a linear SDMM scheme can be seen to be codewords in the code $\mathcal{C}_A \star \mathcal{C}_B$. If this code has error correction capabilities, then we can correct possible errors in each of the codewords by considering each matrix position independently. This approach requires $N = K + 2B$ workers, where $K$ is the recovery threshold of the scheme and $B$ is the number of Byzantine workers. However, a more efficient method would be to consider all of the codewords at once, since we can use the structure of the errors to our advantage.

Figure \ref{fig:response_diagram} depicts how the responses of a linear SDMM scheme can be seen as a collection of codewords from the star product code $\mathcal{C}_A \star \mathcal{C}_B$. Each layer in the diagram depicts one of the responses from one of the workers. By collecting the matching matrix positions to a vector of length $N$ we obtain codewords in the code $\mathcal{C}_A \star \mathcal{C}_B$. If one of the workers returns an incorrect result, say worker 2 in Figure \ref{fig:response_diagram}, then the corresponding coordinate in each of the codewords is erroneous. Similarly, if one of the workers fails to return a response in time, say worker 4 in Figure \ref{fig:response_diagram}, then the corresponding coordinate is an erasure in each of the codewords. In particular, if we consider the $\ell$-interleaved code obtained by taking some collection of $\ell$ codewords in the response matrices, then the errors caused by the Byzantine workers are burst errors in the interleaved code.

There are several algorithms for decoding interleaved codes that can reach error correction capabilities that can correct up to twice as many errors as non-interleaved decoders \cite{schmidt2009collaborative, holzbaur2021success}. This is achieved by collaborative decoding, where the fact that the erroneous symbols are in the same place in each codeword is utilized.

In principle, interleaved coding techniques can be used with any linear SDMM scheme. However, many codes that are used in different SDMM constructions do not have an explicit error correction algorithm. SDMM schemes that are based on polynomial interpolation, such as the secure MatDot or $\mathrm{GASP}_\text{big}$ schemes, can be utilized, since the Reed--Solomon code has well-known error correction algorithms.

The proposed idea for correcting errors from the responses of a linear SDMM scheme is the following.
\begin{itemize}
    \item Choose some subset of $\ell$ matrix entries and collect the corresponding $\ell$ codewords of $\mathcal{C}_A \star \mathcal{C}_B$ as a codeword of an interleaved code.
    \item Find the error locations from the interleaved code using error correction algorithms for the $\ell$-interleaved code.
    \item Treat the erroneous coordinates as erasures and decode as usual.
\end{itemize}
As error correction of the interleaved codewords requires more computational power than just decoding, it is not advantageous to choose $\ell$ to be maximal, \emph{i.e.}, choosing all of the matrix entries to the interleaved codeword. On the other hand, collaborative decoding algorithms do not guarantee success with probability 1, so $\ell$ has to be chosen such that the success probability is suitably high. Furthermore, choosing a small subset of the matrix entries might lead to situations where not all of the errors are detected by the error correction algorithm.

\subsection{Analyzing Error Correction Capabilities}\label{sec:analyzing_error_correction}

\begin{figure}[t]
    \centering
    \includegraphics[width=0.55\columnwidth]{num_errors_vs_P_unsuccessful.png}
    \caption{Probability of unsuccessful decoding as a function of the number of errors $t\leq B$, where $q = 2^{32} - 5$, $D = 19$, $B = D - 2 = 17$ and $\ell = 9$. It is clear from the figure that the probability of unsuccessful decoding is extremely small.}
    \label{fig:number_of_errors}
\end{figure}

Collaborative error correction algorithms have been designed for interleaved Reed-Solomon codes, since they are prevalent in many applications where burst errors are common. In this section we analyze the success probability of some interleaved Reed--Solomon decoders in the context of the secure MatDot and $\mathrm{GASP}_\text{big}$ schemes. The same techniques are applicable for other linear SDMM schemes based on Reed--Solomon codes.

We assume that the errors sent by the Byzantine workers are uniformly distributed, \emph{i.e.}, the Byzantine workers respond with a matrix of appropriate size whose entries are uniformly distributed over the field $\ensuremath{\mathbb{F}}_q$. This is a natural assumption if the errors occur naturally without malice. Additionally, this assumption is popular in the literature, where failure probabilities are analyzed. In the examples in this section we use a field of size $q = 2^{32} - 5$, since it is large enough to use in problems, where the elements are discretized from the real numbers.

Bounded distance decoders for interleaved Reed--Solomon codes are discussed in \cite{schmidt2009collaborative, holzbaur2021success}. These decoding algorithms generalize the Berlekamp--Massey approach of decoding Reed--Solomon codes to interleaved codes. Additionally, \cite{schmidt2009collaborative, holzbaur2021success} give bounds on the success probability of the decoders when the errors are assumed to be uniformly distributed with specified column weights. Theorem 7 of \cite{holzbaur2021success} gives an upper bound on the probability of unsuccessful decoding as
\begin{align}
    \left( \frac{q^\ell - q^{-1}}{q^\ell - 1} \right)^t \cdot \frac{q^{-(\ell + 1)(t_\text{max} - t)}}{q - 1}, \label{eq:holzbaur_upper_bound}
\end{align}
where $t_\text{max} = \frac{\ell}{\ell + 1}(D - 1)$, $D$ is the minimum distance and $t$ is the number of errors.

To analyze the probability of unsuccessful decoding we use the upper bound in \eqref{eq:holzbaur_upper_bound}. This upper bound depends only on the field size $q$, the interleaving order $\ell$, the minimum distance $D$, and the number of Byzantine workers $t$, which is bounded from above by the maximal number of Byzantine workers $B$. The upper bound of unsuccessful decoding as a function of the number of Byzantine workers is depicted in Figure \ref{fig:number_of_errors}.

From Figure \ref{fig:number_of_errors} it is clear that the decoding is successful with a strikingly high probability. This shows that $N = K + B + 1$ workers are needed to correct errors from $B$ Byzantine workers, where the scheme has recovery threshold $K$. By adding $S$ straggling workers, we require an additional $S$ responses, which means that we need at least $N = K + S + B + 1$. This is an improvement over independent decoding, which requires $N = K + 2B + S$ workers, since the number of errors can be up to $\lfloor \frac{N - K}{2} \rfloor$.

\section*{Acknowledgements}

This work has been supported by the Academy of Finland under Grant No.\ 336005.