Сравнение методов решения систем линейных алгебраических уравнений / Хабр

Эта статья является продолжением темы статьи https://habr.com/ru/articles/979542/. Данная статья не имеет цели представить особенности реализации методов решения систем линейных алгебраических уравнений (СЛАУ) или сделать какие-либо открытия в области численных методов. Цель данной статьи - предоставить сравнение методов решения СЛАУ и их эффективности. Оценка эффективности методов решения осуществлялась по двум параметрам: евклидовой норме невязки и времени решения. Евклидова норма невязки решения системы из n уравнений Ax = b вычисляется по формуле:

$\parallel Ax-b\parallel=\sqrt{\sum\limits_{i=0}^n [\sum\limits_{j=0}^n{a_{ij}\cdot x_j} - b_i]^2 }$

Ясно, что чем меньше евклидова норма невязки и чем меньше время решения, тем более эффективен метод для решения СЛАУ.

В статье https://habr.com/ru/articles/979542/ были представлены результаты оценки по времени решения и норме невязки четырьмя методами: методом Гаусса, LU декомпозиции [2], компактной схемой исключения [1] и QR декомпозиции [3,4] только для системы из шести уравнений, что не даёт представления об эффективности того или иного метода. В той же статье было представлено краткое описание этих методов. Поскольку для снижения затрат времени на решение эффективнее было бы обращаться к данным напрямую, то применение классов MATRIX<> и VECTOR<> оказалось лишённым смысла и при программной реализации методов хранение данных осуществлялось в массивах.

Реализация метода Гаусса:

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

	// сведение исходной системы к системе с верхней треугольной матрицей
	double* alf = new double[size * size*sizeof(double)];
	double* bet = new double[size * sizeof(double)];
	memcpy(alf, a, size*size*sizeof(double));
	memcpy(bet, b, size * sizeof(double));
	for (m = 0; m <= size - 2; m++)
	{
		amm = *(alf + m*size + m);
		if (std::abs(amm) <= DBL_MIN)
		{
			delete[] alf;
			delete[] bet;
			return 0.0;
		}
		for (k = m; k <= size - 1; k++)
			*(alf + m*size + k) /= amm;
		*(bet +m) /= amm;
		for (i = m + 1; i <= size - 1; i++)
		{
			aim = *(alf + i*size  + m);
			for (k = m; k <= size - 1; k++)
				*(alf + i*size + k) -= *(alf +m*size + k) * aim;
			*(bet +i) -= *(bet +m) * aim;
		}//end i 
	}//end m 

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

	// вычисление определителя
	double det = 1.0;
	for (i = 0; i < size; i++)
		det *= *(alf + i * size + i);

	delete[] alf;
	delete[] bet;
	return det;
}

Реализация решения СЛАУ прим��нением LU декомпозиции:

/// <summary>
/// Декомпозиция матрицы A = LU
/// L - нижняя треугольная матрица, на главной диагонали которой расположены единицы
/// U - верхняя треугольная матрица
/// </summary>
/// <param name="alfa">Матрицы L и U, ниже гланой диагонали которой расположены внедиагональные элементы L, 
/// на главной диагонали и выше расположены элементы матрицы U</param>
/// <param name="A">исходная матрица</param>
/// <param name="n">порядок матрицы</param>
/// <returns>Определитель матрицы</returns>
double LUDecomposition(const double* A, double* alfa, int n)
{

	// декомпозиция матрицы
	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++)
		{
			if (i <= j)
			{
				*(alfa + i * n + j) = *(A + i * n + j);
				for (int k = 0; k <= i - 1; k++)
					*(alfa + i * n + j) -= *(alfa + i * n + k) * *(alfa + k * n + j);
			}
			else
			{
				if (abs(*(alfa + j * n + j)) <= DBL_MIN) return 0;
				*(alfa + i * n + j) = *(A + i * n + j);
				for (int k = 0; k <= j - 1; k++)
					*(alfa + i * n + j) -= *(alfa + i * n + k) * *(alfa + k * n + j);
				*(alfa + i * n + j) /= *(alfa + j * n + j);

			}
		}
	double det = 1;
	for (int i = 0; i < n; i++)
		det *= *(alfa + i * n + i);
	return det;
}

/// <summary>
/// Решение системы линейных алгебраических уравнений с применением метода LU декомпозиции A = LU
/// L - нижняя тругольная матрица
/// U - верхняя треугольная матрица 
/// </summary>
/// <param name="A">матрица</param>
/// <param name="b">вектор правой части системы</param>
/// <param name="x">решение системы</param>
void LUDecompositionSolve(const double* A, const double* b, double* x, int n)
{

	double *alfa = new double[n*n*sizeof(double)];
	
	double det = LUDecomposition(A, alfa, n);
	if (abs(det) <= DBL_MIN) return;
	// решение СЛАУ
	double* y = new double[n * sizeof(double)];
	for (int k = 0; k < n; k++)
	{
		*(y + k) = *(b + k);
		for (int j = 0; j <= k - 1; j++)
			y[k] -= *(alfa + k*n +j) * *(y + j);
	}

	for (int k = n - 1; k >= 0; k--)
	{
		*(x + k) = *(y + k);
		for (int j = n - 1; j > k; j--)
			*(x + k) -= *(alfa + k*n + j) * *(x + j);
		*(x + k) /= *(alfa + k*n + k);
	}

	delete[] y;
	delete[] alfa;
}

Реализация компактной схемы исключения:

/// <summary>
/// Формирование матриц alpha и gamma в компактной схеме исключения
/// Матрица alpha верхняя треугольная, gamma - нижняя треугольная
/// </summary>
/// <param name="alpha">формируемые матрицы</param>
/// <param name="n">порядок матриц</param>
/// <returns>значение определителя матрицы</returns>
double FormMatrixCompactScheme(const double *A, double *alpha, int n)
{
	for (int i = 0; i < n; i++)
	{
		*(alpha + i * n) = *(A + i * n);
		if (i > 0) *(alpha + i) = *(A + i) / *A;
	}

	int k = 1, i = 1;
	while (i < n)
	{
		if (k >= n)
		{
			k = 1; i++;
			if (i >= n) break;
		}
		if (i >= k)
		{
			*(alpha + i * n + k) = *(A + i * n + k);
			for (int j = 0; j <= k - 1; j++)
				*(alpha + i * n + k) -= *(alpha + i * n + j) * *(alpha + j * n + k);
		}
		else
		{
			if (std::abs(*(alpha + i * n + i)) <= DBL_MIN) return 0.0;
			*(alpha + i * n + k) = *(A + i * n + k);
			for (int j = 0; j <= i - 1; j++)
			{
				*(alpha + i * n + k) -= *(alpha + i * n + j) * *(alpha + j * n + k);
			}
			*(alpha + i * n + k) /=  *(alpha + i * n + i);
		}
		k++;
	}

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

/// <summary>
/// Решение системы линейных уравнений компактной схемой исключения
/// </summary>
/// <param name="A">матрица системы уравнений</param>
/// <param name="b">вектор правой части системы уравнений</param>
/// <param name="x">решение системы уравнений</param>
/// <param name="n">порядок матрицы системы</param>
void CompactSchemeSolve(const double* A, const double* b, double* x, int n)
{
	double* alpha = new double[n*n*sizeof(double)];
	double det = FormMatrixCompactScheme(A, alpha, n);
	// решение только для неособенной матрицы
	if (abs(det) <= DBL_MIN)
	{
		delete[] alpha;
		return;
	}
	double *beta = new double[n*sizeof(double)];
	*beta = *b / *A;
	for (int i = 1; i < n; i++)
	{
		*(beta + i) = *(b + i);
		for (int j = 0; j <= i - 1; j++)
			*(beta + i) -= *(alpha + i*n + j) * *(beta +j);
		*(beta +i) /= *(alpha + i*n + i);
	}

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


	delete[] alpha;
	delete[] beta;
}

Реализация решения СЛАУ применением QR декомпозиции:

/// <summary>
/// QR разложение квадратной матрицы A, A = QR
/// где Q - ортогональная, а R - верхняя треугольная матрица
/// </summary>
/// <param name="Q">ортогональная матрица Q</param>
/// <param name="R">верхняя треугольная матрица R</param>
/// <param name="A">исходная матрица</param>
/// <param name="n">порядок матрицы A</param>
/// <returns>false - если на главной диагонали матрицы R 
/// хотя бы один элемент равен нулю (A - вырожденная матрица), иначе - true</returns>
bool QRDecomposition(const double* A, double* Q, double* R, int n)
{
	double sum = 0.0;

	for (int j = 0; j < n; j++)
	{
		// q(j) = a(j)
		for (int k = 0; k < n; k++)
			*(Q + k * n + j) = *(A + 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 + k * n + i) * *(A + k * n + j);
			*(R + i * n + j) = sum;

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

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

		if (abs(*(R + j * n + j)) <= DBL_MIN) return false;

		for (int k = 0; k < n; k++)
			*(Q + k * n + j) /= *(R + j * n + j);

	}

	return true;
}

/// <summary>
/// Решение СЛАУ с применением QR декомпозиции матрицы системы A
/// </summary>
/// <param name="A">матрица системы линейных уравнений</param>
/// <param name="b">вектор правой части системы линейных уравнений</param>
/// <param name="x">вектор решения системы линейных уравнений</param>
/// <param name="n">порядок матрицы A</param>
void QRDecompositionSolve(const double* A, const double* b, double* x, int n)
{
	double* Q = new double[n * n * sizeof(double)];
	double* R = new double[n * n * sizeof(double)]; 
	if (!QRDecomposition(A, Q, R, n))
	{
		delete[] Q;
		delete[] R;
		return;
	}
	double * beta = new double[n * sizeof(double)];
	//  beta = Qt*b
	for (int i = 0; i < n; i++)
	{
		*(beta + i) = 0.0;
		for (int k = 0; k < n; k++)
			*(beta + i) += *(Q + k * n + i) * *(b + k);

	}
	delete[] Q;

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

	delete [] R;
	delete[] beta;
}

В данной статье приводится оценка решения СЛАУ по затратам времени и точности решения по норме невязки. Число уравнений системы составляет от нескольких сотен до нескольких тысяч .

Кроме реализации четырёх вышеупомянутых методов была представлена реализация методов:

декомпозиции (разложение Холецкого, метод квадратных корней);
Метод вращений;
Метод релаксации
Метод градиентного спуска

Разложение Холецкого

Разложение Холецкого ( декомпозиция, метод квадратных корней) [6] применяется для симметричных положительно определённых матриц и заключается в представлении матрицы СЛАУ A в виде произведения , где означает транспонированную L матрицу. Тогда исходная система уравнений сводится к решению двух систем уравнений с нижней треугольной матрицей и c верхней треугольной матрицей. Элементы матрицы L вычисляются по формулам:

$\begin{matrix} l_{11}=\sqrt{a_{11}} \quad l_{j1} = \frac{a_{j1}}{l_{11}} \\ l_{ii} = \sqrt{ a_{ii}-\sum\limits_{k=1}^{i-1} l_{ki}^2 } \quad i > 1 \\ l_{ij} = \frac{(a_{ij} - \sum\limits_{k=1}^{i-1} l_{ki}l_{kj})}{l_{ii}} \quad j>i \\ l_{ij} = 0 \quad i > j \end{matrix}$

Для того, чтобы использовать разложение Холецкого для решения СЛАУ с матрицей общего вида, исходную систему уравнений можно преобразовать к виду $\widetilde{A}x=\widetilde{b}$ , где $\widetilde{A} = A^TA$ и $\widetilde{b} = A^Tb$ . Матрица $\widetilde{A} = A^TA$ является симметричной и положительно определённой.

Реализация решения СЛАУ применением разложения Холецкого (метода квадратных корней):

/// <summary>
///  Разложение Холецкого симметричной положительно определённой матрицы A = L*L^t
/// </summary>
/// <param name="L">матрица в разложении Холецкого</param>
/// <param name="A">исходная матрица</param>
/// <param name="n">порядок матрицы A</param>
/// <returns></returns>
bool CholeskyDecomposition(const double*A, double* L, int n)
{
	for (int i = 0; i < n; i++)
		for (int j = 0; j <= i; j++) {

			double sum = 0;
			for (int k = 0; k < j; k++)
				sum += *(L + i * n + k) * *(L + j * n + k);

			if (i == j)
			{
				*(L + i * n + j) = sqrt(*(A + i * n + i) - sum);
				if ( isnan(*(L + i * n + j)) ) return false;
			}
			else
				*(L + i * n + j) = (*(A + i * n + j) - sum) / (*(L + j * n + j));
		}

	return true;
}

/// <summary>
/// Решение системы уравнений с применением метода LLT декомпозиции A = LLT
/// Метод применим только для положительно определённых симметричных матриц
/// </summary>
/// <param name="A">матрица</param>
/// <param name="b">вектор правой части</param>
/// <param name="x">вектор решения СЛАУ</param>
void LLTDecompositionSolve(const double* A, const double* b, double* x, int n)
{	
	
	double* L = new double[n * n * sizeof(double)];
	memset(L, 0, n * n * sizeof(double));
	CholeskyDecomposition(A, L, n);

	// решение системы с нижней треугольной матрицей Ly = b
	double* y = new double[n * sizeof(double)];
	for (int i = 0; i < n; i++)
	{
		*(y + i) = *(b + i);
		for (int k = 0; k < i; k++)
			*(y + i) -= *(L +i*n + k) * *(y + k);
		*(y + i) /= *(L + i * n + i);

	}

	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++)
			if (i < j) *(L + i * n + j) = *(L + j * n + i);

	// решение системы уравнений с верхней труегольной матрицей L^t*x = y   
	for (int i = n - 1; i >= 0; i--)
	{
		*(x + i) = *(y + i);
		for (int j = i + 1; j < n; j++)
			*(x + i) -= *(L + i * n + j) * *(x + j);
		*(x + i) /= *(L + i * n + i);
	}


	delete[] y;
	delete[] L;
}

Метод вращений

Метод вращения применяется для нахождения собственных значений симметричной матрицы [5,6]. В данной статье метод вращений был модифицирован таким образом, чтобы за k шагов последовательных вращений матрица общего вида A сводилась к верхней треугольной матрице:

$\begin{matrix} T^{(0)} = A \\ T^{(k)} = U^{(k)}T^{(k-1)} \\ T^{(k)} = U^{(k)}...U^{(1)}U^{(0)}A \\ \end{matrix}$

где $T^{(k)}$ - итерации исходной матрицы, сводящие исходную матрицу A к верхней треугольной матрице, $U^{(k)}$ - матрица вращения:

$U^{(k)} = \begin{cases} 1 & i = j , i \neq i_0, j \neq j_0\\ cos(\phi^{(k)}) & i = i_0, j= i_0 \\ sin(\phi^{(k)}) & i = i_0, j= j_0 \\ -sin(\phi^{(k)}) & i = j_0, j= i_0 \\ cos(\phi^{(k)}) & i = j_0, j= j_0 \\ 0 & i \neq j, i \neq i_0, j \neq j_0 \end{cases}$

где на k-ом шаге итерации строка , столбец и угол вращения $\phi^{(k)}$ выбираются из условия $t^{(k-1)}_{i_0j_0} = 0$ , где $t^{(k-1)}_{i_0j_0}$ - максимальный по модулю элемент, расположенный ниже главной диагонали матрицы $T^{(k-1)}$ . Из этого условия можно получить значение угла вращения $\phi^{(k)}$ для вычисления значений матрицы $U^{(k)}$ :

$\phi^{(k)} = arctg(-\frac{t^{(k-1)}_{i_0j_0}}{t^{(k-1)}_{i_0i_0}})$

Элементы матрицы $T^{(k)}$ на k-ом шаге итерации отличаются от элементов матрицы $T^{(k-1)}$ на предыдущем шаге итерации только двумя строками:

$\begin{matrix} t^{(k)}_{i_0m} = t^{(k-1)}_{i_0m}\cdot cos(\phi^{(k)}) + t^{(k-1)}_{j_0m}\cdot sin(\phi^{(k)}) \\ t^{(k)}_{j_0m} = -t^{(k-1)}_{i_0m}\cdot sin(\phi^{(k)}) + t^{(k-1)}_{j_0m}\cdot cos(\phi^{(k)}) \end{matrix}$

Вектор правой части СЛАУ b преобразуется приводится к вектору $\beta$ :

$\begin{matrix} \beta^{(0)} = b \\ \beta^{(k)} = U^{(k)}\beta^{(k-1)} \\ \beta^{(k)} = U^{(k)}...U^{(1)}U^{(0)}b \end{matrix}$

Только две компоненты вектора $\beta^{(k)}$ на k-ом шаге итерации отличаются от компонент $\beta^{(k-1)}$ на предыдущем шаге итерации:

$\begin{matrix} \beta^{(k)}_{i_0} = \beta^{(k-1)}_{i_0}\cdot cos(\phi^{(k)}) + \beta^{(k-1)}_{j_0}\cdot sin(\phi^{(k)}) \\ \beta^{(k)}_{j_0} = -\beta^{(k-1)}_{i_0}\cdot sin(\phi^{(k)}) + \beta^{(k-1)}_{j_0}\cdot cos(\phi^{(k)}) \end{matrix}$

Итерационный процесс завершается при условии, когда значение максимального по модулю элемента $t^{(k)}_{i_0j_0}$ матрицы $T^{(k)}$ становится меньше определённого малого значения $\epsilon$ : $\mid t^{(k)}_{i_0j_0}\mid < \epsilon$ .Таким образом, исходная система уравнений Ax = b сводится к системе уравнений $T^{(k)}x = \beta^{(k)}$ , где матрица $T^{(k)}$ близка к верхней треугольной матрице. Реализация метода приводится ниже:

/// <summary>
/// Решение СЛАУ методом вращения, решение сводится
/// к решению СЛАУ с верхней треугольной матрицей
/// </summary>
/// <param name="A">матрица СЛАУ</param>
/// <param name="b">вектор правой части СЛАУ</param>
/// <param name="x">вектор решения СЛАУ</param>
/// <param name="n">порядок матрицы</param>
/// <returns>true - если сходимость была достигнута, false - иначе</returns>
bool RotationSolve(const double *A, const double *b, double *x, int n)
{

	// верхняя тругольная матрица
	double* T = new double[n * n * sizeof(double)]; // на текущей итерации

	double* T0 = new double[n * n * sizeof(double)]; // на предыдущей итерации
	memcpy(T0, A, n * n * sizeof(double));

	// преобразованный вектор правой части
	double* bet0 = new double[n * sizeof(double)]; // на предыдущей итерации
	memcpy(bet0, b, n * sizeof(double));

	double* bet = new double[n * sizeof(double)]; // на текущей итерации
	int iter = 0;
	while (true)
	{
		// поиск позиции максимального по модулю элемента ниже главной диагонали
		int i0 = 0, j0 = 0;
		double val = 0.0, el = 0.0;
		for (int i = 0; i < n; i++)
			for (int j = 0; j < i; j++)
			{
				el = abs(*(T0 + i * n + j));
				if (el > val)
				{
					val = el;
					i0 = i;
					j0 = j;
				}
				
			}
		if (val < DBL_EPSILON) break;
		
		// угол матрицы вращения, находится по условию T(k)(i0, j0) = 0
		double fi = atan((-1.0) * *(T0 + i0 * n + j0) / *(T0 + j0 * n + j0));
		double cs = cos(fi);
		double ss = sin(fi);

		// T(k) = U*T(k-1), U - матрица вращения
		// bet(n) = U*bet(n-1)
		memcpy(T, T0, n * n * sizeof(double));
		memcpy(bet, bet0, n * sizeof(double));
		for (int i = 0; i < n; i++)
		{
			*(T + i0 * n + i) = *(T0 + i0 * n + i) * cs +
				*(T0 + j0 * n + i) * ss;
			*(T + j0 * n + i) = (-1.0) * *(T0 + i0 * n + i) * ss +
				*(T0 + j0 * n + i) * cs;
		}

		*(bet + i0) = *(bet0 + i0) * cs +
			*(bet0 + j0) * ss;
		*(bet + j0) = (-1.0) * *(bet0 + i0) * ss +
			*(bet0 + j0) * cs;

		// результаты для следующей итерации
		memcpy(T0, T, n * n * sizeof(double));
		memcpy(bet0, bet, n * sizeof(double));
		if (++iter > MAX_ITERATION_NUMBER)
		{
			delete[] T0;
			delete[] bet0;
			delete[] T;
			delete[] bet;
			return false;
		}
	}

	delete[] T0;
	delete[] bet0;

	// вычисление вектора решения
	// решение системы уравнений с верхней труегольной матрицей A(n)*x = bet   
	for (int i = n - 1; i >= 0; i--)
	{
		*(x + i) = *(bet + i);
		for (int j = i + 1; j < n; j++)
			*(x + i) -= *(T + i * n + j) * *(x + j);
		*(x + i) /= *(T + i * n + i);
	}

	delete[] T;
	delete[] bet;
	return true;
}

Метод релаксации

Метод релаксации применяется для положительно определённых симметричных матриц, описание представлено в источнике [1]. Метод заключается в последовательности k итераций компонент вектора решения x:

${x_i}^{(k)} = {x_i}^{(k-1)} -\frac{\omega}{a_{ii}}\cdot (\sum\limits_{j=0}^{i-1}a_{ij}x_j^{(k)}) +\sum\limits_{j=i}^{n}a_{ij}x_j^{(k-1)} - b_i)$

Итерационный процесс прекращается при условии когда евклидова норма разности между вектором решения на текущей и на предыдущей итерации достигает определённого малого положительного значения $\epsilon$ :

$\parallel x^{(k)} - x^{(k-1)}\parallel < \epsilon$

Реализация метода релаксации:

/// <summary>
/// Скалярное произведение векторов (v1, v2)
/// </summary>
/// <param name="v1">вектор</param>
/// <param name="v2">вектор</param>
/// <param name="n">размерность векторов v1 и v2</param>
/// <returns>скалярное произведение векторов v1 и v2</returns>
double ScalarProduct(const double *v1, const double *v2, int n)
{ 
	double result = 0.0;
	for (int i = 0; i < n; i++)
		result += *(v1 + i) * *(v2 + i);
	return result;
}

/// <summary>
/// Вычисление сферической нормы n-мерного вектора v
/// </summary>
/// <param name="v">вектор</param>
/// <param name="n">размерность вектора</param>
/// <returns>значение нормы</returns>
double norm(const double* v, int n)
{
	return sqrt(ScalarProduct(v, v, n));
}

/// <summary>
/// Решение системы уравнений A*x = b методом  релаксации
/// для симметричной положительно определённой матрицы A
/// </summary>
/// <param name="A">матрица СЛАУ</param>
/// <param name="b">вектор правой части</param>
/// <param name="x">вектор решения СЛАУ</param>
/// <param name="n">порядок матрицы A</param>
/// <param name="omega">релаксационный множитель</param>
/// <returns>true - если сходимость была достигнута, false - иначе</returns>
bool Relaxation(const double* A, const double* b, double* x, int n, double omega)
{
	double* x0 = new double[n * sizeof(double)];
	double *errv = new double[n * sizeof(double)];
	memcpy(x0, b, n * sizeof(double));

	int iter = 0; // число итераций
	double err_norm = 0.0; // погрешность
	do
	{		
		for (int i = 0; i < n; i++)
		{
			double sum = 0.0;
			for (int k = 0; k <= i - 1; k++)
				sum += *(A + i * n + k) * *(x + k);

			for (int k = i; k < n; k++)
				sum += *(A + i * n + k) * *(x0 + k);

			*(x + i) = *(x0 + i) - omega * (sum - *(b + i)) / *(A + i * n + i);						
		}
		// норма разности между значением на текущей и предыдущей итерации
		for(int k =0; k < n; k++)
			*(errv + k) = *(x + k) - *(x0 + k);
		err_norm = norm(errv, n);
		if (++iter > MAX_ITERATION_NUMBER)
		{
			delete[] errv;
			delete[] x0;
			return false;
		}
		memcpy(x0, x, n * sizeof(double));
	} while (err_norm > DBL_EPSILON);
	
	delete[] errv;
	delete[] x0;
	return true;
}

Метод градиентного спуска

Описание метода градиентного спуска представлено в монографии [5]. Суть метода заключается в поиске минимума функционала $F(x) = (Ax,x) - 2(b,x) \quad x \in R^n$ , (x,y) - скалярное произведение векторов x и y:

$(x,y) = \sum\limits_{i=1}^n x_iy_i$

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

$x^{(k)} = x^{(k-1)} - t\cdot grad(F(x^{(k-1)})\quad t>0$

Значение t на k-ом шаге итерации находится из условия:

$\begin{matrix} \phi(t) = x^{(k-1)} - t\cdot grad(F(x^{(k-1)}),\quad \frac{\text{d}{\phi(t)}}{\text{d}t} =0 \\ grad(F(x^{(k-1)}) = 2\cdot (Ax^{(k-1)}-b), \end{matrix}$

откуда:

$\begin{matrix} t = \frac{(r^{(k-1)},r^{(k-1)})}{2\cdot (r^{(k-1)},Ar^{(k-1)})}, \\ где\quad r^{(k-1)} = A\cdot x^{(k-1)} - b \end{matrix}$

$\parallel x^{(k)} - x^{(k-1)}\parallel < \epsilon$

Реализация метода градиентного спуска:

/// <summary>
/// Скалярное произведение векторов (v1, v2)
/// </summary>
/// <param name="v1">вектор</param>
/// <param name="v2">вектор</param>
/// <param name="n">размерность векторов v1 и v2</param>
/// <returns>скалярное произведение векторов v1 и v2</returns>
double ScalarProduct(const double *v1, const double *v2, int n)
{ 
	double result = 0.0;
	for (int i = 0; i < n; i++)
		result += *(v1 + i) * *(v2 + i);
	return result;
}

/// <summary>
/// Вычисление сферической нормы n-мерного вектора v
/// </summary>
/// <param name="v">вектор</param>
/// <param name="n">размерность вектора</param>
/// <returns>значение нормы</returns>
double norm(const double* v, int n)
{
	return sqrt(ScalarProduct(v, v, n));
}

/// <summary>
/// Решение системы уравнений A*x = b методом градиентного спуска
/// для симметричной положительно определённой матрицы A
/// </summary>
/// <param name="A">матрица СЛАУ</param>
/// <param name="b">вектор правой части</param>
/// <param name="x">вектор решения СЛАУ</param>
/// <param name="n">порядок матрицы A</param>
/// <param name="omega">релаксационный множитель</param>
/// <returns>true - если сходимость была достигнута, false - иначе</returns>
bool GradientDescent(const double* A, const double* b, double* x, int n)
{

	// градиент: 2*(Ax0 - b)
	auto Gradient = [A, b, n](const double* x, double *result)
	{
		double* ax = new double[n * sizeof(double)];
		memset(ax, 0, n * sizeof(double));
		for (int i = 0; i < n; i++)
		{
			for (int j = 0; j < n; j++)
				*(ax + i) += *(A + i * n + j) * *(x + j);
			*(result + i) = 2.0 * (*(ax + i) - *(b + i));
		}
		delete[] ax;
	};

	// A*(Ax0 - b)
	auto VectorRight = [A, b, n](const double* x, double* result)
	{
		double* ax = new double[n * sizeof(double)];
		memset(ax, 0, n * sizeof(double));
		for (int i = 0; i < n; i++)
		{
			for (int j = 0; j < n; j++)
				*(ax + i) += *(A + i * n + j) * *(x + j);
			*(ax + i) -= *(b + i);
		}

		memset(result, 0, n * sizeof(double));
		for (int i = 0; i < n; i++)
			for(int j = 0; j < n; j++)
				*(result + i) += *(A + i * n + j) * *(ax + j);
		delete[] ax;
	};

	double *x0 = new double[n * n * sizeof(double)];
	memset(x0, 0, n * sizeof(double));
	*x0 = 1.0;
	double* errv = new double[n * sizeof(double)];
	double *grad = new double[n * sizeof(double)];
	double* agrad = new double[n * sizeof(double)];

	int iter = 0;
	double err_norm = 0.0, t = 0.0;
	do
	{ 
		Gradient(x0, grad);
		double sc1 = ScalarProduct(grad, grad, n) * 0.25;

		VectorRight(x0, agrad);
		double sc2 = 2.0 * ScalarProduct(grad, agrad, n);
		t = sc1 / sc2;

		for (int i = 0; i < n; i++)
			*(x + i) = *(x0 + i) - t * *(grad + i);
		// норма разности между значением на текущей и предыдущей итерации
		for (int i = 0; i < n; i++)
			*(errv + i) = *(x + i) - *(x0 + i);
		err_norm = norm(errv, n);
		if (++iter > MAX_ITERATION_NUMBER)
		{
			delete[] errv;
			delete[] x0;
			delete[] grad;
			delete[] agrad;
			return false;
		}
		memcpy(x0, x, n * sizeof(double));

	} while (err_norm > DBL_EPSILON);
	delete[] errv;
	delete[] x0;
	delete[] grad;
	delete[] agrad;
	return true;
}

Результаты вычислений

Матрица СЛАУ была заполнена на 100%, т. е. не содержала нулевых элементов; коэффициент релаксации $\omega = 0.9$ ; максимальное число итерации для методов вращения, релаксации и градиентного спуска равнялось 60 000 и $\epsilon = 2.2204\cdot 10^{-16}$ (DBL_EPSILON). Результаты вычислений представлены в нижеследующей таблице:

Метод	Время, сек	Норма невязки
Число уравнений 50
Метод Гаусса	0.0001875	$1.02721\cdot 10^{-9}$
LU декомпозиция	0.0001949	$8.14592\cdot 10^{-10}$
Компактная схема исключения	0.0002154	$1.46857\cdot 10^{-9}$
QR декомпозиция	0.0004497	$1.07459\cdot 10^{-10}$
декомпозиция	0.0005565	$9.07323\cdot 10^{-11}$
Метод вращений	0.228928	$5.56271\cdot 10^{-11}$
Метод релаксации и метод градиентного спуска	плохая сходимость
Число уравнений 100
Метод Гаусса	0.0010701	$8.73551\cdot 10^{-10}$
LU декомпозиция	0.0012279	$1.33149\cdot 10^{-9}$
Компактная схема исключения	0.0010931	$8.73551\cdot 10^{-10}$
QR декомпозиция	0.0028991	$5.61553\cdot 10^{-11}$
декомпозиция	0.0035761	$6.66858\cdot 10^{-12}$
Метод вращений и метод релаксации	плохая сходимость
Метод градиентного спуска	3.07913	$3.21141\cdot 10^{-10}$
Число уравнений 200
Метод Гаусса	0.0066878	$1.49908\cdot 10^{-8}$
LU декомпозиция	0.007915	$8.12037\cdot 10^{-9}$
Компактная схема исключения	0.0085465	$1.49908\cdot 10^{-08}$
QR декомпозиция	0.0207982	$1.5519\cdot 10^{-10}$
декомпозиция	0.0271643	$1.74292\cdot 10^{-10}$
Число уравнений 500
Метод Гаусса	0.105087	$1.50375\cdot 10^{-7}$
LU декомпозиция	0.131105	$1.01746\cdot 10^{-7}$
Компактная схема исключения	0.156237	$8.80109\cdot 10^{-8}$
QR декомпозиция	0.413718	$7.76066\cdot 10^{-10}$
декомпозиция	0.44497	$1.52446\cdot 10^{-9}$
Число уравнений 1000
Метод Гаусса	0.752816	$3.194\cdot 10^{-7}$
LU декомпозиция	1.02028	$1.39641\cdot 10^{-6}$
Компактная схема исключения	1.02911	$9.31886\cdot 10^{-7}$
QR декомпозиция	3.61937	$3.50069\cdot 10^{-9}$
декомпозиция	3.81145	$3.07008\cdot 10^{-9}$
Число уравнений 2000
Метод Гаусса	6.14911	$1.92358\cdot 10^{-5}$
LU декомпозиция	10.7675	$3.8148\cdot 10^{-5}$
Компактная схема исключения	11.4573	$1.92358\cdot 10^{-5}$
QR декомпозиция	51.5436	$1.48766\cdot 10^{-8}$
декомпозиция	49.5564	$5.38007\cdot 10^{-9}$
Число уравнений 3000
Метод Гаусса	20.6088	0.000141378
LU декомпозиция	37.7988	0.000180529
Компактная схема исключения	37.5891	0.000143626
QR декомпозиция	192.067	$3.61746\cdot 10^{-7}$
декомпозиция	189.744	$4.51635\cdot 10^{-6}$

По результатам расчётов можно выделить три группы методов:

1) Метод Гаусса, компактная схема исключения и решение с применением LU декомпозиции.

2) Методы решения применением QR и декомпозиции.

3) Итерационные методы: метод вращений, метод релаксации и метод градиентного спуска.

Методы первой группы в сравнении с методами второй группы отличаются меньшим временем решения и меньшей точностью, которая определяется нормой невязки. Причем время решения первой группы в разы меньше, но при этом точность решения второй группы на порядок или два выше точности решения первой группы методов. Наименьшее время решения среди методов первой и второй группы демонстрирует метод Гаусса, но уже начиная с 2000 уравнений точность метода Гаусса становится неудовлетворительной.

Итерационные методы третьей группы непригодны для заполненных матриц и плохо сходятся для заданной точности $\epsilon = 2.2204\cdot 10^{-16}$ . Их применение эффективно на разреженных матрицах, о чём пойдёт речь в следующем разделе статьи.

Реализация и сравнение решения методами вращений, релаксации и градиентного спуска

Методы вращений, релаксации и градиентного спуска не подходят для решения СЛАУ с полностью заполненной матрицей. В этом разделе приводится реализация этих методов для разреженных матриц СЛАУ. Разреженные матрицы характеризуются степенью заполненности матрицы - отношением числа ненулевых элементов к числу элементов матрицы.

Для разреженной матрицы не имеет смысла хранить нулевые элементы, поэтому для хранения и обработки данных использовались не массивы, а другие структуры хранения данных. Для хранения элементов матрицы или компонент вектора была использована структура:

// элемент матрицы
struct SparseElement
{
	int row; // строка
	int column; // столбец
	double value; // значение
};

Для представления данных матрицы или вектора использовался STL контейнер vector<>, в частности vector<SparseElement>. Такая структура данных позволяет единообразно осуществлять алгебраические операции как с матрицами, так и с векторами, представляющими частный случай матрицы (каждый компонент вектора имеет значение SparseElement::column = 1). Использование контейнера vector<> также позволяет задействовать инструментарий STL для поиска и сортировки массива данных.

Краткое описание методов вращений, релаксации и градиентного спуска было приведено в предыдущих разделах, в этом разделе приводится реализация методов для разреженных матриц:

#include <iostream>
#include <vector>
#include <algorithm>
#include <iterator>
#include <map>

#define MAX_ITERATION_NUMBER 60000
#define EPS 1.0e-12

// элемент матрицы
struct SparseElement
{
	int row; // строка
	int column; // столбец
	double value; // значение
};

using namespace std;

/// <summary>
/// Поиск значения элемента в строке i и столбце j
/// </summary>
/// <param name="matrix">разреженная матрица</param>
/// <param name="_row">строка</param>
/// <param name="_column">столбец</param>
/// <returns>ненулевое значение элемента, 0 - если элемент не найден</returns>
double GetValue(const vector<SparseElement>& matrix,int _row, int _column = 1)
{
	auto iter = find_if( matrix.begin(), matrix.end(),
		[_row,_column](const SparseElement &el) { return el.row == _row && el.column == _column; }
	);
	return iter != matrix.end() ? iter->value : 0.0;
} 

/// <summary>
/// Присвоить значение элемента матрицы: если найден - присвоить значение,
/// если нет и значение ненулевое - создать элемент
/// </summary>
/// <param name="matrix">матрица или вектор</param>
/// <param name="_row">строка</param>
/// <param name="_column">столбец</param>
/// <param name="_value">значение</param>
void SetValue(vector<SparseElement>& matrix, int _row, int _column, double _value)
{
	auto iter = find_if(matrix.begin(), matrix.end(),
		[_row, _column](const SparseElement& el) { return el.row == _row && el.column == _column; }
	);
	if (iter != matrix.end())
		iter->value = _value;
	else
	{
		if (abs(_value) <= DBL_MIN) return;
		matrix.push_back({ _row,_column, _value });
	}
}

/// <summary>
/// Транспонирование матрицы
/// </summary>
/// <param name="matrix">исходная матрица</param>
/// <param name="n">порядок матрицы</param>
/// <returns>транспонированну  matrix матрицу</returns>
vector<SparseElement> SparseTranspose(const vector<SparseElement>& matrix)
{
	vector <SparseElement> result;
	for (const auto& elem : matrix)
		result.push_back({ elem.column, elem.row, elem.value });

	return result;
}

/// <summary>
/// Умнoжение матриц first и second
/// </summary>
/// <param name="first">матрица слева</param>
/// <param name="second">матрица справа</param>
/// <returns>произведение матриц</returns>
vector <SparseElement> SparseMultiply(const vector<SparseElement>& first,
	const vector<SparseElement>& second)
{
	map<int, map<int, double>> matrix;
	vector <SparseElement> result;

	for (const auto& el_first : first)
		for(const auto& el_second : second)
			if (el_first.column == el_second.row)
				matrix[el_first.row][el_second.column] += el_first.value * el_second.value;

	for (const auto& matrix_el : matrix)
		for (const auto& mrow_el : matrix_el.second)
			result.push_back({ matrix_el.first, mrow_el.first, mrow_el.second });

	return result;
}

/// <summary>
/// Разность двух матриц first и second
/// </summary>
/// <param name="first">разреженная матрица</param>
/// <param name="second">разреженная матрица</param>
/// <returns>разреженную матрицу first - second</returns>
vector<SparseElement> SparseDifference(const vector<SparseElement>& first,
	const vector<SparseElement>& second)
{
	vector<SparseElement> result;
	
	for (const auto& elem_first : first)
		result.push_back({ elem_first.row, elem_first.column,elem_first.value });
	for (const auto& elem_second : second)
	{
		int i = elem_second.row;
		int j = elem_second.column;
		auto iter = find_if(result.begin(), result.end(),
			[i, j]
			(const SparseElement& el) { return el.row == i && el.column == j; }
		);
		if(iter != result.end())
			iter->value -= elem_second.value;
		else 
			result.push_back({ elem_second.row, elem_second.column, -elem_second.value });
	}

	return result;
}

/// <summary>
/// Скалярное произведение векторов vect1 и vect2
/// </summary>
/// <param name="vect1">вектор</param>
/// <param name="vect2">вектор</param>
/// <returns>значение скалярного произведения</returns>
double SparseScalarProduct(const vector<SparseElement>& vect1,
	const vector<SparseElement>& vect2)
{
	vector<SparseElement> value = SparseMultiply(SparseTranspose(vect1), vect2);
	return value[0].value;
}
/// <summary>
/// Евклидова (сферическая норма) вектора
/// </summary>
/// <param name="vect">разреженный вектор</param>
/// <returns>значение нормы</returns>
double SparseNormv(const vector<SparseElement>& vect) 
{
	return sqrt(SparseScalarProduct(vect, vect));
}
/// <summary>
/// для сортировки матриц и векторов
/// </summary>
/// <param name="one"></param>
/// <param name="two"></param>
/// <returns></returns>
bool VectorCompare(const SparseElement& one, const SparseElement& two)
{
	return one.row < two.row;
}

/// <summary>
/// Погрешность решения - норма невязки A*x - b решения СЛАУ
/// </summary>
/// <param name="A">разреженная матрица СЛАУ</param>
/// <param name="b">вектор правой части</param>
/// <param name="x">полученное решение СЛАУ</param>
/// <returns>погрешность</returns>
double ErrorMeasure(const vector<SparseElement>& A, 
					const vector<SparseElement>& b, 
					const vector<SparseElement>& x)
{
	// A*x
	vector <SparseElement> ax = SparseMultiply(A, x);

	return SparseNormv(SparseDifference(ax, b));
}

/// <summary>
/// Решение СЛАУ методом вращения c разреженной матрицей, решение сводится
/// к решению СЛАУ с верхней треугольной матрицей
/// </summary>
/// <param name="A">матрица СЛАУ</param>
/// <param name="b">вектор правой части СЛАУ</param>
/// <param name="x">вектор решения СЛАУ</param>
/// <param name="n">порядок матрицы и размерность вектора правой части</param>
/// <returns>
/// -1 - деление на ноль, нулевой элемент на главной диагонали
/// -2 - число элементов в строке матрицы не совпадает с размерностью вектора
/// или матрица не квадратная
/// -3 - элемент не найден
///   0 - превышено максимальное число итераций, 
/// 1 - решение завершено успешно
/// </returns>
int SparseRotationSolve(const vector<SparseElement>& A, 
						const vector<SparseElement>& b, 
						vector<SparseElement>& x,
						int n)
{
	// верхняя тругольная матрица
	vector<SparseElement> T; // на текущей итерации

	vector<SparseElement> T0 = A; // на предыдущей итерации

	// преобразованный вектор правой части
	vector<SparseElement> bet0 = b; // на предыдущей итерации

	vector<SparseElement> bet; // на текущей итерации

	double val = 0;
	int i0 = 0, j0 = 0;
	int iter_num = 0; // число итерации, ограничивается максимальным числом итераций
	while (true)
	{
		vector<SparseElement>::iterator iter;
		// любой элемент ниже главной диагонали
		for (iter = T0.begin(); iter != T0.end(); ++iter)
			if (iter->column < iter->row) break;
		// поиск позиции максимального по модулю элемента ниже главной диагонали
		 iter = max_element(iter, T0.end(),
			[](const SparseElement& one, const SparseElement& two)
			{
				if (one.column < one.row && two.column < two.row)
					return abs(one.value) < abs(two.value);
				else
					return false;
			}
		); 
		
		j0 = iter->column;
		i0 = iter->row;
		val = iter->value;

 		if (abs(val) < EPS) break;
		// угол матрицы вращения, находится по условию T(k-1)(i0, j0) = 0
		double t2 = GetValue(T0, j0, j0);
		if (t2 == 0.0) return -1;

		double fi = atan(-val / t2);
		double cs = cos(fi);
		double ss = sin(fi);

		// T(k) = U*T(k-1), U - матрица вращения
		T = T0;

		// заполнение строк T(i0,k) и T(j0,k)	
		double t_i0 = 0.0, t_j0 = 0.0, val_j0 = 0.0;;
		bool found_i0 = false, found_j0 = false;
		for (int k = 1; k <= n; k++)
		{
			// T0(i0,k), T0(j0,k)
			t_i0 = t_j0 = 0.0;
			for (const auto& elem : T0)
			{
				if (elem.row == i0 && elem.column == k)
					t_i0 = elem.value;
				if (elem.row == j0 && elem.column == k)
					t_j0 = elem.value;

			}

			found_i0 = found_j0 = false;
			val_j0 = 0.0;
			// T(i0,k) = T0(i0,k)*cos(phi) + T0(j0,k)*sin(phi)
			val = t_i0 * cs + t_j0 * ss;

			// T(j0,k) = -T0(i0,k)*sin(phi) + T0(j0,k)*cos(phi)
			val_j0 = -t_i0 * ss + t_j0 * cs;

			for (auto& elem : T)
			{
				if (elem.row == i0 && elem.column == k)
				{
					elem.value = val; 
					found_i0 = true;
				}
				if (elem.row == j0 && elem.column == k)
				{
					elem.value = val_j0; 
					found_j0 = true;
				}
			}

			if (!found_i0)
				T.push_back({ i0, k, val });
			if(!found_j0)
				T.push_back({ j0, k, val_j0 });
		}

		// bet(k) = U*bet(k-1)
		bet = bet0;

		// b0(i0), b0(j0)
		double b_i0 = 0.0, b_j0 = 0.0;
		for (const auto& elem : bet0)
		{
			if (elem.row == i0 && elem.column == 1)
				b_i0 = elem.value;
			if (elem.row == j0 && elem.column == 1)
				b_j0 = elem.value;
		}

		val = b_i0 * cs + b_j0 * ss;

		val_j0 = -b_i0 * ss + b_j0 * cs;

		found_i0 = found_j0 = false;
		for (auto& elem : bet)
		{
			if (elem.row == i0 && elem.column == 1)
			{
				elem.value = val;
				found_i0 = true;
			}
			if (elem.row == j0 && elem.column == 1)
			{
				elem.value = val_j0;
				found_j0 = true;
			}
		}

		if (!found_i0)
			bet.push_back({ i0, 1, val });
		if (!found_j0)
			bet.push_back({ j0, 1, val_j0 });

		// результаты для следующей итерации
		T0 = T;
		bet0 = bet;
		if (++iter_num > MAX_ITERATION_NUMBER) return 0;
	}
	
	T0.clear();
	bet0.clear();

	// вычисление вектора решения
	// решение системы уравнений с верхней треугольной матрицей T(n)*x = bet(n)   
	if (!x.empty()) x.clear();
	
	for (int i = n; i >= 1; i--)
	{ 
		val = GetValue(bet, i);
	
		double tval = 0.0, xval = 0.0;
		for (int j = i + 1; j <= n; j++)
		{
			tval = GetValue(T, i, j);
			if (tval == 0) continue;
			xval = GetValue(x, j);
			if (xval == 0) continue;
			val-= tval * xval;
		}

		tval = GetValue(T, i, i);
		if (abs(tval) <= DBL_MIN) return -1;

		if (abs(val) <= DBL_MIN) continue;

		x.push_back({i, 1, val/tval});
	}
	
	reverse(x.begin(), x.end());
	T.clear();
	bet.clear();	
	return 1;
}

/// <summary>
/// Решение системы уравнений A*x = b методом  релаксации
/// для симметричной положительно определённой матрицы A
/// </summary>
/// <param name="A">матрица СЛАУ</param>
/// <param name="b">вектор правой части</param>
/// <param name="x">вектор решения СЛАУ</param>
/// <param name="omega">релаксационный множитель</param>
/// <returns>
///  1 - если сходимость была достигнута,
/// -1 - элемент матрицы на главной диагонали равен нулю,
///  0 - превышено максимальное число итераций
/// </returns>
int SparseRelaxation(const vector<SparseElement>& A,
					const vector<SparseElement>& b, 
					vector<SparseElement>& x, 
					double omega)
{
	vector<SparseElement> x0;
	x0.push_back({ 1, 1, 1.0 });
	double err_norm0 = 1.0;
	map<int, map<int, double>> Amap;

	for (const auto& elem : A)
		Amap[elem.row][elem.column] = elem.value;

	int iter = 0; // число итераций
	double err_norm = 0.0; // погрешность
	do
	{		
		double sum = 0.0, val_d = 0.0;

		for (const auto& elem_a : Amap)
		{	
			map<int,double> Arow = elem_a.second;
			double sum = 0.0;
			int row = elem_a.first;
			
			for (const auto& elem_x : x)
			{
				if (elem_x.row <= row - 1 && Arow.count(elem_x.row) > 0)
					sum += Arow[elem_x.row] * elem_x.value;
			}

			for (const auto& elem_x : x0)
			{ 
				if (elem_x.row >= row && Arow.count(elem_x.row) > 0)
					sum += Arow[elem_x.row] * elem_x.value;
			}

			if (Arow.count(row) == 0) return -1;
			double val_d = Arow[row];	
			double val = GetValue(x0, row) - omega * (sum - GetValue(b, row)) / val_d;
			SetValue(x, row, 1, val);
		}


		// норма разности между значением на текущей и предыдущей итерации
		err_norm = SparseNormv(SparseDifference(x, x0));
		if (abs(err_norm - err_norm0) <= DBL_MIN) return 1;
		if (++iter > MAX_ITERATION_NUMBER) return 0;

		if (!is_sorted(x.begin(), x.end(), VectorCompare))
			sort(x.begin(), x.end(), VectorCompare);
		x0 = x;
		err_norm0 = err_norm;

	} while (err_norm > EPS);

	return 1;
}

/// <summary>
/// Решение системы уравнений A*x = b методом градиентного спуска
/// для симметричной положительно определённой матрицы A
/// </summary>
/// <param name="A">матрица СЛАУ</param>
/// <param name="b">вектор правой части</param>
/// <param name="x">вектор решения СЛАУ</param>
/// <returns>
///  1 - если сходимость была достигнута,
///  0 - превышено максимальное число итераций
/// </returns>
int SparseGradientDescent(const vector<SparseElement>& A, 
	const vector<SparseElement>& b, 
	vector<SparseElement>& x)
{

	// градиент: 2*(Ax0 - b)
	auto Gradient = [A, b](const vector<SparseElement>& x)
		{		
			vector<SparseElement> result = SparseDifference(SparseMultiply(A, x), b);
			for (auto& elem : result)
				elem.value *= 2.0;
			return result;
		};

	// A*(Ax0 - b)
	auto VectorRight = [A, b](const vector<SparseElement>& x)
		{
			vector<SparseElement> ax = SparseDifference(SparseMultiply(A, x), b);
			vector<SparseElement> result = SparseMultiply(A, ax);
			return result;
		};

	vector<SparseElement> x0;
	x0.push_back({ 1, 1, 1.0 });
	double err_norm0 = 1.0;
	vector<SparseElement> grad;

	int iter = 0;
	double err_norm = 0.0, t = 0.0;
	do
	{
		grad = Gradient(x0);
		double sc1 = SparseScalarProduct(grad, grad) * 0.25;

		double sc2 = SparseScalarProduct(grad, VectorRight(x0));
		t = sc1 / sc2;

		for (auto& elem : grad)
			elem.value *= t;

		x = SparseDifference(x0, grad);

		// норма разности между значением на текущей и предыдущей итерации
		err_norm = SparseNormv(SparseDifference(x, x0));
		if (abs(err_norm - err_norm0) <= DBL_MIN) return 1;
		if (++iter > MAX_ITERATION_NUMBER) return 0;
		x0 = x;
		err_norm0 = err_norm;
	} while (err_norm > EPS);
	return 1;
}

Расчёты производились для различных степеней заполнения матрицы. Коэффициент релаксации $\omega = 0.9$ , максимальное число итерации 60 000, $\epsilon = 1.0\cdot 10^{-15}$ для менее 200 уравнений и $\epsilon = 1.0\cdot 10^{-12}$ для 200 и более уравнений. Результаты расчётов представлены в нижеследующей таблице:

Метод решения	Степень заполнения матрицы, %	Время решения, сек	Норма невязки
Число уравнений 50
Метод вращений	37.6	6.1541	$9.85611\cdot 10^{-13}$
Метод релаксации		0.0202589	$1.98615\cdot 10^{-14}$
Метод градиентного спуска		2.20036	$5.84411\cdot 10^{-13}$
Метод вращений	65.2	10.129	$1.30145\cdot 10^{-12}$
Метод релаксации		0.0291353	$2.76372\cdot10^{-14}$
Метод градиентного спуска		3.07949	$6.58427\cdot 10^{-13}$
Метод вращений	96.4	13.4961	$1.48361\cdot 10^{-12}$
Метод релаксации		0.039939	$3.4582\cdot 10^{-14}$
Метод градиентного спуска		4.15467	$5.30127\cdot 10^{-13}$
Число уравнений 100
Метод вращений	36.8	85.0761	$1.79234\cdot 10^{-12}$
Метод релаксации		0.0823859	$5.58988\cdot 10^{-14}$
Метод градиентного спуска		6.80603	$4.51663\cdot 10^{-13}$
Метод вращений	64.6	121.906	$1.90959\cdot 10^{-12}$
Метод релаксации		0.149928	$6.52918\cdot 10^{-14}$
Метод градиентного спуска		13.2371	$4.86747\cdot 10^{-13}$
Метод вращений	84.4	138.603	$2.90557\cdot 10^{-12}$
Метод релаксации		0.125472	$4.11035\cdot 10^{-14}$
Метод градиентного спуска		13.6735	$4.49736\cdot 10^{-13}$
Число уравнений: 200
Метод релаксации	28.2	0.272274	$7.07554\cdot 10^{-14}$
Метод градиентного спуска	28.2	23.7676	$6.23567\cdot 10^{-13}$
Метод релаксации	51.4	0.304953	$5.67969\cdot 10^{-12}$
Метод градиентного спуска	51.4	38.7658	$4.87242\cdot 10^{-10}$
Метод релаксации	75.3	0.408678	$1.0519\cdot 10^{-12}$
Метод градиентного спуска	75.3	47.8197	$4.8624\cdot 10^{-10}$
Число уравнений: 500
Метод релаксации	22.7	1.13415	$2.95622\cdot 10^{-12}$
Метод градиентного спуска	22.7	164.936	$4.8934\cdot 10^{-10}$
Метод релаксации	51.1	1.92343	$6.75017\cdot 10^{-12}$
Метод градиентного спуска	51.1	438.159	$4.86602\cdot 10^{-10}$
Метод релаксации	75.1	2.58393	$9.2081\cdot 10^{-12}$
Метод градиентного спуска	75.1	594.39	$4.85809\cdot 10^{-10}$
Число уравнений: 1000
Метод релаксации	36.1	6.42248	$5.14759\cdot 10^{-12}$
Метод релаксации	64.1	9.33765	$5.40416\cdot 10^{-12}$
Метод релаксации	79.8	10.9605	$6.3317\cdot 10^{-12}$
Число уравнений: 2000
Метод релаксации	27.8	23.2207	$5.30702\cdot 10^{-12}$
Метод релаксации	64	39.8032	$8.12667\cdot 10^{-12}$
Метод релаксации	84	55.3242	$8.7574\cdot 10^{-12}$
Число уравнений: 3000
Метод релаксации	36	69.5115	$6.70979\cdot 10^{-12}$
Метод релаксации	64	99.1673	$7.42065\cdot 10^{-12}$
Метод релаксации	89	132.051	$8.5561\cdot 10^{-12}$

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

Из приведенных результатов видно, что безусловным лидером как по точности, так и по времени решения является метод релаксации. Все методы демонстрировали высокую точность решения.

Выводы

Там, где не требуется высокая точность, вполне подходят методы Гаусса, LU декомпозиции, компактной схемы исключения, при этом метод Гаусса более предпочтителен, так как требует меньше временных затрат на решение. Для достижения более высокой точности предпочтительнее использовать для решения методы QR декомпозиции или декомпозиции, демонстрирующие на один или два порядка более высокую точность, но при этом приходится расплачиваться большими затратами времени на решение. Методы вращения, релаксации и градиентного спуска не подходят для решения СЛАУ со 100% степенью заполнения матрицы в связи с их плохой сходимостью.

Для разреженных матриц метод релаксации наиболее предпочтителен как в плане затрат времени, так и точности решения.

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

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