# Presentation: 045446

## V. Nonlinear Regression By Modified Gauss-Newton Method: Theory

- Method to calculate model parameter estimates that result in the best fit, in a least squares sense, to the calibration data.
- Implemented as an iterative form of linear regression
- Chapter 5 and Appendix A (Hill and Tiedeman, 2007)

## Outline

- Linear regression and normal equations
- Nonlinear regression and normal equations
- Modifications that make the nonlinear regression normal equations work in practice
- Stopping criteria
- Limits on estimated parameters
- Log transformed parameters
- Exercises 5.2a and 5.2b
- Addition of prior information
- Exercise 5.2c

## Simple Linear Regression

- Suppose we collect some data and want to determine the relation between the observed values, y, and the independent variable, x:

- observed
- response

- true,
- unknown
- intercept

- true,
- unknown
- slope

- true,
- unknown
- random error

- If we think the relation between y and x is linear, we can model the data using a **linear model**:

- 0 and 1 are the **true parameters** of this linear model.

## Simple Linear Regression

- Don’t know the true values of the parameters.
- Estimate them using the assumed model and the observations, by expressing the linear model in terms of estimated parameter values:

- observed
- response

- Estimate b0 and b1 to obtain the best fit of the simulated values to the observations.
- One method: **Minimize sum of squared errors, or residuals**.

- estimate  of  0

- estimate  of  1

- residual

- simulated values, *y**i*

## Simple Linear Regression

- Sum of squared residuals:

- To minimize:

- This results in the **normal equations**:

- Solve these equations to obtain expressions for b0 and b1, the parameter estimates that give the best fit of the simulated and observed values.
- If have more than 2 parameters, need to replace summations with matrix notation. First, look at using 2 parameters and matrix notation.

**Notes:**
Remember from calculus: How do we find the minimum of a function?

Reason for ‘normal’ terminology: At the least squares estimates b’, vectors y’=Xb’ and e=y-Xb’ are ORTHOGONAL, or ‘normal’ to eachother.

Normal equations: 2 unknowns (b0,b1) and 2 equations. Thus, can solve by hand for b0 and b1 just using a calculator and some algebra.

## Linear Regression in Matrix Form

- Linear regression model:                                   , i=1.n   

- In general, _*X*_ is of dimension ND x NP, and _*b*_ is of dimension NP, where ND is the number of observations and NP is the number of parameters.
- The **normal equations** (_*b*_’ is the vector of least-squares estimates of _*b*_):

- vector of
- observed
- values

- matrix of
- coefficients

- vector of
- residuals

- vector of
- parameters

- 

- Using summa-tions

- Using matrix notation:

**Notes:**
Now we will generalize the previous discussion of linear regression.

Y=Xb+e: more general form of linear model.

Definitions of y,X,etc: just writing out equation for many values of i, using vectors and matrices. Vectors are of length ND = # obs.

If we calculate the sum of squared residuals using the expression y=Xb+e, and minimize it, we get the normal equations and their solution shown here.

Show what Transpose is.

## Linear Regression with Weighting

- When using weights, we minimize the sum of squared weighted residuals:

- We again set the derivative of the objective function to zero:

- This leads to the **normal equations****:**

- Solving for _*b*_’:

**Notes:**
Obj func written as (w^(1/2) e )^2 assumes using a DIAGONAL weight matrix and that the different observations are uncorrelated.

Writing as (eT w e) allows for a full weight matrix, with terms on the off-diagonals. This allows for a weighting scheme to account for correlations among different observations.

Example: streamflow measured at 3 points. The gain/loss over the 2 reaches are correlated because each is calculated using the middle point measurement.

One study in gw lieterature that has incorporated the off diagonal terms found that including correlations among the observations didn’t matter in terms of the regression results.

In mf2k, can have full weight matrix on flows, adv, prior.

## Linear versus Nonlinear Models

**Linear models**: Sensitivities of matrix _*X*_ are **not** a function of the model parameters:

- Linear models have elliptical objective function surfaces. With two parameters:

*One step to get to the minimum.*

**Notes:**
First eqn: Simulated values for a 2-parameter linear model

Terms of X are independent of values of b0 and b1.
These are sensitivities: derivative of a simulated value w/r/t a model parameter. Same quantities you worked with in ex. 6.
dS/db = f (dy/db), and dy/db are constant for linear model; therefore dS/db is constant no matter what the values of b are.

Obj func surface shown: Nice shape – well defined global minimum, no local minima.
For linear problems, get to minimum of this surface (to optimal values of b0,b1) in just 1 step of solving normal equations.

## Linear versus Nonlinear Models

**Nonlinear models**: Sensitivities of matrix _*X*_ **are** a function of the model parameters.
- Ground-water models are commonly nonlinear. This nonlinearity comes from Darcy’s Law:

**Q = -KA (********h/********x)**
**h = h****1**** – x (Q/KA)**

- Derivatives:

******h/********x**** = -Q/KA ****	**Not a function of x. Linear in x.
******h/********Q**** = -x/****K****A ****	**Function of K. Linear in Q, but nonlinear in K.
******h/********K**** = x****Q****/****K****2****A****	**Function of K *and* Q. Nonlinear in both K and Q.

**Notes:**
Used Darcy’s Law before for illustrating parameter correlations, by solving for h.
Now, use to illustrate model nonlinearity, by taking derivatives w/r/t terms of the equations.

dh/dQ result – Q doesn’t appear in RHS.

## Nonlinear Regression

**An iterative form of linear regression, with some modifications to the normal equations to make them work in practice**. General procedure:
- Linearize the model around the current parameter values. This results in a linearized objective-function surface.
- Using the normal equations, calculate new parameter values that are closer to the minimum of the linearized objective-function surface, and therefore, hopefully closer to the minimum of the nonlinear objective-function surface.
- Repeat from step 1.

**Notes:**
Now sensitivities are different for different parameter values. Thus objective function surface is different. Not just value of S being different for different parameters, but whole surface in all of parameter space.

## Nonlinear Regression

- Nonlinear model:

_*y*_’(_*b*_) is nonlinear function of _*b*_.
- To linearize _*y*_’(_*b*_), we use a Taylor Series expansion about the current values of b.
- For a one-parameter model:

**Notes:**
Can’t write in terms of X anymore because X changes with b.

Superscript l: linearized.
This equation calculates the linearized value of y’i at b, by linearizing the function y’ about b0.

First term on RHS: actual simulated value at b0
Second term: is delta-y – [ (delta y)/(delta b) ] (delta b)

## Nonlinear Regression

- One-parameter model example: **The Theim equation (pg 70)**.
- drawdown=[Q/(2**T**)] ln(r/r0)

- Models
- Nonlinear and
- Linearized

- Objective Functions
- Nonlinear and
- Linearized

**Notes:**
Theim model: steady flow to pumped will in confined aquifer.
h=f(T), recall from Darcy’s law that h = nonlinear function of K, or T here. Thus, nonlinear model.

Upper graph, black curve: drawdown (head) plotted against T. See how is nonlinear.
Red curve: linearized model. Gives exact value at T=.005, deviates more & more from exact model as we get further and further from T=.005

Lower graph: Assumed some drawdown data, to construct an objective function. 
Red curve: is exact at T=.005, deviates significantly away from T=.005.

As move closer from T=.005 to minimum of linearized objective function, we are also moving closer to minimum of nonlinear function. But, bad to move too far!.

## Nonlinear Regression

- Two-parameter model example: **The Theis equation (pg 75)**.

- True minimum

- Points about which the surface is linearized

**Notes:**
Dot is true minimum.

X’s are points we will linearize around.

Lower left figure: We would overshoot if we went all the way to the minimum of the linearized surface.

## Nonlinear Regression

- Objective function for nonlinear model:

- To minimize S, we first substitute the linearized approximation of _*y*_’(_*b*_) into S. The linear approximation written in terms of **sensitivity matrix **_***X***_:

- Substituting this into S, and setting

- leads once again to

- The **NORMAL EQUATIONS!!!**

** ****(**_**X**_**T **_****_** **_**X**_**)   **_**d**_**  =****  **_**X**_**T**** **_****_** (**_**y**_** - **_**y**_******* *****)**** **** **

**Notes:**
In S(b) eqn, we write y’(b) because Xb is not constant.

b0 are starting values for first iteration, then are updated values in subsequent iterations.

In these normal equations: y is vector of obs. y’ is calculated at vector b=b0.

This shouldn’t look too mysterious!

## Nonlinear Regression

** ****(**_**X**_**T **_****_** **_**X**_**)   **_**d**_**  =****  **_**X**_**T**** **_****_** (**_**y**_** - **_**y**_******* *****)**** **** **

- The **NORMAL EQUATIONS:**

- Same form as the normal equations for linear regression, except:
_**d**_ replaces _b_
_**y**_** **-** **_**y**_**** replaces _y_
- These differences are because of the iterative process used to solve the nonlinear regression.
-  _**d**_ is the parameter change vector: _**b**_**r+1**** **=** **_**b**_**r**** **+** **_**d**_. Solve normal equations for _**d**_, then calculate _**b**_**r+1**, which are the parameter values at the next iteration, r+1.
_**X**_ is the matrix of sensitivities calculated at _**b**_**r**.
- Derived independently by Gauss and Newton in the early 1800’s.
- Performed poorly; unstable

**Notes:**
As written, this is the Gauss-Newton method. Both Gauss and Newton independently came up with method of linearizing & iterating to solve nonlinear regression problems.

## Geometry of the Normal Equations

- These normal equations were developed in the early 1800’s. However, in the form presented above, they typically have convergence problems. Modifications are needed to make them work well.

**Notes:**
First term: alters d to point to minimum of linearized obj function.

Took until 1960’s to get equations to work well, to solve problems such as overshoot.

## Making the Normal Equations Work: Scaling

**Scaling** is needed when the parameter values are very different in magnitude, so that sensitivities are also very different. Improves accuracy of the calculated change vector, _**d**_.
**Scaling** does NOT change the magnitude or direction of the parameter change vector, _**d**_.

**Notes:**
Cii = [(XTwX)ii]^(-1/2)

Scaling is conceptually very simple

## Making the Normal Equations Work:The Marquardt Parameter (1950’s)

- Direction change is needed when the scaled parameter change vector, _**d**_, is in a direction that is not likely to help.
- Addition of the **Marquardt term** results in a change vector, _**d**_, that points in a direction that is **closer** to the **steepest descent direction**.

- b2

- b1

**Notes:**
Diagonals of blue term are all 1.0. So, mr is a small number added to each of the diagonals to make them slightly different from 1.0

## Calculating the Marquardt Parameter

- Marquardt parameter used to improve regression performance for ill posed problems. Initially _**m**_**r** =0 for all iterations.
- If _**d**_ is too close to being **orthogonal** to **steepest descent direction**, then Marquardt parameter is used:

- Criteria for implementing Marquardt parameter suggested by Cooley and Naff (1990): If cosine of angle is <0.08  (angle > ~85), then _**m**_**r,new**= *a* x _**m**_**r,old** + *b*
- Cooley and Naff (1990) suggest a=1.5 and b=0.001. In MODFLOW-2000 (PES file) and UCODE_2005, user can specify cosine, a, and b.

**Notes:**
In this example, marquardt parameter used because minimum of linearized obj. func. Is way off to the left.

Blue arrow: direction that d points.
Green arrow: steepest descent direction. 
Use of Marquardt parameter causes blue arrow to become closer to the green arrow.

There have been cases where changing a & b helped, but usually this is not necessary.

## Making the Normal Equations Work: Damping

**Damping** – allow the parameter values to change less than indicated by _**d**_. Damping helps remedy overshoot problems.
- Implemented by inclusion of **damping factor** ******r** in calculation of _**b**_**r+1**: 		_**b**_**r+1**** **=** ********r**** **_**d**_** + **_**b**_**r**.
- Changes the magnitude but not the direction of _**d**_.
- The value of ******r** is calculated internally by MF-2000 or UCODE_2005.
- The value calculated equals the damping required so that no parameter changes by more than user-specified factor MaxChange. A common value is 2.0.

**Notes:**
Obtain d from solution of normal equations.

Pr is a scalar that scales values of d for ALL parameters.

## Making the Normal Equations Work: Damping

- The _**factor**_ by which the regression ‘wants’ to change the jth parameter is:

- where

- Is calculated with ******r** = 1.0.

- If *DMX* is greater than user-specified MaxChange, then ******r** is calculated as:

- Example: Suppose MaxChange = 2.0, and DMX = 10.0:
**	********r** = 2.0 / 10.0 = 0.2
- In this example, each model parameter will be actually changed by only 0.2 times the value of DMX calculated by the regression for that parameter.
- UCODE_2005 allows a different MaxChange values to be defined for each parameter.

**Notes:**
djr: What the regression calculates when solving the normal equations. 

## Stopping Criteria: TolPar and TolSOSC

- Gauss-Newton process is iterative – so, when do you stop iterating?
- Need a **convergence criteria**.
- Two convergence criteria in MODFLOW-2000 and UCODE_2005: **TolPar** and TolSOSC.
**TolPar** (Tolerance fro Parameters): The largest fractional change in a parameter value. For regression to converge:

**TolPar** should ideally be  0.01; larger values may be needed in initial regression runs or for insensitive parameters
- TolSOSC: Convergence criterion for the sum of squared weighted residuals objective function. This criterion stops regression when the model fit isn’t changing much.
- In the final regression runs, the TolPar convergence criterion should be met.

**Notes:**
Small fractional change: means we are at or very close to the minimum of the linearized objective function – Regression ‘sees’ a flat spot.

## MF2KFlow Chart

- Flowchart showing the major steps of the Flow, Observation (OBS), Sensitivity (SEN), and Parameter-Estimation (PES) Processes. (Figure 1 of Hill and others, 2000) .

- MF2KFlow Chart

**Notes:**
Shows how modified Gauss Newton method fits into MF2K flow chart. 
Where is the most computation time in this procedure?

Yellow boxes: assemble needed values for X matrix
Middle green box: solve the normal equations.

## UCODE_2005 Flow Chart

- Flowchart showing the major steps of UCODE_2005. (Figure 1 of Poeter, Hill, Banta, Mehl, and Christensen, 2005) .

## UCODE Flow Chart

- Flowchart showing the major steps of UCODE_2005. (Figure 1 of Poeter, Hill, Banta, Mehl, and Christensen, 2005)
- Lists the input that control each major step.

## Use of Limits on Estimated Parameter Values

- Often used to constrain estimated parameter values to avoid unrealistic values. But unrealistic values can be a valuable indicator of model error!!
- But what about insensitive parameters?
- Applying **limits** results in the estimated parameter being at the edge of reasonable range of values.
- Using **prior information** instead results in parameter values that are in the middle of the range of reasonable parameter values.
- Applying **limits** results in difficulties in propagating uncertainty in limited parameters to uncertainty in predictions.
- Using **prior information** provides a clear framework for propagating uncertainty in the parameters to uncertainty in the predictions.
- Limits are allowed in UCODE_2005, because it is applicable to ANY model, and limits may be needed to maintain parameter values for which solutions can be obtained. See the Parameter_Data input block, documentation p. 70.
- In MODFLOW-2000, this is achieved internally by, for example, not letting hydraulic conductivity parameters be negative unless the user says otherwise.

**Notes:**
If a K or S parameter goes negative, mf2K sets it to a value 2 orders of magnitude < starting value.

## Log-Transformed Parameters

**Log-transforming** parameters can sometimes make a nonlinear regression problem behave more linearly:

- From Carrera and Neuman, WRR, 1986, 22(2), p. 211-227

- Log-transforming also prevents the regression from calculating negative values for the parameters in their native space:

- Log-normal distribution

- Normal distribution

- In MODFLOW-2000 and UCODE, user generally sees values in native space, even when log-transforming is used.  Parameter estimates and statistics are printed in the output files as native and log10 transformed values. (Codes do calculations in natural log space).

## Exercises 5.2a and 5.2b

- DO EXERCISE 5.2a: Define range of reasonable parameter values.
- DO EXERCISE 5.2b: First attempt at regression.
- First, use native parameter values
- Second, try log transforming the K-related parameters

**Notes:**
Reasonable values are in input files. These do not enter into any regression calculations – are just so parameter values can be compared to them in output files.

8b: with native parameter values FAILS – no converge in 10 iters, some K’s go negative.
With log trans, still fails – crashes because an exponent gets too large.

Let’s now look at output.

## Exercise 8a – Reasonable values

- MFI2K screen

## Exercise 8b – Parameter Estimation variables

- MFI2K screen

## Exercise 5.2b – Looking at Results

- Can use GW Chart to look at some results.
- Choose UCODE_2005.
- File  Open File, then choose ex5.2._b. This file contains the parameter estimates for all parameter estimation iterations. Click ‘yes’ in the following dialogue boxes.
- Before proceeding with the exercise, close GW_Chart. MODFLOW-2000 will not run if a needed file is open in GW_Chart.

## Exercise 5.2b – No log transform

## Exercise 5.2b – K’s log transformed

## Analysis of non convergence

- Several K parameters being changed to orders of magnitude smaller than start values.
- What might we do now? Regression did not work with or without log transforming.
- Recall CSS – indicated might not be able to estimate all parameters.
- Should we fix some parameters? Another option: Add prior.
- If so, which parameters?

**Notes:**
Weiss & Smith (GW, 1998) have a nice discussion of which parameters should have prior.  Generally, if a parameter is sensitive, it doesn’t need prior.  If a parameter is insensitive, it does.  Also, errors in sensitive parameters can strongly affect simulated results, while error errors in insensitive parameters have less impact.  Therefore, generally prior should be applied to insensitive parameters. 

## Regression Runs of Exercise 5.2b

- Exercise 5.2b parameter value changes over 10 parameter-estimation iterations when LN=0 for all parameters, and over 6 parameter-estimation iterations when LN=1 for all K parameters.

**Notes:**
Several K parameters being changed to orders of magnitude smaller than start values.
What might we do now? Regression did not work with or without log transforming.

Recall CSS – indicated might not be able to estimate all parameters.
Should we fix some parameters? Another option: Add prior.

## Prior Information

- Allows **direct, independent measurements of model input values** to be included in the regression:

**Observations**	           **Prior**

**Use prior carefully**. Before adding prior, first try performing regression without prior, and assess how much information the dependent observations alone provide towards estimating the model parameters.
- Using prior instead of fixing the value of a parameter allows uncertainty in the parameter to be included in the calculated measures of parameter and prediction uncertainty.
- When model execution times are long, can fix parameter value initially, then use prior information for final regression runs and(or) to evaluate uncertainty

**Notes:**
2nd bullet: this is what we’ve done already in exercise 8b.

Prior information is like a penalty function.  The sum of squared residuals is increased as the parameters depart from our prior information.

## Prior Information

- In MODFLOW-2000 and UCODE, the simulated values related to the prior information are of the form:

- where *b**j* is a parameter value and *a**pj* is a coefficient

- Commonly, the prior information relates to a single parameter value.
- Example: prior information is the field estimate of K from a large-scale aquifer test. *Pp* is the field estimate, and *P’p* is the regression estimate of the K parameter.
- Prior information can also relate to more than one parameter.
- Example: prior information is an annual recharge estimate, but we estimate seasonal recharge in the model. *Pp* is the annual recharge estimate, and *P’p* is the sum of the regression estimates of seasonal recharges.

## Weighting Prior Information

**If the weight reflects the actual uncertainty** in the value specified, that is,

- then this is truly prior information, as used in Bayesian theory.

**If the weight is larger than is consistent with the actual uncertainty**, that is,

- so that the value of STATISTIC is less than that which would accurately reflect the uncertainty, then what is called prior information needs to be classified instead as ***regularization***.

## Weighting Prior Information

- It is sometimes necessary to use regularization to make the problem tractable, but the result is that measures of uncertainty produced by the model will not be correct.
- One approach to this problem is to calibrate the model with the regularization needed to produce a tractable problem, and then change the weighting when calculating prediction uncertainty.

**Notes:**
Can obtain weights on prior in same way as was discussed for observations. Say 95% confident that prior for K of a unit is good to +/- 1 order of magnitude. If pose this way:

1e-5 +/- order of magnitude: 1e-4 to 1e-6
In log10 space: -4 to -6
Std dev in log10 space: [-4 – (-6)] / [2(1.96)] =~ 0.5

## Entering Prior Equations in MODFLOW-2000

- Prior equation entered in Parameter-Estimation Process input file using the format:
- 	**EQNAM PRM "=" ****[SIGN] [COEF "*"] PNAM**** ****[SIGN] [COEF "*"] PNAM**** .........**** ****"STAT" STATP STAT-FLAG PLOT-SYMBOL**
- For the simple example, the prior equation is entered as:
- 	PR_K1   4.0   =   86400 * K1   STAT   0.5e-4   1   4
- In MODFLOW-2000, STATP can relate to the native value or to the log-transformed value of the prior estimate. STAT-FLAG identifies what STATP is (variance, standard deviation, or coefficient of variation) and whether STATP is related to the native or untransformed value.

**Notes:**
Bring up Par Est dialogue in mfi here. 

## Entering Prior Equations in MODFLOW-2000

## Entering Prior Equations in UCODE

- Prior equation in UCODE_2005: Linear_prior_information using the format table:
- In UCODE, if the parameter is log-transformed,
- the value is in native space
- the STATISTIC must relate to the log-transformed value of the prior estimate. STATFLAG identifies what STAT is (variance, standard deviation, or coefficient of variation).
- The equation needs to include the log10 of the parameter.

- BEGIN LINEAR_PRIOR_INFORMATION TABLE
- nrow=2 ncol=5  columnlabels  groupname=prior
- priorname    equation   PriorInfoValue    Statistic  StatFlag
- PRIOR_VK_CB  VK_CB             1.0E-7         0.3        CV
- Eqn_Rch_Ann  0.5*Rch1+0.5*Rch2   37.0         4.0        VAR
- END LINEAR_PRIOR_INFORMATION

## Exercise 5.2c

- DO EXERCISE 5.2c: Assign prior information on parameters, and re-run the regression.
- Do not log transform
- Add prior
- Which parameters should we add prior to?
- Use the starting values as prior value with a coefficient of variation of 0.3
- PROBLEM: Compare ex5.2c.#uout (from 5.2c) and to ex5.2.#uout (from 5.2b)

## Regression Run of Exercise 5.2c

- Exercise 5.2c parameter value changes over the 6 parameter-estimation iterations.

## Prior Information Summary

- Prior information
- Can help stabilize and improve the inversion procedure.
- But use realistic weights when analyzing uncertainties
- Can be thought of as:
- Additional observations, or
- A penalty function
- Is a way for the modeler to incorporate their own judgment about the parameter values.
- Use carefully! (see Weiss & Smith, GW 1998)
- Apply to insensitive parameters, not sensitive parameters.