Simplex, Non Linear Least Squares

 

References

  1. Wolfram, Mathematica
  2. Press, et al., Numerical Recipes, 2nd Edition
  3. Billo, Excel for Chemists, Ch. 12
  4. Bevington, Data Reduction and Error Analysis for the Physical Sciences, Ch. 11
  5. DeVries, A First Course in Computational Physics, p. 138
  6. Diamond and Hanratty, Spreadsheet Applications in Chemistry, p. 61
  7. Jurs, Computer Software Applications in Chemistry, Ch. 9
  8. Shaw and Tigg, Applied Mathematica, Section 5.3

 

How do we treat the problem of fitting data with a non-linear function, as occurs frequently in the sciences?  The first approach, which is used by some but NOT RECOMMENDED, is to apply a linearizing transformation to convert the problem into a linear least squares problem, which we know how to solve.  The problem with this approach is that frequently, we screw up the weighting factors and distort the fit in the process.  If you persist in this approach , the following graphic demonstrates the general approach that is needed.

 

Despite the possibility of using a linearizing transformation to change a non-linear function into a linear function, the PREFERRED ALTERNATIVE is to use a non-linear least squares fitting routine.  We might begin by asking where we encounter non-linear fitting  problems.  Some examples are:

·         Kinetics

·         Gas absorption

·         Fitting with arbitrary functions

·         Peak fitting

In addition, we will be discussing the technique of Simplex which aslos serves as a general introduction to the subject of optimization.  Optimization is used is:

Simplex

 

Simplex is an optimization technique in which one fins the best response in a multidimensional space, where the variables may be interrelated.  Simple is computationally slow but is quick and easy to implement.  The essence of the method is to observe the response of a system, change variables intelligently, not the new response, change variables, etc.  One is searching for a minimum or a maximum in a response surface.  The response surface is in n-dimensional space where the axes correspond to the variables of the problem.  A simplex is a geometric figure with N+1 vertices in an n-dimensional space.  To illustrate this concept, consider a two parameter space.

 

In this drawing, B represents the best vertex (as it is closest to the minimum), W the worst vertex and N the next best vertex.  The trick is to walk the simplex towards a minimum.  This is done by discarding W, finding a rotation point C where C is defined as (B+N)/2 and getting new vertex R where R=C + (C-W) on a line connecting W and C.  The procedure is shown graphically below.

 

 

Here are some of the tricks in doing this procedure>

·        If new vertex R has the worst response, discard N and turn the simplex.

·        Alter the size of the simplex to zero in on the minimum or to get “out of a hole”.

 

Chemical Applications

 

  1. Organic Synthesis  (J. Org. Chem 45, 2315 (1980).  In this synthesis the reaction yield was optimized with respect the the reactan concentrations, pH, T, solvent and time.
  2. Chromatography  (Anal. Chem 45, 278A (1973).  The resolution was optimized as a function of column temperature, and gas flow.
  3. Instrument Operation (Appl. Spectr. 36, 37 (1982).  The operation of an ICP spectrometer was optimized for 5 elements.
  4. Separation Procedure  (Anal. Chem 55, 1553 91983).  The extraction of trace lements was optimized with respect to temperature, acid ratios, bath time, etc.
  5. Non-linear least squares by Simplex.  Minimize  the following quantity

The following graphic shows a sample problem in radioactive decay

 

Computer Methods

 

FORTRAN 90

 

      program simplx

c.....

c.....simplex minimization of a function

c.....

      dimension c(10),e(10),p(10,10),r(10),x(10)

      dimension data(100,10)

c....enter initial information

      maxcnt=500

      errmin=1.0e-03

      nout=6

      ninp=5

      nsim=1

      write(nout,*)'  this is program simplex'

      write(nout,*)' enter number of observations'

      read(ninp,*)nobs

      write(nout,*)' enter number of variables (x+y)'

      read(ninp,*)nv

      write(nout,*)' enter x1,x2,....xp and y'

      do 6 i=1,nobs

 6    read(ninp,*)(data(i,j),j=1,nv)

      if (nobs.gt.nv) go to 11

      write(nout,9)

 9    format(' no. of obs. must be gt # of obs')

      stop

 11   write(nout,*)'enter number of parameters'

      read(ninp,*)np

      np1=np+1

      write(nout,*)' enter initial estimates of parameters'

      read(ninp,*) (x(i),i=1,np)

      e(1)=error(x,data,nv,nobs,kount)

      write(nout,15)e(1)

 15   format(' starting error function value',g12.4)

      write(nout,*)' want debug level output (y or n)'

      read(ninp,18)ans

 18   format(a1)

      idb=0

      if(ans.eq.'y')idb=1

c...   initialize the simnplex

      kount=0

      do 22 j=1,np

 22   p(1,j)=x(j)

      do 28 i=2,np1

      do 26 j=1,np

 26     p(i,j)=x(j)

      p(i,i-1)=1.1*x(i-1)

      if(abs(x(i-1)).lt.1.0e-12) p(i,i-1)=0.001

 28   continue

c...find plow and phigh  / best=plos/worst=phigh

 31   ilo=1

      ihi=1

      do 34 i=1,np1

      do 32 j=1,np

 32     x(j)=p(i,j)

      e(i)=error(x,data,nv,nobs,kount)

      if(e(i).lt.e(ilo)) ilo=i

      if(e(i).gt.E(ihi)) ihi=i

 34   continue

      write(nout,*) 'initial simnplex'

      do 40 k=1,np1

      write(nout,39)k,e(k), (p(k,j),j=1,np)

 39   format(3x,' vertex',i2,' eror and paramerters:',5f8.3)

 40   continue

c...find pnhi the next highest next=pnhi     

 41   nhi=ilo

      do 43 i=1,np1

      if(e(i).ge.e(nhi).and.i.ne.ihi)nhi=i

 43   continue

c..compute the centroid

      do 46 j=1,np

      c(j)=-p(ihi,j)

      do 44 i=1,np1

                  c(j)=c(j)+p(i,j)

 44     continue

      c(j)=c(j)/np

 46   continue

 51   continue

c... print current best vertex

      write(nout,53)kount,nsim

 53   format(' after', i3, 'error evaluations and',i3,' simplexes')

      write(nout,54)(p(ilo,j),j=1,np)

 54   format('  parameterm estimates:',5g12.4)

      write(nout,55)e(ilo)

 55   format('  error function:',g12.4)

c...stopping criterion

      if(kount.gt.maxcnt)stop

      if(abs(e(ilo)-e(ihi))/e(ilo).lt.errmin) go to 56

      go to 61

 56   write(nout,*)'  ===> error criterion satisfied'

      write(nout,54)(p(ilo,j),j=1,np)

      stop

c...reflection     

 61   do 62 j=1,np

      r(j)=1.9985*c(j)-0.9985*p(ihi,j)

 62   continue

      er=error(r,data,nv,nobs,kont)

      if(idb.gt.0)write(nout,65)er,(r(j),j=1,np)

 65   format(' reflection vertex', 3f10.5)

C...REFLECT AGIN IF SUCCESSFUL

      if(er.lt.e(ilo)) go to 91

      if(er.ge.e(ihi)) go to 122

c..replace worst vertex with new one

 79   do 80 j=1,np     

      p(ihi,j)=r(j)

 80   continue

      nsim=nsim+1

      e(ihi)=er

      if(er.gt.e(nhi)) go to 51

      ihi=nhi

      go to 41

c..expand the simplex

 91   ilo=ihi

      ihi=nhi

      do 93 j=1,np

      x(j)=1.95*r(j)-0.95*c(j)

 93   continue

      ex=error(x,data,nv,nobs,kount)

      if(ex.lt.er)go to 104

c..r better than x

      do 99 j=1,np

      p(ilo,j)=r(j)

 99   continue

      nsim=nsim+1

      e(ilo)=er

      go to 110

c.. x is better than r

 104  do 105 j=1,np     

      p(ilo,j)=x(j)

 105  continue

      if(idb.gt.0)write(nout,106)ex,(x(j),j=1,np)

 106  format(' expansion vertex',3f10.5)

      nsim=nsim+1

      e(ilo)=ex

 110  continue

      go to 41

c...  contract the simplex

 122  do 123 j=1,np     

      r(j)=0.5015*c(j)+0.4985*p(ihi,j)

 123  continue

      er=error(r,data,nv,nobs,kount)

      if(idb.gt.0) write(nout,124) er, (r(j),j=1,np)

 124  format(' contraction vertex', 3f10.5)

      if(er.lt.e(ilo)) go to 91

      if(er.lt.e(ihi)) go to 79

c...scale

      write(nout,135)

 135  format(' enter scale (<0 expands, >0 shrinks, 0=stop):')

      read(ninp,*)scal

      if(scal.eq.0.0)go to 999

 137  do 138 i=1,np1

      do 138 j=1,np

                  p(i,j)=p(i,j)+scal*(p(ilo,j)-p(i,j))

 138  continue

      go to 31

 999  stop

      end

c-----------------------------------------------------

      function error(x,data,nv,nobs,kount)

c---computes the error function for the data set

c------smaller value is better

      dimension x(10),data(100,10)

      error=0.0

      do 10 i=1,nobs

      yobs=data(i,nv)

c.....chnage the next statment to change the function being fit

       ycalc=x(1)*(1.0-exp(-x(2)*data(I,1)))

       resi=yobs-ycalc

       error=error+resi*resi

 10   continue

      kount=kount+1

      return

      end

 

 

 

 

 

Here is a sample execution of this program:

 

The advantages of the Simplex method include:

·        “short” code

·        will not diverge

·        no derivatives needed

The disadvantages of the Simplex method are:

In any case the results of one study ought to be tested for a local minima by restarting the simplex in a different place.

 

 

MORE EFFICIENT METHODS (THAN SIMPLEX) TO SEARCH PARAMETER SPACE.

 

 

 

Computer methods of searching parameter space.

 

 

 

The following example shows the use of these Mathematica functions.

Click here to see Mathematica example

 

 

 

 

The following demo shows the use of Excel to solve a nonlinear fitting problem with Solver and with a macro named SOLVSTAT.XLM

Click here to see Excel demo

 

 

 

FORTRAN 90

 

The typical Fortran non-linear least squares fitting procedure is implemented using the program MRQMIN and its driver xmrqmin.  Note you must supply the following information:

  1. the fitting function, identifying fitting parameters
  2. the number of parameters
  3. the derivatives of the function with respect to the parameters.