Fast Computation of the Rank Profile Matrix and the Generalized Bruhat Decomposition Jean-Guillaume Dumas Universit´e Grenoble Alpes, Laboratoire LJK, umr CNRS, BP53X, 51, av. des Math´ematiques, 6 F38041 Grenoble, France 1 0 2 Cl´ement Pernet n a Universit´e Grenoble Alpes, Laboratoire de l’Informatique du Parall´elisme, Universit´e de Lyon, J France. 8 Ziad Sultan ] C S Universit´e Grenoble Alpes, Laboratoire LJK and LIG, Inria, CNRS, Inovall´ee, 655, av. de . l’Europe, F38334 St Ismier Cedex, France s c [ 1 Abstract v 8 9 The row (resp. column) rank profile of a matrix describes the stair-case shape of its row (resp. 7 column) echelon form. We here propose a new matrix invariant, the rank profile matrix, sum- 1 marizing all information on the row and column rank profiles of all the leading sub-matrices. 0 We show that this normal form exists and is unique over any ring, provided that the notion . 1 of McCoy’s rank is used, in the presence of zero divisors. We then explore the conditions for a 0 Gaussianeliminationalgorithmtocomputeallorpartofthisinvariant,throughthecorrespond- 6 ingPLUQdecomposition.ThisenlargesthesetofknownEliminationvariantsthatcomputerow 1 orcolumn rank profiles.AsaconsequenceanewCrout basecase variant significantly improves : thepracticalefficiencyofpreviouslyknownimplementationsoverafinitefield.Withmatricesof v i very small rank, we also generalize the techniques of Storjohann and Yang to the computation X of the rank profile matrix, achieving an (rω+mn)1+o(1) time complexity for an m×n matrix r of rank r,where ω is theexponentof matrix multiplication. Finally, bygive connections tothe a Bruhat decomposition, and several of its variants and generalizations. Thus, our algorithmic improvements for the PLUQ factorization, and their implementations, directly apply to these decompositions. In particular, we show how a PLUQ decomposition revealing the rank profile matrix also reveals both a row and a column echelon form of the input matrix or of any of its leading sub-matrices, by a simple post-processing made of row and column permutations. Key words: Gaussian elimination, Rank profile, Echelon form, PLUQ decomposition, Bruhat decomposition, McCoy’s rank. Preprint submitted to Journal of Symbolic Computation 4 December 2015 Contents 1 Introduction 2 2 The rank profile matrix 5 2.1 Definition over a field 5 2.2 Generalization overa ring 7 3 When does a PLUQ algorithm reveal the rank profile matrix? 11 3.1 Ingredientsof a PLUQ decomposition algorithm 11 3.1.1 Pivot search 12 3.1.2 Pivot permutation 13 3.2 How to reveal rank profiles 14 4 Algorithms for therank profile matrix 17 4.1 Iterativealgorithms 17 4.1.1 Rowand Column order Search 17 4.1.2 Lexicographic order based pivot search 18 4.1.3 Product order based pivot search 18 4.2 Recursivealgorithms 19 4.2.1 Slab recursive algorihtms 19 4.2.2 Tile recursive algorithms 19 5 Improvementsin practice 19 6 Relations with other triangularizations 23 6.1 The LEU decomposition 23 6.2 The Bruhat decomposition 24 6.3 Relation to LUPand PLU decompositions 25 6.4 Computing Echelon forms 25 6.5 The generalized Bruhat decomposition 27 7 Improvementfor low rank matrices 27 7.1 Storjohann and Yang’s algorithm 27 7.2 OnlineLU decomposition 29 References 30 1. Introduction Triangular matrix decompositions are widely used in computational linear algebra. Besides solving linear systems of equations, they are also used to compute other objects more specific to exact arithmetic: computing the rank, sampling a vector from the null- space, computing echelon forms and rank profiles. The row rank profile (resp. column rank profile) of an m×n matrix A with rank r, denoted by RowRP(A) (resp. ColRP(A)), is the lexicographicallysmallest sequence of r indices of linearly independent rows (resp. columns) of A. An m×n matrix has generic row(resp.column)rankprofileifitsrow(resp.column)rankprofileis(1,..,r).Lastly,an m×n matrix has generic rank profile if its r first leading principal minors are non-zero. Note that if a matrix has generic rank profile, then its row and column rank profiles are ⋆ This research was partly supported by the HPAC project of the French Agence Nationale de la Recherche(ANR11BS02013). Email addresses: [email protected] (Jean-GuillaumeDumas), [email protected] (Cl´ementPernet),mailto:[email protected] (ZiadSultan). URLs: http://www-ljk.imag.fr/membres/Jean-Guillaume.Dumas/ (Jean-GuillaumeDumas), http://lig-membres.imag.fr/pernet/ (Cl´ementPernet), http://moais.imag.fr/membres/ziad.sultan (ZiadSultan). 2 generic,buttheconverseisfalse:thematrix[01]doesnothavegenericrankprofileeven 10 if its row and column rank profiles are generic. The row support (resp. column support) of a matrix A, denoted by RowSupp(A) (resp. ColSupp(A)), is the subset of indices of its non-zero rows (resp. columns). We recall that the row echelon form of an m×n matrix A is an upper triangular matrixE =TA,for anon-singularmatrixT,withthe zerorowsofE atthe bottomand the non-zero rows in stair-case shape: min{j : a 6= 0} < min{j : a 6= 0}. As T is i,j i+1,j nonsingular,the columnrank profile of A is that of E,and therefore correspondsto the column indices of the leading elements in the staircase. Similarly the row rank profile of A is composed of the row indices of the leading elements in the staircase of the column echelon form of A. Rank profile and triangular matrix decompositions. The rank profiles of a matrix and the triangular matrix decomposition obtained by Gaussian elimination are strongly re- lated.Theeliminationofmatriceswitharbitraryrankprofilesgivesrisetoseveralmatrix factorizations and many algorithmic variants. In numerical linear algebra one often uses the PLUQ decomposition, with P and Q permutation matrices, L a lower unit trian- gular matrix and U an upper triangular matrix. The LSP and LQUP variants of [1] are used to reduce the complexity of rank deficient Gaussian elimination to that of ma- trix multiplication. Many other algorithmic decompositions exist allowing fraction free computations [2], in-place computations [3,4] or sub-cubic rank-sensitive time complex- ity [5,4]. The reader may refer to [4] for a detailed comparison between these matrix factorizations, and further details on the CUP (resp. PLE) variants, revealing the row (resp.column)rankprofiles.Allthesealgorithms,togetherwiththeschoolbookGaussian elimination algorithm share the property that, for a row rank profile computation, the pivotsearchprocessesrowsin order,andsearchesa pivotinallpossible columnposition before declaring the row linearly dependent with the previous ones. As a consequence, blockingislimitedtoonlyonedimension(inthiscasetherowdimension)leadingtoslab algorithms [6] operating on rectangular blocks of unbalanced dimensions. This reduces the data locality of the algorithm, and therefore penalizes the efficiency of implementa- tions in practice. In parallel, this blocking also puts more constrains on the dependency between tasks [7]. Contribution with respect to the state of the art. In [8] we proposed a first Gaussian elimination algorithm, with a recursive splitting of both row and column dimensions, whichsimultaneouslycomputestherowandcolumnrankprofilewhilepreservingthesub- cubic rank-sensitive time complexity and keeping the computation in-place. It showed that slab blocking is not a necessarycondition for a Gaussian elimination to revealrank profiles. Consequently, we have further analyzed the conditions on the pivoting that reveal the rank profiles in [9], where we introduced a new matrix invariant, the rank profile matrix. This normal form contains the row and column rank profile information of the matrix and that of all its leading sub-matrices. This normal form is closely related to a permutation matrix appearing in the Bruhat decomposition [10] and in related variants [11,12,13,14,15]. Still, no connection to the rankprofilesweremade.Inanothersetting,theconstructionofmatrixSchubertvarieties in [16, Ch. 15] defines a similar invariant, but presents neither a matrix decomposition nor any computational aspects. More precisely, the present paper gathers the key contributions of [8] and [9]: 3 • wedefineofanewmatrixinvariantoverafield,therankprofilematrix,summarizing allinformationontherowandcolumnrankprofilesofalltheleadingsub-matrices; • we study the conditions for a Gaussian elimination algorithm to compute all or part of this invariant, through the corresponding PLUQ decomposition; • asaconsequence,weshowthattheclassicaliterativeCUPdecompositionalgorithm can actually be adapted to compute the rank profile matrix. Used, in a Crout variant,asabase-casetoourimplementation,itdeliversasignificantimprovement in efficiency; • wealsoshowthatboththerowandthecolumnechelonformsofamatrixcanbere- coveredfromsomePLUQdecompositionsthanksto anelementarypost-processing algorithm. Further, we here develop three novel aspects: • we extend the notion of rank profile matrix over arbitrary rings; • wemakefurtherconnectionswithexistingmatrixdecompositions,inparticularthe Bruhat and generalized Bruhat decompositions; • lastly, we extend the recent algorithmic improvements of [17,18,19] for low rank matrices:indeed,weshowherethatthealgorithmin[18]computestherankprofile matrix; we also propose an algorithmic variant reducing the leading constant by a factorofthree;andweproposeanalgorithmicreductioncomputingtherankprofile matrix in time bounded by (rω +mn)1+o(1). Organization of the article. We first introduce in Section 2 the rank profile matrix R , A a normal form summarizing all rank profile information of a matrix and of all its lead- ing sub-matrices. This definition is extended over arbitrary rings in Section 2.2. There, two notions of rank are considered: the spanning rank and McCoy’s rank and we show that only the latter is suited to the definition of the rank profile matrix. We then study the conditions that PLUQ decomposition algorithms must satisfy in order to reveal the rank profile structure of a matrix. For this, we decompose, in Section 3, the pivoting strategy of any PLUQ algorithm into two types of operations: the search of the pivot and the permutation used to move it to the main diagonal. We propose new search and new permutation strategies and show what rank profiles are computed using any possi- ble combination of these operations. In particular we show three new pivoting strategy combinations that compute the rank profile matrix. As an illustration, we show in Section 4 how these pivoting strategies instantiate in iterativeorrecursivealgorithms,usingslabortileblocking.Connectionsaremadetothe mostcommoneliminationalgorithmsandwestateinfull detailsthe recenttile recursive algorithm of [8], implementing the new pivoting strategy. Section 5 shows how this better understanding of the pivoting conditions helped to design an iterative Crout CUP with rotations, to be used as a base case for the tile recursivealgorithm.UsingtheanalysisofSection3,itisshownthatitcomputestherank profile matrix while improving the previous base case implementation by a significant speed-up. We then show in Section 6 how a PLUQ decomposition revealing the rank profile matrix relates with other triangular matrix decompositions, such as the LEU and the classical,modifiedorgeneralizedBruhatdecompositions,orthecomputationofbothrow and column echelon forms, from a single PLUQ decomposition. Lastly, we focus in section 7 on the recent algorithmic improvements for the com- putation of the row or column rank profile for matrices with low rank. We show that 4 the algorithm in [18] actually computes the rank profile matrix, but with a classical complexity. Then we show how to adapt its fast improvement [19] to also compute the rank profile matrix, before discussing algorithmic variants that improve on the leading constant factor in the latter algorithm’s complexity bound. Notations. In the following,0 denotes the m×n zero matrix.For two list of indices m×n P and Q, A denotes the sub-matrix of A formed by the rows of index in P and the P,Q columnsofindex inQ.Inparticular,A denotesthe contiguousblockofcoefficients i..j,k..l in A of rows position between i and j and columns position between k and l. We may denote by ∗ the sequence of all possible row or column indices: e.g. A denotes the i-th i,∗ row of A. To a permutation σ : {1,...,n} → {1,...,n} we define the associated permutation matrixP(σ),permutingrowsbyleftmultiplication:therowsofP(σ)AarethatofAper- mutedbyσ.Reciprocally,forapermutationmatrixP,wedenotebyσ(P)theassociated permutation. 2. The rank profile matrix We propose in Theorem 3 the definition of the rank profile matrix, an invariant sum- marizing all information on the rank profiles of a matrix. As will be discussed in this sectionandinsection6,thisinvariantiscloselyrelatedtotheBruhatdecomposition[10] and its generalizations [12,20,16]. 2.1. Definition over a field We consider first matrices over a field K and a valid pivot in Gaussian elimination is any non-zero element of K. Definition 1. An r-sub-permutation matrix is a matrix of rank r with only r non-zero entries equal to one. Lemma 2. An m×n r-sub-permutation matrix has at most one non-zero entry per row andpercolumn,andcanbewrittenP Ir QwhereP andQarepermutation 0(m−r)×(n−r) matrices. h i Theorem 3. Let A∈Km×n of rank r. There exists a unique m×n r-sub-permutation matrix R of which every leading sub-matrix has the same rank as the corresponding A leading sub-matrix of A. This sub-permutation matrix is called the rank profile matrix of A. Proof. Weproveexistencebyinductionontherowdimensionoftheleadingsubmatrices. IfA =0 ,setting R(1) =0 satisfiesthe defining condition.Otherwise,letj 1,1..n 1×n 1×n be theindexofthe invertibleelementinA andsetR(1) =eT the j-thn-dimensional 1,1..n j row canonical vector, which satisfies the defining condition. Nowforagiveni∈{1,...,m},supposethatthereisauniquei×nrankprofilematrix R(i)suchthatrank(A )=rank(R )foreveryj ∈{1..n}.Ifrank(A )= 1..i,1..j 1..i,1..j 1..i+1,1..n rank(A ), then R(i+1) = R(i) . Otherwise,consider k, the smallest column index 1..i,1..n 01×n h i 5 such that rank(A ) = rank(A )+1 and set R(i+1) = R(i) . Any leading 1..i+1,1..k 1..i,1..k eT k sub-matrix of R(i+1) has the same rank as the corresponding leadihng suib-matrix of A: first,foranyleadingsubsetofrowsandcolumnswithlessthanirows,thecaseiscovered B u by the induction; second define = A , where u,v are vectors and x is a 1..i+1,1..k vT x scalar.Fromthedefinitionofk,vislinearlydependentwithB andthusanyleadingsub- matrix of B has the same rank as the corresponding sub-matrix of R(i+1). Similarly, vT fromthedefinitionofk,thesamereasoningworkswhenconsideringmorethankcolumns, (cid:2) (cid:3) with a rank increment by 1. Lastly we show that R(i+1) is an r -sub-permutation matrix. Indeed, u is linearly i+1 dependent with the columns of B: otherwise, rank( B u ) = rank(B)+ 1. From the definition of k we would then have rank( B u ) = rhank(iB u )+1 = rank(B)+2 = vT x rank( vBT )+2, but this is a contradiction(cid:2). Con(cid:3)sequently,hthe ki-th column of R(i) is all zero, and R(i+1) is an r -sub-permutation matrix. (cid:2) (cid:3) i+1 To prove uniqueness, suppose there exist two distinct rank profile matrices R(1) and R(2) for a given matrix A and let (i,j) be the lexicographically minimal coordinates where R(1) 6= R(2). As a consequence the rank of the (i,j)-leading submatrices of R(1) i,j i,j and R(2) differ but should both be equal to rank(A ), a contradiction. ✷ 1..i,1..j 2030 1000 Example 4. A= 1000 has R = 0010 for rank profile matrix over Q. 0040 A 0000 (cid:20)0201(cid:21) (cid:20)0100(cid:21) Remark 5. The permutation matrix introduced in the modified Bruhat decomposition of [20], and defined there only for invertible matrices, is also the matrix E introduced in Malaschonok’sLEUdecomposition[14, Theorem1].In the latter paper,analgorithm for this decomposition was only shown over a field for m = n = 2k, and no connection wasmadetothe relationwithranksandrankprofiles.We haveshownin[8,Corollary1] that E is in fact the rank profile matrix and we made the explicit connection in in [8, Corollary 1], as recalled in Section 6. We here generalize the existence to arbitrary rank r and dimensions m and n and after proving its uniqueness, we propose this definition as a new matrix normal form. We next also generalize the construction over a ring in Section 2.2. The rank profile matrix has the following properties: Lemma 6. Let A be a matrix. (1) R is diagonal if A has generic rank profile. A (2) R is a permutation matrix if A is invertible A (3) RowRP(A)=RowSupp(R ); ColRP(A)=ColSupp(R ). A A Moreover, for all 1≤i≤m and 1≤j ≤n, we have: (4) RowRP(A )=RowSupp((R ) ) 1..i,1..j A 1..i,1..j (5) ColRP(A )=ColSupp((R ) ), 1..i,1..j A 1..i,1..j These properties show how to recover the row and column rank profiles of A and of any of its leading sub-matrix. 6 2.2. Generalization over a ring Over a ring, several notions of rank exist: the spanning rank, McCoy’s rank, Smith’s rank, etc. Definition 7. Let A be an m×n matrix over a ring R: • McCoy’s rank, denoted by r , is the largest r such that no size s minor of A M with s > r is a unit. Equivalently, it is also the largest r such that there exists a selection of r columns whose right nullspace is the zero vector [21, Theorem 51]. Over a principal ideal ring (PIR), this is the number of invertible determinantal divisors. • The spanning rank, denoted by s , the smallest r such that A = BC, where B r is m×r and C is r×n. Over a principal ideal ring (PIR), this is equivalent to Smith’srank[22,Definiton2.2],thenumberofnon-zerodeterminantaldivisors[23, Proposition 1.1.(d)]. First, the rank profile matrix over a field gives rise to a rank profile matrix over any principal ideal ring. Corollary 8. Let D be a principal ideal domain (PID) and let A∈Dm×n with McCoy’s rank r. There exists a unique m×n r-sub-permutation matrix R of which every leading A sub-matrix has the same rank as the corresponding leading sub-matrix of A. This sub- permutation matrix is called the rank profile matrix of A. Proof. From[23,Proposition1.6],overaPIDwithfieldoffractionsK,s =r =rank . r M K Thus R over K satisfies the requirements over D and is the unique such matrix. ✷ A However, when the ring R contains zero divisors, the notion of spanning rank differs from that of McCoy’s rank. Lemma 9 shows that the notion of rank profile matrix can not be defined from the spanning rank. Lemma 9. Over Z/4Z, no matrix with exactly 1 or 2 non-zero entries is such that its leading sub-matrices have the same spanning ranks as that of A=[02]. 21 Proof. The spanning rank of A is 1, as A=[2][21]. However, note that the spanning 1 rank of the leading 1×1 sub-matrix of A is 0, while that of the leading 1×2 and the leading 2×1 sub-matrices of A are equal to one. Hence, if a matrix R = ab is such c d that any of its leading sub-matrices has the same spanning rank as the corresponding (cid:2) (cid:3) sub-matrix of A, then a=0 and b,c6=0. Now as R must have spanning rank one, there must existx,y,z,t,such that[xy][z t]= 0c db . As b,c6=0,then x,y,z,tareall nonzero and since a = 0, we have x = z = 2, but then, z,t ∈ {−1,1} otherwise b = 0 or c = 0. Consequently d6=0. ✷ (cid:2) (cid:3) WithMcCoy’srank,onthecontrary,therankprofilematrixofthematrixinLemma9 could be defined as the matrix [00]. Indeed, remark that the McCoy’s rank of the sub- 01 matrix [02] is 0. Hence, we focus in the remaining of this section on McCoy’s rank and only assume that R is a commutative ring R with a unit element 1. We will need the following lemmata. 7 Lemma 10. Consider a matrix B ∈ Rm×n, and two vectors u and v. Then δ = 1 B rM B u −rM(B) and δ2 =rM −rM(B) satisfy δi ∈{0,1}. vT (cid:16)h i(cid:17)   Proof. AppendingacolumntoBcannotreducethemaximalnumberofcolumnshaving arightnull-spacereducedtothezerovector,henceδ ≥0.Letr =r (B).Ifr=n,then 1 M rM([B u]) ≤ n+1 and δ ≤ 1. Otherwise, this means that any subset of r+1 columns ofB vanishesinalinearcombinationwithsomenon-zerocoefficients.Consequently,any subset of r +2 columns of [B u], which necessarily includes at least r +1 columns of B, vanishes in a linear combination with some non-zero coefficients, thus showing that δ ≤1. Applying the same reasoning on B T, proves that δ ∈{0,1}. ✷ 1 vT 2 (cid:2) (cid:3) Lemma 11. Consider an m×n matrix B, two vectors, u and v, and a scalar α. Denote by δ ,δ ,δ ,δ the discrepancies in rank of the following matrices: 1 2 3 4 (1) rM(B)+δ1 =rM B u δ (cid:16)h i(cid:17) 4 B (2) r (B)+δ =r M 2 M vT δ1 B u   B u (3) rM B u +δ3 =rM δ2 vT α (cid:16)h i(cid:17) δ3 vT α   B B u (4) r +δ =r M 4 M vT vT α Then:     (i) (δ =0 and δ =1) is equivalent to (δ =0 and δ =1). 1 4 2 3 (ii) (δ =1 and δ =1) is equivalent to (δ =1 and δ =1). 1 3 2 4 (iii) (δ =1 and δ =0) implies (δ =0 and δ =1). 1 3 2 4 (iv) (δ =1) implies (δ =1). 1 4 (v) (δ =1) implies (δ =1). 2 3 Proof. Let r = r (B). From Lemma 10, all the discrepancies are either zero or one, M since the proposed ranks involve matrices with only one extra row or column. From the definitions of δ and δ , on the one hand, those of δ and δ on the other hand, we have 1 3 2 4 that r+δ +δ =r+δ +δ . (1) 2 4 1 3 Thus: (i) if (δ = 0 and δ = 1), then δ +1 = δ and thus the only possibility is (δ = 1 4 2 3 2 0 and δ = 1). Reciprocally, δ = δ + 1 yields (δ = 0 and δ = 1) as only 3 4 1 1 4 possibility. (ii) Similarly there is only way for the rank to be augmented by 2. (iii) Suppose that δ1 = 1 and δ3 = 0. Then rM([B u]) = rM( vBT αu ) = r+1. There existsubsetsofr+1indices I andJ suchthatthe sub-matrix[B u] hasaunit (cid:2) (cid:3) I,J determinant.Sinceδ =1thennecessarilyn+1∈J,otherwiser (B)=r+1.As 1 M 8 r (B ) = r, there is a k ∈ I, such that X = B has a unit M I,J\{n+1} I\{k},J\{n+1} determinant and is therefore invertible over R. Consider the permutation matrices P and Q such that PBQ = [X Y ]. Let b1 = Pu where b has dimension r. We Z T b2 1 have (cid:2) (cid:3) X−1 −X−1Y I 0 r PBQ = .     I M N n−r From[24,4.11.(c)],weknowthatmultiplicationbymatriceswithunitdeterminant preservesr .Note thatN does notcontaina unit, otherwiser (B)=r+1.Now M M we also have X−1 −X−1Y −X−1b 1 Q I 0 0 P B u  I 0 = r ,   n−r   1 M N β h i    1         where β =−Mb +b is a vector containing at least one unit u at index i and 1 2 β X−1 −X−1Y −X−1b I 0 0 1 r P B u Q  I 0 =M N β=C.  1vT α 1 n−r  1  wT xT γ            If xT would contain a unit at index j , we could form the subsets of indices x I = {1,...,r,i ,m+1} and J = {1,...,r,j ,n+1} such that C has a unit β x I,J determinant,showingthatr ( B u )≥r+2whichis acontradiction.Therefore, M vT α since (cid:2) (cid:3) I 0 r P B X−1 −X−1Y Q =M N,  1vT  I  n−r wT xT        we deduce that r ( B )=r. Hence, δ =0, and from (1), we have δ =1. M vT 2 4 (iv) Combining (ii) and (iii) leads directly to (iv). (cid:2) (cid:3) (v) Applying (iv) on the transpose of A yields (v). ✷ Lemma 11, point (iii), is typically what is not true with the spanning rank: on the counter-example matrix A =[02] of Lemma 9, we have δ =1, δ =0, but δ =1 and 1 21 1 3 2 δ = 0. We are now ready to extend the proof of existence of the rank profile matrix to 4 matrices over a ring, using the notion of McCoy’s rank. Definition 12. Anr-sub-permutation matrixisamatrixofMcCoyrankr withonlyr M non-zero entries equal to one. Lemma 13. Anm×nr-sub-permutationmatrix has at mostonenon-zeroentryperrow andpercolumn,andcanbewrittenP Ir QwhereP andQarepermutation 0(m−r)×(n−r) matrices. h i 9 Proof. Let Abe anm×n r -sub-permutationmatrix.Let µand ν be respectively the M numberofnon-zerocolumnsandrowsinA.AsAhasonlyr =r (A)non-zeroelements, M we have that s = min(µ,ν) ≤ r. If s = µ, then A writes A = B I 0 P where B is s the m×s matrix formed by the s non-zero columns of A and P a hpermiutation matrix. From [24, Lemma 4.14], we have that rM(A) ≤ min{rM(B),s = rM Is 0 ,rM(P)}. This decomposition thus shows that the rank of A is at most s and t(cid:16)hhe onliy(cid:17)way for A tohavearankr istohaveonenonzeroelementpercolumn.Thesamereasoningapplies on the non-zero rows when s=ν. ✷ Theorem 14. Let A ∈ Rm×n of McCoy rank r. There exist a unique m ×n r-sub- permutation matrix R of which every leading sub-matrix has the same McCoy rank A as the corresponding leading sub-matrix of A. This sub-permutation matrix is called the (McCoy) rank profile matrix of A. Proof. Weproveexistencebyinductionontherowdimensionoftheleadingsubmatrices. Consider the first row of A. If it does not contain any unit, then setting R(1) =0 1×n satisfiesthedefiningcondition.Otherwise,letj betheindexoftheleftmostunitelement in A and set R(1) =eT the j-th n-dimensional canonical row vector, which satisfies 1,1..n j the defining condition. Now suppose that there exists a rank profile matrix R(s) of which every leading sub- matrixhas the same McCoyrankasthe correspondingleading sub-matrixofA .Let 1..s,∗ r =r (A ). Then, from Lemma 10, either r =r or r =r +1: s M 1..s,∗ s+1 s s+1 s (i) Ontheonehand,ifr =r ,letj bethelargestindexsuchthatr (A )< s+1 s 0 M 1..s,1..j0 rs. We then use Lemma 11, point (iii), with B = A1..s,1..j0, [B u] = A1..s,1..(j0+1), B = A , B u = A : we thus have that r (A ) = vT 1..s+1,1..j0 vT α 1..s+1,1..(j0+1) M 1..s,1..j0 r (A ) = r −1 and r (A ) = r (A ) = r (cid:2)M (cid:3) 1..(s+1),1..j0 (cid:2)s (cid:3) M 1..s,1..j0+1+k M 1..s+1,1..j0+1+k S for all k =0..(n−j −1). 0 Let now j be the largest index such that r (A ) < (r −1), the same 1 M 1..s,1..j1 s reasoning shows that the McCoy ranks of the leading matrices with j +1 to j 1 0 columns ands or s+1rowsareidentical.Continue with j until there areno more i R(s) columns to consider. This shows that is a suitable rank profile matrix for   0 A1..(s+1),∗.   (ii) On the other hand, if r = r +1, it means in particular, that r (A ) 6= s+1 s M 1..s,1..n r (A ). We thus againconsider the leading matrices with s or s+1 rows. M 1..s+1,1..n Letnowj bethesmallestindexsuchthatr (A )6=r (A ).Asthere M 1..s,1..j M 1..s+1,1..j is only one extra row in the latter, we must have exactly r (A )+1=r (A ) (2) M 1..s,1..j M 1..s+1,1..j R(s) We show now that if e is the j-th canonical vector then T = is a suitable j  eT  j rankprofilematrixforA1..s+1,∗:weuse Lemma11,point(i),withB=A1..s,1..j−1, [B u] = A1..s,1..j, vBT = A1..s+1,1..j−1, vBT αu = A1..s+1,1..j. From the definition of j in Equation ((cid:2)2), (cid:3)we have that rM([(cid:2)B u])(cid:3)6= rM vBT αu and therefore that (cid:0)(cid:2) (cid:3)(cid:1) 10