寻找迭代 Schultz 方法

Question

请帮我解决寻找迭代 Schultz 方法的问题。发生溢出，但我不知道如何解决

给出错误信息 RuntimeWarning：matmul 遇到溢出 B=2*X-X@A@X

import numpy as np

def schultz_inverse(A, tolerance=1e-6, max_iterations=100):
    A = A.astype(float) # matrix data type conversion
    n = A.shape[0]
    X = np.eye(n) # initial approximation for X

    for i in range(max_iterations):
        B = 2 * X - X @ A @ X
        if np.allclose(B, X, rtol=tolerance):
            return B
        X = B

    raise Exception("The method did not converge after the maximum number of iterations")

# Usage example
A = np.array([['4', '1'], ['2', '3']])
A_inv = schultz_inverse(A)
print(A_inv)

我试过把float64改成float32，设置更少的迭代次数，改了错误，但是没用

在出口或 [[-信息-信息] [-inf -inf]] 要么是错误

Answer 1

我认为问题不在于实施。

应该是X的init.

对于案例 A=I, X=I, 2X-XAX = I。这是您编写代码的最简单案例。

对于您的测试用例....

你可以看到它发散......相当快（比我预期的快）。

因此，我做了一个毫不费力的改变。从较小的 X 开始。

为了更好的东西......我建议你 init X as I/det(A)... 这是

np.eye(n)/np.linalg.det(A)

; 甚至...作为

np.eye(n)/np.sum(A)