Python - 用 LU 分解求解线性方程组的数值问题
Python - Numerical issues for solving systems of linear equations with LU decomposition
信息
OS: 曼扎罗 Linux
开发环境：Spyder4 (Python3.7)
库：Numpy
问题
嗨，
我写了几个函数来按照这三种方法求解线性方程组：
- LU分解
- 雅可比方法
- 高斯-赛德尔法
程序运行完美。但是，LU 分解的结果让我很困扰。
例如，如果我的矩阵 A 和向量 b 是
A = np.array([[10., -1., 2., 0.],
[-1., 11., -1., 3.],
[2., -1., 10., -1.],
[0., 3., -1., 8.]])
b = np.array([6., 25., -11., 15.])
那么我的结果是
- LU:[ 0.95034091 1.67840909 -0.9125 1.875 ]
- 雅各比：[ 1. 2. -1. 1.]
- 高斯-赛德尔：[ 1. 2. -1. 1.]
如您所见，LU 给我的结果略有不同。这是舍入或截断错误的问题吗？任何建议或帮助将不胜感激。
干杯！
编辑 1：
我解决了这个问题，它只是 LU 分解期间代码中的错误。感谢大家的反馈。
下面的共享代码：
def lu_decomposition(matrix_in, b_in, n):
lower = [[0 for x in range(n)]
for y in range(n)];
upper = [[0 for x in range(n)]
for y in range(n)];
# Doolittle's Method
# Decomposing matrix into upper and lower matrices
for i in range(n):
# Upper triangle
for k in range(i, n):
# Sigma from j to i for each row k of U
sum = 0
for j in range(i):
sum += lower[i][j] * upper[j][k]
# Evaluate U for row k
upper[i][k] = matrix_in[i][k] - sum
# Lower Triangle
for k in range(i,n):
if(i == k): # Entry of a diagonal element
lower[i][i] = 1
else:
# Sigma from j to i for each column k of L
sum = 0
for j in range(i):
sum += (lower[k][j] * upper[j][i])
# Evaluate L for column k
lower[k][i] = int( (matrix_in[k][i] - sum)/ upper[i][i])
# Perform forward substitution Ly=b
y = [0 for x in range(n)]
lower_inv = np.linalg.inv(lower)
y = np.dot(lower_inv, b_in)
# Perform back substitution Ux=y
x_sol = [0 for x in range(n)]
upper_inv = np.linalg.inv(upper)
x_sol = np.dot(upper_inv, y)
# printing results
# setw is for displaying nicely
print("Lower Triangular\t\tUpper Triangular");
# Displaying the result :
for i in range(n):
# Lower
for j in range(n):
print(lower[i][j], end = "\t");
print("", end = "\t");
# Upper
for j in range(n):
print(upper[i][j], end = "\t");
print("");
print("Here's the solution for your matrix: ")
print(x_sol)
您的代码中有很多问题需要解决：
- Python 不需要；要结束一行，python 知道行何时结束。 Python 只使用 ;将两行代码放在同一物理行上。喜欢"x=7;print(x)"。
- 如果您不在列表推导中使用变量，那么习惯上使用 _。例如，[0 for _ in range(10)]。当然，python还有更好的办法“[0]*10”。
- 我注意到您导入了 numpy，但没有使用 numpy 数组的操作。这将极大地改进您的代码（并使其更快）并使其更易于阅读。例如，您可以同时写入和读取整行或整列。 "matrix[0,:]"（第一行），"matrix[:,0]"（第一列）。
- 无需在 python 中包含对象的长度，因为 python 会自动存储其对象的长度，始终可以使用内置的 len() 检索。
- 对您的代码来说并不重要（因为您可能正在练习这样做），但您可能已经知道，lu 分解已经存在，例如 scipy.linalg.lu()。事实上，快速检查 linalg.lu(A,True) 就会发现你的 L 和 U 都有错误。
- 看到您手动生成 L 和 U 矩阵然后继续在 L 和 U 上使用 np.linalg.inv() 真的很奇怪。如果您愿意使用 np.linalg.inv( ) 函数，那么你的答案是一行，"np.linalg.inv(A) @ b"。通常，人们更容易找到 L 和 U 来手动求解 X。要找到 L 和 U 然后使用 numpy 的反函数，有点违背了目的。
- 虽然它有时会有帮助，但 python 并不要求您在创建对象之前在内存中划出 space。 Python 自动管理您的内存创建和删除（无需创建空的零列表）。
- np.dot() 只是手动访问 python 中的“@”运算符的方法。
一些示例：
import numpy as np
from scipy import linalg
def lu_dec_solver_v1(A,b):
return np.linalg.inv(A) @ b
def lu_dec_solver_v2(A,b):
L, U = linalg.lu(A,True)
L_inv = np.linalg.inv(L)
U_inv = np.linalg.inv(U)
return U_inv @ (L_inv @ b)
def lu_dec_solver_v3(A,b):
U = A.copy()
L = np.identity(len(A))
for n in range(0,len(A)-1):
for m in range(n+1,len(A)):
L[m,n] = U[m,n]/U[n,n]
U[m,:] += -L[m,n]*U[n,:]
L_inv = np.linalg.inv(L)
U_inv = np.linalg.inv(U)
return U_inv @ (L_inv @ b)
如果您想自己实现 LU 分解，this 是一种不需要任何外部库的代码。我把它附加在这里，以防你只是想复制它
def Abs(x):
'''
Input x and you'll get the abs.
'''
#Not the best way, but it is an easy way to not use numpy
return (x**2)**0.5
def ind_max(row,N):
'''
Find the index of the maximum of a list (row) of lentgth N.
'''
_in=0
_max=row[0]
i=0
while i<N:#the end of the row should be included (convension in how I use LUP..)
if row[i]>_max:
_max=row[i]
_in=i
i+=1
return _in
def Sum(List,N):
'''
Calculates the sum of a List of size N
'''
s=0
for i in range(N):
s+=List[i]
return s
def index_swap(A,index_1,index_2):
'''
index_swap takes an array and interchanges
A[index_1] with A[index_2].
Example:
A=[0,1,2,3]
index_swap(A,0,2)
A
>>[2, 1, 0, 3]
'''
tmp=A[index_1]
A[index_1]=A[index_2]
A[index_2]=tmp
def apply_permutations_vector(A,P,N):
'''
Applies the permutations given by P from LUP
to a list A of length N, and returns the result.
Example:
A=[1,2,5,8,3]
P=[2,4,0,3,1]
apply_permutations_vector(A,P,5)
>>[5, 3, 1, 8, 2]
'''
#that is, you make a list like this (P basically gives you the indices of A):)
Ap=[A[ P[i] ] for i in range(N)]
return Ap
def apply_permutations_matrix(M,P,N):
'''
Applies the permutations given by P from LUP
to a N*N array M of length N, and returns the result.
M=[
[ 0, 2, 2 , 3 , 5],
[-3, -1, 1 , 5 , 9],
[ 1, -1, 1 , 4 , 7],
[ 1, -1, 1 , 0 , 2],
[ 1, -1, 1 , 0 , 3]
]
P=[2,0,1,4,3]
apply_permutations_matrix(M,P,5)
>>[
[ 1, -1, 1, 4, 7],
[ 0, 2, 2, 3, 5],
[-3, -1, 1, 5, 9],
[ 1, -1, 1, 0, 3],
[ 1, -1, 1, 0, 2]
]
'''
Mp=[ [M[ P[i] ][j]for j in range(N)] for i in range(N) ]
return Mp
def LUP(M,N,_tiny=1e-20):
U=[ [ M[i][j] for j in range(N)] for i in range(N) ]
L=[ [ 0 if i!=j else 1 for j in range(N)] for i in range(N) ]
#this is the "permutation vector". if it is e.g. [2 1 0 3] it means you make 0<->2
P=[ i for i in range(N) ]
for k in range(1,N):
for i in range(k,N):
#find the index of the maximum in column
_col=[Abs(U[_r][k-1]) for _r in range(k-1,N)]
#find the index of the maximum of _col
# notice that the length of _col is N-(k-1)
len_col=N-(k-1)
pivot=ind_max( _col ,len_col) + k - 1 #convert the index of _col (it has a length of len_col) to the index of a row of U
##################################################
#if you remove it, then you get a lot of infinities
#it has to do with the fact that if U[pivot][k-1] <_tiny , then U[k-1][k-1] will be a zero,
#L[i][k-1] explodes.
#You are allowed to skip this i, then, because if U[pivot][k-1] <_tiny , then all U[i][k-1] are small!
#Check that this is true by uncommenting print(_col)
if Abs(U[pivot][k-1]) < _tiny :
#print(_col)
break
###################################################
#if the maximum is not at k-1, swap!
if pivot != k-1 :
# Permute rows k-1 and pivot in U
index_swap(P,k-1,pivot)
tmpU=[U[k-1][_r] for _r in range(k-1,N)]
#print(U)
for _r in range(k-1,N):
U[k-1][_r]=U[pivot][_r]
#print(U)
for _r in range(k-1,N):
U[pivot][_r]=tmpU[_r-(k-1)]#again we have to convert the index of tmpU
#print(U)
#print("=========================")
tmpL=[L[k-1][_r] for _r in range(k-1)]
#print(L)
for _r in range(k-1):
L[k-1][_r]=L[pivot][_r]
#print(L)
for _r in range(k-1):
L[pivot][_r]=tmpL[_r]
#print(L)
#print("========================")
L[i][k-1]=U[i][k-1]/U[k-1][k-1]
for j in range(k-1,N):
U[i][j]=U[i][j]-L[i][k-1]*U[k-1][j]
return L,U,P
信息
OS: 曼扎罗 Linux
开发环境：Spyder4 (Python3.7)
库：Numpy
问题
嗨，
我写了几个函数来按照这三种方法求解线性方程组：
- LU分解
- 雅可比方法
- 高斯-赛德尔法
程序运行完美。但是，LU 分解的结果让我很困扰。 例如，如果我的矩阵 A 和向量 b 是
A = np.array([[10., -1., 2., 0.],
[-1., 11., -1., 3.],
[2., -1., 10., -1.],
[0., 3., -1., 8.]])
b = np.array([6., 25., -11., 15.])
那么我的结果是
- LU:[ 0.95034091 1.67840909 -0.9125 1.875 ]
- 雅各比：[ 1. 2. -1. 1.]
- 高斯-赛德尔：[ 1. 2. -1. 1.]
如您所见，LU 给我的结果略有不同。这是舍入或截断错误的问题吗？任何建议或帮助将不胜感激。
干杯！
编辑 1：
我解决了这个问题，它只是 LU 分解期间代码中的错误。感谢大家的反馈。
下面的共享代码：
def lu_decomposition(matrix_in, b_in, n):
lower = [[0 for x in range(n)]
for y in range(n)];
upper = [[0 for x in range(n)]
for y in range(n)];
# Doolittle's Method
# Decomposing matrix into upper and lower matrices
for i in range(n):
# Upper triangle
for k in range(i, n):
# Sigma from j to i for each row k of U
sum = 0
for j in range(i):
sum += lower[i][j] * upper[j][k]
# Evaluate U for row k
upper[i][k] = matrix_in[i][k] - sum
# Lower Triangle
for k in range(i,n):
if(i == k): # Entry of a diagonal element
lower[i][i] = 1
else:
# Sigma from j to i for each column k of L
sum = 0
for j in range(i):
sum += (lower[k][j] * upper[j][i])
# Evaluate L for column k
lower[k][i] = int( (matrix_in[k][i] - sum)/ upper[i][i])
# Perform forward substitution Ly=b
y = [0 for x in range(n)]
lower_inv = np.linalg.inv(lower)
y = np.dot(lower_inv, b_in)
# Perform back substitution Ux=y
x_sol = [0 for x in range(n)]
upper_inv = np.linalg.inv(upper)
x_sol = np.dot(upper_inv, y)
# printing results
# setw is for displaying nicely
print("Lower Triangular\t\tUpper Triangular");
# Displaying the result :
for i in range(n):
# Lower
for j in range(n):
print(lower[i][j], end = "\t");
print("", end = "\t");
# Upper
for j in range(n):
print(upper[i][j], end = "\t");
print("");
print("Here's the solution for your matrix: ")
print(x_sol)
您的代码中有很多问题需要解决：
- Python 不需要；要结束一行，python 知道行何时结束。 Python 只使用 ;将两行代码放在同一物理行上。喜欢"x=7;print(x)"。
- 如果您不在列表推导中使用变量，那么习惯上使用 _。例如，[0 for _ in range(10)]。当然，python还有更好的办法“[0]*10”。
- 我注意到您导入了 numpy，但没有使用 numpy 数组的操作。这将极大地改进您的代码（并使其更快）并使其更易于阅读。例如，您可以同时写入和读取整行或整列。 "matrix[0,:]"（第一行），"matrix[:,0]"（第一列）。
- 无需在 python 中包含对象的长度，因为 python 会自动存储其对象的长度，始终可以使用内置的 len() 检索。
- 对您的代码来说并不重要（因为您可能正在练习这样做），但您可能已经知道，lu 分解已经存在，例如 scipy.linalg.lu()。事实上，快速检查 linalg.lu(A,True) 就会发现你的 L 和 U 都有错误。
- 看到您手动生成 L 和 U 矩阵然后继续在 L 和 U 上使用 np.linalg.inv() 真的很奇怪。如果您愿意使用 np.linalg.inv( ) 函数，那么你的答案是一行，"np.linalg.inv(A) @ b"。通常，人们更容易找到 L 和 U 来手动求解 X。要找到 L 和 U 然后使用 numpy 的反函数，有点违背了目的。
- 虽然它有时会有帮助，但 python 并不要求您在创建对象之前在内存中划出 space。 Python 自动管理您的内存创建和删除（无需创建空的零列表）。
- np.dot() 只是手动访问 python 中的“@”运算符的方法。
一些示例：
import numpy as np
from scipy import linalg
def lu_dec_solver_v1(A,b):
return np.linalg.inv(A) @ b
def lu_dec_solver_v2(A,b):
L, U = linalg.lu(A,True)
L_inv = np.linalg.inv(L)
U_inv = np.linalg.inv(U)
return U_inv @ (L_inv @ b)
def lu_dec_solver_v3(A,b):
U = A.copy()
L = np.identity(len(A))
for n in range(0,len(A)-1):
for m in range(n+1,len(A)):
L[m,n] = U[m,n]/U[n,n]
U[m,:] += -L[m,n]*U[n,:]
L_inv = np.linalg.inv(L)
U_inv = np.linalg.inv(U)
return U_inv @ (L_inv @ b)
如果您想自己实现 LU 分解，this 是一种不需要任何外部库的代码。我把它附加在这里，以防你只是想复制它
def Abs(x):
'''
Input x and you'll get the abs.
'''
#Not the best way, but it is an easy way to not use numpy
return (x**2)**0.5
def ind_max(row,N):
'''
Find the index of the maximum of a list (row) of lentgth N.
'''
_in=0
_max=row[0]
i=0
while i<N:#the end of the row should be included (convension in how I use LUP..)
if row[i]>_max:
_max=row[i]
_in=i
i+=1
return _in
def Sum(List,N):
'''
Calculates the sum of a List of size N
'''
s=0
for i in range(N):
s+=List[i]
return s
def index_swap(A,index_1,index_2):
'''
index_swap takes an array and interchanges
A[index_1] with A[index_2].
Example:
A=[0,1,2,3]
index_swap(A,0,2)
A
>>[2, 1, 0, 3]
'''
tmp=A[index_1]
A[index_1]=A[index_2]
A[index_2]=tmp
def apply_permutations_vector(A,P,N):
'''
Applies the permutations given by P from LUP
to a list A of length N, and returns the result.
Example:
A=[1,2,5,8,3]
P=[2,4,0,3,1]
apply_permutations_vector(A,P,5)
>>[5, 3, 1, 8, 2]
'''
#that is, you make a list like this (P basically gives you the indices of A):)
Ap=[A[ P[i] ] for i in range(N)]
return Ap
def apply_permutations_matrix(M,P,N):
'''
Applies the permutations given by P from LUP
to a N*N array M of length N, and returns the result.
M=[
[ 0, 2, 2 , 3 , 5],
[-3, -1, 1 , 5 , 9],
[ 1, -1, 1 , 4 , 7],
[ 1, -1, 1 , 0 , 2],
[ 1, -1, 1 , 0 , 3]
]
P=[2,0,1,4,3]
apply_permutations_matrix(M,P,5)
>>[
[ 1, -1, 1, 4, 7],
[ 0, 2, 2, 3, 5],
[-3, -1, 1, 5, 9],
[ 1, -1, 1, 0, 3],
[ 1, -1, 1, 0, 2]
]
'''
Mp=[ [M[ P[i] ][j]for j in range(N)] for i in range(N) ]
return Mp
def LUP(M,N,_tiny=1e-20):
U=[ [ M[i][j] for j in range(N)] for i in range(N) ]
L=[ [ 0 if i!=j else 1 for j in range(N)] for i in range(N) ]
#this is the "permutation vector". if it is e.g. [2 1 0 3] it means you make 0<->2
P=[ i for i in range(N) ]
for k in range(1,N):
for i in range(k,N):
#find the index of the maximum in column
_col=[Abs(U[_r][k-1]) for _r in range(k-1,N)]
#find the index of the maximum of _col
# notice that the length of _col is N-(k-1)
len_col=N-(k-1)
pivot=ind_max( _col ,len_col) + k - 1 #convert the index of _col (it has a length of len_col) to the index of a row of U
##################################################
#if you remove it, then you get a lot of infinities
#it has to do with the fact that if U[pivot][k-1] <_tiny , then U[k-1][k-1] will be a zero,
#L[i][k-1] explodes.
#You are allowed to skip this i, then, because if U[pivot][k-1] <_tiny , then all U[i][k-1] are small!
#Check that this is true by uncommenting print(_col)
if Abs(U[pivot][k-1]) < _tiny :
#print(_col)
break
###################################################
#if the maximum is not at k-1, swap!
if pivot != k-1 :
# Permute rows k-1 and pivot in U
index_swap(P,k-1,pivot)
tmpU=[U[k-1][_r] for _r in range(k-1,N)]
#print(U)
for _r in range(k-1,N):
U[k-1][_r]=U[pivot][_r]
#print(U)
for _r in range(k-1,N):
U[pivot][_r]=tmpU[_r-(k-1)]#again we have to convert the index of tmpU
#print(U)
#print("=========================")
tmpL=[L[k-1][_r] for _r in range(k-1)]
#print(L)
for _r in range(k-1):
L[k-1][_r]=L[pivot][_r]
#print(L)
for _r in range(k-1):
L[pivot][_r]=tmpL[_r]
#print(L)
#print("========================")
L[i][k-1]=U[i][k-1]/U[k-1][k-1]
for j in range(k-1,N):
U[i][j]=U[i][j]-L[i][k-1]*U[k-1][j]
return L,U,P