Matrix decomposition into square positive integer matrices This is an attempt at an analogy with prime numbers. Let's consider only square matrices with positive integer entries. Which of them are 'prime' and how to decompose such a matrix in general?
To illustrate, there is a product of two general $2 \times 2$ matrices:
$$AB=\left[ \begin{matrix} a_{11} &  a_{12} \\  a_{21} &  a_{22} \end{matrix} \right] \left[ \begin{matrix} b_{11} &  b_{12} \\  b_{21} &  b_{22} \end{matrix} \right]=\left[ \begin{matrix} a_{11} b_{11}+a_{12} b_{21} &  a_{11} b_{12}+a_{12} b_{22} \\  a_{21} b_{11}+a_{22} b_{21} &  a_{21} b_{12}+a_{22} b_{22} \end{matrix} \right]$$
Exchanging $a$ and $b$ we obtain the expression for the other product $BA$.
Now, if we allow zero, negative and/or rational entries we can probably decompose any matrix in an infinite number of ways.
However, if we restrict ourselves:
$$a_{jk},~b_{jk} \in \mathbb{N}$$
The problem becomes well defined.

Is there an algorithm to decompose an arbitrary square positive integer matrix into a product of several positive integer matrices of the same dimensions?

There is a set of matrices which can'be be decomposed, just like the prime numbers (or irreducible polynomials, for example). The most trivial one is (remember, zero entries are not allowed):
$$\left[ \begin{matrix} 1 & 1 \\  1 &  1 \end{matrix} \right]$$
There are no natural numbers $a_{11},b_{11},a_{12},b_{21}$, such that:
$$a_{11} b_{11}+a_{12} b_{21}=1$$
The same extends to any dimension $d$. Any 'composite' $d \times d$ matrix will have all entries $ \geq d$. Thus, for square matrices we can name several more 'primes':
$$\left[ \begin{matrix} 2 & 1 \\  1 &  1 \end{matrix} \right],~~~\left[ \begin{matrix} 1 & 2 \\  1 &  1 \end{matrix} \right],~~~\left[ \begin{matrix} 1 & 1 \\  2 &  1 \end{matrix} \right],~~~\left[ \begin{matrix} 1 & 1 \\  1 &  2 \end{matrix} \right],~~~\left[ \begin{matrix} 2 & 2 \\  1 &  1 \end{matrix} \right],~~~\left[ \begin{matrix} 1 & 1 \\  2 &  2 \end{matrix} \right], \dots$$
And in general, any matrix which has at least one entry equal to $1$.
It makes sense, that most entries in 'composite' matrices will be large, since we are multiplying and adding natural numbers. For example:
$$\left[ \begin{matrix} 1 &  2 &  4 \\  3 &  3 &  1 \\  3 &  4 &  4 \end{matrix} \right] \left[ \begin{matrix} 2 &  5 &  5 \\  4 &  5 &  5 \\  5 &  1 &  4 \end{matrix} \right]=\left[ \begin{matrix} 30 &  19 &  31 \\  23 &  31 &  34 \\  42 &  39 &  51 \end{matrix} \right]$$
$$\left[ \begin{matrix} 2 &  5 &  5 \\  4 &  5 &  5 \\  5 &  1 &  4 \end{matrix} \right] \left[ \begin{matrix} 1 &  2 &  4 \\  3 &  3 &  1 \\  3 &  4 &  4 \end{matrix} \right] =\left[ \begin{matrix} 32 &  39 &  33 \\  34 &  43 &  41 \\  20 &  29 &  37 \end{matrix} \right]$$

If no decomposition algorithm for this case exists, is it at least possible to recognize a matrix that can't be decomposed according to the above rules?


I now understand that an algorithm for the general case is unlikely. Searching the web, I only found several articles, concerned with positive integer matrix decomposition into binary matrices.

I tried to consider a particular case, which seems to be the most simple:
$$C=\left[ \begin{matrix} 2 &  c_b \\  c_a &  c \end{matrix} \right]$$
Where $c_a,c_b,c \in \mathbb{N}$ - arbitrary natural numbers. This matrix obviously is either 'prime' or can be decomposed into two 'prime' matrices.
$$\left[ \begin{matrix} 2 &  c_b \\  c_a &  c \end{matrix} \right]=\left[ \begin{matrix} 1 &  1 \\  a_{1} &  a_{2} \end{matrix} \right] \left[ \begin{matrix} 1 &  b_{1} \\  1 &  b_{2} \end{matrix} \right]$$
$$ \begin{cases} a_1+a_2=c_a \\ b_1+b_2=c_b \\ a_1b_1+a_2b_2=c \end{cases} $$
Again  $a_1,b_1,a_2,b_2 \in \mathbb{N}$. From which we immediately obtain the most important condition for $C$ to be 'composite':

$$ \max(c_a,c_b)\leq c<c_ac_b$$

The equality on the left is possible only in the trivial cases of $a_1=a_2=1$ or $b_1=b_2=1$, which I'm not going to consider.
If this conditions doesn't hold, then $C$ is 'prime'. Which gives us new examples of 'prime' matrices:
$$\left[ \begin{matrix} 2 &  5 \\  7 &  39 \end{matrix} \right],~~~\left[ \begin{matrix} 2 &  15 \\  8 &  11 \end{matrix} \right],~~~\dots$$

Another property is this - if we simultaneously permute $a_k$ and $b_k$, then $C$ doesn't change. Which means, that if $C$ is composite, it has at least two factorizations:
$$\left[ \begin{matrix} 2 &  c_b \\  c_a &  c \end{matrix} \right]=\left[ \begin{matrix} 1 &  1 \\  a_{1} &  a_{2} \end{matrix} \right] \left[ \begin{matrix} 1 &  b_{1} \\  1 &  b_{2} \end{matrix} \right]=\left[ \begin{matrix} 1 &  1 \\  a_{2} &  a_{1} \end{matrix} \right] \left[ \begin{matrix} 1 &  b_{2} \\  1 &  b_{1} \end{matrix} \right]$$
Thus, without loss of generality, we can impose:
$$a_1 >a_2,~~~~\frac{c_a}{2}<a_1<c_a \tag{*}$$
The trivial case $a_1=a_2$ gives nilpotent matrix $C$, which I'm not going to consider here.
Solving the system of equations above for $b_1$, we obtain:
$$b_1=\frac{c-c_b(c_a-a_1)}{2a_1-c_a} \tag{**}$$
$$c_a-\frac{c}{c_b}<a_1<c_a$$
The above means that both the intervals $(\frac{c_a}{2},c_a)$ and $(c_a-\frac{c}{c_b},c_a)$ should contain at least one integer:

$$\frac{c_a}{2}>1,~~~\frac{c}{c_b}\geq 2$$

With the condition $(*)$ we have two distinct cases:
$$b_1>b_2,~~~\frac{c_b}{2}<b_1<c_b,~~~~c_b-\frac{c}{c_a}<b_1  \tag{1}$$

$$c>\frac{c_ac_b}{2},~~~~~\frac{c}{c_a} \geq 2,~~~~~\frac{c_b}{2}>1$$

$$b_1<b_2,~~~b_1<\frac{c_b}{2},~~~~b_1<c_b-\frac{c}{c_a} \tag{2}$$

$$c<\frac{c_ac_b}{2},~~~~~\frac{c_b}{2}>1$$

The conditions on $c$ follow from the rearrangement inequality.
There seems to be no additional condition for the case $(2)$, but I'm not sure.
As an example of a composite matrix:
$$C=\left[ \begin{matrix} 2 &  5 \\  7 &  19 \end{matrix} \right]$$
This is case $(1)$ and all the necessary conditions hold.
$$\frac{c_a}{2}=3.5,~~~~c_a-\frac{c}{c_b}=3.2 \rightarrow a_1 = 4,5,6$$
$$\frac{c_b}{2}=2.5,~~~~c_b-\frac{c}{c_a}=2\frac{4}{7} \rightarrow b_1 = 3,4$$
Using $(**)$ we obtain:
$$a_1=4 \to b_1=4$$
$$a_1=5 \to b_1=3$$
Which gives us (I think) all the solutions:
$$\left[ \begin{matrix} 2 &  5 \\  7 &  19 \end{matrix} \right]=\left[ \begin{matrix} 1 &  1 \\  4 &  3 \end{matrix} \right] \left[ \begin{matrix} 1 &  4 \\  1 &  1 \end{matrix} \right]=\left[ \begin{matrix} 1 &  1 \\  3 &  4 \end{matrix} \right] \left[ \begin{matrix} 1 &  1 \\  1 &  4 \end{matrix} \right]$$
$$\left[ \begin{matrix} 2 &  5 \\  7 &  19 \end{matrix} \right]=\left[ \begin{matrix} 1 &  1 \\  5 &  2 \end{matrix} \right] \left[ \begin{matrix} 1 &  3 \\  2 &  1 \end{matrix} \right]=\left[ \begin{matrix} 1 &  1 \\  2 &  5 \end{matrix} \right] \left[ \begin{matrix} 1 &  2 \\  1 &  3 \end{matrix} \right]$$


This case was considered to illustrate how I would approach this problem. Which is very long and complicated way. If someone knows better and faster ways, I would be grateful. 

