Finding Permutation Matrix with NumPy

December 9, 2022

I am looking for the correct permutation matrix that would take matrix a and turn it into matrix b given

a = np.array([[1,4,7,-2],[3,0,-2,-1],[-4,2,1,0],[-8,-3,-1,2]])
b = np.array([[-4,2,1,0],[3,0,-2,-1],[-8,-3,-1,2],[1,4,7,-2]])

I tried

x = np.linalg.solve(a,b)

However, I know this is incorrect and it should be

np.array([[0,0,1,0],[0,1,0,0],[0,0,0,1],[1,0,0,0]])

What numpy code would deliver this matrix from the other two?

>Solution :

Generally, if you have some PA = B and you want P then you need to solve the equation for P. Matrix multiplication is not commutative, so you have to right multiply both sides by the inverse of A. With numpy, the function to get the inverse of a matrix is np.linalg.inv().

Using the matrix multiplication operator @, you can right-multiply with b to get P, taking note that you will end up with floating point precision errors:

In [4]: b @ np.linalg.inv(a)
Out[4]:
array([[-1.38777878e-16, -3.05311332e-16,  1.00000000e+00,
        -1.31838984e-16],
       [ 0.00000000e+00,  1.00000000e+00,  6.93889390e-17,
         0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00, -1.11022302e-16,
         1.00000000e+00],
       [ 1.00000000e+00, -4.44089210e-16,  5.55111512e-17,
         0.00000000e+00]])

Rounding to a single decimal place reveals that this is indeed the correct permutation matrix:

In [5]: np.round(b @ np.linalg.inv(a))
Out[5]:
array([[-0., -0.,  1., -0.],
       [ 0.,  1.,  0.,  0.],
       [ 0.,  0., -0.,  1.],
       [ 1., -0.,  0.,  0.]])

If you’re like me, negative zero makes you uneasy, so cast to int if you want:

In [6]: np.round(b @ np.linalg.inv(a)).astype(int)
Out[6]:
array([[0, 0, 1, 0],
       [0, 1, 0, 0],
       [0, 0, 0, 1],
       [1, 0, 0, 0]])

As @Mad Physicist points out, you can also use > 0.5 to convert the matrix to a boolean matrix:

In [7]: bool_P = (b @ np.linalg.inv(a)) > 0.5

In [8]: bool_P @ a
Out[8]:
array([[-4,  2,  1,  0],
       [ 3,  0, -2, -1],
       [-8, -3, -1,  2],
       [ 1,  4,  7, -2]])

In [9]: bool_P @ a == b  # our original equation, PA = B
Out[9]:
array([[ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True]])