Получение разных результатов при использовании разных (но правильных) алгоритмов для поиска LU-декомпозиции в Python. Оба алгоритма отлично работают в Java

#java #python-3.x #algorithm #graph-algorithm

#java #python-3.x #алгоритм #граф-алгоритм

Вопрос:

 Matrix:
2  -1  -1
-1   3  -1
-1  -1  3
  

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

Я использую два алгоритма:

  1. Алгоритм Дулиттла (этот дает правильный ответ)
  2. Алгоритм, приведенный в «книге Введение в алгоритмы» (этот не является. Но отлично работают в Java).

Реализация алгоритма Дулиттла

 def luDecomposition(mat, n):
    lower = [[0.0 for _ in range(n)]
             for _ in range(n)]
    upper = [[0.0 for _ in range(n)]
             for _ in range(n)]

    # Decomposing matrix into Upper
    # and Lower triangular matrix
    for i in range(n):

        # Upper Triangular
        for k in range(i, n):

            # Summation of L(i, j) * U(j, k)
            total = 0.0
            for j in range(i):
                total  = (lower[i][j] * upper[j][k])

                # Evaluating U(i, k)
            upper[i][k] = mat[i][k] - total

            # Lower Triangular
        for k in range(i, n):
            if i == k:
                lower[i][i] = 1.0  # Diagonal as 1
            else:

                # Summation of L(k, j) * U(j, i)
                total = 0.0
                for j in range(i):
                    total  = (lower[k][j] * upper[j][i])

                    # Evaluating L(k, i)
                lower[k][i] = (mat[k][i] - total) / upper[i][i]

    return lower, upper
  

Вывод:

 Lower Triangular Matrix:      Upper Triangular Matrix:
[[ 1.   0.   0. ]                [[ 2.  -1.  -1. ]
 [-0.5  1.   0. ]                 [ 0.   2.5 -1.5]
 [-0.5 -0.6  1. ]]                [ 0.   0.   1.6]]
  
  1. Реализация в соответствии с книгой.

     def luDecomposition(mat, n):
    
        # Initializing the matrices
        lower = np.array([[0.0]*n for _ in range(n)])
        for i in range(n):
            lower[i][i] = 1.0
    
        upper = np.array([[0.0]*n for _ in range(n)])
    
        # Decomposing matrix into Upper
        # and Lower triangular matrix
        for k in range(n):
            upper[k][k] = mat[k][k]
    
            for a in range(k   1, n):
                lower[a][k] = mat[a][k] / upper[k][k]
                upper[k][a] = mat[k][a]
    
            for x in range(k   1, n):
                for y in range(k   1, n):
                    mat[x][y] -= lower[x][k] * upper[k][y]
    
        return lower, upper
      

Вывод:

 Lower Triangular Matrix:       Upper Triangular Matrix:
[[ 1.   0.   0. ]              [[ 2. -1. -1.]
 [-0.5  1.   0. ]               [ 0.  2. -1.]
 [-0.5 -0.5  1. ]]              [ 0.  0.  1.]]
  

Я проверял несколько раз, и коды реализованы правильно.

Только если вам интересно: Это часть задания по нахождению количества остовных деревьев графа с использованием теоремы Кирхгофа. На шаге 4 теоремы кофактор вычисляется с использованием LU-декомпозиции для сокращения вычислений. Полный код можно найти здесь

Комментарии:

1. Я не могу воспроизвести это. Вы уверены, что приняли во внимание, что массив (который вы передаете в качестве аргумента) видоизменен и, следовательно, уже не тот после вызова? Может быть, вы последовательно запустили два алгоритма без повторной инициализации исходного массива?

2. @trincot Нет, я убедился, что массив был изменяемым. Я сам ответил на вопрос, указав на проблему.

Ответ №1:

Я получил это решение от reddit. Оригинальный ответ и последующие вопросы здесь


Я передавал матрицу как массив numpy, поэтому я получил неправильный результат.

 luDecomposition(np.array([[2,-1,-1],[-1,3,-1],[-1,-1,3]]), 3)
  

выдает неверный результат, но

 luDecomposition([[2,-1,-1],[-1,3,-1],[-1,-1,3]], 3)
  

выдает правильный. Это потому, что numpy принимает np.array([[2,-1,-1],[-1,3,-1],[-1,-1,3]]) быть матрицей целых чисел, и оператор -= не изменит этот тип данных:

Чтобы справиться с этим, вы могли бы сделать

 mat = np.asarray(mat, dtype=float)  # or complex or whatever
  

однако изменение матрицы, которая была передана в качестве аргумента, подобного этому, как правило, является плохой практикой (неожиданные побочные эффекты, как правило, плохие), поэтому, вероятно, лучше всего сделать копию с чем-то вроде:

 mat = np.array(mat, dtype=float)
mat = np.asarray(mat, dtype=float).copy()
  

На самом деле я не уверен, какой из них лучше, но я бы рискнул предположить, что второй быстрее, если mat уже является массивом соответствующего типа (но в противном случае он копирует дважды).


Причина, по которой алгоритм Дулиттла дал правильный ответ, заключается в том, что он не изменял матрицу, а другой алгоритм изменил.