Recursive Matrix Algorithms in Commutative Domain for Cluster with Distributed Memory

TL;DR

This paper reviews theoretical matrix algorithms in commutative domains, presents experimental results on parallel implementations on a supercomputer, and demonstrates their scalability and efficiency for distributed memory clusters.

Contribution

It introduces new parallel matrix algorithms in commutative domains and evaluates their scalability and efficiency on a supercomputer cluster.

Findings

Algorithms achieve high efficiency on distributed memory clusters

Parallel programs scale well with increasing processors

Applications include various computational tasks in algebra and matrix computations

Abstract

We give an overview of the theoretical results for matrix block-recursive algorithms in commutative domains and present the results of experiments that we conducted with new parallel programs based on these algorithms on a supercomputer MVS-10P at the Joint Supercomputer Center of the Russian Academy of Science. To demonstrate a scalability of these programs we measure the running time of the program for a different number of processors and plot the graphs of efficiency factor. Also we present the main application areas in which such parallel algorithms are used. It is concluded that this class of algorithms allows to obtain efficient parallel programs on clusters with distributed memory.

Equations

$$\left(\begin{array}{cc}{A}_0& A_1\\ A_2& A_3\end{array}\right)\times \left(\begin{array}{cc}{B}_0& B_1\\ B_2& B_3\end{array}\right)+\left(\begin{array}{cc}{C}_0& C_1\\ C_2& C_3\end{array}\right)=\left(\begin{array}{cc}{D}_0& D_1\\ D_2& D_3\end{array}\right)$$

$${\mathcal{A}}^{-1}=\left(\begin{array}{cc}\text{I}& -A_0^{-1}{A}_1\\ 0& \text{I}\end{array}\right)\left(\begin{array}{cc}\text{I}& 0\\ 0& (A_3-A_2{A}_0^{-1}{A}_1{)}^{-1}\end{array}\right)\times$$

$$\left(\begin{array}{cc}\text{I}& 0\\ -A_2& \text{I}\end{array}\right)\left(\begin{array}{cc}{A}_0^{-1}& 0\\ 0& \text{I}\end{array}\right)=\left(\begin{array}{cc}{M}_1{M}_5-M_0& M_1{M}_4\\ M_5& M_4\end{array}\right).$$

$${\mathcal{A}}^{-1}=\left(\begin{array}{cc}{A}^{-1}& 0\\ -C^{-1}B{A}^{-1}& C^{-1}\end{array}\right).$$

$$Chol:R^{n\times n}\to (R^{n\times n},R^{n\times n}),$$

$$Chol(\mathcal{A})=(H,H^{-1})$$

$$I_{ij}=E_{ij}{E}_{ij}^T,{\bar{I}}_{ij}=\text{I}-I_{ij},Y_{ij}=E_{ij}^T{S}_{ij}-d_{ij}\text{I},i,j\in 1,2.$$

$$A_{ext}:R^{n\times n}\times (R\backslash 0)\to (R^{n\times n}{)}^3\times (R\backslash 0),$$

$$A_{ext}(M,d_0)=(A,S,E_S,d)$$

$$((\frac{S_{11}{E}_{21}^T{A}_{21}}{d_{11}}{M}_{22}^1-I_{11}{M}_{12}^1{d}_{21})/d_{11}\ast Y_{12}+S_{12}{d}_{21})/d_{11},$$

$$A^1=A_{12}{A}_{11},L=(A^1-(I_{11}{M}_{12}^1{E}_{12}^T{A}^1)/d_{11})/d_{11}\ast d_{22},$$

$$A^2=A_{22}{A}_{21},P=(A^2-(I_{21}{M}_{22}^2{E}_{22}^T{A}^2)/d_s)/d_{21},$$

$$A=\left(\begin{array}{cc}(L+FG)/d_{12}& F\\ (PG)/d_{12}& P\end{array}\right),S=\left(\begin{array}{cc}(M_{11}^2{d}_{22})/d_{21}& M_{12}^3\\ (S_{21}{d}_{22})/d_{21}& M_{22}^3\end{array}\right),$$

$$Then{A}_{ext}(M,d_0)=(A,S,E,d).$$

$$t_{n}n=const.$$

$$f=\frac{t_{n}n}{t_{k}k}100\%$$

Full text

I Introduction

J. Dongarra at his talk at International Congress ICMS-2016 [1] put attansion on the several difficult challenges. He noted that the task of managing calculations on a cluster with distributed memory for algorithms with sparse matrices is today one of these the most difficult challenges.

We have to add also one more problem. It is a high computational complexity, that can be connected with the type of the basic algebra: you can take a matrix over field or over commutative ring.

For sparse matrices, it is not true that all computations over polynomials or integers can be effectivelly reduced due to the technic of modular computations. It was proved in theoretical investigatios of computational complexity for some algorithms with sparse matrices in [2]. Below we give experiments with large sparse matrices, which confirm these theoretical reults.

We consider the class of block-recursive matrix algorithms. The most famous of them are standard and Strassen’s block matrix multiplication, Strassen’s block-matrix inversion [3].

Block-recursive algorithms were not so important as long as the calculations were performed on computers with shared memory. Only in the nineties it became clear that block-recursive matrix algorithms are required to operate with sparse large matrices on a supercomputer with distributed memory.

Note that the generalization of Strassen’s matrix inversion algorithm [3] with additional permutations of rows and columns by J. Bunch and J. Hopkroft [4] is not a block-recursive algorithm.

The block recursive algorithm for the solution of systems of linear equations and for adjoint matrix computation which is some generalisation of Strassen’s inversion in commutative domains was suggested in the papers [8], [9] and [11]. See also at the book [10]. However, in all these algorithms, except matrix multiplication, a very strong restriction are imposed on the matrix: the leading corner minors should not be zero.

This restriction was removed later. The algorithm that computes the adjoint matrix, the echelon form, and the kernel of the matrix operator for the commutative domains was proposed in [12]. The block-recursive algorithm for the Bruhat decomposition and the LEU decomposition for the matrix over the field was obtained in [13], and these algorithms were generalized for the matrices over commutative domains in [15] and in [16].

In this article we review the main achievements in this class of algorithms and present the results of experiments that we conducted with these algorithms on a supercomputer MVS-10P at the Joint Supercomputer Center of the Russian Academy of Science.

In the next section, we present the main application areas in which such algorithms are used.

II Some important areas for applications of sparse matrices algorithms

II-A * Computations of functions of electronic circuits*

The behavior of electronic circuits can be described by Kirchhoff’s laws. The three basic approaches in this theory are direct current, constant frequency current and a current that varies with time. All these cases require the compilation and solution of sparse systems of equations (numerical, polynomial or differential). The solution of such differential equations by the Laplace method also leads to the solution of polynomial systems of equations [17].

II-B * Control systems*

In 1967 Howard H. Rosenbrock introduced a useful state-space representation and transfer function matrix form for control systems, which is known as the Rosenbrock System Matrix [18]. Since that time, the properties of the matrix of polynomials being intensively studied in the literature of linear control systems.

II-C Computation of Gröbner bases

Another important application is the calculation of Gröbner bases. A matrix composed of Buchberger S-polynomials is a strongly sparse matrix. Reduction of the polynomial system is performed when calculating the echelon and diagonal forms of this matrix. The algorithm F4 [19] was the first such matrix algorithm.

II-D Solving ODE’s and PDE’s.

Solving ODE’s and PDE’s is often based on solution of leanear systems with sparse matrices over numbers or over polynomials. One of the important class of sparse matrix is called quasiseparable. Any submatrix of quasiseparable matrix entirely below or above the main diagonal has small rank. These quasiseparable matrices arise naturally in solving PDE’s for particle interaction with the Fast Multi-pole Method (FMM). The efficiency of application of the block-recursive algorithm of the Bruhat decomposition to the quasiseparable matrices is studied in [21].

III Development of the recursive matrix agorithms in integral domain

We can trace how developed the matrix recursive agorithms in integral domain, which eventually led to the creation of modern algorithms. There are several separate periods.

III-A Algorithms for solution of a system of linear equations of size $n$ in an integral domain, which served as the basis for recursive algorithms

(1983) Forward and backward algorithm ($\sim n^{3}$) [5].

(1989) One pass algorithm ($\sim\frac{2}{3}n^{3}$) [6].

(1995) Combined algoritm with upper left block of size $r$ ($\sim\frac{7}{12}n^{3}$ for $r=\frac{n}{2}$) [7].

Really, this was already the first step of a recursive algorithm. It was first discovered that when the matrix is divided into equal four blocks ($r=\frac{n}{2}$), the least computational complexity is achieved. Consequently, further dichotomous division of blocks can give the best algorithm. It remains to prove several determinant identities that would allow us to do recursive calculations.

III-B * Recursive algorithms for solution of a system of linear equations and for

adjoint matrix computation in an integral domain without permutations *

(1997) Recursive algorithm for solution of a system of linear equations [8].

(2000) Adjoint matrix computation (with 6 levels) [9].

(2006) Adjoint matrix computation alternative algorithm (with 5 levels) [11].

Now it remained to solve the problem of permutation of blocks and ensure the fulfillment of determinant identities.

III-C * Main recursive algorithms for matrices in a domain*

(2008) Computation of adjoint and inverse matrices and the operator kernel in a domain [12].

(2010) Bruhat and LEU decompositions in a feild [13].

(2012) Bruhat and LDU decompositions in a domain [14], [15].

Recursive algorithms for sparse matrices in commutative domains with the complexity of matrix multiplication are obtained. The complexity of computing the matrix product for matrices of size $n$ we denote by $\sim n^{\beta}$.

III-D New achivements and new applications

(2013) It is proved that the LEU algorithm has the complexity $O(n^{2}r^{\beta-2})$ for rank $r$ matrices [20].

(2015) New algorithms for Bruhat and LDU decompositions in a domain (alternative algorithm) [16].

(2017) It is proved that the LEU algorithm has the complexity $O(n^{2}s^{\beta-2})$ for quasiseparable matrix. A matrix is called quasiseparable if any it’s submatrix which entirely disposed below or above the main diagonal has small rank $s$, $s<<n$ [21].

IV Recursive standard and strassen’s matrix multiplication

The graph of recursive algorithm for standard matrix multiplication is shown at Figure 1.

[TABLE]

Number of operations for the standard algorithm is $\sim n^{3}$.

The graph of the Strassen multiplication algorithm can be easily represented similarly. The number of operations for this algorithm is $\sim n^{\log_{2}7}$.

The algorithm for multiplying matrices on leaf tops should take into account the sparse matrix structure and compact storage form.

We note that there exists a boundary with respect to the density of the matrix, which separates the region of applicability of the Strassen multiplication. We note that there exists a theoretical boundary for the density of sparse matrix, which separates the region of efficient application of the Strassen algorithm of multiplication. If the density of the matrix is below this boundary, then only standard multiplication is effective. This is due to the fact that the addition of blocks, which is performed in the Strassen algorithm, leads to an increase in the density of the matrix blocks (see details in [2]).

V Recursive Strassen’s matrix inversion

If $\mathcal{A}=\begin{pmatrix}A_{0}&A_{1}\\ A_{2}&A_{3}\end{pmatrix}$, $\det({\mathcal{A}})\neq 0$ and $\det(A_{0})\neq 0$ then

[TABLE]

[TABLE]

We have denoted here $M_{0}=-A_{0}^{-1}$,  $M_{1}=M_{0}A_{1}$,  $M_{2}=A_{2}M_{0}$, $M_{3}=M_{2}A_{1}$,  $M_{4}=(A_{3}+M_{3})^{-1}$,  $M_{5}=-M_{4}M_{2}$.

The graph of recursive Strassen’s matrix inversion is shown on Figure 2.

VI Recursive inversion of triangular matrix

If $\mathcal{A}=\begin{pmatrix}A&0\\ B&C\end{pmatrix}$ is triangular matrix of order $2^{k}$ and $\det({\mathcal{A}})\neq 0$ then

[TABLE]

VII Recursive Cholesky decomposition

Let $\mathcal{A}=\begin{pmatrix}A_{1}&A_{2}\\ A_{2}^{T}&A_{3}\end{pmatrix}$ be a positive definite symmetric matrix and $H=\begin{pmatrix}B&0\\ C&D\end{pmatrix}$ be a low triangle matrix with the property ${\mathcal{A}}=HH^{T}$. The mapping

[TABLE]

[TABLE]

is called an Cholesky decomposition. Let $n=2^{k}$, then you can used following recursive algorithm.

VII-1

Let $Chol({A_{1}})=(B,B^{-1})$.

We can compute C=$A_{2}^{T}(B^{-1})^{T}$ and $F=A_{3}-CC^{T}$

VII-2

Let $Chol({F})=(D,D^{-1})$.

Then $H=\begin{pmatrix}B&0\\ C&D\end{pmatrix}$ and $H^{-1}=\begin{pmatrix}B^{-1}&0\\ -D^{-1}CB^{-1}&D^{-1}\end{pmatrix}$.

VIII Recursive computation of the adjoint matrix, kernel and determinant

We consider matrices over a commutative domain.

Semigroup $P_{n}$ is formed by $n\times n$ matrices, which have the number of unit elements coincides with its rank, and the remaining elements are zero. The semigroup $D_{n}$ is formed by the diagonal matrices: $D_{n}\subset P_{n}$, $|D_{n}|$ = $2^{n}$. The identity matrix ${\mathbf{I}}$ is a unit in $D_{n}$ and in $P_{n}$.  For each matrix $E\in P_{n}$ we define diagonal matrices $I_{E}=EE^{T}\in D_{n},J_{E}=E^{T}E\in D_{n}.$ Also, we used the involution fanction on $D_{n}:{\bar{I}}={\mathbf{I}}-I$, $\forall I\in D_{n}$ with the property ${\bar{\bar{I}}}=I$.

For the matrix $E$, the matrix ${\bar{I}_{E}}$ is a left annihilator, and the matrix ${\bar{J}_{E}}$ is a right annihilator.  So we can denote the set of echelon matrix of order $n$: $S_{n}=\{S\ |\exists E\in P_{n},\ \exists d\in R\setminus 0:\ S=I_{E}S,\ dE=SJ_{E}\}.$ In other words, $dE$ ($E\in P_{n},\ d\in R\setminus 0$) is a block of echelon matrix $S\in S_{n}$ with ${\rm rank}\ E={\rm rank}\ S$, such that the sets of zero rows of the matrices $S$ and $E$ coincide, and each nonzero column of the matrix $dE$ coincides with the same column of the matrix $S$. We write: $E=E_{S},S\in S_{n}$.

Below we will use such notation for any matrix $S_{ij}\in S_{n}$ and $E_{ij}=(E_{ij})_{S_{ij}}$:

[TABLE]

The mapping

[TABLE]

[TABLE]

for $n=2^{k}$ is called an extended adjoint mapping of the pair ($M,d_{0}$) if it is defined recursively as follows.

For $M=0$ we define $A_{ext}(0,d_{0})=(d_{0}I,0,0,d_{0}).$

For $k=0$ and $M=a\neq 0$ we define $A_{ext}(a,d_{0})=(d_{0},a,a,a).$

In all other cases, we split the matrix $M$ into four equal blocks $M=\begin{pmatrix}M_{11}&M_{12}\\ M_{21}&M_{22}\end{pmatrix}$

VIII-1

Let $A_{ext}(M_{11},d_{0})=(A_{11},S_{11},E_{11},d_{11}).$

We denote $M^{1}_{12}={A_{11}M_{12}}/{d_{0}}$, $M^{1}_{21}=-{M_{21}Y_{11}}/{d_{0}}$, $M^{1}_{22}={M_{22}d_{11}-M_{21}E^{T}_{11}M^{1}_{12}}/{d_{0}}$.

VIII-2

Let $A_{ext}(\bar{I}_{11}M^{1}_{12},d_{11})=(A_{12},S_{12},E_{12},d_{12})$.

VIII-3

Let $A_{ext}(M^{1}_{21},d_{11})=(A_{21},S_{21},E_{21},d_{21}).$

We denote $M^{2}_{22}=-{A_{21}M^{1}_{22}Y_{12}}/{(d_{11})^{2}}$, $d_{s}={d_{21}d_{12}}/{d_{11}}.$

VIII-4

Let $A_{ext}(\bar{I}_{21}M^{2}_{22},d_{s})=(A_{22},S_{22},E_{22},d_{22}).$

We denote $M^{2}_{11}=-{S_{11}Y_{21}}/{d_{11}}$, $M^{2}_{12}=$

[TABLE]

$M^{3}_{12}=-{M^{2}_{12}Y_{22}}/{d_{s}}$, $M^{3}_{22}=S_{22}-{I_{21}M^{2}_{22}Y_{22}}/{d_{s}}$,

[TABLE]

[TABLE]

$F=-({({S_{11}E^{T}_{21}A_{21}})/{d_{11}}*d_{22}+({M^{2}_{12}E^{T}_{22}A^{2}})/{d_{s}}})/{d_{21}}$,    $G=-({({M_{21}E^{T}_{11}A_{11}})/{d_{0}}*d_{12}+({M^{1}_{22}E^{T}_{12}A^{1}})/{d_{11}}})/{d_{11}},$

[TABLE]

$E=\begin{pmatrix}E_{11}&E_{12}\\ E_{21}&E_{22}\end{pmatrix},\ d=d_{22}.$

[TABLE]

Sentence. The map $A_{ext}(M,1)=(A,S,E,d)$ defines an extended adjoint matrix $A$, an echelon matrix $S$, and a matrix $E_{S}$ such that $AM=S$ and $dE=SJ_{E}$ [12]. The graph of extended adjoint map is shown at Figure 3.

IX Results of experiments with matrix recursive algorithms on a cluster with distributed memory

The block-recursive matrix algorithms require a special approachs to managing parallel programs. One approach to the cluster computations management is a scheme with one dispatcher.

We consider another scheme of cluster menagement. It is a scheme with multidispatching, when each involved computing module has its own dispatch thread and several processing threads. Each processor, along with its subtask, receives a list of slave processors. During the computation, this list changes when new processors are added or when some of these processors complete their subtasks [22], [23].

IX-A Scalability

We have done the experiments using the supercomputer MVS-10P based on RSC Tornado architecture. It is a 10 Petaflops supercomputing system at the Joint Supercomputer Center (JSCC) of the Russian Academy of Science (RAS).

We demonstrate the results of experiments with parallel programms on the base of a scheme with multidispatching, using C++ and OpenMPI. We demonstrate a scalability of these programs. To do this, we plot the graph of efficiency factor.

For an “ideal parallel program” the product of the time $t_{n}$ for solving the computational problem by the number $n$ of cores in the computational cluster must remain constant:

[TABLE]

So the value

[TABLE]

can be taken as the “efficiency factor” of the n-cores with respect to the k-cores.

In order to investigate how the efficiency factor changes with increasing number of cores for a given program, it is possible to conduct a series of experiments with different number of cores in the cluster. If the program has this coefficient above 50% for some range of number of cores, then we conside that it has good scalability in this range. We conducted several series of experiments and obtained graphs that show how the efficiency factor varies. In all experiments, except Strassen’s matrix inversion, we took matrices with 15 bit integers.

For the algorithm of recursive Strassen’s matrix inversion, we took a dense matrix with double-precision numbers. Figures 4 and 5 show the results of experiments with dence matrices of sizes 8000x8000 and 16000x16000, correspondingly. For a cluster having 200 cores, the efficiency coefficient is equal to 51% for a matrix size of 8000x8000 and 73% for a matrix size of 16000x16000.

In Figure 6, the efficiency is shown for a recursive algorithm for computing the adjoint matrix and the kernel, when the number of cores in a cluster changes from 8 to 400. We took arbitrary dense matrices with size 8000x8000. For a cluster having 200 cores, the efficiency coefficient is equal to 66%. Efficiency coefficient drops to 44% when the number of cores reaches 400.

Figures 4– 7 demonstrate that the larger the size of the matrix, the better the efficiency coefficient remains with increasing cluster sizes.

IX-B Comparison of the calculation time

The next three figures show how the calculation time varies with   the growth of the matrix sizes for different algorithms and for different sparseness of the matrices.

We investigate the algorithm of computing the adjoint matrix and the kernel (see the algorithm in Figure 3). Experiments with two different versions of this algorithm are compared in Figures 7 and 8. These experiments were carried out for integer matrices on a cluster with 100 cores. We compared a standard algorithm and an algorithm in which the Chinese remainder theorem (CRT) was applied. It is well known that the application of the Chinese remainder theorem makes it possible to reduce the number of operations in such algorithms by $\sim n$ times, where $n$ is the size of the matrix, if the standard algorithm for multiplying integers is used. This is true for dense matrices.

Figure 7 shows the results of experiments for dense matrices. For a matrix of size 2500x2500, the CRT algorithm is about twice as fast as a standard algorithm.

Figure 8 shows the results of experiments for sparse matrices that have a density of 1 percent. For a matrix of size 2500x2500, the CRT algorithm is approximately twice as slow as the standard algorithm.

We see that for very sparse matrices, the CRT algorithm should not be used.

X Results of experiments with sequential programs: comparison with Mathematica and MAPLE

An experimental comparison of the sequential recursive algorithm for computing the adjoint matrix with similar programs in Mathematica and MAPLE systems is shown in Figure 9.

We compared our algorithms with Mathematica 11 and MAPLE 2015. For comparison, we took random dense integer matrices and performed calculations with identical matrices in 4 programs. The best time of calculations was demonstrated by sequential recursive algorithm. For example, for 600x600 matrices, it is twice as fast as Mathematica 11 and 7 times faster than MAPLE 2015. A slightly worse calculation time was shown by the fourth program. This is the parallel program based on recursive algorithm that was run on a single processor.

XI Conclusion

The algorithms we discussed were used in the cloud computing system of the computer algebra Math Partner. You can used this cloud system on the websites http://mathpar.cloud.unihub.ru.

Functions that correspond to these algorithms can be called by the operators adjoint, kernel, det, LDU, BruhatDecomposition.

The algorithms that have been discussed above have a wide application. Therefore, it would be important to have such a public site where users could perform parallel computing to solve their specific tasks using a large number of processors.

XII Acknowledgment

The authors are grateful to the Ivannikov Institute for System Programming of the Russian Academy of Science for hosting the cloud computer algebra system Math Partner and to the Joint Supercomputer Center of the Russian Academy of Science   for providing the ability to perform calculations on the supercomputer MVS-10P.