Implementation of an algorithm for finding the resolvent of a matrix using the adject matrix and characteristic polynomial

Аннотация

The article describes the implementation of an algorithm for calculating the resolvent of a matrix using the adjoint matrix and the characteristic polynomial of the matrix in Python. The graphs show the speed of the algorithm for various matrix dimensions.

Толық мәтін

Введение. При решении многих задач требуется вычислять резольвенту матрицы [1], не зная ее собственных значений. В частности, такая задача возникает при решении систем линейных алгебраических уравнений с малым параметром методом Ляпунова-Шмидта [2; 3]. В работе [1] приведен подход вычисления резольвенты матрицы с использованием присоединенной матрицы и характеристического многочлена матрицы. Для вычисления последних Д. К. Фаддеевым предложен метод одновременного вычисления коэффициентов

характеристического многочлена и присоединенной матрицы [1].

Изложим алгоритм вычисления резольвенты матрицы на основе подхода, изложенного в [1] и метода Д. К. Фаддеева одновременного вычисления коэффициентов характеристического многочлена и присоединенной матрицы.

Алгоритм вычисления резольвенты матрицы. Для вычисления резольвенты постоянной (m × m) −матрицы A

$R\left(\lambda \right)={\left[\lambda E-А\right]}^{-1}$ воспользуемся формулой [1, с. 112]

$R\left(\lambda \right)=\frac{1}{x_{{\min }^{\left(\lambda \right)}}}С\left(\lambda \right)$

где C(𝜆) – приведенная присоединенная матрица для матрицы [𝜆𝐸 − A] (см. [1, с. 99]);

𝜒_min (𝜆) – минимальный характеристический многочлен матрицы A.

Приведенную присоединенную матрицу C(𝜆) и минимальный характеристический многочлен матрицы A будем вычислять по формулам [1, с. 99-100]

$\begin{array}{l}C\left(\lambda \right)=\frac{1}{d\left(\lambda \right)}{B}_A\left(\lambda \right)\\ x_{min}\left(\lambda \right)=\frac{\Delta \left(\lambda \right)}{d\left(\lambda \right)}\end{array}$

где B_𝐴(𝜆) – присоединенная матрица для матрицы A [1, с. 92], 𝑑(𝜆) – наибольший общий делитель всех элементов матрицы B_𝐴(𝜆),

$\Delta \left(\lambda \right)={\lambda }^n-p_1{\lambda }^{n-1}-p_2{\lambda }^{n-2}-\mathrm{...}-p_{{}_{n-1}}\lambda -p_n$ – характеристический многочлен матрицы A.

Для нахождения присоединенной матрицы B_𝐴(𝜆) воспользуемся формулой [1, с. 94]

$B_A\left(\lambda \right)=B_0{\lambda }^{n-1}+B_1{\lambda }^{n-2}+B_2{\lambda }^{n-3}+\mathrm{...}+B_{n-1},$

где B₀– единичная (m× m) −матрица, а (m× m) −матрицы B_k, k= 1, … , n − 1, находятся методом Д. К. Фаддеева по формулам [1, с. 97]

$\begin{array}{l}{A}_k=A{B}_{k-1},\\ p_k=\frac{1}{k}Sp\text{(}{A}_k\text{),}\\ B_k=A_k-p_{k}E,\\ k=\mathrm{1,}\mathrm{...},m.\end{array}$

Здесь 𝐸₀ – единичная (m× m) −матрица.

Программная реализация алгоритм вычисления резольвенты матрицы. Приведем фрагменты алгоритма вычисления резольвенты матрицы A на языке Python с использованием пакета SymPy.

Фрагмент кода инициализации матрицы A с целыми случайными элементами

𝑎_ij, i, j = 1, … , m из отрезка [−10, 10]:

max_aij = 10

A = np.random.randint(-max_aij, max_aij + 1, (m, m))

Фрагмент кода реализации метода Д. К. Фаддеева одновременного вычисления характеристического многочлена ∆(𝜆) и присоединенной матрицы B_A(𝜆):

def fadeev_method(A):

m = len(A)

E = sympy.eye(m)

B = [0 for i in range(m + 1)]

B[0] = sympy.eye(m)

tmpA = [0 for i in range(m + 1)]

p = [0 for i in range(m +  1)]

for k in range(1, m + 1):

tmpA[k] = A * B[k -  1]

p[k] =  1 / k * matrix_trace(tmpA[k])

B[k] = tmpA[k] - p[k] * E

deltaLambda = lyamda ** m

for i in range(1, m + 1):

deltaLambda -= p[i] * lyamda ** (m - i)

                B_A = sympy.zeros(m)

for k in range(m):

B_A += B[k] * lyamda ** (m - 1 - k)

return B_A, deltaLambda

Фрагмент кода вычисления наибольшего общего делителя всех элементов присоединенной матрицы B_A(𝜆):

d = B_A[0, 0]

for i in range(m):

for j in range(m):

d = sympy.gcd(d, B_A[i, j])

Фрагмент кода вычисления приведенной присоединенной матрицы C(𝜆), минимального характеристического многочлена 𝜒_min(𝜆) и резольвенты R(𝜆) матрицы A:

C =  1 / d * B_A

x_min = deltaLambda.as_expr() / d

R =  1 / x_min * C

Вычислительный эксперимент. Вычисления проводились на ноутбуке с процессором 11th Gen Intel(R) Core(TM) i7-11800H @ 2.30GHz и операционной системой Windows 11 Pro.

На графиках представлены зависимости скорости работы алгоритма вычисления резольвенты матрицы A с целыми случайными элементами при различных размерностях матрицы.

Рисунок. Графики зависимостей скорости работы алгоритма от размерности матрицы A при условиях:
a) |𝑎_ij | ≤ 10, b) |𝑎_ij | ≤ 10⁵, c) |𝑎_ij | ≤ 10⁶, d) |𝑎_ij | ≤ 10⁷, i, j = 1, … , m

Из графиков видно, что при увеличении значений коэффициентов матрицы A по абсолютной величине скорость работы алгоритма увеличивается незначительно.

Для сравнительного анализа использовалась функция пакета Sympy нахождения обратной матрицы, зависящей от параметра. Вычисления показали, что при m = 30 и |𝑎_ij| ≤ 10 эта функция уже затрачивает порядка 10 мин., что значительно превышает скорость работы    разработанного алгоритма вычисления резольвенты матрицы с использованием приведенной присоединенной матрицы и минимального характеристического многочлена.