最小二乘计算最优解不管是哪个行业肯定都用到的非常多。对于遥感图像处理中，尤其是对图像进行校正处理，关于控制点的几种校正模型中，都用到最小二乘来计算模型的系数。比如几何多项式，或者通过GCP求解RPC系数，以及RPC的间接优化等都是离不开最小二乘的。下面使用MTL库编写的最小二乘求解AX=B形式的X最优解。
关于MTL库的类型定义可以参考之前写的求解特征值和特征向量那篇博客。地址为：http://blog.****.net/liminlu0314/article/details/8957155。
首先是函数定义：
/**
* @brief 求解矩阵Ax=B的最小二乘解集
* 注意：矩阵A的行数与B的行数必须相等，否则不能计算，此外X的大小与A的列数必须一致
* @param A 系数矩阵
* @param B 增广矩阵中的B
* @param X 解集
* @return 是否计算成功
*/
bool SloveMartix(const Matrix A, const Matrix B, Vector &X);
接下来就是函数实现，其实用MTL库很简单，一共也就十几行的代码。实现代码如下。
bool SloveMartix(const Matrix A, const Matrix B, Vector &X)
{
/*/
/* 求解矩阵 AX = B 的最小二乘解集
/* 注意：矩阵A的行数与B的行数必须相等，否则不能计算，此外X的大小与A的列数必须一致
/* 按最小二乘法组成法方程：
/* (A'*A)*X = (A'*B) 其中A'为A的转置矩阵
/* 设：N=(A'*A), W=(A'*B) 法方程为 N*X=W
/*/ int irows = A.nrows();
int icols = A.ncols(); // 矩阵A的行数与B的行数必须相等，否则不能计算，此外X的大小与A的列数必须一致
if (irows != B.nrows() || icols != X.size())
return false; Matrix At(icols, irows); //A的转置矩阵
transpose(A,At); Matrix N(icols, icols);
mult(At, A, N); //计算N=(A'*A) Matrix W(icols, 1);
mult(At, B, W); //计算W=(A'*B) Vector pvector(icols);
lu_factor(N, pvector); Vector b(icols);
for(size_t i=0; i<b.size(); i++)
b[i] = W(i,0); lu_solve(N, pvector, b, X); return true;
}
上面的代码与Matlab比较过，计算的结果是完全一样，除了小数点十几位后有些精度上的问题。不过一般肯定用不到那么高的精度。
参考资料：
[1]http://zh.wikipedia.org/wiki/%E6%9C%80%E5%B0%8F%E4%BA%8C%E4%B9%98%E6%B3%95
[2] http://baike.baidu.com/view/139822.htm
[3]http://blog.sciencenet.cn/blog-430956-621997.html