Refactor lup_decompose() into plu_decompose()

README.md

@@ -268,7 +268,7 @@ True
 
 Galois field arrays also contain matrix decomposition routines not included in NumPy. These include:
 
-- `FieldArray.row_reduce`, `FieldArray.lu_decompose`, `FieldArray.lup_decompose`
+- `FieldArray.row_reduce`, `FieldArray.lu_decompose`, `FieldArray.plu_decompose`
 
 #### NumPy ufunc methods
 

docs/basic-usage.rst

@@ -159,7 +159,7 @@ operation in :math:`\mathrm{GF}(p^m)` not in :math:`\mathbb{R}`. Some of these f
 
 Galois field arrays also contain matrix decomposition routines not included in NumPy. These include:
 
-- :func:`galois.FieldArray.row_reduce`, :func:`galois.FieldArray.lu_decompose`, :func:`galois.FieldArray.lup_decompose`
+- :func:`galois.FieldArray.row_reduce`, :func:`galois.FieldArray.lu_decompose`, :func:`galois.FieldArray.plu_decompose`
 
 NumPy ufunc methods
 ...................

galois/_fields/_linalg.py

@@ -188,15 +188,15 @@ def lu_decompose(A):
     return L, U
 
 
-def lup_decompose(A):
+def plu_decompose(A):
     if not (A.ndim == 2 and A.shape[0] == A.shape[1]):
         raise ValueError(f"Argument `A` must be a square matrix, not {A.shape}.")
 
     field = type(A)
     n = A.shape[0]
     Ai = A.copy()
     L = field.Zeros((n,n))
-    P = field.Identity(n)
+    P = field.Identity(n)  # Row permutation matrix
 
     for i in range(0, n-1):
         if Ai[i,i] == 0:
@@ -219,7 +219,8 @@ def lup_decompose(A):
     L[-1,-1] = 1  # Set the final diagonal to 1
     U = Ai
 
-    return L, U, P
+    # NOTE: Return column permutation matrix
+    return P.T, L, U
 
 
 ###############################################################################
@@ -271,7 +272,8 @@ def det(A):
     elif n == 3:
         return A[0,0]*(A[1,1]*A[2,2] - A[1,2]*A[2,1]) - A[0,1]*(A[1,0]*A[2,2] - A[1,2]*A[2,0]) + A[0,2]*(A[1,0]*A[2,1] - A[1,1]*A[2,0])
     else:
-        L, U, P = lup_decompose(A)
+        P, L, U = plu_decompose(A)
+        P = P.T  # Convert row permutation matrix into column permutation matrix
         idxs = np.arange(0, n)
         nrows = n - np.count_nonzero(P[idxs,idxs]) # The number of moved rows
         S = max(nrows - 1, 0)  # The number of row swaps

galois/_fields/_main.py

@@ -17,7 +17,7 @@
 from .._poly_conversion import integer_to_poly, poly_to_integer, str_to_integer, poly_to_str, sparse_poly_to_integer, sparse_poly_to_str, str_to_sparse_poly
 
 from ._dtypes import DTYPES
-from ._linalg import dot, row_reduce, lu_decompose, lup_decompose
+from ._linalg import dot, row_reduce, lu_decompose, plu_decompose
 from ._functions import FunctionMeta
 from ._ufuncs import UfuncMeta
 
@@ -1917,34 +1917,41 @@ def lu_decompose(self) -> "FieldArray":
         """
         return lu_decompose(self)
 
-    def lup_decompose(self) -> "FieldArray":
+    def plu_decompose(self) -> "FieldArray":
         r"""
         Decomposes the input array into the product of lower and upper triangular matrices using partial pivoting.
 
         Returns
         -------
+        galois.FieldArray
+            The column permutation matrix.
         galois.FieldArray
             The lower triangular matrix.
         galois.FieldArray
             The upper triangular matrix.
-        galois.FieldArray
-            The permutation matrix.
+
+        Notes
+        -----
+        The PLU decomposition of :math:`\mathbf{A}` is defined as :math:`\mathbf{A} = \mathbf{P} \mathbf{L} \mathbf{U}`. This is equivalent to
+        :math:`\mathbf{P}^T \mathbf{A} = \mathbf{L} \mathbf{U}`.
 
         Examples
         --------
         .. ipython:: python
 
             GF = galois.GF(5)
             A = GF([[1, 3, 2, 0], [3, 4, 2, 3], [0, 2, 1, 4], [4, 3, 3, 1]])
-            L, U, P = A.lup_decompose()
+            P, L, U = A.plu_decompose()
+            P
             L
             U
-            P
 
-            # P A = L U
-            np.array_equal(P @ A, L @ U)
+            # A = P L U
+            np.array_equal(A, P @ L @ U)
+            # P.T A = L U
+            np.array_equal(P.T @ A, L @ U)
         """
-        return lup_decompose(self)
+        return plu_decompose(self)
 
     def field_trace(self) -> "FieldArray":
         r"""

scripts/generate_field_test_vectors.py

@@ -385,6 +385,27 @@ def make_luts(field, sub_folder, seed, sparse=False):
     d = {"X": X, "L": L, "U": U}
     save_pickle(d, folder, "lu_decompose.pkl")
 
+    set_seed(seed + 204)
+    shapes = [(2,2), (3,3), (4,4), (5,5), (6,6)]
+    X = []
+    L = []
+    U = []
+    P = []
+    for i in range(len(shapes)):
+        x = randint_matrix(0, order, shapes[i])
+        dtype = x.dtype
+        X.append(x)
+        x = matrix(FIELD, [[F(e) for e in row] for row in x])
+        p, l, u = x.LU()
+        p = np.array([[I(e) for e in row] for row in p], dtype)
+        l = np.array([[I(e) for e in row] for row in l], dtype)
+        u = np.array([[I(e) for e in row] for row in u], dtype)
+        P.append(p)
+        L.append(l)
+        U.append(u)
+    d = {"X": X, "P": P, "L": L, "U": U}
+    save_pickle(d, folder, "plu_decompose.pkl")
+
     ###############################################################################
     # Polynomial arithmetic
     ###############################################################################

tests/fields/conftest.py

@@ -307,6 +307,20 @@ def field_lu_decompose(field_folder):
     return d
 
 
+@pytest.fixture(scope="session")
+def field_plu_decompose(field_folder):
+    GF, folder = field_folder
+    with open(os.path.join(folder, "plu_decompose.pkl"), "rb") as f:
+        print(f"Loading {f}...")
+        d = pickle.load(f)
+    d["GF"] = GF
+    d["X"] = [GF(x) for x in d["X"]]
+    d["P"] = [GF(p) for p in d["P"]]
+    d["L"] = [GF(l) for l in d["L"]]
+    d["U"] = [GF(u) for u in d["U"]]
+    return d
+
+
 ###############################################################################
 # Fixtures for arithmetic methods over finite fields
 ###############################################################################

tests/fields/test_linalg.py

@@ -223,15 +223,6 @@ def test_matmul_2d_2d(field):
 #     assert array_equal(A @ B, np.matmul(A, B))
 
 
-def test_lup_decomposition():
-    GF = galois.GF(3)
-    I = GF.Identity(3)
-    for trial in range(100):
-        A = GF.Random((3,3))
-        L, U, P = A.lup_decompose()
-        assert array_equal(L @ U, P @ A)
-
-
 def test_det_2x2(field):
     dtype = random.choice(field.dtypes)
     a = field.Random(dtype=dtype)
@@ -342,3 +333,18 @@ def test_lu_decompose(field_lu_decompose):
         assert np.array_equal(u, U[i])
         assert type(l) is GF
         assert type(u) is GF
+
+
+def test_plu_decompose(field_plu_decompose):
+    GF, X, P, L, U = field_plu_decompose["GF"], field_plu_decompose["X"], field_plu_decompose["P"], field_plu_decompose["L"], field_plu_decompose["U"]
+
+    for i in range(len(X)):
+        dtype = random.choice(GF.dtypes)
+        x = X[i].astype(dtype)
+        p, l, u = x.plu_decompose()
+        assert np.array_equal(p, P[i])
+        assert np.array_equal(l, L[i])
+        assert np.array_equal(u, U[i])
+        assert type(p) is GF
+        assert type(l) is GF
+        assert type(u) is GF