MyMath数学库(三) - 矩阵操作

我的无人车之路 mymath

---

$\quad$矩阵操作使用了著名的C++矩阵库Eigen在其基础之上封装了模板函数，主要是提供两个功能，一个是求矩阵的伪逆矩阵，一个是计算状态空间下连续形式到离散形式的线性变换。源文件在matrix_operations.hpp中，直接包含即可使用，但是在编译程序时，要使用-I参数包含Eigen的路径。
$\quad$虽然此库只涉及了两个函数，但是确蕴含了很丰富的线性代数与优化知识，在介绍完模块的基本使用后，后边小节补充了一些关于矩阵的逆以及双线性变换的知识，不是非常系统。其中涉及的几个知识点都有讨论。

3.1 函数表

函数名	说明
pseudo_inverse(方阵)	求一方阵的Moore-Penrose伪逆。
pseudo_inverse(任意矩阵)	求一任意矩阵的Moore-Penrose伪逆。
continuous_to_discrete(类型模板)	计算状态空间表示形式的连续形式到离散形式的双线性变换。
continuous_to_discrete(以double类型基础类型的矩阵Matrixd为参数)	计算状态空间表示形式的连续形式到离散形式的双线性变换。

3.2 使用例子

3.3 Eigen开源库

$\quad$Eigen是一个高层次的C ++库，有效支持线性代数，矩阵和矢量运算，数值分析及其相关的算法。Eigen是一个开源库。这里是官网地址。
$\quad$不用编译直接包含即可使用，这里就不做太多叙述了。强力推荐。

3.4 行列式求解

3.5 $n \times n$矩阵的逆

$\quad$若存在一个矩阵$B$使得$AB=BA=I$，则称$n \times n$矩阵$A$为非奇异的或可逆的。矩阵$B$称为$A$的乘法逆元。其中$I=\begin{bmatrix} 1 & 0 &\cdots &0\\0 & 1 & \cdots & 0 \\ \vdots & & \\ 0 & \cdots & & 1 \end{bmatrix}$。我们将非奇异矩阵$A$的乘法逆元简称为$A$的逆元，并记作$A^{-1}$。
$\quad$一个矩阵有逆，首先条件是一个方阵，其次每个向量都是线性无关的。这里列出线性方程组的一些情况。

1. $n$个方程，$m$个未知数，$n=m$具备有唯一解。(非奇异方阵)
2. $n$个方程, $m$个未知数，$n<m$方程有无数解。(亚定方程组)
3. $n$个方程, $m$个未知数，$n>m$方程可能有唯一解（$m-n$个线性相关的方程）,也可能无解，只能求最小二乘解。(超定方程组)

3.5.1 非奇异方阵的求逆算法

$\quad$设非奇异方阵$A$，由于$A$与单位矩阵$I$是行等价的。所以可以通过一系列的线性初等变换得到$I$。

$$A=E_kE_{k-1} \cdots E_1I=E_kE_{k-1} \cdots E_1$$
由于$E_i(i=1,\ldots,k)$是可逆的，其乘积也是可逆的，则$A^{-1}=(E_kE_{k-1} \cdots E_1)^{-1}=E_1^1 \cdots E_{k-1}^{-1}E_k^{-1}$。
$\quad$在具体实施方面可以通过$A$的增广矩阵$(A|I)$求行最简形$(I|A^{-1})$得到$A$的逆矩阵。矩阵的行初等变换可以通过高斯消元算法实现。具体算法定义如下：

【算法】高斯消元法

$\quad$在算法进行之前首先定义两个辅助算法，一个是对矩阵进行排序算法，其次检查矩阵是否为行阶最简型号。

矩阵排序算法

$\quad$假设函数名为matrix_sort(M)。采用什么样的排序算法都可以，当两个向量比较时自定义算法，从第$1$列到第$n$列的元素，一个一个元素进行检测直到遇到一个非$0$元素则计数完毕，计数越多则表示当前向量越大，全$0$向量最大排在矩阵的最后一列。

检查是否行阶段梯形矩阵

$\quad$假设函数名为is_row_echelon_form(M)在调用此函数之前，首先使用matrix_sort对矩阵进行排序。因为有排序算法的原因，所以从最后一行开始检测，如果是全$0$行则检测它的上一行，如果不是全$0$行则检测当前行的倒数第二列是否是$0$，如果是$0$则说明是行最简型的。

行最简行型算法

目标矩阵为一个$i \times j$的矩阵$M$。

1. 使用matrix_sort(M)函数对目标矩阵$M$进行排序。目的是将首变量为非$0$的行排在前面。
2. 使用is_row_echelon_form(M)判断当前矩阵$M$是否是一个行阶段梯形矩阵，如果是则退出，不是则继续流程3。
3. 从第$n$行开始，选取$a_{n}$作为首变量。排序算法保证了这里的首变量是此行的第$n$个元素。
4. 遍历其余的$n+(j-n)$行，并且选取$a_{n+(j-n)}$作为要操作的元素，检测是否为$0$。
5. 如果$a_{n+(j-n)}$为$0$则跳过，并跳过剩下的行。因为排序的原因，在它之后的所有行数对应列的元素都应该为$0$。
6. 如果$a_{n+(j-n)}$不为$0$，则使用$\frac{a_{n+(j-n)}}{a_n}$乘以第$n$行的结果减去第$n$行的结果作为新的第$n$行。
7. 这里由于数值计算的原因，可能有些向量元素只是趋近于$0$，将这些趋近于$0$的值都当作$0$来处理。
8. $n = n+1$。
9. 回到流程1。

3.6 矩阵的值域与零空间

3.7 矩阵的正交空间

3.8 矩阵的最小二乘解

$\quad$一个超定线性方程组的解一般可以化为最小二乘问题。给定一个$m \times n$的方程组$Ax=b$，其中$m \gt n$，一般我们不能期望找到一个向量$x \in \mathbf R^n$，使得$Ax$等于$b$。事实上，可以寻找一个向量$x$，使得$Ax$最接近$b$。
$\quad$给定一个方程组$Ax=b$，其中$A$为一个$m\times n$($m \gt n$)矩阵，并且$b \in \mathbf R^m$，则对每一$x \in \mathbf R^n$，可以构造一个残差。

$$r(x)=b-Ax$$ 其中$b$和$Ax$间的距离为 $$||b-Ax||=||r(x)||$$ 我们希望寻找一个向量$x \in \mathbf R^n$，使得$||r(x)||$是最小的。最小化$||r(x)||$等价于最小化$||r(x)||^2$。达到最小值的向量$\hat x$称为方程组$Ax=b$的最小二乘解。
$\quad$若$\hat x$为方程组$Ax=b$的最小二乘解，且$p=A \hat x$，则$p$就是$A$的列空间中和$b$最接近的向量。向量$p$称为$b$在$R(A)$上个的投影。 $$b-p = b-A \hat x = r(\hat x)$$ 其中$\hat x$必为$R(A)^{\bot}$（$R(A)$正交补）中的元素。因此$\hat x$为最小二乘问题的解的充要条件是$r(\hat x) \in R(A)^{\bot}$。我们知道$R(A)^{\bot}=N(A^T)$，$N(A^T)$为$A$的零域。这也可说明一个向量$\hat x$是方程组$Ax=b$的最小二乘解，当且仅当$r(\hat x) \in N(A^T)$。残差属于$A$的零域。
$$\mathbf 0 = A^Tr(\hat x) = A^T(b-A \hat x)$$
$\quad$因此，为求解最小二乘问题$Ax=b$，我们必须求解 $$A^TAx=A^Tb$$ 此方程表示的$n \times n$线性方程组成为正规方程组。一般地，正规方程组可能存在多个解；然而，若$\hat x$和$\hat y$均为解，则由于$b$在$R(A)$中的投影$p$是唯一的，故 $$A \hat x=A \hat y = p$$ 有唯一解。其解 $$\hat x = (A^TA)^{-1}A^Tb$$

3.9 矩阵的内积空间

3.10 矩阵的特征值与特征向量

$\quad$考虑方程$Ax=\lambda x$。这个方程出现在很多线性代数的应用问题中。如果方程有非零解$x$，则$\lambda$称为$A$的特征值，且$x$称为属于$\lambda$的特征向量。一般地，我们可以将特征值看成与线性变换相关联的自然频率。若$A$为一$n \times n$矩阵，则可将$A$看成一个$R^n$上的线性算子。
$\quad$很多应用问题都涉及将一个线性变换重复作用到一个向量上。求解这类问题的关键是针对算子选择一个在某种意义下很自然的坐标系或基，并使得包含该算子的计算得以简化。对应于这一组新的基向量（特征向量），我们关联一个缩放因子（特征值）表示该算子的自然频率。
$\quad$我们设一$m \times n$的矩阵$A$, $w_0 \in R^n$，我们反复应用$A$到$x$上。$A^n w_0 = w_n$，最后这个$w_n$稳定下来保持不变，那么这个向量$w_n$成为这个过程的稳态向量。线性算子$A$只是将向量$w_0$缩小或扩大了$n$倍而已，所以将这个倍数提取出来为$\lambda$，向量

3.11 矩阵的特征值分解

3.12 连续到离散的变换

【算法】双线性变换