最小二乘法拟合一些基本空间元素的C#实现 —— （一）术语、基本算法与思路

本系列的主要目的是将拟合算法用C#实现，相关的操作用到了Math.NET数值运算库，所以将一些基本概念，算法和对应的C#代码放到前面。

对于线性模型，例如直线和平面，算法大约分三步：

1. 参数化：将模型用参数化方式表达出来。
2. 构造Design（非线性中的Jacobian）矩阵：根据参数公式求所有参数的偏导数作为列(Column)，将实测值作为行(Row)。将偏差值作为等式右值。
3. 求解。

对于非线性模型，例如圆柱，Gauss-Newton算法类似：

1. 参数化：
2. 求（或猜）初始解。
3. 找到Jacobian矩阵及右值
4. 求解
5. 通过4计算逼近值，重复3和4
6. 迭代结束

线性模型的解法是基础，详细的分析请看这篇文章，这里只解释代码部分。在Math.Net中

对于矩阵\(A=\begin{bmatrix}A_{11} & A_{12}\\A_{21} & A_{22}\end{bmatrix}\)
A.Transpose() 代表A的转置矩阵\(A^T=\begin{bmatrix}A_{11} & A_{21}\\A_{12} & A_{22}\end{bmatrix}\)

这在坐标运算中是非常重要的，因为坐标运算中坐标系都表示为正交矩阵（抱歉，没找到证明），当计算一个点（用二维举例）\(P=\left(x,y\right)\)在新坐标系A中的坐标\(P^{‘}=\left(x^{‘},y^{‘}\right)\)时。\[P^{‘}=PA\]

而已知坐标系A中的坐标点\(P^{‘}\)求原坐标值可用\[P=P^{‘}A^T\]

一些类似的代码为：

A.Inverse()代表A的逆矩阵\(A^{-1}\)

A.Multiply(B)代表矩阵乘法\(A*B\)或者矩阵乘以向量\(A*b\)

对于常用的算法：

线性模型

已知Jacobian矩阵\(J\)及等式右值（偏差值向量）\(y\)，另\[Ax=y\]求解向量x

1. 正规方程（normal equations）
   \[A^TAx=A^Ty\]
   逆矩阵解法
   \[x=\left(A^TA\right)^{-1}A^Ty\]
   实现代码 

   x = A.Transpose().Multiply(A).Inverse().Multiply(A.Transpose().Multiply(y));

   更快和更稳定的方法将是对左侧的矩阵进行乔姆斯基分解（Cholesky decomposition），实现代码

   x = A.Transpose().Multiply(A).Cholesky().Solve(A.Transpose().Multiply(y));

   这两个在一些限定条件下是等价的，而这些限定条件在点数远大于参数的情况下是满足的。

2. QR分解
   推论见原文，代码更简单： 

   x=A.QR().Solve(y);

3. 奇异值分解（Singular Value Decomposition, SVD）
   A是方阵的情况下（即点数等于参数）代码为 

   x = A.Svd(true).Solve(y);

   A不是方阵的情况下代码为

// compute the SVD 
Svd svd = A.Svd(true);
// get matrix of left singular vectors with first n columns of U 
Matrix U1 = svd.U().SubMatrix(0, m,0, n); 
// get matrix of singular values 
Matrix S = new DiagonalMatrix(n, n, svd.S().ToArray()); 
// get matrix of right singular vectors 
Matrix V = svd.VT().Transpose(); 
x = V.Multiply(S.Inverse()).Multiply(U1.Transpose().Multiply(y));

非线性模型

这部分对于代码并没有什么区别，只是需要迭代（循环）及判断终止条件。