高斯法求解线性方程组

最新推荐文章于 2023-03-07 08:14:05 发布

我真的不是cjc

最新推荐文章于 2023-03-07 08:14:05 发布

阅读量1.6k

点赞数 2

分类专栏： # 算法笔记文章标签：算法 c++ 微积分高数线性代数

本文链接：https://blog.csdn.net/qq_30142297/article/details/124033621

版权

算法笔记专栏收录该内容

14 篇文章 0 订阅

订阅专栏

求解线性方程组的高斯方法（Gauss method for solving system of linear equations）

给出一个有m个未知数的n个线性代数方程组（SLAE, System of Linear Algebraic Equations）。你需要解决的是：确定它是否无解，是否仅有一个解或有无限多的解。如果它至少有一个解，请找出其中任何一个。

$\begin{aligned} a_{11} x_1 + a_{12} x_2 + &\dots + a_{1m} x_m = b_1\\ a_{21} x_1 + a_{22} x_2 + &\dots + a_{2m} x_m = b_2\\ &\vdots \\ a_{n1} x_1 + a_{n2} x_2 + &\dots + a_{nm} x_m = b_n \end{aligned}$

其中系数 $a_{ij}$ （ $i$ 从 $1$ 到 $n$ ， $j$ 从 $1$ 到 $m$ ）和 $b_i$ （ $i$ 从 $1$ 到 $n$ ）是已知的，变量 $x_i$ （ $i$ 从 $1$ 到 $m$ ）是未知的。

这个问题也有一个简单的矩阵表示：

$A x = b$

其中 $A$ 是系数 $a_{ij}$ 的大小为 $n \times m$ 的矩阵， $b$ 是大小为 $n$ 的列向量。

值得注意的是，本文介绍的方法也可以用来解以任意数 $p$ 为模的方程，即：

$\begin{aligned} a_{11} x_1 + a_{12} x_2 + &\dots + a_{1m} x_m \equiv b_1 \pmod p \\ a_{21} x_1 + a_{22} x_2 + &\dots + a_{2m} x_m \equiv b_2 \pmod p \\ &\vdots \\ a_{n1} x_1 + a_{n2} x_2 + &\dots + a_{nm} x_m \equiv b_n \pmod p \end{aligned}$

高斯（Gauss）

严格来说，下面描述的方法应该称为“Gauss-Jordan”，或 Gauss-Jordan 消元法，因为它是 Jordan 在 1887 年描述的 Gauss 方法的一种变体。

概述

该算法是对每个方程中的变量进行 “顺序消除”，直到每个方程只剩下一个变量。如果 $n = m$ ，你可以认为它是将矩阵 $A$ 转化为单位矩阵，并在这种明显的情况下解决方程，其解是唯一的且等于系数 $b_i$ 。

高斯消元法基于两个简单的变换：

可以交换两个方程
任何方程都可以用该行（具有非零系数）和其他一些行（具有任意系数）的线性组合代替

在第一步， Gauss-Jordan 算法将第一行除以 $a_{11}$ 。然后，该算法将第一行与其余各行相加，使第一列中的系数成为全部为零。为此，在第 $i$ 行，我们必须把第一行乘以 $a_{i1}$ 。注意，这个操作也必须在向量 $b$ 上执行。在某种意义上，它的表现就像向量 $b$ 是矩阵 $A$ 的第 $m$ 列一样。

因此，第一步之后，矩阵 $A$ 的第一列将由第一行的 $1$ 和其他行的 $0$ 组成。

同样地，我们执行算法的第二步，我们考虑第二行的第二列。首先，该行被 $a_{22}$ 除以，然后从其他行中减去，这样所有第二列就变成了 $0$ （除了第二行）。

我们对矩阵 $A$ 的所有列继续这个过程。如果 $n = m$ ，那么 $A$ 将成为单位矩阵。

寻找基准元素（Search for the pivoting element）

所描述的方案遗漏了许多细节。在第 $i$ 步，如果 $a_{ii}$ 为零，我们不能直接应用描述的方法。相反，我们必须首先 选择一个基准行 ：找到矩阵中第 $i$ 列非零的一行，然后交换这两行。

请注意，这里我们交换行而不是列。这是因为如果您交换列，那么当您找到解决方案时，您必须记住交换回正确的位置。因此，交换行要容易得多。

在许多实现中，当 $a_{ii} \neq 0$ 时，你可以看到仍然用一些启发式行来交换第 $i$ 行，使用一些启发式的方法，如选择绝对值最大的 $a_{ji}$ 的启发式行。这种启发式方法是用来在后面的步骤中减少矩阵的取值范围。如果没有这个启发式方法，即使是大小为 $20$ 的矩阵，误差也会过大，并可能导致 C++ 的浮点数据类型溢出。

降阶案例（Degenerate cases）

在 $m = n$ 和系统非降阶的情况下（即它有非零行列式，并有唯一的解决方案），上述算法将把 $A$ 转化为单位矩阵。

现在我们考虑 一般情况 ，即 $n$ 和 $m$ 不一定相等，而且方程组是可能降阶的。在这些情况下，第 $i$ 步中的基准元素可能找不到。这意味着在第 $i$ 列，从当前行开始，全部包含零。在这种情况下，要么变量 $x_i$ 没有可能的值（意味着 SLAE 没有解决方案），要么 $x_i$ 是一个独立变量，可以取任意值。在实现 Gauss-Jordan 算法时，你应该继续后续变量的工作，只是跳过第 $i$ 列（这相当于删除矩阵的第 $i$ 列）。

因此，可以发现过程中的一些变量是独立的。当变量的数量 $m$ 大于方程的数量 $n$ 时，那么至少可以找到 $m - n$ 个独立变量。

一般来说，如果你发现至少有一个自变量，它可以取任何任意的值，而其他（因）变量则通过它来表达。这意味着，当我们在实数领域工作时，方程组可能有无无穷多的解。但你应该记住，当有独立变量时，SLAE可能根本没有解。当剩下的未处理的方程至少有一个非零常数项时，就会发生这种情况。你可以通过给所有独立变量赋零来检查这一点，计算其他变量，然后插入到原来的SLAE中，检查它们是否满足。

实现

以下是 Gauss-Jordan 的实现。选择基准行是通过启发式完成的：在当前列中选择最大值。

函数 Gauss 的输入是矩阵组 $a$ ，该矩阵的最后一列是向量 $b$ 。

该函数返回系统的解的数量 $1,\textrm{or } \infty)$ 。如果至少有一个解存在，那么它将被返回到向量 $A$ 中。

const double EPS = 1e-9;
const int INF = 2; // it doesn't actually have to be infinity or a big number

int gauss (vector < vector<double> > a, vector<double> & ans) {
	int n = (int) a.size();
	int m = (int) a[0].size() - 1;

	vector<int> where (m, -1);
	for (int col=0, row=0; col<m && row<n; ++col) {
		int sel = row;
		for (int i=row; i<n; ++i)
			if (abs (a[i][col]) > abs (a[sel][col]))
				sel = i;
		if (abs (a[sel][col]) < EPS)
			continue;
		for (int i=col; i<=m; ++i)
			swap (a[sel][i], a[row][i]);
		where[col] = row;

		for (int i=0; i<n; ++i)
			if (i != row) {
				double c = a[i][col] / a[row][col];
				for (int j=col; j<=m; ++j)
					a[i][j] -= a[row][j] * c;
			}
		++row;
	}

	ans.assign (m, 0);
	for (int i=0; i<m; ++i)
		if (where[i] != -1)
			ans[i] = a[where[i]][m] / a[where[i]][i];
	for (int i=0; i<n; ++i) {
		double sum = 0;
		for (int j=0; j<m; ++j)
			sum += ans[j] * a[i][j];
		if (abs (sum - a[i][m]) > EPS)
			return 0;
	}

	for (int i=0; i<m; ++i)
		if (where[i] == -1)
			return INF;
	return 1;
}

实现说明：

该函数使用两个指针——当前列 $c o l$ 和当前行 $r o w$ 。
对于每个变量 $x_i$ ，值 $w h e r e (i)$ 是该列不为零的行。需要这个向量是因为有些变量可以是独立的。
在这个实现中，当前的第 $i$ 行并没有像上面描述的那样除以 $a_{ii}$ ，所以最后的矩阵并不是单位矩阵（尽管显然除掉第 $i$ 行可以帮助减少错误）。
找到解决方案后，将其重新插入矩阵以检查系统是否至少有一组解。如果测试解决方案成功，则函数返回 $1$ 或 $i n f$ ，取决于是否存在至少一个自变量。

时间复杂度

现在我们应该估计一下这个算法的复杂性。该算法由 $m$ 个阶段组成，在每个阶段：

搜索并重新排列基准行。使用上面提到的启发式方法时需要 $O (n + m)$
如果找到当前列中的基准元素，那么我们必须将此等式添加到所有其他等式中，这需要 $O (n m)$ 的时间

因此，算法的最终复杂度为 $m)\cdot nm)$ 。如果 $n = m$ , 复杂度即为 $O(n^3)$ 。

请注意，当SLAE不在实数上，而是在模二中时，那么系统的求解速度就会快得多，这将在下面说明。

算法加速

通过将算法分为两个阶段：正向和反向，之前的实现可以加快2倍。

前进阶段。与之前的实现类似，但当前行只加到它之后的行上。因此，我们得到一个三角矩阵，而不是对角矩阵。
反向阶段。当矩阵为三角矩阵时，我们首先计算出最后一个变量的值。然后将此值插入，找到下一个变量的值。然后插入这两个值以查找下一个变量…

反相只需要 $O (n m)$ ，这比正向阶段快得多。在前向阶段，我们将操作数量减少了一半，从而减少了实现的运行时间。

求解模块化 SLAE

为了解决某些模块中的 SLAE，我们仍然可以使用所描述的算法。但是，如果模块等于 $2$ ，我们可以使用按位运算和 C++ 位集数据类型更有效地执行 Gauss-Jordan 消除：

int gauss (vector < bitset<N> > a, int n, int m, bitset<N> & ans) {
    vector<int> where (m, -1);
    for (int col=0, row=0; col<m && row<n; ++col) {
        for (int i=row; i<n; ++i)
            if (a[i][col]) {
                swap (a[i], a[row]);
                break;
            }
        if (! a[row][col])
            continue;
        where[col] = row;

        for (int i=0; i<n; ++i)
            if (i != row && a[i][col])
                a[i] ^= a[row];
        ++row;
    }
        // The rest of implementation is the same as above
}