作业: 数值线性代数 03

Problem 3

设为可逆矩阵, 考虑如何利用 Gauss 消去法或 LU 分解计算的逆矩阵, 并实现相应算法. 将计算结果与MATLAB、Octave、北太天元等软件中的内置线性方程求解函数或 inv 函数进行比较, 以验证算法的正确性.

Proof

算法的实现:

import numpy as np


def lu(A):
    A = A.astype(float)
    n = A.shape[0]
    L = np.eye(n)
    U = A.copy()
    for k in range(n - 1):
        for i in range(k + 1, n):
            L[i, k] = U[i, k] / U[k, k]
            U[i, k:] -= L[i, k] * U[k, k:]
    return L, U


def f_calc(L, b):
    n = len(b)
    x = np.zeros(n)
    for i in range(n):
        x[i] = b[i] - L[i, :i] @ x[:i]
    return x


def b_calc(U, b):
    n = len(b)
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - U[i, i + 1:] @ x[i + 1:]) / U[i, i]
    return x


def inverse_mat(A):
    n = A.shape[0]
    L, U = lu(A)
    A_inv = np.zeros_like(A)
    I = np.eye(n)
    for j in range(n):
        y = f_calc(L, I[:, j])
        x = b_calc(U, y)
        A_inv[:, j] = x
    return A_inv

测试

A = np.array([
    [4, 2, 1],
    [2, 5, 2],
    [1, 2, 4]
], dtype=float)

A_inv = inverse_mat(A)

print("LU inv:")
print(A_inv)

print("NP inv:")
print(np.linalg.inv(A))

print("Err:")
print(np.linalg.norm(A_inv - np.linalg.inv(A)))

print("A @ A_inv:")
print(A @ A_inv)

结果

LU inv:
[[ 0.31372549 -0.11764706 -0.01960784]
 [-0.11764706  0.29411765 -0.11764706]
 [-0.01960784 -0.11764706  0.31372549]]
NP inv:
[[ 0.31372549 -0.11764706 -0.01960784]
 [-0.11764706  0.29411765 -0.11764706]
 [-0.01960784 -0.11764706  0.31372549]]
Err:
0.0
A @ A_inv:
[[1.00000000e+00 2.77555756e-17 0.00000000e+00]
 [0.00000000e+00 1.00000000e+00 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 1.00000000e+00]]

误差很小, 推测是由于实现不同导致的浮点数误差. Numpy 的实现总体也基于LU分解, 但会进行一定的主元选择, 且底层计算过程与上述实现不同. 上述的实现没有包含对错误输入的检查与报错, 且假定了输入矩阵为浮点数矩阵. 在输入是浮点数可逆矩阵, 且不选主元消去能够顺利进行的条件下计算结果是正确的.

Problem 6

如果一个矩阵满足

则称其为带宽为的带状矩阵.

假设可以不选主元直接计算的 LU 分解, 则矩阵和矩阵的稀疏性分别有什么特点?
假设为不依赖于的常数, 证明不选主元的 LU 分解的计算复杂度为 .
如果计算的列主元 LU 分解, 则矩阵和矩阵的稀疏性又分别有什么特点?

Proof

同样保持的带宽, 即满足, 满足. 这是由于在消去的过程中, 第行只会消去行, 不会影响行在第列上的. 同样, 第行也只会被行消去, 不会影响右侧的.
每一行进行消去时, 需要计算次消去系数, 每次消去只需要更新个值, 因此复杂度为. 而一共有行, 因此总体的复杂度为. 为常数, 因此复杂度为.
在列主元 LU 分解的过程中, 最坏情况为第行与行交换, 因此会使的带宽扩大两倍. 即满足. 而的元素同时也由置换矩阵的乘积决定, 因此不一定会保持带状结构, 但置换矩阵不会改变非零元的个数.

Problem 8

设 , 且是可逆矩阵. 描述一个使用的 LU 分解求解如下线性方程组的快速方法:

Proof

. 先考虑. 记, 先对用前代法求出. 然后对用回代法求出.

然后考虑. 有, 记, . 先对用前代法求出. 然后对用回代法求出.

Problem 10

在求解某个计算问题时, 如果存在数据更新, 我们自然希望能够充分利用更新前问题的解, 而不是完全从头开始求解更新后的问题. 考虑线性方程组

其中非奇异, . 当时, 可以看作是的低秩更新.

设. 计算矩阵的分块 LU 分解以及分块 UL 分解.
利用 1. 证明非奇异当且仅当非奇异, 并进一步证明 Sherman-Morrison-Woodbury 公式
假设已经有了矩阵的 LU 分解, 利用 Sherman-Morrison-Woodbury 公式设计一个求解线性方程组

且运算量为的算法. 相比之下, 重新计算矩阵的 LU 分解的运算量为 .

Proof

第一行乘以消去第二行

同样, 使用第二行消去第一行
考虑该分块矩阵的LU分解, 非奇异等价于非奇异. 考虑其UL分解, 非奇异等价于非奇异. 因此非奇异当且仅当非奇异.

记. 考虑的逆

考虑有
即求解
1. 对于, 利用LU分解可以求得的解, 复杂度为.
2. 计算即求解. 考虑, . 即求解个方程, 复杂度为.
3. 计算, 由于已有, 通过矩阵乘法求得结果的复杂度为.
4. 计算, 由于已有, 复杂度为.
5. 通过求. 使用LU求解的复杂度为.
6. , 复杂度为
因此总体的复杂度为

, 因此复杂度为.

Chapter 2