\section{The Levenberg-Marquardt Algorithm}\label{S:LevenbergMarquardt}
\setcounter{progress_ck_counter}{0}

\noi The \emph{Levenberg-Marquardt Algorithm}\index{Levenberg-Marquardt Algorithm} combines the Steepest Descent and Newton Methods in a manner
that emphasizes the first of the two when $\vecx_k$ is far from a local minimum
of $f$ and the second when $\vecx_k$ is close. The method was first developed in 1944 by K. Levenberg\index{Levenberg, K.} and later refined in the early 1960s by D. W. Marquardt\index{Marquardt, D.W.} \cite{Levenberg},\cite{Marquardt}.

%\subsection{Interpolating between the Steepest \mbox{Descent} and \mbox{Newton Methods}}
\subsection{Interpolating between the Steepest Descent and Newton
\newline Methods}

\noi Suppose $S \subseteq \R^n$ is open, $f:S \to \R$ is twice-differentiable, and that we wish to solve the unconstrained
NLP,
\begin{equation*}
\mbox{minimize } f(\vecx), \qquad \vecx \in S.
\end{equation*}
Assume $\vecx_0$ belongs to $S$. The \emph{Levenberg Method}\index{Levenberg Method} starts with the initial value $\vecx_0$ and, at each iteration, uses a \emph{trial iterate} formula given by
\begin{equation}\label{E:LevenbergMarquardtTrialIterateFormula}
\vecw=\vecx_k-\left[ H_f\left(\vecx_k \right)
+\lambda I_n\right]^{-1}\nabla f(\vecx_k) \qquad k=0,1,2,\ldots.
\end{equation}

\noi
In this formula, $I_n$ denotes the $n$-by-$n$ identity matrix, $\lambda$ is a
positive \emph{damping parameter}\index{Levenberg-Marquardt Algorithm!damping parameter}, whose
role is best understood by comparing extreme cases. If $\lambda=0$, then
(\ref{E:LevenbergMarquardtTrialIterateFormula}) reduces to the iteration formula for Newton's Method. If $\lambda$ is extremely large, then
\begin{equation*}
\vecw \approx \vecx_k - \lambda^{-1}\nabla f(\vecx_k),
\end{equation*}
so that $\vecw$ is essentially the first Steepest Descent Method iterate of $f$
starting from $\vecx_k$. Outside of these two extreme cases, $\lambda$
controls the extent to which each method contributes to the interpolation of the two.
Of course, we also assume $\lambda$ yields a matrix sum, $\ds H_f\left(\vecx_k \right)
+\lambda I_n$, that is invertible. This will certainly be the case if $\lambda$ is sufficiently
large, as then $\ds H_f\left(\vecx_k \right)
+\lambda I_n \approx \lambda I_n$.

\subsection{The Levenberg Method}

\noi The key to using (\ref{E:LevenbergMarquardtTrialIterateFormula}) effectively is to adjust
the damping parameter if the trial iterate has an objective value no less than that of $\vecx_k$.
We start with a relatively small value of $\lambda$ so that our interpolation formula
is biased toward Newton's Method and hence, in light of Theorem
\ref{Th:NewtonsMethodConvergenceRateTheorem}, is more likely to yield a sequence
converging to a minimum at a quadratic rate. If our
choice of $\lambda$ results in a value of $\vecw$, whose objective value is less than that of
$\vecx_k$, i.e., if $\ds f(\vecw) < f(\vecx_k)$, then we retain the trial iterate
and set $\vecx_{k+1}=\vecw$. From there, we decrease the damping parameter, $\lambda$, and calculate a new trial iterate, $\vecw$, where $\vecx_{k+1}$ replaces $\vecx_k$ in (\ref{E:LevenbergMarquardtTrialIterateFormula}). On the other hand, if
$\ds f(\vecw) \geq f(\vecx_k)$, then we increase $\lambda$ until we
obtain a new
trial iterate, $\vecw$, for which $\ds f(\vecw) < f(\vecx_k)$. This process is repeated
as necessary and leads to the method of Levenberg.  We now summarize the steps of the algorithm.

\newpar
{\underbar{The Levenberg Method}} \newline
To obtain an approximate solution of NLP (\ref{SteepestDescentStandardNLP}) under
the assumption $f$ is twice-differentiable on $S$:
\be
\item Start with an initial value $\vecx_0$ and a small, positive damping parameter $\lambda$.
\item For $k \geq 0$, calculate the trial iterate\index{Levenberg-Marquardt Algorithm!trial iterate},
\begin{equation}\label{E:LevenbergMarquardtLevenbergFormula}
\vecw=\vecx_k-\left[H_f\left(\vecx_k \right)+\lambda I_n\right]^{-1}\nabla f(\vecx_k).
\end{equation}
\item If $\ds f(\vecw) < f(\vecx_k)$, then set $\vecx_{k+1}=\vecw$, decrease $\lambda$, and return to
(2) with $\vecx_k$ replaced by $\vecx_{k+1}$.
\item If $\ds f(\vecw) \geq f(\vecx_k)$, then increase $\lambda$ and return to (2).
\item Repeat the process until a desired level of tolerance has been achieved, e.g., $\ds \left\| \nabla f(\vecx_k)\right\|$
is smaller than a specified value $\epsilon$.
\ee

\noi
The algorithm does not stipulate the extent to which we should decrease or increase $\lambda$ in (3) and (4), respectively. We return to this issue shortly. For now, mere experimentation with various positive values will suffice.

\newpar  %Maple: Cobb_Douglass_Marquardt_illustration_of_steps.mws
As an example, we now apply one iteration to the unconstrained \emph{ConPro Manufacturing Company} NLP,
\begin{multline}\label{E:LevenbergMarquardtCMC}
 \mbox{minimize } f(x_1,x_2 )= -1400 x_1^{\frac{1}{2}}x_2^{\frac{1}{3}} + 350x_1 +200x_2, \\
	\quad \mbox{ where }\vecx=\begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \in S=\R_+^2=\left\{(x_1,x_2) \mid x_1, x_2 >0 \right\},
\end{multline}
with initial value
$\ds \vecx_0=\begin{bmatrix} 10 \\ 10 \end{bmatrix}$ and $\lambda=.01$. Calculations establish that
$\ds f(\vecx_0)=-4038.1$, $\ds \nabla f (\vecx_0)=\begin{bmatrix} -126.904 \\ -117.93 \end{bmatrix}$,
and $\ds H_f(\vecx_0) = \begin{bmatrix} 23.845 & -15.897 \\ -15.897 & 21.196 \end{bmatrix}$. When these values are substituted
into the interpolation formula (\ref{E:LevenbergMarquardtLevenbergFormula}), we obtain
$\ds \vecw=\begin{bmatrix} 23.470 \\ 24.242 \end{bmatrix}$. Since $f(\vecw)=-7097.364 <f(\vecx_0)$, we set $\vecx_1=\vecw$,
decrease $\lambda$, and repeat the process with $\vecx_0$ replaced by $\vecx_1$.

\begin{progress_ck}
Perform the second iteration of the Levenberg Method for the unconstrained \emph{ConPro Manufacturing Company} NLP, using the previously obtained value
of $\vecx_1$, together with $\lambda=.005$.
\end{progress_ck}

\subsection{The Levenberg-Marquardt Algorithm}

\noi For large values of $\lambda$, the
Levenberg Method iterates depend primarily upon the gradient $\nabla f$. We can incorporate useful second derivative
information without introducing an excessive amount of extra work, insofar
as matrix inversion is concerned, if we replace
$I_n$ in (\ref{E:LevenbergMarquardtLevenbergFormula})
with the diagonal matrix
\begin{equation}\label{E:LevenbergMarquardtHessianDiagonal}
\mbox{Diag}\left(H_f\left(\vecx_k \right)\right)=
\begin{bmatrix}
\left.\frac{\partial^2 f}{\partial x_1^2 }\right|_{\vecx=\vecx_k} & 0 & \ldots & 0 \\
0 & \left.\frac{\partial^2 f}{\partial x_2^2}\right|_{\vecx=\vecx_k} & \ldots  & 0  \\
\vdots & & \ddots & \vdots  \\
0 & 0 & \ldots & \left.\frac{\partial^2 f}{\partial x_n^2}\right|_{\vecx=\vecx_k}
\end{bmatrix}.
\end{equation}

\noi
Note that $\ds \mbox{Diag}\left(H_f\left(\vecx_k \right)\right)$ utilizes
only those entries along the diagonal of $\ds H_f\left(\vecx_k \right)$.

\newpar
Substitution of (\ref{E:LevenbergMarquardtHessianDiagonal}) for
$I_n$ in (\ref{E:LevenbergMarquardtLevenbergFormula}) leads
to a new interpolation formula for computing
the trial iterate:

\begin{equation}\label{E:LevenbergMarquardtHessianTrialIterateFormula}
\vecw=\vecx_k-\left[ H_f\left(\vecx_k \right)
+\lambda\mbox{Diag}\left(H_f\left(\vecx_k \right)\right)\right]^{-1}\nabla  f(\vecx_k).
\end{equation}
 We now introduce a decision process for accepting the trial iterate, one that is similar to that used for the Levenberg Method, but that depends upon both an \emph{initial damping parameter}, $\lambda_0$, together
 with a \emph{scaling factor}\index{Levenberg-Marquardt Algorithm!scaling factor}, $\rho>1$. At each iteration, we compute the trial iterate
 using the previous iterate, $\vecx_k$, along with a value of $\lambda$ that depends upon $\lambda_0$, $\rho$, and $k$. For many applications, $\lambda_0=.001$ and $\rho=10$
 will suffice \cite{Bevington}.
This new decision process leads to the Levenberg-Marquardt Algorithm.

\newpar
{\underbar{The Levenberg-Marquardt Algorithm}}\index{Levenberg-Marquardt Algorithm!summary of steps} \newline
To obtain an approximate solution of NLP (\ref{SteepestDescentStandardNLP}) under
the assumption $f$ is twice-differentiable on $S$:
\be
%1
\item Start with an initial value, $\vecx_0$, in $S$, an initial damping parameter, $\lambda_0$,
and a scaling parameter, $\rho$. For $k \geq 0$ do the following:
%2
\item Determine a trial iterate, $\vecw$, using
(\ref{E:LevenbergMarquardtHessianTrialIterateFormula}) with
$\lambda=\lambda_k \rho^{-1}$.
%3
\item If $f(\vecw)< f(\vecx_k)$, where $\vecw$ is determined in (2), then set $\vecx_{k+1}=\vecw$ and
$\lambda_{k+1}=\lambda_k \rho^{-1}$. Return to (2), replace $k$ with $k+1$, and compute a new trial iterate.
%4
\item If $f(\vecw)\geq f(\vecx_k)$ in (3), determine a new trial
iterate, $\vecw$, using (\ref{E:LevenbergMarquardtHessianTrialIterateFormula}) with
$\lambda=\lambda_k$.
%5
\item If $f(\vecw)< f(\vecx_k)$, where $\vecw$ is determined in (4), then set $\vecx_{k+1}=\vecw$ and
$\lambda_{k+1}=\lambda_k$. Return to (2), replace $k$ with $k+1$, and compute a new trial iterate.
%6
\item If $f(\vecw) \geq f(\vecx_k)$ in (5), then determine the smallest value of $m$ so that
when a trial iterate, $\vecw$, is computed using (\ref{E:LevenbergMarquardtHessianTrialIterateFormula}) with $\lambda=\lambda_k \rho^m$, then $f(\vecw)< f(\vecx_k)$. Set $\vecx_{k+1}=\vecw$
and $\lambda_{k+1}=\lambda_k \rho^m$. Return to (2), replace $k$ with $k+1$, and compute
a new trial iterate.
\item Terminate the algorithm when $\ds \left\|  \nabla f\left(\vecx_k\right)\right\| <\epsilon$,
where $\epsilon$ is a specified tolerance.
\ee

\newpar An important observation to make
is that in step (3), by setting $\lambda_{k+1}=\lambda_k \rho^{-1}$, we are emphasizing to a greater extent the Newton Method component in the
interpolation formula (\ref{E:LevenbergMarquardtHessianTrialIterateFormula}), with the hope that the sequence of iterates converges to the minimum
at a quadratic rate. If our decision-making process takes us to (6), then
setting $\lambda_{k+1}=\lambda_k \rho^{m} > \lambda_k$ signifies we are forced to emphasize the Steepest
Descent Method component in (\ref{E:LevenbergMarquardtHessianTrialIterateFormula}).

\newpar
An application of the Levenberg-Marquardt Algorithm to the
unconstrained \emph{ConPro Manufacturing Company} NLP
(\ref{E:LevenbergMarquardtCMC}) with $\ds \vecx_0=\begin{bmatrix} 10 \\ 10 \end{bmatrix}$, $\lambda_0=.001$, $\rho=10$, and tolerance $\epsilon =.01$ yields results
shown in Table \ref{Ta:LevenbergMarquardtCMCResults}.

\begin{table}[h]
\caption{Results of Levenberg-Marquardt Algorithm applied to the unconstrained \emph{ConPro Manufacturing Company} NLP with
a tolerance of $\epsilon =.01$, $\lambda=.0001$, and $\rho=10$}
\label{Ta:LevenbergMarquardtCMCResults}
\begin{tabular}{|c|c|c|}
\hline
$k$ & $\vecx_k^t$ &   $\ds \left\|  \vecx_k-\vecx_{\bigstar} \right\|$ \\
\hline
					0 & $[10,10]$ &  119.759 \\ \hline
                1  & $[28.057, 29.105]$  & 93.527 \\ \hline
                 2 & $[56.811,67.711]$  & 51.116396 \\ \hline
                 3 & $[80.072,91.034]$  & 12.721 \\ \hline
                4 & $[86.736,100.920]$  &  .802 \\ \hline
                5 & $[87.111,101.626]$  &  0.004 \\ \hline
\end{tabular}
\end{table}

\noi
Table \ref{Ta:SteepestDescent_CMC_SteepestDescent} indicates that eight
iterations of the Steepest Descent Method are required to achieve the same level of tolerance. Newton's Method only requires five. Thus, for this example the Levenberg-Marquardt Algorithm requires fewer iterations than the Steepest Descent Method and a number equal to that of Newton's Method. It should not come as a surprise then that, under appropriate conditions, this new algorithm exhibits quadratic
convergence\index{Levenberg-Marquardt Algorithm!quadratic convergence rate} \cite{Yamashita}.

%Specifically,
%if the sequence of Levenberg-Marquardt iterates $\ds \left\{ \vecx_k \right\}$
%converges to a local minimum, $\vecx_{\bigstar}$, and
%if $\ds H_f\left(\vecx_{\bigstar}\right)$ is invertible, then this sequence
%converges to $\vecx_{\bigstar}$ at a quadratic rate. (See for example,
%\cite{Yamashita}).

\subsection{A Maple Implementation of the \mbox{Levenberg-Marquardt} \mbox{Algorithm}}\index{Levenberg-Marquardt Algorithm!Maple implementation}

\noi The \mf{LevenbergMarquardt} procedure found in the worksheet \textbf{Levenberg-Marquardt Algorithm.mw} functions in a manner similar to the
\mf{SteepestDescentMethod} and \mf{NewtonsMethod} procedures
from Sections \ref{S:SteepestDescent} and \ref{S:NewtonsMethod},
respectively. However, we must take into account different choices of
the damping parameter and scaling factor. Our procedure call takes the form
\begin{equation*}
\mbox{\mf{LevenbergMarquardt(\ma{function,initial,N,tolerance,damp,scale})},}
\end{equation*}
where \ma{function} is a twice-differentiable expression to be minimized, \ma{initial} is the initial value, $\ds \vecx_0$, (written as a
list), \ma{N} is the maximum number of iterations, \ma{tolerance} is the specified tolerance,
\ma{damp} is the initial damping parameter, and \ma{scale} is the scaling factor. The algorithm terminates when either the number of iterations is reached or the gradient norm is less than the desired tolerance. The returned value is the approximate minimum value. Intermediate
iterates leading to the approximate minimum are stored in a list of points, \mf{iterates}, a global variable accessible outside of the procedure. In the procedure, we utilize an auxiliary function, $d$,
that returns formula (\ref{E:LevenbergMarquardtHessianTrialIterateFormula}) in terms
of an input value, $\lambda$, at each iteration.

\newpar
The following worksheet constructs the procedure and approximates the minimum of the \CM objective,
\newline
$\ds f(x_1,x_2 )= -1400 x_1^{\frac{1}{2}}x_2^{\frac{1}{3}} + 350x_1 +200x_2$
using $\ds \vecx_0=\begin{bmatrix} 10 \\ 10 \end{bmatrix}$, $N=5$, a tolerance
of $\epsilon=.01$, a damping parameter of .001, and a scaling factor of 10.

\begin{maplelist}
\item \mf{restart: with(VectorCalculus):with(LinearAlgebra):}
\item \mf{LevenbergMarquardt := proc(function, initial, N, tolerance,
\newline damp, scale)
local x,f,d,Delf,Hf, DiagH,lambda,rho,w,j,m,epsilon; \newline
global iterates: \newline
\# Create local and global variables. \newline
x:=array(0 .. N): \# Define array of iterates. \newline
f:=unapply(function, [x1, x2]): \newline
 \# Create function. \newline
Delf:=unapply(Gradient(function, [x1, x2]), [x1, x2]): \newline
\# Create corresponding gradient function. \newline
Hf:=unapply(Hessian(function, [x1, x2]), [x1, x2]): \newline
\# Create corresponding Hessian matrix function. \newline
x[0]:=evalf(initial):
\newline \# Set initial array value.\newline
epsilon := 10:
\newline \# Set initial epsilon to a large value. \newline
lambda[0]:=damp: rho:=scale: \newline
\# Set initial damping parameter and scaling factor. \newline
j:=0: while (tolerance <= epsilon and j<=N-1) do \newline
\# Create loop structure to perform algorithm. \newline
DiagH := DiagonalMatrix(Diagonal(Hf(x[j][1], x[j][2]))):\newline
\# Create diagonal matrix using Hessian entries. \newline
d :=r-> convert(-MatrixInverse(Hf(x[j][1], x[j][2]) \newline
+r*DiagH).Delf(x[j][1], x[j][2]), list) \newline
\# Formulate auxiliary function. \newline
w := x[j]+d(lambda[j]/rho); \newline
\# Create trial solution.\newline
Test for acceptance of first trial solution.\newline
if f(w[1], w[2]) < f(x[j][1], x[j][2]) then x[j+1] := w:\newline
lambda[j+1] := lambda[j]/rho \newline
else w := x[j]+d(lambda[j]); \newline
\# Test for acceptance of second trial solution. \newline
if f(w[1], w[2]) < f(x[j][1], x[j][2]) then x[j+1] := w: \newline
lambda[j+1]:=lambda[j] \newline
else w := x[j]; \newline
\# If necessary, scale damping factor \newline
until a suitable trial solution found. \newline
m:=1:while f(x[j][1], x[j][2]) <= f(w[1], w[2]) do \newline
w:=x[j]+d(lambda[j]*rho\^{}m): \newline
m:= m+1 end do; x[j+1] := w: lambda[j+1]:=lambda[j]*rho\^{}m \newline
end if: end if; \newline
epsilon := evalf(Norm(Delf(x[j+1][1], x[j+1][2]), Euclidean)); \newline
\# Update epsilon using gradient.\newline
j:=j+1: \newline
\# Increase loop index. \newline
end do: \newline
iterates:=[seq(x[i],i=0..j)]:RETURN(x[j]):end:
\newline
\# Construct iterates and return last array value.}
\item \mf{f:=(x1,x2)->-1400x1\^{}(1/2)x2\^{}(1/3)+350x1+200x2;
\newline \# Enter function}
\mo{f:=(x_1,x_2) \rightarrow -1400x_1^{\frac{1}{2}}x_2^{\frac{1}{3}}
+350x_1+200x_2}
\item \mf{LevenbergMarquardt(f(x1, x2), [10,10],5,.01,.001,10);\newline
\# Procedure call using initial value [10,10], bound of 
\newline 5 iterations, tolerance, .00, damping parameter of .001, 
\newline and scaling factor of 10.}
\mo{ [87.10975285,101.6261946]}
\item \mf{iterates;
\newline \# Print iterates.}
\begin{align*}
&[[10., 10.], [28.00148014, 29.04621145], [56.73574951, 61.62388852],\\
 &[80.03451430, 90.98472144], [86.73264262, 100.9143105], \\
 &[87.10972680, 101.6261320]]
\end{align*}
\end{maplelist}

\subsection{Nonlinear Regression}\index{regression!nonlinear}

\noi Levenberg's original intent was to develop an algorithm well-suited
for solving nonlinear least-squares problems. Specifically, given a set of data pairs,
$\ds
\left\{ (\vecx_1,y_1),(\vecx_2,y_2),\ldots,(\vecx_N,y_N) \right\}$, where each $y_k$ is real and each $\vecx_k$
belongs to $\R^n$, we wish to determine a function, $f$, of a certain type that minimizes the \emph{sum of squared-errors}\index{sum of squared-errors},
\begin{equation}\label{E:LevenbergMarquardtNonlinearRegressionFormula}
\sum_{k=1}^N \left(f(\vecx_k)-y_k\right)^2.
\end{equation}
Generally, we assume that $f$ belongs to a family of functions that are related
by one or more parameters. Multiple linear regression, as discussed in Section
\ref{SS:MultipleLinearRegression}, is a particular example.

\newpar
Unlike the case of multiple linear regression, where (\ref{E:LevenbergMarquardtNonlinearRegressionFormula}) results in an unconstrained
NLP whose critical points are determined analytically, minimization of
(\ref{E:LevenbergMarquardtNonlinearRegressionFormula})
in the nonlinear setting frequently requires a numeric algorithm\index{Levenberg-Marquardt Algorithm!nonlinear regression}.

\newpar
Such was the case for the \emph{ConPro Manufacturing Company} when it constructed its production function, which was based upon the data shown in Table \ref{Ta:LevenbergMarquardtCMCProductionData}.

\begin{table}[h]
\caption{Production data for \CMC}
\label{Ta:LevenbergMarquardtCMCProductionData}
\begin{tabular}{|c|c|c|}
\hline
$x_1$ &  $x_2$ &  $P(x_1,x_2)$  \\
\hline
60&60&29\\ \hline
60&80&36\\ \hline 60&100&34\\ \hline
60&120&37\\ \hline 80&60&36\\ \hline
80&80&38\\ \hline 80&100&41\\ \hline
80&120&47\\ \hline 100&60&35\\ \hline
100&80&45\\ \hline 100&100&47\\ \hline
100&120&47\\ \hline 120&60&46\\ \hline
120&80&48\\ \hline 120&100&47\\ \hline
120&120&54\\ \hline
\end{tabular}
\end{table}

\newpar
Recall $x_1$ and $x_2$ denote the number of units of material and labor, respectively, used to produce $P(x_1,x_2)$ units of pipe.

\newpar
A traditional choice from economic theory for modeling production is the
\emph{Cobb-Douglas Production Function}\index{Cobb-Douglas Production Function},
\begin{equation}\label{E:LevenbergMarquardtCobbProductionFunction}
P(x_1, x_2) = x_1^{\alpha}x_2^{\beta}, \mbox{ where } \alpha,\beta >0.
\end{equation}

\noi
To determine the Cobb-Douglas function that best fits the data in Table
\ref{Ta:LevenbergMarquardtCMCProductionData}, we minimize the \emph{sum of squared-errors}\index{sum of squared-errors}
\begin{equation}\label{E:LevenbergMarquardtProductionDataObjectiveFunction}
f(\alpha,\beta) = \sum_{k=1}^{16}\left(x_1^{\alpha}x_2^{\beta}-P(x_1,x_2) \right)^2, \qquad \alpha, \beta >0,
\end{equation}
where the sum is over the sixteen ordered triples
$\left( x_1,x_2,P(x_1,x_2) \right)$ in Table
\ref{Ta:LevenbergMarquardtCMCProductionData}.
If we apply the Levenberg-Marquardt
Algorithm, using an initial value
$\ds \vecx_0=\begin{bmatrix} \alpha_0 \\ \beta_0 \end{bmatrix}=\begin{bmatrix} .5 \\ .5 \end{bmatrix}$,
initial damping factor $\lambda_0=.001$, a scaling factor of $\rho=10$, and a tolerance of
$\epsilon=.001$, we obtain after five
iterates, $\alpha \approx .510$ and $\beta \approx .322$. These values are quite close to those used by \CM in its
actual production formula.

\newpar
Figure \ref{Fi:f08x03x01} illustrates a plot of the data
triples $(x_1,x_2,P)$ from Table \ref{Ta:LevenbergMarquardtCMCProductionData} along with
the ``best-fit'' production function $\ds P(x_1,x_2)=x_1^{.510}x_2^{.320}$.

\begin{figure}[h]
\centering
\psfrag{x1_axis_label}[Bc][Bc]{$x_1$}
\psfrag{x2_axis_label}[Bc][Bc]{$x_2$}
\psfrag{f_axis_label}[Bc][Bc]{$P$}
\includegraphics[width=.7\linewidth]{images/f08x03x01.eps}
\caption{Plot of \emph{ConPro Manufacturing Company} production data, together with best-fit Cobb-Douglass function $\ds P(x_1,x_2)=x_1^{.510}x_2^{.320}$.}
\label{Fi:f08x03x01}
\end{figure}

\newpar
An important quantity associated with any regression problem is the
\emph{coefficient of determination}\index{coefficient of determination}. This quantity lies between 0 and 1 and
measures the proportion of variation of the data accounted for by the best-fit function.
If the data points used for the regression are given
by $\ds
\left\{ (\vecx_1,y_1),(\vecx_2,y_2),\ldots,(\vecx_m,y_m) \right\}$ and if ${\overline{y}}$ denotes the
mean of $\{y_1,y_2,\ldots, y_m \}$, this proportion then becomes
 \begin{equation}\label{E:LevenbergMarquardtCoefficientOfDetermination}
\frac{ \sum_{k=1}^m \left(f(\vecx_k)-y_k\right)^2 }{\sum_{k=1}^m \left({\overline{y}}-y_k\right)^2} .
\end{equation}

\noi
Direct calculation using (\ref{E:LevenbergMarquardtCoefficientOfDetermination}) shows that the production function
\newline
$\ds P(x_1,x_2)=x_1^{.510}x_2^{.320}$ fits the data in Table \ref{Ta:LevenbergMarquardtCMCProductionData}
with a coefficient of determination equal to .94, thus indicating a strong goodness of fit.


\subsection{Maple's \mf{Global Optimization Toolbox}}

\noi
The Steepest Descent, Newton, and Levenberg-Marquardt Methods are only three of many
numeric algorithms used to approximate solutions of unconstrained NLPs. Their use is limited, however,
in that they offer
no guarantee of locating global solutions, local solutions, or even critical points for that matter.

\newpar
There is a great deal of ongoing research in the field of global optimization, and several techniques, which are beyond the scope of this text,
have been developed for solving unconstrained NLPs.

\newpar
One specific tool developed during the last decade
is known as \emph{Lipschitz Global Optimization}\index{Lipschitz Global Optimization}, or LGO \cite{Pinter}. Implementation of LGO for various software programs, e.g., Maple and Excel, is now
available as an add-on package. In Maple, for example, the \mf{Global Optimization Toolbox} package is loaded using the
command \mf{with(GlobalOptimization)}, and the global minimum is computed using the command
\mf{GlobalSolve}\index{Maple!\mf{GlobalSolve}}, whose arguments essentially take the same form as
those used for the \mf{LPSolve} and \mf{NLPSolve} commands.

\newpar
A brief review of other various global optimization software packages may be found in \cite{Pinter2}. Three texts devoted
to the general topic of global optimization include \cite{Floudas}, \cite{Horst}, \cite{Neumaier}, and \cite{Pinter}.

\section*{Exercises Section \protect\ref{S:LevenbergMarquardt}}

\addcontentsline{toc}{subsection}{Exercises}

\input{Exercises/Exercises_Levenberg_Marquardt}
\endinput
