GSL中的稀疏线性代数

稀疏线性代数

    本章描述求解稀疏线性方程组的函数。本库提供了直接操作gsl_spmatrix和gsl_vector对象的线性代数函数。

本章描述的函数声明在头文件gsl_splinalgh中。

45.1 概述

    本章主要讨论如下线性系统的解，

                                                       Ax=b

    其中A是一个一般矩形n乘n非奇异稀疏矩阵，x是一个未知的n乘1向量，b是一个给定的n乘1右边的向量。求解稀疏线性系统的方法有很多，大体上可分为直接或迭代两类。直接方法包括LU分解和QR分解，迭代方法从向量x的初始猜测开始，通过迭代更新猜测直到收敛。GSL目前不提供任何直接稀疏解算器。

45.2 稀疏迭代求解器

45.2.1 概述

    许多实际的求解n乘n大型稀疏线性系统的迭代方法，涉及到将x的近似解投影到Rn的子空间上。如果我们定义一个m维的子空间K为解x的近似子空间，那么必须施加m个约束来决定下一个近似。这些m个约束定义了另一个m维子空间，用L表示。为了减少生成下一个近似解向量所需的计算工作量，子空间的维数m通常被选择比n小得多。目前存在的许多迭代算法的主要区别在于对K和L的选择。

45.2.2 稀疏迭代求解器的类型

稀疏线性代数库提供了以下类型的迭代求解器：

gsl_splinalg_itersolve_type

gsl_splinalg_itersolve_gmres

这就是广义最小残差法(GMRES)。这是一种使用K=Km，L=AKm的投影方法其中Km是第m个Krylov子空间

     r0=b-Ax0 是初始猜想x0 的残差向量。如果m被设为n，那么Krylov子空间是Rn，GMRES将提供精确解x。然而，该方法的目标是使用一个小得多的子空间Km得到一个非常好的近似x。默认情况下，GMRES方法选择m=MIN(n, 10)，但是用户可以为m指定一个不同的值。GMRES存储需求随着O(n(m+1))的增长而增长，而失败的数量随着的增长而增长。

 在下面的函数gsl_splinalg_itersolve_iterate()中，一次GMRES迭代定义为将近似解向量投影到每个Krylov子空间K1,…,Km上，所以m可以通过“重新启动”该方法并多次调用gsl_splinalg_itersolve_iterate()来保持更小，并为每个新调用提供更新的近似x。如果方法没有充分收敛，用户可以尝试增加参数m。

GMRES被认为是一种鲁棒的通用迭代求解器，但如果矩阵不是正定的，该方法也会停滞不前，直到最后一个投影到子空间Kn=Rn上时才会减少残差。在这些情况下，预处理线性系统可以帮助，但GSL目前不提供任何预处理。

45.2.3 迭代稀疏线性系统

    提供以下函数为稀疏线性求解器分配存储空间，并将系统迭代到一个解决方案。

gsl_splinalg_itersolve * gsl_splinalg_itersolve_alloc(const gsl_splinalg_itersolve_type * T, const size_t n, const size_t m)

    本函数为n乘n稀疏矩阵系统的迭代解分配了一个工作空间。迭代求解器类型由T指定，参数m指定解候选子空间Km的大小。维度m可以设置为0，在这种情况下使用合理的默认值。

void gsl_splinalg_itersolve_free(gsl_splinalg_itersolve * w)

    本函数释放与工作空间w相关的内存。