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 мин., что значительно превышает скорость работы    разработанного алгоритма вычисления резольвенты матрицы с использованием приведенной присоединенной матрицы и минимального характеристического многочлена.