Table of Contents Author Guidelines Submit a Manuscript
International Journal of Reconfigurable Computing
Volume 2010 (2010), Article ID 313479, 14 pages
Research Article

Reconfigurable Hardware Implementation of a Multivariate Polynomial Interpolation Algorithm

1Department of Computer Science, University of Puerto Rico, Río Piedras, PR 00924, Puerto Rico
2Department of Mathematical Sciences, University of Puerto Rico, Mayagüez, PR 00681, Puerto Rico

Received 2 March 2010; Revised 26 July 2010; Accepted 26 October 2010

Academic Editor: Viktor K. Prasanna

Copyright © 2010 Rafael A. Arce-Nazario et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Multivariate polynomial interpolation is a key computation in many areas of science and engineering and, in our case, is crucial for the solution of the reverse engineering of genetic networks modeled by finite fields. Faster implementations of such algorithms are needed to cope with the increasing quantity and complexity of genetic data. We present a new algorithm based on Lagrange interpolation for multivariate polynomials that not only identifies redundant variables in the data and generates polynomials containing only nonredundant variables, but also computes exclusively on a reduced data set. Implementation of this algorithm to FPGA led us to identify a systolic array-based architecture useful for performing three interpolation subtasks: Boolean cover, distinctness, and polynomial addition. We present a generalization of these tasks that simplifies their mapping to the systolic array, and control and storage considerations to guarantee correct results for input sequences longer than the array. The subtasks were modeled and implemented to FPGA using the proposed architecture, then used as building blocks to implement the rest of the algorithm. Speedups up to and were obtained for the subtasks and complete application, respectively, when compared to a software implementation, while achieving moderate resource utilization.

1. Introduction

Recent years have seen a significant increase in methods and tools to collect genetic data from which important information can be extracted using a number of techniques [1]. For instance, microarray data collected at various steps in organism development can help geneticists understand its developmental process and response to environmental stimuli [2]. Several models such as ordinary differential equation models [3], continuous models [4], stochastic models [5], and discrete models, of which Boolean models [6] are special cases, have been proposed. Essentially, each node represents the expression level of a gene (or genes) of interest at a time point. A directed edge from node to node symbolizes the influence of gene(s) in to the gene(s) in .

Our research group focuses on multivariate finite field gene network (MFFGN) models, in which multiple genes are monitored at each time step and their expression levels are discretized to a predefined set of values , where is a prime number, that is, the expression level of each gene is an element of the prime field [7, 8]. One way to deduce gene interaction from these models is to solve the reverse engineering problem, that is, find functions that describe the network's state transitions. An important step in the solution of this problem is determining an interpolating polynomial for the points by using such methods as a multivariate version of Lagrange's interpolation formula [8]. Polynomial interpolation over finite fields also has applications in error correcting codes and is a major building block of many numerical methods [9, 10].

Fast algorithms and implementations are needed to sustain rapidly increasing bioinformatics computational demands [11]. Reconfigurable logic represents an attractive platform for the high performance implementation of finite field algorithms, with many designs achieving significant speedups over their software equivalents. However, despite the availability of tools to ease the hardware design flow, there is still a notable gap between the bioinformatics and hardware experts [12]. Our group developed a reconfigurable logic implementation of a multivariate polynomial interpolation algorithm suitable for reverse engineering in genetic networks. Rather than striving for a specialized punctual design, we used this opportunity to identify and develop common structures that might be of use to other researchers and applications. Well-documented and parameterizable components that perform functionalities that are common in many algorithms will alleviate the effort spent on new implementations.

In this paper, we discuss the interpolation algorithm and our proposed architecture. We emphasize several computational substructures that appear repeatedly throughout the design and which can be implemented effectively using a hardware systolic array-based architecture. These tasks are of the type where we perform a certain reduction or rearrangement of a sequence of elements from multivariate polynomials or boolean expressions. The systolic array concurrently manages data receipt and parallel processing, making it an amenable structure when dealing with streamed data. The simplicity of the array cells, storage, and control unit, allows the instantiation of multiple cells while maintaining competitive clock frequencies, thus achieving high performance. Several tasks critical to interpolation were modeled and implemented to FPGA using the proposed architecture, obtaining speedups up to 172× when compared to a software implementation, while achieving low resource utilization. These implementations were used as components to develop a complete, high-performance multivariate polynomial interpolation methodology on hardware.

Paper Outline
Section 2 discusses the interpolation methodology, while Section 3 details data representation for the implementation. Section 4 describes the hardware implementation in general. Sections 5 and 6 describe an algorithmic generalization for the reduction tasks and present the proposed hardware architecture. Section 7 explains the rest of the implementation blocks. Section 8 reports the results of FPGA implementations, both of reduction tasks and the complete implementation. Section 9 provides our conclusions.

2. Interpolation Methodology

Discrete models, in particular finite field models, have been proposed for regulatory processes such as genetic networks [7] and biochemical networks [13]. In the setting of a finite field model a genetic network can be seen as a finite system whose states are governed by the iterations of a tuple of polynomials , where is the number of genes.

Given a sequence of   -tuples of values from a set , representing the states of genes at times through , the reverse engineering problem seeks to find a function for which , .

When is a finite field, function can be represented by a tuple of polynomials where for each . As pointed out in [8] such polynomials can be found by the following variation of Lagrange's interpolation formula: where each denotes an -tuple and denotes the th component of . Given and , where , is the first component in which and differ.

In the context of genetic networks modeled by polynomials over a finite field, the resulting polynomials give a sense to the biologist about the interactions between genes. Nevertheless, results from (1) may contain redundant variables even after algebraic simplifications.

Example 1. Let us consider the function which is incompletely specified by Table 1.
Using (1), an interpolating polynomial for is given by . This polynomial depends on ,, and and cannot be further simplified by any algebraic means. However, the polynomial also interpolates at the given points but depends only on variables and . In this context, is a redundant variable, since we found another interpolating polynomial for that does not include in its expression.

Table 1: Table for .

Redundant variables are undesirable as they introduce complexity into the polynomials without adding information valuable to the biologist. Furthermore, empirical data suggest that genetic networks are sparsely connected [14]. Redundant variables can be identified using Sasao's algorithm for the detection of dependent variables in incompletely specified multiple valued functions [15].

2.1. Essential Variables and Bases

For any set and any positive integer , we denote each by . Let be a set of variables. For any set of variables and any , we define . For instance, if , , and , the projection of on , is .

Let and where and are integers, and let , where . We say that a variable is essential in or depends on if there exist such that for all , and . A set of variables is a basis for if (1) for every , if and only if and (2) there is no other set of variables that properly contains with this property.

Lemma 1. If is an essential variable and is a basis, then .

Proof. Let be an essential variable. Then there exists and such that , and . Let be any basis. Then , but . Hence .

For any function and any basis , can be expressed solely in terms of the variables of . For any basis , any variable is redundant (with respect to .)

Our goal is to determine interpolating polynomials over finite fields in terms of the variables of any of its bases. To this end, let us first review an algorithm of Sasao [15] which determines all bases for any multiple valued function, not just those over finite fields.

Let , where . For each , let .

The set of variables appearing in each disjunct of constitutes a basis.

Thus, the algorithm consists of two stages. First, determine , which requires steps and second, convert to disjunctive normal form, which requires steps.

The following lemma, whose proof is immediate, is useful in simplifying expressions such as and in the computation of .

Lemma 2. Let be a subset of the variables , let be the disjunction of the variables in and let be a disjunction of some subset of . Also let be the conjunction of some subset of . Then (a), (b).

In what follows we abbreviate by .

Example 2. Let be a function whose values are included in Table 2.
Applying Algorithm 1 and the reduction rules of Lemma 2, we obtain Thus, the bases are , , and .

Table 2

Algorithm 1: Algorithm to determine the collection of bases.

2.2. Interpolating Polynomials with No Redundant Variables

Given any partially defined function over a finite field and any basis for , we would like to determine interpolating polynomials involving only those variables of , that is, polynomials containing no redundant variables.

Definition 1. Let where and let be a basis for . For all define the relation if and only if , for each , , and are defined, and .

In other words, if and only if . It is straightforward to verify the following.

Lemma 3. is an equivalence relation.

Let be a set of representatives of each equivalence class under .

Theorem 1. Let be a field, let , , , be a set of points in , and let be a basis for the function defined by , . Then an interpolating polynomial for whose only variables are those of is where is the smallest index such that and and differ.

Proof. For each , , let be the representative of the equivalence class containing . Then and so , where is the value of each tuple of .

Note that for each of the factors in (3), and and so for all and all

Equation (3) gives us a new algorithm for polynomial interpolation that involves only those variables in the chosen basis , and furthermore, it also avoids redundant tuples by computing only on the representatives of the equivalence classes given by the relation .

Using (3), operations are needed to interpolate points for each function . Given points of a function , the bases for the component functions need not be the same. In this case, we can apply (3) times, thus giving an algorithm with complexity .

Example 3. Let us assume that the of Example 2 represent the corresponding elements of . We found that the bases are , , and . Now let us use the new algorithm given by (3) to find an interpolating polynomial over in terms of the basis . The equivalence classes under are and the set of class representatives is
Applying (3) and making use of (6), we have
Thus, a polynomial which depends only on the variables , , and and interpolates the tuples given by Example 2 is

3. Data Representation

As can be deduced from the previous section, an implementation of the multivariate polynomial interpolation algorithm must be able to represent and compute on disjunctive normal form (DNF), conjunctive normal form (CNF), and multivariate polynomial expressions. The following subsections establish our method of representation, which ultimately determines the techniques used for implementing the algorithms.

3.1. DNF and CNF Expressions

A Boolean expression in DNF using only uncomplemented variables is a disjunction of conjunctive clauses , where , . We shall refer to each term as a . For example, the following is a Boolean expression in DNF: . In our scheme, a DNF expression is represented as a set of cubes where each is represented as a sequence . Observe that complemented variables are not used in the algorithm for finding essential variables and bases, hence their representation is not necessary in this context.

Example 4. The DNF expression is represented as

A Boolean expression in CNF using only uncomplemented variables is represented as a set of disjunctions of literals (DOL) where each DOL is represented as a sequence . Within our implementation there is no need to distinguish the CNF from the DNF representation since no individual component needs to operate on both simultaneously.

Example 5. The CNF expression is represented as

3.2. Polynomial Representation

A multivariate polynomial over , where is prime is defined as , where and . For example, one multivariate polynomial over is . In our scheme, a polynomial is represented as a set of monomials where each is represented as a tuple . Binary representation requires bits for each and .

Example 6. The multivariate polynomial over is represented as . For the three variable polynomial over each monomial requires bits for binary representation.

4. General Description of Implementation

Our algorithm for determining interpolating polynomials that do not contain any redundant variables consists of two stages. The first stage identifies the dependent variables using Sasao's algorithm, while the second uses this information along with (3) to determine . Thus, a key part of the reconfigurable logic implementation is choosing efficient hardware structures to perform the critical parts of both stages. For the first stage, the most time-consuming task is the conversion of a Boolean expression in conjunctive normal form (CNF) to disjunctive normal form (DNF). The second step relies heavily on the identification of distinct binomials as they are generated from (3) and multivariate polynomial multiplication/addition.

To take advantage of the FPGA's fine-grained parallelism and considering their limitations in I/O bandwidth, these computations were implemented in a pipelined manner. In other words, computational blocks were designed to sustain stream processing as allowed by the FPGA's resources, rather than accumulate all needed data and then perform block processing.

Figure 1 shows a block diagram of the implementation. The CNF Generator receives the -tuples and performs Algorithm 1 to generate the disjuntions of literals, that is, the terms. These are streamed to the Cover DOLs block which eliminates redundant disjuntions of literals by using Lemma 2. The nonredundant DOLs are then fed to block CNF2DNF which computes the conjunctions of literals and uses Lemma 2 to eliminate redundant conjunctions. The output of CNF2DNF is the set of bases. A priori biological knowledge about the organism under study is required to choose the most adequate base. The -tuples are also fed to the second stage which uses the chosen base to determine representatives of the equivalence classes from the -tuples. The representatives are then streamed to the modified Lagrange interpolation block which implements (3) by generating the binomials for each pair of class representatives (Binomial Generator), eliminates duplicate binomials (Filter Duplicates), then multiplies and adds the product results.

Figure 1: Block diagram of the implementation, highlighting the subtasks that are implemented using the systolic array architecture.

The highlighted blocks in Figure 1, that is, Cover DOLs, Cover COLs, Determine Class Representatives, and Polynomial Addition, have two characteristics that make them suitable candidates for systolic array processing: (1) each performs an operation on a sequence of elements that are streamed from the preceeding block, (2) the operation is such that it performs reduction of a sequence of elements from multivariate polynomials or Boolean expressions. For example, the Cover Cubes receives preliminary results from the CNF to DNF conversion and eliminates redundant cubes using Lemma 2. This is an operation that can be described as given a sequence of cubes () determine a subsequence such that it contains only cubes where there exists no other such that .

The following sections describe a generalized algorithm and systolic array-based architecture for the type of reduction operations common throughout our interpolation methodology. This is followed by a description of how they are utilized as part of the architectural data path.

5. Generalization of Reduction Tasks

The subtasks that deal with reduction in our interpolation algorithm can be generally described as follows. Given the sequence of elements from , compute the sequence where by eliminating or combining elements based on binary comparisons between them. The empty element () is introduced for the representation of operations where groups of elements can be reduced to a single element. Algorithm 2 shows a general form of these problems, where and are binary functions of the form , and the symbol denotes parallel computation.

Algorithm 2: Nested loop algorithm for reduction tasks.

The tasks of distinctness, polynomial addition, and redundant Boolean term elimination (hereon referred to as cover) can be implemented by specifying and , as explained in the following subsections.

5.1. Distinctness

The task of identifying the distinct elements of a sequence is needed within our interpolation algorithm to determine class representatives and filter out binomial duplicates. It can be implemented by Algorithm 2 with

Example 7. Assume . After computation with Algorithm 2 using (11), the result is .

We provide a correctness proof for this operation. Proofs of the other reduction operations given in this section are similar.

Lemma 4. For an input sequence , where , the output of Algorithm 2 with functions and defined as in (11) outputs where is the subsequence of distinct elements of the input sequence.

Proof. We show by induction that after iterations of the outer loop, where for each ,.
After iterations, . Next suppose that after iterations, for some where for each . Then either (1) or (2) . In the first case, since is not equal to any , and is not changed by the st iteration. In the second case, will be replaced by during the st iteration. Hence in either case, where for each , is either or an element , where .
Hence after iterations, and remains the same after iterations since .

5.2. Polynomial Addition

Assume that each of the elements in is a monomial. A polynomial addition operation can be implemented by Algorithm 2 with where , that is, the monomial without the coefficient.

Example 8. Assume . After computation with Algorithm 2 using (16), the result is .

5.3. Cover

Assume that each of the elements in is a cube. A cube covers another cube (represented as ) if all Boolean variables present in are present in (e.g., , but ). The cover operation can be implemented by Algorithm 2 with where for any .

Example 9. Assume . After computation with Algorithm 2 using (17), the result is . Notice that there may be duplicate terms in the result but these can be easily eliminated using a distinctness operation.

Algorithm 2 is also amenable to other problems that can be expressed as the result of binary comparisons between every two elements of a sequence. For example, the problem of sorting a sequence can be implemented with

6. Hardware Structure for Reduction Tasks

This section explains our proposed hardware structure and outlines the mapping for the reduction tasks. We first discuss the systolic array and cells with the assumption that the array is deep enough to process the input sequence without overflowing. We proceed by discussing the additional components that must be added to guarantee correct processing even if the overflow condition does not hold.

6.1. Systolic Array and Cells

The proposed structure is a linear systolic array of computational cells, each performing an operation deduced from and of the above discussion. Figure 2 illustrates the array and contents of the basic cell. Each cell contains two registers and , where each has enough resources to hold any single element from the alphabet . Each cell also contains the logic to perform functions and . All cells are initialized to the empty symbol. At each clock step a new element from the sequence is input to the first cell and every cell performs a comparison between its stored value and the value from the previous cell and assigns (possibly) new values to the registers according to and . After all elements have been processed the registers contain the resulting sequence and can be shifted out of the array using the connections between the registers.

Figure 2: Block diagram of systolic array.

More formally, our linear systolic array computation can be modeled as a simplified version of the generic systolic array presented in [16], as shown in Algorithm 3. The notation symbolizes the content of register in cell at iteration .

Algorithm 3: Model of the linear systolic array computation.

Figure 3 illustrates the operation of the array using cells that implement the distinctness task by defining and as in (11). Observe that the operations implied by and and the connections between neighboring cells accomplish the comparisons established in Algorithm 2, albeit in a different order and in a parallel manner. In fact, the computation implemented by the systolic array can be interpreted as a reindexing of Algorithm 2 as presented in Algorithm 4. Table 3 shows the comparison indexes generated by Algorithm 4 for the example in Figure 3. Notice that even though Algorithm 4 has reordered the indexes, the same data dependence is maintained as in Algorithm 2. In other words, if is performed before a (where ) in Algorithm 2 the same order is preserved in Algorithm 4. The same condition holds for the order of and .

Table 3: Illustration of index generation by Algorithm 4 for the example shown in Figure 3.

Algorithm 4: Reindexed algorithm for reduction tasks.

Figure 3: Systolic array implementing the distinctness operation for . Each row represents a cell. Within each cell, the left and right columns represent the and registers, respectively. Elements are input to the top row and flow downwards.
6.2. Processing Sequences Longer than the Array Depth

The functionality of the basic cells mandates their implementation on user logic inside the FPGA (rather than on embedded units, such as block RAMS). Depending on the application, the characteristics of the data sequence, and FPGA model, the FPGA might not have enough resources to instantiate a systolic array whose depth is greater than or equal to (the sequence length). In such cases, there is no guarantee that the array by itself will be able to compare each element against each other to obtain the correct solution. To illustrate this anomaly, assume that in the example presented in Figure 3 the input is a sequence of distinct elements . The array would be filled once is registered in the last cell, hence none of the cells would have operated on the combination of and . A control and temporary storage mechanism must be provided to reiterate the data through the pipeline if needed and guarantee that all elements are compared against each other. This can be accomplished by a scheme such as the one shown in Figure 4 through the addition of a FIFO, control unit and several datapath control elements (i.e., multiplexers).

Figure 4: Pipeline with added overflow FIFO and control unit to support sequences longer than .

The Overflow FIFO (OF) stores elements that make it through the array without having been eliminated or registered. These elements will have to be reiterated through the array to test their relation to other elements in the OF. For example, suppose we have an array of depth that implements the distinctness operation. The sequence will fill the pipeline with so that the final overflows the pipeline. These 3 elements are stored in the OF for a later pass through the array. The FIFO depth where is the number of elements in the sequence and is the array depth. Although both systolic array cells and FIFO spend mostly on register resources, a reason why one may be able to increase FIFO depth and not necessarily array depth is that in FPGAs FIFOs are easily mapped to the embedded block RAM, whereas processing cells map to user logic.

The control unit controls the data flow between the array and the OF and determines how to reiterate the overflowed data through the array. For the reduction tasks presented in Section 5 we can deduce three basic modes of operation for the control unit. Figures 5 and 6 illustrate the need for these modes. Figure 5 shows the end of a first pass of a sequence through systolic arrays that implement (a) a sorting algorithm, and (b) sum of polynomials.

Figure 5: Contents of the systolic array and OF after a first pass of the sequences through (a) sorting and (b) polynomial addition.
Figure 6: Contents of systolic array and overflow FIFO for the shown sequence (a) after first pass, (b) after second pass (OF to array refeed), (c) systolic array is emptied for new OF refeed, (d) after final refeed from OF. Dotted lines illustrate the paths of data refeed.
(i)In the case of a nonreducing operation like sorting (Figure 5(a)), the array will contain sorted elements. Thus, elements in the array can be shifted out while elements in the OF are reinputted into the array. This must be repeated until the OF is empty after a pass, that is, at most times.(ii)When adding polynomials (Figure 5(b)), the array will contain at most monomials with their final coefficients. Thus, elements in the array can be shifted out while elements in the OF are reinputted into the array. This must be repeated until the OF is empty after a pass, that is, at most times.

Figure 6 illustrates the various passes needed for the cover operation. At the end of a first pass the array will contain at most cubes. However, in this case a first pass through the array does not guarantee that all cubes in the OF have been compared to all the cubes registered in the array. For example, in Figure 6(a) cube passed through the array without being registered. Meanwhile, the last cube covered and in two cells of the array. Thus, we have at least one cube () in the OF which was not compared to a cube in the array (). The elements in the OF must be refed through the (unflushed) array to allow every stored in the FIFO to be compared against every registered in the array. The refeed is repeated until a complete pass of the OF cubes through the array does not modify the cubes registered in the array. After this (Figure 6(b)), if any in the OF is going to be part of the result, for example, , it must only be still compared to other elements in the OF. The contents of the array may be shifted out to a memory that stores the results, the array is reset and the OF cubes are refed (Figure 6(c)). The three-step process continues until the OF is empty after an OF refeed.

7. Further Implementation Considerations

Aside from the reduction tasks, the computational blocks in our implementation can be classified as performing one of two functionalities: element pair generation or multiple term multiplication. Element pair generators, such as the CNF Generator and the Binomial Generator in Figure 1, input a set of tuples and output a distinct pair of tuples at each computational step. This is accomplished by using the architecture illustrated in Figure 7. As tuples are input (each accompanied by a load signal) both RAMs store the tuples and an Address Control Generation unit counts the number of tuples. When a status signal indicates that all tuples have been received, the ACG begins generating all possible address combinations for , . The distinct pairs of -tuples produced by the CNF Generator are fed to   -comparators to compute the terms in Algorithm 1, as shown Figure 8. In the case of the Binomial Generator, the corresponding terms within the distinct pairs of -tuples are subtracted from each other and the results are used to determine the terms of the binomial for (3), as shown in Figure 9.

Figure 7: Element pair generator.
Figure 8: Computation of the disjunction of literals from a pair of -tuples.
Figure 9: Computation of the terms of a binomial for (3) from a pair of -tuples.

Multiple-term multiplication/conjunction, such as needed in blocks CNF2DNF and Lagrange, is performed using the architectures shown in Figures 10 and 11, respectively. The architecture depicted in Figure 10 receives DOLs from the CNF Generator, the DOL2Literals block generates a literal expression for each found in the DOL. The literals are queued in a circular FIFO and are conjuncted to each of the results in the Preliminary Results FIFO. After each DOL is conjuncted, the resulting cubes are passed through the Cover Cubes block which eliminates redundant cubes. A control unit monitors and generates signals to coordinate the data path activities.

Figure 10: Illustration of a step in the conversion of CNF to DNF. The first two conjunctions have been computed and their redundant terms eliminated and stored in the Preliminary Results FIFO.
Figure 11: Illustration of a step in the multiplication of binomials in . The multiplication of the first two binomials has completed and the results are stored in the Preliminary Results FIFO.

The architecture depicted in Figure 11 performs the computations of (3). Each binomial that is received from the Binomial Generator is registered. Then each of its two monomials is multiplied by each of the monomials that has resulted from previous binomial multiplications. The monomials that result from the multiplication are input to the PrelimPolyAdd block which adds similar monomials to reduce the size of the preliminary product. PrelimPolyAdd outputs its result to the PrelimRes FIFO. The iterative process continues until the last binomial of a product is received, in which case each monomial of the preliminary product is output to the final Polynomial Addition block (see Figure 1). Once the last binomial of the last product has been multiplied, the (final) Polynomial Addition block is signaled and begins to output the monomials of the final result.

8. Results and Discussion

This section presents and discusses results for the individual reduction tasks as well as the complete interpolation algorithm implementation.

8.1. Reduction Tasks

The systolic array cells for the subtasks of distinctness, polynomial addition, and Boolean cover were modeled using Verilog HDL, based on their function definitions presented in Section 5. The array and control units were also modeled using Verilog based on the descriptions provided in Section 6. Cell, systolic array and control unit models were designed as parameterizable modules in terms of width of elements, array depth and control mode, respectively. Models were behaviorally simulated using ModelSim SE 6.0 to validate their results against software versions of the algorithms. Xilinx ISE Design Suite 11.1 was used to synthesize the architectures for a Virtex-4 XC4VLX200ff1513-11 in the DRC Dev Systems DS2000. Default synthesis parameters were used, such as mapping FIFOs to block RAMs and speed as the optimization goal. All the reported results for the individual tasks are before Place and Route.

For each of the three subtasks, several systolic array depths () were implemented, then randomly generated sequences of various bit widths and lengths were input. Table 2 shows timing and resource results for (1) the cover operation on sequences of length 4096 and width of 32 bits, (2) the polynomial addition operation on sequences of length 4096, where each element represents a monomial of 8 variables in , and (3) the distinctness operation on the same sequences as (3). Results are compared against a software implementation compiled using the GNU project c++ compiler with optimization level 3 and running on a 3.40 GHz Intel Pentium D CPU with 2 MB cache and 3 GB RAM. CPU execution times were measured by accessing the processor’s cycle counter. In all cases, we report the smallest measured CPU execution time. The freq column indicates the FPGA clock frequency, and are the FPGA and CPU run times, and speedup is computed as . The numbers shown in parenthesis next to each resource quantity are the resource percentages for the target device.

Several important observations can be highlighted from the results in Table 4. First, the maximum operation frequency is independent of the array depth, allowing performance to scale adequately with depth. In fact, given the simplicity of the required control unit and storage, the longest path was mostly determined by the implementation of functions and inside the basic cells. Second, the parallelism achieved through the use of the systolic array greatly compensates for the refeeds that must be performed when the result is longer than the array depth. For instance, the result for the cover operation was of length 3386 on average, yet there is a speedup of more than 40× for . Finally, the resource utilization to achieve competitive speedups is minimal and easy to estimate when scaling since it is almost completely attributable to the systolic array cells and overflow FIFO. Hence when using this structure as part of a bigger design, the designer could readily determine which parameters are best suited for the area/performance constraints and objectives.

Table 4: Experimental results for subtask blocks.

Part of our purpose here is to argue that the reduction tasks can be implemented efficiently using the systolic array design. With Table 4, we intend to demonstrate the effect of each particular operation for users that might want to use any of them as part of a bigger design. Clearly, the ultimate speedup may be dictated by other components of the final design, as shown by Tables 5 and 6.

Table 5: Experimental results for stage 1.
Table 6: Experimental results for stage 2.
8.2. Interpolation Algorithm

We modeled the proposed interpolation algorithm in Verilog HDL using the reduction tasks and the architectures described in previous sections as building blocks. The behavioral simulation and synthesis tools as well as the software compilation parameters and platform were the same as in the reduction task experiments.

The results are shown separately for each of the two interpolation stages since, once the bases have been computed by the first stage, a scientist's intervention would be necessary to choose the most appropriate before proceeding to the second stage. For both stages, the synthesis tool determined maximum operating frequencies above 150 MHz, thus we chose 150 MHz for the experiments. The Xilinx ISE 11.1 place and route tools (with default effort level) were able to meet the timing requirements reported for Stages 1 and 2. We attribute this to the fact that the great majority of connections in our design are local and grow only linearly with .

8.2.1. Stage 1

Table 5 reports timing and resource results for the first stage. Randomly generated sets of 100 -tuples () were input to implementations with diverse systolic array depths. The table shows results for the case with the shortest execution time for each . The smaller values of generated fewer intermediate terms during computation, hence smaller values of produced better results since they incurred in less unnecessary overhead for filling and emptying. As the number of variables increases, the number of terms produced by Algorithm 1 increases exponentially, which is clearly evidenced by the run times, that is, columns and . The FPGA implementation achieves speedups of 13× to 67×, which mostly increase with . The increase in speedup at the higher values of can be attributed to a higher percentage of time of full pipeline utilization as well as the increased benefit of parallelism in the comparisons between the -tuple variables, that is, the computation of the terms in Algorithm 1. As expected, the use of necessarily serial processing tasks such as the multiplication of the DOLs to obtain the DNF (see Example 2) limit the speedup achieved by the FPGA implementation when compared to the results obtained using the cover subtask by itself.

Although the ultimate goal of reverse engineering might be understanding complete complex biological systems, recent reported models limit their exploration over specific subsystems or punctual mechanisms, for example, the study cell-cycle regulatory network of fission yeast and apoptosis (programmed cell death) [17, 18], using a moderate number of genes between 8 and 20. Our results lead us to believe that the advantages of using the proposed architecture in reconfigurable logic will be sustained even at bigger sizes. As the sizes increase our architecture will have to rely increasingly on iterations rather than parallel computations. Nevertheless, due to the custom nature (finite field) of the data, even at these sizes the parallelism over each iteration will be enough to offer significant speedups over the alternative, fully serial implementation in general purpose processors.

Speedup is maintained at a cost of increasing resource utilization. For the considered cases, BRAMs, which are used to implement the various FIFOs and RAMs of stage 1, is the fastest growing resource. This could pose a challenge to maintain performance at even higher values. However, this issue can potentially be resolved by complementing the BRAM capacity with a dedicated external memory bank, as is commonly included in modern reconfigurable computer platforms [19, 20]. A study of the impact of external memory on the performance of the proposed architecture is future work. An advantage of our proposed architecture is that all the required memories behave as FIFOs. Thus, their accesses patterns are predictable and (as part of future work) we can use this information to device a buffering scheme that hides the latency of accessing the external memories.

The increase in user logic (Slices, LUTs) in higher can be circumvented by using shallower systolic arrays (smaller ), which still maintain competitive speedups. For example, although not shown in the table, for speedups of 30× and 45× were obtained by using values of 64 and 128.

8.2.2. Stage 2

As explained in Section 2.2, our algorithm uses a base computed in Stage 1 to determine a set of equivalence classes from the original set of -tuples. The equivalence classes are then used to compute an interpolating polynomial using (3). A greater number of variables in the chosen and a high variability among the -tuples translates to larger cardinalities of . The greater the cardinality of the higher the number of binomials computed by (3), which constitutes the bulk of computation for Stage 2.

The timing results shown in Table 6 support the previous statements. The sets of 100 -tuples from Stage 1 along with bases of 2 to 4 variables were input to Stage 2. Observe that since the data sets were randomly generated there was high variability in the sets of -tuples, which resulted in similar cardinalities (column Class Reps) and run times (columns and ) for the various sizes of (for a given number of variables in the chosen basis). The FPGA implementation achieves speedups of 16× to 25×, which mostly increase with . The moderate increment in speedup as increases is attributed to the increased benefit of parallelism for the generation of binomials, and the comparison and addition of monomials.

Resource utilization distribution in Stage 2 is roughly uniform between the user logic, that is, Slices, and BRAMs. Furthermore, the increment in resource utilization for increasing values of and numbers of variables in the basis is quite moderate, for example, the percentage of slices goes from 8.69% for to 12.87% for .

9. Conclusion

This paper presents a new methodology based on Lagrange interpolation with two important properties: (1) it identifies redundant variables and generates polynomials containing only nonredundant variables, and (2) it computes exclusively on a reduced data set. The analysis of the methodology for its hardware implementation led us to the identification of several reduction tasks which were generalized to a simple algorithm. The generalized algorithm can be efficiently mapped to a systolic array in which each processing cell implements a pair of binary operations between an incoming and a stored value. The tasks of Boolean cover, distinctness, and multivariate polynomial addition were implemented and served as building blocks to the rest of the application. The FPGA implementation of the reduction operations and the complete application achieved speedups of up to 172× and 67×, respectively, as compared to software implementations run on a contemporary CPU, with moderate resource utilization.


An earlier version of this paper appeared as “A systolic array based architecture for implementing multivariate polynomial interpolation tasks” in the Proceedings of the 2009 International Conference on ReConFigurable Computing and FPGAs (ReConFig’09) [21]. Dr. E. Orozco was partially supported by Grant no. P20RR016470 from the National Center for Research Resources (NCRR), a component of the National Institutes of Health (NIH). The authors thank Dr. Humberto Ortiz-Zuazaga for his helpful discussions on microarray technology and bioinformatics data trends.


  1. U. Muller and D. Nicolau, Microarray Technology and Its Applications, Springer, New York, NY, USA, 2004.
  2. P. A. Ortiz-Pineda, F. Ramírez-Gómez, J. Pérez-Ortiz et al., “Gene expression profiling of intestinal regeneration in the sea cucumber,” BMC Genomics, vol. 10, article 262, 2009. View at Publisher · View at Google Scholar · View at Scopus
  3. H. De Jong, “Modeling and simulation of genetic regulatory systems: a literature review,” Journal of Computational Biology, vol. 9, no. 1, pp. 67–103, 2002. View at Publisher · View at Google Scholar · View at Scopus
  4. J. F. Knabe, M. J. Schilstra, and C. L. Nehaniv, “Evolution and morphogenesis of differentiated multicellular organisms: autonomously generated diffusion gradients for positional information,” in Artificial Life XI : Proceedings of the 11th International Conference on the Simulation and Synthesis of Living Systems, pp. 321–328, MIT Press, 2008.
  5. W. J. Blake, M. Kærn, C. R. Cantor, and J. J. Collins, “Noise in eukaryotic gene expression,” Nature, vol. 422, no. 6932, pp. 633–637, 2003. View at Publisher · View at Google Scholar · View at Scopus
  6. S. A. Kauffman, The Origins of Order: Self-Organization and Selection in Evolution, Oxford University Press, Oxford, UK, 1993.
  7. D. Bollman, O. Colón-Reyes, and E. Orozco, “Fixed points in discrete models for regulatory genetic networks,” EURASIP Journal on Bioinformatics and Systems Biology, vol. 2007, Article ID 97356, 8 pages, 2007. View at Publisher · View at Google Scholar · View at Scopus
  8. R. Laubenbacher and B. Pareigis, “Equivalence Relations on Finite Dynamical Systems,” Advances in Applied Mathematics, vol. 26, no. 3, pp. 237–251, 2001. View at Publisher · View at Google Scholar · View at Scopus
  9. Z. Zilic and Z. G. Vranesic, “A deterministic multivariate interpolation algorithm for small finite fields,” IEEE Transactions on Computers, vol. 51, no. 9, pp. 1100–1105, 2002. View at Google Scholar · View at Scopus
  10. Y. Luo, “A local multivariate Lagrange interpolation method for constructing shape functions,” International Journal for Numerical Methods in Biomedical Engineering, vol. 26, no. 2, pp. 252–261, 2010. View at Publisher · View at Google Scholar
  11. M.-L. T. Lee, Analysis of Microarray Gene Expression Data, Kluwer Academic Publishers, Dordrecht, The Netherlands, 2004.
  12. K. Benkrid, Y. Liu, and A. Benkrid, “A highly parameterized and efficient FPGA-Based skeleton for pairwise biological sequence alignment,” IEEE Transactions on Very Large Scale Integration (VLSI) Systems, vol. 17, no. 4, pp. 561–570, 2009. View at Publisher · View at Google Scholar · View at Scopus
  13. A. S. Jarrah, R. Laubenbacher, B. Stigler, and M. Stillman, “Reverse-engineering of polynomial dynamical systems,” Advances in Applied Mathematics, vol. 39, no. 4, pp. 477–489, 2007. View at Publisher · View at Google Scholar · View at Scopus
  14. R. D. Leclerc, “Survival of the sparsest: robust gene networks are parsimonious,” Molecular Systems Biology, vol. 4, article 213, 2008. View at Publisher · View at Google Scholar · View at Scopus
  15. T. Sasao, “On the number of dependent variables for incompletely specified multiple-valued functions,” in Proceedings of the 30th IEEE International Symposium on Multiple-Valued Logic (ISMVL '2000), pp. 91–97, May 2000. View at Scopus
  16. E. P. Gribomont and V. V. Dongen, “Generic systolic arrays: a methodology for systolic design,” in Proceedings of the 4th International Joint Conference CAAP/FASE on Theory and Practice of Software Development (TAPSOFT '93), vol. 668 of Lecture Notes in Computer Science, pp. 746–761, Springer, Orsay, France, April 1993. View at Publisher · View at Google Scholar
  17. M. I. Davidich and S. Bornholdt, “Boolean network model predicts cell cycle sequence of fission yeast,” PLoS ONE, vol. 3, no. 2, Article ID e1672, 2008. View at Publisher · View at Google Scholar · View at Scopus
  18. L. Tournier and M. Chaves, “Uncovering operational interactions in genetic networks using asynchronous Boolean dynamics,” Journal of Theoretical Biology, vol. 260, no. 2, pp. 196–209, 2009. View at Publisher · View at Google Scholar · View at Scopus
  19. “Convey HC-1 Datasheet,”
  20. “DRC Dev System DS2000 Datasheet,”
  21. R. A. Arce-Nazario, E. Orozco, and D. Bollman, “A systolic array based architecture for implementing multivariate polynomial interpolation tasks,” in Proceedings of the International Conference on ReConFigurable Computing and FPGAs (ReConFig '09), pp. 77–82, December 2009. View at Publisher · View at Google Scholar