2007

LIBRARY
Michigan State
University

This is to certify that the
dissertation entitled

A Rank-Revealing Method for Low Rank Matrices with
Updating, Downdating, and Applications

presented by

Tsung-Lin Lee

has been accepted towards fulfillment
of the requirements for the

Ph D degree in Dggartment of Mathematics

 

 

Pg» ’” / / k ( )

 

Major Professor’s Signature
7//7/07

Date

MSU is an afi‘innative-acfion, equal-opportunity employer

 

SlL-.--u-n—.-.-a-n-l--0_n-I-I-I-b--u-o-a-l-n-o---u. _ 4

.u-I-o-l-u--.--D-I-l-I-l-I—l-o-I-I-.-a-o-l-u-I-I-b-I-I-o-I-

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE

DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/07 p:/CIRC/DateDue indd-p.1

 

A Rank-Revealing Method for Low Rank Matrices with
Updating, Downdating, and Applications

By

Tsung—Lin Lee

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Mathematics

2007

ABSTRACT

A Rank-Revealing Method for Low Rank Matrices with Updating,
Downdating, and Applications

By

Tsung-Lin Lee

As one of the basic problems in matrix computation, rank-revealing has a wide
variety of applications in scientific computing. Although singular value decomposition
is the standard rank-revealing method, it is costly in both computing time and storage
when the the rank or the nullity is low, and it is inefficient in updating and downdating
when rows and columns are inserted or deleted. Consequently, alternative methods
are in demand in those situations. Following up on a recent rank-revealing algorithm
by Li and Zeng in the low nullity case, we present a new rank-revealing algorithm
for low rank matrices with efficient and reliable updating/downdating capabilities. A
comprehensive computing experiment shows the new method is accurate, robust, and

substantially faster than existing rank-revealing algorithms.

To my parents.

iii

ACKNOWLEDGMENTS

With great pleasure, I would like to express my sincere gratitude to Professor Tien—
Yien Li, my dissertation advisor, for his guidance and support during my graduate
study at Michigan State University. As his student, I have been privileged to see and
learn many mathematical insights from him.

I would like to thank Professor Zhonggang Zeng for his inspiring discussions and
suggestions.

I would also like to thank Professor Gang Bao, Professor Chichia Chiu, Professor

Guowei Wei, and Professor Di Liu for their time and concern.

iv

TABLE OF CONTENTS

LIST OF FIGURES vii
LIST OF TABLES viii
Introduction 1
1 A Rank-Revealing Method 4
1.1 Approximate ranks ............................ 4

1.2 The convergence theory .......................... 6

1.3 Computing the approxi-range and the approxi-rowspace ........ 10

1.4 The overall algorithm ........................... 16
1.5 Numerical experiments .......................... 18
1.5.1 Type 1: Low approxi-rank with median noise level ....... 18

1.5.2 Type 2: Increasing approxi-rank, fixed size and median noise level 19

1.5.3 Type 3: Decreasing gap, fixed size and median noise level . . . 20

1.5.4 Type 4: High noise level with increasing size .......... 20

1.5.5 Type 5: Near zero noise level, low approxi-rank, large gap . . 21

2 Updating and Downdating Problems 23
2.1 The USV-plus decomposition ...................... 23
2.2 Updating and downdating ........................ 25
2.2.1 Updating ............................. 26

2.2.2 Downdating ............................ 27

2.3 Numerical results on updating and downdating ............ 32
2.3.1 Row-updating with increasing approxi-ranks .......... 33

2.3.2 Row-updating without changing approxi-ranks ......... 34

2.3.3 Row-updating with a small approxi-rank gap .......... 34

2.3.4 Row-downdating without changing approxi-ranks ....... 35

2.3.5 Row-downdating with decreasing approxi-ranks ........ 35

3 Applications 37

3.1 Information retrieval: latent semantic indexing ............. 37
3.2 Image processing: saving storage of photographs ............ 43
BIBLIOGRAPHY 46

vi

LIST OF FIGURES

1.1 Algorithm larank ................................................................................................... 22
2.1 Algorithm lrowup .................................................................................................. 28
2.2 Algorithm lrowdown .............................................................................................. 32
3.1 The photograph formation process: lattice partition and assignment ................... 44
3.2 Rank 21 approximation of Figure 3.1 .................................................................... 44
3.3 The original image and three lower rank approximation images .......................... 45

Vii

1.1

1.2

1.3

1.4

1.5

2.1

2.2

2.3

2.4

2.5

3.1

3.2

3.3

3.4

3.5

3.6

LIST OF TABLES

Results for Type 1 matrices .................................................................................... 19
Results for Type 2 matrices .................................................................................... 19
Results for Type 3 matrices ................................................................................... 20
Results for Type 4 matrices .................................................................................... 21
Results for Type 5 matrices .................................................................................... 21
Results for row-updating with increasing approxi-ranks ...................................... 33
Results for row-updating without changing approxi-ranks ................................... 34
Results for row-updating with a small approxi-rank gap ...................................... 35
Results for row-downdating without changing approxi-ranks ............................... 36
Results for row-downdating with decreasing approxi-ranks ................................ 36
Term-by-document matrix .................................................................................... 39
Comparisons for a 3000 x 1400 term-by-document matrix from CRAN .............. 41
The retrieval results for three lower rank databases and the raw databases ......... 42
Results for updating databases ............................................................................... 42
Results for downdating databases ......................................................................... 43
Comparison results for a 480 x 640 fingerprint image matrix .............................. 43

viii

Introduction

Rank-revealing appears frequently in scientific computing such as signal process-
ing [13, 28, 36], information retrieval [3, 12, 38] and numerical polynomial algebra
[10, 37]. While the singular value decomposition (SVD) is undoubtedly the most re—
liable method for determining the numerical rank of a matrix, it has drawbacks in
certain situations. In particular, it is expensive when matrix size becomes large but
either the rank or the nullity is low, and it is difficult to update or downdate when
rows or columns are inserted or deleted. Alternative methods have been proposed for
those situations, such as rank-revealing QR decomposition (RRQR) [5, 6, 7], rank-
revealing two-sided orthogonal decompositions (UTV, or URV/ULV) [16, 34, 35],
and rank-revealing LU decomposition (RRLU) [21, 26, 29]. In low-nullity cases, a
new rank—revealing algorithm has been developed by Li and Zeng [24]. We follow up
with a new rank-revealing algorithm for low rank matrices.

For a given m x n matrix A, our method determines the approximate rank
of A by calculating the approximate range of A. We briefly outline the method as fol-
lows. With m 2 n and rank(A) = k, let 01 Z 02 2 ~ - 2 0k > 0 be non-zero singular
values of A along with u,- and v,- being the corresponding unit left singular vectors
and unit right singular vectors associated with 02-, respectively, for i = 1, - - - , k. Un-

less otherwise mentioned, we shall always use “singular vector” to represent the right

singular vector. Since uJT-A = a-vT- forl S j g k,
A = alulv] + - - - + akukv; = uluIA + ---+ ukugA.

Clearly,

A1 := A—alulv] = A—ulu]A

has the same set of singular values along with associated singular vectors as those
of A except the largest singular value 01 of A is replaced with 0 as a singular value
of A1, and the second largest singular value 02 of A becomes the largest singular

value of A1. Thus, the rank of A1 becomes k — 1. Similarly, if k _>_ 2, matrix
A2 := A1 — uQuEA = A — uluIA - ugugA

has the same set of singular values of A except 01 and 02 are replaced with 0 and
the rank of A2 is reduced to k — 2. For the problem of finding the approximate rank,
namely the number of singular values larger than a prescribed threshold 6 > 0 (see
Definition 1 in §1.1), we begin by finding a unit vector fil in the approzi-mnge, namely,
the subspace spanned by the left singular vectors of A associated with singular values
larger than the threshold 0 > 0. This task can be accomplished efficiently by applying
the power iteration on AAT with a proper stopping criterion. We must emphasize here
that we do not require fil to be any of the left singular vectors of A. It can be shown

that (Theorem 4 in §1.3) the rank of the matrix
5471 := A — {1113'er

is one less than the rank of A. Similarly, a unit vector fig in the approxi-range of A1 is

also in the approxi-range of A, and the rank of the matrix
~ ,_ ~ ~T ~ ~T
A2 .— A — ululA — u2u2A

is reduced by another one, making it two below the rank of A. Moreover, Til and fig are
orthogonal since fil is in the left kernel of A2. This process continues recursively and
terminates with the approximate rank identified as well as an orthonormal basis
obtained for the approxi—range.

Our method has been implemented as a Matlab package LOWRANK. Comprehen-
sive numerical results of our code comparing with UTV Tools [16] and Matlab SVD
function are exhibited in §1.5. For low rank matrices, our code is consistently faster
than UTV Tools and the full SVD by a large margin.

Moreover, row /column updating and downdating in our method are quite simple
and straightforward. In §2.3, numerical results on both cases are presented to compare
our method with UTV Tools in this respect. While UTV Tools may sometimes return
incorrect ranks or inaccurate ranges, our method reliably yields accurate results on
all the matrices we tested.

Practical applications of our algorithms on information retrieval and image pro-

cessing are presented in §3.

CHAPTER 1

A Rank-Revealing Method

1 .1 Approximate ranks

The terms rank, nullity, and kernel are used in the exact sense as in common linear
algebra textbooks. The approximate rank, also known as the numerical rank, has a
specific meaning as in Definition 1 below. We use specific terms approxi-mnk, approxi-
range and approxi-rowspace for them respectively, but rank(A) is still the exact rank
of matrix A.

We shall denote matrices by upper case letters such as A , B, R, and column
vectors by lower case boldface letters like u, v and y. Notation ()T stands for the
transpose of a matrix or a vector, (-)-L represents the orthogonal complement of a
subspace, and [H] denotes the 2-norm of a matrix or a vector. The symbol 0,-(M) will

denote the i-th largest singular value of matrix M.

Definition 1. [18] For a given threshold 0 > 0, a matrix A E Rm” is of approxi-
rank k within 0, denoted by rankg (A) = k, if k is the smallest rank of all matrices

within a 2-n01‘m distance 6 of A. Namely,

rank9(A) min {rank(B)] = k. (1.1.1)

_ IlA—Bllso

Hereby, we also say the appromi-nullity of A within 6 is n — k.

The exact rank may be regarded as a special case of the approxi—rank
since rank(A) = rank9(A) for any matrix A within sufficiently small 6.

The minimum in (1.1.1) is attainable [18, 27]: Let the singular value decomposition

of A be

A = UZVT = olulv] + 02u2vg + ---+ onunv; (1.1.2)
where U = [u1, - -- ,um] and V = [v1, - -- ,vn] are orthogonal matrices along with
diagonal matrix 2 = diag{ol, - - - ,on} formed by singular values 01, - - ~ ,on satisfy-
ing

012"'20'k>920],3+1Z"'20n20- (1.1.3)
Then [[A — Ak[] = 0k +1 is the minimum 2-norm distance
from A to any rank k matrix for Ak = U ZkVT with diagonal ma-

trix Z = diag 0 ,-~- ,0 .,0,--- ,0 . Moreover, rank(A ) = rank (A) = k. In
It 1 k k 6

other words, for
5 = inf {a | rank#(A) = k} and 3 = sup {n | rankn(A) = k},

we have 6 = 0k + 1 and 3 = 0k' We call the ratio '7 = 3 / E the approxi-rank gap.

The fundamental subspaces range 7?.(A), kernel [C(A), left kernel [C(AT) and row
space ’R(AT) associated with matrix A can be naturally generalized in the approxi-
mate sense. In terms of the SVD of A in (1.1.2) with singular values satisfying (1.1.3),
the approximate subspaces of A along with their notations are listed as follows.

0 R9(A = span{u1, - -- ,uk}: The approxi-range of A within 6.

0 1C6(A) = span{vk + 1, - - - ,vn}: The approxi-kemel of A within 6.

o R9(AT) = span{v1, - -- ,vk}: The approxi-rowspace of A within 6.

(

0 (C6 AT) = span{uk + 1, - -- ,um}: The approzi-leftkernel of A within 6.

1.2 The convergence theory

For a given m x n matrix A and a rank threshold 6 > 0, we can assume m 2

n, rankg (A) = k and the SVD of A is given in (1.1.2) with

UlZo-Zak=8>025=0k+12...

IV

For a vector 2 79 0 and a subspace W in R’, the distance between z and W, denoted

by dist(z, W), is defined as the distance between supspaces span{z} and W. Namely,

_ z — WWTz
dist(z,W) = H ”z” N

 

if columns of matrix W form an orthonormal basis for subspace W. We
say a sequence {zj}§o= 1 of nonzero vectors converges into sub—
space W if . lim dist(zj,W) = 0.
.7 "* 00
Our strategy of revealing the approxi—rank of A is to construct an orthonormal
basis for the approxi-range R9(A). For this purpose, we use the power iteration
on AAT as follows: For a randomly generated unit vector yo 6 Rm, define se-

quences {xj} and {yj} as

x,- = ATyj_1/IIATy,_1n. y,- = ij/Hiju for j=1.2.--- (1.22)

that converge into the approxi-rowspace R6(AT) and the approxi—range 729(A), re—

spectively, at convergence rates given in the following proposition.

Proposition 2. For 6 > 0 and A E Rmxn with SVD in (1.1.2)
and singular values satisfying (1.2.1), let yo 6 Rm such that yo 9!

R9(A)i, then the sequences {xj} and {yj} generated by iteration (1.2.2) converge

into 726(AT) and R9(A) respectively at linear rate
dist(x ,TR9(A )) g ¢2j_1 and dist(yj,R9(A)) 3 a2], j=1,2,-~ (1.2.3)

with

 

V l .
451 E (g) dist-(y0,’Rg(A)) 2 forl 6 {1,2, - - - }.
\/1 — dist (yo, R9(A))

Proof. Write yo = 61111 +- - -+cmum. Without loss of generality, we assume m 2

ii. Let matrix 0: [ul, - -- ,uk]. From AAT = UEZTUT and for some n E R,

Y1 77 (clogul + - ~ ~ + enogun) (1.2.4)

 

0i 012: ”6+1 02
=61 751114- +ckjuk+ck+1 A2 uk+1+"°+Cn:72}-un
U U

0'

for oz = 77 32. Then

 

 

 

 

 

TAT 02 02
[[UU y1]]= a<cl:2 u1+~ +96? Jig 11k) 2 |a|[[UU Tyo”
and
2 2
Ila-00TH] -—— (ma—,1.)
a o

 

 

 

 

l/\

c 2
IaI(-:- ) nyo— w you

Since yo 9.1 R9(A)i, we have UUT YO 75 0 and

ART

. —UTUT UTU
dlst(y1,R6(A))< _ “3’1 3’1ll< _ <(:)2 llyo- yo||_
HUUT Y1“ HUUT yoll

 

 

and inequality dist(yj,’R9(A)) S <ij in (1.2.3) follows a straightforward induc-
tion. Inequality dist (xj, R9(AT)) S $2 j _ 1 can be proved similarly.
Cl
Most existing rank-revealing methods for a low rank matrix A begin with calcu-
lating the singular vector corresponding to 01 (A). For those methods, the accuracy of
computed subspaces and approxi-ranks relies heavily on the quality of singular vector
estimation [14]. A distinctive feature of our method is that it only computes vectors
in the approxi-range and the approxi-rowspace, and none of those computed vectors
needs to be singular vectors.

For finding a unit vector in the approxi-range, we use iteration (1.2.2) with stop-

6 23'
—— < 6 1.2.5
(HATyjn) m ( )

where cm is chosen on the order of the machine precision. This stopping criterion is

ping criterion

established for the following considerations: From equation (1.2.4), and by a simple

induction, we have

23' c a 23'

 

 

Cl 01 cl 01
23‘ 2i
C 0 C
+ k+1(k+1) Uk+1+...+_n_(271) un] (1.2.6)
C1 0'1 CI 01

where aj is the scalar that normalizes yj. Obviously, sequence {yj] ap—
proaches 726(A) first, and will ultimately converge to 111 if the largest singular
value 01 is strictly larger than the others. Since rankg (A) = k, 0k + 1 S 6 and 773' E
[[ATyJ-[l 3 [[AT” = 01, hence (ok+1/ol)2j S (6/nj)2j forj = 1,2,.... As-
sume (6/17h)2h < em after h steps of iteration (1.2.2). Then (oi/ol)2h < em for i =

k + 1, - -- ,n. Set p E [Ci/c1]. The distance from yh to the approxi-

max
k-l-lSiSn

range 729(A) is

 

2h 2h
C. 0'.
t+1<t+1> uk+1+,_,+§g(gg) u”

d"t A = ,
16(th736( )) 0h Cr 01 01 01

 

 

 

 

< [[pémIIk+1+p€mUk+2+m+p€mUnll = vn—kpem.

Since yo is randomly generated, p is of order 1. Thus, yh can be taken as a unit
vector in the approxi-range 729(A).
Iteration (1.2.2) can also be viewed as a power iteration on ATA. Similar to

(1.2.6), we may obtain

2j—1 C 2j—1

 

 

x] _ fl] 61 01 01 0'1
2'—1 2'—1
+Ck+1<0k+1)] v +m+£71(0_n)] W
01 01 k“ 01 01 ’

where Vi is the (right) singular vector of A associated with o,- and [83- is the

scalar that normalizes the vector xj. Similarly, {xj} converges into the approxi-
' - 1

rowspace R9(AT), and as before, condition (6/||ij]|) ‘7 < em can be used as

an stopping criterion.

1.3 Computing the approxi-range and the
approxi-rowspace

Iteration (1.2.2) produces a vector Z1 in the approxi-range 726(A). We shall show
in this section that the approxi-rank rankg(A - zlzg-A) = rank6(A) — 1. Moreover,
when the approxi-rank rank9(A) is higher than one, applying iteration ( 1.2.2) to A -
zlzirA yields another vector z2 E R6(A) that is orthogonal to z1. This deflation
process can be continued recursively to produce an orthonormal basis for the approxi-

range R9 (A).

Lemma 3. For matrix M E Rm X n with m 2 n and vector h E R", let M =
hT

IV

. [fol 2 02 _>_ 2 on are singular values of .M and 0] 2 of? Z
M
0;, are singular values of M, then 0’1 2 01 2 05 Z 02 2 2 a], 2 on.

Proof. This interlacing property is an analogical consequence of Theorem 7.3.9
in [20]. El
With the same notations as in the last section, we have the following theorem.
Theorem 4. For integerj 6 {1,2, - -- ,lc} and matrizW = [21, z2, , zj] whose
columns form an orthonormal basis for a j-dimensional subspace W of R9(A), ma-

trix A — WWTA has singular values {0,1, 0,2, - ' . ,og] satisfying

I I I
012 012022 20k_]__0'k,
I _ I _...— _
Uk—j+1 " ak—j+2" _Uk—’
and o]=o,~f0ri=k+1,~-- ,n.
Moreover, let W, = span{u’1,-~ ,uz _ j} where u] is the left singular vector of

matrix A — WWTA associated with a; for 1 S i S k — j. Then W, C R9(A) fl WL.

10

Proof. Forj = 1 and B = A — zlzirA, let {21,fig,--- ,fik} be an orthonor-

 

 

 

 

mal basis for 729(A). Also let U = [z1, fig, , fik, Uk+1l’ where Uk+1 =
[uk+1, , um]. Then
AT A A T
U AV = [Z1,U2, ,uk, Uk+1l [Avl Avk Avk+1 AVn
zil-Avl zil—Avk
lag—AV]. fig—AV]:
fill—Avl figAvk
= 0k+1
071
If
0k+1
0n
A hT
where M: with hT = [zTAv , , zTAv ]
1 1 1 k
M
and M =
rig/1v, figAvk

 

 

11

Since U and V are orthogonal matrices, the singular values
of Mare {01,og,--- ,ok}.
On the other hand, UTBV = UT (I—zlle) AV

 

 

 

 

A A T T
=[zl,ug,---,uk,Uk+1] (I—lel)
Avl Av]c Avk+1 Avn]
. 0 -
AT
“2
= 5 [AVI Avk Avk+1 Avn]
fiT
TIC
LUk+1
_M’
0k+1
0n
0
with M, = .By Lemma 3, singular values 0] Z 05 2 2 a]: of M, satisfy
IV!
I I I I
0120120g2022---20k_1ZakZUk.

Since rank(M') = k — 1, a; = 0, hence only It — 1 singular values of B are larger
than 6, and the rest of them satisfy Uh+1= 0k + 1, 056 + 2 = 0k + g, ---, 07,1: on.

Now, left singular vectors u’1,u'2,--- ,ug _ 1 of matrix B corresponding to
singular values 0] 2 0’2 2 2 0L _1 form a basis for the approxi-range
of B within 6 and zl is in the approxi-leftkernel of B = (I — zlzir) A, there-

fore, Z1 6 W"L with W’ = span {u’P-n ,uz _ 1].

12

The assertion for general 1 < j S It follows from a straightforward induction. El

Applying iteration (1.2.2) on matrix A —- 2121A yields a unit vector Zg E
729(A) that is orthogonal to zl. Continuing this process recursively, an orthonormal
basis for R9(A) will be obtained. Likewise, an orthonormal basis for the approxi-
rowspace can be obtained recursively by finding a sequence of vectors in the approxi-
rowspace.
Corollary 5. For integer j 6 {1,2, - - - ,k} and matrix Y =
[y1, yg, , yj] whose columns form an orthonormal basis for a j-dimensional

subspace y of the approxi-rowspace R9(AT), matrix A — AYYT has singular

values {0,1,05, - -- ,oiz} satisfying
0 >o’>o'>~->o' >0
1 — 1 —— 2 — — k— — k,
I _ ’ _ _ ’ — 0
Uk—j+1 _ 0k—j+2— ‘0k_ ’
and a; = oifori=k+1,~-,n.

’.

Moreover, let y’ = span {V’P-u ,vz -j} where V2

is the right singular vector of
matrix A - AYYT associated with a; for 1 g i S k — j. Then 37’ C R9(AT) 0 32¢.
From Theorem 4, one may deduce the following general rule: When a unit vec-

tor z in the approxi-range of A is obtained, let B = A — zzTA, then
0 one of the singular values of A above the rank threshold becomes zero for ma-

trix B;

o the remaining singular values of A above the threshold may shift but stay in
the same interval; and

o singular values of A below the threshold stay the same as singular values of B.

Therefore, the approxi-rank gap of the new matrix after each deflation process will
not become smaller, yielding an essential ingredient for achieving robustness in our

algorithm.

13

The unit vector produced by iteration (1.2.2) is only close to, not exactly in, the
approxi-range of A. The following proposition shows that deflation with such a vector,
the approxi-rank of the resulting matrix within the same threshold remains the same
as long as the approxi-rank gap of A is not too small.

Proposition 6. Let A E Rm X n and 6 > 0. For unit vector z E Rm and 2 being
its orthogonal projection on the approxi-leftkernel ICO(AT). Assume [[ 2 [I = e <
1. Then rank6(A — zzTA) = rank9(A) — 1 ifmin{ok - 6, 6 — 0k +1} > “A” 6.
Proof. Let E be the orthogonal projection of z on 729(A). Set (1 = E/|| E l], B = A—
zzTA, and f3 = A—ddTA. Write z = ’z‘ + E = ’z‘ +(sz)d = ’z‘ +\/1-—€2d. Let h =
am 2 n and U = [h, d] 6 am X 2. Then

zzT—ddT = (2+ 1—62d)(’z‘+ 1—52d)T—ddT

= 22T+1—62’z‘dT+ 1—62d ET—EQddT

= €2hhT+V1—€2€th+VI—€2€th—€2ddT

(2 €\/1 — £2 UT.

 

 

= U
ch—62 —62
Therefore,
[[B-B]] = (zzT—ddT) Al]
2 ,/ _ 2
g zzT—ddTHIIAII s 6 6 1 6 IIAII

 

 

e I—E2 —€2
6 V1—62

= 6 “All = 6“All

\/1—62 ——6

14

6 \/1—452

since matrix is orthogonal.
\/1—e2 -—6
Let {61,6g, ,6n} be singular values of B. As shown in Theorem 4,
those singular values satisfy 61 2 6g 2 2 (Ur—1 2 0k, &k =
0, and 5k+j = 0k+j for j = 1,2,...,n — k. By reindexing, we
write 51 2 a2 2 2 an with an = 0. Let 012 0% 2 2 at, be

the singular values of B. Then the perturbation theorem for singular values [19]

yields [5,- — 02,-] S ”A“ e for i = 1,2, - -- ,n. Consequently,

0]»,_1 2 6k_1—||A||c 2 Uk—llAl]€ > 6 > ok+1+||A||e Z 5k+llAll€ Z 02:.

C]

When the approxi-rank gap 7 = ok/ok +1 is significantly larger than one,
say 7 z 103, and the threshold 6 is not too close to the boundary of the inter-
val (‘71: + 1,0,9), Proposition 6 ensures the deflation process to be safe in our rank-

revealing algorithm.

15

1.4 The overall algorithm

Our algorithm has two main steps. We first find an approxi—range vector by the power
iteration on AAT, followed by implicit deflation via subtracting an outer product
of two vectors from matrix A. The power iteration on AAT for approximating an
approxi-range vector requires 4nmu flops, where u is the average number of iterations
per deflation step. We must emphasize here that our algorithm needs only a unit vector
in the approxi-range instead of a singular vector. From equation (1.2.6), the average
number p of power iteration steps is small for our algorithm. This may help explain
the high efliciency of our code.

Notice that if matrix A E Rm X n is too “tall” (i.e. m > n), we shall calculate
the QR factorization of A (= QR) first, and apply our algorithm on matrix R.

Our rank-revealing algorithm larank for low rank matrices can be outlined as

follows. Let matrix A E Rm X n be given along with threshold 6 > 0.

Step 0. Initialize A1 = A and i = 1.

Step 1. Find unit vectors x,- 6 729(Ag) and yi E R6(A,;) by iteration (1.2.2) on Ai'
If R6(Ai) is empty, that is, IIAzTyZH S 0, exit the algorithm.

Step 2. Set Ai + 1 = Az- — yiyiTA implicitly.

Step 3. Increase i by one and go to Step 1.

On exit, this process returns the approxi-rank k = i — 1,
bases {y1, - -- ,yk} and {x1, - -- ,xk} for 720(A) and R9(AT) respectively.

At Step 2, matrix A, +1 is implicitly obtained. It does not need to be con-
structed or stored. Matrix Az- + 1 = A — ylyIA — — yiyz-TA is stored as ma-

trix [A,Yz-] where Y,- = [y1, - -- ,yi]. When applying iteration (1.2.2) on Ai+1v we

16

use the identities

A,-T+1y = AT(y-y1ny—-~-yiy{y) = AT[y-Y7;(Y,-Ty)].
Ai+1x = Ax—yly](Ax)—---—yiyz-T(Ax) = (Ax)—Y,~[Yz-T(Ax)]

without forming Az- + 1 explicitly. The detailed pseudo—code is given in Figure 1.1.
Vector sets {y1, - -- ,yk} and {x1, - -- ,xk} produced by Algorithm larank are
almost orthonormal. The modified Gram-Schmidt method safely applies for their re-
orthogonalization, yielding orthonormal bases for approxi—range 729(A) and approxi-
rowspace 726(AT) respectively. The flop counts for the pseudo—code is (4nm —— n —

m)u(k + 1) + (4m —1)k(k +1)u assuming rank6(A) = k.

17

1.5 Numerical experiments

Our rank-revealing algorithm for low rank matrices is implemented as a Matlab pack-
age LOWRANK. We shall compare the efficiency, robustness and accuracy of our code
larank in LOWRANK with Matlab built-in svd function as well as lurv and lulv in
the UTV Tools implemented by Fierro, Hansen and Hansen [16]. All tests are carried
out in Matlab 7.0 on a Dell PC with a Pentium D CPU of 3.2 GHz, 1GB of memory,

and machine precision e z 2.2 x 10—16. The main objective of our code

machine
larank is to calculate the approxi-rank, the approxi-range, and the approxi-rowspace
of a matrix that has a low approxi-rank within a user-specified threshold.

If A E Rm x n is of approxi-rank k with threshold 6 > 0, then there is a
“noise” matrix E with ”E“ S 6 where A — E has exact rank Is. Relative pertur-
bation [IE I] / “A” is often referred to as noise level. Usually, the magnitude of relative
perturbation near machine precision, say 10—12, is taken as a low noise level, rela-
tive perturbation near 1, say 10—3, a high noise level, and the median noise level is

around 10—8. In general, the threshold 6 > 0 one chooses reflects the noise level the

matrices may have encountered.

1.5.1 Type 1: Low approxi-rank with median noise level

Matrices for this test are of size 2n x n with approxi-rank fixed at 10 within thresh-

old 6 = 10-8. The singular values range from c to “A” = 20 with approxi-

machine
rank gap 010/011 = 103. Each matrix A is constructed by using those specified
singular values to form a diagonal matrix 2 and setting A = U ZVT with randomly
generated orthogonal matrices U and V with proper sizes. We use this type of ma-
trices to test the efliciency and accuracy of our larank compared with svd, lurv,

and lulv for increasing n. For approxi—ranks, the outputs of all four algorithms are

accurate.

18

Table 1.1. Results for Type 1 matrices.

 

Matrix sizes
400 x 200 800 x 400 1600 x 800 3200 x 1600
time error time error time error time error
SVD 0.31 3e—09 2.19 4e-09 16.6 3e-09 144. 7e-09
lurv 0.66 3e-09 1.52 4e—09 5.97 3e—09 32.5 7e—09
lulv 0.56 4e-09 1.52 6e-09 6.03 5e—09 31.9 5e—09
larank 0.05 3e—09 0.11 4e—09 0.39 3e—09 1.81 4e-09

 

 

 

 

 

 

 

 

 

 

Table 1.1 only lists the time and subspace errors in executing svd, lurv, lulv
and larank. The time measures in seconds and the error measures the distance of
the computed approxi-range to the exact approxi-range. The results show that our

larank is more than ten times faster than lurv and lulv with the same accuracy.

1.5.2 Type 2: Increasing approxi-rank, fixed size and median

noise level
Matrices for this test are of size 1000 x 500. The singular values range
from emachine to ”A” = 20 with approxi-rank gap 103 , and the approxi-ranks

are set to be 10 + 203', for j = 0,1, --- , 5, within a threshold 10—8. We use this type

of matrices to test the efficiency of larank compared with lurv and lulv.

Table 1.2. Results for Type 2 matrices.

 

 

 

Average Approximate rank
Code subspace error 10 30 50 70 90 110
lurv 5e-9 2.16 5.11 8.09 11.9 15.5 16.4
lulv 6e-9 2.17 5.72 8.22 11.8 15.4 17.5
larank 4e-9 0.14 0.34 0.53 0.75 1.02 1.22

 

 

 

 

 

 

 

 

 

 

19

Results in Table 1.2 shows our larank is over ten times faster than lurv and lulv

on all cases with the same accuracy.

1.5.3 Type 3: Decreasing gap, fixed size and median noise

level

Matrices for this test are of size 1000 x 500 with approxi-rank fixed at 10 within a

threshold 10‘8. The singular values stretch from 6 to “A“ = 20. However,

machine
the approxi-rank gaps are set at 12 — 23', for j = 0, 1, - - - , 5 respectively. We use this
type of matrices to test the accuracy of larank compared with lurv and lulv by

comparing the approxi-range error which is the distance of the computed approxi-

range to the corresponding approxi-range.

Table 1.3. Results for Type 3 matrices.

 

Average Approximate rank gaps
Code time 12. 10. 8. 6. 4. 2.
lurv 2.39 4e-8 3e-8 3e-8 4e—8 1e-7 7e-4
lulv 2.21 4e-8 4e-8 6e-8 6e—8 2e-6 1e-3
larank 0.17 3e-8 3e-8 3e-8 4e-8 58-8 3e-5

 

 

 

 

 

 

 

 

 

 

 

 

The results show that even when the approxi-rank gaps are as small as 6.0, these
three codes can still produce quite accurate approxi-ranges. When the gap is 4.0, all
codes become worse while lurv and lulv have encountered tiny errors. When the

gap is 2.0, our larank still enjoys a better accuracy.

1.5.4 Type 4: High noise level with increasing size

The series of matrices in this test are of 2n x n with approxi-rank fixed at 10 within a

threshold 10_2. The singular values range from c to ”A” = 20 with approxi-

machine

20

Table 1.4. Results for Type 4 matrices.

 

Matrix sizes
400 x 200 800 x 400 1600 X 800 3200 x 1600
time error time error time error time error
lurv 0.50 5e—15 1.42 5e-15 8.27 5e—15 34.9 6e—15
lulv 0.49 5e—15 1.50 6e—15 7.80 7e-15 35.0 1e-14
larank 0.02 4e—15 0.13 4e—15 0.48 4e-15 2.08 5e-15

 

 

 

 

 

 

 

 

 

rank gap 010/011 = 103. The results show that all codes achieve accurate approxi-

ranks and approxi-ranges, while our code has a substantial advantage in efficiency.

1.5.5 Type 5: Near zero noise level, low approxi-rank, large
gap

This series of test matrices have singular values in the magnitude of machine precision
except ten singular values are in the interval [1,2]. We use threshold 10_12 to compute

approxi-ranks and approxi-ranges.

Table 1.5. Results for Type 5 matrices.

 

Matrix sizes
400 x 200 800 x 400 1600 x 800 3200 x 1600
time error time error time error time error
lurv 0.48 2e-15 1.39 2e-15 8.33 4e—15 36.4 6e—15
lulv 0.50 2e—15 1.52 2e—15 8.38 5e—15 35.8 5e—15
larank 0.02 1e—15 0.06 2e-15 0.27 3e—15 1.41 3e—15

 

 

 

 

 

 

 

 

 

The results are similar to the results of Type 4. All codes obtain accurate approxi-

ranks and approxi-ranges, while our code maintains the advantage in efficiency.

21

 

 

Algorithm larank

Input: Matrix A E Rm X n, approxi-rank threshold 6 > 0
o Initialize em = [[AIIOO (machine along with empty matrices U and V
o fork: 1,--- ,n do

0 generate a random unit vector yO, set 770 = (0 = 0
o forj=1,2,--- do

a set x = AT[yJ-_1- U(UT)’j — 1)], 773' = ”X”

 

2j — 1 ._ .
o if (%) or [7)] 77] — 1] < em then break j-loop, end

. "J
1f
. set x,- = ,ij, p = Ax -. y = p — U<UTp), (j = llyll
21' ._ .
o if 9 or C C _1 < em then break j-loop, end if
33 <-
.7
D - — 1
y] — 233'
end do

0 if 773‘ or C]- S 6 then break k-loop, end if
0 update U = [U, yj] and V = [V, xj]
o orthogonalize U and V by modified Gram-Schmidt method

end do

Output: Approxi-rank k, orthonormal bases {y1, . - - ,yk} and {X1, - - - ,xk}
for 720(A) and R9(AT) respectively.

(Operations applied on an empty vector or matrix is defined to be an empty
operation.)

Figure 1.1. Algorithm larank
22

CHAPTER 2

Updating and Downdating

Problems

2.1 The USV-plus decomposition

For a given threshold 6 > 0, we assume the singular values {Oil of matrix A E

Rm x n satisfy

012"'20k>920k+12"‘20n-
LetA=U2VT be the SVD of A. Write
A = Ak+E, (2.1.1)
where Ak = UZkVT with 2k = diag{01,---,ok,0,---,0} and E =

ospvT with 2,, = diag {0, . -- ,0, 0,, +1,” ,on}. Clearly, rank0(A) = rank(Ak) =
k and “E“ = 0k + 1 S 6. We shall call Ak the dominant part of A and call E the
noise part of A.

Matrix Ak can be considered a by-product of Algorithm larank given

23

in §1.4. While finding the approxi-rank of A, the algorithm produces a ma—
trix U = [ii1, - - - ,fik] whose columns form an orthonormal basis for the approxi-

range 729(A), and, by Theorem 4,

Thus, U U TA 2 Ah and A = U U TA+E. Similarly, if columns of matrix V form an
orthonormal basis for the approxi-rowspace R9(AT), then A = AVVT + E. Here,
we shall consider several decompositions of A induced by (practical) factorizations
of A k'

Let LQT be the transpose of the “skinny” QR—factorization of B =
U TA, where L E Rk X k is lower triangular and Q E R” x k has orthonor-
mal columns. By a straightforward argument one can easily see that the row
space R(BT), spanned by the orthonormal columns of Q, agrees with the approxi-

rowspace R9(AT). Now substituting Ak = U LQT into (2.1.1) yields
_ ~ T
A — ULQ + E, (2.1.2)

which we call a “ULV—plus decomposition” of A within 6. If the SVD of L in (2.1.2)
is L = XBYT, then
A = (735?” + E,

where U = UX and V = YTQT. We call this an “SVD-plus decomposition”
of A within 6. Let L = OR be the QR—factorization of L where R E Rk x k is

upper triangular, then (2.1.2) becomes

A = URQT+E

24

with U = U Q We shall call this a “URV-plus decomposition” of A within 6.

Those ULV/URV/SVD-plus decompositions of A defined above are convenient
tools when we deal with updating and downdating problems in the next section. The
lower-triangular matrix L, the upper-triangular matrix R and diagonal matrix B are
small for low rank matrix A. We further assign a general name for these three types
of decomposition as the “USV—plus decomposition” within 6 where “ S ” suggests
small size. Of particular importance is that when the approxi-rank k of A is small,

the computation cost of those decompositions will be low.

2.2 Updating and downdating

Suppose the approxi-rank of A has been calculated along with orthonormal bases
for the approxi—range and the approxi—rowspace. When a row/ column is inserted
in A, we wish to update all those results by taking all the available information
into account. This process is called updating [32, 33]. It is called downdating [30] if a
row / column is deleted from A. One of the main reasons for seeking alternatives to SVD
is its difliculties in benefitting from the known information when updating and down-
dating are required. Like SVD, the USV-plus decomposition of A we introduced in
the last section also contains, in addition to the approxi—rank of A, orthonormal bases
for approxi-range and approxi-rowspace of A. Therefore, in updating/downdating we
may update/downdate those results by updating/downdating the USV-plus decom-
position of A. Both of these updating and downdating procedures turn out to be
straightforward in our rank revealing method. Our extensive computing test shows
they are accurate, stable and efficient. While the existing UTV decomposition pro—
cesses robust updating capabilities, its downdating procedure seems somewhat com-

plicated [16].

25

2.2.1 Updating

For A E Rm X n with rank6(A) = k > 0, suppose one of its USV-plus de—
composition is available, say, A = U LVT + E, where U E Rm x k and V E
R" x k whose columns form an orthonormal basis for approxi-range R6(A) and

approxi-rowspace R9(AT) respectively, L E Rk x k is lower triangular with ok(L) >

6, and E E Rmxn with HE” S 6. For a E R”, let A = [1:11.] and 01 =
a

[[a — VVTa“. If a S 6, then a may be taken as a vector in the approxi-rowspace,

making the approxi-rowspaces RM?) = R9(AT), so rank9(A) = It. To update

USV-plus decomposition of A, let B = AV. Then, for b = a — VVTa,

AVVT
aTVVT

A

aT

E
bT

BVT = va =

 

 

 

 

 

A—E A
= A—
aT—bT]

 

 

Hence, A = BVT +

 

E
T ], and the “skinny” QR—factorization B = QR provides

a URV—plus decomposition of A = QRVT +

 

E
[.Ifa > 6,theni'r = 1(a—
bT a

VVTa) is a unit vector of the projection of a on approxi-kernel IC9(A). Now, let Ae =

 

 

 

 

 

 

 

 

A — E. Then

A A A6 £7 U o L o VT E

A = T = T + T = T ~T + T

a a 0 0 1 a V a v 0
is a ULV-plus decomposition of A.
L
When rank9(A) = k is small, finding SVD of TV [ is inexpensive, and
a a

 

importantly, it provides the singular value decomposition of the dominant part of A.
If a z 6, one of the singular values of A may become close to the thresh-
old 6 which may result in a smaller approxi-rank gap. Consequently, we may lose the

accuracy of approxi-range and approxi-rowspace as estimated before. In this case, we

26

apply Algorithm larank in §1.4 with input matrix A and use the vector v and the
columns of V as the initial vectors individually for power iterations to obtain a new
orthonormal base of approxi-rowspace of A. Certainly, this procedure may be used in
general when more accurate approxi-rowspace or approxi-range are required.

When a new vector is inserted, we may always consider it is inserted in the last
row by multiplying a permutation P first. On the other hand, the case of a URV-plus
decomposition of A as well as the column updating can be computed in a similar way.

We summarize the row updating algorithm lrowup in Figure 2.1 as a pseudo-code.

2.2.2 Downdating

To elaborate our downdating procedure, we need a singular value extracting strategy

Rk X k is upper triangular and Rv = on, where o is a sin-

given below. Suppose R 6
gular value of R along with corresponding unit left / right singular vectors u and v re-
spectively. We shall construct orthogonal matrices G and G as products of Givens

rotations such that

GRG =

 

R 0 ], (2.2.1)
0 o

where R 6 IR“3 _ 1)x(k _ 1) remains upper triangular. Similarly, for a lower trian-
gular matrix L E Rk X k with Lv = ou, orthogonal matrices G and G can be con-

structed such that GLG is in the form of I: L , where L E R09 _ 1))“: T 1) is

 

0’
lower triangular.

The process for constructing G and G is recursive. Let G 1 be the Givens rotation

which eliminates the first entry of v. That is, if we write 01 = [g1, - - - , gk]T then

3] '0'
T

X

G1v= g? v: .
T

tgk- -X-

 

 

 

 

27

 

Algorithm lrowup

Input: matrix A, approxi-rank k, rank threshold 6, new row aT, row in-

dex p at which aT to be inserted in A, matrices U, S, V for the USV-plus
decomposition

0 set a = “a - VVTa", i7 = 213(3 — VVTa).

o if a z 0, then apply Algorithm larank on A and use 9 and the columns
of V individually as the initial vectors for the power iterations, end if

c if a > 0, then
0 update approxi-rank k = k + 1

0 form U by inserting eT above the p—th row of U 0 .
P k + 1

S

aTV 2],V= [V v]andU=Up.

osetnewS=[

else

T above the p—th row of A.

0 form the new matrix A by inserting a
0 set W = AV, find the skinny QR factorization W = QR.

o setS=R,U=Q.

end if

Output: k , S, U, V

 

 

Figure 2.1. Algorithm lrowup

and gyv = 0. Because v and u are left and right singular vectors of RT associated
with singular value a respectively, we have RTu = av and therefore (Rgl)T u =
girRTu = agirv = 0. Only first two entries of Rgl can be nonzero since R is

upper triangular and all entries of g1 are zeros except the first two. Let Rgl =

28

[r1,r2,0, - -- ,0]T and u = [211,212, - -- ,ule. This yields r1u1+ 7‘2u2 = 0 and

 

 

’1‘1 x x x-
7‘2 X X X
301T: R[g1,...,gk] = 0 0 x x
O 0 0 x‘
c 3
Let G1 = —s c 0 be the Givens rotation withc=r1/‘/r21+r22 ands:
0 0 [ls—2

r2/ “121 + 7%. Clearly, G1 RC}— becomes upper triangular and the first entry of G 1“ is
zero. In summary, Rv = on implies G1 RGFG IV = 0G 1n and with upper triangular

matrix R1 = GlRGir, we have

0-
x x

R1 =0
-x 5x4

 

 

 

 

(5251RGIG;)(G2GIV) = R2 x = a X

 

 

 

 

where R2 = GQGI RG-ll-G; is upper triangular, and ultimately we have

~ ~ T T
Gk—1"'GIRG1"'Gk_1ek=09k-

29

The last column of the upper triangular matrix Rk _ 1 =
Gk_1---G1RGir-~G;_1 is thus [0,--- ,0,0]T as in (2.2.1). The assertion
for the lower triangular case follows the similar argument.

Now, let A be the matrix resulting from deleting the top row aT of A E
Rm X n_ The approxi-rank rank9(A) may or may not decrease. If rank6(A) =
0, then rank6(A) remains zero, requiring no further computation. For rank9(A) =
k > 0, write A = URVT+E, where columns ofU E Rm X k and V E R" X k form
an orthonormal basis for approxi-range 729(A) approxi-rowspace 729(AT) respec-
tively, R E kak is upper triangular with ak(R) > 6, and E 6 1Ran with ”E“ S
6. Since A — AVVT = E, we have

aT aT
X X

where E E R(m—1)xn is the matrix resulting from deleting the top row

VVT=E=

 

 

 

aT — aTVVT
E )

of E and “E” S “E“ S 6. Consequently, A = AVVT + E. Let AV = QR be
the “skinny” QR-factorization of AV, then
54‘ = QfiévT + E. (2.2.2)

Let 0min(R) be the smallest singular value of R. If amin(R) > 6, then rank0(A) =

k and (2.2.2) is a URV-plus decomposition of A. If amin(R) S 6, we shall ex-

tract the singular value amin(R) from R. By (2.2.1), two products of Givens ro-

 

 

~ ii 0
tations G1 and 02 exist such that GlRG2 = ~ for certain upper
0 0min(R)
triangular matrix R. It follows that
R 0 A 35 0 VT A
.4:ch ~ GJVT+E=[Ud,d] ~ “Yr +E,
0 0min(R) 0 0min(R) W

 

 

 

 

 

 

30

VT

w
wT

where [Ud, d] = (20?, Ud e RWUC— 1), d e W”,

 

]=G;VT,VWE

RnXUs - 1), and w E R". Hence,
A _ A T _ A _ ~ T
A — UdRVw + F, where F — E + 0mm(R)dw . (2.2.3)

Since w is in the approxi-rowspace of A and d is in the approxi-range of A, we
have ”F” S 6. Therefore, rank9(A) = k—l and (2.2.3) is a URV-plus decomposition
of A.

Since rank9(A) = k is small, we find the SVD of R directly which gives the left and
right singular vectors of R associated with the smallest singular value. The argument

is similar when any other row or column of A is deleted.

Our row downdating algorithm lrowdown is summarized in Figure 2.2.

31

 

Algorithm lrowdown

Input: matrix A, approxi-rank k, rank threshold 6, index p of the row to be
deleted, matrix V in the USV—plus decomposition.

0 form the new matrix A by deleting the p—th row of A, set W = AV.
0 find the skinny QR factorization of W = QR.

a find 0min(R) and the corresponding singular vector Vmin by
RANKREV [24]

o if 0min(R) > 6, then
0 set S = R, U = Q
(The approxi-rank stays at k and V does not change).
else

0 set k = k — 1 (approxi-rank reduced by one)

0 get Ud, R, and Vw as in (2.2.3) using the singular value extracting
strategy

osetS=R,U=UdandV=Vw.

end if

Output k, U, S, V.

 

 

Figure 2.2. Algorithm lrowdown

2.3 Numerical results on updating and downdat-
ing

Our updating and downdating algorithms have been extensively tested for cases of
inserting/ deleting rows or columns. Since UTV Tools [16] contains only row-updating
and row-downdating modules, we shall restrict our comparison with UTV Tools to
those situations only. The results of our method for column updating and downdating

are quite similar.

32

The two modules in UTV Tools for updating are urv_up and u1v_up accounting
for inserting a row at the bottom and two modules for downdating are urv_dw and

u1v_dw applying to deleting the top row.

2.3.1 Row-updating with increasing approxi-ranks

The test matrix is initially a 1000 x 500 matrix having an approxi-rank 10 with
threshold 10—8. The approxi-rank gap is '7 = 103. After executing larank in our
LOWRANK package and modules lurv and lulv in UTV Tools on this matrix sepa-
rately for rank-revealing, a random vector is inserted at the bottom at each updating
step. Therefore, each update will increase the approxi—rank by one. Our tests show
that our lrowup code and two counterparts u1v_up/urv.up in UTV Tools are all
accurate in identifying the increasing approxi-ranks. Table 2.1 shows the execution
time, subspace errors and the computed ranks after inserting 30 random rows. Our
lrowup appears to be considerably faster than modules in UTV Tools, while urv_up
achieves better accuracy in the updated approxi-range. For better accuracy of our
code we add one refinement step in each updating step which helps our code lrowup

achieve leading accuracy. Nonetheless, it is still faster than u1v_up and urv_up.

Table 2.1. Results for row-updating with increasing approxi—ranks

 

 

 

time (seconds) range error computed rank
ulv-up 7.08 7e-5 4O
urv_up 8.92 2e—8 4O
lrowup 0.33 7e-5 40
lrowup_1 4.64 2e-9 40

 

 

 

 

 

1rowup-1 is lrowup with one refinement step.

33

2.3.2 Row-updating without changing approxi-ranks

When approxi-rank does not change in row updating, module urv_up in UTV Tools
seems to have difficulties in identifying the approxi-ranks during the recursive up-
dating. In contrast, our code lrowup always outputs accurate approxi—ranks in all
occasions and the speed is about twice as fast on a typical example shown in Ta-
ble 2.2. The initial matrix has the same features as the example in §2.3.1 except the
approxi-rank is set at 130. A sequence of rows consisting of linear combinations of
the existing rows are inserted at the bottom one at a time. The approxi-rank stays
at 130. However, after certain steps in the recursive updating, urv_up outputs inac-

curate approxi-ranks.

Table 2.2. Results for row-updating without changing approxi-ranks

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Number of linearly dependent rows inserted
1 2 - - - 5 6 7 - - - 10

Time ulv_up 0.42 0.23 ' - - 0.30 0.25 0.28 - - - 0.24
(seconds) urv-up 0.47 0.30 - ~ - 0.38 0.23 0.33 - - - 0.28
lrowup 0.14 0.14 - -- 0.13 0.14 0.14 - - - 0.16

Approxi—range ulv-up 3e-8 4e—8 - - - 5e-8 5e-8 58-8 - - - 5e-8
error urv_up 2e—7 3e-7 ~ - - 2e—7 (0.65) (0.85) - - - (0.99)
lrowup 3e-8 4e—8 - - - 6e-8 5e-8 5e-8 - -- 6e—8

Approxi—rank ulv_up 130 130 . - ~ 130 130 130 - - - 130
output urv_up 130 130 - - - 130 (131) (131) - - - (131)

lrowup 130 130 - - - 130 130 130 - - - 130

 

 

Data in parentheses indicate inaccurate computation.

2.3.3 Row-updating with a small approxi-rank gap

This experiment compares the updating performance of our lrowup and UTV Tools

when the updated matrix has a small approxi-rank gap. The initial matrix is the same

34

as the one in §2.3.1. The inserted vector is a linear combination of existing rows plus
a random vector with a norm close to the threshold 10-8. As shown in Table 2.3,
all codes produce correct approxi-ranks, while our code lrowup takes less execution

time and obtains more accurate approxi-range.

Table 2.3. Results for row-updating with a small approxi-rank gap

 

The updated matrix has approxi-rank 11 with gap 53.8

 

 

 

 

time (seconds) range error computed rank
ulv_up 0.25 3e-3 11
urv_up 0.34 5e—7 11
lrowup 0.14 6e-8 1 1

 

 

 

 

 

2.3.4 Row-downdating without changing approxi-ranks

For the case where the approxi—rank remains invariant when a row is deleted, we
construct an initial matrix A 6 R1000 x 500 with approxi—rank 30 within thresh-
old 10_8. The approxi-rank gap is 7 = 103. Then 30 rows consisting of linear com-
binations of the existing rows of A are generated and stacked on top of A. Deleting
those rows one-by-one does not change the approxi-rank. Table 2.4 shows the results
of downdating these 30 rows recursively. The results of our code lrowdown and its
counterparts u1v_dw and urv_dw in UTV Tools are quite similar in both robustness
and accuracy, while our code lrowdown runs more than five times as fast as ulv_dw

and urv_dw.

2.3.5 Row-downdating with decreasing approxi-ranks

As mentioned in [16], UTV decomposition may have difficulties in downdating

especially when applied to the cases where the approxi-ranks are reduced. This

35

Table 2.4. Results for row-downdating without changing approxi-ranks

 

time (seconds)

range error

computed rank

 

 

ulv-dw

urv-dw

L lrowdown

 

14.6
10.2
1.80

 

1e-8
2e-9
2e—9

 

30
30
30

 

 

phenomenon does occur in the experiment shown below. We downdate a matrix
of 1010 x 500 obtained by stacking 10 random rows at the top of matrix A of
size 1000 X 500 with approxi-rank 50 within threshold 10-8. The approxi—rank gap is
set at 103. During the test, the 10 random rows are deleted one-by—one. The approxi-
rank should decrease by one at every downdating step. .

Table 2.5 shows that in step 1 to 4, downdating of the approxi-ranks were all
accurate. While all codes exhibit similar accuracy, our code lrowdown runs more
than twice as fast as ulv.dw and urv_dw. At step 5, ulv.dw miscalculates the approxi-

rank by one and this error is carried on to remaining downdating steps. Whereas, our

code lrowdown always produces the correct approxi-ranks.

Table 2.5. Results for row-downdating with decreasing approxi-ranks

 

Number of linearly dependent rows deleted

 

 

 

 

 

 

 

 

 

1 2 3 4 5 6 10
Time ulv-dw 0.48 0.33 0.33 0.33 0.28 0.30 0.44
(seconds) urv_dw 0.27 0.20 0.23 0.22 0.19 0.22 0.30
lrowdown 0. 1 1 0.09 0.08 0.08 0.06 0.09 0.08 4
Approxi-range ulv.dw 1e-8 2e-8 4e-8 4e-5 (1.0) (1.0) (1.0)
error urv_dw 1e—8 8e-9 8e-9 1e—8 1e—8 1e-8 1e-8
lrowdown 8e-9 4e—9 6e—9 4e—9 4e—9 4e—9 3e-9
Approxi—rank ulv-dw 59 58 57 56 (56) (56) (56)
output urv-dw 59 58 57 56 55 54 50
lrowdown 59 58 57 56 55 54 50

 

 

 

 

 

Data in parentheses indicate inaccurate computation.

36

 

CHAPTER 3

Applications

There are many scientific computing problems where only the dominant part of a
matrix is needed. Those problems include information retrieval and image storage to
be presented in this section. For matrix A E Rmxn’ write its USV-plus decomposition

with approxi-rank k as
A = USVT + E where U e Rmxk, S e kak, V 6 RM".

When k < n, using the dominant part U S VT as a low rank approximation to A may
reduce the memory cost by an order of magnitude and substantially cut the subse-

quential computing time, for instance, matrix—vector product y = Ax can be ap-

proximated by U [S (VTx)] with 0(n) flops instead of 0(n2) if k = 0(1).

3.1 Information retrieval: latent semantic index-
ing

A novel method called latent semantic indexing, which uses key words to find rel-
evant documents from a library database, relies critically on the computation of

low rank approximation for large matrices [4, 38]. This method can also be applied

37

to webpage search engines [3]. Assume there are m terms T1, , Tm extracted
from n documents D1, - - - , Dn. The database can be stored by a term-by-document

matrix A = (aij) E Rm x n, where
a,- j = the number of times term T,- occurs in document Dj.

In other words, the j-th column of A represents the document Dj and the i-th element
of the column is the frequency of the term Ti appears in the document. Hence, we call
the j-th column of A the document vector associated with Dj. For more sophisticated
techniques, weighted frequency strategies may be imposed on “ij [2, 12].

The matrix A is usually contaminated with a certain level of noise caused by the
presentation style, ambiguity in the use of vocabulary [25], etc. In such situations,
using the dominant part U SVT of A will achieve almost the same objective as us-
ing A itself, evidenced by our numerical test given below. In practice, k of Ak is
much smaller than min{m,n}. When a set of key words is submitted, a query vec-

T=<

tor q q1,- - - ,qn) is formed by letting

l 1, if T ,- appears in the set of key words,
(12' = ( (3.1.1)
[ 0, otherwise.

38

Table 3.1. Term-by-document matrix

 

The title of article

 

A1 Updating and downdating an upper trapezoidal sparse orthogonal factorization

A2 A rank revealing method with updating, downdating and applications

 

 

 

A3 A homotopy for solving polynomial systems
A4 UTV tools: MATLAB templates for rank revealing UTV decompositions

 

A5 Discrete orthogonal polynomials: polynomial modification of a classical functional

 

A6 Regularity results for solving systems of polynomials by homotopy method

 

A7 The polynomial rank of a commutative ring

 

 

 

A8 Orthogonal polynomials: applications and computation

 

The underlined terms are extracted to form the following 12 x 8 term—by-document matrix

 

Document
A4 A5

 

E
>
NJ
.‘>
w
:>
05
:>
q
:>
m

Term

 

application
decomposition
downdating
factorization
homotopy
method
orthogonal
polynomial
rank
revealing
system
updating

 

HOOOOv—IOOHHOO
p—IQI—IHOOHOOHO"
OHOOHOOI—‘OOOO
OCHHOOOOOOD—‘O
OOOONt—‘OOOOOO
ot—aoor—IoHv—IOOOO
OOOHI—‘OOOOOOO
OOOOHHOOOOOy—I

 

 

 

39

The query q is compared with document D ', or the document vector Ake -, by
measuring the magnitude of
qTAkej

cos6

.= __ 3.1.2
9 uqullAkeJ-l] ( l

The larger this magnitude is the more relevant the document Dj relates to the
query q.

Table 3.1 demonstrates a small size database consists of 12 terms chosen from
titles of 8 articles and the corresponding term-by-document matrix. When a user
submits a set of key words: rank, revealing, updating, downdating, application, the
associated query vector is

T
q=(101000001101)
By using rank k = 3 approximation to the term-by-document matrix, the first three
most relevant articles will be A2 (0.9136), A4 (0.7844), and A1 (0.5917), where the
number in each parenthesis indicates the cosine of the angle of the query vector and
the corresponding document vector.

From our rank-revealing method, the decomposition of Ak = U S VT, where U E
Rm X k and V E R" x k have orthonormal columns and S E Rk x k, reduces the

storage of database as well as the amount of the computations of the magnitude in

(3.1.2) as shown below. Set W = SVT E Rk x n and rewrite (3.1.2) by

cos6- = qTUSVTej , = qTU (SVTej) = qTU (Wej).
J “‘1“ llUSVTejll "Q“ llSVTejll ”‘1” “Weill

 

 

(3.1.3)

Consequently, the storage of database is reduced from mn to (m + n)k by saving

matrices U and W instead of saving the whole matrix A k' For computing cos 6?- for

40

all documents with respect to a given query, we normalize all columns of W which
requires less computation on normalizing all columns of AA? In addition, when a
term or a document is added in or removed from the database, our updating and

downdating methods can be applied.

Table 3.2. Comparisons for a 3000 X 1400 term-by-document matrix from CRAN.

 

 

 

 

threshold approxi-rank compression ratio running time (seconds)
mn .
6 k m . 1 larank lurv llllv SVD
12% x “A” 56 17.0 : 1 7.92 75.0 66.3
10% x “All 70 13.6 : 1 11.6 67.4 65.8 61.0
8% x ”A“ 90 10.6 : 1 16.5 85.6 80.3

 

 

 

 

 

 

 

 

We use a standard document collection CRAN [9, 31] to be our test sam-
ple. The collection provides about 30000 terms selected from 1400 documents. We
choose first 3000 terms from CRAN to form a term-by—document matrix A E
R3000 X 1400 and execute our algorithm larank, Matlab built-in svd function as
well as two codes lurv and lulv in UTV tools for three different prescribed thresh-
olds. The approxi-rank k’s, the compression ratios and the running time of computing
the decomposition of Ak’s are shown in Table 3.2, which illustrates the considerable
efficiency of our algorithm larank.

Using those three databases calculated by our algorithm larank and the raw
database A, we submit a query with seven key words: thick, ring, part, slight, down-
stream, yaw, and clamp to compare the retrieval results. Table 3.3 lists the indices
of the first eight relevant documents for each database. It shows that three retrieval
results from low rank databases have at least six same documents as the retrieval

result from the raw database. However, as shown in Table 3.2, the decomposition of

low rank databases requires much less storage. Table 3.4 compares the running time of

41

Table 3.3. The retrieval results for three lower rank databases and the raw database.

 

database The first 8 relevant documents

lst 2nd 3rd 4th 5th 6th 7th 8th
Ak, k = 56 766 1031 364 512 680 857 733 26
(.5507) (.5026) (.4706) (.4696) (.4541) (.4469) (.4363) (.4340)
Ak, k = 70 766 1031 512 733 680 857 943 513
(.5495) (.4908) (.4779) (.4585) (.4580) (.4365) (.4325) (.4309)
Ak, k = 90 766 1031 512 680 926 7 33 943 857
(.5536) (.4845) (.4780) (.4566) (.4359) (.4348) (.4344) (.4315)
A 766 1031 512 680 943 201 733 857
(.4880) (.4725) (.4558) (.4558) (.4364) (.4226) (.4193) (.4193)

 

 

 

 

 

 

 

 

 

The number above the parenthesis is the index j of document D -. The number in
parenthesis represents cos6- with respect to the query and the corresponding docu-
ment D -. Boldfaced numbers are the indices of the common documents in the retrieval
result from the raw database A.

adding 10 rows (terms) on the bottom of database A k with k = 56 by using our updat-
ing algorithm lrowup and two updating codes urv_up and u1v_up in UTV tools. The
comparisons for removing top 5 rows (terms) from database Ak with k = 56 are

shown in Table 3.5.

Table 3.4. Results for updating database

 

code [I lrowup urv_uP HIV—UP I]
time (seconds) I] 1.95 14.2 12.8 ]]

 

 

 

 

 

 

Comparison results for updating 10 rows on the bottom of the database A]: with k = 56.

42

Table 3.5. Results for downdating database

 

H code lrowdown urv-dw ulv-dw ]]
H time (seconds) 3.00 6.59 7.09 H

 

 

 

 

 

 

Comparison results for doumdating 5 top rows from the database Ak with k = 56.

3.2 Image processing: saving storage of pho-
tographs

An image can be stored in a matrix whose entries correspond to the levels of color
intensity [1] at pixels. In certain situations, a huge number of images need to be
archived while high resolution is not essential, like fingerprints. Our USV-plus de-
composition can greatly reduce the storage while maintaining an acceptable quality
of the images.

With a certain color map which associates a number with a level of color intensity
built in a photograph formation device such as cameras and scanners, the device
partitions an image by an m X n lattice and fits a number for each cell (pixel) according
to the color map, resulting in an m X n matrix [11, 22]. Figure 3.1 demonstrates that
an image of a baseball is partitioned by a 480 X 640 lattice and each cell corresponds

to a gray level which ranges from 0 to 255 to form a matrix A E R480 X 640.

Table 3.6. Comparison results for a 480 X 640 fingerprint image matrix.

 

 

 

threshold approxi-rank compression ratio running time (seconds)
mn .
6 k W . 1 larank lurv lulv SVD
2.1% x “All 18 15.2 : 1 0.17 2.78 2.78
1.3% X “A” 34 8.07: 1 0.34 5.31 5.28 1.97
0.8% x ”All 51 5.38 : 1 0.67 7.66 7.61

 

 

 

 

 

 

 

 

 

43

Figure 3.1. The photograph formation process: lattice partition and assignment

 

Generally, most of the singular values of photograph matrices are rela-
tively small [23]. When we truncate those terms with small singular val-
ues 0k + 1, - -- ,an from the SVD of A = 01u1v¥+--- + Ununvg, the image from
the resulting matrix Ale still maintains the main feature of the image from A. Be-
cause, by writing Ale 2 A — E, the 2—norm of E is relatively small as shown in §1.1,
thus Ak a: A. Figure 3.2 shows a lower rank approximation image of the picture in
Figure 3.1 by truncating those singular values less than 1.3% of the largest singular
value 01. We can see that the main objects and contours still can be recognizable.

Using the dominant part of matrix A reduces the storage space from mn to (m +

n)k by saving uj and ajvJ- for j = 1, - - - ,k instead of the whole m x 71. matrix.

44

Figure 3.3. The original image and three lower rank approximation images

The Original 480x640 Image Rank 18 Approximation Image

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 3.3 illustrates a 480 x 640 fingerprint photograph from FVC2004 [17] along
with three lower rank approximation images. The color map is the same as the one
used in Figure 3.1. Table 3.6 shows the compression ratio for each image and the
comparisons of the running time for computing the decomposition of matrices by
using larank, lurv, lulv, and svd. Again, the efficiency of our algorithm larank

seems to dominate the existing codes.

45

BIBLIOGRAPHY

[1] T. ACHARYA AND A. K. RAY, Image Processing Principles and Applications,
Wiley, Hoboken, NJ, 2005.

[2] M. W. BERRY, Using linear algebra for intelligent information retrieval, SIAM
Rev., 37(1995), pp. 573—595.

[3] M. W. BERRY AND M. BROWNE, Understanding Search Engines: Mathematical
Modeling and Text Retrieval, SIAM, Philadelphia, 1999.

[4] M. W. BERRY AND R. D. FIERRO, Low-rank orthogonal decompositions for

information retrieval applications, Numerical Linear Algebra with Applications.
3(1996), pp. 301-328.

[5] C. H. BISCHOF AND G. QUINTANA-ORTI, Algorithm 782: Codes for rank-
revealing QR factorizations of dense matrices, ACM TTans. Math. Software,
24(1998), pp. 254—257.

[6] A. BJDRCK, Numerical Methods for Least Squares Problems, SIAM, Philadelphia,
1996.

[7] T. R. CHAN, Rank revealing QR factorizations, Linear Algebra Appl., 88/89
(1987), pp. 67—82.

[8] T. R. CHAN AND P. C. HANSEN, Low-rank revealing QR factorizations, Numer-
ical Linear Algebra with Applications, 1(1994), pp. 33—44.

[9] CORNELL SMART SYSTEM, ftp://ftp.cs.cornelledu/pub/smart.

[10] B. H. DAYTON AND Z. ZENG, Computing the multiplicity struture in solving
polynomial systems, Proceedings of ISSAC ‘05, ACM Press, pp. 116—123, 2005.

[11] J. W. DEMMEL, Applied Numerical Linear Algebra, SIAM, Philadelphia, 1997.

46

[12] S. DEERWESTER, S. T. DUMAIS, G. W. FURNAS, T. K. LANDAUER AND

R. HARSHMAN, Indexing by latent semantic analysis, J. Amer. Soc. Inform. Sci.,
41(1990), pp. 391—407.

[13] F. DEPRETTERE, SVD and Signal Processing, Algorithms, Applications, and
Architectures, North—Holland, Amsterdam, 1988.

[14] R. D. FIERRO, AND J. R. BUNCH, Bounding the subspaces from rank revealing
two-sided orthogonal decompositions, SIAM J. Martix Analysis and Applicatons,
16(3), pp. 743—759, 1995.

[15] R. D. FIERRO, AND P. C. HANSEN, Low-rank revealing UTV decompositions,
Numer. Algorithms, 15(1997), pp. 37—55.

[16] R. D. FIERRO, P. C. HANSEN, AND P. S. K. HANSEN, UTV tools: MATLAB
templates for rank-revealing U TV decompositions, Numer. Algorithms, 20(1999),
pp. 165—194.

[17] FVC 2004 DATABASE, http://biometrics.cse.msu.edu/fvc04db.

[18] G. H. GOLUB, V. KLEMA, AND G. W. STEWART, Rank Degeneracy and Least
Squares Problems, Tech. rep. TR 456, University of Maryland, Baltimore, MD,
1976.

[19] G. H. GOLUB AND C. F. VAN LOAN, Matrix Computations, 3rd ed., Johns
Hopkins University Press, Baltimore, MD, 1996.

[20] R. A. HORN, AND C. R. JOHNSON, Matrix Analysis, Cambridge University
Press, 1985.

[21] T.-M. HWANG, W.-W. LIN AND E.K. YANG, Rank-revealing LU Factorization,
Linear Algebra and its Applications, 175(1992), pp. 115—141.

[22] A. K. JAIN, Fundamentals of Diginal Image Processing, Prentice-Hall, Engle—
wood Cliffs, NJ, 1989.

[23] S. J. LEON, Linear Algebra with Applications, Pearson Prentice Hall, NJ, 2006.

[24] T. Y. LI AND Z. ZENG, A rank-revealing method with updating, downdating, and
applications, SIAM J. Martix Analysis and Applicatons, 26(2005), pp. 918—946.

[25] C. D. MEYER, Matrix Analysis and Applied Linear Algebra, SIAM, Philadelphia,
2000.

47

[26] L. MIRANIAN AND M. GU, Strong rank revealing LUfactorization, Linear Al-
gebra and its Applications, 367(2003), pp 1—16.

[27] L. MIRSKY, Symmetric gauge functions and unitarily invariant norms, Quart.
J. Math. Oxford Ser. 11(1960), pp. 50—59.

[28] M. MOONEN AND B. DE MOOR, S VD and Signal Processing, III, Algorithms,
Applications, and Architectures, Elsevier, Amsterdam, 1995.

[29] C.-T. PAN, 0n the existence and computation of rank-revealing LU factoriza—
tions, Linear Algebra and its Applications, 316(2000), pp. 199—222.

[30] H. PARK AND L. ELDEN, Downdating the Rank Revealing URV Decomposition,
SIAM J. Martix Analysis and Applicatons, 16, pp. 138—155, 1995.

[31] G. SALTON AND M. MCGILL, Introduction to Modern Information Retrieval,
McGraw—Hill, New York, 1983.

[32] G. W. STEWART, An Updating Algorithm for Subspace Tracking, IEEE Trans.
Signal Processing, 40, pp. 1535—1541, 1992.

[33] G. W. STEWART, Updating a Rank-Revealing ULV Decompositions, SIAM J.
Martix Analysis and Applicatons, 14, pp. 494—499, 1993.

[34] G. W. STEWART, U TV decompositions, in Numerical Analysis 1993, D. F. Grif-
fith and G. A. Watson, eds, Pitman Res. Notes Math. Ser. 303, Longman, Harlow,
UK 1994, pp. 225—236.

[35] G. W. STEWART, Matrix Algorithms: Basic Decompositions, SIAM, Philadel-
phia, 1998.

[36] R. VACCARO, SVD and Signal Processing, II, Algorithms, Applications, and
Architectures, Elsevier, Amsterdam, 1991.

[37] Z. ZENG, Computing multiple roots of inexact polynomials, Mathematics of Com-
putation, 74, pp. 869—903, 2005.

[38] H. ZHA, AND H. D. SIMON, Timely communication on updating problems in
latent semantic indexing, SIAM J. Sci. Comput., 21(2), pp. 782—791, 1999.

48

 

u[[[[[[[[[[[[.