LUP (PLU) не удалось разложить со случайной матрицей

Question

Я пытаюсь закодировать факторизацию LUP (или PLU, это то же самое) в python.У меня есть код, который работает для небольшой матрицы (под размером 4х4).Однако, когда я попробовал это с случайно сгенерированной матрицей, декомпозиция не удалась.

import numpy as np


def LUP_factorisation(A):
    """Find P, L and U : PA = LU"""
    U = A.copy()
    shape_a = U.shape
    n = shape_a[0]
    L = np.eye(n)
    P = np.eye(n)
    for i in range(n):
        print(U)
        k = i
        comp = abs(U[i, i])
        for j in range(i, n):
            if abs(U[j, i]) > comp:
                k = j
                comp = abs(U[j, i])
        line_u = U[k, :].copy()
        U[k, :] = U[i, :]
        U[i, :] = line_u
        print(U)
        line_p = P[k, :].copy()
        P[k, :] = P[i, :]
        P[i, :] = line_p
        for j in range(i + 1, n):
            g = U[j, i] / U[i, i]
            L[j, i] = g
            U[j, :] -= g * U[i, :]
    return L, U, P


if __name__ == "__main__":
    A = np.array(
        [[1.0, 2.2, 58, 9.5, 42.65], [6.56, 58.789954, 4.45, 23.465, 6.165], [7.84516, 8.9864, 96.546, 4.654, 7.6514],
         [45.65, 47.985, 1.56, 3.9845, 8.6], [455.654, 102.615, 63.965, 5.6, 9.456]])
    L, U, P = LUP_factorisation(A)
    print(L @ U)
    print(P @ A)

С примером, который я дал, это работает: у нас PA = LU.Но когда я делаю, например:

A = np.random.rand(10, 10)

Тогда я не получаю хорошего результата, потому что PA отличается от LU.Есть идеи ?Спасибо.

Answer 1

Как пишет @MattTimmermans, вы должны поменять местами строки как в L, так и в U.

Обычно это неявно обрабатывается путем сохранения LU в A, а затем перестановки автоматически применяются как к L, так и к U. См. https://en.wikipedia.org/wiki/LU_decomposition#C_code_examples

Но вы разделили их, поэтому вам нужно добавить

    line_l = L[k, :].copy()
    L[k, :] = L[i, :]
    L[i, :] = line_l

Только тестирование по диагонали с доминирующими матрицами действительно плохо; и только тестирование процедур линейной алгебры со случайными матрицами, как известно, является плохим, поскольку их свойства очень специфичны, а не «случайны». См. Работу Trefethen и его учеников, например, http://dspace.mit.edu/handle/1721.1/14322

Цель тестирования должна заключаться в том, чтобы находить ошибки, а не в том, чтобы тестовые случаи были настолько простыми, чтобы это работало.

Answer 2

Убедитесь, что диагональ входной матрицы A является доминирующей.Поэтому добавьте некоторое значение к диагонали А, например,

A = A + np.eye(A.shape)

или

A = A + 100* np.eye(A.shape)

Надеюсь, это поможет!