4

I want to implement my own LU decomposition P,L,U = my_lu(A), so that given a matrix A, computes the LU decomposition with partial pivoting. But I only know how to do it without pivoting. Can anyone help to do the partial pivoting?

def lu(A):

    import numpy as np

    # Return an error if matrix is not square
    if not A.shape[0]==A.shape[1]:
        raise ValueError("Input matrix must be square")

    n = A.shape[0] 

    L = np.zeros((n,n),dtype='float64') 
    U = np.zeros((n,n),dtype='float64') 
    U[:] = A 
    np.fill_diagonal(L,1) # fill the diagonal of L with 1

    for i in range(n-1):
        for j in range(i+1,n):
            L[j,i] = U[j,i]/U[i,i]
            U[j,i:] = U[j,i:]-L[j,i]*U[i,i:]
            U[j,i] = 0
    return (L,U)
Jacqueline Yishu Liu
  • 41
  • 1
  • 1
  • 2
  • 4
    Isn't this more adequately asked on http://math.stackexchange.com/ ? – h7r Feb 10 '15 at 21:06
  • 1
    [LU factorization](http://www4.ncsu.edu/~kksivara/ma505/handouts/lu.pdf) and [LU factorization with pivoting](http://www4.ncsu.edu/~kksivara/ma505/handouts/lu-pivot.pdf) from Trefethen and Bau, with crystal clear pseudocode – gboffi Feb 10 '15 at 23:19

2 Answers2

5
def naive_lu_factor(A):
    """
        No pivoting.

        Overwrite A with:
            U (upper triangular) and (unit Lower triangular) L 
        Returns LU (Even though A is also overwritten)
    """
    n = A.shape[0]
    for k in range(n-1):                
        for i in range(k+1,n):          
            A[i,k] = A[i,k]/A[k,k]      # " L[i,k] = A[i,k]/A[k,k] "
            for j in range(k+1,n):      
                A[i,j] -= A[i,k]*A[k,j] # " U[i,j] -= L[i,k]*A[k,j] "

    return A # (if you want)

def lu_factor(A):
    """
        LU factorization with partial pivorting

        Overwrite A with: 
            U (upper triangular) and (unit Lower triangular) L 
        Return [LU,piv] 
            Where piv is 1d numpy array with row swap indices 
    """
    n = A.shape[0]
    piv = np.arange(0,n)
    for k in range(n-1):

        # piv
        max_row_index = np.argmax(abs(A[k:n,k])) + k
        piv[[k,max_row_index]] = piv[[max_row_index,k]]
        A[[k,max_row_index]] = A[[max_row_index,k]]

        # LU 
        for i in range(k+1,n):          
            A[i,k] = A[i,k]/A[k,k]      
            for j in range(k+1,n):      
                A[i,j] -= A[i,k]*A[k,j] 

    return [A,piv]

def ufsub(L,b):
    """ Unit row oriented forward substitution """
    for i in range(L.shape[0]): 
        for j in range(i):
            b[i] -= L[i,j]*b[j]
    return b

def bsub(U,y):
    """ Row oriented backward substitution """
    for i in range(U.shape[0]-1,-1,-1): 
        for j in range(i+1, U.shape[1]):
            y[i] -= U[i,j]*y[j]
        y[i] = y[i]/U[i,i]
    return y

Solve for x (with and without partial pivoting) using unit forward and backward substitution:

# No partial pivoting
LU = naive_lu_factor(A)
y = ufsub( LU, b )
x = bsub(  LU, y )

# Partial pivoting
LU, piv = lu_factor(A)                      
b = b[piv]
y = ufsub( LU, b )
x = bsub(  LU, y )
LJ Brown
  • 59
  • 1
  • 2
  • The pivoting has a mistake. I'm not sure what it is, but I'm making the same one. The permutations play a role in the Lk. – RicardoMM May 19 '22 at 23:52
-2

You can use Scipy's scipy.linalg.lu for this.

http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.linalg.lu.html

Please check this link also:

http://www.quantstart.com/articles/LU-Decomposition-in-Python-and-NumPy

Gunjan naik
  • 501
  • 1
  • 7
  • 18