三角形方程组和三角分解
前代法
求解下三角形方程组
\[Ly = b \]其中 \(b=(b_1,\cdots,b_n)^T\in\mathbb{R}^n\) 已知， \(y=(y_1,\cdots,y_n)^T\in\mathbb{R}^n\) 未知，而
\[L = \left( \begin{matrix} l_{11}\\ l_{21} & l_{22}\\ \vdots & \vdots & \ddots\\ l_{n1} & l_{n2} & \cdots & l_{nn} \end{matrix} \right) \]是已知的非奇异下三角阵，且 \(l_{ii}\neq 0,\ i=1,\cdots,n\) .
前代法通过公式
\[y_i = \left(b_i-\sum_{j=1}^{i-1}l_{ij}y_j\right)\bigg{/}l_{ii} \]从上到下求解方程组.
代码实现
vector<double> lower(vector<vector<double>> &A, vector<double> &b)
{
int n = b.size();
vector<double> x(n);
double sum;
x[0] = b[0];
for (int i = 1; i < n; i++)
{
sum = 0;
for (int j = 0; j < i; j++)
{
sum += A[i][j] * x[j];
}
x[i] = (b[i] - sum) / A[i][i];
}
return x;
}
回代法
求解上三角形方程组
\[Ux = y \]其中 \(y=(y_1,\cdots,y_n)^T\in\mathbb{R}^n\) 已知， \(x=(x_1,\cdots,x_n)^T\in\mathbb{R}^n\) 未知，而
\[U = \left( \begin{matrix} u_{11} & u_{12} & \cdots & u_{1n}\\ & u_{22} & \cdots & u_{2n}\\ & & \ddots & \vdots\\ & & & u_{nn}\\ \end{matrix} \right) \]是已知的非奇异上三角阵，且 \(u_{ii}\neq 0,\ i=1,\cdots,n\) .
回代法通过公式
从下到上求解方程组.
代码实现
vector<double> upper(vector<vector<double>> &A, vector<double> &b)
{
int n = b.size();
vector<double> x(n);
double sum;
for (int i = n - 1; i >= 0; i--)
{
sum = 0;
for (int j = i + 1; j < n; j++)
{
sum += A[i][j] * x[j];
}
x[i] = (b[i] - sum) / A[i][i];
}
return x;
}
Gauss 变换
通过一系列初等变换，逐步用下三角阵对 \(A\) 进行约化
\[L_k = I - l_ke_k^T,\quad l_k = (0,\cdots,0,l_{k+1,k},\cdots, l_{nk})^T \]左乘 \(L_k\) 进行行变换
\[L_k = \left( \begin{matrix} 1\\ & \ddots\\ & & 1\\ & & -l_{k+1,k} & 1\\ & & \vdots & & \ddots\\ & & -l_{nk}& & & 1 \end{matrix} \right) \]这种类型的初等下三角阵称为 Gauss 变换，向量 \(l_k\) 称为 Gauss 向量.
高斯向量具有重要性质，因为 \(e_k^Tl_k = 0\) ，从而
\[(I-l_ke_k^T)(I+l_ke_k^T) = I - l_ke_k^Tl_ke_k^T = I \]即有
\[L_k^{-1} = I + l_ke_k^T \]高斯变换的逆容易得到.
三角分解的计算
设 \(A\in\mathbb{R}^{n\times n}\) ，则 \(A\) 的三角分解指分解 \(A=LU\) ，其中 \(L\in\mathbb{R}^{n\times n}\) 为下三角阵， \(U\in\mathbb{R}^{n\times n}\) 为上三角阵，此分解也称为 LU 分解.
我们通过左乘一系列高斯变换来消去 \(A\) 的下三角部分
\[A^{(k-1)} = L_{k-1}\cdots L_1A = \left( \begin{matrix} A_{11}^{(k-1)} & A_{12}^{(k-1)}\\ 0 & A_{22}^{(k-1)} \end{matrix} \right) \]其中 \(A_{11}^{(k-1)}\) 是 \(k-1\) 阶上三角阵，而
\[A_{22}^{(k-1)} = \left( \begin{matrix} a_{kk}^{(k-1)} & \cdots & a_{kn}^{(k-1)}\\ \vdots & \ddots & \vdots\\ a_{nk}^{(k-1)} & \cdots & a_{nn}^{(k-1)} \end{matrix} \right) \]如果 \(a_{kk}^{(k-1)}\neq 0\) ，则可以再确定一个高斯变换 \(L_k = I-l_ke_k^T,\ l_{ki} = a_{ki}^{(k-1)}/a_{kk}^{(k-1)},\ i=k+1,\cdots,n\) .
最终我们得到
\[L= (L_{n-1}L_{n-2}\cdots L_1)^{-1},\quad U = A^{(n-1)},\quad A = LU \]利用高斯变换的特点
\[L = (L_{n-1}L_{n-2}\cdots L_1)^{-1} = L_1^{-1}\cdots L_{n-1}^{-1} = (I+l_1e_1^T)\cdots (I+l_{n-1}e_{n-1}^T) = I + l_1e_1^T + \cdots l_{n-1}e_{n-1}^T \]即 \(L\) 有如下形状
\[L = \left( \begin{matrix} 1\\ l_{21} & 1\\ \vdots & \vdots & \ddots\\ l_{n1} & l_{n2} & \cdots & 1 \end{matrix} \right) \]此方法称为 Gauss 消去法；通常称 \(a_{kk}^{(k-1)}\) 为主元.
高斯消去法的存储
由于每次消元都与之前的列无关，并且下三角阵对角元为 \(1\) ，可以省略，因此可以将每一步高斯变换的结果直接存放在原先的列中，我们只展示最后的结果
\[A^{(n-1)} = \left( \begin{matrix} u_{11} & u_{12} & \cdots & u_{1n}\\ l_{21} & u_{22} & \cdots & u_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ l_{n1} & l_{n2} & \cdots & u_{nn} \end{matrix} \right) \]高斯消去法是原地消去，不需要额外的存储空间.
代码实现
void LU(vector<vector<double>> &A)
{
int n = A.size();
for (int i = 0; i < n; i++)
{
double key = A[i][i];
// 无法进行高斯变换
if (key == 0)
{
return;
}
// 第 i 列高斯变换
for (int j = i + 1; j < n; j++)
{
A[j][i] = A[j][i] / key;
for (int k = i + 1; k < n; k++)
{
A[j][k] -= A[i][k] * A[j][i];
}
}
}
}
定理 1.1.1 主元 \(a_{ii}^{(i-1)},\ i=1,\cdots,k\) 均不为零的充要条件是 \(A\) 的 \(i\) 阶顺序主子式 \(A_i\) 非奇异.
证明 容易应用归纳法证明.
定理 1.1.2 若 \(A\in\mathbb{R}^{n\times n}\) 的顺序主子式 \(A_k\in\mathbb{R}^{k\times k},\ k=1,\cdots,n-1\) 非奇异，则存在唯一分解 \(A=LU\) .
选主元三角分解
由于计算机精度的限制，如果主元太小，则会导致高斯变换的误差增大，因此我们需要选取尽可能小的元素作为主元.
\[A^{(k-1)} = \left( \begin{matrix} A_{11}^{(k-1)} & A_{12}^{(k-1)}\\ 0 & A_{22}^{(k-1)} \end{matrix} \right) \]我们在每一步高斯变换前，进行初等变换 \(I_{pq}\) ，左乘 \(I_{pq}\) 会交换 \(p,q\) 行，右乘 \(I_{pq}\) 会交换 \(p,q\) 列。选取
\[\left|a_{pq}^{(k-1)}\right| = \max\left\{\left|a_{ij}^{(k-1)}\right|:k\le i,j\le n\right\} \]如果 \(a_{pq}^{(k-1)} = 0\) ，说明 \(A\) 秩为 \(k-1\) ，结束过程；否则利用 \(I_{kp}\) 与 \(I_{kq}\) 交换 \(a_{kk}\) 与 \(a_{pq}\) ，最终得到
\[L_rP_r\cdots L_1P_1AQ_1\cdots Q_r = U \]其中 \(P_k = I_{kp},\ Q_k={I_{kq}}\) ，此过程称为全主元 Gauss 消去法.
令
\[\begin{aligned} Q &= Q_1\cdots Q_r\\ P &= P_r\cdots P_1\\ L &= P(L_rP_r\cdots L_1P_1)^{-1} \end{aligned} \]则有
\[PAQ = LU \]称为全主元三角分解；注意到 \(P_k^{-1}=P_k\) ，有
\[L = P_r\cdots P_1\cdot (L_rP_r\cdots L_1P_1)^{-1} = P_r\cdots P_2L_1^{-1}P_2L_2^{-1}\cdots P_rL_r^{-1} \]如果我们以 \(L_1^{-1}\) 为中心，令
\[L^{(1)} = L_1^{-1},\quad L^{(k)} = P_kL^{(k-1)}P_kL_k^{-1} \]其中的每一步的 \(P_kL^{(k-1)}P_kL_k^{-1}\) 只改变 \(I_{n-k+1}\) 而不会改变下三角结构
\[L^{(k-1)} = \left( \begin{matrix} L_{11}^{(k-1)} & 0\\ L_{21}^{(k-1)} & I_{n-k+1} \end{matrix} \right) \]因此结果 \(L\) 仍为单位下三角阵；并且由于每一步选取最大的主元，由高斯变换的构造得到 \(L\) 中元素的模不会超过 \(1\) .
定理 1.2.1 若 \(A\in\mathbb{R}^{n\times n}\) ，则有排列方阵 \(P,Q\in\mathbb{R}^{n\times n}\) 使得
\[PAQ = LU \]其中 \(L\) 中所有元素的模不大于 \(1\) ， \(U\) 中非零对角元的个数恰为 \(A\) 的秩.
虽然全主元高斯消去法弥补了不选主元的高斯消去法的不足，但仍需要付出极大的代价，因为在 \(A\) 非奇异的情况下，选主元需进行
\[\sum_{k=1}^{n-1}(n-k+1)^2 = \dfrac{1}{3}n^3+O(n^2) \]次比较判断。于是我们考虑列主元高斯消去法
\[\left|a_{pk}^{(k-1)}\right| = \max\left\{\left|a_{ik}^{(k-1)}\right|:k\le i\le n\right\} \]这样只进行每一列的选主元，能够大量减少时间.
代码实现
// 我们同时进行对 A, b 的变换
void column_lu(vector<vector<double>> &A, vector<double> &b)
{
int n = A.size();
double max, tmp; // 记录最大值
int index; // 记录要交换的主元位置
for (int i = 0; i < n; i++)
{
max = 0;
index = i;
// 选取列主元
for (int p = i; p < n; p++)
{
tmp = fabs(A[p][i]);
if (tmp > max)
{
max = tmp;
index = p;
}
}
// 在这里奇异，说明 A 的秩为 i
if (max == 0)
{
cout << "Singular Matrix!" << endl;
return;
}
// 交换主元
for (int q = 0; q < n; q++)
{
tmp = A[i][q];
A[i][q] = A[index][q];
A[index][q] = tmp;
}
// 同时交换 b
tmp = b[i];
b[i] = b[index];
b[index] = tmp;
// 正常的高斯消去法
for (int j = i + 1; j < n; j++)
{
A[j][i] = A[j][i] / A[i][i];
for (int k = i + 1; k < n; k++)
{
A[j][k] = A[j][k] - A[i][k] * A[j][i];
}
}
}
}
平方根法
平方根法又叫做 Cholesky 分解法，是求解对称正定线性方程组最常用的方法之一.
定理 1.3.1 (Cholesky 分解定理) 若 \(A\in\mathbb{R}^{n\times n}\) 对称正定，则存在对角元均为正数的下三角阵 \(L\in\mathbb{R}^{n\times n}\) ，使得
\[A = LL^T \]称为 Cholesky 分解，其中 \(L\) 称为 \(A\) 的 Cholesky 因子.
证明 由于对称正定，因此其全部主子阵均正定，则存在 \(\widetilde{L},U\) 使得 \(A = \widetilde{L}U\) ，令
\[D = \mathrm{diag}(u_{11},\cdots,u_{nn}),\quad \widetilde{U} = D^{-1}U \]则有
\[\widetilde{U}^TD\widetilde{L}^T = A^T = A = \widetilde{L}D\widetilde{U}\ \Rightarrow\ \widetilde{L}^T\widetilde{U}^{-1} = D^{-1}\widetilde{U}^{-T}\widetilde{L}D \]注意到左边是单位上三角阵，右边是下三角阵，因此两边均为单位矩阵，于是 \(\widetilde{U} = \widetilde{L}^T,\ A=\widetilde{L}D\widetilde{L}^T\) ；由于 \(D\) 对角元为正数，则取
\[L = \widetilde{L}\mathrm{diag}(\sqrt{u_{11}},\cdots,\sqrt{u_{nn}}) \]有 \(A = LL^T\) ，此即为满足要求的分解.
平方根法通过比较元素来计算分解，对比 \(A\) 和 \(LL^T\) 两边的对应元素得到
\[l_{kk} = \left(a_{ik}-\sum_{p=1}^{k-1}l_{kp}^2\right)^{\frac{1}{2}},\quad l_{ik} = \left(a_{ik}-\sum_{p=1}^{k-1}l_{ip}l_{kp}\right)\bigg{/}l_{kk} \]从左向右，每次计算一列元素，先算出对角元素，然后推导其余元素.
代码实现
void cholesky(vector<vector<double>> &A)
{
double sum = 0;
int n = A.size();
for (int i = 0; i < n; i++)
{
sum = 0;
for (int p = 0; p < i; p++)
{
sum += A[i][p] * A[i][p];
}
A[i][i] = sqrt(A[i][i] - sum);
// 奇异矩阵
if (A[i][i] == 0)
{
return;
}
for (int j = i + 1; j < n; j++)
{
sum = 0;
for (int k = 0; k < i; k++)
{
sum += A[j][k] * A[i][k];
}
A[j][i] = (A[j][i] - sum) / A[i][i];
A[i][j] = A[j][i];
}
}
}
为了避免开方运算，我们使用改进的平方根法
\[A = LDL^T \]其中 \(L\) 是单位下三角阵， \(D\) 是对角元均为正的对角阵.
仍然通过对比元素计算分解
\[d_j = a_{jj} - \sum_{k=1}^{j-1}l_{jk}d_kl_{jk},\quad l_{ij} = \left(a_{ij}-\sum_{k=1}^{j-1}l_{ik}d_{k}l_{jk}\right)\bigg{/}d_{j} \]在实际计算时，将 \(L\) 的严格下三角元素存储在 \(A\) 的对应位置， \(D\) 的对角元存储在 \(A\) 的对应对角位置.
代码实现
void cholesky_d(vector<vector<double>> &A)
{
double sum = 0;
int n = A.size();
for (int i = 0; i < n; i++)
{
sum = 0;
for (int p = 0; p < i; p++)
{
sum += A[i][p] * A[i][p] * A[p][p];
}
A[i][i] -= sum;
// 奇异矩阵
if (A[i][i] == 0)
{
return;
}
for (int j = i + 1; j < n; j++)
{
sum = 0;
for (int k = 0; k < i; k++)
{
sum += A[j][k] * A[i][k] * A[k][k];
}
A[i][j] = A[j][i] - sum;
A[j][i] = A[i][j] / A[i][i];
}
}
}
追赶法
三对角阵的分解
\[\left( \begin{matrix} a_{1} & b_2\\ c_{2} & a_{2} & b_3\\ & \ddots & \ddots & \ddots\\ & & c_{n-1} & a_{n-1} & b_{n}\\ & & & c_{n} & a_n \end{matrix} \right) = \left( \begin{matrix} \alpha_{1}\\ \gamma_{2} & \alpha_{2}\\ & \ddots & \ddots\\ & & \gamma_{n-1} & \alpha_{n-1} \\ & & & \gamma_{n} & \alpha_n \end{matrix} \right) \left( \begin{matrix} 1 & \beta_2\\ & 1 & \beta_3\\ & & \ddots & \ddots\\ & & & 1 & \beta_{n}\\ & & & & 1 \end{matrix} \right) \]对比元素得到 \(\alpha_1 = a_1\) ，而
\[\gamma_i = c_i,\quad \beta_i = b_i/\alpha_{i-1},\quad \alpha_i = a_i-\gamma_i\beta_i,\quad i=2,\cdots,n \]此方法称为追赶法.