Decomposing a matrix into a product of sparse matrices

Question

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.

Answer 1

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.

Answer 2

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.$