Title: DEEP EXPECTATION-MAXIMIZATION IN HIDDEN MARKOV MODELS VIA SIMULTANEOUS PERTURBATION STOCHASTIC APPROXIMATION

Abstract: We propose a novel method to estimate the parameters of a collection of Hidden Markov Models (HMM), each of which corresponds to a set of known features. The observation sequence of an individual HMM is noisy and/or insufficient, making parameter estimation solely based on its corresponding observation sequence a challenging problem. The key idea is to combine the classical Expectation-Maximization (EM) algorithm with a neural network, while these two are jointly trained in an end-to-end fashion, mapping the HMM features to its parameters and effectively fusing the information across different HMMs. In order to address the numerical difficulty in computing the gradient of the EM iteration, simultaneous perturbation stochastic approximation (SPSA) is employed to approximate the gradient. We also provide a rigorous proof that the approximated gradient due to SPSA converges to the true gradient almost surely. The efficacy of the proposed method is demonstrated on synthetic data as well as a real-world e-Commerce dataset.

1. INTRODUCTION:

Hidden Markov Model (HMM) is a powerful statistical framework in sequential data analysis. One advantage of HMM lies in its strong model interpretability (Rabiner, 1989) . Estimating the parameters of a HMM given its observation sequence is a well-studied problem. The celebrated Expectation-Maximization (EM) algorithm (Dempster et al., 1977) finds a locally optimal estimation in the maximum likelihood sense with provable optimality (Wu, 1983) . Numerous modifications and enhancements to the EM algorithm have been proposed (Jamshidian & Jennrich, 1997; Capp, 2011) . Notably, explicitly modeling the duration of the states in a HMM can be enabled by an extension named Hidden Semi-Markov Models (HSMM) (Ostendorf et al., 1996) . This paper concerns estimating the parameters of a collection of HMMs, each of which corresponds to a set of known features. The observation sequence of each individual HMM is noisy and/or insufficient, making parameter estimation solely based on its corresponding observation sequence using the traditional EM algorithm impractical. The problem arises in recommender systems of online E-Commerce platforms. The objective is to train a HMM (HSMM) model to predict the category of product that the user would purchase based on the user's feature (demographic, location, income level) and observation sequence (past click history). The observation sequence of a particular user is generally noisy (click not reflective of true purchase intention) and insufficient (user only clicked a small number of items), making parameter estimation highly challenging.

The intuition is that while the information regarding each individual user is quite limited, we have access to a much larger amount of information concerning all the users collectively. If we were able to quantify some kind of correlation among the users using their features, we can combine their information for parameter estimation. Formally, we propose to deduce correlations among the HMMs using their features, so that we can fuse the information across different HMMs to estimate the parameters. The basic idea is to combine the classical Expectation-Maximization (EM) algorithm with a neural network which maps a HMM's feature to its parameters. The weight of the neural network captures the correlation among the HMMs and is shared across all the users. The EM iteration, with the estimated HMM parameter at step k (denoted by λ k ) as input and the refined parameter λ k+1 as output, can be considered as a building block (or layer) in the resulted neural network. In order to enable end-to-end training, it is necessary to compute the gradient of the EM iteration ∂λ k+1 /∂λ k . While this gradient exists, we will show in Section 3 that exactly computing this gradient is computationally expensive. A key technical contribution of this paper is in proving the gradient of the EM iteration can be reliably approximated via simultaneous perturbation stochastic approximation (SPSA) (Spall, 1992; Sadegh & Spall, 1998; Spall, 1998; Sadegh, 1997) , which only involves function evaluation and can be easily computed. With such gradient approximation technique, we expect a large body of functions whose gradient is difficult to compute can now be incorporated as layers in deep neural networks, dramatically expanding the choice of network layers and increasing the flexibility of network architecture design. We name our method DENNL (Differentiable EM as Neural Network Layer).

It is worth noting that a user's feature consists of a number of categorical and numerical fields. Instead of hand-crafting a dedicated neural network to encode the user feature, we used a pre-trained low-dimensional user embedding to represent such information. The user embedding was trained on a link prediction task using the graph convolutional network (GraphSage) (Hamilton et al., 2017; Cheng et al., 2016) . Much similar to the success of pre-trained word embedding in downstream natural language processing tasks (Peters et al., 2018; Devlin et al., 2018) , we demonstrate that the embedding trained using GraphSage is effective in summarizing information with complex structures and can be highly instructive in seemingly unrelated downstream tasks.

2. BACKGROUNDS:

Hidden Markov Model and EM Algorithm. In this section, we only introduce the concepts and notations related to HMM and EM algorithm that are necessary for further discussion, and refer the readers to (Rabiner, 1989) and (Bilmes et al., 1998) for a more comprehensive treatment of the topics.

A HMM is defined by a set of hidden states, the transition probability of the hidden states, and the emission distribution. The Markov property mandates that the state transition and the emission distribution is only dependent on the current hidden state. The realizations of the emission distribution are what we can observe as an observation sequence. Let the distribution of the initial hidden state be π and Q t , {t = 1, • • • , T } be the hidden state at time step t. Then transition probability is defined by a time-homogeneous transition matrix A with elements a i,j = P (Q t = j|Q t-1 = i), where i, j ∈ {1, • • • , N } are the indices of the hidden states. Let O t be observable at time step t and o t to be a realization of O t . The emission distribution defines the probability of a particular observation conditioned on the current hidden state

b j (o t ) = P (O t = o t |Q t = j).

(1)

The emission probability of all hidden states can be represented by B = {b j (.)}.

Now we follow the notation in (Bilmes et al., 1998) to present a brief introduction of EM algorithms applied in HMMs. A key concept is the forward and backward probability, which is used to compute the likelihood efficiently. Forward probability is the probability of seeing the partial sequence o 1 , • • • , o t and ending up in hidden state i at time t, Define the forward probability as:

α i (t) = p(O 1 = o 1 , • • • , O t = o t , Q t = i|λ),

where parameters λ = (A, B, π) = {λ 1 , • • • , λ M } and M is the number of parameters. The forward probability can be computed efficiently in an iterative manner,

EQUATION

Similarly, Backward probability is the probability of the ending partial sequence o t+1 , • • • , o T given the hidden state i at time t and is defined as:

EQUATION

Similarly backward probability can also be computed iteratively:

EQUATION

Given the observation sequence o 1 , • • • , o T and the parameter estimation at the k-th iteration λ (k) , the parameter estimation at the k + 1-th iteration can be computed as follows:

EQUATION

EQUATION

EQUATION

where δ ot,v h is an indicator function which equals to 1 if o t = v h and 0 otherwise.

Hidden Semi-Markov Model. HSMM (Ostendorf et al., 1996) explicitly models the duration of a state, by introducing a duration distribution, with which multiple observations could be emitted from a hidden state before state transition occurs. A HSMM can be parameterized as HMM by introducing a counter of the time steps in the current state (Yildirim et al., 2013; Bietti et al., 2015) . For this reason, the computation methods as well as the analysis of HMM naturally apply to HSMM.

3. PROBLEM STATEMENT AND PROPOSED SOLUTION:

In this section we formally state the problem and describe our proposed method. The objective is to estimate the parameters {λ (u) } of a collection of HMM models, using the corresponding observation sequence {o t,1:T (u) } and feature embedding {e (u) } as input. Here (u) is the index of a HMM in the collection of HMMs, and T (u) is the length of the observation sequence. The observation sequence {o t,1:T (u) } is noisy or insufficient for the conventional EM algorithm to apply. The feature embedding e (u) is a low-dimensional vector that encodes our prior knowledge of {λ (u) }. We propose to use neural network to combine the information across all the HMMs for more effective parameter estimation. λ (u) = Ψ e (u) , o t,1:T (u) , Φ (8) where Φ is the parameter of the neural network that is shared by the HMM models. Formally, the loss function to be minimized is

EQUATION

where L(•) is the log-likelihood of observing sequence o t,1:T (u) given HMM parameter λ (u) .

The overview of our proposed method is shown in Figure 1 . Our architecture can be considered as an unrolled Recurrent Neural Network (RNN), where each recurrence (or stage) of the RNN corresponds to an EM iteration. Our key novelty is in rendering EM iteration as a network layer. The resulted "EM layer" serves as the output function of the RNN. This is in contrast to the traditional RNN where the output function is a fully connected layer with non-linearity. Note that the EM layer is parameter-less and directly reflects the domain knowledge and design intent of network designer. As we will show in Section 6, such recurrent structure that mimics the iterative EM iterations in a conventional setting is critical to the success of our proposed method.

As in the conventional EM algorithm, an initial parameter is provided and is iteratively refined. After the kth EM iteration is performed, a multilayer perceptron (MLP) which is labeled as DNN in the figure transforms λ (k) , the latent state from the previous iteration h k , and the user embedding to an updated λ (k+1) and an updated latent state h k+1 . The parameters of the MLP are shared across EM iterations. It is worth emphasizing that the latent state h captures the state of the DNN and shall not be confused with the hidden state of the HMMs. We use the Softplus function and L 1 normalization to ensure the elements of the resulted λ are in proper range.

Differentiability of the EM iteration. The EM iterations shown in Figure 1 can be considered as layers of a neural network, with λ (k) and the observation sequence as input and λ (k+1) as output.

To enable end-to-end training, it is necessary to compute the partial gradient ∂λ (k+1) /∂λ (k) . Even though the EM iteration for HMM can be explicitly stated (Equation 5), directly computing the partial gradient ∂λ (k+1) /∂λ (k) using a automatic differentiation engine faces two numerical challenges, both of which are related to the forward and backward probability (Equation 2 and 4). End-to-end training is enabled by approximating the gradient of the "EM layer" using SPSA, as shown in the low half of the figure.

EQUATION

• While evaluating the alpha-beta equation is a O(T ) operation where T is the length of the observation sequence, the number of nodes required to build the forward and backward probability in the computation graph of an automatic differentiation engine is on the order of O(T 2 ). Empirically we found this leads to intractable computation cost. • The length of the observation sequence T varies across different HMM models, thus a separate computation graph for the forward and backward probability has to be rebuilt for every distinct T value. This is particularly problematic when the HMMs in a mini-batch comes with different T values.

To address the aforementioned numerical challenges, we propose to approximate the gradient ∂λ (k+1) /∂λ (k) using SPSA. Much similar to how we can approximate the gradient of a onedimensional function using the finite difference method, SPSA estimates the gradient of a function by evaluating the function with stochastically perturbed input vectors. A series of work by Spall and colleagues (Spall, 1992; Sadegh & Spall, 1998; Spall, 1998; Sadegh, 1997) established theoretically rigorous error bound for certain scalar functions. Note that such gradient approximation only involves evaluating the EM iteration and can be efficiently computed.

A important question to answer is how well can SPSA approximate the gradient ∂λ (k+1) /∂λ (k) in our specific application. In the next section, we formally describe how SPSA is performed and establish its convergence property for EM iterations.

4. SPSA AND CONVERGENCE ANALYSIS:

As the name suggests, Simultaneous Perturbation Stochastic Approximation approximates the gradient (Jacobian) of a function by simultaneously applying a small perturbation to all the dimensions of the argument of the function. Let ∆ ∈ R M be a M -dimensional random vector, where elements ∆ 1 , • • • , ∆ M are independent with zero mean. We do not assume specific distribution of ∆ as long as E(1/∆) is bounded. We treat the EM iteration stated in equation 5, equation 6 and equation 7 as a function, with the HMM parameter at the k-th iteration λ (k) being the function argument and the refined parameter at (k + 1)-th iteration λ (k+1) as the output:

EQUATION

) where f : R M → R M and o = {o 1 , • • • , o T } is the observation sequence.

The SPSA estimates the gradient (Jacobian) of a function through two function evaluations. The arguments of the two function evaluation are positively and negatively perturbed respectively as following:

EQUATION

where c k is a positive scalar and (k,+) , (k,-) are measurement noise. Let F k be the random field generated by λ (0) , • • • , λ (k) . We assume the measure noise satisfy the following conditions:

EQUATION

)

where 0 is the M -dimensional zero column vector.

With λ (k,+) and λ (k,-) , the approximation of the true gradient (Jacobian) J (k) due to SPSA is defined as

EQUATION

)

where λ -) , ∆ (k) , respectively.

(k,+) i , λ (k,-) i , ∆ (k) i are the i-th element of λ (k,+) , λ (k,

We now quantify the approximation error of Ĵ(k) and establish the convergence property of the gradient approximation due to SPSA for the EM iteration. Our primary theoretical contribution are Lemma 1 and Theorem 1, which are stated below. Note that (Spall, 1992) offered similar theoretical results but their results only apply to functions with scalar output. Our proof took a very different approach from the proof in (Spall, 1992) to bound the Frobenius Norm of the approximated Jacobian to account for the fact that EM iteration is a function with vector input and vector output (f : R M → R M ). Please see Appendix A for the detailed proof.

Lemma 1. For i-th function f i of f in equation 10,

∂ 3 fi ∂λj 1 ∂λj 2 ∂λj 3

exist and are uniformly bounded on the ranges of parameter λ such that for i, j 1 , j 2 , j

EQUATION

)

where o is the observation sequence and τ 0 is a positive scalar.

Theorem 1. Suppose that for each k,

EQUATION

M are independent to each other and have zero mean such that E(∆

EQUATION

where . F is Frobenius norm.

Sketch of Proof::

The difference between the i, j-th element of the approximated Jacobian Ĵ(k) and the true Jacobian J (k) is denoted by

Ĵ(k) i,j -J (k)

i,j . We compute the third order Taylor series expansion (Taylor & Lay, 1958) of this difference, and establish that the first two order terms vanish. The constant term is obviously zero. The first order term disappears because the elements of ∆ k are independent with zero mean. The second order term is also zero because the positive perturbation and the negative perturbation cancel out. According to Lemma 1, the third order term is uniformly bounded, which guarantees the convergence rate of SPSA for estimating the gradient of the EM iteration.

5. RELATED WORK:

It is shown in (Wen et al., 2017) that even a crude approximation of the gradient is sufficient for stochastic gradient descent to converge. Such theoretical result is among our motivations to explore gradient approximation techniques for layers whose gradient is expensive to compute. There are a number of works on enhancing classical sequential data analysis methods with deep neural networks. Krishnan (Krishnan et al., 2017) employs a RNN to parameterize a variational approximation of the posterior distribution in a Gaussian state space model. Other works along this line include (Bui et al., 2016) , (Salinas et al., 2017) and (Linderman et al., 2017) .

Perhaps the most related existing work is (Rangapuram et al., 2018) , where the authors used a deep neural network to parameterize a collection of linear state space models so the information across multiple time series can be shared. The underlying linear state space model is updated via Kalman Filtering, which only involves matrix-to-matrix multiplication and is differentiable. The price to pay for such differentiability is the underlying sequential data analysis method has to be a relative simple form (e.g. linear state space model) with limited expressiveness and flexibility. In contrast, in this work we do not have to limit ourselves to methods whose gradient are easy to compute.

Another closely related work is Gref et al. (Greff et al., 2017) , where the author incorporates the so-called general EM algorithm (Wu, 1983) in the neural network for unsupervised clustering. The key difference is that instead of performing a complete EM update, the general EM algorithm (Greff et al., 2017) only takes one step of gradient descent in the M step of the EM algorithm to refine the parameter estimation. It may take a much larger number of iterations for the general EM algorithm to converge (Wu, 1983) . Larger number of iterations leads to a deeper neural network, as an additional layer is stacked in the neural network per iteration. A network that is too deep may cause difficulty in training (He et al., 2016) . Hinton et al. (Hinton et al., 2018) directly computes the gradient of an EM based routing procedure. With a small number of EM iterations, (Hinton et al., 2018) does not face the numerical challenge outlined in Section 3 because the sizes of their inputs in a mini-batch are fixed and the problem is of much smaller scale.

Relationship with other recommender systems algorithms. Graph embedding (Hamilton et al., 2017; Grbovic & Cheng, 2018) is one of the backbone algorithms for today's industrial recommender systems. While graph embedding is highly effective in capturing complex interactions among hundreds of millions of users and items, most of the algorithms fail to consider the temporal evolvement of the user's interest (with a few exceptions such as (Kumar et al., 2019) , (Trivedi et al., 2017) and (Zhou et al., 2018) , none of which provides an interpretable model of user interest evolvement). We do not consider our proposed method as a general solution for item level recommendation, as it is unlikely that our proposed method in its current form can be scaled to millions of items. Rather, the focus of our proposed method is to capture highly detailed temporal interest evolvement of the users across a limited number of categories.

6. EXPERIMENT:

In this section, we evaluate the effectiveness of our proposed method on synthetic data as well as real-world e-Commerce data.

The primary contribution of our proposed method is to share the information across multiple HMMs for parameter estimation in cases where the observation sequence is noisy or insufficient. To the best of our knowledge, there is no directly comparable method in the existing literature. We perform ablation studies against the following baselines to justify our design choices:

1. Conventional EM. The HSMM parameters are estimated by the conventional EM algorithm, with no information sharing across the HSMM models. 2. EM Mapping. We first estimate the HSMM parameters using the conventional EM algorithm, then train a neural network (a stack of fully connected layers) that maps the user embedding to the HSMM parameters due to conventional EM algorithm. The training objective is to minimize the square error. Though the neural network establishes some correlation across the HSMM models, the HSMM parameters due to the conventional EM algorithm solely relies on the user's corresponding observation sequence that is noisy and/or insufficient. 3. MLP Mapping. In this baseline, a multilayer perceptron (MLP) maps the user embedding to the HMM parameters, then the likelihood is computed based on the HMM parameters and the observation sequence. We use gradient ascent to maximize the likelihood. Note that this baseline is purely based on gradient ascent and does not involve the EM algorithm.

4. Initial Fusion. A DNN maps user embedding to HMM parameters, which are refined by the EM iterations. Here we only fuse the observations across different HMMs at the initial stage. In a sense, this architecture identifies an optimal initial parameter guess. In contrast, our proposed method fuses observations across different HMMs at every EM iteration.

Figure 2 : Overview of baselines. Left: MLP Mapping. Right: Initial Fusion.

6.1. SYNTHETIC DATA:

The purpose of evaluating the proposed method on synthetic data is to compare the decoded hidden state with known ground truth state. Being able to correctly infer the hidden state is a key metric of HMM models (Sarkar & Dunson, 2018) . We pre-designed the parameters of a collection HSMM, each of which has 10 hidden states. The emission distribution is a 5-dimension Gaussian distribution. The duration distribution is Poisson distribution. To generate an informative feature embedding for each HSMM, we serialize the HSMM parameters as a vector and apply a non-linear transformation (in the form of Ax + b followed by ReLU, where x is the parameter vector, A is a known full rank matrix and b is a known vector). The resulted feature embedding is strongly related to its underlying parameter, but one cannot recover the HMM parameters from its feature embedding. The evaluation metric is the percentage of the decoded states (Viterbi, 1967; Munkres, 1957) that match the ground truth state, or the labeling accuracy.

We first evaluate our method in cases where the observation sequence is insufficient. We also compared DENNL to baseline methods when a certain percentage of the observations in the observation sequence generated by the ground truth model is contaminated by noise. If a particular observation is chosen to be noisy, its value is replaced by the value of another randomly chosen observation in the sequence.

The results are reported in Table 1 . 1 , DENNL consistently outperforms the baseline methods with noisy or insufficient observations. When the observation sequence is of sufficient length or noise-less, the benefit of the proposed method diminishes as expected. Even though Direct Mapping method enables some information sharing across HMMs, its performance was poor because the regression target is due to the conventional EM and can be inaccurate to begin with. This experiment illustrates the necessity to fuse information across HMM in every EM iteration. While the Initial Fusion method seems to enable information sharing, its performance is very close to that of conventional EM. This method essentially attempts to identify an optimal initial solution with shared information across HMMs. However, at least in our application, the EM algorithm appears to be insensitive to the initial solution. The MLP Mapper method enables information sharing with simple architecture, but is slow to converge and is not competitive even against conventional EM. To explain this, note that in a conventional setting, we seldom use gradient ascent for HMM parameter estimation. Rather, the common practice is to use the EM algorithm which is carefully designed for the particular application.

This is the key message that we would like to convey to the community: operations that are inspired by well-studied classical methods (such as EM iteration) can be significantly more effective than operations that are synthesized using generic layers (convolution, fully connected, etc). If the gradient of the classical method is incomputable, it could potentially be numerically approximated. This formulation dramatically expands the choices of network layers and allows network designers to directly apply their domain knowledge in the network architecture.

6.2. E-COMMERCE DATA:

As previously discussed in Section 5, DENNL is not considered as a general solution for recommendation systems in its current form. Rather, it is suitable to model the shift of user's interest across a small number of categories. We apply DENNL in a clearly defined and fast-growing sector on our e-Commerce platform, namely Home Decoration. The Home Decoration sector comes with certain well-defined shopping patterns that follow the progress of the decoration of a newly purchased condo. Formally, we model the behavior of the user in the Home Decoration category by HSMM with 17 states, each of which corresponds to a sub-category in the Home Decoration sector. Please see Appendix A.2 for more information of the e-Commerce dataset.

In this setting we no longer have access to the ground truth states, thus we can no longer use labeling accuracy as the metric. Instead, we use a link prediction metric. Using the trained HSMM model, we could compute the likelihood of a user interacting with each of the 17 sub-categories. The prediction is considered to be successful if the user interacted with the predicted sub-category with the highest likelihood, top-3 and top-5 likelihood. Our proposed method outperforms the conventional EM as well as the GraphSage (Hamilton et al., 2017) in this experiment. Note that the user feature we used is the user embedding due to GraphSage, thus DENNL has all the information that GraphSage can offer. One way to interpret the result is that our method is enhancing GraphSage with an HSMM that captures the temporal interest shift of the users.

7. CONCLUSION:

In this work, we proposed a novel method in combining EM algorithm with neural network to predict the parameters of a collection of Hidden Markov Models. The basic idea is to use the neural network to capture correlation among the HMMs so the information across different HMMs can be shared and re-used. The key technical contribution of this work is in providing a proof that the gradient of the EM iteration can be efficiently and reliably approximated using the SPSA. We evaluated our proposed on synthetic data as well as real-world e-Commerce data. An interesting observation is that the user embedding trained in an unrelated task using graph convolutional neural network can be highly informative in HMM parameter estimation. Parallels can be drawn between this observation and the success of pre-trained word embedding in downstream natural language processing tasks.

A APPENDIX:

A.1 PROOF OF LEMMA 1 AND THEOREM 1

Lemma 1. For i-th dimension f i of f in equation 10,

∂ 3 fi ∂λj 1 ∂λj 2 ∂λj 3

exist and are uniformly bounded on the range of parameter λ such that for i, j 1 , j 2 , j

EQUATION

where o is the observation sequence and τ 0 is a positive scalar.

Proof. We first show that that the numerator and denominator of f i is the polynomial function of M -parameters. Since the range of parameters are closed and bounded set, then it follows (Rudin et al., 1964) that the third partial gradient of f i is uniformly continuous and up bounded and then Lemma 1 is established.

Rewrite equation 2 as

EQUATION

where

EQUATION

Then it follows that

EQUATION

where t = 1, • • • , T . For t = 1, let t l=2 Σ (l) = I N where I N is a N × N identity matrix. According to equation 18, we have

EQUATION

where e i is N -column vector, the i-th element equals to one and the rest elements equal to zeros.

Rewrite equation 4 as:

EQUATION

where 1 N is N -dimensional column vector with all elements equaling to one, t = 1, • • • , T . For t = T , let T l=t+1 Σ (l) = I N . According to equation 22, we have

β i (t) = e T i × T l=t+1 Σ (l) × 1 N ,

which together with equation 21 yields that

EQUATION

)

Define N -dimensional column vector as

EQUATION

)

Let m i l ,j l be the i l , j l -th element of Σ (l) . Then the -th element of

EQUATION

Combining equation 23, equation 24 and equation 25, we have

EQUATION

According to equation 19, given that observation o l = h, m i l ,j l = a i l ,j l b j l (o l ) = a i l ,j l b j l (h). From equation 26, it is clearly that α i (t)β i (t) is the polynomial function of parameters {A, B, π}. From equation 5 and equation 26, we have

EQUATION

)

The numerator and denominator of π (k+1) i are polynomial function of parameters λ = {A, B, π}. According to (Rudin et al., 1964) , the third partial gradient of π (k+1) i to λ is uniformly bounded on the range of parameters, which is bounded and closed set on R M . Since a i,j b j (o t+1 ) = m it,jt , where i t = i, j t = j, then similar as equation 26, we have

EQUATION

)

where index i, j, t of coefficients c(i,j,t, ) i2,j2,••• ,i T ,j T corresponds to i, j, t in α i (t)a i,j b j (o t+1 )β j (t + 1). Combining equation 6, equation 26 and equation 28, we have

EQUATION

Similarly, it follows from equation 7 and equation 26 that

EQUATION

Note that the numerator and denominator of a 

EQUATION

M are independent to each other and have zero mean such that E(∆

EQUATION

where . F is Frobenius norm.

Proof. According to equation 11, equation 12, and equation 14, the i-th row, j-th column element difference between Ĵ(k) and J

EQUATION

Taylor series expansion (Taylor & Lay, 1958) of f i (.) is represented as

EQUATION

f i λ k -c k ∆ (k) = f i (λ (k) ) - M j1=1 c k ∆ (k) j1 ∂f i (λ (k) ) ∂λ j1 + 1 2 M j1=1 M j2=1 c 2 k ∆ (k) j1 ∆ (k) j2 ∂ 2 f i (λ (k) ) ∂λ j1 ∂λ j2 - 1 6 M j1=1 M j2=1 M j3=1 c 3 k ∆ (k) j1 ∆ (k) j2 ∆ (k) j3 ∂ 3 f i λ (k) -ξ (-) ∆ (k) ∂λ 1 ∂λ 2 ∂λ 3 ,

where ξ (+) , ξ (-) ∈ (0, c k ). Then it follows that

EQUATION

where 

EQUATION

EQUATION

In addition, from equation 13, we have

E (k,+) i - (k,-) i 2c k ∆ (k) j F k = 0

which together with equation 32 and equation 36 yields that as k → ∞, almost surely

E Ĵ(k) i,j -J (k) i,j F k = O(c 2 k ), It follows that E Ĵ(k) -J (k) F k F = O(c 2 k ).

A.2 DESCRIPTION OF E-COMMERCE DATASET:

The operation specialist who has deep domain knowledge in the Home Decoration sector have concluded that after a condo is purchased, the customer will likely go through 4 distinct phases: Planning , Indoor Construction, Furniture, and Accessories. In the Planning phase the customer mainly browses listing of designing services and contractors. Representative items in the Indoor Construction includes electrical and plumbing supplies. The user's interest changes according to the actual decoration progress of the condo. Once we have identified that a user is likely to switch from one phase to another, we can begin to recommend the items that belong to the next phase. The ability to incorporate such detailed temporal information is a key improvement over existing recommender system algorithms.

The hidden states include the four phases discussed above, as well as a NotDecorating state to indicate that the user is not actively considering items in the Home Decoration. We summarize the user's behavior on the platform on a daily basis. Each day, we assign a strength value to a total of 17 sub-categories (each of the four phases consists of four sub-categories, plus one for no action). The strength value is calculated based on the number of clicks/purchases and serves an indicator of the level of interest of the user in a particular sub-category. The strength value across the 17 sub-categories is normalized to 1. With this formulation, the emission probability is the Multinomial distribution. The duration distribution is Poisson distribution.

The user feature that we provided to the neural network is a 128 dimension embedding trained in a link prediction task using GraphSage. The task is to predict whether a user would click an item solely based on the user embedding and item embedding. Such embedding not only encodes the numerical and categorical information of the users (such as location and income level), but also captures the user-to-item interaction on the platform. We chose 10000 users on the platform that has a large number of clicks/purchases in Home Decoration, and used 180 days of their behavior between January 2018 and July 2018 as training data. The user's behavior between July 1st and July 2nd 2018 is used as validation data.

:



Figures and Tables:

Figure 1: Overview of DENNL. The EM iteration are incorporated in a deep neural network as layers. End-to-end training is enabled by approximating the gradient of the "EM layer" using SPSA, as shown in the low half of the figure.

are also polynomial function of parameters λ = {A, B, π}. Then the third partial gradient a are also bounded on the range of parameters. Since f in equation 10 is equivalent to equation 27, equation 29 and equation 30, the first partial gradient of f is uniformly bounded.Theorem 1. Suppose that for each k, ∆

< τ 2 , as k → ∞, then given F k , the expectation of RES in equation 33 is bounded such that as k → ∞, almost surely |E(RES|F k )| ≤ τ

Labeling accuracy comparison of DENNL to the baseline methods on the synthetic data.

Table content: <table><tr><td/><td>Sequence</td><td>Sequence</td><td>Sequence</td><td>Noisy</td><td>Noisy</td></tr><tr><td>Method</td><td>Length:</td><td>Length:</td><td>Length:</td><td>Sample:</td><td>Sample:</td></tr><tr><td/><td>50</td><td>100</td><td>400</td><td>20%</td><td>10%</td></tr><tr><td>Conventional EM</td><td>82.0%</td><td>85.5%</td><td>94.1%</td><td>83.1%</td><td>91.5%</td></tr><tr><td>Direct Mapping</td><td>71.4%</td><td>72.2%</td><td>73.0%</td><td>67.9%</td><td>71.4%</td></tr><tr><td>MLP Mapper</td><td>51.5%</td><td>50.3%</td><td>52.1%</td><td>46.2%</td><td>48.5%</td></tr><tr><td>Initial Fusion</td><td>82.4%</td><td>83.7%</td><td>92.6%</td><td>82.8%</td><td>89.7%</td></tr><tr><td>DENNL</td><td>88.3%</td><td>90.0%</td><td>93.5%</td><td colspan="2">88.1% 92.3 %</td></tr><tr><td>Discussion: As shown in Table</td><td/><td/><td/><td/><td/></tr></table>

Sub-cateogry prediction accuracy in the e-Commerce dataset. Method

Table content: <table><tr><td/><td>Top-1</td><td>Top-3</td><td>Top-5</td></tr><tr><td>Conventional EM</td><td colspan="3">10.3% 23.9% 50.5%</td></tr><tr><td colspan="4">GraphSage (Hamilton et al., 2017) 17.7% 37.4% 65.5%</td></tr><tr><td>MLP Mapper</td><td>7.4%</td><td colspan="2">19.6% 35.1%</td></tr><tr><td>Initial Fusion</td><td colspan="3">11.5% 21.6% 48.2%</td></tr><tr><td>DENNL</td><td colspan="3">21.8% 48.5% 75.2%</td></tr></table>

Peer review (iclr2020):
The authors propose to use numerical differentiation (using random perturbation) to approximate the Jacobian of a particular update (essentially equations 5~7) which plays an important role in the estimation of HMMs. To do so, the authors provide first a concise intro to HMM models (well known stuff in S2), presenting the iteration in detail, jump into their model (cryptically presented in my opinion in S3) and then propose a numerical approximation scheme using SPSA (building upon literature from the 90*s, with Theo 1 being the main contribution), before moving onto experiments. I have found the paper poorly presented. Its general motivation stands on a shaky ground (as illustrated by the choice of words by the authors, see below). In terms of presentation, reminders on HMM are welcome, but unfortunately the authors have not kept the same standard for clarity of notation in Section 3, which makes reading and understanding what the authors are doing quite difficult. Not being a specialist in this field, I have struggled a bit to understand the model itself, and the practical motivation of adding a DNN in the middle of what is otherwise an unrolled back-and-forth between k steps of EM estimation of transition parameters and the addition of a DNN layer. Despite the complexity in the story, what the authors propose is essentially to apply a numerical approximation scheme for Jacobians of these EM updates instead of backprop. Since this is the crux of the contribution, I feel some more numerical evidence that their approach works compared to baselines (e.g. Hinton et al 2018) is needed. For these reasons my assessment is a bit on the lower side. - parenthesis bug in b_j(... in Eq.4 - in Eq. 5, index i appears both in numerator (as regular index) and denominator (as sum index) - what is Psi in Eq.8 ? - *While HMM is arguably less prevalent in the era of deep learning*: odd way to start an intro. All papers cited date back to more than 2011, 2 in 2006, all the rest in 20th century. This is particularly strange given the few citations to papers >2015 in Section 5. - the observation sequence o_{t,1:T(u)} is *weakly* indexed by u (since T(u) is just a length) - What is the forall u notation below Eq. 9? - *the number of nodes required to build the forward and backward probability in the computation graph of an automatic differentiation engine is on the order of O(T^2). Empirically we found this leads to intractable computation cost.* since this is critical, where is this empirical evidence? this seems to be a storage problem and cannot be a complexity issue. There are ways to mitigate this problem by only storing partially information, I feel this comparison would add a lot of value to the authors* claim. - Where is J^{(k)} defined (as opposed to hat{J}^{(k)}) defined in Eq.14?