下一节:迭代线性求解器简述 上一级:直接线性求解器简述 上一节:Direct Solvers for Sparse

结构化矩阵的直接求解器

在第§10.2.3节中描述的涉及结构化矩阵的线性方程组具有一个吸引人的特性,即可以在O(n^2)时间内求解,而不是O(n^3)时间。

使得这些系统能够快速(在O(n^2)时间内)求解的共同特性是它们具有低位移秩[257]。

位移秩可以定义如下:对于给定的矩阵AF以及矩阵T,我们定义运算符

\Delta_{A,F}(T)=A\cdot T-T\cdot F
并记\alpha={\rm rank}(\Delta_{A,F}(T))。我们称\alphaT相对于运算符\Delta_{A,F}的“位移秩”。如果\alpha很小(通常为O(1)),我们称矩阵T具有“低位移秩”。由于 {\rm rank}(\Delta_{A,F}(T))=\alpha,存在 n\times \alpha矩阵GB,使得 \Delta_{A,F}(T)=G\cdot B^T。矩阵GB被称为T生成元。我们将具有低位移秩的矩阵称为具有位移结构

以Cauchy矩阵为例 C = (C_{ij})=(1/(x_i-y_j)),我们有

\Delta_{\mathrm{diag}(x),\mathrm{diag}(y)}(C)=\mathrm{diag}(x) \cdot C-C \cdot \mathrm{diag}(y)=(1,1,\ldots,1)^T\cdot (1,1,\ldots,1).
因此,Cauchy矩阵的位移秩为1,生成元为 G=B=(1,1,\ldots,1)^T

对于具有低位移秩\alpha(不一定等于1)的矩阵,相对于运算符 \Delta_{{\rm diag}(x),{\rm diag}(y)},称为Cauchy-like矩阵。类似地,可以定义Vandermonde-like矩阵、Toeplitz-like矩阵等。这些高位移秩系统的求解时间为O(\alpha n^2)

快速算法利用位移方程来生成矩阵F的LU分解。这些算法在矩阵的生成元上工作,而不是矩阵本身,这使得它们能够更快地运行。

块Toeplitz矩阵和类型为T^TT的矩阵(其中T是Toeplitz矩阵)也具有低位移秩[257]。

如果AB具有位移结构,那么类型为 (A+\sigma B)x=b的系统可以使用快速低位移秩方法求解,前提是A+\sigma B也具有位移结构(相对于可能不同的位移运算符)。例如,如果A是Hankel矩阵且B=I,则A+\sigma B是Toeplitz-plus-Hankel矩阵,也具有位移结构。如果A+\sigma B不具有位移结构,那么只能使用零移位\sigma=0或仅需快速矩阵向量乘积的快速迭代方法[257]。

一些方法对矩阵施加了额外的限制,要求其对称(Toeplitz)或正定。这可能导致对移位选择的额外限制。

在使用结构化矩阵的快速算法时应特别注意,因为速度通常以精度为代价。一些位移秩算法通过在生成元上工作来模拟高斯消元法(带部分、完全或无枢轴)[193]。即便如此,这些算法由于可能的“生成元增长”[257],不一定具有与高斯消元法相同的数值特性,这种情况虽不常见,但即使在非常良态的矩阵中也可能出现。

有许多方法可以稳定快速算法,包括生成元增长问题[257, p. 111],尽管偶尔会出现不稳定性,但这些方法非常吸引人。

快速算法在文献中有很好的描述,但可靠的软件库尚不可用。人们还应注意O(n^2)符号中隐藏的常数,这些常数使得快速方法仅在n足够大(通常至少几百,取决于结构)时比传统的O(n^3)方法更快。传统的O(n^3)算法通常经过优化以使用BLAS 3,并有效利用计算机的内存层次结构,以使代码接近处理器的全速运行(参见第§10.2.1节)。快速算法通常只能使用BLAS 2,并且可能难以或无法优化快速算法以有效利用计算机的缓存和快速内存。这意味着,即使快速O(n^2)算法和慢速O(n^3)算法执行相同数量的浮点运算,由于这些优化问题,“慢”算法很可能更快。

我们还应提及用于求解Vandermonde和三项Vandermonde系统的Bjorck-Pereyra型算法[51,230]。这些O(n^2)方法具有显著的数值特性。在Vandermonde矩阵中节点顺序和右侧符号模式的某些额外限制下,可以完全达到机器精度求解这些系统。



下一节:迭代线性求解器简述 上一级:直接线性求解器简述 上一节:Direct Solvers for Sparse
Susan Blackford 2000-11-20