QR AND

GENERALIZATIONS

OF

THE SINGULAR

VALUE

AND

QR

DECOMPOSITIONS*

BART DE MOOR AND PAUL VAN DOOREN$

We dedicate thispaper to Gene Golub, a true source

_of

inspiration

for

ourwork, but alsoagenuinefriend, onthe occasion

_of

his60th birthday

Abstract. This paper discusses multimatrix generalizations of twowell-knownorthogonalrank factorizations ofa matrix: the generalized singular valuedecomposition and thegeneralized QR-(or URV-) decomposition. These generalizations can be obtained for any number ofmatrices of compatible dimensions. This paper discussesin detail the structure of these generalizations and theirmutual relations and givesaconstructiveprooffor the generalized QR-decompositions.

Key words, singularvalue decomposition, QR-factorization, URV-decomposition, complete orthogonaldecomposition

AM$(MOS) subject classifications. 15A09, 15A18, 15A21, 15A24,65F20

1. Introduction.

In

this paper,wepresent multimatrixgeneralizations ofsome

well-knownorthogonalrank factorizations.

We

showhow the ideaofa

QR-decompo-sition

(QRD),

a URV-decomposition

(URVD),

and a singular value decomposition

(SVD)

for one matrix can be generalized to any number of matrices. While

gener-alizations of the

SVD

for any number ofmatrices have been derived in

[9],

one of

the main contributionsofthispaperis the constructive derivationofageneralization for the

QRD

(or URVD)

for any number ofmatrices of compatible dimensions. The idea is to reduce the set ofmatrices

A1, A2,...,

Ak

to a simpler form using unitary

transformationsonly. Hereby, weavoidexplicitproducts andinverses ofthe matrices that are involved.

We

show that these generalized QR-decompositions

(GQRD)

can

beconsidered as apreliminaryreduction for anygeneralized singular value

decompo-sition

(GSVD).

The reason is thatthere is a certainone-to-one relationbetweenthe structure ofa

GQRD

and the "corresponding"

GSVD,

which is explained in detail below.

This paper is organized as follows.

In

2,

we provide a summary of orthogonal

rank

_{factorizations}

for one matrix.

We

briefly review the

SVD,

the

QRD,

and the

URVD

asspecial cases.

In

3,

wegiveasurvey o1existinggeneralizations of the

SVD

and

QRD

for two or three matrices.

In

4,

wesummarize the results on

GSVDs

for

anynumber ofmatricesofcompatible dimensions. Section

5,

which containsthemain newcontributionofthis paper, describes ageneralization ofthe

QRD

and the

URVD

for any number of matrices.

We

derive a constructive, inductive proofwhich shows

thata

GQRD

can be usedas apreliminaryreduction for acorresponding

GSVD.

In

6,

weanalyzeindetail thestructureof the

GQRDs

and

GSVDs

and show thatthere is aone-to-one relation between the two generalizations. Thisrelation is elaborated

in more detail in

7,

where weillustrate howa

GQRD

can be used as apreliminary

Step

inthederivationofacorresponding

GSVD.

Receivedbythe editorsApril 16, 1991; accepted forpublication (inrevisedform) October 18,

1991. This researchwas partially supported bytheBelgian Programon Interuniversity Attraction Poles andtheEuropeanCommunityResearchProgram ESPRIT, BRA3280.

Departmentof Electrical Engineering, Katholieke UniversiteitLeuven, B-3001Leuven,Belgium (demoor@esat.kuleuven.ac.be). Thisauthoris aResearchAssociate of the Belgian National Fund for ScientificResearch_(NFWO).

Coordinated Science Laboratory, University of Illinois, Urbana, Illinois 61801 (vdoorenOuicsl.csl.uiuc.edu).

While all results inthis paper are stated for complex matrices, they canbe spe-cializedto thereal casewithout much difficulty. Thiscanbe done inmuch the same way as with the

SVD

for complex and real matrices.

In

particular, it suffices to re-state most results using the term real orthonormalinsteadof unitaryand to replace

a superscript "."

_(which

denotes thecomplex conjugate transpose ofa

matrix)

by a superscript

"t" (which

is the transpose ofa

_matrix).

2. Orthogonal rankfactorizations.

Any

matrix

A

E

C

mn _{canbe factorized}

as

where

R

C

nxn _{is upper trapezoidal and}

_H

_{is a real n x n permutation matrix} that permutes the columns of

A

so that the first

ra

rank(A)

columns are linearly independent. The matrix

Q

C

mxm _{isunitary andcanbe partitionedas}

IFa m

ra

).

If we partition

R

accordingly as

R

(Rll

R12),

where

_Rll

C

rxra _{is upper}

triangularandnonsingular, weobtain

A=QI(RI Rle)H

whichis sometimes cMled the

QR-factorization

of

A.

Ifwe rewrite

₍₁₎

as

Q’A-(

R

weseethat

Q

isanorthogonal transformationthatcompressestherows of

A.

There-fore,

it is called a row compression.

A

similar construction exists, of course, for a

columncompression.

A

complete orthogonal

_{factorization}

ofanmx matrix

A

is any factoriation of the form

(2)

A

U

(

T

0

V*

where

T

is

_rax

_ra

square nonsingular nd

ra

rank(A).

One

particular case isthe

SVD,

which has become an important tool in the nalysis and numericM solution of numerous problems, especially since the development of numericMly robust

algo-rithms by Golub ndhiscoworkers

[15], [16],

[17].

The

SVD

is complete orthogonM

factorizationwhere thematrix

T

is diagonalwith positive diagonM elements:

A

UEV*.

Here U

C

mxm _and

_V

_C

nxn

are unitary and E

mxn

is oftheform

(71 0 0 0

0 _a2 0 0

0 0

_rra

0

0 0 0 0

In this paper, we use the convention that zero blocks may be "empty" matrices, i.e., certain block dimensions may be 0.

GENERALIZED SVD AND QR 995 The positive numbers 0-1

_

0"2

_

_

0-ra

>

0 are called the singular values of

A,

whilethe columns of

U

and

V

arethe

_left

and right singular vectors.

In

applications where m

>>

n,it isoftenagoodideato use the

QRD

ofthematrix as a preliminary step in the computation ofits

SVD.

The

SVD

of

A

isobtained via

the

SVD

ofits triangular factor as

A

QR-

Q(UrE.V?)

(QUr)ErV?.

This idea of combining the

QRD

and the

SVD

of the triangular matrix, in order to compute the

SVD

ofthe full matrix, is mentioned in

_[22,

p.

_119]

and more

fully

analyzedin

[3].

In

_[18]

the methodisreferredtoasR-bidiagonalization.

Its

flopcount is

(mn2+n3),

ascomparedto

(2mn2-2/3n3)

forabidiagonalization ofthefullmatrix.

Hence,

whenever m

_>

5/3n,

it is more advantageousto use the R-bidiagonalization algorithm.

There exist still other complete orthogonal factorizations ofthe form

(2)

where

only

T

is required to be triangular

(upper

or

_{lower) (see,}

e.g.,

[18]).

Sucha factoriza-tionwas called aURV-decompositionin

[27].

Here

where

U

E

C

m m

_V

_C

gularuppertriangular.

are unitarymatrices and

R

C

raxr _{is square}

nonsin-It

is well known that the QR-factorization of a singular matrix

A

and of its

transpose A* canbeused for findingthe image and kernelof

A

(URV-decompositions

actually give both at

_once).

In

this paper, we try to extend these ideas to several matrices.

Suppose

we have a sequence ofmatrices

Ai,

1,...,

k,

and we want to

knowthekernels

(or

null

spaces)

of eachpartialproduct

_{A1. A2...}

Aj.

We

could form these products andcompute QR-decompositions ofeach ofthem. That can, in

fact,

be avoided, as shown below.

Let

us take the "special" example

Ai

A,

1, 2, 3,

with 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0

It

iswellknown that the null spaces of

A

infactgive the Jordan structureof

A,

and thisstructure is alreadyobviousfrom the form of

A.

But

letusreconstructit from a sequence ofQR-decompositions

(in

factweneed hereRQ-decompositions of

A).

The first oneis, ofcourse, acolumn compression of

A1,

for whichwe usethe permutation

of columns 2 and 4

_(denoted

by thematrix

P24):

0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 1 0 0 0

The separation line here indicatesthat the first two columns of

P24

(i.e.,

el and

ea)

span thekernelof

A

A.

For

the kernelof

A

2

butapply theinverseof the orthogonal transform

_{P24 (which}

is again

P24)

tothe rows of

_A2

A:

0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0

Since

_A1A2

_{(AIP24)(P24A2),}

it is clear that the kernel of

_AIA2

is also the kernel

ofthe bottom part of

_P24A2.

The following column compression of

_P24A2

actually yields the kernel of both

_A2

and theproduct

_AA2.

Perform indeedthe orthogonal transformation

P24P35:

P24A2P24P35

0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0

We

see thatthekernelof

A2

comprisesthefirst twocolumns of

_P24P35

(i.e.,

el ande4

as

before)

and the kernel of

A1A2

comprises the firstfour columns of

P24P35,

i.e., e2, e4, and e5.

An

additional step ofthis procedure finally shows that the kernel of

the product

A1A2A3

(AIP24)(P24A2P24P35)(P35P24A3)

is thatofthe bottom part of thematrix 0 1 0 0 0 0 0 0 0 1

P35P24A3--

0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 which is a zero matrix.

Hence

the kernel of

A

3

is the whole space, as expected.

The interestingpart ofthis simple example"is the fact that we have not formed the

intermediate products to get their corresponding kernels. The case treated here of equal matrices

Ai

is asimple one

_(and

could besolvedusing theresults of

[19]),

but in the nextfewsections we show how thiscanalso be done forarbitrarysequencesof

matrices. The key idea is that at each step wedo anumber of QR-factorizations on

the blocks ofapartitioned matrix

(column

blocks in our

case).

This then induces a new partitioning onthe rows ofthis matrix,on the columns of thenext matrix, and

so on.

3. Generalizations for two or three matrices.

In

the last decade or so, several generalizations for the

SVD

have been derived. The motivation is basically the necessity to avoidthe explicit formationof products and matrix quotients inthe

computation of the

SVD

of products and quotients of matrices.

Let

A

and

B

be

nonsingular square matrices and assume that weneed the

SVD

of

AB-*

USV*.

2

It

is well known that the explicit calculation of

B

-1,

followed by the computation

of the product, may result in loss of numerical precision

(digit

cancellation),

even

before any factorization is

attempted!

The reason is the finite machine precision of Thenotation B-* refers tothe complexconjugatetransposeof theinverseofthematrixB.

any calculator.

Therefore,

it seems more appropriate to come up with an implicit combined factorizationof

A

and

B

separately, such as

A

UD1X

-1,

(3)

B=X-*D2V*,

where

U

and

V

are unitaryand

X

nonsingular. The matrices

_D1

and

D2

arereal but

"sparse" (quasi-diagonal),

and theproduct

_DID

isdiagonalwith positivediagonal elements. Thenwefind

AB-*

_UDIX-1XDtV

_U(DDt)v

*.

A

factorization as in

₍₃₎

is always possible for two square nonsingular matrices.

In

fact,

it is always possible for two matrices

AE

C

mn _and

_B

E

C

nv

(as

long as the numberof columns of

A

isthesame asthe numberofrowsof

B,

which werefertoas a

compatibility

condition).

In

general, thematrices

A

and

B

mayevenberank deficient. The combined factorization

(3)

is called the quotient singular value decomposition

(QSVD)

andwas firstsuggestedin

_[32]

andrefined in

_[23]

(it

wasoriginally calledthe

generalized

SVD,

but wehave suggested astandardizednomenclature in

_[6]).

A

similar idea might be exploited for the

SVD

ofthe product of two matrices

AB

USV*,

viatheso-called product singular value decomposition

(PSVD)

A

UD1

X-l,

(4)

B

XD2V*,

sothat

AB

U(DID2)V*,

which is an

SVD

of

AB.

The combined factorization

(4)

was proposed in

[13]

as aformalization ofideas in

[21].

In

thegeneral case, for two

compatiblematrices

A

and

B

(which

mayberank

deficient),

the

PSVD

of

(4)

always

existsandprovidesthe

SVD

of

AB

without theexplicit construction ofthe product. Similarly, if

A

and

B

are compatible, the

QSVD

always

exists.

However,

it does not always deliver the

SVD

of

ABt

when

B

is rank deficient

(Bt

is the pseudoinverse

of

_B).

Another generalization, this time for three matrices, is the restricted singular

value decomposition

_(RSVD).

It

was proposed

n

[35],

and numerous applications were reviewedin

_[7].

It

was soon found that allofthesegeneralized

SVDs

for twoor three matricesarespecialcasesofageneral

theorem,

presentedin

_[9].

Themainresult

is that there exist

GSVDs

for any number ofmatrices

AI,A.,...,Ak

of compatible dimensions. The general structure of these

GSVDs

wasfurther

analyzed

in

_{10].

The

dimensions ofthe blocks that occur in any

GSVD

can be expressed asranks of the matricesinvolved and ascertain products and concatenationsof these.

We

present a summary ofthe resultsbelow.

As

for generalizations of the

QRD,

it is mainly Paige

{25]

who pointed out the importance of generalized

QRDs

for two matrices as a basic conceptual and math-ematical tool. The motivation is that in some applications, we need the

QRD

ofa

product oftwo matrices

AB

where

A

mx

and

B

nxp.

For

general matrices

A

and

B

such a computation avoids forming the product explicitly, and transforms

A

and

B

separatelyto obtain the desired results. Paige

[25]

refers to such a factor-ization as a product

QR

_{factorization.}

Similarly, in some applications we need the

QR-factorization of

AB

-

where

B

is square andnonsingular.

A

general numerically robust algorithm would not compute theinverse of

B

northe product explicitly, but would transform

A

and

B

separately. Paige

[25]

proposed calling such a combined

decompositionoftwo matricesageneralized

QR

factorization,

following

[20].

We

pro-pose here toreservethename generalized

QRD

for the completeset ofgeneralizations of the QR-decompositions,which aredevelopedin thispaper.

We

also propose anovel

nomenclature,

as wedidfor the generalizations ofthe

SVD

in

[6].

Stoer

[28]

appears tobethefirst to havegivenareliable computationofthistype of generalized QR-factorization for two matrices

_(see

_[14]).

Computational methods for producingthetwo types of generalized

QR

factorizations for twomatrices, as de-scribed

above,

haveappearedregularlyin theliteratureas

_{(intermediate)}

stepsinthe

solutionofsomeproblems.

In

this paper, wederiveaconstructiveproofof

generaliza-tionsofthe

QRD

for any number ofmatrices.

As

wesee

below,

ourgeneralized

QRDs

can also beconsideredthe appropriate generalization ofthe

URVD

ofa matrix.

4. Generalizedsingularvaluedecompositions.

In

thissection, wepresenta

general theoremthatcanbeconsideredthe appropriategeneralization forany number

ofmatrices of the

SVD

ofonematrix.

It

contains the existing generalizations of the

SVD

for two matrices

_(i.e.,

the

PSVD

and the

QSVD)

and three matrices

(i.e.,

the

RSVD)

asspecial cases.

A

constructiveproofcan befoundin

[9].

THEOREM 4.1

(generalized

singular value decompositions for k

matrices).

Con-sideraset

_of

k matrices with compatible dimensions:

A1

(no

nl),A2 (nl

n2),...,

Ak-1

(nk-2

nk-1),Ak (nk-1

nk).

Then there exist

--Unitarymatrices

_UI

_(no

_no)

and

Vk (nk

nk).

--Matrices

_Dj,

j 1,

2,...,

(k

1)

of

the

_form

nj_ nj where

(6)

(7)

rj_ rj 2 rj 2 2 rj_ rj rj J rj J nj-1 rj-1 rj 2 3 J nj -rj [rj rj rj rj

I

0 0 0 0 0 0 0 0 0 0

I

0 0 0 0 0 0 0 0 0 0

I

0 0 0 0 0 0

I

0 0 0 0 0

ro

O,

A

matrix

Sk

of

the

_form

rk rk- rk

r

₂

r_l

rk

r_k J

rank(Aj).

rj

E

_rj i=1 k

r r r

r_k nk rk

o

o

o

o

0 0 0 0 0 0

S

0 0 0 0 0 0 0 0 0 0

S

0 0 0 0 0 nk-l rk-l r_k

_\

0 0

where

(8)

rk

E r

rank(A)

i--1

and the

_{rk rk}

matrices

S

are diagonalwith positivediagonal elements. Expressions

for

the integers rj aregiven in

₆

a

Nonsingularmatrices

Xj

(n

n)

and

Zi,

j

1,

2,...,

(k-

1)

where

_Z

is either

Zj

_X*

or

Z

_Xi

_(i.e.,

both choices are always

possible),

such that the

given matrices can be

_factorized

as

A1

U1D1X 1,

A2

ZID2X

1,

A3-

_Z2D3X

1,

Ai-

Zi-IDiX-1,

Ak

Zk- Sk

V;

Observe that the matrices

Dj

in

₍₅₎

and

Sk

in

₍₇₎

are generally not diagonal.

Theironlynonzero

blocks, however,

are diagonalblockmatrices.

We

propose tolabel

themasquasi-diagonalmatrices. Thematrices

Dj,

j 1,..., k-1arequasi-diagonal,

their only nonzero blocks being identity matrices. The matrix

Sk

is quasi-diagonal

andits nonzeroblocksarediagonalmatriceswith positivediagonal elements. Observe

thatwealwaystakethelastfactorinevery factorizationastheinverseofanonsingular matrix, which is only a matter ofconvention

(another

convention would result in a modified definition of the matrices

Zi). As

for the name of a certain

GSVD,

we propose to adopt the followingconvention

(see

also

[9]).

DEFINITION 4.2

_(the

nomenclaturefor

_GSVDs).

If k 1 in Theorem

4.1,

then

the correspondingfactorizationofthe matrix

A1

willbe called the

(ordinary)

singular value decomposition. If for a matrix pair

Ai, Ai+l,

1

<_

_<

k- 1 in Theorem 4.1,

we have

Zi

Xi, then the factorization ofthe pair is said to be of

P

type.

If,

on the other

hand,

for a matrix pair

Ai, Ai+I,

1

<_

_<

k- 1 in Theorem 4.1, we have

Zi

X-*,

then the factorizationofthe pair is said to be of

Q

type. The nameofa

GSVD

ofthe matrices

Ai,

1, 2,... k

>

1 as in Theorem 4.1, is then obtained by simplyenumerating the different factorizationtypes.

Let

usgive some examples.

Example. Considertwo matrices

_A

(no

n)

and

_A2

_(n

n2).

Then,

we have

two possible

GSVDs"

A1

P

type

Q

type

U1D1X

UID1X

X S:

V

X*

S: V;

The

P-type

factorization is called the

PSVD

(see

[8]

and references

therein),

while the

Q-type

factorizationis called the

QSVD.

Example.

Let

uswrite a

PQQP-SVD

for five matrices"

A1

U1DIX

1,

A2

_XID2Xf

1,

A3

Xf*

D3X

1,

A4

_{Xf D4X}

1,

A

X4SV.

We

also introducea notation using powers thatsymbolize acertainrepetitionof

aletter orofasequence of letters:

p3Q2-SVD

PPPQQ-SVD,

(PQ)2Q3(PPQ)2-SVD

PQPQQQQPPQPPQ-SVD.

Despite thefact thatthereare 2k-1 _{different sequencesof letters}

_P

_and

_Q

_atlevel

k

>

1, not all ofthesesequences correspondto different

GSVDs.

The reasonfor this is

that,

for instance, the

QP-SVD

of

(A

1,

A

2,

A

3)

canbe obtainedfrom the

PQ-SVD

of

_((A3)

*,

(A2)

*,

(A1)*).

Similarly, the

P2(Qp)3-SVD

of

(A1,...

,A

9)

isessentiallythe same as the

(PQ)3p2-SVD

of

((A9)*,..., (A1)*).

The number ofdifferent

factoriza-(2k-1

2k/2

_(2k-1

(k

tions for k matrices is, in

fact,

/ for k even and /:2

-1)/2)

for k odd.

A

possible way to visualize Theorem 4.1 is to builda tree with all different

fac-torizationsfor

1,

2, 3, etc matrices asfollows:

O

P

Q

p2

_pQ

_Q2

p3

pQ

pQp

pQ

QpQ

Q3

5. Generalized

URVDs.

In

thissection, wederiveageneralization for several matrices, ofthe

URVD

ofone matrix.

We

proceed inseveral stages. First, we show how k matrices can be reduced to block

tria.ngular

matrices using unitary

transfor-mations only.

Next,

weshow how the block triangular factors can betriangularized furthertotriangular factors.

THEOREM 5.1 Given k complex matrices

AI

_(no

_hi),

_A2

_(hi

n2),

Ak

(nk-1

nk),

there always exist unitarymatrices

Qo,

Q1,..., Qk

such that

where

_Ti

is a block lower triangular orblock upper triangularmatrix

_(both

cases are

always

possible)

with thefollowing structures:

--Lower

block triangular

(denoted

by a superscript

l)

(9)

T/t

2 i--1

r+l

ri ri r ri ri_

Ti,1

0 0 0 0 2 ri_

_*

Ti,2

0 0 0

r_

T,

0

--Upper

block triangular

(denoted

by asuperscript

u)"

(10)

T/u

2 i-1

_r+l

ri ri r ri ri-1

*

0 2

Ti,2

*

*

0 ri_ 0 0

T,

0 ri_

where

Ti,j

j

1,...,

are

_full

columnrankmatricesand each represents anonzero

block. The block dimensions coincide with those

of

Theorem4.1.

In

particular,

r

no,

r

+1

_nullity(Ai)

ni ri, and

E

j

rank(A)

r r j--1 ri_ hi-1. j=l

Our

proof ofTheorem 5.1 is inductive:

We

obtain the required factorization of

Ai

fromthat of

Ai-1.

Proof.

The inductionis initialized for 1 as follows. First, take thecasewhere

T1

is to be lower block triangular.

Use

a unitary column compression matrix

Q1

to

reduce thematrix

A

to

where and

T=AIQ=r

(TI,1

0),

r

rank(T1,1)

rank(A),

rl2

nullity(A1)

n

rl,

r

o.

The casewhere

T1

is requiredto beupper block triangularis similar:

T*

AIQ1

=r

_(T1,1

0

).

Observethatwe havetaken

Qo

_Ino.

Now,

we canstartourinduction.

Assume

that wehave the requiredfactorization

for the first i-1 matrices

=QoAQ,

where the matrices

Tj,j

1,...

,i- 1 have the block structure as in Theorem 5.1.

We

now want to find a unitary matrix Qi such that

Ti

Q_IAQi

is either lower or upper block triangular. First, consider the case where

_Ti

is to be lower block triangular. The matrix

_Q-IA

canbe partitioned accordingtothe dimensionsof the block columns of

T_I

as

(11)

ni ri_ 2 r

,

Q_

n

ri_

It

isalwayspossible to constructaunitarymatrix

Q

to compress the columns of each

ofthe block rowsto the left as

(12)

1-1 /and where the subblocks

Ti,j

are of full column

rank,

denoted by

r,

r

+1

_nullity(A).

_Hereto,

_{we first compress the first block row of}

₍₁₁₎

_{to the left} with unitarycolumntransformationsappliedtothe fullmatrix. Thenweproceedwith thesecond blockrowin thedeflatedmatrix

(i.e.,

withoutmodifyingthe previousblock

column).

By

repeating this procedure times, wefind therequired form

(12).

Obviously,

1_<

l-l,

i.

(13)

r ri_l,

The construction of

Ti

when it is required to be upper block triangular is similar.

Construct

a unitary matrix

Q

that compresses the columns of the block rows of

Q-IAi

tothe right. The onlydifference isthatwe nowstartfrom thebottomto find that

(14)

r+l

2 ri /’i ri 0

Til

*

ri-1 2 _{0 0}

T2

*

T

Q_IAiQi

ri-1 0 0

_Ti,i

l’i-We

can nowapplyanadditional

(block)

column permutation to theright ofthematrix

T/

so asto findthe matrix of

(10).

This completes theproof.

We

now demonstrate that the matrices

T,j

can always be further reduced to

inthe casewhen

_Ti

is lower block triangular.

Here,

R,j

is alower triangularmatrix

Similarly, we can always reduce

Ti,j

to

inthecasewhere

_Ti

isupper blocktriangular.

Here

R

.

is anuppertriangularmatrix 3

In

order todemonstrate this, weneed the following result.

LEMMA

5.2.

Let

PI,...,

Pk

bek givencomplexmatriceswhere

_Pi

hasdimensions pi-1 pi, pi-1

>_

pi and

rank(Pi)

pi. Then there always exist unitary matrices

Q0, Q1,...,Qk

such that

i--1

where

_Ri

is either

_of

the

_form

(15)

Pi

Ri-Pi-l-pi

)pi

(

0

R

with

R

a lowertriangular matrix, or

(16)

R-

pi

p-pi Pi

with

_R

upper triangular.

For

every

1,...,k,

both choices,

(15)

and

_(16),

are

alwayspossible.

Proof.

Again, the proofis by induction, but now for decreasingindex i.

For

the

initialization, start with k and obtain a QR-decomposition of

Pk

with either an upperor a lower triangular factorasrequired. This defines the unitary matrix

Qk-.

We

takeQk

Ipk. Hence,

wefind

--Lower

triangular:

(0)

--Upper

triangular:

Pk

lk-lRk

Qk-1

We

cannowstart

the.

inductionfor k 1, k

2,...,

1.

Therefore,

assumethat wehave the requiredfactorizations for the matrices

Pk, Pk-1,..-,

_Pi+"

R

k

(*k_ Pk

k

R

k-1

l*k_

2

Pk-

l

k-1,

Then,

if

_Ri

istobe lower triangular, obtain aQR-decomposition ofthe product PiQi

as

(o)

PiQi Qi-lRi Qi-1

R

sothat

R

Q-IPQ"

If

Ri

is required tobe uppertriangular, obtain a QR-decompositionas

PiQi Qi-IRi Qi-1

sothat

R

Q-I

PQi.

Thiscompletes the construction.

We

now repeatedly apply

Lemma

5.2 on the full rank blocks inthe matrices

_Ti

in

₍₉₎

and

(10).

First, weapply

Lemma

5.2 to the sequence ofk subblocks

Next,

weapplyit tothe sequence of the k 1 subblocks

In

general,we apply

Lemma

5.2 ktimes to the k sequences of subblocks

Tj,j,Tj+I,j,...,Tk,j

for j

1,...,k.

In

applying

Lemma

5.2tothe jth of these sequences,we canfindasequenceofunitary matrices

QJ] QJ]

[J]_Udk-j+l and matrices

Ri,j

such that

Ti,j

o[J]

i-j

Ri,j’i-j+l

i=j,...,k, where

or

We

nowdefine the unitarymatrices

(i

for 0,...,

k,

which areblock diagonalwith blocks

with

Qk+l]

I.

Next

wedefine

Qi_

TiOi,

i=O,...,k.

Then,

it canbe verified

that

for the lower triangular case we obtain

2

r+l

r r r ri_l ,1 0 0 0 2

Ri,2

0 0 I’i_ ri_l

_*

_*

Ri,i

0

and for the upper triangularcase wefind that

(18)

i--1

r+l

r r r ri_l ,1

_*

_*

0 9.

R,9.

0 ri_ 0 0

_R,

0 ri_

Ifwe nowcombine

_(9)-(17)

and

(10)-(18),

weobtaina combined factorizationofthe form

Hence,

wehave proved the followingtheorem.

THEOREM

5.3

(generalized

URVDs).

Givenk complexmatrices

_A1

_(no

x

nl),

A2

(n,

x

n9),

Ak

(nk-

x

nk),

there always exist unitary matrices

Qo,Q,...,Qk

such that

where is a lower triangular or upper triangular matrix

_(both

cases are always

possible)

with the following structures:

--Lower

triangular

(denoted

by a superscript

l)"

where

2

_r+l

r r r

11

Ri,1

0 0 0

2

,li

ri-1

*

Ri,2

0 0

and

_R,j

is a square nonsingular lower triangularmatrix.

--Upper

triangular

(denoted

by a superscript

u)"

i--1

_ri+l

r r r ri_l ,1

*

*

0 2

Ri,2

0 ri_ 0 0

Ri,i

0 ri_ where

and

R

is a square nonsingular upper triangular matrix. The block dimensions co-incide with those

_of

Theorem 4.1.

As

for the nomenclature of these generalized

URVDs,

we propose the following

definition.

DEFINITION 5.4

_{(nomenclature}

for generalized

URV).

The nameofageneralized

URVD

of kmatricesofcompatible dimensionsisgeneratedbyenumerating the letters

L

(for lower)

and

U

(for

upper),

according to the lower or upper triangularity ofthe matrices

T,

1,...,

k inthe decompositionofTheorem 5.3.

For

k matrices, thereare 2k different sequences with two letters.

For

instance, for k 3, there are eight generalized

URVD

(LLL, LLU,

LUL, LUU, ULL,

ULU,

UUL,

VVV).

Remarks. The decompositions in Theorems 5.1 and 5.3 both use column and

row compressions of a matrix as a cornerstone for the rank determination of the individualblocks.

As

already pointed out in

2,

the rank determination can bedone via anordinary

SVD

(OSVD),

buta moreeconomicalmethodusesthe

QRD

asinitial

step, since typically the matrices involved here have many more columns than rows or vice versa.

A

further alternative would be to replace the

OSVD

ofthe triangular

matrix resulting from the initial

QRD

by a rank-revealing

QRD.

Since the time of

the initial paper drawing attention to this

[5],

much progress has been made in this

area, and we only want to stress here that such alternatives can only benefit our decomposition.

The overall complexity of this

GQRD

is easilyseen to be comparableto that of performing two

QRDs

ofeach matrix

_A

involved.

For

each

_A

we indeed applythe

left transformation

Q-I

derivedfrom theprevious matrix andthenapplya "special"

compression

Q

of the resulting matrix while respecting its block structure. Both

steps have a complexity comparable to a

QRD

ofa matrix of the same dimensions.

For

parallel machines we can check that the "block" algorithms

[18]

for one-sided

orthogonal transformations such as the

QRD

can also be applied to the present de-composition, and thattheywillyield satisfactory speedups. The mainreasonfor this is that the two-sided orthogonal transforms applied to each

Ai

are done separately, and hencethey canessentially be considered one-sidedfor parallellization purposes.

6.

On

the structure of the

GSVD

and the

GQRD.

In

this section, we first point out how for each

GSVD

there are twogeneralized

URVDs,

and weclarify

the correspondence between the two types of generalized decompositions.

Next,

we J _{in Theorems 4.1 and 5.1} give asummary ofexpressions for the block dimensions r

in terms of the ranks ofthe matrices

A

1,...,

Ak

and concatenations and products

thereof. These expressions werederived in

[10].

Recall the nomenclature for the generalized

URVDs

(Definition

5.4)

and the

GSVDs

(Definition

4.2).

The relationshipbetweenthese two definitionsis asfollows.

A

pairofidentical

letters,

i.e.,

L-L

or

U-U,

thatoccursinthe factorization of

Ai,

_Ai+l

corresponds to a

P-type

factorizationofthe pair.

A

pair ofalternating

letters,

i.e.,

L-U

or

U-L,

thatoccursin the factorizationof

A,

_Ai+l

corresponds to a

Q-type

fac-torization of the pair.

As

an example, for a

PQP-SVD

of four matrices, there are twopossible corresponding generalized

URVDs,

namelyan LLUL-decompositionand

aUULU-decomposition.

As

withthe

GSVD,

we canalso introduce theconventionto use powers of

_(a

sequence

of)

letters.

For

instance, for a

p3Q2-SVD,

there are two

GURVs,

namely, an

LnUL-URV

and a

U4LU-URV.

j 4

_Let

_{usfirst considerthe}

We

nowderive expressions forthe block dimensions

_r.

case ofa

GSVD

that consists only of

P-type

factorizations.

Denote

the rank ofthe

product ofthe matrices

A,A+I,...,

Aj

with

_<

j by

ri(i+l)...(j-1)j

rank(AA+l

Aj_IAj).

THEOREM

6.1

_(on

the structure of

Pk-I-SVD, Lk-URV,

and

Uk-URV).

Con-sider any

of

the

_{factorizations}

above

_for

the matrices

A1,

A2,...,

Ak.

Then,

the block dimensions

_rji

that appearin Theorems

4.1,

5.1, and5.3 are given by:

(19)

(20)

--r( rj 1)(2)...(j), rj ri(i+l)...(j r(i-1)(i)...(j), with rj ri

if

j.

Next,

consider the case of a

GSVD

that only consists of

Q-type

factorizations.

Denote

the rank ofthe block bidiagonal matrix

(21)

A

0 0 0 0 0

Ai*+l A+e

0 0 0 0 0

_Ai*+3

_Ai+4

0 0 0 0 0 0

_{A_}

₃

Aj_2

0 0 0

_{A_}

_Aj

(by

rili+ll...lj_llj).

THEOREM

6.2

_(on

the structure of

Qk-I-SVD,

(LU)k/2-URV

(k even),

(UL)

k/2-URV

(k even),

(UL)(k-1)/2U-URV

(k odd),

and

(LU)(k-1)/2L-URV

(k odd)).

Con-sider any

_of

the above

_{factorizations}

_for

the matrices

A1,

Ae,..., Ak.

Then,

Ifj-iiseven,

2

it}+2

r}+4

j--2 j

ril...Ij ril...Ij-1 -t-

(r

+

rj

+...

+

rj)

-t- -t- -t-

+

rj -t-rj;

Ifj-i is

odd,

4 _{Recall that the subscript refersto the ith matrix, while the superscript j refers to the jth}

For

thegeneral case, we needamixture of the twoprecedingnotations for block bidiagonal matrices, the blocks ofwhich canbe products of matrices, such as

(A

A.

..A3_1

0 0

(A3

Ai4-1)

Ai4

A-I

0

0

_A,

Aj

where 1

_{_< io}

<

il

<

i2

<

i3

<

<

it

_<

j

_<

k. Theirrankis denoted by

For

instance, the rank of thematrix

A2A3 0 0

A?

A5A6A7

0

0

_(A8A9)*

_AlO

is _{represented by r(2)(3)141(5)(6)(7)1(8)(9)l(10).}

THEOREM

6.3

(on

the structure of a

GSVD

and a

_GURV).

r(io)(io+l)...(i-l)li...(i.-1)l...li...j can be derived as

follows:

i=1,2,.,

l+l"

1. Calculate the following

₊

1 integers

_By,

The rank 2

r}O

sj rj

-[-rj

2

ro+l

rO+2

rl

sj

+

+...+

2. Depending on even orodd there are two cases:

--I

even:

t

odd:

Observe that Theorems 6.1 and 6.2 are special cases of Theorem 6.3. While interms ofdifferences Theorem 6.1 provides a_{direct expression of the dimensions rj}

of ranks of products, Theorems 6.2 and 6.3 do soonly implicitly. This is illustrated inthe following examples.

Example.

Let

us determinethe blockdimensions of the quasi-diagonalmatrix$4 in a

QPP-SVD

ofthe matrices

A1, A2,

A3,

A4

(which

arealso the block dimensions of

an

LUUU-

or a

ULLL-decomposition). From

Theorem 6.2 wefindthat

and

r(1)1(2)(3)(4 rl

+

s4

2,

sothat

r42

r11(2)(3)(4 rl.

Finally, since r4

r

+

r

+

r43

+

r4

4,

wefindthat

r

rl

+

r(2)(3)(4) r11(2)(3)(4).

Observe that this last relation can be interpreted geometrically as the dimension of

theintersectionbetweenthe row spacesof

A1

and

A2A3A4:

r41

dim

spanrow(A1)

+

dim

spanrow(A2A3A4

dim

spanro

A2A3A3

Example. Consider the determination of

r, r,

r53,

r,

r

in a

PQ3-SVD

of five matrices

A1,A2, A3, A4,

A5

withTheorem

6.3,

which coincides withthe structure of

a

UULUL-URV

or an

LLULU-URV

_(see

Table

1).

TABLE r415 r2J314[5 r(1)(2)J3J4J5

8

r

+r +r

+r

4

8

r55

8 _=-r5

+

_r5

+

_r5 85 r5 1+}_2 8 r r 82 r3 3_{__r}4 85

8

r55

3

s

--r

+

r₅

These relationscanbe used to set up aset ofequationsfortheunknowns

r,

r,

r53,

r5

,

r55,

usingTheorem 6.3 as 1 1 1 1 1

r

_r5 0 0 0 0 1

r52

_r415

r4 1 1 1 0 1

_r53

_{r3j4j5 r3]}4 0 0 1 0 1

_r54

_r2131415

_r2J314

0 1 1 0 1

r55

_{r(1)(2)J3J4J5 r(1)(2)J3J4} the solutionofwhich is

1010 BARTDE MOOR AND PAULVANDOOREN

r

r(1)(2)131415 r(1)(2)1314

r2131415

--

r21314,

r

r2131415

r21314

r415

+

ra,

r54

r5

_r31415

-

r314,

r55

r415

r4.

7.

A

further block diagonalization of the

GQRD. In

thissection, wenote

that afurther block diagonalization ofa

GQRD

can be interpretedas apreliminary step towards the corresponding

GSVD. We

proceedin two stages. First, we observe that each upper or lower triangular matrix in the generalized

URVD

of Theorem 5.3 can be block diagonalized.

Next,

we show how these block diagonalizations can

be propagated backward through the

GQRD.

The first step is the factorization of the upper and lower triangular matrices

_Ti

of Theorem 5.3 into an upper or lower triangularmatrixandablock diagonalmatrix.

For

lowertriangularmatrices

i

/,

we can obtain afactorization ofthe form

where

2 i--1

ri_l ri_l ri_ ri_l

I

0 0 0 ri_ 2 ri_ *

I

0 0 ?i-1 * * *

I

2

r+l

/’i /’i ri_l ,1 0 0 0 2

Ri

2 0 0 0 0

_Ri,i

0

i--Since the diagonal blocks

Ri,j

are of full column

rank,

such a factorization is

al-ways possible.

In

a similar way, for upper iangular matrices

i

/u,

we find a factorizationof the form

with

Ui

anuppertriangular blockmatrixwithidentitymatricesonthe blockdiagonal:

2 i--1

ri--1 i--1 ri--1 i--1

I

ri_ 2 ₀

_I

ri_ 0 0 0

I

ri_ 2

r+l

ri ?i l’i ri-1 ,1 0 0 0

Ri,2

0 0

r_

0 0

Ri,i

0

Now

suppose that wehave done thisfor all matrices

i,

1,..., kin a

GQRD

of Theorem 5.3.

We

show how we can propagate a further block diagonalization backward throughthe

GQRD,

inaway thatiscompletelyconsistent withthe

corre-sponding

GSVD

of Theorem4.1.

To

simplify the notation, wesimply replace

_Ti

by

Ti

and

Di

by

D

inthe following.

First, assumethat

Tk

is lowerblock triangular.

It

thenfollows from theprevious sectionthatwe can factorize

_Tk

as

Tk

LD.

Dependingonwhether

T-I

isupper orlower triangular, wehave twocases"

.T-I

T_

lower triangular.

In

this

case,

theproduct

Tk_IL

k islower

trian-gularas

well,

andwe can obtain asimilardecomposition

Tk

t-lLa

La_Da_,

where

La_

is againlower triangularand

D_

has the same diagonal blocks

Ri,

as

Tk-1

Tkl

upper triangular.

In

this case, the product

_T%L*

is upper

triangular, andwe can obtain afactorization

-Uk-lDk

T_I

L;*

where

_U-I

is upper triangularand

_Dk-1

hasthe samediagonal blocks

R,y

as

_Tkl.

It

is easily verified that when

_Tk

is upper triangular, similar conclusions can be obtained.

In

generM,

let

_T

be lower triangular and assume thatit is factorized as

T

L

D Z.

Assume

that

T_

is lowertriangular. Then

_T-I

canbe fctored as

T_

L_D_L 1.

If

T_

is uppertriangular, it can be factoredas

TiL

U_

D_

1Li,

where

U-I

is uppertriangular. Thecases with

_T

uppertriangulararesimilar. Table

2 summarizesall possibilities.

Ti Lowertriangular

Ti LiDZ

T

Uppertriangular

T

Ui

D Z

TABLE 2

Ti- Lowertriangular Ti-1 Uppertriangular Ti-1 Lowertriangular Ti_ Uppertriangular

Ti-1 Li-1 Di-1

_U

Ti-i Vi-1Di_

_U

-1

Example.

Let

us applythis result to a sequenceof four matrices

A1,A2,A3,A4

with compatibledimensions. Iftherequired sequence is

ULUU,

then

A1

QoTQ*

Qo(UD1L:)Q

(QoU1)DI(QIL2)*

A:

QIT2Q

QI(L:DU)Q

(QIL2)D2(Q2Ua)*,

A3

Q2TQ

Q2(U3D3ui)Q

(Q2U3)D3(Q3U4)

-1

Note

that

U1

_Ino.

Thisfollows immediately fromtheblockstructureof

_U

for 1.

Observe that the relationships between the common factors in theleft-hand side of these expressions conform with the requirements for a

QQP-SVD.

Only

the middle

factors

D,

1, 2, 3,

4 are notquasi-diagonal.

8. Conclusions.

In

this paper, aconstructive proofwas givenofamultimatrix

generalization of the concept of rank factorization. The connection ofthis new

de-composition with the analogous

GSVD

was also shown. The blockstructure of both generalizations and the ranks of theindividualdiagonal blocksinboth decompositions

were indeed shown to be identical.

As

is shown in a forthcoming paper, the spaces

spanned by certainblock columns of the orthogonaltransformation matrices

Q

are,

in

fact,

identicaltothoseofthe

GSVD.

The difference lies only inaparticular choice

ofbasisvectorsfor thesespaces. Theconsequencesof theseconnections arestillunder

investigation.

We

mentionthe following results here:

Updating the above decompositiontoyield the

GSVD

requiresnonorthogonal

transformation. Theseupdatingtransformations canbe chosenblock triangularwith

diagonal block sizescompatiblewiththe indexsets derived in Theorem4.1.

A

modified orthogonal decomposition can be defined where the compound

matrix is not triangularized but diagonalized. This new factorization is a variant of

the above decomposition where now a special coordinate system is chosen for each

of theindividual orthogonal transformations

Q.

The resultis an orthogonal

decom-position ofthe type ofTheorem 5.3 where now the generalized singular values can be extracted from the diagonal elements of some triangular blocks. The orthogo-hal updating needed to obtain this new decomposition can be done with techniques

described in

_[2].

A

geometric interpretation can be givenofthe bases obtainedfromthe trans-formation matrices

Q

in Theorem 5.1.

As

particular examples of these spaces we retrieve the followingwell-knownconcepts.

(a)

For

the case

_A

_(A-I)

the

GQRD

infactreconstructs the nested null spaces

ofthe matrices

_(A

_aI)

,

which reveal theJordanstructureofthe matrix

A

at theeigenvalue a

(see

also theexample in

2).

(b)

For

the cases

_A2

_(A-aB)

and

_A2+1

B

the decomposition reconstructs the

nested null spacesof thesequences

[B-.I(A

aB)]

and

[(A

aB)B-1]

,

which

revealthe Kronecker structureof thepencil

AB- A

atthe generalized eigenvalue

c

_(see

_[30]

and

[31]).

(c)

For

the cases

A

D

and

_A

C.

A

-.

B,i 1,..., the decomposition reconstructs theinvertibility subspaces ofthe discrete time

system

Xk+l

Axk

+

Buk,

Yk

Cxk

+

Duk.

These are in fact also the spaces constructed by the structure algorithm of Sil-verman

[29],

andthey play arole inseveralkey problems of geometrical systems theory

[34].

Other applications of

GSVDs

have beendescribed in

_{[7], [8], [11], [13],}

and

[35],

while applications of the generalized QR-decompositions are described in

_[25]

and Acknowledgment.

Gene

Golub’s hospitality inthesummer of 1989 liesat the

generatingthe

MATLAB

codeforthe

GSVD

and the

GQRD

basedonthe constructive

proofs of the

GSVD

in

[9]

and ofthe

GQRD

ofTheorem 5.3ofthis paper.

REFERENCES

[1] A. BJRCK,Solvinglinearleastsquaresproblems byGram-Schmidt orthogonalization,BIT,7

(1967),pp. 1-21.

[2] A. BOJANCZYK, M. EWERBRING, F. LUK, AND P. VAN DOOREN, An accurateproductSVD

algorithm, in Proc. Second Internat. Sympos. on SVD and Signal Processing, Kingston,

RI, pp.217-228, June1990; SignalProc.,to appear.

[3] T. CHAN, Animproved algorithmforcomputing the singular valuedecomposition,ACMTrans.

Math. Software, 8(1982),pp. 72-83.

[4]

,

Algorithm581: Animproved algorithmforcomputingthe singular value decomposition,

ACM Trans. Math. Software,8(1982), pp. 84-88.

[5]

,

RankrevealingQRfactorizations,LinearAlgebra Appl.,88/89 (1987),pp. 67-82. [6] B. DE MOOR AND G. H. GOLUB, Generalized singular value decompositions: A proposal

for astandardizednomenclature,ESAT-SISTA Report 1989-10,Departmentof Electrical

Engineering,Katholieke Universiteit Leuven,Leuven,Belgium, April 1989.

[7]

,

The restricted singular value decomposition: Properties and applications, SIAM J.

MatrixAnal.Appl., 12 (1991), pp. 401-425.

[8] B. DE MOOR, On the structure and geometry ofthe product singular value decomposition, LinearAlgebra Appl., 168(1992), pp. 95-136.

[9] B. DE MOOR AND H. ZHA, A tree _of generalizations _ofthe singular value decomposition,

ESAT-SISTA Report 1990-11, Departmentof Electrical Engineering, Katholieke

Univer-siteitLeuven, Leuven, Belgium, June1990; LinearAlgebra Appl.,to appear.

[10] B. DE MOOR, On the structure _ofgeneralized singular value decompositions, ESAT-SISTA Report 1990-12, Department of Electrical Engineering, Katholieke Universiteit Leuven, Leuven,Belgium, June1990.

[11]

,

Generalizations oftheOSVD: Structure, properties, applications, in Proc. Second

In-ternat. Workshop on SVDandSignal Processing, University of Rhode Island, Kingston, RI, June1990,pp. 209-216; SignalProc.,to appear.

[12]

,

A history _ofthe singular value decomposition, ESAT-SISTA Report, Department of Electrical Engineering, Katholieke UniversiteitLeuven, Leuven,Belgium, February1991. [13] K. V.FERNANDOANDS. J. HAMMARLING, Aproductinducedsingular valuedecomposition_for

two matrices and balanced realisation, in Linear Algebrain Signal Systems andControl, B. N.Datta, C. R. Johnson,M. A. Kaashoek, R. Plemmons,andE. Sontag, eds., Society for Industrial andApplied Mathematics,Philadelphia,PA, 1988,pp. 128-140.

[14] P. E. (ILL, W. MURRAY, AND M.

n.

WRIGHT, Practical Optimization, Academic Press,

London, 1981.

[15] G H. GOLUBANDW.KAHAN, Calculating the singular values and pseudo-inverseofamatrix,

SIAMJ. Numer. Anal.,2(1965), pp. 205-224.

[16] P A. BUSINGERANDG.H. GOLUB,Algorithm358: Singular value decomposition_ofacomplex

matrix,Comm.Assoc. Comp. Mach., 12 (1969),pp. 564-565.

[17] G H. GOLUB AND C. REINSCH, Singular value decomposition and least squares solutions,

Numer. Math., 14 (1970),pp. 403-420.

[18] ( H. GOLUBAND C. VAN LOAN, MatrixComputations, Second Edition,The JohnsHopkins

UniversityPress,Baltimore, MD,1989.

[19] G GOLUBANDJ. WILKINSON,Ill-conditioned eigensystems andthecomputationoftheJordan

canonicalform, SIAMRev., 18(1976),pp. 578-619.

[20] S HAMMARLING, The numerical solution _ofthe general Gauss-Markov linearmodel, NAG

Tech. ReportTR2/85,NumericalAlgorithms GroupLimited,Oxford, 1985.

[21] M. T. HEATH, A. J. LAUB, C. C. PAIGE, AND R. C. WARD, Computing the singularvalue decomposition of a product _of two matrices, SIAM J. Sci. Statist. Comput., 7 (1986), pp. 1147-1159.

[22] R. J. HANSONANDC. L.LAWSON,Extensions andapplications_ofthe Householder algorithms forsolving linear least squaresproblems, Math. Comp.,23 (1969),pp. 787-812.

[23] C.C. PAIGEANDM.A. SAUNDERS, Toward8ageneralized singular value decomposition,SIAM J. Numer.Anal.,18 (1981), pp. 398-405.

[24] C. C. PAIGE, Computing thegeneralized singular value decomposition, SIAM J. Sci. Statist.

Comput., 7(1986),pp. 1126-1146.

CoxandS. Hammarling,eds., OxfordUniversityPress, Oxford,U.K., 1990, pp. 73-91. [26] G. W. STEWART, A method_for computing the generalized singular value decomposition, in

Proc. MatrixPencils, PiteHavsbad, 1982, B.KagstrSmandA. Ruhe,eds., Lecture Notes inMathematics973, Springer-Verlag, Berlin, NewYork, 1983.

[27]

,

An updating algorithm_forsubspace tracking, UMIACS-TR-90-86, CS-TR2494,

Com-puterScience Tech.ReportSeries, University ofMaryland, CollegePark,MD,July 1990. [28] J. STOER, Onthe numerical solution _ofconstrainedleast-squares problems, SIAM J. Numer.

Anal.,8(1971),pp. 382-411.

[29] L. M. SILVERMAN, Discrete Riccati Equations: AlternativeAlgorithms, AsymptoticProperties andSystem TheoryInterpretations, inControland DynamicalSystems, AcademicPress, New York,1976.

[30] P. VANDOOREN, Thecomputation_ofKronecker’s canonical_{form of}asingular pencil, Linear

AlgebraAppl., 27(1979), pp. 103-140.

[31]

,

Reducing subspaces: Definitions, properties and algorithms, in Proc. MatrixPencils,

Lecture NotesinMathematics973, Springer-Verlag, Berlin,NewYork, 1983,pp. 58-73. [32] C. F. VAN LOAN, Generalizingthe singular value decomposition, SIAM J. Numer. Anal., 13

(1976),pp. 76-83.

[33] J. H. WILKINSON, The Algebraic Eigenvalue Problem, The Clarendon Press, Oxford, U.K., 1965.

[34] M. WONHAM, Linear Multivariable Control, a Geometric Approach, Springer-Verlag, New York, 1974.

[35] H. ZHA, The restrictedSVD_formatrixtriplets and rank determination_ofmatrices, Scientific

Report89-2, Konrad-ZuseZentrumfiirInformationstechnik, Berlin,Germany, 1989;SIAM

J.Matrix Anal. Appl., 12 (1991),pp. 172-194.

[36]

,

The implicit QR decomposition and its applications, ESAT-SISTA Report 1989-25,