SparseLM Demo示例分析

编译好SparseLM后，其中的第一个Demo是Chained Rosenbrock function，其方程式如下：

f1=10*(x[i]*x[i]-x[i+1])f2=x[i]-1

这个方程在JORGE的文章中也写成：

在Andrew的文章中又是：

而Jan Vlcek写的是：

我不是学数学的，不知道怎么还写了这么多种不同的样子的方程；此处暂且按照demo中的样子来分析吧。

首先还是看一下运行结果：

观察其代码，和levmar使用其实是非常类似；都从一个迭代主函数Sparse sparselm_dercrs()
调用另外两个函数，一个是描述测量值和参数直接的代数关系，一个是用来描述雅克比矩阵；唯一的区别是，雅克比矩阵在此处是一个稀疏矩阵，其存储方式和之前的稠密矩阵有极大的不同。

It accepts sparse Jacobians encoded in either compressed row storage (CRS) or compressed column storage (CCS) (a.k.a. Harwell-Boeing) format, allowing user applications to choose the representation that is most natural to them. 

从这段话可以看出，一共有行优先（CRS）和列优先（CCS）两种存储方式（ps.突然想起来Eigen的数据就是列优先格式进行存储的）。

现在开始一行一行的看demo的代码

/* Chained Rosenbrock function *///必须是偶数，这里如此写是因为程序中判定偶数的方式导致的#define   MCHROSEN   10 #define   NCHROSEN    2*(MCHROSEN-1)int main(){    register int i;    //循环迭代中一些参数，全部默认就好了    double opts[SPLM_OPTS_SZ], info[SPLM_INFO_SZ];    //待优化参数的数组    double p[MCHROSEN];    //m：参数的个数，n: 测量值的个数，ret：迭代次数    int m, n, ret;    //雅克比矩阵中非零项个数    int nnz;    //设置参数，默认    opts[0]=SPLM_INIT_MU;     opts[3]=SPLM_STOP_THRESH;    opts[4]=SPLM_DIFF_DELTA;    opts[5]=SPLM_CHOLMOD;    m=MCHROSEN;    //从宏中看到NCHROSEN=(MCHROSEN-1)*2其实没有特别的要求，只需求测量值的数目不能太少    n=NCHROSEN;    //这里刚开始理解起来可能很费解，我们这样想吧。    //一组观测值，含有f1,f2两个方程，分别对x[i],x[i+1]求偏导，那么就可以得到一个2×2的矩阵，其中还有一个零，因此剩下三个就是非零的值。    //因此这里地方大概意思就是：测量值对数*3    nnz=3*NCHROSEN/2;    //赋初值    for(i=0; i<MCHROSEN; i+=2)    {      p[i]=-1.2;      p[i+1]=1.0;    }     //调用迭代函数     ret=sparselm_dercrs(chainedRosenbrock,//描述关系的函数指针     chainedRosenbrock_anjacCRS, //生成CRS格式代数雅克比矩阵     p, //参数初值和结果保存指针     NULL, //测量值，若是零向量就设为NULL     m, //参数个数     0, //参数中不能调整的值的个数     n, //观测值个数     nnz,//雅克比矩阵中非零个数     -1, //J^t*J中非零个数，如果不知道就设为-1     1000, //最大迭代次数     opts, info, //配置参数     NULL//这个参数是用来向程序中传递额外的数据的，没有就设置为空); 

//初始化的时候当i为奇数，p[i]=-1.2；当i为偶数，p[i]=1.0;void chainedRosenbrock(double *p, double *hx, int m, int n, void *adata){  register int k, k1, i;  for(k=0; k<n; ++k)  {    //下边这个操作其目的是为了对应上观测值和参数，因为默认是第i对数据，对应p[i]和p[i+1]两个参数    k1=k+1; // k从0开始计数,转化到以1开始计数    i=DIV(k1+1, 2) - 1; // 再将i转化到以0开始计数    if(k1%2==1) // k1是奇数时      hx[k]=10.0*(p[i]*p[i]-p[i+1]);    else // k1是偶数时      hx[k]=p[i]-1.0;  }}//雅克比矩阵//这个应该是默认的一个函数void chainedRosenbrock_anjacCRS(double *p, struct splm_crsm *jac, int m, int n, void *adata){  register int k, k1, i;  int l;  for(k=l=0; k<n; ++k)  {    //这个地方的理解放到最后边，这个地方比较重要!!!!!    //大致理解就是记住第k行jac矩阵中有几个非0值    jac->rowptr[k]=l;    //同上    k1=k+1;    i=DIV(k1+1, 2) - 1;    if(k1%2==1)    {      //hx[k]=10*(p[i]*p[i]-p[i+1])      //赋值在（k,l)      jac->val[l]=20.0*p[i];       //列标移动并未当前的赋值      //保存列标签column index      jac->colidx[l++]=i;      jac->val[l]=-10.0;       jac->colidx[l++]=i+1;    }    else    {      // k1 even, hx[k]=p[i]-1.0      jac->val[l]=1.0;      jac->colidx[l++]=i;    }  }  //移动到雅克比矩阵末尾  jac->rowptr[n]=l;}

以上大概就是整个的一个代码，Demo中还提供了另外三个雅克比矩阵的函数，功能类似，变化在于存储方式（CRS、CSS)和稀疏矩阵赋值方式。
赋值方式还可以调用

extern int splm_stm_nonzeroval(struct splm_stm *sm, int i, int j, double val);

相对于来说更为直观。

关于其他细节，以后再学习。

关于上边提到的非常重要的细节，首先拿到struct splm_crsm的定义：

struct splm_crsm{    int nr, nc;   /* #rows, #cols 分别是稀疏矩阵行列值 */    int nnz;      /* 稀疏矩阵中非零值个数nnz */    double *val;  /* 指向用于保存稀疏矩阵中非零值的指针，其大小为nnz*/    int *colidx;  /* 保存val中对应值在稀疏矩阵中某一列信息的指针，其大小和val是一模一样的，为nnz*/    int *rowptr;  /* 这个数据长度为nr+1,记录从当前行到第一行总共有多少个非零值，并且rowptr[nr]=nnz */};

举个例子：

结构体中存储的数据应该是如下的：
nr = 5,nc = 5；
nnz = 9;
val={1 1 1 1 1 1 1 1 1}
colidx={1 2 3 1 2 3 2 1 2}
rowptr={0 3 4 6 7 9}