Untitled

import numpy as np

def pivotedOuterProductLU(A):
    """
    Partially pivoted LU factorization (row pivoting) using an
    outer-product (rank-1) update approach.

    Given an n x n matrix A, returns (L, U, P) such that:
        P * A = L * U,
    where L is lower-triangular (1 on the diagonal), U is upper-triangular,
    and P is the permutation matrix that encodes the row swaps.
    """
    A = A.astype(float)
    n = A.shape[0]

    # Initialize U as a copy of A; L as identity; P as identity
    U = A.copy()
    L = np.eye(n, dtype=float)
    P = np.eye(n, dtype=float)

    # Loop over each "pivot" column k
    for k in range(n - 1):

        # 1) Find the pivot row in column k, among rows k..(n-1)
        pivot_row = np.argmax(np.abs(U[k:, k])) + k

        # 2) Swap row k with pivot_row in U, P, and the multipliers of L
        if pivot_row != k:
            U[[k, pivot_row], :] = U[[pivot_row, k], :]
            P[[k, pivot_row], :] = P[[pivot_row, k], :]
            # Swap only the first k columns in L (the multipliers so far)
            L[[k, pivot_row], :k] = L[[pivot_row, k], :k]

        # 3) The pivot is now U[k, k]; check for zero pivot
        pivot_val = U[k, k]
        if abs(pivot_val) < 1e-15:
            raise ValueError("Encountered near-zero pivot; factorization fails.")

        # 4) Define the multipliers below the pivot (outer-product step)
        for i in range(k+1, n):
            L[i, k] = U[i, k] / pivot_val
            # Subtract the rank-1 update from row i of U
            U[i, k+1:] = U[i, k+1:] - L[i, k] * U[k, k+1:]

    return L, U, P

# --------------- TEST on the given 4x4 matrix ---------------
if __name__ == "__main__":
    A_test = np.array([
        [ 4,  3,  2,  1],
        [ 8,  8,  5,  2],
        [16, 12, 10,  5],
        [32, 24, 20, 11]
    ], dtype=float)

    L, U, P = pivotedOuterProductLU(A_test)

    print("A =\n", A_test, "\n")
    print("L =\n", L, "\n")
    print("U =\n", U, "\n")
    print("P =\n", P, "\n")

    # Check that P * A == L * U (up to floating-point rounding)
    lhs = P @ A_test
    rhs = L @ U
    print("P * A =\n", lhs, "\n")
    print("L * U =\n", rhs, "\n")
    print("Difference norm =", np.linalg.norm(lhs - rhs))