extfieldsofstudy:
Computer Science
Mathematics

text:
Feasibility Conditions for Interference Alignment

The degrees of freedom of MIMO interference networks with constant channel coefficients are not known in general. Determining the feasibility of a linear interference alignment solution is a key step toward solving this open problem. Our approach in this paper is to view the alignment problem as a system of bilinear equations and determine its solvability by comparing the number of equations and the number of variables. To this end, we divide interference alignment problems into two classes - proper and improper. An interference alignment problem is called proper if the number of equations does not exceed the number of variables. Otherwise, it is called improper. Examples are presented to support the intuition that for generic channel matrices, proper systems are almost surely feasible and improper systems are almost surely infeasible.

I. INTRODUCTION
The degrees of freedom (DoF) of wireless interference networks represent the number of interference-free signalingdimensions in the network. In a network with K transmitters and K receivers and non-degenerate channel conditions, it is well known that K non-interfering spatial signaling dimensions can be created if the transmitters or the receivers are able to jointly process their signals. Cadambe and Jafar [1] recently introduced the idea of interference alignment for the K-user wireless interference network with time-varying/frequencyselective channel coefficients, and showed that K/2 spatial signaling dimensions are available inspite of the distributed nature of the network which precludes joint processing of signals at transmitters or receivers. While a number of interference alignment solutions have appeared since [1] for different channel settings, many fundamental questions remain unanswered. One such problem is to determine the feasibility of linear interference alignment for MIMO interference networks with constant channel coefficients. It is this open problem that we address in this paper.
Suppose we denote by (M × N, d) K , the K-user interference network where every transmitter has M antennas, every receiver has N antennas and each user wishes to achieve d DoF. We call such a system a symmetric system. Consider the following examples.
• (2 × 2, 1) 3 -It is shown by [1] that the in the 3user interference channel with 2 antennas at each node, each user can achieve 1 degree of freedom by linear † The author is supported in part by The Scientific and Technological Research Council of Turkey (TUBITAK). The author is on leave at University of California, Irvine.
interference alignment, i.e. by linear beamforming at the transmitters and linear combining at the receivers. • (5×5, 2) 4 -Consider the 4-user interference channel with 5 antennas at each user. Suppose we wish to achieve 2 DoF per user, for a total of 8 network DoF. An analytical solution to this problem is not known but numerical evidence in [2] clearly indicates that a linear interference alignment solution exists. Numerical algorithms are one way to determine the feasibility of linear interference alignment. But, is there a way to analytically predict the feasibility of alignment? In other words, without running the numerical simulation could we have predicted whether a linear interference alignment solution will exist for the (5 × 5, 2) 4 system? • Now consider three distinct systems -(6 × 4, 2) 4 , (7 × 3, 2) 4 , (8×2, 2) 4 . Are these systems feasible? Clearly, the last one, (8 × 2, 2) 4 is feasible, because simple transmit zero-forcing is enough to eliminate the interference at every receiver. Is the feasibility of (8 × 2, 2) 4 case related to the feasibility of other cases where each transmitter successively donates one antenna to the receiver? We will show that these three systems and the (5 × 5, 2) 4 case belong to the same group.
In this paper, we address all these questions, as well as more complex asymmetric cases where each node may have different number of antennas and each user may demand different number of degrees of freedom. The basic approach is to consider the linear interference alignment problem formulation as a system of bilinear equations. We determine the correct way to count the number of variables and equations for a general MIMO interference alignment problem. Then, based on the number of variables and equations we classify the system as either proper (number of equations does not exceed number of variables), or improper (number of equations exceeds the number of variables), with the intuitive understanding that proper systems are almost surely feasible, and improper systems are almost surely infeasible.
Aside from the detailed conditions that classify a general (asymmetric) system as proper or improper, our analysis yields a useful rule of thumb for symmetric systems (M × N, d) K . We find that donating an antenna from each transmitter to the corresponding receiver, or vice versa, does not change the nature of the system (proper or improper), provided every node still has at least d antennas. For example, consider the system (2 × 3, 1) 4 . Is this system proper? While this question may be difficult to answer at first, now suppose we transfer one antenna from each transmitter to its corresponding receiver, to obtain the system (1 × 4, 1) 4 . This system is clearly proper, because simple zero-forcing achieves 1 degree of freedom for each user in this system. This tells us that the system (2 × 3, 1) 4 is also proper. Thus, the following are groups of proper systems, related by the rule of thumb defined above.
The first and last members of each group are easily seen to be proper (because a simple zero-forcing solution exists), thereby also determining the status of the rest of the members of the group. Improper systems can be similarly grouped as well. Completely asymmetric cases require a more sophisticated set of feasibility conditions. The analysis is supported by numerical results for a wide variety of cases, including the specific examples listed above. We begin with the system model.

II. SYSTEM MODEL
We consider the same K-user MIMO interference channel as considered in [2]. The received signal at the n th channel use can be written as follows: ) denotes the Kuser MIMO interference network, where the k th transmitter and receiver have M k and N k antennas, respectively and the k th user demands d [k] DoF. As defined earlier, (M × N, d) K denotes the K-user symmetric MIMO interference network, where each transmitter and receiver has M and N antennas, respectively, and each user demands d DoF, so that the total DoF demand is Kd. Some sample asymmetric and symmetric systems are shown in Fig. 1.

III. LINEAR INTERFERENCE ALIGNMENT SCHEME
In interference alignment precoding, the transmitted signal from the k th user is vector that denotes the d [k] independently encoded streams transmitted from the k th user. The M [k] ×d [k] precoding filters V [k] are designed to maximize the overlap of interference signal subspaces at each receiver while ensuring that the desired signal vectors at each receiver are linearly independent of the interference subspace. Therefore, each receiver can zeroforce all the interference signals without zero-forcing any of the desired signals. The zero-forcing filters at the receiver are denoted by U [k] . In [2], it is shown that an interference alignment solution requires the simultaneous satisfiability of the following conditions: where † denotes the conjugate transpose operator. Very importantly, [2] explains how the condition (2) is automatically satisfied almost surely if the channel matrices do not have any special structure, ) and are designed to satisfy (1), which is independent of all direct channels H [kk] . The scenarios where it is difficult to theoretically determine the feasibility of interference alignment can be evaluated numerically using an iterative algorithm proposed in [2]. In this work, we develop analytical criteria to predict the feasibility of interference alignment. Our approach is to count the number of equations and variables in (1). We assume generic MIMO channels with no structure and force the required ranks of the transmit and receive filters by design. Thus, (2) is satisfied automatically and we only need to count the number of variables and equations for (1).

IV. PROPER SYSTEM
The notion of a proper network is based on counting the number of equations and the number of variables involved in the bilinear equations (1). While the formal definition appears later, put simply, the system network is proper if and only if the cardinality of every subset of the equations is less than or equal to the number of variables involved in that subset of equations. Otherwise, the network is considered improper. The reason for this classification is the following intuition that forms the basis of our approach in this paper.
Key Insight : The interference alignment problem is almost surely feasible for proper systems and almost surely infeasible for improper systems.
Instead of a formal proof, we refer to Bezout's theorem, which states that a system of m generic polynomial equations in m variables will almost surely have as many common solutions as the product of the degrees of the polynomials. For our case, it suffices if there is one solution. The feasibility conditions in our case correspond to a system of polynomial equations, whose generic nature is due to the generic channel matrices. Due to the special bilinear form of the equations it needs to be shown that the polynomials are sufficiently generic. We are able to construct the rigorous proof for the (2×3, 1) 4 case which is omitted here due to lack of space. We expect generalizations to all cases should be possible, albeit cumbersome.
Next we explicitly account for all equations and variables. Let us start with the total number of equations N e and the total number of variables N v .

A. Counting the Number of Equations N e and Variables N v
To obtain N e and N v , we rewrite the condition in (1) as follows: m are the transmit and the receive beamforming vectors (columns of precoding and interference suppression filters, respectively). N e is directly obtained from (3) as follows: However, calculating the number of variables N v is less straightforward. In particular, we have to be careful to not count any superfluous variables that do not help with interference alignment. At Thus, the columns of V [k] are the basis for the transmitted signal space. However, the basis representation is not unique for a given subspace. In particular, consider any full rank Then vectors. It is easy to argue that there is no other basis representation for the transmitted signal space with fewer variables.
Therefore, by eliminating all superfluous variables for the interference alignment problem, the number of variables to be designed for the precoding filter of the k th transmitter, [k] . Likewise, the actual number of variables to be designed for the interference suppression filter of the k th receiver, . As a result, the total number of variables in the network to be designed is:

B. Proper System Characterization
To formalize the definition of a proper system, we introduce some notation. We use the notation E mn kj to represent the equation The set of variables involved in an equation E are indicated by the function var(E). Clearly where | · | is the cardinality of a set. Using this notation, we denote the set of N e equations as follows: This leads us to the formal definition of a proper system.
In other words, for all subsets of equations, the number of variables involved must be at least as large as the number of equations in that subset. Condition (14) provides us a way to predict the feasibility of interference alignment in a general ) system. However, note that it can be computationally cumbersome because we have to test all subsets of the set of all equations. In most cases however, especially if the system is improper, simply comparing the total number of equations and the total number of variables may suffice.
Theorem 1: Example 1: Consider the (5 × 5, 3)(5 × 5, 2) 3 system, i.e. a 4-user interference networks where all nodes have 5 antennas, user 1 demands 3 DoF and users 2, 3, 4 demand 2 DoF each. There are a total of 60 equations, so that there are 2 60 − 1 subsets of the set of all equations. Testing each of them could be very challenging. However, since the total number of variables N v = 48 is less than the number of equations, the system is easily seen to be improper.
Similarly, one can sometimes identify the bottleneck equations in the system by checking the equations with the fewest number of variables, i.e., the equations involving the transmitters and receivers with the fewest antennas.
Example 2: As a simple example, consider the system (2 × 1, 1) 2 which is clearly feasible (proper) because simple zeroforcing is enough for achievability. However, now consider the (2 × 1, 1)(1 × 2, 1) system, which also has the same total number of equations N e and variables N v as the (2 × 1, 1) 2 system. Thus, only comparing N v , N e would lead one to believe that this system is proper. However, suppose we check E 11 12 which connects transmitter 2 and receiver 1, both of which have only one antenna each, i.e., let S = {E 11 12 } so that |S| = 1. However var(E 11 12 ) = 0. Thus, this system has an equation with zero variables, which makes it improper, and therefore infeasible.
Example 3: Several interesting cases emerge from applying the condition (15). For example, consider the 2 user interference channel (2 × 3, 1)(3 × 2, 1) where a total of 2 degrees of freedom are desired. It is easily checked that this system is proper and the achievable scheme is described in [3]. However, now consider the 4-user interference channel consisting of two sets of these channels, all interfering with each other (2 × 3, 1) 2 (3 × 2, 1) 2 which is a 4-user interference channel and a total of 4 degrees of freedom are desired. It is easily verified that this is a proper system, i.e. it satisfies (15). It turns out that a closed form solution for alignment can be found in this case. It is somewhat surprising that going from two users to four users simply doubled the degrees of freedom in this case, without any penalty for interference between the users in terms of degrees of freedom. Since in this paper our focus is only on feasibility and not on closed form solutions, we omit the details of this example.
In general, testing if a large asymmetric system is proper can be cumbersome. However, for the case of symmetric systems of the type (M × N, d) K , the process is much simpler. Remark: Theorem 2 implies that for every user to achieve d degrees of freedom in a K user interference channel, it suffices to have a total of M + N ≥ (K + 1)d antennas between the transmitter and receiver of a user. The antennas can be distributed among the transmitter and receiver arbitrarily, as long as each of them has at least d antennas. In particular, to achieve K degrees of freedom in a K user symmetric network we only need a total of K +1 antennas between the transmitter and receiver of each user. The system (2 × 3, 1) 4 is such an example.
Proof: Because of the symmetry each equation is involved with the same number of variables and any deficiency in the number of variables shows up in the comparison of the total number of variables versus the total number of equations. Plugging in the values of N v , N e computed earlier, we have the result of Theorem 2.
The following corollary shows the limitations of linear interference alignment over constant MIMO channels (with no symbol extensions).
Corollary 1: The ratio of the sum degrees of a proper (M × N, d) K system, normalized by a single user's degree of freedom in the absence of interference is bounded by: Proof: The proof is straightforward from the condition of Theorem 2.
Remark: For the case that M = N note that the total number of degrees of freedom for a proper system is no more than twice the number of degrees of freedom achieved by each user in the absence of interference. Note that for diagonal (time-varying) channels it was shown in [1] that the total number of degrees of freedom is K/2 times the number of degrees of freedom achieved by each user in the absence of interference. This result shows that the diagonal structure of the channel matrix is very helpful. Going from the case of no structure (arbitrary MIMO channels) to diagonal structure the ratio of total degrees of freedom to the single user degrees of freedom increases from a maximum value of 2 to K/2.
The following corollary identifies the groups of symmetric systems, which are either all proper or all improper. Example 7: The systems (1 × 4, 1) 4 , (2 × 3, 1) 4 , (3 × 2, 1) 4 , (4 × 1, 1) 4 are in the same group, formed by switching between transmit and receive antennas. It is easy to see that the (4 × 1, 1) 4 system is proper, because simple zero-forcing suffices to achieve the DoF demand. By virtue of being in the same group, the rest are proper as well.
Example 8: Consider a (2 × 8, 2) 4 symmetric system. Again, the DoF for this interference channel can be trivially obtained by zero-forcing at each receiver. By switching antennas from each transmitter to receiver, an equivalent (5 × 5, 2) 4 symmetric system of Example 6 is obtained, which is also proper.

V. NUMERICAL RESULTS
We tested numerous interference alignment problems, both symmetric and asymmetric, and especially including each of the examples presented in this paper, using the numerical algorithms of [2]. In every case so far, we have found the results to be consistent with the guiding intuition of this work, proper systems are feasible and improper systems are not.
In this section we provide numerical results for a few interesting and representative cases. The results are in terms of the leakage interference, as defined in [2] -i.e. the fraction of the interference power that is present in the dimensions reserved for the desired signal. The interference percentage at the k th receiver is evaluated as follows: where λ j denotes the smallest eigenvalue of a matrix, Tr denotes the trace of a matrix and Q [k] denotes the interference covariance matrix at the k th receiver: The numerator and the denominator of (18) are the interference and desired signal space powers at the k th receiver, respectively.
In Fig. 2, the interference percentages versus the total number of beams are presented. The total number of beams are started from the DoF of each network. Therefore, after the first point on the x-axis, interference percentage of each network is not zero. The non-zero interference percentage indicates that interference alignment is not possible.

VI. CONCLUSION
We propose an analytical method to predict the feasibility of interference alignment via linear schemes over constant MIMO channels that have no structure (i.e. no symbol extensions). Our approach to determine feasibility is to count the number of equations and variables. We define a system as proper if the number of variables is not smaller than the number of equations, and as improper otherwise, with the understanding that proper systems are feasible almost surely, while improper systems are almost surely infeasible. This observation follows from Bezout's theorem and the generic nature of the channel coefficients, and can be rigorously shown in some cases, to be included in the full paper. The full paper will also include closed form solutions for interference alignment for several cases, such as (2 × 3, 1) 2 (3 × 2, 1) 2 , (3 × 3, 1)(2 × 3, 1) 3 , (2 × 4, 1)(2 × 3) 3 , (2 × 4, 1) 2 (4 × 2, 1) 2 (3 × 3, 1) etc.

ACKNOWLEDGMENT
The first author would like to thank V. R. Cadambe at University of California, Irvine for valuable discussions on interference alignment.


extfieldsofstudy:
Computer Science
Medicine

text:
Functional identification of biological neural networks using reservoir adaptation for point processes

The complexity of biological neural networks does not allow to directly relate their biophysical properties to the dynamics of their electrical activity. We present a reservoir computing approach for functionally identifying a biological neural network, i.e. for building an artificial system that is functionally equivalent to the reference biological network. Employing feed-forward and recurrent networks with fading memory, i.e. reservoirs, we propose a point process based learning algorithm to train the internal parameters of the reservoir and the connectivity between the reservoir and the memoryless readout neurons. Specifically, the model is an Echo State Network (ESN) with leaky integrator neurons, whose individual leakage time constants are also adapted. The proposed ESN algorithm learns a predictive model of stimulus-response relations in in vitro and simulated networks, i.e. it models their response dynamics. Receiver Operating Characteristic (ROC) curve analysis indicates that these ESNs can imitate the response signal of a reference biological network. Reservoir adaptation improved the performance of an ESN over readout-only training methods in many cases. This also held for adaptive feed-forward reservoirs, which had no recurrent dynamics. We demonstrate the predictive power of these ESNs on various tasks with cultured and simulated biological neural networks.

Introduction
One central goal in neuroscience is to understand how the brain represents, processes, and conveys information. Starting from Hebb's cell assemblies (Brown and Milner 2003;Hebb 1949), many neurobiologically founded theories and hypotheses have been developed towards this goal. It stands clear now that spikes are the elemental quanta of information processing in the mammalian cortex (Hodgkin and Huxley 1952;Kandel et al. 2000). As a result of extensive experiments of cortical recordings, it has been widely postulated and accepted that function and information of the cortex are encoded in the spatio-temporal network dynamics (Abeles et al. 1993;Lindsey et al. 1997;Mainen and Sejnowski 1995;Prut et al. 1998;Riehle et al. 1997;Villa et al. 1999). The right level of describing the dynamics, however, is a matter of intensive discussions (Prut et al. 1998;Rieke et al. 1999). Are spike rates or spike timings more relevant? What is the right temporal precision if the latter proves significant? What should be the spatial resolution of this description? How far can the population activity of neurons be related to function or behavior? Does the correlated activity of multiple neurons indicate a functionally relevant state?
Depending on the answers to the above questions one would preferably apply different models to relate the network activity to function. Another approach is to employ a generic network model, which can be assumed to be universal for problems of neural encoding. The parameters of the model would be learned by adaptive algorithms. Obviously, such a model should be able to deal with single spikes with high temporal precision as well as population rates. It should also be able to catch, with the appropriate parameters, network synchrony and polychrony (Izhikevich 2006).

Reservoir computing
Liquid State Machines (LSM) and Echo State Networks (ESN) have been introduced as efficiently learning recurrent neural networks (Jaeger 2001;Maass et al. 2002). The common key contribution of these approaches is the proposal of a recurrent neural network with a fixed connectivity, i.e. a reservoir, which does not have stable states and has a fading memory of the previous inputs and network states. In response to an input stream, the reservoir generates a higher-dimensional spatio-temporal dynamics reflecting the structure in the input stream. The higher dimensional reservoir state can be mapped to a target output stream online, with a second module, namely a readout. LSMs are networks of spiking integrate-and-fire (IAF) neurons whereas ESNs use continuous valued sigmoid neurons and a single layer of readout neurons (see Section 2).
With the appropriate parameters, reservoir dynamics can be sensitive to different features of the input such as correlations, polychronous and synchronous spikes, different frequency bands or even temporally precise single spikes. For instance, LSMs can approximate any time invariant filter with fading memory if their traces of internal states differ at least for one time point in response to any two different input streams (Separation Property, SP) and if the readout modules can approximate any continuous function from R m → R, where m ∈ N (Approximation Property, AP) (Maass et al. 2002). Satisfying the separation property depends on whether the reservoir is composed of sufficient basis filters. A random reservoir needs to have a rich repertoire of basis filters in order to approximate the target time-invariant filter. This could be achieved for several tasks with sufficiently large random reservoirs (Maass et al. 2002). Furthermore, reservoirs have been shown to simulate a large class of higher order differential equations with appropriate feedback and readout functions (Maass et al. 2007). These findings suggest that reservoir computing can be used as a generic tool for the problems of neural encoding.
Biological neural networks, in vivo, process a continuous stream of inputs in real time. Moreover, they prove successful to react and compute fast, independent of their instantaneous state and ongoing activity, when prompted by sudden changes in stimuli. In other words, they perform any time computing. Reservoir computing has also been suggested as a model of cortical information processing for their capability of online and anytime computing, for their fading memory and for their separation properties. It has been argued that specific cortical circuitry is possible to build into generic LSM framework (Maass et al. 2004). Bringing reservoir computing into the problem will not only deliver expressive models that can distinguish a rich set of input patterns, but also may provide more biological relevance to the theoretical tool.

Neuronal cell cultures as biological neural networks
While brain tissue has highly specialized architecture and developmental history, generic biological networks can be created as cell cultures of mammalian cortical neurons that have been dissociated and regrown outside an organism. They are closed system in vitro living networks, which are frequently used to study physiological properties of neural systems (Marom and Shahaf 2002). Although the anatomical structure of the brain is not preserved in cultures, their inherent properties as networks of real neurons, their accessibility and small size make them very attractive for investigating information processing in biological networks. Using cultured networks also eliminates the problem of interference with ongoing activity from different parts of the brain. Compared to in vitro brain slices, cultured networks are more accessible in terms of ratio of the recorded neurons to the whole network. In other words, they pose a less serious under-sampling problem. Studying such networks can provide insight into the generic properties biological neural networks, independent of a specific anatomy (Marom and Shahaf 2002 to be able to learn. Shahaf and Marom (2001) demonstrated how the response of the network to a stimulus can be 'improved' by systematically training the culture. Jimbo et al. (1999) showed how neuronal cultures can increase or decrease their overall response to pulse inputs stimulus by tetanic stimuli training. Micro-Electrode Arrays (MEA) have commonly been employed to both stimulate and record from neuronal cultures and slices (Egert et al. 2002;Marom and Shahaf 2002) (Fig. 1). Standard MEAs allow for simultaneous recordings and stimulations from up to 60 electrodes from a surface area of around 2 mm 2 . Each electrode picks up the extracellular electrical field of one or several (1 to 3) neurons. Neurons transmit information by action potentials, or spikes, which can be extracted by adaptive high-pass filtering of the extracellular signals (Egert et al. 2002;Wagenaar et al. 2005).
The activity in neuronal cultures is composed of irregular network-wide bursts of spikes, even in absence of an external stimulation (Marom and Shahaf 2002). These networks display little or no spiking activity most of the time (that is, between bursts) and very high spiking activity during the bursts (Fig. 2). Although the inter-burst intervals are not necessarily regular, there seem to be spatio-temporal patterns that occur within bursts (Feber le et al. 2007;Rolston et al. 2007). For example, a burst might always start with the activity of the same channel and continue with the activity of another particular channel.

Problem statement
Although we cannot relate the activity dynamics to a physiological function in random in vitro BNNs, by studying their activity dynamics, we gain experience and information about generic network properties forming the basis of in vivo networks. Moreover, one can assign pseudo-functions to random BNNs by artificially mapping network states to a predefined set of actions (Chao et al. 2008). One can also regard the response spike train of a BNN to a stream of various stimuli as a very detailed characterization of its pseudofunction and aim at modeling stimulus-response relations. In the present work, we take this approach. We record the output responses of simulated and cultured BNNs to random multivariate streams of stimuli. We tackle the question whether it is possible to train an artificial neural network that predicts the response of a reference biological neural network under the applied stimulus range. In other words, we aim at generating an equivalent network of a BNN in terms of stimulus-response relations. Given the same stimulus, the equivalent network should predict the output of the biological neural network.  A model or a predictor for biological neural networks can be useful for relating the physiological and physical determinants of its activity and thereby, can be a tool for analyzing information coding in these networks. It can also be helpful for interacting with BNNs by means of electrical stimulation. Here, we employ Echo State Networks (ESN) as a reservoir computing tool and propose an algorithm to find appropriate models for the relations between continuous input streams and multi-unit recordings in biological neural networks. The algorithm uses point process log-likelihood as an optimization criterion and adapts the readout and reservoir parameters of ESNs accordingly. Moreover, we investigate the performance of our approach on different feed-forward and recurrent reservoir structures and demonstrate its applicability to stimulus-response modeling problem in BNNs.
We shortly review ESNs in Section 2 and point process modeling of spike data in Section 3. An elaboration of ESN adaptation for point processes is presented in Section 4. In Section 5 we present our evaluation methods. A detailed experimental section and their implications can be found in Sections 6 and 7, respectively.

Echo State Networks
Echo State Networks (ESN) (Jaeger 2001) consist of a recurrent neural network (RNN) reservoir and a readout module (Fig. 3). Recurrent neural networks (RNN) are in general difficult to train with gradient descent methods (Bengio et al. 1994). The echo state approach to RNNs is motivated by the observation that a fixed (untrained) RNN can be very useful for discriminating multivariate time series if the RNN is gradually forgetting internal states and thereby also the inputs at previous time steps. This property is described as generating the echo states of the input. Existence of  states is independent of the specific structure of the reservoir, but strongly determined by the extension of the eigenvalues of its connectivity matrix. The input stream is used to stimulate a higher dimensional dynamical system, which can be mapped to a target output stream by adaptive linear or sigmoidal readout units. Even with random fixed reservoirs, ESNs make up a successful reservoir computing framework with applications in several engineering and nonlinear modeling tasks (Jaeger 2003;Jaeger and Eck 2006;Jaeger and Haas 2004).

Dynamics
The ESN dynamics is described as N × 1 state vector for time step n, u n K × 1 input vector for time step n, y n L × 1 output vector for time step n, ';' vertical vector concatenation.
Note that ESNs can have feedback projections from the readout module to the reservoir (Jaeger 2001). Throughout this note, however, we employ ESN architectures without such feedback projections. We choose exp(x)+exp(−x) , which allows for existence of echo states.

ESN with leaky integrators
Fading memory in the reservoir is due to recurrent connectivity if the individual neurons are memoryless. Using leaky integrator neurons results in a reservoir made from a network of low pass filters, which can be referred to as Leaky Integrator Echo State Networks (LI-ESN). LI-ESNs have been shown to work well on noisy systems with slowly changing dynamics (Jaeger et al. 2007). The leakage time constant is another aspect of fading memory. For instance, even a purely feedforward reservoir will have a fading memory due the memory of the individual neurons. We describe the LI-ESN dynamics with y n = f out W out u n ; x n , where τ is a global parameter for all reservoir neurons.

Feed-Forward Echo State Networks
Acyclic random neural networks with leaky integrator neurons can be used as reservoirs, as they also have fading memory. We refer an ESN with such a reservoir as a Feed-Forward Echo State Network (FF-ESN). Naturally, it can be asked whether FF-ESNs can have the same expressive power as the recurrent ESNs. Intuitively, it can be expected that a recurrent ESN with the same number of neurons discriminates between more input patterns than a FF-ESN, i.e. the recurrent ESN has a much more powerful separation property. In case of adaptive reservoirs, however, it is worthwhile to test whether the reservoir adaptation algorithm works especially better on feed-forward reservoirs than on recurrent ones. Here, we also propose the use of FF-ESNs for functional identification of biological neural networks. A recurrent reservoir can easily be made acyclic by inverting the connections that cause recurrence (Fig. 4).

ESN learning
Learning in Echo State Networks will typically include learning readout parameters by linear regression, e.g. using the pseudo inverse matrix where L×T matrix Y is the collection of target vectors y for T time steps, X , U are the same collection for reservoir states and input vectors. + is pseudo inverse operation. Much more efficient is the use of the Wiener-Hopf equation for ESN readout learning (Jaeger 2001). Steil (2004) proposed the backpropagationdecorrelation algorithm, which seeks a trade-off between error minimization and decorrelation of the reservoir activity. Although ESN learning is restricted to readout learning in many cases, there have recently been several approaches to adapt reservoirs, among which a significant improvement over untrained reservoirs was reported by Steil (2007). Adopting the intrinsic plasticity rule from real biological systems improved the performance of the backpropagation-decorrelation algorithm.
In the current approach, we adapt the reservoir connectivity and individual neuronal time constants with a one-step propagation of the log-likelihood into the reservoir. This is basically a gradient descent approach, where the log-likelihood of a point process is used as an optimization criterion. We elaborate our approach in Section 4.

Point process modeling of spike data
Signals that are comprised of point-like events in time, e.g. trains of action potentials generated by neurons, can be characterized in terms of stochastic point processes (Brown et al. 2001;Chornoboy et al. 1988;Okatan et al. 2005;Rajaram et al. 2005). A spiking process, modeled as a point process, is in turn fully characterized by its conditional intensity function (Cox and Isham 1980;Daley and Vere-Jones 2003). The instantaneous firing rate is defined as a conditional intensity function where N i (t) is the cumulative spike count of unit i. λ i (t|I t , H t ) represents the conditional probability density that a spiking event occurs at unit i. I t and H t stand for input history and spike history, respectively. For small , the probability of a spike can be computed by Bernoulli approximation to the point process (Daley and Vere-Jones 2003;) H t and parameter set θ, the log-likelihood L of a sample count path for a node (electrode) i is expressed as in Chornoboy et al. (1988): where p(N i |θ) is the conditional probability density for the count path. The instantaneous log-likelihood (Brown et al. 2001) of a point process trajectory (spike train) is defined as Note that With a conversion to discrete time, i.e. t → n, → 1, instantaneous log-likelihood becomes, where n is the discrete time step index. Any parameter θ j of the model is then learned then by the standard gradient descent learning.
where η is the learning rate. For the sake of readability, we will leave the calculation points of partial derivatives out and simply use the notation θ n+1 The quality of point process modeling strongly depends on the expressive power of function λ(t|I t , H t , θ). Linear or log-linear models might be preferable for successful adaptation by gradient descent, whereas nonlinear functions allow for more expressive models. In the current work, we model this function by an Echo State Network, which inherently incorporates input and network history into the instantaneous network states.

Point process modeling with ESN and log likelihood propagation
We relate input and reservoir state to the conditional intensity function in discrete time as where λ n is an L-dimensional vector of conditional intensity estimations for time step n. We use an exponential f out where 0 ≤ A ≤ 1 is a constant. Exponential functions were already applied in point process modeling and shown to have desirable properties such as avoiding local minima in Generalized Linear Models (GLM) (Paninski 2004). ESN readout learning can be trivially adapted to point process data, where the output of the ESN estimates the conditional intensity function in batch mode and also in online mode We are going to provide the expansion of the above learning rule and its extension to reservoir adaptation assuming an online mode. Combining Eqs. (2), (3), and (4) yields where ξ n = W out [u n ; x n ].

Reservoir adaptation
One interesting question is whether it is sufficient to learn the output weights or whether one needs to adapt the reservoir using the point process log-likelihood as the fitness criterion. In the current work, we adapt recurrent and feed-forward reservoirs with one step propagation of the point process log-likelihood into the reservoir and compare the results to non-adaptive reservoirs. Any connectivity weight W kl within the reservoir can also be adapted with the gradient descent rule, As we assume the output channels are mutually independent given the parameters, Obviously, n i has a long term dependency on W kl due to reservoir memory, i.e. n i not only depends on W n−1 kl but also on W n−h kl , where h > 1. However, for the sake of computational efficiency, and in view of the fact that long term partial derivatives tend to vanish exponentially with respect to time (Bengio et al. 1994), we utilize a one-step propagation of the instantaneous log-likelihood into the reservoir.
Note that if neuron ψ projects to the readouts, any connection to ψ is updated with one-step propagation. As every reservoir neuron is projecting to the readout units, every reservoir connection is updated with the learning rule.
W n−1 kl contributes to n i for all i, k and l through x n k . We can now utilize the chain rule, Partial derivatives with respect to state variables can be computed by combining Eqs. (2), (3), and (4).
Again following Eq. (11), where a n k is the k − th element of the vector W inp u n + W n−1 x n−1 and f k (a n k ) is the k-th element of the resulting state vector. We define α k as the inverse of the time constant for reservoir unit k, 0 < α k = 1/τ k < 1. With Eqs. (9), (10), (12), (13) and (14), the update rule in Eq. (8) is completed.

Adapting neuronal time constants
Gradient descent is also used for adapting the time constants of the reservoir neurons. Note that in this case time constant is not a global parameter anymore. Let to have a continuity over all α n j and to keep them in the desired range.
The complexity of the readout-only learning is O(L N) for each time step, where L is the number of readout units. Running the trained ESN on the test input stream has also a complexity of O(L N) for sparse reservoirs, where the number of connections linearly scales with the number of reservoir units. All the reservoirs we use belong to this sparse type. In this case, the complexity of the one-step log-likelihood propagation is also O(L N).

Existence of local maxima and confidence intervals
Due to deep architecture and reservoir transfer function, we cannot guarantee that gradient descent in the reservoir parameters yields uniquely true, i.e. globally optimal, parameters. Point process log-likelihood is not concave with respect to the whole set of parameters.
Reservoir adaptation, as a result, is not a tool for finding a globally optimal equivalent of a given BNN. The quality of reservoir adaptation is evaluated based on its improvement on the predictive performance over fixed reservoirs. For a fixed reservoir, on the other hand, learning reduces to redout-only learning, and hence, to adapting the parameters of a Generalized Linear Model (Mccullagh and Nelder 1989;Paninski 2004), which maps reservoir states and input to the conditional intensity. In this case, point-process loglikelihood is a concave function of readout parameters (see Appendix 2) and does not have non-global local maxima with respect to them (Paninski 2004). This shows that gradient descent in the readout parameters will asymptotically result in a global maximum. The absence of local minima holds for readout parameters under a fixed reservoir or if an adaptive reservoir is fixed after some epochs of training. Note that we take an online (stochastic) gradient descent approach in this work. Although online gradient descent takes a stochastic path in the parameter space for maximization of the log-likelihood, empirical evidence suggests that the average behavior of gradient descent is not affected by online updates (Bottou 2004). Upon training, the observed Fisher Information Matrix can be used to approximate the covariance matrix and confidence intervals on the readout parameters, where out is the estimation of the covariance matrix and ∇ 2 W out L(W out |N) is the Hessian of the loglikelihood. The Hessian matrix is computed at the final point estimates of the readout parameters. Confidence intervals are computed using the standard deviations obtained from the diagonal of the covariance matrix. It should be noted that we do not perform a full probabilistic learning of the parameters. We approximate the distribution on the readout parameters by the asymptotic normal distribution of the maximum likelihood estimator (Davison 2003;Pawitan 2001).

Evaluation of the learned models
We compared our results with feed-forward and recurrent architectures to the baseline method where only the readout weights were adapted. For different ESN types and architectures, we comparatively evaluated their capabilities on modeling BNN stimulus-response relations by testing the predictive power of the ESN on the observed spike trains. The continuous ESN output signal was tested for compatibility with the actual observed spikes.
Receiver Operating Characteristic (ROC) curve analysis was employed to test the quality of prediction. ROC curves are extensively used in signal detection theory and machine learning to evaluate prediction performance on binary variables. In a binary classification problem (here spike or no spike), a ROC curve is a plot of true positive rate vs. false positive rate. True positive rate is defined as Similarly, false positive rate is defined as As the ESN output is continuous, true positive rate and false positive rate will vary with respect to a moving threshold. An example ROC curve is shown in Fig. 6 (bottom). In contrast to classification accuracy, the area under a ROC curve (AUC) is more robust with respect to prior class probabilities (Bradley 1997). An AUC around 0.5 indicates a random prediction whereas an AUC for a perfect prediction will equal 1.0.

Experimental results
We employed Echo State Networks with different types and sizes of reservoirs on spike prediction tasks, where the spikes were recorded from simulations of random cortical networks and from cultured networks of cortical neurons. We investigated whether ESNs successfully predict output spikes of BNNs when they are presented with the same input streams.

Simulations of random cortical networks
We simulated 10 surrogate cortical neural networks of 1000 neurons each with 0.1 random connectivity. We used the Izhikevich neuron model and included spike-timing dependent plasticity (STDP) in our simulations (Izhikevich 2006). At each 5 ms, one of the 100 input channels was randomly selected from a uniform distribution and a pulse of 5 ms width was sent to the network. Each input channel had excitatory projections Table 1 Areas under ROC curves (in %) for the prediction of activity in 10 simulated biological neural networks with fixed recurrent (R-fxd), feed-forward adaptive (FF-adp) and recurrent adaptive (R-adp) reservoirs to 10 % of the neurons in the network. The networks were simulated for 2 h in real time with STDP and 0.5 h in real time with frozen synapses thereafter. In each network, the spikes from a randomly selected excitatory neuron were recorded with 5 ms discretization. We then trained ESNs to estimate the conditional intensity of the selected neuron's spiking process from the history of the input pulses. LI-ESNs with three different reservoir types of various sizes were used for this task, namely recurrent fixed recurrent reservoirs (R-fxd), recurrent adaptive reservoirs (R-adp) and feed-forward adaptive (FF-adp) Confidence intervals on the readout parameters of 3 reservoir types of size 100 for a selected biological neural network. The learned readout weights and 99% confidence intervals (±2.576σ ) are shown. Confidence intervals are computed by assuming a multivariate normal distribution and an infinite prior covariance matrix. Hence, inverse of the observed Fisher information matrix is used to compute confidence intervals. The weights indexed from 1 to 100 are reservoir-to-output parameters, whereas the weights indexed from 101 to 200 are direct input-to-output parameters. Greater uncertainty for input-to-output parameters are noticeable, which results from low density of the input pulses per channel reservoirs. Reservoir connectivity and each individual time constant were adapted using the method described in Section 4. For each of the 10 surrogate biological neural networks, we used reservoirs of sizes 100, 500 and 1000. A separate random sparse reservoir for each size and for each surrogate biological network was generated, where each reservoir unit was connected to 10 other units. For each random reservoir, a feedforward reservoir was generated by inverting the edges that cause recurrence as described in Fig. 4. Note that the same sized feed-forward and the recurrent reservoirs for the same biological neural network had the same topology apart from edge directions. This yielded 60 different reservoirs for the whole experiment. From each of the 10 BNN simulations we recorded data for 30 min in real time. Using a discretization time bin of 5 ms, the simulation yielded time series data for 360,000 time steps, of which 328,000 time steps for were used for training and 30,000 time steps were used for testing in each sub-experiment. For the adaptive reservoirs, the training phase included 20 full adaptation iterations of the reservoirs and 60 readout adaptation iterations. One iteration covered one forward traversal of the training data in time. No learning was performed on the test data.
Results of the ROC curve analysis show that the ESN approach performs well in conditional intensity estimation (Table 1). For all the experiments on 10 different BNNs, the estimated conditional intensity predicts output firing far better than random prediction. For BNN2, point estimates and uncertainty measures of readout parameters after readout-only-training were also presented. (Fig. 5, see Section 4.2 for computation of uncertainty measures).
We further investigated whether the predictive power of ESNs with feed-forward adaptive reservoirs readout only learning res. adaptation of leaky integrator neurons is comparable to that of recurrent reservoirs. Adaptive reservoirs are in general better than non-adaptive reservoirs. This difference, however, vanishes with increasing reservoir size (Table 1). In these experiments, feed-forward adaptive reservoirs significantly outperformed the non-adaptive recurrent reservoirs of all sizes. Furthermore, they significantly outperformed recurrent adaptive reservoirs for sizes of 100 and 500. Although adaptation brought significant improvement also to recurrent reservoirs for sizes 100 and 1000, feed-forward adaptation performed in general better than recurrent adaptation in estimating the conditional intensity. Small reservoirs (i.e. for N = 100) of the feed-forward architecture were drastically superior to the recurrent reservoirs.
One can gain some insight on the comparative performances of different reservoirs by looking at the development of the log-likelihood during training (Fig. 6). Note that we adapted the reservoir connectivity and the time constants only in the first 20 iterations. The remaining 60 iterations included only learning of readout parameters. This allowed for observing the effects of two learning phases separately. The log-likelihood clearly increases with fewer iterations in adaptive reservoirs. A sudden drop of training log-likelihood is remarkable at the transition from the full training phase to the readout training phase. Although this transition significantly reduces the fitness and might make the adaptive reservoirs worse than non-adaptive ones, previous full reservoir training pays out in the upcoming iterations.
Note that Table 1 refers only to the prediction of output spikes from input pulse trains. If the spike history of the output neuron is also taken into account, the prediction performance increases considerably. In this case, the output spike history is also fed into the ESN through an additional input channel. For BNN2, we visualize the conditional intensity estimation that takes also the output spike history into account (Fig. 7). Note the increased area under the ROC curve when compared to the second row in Table 1. Cross-correlation coefficients around 0 time lag between the conditional intensity estimation and the target spike train also improve when the output spike history is used (Fig. 8), indicating an increased similarity of estimations and observations. Cross-correlation coefficients are computed as where m is the time lag, s is the binary valued discretetime signal for the target spike train and λ(n) is the  Fig. 8 Cross-correlation coefficients between the conditional intensity estimation and the observed spike patterns for the output neuron of a selected biological neural network (BNN2) with feedforward ESN adaptation for a reservoir of 1000 neurons. Using only the input stream, AUC = 85.5% (top). Using the input stream and the spike history of the output neurons, AUC = 91.6% (bottom). Each bar shows the correlation coefficient for 1 time step of the ESN (5 ms). Around 0 time lag, the conditional intensity that uses the spike history has a higher correlation with the observed spikes. Note that there is negative correlation for some negative time lags (−90 < timelag < 0). This results from the fact that conditional intensity estimation is very low after an observed spike since the algorithm learns the negative autocorrelation of the output spikes for this time lag interval. In other words, the learned model is reflecting the inter-spike interval in its output conditional intensity estimation at time step n.s and λ stand for mean values, whereas σ s and σ λ denote standard deviations. T is the data length in time steps.
As for ROC curve areas, the peak around zero lag indicates the similarity between the conditional intensity estimation and the reference spike train.

Approximating gradient vectors with one-step propagation
In order to investigate the consequences of one-step approximation to gradient descent, we compared the gradient vectors obtained by full gradient computation and those obtained by one-step approximation. Full gradient computation would replace Eqs. (14) and (16) as respectively. The resulting learning rule would be an adaptation of the Real Time Recurrent Learning (Williams and Zipser 1989) for point processes, which is computationally very expensive for our data and reservoir sizes. It is, however, still viable to compare the gradients obtained by one-step propagation to those obtained by full gradient computation for a limited results also indicate that output spikes are encoded in reservoir activity. Point process readout learning proves more efficient in relating the reservoir activity to output spikes than CCA number of time steps. Here, we compared gradient vectors of the log likelihood with respect to all reservoir parameters for 19,000 time steps. We computed the correlation and the cosine between full and one-step gradients for each time step separately, where ∇ n−1 n and ∇ n are the gradient vectors at time step n obtained by one-step propagation and by full gradient computation, respectively; ∂ n ∂ W n−1 r and ∂ n ∂ W r are their r-th components; μ n one/ f ull and σ n one/ f ull are the mean values and the standard deviations for the components of the gradient vector; c n and n are the correlation and the angle between the full and one-step gradient vectors at time step n, respectively.
The distribution of the correlations and cosines provides insight to the consequences of one-step approximation in feed-forward and recurrent reservoirs. Over 19,000 time steps, the correlations between the one-step  (Fig. 9). Cosines constitute a very similar distribution, as the gradient vector components for each time step have approximately a mean value of 0. Strong correlations between the full and one-step gradients deliver additional explanation for the better performance adaptive reservoirs. The fact that the loglikelihood propagates to every reservoir parameter in one-step (due to full reservoir-to-output connectivity) and that one-step gradient is the dominant component of the full gradient information are possible reasons for the detected similarity. Stronger correlations for feed-forward reservoirs are in agreement with the rel-ative better performance of the feed-forward reservoir adaptation.

Information encoded in reservoir activity
A canonical correlation analysis of the reservoir activity and the output spiking activity also revealed that output spikes are encoded in the reservoir activity. Prior to readout-only learning, we recorded the reservoir states for one traversal of the training data and performed a canonical correlation analysis (CCA) (Chatfield and Collins 2000) with the observed training spikes. We then used parameters obtained by CCA for relating reservoir states to output spikes in the test data and detected maximum correlations for all of the 10 experiments in Table 1. Although not as good as in readout  learning, CCA also reveals the correlation of the reservoir activity with output spikes (Fig. 10), indicating the spike encoding in the reservoir activity.
In addition to CCA of reservoir states and output spikes, we further visualized the memory traces of the input pulses in reservoir states. We compared how past inputs are reflected in the reservoir activity for different reservoir types by considering input triggered averages (ITA). We define the input triggered average for an input-neuron pair as where m is the time lag, u n h is the input pulse from channel h at time step n (1 or 0), x n+m j is the state of the j-th neuron at time step n + m. The normalized ITA is then wherex j and σ x j are the mean value and the standard deviation for x j , respectively. Figure 11 depicts the ITA hj curves for 800 input neuron pairs in 3 different reservoir types for a selected biological neural network. Memory traces for adaptive reservoirs, especially for Fig. 14 Cross-correlation coefficients between the conditional intensity estimation and the observed output events for a selected (1) active output electrode without using the output spike history (res. size of 500). Coefficients are shown for fixed recurrent (a), feed-forward adaptive (b), and recurrent adaptive reservoirs (c). Each bar shows the correlation coefficient for 5 time steps of the ESN (5 ms). As a baseline a prediction based on event rate (without an ESN, conditional intensity proportional to event rate) is also shown (d). The recurrent adaptive reservoir performs better here, considering the height of the peak and the width of the cross-correlogram  Results are aggregated over 10 random reservoirs with size 500 for all active output electrodes, for the best active electrode and for the worst active electrode recurrent adaptive reservoirs, are more scattered over time. This increased heterogeneity in the responses is very likely a reflection of time constants in adaptive reservoirs.

Prediction of spontaneous events in neuronal cultures
To test our approach on living neural networks, we aimed at predicting spontaneous events in neuronal cultures. We defined an event as a group of spikes recorded by the same MEA electrode, whose maximum inter spike interval is less than 60 ms. The time of the first spike in the group is regarded as the event time. It should be noted again that the activity in neuronal cultures is typically a sequence of network bursts. If there were at least 100 ms between two successive events, we regarded them as parts of different bursts. These criteria mostly excluded isolated spikes from network bursts and clustered fast repetitive firings at an electrode, e.g. a cellular burst, into a single event. We recorded spikes from a neuronal cell culture for 30 min, detecting bursts and events according to the above definition. We used the data for training except for the last 200 s, which were reserved for testing. We selected approximately 3/4th of the 60 MEA electrodes and treated their events as the input stream; and the remaining electrodes as output. This selection was based on spatial order of the electrodes, i.e. input and output electrodes occupied two distinct areas in the culture. If an electrode had never recorded spikes in the training event train, it was regarded as inactive and was excluded from the learning and testing processes. The evaluation task was to estimate the conditional intensity of the output events for each output electrode at time step n using the total input event stream until time step n (1 ms bin size). Note that only very few events occur outside the bursts. Therefore, the algorithm performs learning and prediction only during bursts. Based on estimates for all time steps during network bursts in the last 200 s, we evaluated the predictive performance of the learned ESN models using ROC curve analysis on the target events and estimated conditional intensity. We selected the electrodes that recorded at least 15 events within the evaluation window of 200 s. Figure 12 (left) shows the AUC for each electrode using reservoirs of different types all with 500 neurons. Note that the prediction results in Fig. 12 is based only on input event streams, i.e. the event history of the output electrodes is excluded. Except for output electrode 3, all AUC measures are far above the random level. This supports the notion that most output neurons in the culture do not fire randomly but they rather follow patterns (Shahaf et al. 2008). Even electrode 3 is above the chance level but it has low predictability compared to other electrodes. Aggregated test data AUC over active electrodes (except electrode 3) for different ratios of training data is shown in Fig. 12 (right). As expected, AUC improved with more training data. The  improvement, however, was not dramatic. For a selected output electrode the details of the ROC analysis is illustrated in Fig. 13. Table 2 shows aggregated areas under ROC curves over all active electrodes except electrode 3 for a selected culture (Culture 1), without (top) and with (middle) using the output event history. Each entry in the table was obtained with a single reservoir of the corresponding size. As expected, the prediction performance usually increased with the reservoir size. Using the event history of the output electrodes did not improve the prediction performance, thus, no additional information could be found in the output event history. Note that utilizing the event history results in an approximately 25% increase in the input dimensionality. Possible redundancy of the event history with the already existing input might have caused a suboptimal performance of the learning algorithm, decreasing the prediction accuracy in several reservoirs. Feed-forward adaptive and recurrent adaptive reservoirs performed better than the fixed reservoirs for most reservoir sizes.
As a baseline for comparison purposes, we used the input event rate (regardless of the spatial information, i.e. without using the electrode number) for the conditional intensity estimation (Table 2, below). Although this method proved to have some predictive power especially for a rate kernel of 20 ms, it is obviously not as good as ESN predictors. Cross-correlation analysis reveals further information about comparative performances of different methods (Fig. 14). The higher peak and the smaller width of the cross-correlogram is an indicator for the better performance of recurrent adaptive reservoirs. To further compare different reservoir types, we experimented with 3 different cultures and employed 10 randomly generated reservoirs of size 500. We used the acyclic forms of the same reservoirs for feed-forward adaptation. In 2 of these 3 cultures (1 and 3), the recurrent adaptive reservoirs had the best average predictive power (Table 3). In the same 2 cultures, the prediction of the recurrent adaptive reservoirs significantly outperformed the others for several output electrodes (Table 4).

Fig. 15
Estimated conditional intensity for a selected output electrode (5) of the neuronal culture in the next-event-prediction task using a feed-forward adaptive reservoir with 500 neurons.  Results are aggregated over 10 random reservoirs with size 300 for all active output electrodes, for the best active electrode and for the worst active electrode 6.3 Next-event prediction in neuronal cultures Recent findings indicate that the temporal order of the events in neuronal cultures carries significant information about the type of stimuli (Shahaf et al. 2008). In this experimental setting, we investigated whether the ESNs can model the structure in the temporal order of the events. We deleted timing information of the events from the previous setting and, as a result, we obtained data of temporally ordered events. More precisely, we used operational time steps for ESN analysis, i.e. each event appeared only in 1 time step and each time step contained only 1 event.
In this setting, we treated every electrode as output and used the whole event history for prediction. The task was again to estimate the conditional intensity of the output events for each electrode at time step n using the total event stream until time step n. Note that point process modeling is not necessarily the optimal tool for this data, where each abstract time step can contain only 1 event. The general framework taken in the current work, however, can still be employed for predictive modeling of temporal event orders. Again, ROC curve analysis indicates a good average prediction quality with larger reservoirs, feed-forward adaptive reservoirs outperforming the recurrent reservoirs up to sizes of 500 (Table 5, Fig. 15).
To further compare different reservoir types, 10 random reservoirs of 300 neurons together with their acyclic forms were tested on 3 different cultures ( Table 6). The predictability varied with respect to the electrode that recorded the event, with AUC ranging from 50.5% to 93.4%. In pairwise comparison of the 3 reservoir types, FF-adp reservoirs significantly (signedrank test, p < 0.05) outperformed other methods for many electrodes (Table 7). Conversely, the number of electrodes, for which the other methods performed better, was much smaller.

Conclusion
The implications of our results are manifold. Firstly, our results indicate that reservoir computing is a potential candidate for modeling neural activity including neural encoding and decoding. With their expressive power for different activity measures (e.g. spike rates, correlations etc.), reservoir computing tools might help for analysis of neural data. In our experiments, ESNs of leaky integrator neurons proved successful for modeling response-dynamics (e.g. stimulus-response relations and spatio-temporal dynamics) of simulated and cultured biological neural neural networks.
On the methodological side, we showed that ESN learning algorithms can be adapted for event data, such as spikes or spike groups, using a point process framework. We proposed a reservoir adaptation method for event data, which can be used to adapt connectivity and time constants of the reservoir neurons. The experimental results indicate that reservoir adaptation can significantly improve the ESN performance over readout-only training. We utilized feed-forward networks with leakyintegrator neurons as reservoirs with a comparable predictive power to recurrent reservoirs when trained with log-likelihood propagation. In modeling stimulusresponse relations of simulated BNNs, feed-forward reservoir adaptation outperformed other methods up to 500 neurons. This outperformance was statistically significant. For the event timing prediction task in neuronal cultures, adaptive recurrent reservoir adaptation outperformed the other methods (in 2 of 3 cultures), whereas feed-forward adaptation were better in the next-event type prediction task in all 3 cultures. This might indicate that the type of encoding in neural systems (order or timing) favors different architectures for decoding. An analysis of the structure-coding relations, however, is beyond the scope of this note. Why feed-forward reservoir adaptation can work better than recurrent reservoir adaptation necessitates also more theoretical analysis. We manually optimized global parameters (spectral radius, learning rates and A) for recurrent fixed reservoirs. Recurrent adaptive reservoirs started from these values. Although we can think of no obvious disadvantage for recurrent reservoirs, feedforward adaptation outperformed recurrent adaptation in some tasks. We experimentally showed that one-step propagation approximates the gradients better in feedforward reservoirs than in recurrent ones. In our opinion, better structuring of the reservoir parameters, i.e. separation of the memory parameters from the reservoir connectivity, might further underlie the relative better performance of feed-forward adaptation. For instance, a small connectivity change in the recurrent adaptation might have a more dramatic effect on the reservoir memory than in case of feed-forward adaptation. Our findings on feed-forward networks are also in accordance with the recent work by Ganguli et al. (2008), Goldman (2009), Murphy and Miller (2009), who show that stable fading memory can be realized by feed-forward or functionally feed-forward networks and that feed-forward networks have desirable properties in terms of trade-off between noise amplification and memory capacity.
Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.

Appendix 1: Implementation details
For all the experimental settings, we generated recurrent sparse random reservoirs. Regardless of the reservoir size, each reservoir unit was connected to 10 other units. The initial connection strengths were drawn from a uniform distribution over the interval [−0.5 0.5]. The entries of the input weight matrix W in were drawn from a uniform distribution over [−1 1]. We normalized the connectivity so that the maximum absolute eigenvalue, i.e. the spectral radius, was 1, with which the ESNs were experimentally detected to work well on all experiments. We additionally generated the feed-forward form of the each recurrent reservoir as described in Fig. 4. After changing from recurrent to feedforward architecture, we did not manually modify the connectivity parameters. The leakage time constants, τ j = 1/α j , were initialized according to where α j were uniformly drawn from [−1.5 1.5].
Note that the initial value of the spectral radius and the initial range for time constants were manually arranged to work well with recurrent fixed reservoirs. Adaptive reservoirs used this fixed settings as starting parameters.
For fixed reservoirs, we used 80 epochs of readoutonly training. Initially, we set η out = 0.7 and gradually decreased the output learning rate with where E was the epoch number. This was done for fast learning and convergence of the algorithm. For adaptive reservoirs, we used 20 epochs for reservoir and readout training and 60 epochs for readoutonly training. The readout-only learning employed the same time dependent η out as the fixed reservoirs. For reservoir training, we set the initial values η out = η res = η α = η = 0.2 and decreased it with η ← 0.5 η, when the increase in the training log-likelihood dropped below 0.0003 per output channel and time step.
The parameter A in f out (ξ ) = exp(A ξ) was set to 0.2, which was manually found to work well for the selected learning rates.
The ESN learning algorithm described in this note was implemented in MATLAB and C++ languages.
The results reported here were produced by the C++ implementation. The source code is available under http://www.informatik.uni-freiburg.de/∼guerel/PP_ ESN.zip (currently only MATLAB).


Incremental HMM Alignment for MT System Combination

Inspired by the incremental TER alignment, we re-designed the Indirect HMM (IHMM) alignment, which is one of the best hypothesis alignment methods for conventional MT system combination, in an incremental manner. One crucial problem of incremental alignment is to align a hypothesis to a confusion network (CN). Our incremental IHMM alignment is implemented in three different ways: 1) treat CN spans as HMM states and define state transition as distortion over covered n-grams between two spans; 2) treat CN spans as HMM states and define state transition as distortion over words in component translations in the CN; and 3) use a consensus decoding algorithm over one hypothesis and multiple IHMMs, each of which corresponds to a component translation in the CN. All these three approaches of incremental alignment based on IHMM are shown to be superior to both incremental TER alignment and conventional IHMM alignment in the setting of the Chinese-to-English track of the 2008 NIST Open MT evaluation.

Introduction
Word-level combination using confusion network (Matusov et al. (2006) and Rosti et al. (2007)) is a widely adopted approach for combining Machine Translation (MT) systems' output. Word alignment between a backbone (or skeleton) translation and a hypothesis translation is a key problem in this approach. Translation Edit Rate (TER, Snover et al. (2006)) based alignment proposed in Sim et al. (2007) is often taken as the baseline, and a couple of other approaches, such as the Indirect Hidden Markov Model (IHMM, He et al. (2008)) and the ITG-based alignment (Karakos et al. (2008)), were recently proposed with better results reported. With an alignment method, each hypothesis is aligned against the backbone and all the alignments are then used to build a confusion network (CN) for generating a better translation.
However, as pointed out by Rosti et al. (2008), such a pair-wise alignment strategy will produce a low-quality CN if there are errors in the alignment of any of the hypotheses, no matter how good the alignments of other hypotheses are. For example, suppose we have the backbone "he buys a computer" and two hypotheses "he bought a laptop computer" and "he buys a laptop". It will be natural for most alignment methods to produce the alignments in Figure 1a. The alignment of hypothesis 2 against the backbone cannot be considered an error if we consider only these two translations; nevertheless, when added with the alignment of another hypothesis, it produces the low-quality CN in Figure 1b, which may generate poor translations like "he bought a laptop laptop". While it could be argued that such poor translations are unlikely to be selected due to language model, this CN does disperse the votes to the word "laptop" to two distinct arcs. Rosti et al. (2008) showed that this problem can be rectified by incremental alignment. If hypothesis 1 is first aligned against the backbone, the CN thus produced (depicted in Figure 2a) is then aligned to hypothesis 2, giving rise to the good CN as depicted in Figure 2b. 1 On the other hand, the Figure 1: An example bad confusion network due to pair-wise alignment strategy correct result depends on the order of hypotheses. If hypothesis 2 is aligned before hypothesis 1, the final CN will not be good. Therefore, the observation in Rosti et al. (2008) that different order of hypotheses does not affect translation quality is counter-intuitive. This paper attempts to answer two questions: 1) as incremental TER alignment gives better performance than pair-wise TER alignment, would the incremental strategy still be better than the pairwise strategy if the TER method is replaced by another alignment method? 2) how does translation quality vary for different orders of hypotheses being incrementally added into a CN? For question 1, we will focus on the IHMM alignment method and propose three different ways of implementing incremental IHMM alignment. Our experiments will also try several orders of hypotheses in response to question 2. This paper is structured as follows. After setting the notations on CN in section 2, we will first introduce, in section 3, two variations of the basic incremental IHMM model (IncIHMM1 and IncIHMM2). In section 4, a consensus decoding algorithm (CD-IHMM) is proposed as an alternative way to search for the optimal alignment. The issues of alignment normalization and the order of hypotheses being added into a CN are discussed in sections 5 and 6 respectively. Experiment results and analysis are presented in section 7.

Preliminaries: Notation on Confusion Network
Before the elaboration of the models, let us first clarify the notation on CN. A CN is usually described as a finite state graph with many spans. Each span corresponds to a word position and contains several arcs, each of which represents an alternative word (could be the empty symbol , ) at that position. Each arc is also associated with M weights in an M -way system combination task. Follow Rosti et al. (2007), the i-th weight of an arc is r 1 1+r , where r is the rank of the hypothesis in the i-th system that votes for the word represented by the arc. This conception of CN is called the conventional or compact form of CN. The networks in Figures 1b and 2b are examples.
On the other hand, as a CN is an integration of the skeleton and all hypotheses, it can be conceived as a list of the component translations. For example, the CN in Figure 2b can be converted to the form in Figure 3. In such an expanded or tabular form, each row represents a component translation. Each column, which is equivalent to a span in the compact form, comprises the alternative words at a word position. Thus each cell represents an alternative word at certain word position voted by certain translation. Each row is assigned the weight 1 1+r , where r is the rank of the translation of some MT system. It is assumed that all MT systems are weighted equally and thus the Figure 3: An example of confusion network in tabular form rank-based weights from different system can be compared to each other without adjustment. The weight of a cell is the same as the weight of the corresponding row. In this paper the elaboration of the incremental IHMM models is based on such tabular form of CN.
Let E I 1 = (E 1 . . . E I ) denote the backbone CN, and e J 1 = (e 1 . . . e J ) denote a hypothesis being aligned to the backbone. Each e j is simply a word in the target language. However, each E i is a span, or a column, of the CN. We will also use E(k) to denote the k-th row of the tabular form CN, and E i (k) to denote the cell at the k-th row and the i-th column. W (k) is the weight for E(k), and Note that E(k) contains the same bag-of-words as the k-th original translation, but may have different word order. Note also that E(k) represents a word sequence with inserted empty symbols; the sequence with all inserted symbols removed is known as the compact form of E(k).

The Basic IncIHMM Model
A naïve application of the incremental strategy to IHMM is to treat a span in the CN as an HMM state. Like He et al. (2008), the conditional probability of the hypothesis given the backbone CN can be decomposed into similarity model and distortion model in accordance with equation 1 The similarity between a hypothesis word e j and a span E i is simply a weighted sum of the similarities between e j and each word contained in E i as equation 2: The similarity between two words is estimated in exactly the same way as in conventional IHMM alignment.
As to the distortion model, the incremental IHMM model also groups distortion parameters into a few 'buckets': The problem in incremental IHMM is when to apply a bucket. In conventional IHMM, the transition from state i to j has probability: It is tempting to apply the same formula to the transitions in incremental IHMM. However, the backbone in the incremental IHMM has a special property that it is gradually expanding due to the insertion operator. For example, initially the backbone CN contains the option e i in the i-th span and the option e i+1 in the (i+1)-th span. After the first round alignment, perhaps e i is aligned to the hypothesis word e j , e i+1 to e j+2 , and the hypothesis word e j+1 is left unaligned. Then the consequent CN have an extra span containing the option e j+1 inserted between the i-th and (i + 1)-th spans of the initial CN. If the distortion buckets are applied as in equation 3, then in the first round alignment, the transition from the span containing e i to that containing e i+1 is based on the bucket c(1), but in the second round alignment, the same transition will be based on the bucket c(2). It is therefore not reasonable to apply equation 3 to such gradually extending backbone as the monotonic alignment assumption behind the equation no longer holds. There are two possible ways to tackle this problem. The first solution estimates the transition probability as a weighted average of different distortion probabilities, whereas the second solution converts the distortion over spans to the distortion over the words in each hypothesis E(k) in the CN.

Distortion Model 1: simple weighting of covered n-grams
Distortion Model 1 shifts the monotonic alignment assumption from spans of CN to n-grams covered by state transitions. Let us illustrate this point with the following examples.
In conventional IHMM, the distortion probability p (i + 1|i, I) is applied to the transition from state i to i+1 given I states because such transition jumps across only one word, viz. the i-th word of the backbone. In incremental IHMM, suppose the i-th span covers two arcs e a and , with probabilities p 1 and p 2 = 1 − p 1 respectively, then the transition from state i to i + 1 jumps across one word (e a ) with probability p 1 and jumps across nothing with probability p 2 . Thus the transition probability should be p 1 · p (i + 1|i, I) + p 2 · p (i|i, I).
Suppose further that the (i + 1)-th span covers two arcs e b and , with probabilities p 3 and p 4 respectively, then the transition from state i to i + 2 covers 4 possible cases: 1. nothing ( ) with probability p 2 · p 4 ; 2. the unigram e a with probability p 1 · p 4 ; 3. the unigram e b with probability p 2 · p 3 ; 4. the bigram e a e b with probability p 1 · p 3 . Accordingly the transition probability should be The estimation of transition probability can be generalized to any transition from i to i by expanding all possible n-grams covered by the transition and calculating the corresponding probabilities. We enumerate all possible cell sequences S(i, i ) covered by the transition from span i to i ; each sequence is assigned the probability where the cell at the i -th span is on some row E(k). Since a cell may represent an empty word, a cell sequence may represent an n-gram where We denote |S(i, i )| to be the length of n-gram represented by a particular cell sequence S(i, i ). All the cell sequences S(i, i ) can be classified, with respect to the length of corresponding n-grams, into a set of parameters where each element (with a particular value of n) has the probability The probability of the transition from i to i is: That is, the transition probability of incremental IHMM is a weighted sum of probabilities of 'ngram jumping', defined as conventional IHMM distortion probabilities. However, in practice it is not feasible to expand all possible n-grams covered by any transition since the number of n-grams grows exponentially. Therefore a length limit L is imposed such that for all state transitions where |i − i| ≤ L, the transition probability is calculated as equation 4, otherwise it is calculated by: for some q between i and i . In other words, the probability of longer state transition is estimated in terms of the probabilities of transitions shorter or equal to the length limit. 2 All the state transitions can be calculated efficiently by dynamic programming.
A fixed value P 0 is assigned to transitions to null state, which can be optimized on held-out data. The overall distortion model is:

Distortion Model 2: weighting of distortions of component translations
The cause of the problem of distortion over CN spans is the gradual extension of CN due to the inserted empty words. Therefore, the problem will disappear if the inserted empty words are removed. The rationale of Distortion Model 2 is that the distortion model is defined over the actual word sequence in each component translation E(k). Distortion Model 2 implements a CN in such a way that the real position of the i-th word of the kth component translation can always be retrieved. The real position of E i (k), δ(i, k), refers to the position of the word represented by E i (k) in the compact form of E(k) (i.e. the form without any inserted empty words), or, if E i (k) represents an empty word, the position of the nearest preceding non-empty word. For convenience, we also denote by δ (i, k) the null state associated with the state of the real word δ (i, k). Similarly, the real length of E(k), L(k), refers to the number of non-empty words of E(k).
The transition from span i to i is then defined as where k is the row index of the tabular form CN.
Depending on E i (k) and E i (k), p k (i|i ) is computed as follows: where p refers to the conventional IHMM distortion probability as defined by equation 3.
2. if E i (k) represents a real word but E i (k) the empty word, then Like conventional HMM-based word alignment, the probability of the transition from a null state to a real word state is the same as that of the transition from the real word state associated with that null state to the other real word state. Therefore, 3. if E i (k) represents the empty word but E i (k) a real word, then The second option is due to the constraint that a null state is accessible only to itself or the real word state associated with it. Therefore, the transition from i to i is in fact composed of the first transition from i to δ(i, k) and the second transition from δ(i, k) to the null state at i.

if both E i (k)
and E i (k) represent the empty word, then, with similar logic as cases 2 and 3,

Incremental Alignment using Consensus Decoding over Multiple IHMMs
The previous section describes an incremental IHMM model in which the state space is based on the CN taken as a whole. An alternative approach is to conceive the rows (component translations) in the CN as individuals, and transforms the alignment of a hypothesis against an entire network to that against the individual translations. Each individual translation constitutes an IHMM and the optimal alignment is obtained from consensus decoding over these multiple IHMMs. Alignment over multiple sequential patterns has been investigated in different contexts. For example, Nair and Sreenivas (2007) proposed multipattern dynamic time warping (MPDTW) to align multiple speech utterances to each other. However, these methods usually assume that the alignment is monotonic. In this section, a consensus decoding algorithm that searches for the optimal (non-monotonic) alignment between a hypothesis and a set of translations in a CN (which are already aligned to each other) is developed as follows.
A prerequisite of the algorithm is a function for converting a span index to the corresponding HMM state index of a component translation. The two functions δ and δ s defined in section 3.2 are used to define a new function: Accordingly, given the alignment a J 1 = a 1 . . . a J of a hypothesis (with J words) against a CN (where each a j is an index referring to the span of the CN), we can obtain the alignmentã k = δ(a 1 , k) . . .δ(a J , k) between the hypothesis and the k-th row of the tabular CN. The real length function L(k) is also used to obtain the number of non-empty words of E(k).
Given the k-th row of a CN, E(k), an IHMM λ(k) is formed and the cost of the pair-wise alignment,ã k , between a hypothesis h and λ(k) is defined as: The cost of the alignment of h against a CN is then defined as the weighted sum of the costs of the K alignmentsã k : where Λ = {λ(k)} is the set of pair-wise IHMMs, and W (k) is the weight of the k-th row. The optimal alignmentâ is the one that minimizes this cost: A Viterbi-like dynamic programming algorithm can be developed to search forâ by treating CN spans as HMM states, with a pseudo emission probability as and a pseudo transition probability as Note that P (e j |E a j ) and P (j|i) are not true probabilities and do not have the sum-to-one property.

Alignment Normalization
After alignment, the backbone CN and the hypothesis can be combined to form an even larger CN. The same principles and heuristics for the construction of CN in conventional system combination approaches can be applied. Our incremental alignment approaches adopt the same heuristics for alignment normalization stated in He et al. (2008). There is one exception, though. All 1-N mappings are not converted to N − 1 -1 mappings since this conversion leads to N − 1 inser-tion in the CN and therefore extending the network to an unreasonable length. The Viterbi alignment is abandoned if it contains an 1-N mapping. The best alignment which contains no 1-N mapping is searched in the N-Best alignments in a way inspired by Nilsson and Goldberger (2001). For example, if both hypothesis words e 1 and e 2 are aligned to the same backbone span E 1 , then all alignments a j={1,2} = i (where i = 1) will be examined. The alignment leading to the least reduction of Viterbi probability when replacing the alignment a j={1,2} = 1 will be selected.

Order of Hypotheses
The default order of hypotheses in Rosti et al. (2008) is to rank the hypotheses in descending of their TER scores against the backbone. This paper attempts several other orders. The first one is system-based order, i.e. assume an arbitrary order of the MT systems and feeds all the translations (in their original order) from a system before the translations from the next system. The rationale behind the system-based order is that the translations from the same system are much more similar to each other than to the translations from other systems, and it might be better to build CN by incorporating similar translations first. The second one is N-best rank-based order, which means, rather than keeping the translations from the same system as a block, we feed the top-1 translations from all systems in some order of systems, and then the second best translations from all systems, and so on. The presumption of the rank-based order is that top-ranked hypotheses are more reliable and it seemed beneficial to incorporate more reliable hypotheses as early as possible. These two kinds of order of hypotheses involve a certain degree of randomness as the order of systems is arbitrary. Such randomness can be removed by imposing a Bayes Risk order on MT systems, i.e. arrange the MT systems in ascending order of the Bayes Risk of their top-1 translations. These four orders of hypotheses are summarized in Table 1. We also tried some intuitively bad orders of hypotheses, including the reversal of these four orders and the random order.

Evaluation
The   (2008)). In the following sections, the incremental IHMM approaches using distortion model 1 and 2 are named as IncIHMM1 and IncIHMM2 respectively, and the consensus decoding of multiple IHMMs as CD-IHMM. The baselines include the TER-based method in Rosti et al. (2007), the incremental TER method in Rosti et al. (2008), and the IHMM approach in He et al. (2008). The development (dev) set comprises the newswire and newsgroup sections of MT06, whereas the test set is the entire MT08. The 10best translations for every source sentence in the dev and test sets are collected from eight MT systems. Case-insensitive BLEU-4, presented in percentage, is used as evaluation metric.
The various parameters in the IHMM model are set as the optimal values found in He et al. (2008). The lexical translation probabilities used in the semantic similarity model are estimated from a small portion (FBIS + GALE) of the constrained track training data, using standard HMM alignment model (Och and Ney (2003)). The backbone of CN is selected by MBR. The loss function used for TER-based approaches is TER and that for IHMM-based approaches is BLEU. As to the incremental systems, the default order of hypotheses is the ascending order of TER score against the backbone, which is the order proposed in Rosti et al. (2008). The default order of hypotheses for our three incremental IHMM approaches is N-best rank order with Bayes Risk system order, which is empirically found to be giving the highest BLEU score. Once the CN is built, the final system combination output can be obtained by decoding it with a set of features and decoding parameters. The features we used include word confidences, language model score, word penalty and empty word penalty. The decoding parameters are trained by maximum BLEU training on the dev set. The training and decoding processes are the same as described by Rosti et al. (2007).   Table 2 lists the BLEU scores achieved by the three baseline combination methods and IncIHMM2. The comparison between pairwise and incremental TER methods justifies the superiority of the incremental strategy. However, the benefit of incremental TER over pair-wise TER is smaller than that mentioned in Rosti et al. (2008), which may be because of the difference between test sets and other experimental conditions. The comparison between the two pair-wise alignment methods shows that IHMM gives a 0.7 BLEU point gain over TER, which is a bit smaller than the difference reported in He et al. (2008). The possible causes of such discrepancy include the different dev set and the smaller training set for estimating semantic similarity parameters. Despite that, the pair-wise IHMM method is still a strong baseline. Table 2 also shows the performance of IncIHMM2, our best incremental IHMM approach. It is almost one BLEU point higher than the pair-wise IHMM baseline and much higher than the two TER baselines.  Table 3: Comparison between the three incremental IHMM approaches els is to shift the distortion over spans to the distortion over word sequences. In distortion model 2 the word sequences are those sequences available in one of the component translations in the CN. Distortion model 1 is more encompassing as it also considers the word sequences which are combined from subsequences from various component translations. However, as mentioned in section 3.1, the number of sequences grows exponentially and there is therefore a limit L to the length of sequences. In general the limit L ≥ 8 would render the tuning/decoding process intolerably slow. We tried the values 5 to 8 for L and the variation of performance is less than 0.1 BLEU point. That is, distortion model 1 cannot be improved by tuning L. The similar BLEU scores as shown in Table 3 implies that the incorporation of more word sequences in distortion model 1 does not lead to extra improvement.

Comparison among the Incremental IHMM Models
Although consensus decoding is conceptually different from both variations of IncIHMM, it can indeed be transformed into a form similar to IncIHMM2. IncIHMM2 calculates the parameters of the IHMM as a weighted sum of various probabilities of the component translations. In contrast, the equations in section 4 shows that CD-IHMM calculates the weighted sum of the logarithm of those probabilities of the component translations. In other words, IncIHMM2 makes use of the sum of probabilities whereas CD-IHMM makes use of the product of probabilities. The experiment results indicate that the interaction between the weights and the probabilities is more fragile in the product case than in the summation case. Table 4 lists the BLEU scores on the test set achieved by IncIHMM1 using different orders of hypotheses. The column 'reversal' shows the impact of deliberately bad order, viz. more than one BLEU point lower than the best order. The random order is a baseline for not caring about order of hypotheses at all, which is about 0.7 BLEU  Table 4: Comparison between various orders of hypotheses. 'System' means system-based order; 'Rank' means N-best rank-based order; 'BR' means Bayes Risk order of systems. The numbers are the BLEU scores on the test set.

Impact of Order of Hypotheses
point lower than the best order. Among the orders with good performance, it is observed that N-best rank order leads to about 0.2 to 0.3 BLEU point improvement, and that the Bayes Risk order of systems does not improve performance very much. In sum, the performance of incremental alignment is sensitive to the order of hypotheses, and the optimal order is defined in terms of the rank of each hypothesis on some system's n-best list.

Conclusions
This paper investigates the application of the incremental strategy to IHMM, one of the state-ofthe-art alignment methods for MT output combination. Such a task is subject to the problem of how to define state transitions on a gradually expanding CN. We proposed three different solutions, which share the principle that transition over CN spans must be converted to the transition over word sequences provided by the component translations. While the consensus decoding approach does not improve performance much, the two distortion models for incremental IHMM (IncIHMM1 and IncIHMM2) give superb performance in comparison with pair-wise TER, pair-wise IHMM, and incremental TER. We also showed that the order of hypotheses is important as a deliberately bad order would reduce translation quality by one BLEU point.