Symmetric indefinite triangular factorization revealing the rank profile matrix
Jean-Guillaume Dumas Clement Pernet

=1





















plain







 





  
  
  
  
  
  
  














theoremTheorem
remarkRemark
corollaryCorollary
definitionDefinition
lemmaLemma
propositionProposition
problemProblem































⋮
  vdots


⋱
  
    ddots
    





smatrix[[ ]








matrix[1][*@̧MaxMatrixCols c]
  -ifnextchar@ifnextchar
  #1




















Symmetric indefinite triangular factorization revealing the rank profile matrix[This work is partly funded by .]
    Jean-Guillaume Dumas[
  Université Grenoble Alpes.
  Laboratoire Jean Kuntzmann, CNRS, UMR 5224.
  700 avenue centrale, IMAG - CS 40700,
  38058 Grenoble, cedex 9
  France.
mailto:Jean-Guillaume.Dumas@univ-grenoble-alpes.fr,Clement.Pernet@univ-grenoble-alpes.fr{firstname.lastname}@univ-grenoble-alpes.fr]
Clément Pernet[2]
    December 30, 2023
=======================================================================================================================================================================================================================================================================================================================================




We present a novel recursive algorithm for reducing a symmetric matrix to a
triangular factorization which reveals the rank profile matrix.  
That is, the algorithm computes a factorization  P^TAP = LDL^T where P is a
permutation matrix, L is lower triangular with a unit diagonal and  D is
symmetric block diagonal with 1×1 and 2×2
antidiagonal blocks.
The novel algorithm requires O(n^2r^ω-2) arithmetic operations.
Furthermore, experimental results demonstrate that our algorithm can even
be slightly more than twice as fast as the state of the art unsymmetric Gaussian
elimination  in most cases, that is it achieves
approximately the same computational speed. 
By adapting the pivoting strategy developed in the unsymmetric case, we show
how to recover the rank profile matrix from the permutation matrix and the
support of the block-diagonal matrix.
There is an obstruction in characteristic 2 for revealing the rank profile
matrix which requires to relax the shape of the block diagonal by allowing the
2-dimensional blocks to have a non-zero bottom-right coefficient.
This relaxed decomposition can then be transformed into a standard 
PLDL^TP^T
decomposition at a negligible cost.






§ INTRODUCTION


Computing a triangular factorization of a symmetric matrix is a commonly used
kernel to solve symmetric linear systems, or to compute the signature of
symmetric bilinear forms.
Besides the fact that it is expected to save half of the arithmetic cost of a standard
(non-symmetric) Gaussian elimination, it can also recover invariants, such as
the signature, specific to symmetric matrices, and thus, e.g., be used to
certify positive or negative definite or
semidefiniteness <cit.>.

It is a fundamental computation in numerical linear algebra, and is therefore most
often presented in the setting of real matrices. When the matrix is positive
definite, the Cholesky factorization can be defined: A=LL^T, where L is
lower triangular for which square roots of diagonal elements have to be
extracted. Alternatively, gathering the diagonal elements in a central diagonal
matrix yields the LDLT factorization A=LDL^T which no longer requires square
roots. Similarly as for the LU decomposition, it is only defined for matrices
with generic rank profile, i.e. having their r=rank(A) first
leading principal minors non-zero. For arbitrary matrices, symmetric
permutations may lead to the former situations:
PAP^T=LDL^T. However, this is
unfortunately not always the case. For instance there is no permutation
P such that 0   1
1   0 has a LDLT
factorization with a diagonal D.
This lead to a series of generalizations where the matrix D was replaced first
by a tridiagonal symmetric matrix by Parlett and Reid <cit.>,
improved by Aasen <cit.>, achieving half the arithmetic cost
of Gaussian elimination. Bunch and Kaufman then replaced this tridiagonal matrix
by a block diagonal composed of 1 or 2-dimensional diagonal blocks.



  
Pivoting
In numerical linear algebra, the choice of the permutation matrix is mainly
driven by the need to ensure a good numerical quality of the decomposition.
Bunch and Parlett <cit.> use a full pivoting technique,
requiring a cubic number of tests. Bunch and Kaufman pivoting strategy,
implemented in LAPACK, uses a partial pivoting requiring only a quadratic number
of tests.

In the context of exact linear algebra, for instance when computing over a finite
field, numerical stability is no longer an issue. However, the computation of
echelon forms and rank profiles, central in many applications, impose further
constraints on the pivoting. A characterization of the requirements for the
pivoting strategy is given in <cit.> so that a PLUQ
decomposition can reveal these rank profiles and echelon forms in the
non-symmetric case.
In particular, it is shown that pivot selection minimizing the lexicographic order
on the coordinate of the pivot, combined with row and column rotations to move
the pivot to the diagonal, enable the computation of the rank profile matrix,
an invariant from which all rank profile information, the row and the column
echelon form can be recovered.



  
Recursive algorithms
As in numerical linear algebra, we try to gather arithmetic operations in level
3 BLAS operations (matrix multiplication based), for it delivers the best
computation throughput.
Numerical software often use tiled implementations, especially when the pivoting is more
constrained by the symmetry <cit.>,
or in order to define communication avoiding variants <cit.>.
In exact linear algebra sub-cubic matrix multiplication, such as Strassen's
algorithm, can be extensively used with no numerical instability issues. This
led to the design of recursive algorithms, which was proven successful
in the unsymmetric case, including for shared memory parallel computations <cit.>. 

































  
Contribution
The contribution here is to propose a recursive algorithm producing a symmetric
factorization PLDLTPT  over any field,
from which the rank profile matrix of the input can be recovered. This algorithm
is a recursive variant of Bunch and Kaufman's algorithm <cit.> where the pivoting
strategy has been replaced by the one developped previously by the authors in
the unsymmetric case <cit.>. Compared to the 
recursive adaptation of Aasen's algorihtm in <cit.>, our
algorithm leads to a similar data partitionning but does not suffer from an
arithmetic overhead compared to Aasen's algorithm.
Our algorithms has time complexity O(n^2r^ω-2) where ω is an
admissible exponent for matrix multiplication and r is the rank of the input
matrix.
With ω=3, the leading constant in the time complexity is 1/3, matching
that of the best alternative algorithms based on cubic time linear algebra.






In Section <ref> we show that in characteristic two the rank profile
matrix can not always be revealed by a symmetric factorization with antidiagonal
blocks: sometimes antitriangular blocks are also required.  
Then we recall in Section <ref> the main required level 3 linear algebra
subroutines. In Section <ref> we present the main recursive
algorithm. An alternative iterative Crout variant is 
presented in Section <ref> to be used as a base case in the recursion.
We finally show, in Section <ref>, experiments of the resulting
implementation over a finite field. They demonstrate the efficiency of cascading
the recursive algorithm with the base case variant, especially with matrices involving a
lot of pivoting. They finally confirm a speed-up by a factor of about 2 compared to the state of the art
unsymmetric Gaussian elimination.



§ THE SYMMETRIC RANK PROFILE MATRIX



 §.§ The pivoting matrix

Theorem <ref> recalls the definition of the rank profile matrix.

   Let A∈^m× n. There exists a unique m× n
 {0,1}-matrix A with  r 1's in rook placement
 of which every leading sub-matrix has the same
 rank as the corresponding
 
 leading sub-matrix of A.
 This matrix is called the rank profile matrix of A.



A symmetric matrix has a symmetric rank profile matrix.


  Otherwise, the rank of some leading submatrix of A and the same
  leading submatrix of A^T would be different which is absurd.


Also, any symmetric matrix has a triangular decomposition
A=PLDL^T P^T where L is unit lower triangular, D is block diagonal,
formed by 1-dimensional scalar blocks or 2-dimensional blocks of the form
0   x
x   0 and P a permutation matrix.

We here further define Ψ as the support matrix of D: namely, a
block 
diagonal {0,1}-matrix such that D=ΨD, with D a
diagonal matrix. 
  
    The pivoting matrix of a PLDLTPT decomposition is the matrix Π=PΨP^T.
  

  
    A PLDLTPT reveals the rank profile matrix of a symmetric matrix A if its
    pivoting matrix equals the rank profile matrix of A.
  



 §.§ Antitriangular blocks in characteristic two

In zero or odd characteristic, we show next that one can always find such a
PLDLTPT decomposition revealing the rank profile matrix. In  characteristic two,
however, this is not always possible.
  
    In characteristic 2, there is no symmetric indefinite elimination revealing the rank
  profile matrix of A=[ 0 1; 1 1 ].
  

Let J be the 2×2 anti-diagonal identity matrix. This is also the rank
profile matrix of A.
Now, we let L=[ 1 0; x 1 ],
D=[ y 0; 0 z ]. As the permutation matrices
involved, P and Ψ, can only be either the identity matrix or J, there
are then four cases: 


  * A=L·D·L^T=[       y     x y;     x y x^2 y+z ],
  but y=0 and (D)=(A)=2 are incompatible.

  * A=J·L·D·L^T·J^T=[ x^2y+z    x y;    x y      y ], but J·I·J≠J=A.

  * A=L·D·J·L^T=[     0     y;     z xy+xz ],
  but we need y=z for the symmetry and then 2xy=0≠ 1 in
  characteristic 2. 

  * A=J·L·D·J·L^T·J^T=[ xy+xz     z;     y     0 ] but the bottom right coefficient of A is non zero. 



However, one can generalize the PLDLTPT decomposition to a block diagonal
matrix D having 2-dimensional blocks of the form
0   c
c   d (lower antitriangular).
Then the support matrix Ψ of
D is the block diagonal {0,1} matrix such that D=ΨD, with
D an upper triangular bidiagonal matrix (or equivalently such that
D=DΨ, with D lower triangular bidiagonal). 

With these generalized definitions, we show in Section <ref>, that
there exists RPM-revealing PLDLTPT decompositions.




 §.§ Antitriangular decomposition

Then, such a generalized decomposition can always be further reduced to a strict
PLDLTPT decomposition by eliminating each of the antitriangular blocks.  
For this, the observation is that in
characteristic two, a symmetric lower antitriangular 2×2 block is
invariant under any symmetric triangular transformation:
[ 1; x 1 ][ c; c d ][ 1 x; 1 ]=
[     c;     c 2cx+d ]≡[ c; c d ] 2=
[ 1; 1 ][ c; c d ][ 1; 1 ].
Thus for each 2×2 block in a tridiagonal decomposition, the
corresponding 2×2 diagonal block in L can be replaced by I_2, via a
multiplication by 1   
-x   1.

Further, we have that: 
c
c   d=
J1   
c/d   1d   
   -c^2/d1   c/d
   1J.

Now J commutes with the identity I_2 matrix. 
Therefore we have that:
1   
x   11   
-x   1J1   
c/d   1
=
J1   
x   11   
-x   11   
c/d   1.

Thus, to eliminate the antitriangular blocks, 
create a triangular matrix L_J that starts as the identity
and where its i,i+1 blocks corresponding to a
1   
x   1 block in L  is a
1   
c/d-x   1 block (associated to an
antitriangular c
c   d block, with d≠0, in
D).

Then replace the triangular
matrix L by L̃=L·L_J.

Also, modify the diagonal matrix D, to D̃ such that 
the c
c   d blocks of D are replaced by
d   
   -c^2/d blocks in D̃.

Finally, create a permutation matrix P_J, starting from the identity matrix,
where each identity block at position i,i+1 corresponding to an antitriangular
block in D is replaced by J. Then P̃=P·P_J.


From this we have now a symmetric PLDLTPT factorization,
A=P̃L̃D̃L̃^TP̃^T, with purely 1×1
and 2×2 antidiagonal blocks in D̃
(but then a direct access to the rank profile matrix, PΨP^T, might not be
possible from P̃ and D̃). 

In the following we present some building blocks and then algorithms
computing RPM-revealing symmetric indefinite triangular factorization.



  

§ BUILDING BLOCKS


  We recall here some of the standard algorithms from the BLAS3 <cit.> and
  LAPACK <cit.> interfaces and generalization thereof <cit.>, which will be used to
  define the main block  recursive symmetric eliminating algorithm.

  
    
  (C,A,B): general matrix multiplication. Computes C←C - AB. 
    
  (U,B):   multiply a triangular and a rectangular matrix in-place. Computes B←UB where B is m× n and U
      is upper or lower triangular.
    
  (C,U,B):  multiply a triangular and a rectangular
      matrix. Computes C←C- UB where B and C are m× n and U
      is upper or lower triangular. This is an adaptation of the BLAS3 to
      leave the B operand unchanged.
    
  (U,B):  solve a triangular system with matrix right
      hand-side. Computes B←U^-1B where B is m× n and U
      is upper or lower triangular.
    
  (C,A,D): symmetric rank k update with diagonal
      scaling. Computes the upper or lower triangular part of the symmetric
      matrix C←C - ADA^T where A is n× k and D is
      diagonal or block diagonal.
    
  (C,A,D,B): symmetric rank 2k update with diagonal
      scaling. Computes the upper or lower triangular part of the symmetric
      matrix C←C - ADB^T-BDA^T where A and B are n× k, B and D is
      diagonal or block diagonal. 
  

  In addition, we need to introduce the routine solving  Problem <ref>.
  
    
    Let  be a field of characteristic different than 2. Given a symmetric
    matrix C∈^n× n and a unit upper triangular matrix
    U∈^n× n, find an upper triangular matrix X∈^n× n
    such that X^TU+U^TX=C.
  
In characteristic 2, the diagonal of X^TU+U^TX is always zero for any matrix X and
U, hence Problem <ref> has no solution as soon as C has a non-zero
diagonal element. 

However in characteristic zero or odd, Algorithm <ref> presents a
recursive implementation of this 
routine, and is in the same time a constructive proof of the existence of such a
solution. Note that it performs a division by 2 in line <ref>, and
therefore requires that the base field has not characteristic two.


 Note that algorithm <ref> computes the solution X in
  place on the symmetric storage of C: by induction X_1 and X_3 overwrite
  C_1 and C_3, and X_2 overwrites C_2 according to the specifications of
  the generalized  routine.



  Algorithm is correct and runs in O(n^ω) arithmetic operations.
  















  Using the notations of Algorithm <ref>, let X=
  [ X_1 X_2; X_3 ]. Then exanding X^TU+U^TX gives
  
    [                   X_1^TU_1+U_1^TX_1                 X_1^TU_2 + U_1^TX_2;            (X_1^TU_2 + U_1^T X_2)^T X_2^TU_2+U_2^TX_2+X_3^TU_3+U_3^TX_3 ]
    
    =    [                        C_1                        C_2;                      C_2^T C_3-D_3 +X_3^TU_3+U_3^TX_3 ]     =C.

  which proves the correctin by induction.
  The arithmetic cost satisfy a recurrence of the form T(n)=2T(n/2)+Cn^ω
  and is therefore T(n)=O(n^ω).
  
  
  
  
  
  
  
  
  
  




§ A BLOCK RECURSIVE ALGORITHM




 §.§ Sketch of the recursive algorithm


The design  of a block recursive algorithm is based on the generalization of
the 2× 2 case into a block 2× 2 block algorithm. While scalars
could be either 0 or invertible, the difficulty in elimination algorithms, is
that a submatrix could be rank defficient but non-zero. We start here an
overview of the recursive algorithm by considering that the leading principal
block is either all zero or invertible. We will later give the general
presentation of the algorithm where its rank could be arbitrary.






Let  M∈^(m+n)×(m+n) be the symmetric matrix to be
factorized. Consider its block decomposition 
M=[   A   B; B^T   C ] where A∈^m× m and
C∈^n× n are also symmetric. 

If A is full rank, then a recursive call will produce 
A=PLDL^TP^T, and M can thus be decomposed as:

    M=
    [ P 0; 0 I ][ L 0; G I ][ D 0; 0 Z ][ L^T G^T;   0   I ][ P^T   0;   0   I ],

where G is such that PLDG^T=B and Z=C-GDG^T.
Thus G can be computed as the transpose of
D^-1L^-1P^-1B which can be
obtained by a call to , some permutations and a diagonal scaling. Then
Z is computed by a call to . A second recursive call will then
decompose Z and lead to the final factorization of M.

Now if A is the zero matrix, one is reduced to factorize the matrix 
N=0   B
B^T    C. In order to recover the rank
profile matrix, one has to first look for pivots in B before considering the
block C. Therefore diagonal pivoting is not an option here.
Then the matrix B, which we assume has full rank for the moment, can be decomposed in a
PLDUQ factorization (P and Q permutation matrices, L and U 
respectively unit lower and unit upper triangular, D is diagonal).
We then need to distinguish two cases depending on whether the field
characteristic is two or not.



  §.§.§ Zero or odd characteristic case

If the characteristic zero or odd, N can thus be decomposed as:

    N=
    [   P   0;   0 Q^T ][   L   0;   G U^T ][ 0 D; D 0 ][ L^T G^T;   0   U ][ P^T   0;   0   Q ],

where G is such that Q^T(GDU+U^TDG^T)Q=C.
To compute G, one can first permute C to get C'=QCQ^T (which remains
symmetric) and then use a call to .



  §.§.§ Characteristic two case

In characteristic two, the equation GDU+U^TDG^T=QCQ^T in unknown G has in
general  no solution (as soon as C has a non-zero diagonal element).






However, one can still relax Problem <ref> and allow the elimination to
leave a diagonal of elements not zeroed out.
Following Lemma <ref>, the idea is then to
decompose N into a block tridiagonal form: 

    N=
    [   P   0;   0 Q^T ][   L   0;   G U^T ][ 0 D; D Δ ][ L^T G^T;   0   U ][ P^T   0;   0   Q ],

where Δ is a diagonal matrix and now G is such that
Q^T(GDU+U^TDG^T+U^TΔU)Q=C.
Therefore Δ can be chosen such that the diagonal of 
C”=QCQ^T-U^TΔU=C'-U^TΔU is zero. 


As U is unit upper triangular, a simple pass over its
coefficients is sufficient to find such a Δ: let
Δ_ii=C'_ii-∑_j=1^i-1Δ_jjU_j,i^2.
The algorithm is thus to permute C to get C';
then compute Δ with the recursive relation above
and update C”=C'-U^TΔU with a . 
C” remains symmetric but with a zero diagonal and now can be applied.




 §.§ The actual recursive algorithm





  §.§.§ First phase: recursive elimination

In the general case, the leading matrices are not full rank, and we have to
consider intermediate steps.
For the symmetric matrix M∈^(m+n)×(m+n) of
Section <ref>, its leading principal block A∈^m×m is
of rank r≤m.  
Thus its actual recursive decomposition is of the form:

    A=P_1
    [ L_1; M_1 ][ D_1 ][ L_1^T M_1^T ]P_1^T,

where L_1∈^r×r is full rank
unit lower triangular, D_1∈^r×r is block diagonal with
1 or 2-dimensional diagonal blocks, and M_1∈^(m-r)×r.

Therefore, forgetting briefly the permutations, the decomposition of M
becomes: 
    M=
    [ L_1   0   0; M_1   I   0;   G   0   I ][ D_1   0   0;   0   0   Y;   0 Y^T   Z ][ L_1^T M_1^T   G^T;     0     I     0;     0     0     I ],

where Y is such that 
B=[ L_1; M_1 ]D_1G^T+[ 0; Y ].

From this point on, there remains to factorize the submatrix
0   Y
Y^T   Z. This will be carried out by the
algorithm described in the next section, working on a matrix with a zero 
leading principal submatrix.  



Supposing for now that this is possible, Algorithm <ref> summarizes
the whole procedure.







  §.§.§ Second phase: off-diagonal pivoting


Consider N=0   B
B^T    C, where
B∈^m×n, with m≤n, has now an arbitrary rank
r≤m. Then its
PLDUQ decomposition is of the form

    B=P[ L_1; M_1 ][ D_1 ][ U_1 V_1 ]Q,

with D_1 diagonal, and L_1 and U_1 unit square triangular matrices, all
three of order r.
Then consider a conformal block decomposition of
QCQ^T=C_1   C_2
C_2^T   C_3  where C_1 is r× r.
It remains to eliminate C_1 and C_2 with the pivots found in B, which
leads to the following factorization:

    N=
    [   P_1; Q_1^T ][   L_1     0     0;   M_1     0     0;   G_1 U_1^T     0;   G_2 V_1^T     I ][   0 D_1   0; D_1   0   0;   0   0   Z ]×
    [ L_1^T M_1^T G_1^T G_2^T;     0     0   U_1   V_1;     0     0     0     I ][ P_1^T;   Q_1 ]


where G_1 satisfies 
    U_1^TD_1G_1^T+G_1D_1U_1=C_1

  and
G_2=(C_2-V_1^TD_1G_1^T)U_1^-1D_1^-1 and Z=C_3-(V_1^TD_1G_2^T+G_2D_1V_1).



In order to produce a LDLT decomposition, there still remains to perform
permutations to


  * compact the leading elements of the lower triangular matrix into a
  2r× 2r invertible leading triangular submatrix,

  * make the D_1
D_1 matrix block
  diagonal with 1 or 2-dimensional diagonal blocks.


The permutation matrix 
    P_c=I_r    0    0   0
    0   0   I_m-r   0
    0   I_r   0   0
    0   0   0   I_n-r,

corresponding to a block circular rotation, takes care of
condition <ref>, while preserving precedence in the
non-pivot rows. This is a requirement for the factorization to reveal the rank
profile matrix <cit.>.
The decomposition becomes

    N=
    [   P_1; Q_1^T ]P_c
    [   L_1     0     0;   G_1 U_1^T     0;   M_1     0     0;   G_2 V_1^T     I ][   0 D_1   0; D_1   0   0;   0   0   Z ]×
    [ L_1^T G_1^T M_1^T G_2^T;     0   U_1     0   V_1;     0     0     0     I ]P_c^T  [ P_1^T;   Q_1 ]

























In order to achieve Condition <ref>, we will transform the
matrix 0   D_1
D_1   0 into the block diagonal
  matrix Diag( 0   d_i
d_i   0) where d_i is
    the ith diagonal element in D_1.
To describe the process, we will focus on the matrix
 
    N_2=
    [   L_1     0;   G_1 U_1^T ][   0 D_1; D_1   0 ][ L_1^T G_1^T;     0   U_1 ]
    =L·Δ·L^T,

and consider a splitting in halves of the matrix
D_1=D_11
   D_12 where D_11 has order
r_1 and D_12 order r_2. This leads to the conformal decompostion

    [   L_11      0      0      0;   L_12   L_13      0      0;   G_11   G_14 U_11^T      0;   G_12   G_13 U_12^T U_13^T ][    0    0 D_11    0;    0    0    0 D_12; D_11    0    0    0;    0 D_12    0    0 ][ L_11^T L_12^T G_11^T G_12^T;      0 L_13^T G_14^T G_13^T;      0      0   U_11   U_12;      0      0      0   U_13 ]


Then considering the permutation matrix

    P_d=I_r_1    0    0   0
    0   0   I_r_2   0
    0   I_r_1   0   0
    0   0   0   I_r_2,

one  can form 
P_d^TΔP_d=0   D_11   0   0
D_11   0   0   0
0   0   0   D_12
0   0   D_12   0
and
P_d^TLP_d=
L_11   0   0   0
G_11   U_11^T   G_14   0
L_12   0   L_13   0
G_12   U_12^T   G_13   U_13^T.

  Applying this process recursively changes Δ into the desired
  block diagonal form. Then the transformation of L will remain
  lower triangular if and only if all G_14 matrices are zero: this means
  that  G_1 must be lower triangular in the first place.

Finding G_1 lower triangular satifying Equation (<ref>), is an instance
of Problem <ref> for which the routine provides a solution.


Note that the actual permutation to transform 
0   D_1
D_1   0 into a 2×2-blocks diagonal
matrix is a permutation matrix, P_i, resulting from the one by one
interleaving of the rows of I_r   0 and
0   I_r. If
e_i=0…0   1   0…0^T is the i-th canonical vector, then:

    P_i=[   e_1 e_r+1   e_2 e_r+2     …   e_r  e_2r ].

Similarly the triangular factor of the factorization is thus a one by one
interleaving of the rows of L_1   0 and
G_1   U_1^T as well as a one by one interleaving of
the columns L_1
G_1 and
0
U_1^T, which overall remains triangular.

Finally, a call to Algorithm <ref> produces a factorization for the
remaining Z block and a final block rotation,  

    [  I_2r     0     0;     0     0 I_m-r;     0 I_n-r     0 ],
 
moves the intermediate zero rows and columns to the bottom right.
The full algorithm is presented in details in Algorithm <ref> (for zero
or odd characteristic, the characteristic two case being presented afterwards in
Section <ref>).






 §.§ Characteristic two

The case of the characteristic two can be handled similarly, just computing the extra
diagonal and updating after the decomposition, as sketched in
Section <ref>. 
Indeed, the only issue is the division by 2 in , which is removed if the
diagonal of C'_1 is zero. 
Therefore, Algorithm <ref> is unchanged, the block diagonal matrix
just has lower symmetric antitriangular 2×2 blocks instead of only
antidiagonal ones. The only few additional operations appear in 
Algorithm <ref>  and are the contents of the
"characteristic()=2…" branchings.














































Then the tridiagonal form with symmetric antitriangular
2×2 blocks thus obtained by Algorithm <ref> can be used to
either reveal the rank profile matrix
(via computing Ψ, the support matrix of D, and the pivoting matrix
=PΨ P^T) or a PLDLTPT factorization, both at an
extra linear cost, as shown in Section <ref>.

Overall, we have proven:
 Algorithm <ref> correctly computes a symmetric
  indefinite PLDLTPT factorization revealing the rank profile matrix. 




§ BASE CASE ITERATIVE VARIANT


The recursion of Algorithm <ref> should not be performed all the way
to a dimension 1 in practice.
For implementations over a finite field, it would induce an unnecessary large
number of modular reductions and a significant amount of data movement for the
permutations. Instead, we propose in Algorithm <ref> an iterative algorithm computing a PLDLTPT
revealing the rank profile matrix to be used as a base case in the recursion.

This iterative algorithm has the following features:

  
  * it uses a pivot search minimizing the lexicographic order
(following the caracterization in <cit.>): if the diagonal element
    of the current row is 0, the pivot is chosen as the first non-zero element
    of the row, unless the row is all zero, in which case, it is searched in the
    following row;
  
  * the pivot is permuted with cyclic shifts on the row and columns, so as
    to leave the precedence in the remaining rows and columns unchanged.
  
  * the update of the unprocessed part in the matrix is delayed following
    the scheme of a Crout elimination 
    schedule <cit.>. It does not only improves efficiency thanks to a
    better data locality, but it also reduces the amount of modular reductions,
    over a finite field, as shown for the unsymmetric case
    in <cit.>. 


We denote by ρ_i,n the cyclic shift permutation of order n moving
element i to the first position: ρ_i,n=(i,0,1,…, i-1,i+1,… n-1).
Indices are 0 based, index ranges are excluding their upper bound. For instance,
A_i,0..r denotes the r first elements of the i+1st row of A, and
A_0..r,0..r is the 0-dimensional matrix when r=0.







§ EXPERIMENTS

We now report on experiments of an implementation of these algorithms in the
FFLAS-FFPACK library <cit.>, dedicated to dense linear algebra over finite
fields. We used the version committed under the reference
https://github.com/linbox-team/fflas-ffpack/commit/e12a9989f07a9e128d9b4dd59681e8f62e8fe3b1e12a998
of the master branch. It was compiled with gcc-5.4 and was linked with the
numerical library OpenBLAS-0.2.18. Experiments are run on a single core of an an
Intel Haswell i5-4690, 643.5GHz.

Computation speed are normalized as effective Gfops, an estimate of the
number of field operations that an algorithm with classic matrix arithmetic
would perform per second, divided by the computation time. For a matrix of order n
and rank r, we defined this as:

    Effective Gfops = (r^3/3+n^2r-r^2n) / (10^9 ×time).

All experiments are over the 23-bits finite field /8388593.

Figure <ref> compares the computation speed of the pure recursive
algorithm, the base case algorithm and a cascade of these two, with a threshold
set to its optimum value from experiments on this machine.
Remark that the pure recursive variant performs rather well with generic rank
profile matrices, while matrices with uniformly random rank profile matrix make
this variant very slow, due to an excessive amount of pivoting.
As expected, the base case Crout variant speeds up these instances for small dimensions, but
then its performance stagnate on large dimensions, due to poor cache efficiency. Lastly the
cascade algorithm combines the benefits of the two variants and therefore
performs best in all settings.
We here used a threshold n=128 for the experiments with random RPM matrices, but of
only n=48 for generic rank profile matrices, since the recursive variant becomes
competive much earlier. In most cases, the rank profile structure of given
matrices is unknown a priori, making the setting of this threshold
speculative. One could instead implement an introspective strategy, updating
the threshold from experimenting with running instances.














































Table <ref> compares the computation time of the symmetric
decomposition algorithm with that of the unsymmetric case (running the PLUQ
algorithm of <cit.>). These experiments confirm a speed-up factor
of about 2 between these routines, which is the expected gain in the constant in
the time complexity. Note that on large instances, the PLUQ elimination performs
better with random RPM instances than generic rank profiles, contrarily to the
LDLT routine. This is due to the lesser amount of arithmetic operations when the
RPM is random (some intermediate submatrices being rank deficient). On the other
hand, these matrices generate more off-diagonal pivots, which cause more
pivoting in LDLT than in PLUQ, explaining the slow down for the symmetric case.
















plainurl ]