Decoding in autoregressive models (ARMs) consists of searching for a high scoring output sequence under the trained model. Standard decoding methods, based on unidirectional greedy algorithm or beam search, are suboptimal due to error propagation and myopic decisions which do not account for future steps in the generation process. In this paper we present a novel decoding approach based on the method of auxiliary coordinates (Carreira-Perpinan & Wang, 2014) to address the aforementioned shortcomings. Our method introduces discrete variables for output tokens, and auxiliary continuous variables representing the states of the underlying ARM. The auxiliary variables lead to a factor graph approximation of the ARM, whose maximum a posteriori (MAP) solution is found exactly using dynamic programming. The MAP solution is then used to recreate an improved factor graph approximation of the ARM via updated auxiliary variables. We then extend our approach to decode in an ensemble of ARMs, possibly with different generation orders, which is out of reach for the standard unidirectional decoding algorithms. Experiments on the text infilling task over SWAG and Daily Dialogue datasets show that our decoding method is superior to strong competing decoding methods.
I think this is a good study, unless I miss something. It proposes a new solution to one of the fundamental problems in the field, which significantly improves previous solutions in terms of accuracy. I will recommend it for acceptance unless I miss something. (I’m not an expert in the field and the problem seems to be so fundamental, making me cautious about judging the novelty of the proposed solution.)
There seems to be some unclarity about the optimization algorithm. In short, I suspect that the proposed optimization has some difficulty with its convergence. 
-	In Sec.4, the authors suggest “a two-step block coordinate descent algorithm to alternate between (i) optimizing y’s while g’s are fixed, and (ii) optimizing g’s while y’s are fixed”. They further suggest `a variant of the Viterbi algorithm’ for (i) and gradient-based algorithms for (ii). These seem to make sense to me so far. 
-	Then, in Sec.6, they state that “We use the Nesterov optimizer with a learning rate of 0.1. All \mu’s are initialised with 0.5 and are multiplied by 1.2 after 5 epochs”. I think this needs some explanation. 
-	For instance, what is the thought behind the choice of the optimizer? Why was Nesterov’s acceleration necessary instead of plain GD? Is it essential to use the scheduled adjustment of \mu?
-	It states “All \mu’s” but there seems only a single \mu in Eq.(5). Am I missing something?
-	 Some explanation on what “Nesterov optimizer” and \mu would also be necessary for a wide range of readers. 
-	There is a statement on the initialization of the optimization in p.7: “h corresponding to the beam search solution to initialize the g variables”. How sensitive is the optimization to the initial values? For instance, will the results change if g’s are initialized by the solution of the greedy method. 
-	There is a plot without a figure caption in p.8, which shows how the objective cost decreases with parameter updates. Why are their values at ‘epoch’=0 lower than those at the subsequent epochs? Does this mean the initial values give more optimal parameters?
Additional comments:
-	It would be more friendly to the readers to show the definition of the notations like y^n_1. 
-	There is a statement like “Decoding was run for 10 epochs.” I suppose epochs here mean iteration count of the alternated parameter updates in the proposed algorithm. Why do the authors call it `epoch’? 
-	No figure caption for the plot in p.8. It also uses `epochs’ for the horizontal axis title. No axis title for the vertical axis.