# Python and NumPy: Array Division

## Array Division of Square Matrices Using the Inverse

Lets create the following two matrices and multiply them using the dot product to get C:

Python
[[19 22]
[43 50]]

Matrix Division is always gently brushed over at school with only the special case of a 2 by 2 square matrix being solved using the inverse of the square matrix.

For:

A⁕B=C

With AInv, the Inverse Matrix of A is defined as:

$\displaystyle \text{AInv}=\frac{1}{{ad-bc}}\left( {\begin{array}{*{20}{c}} d & {-b} \\ {-c} & a \end{array}} \right)$

Substituting in the numbers:

$\displaystyle \text{AInv=}\frac{1}{{\left( {\left( 1 \right)*\left( 4 \right)-\left( 2 \right)*\left( 3 \right)} \right)}}*\left( {\begin{array}{*{20}{c}} 4 & {-2} \\ {-3} & 1 \end{array}} \right)=-\frac{1}{2}\left( {\begin{array}{*{20}{c}} 4 & {-2} \\ {-3} & 1 \end{array}} \right)=\left( {\begin{array}{*{20}{c}} {-2} & 1 \\ {\frac{3}{2}} & {-\frac{1}{2}} \end{array}} \right)$

In Python to find the Inv Matrix we can use the function inv in the linear algebra (linalg) section of the numpy library:

Python
[[-2.   1. ]
[ 1.5 -0.5]]

Now when a Matrix is multiplied by it’s inverse will give the Identity matrix I:

Ainv⁕A=A⁕Ainv=I

Lets check this in Python:

Python
[[1.0000000e+00 4.4408921e-16]
[0.0000000e+00 1.0000000e+00]]
Python
[[1.00000000e+00 1.11022302e-16]
[0.00000000e+00 1.00000000e+00]]

The numbers in the anti-diagonal are effectively zero, showing only to the 16th decimal place. We could use np.round to round to 5 decimal places:

Python
[[1. 0.]
[0. 1.]]

If we multiply the expression

Python
[[1. 0.]
[0. 1.]]

A⁕B=C

Throughout by the inverse of A:

AInv⁕A⁕B=AInv⁕C

I⁕B=AInv⁕C

When a Matrix is multiplied by the Identity Matrix it is unchanged

We can check this in Python using:

Python
[[5. 6.]
[7. 8.]]
Python
[[5. 6.]
[7. 8.]]

Therefore:

AInv⁕A⁕B=AInv⁕C

Becomes:

B=AInv⁕C

This gives the correct value of B. Lets check this in Python:

Python
[[5. 6.]
[7. 8.]]

This can be thought of as:

A(1,1)⁕B(1,1)+A(1,2)⁕B(2,1)=C(1,1) A(1,1)⁕B(1,2)+A(1,2)⁕B(2,2)=C(1,2) A(2,1)⁕B(1,1)+A(2,2)⁕B(2,1)=C(2,1) A(2,1)⁕B(1,2)+A(2,2)⁕B(2,2)=C(2,2)

Assuming we don’t know the values of B, but we know the values of A and C. This is 4 sets of equations, containing 4 unknowns:

1. B(1,1)+2⁕B(2,1)=19
2. B(1,2)+2⁕B(2,2)=22
3. 3⁕B(1,1)+4⁕B(2,1)=43
4. 3⁕B(1,2)+4⁕B(2,2)=50

Rearranging 1:

B(1,1)+2⁕B(2,1)=19 B(1,1)=19-2⁕B(2,1)

Substituting into 3:

3⁕(19-2⁕B(2,1))+4⁕B(2,1)=43 57-6⁕B(2,1)+4⁕B(2,1)=43 -2⁕B(2,1)=43-57 -2⁕B(2,1)=-14 B(2,1)=-14/-2 B(2,1)=7

Substituting this into 1 (which was rearranged):

B(1,1)=19-2⁕B(2,1)
B(1,1)=19-2⁕7
B(1,1)=19-14
B(1,1)=5

Now rearranging 2:

B(1,2)+2⁕B(2,2)=22
B(1,2)=22-2⁕B(2,2)

Now substituting into 4:

3⁕B(1,2)+4⁕B(2,2)=50
3⁕(22-2⁕B(2,2))+4⁕B(2,2)=50
66-6⁕B(2,2)+4⁕B(2,2)=50
-6⁕B(2,2)+4⁕B(2,2)=50-66
-2⁕B(2,2)=-16
B(2,2)=-16/-2
B(2,2)=8

Substituting this into the rearranged 2:

B(1,2)=22-2⁕B(2,2)
B(1,2)=22-2⁕8
B(1,2)=22-16
B(1,2)=6

This gives:

B(1,1)=5
B(1,2)=6
B(2,1)=7
B(2,2)=8

As expected.

# Square Matrix Division

Python can perform Matrix Division using the least squares function solve which carries out matrix division using the inverse and hence requires a square matrice or lstsq of the linear algebra library which uses least squares to find the coefficients. Note I have found that NumPy’s lstsq method to be substantially inferior to MATLAB/Octave.

Let’s look at the case:

A⁕B=C

Solving for B we get

Python
[[5. 6.]
[7. 8.]]

[[5. 6.]
[7. 8.]]

[]

2

[5.4649857  0.36596619]

Starting with:

$\displaystyle A*B=C$

$latex \displaystyle \frac{1}{A}AB=\frac{1}{A}*C$

$\displaystyle B=\frac{1}{A}*C$

To find B opposed to A it is slightly different:

$\displaystyle A*B=C$

$latex \displaystyle AB\frac{1}{B}=C*\frac{1}{B}$

$\displaystyle A=C*\frac{1}{B}$

This gives us the same form we used previously but with both sides as Transposes.

Python
[[1. 3.]
[2. 4.]]

[[1. 3.]
[2. 4.]]

[]

2

[5.4649857  0.36596619]

We can then Transpose the result:

Python
[[1. 2.]
[3. 4.]]

## Non-Square Matrix Division

We can use the data we multiplied together to work back from to illustrate division:

Assuming we know S and A and want to calculate E:

Python
array([[ 3,  5,  4],
[15, 25, 20],
[ 3,  5,  4],
[ 6, 10,  8]])

[[0.06 0.1  0.08]]

[[0.03225806]
[0.16129032]
[0.03225806]
[0.06451613]]

Here the results are not precisely as expected. S1 is S*1/50 and E1 is approximately E*1/31.