Fast Computation of the Rank Profile Matrix and the Generalized Bruhat Decomposition

01/08/2016
by   Jean-Guillaume Dumas, et al.
0

The row (resp. column) rank profile of a matrix describes the stair-case shape of its row (resp. column) echelon form.We here propose a new matrix invariant, the rank profile matrix, summarizing all information on the row andcolumn rank profiles of all the leading sub-matrices.We show that this normal form exists and is unique over any ring, provided that the notion of McCoy's rank is used, in the presence of zero divisors.We then explore the conditions for a Gaussian elimination algorithm to compute all or part of this invariant,through the corresponding PLUQ decomposition. This enlarges the set of known Elimination variants that compute row or column rank profiles.As a consequence a new Crout base case variant significantly improves the practical efficiency of previously known implementations over a finite field.With matrices of very small rank, we also generalize the techniques ofStorjohann and Yang to the computation of the rank profile matrix, achieving an (r^ω+mn)^1+o(1) time complexity for an m × n matrix of rank r, where ω is the exponent of matrix multiplication. Finally, by give connections to the Bruhat decomposition, and severalof its variants and generalizations. Thus, our algorithmicimprovements for the PLUQ factorization, and their implementations,directly apply to these decompositions.In particular, we show how a PLUQ decomposition revealing the rankprofile matrix also reveals both a row and a column echelon form ofthe input matrix or of any of its leading sub-matrices, by a simplepost-processing made of row and column permutations.

READ FULL TEXT VIEW PDF

page 1

page 2

page 3

page 4

02/18/2022

Rank-Sensitive Computation of the Rank Profile of a Polynomial Matrix

Consider a matrix 𝐅∈𝕂[x]^m × n of univariate polynomials over a field 𝕂....
09/11/2019

Elimination-based certificates for triangular equivalence and rank profiles

In this paper, we give novel certificates for triangular equivalence and...
02/13/2017

Certificates for triangular equivalence and rank profiles

In this paper, we give novel certificates for triangular equivalence and...
07/05/2022

Numerical considerations and a new implementation for ICS

Invariant Coordinate Selection (ICS) is a multivariate data transformati...
02/26/2018

Symmetric indefinite triangular factorization revealing the rank profile matrix

We present a novel recursive algorithm for reducing a symmetric matrix t...
08/21/2019

Perturbations of CUR Decompositions

The CUR decomposition is a factorization of a low-rank matrix obtained b...
03/11/2016

Matrix factoring by fraction-free reduction

We consider exact matrix decomposition by Gauss-Bareiss reduction. We in...

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 matrix with rank , denoted by RowRP(A) (resp. ColRP(A)), is the lexicographically smallest sequence of indices of linearly independent rows (resp. columns) of . An matrix has generic row (resp. column) rank profile if its row (resp. column) rank profile is . Lastly, an matrix has generic rank profile if its first leading principal minors are nonzero. Note that if a matrix has generic rank profile, then its row and column rank profiles are generic, but the converse is false: the matrix does not have generic rank profile even if its row and column rank profiles are generic. The row support (resp. column support) of a matrix , denoted by (resp. ), is the subset of indices of its nonzero rows (resp. columns).

We recall that the row echelon form of an

matrix is an upper triangular matrix , for a nonsingular matrix , with the zero rows of at the bottom and the nonzero rows in stair-case shape: . As is nonsingular, the column rank profile of is that of , and therefore corresponds to the column indices of the leading elements in the staircase. Similarly the row rank profile of is composed of the row indices of the leading elements in the staircase of the column echelon form of .

Rank profiles and triangular matrix decompositions.

The rank profiles of a matrix and the triangular matrix decompositions obtained by Gaussian elimination are strongly related. The elimination of matrices with arbitrary rank profiles gives rise to several matrix factorizations and many algorithmic variants. In numerical linear algebra one often uses the PLUQ decomposition, with and permutation matrices, a lower unit triangular matrix and an upper triangular matrix. The LSP and LQUP variants of Ibarra et al. (1982) have been introduced to reduce the complexity of rank deficient Gaussian elimination to that of matrix multiplication. Many other algorithmic decompositions exist allowing fraction free computations Jeffrey (2010), in-place computations Dumas et al. (2008); Jeannerod et al. (2013) or sub-cubic rank-sensitive time complexity Storjohann (2000); Jeannerod et al. (2013). The reader may refer to Jeannerod et al. (2013) for a detailed comparison between these matrix factorizations, and further details on the CUP (resp. PLE) variants, revealing the row (resp. column) rank profiles. All these algorithms, together with the schoolbook Gaussian elimination algorithm share the property that, for a row rank profile computation, the pivot search processes rows in order, and searches a pivot in all possible column position before declaring the row linearly dependent with the previous ones. As a consequence, blocking is limited to only one dimension (in this case the row dimension) leading to slab algorithms Klimkowski and van de Geijn (1995) operating on rectangular blocks of unbalanced dimensions. This reduces the data locality of the algorithm, and therefore penalizes the efficiency of implementations in practice. In parallel, this blocking also puts more constrains on the dependencies between tasks Dumas et al. (2015a).

Contribution with respect to the state of the art.

In Dumas et al. (2013) we proposed a first Gaussian elimination algorithm, with a recursive splitting of both row and column dimensions, which simultaneously computes the row and column rank profile while preserving the sub-cubic rank-sensitive time complexity and keeping the computation in-place. It showed that slab blocking is not a necessary condition for a Gaussian elimination to reveal rank profiles. Consequently, we have further analyzed the conditions on the pivoting that reveal the rank profiles in Dumas et al. (2015b), 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 Bruhat (1956) and in related variants Della Dora (1973); Grigoriev (1981); Bourbaki (2008); Malaschonok (2010); Manthey and Helmke (2007). Still, none of these did connect it to the notion of rank profile. In another setting, the construction of matrix Schubert varieties in (Miller and Sturmfels, 2005, 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 Dumas et al. (2013) and Dumas et al. (2015b):

  1. we define a new matrix invariant over a field, the rank profile matrix, summarizing all information on the row and column rank profiles of all the leading sub-matrices;

  2. we study the conditions for a Gaussian elimination algorithm to compute all or part of this invariant, through the corresponding PLUQ decomposition;

  3. as a consequence, we show that the classical iterative CUP decomposition algorithm can actually be adapted to compute the rank profile matrix. Used, in a Crout variant, as a base-case to our recursive implementation over a finite field, it delivers a significant improvement in efficiency;

  4. we also show that both the row and the column echelon forms of a matrix can be recovered from some PLUQ decompositions thanks to an elementary post-processing algorithm.

Further, we develop three novel aspects:

  1. we study how the notion of rank profile matrix can be generalized over an arbitrary ring. We show that none of the common definition of rank over an arbitrary commutative ring allow to define a rank profile matrix in general. However, over a principal ideal domain and over a finite chain ring, we can produce such a definition.

  2. we make further connections with existing matrix decompositions, in particular the Bruhat and the generalized Bruhat decompositions, exhibiting both the row and the column echelon in a single decomposition;

  3. lastly, we extend the recent algorithmic improvements of Cheung et al. (2013); Storjohann and Yang (2014, 2015) for low rank matrices: indeed, we show here that the algorithm in Storjohann and Yang (2014) computes the rank profile matrix; propose an algorithmic variant thereof, reducing the leading constant by a factor of three; and lastly show an algorithmic reduction computing the rank profile matrix in time bounded by with a Las-Vegas probabilistic algorithm.

Organization of the article.

We first introduce in Section 2 the rank profile matrix , and study in Section 2.2 its generalization over an arbitrary ring. We then study in Section 3 under which condition a PLUQ decomposition algorithms reveals the rank profile structure of a matrix. For this, we investigate existing and new pivoting strategies, based on all combination of elementary search and permutation operations, showing for each of them what part of the rank profile information is being computed. In particular we propose three new pivoting strategies that compute the rank profile matrix. As an illustration, we show in Section 4 how these pivoting strategies instantiate in iterative or recursive algorithms, using slab or tile blocking. Connections are made to the most common elimination algorithms and we state in full details the recent tile recursive algorithm of Dumas et al. (2013), implementing one of the new pivoting strategy. Section 5 shows how this better understanding on the pivoting strategies has resulted in the design of an iterative Crout CUP decomposition with rotations, to be used as a base case for the tile recursive algorithm, speeding up the computation efficiency, while still recovering the rank profile matrix invariant. 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, modified or generalized Bruhat decompositions, or the computation of both row and column echelon forms, from a single PLUQ decomposition. Lastly, we extend in Section 7 the recent algorithms of Storjohann and Yang (2014, 2015) for the row or column rank profile of matrices with low rank to computed the rank profile matrix within the same complexities.

Notations.

In the following, denotes the zero matrix. For two list of indices and , denotes the sub-matrix of formed by the rows of index in and the columns of index in . In particular, denotes the contiguous block of coefficients in of rows position between and and columns position between and . We may denote by the sequence of all possible row or column indices: e.g. denotes the -th row of . To a permutation we define the associated permutation matrix , permuting rows by left multiplication: the rows of are that of permuted by . Reciprocally, for a permutation matrix , we denote by the associated permutation.

2 The rank profile matrix

We propose in Theorem 3 the definition of the rank profile matrix, an invariant summarizing all information on the rank profiles of a matrix. As will be discussed in this section and in Section 6, this invariant is closely related to the Bruhat decomposition Bruhat (1956) and its generalizations Grigoriev (1981); Tyrtyshnikov (1997); Miller and Sturmfels (2005).

2.1 Definition over a field

We first consider matrices over an arbitrary commutative field .

Definition 1.

An -sub-permutation matrix is a matrix of rank with only non-zero entries equal to one.

Lemma 2.

An -sub-permutation matrix has at most one non-zero entry per row and per column, and can be written where and are permutation matrices.

Theorem 3.

Let of rank . There exists a unique -sub-permutation matrix of which every leading sub-matrix has the same rank as the corresponding leading sub-matrix of . This sub-permutation matrix is called the rank profile matrix of .

Proof..

We prove existence by induction on the row dimension. of the leading submatrices.

If , setting satisfies the defining condition. Otherwise, let be the index of the first invertible element in and set the j-th -dimensional canonical row vector, which satisfies the defining condition.

Now for a given , suppose that there is a unique rank profile matrix such that for every and . If , then . Otherwise, consider , the smallest column index such that and set . Define , where are vectors and is a scalar. By definition of , we have .

First we show that is an -sub-permutation matrix. If not, then the -th column of would contain a which, by induction hypothesis, would imply that . Hence we would have , a contradiction.

Then we show that any leading sub-matrix of with has the same rank as the corresponding leading sub-matrix of . The case is covered by the induction; second, since , any leading sub-matrix of has the same rank as the corresponding sub-matrix of , which covers the case . Lastly, for and , the definition of implies .

To prove uniqueness, suppose there exist two distinct rank profile matrices and for a given matrix and let be the lexicographically minimal coordinates where . The rank of the -leading submatrices of and differ but should both be equal to , a contradiction. ∎

Example 4.

has for rank profile matrix over .

Remark 5.

The permutation matrix introduced in the modified Bruhat decomposition of Tyrtyshnikov (1997), and defined there only for invertible matrices, is also the matrix introduced in Malaschonok’s LEU decomposition (Malaschonok, 2010, Theorem 1). In the latter paper, an algorithm for this decomposition was only shown over a field for , and no connection was made to the relation with ranks and rank profiles. We have shown in (Dumas et al., 2013, Corollary 1) that is in fact the rank profile matrix and made the connection to the PLUQ decomposition explicit, as recalled in Section 6. We here generalize the existence to arbitrary rank and dimensions and and after proving its uniqueness, we propose this definition as a new matrix normal form.

The rank profile matrix has the following properties:

Lemma 6.

Let be a matrix.

  1. is diagonal if and only if has generic rank profile.

  2. is a permutation matrix if and only if is invertible

  3. ; .

Moreover, for all and , we have:

  1. ,

These properties show how to recover the row and column rank profiles of and of any of its leading sub-matrix.

2.2 Generalization over a commutative ring with unity

We now explore if and under which condition, the notion of rank profile matrix can be generalized over an arbitrary commutative ring with unity. As the rank profile matrix over a field relies on the notion of rank, which generalization over an arbitrary ring leads to several definitions, we review the most common of these definitions and explore their capacity to define a rank profile matrix invariant.

Let be a commutative ring with unity. Following (Brown, 1992, §4), we will denote by the -th determinantal ideal, namely the ideal in generated by all minors of , and set by definition . We have . We will also denote by the annihilator ideal of a set .

Definition 7.

For an matrix over a communtative ring , the rank of can be defined as:

  1. , called spanning rank by Brown (1998) or Schein rank by (Bhaskara Rao, 2002, §2.4),

  2. introduced in (McCoy, 1948, Theorem 51),

  3. (Bhaskara Rao, 2002, §2.4) (this is also defined in (Brown, 1992, §4, Ex. 11)),

  4. used in (Norton and Salagean, 2000, Def. 2.5).

We illustrate the value of these ranks in Table 1 for three matrices (columns 3, 6, 9) and their leading and sub-matrices (columns 2, 5, 8). These ranks will be used to argue whether a rank profile matrix can be defined in each case (columns 4, 7 and 10 if it exists).

over over
, , ,
1 1 none 1 2 1 2
0 1 1 1 none 1 2
1 1 none 1 2 1 2
0 0 0 0 1,0 2 none
Table 1: Ranks of matrices according to the various definitions. These three matrices are counterexamples showing that none of the four rank definitions of Definition 7 can be used to define a rank profile matrix in general.

For instance since over . This matrix shows that can not be used to define a rank profile matrix. Indeed, if it would exist, its first row would be (as the leading submatrix has rank 0 and the leading submatrix has rank 1). Similarly the first column of this rank profile matrix would be . But the rank of the permutation matrix would be . We will use the same reasoning pattern, to show that can not define a rank profile matrix for , and neither does for over : in these two cases, the rank profile matrix should be to satisfy the rank condition on the and sub-matrices, but the rank is only 1. Laslty, and , therefore the rank profile matrix of should be of the form , but it can then never be of rank 2, as the matrix .

Remark also that the rank profile matrix invariant is strongly connected with elimination as will be presented in the next sections. It therefore needs to be based on a notion of rank that is stable with multiplication by invertible matrices. This is the case with McCoy’s rank (Brown, 1992, 4.11.(c)), but not with . Indeed the rank of is 1 whereas the rank of is 0.

Consequently there is no notion of rank over an arbitrary commutative ring with unity supporting the definition of a rank profile matrix. We will now show that with additional assumptions on the ring, some of the above definitions for the rank coincide and make the existence of the rank profile matrix possible.

2.2.1 Over a principal ideal domain

Over a principal ideal domain, the existence of an underlying field of fractions guarantees that the McCoy rank, the spanning rank and the rank over the field of fractions coincide.

Corollary 8.

Let D be a principal ideal domain (PID) and let with McCoy’s rank . There exists a unique -sub-permutation matrix of which every leading sub-matrix has the same rank as the corresponding leading sub-matrix of . This sub-permutation matrix is called the rank profile matrix of .

Proof..

From (Brown, 1998, Proposition 1.6), over a PID with field of fractions , . Thus over satisfies the requirements over D and is the unique such matrix. ∎

2.2.2 Over a finite chain ring

A finite chain ring is a finite commutative ring with identity which ideals are ordered by inclusion (Clark and Liang, 1973). Equivalently, it is a finite local ring which maximal ideal is principal. These rings include and for prime, as well as all Galois rings of characteristic .

Over a finite chain ring, McCoy rank and the UnitMinor rank coincide, and therefore allow to define a rank profile matrix, as shown next.

Lemma 9.

Over a finite chain ring, McCoy’s rank is the largest positive integer such that there exist a unit minor of A.

Proof..

Krull’s theorem (Krull, 1938, Theorem 2) (see also (Hungerford, 1968, Proposition 4)) states that in a local ring, an element is a unit if and only if it is not in the maximal ideal. As the maximal ideal of the ring is principal, let be a generator thereof. Then is nilpotent. Indeed the ring is finite so there exists indices such that . Therefore . But as is not a unit, is not a unit either, but then must be a unit (otherwise the maximal ideal would contain and would not be proper). Therefore . The nilpotency index of is then the smallest positive integer such that .

Let and suppose that all minors are non-unit. They belong to the maximal ideal and for some . Therefore, , a contradiction. ∎

Corollary 10 ((Norton and Salagean, 2000, Corollary 2.7)).

Over a finite chain ring , with maximal ideal , McCoy’s rank is the rank over the field .

Corollary 11.

Let be a finite chain ring and let with McCoy’s rank . There exists a unique -sub-permutation matrix of which every leading sub-matrix has the same McCoy’s rank as the corresponding leading sub-matrix of . This sub-permutation matrix is called the rank profile matrix of .

Proof..

From Corollary 10, we consider the residue field . Then over satisfies the requirements over and is the unique such matrix. ∎

3 When does a PLUQ algorithm reveal the rank profile matrix?

From now on, for the sake of simplicity, we consider algorithms over a field.

3.1 Ingredients of a PLUQ decomposition algorithm

Over a field, the LU decomposition generalizes to matrices with arbitrary rank profiles, using row and column permutations (in some cases such as the CUP, or LSP decompositions, the row permutation is embedded in the structure of the or matrices). However such PLUQ decompositions are not unique and not all of them will necessarily reveal rank profiles and echelon forms. We will characterize the conditions for a PLUQ decomposition algorithm to reveal the row or column rank profile or the rank profile matrix.

We consider the four types of operation of a Gaussian elimination algorithm in the processing of the -th pivot:

Pivot search:

finding an element to be used as a pivot,

Pivot permutation:

moving the pivot in diagonal position by column and/or row permutations,

Update:

applying the elimination at position : ,

Normalization:

dividing the -th row (resp. column) by the pivot.

Choosing how each of these operation is done, and when they are scheduled results in an elimination algorithm. Conversely, any Gaussian elimination algorithm computing a PLUQ decomposition can be viewed as a set of specializations of each of these operations together with a scheduling.

The choice of doing the normalization on rows or columns only determines which of or will be unit triangular. The scheduling of the updates vary depending on the type of algorithm used: iterative, recursive, slab or tiled block splitting, with right-looking, left-looking or Crout variants (see Dongarra et al. (1998)). Neither the normalization nor the update impact the capacity to reveal rank profiles and we will thus now focus on the pivot search and permutation.

Choosing a search and a permutation strategy sets the matrices and of the PLUQ decomposition obtained and, as we will see, determines the ability to recover information on the rank profiles. Once these matrices are fixed, the and the factors are unique. We therefore introduce the pivoting matrix.

Definition 12.

The pivoting matrix of a PLUQ decomposition of rank is the -sub-permutation matrix

The nonzero elements of are located at the initial positions of the pivots in the matrix . Thus summarizes the choices made in the search and permutation operations.

3.1.1 Pivot search

The search operation vastly differs depending on the field of application. In numerical dense linear algebra, numerical stability is the main criterion for the selection of the pivot. In sparse linear algebra, the pivot is chosen so as to reduce the fill-in produced by the update operation. In order to reveal some information on the rank profiles, a notion of precedence has to be used: a usual way to compute the row rank profile is to search in a given row for a pivot and only move to the next row if the current row was found to be all zeros. This guarantees that each pivot will be on the first linearly independent row, and therefore the row support of will be the row rank profile. The precedence here is that the pivot’s coordinates must minimize the order for the first coordinate (the row index). As a generalization, we consider the most common preorders of the cartesian product inherited from the natural orders of each of its components and describe the corresponding search strategies, minimizing this preorder:

Row order:

iff : search for any invertible element in the first nonzero row.

Column order:

iff . search for any invertible element in the first nonzero column.

Lexicographic order:

iff or and : search for the leftmost nonzero element of the first nonzero row.

Reverse lexicographic order:

iff or and : search for the topmost nonzero element of the first nonzero column.

Product order:

iff and : search for any nonzero element at position being the only nonzero of the leading sub-matrix.

Example 13.

Consider the matrix , where each literal is a nonzero element. The nonzero elements minimizing each preorder are the following:

Row order Column order
Lexicographic order Reverse lexic. order
Product order

3.1.2 Pivot permutation

The pivot permutation moves a pivot from its initial position to the leading diagonal. Besides this constraint all possible choices are left for the remaining values of the permutation. Most often, it is done by row or column transpositions, as it clearly involves a small amount of data movement. However, these transpositions can break the precedence relations in the set of rows or columns, and can therefore prevent the recovery of the rank profile information. A pivot permutation that leaves the precedence relations unchanged will be called -monotonically increasing.

Definition 14.

A permutation of is called -monotonically increasing if its last values form a monotonically increasing sequence.

In particular, the last rows of the associated row-permutation matrix are in row echelon form. For example, the cyclic shift between indices and , with defined as , that we will call a -rotation, is an elementary -monotonically increasing permutation.

Example 15.

The -rotation is a -monotonically increasing permutation. Its row permutation matrix is .

Monotonically increasing permutations can be composed as stated in Lemma 16.

Lemma 16.

If is a -monotonically increasing permutation and a -monotonically increasing permutation with then the permutation is a -monotonically increasing permutation.

Proof..

The last values of are the image of a sub-sequence of values from the last values of through the monotonically increasing function . ∎

Therefore an iterative algorithm, using rotations as elementary pivot permutations, maintains the property that the permutation matrices and at any step are -monotonically increasing. A similar property also applies with recursive algorithms.

3.2 How to reveal rank profiles

A PLUQ decomposition reveals the row (resp. column) rank profile if it can be read from the first values of the permutation matrix (resp. ). Equivalently, by Lemma 6, this means that the row (resp. column) support of the pivoting matrix equals that of the rank profile matrix.

Definition 17.

The decomposition reveals:

  1. the row rank profile if ,

  2. the column rank profile if ,

  3. the rank profile matrix if .

Example 18.

has for rank profile matrix over . Now the pivoting matrix obtained from a PLUQ decomposition with a pivot search operation following the row order (any column, first nonzero row) could be the matrix . As these matrices share the same row support, the matrix reveals the row rank profile of .

Remark 19.

Example 18 suggests that a pivot search strategy minimizing row and column indices could be a sufficient condition to recover both row and column rank profiles at the same time, regardless the pivot permutation. However, this is unfortunately not the case. Consider for example a search based on the lexicographic order (first nonzero column of the first nonzero row) with transposition permutations, run on the matrix: . Its rank profile matrix is whereas the pivoting matrix would be , which does not reveal the column rank profile. This is due to the fact that the column transposition performed for the first pivot changes the order in which the columns will be inspected in the search for the second pivot.

We will show that if the pivot permutations preserve the order in which the still unprocessed columns or rows appear, then the pivoting matrix will equal the rank profile matrix. This is achieved by the monotonically increasing permutations. Theorem 20 shows how the ability of a PLUQ decomposition algorithm to recover the rank profile information relates to the use of monotonically increasing permutations. More precisely, it considers an arbitrary step in a PLUQ decomposition where pivots have been found in the elimination of an leading sub-matrix of the input matrix .

Theorem 20.

Consider a partial PLUQ decomposition of an matrix :

where is lower triangular and is upper triangular, and let be some leading sub-matrix of , for . Let be a PLUQ decomposition of . Consider the PLUQ decomposition

Consider the following clauses:

  1. is -monotonically increasing or ( is -monotonically increasing and )

  2. is -monotonically increasing or ( is -monotonically increasing and )

Then,

  1. if (i) or (ii) or (iii) then

  2. if (vii) then ((i) and (iv)) ;

  3. if (viii) then ((ii) and (v)) ;

  4. if (vii) and (viii) then (iii) and (vi) .

Proof..

Let and where is and is . On one hand we have

(1)

On the other hand,

Let and denote by the leading sub-matrix of .

  1. The clause (i) or (ii) or (iii) implies that all pivots of the partial elimination were found within the sub-matrix . Hence and we can write and , and the matrix writes

    (2)

    Now as a sub-matrix of of rank and since

    where the first term, , has rank and the second term has a disjoint row support.

    Finally, consider the term of equation (2). As its row support is disjoint with that of the pivot rows of , it has to be composed of rows linearly dependent with the pivot rows of to ensure that . As its column support is disjoint with that of the pivot columns of , we conclude that it must be the zero matrix. Therefore the leading sub-matrix of is zero.

  2. From (a) we know that . Thus . Recall that . No pivot row of can be made linearly dependent by adding rows of , as the column position of the pivot is always zero in the latter matrix. For the same reason, no pivot row of can be made linearly dependent by adding rows of . From (i), the set of pivot rows of is , which shows that

    (3)

    Let be the map representing the sub-permutation (i.e. such that ). If is -monotonically increasing, the matrix has full column rank and is in column echelon form, which implies that

    (4)

    since has full row rank. If is monotonically increasing, we can write , where the matrix is in column echelon form. If , the matrix writes . Hence we have which also implies

    From equation (3.2), the row support of is that of , which is the union of the row support of these two terms as they are disjoint. Under the conditions of point (b), this row support is the union of and , which is, from (4) and (3), .

  3. Similarly as for point (b).

  4. From (a) we have still . Now since , there is no other nonzero element in than those in and . The row and column support of and that of are disjoint. Hence

    (5)

    If both and are -monotonically increasing, the matrix is in column echelon form and the matrix in row echelon form. Consequently, the matrix is a copy of the matrix with zero-rows and zero-columns interleaved, which does not impact the linear dependency relations between the nonzero rows and columns. As a consequence

    (6)

    Now if is -monotonically increasing, is -monotonically increasing and , then, using notations of point (b), where is in column echelon form. Thus for the same reason. The symmetric case where is -monotonically increasing and works similarly. Combining equations (3.2), (5) and (6) gives .

4 Algorithms for the rank profile matrix

Using Theorem 20, we deduce what rank profile information is revealed by a PLUQ algorithm by the way the search and the permutation operations are done. Table 2 summarizes these results, and points to instances known in the literature, implementing the corresponding type of elimination. More precisely, we first distinguish in this table the ability to compute the row or column rank profile or the rank profile matrix, but we also indicate whether the resulting PLUQ decomposition preserves the monotonicity of the rows or columns. Indeed some algorithm may compute the rank profile matrix, but break the precedence relation between the linearly dependent rows or columns, making it unusable as a base case for a block algorithm of higher level.

Search Row Perm. Col. Perm. Reveals Monotonicity Instance
Row order Transposition Transposition RowRP Ibarra et al. (1982); Jeannerod et al. (2013)
Col. order Transposition Transposition ColRP Keller-Gehrig (1985); Jeannerod et al. (2013)
Lexico. Transposition Transposition RowRP Storjohann (2000)
Transposition Rotation RowRP, ColRP, Col. here
Rotation Rotation RowRP, ColRP, Row, Col. here
Rev. lexico. Transposition Transposition ColRP Storjohann (2000)
Rotation Transposition RowRP, ColRP, Row here
Rotation Rotation RowRP, ColRP, Row, Col. here
Product Rotation Transposition RowRP Row here
Transposition Rotation ColRP Col here
Rotation Rotation RowRP, ColRP, Row, Col. Dumas et al. (2013)
Table 2: Pivoting Strategies revealing rank profiles

4.1 Iterative algorithms

We start with iterative algorithms, where each iteration handles one pivot at a time. Here Theorem 20 is applied with , and the partial elimination represents how one pivot is being treated. The elimination of is done by induction.

4.1.1 Row and Column order Search

The row order pivot search operation is of the form: any nonzero element in the first nonzero row. Each row is inspected in order, and a new row is considered only when the previous row is all zeros. With the notations of Theorem 20, this means that is the leading sub-matrix of , where is the index of the first nonzero row of . When permutations and , moving the pivot from position to are transpositions, the matrix is the element of the canonical basis. Its row rank profile is which is that of the leading sub-matrix . Finally, the permutation is -monotonically increasing, and Theorem 20 case (b) can be applied to prove by induction that any such algorithm will reveal the row rank profile: . The case of the column order search is similar.

4.1.2 Lexicographic order based pivot search

In this case the pivot search operation is of the form: first nonzero element in the first nonzero row. The lexicographic order being compatible with the row order, the above results hold when transpositions are used and the row rank profile is revealed. If in addition column rotations are used, which is -monotonically increasing. Now which is the rank profile matrix of the leading sub-matrix of . Theorem 20 case (d) can be applied to prove by induction that any such algorithm will reveal the rank profile matrix: . Lastly, the use of row rotations, ensures that the order of the linearly dependent rows will be preserved as well. Algorithm 2 is an instance of Gaussian elimination with a lexicographic order search and rotations for row and column permutations.

The case of the reverse lexicographic order search is similar. As an example, the algorithm in (Storjohann, 2000, Algorithm 2.14) is based on a reverse lexicographic order search but with transpositions for the row permutations. Hence it only reveals the column rank profile.

4.1.3 Product order based pivot search

The search here consists in finding any nonzero element such that the leading sub-matrix of is all zeros except this coefficient. If the row and column permutations are the rotations and , we have . Theorem 20 case (d) can be applied to prove by induction that any such algorithm will reveal the rank profile matrix: . An instance of such an algorithm is given in (Dumas et al., 2013, Algorithm 2). If (resp. ) is a transposition, then Theorem 20 case (c) (resp. case (b)) applies to show by induction that the columns (resp. row) rank profile is revealed.

4.2 Recursive algorithms

A recursive Gaussian elimination algorithm can either split one of the row or column dimension, cutting the matrix in wide or tall rectangular slabs, or split both dimensions, leading to a decomposition into tiles.

4.2.1 Slab recursive algorithms

Most algorithms computing rank profiles are slab recursive Ibarra et al. (1982); Keller-Gehrig (1985); Storjohann (2000); Jeannerod et al. (2013). When the row dimension is split, this means that the search space for pivots is the whole set of columns, and Theorem 20 applies with . This corresponds to either a row or a lexicographic order. From case( b), one shows that, with transpositions, the algorithm recovers the row rank profile, provided that the base case does. If in addition, the elementary column permutations are rotations, then case (d) applies and the rank profile matrix is recovered. Finally, if rows are also permuted by monotonically increasing permutations, then the PLUQ decomposition also respects the monotonicity of the linearly dependent rows and columns. The same reasoning holds when splitting the column dimension.

4.2.2 Tile recursive algorithms

Tile recursive Gaussian elimination algorithms Dumas et al. (2013); Malaschonok (2010); Dumas and Roch (2002) are more involved, especially when dealing with rank deficiencies. Algorithm 1 recalls the tile recursive algorithm that the authors proposed in Dumas et al. (2013).

a matrix over a field
: and permutation matrices, , the rank of
where and are resp. unit lower and upper triangular, and
if  then
     Call a base case implementation
end if
Split where is .
Decompose
Here .