Learning a Phrase-based Translation Model from Monolingual Data with Application to Domain Adaptation

Currently, almost all of the statistical machine translation (SMT) models are trained with the parallel corpora in some specific domains. However, when it comes to a language pair or a different domain without any bilingual resources, the traditional SMT loses its power. Recently, some research works study the unsupervised SMT for inducing a simple word-based translation model from the monolingual corpora. It successfully bypasses the constraint of bitext for SMT and obtains a relatively promising result. In this paper, we take a step forward and propose a simple but effective method to induce a phrase-based model from the monolingual corpora given an automatically-induced translation lexicon or a manually-edited translation dictionary. We apply our method for the domain adaptation task and the extensive experiments show that our proposed method can substantially improve the translation quality.

Introduction
During the last decade, statistical machine translation has made great progress. Novel translation models, such as phrase-based models (Koehn et a., 2007), hierarchical phrase-based models (Chiang, 2007) and linguistically syntax-based models (Liu et a., 2006;Huang et al., 2006;Galley, 2006;Zhang et al, 2008;Chiang, 2010;Zhai et al., , 2012 have been proposed and achieved higher and higher translation performance. However, all of these state-of-the-art translation models rely on the parallel corpora to induce translation rules and estimate the corresponding parameters.
It is unfortunate that the parallel corpora are very expensive to collect and are usually not available for resource-poor languages and for many specific domains even in a resource-rich language pair.
Recently, more and more researchers concentrated on taking full advantage of the monolingual corpora in both source and target languages, and proposed methods for bilingual lexicon induction from non-parallel data (Rapp, 1995(Rapp, , 1999Koehn and Knight, 2002;Haghighi et al., 2008;Daumé III and Jagarlamudi, 2011) and proposed unsupervised statistical machine translation (bilingual lexicon is a byproduct) with only monolingual corpora (Ravi and Knight, 2011;Nuhn et al., 2012;Dou and Knight, 2012).
In the bilingual lexicon induction (Koehn and Knight, 2002;Haghighi et al., 2008;Daumé III and Jagarlamudi, 2011), with the help of the orthographic and context features, researchers adopted an unsupervised method, such as canonical correlation analysis (CCA) model, to automatically induce the word translation pairs between two languages from non-parallel data only requiring that the monolingual data in each language are from a fairly comparable domain.
The unsupervised statistical machine translation method (Ravi and Knight, 2011;Nuhn et al., 2012;Dou and Knight, 2012) viewed the translation task as a decipherment problem and designed a generative model with the objective function to maximize the likelihood of the source language monolingual data. To tackle the large-scale vocabulary, they mainly considered the word-based model (e.g. IBM Model 3) and applied the Bayesian method with Gibbs sampling or slice sampling. Finally, they used the learned translation model directly to translate unseen data (Ravi and Knight, 2011;Nuhn et al., 2012) or incorporated the learned bilingual lexicon as a new in-domain translation resource into the phrase-based model which is trained with out-of-domain data to improve the domain adaptation performance in machine translation (Dou and Knight, 2012).
In this paper, we focus on exploring this direction and propose a simple but effective meth-od to induce the phrase-level translation rules from monolingual data. The main idea of our method is to divide the phrase-level translation rule induction into two steps: bilingual lexicon induction and phrase pair induction.
Since many researchers have studied the bi-lingual lexicon induction, in this paper, we mainly concentrate ourselves on phrase pair in-duction given a probabilistic bilingual lexicon and two in-domain large monolingual data (source and target language). In addition, we will further introduce how to refine the induced phrase pairs and estimate the parameters of the induced phrase pairs, such as four standard translation features and phrase reordering feature used in the conventional phrase-based models (Koehn et al., 2007). The induced phrase-based model will be used to help domain adaptation for machine translation.
In the rest of this paper, we first explain with examples to show what new translation knowledge can be learned with our proposed phrase pair induction method (Section 2), and then we introduce the approach for probabilistic bilingual lexicon acquisition in Section 3. In Sec-tion 4 and 5, we respectively present our method for phrase pair induction and introduce an ap-proach for phrase pair refinement and parameter estimation. Section 6 will show the detailed ex-periments for the task of domain adaptation. We will introduce some related work in Section 7 and conclude this paper in Section 8.

What Can We Learn with Phrase Pair Induction?
Readers may doubt that if phrase pair induction is performed only using bilingual lexicon and monolingual data, what new translation knowledge can be learned?
The bilingual lexicon can only express the translation equivalence between source-and tar-get-side word pair and has little ability to deal with word reordering and idiom translation. In contrast, phrase pair induction can make up for this deficiency to some extent. Furthermore, our method is able to learn some unknown word translations.
From the induced phrase pairs with our meth-od, we have conducted a deep analysis and find that we can learn three kinds of new translation knowledge: 1) word reordering in a phrase pair; 2) idioms; and 3) unknown word translations. Table 1 gives examples for each of the three kinds. For the first example, the source and tar-get phrase are extracted respectively from mono-lingual data, each word in the source phrase has a translation in the target phrase, but the word order is different. The word order encoded in a phrase pair is difficult to learn in a word-based SMT. In the second example, the italic source word corresponds to two target words (in italic), and the phrase pair is an idiom which cannot be learned from word-based SMT. In the third ex-ample, as we learn from the source and target monolingual text that the words around the italic ones are translations with each other, thus we cannot only extract a new phrase pair but also learn a translation pair of unknown words in italic.

Probabilistic Bilingual Lexicon Acquisition
In order to induce the phrase pairs from the in-domain monolingual data for domain adaptation, the probabilistic bilingual lexicon is essential.
In this paper, we acquire the probabilistic bi-lingual lexicon from two approaches: 1) build a bilingual lexicon from large-scale out-of-domain parallel data; 2) adopt a manually collected in-domain lexicon. This paper uses Chinese-to-English translation as a case study and electronic data is the in-domain data we focus on.
In Chinese-to-English translation, there are lots of parallel data on News. Here, we utilize about 2.08 million sentence pairs 1 in News domain to learn a probabilistic bilingual lexicon. Basically, we can use GIZA++ (Och, 2003) to get the probabilistic lexicon. However, the problem is that each source-side word associates too many possible translations which contain much noise. For instance, in the lexicon obtained with GIZA++, each source-side word has about 13 translations on average. The noise of the lexicon can influence the accuracy of the induced phrase pairs to a large extent. To learn a lexicon with a high precision, we follow Munteanu and Marcu (2006) to apply Log-Likelihood-Ratios (Dunning, 1993;Melamed, 2000;Moore, 2004aMoore, , 2004b to estimate how strong the association is between a source-side word and its aligned target-side word. We employ the same algorithm used in (Munteanu and Marcu, 2006) which first use the GI-ZA++ (with grow-diag-final-and heuristic) to obtain the word alignment between source and target words, and then calculate the association strength between the aligned words. After using the log-likelihood-ratios algorithm 2 , we obtain a probabilistic bilingual lexicon with bidirectional translation probabilities from the out-of-domain data. In the final lexicon, the number of average translations is only 5. We call this lexicon LLRlex.
In the electronic domain, we manually collected a lexicon which contains about 140k entries. It should be noted that there is no translation probability in this lexicon. In order to assign probabilities to each entry, we apply the Corpus Translation Probability which used in (Wu et al., 2008): given an in-domain source language monolingual data, we translate this data with the phrase-based model trained on the out-of-domain News data, the in-domain lexicon and the indomain target language monolingual data (for language model estimation). With the source language data and its translation, we estimate the bidirectional translation probabilities for each entry in the original lexicon. For the entries whose translation probabilities are not estimated, we just assign a uniform probability. That is if a source word has n translations, then the translation probability of target word given the source word is 1/n. We call this lexicon Domain-lex.
We combine LLR-lex and Domain-lex to obtain the final probabilistic bilingual lexicon for phrase pair induction.

Phrase Pair Induction Method
Given a probabilistic bilingual lexicon and two monolingual data, we present a simple but effective method for phrase pair induction in this section.

A Naï ve Method
We first introduce a relatively naï ve way to extract phrase pairs from the given resources. For a source phrase (word sequence), we can reorder the words in the phrase (permutation) first, and then obtain the target phrases with the bilingual lexicon (translation), and finally check if the target phrase is in the target monolingual data. The algorithm is given in Figure 1. Figure 1 shows that the naï ve algorithm is very easy to implement. However, the time complexity is too high.
It is very time consuming. One may design smarter algorithms. For example, one can collect distinct n-grams from source and target monolingual data. Then, for a source-side phrase with length L, one can find the best translation candidate using the probabilistic bilingual lexicon from the target-side phrases with the same length L. The biggest disadvantage of these algorithms is that they can only induce phrase pair (with the Input: Probabilistic bilingual lexicon V (each source word s maps a translation set V[s]) Source language monolingual data S={s n } n=1...N Target language monolingual data T={t m } m=1...M Output: Phrase pairs P 1: For each distinct source-side phrase j i s in S:

2:
If each j ki ss  in V: 3: same length) encoding word reordering, but cannot learn phrase pairs in different length. Furthermore, they cannot learn idioms and unknown word translations from monolingual data. Obviously, these kind of approaches is not optimal.

Phrase Pair Induction with Inverted Index
In order to make the phrase pair induction both effective and efficient, we propose a method using inverted index data structure which is usually a central component of a typical search engine.
The inverted index is employed to represent the target language monolingual data. For a target language word, the inverted index not only records the sentence position in monolingual data, but also records the word position in a sentence. Some examples are shown in Table 2. By doing this, we do not need to iterate all the permutations of source language phrase j i s to explore possible phrase pairs encoding word reordering. Furthermore, it is possible to learn idiom translation and unknown word translations. We will elaborate how to induce phrase pairs with the help of inverted index.  Table 2: Some examples of inverted index for target language words, (2,5) means that "communication" occurs at the 5 th word of the 2 nd sentence in the target monolingual data.
The new algorithm for phrase pair induction is presented in Figure 2. Line 1 iterates all the distinct phrases in the source-side monolingual data. It can be implemented by collecting all the distinct n-grams in which n is the phrase length we are interested in (3 to 7 in this paper). For each distinct source-side phrase, Line 2-5 efficiently collects all the positions in the target monolingual data for the translations of each word in the source phrase. Line 6 sorts the positions so that we can easily find the position sequence belonging to a same sentence. Line 8-9 discards all the position sub-sequences that lack translations for more than one source-side words. That is to say we allow at most one unknown word in an induced phrase pair in order to make the induction more accurate. Line 10 and Line 12 is the core of this algorithm. We first define a constraint before detailing the algorithm. Constraint: we require that there exists at most one phrase in a target sentence that is the translation of the source-side phrase.
According to our analysis, it is not often to find that two phrases (length larger than 2) in a same sentence have the same meaning. Even if it happens, it is reasonable to keep the one with the highest probability. Given a position sequence belonging to a same sentence, Line 10 smoothes the probability of the single word gap according to the probabilities of the around words. Single word gap means that this word is not aligned but its left and right words are aligned with the words of the source-side phrase. Suppose the target sub-sequence is and ir t  is the only word that is not aligned with source-side words. We smooth the probability   | ir p t null  as follows: The above formula means that if the left or the right side only has one word, then the smoothed probability is one half of the minimum of the probabilities of the two neighbors, otherwise the smoothed probability is the average of the probabilities of the two neighbors. This smoothing strategy encourages that if more words around the un-aligned word are translations of the source-side phrase, then the gap word is more likely to belong to the translations of the sourceside phrase. After probability smoothing of the single gap word, we are ready to extract the candidate translation of the source-side phrase. Similar with Line 9 in Figure 2, we further filter the target continuous phrase if more than one word in source-side phrase has no translation in this target phrase. After that, we just choose the continuous target phrase with the largest probability if two or more continuous target phrases exist in the same target sentence. The probability of a target-side phrase given the source-side phrase is computed similar to that of (Koehn et al., 2003) except that we impose length normalization: where the alignment a is produced using probabilistic bilingual lexicon. If a target word in t is a gap word, we suppose there is a word alignment between the target gap word and the source-side null.
Similarly, we can compute the probability of source-side phrase given the target-side phrase . Line 13 in Figure 2 collects the induced phrase pairs. However, we find in the target monolingual data (1 million sentences) that each distinct word happens 110 times on average. Then, for a sources-side phrase with 7 words, the average length of positionArray will be 3850, since each source word has averagely 5 target translations (mentioned in Section 3). Therefore, the algorithm is relatively efficient in the average case.

Phrase Pair Refinement
Some of the phrase pairs induced in Section 4 may contain noise. According to our analysis, we find that the biggest problem is that in the target-side of the phrase pair, there are two or more identical words aligned to the same source-side word. For example, we extract a phrase pair as follows: 的 商业 信息 of business information of In the above phrase pair, there are two words "of" in the target side and the first one is redundant. The phrase pair induction algorithm presented in Section 4 cannot deal with this situation. In this section, we propose a simple approach to handle this problem. For each entry in LLR-lex, such as (的, of), we can learn two kinds of information from the out-of-domain wordaligned sentence pairs: one is whether the target translation is before or after the translation of the preceding source-side word (Order); the other is whether the target translation is adjacent with the translation of the preceding source-side word (Adjacency). If the source-side word is the beginning of the phrase, we calculate the corresponding information with the succeeding word instead of the preceding word. For the entries in Domain-lex, we constrain that the target translation should be adjacent with the translations of its source-side neighbors and translation order is the same with the source-side words.
With the Order and Adjacency information, we first check the order information, and then check the adjacency information if the duplicates cannot be handled using order information. For example, since (的, of) is an entry in LLRlex and we have learned that "of" is much more likely to be behind the translation of the succeeding word. Thus, the first word "of" can be discarded. This refinement can be applied before finding the phrase pair with maximum probability (Line 12 in Figure 2) so that the duplicate words do not affect the calculation of translation probability of phrase pair.

Translation Probability Estimation
It is well known that in the phrase-based SMT there are four translation probabilities and the reordering probability for each phrase pair.
The translation probabilities in the traditional phrase-based SMT include bidirectional phrase translation probabilities and bidirectional lexical weights. For the lexical weights, we can use the   |, lex p s t a and   |, lex p t s a computed in the above section without length normalization. However, for the phrase-level probability, we cannot use maximum likelihood estimation since the phrase pairs are not extracted from parallel sentences.
In this paper, we borrow and extend the idea of (Klementiev et al., 2012)  occurs in the context of s and in the entire source language monolingual data, and max n is the maximum number of occurrence of any source-side word in the source language monolingual data. The k-th element of t's vector can be computed with the same method. We finally normalize these vectors with L 2 -norm.
With the s's and t's contextual vector representations, we calculate two similarities: 1) project s's vector into target side t with the lexical mapping p(t|s), and then get the similarity by computing the cosine of two angles between t's and t 's vectors; 2) project t's vector into source side s with the lexical mapping p(s|t), and then obtain the similarity between s's and s 's vectors. These two contextual similarities will serve as two phrase-level translation probabilities.

Reordering Probability Estimation
For the reordering probabilities of newly induced phrase pairs, we can also follow Klementiev et al. (2012) to estimate these probabilities using source and target monolingual data. The method is to calculate six probabilities for monotone, swap or discontinuous cases. For the phrase pair (的 商业 信息, business information of), we find a source sentence containing 的 商业 信息, and find a target sentence containing business information of. If there is another phrase pair   , st , t exactly follows business information of and s occurs in the same source sentence with 的 商业 信息, then we compare the position relationship between s and 的 商业 信息. We increment the swap count if s is just before 的 商业 信息. After counting, we finally use maximum likelihood estimation method to compute the reordering probabilities.

Related Work
As far as we know, few researchers study phrase pair induction from only monolingual data.
There are three research works that are most related with ours. One is using an in-domain probabilistic bilingual lexicon to extract subsentential parallel fragments from comparable corpora (Munteanu and Marcu, 2006;Quirk et al., 2007;Cettolo et al., 2010). Munteanu and Marcu (2006) first extract the candidate parallel sentences from the comparable corpora and further extract the accurate sub-sentential bilingual fragments from the candidate parallel sentences using the in-domain probabilistic bilingual lexicon. Compared with their work, our focus is to induce phrase pairs directly from monolingual data rather than comparable data. Thus, finding the candidate parallel sentences is not possible in our situation.
Another is to make full use of monolingual data with transductive learning (Ueffing et al., 2007;Schwenk, 2008;Wu et al., 2008;Bertoldi and Federico, 2009). For the target-side monolingual data, they just use it to train language model, and for the source-side monolingual data, they employ a baseline (word-based SMT or phrasebased SMT trained with small-scale bitext) to first translate the source sentences, combining the source sentence and its target translation as a bilingual sentence pair, and then train a new phrase-base SMT with these pseudo sentence pairs. This method cannot learn idiom translations and unknown word translations.
The third is to estimate the translation parameters and reordering parameters using monolingual data given the phrase pairs (Klementiev et al., 2012). Their work supposes the phrase pairs are already given and then corresponding parameters can be learned with monolingual data. Different from their work, we concentrate ourselves on inducing phrase pairs from monolingual data and then borrow some ideas from theirs for parameter estimation. Furthermore, we extend their contextual similarity between source and target phrases to both directions.

Experimental Setup
Our purpose is to induce phrase pairs to improve translation quality for domain adaptation. We have introduced the out-of-domain data and the electronic in-domain lexicon in Section 3. Here we introduce other information about the in-domain data. Besides the in-domain lexicon, we have collected respectively 1 million monolingual sentences in electronic area from the web. They are neither parallel nor comparable because we cannot even extract a small number of parallel sentence pairs from this monolingual data using the method of (Munteanu and . We further employ experts to translate 2000 Chinese electronic sentences into English. The first half is used as the tuning set (elec1000tune) and the second half is employed as the testing set (elec1000-test).
We construct two kinds of phrase-based models using Moses (Koehn et al., 2007): one uses out-of-domain data and the other uses in-domain data. For the out-of-domain data, we build the phrase table and reordering table using the 2.08 million Chinese-to-English sentence pairs, and we use the SRILM toolkit (Stolcke, 2002) to train the 5-gram English language model with the target part of the parallel sentences and the Xinhua portion of the English Gigaword. For the in-domain electronic data, we first consider the lexicon as a phrase table in which we assign a constant 1.0 for each of the four probabilities, and then we combine this initial phrase table and the induced phrase pairs to form the new phrase table. The in-domain reordering table is created for the induced phrase pairs. An in-domain 5gram English language model is trained with the target 1 million monolingual data.
We use BLEU (Papineni et al., 2002) score with shortest length penalty as the evaluation metric and apply the pairwise re-sampling approach (Koehn, 2004) to perform the significance test.

Experimental Results
In this section, we first conduct experiments to figure out how the translation performance degrades when the domain changes. To better illustrate the comparison, we first use News data to evaluate the NIST evaluation tests and then use the same News data to evaluate the electronic test sets. For the NIST evaluation, we employ Chinese-to-English NIST MT03 as the tuning set and NIST MT05 as the test set. Table 3 gives the results. It is obvious that, it is relatively high when using the News training data to evaluate the same News test set. However, when the test domain is changed, the translation performance decreases to a large extent.
Given the in-domain bilingual lexicon and two monolingual data, previous works also proposed some good methods to explore the potential of the given data to improve the translation quality. Here, we implement their approaches and use them as our strong baseline. Wu et al. (2008) regards the in-domain lexicon with corpus translation probability as another phrase table and further use the in-domain language model besides the out-of-domain language model. Table 4 gives the results. We can see from the table that the domain lexicon is much helpful and significantly outperforms the baseline with more than 4.0 BLEU points. When it is enhanced with the in-domain language model, it can further improve the translation performance by more than 2.5 BLEU points. This method has made good use of in-domain lexicon and the target-side indomain monolingual data, but it does not take full advantage of the in-domain source-side monolingual data.
In order to use source-side monolingual data, Ueffing et al. (2007), Schwenk (2008), Wu et al. (2008) and Bertoldi and Federico (2009) employed the transductive learning to first translate the source-side monolingual data using the best configuration (baseline+in-domain lexicon+indomain language model) and obtain 1-best translation for each source-side sentence. With the source-side sentences and their translations, the new phrase table and reordering table are built. Then, these resources are added into the best configuration. The experimental results are presented in the last low of Table 4. From the results, we see that transductive learning can further improve the translation performance significantly by 0.6 BLEU points.
In tranductive learning, in-domain lexicon and both-side monolingual data have been explored. However, this method does not take full advantage of both-side monolingual data because it uses source and target monolingual data individually. In our method, we explore fully the source and target monolingual data to induce translation equivalence on the phrase level. In order to make the phrase pair induction more efficient, we first sort all the sentences in the both-side monolingual data according to the word hit rate in the bilingual lexicon. Then, we conduct six sets of experiments respectively on the first 100k, 200k, 300k, 500k and whole 1m sentences. All the experiments are run based on the configuration with BLEU 13.41 in Table 4, and we call this configuration BestConfig. Note that the unknown words are only allowed if the source-side of a phrase pair has more than 3 words. Table 5 shows the results.    We can see from the table that our method obtains the best translation performance. When using the first 100k sentences for phrase pair induction, it obtains a significant improvement over the BestConfig by 0.65 BLEU points and can outperform the transductive learning method. When we use more monolingual data, the performance becomes even better. The method of phrase pair induction using 300k sentences performs quite well. It outperforms the BestConfig significantly with an improvement of 1.42 BLEU points and it also performs much better than transductive learning method with gains of 0.82 BLEU points. With the monolingual data larger and larger, the gains become smaller and smaller because the word hit rate gets lower and lower. These experimental results empirically show the effectiveness of our proposed phrase pair induction method. A question remains that how many new phrase pairs are induced with different size of monolingual data. Here, we give respectively the statistics before and after filtering with the 1000 test sentences. Table 6 shows the statistics. We can see from the table that lots of new phrase pairs can be induced since the source and target monolingual data is in the same domain. However, since the source and target monolingual data is far from parallel, most of the phrase pairs are not long. For example, in the 108,948 distinct phrase pairs, we find that the phrase pair distribution according to source-side length is (3:50.6%, 4:35.6%, 5:3.3%, 6:9.8%, 7:0.7%). It is easy to see that the phrase pairs whose source-side length longer than 4 account for only a very small part.

Conclusion and Future Work
This paper proposes a simple but effective method to induce phrase pairs from monolingual data. Given the probabilistic bilingual lexicon and both-side monolingual data in the same domain, the method employs inverted index structure to represent the target-side monolingual data, and induce the translations for each distinct sourceside phrase with the help of the bilingual lexicon. We further propose an approach to refine the result phrase pairs to make them more accurate. We also introduce how to estimate the translation and reordering parameters for the induced phrase pairs with monolingual data. Extensive experiments on domain adaptation have shown that our method can significantly outperform previous methods which also focus on exploring the indomain lexicon and monolingual data.
However, through the analysis we find that our induced phrase pairs still contain some noise, such as the words in source-and target-side of the phrase pair are all aligned but the target-side phrase expresses the different meaning. Furthermore, our proposed method cannot learn expressions which are not lexical translations but are semantic ones. In the future, we will study further on these phenomena and propose new methods to handle these problems.


extfieldsofstudy:
Mathematics
Computer Science

text:
Parallel inversion of huge covariance matrices

An extremely common bottleneck encountered in statistical learning algorithms is inversion of huge covariance matrices, examples being in evaluating Gaussian likelihoods for a large number of data points. We propose general parallel algorithms for inverting positive definite matrices, which are nearly rank deficient. Such matrix inversions are needed in Gaussian process computations, among other settings, and remain a bottleneck even with the increasing literature on low rank approximations. We propose a general class of algorithms for parallelizing computations to dramatically speed up computation time by orders of magnitude exploiting multicore architectures. We implement our algorithm on a cloud computing platform, providing pseudo and actual code. The algorithm can be easily implemented on any multicore parallel computing resource. Some illustrations are provided to give a flavor for the gains and what becomes possible in freeing up this bottleneck.

Introduction
We consider a symmetric positive definite matrix K of order n where n is very large; on the order of 10,000s to millions or more. Our interest is in evaluation of K −1 , which cannot be computed sufficiently quickly using current algorithms. Even if we could compute the inverse, substantial numeric instability and inaccuracies would be expected due to propagation of errors arising from finite machine precision arithmetic (Trefethen and Bau III, 1997).
Typically in statistical applications, one needs to evaluate expressions such as K −1 A, where A is some matrix of appropriate order. Instead of directly computing K −1 , one popular approach is to consider the QR decomposition for K = QR and evaluate the expression, R −1 Q T A (Press et al., 2007). The matrix R is lower triangular and therefore R −1 can be evaluated by backward substitution, which requires O(n 2 ) flops instead of the O(n 3 ) flops for usual inversion. QR decomposition is known to be relatively more stable and efficient than other standard algorithms (Cox and Higham, 1997), a close competitor being the Cholesky decomposition.
The problem is that QR decomposition for K requires O(n 3 ) computations, which is of the same order as that of inversion and therefore prohibitively expensive for large n. Even with QR decomposition being relatively stable, for very large n, finite precision numerical errors are large. This is accentuated when the matrix K has a small spectrum, in the sense of fast decaying eigenvalues, as is typically obtained for covariance matrices obtained from smooth covariance kernels and in large least square problems among a host of other areas. However, this small spectrum, while being the bane of full rank algorithms, makes low-rank algorithms highly accurate.
We propose a new class of algorithms combining ideas from recent developments in linear algebra to find approximate low-rank decompositions. These low-rank decompositions provide several orders of magnitude improvement in speed while also improving accuracy. We propose a general blocking method to enable implementation in parallel on distributed computing architectures. We also consider accuracy of these approximations and provide bounds which are obtained with high probability. Our approach amalgamates ideas from three recent but apparently unrelated developments in numerical linear algebra and machine learning, which we briefly outline below.
There has been an increasing literature on approximating a matrix B of order n × n by BΩ where Ω is a matrix with random entries of order n × d with d << n (Halko et al., 2011). Originally motivated by Johnson Lindenstrauss transforms (Johnson and Lindenstrauss, 1984), these results have been used in a host of application areas, including compressive sensing (Donoho, 2006;Candès et al., 2006) and approximate Gaussian process regression (Banerjee et al., 2013). Typically Ω is populated by random entries, such as draws from Gaussian or Rademacher distributions. The product BΩ involves O(n 2 d) flops and can be itself expensive. Recent developments have shown that structured random matrices, such as random Hadamard transforms, may improve the efficiency significantly, bringing down the number of flops required to the order of n log d (Woolfe et al., 2008). In our case, we show that with certain classes of structured random matrices, Ω of order n × d, we have that KΩ has approximately the same column space as that of K (in a sense to be made precise later), and that the product KΩ may be computed efficiently. In addition, such structured random matrices spread the information of the matrix K, so that we obtain bounds on the minimum and maximum eigenvalues of KΩ, implying that approximate decompositions using KΩ are almost entirely devoid of usual inaccuracies arising from numerical conditioning.
Having efficiently obtained a tall skinny matrix KΩ from the square positive definite matrix K, we now consider recent developments which show that QR decompositions of such tall skinny matrices may be done efficiently in parallel (Constantine and Gleich, 2011;Agullo et al., 2010). The key idea is to partition the tall matrix into blocks and efficiently combine QR decomposition of each of these blocks. The main considerations are the choice of the number of blocks and column size d of the projection matrix Ω, which regulate computation time and accuracy. We provide some theoretical pointers while empirically demonstrating gains achieved on modest architectures. Approximate QR decompositions are useful in obtaining other approximate matrix decompositions, such as approximate spectral decompositions or matrix products. Several recent articles have focused on these issues and Halko et al. (2011) provides an excellent review for algorithms which switch between approximate decompositions. For illustration, we focus on the scenario when the principal interest is in evaluating a Gaussian type likelihood and fine tune the algorithms in this context. We finally consider performance versus other competitors.

Notational Preliminaries
We begin with some notations which we shall use through the rest of the article. In general we will represent matrices in the upper case, A, B etc and row or column vectors in the lower case, a, b etc. We shall use Frobenius and spectral matrix norms: A F = i,j a(i, j) 2 and A 2 = max{ Ax : x = 1} respectively, with the notation A ψ to imply the conclusion holds for both norms. We begin with a realvalued positive definite matrix K ∈ R n×n and k(i, j) will denote the (i, j)th element of K. A spectral decomposition of K = U DU T can be partitioned as, Since K is positive definite, the eigenvalues are positive and we assume without loss of generality that the diagonal matrix D of eigenvalues contains them in descending order of magnitude and the eigenvector matrix U is orthonormal. D m is therefore the matrix containing the m largest eigenvalues and U m the corresponding eigenvector matrix. By the Eckart Young theorem (Stewart, 1993), the best rank m approximation K m to K in terms of minimizing K − K m ψ is given by K m = U m D m U T m . We shall focus on the case where the matrix K has fast decaying eigenvalues. Some examples of problems where this happens is in large least squares optimization problems, Ax−b , with K = A T A, or with K being a covariance matrix generated for a Gaussian process at a dense set of locations with a smooth covariance function (Banerjee et al., 2013;Drineas et al., 2011). For some of our results, we will assume that eigenvalues of K decay at an exponential rate, by which we assume that there exists positive constants λ 1 , λ 2 such that d i,i ≤ λ 1 e λ2i ∀i. λ 1 , λ 2 shall be referred to as the proportionality and rate constants of exponential decay respectively.
An important consideration in this article will be the numerical stability of the algorithms and decompositions. One way in which the stability of a matrix decomposition can be measured is by how well the matrix in question is conditioned. The condition number of a positive semidefinite matrix A is given by c(A) = σ l /σ s , where σ l , σ s are the largest and smallest eigenvalues of the matrix respectively. For the positive definite matrix K it would be c(K) = d(1, 1)/d(n, n).

The Main Algorithm
Our fast approximate inversion algorithm uses three different ideas from numerical techniques and randomized linear algebra. The first is to approximate the spectrum of the large matrix K by post-multiplying it with a random matrix Ω, the second is to consider incorporating special structure in the random matrix Ω for faster evaluation of the product and lastly, blocking strategies, to enable computation of decomposition of the product KΩ on multicore architectures.

Column space approximation
As mentioned, we shall be considering positive definite matrices K with very fast decaying spectrums. In section §4, we provide theoretical justification to show that eigenvalues of positive definite matrices, produced as discrete realizations of a wide class of positive definite kernels, decay very fast. For a positive definite matrix K with fast decaying spectrums, it is reasonable to expect that if Ω is a n × r matrix with r n, populated with independent entries, then the product KΩ will capture most of the information in K, or approximate the column space of K. Such approximation schemes have been in focus in the rapidly expanding field of randomized linear algebra (refer to Halko et al. (2011) for a review). Banerjee et al. (2013) consider this method in the context of approximations for Gaussian processes for large data sets. In general, it is difficult to measure the amount of information captured in KΩ from K. One way to do this is to consider the generalized projection of K onto KΩ and then consider differences. Letting Q be an orthonormal basis for the column space of KΩ, we may consider the error, K − QQ T K ψ , where Q T Q = I, by virtue of orthonormality. In general, we may consider an approximation scheme with two objectives in mind.
1. Fix a target rank m and try to minimize the matrix norm error when using Ω having the target rank.
2. Fix a target matrix norm error and try to achieve that error with high probability with the smallest possible Ω.
We shall consider each of objectives in our projection approximation designs. The simplest possible choice for the random matrix Ω is a submatrix of columns from the n × n identity matrix (this submatrix is sometimes called a permutation submatrix, denoted by P ). This amounts to choosing r columns at random from the matrix K. Low rank approximations using a submatrix have been explored in a variety of contexts, including subset of regressors (Quinonero Candela and Rasmussen, 2005) and least squares (Drineas et al., 2011). With the advent of compressive sensing (Donoho, 2006;Candès et al., 2006), it has been shown that instead of restricting Ω to permutation submatrices, more general random matrices often have better performance (Drineas et al., 2011). Common choices of the random matrix Ω, which have been used in the compressive sensing and matrix approximation literature are: 1. elements of Ω are independent and identically distributed Gaussian random variables appropriately scaled, 2. elements of Ω are independent and identically distributed Rademacher variables, 3. elements of Ω are independently sampled from the uniform distribution on the unit sphere.
In general it has been shown that most choices, based on sampling from a centered distribution and then appropriately scaling, work in achieving accurate error bounds with high probability (Candès et al., 2006).

Approximations via structured random matrices
The biggest drawback of the approximation techniques in §3.4.1 is that the matrix product KΩ may itself be prohibitively expensive O(n 2 r) for very large K. Instead we may form the random matrix Ω such that the product KΩ can be evaluated faster, exploiting properties of the structured distribution.
Definition 1. We shall call a random matrix a structured random matrix if it is of the form Ω = cRT P , where • c is an appropriate scaling constant, such that the columns of Ω are orthonormal, • R is an n × n diagonal matrix whose diagonal elements are independent Rademacher entries, • T is an n × n appropriate orthogonal transform, facilitating the fast multiply, • P is an n × r permutation submatrix.
There are different choices of orthogonal transforms we may use, all of which facilitate the fast multiply, examples for T being the discrete Fourier matrix, discrete cosine matrix, discrete Hartley matrix or the Walsh-Hadamard matrix. The Walsh-Hadamard matrix in particular has been in focus in recent literature (Tropp, 2011) and tight bounds on the approximation accuracy have been obtained. The Walsh Hadamard matrix of order n × n is defined as, The Walsh-Hadamard and discrete Fourier transforms have some optimality with respect to preserving matrix coherence (Tropp, 2011;Boutsidis and Gittens, 2012), but a disadvantage is Walsh-Hadamard transforms exist only for positive integers which are powers of 2, while the discrete Fourier transforms may involve complex numbers. The alternatives, Hartley and discrete cosine transforms while being marginally suboptimal do not have these problems (Avron et al., 2010). In the following section, we formalize the theoretical setting and describe probabilistic error bounds for column space approximations using transforms with general orthonormal matrices (therefore valid for any of the aforementioned transforms).

Blocked decompositions for parallelization
Our focus has been on obtaining a low rank rectangular matrixK of order n × r, with r n, as an approximation to K, such that the column space ofK is very close to the column space of K. To measure the approximation error, we used the error in column space approximation via projection, measured as, ||K − QQ T K ψ , where Q is an orthogonal basis for the column space ofK. Such a Q is typically obtained from the QR decomposition ofK, where Q is an n × r matrix with orthonormal columns and R is an upper triangular matrix (lower triangular in some conventions). Obtaining a QR decomposition of an n × r matrix has computational cost O(nr 2 ). With the matrix multiplication for forming KΩ, the eventual computational cost is O(n 2 r) (for n r), which may be too expensive for very large n. To further reduce computational burden, we may consider parallelizing the computations. Matrix multiplication is trivially parallelizable, while most matrix decompositions are not (Choi et al., 1994). State-ofthe-art linear algebra algorithms, for example in ScaLAPACK routines (Blackford et al., 1997), QR, Cholesky or SV D decompositions, involve dense manipulations of the whole matrix, which are not parallelizable or cannot be blocked. In addition to computational cost, storing the full matrixK may be problematic in terms of memory requirements for very large n. However, in our case, we are interested in the QR decomposition of a tall matrix, number of rows being much greater than the number of columns, for which blocking strategies have been recently developed, exploiting modern parallel computing platforms like MapReduce. We elaborate on the blocking strategies next.
For notational clarity assume that n = 8r and consider a partition ofK, In the above each K i is of size 2r × r. Denote the QR factorization of the small matrix K i = Q i R i . Along the lines of the T SQR factorization presented in Constantine and Gleich (2011) we consider the following blocked factorization, Each of the small matrices Q i being orthogonal, the left factor of the decomposition in (1), Q b 1 , is formed by stacking the small matrices Q i as diagonal blocks, in an orthogonal matrix. It is also an orthogonal basis for the column space ofK and in many application areas it suffices to work with Q b 1 . In some other applications, we may require evaluation of the product Q b 1 x, where x is some vector of order r × 1, which can be evaluated in a similar blocked fashion, utilizing the small Q i 's. The matrix R b 1 , formed by stacking the small upper triangular matrices R i , is not upper triangular and hence the decomposition Q b 1 R b 1 is not a QR decomposition ofK. To achieve the QR decomposition ofK consider the following sequence of blocked and finally the QR decomposition of the 2r × r matrix R b 2 , where we can give the final QR decomposition ofK as, R b 3 by virtue of the construction is an upper triangular matrix, whileQ, a product of orthogonal matrices is itself orthogonal. A similar iteration can be carried out for any n, r by appropriately adjusting the block sizes. In some instances, if only (R b 3 ) −1 is required, the matrix product for formingQ need not be evaluated at all. On the other hand, if a product of the formQ T x is required for some some vector x, the entire product maybe evaluated via blocking. We give more details about such strategies in our examples.
There are other ways of constructing a sequence of small blocked QR decompositions leading to a QR decomposition of the tall matrixK. The scheme presented in (4) is the one allowing for maximum assimilation of computations at the same time. Depending on the architecture being used, this may not be the most optimal choice. Below we present another analogous blocked QR scheme, which allows for lesser number of computations to be carried on simultaneously, but possibly uses lesser communication between units of a multicore architecture. Assume, again for notational clarity that n = 8r. The matrixK is now partitioned into 8 blocks, where, Q 1 R 1 is the orthogonal matrix from the QR decomposition of K 1 , Q 2 is from the QR decomposition of R 1 K 2 , etc and I denotes identity matrix of appropriate size. The matrices R 1 , R 2 . . . are intentionally omitted from (5)  . This alternative blocking scheme involves minimum communication between units of the multicore architecture. The optimal schemes will usually be a mixture of the strategies (2) and (5). For more variants of blocking algorithms to get QR decompositions for tall matrices, we refer the readers to Constantine and Gleich (2011); Agullo et al. (2010). We now consider the theoretical gains possible from the blocked algorithm, while considering block size. Consider the ideal scenario, when there is no overhead or communication cost and each unit of the multicore architecture has equal speed. We consider the strategy corresponding to (2). Each of the small QR decompositions involve an 2r × r matrix, which has O(r 3 ) cost. At the first iteration, as in (2), using the same number of units as the number of blocks, say b, we can perform all the computations in O(r 3 ). In fact, it is easy to see that the computational order of O(r 3 ) holds true for each of the subsequent iterations. The total number of iterations needed, following the first strategy would be log 2 (b) + 1, which brings the total computational cost to O(log 2 (b)r 3 ) for the QR ofK. Forming the product KΩ and additional matrix multiplications (which will almost always be needed in applications after the QR decomposition), bring the cost to O(n 2 r) without parallelization, whereas with the blocking algorithm, provided we perform the matrix multiplications with the same blocked structure, we achieve the whole computation in O(log 2 (b)r 3 ). Following the first strategy and from the above discussion, in the ideal scenario, the number of blocks to be used is b = n 2r + 1. The theoretical possible maximum speed-up, assuming no communication cost, is where we measure the speed-up from Amdahl's equation (Mattson et al., 2005), as proportion of speed of the unparallelized computations to the speed of the parallelized computations. In this subsection we derive some motivating results justifying our assertion in §3.1 that positive definite matrices obtained as dense discrete realizations of positive definite kernels will have a very fast decaying spectrum, enabling us to obtain good approximations using much lower rank matrices. The ideas are not entirely new, we adapt abstract results from theory of integral equations and stochastic series expansions to our context. To prepare the background, we start with the well-known Mercer's theorem for positive definite functions.
Definition 2. Let D be a compact metric space and C a function, C : D × D → R + ∪ {0}. C is said to be positive definite if for all n, scalars c 1 , . . . , c n ∈ R and x 1 , . . . , x n ∈ D, we have i j c i c j C(x i , x j ) > 0. It follows trivially from the definition that if K is an n × n matrix with k(i, j) = c(x i , x j ) for x 1 , . . . , x n ∈ D, then K is a positive definite matrix.
With this definition, we give the following version of Mercer's theorem for the positive definite function C(·, ·) (Kühn, 1987), where D is a compact metric space, there exists a real valued scalar sequence, {λ i } ∈ l 1 , λ 1 ≥ λ 2 ≥ . . . ≥ 0 and an orthonormal sequence of functions, {φ i (·)} ∈ L 2 (D), φ i (·) : D → R, such that, ∀x 1 , x 2 ∈ D and this sequence converges in the mean square sense and uniformly in D.
We omit the the proof of this theorem and refer the readers to (Kühn, 1987).Mercer's theorem is a generalization of the spectral decomposition for matrices extended to functions. {λ i }, {φ(·)} are referred to as the eigenvalues and eigenfunctions respectively of C(·, ·). Alternatively, suppose D ⊂ R d and considering Lebesgue measure, we have the following integral equation for the eigenvalues and eigenfunctions, analogous to the definition of matrix eigenvalues and eigenfunctions, ∀x ∈ D, i ∈ N. The eigenvalues and eigenfunctions of C(·, ·) can be approximated using the eigenvalues and eigenvectors of the positive definite matrix K, obtained by evaluating C(·, ·) at any set of n locations x 1 , . . . x n ∈ D. To see this, consider a discrete approximation of the integral equation (6), using the points x 1 , . . . , x n , with equal weights in the quadrature rule, Substituting y = x 1 , . . . , x n in (7) we get a system of linear equations which is equivalent to the matrix eigenproblem, KU i = d(i, i)U i , i = 1, . . . , n, where U i is the i th row of the eigenvector matrix U . This approximation is related to the Galerkin method for integral equations and the Nystrom approximation (Baker and Baker, 1977;Delves and Mohamed, 1988). It also corresponds to finite truncation of the expansion from the Mercer's theorem. In general it can be shown that d(i, i)/n and √ n u(i, j) converge to λ i and φ i (x j ) respectively as n → ∞ (Baker and Baker, 1977).
The accuracy of the finite truncation of the Mercer expansion has been in recent focus and error bounds have been obtained depending on the degree of smoothness of the positive definite function (Todor, 2006;Schwab and Todor, 2006), in the context of stochastic uncertainty quantification. To borrow from these ideas and adapt the abstract results in our case, we begin with definitions quantifying the smoothness of the covariance functions.
Definition 3. Let C(·, ·) be a positive definite function on D × D, where D is a compact metric space. We call C(·, ·) piecewise Sobolev (p, q) smooth, for some p, q ∈ N, if there exists a finite disjoint partition {D j , j ∈ J} of D, with D ⊂ ∪ j∈J D j , such that for any pair (j 1 , j 2 ) ∈ J, C(·, ·) is Sobolev (p, q) on D j2 ×D j2 .
Definition 4. Let C(·, ·) be a positive definite function on D × D, where D is a compact metric space. We call C(·, ·) piecewise analytic smooth, for some p, q ∈ N, if there exists a finite disjoint partition {D j , j ∈ J} of D, with D ⊂ ∪ j∈J D j , such that for any pair (j 1 , j 2 ) ∈ J, C(·, ·) is analytic on D j2 × D j2 .
Commonly used covariance functions such as the squared exponential, C(x, y) = c 1 exp(−c 2 x − y 2 ) and the Matern function, depending on the value of the parameter ν, fall into one of the above categories. In fact the squared exponential function is Sobolev smooth on any compact domain and we can give stronger results for decay of its eigenvalues.
Using the couple of definitions categorizing smoothness, we give the following result, which is analogous to a functional version of the Eckart-Young theorem.
Lemma 2. Let D ⊂ R d . For any m ∈ N, and let F m be an m dimensional closed subspace of positive definite functions on D × D, where by m dimensional we mean that ∃ sequence of orthonormal function {ψ j (·)} m j=1 on D such that{ψ 2 j (·)} m j=1 spans F m , with positive coefficients by virtue of positive definiteness. Let C Fm (·, ·) be the projection of C(·, ·) onto F m , and infimum of the errors in approximating C(·, ·) by m dimensional functions ∈ F m , ie E m = inf Fm { C − C Fm 2 , C m ∈ F m }. Then, (i) If C(·, ·) is piecewise Sobolev (p, q) ∀(p, q) ∃ for each s > 1, a positive constant c s depending only on C(·, ·), s, such that E m ≤ c s j≥m j −s . Henceforth we call this bound E s m .
(ii) If C(·, ·) is piecewise analytic ∃ , positive constants c 1 , c 2 depending only on C(·, ·), such that E m ≤ c s j≥m c 1 exp c2m 1/d . Henceforth we call this bound E a m .
Proof. Let the Mercer expansion for C(x, y) be j≥1 λ j φ j (x)φ j (y). Define a random function, f (x) = j≥1 η j λ j φ j (x), where η j are independent standard normal random variables. Let C Fm (x, y) = m j=1 θ j ψ j (x)ψ j (y) corresponding to the basis {ψ 2 (·)} m j=1 for F m . Similar to f , corresponding to C Fm (·, ·), define a random function, f Fm (x) = j≥1 ζ j θ j ψ j (x), where ζ j 's are independent random normal variables. Defining f − f Fm = E(f − f Fm ) 2 , and applying theorem 2.7 in Schwab and Todor (2006), we get that, E m = j≥m λ j . Then the result follows by applications of Corollary 3.3 and Proposition 3.5 in Todor (2006).
The above lemma quantifies the finite truncation accuracy of the Mercer theorem for smooth kernels. The bounds obtained are optimal and in general cannot be improved. The spaces F m 's can be made more general to encompass all square integrable bivariate functions, but the proof becomes more involved in that case and we omit it for the sake of brevity. In case of the squared exponential kernel, which is smooth of all orders, the finite truncation is even sharper and we give the following corollary for it: Corollary 3. Let C(x, y) = θ 1 exp −θ2 x−y 2 , for x, y ∈ D, compact ⊂ R d and θ 1 , θ 2 > 0. Using all notations as in lemma 2, ∃ positive constant c θ1,θ2 such that E m ≤ c θ1,θ2 j≥m θ j 1/d 2 Γ(j 1/d /2) . Henceforth we shall call this bound E e r .
The proof is exactly similar to the proof of Lemma 2 and an application of Proposition 3.6 in Todor (2006).
We now describe some results relating to the matrix eigenvalues, using the strength of Lemma 2 above. Proof. We begin with the discrete approximate solution of the integral equation, using the Galerkin technique, described in equation (7). Note that this method coincides with the Raleigh-Ritz method with the identity matrix in Baker and Baker (1977), since C(·, ·) is symmetric positive definite. Applying then theorem 3.31 on page 322 of Baker and Baker (1977), we have d(m, m) ≤ nλ m , where λ m is the m th eigenvalue from the Mercer expansion of C(·, ·). The result then follows by a straight-forward application of Lemma 2 above and Corollary 3.3, Proposition 3.5 in Todor (2006). For the squared exponential kernel, it is straightforward application of Corollary 3 above and Proposition 3.6 in Todor (2006).
This theorem shows that for most covariance functions, the positive definite matrices that are generated by their finite realizations have eigenvalues that decay extremely fast, which was our assertion at the start of the section. In the simulated experiments we compute eigenvalues of some such covariance kernels and show that empirical results support our theory.

Approximation accuracy and condition numbers
We now present some results regarding the accuracy of our approximation algorithms. We shall be concerned with the column space approximation error, when we useK = Ω instead of K, measured as K − QQ T K ψ , where Ω is a structured random matrix as in definition 1 and Q is orthonormal basis for the columns ofK. Such a Q can be obtained from the QR decomposition ofK.
Let Ω be an n × r structured random matrix as formulated in definition 1 and Q be the left factor from the QR decomposition of KΩ. Choose r, k such that 4[ √ k + √ 8 ln kn] 2 ln K ≤ r ≤ n. With probability at least (1 − O(1/k)), the following hold true, (i) If C(·, ·) is piecewise Sobolev smooth ∀(p, q), then (a) K − QQ T K 2 ≤ (1 + 7n/r)c s r −s , where c s is a positive constant depending on s for any s > 1; (b) K − QQ T K F ≤ (1 + 7n/r)E s r . (ii) If C(·, ·) is piecewise analytic, then (a) K − QQ T K 2 ≤ n(1 + 7n/r)c 1 exp c2m 1/d , where c 1 , c 2 are positive constants ; (b) K − QQ T K F ≤ n(1 + 7n/r)E a r . (iii) In particular, for a squared exponential kernel, C(x, y) = θ 1 exp −θ2 x−y 2 , then (a) K − QQ T K 2 ≤ (1 + 7n/r)c θ1,θ2 θ m 1/d 2 Γ(m 1/d /2) , where c θ1,θ2 is a positive constant depending on θ 1 , θ 2 ; (b) K − QQ T K F ≤ (1 + 7n/r)E e r , where E s r , E a r , E e r are bounds as in Theorem 2 and Corollary 3. Proof. First note that, where ((K) TK ) + denotes the Moore-Penrose inverse of ((K) TK ). Let Q k R k denote the QR decomposition ofK, so that we have, Let K r be the best rank r approximation to K and let the spectral decomposition of K = U DU T . Then from the Eckart-Young theorem, K − K m 2 = d(r, r) and K − K m 2 = n j=r d(j, j). The result then follows by an application of our Theorem 4 and Theorem 11.2 in Halko et al. (2011).
In addition to providing accurate column space approximations, these orthogonal transforms spread the information of the matrix K and improve its conditioning, while preserving its geometry. By preserving the geometry, we mean preserving the norms of the eigenvectors -we explore more of this aspect empirically in the simulations. Any low rank approximation improves the conditioning, but it has been shown (Boutsidis and Gittens, 2012;Avron et al., 2010) that projections using the orthogonal transforms as above, improve conditioning substantially beyond what is achieved by just using any low rank approximation with high probability in special cases. Halko et al. (2011); Boutsidis and Gittens (2012) obtain bounds on the eigenvalues of AΩ where Ω is a structured random matrix and A is orthogonal. We are in interested in the stability of the numerical system,KR −1 k , as explained in the introduction of this article, where the QR decomposition ofK = Q k R k . We present the following result, without proof, trivially modifying results from Boutsidis and Gittens (2012); Avron et al. (2010), Theorem 6. Fix 0 < < (1/3) and choose 0 < δ < 1, r such that 6 2 [ (r) + 8 ln n/δ] 2 ln 2n/δ ≤ r ≤ n. Then with probability at least 1 − 2δ, we have the condition number of the linear system of interest, This improvement in conditioning number causes huge improvement in learning algorithms over other competing approaches as we demonstrate later.

Simulation Examples
We begin with empirical investigation into the eigenvalue decays of matrices realized from commonly used covariance kernels. We consider an equispaced grid of 100 points in [0, 1], to generate the 100×100 covariance matrices with 100 positive eigenvalues. We start with such a moderate size, 100, to illustrate that even for this small size, eigenvalues decay extremely fast, as was hypothesized and theoretically proven earlier.
This first one we use is the squared exponential kernel, The parameter θ 1 is a scaling parameter and does not change the eigendirections. In our present simulation, we set θ 1 = 1. The parameter θ 2 controls the smoothness of the covariance, as it decays with the distance x − y . We consider a range of θ 2 values, 0.05, 0.5, 1, 1.5, 2, 10, respectively, considering negligible decay with distance and very fast decay with distance. The plot of the eigenvalues is shown in figure 1. The figure shows similar rates of decay across the range of values of the smoothness parameter θ 2 . Depending on the smoothness parameter, the value of the largest eigenvalue changes by several orders, but in all cases, by the 10th largest eigenvalue, the sequence has approximately converged to 0, which indicates that for a low rank approximation, a choice of rank 10 would be sufficient for our purposes. This second one we use is the Matern covariance kernel (Williams and Rasmussen, 1996), where B ν is the modified Bessel function of the second kind. The parameter θ 1 is a scaling parameter and does not change the eigendirections, analogous to θ 1 of the squared exponential kernel. In our present simulation, we set θ 1 = 1, θ 2 = 1. The parameter ν controls the smoothness of the covariance, analogous to the parameter θ 2 of the squared exponential kernel. As a special case, as ν → ∞, we obtain the squared exponential kernel. In typical spatial applications, ν is considered to be within [0.5, 3], as the data are rarely informative about the smoothness levels beyond these. We consider a range of ν values, 0.5, 1, 1.5, 2, 2.5, 3, respectively, considering negligible decay with distance and very fast decay with distance. The plot of the eigenvalues is shown in figure 2. The figure shows similar rates of decay across the range of values of the smoothness parameter ν, exactly as exhibited with the squared exponential kernel. Once again, a low rank approximation of rank 10 would suffice in this case for approximating the 100 × 100 matrix.
Using the same set of simulations, we compute conditioning numbers, in turn for each of the covariance kernels. We compute the conditioning number of the the full 100 × 100 covariance matrix and its best rank m approximations for m = 5, 10, 15, 20, 50 respectively. Results for the squared exponential kernel are tabulated in table 1 and for the Matern covariance kernel in table 2. The full covariance matrix is extremely ill-conditioned in each case, the best rank m approximations improve the conditioning, as is to be expected. Even with the best rank m approximations, the conditioning numbers are still very large indicating very unstable linear systems. In the table we omit the exact digits for the very large numbers and just indicate their orders in terms of powers of 10. In case of the squared exponential kernel, condition numbers decrease from left to right and from top to bottom. This pattern is in general not true for the condition numbers of the Matern kernel, whose condition numbers are several orders smaller than that for the squared exponential but still to large for it to be stable linear systems. In general for the Matern kernel, as ν becomes larger, the kernel becomes smoother, the condition numbers are larger and the eigenvalues decay faster.
We next move to the simulations to consider the effect of blocking and gain in efficiency from parallelization. For large sample sizes, we will get very poor and time consuming inverse estimates, as demonstrated by the large conditioning numbers from the simulations of the previous section. To circumvent this, we start therefore with a known spectral decomposition and apply our approximations, pretending that the decomposition was unknown. We consider the known spectral decomposition as K = EDE T , where E is an orthonormalized matrix with randomly generated iid Gaussian entries. D is the diagonal matrix of eigenvalues in decreasing order or magnitude and for this simulation example we assume exponential decay with scale = 1 and rate 0.01. We generate this covariance matrix K for sample sizes n = 1000, 5000, 10000, 50000 respectively and apply our blocked approximate inversion algorithm. For the random matrix Ω used for the projection, we use three different choices, a matrix with scaled iid standard Gaussian entries, a structured random Hartley transform and a structured random discrete cosine transform. We have a total 64 cores at our disposal and to get a flavor of the gains possible by parallelization, we use 8, 16, 32, 64 cores in turn and measure the gain in efficiency as the ratio of time taken to the time for the algorithm without parallelization.
We report the gain in efficiency in figure 3. The figure reveals the efficiency gain by using the blocked algorithm. It shows that efficiency gain increases with number of cores in the larger sample sizes. Specifically when the sample size is 50, 000, it seems that more than 64 cores could have further increased the efficiency. The gains reported are obviously lower than the theoretical maximum possible gains, but are still substantial and have potential of further increase with larger number of cores in huge sample sizes.

Real data examples
We consider a real data examples with a very large data set, to demonstrate the large data handling power and efficiency gain via our approach. We consider subsets of the actual data available and variables such that the data fit into a Gaussian process framework.
The example we consider is a subset of the Sloan sky digital survey (Stoughton et al., 2007). In the survey, redshift and photometric measurements were made for a very large number of galaxies. The photometric measurements consist of 5 band measurements and are available for most of the galaxies, while the redshift measurements are not available for a large number of galaxies due to spectroscopy constraints. The main interest is therefore in predicting the red-shifts, given the other 5 measurements. For this example we consider a subset comprising of 180, 000 galaxies, whose photometric measurements and red-shift measurements are known. We hold out 30, 000 galaxies and pretend that redshifts are unknown for these galaxies, while using the remaining 150, 000 to train our Gaussian process model.
For the computations we use both a squared exponential kernel and a Matern kernel. The data are scaled and centered before calculations so that we are not concerned with inference regarding the θ 1 parameters for the covariance kernel. For the Matern kernel, we use a discrete uniform prior for the parameter ν based on 100 equispaced points on the grid [0.5, 3]. For the squared exponential kernel, we use a 100 trials, each having a random selection 1, 000 points from the training set to set a band for θ 2 , to estimate range of covariance of the data. It appears that it is sufficient to consider θ 2 in the range [0.05, 1]. We therefore use a grid on 100 uniformly spaced points in [0.05, 1] and analogous to the prior for ν, place a discrete uniform prior on θ 2 of the squared exponential kernel.
For data sets of this size it becomes extremely inefficient to use predictive processes or knot based approaches per se, without blocking or parallelization and knot selection is almost practically infeasible to attempt. We work with the modified predictive processes approach of Finley et al. (2009), without knot selection as competitors to our approach. The modified predictive processes is conceptually equivalent to the fully independent training conditional approach (Quinonero Candela and Rasmussen, 2005) or the Gaussian process approximation with pseudo inputs, with (Snelson and Ghahramani, 2006), with the same priors for the parameters as our approach, so that the only difference is in the way the covariance inversion is done.
For each of the approaches, including ours, one issue is how to select rank of the approximate projections. We use an approximate estimate of the covariance decay to come up with the low rank projection to be used. To calculate this we select 1, 000 data points spread out across the training set, in the following approximate manner, first select the pair of points which are most distant from each other, then select the point which is almost equidistant from the pair selected, and so on, till we have a well spread out selection. Here distance is measured as the Euclidean distance between the 5 variate photometric measurement vectors. Using these 1, 000 data points, we calculate the prior covariance matrix, using the Gaussian covariance kernel and the Matern kernel for each value of the smoothness parameter on its grid of values (Our simulation examples have revealed that the decay of eigenvalues depends on this smoothing parameter). We then estimate the rank such that Frobenius norm error in using the best low rank approximation from the Eckart Young theorem for the full covariance matrix would be no more than 0.001 (bearing in mind that this is an estimate of the best possible error, actual error in different projections will be more). In the computations for our approach, we use a different random projection to the targeted rank for each grid value of the smoothing parameter, while for the knot based approaches, a different random selection of knots is used.
For our approach, we use 3 different choices of the projection matrix Ω, a matrix with scaled iid standard Gaussian entries, a structured random matrix with the Hartley transform and a structure random matrix with the discrete cosine transform. We use blocking as in our main algorithm for each of the three types of projection matrices on all 64 cores of the distributed computing structure. The knot based approach has no obvious method of blocking or for it to be used on a grid computing environment. The MCMC algorithm in each case is run for 10, 000 iterations, with the first 2, 000 discarded as burn-in. Usual diagnostics do not show evidence of non-convergence of the MCMC algorithms and we report the effective sample size for the smoothing parameter as obtained from the CODA package in R Plummer et al. (2006). Time taken by each method is computed and reported. The time calculations are made by recording the exact CPU times for iterations, excluding the time for Frobenius norm accuracy and conditioning number calculations and also excluding the time taken to set up parallel Matlab workers via matlabpool. We measure predictive accuracy, of the predictions in the hold out set, based on the relative mean squared error.
We tabulate the results from the experiment in the table 3. The performance of any of the projection based approaches is substantially better than the knot based PP in terms of any of the parameters reported, predictive accuracy, time taken or the effective sample size of the smoothing parameter. Time taken by the PP in case of this large dataset was several days, as compared to a few hours for the blocked projection approaches. The improvement in predictive accuracy can probably be attributed to much better conditioned stable linear systems being used inherently and hence improved inference. The discrete cosine transform and Hartley transform perform marginally better than the random Gaussian projection, however the structured random transforms show marked improvement in the time taken, possibly due to the faster matrix multiplies.
There is little to choose between the discrete cosine transform or the discrete Hartley transform.

Discussion
In this article we present a new method of blocked approximate inversion for large positive definite matrices, using subsampled random orthogonal transforms, motivated by several recent developments in random linear algebra. The blocking strategy comes from development of parallel QR algorithms for tall matrices and in our case tall matrices are obtained from large positive definite matrices as the first projection step of a low rank approximation. We have also presented a comprehensive set of theoretical results quantifying the decay of eigenvalues of covariance matrices, which in turns helps bound the errors from the low rank projections. The examples show marked improvement in terms of numerical stability, prediction accuracy and efficiency and should be routinely applicable in a wide variety of scenarios.
One interesting future direction is the investigation of several structured random projections taken together. One may break up a target rank into several parts, the first of which is spent on a random projection, and the remaining pieces are learnt through the direction of maximum accuracy. Both theoretical (finding projections in the directions of maximum accuracy, which may be the direction of maximum information descent) and empirical investigation of this issue is currently being investigated. Table 1: Table of condition numbers for the squared exponential kernel for different values of the smoothing parameter θ 2 versus truncation levels, for truncating to the best possible approximation according to the Eckart-Young theorem. The rows represent the levels of truncation and the columns, the values of the smoothing parameter.


extfieldsofstudy:
Computer Science
Biology
Medicine

s2fieldsofstudy:
Computer Science
Biology

text:
Manually curated database of rice proteins

‘Manually Curated Database of Rice Proteins’ (MCDRP) available at http://www.genomeindia.org/biocuration is a unique curated database based on published experimental data. Semantic integration of scientific data is essential to gain a higher level of understanding of biological systems. Since the majority of scientific data is available as published literature, text mining is an essential step before the data can be integrated and made available for computer-based search in various databases. However, text mining is a tedious exercise and thus, there is a large gap in the data available in curated databases and published literature. Moreover, data in an experiment can be perceived from several perspectives, which may not reflect in the text-based curation. In order to address such issues, we have demonstrated the feasibility of digitizing the experimental data itself by creating a database on rice proteins based on in-house developed data curation models. Using these models data of individual experiments have been digitized with the help of universal ontologies. Currently, the database has data for over 1800 rice proteins curated from >4000 different experiments of over 400 research articles. Since every aspect of the experiment such as gene name, plant type, tissue and developmental stage has been digitized, experimental data can be rapidly accessed and integrated.

INTRODUCTION
Decades of research on rice has generated several knowledge resources (1) such as genome sequences (2), annotations (3,4), transcript data (5,6), mutant resources (7,8) and germplasm collections (9,10) available through several excellent databases (4,6,7,(11)(12)(13)(14)(15)(16) or published literature, with a sole aim to understand every aspect of rice biology. Although, these databases have exhaustive data on rice, they do not precisely catalogue and integrate every detail of the experimental data published in scientific literature. Characterization of every genetic component in rice (17) requires seamless integration of data from various sources and thus, data curation has emerged as a major challenge (18,19). Bulk of scientific data are available as peer-reviewed articles accessible either through manual reading or through databases that mine the text of the articles (20,21). Extensive efforts have been done to develop efficient text-mining approaches (22)(23)(24) to enrich and integrate scientific data (25). Primarily, text-mining is done with a certain perspective and may not capture all the aspects of the data presented in the article. Moreover, text mining is a slow post-publication exercise, it is almost impossible to keep pace with the accumulation of published literature. Thus, there is a need to develop novel curation concepts which can deliver precise and fast integration of data (17). In the current study, we have demonstrated the feasibility of digitizing the experimental data instead of mining the text of the article. Digitization of experimental data itself offers several advantages: (1) Digitization would render the experimental data amenable to computerized search such that exact data from multiple articles can be searched/retrieved rapidly. (2) Digitization into a standard computer readable format would facilitate a seamless and semantic integration of data. (3) Data models for digitization of experimental data may be extended further to enable pre-publication integration of data. Thus, from the very beginning, the data would be in a format where it can be integrated into any database with minimal manual intervention. Such approach would help close the gap between the curated databases and the accumulated scientific literature.
*To whom correspondence should be addressed. Tel: +91 9811574152; Fax: +91 11 24119430; Email: saurabh@genomeindia.org The authors wish it to be known that, in their opinion, the first three authors should be regarded as Joint First Authors.

Digitization of the experimental data via manual curation: concept
The manual curation workflow adopted in developing the database is designed to digitize the actual experimental data published in peer-reviewed articles. Such experimental data are mostly presented as images or graphs and is thus obviously not suitable for any computerized search. Moreover, there are no universal data representation standards for such experimental data. Thus, we developed the data curation workflow that was guided by the need to digitize experimental data as well as to lay foundation for development of standard representation of experimental data. The basic concept is depicted in Figure 1a and is dealt in detail elsewhere (Raghuvanshi et al. 2013, manuscript under preparation). In simple terms, every data point in a graph or image represents data of several dimensions such as gene name, plant type, tissue and growth conditions. Over the years, several ontologies and standard notations have been developed to represent at least some of these data types (26)(27)(28). Thus, a combined systematic use of these notations can represent information contained in every data point such that every data point is represented by a collection of these predefined terms. Since the notations would be alphanumeric in nature with pre-defined definition, the data can now be stored in any relational database and can be easily searched and correlated. Based on this concept, we have extensively used existing ontologies as well as in-house developed notations. To represent rice gene/protein, 'Rice Genome Annotation Project (RGAP)' locus ids have been primarily used (29). In cases where it was not possible to assign a RGAP id, the GenBank id has been used. Moreover, in several studies rice gene is studied in a heterologous plant system. In such cases, the experimental data may represent the expression of the heterologous gene. Thus, if the gene is from Arabidopsis, the 'The Arabidopsis Information Resource' (TAIR) (30) gene ids are used otherwise GenBank id of the heterologous gene has been used. The tissue and plant developmental stage is denoted by 'Plant Ontology (PO)' terms (31). Distinct PO terms have been used to represent plant developmental stage and plant tissue used in the analysis. In several cases, more than one PO term had to be used to represent a particular tissue such as 'PO:0025034-PO:0005352-PO:0000074' for 'leaf-xylem-parenchyma'. The growth conditions such as temperature, water status, light or presence of any chemical such as hormones, metals, etc. have been denoted by 'Environmental Ontology (EO)' terms (27). In order to accurately denote environmental conditions, the 'EO' term is appended with additional data such as time of treatment, concentration, etc. The 'Trait Ontology' terms have been used to represent phenotypic or biochemical traits while the functional aspect of the protein is depicted with the help of 'Gene Ontology' (32). It may be noted that due to the vast complexity of the published experimental data, none of the above mentioned ontologies were sufficient to encode every aspect of data. Thus, a large number of new terms had to be defined for most ontologies to efficiently digitize the experimental data. The ids of all new ontology terms start with the numerical digit '1' (such as TO:1000413-ascorbic acid content) to differentiate them from the pre-existing terms.
Besides the use of existing ontologies several other notations had to be developed to represent data for which no ontology exists. One of the most important is the notation to represent various experimental techniques. There are numerous different experimental techniques that are used in biological sciences and it is important to be aware of the experimental technique before any inference can be drawn from the data. Moreover, since all different type of experimental data is being coded with the help of similar set of ontologies and notations, it is very important to record the nature of the experiment. For example, a data point represented by '[LOC_Os05g46480/ PO:0009005 (roots)]' would either represent the transcript levels of LOC_Os05g46480 in roots in a RT-PCR experiment or protein levels in a Western analysis. Thus, the context in which each term is used is essential for interpretation of the data. Notations were also defined to represent the type of rice plant that was used in the experiment, i.e. whether it was a wild type, transgenic over-expression, transgenic-RNAi line, mutant, etc. as well as to represent the promoters used to drive the transgene.

The current release of the database
The current release of the database catalogues data for over 1800 rice genes from 400 peer-reviewed research articles (Figure 1b). Data from more than 4000 different experiments have been digitized using in-house developed manual data curation models. These experiments are based on a total of 140 different experimental techniques mostly related to gene expression analysis, biochemical activity analysis, protein-protein interaction, DNAprotein interaction and cellular/sub-cellular localization. Majority of the curated articles are from the year 2007 onwards and published in a variety of peer-reviewed plant biology journals. Figure 1c shows the distribution of the 'top five' Environmental Ontology, Gene Ontology (biological process) and Trait Ontology terms. Based on the curated experimental data, these genes have been associated to more than 600 different traits and over 300 Gene Ontology terms. The most frequent traits include physical traits such as 'survival rate', 'plant height' and 'root length'. Since there is a large collection of research articles on rice and the curation process is pretty extensive, we have tried to concentrate on studies related to abiotic stress in rice. This is clear from the distribution of the most frequent Gene Ontology terms that describe response to various abiotic stress conditions. The database also has protein-protein interaction data for 199 rice proteins as well as DNA-protein interaction data for 51 rice proteins.
How to access the database?
Every aspect (gene id, tissue, growth conditions, etc.) of the experimental data has been encoded with the help of defined notations (ontologies, etc.) and thus it is possible to retrieve the same set of data from several different perspectives. The database can be accessed either by browsing or searching for a particular keyword or term (Figure 2). A global overview of the data from different perspectives can be acquired by 'browsing' the database. Specific queries can be made by searching the database with the help of any of the ontology terms or with the help of a keyword. Both basic and advanced search options are available.

Browsing
The entire content of the database is summarized from six different aspects for browsing. The data can be browsed by either PubMed id, rice gene locus id, plant tissue/developmental stage (Plant Ontology term), growth conditions (Environmental Ontology), phenotypic or biochemical trait (Trait Ontology) and gene function/ localization (Gene Ontology term). While browsing by 'PubMed id', a list of all the articles curated in the database along with basic information such as title, authors, journal as well as the abstract is shown (Figure 3a). This list can be sorted either by PubMed id, journal or year of publication. A search option based on 'PubMed id', journal, author, year as well as any key word that may appear in the title of the article is also provided. On 'clicking' the PubMed id of interest, a list of all the rice gene loci that have been studied in that article is shown followed by a list of experiments that have been digitized (Figure 3b). Further, selecting any rice gene locus opens the 'Rice gene details' page that summarizes the information about the rice gene as presented in the selected article ( Figure 4). Initially, the data shown would be restricted to the selected article, however, it is also possible to access data for the rice gene locus gathered from all the articles that have been curated in the database by selecting the option 'Search entire database for LOC . . . .'. The information in the 'Rice gene details' page is divided into several sections. The section titled 'Basic information' presents the information about the protein domain as available in the Pfam database (33). The 'Functional details' sections list all the GO terms that have been mapped to the rice gene based on the experimental data. 'Clicking' on any of the GO term id provides the details of the experiment on the basis of which the GO term has been assigned. GO terms with (RGAP) tag have been acquired from the RGAP database (30). Similarly, the 'Plant developmental stage/ tissue details' section lists all the tissue or developmental stages where the selected rice gene has been studied. The presence or absence of expression/protein activity in a particular tissue/developmental stage is indicated by '(+)' or '(-)' sign, respectively. Thus, '(-)' means that although expression level or gene activity of the gene has been studied in that particular tissue but no detectable expression or activity was found. Selecting any of the PO terms presents the digitized data of the related experiment. In most cases, the data pertaining to a gene in a particular PO would only be part of the total experiment (e.g. one bar of the entire real-time RT-PCR graph). In such cases, initially data for only the selected data point would be shown, however, one can see the data for the entire experiment by selecting the 'Show data for entire experiment' option. The information about the environmental conditions under which the expression and activity of the gene has been studied is given as a list of Environmental Ontology terms in the 'Environment details' section. Thus, for example the affect of abiotic stress conditions such as drought, heat, salinity or hormones and other chemicals is recorded with the help of the EO terms. Selecting any particular EO term shows the digitized experimental data. Again 'Show data for entire experiment' option can be used if one wants to see the data of the entire experiment. The section entitled 'Trait details' summarizes all the phenotypic or biochemical traits that have been associated with the selected gene. The digitized experimental details can be accessed as explained above in the case of 'Plant developmental stage/tissue details'. The last section summarizes the physical interaction data (protein-protein and protein-DNA) for the rice protein.
The database can also be browsed through rice gene ids. This gives a list of all the rice genes that have been curated in the database. The first column of the table contains rice gene id that has been curated which in most cases is the RGAP locus id. A significant number of papers do not cite the relevant RGAP locus id. In such cases, RGAP locus id was assigned on the basis of sequence alignment at DNA level. If there was 100% identity and coverage between the gene sequence (mostly GenBank id) mentioned in the article and any RGAP locus id then the data were encoded using RGAP locus id. In case of lower but significant similarity the data were encoded using the original id (or the corresponding protein id, in case the mentioned id is of a genomic sequence or cDNA) as mentioned in the article and the closest RGAP locus id is mentioned in the column titled 'Nearest RGAP protein id'. The column 'Associated proteins' gives information of all the rice genes that have been experimentally shown to be associated with the concerned rice protein. Rice proteins having physical association are indicated with tag '(I)' whereas proteins that may affect the expression are indicated by tag '(R)'. If a protein affects the function of the concerned protein then it is indicated by tag '(F)'. Selecting any gene id in the first column opens the 'Rice gene details' which has been described as above. Similarly, the database can also be browsed by Plant Ontology term. This gives the list of all the Plant Ontology terms that have been curated in the database. Selecting any one term lists all the rice genes that have been studied in the selected plant tissue or developmental stage. Further selecting the rice gene id opens the 'Rice gene details' page. Browsing by Environmental Ontology, Trait Ontology and Gene Ontology follows the similar trend.

Searching the database
The database can be searched on the basis of gene id, plant developmental stage, plant part/tissue, environmental conditions, associated traits, associated molecular function, associated biological process or cellular localization. The search can be done by either specifying the exact 'term' id or the keyword. In order to facilitate this the search function operates in two tiers. If an exact 'term id' is defined (such as protein id, ontology term id) the relevant results are shown immediately. In case a keyword is given, in the first stage a list of all the related terms is displayed. The user may then select one or more of these terms and then proceed to the second stage where data relevant to the selected term is shown. It may be noted that search would be much faster when an exact id is specified such as 'PO:OOO9005' instead of the keyword 'root'. It is also possible to formulate search with multiple terms and boolean operators AND and OR. The output of the search is a list of rice gene ids as well as the link to the exact experiment (PubMed id and experiment no.) where the search term has been used.

Major highlights
Digitization of experimental data of published peerreviewed studies is a relatively less explored dimension of data curation. Nevertheless, it offers several advantages. Some of the major aspects that stand out are as follows: (a) Seamless semantic integration of data: The digitization of the experimental data has been done using the same basic elements (i.e. ontologies and notations) thus providing a natural connectivity of data coming from different experimental set-ups. (b) Universality of the basic concept: The fundamental models developed to digitize the experimental data on rice genes are universal in nature and can be easily adopted for other organisms as well. The basic requirement is the availability of organism specific ontologies. The current exercise has led to a significant enrichment of plant ontologies as a large number of novel terms had to be defined in order to digitize the published experimental data. (c) In-depth and precise functional annotation: One of the major advantages of indexing the experimental data itself is the preciseness of the functional annotation of the genes. The assigned GO terms lie very low in the hierarchy and are thus more precise in nature.
There is a direct connection of the associated 'GO' term with the experimental data that was used to assign the GO term. Another aspect is the cataloguing of the dependencies of a particular GO term assigned to the protein. For example a 'kinase' activity of the protein may be dependent on the presence of other protein, certain environmental conditions or restricted to a plant tissue. Such dependencies have also been catalogued. (d) Learning resource: Such databases can be great learning resources for young researchers. Naı¨ve questions such as 'What concentration of a particular chemical (hormone/metal ion/drug, etc.) should be used while designing experiments?' A simple search with a suitable 'EO' term would line up all the experiments conducted with the concerned chemical in a matter of seconds.

Future perspective
The database will be updated every 6 months and as the curation progresses the data would be enriched by inclusion of studies from diverse aspects such as biotic stress, yield, etc. The ultimate aim is to digitize all the relevant articles but since the task is enormous we are proceeding in a methodical manner so that important aspects such stress biology (abiotic or biotic), yield are better represented. The concept of digitization of the experimental data from peer-reviewed research articles is relatively new and is bound to evolve. Better user-friendly portals are being developed to facilitate easy and fast digitization of experimental data. Such portals would also enable active participation of third party curators to meet the challenge of curating all the articles. Lab data management implementation may integrate these curation models so that they can be used not only to store data but also in designing of the experiment from the very beginning. Thus, at the time of publication, the data would already be in a format compatible for integration in the database.

FUNDING
Funding for open access charge: The Department of Biotechnology, Government of India, India. The work is supported by grant from The Department of Biotechnology, Government of India, India.