Unified Conversational Recommendation Policy Learning via Graph-based Reinforcement Learning

Conversational recommender systems (CRS) enable the traditional recommender systems to explicitly acquire user preferences towards items and attributes through interactive conversations. Reinforcement learning (RL) is widely adopted to learn conversational recommendation policies to decide what attributes to ask, which items to recommend, and when to ask or recommend, at each conversation turn. However, existing methods mainly target at solving one or two of these three decision-making problems in CRS with separated conversation and recommendation components, which restrict the scalability and generality of CRS and fall short of preserving a stable training procedure. In the light of these challenges, we propose to formulate these three decision-making problems in CRS as a unified policy learning task. In order to systematically integrate conversation and recommendation components, we develop a dynamic weighted graph based RL method to learn a policy to select the action at each conversation turn, either asking an attribute or recommending items. Further, to deal with the sample efficiency issue, we propose two action selection strategies for reducing the candidate action space according to the preference and entropy information. Experimental results on two benchmark CRS datasets and a real-world E-Commerce application show that the proposed method not only significantly outperforms state-of-the-art methods but also enhances the scalability and stability of CRS.

INTRODUCTION
Conversational recommender systems (CRS) aim to learn user's preferences and make recommendations through interactive conversations [12,16,41]. Since it has the natural advantage of explicitly acquiring user's preferences and revealing the reasons behind recommendation, CRS has become one of the trending research topics for recommender systems and is gaining increasing attention. Unlike traditional recommender systems [8,23,43] or interactive recommender systems (IRS) [49,51], which mainly focus on solving the problem of which items to recommend, there exists generally the other two core research questions for CRS [11], namely what questions to ask and when to ask or recommend. Recent works have demonstrated the importance of interactivity of asking clarifying questions in CRS [4,41,50]. More importantly, deciding when to ask or recommend is the key to coordinating conversation and recommendation for developing an effective CRS [12,14,26].
Different problem settings of CRS have been proposed, either from the perspective of dialog systems, being a variation of taskoriented dialog [13,16,48], or from the perspective of recommender systems, being an enhanced interactive recommender system [3,12,26]. In this work, we study the multi-round conversational recommendation (MCR) setting [12,26], where the system asks questions about users' preferences on certain attributes or recommends items multiple times, and the goal is to make successful recommendation with the minimum number of interactions.
In MCR scenario, the CRS is typically formulated as a multi-step decision-making process and solved by reinforcement learning (RL) methods for policy learning [12,14,26]. As shown in Figure 1(a), RL-based IRS is only required to learn the policy to decide which items to recommend. However, the situation is more complicated in CRS, since there are two components that need to be coherently considered [19], i.e., conversation and recommendation components. Existing methods CRM [26] and EAR [12] employ policy gradient [27] to improve the strategies of when and what attributes to ask, while the recommendation decision is made by an external recommendation model. In order to reduce the action space in policy learning, another state-of-the-art method SCPR [14] only considers learning the policy of when to ask or recommend, while two isolated components are responsible for the decision of what to ask and which to recommend. The policy learning frameworks of these two kinds of CRS are presented in Figure 1(b) and 1(c).
Despite the effectiveness of these methods, there are some issues remained to be tackled for real-world applications: (i) The models trained with existing CRS methods lack generality to different Work done when Yang Deng was an intern at Alibaba. This work was supported by Alibaba Group through Alibaba Research Intern Program, and a grant from the Research Grant Council of the Hong Kong Special Administrative Region, China (Project Codes: 14200719). Figure 1: Illustration of policy learning frameworks for IRS and CRS, including CRM [26], EAR [12], and SCPR [14].
domains or applications, since there are three different decisionmaking processes to be considered in CRS, including what attributes to ask, which items to recommend, and when to ask or recommend. It requires extra efforts to train an offline recommendation model [12,14] or pretrain the policy network with synthetic dialogue history [12,26]. (ii) The policy learning is hard to converge, since the conversation and recommendation components are isolated and lack of mutual influence during the training procedure.
To address these issues, in this work, we formulate the aforementioned three separated decision-making processes in CRS as a unified policy learning problem to harness the ultimate goal of CRS as well as fill the gap between the recommendation and conversation components during the training procedure. Such unified conversational recommendation policy learning (UCRPL) aims at learning a unified policy to decide the action, either ask an attribute or recommend items, at each conversation turn to maximize the cumulative utility over the whole MCR process. The overview of the proposed unified policy learning for CRS is depicted in Figure 2.
However, the UCRPL problem comes along with two challenges: (i) How to systematically combine conversation and recommendation components for unified policy learning? (ii) How to deal with sample efficiency issues? As the action space becomes overwhelmingly large in UCRPL, including all available attributes and items, it requires tremendous amount of interaction data to learn the optimal policy. Fortunately, the graph structure, capturing the rich correlated information among different types of nodes (i.e., users, items, and attributes), enables us to discover collaborative user preferences towards attributes and items. Thus, we can leverage the graph structure to integrate the recommendation and conversation components as an organic whole, where the conversation session can be regarded as a sequence of nodes maintained in the graph to dynamically exploit the conversation history for predicting the action at the next turn. On the other hand, although the connectivity of the graph can also be used to eliminate invalid actions by path reasoning [14], there are still a large number of candidates left for action searching. Since users are not likely to be interested in all items and attributes, we can focus on the potentially important candidates to improve the sample efficiency of UCRPL.
To this end, we propose a novel and adaptive graph-based reinforcement learning framework, namely UNIfied COnversational RecommeNder (UNICORN). Specifically, due to the evolving nature of CRS, we leverage a dynamic weighted graph to model the changing interrelationships among users, items and attributes during the conversation, and consider a graph-based Markov Decision Process (MDP) environment for simultaneously handling the  decision-making of recommendation and conversation. Then we integrate graph-enhanced representation learning and sequential conversation modeling to capture dynamic user preferences towards items and attributes. In addition, two simple yet effective action selection strategies are designed to handle the sample efficiency issue. Instead of enumerating the whole candidate item and attribute set, we adopt preference-based item selection and weighted entropy-based attribute selection strategies to only consider potentially important actions.
In summary, the contributions of this paper are as follows: • We formulate three separated decision-making processes in conversational recommender systems, namely when to ask or recommend, what to ask, and which to recommend, as a unified conversational recommendation policy learning problem. • To tackle the challenges in UCRPL problem, we propose a novel and adaptive reinforcement learning framework, namely UNIfied COnversational RecommeNder (UNICORN), based on dynamic weighted graph. To handle the sample efficiency issue, we further design two simple yet effective action selection strategies. • Experimental results show that the proposed method significantly outperforms state-of-the-art CRS methods across four public benchmark datasets and a real-world E-Commerce application.

RELATED WORKS
Conversational Recommendation. Based on the problem settings, current CRS studies can be categorized into four directions [5,11]: (1) Exploration-Exploitation Trade-offs for Cold-start Users [4,17,40]. These approaches leverage bandit approaches to balance the exploration and exploitation trade-offs for cold-start users in conversational recommendation scenarios. (2) Question-driven Approaches [3,41,50] aim at asking questions to users to get more information about their preferences, which is often addressed as "asking clarification/clarifying question". (3) Dialogue Understanding and Generation [16,18,47,48]. These studies focus on how to understand users' preferences and intentions from their utterances and generate fluent responses so as to deliver natural and effective dialogue actions. (4) Multi-round Conversational Recommendation [12,14,26]. Under this problem setting, the system asks questions about the user's preferences or makes recommendations multiple times, with the goal of achieving engaging and successful recommendations with fewer turns of conversations. Among these problem settings, we focus on the MCR problem. RL in Recommendation. Reinforcement learning (RL) has been widely introduced into recommender systems due to its advantage of considering users' long-term feedbacks [44,45]. RL-based recommendation formulates the recommendation procedure as an MDP of the interactions between the user and a recommendation agent, and employs RL algorithms to learn the optimal recommendation strategies [22,25,44,45]. Recent works on sequential recommendation [30,37] and interactive recommendation [39,49,51] adopt RL to capture users' dynamic preferences for generating accurate recommendations over time. The goal of these approaches typically is to learn an effective policy for determining which items to recommend. As for CRS, RL-based methods are adopted to improve the strategies of the other two decision processes, including (i) what attributes to ask [12,26] and (ii) when to ask or recommend [14].
In order to simplify the overall framework of MCR with better scalability and generality, we formulate these three core decision processes in CRS as a unified policy learning problem.
Graph-based Recommendation. Graph-based recommendation studies mainly leverage the graph structure for two purposes. The first one is to enhance the recommendation performance by graphbased representation learning, including exploiting the structure information for collaborative filtering [7,31,46], and adopting knowledge graph embeddings as rich context information [9,38].
The other group of studies models recommendation as a path reasoning problem for building explainable recommender systems [20,32]. Recent years have witnessed many successful applications of graphbased RL methods on different scenarios of recommender systems [14,15,33,36,42,49]. For example, Xian et al. [36] employ a policy-guided graph search method to sample reasoning paths for recommendation, which is enhanced with adversarial actor-critic for demonstration-guided path reasoning [42]. Lei et al. [15] and Zhou et al. [49] employ graph convolutional network (GCN) [10] for state representation learning to enhance the performance of traditional RL methods on recommendation policy learning. In this work, we study graph-based RL method for the conversational recommendation scenario based on a dynamic weighted graph.

PROBLEM DEFINITION
Multi-round Conversational Recommendation. In this work, we focus on the multi-round conversational recommendation (MCR) scenario [12,14], the most realistic conversational recommendation setting proposed so far, in which the CRS is able to ask questions about attributes or make recommendations multiple times. Specifically, on the system side, the CRS maintains a large set of items V to be recommended, and each item is associated with a set of attributes P . In each episode, a conversation session is initialized by a user specifying an attribute 0 . Then, the CRS is free to ask the user's preference on an attribute selected from the candidate attribute set P cand or recommend a certain number of items (e.g., top-) from the candidate item set V cand . Following the assumptions from Lei et al. [12], the user preserves clear preferences towards all the attributes and items. Thus, the user will respond accordingly, either accepting or rejecting the asked attributes or the recommended items. The CRS updates the candidate attribute and item sets, and decides the next action based on the user response. The system-ask and user-respond process repeats until the CRS hits the target item or reaches the maximum number of turn . Unified Conversational Recommendation Policy Learning. Multi-round conversational recommendation aims to make successful recommendations within a multi-round conversation session with the user. At each timestep , according to the observation on past interactions, the CRS selects an action , either asking an attribute or recommending items. In return, the user expresses his/her feedback (accept or reject). This process repeats until the CRS hits the user-preferred items or reaches the maximum number of turn . Such MCR task can be formulated as a Markov Decision Process (MDP). The goal of the CRS is to learn a policy maximizing the expected cumulative rewards over the observed MCR episodes as where is the state representation learned from the system state and the conversation history, is the action that the agent takes at timestep , and (·) is the intermediate reward, abbreviated as .

METHODOLOGY
The overview of the proposed method, UNICORN, is depicted in Figure 3, which consists of four main components: Graph-based MDP Environment, Graph-enhanced State Representation Learning, Action Selection Strategy, and Deep Q-Learning Network.

Graph-based MDP Environment
The MDP environment is responsible for informing the agent about the current state and possible actions to take, and then rewards the agent based on how the current policy fits the observed user interactions. Formally, the MDP environment can be defined by a tuple (S, A, T , R), where S denotes the state space, A denotes the action space, T : S × A → S refers to the state transition function, and R : S × A → R is the reward function.

State.
As for the graph-based MDP environment, the state ∈ S at timestep is supposed to contain all the given information for conversational recommendation, including the previous conversation history and the full graph G that includes all the users, items, and attributes. Given a user , we consider two major elements: where rej ] denotes the conversation history until timestep , and G ( ) denotes the dynamic subgraph of G for the user at timestep (The graph construction will be introduced in Section 4.2). P denotes the user-preferred attribute. P rej and V rej are the attributes and items rejected by the user, respectively. The initial state 0 is initialized by the user-specified attribute 0 where V P ( ) is the set of item vertices directly connecting all P ( ) (i.e., items satisfying all the preferred attributes), and P V ( ) cand is the set of attribute vertices directly connecting to one of V ( ) cand (i.e., attributes belonging to at least one of the candidate items).

Transition.
We consider that the current state will transition to the next state +1 when the user responds to the action . In specific, if CRS asks an attribute and the user accepts it, the next state +1 will be updated by P ( +1) = P ( ) ∪ . Conversely, if the user rejects the action , +1 will be updated by P , G ( +1) ].

Reward.
Following previous MCR studies [12,14], our environment contains five kinds of rewards, namely, (1) rec_suc , a strongly positive reward when the user accepts the recommended items, (2) rec_fail , a negative reward when the user rejects the recommended items, (3) ask_suc , a slightly positive reward when the user accepts the asked attribute, (4) ask_fail , a negative reward when the user rejects the asked attribute, and (5) quit , a strongly negative reward when reaching the maximum number of turns.

Graph-enhanced State Representation
As we formulate conversational recommendation as a unified policy learning problem over a graph-based MDP environment, it is required to encode both the conversational and graph structural information into the latent distributed representations. In order to make use of the interrelationships among users, items, and attributes, we first adopt graph-based pre-training methods [2,35] to obtained node embeddings for all the nodes in the full graph G. Figure 3, we represent the current state of the graph-based MDP environment as a dynamic weighted graph. Formally, we denote an undirected weighted graph as G = (N, ), with the node ∈ N , the adjacency matrix element , denoting the weighted edges between nodes and . In our case, given the user , we denote the dynamic graph at timestep as G ( ) = (N ( ) , ( ) ):

Dynamic Weighted Graph Construction. As shown in
where ( ) is a scalar indicating the recommendation score of the item in the current state. In order to incorporate the user preference as well as the correlation between the asked attributes and the items, such weight ( ) is calculated as where (·) denotes the sigmoid function, , , and are the embeddings of the user, item, and attribute, respectively.

4.2.2
Graph-based Representation Learning. In order to comprehensively take advantage of the correlation information among the involved user, items, and attributes from the connectivity of the graph, we employ a graph convolutional network (GCN) [10] to refine the node representations with structural and relational knowledge. The representations of the node in the ( + 1)-th layer can be computed by: where N denotes the neighboring indices of the node , and are trainable parameters representing the transformation from neighboring nodes and the node itself, and is a normalization adjacent matrix as = − 1 2 − 1 2 with = , . 4.2.3 Sequential Representation Learning. Apart from the interrelationships among the involved user, items, and attributes, the CRS is also expected to model the conversation history in the current state. Unlike previous studies [12,14] that adopt heuristic features for conversation history modeling, we employ Transformer encoder [29] for capturing the sequential information of the conversation history as well as attending the important information for deciding the next action. As described in [29], each Transformer layer consists of three components: (i) The layer normalization is defined as LayerNorm(·). (ii) The multi-head attention is defined as MultiHead( , , ), where , , are query, key, and value, respectively. (iii) The feed-forward network with ReLU activation is defined as FFN(·). Take the -th layer for example: where ∈ R × denotes the embeddings, and is the sequence length. In our case, the input sequence (0) is the accepted attributes P ( ) in the current conversation history with the learned graph-based representation { ( ) : ∈ P ( ) }, where is the number of layers in GCN. After the sequential learning with Transformer layers, we can aggregate the information learned from both the graph and the conversation history by a mean pooling layer to obtain the state representation of : For simplicity, we denote the learned state representation of as ( ), where is the set of all network parameters for state representation learning, including GCN and Transformer layers.

Action Selection Strategy
A large action search space will harm the performance of the policy learning to a great extent [14]. Thus, it attaches great importance to handle the overwhelmingly large action space in UCRPL. To this end, we propose two simple strategies to improve the sample efficiency for candidate action selection. Preference-based Item Selection. In general, for candidate items to be recommended, we can consider only the action of making recommendations from a small number of candidate items that fit the user preference the most, since users are not likely to be interested in all items. To achieve this, we select top-candidate items from V ( ) cand into the candidate action space A at each timestep , which is ranked by the recommendation score ( ) in Eq. (6).
Weighted Entropy-based Attribute Selection. Whereas for candidate attributes to be asked, the expected one is supposed to be able to not only better eliminate the uncertainty of candidate items, but also encode the user preference. Inspired by [14], we adopt weighted entropy as the criteria to prune candidate attributes: where V denotes the items that have the attribute . Similar to item selection, we also select top-candidate attributes from P ( ) cand into A based on the weighted entropy score ( ) .

Deep Q-Learning Network
After obtaining the graph-enhanced state representation and the candidate action space, we introduce the deep Q-learning network (DQN) [21] to conduct the unified conversational recommendation policy learning. We further implement some techniques to enhance and stabilize the training of DQN. The training procedure of the Unified Conversational Recommender is presented in Algorithm 1.

4.4.1
Dueling Q-Network. Following the standard assumption that delayed rewards are discounted by a factor of per timestep, we define the Q-value ( , ) as the expected reward based on the state and the action . As shown in the rightmost part of Figure 3, the dueling Q-network [34] employs two deep neural networks to compute the value function (·) and advantage function (·), respectively. Then the Q-function can be calculated by: where (·) and (·) are two separate multi-layer perceptions with parameters and , respectively, and let = { , }. The optimal Q-function * ( , ), which has the maximum expected reward achievable by the optimal policy * , follows the Bellman equation [1] as: 4.4.2 Double Q-Learning with Prioritized Experience Replay. During each episode in the MCR process, at each timestep , the CRS agent obtains the current state representation ( ) via the graph-enhanced state representational learning described in Section 4.2. Then the agent selects an action from the candidate action space A , which is obtained via the action selection strategies described in Section 4.3. Here we incorporate -greedy method to balance the exploration and exploitation in action sampling (i.e., select a greedy action based on the max Q-value with probability 1 − , and a random action with probability ).
Then, the agent will receive the reward from the user's feedback. According to the feedback, the current state transitions to the next state +1 , and the candidate action space A +1 is updated accordingly. The experience ( , , , +1 , A +1 ) is then stored into the replay buffer D. To train DQN, we sample mini-batch of experiences from D, and minimize the following loss function: where is the target value based on the currently optimal * . To alleviate the overestimation bias problem in conventional DQN, we adopt Double Q-learning [28], which employs a target network ′ as a periodic copy from the online network. The target value of the online network is then changed to: where ′ denotes the parameter of the target network, which is updated by the soft assignment as: where is the update frequency.
In addition, the conventional DQN samples uniformly from the replay buffer. In order to sample more frequently those important transitions from which there is much to learn, we employ prioritized replay [24] as a proxy for learning potential, which samples transitions with probability relative to the absolute TD error:

Model Inference.
With the learned UNICORN model, given a user and his/her conversation history, we follow the same process to obtain the candidate action space and the current state representation, and then decide the next action according to max Q-value in Eq. (13). If the selected action points to an attribute, the system will ask the user's preference on the attribute. Otherwise, the system will recommend top-items with the highest Q-value to the user.

EXPERIMENT 5.1 Dataset
We evaluate the proposed method, UNICORN, on four existing multi-round conversational recommendation benchmark datasets, and further conduct evaluation on a real-world E-Commerce platform. The statistics of these datasets are presented in Table 1.
• LastFM and Yelp. The LastFM dataset is used for evaluation on music artist recommendation, while the Yelp dataset is for business recommendation. Lei et al. [12] manually categorize the original attributes in LastFM into 33 coarse-grained groups, and build a 2-layer taxonomy with 29 first-layer categories for Yelp. • LastFM* and Yelp*. Lei et al. [14] consider that it is not realistic to manually merge attributes for practical applications, so they adopt original attributes to reconstruct these two datasets. For a fair comparison, we adopt both versions for experiments. • E-Commerce. In order to evaluate the proposed method in realworld E-Commerce scenario, we collect the dataset by sampling logs from Taobao, the largest E-Commerce platform in China, in the same way as LastFM and Yelp. We randomly sample 30,000 items from the Women_Dress category and a certain number of users who have interacted with these items. And we treat the concatenation of the product property and the property value as the attribute to be asked, e.g., "Color=Red", where "Color" is a product property and "Red" is a property value. Following the MCR data construction described in Lei et al. [12], we obtain the E-Commerce dataset with the statistics shown in Table 1. As for graph construction, we consider two kinds of relations, including "purchased_by (item→ user)" and "belong_to (attribute→ item)". Compared with the imbalanced ratio of #items:#attributes at LastFM (0.88:1) and Yelp (119:1), the E-Commerce dataset preserves a large number of items as well as attributes, which is more in line with real-world applications in E-Commerce.

User Simulator.
Due to the interactive nature of MCR, it needs to be trained and evaluated by interacting with users. Following the user simulator adopted in [12,26], we simulate a conversation session for each observed user-item interaction pair ( , ).
We regard the item as the ground-truth target item and treat its attribute set P as the oracle set of attributes preferred by the user in this conversation session. The session is initialized by the simulated user who specifies a certain attribute randomly chosen from P . Then the session follows the process of "System Ask, User Respond" [41] as described in Section 4.

Baselines.
We compare the proposed method with several state-of-the-art methods on MCR as follows: • Max Entropy. This is a rule-based strategy to decide the next action, where the CRS selects an attribute to ask based on the maximum entropy within the current state, or chooses to recommend the top-ranked items with a certain probability [12]. • Abs Greedy [4]. This method only performs recommendation actions and updates itself until the CRS makes a successful recommendation or exceeds the maximum turns of conversation. • CRM [26]. This method is originally designed for single-round conversational recommendation, which employs policy gradient [27] to learn the policy deciding when and which attributes to ask. We follow [12] to adapt this method into MCR scenario. • EAR [12]. A three-stage method is proposed to enhance the interaction between the conversation and recommendation components with a similar RL framework as CRM. • SCPR [14]. This is the state-of-the-art method on MCR setting, which leverages interactive path reasoning on the graph to prune off candidate attributes and adopts the DQN [21] framework to determine when to ask or recommend.

Evaluation Metrics.
Following previous studies on multiround conversational recommendation [12,14], we adopt success rate at the turn (SR@ ) [26] to measure the cumulative ratio of successful conversational recommendation by the turn . Besides, average turn (AT) is adopted to evaluate the average number of turns for all sessions (if the conversation session reaches the maximum turn , the turn for such session is also counted as ). The higher SR@ indicates a better performance of the CRS at a turn , while the lower AT means an overall higher efficiency.
Although SR@ and AT are widely adopted for the evaluation of CRS, both of them are not sensitive to the rank order in the recommendation results and only consider a certain aspect of the performance of CRS. In order to conduct a comprehensive evaluation of CRS, we propose to extend the normalized discounted cumulative gain (NDCG@ ) to be a two-level hierarchical version, namely hNDCG@( , ), with the definition as follows: where and represent the number of conversation turns and recommended items in each turn, ( , ) denotes the relevance of the result at the turn and position . Then hNDCG@( , ) can be drawn from hDCG@( , ) with the same way as the original NDCG@ . Since there is only one target item for each session in MCR, we simply use hDCG@( , ) for evaluation. The intuition behind hNDCG@( , ) is that the less number of turns of a successful session is favorable for the CRS, while the target item is expected to be ranked higher in the recommendation list at the success turn.

Implementation Details.
Following [14], we split the E-Commerce dataset by 7:1.5:1.5 for training, validation, and testing, and set the size of the recommendation list as 10, the maximum turn as 15. We adopt TransE [2] from OpenKE [6] to pretrain the node embeddings in the constructed graph with the training set. We adopt the user simulator described in Section 5.2.1 to interact with the CRS for online training the model using the validation set. For all implemented methods, we conduct the online training for 10,000 episodes. To maintain a fair comparison with other baselines [12,14], we adopt the same reward settings to train the proposed model:  GCN layers and Transformer layers are set to be 2 and 1, respectively. The numbers of selected candidate attributes and items are set to be 10. During the training procedure of DQN, the size of experience replay buffer is 50,000, and the size of mini-batch is 128. The learning rate and the 2 norm regularization are set to be 1e-4 and 1e-6, with Adam optimizer. The discount factor and the update frequency are set to be 0.999 and 0.01. Table 2 shows the performance comparison between the proposed method, UNICORN, and all baselines across five datasets. In general, UNICORN outperforms all the baselines by achieving a significantly higher success rate and less average turn, which is also comprehensively validated by the improvements on hDCG. As for the real-world E-Commerce dataset, SCPR outperforms EAR and CRM, whose performance is largely affected by the large action space in E-Commerce dataset. Despite the effectiveness of SCPR to handle the issue of large action space in CRS, its performance in such a real-world application is still restricted by the separation of its conversation and recommendation components. Our UNICORN not only enables the conversation and recommendation to be mutually enhanced during the training process, but also attains an effective sample efficiency with the proposed action selection strategies. This leads to a substantial margin  from these baselines, about 18% for SR@15, 2 turns for AT, and 30% for hDCG. Detailed analyses can be found as follows.

5.3.2
Training Efficiency. Figure 4 shows the test performance curves of different methods. Since Max Entropy and Abs Greedy are unsupervised methods, their curves are presented as two lines for comparisons. It can be clearly observed that UNICORN is trained more stably and requires fewer training episodes (i.e., interaction data) to achieve a better performance than other strong baselines. Among these baselines, the curve of SCPR is the most vibrant, since it only considers the policy of when to ask or recommend, while the decisions of question-asking and recommendation are made by two separated components. As for EAR and CRM, due to the large action space in the last three datasets, there is no much performance increase during the online training process, even getting worse. These results demonstrate the efficiency and effectiveness of the proposed unified policy learning for CRS.

Comparison at Different Conversation Turns.
Besides SR@15, we also present the performance comparison of success rate at each turn (SR@ ) in Figure 5. In order to better observe the differences among different methods, we report the relative success rate compared with the state-of-the-art baseline SCPR. For example, the line of = 0 represents the curve of Success Rate* for SCPR against itself. There are several notable observations as follows: (i) The proposed UNICORN substantially and consistently outperforms these baselines across all the datasets and almost each turn in the conversation session.
(ii) Due to the nature of greedy recommendation approach (Abs Greedy), it may successfully hit target items at the early stage of the conversation, leading to a relatively strong performance at the first few turns, but the performance falls quickly as the turn increases.
(iii) UNICORN achieves an outstanding performance in the middle stage of the conversation, where there are still a large number of candidate items and attributes to be pruned. This shows the strong scalability of UNICORN to effectively handle large candidate action space in different situations.
(iv) The performance of SCPR is getting closer to UNICORN at the latter stage of the conversation, as the candidate item and attribute set is getting smaller and the task becomes easier.
(v) EAR and CRM share similar performance as Abs Greedy in those datasets with a large candidate attribute set, i.e., Yelp* and E-Commerce, indicating their policy learning is merely working when encountering a large action space.

Model Components.
We first evaluate the model components for the state representation learning, including the graphbased representation learning and the sequential representation learning. The ablation studies are presented in the first part in Table 3 (row (a-b)). We observe that the performance on 4 out of 5 datasets suffers a larger decrease by discarding sequential learning (row (a)) than discarding graph-based learning (row (b)). We attribute this to two reasons: (i) The pretrained TransE embeddings have already encoded certain graph-based knowledge into the node representations. (ii) Since there are relatively larger number of attributes in these four datasets than LastFM, it is more important to model the interrelationships among diverse user-accepted attributes for capturing user preferences.

5.4.2
Pre-training Embeddings. The second part in Table 3 (row (c-d)) presents the results that we replace pretrained TransE embeddings with randomly initialized embeddings (row (c)) and pretrained embeddings from FM model described in Lei et al. [12] (row (d)). We observe a substantial performance gap between UNI-CORN with randomly initialized embeddings and that with pretrained embeddings, either TransE or FM. This result shows that  pretraining is necessary and simply online training with interactive conversations is far from efficient and sufficient for training a CRS with limited interaction data. In addition, using TransE embeddings further improves the performance from FM embeddings to a great extent, which also validates the advantages of graph-enhanced RL. Table 3 (row (e-h)) presents the results that we replace or discard the proposed action selection strategies. One alternative attribute selection strategy is to adopt the same strategy as the preference-based item selection by changing the object from items to attributes. Another one is to use the original maximum entropy function [12]. The results (row (e,f)) show that the performance of UCRPL suffers a noticeable decrease when adopting the preference-based or entropy-based strategy, indicating that it is required to consider both the user preference and the capability of reducing candidate uncertainty when deciding the asked attribute. Without attribute selection, we observe that the impact on applications with small action space (e.g., LastFM and Yelp) is less than those with large action space (e.g., LastFM*, Yelp*, and E-Commerce). UNICORN is merely working without item selection, since there are no pretrained recommendation components in the framework and the preference-based item selection serves as an auxiliary item recall process.

Parameter Sensitivity Analysis
The upper part in Table 4 summarizes the experimental results (hDCG) by varying the number of selected candidate actions. Since similar conclusions can also be drawn on other datasets, we only report the results on LastFM* and E-Commerce datasets due to space limit. As for the number of selected attributes, it is likely to discard the important attributes when only selected the attribute with the highest weighted entropy (e.g., =1). However, within a certain training interaction period (10,000 episodes in our case), UNICORN generally achieves the best performance when only selecting a small number of candidate items for policy learning.
The results also demonstrate the necessity of pruning the available actions when there is a large action search space in UCRPL. The lower part shows the experimental results by varying the number of network layers. As for GCN, 1-hop aggregation ( =1) is not enough for capturing all the useful information for graph representational learning. However, there is no much difference between 2-hop ( =2) and 3-hop ( =3). Besides, for Transformer, we can see that 1-layer ( =1) is sufficient for a good performance.

Qualitative Analysis
In order to intuitively study the difference between the proposed UNICORN and other state-of-the-art CRS methods, we randomly sample a real-world interaction from the E-Commerce dataset. The generated conversations by UNICORN, SCPR, EAR, and CRM with the user simulator are presented in Figure 6. Facing the large candidate action space, CRM tends to only trigger the recommendation component to make recommendations, and EAR continuously asks the questions that are not preferred by the user. Despite the success of SCPR in predicting user-preferred attributes, the policy learning in SCPR only decides when to ask or recommend, based on the number of candidate items, which leads to some unnecessary or redundant question-asking turns. UNICORN systematically addresses these issues by making a comprehensive decision of the next action.

CONCLUSIONS
In this work, we formulate three separated decision-making processes in CRS, including when to ask or recommend, what to ask and which to recommend, as a unified policy learning problem. To tackle the unified conversational recommendation policy learning problem, we propose a novel and adaptive RL framework, which is based on a dynamic weighted graph. In addition, we further design two simple yet effective action selection strategies to handle the sample efficiency issue. Experimental results show that the proposed method significantly outperforms state-of-the-art CRS methods across four benchmark datasets and the real-world E-Commerce application with remarkable scalability and stability.


Learning model discrepancy: A Gaussian process and sampling-based approach Mechanical Systems and Signal Processing

Predicting events in the real world with a computer model ( simulator ) is challenging. Every simulator, to varying extents, has model discrepancy , a mismatch between real world observations and the simulator (given the ‘true’ parameters are known). Model discrepancy occurs for various reasons, including simpliﬁed or missing physics in the simulator, numerical approximations that are required to compute the simulator outputs, and the fact that assumptions in the simulator are not generally applicable to all real world contexts. The existence of model discrepancy is problematic for the engineer as performing calibration of the simulator will lead to biased parameter estimates, and the resulting simulator is unlikely to accurately predict (or even be valid for) various contexts of interest. This paper proposes an approach for inferring model discrepancy that overcomes non-identiﬁability problems associated with jointly inferring the simulator parameters along with the model discrepancy. Instead, the proposed procedure seeks to identify model discrepancy given some parameter distribution, which could come from a ‘likelihood-free’ approach that con- siders the presence of model discrepancy during calibration, such as Bayesian history matching. In this case, model discrepancy is inferred whilst marginalising out the uncertain simulator outputs via a sampling-based approach, therefore better reﬂecting the ‘true’ uncertainty associated with the model discrepancy. Veriﬁcation of the approach is per- formed before a demonstration on an experiential case study, comprising a representative ﬁve storey building structure. (cid:1) 2020 The Authors. Published by Elsevier Ltd. ThisisanopenaccessarticleundertheCCBY

where z 2 R NÂ1 are real world observations for a set of inputs X 2 R NÂdx . The simulator g Á; Á ðÞ depends on a set of inputs X and parameters h 2 R MÂ1 , whereas the model discrepancy term d Á ðÞ is assumed to only be dependent on the inputs X. Finally, e 2 R NÂ1 is assumed Gaussian additive noise. The decision about what the inputs X and parameters h of a simulator are will be application specific, where generally, the inputs X are chosen to reflect the measured outputs z and the parameters h are additional variables in the simulator that could be tuned. For example, the inputs X could be variables that drive a physical process, such as a force or crack length, or independent variables like spectral lines or frequency bins, with the parameters h being, for example, material properties, or even mass, stiffness and damping coefficients, such as in direct model updating [18].
The main difference between the proposed method and existing techniques, such as in [1], is that the parameters h are identified before the functional model discrepancy term d Á ðÞ is inferred; rather than jointly inferring both h and d Á ðÞ. This decoupling assumption is made as the alternative joint inference problem is susceptible to non-identifiability issues, caused by modelling model discrepancy d Á ðÞwith a Gaussian process whilst inferring the parameters [4,6]. In fact, this modelling choice in the joint approach, leads to a rather poor likelihood, where even 'bad' parameter samples are given a high probability in the likelihood function [6], due to the ability of the GP to model any arbitrary function well [9], making the likelihood insensitive. Instead, by decoupling the problem, inferring the parameters and then the model discrepancy, nonidentifiability issues and problems with the likelihood can be overcome.
In order to decouple the parameter and model discrepancy inferences, the calibration method must be able to account for model discrepancy in another way. 'Likelihood-free' approaches, such as Bayesian history matching, offer such a technique. These methods incorporate model discrepancy though a notion of distance, removing issues associated with defining a specific likelihood, whilst approximating the parameter posterior distribution p hjz ðÞ . Once obtained, model discrepancy can be inferred using a Gaussian process model -without affecting the parameter posterior distribution. However, the Gaussian process model must be constructed from the uncertain simulator outputs p yjX; h ðÞ to the real world observations z. Unfortunately, it is not possible in closed-form to create a Gaussian process from uncertain inputs, meaning that a sampling-based solution is required.
It is noted that in the joint inference method, the parameters h are inferred, given an empirical Bayes estimate of the model discrepancy hyperparameters, requiring the Gaussian process to be conditioned on the parameter prior distributions, e.g.
R p zjy; /; h ðÞ p h ðÞ dh [1][2][3][4][5][6][7] -this can bias the inferred model discrepancy. This conditioning, typically leads to a restricted choice in prior distributions p h ðÞthat are conjugate with a Gaussian process, such as a Gaussian [1] or uniform [4] distribution; where non-conjugate priors require an additional expensive sampling procedure, on top of the parameter estimation, which is performed in low dimensions by a quadrature approach, and in high dimensions by another sampling step [1]. These issues are removed by considering the decoupled approach proposed in this paper.
By decoupling the inference procedure, the model discrepancy method can also be applied in scenarios where the parameter distribution is obtained by some elicitation process or experimentation. This makes the technique more generally applicable to a wider range of problems outside of those originally considered by the joint approach.
The proposed method assumes that some parameter distribution p hjz ðÞ has been obtained from a 'likelihood-free' calibration method, or p h ðÞhas been acquired from some elicitation process; where, for simplicity, p h ðÞwill be used to denote a generic parameter distribution. The approach then seeks to find the additive model discrepancy (and noise) term modelled using Gaussian process (GP) regression. As the simulator outputs, given the parameter distribution, are uncertain, a sampling-based approach is used to marginalise out the simulator outputs, meaning (potentially calibrated and) biascorrected model predictions can be made, reflecting the uncertainty from the parameter distribution. A brief outline of the approach is as follows: Obtain N s samples from the parameter distribution i.e. for the jth sample, h j ðÞ $ p h ðÞ . Propagate those N s samples through the simulator to obtain N s simulator output (denoted y 2 R NÂ1 ) samples i.e. y j Learn a GP mapping for each of the N s output samples y j ðÞ to a set of training observations z i.e. GP j ðÞ : y j ðÞ ; X ÈÉ ! z and obtain a weight w j for each regression model -the weights will be formally introduced in Section 2.2.
Calculate the weighted average of the set of GP regression models generating a bias-corrected model prediction It is noted that in the case where the simulator is computationally expensive to evaluate, a more computationally efficient emulator, or surrogate model [19], can be constructed. This efficient approximation can be sampled instead of the simulator in step two, where any emulator technique within the literature could be implemented [19][20][21]; in this paper a Gaussian process emulator is utilised.

Gaussian process regression
Model discrepancy is modelled in this paper by GP regression as it is a flexible, nonparametric tool, and because it has a Bayesian formulation allowing the uncertainty associated with the inferred functional form to be estimated [22,9]. These properties, the ability to approximate any unknown function well whilst quantifying the uncertainty in the prediction, are useful as the functional form of the model discrepancy is unknown a priori and quantifying the uncertainty associated with this form may aid simulator developers in targeting improvements to their computer model. In addition, by decoupling the inference problem, the choice of modelling model discrepancy with a Gaussian process no longer affects the likelihood in the parameter inference stage, making it a more suitable assumption.
The model discrepancy term is assumed in Eq. (1) to be additive, meaning it can be formed as a map from the simulator outputs y and inputs X to the observational data z; GP : y; X fg ! z, where the inferred model discrepancy GP can be related back to the inputs as d X ðÞ . For this reason GP regression is introduced in this section with the simulator outputs y being part of the inputs to the GP along with X, and where the noisy observations z are the outputs of the GP.
A Gaussian process states a prior distribution over a latent function f y; X ðÞ (of the noisy function z y; X where GP Á; Á ðÞ is a Gaussian process, with a mean function m Á ðÞand covariance function k Á; Á ðÞ which define the prior belief about the types of possible functions that could model the function f 2 R NÂ1 . Here a zero mean function is assumed, i.e. m y; X ðÞ ¼ 0, although it is trivial to add a non-zero mean function. The covariance function defines the correlation between any two points in the input space (hence being a function of y; X ðÞ and y 0 ; X 0 ÀÁ ) in a Reproducing Kernel Hilbert Space (RKHS) and is fully specified by a set of hyperparameters /, i.e. K ¼ k Á; Á; / ðÞ . The covariance function utilised in this paper is a Matérn 3/2 covariance function, as it is ideal for modelling relatively 'smooth' real world functions, 1 as it is (3/2-1) times mean square differentiable [23], and is defined as, where, and, where K f ;f 2 R NÂN is the covariance matrix for inputs 2 y 2 R NÂ1 and X 2 R NÂdx , r 2 f is the signal variance hyperparameter, and L y ¼ diag l y1 ; ...; l yd ÀÁ and L x ¼ diag l x1 ; ...; l xd ðÞ are lengthscale hyperparameters (making the covariance function an automatic relevance determination prior, i.e. it reduces the effect of redundant inputs). The covariance structure here separates out y and X allowing them to have an independent relationship with the outputs z. The hyperparameter vector for the Matérn 3/2 covariance function is therefore / ¼ r f ; L y ; L x

ÈÉ
. It is noted that the notation K f ;Ã ¼ k y; X ðÞ ; y Ã ; X Ã ðÞ ðÞ is used, where f indicates training and Ã test data. In order to make predictions, the joint Gaussian distribution is formed between a set of training D ¼ y; X fg ; z fg and testing data y Ã ; X Ã fg ; z Ã fg , assuming a Gaussian likelihood, where I f 2 R NÂN and I Ã 2 R NÃÂNÃ are identity matrices and r 2 n is a Gaussian noise variance. Following standard Gaussian conditionals, the posterior for a Gaussian process regression model can be formed as, Conventionally, GP models are inferred by taking a type-II maximum likelihood approach [9], i.e. finding the hyperparameters that maximise the marginal likelihood, leading to an empirical Bayes estimate of the hyperparameters/. By combining the noise variance with the set of covariance function hyperparameters, i.e. / ¼ r f ; L y ; L x ; r 2 n È É , the empirical Bayes estimates for the set of hyperparameters may be found through optimisation, -here a global optimisation approach is used, specifically via quantum particle swarm [24] -by minimising the negative log marginal likelihood, where, À log p zjy; X; / ðÞ ¼ It is noted that a fully Bayesian analysis would require marginalisation of the hyperparameters, which is not solvable in closed-form due to the dependence of the hyperparameters in the covariance function. However, the fully Bayesian solution may be inferred from a sampling-based approach [25][26][27] and is explored in the following section.

Sampling-based approach
The method outlined in this paper utilises GP regression in identifying the map GP : y; X fg ! z, and therefore inferring the model discrepancy term d X ðÞ . However, the output from a simulator y will typically be uncertain, i.e. p yjX; h ðÞ , arising from parametric uncertainty in p h ðÞ . However, Gaussian process regression cannot be solve in closed-form for uncertain inputs, and even though the simulator inputs X are deterministic, the simulator outputs y are uncertain. To create biascorrected predictions that account for this parametric uncertainty, the simulator outputs y Ã must be integrated out, forming the following integral, where p z Ã jX Ã ; h; D; / ð Þ is the bias-corrected predictive output, p y Ã jX Ã ; h ðÞ is the simulator prediction at test inputs X Ã which is conditioned on the parameter distribution p h ðÞ . It is noted that in previous work [8] the bias-corrected outputs have been approximated using the maximum a posteriori (MAP) estimate of the parameters (and an empirical Bayes estimate of the hyperparameters/) meaning p z Ã jX Ã ; h MAP ; D;/ . However, this will not account for the complete parametric uncertainty from p h ðÞand may result in a biased estimate of the model discrepancy. The proposed method seeks to approximate Eq. (10) via a sampling-based approach -specifically from an importance sampling viewpoint. Importance sampling is a technique for obtaining unbiased estimates of expectation integrals [28], such as in Eq. (10), and can be generalised as, where fx ðÞis a function, px ðÞthe nominal distribution over the variable x and qx ðÞis the proposal distribution, where X $ q are independent draws from the proposal distribution. The expectation in Eq. (11) can be formed as, where N are the number of samples and wX ðÞ ¼ pX ðÞ =qX ðÞ are the importance weights [28]. Eq. (10) can be approximated by setting the nominal and proposal distributions equal to the simulator output predictive distribution p y Ã jX Ã ; h ðÞ . This means that N s samples can be obtained from the parameter distribution, i.e. h j ðÞ $ p h ðÞand propagated through the simulator to obtain output samples y j , meaning that the weight for each sample equals one, i.e. w y Ã j ðÞ ÀÁ ¼ 1 (this effectively is the same as approximating the integral via Monte Carlo sampling, however the language of importance sampling will useful later in this section). The predictive equation p z Ã jX Ã ; h; D; / ð Þ can now be approximated using the set of weights and Gaussian process predictions for each sample, and the laws of total expectation and variance, where E z Ã ðÞ and V z Ã ðÞ are the GP predictive mean and covariance from Eq. (7). The bias-corrected predictions are approximately Gaussian, given that they are formed from weighted averaged Gaussian processes. The method is outlined in Algo-rithm 1. The main computational expense is in sampling the simulator outputs, which can be reduced by replacing the simulator with a computationally efficient emulator [19][20][21]. One problem with this approach is that the predictions are still dependent on a set of hyperparameters (that have been inferred from the GP associated with the parameter MAP estimates). However, these hyperparameters can also be marginalised out of the predictive equations using importance sampling, as discussed below.

Algorithm 1. Model discrepancy inference dependent on empirical Bayes estimates of /
The following integrals can be solved to generate a bias-corrected prediction not dependent on either the simulator outputs or GP hyperparameters, requiring the posterior of the hyperparameters p /jD ðÞ , which can also be approximated via importance sampling. By solving Eq. (14), rather than optimising the GP via a type-II maximum likelihood technique, the fully Bayesian solution of the hyperparameters can be acquired.
The posterior distribution of the hyperparameters p /jD ðÞ / p zjy; X; / ðÞ p / ðÞ (where D¼ y; X fg ; z fg ), can be approximated with importance sampling by setting the unnormalised nominal distribution as p zjy; X; / ðÞ p / ðÞ . By keeping the proposal and nominal distributions for the parameters the same as in the first approach, and by setting the proposal distribution for the hyperparameters equal to the prior distribution, i.e. / k ðÞ $ p / ðÞ ¼ q / ðÞ , the weights for each of the N / hyperparameter samples are equal to p zjy; X; / k ðÞ , i.e. the marginal likelihood of the GP model (in negative log-form in Eq. (9)). This formulation now allows both y Ã and / to be integrated out with importance sampling forming p z Ã jX Ã ; h; D ð Þ as, where the overall procedure is outlined in Algorithm 2. Is it noted that the approaches rely on importance sampling, which will suffer from the curse of dimensionality as the dimension of Y and / increases, as it will be less likely that a given sample will carry some meaningful weight. This issues can be mitigated by an adaptive approach [26], which is left for further research.

Obtaining the parameter distribution
As aforementioned, the proposed model discrepancy inference approach is primarily designed to accompany a 'likelihood-free' calibration process, where model discrepancy is account though a distance measure. These approaches allow for p h ðÞto be obtained without the GP model discrepancy influencing the posterior parameter distribution, removing problems associated with non-identifiability. This section briefly introduces Bayesian history matching as one such approach that can be used in combination with the proposed procedure (for more details on Bayesian history matching the interested reader is referred to [8]).

Bayesian history matching
Bayesian history matching (BHM) is an approximate Bayesian approach for calibrating statistical models of the form in Eq.
(1). The method seeks to determine whether parameter combinations are 'implausible', i.e. they are not likely to have produced the observations z, based on a criteria such that the remaining non-implausible parameter space is identified; leading to an approximation of the posterior distribution p hjz ðÞ . The criteria for discarding implausible samples is a combination of an implausibility metric and a threshold T, where the implausibility metric accounts for model discrepancy though a notion of distance.
BHM assumes that the simulator is computationally expensive to evaluate, and hence replaces the simulator with a computationally efficient GP emulator g X; h ðÞ $ G P : X; h fg ! y. This replacement of the simulator for an emulator is possible as a GP estimates the uncertainty associated with the approximation through the predictive variance V yX ; h ðÞ ðÞ from  (7)). This means that parameter combinations are only discarded when the approximation is certain enough, given the other uncertainty in the implausibility metric. The implausibility metric assess the distance between the mean emulator prediction E g Ã ðÞ and the observed data z, weighted by several uncertainties, The threshold is set given the statistical properties of the implausibility metric [8], e.g. for Eq. (16) the threshold can be set as T ¼ 3 given Pukelsheim's 3r rule [29]. By updating the emulator approximation at each iteration, the criteria can discard more parameter combinations with increased confidence. Finally, convergence is reached when either the uncertainty in the emulator is lower than the remaining uncertainties i.e. V o þ V m > V yX ; h ðÞ ðÞ or all the parameters are discarded.

Case study: numerical verification problems
The proposed two stage calibration and model discrepancy inference was verified on two numerical case studies; one where x 2 R NÂ1 and the other x 2 R NÂ2 . In addition, the first numerical case study is used to benchmark the proposed decoupled approach, using Bayesian history matching and the sampling-based model discrepancy procedure, against the hierarchical Bayesian model formulation, where the model discrepancy and model parameters are jointly inferred together. In both case studies the simulator modelled the tip deflection of a cantilever beam subject to an open crack with a point force of 10kN at the tip. The stiffness reduction model for an open crack used in this case study was that proposed by Christides and Barr [30], where the stiffness along the beam EI Á ðÞis a function of the length along the beam x, Young's modulus E, the second moment of area for the undamaged beam I 0 , the beam thickness t, the crack location l loc and a, a coefficient experimentally defined by Christides and Barr as 0:667. The constant C ¼ I 0 À I c ðÞ =I c is a function of the undamaged I 0 and damaged second moments of area I c , which for a rectangular beam are I 0 ¼ wt 3 ÀÁ =12 and I c ¼ wtÀ l cr ðÞ 3 =12; where w is the beam width and l cr is the crack length. The tip deflection was numerically estimated via Euler-Bernoulli bending beam equation, where M Á ðÞis the moment along the beam. In both case studies the beam used in the analysis was rectangular with the following dimensions: l ¼ 1m, w ¼ 0:5m and t ¼ 0:1m.

Numerical case study: one input problem
The first illustrative case study considers a scenario where the input was crack location x ¼ l loc , the parameters were Young's modulus E and crack length l cr , i.e. h ¼ l cr ; E fg , and the output (both from the simulator and experiments) was the tip deflection. In this analysis the true parameters were defined asl cr ¼ 38mm and b E ¼ 68GPa, i.e.ĥ ¼ 38; 68 fg ; the simulator evaluated at these parameters is depicted in Fig. 1. The training inputs for both the simulator and experimental data were 13 equally spaced points from 0:1m to 0:9m, and the simulator parameters were evaluated between 0mm and 50mm, and 50GPa and 90GPa, resulting in 25 equally spaced data points. These training inputs were used to construct an emulator with a linear mean and Matérn 3/2 covariance. The model discrepancy was defined as, d x ðÞ¼0:31 :5 À x ðÞ Â sin 1:8 x À 0:2 ðÞ Â 2p ðÞ ð 19Þ shown in Fig. 1. The experimental data z was formed from the simulator output plus the additive model discrepancy where the observation noise was Gaussian distributed with variance 0:001, where the experimental data is displayed in Fig. 1. The prior model discrepancy uncertainty (used in BHM) was V m ¼ 0:05, reflecting the expected magnitude of the model discrepancy, where the error bars on the experimental data in Fig. 1 show the total prior uncertainties. BHM was used to find the approximate posterior distribution, shown in Fig. 2 where the true parameters (shown in red) are close to the mode of the joint posterior distribution. Samples from the posterior distribution are shown in Fig. 1, where the output from the mode of the posterior distribution is visually in good agreement with the output at the true parameter values.  The proposed model discrepancy inference procedure was subsequently run for three scenarios: 1. Using a MAP estimate of the simulator parameters h MAP and an empirical Bayes estimate of the GP hyperparameters/ (when the simulator output for h MAP is used). 2. Marginalising out the simulator outputs y Ã using importance sampling, with an empirical Bayes estimate of the GP hyperparameters/ (when the simulator output for h MAP is used); N s ¼ 1000. 3. Marginalising out both the simulator outputs y Ã and GP hyperparameters / via importance sampling; N s ¼ 1000 and For each scenario a zero mean and Matérn 3/2 covariance function were used. The hierarchical Bayesian approach was also applied to the same training dataset where the prior parameter distributions were l cr $N 35; 10 ðÞ and E $N 70; 36 ðÞ and the model discrepancy Gaussian process was also modelled with a Matérn 3/2 covariance function. These prior parameter distributions are more informative than those consider in BHM (which can be understood as a uniform prior), where a comparison is shown in Fig. 3. The reason for using more informative priors in the hierarchical Bayesian approach than in BHM, is due to the findings of Brynjarsdóttir and O'Hagan [6], where to reduce the problem of non-identifiability, one solution is to apply more informative priors, constraining the posterior when the likelihood is relatively flat as a result of the model discrepancy Gaussian process; although it is noted that obtaining more informative priors is challenging in practical applications. Other aspects of the hierarchical Bayesian analysis were kept the same as in the BHM approach, such that objective comparisons could be made. The posterior distributions from the hierarchical Bayesian approach were obtained via an adaptive Markov chain Monte Carlo scheme [31,32], such that 100; 000 posterior samples were obtained with a 50; 000 sample burn-in period. The autocorrelations of the chains were checked for stationarity in order to confirm convergence. The posterior distribution is shown in Fig. 4, where it can be seen that the parametric uncertainties are much larger than in the posterior distribution from Bayesian History Matching (BHM). The effect arises as the likelihood has dominated in the posterior distribution, enlarging the area of probably parameter values due to the insensitive likelihood function [4,6] -which arises as a result of modelling the model discrepancy as a Gaussian process during joint inference approach.
The results from the three decoupled approaches and the hierarchical Bayesian approach are shown in Fig. 5, were it can be seen that all of the methods have managed to accurately predict the tip deflection well, reflected in low normalised mean squared errors (NMSEs) in Table 1 for a 200 point independent test dataset. However, due to the large parametric uncertainty  in the posterior distribution from the hierarchical Bayesian approach, the inferred model discrepancy not only has a large variance, but also a less accurate mean prediction with a NMSE of 39.349, over 17 times larger than the highest NMSE from the decoupled approach. This shows the challenges an insensitive likelihood causes on the inference process, and why a decoupled solution is one approach that can be used to overcome these challenges. Furthermore, the hierarchical Bayesian approach has an underestimated predictive variance for the tip deflection d tip , with a large number of data points, particularly around the first peak at 0.36 m, exceeding a three standard deviation interval. This relatively simple numerical case study shows the problems with a hierarchical Bayesian approach and further motives the need for alternative solutions to the model discrepancy inference problem, such as the decoupled approach proposed in the paper. In terms of comparing the three decoupled-based approaches, the main difference, as expected, is in the estimated uncertainty for the model discrepancy. Scenario one has the smallest uncertainty, with a larger number of experimental test data points outside a 3r range when compared to the other two scenarios. The first scenario is also overconfident in the model discrepancy predictions, especially around 0.1m, where the true model discrepancy is outside of the 3r range. Scenarios two and three increase the uncertainty in the model discrepancy, reflecting the parameter uncertainty in the posterior distribution, meaning the true model discrepancy remains within the 3r variance range. The NMSE is lowest for the model discrepancy in scenario three, with scenario two producing the largest error in its mean prediction. It can be argued from the results, that scenario one is overconfident and although the mean prediction is better than scenario two, its distribution could be misleading and less helpful to the engineer by not reflecting the true uncertainty in the analysis. Finally, the posterior of the Gaussian process hyperparameters is obtainable as part of scenario three, and presented in Fig. 6.

Numerical case study: two input problem
The second case study considers a scenario with multiple inputs, where X ¼ l loc ; l len fg i.e. crack location and length, where the output is tip deflection. The parameter in this analysis is the Young's modulus h ¼ E, where the true parameter value is 68GPa. The training inputs for both the simulator and experimental data were 64 data points evenly spaced between 0:1m and 0:5m, and 0mm and 50mm (for the crack location and length respectively) where the outputs are shown in Fig. 7. The simulator parameters were evaluated at four points between 50GPa and 90GPa and an emulator was constructed using a linear mean and Matérn 3/2 covariance function. The model discrepancy was defined as,  and displayed in Fig. 7. Again the experimental data, shown in Fig. 7, was formed from the simulator at the true parameters plus the model discrepancy with Gaussian additive noise with a variance of 0.001. The prior model discrepancy variance was V m ¼ 0:05. The approximate posterior from BHM is presented in Fig. 7 where the difference between the mode and true parameter value is 0.9%. Samples from the posterior are shown in Fig. 7 showing the simulator has been adequately calibrated.
The model discrepancy inference procedure was run for three scenarios:  1. Using a MAP estimate of the simulator parameters h MAP and an empirical Bayes estimate of the GP hyperparameters/ (when the simulator output for h MAP is used). 2. Marginalising out the simulator outputs y Ã using importance sampling, with an empirical Bayes estimate of the GP hyperparameters/ (when the simulator output for h MAP is used); N s ¼ 1000. 3. Marginalising out both the simulator outputs y Ã and GP hyperparameters / via importance sampling; N s ¼ 1000 and For each scenario a zero mean and Matérn 3/2 automatic relevance determination covariance function were used. The results in Fig. 8 show adequate predictions for all the scenarios with scenario 3 producing the lowest predictive NMSE on a 400 point independent test dataset (see Table 2). Interestingly, the second approach achieves a lower NMSE on the model discrepancy when compared to scenario three (with both performing better than scenario one). It can be seen in Fig. 8 that the model discrepancy has been correctly captured by all three approaches. Finally, scenario three obtained the posterior hyperparameter distribution, useful in understanding the extracted model discrepancy. (see Fig. 9).

Case study: five storey shear structure
In order to demonstrate the effectiveness of the proposed approach an experimental case study is presented. This case study seeks to infer the model discrepancy of a modal finite element model used to predict the change in natural frequency when different masses are applied to the fourth floor (simulating a damage scenario). Estimation of the parameter distribution was performed using Bayesian history matching as outlined in [8]. A brief overview of the calibration process is introduced below, where the reader is referred to [8] for more details.

Calibration using Bayesian history matching
Bayesian history matching was applied to infer the material properties h ¼ E; m; q fg (Young's modulus, Poisson's ratio and density) of a finite element model of a five storey shear structure (g Á; Á ðÞ ) known to have model-form errors due to modelling  simplifications (particularly of the boundary condition between the structure and fixing). Observational data were the first five bending modes of a representative building structure, z ¼ x 1 ; ...; x 5 f g , constructed from aluminium 6082, depicted in Fig. 10. These data were obtained via modal testing, where an electrodynamic shaker applied a Gaussian noise excitation with a bandwidth of 409.6 Hz, and five uniaxial accelerometers were used to capture the acceleration response at each of the five floors (where sample rate and time were chosen to allow frequency resolution of 0.05 Hz). Masses were incrementally added to the fourth floor of the structure m ¼ 0; 0:1; ...; 0:5 f g kg, representing pseudo-damage, and were treated as the inputs in this analysis i.e. m ¼ x. Ten repeat estimates of the natural frequencies were obtained for each mass providing a representation of observational uncertainty.
Calibration was performed on training data, which were the ten repeat observations of the bending natural frequencies when x ¼ 0; 0:3; 0:5 f g . The testing data were the ten repeat observations of the bending natural frequencies when x Ã ¼ 0:1; 0:2; 0:4 f g . The prior parameter bounds were AE15% of typical material properties for aluminium 6082; E ¼ 71GPa, m ¼ 0:33; q ¼ 2770kg/m 3 . These parameter bounds behave in a similar way to a uniform prior over the space. The approximate posterior distribution of the parameters, identified from the Bayesian history matching analysis, is displayed in Fig. 11.
Samples of the simulator output distribution (for each of the five natural frequencies) y  Fig. 12, where the error bars define to the prior model discrepancy and observational uncertainties. It is clear from Fig. 12 that there is a large amount of model discrepancy for the first natural frequency.

Model discrepancy inference
Inference of the model discrepancy from the Bayesian history matching analysis in Section 5.1 was performed using three approaches: 1. Using a MAP estimate of the simulator parameters h MAP and an empirical Bayes estimate of the GP hyperparameters/ (when the simulator output for h MAP is used). 2. Marginalising out the simulator outputs y Ã using importance sampling, with an empirical Bayes estimate of the GP hyperparameters/ (when the simulator output for h MAP is used); N s ¼ 1000. 3. Marginalising out both the simulator outputs y Ã and GP hyperparameters / via importance sampling; N s ¼ 1000 and For each of the three methods the model discrepancy was inferred as a map GP i : Y; X fg ! z i 8i 2 1 : 5 fg . It is noted that a multiple output Gaussian process could be implemented [33], meaning only one map would need to be inferred from Y; X fg to Z. This would not change the general formulation of the approach and is therefore left for further research. The priors for each of the five GP models were zero mean functions with Matérn 3/2 automatic relevance determination covariance functions (specified by Eqs. (3)-(5)).
For the third approach -the marginalisation of both the simulator outputs and GP hyperparameters -independent Gaussian priors were defined for each hyperparameter in the set. The priors for lengthscale hyperparameters assumed that the process will change slowly with the input (i.e. large lengthscales), log l The calibrated and bias-corrected predictions for each of the three approaches are displayed in Fig. 13, with the inferred model discrepancies shown in Fig. 14. Firstly, it can be seen that the mean predictions of all three approaches, both in terms of their output predictions' and inferred model discrepancies', are visually similar. The main difference between all three approaches is their estimation of the uncertainty, with significant differences in the inferred model discrepancy uncertainty.  The second and third methods have propagated the posterior parameter uncertainty through to the model discrepancy and the output predictions, unlike the first approach which collapses this uncertainty down to the parameter MAP estimates. This increase in uncertainty in the model discrepancy from methods two and three is useful to the engineer as it provides a better reflection of the underlying model discrepancy and will be helpful in identifying simulator improvements, as the 'true' model discrepancy is more likely to be contained within the confidence intervals. This is particularly clear given the results in the numerical case studies, where method one resulted in model discrepancy predictions where the 'true' model discrepancy  occasionally exceeded 3r. In terms of output prediction, each of the methods visually appear to have captured the noise process, with the third method showing increases in uncertainty outside of the training observations. This effect is likely to be caused by the small number of training observations, causing the uncertainty to increase away from the training data, indicating that the prior has a large effect on the posterior due to the small number of observations in the likelihood. However, the extra uncertainty quantified by marginalising out the posterior hyperparameter distributions is useful for gaining an insight into the level of trust in the identified model discrepancy given the limited training data used to estimate the discrepancy. Several validation metrics have been applied in order to quantitatively assess the performance of the inferred models for each of the three scenarios. The first metric, the normalised mean squared error (NMSE) (the sum of the squared errors divided by the variance and the number of data points), assesses the performance of the mean prediction. The second metric is the maximum mean discrepancy (MMD) distance, a measure of the distance between two distributions [34]. The distance is the difference between the means of two kernel embeddings of the data (where here a Gaussian kernel is used, where the scale parameter is inferred via a median heuristic [35]). The third metric is the posterior likelihood, and is a measure of the probability of the data coming from the inferred GP model. The validation metrics are applied to the predictions from each three scenario and are shown in Fig. 15 and Table 3.
In terms of the mean predictions, the NMSEs indicate that on average the third approach provided the best mean performance. In fact, both the second and third methods outperformed the first method in all of their mean predictions (apart from method three's predictions of the fourth natural frequency). This shows that including and propagating these sources of uncertainties are beneficial for the overall mean predictive performance, with the results supporting the conclusions from the numerical case studies. Helpfully, the models for all three scenarios show the same general mean predictive behaviour across the five natural frequencies, with predictions being poorest for the second natural frequency. Comparing the output distributions, the MMD distances for each scenario are relatively comparable, with the third method performing best on average. The reason for similar MMD distances is that the data distribution is being inferred from 10 observations at every input, and that this has a greater effect on the distance than each scenarios change in predictive distribution. This demonstrates the challenges in validating predictive distributions when only a few number of validation data points are available. In comparison to the NMSEs, the posterior likelihood indicates a different assessment of which natural frequency is predicted best -the second natural frequency is most likely to have produced the observational data. The posterior likelihoods indicate that on average the first method was more likely to have produced the observational data than the other approaches. However, the second method produces the highest posterior likelihoods for the third and fourth natural frequencies. Furthermore, the second and third methods have comparable posterior likelihoods for the first to third natural frequencies, which are very similar to method one. It is noted that the first method does not reflect the uncertainty associated with the parameters, and therefore may be overconfident in it's predictions, given the training observations.
From the third approach it is possible to obtain the posterior distribution of the GP hyperparameters. Fig. 16 depicts two of the posterior hyperparameters distributions; the first and fifth natural frequencies. The posteriors show that the lengthscales for the simulator output L y are all uncorrelated (as expected by the automatic relevance determination and distance metric assumptions in the covariance function, i.e. L y is diagonal). Furthermore, both show that the modes of the signal variance and noise variance are fairly constant when compared to the lengthscales. This means that the noise and signal variance have been well-identified, and the output uncertainty is mainly attributed the uncertainty in the lengthscales. Resultantly, the posterior hyperparameter distributions provide a high level of insight into the inferred model discrepancy. With more observational data, these posterior distributions will provide insight into the type of missing functional-form in the simulator, as the lengthscale distributions will be expected to decrease in uncertainty.

Conclusions
Every computer model (here defined as a simulator) will imperfectly reflect the real world due to some level of model discrepancy (whether due to missing physics, simplifications or approximations etc.). Without identifying the level of model discrepancy within a simulator, predictions are likely to be inaccurate. This paper proposes a method based on Gaussian process (GP) regression and a sampling-based approach for identified model discrepancy given some parameter distribution. This method has been demonstrated to be effective on numerical examples and an experimental case study of a representative five storey building structure.
The approach in this paper allows for bias-corrected predictions to be constructed that marginalise out the simulator outputs, with the additional ability to marginalise out the GP hyperparameters. By performing this process, the bias-corrected predictive distributions better reflect the parameter and hyperparameter uncertainty in the predictions and help the engineer identify 'true' improvements to a simulator, rather than those based on overconfident estimations of the model discrepancy. The technique relies on generating a set of Gaussian process maps from the uncertain simulator outputs and deterministic inputs to observational data and performing weighted averages to form the bias-corrected predictive distributions.
Three scenarios were investigated: using a MAP estimate of the simulator parameters and simulator outputs, and an empirical Bayes estimate of the GP hyperparameters; marginalising out the simulator outputs using importance sampling, and an empirical Bayes estimate of the GP hyperparameters; and marginalising out both the simulator outputs and the GP hyperparameters via importance sampling.
The numerical case studies show that a two stage decoupled process, utilising Bayesian history matching and the proposed model discrepancy procedure, is appropriate for calibrating a simulator and extracting model discrepancy. In addition, the first numerical case study demonstrated issues associated with the hierarchical Bayesian approach that seeks to jointly infer the parameters and model discrepancy. In accordance with the literature [4,6], the approach leads to an insensitive likelihood, and can cause non-identifiability issues. In the numerical case study the hierarchical Bayesian approach inferred a model discrepancy distribution that was more uncertain than the proposed decoupled approaches, and had a worse mean predictive performance. Furthermore, both numerical case studies showed that considering both the simulator output  and hyperparameter uncertainties provides more information about the model discrepancy function, and in improves the performance of the mean prediction.
Finally, an experimental case study was provided, where it was shown that the third approach (marginalising both the simulator outputs and the GP hyperparameters) produced the best mean predictions. In addition, the uncertainty associated with the parameters was propagated onto the model discrepancy and output predictions, better reflecting the uncertainty quantified by the BHM parameter posterior distribution. The inclusion of the parametric uncertainty is valuable in understanding the 'true' model discrepancy, and is beneficial in determining what is known about the model discrepancy from the analysis. In addition, the third approach provides posterior distributions of the hyperparameter, which provide further insight into the model discrepancy process.

Declaration of Competing Interest
None.