BasicConceptsNLA.muse

#title 数值线性代数基础

<contents>

线性代数在现代数学与科学计算中起着十分基础的作用.按照计算的目的与用途,可以分为数值线性代数与精确线性代数.

数值线性代数服务于工程计算,运用数值方法主要求解三大类问题:线性方程组,最小二乘问题与特征值问题.所谓数值方法,是指由于工程计算中多不需要严格解,而只要依据计算的需要,实施精确到一定程度的快速计算即可,因而计算结果往往并非严格解,常常有一定的误差,称为数值解.由于这一特点,在进行算法设计时,往往不需考虑算法能否"真正地"严格求解,即使对一些有严格算法的问题 (如线性方程组)也可以考虑通过迭代法"逼近"严格解到一定精度.而在评价一个算法的优劣时,最重要的则是如下两条:

 1. 这一算法是否稳定?即,能否在全部或绝大部分情况下得到我们所需要的精度?与数值线性代数有关的算法稳定性分析的基本概念,我们会在"数值型算法分析基础"一节中简单介绍.
 2. 算法效率如何?即,能否在较短时间内计算出我们所需要的结果?在数值算法分析中,常用flop (floating point operation,浮点运算),即完成该算法所需进行的浮点运算次数来表征算法效率.由于通常计算机进行浮点数加减法的速度要快于浮点数乘除法,故在数值线性代数部分中我们以flop只表示浮点除乘除法的次数.例如,矩阵阶数为$n$时,Gauss消元法求解线性方程组需要$O(n^3)$次浮点数乘法(flop),其中$O$表示忽略了一个渐近的常数因子,即在$n$趋于无穷大时,flop数与$n^3$同阶.

本章将介绍数值线性代数的如下基本问题:

 1. 一般线性代数与矩阵分析的一些基本概念.
 2. 数值型算法分析基础.我们将讨论数值运算的误差与稳定性的基本概念.
 3. 矩阵乘法的快速算法.由于矩阵乘法构成几乎所有线性代数算法的基础,因此对它进行加速是很重要的.

* 线性代数与矩阵分析的基本概念

我们假定读者已经熟悉矩阵的概念与基本运算,以及线性空间的基本概念,这些内容都可以参考一般线性代数的教材,如<cite>ZhaXu04</cite>.这里我们将一些今后常用的概念和记号交代一下.

** 矩阵与线性空间

在数值分析中,我们总假定线性空间的基域是实数域$\mathbb{R}$或复数域$\mathbb{C}$.在如下的介绍中,我们如无特殊说明,都只在$\mathbb{C}$上考虑,但所有定义都容易推广到$\mathbb{R}$上.

我们通常将矩阵用大写拉丁字母或希腊字母表示,向量则以小写字母表示.如$A\in \mathbb{C}^{m\times n}$,$x\in \mathbb{C}^n$.用$A^T$表示矩阵$A\in \mathbb{C}^{m\times n}$的转置,即$(A^T)_{ji}=A_{ij},1\leqslant i \leqslant m,1\leqslant j\leqslant n.$用$A^H$表示复矩阵$A\in \mathbb{C}^{m\times n}$的共轭转置 (Hermite转置)，即$(A^H)_{ji}=\bar{A}_{ij},1\leqslant i \leqslant m,1\leqslant j\leqslant n.$

设$V$是一个线性空间,$S$是$V$中的一个向量组,$\alpha\in V$称为$S$的线性组合 (或由$S$线性表出,线性生成)是指$\alpha$可由$S$中的某有限多个向量线性表出.若$V$的子集合$W$在$V$中加法和数乘运算下是$F$上的线性空间,则称$W$是$V$的子空间.对$V$的子集合$S$,由$S$线性生成的子空间 (即$S$中任意有限个向量线性组合全体)称为由$S$生成 (或张成)的子空间,记为$\mathrm{span} (S)$.特别地,当$S=\{\alpha_1,\cdots,\alpha_s\}$为有限集时,$$\mathrm{span} (S)=\mathrm{span}\{\alpha_1,\cdots,\alpha_s\}=\mathbb{C}\alpha_1+\cdots+\mathbb{C}\alpha_s.$$

设矩阵$A\in \mathbb{C}^{m\times n}$,$A$的列空间是指由$A$的所有列向量张成的空间,记为$\mathrm{span} (A)$,它是$\mathbb{C}^m$中所有能表示为$Ax$形式的向量的集合,因此也称为矩阵$A$的值域 (range).$A$的列空间的维数称为矩阵$A$的列秩,记作$\mathrm{rank}_c (A)$. 与此相对应,$A$的行空间是指$A$的所有行向量张成的空间,记为$\mathrm{span}_r (A)$,其维数称为$A$的行秩,记作$\mathrm{rank}_r (A)$.但可证明,对于任意矩阵$A$,其行秩与列秩总是相等,通常称为矩阵$A$的秩,记作$\mathrm{span} (A)$. 矩阵$A$的零空间 (nullspace)是指满足$Ax=0$的向量$x\in \mathbb{C}^n$的集合,记作$\mathrm{null} (A)$或$\mathrm{ker} (A)$.

有如下定理,它反映了线性映射的某种"守恒"性质:

<theorem label="Thm:Dim" name="维数定理">
设$A\in \mathbb{C}^{m\times n}$,则$$n=\mathrm{dim} (\mathrm{null} (A))+rank (A).$$
</theorem>

通常我们研究的线性空间具有内积空间的结构.定义向量$x,y\in \mathbb{C}^n$的内积为$$(x,y)\equiv x^Hy\in \mathbb{C}.$$可以看出,对任意非零向量$x$,$(x,x)=x^Hx>0$,称$||x||\equiv\sqrt {x^Hx}$为$x$的范数 (norm).有时也称$$(x,Ay)\equiv x^HAy$$为广义内积,其中$A\in \mathbb{C}^{n\times n}$为Hermite矩阵,即$A^H=A$.若对于任意非零向量$x\in \mathbb{C}^n,(x,Ax)>0$(或$ (x,Ax)\leqslant 0$),则称$A$为正定矩阵 (或半正定矩阵).对于实向量,只要把内积定义中的共轭转置改为转置即可.特别地,若$(x,y)=0$,则称$x$与$y$正交 (垂直),这与解析几何中向量垂直的概念是一致的.若$V$的两个子空间$U$和$W$满足对任意$u\in U,w\in W$,$u$与$w$正交且$U+W=V$,则$U$和$W$互称为正交补,记为$U=W^\perp$.

对于矩阵$A\in \mathbb{C}^{n\times n}$,若$\lambda\in \mathbb{C}$及非零向量$\alpha\in \mathbb{C}^n$使得$$A\alpha=\lambda a,$$则称$\lambda$是矩阵$A$的特征值或特征根,$\alpha$称为$A$的属于$\lambda$的特征向量.我们注意到,$A\alpha=\lambda\alpha$相当于$(\lambda I-A)\alpha=0$有非零解$\alpha$,故$\lambda$是$A$的特征根当且仅当$$\det (\lambda I-A)=0,$$其中,多项式$f (\lambda)=\det (\lambda I-A)$称为矩阵$A$的特征多项式.计算矩阵的特征值与特征多项式是数值线性代数中的一个重要的问题,我们将在"特征值问题"两章中讨论.

对于一个任意的矩阵$A$,我们常可以通过变换矩阵$P$和$Q$使$PAQ$成为"特定的"(其精确定义以后给出) 形式,称为矩阵的标准形 (canonical form),矩阵$P$和$Q$称为变换矩阵.若要求$P$和$Q$是可逆矩阵,则称之为相抵变换.特别地,如果$P$和$Q$互为逆矩阵,则这种变换称为相似变换;如果$P$和$Q$互为转置矩阵,则这种变换称为相合变换.如果复矩阵$P$和$Q$满足$P^HP=I$以及$Q^HQ=I$,则$P$和$Q$称为酉矩阵 (unitary matrix),相应的变换也称为酉变换;与之相对应实数域上的矩阵$P$和$Q$如果满足$P^TP=I$以及$Q^TQ=I$,则称为正交矩阵 (orthogonal matrix),相应变换称为正交变换.容易看到,酉矩阵和正交矩阵的每一列 (行) 都是单位向量,且不同的列 (行) 之间相互正交.这几种变换是我们今后常用的变换,由它们可以给出矩阵的多种标准形式 (或分解),是多种算法的基础.

为了接下来讨论方便起见,我们先给出所谓奇异值分解 (singular value decomposition,SVD)的概念.

<theorem label="SVD" name="奇异值分解,或矩阵的酉相抵标准形">
设$n\times m$复矩阵$A$的秩为$r$,则有酉方阵$U_1$和$U_2$使
<latex>
\begin{equation*}
  A=U_1
  \begin{bmatrix}
    \lambda_1& & & & & \\
    & \ddots& & & & \\
    & & \lambda_r& & & \\
    & & & 0& & \\
    & & & & \ddots& \\
    & & & & & 0
  \end{bmatrix}U_2，
\end{equation*}
</latex>其中$\lambda_1\geqslant \lambda_2\geqslant\cdots\geqslant\lambda_r>0$,而$\lambda_1^2,\cdots,\lambda_r^2$为$A^HA$的所有非零特征根.
</theorem>

以上定理的证明可参考<cite>ZhaXu04</cite>.可对于实数域上的矩阵证明类似的奇异值分解 (或称为正交相抵标准形).

** 向量范数与矩阵范数

以上我们已经给出了向量范数的概念,我们进一步给出它的基本性质,其证明都很容易:
<latex>
\begin{align}
   &||x||\geqslant 0,&& x\in \mathbb{C}^n,& (||x||=0 \text{~iff~}x=0.)\notag\\ @@eq:norm
  &||x+y||\leqslant ||x||+||y||,&&x,y\in \mathbb{C}^n.&\\
  &||\alpha x||=|\alpha | ||x||,&&\alpha\in \mathbb{C},x\in \mathbb{C}^n.&\notag
\end{align}
</latex>事实上,以上性质是一般意义下的范数所具有的性质.据此,若$A$是正定Hermite阵,则广义内积$(x,Ay)$同样"诱导"出向量的范数$$||x||_A\equiv (x,Ax),$$脚标$A$用于标识不同的范数定义.除内积诱导的范数外,还有一类重要的$p-$范数经常使用,其定义如下:$$||x||_p= (|x_1|^p+\cdots+|x_n|^p)^\frac{1}{p},~p>1.$$ $p=1,2,\infty$的情形特别重要:
<latex>
\begin{align*}
&||x||_1 =  |x_1|+\cdots+|x_n|,\\
&||x||_2 = (|x_1|^2+\cdots+|x_n|^2)^\frac{1}{2},\\
&||x||_\infty =\max\limits_{1\leqslant i\leqslant n}|x_i|.
\end{align*}
</latex>由H\"older不等式$$|x^Hy|\leqslant||x||_p||y||_q,~\frac{1}{p}+\frac{1}{q}=1$$可以证明,以上定义确实满足范数定义要求的性质##eq:norm .

若把复数域上$n\times m$维矩阵全体看做一个向量空间,也可以有范数的多种定义.最常用的是以下两种:

 1. Frobenius范数:将$n\times m$维矩阵全体看做$m\times n$向量空间,定义2-范数为$$||A||_F= (\sum\limits_{ij}|A_{ij}|^2)^{\frac{1}{2}},$$这称为Frobenius范数.
 2. 算子范数:由$m$维和$n$维向量空间的$p-$范数诱导的算子范数定义如下:$$||A||_p=\sup\limits_{x\neq 0}\frac{||Ax||_p}{||x||_p},$$由上确界的定义可以验证这确实满足范数的定义##eq:norm ,并且对于有限维向量空间有$$||A||_p=\sum\limits_{x\neq 0}\left\| A \left(\frac{x}{||x||_p}\right)\right\|_p=\max\limits_{||x||_p=1}||Ax||_p.$$

向量范数与矩阵范数用于描述两个向量或矩阵的"接近"程度,在考察算法的精度时经常使用.

下面给出两个同维子空间的"距离"的定义,它用于描写这两个子空间的分离程度,在研究特征值算法的收敛速率时会用到.

设$S$为$V$的子空间,若矩阵$P$满足以下条件,则称$P$为$V$到$S$的正交投影映射:

 1. $\mathrm{span} (P)=S$.
 2. $P^2=P$.
 3. $P^T=P$.

由此定义可以得到,若$x\in V$,则$Px\in S$,且$(I-P)x\in S^\perp$.若$P_1$和$P_2$均为$V$到$S$的正交投影,则$$\| (P_1-P_2)x\|_2^2= (P_1z)^T (I-P_2)z+ (P_2z)^T (I-P_1)z=0,$$由此得到$P_1=P_2$,即向同一个子空间的正交投影映射是唯一的.若$S$有一组单位正交基$[\alpha_1,\cdots,\alpha_s]\equiv A$,则容易证明$P=AA^T$正是$V$到$S$的正交投影.

可以根据线性子空间与正交投影映射的一一对应关系定义两个同维子空间之间的距离.设$S_1$和$S_2$为$V$的线性子空间,$\mathrm{dim}S_1=\mathrm{dim}S_2$,$P_i$为$V$到$S_i$的正交投影 ($i=1,2$).则两线性子空间之间的距离定义为$$\mathrm{dist} (S_1,S_2)=\|P_1-P_2\|_2.$$我们用$\mathrm{dist} (S_1,S_1)$来表征两个同维子空间的分离程度.注意到$$0\leqslant \mathrm{dist} (S_1,S_2)\leqslant 1,$$且该距离为0当且仅当$S_1=S_2$,距离为1当且仅当$S_1\cap S_2^\perp \neq 0$.

** 特殊形式的矩阵

具有特殊结构的矩阵在数值线性代数中有着重要的地位.事实上,这是一个相当普遍的原则,即任何算法设计均应充分利用计算对象的特殊结构,而对于一般形式的矩阵则应进行约化,使之变成利于快速精确求解的矩阵形式.下面我们介绍一些今后会用到的特殊形式的矩阵.

上文中我们已经介绍了酉矩阵和Hermite矩阵.

置换矩阵$P$是指它的每一行和每一列都只有一个元素是1,其他元素都是0.从而$PX$与$XP$的结果是分别将$X$的行与列做一个置换.

反向置换阵$E$定义为
<latex>
\begin{equation*}
  E=
  \begin{bmatrix}
    & & 1\\
    & \adots & \\
    1& & 
  \end{bmatrix},
\end{equation*}
</latex>它左乘向量的效果是使向量各分量的顺序颠倒过来,即$$Ex=E (x_1,\cdots,x_n)^T= (x_n,\cdots,x_1).$$

稀疏矩阵是一大类非零元素的数目远小于矩阵元素总数的矩阵.可以看出置换矩阵就是一类特殊的稀疏矩阵.今后我们主要应用的稀疏矩阵是所谓带状阵,形象地说,就是其非零元素在矩阵中排成带状.若矩阵$A$满足$\forall i>j+s_1,A_{ij}=0$,则称$A$下带宽为$s_1$.类似地,若$\forall j>i+s_2,A_{ij}=0$,则称$A$上带宽为$s_2$.以下几类带状阵是值得注意的:

 1. 对角阵:若方阵的上下带宽均为0,则称为对角阵.直观地说,它的非零元素仅排列在对角线上.
 2. 双对角阵与Hessenberg阵:若方阵的上 (或下)带宽为1,则称为下 (或上)Hessenberg阵.若进一步要求上Hessenberg阵的上带宽为0,则成为下双对角阵,类似定义上双对角阵.
 3. 三对角阵:若方阵的上下带宽均为1,则称之为三对角阵.若它还是对称矩阵,则称为对称三对角阵.
 4. 块对角阵:若方阵经过合适的分块,可以看做对角阵,则称为块对角阵.

反向对称阵:若矩阵沿次对角线进行翻转后不变,则称之为反向对称阵.

我们还会用到Vandermonde阵与Toeplitz阵,将在应用时介绍.

* 数值型算法分析基础

由于数值计算通常是有误差的,因此对误差的分析构成算法分析的重要方面.按照误差的来源,可分为模型误差,观测误差与数值计算误差.前二者是数值分析中不考虑的.计算过程中产生的误差,主要起源与两个方面:

 1. 算法误差,即由于问题不可精确求解 (如高于5阶的一般矩阵的特征值问题)或算法的特殊设计 (如迭代法),算法本身就包含了误差的因素.典型的算法误差如截断误差,在运用级数对函数值进行逼近时,根据精度需要,只截断到某一项而舍去其余.(可参考"矩阵函数"一章中中"逼近算法"一节.)
 2. 舍入误差.这是由于计算机硬件存储数位的长度限制,在对原始数据与中间数据进行存储与计算时,可能产生误差.由于数值计算的特殊性,对这一误差来源往往要慎重考虑.为了消减舍入误差,在计算中有一些原则需要遵守,如避免相近数相减或大数与小数相加减等,以防精度严重丢失.

在一定条件下进行数值运算,对运算过程中产生误差的上界进行估计,从而评价算法的优劣,避免精度的过度损失是误差分析的目的.由于舍入误差的存在与不可避免,对于一个特定算法,我们往往要评价其结果对于中间数据的微小扰动的稳定性,这导致"条件数 (conditioning)"概念的出现.对一个数值问题,如果输入数据仅有微小扰动时,引起计算结果的相对误差很大,这就是病态 (ill-conditioned) 问题.例如,计算函数值$f (x)$时,若$x$有扰动$\Delta x-x^*$,其相对误差为$\frac{\Delta x}{x}$,函数值$f (x^*)$的相对误差为$\frac{f (x)-f (x^*)}{f (x)}$.相对误差比值$$\left|\frac{f (x)-f (x^*)}{f (x)}\right|\Big/\left|\frac{\Delta x}{x}\right|\approx\left|\frac{xf' (x)}{f (x)}\right|=C_p,$$ $C_p$称为函数计算问题的条件数.自变量相对误差一般不会很大,如果条件数$C_p$很大,将引起函数值相对误差很大,出现这种情况的问题就是病态问题.

对于矩阵计算的许多算法进行误差分析,发现问题的条件数往往与以下定义的矩阵的条件数有关 (多数是与变换矩阵的条件数成正比).因此,将其单独拿出来进行一些讨论是有意义的.

<definition label="conditioning" name="矩阵的条件数">
设$A\in \mathbb{C}^{n\times n}$为可逆矩阵,称乘积$\|A\|\|A^{-1}\|$为矩阵$A$(对于范数$\|\cdot\|$)的条件数,记为$\kappa (A)$.
</definition>

这里,"条件数"与一个矩阵相关,因而与该矩阵参与的所有数值计算均相关.如果$\kappa (A)$小,$A$称为良态的;如果$\kappa (A)$大,则$A$称为病态的.若$A$奇异,习惯写成$\kappa (A)=\infty$,以说明用奇异矩阵进行变换时,求解精度往往变得很差.

特别地,当以上定义中范数取矩阵的2-范数时,$A$的条件数与其SVD有着重要的联系.记非奇异矩阵$A$的SVD为$$A=U\mathrm{diag} (\sigma_1,\cdots,\sigma_n)V,$$则$$\|A\|_2=\sigma_1,\|A^{-1}\|_2=1/\sigma_n.$$因此2-范数意义如下$$\kappa (A)=\frac{\sigma_1}{\sigma_n}.$$这个比值可以解释为$\mathbb{C}^n$中单位球面在$A$作用下的象的哪个超椭球的离心率.注意到,这一离心率不小于1.当且仅当$A$为酉矩阵 (或正交矩阵) 时,$\kappa (A)=1$达到最小值.这正是数值算法设计中往往采用酉变换 (或正交变换)的原因.

根据以上关于矩阵条件数与SVD关系的讨论,我们还可以将矩阵的条件数推广到任意维数$n\times m$的矩阵上去.设$A\in \mathbb{C}^{n\times m}$的SVD为$$A=U\mathrm{diag} (\sigma_1,\cdots,\sigma_r,0,\cdots,0)V,$$定义$A$的 (广义)条件数为$$\kappa_2A=\sigma_1/\sigma_r.$$这样定义的条件数也能反映$A$用作变换矩阵时数值运算的精度.

* 矩阵乘法

在有关矩阵的计算中,矩阵乘法具有基础性的意义.对于$n$阶矩阵的乘法,常规算法具有$O(n^3)$的复杂度,因此加速矩阵计算是很重要的.

** 基于向量内积算法的Winograd加速算法

以下讨论主要来自文献<cite>Win68</cite>.
<algorithm label="alg:WinogradInnerProduct" name="Winograd内积算法">
设$x=(x_1,\cdots,x_n)^T$,$y=(y_1,\cdots,y_n)^T$,记$\xi=\sum\limits_{j=1}^{\lfloor n/2\rfloor}x_{2j-1}x_{2j}$,$\eta=\sum\limits^{\lfloor n/2\rfloor}_{j=1}y_{2j-1}y_{2j}$,则内积$(x,y)$可由下式给出:

<latex>
\begin{equation*}
  (x,y)=
  \begin{cases}
    \sum\limits^{\lfloor n/2\rfloor}_{j=1}(x_{2j-1}+y_{2j})(x_{2j}+y_{2j-1})-\xi-\eta,& \text{$n$是偶,} \\
    \sum\limits^{\lfloor n/2\rfloor}_{j=1}(x_{2j-1}+y_{2j})(x_{2j}+y_{2j-1})-\xi-\eta+x_ny_n,& \text{$n$是奇.}
  \end{cases}
\end{equation*}
</latex>
</algorithm>

将这种算法用于$C=AB$的矩阵元素运算时,由于减少重复计算$\xi$,$\eta$,可使计算所需的乘法次数减半,但同时使所需的加法运算增加.Winograd算法也是$O(N^3)$的算法,仅适用于小规模的矩阵求积运算,且由于该算法破坏了向量内积的整体间运算,同时增加了内存开销,因而其算法改进价值并不很大.

** Strassen算法

Strassen算法(1968)是一种分治策略的算法.它以分块矩阵运算为基础.

下面介绍改进型Strassen算法,它较原始算法<cite>Str68</cite>需要更少的矩阵加法运算<cite>mca</cite>.

<algorithm label="alg:Str" name="Strassen算法">
设$A$,$B$为$n$阶矩阵,必要时通过补充零行零列加以扩充,将其分块:

<latex>
\begin{gather*}
  A=
  \begin{bmatrix}
    A_{11}& A_{12}\\
    A_{21}& A_{22}
  \end{bmatrix}
  ,B=
  \begin{bmatrix}
    B_{11}& B_{12}\\
    B_{21}& B_{22}
  \end{bmatrix}
 ,C=AB=
  \begin{bmatrix}
    C_{11}& C_{12}\\
    C_{21}& C_{22}
  \end{bmatrix}.
\end{gather*}
</latex>
进行如下递归运算:
 1. 若$n\leq l$($l$为递归下界),做直接乘法.
 2. 计算
<latex>
\begin{gather*}
  S_1=A_{21}+A_{22},S_2=S_1-A_{11},S_3=A_{11}-A_{21},S_4=A_{12}-S_2,\\
  T_1=B_{12}-B_{11},T_2=B_{22}-T_1,T_3=B_{22}-B_{11},T_4=T_2-B_{21}.
\end{gather*}
</latex>
 3. 计算
<latex>
\begin{gather*}
  P_1=A_{11}B_{11},P_2=A_{12}B_{21},P_3=S_4B_{22},P_4=B_{22}T_4,\\
  P_5=S_1T_1,P_6=S_2T_2,P_7=S_3T_3.
\end{gather*}
</latex>
 4. 计算
<latex>
\begin{gather*}
  U_1=P_1+P_2,U_2=P_1+P_6,U_3=U_2+P_7,U_4=U_2+P_5,\\
  U_5=U_4+P_3,U_6=U_3-P_4,U_7=U_3+P_5.
\end{gather*}
</latex>
 5. 返回
<latex>
\begin{equation*}
\begin{bmatrix}
   U_1& U_5\\
   U_6& U_7
\end{bmatrix}.
\end{equation*}
</latex>
</algorithm>

以上算法的正确性直接代入即可验证.可以看出,每次递归需要7次乘法与15次加法,从而其算法复杂度是$O(N^{\log_27})\simeq O(N^{2.808})$.

Strassen算法在之后有许多推广,最优渐进复杂度可以降到$O(N^{2.376})$.但在实际中,仅当$N$极大时才有价值,故通常并不采用.可参考<cite>Pan84</cite>.