使用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

使用MTL库求解最小二乘解相关推荐

使用MTL库求解矩阵特征值和特征向量
关于矩阵的特征值和特征向量求解,大部分的数学运算库都进行了提供,下面是使用MTL库的接口进行封装. #include <mtl/matrix.h> #include <mtl/mtl ...
关于MFC下使用MTL库编译错误的问题
在使用Matrix Template Library(MTL)库进行矩阵运算还是很给力的,但是遇到了一个比较悲剧的问题就是,在控制台程序中一切完美,而在MFC下使用该库会编译不过去,(编译环境VS20 ...
mtl库在GCC编译器下的使用
最近一直在改造算法库,将其移植到Linux平台下.使用GCC编译器是发现MTL库中出现一大堆的问题.使用的MTL库下载地址为:http://osl.iu.edu/research/mtl/downlo ...
给MTL库添加求行列式值
在使用MTL库的时候,发现mtl库没有求行列式的值的函数,google了一把,找到下面的网页 How to Find Determinant of nxn matrix? 参考里面的说明,给mtl加上 ...
利用sympy库求解常微分方程：dsolve()函数
[小白从小学Python.C.Java] [计算机等级考试+500强双证书] [Python-数据分析] [sympy库的使用] 使用sympy库求解常微分方程 .dsolve()方法选择题下列说 ...
基于 python pulp 库求解船舶泊位调度线性规划问题
目录基于 python pulp 库求解船舶泊位调度线性规划问题泊位调度问题建模代码实现准备包代码讲解绘制图像完整代码题外话基于 python pulp 库求解船舶泊位调度线性规划问 ...
fortran使用MKL函数库求解普通稀疏矩阵与向量的乘积
使用MKL函数库中的mkl_dcsrgemv可以实现普通稀疏矩阵与向量的乘积 program mainimplicit noneinteger :: n, n1, ja(9), ia(7)real*8 ...
C++给蜗蜗留口气代码源oj题库求解
蜗蜗最近爱上了下蜗蜗棋,蜗蜗棋和地球上的围棋有相似之处.他最近在分析自己以前下过的一些残局(即将结束的棋局).现在有一个 n×n的棋盘.每个位置有一个数,这个数为 0 ,则表示这里没有棋子,是空的:这 ...
使用OpenCV库快速求解相机内参
本文主要介绍如何使用OpenCV库函数求解相机内参.具体可查阅官网:https://docs.opencv.org/master/dc/dbb/tutorial_py_calibration.html ...

使用MTL库求解最小二乘解

使用MTL库求解最小二乘解相关推荐

最新文章

热门文章