在使用不同但正确的算法查找Python中的LU分解时获得不同的结果两种算法在Java中都可以正常工作

问题描述

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

我知道有用于找到LU分解的内置模块，但这是算法分配，所以我不打算使用它们。

我正在使用两种算法：

1. Doolittle's Algorithm（此人给出正确答案）
2. Algorithm provided in "Introduction to Algorithms book"（这不是。但是在Java中工作正常）。

Doolittle的算法实现

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

我已经检查了多次，并且代码正确实现。

只有在您感兴趣的情况下： 这是使用Kirchoff's Theorem查找图的生成树数量的分配的一部分。在定理的第4步中，使用LU分解来计算辅因子以减少计算量。完整的代码可以在here

中找到

解决方法

我从reddit得到了这个解决方案。原始答案和后续问题here

---

我将矩阵作为numpy数组传递，这就是为什么我得到错误结果的原因。

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

给出错误的结果，但是

luDecomposition([[2,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已经是适当类型的数组，则第二个会更快（但否则会复制两次）。

---

Doolittle算法之所以给出正确答案的原因是因为它没有改变矩阵，而另一种算法却没有改变。