finding the iterative Schultz method

Question

Please help me solve problems with finding the iterative Schultz method. Overflow occurs, but I don't know how to fix it

Gives an error message RuntimeWarning: overflow encountered in 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)

I tried to change float64 to float32 and set fewer iterations and change the error, but it didn't help

at the exit or [[-inf -inf] [-inf -inf]] either an error

Answer 1

I think the problem is not the implementation.

It should be the init of X.

For case A=I, X=I, 2X-XAX = I. This is the simplest case that you code works.

For your test case....

You can see it diverges... quite fast tbh (faster than I expected).

Therefore, I did an effortless change. Which is start with a smaller X.

For something a little better... I suggest you init X as I/det(A)... Which is np.eye(n)/np.linalg.det(A); Or even... as np.eye(n)/np.sum(A)