赖敬文

Msn: 252523163

Mail: jingwenlai@163.com

Website: http://blog.csdn.net/jingwenlai_scut

Description:

Posts by 赖敬文:

Matlab求解nonlinear least squares problems及相互验证
Jingwenlai contact: jingwenlai@163.com
1 简介 1
2 用实例来了解Matlab求解非线性最小二乘问题的步骤 1
2.1 第一个例子: 1
2.2 例子２: 3
2.3 可相互印证的另一种方法:利用Optimization Toolbox GUI: 5
3. 如何将Matlab求解的思路转换成C++代码 (Gauss-newton的C++实现) 6
3.1 说明 6
3.2 程序源代码 6
3.3 运行结果 10
3.4 说明 10
4. 总结 11

1 简介
Matlab 作为成熟的软件，其数值求解部分可以作为我们在编程中对结果进行检验的参照。本文通过两个例子，介绍如何利用Matlab来求解非线性最小二乘问题。以及如何将其转化成C++代码，最后给出部分讨论供参考。
2 用实例来了解Matlab求解非线性最小二乘问题的步骤
2.1 第一个例子:
F(x1,x2,x3)=(f1(x1,x2,x3),f2(x1,x2,x3),f3(x1,x2,x3))^t;
其中,
F1(x1,x2,x3)=3*x1-cos(x2*x3)-1/2;
F2(x1,x2,x3)=x1^2-81*(x2+0.1)^2+sinx3+1.06;
F3(x1,x2,x3)=exp(-x1*x2)+20*x3+(10*pi()-3)/3;
先编辑一个M-file用于定义F函数
function Z = Example1F(X) %save as file Example1F.m
x1 = X(1); x2 = X(2); x3 = X(3);
Z = zeros(1,3);
Z(1) = 3*x1-cos(x2*x3)-1/2;
Z(2) = x1*x1 – 81*(x2+0.1)^2 + sin(x3) + 1.06;
Z(3) = exp(-x1*x2) + 20*x3 + (10*pi()-3)/3;
end
再编辑其相应的jacobian函数（即在对于F中的每个f,分别对x1,x2,x3求导）。
function W = JF(X) %saved as file Example1JF.m
x1 = X(1); x2 = X(2); x3 = X(3);
W = [3 x3*sin(x2*x3) x2*sin(x2*x3);
2*x1 -162*(x2+0.1) cos(x3);
-x2*exp(-x1*x2) -x1*exp(-x1*x2) 20];
end
然后编辑这个是用于GaussNewton法求解的Matlab函数.
function [P,iter,err]=newdim(F,JF,P,delta,epsilon,max1)
%save as file netdim.m
%Input — F is the system saved as the M-file F.m
% – JF is the Jacobian of F saved as the M-file JF.m
% – P is the initial approximation to the solution
% – delta is the tolerance for P
% – epsilon is the tolerance for F(P)
% – max1 is the maximum number of iterations
%Output – P is the approximation to the solution
% – iter is the number of iterations required
% – err is the error estimate for P
Y = feval(F,P); %用P作为输入，F作为函数求值，根据前述，这里得出的值即为F的值
for k=1:max1
J = feval(JF,P); %计算F的Jacobian矩阵
Q = P-(J\Y’)'; %计算新的估值
Z=feval(F,Q);
err=norm(Q-P);
relerr=err/(norm(Q)+eps);
P=Q;
Y=Z;
iter = k;
if(err> [P,iter,err]=newdim(@Example1F,@Example1JF,X,0.001,0.01,100)
P = 0.5000 0.0000 -0.5236
iter = 3
err = 0.0016
这个结果与<<数值分析(第七版 影印版)Numerical Analysis>>P615的结果是一致的。

2.2 例子２:
该例子参见<> page 177
F(x,y) = (f1(x,y),f2(x,y))
Where,
f1(x,y)=x^2-2*x-y+0.5;
f2(x,y)=x^2+4*y^2-4;
计算其相应的Jacobian矩阵
J(x,y) =[ 2*x-1 -1 ;
2*x 8*y ];
其F函数对应的m-file如下：
function Z = F(X) %save as file F
x = X(1);y=X(2);
Z = zeros(1,2);
Z(1)=x^2-2*x-y+0.5;
Z(2)=x^2+4*y^2-4;
其相应的Jacobian函数对应的M-file如下：
function W=JF(X) % save as file JF
x=X(1);y=X(2);
W=[2*x-2 -1;2*x 8*y];
同样地，我们也编辑一个newdim.m函数，这个是从＜＜数值方法(第四版)Numerical Methods using Matlab,Fourth Edition>>的page 179中对非线性系统求解的Newton-RaphsonMethod中得来。
function [P,iter,err]=newdim(F,JF,P,delta,epsilon,max1)
% save as file newdim.m
%Input — F is the system saved as the M-file F.m
% – JF is the Jacobian of F saved as the M-file JF.m
% – P is the initial approximation to the solution
% – delta is the tolerance for P
% – epsilon is the tolerance for F(P)
% – max1 is the maximum number of iterations
%Output – P is the approximation to the solution
% – iter is the number of iterations required
% – err is the error estimate for P
Y = feval(F,P);
for k=1:max1
J = feval(JF,P);
Q = P-(J\Y’)';
Z=feval(F,Q);
err=norm(Q-P);
relerr=err/(norm(Q)+eps);
P=Q;
Y=Z;
iter = k;
if(err break
end
end
由此，则可以在Matlab的控制台中调入如下命令，对其进行求解。
X=[2.00 0.25]
[P,iter,err]=newdim(@F,@JF,X,0.001,0.001,100);
由此得到的答案与P177中例子所求得的答案一致。

如果想一步步的查看结果，则可以在函数的入口中加入一个kmax,用于控制其循环次数。
还可以对上面的newdim.m进行修改，改成更易理解的GaussNewton.m.
function [P,iter,err]=GaussNewton(F,JF,P,delta,epsilon,max1)
% save as file GaussNewton.m
%Input – F is the system saved as the M-file F.m
% – JF is the Jacobian of F saved as the M-file JF.m
% – P is the initial approximation to the solution
% – delta is the tolerance for P
% – epsilon is the tolerance for F(P)
% – max1 is the maximum number of iterations
%Output – P is the approximation to the solution
% – iter is the number of iterations required
% – err is the error estimate for P
for k=1:max1
Fx = feval(F,P);% use P and function F to evaluate the function’s values
Jx = feval(JF,P);% evaluate Jacobian
Q = P – (Jx\Fx’)'; % x = x+y
err=norm(Q-P);

P = Q; % for next iteration
iter = k;
if(err < delta)|(abs(Fx) 1
W = [3 x3*sin(x2*x3) x2*sin(x2*x3);
2*x1 -162*(x2+0.1) cos(x3);
-x2*exp(-x1*x2) -x1*exp(-x1*x2) 20];
end
end

再用optimtool打开Optimization Tool,设定其一些参数，lsqnonlin,Gauss-Newton,
其中最重要的是Objective function,而其写法参照上面。然后，若程序中含有对Jacobian的求解，则选择由程序给定，否则，则Solver会估计，然后选择初始值，则可以点Start求解，最后求解的结果将显示，最终的点也会显示，且可以看到进行了多少次迭代。这里面使用的方法其实跟前面的方法是一致的，只是这里用了内部的lsqnonlin函数，而前面我们通过自己构造Gauss-Newton方法来求解。

从中也可以看出其与第一个例子中的结果是一致的。(左下角为求解结果).
3. 如何将Matlab求解的思路转换成C++代码 (Gauss-newton的C++实现)
3.1 说明
这里使用的例子是<>p177的例子。其实现参照了P179的matlab代码。（TODO：现在还少一步对结果的判断）.
3.2 程序源代码
#include
#include

#include
using namespace std;

typedef CGAL::Taucs_solver_traits TaucsSolver;
typedef TaucsSolver::Matrix Matrix;
typedef TaucsSolver::Vector Vector;

//testing, n = 2, m = 2;
void func(float * p, float * hx,int m,int n)
{
float x = p[0]; float y = p[1];
hx[0] = x*x-2*x-y+0.5;
hx[1] = x*x+4*y*y-4;
}
//testing, n = 2, m = 2; sizeof(j) is nxm;
void objfunc(float * p ,float * j, int m,int n)
{
float x = p[0];
float y = p[1];
j[0] = 2*x-2; j[1] = -1;
j[2] = 2*x; j[3] = 8*y;
}

bool GaussNewton(void (*func)(float * p,float * hx, int m,int n),
void (*jacf)(float * p, float * j, int m,int n),
float * p,
float * x,
int m,
int n,
int itmax,
// float opts[4],
float eps3)
{
int k = 0; //iterater
bool found = false;
float * hx;
float * jac;

hx = new float [2];
jac = new float [4];

//compute the J and -f
(*func)(p,hx,m,n); //compute the f
cout<<"hx is "< cout<

//check whether the f is small enough
float p_eL2 = 0.0;
float tmp;
for(int i = 0; i < n; i++)
{
tmp = x[i]-hx[i];
p_eL2 += tmp*tmp;
}

TaucsSolver solver;
float d =1;
while( !found &&(k {
if(p_eL2 < eps3)
{
found = true;
std::cout<<"p_eL2 < eps3 is the stop critera"< }
//compute jac
cout<<"P before compute jac is "< cout<

(*jacf)(p,jac,m,n);

cout<<"the jac is "< cout< cout<

Matrix A(n,m);// 180x84 as Jacobian
Vector B(m); // 84 as f
Vector X(m);//m=84
//setting Jacobian Matrix
for(int i = 0; i < n; i++)
{
for(int j = 0; j < m; j++)
{
A.set_coef(i,j,jac[i*m+j]);
}
}
cout<<"The A is "< for(int i = 0; i < n; i++)
{
for(int j = 0; j < m; j++)
{
cout<<" "< }
cout< }
//setting -f for solving JX = -f;
for(int i = 0; i < m; i++)
{
B[i] = -hx[i];
}
cout<<"The -f is "< for(int i = 0; i < m; i++)
cout<<" "<

//solve
solver.linear_solver(A,B,X,d);

//for the next iteration,evaluate function at Pdp
for(int i = 0; i < m; i++)
{
p[i] += X[i];
}
cout< cout<

//compute the Y for next iteration
(*func)(p,hx,m,n); //compute the f

k++; //for next iteration
}

return true;
}

int main()
{
float * p;
float * x;
float * j;
int n = 2;
int m = 2;
j = new float [n*m];
p = new float [m];
x = new float [m];

p[0] = 2.0; p[1] = 0.25;
int maxIter = 5;

bool ret = GaussNewton(func,objfunc,p,x,m,n,maxIter,0.001);
if(ret = true)
{
cout<<"The result is:"< cout< }
return 0;
}

3.3 运行结果

可以看到，此结果与前面的例子2的求解结果是一致的。
3.4 说明
上面的用C++实现的版本是用CGAL内置的封装后的TaucsSolver_Trait来求解的，若需要更多的自由度，可以使用Taucs、CHOLMOD、UMFPACK或MKL等库提供的API来做，因为其暴露了更多的API，可以提供更细粒度的操作。
4. 总结
从上面的求解过程我们可以看到，利用gauss-newton方法来求解最主要的是写出合适的J以及JF,有了这两个函数之后，在C++求解的时候，实际上传给第三方数值求解库的都是这两个函数的function pointer. 另外，Matlab在初期的算法的正确性的难证上帮助很大，还有另一种方法是直接在C++中嵌入Matlab的求解引擎，然后在c++中调用matlab的m-file即可，这种做法有点类似于做脚本的做法。
另外，有一个小经验，做编辑的时候，Cube是一个最简单的图形，便于手工验证算法的正确性。

演示视频地址: http://v.youku.com/v_show/id_XMTY4MjY1Nzk2.html

paper地址 :http://gl.ict.usc.edu/Research/3DDisplay/

如何开始研究？－－ 图形学相关的库介绍
jingwenlai 2010-4-29

以下列出一些我在学习过程中碰到过的库及其相应的介绍，希望对刚入门的师弟师妹有一个大体的印象。 如果你觉得我有些说得不准确的地方，请及时帮忙指正，谢谢。

图形学大体上可以分为两类：一类是建模，一类是渲染

* 建模
现在我们所做的数字几何处理(Digital Signal Processing)大体上都是在做建模这一类，不管是做模型分析还是模型编辑。对于模型，我们首先需要了解它的数据结构。
** 数据结构：
我们所处理的模型大部分情况下是大家通常所说的mesh，实际上其是由点、线、面所组成的一个物体（如果大家有用过maya/3dsmax之类的软件的话，在framewire渲染模式下，可以看得比较清楚）。如何对这些元素进行组织，就得需要底层的数据结构的支持，而这些数据结构有很多种，现在大家最常用的是一个叫做“半边结构”（Halfedge）的数据结构，在这里有关于它的最初始的介绍http://www.holmes3d.net/graphics/dcel/， 在CGAL库中实现的Mesh结构也是半边结构，大家可以看它的Manual中的Polyhedral那一章。OpenMesh中实现的Mesh数据结构同样也是半边结构。 CGAL或OpenMesh对于Mesh部分提供的功能都差不多，所不同的只是其API的使用方式。 还有一些比较小众的库，如trimesh2

有了数据结构，就可以在上面做很多分析或者操作，比如简化，细分，分析模型的特征等。在这里，我们需要对数据处理的结构或模型用可视化的方式直观地显示出来，所以，也就引出了我们需要使用的图形API： OpenGL/DirectX
**图形API: OpenGL/DirectX
我们需要图形API的原因是想把我们对模型的处理结构及模型用可视化的方式直观显示出来，这就需要我们学习OpenGL或DirectX的使用，现在我们组使用的大部分是OpenGL,好处是可以跨平台，而且学习资源丰富。当然，现在有一些库对两者进行封装，提供统一的接口以供使用，其中最出名的算是OGRE. 但是我们所需要的大部分是对点线面的原子操作，大高层的库反倒不好用，所以大部分情况下还是用底层的API

好了，有了API，但是API只是接口，我们需要一些东西来承载它，因此也就需要界面库，因为整个OpenGL/DirectX其需要Context(上下文)来对它进行初始化。
** GUI
这里的界面库包括以下几类：一类是大部分介绍OpenGL编程的时候都会介绍的GLUT,一类是更丰富一点的MFC, Qt之类。CGAL官方出的subdivision tutorial中使用的界面是MFC, 熊老师在它的基础上扩展了很多，已经加入了很多他的算法。另外，Kobbelt那一组的Mario等人在SIGGRAPH 2007/2008, Eurographics 2007都出过course notes,包含十多章的PDF（基本上每一章是一个算法）以及一个源代码的打包文件，它里面使用的是GLUT，灿江当时做的时候是从那个上面开始扩展的。而我挑的是一个比较偏的叫FLTK的GUI库，上手会比较快，但当程序大了之后，比较难维护，现在开始往Qt那边转了。
我个人觉得比较好的学习方法是挑一个自己喜欢的不算太大的程序（如CGAL的subdivision tutorial或OpenMesh的SIGGRAPH 2007/2008 course notes中附带的程序）， 先把它整个程序流程弄懂，然后再在那个基础上慢慢扩展自己的算法。

熟悉了上面这三种类型的库之后，基本上就可以开始做事了，你可以针对你自己所想要实现的算法查找相应的资料，这里我主要给出我接触过的有关数值求解那部分的库。
** 数值求解
需要数值求解的原因大部分是因为近年来有关模型编辑（mesh editing）方面的内容，都会涉及到大的稀疏矩阵的求解，它底层的主要是一些向量与向量，向量与矩阵，矩阵与矩阵之间的操作
＊＊＊CPU方面：
最早的这方面的库是BLAS,LINPACK这一类库。这一类型的库比较底层，所以后面有许多人给出类似的实现，但是更高层一些（即面向程序员更易用）：
1.Intel 的MKL: 这个库在浙大有人用（08年大连开会认识的浙大的赵勇用得是这个）， AMD也有相应的实现叫acml，但是好像使用的人不多。
2.另一阵营则是学术界自行实现的了，其中一个是现在在论文中引用数较多的taucs, 另一个是http://www.cise.ufl.edu/research/sparse/ 介绍的其它库，包括UMFPACK, CHOLMOD
其性能的对比可参见《Experiences of Sparse Direct Symmetric Solvers》，《A Fast Multigrid Algorithm for Mesh Deformation》，这类型的数值求解库比较多，但是只需要选一个就可以了。
3.还有一种方法我值得可以推荐的是借用matlab, matlab可以用于难证上述的库求解的正确性，或者如李琳那样使用matlab提供的api, 直接把需要求解的过程交给matlab,然后把matlab求解出来的结果再在自己的程序中使用。熊老师在这方面做得工作更多，有一次听他说他后续的有关数值求解（如特征值特征向量等）的求解全部交由matlab处理，这种方式的优点是使原型开发更快，而且专注于更关心的内容，缺点可能是性能方面。这种方式我个人现在比较推荐，可以等算法实现稳定后，再考虑性能优化的问题。

***GPU方面:
现在随着CUDA的流行，CUDA中提供了cuBlas,提供一些基本的与BLAS相对应的操作，现在后续地基于它的数值求解算法也有人在做，速度会比CPU的快。
这部分目前了解不多

* 渲染
渲染这一部分，我看得不多，一般来说我觉得分为两派，一派是追求尽可能的真实还原真实世界的场景，这一部分包含的算法包括Raytracing, radiosity等，另一部分则追求non-photorealistic（非真实）, 其中有一些是希望利用算法判断模型特征然后勾勒出一些轮廓，有一些是希望做出漫画效果，上一次灿江发的2010做视频实时渲染成水彩画的我想也属于这一类。

学习的过程中走过一些弯路，希望本文对你有所帮助。

1。选择研究方向要看学术影响力和工业上的应用价值。

2。做SIGGRAPH的三个标准：

Motivation：要有好的“故事”，大部分人同意从客观上值得做。

新的（好的）视觉特效：尤其是视觉艺术家能认可。

高性能：起码在2个指数级。

3。做SIGGRAPH的6招：

第一招：补洞。 Read more