Earth Mover's Distance Pooling over Siamese LSTMs for Automatic Short Answer Grading

1 Introduction

Grading student work is critical for assessing their understanding, and providing teachers instructive feedback. However, answer grading can become monotonous and tedious for teachers. Computer assisted assessment has been in practice in schools and colleges for many years now, but primarily for recognition questions with constrained answers such as multiple choice questions. Prior research has shown that such recognition questions are deficient in that they do not capture multiple aspects of acquired knowledge, such as reasoning and self-explanation [Conole and Warburton, 2005]. Consequently, open-ended recall questions that seek responses constructed by students are more commonly used in academia. The focus of this paper is automatic grading of such constructed answers with reference to instructor-provided model answers. In particular, we are interested in short answers between a few words and a few sentences long, i.e., everything in between fill-in-the-gap and essay-type answers [Burrows et al., 2015; Roy et al., 2015]. Table 1 shows an example of a short answer grading task.

Grading student-constructed short answers, given instructor-provided model answers, is a complex natural language understanding task owing to linguistic variations (the same answer could be articulated in different ways), the subjective nature of grading (multiple possible correct answers or no correct answer) and lack of consistency in

Shourya Roy

American Express Big Data Labs shourya.roy@aexp.com

human rating (non-binary scoring on an ordinal scale within a range). Despite state of the art results in most natural language processing tasks, neural models have not been applied extensively for ASAG, barring the use of word embeddings for supervised ASAG [Sakaguchi et al., 2015] and neural networks for essay grading [Alikaniotis et al., 2016]. In this paper, we present a new ASAG system comprising a novel pooling method over Siamese LSTMs, followed by a flexible regression layer for generating a score. Specifically, we make the following contributions.

1.1 EMD Pooling over LSTMStates

Although natural language inference (NLI) has some critical differences from ASAG (see Section 2), we first tried two network architectures established in NLI. In the Siamese architecture, two networks are applied to two texts, and their final states combined to make predictions [Severyn and Moschitti, 2015]. This has been superseded by a more global comparison across all pairs of states, using various attention mechanisms [Rocktäschel et al., 2015; Yin et al., 2015; He and Lin, 2016]. Surprisingly, attention mechanisms did not perform well for ASAG (see Section 6.4).

Instead, we solve an Earth Mover Distance (EMD) problem on a matrix of pairwise distances between each state vector of the model and student answers. EMD is a metric and more generally known as Wasserstein distance [Levina and Bickel, 2001]. EMD has been used between distributional word vectors [Kusner et al., 2015], but not to recurrent states, to our knowledge. Part of the technical challenge is to backprop from ASAG-appropriate losses through the EMD pooling to the LSTMs, in the face of certain aggregate constraints in EMD. This pooling layer may be of independent interest beyond the ASAG task.

1.2 Sinkhorn-KnoppMatrix Scaling

To alleviate the cubic time complexity of EMD calculation, we replace the linear program EMD solver with the SinkhornKnopp matrix scaling [Sinkhorn and Knopp, 1967; Cuturi,

2013] procedure, which replaces EMD with the Sinkhorn distance [Knight, 2008], an excellent approximation. Another potential benefit is that matrix scaling is a differentiable operator [Huang et al., 2016].

1.3 Ordinal Loss Optimization

Another important difference between ASAG and NLI is that the desired output in ASAG is a real or ordinal score, not a categorical label (entailment, contradiction, irrelevance) as in NLI. The primary goal in ASAG is fair grading, i.e., to minimize numeric difference between gold and system scores. Large correlation usually follows, but is of secondary importance. Therefore we must append a final ordinal prediction layer to the EMD layer. Whereas softmax is most popular for categorical prediction, we found that a two-way hinge loss used in support-vector ordinal regression (SVOR) [Chu and Keerthi, 2007] is superior to a natural non-convex choice [Cheng, 2007].

1.4 Data Augmentation

Training neural models often require large amounts of training data, which is rare among public ASAG datasets. Our primary dataset from [Mohler et al., 2011] contains only about 2,200 (model answer, student answer) pairs. In computer vision, images are readily augmented by scaling, rotation, and changes in brightness and contrast. Here we use a simple but effective method to increase the number of training pairs.

Organization: After reviewing related work in section 2, we present the components of our system in section 3. Details of the earth mover pooling layer are in section 4. Section 5 describes how the EMD pooling output is turned into a score. Section 6 report on experiments on public ASAG benchmarks, and section 7 concludes the paper. We intend to place our code in the public domain.

2 RelatedWork

In this section, we review prior work in two relevant research streams viz. supervised ASAG and neural models for related NLP tasks with emphasis on NLI.

SupervisedASAG: Early supervisedASAGsystems incorporated task- and dataset-specific features into standard classification and regression algorithms [Pulman and Sukkarieh, 2005]. Ensembles of classifiers using different subsets of features [Heilman and Madnani, 2013] are widely used. Features are extensively tuned towards specific datasets, and new techniques are rarely evaluated on datasets published earlier. A comprehensive overview of prior supervised ASAG systems can be found in two recent survey papers and references therein [Burrows et al., 2015; Roy et al., 2015]. These surveys have also emphasized the importance of standardized evaluation across datasets. [Ramachandran et al., 2015] proposed a graph based approach to extract patterns from groups of questions and their answers towards constructing regular expression alike patterns and showed better performance than the Tandalla's winning entry in [Kaggle, 2012b].

Recent ASAG efforts have leveraged advancements in distributional word semantics and neural methods for embedding words. Word level similarities between student and model answers based on pre-trained word vectors [Mikolov et al., 2013b] have been used as features to train regression and classification models [Sakaguchi et al., 2015; Sultan et al., 2016]. [Alikaniotis et al., 2016] introduced score-specific word embeddings in LSTMs to embed student and model essays to single state vectors. These were turned into an output score via linear regression with square loss.

Natural Language Inference (NLI): In the NLI task [MacCartney, 2009], we are given a premise sentence (e.g., The cat is playing Frisbee on the beach) and have to judge if a hypothesis (e.g., A pet is using a plastic toy on sand; or An elephant is bathing in the river) is entailed by the premise, or contradicts it, or is unrelated. Thus, a special case is recognizing textual entailment (RTE) [Sammons et al., 2011]. A public NLI corpus was recently released [Bowman et al., 2015]. Their baseline approach represented the premise and hypothesis as concatenated word embedding vectors, which were input to a multi-layer LSTM-RNN architecture that outputs one of the above classes.

Various Siamese networks with tied weights have been used to compare or label pairs of short texts. [Severyn and Moschitti, 2015] used Siamese convnets to match candidate answer passages to queries. [Mueller and Thyagarajan, 2016] used Siamese LSTMs for NLI. In recent proposals [Rocktäschel et al., 2015; Wang and Jiang, 2015; Liu et al., 2016], each hypothesis token attended to specific focus tokens in the premise sentence. Some NLI approaches [Yin et al., 2015; He and Lin, 2016] compute attention from a matrix of pairwise interactions like us, but these are then used to predict categorical NLI labels.

Critical differences between NLI/RTE and ASAG: ASAG has often been connected to NLI/RTE by regarding the student answer as premise and the model answer as hypothesis [Ostermann et al., 2015; Mohler et al., 2011]. However the differences between the tasks are stronger than their apparent commonality. [Ostermann et al., 2015] conducted an entailment annotation exercise with fine grained labels and compared with scores assigned by teachers. They found:

• Not all RTE tags could be mapped to ASAG scores. ASAG scores are typically ordinal, e.g., on a 5-point scale whereas RTE tasks have categorical labels such as 'entailed', 'not entailed', and 'contradictory'
• Partial entailment instances had poor agreement or correspondence with human-provided ASAG scores.

In ASAG, key concepts in the model answer must be covered by the student answer, but the concepts may be diffused over many tokens, and extra material in the answer is not necessarily a disqualification. A student answer which is more specific than the model answer may not entail the latter, but nonetheless is likely to be scored as correct. Finally, the organizers of the "student response analysis" (SRA) task in SemEval-2013 workshop concluded that the correlation between answer assessment judgments and entailment judgments is not perfect. 1

Finally, there are two other NLP tasks which are related to ASAG. In automatic essay scoring (AES), the objective is to score essays on the basis of composition, fluency, grammatical correctness etc. Unlike ASAG, AES evaluates essays as per "absolute" criteria, not relative to instructor-provided model answers [Kaggle, 2012a]. In paraphrase detection (PD), the objective is to detect if two sentences or passages have the same meaning [Socher et al., 2011] can also be seen as similar to ASAG. However, unlike PD, ASAG is asymmetric and there may not be a one-to-one correspondence between concepts in model and student answers.

1 https://www.cs.york.ac.uk/semeval-2013/task7/

3 SystemOverview

Figure 2 gives a high-level view of our system. The components are described in detail in sections 4 and 5. During inference, the system gets a question, associated with a model answer M and a student answer A, and has to estimate a nonnegative real or ordinal score yA. During training, gold yA are provided. M and A are modeled as sequences of tokens. Tokens are mapped to pre-trained (300-dimensional) word embeddings [Mikolov et al., 2013b] and input to a conventional Siamese bi-LSTM [Hochreiter and Schmidhuber, 1997] with tied parameters, which we call LSTMM and LSTMA.

Suppose there are ℓM tokens in model answer M and ℓA tokens in student answer A. Bi-LSTM LSTMM emits hM = (hM,i : i ∈ [1, ℓM]) where hM,i = [-→h M,i, ←- h M,i] is a concatenation of the forward and backward states at each token position i ∈ [1, ℓM]. Similarly, LSTMA emits hA = (hA,j = [-→h A,j , ←- h A,j]) at each token position j ∈ [1, ℓA]. In our implementation each h vector is in R50+50 . As Figure 2 shows, hM and hA are input into a pooling layer described in section 4. The output of the pooling layer goes to a flexible regression layer, discussed in section 5.

4 EarthMover Distance (EMD) Pooling Layer

At this stage, many Siamese LSTMs would compare/combine state vectors -→h M,ℓM and -→h A,ℓA (and possibly ←- h M,1 and ←- h A,1) to predict a label. But, building on our intuition about ASAG thus far, we inject all of hM and hA into a "crossbar" type pooling layer, defined by the earth mover distance (EMD). We first define the distance between two hidden states: dij = ‖hM,i - hA,j‖2 .

Next, we define earth mover distance as a minimization over an auxiliary transport matrix T ∈ RℓM×ℓA :

where ∑i Tij = 1/ℓA and ∑j Tij = 1/ℓM. (2) Intuitively, EMD is defined by a minimal fractional transportation between the state vectors in hM and hA. [Kusner et al., 2015] calculated EMD between word embedding sequences, but not to LSTM states that capture long distance semantics in the answers. In fact, in section 6, we show that EMD between words is inferior to that between LSTM states.

4.1 Alternating Optimization

We need to backprop errors through the EMD pooling network to the LSTM weights Θ, because hM,hA, and, by extension, dij change as training progresses. However, in (1) and (2), T 's dependence on Θ is not in a closed form. We address this problem using alternating optimization. Every half-step, we fix Θ, their hidden states hM,hA, and therefore dijs. Then we optimize T using LP. In the other half-step, we fix T and backprop the loss into Θ. For reasonable choices of the regression and loss layers in Figure 2 (see section 5.3), the loss becomes differentiable wrt Θ.

4.2 Sinkhorn Approximation

Unfortunately, the time complexity of solving the EMD optimization via LP is O(d3 log d) where d = max{ℓM, ℓA} [Kusner et al., 2015]. This becomes a bottleneck for training the network, because this computation needs to be performed for every single update.

To alleviate the cubic time complexity of EMD, we follow the approach of [Cuturi, 2013] to add an entropy regularization term to the transport objective (1). For any fixed λ > 0, the regularized transportation problem is defined as

subject to the constraints (2). The larger the λ, the closer this relaxation is to the original EMD. This modified optimization is strictly convex, and an efficient matrix scaling algorithm [Sinkhorn and Knopp, 1967; Cuturi, 2013] can be used to solve it in O(d2 ) time, giving the Sinkhorn distance, a good approximation to EMD. It may be possible to backprop through the matrix scaling steps to the LSTMs as well [Huang et al., 2016], but this is left for future work.

5 Output Regression Layer and Loss Design

The EMD pooling layer outputs a real scalar, but this does not directly correspond to answer grades/score. In this section we discuss alternatives to obtain a score from the LSTM or EMD outputs. The best choice may depend on how the ASAG system will be evaluated.

5.1 Regression to Continuous Score (Baseline)

A standard Siamese LSTM system would compare -→h M,ℓM with -→h A,ℓA . These can be combined into a realvalued answer score predictor. E.g., [Mueller and Thyagarajan, 2016] use the following transformation: exp (- ∥∥∥-→h M,ℓM - -→h A,ℓA ∥∥∥1) which is large if the LSTM states are similar, and small if they are different. The final loss can be L1 or L2 (mean square) error between the above prediction and ground truth. (More about the choice of losses in section 6.3.)

5.2 Ordinal Labels with Non-convex Loss (Logits)

In a standard multi-class classification problem with categorical labels 1, . . . ,K (no order between them), neural networks output a multinomial distribution (pk : k ∈ [1,K]) from a last softmax stage, whose cross-entropy KL(p||· · ·) from a one-hot gold distribution (. . . , 0, 1, 0, . . .) is used as the loss. In case of ordinal regression [Cheng, 2007; Gutiérrez et al., 2016], the gold label representation changes from one-hot to "prefix-hot": if the label is k, the first k elements of the vector are ones, and the suffix all zeros, i.e., (1, . . . , 1, 0, . . . , 0). The goal of the regression stage is to learn a function that maps the input (EMD, say) to a vector o = (o1, o2, ..., ok, ...oK), where oi ≈ 1 if i ≤ k and oi ≈ 0 if i ≥ k. ‖o‖1 is an estimate of k itself, rather than 1.

Suppose the pooled output in Figure 2 is z ∈ R. For EMD pooling, larger z implies smaller k and vice versa, so we will negate z in the following. We include in the trainable parameters b1 ≤ b2 ≤ · · · ≤ bK, and predict output vector pk = σ(-z - bk), where σ(·) = 1/(1 + e-· ) is the sigmoid function. As loss between p and ground truth, we can use L1 or L2 or KL divergence. During inference, to make a prediction, our method scans p1, p2, . . . , pK in order. It stops when pk drops below some predefined threshold, or k reaches K.

5.3 Ordinal Labels with Convex Loss (SVOR)

Instead of using b1, . . . , bK as the transition points for sigmoids, we can adopt standard convex formulations for ordinal regression [Herbrich et al., 1999; Chu and Keerthi, 2007] and use them to define consecutive ranges within which the EMD output must fall to emit successive ordinal scores.

To ease notation, let -∞ = b0 ≤ b1 ≤ · · · ≤ bK = +∞, where only b1, . . . , bK-1 are learnable parameters, inducing K intervals in R. If the gold label is k, then -z should lie in [bk-1, bk]. Tacking on a margin, we define the "tub loss" as a two-sided hinge loss around [bk-1 + 1, bk - 1]:

The alternating training protocol of section 4.1 can now be formally specified as the LSTM-EMD-SVOR algorithm.

• 1: Initialize tied LSTM parameters Θ and regression parameters b
• 2: while validation error reduces do
• 3: for each labeled instance 〈M, A, k〉 do
• 4: Input M, A to LSTMs and calculate dij s
• 5: Fix {dij}, find optimal T via matrix scaling
• 6: Fix T , calculate z, then tub(z, k; b) loss
• 7: Backprop to update Θ and b

During inference, Θ, b are fixed. We input M, A as before, compute z, and then report k such that -z ∈ [bk-1, bk]. Depending on the optimizer, we can enforce bk-1 ≤ bk directly, or by a suitable change of variable, like bk-1 + eck = bk.

6 Experiments

We evaluate our system on two grading tasks that we call CS Dataset and SemEval Dataset, described in section 6.1. In section 6.2 we describe how we augmented them to better train our networks. Subsequently we describe different evaluation metrics used, in section 6.3. The rest of this section presents accuracy numbers comparing with relevant prior art.

6.1 Datasets

We use two datasets for reporting experimental results which we describe next.

CS dataset: This is one of the earliest ASAG datasets consisting of 80 questions from 12 assignments of an undergraduate "Data Structure" course leading to total 2, 273 student responses [Mohler et al., 2011]. 2 Each student answer was evaluated on a scale of [0, 5] with 0.5 interval leading to 11 permissible scores. Student answers were independently evaluated by two annotators, with their average designated as the gold label. We train the proposed network with the ordinal-loss based final layer on this dataset for each assignment by training on the remaining 11 assignments. Before data augmentation (section 6.2), the average size of training set (in 12-fold cross validation) is about 2, 100 student responses.

SemEval dataset: This dataset is a part of the "Student Response Analysis" (SRA) in the Semantic Evaluation (SemEval) workshop in 2013 [Dzikovska et al., 2013]. 3 We use the SCIENTSBANK subset of the data which contains approximately 10, 000 answers to 197 assessment questions from 15 different science domains. The answers were graded by multiple annotators on a nominal scale viz. 'Correct', 'Partially correct/incomplete', 'Contradictory' (student answer contradicts the reference answer), 'Irrelevant' and 'non domain'.

Owing to its ordinally incomparable labels, we cannot report numbers that are directly comparable to published numbers on this task. Nevertheless, we adapt it to compare between the alternative network architectures we have proposed. We use the ordinal labels 'Correct', 'Partially correct/incomplete', and 'Irrelevant' and assign them ordinal labels 2, 1 and 0, to be able to report MAE, RMSE and correlation as in the Mohler data set.

SemEval test protocols: The test set is divided into three subsets with varying degrees of similarity with the training examples. The Unseen Answers (UA) dataset consists of responses to questions that are present in the training set. Unseen Questions (UQ) contains responses to in-domain but previously unseen questions. Three of the fifteen domains were held out for a final Unseen Domains (UD) test set, containing completely out-of-domain question-response pairs.

2 http://web.eecs.umich.edu/~mihalcea/downloads/ ShortAnswerGrading_v2.0.zip

3 https://www.cs.york.ac.uk/semeval-2013/task7/index.php\ %3Fid=data.html

6.2 Training Data Augmentation

Public ASAG datasets are rarely large enough to effectively train a large number of model parameters in complex neural networks. Vision and speech [Cui et al., 2015] labeling tasks augment training data using label-preserving transformations. [Zhang and LeCun, 2015] used a thesaurus derived fromWordnet to replace some words from text with their synonyms and thereby increasing training data size. Here we describe our data augmentation method. A substantial number of the answers in our data-set are given perfect scores as per the instructors. We hypothesize that the student answers which received perfect scores by the instructor are equivalent to instructor-provided model answers. For example, if out of n students m received perfect scores (m < n) then by above strategy we can generate m × (n - 1) new (model answer, student answer) training pairs. Augmenting the Mohler data set results in an increase from about 2500 to 35000 training pairs. Augmenting the SemEval data set increases training pairs from about 5000 to 78000. During testing, we use only the model answer provided by the instructor.

6.3 EvaluationMetrics

ASAG systems that predict ordinal labels have been evaluated using a variety of metrics, in two broad categories: value deviation and correlation.

Mean absolute error (MAE): If y∗ is the gold score and y is the system-generated score, this is simply |y-y∗ |, averaged over instances.

Root mean square error (RMSE): If i ∈ [1, I] indexes instances, this is √(1/I) ∑i (yi - y∗ i )2 .

Correlations: Here the goal is to compare the orders imposed by gold and system scores over all answers. Pearson's r is the most popular correlation coefficient in ASAG.

While both absolute and correlation measures have been used in ASAG, we argue that absolute deviation makes more sense, because we are interested in fair assessment, rather than merely ordering students. Further, MAE is believed to be superior to RMSE [Willmott and Matsuura, 2005], in part because of nonuniform scaling: errors < 1 shrink, while errors > 1 expand when squared. Nevertheless, we report all the measures for easy comparison with relevant prior art.

The upper block of Table 3 gives baseline numbers. Tfidf is the simple TFIDF similarity between M and A. Lesk is a standard WordNet-based similarity (see [Mohler et al., 2011]). The numbers for [Mohler et al., 2011] are from that paper. We ran code provided by [Yin et al., 2015] and [He and Lin, 2016] on our ASAG data. We also tried the attention model of [Rocktäschel et al., 2015], but failed to get any better baseline results than the ones shown. The second block shows ablation studies on our network architecture. The lower block shows the best design choices.

6.4 Ablation Studies

Although exhaustively covering all possible combinations in {Word2vec, ConvNet, LSTM} × {last state, all states} ×

{Avg, Max, EMD} × {SoftMax, L2, Logits, SVOR} would be too tedious, we compare some important points in this design space.

Choice of Output Stage: In table 3, using L2 loss at the output generally results in poorer performance. Results are mixed among the best two choices in rows 13 and 14. SVOR is much better than Logits for MAE and RMSE, but not as good for correlation, which may be because our loss function is not ranking-oriented. Either may be used, depending on whether the assessment goal is accurate scoring or ranking.

Last-state vs. All states: LSTM-Last-L2 is the only row where the last state vector was used, and is generally inferior to most methods where all states were used (except with L2 output loss). To better understand the value of information from all states, we plotted Tij from sample M,A pairs as heatmaps in Figure 4, for one good and one poor answer. The good heatmap has large bright patches representing high transportation values (i.e., LSTM state matches), where as the poor heatmap has large dark patches. Last-state regression score roughly corresponds to the brightness of the bottom right corner cell TℓM,ℓA . The heatmaps show that even a good (respectively, poor) answer may have a relatively low (respectively, high) value of this cell. Therefore, incorporating signals from all intermediate states is important. Attentionbased networks also try to capture this, but they have many more parameters which are more delicate to optimize compared to our frugal model (rows 4 and 5 in Table 3).

EMD over words vs. LSTM states: At this point one may be convinced that a whole sequence-to-sequence match is important, but it would be tempting to try a simpler alternative to LSTM, namely, the word mover distance of [Kusner et al., 2015] which computes EMD across stateless sequences of word embeddings [Mikolov et al., 2013a]. But, comparing row 10 against the most closely comparable rows 13 and 14 of Table 3, we see clear evidence that EMD pooling over stateless word embeddings is not nearly as good as EMD pooling over LSTM states.

6.5 Comparison with Current Systems

Comparing rows 3 and 13 in Table 3 shows that we perform better than [Mohler et al., 2011]. They employ a support vector regression machine that predicts scores using a set of dependency graph alignment and lexical similarity measures while our method uses no handcrafted features.

[Ramachandran et al., 2015] adopt a different evaluation setup. For each assignment/test, they use 80% of the data for training and the rest as test. This setup thus enables in-domain model training. Their system trains a random forest regressor based on features derived from automatically generated regexp patterns to capture semantic variations and syntactic structures of good answers. Results in this setup are shown in Table 6. Our model performs better on RMSE and matches on correlation (they did not report MAE).

6.6 SemEval Performance

Table 5 compares some of the system variations in Table 4.1 for the SemEval data. This data set has many tied gold and systems scores, which explains the low correlation across the board. However, for MAE and RMSE, the trends seen in the Mohler CS data are more-or-less preserved here, with LSTMEMD-SVOR emerging superior to other combinations.

6.7 Effect of λ

Table 7 shows how MAE and time (for one EMD or Sinkhorn calculation) varies with λ. We chose λ = 10 for all experiments as a robust compromise.

7 Conclusion

In this paper, we proposed a novel neural architecture for automatic short answer grading (ASAG). Our system combines Siamese bi-LSTMs, a novel pooling layer based on the Sinkhorn distance between LSTM state sequences, and a support vector ordinal output layer. Training is enhanced via a task-specific data augmentation strategy. Experiments on two publicly available ASAG datasets established that our system has scoring accuracy superior to recent baselines.