Python-код, который реализует LU-разложение Я нашел следующее

Question

В моем поиске кода Python, который реализует декомпозицию LU, я обнаружил следующее.У меня два вопроса:

1. Мне интересно, использует ли этот код частичное вращение или нет;Я ищу тот, который не использует частичное вращение.
2. Когда я запускаю этот код, я получаю следующие ошибки:

Traceback:

IndexError                                Traceback (most recent call last)
<ipython-input-13-88eb5f3643e3> in <module>()
 60 
 61 A = [[7, 3, -1, 2], [3, 8, 1, -4], [-1, 1, 4, -1], [2, -4, -1, 6]]
---> 62 P, L, U = lu_decomposition(A)
 63 
 64 print ("A:")

<ipython-input-13-88eb5f3643e3> in lu_decomposition(A)
 54         for i in range(j, n):
 55             s2 = sum(U[k][j] * L[i][k] for k in range(j))
---> 56             L[i][j] = (PA[i][j] - s2) / U[j][j]
 57 
 58     return (P, L, U)

 IndexError: list index out of range

А вот код, который я использую

import numpy as np
import pprint

def mult_matrix(M, N):
    """Multiply square matrices of same dimension M and N"""

    # Converts N into a list of tuples of columns                                                                                                                                                                                                      
    tuple_N = zip(*N)

    # Nested list comprehension to calculate matrix multiplication                                                                                                                                                                                     
    return [[sum(el_m * el_n for el_m, el_n in zip(row_m, col_n)) for col_n in tuple_N] for row_m in M]

def pivot_matrix(M):
    """Returns the pivoting matrix for M, used in Doolittle's method."""
    m = len(M)

    # Create an identity matrix, with floating point values                                                                                                                                                                                            
    id_mat = [[float(i ==j) for i in range(m)] for j in range(m)]

    # Rearrange the identity matrix such that the largest element of                                                                                                                                                                                   
    # each column of M is placed on the diagonal of of M                                                                                                                                                                                               
    for j in range(m):
        row = max(range(j, m), key=lambda i: abs(M[i][j]))
        if j != row:
            # Swap the rows                                                                                                                                                                                                                            
            id_mat[j], id_mat[row] = id_mat[row], id_mat[j]

    return id_mat

def lu_decomposition(A):
    """Performs an LU Decomposition of A (which must be     square)                                                                                                                                                                                        
    into PA = LU. The function returns P, L and U."""
    n = len(A)

    # Create zero matrices for L and U                                                                                                                                                                                                                 
    L = [[0.0] * n for i in range(n)]
    U = [[0.0] * n for i in range(n)]

    # Create the pivot matrix P and the multipled matrix PA                                                                                                                                                                                            
    P = pivot_matrix(A)
    PA = mult_matrix(P, A)

    # Perform the LU Decomposition                                                                                                                                                                                                                     
    for j in range(n):
    # All diagonal entries of L are set to unity                                                                                                                                                                                                   
        L[j][j] = 1.0

        # LaTeX: u_{ij} = a_{ij} - \sum_{k=1}^{i-1} u_{kj} l_{ik}                                                                                                                                                                                      
        for i in range(j+1):
            s1 = sum(U[k][j] * L[i][k] for k in range(i))
            U[i][j] = PA[i][j] - s1

        # LaTeX: l_{ij} = \frac{1}{u_{jj}} (a_{ij} - \sum_{k=1}^{j-1} u_{kj} l_{ik} )                                                                                                                                                                  
        for i in range(j, n):
            s2 = sum(U[k][j] * L[i][k] for k in range(j))
            L[i][j] = (PA[i][j] - s2) / U[j][j]

    return (P, L, U)


A = [[7, 3, -1, 2], [3, 8, 1, -4], [-1, 1, 4, -1], [2, -4, -1, 6]]
P, L, U = lu_decomposition(A)

print ("A:")
pprint.pprint(A)
print ("P:")
pprint.pprint(P)

print ("L:")
pprint.pprint(L)

print ("U:")
pprint.pprint(U)

Answer 1

Проблема очевидна при небольшом диагностическом использовании print:

    for i in range(j, n):
        s2 = sum(U[k][j] * L[i][k] for k in range(j))
        print("size PA:", len(PA), len(PA[i]), i, j)
        print("PA", PA[i][j])
        print("size L :", len(L),  len(L[i]),  i, j)
        print("L",  L[i][j])
        print("size U :", len(U),  len(U[i]),  i, j)
        print("U",  U[j][i])
        L[i][j] = (PA[i][j] - s2) / U[j][j]

Вывод:

size PA: 4 4 0 0
PA 7.0
size L : 4 4 0 0
L 1.0
size U : 4 4 0 0
U 7.0
size PA: 4 0 1 0
Traceback (most recent call last):
  File "so.py", line 68, in <module>
    P, L, U = lu_decomposition(A)
  File "so.py", line 57, in lu_decomposition
    print("PA", PA[i][j])
IndexError: list index out of range

Вы попадаете в точку, где PA имеет 4 строки,но первый ряд пустВы не можете индексировать местоположение 0 пустого списка.

Что касается решения проблемы, вам необходимо диагностировать, как работает алгоритм, и где ваша реализация не соответствует данному алгоритму.Я, например, не собираюсь перепроектировать алгоритм из кода, который едва документирован и имеет дело с одно- и двухбуквенными сокращениями для анонимных понятий.Как говорится в правилах размещения сообщений, нам будет проще помочь.