newmatnl.cpp

矩阵计算库

//$$ newmatnl.cpp         Non-linear optimisation

// Copyright (C) 1993,4,5: R B Davies


#define WANT_MATH
#define WANT_STREAM

#include "newmatap.h"
#include "newmatnl.h"




void FindMaximum2::Fit(ColumnVector& Theta, int n_it)
{
   Tracer tr("FindMaximum2::Fit");
   Real z,w,x,x2,g,l1,l2,l3,d1,d2,d3;
   ColumnVector Theta1, Theta2, Theta3;
   int np = Theta.Nrows();
   ColumnVector H1(np), H3, HP(np), K, K1(np);
   Boolean oorg, conv;
   int counter = 0;

   {                              // subblock so won't have labels in same
                                  // block as destructors.

                                  // I know it is supposed to be evil to
                                  // use "goto". Nevertheless this seems
                                  // to be the best way of coding the
                                  // algorithm, even in C++.

      Theta1 = Theta; HP = 0.0; g = 0.0;

   Start:
      tr.ReName("FindMaximum2::Fit/Start");
      Value(Theta1, TRUE, l1, oorg); if (oorg) Throw(DataException(Theta));

   Restart:
      tr.ReName("FindMaximum2::Fit/ReStart");
      conv = NextPoint(H1, d1); if (conv) goto Convergence;
      if (counter++ > n_it) goto Fail;

      z = 1.0 / sqrt(d1);
      H3 = H1 * z; K = (H3 - HP) * g; HP = H3;
      g = 0.0;                     // de-activate to use curved projection
      if (g==0.0) K1 = 0.0; else K1 = K * .2 + K1 * .6;
      // (K - K1) * alpha + K1 * (1 - alpha) = K * alpha + K1 * (1 - 2 * alpha)
      K = K1 * d1; g = z;

   Continue:
      tr.ReName("FindMaximum2::Fit/Continue");
      Theta2 = Theta1 + H1 + K;
      Value(Theta2, FALSE, l2, oorg);
      if (counter++ > n_it) goto Fail;
      if (oorg)
         { H1 = H1 * 0.5; K = K * 0.25; d1 *= 0.5; g *= 2.0; goto Continue; }
      d2 = LastDerivative(H1 + K * 2.0);

   Interpolate:
      tr.ReName("FindMaximum2::Fit/Interpolate");
      z = d1 + d2 - 3.0 * (l2 - l1);
      w = z * z - d1 * d2;
      if (w < 0.0) goto Extrapolate;
      w = z + sqrt(w);
      if (1.5 * w + d1 < 0.0) goto Extrapolate;
      if (d2 > 0.0 && l2 > l1 && w > 0.0) goto Extrapolate;
      x = d1 / (w + d1); x2 = x * x; g /= x;
      Theta3 = Theta1 + H1 * x + K * x2;
      Value(Theta3, TRUE, l3, oorg);
      if (counter++ > n_it) goto Fail;
      if (oorg)
      {
         if (x <= 1.0)
            { x *= 0.5; x2 = x*x; g *= 2.0; d1 *= x; H1 = H1 * x; K = K * x2; }
         else
         {
            x = 0.5 * (x-1.0); x2 = x*x; Theta1 = Theta2;
            H1 = (H1 + K * 2.0) * x;
            K = K * x2; g = 0.0; d1 = x * d2; l1 = l2;
         }
         goto Continue;
      }

      if (l3 >= l1 && l3 >= l2) { Theta1 = Theta3; l1 = l3; goto Restart; }

      d3 = LastDerivative(H1 + K * 2.0);
      if (l1 > l2)
         { H1 = H1 * x; K = K * x2; Theta2 = Theta3; d1 = d1*x; d2 = d3*x; }
      else
      {
         Theta1 = Theta2; Theta2 = Theta3;
         x -= 1.0; x2 = x*x; g = 0.0; H1 = (H1 + K * 2.0) * x;
         K = K * x2; l1 = l2; l2 = l3; d1 = x*d2; d2 = x*d3;
         if (d1 <= 0.0) goto Start;
      }
      goto Interpolate;

   Extrapolate:
      tr.ReName("FindMaximum2::Fit/Extrapolate");
      Theta1 = Theta2; g = 0.0; K = K * 4.0; H1 = (H1 * 2.0 + K);
      d1 = 2.0 * d2; l1 = l2;
      goto Continue;

   Fail:
      Throw(ConvergenceException(Theta));

   Convergence:
      Theta = Theta1;

   }
}



void NonLinearLeastSquares::Value
   (const ColumnVector& Parameters, Boolean, Real& v, Boolean& oorg)
{
   Tracer tr("NonLinearLeastSquares::Value");
   Y.ReDimension(n_obs); X.ReDimension(n_obs,n_param);
   // put the fitted values in Y, the derivatives in X.
   Pred.Set(Parameters);
   if (!Pred.IsValid()) { oorg=TRUE; return; }
   for (int i=1; i<=n_obs; i++)
   {
      Y(i) = Pred(i);
      X.Row(i) = Pred.Derivatives();
   }
   if (!Pred.IsValid()) { oorg=TRUE; return; }  // check afterwards as well
   Y = *DataPointer - Y; Real ssq = Y.SumSquare();
   errorvar =  ssq / (n_obs - n_param);
   cout << "\n" << setw(15) << setprecision(10) << " " << errorvar;
   Derivs = Y.t() * X;          // get the derivative and stash it
   oorg = FALSE; v = -0.5 * ssq;
}

Boolean NonLinearLeastSquares::NextPoint(ColumnVector& Adj, Real& test)
{
   Tracer tr("NonLinearLeastSquares::NextPoint");
   QRZ(X, U); QRZ(X, Y, M);     // do the QR decomposition
   test = M.SumSquare();
   cout << " " << setw(15) << setprecision(10)
      << test << " " << Y.SumSquare() / (n_obs - n_param);
   Adj = U.i() * M;
   if (test < errorvar * criterion) return TRUE;
   else return FALSE;
}

Real NonLinearLeastSquares::LastDerivative(const ColumnVector& H)
{ return (Derivs * H).AsScalar(); }

void NonLinearLeastSquares::Fit(const ColumnVector& Data,
   ColumnVector& Parameters)
{
   Tracer tr("NonLinearLeastSquares::Fit");
   n_param = Parameters.Nrows(); n_obs = Data.Nrows();
   DataPointer = &Data;
   FindMaximum2::Fit(Parameters, Lim);
   cout << "\nConverged\n";
}

void NonLinearLeastSquares::MakeCovariance()
{
   if (Covariance.Nrows()==0)
   {
      UpperTriangularMatrix UI = U.i();
      Covariance << UI * UI.t() * errorvar;
      SE << Covariance;                 // get diagonals
      for (int i = 1; i<=n_param; i++) SE(i) = sqrt(SE(i));
   }
}

void NonLinearLeastSquares::GetStandardErrors(ColumnVector& SEX)
   { MakeCovariance(); SEX = SE.AsColumn(); }

void NonLinearLeastSquares::GetCorrelations(SymmetricMatrix& Corr)
   { MakeCovariance(); Corr << SE.i() * Covariance * SE.i(); }

void NonLinearLeastSquares::GetHatDiagonal(DiagonalMatrix& Hat) const
{
   Hat.ReDimension(n_obs);
   for (int i = 1; i<=n_obs; i++) Hat(i) = X.Row(i).SumSquare();
}