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.

Advertisements