Communication-efficient parallel Bruhat decomposition

Communication-efficient parallel Bruhat decomposition Ioanna Evtushevskaya-Konovalova and Alexander Tiskin

arXiv:2608.07724v1 [cs.DS] 7 Aug 2026

St Petersburg State University

Abstract. The model of bulk-synchronous parallel (BSP) computation is an emerging paradigm of general-purpose parallel computing. Bruhat decomposition is an important method in numerical linear algebra, generalising ordinary LU decomposition while providing a symmetric structured way of expessing pivoting. Block-recursive algorithms have been developed in the past for various numerical linear algebra problems, including Bruhat decomposition; however, expressing recursive algorithms in the BSP model still remains a challenge. In this paper, we consider the communication and synchronisation complexity of Bruhat decomposition in the BSP model. We develop and analyze two communication-efficient algorithms based on different recursion schemes.

Keywords: bulk-synchronous parallel, Bruhat decomposition, communication, synchronisation, numerical linear algebra.

1

Introduction

Bruhat decomposition is a representation of a non-singular square matrix A as a product A = U1 P U2 , where U1 , U2 are upper-triangular matrices and P is a permutation matrix. Matrix P is uniquely determined. Bruhat decomposition is a natural generalization of the LU decomposition: indeed, if all leading principal minors of matrix A are non-singular, the matrix P becomes anti-diagonal, and then the product U1 P is precisely the lower-triangular matrix L of the decomposition A = LU . Classical algorithms for parallel computation of the LU decomposition are described e.g. in [3]. In this work, we will construct a more general algorithms for Bruhat decomposition. For convenience, we reformulate the problem slightly. Let a non-singular square matrix A be given. We want to find non-singular matrices U (uppertriangular), L (lower-triangular), and P (permutation) such that the following holds: LAU = P. Indeed, given these matrices, it is easy to compute matrices for Bruhat decomposition in its classical form. In this work, we develop and analyze two algorithms based on different recursion schemes. We describe them from the BSP model perspective, including parallelization and cost analysis.

2

Ioanna Evtushevskaya-Konovalova and Alexander Tiskin

Throughout the paper, we ignore small irregularities arising from imperfect matching of parameters. For example, when we say that p processors are partitioned into two equal groups, we assume implicitly that the groups’ sizes may differ by ±1.

2

Existing work

Malaschonok in [4] describes parallel Bruhat decomposition algorithms that rely on the generalized Bruhat decomposition. His algorithm was considered outside any particular parallel computing model. We will consider a similar algorithm in the BSP model and we will estimate this algorithm’s costs. Our second algorithm relies heavily on banded matrices, and the connection of the Bruhat decomposition to such matrices is described by Strang [7]. However, Strang does not discuss his algorithm from the standpoint of parallel computation. Generalized Bruhat decomposition was introduced by Grigoriev in [2] and Malaschonok in [4]. This decomposition preserves, in a certain sense, the rank stucture of matrix A. Definition 1. A subpermutation matrix is a matrix, where every row and every column contains at most one non-zero entry, and every non-zero entry is equal to 1. Definition 2. Given an n × n matrix A, let i1 , . . . , im , j1 , . . . , jk be the indices of its zero rows and columns. An LEU -decomposition is A = LEU , where U , L are non-singular n × n upper- (respectively lower-)triangular matrices, and E is an n × n subpermutation matrix, i1 , . . . , im , j1 , . . . , jk are the indices of zero rows and columns of matrix E, and rk(E) = rk(A). Columns i1 , . . . , im of L and rows j1 , . . . , jk of U coincide with the columns and rows of the identity matrix. Every square matrix has an LEU -decomposition. (This is proved e.g. in [4].) Various generalizations of the Bruhat decomposition have also been researched by Grigoriev [2], Dumas et al. [1]. Odeh et al. [6]. Tiskin [9] has proposed BSP algorithms for the related problems of Gaussian elimination with pairwise pivoting and QR decomposition.

3

Bulk Synchronous Parallelism

3.1

The BSP model

We consider parallel computations in the BSP (Bulk Synchronous Parallel) model. This model was proposed by Valiant [10]; it is also described e.g. by Tiskin [8]. The BSP model assumes several processors that can perform computations in parallel and exchange information with each other. Each processor has a local memory, and there is also a global memory accessible to all processors. We will denote the number of processors by p. For information exchange, a network is assumed that routes communication between pairs of processors. Parameter g is introduced as the reciprocal of network bandwidth; that is, sending/receiving

Communication-efficient parallel Bruhat decomposition

3

a unit of information takes g conventional time units. Thus, in this model, it is important to optimize not only the distribution of computation blocks but also the communication between them. The computation in the BSP model is a sequence of supersteps. At the end of each superstep, a barrier synchronization occurs, ensuring that all processors have completed the current stage of local computation and network communication, and are ready to proceed to the next stage of computation. Thus, during a superstep, a processor cannot use the results of computations by other processors obtained in the same superstep, and is guaranteed to receive them only after the synchronization at the end of that superstep. Synchronization also signals that the superstep is complete and the next one can begin. Another parameter l is introduced as the cost of synchronization. Thus, the number of supersteps (each involving a synchronization) also needs to be optimized. Thus, the system has three main parameters: p, g, l. We will use the following notation: – W (k, r) – computation cost on processor r at superstep k, – H(k, r) – communication cost on processor r at superstep k, – S – number of supersteps. The computation and communication costs at superstep k are: W (k) = max(W (k, r)) r

H(k) = max(H(k, r)) r

Thus, the total cost is calculated as: S X

(W (i) + g · H(i) + l)

i=0

In the BSP model, the computation cost must be evenly balanced. Specifically, if a sequential algorithm for a computational problem has complexity T (n), then the corresponding parallel algorithm must achieve a computational cost of T (n)/p . Subject to this condition, the communication and synchronization costs are to be minimized. It is also standard to assume the slackness condition n ≫ p, which we adopt in this paper. 3.2

BSP matrix multiplication

When computing the Bruhat decomposition and the LEU decomposition, it is necessary to compute matrix products in parallel. Matrix multiplication can be carried out in parallel as follows, it can also be computed in parallel using the algorithm developed in [5]. Let A, B be n × n matrices, C = AB, and p be the number of processors. We partition the matrices into p2/3 square blocks of size n/p1/3 × n/p1/3 . Denote the block size as d = n/p1/3 . The block of matrix A formed by the intersection of rows d(i − 1) + 1, d(i − 1) + 2, . . . , di and columns d(j − 1) + 1, d(j − 1) + 2, . . . , dj is denoted by A[i, j], and similarly for matrices B and C. Then it is clear that

4

Ioanna Evtushevskaya-Konovalova and Alexander Tiskin

C[i, j] =

1/3 p X

A[i, k]B[k, j].

k=1

Thus, it is necessary to compute A[i, k]B[k, j] for all 1 ≤ i, k, j ≤ p1/3 . Each processor is assigned a triple (i, k, j) and computes the product C[i, k, j] = A[i, k]B[k, j], writing the result to external memory. Matrix C is obtained by Pp1/3 summing the corresponding blocks: C[i, j] = k=1 C[i, k, j].
Each processor performs O (n/p1/3 )3 = O(n3 /p) elementary multiplications 1/3 when computing the product of two blocks of size n/p , and sends one com
1/3 2 puted square block of the same size i.e. O (n/p ) = O(n2 /p2/3 ). Thus, the costs for computation, communication, and synchronization are: W (n, p) = O(n3 /p), H(n, p) = O(n2 /p2/3 ), S(n, p) = O(1).

4

Block-recursive algorithm

We will solve the problem in the following form: it is required to find matrices L, U such that LAU = P. Then, by computing the inverses L−1 , U −1 , we obtain the Bruhat decomposition. We will solve the problem iteratively. At iteration k of the algorithm, we eliminate some elements of the current matrix A(k) by multiplication on the left with lower-triangular L(k) and on the right with upper-triangular U (k) . Thus, we obtain a matrix A(k+1) that will gradually become a permutation matrix. Matrix A(k) is of the form ! (k) e(k) a01 (k) (k) a10 a11 Here e(k) is a block of size m × m (where m is to be specified later.) This block is a subpermutation matrix. It may contain entirely zero rows or columns; we will call them special and handle them in a special way. Note that the blocks (k) (k) a01 and a10 are rectangular, in general not square. In the following description of an individual algorithm iteration, we will reason about square blocks of size m and 2m. For convenience, blocks of size m are denoted by lowercase letters, and blocks of size 2m by uppercase letters. Each block is processed in parallel using recursive calls on smaller blocks, down to a certain sufficiently small block size n0 . Blocks of this size will be processed sequentially.

Communication-efficient parallel Bruhat decomposition

5

Fig. 1. Block-recursive algorithm

4.1

Recursive scheme

The main idea is that when processing a certain block, we can reduce it to several matrix multiplications and the computation of an LEU -decomposition on smaller blocks. Computation will be distributed among processors, and at the bottom of the recursion tree, each processor will compute the LEU -decomposition for a separate block that is sufficiently small so that the LEU -decomposition can be computed sequentially. Our algorithm is described as a sequence of iterations. From the point of view of parallel computations, an iteration looks as follows: first, the LEU decomposition for the upper-left block of the matrix is computed recursively; then the processors receive the bottom-left and upper-right blocks and compute certain auxiliary matrices. Each of these matrices is computed using the BSP matrix multiplication algorithm from Subsection 3.2; moreover, different matrices can also be computed in parallel by different groups of processors. Then the bottom-left and upper-right blocks are updated. Next, the LEU -decomposition is computed recursively for the bottom-left and upper-right blocks. After computing the LEU -decomposition for the bottom-left and upper-right blocks and updating the lower-right block via matrix multiplications, we make one more recursive call and compute the LEU -decomposition for this block. Finally, we obtain the LEU -decomposition for the entire matrix from block LEU -decompositions, again using the BSP matrix multiplications. The algorithm is illustrated in Figure 1. Consider the process of reducing a matrix to subpermutation form using LEU -decomposition. Suppose that at some point we have obtained a matrix B of size m × m partitioned into blocks of size m/2 × m/2:   b00 b01 B= b10 b11 1. LEU -decomposition for b00 This is obtained recursively as l00 b00 u00 = e00 where e00 is a subpermutation matrix. Transformation matrices are formed as     l 0 u00 0 L1 = 00 U1 = 0 I 0 I

6

Ioanna Evtushevskaya-Konovalova and Alexander Tiskin

where I is the identity matrix of conforming size. After applying these transformations, the current matrix takes the form   e00 b′01 L1 BU1 = B1 = ′ b10 b11 2. Partial elimination of rows in b′01 and columns in b′10 After the previous step, we have some non-zero rows/columns in e00 . If e00 has a non-zero column indexed i, we can obtain a zero column indexed i in b′10 . Similarly, we can obtain a zero row indexed i in b′01 . Thus, elements in blocks b′01 and b′10 are eliminated using lower-triangular and upper-triangular matrices:    ′ I 0 Iu L′ = ′ U′ = l I 0 I where l′ = −b′10 e⊤ 00 ,

u′ = −b′01 e⊤ 00 .

After applying these transformations, the current matrix takes the form   e00 b∗01 ′ ′ L B1 U = B2 = ∗ ′ b10 b11 3. LEU -decomposition for b∗01 and b∗10 This is obtained recursively as ∗ ∗ ∗ l01 b01 u01 = e∗01

Transformation matrices are formed as ∗  l01 0 L2 = ∗ 0 l10

∗ ∗ ∗ l10 b10 u10 = e∗10

 ∗  u01 0 U2 = 0 u∗10

After applying these transformations, the current matrix takes the form   e00 e01 L2 B2 U2 = B3 = e10 b′′11 4. Partial elimination in rows and columns of block b′′11 Similarly to step 2, we obtain matrices L′′ , U ′′    ′′  I 0 Iu L′′ = ′′ U ′′ = l I 0 I l′′ = −b′′11 p⊤ 01 ,

u′′ = −b′′11 p⊤ 10

After applying these transformations, the current matrix takes the form   e e L′′ B3 U ′′ = B3′ = 00 ∗01 e10 b11

Communication-efficient parallel Bruhat decomposition

7

5. LEU -decomposition for the b∗11 This is obtained recursively as ∗ ∗ ∗ l11 b11 u11 = e∗11

with transformation matrices   I 0 L3 = ∗ 0 l11

 U3 =

I 0 0 u∗11



After applying these transformations, the current matrix takes the form   e00 e01 L3 B3 U3 = B4 = =E e10 e11 Summary Overall, the original matrix B is reduced to sub-permutation form by successively multiplying on the left by lower-triangular matrices and on the right by upper-triangular matrices: L1 L′ L2 L′′ L3 B U3 U ′′ U2 U ′ U1 = E, where E is the resulting sub-permutation matrix. This is equivalent to finding the LEU decomposition of B. If B is non-singular, than E is a permutation matrix. 4.2

Cost analysis

We can derive the following recurrence relations for our algorithm’s costs of computation, communication, and synchronization. Let m be the matrix size and q the number of processors in the current recursive call. As mentioned earlier, at the bottom of the recursion tree, each processor will process a separate block that is sufficiently small to be processed sequentially. The size of this block can be varied. On one hand, it must be small enough for uniform division of computation. On the other hand, partitioning the matrix into overly small blocks increases synchronization cost. Let n0 = n/pα be the block size, where α is a parameter to be defined later. By varying α, we select the optimal block size n0 . Obviously, the recursion depth depends linearly on α. The computation, communication and synchronization costs in this scheme can be described by the following recurrences: W (m, q) = 2W (m/2, q) + W (m/2, q/2) + O(m3 /q) H(m, q) = 2H(m/2, q) + H(m/2, q/2) + O(m2 /q 2/3 ) S(m, q) = 2S(m/2, q) + S(m/2, q/2) + O(1) Indeed, processing a block of size m with q processors proceeds as follows: first, all q processors process the upper-left block of size m/2; then the processors are split into two equal groups of q/2 processors and process the upper-right and the bottom-left blocks of size m/2 in parallel; then all q processors process the bottom-right block of size m/2. Moreover, in addition to recursive calls, matrix multiplications are needed to update the blocks. The base values for the recurrences are provided by the respective costs of sequential computation: W (n0 , q) = O(n30 )

H(n0 , q) = O(n20 )

S(n0 , q) = O(1)

8

Ioanna Evtushevskaya-Konovalova and Alexander Tiskin

λ−2

λ−1

λ

µ−2

µ−1

µ

Fig. 2. Recursion dag for the block-recursive algorithm

From the above recurrences, we can obtain explicit expressions for W (n, p), H(n, p) and S(n, p). When analyzing algorithms, we will always omit the base of the logarithm, implying the binary logarithm. Let µ = log(m), λ = log(q) and introduce new functions: W̃ (µ, λ) = W (m, q)

H̃(µ, λ) = H(m, q)

S̃(µ, λ) = S(m, q)

Let k = log(n). Then the recurrences become:  3µ  2 W̃ (µ, λ) = 2W̃ (µ − 1, λ) + W̃ (µ − 1, λ − 1) + O 2λ  2µ  2 H̃(µ, λ) = 2H̃(µ − 1, λ) + H̃(µ − 1, λ − 1) + O 22λ/3 S̃(µ, λ) = 2S̃(µ − 1, λ) + S̃(µ − 1, λ − 1) + O(1) These recurrences gives rise to a dag (directed acyclic graph) whose vertices are points in the 2D integer lattice. Their coordinates are interpreted as follows. The horizontal coordinate of a vertex equals the logarithm of the number of processors; the vertical coordinate equals the logarithm of the size of the block being processed. Each recursion level consists of vertices with a fixed vertical coordinate. Deepening the recursion corresponds to moving from the root up and to the left, as shown in Figure 2. Each vertex of the dag has two child vertices. In particular, vertex (µ, λ) has two child vertices: (µ − 1, λ), connected to the parent by a vertical edge, and (µ − 1, λ − 1), connected to the parent by a diagonal edge. All the vertical edges are doubled, while diagonal ones are not. (This corresponds to coefficients 1 and 2 in the recurrence.)

Communication-efficient parallel Bruhat decomposition

9

To express the costs W̃ (µ, λ), H̃(µ, λ), S̃(µ, λ) we write for a vertex v with coordinates (µ, λ)  2µ   3µ  2 2 h(v) = O s(v) = O(1) w(v) = O 2λ 22λ/3 r(v) = the number of paths from the root to v We have W̃ (ν, π) =

P

v w(v)r(v)

H̃(ν, π) =

P

v h(v)r(v)

S̃(ν, π) =

P

v s(v)r(v)

where the sums range across all the vertices of the recursion dag. Consider level i of this dag: vertices at this level have coordinates of the form (k − i, i − j) for all 0 ≤ j ≤ i. A path to such a vertex is of length i with exactly j vertical edges. Since vertical edges are doubled, traversing a vertical
edge doubles the number of paths. Thus, the total number of paths equals ji 2j . Next, we will separately sum up the terms corresponding to vertices of each level. First, we obtain the computation cost. Summing up the terms w(v) corresponding to vertices of level i:
P Pi W̃i = v=(k−i,j),j=0,...,i w(v)r(v) = j=0 ji 2j+3(k−i)+(i−j) p1 =

Pi Pi i 3k−2i 1 = 23k−2i p1 j=0 ji = 23k−i p1 j=0 j 2 p Then, summing up over all levels, we obtain:  3 Pk Pk n 3k −i 1 3k 1 W̃ = 2 2 = O(2 ) = O i i=0 i=0 p p p Now we similarly obtain the communication cost. Summing the terms corresponding to vertices of level i, we obtain
P Pi 1 H̃i = v=(k−i,j),j=0,...,i h(v)r(v) = j=0 ji 2j+2(k−i)+2(i−j)/3 · p2/3

Pi Pi i j/3 1 1 1 = j=0 ji 22k−4i/3+j/3 · p2/3 = 22k−4i/3 p2/3 = 22k−4i/3 p2/3 (1 + 21/3 )i j=0 j 2 Then, summing over all levels, we obtain:
i Pk Pk Pk 1 2k 1 −4i/3 −4/3 (1 + 21/3 )i = 22k p2/3 (1 + 21/3 ) ≈ i=0 H̃i = 2 p2/3 i=0 2 i=0 2  2  P∞ 1 n i 2k+3 1 ≈ 22k p2/3 = O 0.89 ≈ 2 2/3 i=0 p p2/3 Now we similarly obtain the synchronization cost:  P P P
 k i k j i i S(n, p) = O =O = O(3k ) i=0 3 i=0 j=0 2 j Overall we have obtained the following asymptotic costs:  3  W (n, p) = O pn2/3  2   2  H(n, p) = O pn2/3 + O(3k ) = O pn2/3 S(n, p) = O(3k ) = O(nlog 3 ) ≈ O(n1.58 )

10

Ioanna Evtushevskaya-Konovalova and Alexander Tiskin

Recall that the recursion tree is evaluated in parallel while the block size is large, and switches to sequential evaluation once a critical block size n0 = n/pα is reached. Denote k0 = log(n/n0 ). For every vertex u = (k0 , λ − j), we delete all its descendants and make it a leaf. We have w(u) = O(23k0 )

h(u) = O(22k0 )

s(u) = O(1)

r(u) = r(u)

We denote the set of inner vertices (vertices at levels 0 to k0 − 1) and leaves (level k0 ) by Inner , Leaf , respectively. Using the previous analysis of the recurrence dag, we can bound the contribution of the Inner terms as follows:  2   3 H̃inner ≤ O pn2/3 W̃inner ≤ O np The bound on synchronization is obtained from the previous analysis by substituting k0 for k:    P P log 3
 Pk0 i k0 i n j i k0 S̃inner = O = O( i=0 3 ) = O(3 ) = O i=0 j=0 2 j n0 Independently, we sum up the Leaf terms: P Pk0 k0
3(k−k0 ) j W̃leaf = v∈Leaf w(v)r(v) = j=0 2 j 2
P k0 k0 j 3(k−k0 ) k0 = 23(k−k0 ) j=0 3 j 2 =2  log 3 = 23(log n−log(n/n0 )) 3log(n/n0 ) = n30 nn0      3 n3 n α log 3 p = O =O α α(3−log 3) p p H̃leaf =

P

v∈Leaf h(v)r(v) =

= O(n20 )3k0 = O



n pα

Pk0

j=0

2 

k0 j

j+2(k−k ) Pk0 k0
j 0 2 = 22(k−k0 ) j=0 j 2

pα log 3

S̃leaf = O(3k0 )  3 Since we need W̃leaf = O np , it must hold α ≥ αmin = 1/(3 − log 3) ≈ 0.70. Finally, we sum up the Inner and Leaf expressions together:  3 W (n, p) = W̃inner + W̃leaf = O np H(n, p) = H̃inner + H̃leaf =  2   2 !  2    n n n n2 α log(3) =O + O p = O + O pα p2/3 p2/3 pα(2−log 3)   2  2 O n2/3 α ≥ 3(2−log 3) = αmax ≈ 1.59 p 2  H(n, p) = n O α(2−log 3) otherwise p S(n, p) = O(3k0 )

Communication-efficient parallel Bruhat decomposition

11

Thus, we only need to consider αmin ≤ α ≤ αmax , since α ≥ αmin is needed for the optimal local computation cost, and for α > αmax , the communication cost of matrix multiplication dominates, while the syncronization cost is higher than that of α = αmax . The parameter α describes a trade-off between communication and synchronization costs.

5

Strip-recursive algorithm

In this section, we describe an algorithm based on a different recursion scheme, which we call antidiagonal recursion. We partition the matrix into antidiagonal strips. The algorithm proceeds by splitting a strip into two half-sized strips and processing them recursively; between the two recursive calls, the matrix needs to be updated in a carefully controlled manner in order to optimize communication. 5.1

Banded Matrices

The algorithm operates on matrices of aim special form. Definition 3. Matrix A is called lower banded of bandwidth k if Aij ̸= 0 ⇒ 0 ≤ i − j ≤ k. An upper banded matrix is defined symmetrically. Lemma 1. Let matrices A, B be lower banded of bandwidth k. Then the matrix AB is lower banded of bandwidth 2k. Lemma 2. Let A, B be n × n matrices, lower banded of bandwidth m. Then the product AB can be computed in parallel on p processors with the following costs:  4/3 2/3   2  n H = O m p2/3 , S = O(1). W = O mp n , Proof. Each of the matrices contains O(mn) nonzero elements. Each element of matrix A is multiplied by O(m) elements of matrix B, so the computation contains O(m2 n) elementary multiplications. Our goal is to partition these operations into equal near-cubic cuboid blocks. We do this as follows: partition the nonzero part of matrix A (i.e., the strip of width m below main diagonal) into strips parallel to main diagonal of width O(m2/3 n1/3 p1/3 ); then partition each strip horizontally into regular parallelograms of equal width (Figure 3). Proceed similarly with the second matrix, but partition the strips into parallelograms vertically. The product of each pair of blocks requires O(m2 n) operations. After computing the product of a pair of blocks, O((m2/3 n1/3 )2 ) = O(m4/3 n2/3 ) values are obtained, which must be transferred between processors to compute the values of the product matrix. Thus, multiplying matrices of this type requires computation cost W = O(m2 n/p) and communication cost H = O(m4/3 n2/3 /p2/3 ). The number of supersteps is obviously constant: S = O(1).

12

Ioanna Evtushevskaya-Konovalova and Alexander Tiskin m2/3 n1/3 p1/3

m2/3 n1/3 p1/3

Fig. 3. Multiplication of lower banded matrices

5.2

Antidiagonal strips

Next, we consider an algorithm based on the following recursive scheme. The matrix is partitioned into anti-diagonal strips of width n0 = n/pα . The algorithm processes these strips sequentially, from left-to-right and top-down. In lemma 3 we describe parallel algorithm for reducing in BSP model. We call a strip reduced if every row (respectively, column) of the matrix contains at most one nonzero entry within the strip, and that entry is 1. A strip of width n0 can be reduced by multiplying the matrix on the left and right by banded matrices of bandwidth n0 . Before reduction, preprocessing of the current strip is required in order to eliminate certain rows and columns. This elimination is performed by multiplying on the left and right by banded matrices of bandwidth 2n0 . Lemma 3. An antidiagonal strip of width n0 of an n × n matrix can be reduced by multiplying the matrix on the left and right by banded matrices, which are also of size n × n and bandwidth n0 . This reduction can be performed in the BSP model with computation cost W = O(n30 ) and communication cost H = O(n20 ). Proof. To reduce the strip, partition it horizontally into blocks of n0 rows each (there will be between 1 and p blocks). Then each processor processes one block assigned to it. We perform the reduction using a sequence of transvections on each row and each column with rows below and columns to the right within the same block. Now each block contains at most one nonzero element in each row and each column. Each column intersects at most two blocks, hence contains at most two nonzero elements. Then each processor sends the resulting block to the processor handling the next block, and if necessary, row transvection is performed to eliminate the second nonzero element in the column. Note that each row is subtracted from at most n0 rows below, and each column from at most n0 columns to the right; therefore, the transvection matrices contain nonzero elements only in the strip of width n0 below (above) the main diagonal. Thus, reduction of the strip is accomplished by multiplying on the left and right by lower-banded and upper-banded matrices of bandwidth n0 , respectively. Processing a block of the indicated size requires computation cost O(n30 ) and communication cost O(n20 ); the communication and computation costs of the ⊔ ⊓ second stage are O(n20 ) and O(n0 ), respectively.

Communication-efficient parallel Bruhat decomposition

13

We introduce a numbering of the antidiagonals from left to right and top to bottom. Antidiagonal r consists of elements Ai,j , where i + j − 1 = r. The strip consisting of consecutive antidiagonals numbered m, m + 1, . . . , m + k inclusive will be denoted by A⟨m, m + k⟩. Lemma 4. Suppose that strip A⟨k, k + m⟩ has been reduced and let i1 , . . . , ik , j1 , . . . , jk be the indices of rows and columns, containing nonzero elements in this strip. By multiplying A on the left by a lower-banded matrix L and on the right by an upper-banded matrix U , both of bandwidth 2m, we can obtain a matrix A′ , where A′ ⟨k + m, k + 2m⟩ consists only of zero elements in rows i1 , . . . , ik and columns j1 , . . . , jk . Proof. Let us construct the elimination matrix for the columns (the elimination matrix for the rows is constructed similarly). For all i, j, where Aij = 1 and k + 1 ≤ i + j ≤ k + m + 1, i.e., the element Aij belongs to strip A⟨k, k + m⟩, let Lr,i = −Ar,j for k ( + m + 1 − j ≤ r ≤ k + 2m + 1 − j. Other elements of L are 0, i ̸= j defined as Li,j = . Direct computation verifies that matrix L does 1, i = j indeed eliminate the required columns. It is also clear that it is lower-banded of bandwidth 2m. ⊔ ⊓ 5.3

Recursive scheme

The algorithm takes strip A⟨k, k + m⟩ as input and reduces it as follows: 1. Reduction of A⟨k, k + m/2⟩ This is performed recursively, thus obtaining elimination matrices Lm and Um of bandwidth 2m; Lm AUm = A1 2. Update on A1 ⟨k + m/2 + 1, m⟩ Thus obtaining elimination matrices L′m and ′ Um , with which we eliminate certain rows and columns in the lower-right strip: more precisely, those that already contain a nonzero element in A⟨k, k + m/2⟩. ′ are lower-banded and upper-banded of bandwidth 2m, Matrices L′m and Um ′ respectively; Lm A1 Um = A2 3. Reduction of A2 ⟨k +m/2+1, m⟩ This is performed recursively, thus obtaining elimination matrices L̃m and Ũm (again a lower-banded and an upper-banded matrix of bandwidth 2m); L̃m A2 Ũm = AF inal Summary Overall, A⟨k, k + m⟩ is reduced by multiplying on the left by lowertriangular elimination matrices and on the right by upper-triangular elimination ′ matrices L̃m L′m Lm AUm Um Ũm = AF inal , where AF inal ⟨k, k + m⟩ is reduced. 5.4

Cost analysis

A recursive call on a strip of bandwidth m performs two recursive calls on strips of bandwidth m/2 and several banded matrix products, as described in steps 2

14

Ioanna Evtushevskaya-Konovalova and Alexander Tiskin

and 4. These matrix products are performed in parallel similarly to the proof of Lemma 4 with the following costs:  4/3 2/3   2  n , S = O(1). W = O mp n , H = O m p2/3 ′ Next, we need to compute L2m = Lm L′m Lm and U2m = Um Um Um ; this is done with the same costs, according to Lemma 4. Thus, we can derive the following recurrence relations:  2  W (m, p) = 2 W (m/2, p) + O mp n ,  4/3 2/3  n H(m, p) = 2 H(m/2, p) + O m p2/3 ,

S(m, p) = 2 S(m/2, p) + O(1). The base values for the recurrences are provided by the respective costs of reducing a strip of width n0 , as described in Lemma 3: W (n0 , q) = O(n30 )

H(n0 , q) = O(n20 )

S(n0 , q) = O(1)

Since the recursion depth is α log p, and n0 = n/pα , the recurrence relations are solved as follows:  3  3  3 Pα log p  i 2  W (n, p) = i=0 O (2 np0 ) n + 2α log p O(n0 3 ) = O np + O pn2α = O np , Pα log p  i 43 23  H(n, p) = i=0 O (2 n0 )2 n + 2α log p O(n0 2 ) p3

=

αX log p

4i

23 O

 n 34 n 23 

i=0

=O

0

p

2 3

 2 n

2

p3

4α

+ 2α log p O(n0 2 ) = p 3 O



n2 2 4α p3+ 3



+ pα O



n2 p2α



  2 O n2 α ≥ αmax = 23  p 23  + O( pnα ) = O nα otherwise p 2

S(n, p) = O(2α log p ) = O(pα ).  3 Since we need W (n, p) = O np it must hold α ≥ αmin = 1/2. Similarly to the block-recursive algorithm, parameter α describes a trade-off between communication and synchronization costs. Thus, we only need to consider αmin ≤ α ≤ αmax , since α ≥ αmin is needed for the optimal local computation cost, and for α > αmax = 32 , the communication cost of matrix multiplication dominates, while the synchronization cost is higher than that of α = αmax .

6

Conclusion

We have developed two algorithms for Bruhat decomposition in the BSP model: the block-recursive algorithm and the strip-recursive algorithm, and compared

Communication-efficient parallel Bruhat decomposition

15

their effectiveness. The block-recursive algorithm uses a more familiar partition into blocks and a simpler recursive scheme. The strip-recursive algorithm requires some technical preparation (introducing special types of matrices and adapting existing parallel algorithms to such matrices). On the positive side, this algorithm uses a simpler recursion scheme and has lower cost. Both algorithms rely on efficient BSP matrix multiplication, and exhibit a trade-off between communication and synchronization costs. In the strip-recursive algorithm this trade-off is more efficient than in the block-recursive algorithm, matching the analogous trade-off by Tiskin [9] for LU -decomposition and generic pairwise elimination.

References 1. Jean-Guillaume Dumas, Clément Pernet, and Ziad Sultan. Fast computation of the rank profile matrix and the generalized Bruhat decomposition. Journal of Symbolic Computation, 83:187–210, 2017. 2. D.Yu. Grigor’ev. Additive complexity in directed computations. Theoretical Computer Science, 19(1):39–67, 1982. 3. N J Higham. Accuracy and Stability of Numerical Algorithms. SIAM, second edition, 2002. 4. Gennadi Malaschonok. Fast Generalized Bruhat Decomposition. In Proceedings of CASC, volume 6244 of Lecture Notes in Computer Science, pages 194–202, 2010. 5. W. F. McColl and A. Tiskin. Memory-Efficient Matrix Multiplication in the BSP Model. Algorithmica, 24(3-4):287–297, 1999. 6. O. H. Odeh, D. D. Olesky, and P. van den Driessche. Bruhat Decomposition and Numerical Stability. SIAM Journal on Matrix Analysis and Applications, 19(1):89– 98, jan 1998. 7. Gilbert Strang. The Algebra of Elimination. In Excursions in Harmonic Analysis, volume 3, pages 3–22. 2015. 8. Alexander Tiskin. The bulk-synchronous parallel random access machine. Theoretical Computer Science, 196(1-2):109–130, 1998. 9. Alexander Tiskin. Communication-efficient parallel generic pairwise elimination. Future Generation Computer Systems, 23:179–188, 2007. 10. L G Valiant. A Bridging Model for Parallel Computation. Communications of the ACM, 33(8):103–111, 1990.

Related documents

• FastPair: GPU-Optimized String Decoding#1084928
• FastPair: GPU-Optimized String Decoding#919348
• Terastate-per-second QUBO Brute-Force on a Single GPU: A Matrix Prefix-Suffix Decomposition#343428
• Fractional Verkle Trees: A Hypertree Decomposition and Verified Proof Serialization Architecture for High-Performance Blockchain State Accumulators#282772
• A first introduction to Matrix Product State algorithms for the integration of Lindblad equation#1012999