東京大学 情報理工学系研究科 コンピュータ科学専攻 2022年2月実施 問題1

标签：

• Tokyo-University
• Lower-Upper-Decomposition
• Systems-of-Linear-Equations
• Algorithm-Complexity

Author

zephyr

Description

For a matrix denoted by $\mathbf{M}$, we write $\mathbf{M}_{i:j,k:l}$ for the submatrix formed from the rows $i$ through $j$ and the columns from $k$ through $l$ of the matrix. We also denote $\mathbf{M}_{i:i,k:l}$ and $\mathbf{M}_{i:j,k:k}$ by $\mathbf{M}_{i,k:l}$ and $\mathbf{M}_{i:j,k}$, respectively. (Note that they are row and column vectors, respectively.)

Assume below that $\mathbf{A}$ is an $n \times n$ nonsingular matrix which can be LU factorized. In addition, assume that a positive integer $w$ exists such that all entries in the $k$-th sub- and super-diagonal of $\mathbf{A}$ are zero if $k > w$. Namely, if $|i - j| > w$, $A_{i,j}$ is zero.

We denote the LU factorization of $\mathbf{A}$ by $\mathbf{A} = \mathbf{LU}$, where $\mathbf{L}$ is a lower triangular matrix with unit diagonal elements, and $\mathbf{U}$ is an upper triangular matrix.

The following algorithm P computes $\mathbf{L}$ and $\mathbf{U}$ from the input $\mathbf{A}$ in-place in such a way that the strictly lower triangular part and the upper triangular part of $\mathbf{A}$ will be $\mathbf{L}$ and $\mathbf{U}$, respectively. In other words, $A_{i,j}$ will be $L_{i,j}$ for $i > j$, and $A_{i,j}$ will be $U_{i,j}$ for $i \leq j$, when it terminates.

Algorithm P

for (k = 1; k < n; k = k + 1) {
    A[k+1:n,k] <- A[k+1:n,k] * (1 / A[k,k]);
    A[k+1:n,k+1:n] <- A[k+1:n,k+1:n] - A[k+1:n,k] * A[k,k+1:n];
}

Answer the following questions.

(1) Compute the LU factorization of the following matrix.

$$\begin{pmatrix} 1 & 1 & 0 & 0 \\ -1 & 1 & 1 & 0 \\ 0 & -2 & 2 & 1 \\ 0 & 0 & -3 & 3 \end{pmatrix}$$

(2) Prove that $L_{i,j} = 0$ and $U_{i,j} = 0$ if $|i - j| > w$.

(3) Assume that $n \gg 1$, $w \ll n$, and $m \gg n$. We wish to solve a set of linear systems,

$$\mathbf{A}\mathbf{x}_i = \mathbf{b}_i, \quad i = 1, 2, \ldots, m,$$

where $\mathbf{b}_i \, (i = 1, 2, \ldots, m)$ are the constant vectors, and $\mathbf{x}_i \, (i = 1, 2, \ldots, m)$ are the unknown vectors. Explain the relative merits of the following methods (I) and (II) with respect to the amount of arithmetic operations.

• (I) Compute $\mathbf{A}^{-1}$ first, and then compute each $\mathbf{x}_i$ as $\mathbf{A}^{-1} \mathbf{b}_i$.
• (II) Compute the LU factorization of $\mathbf{A}$ first, and then solve $\mathbf{L} \mathbf{y}_i = \mathbf{b}_i$ for $\mathbf{y}_i$. Solve $\mathbf{U} \mathbf{x}_i = \mathbf{y}_i$ for $\mathbf{x}_i$ lastly.

---

对于矩阵 $\mathbf{M}$，我们写 $\mathbf{M}_{i:j,k:l}$ 表示从矩阵的第 $i$ 行到第 $j$ 行以及第 $k$ 列到第 $l$ 列构成的子矩阵。我们还将 $\mathbf{M}_{i:i,k:l}$ 和 $\mathbf{M}_{i:j,k:k}$ 分别表示为 $\mathbf{M}_{i,k:l}$ 和 $\mathbf{M}_{i:j,k}$。（注意它们分别是行向量和列向量。）

假设 $\mathbf{A}$ 是一个 $n \times n$ 的非奇异矩阵，可以进行 LU 分解。另外，假设存在一个正整数 $w$，使得当 $|i - j| > w$ 时，$\mathbf{A}$ 在第 $k$ 条下对角线和超对角线上的所有元素都为零。即，如果 $|i - j| > w$，则 $A_{i,j}$ 为零。

我们用 $\mathbf{A} = \mathbf{LU}$ 来表示 $\mathbf{A}$ 的 LU 分解，其中 $\mathbf{L}$ 是具有单位对角元素的下三角矩阵，$\mathbf{U}$ 是上三角矩阵。

下面的算法 P 从输入 $\mathbf{A}$ 中就地计算 $\mathbf{L}$ 和 $\mathbf{U}$，以这样的方式，$\mathbf{A}$ 的严格下三角部分和上三角部分将分别是 $\mathbf{L}$ 和 $\mathbf{U}$。换句话说，当它终止时，$A_{i,j}$ 将为 $L_{i,j}$（当 $i > j$）和 $U_{i,j}$（当 $i \leq j$）。

算法 P

for (k = 1; k < n; k = k + 1) {
    A[k+1:n,k] <- A[k+1:n,k] * (1 / A[k,k]);
    A[k+1:n,k+1:n] <- A[k+1:n,k+1:n] - A[k+1:n,k] * A[k,k+1:n];
}

回答以下问题。

(1) 计算以下矩阵的 LU 分解。

$$\begin{pmatrix} 1 & 1 & 0 & 0 \\ -1 & 1 & 1 & 0 \\ 0 & -2 & 2 & 1 \\ 0 & 0 & -3 & 3 \end{pmatrix}$$

(2) 证明 $L_{i,j} = 0$ 和 $U_{i,j} = 0$ 如果 $|i - j| > w$。

(3) 假设 $n \gg 1$，$w \ll n$，并且 $m \gg n$。我们希望求解一组线性系统，

$$\mathbf{A} \mathbf{x}_i = \mathbf{b}_i, \quad i = 1, 2, \ldots, m,$$

其中 $\mathbf{b}_i \, (i = 1, 2, \ldots, m)$ 是常量向量，$\mathbf{x}_i \, (i = 1, 2, \ldots, m)$ 是未知向量。解释以下方法 (I) 和 (II) 在算术运算量方面的相对优点。

• (I) 首先计算 $\mathbf{A}^{-1}$，然后将每个 $\mathbf{x}_i$ 计算为 $\mathbf{A}^{-1} \mathbf{b}_i$。
• (II) 首先计算 $\mathbf{A}$ 的 LU 分解，然后解 $\mathbf{L} \mathbf{y}_i = \mathbf{b}_i$ 得到 $\mathbf{y}_i$。最后解 $\mathbf{U} \mathbf{x}_i = \mathbf{y}_i$。

题目描述

对矩阵 $M$，以 $M_{i:j,k:l}$ 表示由第 $i$ 至 $j$ 行、第 $k$ 至 $l$ 列组成的子矩阵；$M_{i,k:l}$ 和 $M_{i:j,k}$ 分别表示相应行向量、列向量。

设 $A$ 是可作 LU 分解的 $n\times n$ 非奇异带状矩阵：存在正整数 $w$，使 $|i-j|>w$ 时 $A_{i,j}=0$。记 $A=LU$，其中 $L$ 为单位下三角矩阵，$U$ 为上三角矩阵。题中原地算法 P 对 $k=1,\ldots,n-1$ 依次执行

$$A_{k+1:n,k}\leftarrow A_{k+1:n,k}/A_{k,k},$$

$$A_{k+1:n,k+1:n}\leftarrow A_{k+1:n,k+1:n} -A_{k+1:n,k}A_{k,k+1:n},$$

终止后 $A$ 的严格下三角部分存放 $L$ 的相应元素，上三角部分存放 $U$。

（1）计算矩阵

$$\begin{pmatrix} 1&1&0&0\\ -1&1&1&0\\ 0&-2&2&1\\ 0&0&-3&3 \end{pmatrix}$$

的 LU 分解。

（2）证明当 $|i-j|>w$ 时，$L_{i,j}=0$ 且 $U_{i,j}=0$。

（3）设 $n\gg1$、$w\ll n$、$m\gg n$，需要求解 $Ax_i=b_i\ (i=1,\ldots,m)$。从算术运算量比较：

• 先求 $A^{-1}$，再计算每个 $x_i=A^{-1}b_i$；
• 先对 $A$ 作 LU 分解，再分别前代求 $Ly_i=b_i$、回代求 $Ux_i=y_i$。

说明两种方法的相对优劣。

考点

• LU 分解与原地消元：按高斯消元过程计算单位下三角和上三角因子。
• 带状矩阵的零结构：证明无主元交换的消元不会在带宽之外产生填充。
• 多右端线性方程组：比较显式求逆与一次分解、多次三角求解的运算量，并利用带宽降低复杂度。

Kai

(1)

We start by applying the given algorithm P to compute the LU factorization of $\mathbf{A}$.

Step-by-step computation:

1. Initial matrix $\mathbf{A}$:

$$\mathbf{A} = \begin{pmatrix} 1 & 1 & 0 & 0 \\ -1 & 1 & 1 & 0 \\ 0 & -2 & 2 & 1 \\ 0 & 0 & -3 & 3 \end{pmatrix}$$

1. First iteration ($k = 1$):
   • Update $A[2:4,1]$ by dividing by $A[1,1]$:

$$A[2:4,1] \leftarrow A[2:4,1] \times (1 / A[1,1]) = \begin{pmatrix} -1 \\ 0 \\ 0 \end{pmatrix} \times 1 = \begin{pmatrix} -1 \\ 0 \\ 0 \end{pmatrix}$$

• Update $A[2:4,2:4]$:

$$A[2:4,2:4] \leftarrow A[2:4,2:4] - A[2:4,1] \times A[1,2:4] = \begin{pmatrix} 1 & 1 & 0 \\ -2 & 2 & 1 \\ 0 & -3 & 3 \end{pmatrix} - \begin{pmatrix} -1 \\ 0 \\ 0 \end{pmatrix} \times \begin{pmatrix} 1 & 0 & 0 \end{pmatrix}$$

Simplifying:

$$A[2:4,2:4] = \begin{pmatrix} 2 & 1 & 0 \\ -2 & 2 & 1 \\ 0 & -3 & 3 \end{pmatrix}$$

1. Second iteration ($k = 2$):
   • Update $A[3:4,2]$ by dividing by $A[2,2]$:

$$A[3:4,2] \leftarrow A[3:4,2] \times (1 / A[2,2]) = \begin{pmatrix} -2 \\ 0 \end{pmatrix} \times \frac{1}{2} = \begin{pmatrix} -1 \\ 0 \end{pmatrix}$$

• Update $A[3:4,3:4]$:

$$A[3:4,3:4] \leftarrow A[3:4,3:4] - A[3:4,2] \times A[2,3:4] = \begin{pmatrix} 2 & 1 \\ -3 & 3 \end{pmatrix} - \begin{pmatrix} -1 \\ 0 \end{pmatrix} \times \begin{pmatrix} 1 & 0 \end{pmatrix}$$

Simplifying:

$$A[3:4,3:4] = \begin{pmatrix} 3 & 1 \\ -3 & 3 \end{pmatrix}$$

1. Third iteration ($k = 3$):
   • Update $A[4,3]$ by dividing by $A[3,3]$:

$$A[4,3] \leftarrow A[4,3] \times (1 / A[3,3]) = -3 \times \frac{1}{3} = -1$$

• Update $A[4,4]$:

$$A[4,4] \leftarrow A[4,4] - A[4,3] \times A[3,4] = 3 - (-1) \times 1 = 4$$

The resulting matrices $\mathbf{L}$ and $\mathbf{U}$ are:

$$\mathbf{L} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 0 & -1 & 1 & 0 \\ 0 & 0 & -1 & 1 \end{pmatrix}, \quad \mathbf{U} = \begin{pmatrix} 1 & 1 & 0 & 0 \\ 0 & 2 & 1 & 0 \\ 0 & 0 & 3 & 1 \\ 0 & 0 & 0 & 4 \end{pmatrix}$$

(2)

Given that $\mathbf{A}$ is banded with bandwidth $w$, meaning $A_{i,j} = 0$ for $|i - j| > w$. The LU factorization $\mathbf{A} = \mathbf{LU}$ maintains this band structure:

• Lower Triangular Matrix $\mathbf{L}$:

  • By definition, $\mathbf{L}$ has non-zero entries only in the lower triangular part and on the diagonal. For $|i - j| > w$, $L_{i,j} = 0$ because $A_{i,j} = 0$ for these entries. Hence, $\mathbf{L}$ also maintains the band structure.

• Upper Triangular Matrix $\mathbf{U}$:

  • Similarly, $\mathbf{U}$ has non-zero entries only in the upper triangular part. For $|i - j| > w$, $U_{i,j} = 0$ because $A_{i,j} = 0$ for these entries. Hence, $\mathbf{U}$ maintains the band structure.

Therefore, $L_{i,j} = 0$ and $U_{i,j} = 0$ if $|i - j| > w$.

(3)

Method (I): Compute $\mathbf{A}^{-1}$ first, then $\mathbf{x}_i = \mathbf{A}^{-1} \mathbf{b}_i$

• Steps:
  1. Compute $\mathbf{A}^{-1}$, which requires $O(n^3)$ operations.
  2. Compute each $\mathbf{x}_i = \mathbf{A}^{-1} \mathbf{b}_i$, which requires $O(n^2)$ operations per $\mathbf{x}_i$.
  3. Total operations for $m$ systems: $O(n^3) + m \cdot O(n^2) = O(n^3 + mn^2)$.

Method (II): Compute the LU factorization, then solve $\mathbf{L} \mathbf{y}_i = \mathbf{b}_i$ and $\mathbf{U} \mathbf{x}_i = \mathbf{y}_i$

• Steps:

  1. Compute the LU factorization of $\mathbf{A}$, which requires $O(n^3)$ operations.

  2. Solve $\mathbf{L} \mathbf{y}_i = \mathbf{b}_i$ and $\mathbf{U} \mathbf{x}_i = \mathbf{y}_i$ for each $\mathbf{x}_i$, which requires $O(n^2)$ operations per $\mathbf{x}_i$.

  3. Total operations for $m$ systems: $O(n^3) + m \cdot O(n^2) = O(n^3 + mn^2)$.

Comparison

• Arithmetic Operations:

  • Both methods have the same asymptotic complexity $O(n^3 + mn^2)$.

• Efficiency in Practice:

  • Method (II) is generally more efficient and numerically stable because it avoids directly computing the inverse of $\mathbf{A}$, which can be computationally expensive and prone to numerical errors.

• Storage Requirements:

  • Method (I) requires storing $\mathbf{A}^{-1}$, which can be memory-intensive.
  • Method (II) only requires storing the LU factors, which is typically more storage-efficient.

Conclusion: Method (II) is preferred due to its numerical stability and potentially lower storage requirements, even though both methods have the same theoretical computational complexity.

Knowledge

LU分解 稀疏矩阵 复杂度分析 递归 内存效率

难点解题思路

对大型稀疏矩阵，直接计算逆矩阵通常会引入数值不稳定性，且计算开销较大。LU 分解通过三角矩阵简化了问题，解决了计算复杂度和数值稳定性的问题。

解题技巧和信息

对于稀疏矩阵的 LU 分解，注意带状结构的保持，有助于减少计算量。在求解线性方程组时，优先考虑分解矩阵的方法以优化计算效率。

重点词汇

• LU factorization LU 分解
• band matrix 带状矩阵
• numerical stability 数值稳定性
• computational complexity 计算复杂度

参考资料

1. Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations. Johns Hopkins University Press. Chap. 3.
2. Trefethen, L. N., & Bau, D. (1997). Numerical Linear Algebra. SIAM. Chap. 21.