A Matrix Ensemble Kalman Filter-based Multi-arm Neural Network to
  Adequately Approximate Deep Neural Networks
Deep Learners (DLs) are the state-of-art predictive mechanism with
applications in many fields requiring complex high dimensional data processing.
Although conventional DLs get trained via gradient descent with
back-propagation, Kalman Filter (KF)-based techniques that do not need gradient
computation have been developed to approximate DLs. We propose a multi-arm
extension of a KF-based DL approximator that can mimic DL when the sample size
is too small to train a multi-arm DL. The proposed Matrix Ensemble Kalman
Filter-based multi-arm ANN (MEnKF-ANN) also performs explicit model stacking
that becomes relevant when the training sample has an unequal-size feature set.
Our proposed technique can approximate Long Short-term Memory (LSTM) Networks
and attach uncertainty to the predictions obtained from these LSTMs with
desirable coverage. We demonstrate how MEnKF-ANN can "adequately" approximate
an LSTM network trained to classify what carbohydrate substrates are digested
and utilized by a microbiome sample whose genomic sequences consist of
polysaccharide utilization loci (PULs) and their encoded genes.
Ved Piyush, Yuchen Yan, Yuzhen Zhou, Yanbin Yin, Souparno Ghosh
A Matrix Ensemble Kalman Filter-based Multi-arm Neural Network to Adequately Approximate Deep Neural Networks

###### Abstract

Deep Learners (DLs) are the state-of-art predictive mechanism with applications in many fields requiring complex high dimensional data processing. Although conventional DLs get trained via gradient descent with back-propagation, Kalman Filter (KF)-based techniques that do not need gradient computation have been developed to approximate DLs. We propose a multi-arm extension of a KF-based DL approximator that can mimic DL when the sample size is too small to train a multi-arm DL. The proposed Matrix Ensemble Kalman Filter-based multi-arm ANN (MEnKF-ANN) also performs explicit model stacking that becomes relevant when the training sample has an unequal-size feature set. Our proposed technique can approximate Long Short-term Memory (LSTM) Networks and attach uncertainty to the predictions obtained from these LSTMs with desirable coverage. We demonstrate how MEnKF-ANN can "adequately" approximate an LSTM network trained to classify what carbohydrate substrates are digested and utilized by a microbiome sample whose genomic sequences consist of polysaccharide utilization loci (PULs) and their encoded genes. The scripts to reproduce the results in this paper are available at [https://github.com/Ved-Piyush/MEnKF-ANN-PUL](https://github.com/Ved-Piyush/MEnKF-ANN-PUL).

## 1 Introduction

Deep Learners (DLs) have achieved state-of-art status in empirical predictive modeling in a wide array of fields. The ability of DLs to synthesize vast amounts of complex data, ranging from high dimensional vectors to functions and images, to produce accurate predictions has made them the go-to models in several areas where predictive accuracy is of paramount interest. Bioinformatics has also seen a steep increase in articles developing or deploying DL techniques in recent years
[Min et al., 2017]. However, conventional DLs trained via gradient descent with back-propagation require tuning of a large number of hyperparameters. Additionally, given the vast number of weights that DLs estimate, they are prone to overfitting when the training sample size is relatively small. Since the gradient descent algorithms compute the weights deterministically, DLs in their vanilla form do not yield any uncertainty estimate.

Several techniques have been proposed to alleviate the foregoing issues in DLs. For instance, Bayesian Neural Network (BNN) [Kononenko, 1989, Neal, 2012] was explicitly devised to incorporate epistemic and aleatoric uncertainty in the parameter estimation process by assigning suitable priors on the weights. The Bayesian mechanism can process these priors and generate uncertainty associated with DL predictions. Additionally, with judicious choice of priors, BNNs can be made less prone to overfitting [Fortuin et al., 2021, Srivastava et al., 2014]. Variational inference is another popular technique for uncertainty quantification in DLs. In particular, Hinton and Van Camp [1993] showed that a posterior distribution for the model weights could be obtained by minimizing the Kullback-Leibler distance between a variational approximation of the posterior and the true posterior of the weights. The Bayes by Backprop [Blundell et al., 2015] is another technique that uses variational formulation to extract uncertainty associated with the weights in DLs.

However, the Monte Carlo dropout technique [Srivastava et al., 2014], wherein each neuron (and all of its connections) is randomly dropped with some probability during the model training process, has arguably turned out to be the most popular method to regularize DLs and extract predictive uncertainty. In addition to the conceptual simplicity of the dropout technique, it was shown that the models trained using dropout are an approximation to Gaussian processes and are theoretically equivalent to variational inference [Gal and Ghahramani, 2016]. Regardless of its conceptual simplicity and theoretical underpinning, dropout methods require gradient computation. Hence, it can be quite computationally intensive in DLs with millions of parameters.

Another suite of methods for approximating DLs uses Kalman Filters (KF), or its variants, to obtain approximate estimates of the DL parameters [Yegenoglu et al., 2020, Rivas and Personnaz, 1998, Wan and Van Der Merwe, 2000, Julier and Uhlmann, 2004, Chen et al., 2019]. In particular, the Ensemble Kalman Filter (EnKF) technique offers a computationally fast approximation technique to DLs. For instance, Chen et al. [2019] train a single hidden layer neural network using the EnKF updating equations outlined in Iglesias et al. [2013] and show how using the augmented state variable, one can estimate the measurement error variance. In the DL setting, Chen et al. [2018] demonstrate the utility of EnKF in approximating a Long Short Term Memory (LSTM) model. Yegenoglu et al. [2020] use the EnKF to train a Convolutional Neural Network directly using the Kalman Filtering equations derived in Iglesias et al. [2013]. All these methods approximate single-arm DLs and, therefore, cannot be used in situations where input features can be represented in multiple ways. Our target is to develop a multi-arm approximator to the DL. In principle, we can train different DL approximators for each representation and perform a post-training model averaging. However, that would increase the computation cost substantially. We argue that since multiple different feature representations do not necessarily offer complimentary information, developing a multi-arm approximator that performs model averaging while training would perform "adequately". To that end, we develop a Matrix Ensemble Kalman Filter (MEnKF)-based multi-arm ANN that approximates a deep learner and simultaneously performs model averaging.

We apply our method to approximate an LSTM model trained to classify what carbohydrate substrates are digested and utilized by a microbiome sample characterized by genomic sequences consisting of polysaccharide utilization loci (PULs) [Bjursell et al., 2006] and their encoded genes. We use two different representations of the genomic sequences consisting of the PULs in two different arms, our MEnKF-ANN approximator, and demonstrate that our approximator closely
follows the predicted probabilities obtained from the trained LSTM. We also generate prediction intervals around the LSTM-predicted probabilities. Our results show that the average width of the prediction interval obtained from the MEnKF-ANN approximator is lower than that obtained from the original LSTM trained with MC dropout. We also perform extensive simulations, mimicking the focal dataset, to demonstrate that our method has desirable coverage for test samples compared to the MC dropout technique. Finally, we emphasize that even though the original problem is binary classification, our MEnKF-ANN approximator is designed to emulate the probabilities obtained from the original LSTM model and quantify the uncertainties in the LSTM-predicted probabilities.

The remainder of the article is organized as follows: In section 2, we describe the aspects of an experimentally obtained microbiome dataset that motivated us to design this approximator. In section 3, for the sake of completion, we offer a brief review of KF, EnKF, and how these techniques have been used to train DLs. Section 4 details the construction of our MEnKF-ANN method. In section 5, we offer extensive simulation results under different scenarios and follow it up with the application on real data in section 6. Finally, section 7 offers concluding remarks and future research directions.

## 2 Motivating Problem

The human gut, especially the colon, is a carbohydrate-rich environment (Kaoutari et al., 2013). However, most of the non-starch polysaccharides (for example, xylan, pectin, resistant glycans) reach the colon undegraded (Pudlo et al., 2022) because human digestive system does not produce the enzymes required to degrade these polysaccharides (Flint et al., 2012). Instead, humans have developed a complex symbiotic relationship with gut microbiota, with the latter providing a large set of enzymes for degrading the aforementioned non-digestible dietary components (Valdes et al., 2018). Consequently, an essential task in studying the human gut microbiome is to predict what carbohydrate substrates a microbiome sample can digest from the genetic characterization of the said microbiome (Koropatkin et al., 2012).

In order to generate a focused genetic characterization of the microbes that relates to their carbohydrate utilization property, one often investigates the genes encoding the Carbohydrate Active Enzymes (CAZymes) and other proteins that target glycosidic linkages and act to degrade, synthesize, or modify carbohydrates (Lombard et al., 2014; Zhang et al., 2018). This set of genes tend to form physically linked gene clusters in the genome known as polysaccharide utilization loci (PULs) (Bjursell et al., 2006). Consequently, the gene sequences associated with PULs of microbes could be used as a predictor to ascertain the carbohydrate substrate the microbe can efficiently degrade. However, these gene sequences are string-valued quantities (Huang et al., 2018; Stewart et al., 2018) and hence their naive quantitative representations (for instance, one-hot-encoding or count vectorization) often do not produce classifiers with acceptable accuracy (Badjatiya et al., 2017). Instead, we can use LSTM to process the entire sequence of string-valued features and then implement a classifier with a categorical loss function. The trained LSTM, by default, produces an embedding of the gene sequences in a vector space. Alternatively, we can also use a Doc2Vec embedding of the entire sequence associated with the PUL or an arithmetic average of Word2Vec embedding of each gene in the sequence and train a shallow learner (an ANN, for example). Since various representations of the predictors are available, we can train a multi-arm DL that takes different representations of features in different arms and performs concatenation/late integration of the embeddings before the prediction layer (Liu et al., 2020,
hanifi-Noghabi et al., 2019]. However, such multi-arm DLs require a relatively large number of training samples - typically tens of thousands.

Since the experimental characterization of new PULs for carbohydrate utilization is an expensive process [Ausland et al., 2021], we do not have large enough labeled samples to train complex multi-arm DLs. This predicament motivates us to develop a multi-arm approximator to a DL with the following capabilities: (a) it must "adequately" approximate the focal single-arm DL, (b) it should be able to ingest different feature representations in different arms and perform model averaging, (c) it should be able to detect if the set of representations supplied to it is substantially different from the representations used to train the original DL, i.e., sensitive to major misspecification. Since the original DL is trained on a single representation of features and the approximator is ingesting multiple representations, the latter is misspecified in a strict sense. However, this misspecification is deliberately introduced to circumvent training a multi-arm DL and assess, via model averaging, whether there is any benefit in using multiple representations of the same feature set.

We extract the dataset from the dbCAN-PUL database [Ausland et al., 2021] that contains experimentally verified PULs and the corresponding GenBank sequences of these PULs along with known target carbohydrate substrates. Figure 1 shows an example of a gene sequence associated with a PUL for the substrate Pectin. We have a total of approximately 411 data points. Figure 2 shows the dataset's frequency distribution of various target substrates. We do not have sufficient samples to train a complex DL to classify all the available substrates. Hence we propose to classify the two most frequently occurring target substrates - Xylan and Pectin - and train an LSTM binary classifier. Seventy-four samples belong to these two classes of substrates in a reasonably balanced way. One way to attach uncertainty to the probabilities estimated by the LSTM architecture is to activate the dropout layers during the prediction phase. This will generate multiple copies of the prediction depending on the architecture of the dropout layers. However, we need to decide on how many dropout layers to include and where to place them. Often the prediction intervals are pretty sensitive to the number and placements of dropout layer(s). For instance, the top left and bottom left panels of Figure 3 show the prediction intervals associated with eight held-out test samples obtained when two dropout layers were included - one inside the LSTM and one just before the final prediction layer. In contrast, the top right and bottom right panels of Figure 3 show the prediction intervals associated with the same test samples obtained when the second dropout layer was removed from the foregoing LSTM architecture. Observe how the number and placement of dropout layers influence the variability in the width of these intervals. If we wish to control the width variability, the placement of the dropout layer becomes a tuning parameter and further increases the hyperparameter search space.

We empirically show that our MEnKF-ANN approximator, trained on the logit transformed LSTM-estimated probabilities as the response and the embedding of the sequences obtained from LSTM and Doc2Vec operations as two types of features produce more stable prediction intervals regardless of the location of the dropout layer in the original LSTM.

Figure 1: Pectin PUL
Figure 3: Boxplots showing the Predictions superimposed with the ground truth value from Heavy and Low Just LSTMs

Figure 2: Frequency distribution for the various substrates
Background

This section offers a brief overview of Kalman Filters and Ensemble Kalman Filters and discusses how these methods have been used to train NNs and approximate DLs. For an extensive discussion on KF and EnKF technique, we direct the audience to Katzfuss et al. (2016) and Evensen (2003), respectively.

### Linear State Space Model & Kalman Filter

Consider a linear Gaussian state-space model given by

\[y_{t}=H_{t}x_{t}+\epsilon_{t},\ \ \epsilon_{t}\sim\mathcal{N}_{m_{t}}( 0,R_{t}) \tag{1}\] \[x_{t}=M_{t}x_{t-1}+\eta_{t},\ \ \eta_{t}\sim\mathcal{N}_{n}( 0,Q_{t}) \tag{2}\]

where \(y_{t}\) is the \(m_{t}\) dimensional observation vector at time step \(t\), \(x_{t}\) is the \(n\) dimensional state variable at that time, \(H_{t}\) and \(M_{t}\) denote the observation and the state transition matrices.

Assume that the filtering distribution of the state vector at \(t-1\) is given by

\[x_{t-1}|y_{1:t-1}\sim\mathcal{N}(\hat{\mu}_{t-1},\hat{\Sigma}_{t-1}), \tag{3}\]

KF computes the forecast distribution at \(t\) using (2) as

\[x_{t}|y_{1:t-1}\sim\mathcal{N}(\hat{\mu}_{t},\hat{\Sigma}_{t})\] \[\hat{\mu}_{t}:=M_{t}\hat{\mu}_{t-1},\] \[\tilde{\Sigma}_{t}:=M_{t}\hat{\Sigma}_{t-1}M_{t}^{{}^{\prime}}+Q_ {t} \tag{4}\]

Once the measurement at time step \(t\) becomes available, the joint distribution of \((x_{t},y_{t})\) is given by

\[\begin{pmatrix}x_{t}\\ y_{t}\end{pmatrix}\bigg{|}y_{1:t-1}\sim N\left(\begin{pmatrix}\tilde{\mu}_{t}\\ H_{t}\tilde{\mu}_{t}\end{pmatrix},\begin{pmatrix}\tilde{\Sigma}_{t}&\tilde{ \Sigma}_{t}H_{t}^{{}^{\prime}}\\ H_{t}\tilde{\Sigma}_{t}&H_{t}\tilde{\Sigma}_{t}H_{t}^{{}^{\prime}}+R_{t}\end{pmatrix}\right) \tag{5}\]

Then the updated filtering distribution is \(x_{t}|y_{1:t}\sim\mathcal{N}(\hat{\mu}_{t},\hat{\Sigma}_{t})\) where \(\hat{\mu}_{t}\) and \(\hat{\Sigma}_{t}\) are given by

\[\hat{\mu_{t}}:=\tilde{\mu_{t}}+K_{t}(y_{t}-H_{t}\tilde{\mu_{t}}),\] \[\hat{\Sigma_{t}}:=(I_{n}-K_{t}H_{t})\tilde{\Sigma}_{t}, \tag{6}\]

with \(K_{t}:=\tilde{\Sigma}_{t}H_{t}^{{}^{\prime}}(H_{t}\tilde{\Sigma}_{t}H_{t}^{{} ^{\prime}}+R_{t})^{-1}\) being the Kalman Gain Matrix. For large \(n\) and \(m_{t}\), computing the matrices in (6) is computationally expensive and often leads to numeric instability.

### Ensemble Kalman Filter

The idea of EnKF is to take an ensemble of size \(N\) from the filtering distribution at \(t-1\). This ensemble is denoted as \(\hat{x}_{t-1}^{(0)},\ \ \hat{x}_{t-1}^{(1)},...,\hat{x}_{t-1}^{(N)}\sim \mathcal{N}_{n}\left(\hat{\mu}_{t-1},\hat{\Sigma}_{t-1}\right)\). In the forecast step of EnKF, (2)
is applied to the ensemble members to obtain their evolution from \(t-1\) to \(t\). That is

\[\tilde{x}_{t}^{(i)}=M_{t}\hat{x}_{t-1}^{(i)}+\eta_{t}^{(i)},\ \ \eta_{t}^{(i)}\sim \mathcal{N}(0,Q_{t}),\ \ \ \ i=1,\ldots,N \tag{7}\]

It can also be shown that \(\tilde{x}_{t}^{(i)}\sim\mathcal{N}(\tilde{\mu}_{t},\tilde{\Sigma}_{t})\). Similar to the update step of the Kalman Filter, all the members of this ensemble must be updated when the measurement at time step t becomes available. To update these ensemble members, first, a sample is obtained for the measurement error sequences that is \(\epsilon_{t}^{(1)},\epsilon_{t}^{(2)},...,\epsilon_{t}^{(N)}\sim\mathcal{N}_{ m_{t}}(0,R_{t})\). Then using these simulated measurement errors \(N\) perturbed observations \(\tilde{y}_{t}^{(1)},\tilde{y}_{t}^{(2)},\ldots,\tilde{y}_{t}^{(N)}\) are obtained using \(\tilde{y}_{t}^{(i)}=H_{t}\tilde{x}_{t}^{(i)}+\epsilon_{t}^{(i)}\). Since the joint distribution of \((\tilde{x}_{t}^{(i)}\), \(\tilde{y}_{t}^{(i)})\) is the same as in (5), the updating equations are obtained by shifting the forecasted ensemble in (7) as follows

\[\hat{x}_{t}^{(i)}=\tilde{x}_{t}^{(i)}+K_{t}(y_{t}-\tilde{y}_{t}^{(i)}),\ \ \ \ i=1, \ldots,N \tag{8}\]

It can be easily shown that \(\hat{x}_{t}^{(i)}\sim\mathcal{N}_{n}(\hat{\mu}_{t},\hat{\Sigma}_{t})\). The computational benefit comes from the fact that instead of computing the Kalman gain matrix in (8) explicitly, the sample covariance matrix of the forecasted ensemble (\(\tilde{S}_{t}\), say) is used to estimate the Kalman Gain matrix as \(\hat{K}_{t}:=\tilde{S}_{t}H_{t}^{{}^{\prime}}(H_{t}\tilde{S}_{t}H_{t}^{{}^{ \prime}}+R_{t})^{-1}\).

### KF and EnKF for Deep Learners

Although conventional KF is only suitable for estimating parameters in linear state-space models, several extensions have been proposed to generalize KF in nonlinear settings. For instance, Rivals and Personnaz (1998) used extended KF to train feed-forward NN. Wan and Van Der Merwe (2000) introduced the unscented KF that better approximates nonlinear systems while making it amenable to the KF framework. Anderson (2001) used the concept of state augmentation that offered a generic method to handle nonlinearity in state-space models via the KF framework. Iglesias et al. (2013) utilized this state augmentation technique to develop a generic method to train ANNs. They derived the state-augmented KF's forecast and updated equations in ANNs, thereby providing the algebraic framework to train DLs using the Ensemble Kalman Filters approach. These equations were subsequently used by Yegenoglu et al. (2020) to train a Convolutional Neural Network using EnKF. Furthermore, Chen et al. (2019) also used the updating equations in Iglesias et al. (2013) to train a single hidden layer ANN and demonstrated how using state augmentation one can estimate the measurement error variance. State-augmented EnKF formulation also estimated parameters in LSTMs (Chen et al., 2018).

All the foregoing models offer techniques to estimate parameters of complex nonlinear DLs using the EnKF framework. However, they are unsuitable when we have multiple feature representations. We want to approximate a DL with a multi-arm ANN trained via EnKF, as discussed in section 2.

## 4 Methodology

First, we offer a generic construction of the proposed MEnKF-ANN procedure and describe how this method could be deployed to solve the problem in section 2. We will use the following notations. \(Y\in\mathcal{R}\) is our target response. We have a total of \(m=\sum_{t=1}^{T}m_{t}\) training instances, with \(m_{t}\) being the number of training data points in the \(t^{th}\) batch. \(v_{t}^{f}\in\mathcal{R}^{p}\) and \(v_{t}^{g}\in\mathcal{R}^{q}\) denote two different representations of the features (possibly of different dimensions) for the \(t^{th}\) batch of data.
Consider two ANNs, denoted by \(f\) and \(g\). Assume that the neural networks f and g have \(n_{f}\), \(n_{g}\) number of learnable parameters. For illustrative purposes, we will assume \(n_{f}=n_{g}\). If \(n_{f}\neq n_{g}\), we can use suitable padding when updating the weights. In the \(t^{th}\) batch of data, we assign the feature sets \(v^{f}_{t}\) and \(v^{g}_{t}\) to networks \(f\) and \(g\), respectively. We denote \(w^{f}_{t}\) and \(w^{g}_{t}\) to be the updated weights for the neural network \(f\) and \(g\), respectively.

### Matrix Kalman Filter based Multi-arm ANN

Consider the state matrix, \(X_{t}\), associated with the \(t^{th}\) batch of data given by

\[X_{t}^{(m_{t}+n_{g}+1)\times 2}=\begin{bmatrix}f(v^{f}_{t},w^{f}_{t}),&g(v^{g}_{ t},w^{g}_{t})\\ w^{f}_{t},&w^{g}_{t}\\ 0,&a_{t}\end{bmatrix} \tag{9}\]

where \(a_{t}\) and \(b_{t}\) are real-valued scalar parameters. Define \(H_{t}^{m_{t}\times(m_{t}+n_{g}+1)}=[I_{m_{t}},0_{m_{t}\times(n_{g}+1)}]\) and \(G_{t}^{2\times 1}=[1-\sigma(a_{t}),\ \ \sigma(a_{t})]^{T}\) where \(\sigma(.):\mathcal{R}\rightarrow[0,1]\), with the sigmoid function being a popular choice of \(\sigma(.)\). Additionally, define \(\Theta_{t-1}=I_{m_{t}+n_{g}+1}\) and \(\psi_{t-1}=I_{2}\). We are now in a position to define the Matrix Kalman Filter. The measurement equation is given by:

\[Y_{t}=H_{t}X_{t}G_{t}+\epsilon_{t} \tag{10}\]

with the state evolution equation being

\[X_{t}=\Theta_{t-1}X_{t-1}\psi_{t-1}+\eta_{t} \tag{11}\]

Writing in \(vec\) format, (11) becomes

\[x_{t}=vec(X_{t})=(\psi^{T}_{t-1}\otimes\Theta_{t-1})vec(X_{t-1})+vec(\eta_{t}) \tag{12}\]

Now letting \(\phi_{t-1}=\psi^{T}_{t-1}\otimes\Theta_{t-1}\) and \(\tilde{\eta}_{t}=vec(\eta_{t})\) we get from (12)

\[x_{t}=\phi_{t-1}x_{t-1}+\tilde{\eta}_{t} \tag{13}\]

(10) can similarly be compactified as

\[y_{t}=\mathcal{H}_{t}x_{t}+\epsilon_{t} \tag{14}\]

where \(\mathcal{H}_{t}=G_{t}^{T}\otimes H_{t}\). Observe that(14) and (13) have the same form as the standard representation of linear state space model described in (1) and (2). Therefore, we can get the matrix state space model's solution by converting it to the vector state space model and then using EnKF to approximate the updating equations. We direct the audience to [Choukroun et al., 2006] for more details on Matrix Kalman Filters.

### Interpreting MEnKF-ANN and a Reparametrization

The above construction of \(X_{t}\), \(H_{t}\), and \(G_{t}\) performs automatic model averaging while training. First, consider the matrix multiplication of \(H_{t}X_{t}\) from (10). This would be a \(m_{t}\times 2\) dimensional matrix in which the first column is the prediction, for the \(t^{th}\) batch, from the neural network \(f\) and the second column is the prediction from the neural network \(g\). Post multiplication by \(G_{t}\) would
take the weighted average of each row in \(H_{t}X_{t}\) where the weights are defined inside the \(G_{t}\) matrix. Now consider the matrix multiplication of \(H_{t}X_{t}G_{t}\) from (10)

\[H_{t}X_{t}G_{t} = \left[f(v_{t}^{f},w_{t}^{f}),\ \ g(v_{t}^{g},w_{t}^{g})\right]\begin{bmatrix}1- \sigma(a_{t})\\ \sigma(a_{t})\end{bmatrix} \tag{15}\] \[= \left[(1-\sigma(a_{t}))f(v_{t}^{f},w_{t}^{f})+\sigma(a_{t})g(v_{ t}^{g},w_{t}^{g})\right]\]

(15) clearly demonstrates how our construction explicitly performs model averaging across the batches with \(1-\sigma(a_{t})\) and \(\sigma(a_{t})\) being the convex weights allocated to the ANNs \(f\) and \(g\), respectively.

Although the foregoing construction connects Matrix KF formulation with multi-arm ANN and performs explicit model averaging, it suffers from a computational bottleneck. Using (13) and (14) the estimated Kalman Gain Matrix would be \(K_{t}=\widetilde{S}_{t}\mathcal{H}_{t}^{T}(\mathcal{H}_{t}\widetilde{S}_{t} \mathcal{H}_{t}^{T}+\sigma_{y}^{2}I_{m_{t}})^{-1}\). However, in the above parameterization we have \(G_{t}=[1-\sigma(a_{t}),\ \ \sigma(a_{t})]^{T}\) and \(\mathcal{H}_{t}=G_{t}^{T}\otimes H_{t}\). This would require computation of the estimated Kalman Gain matrix for each member in EnKF since, at any given iteration of our MEnKF-ANN, we have an \(a_{t}\) for each member of the ensemble. Thus computation complexity associated with the Kalman Gain computation increases linearly with the size of the ensemble in the above parametrization of the MEnKF-ANN.

To alleviate this computational bottleneck, consider the following parametrization:

\[X_{t}=\begin{bmatrix}(1-\sigma(a_{t}))f(v_{t}^{f},w_{t}^{f}),\ \sigma(a_{t})g(v_{t}^{g},w_{t}^{g})\\ w_{t}^{f},\ \ w_{t}^{g}\\ 0,\ \ a_{t}\end{bmatrix} \tag{16}\]

and \(G_{t}=[1,\ \ 1]^{T}\). We still have explicit model averaging in the measurement equation, i.e.,

\[H_{t}X_{t}G_{t}=\left[(1-\sigma(a_{t}))f(v_{t}^{f},w_{t}^{f})+\sigma(a_{t})g(v _{t}^{g},w_{t}^{g})\right] \tag{17}\]

but \(\mathcal{H}_{t}\) does not depend on \(a_{t}\). Therefore the matrix products for the Kalman Gain computation can now be computed once for each batch.

Turning to the variance parameter in the measurement equation (14). Assume \(\epsilon_{t}\sim\mathcal{N}_{m_{t}}(0,\nu_{y}^{2}I_{m_{t}})\) To estimate \(\nu_{y}^{2}\), we augment the state vector as follows:

\[X_{t}^{(m_{t}+n_{g}+2)\times 2}=\begin{bmatrix}(1-\sigma(a_{t}))f(v_{t}^{f},w_{t} ^{f}),\ \sigma(a_{t})g(v_{t}^{g},w_{t}^{g})\\ w_{t}^{f},\ \ w_{t}^{g}\\ 0,\ \ a_{t}\\ 0,\ \ b_{t}\end{bmatrix} \tag{18}\]

where \(\nu_{y}^{2}=\log(1+e^{b_{t}})\) and \(H_{t}\) in (10) now becomes \([I_{m_{t}},0_{m_{t}\times(n_{g}+2)}]\). We used a softplus transformation of \(\nu_{y}^{2}\), instead of the usual log transformation for computational stability.

### Connecting MEnKF-ANN with DL

Recall that our dataset consists of string-valued gene sequences associated with experimentally determined PULs, with the response being the carbohydrate substrates utilized by the said microbe. Since we consider only two categories of PULs, we have a binary classification problem.
An LSTM trained with a binary cross-entropy loss is the approximand DL in our case. Suppose \(p\) is the probability of observing a sample of a particular category. In that case, the trained LSTM produces \(\hat{p}\) for each training instance, along with an embedding of the associated gene sequences. Our MEnKF-ANN approximator uses \(logit(\hat{p})\) as the target response. The LSTM embedding of the gene sequences is fed into one arm of the approximator while the other arm ingests Doc2Vec encoding of the gene sequences. Thus, our MEnKF-ANN approximates the probabilities estimated by an LSTM. The convex weights \(\sigma(a)\) ascertain which embedding has more predictive power. Clearly, MEnKF-ANN operates as a model stacker, and the predictive uncertainty interval that it produces, by default, around its target approximand quantifies how well simpler ANNs, fitted without backpropagation, can approximate a deep learner.

To initialize the ensemble in the approximator, we draw the members in the state vector (18) from \(\mathcal{N}_{2(m_{t}+n_{g}+2)}(\mathbf{0},\nu_{x}^{2}I)\), where \(\nu_{x}^{2}\) is a tuning parameter that plays a key role in controlling the spread of the ensemble members and the dimension of \(I\) matches with the dimension of normal distribution. Following Chen et al. (2018, 2019), we assume the state transition is deterministic, i.e., \(x_{t}=\phi_{t-1}x_{t-1}\) and hence we do not have the variance parameter corresponding to \(\tilde{\eta}\) in the augmented state vector. When we reach the \(t^{th}\) batch of data, for the \(i^{th}\) member in the ensemble (\(i=1,2,...,N\)), we update each element in the augmented state vector \(w_{t}^{f,(i)}\), \(w_{t}^{g,(i)}\), \(a_{t}^{(i)}\), \(b_{t}^{(i)}\) using the updating equation (8) suitably modified to handle deterministic state transition.

## 5 Simulations

We conducted extensive simulations to assess how well our MEnKF-ANN can approximate an LSTM binary classifier. This simulation exercise aims to demonstrate that our MEnKF-ANN is not only "adequate" in approximating the probabilities produced by LSTM but can also capture the "true" probabilities that generate binary labels. We compute the coverage and width of the prediction intervals of the target probabilities in the test set to assess the "adequacy" of the approximator. Then, we compare this coverage and width with those computed directly via an LSTM trained with MC dropout. Admittedly, the prediction intervals obtained from the latter are different from those computed from MEnKF-ANN. However, if the ground truth probabilities are known, an adequate approximator should be able to achieve near-nominal coverage when the approximand is not misspecified.

Our simulation strategy mimics the focal dataset and uses the gene sequences associated with the original PULs to generate labels. As mentioned above, we extracted \(\hat{p}\) from the LSTM trained on the original dbCAN-PUL data. We call this LSTM the _true LSTM_. We consider \(\hat{p}\) the true probabilities for synthetic data generation. We then use noisy copies of \(\hat{p}\) to generate a synthetic label in the following way: generate \(logit(\tilde{p}_{i}^{(j)})=logit(\hat{p}_{i})+\epsilon_{i}^{*(j)},\ i=1,2,...,m,j=1,2,...,J\), where \(J\) is the number of the simulated dataset and \(m\), is the number of data points in each simulated set, the perturbation \(\epsilon_{i}^{*(j)}\) are iid Normal(0,0.01\({}^{2}\)). We generate synthetic labels \(\tilde{Y}\) by thresholding \(\tilde{p}_{i}^{(j)}\) at 0.5, i.e \(\tilde{Y}_{i}^{(j)}=I(\tilde{p}_{i}^{(j)}>0.5)\). Then the simulated dataset consists of \(D^{(j)}=\{\mathbf{F},\tilde{Y}^{(j)},\ j=1,2,...,J\}\), where \(\mathbf{F}\) is the set of original gene sequences from dbCAN-PUL.

Now in each \(D^{(j)}\), we train a second LSTM (with two dropout layers) and extract \(\tilde{\tilde{p}}_{i}^{(j)},i=1,2,...,m\) and the embedding of the gene sequences. We call these LSTMs, trained on \(D^{(j)}\), the _fitted LSTMs_. Note that the embeddings from _fitted LSTMs_ could potentially be different from those obtained from the _true LSTM_. We denote the embedding from _fitted LSTMs_ by \(v_{i}^{(j),f},\ j=1,2,...J\). Our MEnKF-ANN is constructed to approximate the _fitted LSTMs_. To that end, the approximator uses \(logit(\tilde{\tilde{p}}_{i}^{(j)})\) as the target response. \(v_{i}^{(j),f}\) are supplied as features to one arm of the ANN, the
other arm ingests \(v_{i}^{(j),g}\) - the Doc2Vec embedding of \(\mathbf{F}\). Once the MEnKF-ANN is trained, we use a hold-out set in each simulated dataset to generate predictive probabilities from the forecast distribution for each member in the KF ensemble and compute the empirical 95% predictive interval at \(logit^{-1}\) scale. To measure the adequacy of MEnKF-ANN, we compute the proportion of times the foregoing predictive interval contains \(\hat{p}\) in held-out test data. We expect this coverage to be close to the nominal 95%, and the average width of these intervals should not be greater than 0.5. Additionally, observe that the data-generating model uses LSTM embedding of \(F\); hence, using Doc2Vec embedding as input is misspecification. Consequently, we expect the average model weight associated with \(v^{f}\) to be larger than \(v^{g}\). Table 1 shows the performance of MEnKF-ANN in terms of coverage, the average width of prediction intervals, and average LSTM weight under two specifications of ensemble size (\(N\)) and initial ensemble variance (\(\nu_{x}^{2}\)). To compare these results, we offer the coverage and average width of the prediction intervals when both the dropout layers are activated in the _fitted LSTM_ during the prediction phase in Table 2. Observe how MEnKF-ANN recovered the _true probabilities_ even better than the correctly specified LSTM with dropout. The average interval widths obtained from MEnKF-ANN are also lower than those from the _fitted LSTM_. These demonstrate the adequacy of MEnKF-ANN in approximating the target DL. Additionally, we observe that the average LSTM model weight is \(\approx 1\) indicating the ability of our approximator to identify the correctly specified data-generating model. Figure 4 shows the histogram of the predictive samples obtained from the ensemble members for eight test samples in a randomly chosen replicate. The red vertical line denotes the true logits, and the green vertical lines show the fences of the 95% prediction interval.

Now, to demonstrate a situation where MEnKF-ANN is "inadequate", we supply the approximator with a completely different feature set representation. Instead of using the LSTM embedding \(v^{f}\), we use word2vec embedding of each gene in the predictor string and take the arithmetic average of these word2vec embeddings to represent the entire sequence. We denote this feature set by \(\tilde{v}^{f}\) and then train the MEnKF-ANN using \(\tilde{v}^{f}\) and \(v^{g}\) as the features and \(logit(\tilde{\tilde{p}}^{(j)})\) as the target response. Evidently, MEnKF-ANN is highly misspecified. Table 3 reports the coverage and average width of the prediction interval obtained from this model. Observing the huge width of the intervals essentially invalidates the point prediction. Such large width indicates that MEnKF-ANN may not approximate the target DL. Therefore, we caution against using the coverage and width metrics to assess the "adequacy" of the _fitted LSTM_ itself.

\begin{table}
\begin{tabular}{l l l l l} \hline \(N\) & \(\nu_{x}^{2}\) & Coverage & Width & LSTM weight \\ \hline
216 & 16 & 90.25\% & 0.33 & 0.9997 \\
216 & 32 & 89.25\% & 0.32 & 0.9999 \\ \hline \end{tabular}
\end{table}
Table 1: Performance of MEnKF-ANN using LSTM embedding and Doc2Vec

\begin{table}
\begin{tabular}{l l l l} \hline Rate & Reps & Coverage & Width \\ \hline
0.5 & 50 & 81.25\% & 0.53 \\
0.5 & 200 & 84.50\% & 0.56 \\ \hline \end{tabular}
\end{table}
Table 2: Coverage and width of prediction intervals obtained from _fitted LSTM_ with two dropout layers
## 6 Application

Recall that our focal dataset consists of \(n=74\) samples belonging to Xylan and Pectin. However, training an LSTM on a small sample size would require aggressive regularization, even with this reduced label space. Therefore, we draw on an extensive collection of unlabelled data containing gene sequences associated with CAZyme gene clusters (CGC) computationally predicted from genomic data [20, 17]. Although this unlabelled data contain approximately 250K CGC gene sequences, unlike experimentally characterized PULs, these sequences do not have known carbohydrate substrate information. They hence cannot be directly used for classification purposes. We, therefore, use this unlabelled dataset to learn the word2vec embeddings of each gene appearing in the unlabelled dataset. These embeddings are then used to initialize the embedding layer of the target LSTM classifier.

Turning to the labeled dataset, instead of performing full cross-validation, we resort to subsampling procedure [13]. We take a subsample of sixty-six instances for training and hold eight instances for testing purposes. The subsample size (\(b\)) is chosen such that \(b(n)/n\approx 8\sqrt{n}/n\to 0\), as \(n\rightarrow\infty\). Although the subsampling theory requires generating \(\begin{pmatrix}n\\ b\end{pmatrix}\) replicates, the computational cost for generating \(\approx 10^{11}\) replicates, in our case, is prohibitive. Instead, we generate 50 independently subsampled replicates comprising training and testing sets of sizes 66 and 8, respectively. In each replication, an LSTM with the same architecture is trained on the foregoing 66 training instances. Under this scheme, the probability that the \(i^{th}\) instance in our dataset appears at least once in the test set is \(\approx 99.6\%\).

The LSTM-estimated probabilities of observing a _Pectin_ substrate are extracted from each replicate. These probabilities are logit transformed and used as the target response for our MEnKF-ANN approximator. We feed the LSTM embedding and Doc2Vec embedding of the gene sequences

Figure 4: True Logits Superimposed on Predicted Logits from MEnKF-ANN using LSTM and Doc2Vec embedding

\begin{table}
\begin{tabular}{l l l l l} \hline \(N\) & \(\nu_{x}^{2}\) & Coverage & Width & Word2Vec weight \\ \hline
216 & 16 & 96.25\% & 0.83 & 0.9155 \\
216 & 32 & 94.25\% & 0.84 & 0.9787 \\ \hline \end{tabular}
\end{table}
Table 3: Performance of MEnKF-ANN using Word2Vec and Doc2Vec
into the two arms of the approximator along with the foregoing logit-transformed estimated probabilities. We then generate predictions on the held-out test data points in each replicate. Finally, we compare the LSTM-prediction of probabilities with those generated by MEnKF-ANN generated predictions. The average MAE and the proportion of times a 95% prediction interval contains the LSTM-generated predictions in the held-out data set, under two different MEnKF-ANN hyperparameter choices are shown in Table 4 indicating that our approximator can be adequately used to generate the predictions. We do not report the LSTM weights estimated by MEnKF-ANN because, as we observed in the simulation (Table 1), the approximator overwhelmingly prefers the LSTM embeddings. Figure 5 shows the scatter plot of MEnKF-ANN-predicted and LSTM-predicted probabilities for the held-out data across 50 replicates. Figure 6 shows the boxplots associated with MEnKF-ANN predictions for the same set of test samples for which LSTM generated prediction boxplots were shown in the left column of Figure 3. Evidently, MEnKF-ANN can adequately approximate the target LSTM.

Turning to the stability of prediction intervals, Table 5 shows the average width of the 95% prediction intervals obtained under two configurations of LSTM and their respective MEnKF-ANN approximators. LSTM\({}_{1}\) has two dropout layers (one in the LSTM layer and one before the final prediction layer) with a 50% dropout rate and 200 replicates. LSTM\({}_{2}\) has one dropout layer (in the LSTM layer) with a 50% dropout rate and 200 replicates. MEnKF-ANN\({}_{11}\) approximates LSTM\({}_{1}\) with 216 ensemble members and \(\nu_{x}^{2}=16\), MEnKF-ANN\({}_{12}\) also approximates LSTM\({}_{1}\), but now with 216 ensemble members and \(\nu_{x}^{2}=32\). Similarly, MEnKF-ANN\({}_{21}\) and MEnKF-ANN\({}_{22}\) approximates LSTM\({}_{2}\) with 216 ensemble members and \(\nu_{x}^{2}=16\) and \(\nu_{x}^{2}=32\), respectively. Observe that the variation in the average width between LSTM\({}_{1}\) and LSTM\({}_{2}\) is considerably higher than the variation between MEnKF-ANN\({}_{11}\) and MEnKF-ANN\({}_{21}\) or between MEnKF-ANN\({}_{12}\) and MEnKF-ANN\({}_{22}\). This indicates that the approximator produces more stable prediction intervals than obtaining prediction by activating the dropout layer during prediction.

Finally, we demonstrate how MEnKF-ANN can naturally handle two predictive models with potentially different feature sets. This situation is relevant because, owing to the small sample size, we can train a shallow learner (ANN with backpropagation, for instance) that takes Doc2Vec representation of gene sequences as predictors to estimate the probabilities of observing the _Pectin_ substrate. Now, we can average the probabilities estimated by the LSTM (\(\hat{p}_{LSTM}\), say) and ANN (\(\hat{p}_{ANN}\), say) to produce a model-averaged estimated probability of observing _Pectin_ (\(\hat{\bar{p}}\), say). However, how would we attach uncertainty to \(\hat{\bar{p}}\)? The multi-arm construction of MEnKF-ANN provides a natural solution in this situation. We supply, as described in the foregoing sections, LSTM embeddings and Doc2Vec embeddings to the two arms of MEnKF-ANN but use \(logit(\hat{\bar{p}})\) as the target response here. Thus MEnKF-ANN is now approximating the average output of two primary models. These primary models are trained on the same response variable but use two different representations of features. Table 6 shows the performance of MEnKF-ANN in this situation for some combinations of \(N\) and \(\nu_{x}^{2}\). The coverage is measured with respect to \(\hat{\bar{p}}\) on the test sets. Although the average width and MAE are larger than those reported in Table 4, we observe that the LSTM weights \(\approx\) 0.5, which is what we would expect because MEnKF-ANN is _seeing_ equally weighted outputs from LSTM and ANN.

## 7 Discussion

State-augmented Kalman Filter and its variants provide a gradient-free method that can be extended to approximate popular neural network-based deep learners for regression and classification
tasks. In this article, we have developed a Matrix Ensemble Kalman Filter-based multi-arm neural network to approximate an LSTM. We have demonstrated that this technique adequately approximates the target DL regarding coverage and the average width of the prediction interval. We have also demonstrated how the in-built model averaging capability can be leveraged to attach

\begin{table}
\begin{tabular}{l l l l l l} \hline \hline Target model & Average Width & Approximator & Average Width \\ \hline LSTM\({}_{1}\) & 0.492 & MEnKF-ANN\({}_{11}\) & 0.102 \\  & & MEnKF-ANN\({}_{12}\) & 0.085 \\ LSTM\({}_{2}\) & 0.371 & MEnKF-ANN\({}_{21}\) & 0.119 \\  & & MEnKF-ANN\({}_{22}\) & 0.108 \\ \hline \hline \end{tabular}
\end{table}
Table 6: Performance of MEnKF-ANN trained on the averaged probability of LSTM and shallow ANN

\begin{table}
\begin{tabular}{l l l l l l} \hline \hline \(N\) & \(\nu_{x}^{2}\) & Coverage & Width & MAE & CPU Time \\ \hline
216 & 16 & 90.50\% & 0.1024 & 0.0200 & 2.39 mins \\
216 & 32 & 85.50\% & 0.0850 & 0.0161 & 3.67 mins \\ \hline \hline \end{tabular}
\end{table}
Table 4: Performance of MEnKF-ANN using LSTM embedding and Doc2Vec for dbCAN-PUL data

Figure 5: Scatterplot of First LSTM Predicted Probabilities vs. EnKF Predicted Probabilities

\begin{table}
\begin{tabular}{l l l l l l} \hline \hline Target model & Average Width & Approximator & Average Width \\ \hline LSTM\({}_{1}\) & 0.492 & MEnKF-ANN\({}_{11}\) & 0.102 \\  & & MEnKF-ANN\({}_{12}\) & 0.085 \\ LSTM\({}_{2}\) & 0.371 & MEnKF-ANN\({}_{21}\) & 0.119 \\  & & MEnKF-ANN\({}_{22}\) & 0.108 \\ \hline \hline \end{tabular}
\end{table}
Table 5: Comparison of the average width of prediction interval LSTM + MC Dropout and MEnKF-ANN approximator for each LSTM
uncertainty to the averaged predictions generated by two different models. Our simulations suggest that by using an explicit model averaging construction, our approximator can also identify its target approximand. We have also observed that the prediction intervals generated by the approximator are less sensitive to the location of dropout layers and hence provide more stable prediction intervals than obtaining predictions by activating the dropout layers within the DL itself. Admittedly, our procedure requires an additional round of training, but its fast computation time (see Table 4), along with its ability to emulate the approximand, adequately compensate for that. We have also deployed our approximator on a highly accessed database, dbCAN-PUL, to attach uncertainty to the predicted probabilities produced by (a) the primary LSTM model and (b) an ensemble of LSTM and ANN models. The primary LSTM and ANN models were trained to classify two carbohydrate substrates using the gene sequences characterized by the PULs of the gut microbiome. We anticipate this technique will be helpful to domain experts in assessing the reliability of predictions generated by deep learners or an ensemble of learners.

In the future, we propose to expand our model to handle more than two classes. This would enable us to utilize the information in the dbCAN-PUL database better. Another possible direction is to develop an analog of MEnKF-ANN that can directly handle binary data. Although the KF technique crucially requires Gaussianity assumption, but Fasano et al. (2021) recently developed an extension of the KF method that can handle binary responses. We are actively investigating

Figure 6: Boxplots showing the MEnKF-ANN Predictions superimposed with the ground truth value for Heavy and Low Dropout
how this technique can be adapted to our MEnKF-ANN framework.

## 8 Competing interests

No competing interest is declared.

## 9 Author contributions statement

V.P., Y.Z., and S.G. conceived the models and experiment(s), and V.P. conducted the experiment(s). V.P. and S.G. analyzed the results. Y.Yan and Y.Yin contributed to the training data. V.P. drafted the manuscript. All authors reviewed the manuscript. Y. Yin secured the funding.

## 10 Acknowledgments

The authors thank the anonymous reviewers for their valuable suggestions. This work is supported in part by funds from the National Institute of Health (NIH: R01GM140370, R21AI171952) and the National Science Foundation (NSF: CCF2007418, DBI-1933521). In addition, we thank the lab members for their helpful discussions. This work was partially completed utilizing the Holland Computing Center of the University of Nebraska-Lincoln, which receives support from the Nebraska Research Initiative.
