Linear and Weighted Least Squares

4.3.2 Linear and Weighted Least Squares

Once an estimation problem has been formulated, the next question is how to determine the unknown state or parameters from the available observations.

The simplest important case is the linear estimation problem:

\[ \boxed{ \mathbf z=\mathbf H\mathbf x+\mathbf v } \]

where

• \mathbf z is the measurement or observation vector,
• \mathbf x is the vector of unknown states or parameters,
• \mathbf H is the design, observation, or sensitivity matrix,
• \mathbf v is the measurement error.

If the measurements were perfect and the equations exactly determined, the problem could sometimes be solved directly. Real estimation problems, however, usually contain more measurements than unknowns, and those measurements contain noise.

Therefore, in general,

\[ \boxed{ \mathbf z \neq \mathbf H\mathbf x } \]

for every individual observation simultaneously.

Instead of searching for a state that satisfies every noisy equation exactly, least squares searches for the state that makes the overall disagreement between the observations and model predictions as small as possible.

Conceptually,

\[ \boxed{ \text{Measurements} + \text{Measurement Model} \rightarrow \text{Residuals} \rightarrow \text{Residual Minimization} \rightarrow \text{Best-Fit Estimate} } \]

Linear least squares is therefore one of the fundamental mathematical foundations of state estimation, parameter estimation, sensor calibration, trajectory reconstruction, orbit determination, and many other estimation problems.

4.3.2.1 From the Estimation Problem to Least Squares

From Section 4.3.1, the general nonlinear observation model was

\[ \mathbf z=h(\mathbf x,\boldsymbol{\theta})+\mathbf v. \]

For the present section, consider the linear case:

\[ \boxed{ \mathbf z=\mathbf H\mathbf x+\mathbf v } \]

Suppose there are m measurements and n unknown quantities:

\[ \mathbf z\in\mathbb R^m, \qquad \mathbf x\in\mathbb R^n, \qquad \mathbf H\in\mathbb R^{m\times n}. \]

In many practical estimation problems,

\[ m>n. \]

The system is then overdetermined.

Because the measurements contain noise, a value of \mathbf x that satisfies every equation exactly generally does not exist.

Therefore, the problem changes from

\[ \text{Solve exactly} \]

to

\[ \boxed{ \text{Find the state that best fits all observations simultaneously.} } \]

This is the central idea of least-squares estimation.

4.3.2.2 Linear Observation Model

Consider the observation equations

\[ \begin{aligned} z_1 &= H_{11}x_1+H_{12}x_2+\cdots+H_{1n}x_n+v_1,\\ z_2 &= H_{21}x_1+H_{22}x_2+\cdots+H_{2n}x_n+v_2,\\ &\vdots\\ z_m &= H_{m1}x_1+H_{m2}x_2+\cdots+H_{mn}x_n+v_m. \end{aligned} \]

These equations can be collected into matrix form:

\[ \boxed{ \mathbf z=\mathbf H\mathbf x+\mathbf v. } \]

The matrix

\[ \mathbf H= \begin{bmatrix} H_{11}&H_{12}&\cdots&H_{1n}\\ H_{21}&H_{22}&\cdots&H_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ H_{m1}&H_{m2}&\cdots&H_{mn} \end{bmatrix} \]

contains the relationship between the unknown quantities and observations.

Each row corresponds to an observation equation.

Each column describes how the complete observation set responds to one unknown quantity.

Thus,

\[ \boxed{ \text{Rows of }\mathbf H \rightarrow \text{Measurements} } \]

and

\[ \boxed{ \text{Columns of }\mathbf H \rightarrow \text{Unknown State or Parameter Components}. } \]

This interpretation becomes extremely important when studying rank, geometry, conditioning, observability, and parameter correlation.

4.3.2.3 Predicted Measurements and Residuals

Given a candidate state estimate \hat{\mathbf x}, the corresponding predicted observations are

\[ \boxed{ \hat{\mathbf z}=\mathbf H\hat{\mathbf x}. } \]

The residual vector is

\[ \boxed{ \mathbf r = \mathbf z-\hat{\mathbf z} = \mathbf z-\mathbf H\hat{\mathbf x}. } \]

Component-wise,

\[ r_i=z_i-\hat z_i. \]

Therefore,

\[ \boxed{ \text{Residual} = \text{Observed Measurement} - \text{Predicted Measurement}. } \]

If the candidate state changes, the predicted observations change, and therefore the residuals also change.

Least squares chooses the candidate state for which the overall residual magnitude is minimized.

Remember that the residual is a measurement-space quantity. It should not be confused with the state-estimation error.

4.3.2.4 Why an Exact Solution Usually Does Not Exist

Suppose three measurements are used to estimate one unknown quantity:

\[ \begin{aligned} z_1 &= x+v_1,\\ z_2 &= x+v_2,\\ z_3 &= x+v_3. \end{aligned} \]

Suppose the measured values are

\[ z_1=10.1,\qquad z_2=9.9,\qquad z_3=10.2. \]

There is no single value of x that satisfies

\[ x=10.1,\qquad x=9.9,\qquad x=10.2 \]

simultaneously.

The inconsistency is not necessarily a problem with the model. It is the natural consequence of measurement uncertainty.

Instead, we seek the value of x that gives the smallest overall disagreement.

Thus,

\[ \boxed{ \text{Noisy Overdetermined System} \rightarrow \text{Best-Fit Solution Rather Than Exact Solution}. } \]

4.3.2.5 Ordinary Least-Squares Cost Function

For the linear model,

\[ \mathbf r=\mathbf z-\mathbf H\mathbf x. \]

Ordinary least squares defines the cost function

\[ \boxed{ J(\mathbf x)=\mathbf r^T\mathbf r. } \]

Substituting the residual,

\[ \boxed{ J(\mathbf x) = (\mathbf z-\mathbf H\mathbf x)^T (\mathbf z-\mathbf H\mathbf x). } \]

Since

\[ \mathbf r^T\mathbf r = \sum_{i=1}^{m}r_i^2, \]

we may also write

\[ \boxed{ J(\mathbf x) = \sum_{i=1}^{m} (z_i-\hat z_i)^2. } \]

The least-squares estimate is

\[ \boxed{ \hat{\mathbf x}_{LS} = \arg\min_{\mathbf x} J(\mathbf x). } \]

This is why the method is called least squares: it minimizes the sum of the squared measurement residuals.

4.3.2.6 Why Square the Residuals?

Suppose the residuals are

\[ r_1=+2, \qquad r_2=-2. \]

If they were simply added,

\[ r_1+r_2=0, \]

which could misleadingly suggest perfect agreement.

Squaring gives

\[ r_1^2+r_2^2=4+4=8. \]

Therefore positive and negative errors cannot cancel.

Squaring also produces a smooth differentiable objective function:

\[ J(\mathbf x)=\mathbf r^T\mathbf r\ge0. \]

This allows the minimum to be found analytically for the linear case.

Thus,

\[ \boxed{ \text{Residual Squaring} \rightarrow \text{No Sign Cancellation} + \text{Nonnegative Cost} + \text{Convenient Optimization}. } \]

4.3.2.7 Derivation of the Ordinary Least-Squares Solution

Start with

\[ J(\mathbf x) = (\mathbf z-\mathbf H\mathbf x)^T (\mathbf z-\mathbf H\mathbf x). \]

Expand:

\[ J = \mathbf z^T\mathbf z - 2\mathbf x^T\mathbf H^T\mathbf z + \mathbf x^T\mathbf H^T\mathbf H\mathbf x. \]

At the minimum,

\[ \frac{\partial J}{\partial\mathbf x}=0. \]

Therefore,

\[ -2\mathbf H^T\mathbf z + 2\mathbf H^T\mathbf H\mathbf x = 0. \]

Hence,

\[ \boxed{ \mathbf H^T\mathbf H\hat{\mathbf x} = \mathbf H^T\mathbf z. } \]

These are the normal equations.

If \mathbf H^T\mathbf H is nonsingular,

\[ \boxed{ \hat{\mathbf x}_{LS} = (\mathbf H^T\mathbf H)^{-1} \mathbf H^T\mathbf z. } \]

This is the classical ordinary least-squares solution.

4.3.2.8 Meaning of the Normal Equations

The normal equations are

\[ \mathbf H^T (\mathbf z-\mathbf H\hat{\mathbf x}) = 0. \]

Since

\[ \hat{\mathbf r} = \mathbf z-\mathbf H\hat{\mathbf x}, \]

we obtain

\[ \boxed{ \mathbf H^T\hat{\mathbf r}=0. } \]

This has an important geometric interpretation.

The residual vector at the least-squares solution is orthogonal to every column of \mathbf H.

Therefore,

\[ \boxed{ \hat{\mathbf r} \perp \operatorname{Col}(\mathbf H). } \]

Least squares therefore does more than simply produce a formula. It finds the point in the measurement-model subspace that is closest to the actual observation vector.

Engineering interpretation: The orthogonality condition means there is no remaining component of the residual that can be reduced by moving the estimate within the column space of the design matrix. The remaining residual lies outside the set of measurements that the assumed linear model can reproduce.

4.3.2.9 Geometric Interpretation of Least Squares

The model

\[ \mathbf H\mathbf x \]

can produce only vectors lying in the column space of \mathbf H:

\[ \operatorname{Col}(\mathbf H). \]

The measured vector \mathbf z, because of noise, may not lie exactly in that space.

Least squares projects \mathbf z onto the column space:

\[ \boxed{ \mathbf z \rightarrow \hat{\mathbf z} = \mathbf H\hat{\mathbf x}. } \]

The residual is

\[ \hat{\mathbf r} = \mathbf z-\hat{\mathbf z}. \]

At the optimum,

\[ \boxed{ \mathbf H^T\hat{\mathbf r}=0. } \]

Conceptually,

\[ \boxed{ \text{Observation Vector} = \text{Best Model-Compatible Component} + \text{Orthogonal Residual}. } \]

This geometric viewpoint is especially useful for understanding rank deficiency and poorly observable directions.

4.3.2.10 Projection Matrix

For full-column-rank \mathbf H,

\[ \hat{\mathbf z} = \mathbf H (\mathbf H^T\mathbf H)^{-1} \mathbf H^T\mathbf z. \]

Define

\[ \boxed{ \mathbf P_H = \mathbf H (\mathbf H^T\mathbf H)^{-1} \mathbf H^T. } \]

Then

\[ \boxed{ \hat{\mathbf z} = \mathbf P_H\mathbf z. } \]

The residual becomes

\[ \boxed{ \hat{\mathbf r} = (\mathbf I-\mathbf P_H)\mathbf z. } \]

Thus \mathbf P_H projects the measurement vector onto the column space of \mathbf H.

This gives a direct mathematical representation of the geometric interpretation of least squares.

4.3.2.11 Simple Scalar Least-Squares Example

Consider again

\[ z_i=x+v_i \]

with

\[ z_1=10.1,\qquad z_2=9.9,\qquad z_3=10.2. \]

The cost function is

\[ J(x) = (10.1-x)^2 + (9.9-x)^2 + (10.2-x)^2. \]

Differentiate:

\[ \frac{dJ}{dx} = -2(10.1-x) -2(9.9-x) -2(10.2-x). \]

Setting this equal to zero gives

\[ 3x=30.2. \]

Therefore,

\[ \boxed{ \hat x=10.0667. } \]

For equally accurate independent measurements,

\[ \boxed{ \hat x = \frac{1}{m} \sum_{i=1}^{m}z_i. } \]

The arithmetic mean is therefore a special case of least-squares estimation.

4.3.2.12 Why Ordinary Least Squares Is Not Always Appropriate

Ordinary least squares treats every residual equally:

\[ J=\sum_i r_i^2. \]

But suppose two sensors measure the same quantity and have standard deviations

\[ \sigma_1=0.1 \]

and

\[ \sigma_2=10. \]

The first measurement is much more precise.

Treating both measurements equally ignores known information about their quality.

This motivates weighted least squares.

The principle is

\[ \boxed{ \text{Measurements Should Influence the Estimate According to Their Uncertainty}. } \]

4.3.2.13 Measurement Covariance

Let

\[ \mathbf v \]

represent measurement error.

Its covariance is

\[ \boxed{ \mathbf R = E[ \mathbf v\mathbf v^T ]. } \]

If the measurement errors are independent,

\[ \mathbf R = \begin{bmatrix} \sigma_1^2&0&\cdots&0\\ 0&\sigma_2^2&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&\sigma_m^2 \end{bmatrix}. \]

Smaller variance means greater measurement precision.

Therefore,

\[ \boxed{ \sigma_i^2\downarrow \Rightarrow \text{Measurement Confidence}\uparrow. } \]

The covariance matrix allows this information to be incorporated mathematically into the estimator.

4.3.2.14 Weight Matrix

The conventional statistical weighting matrix is

\[ \boxed{ \mathbf W=\mathbf R^{-1}. } \]

For independent measurements,

\[ \mathbf W = \begin{bmatrix} 1/\sigma_1^2&0&\cdots\\ 0&1/\sigma_2^2&\cdots\\ \vdots&\vdots&\ddots \end{bmatrix}. \]

Therefore,

\[ \boxed{ w_i=\frac{1}{\sigma_i^2}. } \]

A measurement with smaller uncertainty receives greater weight.

A measurement with larger uncertainty receives less weight.

Thus,

\[ \boxed{ \text{Measurement Accuracy} \rightarrow \text{Statistical Weight} \rightarrow \text{Influence on Estimate}. } \]

4.3.2.15 Weighted Least-Squares Cost Function

Weighted least squares defines

\[ \boxed{ J(\mathbf x) = \mathbf r^T\mathbf W\mathbf r. } \]

Using

\[ \mathbf W=\mathbf R^{-1}, \]

we obtain

\[ \boxed{ J(\mathbf x) = (\mathbf z-\mathbf H\mathbf x)^T \mathbf R^{-1} (\mathbf z-\mathbf H\mathbf x). } \]

For independent measurements,

\[ \boxed{ J = \sum_{i=1}^{m} \frac{r_i^2}{\sigma_i^2}. } \]

Each residual is therefore evaluated relative to the expected uncertainty of its measurement.

This is more meaningful than comparing raw residual magnitudes when the measurements have different accuracy or units.

4.3.2.16 Normalized Residual Interpretation

For an independent scalar measurement,

\[ \frac{r_i}{\sigma_i} \]

represents the residual relative to its standard deviation.

The weighted cost becomes

\[ J = \sum_i \left( \frac{r_i}{\sigma_i} \right)^2. \]

Therefore weighted least squares effectively asks:

How large is each measurement mismatch relative to the uncertainty expected for that measurement?

A residual of 1 m may be extremely significant for a millimetre-accuracy sensor but almost insignificant for a measurement with tens of metres of uncertainty.

Thus,

\[ \boxed{ \text{Residual Magnitude Alone Does Not Determine Its Statistical Importance}. } \]

4.3.2.17 Derivation of Weighted Least Squares

Start from

\[ J(\mathbf x) = (\mathbf z-\mathbf H\mathbf x)^T \mathbf W (\mathbf z-\mathbf H\mathbf x). \]

Assuming \mathbf W is symmetric,

\[ \frac{\partial J}{\partial\mathbf x} = -2\mathbf H^T\mathbf W (\mathbf z-\mathbf H\mathbf x). \]

At the minimum,

\[ \mathbf H^T\mathbf W (\mathbf z-\mathbf H\hat{\mathbf x}) = 0. \]

Therefore,

\[ \boxed{ \mathbf H^T\mathbf W\mathbf H \hat{\mathbf x} = \mathbf H^T\mathbf W\mathbf z. } \]

These are the weighted normal equations.

If the matrix is nonsingular,

\[ \boxed{ \hat{\mathbf x}_{WLS} = (\mathbf H^T\mathbf W\mathbf H)^{-1} \mathbf H^T\mathbf W\mathbf z. } \]

Using

\[ \mathbf W=\mathbf R^{-1}, \]

the solution becomes

\[ \boxed{ \hat{\mathbf x}_{WLS} = (\mathbf H^T\mathbf R^{-1}\mathbf H)^{-1} \mathbf H^T\mathbf R^{-1}\mathbf z. } \]

This is one of the most important equations in estimation theory.

4.3.2.18 Information Matrix

The matrix

\[ \boxed{ \boldsymbol{\Lambda} = \mathbf H^T\mathbf R^{-1}\mathbf H } \]

is often interpreted as an information matrix.

It describes how much measurement information is available about the unknown state.

The contribution of measurement i can be understood through its sensitivity and uncertainty.

A measurement provides more useful information when it has

• strong sensitivity to the unknown quantity,
• low measurement uncertainty,
• useful geometry,
• independent information relative to other observations.

Conceptually,

\[ \boxed{ \text{Information} \sim \frac{\text{Sensitivity}^2} {\text{Measurement Uncertainty}}. } \]

This relationship helps connect least squares with covariance, observability, experiment design, and recursive estimation.

Engineering interpretation: The information matrix is also useful for diagnosis. Its eigenvectors identify combinations of states that are strongly or weakly constrained, while small eigenvalues indicate weak-information directions that typically correspond to large uncertainty after inversion.

4.3.2.19 Covariance of the Least-Squares Estimate

Suppose

\[ \mathbf z = \mathbf H\mathbf x+\mathbf v \]

with

\[ E[\mathbf v]=0 \]

and

\[ E[\mathbf v\mathbf v^T]=\mathbf R. \]

The weighted least-squares estimate is

\[ \hat{\mathbf x} = (\mathbf H^T\mathbf R^{-1}\mathbf H)^{-1} \mathbf H^T\mathbf R^{-1}\mathbf z. \]

Under the standard linear model assumptions, the covariance of the estimation error is

\[ \boxed{ \mathbf P_{\hat x} = (\mathbf H^T\mathbf R^{-1}\mathbf H)^{-1}. } \]

Therefore,

\[ \boxed{ \text{Information Matrix} = \mathbf P_{\hat x}^{-1}. } \]

or

\[ \boxed{ \mathbf P_{\hat x}^{-1} = \mathbf H^T\mathbf R^{-1}\mathbf H. } \]

This is a powerful interpretation:

\[ \boxed{ \text{More Independent Measurement Information} \rightarrow \text{Smaller Estimate Uncertainty}. } \]

4.3.2.20 Meaning of Estimate Covariance

Suppose

\[ \mathbf P_{\hat x} = \begin{bmatrix} \sigma_{x_1}^2 & P_{12}\\ P_{21} & \sigma_{x_2}^2 \end{bmatrix}. \]

The diagonal terms represent estimated variances:

\[ P_{ii}=\sigma_{x_i}^2. \]

The off-diagonal terms describe covariance between estimation errors.

The correlation coefficient is

\[ \boxed{ \rho_{ij} = \frac{P_{ij}} {\sqrt{P_{ii}P_{jj}}}. } \]

If

\[ |\rho_{ij}|\approx1, \]

the estimated quantities are strongly correlated.

This may indicate that the observations have difficulty distinguishing their individual effects.

Therefore covariance provides information not only about how uncertain the solution is, but also about which unknowns are coupled.

Engineering interpretation: For multidimensional states, covariance also has a geometric interpretation. A position covariance defines an uncertainty ellipsoid whose principal axes follow the covariance eigenvectors, so a small scalar residual can coexist with a large uncertainty direction.

4.3.2.21 Unbiasedness of Linear Least Squares

For

\[ \mathbf z=\mathbf H\mathbf x+\mathbf v, \]

suppose

\[ E[\mathbf v]=0. \]

Then

\[ E[\hat{\mathbf x}] = \mathbf x \]

for the standard full-rank least-squares estimator.

Therefore,

\[ \boxed{ E[\hat{\mathbf x}-\mathbf x]=0. } \]

The estimator is unbiased under the assumed model.

However, this conclusion depends on the model assumptions.

If the measurement contains an unmodelled bias,

\[ \mathbf z = \mathbf H\mathbf x+\mathbf b+\boldsymbol{\eta}, \]

then the estimated state can also become biased.

Thus,

\[ \boxed{ \text{Unbiased Noise Assumption} \neq \text{Guarantee of an Unbiased Real System}. } \]

Model errors and systematic sensor errors still matter.

4.3.2.22 Gauss-Markov Interpretation

Under the classical assumptions that

• the model is linear,
• the observation errors have zero mean,
• the covariance is correctly represented,
• the design matrix has sufficient rank,

weighted least squares has an important optimality property.

Among linear unbiased estimators, it provides a minimum-variance solution under the corresponding covariance model.

This is closely related to the Gauss-Markov theorem.

Thus least squares is not simply a convenient curve-fitting technique. Under appropriate assumptions, it has a rigorous statistical interpretation.

4.3.2.23 Weighted Least Squares and Gaussian Maximum Likelihood

Suppose

\[ \mathbf v \sim \mathcal N(\mathbf0,\mathbf R). \]

Then

\[ p(\mathbf z|\mathbf x) \propto \exp \left[ -\frac12 (\mathbf z-\mathbf H\mathbf x)^T \mathbf R^{-1} (\mathbf z-\mathbf H\mathbf x) \right]. \]

Maximizing the likelihood is equivalent to minimizing

\[ (\mathbf z-\mathbf H\mathbf x)^T \mathbf R^{-1} (\mathbf z-\mathbf H\mathbf x). \]

Therefore,

\[ \boxed{ \text{Weighted Least Squares} \Longleftrightarrow \text{Maximum Likelihood} } \]

for linear measurements with Gaussian errors.

This important relationship will be developed more generally in 4.3.6 Maximum Likelihood and MAP Estimation.

4.3.2.24 Correlated Measurement Errors

The measurement covariance need not be diagonal.

For example,

\[ \mathbf R= \begin{bmatrix} \sigma_1^2 & R_{12}\\ R_{21} & \sigma_2^2 \end{bmatrix}. \]

If

\[ R_{12}\neq0, \]

the measurement errors are correlated.

The WLS objective remains

\[ \boxed{ J = \mathbf r^T\mathbf R^{-1}\mathbf r. } \]

The complete covariance must be used.

Ignoring the off-diagonal terms effectively treats correlated observations as if they contained more independent information than they actually do.

Engineering interpretation: collecting more measurements does not automatically provide proportionally more information if their errors are strongly correlated.

Engineering interpretation: Correlation can arise when measurements share a clock, calibration process, atmospheric correction, image-processing chain, or other common error source. Off-diagonal covariance terms therefore represent real shared uncertainty and should not be discarded merely to simplify the algebra.

4.3.2.25 Whitening the Measurements

Weighted least squares can also be converted into an ordinary least-squares problem.

Suppose a matrix \mathbf L satisfies

\[ \mathbf R=\mathbf L\mathbf L^T. \]

For example, \mathbf L may be obtained from a Cholesky factorization.

Define whitened quantities:

\[ \boxed{ \tilde{\mathbf z} = \mathbf L^{-1}\mathbf z } \]

and

\[ \boxed{ \tilde{\mathbf H} = \mathbf L^{-1}\mathbf H. } \]

Then

\[ \tilde{\mathbf z} = \tilde{\mathbf H}\mathbf x + \tilde{\mathbf v}, \]

where the transformed measurement noise has identity covariance:

\[ E[ \tilde{\mathbf v}\tilde{\mathbf v}^T ] = \mathbf I. \]

Therefore,

\[ \boxed{ \text{Weighted Least Squares} \rightarrow \text{Ordinary Least Squares in Whitened Coordinates}. } \]

This is both mathematically useful and numerically important.

4.3.2.26 Measurements with Different Units

Suppose the measurement vector contains

\[ \mathbf z= \begin{bmatrix} \rho\\ \dot\rho\\ \alpha \end{bmatrix}, \]

where

• \rho is range in metres,
• \dot\rho is range rate in metres per second,
• \alpha is an angle in radians.

Directly minimizing

\[ r_\rho^2+r_{\dot\rho}^2+r_\alpha^2 \]

can be physically misleading because the residuals have different units and scales.

Weighting by covariance gives

\[ J = \frac{r_\rho^2}{\sigma_\rho^2} + \frac{r_{\dot\rho}^2}{\sigma_{\dot\rho}^2} + \frac{r_\alpha^2}{\sigma_\alpha^2} \]

when the measurements are independent.

Each residual is then measured relative to its expected uncertainty.

Thus,

\[ \boxed{ \text{Covariance Weighting Makes Heterogeneous Measurements Statistically Comparable}. } \]

4.3.2.27 A Priori Information

Suppose an approximate state is already available:

\[ \bar{\mathbf x}_0 \]

with covariance

\[ \mathbf P_0. \]

The measurement cost is

\[ J_m = (\mathbf z-\mathbf H\mathbf x)^T \mathbf R^{-1} (\mathbf z-\mathbf H\mathbf x). \]

The prior cost is

\[ J_p = (\mathbf x-\bar{\mathbf x}_0)^T \mathbf P_0^{-1} (\mathbf x-\bar{\mathbf x}_0). \]

The combined objective becomes

\[ \boxed{ \begin{aligned} J(\mathbf x) ={}& (\mathbf z-\mathbf H\mathbf x)^T \mathbf R^{-1} (\mathbf z-\mathbf H\mathbf x) \\ &+ (\mathbf x-\bar{\mathbf x}_0)^T \mathbf P_0^{-1} (\mathbf x-\bar{\mathbf x}_0). \end{aligned} } \]

This combines information from the prior estimate and new observations.

4.3.2.28 Least Squares with A Priori Constraint

Differentiating the combined cost gives

\[ \left( \mathbf H^T\mathbf R^{-1}\mathbf H + \mathbf P_0^{-1} \right) \hat{\mathbf x} = \mathbf H^T\mathbf R^{-1}\mathbf z + \mathbf P_0^{-1}\bar{\mathbf x}_0. \]

Therefore,

\[ \boxed{ \hat{\mathbf x} = \left( \mathbf H^T\mathbf R^{-1}\mathbf H + \mathbf P_0^{-1} \right)^{-1} \left( \mathbf H^T\mathbf R^{-1}\mathbf z + \mathbf P_0^{-1}\bar{\mathbf x}_0 \right). } \]

The posterior covariance is

\[ \boxed{ \mathbf P = \left( \mathbf P_0^{-1} + \mathbf H^T\mathbf R^{-1}\mathbf H \right)^{-1}. } \]

This equation gives a particularly important information interpretation:

\[ \boxed{ \text{Posterior Information} = \text{Prior Information} + \text{Measurement Information}. } \]

That is,

\[ \boxed{ \mathbf P^{-1} = \mathbf P_0^{-1} + \mathbf H^T\mathbf R^{-1}\mathbf H. } \]

This relationship becomes central when connecting batch estimation to Bayesian and recursive estimation.

4.3.2.29 Strength of the Prior

If

\[ \mathbf P_0 \]

is small, the prior is precise.

Then

\[ \mathbf P_0^{-1} \]

is large, giving the prior strong influence.

Conversely,

\[ \mathbf P_0 \]

large implies weak prior information.

Therefore,

\[ \boxed{ \mathbf P_0\downarrow \Rightarrow \text{Prior Weight}\uparrow. } \]

But a prior should not be made artificially precise simply to stabilize an otherwise weak estimation problem.

An incorrect overconfident prior can dominate valid measurement information and bias the solution.

4.3.2.30 A Priori Information as Additional Observations

Prior information can also be interpreted as another observation equation:

\[ \bar{\mathbf x}_0 = \mathbf I\mathbf x+\mathbf e_0, \]

where

\[ E[\mathbf e_0\mathbf e_0^T] = \mathbf P_0. \]

The complete estimation problem can therefore be written as

\[ \begin{bmatrix} \mathbf z\\ \bar{\mathbf x}_0 \end{bmatrix} = \begin{bmatrix} \mathbf H\\ \mathbf I \end{bmatrix} \mathbf x + \begin{bmatrix} \mathbf v\\ \mathbf e_0 \end{bmatrix}. \]

This demonstrates that prior information and measurements can be treated within the same mathematical framework.

4.3.2.31 Rank Requirement

For a unique unconstrained least-squares solution,

\[ \mathbf H \]

must have sufficient independent columns.

For n unknowns,

\[ \boxed{ \operatorname{rank}(\mathbf H)=n. } \]

If

\[ \operatorname{rank}(\mathbf H)

then some combination of unknown quantities cannot be uniquely determined from the observations.

Consequently,

\[ \mathbf H^T\mathbf H \]

is singular in ordinary least squares.

Similarly,

\[ \mathbf H^T\mathbf R^{-1}\mathbf H \]

will be singular when the weighted problem lacks sufficient independent information.

Thus,

\[ \boxed{ \text{Full Measurement Count} \neq \text{Full Information Rank}. } \]

4.3.2.32 Rank Deficiency Example

Consider

\[ \mathbf H= \begin{bmatrix} 1&1\\ 2&2\\ 3&3 \end{bmatrix}. \]

The second column is exactly equal to the first.

Therefore,

\[ \operatorname{rank}(\mathbf H)=1 \]

even though there are two unknowns.

If

\[ \mathbf x= \begin{bmatrix} x_1\\ x_2 \end{bmatrix}, \]

the measurements depend only on

\[ x_1+x_2. \]

The estimator may determine that combination but cannot distinguish x_1 and x_2 individually.

This is the linear algebra version of the identifiability problem discussed in Section 4.3.1.

4.3.2.33 More Measurements Do Not Guarantee Better Geometry

Suppose one observation row is

\[ \mathbf h^T. \]

Repeating exactly the same observation many times gives

\[ \mathbf H= \begin{bmatrix} \mathbf h^T\\ \mathbf h^T\\ \mathbf h^T\\ \vdots \end{bmatrix}. \]

The number of measurements increases, but the measurement geometry does not fundamentally change.

If some state direction was unobservable from \mathbf h, repeating the same geometry does not create sensitivity to that direction.

Therefore,

\[ \boxed{ \text{More Measurements} \neq \text{More Independent Measurement Directions}. } \]

Diverse measurement geometry is often more valuable than simply increasing measurement count.

4.3.2.34 Conditioning

A matrix may have full rank and still be difficult to solve numerically.

Suppose two columns of \mathbf H are nearly linearly dependent.

Then

\[ \mathbf H^T\mathbf H \]

may be nearly singular.

The problem is then ill-conditioned.

A common measure is the condition number:

\[ \boxed{ \kappa(\mathbf H) = \frac{\sigma_{\max}} {\sigma_{\min}}, } \]

where \sigma_{\max} and \sigma_{\min} are the largest and smallest singular values.

If

\[ \kappa(\mathbf H)\gg1, \]

small changes in the observations can produce large changes in the estimated state.

Thus,

\[ \boxed{ \text{Full Rank} \not\Rightarrow \text{Well-Conditioned Estimation}. } \]

Engineering interpretation: Rank asks whether a unique solution exists in principle; conditioning asks how sensitive that solution is to small perturbations. A full-rank matrix can therefore still produce a fragile estimate when its smallest singular value is very small.

4.3.2.35 Why Forming the Normal Equations Can Be Numerically Risky

The analytical expression

\[ \hat{\mathbf x} = (\mathbf H^T\mathbf H)^{-1} \mathbf H^T\mathbf z \]

is useful for understanding least squares.

However, explicitly forming and inverting

\[ \mathbf H^T\mathbf H \]

is not always the best numerical implementation.

Approximately,

\[ \boxed{ \kappa(\mathbf H^T\mathbf H) \approx \kappa(\mathbf H)^2. } \]

Therefore, forming the normal equations can worsen conditioning.

A problem that is moderately ill-conditioned in \mathbf H can become substantially more difficult when expressed through \mathbf H^T\mathbf H.

Engineering interpretation: the closed-form inverse is an important theoretical equation, but production estimation software should not automatically implement the equation literally using an explicit matrix inverse.

4.3.2.36 Never Invert a Matrix Just Because the Equation Contains an Inverse

Mathematical derivations often write

\[ \mathbf A^{-1}\mathbf b. \]

Numerically, it is generally preferable to solve

\[ \boxed{ \mathbf A\mathbf x=\mathbf b } \]

directly rather than explicitly computing

\[ \mathbf A^{-1}. \]

Thus the normal equations should generally be treated as

\[ \boxed{ (\mathbf H^T\mathbf W\mathbf H)\hat{\mathbf x} = \mathbf H^T\mathbf W\mathbf z } \]

and solved using an appropriate linear-system method.

This distinction between mathematical notation and numerical implementation is important throughout estimation theory.

4.3.2.37 QR Factorization

A more numerically robust solution can be obtained using QR factorization.

Suppose

\[ \boxed{ \mathbf H=\mathbf Q\mathbf R_H, } \]

where

\[ \mathbf Q^T\mathbf Q=\mathbf I \]

and \mathbf R_H is upper triangular.

Then

\[ \|\mathbf z-\mathbf H\mathbf x\|^2 = \|\mathbf Q^T\mathbf z-\mathbf R_H\mathbf x\|^2. \]

For a full-column-rank problem, the least-squares solution can be obtained from the triangular system

\[ \boxed{ \mathbf R_H\hat{\mathbf x} = \mathbf Q^T\mathbf z. } \]

QR avoids explicitly forming \mathbf H^T\mathbf H and is generally more numerically stable than solving the normal equations directly.

4.3.2.38 Cholesky Factorization

If the weighted information matrix

\[ \mathbf A = \mathbf H^T\mathbf R^{-1}\mathbf H \]

is symmetric positive definite, it can be factored as

\[ \boxed{ \mathbf A=\mathbf L\mathbf L^T. } \]

The normal equations

\[ \mathbf A\hat{\mathbf x}=\mathbf b \]

can then be solved using two triangular systems:

\[ \mathbf L\mathbf y=\mathbf b \]

followed by

\[ \mathbf L^T\hat{\mathbf x}=\mathbf y. \]

Cholesky factorization is computationally efficient when its assumptions are satisfied.

However, if the problem is nearly rank deficient, QR or SVD may provide better numerical diagnostics.

4.3.2.39 Singular Value Decomposition

The Singular Value Decomposition writes

\[ \boxed{ \mathbf H = \mathbf U\boldsymbol{\Sigma}\mathbf V^T. } \]

The diagonal matrix \boldsymbol{\Sigma} contains singular values

\[ \sigma_1\ge\sigma_2\ge\cdots\ge\sigma_n\ge0. \]

Very small singular values indicate weakly determined directions.

If a singular value is zero, the corresponding direction is unobservable in the linear estimation problem.

Using the pseudoinverse,

\[ \boxed{ \hat{\mathbf x} = \mathbf H^{+}\mathbf z. } \]

SVD is particularly valuable for

• detecting rank deficiency,
• examining weak estimation directions,
• computing minimum-norm solutions,
• diagnosing parameter correlation,
• understanding numerical conditioning.

Thus SVD is not only a solver; it is also a powerful estimation-analysis tool.

4.3.2.40 Pseudoinverse and Minimum-Norm Solution

When \mathbf H is not square or lacks full rank, the Moore-Penrose pseudoinverse may be used:

\[ \boxed{ \mathbf H^+ = \mathbf V\boldsymbol{\Sigma}^+\mathbf U^T. } \]

Then

\[ \boxed{ \hat{\mathbf x} = \mathbf H^+\mathbf z. } \]

For an underdetermined system, multiple states may produce the same minimum residual.

The pseudoinverse selects the minimum-Euclidean-norm solution under the standard formulation.

However,

\[ \boxed{ \text{A Pseudoinverse Producing a Number} \neq \text{Physical Identifiability}. } \]

If the data do not independently determine a parameter, numerical machinery cannot create that missing physical information.

4.3.2.41 Scaling of the Estimation Problem

Suppose the state contains quantities with very different magnitudes:

\[ \mathbf x= \begin{bmatrix} r\\ v\\ b \end{bmatrix} \]

where position may be of order 10^6, velocity 10^3, and a sensor bias 10^{-5}.

Poor scaling can make numerical conditioning worse.

Introduce scaled variables:

\[ \boxed{ \mathbf x=\mathbf S_x\tilde{\mathbf x}. } \]

Then

\[ \mathbf z = \mathbf H\mathbf S_x\tilde{\mathbf x} + \mathbf v. \]

Define

\[ \tilde{\mathbf H} = \mathbf H\mathbf S_x. \]

The scaled estimation problem may have much better numerical behaviour.

Scaling does not create new information, but it can prevent numerical representation from unnecessarily degrading an otherwise valid problem.

Engineering interpretation: Scaling is especially important when one solve contains heterogeneous quantities such as position, velocity, angles, sensor biases, and dimensionless coefficients. Representative scaling can make design-matrix columns comparable in magnitude and improve numerical behaviour without changing the underlying physics.

4.3.2.42 Regularization

If an estimation problem is poorly conditioned, additional information or constraints may sometimes be introduced through regularization.

A common form is

\[ \boxed{ J(\mathbf x) = \|\mathbf z-\mathbf H\mathbf x\|^2 + \lambda\|\mathbf L\mathbf x\|^2. } \]

Here,

• \lambda controls regularization strength,
• \mathbf L determines which combinations of states are penalized.

The solution satisfies

\[ \boxed{ ( \mathbf H^T\mathbf H + \lambda\mathbf L^T\mathbf L ) \hat{\mathbf x} = \mathbf H^T\mathbf z. } \]

For

\[ \mathbf L=\mathbf I, \]

this becomes a common Tikhonov or ridge-type regularization.

Regularization can stabilize weak problems, but it introduces additional assumptions.

Therefore,

\[ \boxed{ \text{Regularization Stabilizes a Problem by Adding Information or Preference}; \quad \text{It Does Not Create Measurement Information}. } \]

4.3.2.43 Hard Constraints versus Soft Constraints

Suppose a quantity must satisfy

\[ \mathbf C\mathbf x=\mathbf d. \]

This is a hard equality constraint.

It can be enforced directly using constrained optimization methods.

Alternatively, suppose the constraint is approximate:

\[ \mathbf C\mathbf x \approx \mathbf d. \]

Then it may be included as a weighted residual:

\[ \boxed{ J_c = (\mathbf C\mathbf x-\mathbf d)^T \mathbf R_c^{-1} (\mathbf C\mathbf x-\mathbf d). } \]

This behaves like a soft constraint.

Thus,

\[ \boxed{ \text{Hard Constraint} \rightarrow \text{Must Be Satisfied} } \]

whereas

\[ \boxed{ \text{Soft Constraint} \rightarrow \text{Penalizes Departure According to Uncertainty}. } \]

The distinction is important because uncertain engineering knowledge should not automatically be treated as exact truth.

4.3.2.44 Combining Multiple Sensor Types

Suppose one sensor produces

\[ \mathbf z_A = \mathbf H_A\mathbf x+\mathbf v_A \]

and another produces

\[ \mathbf z_B = \mathbf H_B\mathbf x+\mathbf v_B. \]

The observations can be stacked:

\[ \boxed{ \mathbf z= \begin{bmatrix} \mathbf z_A\\ \mathbf z_B \end{bmatrix}, \qquad \mathbf H= \begin{bmatrix} \mathbf H_A\\ \mathbf H_B \end{bmatrix}. } \]

If the sensor errors are independent,

\[ \mathbf R= \begin{bmatrix} \mathbf R_A&0\\ 0&\mathbf R_B \end{bmatrix}. \]

Then the same WLS equation estimates the common state:

\[ \boxed{ \hat{\mathbf x} = (\mathbf H^T\mathbf R^{-1}\mathbf H)^{-1} \mathbf H^T\mathbf R^{-1}\mathbf z. } \]

This shows that multisensor estimation is mathematically not simply “averaging sensors.” Each sensor contributes information according to its measurement model, sensitivity, geometry, and uncertainty.

4.3.2.45 Information Contribution from Multiple Independent Sensors

For independent sensor groups,

\[ \mathbf H= \begin{bmatrix} \mathbf H_A\\ \mathbf H_B \end{bmatrix} \]

and

\[ \mathbf R= \begin{bmatrix} \mathbf R_A&0\\ 0&\mathbf R_B \end{bmatrix}. \]

Then

\[ \mathbf H^T\mathbf R^{-1}\mathbf H = \mathbf H_A^T\mathbf R_A^{-1}\mathbf H_A + \mathbf H_B^T\mathbf R_B^{-1}\mathbf H_B. \]

Therefore,

\[ \boxed{ \boldsymbol{\Lambda} = \boldsymbol{\Lambda}_A + \boldsymbol{\Lambda}_B. } \]

This gives an elegant interpretation:

\[ \boxed{ \text{Independent Sensor Information Adds}. } \]

A sensor can be valuable even when it is individually less precise if it provides information in a state direction that other sensors constrain poorly.

4.3.2.46 General Relative-Navigation Example

Consider a general relative-navigation problem with translational state

\[ \mathbf x_{\mathrm{rel}} = \begin{bmatrix} x\\y\\z \end{bmatrix}. \]

Suppose several linearized measurements around a nominal relative state satisfy

\[ \delta\mathbf z = \mathbf H \delta\mathbf x_{\mathrm{rel}} + \mathbf v. \]

Different observations may provide different directional sensitivities.

For example, one measurement may constrain primarily the line-of-sight direction while another contributes information about range or another independent geometric direction.

The weighted correction is

\[ \boxed{ \delta\hat{\mathbf x}_{\mathrm{rel}} = ( \mathbf H^T\mathbf R^{-1}\mathbf H )^{-1} \mathbf H^T\mathbf R^{-1} \delta\mathbf z. } \]

The updated relative-state estimate is

\[ \boxed{ \hat{\mathbf x}_{\mathrm{rel}}^{+} = \hat{\mathbf x}_{\mathrm{rel}}^{-} + \delta\hat{\mathbf x}_{\mathrm{rel}}. } \]

The important principle is general:

\[ \boxed{ \text{Complementary Measurement Geometry} \rightarrow \text{Improved Relative-State Estimation}. } \]

This is not tied to any particular rendezvous or proximity-operations implementation.

4.3.2.47 Parameter-Estimation Example

Least squares can estimate parameters as well as states.

Suppose a measured quantity satisfies

\[ z_i = a_i\theta+v_i, \]

where \theta is an unknown model parameter.

Stacking the observations gives

\[ \mathbf z = \mathbf H\theta+\mathbf v, \]

with

\[ \mathbf H= \begin{bmatrix} a_1\\ a_2\\ \vdots\\ a_m \end{bmatrix}. \]

The WLS estimate is

\[ \boxed{ \hat\theta = \frac{ \mathbf H^T\mathbf R^{-1}\mathbf z }{ \mathbf H^T\mathbf R^{-1}\mathbf H }. } \]

The parameter is estimable only if the observations are sufficiently sensitive to it.

If all

\[ a_i\approx0, \]

then the parameter has little influence on the measurements and its uncertainty will remain large.

4.3.2.48 Estimating Multiple Parameters

Suppose

\[ \boldsymbol{\theta} = \begin{bmatrix} \theta_1\\ \theta_2 \end{bmatrix} \]

and

\[ \mathbf z = \mathbf H_\theta\boldsymbol{\theta} + \mathbf v. \]

If the two columns of \mathbf H_\theta are almost parallel,

\[ \mathbf h_{\theta_1} \approx c\mathbf h_{\theta_2}, \]

the two parameters produce nearly indistinguishable effects on the observations.

The information matrix

\[ \mathbf H_\theta^T \mathbf R^{-1} \mathbf H_\theta \]

then becomes poorly conditioned.

The resulting parameter covariance may contain large variances and strong correlations.

Thus,

\[ \boxed{ \text{Good Residual Fit} \neq \text{Well-Identified Individual Parameters}. } \]

4.3.2.49 Residual Analysis After Estimation

After solving for

\[ \hat{\mathbf x}, \]

the post-fit residual is

\[ \boxed{ \hat{\mathbf r} = \mathbf z-\mathbf H\hat{\mathbf x}. } \]

Residual analysis can reveal useful information.

Questions include:

• Are residuals centred approximately around zero?
• Are there persistent biases?
• Are some sensors producing systematically larger residuals?
• Are residual magnitudes consistent with the assumed measurement uncertainty?
• Do residuals show time dependence?
• Are there outliers?
• Do residuals correlate with geometry, temperature, operating mode, or another variable?
• Does a particular measurement type dominate the cost?

A small total least-squares cost does not by itself prove that the physical model is correct.

Residual structure should therefore be examined, not only the final scalar cost.

Engineering interpretation: Residual analysis should be performed by measurement type as well as globally. A single total cost can hide systematic behaviour in one sensor class, whereas grouped residual statistics can reveal a biased sensor, an incorrect covariance model, geometry-dependent mismatch, or an unmodelled effect.

4.3.2.50 Pre-Fit and Post-Fit Residuals

If an initial estimate

\[ \mathbf x_0 \]

is available, the pre-fit residual is

\[ \boxed{ \mathbf r_{\text{pre}} = \mathbf z-\mathbf H\mathbf x_0. } \]

After estimating a correction

\[ \Delta\hat{\mathbf x}, \]

the updated state is

\[ \hat{\mathbf x} = \mathbf x_0+\Delta\hat{\mathbf x}. \]

The post-fit residual is

\[ \boxed{ \mathbf r_{\text{post}} = \mathbf z-\mathbf H\hat{\mathbf x}. } \]

Comparing pre-fit and post-fit residuals shows how much of the observation mismatch the estimated state correction explains.

However, forcing residuals to become extremely small can also indicate overfitting if too many parameters are estimated.

4.3.2.51 Outliers and the Limitation of Squared Residuals

Ordinary and weighted least squares square the residuals.

Therefore, a very large residual receives a disproportionately large penalty.

For example,

\[ r=10 \]

contributes

\[ r^2=100. \]

This means least squares can be sensitive to outliers.

Outliers may arise from

• corrupted measurements,
• tracking loss,
• false feature association,
• multipath,
• sensor saturation,
• timing errors,
• data-processing faults.

Classical least squares assumes the adopted residual model is appropriate.

When outliers are significant, robust estimation methods may be required.

Examples include robust loss functions, data screening, gating, or iterative reweighting.

These methods should not be confused with ordinary statistical weighting by measurement covariance.

4.3.2.52 Weighted Least Squares Is Not Arbitrary Weight Tuning

A common misunderstanding is to treat WLS weights as tuning gains.

In statistically formulated WLS,

\[ \boxed{ \mathbf W=\mathbf R^{-1}. } \]

Therefore, the weights should represent the measurement uncertainty model.

Artificially increasing a sensor's weight means mathematically claiming that its measurements are more precise.

Artificially reducing its weight means claiming greater uncertainty.

Thus,

\[ \boxed{ \text{Weight Selection} \leftrightarrow \text{Measurement-Uncertainty Assumption}. } \]

If weights are adjusted empirically, that choice should be documented and validated rather than presented as a physical covariance without justification.

4.3.2.53 Measurement Quantity versus Measurement Information

Two measurements with similar numerical accuracy may provide very different state information.

Suppose measurement A has very small noise but weak sensitivity to x_2.

Measurement B has somewhat larger noise but strong sensitivity to x_2.

Measurement B may contribute more information about x_2.

The information matrix captures both effects:

\[ \boxed{ \boldsymbol{\Lambda} = \mathbf H^T\mathbf R^{-1}\mathbf H. } \]

Therefore,

\[ \boxed{ \text{Useful Estimation Information} = \text{Measurement Precision} + \text{State Sensitivity} + \text{Independent Geometry}. } \]

Sensor accuracy alone does not determine estimator performance.

4.3.2.54 Adding a Measurement

Suppose the current information matrix is

\[ \boldsymbol{\Lambda}_{old}. \]

A new independent scalar measurement has sensitivity row

\[ \mathbf h^T \]

and variance

\[ \sigma^2. \]

Its information contribution is

\[ \boxed{ \Delta\boldsymbol{\Lambda} = \frac{1}{\sigma^2} \mathbf h\mathbf h^T. } \]

Therefore,

\[ \boxed{ \boldsymbol{\Lambda}_{new} = \boldsymbol{\Lambda}_{old} + \frac{1}{\sigma^2} \mathbf h\mathbf h^T. } \]

This equation clearly shows why a new measurement is useful only in the state directions to which it is sensitive.

It also begins to reveal the connection between batch least squares and recursive estimation.

4.3.2.55 Connection to Recursive Estimation

Suppose measurements arrive one at a time.

Instead of rebuilding the entire batch problem from the beginning, the estimator can update the accumulated information.

Conceptually,

\[ \boxed{ \text{Old Information} + \text{New Measurement Information} \rightarrow \text{New Estimate}. } \]

This is one of the mathematical foundations of recursive least squares and Gaussian recursive estimation.

The detailed connection will be developed in 4.3.4 Recursive Estimation from Least-Squares and Bayesian Viewpoints.

4.3.2.56 Static Least Squares versus Dynamic Estimation

The least-squares formulation developed so far assumes a common unknown vector \mathbf x appearing directly in the observation equations:

\[ \mathbf z=\mathbf H\mathbf x+\mathbf v. \]

But spacecraft observations are frequently collected at different times:

\[ t_1,t_2,\ldots,t_N. \]

The state at each measurement epoch is related to an earlier state through the dynamics.

For example,

\[ \delta\mathbf x(t_k) = \boldsymbol{\Phi}(t_k,t_0) \delta\mathbf x(t_0). \]

Then the measurement sensitivity to the initial state becomes

\[ \boxed{ \frac{\partial h_k} {\partial\mathbf x_0} = \frac{\partial h_k} {\partial\mathbf x_k} \boldsymbol{\Phi}(t_k,t_0). } \]

This converts measurements distributed along a trajectory into information about a common initial state.

That is the subject of 4.3.3 Batch Dynamic State Estimation.

4.3.2.57 Static versus Dynamic Design Matrix

For static estimation,

\[ \mathbf H = \frac{\partial h} {\partial\mathbf x}. \]

For a dynamic batch problem, the effective design matrix contains sensitivities propagated through time:

\[ \boxed{ \mathbf H_k^{(0)} = \mathbf H_k \boldsymbol{\Phi}(t_k,t_0). } \]

Measurements can then be stacked:

\[ \begin{bmatrix} \delta\mathbf z_1\\ \delta\mathbf z_2\\ \vdots\\ \delta\mathbf z_N \end{bmatrix} = \begin{bmatrix} \mathbf H_1\boldsymbol{\Phi}_{1,0}\\ \mathbf H_2\boldsymbol{\Phi}_{2,0}\\ \vdots\\ \mathbf H_N\boldsymbol{\Phi}_{N,0} \end{bmatrix} \delta\mathbf x_0 + \mathbf v. \]

This has exactly the same algebraic structure as linear least squares.

Therefore,

\[ \boxed{ \text{Dynamic Batch Estimation} = \text{Least Squares} + \text{Dynamic Sensitivity Propagation}. } \]

4.3.2.58 Linear Least Squares versus Nonlinear Least Squares

Linear least squares assumes

\[ \mathbf z=\mathbf H\mathbf x+\mathbf v \]

with a fixed \mathbf H.

For a nonlinear model,

\[ \mathbf z=h(\mathbf x)+\mathbf v, \]

the measurement equation can be linearized around the current estimate:

\[ h(\mathbf x+\delta\mathbf x) \approx h(\mathbf x) + \mathbf H\delta\mathbf x. \]

Then

\[ \boxed{ \delta\mathbf z \approx \mathbf H\delta\mathbf x+\mathbf v. } \]

A linear least-squares problem is solved for the correction \delta\mathbf x, the nominal state is updated, and the process is repeated.

Conceptually,

\[ \boxed{ \text{Nonlinear Estimation} \rightarrow \text{Linearize} \rightarrow \text{Solve LS Correction} \rightarrow \text{Update} \rightarrow \text{Relinearize}. } \]

This becomes the foundation of 4.3.5 Nonlinear Least Squares and Differential Corrections.

4.3.2.59 When Linear Least Squares Is Appropriate

Linear least squares is appropriate when

• the observation equations are genuinely linear in the unknowns,
• or a local linear model is sufficiently accurate,
• the measurement uncertainty model is reasonably characterized,
• the unknowns are sufficiently constrained by the observations,
• numerical conditioning is acceptable.

Typical applications include

• static calibration,
• linear parameter estimation,
• bias estimation,
• fitting linearized measurement corrections,
• local state corrections,
• components of orbit-determination algorithms,
• components of navigation and sensor-calibration algorithms.

The mathematics is simple, but the assumptions must still be checked.

4.3.2.60 When Least Squares Can Give a Misleading Solution

A numerical least-squares solution does not automatically mean the estimation problem is physically valid.

Potential problems include

• incorrect measurement models,
• incorrect covariance assumptions,
• measurement bias,
• outliers,
• poor geometry,
• rank deficiency,
• nearly dependent parameters,
• incorrect reference frames,
• inconsistent units,
• weak sensitivities,
• model mismatch,
• nonlinear effects outside the linearization region,
• overconfident priors,
• numerical ill-conditioning.

Therefore,

\[ \boxed{ \text{Successful Matrix Solution} \neq \text{Successful Estimation Problem}. } \]

The physical meaning of the solution must always be assessed.

4.3.2.61 Practical Linear/WLS Estimation Workflow

A disciplined workflow is:

1. Define the unknown state or parameter vector \mathbf x.
2. Assemble the measurement vector \mathbf z.
3. Define the observation equations.
4. Construct the design matrix \mathbf H.
5. Verify units, reference frames, epochs, and sign conventions.
6. Define measurement covariance \mathbf R.
7. Include measurement correlations where required.
8. Determine whether prior information (\bar{\mathbf x}_0,\mathbf P_0) exists.
9. Examine dimensions and rank.
10. Examine column scaling.
11. Examine singular values and condition number when necessary.
12. Form residuals.
13. Whiten the problem if useful.
14. Select an appropriate numerical solver.
15. Solve for \hat{\mathbf x} without unnecessary explicit matrix inversion.
16. Compute or assess solution covariance.
17. Calculate post-fit residuals.
18. Examine parameter correlations.
19. Check physical plausibility.
20. Investigate outliers and residual structure.
21. Verify that the solution is supported by actual measurement information.

Conceptually,

\[ \boxed{ \begin{array}{c} \text{Define Unknowns}\\ \downarrow\\ \text{Assemble Observations}\\ \downarrow\\ \text{Construct }\mathbf H\\ \downarrow\\ \text{Define }\mathbf R\text{ and Prior Information}\\ \downarrow\\ \text{Check Rank / Geometry / Conditioning}\\ \downarrow\\ \text{Form Weighted LS Problem}\\ \downarrow\\ \text{Solve Numerically}\\ \downarrow\\ \text{Compute Covariance}\\ \downarrow\\ \text{Analyse Residuals and Correlations}\\ \downarrow\\ \text{Validate Physical Meaning} \end{array} } \]

Engineering interpretation: The workflow is iterative in practice. If residuals, covariance, singular values, or physical checks reveal a weakness, return to the measurement model, state definition, covariance model, or data selection rather than only changing the numerical solver.

4.3.2.62 Complete Linear Least-Squares Architecture

The complete structure can be summarized as

\[ \boxed{ \begin{array}{c} \text{Unknown State / Parameters }\mathbf x\\ \downarrow\\ \text{Linear Observation Model}\\ \mathbf z=\mathbf H\mathbf x+\mathbf v\\ \downarrow\\ \text{Predicted Observation}\\ \hat{\mathbf z}=\mathbf H\hat{\mathbf x}\\ \downarrow\\ \text{Residual}\\ \mathbf r=\mathbf z-\hat{\mathbf z}\\ \downarrow\\ \text{Measurement Covariance }\mathbf R\\ \downarrow\\ \text{Weighted Cost}\\ J=\mathbf r^T\mathbf R^{-1}\mathbf r\\ \downarrow\\ \text{Information Matrix}\\ \boldsymbol{\Lambda} = \mathbf H^T\mathbf R^{-1}\mathbf H\\ \downarrow\\ \text{Solve}\\ \hat{\mathbf x} = \boldsymbol{\Lambda}^{-1} \mathbf H^T\mathbf R^{-1}\mathbf z\\ \downarrow\\ \text{Estimate Covariance}\\ \mathbf P_{\hat x} = \boldsymbol{\Lambda}^{-1}\\ \downarrow\\ \text{Residual / Correlation / Conditioning Assessment} \end{array} } \]

The inverse notation above expresses the analytical result. In numerical software, direct factorization and linear solves are generally preferable to explicitly computing matrix inverses.

4.3.2.63 Ordinary LS versus Weighted LS versus Prior-Constrained LS

The three major forms can be summarized as follows.

Their conceptual progression is

\[ \boxed{ \text{Residual Fit} \rightarrow \text{Uncertainty-Weighted Residual Fit} \rightarrow \text{Measurement + Prior Information}. } \]

4.3.2.64 Key Equations

Linear measurement model

\[ \boxed{ \mathbf z=\mathbf H\mathbf x+\mathbf v } \]

Predicted observation

\[ \boxed{ \hat{\mathbf z}=\mathbf H\hat{\mathbf x} } \]

Residual

\[ \boxed{ \mathbf r=\mathbf z-\mathbf H\mathbf x } \]

Ordinary least-squares cost

\[ \boxed{ J=\mathbf r^T\mathbf r } \]

Ordinary normal equations

\[ \boxed{ \mathbf H^T\mathbf H\hat{\mathbf x} = \mathbf H^T\mathbf z } \]

Ordinary least-squares solution

\[ \boxed{ \hat{\mathbf x}_{LS} = (\mathbf H^T\mathbf H)^{-1} \mathbf H^T\mathbf z } \]

Measurement covariance

\[ \boxed{ \mathbf R = E[\mathbf v\mathbf v^T] } \]

Weight matrix

\[ \boxed{ \mathbf W=\mathbf R^{-1} } \]

Weighted least-squares cost

\[ \boxed{ J = \mathbf r^T\mathbf R^{-1}\mathbf r } \]

Weighted normal equations

\[ \boxed{ \mathbf H^T\mathbf R^{-1}\mathbf H \hat{\mathbf x} = \mathbf H^T\mathbf R^{-1}\mathbf z } \]

Weighted least-squares solution

\[ \boxed{ \hat{\mathbf x}_{WLS} = (\mathbf H^T\mathbf R^{-1}\mathbf H)^{-1} \mathbf H^T\mathbf R^{-1}\mathbf z } \]

Information matrix

\[ \boxed{ \boldsymbol{\Lambda} = \mathbf H^T\mathbf R^{-1}\mathbf H } \]

Estimate covariance

\[ \boxed{ \mathbf P_{\hat x} = (\mathbf H^T\mathbf R^{-1}\mathbf H)^{-1} } \]

Prior-constrained information matrix

\[ \boxed{ \mathbf P^{-1} = \mathbf P_0^{-1} + \mathbf H^T\mathbf R^{-1}\mathbf H } \]

Prior-constrained estimate

\[ \boxed{ \hat{\mathbf x} = \mathbf P \left( \mathbf P_0^{-1}\bar{\mathbf x}_0 + \mathbf H^T\mathbf R^{-1}\mathbf z \right) } \]

Rank requirement

\[ \boxed{ \operatorname{rank}(\mathbf H)=n } \]

for a unique unconstrained estimate of n unknown quantities.

SVD

\[ \boxed{ \mathbf H = \mathbf U\boldsymbol{\Sigma}\mathbf V^T } \]

Pseudoinverse solution

\[ \boxed{ \hat{\mathbf x} = \mathbf H^+\mathbf z } \]

4.3.2.65 Key Takeaways

Linear least squares begins with

\[ \boxed{ \mathbf z=\mathbf H\mathbf x+\mathbf v. } \]

When noisy measurements make the system inconsistent, the goal is not to satisfy every observation exactly.

Instead,

\[ \boxed{ \hat{\mathbf x} = \arg\min_{\mathbf x} \|\mathbf z-\mathbf H\mathbf x\|^2. } \]

Ordinary least squares treats all residuals equally.

Weighted least squares incorporates measurement uncertainty:

\[ \boxed{ J = \mathbf r^T\mathbf R^{-1}\mathbf r. } \]

Its solution is

\[ \boxed{ \hat{\mathbf x} = ( \mathbf H^T\mathbf R^{-1}\mathbf H )^{-1} \mathbf H^T\mathbf R^{-1}\mathbf z. } \]

The matrix

\[ \boxed{ \mathbf H^T\mathbf R^{-1}\mathbf H } \]

represents measurement information about the unknown state, while

\[ \boxed{ \mathbf P_{\hat x} = ( \mathbf H^T\mathbf R^{-1}\mathbf H )^{-1} } \]

represents the corresponding estimate uncertainty under the assumed linear statistical model.

Prior information adds naturally:

\[ \boxed{ \text{Posterior Information} = \text{Prior Information} + \text{Measurement Information}. } \]

A large number of observations does not guarantee a good solution.

\[ \boxed{ \text{Measurement Count} \neq \text{Independent Information}. } \]

The design matrix must contain sufficient independent sensitivity:

\[ \boxed{ \operatorname{rank}(\mathbf H)=n. } \]

Even full-rank problems may be poorly conditioned.

Therefore rank, singular values, geometry, scaling, parameter correlation, and numerical conditioning must also be examined.

The analytical inverse formulas are useful for understanding the mathematics, but numerical implementations should generally use stable matrix factorizations or linear solves such as QR, Cholesky, or SVD, rather than blindly forming explicit inverses.

The central engineering principle is

\[ \boxed{ \text{Least Squares Finds the Best Fit Allowed by the Model and Data;} \quad \text{It Does Not Guarantee That the Model or Data Are Correct}. } \]

4.3.2.66 What Comes Next

The formulation developed in this section assumed that all measurements could be related directly to a common unknown vector through

\[ \mathbf z=\mathbf H\mathbf x+\mathbf v. \]

Spacecraft estimation problems are often more complicated because observations occur at different epochs while the spacecraft state evolves continuously.

A measurement at time t_k depends on

\[ \mathbf x(t_k), \]

but the quantity being estimated may be the state at an earlier reference epoch:

\[ \mathbf x(t_0). \]

The dynamics provide the connection:

\[ \boxed{ \delta\mathbf x(t_k) = \boldsymbol{\Phi}(t_k,t_0) \delta\mathbf x(t_0). } \]

Consequently,

\[ \boxed{ \frac{\partial h_k}{\partial\mathbf x_0} = \frac{\partial h_k}{\partial\mathbf x_k} \boldsymbol{\Phi}(t_k,t_0). } \]

Measurements distributed over an observation arc can therefore all contribute information about the same initial state.

This leads directly to:

\[ \boxed{ \textbf{4.3.3 Batch Dynamic State Estimation} } \]

where least squares is extended from a static linear problem to a dynamic trajectory-estimation problem involving multiple measurement epochs, trajectory propagation, initial-state estimation, state-transition matrices, dynamic sensitivities, and batch correction.

Next: 4.3.3 Batch Dynamic State Estimation

Linear and weighted least squares provide the static mathematical foundation for estimating unknown quantities from noisy observations. Spacecraft estimation problems, however, commonly involve measurements acquired at different times while the state evolves according to the equations of motion.

The next section extends least squares to a dynamic observation arc. It develops trajectory propagation, reference-epoch state estimation, measurement stacking, State Transition Matrix propagation, dynamic design matrices, sensitivity propagation, batch normal equations, iterative correction, convergence, and post-fit assessment.