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

Подавляющая часть прикладных задач математической физики решается численными методами, неотъемлемой частью которых является решение системы линейных алгебраических уравнений (СЛАУ) $A\cdot x = b$ , где A - матрица коэффициентов системы, b - вектор правой части, x - вектор решений.

Поэтому, для получения наиболее точного решения прикладной задачи необходимо научиться наиболее точно решать СЛАУ. В данной статье будет рассмотрено решение СЛАУ несколькими методами:

методом Гаусса;
LU декомпозицией исходной матрицы;
QR декомпозицией исходной матрицы;
компактной схемой исключения.

Шаблоны классов MATRIX<T> (Матрица) и VECTOR<T>

Для упрощения работы с матрицами и векторами были созданы шаблоны классов MATRIX<T> (Матрица) и VECTOR<T> (Вектор). $n\cdot m$ элементов матрицы типа T, где n - число строк, а m - число столбцов, размещаются в одномерном массиве данных по строкам. Элементы n-мерного вектора располагаются в одномерном массиве типа T. Шаблоны классов MATRIX<T> и VECTOR<T> содержат перегруженные операции =, *, +, -. Шаблоны классов MATRIX<T> и VECTOR<T> представлены в файле MatrixVectorTemplate.hpp проекта MatrixVector. Функции, осуществляющие решение СЛАУ, представлены дружественными функциями классов MATRIX<T> и VECTOR<T>

	template <typename T> friend bool Gauss(const MATRIX<T>& a, const VECTOR<T>& b, VECTOR<T>& x);
	template <typename T> friend void CompactSchemeSolve(MATRIX<T>& A, VECTOR<T>& b, VECTOR<T>& x);
	template <typename T> friend void QRDecompositionSolve(MATRIX<T>& A, VECTOR<T>& b, VECTOR<T>& x);
	template <typename T> friend void LLTDecompositionSolve(MATRIX<T>& A, VECTOR<T>& b, VECTOR<T>& x);
	template <typename T> friend void LUDecompositionSolve(MATRIX<T>& A, VECTOR<T>& b, VECTOR<T>& x);

Решение СЛАУ

Проверка решения СЛАУ различными методами осуществлялась на примере матрицы коэффициентов СЛАУ A и вектора правой части b.

$A = \begin{bmatrix}1789.0 & 17.0 & 45.0 & 1205.0 & 13.0 & 23.09 \\16.0 & 22.0 & 48.0 & 34.0 & 56088.0001 & 3456.0 \\99.0 & 21.0 & 14.0 & 7.0 & 2.0 & 782.0 \\112.0 & 117.0 & 29.0 & 456.0234 & 22.0 & 435.0 \\0.345 & 0.0089 & 1.004 & 2.67 & 16.42 & 0.009 \\0.009 & 0.003 & 240.34 & 0.998 & 0.87 & 0.22 \end{bmatrix} b = \begin{bmatrix} 19.0 \\77.0 \\66.0 \\11.0 \\22.0 \\234.0\end{bmatrix}$

Производилась оценка точности решения в виде евклидовой нормы невязки $\parallel A\cdot x - b\parallel$ , а также вычислялось время решения.

Метод Гаусса

Метод Гаусса является, пожалуй, самым известным методом решения СЛАУ. Он хорошо описан в методической, учебной и научной литературе. Код представлен в виде дружественной классам MATRIX<T> и VECTOR<T> функции:

/// <summary>
/// Решение системы линейных алгебраических уравнений (СЛАУ) методом Гаусса с вычислением определителя матрицы системы 
/// </summary>
/// <param name="a">матрица коэффициентов СЛАУ</param>
/// <param name="b">вектор правой части СЛАУ</param>
/// <param name="x">решение СЛАУ</param>
/// <returns></returns>
template <typename T>
bool Gauss(const MATRIX<T>& a, const VECTOR<T>& b, VECTOR<T>& x)
{
	int i, k, m;
	T amm, aim;

	// матрица должна быть квадратной и размерность вектора должна совпадать 
	// с размерностью матрицы
	if (a.m_columns != b.m_size || a.m_columns != a.m_rows) return false;
	int size = a.m_rows;

	// сведение исходной системы к системе с верхней треугольной матрицей
	MATRIX<T> alf(a.m_rows, a.m_columns);
	VECTOR<T> bet(b.m_size);
	alf = a;
	bet = b;
	for (m = 0; m <= size - 2; m++)
	{
		amm = *(alf.m_data + m * size + m);
		for (k = m; k <= size - 1; k++)
			*(alf.m_data + m * size + k) /= amm;
		*(bet.m_data + m) /= amm;
		for (i = m + 1; i <= size - 1; i++)
		{
			aim = *(alf.m_data + i * size + m);
			for (k = m; k <= size - 1; k++)
				*(alf.m_data + i * size + k) -= *(alf.m_data + m * size + k) * aim;
			*(bet.m_data + i) -= *(bet.m_data + m) * aim;
		}//end i 
	}//end m 

	// нахождение решения СЛАУ с верхней треугольной матрицей
	*(x.m_data + size - 1) = *(bet.m_data + size - 1) / *(alf.m_data + size * (size - 1) + size - 1);
	for (i = size - 2; i >= 0; i--)
	{
		*(x.m_data + i) = *(bet.m_data + i);
		for (k = i + 1; k < size; k++)
			*(x.m_data + i) -= *(alf.m_data + i * size + k) * *(x.m_data + k);
	}//end i
	return true;

}

LU декомпозиция

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

$\begin{matrix} L\cdot y = b \\ U\cdot x = y \end{matrix}$

Код представлен в виде дружественной классам MATRIX<T> и VECTOR<T> функции:

/// <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>
/// Решение системы линейных алгебраических уравнений с применением метода LU декомпозиции A = LU
/// L - нижняя тругольная матрица
/// U - верхняя треугольная матрица 
/// </summary>
/// <param name="A">матрица</param>
/// <param name="b">вектор правой части системы</param>
/// <param name="x">решение системы</param>

template <typename T>
void LUDecompositionSolve(MATRIX<T>& A, VECTOR<T>& b, VECTOR<T>& x)
{
	if (A.m_rows != A.m_columns || A.m_rows != b.m_size) return;

	MATRIX<T> alfa(A.m_rows, A.m_columns);
	int n = A.m_rows;

	T det = A.LUDecomposition(alfa);
	if (abs(det) > 1.0e-36)
	{
		// решение СЛАУ
        // Решение системы Ly = b
		VECTOR<T> y(n);
		for (int k = 0; k < n; k++)
		{
			T val = 0;
			for (int j = 0; j <= k - 1; j++)
				val += *(alfa.m_data + k * alfa.m_columns + j) * *(y.m_data + j);
			*(y.m_data + k) = *(b.m_data + k) - val;
		}

        // Решение системы Ux = y
		for (int k = n - 1; k >= 0; k--)
		{
			T sum = 0;
			for (int j = n - 1; j > k; j--)
				sum += *(alfa.m_data + k * alfa.m_columns + j) * *(x.m_data + j);
			*(x.m_data + k) = (*(y.m_data + k) - sum) / *(alfa.m_data + k * alfa.m_columns + k);
		}
	}
}

Компактная схема исключения

Компактная схема исключения описана в справочнике [1]. Компактная схема исключения заключается в сведении решения СЛАУ к решению СЛАУ с верхней треугольной матрицей. Код представлен в виде дружественной классам MATRIX<T> и VECTOR<T> функции:

/// <summary>
/// Формирование матрицы alpha. Верхняя треугольная матрица
/// раполагается на главной диагонали и выше главной диагонали.
/// Значение определителя - произведпение диагональных элементов
/// </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>
/// Решение системы линейных уравнений компактной схемой исключения
/// </summary>
/// <param name="A">матрица системы уравнений</param>
/// <param name="b">вектор правой части системы уравнений</param>
/// <param name="x">решение системы уравнений</param>
template <typename T>
void CompactSchemeSolve(MATRIX<T>& A, VECTOR<T>& b, VECTOR<T>& x)
{
	if (A.m_rows != A.m_columns) return;
	MATRIX<T> alpha(A.m_rows, A.m_columns);
	T det = A.FormMatrixCompactScheme(alpha);
	// решение только для неособенной матрицы
	if (det != 0.0)
		TriangleSolve(alpha, b, x);
}

QR декомпозиция

Суть метода заключается в том, что исходная матрица представляется в виде произведения ортогональной матрицы Q и верхней треугольной матрицы R, A = QR [3,4]. Декомпозиция матрицы осуществляется применением процедуры ортогонализации Грама-Шмидта к столбцам матрицы СЛАУ A. Таким образом, исходная система уравнений сводится к системе уравнений с верхней треугольной матрицей:

$\begin{matrix} R\cdot x = Q^{T}\cdot b \\ где\: Q^{T} - транспонированная\: к\: Q\: матрица \end{matrix}$

Код представлен в виде дружественной классам MATRIX<T> и VECTOR<T> функции:

/// <summary>
/// Получение транспонированной матрицы 
/// </summary>
/// <returns>транспонированную матрицу</returns>
template <typename T>
MATRIX<T> MATRIX<T>::Transpose()
{
	MATRIX<T> tr(m_columns, m_rows);
	for (int i = 0; i < m_columns; i++)
		for (int j = 0; j < m_rows; j++)
			*(tr.m_data + i * tr.m_columns + j) = *(m_data + j * m_columns + i);
	return tr;
}

/// <summary>
/// QR разложение квадратной матрицы
/// </summary>
/// <param name="Q">матрица Q, ортогональная матрица</param>
/// <param name="R">матрица R, верхняя треугольная матрица</param>
template <typename T>
bool MATRIX<T>::QRDecomposition(MATRIX& Q, MATRIX& R)
{
	if (m_columns != m_rows) return false;
	int n = m_rows;
	T sum = 0;

	for (int j = 0; j < n; j++)
	{
		// q(j) = a(j)
		for (int k = 0; k < n; k++)
			*(Q.m_data + k * n + j) = *(m_data + k * n + j);

		for (int i = 0; i <= j - 1; i++)
		{
			// rij = q(i)^T*a(j)
			sum = 0.0;
			for (int k = 0; k < n; k++)
				sum += *(Q.m_data + k * n + i) * *(m_data + k * n + j);
			*(R.m_data + i * n + j) = sum;

			// r(i,j)*q(i)
			for (int k = 0; k < n; k++)
				*(Q.m_data + k * n + j) -= *(R.m_data + i * n + j) * *(Q.m_data + k * n + i);
		}

		// r(j,j) = || q(j)) ||2
		sum = 0.0;
		for (int k = 0; k < n; k++)
			sum += *(Q.m_data + k * n + j) * *(Q.m_data + k * n + j);
		*(R.m_data + j * n + j) = sqrt(sum);

		if (*(R.m_data + j * n + j) == 0.0)
		{
			return false;
		}
		for (int k = 0; k < n; k++)
			*(Q.m_data + k * n + j) /= *(R.m_data + j * n + j);

	}

	return true;
}

/// <summary>
/// Решение СЛАУ с применением QR декомпозиции матрицы системы A
/// </summary>
/// <param name="A">матрица системы линейных уравнений</param>
/// <param name="b">вектор правой части системы линейных уравнений</param>
/// <param name="x">вектор решения системы линейныз уравнений</param>
template <typename T>
void QRDecompositionSolve(MATRIX<T>& A, VECTOR<T>& b, VECTOR<T>& x)
{
	if (A.m_rows != A.m_columns || A.m_rows != b.m_size) return;
	MATRIX<T> Q(A.m_rows, A.m_columns), R(A.m_rows, A.m_columns);
	int n = A.m_rows;
	A.QRDecomposition(Q, R);
	VECTOR <T>beta(b.m_size);
	beta = Q.Transpose() * b;

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


}

Результаты решения системы

Был получен результат решения СЛАУ, вектор решения x:

$x = \begin{bmatrix} -5.99353 \\ -37.4183 \\ 0.934134 \\ 9.3746 \\ -0.0964064 \\ 1.74762 \end{bmatrix}$

Время решения СЛАУ и полученная норма невязки представлены в нижеследующей таблице:

Метод решения	Время решения, сек	Норма невязки
Метод Гаусса	$1.06\cdot10^{-5}$	$8.497\cdot10^{-12}$
LU декомпозиция	$1.49\cdot10^{-3}$	$1.02\cdot10^{-11}$
Компактная схема исключения	$1.19\cdot10^{-5}$	$7.048\cdot10^{-12}$
QR декомпозиция	$1.81\cdot10^{-5}$	$3.693\cdot10^{-12}$

Выводы

Как следует из результатов решения, приведённых в таблице выше, наибольшую точность даёт QR декомпозиция, а наилучший результат по времени решения даёт метод Гаусса. Самые худшие результаты и по времени решения, и по точности даёт LU декомпозиция, а наиболее оптимальным методом решения СЛАУ является компактная схема исключения.

Литература

Корн, Г. Справочник по математике для научных работников и инженеров : определения, теоремы, формулы / Г. Корн, Т. Корн ; пер. с англ. под общ. ред. И. Г. Арамановича. — 3‑е изд. — Москва : Наука, 1973. — 832 с.
Википедия. LU - разложение/ URL: https://ru.wikipedia.org/wiki/LU-разложение
Wikipedia. QR decomposition / URL: https://en.wikipedia.org/wiki/QR_decomposition (дата посещения 18.12.2022)
Ланкастер, П. Теория матриц / П. Ланкастер; пер. с англ. Демушкина С.П. — Москва : Наука, 1978. — 280 с.