Abstract:
Circuitry speeds up the QR decomposition of a matrix. The circuitry can be provided in a fixed logic device, or can be configured into a programmable integrated circuit device such as a programmable logic device. This implementation performs Gram-Schmidt orthogonalization with no dependencies between iterations. QR decomposition of a matrix can be performed by processing entire columns at once as a vector operation. Data dependencies within and between matrix columns are removed, as later functions dependent on an earlier result may be generated from partial results somewhere in the datapath, rather than from an earlier completed result. Different passes through the matrix are timed so that different computations requiring the same functional units arrive at different time slots. After the Q matrix has been calculated, the R matrix may be calculated from the Q matrix by taking its transpose and multiplying the transpose by the original input matrix.

Description:
BACKGROUND OF THE INVENTION 
     This invention relates to performing QR decomposition operations in integrated circuit devices, and particularly in programmable integrated circuit devices such as programmable logic devices (PLDs). 
     QR decomposition (also called a QR factorization) of a matrix is a decomposition of the matrix into an orthogonal matrix Q and a right triangular matrix R. QR decomposition may be used, for example, to solve the linear least squares problem. QR decomposition also is the basis for a particular eigenvalue algorithm called the QR algorithm 
     One known technique for performing QR decomposition is the Gram-Schmidt technique, which calculates the Q matrix as follows (where A is the input matrix, having columns a i  and elements a ij ): 
     for i=1 to n do
         v i =a i          

     end for
         for i=1 to n do   r ii =∥v i ∥   q i =v i /r ii      for j=i+1 to n do
           r ij =q T   i v j      v j =v j −r ij q i      
           end for       

     end for 
     As can be seen, there are two data dependencies, one in the outer loop, where q i  cannot be computed until r ii  has been computed, and one in the inner loop, where r ij  cannot be computed until q T   i  has been computed. Such data dependencies can cause delays when the computation is performed in hardware, and also may be of concern in a software implementation in a multicore processor environment, or even in a single core processor environment if the processor is deeply pipelined and the pipeline is optimized for more functions more common than division. 
     SUMMARY OF THE INVENTION 
     The present invention relates to circuitry for speeding up the QR decomposition of a matrix. The circuitry can be provided in a fixed logic device, or can be configured into a programmable integrated circuit device such as a programmable logic device (PLD). This implementation performs Gram-Schmidt orthogonalization with no dependencies between iterations. 
     In accordance with embodiments of the invention, QR decomposition of a matrix can be performed by processing entire columns at once as a vector operation. Data dependencies between matrix columns are removed, as later functions dependent on an earlier result may be generated from partial results somewhere in the datapath, rather than from an earlier completed result. There is no added hardware cost, as long as the different passes through the matrix are timed so that different computations requiring the same functional units arrive at different time slots. 
     By structuring the Gram-Schmidt algorithm differently, the result is that there are no immediate data dependencies between columns of the processed matrix. Any data dependencies that exist are separated by many clocks, allowing all calculations time to flow through a pipelined datapath with a very long latency, so that no wait states are introduced in the calculation of the Q matrix. In addition, the R matrix can be directly computed from the Q matrix. In an embodiment according to the present invention, direct access to the required Q elements allows re-use of one of the Q calculation datapaths to immediately generate an R matrix output element on each clock cycle once calculation of the R matrix begins. 
     Circuitry according to the invention implements the following modified Gram-Schmidt algorithm: 
     for i=1 to n do
         v i =a i          

     end for 
     for i=1 to n do
         r ii =∥v i ∥   for j=i+1 to n do
           r ij =v j v i /∥v i ∥   
           end for   q i =v i /r ii      for j=i+1 to n do
           v i =v j −r ijq q i      
           end for       

     end for 
     As described later, the r ij q i  term will be replaced by a set of values that already exist at some point in the datapath, completely removing q i  as a data dependency. 
     As can be seen, data dependencies have been removed in this implementation. All of the r ij  terms are computed first. There is no data dependency on ∥v i ∥ (i.e., the norm of the v i  vector, which is the square root of the sum of the squares of the elements of the vector) because the v j v i  terms are computed immediately following the v i v i  term used for the norm. The second datapath then streams all of these functions in, selectively computing either the ∥v i ∥ norm or the following r ij  terms (once viva has been computed—e.g., during computation of the norm—it preferably is latched for use in the r ij  calculations). 
     Both the input matrix A and the processing matrix Q preferably are stored in a plurality of row memories. This allows an entire column to be accessed in a single cycle by reading the corresponding elements across the row memories. 
     At the start of each iteration, an entire starting column a i  is read, and the norm of that column in calculated in two steps. The first column is latched into one of the inputs to a first datapath during the calculation of a first pass. The second input to the first datapath is the current column from the memory. First, when the current column is the first column, the inner product of the column with itself is calculated through the first datapath, and latched to the input of a second datapath. The second datapath generates the inverse norm of the first inner product, and then uses the latched first inner product to generate the subsequent values, which are the inner product of the current column with the first column, divided by the inner product of the first column with itself. The columns may be read, one per clock, to create a stream of continuous norm combinations. This sequence of operations generates the results of the first inner loop shown above: 
     r ii =∥v i ∥ 
     for j=i+1 to n do
         r ij =v j v i /∥v i ∥       

     end for 
     The same sequence of columns is read a second time, timed so that the columns arrive at the first datapath at the same time as the norm combinations for that column. The first column (in this case v i ) is latched on the input of the first datapath, which then calculates the new column iterations. Specifically, the iteration for each column is the difference between that column and the product of that column&#39;s norm combination (r ij ) multiplied by the latched first column. This sequence of operations generates the results of the second inner loop: 
     q i =v i /r ii    
     for j=i+1 to n do
         v j =v j −r ij q i          

     end for 
     Taken as a whole, the computation may be written as follows: 
     for i=1 to n do
         v i =a i          

     end for 
     for i=1 to n do
         r ii =∥v i ∥   for j=i+1 to n do
           r ij =v j v i /∥v i ∥   
           end for   q i =v i /r ii      for j=i+1 to n do
           v j =v j −r ij q i      
           end for       

     end for 
     Remembering that r ii =∥v i ∥, r ij =v j v i /∥v i ∥, and q i =v i /r ii , then r ij q i =(v j v i /∥v i ∥)×v i /r ii =(v j v i /∥v i ∥)×v i /v i ∥=(&lt;v j v i &gt;/∥v i ∥ 2 )*×v i . This means that the computation can be written as: 
     for i=1 to n do
         v i =a i          

     end for 
     for i=1 to n do
         r ii =∥v i ∥   for j=i+1 to n do
           r ii =v j v i /∥v i ∥   
           end for   q i =v i /r ii      for j=i+1 to n do
           v j =v j −(&lt;v j ,v i &gt;/∥v i ∥ 2 )×v i      
           end for       

     end for 
     The v j  vectors are read in one per clock. The (&lt;v j ,v i &gt;/&lt;v i ,v i &gt;) norm combinations are streamed in from the second datapath, calculated in the first pass. The v i  term is latched into the first datapath for the entire second pass. In other words the r ij  values do not have to be explicitly calculated and stored, as their contribution to the second pass is incorporated into the (&lt;v j ,v i &gt;/&lt;v i ,v i &gt;) norm combination. 
     As discussed in more detail below, after the Q matrix has been calculated, the R matrix may be calculated from the Q matrix by taking its transpose and multiplying the transpose by the original input A matrix, as follows:
 
R=Q T A
 
     Therefore, in accordance with the present invention, there is provided circuitry for performing QR decomposition of an input matrix. The circuitry includes a first datapath for performing multiplication and addition operations on columns of the input matrix, where the first datapath includes a plurality of multipliers, a corresponding plurality of adders each having an input connected to an output of one of the multipliers, and a summer having inputs connected to the outputs of the multipliers, and a second datapath for performing an inverse square root operation and a multiplication operation on an output of the summer of the first datapath. On a first pass, the first datapath computes respective inner products of the first indexed column of the input matrix for that iteration with each column of the input matrix, and the second datapath computes an inverse norm of the first column and multiplies that inverse norm of the first column by respective inner products of the first column with each other column of the input matrix to form respective norm combinations. On a second pass, the adders of the first datapath compute a respective difference between each other column and a product of the first column and a respective one of the norm combinations. 
     A method of configuring such circuitry on a programmable device, a programmable device so configurable, and a machine-readable data storage medium encoded with software for performing the method, are also provided. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       Further features of the invention, its nature and various advantages will be apparent upon consideration of the following detailed description, taken in conjunction with the accompanying drawings, in which like reference characters refer to like parts throughout, and in which: 
         FIG. 1  shows one embodiment, according to the invention, of a datapath arrangement for QR decomposition; 
         FIG. 2  is a cross-sectional view of a magnetic data storage medium encoded with a set of machine-executable instructions for performing the method according to the present invention; 
         FIG. 3  is a cross-sectional view of an optically readable data storage medium encoded with a set of machine executable instructions for performing the method according to the present invention; and 
         FIG. 4  is a simplified block diagram of an illustrative system employing a programmable logic device incorporating the present invention. 
     
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     A circuit implementation  100  for QR decomposition according to an embodiment of the present invention is shown in  FIG. 1 . 
     Circuit  100  includes two input matrix stores  101 ,  111 , a first datapath  102 , and a second datapath  103 . An input matrix A is input to both matrix store  101  and matrix store  111 . The calculations of orthogonal matrix Q are made on the copy in matrix store  111 . The copy in matrix store  101  is kept for the computation of right triangular matrix R after matrix Q has been computed. Each of matrix stores  101 ,  111  may include a plurality of row memories  121 . This enables entire columns to be read out in a single clock cycle by reading across the row memories—i.e., by reading corresponding elements from each row memory in the same clock cycle—which facilitates the row-based operations described above. For example, if an embodiment of the invention is implemented in programmable logic devices from Altera Corporation, of San Jose, Calif., many models of those devices provide a plurality of 9 kb user memories which may be employed for this purpose. Other memory structures also may be used. 
     Computation of matrix Q is performed in two passes. The first pass, which computes the r ij  terms, uses both datapaths  102  and  103 , while the second pass, which computes the q i  columns, uses only datapath  102 . 
     In the first pass, the inverse norm of the first column is calculated, followed by the inner product of the first and current columns divided by the inner product of the first column by itself. The norm of the first column is the square root of the inner product of the first column with itself. These operations therefore amount to taking the respective inner products of the first column with itself and each other column. The inner product of the first column with itself is used to provide the norm of the first column by taking the inverse square root of that inner product. That inner product of the first column with itself is then latched or stored at the entry to datapath  103  and used with each subsequent inner product to obtain the quotient of the inner product of each column and the norm of the first column, as described below. 
     The inner products are calculated in first datapath  102  containing multipliers  122  and adders  132 . The inverse square root is calculated at  113  in second datapath  103 . Multiplier  183  has as one input inverse square root  113 , and as its other input may have either “1” or the inverse square root  113 , so that, depending on what is required, it outputs either inverse square root  113  or it square (i.e., the inverse of the input to inverse square root  113 ). Multiplier  123  thus multiplies the input to datapath  103  either by inverse square root  113  or by the square of inverse square root  113  (i.e., divides the input by the square root or by the square of the square root). 
     The inputs to first datapath  102  during the first pass, as controlled by a bank of multiplexers  109 , are only inputs  104 ,  105  to multipliers  122 , with one of the inputs  104  always being the first column v i , while the inputs  106  to the individual adders  132  are zeroed. 
     In the first pass, multipliers  122  thus calculate the products of the elements of the first column v i  with its own elements or the elements of the current column v j . Those products  142  are then summed by summer  152  to yield the inner products v i v i  or v i v j . v i v i  is latched at  133  at the input to inverse square root  113 , which therefore always outputs 1/∥v i ∥. On the first cycle of the second pass, that output  143  is selected by output multiplexer  153  and input to first datapath  102  to calculate q i . Otherwise, output multiplexer  153  outputs a stream of norm combinations  163  (v i v j /∥v i ∥ 2 ) which are fed back to the input of first datapath  102  for use during the second pass as described below. 
     During the second pass, inputs  106  to the individual adders  132 , again selected by multiplexers  109 , are the elements of the current column v j . The other inputs  116  to adders  132  are negated (not shown) so that adders  132  subtract the outputs of multipliers  122  from the current column v j . One of the inputs to multipliers  122  is the stream of norm combinations  163 , and the other is the elements of the first column v i , latched at  162 . First, q i  is calculated by dividing each element of column v i  by the value  143  of 1/∥v i ∥. The remainder of the second pass is timed so that each column arrives at the first datapath at the same time as the norm combinations  163  for that column. The outputs of the individual adders  132  are written back to the Q memory  101 . After all j iterations have been completed, the Q matrix processing for the ith column is complete. This is repeated for all i columns. 
     The R matrix can then be calculated using R=Q T A as set forth above. In this calculation, dot products are formed from a respective row of the Q T  matrix and corresponding column of the A matrix. Because the rows of the Q matrix are stored in individual memories as described above, allowing columns to easily be read by reading across all rows at once, and because the rows of the Q T  matrix are the columns of the Q matrix, the desired dot products are those of corresponding columns of the Q and A matrices. Multiplexers  109  into the first datapath are changed again, to provide one multiplier input  104  from the A memory  101 , and the other multiplier input  105  from the Q memory  111 . Inputs  106  to the individual adders are zeroed. Each dot product will produce an element of the R matrix, which will be burst out at  173 , row by row. 
     Thus it is seen that circuitry and methods for performing QR decomposition with reduced data dependencies has been provided. This invention may have use in hard-wired implementations of QR decompositions, as well as in software implementations in multi-core processors where data dependencies across processor cores may be a factor. 
     Another potential use for the present invention may be in programmable integrated circuit devices such as programmable logic devices, where programming software can be provided to allow users to configure a programmable device to perform matrix operations. The result would be that fewer logic resources of the programmable device would be consumed. And where the programmable device is provided with a certain number of dedicated blocks for arithmetic functions (to spare the user from having to configure arithmetic functions from general-purpose logic), the number of dedicated blocks needed to be provided (which may be provided at the expense of additional general-purpose logic) can be reduced (or sufficient dedicated blocks for more operations, without further reducing the amount of general-purpose logic, can be provided). 
     The calculations described above may be done as fixed-point calculations or floating-point calculations. If floating point calculations are used, the various operators used for the calculations described above can be configured in a programmable device using, e.g., the techniques described in copending, commonly-assigned U.S. patent application Ser. No. 11/625,655, filed Jan. 22, 2007, which is hereby incorporated by reference herein in its entirety, which conserves resources by limiting the normalization of intermediate results, as described therein. 
     Instructions for carrying out a method according to this invention for programming a programmable device to perform matrix decomposition may be encoded on a machine-readable medium, to be executed by a suitable computer or similar device to implement the method of the invention for programming or configuring PLDs or other programmable devices to perform addition and subtraction operations as described above. For example, a personal computer may be equipped with an interface to which a PLD can be connected, and the personal computer can be used by a user to program the PLD using a suitable software tool, such as the QUARTUS® II software available from Altera Corporation, of San Jose, Calif. 
       FIG. 2  presents a cross section of a magnetic data storage medium  800  which can be encoded with a machine executable program that can be carried out by systems such as the aforementioned personal computer, or other computer or similar device. Medium  800  can be a floppy diskette or hard disk, or magnetic tape, having a suitable substrate  801 , which may be conventional, and a suitable coating  802 , which may be conventional, on one or both sides, containing magnetic domains (not visible) whose polarity or orientation can be altered magnetically. Except in the case where it is magnetic tape, medium  800  may also have an opening (not shown) for receiving the spindle of a disk drive or other data storage device. 
     The magnetic domains of coating  802  of medium  800  are polarized or oriented so as to encode, in manner which may be conventional, a machine-executable program, for execution by a programming system such as a personal computer or other computer or similar system, having a socket or peripheral attachment into which the PLD to be programmed may be inserted, to configure appropriate portions of the PLD, including its specialized processing blocks, if any, in accordance with the invention. 
       FIG. 3  shows a cross section of an optically-readable data storage medium  810  which also can be encoded with such a machine-executable program, which can be carried out by systems such as the aforementioned personal computer, or other computer or similar device. Medium  810  can be a conventional compact disk read-only memory (CD-ROM) or digital video disk read-only memory (DVD-ROM) or a rewriteable medium such as a CD-R, CD-RW, DVD-R, DVD-RW, DVD+R, DVD+RW, or DVD-RAM or a magneto-optical disk which is optically readable and magneto-optically rewriteable. Medium  810  preferably has a suitable substrate  811 , which may be conventional, and a suitable coating  812 , which may be conventional, usually on one or both sides of substrate  811 . 
     In the case of a CD-based or DVD-based medium, as is well known, coating  812  is reflective and is impressed with a plurality of pits  813 , arranged on one or more layers, to encode the machine-executable program. The arrangement of pits is read by reflecting laser light off the surface of coating  812 . A protective coating  814 , which preferably is substantially transparent, is provided on top of coating  812 . 
     In the case of magneto-optical disk, as is well known, coating  812  has no pits  813 , but has a plurality of magnetic domains whose polarity or orientation can be changed magnetically when heated above a certain temperature, as by a laser (not shown). The orientation of the domains can be read by measuring the polarization of laser light reflected from coating  812 . The arrangement of the domains encodes the program as described above. 
     A PLD  90  programmed according to the present invention may be used in many kinds of electronic devices. One possible use is in a data processing system  900  shown in  FIG. 4 . Data processing system  900  may include one or more of the following components: a processor  901 ; memory  902 ; I/O circuitry  903 ; and peripheral devices  904 . These components are coupled together by a system bus  905  and are populated on a circuit board  906  which is contained in an end-user system  907 . 
     System  900  can be used in a wide variety of applications, such as computer networking, data networking, instrumentation, video processing, digital signal processing, or any other application where the advantage of using programmable or reprogrammable logic is desirable. PLD  90  can be used to perform a variety of different logic functions. For example, PLD  90  can be configured as a processor or controller that works in cooperation with processor  901 . PLD  90  may also be used as an arbiter for arbitrating access to a shared resources in system  900 . In yet another example, PLD  90  can be configured as an interface between processor  901  and one of the other components in system  900 . It should be noted that system  900  is only exemplary, and that the true scope and spirit of the invention should be indicated by the following claims. 
     Various technologies can be used to implement PLDs  90  as described above and incorporating this invention. 
     It will be understood that the foregoing is only illustrative of the principles of the invention, and that various modifications can be made by those skilled in the art without departing from the scope and spirit of the invention. For example, the various elements of this invention can be provided on a PLD in any desired number and/or arrangement. One skilled in the art will appreciate that the present invention can be practiced by other than the described embodiments, which are presented for purposes of illustration and not of limitation, and the present invention is limited only by the claims that follow.