Gaussian Elimination with Partial Pivoting

In [1]: 
import numpy as np
In [38]: 
n = 3

np.random.seed(22)
A = np.round(5*np.random.randn(n, n))
A[0,0] = 0
A
Out[38]:
array([[ 0., -7.,  5.],
       [-1., -2., -5.],
       [ 5., -6.,  3.]])

Now define a function row_swap_mat(i, j) that returns a permutation matrix that swaps row i and j:

In [39]: 
def row_swap_mat(i, j):
    P = np.eye(n)
    P[i] = 0
    P[j] = 0
    P[i, j] = 1
    P[j, i] = 1
    return P

What do these matrices look like?

In [40]: 
row_swap_mat(0,1)
Out[40]:
array([[ 0.,  1.,  0.],
       [ 1.,  0.,  0.],
       [ 0.,  0.,  1.]])

Do they work?

In [41]: 
row_swap_mat(0,1).dot(A)
Out[41]:
array([[-1., -2., -5.],
       [ 0., -7.,  5.],
       [ 5., -6.,  3.]])

U is the copy of A that we'll modify:

In [63]: 
U = A.copy()

Create P1 to swap up the right row:

In [64]: 
P1 = row_swap_mat(0, 2)
U = P1.dot(U)
U
Out[64]:
array([[ 5., -6.,  3.],
       [-1., -2., -5.],
       [ 0., -7.,  5.]])
In [65]: 
M1 = np.eye(n)
M1[1,0] = -U[1,0]/U[0,0]
M1
Out[65]:
array([[ 1. ,  0. ,  0. ],
       [ 0.2,  1. ,  0. ],
       [ 0. ,  0. ,  1. ]])
In [66]: 
U = M1.dot(U)
U
Out[66]:
array([[ 5. , -6. ,  3. ],
       [ 0. , -3.2, -4.4],
       [ 0. , -7. ,  5. ]])

Create P2 to swap up the right row:

In [67]: 
P2 = row_swap_mat(1,2)
U = P2.dot(U)
U
Out[67]:
array([[ 5. , -6. ,  3. ],
       [ 0. , -7. ,  5. ],
       [ 0. , -3.2, -4.4]])
In [68]: 
M2 = np.eye(n)
M2[2,1] = -U[2,1]/U[1,1]
M2
Out[68]:
array([[ 1.        ,  0.        ,  0.        ],
       [ 0.        ,  1.        ,  0.        ],
       [ 0.        , -0.45714286,  1.        ]])
In [69]: 
U = M2.dot(U)
U
Out[69]:
array([[ 5.        , -6.        ,  3.        ],
       [ 0.        , -7.        ,  5.        ],
       [ 0.        ,  0.        , -6.68571429]])

So we've built $M_2P_2M_1P_1A=U$.

In [75]: 
M2.dot(P2).dot(M1).dot(P1).dot(A)
Out[75]:
array([[ 5.        , -6.        ,  3.        ],
       [ 0.        , -7.        ,  5.        ],
       [ 0.        ,  0.        , -6.68571429]])

Can't simply sort the permutation matrices out of the way:

In [77]: 
M2.dot(P2).dot(M1).dot(P1)
Out[77]:
array([[ 0.        ,  0.        ,  1.        ],
       [ 1.        ,  0.        ,  0.        ],
       [-0.45714286,  1.        ,  0.2       ]])
In [80]: 
M2.dot(M1).dot(P2).dot(P1)
Out[80]:
array([[ 0.        ,  0.        ,  1.        ],
       [ 1.        ,  0.        ,  0.2       ],
       [-0.45714286,  1.        , -0.09142857]])

But: We now know the right reordering!

In [87]: 
P = P2.dot(P1)
P
Out[87]:
array([[ 0.,  0.,  1.],
       [ 1.,  0.,  0.],
       [ 0.,  1.,  0.]])
In [90]: 
PA = P.dot(A)
PA
Out[90]:
array([[ 5., -6.,  3.],
       [ 0., -7.,  5.],
       [-1., -2., -5.]])