Finding Permutation Matrix with NumPy-CodePudding

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?

CodePudding user response：

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]])

CodePudding user response：

solve computes x in the equation Ax = B. You can use it, but the arguments have to match.

Your equation is XA = B. You can take the transpose of both sides to get A^TP^T = B^T. So you can get the result from

np.linalg.solve(A.T, B.T).T

To get a boolean result that's guaranteed to be zeros and ones, you can apply a threshold:

np.linalg.solve(A.T, B.T).T > 0.5