Beyond Sentential Semantic Parsing: Tackling the Math SAT with a Cascade of Tree Transducers

We present an approach for answering questions that span multiple sentences and exhibit sophisticated cross-sentence anaphoric phenomena, evaluating on a rich source of such questions – the math portion of the Scholastic Aptitude Test (SAT). By using a tree transducer cascade as its basic architecture, our system propagates uncertainty from multiple sources (e.g. coreference resolution or verb interpretation) until it can be confidently resolved. Experiments show the first-ever results 43% recall and 91% precision) on SAT algebra word problems. We also apply our system to the public Dolphin algebra question set, and improve the state-of-the-art F1-score from 73.9% to 77.0%.

Introduction
Math word problems pose questions that are challenging for current question answering (QA) systems to solve. Consider the following question originating from a study guide for the Math SAT 1 : Example 1: Suppose 3x + y = 15, where x is a positive integer. What is the difference between the largest possible value of y and the smallest possible value of x, assuming that y is also a positive integer?
The correct response is 11; however its relationship with the other numbers in the question (3 and 15) is oblique and not easily mapped to an operator tree or equation template. This encourages us to build a semantic parser that produces an explicit representation of what the question is asking, if we want to make quantitative progress on the question set. However, while it is not hard to formalize the semantics: it is not clear how to devise a compositional transformation from the original question to the formal semantics, since the meaning is dispersed throughout the discourse, such that neither the maximization nor the minimization can be locally derived from some subtree of the syntactic structure.
Moreover, SAT questions quickly reach the limits of preprocessing tools like anaphora resolution: Example 2: r, s, t In the sequence above, if each term after the first is x more than the previous term, what is the average of r, s, and t in terms of r and x?
Understanding this question requires a nuanced resolution of each term after the first to the subsequence s, t , a coreference resolution beyond the grasp of the current state of the art.
Generally speaking, question discourse (with its complex cross-sentence semantics and anaphora) has not been a major focus of QA research. In this paper, we use math SAT questions to develop an approach to handling question discourse. Our parser uses an intermediate semantic language that allows complex semantics (like those of Example 1) to be compositionally constructed from a multisentence question passage (Section 5.1). By architecting our semantic parser as a cascade of nondeterministic tree transducers (Gécseg and Steinby, 1997), we can propagate uncertainty until it can be confidently resolved -sometimes as late as during program interpretation (Section 6). The integrated approach also allows us to handle novel classes of anaphoric phenomena by framing anaphora resolution as an operation on a parse forest decorated with implicits (Section 5.2).
Ultimately we produce an end-to-end system (called EUCLID) that achieves 43% recall and 91% precision on SAT closed-vocabulary algebra questions, a subset (described in more detail in the next section) that constitutes approximately 45% of a typical math SAT. We also achieve state-of-the-art results on the publicly released Dolphin question set (Shi et al., 2015), a set of more than 1500 algebra questions released by Microsoft Research. Finally, we provide a look at our early progress on extending the system to tackle the math SAT in its entirety.

Anatomy of a Math SAT
To assess our semantic parser, we compiled three question sets. Two question sets were created from sample SAT exams found in study guides (published by Kaplan and McGraw-Hill). We used the Kaplan set (12 exams, 648 total questions) for training/development and the McGraw-Hill set (13 exams 2 , 686 total questions) for devtest. We reserved official practice exams (8 exams, 396 total questions) released by the College Board for final testing. We did not subselect questions from the exams, rather we used them in their entirety. 3 We encoded mathematical formatting using LATEX.
During the compilation of these questions, they were split into 4 broad categories: 1. Algebra (closed vocabulary) (e.g. Examples 1 and 2) : Algebra questions drawn from a limited mathematical vocabulary.

2.
Algebra (open vocabulary) (e.g. "At a basketball tournament involving 8 teams, each team played 4 games with each of the other teams. How many games were played at this tournament?") : Algebra questions drawn from an open-ended vocabulary.

4.
Other A catch-all for questions that do not fall neatly into the above categories.
In this paper, we focus our attention on closedvocabulary algebra, which constitutes approximately 45% of the questions.

Related Work
Most of the recent work on math questions has focused on open-vocabulary algebra problems, also known as math story problems. Benchmark datasets include Alg514 , AI2 (Hosseini et al., 2014), Illinois and Commoncore (Roy and Roth, 2015), and DRAW (Upadhyay and Chang, 2016). A common property of these datasets is that they have been curated such that any given question can be solved by a limiteddepth operator tree (AI2, Illinois, Commoncore) or a limited set of equation templates (Alg514 and DRAW). Because of this, it is feasible to use discriminative approaches Hosseini et al., 2014;Roy and Roth, 2015;Zhou et al., 2015;Koncel-Kedziorski et al., 2015;Mitra and Baral, 2016) that extract the quantities, featurize the question, and then perform a weighted search over the space of instantiated operator trees or equation templates. However it is not clear how one can extend these discriminative techniques to handle the complex semantics found in Examples 1 and 2. Very recently, (Matsuzaki et al., 2017) published a paper about their semantic parsing approach to pre-university math problems (harvested from Japanese exams rather than the Math SAT). It is challenging to do a direct comparison, since they report results only on the Japanese-language exams. They report end-to-end system results of 11% recall and 50% precision. (Shi et al., 2015) harvested a fairly diverse set of closed-vocabulary algebra problems (called Dolphin) from the web and provided the first results on that dataset. Here, we demonstrate how to handle the more complex discourse semantics and anaphoric phenomena found in Math SAT questions, and establish a new state-of-the-art result on the Dolphin benchmark.
4 System Overview Figure 1 shows a high-level view of our QA system. We will give a general overview in this section, and then explore more advanced concepts and examples in the subsequent section.

Intermediate Languages
Our QA system has two basic languages that mediate the transformation from the question passage to the answer: a syntactic language A and a semantic language B.
Syntactic language A has a constituent-style syntax convenient 4 for the tree transducers in our cascade. In Figure 2, we show an example. We have three basic node types: clauses, entities (these correspond to noun phrases), and details (these correspond to adjectival and adverbial modifiers). Each node has a table of fields (key-value pairs) that store child relationships and auxiliary information like tense and number. For brevity, this additional structure is omitted from Figure 2, but a more explicit visualization can be found in Figure 4 (top).
A program in semantic language B is a set of constraint declarations. For instance, the question from the form of the tree is unimportant, it will be convenient to use a more legible LISP-style format, e.g.
Every constraint in this program should be easily understandable, except for (proto m M ), which loosely means that M is the set of all possible values of m. In Section 5.1, we discuss the proto directive in more detail.

Syntactic Parsing
The first stage of our QA system parses the question passage into language A 5 . We implemented the parser as the backward application of an extended top-down tree-to-string (XTOPs) transducer 6 .
We refer the reader to (Maletti et al., 2009) for a theoretical presentation of XTOPs, and instead give a brief intuitive presentation of the device. An XTOPs transducer defines a top-down transformation from a tree language to a string language, via a set of stateful rewrite rules. For instance, rules (i) through (v) of Figure 4 can generate the string "a positive integer" from the A-tree pictured at the top of the figure, given start state qNP.
Given an XTOPs transducer M , we can parse string s through backward application of the transducer, i.e. compute the set of trees M −1 (s) that could have generated string s from the start 5 Recall: language A is the syntactic language described in the previous section. An example is shown in Figure 2. 6 We chose to implement the parsing step by engineering a transducer rather than using an off-the-shelf statistical parser. While we tried to retrofit a parser -e.g. as done by (Seo et al., 2015) -to serve our needs, it turned out to be somewhat more robust (and relatively simple) to engineer our own. state. Efficient backward application of XTOPs transducers is supported by packages like Tiburon (May and Knight, 2006).
Our XTOPs transducer has approximately 140 states and 550 engineered rules (approximately 200 of these rules are used for parsing formal mathematics and a subset of LaTeX). Most lexical rules are automatically generated on-the-fly from WordNet (Miller, 1995).

Compilation
We then compile the parses of the question passage, by running them forward through a cascade of bottom-up tree transducers (Engelfriet, 1975). Again we refer the reader to the literature (Maletti, 2011(Maletti, , 2014 for a theoretical presentation of bottom-up tree transducers, and use Figure 5 to provide intuition about the device. A bottomup tree transducer defines a transformation from a tree language to a (possibly different) tree language, via a set of stateful bottom-up rewrite rules.
In Figure 5, we show how this transformation works in the context of the semantic translation step, which uses a multi bottom-up transducer (MBOT) to map our syntactic language A into our semantic language B. There is a single state (indicated by a gray shaded rectangle) that has two children: (i) a return value, and (ii) a set of sideeffect statements.
The first rule application transforms "all values of m" into a return value of M (a new variable introduced to indicate the set of all values of variable m) and a side-effect (proto m M ), indicating that M equals the set of all possible values of m. The second rule application transforms "the sum of M " into a return value of (sum M ), and propagates upward the accumulated side-effects.
We implemented all three compilation steps from Figure 1 (anaphora resolution, semantic translation, and semantic analysis) as the forward application of a bottom-up tree transducer. Anaphora resolution resolves any nodes that refer to other nodes in the tree. Semantic translation translates syntactic language A into semantic language B. Semantic analysis type-checks the trees for internal consistency.

Interpretation
Finally, each derived B-tree is sent to an evaluator to obtain an answer. Our main evaluator is a wrapped version of Z3 (de Moura and Bjorner, 2008), a widely used Satisfiability Modulo Theories (SMT) solver. If it does not find an answer, we fall back to a numeric optimization solver similar to one used by (Seo et al., 2015).

Spotlights
Having provided a bird's eye view in the last section, we now spotlight some key details of our QA system.

Spotlight: Complex Aggregations
A core challenge of semantic parsing is how best to read complex semantic phenomena from a syntactic representation. Two such phenomena are superlatives and counting. GeoQuery (Zelle and Mooney, 1996) has examples 7 of these, as does 8 WebQuestions (Berant et al., 2013). Unfortunately, it is not clear how existing strategies for dealing with aggregative constructs (e.g. (Liang et al., 2011)) can be extended to the more complex multi-sentence questions found on the SATs. For instance, the basic semantics of Example 1 (enumerated in Section 1) is dispersed throughout the question passage, such that neither the maximization nor the minimization can be locally derived from some subtree of the dependency structure.
To deal with this challenge, we designed our semantic language B to decompose the semantics of aggregative constructs into order-independent atoms. Consider the following restatement of the semantics of Example 1: where proto(ż, Z) designates that a variableż should be treated as the prototype variable of a statement in set-builder notation, i.e. Z = {ż | ...}. We treat any other statement featuring prototype variableż as a constraint appearing on the right side of the set-builder statement. If there are multiple prototype statements, they are grouped into a single set-builder statement (as occurs witḣ x andẏ in our example).
The power of this decomposition is that it can be reconstructed piecemeal from an arbitrarily complex passage. The atomic statements can be interpreted locally in an arbitrary order, as in Figure 6, then synthesized into set-builder notation during evaluation.

Spotlight: Complex Anaphoric Phenomena
The anaphora resolution task (Ge et al., 1998) is typically defined at the lexical level. For instance, in the sentence "c is equal to its square," a traditional evaluation like the CoNLL-2011 Shared Task (Pradhan et al., 2011) would ask whether the string "its" is aligned to the string "c". These evaluations also assume that both the reference and the referent (a.k.a. antecedent) are contiguous strings in the text. Math SAT problems exhibit a host of new challenges that fall outside traditionally studied definitions of anaphora resolution: • One-to-many coreference 9 (One integer is 5 more than another. What is the sum of the numbers?): "The numbers" refers to two discontiguous strings: "one integer" and "another".
• Implicit set reference (Two numbers sum to 5. If the first is 2, what is the second?): "The second" implies a latent set that needs to be resolved (to "two numbers") in order to understand the sentence. This phenomenon is shown in Figure 7.
• Implicit clausal reference (If 7 is divided by 3, what is the remainder?): "The remainder" implies a latent clause that needs to be resolved (to "7 is divided by 3") in order to understand the sentence.
We address this broader class of anaphora by a two-pass process: 1. First, we introduce implicit sets and clauses when appropriate. For instance, implicit sets are introduced for superlative and ordinal constructions, while implicit clauses are introduced for functional nouns like "remainder." In Figure 7, these implicits are depicted as bracketed phrases (i.e. [of a set]).
2. Anaphora resolution then proceeds as a bottom-up tree-labeling process, shown in Figure 7. For each subtree, a resolution function ρ partially maps subtree entities to subtree nodes. Note that ancestors can overwrite the resolutions of their descendants. This occurs in the second example of Figure 7, where the implicit set E7 is initially resolved to implicit set E4, but is later resolved to the entity E1 ("two numbers") once it comes into scope.
In our initial system, the resolution function ρ was Figure 7: Bottom-up anaphora resolution in our QA system. For convenience, we show the correspondence of the nodes of our syntactic parse (top) to the original question passage (bottom). In the parse tree, "E" is an abbreviation for ENTITY. engineered heuristically. We later replaced this with a learned version (by using our system to generate training data). Due to space considerations, details are omitted.

Results with the Unweighted Nondeterministic Cascade
In the basic cascade from Section 4, the number of trees passed from module to module can expand, but it can also contract (for instance, in the semantic translation step, there can be multiple ways of translating a parse, or none at all). This allows the QA system to disambiguate question passages by eliminating parses for which there is no consistent semantics. On the subset of the Kaplan questions for which at least one parse exists, the average number of trees after the parsing step is 7.5. The average number of trees after the semantic analysis step goes down to only 2.4. At that point, obviously we need to choose some priority in which to feed these finalized programs to the evaluation module. Using a simple heuristic (process smaller programs first), we obtain 70.2% recall and 95.8% precision on the Kaplan closed-vocabulary algebra questions 10 . This high precision can be partially attributed to the fact that most SAT questions are multiplechoice (thus we can sequentially evaluate the finalized programs until we find a viable answer). We do not have that luxury on the Dolphin dataset, a set of direct-answer algebra questions curated by Microsoft Research (split into a development set of 374 questions and a test set of 1504 ques-10 Recall and precision numbers are computed over the entire set of questions, regardless of whether they have a valid parse. tions). On the subset of the development questions for which at least one parse exists (90.3% of the questions), the average number of trees after the parsing step is 4.3. The average number of trees after the semantic analysis step goes down to 1.5. Our basic system obtains 66.3% recall on the development questions. Naturally the precision is not as high as on the multiple choice questions, but surprisingly we still obtain 85.5% precision, even with an unweighted cascade.

Introducing a Parse Ranker
Most of this precision loss is due to legitimate parse ambiguity that cannot be resolved through semantic interpretability alone. Rather, the disambiguation requires some additional pragmatic convention. Consider the example: "When the reciprocal of three times a number is subtracted from the reciprocal of the number, the result is one sixth. Find the number." By interpreting "the reciprocal of three" as 1 3 , the meaning of this question becomes "When 1 3 times a number is subtracted from the reciprocal of the number, the result is one sixth. Find the number." This is not however the most human-intuitive interpretation of the question. Somehow the system must identify the pragmatic cues that cause humans to disprefer this interpretation.
To identify these cues, we insert a parse ranking module between the parsing module and the anaphora resolution module (see Figure 1 for a reminder of the system components). The goal of the parse ranker is to associate a lower cost to "more intuitive" interpretations when there are multiple plausible syntactic interpretations. The  Table 2: Results on the Dolphin question sets. The increase in recall is statistically significant with a P -value < 0.01.
rest of the cascade propagates these costs. Similar to existing work, e.g. (Charniak and Johnson, 2005), we implement the cost function as an L 1regularized logistic regression model.
Adding the trained parse ranker module improves performance on the Dolphin development set to 75.7% recall and 97.3% precision (from 66.3% recall and 85.5% precision).

Final Results
Results from our final system are shown in Table 1 (for the closed-vocabulary algebra subsets of our math SAT question sets) and Table 2 (for the Dolphin question sets). EUCLID generalizes reasonably well to the blind SAT questions, achieving approximately 60% of the system's recall on the training questions, at a precision of approximately 91%. To give a sense of the extent of the generalization from training to test, Table 3 offers a couple of correctly answered questions from the blind 11 McGraw-Hill set, plus their closest analog in the training questions (by edit distance). The performance on the blind test sets (including all questions, not just closed-vocabulary algebra) corresponds to an SAT score of approximately 350 (out of 800). A random-guessing baseline has an expected score of 200. If x and y are positive even integers different positive such that a < b < c integers and and a + b + c = 60, 3x + y = 17, the then the greatest difference between possible value of c is the largest possible value of y and the smallest possible value of x is development (blind) training failed to parse 50% failed to map parse into a semantic program 24% failed to produce an answer from any semantic program 18% produced an incorrect answer 8% valid parse. Roughly a quarter of the questions had at least one valid parse, but none of these resulted in a semantic program. 18% of the questions compiled into at least one semantic program, but none of these produced an answer when fed to the interpreter. 8% of the questions compiled into a semantic form that produced an incorrect answer.
Besides the Math SAT datasets, EUCLID also has state-of-the-art performance on the public Dolphin question set, achieving an absolute recall improvement 12 of nearly 5% with a small loss in precision. This raises the state-of-the-art F 1 -score on this data set from 73.9% to 77.0%.  Table 5: Snapshot of early progress across several subgenres of the Math SAT. In our early stages, we are hillclimbing on the hosseini datasets from (Hosseini et al., 2014), which are simpler than the open-vocabulary algebra questions from the Math SAT.

Towards a Broad-Coverage SAT solver
This paper reports on the first steps of a longerterm initiative to build a unified system that can pass the math SAT. We have made some preliminary forays into extending the system to handle the more complex subdomains described in Section 2, namely open-vocabulary algebra and geometry. Key research challenges presented by these new domains are: • Mapping into richer semantic languages: The math story problems of open-vocabulary algebra require languages that reason about state change and can introduce assumptions not explicitly represented in the text.
• Robustly synthesizing diagram and text information: For geometry questions, we are building on key early work in this area performed by (Seo et al., 2014(Seo et al., , 2015.
• Extending the system in a scalable way: Writing new transducer rules for each new domain is not a sustainable way to extend our system. We are exploring how to use natural language to "program" our system, e.g. by automatically inducing transducer rules for paraphrase text.
A snapshot of our current progress is shown in Table 5.

Discussion
In the process of creating a QA system for math SAT questions, this project has yielded several general strategies for beyond-sentential semantic parsing. For instance: • One can modularize the parser as a cascade of tree transducers, allowing uncertainty about anaphora resolution and lexical interpretation to be propagated until it can be confidently resolved, sometimes as late as program interpretation (see Section 6).
• One can atomize complex semantic phenomena (e.g. aggregrative constructs) into small order-independent pieces. This allows a simpler transformation from a syntactic form, because these atoms can be locally recognized, incrementally composed, and globally reconstituted into a structured semantics (see Section 5.1).
• One can reframe bread-and-butter NLP tasks (e.g. anaphora resolution) to fit better within (and take advantage of) the context of the transducer cascade (see Section 5.2) An important focus of this paper has been issues of representation, namely how to develop intermediate structured languages that facilitate the automatic transformation of question discourse into a response. Because we can use the resulting QA system to generate gold intermediate trees for any correctly answered question in our dataset, one way to view this work is as a data annotation project. One distinguishing advantage is that our intermediate languages come with a "proof of usefulness." They are not designed based on speculative utility -rather, they have already proven useful in the context of a functioning QA system.


An End-to-End Trainable Neural Network Model with Belief Tracking for Task-Oriented Dialog

We present a novel end-to-end trainable neural network model for task-oriented dialog systems. The model is able to track dialog state, issue API calls to knowledge base (KB), and incorporate structured KB query results into system responses to successfully complete task-oriented dialogs. The proposed model produces well-structured system responses by jointly learning belief tracking and KB result processing conditioning on the dialog history. We evaluate the model in a restaurant search domain using a dataset that is converted from the second Dialog State Tracking Challenge (DSTC2) corpus. Experiment results show that the proposed model can robustly track dialog state given the dialog history. Moreover, our model demonstrates promising results in producing appropriate system responses, outperforming prior end-to-end trainable neural network models using per-response accuracy evaluation metrics.

Introduction
Task-oriented spoken dialog system is a prominent component in today's virtual personal assistants, which enable people to perform everyday tasks by interacting with devices via voice input. Traditional task-oriented dialog systems have complex pipelines, with a number of independently developed and modularly connected components. There are usually separated modules in a pipeline for natural language understanding (NLU), dialog state tracking (DST), dialog management (DM), and natural language generation (NLG) [1,2,3,4]. One limitation with such pipeline approach is that it is inherently hard to adapt a system to new domains, as all these modules are trained and finetuned independently. Moreover, errors made in upper stream modules may propagate to downstream components, making it tedious to identify and track the source of error [5].
To address these limitations, efforts have been made recently in designing end-to-end frameworks for task-oriented dialogs. Wen et al. [6] proposed an end-to-end trainable neural network model with modularly connected neural networks for each system component. Zhao and Eskenazi [5] introduced an end-to-end reinforcement learning framework that jointly performs dialog state tracking and policy learning. Li et al. [7] proposed an end-to-end learning framework that leverages both supervised and reinforcement learning signals and showed promising dialog modeling performance. Such end-toend trainable neural network models can be optimized directly towards the final system objective functions (e.g. task success rate) and thus ameliorate the challenges of credit assignment and online adaptation [5].
In this work, we present an end-to-end trainable neural network model for task-oriented dialog that applies a unified net-work for belief tracking, knowledge base (KB) operation, and response creation. The model is able to track dialog state, interface with a KB, and incorporate structured KB query results into system responses to successfully complete task-oriented dialogs. We show that our proposed model can robustly track belief state given the dialog history. Our model also demonstrates promising performance in providing appropriate system responses and conducting task-oriented dialogs compared to prior end-to-end trainable neural network models.

Dialog State Tracking
In spoken dialog systems, dialog state tracking, or belief tracking, refers to the task of maintaining a distribution over possible dialog states which directly determine the systems actions. Dialog state tracker is a core component in many state-of-theart task-oriented spoken dialog systems [6,8]. Conventional approaches for DST include using rule-based systems and generative methods that model the dialog as a dynamic Bayesian network [9]. Discriminative approaches using sequence models such as CRF [10] or RNN [11,12] address the limitation of generative models with the flexibility in exploring arbitrary features [13] and achieve state-of-the-art DST performance.

End-to-End Task-Oriented Dialog Models
Conventional task-oriented dialog systems typically require a large number of domain-specific rules and handcrafted features, which make it hard to extend a good performing model to new application domains. Recent approaches to task-oriented dialog cast the task as a partially observable Markov Decision Process (POMDP) [4] and use reinforcement learning for online policy optimization by interacting with users [14]. The dialog state and system action spaces have to be carefully designed in order to make the reinforcement policy learning tractable [4].
With the success of end-to-end trainable neural network models in non-task-oriented chit-chat dialog settings [15,16], efforts have been made in carrying over the good performance of end-to-end trainable models to task-oriented dialogs. Wen et al. [6] proposed a neural network based model that is endto-end trainable yet still modularly connected. The model has separated modules for intent estimation, belief tracking, policy learning, and response generation. Our model, on the other hand, use a unified network for belief tracking, KB operation, and response generation, to fully explore knowledge that can be shared among different tasks. Bordes and Weston [17] recently proposed modeling the dialog with a reasoning approach using end-to-end memory network. Their model selects a best system response directly from a list of response candidates without explicitly tracking dialog state. Comparing to this approach, our

Food slot Price slot
LSTMD, sk+1 LSTMD, sk-1 Figure 1: System architecture of the proposed end-to-end trainable neural network model for task-oriented dialog. model tracks dialog state over the sequence of turns explicitly, as it is shown in [18] that robust dialog state tracking is likely to boost success rate in task completion. Moreover, when generating final system response, instead of letting the model to select a final response directly from a large pool of candidate responses, we let our model to select skeletal sentence structure from a short list of candidates and then replace the delexicalised tokens with the state tracking outputs. This method will help to reduce the number of training samples required [11] and make the model more robust to noises in dialog state.

Proposed Method
We model task-oriented dialog as a multi-task sequence learning problem, with components for encoding user input, tracking belief state, issuing API calls, processing KB results, and generating system responses. The model architecture is as shown in Figure 1. Sequence of turns in a dialog is encoded using LSTM [19] recurrent neural networks. Conditioning on the dialog history, state of the conversation is maintained in the LSTM state. The LSTM state vector is used to generate: (1) a skeletal sentence structure by selecting from a list of delexicalised system response candidates, (2) a probability distribution of values for each slot in belief tracker, and (3) a pointer to an entity in the retrieved KB results that matches the user's query. The final system response is generated by replacing the delexicalised tokens with the predicted slot values and entity attribute values. Each model component is described in detail in below sections.

Utterance Encoding
Utterance encoding here refers to encode a sequence of words into a continuous dense vector. Popular methods include using bag-of-means on word embeddings and RNNs [20,21]. We use bidirectional LSTM to encode the user input to an utterance vector. Let U k = (w1, w2, ..., wT k ) be the user input at the kth turn with T k words. The user utterance vector U k is represented 1 are the last forward and backward utterance-level LSTM states at kth turn.

Belief Tracking
Belief tracking, or dialog state tracking, maintains and adjusts the state of a conversation, such as user's goals, by accumu-late evidence along the sequence of a dialog. After collecting new evidence from a user's input at turn k, the neural dialog model updates the probability distribution P (S m k ) over candidate values for each slot type m ∈ M . For example, in restaurant search domain, the model maintains a multinomial probability distribution over each of user's goals on restaurant area, food type, and price range. At turn k, the dialog-level LSTM (LSTMD) updates its hidden state s k and use it to infer any updates on user's goals after taking in the user input encoding U k and KB indicator I k (to be described in section below).
where SlotDistm is a multilayer perceptron (MLP) with softmax activation function over the slot type m ∈ M .

Issuing API Calls
Conditioning on the state of the conversation, the model may issue an API call to query the KB based on belief tracking outputs. A simple API call command template is firstly generated by the model. The final API call command is produced by replacing the slot type tokens in the command template with the best hypothesis for each of the goal slot from the belief tracker. In restaurant search domain, a simple API call command template can be "api call area f ood pricerange ", and the slot type tokens are to be replaced with the belief tracker outputs to form the final API call command "api call west italian dontcare".

KB Results Processing
Once the neural dialog model receives the KB query results, it suggests options to users by selecting entities from the returned list. Instead of treating KB results as unstructured text (more specifically, as user utterances as in [17,22,23]) and processing them with machine reading comprehension approach, we treat KB results as a list of structured entities and let the model to select appropriate entity pointers. Outputs from KB search or database query typically have well defined structures, with entity attributes associated with entity index. Other than letting the model to learn such entity-attribute association purely from the training dialog corpus as in [17,22,23], we keep such structural information in our system and let the model to learn to select proper entity pointers from a ranked list.
At turn k of a dialog, a binary KB indicator I k is passed to the neural dialog model. This indicator is decided by the number of retrieved entities from the last API call and the current entity pointer. When the system is in a state to suggest an entity to user, if a zero value I k is received, the model is likely to inform user the unavailability of entity matching the current query. Otherwise if I k has a value of one, the model will likely pick an entity from the retrieved results based on the updated probability distribution of the entity pointer P (E k ): where EntityPointerDist is an MLP with softmax activation.

System Response Generation
At kth turn of a dialog, a skeletal sentence structure R k is selected from a list of delexicalised response candidates. The final system response is produced by replacing the delexicalised tokens with the predicted slot values and entity attribute values. For example, replacing f ood to italian, and replacing R name to prezzo as in Figure 1.
where ResponseDist is an MLP with softmax activation.

Model Training
We train the neural dialog model by finding the parameter set θ that minimize the cross-entropy of the predicted and true distributions for goal slot labels, entity pointer, and delexicalised system response jointly: where λs are the linear interpolation weights for the cost of each system output. S m k * , E * k , and R * k are the ground truth labels for each task at the kth turn.

Alternative Model Designs
The model architecture ( Figure 1) described above assumes that the hidden state of the dialog-level LSTM implicitly captures the complete state of the conversation, i.e. the user goal estimation and the previous system actions. Intuitively, the model is likely to provide a better response if it is informed about the goal slot value estimations explicitly and is aware of its previous responses made to the user. Thus, we design and evaluate a few alternative model architectures to verify such assumption: (1) Model with previously emitted delexicalised system response connected back to dialog-level LSTM state: (2) Model with previously emitted slot labels connected back to dialog-level LSTM state:

Dataset
We use data from DSTC2 [24] for our model evaluation. This challenge is designed in the restaurant search domain. Bordes and Weston [17] transformed the original DSTC2 corpus by adding system commands and removing the dialog state annotations. This transformed corpus contains additional API calls that the system would make to the KB and the corresponding KB query results. In this study, we combine the original DSTC2 corpus and this transformed version by keeping the dialog state annotations and adding the system commands (API calls). We can thus perform more complete evaluation of our model's capability in tracking the dialog state, processing KB query results, and conducting complete dialog. Statistics of this dataset is summarized in the Table 1.

Model Configuration and Training
We perform mini-batch model training with batch size of 32 using Adam optimization method [25]. Regularization with dropout is applied to the non-recurrent connections [26] during model training with dropout rate of 0.5. We set the maximum norm for gradient clipping to 5 to prevent exploding gradients. Hidden layer sizes of the dialog-level LSTM and the utterance-level LSTM are set as 200 and 150 respectively. Word embeddings of size 300 are randomly initialized. We also explore using pre-trained word vectors [27] that are trained on Google News dataset to initialize the word embeddings.

Results and Analysis
Similar to the evaluation methods used in [17,22,23], we evaluate the task-oriented dialog model in a ranking setting. We report the prediction accuracy for goal slot values, entity pointer, delexicalised system response, and final system response which has the delexicalised tokens replaced by predicted values.
We first experiment with different text encoding methods and recurrent model architectures to find best performing model. Table 2 shows the evaluation results of models using different user utterance encoding methods and different word embedding initialization. Bidirectional LSTM (Bi-LSTM) shows clear advantage in encoding user utterance comparing to bagof-means on word embedding (BoW Emb) method, improving the joint goal prediction accuracy by 4.6% and the final system response accuracy by 1.4%. Using pre-trained word vectors (word2vec) boosts the model performance further. These results show that the semantic similarities of words captured in the pre-trained word vectors are helpful in generating a better representation of user input, especially when the utterance contains words or entities that are rarely observed during training. Table 3 shows the evaluation results of different recurrent model architectures. The Hierarchical LSTM model in Table  3 refers to the last model in Table 2. Models in row 2 to 4 of Table 3 refer to the three recurrent model architectures dis- cussed in section 3.7. As illustrated in these results, models that conditioned on previous emitted labels in generating system response achieve lower prediction accuracy across all of the four evaluation metrics. These observations are contrary to our intuition and analysis made in previous section. We believe the degraded performance is mainly due to the data sparsity issue of the dataset that used in the experiment. Given a certain dialog context, there might be multiple system response candidates that can be used to generate a suitable final response. With limited number of training samples, the model is likely to overfit the training set and not to generalize the strong modeling capacity well during inference. The overfitting issue might be less of a problem in the word-by-word response generation setting, and this is to be studied further in our future work. As observed in Table 3, the joint goal tracking performance is directly related to final response prediction accuracy. To further investigate our model's capability on belief tracking, we conduct error breakdown analysis for each goal slot. We compare our model to two other recently proposed belief tracking models, an RNN based model [11] and the Neural Belief Tracker [28], in the setting of only using live ASR hypothesis as model input. As the results show in Table 4, our system achieves promising belief tracking performance comparable to the state-of-the-art systems. Finally, we report the model performance in producing final system responses and compare it to other published results following the per-response accuracy metric used in prior work. Even though using the same evaluation measurement, our model is designed with slightly different settings comparing to other published models in Table 5. Instead of using additional matched type features (i.e. KB entity type feature for each word, e.g. whether a word is a food type or area type, etc.) as in [17,23], we use user's goal slots at each turn that are mapped from the original DSTC2 dataset as additional supervised signals in our model. Moreover, instead of treating KB query results as unstructured text, we treat them as structured entities and let our model to pick the right entity by selecting the most appropriate entity pointer. Our proposed model successfully predicts 52.8% of the true system responses, outperforming prior end-to-end trainable neural dialog systems. Table 5: Performance of the proposed model in per-response accuracy comparing to previous approaches.

Model
Per-res Accuracy Memory Networks [17] 41.1 Gated Memory Networks [29] 48.7 Sequence-to-Sequence [22] 48.0 Query-Reduction Networks [23] 51.1 Hierarchical LSTM 52.8 To further understand the prediction errors made by our model, we conduct human evaluation by inviting 10 users to evaluate the appropriateness of the responses generated by our system. While some of the errors are made on generating proper API calls due to the errors in dialog state tracking results, we also find quite a number of responses that are considered appropriate by our judges but do not match to the reference responses in the test set. For example, there are cases where our system directly issues the correct API call (e.g. "api call south italian expensive") based on user's inputs, instead of asking user for confirmation of a goal type (e.g. "Did you say you are looking for a restaurant in the south of town?") as in the reference corpus. By taking such factors into consideration, our system was able to generate appropriate responses in 73.6% of the time based on feedback from the judges. These results show that the per-response accuracy evaluation metric may not well correlate with human judgments [30], and better dialogs evaluation measurements should to be further explored.

Conclusions
In this work, we propose a novel end-to-end trainable neural network model for task-oriented dialog systems. The model is able to track dialog belief state, interface with knowledge bases by issuing API calls, and incorporate structured query results into system responses to successfully complete task-oriented dialogs. In the evaluation in a restaurant search domain using a converted dataset from the second Dialog State Tracking Challenge corpus, our proposed model shows robust performance in tracking dialog state over the sequence of dialog turns. The model also demonstrates promising performance in generating appropriate system responses, outperforming prior end-to-end trainable neural network models.


extfieldsofstudy:
Computer Science
Biology

text:
Chemical Boltzmann Machines

How smart can a micron-sized bag of chemicals be? How can an artificial or real cell make inferences about its environment? From which kinds of probability distributions can chemical reaction networks sample? We begin tackling these questions by showing four ways in which a stochastic chemical reaction network can implement a Boltzmann machine, a stochastic neural network model that can generate a wide range of probability distributions and compute conditional probabilities. The resulting models, and the associated theorems, provide a road map for constructing chemical reaction networks that exploit their native stochasticity as a computational resource. Finally, to show the potential of our models, we simulate a chemical Boltzmann machine to classify and generate MNIST digits in-silico.

Introduction
To carry out complex tasks such as finding and exploiting food sources, avoiding toxins and predators, and transitioning through critical life-cycle stages, singlecelled organisms and future cell-like artificial systems must make sensible decisions based on information about their environment [1,2]. The small volumes of cells makes this enterprise inherently probabilistic: environmental signals and the biochemical networks within the cell are noisy, due to the stochasticity inherent in the interactions of small, diffusing molecules [3,4,5]. The small volumes of cells raises questions not only about how stochasticity influences circuit function, but also about how much computational sophistication can be packed into the limited available space.
Perhaps surprisingly, neural network models provide an attractive architecture for the types of computation, inference, and information processing that cells must do. Neural networks can perform deterministic computation using circuits that are smaller and faster than boolean circuits composed of AND, OR, and NOT gates [6], can robustly perform tasks such as associative recall [7], and naturally perform Bayesian inference [8]. Furthermore, the structure of biochemical networks, such as signal transduction cascades [9,1,10] and genetic regulatory networks [11,12,13,14,15], can map surprisingly well onto neural network architectures. Chemical implementations of neural networks and related machine In a micron-scale environment, molecular counts are low and a real (or synthetic) cell will have to respond to internal and environmental cues. Probabilistic inference using chemical Boltzmann machines provides a framework for how this may be achieved.
learning models have also been proposed [16,17,18,19,20], and limited examples demonstrated [21,22,23,24], for synthetic biochemical systems. Previous work on biochemical neural networks and biochemical inference invoked models based on continuous concentrations of species representing neural activities. Such models are limited in their ability to address questions of biochemical computation in small volumes, where discrete stochastic chemical reaction network models must be used to account for the low molecular counts. The nature of biochemical computation changes qualitatively in this context. In particular, stochasticity has been widely studied in genetic regulatory networks [25], signaling cascades [26], population level bet hedging in bacteria [27], and other areas [28,29] -where the stochasticity is usually seen as a challenge limiting correct function, but is occasionally also viewed as a useful resource [30]. Our work falls squarely in the latter camp: we attempt to exploit the intrinsic stochastic fluctuations of a formal chemical reaction network (CRN) to build natively stochastic samplers by implementing a stochastic neural network. This links to efforts to build natively stochastic hardware for Bayesian inference [31,32] and to the substantial literature attempting to model, and find evidence for, stochastic neural systems capable of Bayesian inference [33,34].
Specifically, we propose CRNs that implement Boltzmann machines (BMs), a flexible class of Markov random fields capable of generating diverse distributions and for which conditioning on data has straightforward physical interpretations [8,35]. BMs are an established model of probabilistic neural networks due to their analytic tractability and identity with spin systems in statistical physics [36] and Hopfield networks in computer science [7]. These networks have been studied extensively and used in a wide range of applications including image classification [37] and video generation [38]. We prove that CRNs can implement BMs and that this is possible using detailed balanced CRNs. Moreover, we show that many of the attractive features of BMs can be applied to our CRN constructions such as inference, a straightforward learning rule and scalability to real-world data sets. We thereby introduce the idea of a chemical Boltzmann machine (CBM), a chemical system capable of exactly or approximately performing inference using a stochastically sampled high-dimensional state space, and explore some of its possible forms.

Boltzmann Machines (BMs):
Boltzmann machines are a class of binary stochastic neural networks, meaning that each node randomly switches between the values 0 and 1 according to a specified distribution. They are widely used for unsupervised machine learning because they can compactly represent and manipulate high-dimensional probability distributions. Boltzmann machines provide a flexible machine learning architecture because, as generative models, they can be used for a diverse set of tasks including data classification, data generation, and data reconstruction. Additionally, the simplicity of the model makes them analytically tractable. The use of hidden units (described below) allows Boltzmann machines to represent high order correlations in data. Together, these features make Boltzmann machines an excellent starting point for implementing stochastic chemical computers.
Fix a positive integer N ∈ Z >0 . An N -node Boltzmann machine (BM) is specified by a quadratic energy function E : where θ i ∈ R is the bias of node i, and w ij = w ji ∈ R is the weight of the unordered pair (i, j) of nodes, with w ii = 0. One may specify a BM architecture, or graph topology, by choosing additional weights w ij that are to be set to 0. In this paper, we will use N (i) = {j s.t. w ij = 0} to denote the neighborhood of i. From a physical point of view, we are implicitly using temperature units k B T for energy, which we will continue to do throughout this paper. A BM describes a distribution P (x) over state vectors x = (x 1 , . . . , x N ) ∈ {0, 1} N , Nodes of a BM are often partitioned into sets V and H of visible and hidden, respectively. Nodes in V represent data, and auxiliary nodes in H allow more complex distributions to be represented in the visible nodes. An implementation of a BM is a stationary stochastic process that samples from this distribution in the steady state. A BM can be implemented in silico using the Gibbs sampling algorithm [39], which induces a discrete time Markov chain (DTMC) on the state space {0, 1} N in such a way that the stationary distribution of this Markov chain corresponds to the distribution P (x). In each round, one node i ∈ {1, . . . , n} is chosen at random for update. For any two adjacent configurations x and x which differ only at node i -i.e., x i = 1−x i and x j = x j for all j = i -we set the transition probabilities T x→x of the DTMC so that Any function T x→x can be chosen so long as (3) is satisfied. One common choice A Boltzmann machine is also an inference engine. One can do inference on P (x) by conditioning on the values of a subset of the nodes. Suppose nonempty node subsets U and Y form a partition of the nodes {1, 2, . . . , N }, and fix u ∈ {0, 1} U . To obtain samples from P (y | u) where y ∈ {0, 1} Y , one clamps every node i ∈ U to the state u i while running Gibbs sampling, i.e., one does not allow these nodes to update. Clamping nodes to an input state is the same as specifying the input data for a statistical model. Steady state samples y ∈ {0, 1} Y of this procedure are draws from the distribution P (y | u).
Boltzmann machines can be used to learn a generative model from unlabeled data. After specifying the architecture, one then proceeds to find the weights, w ij , and biases, θ i , that maximize the likelihood of the observed data according to the model, using gradient descent from a random initial parameterization. This reduces to a very simple two-phase learning rule where weights on active edges are strengthened in a "wake phase" during which the network is clamped to observed data and are weakened in a "sleep phase" during which the network runs free [8,35]. Given a target distribution Q(x), this gradient descent corresponds to calculating the gradient of the Kullback-Leibler divergence from P to , with respect to the parameters θ i and w ij : (4) where · P and · Q denote expected values with respect to the distributions P and Q respectively. When hidden units are present, the distribution Q (which is defined on visible units only) is extended to hidden units based on clamping the visible units according to Q and using the conditional distribution P (y|u) to determine the hidden units.

Chemical Reaction Networks (CRNs):
Fix a finite set S = (S 1 , S 2 , . . . , S M ) of M species. A reaction r is a formal chemical equation abbreviated as r = µ r → ν r where µ r , ν r ∈ N S are the stoichiometric coefficient vectors for the reactant and product species respectively, and N = Z ≥0 . A reaction rate constant, k r ∈ R >0 , is associated with each reaction. In this paper, we define a chemical reaction network (CRN) as a triple C = (S, R, k) where S is a finite set of species, and R is a set of reactions, and k is the associated set of reaction rate constants.
We will denote chemical species by capital letters, and their counts by lower case letters; e.g., s 1 denotes the number of species S 1 . Thus the state of a stochastic CRN is described by a vector on a discrete lattice, s = (s 1 , s 2 . . . s M ) ∈ N S . The dynamics of a stochastic CRN are as follows [40]. The probability that a given reaction occurs per unit time, called its propensity, is given by Each time a reaction fires, state s changes to state s + ∆ r , where ∆ r = ν r − µ r is called the reaction vector, and the propensity of each reaction may change.
Viewed from a state space perspective, stochastic CRNs are continuous time Markov chains (CTMCs) with transition rates and thus their dynamics follow dP (s, t) dt = s =s where P (s, t) is the probability of a state with counts s at time t. Equivalently, they are governed by the chemical master equation, A stationary distribution π(s) may be found by solving dP (s,t) dt = 0 simultaneously for all s; in general, it need not be unique, and even may not exist. Given an initial state s 0 , π(s) = P (s, ∞) is unique if it exists. For that initial state, the reachability class Ω s0 ⊆ N M is the maximal subset of the integer lattice accessible to the CRN via some sequence of reactions in R. We will specify a CRN and a reachability class given an initial state as a shorthand for specifying a CRN and a set of initial states with identical reachability classes.

Detailed Balanced Chemical Reaction Networks:
A CTMC is said to satisfy detailed balance if there exists a well-defined function of state s, E(s) ∈ R, such that for every pair of states s and s , the transition rates R s→s and R s →s are either both zero or If the state space Ω is connected and the partition function Z = s∈Ω e −E(s) is finite, then the steady state distribution π(s) = 1 Z e −E(s) is unique, and the net flux between all states is zero in that steady state.
There is a related but distinct notion of detailed balance for a CRN. An equilibrium chemical system is constrained by physics to obey detailed balance at the level of each reaction. In particular, for a dilute equilibrium system, each species S i ∈ S has an energy G[S i ] ∈ R that relates to its intrinsic stability, and where ∆ i r + is the ith component of ∆ r + = ν r + − µ r + , and we have defined Any CRN for which there exists a function G satisfying (11) is called a detailed balanced CRN. To see that the CTMC for a detailed balanced CRN also itself satisfies detailed balance in the sense of (10), let s = s + ∆ r+ and note that (6) and (11) imply that for all reactions r + . Here, G(s) is a well-defined function of state s (the free energy) that can play the role of E in (10). If there are multiple reactions that bring s to s , they all satisfy (12), and therefore the ratio R s→s /R s →s satisfies (10) and the CTMC satisfies detailed balance. It is possible to consider non-equilibrium CRNs that violate (11). Such systems must be coupled to implicit reservoirs of fuel molecules that can drive the species of interest into a non-equilibrium steady state [41,42,43]. Usually -but not always [44,45] -the resultant Markov chain violates detailed balance. In Section 3.1, we shall consider a system that exhibits detailed balance at the level of the Markov chain, but is necessarily non-equilibrium and violates detailed balance at the detailed chemical level.
Given an initial condition s 0 , a detailed balanced CRN will be confined to a single reachability class Ω s0 . Moreover, from the form of G(s), the stationary distribution π(s) on Ω s0 of any detailed balanced CRN exists, is unique, and is a product of Poisson distributions restricted to the reachability class [46], with the partition function Z = s ∈Ωs 0 e −G(s ) dependent on the reachability class. Note that this implies that the partition function is always finite, even for an infinite reachability class.

Clamping and Conditioning with Detailed Balanced CRNs:
In a Boltzmann machine that has been trained to generate a desired probability distribution when run, inference can be performed by freezing, also known as clamping, the values of known variables, and running the rest of the network to obtain a sample; this turns out to exactly generate the conditional probability. A similar result holds for a subclass of detailed balanced CRNs that generate a distribution, for an appropriate notion of clamping in a CRN. Imagine a "demon" that, whenever a reaction results in a change in the counts of one of the clamped species, will instantaneously change it back to its previous value. If every reaction is such that either no clamped species change, or else every species that changes is clamped, then the demon is effectively simply "turning off" those reactions. More precisely, consider a CRN, C = (S, R, k). We will partition the species into two disjoint groups Y = S f ree and U = S clamped , where S f ree will be allowed to vary and S clamped will be held fixed. We will define free reactions, R f ree , as reactions which result in neither a net production nor a net consumption of any clamped species. Similarly, clamped reactions, R clamped are reactions which change the counts of any clamped species. The clamped CRN will be denoted C| U =u to indicate the the species U i ∈ U have been clamped to the values u i . The clamped CRN is defined by C| U =u = (S, R f ree , k f ree ), that is, the entire set of species along with the reduced set of reactions and their rate constants. In the clamped CRN it is apparent that the clamped species will not change from their initial conditions because no reaction in R f ree can change their count. However, these clamped species may still affect the free species catalytically. If the removed reactions, R clamped , never change counts of non-clamped species, then C| U =u is equivalent to the action of the "demon" imagined above.
We use equation 13 to prove that clamping a detailed balanced CRN is equivalent to calculating a conditional distribution, and to show when the conditional distributions of a detailed balanced CRN will be independent. Together, these theorems provide guidelines for devising detailed balanced CRNs with interesting (non-independent) conditional distributions and for obtaining samples from these distributions via clamping.
We will need one more definition. Let C be a detailed balanced CRN with reachability class Ω s0 for some initial condition s 0 = (u 0 , y 0 ). Let Γ s0 be the reachability class of the clamped CRN C| U =u0 with species U clamped to u 0 and species Y free. We say clamping preserves reachability if where In other words, clamping preserves reachability if, whenever a state s = (u 0 , y) is reachable from s 0 by any path in C, then it is also reachable from s 0 by some path in C| U =u0 that never changes u. Theorem 1. Consider a detailed balanced CRN C = (S, R, k) with reachability class Ω s0 from initial state s 0 . Partition the species into two disjoint sets Let the projection of s 0 onto U and Y be u 0 and y 0 . The conditional distribution P (y | u) implied by the stationary distribution π of C is equivalent to the stationary distribution of a clamped CRN, C| U =u starting from initial state s 0 with u 0 = u, provided that clamping preserves reachability.
Then, the conditional probability distribution of the unclamped CRN is given by Simply removing pairs of forwards and backwards reactions will preserve detailed balance for unaffected transitions, and hence the clamped system remains a detailed balanced CRN with the same free energy function. We then readily see that the clamped CRN's stationary distribution, π c (y|u) is given by The original CRN and the clamped CRN do not need to have the same initial conditions as long as the initial conditions have the same reachability classes. However, even if the two CRNs have the same initial conditions, it is possible that the clamping process will make some part of Ω Y s0|U =u0 inaccessible to C| U =u , in which case this theorem will not hold.
Theorem 2. Assume the reachability class of a detailed balanced CRN can be expressed as the product of subspaces, Ω s0 = L j=1 Ω j s0 . Then the steady-state distributions of each subspace will be independent for each product space: π(s) = L j=1 π j (s j ), where s = (s 1 , . . . , s L ) and π j is the distribution over Ω j s0 . Proof.
If Ω s0 is decomposable into a product of subspaces Ω j s0 , with j = 1...L, then each subspace involves disjoint sets of species Y j = {S j 1 , . . . , S j Mj }. In this case the steady-state distribution of a detailed balanced CRN can be factorized due to the simple nature of G(s) given by eq. (12): where s j = (s j 1 , s j 2 , . . . , s j Mj ) is state of the set of species within the subspace j.
The product form π(s) means that species from separate subspaces Ω j s0 are statistically independent. To develop non-trivial conditional probabilities for the states of different species, therefore, it is necessary either to use a non-detailed balanced CRN by driving the system out of equilibrium, or to generate complex interdependencies through conservation laws that constrain reachability classes and "entangle" the state spaces for different species. We explore both of these possibilities in the following sections.

Direct Implementation of a Chemical Boltzmann Machine (DCBM):
We first consider the most direct way to implement an N -node Boltzmann machine with a chemical system. Recall that a BM has a state space Ω BM = {0, 1} N and an energy function  . Although there are 2 2N states in which where s is any state with x OF F i = 1, and s has x ON i = 1 but is otherwise the same. Such a choice would enforce detailed balance of the CTMC, with the desired steady-state distribution To implement such a CRN, we define a reaction set R that contains a distinct pair of reactions for each possible state of the neighbors of i for which w ij = 0. Let α i ∈ {ON, OF F } |N (i)| denote a state of neighboring species. Then, the necessary reactions and rate constants are for each i and every possible state α. In physical terms, the species representing the neighbors of node i collectively catalyze X OF F i X ON i , with a separate pair of reactions for every possible α i . While this entails a large number of reactions (2 |N (i)|+1 for each node i), it allows the rate constants for each configuration of neighbors to be distinct, and thus to satisfy the ratio of rate constants given in (19). For CRN states that satisfy the conservation laws x ON i + x OF F i = 1, there will be a unique reaction that can flip any given bit, and thus the CTMC detailed balance (17) also holds, yielding the correct π(s). The construction is illustrated by example in Figure 2 and compared to other constructions in Figure 3.
The distribution π(s) is identical to that of the BM, both with and without clamping. Reachability is preserved by clamping, as all states satisfying the conservation laws and clamping can be reached in the clamped CRN. All results derived for traditional BMs therefore apply, including conditional inference and the Hebbian learning rule. The construction can be generalized to any graphical model and indeed to any finite Markov chain defined on a positive integer lattice.
With the DCBM, we have shown that CRNs can express the same distributions as BMs, and are thus very expressive. However, since each possible state α i of N (i) is associated with two reactions, the number of reactions of the CRN is exponentially large in the typical node degree d in the original BM. Moreover, the scheme requires high molecularity reactions in which multiple catalysts effect a single transition (the molecularity grows linearly with d). Physical implementations are therefore likely to be challenging. We further note that as a consequence of Theorem 2, the DCBM cannot be detailed balanced at the level of the underlying chemistry, due to its simple conservation laws. Physically, this means that the DCBM must use a fuel species to drive each reaction. Details of this argument are given in the appendix (7.1).

The Edge Species CBM Construction (ECBM):
Can a detailed balanced CRN also implement a Boltzmann machine, or is it necessary to break detailed balance at the level of the CRN reactions, as in the DCBM? Here we show that it is not necessary by introducing a detailed balanced CRN that uses species to represent both the nodes and edges of a BM.
(20) This reaction scheme is visualized in Figure 3. Just like in the DCBM, there is a separate pair of reactions for each node i for each state of its neighbors α i . In this case, however, the backwards reaction in (20) does represent a transition that is a true chemical inversion of the forwards reaction. Further, given a valid initial state, clamping any subset of the X ON/OF F i species preserves reachability. The ECBM is detailed balanced and therefore its stationary distribution has the form (13). Substituting the conservation law w ij = x i x j and omitting species with 0 energy results in Comparing this expression to the distribution of a BM, equation (2), the above expressions are equivalent provided that their partition functions are equivalent. To see this is the case, notice that: 1) the partition function is just a sum over the Gibbs factors across the entire state space.
2) The Gibbs factors take the same form between the ECBM and BM (as shown above). And 3) the reachable state spaces spaces are equivalent. Thus a sum over all possible Gibbs factors will be equal. Therefore, Z BM = Z π and the theorem is proven.
Via the ECBM, we have shown that even detailed balanced CRNs can represent rich distributions and are able to calculate conditional distributions through clamping as proven in Theorem 1. Due to being detailed balanced, this construction requires no fuel molecules and performs sampling via the intrinsic equilibrium fluctuations of the CRN. Moreover, it is only necessary to tune molecular energies in this construction, since appropriate relative rate constants follow by definition. This construction is possible due to the complex set of conservation laws that ensure that the reachability classes of all the X ON/OF F i species are tightly coupled via the W ON/OF F ij species. One implication is that this construction does not generalize easily to non-binary species counts. Additionally, issues related to high molecularity reactions and large number of reactions remain.

Approximate Bimolecular Implementations
The DCBM and the ECBM both require reactions of high molecularity. High molecularity reactions and systems involving many species are physically challenging to implement and also potentially suffer from long mixing times. In this section, we discuss an approximation scheme to create CBMs with lower molecularity reactions and thus overcome these issues.

Taylor Series Chemical Boltzmann Machine (TCBM):
Here, we demonstrate a compact CBM that approximates a BM. It is not detailed balanced on either the Markov chain or the CRN level, but uses only 2N species and unimolecular and bimolecular reactions. The TCBM is a non-equilibrium CBM of the kind discussed in Section 3.1 that uses a dual-rail representation and single-node transitions to approximately implement a BM. The reactions are given by: which, with appropriate initial conditions, preserve the conservation law that This model's parameters can be taken directly from the weights of a BM, w ij . First, define a symmetric matrix W . Decompose this matrix into the difference of two positive matrices, W = A − B, where a ij ∈ A are all w ij > 0 and b ij ∈ B are the absolute values of all w ij < 0. Finally, k is an arbitrary overall rate. This construction can be understood as an approximation of equation (17), which dictates that for two states s and s that differ only in bit i with x ON i = 1 in state s , the CTMC transition rates must satisfy The bias θ i has been absorbed into w ij for notational clarity by assuming there is some x ON 0 = 1 whose weights act as biases. The TCBM is a bimolecular CRN obeying the same conservation laws as the DCBM in which each species j acts as an independent catalyst for transitions in i with reaction rates determined by a ij and b ij . The relative propensities of this network are exactly equal to the linear expansion of the relative propensities shown in the last term in (23). Specifically, the numerator is the sum of the reaction propensities for X OF F i → X ON i and the denominator is the sum of the reaction propensities for X ON i → X OF F i , in each case plus a constant term due to the unimolecular reactions. We thus propose the simple scheme in (22) as an approximate CBM; Figure 4A depicts this TCBM schematically. This model bears some resemblance to protein phosphorylation networks where adding or removing a phosphate group is analogous to turning a species on or off; both are driven, catalytic processes capable of diverse computation.

Approximate BCRN Inference:
Remarkably, this simple approximate CBM can reasonably approximate the inferential capabilities of a BM. We demonstrate this by using (22) to convert a BM trained on the MNIST dataset [47] to a TCBM ( Figure 4). We then compare the BM and the TCBM side by side. Digit classification is shown in Figure 4 panels C and D for a BM and in Figure 4 panels F and G for a CBM as confusion heatmaps. Classification is carried out by clamping the image nodes to MNIST images and averaging the values of the classification nodes. As is apparent from these plots, the BM does a fairly reasonable job classifying these digits, but struggles on the number 5. The CBM functions as a very noisy version of the BM with nodes in general much more likely to be on. The CBM has also faithfully inherited the capabilities and limitations of the BM and similarly struggles to classify the digit 5. Digit generation is shown in Figure 4E for a BM and 4H for a CBM. Generation was carried out by clamping a single class node to 1 and all other class nodes to 0, then averaging the output of the image nodes after the network had equilibrated. For each generated image, we show the raw output and the top 85th percentile of nodes, a thresholding which helps visualize the noisy output. As is apparent from the raw output, the CBM approximation scheme does not generate images nearly as distinctly as the BM. However, this approximation does faithfully reproduce plausible digits when filtered for the top 85th percentile. Additional training and simulation details can be found in the appendix (7.2) The overall performance of the CRN is reasonable, given the fact that weights were simply imported from a BM without re-optimization. The TCBM only approximates the distribution implied by these weights and, in the absence of detailed balance, does not have an established formal relationship between clamping and conditioning.

Detailed Balanced CRN Learning Rule
A broad class of detailed balanced chemical reaction networks can be trained with a Hebbian learning rule between a waking phase (clamped) and sleeping phase (free) that is reminiscent of the classic gradient descent learning algorithm for a BM [8,35]. A similar algorithm for detailed balanced CRNs follows.
First we state a simple case of Theorem 4 where we just want a CRN with stationary distribution π over Ω s0 to learn a target distribution Q also defined on Ω s0 . Then, the learning rule is given by Here, s i π and s i Q denote the expected count of the species S i with respect to the probability distributions π and Q, respectively, and g i = G[S i ] is the energy of species S i . Theorem 4 generalizes this procedure to cases where there are hidden species.
Theorem 4. Let C = (S, R, k) be a detailed balanced chemical reaction network with stationary distribution π(s) on Ω s0 . Consider a partition (V, H) of the set S of species into visible and hidden species such that π(s) = π(v, h). Require that for all visible states v, the clamped CRN C| V =v preserves reachability. Let be the weighted mixture of stationary distributions of the clamped CRNs C| V =v with v drawn from the distribution Q. Then, the gradient of the Kullback-Leibler divergence from π V to Q with respect to the energy, g i = G[S i ], of the species S i is given by Proof. Applying Theorem 1, the clamped CRN ensemble π Q (s) may be written Additionally we will need the partial derivative of a Gibbs factor and the partition function with respect to g i , Using these results, the partial derivative of any detailed balanced CRN's distribution at a particular state s, with respect to an energy g i , is ∂π(s) Noting that Q has no dependence on g i , the gradient of the Kullback-Leibler divergence can then be written, In the special case where there are no hidden species, which is to say the target distribution Q is defined over the whole reachability class Ω s0 , then π V (v) = π(s) and π Q (s) = Q(s) and the gradient has the simple form shown in equation (24). Applying gradient descent via dgi dt = − ∂D KL ∂gi , we thus have a simple in silico training algorithm to train any detailed balanced CRN such that it minimizes the Kullback-Leibler divergence from π V to Q. If H = ∅, simulate the CRN freely to estimate the average counts s i under π(s). Then compare to the average counts under the target Q(s) and update the species' energies accordingly. If H = ∅, clamp the visible species to some v ∈ Ω V s0 with probability Q(v) and simulate the dynamics of the hidden units. Repeat to sample an ensemble of clamped CRNs C| V =q . Because clamping v preserves reachability, Gillespie simulations of the CRN with the V species clamped to the data values v will sample appropriately. This gives the average counts under π Q .
This CBM learning rule is more general than the classical Boltzmann machine learning rule, as it applies to arbitrary detailed balanced CRNs, including those with arbitrary conservation laws and arbitrarily large species counts (but still subject to the constraint that reachability under clamping must be preserved). That said, at first glance the CBM learning rule appears weaker than the classical Boltzmann machine learning rule, as it depends exclusively on mean values s i , whereas the Boltzmann machine learning rule relies primarily on secondorder correlations x i x j . In fact, though, conservation laws within the CRN can effectively transform mean values into higher-order correlations. A case in point would be to apply the CBM learning rule to the ECBM network: For x ON j π , which exactly matches the classical Boltzmann machine learning rule if we assert that the energies of OF F species are fixed at zero.

Discussion
We have given one approximate and two exact constructions that allow CRNs to function as Boltzmann machines. BMs are a "gold standard" generative model capable of performing numerous computational tasks and approximating a wide range of distributions. Our constructions demonstrate that CRNs have the same computational power as a BM. In particular, CRNs can produce the same class of distributions and can compute conditional probabilities via the clamping process. Moreover, the TCBM construction appears similar in architecture to protein phosphorylation networks. Both models are non-equilibrium (i.e., require a fuel source) and make use of molecules that have an on/off (e.g., phosphorylated/unphosphorylated) state. Additionally, there are clear similarities between our exact schemes and combinatorial regulation of genetic networks by transcription factors. In this case, both models make use of combinatoric networks of detailed-balanced interactions (e.g., binding/unbinding) to catalyze a state change in a molecule (e.g., turning a gene on/off). We note that our constructions differ from some biological counterparts in requiring binary molecular  counts. However, in some cases we believe that biology may make use of conservation laws (such as having only a single copy of a gene) to allow for chemical networks networks to perform low-cost computations. In the future, we plan to examine these cases in a biological setting as well as generalize our models to higher counts. Developing these CBMs leads us to an important distinction between equilibrium, detailed-balanced CRNs with steady state distributions determined by molecular energies, and CRNs that do not obey detailed balance in the underlying chemistry. The second category includes those that nonetheless appear detailed balanced at the Markov chain level. Physically, this distinction is especially important: a non-detailed balanced CRN will always require some kind of implicit fuel molecule (maintained by a chemostat) to run and the steady state will not be an equilibrium steady state due to the continuous driving from the fuel molecules. A detailed balanced CRN (at the chemical level) requires no fuel molecules: and thus the chemical circuit can act as a sampler without fuel cost. Despite this advantage, working with detailed balanced CRNs presents additional challenges: to ensure that chemical species do not have independent distributions, species counts must be carefully coupled via conservation laws.
Our constructions also highlight important complexity issues underlying CBM design. The number of species, the number of reactions, and the reaction molecularity needed to implement a particular BM are relevant. Trade-offs appear to arise between these different factors, and the thermodynamic requirements of a given design. A breakdown of the main features of each CBM is given in Table 1. Summarizing, the TCBM is by far the simplest construction, using O(N ) species, at most O(N 2 ) reactions, with molecularity ≤ 2. However, this happens at the expense of not being an exact recreation of a BM, and the requirement of a continuous consumption of fuel molecules. The DCBM is the next simplest in complexity terms, using O(N ) species, O(N 2 N ) reactions, and molecularity of at most N . Like the TCBM, the DCBM requires fuel molecules because it is not detailed balanced at the CRN level. The ECBM is considerably more complex than the DCBM, using quadratically more species, O(N 2 ), the same number of reactions, O(N 2 N ) and double the reaction molecularity. The ECBM makes up for this increased complexity by being detailed balanced at the CRN level, meaning that it functions in equilibrium without implicit fuel species.
Finally, we have shown that a broad class of detailed balanced CRNs can be trained using a Hebbian learning rule between a waking phase (clamped) and sleeping phase (free) reminiscent of the gradient descent algorithm for a BM. This exciting finding allows for straightforward optimization of detailed balanced CRNs' distributions.
This work provides a foundation for future investigations of probabilistic molecular computation. In particular, how more general restrictions on reachability classes can generate other "interesting" distributions in detailed balanced CRNs remains an exciting question. We further wonder if the learning rule algorithm can be generalized to certain classes of non-detailed balanced CRNs. Another immediate question is whether these ideas can be generalized to nonbinary molecular counts. From a physical standpoint, plausible implementations of the clamping process and the energetic and thermodynamic constraints require investigation. Indeed, a more realistic understanding of how a CBM might be implemented physically will help us identify when these kinds of inferential computations are being performed in real biological systems and could lead to building a synthetic CBM.

Application of Theorem 2: The Direct CBM Must Use Implicit Fuel Species
Here, we use Theorem 2 to analyze the direct implementation of a CBM and show that it cannot be detailed balanced and thereby requires implicit fuel molecules. First, notice that the the conservation laws used in this construction are of a simple form. The states accessible by (X ON i , X OF F i ) are independent of (X ON j , X OF F j ) for i = j, and therefore the reachability class is a product over the subspaces of each individual node. As a consequence, by Theorem 2, the system must be out of equilibrium and violate detailed balance at the level of the CRN because, by construction, this system is equivalent to a BM and has correlations between nodes i and j whenever w ij = 0. In physical terms, the presence of catalysts cannot influence the equilibrium yield of a species, and therefore a circuit which uses catalysis to bias distributions of species must be powered by a supply of chemical fuel molecules [41,42,43]. It is also worth noting that, as a consequence, this scheme cannot be implemented by tuning of (free) energies; it is fundamentally necessary to carefully tune all of the rate constants individually (via implicit fuel molecules) to ensure that detailed balance is maintained at the level of the Markov chain for the species of interest.

BM Training and TCBM Simulation Details:
We trained a BM using stochastic gradient descent on the MNIST dataset, down sampled to be 10 pixels by 10 pixels [47]. The BM has 100 visible image units (representing a 10 x 10 image), 10 visible class nodes, and 40 hidden nodes as depicted in Figure 4B. Our training data consisted of the concatenation of down sampled MNIST images and their classes projected onto the 10 class nodes. The weights and biases of the trained BM were converted to reaction rates for a CBM using the Taylor series approximation. This CBM consists of 300 species, 300 unimolecular reactions and 22350 bimolecular reactions. The resulting CBM was then compared side-by-side with the trained BM on image classification and generation. The BM was simulated using custom Gibbs sampling written in Python. The CRN was simulated on a custom Stochastic Simulation Algorithm (SSA) [40] algorithm written in Cython. All simulations, including network training, were run locally on a notebook or on a single high performance Amazon Cloud server.
Classification was carried out on all 10000 MNIST validation images using both the BM and the CBM. Each 10 by 10 gray-scale image was converted to a binary sample image by comparing the gray-scale image's pixels (which are represented as real numbers between 0 and 1) to a uniform distribution over the same range. The network's image units were then clamped to the binary sample and the hidden units and class units were allowed to reach steady state. This process was carried out 3 times for each MNIST validation image, resulting in 30000 sample images being classified. Raw classification scores were computed by averaging the class nodes' outputs for 20000 simulation steps after 20000 steps of burn-in (Gibbs sampling for the BM, SSA for the CBM). Max classification was computed by taking the most probable class from the raw classification output. Raw classification and max classification confusion heatmaps, showing the average classification across all sample images as a function of the true label are shown in Figure 4 panels C and D for a BM and in Figure 4 panels F and G for a CBM.
Image generation was carried out by clamping the class nodes with a single class, 0...9, taking the value of 1 and all other classes being 0, and then allowing the network to reach steady state. Generated images were computed by averaging the image nodes over 50000 simulation steps (Gibbs sampling for the BM, SSA for the CBM) after 25000 steps of burn-in. Generation results are shown in Figure 4E for a BM and Figure 4H for a CBM.