Least Squares

4.1 Least Squares and the Normal Equations

1) Inconsistent System

Inconsistent System: A system

$$Ax=b$$

that has no exact solution.

Inconsistent systems commonly occur when the number of equations exceeds the number of unknowns or when measured data contains errors.

2) Overdetermined System

Overdetermined System: A system with more equations than unknowns.

For

$$A\in\mathbb{R}^{m\times n},$$

the system is overdetermined when

$$m>n.$$

An overdetermined system generally has no exact solution.

3) Euclidean Norm

Euclidean Norm: The length of a vector.

For

$$r=(r_1,\ldots,r_m)^T,$$

the Euclidean norm, or $2$-norm, is

$$\|r\|_2 = \sqrt{r_1^2+\cdots+r_m^2}.$$

4) Residual

Residual: The difference between the observed vector and the vector predicted by an approximation $x$.

$$r=b-Ax.$$

An exact solution has

$$r=0.$$

5) Least-Squares Solution

Least-Squares Solution: A vector $x$ that minimizes the Euclidean norm of the residual.

$$\min_x\|b-Ax\|_2.$$

Equivalently, it minimizes

$$\|b-Ax\|_2^2 = \sum_{i=1}^m r_i^2.$$

6) Column Space

Column Space: The set of all linear combinations of the columns of $A$.

$$\operatorname{Col}(A) = \{Ax\mid x\in\mathbb{R}^n\}.$$

The least-squares approximation $Ax$ is the point in $\operatorname{Col}(A)$ closest to $b$.

7) Orthogonality Condition

At the least-squares solution, the residual is orthogonal to every column of $A$.

$$A^T(b-Ax)=0.$$

Equivalently,

$$A^Tr=0.$$

Thus, $Ax$ is the orthogonal projection of $b$ onto $\operatorname{Col}(A)$.

8) Normal Equations

Normal Equations: The system obtained from the least-squares orthogonality condition.

$$A^TAx=A^Tb.$$

If the columns of $A$ are linearly independent, then $A^TA$ is nonsingular and the least-squares solution is unique.

Formally,

$$x=(A^TA)^{-1}A^Tb.$$

The inverse should not normally be computed explicitly.

9) Squared Error

Squared Error: The sum of the squared residual components.

$$\operatorname{SE} = r_1^2+\cdots+r_m^2 = \|r\|_2^2.$$

10) Root Mean Squared Error

Root Mean Squared Error: The square root of the average squared residual.

$$\operatorname{RMSE} = \sqrt{ \frac{r_1^2+\cdots+r_m^2}{m} }.$$

Equivalently,

$$\operatorname{RMSE} = \frac{\|r\|_2}{\sqrt{m}}.$$

A smaller RMSE indicates a better fit to the data.

11) Least-Squares Data Fitting

Given data points

$$(t_1,y_1),\ldots,(t_m,y_m),$$

choose a model

$$y=f(t;c_1,\ldots,c_n).$$

If the model is linear in its parameters, substitution of the data produces a matrix problem

$$Ac\approx y.$$

The optimal parameter vector minimizes

$$\|Ac-y\|_2.$$

12) Linear Regression

Linear Regression: Least-squares fitting using the model

$$y=c_0+c_1t.$$

The design matrix is

$$A= \begin{bmatrix} 1&t_1\\ 1&t_2\\ \vdots&\vdots\\ 1&t_m \end{bmatrix},$$

and the parameter vector is

$$c= \begin{bmatrix} c_0\\ c_1 \end{bmatrix}.$$

13) Polynomial Regression

Polynomial Regression: Least-squares fitting using a polynomial model.

For a degree-$d$ model,

$$y = c_0+c_1t+\cdots+c_dt^d.$$

The design matrix is

$$A= \begin{bmatrix} 1&t_1&t_1^2&\cdots&t_1^d\\ 1&t_2&t_2^2&\cdots&t_2^d\\ \vdots&\vdots&\vdots&&\vdots\\ 1&t_m&t_m^2&\cdots&t_m^d \end{bmatrix}.$$

Polynomial regression approximates the data and does not generally pass exactly through every point.

14) Conditioning of the Normal Equations

The condition number of the normal-equation matrix satisfies approximately

$$\operatorname{cond}(A^TA) = \operatorname{cond}(A)^2$$

when the condition number is based on the $2$-norm.

Forming $A^TA$ can therefore square the original conditioning problem.

The normal equations are simple but may be inaccurate for ill-conditioned problems.

4.2 A Survey of Models

1) Model Linear in Parameters

Model Linear in Parameters: A model that can be written as

$$y = c_1\phi_1(t) +\cdots+ c_n\phi_n(t),$$

where the basis functions $\phi_i(t)$ are known and the parameters $c_i$ appear linearly.

The design matrix entries are

$$a_{ij}=\phi_j(t_i).$$

Such a model can be fitted using linear least squares even when the basis functions themselves are nonlinear in $t$.

2) Periodic Model

Periodic Model: A model constructed using sine and cosine functions.

A basic periodic model with period $T$ is

$$y = c_0 + c_1\cos\left(\frac{2\pi t}{T}\right) + c_2\sin\left(\frac{2\pi t}{T}\right).$$

Additional harmonics can be included:

$$y = c_0 + \sum_{k=1}^d \left[ a_k\cos\left(\frac{2\pi kt}{T}\right) + b_k\sin\left(\frac{2\pi kt}{T}\right) \right].$$

This model is nonlinear in $t$ but linear in its parameters.

3) Orthogonal Basis Functions

Orthogonal Basis Functions: Functions whose sampled column vectors have zero dot products.

If the columns $A_i$ and $A_j$ satisfy

$$A_i^TA_j=0 \qquad\text{for }i\neq j,$$

then $A^TA$ is diagonal.

Orthogonal basis functions simplify the least-squares computation and improve numerical conditioning.

4) Data Linearization

Data Linearization: Transforming a nonlinear model into a model that is linear in new parameters.

This allows the transformed problem to be solved using linear least squares.

Linearization changes the error being minimized and does not generally produce the same solution as nonlinear least squares in the original variables.

5) Exponential Model

An exponential model has the form

$$y=c_1e^{c_2t}.$$

Taking logarithms gives

$$\ln y = \ln c_1+c_2t.$$

Define

$$z=\ln y \qquad\text{and}\qquad a=\ln c_1.$$

The transformed model is

$$z=a+c_2t.$$

After fitting $a$ and $c_2$,

$$c_1=e^a.$$

This transformation requires

$$y_i>0.$$

6) Power Model

A power model has the form

$$y=c_1t^{c_2}.$$

Taking logarithms gives

$$\ln y = \ln c_1+c_2\ln t.$$

Define

$$z=\ln y, \qquad s=\ln t, \qquad a=\ln c_1.$$

Then

$$z=a+c_2s.$$

This transformation requires positive $t_i$ and $y_i$.

7) Original-Space Error

For an exponential model, the original least-squares objective is

$$\sum_{i=1}^m \left( c_1e^{c_2t_i}-y_i \right)^2.$$

8) Transformed-Space Error

The log-linearized objective is

$$\sum_{i=1}^m \left( \ln c_1+c_2t_i-\ln y_i \right)^2.$$

These objectives are different and generally produce different parameter values.

The appropriate model depends on whether absolute errors or relative/logarithmic errors are more meaningful.

4.3 QR Factorization

1) Dot Product

For vectors $u,v\in\mathbb{R}^m$, the dot product is

$$u^Tv = \sum_{i=1}^m u_iv_i.$$

2) Orthogonal Vectors

Orthogonal Vectors: Vectors whose dot product is zero.

$$u^Tv=0.$$

3) Orthonormal Vectors

Orthonormal Vectors: Unit vectors that are mutually orthogonal.

$$q_i^Tq_j = \begin{cases} 1, & i=j,\\ 0, & i\neq j. \end{cases}$$

4) Orthogonal Matrix

Orthogonal Matrix: A square matrix $Q$ satisfying

$$Q^TQ=QQ^T=I.$$

Therefore,

$$Q^{-1}=Q^T.$$

Orthogonal matrices preserve Euclidean norms:

$$\|Qx\|_2=\|x\|_2.$$

They also preserve dot products and angles.

5) Reduced QR Factorization

For

$$A\in\mathbb{R}^{m\times n}, \qquad m\geq n,$$

with linearly independent columns, the reduced QR factorization is

$$A=QR,$$

where:

• $Q\in\mathbb{R}^{m\times n}$ has orthonormal columns.
• $R\in\mathbb{R}^{n\times n}$ is upper triangular.

Thus,

$$Q^TQ=I_n.$$

6) Full QR Factorization

The full QR factorization is

$$A=QR,$$

where:

• $Q\in\mathbb{R}^{m\times m}$ is orthogonal.
• $R\in\mathbb{R}^{m\times n}$ is upper trapezoidal.

The first $n$ columns of $Q$ span $\operatorname{Col}(A)$.

7) Projection onto a Vector

The projection of $a$ onto a unit vector $q$ is

$$\operatorname{proj}_q(a) = (q^Ta)q.$$

For a nonunit vector $u$,

$$\operatorname{proj}_u(a) = \frac{u^Ta}{u^Tu}u.$$

8) Classical Gram-Schmidt Orthogonalization

Classical Gram-Schmidt: A method that converts linearly independent vectors into orthonormal vectors.

Given the columns

$$A_1,\ldots,A_n,$$

calculate

$$y_j = A_j - \sum_{i=1}^{j-1} (q_i^TA_j)q_i$$

and normalize:

$$q_j = \frac{y_j}{\|y_j\|_2}.$$

The entries of $R$ are

$$r_{ij}=q_i^TA_j \qquad\text{for }i<j$$

and

$$r_{jj}=\|y_j\|_2.$$

9) Modified Gram-Schmidt

Modified Gram-Schmidt: A numerically improved version of Gram-Schmidt that subtracts projections one at a time.

Start with

$$y=A_j.$$

Then, for $i=1,\ldots,j-1$,

$$r_{ij}=q_i^Ty$$

and

$$y\leftarrow y-r_{ij}q_i.$$

Finally,

$$r_{jj}=\|y\|_2$$

and

$$q_j=\frac{y}{r_{jj}}.$$

Modified Gram-Schmidt is mathematically equivalent to classical Gram-Schmidt in exact arithmetic but usually more accurate in floating-point arithmetic.

10) Least Squares by Reduced QR

Given

$$A=QR,$$

the least-squares problem becomes

$$\min_x\|QRx-b\|_2.$$

Because $Q$ has orthonormal columns, the solution satisfies

$$Rx=Q^Tb.$$

Since $R$ is upper triangular, $x$ is found using back substitution.

This method avoids forming $A^TA$.

11) Least Squares by Full QR

For the full factorization, partition

$$Q= \begin{bmatrix} Q_1&Q_2 \end{bmatrix}$$

and

$$R= \begin{bmatrix} R_1\\ 0 \end{bmatrix}.$$

Then

$$\|Ax-b\|_2^2 = \|R_1x-Q_1^Tb\|_2^2 + \|Q_2^Tb\|_2^2.$$

The least-squares solution satisfies

$$R_1x=Q_1^Tb.$$

The minimum residual norm is

$$\|r\|_2 = \|Q_2^Tb\|_2.$$

12) Householder Reflector

Householder Reflector: An orthogonal matrix that reflects vectors across a hyperplane.

For a nonzero vector $v$,

$$H = I - 2\frac{vv^T}{v^Tv}.$$

It satisfies

$$H^T=H,$$ $$H^{-1}=H,$$

and

$$H^TH=I.$$

13) Constructing a Householder Reflector

To transform a vector $x$ into a multiple of the first coordinate vector, choose

$$\alpha = -\operatorname{sign}(x_1)\|x\|_2$$

and

$$v=x-\alpha e_1.$$

Then

$$H = I-2\frac{vv^T}{v^Tv}$$

satisfies

$$Hx=\alpha e_1.$$

The sign is chosen to reduce cancellation.

14) Householder QR Factorization

Householder reflectors successively introduce zeros below the diagonal:

$$H_k\cdots H_2H_1A=R.$$

Because each reflector is orthogonal and symmetric,

$$A=H_1H_2\cdots H_kR.$$

Thus,

$$Q=H_1H_2\cdots H_k.$$

Householder QR is generally more stable and efficient than Gram-Schmidt for dense matrices.

15) QR versus Normal Equations

Method	Main Equation	Conditioning	Stability
Normal equations	$A^TAx=A^Tb$	Squares $\operatorname{cond}(A)$	Lower
QR factorization	$Rx=Q^Tb$	Uses $A$ directly	Higher

QR factorization is preferred when numerical accuracy matters.

4.4 Generalized Minimum Residual Method

1) Krylov Subspace

Krylov Subspace: The vector space generated by repeated multiplication of a residual by $A$.

Given

$$r_0=b-Ax_0,$$

the order-$k$ Krylov subspace is

$$\mathcal{K}_k(A,r_0) = \operatorname{span} \left\{ r_0,Ar_0,\ldots,A^{k-1}r_0 \right\}.$$

2) GMRES

Generalized Minimum Residual Method: An iterative method for solving large, sparse, nonsymmetric systems.

At step $k$, GMRES selects

$$x_k\in x_0+\mathcal{K}_k(A,r_0)$$

that minimizes

$$\|b-Ax_k\|_2.$$

3) Arnoldi Process

Arnoldi Process: An orthogonalization process that constructs an orthonormal basis of a Krylov subspace.

Starting with

$$q_1=\frac{r_0}{\|r_0\|_2},$$

calculate

$$w=Aq_k.$$

For $j=1,\ldots,k$,

$$h_{jk}=q_j^Tw$$

and

$$w\leftarrow w-h_{jk}q_j.$$

Then,

$$h_{k+1,k}=\|w\|_2$$

and

$$q_{k+1} = \frac{w}{h_{k+1,k}}.$$

4) Arnoldi Relation

Let

$$Q_k= \begin{bmatrix} q_1&\cdots&q_k \end{bmatrix}.$$

The Arnoldi process produces

$$AQ_k = Q_{k+1}H_k,$$

where $H_k$ is an upper-Hessenberg matrix.

5) Upper-Hessenberg Matrix

Upper-Hessenberg Matrix: A matrix whose entries below the first subdiagonal are zero.

$$h_{ij}=0 \qquad\text{when }i>j+1.$$

6) GMRES Least-Squares Problem

Write the correction as

$$x_k=x_0+Q_ky.$$

GMRES determines $y$ by minimizing

$$\left\| H_ky-\|r_0\|_2e_1 \right\|_2.$$

The updated approximation is

$$x_k=x_0+Q_ky.$$

The residual norm decreases monotonically as the Krylov subspace grows.

7) Restarted GMRES

Restarted GMRES: A memory-limited version of GMRES that restarts after a fixed number of iterations.

After $m$ iterations:

1. Keep the current approximation $x_m$.
2. Discard the existing Krylov basis.
3. Restart GMRES using $x_m$ as the new initial guess.

Restarting reduces memory use but may slow or prevent convergence.

8) Preconditioned GMRES

Choose a preconditioner $M$ and solve

$$M^{-1}Ax=M^{-1}b.$$

The initial residual becomes

$$r_0=M^{-1}(b-Ax_0).$$

The Arnoldi multiplication becomes

$$w=M^{-1}Aq_k.$$

The inverse $M^{-1}$ should not be formed explicitly; instead, systems involving $M$ should be solved.

4.5 Nonlinear Least Squares

1) Nonlinear Residual Function

Let

$$r(x) = \begin{bmatrix} r_1(x)\\ \vdots\\ r_m(x) \end{bmatrix}$$

be a nonlinear residual vector.

The nonlinear least-squares problem is

$$\min_x \frac{1}{2}\|r(x)\|_2^2.$$

2) Least-Squares Objective Function

Define

$$E(x) = \frac{1}{2} r(x)^Tr(x) = \frac{1}{2} \sum_{i=1}^m r_i(x)^2.$$

The factor $1/2$ simplifies derivative formulas.

3) Residual Jacobian

Residual Jacobian: The matrix of partial derivatives of the residual functions.

$$J(x) = Dr(x) = \begin{bmatrix} \dfrac{\partial r_1}{\partial x_1} & \cdots & \dfrac{\partial r_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \dfrac{\partial r_m}{\partial x_1} & \cdots & \dfrac{\partial r_m}{\partial x_n} \end{bmatrix}.$$

4) Gradient of the Squared Error

The gradient of the objective is

$$\nabla E(x) = J(x)^Tr(x).$$

A stationary point satisfies

$$J(x)^Tr(x)=0.$$

5) Exact Hessian

The Hessian of $E$ is

$$\nabla^2E(x) = J(x)^TJ(x) + \sum_{i=1}^m r_i(x)\nabla^2r_i(x).$$

Gauss-Newton neglects the second term.

6) Gauss-Newton Method

Gauss-Newton Method: An iterative method for nonlinear least squares based on a local linear approximation of the residual.

At iteration $k$, set

$$J_k=J(x_k).$$

Solve

$$J_k^TJ_kv_k = -J_k^Tr(x_k)$$

and update

$$x_{k+1} = x_k+v_k.$$

Equivalently, solve the linear least-squares problem

$$\min_v \|J_kv+r(x_k)\|_2.$$

Using QR to solve this subproblem is generally more stable than forming $J_k^TJ_k$.

7) Gauss-Newton Convergence

Gauss-Newton works best when:

• The initial guess is close to a solution.
• The Jacobian has full column rank.
• The residual at the solution is small.
• The model is not severely nonlinear near the solution.

It may converge to a local minimum, a maximum, or another stationary point.

8) Nonlinear Parameter Model

Given data points

$$(t_i,y_i), \qquad i=1,\ldots,m,$$

and a nonlinear model

$$y=f(t;c),$$

define the residuals

$$r_i(c) = f(t_i;c)-y_i.$$

The parameter vector is determined by minimizing

$$\sum_{i=1}^m r_i(c)^2.$$

The residual Jacobian entries are

$$J_{ij} = \frac{\partial r_i}{\partial c_j} = \frac{\partial f(t_i;c)}{\partial c_j}.$$

9) Regularization

Regularization: The addition of a stabilizing term to an ill-conditioned problem.

Regularization reduces sensitivity and controls excessively large parameter changes.

10) Levenberg-Marquardt Method

Levenberg-Marquardt Method: A regularized Gauss-Newton method for nonlinear least squares.

At iteration $k$, let

$$J_k=J(x_k).$$

Solve

$$\left( J_k^TJ_k + \lambda\operatorname{diag}(J_k^TJ_k) \right)v_k = -J_k^Tr(x_k)$$

and update

$$x_{k+1} = x_k+v_k.$$

The parameter $\lambda\geq0$ controls the amount of regularization.

11) Effect of the Damping Parameter

When

$$\lambda=0,$$

Levenberg-Marquardt becomes the Gauss-Newton Method.

For a large $\lambda$, the step behaves more like a scaled steepest-descent step.

Increasing $\lambda$ generally:

• Improves conditioning.
• Produces smaller, more conservative steps.
• Expands the set of usable initial guesses.

Decreasing $\lambda$ generally:

• Produces faster Gauss-Newton-like convergence near a solution.
• Reduces the effect of regularization.

12) Adaptive Damping

A common adaptive strategy is:

• Accept the step and decrease $\lambda$ when the residual decreases.
• Reject the step and increase $\lambda$ when the residual increases.

This balances convergence speed and robustness.

Comparison of Least-Squares Methods

Method	Problem Type	Main Equation	Main Property
Normal Equations	Linear	$A^TAx=A^Tb$	Simple but squares conditioning
QR Factorization	Linear	$Rx=Q^Tb$	More stable
GMRES	Large sparse linear system	Krylov least-squares problem	Handles nonsymmetric matrices
Gauss-Newton	Nonlinear least squares	$J^TJv=-J^Tr$	Fast near a good solution
Levenberg-Marquardt	Nonlinear least squares	Regularized Gauss-Newton system	More robust under poor conditioning

Essential Concepts

1. Inconsistent System: A system with no exact solution.
2. Overdetermined System: Has more equations than unknowns.
3. Residual: The vector $r=b-Ax$.
4. Least-Squares Solution: Minimizes $\|b-Ax\|_2$.
5. Orthogonality Condition: The residual satisfies $A^Tr=0$.
6. Normal Equations: The system $A^TAx=A^Tb$.
7. Squared Error: The value $\|r\|_2^2$.
8. RMSE: The residual norm divided by $\sqrt{m}$.
9. Design Matrix: Contains the model basis functions evaluated at the data points.
10. Linear Regression: Fits a straight line using least squares.
11. Polynomial Regression: Fits a polynomial without requiring exact interpolation.
12. Data Linearization: Transforms a nonlinear model into a linear parameter model.
13. Orthogonal Matrix: Satisfies $Q^{-1}=Q^T$.
14. QR Factorization: Expresses $A$ as an orthogonal factor times an upper-triangular factor.
15. Gram-Schmidt: Constructs orthonormal vectors from linearly independent vectors.
16. Modified Gram-Schmidt: Improves floating-point accuracy.
17. Householder Reflector: Stably introduces zeros using an orthogonal reflection.
18. Krylov Subspace: Spans repeated matrix applications to a residual.
19. GMRES: Minimizes the residual over an expanding Krylov subspace.
20. Residual Jacobian: Contains derivatives of nonlinear residual functions.
21. Gauss-Newton Method: Applies linear least squares to a local residual approximation.
22. Regularization: Stabilizes an ill-conditioned minimization problem.
23. Levenberg-Marquardt Method: Combines Gauss-Newton with damping.