rt'\ 

I'\j 

,... 

..() VI ,-j tn 

I 

rt'\ U 0- C 

Z :::> 

a 

o -t Cl u...,... 

Zt!) 0 ..... <C Z I/) • .... w..c wl/)ao. ...... -I I/) :::> • COw~t!» 

<C u u z·. ..... ow .... c -IO::~I/):::> 

wO' ..... 1/) 

0:: J:WIU 

..J U U •• 

UJocac 

--I<ca ... 0--1 0..0 

-.IJ<lCl \t.-

o- 0:: Z. -I '. ,-j<l<lCf_ 

0-0. z"O ..... I/)~u 

'~:r- ..... O::ZJ:V') 

UUJ~ 

, ......... Zl/l 

<l U oc oc .• 

V') .... OWI/l 

<l u.. ,-' a Q) 

Zu....Jo..c 

..... UJ<l:t:~ 

~~;~::,< . .:- L~~.-· -'\ - ,.. J 

: ; .. i ," _ ~; • .: _~ ~ : 

UNIVERSITY OF CALIFORNIA 

/;1/ /! . / 

Los Angeles 

; ' ' :.=_. o.r 

Reliable and Efficient Parallel 

t' -. ,'1 r 

PrOcessIng 

Algorithms and Architectures 

for Modern Signal Processing 

-.IJ 

..() 

0-

LI'I 

rt'\ 

A dissertation submitted in partial satisfaction of the 

~ 

0 requirements for the degree 

~ 

-.IJ 

"'-

Doctor of Philosophy in Electrical Engineering 

rt'\ 

t!) 

By 

KuoJuey Ray Liu 

1990 
© Copyright by 

KuoJuey Ray Liu 

1990 
The dissertation of KuoJuey Ray Liu is approved. 

1 os rcegovac 

i 

, 

~ng~ Committee Chair 

University of California, Los Angeles 

1990 

11 
Dedication 

To my wife Ching-Ling, 

who shared this experience with love, understanding, and support. 

To my great-grand mother Lang-Shuang 

for her kindly love since my childhood. 

To my son Jeffry 

for the joyful time that we share together. 

11l 
Contents 

Dedication III 

List of Figures IX 

List of Tables X 

ACKNOWLEDGEMENTS Xl 

VITA Xll 

ABSTRACT XIV 

1 Introduction 1 

2 QRD RLS Algorithms Using Householder Transformation 8 

2.1 QRD Recursive Least-Squares Algorithm. . 9 

2.2 Systolic Block Householder Transformation. 15 

2.2.1 Vectorized SBHT QRD Systolic Array 19 

2.3 SBHT RLS Algorithm . . . . . . . . 20 

2.3.1 Vectorized SBHT RLS Array 26 

2.4 Two-level Pipelined Implementations 30 

IV 
2.5 Constrained RLS Problems. . . . . . . . . . . . . . . . . . . . . . .. 36 

3 Real-Time Algorithm-Based Fault-Tolerance for QRD RLS Systolic 

Array 

3.1 Algorithm-Based Fault-Tolerance 

3.2 Fault Model . . . . . . . . 

3.3 Real-Time Fault-Tolerance 

3.3.1 Concurrent Error Detection - Residual Method 

3.3.2 Fault Diagnosis .......... 

3.3.3 Order-Degraded Reconfiguration . 

3.4 Error Propagation and Recovery Latency 

3.4.1 Robust Error Detection 

3.4.2 Latencies 

3.5 Conservation Test . 

4 Dynamic Range, Stability, and Fault-tolerant Capability of Finiteprecision QRD RLS Systolic Algorithm 

4.1 Quasi Steady-State and Ensemble Behavior 

4.2 Dynamic Range and Lower Bound on Memory Size 

4.3 Stability and Quantization Effect . . . . . . . . 

4.4 Finite-length Effect of Fault-tolerant Capability 

4.4.1 ,Missing Error Detection 

4.4.2 False Alarm . . . . . . . 

4.4.3 Overall Memory Size Consideration 

5 Order Degraded Performance and Residual Estimations 

v 

41 

42 

46 

47 

47 

57 

62 

63 

64 

67 

69 

72 

73 

78 

86 

89 

89 

92 

96 

98 
5.1 Order Degraded Performance .. 

5.1.1 Geometric Interpretation . 

5.2 Residual Estimation in Faulty Situation 

5.2.1 Faulty Internal Cell .. 

5.2.2 Faulty Boundary Cell . 

98 

103 

104 

105 

108 

6 Multi-phase Systolic Algorithms for Spectral Decomposition 109 

6.1 Recent Developments. . . . . . . 110 

6.2 Systolic Array Matrix Processing 

6.3 QR Algorithm ......... . 

6.4 Multi-phase Systolic Algorithms 

6.4.1 Multi-phase Triangular Systolic Array 

6.4.2 Multi-phase Rectangular Systolic Array. 

6.4.3 The Hessenberg Reduction. 

6.4.4 Computing the Eigenvectors 

6.5 Performance Efficiency . . . . . . 

6.5.1 Comparisons of the arrays 

6.5.2 Rate of Convergence ... 

6.6 Efficient Fault-tolerance Schemes 

7 Conclusions and Future Research 

VI 

114 

120 

124 

125 

130 

133 

136 

138 

138 

139 

141 

146 
List of Figures 

2.1 QRD RLS systolic array .......... . 

2.2 Processing Cells of QRD RLS systolic array 

2.3 The SBHT QRD systolic array. . . . . . . . 

2.4 The processing cells of the SBHT QRD systolic array. 

2.5 The SBHT RLS systolic array obtained by direct generalization of the 

13 

14 

21 

22 

Givens rotation array.. . . . . . 27 

2.6 The SBHT RLS systolic array. . 29 

2.7 Operations of the processing cells by using modified Householder transformation. . . . . . . . . . . . . . . . . . . . . . . . . . 32 

2.8 Processing cells for two-level pipelined implementation. 33 

2.9 Operations of the processing cells . . . . 34 

2.10 The MVDR beamforming systolic array. 40 

3.11 "Matrix-matrix multiplication and LU decomposition with (weighted) 

checksum encoding . . . . . . . . . . . . . . . . . . . . . . 

3.12 Recursive LS systolic array with multiple desired responses 

3.13 Fault-tolerant recursive LS systolic array based on linear combination 

of input data in top row array feeding into column error detection array 

44 

48 

with output fault indication variable eo.. . . . . . . . . . . . . . . .. 52 

VB 
3.14 Processing cells of fault-tolerant systolic array. . . . . . . . . . . . .. 53 

3.15 Plot of leol of an adaptive QRD LS filter with order p = 3 when a fault 

occurred in P E23 from t = 25 to t = 35.. . . . . 55 

3.16 Plot of leol in Fig.3(a) with threshold set at 0.3. 55 

3.17 Plot of the hardware efficiency of the residual method versus the order 

of the LS estimation. . . . . . . . . . . . . . . . . . . . . . . . . . .. 57 

3.18 Reduced (p - 1) x (p - 1) tri-array after the deletion of the row and 

column with a faulty cell. 

3.19 Plot of expected processing latency and expected error propagation 

latency versus the order of the LS estimation. ............ . 

3.20 Comparisons of expected recovery latency for FFL and CSE methods. 

4.21 Plots of the variances in dB 

4.22 Contents of the Cells . . . . 

4.23 The first two rows of the array. 

5.24 Geometric illustration 

6.25 Triangular systolic array for QR decomposition. 

6.26 Computation of R-T x using a triarray. . . . . . 

6.27 i\lultiplication of a triangular matrix and a full dense matrix .. 

6.28 Matrix-matrix multiplication in a rectangular array. 

6.29 A circular multiplexer. . . 

6.30 A first in/first out buffer .. 

6.31 Phase 1: The QR decomposition. 

6.32 Phase 2: Computing the Q matrix. 

6.33 Phase 3: Computing the matrix product RQ .. 

V 111 

63 

68 

70 

78 

82 

93 

104 

115 

119 

120 

121 

124 

125 

127 

128 

129 
6.34 Multi-phase rectangular array for the QR iteration. 131 

6.35 System configuration of the multi-phase triarray. . . 137 

6.36 System configuration of the multi-phase rectangular array. 137 

6.37 The number of iterations for a QR algorithm to converge versus the 

matrix size . . . . . . . . . . . 140 

6.38 Row checksum for Ar = Q Rr. 143 

6.39 Column checksum for R-T A~ = Q~. 144 

IX 
List of Tables 

2.1 Comparisons of the 5BHT and Givens rotation methods 36 

3.2 The output of eo and e of an adaptive QRD L5 filter .. 54 

4.3 Mean distribution for different input data with different ,\ values. 76 

4.4 Variance distributions for different input data with different ,\ values 77 

6.5 Operations of the processing cells for different phases. . . . . . . . .. 116 

6.6 The timing table for the rotation parameters to reach the right edge 

of the QR triarray. . . . . . . . . . . . . . . . . . . . . . . . . . 117 

6.7 Comparisons of the multi-phase triarray and rectangular array. . 138 

x 
Acknow ledgments 

I would like to thank Professor Kung Yao for his counsel, comments, and constant 

encouragement on my research, especially his time and effort on reading and revising 

many papers that we published. I also would like to thank the members of my doctoral 

committee for their time spent in reviewing the dissertation research. 

I am grateful to Professor Ken Martin for his kindness to provide me the working 

environment in ICSL in which many of my research were done. I would like to thank 

S.F. Hsieh for various helpful discussions that led to many research ideas. Also, I 

This work was partially supported by the NASAl Ames under Grant 2-374, 

the National Science Foundation under Grant NCR-8814407, and the UC MICRO 

Grant. 

Finally, I cannot adequately express my deeply appreciation to my wife, ChingLing, for her constant support and love during these years of graduate study, especially 

those days that we can only buy 30e toast with an almost empty bank account. I 

also appreciate my parents for their encouragement and love. 

Xl 
1983 

1983-1985 

1986-1987 

1987 

1987 

1987 

1987-1988 

1988-1990 

1989 

1989 

VITA 

KuoJuey Ray Liu 

Born in 

B.S., Electrical Engineering Department, National Taiwan University, Taipei, Taiwan 

Communication Officer, Signal Corp. Taiwan 

Teaching Assistant, University of Michigan 

Fifth Annual VLSI Design Contest - Fifth prize: VLSI 

implementation of a median filter 

M.S.E., Department of Electrical Engineering and 

Computer Science, University of Michigan, Ann Arbor 

President Research Partnership, Office of the President, 

University of Michigan 

University Fellowship, UCLA 

Teaching Associate, UCLA 

Hortense Fishbaugh Memorial Scholarship, UCLA 

Graduate Student Award in Science and Engineering, 

Taiwanese- American Foundation 

PUBLICATIONS AND PRESENTATIONS 

1. K.J.R. Liu and K. Yao "On uniform one-chip VLSI design considerations for 

some discrete orthogonal transforms" , Proc. IEEE International Conference on 

Acoustics, Speech, and Signal Processing (ICASSP) pp.2136-2139, New York, 

April 1988. 

2. K.J .R. Liu and K. Yao, "A systematic approach to bit recursive systolic array 

design" Proc. IEEE International Conference on Systolic Array, pp.685-694, 

San Diego, May 1988. 

3. K.J.R. Liu, "Optimal Graph-Theoretical Layout of CMOS Functional Cells" 

Tech. Report, UCLA-ENG 88-29, Sep. 1988. 

Xll 
4. K.J.R. Liu and K. Yao, "Gracefully degradable real-time algorithm-based faulttolerant method for QR recursive least-squares systolic array" , in Systolic Array 

Processors, Ed. McCanny, McWhirter, and Swartzlander, ppAO 1-410, Prentice 

Hall (UK), 1989. 

5. S.F. Hsieh, KJ. R. Liu, and K. Yao, "A fast and effective algorithm for sinusoidal frequency estimation", IEEE Int'l Symposium on Information Theory, 

San Diego, Jan. 1990. 

6. K.J.R. Liu and K Yao, "Spectral decomposition via systolic triarray based on 

QR iteration" , Proc. IEEE Int'l Conf. Acoustic, Speech, and Signal Processing 

(ICASSP), pp.l017-1020, Albuquerque, April 1990. 

7. S.F. Hsieh, K.J. R. Liu, and K. Yao, "Applications of truncated QR methods 

to sinusoidal frequency estimation", Proc. IEEE Int'l Conf. Acoustic, Speech, 

and Signal Processing (ICASSP), pp.2571-2574, Albuquerque, April 1990. 

8. K.J.R. Liu, S.F. Hsieh, and K. Yao, "Recursive LS filtering using block Householder transformation" Proc. IEEE Int'l Conf. Acoustic, Speech, and Signal 

Processing (ICASSP), pp. 1631-1634, Albuquerque, April 1990. 

9. KJ.R. Liu, "Dynamic range for finite-precision QRD LS algorithm and its stability", Proc. Int'l Sym. Circuits and Systems (ISCAS), pp.3142-3145, New 

Orleans, May 1990. 

10. S.F. Hsieh, KJ.R. Liu, and K Yao, "Comparisons of truncated QR and SVD 

methods for sinusoidal frequency estimation", Proc. The second Int'l Workshop 

on SVD and Signal Processing, Kingston, June 1990. 

11. KJ.R. Liu and K. Yao, "Multi-phase systolic architectures for spectral decomposition" Proc. Int'l Conf. on Parallel Processing, August, 1990. 

12. KJ.R. Liu, S.F. Hsieh, and K. Yao, "Two-level pipelined implementation of 

systolic block Householder transformations with application to RLS algorithm" 

Proc. Int'l Conf. on Application-Specific Array Processors, Princeton, Sep. 

1990. 

13. K.J.R. Liu, S.F. Hsieh, and K. Yao, " Performance comparisons of parallel SVD 

in VLSI array processors", VLSI Signal Processing IV, San Diego, Nov. 1990. 

Xlll 
ABSTRACT OF THE DISSERTATION 

Reliable and Efficient Parallel Processing 

Algorithms and Architectures 

for Modern Signal Processing 

by 

KuoJuey Ray Liu 

Doctor of Philosophy in Electrical Engineering 

University of California, Los Angeles, 1990 

Professor Kung Yao, Chair 

Least-squares (LS) estimations and spectral decomposition algorithms constitute 

the heart of modern signal processing and communication problems. Implementations of recursive LS and spectral decomposition algorithms onto parallel processing 

architectures such as systolic array with efficient fault-tolerant schemes are the major 

concerns of this dissertation. 

There are four major results in this dissertation. First, we propose the systolic 

block Householder transformation with application to the recursive least-squares minimization. It is successfully implemented on a systolic array with a two-level pipelined 

implementation at the vector level as well as at the word level. 

XIV 
Second, a real-time algorithm-based concurrent error detection scheme based on 

the residual method is proposed for the QRD RLS systolic array. The fault diagnosis, 

order degraded reconfiguration, and performance analysis are also considered. 

Third, the dynamic range, stability, error detection capability under finite-precision 

implementation, order degraded performance, and residual estimation under faulty 

situations for the QRD RLS systolic array are studied in details. 

Finally, we propose the use of multi-phase systolic algorithms for spectral decomposition based on the QR algorithm. Two systolic architectures, one based on 

triangular array and another based on rectangular array, are presented for the multiphase operations with fault-tolerant considerations. Eigenvectors and singular vectors 

can be easily obtained by using the multi-phase operations. Performance issues are 

also considered. 

xv 
Chapter 1 

Introduction 

In order to meet the demand of high throughput and high computational complexity 

of modern signal processing, parallel processing algorithms and architectures have 

been extensively studied and implemented for dedicated applications. Rapid advances in VLSI microelectronics make it practical to build low-cost and high-density 

application-specific integrated circuits (ASIC) to meet the demands of speed and 

performance of modern signal processing. Recent VLSI/WSI technology permits the 

building of million of transistors in a single chip, while a large system may require 

hundreds of these chips to function properly. For a complex system, a single fault 

from any part of the system can make the whole system useless. For various critical 

applications, highly-reliable computations are demanded. Fault-tolerance is therefore 

needed in many of these problems. From these reasons, there is a potential proliferation of activities in the new area of fault-tolerant signal processing which explores 

reliable ways to implement fast and efficient signal processing algorithms. The goal 

of this dissertation is to study the parallel processing algorithms and architectures as 

well as associated efficient fault-tolerant techniques for modern signal processing. 

1 
Least-squares (LS) problems and spectral decomposition algorithms constitute the 

heart of modern signal processing and communications applications such as adaptive 

filtering, spectral estimation, array signal processing, channel equalization, etc .. Implementations of recursive LS and spectral decomposition algorithms onto parallel 

processing architectures such as systolic array with efficient fault-tolerant schemes 

are the major concerns of this dissertation. 

Efficient implementation of the LS algorithm, particularly the recursive LS algorithm (RLS), is needed to meet the high throughput and speed requirements of 

modern signal processing. There are many possible approaches such as fast transversal method and lattice method which can perform RLS algorithm efficiently [6, 29]. 

Unfortunately, these methods can encounter numerical difficulties due to the accumulation of round-off errors under a finite-precision implementation as summarized 

in [18]. This may lead to a divergence of the computations of the RLS algorithm [18]. 

A new type of systolic algorithm based on the QR decomposition (QRD) known as 

the QRD RLS was first proposed by l\lc\Vhirter in [72]. This algorithm is one of the 

most promising algorithms in that it is numerical stable [6, 54] as well as suitable for 

parallel processing implementation on a systolic array [29, 72]. 

Up to now, most of the QRD RLS implementations were based on the Givens 

rotation method and modified Gram-Schmidt method which both are rank-l update 

approaches [19, 26, 31, 56, 67, 72, 43]. It is well-known that the Householder transformation (HT), which is a rank-k update approach, is one of the most computation-

"J 

ally efficient methods to compute QRD. The error analysis carried out by Wilkinson 

[97, 39] showed that the HT outperforms the Givens method under finite precision 

computations. Presently, there is no known technique to implement the HT on a 

2 
systolic array parallel processing architecture, since there is a belief that non-local 

communications in the implementation are necessary due to the vector processing nature of the Householder transformation. One of the purposes of this paper is to show 

that we can implement the HT on a systolic array with only local connections. Thus, 

it is amenable to VLSI implementation and is applicable to real-time high throughput 

applications of modern signal processing. 

In Chapter 2, we first propose a systolic Householder algorithm called a systolic 

block Householder transformation (SBHT) to compute the QRD with an implementation on a vectorized systolic array. Then a RLS algorithm based on the SBHT 

called SBHT RLS algorithm is proposed to perform RLS operations on the array. 

We shall show that the SBHT array and the SBHT RLS array are generalizations 

of Gentleman-Kung's QRD array [26] and McWhirter's QRD RLS systolic array [72] 

respectively. The difficulty in the applications of the above arrays is mainly due to 

the the vectorized operations of the processing cells. This results in a high cell complexity as well as a high I/O bandwidth. By using a modified HT algorithm proposed 

by Tsao [95], a two-level pipelined implementation of the SBHT RLS algorithm can 

be achieved. That is, the algorithm is pipelined at the vector level as well as at the 

word level. The complexity of the processing cell as well as the I/O bandwidth are 

thus reduced. In general, the cell complexity of the SBHT array is higher and the system latency is longer than that of the conventional Givens rotation implementations. 

With the two-level pipelined implementation, the throughput of the SBHT RLS systolic array is as fast as that of McWhirter's Givens rotation array, and it offers better 

numerical property than the Givens method. In addition, an extension of the SBHT 

RLS array to MVDR beamformation, which is a constrainted RLS problem, is also 

3 
considered. 

Fault-tolerance has been defined as the ability of a system to execute specified algorithms correctly regardless of hardware failures and program errors [?]. In order 

to achieve the goal of fault-tolerance, redundancy has to be introduced. When we 

encounter a specific VLSI signal processing problem, an inherent nature of that signal 

processing algorithm can be used to develop a highly efficient specific fault-tolerant 

technique named algorithm-based fault-tolerance. In Chapter 3, we propose a new 

algorithm-based fault-tolerant scheme derived from the inherent nature of the QRD 

LS systolic algorithm called the residual method. For a LS problem, especially in 

communication and signal processing applications, we abstract information from the 

residuals which are the differences of the optimal LS estimations from input data and 

the desired responses. Based on the fact that a QRD LS systolic array can compute 

the residuals of different desired responses simultaneously, an artificial desired response can be designed to detect any error produced by a faulty processor. We show 

that if the artificial desired response is designed as a non-zero linear combination of 

all data inputs, the residual output of this response will be zero if no fault occurred. 

Any fault in the system will cause the residual to be non-zero and thus the fault is 

detected in real-time. Thus, the residual method can be easily incorporated with the 

systolic array antenna beamforming systems such as that considered in [77] and will 

have a great impact on the next generation of radar and sonar systems where the 

fault-tolerance scheme is quite essential for reliable real-time operation. 

Once the fault has been detected, two methods to diagnose the location of the 

faulty row are addressed. The first method, called the flushing fault location method, 

is based on the weight flushing technique to flush the weight vector out by using an 

4 
identity matrix input. The second method is to use the checksum encoding property 

to detect the row which does not meet the checksum condition. When the faulty row 

is determined, this row and the column associated with the same boundary cell are 

eliminated by a reconfiguration operation. Then the system operates in an orderdegraded manner which is acceptable in many least-squares applications. 

Though the QRD RLS algorithm is generally recognized as having good numerical 

properties such as numerical stability under finite-precision implementation [8, 60], 

there is no detailed study of this until a recent paper by Leung and Haykin [60]. 

Presently, it is still not known how to obtain the dynamic range of the algorithm in 

determining the word-length to ensure correct operation of the algorithm. 

In Chapter 4, we first observe that the cosine parameters generated by boundary 

cells will eventually reach the quasi steady-state if ,\ is close to one which is generally the case. We will show that the quasi steady-state and ensemble values of sine 

and cosine parameters are the same for all boundary cells. It is independent of the 

statistics of the input data sequence and the position of the boundary cell which generates the sine and cosine parameters. Simulation results are presented to support 

this observation. These results yield the tools needed to further investigate many 

properties of the QRD LS systolic algorithm. Then, we can obtain upper bounds on 

the dynamic range of the processing cells. Thus, lower bounds on the memory size 

can be obtained from these upper bounds on the dynamic range to prevent overflow 

and to ensure correct operations of the QRD LS algorithm. With these results, we reconsider the stability problem under quantization effects with a more general analysis 

and obtain tighter bounds than given in previous work [60]. Two important factors 

of the fault-tolerant capability, the missing error detection and the false alarm effects 

5 
are also studied in this chapter. 

When the faulty row is found in the fault-tolerant QRD LS systolic array, we enter 

an order degraded operation. In Chapter 5, we study the performance degradation 

when the order of the LS is reduced and give a geometric interpretation. For a 

short-time transient fault, we do not switch to the order degraded operation. Thus 

we propose an approximate method to estimate the optimal residual in this faulty 

situation. 

Computing the spectral decomposition of a matrix is an important issue in many 

modern signal processing and system applications. The feasibility of real-time processing for sophisticated modern signal processing systems, depends crucially on efficient 

implementation of parallel processing of the algorithms and associated architectures 

needed to perform these operations [14, 58]. While many variations exist in the literatures for solving these matrix problems, the heart of all these iterative methods 

are based either on the Jacobi-Hestennes method or the QR algorithm [30, 101, 107]. 

While there are some fundamental differences between these two approaches, both algorithms have good numerical stability and convergence rate properties and thus are 

desirable for possible implementation. Since present VLSI technology is capable of 

building a multiprocessor system on a chip, many researchers have proposed different 

parallel processing architectures to solve eigenvalue and singular value decomposition 

(SVD) problems. 

Presently, there is no known simple efficient systolic array approach for the generation of eigenvectors. The main reason is that there is no single architecture that is 

capable of handling all the steps required in the algorithm such that we can pipeline 

6 
the successive iterations readily. The communication cost among different architectures is high and the interface problem for an efficient data flow is demanding. In this 

paper, we propose two multi-phase systolic algorithms to solve the spectral decomposition problem based on the QR algorithm. By multi-phase operations we mean that 

the processing cells can perform different arithmetical operations in different phase 

of the computations. Two systolic arrays, one is triangular and the other is rectangular, are designed based on the multi-phase concept. A key feature in our method 

for the successfully application of the QR algorithm is that the Q matrix of the 

QR decomposition can be computed explicitly by multiphase operations. With the 

proper feedback of this Q matrix, the QR algorithm can be computed and pipelined 

effectively in a single systolic array. From the accumulation of those Q matrices in 

another array, eigenvectors and singular vectors can be computed without needing 

global communication inside the array. All these issues are considered in Chapter 6. 

Finally, Chapter 7 summarizes the research results of this dissertation and provides 

future research directions in these areas. 

7 
Chapter 2 

QRD RLS Algorithms Using 

Householder Transformation 

The QRD RLS algorithm is one of the most promising RLS algorithms, due to its 

robust numerical stability and suitability for VLSI implementation based on a systolic array architecture. Up to now, among many techniques to implement the QR 

decomposition, only the Givens rotation and modified Gram-Schmidt methods have 

been successfully applied to the development of the QRD RLS systolic array. It is 

well-known that Householder transformation (HT) outperforms the Givens rotation 

method under finite precision computations. Presently, there is no known technique 

to implement the HT on a systolic array architecture. In this chapter, we propose 

a Systolic Block Householder Transformation (SBHT) approach, to implement the 

HT on a systolic array as well as its application to the RLS algorithm. Since the 

data is fetched in a block manner, vector operations are in general required for the 

vectorized array. However, by using a modified HT algorithm, a two-level pipelined 

implementation can be used to pipeline the SBHT systolic array both at the vector 

8 
and word levels. The throughput rate can be as fast as that of the Givens rotation 

method. Our approach makes the HT amenable for VLSI implementation as well 

as applicable to real-time high throughput applications of modern signal processing. 

The constrained RLS problem using the SBHT RLS systolic array is also considered 

in this chapter. 

In section 2.1, a brief review of the QRD RLS algorithm is given. In section 2.2, 

the SBHT is presented while the SBHT RLS algorithm is considered in section 2.3. 

The two-level pipelined implementation of the SBHT RLS systolic array is discussed 

in section 2.4. Finally, in section 2.5, the constrained RLS problem, as applied to 

MVDR beamformation, using an extension of the SBHT RLS array is presented. 

2.1 QRD Recursive Least-Squares Algorithm 

Consider a n x p real-valued data matrix A(n) denoted by 

A(n) = [1!(1)'1!(2),·.· ,1,t(n)f = [f!(1),Q(2),··· d~(P)], (2.1 ) 

a n x 1 desired response vector u(n) = [d(1),d(2),· .. ,d(n)jT, a p x 1 weight vector 

tv( n), and a n x 1 residual vector 

f(n) = [e(1), e(2),.·. ,e(n)]T = A(n)w(n) -1L(n). 

Let the index of performance be defined by the weighted 12 norm of 

where i\(n) = diag[A(n-t),A(n-2), ... ,A,1] with a real-valued forgetting factor 0 < 

A ~ 1. Then the LS solution, satisfies 

~min(n) = min IIA(n)tv(n) - y(n)lI~ = IIA(n)w(n) - y(n)II~. (2.2) .!!L - -

9 
Assume A to be full rank with n ~ p. Then the QR Decomposition of A(n)A(n) 

yields 

Q(n)A(n)A(n) = [R(nf,O]T, 

where R(n) is a p x p upper triangular matrix and Q(n) is an n x n unitary matrix. 

Thus 

where 

with a p x 1 vector P(n) and a n x 1 vector Q(n). The LS solution w(n) can be 

obtained from 

R(n)w(n) = P(n). (2.3) 

Then ~min(n) = IIQ(n)112. 

For some radar/sonar and communication problems, the weight vector w( n) may 

not be of direct interest. For example, the residual is of interest in the multiple 

sidelobe canceller adaptive array problem. Let y.( n) be the nth new incoming row 

vector of data in A( n) of (2.1). Then the LS residual 

e(n) = !?(n)w(n) - d(n), n = p,p + 1,···, (2.4 ) 

can be obtained without computing w{n) explicitly. Thus, the complexity of the LS 

problem is further reduced and the solution of e( n) is feasible with a systolic array 

implementation [78]. 

To compute the linear LS problem recursively, a unitary matrix Q(n-1) is defined 

as 

10 
and we have 

[ 

Q(n-1) 

Q(n -1) = - - - --

!l~-l 

>.R(n - 1) 

Q(n - l)A(n)A(n) = 0 

Suppose the Givens rotation method is used for the QRD, then 1£T(n) can be annihilated by applying a sequence of Givens rotations 

(2.5) 

where the n x n transformation matrix Gi is defined by 

I i

_ 1 a a a - -- ---

a Ci a Si 

Gi = - -- ---

0 a I n - i - 1 0 

- -- - --

0 -Sj a Cj 

with 

a b 

Ci = Ja2 + b2 ' 

Sj = va2 + b2 ' 

where a and b are elements of vectors in the ith and nth rows under rotation. The 

matrix G( n) can be shown to have the form 

]((n) o h.(n) 

G(n) = o I n - p- 1 o 

h.H (n) o ,(n) 

11 
where J((n) is a p x p matrix, k(n) is a p x 1 vector, I n - p _ 1 is the (n-p-1) x (n-p-1) 

identity matrix, and ,(n) is a scalar given by ,(n) = nf=l ct(n), n ~ p with ct{n) is 

the cosine parameter associated with the ith Given rotation. The error vector can be 

transformed by Q( n - 1) to 

AR(n-1) AP(n-1) 

Q(n - l)Af(n) = o !Q{n) - AQ(n -1) 

d(n) 

vVhen the Given rotation is applied on it, we get 

[ 

R(n) 1 [ P{n} 1 G(n)Q(n -l)Af(n) = - - - !Q{n) - - - - . 

o Q(n) 

It can be easily seen that the last element of Q(n), vn(n), can be obtained naturally 

during the triangularization process [78] and is given by 

In [78], McWhirter has shown the LS residual e( n) can be expressed as 

e(n) = vn(nh(n). 

The above results lead to the systolic implementation of QRD RLS algorithm without computing weight vector explicitly [78]. The systolic array is shown in Fig. 2.1. 

It consists of two parts: a triangular array for computing QRD and a linear column 

array called response array (RA) for computing LS residual. When a new data vector 

is updated, a new triangular matrix R sits in the triangular QR array and a new 

vector P sits in the RA. 

12 
1 

-1 

a(1) __ a(2) --!(3)--

RA 

QRDAITay~ 

Figure 2.1: QRD RLS systolIc . array 

13 
(1) Boundary Cell 

iin 

.---+ (C,S) 

"Yout 

(2) Internal Cell 

~, 

(c, s) .. r .. (c, S) 

" 

Iout 

(3) Final Cell 

Xout 

If Xin == 0 then 

C +- 1; s +- o· ...... , lout +- "Yin; 

r == Ar, 

otherwise 

r' == . IA2r2 + x~ . V· In' 

C +- Ar /r'; s +- Xin/7" 

r +- r'; "'rout == C"'rin 

end 

Xout +- CXin - sAr 

r ..- SXin + cAr 

Figure 2.2: Processing Cells of QRD RLS systolic array 

14 
2.2 Systolic Block Householder Transformation 

The Givens rotation method discussed above is a rank-l update approach since each 

input is a row of data. For the systolic block Householder transformation (SBHT), 

we need a block data formulation. Denote the data matrix as 

[

if I [ X(n - 1) 1 X(n) = . = - -- : X T 

XT n n 

and the desired response vector as 

Y y(n - 1) 

[ 

Yl I 

y(n} = ;: = [ -;. - 1 

where xT is the ith data block, 

T 

X y-1}k+1 

X(i_1}k+2 

T X ik 

[X'1 X 2"' X ' ] I, 1, t,p 

E ~nkxp 

E ~nk, 

X(i-1)k+1.1 X(i-1}k+1.2 

X (i-1 )k+2.1 X(i-1)k+2.2 

XU-I)k+I,. 

X(i-1)k+2.p 

I = 

Xik.1 Xik.2 Xik.p 

and Yi is the ith desired response block, 

[ Yi = . ~~:=:~::: I 

Yik 

k is the block size and p is the order (columns) of the system. 

For a rank-k update QR decomposition, suppose we have 

15 

(2.6) 

(2.7) 

(2.8) 

E ~kxp (2.9) 

(2.10) 
[ 

R(n - 1) ] 

Q(n - 1)X(n - 1) = - ~ - . (2.11 ) 

Denote 

[ 

Q(n - 1) 

Qk(n - 1) = - - -

OT 

(2.12) 

then we have 

R(n - 1) 

o (2.13) 

If we can find a matrix H( n) such that 

[ 

R(n) ] 

H(n)Qk(n - l)X(n) = - ~ - , (2.14) 

then the new Q( n) is 

Q(n) = H(n)Qk(n - 1). (2.15 ) 

An n x n Householder transformation matrix T is of the form 

2ZZT 

T=I- W ' (2.16) 

where z E Rn. When a vector x is multiplied by T, it is reflected in the hyperplane 

defined by span{z}l.. Choosing z = x ± Ilxll2el, where el = [1,0,0,···,0] E Rn, then 

x is reflected onto el by T as 

(2.17) 

That is, all of the energy of x is reflected onto unit vector el after the transformation. 

We can zero out X;: by applying successive Householder transformations as follows, 

16 
o o 

O 0 (i-I) (i-I) •.. X· ... X , " nt" 'n,p O ... 0 0 XU). . •• y(i) , '" n.I+I, '--n.V 

for i = 1"" ,p, where x~~~ = Xn.i, R(O)(n - 1) = R(n - 1), and the resultant matrix 

H(n) is 

(2.18) 

where each H(i)(n) represents a Householder transformation which zeros out the ith 

1 f h d d X T' (i-I) co umn 0 t e up ate n' Le., Xni . 

To obtain H(1)( n), denote 

[ . T . T]T Zl = rlI - 0"1 : O(n_l)k_1 : Xn,l ' 

where rlI is the (1,1) element of R(n -1),0"; = rrt + Ilxn ,dI 2. Then from (2.16) 

hg)(n) OT h~;)T(n) 

---- ---

H(1)(n) = 0 I(n-t)k-t 0 (2.19 ) 

---- - --

h~~) (n) 0 H~~) (n) 

a k x k matrix and 

2x xT 

H(l)( ) _ I n.t n.t 

22 n - k - 2 ' 

O"ZI 

(2.20) 

with O"i l = IIztll~ = 2(0"; - O"lrll). Define tPl = O"r - O"lrlI, (2.20) can be rewritten in 

a form without multiplication of the number 2 as 

X xT 

H (I) - 1- n,1 n.1 

22 - . 

tPI 

17 
In general, 

Hi7)(n) 0 Hl';)(n) 

H(m)(n) = 0 I(n-l)k-p 0 (2.21 ) 

H~~)(n) 0 H~';)(n) 

where Hi~)(n) E ~pxp is an identity matrix except for the mth diagonal entry; 

H~';)(n) E ?Rpxk is a zero matrix except for the mth row; H~~){n) = Hi,;)T(n); and 

x(m-I)x(m-l)T 

H (m)( ) _ I n,m n,m E \Okxk 

22 n - k - tPm :n. (2.22) 

seen that H~~(n)HW(n) = O,H~?(n)H~~)(n) = 0, for Vi f. j. Thus we have the 

following lemma, 

Lemma 2.1 The Householder transformation matrix, H(n) E ?Rnkxnk, is orthogonal 

and is of the form 

Hll (11) 0 H12(n) 

H(n) = 0 I(n-l)k-p 0 = H(p)(n)H(p-l)(n) ... H(1)(n), 

H21 (n) 0 H22(n) 

(2.23) 

with 

Hll(n) - HW(n) ... H~;)(n)Hg)(n) 

H22(n) = H~i)(n) ... H~~(n)H~~)(n),O (2.24) 

18 
If the block size k = 1, due to the fact that Givens rotation method is a special case of 

rank-l update Householder transformation [30], the H matrix in Lemma 2.1 becomes 

a Givens rotation matrix G of the form [78] 

K(n) 0 h(n) 

---- ----

G(n) = 0 I n - p - 1 o ---- ----

hH(n) 0 ,(n) 

where K(n) is a p x p matrix, h(n) is a p x 1 vector, and ,(n) is a scalar given by 

,(n) = Df::::l ct(n), n 2 p where c;(n) is the cosine parameter associated with the ith 

Given rotation. 

2.2.1 Vectorized SBHT QRD Systolic Array 

Now we propose a vectorized systolic array to implement the QRD based on the 

SBHT. Similar to the QR triarray of Gentleman-Kung [29], this array has both 

boundary and internal cells. The boundary cell takes an input of block size k from the 

previous processor or directly from the input port, updates its content and generates 

reflection vector, and sends it to the right for the internal cell processing. Define 

_(i_1)T _ [ .. • (i_l)T] . _ 1 x . .. . z··· n,t - Oi_l . Tii . O(n-l)k-i . Xn,i ' -, ,p, 

and Zi = x~i,il) - aie"~ where ei is a zero vector except for a unity at the ith. position. 

When an internal cell receives the reflection vector, instead of forming the matrix 

zizT and performing matrix arithmetics, it performs an inner product operation to 

update its content Tij by doing 

(2.25) 

19 
and sends the reflected data Xn,j downward for further processing. Fig. 2.3 and 

Fig. 2.4 show the SBHT QRD array architecture and the operations of the processing 

cells respectively. When the block size is k = 1, this vectorized array degenerates to 

the Gentleman-Kung's Givens rotation triarray. 

2.3 SBHT RLS Algorithm 

The LS problem is to choose a weight vector w( n) E ~P, such that the block-forgetting 

norm of 

e2(n) el (n) I 

: =X(n)w(n)-y(n) 

en(n) 

is minimized. That is, 

where 

n 

1\f(n)IIAk = IIAk(n)f(n)11 = L A2(n-i) ·llei(n)II~, 0< A S 1, 

i=1 

Ak(n) is a block-diagonal exponential weighting matrix of the form, 

and II . 112 is the Euclidean norm, 

k 

o 

~ I E~nkxnk, ~ 

lIei(n)lI~ = L le(i-l)k+j(nW· 

)=1 

(2.26) 

(2.27) 

(2.28) 

(2.29) 

(2.30) 

The exponential forgetting weighting A is incorporated in the RLS filtering scheme to 

avoid overflow in the processors as well as to facilitate nonstationary data updating. 

20 
Figure 2.3: The SBHT QRD systolic array. 

21 

Data of 

block size k 
x 

u 

x 

k 

0 a 

u u 

~. + 1 k + 1 

y 

If IIxl1 2 = 0 then 

(1 = 0, r = Ar, 

else 

S - A2r2 + IIxm 

s -.;s 

(1 - S + SAr 

U - [ Ar:s 1 

r-s 

end 

If a = a then 

y = X, r = }.r. 

else 

t - ~-I . uT . [ :" 1 

[;j-[:"j-tu 

end 

Figure 2.4: The processing cells of the SBHT QRD systolic array. 

22 
If we know the QRD of the weighted augmented data matrix at time n (in the 

block sense, which is equivalent to the nk snapshots), which is given by 

[] [ 

R(n) 

Ak(n) X(n) = y(n) = [Qf(n)=Qf(n)] 0 

where 

[ 

Ql(n) 1 

Q(n) = - - - , 

Q2(n) 

: u(n) 1 '

: v(n) 

(2.31) 

and Ql(n) E ~pxnk and Q2(n) E ~(nk-p)xnk constitute an orthogonal transformation matrix with the former spanning the column space of the weighted data matrix 

Ak(n)X(n) and the latter the null space, and R(n) E ~pxp is an upper triangular 

matrix. The optimal weight vector can be obtained by solving 

R(n)w(n) = u(n). (2.32) 

Obviously, Adn)X(n) = Q[{n)R(n). As a result, the weighted optimal residual of 

(2.26) is, 

Q[(n)R(n)w(n) - Q{(n)u(n) - Qf{n)v(n) 

= -Qf(n)v(n), 

which lies in the null space of the weighted data matrix. 

(2.33) 

Now, suppose we have the data matrix up to time n - 1 and the QRD of Ak(nl)[X(n - 1) : y(n - 1)], then the recursive LS problem is to efficiently compute the 

optimum residual at time n from the results we have at time n - 1. In particular, we 

are interested in the new nth block of the optimal residual, 

(2.34 ) 

23 
From (2.13), (2.14) and (2.23), (2.31) can be expressed as 

H l1 {n) 0 H 12 {n) ,\R{n-1) '\u(n-1) 

[ R~n) ---- ----

u(n) ] ; 0 I(n-l)k-p 0 0 '\v(n-1) 

v{n) ---- ----

H2I (n) 0 H22(n) XT 

n Yn 

By recursion on n, we relate Q(n) and Q(n - 1) using (2.15) and have 

Hll(n) 0 H12(n) QI(n - 1) 0 

--- ---

Q(n) = 0 I(n-Ilk-p 0 Q2(n - 1) 0 

--- ---

H2I (n) a H22(n) a Ik 

H ll(n)QI(n - 1) H I2{n) 

------- ---

- Q2(n-1) 0 (2.35 ) 

------- ---

H2dn)Q1(n - 1) H22 

\Ve can see that Q2(n) is updated from QI(n - 1) and Q2(n - 1) by 

(2.36) 

Q(n)A,(n) [ Y~n;~~) 1 

AHll {n)u(n - 1) + H 12 (n)Yn 

= Av(n-1) (2.37) 

where 

(2.38) 

24 
Therefore, from (2.33), (2.36), and (2.37), the weighted optimal residual vector can 

be obtained from parameters in time n - 1 by 

Ak(n)f(n) = [~~-_1In) 1 = [ -~~~ = ~:(~ ~ 12 ~ ~r~n_-!~~l~~Vn 1 '

en(n) -Hf2{n)vn 

(2.39) 

where f( min) denotes the estimate of f at time m, m $ n, given all of the data up 

to time n. The new nth block of the optimal residual is then obtained as 

(2.40) 

For the block size of k = 1, all vector parameters in (2.40) become scalars and can 

be expressed as 

p 

en(n) = - IT CiVn, 

i=1 

(2.41) 

which was first shown by McWhirter in [78]. Note that there is little difference in the 

optimal residuals estimated by the SBHT and the Givens rotation methods. To be 

specific, the optimal residual vector in (2.40) is given by 

e(n-l)k+l((n - l)k + link) 

enk_l(nk -link) 

enk(nklnk) 

(2.42) 

while the optimal residual estimated by the Givens rotation method in (2.41) is 

(2.43) 

In this sense, the SBHT RLS gives a better estimate of the residual since it uses 

more data samples to estimate the optimal residual. Take k = 2 as an example; 

25 
the optimal residual obtained from the SBHT RLS and the Givens methods are 

[e2n_l((2n -1)12n),e2n(2nI2n)] and [e2n_l((2n -1)1(2n -1),e2n(2nI2n)] respectively. 

It is clear now that the SBHT RLS method gives a better estimate for the previous 

residual than the Givens rotation method because the former makes use of the future 

data sample at time 2n to estimate the residual at time 2n - 1, while the latter does 

not. 

2.3.1 Vectorized SBHT RLS Array 

To obtain the residual vector for RLS filtering in the systolic array, there are two 

possible approaches. The first one is to generalize the architecture of McWhirter's 

Givens rotation approach [781. A SBHT QRD array is incorporated with a RA as 

shown in Fig. 2.5. Observe that Vn in (2.38) results from the reflection computation in 

(2.37), therefore Vn is obtained naturally from the output of the RA. Each boundary 

cell then forms the matrix H~'J and propagates it down diagonal boundary cells. Since 

H~id is generated earlier than HW for i < j, equation (2.30) has to be computed from 

left to right. Obviously, we prefer to compute (2.40) from right to left such that only 

inner product computations are performed instead of matrix multiplications. As a 

result, each boundary cell performs the matrix multiplication to cumulate H~·J when 

it is propagated down diagonal boundary cells. The matrix multiplications needed in 

the boundary cells in this approach are objectionable since they not only slow down 

the throughput but also increase the complexity of the boundary cells. We note, 

McWhirter's original approach based on Givens rotation worked well since only scalars 

need to be propagated down the diagonal boundary cells and the multiplications for 

scalars can be in the reverse order. 

26 
I Skewed block input data --~I 

Matrix propagation Channel 

Desired 

response 

k 

k 

k 

k 

Residual 

vector 

Figure 2.5: The SBHT RLS systolic array obtained by direct generalization of the 

Givens rotation array. 

27 
Instead of forming the matrix H~l and propagating it down, another approach 

IS to use the facts that H~l can be expressed by using (2.22) and the reflection 

vectors are sent to the right from boundary cells as described in Section 2.2.1. From 

these observations, (2.40) can be computed in similar operations performed by the 

internal cells. A generalized architecture, as shown in Fig. 2.6, is thus introduced to 

circumvent this problem. A column array of internal cells called backward propagation 

array (BPA) is added at the right hand side to perform the backward propagation of 

v n • Each row, say the ith one, needs 2(p - i) delayed buffers as shown in Fig. 2.6. 

The Vn obtained at the output of RA is then backward propagated through the BPA. 

From (2.22), each cell of this array performs the operation 

(i-I) 

H (i)( ) I _ I _ Xn,i ( (i-:-l)T. ')' 2 1 

22 n v n - V n tPi Xn,t V n' l = p, .. " , , (2.44 ) 

where v' n is an updated V n • This is a subset of the operations performed by the 

internal cell shown in (2.25). The residual vector is obtained from the top of the 

newly added column array. 

The costs for this proposed architecture are: an increased latency time from (2p + 

1 )t3 of Mc\Vhirter's Givens method to 3ptv, where t3 represents the processing time 

for the scalar operations used in the Givens rotation method and tv is the processing 

time for vector operations used in the SBHT method; the number of delay elements 

needed increases from p to 2:i 2(p - i) = p(p - 1); and p additional internal processing 

cells. The operations of the boundary and internal cells are given in Fig. 2.5. These 

results clearly show that HT can be implemented simply on a systolic array to achieve 

massive parallel processing with vector operations. This provides an efficient method 

to obtain a high throughput rate for recursive LS filtering by using the HT method. 

28 
'"rj 

oq. 

~ 

""1 

(1) 

tV 

O'l 

1-3 

0'" 

(1) 

en 

to 

tV =:r: co 1-3 

~ 

t'"' 

en 

til 

'< 

til 

<+ 

0 

:::-: 

(') 

PI 

""1 

"1 

PI 

~ 

Skewed block input data ---"'I Desired 

response 

Residual 

vector 

~'k )"'~k ;V~k Vi.- k .... ~k /"1' ~ ; _, ... k 

~ '\/' , ~ .. ~ .. ~. 2(p-1) , 

~~-.~--~--~~. ~~ ~ ~ ~ ~ ~ ... ~ I~ !. ~ ... I. 

, " "- '" ~ '" " \,. I'll: ,7 ,7 ,7 ~fT ,7 ~ k+1 k+1 k+1 k+1 k+1 "' [7,----, 

.... h< .~~ ~rt< I.-'~k l~1' ; r ; J .. 

.. V \.Y .. V \/ 2(p-2) ..... ~k , • • I... ~ • [J]] I ... 

L ~ ~, " ~ ~~ ""I: ,7 " l~ " ~ k+ 1 v k+ 1 k+ , k+ 1 ' ~--.-'" ; L; k 

J' 

I.I'v k ~ ~ k ' ~ ... k 

')(V ~ V ' V 2(p-3) ~" ... .. [I] .. '" ,. 7 t--"-......,...... '" "\:~ --"' ~""III[' ............... k+1 k+1~ k~1 ~17--.-", 

~~ ~~k ~ , 

~ V ~ V ... ~ k ~ .. 

y .. .. 

L'" ~ ~ k+1 ' ~.." 

_ ~ ~k k+1 ----'1'-' 

.... ~ k 

'''w ~ "k ... 

Response Array '" ~ f 

" 7' ~ 

Backward Propagation 

Array 
2.4 Two-level Pipelined Implementations 

The SBHT QRD array and RLS array discussed in the above sections are derived 

from the conventional Householder transformation shown in (2.16). Due to the vector 

processing nature of the conventional method, the cells of both arrays perform vector 

operations such as inner products. This means the complexity of each cell is high 

and the I/O bandwidth is large in order to have an effective communication of vector 

data. Each cell, due to the complexity for vector processing, may lead to a large 

processor. Clearly, this is not a desirable property for VLSI implementation. Thus, 

we are motivated to find a suitable algorithm to pipeline the data down to the word 

level such that the I/O bandwidth as well as the complexity can be reduced. In 

addition, we still wish to achieve a high throughput which is needed in many modern 

signal processing applications. 

The conventional approach in computing Householder transformation, y = Tq, 

based on (2.16) is to form z and IlzW from x first and then zT q/llzl12 and q2Z(zT q/llzI12) as considered before. It can be stated in the following form: 

HT Algorithm (Conventional) 

Step 2. If Sxx = 0, then y = q. 

Step 3. If Sxx =f 0 then 

(1) S =..;5;;, z = x+ [s,O,O,,· ,,0V, 

(2) ¢ = Sxx + SXll Szq = zTq , 

(3) d = S: 1/ ¢, Y = q - dz. 

In [104], Tsao pointed out that by skipping the computation of ¢ and avoiding 

the cumbersome intermediate steps of forming vector z for further computations, a 

modified algorithm with smaller round off error and less operations can be obtained. 

30 
Only step 3 of the conventional algorithm is to be modified as follows, 

Modified HT Algorithm (Tsao, 1975) 

Step 3. If Sxx i- 0 then 

(1) S = vg;;, a = Xl + S, 

(2) Sxq = x T q , 

(3) YI = -Sxq/s, d = (ql - yd/a, Yj = qi - dXj, i = 2,···, n. 

With this algorithm, the operations of the cells of the vectorized systolic arrays can be 

modified as shown in Fig. 2.7. As we can see that, for the boundary cell, the vector u, 

which consists of the forgotten diagonal element of the upper-triangular matrix and 

one column of the input data block (updated or not), can be sent out immediately 

when the input x is available without waiting for any computations as required in 

the implementation using the conventional algorithm, Due to this advantage and the 

modified operations in the internal cells, we can then pipeline the vector operations 

down to the word level such that each cell only performs scalar operations which will 

significantly reduce the complexity of the cell. 

A two-level pipelined implementation of the modified HT algorithm is given in 

FIg, ' 'J _,0 Q an d F' Ig. ? _. 9 , The boundary cell performs three major functions - squareand-cumulate, square root, and addition. For each data block, the boundary cell 

fetches one data sample, cumulates the square of the sample, and sends the data to 

the right for internal cell processing, When all the data of the block are processed, 

the content of S is then sent down for square-root operation. The resultant S is sent 

to the right for internal cell processing as well as sent down to obtain a, which is then 

sent to the right when available. At the same time, when an internal cell receives a 

Uj, it multiplies Uj with an input Xi and cumulates all these products to obtain S. 

When S is available, it is sent down for division operation with s, which arrives at 

31 
( ( 

( ( 

x 

x 

u 

c 

(J 

y 

u 

$ 

(j 

u 

$ 

(j 

If I\Xl\2 = 0 then 

(j = 0, r = ~r, 

else 

5 _ ~2r2 + I\xll~ 

s -..;s 

(j - S + ~r u_[:] 

r-s 

end 

If (j = 0 then 

y = X. T = Ar, 

c:; _ uT. AT else [] 

~ X 

t - -5/ s 

d - (Ar - t)/(1 

r-t 

y - X - d [ ~:~ 1 

Uk+l 

end 

Figure 2.7: Operations of the processing cells by using modified Householder transformation. 

32 
I ( 

--

Figure 2.8: Processing cells for two-level pipelined implementation. 

33 
( I 

( ( 

Uk~].· ··.U:.lJ] 

.~ 

o 

I~ 

.r .. 

. 1' ~ 

.T] 

send.right Ar 

5 _ A2r2 

do i = I, k 

fetch Xi; 5 - 5 + x:; send.right Xi, 

end do 

r' _ Ar; s - VS; send.right s, 

r _ S; (7 _ r' + S; send.right (7. 

5 = O. 

dOl = 1, 1\' + 1 

fetch u, . 

if i = 1 then 

S - S + u,Ar; sendJight u,. 

else 

fetch I,_1; S = S + UtIt_l: sendJ"ight ut • 

end if 

fetch 5: r' - Ar; t -- -Sis. 

fetch 0: r -- t: d -- (r' - t)/o. 

end do 

do 1 = 1, /.: 

y, = It - dU.+l: send_down V"~ 

end do 

Figure 2.9: Operations of the processing cells 

34 
the same time, to obtain t; then t is sent down and (j again arrives at the same time 

to compute d. To compute Yi of (2.8) in Step 3, we need some registers to store Ui 

and Xi temporarily. Since data from the next block is continuously being sent into 

the system, each internal cell needs 2(k + 3) registers to store Ui and Xi as indicated 

in Fig. 2.8. When d is available, Yi is then obtained one by one and sent down for 

further processing. For data from the next block, it goes through the same processing. 

When a new d is available in the internal cell, the corresponding Xi and Ui are already 

waiting in the registers. Therefore the vector operations are successfully pipelined 

down to the word level. This means that by using the modified HT algorithm, we 

have not only pipelined the SBHT arrays at the vector level but also at the word 

level. The input data is now skewed in the word level as shown in Fig. 2.8 rather 

than in the vector level as shown in Fig. 2.6. 

Since the most time consuming operation of this two-level pipelined implementation is the square root operation which is also the critical operation in McWhirter's 

Givens rotation implementation, the throughput of this two-level pipelined implementation is as fast as that of the McWhirter's Givens array. However, with a longer 

pipeline, a longer system latency for the SBHT method is required. This is due to the 

fact that the registers of the internal cells have to be all filled before we can obtain 

the residual vector. For the SBHT RLS systolic array of order p, we have (p2 + 3p )/2 

internal cells, including the BPA. Thus, there are a total of (p2 + 3p)(k + 3) registers for the whole system. The system latency is given by ts = 2p(k + 4), which 

is linearly proportional to p and k. However, for the Givens rotation method, the 

system latency is only t" = 2p + 1. Comparisons of both RLS arrays based on the 

SBHT and Givens rotation are summarized in Table 2.1. In general, the throughput 

35 
Givens rotation SBm 

Number of cells (p2+3p)/2 (p2+5p){l 

Number of delay elements p p(p-l) 

Number of registers 0 U;+3p)(k+3) 

System latency 2p+l 2p(k+4) 

Cell complexity less higher 

Numerical stability good better 

Table 2.1: Comparisons of the SBHT and Givens rotation methods 

of the SBHT RLS systolic array is as fast as Givens rotation method. Of course, 

while the cell complexity of the SBHT array is higher, it does offer better numerical 

property [107]. A detailed backward error analysis carried out by Wilkinson showed 

that for an n x n matrix A, after n(n - 1)/2 Givens rotations, the roundoff error in 

the upper-triangular matrix is in the order of O(Kgn3/

2pIIAID [107, pp.138], while a 

series of n - 1 HT's gives about O(KhnpiIAII) [107, pp.160], with Kg and Kh being 

constants and p a floating machine computation precision constant. 

2.5 Constrained RLS Problems 

In the above sections, we have dealt with an unconstrained RLS problem. The RLS 

systolic array considered there was motivated originally by the sidelobe canceller 

beamformation problem [78]. Other practical motivation could have come from the 

36 
adaptive filtering problem [32]. However, there are other signal processing applications which are modeled by a constrained RLS problem. The MVDR beamforming 

constitutes such an example [80, 85, 97]. It is interesting to determine whether a systolic array for an unconstrained RLS problem can also be used for a constrained RLS 

problem. In [80], McWhirter and Shepherd showed an extension of the unconstrained 

RLS array to the MVDR beamforming problem. Based on their approach, we shall 

also demonstrate the implementation of a MVDR beamformation problem using a 

SBHT RLS array. 

The MVDR beamforming problem is to minimize 

(2.45 ) 

subject to the linear constraints of 

C(l)T w(l)(n) = (3«(), e = 1,···, L, (2.46) 

where L is the number of constraints. We are interested in the a posteriori residual 

vector 

(2.4 7) 

The optimal solution of the weight vector is known [80] to be given by 

• (l) _ (3(l)M- 1(n)c«() _ (3(I)R- 1 (n)a(l)(n) 

w (n) - c(l)TM-l(n)c«() - lIa(l)(n)1I2 ' (2.48) 

where M = XT(n)Ak(n)X(n) is the weighted covariance matrix, R(n) is the upper 

triangular matrix resulted from the QRD of the weighted data matrix AkX(n), and 

(2.49) 

Therefore the optimal residual vector at time n is 

(3(1) 

.(i)( ) _ . XTR- 1 ( ) «()( ) en n - lIa(l)(n)1I2 n nan. (2.50) 

37 
A crucial step needed is for the efficient recursive updating of a(t)( n). A novel approach was proposed for performing this updating [80]. Specifically, from (2.13), 

(2.14), and (2.49), 

(2.51) 

o 

where b(l)(n - 1) is an arbitrary ((n - l)k - p) x 1 vector. Then from Lemma 2.1, 

(2.13), and (2.14), we have 

o 

(2.52) 

Thus, a(e)(n) = ,\-2('\H11(n)a(f)(n - 1)) can be obtained by updating a(l)(n - 1) in 

a way similar to that u(n) is obtained by updating u(n - 1) using (2.37). The only 

differences are the input for updating a(l)(n -1) is a zero vector and a scaling factor 

,\-2. Due to the structure of H in Lemma 1, the vector b(i)(-) plays no role in the 

updating of a(l)(.). Furthermore, from (2.32) and (2.34), we have 

(2.53) 

From (2.37), we see that u(n) results from the update of [y(n - 1) : Yn], where Yn is 

the new input. Now replacing u(n) with a(l)(n) and Yn with a zero vector, we have 

(2.54) 

38 
and from (2.50), we then obtain 

(2.55) 

This equation reveals that by the proper scaling of en(n), which can be obtained 

from the SBHT RLS systolic array, we can obtain the a posteriori residual vector, 

e~t)(n), of the MVDR beamformation. Fig. 2.10 shows an extension of the SBHT 

RLS array for the new problem. Now one more data channel is needed for the RA to 

pipeline cumulation of Ila(l)(n)l/2, and the scaling of the residual vector is done at the 

bottom of the RA when a new l/a(l)(.)1/2 is available. Each RA/BPA pair in Fig. 2.10 

represents one of the ]{ constraints. The optimal a posteriori residual vector of each 

linear constraint is obtained at the output of the corresponding backward propagation 

array. 

As pointed out in [80J, there are two ways to initialize the array. One method is to 

set R(O) = 81, where 8 is a small scalar, and thus from (2.49), a«()(O) = 8-1c«(), £. = 

1"", L. Another method is to obtain R(n) to some time n, then use (2.49) to obtain 

a(t)(n). The details of a two-mode operation required for this initialization procedure 

are also considered in [80J. 

39 
ttk ... ~ aq. 

c::: 

"1 

(b 

'"--'\,. ....... J'\.. ,-, 

k+1 ~ 

..... 

0 

t-3 

t:r' 

(b 

:;: 

< 

0 

~ ~ 0'" 0 (b 

~ 

S 

0' 

"1 

2. 

t:I 

()q 

rn « 

rn 

("+ 

£. 

(") -. 

~ 

"1 

"1 

~ « 

Skcwed hlock jnrlltlbta ---_ o o 

/~ ~k /~ ~k ,"-' "k ~~k ~"k 

'oi ~ , ~ ... , V , - .. - '\. ~ '\. ~ '\. ~ '\. "'If 

"\.. Y "\.. " '\. , y-, ... , ,- -, 

k+1 k+1 "'k k+1 k+1 

h< ,~ "'k 

I ... Y ~ ... .. I.-' .-1~ k "~ k , '" ... V '>I", , ~ V ~~'VV 

.. - - '"- "" J"\.. "\.. "X 

k~-17 , ... , , ... , 

"'k 

k+1 k+1 " k+1 

.1/ ~" ~/ "k i'/ k k 

~/ ~ , ~ I , - '\. T "" ........ ,"\.. "\.. " ,7 ,---, 

k+1 k+1 k~{ 

, ","k "-'~ k ~~k 

').IV '" ... V 

- I"\" " , ' , , k+1 "k+1 , lIalll(.)II' ,/ k lIa('l)(-)II~," k 

t~" '\/ 

o ~ 

Residual 

vector 

/~~ 

2(p-1 ) ':--. k 

~ 

.. .. '\. "It: ~ _'It" 

"\.. "" , , 

k+1 " /" ~ 

2(p-2) '~ "k 

"[1]] '\,. "I[ 

'- > ~ " , ,,. k+1 J " '~ k 

.... 2(p-3) .. .. 

[l] ~ _T "\.. > "\.. " ,---, k+1 ~ i' ... ~ k 

" 

• 

""-

"- ...... 

k+1 I " "", k 

... 

, 

ilk U ~} 0: ~ , 

--

Residual 

vector 

/~~ 

... ~ ~k 

,,- 1It:" 

'-

y-, 

~,. " "'" "k .. X , ,. 

J " 'r-.. k 

.. - " ...... '-

, ... I " '~ k 

" .... 

'-~-,. 

J ~ "", k 

I'-

.. > , 
Chapter 3 

Real-Time Algorithm-Based 

Fault-Tolerance for QRD RLS 

Systolic Array 

In this chapter, we propose a new algorithm-based fault-tolerant method derived from 

the inherent nature of the QR least-squares systolic algorithm. Since the residuals 

of different desired responses can be computed simultaneously, an artificial desired 

response can be designed to detect an error produced by a faulty processor. We 

show that if the artificial desired response is designed as some proper combinations 

of the input data, the output residual of the system will be zero if there is no fault. 

However, any occurring fault in the system will cause the residual to be non-zero and 

the fault can be detected in real-time. Once the fault has been detected, the system 

enters into the fault diagnosis phase from the concurrent error detection phase. Two 

methods, the flushing fault location and the checksum encoding methods, can be 

used to diagnose the location of the faulty row. When the faulty row is determined, 

41 
this row and the associated column with the same boundary cell are eliminated by a 

reconfiguration operation. Then the system degrades in a graceful manner which is 

generally acceptable for many least-squares applications. Those eliminated processors 

enter into a self-checking phase, and when the transient fault condition is removed, a 

reconfiguration is performed to resume the normal full order operation. The analysis 

of error propagation and recovery latency is also considered in this chapter. 

This chapter is organized in the following manner. In section 3.1, we review 

some prevIOus works in algorithm-based fault-tolerance and their advantages and 

disadvantages. In section 3.2, we address the fault model. Various fault-tolerance 

schemes, including concurrent error detection, fault diagnosis, and order degraded 

reconfiguration are presented in section 3.3. In section 3.4, we prove the residual 

method is robust and discuss some important latencies. Finally, in section 3.5, we 

present some conservation tests which may also be used to detect faults. 

3.1 Algorithm-Based Fault-Tolerance 

For a given algorithm, some inherent nature of that algorithm can be used to develop a highly efficient specific fault-tolerant technique denoted as algorithm-based 

fault-tolerance. The term algorithm-based means it is an algorithm-specific and not 

a general scheme that can be applied to all general problems. An inherent property 

of many matrix computations is that checksums preserved after the computations. 

j 

This observation has led to the discovery of algorithm-based fault-tolerance which 

was first studied at University of Illinois at Urbana. A recently reported algorithmbased fault-tolerant technique, called checksum (and weighted checksum) encoding 

scheme proposed by Hwang, Abraham, Jou, Chen et al., has evolved from the study 

42 
of VLSI matrix computations [3,4,5,6]. This scheme belongs to the category of information redundancy [21]. Since few hardware and time redundancies are necessary, it 

is promising for its low-cost and low overhead for VLSI/WSI multiprocessor systems. 

Many applications of the checksum (or weighted checksum) scheme have been successfully applied to various signal processing and linear algebra operations [13]. The 

major drawback of the checksum scheme proposed in [6] is that the system throughput will be slowed down because the system clock has to be extended long enough 

to accommodate the longer signal path of the non-local interconnection caused by 

the checksum scheme. Unfortunately, local connection is one of the basic desirable 

requirements of implementations. Define the checksum vector 

eT 

= [1 1 ... 1] - ", 

and weighted checksum vector 

the column, row and full weighted checksum matrices Acw, Arw, and Ajw of a square 

n-by-n matrix A as 

A 

Acw = ~TA (3.56) 

~~A 

Arw = [ A A~ A~w ], (3.57) 

A A~ A~w 

Ajw = ~TA §:.T A§:. eT Ae - -w (3.58) 

eT A -w eT Ae -w - ~A§:.w 

43 
A x B 

ees 

A e e 

S S -

1 2 

weSl 

weS2 

R 

e -

s 

weSl 

weS2 

c 

R 

C 

S 

'---....;;::;e.....;;;;..;es~ __ 1 D 

e e 

s s 

1 2 

Figure 3.11: Matrix-matrix multiplication and LV decomposition with (weighted) 

checksum encoding 

If there is any fault occurring during the computation, the checksum criterion is 

not met and thus the fault is detected. The checksum and weighted checksum 

fault-tolerance schemes have been proposed by Huang, Jou, Chen, and Abraham 

[14, 15, 33, 41J to various signal processing and linear algebra operations such as matrix multiplication and inversion, QR decomposition, LV decomposition and singular 

value decomposition (SVD) et al .. For examples, Fig. 3.11 shows the applications 

of checksum scheme to matrix-matrix multiplication and LV decomposition. The 

weighted checksum scheme can be further used to correct errors [41, 2]. The number of error to be corrected depends on how many weighted checksum vectors are 

44 
used. The major advantages of (weighted) checksum scheme are low-cost and low 

overhead for error detection and correction of multiple processor systems. However, 

the disadvantages are: 

1. There is the problem to distinguish roundoff errors from those caused by failures 

[47J; 

2. If the intermediate results of each processor are continuously propagated to 

other processors such as QRD LS systolic array, the system throughput will 

be slowed down because the system clock has to be extended long enough to 

accomodate the longer signal path of the non-local interconnections caused by 

the checksum scheme. 

There are some difficulties in applying checksum schemes to system problems such 

as communication data equalization, and radar/sonar adaptive antenna array, where 

the signal arrives continuously and high throughput rate is required. As an example, 

for a recursive QRD LS systolic array, the data is coming in row by row. As pointed 

out in [17], the QRD of a row checksum encoding matrix Ar results in a row checksum 

encoding upper triangular matrix R,.. That is, Ar = QRr. If the checksum scheme 

is used to detect error, during the recursive operation, when a new upper triangular 

matrix R is formed, to obtain a new checksum for each row, we need to sum up each 

elements of the corresponding row. Suppose only local connection, which is one of the 

basic desirable requirements of VLSI implementations, is allowed, to prevent severe 

throughput degradation, a new channel used to pipe out the partial sum of matrix R 

has to be built. While a new content of each cell is available, this content has to be 

added to the partial checksum sent from the previous cells through the new checksum 

45 
channel and then the partial checksum has to be passed to the right for the next cell. 

With these requirements and operations, not only the complexity is increased but 

also the throughput is hindered. Even when global communications of the processors 

are allowed, the checksum scheme for fault detection still reduces the throughput of 

the system as pointed out above. 

Most recently, [4] has proposed a new algorithm-based fault-tolerant technique 

applicable to the recursive LS triangular systolic array. With the addition of one 

extra column of processor array, errors in the tri-array can be detected by observing 

the scalar output of this column array. While there are some similarities between the 

results reported in [4] and that considered in this chapter (as well as that in an earlier 

version of our paper [65]), there are much differences in assumptions, techniques, and 

results between these two approaches. In any case, these two works were performed 

independently of each other. 

3.2 Fault Model 

As the VLSI technology progresses, the geometric features become smaller. Any 

defect affecting a given part of the circuitry may cause an entire module or a logic 

block to become faulty and to produce arbitrary errors. Thus, the traditional gatelevel single stuck-at fault model is no longer appropriate for VLSljWSI system. A 

cell or module is allowed to produce arbitrary errors if any part of the cell is under 

failures [17]. However, we assume that at most one cell can be faulty at a given short 

period of time. This is based on the assumption that the system reliability is such 

that the mean time between failures is long enough and the probability of more than 

one fault occurring is very small. Some basic assumptions we need are as follows: 

46 
1. If any part of the cell become faulty, the whole cell will not function correctly; 

2. The probability of the communication links and registers failing is very small 

and thus negligible [55]. 

The second assumption is reasonable since these components are typically much simpIer and smaller than the processing cells themselves [55]. In addition, they can be 

implemented conservatively with high redundancy or with self-testing circuitry to 

mask a possible fault. 

3.3 Real-Time Fault-Tolerance 

While a demanding system is usually expected to operate under a high performance 

condition, when a fault has occurred, a slightly degraded performance is often acceptable under the circumstance. For a recursive QRD L8 systolic array operating under 

a normal fault-free condition, the optimum L8 solution is attained. A reduced order 

degraded L8 solution can still yield a reasonable performance if the row and column 

with the faulty processor can be eliminated from the computation. This concept leads 

to our proposed design of a gracefully degraded fault-tolerance approach for the QRD 

L8 systolic array. 

3.3.1 Concurrent Error Detection - Residual Method 

An inherent nature of the QRD L8 systolic algorithm is that for a given data matrix 

A, the minimization of IIAwj(n) - y .(n )\\" for many desired response vectors y., i E I, - -. "'-' 

can be performed concurrently by appending some more RA's to the systolic array 

such as shown in Fig. 3.12. The output of the ith RA, ei(n), is the minimal residual 

47 
-1 

QRD Array--./ 

e· J 

Figure 3.l2: Recursive L5 systolic array with multiple desired responses 

of l?(n)tL',(n) - d,(n), where dj(n) is the nth desired response of ith input vector 

lLj' Let '1Lo belong to the column space of A. That is, '1Lo E span{~(i), 1 ~ i ~ p}, 

then Y.n = L~=l Cj~(j). The optimal L5 residual and the associated weight vector 

for '1Lo are eo( n) = 0 and lQi( n) = [ell C2, ••• ,<;,]T, for n ~ p. The actual selection of 

{Cl' C2,' .. ,c,,} will be given later. Various extentions of these fundamental properties 

of the optimum residual of a L5 estimation problem form the basis of the proposed 

residual method approach toward concurrent error detection. Denote ~ to be the 

artificial desired response (ADR) and the associated RA as the error detection 

48 
array (EDA). 

Lemma 3.1 Given the ADR 1!.0 = g(p), the contents of the EDA are identical to 

the contents of the pth column array of the QRD triangular array, and the optimal 

residual eo, the output of the pth cell of the EDA, is always zero. 

Proof: Both arrays are p x 1 vector. The first p - 1 elements of both arrays are 

identical given by the fact that the same data are rotated by the same Ci and Si 

generated by boundary cells P Eii , 1 ::; i ::; p - 1, where P Eii denoted the processor 

at position (i, i). Thus the outputs of the (p - 1 )th cells of both arrays are identical. 

Initially, the contents of the pth cells of both array are zeros. Let the first non-zero 

output of (p - 1 )th cells be x, by the update equations of both cells, the first non-zero 

contents of both pth cells equal x. If the second non-zero output of (p - 1 )th cells 

is z, the updated content of the pth boundary cell equals y'A2x2 + z2 and that of 

the pth cell of the EDA is S1'Z + Cp,.\ = y'',.\2X2 + z2, where S1' = z/J,.\2 x 2 + z2 and 

C1' = "\x / J ,.\2X2 + z2. Therefore, the contents of both array are identical. Since the 

rotation coefficients, c1' and S1" generated by P E1'1' boundary cell are proportional to 

x and z respectively, the output of the pth cell of the EDA, eo, is Cpz - S1'x which is 

always zero.O 

If there is a fault in either the pth column of the array or the EDA, these contents 

are no longer identical and then lead to a non-zero eo. Thus a fault is detected. 

However, if there is any fault outside of these two arrays, then the errors produced 

by that fault will affect both of these arrays in the same manner (i.e., contents of 

both arrays are still identical) and resulting in a zero eo. Thus, these faults will not 

be detected by the 1!.0 = g(p) design. Clearly, we can generalize the above results by 

the following Lemma. 

49 
Lemma 3.2 Given the ADR 1!.0 = Q(k), for 1::5 k ::5 p, the contents of the kth column 

array of the QRD triarray and the first k cells of the EDA are identical. The output 

of the kth cell of the EDA is zero. The contents of celli, k + 1 ::5 I ::5 p, of the error 

detection array are all zeros. 

Corollary 3.1 A fault occurring outside of the kth column of the QRD triarray and 

the EDA will not be detected if the ADR is designed as 1!.0 = Q(k). 

Proof: From the previous discussion, a fault occurring in the ith, 1::5 i ::5 k - 1, 

column of QR array will not be detected. From Lemma 3.2, the output of the kth 

cell and the contents of cell I, k + 1 ::5 I ::5 p, of the EDA are all zeros. Thus, any 

fault occurring to the ith, k + 1 ::5 i ::5 p, column of QR array will be masked by 

these zeros. The optimal residual eo is always zero unless there is any inconsistency 

between the kth column of the QR array and the EDA.D 

From all the above observations, by selecting 1!.0 properly as given in the following 

theorem, we can detect the presence of a fault in any location of the system. 

Theorem 3.1 (Concurrent Error Detection Theorem) Consider the selection 

of the artificial desired response Y..o = 2:f=l g.(i). If there is no fault in the system, 

then the output of the EDA with Y..o as an input yields eo = O. If there is a fault in 

the system, then eo # O. 

Proof: From Lemma 3.2, each Q(i) is "zeroed out" by the ith cell of the EDA. Any 

error produced by a faulty processor, say in the ph column of the QR array, will not 

be zeroed out by the ph cells of the EDA. The output of the ph cell is then non-zero 

and propagates down to the output. Therefore, whenever el # 0, there is a fault in 

the system. 0 

50 
The ADR 'HiJ = Ef=l Q(i) is obtained by implementing a top row encoding array 

(EA) consisting of p summing cells as shown in Fig. 3.13. The response array (RA), 

with the desired response 1L as input and e as output, located at the right of the EDA 

is incorporated with the system to produce the desired residual. Once eo :/:- 0, which 

indicates the system had a fault, then e(n) is considered to be in error and will not 

be used. The error detection is thus achieved in real-time. 

Example 3.1: An adaptive filter using QRD LS systolic array with order p = 3 

is simulated. In between t = 25 and t = 35, a fault occurred in cell P E23 in such a 

way that random noise within range [-1,1] is generated. From Table 3.2, we can see, 

due to the propagation delay which will be considered in Section 3.4.2, from t = 28 to 

t = 38, eo i= 0 results from errors generated by the faulty cell. The optimal residuals 

in e from t = 28 to t = 38 are then considered faulty. After t = 39, eo then decays 

down due to the adaptive nature of the algorithm. Fig. 3.15 plots leol versus t and 

shows the adaptive effect of the algorithm. A threshold device can be used with eo to 

provide a decision on the size of the error that can be tolerated. Fig. 3.16 shows the 

leol in Fig. 3.15 with threshold set at 0.3. A generalization of the proposed scheme is 

stated below. 

Theorem 3.2 (Generalized Error Detection Theorem) Any fault occurring in 

the system can be detected if the ADR is given by'HiJ = E;=l CiQ(iL where Ci :/:- 0, 1 ~ 

i ~ p.D 

The simplest ADR that can detect fault is indeed a checksum encoded data (given 

by Theorem 3.1) which is a special case of the set of ADR given by Theorem 3.2. 

However, unlike the checksum fault detection scheme in [17, 47], Theorem 3.1 and 

51 
( ( 

EA 

_ 4) ____ Yo ---- ~3)---- a 2)----

---r= ~(l)----- RA 

t " ", ( 

EDA 

_~--.~ v '\ .,v '\ ~v '\ .v "" ~~~ ~--~ ~~---- i\...J - ~ / 1\../ 1\../ 

-1---+------~--------+_------~--~ ~ __ +-'~ ~~ __ ~ 

... 

• 

.... 

... 

, . 

., 

'. 

, 

• 

• 

-

" e 

Figure 3.13: Fault-tolerant recursive LS systolic array based on linear combination 

of input data in top row array feeding into column error detection array with output 

fault indication variable EO. 

52 
(1) Encoding Cell 

Xin 

Yin Yout 

Xout 

(2) Internal Cell 

~r 

(c,s) .. r .. (c,s) 

" 

Iout 

(3) Final Cell 

lout 

!lout 

Xout +- .z,,, 

J/out +- .z,,, + J/i" 

If(c= 1 and 8=0) then 

Xout +- Xin, r +- Ar 

else 

Xout +- CX,,, - BAr 

r +- BXin + cAr 

end if 

!lout +- '"rin!lin 

lout +- '"rin 

Figure 3.14: Processing cells of fault-tolerant systolic array. 

53 
t-l .0-0.000000 .-0.000000 

t-2 .0-0.000000 .-0.000000 

t-3 .0-0.000000 .-0.000000 

t-4 .0-0.000000 .-0.000000 

t-S .0-0.000000 .-0.000000 

t-6 .0-0.000000 .-0.000000 

t-7 .0-·0.000000 .-0.000000 

t-e .0-·0.000000 .-·0.000000 

t-9 .0-·0.000000 .-·0.000000 

t-l0 .0-·0.000000 .--0.000000 

t-ll .0-0.000000 .-·0.000000 

t-12 .0-0.000000 .-0.000684 

t-13 .0-0.000000 .-·0.1~3411 

t-14 .0-0.000000 .-·0.138241 

t-1S .0-0.000000 .-0.155885 

t-16 .0-·0.000000 .-0.065978 

t-17 .0-·0.000000 .-0.~07390 

t-18 .0-·0.000000 .-0.000090 

t-19.0-0.000000 .-·0.704896 

t-20 .0-·0.000000 .-0.011910 

t-21 .0-0.000000 .-0.053229 

t-22.0-0.000000 .-·0.787624 

t·23 .0-·0.000000 .-0.015445 

t-24 .0-0.000000 .-·0.287406 

t-25 .0-0.000000 .-·0.677166 

t-26 .0-0.000000 .-·0.050S05 

t-27 .0-0.000000 .-0.117616 

t-a .0-·0.966001 .-·0.553705 

t -29 .0--0.611989 .-0.325937 

t -30 .0-1.083665 .-0.102332 

t-31 .0--0.506217 .--0.100179 

t-32 .0--0.391151 .--0.430207 

t-33 .0-0.538918 .-0. 689~4 9 

t-l4 .0--0.426679 .--0.220855 

t-35 .0-·0.558384 .-0.312747 

t-)6 .0-0.975315 .--0.117668 

t-37 .0--0.748951 .-0.031787 

t -38 .0-0.144618 .-0.219978 

t-39 .0--0.053715 .-0.38649( 

t-40 .0--0.014142 .-0.095640 

Table 3.2: The output of eo and e of an adaptive QRD L5 filter 

54 
A PIOl of leol 

O.SI---fII----+------I------I------I 

Time 

Figure 3.15: Plot of leol of an adaptive QRD LS filter with order p = 3 when a fault 

occurred in P E23 from t = 25 to t = 35. 

A PIOl of leol with threshold sel al 0.3 1Sr------r-----..,.------r------.. 

I~Tr---~----+-----~----~ 

O.sl--+f-----I------+------I--------l 

~I O~~---~----~----~------J o 125 250 375 

Time 

.. l 

Figure 3.16: Plot of leo I in Fig.3(a} with threshold set at 0.3. 

55 
3.2 provide a real-time fault detection scheme using the inherent nature of QRD RLS 

systolic array. 

Optimal Efficiency 

Since a column of linear EDA and a row of linear EA are required, the complexity 

of this fault detection scheme is 2p. That is, 2p redundant processors is required. 

Compare to the complexity of the triangular QR LS array, (p2 + 3p) /2, it is a costeffective real-time fault detection scheme. Here we do not count the final output 

multiplier cells in the EDA and RA. 

Vve define the hardware efficiency nh to be the ratio of the hardware cost of 

implementing the algorithm to the cost of implementing the algorithm with an errordetection capability. We see 

n ( ) _ p2 + 3p 

h P - p2 + 7p' (3.59) 

Thus 1/2 ~ nh(p) ~ 1 since a single error can be detected by duplicating the hardware. \Vhen fh(p) = 1, we say the error-detection scheme is most hardware efficient. 

Define the time efficiency nt to be the ratio of the time to implement the algorithm 

and the time to implement the algorithm incorporating the error-detection scheme. 

Obviously, time efficiency is bounded by 0 ~ nt ~ 1. When nt = 1, we say the errordetection scheme is most time efficient. If an error-detection scheme is both most 

hardware and time efficiency, then it is said to be optima/. For the proposed residual 

method, clearly limp _ oo nh(p) = 1 as shown in Fig. 3.17. That is, it is asymptotically 

most hardware efficient. However, the time efficiency nt{p) = 1, 'rip, so that it is 

also most time efficient. Therefore, the residual method is an asymptotically optimal 

error-detection scheme for the recursive LS systolic array. 

56 
>. 

u 

c 

E ·GI 

LU 

I 

o. 9 

o. 8 

0.7 

0.6 

0.5 

o 

Hardware Efficiency 

/ ~ 

( 

/ 

50 100 150 200 

Order of LC;\SI-Squares (p) 

Figure 3.17: Plot of the hardware efficiency of the residual method versus the order 

of the LS estimation. 

3.3.2 Fault Diagnosis 

When a fault is detected, the system leaves the concurrent error detection phase and 

enters the fault diagnosis phase. The main purpose of this phase is to find the faulty 

processor row. Either of two methods, the flushing fault location (FFL) method or 

the checksum encoding (CSE) method, can be used to diagnose and locate the faulty 

row. The FFL method is developed under the assumption that only the residual 

output eo can be accessed externally, while the CSE method assumes that all the 

cells of the EDA can be accessed. 

57 
Flushing Fault Location Method 

During the concurrent error detection phase, a fault is detected based on unknown 

values of the incoming data in Q(k), 1 ::; k ::; p, and contents of all the cells. However, 

in the FFL method, we will control the desired incoming data as well as the contents 

of the tri-array and the EDA, in order to obtain an appropriate value in eo to locate 

the faulty processor row. In the FFL method, we do not use the operations of the 

RA cells or the EA cells. From (2.3), the weight vector W = (Wl,W2,""WI') can be 

solved by using the back substitution method of 

(3.60) 

where Pi and rij are elements of vector P and matrix R. A linear array to performed 

the back substitution as in [29] can be prevented by using the weight flushing technique 

[106]. 

In a fault-free triangular LS systolic array, by "freezing" the QRD upper triangular 

matrix R{n) and the associated column vector P(n) in (2.3), the optimum solution 

tv(n) can be "flushed" out sequentially with a skewed identity matrix input. In these 

operations, the internal cells (defined in Fig. 3.17) in a given row, act as pure bypass 

elements with no Givens rotations, when the input to the boundary cell of that row 

is zero and c = 1 and s = 0 are being propagated. For example, let p = 3 and 

(3.61 ) 

However, due to errors generated by the faulty processor, various parts of R(n) 

58 
and P stored in the array are no longer correct. A new test triangular matrix T 

and a test vector Et are load to the tri-array and EDA respectively. The values of 

T and Et can be either pre-stored or distributed by an host computer when a fault 

is detected. Specifically, define T = [il,b,'" ,ip ] as a p x p all-1's upper triangular 

matrix, where L is a p X 1 vector, and Pt as a p x 1 all-1's vector. Since Pt = i p , 

the optimal solution vector ~ = [0,0,· .. ,1]' as given by (2.3). One of the reasons 

for the selection of these all l's in T and Pt is to reduce memory requirement. Only 

one-bit register is required for each cell or distribution requirement. 

With T and Pt frozen, consider a skewed p x p identity matrix input to the first 

p columns and all zeros input to the EDA. In the absence of a fault, components of 

the optimum solution vector, ~ = [0,0,.··,1]' are outputted sequentially at eo [106]. 

Denote £' = [0, ... ,0,1,0" .. ,0] as a 1 x p vector of all zeros except for an one at the 

ith position, representing the ith row of the skewed identity matrix input to the array. 

Then the output eo( 1) in response to t~ is processed by all the cells from the first to 

the pth row of the array. In general, eo( i) is the response to t: (i.e. eo( i) = t: ilia), 

and is processed by all the cells from the ith to the pth row. As considered above, in 

the absence of a fault, ea(i) = 0, 1 :S i :S p - 1, and eo(p) = 1. However, with a fault 

in the kth row, then eo( i) =I 0, 1 ~ i ~ k, but responses due to tJ for k + 1 ~ j :S p, 

encountering only fault free processing cells from the ph to the pth row, will yield the 

correct value. This property can be used to locate the faulty row in the FFL method. 

Theorem 3.3 (Flushing Fault Location Theorem) When a fault is detected and 

the system enters the fault diagnosis phase, both the EA and RA arrays are made inoperative and all 1 's are loaded into the contents of the processor cell. A skewed 

identity p x p matrix is flushed into the system with all zeros input to the EDA. The 

59 
EDA output eo( i), 1 :::; i :::; p is obtained sequentially. Assume the first zero output at 

eo(k + 1) = 0, occurs for some k, 1:::; k :::; p - 2, then the faulty processor is in the 

kth row. If there is no such k for eo(k + 1) = 0,1 :::; k :::; p - 2, and eo(p) = 1, then 

the (p - 1 )th row is the faulty row; othenvise, with eo(p) =f. 1, the pth row is the faulty 

row.D 

Example 3.2: A QRD RLS array with order p = 5 is considered. Suppose a 

fault has occurred in P E34 • When a skewed 5 x 5 identity matrix is flushed into the 

system, due to the randomly generated noise from the faulty cell P E34 , the outputs 

from eo are given by [0.2127, -0.5714,0.7453,0.,1.]. In the absence of an error, the 

outputs should be [0.,0.,0.,0.,1.]. Since the first three elements are erroneous, based 

on Theorem 3.3, the faulty cell is in the third row. 

It is obvious the flushing of ¢ is unnecessary since computations involving the 

-1 

entire QR array are definitely incorrect in the fault diagnosis phase. 

Checksum Encoding Method 

The basic assumption of the CSE method is that all the cells of the EDA can be 

accessed. Further more, all the contents of the tri-array can be piped out in the 

diagnosis phase. In this chapter, instead of using the CSE as a fault detector as 

in [17], it is used to diagnose fault location when a fault has been detected. The 

disadvantages of the checksum scheme for real-time application is thus prevented. It 

has been shown in [17, 47] that the QRD of a row checksum matrix Ar results in 

a row checksum upper triangular matrix Rr • Let rjj(n) be the content of processor 

P Ejj of the tri-array at time nand Pj( n) be the content of the ith processor of the 

EDA at time n. 

60 
Theorem 3.4 (Checksum Encoding Theorem) Given the artificial desired response Yo = Lf=l oiQ(i), for OJ i= 0, the checksum 

p-i 

L Oiri(i+k)(n + k) = Pi(n + p - i + 1), (3.62) 

k=O 

holds for i = 1, 2, ... ,p, n = p, p + 1, p + 2, ... I if no fault has occurred. If there is 

an m such that the checksum does not hold for m ::; i ::; P, then there is a fault in the 

system and the faulty processor is in the first row that does not meet the checksum.D 

The time indices are introduced to describe the time difference for a given row 

input of data to the array. Each column of inputs, say Q.( k), is zeroed out by the kth cell 

of ED A. Thus the content of the kth cell of ED A is affected only by Q( i), k::; i ::; p. 

Therefore, if there is a faulty processor, say in row m, then all the rows below do not 

satisfy the checksum because of the error produced by the faulty one which cannot 

be zeroed out by the mth cell of EDA. 

Example 3.3: Consider a [R : Eo] matrix of intermediate results for a QRD RLS 

array with order p = 4, 

0.7930 0.7462 0.4655 0.9774 2.9821 

O. 1.2973 1.0816 1.0379 3.4168 

O. o. 0.2729 0.3643 1.7132 

O. o. o. 0.1675 0.8375 

where the (i, i + k) element, Ti,(i+k), of R, takes the value at time n + k due to time 

skewing of the input. That is, Ti,(i+k)(n + k). The ith element of P takes the value at 

time n + 5 - i. That is, Pi(n + 5 - i). As we can see, starting at the third row, the 

checksums are no longer consistent for each row. From Theorem 3.4, the faulty cell 

is in the third row. 

61 
Unlike the FFL method, the CSE method cannot stop the concurrent error detection phase immediately when a fault is detected. Because of the skewed manner of 

inputting the data, if we stop the operation immediately, the checksum property will 

not hold according to Theorem 3.4. Each processor P Eii' 1 ~ i ~ j ~ p, of the QR 

tri-array has to take j - i more data and the ith cell of the EDA has to take p + 1 - i 

more data so that the checksum is satisfied for each row of the systolic array. If each 

data requires one system clock, we observe that at most p more system clocks are 

needed to process those unfinished data after the moment a fault is detected. The 

last row takes one clock and the first row takes p clocks. Generally, the ith row takes 

p + 1 - i clocks to process the unfinished data. Thus, those rows which take fewer 

clocks can pipe their final results out to the right to check their checksum while others 

rows are still working on their unfinished data. 

3.3.3 Order-Degraded Reconfiguration 

An order-degraded performance is reasonable and often acceptable in many LS applications. A reconfiguration is needed to reroute data paths for order-degraded operation. l\1any models and approaches can be found in the literatures [58,44,57,92] for 

the reconfiguration of VLSI array processors. Here we use a similar model described 

in [58]. When the faulty row, say row k, is determined, the cells in the kth column and 

row become connection elements and enter a dormant state. In the dormant state, 

each cell tests itsdf to check its status repeatedly [58]. The reduced (p - 1) x (p - 1) 

tri-array then operates in an order-degraded LS computational manner. Fig. 3.18 

shows an example of bypassing the faulty row and the associated column to become 

an order-degraded LS array. When the (transient) fault is removed, the dormant cells 

62 
Figure 3.18: Reduced (p - 1) x (p - 1) tri-array after the deletion of the row and 

column with a faulty cell. 

reactivate and generate an interrupt immediately. The reactivation scheme recovers 

all of the cells which become connection elements before and turns them into active 

cells. Then the full-order LS operation is resumed. Detail reviews on various schemes 

and technologies of reconfiguration can be found in [13]. 

3.4 Error Propagation and Recovery Latency 

In this section, we derive and compare the performances of the proposed fault-tolerant 

schemes. 

63 
3.4.1 Robust Error Detection 

Some definitions needed in the analysis are first gIven. A processor is said to be 

faulty if a fault has occurred in it and is said to have been contaminated if it contains 

erroneous data or states. The faulty moment is the time instant a fault first occurs 

within the processor, and the contaminating moment is the time instant the first error 

occurs in the processor [98]. A fault is said to be observable if some of the observable 

states, such as signal data and system output, of the processor can be affected by 

the fault. It may take some latent period to change the observable states incorrectly 

once the fault occurs. An immediately observable fault is a fault which affects some 

of the observable states right after its occurrence. To make the analysis tractable, 

we first assume all faults are immediately observable faults such that when the fault 

occurs, it produces error at the output of the faulty processor immediately. Later, 

we will relax this assumption for more general results. The error is then propagated 

and contaminates those processors affected. The moment the error is first observed 

at the output of the system is called the error observed moment. First, assume a 

fault occurs in an internal cell P Eij , i i- j, at a faulty moment. The output of this 

faulty cell is thus erroneous and can be described by X~ut = Xout + 0, where Xout is 

the fault-free output and 0 is the error generated by the fault. The error propagation 

path can be described by 

P Eij --+ P E(i+I)j --+ ... --+ P E jj , 

and then P Ek1 , k 2 j, I 2 j are all contaminated. 

From the operations executed by the internal cell, the error is modified to Ci+l {y by 

P E(i+I)j and the cumulative modifications of the error before reaching the boundary 

64 
cell, P Ejj, is 

j-1 

11 = b II Ck, 

k=i+1 

(3.63) 

where Cj is the cosine parameter generated by the boundary cell P Eii . Therefore, Cj 

and S j are erroneous and are given by 

>.r 

J>.2 r 2 + (Xin + 11)2' 

(Xin + 11) (3.64) 

In this case, Sj is no longer proportional to Xin, Q(j) will not be zeroed out by the ph 

cell of the EDA. The size of the error generated by this cell is 

>'r11 71· = C·X· - s·>.r = --r=::;r=~===::;;: '1) J m J v' ,2 + 2 . + 2' r 1Jxm 11 

(3.65) 

where r' = J >.2r2 + x;n is the new updated uncontaminated value of the content of 

P Ejj . Although those Ck'S and Sk'S for j ~ k ~ p are contaminated, Q( k), j ~ k ~ p, 

are zeroed out by the kth cell of the EDA because of the consistency of the sine 

and cosine parameters used by the kth column array of the QR tri-array and the 

EDA. \\Then 1Jj propagates down to the output of the EDA, 1Jj is influenced by the 

contaminated cosines c' of each following row. The error output at eo due to an error 

[) generated at P Eij is then given by 

5 (. .) eo 1,) = 

i ~ j - 2, 

(3.66) 

i=j-l. 

where I = 11;::-11 Cj n~=J c~. Next, assume a fault occurs in a boundary cell, P Ejj , 1 ~ 

j ~ p, at the faulty moment. Both erroneous cj and sj produced by P Ej; can be 

65 
written by 

(3.67) 

where be and b6 represent errors in the numerators while r~ represents the erroneous 

content of the denominators of Cj and Sj. The error produced by the ph cell of the 

EDA is then given by 

(3.68) 

and the output error at eo due to a faulty boundary cell is given by 

(3.69) 

A fault can occur either in the internal or boundary cell. First assume the fault 

occurs in an internal cell. The cosine parameter, Cj, produced by P Ejj is bounded 

by 0 :::; Cj :::; 1. 

Lemma 3.3 For some CJl if there exists an n E I such that cj(n) =F 0, then cj(m) > 0 

for all m > n, mEl. 

Proof: The cosine parameter IS gIven by CHI = Ar(k)/r'(k), where r(k + 1) = 

r'(k) = JA2 r 2 (k) + xtn(k). 3n ---+ cj(n) =F 0 is equivalent to say that 3n ---+ r(n) =F 0 

or r(n) > O. Since r(k+1) = r'(k) ~ Ar(k), we have r(k) > 0 for all k ~ n. Therefore, 

cj(k) > 0 for all k > n.D 

For linearly independent input column vectors, all cj(n) =F 0, and from (3.66), we 

see eg =F 0 if an e~ror has been introduced. Thus the fault can always be detected in 

this case. When the fault occurs in a boundary cell, from (3.69) we conclude that the 

only chance of a loss of detection is that eg = O. This happens only if bel b6 = Ar I Xin, 

which is very unlikely since the value of 8e. 86 , r, and Xin are not related in any 

66 
manner. Furthermore, even if he/h~ equals Ar(n)/Xin(n) for some instant of time at 

n, it is even more unlikely that he/h, equal Ar(n + 1)/xin(n + 1) at time n + 1. We 

can summarize these results in the following theorem. Note that the statement of 

Theorem 3.5 is valid only for infinite-precision implementation. 

Theorem 3.5 (Robust Error Detection Theorem) An error produced by a faulty 

processor at the faulty moment will be detected at the EDA output eo and the probability of error detection given a fault occurs equal one. That is, 

Pr( error detected at eo I a fault occurred) = 1. 0 (3.70) 

3.4.2 Latencies 

Now, we consider some basic issues related to latencies in the array. 

Definition 5.1: The system latency, t~, is the time between the moment of data 

input to the system and the moment of the output of this data from the system. 

Definition 5.2: The processing latency, tp , of processor P Eij is the time between 

the moment a data in a wavefront inputs to the system and the moment P Eij is 

processing data from that wavefront. 

Definition 5.3: The error propagation latency, t~, is the time between the faulty 

moment and the error observed moment. It is clear that the system latency of the 

QR recursive LS array depends on the number of processors and delay elements on 

the boundary and is given by ts = 2p + 1. The processing latency of processor 

P Eij , 1 ~ i ~ j ~ p + 1, is given by tp = (i + j) - 1. Since there are a totally of 

p(p + 3)/2 processors, the expected processing latency is 

(3.71 ) 

67 
Eltpected Processing and Error Propagation Latency 

1I0,..------r-------r------,..----_ 

Order or Lc~~I·Sql1~rc< (p) 

Figure 3.19: Plot of expected processing latency and expected error propagation 

latency versus the order of the LS estimation. 

The error propagation latency is given by 

(3.72) 

and the expected value is 

E(t ) = (p+ l)(p+ 2) ~ (p + 3) . (3.73) 

Fig. 3.19 shows the plots of E(t p ) and E(tt}. 

Definition 5.4: The fault diagnosis time, t j, of a faulty processor P Eij is the minimum time required to locate the faulty row right after the error observed moment. 

Definition 5.5: The recovery latency, tr , be the time between the faulty moment 

and the moment the faulty row is determined. Since the system latency is 2p + 1, 

for the FFL method the fault diagnosis time of processor P Eij for the array is 

68 
tfFL = (2p + 1) + i. We can show that the fault diagnosis time for the CSE method is 

t<jSE = p+2-i. The expected value for fault diagnosis time are E(tfFL) = (5/2)p+1 

and E(t<jSE) = (1/2)p + 2 respectively. By the definition of the recovery latency, 

we have tr = t~ + tf. Therefore, the recovery latency are t;FL = 4p - j + 3 and 

t?SE = 3p - 2i - j + 4, while the expected recovery latency are 

7p2 + 23p + 10 

2p+6 

3p2 + 13p + 16 

2p + 6 (3.74) 

respectively. Due to the facts that multiple ports can be accessed externally and we 

can use the parallel pipe-out feature of the CSE method, it is not surprised that the 

performance of the CSE method is better than that of the FFL method as indicated 

in Fig. 3.20. However, for both cases, the order of the expected recovery latency is 

O(p), which is linear with respect to p. In practice, a (transcient) fault may not 

be necessarily an immediately observable fault. Without this assumption, all the 

values obtained in this section become the lower bounds of those parameters. That 

is, performance obtained by the assumption of immediately observable fault is the 

best performance we can achieve. 

3.5 Conservation Test 

Thus far the fault-tolerance system discussed above is designed to detect the fault 

occurring in the QR triarray and the EDA. It is not applicable to the RA which 

computes the real desired response. In this section, we introduce some conservation 

properties to tackle this problem and the data is still assumed to be real. 

A: Cell Level Test 

69 
UJ: VI: 

u: 

~I 

Expected Recovery Latency of FFL & CSE 

352~--------~----------~--------~--------~ 

2~r---------~----------~--------~~--------4 

176.~---------+--------~~---------+--------~ 

................ ~~ .... ............... 88~--------~~--------~.-.. ~ .. ~~----+----------4 

......................... 

2S 7S 100 

Order pf Lc;'~t.Sq"arc.~ (p) 

Figure 3.20: Comparisons of expected recovery latency for FFL and CSE methods. 

AI: Energy Conservation The rotation operation of internal cell is described in 

Fig.I. Denote the updated r as r'. Then it can be seen that 

2 12 2 \ 2 2 

Iout + r = Iin + A r , (3.75) 

which means the energy, (i.e. the 2-norm) is conserved before and after the 

operation. This is due to the fact that the Givens rotation transformation is an 

unitary transformation which conserves the energy. 

A2: Inverse of TJnitary Matrix The above energy conservation test requires the 

square operation which is different from the operation of the internal cell. In a 

VLSIjWSI system, it is desirable to keep the number of different kinds of processors or cells as low as possible [74, 54]. Observe that the inverse of a unitary 

70 
matrix is the Hermitian transpose of it. If there is no fault occurring during the 

computations, the inverse computation should give the original values. That is 

[ 

Xin ]_ [c s] [ Xout ]. 

Ar -s c r' 

(3.76) 

One of the above inverse computation is sufficient to detect a fault. For example, 

the inverse computation of Xin = CXout + sr' will not match the original Xin, if 

Xout or r' is erroneous. This kind of test can be carried out by the internal cell. 

B: System Level Test 

Consider the RA which computes the real desired response. Based on the previous 

discussion, the energy is conserved for unitary transformation. Denote Pi ( n), 1 ~ 

i ~ p, to be the content of the ith cell of the array at time nand d as defined in 

section 2. The system level energy conservation is described by 

p n 

L P/(n + i-I) = L d2(j). (3.77) 

i=1 j=1 

That is, the energy of the input signal is equal to the energy in the RA. Note that the 

introduced time index is used to describe the operations of continuously updating the 

incoming signal. Both cell and system level tests can be built either in the system or 

applied externally depending on the requirements of the application. 

71 
Chapter 4 

Dynamic Range, Stability, and 

Fault-tolerant Capability of 

Finite-precision QRD RLS 

Systolic Algorithm 

In this chapter, we first observe that the cosine parameters generated by boundary 

cells will eventually reach quasi steady-state if ,,\ is close to one which is the usual 

case. We will show that the quasi steady-state and ensemble values of sine and cosine 

parameters are the same for all boundary cells. It is independent of the statistics of 

the input data sequence and the position of the boundary cell which generates the sine 

and cosine parameters. Simulation results are presented to support this observation. 

These results yield the tools needed to further investigate many properties of the 

QRD LS systolic algorithm. Then, we can obtain upper bounds of the dynamic 

range of processing cells. Thus, lower bounds on the memory size can be obtained 

72 
from upper bounds of the dynamic range to prevent overflow and to ensure correct 

operations of the QRD LS algorithm. With these results, we reconsider the stability 

problem under quantization effects with a more general analysis and obtain tighter 

bounds than given in previous work [60]. Two important factors of the fault-tolerant 

capability, the missing error detection and the false alarm effects are also studied in 

this chapter. 

Quasi steady-state of the rotation parameters is discussed in section 4.1. Dynamic 

range and lower bound on memory size are derived in section 4.2. Stability and 

quantization effects are studied in section 4.3. Finally, the fault-tolerant capability is 

presented in section 4.4. 

4.1 Quasi Steady-State and Ensemble Behavior 

From the updated recursive equation of the boundary cell (Fig. 2.2), we have 

k 

r2(k + 1) = -\2r2(k) + x2(k) = I: -\2i x2(k - i). ( 4.78) 

i=O 

Assume {:r} is a bounded random sequence of zero-mean and variance (72, then the 

expected value of r2(k + 1) is 

( 4.79) 

For 1,\1 < 1, then 

( 4.80) 

Since J is a concave function, from Jensen inequality 

(4.81) 

73 
and from (4.78) 

IXminl < limr(k+l)< IXmarl 

v'1 - A2 - k-oo - VI - A2 (4.82) 

where IXmarl and IXminl are the maximum and minimum values of {Ix!}, respectively. 

The cosine parameter of the Givens rotation is computed by c(k+l) = Ar{k)jr(k+ 

1). The steady-state of this parameter exists if limk_oo c( k) exists. For the sequence 

{c(·)} to have a steady-state, we need limk_oo r(k)jr(k+l) = a, where a is a constant. 

If a < 1, then the sequence {r(.)} is unbounded which conflicts with (4.82) that 

indicates {r(.)} should be bounded; if a > 1, then limk_oo r( k) = 0 which, again, 

conflicts with (4.82) (unless {x} is a zero sequence). Therefore, a has to be a unity 

to guarantee the steady-state of {c(·)} exists. That is, 

1

. r( k) 1m -+ 1, 

k-oo r(k + 1) ( 4.83) 

and the steady-state value of cosine, if exists, is 

1

· (k ) l' Ar(k) \ 1m c . + 1 = 1m (k ) = II. k ..... oo k-oo r . + 1 (4.84 ) 

From (4.78), we can see that if A = 1, then limk ..... oo r(k) -+ 00 such that limk_oo r(k)jr(k+ 

1) = 1. In this case, though the steady-state of {c(.)} exists, {r(.)} is unbounded. 

Usually A is chosen between .99 and 1 which is very close to onel . When we update 

r(k) to r(k + 1) using (4.78), a A portion of r(k) is forgotten and an input x(k) is 

added into it. If A is close to one, when k is very large, r( k) will come close to r( k + 1) 

and the input x( k) plays less and less significant role in computing r( k + 1) as the 

) 

case when A = 1. It is obvious that for A close to one, 

lim Er(k) ~ lim Er(k + 1). k ..... oo k-oo 

1 For different expressions as in [32, 60, 78], A is between .98 and 1 

74 
Therefore, from the averaging principle [76] which has been used successfully in many 

situations, the expected cosine can be approximated by 

. Er(k) hm Ec(k+ 1):::). E (k ) ~).. k-oo r + 1 ( 4.85) 

When>. is close to one, from above discussions, we have 

lim c(k) = lim ~(k)) = ). + 8(>., x), k-oo k-oo r + 1 ( 4.86) 

where 8(>., x) represents the small deviation due to the forgotten). portion of rand 

input of x. If 8 is very small such that it is negligible when k is large, we say that the 

sequence {c(.)} reaches the quasi steady-state. 

Generally, it is difficult to quantitatively describe 8()', x). Here we model the input 

signal to the systolic array as a second-order AR process described by 

( 4.87) 

where v(n) is a white Gaussian noise process of zero mean and unit variance. By 

choosing of different AR parameters a} and a2, we obtain different realizations of 

the AR process [32]. In our simulations, three different categories of signal are encountered. The first category consists of three stationary AR processes which are 

ARI (a} = -0.1, a2 = -0.8), AR2 (a} = 0.1, a2 = -0.8) with real roots and AR3 

(a} = -0.975, a2 = 0.95) with complex-conjugate roots. The second catogory is a 

non-stationary AR process, AR4 (a} = -0.6, a2 = -0.5), and the third catogory is a 

white Gaussian noise process, WN, with zero mean and unit variance. All of the AR 

processes are normalized to unit variance. Table 4.3 shows the mean distribution for 

different input data with different). values. This table justifies the result in (4.85). 

Table 4.4 shows the variance distribuation of different input data with different ). 

75 
ARI AR2 AR3 AR4 WN . 

= 980 .9800 .9800 .9802 .9799 .9801 

A=.985 .9849 .9849 .9851 .9848 .9850 

A=.990 .9897 .9897 .9900 .9897 .9899 

A=.991 .9907 .9907 .9910 .9907 .9909 

A=.993 .9927 .9927 .9930 .9927 .9929 

A=.995 .9947 .9947 .9950 .9947 .9949 

A=.997 .9967 .9967 .9970 .9967 .9969 

A= 

Table 4.3: !\1ean distribution for different input data with different ,X values 

value. The values of those variances are in the order of 10-4 to 10-6 which implies 

that (; is indeed very small. They can be closely approximated by using quadratic 

polynomials as follows, 

ARI : (1;(,X) - 1.5938 - 3.182), + 1.5882).2 

AR2 : (1;().) = 1.5991 - 3.1919,X + 1.5928,X2 

AR3 : (1;(,X) 1.5812 - 3.1595'x + 1.5784,X2 

AR4 : (1; (>.) - 1.4492 - 2.8936), + 1.4444).2 

AR5 : (1;(,X) = 1.6437 - 3.2904'x + 1.6431,X2 (4.88) 

We can see, although the st.atistics of the input data are different, the variances 

can be described by ,X in a very similar way (see Fig. 4.21). This means, when ,X is 

close to one and the quasi steady-state is reached, the size of the variation {; is mainly 

governed by ). instead of the statistics of the input data. Fig. 4.21 shows the plots of 

the variances in dB scale. 

76 
AR1 AR2 AR3 AR4 WN 

= 980 7.3885e-4 7.5465e-4 6.8163e-4 6.6721e-4 7.3367e-4 

A=.985 4.3970e-4 4.S144e-4 3.9S77e-4 3.9S17e-4 4.3308e-4 

A=.990 2.0903e-4 2. 1463e-4 1.8376e-4 1.8918e-4 2.0080e-4 

A=.991 1.7154e-4 1.7875e-4 1.4883e-4 1.SS62e-4 1.665ge-4 

A=.993 1.0991e-4 1.1390e-4 9.1016e-5 9.6440e-5 1.0323e-4 

A=.995 5.9724e-S 6.0796e-S 4.678ge-S S.18S6e-S 5.3525e-5 

A=.997 2.3007e-5 2.4735e-5 1.6808e-5 1.9908e-5 2.0S04e-S 

A= 

Table 4.4: Variance distributations for different input data with different). values 

With these results, we conclude that sequence {c(.)} reaches the quasi steady-state 

regardless the input statistics if ). is close to one. Thus, we can write 

( 

lim c( k + 1) "" lim Ec( k + 1) == )., k-oo k-oo 

lim s(k + 1) "" lim Es(k + 1) == ,II _ ).2. k-oo k_oo (4.89) 

The quasi steady-state and ensemble values of sine and cosine parameters are the same 

for all boundary cells. It is independent of the statistics of the input data sequence 

and the position of the boundary cell which generates the sine and cosine parameters. 

These results yield the tools needed to further investigate many properties of the 

QRD L5 systolic algorithm. 

77 
I( 

.3o,-----,------r------r----...... 

~Or_---~-----+-~~~-+----~ 

·50 t-----t------+-----+--l.\--~ 

.60 ~-----1.-----L.------L----..J 

0.98 0.99 

I.mbd. 

1.00 

Figure 4.21: Plots of the variances in dB 

4.2 Dynamic Range and Lower Bound on Memory Size 

The dynamic range of the content of boundary cell P Ell can be upper bounded by 

Therefore, 

lim Ir (k)1 < IXmarl ~~. k-oo 11 - ~h _ ).2 

(4.91) 

For internal cell P Elj , we have 

hj(k + 1)1 = Is(k)x(k) + c(k).r(k)1 

78 
- Is(k)x(k) + c(k),X[s(k - 1)x(k -1) + c(k -1)'xr(k - 1)]1 

k i-I 

< L ,Xilx(k - i)s(k - i)1 II c(k - j) 

i=O j=O 

k i-I 

< IXmarl L,Xi Is(k - i)1 II c(k - j) ( 4.92) 

i=O j=O 

From the basic relation between the geometric mean and the arithmetic mean, we 

know 

( 4.93) 

If n is large enough, then from the law of large number, we know 

Therefore, 

i=1 

when n is large. We can further simplify the bound for k --+ 00 by using this inequality 

as follows, 

k 

lim IrIj(k + 1)1 < Ixmaxllim L ,Xis(k - i)E(c(k - i))i 

k ..... = k ..... = i=O 

IXmaxl t ,X2i . VI - ,X2 = IXmaxl = R. 

i=O VI - P ( 4.94) 

From (4.90) and (4.93), we can see the steady-state dynamic range of the first row 

is upper bounded by ?R for both boundary and internal cells. The dynamic range of 

the second row depends on the output of internal cells of the first row. Denote the 

output of the first row as Xouh we have 

Xout(k + 1) = c(k)x(k) - s(k)'xr(k). (4.95 ) 

The first term of the right-hand side of (4.95) can be bounded by 

lim Ic(k)x(k)1 ~ 'xlxmaxl k-= ( 4.96) 

79 
and from (4.94) the second term is bounded by 

( 4.97) 

There are two possible cases. 

Case 1: Highly fluctuated input 

The value of x(k) may vary differently from time to time such that s(k) may 

have an opposite sign of x(k). For this case 

( 4.98) 

Case 2: Smooth input 

For this case, the input data sequence does not change its value rapidly, therefore 

s(k) may always have the same sign as x(k). The bound is 

( 4.99) 

From (4.91) and (4.94), it is obvious the steady-state dynamic range of the second 

row is bounded by 

lim IT2j{k)1 < '\Ixmaxl = 2,\R, 

k-oo - VI - ,\2 ( 4.100) 

for the highly fluctuated input and 

(4.101) 

for the smooth input. From above results, the steady-state dynamic range of the mth 

row is bounded by 

( 4.102) 

80 
for the highly fluctuated input and 

(4.103) 

for the smooth input. 

Fig. 4.22 shows a simulation of the contents of internal cells of the first row and 

the second row. Denote Em as the word-length of the mth row, to prevent overflow and 

to insure the correct operation of the QRD L8 algorithm, we require 2Bm 2:: (2A)m-l~ 

for fixed point operation, and therefore 

Em 2:: r(m - 1)(1 + log2 A) + log2 ~1· (4.104 ) 

The memory size required for each row to prevent overflow increases and decreases 

for the highly fluctuated and smooth inputs respectively. For the fluctuated input, 

when (2A)m-l = 2, one more bit is needed for the memory of each successive row. 

The number of rows m for each increase is 

1 

m = f1 + 1 + log2 A 1, (4.105) 

which is a monotonically decreasing function of A. If A S 0.5, then there is no such m 

exists. That is, the memory size of the array can be fixed at ~ without the overflow 

problem. For smooth input, when Am - 1 = t, one bit can be discarded from the 

memory of the following rows. The number of rows m for each decrease is 

1 m=f1-- 1 ,1, og2/\ (4.106) 

which is a monotonically increasing function of A. For A S 0.5, m = 2. That is, for 

every two rows we can decrease one bit for the memory. 

81 
A PIO! of Ihc CoIMcnu 01 "lbc rUll au»w V~"t' 

. .a. Plot of the Coments of The Second Row (AR I) 

'I~----------------------------------------, 

A Plot of the Contents of The Third Row (AR I) 

o 

Figure 4.22: Contents of the Cells 

82 
A Plot of die Conlall. of1be FInI Row (AR2) 

.. Ir----------------------., 

/'\ ... 

o .... A Plot orlhe Conlents or The Second Row (ARl) 

o ~'I 

... Piol of lIIe ContenlS or The Third Row (AR2) 

o IOD 

83 
A PlO! of die ConICllU of1bc F" ... Row (ARl) 

o,~--------------------------------------~ 

o A Plot of die ConlcnU of The Second Row (ARJI -

aL-__________________________________________ ~ 

n 

It. Piol ohllc ConICnU of The Third Row (ARJ) 

o 

84 
A PlOi of !he ConICllIS of The FUll Row (AR4) 

~~--------------------------------------------~.m 

A Piol of IItc Conlenu of The Second Row (AR4) 

.;~--------------------------------------------., 

o'~----------------------------------~ o 

A Piol of lite Conlcnu of The Third Row (AR4) 

85 
4.3 Stability and Quantization Effect 

In this section, we consider the stability under quantization effects. Here, the stability 

is defined in the sense of bounded input/bounded output (BIBO) as in [60]. From 

(4.98) and (4.99), the output of the mth row is bounded by 

( 4.107) 

for the highly fluctuated input and 

( 4.lOS) 

for the smooth input. 

The order p of the LS problem is finite and fixed. The output of the last row of 

the QR triarray is bounded, in the worst case, by limk ..... = IXoutp I ~ (2,\)P- 1 Ixmax l. The 

resid ual is then asymptotically bounded by 

( 4.109) 

where ,( k) = llf=l Ci( k) [7S], Thus, for ,\ < 1, if the input data are bounded, that is, 

IXmaxl < 00, the output is always bounded. The QRD LS systolic array constitutes 

a BIBO stable system under unlimitted precision condition. In practice, the memory 

size of each processing cell is finite-length. Leung and Haykin [60] first considered 

the stability under this effect and showed the QRD LS algorithm is stable under 

finite-precision implementation. Here we reconsider this problem and obtain a more 1 

general analysis and a tighter bound. 

Let Q(.) be the quantization operator and x be the quantized value of x. Since 

the quantization error for the additions of quantized parameters is much smaller than 

86 
that of the multiplications of them, to make the analysis simplier, we may express 

the quantization error for additions as 

n n 

Q(Lui) = LUi +bn ( 4.110) 

i=l i=l 

From (4.78), the square of the quantized content of the boundary cell is 

k 

f2(k + 1) = Q(Q()..2f2(k)) + Q(x2(k))) = LQ()..2ix2(k - i)) + bk+l' (4.111) 

i=O 

Since the quantization operator Q is a bounded operator such that IQ(x)1 ~ J(lxl for 

all x and some ]{ [60], (4.111) can be bounded by 

, " 2 " 2 " 2k < Amax' x max (1 + >. + ... + >. ), (4.112) 

where xmax is the maximum quantized value of sequence x. The asymptotic behavior 

can be obtained by taking the limit on both sides. Thus, 

(4.113) 

Therefore, the quantized content is 

(4.114) 

With the same arguments as in section 4.1, we then have 

1

· f( k) 1 

1m --+ 

k-oo f(k + 1) (4.115) 

and the quantized steady-state value of cosine is 

. . 5..r(k) hm c(k + 1) = hm "(k ) = ).., k-oo k-oo r . + 1 (4.116) 

87 
and the quantized steady-state value of sine is 

Analogous to section 4.1, we can further obtain limk_oo E c( k) = ~ and limk_oo Es( k) = 

VI - ~2. 

Now consider the quantized content of the internal cell, from (4.92) 

IfIj(k + 1)1 - IQ(Q(s(k)x(k)) + Q(c(k)~f(k))1 

k i-I L Q(~ilx(k - i)s(k - i)1 II c(k - j)) + 8k+1 

i=O j=O 

k i-I 

< /(" maxlxmaxl L ~ils(k - i)1 II c(k - j), (4.117) 

i=O j=O 

where /(" max results from quantization error including 8k+I' From section 4.2, (4.116), 

and (4.11 i), the quantized steady-state dynamic range of the internal cell is bounded 

by 

I· 1- (k)1 < }'" IXmaxl }'" Ii) 1m rlj • _ \ max = \ max;J\.· 

k ..... oo VI _ ~2 

(4.118) 

The output of the mth row is bounded, under the quantization effect, by 

(4.119) 

for the highly fluctuated input and 

( 4.120) 

for smooth input. 

From these results, the quantized asymptotic value of the residual can be obtained 

as 

(4.121) 

Thus, if >. ~ 1 and the input data are bounded, the QRD LS systolic array constitute 

a BIBO stable system under the quantization effect. 

88 
4.4 Finite-length Effect of Fault-tolerant Capability 

In this section, we discuss the finite-length effects of fault-tolerant capability. The 

first problem is that of missing error detection which results from the cumulative 

multiplications of the cosine value with a small error. Under a finite-precision implementation, this may result in a failure of error detection. The minimum word-length 

to circumvent this problem is then derived. The second problem is called the false 

alarm. \\lith the quantization effects, the system without fault may produce quantization error to cause the false alarm. A threshold device is then introduced to tackle 

this problem. 

4.4.1 Missing Error Detection 

By missing error detection we mean that a small error generated by a faulty processing 

cell is not detected due to the finite-precision computation. Assume a fault occurs 

in an internal cell P Eij , i #- j, at a faulty moment. The output of this faulty cell is 

thus erroneous and can be described by X~ut = Xout + 6, where Xout is the fault-free 

output and 6 is the error generated by the fault. The error propagation path can be 

described by 

P Ei} ---? P E(i+l)j ---? ••• ---? P E}}, 

and then P Ekl , k ~ j, I ~ j are all contaminated. From the operations executed 

by the internal cell, the error is modified to Ci+16 by P E(i+l)j and the cumulative 

89 
modifications of the error before reaching the boundary cell, P Ejj , is 

i-I 

17 = fJ II Ck, 

k=i+1 

( 4.122) 

where Ci is the cosine parameter generated by the boundary cell P Eii . Let cj and sj 

denote the erroneous Cj and Sj respectively. The cj and sj are then given by 

s'. = Xin + 1] 

] J)..2 r 2 + (Xin + 17)2 

( 4.123) 

In this case, sj is no longer proportional to Xin, Q(J) will not be zeroed out by the ph 

cell of the EDA. The size of the error generated by this cell is 

(4.124 ) 

where r' = J }..2r2 + x7n is the new updated uncontaminated value of the content of 

P Ejj . When 1]j propagates down to the output of the EDA, 1]j is influenced by the 

contaminated cosines c' of each following row. The error output at eo due to an error 

8 generated at P Eij is then given by 

p p 

o (. .) fO Z,) = II -I c, II ' m 1]j = -I cm 1] 

m=j+l m=) 

i-I p 

-I II Ck· II c:n8, (4.125) 

k=i+l m=] 

I j-l p 

eg(i,j) = - II CI II cZ II c~8. (4.126) 

1=1 k=i+l m=i 

Next, assume a fault occurs in a boundary cell, P Eji , 1:S j :S p, at the faulty 

moment. Both erroneous cj and sj produced by P Eii can be written by 

, 

c = ) r' ( 

90 

( 4.127) 
where 8e and 88 represent errors in the numerators while r~ represents the erroneous 

content of the denominators of Cj a nd 8j. The error produced by the ph cell of the 

EDA is then given by 

( 4.128) 

and the output error at eo due to a faulty boundary cell is given by 

6(' .) rrP

' Xin8e - >'r8s eo },} = 1 C' 

m r' 

m=j+l ( 

i P rr CI' rr c~. r/j. ( 4.129) 

1=1 m=j+l 

From (4.126) and (4.129), we can see that eg i 0, under unlimited precision condition, 

if there is a fault occurring in the system, except when Uin8e = >.r8s in (4.128). 

However, this is unlikely to happen. From Lemma 3.3, for 0 < Cj ~ 1, the error 

may not be detected after multiple multiplications of Ci in (4.126) and (4.129) under 

finite-precision implementation. It is obvious there is no such problem when 8 is 

large. Since r in (4.123) tends to be a large number asymptotically, it is reasonable 

to assume the error size 8 generated by a fault is much smaller than r when 8 is small. 

Under this circumstance, from (4.123), we have cj ~ Cj. In the quasi steady-state, 

the asymptotic behavior of erroneous cosine is cj ~ Cj = >.. From (4.126) and (4.129)' 

the error output eg due to an error size 8 is then approximated by 

6(' .) C:1 _ \2p-i( eo Z,} - A U (4.130) 

for a faulty internal cell and 

(4.131) 

for a faulty boundary cell. Denote BL::. be the word-length of each memory of fixed 

point arithmetics. That is, each memory size is of B L::. bits and let l::. = min( 8,77]). 

91 
To insure the detection of error of size 6, we need 

(4.132) 

Therefore, the memory size should be at least 

( 4.133) 

such that the small error size can be detected. The second term of the right-hand 

size is obvious since the error size 6 must be detected; the first term is to account 

for the effects that the error propagates through the array. 

4.4.2 False Alarm 

Due to the finite-precision implementation, the residual output of the EDA will not 

be an actual zero if there is no fault in the system. We call this effect a false alarm. 

Here, we are going to model and quantitatively describe the false alarm effect and 

introduce a threshold device to overcome this problem. 

Cancellation Principle 

Suppose now we have a QRD RLS array of order p = 3. Denote the first and second 

where the checksums Xl +X2+X3 and x~ +x;+x; are inputs to the EDA. Fig. 4.23 shows 

the first two rows of the array. After both data pass through the array, according to 

the operations of the processing cells, the contents of the cells of the first row are 

92 
(( 

Figure 4.23: The first two rows of the array 

(4.134) 

where c = xdrll and s = xUrll are the rotation parameters generated by the 

boundary cell and r,j is the content of P Eij • The output of the internal cells are 

Z)3 

(4.135) 

Since sx; + ex) = J x? + X/~ and ex; - SXl = 0, we have r14 = ru + r12 + r13 and 

Z14 = Z12 + Z13' That is, both the contents and the outputs of the first row still meet 

the checksum. The output of the first cell of EDA, Z14, can be rewritten as 

(4.136) 

We can see that the data from the first column got cancelled out by the first cell of 

the EDA. Since the outputs meet the checksum, with the same principle, the data 

93 
from the second column will got cancelled out by the 2nd cell of EDA. Thus, this 

observation can be generalized and stated as bellowed: 

Cancellation Princeple: With the checksum encoding data inputted to EDA, the 

data from the ith column got cancelled out by the ith cell of the EDA.D 

For a finite-precision implementation, due to the roundoff error, the data from the 

ith column will not be totally cancelled out by the ith cell of the EDA. This effect 

results in the false alarm problem. 

Floating Point Arithmetics 

A floating point number f can be represented by [30] 

( 4.137) 

where f3 is the base, t is the precision, and [L, U] is the exponent range. The floating 

point operator fl can be shown to satisfy [30] 

x = fl(x) - x{1 + f), 

fl(a op b) (a op b)(1 + f) kl:::; u, (4.138) 

where u is the unit roundoff defined by 

_ 1 f3 1- t u - - , 2 

for rounded arithmetics, 

and op denote any of the four arithmetic operations +, -, x,"';-. 

Roundoff Analysis 

For a QRD RLS array of order p with floating point arithmetics, denote the first row 

94 
and If I < u is a constant2, and the second row of input vector as (x~, x;, .. " x~, Ef=l xi+ 

f,,). The content of the first boundary cell is given by 

~ 2 ~ ,2( 1 ) Xl + X I + f , ( 4.139) 

and the rotation parameters are e = Jl(xdrn) and s = Jl(xUrn). The contents of 

the internal cells then can be obtained as 

T-ij Jl(Jl(sxj) + Jl(cxj)) 

[sxj(1 + f) + cXj(1 + f)](1 + f) 

~ (1 + 2f)(SXj + eXj), 1 < j S p, ( 4.140) 

and 

" " 1\,1'+1 J1(fl(8C[:. x: + fp)) + Jl(c(2:. Xi + fp))) 

i=l i=l 

I' " 

~ (8 2:. x~ + c 2:. xd + 6f". (4.141) 

i=1 i=l 

The mismatch 7"1 resulting from the finite precision computations of the first row is 

7"} = 6fp - (fJxi + Xl~ + 2f t(8Xj + eXi)) ( 4.142) 

i=2 

and it is bounded by 

hi < 6plfXmax i + 12fXmax i + 4(p -l)lfxmaxl 

= (lOp - 2)1fXmax i S 10plfxmaxl· ( 4.143) 

For the second row, using the same approach, the mismatch is bounded by lO(p -

1)lfxmax l. The total mismatch from the whole array for a given input vector is given 

2To simplify the notation, we do not give indeices to different ('s. 

95 
by 

p-l 

ITI ~ L lO(p - i)lfXmaxl = 5p(p + 1)lfxmax l· (4.144) 

i=O 

Let n denote the number of input data rows. As pointed out in [109], for the Givens 

rotation method, the roundoff error is proportional to O(n1.5). The possible mismatch 

becomes 

(4.145) 

This bound can be interpreted as: for each row of input, each processing cell contributes about Itxmaxl amount of roundoff error. Since there are about p(p + 1) 

processing cells, the total possible roundoff error is then p(p + 1) Itxmax I. 

In order to prevent false alarm, the threshold th has to be set at least at IT{n)l. 

Suppose f3 = 2, t = 16, then u = 2-16 . Given an scaled input data such that 

IXmaxl = 1 and n = 104 , the threshold of a QRD RLS array of order p = 100 is 

( 4.146) 

4.4.3 Overall Memory Size Consideration 

To prevent missing error detection, we want to set the detectable error size l::" 

min( 8, 7]j) as small as possible. While to prevent the false alarm, we also want to 

choose a sufficiently high threshold. Both situations cannot be satisfied simultaneously since both goals are conflicting. 

To detect the error size l::", from (4.130), (4.131), and (4.132), we set the threshold 

th ~ )..2pl::". That is, the minimal detectable error size is )..-2P. tho From (4.133), 

(4.147) 

96 
A criterion to choose B6 is then given by 

B6 = min( r -2plog2 A -log2 6.1, r -log2 th 1)· ( 4.148) 

To prevent overflow, from (4.104), set 

(4.149) 

For a QRD RLS systolic array to detect an error size of 6. without false alarm and 

overflow problems, the memory size B is required to be at least 

(4.150 ) 

97 
Chapter 5 

Order Degraded Performance and 

Residual Estimations 

When the faulty row is found in the fault-tolerant QRD RL5 systolic array, it enters 

into an order degraded operation. In this chapter, the performance degradation of 

the reduced-order QRD L5 is studied. The optimal residual estimation under faulty 

situation is also considered in this chapter. 

5.1 Order Degraded Performance 

Consider a order p L5 problem with a n x p complex-valued data matrix Ap( n) denoted 

by 

Ap(n) = [1!(1),~(2), ... ,1!(n)V = [Q.(1),g(2),··· ,g(p)] = [Ap_dn) : g(p)], (5.151) 

a n x 1 desired response vector 

u(n) = [d(1),d(2), ... ,d(n)1T, 

98 
a p x 1 weight vector 

(5.152) 

and a n x 1 residual vector 

.{" ( n) = [e" ( 1 ), e" (2), ... , e" ( n )] T = A" ( n ) w" ( n) - 1L ( n ) . (5.153) 

Let the index of performance be defined by the weighted 12 norm of 

(5.154) 

where 

with a real-valued forgetting factor 0 < A ~ 1 and A = A 2 . Then the least-squares 

solution, satisfies 

(5.155) 

The optimal weight vector can be obtained by solving the normal equation 

(5.156) 

where 

H - A,,_l (n )AA,,_l (n) 

(5.157) 

99 
and 

(5.158) 

It can be easily shown that 

(5.159) 

For a order p - 1 LS problem, then as before, we want to minimize the weighted 12 

norm of 

(5.160) 

with weihgt vector 

Obviously, the optimal weight vector of order p and p-1 can be obtained as wp(n) = 

¢;l(n)tjJp(n) and t'lJp_l(n) = ¢;~l(n)1jJp_l(n) respectively. The difference vector of the 

optimal residual vectors of different order is defined by 

(5.161) 

and the weighted 12 norm of 6( n) is defined to be the order degraded performance, 

r( n), given by 

r(n) - 116(n)ll~ 

- 6wH(n)¢p_l(n)6iv(n) + IIw;(n)112 '11!l(p)lI~ 

+2· Re[w~(n)6iQH (n)A:_1AQ(p)], 

100 

(5.162) 
To relate wp_llp(n) with .Ylp_l(n), we define a LS problem of order p - 1 with 

Ap- 1 as the data matrix and g,(p) as the desired response. That is, we would like to 

minimize the weighted 12 norm of the residual vector 

From (5.153) and (5.156), the optimal weight vector can be obtained by solving 

(5.164 ) 

where fl(n) is defined in (5.157). Let the optimal index of performance of this LS 

problem be ~a,p-l(n) = 11.k,p_1(n)II~, then </>;1 can be represented by 

aT -1'-1 a 

+) 1 [-~,p-1(n)l[-'H () ~p-l n ~a,p-1(n 1 ' 

Then wp ( n) can be represented by expressions of order p - 1 as 

and therefore, 

</>-1 (n)AH (n) + .l£a,p_l(n)r:,p_l(n) 1'-1 1'-1 {o,p_l!n) 

-~p_l(n) 

~o,p-l (n) 

. 7(n) . W p _1(n) + ~ ()1lla,p-l(n), 

7(n) 

~a,p-1 (n)' 

101 

a,p-l n 

1]. (5.165) 

(5.166) 

(5.167) 

(5.168) 
where I(n) =< la,p-1(n),M.(n) >A and < ~'M. >A=< A~, AM. > is the weighted inner 

product. Thus, 

1\ A () I( n) A () 

ulQ n = t ()llia,P-1 n . "'a,p-1 n 

(5.169) 

Now, we may proceed to calculate r(n) in (5.162). The first term of (5.162) becomes 

b.wH (n)¢p_dn)b.w(n) 

II I ( n ) 112 A () AH ( ) A- () A () - e ().llia,P-1 n p-1 n Ap_1 n .llia,p-1 n a,p-1 n 

III(n)1I2 A H- A 

= e ()(.~a,p_1(n)+!!.(p)) A(~,p_1(n)+!!.(p)) a,p-1 n 

III(n)1I2 A 

- e () (ea,p-1(n) + 7](n) + 2· Re( < ~,p_1(n),!!.(p) >A)). (5.170) a,p-1 n 

The second term of (5.162) becomes 

(5.171 ) 

The third term of (5.162) is 

Combining (5.170), (5.171), and (5.172) together, we have 

f(n) = III(n)W = II < fa,p-1(n),y(nl >/\ 112. 

ea,p-1(n) Illa,p-1(n)IIA (5.173) 

Denote the last row of input data matrix y(n) = [yp_1(n),up(n)], the difference of 

the optimal residual at time n then can be obtained as 

102 
III(n)1I2 II T ( ) A () ()11 2 

- (2 () 1!,,-1 n lQa,,,-l n - u" n 

'>a,,,-l n 

III(n)1I2 A 2 - (2 () Ilea,,,-l(n)1! . '>a,,,-l n 

(5.174) 

5.1.1 Geometric Interpretation 

From (5.173), we can see the order degraded performance is indeed the energy of the 

projection of the desired response 1!.(n) onto the subsapce spanned by the optimal 

residual1",,_I(n). As the vector lL(n) becomes more orthogonal to 1",,_l(n), the 

order degraded performance is also reduced. Denote the column space of A,,_1 (n) by 

{Ap_dn)}. Then the projection operator PAp_1 projects vectors onto space {Ap_l(n)} 

and the orthogonal projection operator Pip _ 1 = I - P Ap - 1 projects vectors onto the 

space {A;_1 (n )} which is orthogonal to space {Ap - 1 (n) }. The entire space S spanned 

by lL( n) can be represented by 

(5.175) 

and all of these subspaces are orthogonal to each other. It is obvious that 1",,_l(n) = 

Pip _ 1Q.(p)· By projecting the desired response 1!.( n) to these subspaces, we obtain 

1!.(n) PAp _11!.(n) + P1a,p_l1!.(n) + P}p1!.(n) 

A () • () ( ) 1",,-1 (n) A ( ) 

= lL n + < £',p-l n ,lL n >A 111"p-l(n)ll~ +!.p n , (5.176) 

and all of these vectors are also orthogonal to each other. If we drop the vector Q(p), 

then the one-dimensional subspace of {Pip- 1Q(p)} cannot be used to represent lL(n). 

Therefore, the components of '}!,.{ n) in this subspace is lost and this introduces an error 

vector P'{ y(n) with energy .... ,1'-1_ 

IIPt (n)1I2 = III(n)112 = f(n), _,p-IlL ((n) ,>a,,,-l 

(5.177) 

103 
((£r 

f-I .~ 

· · · 

· · · · 

· 

Figure 5.24: Geometric illustration 

The energy of the last component of this error vector is given by (5.174). Fig 5.24 

illustrates the geometric interpretation discussed above. 

5.2 Residual Estimation in Faulty Situation 

From Section 4.4, we know that if there is a fault occurring in the system, the error 

output at EO due to an error D generated at an internal cell P Eij is given by 

p 

e~ (i ,j) = - T II c'm 1] , (5.178) 

m=) 

where 1] == D n~::+l Ck, and the error due to a faulty boundary cell is given by 

) p 

e~(j,j) = II Cl' II c~ .1]j, (5.179) 

1=1 m=j+l 

where 1]j = (XinDc - )..rb~)/r: with be and b~ representing errors in the numerators 

while r~ represents the erroneous content of the denominators of Cj and Sj. Whenever 

eg( n) :f. 0, we detect there is a fault occurring in the system. The optimal residual 

104 
e(n) is thus erroneous and is denoted as e6(n). In this section, we devise a method 

to estimate the true optimal residual e( n) given the contaminated residual e6 (n). 

5.2.1 Faulty Internal Cell 

Once the faulty cell P Eij is identified, then T} in (5.178) can be estimated by 

(5.180 ) 

where c~s are obtained at the right hand side of the systolic array. If T} is large, then 

through many non-linear operations in the cells, it is essentially impossible to keep 

track of the effect of T} on e6 (n). Thus, for analytical tractability, we assume the error 

T} is small in the sense that T}X P « r. With this assumption, the erroneous cosine can 

be approximated by 

'xr 

J,X2 r2 + (x + T})2 

'xr 'xr 

-;;::::;;;===;;:== '" --- Vr,2 + 2T}x - r' + T}X 

c) T}x 

1 + !E ~ Cj - -:;:;- = Cj - T}Sj 

r' 

and the erroneous sine is given by 

(5.181) 

(5.182) 

And since T}xP « r, we also have T}Sj « 1 and T}SjCj « Cj < 1 which yield cj < 1, sj < 

1, and C'/ + sj2 ~ 1. 

The erroneous cosine and sine generated by the boundary cell P Ejj are then sent 

to the right for processing by the internal cells P Ej,)+kl 1 :s k :s p - j. The output 

105 
of these cells are erroneous and are given by 

(5.183) 

where x is the input data and Xout and X~ut are the uncontaminated and erroneous 

output respectively. Denote the error size produced by cell P Ei,j as ~i,j. With this 

notation, ~j-1,j = Tj, and the error size produced by cell P Ej,j+k is 

(5.184) 

Recall that O(r) '" O(x/Vl -.\2) and S '" VI - -\2, we have O(~j,j+k) '" O(TJx). As 

discussed above, the contaminated rotation parameters generated by the boundary 

cell P Ej+k,j+k. are given by 

{ 

cj+k = Cj+k - ~j+k-l,j+kSj+k 

sj+k = S)+k + ~j+k-l,)+kSj+kCj+k 1 ~ k ~ p - j. 

(5.185 ) 

By solving this linear system of equations, we can estimate the uncontaminated sine 

and cosine as 

{ 

. '" 1( , 0,-1 + 1(--' 0,-1 )2 +4 ' C)+k - 2 c j +k - ~j+k-l,j+k V Cj+k - ~j+k-l,j+k Sj+k' 

Sj+k ~ )1 - cJ+k 0 ~ k ~ p - j. 

(5.186) 

When (::r;~k-1,j+k « 1, we can further simply (5.186) as 

(5.187) 

To obtain ~j+k,j+" 1 ~ k ~ p - j, k ~ I ~ p - j, from the operations of the internal 

cell 

106 
- (Ci+k - ~i+k-l,j+kSj+k)(X + ~j+k-l,i+I) - (Sj+k + ~j+k-l,j+kSj+kCj+k) . >'rj+k,i+/ 

- XoutJ+k,J+1 - (Sj+k~j+k-l,j+k - Cj+k)~j+k-llj+1 - ~j+k-l,j+kSj+kCj+k>'ri+kJ';>~~88) 

A recursive equation is then obtained to compute ~j+klj+/ from ~i+k-l,j+1 and SJ'j+k-l,j+k 

as 

~j+klj+/ 

1 ~ k ~ p - j, k < 1 ~ p - j, (5.189) 

where Pj+k = Si+k~j+k-l,j+k - Cj+k. The order of ~j+k,j+/ is that of 

O( ~ j+k-l,j+kS j+kCj+k Ar j+k,j+/)' 

Therefore, O(~j+l,j+d '" O(1]X2), We can also show that O(~j+k,J+d '" O(1]Xk) < 

O( 1]XP ). Thus the error size is consistent with the original assumption of the analyizable situation. 

With (5.189), we can recursively compute the error size produced by the cell 

P Ep,p+l and that of ~p,p+l' The erroneous residual can be obtained as 

(5.190) 

The optimal uncontaminated residual is then estimated as 

(5.191) 

To update the content of P Ej+k,j+k, we try to keep it as uncontaminated as possible 

by performing 

rj+k,j+I(n + 1) - Sj+k(x + ~j+k-l,j+/) - Cj+kArj+k,i+I(n), 

1 ~ k ~ p - j, k < 1 ~ p - j. 

107 

(5.192) 
We summarize the recursive estimation procedure as follows, 

1. Estimate error size TJ = -eg(n)/b TI~=i c~) from the EDA. 

2. Estimate the uncontaminated rotation parameters by using 

{ 

Cj+k ~ HC~+k - ~j~k-l,j+k + -/(S+k - ~j~k-l,j+k)2 + 4sj+k 

Sj+k ~ -/1 - CJ+k 0 ~ k ~ p - j. 

3. Recursive estimation of the error size ~j+k,j+1 by computing 

4. Estimate the output residual from 

5. Update the contents of the processing cells by doing 

1 ~ k ~ p - j, k < I ~ p - j. 

5.2.2 Faulty Boundary Cell 

For a faulty boundary cell, from (5.179), we can estimate 

TJj = TIj TIP -J2 

1=1 CI' m=J+1 (;m 

(5.193) 

Unfortunately, there is no way to estimate 8c and 85 from TJj. Therefore, we are unable 

to estimate the optimal residual if the faulty cell is a boundary cell. 

108 
Chapter 6 

Multi-phase Systolic Algorithms 

for Spectral Decomposition 

In this chapter, we propose two multi-phase systolic algorithms to solve the spectral 

decomposition problem based on QR algorithm. The spectral decomposition constitutes one of the most computationally intensive needs of modern signal processing 

applications. \Vhile the QR algorithm is well known to be an effective method to 

solve the eigenvalue problem, there is still no single systolic array architecture that 

can compute the unitary Q matrix readily and perform the QR algorithm efficiently. 

Previous methods based on the Jacobi-like approach required global communication 

or broadcast in computing the eigenvector, and methods using the QR algorithm 

had communication problems among different architectures. In this chapter, we show 

that the Q matrix can be computed easily by using multi-phase systolic algorithms 

and thus the eigenvectors can also be computed without any global communication 

in the array. Two arrays, a triangular and rectangular, are presented to implement 

109 
the multi-phase algorithms. Details on these multi-phase operations of the QR algorithm as well as architectural consequences are discussed in the chapter. Efficient 

fault-tolerant schemes for these multi-phase operations are also considered. 

Recent developments in the parallel processing architectures, especially the systolic architectures, for the spectral decomposition are discussed in section 6.1. In 

section 6.2, some preliminary matrix operations useful for the multi-phase operations 

are discussed. In section 6.3, we review the QR algorithm and show the evaluation of the eigenvector from cumulative multiplications of the Q matrices. Then two 

multi-phase systolic arrays for the QR algorithm and the Hessenberg reduction are 

presented in section 6.4. Their performances, numerical stabilities, and convergence 

rates are studied in section 6.5. Finally, some efficient fault-tolerant schemes that can 

be incorporated with the arrays are discussed in section 6.6. 

6.1 Recent Developments 

Computing the spectral decomposition of a matrix is an important issue in many modern signal processing and system applications. The feasibility of real-time processing 

for sophisticated modern signal processing systems, depends crucially on efficient implementation of parallel processing of the algorithms and associated architectures 

needed to perform these operations [14, 58]. While many variations exist in the literatures for solving these matrix problems, the heart of all these iterative methods 

are based either on the Jacobi-Hestennes method or the QR algorithm [30, 101, 107]. 

While there are some fundamental differences between these two approaches, both algorithms have good numerical stability and convergence rate properties and thus are 

desirable for possible implementation. Since present VLSI technology is capable of 

110 
building a multiprocessor system on a chip, many researchers have proposed different 

parallel processing architectures to solve eigenvalue and singular value decomposition 

(SVD) problems. 

For any complex-valued m x n matrix A, the classical spectral decomposition [102] 

of the n x n Hermetian matrix AHA, is given by 

n 

AH A = 'LAjlljllf = VAVH, (6.194) 

j=l 

where V = U!1"",31nJ is an n x n unitary matrix, A = diag[Al,···,An}, and His 

the complex conjugate transpose operator. The A~8 are the eigenvalues satisfying 

>'1 ~ >'2" . ~ >'n ~ 0 and the v:s are the eigenvectors satisfying AH A31j = Ajllj. The 

decomposition of AH A follows from the SVD [30J of A given by 

(6.195 ) 

where U = [lil"" ,linJ is an m x n matrix with orthogonal column vectors,S = 

diag[sl' ... , sn], and V is an n x n unitary matrix. The 8~8 are the singular values 

satisfying 81 ~ 82' ~ Sn ~ 0 and are the positive square roots of A~8 such that 

A = 52. In this paper, we shall use spectral decomposition in the broad sense of 

not only including the decompositions of (1) and (2), but also include the eigenvalue 

decomposition of an arbitrary complex-valued n x nmatrix A given by 

AX = AX, (6.196) 

where X is an n X n matrix of eigenvectors and A = diag[>'ll ... , >'nJ is the matrix of 

eigenvalues of A. 

Luk [72], Brent [12], and Gao and Thomas [27], have used effectively the J acobilike method to solve these problems for either a multiprocessor system or systolic 

111 
array. The basic problem concerns the diagonalization of a 2 x 2 matrix by the 

rotation matrices J(O) and J((cP) in 

[ 

w x 1 [ d

1 

J(Of y z J((cP) = 0 ;,1 (6.197) 

A two-stage procedure is then used to find 0 and cP [72]. To find the SVD of a square 

matrix A, an appropriate sequence of 2 x 2 matrices is computed by using the basic 

Jacobi transformation in 

(6.198) 

where Jij and J(ij are rotations in the (i, j) plane chosen to annihilate the (i, j) and 

(j, i) elements of A respectively [72]. While the Jacobi-like method, as considered in 

[72], is currently known as one of the most effective parallel SVD algorithm for full 

dense matrices, the computations required to obtain the rotational matrices needed 

in this approach to obtain the singular vectors are not simple and can not be obtained 

without broadcast [72]. Moreno and Lang [82] also considered some alternatives to 

the algorithms in [12]. 

On the other hand, other researchers have used QR algorithm to solve the eigenvalue problems. Heller and Ipsen [33, 43] performed the QR iteration for banded matrix based on orthogonal systolic network and Schreiber [94] combined their network 

with Gentleman-Kung's QR array to compute the QR algorithm. These methods 

required the computation of the unitary matrix Q. However, problems exist in the 

concurrent computation of Q and the pipeline operation of the QR iteration [43]. In 

[81], Moldovan et a1. studied the mapping of a large QR algorithm onto a fixed size 

array. Torralba and Navarro [103] further purposed a size-independent linear array for 

QR iteration and Hessenberg reduction. While this approach can provide an efficient 

112 
computation of one iteration of the QR iteration, it is not obvious how to pipeline 

the iteration. 

For some system applications, such as matrix rank determination and system identification [52], the efficient computation of singular values is sufficient, while in other 

applications such as antenna beamformation [79, 97], spectral estimation [50, 93], 

direct finding [31, 85], etc., the eigenvectors are crucially needed. This makes practical implementation of systolic arrays discussed above difficult for many applications 

since they either cannot compute the eigenvector or cannot obtain the eigenvector 

without broadcast. For example, for the MUSIC algorithm [32], once we determine 

the signal subspace and noise subspace from the eigenvalues, the sample spectrum is 

then determined by 

1 

S(w) = sH(w)XNX/Js(w)' 

where X N is the matrix of eigenvectors which generate the noise subspace and 

s(w) - [1 e-jw ... e-jw(K-l)] -, " , 

with J{ is the dimension of matrix XN. A system which consists of several systolic 

modules to compute the MUSIC algorithm has been proposed in [90]. However, communication problems among the modules and the difficulty of matching the pipeline 

rates and timings among different modules may pose difficulties for practical implementation. 

Presently, there is no known simple efficient systolic array approach for the generation of eigenvectors. The main reason is that there is no single architecture that is 

capable of handling all the steps required in the algorithm such that we can pipeline 

the successive iteration readily. The communication cost among different architectures is high and the interface problem for an efficient data flow is demanding. In this 

113 
paper, we propose two multi-phase systolic algorithms to solve the spectral decomposition problem based on the QR algorithm. By multi-phase operations we mean that 

the processing cells can perform different arithmetical operations in different phase 

of the computations. Two systolic arrays, one is triangular and the other is rectangular, are designed based on the multi-phase concept. A key feature in our method 

for the successfully application of the QR algorithm is that the Q matrix of the 

QR decomposition can be computed explicitly by multiphase operations. With the 

proper feedback of this Q matrix, the QR algorithm can be computed and pipelined 

effectively in a single systolic array. From the accumulation of those Q matrices in 

another array, eigenvectors and singular vectors can be computed without needing 

global communication inside the array. 

6.2 Systolic Array Matrix Processing 

In this section, we consider some preliminary matrix and associated systolic array 

operations needed in the multi-phase systolic algorithms for spectral decompositions. 

A. QR Decomposition 

A non-degenerate m x n rectangular matrix A can be factorized into two matrices 

Q and R such that A = Q R, where Q is an m x m unitary matrix and R is an m x n 

upper triangular matrix. The matrix Q can be computed using sequences of Givens 

rotations. An elementary Givens transformation has the form of 

o 0 

r~ 1 x~ 

(6.199) 

114 
Figure 6.25: Triangular systolic array for QR decomposition. 

where 

ri 

Several different QR arrays have been considered by Gentleman and Kung [26], Heller 

and Ipsen [31], and Luk [67]. In particular, the computation of the Q matrix without 

broadcast is difficult for the array considered in [66, pp.266]. On the other hand, [26] 

has shown that a triangular systolic array can be used to obtain the upper triangular 

matrix R based on sequences of Givens rotations. This approach also leads to an 

efficient method for performing recursive least-squares computation [72], and is also 

useful for finding the singular value of a matrix [25]. This systolic array is shown in 

Fig. 6.25 and the operations of the cells are described in the first column (i.e., phase 

1) of Table 6.5. While the rotation parameters are propagated to the right, the Q 

matrix will not appear directly at the right as originally suggested by [86]. In order to 

115 
Phase 1 Phase 2 Phase 3 

2rl If Z," = 0 then 

c-li .-0; 

otherwise 

r' = J).2r2 + z~"i • - z, .. /r • - z,,,r 

c - ).r/r'; • - z,,,/r' 

r _r'; 

end 

:t,,, 

Zwl - cz,,, - .).r 

.wl - .,,, + zt"r .z," + c).r 

Zwl - z," - .r 

r - :tout 

Table 6.5: Operations of the processing cells for different phases. 

demonstrate this point, denote Gii as the Givens rotation matrix of the (i,j) plane, 

(( 

then matrix Q can be obtained as 

(6.200) 

i=m-l i=i+l 

where U is an ordered matrix product such that UJ=m-l Gi = Gm - 1Gm - 2 • •• Gil while 

n denoted a conventional ordinary product. From Table 6.6, we can see, for an xn QR 

triarray, the first rotation parameter coming out at the right edge occurs at time n+ 1. 

After that, rotation parameters for different plane rotation come out successively. If 

assuming that the operation of n discussed above can be obtained immediately, then 

there are m operations of U when all of n are available. This observation leads to 

the conclusion that we cannot obtain the Q matrix by cumulatively multiplying the 

rotation parameters propagated to the right edge. Thus, This is not an effective way 

to obtain the Q matrix. 

116 
Time n+l n+2 n+3 n+4 n+S n+6 

First (1,2) (1,3) (1,4) (1,5) (1,6) (1,7) row 

Second (2,3) row (2,4) (2,5) (2,6) 

Third (3,4) (3,5) 

row 

Table 6.6: The timing table for the rotation parameters to reach the right edge of the 

QR triarray. 

B. Computation of R-T ~ 

In [7], Comon and Robert presented a rectangular systolic array for the computation of B-1 A, where B and A are square and rectangular matrices respectively. The 

computation takes two phases. First, the matrix B is fed into the array and B-1 

is computed. In the second phase, the matrix A is inputted to produce B-1 A. For 

the special case where B is an upper triangular matrix denoted by R, instead of a 

full dense matrix, Mc\Vhirter and Shepherd [73] used the property that a triangular 

array can compute R-T L in one phase with the matrix R prestored in the triarray. 

Since this property is needed in phase 2 of our work, and no derivation was given in 

[73], we present a brief derivation of this result. Define rij = (R)ij and r:j = (R-l )ij, 

where TI) = 0 and rL = 0 for i > j, then it can be shown that 

t =), 

(6.201) 

i < j $ n. 

117 
Let [Yl,·· . ,Yn]T = R-T ~, then 

j 

Yj = L:xir~j' j = I,···,n. 

i=1 

In particular, Yi can be expressed in terms of rij as 

Yj 

I j-l j-l 

- . (Xj - L: Xi L: r~krkj) rjj 1=1 k=i 

1 j-l k 

= -. (Xj - L: L:xir:krkj). 

rjj k=1 i=1 

From (6.202), Yj is given by 

(6.202) 

(6.203) 

(6.204) 

Thus, Yi can be computed recursively according to the above equation in the following 

algorithm: Fig. 6.26 shows the data flow of the input ~ and the output '#.: 

Recursive Algorithm for Computing Ji. = R-T:r.. 

1 Yl = ;;7 . Xl 

for j = 2 to n 

begin 

end 

Zj = Xj 

for k = 1 to j - 1 

Zj = Zj - Ykrkj 

Yj = z)/rJj 

The corresponding systolic array to implement the above algorithm is the same 

as the one shown in Fig. 6.25. The operations of the cells are shown in the second 

column of Table 6.5. The first part of the (6.204) , (i.e., the division), is performed 

by the boundary cell while the second part of (6.204) is cumulated by the internal 

cells. With R pre-stored in the triarray, Fig. 6.26 shows the data flow of the input :r.. 

and the output Ji.. 

118 
R 

Figure 6.26: Computation of R-T x using a triarray. 

C. Triangular-Matrix Multiplication 

The multiplication of a triangular matrix R and an rectangular full dense matrix 

B is gi ven by 

n 

(C)ij = (RB)ij = 2: rikbkj, (6.205 ) 

k=i 

where rik and bkj are elements of matrices Rand B. Using the same array as in 

Fig. 6.251 with R prestored in the triarray and the operations shown in the third 

column of Table 6.51 this multiplication can be easily obtained if B is inputted row 

by row as in Fig. 6.27. 

D. Matrix Multiplication 

There are many ways to implement a full matrix-matrix multiplication in systolic 

array [52]. In Fig. 6.28 1 we show a typical architecture that can be incorporated with 

119 
I( 

~1 

R 

Figure 6.27: Multiplication of a triangular matrix and a full dense matrix. 

the multi-phase operations to obtain eigenvectors. With input matrices Q and A 

arranged as in Fig. 6.28, the matrix BT, where B = AQ, will sit in the rectangular 

array when the computation is completed. Details on this issue will be discussed in 

later sections. 

6.3 QR Algorithm 

In this section we review briefly the basic operation of the QR algorithm and show 

the evaluation of the eigenvectors from the cumulative multiplication of successive 

Q matrices. For a complex-valued n x n matrix A, it states that there is a unitary 

transform U such that R = U AU H is a upper triangular matrix with diagonal eigenvalues of descending order. This follows from the QR Algorithm [27, 92, 97] where 

120 
Figure 6.28: Matrix-matrix multiplication in a rectangular array. 

with unitary Qk and upper triangular Rk. Furthermore, Ak converges to the upper 

triangular matrix with diagonal eigenvalue elements. However, it is not obvious how 

to compute the eigenvectors from those Qk and Rk we have calculated. With similar 

derivations as in [97], here we shows how to obtain the eigenvector associated with the 

largest eigen value from cumulative multiplications of Q k. From the above discussions, 

we have 

(6.206) 

Define 

k 

th U Qi = Ql Q2 ... Qk, 

i=l 

121 
Then we have 

1 

Rk - URi = RkRk-I ., • RI • 

i=k 

Thus the multiplication of QkRk can be expressed as 

(6.207) 

(6.208) 

(6.209) 

Let the eigenvalues of A satisfy, I>'d > 1>'21 > ... > I>'nl. Denote the matrix 

eigenvectors and eigenvalues of A by X and A respectively. Then Ak is given by 

(6.210) 

Let the QR decomposition of X be X = QR and the LU decomposition of X-I be 

X-I = LV, where L is an unit-lower triangular matrix. Then 

(6.211) 

where 

(6.212) 

and 

(6.213) 

Therefore, when k is large enough, we have limk ...... oo Ek = 0 and thus A k LA -k approachs the identity matrix. Then (6.211) can be rewritten as 

(6.214) 

122 
Since the term RAkU is an upper triangular matrix, comparing to (6.209) we can 

see that Ok ~ Q when k is large. That is, the Q matrix of the QR decomposition 

of Ak approachs to that of the Q matrix of the QR decomposition of the matrix of 

eigenvector X. Define 

(6.215) 

and rij as the (i,j) element of R. From Ok ~ X, we find r1111 ~ :£1 when k is 

large. Since ~l is the eigenvector associated with the largest eigenvalue, we conclude 

that the first column of the matrix Ok approach the eigenvector associated with the 

largest eigenvalue of matrix A when k is large. If the matrix A is symmetric, which 

is often the case for many signal processing applications, the similar transformation 

Ak+1 = Or: AOk is also symmetric. Since Ak+1 approaches the upper triangular matrix 

by the QR algorithm, Ak+l approaches a diagonal matrix. That is 

(6.216) 

and 

(6.217) 

In this case, for large k, the columns of Ok become proportional to the columns of 

eigenvector in X. 

If A is real, then Ak will converge to a real block upper triangular matrix with 

1 x 1 and 2 x 2 main diagonal blocks. The complex conjugate pairs of eigenvectors 

of the 2 x 2 blocks can be solved easily using the quadratic formula. When A is not 

a square matrix, the singular values and vectors are of interest. For a m x n matrix 

B, where m > n, the SVD of B shows B = U1:VT , where U is a m x n matrix of 

123 
Figure 6.29: A circular multiplexer. 

orthogonal columns, V is a n x n unitary matrix, and 1: is a n x n diagonal matrix 

with diagonal elements called singular values given in descending order. For most 

situations where high condition numbers are not encountered, a simple symmetric 

n x n matrix C = BT B can be formed and the matrix V can be found by direct use 

of the QR algorithm. Similarly, U can be found by using D = BBT. 

6.4 Multi-phase Systolic Algorithms 

In this section, we introduced the multi-phase systolic algorithms to compute the QR 

algorithm. Two arrays, triangular and rectangular, can be used to compute the QR 

algorithm with some advantages and disadvantages for each. We shall show that our 

methods compute the Q matrix explicitly without requiring any global communication 

within the array. 3efore we consider the multi-phase algorithms, two communication 

switches are first discussed. A circular multiplexer is a device which takes its inputs 

and distributes them in different output positions as shown in Fig. 6.29. We use a 

skewed row to represent the circular multiplexier. A first in/first out (FIFO) buffer 

124 
I 

· IIIIIIII"-r 

III 

Output ....... ----: __________ :rt---

Figure 6.30: A first in/first out buffer. 

Input 

is a buffer which takes its input to output in a first come first out manner as shown in 

Fig. 6.30. Both devices are commonly used in computer and microprocessor systems 

for data arrangement [36]. The computation of a QR algorithm consists of two basic 

steps. Initially, set Al = A. 

(2) compute Ak+l = RkQk, stop if converge, otherwise go back to step (1). 

6.4.1 Multi-phase Triangular Systolic Array 

The QR Decomposition triarray proposed by Gentleman and Kung [26] is used in our 

approach. The R matrix is stored in the triarray after the computation. To compute 

the matrix Ak+1 in step (2), the Qk matrix has to be computed first. Let us call the 

computations in step (1) and step (2) an iteration. Several iterations are required for 

Ak to converge. For each iteration, we propose a three phase operation on a triarray 

as follows: 

125 
• Phase 1: QR decomposition for Ak 

Compute the QR decomposition of the matrix Ak = QkRk, with the upper 

triangular matrix Rk being stored in the triarray [29]. The data in Ak is inputted 

row by row skewed in time as shown in Fig. 6.31. 

• Phase 2: Computing the Qk matrix 

From the QR Decomposition, we have R;;T Ar = Qr. Let the ith column of 

matrices Ar and Qr be denoted by Qi and 9..i respectively. Then 

(6.218) 

Section 6.2 showed that R;;T ~ can be computed in a triarray same as the one 

used in Phase 1. Since the ith column of Ar is the ith row of Ak, then with Ak 

inputted row by row skewed in time as shown in Fig. 6.32, the operations of 

the processing cells are given in the second column of Table 6.5. The triarray 

computes the Qk matrix of Ak. The matrix Qk is then outputted row by row 

as shown in Fig. 6.32. In order to start Phase 3, the matrix Qk has to be in 

the form of Fig. 6.33. Observe that the output Qk of phase 2 shares the same 

snap-shot order as the desired arrangement of Qk in Phase 3 after a transpose 

operation. A circular multiplexer is used to distribute each column output of 

Qk into row input as indicated in Fig. 6.32. 

• Phase 3: Computing RkQk 

With the operations of the processing cell as shown in the third column of 

Table 6.5 and the Qk obtained in Phase 2, Fig. 6.33 shows the computation 

of Ak+l = RkQk in the triarray. Then the matrix Ak+1 comes out column 

by column from the right side of the triarray. Again, we observe that Ak+l 

126 
\1 

1 

Figure 6.31: Phase 1: The QR decomposition. 

127 
R 

Figure 6.32: Phase 2: Computing the Q matrix. 

128 
R 

Figure 6.33: Phase 3: Computing the matrix product RQ. 

129 
shares the same snap-shot order as the desired arrangement of Ak in Phase 1 

after a transpose operation. If not convergent, a new iteration is repeated by 

feeding back Ak+l into the triarray after using a circular multiplexer as shown 

in Fig. 6.33. Then Phase 1 operation begins as in Fig. 6.31. 

An attractive property of this multi-phase operation is that the feedback requirement of the matrices in different phases are identical. Thus, only a circular multiplexer is needed for each row outside the array. Observe that each column of the 

matrices inputted in all of the phases need n time steps to process and the next phase 

can be started at time n + 1. We find once the data outputted at the right hand 

side of the triarray, after passing through the circular multiplexer, it can be piped 

into the array for the next phase computation without suffering any delay. If we assume the multiplexer is ideal such the delay in the multiplexer can be ignored, it takes 

3n +(2n -1) = 5n -1 system clocks for one iteration. The (2n-l) term represents the 

initial time to feed the data into the array. Suppose the number of iteration required 

for convergence is 5, then the total number of system clocks needed is 35n + (2n - 1). 

Thus, the converge rate of this algorithm is of the order of 0((35 + 2)n). After the 

convergence of the Ak matrix, those values on the boundary cell are the eigenvalues 

of the A matrix. 

6.4.2 Multi-phase Rectangular Systolic Array 

The above method requires the use of the R-T operation in the computations. From 

a numerical stability point of view, we may want to consider an alternative that uses 

a square matrix for cumulative multiplication of the rotation parameters. Fig. 6.34 

shows a square matrix which is an extended version of the Gentleman-Kung's trian130 
,~ ~Ir ~~ ,..- ... '-

0 -- ~ ~ 

[A/II] 

L(' 1, " 1, 

.. .. .. . 

't ~, ,~ l' ..... .. .. 

--. 

4 • ~, • n 

" l' ..... ..... .. ...- -

~ I ~ 

~ I ~ 

" " 1t l' --- .... ~ 

-- ~ 

Figure 6.34: Multi-phase rectangular array for the QR iteration. 

131 
gular array with two delay elements (represented by black dots in Fig. 6.34) in the 

vertical communication links of the lower triangular part of the array. The processors 

in the lower triangular part are identical to the internal cells in the upper triangular 

part. Denote 

(6.219) 

as the parallel combination of matrix A and I, where A = [QllQ2,··· ,Qn] and I is 

a n x n identity matrix with ~ as its ith column. The square array takes the input 

[AIIIJ, while rotating A into an upper triangular matrix, I is used to cumulate the 

rotation parameters by 

Q[AII I] = [RIIQ]· (6.220) 

Noted that processors in the upper triangular part not only rotate the matrix A but 

also cumulatively multiply the rotation parameters with I. Thus, its work load is, 

in general, twice as that of Gentleman-Kung's internal cell. Processors in the lower 

part, on the other hand, only cumulate Q from the propagated rotation parameters. 

A two phases operation for QR iteration is proposed as follows: 

• Phase 1: QR decomposition 

Compute the QR decomposition of matrix Ak = QkRk; both Qk and Rk are 

obtained and stored. Then each row of Qk is piped out and fed back to the 

array through a FIFO buffer as shown in Fig. 6.34 . 

• Phase 2: Computing RkQk 

In this phase, the operation is identical to that of the Phase 3 in the triangular 

array. A circular multiplexer is used to transform Ak+l from row output into 

column input. Continue this iteration until converged. 

132 
Due to the delay elements at the lower triangular part of the array and the work 

load of each processor (except those in lower triangular part) is twice as that of 

triangular array, the time to obtain the ith row of the Qk matrix, t" is 

tj = max(2(n + 3i - 3), 2(2n + i - 2)), 

where 2( n + 3i - 3) is the time for the left-most cell of the ith row to obtain its Q 

element and 2(2n + i - 2) is the time for the right-most cell to finish. Obviously, 

when i ~ r nr 1, tj = 2( n + 3i - 3). Thus, the time required to obtain the whole Q 

matrix is tn = 8n - 6. By assuming that it takes time n to sequentially pipe out the 

Qk matrix, this algorithm takes (9n - 6) + n + (2n - 1) to complete an iteration in 

the worst cast. Again, denote the number of iteration as 5, this algorithm converge 

in the order of 0(5(10n - 6) + 2n -1) = 0((105 + 2)n). Of course, the performance 

can be improved by piping out each row of Q matrix when it is available instead of 

waiting for the whole Q matrix is available. \Vith this, the performance can reach to 

the order of O( (95 + 2)n). 

6.4.3 The Hessenberg Reduction 

In order to perform the QR algorithm efficiently in conventional Von Neumann type 

series computers, we usually transform the data matrix A into an upper Hessenberg 

matrix before the QR iteration is started. \Vith this transformation, the amount of 

work per iteration is reduced from 0(n3 ) to O(n2 ) [30]. However, this argument may 

not be true for parallel processing architectures. The reasons are two folds: 

1. Due to the hardware resource in a parallel processing architecture, the computations can be performed concurrently without hindering the processing time. 

133 
For example, the computation times of the two multi-phase arrays discussed 

above are of the order O(n). 

2. The data matrix is usually not of the Hessenberg form. The pre-processing of 

the data matrix to Hessenberg form may not be able to be incorporated with the 

following computations. That is, the pre-processing must be done separately. 

Both reasons may lead to the conclusion that the Hessenberg form is not of interest and practical for parallel processing of the QR algorithm unless the Hessenberg 

form can be obtained easily by using the same parallel processing architecture. Many 

papers prevented to answer this question by assuming the Hessenberg form (or the 

tridiagonal form) is already available. Fortunately, the Hessenberg form can be obtained easily as an additional part of the multi-phase operations. 

To obtain the Hessenberg form, we can choose an unitary similarity transformation 

U such that Al = UH AU is a upper Hessenberg matrix [30]. The transformation U 

can be obtained from sequences of Givens rotations. Denote Gi as the product of 

the Givens rotation matrices which zero out the proper positions of the ith column. 

Since the first ith rows will not be affected by Gi , the matrix Gi is of the form Gi = 

diag(Ii' G\), where Ii is an identity matrix of dimension i. Suppose the Hessenberg 

form through its first k - 1 columns have been obtained 

Bll B12 Bl3 

(G1 ·•· Gk_d H A(G1 ... Gk-d = B21 B22 B23 

o B32 B33 

134 

(6.221) 
where Bll and B33 are (k - 1) x (k - 1) and (n - k) x (n - k) matrices respectively. 

Then 

(Gl '" Gk)H A(Gl ... Gk) = B2l (6.222) 

o 

is a Hessenberg form through its first k columns. Thus 

(6.223) 

is an upper Hessenberg matrix and 

[ 

1 OT 1 U = G 1 ... G n-l = II ~ . (6.224) 

Denote 

A-- [g?]-UA- -- [1 llT 1· [ g? 1 , 

A II G R 

(6.225 ) 

where A GRand R is of the form of an upper triangular matrix without the 

lowest vertex element. Then Al = AU = UH AU. Obviously, this is similar to the 

computations in the QR iteration. To obtain G and R, let 

(6.226) 

The QRD of A is 

A A A [G llT 1 [ R 1 A = QR = !l 1 . 0,0" .. ,1 . (6.227) 

Now, we can use the multi-phase operations to obtain the upper Hessenberg form. 

We call this the Phase 0 operation. 

Phase 0: The Hessenberg Reduction. 

135 
(1) Phase 1 Operation: QRD of A. 

from A = Q fl, we obtain R in the triarray. 

(2) Phase 2 Operation: Computing G matrix. 

From Q = R-T AT, we obtain matrix G. 

(3) Phase 3 Operation: Computing the Hessenberg matrix AI' 

By forming 

and 

A= [~l 

u=[~ ~l' we obtain the Hessenberg matrix Al = AU. 

6.4.4 Computing the Eigenvectors 

To compute an eigenvector, a matrix multiplication systolic array can be incorporated 

with the multi-phase array such that those matrices QI,'" ,Qk are cumulated to form 

the Qk matrix. Noted that Ch = Qk-IQk and the matrix Qk-I is available at the start 

of the kth iteration. while the matrix Qk coming out at Phase 2 operation of the kth 

iteration. Then Qk is obtained by multiplying Qk-I and Qk as shown in Fig. 6.28. A 

system configuration for triangular array is shown in Fig. 6.35 and that for rectangular 

is shown in Fig. 6.36. As discussed in section 6.3, for a symmetric A matrix, when 

Ak converged, Qk yields the matrix of eigenvectors. For a non-symmetric A matrix, 

the first column of Qk yields the eigenvector associated with the largest eigenvalue. 

136 
(I 

Host Computer 

Bus 

Matrix 

Multiplier 

Figure 6.35: System configuration of the multi-phase triarray. 

"I 

Host Computer 

~ 

~ 

/ " , 

Bus 

'\. j~ 4 / 

... ... 4~ - , ., 

Matrix .... Rectangular - Multiplier Array 

Figure 6.36: System configuration of the multi-phase rectangular array. 

137 
Triangular Rectangular 

array array 

Computation O«3S+2)n) O«lOS+2)n) 

time worst case 

Numerical fair stable stability 

Number of n(n+1) rt plus 

cells 2 (n2 -n) d-element~ 

I/O ports 2n 3n 

Utilization 1 <1 

Communication 1 2 devices 

Table 6.7: Comparisons of the multi-phase triarray and rectangular array. 

6.5 Performance Efficiency 

6.5.1 Comparisons of the arrays 

Although there are three phases of operations, the arithematical operations in Phase 

2 and Phase 3 form a subset of the operations executed in Phase 1. Therefore we do 

not increase the cell complexity in the multi-phase arrays. The performance and characteristics of both triangular and rectangular arrays considered above are summarized 

in Table 6.7. The advantages of the triangular array are: it has less computational 

time as well as less number of cells, I/O ports, and communication devices. Furthermore, all the the processing cells are fully utilized. However, due to the computation 

of R-T in Phase 2 of the operations, it may be numerical unstable for certain highly 

138 
ill-conditioning data. For example, consider the matrix given by 

0.7601 -0.3967 0.6060 

-0.3967 1.7475 -0.1962 

0.6060 -0.1962 0.4924 

with eigenvalues {2.0, 1.0, 1O-12 }. If the triarray algorithm, which uses R-T to obtain 

QT, is used, the eigenvalues are obtained as {2.0, 1.0,3.6818 . 1O-6}. On the other 

hand, eigenvalues are obtained as {2.0, 1.0,9.9999· 1O-13} if the R-T is not explicitly 

computed. As a result, we have a complexity versus numerical stability tradeoff for 

the two multi-phase arrays. 

6.5.2 Rate of Convergence 

Similar to Luk in [72], by convergence we mean that the parameter of f(Ak) defined 

as 

(6.228) 

where N is the number of off-diagonal elements, has fallen below some prechosen 

tolerance value. As indicated in [72], it is difficult to monitor of f(Ak) in the parallel 

computation. Luk then proposed that the iteration be stopped after a sufficiently 

large number S of iterations. In the studies of Brent and Luk [12, 72], they found 

that S ~ 9 for random symmetric matrices of order n ~ 230 and S ~ 6 for n ~ 24. 

Therefore, they chosen S = 10 for n ~ 100 for Jacobi-like method. Similar to their 

approach, we apply the QR algorithm to random n xn symmetric matrices (aij), where 

the elements aij for 1 ~ i ~ j ~ n were uniformly and independently distributed in 

[-1,1]. The tolerance to meet the stopping condition is off(Ak) ~ 10- 1°. We can 

see from Fig. 6.37 that the number of iterations for a QR algorithm to converge is in 

139 
c 

o 

;: 

~ 

... 

C) 

::: 

100 

80 

60 

40 

20 ~ 

o 

o 

./ I ... I.rmion I 

/ 

V 

/ 

10 20 30 

n 

Figure 6.37: The number of iterations for a QR algorithm to converge versus the 

matrix size. 

the order of 10 for matrix size smaller than 20 x 20. Even though we can reduce the 

matrix to Hessenberg form for full dense matrix or tridiagonal form for symmetric 

matrix, and the QR iteration with origin shift can accelerated the convergent rate 

[32, 91, 98], the number of iteration is still in the order of 10. As an example, the 

4 x 4 tridiagonal matrix 

1 200 

234 0 

o 4 5 6 

o 0 6 7 

still requires eight iterations to converge when the most elegant symmetric QR algorithm, the Symmetrix QR Algorithm, was used [27, pp.424] This kind of property 

is not desirable for parallel processing implementation. It is known that Jacobi-like 

140 
method may require more flops as compared to the symmetric QR algorithm. However, due to parallel implementation, many rotations may take place at the same 

time. The computations involved in QR algorithm and Jacobi-like method are generally of the same complexity. From these discussions, the one which requires less 

number of iterations is more attractive from the parallel implementational point of 

view. Furthermore, the convergence rate of an QR iteration depends on the ratio of 

the eigenvalues. In our simulations, in more than 10% of the cases, the randomly 

generated symmetric matrices required significantly larger number of iterations to 

converge. As pointed out before, it is difficult to monitor the quantity off(Ak) to 

decide when the algorithm converges in the parallel computation. Since the convergence rate is highly dependent on the ratio of eigenvalues, there is no general rule 

to choose a sufficient number of iteration 5 like Jacobi-like method to insure convergence. This is also an undesirable intrinsic property of the QR algorithm for parallel 

implementation as compared to that of the Jacobi-like method. 

6.6 Efficient Fault-tolerance Schemes 

Reliable implementation is quite essential in parallel processing architectures. For 

a complex parallel processing system, a single fault from any part of the system 

can make the whole system useless. For various critical applications using spectral 

decomposition, highly-reliable computations are demanded. Fault-tolerance is therefore needed in many of these problems. A simple and cost effective fault-tolerant 

scheme is the checksum and weighted checksum proposed by Abraham et al [2]. This 

scheme is one of the typical examples of the algorithm-based fault-tolerance which has 

been applied to various signal processing and linear algebra operations [65]. Define 

141 
the checksum vector £.T = [1, 1, ... , 1], the column, row and full checksum matrices 

Ae, Ar and AJ of a square n-by-n matrix A are defined as 

If there is any fault occurring during the computation, the checksum criterion is not 

met and thus the fault is detected. The weighted checksum scheme can be further 

used to correct errors [47]. It has been suggested in [17] that checksum/weighted 

checksum scheme can be incorporated into the QR iteration for error detection. These 

properties as well as others are considered here for the multi-phase arrays. 

Since there are different operations in different phases, the inherent natures of 

the operations of each phase is thus different and should be examined for possible 

fault-tolerant implementation individually. The fault-tolerant schemes for each phase 

of the multi-phase triangular array are given as follows: 

• Phase 1: As pointed out in [17, 65], row checksum is invariant for the QR 

decomposition. It can be seen 

(6.229) 

This means that the QR decomposition of a row checksum matrix results in a 

row checksum upper triangulat matrix. Fig. 6.38 shows the implementation of 

this scheme. 

142 
( 

~3 

~2 a23 

~l a22 

a21 

R 

C 

S 

R 

C 

S 

Figure 6.38: Row checksum for Ar = QRr. 

143 
R 

Figure 6.39: Column checksum for R-T A~ = Q~. 

144 
• Phase 2: Due to the nature of computations in this phase, row checksum is no 

longer valid. Fortunately, column checksum is possible as given by 

(6.230) 

An implementation of this scheme is shown in Fig. 6.39 . 

• Phase 3: Although a row checksum upper triangular matrix Rr and a column 

checksum unitary matrix Qc are obtained in the above phases, unfortunately, 

RrQc does not yield any relevant use. By defining the trace operation as the 

sum of the diagonal elements in a square matrix, we obtain Tr[ABl = Tr[BA]' 

where A and B are square matrices. Therefore, 

n 

Tr[Ak+d = Tr[RkQkl = Tr[QkRkl = Tr[Akl = L Ai, (6.231) 

i=l 

where Ai is the eigenvalue of matrix Ak • This invariant property can be used 

to check the result of the Phase 3 operation. Once the trace of Ak is different 

from the trace obtained before, a fault is then detected during the Phase 3 

computation. 

For the rectangular array, the Phase 2 operation is the same as the Phase 3 

operation of the triangular array. For its Phase 1 operation, an interesting feature of 

this computation is given by 

(6.232) 

That is, a row checksum of the parallel combination of matrices A and I gives a 

row checksum of the upper triangular matrix Rr and a row checksum of the unitary 

matrix Qr. 

145 
Chapter 7 

Conclusions and Future Research 

This dissertation has addressed various aspects of the RLS problems and spectral 

decomposition algorithms based on QRD approach. 

In Chapter 2, we have shown that the Householder transformation can be implemented on a systolic array. By using a two-level pipelined implementation, the 

throughput of the SBHT RLS systolic array can be made as fast as that of the original Givens array in [78]. While the system latency is longer for the SBHT, it provides 

a better numerical stability than the Givens method. Clearly, the Givens array is a 

special case of the SBHT array with a block size of one. In general, the block size 

is an important variable. A larger block size results in a better numerical stability, 

while the system latency is increased. Many known properties of the Givens array 

are also applicable to the SBHT array. For example, the real-time algorithm-based 

fault-tolerant scheme proposed in [65] can also be easily incorporated into the SBHT 

RLS array. From the results described in Chapter 2, it appears that the Householder 

transformation method is useful in real-time high throughput applications of modern 

signal processing as well as in VLSI implementation. 

146 
In Chapter 3, special inherent natures of a QRD recursive LS systolic array are 

used to design a real-time, low-cost algorithm-based fault-tolerant system. By proper 

design of the artificial desired response, the residual output serves as a concurrent error detector. The disadvantages of the CSE scheme to detect fault are thus prevented. 

The proposed method requires the same complexity as that of the CSE scheme in 

[17] without hindering the throughput rate of the system. Two methods, FFL and 

CSE, are then introduced to locate the faulty row. For both methods, the recovery 

latency is achieved in O(p) time. However, the recovery latency of the CSE method is 

generally less than that of the FFL method in that parallel execution is possible and 

multiple ports can be accessed. Once the faulty row is determined, an order-degraded 

reconfiguration is performed so that the system is operating in an gracefully degraded 

manner. Any single fault occurring in the system will not cause a catastrophic degradation of the system and thus one critical requirement for many high performance 

and demanding VLSI signal processing tasks is met. 

This fault-tolerance system is general in the sense that it is independent of the 

implementation algorithms. It is applicable to Given rotation method as well as 

those methods such as Modified Gram-Schmidt, square-root free Given methods [19], 

and Householder transformations [66] etc. The LS problems with constraints such 

as MVDR beamforming is also applicable [80]. The residual method for concurrent 

error detection is robust in the sense that the probability of error detection given by a 

fault occurred equals one (for infinite-precision implementation). This scheme could 

have great impact on next generation of radar systems which may use VLSI array 

processors as their central signal processing units. 

In Chapter 4, We present an important observation that rotation parameters of the 

147 
QRD LS algorithm will eventually reach the quasi steady-state. With this observation, 

the dynamic range of each processing cell is then derived to insure correct operation 

of the algorithm. Also, we can prove the stability of the QRD LS algorithm under 

finite-precision implementation with this observation. The missing error detection 

and false alarm problems are considered. We present a design of the memory size 

which is overflow free without missing error detection and false alarm problems. 

In Chapter 5, the order degraded performance is derived and a geometric interpretation is given. A scheme to estimate the optimal residual under faulty situation 

is also presented. 

The multi-phase systolic algorithms proposed in Chapter 6 can be used efficiently 

to solve the eigenvalue and SVD problems based on the QR algorithm. In particular, the eigenvectors can be obtained without global communication within the array 

using the multi-phase operations. We showed that the QR algorithm can achieve a 

parallel implementation on a single architecture. Two systolic arrays, a triangular 

and a rectangular, are proposed for multi-phase implementation. Efficient algorithmbased fault-tolerance schemes can be incorporated with both arrays easily. Since the 

operations in each phase belong to the same types of computations, the cell complexity is thus not increased by multi-phase operations. There is a tradeoff between 

numerical stability and complexity for both arrays. Each iteration takes O(n) time 

unit while the time required for convergence is O(Sn), where S is the number of 

iterations. Unlike the Jacobi-like method, the convergence rate of the QR algorithm 

depends on the ratio of the eigenvalues. As a consequence, S may vary differently for 

the matrix of the same size, with or without origin shift to accelerate the convergence. 

Generally, S is in the order of 10 for the QR algorithm. While we have demonstrated 

148 
the advantage of the QR algorithm that can yield two multi-phase systolic algorithms 

implementable on single architectures without requiring global connections, the intrinsic convergence rate for the QR algorithm makes it less attractive as compared 

to the Jacobi-like method in parallel implementation. Depending on specific system 

and hardware requirements, one approach may be more desirable than the other. 

Of course, it is most meaningful to have two basic approaches to choose from for 

real-time VLSI signal processing based on spectral decomposition. 

There are various continuing problems extended from this dissertation that can 

be done in the future. For the SBHT RLS systolic algorithm, to pursue a possible 

CORDIC implementation would be an interesting problem. An algorithm-based faulttolerance has been proposed for the QRD RLS systolic algorithm. It would be a valid 

question to ask if there is any efficient fault-tolerant scheme for the transversal or 

lattice algorithms. Different from the residual estimation under faulty situration, it 

would be possible for us to ask if we can estimate the full-order optimal residual under 

the order degraded operation. Also, there are still some questions on re-convergence 

of a fault-tolerant QRD RLS systolic array. When there is a fault occurred in the 

system, if it is a transient fault, how long will this fault affect the system? In another 

words, how soon will the adaptive system recover itself from the faulty situation and 

converge. There are two questions. The first one is how long will it recover from the 

short time transient fault without any cure? Does the position of the faulty processor 

play an important role? The next question is when we switch to order degraded 

operation, do we want to reinitialize the system or to restart from the contaminated 

contents? 

Future research in the VLSI signal processing can be divided into the following 

149 
categories: 

1. Algorithms and architectures aspects 

To reconsider and to develop new parallel and pipelined algorithms, to design dedicated application-specific architectures for parallel algorithms (either 

in VLSI or DSP chip sets), and to find efficient fault-tolerant schemes for specialpurposed applications are the major issues in this category. 

2. Software Supported CAD Aspects 

To provide simulation tools for algorithms and architectures development and 

to support automatic design of complicated systems are the major concerns of 

this category. 

3. Application and Practical Implementation Aspects 

There is a large application spectrum. To name a few, the areas such as 

sonar / radar, communication, image processing, HDTV, coding, video, speech 

processing, and music are of potential interest to researcher in the VLSI signal 

proceSSIng. 

150 
Bibliography 

[1] A.H. Abdallah and Y.H. Hu, "Parallel VLSI computing array implementation for 

signal subspace updating algorithm", IEEE Trans. ASSP, Vol 37, pp.742-748, 1989. 

[2] J .A. Abraham et al., "Fault tolerance techniques for systolic array", IEEE Computer, Vol 20, pp.65-76, July 1987. 

[3] C.J. Anfinson and F.T. Luk, "A linear algebraic model of algorithm-based fault 

tolerance", Proc. IEEE Int'l Conf. Systolic Array, pp.483-493, May 1988. 

[4] C.J. Anfinson et al., " Algorithm-based fault-tolerant techniques for MVDR beamforming", Proc. IEEE ICASSP, pp.2417-2420, May, 1989. 

[5] M. Annaratone et al., "The Warp computer: Architecture, implementation, and 

performance", IEEE Trans. Computer, C-36, pp.1523-1538, Dec. 1987. 

[6] A. Avizienis, "Fault tolerant computing - An overview", IEEE Computer, Vol 4, 

pp.5, Jan. 1971. 

[7] P. Banerjee and J.A. Abraham, "Bounds on algorithm-based fault-tolerance in 

multiple processor systems", IEEE Trans. Computer, C-35, pp296-306, 1986. 

[8] M.G. Bellanger, "Computational complexity and accuracy issues in fast least 

squares algorithms for adaptive filtering", Proc. IEEE ISCAS, pp.2635-2639, Finland, 1988. 

[9] G. Bienvenue and H.F. Mermoz, "New principle of array processing in underwater 

passive listening", in VLSI and Modern Signal Processing, S.Y. Kung et al., 

Eds., Prentice-Hall, 1985. 

[10] A.W. Bojanczyk, R.P. Brent, and F.R. de Hoog, "Parallel QR decomposition of 

Toeplitz matrices" , Proc. SPIE Advanced Algorithms and Architectures for Signal 

Processing, pp.39-44, Aug. 1986. 

[11] A.W. Bojanczyk and A.O. Steinhardt, " Stabilized hyperbolic Householder transformations", IEEE Trans. ASSP, Vol 37, pp.1286-1288, Aug. 1989. 

151 
[12] R.P. Brent and F.T. Luk, "The solution of of singular-value singular-value and 

symmetric eigenvalue problems problems problems on multiprocessor array," SIAM 

J. Sci. Stat. Comput., Vol 6, pp. 69-84, Jan. 1985. 

[13] R.P. Brent and F.T. Luk, "A systolic array for the linear-time solution of Toeplitz 

systems of equations", Journal of VLSI and Computer Systems, 1, pp.1-22, 1985. 

[14] K. Bromley and J.M. Speiser, "Signal processing algorithm, architectures, and 

applications", Proc. SPIE, Vol. 431, pp. 2-6, 1983. 

[15] S.-W. Chan and C.-L. Wey, "The design of concurrent error diagnosable systolic 

arrays for band matrix multiplications", IEEE Trans. CAD, Vol 7, pp.21-37, Jan. 

1988. 

[16] M. Chean and J .A.B. Fortes, "A texonomy of reconfiguration techniques for 

fault-tolerant processor arrays", IEEE Computer, Vol. 23, pp.55-69, Jan. 1990. 

[17] C.-Y. Chen and J.A. Abraham, "Fault-tolerant systems for the computation of 

eigenvalues and singular values", Proc. SPIE, Vol 696, Advanced Algorithms and 

Architectures for Signal Processing, 1986. 

[18] C.-Y. Chen and J.A. Abraham, "Current error detection in VLSI processor arrays", Proc. SPIE, Vol 826, Advanced Algorithms and Architectures for Signal 

Processing II, pp. 205-214, 1987. 

[19] M.J. Chen, "On realizations and performances of least-squares estimation and 

Kalman filterring by systolic array", Ph.D. dissertation, Electrical Engineering 

Dept., UCLA, 1987. 

[20] S.1. Chou and C.M. Rader, "Algorithm-based error detection of a Cholesky factor 

updating systolic array using CORDIC processors", Proc. SPIE, Real-time Signal 

Processing XI, pp.104-111, Aug. 1988. 

[21] J.M. Cioffi, "Limited-precision effects in adaptive filtering", IEEE Trans. CAS, 

VOL CAS-34, pp.821-833, July, 1987. 

[22] J.M. Cioffi, "The fast QR adaptive filter", submitted to IEEE Trans. ASSP. 

[23] J.M. Cioffi, "The fast Householder filters RLS adaptive filter" , Proc. IEEE 

ICASSP, pp.161_9-1622, Albuquerque, April 1990. 

[24] P. Comon and Y. Robert, "A systolic array for computing BA- 1'" IEEE Trans. 

ASSP, ASSP-35, pp.717-723, June, 1987. 

[25] R.J. Cosentino, "Concurrent error corrections in systolic architecture", IEEE 

Trans. CAD, Vol 7, pp.l17-125, Jan. 1988. 

152 
[26] J.A.B. Fortes and C.S. Raghavendra, "Gracefully degradable processor arrays", 

IEEE Trans. Computer, Vol C-34, pp.l033-1044, Nov. 1985. 

[27] G.R. Gao and S.J. Thomas, "An optimal parallel Jacobi-like solution method 

for singular value decomposition", Proc. Int'l Conf. Parallel Processing, pp. 47-53, 

1988. 

[28] G.D. de Villiers, " A Gentleman-Kung architecture for finding the singular value 

of a matrix", Proc. Int'l Conf. Systolic Array, pp.545-554, Ireland, 1989. 

[29] W.M. Gentleman and H.T. Kung, "Matrix triangularization by systolic array," 

Proc. SPIE, Vol. 298, pp. 298-303, 1981. 

[30] G.H. Golub and C.F. Van Loan, Matrix Computation" 2nd edition, Johns 

Hopkins, 1989. 

[31] S. Haykin, "Radar array processing for angle of arrival estimation", in Array 

Signal Processing, pp.194-292, Haykin, Ed., Prentice-Hall, 1985. 

[32] S. Haykin, Adaptive Filter Theory, Prentice-Hall, 1986. 

[33] D.E. Heller and I.C.F. Ipsen, "Systolic networks for orthogonal equivalence transformations and their application", 1982 Conf. Advanced Research in VLSI, M.I.T. 

[34] D.E. Heller and I.C.F. Ipsen, "Systolic networks for orthogonal decomposition", 

SIAM J. Sci. Stat. Comput. Vol 4, pp.261-269, June, 1983. 

[35] B. Hochet, P. Quinton and Y. Robert, "Systolic Gaussian elimination over GF(p) 

with partial pivoting", IEEE Trans. Computer, Vol 38, pp.1321-1324, Sep. 1989. 

[36] W. Hoffmann and B.N. Parlett, "A new proof of global convergence for the 

tridiagonal QL algorithm", SIAM J. Numer. Anal. Vol 15, Oct. 1978. 

[37] S.H. Hosseini, "On fault-tolerant structure, distributed fault-diagnosis, reconfiguration, and recovery of the array processors", IEEE Trans. Computer, Vol 38, 

pp.932-942, July 1989. 

[38] K.-H. Huang and J.A. Abraham, "Algorithm-based fault-tolerance for matrix 

operations", IEEE Trans. Computer, Vol C-33, pp.518-528, June, 1984. 

[39] S.F. Hsieh, K.J. R. Liu, and K. Yao, "A fast and effective algorithm for sinusoidal 

frequency estimation", IEEE Int'l Symosium on Information Theory, San Diego, 

Jan. 1990. 

[40] S.F. Hsieh, K.J. R. Liu, and K. Yao, "Applications of truncated QR methods 

to sinusoidal frequency estimation", Proc. IEEE Int'l Conf. Acoustic, Speech, and 

Signal Processing (ICASSP), pp.2571-2574, Albuquerque, 1990. 

153 
[41] M.G.M. Hussain and M. Jaragh, "A triangular systolic array for the discrete-time 

deconvolution", IEEE Trans. CAS, Vol 36, pp.622-628, Apr. 1989. 

[42] K. Hwang and F.A. Briggs, Computer Architecture and Parallel Processing, McGraw-Hill, 1984. 

[43] I. Ipsen, "Singular value decomposition with systolic array," Proc. SPIE, Vol 

495, pp.13-21, 1984. 

(44) S.N. Jean, C.W. Chang and S.Y. Kung, "Graceful degradation schemes for 

static/dynamic wavefront array", Proc. Int'l Conf. Parallel Processing, pp.249-255, 

Aug. 1988. 

[45] 1. Johnsson, "A computational array for the QR-method," 1982 Conference on 

Advanced Research in VLSI, M.I.T., pp. 123-129. 

[46) J.-Y. Jou and J.A. Abraham, "Fault-tolerant matrix operation on multiple processor system using weigted checksums", Proc. SPIE Vol 495 Real Time Signal 

Processing VII, pp.94-101, 1984. 

[47] J.-Y Jou and J.A. Abraham, "Fault-tolerant matrix arithmetic and signal processing on highly concurrent computing structures", Proc. IEEE, Vol 74, pp.732-

741, May, 1986. 

[48) J.-Y. Jou and J.A. Abraham, "Fault-tolerant algorithms and architectures for 

real time signal processing" , Proc. Int'l Conf. Parallel Processing, pp.359-362, Aug. 

1988. 

[49] S. Kalson and K. Yao, "Systolic array processing for order and time recursive 

generalized least-squares estimation," Proc. SPIE, Vol. 564, Real Time Signal Processing VIII, pp. 28-38, 1985. 

[50] S.M. Kay, Modern Spectral Estimation, Prentice-Hall, 1988. 

[51] I. Koren and M.A. Breuer, "On area and yield considerations for fault-tolerant 

VLSI processor arrays", IEEE Trans. Computer, Vol C-33, pp.21-27, Jan. 1984. 

[52] K. Konstantinides and K. Yao, " Statistical analysis of effective singular values in 

matrix rank determination", IEEE Trans. ASSP, ASSP-36, pp.757-736, May 1988. 

[53] I. Koren and D.K. Pradham, "Yield and performance enhancement through redundancy in VLSI and WSI multiprocessor systems", IEEE Proc. Vol 74, pp.699-

711, May, 1986. 

[54] H.T. Kung, "Why systolic architectures?", IEEE Computer, Vol 15, pp.37, Jan. 

1982. 

154 
[55] H.T. Kung and M.S. Lam, "Wafer-scale integration and two-level pipe lined implementaiton of systolic array", J. Parallel Distrib. Comput., Vol 1 , pp.32-63, 1984. 

[56] S.Y. Kung et ai., VLSI and mordern signal processing, Prentice-Hall, 1985. 

[57] S.Y. Kung, C.W. Chang, and C.W. Jen, "Real-time reconfiguration for faulttolerant VLSI array processors", Proc. Real-Time Systems Symposium, pp.46-54, 

1986. 

[58] S.Y. Kung, VLSI Array Processors, Prentice Hall, 1988. 

[59] H. Lev-Ari and B. Friedlander, "On the systematic design of fault-tolerant processor arrays with application to digital filter", Proc. IEEE Workshop on VLSI 

Signal Processing, pp.483-494, Nov. 1988. 

[60] H. Leung and S. Haykin, "Stability of recursive QRD LS algorithms using finiteprecision systolic array implementation", IEEE Trans. ASSP, VOL37 pp.760-763, 

May 1989. 

[61] R.A. Lincoln and K. Yao, "Efficientsystolic Kalman filtering design by dependence graph mapping", VLSI Signal Processing III, pp.396-407, Nov. 1988. 

[62] F. Ling, D. Manolakis, and J.G. Proakis, "A recursive modified Gram-Schmidt 

algorithm for least-squares estimation", IEEE Trans. ASSP, Vol. ASSP-34, pp.829-

836, Aug. 1986. 

[63] KJ.R. Liu and K. Yao "On uniform one-chip VLSI design considerations for 

some discrete orthogonal transforms", Proc. IEEE International Conference on 

Acoustics, Speech, and Signal Processing (ICASSP) pp.2136-2139, New York, April 

1988. 

[64] K.J.R. Liu and K. Yao, "A systematic approach to bit recursive systolic array 

design" Proc. IEEE International Conference on Systolic Array, pp.685-694, San 

Diego, May 1988. 

[65] K.J.R. Liu and K. Yao, "Gracefully degradable real-time algorithm-based faulttolerant method for QR recursive least-squares systolic array", in Systolic Array 

Processors, Ed. McCanny, McWhirter, and Swartzlander, pp. 401-410, Prentice 

Hall (UK), 1989. 

[66] K.J.R. Liu, S.F. Hsieh, and K. Yao, "Recursive LS filtering using block Householder transformations", Proc. IEEE Int'l ConI. Acoustic, Speech, and Signal Processing, pp.1631-1634, Albuguerque, 1990. 

[67] K.J .R. Liu and K. Yao, "Spectral decomposition via systolic triarray based on 

QR iteration", Proc. IEEE Int'l Conf. Acoustic, Speech, and Signal Processing 

(ICASSP), pp.l017-1020, Albuquerque, 1990. 

155 
I 

[68] K.J.R. Liu, "Dynamic range for finite-precision QRD LS algorithm and its stability", Proc. Int'l Sym. Circuits and Systems (ISCAS), pp.3142-3145, New Orleans, 

May 1990. 

[69] K.J.R. Liu and K. Yao, "Multi-phase systolic architectures for spectral decom~ 

position" Int'l Conf. on Parallel Processing, August, 1990. 

[70] K.J.R. Liu, S.F. Hsieh, and K. Yao, "Two-level pipelined implementation of 

systolic block Householder transformations with application to RLS algorithm" 

Int'l Conf. on Application-Specific Array Processors, Princeton, Sep. 1990. 

[71] F.T. Luk, "A parallel method for computing the generalized singular value decomposition", J. Parallel and Distributed Computing 2, pp.250-260, 1985. 

[72] F.T. Luk, "A triangular processor array for computing singular value," Linear 

Algebra and Its Applications, Vol. 77, pp. 259-273, 1986. 

[73J F.T. Luk, "A rotation method for computing the QR-decomposition", SIAM J. 

Sci. Stat. Comput. Vol 7, pp.452-459, Apr. 1986. 

[74J F.T. Luk and H. Park, "Fault-tolerant matrix triangulization on systolic array", 

IEEE Trans. Computer, Vol 37, pp.1434-1438, Nov. 1988. 

[75J F. T. Luk and S. Qiao, "Analysis of a recursive least-squares signal-processing 

algorithm," SIAM J. Sci. Stat. Comput. Vol. 10, No.3, pp. 407-418, May 1989. 

[76J V.J. Mathews and Z. Xie, "Fixed-point error analysis of stochastic gradient adaptive lattice filters", IEEE Trans. ASSP, Vol 38, pp.70-80, Jan. 1990. 

[77J J.V. McCanny and J.G. McWhitter, "Some systolic array developments in the 

United Kingdom", IEEE Computer, Vol 20, pp.51-64, July 1987. 

[78] J.G. McWhirter, "Recursive least-squares minimization using a systolic array", 

Proc. SPIE, Vol 431, Real Time Signal Processing VI, pp.105-112, 1983. 

[79] J.G. McWhirter and T.J. Shepherd, "An efficient systolic array for MVDR beamforming," Proc. Int'l ConL Systolic Array, pp. 11-20, 1988. 

[80] J.G. McWhitter and T.J. Shepherd, "Systolic array processor for MVDR beamforming", lEE Proc. Vol 135, Pt. F, pp.75-80, 1989. ) 

[81] D.1. Moldovan, C.I. Wu and J.A.B. Fortes, "Mapping arbitrary large QR algorithm into a fixed size VLSI array", Proc. Int'l ConL on Parallel Processing, 

pp.365-373, 1984. 

[82] J.H. Moreno and T. Lang, "A multilevel pipelined processor for the singular 

value decomposition", Proc. SPIE Vol 698, Real Time Signal Processing IX, 1986. 

156 
[83] J.H. Moreno and T. Lang, "Comments on "A systolic array for computing 

BA-1'" IEEE Trans. ASSP, Vol 37, pp.1786-1789, Nov. 1989. 

[84] J.G. Nash, K.W. Przytula, and S. Hansen, "Systolic/Cellular processor for linear 

algebraic operations", VLSI Signal Processing II, pp.306-315, 1986. 

(85] N.L. Owsley, "Sonar array processing", in Array Signal Processing, pp.115-

193, Haykin, Ed., Prentice-Hall, 1985. 

[86] C.C. Paige, "Computing the generalized singular value decomposition", SIAM 

J. Sci. Stat. Comput. Vol 7, Oct. 1986. 

[87] B. N. Parlett, "Analysis of algorithms for reflections in bisectors," SIAM Review, 

Vol. 13, No.2, pp. 197-208, Apr. 1971. 

(88] S.K. Rao and T. Kailath, "Regular iterative algorithms and their implementation 

on processor arrays" ,Proc. of IEEE, 76, pp.259-269, March 1988. 

[89] C.P. Rialan and L.L. Scharf, "Cellular architectures for implementing projection 

operators", IEEE Trans. ASSP, ASSP-35, pp.1619-1627, Nov. 1987. 

[90] W. Robertson and W. Phillips, "A systolic MUSIC system for VLSI implementation", Proc. IEEE ICASSP, pp.2577-2580, 1989. 

[91] R. Roncella and R. Saletti, "A VLSI systolic adder for digital filtering of deltamodulated signal", IEEE Trans. ASSP, Vol 37, pp.749-754, 1989. 

[92] M. Sami and R. Stefanelli, "Reconfigurable architectures for VLSI processing 

arrays", Proc. IEEE, pp.712-722, May 1986. 

[93] RO. Schmidt, "A signal subspace approach to multiple emitter location and 

spectral estimation", Ph.D. dissertation, Stanford Univ. 1981. 

(94] R. Schreiber, "Systolic array for eigenvalue computation", Proc. SPIE Vol 341, 

Real Time Signal Processing V, 1982. 

[95] R Schreiber, "Systolic linear algebra machines in digital signal processing", in 

VLSI and Modern Signal Processing, pp.389-405, S.Y. Kung et al., Eds., 

Prentice-Hall, 1985. 

[96] R Schreiber, "Solving eigenvalue and singular value problems on an undersized 

systolic systolic array," SIAM J. Sci. Stat. Comput. Vol 7, pp.441, Apr. 1986. 

[97] R Schreiber, "Implementation of adaptive array algorithms", IEEE Trans. ASSP 

Vol ASSP-34, pp.l038-1045, Oct. 1986. 

157 
[98J K.G. Shin and T.-H. Lin, "Modeling and measurement of error propogation in 

a multimodule computing system", IEEE Trans. Computer, Vol 37, pp.1053-1066, 

Sep. 1988. 

[99J A. Steinhardt, "Householder transformation," IEEE ASSP Mag., 1988. 

[lOOJ C.W. Stewart, "Incorporating origin shifts into the QR algorithm for symmetric 

tridiagonal matrices", Comms. ACM, Vol 13, June, 1970. 

[lOlJ G.W. Stewart, Introduction to Matrix Computations, Academic Press, 

1973. 

[102J R.A. Thisted, Elements of statistical computing, Chapman and Hall, 1988. 

[103J N. Torralba and J.J. Navarro, "Size-independent systolic algorithms for QR 

iteration and Hessenberg reduction", Proc. Int'l Conf. Systolic Array, pp.166-175, 

1989. 

[104J N .-K. Tsao, "A note on implementing the Householder transformation," SIAM 

J. Numer. Anal., Vol. 12, No.1, pp. 53-58, Mar. 1975. 

[105] B. Van Veen, "Systolic processors for linearly constrained beamforming", IEEE 

Trans. ASSP, Vol 37, pp.600-604, Apr. 1989. 

[106J C.R. Ward, P.J. Hargrave and J.C. McWhirter, "A novel algorithm and architecture for adaptive digital beamforming", IEEE Trans. Antennas Propagat., Vol 

AP-34, pp.338-346, March, 1986. 

[107J J.H. Wilkinson, Algebraic Eigenvalue Problem, Oxford University Press, 

1965. 

[108J J.H. \Vilkinson, "Global convergence of tridiagonal QR algorithm with origin 

shift", Linear Algebra and Its Applications 1, pp.409-420, 1968. 

[109J J. H. Wilkinson, "Modern error analysis," SIAM Review, Vol. 13, No.4, pp. 

548-568, Oct. 1971. 

158