Разделить на нулевое предупреждение в работе алгоритма разложения LU - Doolittle - PullRequest
Новая
       12

Разделить на нулевое предупреждение в работе алгоритма разложения LU - Doolittle

0 голосов
/ 30 апреля 2020

Я реализовал стандартные уравнения / алгоритм LU-разложения матрицы, перейдя по этой ссылке: ( 1 ) и ( 2 ). Это возвращает LU-разложение квадратной матрицы вроде ниже отлично. Моя проблема, однако, это также дает Divide by Zero warning.

Код здесь: 

import numpy as np

def LUDecomposition (A):
    L = np.zeros(np.shape(A),np.float64)
    U = np.zeros(np.shape(A),np.float64)
    acc = 0
    L[0,0]=1
    for i in np.arange(len(A)):

        for k in range(i,len(A)):

            for j in range(0,i):
                acc += L[i,j]*U[j,k]
            U[i,k] = A[i,k]-acc


            for m in range(k+1,len(A)):
                if m==k:
                    L[m,k]=1
                else:

                    L[m,k] = (A[m,k]-acc)/U[k,k]
            acc=0
    return (L,U)

A = np.array([[-4, -1, -2],
              [-4, 12,  3],
              [-4, -2, 18]])

L, U = LUDecomposition (A)

Где я иду не так?

1 Ответ

2 голосов
/ 30 апреля 2020

Похоже, что вы допустили некоторые ошибки в отступах относительно первого внутреннего уровня for loop: U необходимо оценить до L; Вы также неправильно вычислили член суммирования acc и неправильно установили диагональные члены L в 1. После некоторых других изменений синтаксиса вы можете переписать свою функцию следующим образом: 

def LUDecomposition(A):

    n = A.shape[0]
    L = np.zeros((n,n), np.float64)
    U = np.zeros((n,n), np.float64)

    for i in range(n):
        # U
        for k in range(i,n):
            s1 = 0  # summation of L(i, j)*U(j, k) 
            for j in range(i):
                s1 += L[i,j]*U[j,k]
            U[i,k] = A[i,k] - s1

        # L
        for k in range(i,n):
            if i==k:
                # diagonal terms of L 
                L[i,i] = 1
            else:
                s2 = 0 # summation of L(k, j)*U(j, i) 
                for j in range(i):
                    s2 += L[k,j]*U[j,i]
                L[k,i] = (A[k,i] - s2)/U[i,i]

    return L, U

, который дает на этот раз правильный выход для матрицы A по сравнению с scipy.linalg.lu в качестве надежного эталона: 

import numpy as np
from scipy.linalg import lu

A = np.array([[-4, -1, -2],
              [-4, 12,  3],
              [-4, -2, 18]])

L, U = LUDecomposition(A)
P, L_sp, U_sp = lu(A, permute_l=False)

P
>>> [[1., 0., 0.],
     [0., 1., 0.],
     [0., 0., 1.]])

L
>>> [[ 1.          0.          0.        ]
     [ 1.          1.          0.        ]
     [ 1.         -0.07692308  1.        ]]

np.allclose(L_sp, L))
>>>  True

U
>>> [[-4.         -1.         -2.        ]
     [ 0.         13.          5.        ]
     [ 0.          0.         20.38461538]]

np.allclose(U_sp, U))
>>>  True

Примечание : в отличие от алгоритма scipy lapack getrf, эта реализация Doolittle не включает в себя поворот, тогда эти два сравнения верны только в том случае, если матрица перестановок P, возвращаемая scipy.linalg.lu, является единичной матрицей, то есть scipy не выполнила никаких перестановки, что действительно имеет место для вашей матрицы A. Матрица перестановок, определенная в алгоритме Сципи, предназначена для оптимизации номеров условий результирующей матрицы с целью уменьшения ошибок округления. Наконец, вы можете просто проверить, что A = LU, что всегда будет, если факторизация сделана правильно: 

A = np.random.rand(10,10)
L, U = LUDecomposition(A)

np.allclose(A, np.dot(L, U))
>>>  True

Тем не менее, с точки зрения численной эффективности и точности, я бы не рекомендовал вам использовать свою собственную функцию для вычисления разложения LU. Надеюсь, это поможет.

...