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}