Обращение матриц методом Шульца

[Schrodinger's cat]

Всем доброго здоровья и хорошего настроения.
Хочу запрограммировать алгоритм обращения матриц методом Шульца.
Как правильно инициализировать начальное приближение?
Если я делаю, как везде написано, исходя из половины спектрального радиуса (я так понимаю он равен максимальному собственную числу) умноженного либо на единичную матрицу в случае симметричной проблемы либо на $A \cdot A^T$ в случае произвольной проблемы то ничего не получается, ничего не сходится.
Да и можно ли как то обойтись без собственных чисел, так как сама по себе она довольно тяжеловесна?

[ipgmvq]

Schrodinger's cat в сообщении #1588670 писал(а):

ничего не сходится

А Ваша матрица точно обратима? Код и/или матрицу покажите?

[Schrodinger's cat]

ipgmvq
Точно обратима, так как она взята из примера отсюда:
Итерационный метод Шульца нахождения обратной матрицы
Пример 10.18
Собственно как получается первое приближение $U^{(0)}$ непонятно.

ПС: Код пока написан на MatCAD, что здешние обитатели не любят. Но если надо конечно выложу.
Код примитивен:
1 Создаем первое приближение.
2 Выполняем пункт 2 решения.
3 Выполняем пункт 4 решения.
4 Повторяем 2 и 3 итеративно до сходимости.
Если взять первое приближение из этого примера то алгоритм сходится. Проблема только в том как найти это первое приближение.

[zykov]

Schrodinger's cat в сообщении #1588670 писал(а):

Как правильно инициализировать начальное приближение?

Выбор начального приближения

[ipgmvq]

Schrodinger's cat в сообщении #1588670 писал(а):

обойтись без собственных чисел

Помимо такой
$U^{(0)} = \frac{A^{\dag}}{\left\lVert A \right\rVert_2^2}$
где $A^{\dag}$ эрмитово-сопряжённая матрица
можно ещё использовать такие
$U^{(0)} = \frac{A}{\left\lVert A \right\rVert_1^2}$
$U^{(0)} = \frac{A}{\left\lVert A \right\rVert_{\infty}^2}$
$U^{(0)} = \frac{A}{\left\lVert A \right\rVert_F^2}$
$U^{(0)} = \frac{A^{\intercal}}{n \left\lVert A \right\rVert_1 \left\lVert A \right\rVert_{\infty}}$

На Mathcad не пишу. Так это будет выглядеть на Python, где передаваемый номер метода (от 0 до 3 включительно) это как раз эти 4 альтернативных варианта выше:

Используется синтаксис Python

import numpy as np
from numpy.linalg import norm

матрица = np.matrix([[1, 2, 1], [0, 1, 0], [0, 2, 2]])

def найти_обратную(A0, error, X0=None, метод=0):
    методы = [lambda X: X/np.power(norm(X, ord=1), 2),
              lambda X: X/np.power(norm(X, ord=np.inf), 2),
              lambda X: X/np.power(norm(X, ord='fro'), 2),
              lambda X: X.T/(norm(X, ord=1)*norm(X, ord=np.inf)*A0.shape[0])]
    if X0 is None:
        X0 = методы[метод](A0)
        R0 = np.eye(A0.shape[0]) - A0 * X0
        if norm(R0.flatten(), ord=1) >= 1.0:
            return 'Метод не позволяет сформировать начальную матрицу, удовлетворяющую условию'
        else:
            return найти_обратную(A0, error=error, X0=X0, метод=метод)
    else:
        R0 = np.eye(A0.shape[0]) - A0 * X0
        if norm(R0, ord='fro') < error:
            return X0
        else:
            X0 = X0 * (np.eye(A0.shape[0]) + R0)
            return найти_обратную(A0, error=error, X0=X0, метод=метод)
        
обратная = найти_обратную(матрица, error=0.001 , метод=3)
print("Обратная матрица:\n\n", обратная, "\n")
print("Произведение с обратной матрицей:\n\n", обратная * матрица, "\n")
 

[Schrodinger's cat]

Все работает, всем спасибо.