跳到主要内容

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

Author

zephyr, 祭音Myyura

Description

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

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

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

The following algorithm P computes L\mathbf{L} and U\mathbf{U} from the input A\mathbf{A} in-place in such a way that the strictly lower triangular part and the upper triangular part of A\mathbf{A} will be L\mathbf{L} and U\mathbf{U}, respectively. In other words, Ai,jA_{i,j} will be Li,jL_{i,j} for i>ji > j, and Ai,jA_{i,j} will be Ui,jU_{i,j} for iji \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.

(1100111002210033)\begin{pmatrix} 1 & 1 & 0 & 0 \\ -1 & 1 & 1 & 0 \\ 0 & -2 & 2 & 1 \\ 0 & 0 & -3 & 3 \end{pmatrix}

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

(3) Assume that n1n \gg 1, wnw \ll n, and mnm \gg n. We wish to solve a set of linear systems,

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

where bi(i=1,2,,m)\mathbf{b}_i \, (i = 1, 2, \ldots, m) are the constant vectors, and xi(i=1,2,,m)\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 A1\mathbf{A}^{-1} first, and then compute each xi\mathbf{x}_i as A1bi\mathbf{A}^{-1} \mathbf{b}_i.
  • (II) Compute the LU factorization of A\mathbf{A} first, and then solve Lyi=bi\mathbf{L} \mathbf{y}_i = \mathbf{b}_i for yi\mathbf{y}_i. Solve Uxi=yi\mathbf{U} \mathbf{x}_i = \mathbf{y}_i for xi\mathbf{x}_i lastly.

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

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

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

下面的算法 P 从输入 A\mathbf{A} 中就地计算 L\mathbf{L}U\mathbf{U},以这样的方式,A\mathbf{A} 的严格下三角部分和上三角部分将分别是 L\mathbf{L}U\mathbf{U}。换句话说,当它终止时,Ai,jA_{i,j} 将为 Li,jL_{i,j}(当 i>ji > j)和 Ui,jU_{i,j}(当 iji \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 分解。

(1100111002210033)\begin{pmatrix} 1 & 1 & 0 & 0 \\ -1 & 1 & 1 & 0 \\ 0 & -2 & 2 & 1 \\ 0 & 0 & -3 & 3 \end{pmatrix}

(2) 证明 Li,j=0L_{i,j} = 0Ui,j=0U_{i,j} = 0 如果 ij>w|i - j| > w

(3) 假设 n1n \gg 1wnw \ll n,并且 mnm \gg n。我们希望求解一组线性系统,

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

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

  • (I) 首先计算 A1\mathbf{A}^{-1},然后将每个 xi\mathbf{x}_i 计算为 A1bi\mathbf{A}^{-1} \mathbf{b}_i
  • (II) 首先计算 A\mathbf{A} 的 LU 分解,然后解 Lyi=bi\mathbf{L} \mathbf{y}_i = \mathbf{b}_i 得到 yi\mathbf{y}_i。最后解 Uxi=yi\mathbf{U} \mathbf{x}_i = \mathbf{y}_i

题目描述

对矩阵 MM,以 Mi:j,k:lM_{i:j,k:l} 表示由第 iijj 行、第 kkll 列组成的子矩阵;Mi,k:lM_{i,k:l}Mi:j,kM_{i:j,k} 分别表示相应行向量、列向量。

AA 是可作 LU 分解的 n×nn\times n 非奇异带状矩阵:存在正整数 ww,使 ij>w|i-j|>wAi,j=0A_{i,j}=0。记 A=LUA=LU,其中 LL 为单位下三角矩阵,UU 为上三角矩阵。题中原地算法 P 对 k=1,,n1k=1,\ldots,n-1 依次执行

Ak+1:n,kAk+1:n,k/Ak,k,A_{k+1:n,k}\leftarrow A_{k+1:n,k}/A_{k,k},
Ak+1:n,k+1:nAk+1:n,k+1:nAk+1:n,kAk,k+1:n,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},

终止后 AA 的严格下三角部分存放 LL 的相应元素,上三角部分存放 UU

(1)计算矩阵

(1100111002210033)\begin{pmatrix} 1&1&0&0\\ -1&1&1&0\\ 0&-2&2&1\\ 0&0&-3&3 \end{pmatrix}

的 LU 分解。

(2)证明当 ij>w|i-j|>w 时,Li,j=0L_{i,j}=0Ui,j=0U_{i,j}=0

(3)设 n1n\gg1wnw\ll nmnm\gg n,需要求解 Axi=bi (i=1,,m)Ax_i=b_i\ (i=1,\ldots,m)。从算术运算量比较:

  • 先求 A1A^{-1},再计算每个 xi=A1bix_i=A^{-1}b_i
  • 先对 AA 作 LU 分解,再分别前代求 Lyi=biLy_i=b_i、回代求 Uxi=yiUx_i=y_i

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

Kai

(1)

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

Step-by-step computation:

  1. Initial matrix A\mathbf{A}:
A=(1100111002210033)\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=1k = 1):
    • Update A[2:4,1]A[2:4,1] by dividing by A[1,1]A[1,1]:
A[2:4,1]A[2:4,1]×(1/A[1,1])=(100)×1=(100)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]:
A[2:4,2:4]A[2:4,2:4]A[2:4,1]×A[1,2:4]=(110221033)(100)×(100)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]=(210221033) A[2:4,2:4] = \begin{pmatrix} 2 & 1 & 0 \\ -2 & 2 & 1 \\ 0 & -3 & 3 \end{pmatrix}
  1. Second iteration (k=2k = 2):
    • Update A[3:4,2]A[3:4,2] by dividing by A[2,2]A[2,2]:
A[3:4,2]A[3:4,2]×(1/A[2,2])=(20)×12=(10)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]:
A[3:4,3:4]A[3:4,3:4]A[3:4,2]×A[2,3:4]=(2133)(10)×(10)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]=(3133)A[3:4,3:4] = \begin{pmatrix} 3 & 1 \\ -3 & 3 \end{pmatrix}
  1. Third iteration (k=3k = 3):
    • Update A[4,3]A[4,3] by dividing by A[3,3]A[3,3]:
A[4,3]A[4,3]×(1/A[3,3])=3×13=1 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]:
A[4,4]A[4,4]A[4,3]×A[3,4]=3(1)×1=4 A[4,4] \leftarrow A[4,4] - A[4,3] \times A[3,4] = 3 - (-1) \times 1 = 4

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

L=(1000110001100011),U=(1100021000310004)\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)

Induct on the elimination step kk. Before step kk, the active matrix is banded. Hence

Ai,k=0(i>k+w),Ak,j=0(j>k+w).A_{i,k}=0\quad(i>k+w),\qquad A_{k,j}=0\quad(j>k+w).

The multipliers stored in column kk therefore vanish below row k+wk+w. Moreover, the rank-one update Ai,jAi,jAi,kAk,jA_{i,j}\leftarrow A_{i,j}-A_{i,k}A_{k,j} can be nonzero only when k<i,jk+wk<i,j\le k+w; such indices satisfy ij<w|i-j|<w. Thus the update creates no entry outside the band. Induction over kk proves Li,j=Ui,j=0L_{i,j}=U_{i,j}=0 whenever ij>w|i-j|>w.

(3)

By (2), one elimination step updates at most w2w^2 entries, so the LU factorization costs O(nw2)O(nw^2). A forward/back substitution with the banded factors costs O(nw)O(nw).

In method (I), forming the generally dense inverse by solving for its nn columns costs O(n2w)O(n^2w), and each multiplication A1biA^{-1}b_i costs O(n2)O(n^2). Its total cost is

O(nw2+n2w+mn2)=O(mn2).O(nw^2+n^2w+mn^2)=O(mn^2).

Method (II) costs

O(nw2+mnw)=O(mnw).O(nw^2+mnw)=O(mnw).

Since wnw\ll n and mnm\gg n, method (II) uses asymptotically fewer arithmetic operations, by a factor of order n/wn/w in the dominant term. It also stores only the banded LU factors instead of a generally dense inverse and avoids explicitly forming A1A^{-1}, which is preferable for storage and numerical accuracy.

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.