Fix mypy errors in lu_decomposition.py (attempt 2)

arithmetic_analysis/lu_decomposition.py

@@ -1,62 +1,101 @@
-"""Lower-Upper (LU) Decomposition.
+"""
+Lower–upper (LU) decomposition factors a matrix as a product of a lower
+triangular matrix and an upper triangular matrix. A square matrix has an LU
+decomposition under the following conditions:
+ - If the matrix is invertible, then it has an LU decomposition if and only
+ if all of its leading principal minors are non-zero (see
+ https://en.wikipedia.org/wiki/Minor_(linear_algebra) for an explanation of
+ leading principal minors of a matrix).
+ - If the matrix is singular (i.e., not invertible) and it has a rank of k
+ (i.e., it has k linearly independent columns), then it has an LU
+ decomposition if its first k leading principal minors are non-zero.
+
+This algorithm will simply attempt to perform LU decomposition on any square
+matrix and raise an error if no such decomposition exists.
 
-Reference:
-- https://en.wikipedia.org/wiki/LU_decomposition
+Reference: https://en.wikipedia.org/wiki/LU_decomposition
 """
 from __future__ import annotations
 
 import numpy as np
-from numpy import float64
-from numpy.typing import ArrayLike
-
 
-def lower_upper_decomposition(
- table: ArrayLike[float64],
-) -> tuple[ArrayLike[float64], ArrayLike[float64]]:
- """Lower-Upper (LU) Decomposition
-
- Example:
 
+def lower_upper_decomposition(table: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+ """
+ Perform LU decomposition on a given matrix and raises an error if the matrix
+ isn't square or if no such decomposition exists
  >>> matrix = np.array([[2, -2, 1], [0, 1, 2], [5, 3, 1]])
- >>> outcome = lower_upper_decomposition(matrix)
- >>> outcome[0]
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
  array([[1. , 0. , 0. ],
  [0. , 1. , 0. ],
  [2.5, 8. , 1. ]])
- >>> outcome[1]
+ >>> upper_mat
  array([[ 2. , -2. , 1. ],
  [ 0. , 1. , 2. ],
  [ 0. , 0. , -17.5]])
 
+ >>> matrix = np.array([[4, 3], [6, 3]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
+ array([[1. , 0. ],
+ [1.5, 1. ]])
+ >>> upper_mat
+ array([[ 4. , 3. ],
+ [ 0. , -1.5]])
+
+ # Matrix is not square
  >>> matrix = np.array([[2, -2, 1], [0, 1, 2]])
- >>> lower_upper_decomposition(matrix)
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
  Traceback (most recent call last):
  ...
  ValueError: 'table' has to be of square shaped array but got a 2x3 array:
  [[ 2 -2 1]
  [ 0 1 2]]
+
+ # Matrix is invertible, but its first leading principal minor is 0
+ >>> matrix = np.array([[0, 1], [1, 0]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ Traceback (most recent call last):
+ ...
+ ArithmeticError: No LU decomposition exists
+
+ # Matrix is singular, but its first leading principal minor is 1
+ >>> matrix = np.array([[1, 0], [1, 0]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
+ array([[1., 0.],
+ [1., 1.]])
+ >>> upper_mat
+ array([[1., 0.],
+ [0., 0.]])
+
+ # Matrix is singular, but its first leading principal minor is 0
+ >>> matrix = np.array([[0, 1], [0, 1]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ Traceback (most recent call last):
+ ...
+ ArithmeticError: No LU decomposition exists
  """
- # Table that contains our data
- # Table has to be a square array so we need to check first
+ # Ensure that table is a square array
  rows, columns = np.shape(table)
  if rows != columns:
  raise ValueError(
- f"'table' has to be of square shaped array but got a {rows}x{columns} "
-+ f"array:\n{table}"
+ f"'table' has to be of square shaped array but got a "
+f"{rows}x{columns} array:\n{table}"
  )
+
  lower = np.zeros((rows, columns))
  upper = np.zeros((rows, columns))
  for i in range(columns):
  for j in range(i):
- total = 0
-for k in range(j):
-total += lower[i][k] * upper[k][j]
+ total = sum(lower[i][k] * upper[k][j] for k in range(j))
+if upper[j][j] == 0:
+raise ArithmeticError("No LU decomposition exists")
  lower[i][j] = (table[i][j] - total) / upper[j][j]
  lower[i][i] = 1
  for j in range(i, columns):
- total = 0
- for k in range(i):
- total += lower[i][k] * upper[k][j]
+ total = sum(lower[i][k] * upper[k][j] for k in range(j))
  upper[i][j] = table[i][j] - total
  return lower, upper
 

arithmetic_analysis/lu_decomposition.py

@@ -1,62 +1,101 @@
-"""Lower-Upper (LU) Decomposition.
+"""
+Lower–upper (LU) decomposition factors a matrix as a product of a lower
+triangular matrix and an upper triangular matrix. A square matrix has an LU
+decomposition under the following conditions:
+ - If the matrix is invertible, then it has an LU decomposition if and only
+ if all of its leading principal minors are non-zero (see
+ https://en.wikipedia.org/wiki/Minor_(linear_algebra) for an explanation of
+ leading principal minors of a matrix).
+ - If the matrix is singular (i.e., not invertible) and it has a rank of k
+ (i.e., it has k linearly independent columns), then it has an LU
+ decomposition if its first k leading principal minors are non-zero.
+
+This algorithm will simply attempt to perform LU decomposition on any square
+matrix and raise an error if no such decomposition exists.
 
-Reference:
-- https://en.wikipedia.org/wiki/LU_decomposition
+Reference: https://en.wikipedia.org/wiki/LU_decomposition
 """
 from __future__ import annotations
 
 import numpy as np
-from numpy import float64
-from numpy.typing import ArrayLike
-
 
-def lower_upper_decomposition(
- table: ArrayLike[float64],
-) -> tuple[ArrayLike[float64], ArrayLike[float64]]:
- """Lower-Upper (LU) Decomposition
-
- Example:
 
+def lower_upper_decomposition(table: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+ """
+ Perform LU decomposition on a given matrix and raises an error if the matrix
+ isn't square or if no such decomposition exists
  >>> matrix = np.array([[2, -2, 1], [0, 1, 2], [5, 3, 1]])
- >>> outcome = lower_upper_decomposition(matrix)
- >>> outcome[0]
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
  array([[1. , 0. , 0. ],
  [0. , 1. , 0. ],
  [2.5, 8. , 1. ]])
- >>> outcome[1]
+ >>> upper_mat
  array([[ 2. , -2. , 1. ],
  [ 0. , 1. , 2. ],
  [ 0. , 0. , -17.5]])
 
+ >>> matrix = np.array([[4, 3], [6, 3]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
+ array([[1. , 0. ],
+ [1.5, 1. ]])
+ >>> upper_mat
+ array([[ 4. , 3. ],
+ [ 0. , -1.5]])
+
+ # Matrix is not square
  >>> matrix = np.array([[2, -2, 1], [0, 1, 2]])
- >>> lower_upper_decomposition(matrix)
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
  Traceback (most recent call last):
  ...
  ValueError: 'table' has to be of square shaped array but got a 2x3 array:
  [[ 2 -2 1]
  [ 0 1 2]]
+
+ # Matrix is invertible, but its first leading principal minor is 0
+ >>> matrix = np.array([[0, 1], [1, 0]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ Traceback (most recent call last):
+ ...
+ ArithmeticError: No LU decomposition exists
+
+ # Matrix is singular, but its first leading principal minor is 1
+ >>> matrix = np.array([[1, 0], [1, 0]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
+ array([[1., 0.],
+ [1., 1.]])
+ >>> upper_mat
+ array([[1., 0.],
+ [0., 0.]])
+
+ # Matrix is singular, but its first leading principal minor is 0
+ >>> matrix = np.array([[0, 1], [0, 1]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ Traceback (most recent call last):
+ ...
+ ArithmeticError: No LU decomposition exists
  """
- # Table that contains our data
- # Table has to be a square array so we need to check first
+ # Ensure that table is a square array
  rows, columns = np.shape(table)
  if rows != columns:
  raise ValueError(
- f"'table' has to be of square shaped array but got a {rows}x{columns} "
-+ f"array:\n{table}"
+ f"'table' has to be of square shaped array but got a "
+f"{rows}x{columns} array:\n{table}"
  )
+
  lower = np.zeros((rows, columns))
  upper = np.zeros((rows, columns))
  for i in range(columns):
  for j in range(i):
- total = 0
-for k in range(j):
-total += lower[i][k] * upper[k][j]
+ total = sum(lower[i][k] * upper[k][j] for k in range(j))
+if upper[j][j] == 0:
+raise ArithmeticError("No LU decomposition exists")
  lower[i][j] = (table[i][j] - total) / upper[j][j]
  lower[i][i] = 1
  for j in range(i, columns):
- total = 0
- for k in range(i):
- total += lower[i][k] * upper[k][j]
+ total = sum(lower[i][k] * upper[k][j] for k in range(j))
  upper[i][j] = table[i][j] - total
  return lower, upper
 

arithmetic_analysis/lu_decomposition.py

@@ -1,62 +1,101 @@
-"""Lower-Upper (LU) Decomposition.
+"""
+Lower–upper (LU) decomposition factors a matrix as a product of a lower
+triangular matrix and an upper triangular matrix. A square matrix has an LU
+decomposition under the following conditions:
+ - If the matrix is invertible, then it has an LU decomposition if and only
+ if all of its leading principal minors are non-zero (see
+ https://en.wikipedia.org/wiki/Minor_(linear_algebra) for an explanation of
+ leading principal minors of a matrix).
+ - If the matrix is singular (i.e., not invertible) and it has a rank of k
+ (i.e., it has k linearly independent columns), then it has an LU
+ decomposition if its first k leading principal minors are non-zero.
+
+This algorithm will simply attempt to perform LU decomposition on any square
+matrix and raise an error if no such decomposition exists.
 
-Reference:
-- https://en.wikipedia.org/wiki/LU_decomposition
+Reference: https://en.wikipedia.org/wiki/LU_decomposition
 """
 from __future__ import annotations
 
 import numpy as np
-from numpy import float64
-from numpy.typing import ArrayLike
-
 
-def lower_upper_decomposition(
- table: ArrayLike[float64],
-) -> tuple[ArrayLike[float64], ArrayLike[float64]]:
- """Lower-Upper (LU) Decomposition
-
- Example:
 
+def lower_upper_decomposition(table: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+ """
+ Perform LU decomposition on a given matrix and raises an error if the matrix
+ isn't square or if no such decomposition exists
  >>> matrix = np.array([[2, -2, 1], [0, 1, 2], [5, 3, 1]])
- >>> outcome = lower_upper_decomposition(matrix)
- >>> outcome[0]
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
  array([[1. , 0. , 0. ],
  [0. , 1. , 0. ],
  [2.5, 8. , 1. ]])
- >>> outcome[1]
+ >>> upper_mat
  array([[ 2. , -2. , 1. ],
  [ 0. , 1. , 2. ],
  [ 0. , 0. , -17.5]])
 
+ >>> matrix = np.array([[4, 3], [6, 3]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
+ array([[1. , 0. ],
+ [1.5, 1. ]])
+ >>> upper_mat
+ array([[ 4. , 3. ],
+ [ 0. , -1.5]])
+
+ # Matrix is not square
  >>> matrix = np.array([[2, -2, 1], [0, 1, 2]])
- >>> lower_upper_decomposition(matrix)
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
  Traceback (most recent call last):
  ...
  ValueError: 'table' has to be of square shaped array but got a 2x3 array:
  [[ 2 -2 1]
  [ 0 1 2]]
+
+ # Matrix is invertible, but its first leading principal minor is 0
+ >>> matrix = np.array([[0, 1], [1, 0]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ Traceback (most recent call last):
+ ...
+ ArithmeticError: No LU decomposition exists
+
+ # Matrix is singular, but its first leading principal minor is 1
+ >>> matrix = np.array([[1, 0], [1, 0]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
+ array([[1., 0.],
+ [1., 1.]])
+ >>> upper_mat
+ array([[1., 0.],
+ [0., 0.]])
+
+ # Matrix is singular, but its first leading principal minor is 0
+ >>> matrix = np.array([[0, 1], [0, 1]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ Traceback (most recent call last):
+ ...
+ ArithmeticError: No LU decomposition exists
  """
- # Table that contains our data
- # Table has to be a square array so we need to check first
+ # Ensure that table is a square array
  rows, columns = np.shape(table)
  if rows != columns:
  raise ValueError(
- f"'table' has to be of square shaped array but got a {rows}x{columns} "
-+ f"array:\n{table}"
+ f"'table' has to be of square shaped array but got a "
+f"{rows}x{columns} array:\n{table}"
  )
+
  lower = np.zeros((rows, columns))
  upper = np.zeros((rows, columns))
  for i in range(columns):
  for j in range(i):
- total = 0
-for k in range(j):
-total += lower[i][k] * upper[k][j]
+ total = sum(lower[i][k] * upper[k][j] for k in range(j))
+if upper[j][j] == 0:
+raise ArithmeticError("No LU decomposition exists")
  lower[i][j] = (table[i][j] - total) / upper[j][j]
  lower[i][i] = 1
  for j in range(i, columns):
- total = 0
- for k in range(i):
- total += lower[i][k] * upper[k][j]
+ total = sum(lower[i][k] * upper[k][j] for k in range(j))
  upper[i][j] = table[i][j] - total
  return lower, upper
 

arithmetic_analysis/lu_decomposition.py

@@ -1,62 +1,101 @@
-"""Lower-Upper (LU) Decomposition.
+"""
+Lower–upper (LU) decomposition factors a matrix as a product of a lower
+triangular matrix and an upper triangular matrix. A square matrix has an LU
+decomposition under the following conditions:
+ - If the matrix is invertible, then it has an LU decomposition if and only
+ if all of its leading principal minors are non-zero (see
+ https://en.wikipedia.org/wiki/Minor_(linear_algebra) for an explanation of
+ leading principal minors of a matrix).
+ - If the matrix is singular (i.e., not invertible) and it has a rank of k
+ (i.e., it has k linearly independent columns), then it has an LU
+ decomposition if its first k leading principal minors are non-zero.
+
+This algorithm will simply attempt to perform LU decomposition on any square
+matrix and raise an error if no such decomposition exists.
 
-Reference:
-- https://en.wikipedia.org/wiki/LU_decomposition
+Reference: https://en.wikipedia.org/wiki/LU_decomposition
 """
 from __future__ import annotations
 
 import numpy as np
-from numpy import float64
-from numpy.typing import ArrayLike
-
 
-def lower_upper_decomposition(
- table: ArrayLike[float64],
-) -> tuple[ArrayLike[float64], ArrayLike[float64]]:
- """Lower-Upper (LU) Decomposition
-
- Example:
 
+def lower_upper_decomposition(table: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+ """
+ Perform LU decomposition on a given matrix and raises an error if the matrix
+ isn't square or if no such decomposition exists
  >>> matrix = np.array([[2, -2, 1], [0, 1, 2], [5, 3, 1]])
- >>> outcome = lower_upper_decomposition(matrix)
- >>> outcome[0]
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
  array([[1. , 0. , 0. ],
  [0. , 1. , 0. ],
  [2.5, 8. , 1. ]])
- >>> outcome[1]
+ >>> upper_mat
  array([[ 2. , -2. , 1. ],
  [ 0. , 1. , 2. ],
  [ 0. , 0. , -17.5]])
 
+ >>> matrix = np.array([[4, 3], [6, 3]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
+ array([[1. , 0. ],
+ [1.5, 1. ]])
+ >>> upper_mat
+ array([[ 4. , 3. ],
+ [ 0. , -1.5]])
+
+ # Matrix is not square
  >>> matrix = np.array([[2, -2, 1], [0, 1, 2]])
- >>> lower_upper_decomposition(matrix)
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
  Traceback (most recent call last):
  ...
  ValueError: 'table' has to be of square shaped array but got a 2x3 array:
  [[ 2 -2 1]
  [ 0 1 2]]
+
+ # Matrix is invertible, but its first leading principal minor is 0
+ >>> matrix = np.array([[0, 1], [1, 0]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ Traceback (most recent call last):
+ ...
+ ArithmeticError: No LU decomposition exists
+
+ # Matrix is singular, but its first leading principal minor is 1
+ >>> matrix = np.array([[1, 0], [1, 0]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
+ array([[1., 0.],
+ [1., 1.]])
+ >>> upper_mat
+ array([[1., 0.],
+ [0., 0.]])
+
+ # Matrix is singular, but its first leading principal minor is 0
+ >>> matrix = np.array([[0, 1], [0, 1]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ Traceback (most recent call last):
+ ...
+ ArithmeticError: No LU decomposition exists
  """
- # Table that contains our data
- # Table has to be a square array so we need to check first
+ # Ensure that table is a square array
  rows, columns = np.shape(table)
  if rows != columns:
  raise ValueError(
- f"'table' has to be of square shaped array but got a {rows}x{columns} "
-+ f"array:\n{table}"
+ f"'table' has to be of square shaped array but got a "
+f"{rows}x{columns} array:\n{table}"
  )
+
  lower = np.zeros((rows, columns))
  upper = np.zeros((rows, columns))
  for i in range(columns):
  for j in range(i):
- total = 0
-for k in range(j):
-total += lower[i][k] * upper[k][j]
+ total = sum(lower[i][k] * upper[k][j] for k in range(j))
+if upper[j][j] == 0:
+raise ArithmeticError("No LU decomposition exists")
  lower[i][j] = (table[i][j] - total) / upper[j][j]
  lower[i][i] = 1
  for j in range(i, columns):
- total = 0
- for k in range(i):
- total += lower[i][k] * upper[k][j]
+ total = sum(lower[i][k] * upper[k][j] for k in range(j))
  upper[i][j] = table[i][j] - total
  return lower, upper
 

arithmetic_analysis/lu_decomposition.py

@@ -1,62 +1,101 @@
-"""Lower-Upper (LU) Decomposition.
+"""
+Lower–upper (LU) decomposition factors a matrix as a product of a lower
+triangular matrix and an upper triangular matrix. A square matrix has an LU
+decomposition under the following conditions:
+ - If the matrix is invertible, then it has an LU decomposition if and only
+ if all of its leading principal minors are non-zero (see
+ https://en.wikipedia.org/wiki/Minor_(linear_algebra) for an explanation of
+ leading principal minors of a matrix).
+ - If the matrix is singular (i.e., not invertible) and it has a rank of k
+ (i.e., it has k linearly independent columns), then it has an LU
+ decomposition if its first k leading principal minors are non-zero.
+
+This algorithm will simply attempt to perform LU decomposition on any square
+matrix and raise an error if no such decomposition exists.
 
-Reference:
-- https://en.wikipedia.org/wiki/LU_decomposition
+Reference: https://en.wikipedia.org/wiki/LU_decomposition
 """
 from __future__ import annotations
 
 import numpy as np
-from numpy import float64
-from numpy.typing import ArrayLike
-
 
-def lower_upper_decomposition(
- table: ArrayLike[float64],
-) -> tuple[ArrayLike[float64], ArrayLike[float64]]:
- """Lower-Upper (LU) Decomposition
-
- Example:
 
+def lower_upper_decomposition(table: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+ """
+ Perform LU decomposition on a given matrix and raises an error if the matrix
+ isn't square or if no such decomposition exists
  >>> matrix = np.array([[2, -2, 1], [0, 1, 2], [5, 3, 1]])
- >>> outcome = lower_upper_decomposition(matrix)
- >>> outcome[0]
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
  array([[1. , 0. , 0. ],
  [0. , 1. , 0. ],
  [2.5, 8. , 1. ]])
- >>> outcome[1]
+ >>> upper_mat
  array([[ 2. , -2. , 1. ],
  [ 0. , 1. , 2. ],
  [ 0. , 0. , -17.5]])
 
+ >>> matrix = np.array([[4, 3], [6, 3]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
+ array([[1. , 0. ],
+ [1.5, 1. ]])
+ >>> upper_mat
+ array([[ 4. , 3. ],
+ [ 0. , -1.5]])
+
+ # Matrix is not square
  >>> matrix = np.array([[2, -2, 1], [0, 1, 2]])
- >>> lower_upper_decomposition(matrix)
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
  Traceback (most recent call last):
  ...
  ValueError: 'table' has to be of square shaped array but got a 2x3 array:
  [[ 2 -2 1]
  [ 0 1 2]]
+
+ # Matrix is invertible, but its first leading principal minor is 0
+ >>> matrix = np.array([[0, 1], [1, 0]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ Traceback (most recent call last):
+ ...
+ ArithmeticError: No LU decomposition exists
+
+ # Matrix is singular, but its first leading principal minor is 1
+ >>> matrix = np.array([[1, 0], [1, 0]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
+ array([[1., 0.],
+ [1., 1.]])
+ >>> upper_mat
+ array([[1., 0.],
+ [0., 0.]])
+
+ # Matrix is singular, but its first leading principal minor is 0
+ >>> matrix = np.array([[0, 1], [0, 1]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ Traceback (most recent call last):
+ ...
+ ArithmeticError: No LU decomposition exists
  """
- # Table that contains our data
- # Table has to be a square array so we need to check first
+ # Ensure that table is a square array
  rows, columns = np.shape(table)
  if rows != columns:
  raise ValueError(
- f"'table' has to be of square shaped array but got a {rows}x{columns} "
-+ f"array:\n{table}"
+ f"'table' has to be of square shaped array but got a "
+f"{rows}x{columns} array:\n{table}"
  )
+
  lower = np.zeros((rows, columns))
  upper = np.zeros((rows, columns))
  for i in range(columns):
  for j in range(i):
- total = 0
-for k in range(j):
-total += lower[i][k] * upper[k][j]
+ total = sum(lower[i][k] * upper[k][j] for k in range(j))
+if upper[j][j] == 0:
+raise ArithmeticError("No LU decomposition exists")
  lower[i][j] = (table[i][j] - total) / upper[j][j]
  lower[i][i] = 1
  for j in range(i, columns):
- total = 0
- for k in range(i):
- total += lower[i][k] * upper[k][j]
+ total = sum(lower[i][k] * upper[k][j] for k in range(j))
  upper[i][j] = table[i][j] - total
  return lower, upper
 

arithmetic_analysis/lu_decomposition.py

@@ -1,62 +1,101 @@
-"""Lower-Upper (LU) Decomposition.
+"""
+Lower–upper (LU) decomposition factors a matrix as a product of a lower
+triangular matrix and an upper triangular matrix. A square matrix has an LU
+decomposition under the following conditions:
+ - If the matrix is invertible, then it has an LU decomposition if and only
+ if all of its leading principal minors are non-zero (see
+ https://en.wikipedia.org/wiki/Minor_(linear_algebra) for an explanation of
+ leading principal minors of a matrix).
+ - If the matrix is singular (i.e., not invertible) and it has a rank of k
+ (i.e., it has k linearly independent columns), then it has an LU
+ decomposition if its first k leading principal minors are non-zero.
+
+This algorithm will simply attempt to perform LU decomposition on any square
+matrix and raise an error if no such decomposition exists.
 
-Reference:
-- https://en.wikipedia.org/wiki/LU_decomposition
+Reference: https://en.wikipedia.org/wiki/LU_decomposition
 """
 from __future__ import annotations
 
 import numpy as np
-from numpy import float64
-from numpy.typing import ArrayLike
-
 
-def lower_upper_decomposition(
- table: ArrayLike[float64],
-) -> tuple[ArrayLike[float64], ArrayLike[float64]]:
- """Lower-Upper (LU) Decomposition
-
- Example:
 
+def lower_upper_decomposition(table: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+ """
+ Perform LU decomposition on a given matrix and raises an error if the matrix
+ isn't square or if no such decomposition exists
  >>> matrix = np.array([[2, -2, 1], [0, 1, 2], [5, 3, 1]])
- >>> outcome = lower_upper_decomposition(matrix)
- >>> outcome[0]
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
  array([[1. , 0. , 0. ],
  [0. , 1. , 0. ],
  [2.5, 8. , 1. ]])
- >>> outcome[1]
+ >>> upper_mat
  array([[ 2. , -2. , 1. ],
  [ 0. , 1. , 2. ],
  [ 0. , 0. , -17.5]])
 
+ >>> matrix = np.array([[4, 3], [6, 3]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
+ array([[1. , 0. ],
+ [1.5, 1. ]])
+ >>> upper_mat
+ array([[ 4. , 3. ],
+ [ 0. , -1.5]])
+
+ # Matrix is not square
  >>> matrix = np.array([[2, -2, 1], [0, 1, 2]])
- >>> lower_upper_decomposition(matrix)
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
  Traceback (most recent call last):
  ...
  ValueError: 'table' has to be of square shaped array but got a 2x3 array:
  [[ 2 -2 1]
  [ 0 1 2]]
+
+ # Matrix is invertible, but its first leading principal minor is 0
+ >>> matrix = np.array([[0, 1], [1, 0]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ Traceback (most recent call last):
+ ...
+ ArithmeticError: No LU decomposition exists
+
+ # Matrix is singular, but its first leading principal minor is 1
+ >>> matrix = np.array([[1, 0], [1, 0]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
+ array([[1., 0.],
+ [1., 1.]])
+ >>> upper_mat
+ array([[1., 0.],
+ [0., 0.]])
+
+ # Matrix is singular, but its first leading principal minor is 0
+ >>> matrix = np.array([[0, 1], [0, 1]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ Traceback (most recent call last):
+ ...
+ ArithmeticError: No LU decomposition exists
  """
- # Table that contains our data
- # Table has to be a square array so we need to check first
+ # Ensure that table is a square array
  rows, columns = np.shape(table)
  if rows != columns:
  raise ValueError(
- f"'table' has to be of square shaped array but got a {rows}x{columns} "
-+ f"array:\n{table}"
+ f"'table' has to be of square shaped array but got a "
+f"{rows}x{columns} array:\n{table}"
  )
+
  lower = np.zeros((rows, columns))
  upper = np.zeros((rows, columns))
  for i in range(columns):
  for j in range(i):
- total = 0
-for k in range(j):
-total += lower[i][k] * upper[k][j]
+ total = sum(lower[i][k] * upper[k][j] for k in range(j))
+if upper[j][j] == 0:
+raise ArithmeticError("No LU decomposition exists")
  lower[i][j] = (table[i][j] - total) / upper[j][j]
  lower[i][i] = 1
  for j in range(i, columns):
- total = 0
- for k in range(i):
- total += lower[i][k] * upper[k][j]
+ total = sum(lower[i][k] * upper[k][j] for k in range(j))
  upper[i][j] = table[i][j] - total
  return lower, upper
 

arithmetic_analysis/lu_decomposition.py

@@ -1,62 +1,101 @@
-"""Lower-Upper (LU) Decomposition.
+"""
+Lower–upper (LU) decomposition factors a matrix as a product of a lower
+triangular matrix and an upper triangular matrix. A square matrix has an LU
+decomposition under the following conditions:
+ - If the matrix is invertible, then it has an LU decomposition if and only
+ if all of its leading principal minors are non-zero (see
+ https://en.wikipedia.org/wiki/Minor_(linear_algebra) for an explanation of
+ leading principal minors of a matrix).
+ - If the matrix is singular (i.e., not invertible) and it has a rank of k
+ (i.e., it has k linearly independent columns), then it has an LU
+ decomposition if its first k leading principal minors are non-zero.
+
+This algorithm will simply attempt to perform LU decomposition on any square
+matrix and raise an error if no such decomposition exists.
 
-Reference:
-- https://en.wikipedia.org/wiki/LU_decomposition
+Reference: https://en.wikipedia.org/wiki/LU_decomposition
 """
 from __future__ import annotations
 
 import numpy as np
-from numpy import float64
-from numpy.typing import ArrayLike
-
 
-def lower_upper_decomposition(
- table: ArrayLike[float64],
-) -> tuple[ArrayLike[float64], ArrayLike[float64]]:
- """Lower-Upper (LU) Decomposition
-
- Example:
 
+def lower_upper_decomposition(table: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+ """
+ Perform LU decomposition on a given matrix and raises an error if the matrix
+ isn't square or if no such decomposition exists
  >>> matrix = np.array([[2, -2, 1], [0, 1, 2], [5, 3, 1]])
- >>> outcome = lower_upper_decomposition(matrix)
- >>> outcome[0]
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
  array([[1. , 0. , 0. ],
  [0. , 1. , 0. ],
  [2.5, 8. , 1. ]])
- >>> outcome[1]
+ >>> upper_mat
  array([[ 2. , -2. , 1. ],
  [ 0. , 1. , 2. ],
  [ 0. , 0. , -17.5]])
 
+ >>> matrix = np.array([[4, 3], [6, 3]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
+ array([[1. , 0. ],
+ [1.5, 1. ]])
+ >>> upper_mat
+ array([[ 4. , 3. ],
+ [ 0. , -1.5]])
+
+ # Matrix is not square
  >>> matrix = np.array([[2, -2, 1], [0, 1, 2]])
- >>> lower_upper_decomposition(matrix)
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
  Traceback (most recent call last):
  ...
  ValueError: 'table' has to be of square shaped array but got a 2x3 array:
  [[ 2 -2 1]
  [ 0 1 2]]
+
+ # Matrix is invertible, but its first leading principal minor is 0
+ >>> matrix = np.array([[0, 1], [1, 0]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ Traceback (most recent call last):
+ ...
+ ArithmeticError: No LU decomposition exists
+
+ # Matrix is singular, but its first leading principal minor is 1
+ >>> matrix = np.array([[1, 0], [1, 0]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ >>> lower_mat
+ array([[1., 0.],
+ [1., 1.]])
+ >>> upper_mat
+ array([[1., 0.],
+ [0., 0.]])
+
+ # Matrix is singular, but its first leading principal minor is 0
+ >>> matrix = np.array([[0, 1], [0, 1]])
+ >>> lower_mat, upper_mat = lower_upper_decomposition(matrix)
+ Traceback (most recent call last):
+ ...
+ ArithmeticError: No LU decomposition exists
  """
- # Table that contains our data
- # Table has to be a square array so we need to check first
+ # Ensure that table is a square array
  rows, columns = np.shape(table)
  if rows != columns:
  raise ValueError(
- f"'table' has to be of square shaped array but got a {rows}x{columns} "
-+ f"array:\n{table}"
+ f"'table' has to be of square shaped array but got a "
+f"{rows}x{columns} array:\n{table}"
  )
+
  lower = np.zeros((rows, columns))
  upper = np.zeros((rows, columns))
  for i in range(columns):
  for j in range(i):
- total = 0
-for k in range(j):
-total += lower[i][k] * upper[k][j]
+ total = sum(lower[i][k] * upper[k][j] for k in range(j))
+if upper[j][j] == 0:
+raise ArithmeticError("No LU decomposition exists")
  lower[i][j] = (table[i][j] - total) / upper[j][j]
  lower[i][i] = 1
  for j in range(i, columns):
- total = 0
- for k in range(i):
- total += lower[i][k] * upper[k][j]
+ total = sum(lower[i][k] * upper[k][j] for k in range(j))
  upper[i][j] = table[i][j] - total
  return lower, upper