Levenberg-Marquardt算法问题 ,请教关于levmar-2.5中dlevmar_dif()的使用

ml1306

这个是Levenberg-Marquardt算法,找了很多资料,都没有发现他要求提供的函数关系的原型是如何得出的
int dlevmar_dif(
      void (*func)(double *p, double *hx, int m, int n, void *adata),
      double *p, double *x, int m, int n, int itmax, double *opts,
      double *info, double *work, double *covar, void *adata);
/* Secant version of the LEVMAR_DER() function above: the Jacobian is approximated with
* the aid of finite differences (forward or central, see the comment for the opts argument)
*/
int LEVMAR_DIF(
  void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata), /* functional relation describing measurements. A p \in R^m yields a \hat{x} \in  R^n */
  LM_REAL *p,         /* I/O: initial parameter estimates. On output has the estimated solution */
  LM_REAL *x,         /* I: measurement vector. NULL implies a zero vector */
  int m,              /* I: parameter vector dimension (i.e. #unknowns) */
  int n,              /* I: measurement vector dimension */
  int itmax,          /* I: maximum number of iterations */
  LM_REAL opts[5],    /* I: opts[0-4] = minim. options [\mu, \epsilon1, \epsilon2, \epsilon3, \delta]. Respectively the
                       * scale factor for initial \mu, stopping thresholds for ||J^T e||_inf, ||Dp||_2 and ||e||_2 and
                       * the step used in difference approximation to the Jacobian. Set to NULL for defaults to be used.
                       * If \delta<0, the Jacobian is approximated with central differences which are more accurate
                       * (but slower!) compared to the forward differences employed by default.
                       */
  LM_REAL info[LM_INFO_SZ],
                                                   /* O: information regarding the minimization. Set to NULL if don't care
                      * info[0]= ||e||_2 at initial p.
                      * info[1-4]=[ ||e||_2, ||J^T e||_inf,  ||Dp||_2, mu/max[J^T J]_ii ], all computed at estimated p.
                      * info[5]= # iterations,
                      * info[6]=reason for terminating: 1 - stopped by small gradient J^T e
                      *                                 2 - stopped by small Dp
                      *                                 3 - stopped by itmax
                      *                                 4 - singular matrix. Restart from current p with increased mu
                      *                                 5 - no further error reduction is possible. Restart with increased mu
                      *                                 6 - stopped by small ||e||_2
                      *                                 7 - stopped by invalid (i.e. NaN or Inf) "func" values. This is a user error
                      * info[7]= # function evaluations
                      * info[8]= # Jacobian evaluations
                      * info[9]= # linear systems solved, i.e. # attempts for reducing error
                      */
  LM_REAL *work,     /* working memory at least LM_DIF_WORKSZ() reals large, allocated if NULL */
  LM_REAL *covar,    /* O: Covariance matrix corresponding to LS solution; mxm. Set to NULL if not needed. */
  void *adata)       /* pointer to possibly additional data, passed uninterpreted to func.
                      * Set to NULL if not needed
                      */
以上是函数说明,需要的是用这个函数拟合曲线y=c/(1+exp(a+bx))+d;(提供四个点)
以下是在Nemerical Recipes in C这本书中函数所需要的原型,不知道如何转化为dlevmar_dif所要求的形式
void funcs(float x, float a[], float *y, float dyda[], int na)
{   
    int i;
    float fac,ex,arg;
    *y=0.0;
    for(i=1;i<=na-1;i+=3)
    {
    fac=a+a[i+1]*x;
    ex=exp(fac);
    arg=(1+ex)*(1+ex);
    *y=a[i+2]/(1+ex)+a[4];
    dyda=(-a[i+2])*ex/arg;
    dyda[i+1]=(-a[i+2])*x*ex/arg;
    dyda[i+2]=1/(1+ex);
    dyda[i+3]=1;
    }
}
那位大侠知道如何使用或者有更好的曲线拟合库的请指教
它的网站 http://www.ics.forth.gr/~lourakis/levmar/
附上安装方法
levmar-2.5安装步骤，测试通过

转载自http://blog.csdn.net/lkbwx/archive/2010/02/19/5311802.aspx



点击http://www.netlib.org/clapack/CLAPACK-3.1.1-VisualStudio.zip 下载clapack

点击http://www.cmake.org/files/v2.8/cmake-2.8.0-win32-x86.exe   下载cmake

点击http://www.ics.forth.gr/~lourakis/levmar/levmar-2.5.tgz       下载LEVMAR

打开CLAPACK-3.1.1-VisualStudio文件夹，先把\LIB\Win32的lib都删了,以免混淆

双击clapack.vcproj打开工程项目文件，下面的各编译步骤都编译成debug模式

依次编译F2CLIBS,tmglib,blas,CLAPACK,结果是在\LIB\Win32下依次生成了libf2cd.lib,tmglibd.lib,BLASd.lib,clapackd.lib

使用cmake作用于levmar-2.5文件夹，得到工程项目文件LEVMAR.vcproj，打开它

工具->选项->项目和解决方案->vc++目录->包含文件处添加     \INCLUDE目录

工具->选项->项目和解决方案->vc++目录->库文件处添加        \LIB\Win32目录

项目属性->连接器->输入->附加依赖项libf2cd.lib BLASd.lib clapackd.lib tmglibd.lib

项目属性->连接器->输入->忽略特定库Libcmtd.lib

欲在LEVMAR中调用CLAPACK，还需要打开LEVMAR的misc_core.c文件，然后在文件起始处写上#include "blaswrap.h"

参考

http://www.cnblogs.com/random_walk/archive/2009/11/13/1602682.html

http://www.wilmott.com/messageview.cfm?catid=34&threadid=51193

ml1306

找到了例子中的原型，发现其实很简单，但是初值的选择还是有一定的难度，网上资料很少，1stOpt具有全局最优的能力，但是看不到，希望哪位高手能提供一个初值的粗略估算方法
以下是实现的几个简单曲线拟合实例，难点就在于初值选择
#include <stdio.h>
#include <process.h>

//#include <stdlib.h>
#include <math.h>
//#include <float.h>

#include "levmar.h"
#include "compiler.h"
#include "misc.h"
#include "test1.h"
//因为程序是C++，而CLAPACK是f2c程序转换的C语言版本，所以在此处用extern关键字调用
extern"C"
{
#include <f2c.h>
#include <clapack.h>
}
#define ROSD 105.0
#define EV 2.71828183
double ava[6]={1,2,3,4,5,6};
double bva[6]={1,2,5,7,40,43};
///* Rosenbrock function, global minimum at (1, 1) */
//void ros(double *p, double *x, int m, int n, void *data)
//{
//        register int i;
//
//        for(i=0; i<n; ++i)
//                x=((1.0-p[0])*(1.0-p[0]) + ROSD*(p[1]-p[0]*p[0])*(p[1]-p[0]*p[0]));
//}
//
//void jacros(double *p, double *jac, int m, int n, void *data)
//{
//        register int i, j;
//
//        for(i=j=0; i<n; ++i){
//                jac[j++]=(-2 + 2*p[0]-4*ROSD*(p[1]-p[0]*p[0])*p[0]);
//                jac[j++]=(2*ROSD*(p[1]-p[0]*p[0]));
//        }
//}
//
//#define SIZE 4
//
//void func(double *p, double *x, int m, int n, void *data)
//{
//        int i;
//        for(i = 0; i < n; ++i)
//                x = (p[0] - 3.1415926) * (p[0] - 3.1415926) +
//                (p[1] - 2.0)  * (p[1] - 2.0);
//}

/* Wood's function, minimum at (1, 1, 1, 1) */
//void wood(double *p, double *x, int m, int n, void *data)
//{
//        register int i;
//
//        for(i=0; i<n; i+=6)
//        {
//                x=10.0*(p[1] - p[0]*p[0]);
//                x[i+1]=1.0 - p[0];
//                x[i+2]=sqrt(90.0)*(p[3] - p[2]*p[2]);
//                x[i+3]=1.0 - p[2];
//                x[i+4]=sqrt(10.0)*(p[1]+p[3] - 2.0);
//                x[i+5]=(p[1] - p[3])/sqrt(10.0);
//        }
//}
//void func(double *p, double *x, int m, int n, void *data)
//{
//        int i;
//        double a[3]={1,2,3};
//        double b[3];
//        for(i=0;i<3;i++)
//        b=2*a*a*a+3*a*a+a;
//        for(i = 0; i < n; i+=3)
//        {
//                x = b-(p[0]*a*a*a+p[1]*a*a+p[2]*a);
//                x[i+1] =b[i+1]-(p[0]*a[i+1]*a[i+1]*a[i+1]+p[1]*a[i+1]*a[i+1]+p[2]*a[i+1]);
//                x[i+2] =b[i+2]-(p[0]*a[i+2]*a[i+2]*a[i+2]+p[1]*a[i+2]*a[i+2]+p[2]*a[i+2]);
//        }
//}
//void jacfunc(double *p, double *jac, int m, int n, void *data)
//{
//        register int i, j;
//
//        for(i=j=0; i<n; ++i){
//                jac[j++]=-1;
//                jac[j++]=-1;
//
//                jac[j++]=-4;
//                jac[j++]=-2;
//        }
//}
//void funcs(float x, float a[], float *y, float dyda[], int na)
//{   
//    int i;
//    float fac,ex,arg;
//    *y=0.0;
//    for(i=1;i<=na-1;i+=3)
//    {
//    fac=a+a[i+1]*x;
//    ex=exp(fac);
//    arg=(1+ex)*(1+ex);
//    *y=a[i+2]/(1+ex)+a[4];
//    dyda=(-a[i+2])*ex/arg;
//    dyda[i+1]=(-a[i+2])*x*ex/arg;
//    dyda[i+2]=1/(1+ex);
//    dyda[i+3]=1;
//    }
//}
void Weibull(double *p, double *x, int m, int n, void *data)
{
        int i;


        //for(i=0;i<4;i++)
        //bva=1-3*exp(-4*pow(ava,1));
        //for(i=0; i<m; ++i) printf("%.7g ", b);printf("\n");
        for(i = 0; i < n; i++)
        {
                double u = pow(ava,p[3]);
                x = bva-(p[0]-p[1]*exp(-p[2]*u));
        }
}
////void jacWeibull(double *p, double *jac, int m, int n, void *data)
////{
////        register int i, j;
////
////        for(i=j=0; i<n; ++i)
////        {
////                jac[j++]=-1;
////                jac[j++]=exp(-p[2]*pow(ava,p[3]));
////                jac[j++]=-p[1]*exp(-p[2]*pow(ava,p[3]))*-p[2]*pow(ava,p[3])*(log(ava)/log(EV));
////        }
////}

void Logistic(double *p, double *x, int m, int n, void *data)
{
        int i;
       

        //for(i=0;i<3;i++)
        //b=2/(1+1*exp(-2*a));
        //for(i=0; i<m; ++i) printf("%.7g ", b);printf("\n");
        for(i = 0; i < n; i++)
        {
                x = bva-p[0]/(1+p[1]*exp(-p[2]*ava));
        }
}
void funa(double *p, double *x, int m, int n, void *data)
{
        int i;


       
        //for(i=0; i<m; ++i) printf("%.7g ", b);printf("\n");
        for(i = 0; i < n; i++)
        {
                x = bva-p[0]/(1+exp(-p[1]-p[2]*ava));
        }
}
double corrcoef(double *d1, double *d2, int dataL)
{
        int i;
        double xy=0, x=0, y=0, xsum=0, ysum=0;
        double corrc;
        for (i=0; i<dataL; i++)
        {
                xsum += d1;
                ysum += d2;
        }
        xsum /= dataL;
        ysum /= dataL;
        for (i=0; i<dataL; i++)
        {
                d1 -= xsum;
                d2 -= ysum;
                x += d1 * d1;
                y += d2 * d2;
                xy += d1 * d2;
        }
        corrc = abs(xy) / (sqrt(x) * sqrt(y));
        return corrc;
}

int main()
{
        register int i;
        int  ret;
        double p[5], // 5 is max(2, 3, 5)
                x[16]; // 16 is max(2, 3, 5, 6, 16)
        int m, n;
        double opts[LM_OPTS_SZ], info[LM_INFO_SZ];
        char *probname[]={
                "Rosenbrock function",
                "modified Rosenbrock problem",
                "Powell's function",
                "Wood's function",
                "Meyer's (reformulated) problem",
                "Osborne's problem",
                "helical valley function",
                "Boggs & Tolle's problem #3",
                "Hock - Schittkowski problem #28",
                "Hock - Schittkowski problem #48",
                "Hock - Schittkowski problem #51",
                "Hock - Schittkowski problem #01",
                "Hock - Schittkowski modified problem #21",
                "hatfldb problem",
                "hatfldc problem",
                "equilibrium combustion problem",
                "Hock - Schittkowski modified #1 problem #52",
                "Schittkowski modified problem #235",
                "Boggs & Tolle modified problem #7",
                "Hock - Schittkowski modified #2 problem #52",
                "Hock - Schittkowski modified problem #76",
        };

        opts[0]=LM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-20;
        opts[4]= LM_DIFF_DELTA;
//        char *reason[]={
//                "",
//                "stopped by small gradient J^T e",
//                "stopped by small Dp",
//                "stopped by itmax",
//                "singular matrix. Restart from current p with increased mu",
//                "no further error reduction is possible. Restart with increased mu",
//                "stopped by small ||e||_2"
//        };
//
//        int i, m, n, ret;
//        n = 2; m = 2;
//        double p[10], x[10];
//        double info[100];
//
//        char JOBU;
//        char JOBVT;
//
////        int i;
//
//        //数据类型integer是fortran里的。这里在C++下可以使用的原因是f2c.h文件中已经作了定义
//        integer M = SIZE;
//        integer N = SIZE;
//        integer LDA = M;
//        integer LDU = M;
//        integer LDVT = N;
//        integer LWORK;
//        integer INFO;
//
//        integer mn = min( M, N );
//
//        integer MN = max( M, N );
//
//        double a[SIZE*SIZE] = { 16.0, 5.0, 9.0 , 4.0, 2.0, 11.0, 7.0 , 14.0, 3.0, 10.0, 6.0, 15.0, 13.0, 8.0, 12.0, 1.0};
//        double s[SIZE];
//        double wk[201];
//        double uu[SIZE*SIZE];
//        double vt[SIZE*SIZE];
//
//        JOBU = 'A';
//
//        JOBVT = 'A';
//
//        LWORK = 201;

        /* Subroutine int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
        doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *
        ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
        integer *info)
        */
        //dgesvd_( &JOBU, &JOBVT, &M, &N, a, &LDA, s, uu, &LDU, vt, &LDVT, wk, &LWORK, &INFO);

        //printf("INFO=%d \n", INFO );         

        //for ( i= 0; i< SIZE; i++ ) {
        //        printf("s[ %d ] = %f\n", i, s[ i ] );
        //}

        //

        //p[0] = 2; p[1] = 3;  // initial values
        //for(i = 0; i < n; i++) x = 0.0;
        //ret=dlevmar_dif(func, p, x, m, n, 1000, NULL, info, NULL, NULL, NULL); // with numerical jacobian
        //printf("\nUnconstrained, with numerical jacobian\n--------------------------------------");
        //printf("\nLevenberg-Marquardt returned in %d iteration(s): %s\nSolution: ", ret, reason[int(info[6])]);
        //for(i=0; i<m; ++i) printf("%.7g ", p);
        //printf("\n");

        ///* Woods's function */
        //m=4; n=6;
        //p[0]=-3.0; p[1]=-1.0; p[2]=-3.0; p[3]=-1.0;
        //for(i=0; i<n; i++) x=0.0;
        ////x[1]=1.0;

        //ret=dlevmar_dif(wood, p, x, m, n, 1000, NULL, info, NULL, NULL, NULL);  // no jacobian
        //for(i=0; i<m; ++i) printf("%.7g ", p);
        //printf("\n");

        //m=2; n=2;
        //p[0]=2.0; p[1]=1.0;
        //for(i=0; i<n; i++) x=0.0;
        //ret=dlevmar_dif(ros, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
        ////ret=dlevmar_dif(ros, p, x, m, n, 1000, opts, info, NULL, NULL, NULL);  // no Jacobian
        //for(i=0; i<m; ++i) printf("%.7g ", p);
        //printf("\n");
        //double ava[3]={1,2,3};
        //double bva[3]={3,6,10};

        //for(i=0;i<3;i++)
        //b=2/(1+1*exp(-2*a));
        //for(i=0; i<m; ++i) printf("%.7g ", b);printf("\n");
        m=4; n=4;
        p[0]=23;p[1]=23;p[2]=0.0008;p[3]=2;
        for(i=0; i<n; i++) x=0.0;
        //for(i=0;i<6;i++)
        //        bva=5/(1+exp(-1-1*ava));
        //ret=dlevmar_dif(funa,p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
        //ret=dlevmar_dif(Logistic,p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
        ret=dlevmar_dif(Weibull,p, x, m, n, 10000, opts, info, NULL, NULL, NULL);

       
        //ret=dlevmar_dif(ros, p, x, m, n, 1000, opts, info, NULL, NULL, NULL);  // no Jacobian
        for(i=0; i<m; ++i) printf("%.7g ", p);
        printf("\n");
        for(i=0; i<n; ++i) printf("%.7g ", bva);printf("\n");
        double cva[6];
        for(i=0;i<n;i++)
//        cva=p[0]/(1+exp(-p[1]-p[2]*ava));
        //cva=p[0]/(1+p[1]*exp(-p[2]*ava));  //Logistic
        cva=p[0]-p[1]*exp(-p[2]*pow(ava,p[3]));
        for(i=0; i<n; ++i) printf("%.7g ", cva);printf("\n");

        double sva=0;

        sva = corrcoef(cva,bva,n);
        printf("%.7g ", sva);
       
        system("pause");
        return 0;
}