矩阵最小二乘法问题求解_最小二乘法求矩阵方程

一、超定方程组

超定方程组是指方程个数大于未知量个数的方程组。对于方程组$Ax=b$，$A$为n×m矩阵，如果R列满秩，且n>m。则方程组没有精确解，此时称方程组为超定方程组。
在实验数据处理和曲线拟合问题中，求解超定方程组非常普遍。比较常用的方法是最小二乘法。
如果有向量$x$使得下式的值达到最小，则称 $x$为上述超定方程的最小二乘解。

二、矩阵形式的最小二乘法

最小二乘法问题：
$$E(\mathbf{x})=\Vert \mathbf{Ax}-\mathbf{b}{\Vert }^2$$求使$E(\mathbf{x})$达到最小值时$x$的值：
E ( x ) = ∥ A x − b ∥ 2 = ( A x − b ) T ( A x − b ) = ( x T A T − b T ) ( A x − b ) = x T A T A x − b T A x − x T A T b + b T b = x T A T A x − ( A T b ) T x − ( A T b ) T x + b T b

$$\begin{array}{l} E(\mathbf{x})=\|\mathbf{A} \mathbf{x}-\mathbf{b}\|^{2} \\ =(\mathbf{A} \mathbf{x}-\mathbf{b})^{T}(\mathbf{A} \mathbf{x}-\mathbf{b}) \\ =\left(\mathbf{x}^{\mathrm{T}} \mathbf{A}^{\mathrm{T}}-\mathbf{b}^{\mathrm{T}}\right)(\mathbf{A} \mathbf{x}-\mathbf{b}) \\ =\mathbf{x}^{\mathrm{T}} \mathbf{A}^{\mathrm{T}} \mathbf{A} \mathbf{x}-\mathbf{b}^{\mathrm{T}} \mathbf{A} \mathbf{x}-\mathbf{x}^{\mathrm{T}} \mathbf{A}^{\mathrm{T}} \mathbf{b}+\mathbf{b}^{\mathrm{T}} \mathbf{b} \\ =\mathbf{x}^{\mathrm{T}} \mathbf{A}^{\mathrm{T}} \mathbf{A} \mathbf{x}-\left(\mathbf{A}^{\mathrm{T}} \mathbf{b}\right)^{\mathrm{T}} \mathbf{x}-\left(\mathbf{A}^{\mathrm{T}} \mathbf{b}\right)^{\mathrm{T}} \mathbf{x}+\mathbf{b}^{\mathrm{T}} \mathbf{b} \end{array}$$ E(x)=∥Ax−b∥2=(Ax−b)T(Ax−b)=(xTAT−bT)(Ax−b)=xTATAx−bTAx−xTATb+bTb=xTATAx−(ATb)Tx−(ATb)Tx+bTb​令$E(\mathbf{x})$对$x$一阶导数为0：
$$\begin{array}{c}\frac{\partial E}{\partial \mathbf{x}}=2{\mathbf{A}}^{\mathbf{T}}\mathbf{Ax}-2{\mathbf{A}}^{\mathrm{T}}\mathbf{b}=0\\ \mathbf{x}={\left({\mathbf{A}}^T\mathbf{A}\right)}^{-1}{\mathbf{A}}^T\mathbf{b}\end{array}$$这样，求解最小二乘问题只需求解一个线性方程组即可，不再需要像梯度下降那样迭代了。

矩阵求导可以参考：
https://zhuanlan.zhihu.com/p/273729929

对于特殊的$Ax=0$型超定方程组：
$$\frac{\partial E}{\partial \mathbf{x}}=2{\mathbf{A}}^{\mathbf{T}}\mathbf{Ax}=0$$通过以下方式获得最小二乘解：
[V,D] = eig(A'*A)
其中D是特征值对角矩阵(特征值沿主对角线降序)，V是对应D特征值的特征向量(列向量)组成的特征矩阵,A’表示A的转置。其最小二乘解为V(1)，即系数矩阵A最小特征值对应的特征向量就是超定方程组$Ax=0$的最小二乘解。

案例：求单应性矩阵的matlab代码【原文连接】

% 返回值 H 是一个3*3的矩阵
% pts1 和 pts2是2*4的坐标矩阵对应特征点的(x,y)坐标
n = size(pts1,2);
A = zeros(2*n,9);
A(1:2:2*n,1:2) = pts1';
A(1:2:2*n,3) = 1;
A(2:2:2*n,4:5) = pts1';
A(2:2:2*n,6) = 1;
x1 = pts1(1,:)';
y1 = pts1(2,:)';
x2 = pts2(1,:)';
y2 = pts2(2,:)';
A(1:2:2*n,7) = -x2.*x1;
A(2:2:2*n,7) = -y2.*x1;
A(1:2:2*n,8) = -x2.*y1;
A(2:2:2*n,8) = -y2.*y1;
A(1:2:2*n,9) = -x2;
A(2:2:2*n,9) = -y2;

[evec,~] = eig(A'*A);
H = reshape(evec(:,1),[3,3])';
H = H/H(end); % make H(3,3) = 1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22