Time and space efficient generators for quasiseparable matrices

0.33 MB·

by  Clement Pernet

1
ffi
Time and space e cient generators for quasiseparable
matrices
Cle´mentPernet
Universite´GrenobleAlpes,LaboratoireJeanKuntzmann,
7
E´NSdeLyon,Laboratoiredel’InformatiqueduParalle´lisme.
1
0
ArneStorjohann
2
DavidR.CheritonSchoolofComputerScience
n
UniversityofWaterloo,Ontario,Canada,N2L3G1
a
J

2

]
C
S Abstract
.
s Theclassofquasiseparablematricesisdefinedbythepropertythatanysubmatrixen-
c
tirelybeloworabovethemaindiagonalhassmallrank,namelybelowaboundcalled
[
theorderofquasiseparability. ThesematricesarisenaturallyinsolvingPDE’sforpar-

1
ticle interaction with the Fast Multi-pole Method (FMM), or computing generalized
v
eigenvalues. From these application fields, structured representations and algorithms
6
9 havebeendesignedinnumericallinearalgebratocomputewiththesematricesintime
3 linear in the matrix dimension and either quadratic or cubic in the quasiseparability
0 order. Motivated by the design of the general purpose exact linear algebra library
0
LinBox, and by algorithmic applications in algebraic computing, we adapt existing
.
1 techniquesintroducenovelonestousequasiseparablematricesinexactlinearalgebra,
0 where sub-cubic matrix arithmetic is available. In particular, we will show, the con-
7
nection between the notion of quasiseparability and the rank profile matrix invariant,
1
thatwehaveintroducedin2015. Itresultsintwonewstructuredrepresentations,one
:
v being a simpler variation on the hierarchically semiseparable storage, and the second
i
X one exploiting the generalized Bruhat decomposition. As a consequence, most basic
operations,suchascomputingthequasiseparabilityorders,applyingavector,ablock
r
a vector, multiplying two quasiseparable matrices together, inverting a quasiseparable
matrix,canbeatleastasfastandoftenfasterthanpreviousexistingalgorithms.
Keywords: Quasiseparable;HierarchicallySemiseparable;Rankprofilematrix;
GeneralizedBruhatdecomposition;Fastmatrixarithmetic.
Emailaddresses:[email protected](Cle´mentPernet),[email protected]
(ArneStorjohann)
URL:http://ljk.imag.fr/membres/Clement.Pernet/(Cle´mentPernet),
PreprintsubmittedtoElsevier January3,2017
1. Introduction
Weconsidertheclassofquasiseparablematrices,definedbyaboundingcondition
ontheranksofthesubmatricesintheirloweranduppertriangularparts. Thesestruc-
turedmatricesoriginatemainlyfromtwodistinctapplicationfields: computinggener-
alizedeigenvalues(Gohbergetal.,1985;Vandebriletal.,2005,2007),andsolvingpar-
tialdifferentialequationsforparticulesimulationwiththefastmultipolemethod(Car-
rier et al., 1988). This class also arise naturally, as it includes the closure under in-
version of the class of banded matrices. Among the several definitions used in the
litterature,wewillusethatofEidelmanandGohberg(1999)fortheclassofquasisep-
arablematrices.
Definition1. Ann×nmatrixMis(r ,r )-quasiseparableifitsstrictlylowerandupper
L U
triangularpartssatisfythefollowinglowrankstructure: forall1≤k≤n−1,
rank(Mk+1..n,1..k) ≤ rL, (1)
rank(M1..k,k+1..n) ≤ rU. (2)
Thevaluesr andr arethequasiseparableordersofM.
L U
OtherpopularclassesofstructuredmatriceslikeToeplitz,Vandermonde,Cauchy,
Hankelmatricesandtheirblockversions,enjoyaunifieddescriptionthroughthepow-
erfulnotionofdisplacementrank(Kailathetal.,1979).Consequentlytheybenefitfrom
spaceefficientrepresentations(linearinthedimensionnandinthedisplacementrank
s), and time efficient algorithms to apply them to a vector, compute their inverse and
solvelinearsystems:mostoperationshavebeenreducedtopolynomialarithmetic(Pan,
1990;BiniandPan,1994),andbyincorporatingfastmatrixalgebra,thiscosthasbeen
reduced from O˜(s2n) to O˜(sω−1n) by Bostan et al. (2008) (assuming that two n×n
matricescanbemutlipliedinO(nω)for2.3728639≤ω≤3(LeGall,2014)).
Howeverquasiseparablematricesdonotfitintheframeworkofrankdisplacement
structures. Taking advantage of the low rank properties, mainly two types of struc-
turedrepresentationshavebeendeveloppedtogetherwithcorrespondingdedicatedal-
gorithmstoperformcommonlinearalgebraoperations: thequasiseparablegenerators
of Vandebril et al. (2005, 2007), sometimes referred to as sequentially semiseparable
(SSS) and the hierarchically semiseparable representations (HSS) of Chandrasekaran
et al. (2006); Xia et al. (2010). We refer to (Vandebril et al., 2005) and (Vandebril
etal.,2007)forabroadbibliographicoverviewonthetopic. Notealsothealternative
approachofGivensandunitaryweightsinDelvauxandVanBarel(2007).
SequentiallySemiseparablerepresentation. Thesequentiallysemiseparablerepresen-
tationusedbyEidelmanandGohberg(1999);Vandebriletal.(2005,2007);Eidelman
et al. (2005); Boito et al. (2016) for a matrix M, consists of (n−1) pairs of vectors
p(i),q(i)ofsizer ,(n−1)pairsofvectorsg(i),h(i)ofsizer ,n−1matricesa(i)of
L U
https://cs.uwaterloo.ca/~astorjoh/(ArneStorjohann)
2
dimensionr ×r ,andn−1matricesb(i)ofdimensionr ×r suchthat
L L U U

Mi,j = pdg(((iii))),TTab>i<jqh((jj)),, 111≤≤≤iij=<<ijj≤≤≤nnn
ij
where
a>ij =a(i−1)...a(j+1)for j>i+1, aj+1,j =1,
b<ij =b(i+1)...b(i−1)fori> j+1, bi,i+1 =1.
For s = max(r ,r ), this representation, of size O(n(r2 + r2)) = O(s2n) makes it
L U L U
possibletoapplyavectorinO(s2n)fieldoperations,multiplytwoquasiseparablema-
tricesintimeO(s3n)andalsocomputetheinverseofastronglyregularmatrixintime
O(s3n)(EidelmanandGohberg,1999).
The Hierarchically Semiseparable representation. The Hierarchically Semiseparable
representation was introduced in Chandrasekaran et al. (2006) and is related to the
structure used in the Fast Multipole Method (Carrier et al., 1988). It is based on the
splitting of the matrix in four quadrants, the use of rank revealing factorizations of
its off-diagonal quadrants and applying the same scheme recursively on the diagonal
blocks. Afurthercompressionisappliedtorepresentalloff-diagonalblocksaslinear
combinations(calledtranslationoperators)ofblocksofafinerrecursiveorder. While
the space and time complexity of the HSS representation is depending on numerous
parameters,theanalysisinChandrasekaranetal.(2006)seemtoindicatethatthesize
ofanHSSrepresentationisO(sn),itcanbeappliedtoavectorinlineartimeinitssize,
andlinearsystemscanbesolvedinO(s2n). FortheproductoftwoHSSmatrices,we
couldnotfindanybetterestimatethanO(s3n)givenbyShengetal.(2007).
Context and motivation. The motivation here is to propose simplified and improved
representationsofquasiseparablematrices(inspaceandtime). Ourapproachdoesnot
focusonnumericalstabilityforthemoment. Ourfirstmotivationisindeedtousethese
structured matrices in computer algebra where computing e.g. over a finite field or
over multiprecision integers and rationals does not lead to any numerical instability.
Hence we will assume throughout the paper that any Gaussian elimination algorithm
mentioned has the ability to reveal ranks. In numerical linear algebra, a special care
needtobetakenforthepivotingofLUdecompositions(Hwangetal.,1992;Pan,2000),
andQRorSVDdecompositionsareoftenpreferred(Chan,1987;Chandrasekaranand
Ipsen,1994). Partofthemethodspresentedhere,namelythatofSection5,relyonan
arbitrary rank revealing matrix factorization and can therefore be applied to a setting
with numerical instability. In the contrary, Section 6 relies on a class of Gaussian
eliminationalgorithmthatrevealtherankprofilematrix,henceapplyingittonumerical
setting is future work. This study is motivated by the design of new algorithms on
polynomialmatricesoverafinitefield,wherequasiseparablematricesnaturallyoccur,
andmoregenerallybytheframeworkoftheLinBoxlibrary(TheLinBoxGroup,2016)
forblack-boxexactlinearalgebra.
3
Contribution. This paper presents in further details and extends the results of Pernet
(2016), while also fixing a mistake 2. It proposes two new structured representations
forquasiseparablematrices,aRecursiveRankRevealing(RRR)representationthatcan
be viewed as a simplified version of the HSS representation of Chandrasekaran et al.
(2006),andarepresentationbasedonthegeneralizedBruhatdecomposition,whichwe
nameCompactBruhat(CB)representation.Thelaterone,ismadepossiblebythecon-
nection that we make between the notion of quasiseparability and a matrix invariant,
the rank profile matrix, that we introduced in Dumas et al. (2015) and applied to the
generalized Bruhat decomposition in Dumas et al. (2016). More precisely, we show
that the lower and upper triangular parts of a quasiseparabile matrix have a General-
izedBruhatdecompositionsoffofwhichmanycoefficientscanbeshaved. Theresult-
ingstructureofthesedecompositionsallowstohandlethemwithinmemoryfootprint
andtimecomplexitythatdoesnotdependontherankbutonthequasiseparableorder
(which can be arbitrarily lower). These two representations use respectively a space
O(snlogn)(RRR)andO(sn)(CB),henceimprovingoverthatoftheSSS,O(s2n),and
s
matchingthatoftheHSSrepresentation,O(sn).
The complexity of applying a vector remains linear in the size of the represen-
tations. The main improvement in these two representations is in the complexity of
applyingthemtomatricesandcomputingthematrixinverse,wherewereplacebysω−1
thes3factoroftheSSSorthes2factoroftheHSSrepresentations.3 Table1compares
SSS HSS RRR CB
Size O(s2n) O(sn) O(snlogn) O(sn)
s
Construction O(s2n2) O(sn2) O(sω−2n2) O(sω−2n2)
QSxVec O(s2n) O(sn) O(snlogn) O(sn)
s
QSxTS O(s3n) O(s2n) O(sω−1nlogn) O(sω−1n)
s
QSxQS O(s3n) O(s3n) O(sω−1nlog2 n) O(sω−2n2)
s
LinSys O(s3n) O(s2n) O(sω−1nlog2 n)
s
Table1:ComparingthesizeandtimecomplexitiesforbasicoperationsoftheproposedRRRandCBrepre-
sentationswiththeexistingoneSSSandHSSonann×nquasiseparablematrixoforders.
the two proposed representations with the SSS and the HSS in their the size, and the
complexityofthemainbasicoperations.
Outline. Section2definesandrecallssomepreliminarynotionsonlefttriangularma-
trices and the rank profile matrix, that will be used in Section 3 and 6. Using the
strongconnectionbetweentherankprofilematrixandthequasiseparablestructure,we
firstproposeinSection3analgorithmtocomputethequasiseparabilityorders(r ,r )
L U
2Equation(9)inPernet(2016)ismissingtheLeftoperators.Theresultingalgorithmsareincorrect.This
isfixedinsection6.2.
3NotethatmostcomplexitiesforSSSandHSSinthelitteraturearegivenintheformO(n2)orO(n),
consideringtheparametersasaconstant.Theestimatesgivenhere,withtheexponentins,canbefoundin
theproofsoftherelatedpapersoreasilyderivedfromthealgorithms.
4
of any dense matrix in O(n2sω−2) where s = max(r ,r ). Section 4 then describes
L U
thetwoproposedstructuredrepresentationsforquasiseparablematrices: theRecursive
RankRevealingrepresentation(RRR),asimplifiedHSSrepresentationbasedonabi-
narytreeofrankrevealingfactorizations,andtheCompactBruhatrepresentation(CB),
basedontheintermediateBruhatrepresentation. Section5thenpresentsalgorithmsto
computeanRRRrepresentation,andperformthemostcommonoperationswithit:ap-
plying a vector, a tall and skinny matrix, multiplying two quasiseparable matrices in
RRRrepresentation,andcomputingtheinverseofastronglyregularRRRmatrix. Sec-
tion6presentsalgorithmstocomputeaCompactBruhatrepresentation,andmultiply
itwithavector,atallandskinnymatrixoradensematrix.
Notations. Throughoutthepaper,A willdenotethesub-matrixofAofrowindices
i..j,k..l
betweeniand jandcolumnindicesbetweenkandl. Thematrixofthecanonicalbasis,
withaoneatposition(i, j)willbedenotedby∆(i,j). Wewilldenotetheidentitymatrix
of order n by I , the unit antidiagonal of dimension n by J and the zero matrix of
n n
dimensionm×nby0 .
m×n
2. Preliminaries
2.1. Lefttriangularmatrices
Wewillmakeintensiveuseofmatriceswithnon-zeroelementsonlylocatedabove
themainanti-diagonal. Wewillrefertothesematricesaslefttriangular,toavoidany
confusionwithuppertriangularmatrices.
Definition2. Anm×nmatrixAislefttriangularifA =0foralli>n− j.
i,j
The left triangular part of a matrix A, denoted by Left(A) will refer to the left
triangular matrix extracted from it. We will need the following property on the left
triangularpartoftheproductofamatrixbyatriangularmatrix.
Lemma3. LetA = BUbeanm×nmatrixwhereUisn×nuppertriangular. Then
Left(A)=Left(Left(B)U).
Proof. LetA¯ = Left(A),B¯ = Left(B). For j ≤ n−i,wehaveA¯ = (cid:80)n B ·U =
i,j k=1 i,k k,j
(cid:80)j B ·U asUisuppertriangular. Nowfork≤ j≤n−i,B =B¯ ,whichproves
k=1 i,k k,j i,k i,k
thatthelefttriangularpartofAisthatofLeft(B)U.
ApplyingLemma3onAT yieldsLemma4
Lemma4. LetA = LBbeanm×nmatrixwhereLism×mlowertriangular. Then
Left(A)=Left(LLeft(B)).
Lastly, we will extend the notion of order of quasiseparability to left triangular
matrices,inthenaturalway: theorderofleftquasiseparabilityisthemaximalrankof
anyleadingk×(n−k)sub-matrix. Whennoconfusionmayoccur,wewillabusethe
definitionandsimplycallittheorderofquasiseparability.
5
2.2. PLUQdecomposition
Werecallthatforanym×nmatrixAofrankr,thereexistaPLUQdecomposition
A=PLUQwherePisanm×mpermutationmatrix,Qisann×npermutationmatrix,
L is an m×r unit lower triangular matrix, and U is an r×n upper triangular matrix.
matrix. It is not unique, but once the permutation matrices P and Q are fixed, the
triangularfactorsLandUareunique,sincethematrixPTAQT hasgenericrankprofile
andthereforehasauniqueLUdecomposition.
2.3. Therankprofilematrix
Wewilluseamatrixinvariant,introducedin(Dumasetal.,2015,Theorem1),that
summarizestheinformationontheranksofanyleadingsub-matricesofagiveninput
matrix.
Definition 5. (Dumas et al., 2015, Theorem 1) The rank profile matrix of an m×n
matrixAofrankristheuniquem×nmatrixR ,withonlyrnon-zerocoefficients,all
A
equaltoone,locatedondistinctrowsandcolumnssuchthatanyleadingsub-matrices
ofR hasthesamerankasthecorrespondingleadingsub-matrixinA.
A
This invariant can be computed in just one Gaussian elimination of the matrix A,
atthecostofO(mnrω−2)fieldoperations(Dumasetal.,2015), providedsomecondi-
tionsonthepivotingstrategybeingused. ItisobtainedfromthecorrespondingPLUQ
decompositionastheproduct
(cid:34) (cid:35)
I
R =P r Q.
A 0
(m−r)×(n−r)
We also recall in Theorem 6 an important property of such PLUQ decompositions
revealingtherankprofilematrix.
Theorem6((Dumasetal.,2016,Th.24),(Dumasetal.,2013,Th.1)). LetA=PLUQ,
(cid:104) (cid:105)
aPLUQdecompositionrevealingtherankprofilematrixofA.Then,P L 0m×(m−r) PT
(cid:34) (cid:35)
U
islowertriangularandQT Qisuppertriangular.
0
(n−r)×n
Lemma7. Therankprofilematrixinvariantispreservedbymultiplication
1. totheleftwithaninvertiblelowertriangularmatrix,
2. totherightwithaninvertibleuppertriangularmatrix.
Proof. LetB=LAforaninvertiblelowertriangularmatrixL. Thenforanyi≤m, j≤
n,rank(B )=rank(L A )=rank(A ). HenceR =R .
1..i,1..j 1..i,1..i 1..i,1..j 1..i,1..j B A
3. Computingtheordersofquasiseparability
LetMbeann×nmatrixofwhichonewantstodeterminethequasiseparableorders
(r ,r ). LetLandUberespectivelythelowertriangularpartandtheuppertriangular
L U
partofM.
6
MultiplyingontheleftbyJ ,theunitanti-diagonalmatrix,inversestheroworder
n
while multiplying on the right by J inverses the column order. Hence both J L and
n n
UJ arelefttriangularmatrices.Remarkthatconditions(1)and(2)statethatallleading
n
k×(n−k)sub-matricesofJ LandUJ haveranknogreaterthanr andr respectively.
n n L U
We will then use the rank profile matrix of these two left triangular matrices to find
theseparameters.
3.1. Fromarankprofilematrix
First,notethattherankprofilematrixofalefttriangularmatrixisnotnecessarily
(cid:20) (cid:21) (cid:20) (cid:21)
110 100
left triangular. For example, the rank profile matrix of 100 is 010 . However,
000 000
only the left triangular part of the rank profile matrix is sufficient to compute the left
quasiseparableorders.
Suppose for the moment that the left-triangular part of the rank profile matrix of
a left triangular matrix is given (returned by a function LT-RPM). It remains to enu-
merate all leading k×(n−k) sub-matrices and find the one with the largest number
ofnon-zeroelements. Algorithm1showshowtocomputethelargestrankofalllead-
ing sub-matrices of such a matrix. Run on J L and UJ , it returns successively the
n n
quasiseparableordersr andr .
L U
Algorithm1QS-order
Require: A,ann×nmatrix
Ensure: max{rank(A ):1≤k≤n−1}
1..k,1..n−k
R←LT-RPM(A) (cid:46)ThelefttriangularpartoftherankprofilematrixofA
rows←(False,...,False)
cols←(False,...,False)
forall(i, j)suchthatR =1do
i,j
rows[i]←True
cols[j]←True
endfor
s,r←0
fori=1...n−1do
if rows[i]thenr←r+1
if cols[n−i+1]thenr←r−1
s←max(s,r)
endfor
returns
This algorithm runs in O(n) provided that the rank profile matrix R is stored in a
compactway,e.g. usingavectorofrpairsofpivotindices([(i , j ),...,(i , j )].
1 1 r r
3.2. Computingtherankprofilematrixofalefttriangularmatrix
Wenowdealwiththemissingcomponent: computingthelefttriangularpartofthe
rankprofilematrixofalefttriangularmatrix.
7
3.2.1. FromaPLUQdecomposition
A first approach is to run any Gaussian elimination algorithm that can reveal the
rankprofilematrix, asdescribedinDumasetal.(2015). Inparticular, thePLUQde-
compositionalgorithmofDumasetal.(2013)computestherankprofilematrixofAin
O(n2rω−2)wherer =rank(A). Howeverthisestimatemaybepessimisticasitdoesnot
takeintoaccountthelefttriangularshapeofthematrix. Moreover,itdoesnotdepend
ontheleftquasiseparableordersbutontherankr,whichcouldbemuchhigher.
Remark 8. The discrepancy between the rank r of a left triangular matrix and its
quasiseparable order arises from the location of the pivots in its rank profile matrix.
Pivotslocatednearthetopleftcornerofthematrixaresharedbymanyleadingsub-
matrices, and are therefore likely to contribute to the quasiseparable order. On the
otherhand,pivotsnearthemainanti-diagonalcanbenumerous,butdonotaddupto
a large quasiseparable order. As an illustration, consider the two following extreme
cases:
1. amatrixAwithgenericrankprofile. Thentheleadingr×rsub-matrixofAhas
rankrandthequasiseparableorderiss=r.
2. thematrixwithn−1onesimmediatelyabovethemainanti-diagonal.Ithasrank
r=n−1butquasiseparableorder1.
Remark8indicatesthatintheunluckycaseswhenr (cid:29) s,thecomputationshould
reduce to instances of smaller sizes, hence a trade-off should exist between, on one
hand,thediscrepencybetweenr and s,andontheotherhand,thedimensionnofthe
problems. Allcontributionspresentedintheremainingofthepaperarebasedonsuch
trade-offs.
3.2.2. Adedicatedalgorithm
Inordertoreachacomplexitydependingon sandnotr,weadaptinAlgorithm2
thetilerecursivealgorithmofDumasetal.(2013),sothatthelefttriangularstructure
oftheinputmatrixispreservedandcanbeusedtoreducetheamountofcomputation.
Algorithm 2 does not assume that the input matrix is left triangular, as it will be
called recursively with arbitrary matrices, but guarantees to return the left triangular
partoftherankprofilematrix. WhilethetopleftquadrantA iseliminatedusingany
1
PLUQ decomposition algorithm revealing the rank profile matrix, the top right and
bottomleftquadrantsarehandledrecursively.
Theorem9. Givenann×ninputmatrixAwithleftquasiseparableorders,Algorithm2
computes the left triangular part of the rank profile matrix of A in O(n2sω−2) field
operations.
Proof. Firstremarkthat
(cid:34) (cid:35) (cid:34) (cid:35) (cid:34) (cid:35)
D L −1 B
P =P 1 P TP 1 =LA .
1 F 1 −M L −1 I 1 1 B 2
(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)1(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)1(cid:32)(cid:32) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)n(cid:32)(cid:32)−(cid:32)(cid:32)(cid:32)r(cid:32)(cid:32)1(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) 2
(cid:124) (cid:123)(cid:122) (cid:125)
L
8
Algorithm2LT-RPM:LeftTriangularpartoftheRankProfileMatrix
Require: A: ann×nmatrix
Ensure: R: thelefttriangularpartoftheRPMofA
1: if n=1thenreturn[0]
(cid:34) (cid:35)
A A
2: SplitA= A1 2 whereA3is(cid:98)n2(cid:99)×(cid:98)n2(cid:99)
3
(cid:34) (cid:35)
L (cid:104) (cid:105)
3: ComputeaPLUQdecompositionA1 =P1 M1 U1 V1 Q1revealingtheRPM
1
(cid:34) (cid:35)
I
4: R1 ←P1 r1 0 Q1wherer1 =rank(A1).
(cid:34) (cid:35)
B
5: B1 ←P1TA2
2
 
6: (cid:104)C1 C2(cid:105)←A3Q1T (cid:46)HereA= L1M\U11 V01 BB12 .
C C
1 2
7: D←L1−1B1
8: E←C1U1−1
9: F←B2−M1D  
10: G←C2−EV1 (cid:46)HereA= L1M\U11 V01 DF .
E G
(cid:34) (cid:35)
11: H←P1 0rF1×2n
(cid:104) (cid:105)
12: I← 0r1×2n G Q1
13: R2 ←LT-RPM(H)
14: R3 ←LT-RP(cid:34)M(I) (cid:35)
R R
15: returnR← 1 2
R
3
Hence
(cid:34) (cid:104) (cid:105) (cid:35)
(cid:104) (cid:105) U V Q D
L A A =P 1 1 1 .
1 2 1
0 F
FromTheorem6,thematrixLislowertriangularandbyLemma7therankprofilema-
(cid:34) (cid:104) (cid:105) (cid:35)
(cid:104) (cid:105) U V Q D
trixof A A equalsthatofP 1 1 1 . NowasU isuppertriangular
1 2 1 1
0 F
(cid:34) (cid:104) (cid:105) (cid:35)
U V Q 0
andnon-singular,thisrankprofilematrixisinturnthatofP 1 1 1 and
1
0 F
(cid:104) (cid:105)
itslefttriangularpartis R R .
1 2
(cid:104) (cid:105)T
By a similar reasoning, R R is the left triangular part of the rank profile
1 3
(cid:104) (cid:105)T
matrixof A A ,whichshowsthatthealgorithmiscorrect.
1 3
Let s be the left quasiseparable order of H and s that of I. The number of field
1 2
9
operationsrequiredtorunAlgorithm2is
T(n,s)=αrω−2n2+T (n/2,s )+T (n/2,s )
1 LT-RPM 1 LT-RPM 2
forapositiveconstantα. WewillprovebyinductionthatT(n,s)≤2αsω−2n2.
Again, since L is lower triangular, the rank profile matrix of LA is that of A
2 2
andthequasiseparableordersofthetwomatricesarethesame. NowHisthematrix
LA withsomerowszeroedout,hence s ,thequasiseparableorderofHisnogreater
2 1
than that of A which is less or equal to s. Hence max(r ,s ,s ) ≤ s and we obtain
2 1 1 2
T(n,s)≤αsω−2n2+4αsω−2(n/2)2 =2αsω−2n2.
4. Newstructuredrepresentationsforquasiseparablematrices
Inordertointroducefastmatrixarithmeticinthealgorithmscomputingwithqua-
siseparablematrices,weintroduceinthissectionthreenewstructuredrepresentations:
the Recursive Rank Revealing (RRR) representation, the Bruhat representation, and
finallyitscompactversion,theCompactBruhat(CB)representation.
4.1. TheRecursiveRankRevealingrepresentation
This a simplified version of the HSS representation. It uses in the same manner
a recursive splitting of the matrix in a quad-tree, and each off-diagonal block at each
recursivelevelisrepresentedbyarankrevealingfactorization.
Definition10(RR:Rankrevealingfactorization). Arankrevealingfactorization(RR)
ofanm×nmatrixAofrankr isapairofmatricesLandRofdimensionsm×r and
r×nrespectively,suchthatA=LR.
For instance, a PLUQ decomposition is a rank revealing factorization. One can
eitherstoreexplicitelythetwofactorsPLandUQoronlyconsiderthefactorsLandU
keepinginmindthatpermutationsneedtobeappliedontheleftandontherightofthe
product.
Definition11(RRR:RecursiveRankRevealingrepresentation). Arecursiverankre-
(cid:34) (cid:35)
A A
vealing (RRR) representation of an n×n quasiseparable matrix A = 11 12 of
A A
21 22
order sisformedbyarankrevealingfactorizationofA andA andappliesrecur-
12 21
sivelyfortherepresentationofA andA .
11 22
ARecursivelyRankRevealingrepresentationformsabinarytreewhereeachnode
correspond to a diagonal block of the input matrix, and contains the Rank Revealing
factorizationofitsoff-diagonalquadrants.
If A is (r ,r )-quasiseparable, then all off-diagonal blocks in its lower part have
L U
rankboundedbyr ,andtheirrankrevealingfactorizationstakeadvantageofthislow
L
rankuntilablockdimensionn/2k ≈r whereadenserepresentationisused. Thesame
L
appliesfortheuppertriangularpartwithquasiseparableorderr . Thisrepresentation
U
usesO(snlogn)spacewheres=max(r ,r ).
s L U
10