FilesExpand file tree

master

/

lu_decomposition.py

Copy path

Blame

More file actions

Blame

More file actions

Latest commit

 

History

History

History

67 lines (57 loc) · 1.86 KB

master

/

lu_decomposition.py

Copy path

Top

File metadata and controls

• Code
• Blame

67 lines (57 loc) · 1.86 KB

Raw

Copy raw file

Download raw file

Open symbols panel

Edit and raw actions

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

"""Lower-Upper (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:

>>> matrix = np.array([[2, -2, 1], [0, 1, 2], [5, 3, 1]])

>>> outcome = lower_upper_decomposition(matrix)

>>> outcome[0]

array([[1. , 0. , 0. ],

[0. , 1. , 0. ],

[2.5, 8. , 1. ]])

>>> outcome[1]

array([[ 2. , -2. , 1. ],

[ 0. , 1. , 2. ],

[ 0. , 0. , -17.5]])

>>> matrix = np.array([[2, -2, 1], [0, 1, 2]])

>>> 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]]

"""

# Table that contains our data

# Table has to be a square array so we need to check first

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}"

)

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]

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]

upper[i][j] = table[i][j] - total

return lower, upper

if __name__ == "__main__":

import doctest

doctest.testmod()