How to study the decomposition of a square matrix into a product of sparse matrices? 


There are no restrictions on the number of matrices in the product, but the fewer the better.



Decomposing a matrix into a product of sparse matrices
Given an invertible $n \times n$ matrix $\mathrm A$, we perform [Gaussian elimination](https://en.wikipedia.org/wiki/Gaussian_elimination) until we obtain a (nonsingular) diagonal matrix. In other words, we left-multiply $\mathrm A$ by [permutation matrices](https://en.wikipedia.org/wiki/Permutation_matrix) (whose inverses are their transposes) and by [elementary matrices](https://en.wikipedia.org/wiki/Elementary_matrix) of the form


$$\mathrm E\_{ij} := \mathrm I\_n + \beta\_{ij} \, \mathrm e\_i \mathrm e\_j^{\top}$$


whose inverse is simply


$$\mathrm E\_{ij}^{-1} = \mathrm I\_n - \beta\_{ij} \, \mathrm e\_i \mathrm e\_j^{\top}$$


until we obtain a (nonsingular) diagonal matrix $\mathrm D$, i.e.,


$$\mathrm E\_{i\_m j\_m} \mathrm P\_m \cdots \mathrm E\_{i\_2 j\_2} \mathrm P\_2 \mathrm E\_{i\_1 j\_1} \mathrm P\_1 \mathrm A = \mathrm D$$


where $m$, the number of $\mathrm E\_{i j}$ elementary matrices to the left of $\mathrm A$, is at most $n (n-1)$. Thus, one decomposition of $\mathrm A$ into a product of sparse matrices is


$$\boxed{\mathrm A = \mathrm P\_1^{\top} \mathrm E\_{i\_1 j\_1}^{-1} \mathrm P\_2^{\top} \mathrm E\_{i\_2 j\_2}^{-1} \cdots \mathrm P\_m^{\top} \mathrm E\_{i\_m j\_m}^{-1} \mathrm D}$$


where the $\mathrm E\_{i j}^{-1}$ matrices are *sparse* and easy to compute. The fraction of nonzero entries in these matrices is $n^{-1} + n^{-2} \approx n^{-1}$ for the $\mathrm E\_{i j}^{-1}$ matrices and $n^{-1}$ for matrix $\mathrm D$ and for each permutation matrix. Hence, the larger $n$, the sparser the matrices.


Multiplying the sparse elementary matrices pertaining to the same column (e.g., $\mathrm E\_{21}, \mathrm E\_{31}, \dots, \mathrm E\_{n1}$ for the first column), we obtain less sparse [Frobenius matrices](https://en.wikipedia.org/wiki/Frobenius_matrix). To buy a smaller number of factors in the decomposition, we pay with sparsity.


Though the $\mathrm E\_{i j}$ and $\mathrm E\_{i j}^{-1}$ matrices have $n+1$ nonzero entries, only one of them, the $(i,j)$-th entry, is in general different from $1$. Therefore, we only need one real number, $\beta\_{ij}$, and one pair of positive integers, $(i,j)$, to fully specify the elementary matrix $\mathrm E\_{i j}$ (or its inverse). Taking into account also the $n$ nonzero entries on the main diagonal of the diagonal matrix $\mathrm D$, much of the data needed to build the factors in the decomposition can be stored in the following $n \times n$ matrix


$$\begin{bmatrix} d\_1 & \beta\_{12} & \dots & \beta\_{1n}\\ \beta\_{21} & d\_2 & \dots & \beta\_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ \beta\_{n1} & \beta\_{n2} & \dots & d\_n\end{bmatrix}$$


Each permutation matrix can be specified using $n$ positive integers. Thus, at most $n^2$ reals and $n^2$ integers are needed to specify the decomposition.




---


Example
=======


Consider the (invertible) $3 \times 3$ magic square


$$\mathrm A := \begin{bmatrix} 2 & 9 & 4\\7 & 5 & 3\\6 & 1 & 8\end{bmatrix}$$


Using [SymPy](http://www.sympy.org),



```
>>> from sympy import *
>>> A = Matrix([[2,9,4],[7,5,3],[6,1,8]])

```

Eliminate two entries below the $(1,1)$-th entry:



```
>>> E21 = eye(3) + Matrix([[0,0,0],[-7.0/2.0,0,0],[0,0,0]])
>>> E31 = eye(3) + Matrix([[0,0,0],[0,0,0],[-3.0,0,0]]) 
>>> E31 * E21 * A
Matrix([
[2,     9,     4],
[0, -26.5, -11.0],
[0, -26.0,  -4.0]])

```

Eliminate the entry below the $(2,2)$-th entry:



```
>>> E32 = eye(3) + Matrix([[0,0,0],[0,0,0],[0,-26.0/26.5,0]])
>>> E32 * E31 * E21 * A
Matrix([
[2,                    9,                4],
[0,                -26.5,            -11.0],
[0, 1.33226762955019e-15, 6.79245283018868]])

```

Numerical noise is now appearing. Eliminate two entries above the $(3,3)$-th entry



```
>>> E23 = eye(3) + Matrix([[0,0,0],[0,0,11.0/6.79245],[0,0,0]])
>>> E13 = eye(3) + Matrix([[0,0,-4.0/6.79245],[0,0,0],[0,0,0]])
>>> E13 * E23 * E32 * E31 * E21 * A
Matrix([
[2.0,                  9.0, -1.66666736145515e-6],
[  0,                -26.5,  4.58333524377963e-6],
[  0, 1.33226762955019e-15,     6.79245283018868]])

```

Eliminate the entry above the $(2,2)$-th entry:



```
>>> E12 = eye(3) + Matrix([[0,9.0/26.5,0],[0,0,0],[0,0,0]])
>>> E12 * E13 * E23 * E32 * E31 * E21 * A
Matrix([
[2.0, 1.99840144432528e-15, -1.10062938096789e-7],
[  0,                -26.5,  4.58333524377963e-6],
[  0, 1.33226762955019e-15,     6.79245283018868]])

```

We extract the diagonal:



```
>>> D = diag(2.0,-26.5,6.79245)
>>> D
Matrix([
[2.0,     0,       0],
[  0, -26.5,       0],
[  0,     0, 6.79245]])

```

Finally, we reconstruct the magic square:



```
>>> (E21**-1) * (E31**-1) * (E32**-1) * (E23**-1) * (E13**-1) * (E12**-1) * D
Matrix([
[2.0, 9.0,              4.0],
[7.0, 5.0,              3.0],
[6.0, 1.0, 7.99999716981132]])

```

The reconstruction is not exact due to numerical noise, but it is close enough. Thus, magic square $\mathrm A$ can be approximately factored as follows


$$\underbrace{\begin{bmatrix} 1 & 0 & 0\\ 3.5 & 1 & 0\\3 & 0 & 1\end{bmatrix}}\_{= \mathrm E\_{21}^{-1} \mathrm E\_{31}^{-1}} \, \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\0 & 0.98 & 1\end{bmatrix} \, \underbrace{\begin{bmatrix} 1 & 0 & 0.59\\0 & 1 & -1.62\\0 & 0 & 1\end{bmatrix}}\_{\approx \mathrm E\_{23}^{-1} \mathrm E\_{13}^{-1}} \, \begin{bmatrix} 1 & -0.34 & 0\\0 & 1 & 0\\0 & 0 & 1\end{bmatrix} \, \underbrace{\begin{bmatrix} 2 & 0 & 0\\ 0 & -26.5 & 0\\0 & 0 & 6.79\end{bmatrix}}\_{\approx \mathrm D}$$


where the 1st and 3rd of the five matrix factors are [Frobenius matrices](https://en.wikipedia.org/wiki/Frobenius_matrix). In this particular example, no permutation matrices (other than $\mathrm I\_n$, that is) were needed.


