Computing the rank of big sparse matrices modulo p using gaussian - - PowerPoint PPT Presentation

▶
computing the rank of big sparse matrices modulo p using
SMART_READER_LITE
LIVE PREVIEW

Computing the rank of big sparse matrices modulo p using gaussian - - PowerPoint PPT Presentation

Jan 06, 2023 407 likes •961 views

Computing the rank of big sparse matrices modulo p using gaussian elimination Charles Bouillaguet 1 Claire Delaplace 2 12 CRIStAL, Universit de Lille 2 IRISA, Universit de Rennes 1 JNCF, 16 janvier 2017 Claire Delaplace Sparse elimination Mod

slide-1
SLIDE 1

Computing the rank of big sparse matrices modulo p using gaussian elimination

Charles Bouillaguet 1 Claire Delaplace 2

12CRIStAL, Université de Lille 2IRISA, Université de Rennes 1

JNCF, 16 janvier 2017

Claire Delaplace Sparse elimination Mod p JNCF 2017 1 / 20

slide-2
SLIDE 2

Sparse Linear Algebra

Background

Sparse Linear Algebra Modulo p (coefficients : int)

Operations

Rank Linear systems Kernel etc...

Claire Delaplace Sparse elimination Mod p JNCF 2017 2 / 20

slide-3
SLIDE 3

Sparse Linear Algebra

Background

Sparse Linear Algebra Modulo p (coefficients : int)

Operations

Rank Linear systems Kernel etc...

Two families of Algorithms

Direct methods (Gaussian Elimination, LU, ...) Iterative methods (Wiedemann, Lanczos...)

Claire Delaplace Sparse elimination Mod p JNCF 2017 2 / 20

slide-4
SLIDE 4

Sparse Linear Algebra

Related Work

Algorithms

Comparison between a sparse gaussian elimination and the Wiedmann algorithm: [Dumas & Villard 02] Direct methods in the numerical world (e.g. [Davis 06]) Pivots selection heuristic for Gröbner Basis Matrices [Faugère & Lacharte 10]

Software

Exact: not much (LinBox, GBLA (Gröbner basis), MAGMA) Numeric: many (SuperLU, UMFPACK, ...)

Claire Delaplace Sparse elimination Mod p JNCF 2017 3 / 20

slide-5
SLIDE 5

Sparse Linear Algebra

PLUQ Factorization

L U = P× A ×Q−1 L has non zero diagonal U has unit diagonal A can be rectangular A can be rank deficient

Claire Delaplace Sparse elimination Mod p JNCF 2017 4 / 20

slide-6
SLIDE 6

Sparse Linear Algebra

Right-looking and Left-looking LU

Usual right-looking Algorithm: L U S Left-looking GPLU Algorithm [Gilbert & Peierls 88] L U A

Claire Delaplace Sparse elimination Mod p JNCF 2017 5 / 20

slide-7
SLIDE 7

Sparse Linear Algebra

Right-looking and Left-looking LU

Usual right-looking Algorithm: Data access L U Left-looking GPLU Algorithm [Gilbert & Peierls 88] L U A

Claire Delaplace Sparse elimination Mod p JNCF 2017 5 / 20

slide-8
SLIDE 8

Sparse Linear Algebra

Right-looking and Left-looking LU

Usual right-looking Algorithm: L U S Left-looking GPLU Algorithm [Gilbert & Peierls 88] L U A

Claire Delaplace Sparse elimination Mod p JNCF 2017 5 / 20

slide-9
SLIDE 9

Sparse Linear Algebra

Right-looking and Left-looking LU

Usual right-looking Algorithm: Data access L U S Left-looking GPLU Algorithm [Gilbert & Peierls 88] U A

Claire Delaplace Sparse elimination Mod p JNCF 2017 5 / 20

slide-10
SLIDE 10

Sparse Linear Algebra

Right-looking and Left-looking LU

Usual right-looking Algorithm: L U S Left-looking GPLU Algorithm [Gilbert & Peierls 88] L U A

Claire Delaplace Sparse elimination Mod p JNCF 2017 5 / 20

slide-11
SLIDE 11

Sparse Linear Algebra

An up-looking variant of the GPLU

Row-by-row version Adapted to Computer Algebra L U A

Claire Delaplace Sparse elimination Mod p JNCF 2017 6 / 20

slide-12
SLIDE 12

Sparse Linear Algebra

An up-looking variant of the GPLU

Row-by-row version Adapted to Computer Algebra Data access L A

Claire Delaplace Sparse elimination Mod p JNCF 2017 6 / 20

slide-13
SLIDE 13

Sparse Linear Algebra

An up-looking variant of the GPLU

Row-by-row version Adapted to Computer Algebra L U A

Claire Delaplace Sparse elimination Mod p JNCF 2017 6 / 20

slide-14
SLIDE 14

Sparse Linear Algebra

GPLU Algorithm: Application to Exact Linear Algebra

Never been used before for exact computations We implemented it We benchmarked it against LinBox (sparse right-looking)

Our Benchmarks show:

GPLU work best when U is very sparse Sometimes GPLU outperform the right-looking algorithm (often) Sometimes the right-looking algorithm outperform GPLU (less often) = ⇒ Can we take advantage of both methods?

Claire Delaplace Sparse elimination Mod p JNCF 2017 7 / 20

slide-15
SLIDE 15

A new hybrid algorithm

Our Work: A New Hybrid Algorithm (CASC 2016)

Description

Find many pivots without performing any arithmetical operations. Compute the Schur complement S, using an up-looking algorithm. Compute the rank of S.

Claire Delaplace Sparse elimination Mod p JNCF 2017 8 / 20

slide-16
SLIDE 16

A new hybrid algorithm

Our Work: A New Hybrid Algorithm (CASC 2016)

Description

Find many pivots without performing any arithmetical operations. Compute the Schur complement S, using an up-looking algorithm. Compute the rank of S.

Rank of S

Recurse Dense rank computation Wiedemann Algorithm ...

Claire Delaplace Sparse elimination Mod p JNCF 2017 8 / 20

slide-17
SLIDE 17

A new hybrid algorithm

Example:

Claire Delaplace Sparse elimination Mod p JNCF 2017 9 / 20

slide-18
SLIDE 18

A new hybrid algorithm

Example:

Set of pivots found without any arith- metical operations

Claire Delaplace Sparse elimination Mod p JNCF 2017 9 / 20

slide-19
SLIDE 19

A new hybrid algorithm

Example:

Schur Complement

Claire Delaplace Sparse elimination Mod p JNCF 2017 9 / 20

slide-20
SLIDE 20

A new hybrid algorithm

Example:

Claire Delaplace Sparse elimination Mod p JNCF 2017 9 / 20

slide-21
SLIDE 21

A new hybrid algorithm

Example:

Claire Delaplace Sparse elimination Mod p JNCF 2017 9 / 20

slide-22
SLIDE 22

A new hybrid algorithm

Example:

Claire Delaplace Sparse elimination Mod p JNCF 2017 9 / 20

slide-23
SLIDE 23

A new hybrid algorithm

Example:

Parallelizable

Claire Delaplace Sparse elimination Mod p JNCF 2017 9 / 20

slide-24
SLIDE 24

A new hybrid algorithm

Initial Pivots Selection Heuristic [Faugère & Lachartre 10]

1 2 3 4 5 6 7            

           4 2 1 3 5 6 7            

          

Description

Each row is mapped to the column of its leftmost coefficient. When several rows have the same leftmost coefficient, select the sparsest. Move the selected rows before the others and sort them by increasing position of the leftmost coefficient.

Claire Delaplace Sparse elimination Mod p JNCF 2017 10 / 20

slide-25
SLIDE 25

A new hybrid algorithm

Schur Complement Computation

P denotes the permutation that pushes the pre-computed pivots in the top of A. Ignoring permutation over the columns of A: PA =

  • A00

A01 A10 A11

  • =

L00 L10 L11

  • ·

U00 U01 U11

  • Claire Delaplace

Sparse elimination Mod p JNCF 2017 11 / 20

slide-26
SLIDE 26

A new hybrid algorithm

Schur Complement Computation

P denotes the permutation that pushes the pre-computed pivots in the top of A. Ignoring permutation over the columns of A: PA =

  • A00

A01 A10 A11

  • =

L00 L10 L11

  • ·

U00 U01 U11

  • Definition

The Schur Complement S of PA with respect to U00 is given by : S = A11 − A10U−1

00 U01

Claire Delaplace Sparse elimination Mod p JNCF 2017 11 / 20

slide-27
SLIDE 27

A new hybrid algorithm

Schur Complement Computation

P denotes the permutation that pushes the pre-computed pivots in the top of A. Ignoring permutation over the columns of A: PA =

  • A00

A01 A10 A11

  • =

L00 L10 L11

  • ·

U00 U01 U11

  • Definition

The Schur Complement S of PA with respect to U00 is given by : S = A11 − A10U−1

00 U01

Denote by (ai0 ai1) the i-th row of (A10 A11), and consider the following system : (x0 x1) · U00 U01 Id

  • = (ai0 ai1)

Claire Delaplace Sparse elimination Mod p JNCF 2017 11 / 20

slide-28
SLIDE 28

A new hybrid algorithm

Schur Complement Computation

P denotes the permutation that pushes the pre-computed pivots in the top of A. Ignoring permutation over the columns of A: PA =

  • A00

A01 A10 A11

  • =

L00 L10 L11

  • ·

U00 U01 U11

  • Definition

The Schur Complement S of PA with respect to U00 is given by : S = A11 − A10U−1

00 U01

Denote by (ai0 ai1) the i-th row of (A10 A11), and consider the following system : (x0 x1) · U00 U01 Id

  • = (ai0 ai1)

We obtain x1 = ai1 − ai0U−1

00 U01. x1 is the i-th row of S

Claire Delaplace Sparse elimination Mod p JNCF 2017 11 / 20

slide-29
SLIDE 29

A new hybrid algorithm

Schur Complement Computation

In a Nutshell:

Each row of S can be computed independently Guess if S is sparse or dense by computing some random rows If S sparse : compute S using the previous method The Schur complement computation is parallelizable

Claire Delaplace Sparse elimination Mod p JNCF 2017 12 / 20

slide-30
SLIDE 30

A new hybrid algorithm

Performance of the Hybrid Algorithm (CASC 2016)

Experiments carried on an Intel Core i7-3770 with 8 GB of RAM Experiments carried on all matrices from SIMC with integer coefficents Only one core used

Hybrid versus Right-looking (Linbox) and GPLU (time in s)

Matrix Right-looking GPLU Hybrid GL7d/GL7d24 34 276 11.6 Margulies/cat_ears_4_4 3 184 0.1 Homology/ch7-8.b4 173 0.2 0.2 Homology/ch7-8.b5 611 45 10.7

Hybrid versus Wiedmann (time in s)

Matrix Wiedmann Hybrid M0,6-D7 20397 0.8 relat8 244 2 relat9 176694 2024

Claire Delaplace Sparse elimination Mod p JNCF 2017 13 / 20

slide-31
SLIDE 31

Recent Improvement: Better Pivots Selection Heuristic

New: Improvement of the Pivots Selection Heuristic

Why

Enable us to reduce the number of the elimination steps Keep U sparse ⇒ fast elimination steps

Claire Delaplace Sparse elimination Mod p JNCF 2017 14 / 20

slide-32
SLIDE 32

Recent Improvement: Better Pivots Selection Heuristic

New: Improvement of the Pivots Selection Heuristic

Why

Enable us to reduce the number of the elimination steps Keep U sparse ⇒ fast elimination steps

How

Find the largest set: NP-Complete Find a large set using heuristics Our idea: First find the F.-L. pivots, then search for more.

Claire Delaplace Sparse elimination Mod p JNCF 2017 14 / 20

slide-33
SLIDE 33

Recent Improvement: Better Pivots Selection Heuristic

How to Find Remaining Pivots:

  • Choice of new pivots:

The pivotal submatrix must form a DAG ⇒ If there is a cycle, the entry can’t be selected.

Claire Delaplace Sparse elimination Mod p JNCF 2017 15 / 20

slide-34
SLIDE 34

Recent Improvement: Better Pivots Selection Heuristic

How to Find Remaining Pivots:

  • One candidate

Choice of new pivots:

The pivotal submatrix must form a DAG ⇒ If there is a cycle, the entry can’t be selected.

Claire Delaplace Sparse elimination Mod p JNCF 2017 15 / 20

slide-35
SLIDE 35

Recent Improvement: Better Pivots Selection Heuristic

How to Find Remaining Pivots:

  • No entries on upper

rows

Choice of new pivots:

The pivotal submatrix must form a DAG ⇒ If there is a cycle, the entry can’t be selected.

Claire Delaplace Sparse elimination Mod p JNCF 2017 15 / 20

slide-36
SLIDE 36

Recent Improvement: Better Pivots Selection Heuristic

How to Find Remaining Pivots:

  • New pivot found !

Choice of new pivots:

The pivotal submatrix must form a DAG ⇒ If there is a cycle, the entry can’t be selected.

Claire Delaplace Sparse elimination Mod p JNCF 2017 15 / 20

slide-37
SLIDE 37

Recent Improvement: Better Pivots Selection Heuristic

How to Find Remaining Pivots:

  • Choice of new pivots:

The pivotal submatrix must form a DAG ⇒ If there is a cycle, the entry can’t be selected.

Claire Delaplace Sparse elimination Mod p JNCF 2017 15 / 20

slide-38
SLIDE 38

Recent Improvement: Better Pivots Selection Heuristic

How to Find Remaining Pivots:

  • Choice of new pivots:

The pivotal submatrix must form a DAG ⇒ If there is a cycle, the entry can’t be selected.

Claire Delaplace Sparse elimination Mod p JNCF 2017 15 / 20

slide-39
SLIDE 39

Recent Improvement: Better Pivots Selection Heuristic

How to Find Remaining Pivots:

  • Three candidates

Choice of new pivots:

The pivotal submatrix must form a DAG ⇒ If there is a cycle, the entry can’t be selected.

Claire Delaplace Sparse elimination Mod p JNCF 2017 15 / 20

slide-40
SLIDE 40

Recent Improvement: Better Pivots Selection Heuristic

How to Find Remaining Pivots:

  • Choice of new pivots:

The pivotal submatrix must form a DAG ⇒ If there is a cycle, the entry can’t be selected.

Claire Delaplace Sparse elimination Mod p JNCF 2017 15 / 20

slide-41
SLIDE 41

Recent Improvement: Better Pivots Selection Heuristic

How to Find Remaining Pivots:

  • Choice of new pivots:

The pivotal submatrix must form a DAG ⇒ If there is a cycle, the entry can’t be selected.

Claire Delaplace Sparse elimination Mod p JNCF 2017 15 / 20

slide-42
SLIDE 42

Recent Improvement: Better Pivots Selection Heuristic

How to Find Remaining Pivots:

  • May be eliminated

Choice of new pivots:

The pivotal submatrix must form a DAG ⇒ If there is a cycle, the entry can’t be selected.

Claire Delaplace Sparse elimination Mod p JNCF 2017 15 / 20

slide-43
SLIDE 43

Recent Improvement: Better Pivots Selection Heuristic

How to Find Remaining Pivots:

  • Can’t be eliminated

Choice of new pivots:

The pivotal submatrix must form a DAG ⇒ If there is a cycle, the entry can’t be selected.

Claire Delaplace Sparse elimination Mod p JNCF 2017 15 / 20

slide-44
SLIDE 44

Recent Improvement: Better Pivots Selection Heuristic

How to Find Remaining Pivots:

  • New pivot found !

Choice of new pivots:

The pivotal submatrix must form a DAG ⇒ If there is a cycle, the entry can’t be selected.

Claire Delaplace Sparse elimination Mod p JNCF 2017 15 / 20

slide-45
SLIDE 45

Recent Improvement: Better Pivots Selection Heuristic

How to Find Remaining Pivots:

  • Choice of new pivots:

The pivotal submatrix must form a DAG ⇒ If there is a cycle, the entry can’t be selected.

Claire Delaplace Sparse elimination Mod p JNCF 2017 15 / 20

slide-46
SLIDE 46

Recent Improvement: Better Pivots Selection Heuristic

How to Find Remaining Pivots:

  • Choice of new pivots:

The pivotal submatrix must form a DAG ⇒ If there is a cycle, the entry can’t be selected.

Claire Delaplace Sparse elimination Mod p JNCF 2017 15 / 20

slide-47
SLIDE 47

Recent Improvement: Better Pivots Selection Heuristic

Our new pivot selection algorithm

  • In a Nutshell:

Find the F.-L. pivots Push the F.-L. pivots on the top of A. For all non-pivotal column j, if the upmost coefficient aij is on a non-pivotal row, select aij Perform a graph search to find more pivots

Claire Delaplace Sparse elimination Mod p JNCF 2017 16 / 20

slide-48
SLIDE 48

Recent Improvement: Better Pivots Selection Heuristic

About Graph Search

Remark

This method is parallelizable. The worst case complexity of the search is O (n|A|) (same as Wiedemann) In practice : much faster than Wiedmann

Claire Delaplace Sparse elimination Mod p JNCF 2017 17 / 20

slide-49
SLIDE 49

Recent Improvement: Better Pivots Selection Heuristic

Example

Claire Delaplace Sparse elimination Mod p JNCF 2017 18 / 20

slide-50
SLIDE 50

Recent Improvement: Better Pivots Selection Heuristic

Example

F.-L. pivots

Claire Delaplace Sparse elimination Mod p JNCF 2017 18 / 20

slide-51
SLIDE 51

Recent Improvement: Better Pivots Selection Heuristic

Example

More pivots after graph search

Claire Delaplace Sparse elimination Mod p JNCF 2017 18 / 20

slide-52
SLIDE 52

Conclusion

Conclusion

Implemented in C and C++ in the SpaSM (SPArse direct Solver Modulo p) library and publicy available at: https://github.com/cbouilla/spasm Will be implemented in LinBox

Claire Delaplace Sparse elimination Mod p JNCF 2017 19 / 20

slide-53
SLIDE 53

Conclusion

Conclusion

Implemented in C and C++ in the SpaSM (SPArse direct Solver Modulo p) library and publicy available at: https://github.com/cbouilla/spasm Will be implemented in LinBox Open questions: Can we still improve the pivots selection? Is it possible to adapt this method on matrices from CADO-NFS?

Claire Delaplace Sparse elimination Mod p JNCF 2017 19 / 20

slide-54
SLIDE 54

Conclusion

Thank you for your time !

Claire Delaplace Sparse elimination Mod p JNCF 2017 20 / 20