For a matrix denoted by M, we write Mi: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 Mi:i,k:l and Mi:j,k:k by Mi,k:l and Mi:j,k, respectively. (Note that they are row and column vectors, respectively.)
Assume below that A is an n×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 A are zero if k>w. Namely, if ∣i−j∣>w, Ai,j is zero.
We denote the LU factorization of A by A=LU, where L is a lower triangular matrix with unit diagonal elements, and U is an upper triangular matrix.
The following algorithm P computes L and U from the input A in-place in such a way that the strictly lower triangular part and the upper triangular part of A will be L and U, respectively. In other words, Ai,j will be Li,j for i>j, and Ai,j will be Ui,j for i≤j, when it terminates.
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.
1−10011−20012−30013
(2) Prove that Li,j=0 and Ui,j=0 if ∣i−j∣>w.
(3) Assume that n≫1, w≪n, and m≫n. We wish to solve a set of linear systems,
Axi=bi,i=1,2,…,m,
where bi(i=1,2,…,m) are the constant vectors, and xi(i=1,2,…,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 A−1 first, and then compute each xi as A−1bi.
(II) Compute the LU factorization of A first, and then solve Lyi=bi for yi. Solve Uxi=yi for xi lastly.
对于矩阵 M,我们写 Mi:j,k:l 表示从矩阵的第 i 行到第 j 行以及第 k 列到第 l 列构成的子矩阵。我们还将 Mi:i,k:l 和 Mi:j,k:k 分别表示为 Mi,k:l 和 Mi:j,k。(注意它们分别是行向量和列向量。)
假设 A 是一个 n×n 的非奇异矩阵,可以进行 LU 分解。另外,假设存在一个正整数 w,使得当 ∣i−j∣>w 时,A 在第 k 条下对角线和超对角线上的所有元素都为零。即,如果 ∣i−j∣>w,则 Ai,j 为零。
我们用 A=LU 来表示 A 的 LU 分解,其中 L 是具有单位对角元素的下三角矩阵,U 是上三角矩阵。
下面的算法 P 从输入 A 中就地计算 L 和 U,以这样的方式,A 的严格下三角部分和上三角部分将分别是 L 和 U。换句话说,当它终止时,Ai,j 将为 Li,j(当 i>j)和 Ui,j(当 i≤j)。
Induct on the elimination step k. Before step k, the active matrix is banded. Hence
Ai,k=0(i>k+w),Ak,j=0(j>k+w).
The multipliers stored in column k therefore vanish below row k+w. Moreover, the rank-one update
Ai,j←Ai,j−Ai,kAk,j can be nonzero only when
k<i,j≤k+w; such indices satisfy ∣i−j∣<w. Thus the update creates no entry outside the band. Induction over k proves
Li,j=Ui,j=0 whenever ∣i−j∣>w.
By (2), one elimination step updates at most w2 entries, so the LU factorization costs O(nw2). A forward/back substitution with the banded factors costs O(nw).
In method (I), forming the generally dense inverse by solving for its n columns costs O(n2w), and each multiplication A−1bi costs O(n2). Its total cost is
O(nw2+n2w+mn2)=O(mn2).
Method (II) costs
O(nw2+mnw)=O(mnw).
Since w≪n and m≫n, method (II) uses asymptotically fewer arithmetic operations, by a factor of order n/w in the dominant term.
It also stores only the banded LU factors instead of a generally dense inverse and avoids explicitly forming A−1, which is preferable for storage and numerical accuracy.