Reconfigurable Hardware Implementation of a Multivariate Polynomial Interpolation Algorithm

Abstract

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 172× and 67× 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 {0,1,2,…,𝑝−1}, 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 𝑡_𝑘−1, the reverse engineering problem seeks to find a function 𝑓∶𝐹^𝑛→𝐹^𝑛 for which 𝑓(𝑠₀)=𝑠₁, 𝑓(𝑠₁)=𝑠₂,…,𝑓(𝑠_𝑘−2)=𝑠_𝑘−1.

When 𝐹 is a finite field, function 𝑓 can be represented by a tuple of polynomials (𝑃₁,𝑃₂,…,𝑃_𝑛) where 𝑃_𝑖∶𝐹^𝑛→𝐹 for each 𝑖=1,2,…,𝑛. As pointed out in [8] such polynomials can be found by the following variation of Lagrange's interpolation formula:𝑃_𝑖(𝑥)=^𝑘−1_𝑚=0𝑓_𝑖𝑎_𝑚^𝑘−1_{𝑗=0,𝑗≠𝑚}𝑥_𝑙𝑗𝑚−𝑎_{𝑗𝑙𝑗𝑚}𝑎_{𝑚𝑙𝑗𝑚}−𝑎_{𝑗𝑙𝑗𝑚}, (1) 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 𝑃=𝑥²₁𝑥₂𝑥₃+𝑥²₁𝑥²₂+𝑥²₁𝑥₃+2𝑥²₁𝑥₂+𝑥₁𝑥₂+2𝑥₂+2𝑥²₁+2𝑥₁+1. This polynomial depends on 𝑥₁,𝑥₂, and 𝑥₃ and cannot be further simplified by any algebraic means. However, the polynomial 𝑃′=2𝑥₂+2𝑥₂𝑥₃ 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 𝑋={𝑥_𝑖1,𝑥_𝑖2,…,𝑥_𝑖𝑚}⊂𝑈 and any 𝑠∈𝑆^𝑛, we define proj_𝑋𝑠={𝑠_𝑖1,…,𝑠_𝑖𝑚}. For instance, if 𝑈={𝑥₁,𝑥₂,𝑥₃,𝑥₄,𝑥₅,𝑥₆}, 𝑋={𝑥₂,𝑥₃,𝑥₆}, and 𝑠=(1,0,3,2,0,1), the projection of 𝑠 on 𝑋,proj_𝑋𝑠, is {0,3,1}.

Let 𝐴={0,1,…,𝑎−1} and 𝐵={0,1,…,𝑏−1} where 𝑎 and 𝑏 are integers, 𝑎,𝑏≥2 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 {𝑥_𝑖1,𝑥_𝑖2,…,𝑥_𝑖𝑚} is a basis for 𝑓 if (1) for every 𝑐,𝑑∈𝐷, proj_𝑋𝑐=proj_𝑋𝑑 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 𝑐_𝑖≠𝑑_𝑖, proj_{𝑈−{𝑥𝑖}}𝑐=proj_{𝑈−{𝑥𝑖}}𝑑 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 𝑖=0,1,…,𝑏−1, 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 𝑂(2^𝑛) steps.

The following lemma, whose proof is immediate, is useful in simplifying expressions such as 𝑥₁∧𝑥₁𝑥₂=𝑥₁ and 𝑥₁∨𝑥₁𝑥₂=𝑥₁ in the computation of 𝑅.

Lemma 2. Let 𝑆={𝑥_𝑖1,𝑥_𝑖2,…,𝑥_𝑖𝑑} 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 𝑓∶{0,1,2,3}⁵→{1,2,3} be a function whose values are included in Table 2.
Applying Algorithm 1 and the reduction rules of Lemma 2, we obtain 𝑅=𝑥₂∨𝑥₅𝑥₃𝑥₁∨𝑥₅𝑥₁∨𝑥₄=𝑥₁𝑥₂𝑥₃∨𝑥₁𝑥₃𝑥₅∨𝑥₃𝑥₄𝑥₅. (2) Thus, the bases are 𝑋₁={𝑥₁,𝑥₂,𝑥₃}, 𝑋₂={𝑥₁,𝑥₃,𝑥₅}, and 𝑋₃={𝑥₃,𝑥₄,𝑥₅}.

Table 2 

for  𝑖∶=0 to 𝑏−1  do     𝑅=1;     for  𝑗∶=0 to 𝑏−1,𝑗≠𝑖  do         𝑟(𝑖,𝑗)=⋁(𝑢,𝑣)∈𝑆𝑓,𝑖×𝑆𝑓,𝑗{𝑥𝑘|𝑢𝑘≠𝑣𝑘}     end     𝑅=𝑅∧𝑟(𝑖,𝑗);     Express 𝑅 in disjunctive normal form (DNF) end

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 (proj_𝑋𝑐,𝑓(𝑐))=(proj_𝑋𝑑,𝑓(𝑑)). 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 (𝑎₀,𝑏₀), (𝑎₁,𝑏₁), …, (𝑎_𝑘−1,𝑏_𝑘−1) be a set of points in 𝐹^𝑛×𝐹, and let 𝑋 be a basis for the function 𝑓 defined by 𝑓(𝑎_𝑖)=𝑏_𝑖, 𝑖=0,1,…,𝑘−1. Then an interpolating polynomial for 𝑓 whose only variables are those of 𝑋 is 𝑃(𝑋)=^𝑘−1_{𝑚=0,𝑎′𝑚∈𝐷/≡𝑋}𝑏_{𝑚𝑗≠𝑚}𝑥_ℓ𝑗𝑚−𝑎_{𝑗ℓ𝑗𝑚}𝑎^′_{𝑚ℓ𝑗𝑚}−𝑎_{𝑗ℓ𝑗𝑚}, (3) where ℓ_𝑗𝑚 is the smallest index such that 𝑥_ℓ𝑗𝑚∈𝑋 and 𝑎^′_𝑗 and 𝑎^′_𝑚 differ.

Proof. For each 𝑎_𝑚, 𝑚=0,1,…,𝑘−1, let 𝑎^′_𝑚 be the representative of the equivalence class 𝐶_𝑚 containing 𝑎_𝑚. Then 𝑥_ℓ𝑗𝑚−𝑎_{𝑗ℓ𝑗𝑚}𝑎_{𝑚ℓ𝑗𝑚}′−𝑎_{𝑗ℓ𝑗𝑚}=0if𝑥=𝑎_𝑗1if𝑥=𝑎_𝑚 (4) and so 𝑃(𝑎_𝑚)=𝑏_𝑚, where 𝑏_𝑚 is the value of each tuple of 𝐶_𝑚.

Note that for each of the factors𝑐_𝑗𝑚(𝑥)=𝑥_ℓ𝑗𝑚−𝑎_{𝑗ℓ𝑗𝑚}𝑎^′_{𝑚ℓ𝑗𝑚}−𝑎_{𝑗ℓ𝑗𝑚} (5) in (3), 𝑐_𝑗𝑚(𝑎_𝑗)=0 and 𝑐_𝑗𝑚(𝑎_𝑚)=1 and so for all 𝑖=0,1,…,𝑘−1 and all 𝑟=1,2,…𝑐^𝑟_𝑗𝑚𝑎_𝑖=𝑐_𝑗𝑚𝑎_𝑖. (6)

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 0,1,2,3 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 𝐶₁=),(,𝐶{(0,2,1,2,30,2,1,3,1)}₂=,𝐶{(2,1,1,0,1)}₃=,𝐶{(0,1,1,2,0)}₄=),(),(,𝐶{(2,1,2,1,02,1,2,0,12,1,2,2,1)}₅=,𝐶{(0,1,2,1,1)}₆={(2,2,2,2,2)}. (7) and the set of class representatives is 𝐷≡_𝑋1=),(),(),(),(),(.{(0,2,1,2,32,1,1,0,10,1,1,2,02,1,2,1,00,1,2,1,12,2,2,2,2)} (8)
Applying (3) and making use of (6), we have 𝑃𝑋₁𝑥=1∗¹−2𝑥0−2²−1𝑥2−1+1∗¹−0𝑥2−0³−2𝑥1−2²−2𝑥1−2+2∗²−2𝑥1−2¹−2𝑥0−2³−2𝑥1−2+2∗¹−0𝑥2−0³−1𝑥2−1²−2𝑥1−2+3∗²−2𝑥1−2¹−2𝑥0−2³−1𝑥2−1+3∗¹−0𝑥2−0²−1𝑥2−1=2₁−2𝑥₂−1+3𝑥₁𝑥₂−2𝑥₃𝑥−2+4₁−2𝑥₂−2𝑥₃−2+4𝑥₁𝑥₃−1𝑥₂𝑥−2+4₁−2𝑥₂−2𝑥₃−1+4𝑥₁𝑥₂−1. (9)
Thus, a polynomial which depends only on the variables 𝑥₁, 𝑥₂, and 𝑥₃ and interpolates the tuples given by Example 2 is 𝑃𝑋₁=4𝑥₁𝑥₂+4𝑥₂𝑥₃+3𝑥₁+2𝑥₃+1. (10)

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 𝐷=⋁^𝑞−1_𝑖=0𝑐_𝑖, where 𝑐_𝑖=𝑥^𝑡𝑖0₀∧𝑥^𝑡𝑖1₁∧⋯∧𝑥^{𝑡𝑖𝑚−1}_𝑚−1, 𝑡_𝑖𝑗∈{0,1}. 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 {𝑐₀,𝑐₁,…,𝑐_𝑞−1} where each 𝑐_𝑖 is represented as a sequence ⟨𝑡_𝑖0,𝑡_𝑖1,…,𝑡_𝑖𝑚−1⟩. 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 {⟨0,1,1,0⟩,⟨0,1,0,1⟩,⟨1,0,0,0⟩}

A Boolean expression 𝐶 in CNF using only uncomplemented variables is represented as a set of disjunctions of literals (DOL) {𝑑₀,𝑑₁,…,𝑑_𝑞−1} where each DOL 𝑑_𝑖 is represented as a sequence ⟨𝑡_𝑖0,𝑡_𝑖1,…,𝑡_𝑖𝑚−1⟩. 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 {⟨1,0,1⟩,⟨1,1,0⟩,⟨0,1,0⟩}

3.2. Polynomial Representation

A multivariate polynomial 𝑃 over ℤ_𝑝={0,1,…,𝑝−1}, where 𝑝 is prime is defined as 𝑃(𝑥₀,𝑥₁,…,𝑥_𝑛−1)=∑^𝑞−1_𝑖=0𝑚_𝑖, where 𝑚_𝑖=𝛼_𝑖𝑥^𝑡𝑖0₀𝑥^𝑡𝑖1₁⋯𝑥^{𝑡𝑖𝑛−1}_𝑛−1 and 𝛼_𝑖,𝑡_𝑖𝑗∈ℤ_𝑝. For example, one multivariate polynomial over ℤ₅ is 𝑃(𝑥₀,𝑥₁,𝑥₂)=4+𝑥³₀𝑥₁+3𝑥₂+𝑥⁴₂. In our scheme, a polynomial is represented as a set of monomials {𝑚₀,𝑚₁,…,𝑚_𝑞−1} where each 𝑚_𝑖 is represented as a tuple (𝛼_𝑖,⟨𝑡_𝑖0,𝑡_𝑖1,…,𝑡_𝑖𝑛−1⟩). Binary representation requires ⌈log₂𝑝⌉ bits for each 𝛼_𝑖 and 𝑡_𝑖𝑗.

Example 6. The multivariate polynomial 𝑃(𝑥₀,𝑥₁,𝑥₂)=4+𝑥³₀𝑥₁+3𝑥₂+𝑥⁴₂ over ℤ₅ is represented as {(4,⟨0,0,0⟩),(1,⟨3,1,0⟩),(3,⟨0,0,1⟩),(1,⟨0,0,4⟩)}. For the three variable polynomial over ℤ₅ each monomial requires 4×⌈log₂5⌉=12 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.

for  𝑖∶=1 to 𝑛−1  do      for  𝑗∶=𝑖+1 to 𝑛  do    𝑥𝑖∶=𝑊(𝑥𝑖,𝑥𝑗)‖𝑥𝑗∶=𝑆(𝑥𝑖,𝑥𝑗);      end end

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𝑊𝑥_𝑖,𝑥_𝑗𝑥∶=^𝑗if𝑥_𝑖𝑥=𝜀_𝑖𝑆𝑥,otherwise,_𝑖,𝑥_𝑗∶=𝜀if𝑥_𝑖=𝑥_𝑗or𝑥_𝑖𝑥=𝜀_𝑗,otherwise. (11)

Example 7. Assume 𝑋=5,8,8,5,8,1. After computation with Algorithm 2 using (11), the result is 𝑋^′=5,8,1,𝜀,𝜀,𝜀,𝜀.

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 𝑥_𝑖1,𝑥_𝑖2,…,𝑥_𝑖𝑘,𝜀,𝜀,…,𝜀 where 𝑥_𝑖1,𝑥_𝑖2,…,𝑥_𝑖𝑘 is the subsequence of distinct elements of the input sequence.

Proof. We show by induction that after 𝑗≤𝑘 iterations of the outer loop, 𝑥_𝑖,𝑥₂,…,𝑥_𝑗,𝑥_𝑗+1,…,𝑥_𝑚=𝑥_𝑖1,𝑥_𝑖2,…,𝑥_𝑖𝑗,𝑦_𝑗+1,…,𝑦_𝑚, (12) where for each 𝑝≥𝑗+1,𝑦_𝑝∈{𝑥_𝑖𝑗+1,…,𝑥_𝑖𝑘}∪{𝜀}.
After 𝑗=1 iterations, 𝑥₁=𝑥_𝑖1. Next suppose that after 𝑗 iterations, 𝑥₁,𝑥₂,…,𝑥_𝑗,𝑥_𝑗+1,…,𝑥_𝑚=𝑥_𝑖1,…,𝑥_𝑖𝑗,𝑦_𝑗+1,…,𝑦_𝑚. (13) for some 𝑗<𝑘 where 𝑦_𝑝∈{𝑥_𝑖𝑗+1,…,𝑥_𝑖𝑘}∪{𝜀} for each 𝑝≥𝑗+1. Then either (1) 𝑦_𝑗+1=𝑥_𝑗+1 or (2) 𝑦_𝑗+1=𝜀. In the first case, 𝑦_𝑗+1=𝑥_𝑖𝑗+1 since 𝑥_𝑗+1 is not equal to any 𝑥_𝑝,𝑝≤𝑗 and is not changed by the 𝑗+1st iteration. In the second case, 𝑦_𝑗+1=𝜀 will be replaced by 𝑥_𝑖𝑗+1 during the 𝑗+1st iteration. Hence in either case, 𝑥₁,𝑥₂,…,𝑥_𝑗,𝑥_𝑗+1,…,𝑥_𝑚=𝑥_𝑖1,𝑥_𝑖2,…,𝑥_𝑖𝑗+1,𝑦_𝑗+2,…,𝑦_𝑚, (14) where for each 𝑝≥𝑗+2,𝑦_𝑝 is either 𝜀 or an element 𝑥_𝑖𝑝, where 𝑝≥𝑗+2.
Hence after 𝑖=𝑘 iterations, 𝑥₁,𝑥₂,…,𝑥_𝑚=𝑥_𝑖1,𝑥_𝑖2,…,𝑥_𝑖𝑘,𝜀,…,𝜀 (15) 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𝑊𝑥_𝑖,𝑥_𝑗𝑥∶=^𝑖+𝑥_𝑗ifmon𝑥_𝑖=mon𝑥_𝑗,or𝑥_𝑖𝑥=𝜀,_𝑖𝑆𝑥,otherwise,_𝑖,𝑥_𝑗𝑥∶=𝜀ifmon_𝑖=mon𝑥_𝑗,or𝑥_𝑖𝑥=𝜀,_𝑗,otherwise, (16) where mon(𝑥_𝑖)=⟨𝑡_𝑖0,𝑡_𝑖1,…,𝑡_𝑖𝑛−1⟩, that is, the monomial without the coefficient.

Example 8. Assume 𝑋=𝑥²,𝑥𝑦²,3𝑥²,𝑥𝑦,2𝑥𝑦. After computation with Algorithm 2 using (16), the result is 𝑋^′=4𝑥²,𝑥𝑦²,3𝑥𝑦,𝜀,𝜀.

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𝑊𝑥_𝑖,𝑥_𝑗𝑥∶=^𝑗if𝑥_𝑗⊃𝑥_𝑖,𝑥_𝑖𝑆𝑥,otherwise,_𝑖,𝑥_𝑗∶=𝜀if𝑥_𝑖⊃𝑥_𝑗,𝑥_𝑗,otherwise, (17) 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𝑊𝑥_𝑖,𝑥_𝑗∶=min𝑥_𝑖,𝑥_𝑗,𝑆𝑥_𝑖,𝑥_𝑗∶=max𝑥_𝑖,𝑥_𝑗. (18)

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 𝑖.

𝑖∶=0 ‖𝑛𝑟=1{𝑀0𝐹[𝑟]∶=𝜀0[𝑟]∶=𝜀 while stop_condition = false do     𝑖∶=𝑖+1     if  𝑖<𝑚  then  𝐹𝑖−1[0]∶=𝑥𝑖     else  𝐹𝑖−1[0]∶=𝜀     end  if     ‖𝑛𝑟=1{𝐹𝑖[𝑟]∶=𝑆(𝐹𝑖−1[𝑟−1],𝑀𝑖−1𝑀[𝑟])𝑖[𝑟]∶=𝑊(𝐹𝑖−1[𝑟−1],𝑀𝑖−1[𝑟]) end

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.

for  𝑗∶=1 to 2𝑛−3  do      parallel  for  𝑖∶= iif (𝑗≤𝑛,1,𝑗−𝑛) to |(𝑗)/2|  do    𝑥𝑖=𝑊(𝑥𝑖,𝑥𝑗−𝑖+1)‖𝑥𝑗−𝑖+1=𝑆(𝑥𝑖,𝑥𝑗−𝑖+1)      end End

Algorithm 4 

Reindexed algorithm for reduction tasks.

Figure 3 

Systolic array implementing the distinctness operation for 𝑋=0,6,2,2,7,6. 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 𝑑_sa 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 𝑑_sa=4 that implements the distinctness operation. The sequence ⟨3,5,1,3,2,6,7,6⟩ will fill the pipeline with ⟨3,5,1,2⟩ so that the final ⟨6,7,6⟩ overflows the pipeline. These 3 elements are stored in the OF for a later pass through the array. The FIFO depth 𝑑_𝑓≥𝑚−𝑑_sa where 𝑚 is the number of elements in the sequence and 𝑑_sa 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.

(a)

(b)

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 𝑑_sa 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 ⌈𝑚/𝑑_sa⌉ times.(ii)When adding polynomials (Figure 5(b)), the array will contain at most 𝑑_𝑠a 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 ⌈𝑚/𝑑_sa⌉ 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 𝑑_sa 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 0≤𝑎₀<𝑛, 𝑎₀+1≤𝑎₁≤𝑛. The distinct pairs of 𝑛-tuples produced by the CNF Generator are fed to 𝑛  ⌈log₂𝑝⌉-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 (𝑥1∨𝑥2∨𝑥3)∧(𝑥1∨𝑥4)∧(𝑥3∨𝑥5)∧(𝑥1∨𝑥2) 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 (𝑥2+1)(2𝑥3+3)(3𝑥1+2) in ℤ5. 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 (𝑑_sa) 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 𝑑_sa=256. 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 (8≤𝑛≤13) were input to implementations with diverse systolic array depths. The table shows results for the 𝑑_SA case with the shortest execution time for each 𝑛. The smaller values of 𝑛 generated fewer intermediate terms during computation, hence smaller values of 𝑑_SA 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 𝑑_SA), which still maintain competitive speedups. For example, although not shown in the table, for 𝑛=13 speedups of 30× and 45× were obtained by using 𝑑_SA 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 𝑛=8 to 12.87% for 𝑛=13.

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.

Acknowledgments

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.