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
=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}