Decomposing a matrix into a product of sparse matrices

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.

• Must all the matrices in the product have the same sparsity pattern? Dec 29 '16 at 18:31
• no...the patterns (and density) could be different...as long as they remain sparse Dec 29 '16 at 19:13

Given an invertible $n \times n$ matrix $\mathrm A$, we perform Gaussian elimination until we obtain a (nonsingular) diagonal matrix. In other words, we left-multiply $\mathrm A$ by permutation matrices (whose inverses are their transposes) and by elementary matrices 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. 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,

>>> 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. In this particular example, no permutation matrices (other than $\mathrm I_n$, that is) were needed.

• @Rodrigo...this is very nice. A minor point : wouldn't there be some permutation matrices in the product too? or are you combining these with the $E_{ij}$ factors? Dec 30 '16 at 17:26
• in the special case where the field is GF(2) and there are no numerical issues, it still seems that we still need permutation matrices in case the pivot is not in the right place...this would add $n-1$ matrices to the $n(n-1)$ matrices above; I think the permutation can be absorbed into one of the $E_{ij}$'s or as a slight modification to the Frobenius matrices but I haven't worked that out yet Dec 30 '16 at 18:21
• It is not true that pivots are needed only when $A$ is ill-conditioned. There are matrices, such as $\begin{bmatrix}0&1\\1&0\end{bmatrix}$, which have condition number $1$ but for which pivoting is necessary. Dec 30 '16 at 20:54
• Yes - Gaussian elimination with pivoting is equivalent to Gaussian elimination on a permuted $PA$, but as far as I know there is no way to determine $P$ in advance, you just have to permute while you're doing the elimination. Anyway, for this problem the simplest solution is what @unknown suggests, absorbing the permutations into the $E_{ij}$. Dec 30 '16 at 21:14
• All seems correct to me. It's a very clear explanation, great job. Very minor remark: most of your "numerical noise" comes from truncating 6.79245283018868 to 6.79245 when you use it as input. Dec 30 '16 at 21:49

The question of deciding if a specific matrix is a product of "a few" "sparse" matrices seems very interesting to me and perhaps rather hard.

Likewise the question of how bad things could be if the underlying field is finite (say $\mathbb{Z}_2.$)

However a random matrix is, in some sense, highly unlikely to have any unusual decompositions.

Here is a rather extreme conjecture which seems plausible to me: define the density of an $n \times n$ matrix to be the number entries which are not $0$ or $1$ and the density of a product of such matrices to be the sum of the densities of the factors. Then with probability $1$ a random real matrix has density $n^2$ and, more generally, no factorization has lower density. However every matrix is a product of $n^2$ matrices of density $1$ and , with probability $1$, there are $n^2!$ such decompositions.

Perhaps a better definition , at least for integer matrices (with huge entries), would be the number of entries in which a matrix differs from the identity matrix or the previous definition plus $1.$

• The map from the nontrivial coefficients in the factors to the product matrix is a polynomial map; so if your conjecture were false then we'd have a polynomial map from $\mathbb{R}^{m}$ (with $m<n^2$) to $\mathbb{R}^{n^2}$ whose image has nonzero Lebesgue measure. My algebraic geometry is very rusty, but this seems unlikely. Dec 29 '16 at 16:28
• as a special case, matrices that are known to be the inverse of a sparse matrix are of interest...also the field could be GF(2). I don't know if these additional restrictions change the problem much Dec 29 '16 at 17:58