使用 LU 分解实现 Ax = b 求解器时遇到问题

Question

我正在尝试使用 A = LU 算法求解 Ax = b 形式的线性系统，其中 A 是 nxn 实数矩阵，b 是 1xn 实数向量。

我已经实现了必要的功能，但我不确定问题出在哪个或哪些功能上。

import numpy as np

 def forward_sub(L, b):
    """Given a lower triangular matrix L and right-side vector b,
    compute the solution vector y solving Ly = b."""

    y = []
    for i in range(len(b)):
        y.append(b[i])
        for j in range(i):
            y[i]=y[i]-(L[i, j]*y[j])
        y[i] = y[i]/L[i, i]

    return y

def backward_sub(U, y):
    """Given a lower triangular matrix U and right-side vector y,
    compute the solution vector x solving Ux = y."""



    x = [0 for i in range(len(U))]

    for i in range(len(U)-1, 0, -1):
        x[i] = U[i, i]/y[i]
        for j in range (i-1, 0, -1):
            U[i, i] += U[j, i]*x[i]

    return x


def lu_factor(A):


    #LU decompostion using Doolittles method

 L = np.zeros_like(A)
 U = np.zeros_like(A)

 N = np.size(A,0)

 for k in range(N):
        L[k, k] = 1
        U[k, k] = (A[k, k] - np.dot(L[k, :k], U[:k, k])) / L[k, k]
        for j in range(k+1, N):
            U[k, j] = (A[k, j] - np.dot(L[k, :k], U[:k, j])) / L[k, k]
        for i in range(k+1, N):
            L[i, k] = (A[i, k] - np.dot(L[i, :k], U[:k, k])) / U[k, k]


  return (L, U)

def lu_solve(L, U, b):
    # Step 1: Solve Uy = b using forward substitution

    # Step 2: Solve Lx = y using backward substitution

    y = forward_sub(L,b)
    x = backward_sub(U,y)

    return x


def linear_solve(A, b):
    # ...

    L, U = lu_factor(A)
    x = lu_solve(L,U,b)
    return x


b = [6,-4,27]
A = np.matrix([[1,1,1],[0,2,5],[2,5,-1]])

print(linear_solve(A,b))

如上选择 A 和 b 给出 x = [0,-0.5,-0.42] 作为我的解向量，但是它应该给出 x = [5,3,-2]

Answer 1

A是一个整数矩阵。这也使 L 和 U 整数矩阵，但正确的结果是：

L:
[[1.  0.  0. ]
 [0.  1.  0. ]
 [2.  1.5 1. ]]
U:
[[  1.    1.    1. ]
 [  0.    2.    5. ]
 [  0.    0.  -10.5]]

需要一些小数值。通常，即使输入由整数组成，LU 分解也是如此。毕竟有一些分裂。

更改数据类型可以解决这个问题。例如：

A = np.matrix([[1.,1,1],
               [0,2,5],
               [2,5,-1]])

backward_sub 坏了，我不确定具体是怎么回事，但无论如何这是一个奇怪的实现。这个似乎有效：

def backward_sub(U, y):
    """Given a lower triangular matrix U and right-side vector y,
    compute the solution vector x solving Ux = y."""

    x = np.zeros_like(y)

    for i in range(len(x), 0, -1):
      x[i-1] = (y[i-1] - np.dot(U[i-1, i:], x[i:])) / U[i-1, i-1]

    return x

结果是[ 5. 3. -2.]，试一试ideone

使用 LU 分解实现 Ax = b 求解器时遇到问题

Trouble implementing Ax = b solver using LU decomposition

python

algorithm

linear-algebra

python-3.x