SOR法の加速パラメータωを最適化してみた

　SOR法とは連立一次方程式を反復的に解く方法の一つです。SOR法はハイパーパラメータとして加速パラメータωというものを必要とします。教科書的にはωは予め定数として与えておくものらしいですが、今回は誤差が最小になるようなωを各ステップごとに求めながら計算していく方法を実装してみました。

方法

　SOR法のアルゴリズムは簡単なので説明は控えます。
SOR法 - Wikipedia
今回の実装では（式が簡単になるので）ガウス・ザイデル法ではなくヤコビ法を元にして実装しました。よって漸化式は

となります。ここで、ωを変数としたままを求めると各要素がωの一次式になるので、誤差の二乗ノルムはωの二次式で表されます。その二次式の最小値をとしての各要素に代入すればよいというわけです。

結果

　結果として、ωを決め打ちするSOR法と比べて反復回数は少なくなりました(当たり前ですが)。ただ、を求める際に行列ベクトル積を2回計算する必要があるのですが、そのデメリットを補えるほど加速されるかは場合によります。私が試したところ、非対角成分を小さな負の値にしたときにが2より遥かに大きくなるステップが見つかり、その場合はこの方法のほうが結果として速くなりました。

#include <stdio.h>
#include <conio.h>
#include <assert.h>
#include <random>

void init_matrix(double* src, const int dim, const bool diag_domi, const double lowerbound, const double upperbound)
{
	//srcは[dim*dim]の配列であると仮定する。これをdim次正方行列とみなして、
	//対角成分が全て正であるような対称行列になるように適当な値を代入して返す。
	//diag_domiがtrueなら狭義対角優位にする。

	assert(lowerbound < upperbound);

	std::mt19937 gen(1729);
	std::uniform_real_distribution<> ra(lowerbound, upperbound);

	for (int i = 0; i < dim; i++)
	{
		for (int j = 0; j < i; j++)src[i * dim + j] = src[j * dim + i];//対称行列なので
		for (int j = i; j < dim; j++)src[i * dim + j] = ra(gen);
	}

	//対角成分の処理。対角優位にしない場合でも、対角成分で割る処理があるので非ゼロの値にしておく。
	for (int i = 0; i < dim; i++)
	{
		double sum = 0.0;
		for (int j = 0; j < dim; j++)sum += abs(src[i * dim + j]);
		src[i * dim + i] += diag_domi ? sum + ((double)dim) * 0.01 : abs(src[i * dim + i]) + 1.5;
	}
}
void init_vector(double* src, const int dim)
{
	std::mt19937 gen(1729);
	std::uniform_real_distribution<> ra(0.0, 10.0);

	for (int i = 0; i < dim; i++)src[i] = ra(gen);
}
void mul(const double* mat, const double* vec, const int dim, double* dest)
{
	for (int i = 0; i < dim; i++)
	{
		dest[i] = 0.0;
		for (int j = 0; j < dim; j++)dest[i] += mat[i * dim + j] * vec[j];
	}
}
void sub(const double* vec1, const double* vec2, const int dim, double* dest)
{
	for (int i = 0; i < dim; i++)dest[i] = vec1[i] - vec2[i];
}
void mad(const double* vec1, const double p, const double* vec2, const int dim, double* dest)
{
	for (int i = 0; i < dim; i++)dest[i] = vec1[i] * p + vec2[i];
}
void Jacobi(const int dim, const bool diag_domi, const double epsilon, const double lowerbound, const double upperbound)
{
	assert(0 < dim);

	double* A = (double*)malloc(sizeof(double) * dim * dim);
	double* x = (double*)malloc(sizeof(double) * dim);
	double* b = (double*)malloc(sizeof(double) * dim);

	init_matrix(A, dim, diag_domi, lowerbound, upperbound);
	init_vector(b, dim);
	for (int i = 0; i < dim; i++)x[i] = b[i];

	double* x1 = (double*)malloc(sizeof(double) * dim);
	double* v1 = (double*)malloc(sizeof(double) * dim);
	double* v2 = (double*)malloc(sizeof(double) * dim);

	//連立一次方程式Ax=bをJacobi法で解く。

	for (int loop = 0; loop < 1000; loop++)
	{
		//この時点で、近似解はxに格納されている。

		//ヤコビ法を適用する。
		for (int i = 0; i < dim; i++)
		{
			x1[i] = b[i];
			for (int j = 0; j < dim; j++)if (i != j)x1[i] -= A[i * dim + j] * x[j];
			x1[i] /= A[i * dim + i];
		}
		for (int i = 0; i < dim; i++)x[i] = x1[i];

		//誤差ノルム||Ax-b||を求めてerrorに格納する。
		mul(A, x, dim, v1);
		sub(v1, b, dim, v2);
		double error = 0.0;
		for (int i = 0; i < dim; i++)error += v2[i] * v2[i];

		printf("反復%d回目  Jacobi  error=%f\n", loop, error);
		_getch();

		//新しい近似解が目標の精度に達しているなら終了。
		if (error < epsilon)break;
	}

	free(A);
	free(x);
	free(b);
	free(x1);
	free(v1);
	free(v2);
}
void SOR(const int dim, const bool diag_domi, const double epsilon, const double lowerbound, const double upperbound, const double omega = -1.0)
{
	assert(0 < dim);

	double* A = (double*)malloc(sizeof(double) * dim * dim);
	double* x = (double*)malloc(sizeof(double) * dim);
	double* b = (double*)malloc(sizeof(double) * dim);

	init_matrix(A, dim, diag_domi, lowerbound, upperbound);
	init_vector(b, dim);
	for (int i = 0; i < dim; i++)x[i] = b[i];

	double* x1 = (double*)malloc(sizeof(double) * dim);
	double* v1 = (double*)malloc(sizeof(double) * dim);
	double* v2 = (double*)malloc(sizeof(double) * dim);
	double* v3 = (double*)malloc(sizeof(double) * dim);

	//連立一次方程式Ax=bをSOR法で解く。

	for (int loop = 0; loop < 1000; loop++)
	{
		//この時点で、近似解はxに格納されている。

		//加速パラメータωとしてSOR法を適用後の新しい近似解ベクトルx'の値を求める。x'=x1*ω+xとする。
		for (int i = 0; i < dim; i++)
		{
			x1[i] = b[i];
			for (int j = 0; j < dim; j++)x1[i] -= A[i * dim + j] * x[j];
			x1[i] /= A[i * dim + i];
		}

		double error = 0.0;
		double omega_using = omega;

		if (0.0 < omega && omega < 2.0)//引数omegaに有効な値が指定されているならそれを使う。
		{
			for (int i = 0; i < dim; i++)x[i] += omega * x1[i];

			//誤差ノルム||Ax-b||を求めてerrorに格納する。
			mul(A, x, dim, v1);
			sub(v1, b, dim, v2);
			for (int i = 0; i < dim; i++)error += v2[i] * v2[i];
		}
		else//最適なomegaを求める。
		{
			//ベクトルx'の各成分はωの1次式であるが、ここで||Ax-b||を求める。
			//||Ax-b||はスカラーでありωの2次式で表される。||Ax-b||=p2*ω^2+p1*ω+p0とする。
			mul(A, x1, dim, v1);
			mul(A, x, dim, v2);
			sub(v2, b, dim, v3);
			double p2 = 0.0, p1 = 0.0, p0 = 0.0;
			for (int i = 0; i < dim; i++)
			{
				p2 += v1[i] * v1[i];
				p1 += v1[i] * v3[i] * 2.0;
				p0 += v3[i] * v3[i];
			}

			//||Ax-b||が最小となるようなωの値を求める。
			omega_using = -p1 / (2.0 * p2);

			//そのωのときの||Ax-b||を求める。
			error = p0 - ((p1 * p1) / (4.0 * p2));

			//そのωを使って新しい近似解を求め、xに格納する。
			mad(x1, 1, x, dim, v1);
			for (int i = 0; i < dim; i++)x[i] = v1[i];

		}

		printf("反復%d回目  ω=%f  error=%f\n", loop, omega_using, error);
		_getch();

		//新しい近似解が目標の精度に達しているなら終了。
		if (error < epsilon)break;
	}

	free(A);
	free(x);
	free(b);
	free(x1);
	free(v1);
	free(v2);
	free(v3);
}
int main()
{
	printf("JacobiとSOR(ω=1.0)が同じ結果になることを確認する。\n");
	Jacobi(100, true, 0.00001, -1.0, 1.0);
	SOR(100, true, 0.00001, -1.0, 1.0, 1.0);

	printf("ωを最適化させれば速く収束することを確認する。\n");
	SOR(100, true, 0.00001, -1.0, 1.0);

	printf("ω最適化がうまく働く例。\n");
	SOR(100, true, 0.00001, -1.0, 0.0);

	printf("おわり");
	_getch();
	return 0;
}