使用MTL库求解最小二乘解

最小二乘计算最优解不管是哪个行业肯定都用到的非常多。对于遥感图像处理中，尤其是对图像进行校正处理，关于控制点的几种校正模型中，都用到最小二乘来计算模型的系数。比如几何多项式，或者通过GCP求解RPC系数，以及RPC的间接优化等都是离不开最小二乘的。下面使用MTL库编写的最小二乘求解AX=B形式的X最优解。

关于MTL库的类型定义可以参考之前写的求解特征值和特征向量那篇博客。地址为：http://blog.csdn.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

[4]http://www.orbitals.com/self/least/least.htm

[5]http://blog.csdn.net/liminlu0314/article/details/8957155