使用 Neumann 级数的矩阵求逆给出有趣的损失函数
Matrix inversion using Neumann Series giving funny loss function
根据 (steward,1998)。可逆矩阵 A 可以用公式 A^{-1} = \sum^{inf}_{n=0 来近似} (I- A)^{n}
我尝试实现一种算法来近似简单矩阵的逆,损失函数显示出有趣的结果。请看下面的代码。可以找到有关 Neumann 系列的更多信息 here and here
这是我的代码。
A = np.array([[1,0,2],[3,1,-2],[-5,-1,9]])
class Neumann_inversion():
def __init__(self,A,rank):
self.A = A
self.rank = rank
self.eye = np.eye(len(A))
self.loss = []
self.loss2 =[]
self.A_hat = np.zeros((3,3),dtype = float)
#self.loss.append(np.linalg.norm(np.linalg.inv(self.A)-self.A_hat))
def approximate(self):
# self.A_hat = None
n = 0
L = (self.eye-self.A)
while n < self.rank:
self.A_hat += np.linalg.matrix_power(L,n)
loss = np.linalg.norm(np.linalg.inv(self.A) - self.A_hat)
self.loss.append(loss)
n+= 1
plt.plot(self.loss)
plt.ylabel('Loss')
plt.xlabel('rank')
# ax.axis('scaled')
return
Matrix = Neumann_inversion(A,200)
Matrix.approximate()
仅当 $A^n$ 随着 $n$ 的增加趋于零时,该公式才有效。所以你的矩阵必须满足
np.all(np.abs(np.linalg.eigvals(A)) < 1)
尝试
Neumann_inversion(A/10, 200).approximate()
你可以认真对待损失:)
公式的由来有点关系
(1-x) * (1 + x + x^2 + ... x^n) = (1 - x^(n+1))
当且仅当矩阵的所有特征值的大小都小于 1 时,x^(n+1) 项将接近于零,因此总和将近似为 (1-x ).
根据 (steward,1998)。可逆矩阵 A 可以用公式 A^{-1} = \sum^{inf}_{n=0 来近似} (I- A)^{n}
我尝试实现一种算法来近似简单矩阵的逆,损失函数显示出有趣的结果。请看下面的代码。可以找到有关 Neumann 系列的更多信息 here and here
这是我的代码。
A = np.array([[1,0,2],[3,1,-2],[-5,-1,9]])
class Neumann_inversion():
def __init__(self,A,rank):
self.A = A
self.rank = rank
self.eye = np.eye(len(A))
self.loss = []
self.loss2 =[]
self.A_hat = np.zeros((3,3),dtype = float)
#self.loss.append(np.linalg.norm(np.linalg.inv(self.A)-self.A_hat))
def approximate(self):
# self.A_hat = None
n = 0
L = (self.eye-self.A)
while n < self.rank:
self.A_hat += np.linalg.matrix_power(L,n)
loss = np.linalg.norm(np.linalg.inv(self.A) - self.A_hat)
self.loss.append(loss)
n+= 1
plt.plot(self.loss)
plt.ylabel('Loss')
plt.xlabel('rank')
# ax.axis('scaled')
return
Matrix = Neumann_inversion(A,200)
Matrix.approximate()
仅当 $A^n$ 随着 $n$ 的增加趋于零时,该公式才有效。所以你的矩阵必须满足
np.all(np.abs(np.linalg.eigvals(A)) < 1)
尝试
Neumann_inversion(A/10, 200).approximate()
你可以认真对待损失:)
公式的由来有点关系
(1-x) * (1 + x + x^2 + ... x^n) = (1 - x^(n+1))
当且仅当矩阵的所有特征值的大小都小于 1 时,x^(n+1) 项将接近于零,因此总和将近似为 (1-x ).