SSA-Caterpillar in Group Anonymity

Nowadays, it is a common practice to protect various types of statistical data before publishing them for different researches. For instance, when conducting extensive demographic surveys such as national census, the collected data should be at least depersonalized to guarantee proper level of privacy preservation. In practice, even more complicated methods of data protection need to be used. All these methods can be generally divided into two classes. The first ones aim at providing individual data anonymity, whereas the other ones are focused on protecting information about a group of respondents. In this paper, we propose a novel technique of providing group anonymity in statistical data using singular spectrum analysis (SSA).Also, we apply SSA to defining hidden patterns in demographic data distribution.

I. INTRODUCTION
"He who controls the information, controls the world." This famous quote vividly expresses one of the main features of a post-industrial world. Today, necessity of collecting, storing, and especially processing different kinds of information is understood well by those who want to keep up with the pace of contemporary society.
Without a doubt, national population census is the most extensive activity of collecting large amounts of data. The aggregated census data often are accessible in suitable forms like statistical tables or graphs, and on suitable convenient media (paper publications or OLAP databases). Also, they are disseminated by means of suitable channels such as publishers and the Internet [1]. At the same time, big samples of the primary census data usually are reduced to one microfile containing gathered information about respondents, possibly, in a coded form (see Table 1). Obviously, necessary precautions should be taken beforehand in order to guarantee acceptable level of data privacy.  By data privacy we understand a user protection against discovery and further misuse of identity by other users [2, p. 69]. According to ISO/IEC 15408-2:2008, privacy consists of four family functional requirements, anonymity being the principal one.
In the following subsection, we will cover some approaches to providing anonymity in microfiles.

A. Anonymity
In a consolidated proposal for terminology [3], anonymity of a subject means that the subject is not identifiable (uniquely characterized) within a set of subjects. So, in this work, we will treat both individual and group anonymity as a particular subject's property of being unidentifiable within a microfile. But, as individual anonymity concerns individual microfile records (such as citizens, households, and enterprises), group anonymity has to do mainly with certain sets of records defined in a specific way (to be discussed later on in this paper).
All the methods of providing individual anonymity are usually referred to as part of privacy preserving data publishing (PPDP) methods. The name clearly states that some kind of modification should be applied to the data before making them widely accessed. A state-of-the-art survey of PPDP is given in [4]. Here, we will only outline the most frequently used individual anonymity methods.
So-called perturbative PPDP methods achieve data anonymity by altering dataset records in a particular way. For instance, randomization [5] implies adding noise to initial data so that it is hard to identify records with outstanding attribute values. Another approach lies in gaining k-anonymity [6], which means that at least k dataset records correspond to a specific attribute values' combination. And, of course, it is possible to swap attribute values between data records [7].
Among recently developed perturbative PPDP methods we could mention using Fourier transform [8] and singular value decomposition [9].
All the methods mentioned earlier are mainly used in noninteractive privacy mechanisms when the trusted data collector publishes a cleansed version of the collected data. The other situation arises when applying interactive privacy mechanisms when the data collector provides means for the users through which they can make queries and obtain answers about the data. The latter case is somewhat easier to deal with, mainly due to the fact that the data collector already knows what kind of utility the user wants to receive. Nevertheless, there are some specific problems to face (for instance, even if information on a particular individual isn't present in the dataset at all, privacy breach can occur provided that additional sources of information are available). To cope with such issues, the concept of differential privacy has been introduced and successfully developed (e.g., refer to [10]).
Along with perturbative methods, one might also apply various non-perturbative methods such as data recoding and data suppression [11]. These methods, unlike perturbative ones, anonymize the data without altering them.
As opposed to PPDP methods, those ones for providing group anonymity should be capable of modifying not only particular respondents' attribute values but their distribution as well. Moreover, they have to ensure that data distortion being introduced is satisfactory. That is why completely new techniques need to be developed.
At the time this paper is being written, maybe the most complete source of how to provide group anonymity is [12]. The paper proposes a generic scheme of providing group anonymity and illustrates it with some practical examples.
One of the ways to provide group anonymity described in [12] is based on normalization process presented in [13]. This makes it possible to preserve major data statistical features, e.g. mean value and standard deviation. Though, it does not seem to be quite sufficient.
Another approach is based on applying wavelet transform (WT) [14] to initial data so that it helps to conceal important data features not likely to be given away. At the same time, it enables preserving data frequency peculiarities. But, this technique has some slight disadvantages, some of which are: 1) Insufficient Utility Preserving: The way utility is being preserved isn't clear enough.
2) Strong Dependence on the Algorithm: Choosing different wavelet bases yields completely different outcomes. At the same time, there are no well-defined means of defining which base should be chosen.
3) Tight Constraints while Modifying the Data: When performing data perturbation, the researcher is highly forced to modify dataset records in a way implied by the method itself.
All these disadvantages lead to a conclusion that using WT for providing group anonymity is a very much methodoriented technique. Fortunately, there exists quiet dissimilar approach which is free of all three downsides mentioned above. Indeed, applying SSA, as we propose it in this paper, gives a researcher powerful means of modifying initial data almost deliberately and simultaneous preserving data utility.

B. Singular Spectrum Analysis
Singular spectrum analysis (SSA) (traditionally called "Caterpillar" in papers written by Russian researchers) is a relatively novel technique for analyzing time series of quite dissimilar origin. SSA-Caterpillar method finds its applications in fields ranging from climate research [15] and geophysics [16] to processing time series with missing values [17] and detecting spatial data patterns [18].
Generally speaking [17], the method's aim is to split original time series into a sum of trend, periodic, and/or noise components. What makes this technique especially powerful is that the original series do not have to be stationary, the model of the trend can be completely indeterminate, and we do not even have to know anything about periodic components within the series.
Due to main peculiarities of the method, it cannot be fully automatized. This means that the last decision is always left to the researcher. Moreover, principle parameters of the method should be also picked by a human in each and every case independently. This leads to some complications in applying the method to analyzing large amounts of information. Though, despite this fact, it is mostly efficient in performing such tasks as series smoothing, predicting, and filling in missing signal values.
SSA-Caterpillar technique was independently developed in the UK and the USA, from one side, and Russia, from the other one. At present, there exist rather thorough and complete studies of the method full of practical examples, among which [19] and [20] can be mentioned.

A. Group Anonymity
In this subsection, we will present main steps of providing group anonymity in microfiles, though we will omit some facts not necessary for the further discussion. For more complete explanation and real-life illustrations, refer to [12].
Let us propose some necessary definitions. Definition 1. An identifier is a microfile attribute which unambiguously defines a certain respondent in a microfile.
Since keeping identifiers in the microfile when publishing it is the best way to violate individual privacy, they should be always removed from the dataset. Further on, we will presume that any microfile does not contain any identifier. Definition 2. We will call certain attributes Vital attributes help define subsets of microfile respondents to be protected when providing group anonymity. Definition 3. We will call an element () Specific vital value combinations distinguish those respondent categories whose distribution is supposed to be protected.
Definition 4. We will call an element () , S being a subset of microfile data elements corresponding to the p th attribute (called a parameter set), a parameter value. The attribute itself will be called a parameter attribute.
We usually use parameter values to arrange microfile data in some particular order, so that it becomes possible to treat them as certain series.
Definition 5. A group ( , ) G V P is an attribute set consisting of several vital attributes   12 , , ..., So, providing group anonymity actually lies in performing such microfile modifications (with respect to each group ik  separately) that sensitive data become completely confided. Or, speaking more formally, we can divide providing group anonymity into the following stages: 1) Data Preparation: Construct a (depersonalized) microfile M representing statistical data to be processed.
2) Defining Goals and Targets: Define one or several groups ( , ), iii G V P 1,..., ik  standing for respondents categories to be protected.
3) Group Anonymity Calculations: For each i from 1 to k: a) Choosing data representation: Pick a proper data representation (called goal representation) c) Performing goal representation's modification: (2) (also called modifying functional). Obtain a modified goal representation.
d) Obtaining the modified microfile. Define an inverse goal mapping function (3) Obtain a modified microfile. 4) Post-Processing: Prepare the modified microfile for publishing.
We introduced some new terms in the scheme above, so we would like to clarify them a bit before moving forward. Definition 6. A goal representation  (M, G) of a dataset M with respect to a group G is a dataset reflecting specific features of a group within the microfile in a way appropriate for providing group anonymity.
In general, this dataset could be of any possible form convenient for the researcher. In this paper, we will work with one particular goal representation called a quantity signal 12 ( , ,..., ) m q q q q  . It can be constructed by counting up all the respondents in a group with a certain pair of vital value combination and a parameter value, and arranging them afterwards in any order proper for a parameter attribute.
A quantity signal provides a quantitative statistical distribution of group members from initial microfile. In some situations, extremums of this distribution could lead to disclosing restricted information, which makes it a suitable goal representation.
Definition 7. A modifying functional  : of a dataset M with respect to a group G is any function, algorithm, or procedure which transforms initial goal representation into another one ensuring that group anonymity is reached.
Frankly saying, creating modifying functional is possibly the toughest, though yet utterly important part of the whole providing group anonymity process. In this paper, we will pay most attention to one version of modifying functional based on applying SSA-Caterpillar technique.

B. SSA-Caterpillar
In the current subsection, we will discuss main steps of the classical SSA-Caterpillar algorithm. Also, we will propose some pieces of advice on how to correctly pick method's parameters to receive needed outcome. We will heavily base our explanation on [21].
So, let us be given initial non-zero time series   01 ,..., . Then, to split the series into trend, periodic, and noise components, we need to carry out following steps: 1) Embedding: Transform the series into a so-called trajectory matrix: In (4) id  , we can present decomposition of (4) the following way: In (5), Singular value decomposition (5) is optimal in the sense that it provides the best approximation of (4) among all the other matrices of rank lower than d.
From (6), it is evident that the original trajectory matrix (4) can now be viewed as the following sum: The last thing to complete is to return back from the matrix representation of our data to the series one. 4) Diagonal Averaging: Transfer each resulting matrix (6) into a time series.
To carry this operation out, let us first introduce the following notation: Finally, applying (8) to resulting matrices from (7), we obtain time series representation of our initial dataset:  I consists of only one element, then its form will strongly resemble the one of the corresponding eigenvector. This means that if any eigenvector is periodic than the corresponding time series might be periodic as well. Moreover, the series trend is highly supposed to correspond to the eigenvector with relatively low frequency.
2) Two rather close singular values could be the sign of a periodic time series component. This means that we can extract periodic components by grouping resulting matrices corresponding to such pairs of singular values.
3) The greater is the eigenvalue, the greater is the contribution of the corresponding time series component to the time series in general.
4) Relatively small eigenvalue tending to zero are likely to point at noise time series components.

5)
There is no need in choosing window length L greater than N/2. This is true because singular value decompositions of trajectory matrices obtained with window lengths L and 1 NL  are equivalent.
6) In most cases, the greater is L, the more accurate is time series decomposition. That's why it seems reasonable to take L as big as possible (and, according to the previous hint, the most appealing option is N/2). Though, in some particular cases like extracting periodic component it is advisable to pick window length being multiple of the component's period.

С. Applying SSA-Caterpillar to Group Anonymity
In this subsection, we will discuss how exactly SSA-Caterpillar technique might be applied to providing group anonymity in a microfile.
According to the group anonymity scheme presented at the beginning of this Section, the only possible stage to make use of SSA-Caterpillar is constructing the modifying functional (2). But, first we need to define what restrictions does this functional have to satisfy.
On one hand, the modifying functional has to be efficient for altering the quantity signal obtained from the original dataset. In this paper, we propose to alter the quantity signal by replacing its trend component from (9) with another one ensuring that original extremums of the signal are changed. In this case, we are able to conceal any extreme quantities from vital values' distribution which might occur to be securityintensive.
On the other hand, it is vitally necessary to make sure that data utility isn't reduced much. It is obvious that SSA-Caterpillar method is especially convenient for extracting periodic components from the original signal. Such components may be very important for statistical, demographic, and sociological researches (and this is why microfile data are used for in the first place).
So, we come to conclusion, that modifying trend component of the quantity signal with simultaneous preserving all the other components (including periodic ones) seems to be a rather adequate way of providing group anonymity in a microfile. In the next Section, we will illustrate this concept with a practical example based on real-life data.

III. PRACTICAL EXAMPLE
According to the group anonymity scheme from Section II, the first step we have to take is to build up a microfile with the data needed to be protected. We decided to work with the 5-Percent Public Use Microdata Sample File provided by the U.S. Census Bureau [22] concerning the 2000 U.S. census of population and housing microfile data in California.
The second step requires defining groups to be protected. For simplicity of calculations, we took only one group corresponding to the military personnel distributed by their age.
It is important to note here that due to the nature of the group under protection, we limited our parameter set down to values from 17 to 56 years of age. The graphical representation of the signal can be found in Fig. 4a.
As we can clearly notice, there are huge extremums at the 3 rd , the 4 th , the 5 th , and the 6 th signal positions. To provide group anonymity, we need to replace the quantity signal's trend with another one, so that these extremums are not that outstanding.
In general, we could try to transfer these extremums to some other signal elements, or make an attempt at creating new alleged maximums. But, both of these approaches either yield not very truthful data, or cannot be completed in terms of trend substituting. That is why our only option left is to somehow smooth the signal so that these extremums cannot be marked out.
At the same time, as we mentioned in the previous Section, we need to preserve periodic components of the signal to reduce loss of data utility. To be able to accomplish both tasks, we need to decompose the signal first.
For that matter, we decided to take the value of 20 as the window length (mainly because the length of the whole signal is 40). After performing quite obvious steps of constructing trajectory matrix X according to (4) and obtaining resulting matrices (5) (which we omit in this paper), we need to group the components to extract valuable periodic and trend ones.
To make correct choice, we need to follow the instructions given at the end of Subsection II-C, especially one saying about tracking down pairs of almost equal T XX eigenvalues likely to correspond to periodic signal components. That is why we present the square roots of these eigenvalues in Fig. 1a.
Unfortunately, there isn't any pair of similar eigenvalues present there. This means there is no evident periodic component within the original quantity signal. This situation is not quite preferable, because in this case there is no need in using such powerful technique as SSA-Caterpillar only to substitute the signal's trend. Though, it seems to be not very believable that the time series representing military personnel of the state of California does not contain any periodic component at all. Indeed, if we take a closer look at the data we are working with, we will find out that the main reason of such an incident is that California is the home state for San Diego military bases being the largest concentration of naval facilities in the world. So, we are strongly convinced that periodic components, if there is any, are simply suppressed by the outweighing number of young military men working at San Diego naval bases.
This signal is present in Fig. 4b. Once again, let us take a look at Fig. 1b representing T XX eigenvalues (though, now the trajectory matrix corresponds to signal q ). In this case, we can mark out two clearly visible pairs of eigenvalues likely to correspond to the periodic components. The first one is values 3 and 4, and the second one is values 5 and 6. Also, concluding from the magnitude of the first two eigenvalues, we can say that they correspond to the trend signal component. Taking this all into consideration, we can finally group all 20 signal components into the following index subsets: 1  After having obtained resulting matrices (6) and performing diagonal averaging (8), we can present our quantity signal as a sum of the following components: (1) ( In (10), (1) q stands for signal's trend, (2) q and (3) q both reflect periodic components, and (4) q will be treated as noise. Appropriate graphical representations can be seen in Fig. 2, and numerical values go as follows (we present here only three decimal numbers, though all the calculations had been carried out with much higher proximity): (1)  Judging from Fig. 2b and Fig. 2c, we can confidently admit that there are two periodic components in the quantity signal under analysis. The first one has a period of 20, and the second one seems to possess a period of 4. This is a very interesting result, though we are not qualified sociologists to make any proper conclusions about what these periods stand for. Nevertheless, such information has to be preserved for further studies.
It is important to emphasize that though signal components have slightly changed, main periodic components' features (like their period, phase and so on) have persisted. This is exactly what we understand by utility preserving.

IV. CONCLUSIONS AND FURTHER STUDIES
In this paper, we tried to attract attention to some novel techniques in the field of providing anonymity in statistical data, and also attempted to introduce new ways of preserving data utility while altering them in predefined manner.
For the first time, we proposed to apply SSA-Caterpillar to providing group anonymity in population census microfiles, and addressed its main peculiarities. At the same time, we showed that this technique can be efficiently used not only for concealing data distribution features, but also for analyzing them simultaneously. Of course, comprehensive interpreting should be done by qualified sociologists and demographers. And, last but not least, we illustrated the strong importance of correct sampling procedure. Having excluded from consideration an evident cause of incorrect data distribution, we got an opportunity of obtaining meaningful results in the output.
Still, there are some goals to be achieved in the future: 1) Modifying the Signal: There has to be developed more formal way of picking up new trend component of the signal.
2) Data Interpreting: It would be more helpful to be able to interpret data components from the sociological point of view to define those ones to be preserved in any case, and those ones to be lost without major consequences.
3) Data Adjusting: After having altered appropriate goal representations, it is necessary to construct modified microfiles. There needs to be proposed a convenient way of evaluating loss of data utility implied by such constructing procedures.


extfieldsofstudy:
Computer Science
Mathematics

text:
Modular exponentiation of matrices on FPGA-s

We describe an efficient FPGA implementation for the exponentiation of large matrices. The research is related to an algorithm for constructing uniformly distributed linear recurring sequences. The design utilizes the special properties of both the FPGA and the used matrices to achieve a very significant speedup compared to traditional architectures.

Introduction
Field-programmable gate arrays (FPGA) offer a number of special options in computation. Utilizing the unique properties of an FPGA, some algorithms that are impractical to implement on a more traditional architecture can become both convenient to create and resource-efficient. The programmable array of look-up tables commonly found on an FPGA provide both flexibility in creating logic to suit specific needs and naturally lend themselves to great parallelism in computations.
Fast operations on matrices are of great practical interest. Ways to speed up certain matrix calculations still find their way into numerous applications.
Faster implementations of matrix algorithms can be achieved either from a "software" point of view, by improving upon the algorithm itself, or from a "hardware" point of view, by using faster or differently structured architectures.
Theoretical improvements on matrix algorithms include Strassen's algorithm [12] and the Coppersmith-Winograd algorithm [2]. The naive algorithm for matrix multiplication is a well-known Θ(n 3 ) algorithm. Strassen's algorithm uses an idea similar to the Karatsuba-multiplication. It has a time complexity of O(n lg 7 ) by dividing the matrices into sub-matrices. Then by multiplying them in a different arrangement, it manages an overall lower multiplication count compared to the classical algorithm. Research implementing it on the Cell Broadband Engine can be found in [5]. Strassen's algorithm and its applicability to the project is briefly discussed in Section 7. The Coppersmith-Winograd algorithm further improves the complexity to O(n 2.376 ) by combining the idea of Strassen with the Salem-Spencer theorem. [9] discusses and compares the performance of implementations of these algorithms.
Numerous research has been done on creating efficient realizations of different matrix operations on different architectures. [8] and [10] both use FPGAs to perform matrix inversion.
The design presented here is an implementation of matrix multiplication on an FPGA. Works of similar nature can be found in [1] and [4], dealing with FPGA configurations used for floating point matrix multiplication. [11] uses an FPGA design for digital signal processing. [3] discusses another FPGA implementation for accelerating matrix multiplication.
The research in this paper is related to an algorithm for the construction of pseudo random number generators. It requires the exponentiation of large matrices to an extremely high power. This allows for numerous optimizations to be made on the FPGA implementation, resulting in an extremely fast design. A speedup factor of ∼200 is achieved compared to a highly optimized program on a more traditional architecture.
We give the details of a design implemented on a Virtex-5 XC5VLX110T FPGA that multiplies two 896 × 896 sized matrices. The matrices are defined over the mod 4 residue class ring. Using this property and the fact that the hardware uses 6-LUTs (Lookup Tables), we describe first a module that computes the dot product of vectors taken from Z 28 4 in a single clock cycle at 100MHz clock speed. With these modules we construct a matrix multiplier module that computes the C ∈ Z 20×20 4 product matrix of A ∈ Z 20×28d 4 and B ∈ Z 28d×20 4 in d clock cycles at 100MHz. The significance of the value 28 in the implementation and its experimental determination is also discussed. Finally, we describe how to use these modules for multiplying matrices taken from Z 896×896 4 . The proposed algorithm deals with the management of stored data in such a way that it can be accomplished completely in parallel with the computations. The resulting design completes the multiplication in 64800 clock cycles at 100MHz.
Future work for increasing the size of the used matrices, and further optimizing the design's performance using Strassen's algorithm is also described.

Mathematical background
The present work is initiated by a method for the construction of uniformly distributed pseudo random number generators. (See [7].) The generator uses recurring sequences modulo powers of 2 of the form u n ≡ a d−1 u n−1 + a d−2 u n−2 + · · · + a 0 u n−d mod 2 s , a i ∈ {0, 1, 2, 3}, s ∈ Z + The theoretical background can be found in [6].
The construction assumes that the values a 0 , a 1 , . . . , a d−1 are such that holds for some P(x) irreducible polynomial. It is practical to choose P(x) to have maximal order, since the order of P is closely related to the period length of the corresponding recurring sequence. The sequence u n obtained this way does not necessarily have uniform distribution, however exactly one of the following four sequences does: For the details see [7]. Finding the sequence with uniform distribution is of interest. Let  be the companion matrix of sequence u. To find which of the above sequences has a uniform distribution, we have to compute M(u) 2 d+1 −2 mod 4. If M(u) 2 d+1 −2 mod 4 equals the identity matrix, then the period length of u n is 2 d+1 − 2, which means it is not the sequence we are searching for.
The exponentiation of matrices to high powers can quickly become time consuming on traditional computers. The aim of the project was to utilize the special properties of an FPGA to achieve a significant upgrade in speed compared to implementations on more traditional architectures.

Hardware used in the implementation
The project was implemented on a Xilinx XUPV505-LX110T development platform. The board features a variety of ports for communication with the device. As a first approach the RS-232 serial port was used to send data between the board and a PC. A high-speed PCI Express connection is also available if the amount of data transferred would necessitate its use.
The board's most prominent feature is the Virtex-5 XC5VLX110T FPGA. The FPGA's main tool for computation is the array of 6-input look-up tables, arranged into 17280 Slices, with four look-up tables found in each Slice, adding up to a total of 69120 LUTs. A single 6-input LUT can store 64 bits of data, where its six input bits are used as an address to identify the single bit of data that is to be outputted. By manipulating the 64 bit content of the look-up table, it can be configured to carry out arbitrary Boolean functions with at most six input bits. In our design they are used to create LUTs performing a multiply-accumulate function, which are hierarchically arranged into larger and more complex modules. One out of four LUTs on the device can also be used as a 32 bit deep shift register; these are the basis to implement containers storing the data, which is directly fed to the computational module.
Attached to the board, there is a 256MB DDR2 SODIMM module, which is used for storing data exceeding the amount that can be practically stored on the FPGA.

Structure of modules used in the computation
The basic elements of the design are the LUTs denoted by L(a, b, s) = c, where a, b, c and s are two-digit binary numbers. The function carried out by L is a multiply-accumulate (for short: MA) function, i.e.: c ≡ (a · b) + s mod 4 .
The LUTs l 0 and l 1 were configured to the values shown in Table 1 and Table 2 to perform the multiply-accumulate function. P P P P P P P P P  Therefor m is a function that accepts a pair of vectors u, v of two-digit numbers of length n and outputs on c n−1 the two-digit dot-product of the two vectors, i.e. m(u, v) = w. P P P P P P P P P  In total, the number of LUTs used in m is 2n. Note that vectors of arbitrary length can be used in the computation if we connect the output of module m to the sum input of L[0] (c n−1 = s 0 ), and then iteratively shift u and v onto the module's input by n elements at a time: Here u ′ and v ′ are the extensions of u and v by 0's. We shall see that the number chosen for n is critical in setting many characteristics of the entire project. The experiment used for determining n will be discussed in the following chapter.
Our aim is to obtain a module that performs the matrix multiplication of A, B ∈ Z k×k 4 , where Z 4 is the mod 4 residue class ring. In the following, let C ∈ Z k×k 4 be the output matrix, such that C = A × B. Furthermore, let a i be the ith row of matrix A and let b j be the jth column of matrix B.
The multiplier units denoted by m are used to create more complex modules in a hierarchical manner. First, by taking ten m multiplier blocks we create a row of multipliers R = (m 0 , m 1 , . . . , m 9 ). This is used to compute ten consecutive elements of a single row of the output matrix: where c i,j = a i · b j . The input vector a i is used by all ten multiplier units of R. The length of these vectors, as mentioned above, can be arbitrary, but vectors of length greater than n will need to be iteratively shifted to the input of R.
By taking ten row multipliers we can create a unit M 10×10 = (R 0 , R 1 , . . . , R 9 ) which outputs a 10 × 10 sub-matrix of C: Finally, four such units are arranged so that a 20 × 20 sub-matrix of C could be obtained as output: The M 20×20 's inputs are twenty vectors from both matrices A and B. Because of hardware constraints -in particular the number of LUTs on the used device -a larger arrangement of multipliers would be impractical to implement. The module M 20×20 is comprised of 400 m multiplier units. Figure 3 shows the hierarchy of units used to build M 20×20 . The M 20×20 unit can be used iteratively to multiply matrices of arbitrary size, producing 20×20 sub-matrices of the output matrix C with each iteration. After inputting twenty rows from matrix A and twenty columns from matrix B and obtaining the desired output, we can simply repeat the process for a set of rows and columns of A and B respectively, until we obtain the entire output matrix C: Note that in the naive algorithm large matrix mult(A, B), during the main loop (lines 7-10), for each twenty rows read from A, the entire matrix B is read. During the whole procedure, matrix A will be read entirely exactly once, while matrix B will be read κ times. Methods improving on this number are described in section 6.
Since for almost all practical cases the size k of matrices A, B ∈ Z k×k 4 will be greater than the parameter n, the vectors taken from these matrices will need to be iteratively shifted onto the input of the multiplier M 20×20 , n elements at a time. Therefore, an efficient way to both store and then use the vectors taken from the matrices is the creation of FIFO type containers made of shift registers.
Let t d n be a shift register of width n and depth d. It means that t d n can store at most d vectors of length n, or equivalently a single vector of length at most nd. We choose d such that nd ≥ k, thus it can store one row or column from the input matrices A or B. Let the vector filling t d n be f = (f 0 , f 1 , . . . , f d−1 ), where f i ∈ Z n 4 , i = 0, 1, . . . , d − 1. In practice, t d n is a queue data structure. In a single step, t d n outputs a vector of length n and shifts its content by n places. For the i th activation, the container will output f i . After d activations, the container becomes empty.
One container t d n is used to store a single row or column of matrices A or B respectively. Connecting twenty of them in parallel, denoted by , we obtain a container that stores twenty rows or columns. This is exactly the amount of data the M 20×20 multiplier structure requires as input in d iteration steps. After d activations T d 20n has shifted all its stored data to M 20×20 , broken up into pieces of length n for each activation. Two such T d 20n containers are connected to M 20×20 , one for the rows taken from matrix A and one for the columns taken from matrix B. Using M 20×20 and T d 20n in a proper structure, we can execute one iteration cycle of the computation. After filling one T d 20n container with the desired twenty rows from matrix A and one T d 20n container with the desired twenty columns from matrix B, we simply send d activation signals to the containers. This will shift the data onto M 20×20 , which computes the 20 × 20 product matrix in the way described in function iterated m(u, v). The number of steps in one iteration cycle is d.

Experimental determination of parameters
Now, we turn to the determination of n (how many MA modules should be connected into a single multiplier m). This sets the length of the vectors that we use in the computation in a single step and thus has an effect on many other technical parameters of the design. The goal was to find the greatest number such that the multiplier would still reliably produce the correct dot product in a single clock cycle. Clearly, this number dependents on the used hardware and the clock frequency. For the device used, the chosen clock frequency was 100 MHz, the default frequency provided by the board.
The following experiment was devised to determine the value of n: Let S be a multiplier m, called the "Subject", and let E 0 , E 1 , . . . , E 9 be ten more m multipliers, called the "Examiners". Informally, the Examiners' duty was to verify the answers given by the Subject to questions they already knew the answer to. The "questions" here are test data: two vectors v, u of length n generated by the following sequence to obtain suitable pseudo-random values: More formally, let p ∈ Z 10 be a counter that cycles between values 0, 1, . . . , 9, incrementing its value by one with each clock cycle, and returning to value 0 after 9. For each clock cycle during the experiment, the following happens depending on the value of p: • The output of S is checked for equality with the output of E p . If inequality is detected, then an error is noted.
• The test data E p+1 is currently working on is given to S.
• New test data is given to E p−1 .
Procedure testing 1. Let S, E 0 , E 1 , . . . , E 9 be m multipliers 2. Let D be the test data generator 3. Let i ∈ N, p ∈ Z 10 4. forever do 5. i = i + 1 6. p ≡ i mod 10 7. if S out = E p out then return ERROR 8. E p−1 in ← D(i) 9. S in ← E p+1 in 10. end forever end Procedure Note that the output of S is checked every clock cycle, which yields that S has only a single cycle to calculate its answer to the question it was given Figure 5: Activity of testing module when counter's value is p in the preceding clock cycle. A given Examiner, however, has ten times more time to work on its test data. Once in every ten clock cycles, new data is given to the Examiner to work on, and its output is only checked nine clock cycles later, just before it is given new input again. This way the Examiners have enough time to compute the correct answer to the question by the time it is needed.
As the initial value for n, we have chosen 16, a number small enough to be reasonably expected to pass the criteria set for n, but large enough to be of interest. If the experiment reported no error, meaning the Subject was flawlessly able to calculate the dot product for a sufficiently long time, then the value of n was increased and the experiment repeated. After the first error was encountered, meaning the Subject was not able to keep up with the calculations, the largest value was chosen for n for which there were no errors.
On the used device, the largest such value was found to be 28 at a clock speed of 100 MHz and setting the length of m multipliers to n = 28 were able to work error-free for days without interruption.

Computation of large matrices
In the Section 4 we gave an algorithm for using the described modules for computing the product of large matrices. Following the description, the implemented design would make use of the parallelism offered by the FPGA only in the computation of dot products. Making further use of parallel operations, the design's performance can be significantly improved. In this section we describe the implementation choices made to raise the overall performance.
The biggest factor to consider is the management of data. When computing the product of large matrices, the amount of data to store and to move between the computation modules can easily exceed the size which can be practically stored on the FPGA. Fortunately, as mentioned before, a 256MB DDR2 SODIMM is connected to the board as the main data storage device. A module is generated using the Memory Interface Generator v3.5 intellectual property core provided by Xilinx to implement the logic needed to communicate with the DDR2 RAM. The module is structured hierarchically, connecting the memory device to a user interface. All communication with the device is done through two FIFO queues: one queue to send the command and address signals, while the other queue is used for write data and write data mask (when masking is allowed).
A naive utilization of the memory would be to simply read the required data before each iteration of the computation, and writing the output back after it is finished. An undesirable effect of this approach would be that the design would spend significantly more time with memory management than with the actual computation. The desired result would be that memory management (and all other auxiliary operations) were done during the time interval of the computation. Note that since both the size of the matrices and the multiplier module is fixed, the time the multiplication consumes is a fixed constant, which cannot be lowered. Optimally, the time of the computation should be an upper bound for the running time of the entire design. The difficulty of reaching this optimum lies in the high speed of the multiplier modules compared to the memory module.
One way to resolve the problem caused by slow transmission speed is to increase the amount of data stored on the FPGA. Informally, the main idea is to keep enough data in a prepared state, i.e. by the time the multiplier module finishes all of its computations, we have enough new data to continue working. More formally, let us define the following quantities: • Let d be the time necessary to complete one iteration of the computation.
As described in the previous sections, this is equal to the depth of the containers T d 20n .
• Let κ = ⌈ k 20 ⌉, where k is the size of the matrices. (A, B ∈ Z k×k 4 ) This quantity is already used in algorithm large matrix mult. For the rest of the section, it is practical to think of A and B as κ × κ sized block matrices, where each element is a 20 × 20 matrix.
• Let f(A, B) be an arbitrary algorithm executing matrix multiplication on A and B, including the memory management needed for the computation. Let K(f) be the number of times the algorithm needs to fill a T d 20n container, i.e. the number of times it has to read twenty rows or columns from the matrices. Note that completely reading either input matrices once means filling T d 20n containers κ times, since one T d 20n can store twenty rows or columns at a time. Algorithm large matrix mult's main loop (starting at line 7) reads twenty rows from matrix A (filling a T d 20n once) and reads matrix B entirely for each step. Since the loop has κ steps, it follows that K(large matrix mult) = κ 2 + κ.
• Let δ be the time it takes to fill a T d 20n container. This quantity depends on both the width and depth of the container. The total time f(A, B) spends on reading from memory to fill the containers is K(f)δ.
• Let Φ(f) be the total time the design has to spend with memory management. This is the sum of the time it spends on reading matrices A and B from the memory and the time it spends on writing the product matrix C into the memory. The number of times f has to read A and B from the memory depends on f. Note that since the size of the total output matrix C is the same as the size of A and B, writing C into the memory takes time equal to reading either matrices once from the memory. In other words, it takes κδ time. The total time the design has to spend with memory management is Φ(f) = K(f)δ + κδ.
• Let Γ (f) be the time f(A, B) spends on the computation itself. From the definition of d and κ it follows that C(large matrix mult) = dκ 2 .
The goal here is to reduce K(f) in such a way that the data required for the next iteration of the computation is always ready by the time the previous iteration ends. If this arrangement is achieved then C(f) becomes the upper bound for the running time of the design.
Storing more data on the FPGA can be done by adding more T d 20n containers to the design. During an iteration only two such containers are used directly. The rest can be used to load data necessary for the forthcoming iteration steps.
Suppose the design has z + 2 pieces of T d 20n containers. We assign z of the containers to store rows from matrix A, called "row-stores", and two of them to store columns from matrix B, called "column-stores". With this arrangement, we can carry out z − 1 iterations of the computation, using up the data stored in z − 1 row-stores and one column-store. This leaves one row-store and one column-store to load new data into during the computation. Using the above definitions, the allover computation takes (z − 1)d time. If we use all z row-stores and one column-store for the computation while the remaining column-store is devoted to loading new columns into, then we would have to load all z row-stores with new rows once we read all the columns before we can continue the computation. This would take zδ time for each case where we read all the columns but haven't read all the rows yet, which happens ⌊ κ z ⌋ times. In total, it would add ⌊ κ z ⌋zδ to the running time. Instead, the computation of the output matrix moves slightly diagonally. See Figure 7. The z − 1 row-stores used in the computations store a total of (z − 1) · 20 rows. Initially, the row-stores are filled with rows a 0 → a (z−1)·20−1 .
New rows are loaded in at a slower pace than columns are. By the time all columns are read once, the contents of the row-stores have shifted exactly to the next segment of data needed, the next (z − 1) · 20 rows. After matrix B is completely read once, the row-stores are filled with rows a (z−1)20 → a 2(z−1)·20−1 . Reading rows and columns proceeds in this manner until we've completely read matrix A once. For this reason, it is practical to choose z such that (z − 1) | κ. All together we read matrix B κ z−1 times and matrix A once. During each z−1 iterations shown in Figure 6, twenty new columns and (z−1)·20 κ new rows are loaded into the column-store and row-store currently unused by the computation. When the unused row-store is filled with twenty new rows, it becomes active, to be used in the following iterations. The row-store containing the rows with the least index becomes inactive in the computation and starts accepting the new rows read. For convenience, time quantities are measured in clock cycles at 100MHz, the clock speed of the M 20×20 multiplier.
The value of δ depends on the DDR2 RAM used. The device was used at 200MHz, and has a 64 bit wide physical data bus.
The goal of Γ (improved matrix mult) > Φ(improved matrix mult) is achieved, meaning that the running time of the design is equal to the time used by the computation. The speedup provided by the configuration can be shown by comparing its performance to a similar implementation created on a more traditional architecture. A highly optimized C++ program was created for a machine using an Intel E8400 3GHz Dual Core processor with 2GB RAM. The algorithm is strongly specialized for the task, making use of all available options for increasing performance. It uses 64 bit long variables to perform multiplication on 16 pairs of two-digit elements at once in parallel on both processor cores.
The running time of the multiplication of matrices of the same size is over 100 ms. The FPGA implementation, as mentioned above, achieves a runtime of ∼0.6 ms. On average, a speedup factor of 200 is reached using the described FPGA design.

Future work
The future course of research will focus on increasing the size of the used matrices.
As mentioned in the previous section, simply increasing the depth d of the T d 20n containers would be impractical. Since a single LUT on the device can only be configured as a 32 bit deep shift register, setting d > 32 would double the number of LUTs needed for a T d 20n , and the design is already using well over half of the device's LUTs that can be configured this way (13440 out of 17280, to be exact). Increasing the size of the matrices this way would require the restructuring of both the multiplier module and the algorithm used for memory management.
Instead, the currently implemented module can be used as a basic unit for the multiplication of larger matrices. Then the entries of the large matrices are 896 × 896 blocks.
This also allows for further optimization using Strassen's algorithm. Suppose we double the matrix sizes, interpreting them as matrices with four blocks. Using the classical algorithm, multiplying two 1792×1792 sized matrices would take eight multiplication of the blocks. Using a divide-and-conquer strategy, we can exchange one multiplication for a few extra additions. This algorithm, with its O(n lg 7 ) time complexity, could speed up the design on large matrices. We should note however, that the speed of the extra additions have to be carefully considered. Since the multiplication is already extremely fast, a similar improvement may also be necessary for additions if the overall performance upgrade is to remain significant.


Minimum Imputed-Risk: Unsupervised Discriminative Training for Machine Translation

Discriminative training for machine translation has been well studied in the recent past. A limitation of the work to date is that it relies on the availability of high-quality in-domain bilingual text for supervised training. We present an unsupervised discriminative training framework to incorporate the usually plentiful target-language monolingual data by using a rough "reverse" translation system. Intuitively, our method strives to ensure that probabilistic "round-trip" translation from a target-language sentence to the source-language and back will have low expected loss. Theoretically, this may be justified as (discriminatively) minimizing an imputed empirical risk. Empirically, we demonstrate that augmenting supervised training with unsupervised data improves translation performance over the supervised case for both IWSLT and NIST tasks.

Introduction
Missing data is a common problem in statistics when fitting the parameters θ of a model. A common strategy is to attempt to impute, or "fill in," the missing data (Little and Rubin, 1987), as typified by the EM algorithm. In this paper we develop imputation techniques when θ is to be trained discriminatively.
We focus on machine translation (MT) as our example application. A Chinese-to-English machine translation system is given a Chinese sentence x and * Zhifei Li is currently working at Google Research, and this work was done while he was a PHD student at Johns Hopkins University. asked to predict its English translation y. This system employs statistical models p θ (y | x) whose parameters θ are discriminatively trained using bilingual sentence pairs (x, y). But bilingual data for such supervised training may be relatively scarce for a particular language pair (e.g., Urdu-English), especially for some topics (e.g., technical manuals) or genres (e.g., blogs). So systems seek to exploit additional monolingual data, i.e., a corpus of English sentences y with no corresponding source-language sentences x, to improve estimation of θ. This is our missing data scenario. 1 Discriminative training of the parameters θ of p θ (y | x) using monolingual English data is a curious idea, since there is no Chinese input x to translate. We propose an unsupervised training approach, called minimum imputed risk training, which is conceptually straightforward: First guess x (probabilistically) from the observed y using a reverse Englishto-Chinese translation model p φ (x | y). Then train the discriminative Chinese-to-English model p θ (y | x) to do a good job at translating this imputed x back to y, as measured by a given performance metric. Intuitively, our method strives to ensure that probabilistic "round-trip" translation from a targetlanguage sentence to the source-language and back again will have low expected loss.
Our approach can be applied in an application scenario where we have (1) enough out-of-domain bilingual data to build two baseline translation systems, with parameters θ for the forward direction, and φ for the reverse direction; (2) a small amount of in-domain bilingual development data to discriminatively tune a small number of parameters in φ; and (3) a large amount of in-domain English monolingual data.
The novelty here is to exploit (3) to discriminatively tune the parameters θ of all translation model components, 2 p θ (y|x) and p θ (y), not merely train a generative language model p θ (y), as is the norm.
Following the theoretical development below, the empirical effectiveness of our approach is demonstrated by replacing a key supervised discriminative training step in the development of large MT systems -learning the log-linear combination of several component model scores (viewed as features) to optimize a performance metric (e.g. BLEU) on a set of (x, y) pairs -with our unsupervised discriminative training using only y. One may hence contrast our approach with the traditional supervised methods applied to the MT task such as minimum error rate training (Och, 2003;, the averaged Perceptron (Liang et al., 2006), maximum conditional likelihood (Blunsom et al., 2008), minimum risk (Smith and Eisner, 2006;Li and Eisner, 2009), and MIRA (Watanabe et al., 2007;Chiang et al., 2009).
We perform experiments using the open-source MT toolkit Joshua (Li et al., 2009a), and show that adding unsupervised data to the traditional supervised training setup improves performance.

Supervised Discriminative Training via Minimization of Empirical Risk
Let us first review discriminative training in the supervised setting-as used in MERT (Och, 2003) and subsequent work. One wishes to tune the parameters θ of some complex translation system δ θ (x). The function δ θ , which translates Chinese x to English y = δ θ (x) need not be probabilistic. For example, θ may be the parameters of a scoring function used by δ, along with pruning and decoding heuristics, for extracting a high-scoring translation of x.
The goal of discriminative training is to minimize the expected loss of δ θ (·), under a given taskspecific loss function L(y , y) that measures how bad it would be to output y when the correct output is y. For an MT system that is judged by the BLEU metric (Papineni et al., 2001), for instance, L(y , y) may be the negated BLEU score of y w.r.t. y. To be precise, the goal 3 is to find θ with low Bayes risk, where p(x, y) is the joint distribution of the inputoutput pairs. 4 The true p(x, y) is, of course, not known and, in practice, one typically minimizes empirical risk by replacing p(x, y) above with the empirical distributionp(x, y) given by a supervised training set The search for θ * typically requires the use of numerical methods and some regularization. 5 3 Unsupervised Discriminative Training with Missing Inputs

Minimization of Imputed Risk
We now turn to the unsupervised case, where we have training examples {y i } but not their corresponding inputs {x i }. We cannot compute the summand L(δ θ (x i ), y i ) for such i in (2), since δ θ (x i ) requires to know x i . So we propose to replace 3 This goal is different from the minimum risk training of Li and Eisner (2009) in a subtle but important way. In both cases, θ * minimizes risk or expected loss, but the expectation is w.r.t. different distributions: the expectation in Li and Eisner (2009) is under the conditional distribution p(y | x), while the expectation in (1) is under the joint distribution p(x, y). 4 In the terminology of statistical decision theory, p(x, y) is a distribution over states of nature. We seek a decision rule δ θ (x) that will incur low expected loss on observations x that are generated from unseen states of nature. 5 To compensate for the shortcut of using the unsmoothed empirical distribution rather than a posterior estimate of p(x, y) (Minka, 2000), it is common to add a regularization term ||θ|| 2 2 in the objective of (2). The regularization term can prevent overfitting to a training set that is not large enough to learn all parameters.
where p φ (· | ·) is a "reverse prediction model" that attempts to impute the missing x i data. We call the resulting variant of (2) the minimization of imputed empirical risk, and say that is the estimate with the minimum imputed risk 6 . The minimum imputed risk objective of (4) could be evaluated by brute force as follows.
1. For each unsupervised example y i , use the reverse prediction model p φ (· | y i ) to impute possible reverse translations X i = {x i1 , x i2 , . . .}, and add each (x ij , y i ) pair, weighted by p φ (x ij | y i ) ≤ 1, to an imputed training set .
2. Perform the supervised training of (2) on the imputed and weighted training data.
The second step means that we must use δ θ to forward-translate each imputed x ij , evaluate the loss of the translations y ij against the corresponding true translation y i , and choose the θ that minimizes the weighted sum of these losses (i.e., the empirical risk when the empirical distributionp(x, y) is derived from the imputed training set). Specific to our MT task, this tries to ensure that probabilistic "roundtrip" translation, from the target-language sentence y i to the source-language and back again, will have a low expected loss. 7 The trouble with this method is that the reverse model p φ generates a weighted lattice or hypergraph X i encoding exponentially many translations of y i , and it is computationally infeasible to forwardtranslate each x ij ∈ X i . We therefore investigate several approximations to (4) in Section 3.4. 6 One may exploit both supervised data {(xi, yi)} and unsupervised data {yj} to perform semi-supervised training via an interpolation of (2) and (4). We will do so in our experiments. 7 Our approach may be applied to other tasks as well. For example, in a speech recognition task, δ θ is a speech recognizer that produces text, whereas p φ is a speech synthesizer that must produce a distribution over audio (or at least over acoustic features or phone sequences) (Huang et al., 2010). (4) is the reverse prediction model p φ (·|·) that attempts to impute the missing x i . We will train this model in advance, doing the best job we can from available data, including any outof-domain bilingual data as well as any in-domain monolingual data 8 x.

The Reverse Prediction Model
In the MT setting, δ θ and p φ may have similar parameterization. One translates Chinese to English; the other translates English to Chinese.
Yet the setup is not quite symmetric. Whereas δ θ is a translation system that aims to produce a single, low-loss translation, the reverse version p φ is rather a probabilistic model. It is supposed to give an accurate probability distribution over possible values x ij of the missing input sentence x i . All of these values are taken into account in (4), regardless of the loss that they would incur if they were evaluated for translation quality relative to the missing x i .
Thus, φ does not need to be trained to minimize the risk itself (so there is no circularity). Ideally, it should be trained to match the underlying conditional distribution of x given y, by achieving a low conditional cross-entropy 2σ 2 φ 2 2 on some bilingual data, with the regularization coefficient σ 2 tuned on held out data.
It may be tolerable for p φ to impute mediocre translations x ij . All that is necessary is that the (forward) translations generated from the imputed x ij "simulate" the competing hypotheses that we would see when translating the correct Chinese input x i .

The Forward Translation System δ θ and
The Loss Function L(δ θ (x i ), y i ) The minimum empirical risk objective of (2) is quite general and various popular supervised training methods (Lafferty et al., 2001;Collins, 2002;Och, 2003;Crammer et al., 2006;Smith and Eisner, 2006) can be formalized in this framework by choosing different functions for δ θ and L(δ θ (x i ), y i ). The generality of (2) extends to our minimum imputed risk objective of (4). Below, we specify the δ θ and L(δ θ (x i ), y i ) we considered in our investigation.

Deterministic Decoding
A simple translation rule would define If this δ θ (x) is used together with a loss function L(δ θ (x i ), y i ) that is the negated BLEU score 9 , our minimum imputed risk objective of (4) is equivalent to MERT (Och, 2003) on the imputed training data. However, this would not yield a differentiable objective function. Infinitesimal changes to θ could result in discrete changes to the winning output string δ θ (x) in (6), and hence to the loss L(δ θ (x), y i ). Och (2003) developed a specialized line search to perform the optimization, which is not scalable when the number of model parameters θ is large.

Randomized Decoding
Instead of using the argmax of (6), we assume during training that δ θ (x) is itself random, i.e. the MT system randomly outputs a translation y with probability p θ (y | x). As a result, we will modify our objective function of (4) to take yet another expectation over the unknown y. Specifically, we will replace L(δ θ (x), y i ) in (4) with y p θ (y | x) L(y, y i ).
Now, the minimum imputed empirical risk objective of (4) becomes If the loss function L(y, y i ) is a negated BLEU, this is equivalent to performing minimum-risk training described by (Smith and Eisner, 2006;Li and Eisner, 2009) on the imputed data. 10 9 One can manipulate the loss function to support other methods that use deterministic decoding, such as Perceptron (Collins, 2002) and MIRA (Crammer et al., 2006). 10 Again, one may manipulate the loss function to support other probabilistic methods that use randomized decoding, such as CRFs (Lafferty et al., 2001). (8) is now differentiable, since each coefficient p θ (y | x) is a differentiable function of θ, and thus amenable to optimization by gradient-based methods; we use the L-BFGS algorithm (Liu et al., 1989) in our experiments. We perform experiments with the syntax-based MT system Joshua (Li et al., 2009a), which implements dynamic programming algorithms for second-order expectation semirings (Li and Eisner, 2009) to efficiently compute the gradients needed for optimizing (8).

Approximating
As mentioned at the end of Section 3.1, it is computationally infeasible to forward-translate each of the imputed reverse translations x ij . We propose four approximations that are computationally feasible. Each may be regarded as a different approximation of p φ (x | y i ) in equations (4) or (8).
k-best. For each y i , add to the imputed training set only the k most probable translations {x i1 , . . . x ik } according to p φ (x | y i ). (These can be extracted from X i using standard algorithms (Huang and Chiang, 2005).) Rescale their probabilities to sum to 1.
Sampling. For each y i , add to the training set k independent samples {x i1 , . . . x ik } from the distribution p φ (x | y i ), each with weight 1/k. (These can be sampled from X i using standard algorithms (Johnson et al., 2007).) This method is known in the literature as multiple imputation (Rubin, 1987).
Lattice. 11 Under certain special cases it is be possible to compute the expected loss in (3) exactly via dynamic programming. Although X i does contain exponentially many translations, it may use a "packed" representation in which these translations share structure. This representation may furthermore enable sharing work in forward-translation, so as to efficiently translate the entire set X i and obtain a distribution over translations y. Finally, the expected loss under that distribution, as required by equation (3), may also be efficiently computable.
All this turns out to be possible if (a) the posterior distribution p φ (x | y i ) is represented by an un-ambiguous weighted finite-state automaton X i , (b) the forward translation system δ θ is structured in a certain way as a weighted synchronous context-free grammar, and (c) the loss function decomposes in a certain way. We omit the details of the construction as beyond the scope of this paper.
In our experimental setting described below, (b) is true (using Joshua), and (c) is true (since we use a loss function presented by Tromble et al. (2008) that is an approximation to BLEU and is decomposable). While (a) is not true in our setting because X i is a hypergraph (which is ambiguous), Li et al. (2009b) show how to approximate a hypergraph representation of p φ (x | y i ) by an unambiguous WFSA. One could then apply the construction to this WFSA 12 , obtaining an approximation to (3). Rule-level Composition. Intuitively, the reason why the structure-sharing in the hypergraph X i (generated by the reverse system) cannot be exploited during forward translating is that when the forward Hiero system translates a string x i ∈ X i , it must parse it into recursive phrases.
But the structure-sharing within the hypergraph of X i has already parsed x i into recursive phrases, in a way determined by the reverse Hiero system; each translation phrase (or rule) corresponding to a hyperedge. To exploit structure-sharing, we can use a forward translation system that decomposes according to that existing parse of x i . We can do that by considering only forward translations that respect the hypergraph structure of X i . The simplest way to do this is to require complete isomorphism of the SCFG trees used for the reverse and forward translations. In other words, this does round-trip imputation (i.e., from y to x, and then to y ) at the rule level. This is essentially the approach taken by .

The Log-Linear Model p θ
We have not yet specified the form of p θ . Following much work in MT, we begin with a linear model where f (x, y) is a feature vector indexed by k. Our deterministic test-time translation system δ θ simply 12 Note that the forward translation of a WFSA is tractable by using a lattice-based decoder such as that by Dyer et al. (2008). outputs the highest-scoring y for fixed x. At training time, our randomized decoder (Section 3.3.2) uses the Boltzmann distribution (here a log-linear model) p θ (y | x) = e γ·score(x,y) Z(x) = e γ·score(x,y) y e γ·score(x,y ) (10) The scaling factor γ controls the sharpness of the training-time distribution, i.e., the degree to which the randomized decoder favors the highest-scoring y. For large γ, our training objective approaches the imputed risk of the deterministic test-time system while remaining differentiable.
In a task like MT, in addition to the input x and output y, we often need to introduce a latent variable d to represent the hidden derivation that relates x to y. A derivation d represents a particular phrase segmentation in a phrase-based MT system (Koehn et al., 2003) and a derivation tree in a typical syntaxbased system (Galley et al., 2006;Chiang, 2007). We change our model to assign scores not to an (x, y) pair but to the detailed derivation d; in particular, now the function f that extracts a feature vector can look at all of d. We replace y by d in (9)-(10), and finally define p θ (y|x) by marginalizing out d, where D(x, y) represents the set of derivations that yield x and y.

Minimum Imputed Risk vs. EM
The notion of imputing missing data is familiar from other settings (Little and Rubin, 1987), particularly the expectation maximization (EM) algorithm, a widely used generative approach. So it is instructive to compare EM with minimum imputed risk. One can estimate θ by maximizing the loglikelihood of the data {(x i , y i ), i = 1, . . . , N } as If the x i 's are missing, EM tries to iteratively maximize the marginal probability: 924 The E-step of each iteration comprises computing x p θt (x | y i ) log p θ (x, y i ), the expected loglikelihood of the complete data, where p θt (x | y i ) is the conditional part of p θt (x, y i ) under the current iterate θ t , and the M-step comprises maximizing it: (14) Notice that if we replace p θt (x|y i ) with p φ (x | y i ) in the equation above, and admit negated loglikelihood as a loss function, then the EM update (14) becomes identical to (4). In other words, the minimum imputed risk approach of Section 3.1 differs from EM in (i) using an externally-provided and static p φ , instead of refining it at each iteration based on the current p θt , and (ii) using a specific loss function, namely negated log-likelihood.
So why not simply use the maximum-likelihood (EM) training procedure for MT? One reason is that it is not discriminative: the loss function (e.g. negated BLEU) is ignored during training.
A second reason is that training good joint models p θ (x, y) is computationally expensive. Contemporary MT makes heavy use of log-linear probability models, which allow the system designer to inject phrase tables, linguistic intuitions, or prior knowledge through a careful choice of features. Computing the objective function of (14) in closed form is difficult if p θ is an arbitrary log-linear model, because the joint probability p θ (x i , y i ) is then defined as a ratio whose denominator Z θ involves a sum over all possible sentence pairs (x, y) of any length.
By contrast, our discriminative framework will only require us to work with conditional models. While conditional probabilities such as p φ (x | y) and p θ (y | x) are also ratios, computing their denominators only requires us to sum over a packed forest of possible translations of a given y or x. 13 In summary, EM would impute missing data using p θ (x | y) and predict outputs using p θ (y | x), both being conditional forms of the same joint model p θ (x, y). Our minimum imputed risk training method is similar, but it instead uses a pair of

Experimental Results
We report results on Chinese-to-English translation tasks using Joshua (Li et al., 2009a), an open-source implementation of Hiero (Chiang, 2007).

IWSLT Task
We train both reverse and forward baseline systems. The translation models are built using the corpus for the IWSLT 2005 Chinese to English translation task (Eck and Hori, 2005), which comprises 40,000 pairs of transcribed utterances in the travel domain. We use a 5-gram language model with modified Kneser-Ney smoothing (Chen and Goodman, 1998), trained on the English (resp. Chinese) side of the bitext. We use a standard training pipeline and pruning settings recommended by (Chiang, 2007).

NIST Task
For the NIST task, the TM is trained on about 1M parallel sentence pairs (about 28M words in each language), which are sub-sampled from corpora distributed by LDC for the NIST MT evaluation using a sampling method implemented in Joshua. We also used a 5-gram language model, trained on a data set consisting of a 130M words in English Gigaword (LDC2007T07) and the bitext's English side.

Feature Functions
We use two classes of features f k for discriminative training of p θ as defined in (9).

Regular Hiero Features
We include ten features that are standard in Hiero (Chiang, 2007). In particular, these include one baseline language model feature, three baseline translation models, one word penalty feature, three features to count how many rules with an arity of zero/one/two are used in a derivation, and two features to count how many times the unary and binary glue rules in Hiero are used in a derivation.

Target-rule Bigram Features
In this paper, we do not attempt to discriminatively tune a separate parameter for each bilingual rule in the Hiero grammar. Instead, we train several hundred features that generalize across these rules.
For each bilingual rule, we extract bigram features over the target-side symbols (including nonterminals and terminals). For example, if a bilingual rule's target-side is "on the X 1 issue of X 2 " where X 1 and X 2 are non-terminals (with a position index), we extract the bigram features on the, the X, X issue, issue of, and of X. (Note that the position index of a non-terminal is ignored in the feature.) Moreover, for the terminal symbols, we will use their dominant POS tags (instead of the symbol itself). For example, the feature the X becomes DT X. We use 541 such bigram features for IWSLT task (and 1023 such features for NIST task) that fire frequently.

IWSLT Task
In addition to the 40,000 sentence pairs used to train the baseline generative models (which are used to compute the features f k ), we use three bilingual data sets listed in Table 1, also from IWSLT, for discriminative training: one to train the reverse model p φ (which uses only the 10 standard Hiero features as described in Section 5.2.1), 14 one to train the forward model δ θ (which uses both classes of features described in Section 5.2, i.e., 551 features in total), and one for test.
Note that the reverse model φ is always trained using the supervised data of Dev φ, while the forward model θ may be trained in a supervised or semisupervised manner, as we will show below.
In all three data sets, each Chinese sentence x i has 16 English reference translations, so each y i is actually a set of 16 translations. When we impute data from y i (in the semi-supervised scenario), we 14 Ideally, we should train φ to minimize the conditional cross-entropy (5) as suggested in section 3.2. In the present results, we trained φ discriminatively to minimize risk, purely for ease of implementation using well versed steps.

Data set Purpose
# of sentences Chinese English Dev φ training φ 503 503×16 Dev θ training θ 503 * 503×16 Eval θ testing 506 506×16 Table 1: IWSLT Data sets used for discriminative training/test. Dev φ is used for discriminatively training of the reverse model φ, Dev θ is for the forward model, and Eval θ is for testing. The star * for Dev θ emphasizes that some of its Chinese side will not be used in the training (see Table 2 for details).
actually impute 16 different values of x i , by using p φ to separately reverse translate each sentence in y i . This effectively adds 16 pairs of the form (x i , y i ) to the training set (see section 3.4), where each x i is a different input sentence (imputed) in each case, but y i is always the original set of 16 references.

NIST Task
For the NIST task, we use MT03 set (having 919 sentences) to tune the component parameters in both the forward and reverse baseline systems. Additionally, we use the English side of MT04 (having 1788 sentences) to perform semi-supervised tuning of the forward model. The test sets are MT05 and MT06 (having 1082 and 1099 sentences, respectively). In all the data sets, each source sentence has four reference translations.

Main Results
We compare two training scenarios: supervised and semi-supervised. The supervised system ("Sup") carries out discriminative training on a bilingual data set. The semi-supervised system ("+Unsup") additionally uses some monolingual English text for discriminative training (where we impute one Chinese translation per English sentence).
Tables 2 and 3 report the results for the two tasks under two training scenarios. Clearly, adding unsupervised data improves over the supervised case, by at least 1.3 BLEU points in IWSLT and 0.5 BLEU in NIST.

Results for Analysis Purposes
Below, we will present more results on the IWSLT data set to help us understand the behavior of the 926   Table 3: BLEU scores for semi-supervised training for NIST task. The "Sup" system is trained on MT03, while the "+Unsup" system is trained with additional 1788 English sentences from MT04. (Note that while MT04 has 1788×4 English sentences as it has four sets of references, we only use one such set, for computational efficiency of discriminative training.) A star * indicates a result that is signicantly better than the "Sup" baseline (paired permutation test, p < 0.05). methods proposed in this paper.

Imputation with Different Reverse Models
A critical component of our unsupervised method is the reverse translation model p φ (x | y). We wonder how the performance of our unsupervised method changes when the quality of the reverse system varies. To study this question, we used two different reverse translation systems, one with a language model trained on the Chinese side of the bitext ("WLM"), and the other one without using such a Chinese LM ("NLM").  Table 4: BLEU scores for unsupervised training with/without using a language model in the reverse system. A data size of 101 means that we use only the English sentences from a subset of Dev θ containing 101 Chinese sentences and 101×16 English translations; for each English sentence we impute the 1-best Chinese translation. "WLM" means a Chinese language model is used in the reverse system, while "NLM" means no Chinese language model is used. In addition to reporting the BLEU score on Eval θ, we also report "Imputed-CN BLEU", the BLEU score of the imputed Chinese sentences against their corresponding Chinese reference sentences.
BLEU point in the forward translations. Still, even with the worse imputation (in the case of "NLM"), our forward translations improve as we add more monolingual data.

Imputation with Different k-best Sizes
In all the experiments so far, we used the reverse translation system to impute only a single Chinese translation for each English monolingual sentence. This is the 1-best approximation of section 3.4. Table 5 shows (in the fully unsupervised case) that the performance does not change much as k increases. 16 This may be because that the 5-best sentences are likely to be quite similar to one another (May and Knight, 2006). Imputing a longer k-best list, a sample, or a lattice for x i (see section 3.4) might achieve more diversity in the training inputs, which might make the system more robust.  Table 5: BLEU scores for unsupervised training with different k-best sizes. We use 101×16 monolingual English sentences, and for each English sentence we impute the k-best Chinese translations using the reverse system. the expected loss of the forward translations with respect to the observed output) is minimized. This matches the intuition that the probabilistic "roundtrip" translation from the target-language sentence to the source-language and back should have low expected loss. We applied our method to two Chinese to English machine translation tasks (i.e. IWSLT and NIST). We showed that augmenting supervised data with unsupervised data improved performance over the supervised case (for both tasks).
Our discriminative model used only a small amount of training data and relatively few features. In future work, we plan to test our method in settings where there are large amounts of monolingual training data (enabling many discriminative features). Also, our experiments here were performed on a language pair (i.e., Chinese to English) that has quite rich bilingual resources in the domain of the test data. In future work, we plan to consider lowresource test domains and language pairs like Urdu-English, where bilingual data for novel domains is sparse.