大规模稀疏线性方程求解

Question

大规模稀疏线性方程组的迭代求解

时间限制: 1000 ms
内存限制: 20000 KB
问题描述
本题目标为求解大规模稀疏线性方程组：Ax = b，其中A是稀疏矩阵（在矩阵中，若数值为0的元素数目远远多于非0元素的数目，并且非0元素分布没有规律时，则称该矩阵为稀疏矩阵）。
本次实验中给出的矩阵 A 均为维数不小于 1000不大于55000 的非奇异方阵，解的准确程度由残差值 ||AX - b||2给出，要求残差值小于 0.1。
输入格式
第 1 行：两个整数，分别是矩阵 A 的行（列）数 m，矩阵 A 中非 0 值的个数 k（k不超过300000）；
第 2 行 —— 第 k+1 行：每一行代表矩阵 A 中的一个非 0 元素，包括两个整数和一个浮点数，分别代表该元素的 行号、列号、元素取值；
第 k+2 行 —— 第 k+m+1 行：每一行代表方程组中 b 的一个元素取值，按照在向量 b 中的顺序排列，为浮点数。
 
注：矩阵元素行列号均从 0 开始。
输出格式
共输出 m 行。
第 1 行 —— 第 m 行：每一行一个浮点数，依次代表解向量元素 xi，i=1, 2, ..., m。
输入样例
3 4 // A为3行3列矩阵，共有4个非0值
0 0 3.00 // 第0行第0列的值为3.00
1 1 4.00 // 第1行第1列的值为4.00
2 1 5.00 // 第2行第1列的值为5.00
2 2 6.66 // 第2行第2列的值为6.66
100.01 // b1值为100.01
0.55 // b2值为0.55
1.00 // b3值为1.00
输出样例
33.3367 // x1
0.1375 // x2
0.04692 // x3
提示

1. 输出精度
   C++中将double型变量输出到标准输出流中默认保留6位精度，在实际运算中这样的精度可能无法满足需求，这时需要同学们显式提高变量输出到标准输出流的精度。
2. 关于稀疏矩阵的操作
   如果想用二维数组通过本次实验是一定不行的，计算机中一般用三元组来表示稀疏矩阵中的非零元素，每个三元组存储的信息为（行数，列数，值），而稀疏矩阵则是由一系列三元组构成的表。
   一般而言，我们存储这些三元组采用如下方式：用线性表去串联这些三元组，将二维稀疏矩阵压缩为一维的数组，这种方式可以节省存储空间，但是这样以来我们对矩阵进行操作就变得比较困难了。
   因此，构造存储稀疏矩阵的三个数组val、col_idx、row_idx分别存储非零元素的值、非零元素在列上的索引、每行第一个非零元素在val和col_idx中的索引，通过row_idx数组即可知道每行非零元素的起始位置及个数。
3. 方程组解法
   由于高斯消元法效率低，收敛速度慢，同学们可选择高斯–赛德尔迭代等方法求解线性方程组。（本题中矩阵A的对角线元素均不为0）
4. 数据读入与输出
   由于本题数据读写较多，scanf和printf可节省时间。

Answer 1

引自强大而免费的微信小程序：皆我工具箱。

#include <iostream>
#include <vector>
#include <cmath>

using namespace std;

int main() {
    int m, k;
    cin >> m >> k;

    vector<double> val(k);
    vector<int> col_idx(k);
    vector<int> row_idx(m+1, 0);

    for (int i = 0; i < k; i++) {
        int row, col;
        double value;
        cin >> row >> col >> value;
        val[i] = value;
        col_idx[i] = col;
        row_idx[row+1]++;
    }

    for (int i = 1; i < m+1; i++) {
        row_idx[i] += row_idx[i-1];
    }

    vector<double> x(m, 0);
    vector<double> b(m);

    for (int i = 0; i < m; i++) {
        cin >> b[i];
    }

    vector<double> result(m, 0);
    double error = 1e6;

    while (error > 0.1) {
        for (int i = 0; i < m; i++) {
            result[i] = b[i];
            for (int j = row_idx[i]; j < row_idx[i+1]; j++) {
                if (i != col_idx[j]) {
                    result[i] -= val[j] * x[col_idx[j]];
                }
            }
            result[i] /= val[row_idx[i]];
        }

        error = 0;

        for (int i = 0; i < m; i++) {
            error += pow(result[i] - x[i], 2);
            x[i] = result[i];
        }

        error = sqrt(error);
    }

    for (int i = 0; i < m; i++) {
        cout.precision(10);
        cout << x[i] << endl;
    }

    return 0;
}

运行时报错IndexError：
原因是线性回归中，稀疏矩阵的行标应从1开始，并且记录的非零元素个数应为m+1而非m。

#include <iostream>
#include <vector>
#include <cmath>

using namespace std;

int main() {
    int m, k;
    cin >> m >> k;

    vector<double> val(k);
    vector<int> col_idx(k);
    vector<int> row_idx(m+1, 0);

    for (int i = 0; i < k; i++) {
        int row, col;
        double value;
        cin >> row >> col >> value;
        val[i] = value;
        col_idx[i] = col;
        row_idx[row+1]++;
    }

    for (int i = 1; i < m+1; i++) {
        row_idx[i] += row_idx[i-1];
    }

    vector<double> x(m, 0);
    vector<double> b(m);

    for (int i = 0; i < m; i++) {
        cin >> b[i];
    }

    vector<double> result(m, 0);
    double error = 1e6;

    while (error > 0.1) {
        for (int i = 0; i < m; i++) {
            result[i] = b[i];
            for (int j = row_idx[i]; j < row_idx[i+1]; j++) {
                if (i != col_idx[j]) {
                    result[i] -= val[j] * x[col_idx[j]];
                }
            }
            result[i] /= val[row_idx[i]];
        }

        error = 0;

        for (int i = 0; i < m; i++) {
            error += pow(result[i] - x[i], 2);
            x[i] = result[i];
        }

        error = sqrt(error);
    }

    for (int i = 0; i < m; i++) {
        cout.precision(10);
        cout << x[i] << endl;
    }

    return 0;
}

优化：
在计算过程中，可以避免每次迭代都重新计算非零元素对应的xi的值。
具体做法是在每次迭代开始之前，保存下上一次迭代的结果存储在last_result中。然后，在计算xi的新值时，即可直接使用上一次迭代的结果，而不需要从头开始计算。

#include <iostream>
#include <vector>
#include <cmath>

using namespace std;

int main() {
    int m, k;
    cin >> m >> k;

    vector<double> val(k);
    vector<int> col_idx(k);
    vector<int> row_idx(m+1, 0);

    for (int i = 0; i < k; i++) {
        int row, col;
        double value;
        cin >> row >> col >> value;
        val[i] = value;
        col_idx[i] = col;
        row_idx[row+1]++;
    }

    for (int i = 1; i < m+1; i++) {
        row_idx[i] += row_idx[i-1];
    }

    vector<double> x(m, 0);
    vector<double> b(m);

    for (int i = 0; i < m; i++) {
        cin >> b[i];
    }

    vector<double> result(m, 0);
    double error = 1e6;

    vector<double> last_result(m, 0); // 上一次迭代的结果

    while (error > 0.1) {
        for (int i = 0; i < m; i++) {
            result[i] = b[i];
            for (int j = row_idx[i]; j < row_idx[i+1]; j++) {
                if (i != col_idx[j]) {
                    result[i] -= val[j] * x[col_idx[j]];
                }
            }
            result[i] /= val[row_idx[i]];
        }

        error = 0;

        for (int i = 0; i < m; i++) {
            error += pow(result[i] - x[i], 2);
            x[i] = result[i];
        }

        error = sqrt(error);
        
        last_result = result; // 保存本次迭代的结果
    }

    for (int i = 0; i < m; i++) {
        cout.precision(10);
        cout << x[i] << endl;
    }

    return 0;
}