2021-02-28-数值分析-Day03-范数-迭代改善算法

2.9-11 范数

为研究线性方程组近似解的误差估计和迭代法的收敛性，对 \(R……n\) （ \(n\) 维向量空间）中的向量（ \(R^{n\times n}\) 中的矩阵）的“大小”引进度量——向量（或矩阵）的范数。

向量的范数

定义1：设 \(x=(x_1,x_2,\dots,x_n)^T,y=(y_1,y_2,\dots,y_n)^T\in R^n (or\ C^n)\)
将实数 \((x,y)=y^Tx=\sum_{i=1}^{n}x_i\overline{y_i}\)（或复数 \((x,y)=y^Hx=\sum_{i=1}^{n}x_iy_i\) ）称为向量x,y的数量积（内积）。
将非负实数 \(\|x\|_2=(x,x)^{1\over2}=(\sum_{i=1}^{n}x_i^2)^{1\over2}\) 或 \(\|x\|_2=(x,x)^{1\over2}=(\sum_{i=1}^{n}|x_i|^2)^{1\over2}\) 称为向量 \(x\) 的欧式范数。

引理：向量的内积运算支持交换律、分配律、与实数乘时的结合律。
Cauchy-Schwarz不等式 \(|(x,y)|\le \|x\|_2\cdot\|y\|_2\) 当且仅当 \(x\) 与 \(y\) 线性相关时取等。
三角不等式 \(\|x+y\|_2\le\|x\|_2+\|y\|_2\)

定义2-向量的范数：若向量 \(x\in R^n\) （或 \(C^n\) ）的某个实值函数 \(N(x)=\|x\|\) ，满足条件： \[ \begin{aligned} & (1)\ \|x\|\ge 0\ (\|x\|=0当且仅当x=0)\quad\textbf{正定条件}\\ & (2)\ \|\alpha x\|=|\alpha|\|x\|,\ \forall\alpha\in R(或\alpha\in C)\quad\textbf{齐次性}\\ & (3)\ \|x+y\|_2\le\|x\|_2+\|y\|_2\quad\textbf{三角不等式}\\ & (3\to 4)\ \|x-y\|_2\ge\|x\|_2-\|y\|_2 \end{aligned} \] 则称 \(N(x)\) 是 \(R^n\) （或 \(C^n\) ）上的一个向量范数（或模）

几种常用的向量范数：

向量的 \(\infty\) -范数（最大范数）： \(\displaystyle{\|x\|_\infty=\max\limits_{1\le i\le n}|x_i|}\)
向量的 \(1\) -范数： \(\displaystyle{\|x\|_1=\sum_{i=1}^{n}|x_i|}\)
向量的 \(2\) -范数： \(\displaystyle{\|x\|_2=(x,x)^{1\over2}=(\sum_{i=1}^{n}x_i^2)^{1\over2}}\)
向量的 \(p\) -范数： \(\displaystyle{\|x\|_p=(\sum_{i=1}^{n}|x_i|^p)^{1\over p}}\)
向量的 \(\infty\) -范数（最大范数）： \(\displaystyle{\|x\|_{-\infty}=\min\limits_{1\le i\le n}|x_i|}\)
向量的 \(0\) -范数：向量的非零元素个数

向量的收敛性：设 \(\{x^{(k)}\}\) 为 \(R^n\) 中的一个向量序列， \(x^*\in R^n\) ，记 \(x^{(k)}=(x_1^{(k)}=(x_1^{(k)},x_2^{(k)},\dots,x_n^{(k)})^T\, x^*=(x_1^*,x_2^*,\dots,x_n^*)^T\) 。若 \(\lim\limits_{k\to\infty}x_i^{(k)}=x_i^*(i=1,2,\dots,n)\) ，则称 \(x^{(k)}\) 收敛于 \(x^*\) ，记为 \(\lim\limits_{k\to\infty}x_i^{(k)}=x^*\) 。

范数的连续性：设非负函数 \(N(x)=\|x\|\) 为 \(R^n\) 上任一向量范数，则 \(N(x)\) 是 \(x\) 的分量 \(x_1,x_2,\dots,x_n\) 的连续函数。 \(x\) 有很小变化时，同时 \(N(x)\) 亦变化不大。

向量范数的等价性：设 \(\|x\|_s,\|x\|_t\) 为 \(R^n\) 上向量的任意两种范数，则存在常数 \(c_1,c_2>0\) ，使得对一切 \(x\in R^n\) 有 \(c_1\|x\|_s\le\|x\|_t\le c_2\|x\|_s\) 。

矩阵的范数

方阵

定义1：若矩阵 \(A\in R^{n\times n}\) 的某个非负实值函数 \(N(A)=\|A\|\) 满足条件：
\[ \begin{aligned} & (1)\ \|A\|\ge 0\ (\|A\|=0\Leftrightarrow A=0)\quad\textbf{正定条件}\\ & (2)\ \|c A\|=|c|\|A\|,\ c为实数\quad\textbf{齐次性}\\ & (3)\ \|A+B\|\le\|A\|+\|B\|\quad\textbf{三角不等式}\\ & (4)\ \|AB\|\le\|A\|\|B\| \end{aligned} \] 则称 \(N(A)\) 是 \(R^{n\times n}\) 上的一个矩阵范数（或模）。

定义2-矩阵的算子范数：设 \(x\in R^n, A\in R^{n\times n}\) ，给出一种向量范数 \(\|x\|_\nu,\ (如\nu=1,2,\infty,\dots)\) ，相应地定义一个矩阵的非负函数： \(\displaystyle{\|A\|_\nu=\max\limits_{x\neq0}\frac{\|Ax\|_\nu}{\|x\|_\nu}}\) ，可验证该非负函数为矩阵 \(A\) 的范数，称之为矩阵 \(A\) 的算子范数。

定理： \(\|Ax\|_\nu\le\|A\|_\nu\|x\|_\nu\)

几种常用的矩阵范数：

Frobenius范数（F范数）： \(\displaystyle{F(A)=\|A\|_F=\left(\sum\limits_{i,j=1}^{n}a_{i,j}^2\right)^{1\over2}}\) ，由向量的 \(2\) -范数推广而来
行范数（\(\infty\) -范数）： \(\displaystyle{\|A\|_\infty=\max\limits_{1\le i\le n}\sum\limits_{j=1}^{n}|a_{i,j}|}\)
列范数（ \(1\) -范数）： \(\displaystyle{\|A\|_1=\max\limits_{1\le j\le n}\sum\limits_{i=1}^{n}|a_{i,j}|}\)
谱范数（ \(2\) -范数）： \(\displaystyle{\|A\|_2=\sqrt{\lambda_{\max}(A^TA)},\ (\lambda为特征值)}\)

非方阵

矩阵 \(A_{mn}\)

Frobenius范数（F范数）： \(\displaystyle{F(A)=\|A\|_F=\left(\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{n}a_{i,j}^2\right)^{1\over2}}\) ，由向量的 \(2\) -范数推广而来
行范数（\(\infty\) -范数）： \(\displaystyle{\|A\|_\infty=\max\limits_{1\le i\le n}\sum\limits_{j=1}^{n}|a_{i,j}|}\)
列范数（ \(1\) -范数）： \(\displaystyle{\|A\|_1=\max\limits_{1\le j\le n}\sum\limits_{i=1}^{m}|a_{i,j}|}\)
谱范数（ \(2\) -范数）： \(\displaystyle{\|A\|_2=\sqrt{\lambda_{\max}(A^TA)},\ (\lambda为特征值)}\)

2.12 谱半径

特征值和特征向量 \(Ax=\lambda x\) 常数 \(\lambda\) 为特征值，矩阵 \(A\) 为特征向量，且为方阵。计算方法： \(\det(A-\lambda I)=0\)

定义：设 \(A\in R^{n\times n}\) 的特征值为 \(\lambda_i(i=1,2,\dots,n)\) ，称 \(\rho(A)=\max\limits_{1\le i\le n}|\lambda_i|\) 为 \(A\) 的谱半径。

定理：

设 \(A\in R^{n\times n}\) ，则 \(\rho(A)\le\|A\|\) ，即矩阵的谱半径不超过矩阵的任何一种算子范数（对F范数也成立）。
若 \(A\in R^{n\times n}\) 为对称矩阵，则 \(\rho(A)=\|A\|_2\) 。
若 \(\|B\|<1\) ，则 \(I+B\) 为非奇异矩阵，且 \(\left\|(I\pm B)^{-1}\right\|\le\displaystyle\frac{1}{1-\|B\|}\) ，此处 \(\|\ \|\) 指算子范数。

2.13 - 14 方程组的条件数

定义：如果矩阵 \(A\) 或常数项 \(b\) 的微小变化，引起方程组 \(Ax=b\) 的解的巨大变化，则称该方程组为“病态”方程组，若矩阵 \(A\) 为原因则称之为该方程组的“病态”矩阵。否则称“良态”。

设 \(A\) 精确， \(b\) 有误差 \(\delta b\) ，解为 \(x+\delta x\) ，则： \[ A(x+\delta x)=b+\delta b\ \Rightarrow\ \delta x=A^{-1}\delta b\\ \|\delta x\|\le\|A^{-1}\|\|\delta b\|\\ \] 由 \(Ax=b\) 得： \[ \|b\|\le\|A\|\|x\|,\ \frac{1}{\|x\|}\le\frac{\|A\|}{\|b\|}\ (b\neq0) \] 两式相乘得到定理1：设 \(A\) 是非奇异矩阵， \(Ax=b\neq0\) ，且 \(A(x+\delta x)=b+\delta b\) ，则 \(\displaystyle\frac{\|\delta x\|}{\|x\|}\le\|A^{-1}\|\|A\|\frac{\|\delta b\|}{\|b\|}\) 。该定理给出解的相对误差的上街，即常数项 \(b\) 的相对误差在解中至多放大 \(\|A^{-1}\|\|A\|\) 倍。
设 \(b\) 精确， \(A\) 有微小误差（扰动） \(\delta A\) ，解为 \(x+\delta x\) ，则： \[ (A+\delta A)(x+\delta x)=b\ \Rightarrow\ (A+\delta A)\delta x=-(\delta A)x\\ \] 若 \(\delta A\) 不受限制， \(A+\delta A\) 可能奇异，而： \[ (A+\delta A)=A(I+A^{-1}\delta A) \] 当 \(\|A^{-1}\delta A\|<1\) 时， \((I+A^{-1}\delta A)^{-1}\) 存在。则可得： \[ \delta x=-(I+A^{-1}\delta A)^{-1}A^{-1}(\delta A)x\\ \|\delta x\|\le\frac{\|A^{-1}\|\|\delta A\|\|x\|}{1-\|A^{-1}(\delta A)\|} \] 设 \(1-\|A^{-1}(\delta A)\|<1\) ，则有： \[ \frac{\|\delta x\|}{\|x\|}\le\frac{\|A^{-1}\|\|A\|\frac{\|\delta A\|}{\|A\|}}{1-\|A^{-1}\|\|A\|\frac{\|\delta A\|}{\|A\|}} \]

综上， \(\|A^{-1}\|\|A\|\) 刻画了方程组的“病态”程度，称之为矩阵 \(A\) 的条件数： \[ cond(A)_\nu=\|A^{-1}\|_\nu\|A\|_\nu\ (\nu=1,2,\infty) \] \(cond(A)\) 过大时，方程组“病态”。（相对于矩阵中的元比较大小。）

通常使用的条件数： + 无穷条件数 \(cond(A)_\infty=\|A^{-1}\|_\infty\|A\|_\infty\) + 谱条件数 \(cond(A)_2=\|A^{-1}\|_2\|A\|_2=\displaystyle\sqrt{\frac{\lambda_{\max}(A^TA)}{\lambda_{\min}(A^TA)}}\) 当 \(A\) 为对称矩阵时， \(cond(A)_2=\displaystyle\frac{\lambda_1}{\lambda_2}\) ，其中 \(\lambda_1,\lambda_2\) 为 \(A\) 的绝对值最大和绝对值最小的特征值。

条件数的性质： + \(cond(A)_\nu=\|A^{-1}\|_\nu\|A\|_\nu\ge\|A^{-1}A\|=1\) + \(cond(cA)_\nu=cond(A)_\nu, c为常数\) + 若 \(A\) 为正交矩阵（ \(A^{-1}=A^T\) ），则 \(cond(A)_2=1\) 若 \(A\) 为非奇异矩阵， \(R\) 为正交矩阵，则 \(cond(RA)_2=cond(AR)_2=cond(A)_2\)

2.15 事后误差估计与迭代改善算法

事后误差估计

设 \(A\) 为非奇异矩阵， \(x\) 是 \(Ax=b\neq0\) 的精确解。再设 \(\overline x\) 是此方程组的近似解， \(r=b-A\overline x\) ，则： \[ \frac{\|x-\overline x\|}{\|x\|}\le cond(A)\cdot\frac{\|r\|}{\|b\|} \]

迭代改善算法

理论：设 \(Ax=b\) ，其中 \(A\in R^{n\times n}\) 为非奇异矩阵，且为病态方程组（但不过分病态）。 + 用选主元三角分析法实现分解计算 \(PA=LU\) ，其中 \(P\) 为置换阵， \(L\) 为单位下三角阵， \(U\) 为上三角阵，求得方程组近似解 \(x_1\) 。 + 计算剩余向量： \(r_1=b-Ax_1\) + 求解 \(Ad=r_1\) ，得到解 \(d_1\) ，然后改善原解： \(x_2=x_1+d_1\) 。 + 递归推：\(d_2,d_3,\dots,\ x_3,x_4,\dots\)

实现： + 用选主元三角分解实行分解计算 \(PA=LU\) ，且求计算解 \(x_1\) 。 + 对于 \(k=1,2,\dots,N_0\) : + 计算 \(r_k=b-Ax_k\) + 求解 \(LUd_k=Pr_k=PAd_k\) ，\(\left\{\begin{array}{}Ly=Pr_k\\Ud_k=y\end{array}\right.\) + 若 \(\displaystyle\frac{\|d_k\|_\infty}{\|x_k\|_\infty}\le10^{-t}\) 则输出 + 改善 \(x_{k+1}=x_k+d_k\)

补充：选主元三角分解法

LU分解过程中的问题：有解问题无法LU分解（主元为0）；主元过小时引入较大误差

解决方法：选主元

第 \(k\) 列主元及主元所在列下方元素挑选绝对值最大的数
通过初等变换使该元位于主对角线上，最终将 \(A\) 转换为 \(PA\)
则 Doolittle分解转换为 \(PA=LU\)

代码：迭代改善算法

//
// Created by xa on 2021-02-28.
//
#include <iostream>
#include <vector>
#include <algorithm>

using std::vector;

vector<vector<double>> A;
vector<vector<double>> C;
vector<double> b;
vector<double> x;
vector<double> r;

int size = 0;

void iterate(vector<vector<double>> A_, vector<double> b_);
void solve();

void getPA();
void getLU();
void refreshR();

template<typename T>
void swap(T& a, T& b);

void printX();

int main(){
    vector<vector<double>> A_{{1.0303, 0.99030}, {0.99030, 0.95285}};
    vector<double> b_{2.4944, 2.3988};
    iterate(A_, b_);
    for (int i = 0; i < size; ++i)
        std::cout << x[i] << " ";
}

void iterate(vector<vector<double>> A_, vector<double> b_) {
    A.assign(A_.begin(), A_.end());
    b.assign(b_.begin(), b_.end());
    size = A.size();
    getPA();
    C.assign(A.begin(), A.end());
    r.assign(b.begin(), b.end());
    getLU();
    solve();
    x.assign(r.begin(), r.end());

    while (true) {
        refreshR();
        solve();
        if (*(std::max_element(std::begin(r), std::end(r))) /
            *(std::max_element(std::begin(x), std::end(x))) <= 1e-5)
            break;
        for (int i = 0; i < size; ++i)
            x[i] += r[i];
        printX();
    }
}

void solve() {
    vector<double> d(size);
    vector<double> y(size);
    // solve y
    y[0] = r[0];
    for (int k = 1; k < size; ++k) {
        y[k] = r[k];
        for (int i = 0; i < k; ++i)
            y[k] -= C[k][i] * y[i];
    }
    // solve d
    d[size-1] = y[size-1] / C[size-1][size-1];
    for (int i = size-2; i >= 0; --i) {
        double sum = 0;
        for (int j = i+1; j < size; ++j)
            sum += C[i][j] * r[j];
        d[i] = (y[i] - sum) / C[i][i];
    }
    r.assign(d.begin(), d.end());
}

void getPA() {
    for (int k = 0; k < size; ++k) {
        int max = k;
        if (k < size-1) {
            for (int i = k; i < size; ++i)
                if (A[i][k] > A[max][k]) max = i;
            swap(A[k], A[max]);
            swap(b[k],b[max]);
        }
    }
}

void getLU() {
    for (int k = 0; k < size; ++k) {
        if (k > 0) for (int j = k; j < size; ++j)
                for (int m = 0; m < k; ++m) {
                    C[k][j] -= C[k][m] * C[m][j];
                }
        if (k < size-1) for (int i = k+1; i < size; ++i) {
            double sum = 0;
            for (int m = 0; m < k; ++m)
                sum += C[i][m] * A[m][k];
            C[i][k] = (C[i][k] - sum) / C[k][k];
        }
    }
}

void refreshR() {
    for (int i = 0; i < size; ++i) {
        r[i] = b[i];
        for (int j = 0; j < size; ++j)
            r[i] -= A[i][j] * x[j];
    }
}

template<typename T>
void swap(T& a, T& b) {
    double tmp = a;
    a = b;
    b = tmp;
}

void printX() {
    for (int i = 0; i < size; ++i) {
        std::cout << x[i] << ' ';
    }
    std::cout << std::endl;
}

代码：选主元Doolittle算法

//
// Created by xa on 2021-02-28.
//
#include <iostream>
#include <vector>

using std::vector;

vector<vector<double>> A;
vector<double> b;

int size = 0;

vector<double> solve(vector<vector<double>> A_, vector<double> b_);

void getPA();
void getLU();

template<typename T>
void swap(T& a, T& b);

void printA();

int main(){
    vector<vector<double>> A_ {{8, 1, 7}, {3, 7, 9}, {9, 1, 5}};
    vector<double> b_ {3, 8, 6};
    vector<double> x = solve(A_, b_);

    for (int i = 0; i < size; ++i)
        std::cout << x[i] << " ";
}

vector<double> solve(vector<vector<double>> A_, vector<double> b_) {
    A.assign(A_.begin(), A_.end());
    b.assign(b_.begin(), b_.end());
    size = A.size();
    vector<double> x(size);
    vector<double> y(size);

    getPA();
    getLU();

    // solve y
    y[0] = b[0];
    for (int k = 1; k < size; ++k) {
        y[k] = b[k];
        for (int i = 0; i < k; ++i)
            y[k] -= A[k][i] * y[i];
    }

    // solve x
    x[size-1] = y[size-1] / A[size-1][size-1];
    for (int i = size-2; i >= 0; --i) {
        double sum = 0;
        for (int j = i+1; j < size; ++j)
            sum += A[i][j] * x[j];
        x[i] = (y[i] - sum) / A[i][i];
    }
    return x;
}

void getPA() {
    for (int k = 0; k < size; ++k) {
        int max = k;
        if (k < size-1) {
            for (int i = k; i < size; ++i)
                if (A[i][k] > A[max][k]) max = i;
            swap(A[k], A[max]);
            swap(b[k],b[max]);
        }
    }
    printA();
}

void getLU() {
    for (int k = 0; k < size; ++k) {
        if (k > 0) for (int j = k; j < size; ++j)
                for (int m = 0; m < k; ++m) {
                    A[k][j] -= A[k][m] * A[m][j];
                }
        if (k < size-1) for (int i = k+1; i < size; ++i) {
            double sum = 0;
            for (int m = 0; m < k; ++m)
                sum += A[i][m] * A[m][k];
            A[i][k] = (A[i][k] - sum) / A[k][k];
        }
    }
    printA();
}

template<typename T>
void swap(T& a, T& b) {
    double tmp = a;
    a = b;
    b = tmp;
}

void printA() {
    for (int i = 0; i < size; ++i) {
        for (int j = 0; j < size; ++j)
            std::cout << A[i][j] << ' ';
        std::cout << std::endl;
    }
    std::cout << std::endl;
}

2.9-11 范数¶

向量的范数¶

矩阵的范数¶

方阵¶

非方阵¶

2.12 谱半径¶

2.13 - 14 方程组的条件数¶

2.15 事后误差估计与迭代改善算法¶

事后误差估计¶

迭代改善算法¶

补充：选主元三角分解法¶

代码：迭代改善算法¶

代码：选主元Doolittle算法¶