Матрицы и векторы: вычисление обратной матрицы / Хабр

Вычисление обратной матрицы квадратной матрицы высокого порядка, а именно, вычисление алгебраических дополнений и определителя матрицы займёт большое количество машинных ресурсов. Поэтому обращение матрицы при помощи алгебраических дополнений нерационально. В этой статье для вычисления обратной матрицы были использованы нижеследующие методы:

Решение системы AX = E, где A, X, E - квадратные матрицы порядка n, X - обратная A матрица, E - единичная матрица, $E_{ij} = \begin{cases}1 & i = j\\0 & i \neq j\end{cases}$
Использование LU разложения матрицы A = LU [1,2,3]
Обращение матриц было реализовано в виде шаблонов функций шаблона класса MATRIX<T>. В данной статье обсуждается решение вышеуказанными методами обращения квадратной матрицы порядка n и полученные результаты.
Решение системы AX = E
Решение матричного уравнения AX = E для матрицы порядка n A сводится к решению n систем линейных уравнений:

$\begin{matrix} A\cdot \bar{X_{j}} = \bar{e_{j}},\quad j = 1...n \\ \bar{X_{j}} = \begin{bmatrix}X_{1j} \\X_{2j} \\ ... \\ X_{nj} \end{bmatrix},\quad j-ый\;столбец \;матрицы\;X_{ij} \\ \bar{e_{j}},\quad j-ый\;столбец\;матрицы\; E \end{matrix}$

Код реализации решения приводится ниже:

/// <summary>
/// Формирование матриц alpha и gamma в компактной схеме исключения
/// Матрица alpha верхняя треугольная, gamma - нижняя треугольная
/// </summary>
/// <param name="alpha">формируемые матрицы</param>
/// <returns>значение определителя матрицы</returns>
template <typename T>
T MATRIX<T>::FormMatrixCompactScheme(MATRIX& alpha)
{
	if (this->m_columns != this->m_rows) return NAN;
	int n = this->m_rows;

	for (int i = 0; i < n; i++)
	{
		*(alpha.m_data + i * n) = *(this->m_data + i * n);
		if (i > 0) *(alpha.m_data + i) = *(this->m_data + i) / *this->m_data;
	}

	T sum = 0;
	int k = 1, i = 1;
	while (i < n)
	{
		if (k >= n)
		{
			k = 1; i++;
			if (i >= n) break;
		}
		if (i >= k)
		{
			sum = 0.0;
			for (int j = 0; j <= k - 1; j++)
			{
				sum += *(alpha.m_data + i * n + j) * *(alpha.m_data + j * n + k);
			}
			*(alpha.m_data + i * n + k) = *(this->m_data + i * n + k) - sum;
		}
		else
		{
			sum = 0.0;
			for (int j = 0; j <= i - 1; j++)
			{
				sum += *(alpha.m_data + i * n + j) * *(alpha.m_data + j * n + k);
			}
			if (abs(*(alpha.m_data + i * n + i)) < 1.0e-18)
			{
				return 0.0;
			}
			*(alpha.m_data + i * n + k) = (*(this->m_data + i * n + k) - sum) / *(alpha.m_data + i * n + i);
		}
		k++;
	}

	// вычисление определителя
	T det = 1;
	for (int i = 0; i < n; i++)
		det *= *(alpha.m_data + i * n + i);
	return det;
}

/// <summary>
/// Решение системы линейных алгебраических уравнений A*x = b с 
/// верхней треугольной матрицей
/// </summary>
/// <param name="A">верхняя треугольная матрица системы</param>
/// <param name="b">вектор правой части системы</param>
/// <param name="x">решение системы</param>
template <typename T>
void TriangleSolve(MATRIX<T>& A, VECTOR<T>& b, VECTOR<T>& x)
{
	int n = b.m_size;
	VECTOR<T> beta(n);
	T sum = 0;
	*beta.m_data = *b.m_data / *A.m_data;
	for (int i = 1; i < n; i++)
	{
		sum = 0.0;
		for (int j = 0; j <= i - 1; j++)
			sum += *(A.m_data + i * n + j) * *(beta.m_data + j);
		*(beta.m_data + i) = (*(b.m_data + i) - sum) / *(A.m_data + i * n + i);
	}

	// решение системы уравнений с труегольной матрицей		   
	*(x.m_data + n - 1) = *(beta.m_data + n - 1);
	for (int i = n - 2; i >= 0; i--)
	{
		sum = 0.0;
		for (int j = n - 1; j > i; j--)
			sum += *(A.m_data + i * n + j) * *(x.m_data + j);
		*(x.m_data + i) = *(beta.m_data + i) - sum;
	}

}
/// <summary>
/// Копировать вектор в  j-ю колонку матрицы
/// </summary>
/// <param name="v">вектор</param>
/// <param name="j">номер колонки</param>
template <typename T>
void MATRIX<T>::CopyColumn(VECTOR<T>& v, int j)
{
	if (m_columns < j)
	{
		return;
	}
	if (m_rows != v.size())
	{
		return;
	}
	for (int i = 0; i < m_rows; i++)
		*(m_data + i * m_columns + j) = v[i];
}

/// <summary>
/// Вычисляет первую норму матрицы
/// </summary>
/// <returns>значение первой нормы матрицы</returns>
template <typename T>
T MATRIX<T>::normI()
{
	T norm = 0;
	for (int i = 0; i < m_rows; i++)
	{
		T sum = 0;
		for (int j = 0; j < m_columns; j++)
			sum += abs(*(m_data + i * m_columns + j));
		if (sum > norm) norm = sum;
	}
	return norm;
}

/// <summary>
/// Обращение марицы при помощи решения систем уравнений A*X = E,
/// столбцы матрицы X - столбцы обратной к A матрицы,
/// E - единичная матрица
/// </summary>
/// <returns>обратную матрицу</returns>
template <typename T>
MATRIX<T> MATRIX<T>::Invert()
{
	MATRIX A(m_rows, m_columns);
	if (m_rows != m_columns)
		return A;

	MATRIX alpha(m_rows, m_columns);
	T det = FormMatrixCompactScheme(alpha);

	if (abs(det) < 1.0e-36)
		return A;

	for (int i = 0; i < m_rows; i++)
	{
		VECTOR<T> e(m_rows), x(m_rows);
		e[i] = 1;
		TriangleSolve(alpha, e, x);
		A.CopyColumn(x, i);
	}

	// итерационное уточнение обратной матрицы
	if (m_rows >= MIN_SIZE_INVERSION)
	{
		MATRIX E(m_rows, m_columns);
		for (int i = 0; i < m_rows; i++)
			*(E.m_data + i * m_columns + i) = 1;
		int iter = 0;
		T norm;
		do       
        {
			MATRIX R(E - *this * A);
			A = A * (E + R);
			norm = R.normI();
			if (++iter > MAX_ITER_NUMBER) break;
		} while (norm < 1.0e-17);
	}

	return A;
}

Перегрузка операторов *, =, +, - реализована в шаблоне класса MATRIX<T>. Решение системы уравнений осуществлялось с использованием компактной схемы исключения [4].

Обращение матрицы с использованием LU разложения матрицы

Разложение матрицы A = LU на произведение нижней треугольной L и верхней треугольной U может быть использовано не только для решения системы линейных алгебраических уравнений, но и для обращения матрицы. Описание метода обращения матрицы при помощи LU разложения приведён в источниках [2,3] и заключается в следующем:

$\begin{matrix} A = L\cdot U \Leftrightarrow A^{-1} = U^{-1}\cdot L^{-1} \\ U\cdot A^{-1} = L,\;A^{-1}\cdot L=U^{-1} \end{matrix}$

Обозначая искомые элементы матрицы $A^{-1}$ через $x_{ij}$ и учитывая, что треугольные матрицы при обращении сохраняют свою структуру, последние соотношения можно записать в виде:

$\begin{matrix} \begin{bmatrix}u_{11} & u_{12} & ... & u_{1n} \\ 0 & u_{22} & ... & u_{2n} \\ ... & .... & ... & ... \\ 0 & 0 & 0 & u_{nn} \end{bmatrix}\cdot \begin{bmatrix}x_{11} & x_{12} & ... & x_{1n} \\ x_{21} & x_{22} & ... & x_{2n} \\ ... & .... & ... & ... \\ x_{n1} & x_{n2} & ... & x_{nn} \end{bmatrix} = \begin{bmatrix}1 & 0 & ... & 0 \\ * & 1 & ... & 0 \\ ... & .... & ... & ... \\ * & * & ... & 1 \end{bmatrix} \\ \begin{bmatrix}x_{11} & x_{12} & ... & x_{1n} \\ x_{21} & x_{22} & ... & x_{2n} \\ ... & .... & ... & ... \\ x_{n1} & x_{n2} & ... & x_{nn} \end{bmatrix}\cdot \begin{bmatrix}1 & 0 & ... & 0 \\ l_{21} & 1 & ... & 0 \\ ... & .... & ... & ... \\ l_{n1} & l_{n2} & ... & 1 \end{bmatrix} = \begin{bmatrix}* & * & ... & * \\ 0 & * & ... & * \\ ... & .... & ... & ... \\ 0 & 0 & ... & * \end{bmatrix} \end{matrix}$

где символом '*' обозначены элементы, знание значений которых для нахождения обратной матрицы не требуется.

Из этих соотношений можно получить выражения для коэффициентов обратной матрицы $x_{ij}$ :

$\begin{matrix} x_{ij} = -\frac{1}{u_{ii}}\sum\limits_{k=i+1}^n u_{ik}x_{kj}\;(i<j) \\ x_{ji}= \frac{1}{u_{jj}}\sum\limits_{k=j+1}^n (1 - u_{jk}x_{kj})\;(i=j) \\ x_{ij} = -\sum\limits_{k=j+1}^n x_{ik}l_{kj}\;(i > j) \end{matrix}$

Код реализации решения приводится ниже:

/// <summary>
/// Декомпозиция матрицы A = LU
/// L - нижняя треугольная матрица, на главной диагонали которой расположены единицы
/// U - верхняя треугольная матрица
/// </summary>
/// <param name="alfa">Матрицы L и U, ниже главной диагонали которой 
/// расположены внедиагональные элементы L, 
/// на главной диагонали и выше расположены элементы матрицы U</param>
/// <returns>Определитель матрицы</returns>
template <typename T>
T MATRIX<T>::LUDecomposition(MATRIX& alfa)
{
	T det = 0;
	if (m_rows != m_columns) return det;

	// декомпозиция матрицы
	for (int i = 0; i < m_rows; i++)
		for (int j = 0; j < m_columns; j++)
		{
			T val = 0;
			if (i <= j)
			{
				for (int k = 0; k <= i - 1; k++)
					val += *(alfa.m_data + i * alfa.m_columns + k) * *(alfa.m_data + k * alfa.m_columns + j);
				*(alfa.m_data + i * alfa.m_columns + j) = *(m_data + i * m_columns + j) - val;
			}
			else
			{
				for (int k = 0; k <= j - 1; k++)
					val += *(alfa.m_data + i * alfa.m_columns + k) * *(alfa.m_data + k * alfa.m_columns + j);
				if (abs(*(alfa.m_data + j * alfa.m_columns + j)) < 1.0e-36) return 0;
				*(alfa.m_data + i * alfa.m_columns + j) = (*(m_data + i * m_columns + j) - val) / *(alfa.m_data + j * alfa.m_columns + j);

			}
		}
	det = 1;
	for(int i = 0; i < alfa.m_rows; i++)
		det *= *(alfa.m_data + i * alfa.m_columns + i);
	return det;
}

/// <summary>
/// Вычисляет первую норму матрицы
/// </summary>
/// <returns>значение первой нормы матрицы</returns>
template <typename T>
T MATRIX<T>::normI()
{
	T norm = 0;
	for (int i = 0; i < m_rows; i++)
	{
		T sum = 0;
		for (int j = 0; j < m_columns; j++)
			sum += abs(*(m_data + i * m_columns + j));
		if (sum > norm) norm = sum;
	}
	return norm;
}
/// <summary>
/// Обращение матриы с пощью LU разложения
/// </summary>
/// <returns>обратную матрицу</returns>
template <typename T>
MATRIX <T> MATRIX <T>::InvertLT()
{
	MATRIX D(m_rows, m_columns), alpha(m_rows, m_columns);
	if (m_rows != m_columns) // матрица должна быть квадратной
		return D;
	T det = LUDecomposition(alpha); // декомпозиция A = LU
	if (abs(det) < 1.0e-36)
		return D;
	// формирование элементов обратной матрицы
	int n = m_rows;
	for (int i = n-1; i >= 0; i--)
		for(int j = n-1;j >= 0; j--)
		{
			T sum = 0;
			if (i < j)
			{
				for (int k = i + 1; k < n; k++)
					sum -= *(alpha.m_data + i * alpha.m_columns + k) * *(D.m_data + k * D.m_columns + j);
				*(D.m_data + i * D.m_columns + j) =  sum / *(alpha.m_data + i * alpha.m_columns + i);

			}
			else if (i == j)
			{
				for (int k = j + 1; k < n; k++)
					sum += *(alpha.m_data + j * alpha.m_columns + k) * *(D.m_data + k * D.m_columns + j);
				*(D.m_data + j * D.m_columns + i) = (1-sum) / *(alpha.m_data + j * alpha.m_columns + j);

			}
			else
			{
				for (int k = j + 1; k < n; k++)
					sum -= *(alpha.m_data + k * alpha.m_columns + j) * *(D.m_data + i * D.m_columns + k);
				*(D.m_data + i * D.m_columns + j) = sum;

			}
		}
		// итерационное уточнение обратной матрицы
		if (m_rows >= MIN_SIZE_INVERSION)
		{
			MATRIX E(m_rows, m_columns);
			for (int i = 0; i < m_rows; i++)
				*(E.m_data + i * m_columns + i) = 1;
			int iter = 0;
			T norm;
			do
			{
				MATRIX R(E - *this * D);
				D = D * (E + R);
				norm = R.normI();
				if (++iter > MAX_ITER_NUMBER) break;
			} while (norm < 1.0e-17);
		}

	return D;
}

Результаты расчётов

Расчёты по обращению квадратной матрицы порядка n проводились обоими методами для разных значений n. Исходная матрица A порядка n заполнялась при помощи генератора случайных чисел в интервале от -1000 до 1000. В качестве меры погрешности вычислений была использована I-норма

$\parallel A\parallel_{I} = \begin{matrix}max & \sum\limits_{j=1}^n \mid a_{ij}\mid \\ 0\leq i \leq n\end{matrix}$

матрицы $E - A^{-1}A$ , где E - единичная матрица.

Расчёты проводились без итерационного уточнения и с итерационным уточнением обратной матрицы.

Итерационное уточнение осуществлялось следующим образом [2]. Пусть $D_{0}$ первое приближение матрицы $A^{-1}$ . Если при этом $\parallel{R_{0}}\parallel_{I}\leq k < 1,\quad где\; R_{0} = E-AD_{0}$ , то можно построить итерационный процесс:

$\begin{matrix} D_{1} = D_{0}(E-R_{0}),\quad R_{1} = E- AD_{1} \\ D_{2} = D_{1}(E-R_{1}),\quad R_{2} = E- AD_{2} \\ ...................................\\ D_{m} = D_{m-1}(E-R_{m-1}),\quad R_{m} = E- AD_{m} \end{matrix}$

Непосредственной проверкой можно убедиться в том, что [2]:

$R_{m}=R_{m-1}^{2}=R_{m-2}^{4}=...=R_{0}^{2^{m}}$

Откуда далее следует оценка погрешности [2]:

$\parallel{D_{m} - A^{-1}}\parallel_{I}=\parallel{D_{0}}\parallel_{I}\cdot \frac{k^{2^{m}}}{1-k}$

Расчёт без итерационного уточнения

Результаты расчётов без итерационного уточнения матрицы $A^{-1}$ приведены в нижеследующей таблице:

Порядок матрицы A	Решение уравнения AX=E		Декомпозиция A=LU
Порядок матрицы A	Время, сек	Погрешность	Время, сек	Погрешность
50	0,0007115	$1,82559\cdot 10^{-11}$	0,0003418	$1,20276\cdot 10^{-9}$
100	0,0030217	$1,38811\cdot 10^{-10}$	0,0024037	$2,69221\cdot 10^{-8}$
200	0,0238938	$8,36936\cdot 10^{-10}$	0,0193017	$2,20405\cdot 10^{-6}$
500	0,349852	$5,84951\cdot 10^{-9}$	0,288407	$6,42369\cdot 10^{-6}$
1000	2,83838	$4,17686\cdot 10^{-7}$	2,36379	0,0509241
2000	27,9798	$2,53471\cdot 10^{-7}$	33,5439	0,00136141

Расчёт с итерационным уточнением

По приведённым результатам расчёта без итерационного уточнения было сделано заключение о нецелесообразности итерационного уточнения обращения квадратной матрицы порядка ниже 200. Результаты расчётов с итерационным уточнением матрицы $A^{-1}$ приведены в нижеследующей таблице:

Порядок матрицы A	Решение уравнения AX=E		Декомпозиция A=LU
Порядок матрицы A	Время, сек	Погрешность	Время, сек	Погрешность
200	0,0981653	$5,85991\cdot 10^{-12}$	0,0836408	$5,76079\cdot 10^{-12}$
500	1,44216	$1,13022\cdot 10^{-10}$	1,31034	$1,13075\cdot 10^{-10}$
1000	11,4696	$1,70632\cdot 10^{-10}$	10,9775	$2,17188\cdot 10^{-10}$
2000	106,481	$6,39272\cdot 10^{-10}$	109,584	$6,52071\cdot 10^{-10}$

Выводы

Обращение матрицы без итерационного уточнения метод LU декомпозиции показывает лучшие результаты по времени, а метод решения системы AX = E демонстрирует лучшие результаты по точности. При этом, с ростом порядка квадратной матрицы разница в точности расчёта методом LU декомпозиции и методом решения системы AX = E возрастает на несколько порядков, а при обращении матрицы порядка 2000 метод решения системы AX = E демонстрирует лучшие результаты как по времени, так и по точности.

Сравнительный анализ обращения матрицы с итерационным уточнением методом LU декомпозиции и методом решения системы AX = E показывает, что разница как по точности расчёта, так по времени несущественна, а итерационное уточнение вычисления обратной матрицы даёт хорошие результаты по точности методом LU в сравнении с решением без итерационного уточнения.

Список источников

Википедия. LU - разложение/ URL: https://ru.wikipedia.org/wiki/LU-разложение (дата посещения 18.12.2022)
Фаддеев, Д.К. Вычислительные методы линейной алгебры: Учебник/ Фаддеев Д.К, Фаддеева В.Н. — 4‑е изд.— СПб : Лань, 2009. — 736 с.
Численные методы на Python: LU разложение/ URL: https://orion1401.gitbooks.io/numerical_analysys_python/content/lu_razlozhenie.html (дата посещения 22.12.2025)
Корн, Г. Справочник по математике для научных работников и инженеров : определения, теоремы, формулы / Г. Корн, Т. Корн ; пер. с англ. под общ. ред. И. Г. Арамановича. — 3‑е изд. — Москва : Наука, 1973. — 832 с.

Матрицы и векторы: вычисление обратной матрицы

Решение системы AX = E