Numerical Conditioning and Robust Solution Methods

Normal equations, QR, Cholesky, SVD, square-root least squares, condition number, scaling, regularization, rank deficiency, covariance stability, robust factorization choices, numerical diagnostics, and practical solver selection.

4.3.17 Numerical Conditioning and Robust Solution Methods

Numerical conditioning and robust solution methods deal with an important practical question in estimation:

Even if the estimation problem is theoretically solvable, can it be solved reliably and accurately on a computer?

A mathematically valid estimator can still produce poor results if the underlying matrices are badly scaled, nearly singular, rank deficient, or solved using numerically fragile methods.

The central idea is

\[ \boxed{ \text{Good Estimation} = \text{Good Mathematical Model} + \text{Good Numerical Conditioning} + \text{Robust Linear Algebra} } \]

This section focuses on:

These ideas apply to:

4.3.17.1 Why Numerical Conditioning Matters

Consider the linearized estimation problem

\[ \boxed{ \mathbf r = \mathbf G \delta\mathbf x + \mathbf v } \]

where

Even if \mathbf G has full rank, the solution can be very sensitive if its columns are nearly linearly dependent.

Then small errors in

can produce large changes in the estimated correction.

Thus,

\[ \boxed{ \text{Full Rank} \neq \text{Numerically Reliable} } \]

4.3.17.2 Well-Conditioned versus Ill-Conditioned Problems

A problem is well conditioned if small perturbations in the input produce small perturbations in the solution.

It is ill conditioned if small input perturbations can produce large changes in the solution.

Conceptually,

\[ \boxed{ \text{Small Measurement Error} \rightarrow \text{Small State Error} } \]

for a well-conditioned system,

whereas

\[ \boxed{ \text{Small Measurement Error} \rightarrow \text{Potentially Large State Error} } \]

for an ill-conditioned one.

4.3.17.3 Numerical Conditioning versus Physical Observability

These are related but not identical.

Poor physical geometry often produces poor numerical conditioning.

However, numerical problems can also arise from poor variable scaling even when the physical information content is adequate.

Therefore,

\[ \boxed{ \text{Poor Conditioning} = \text{Possible Geometry Problem} + \text{Possible Scaling Problem} + \text{Possible Solver Problem} } \]

The cause should be diagnosed rather than assumed.

4.3.17.4 Linear Least-Squares Problem

The basic weighted least-squares problem is

\[ \boxed{ \min_{\delta\mathbf x} \; J = ( \mathbf r-\mathbf G\delta\mathbf x )^T \mathbf W ( \mathbf r-\mathbf G\delta\mathbf x ) } \]

where

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

If measurements are whitened, this becomes

\[ \boxed{ \min_{\delta\mathbf x} \| \tilde{\mathbf r} - \tilde{\mathbf G} \delta\mathbf x \|^2. } \]

The choice of how this system is solved matters greatly for numerical robustness.

4.3.17.5 Normal Equations

The classical normal equations are

\[ \boxed{ \mathbf G^T \mathbf W \mathbf G \, \delta\mathbf x = \mathbf G^T \mathbf W \mathbf r. } \]

Define

\[ \boxed{ \mathbf N = \mathbf G^T \mathbf W \mathbf G } \]

and

\[ \boxed{ \mathbf b = \mathbf G^T \mathbf W \mathbf r. } \]

Then

\[ \boxed{ \mathbf N \delta\mathbf x = \mathbf b. } \]

4.3.17.6 Why Normal Equations Are Attractive

Normal equations are popular because they

For many moderate-size estimation problems, they can work well.

4.3.17.7 Main Numerical Weakness of Normal Equations

If

\[ \kappa(\mathbf G) \]

is the condition number of the design matrix, then approximately

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

Thus the normal equations square the conditioning problem.

For example,

\[ \kappa(\mathbf G) = 10^6 \]

can produce

\[ \kappa( \mathbf G^T\mathbf G ) \sim 10^{12}. \]

This can severely amplify floating-point error.

4.3.17.8 Avoid Explicit Matrix Inversion

The formal expression

\[ \delta\mathbf x = ( \mathbf G^T \mathbf W \mathbf G )^{-1} \mathbf G^T \mathbf W \mathbf r \]

is useful mathematically.

However, in computation it is usually preferable to solve

\[ \boxed{ \mathbf N\delta\mathbf x = \mathbf b } \]

directly.

Thus,

\[ \boxed{ \text{Solve the Linear System} \quad\text{instead of explicitly computing}\quad \mathbf N^{-1}. } \]

4.3.17.9 Why Explicit Inversion Is Discouraged

Explicit matrix inversion can

Most numerical libraries provide stable linear solvers.

Therefore code such as

\[ \texttt{x = inv(A)*b} \]

is generally inferior to solving

\[ \texttt{A\textbackslash b} \]

in MATLAB.

4.3.17.10 MATLAB Backslash Operator

In MATLAB,

x = A\b;

selects an appropriate linear solver based on the matrix structure.

For a least-squares problem,

x = G\r;

can be preferable to manually forming

inv(G'*G)*G'*r

because MATLAB can use numerically stable factorization internally.

4.3.17.11 Weighted Whitening

Given

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

suppose

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

Define whitened quantities:

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

and

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

Then

\[ \boxed{ \min \| \tilde{\mathbf r} - \tilde{\mathbf G} \delta\mathbf x \|^2 } \]

is equivalent to weighted least squares.

4.3.17.12 Why Whitening Helps

Whitening

For independent measurements with standard deviations \sigma_i,

\[ \boxed{ \tilde r_i = \frac{r_i}{\sigma_i}. } \]

4.3.17.13 QR Decomposition

For the whitened design matrix,

\[ \tilde{\mathbf G}, \]

QR decomposition gives

\[ \boxed{ \tilde{\mathbf G} = \mathbf Q \mathbf R } \]

where

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

and \mathbf R is upper triangular.

4.3.17.14 QR Least-Squares Solution

The problem

\[ \min \| \tilde{\mathbf r} - \tilde{\mathbf G} \delta\mathbf x \|^2 \]

becomes

\[ \min \| \mathbf Q^T \tilde{\mathbf r} - \mathbf R \delta\mathbf x \|^2. \]

For full column rank,

\[ \boxed{ \mathbf R \delta\mathbf x = \mathbf Q^T \tilde{\mathbf r}. } \]

This triangular system can be solved efficiently by back substitution.

4.3.17.15 Why QR Is More Stable than Normal Equations

QR avoids forming

\[ \mathbf G^T\mathbf G. \]

Therefore it avoids squaring the condition number.

Conceptually,

\[ \boxed{ \text{QR} \rightarrow \text{Better Numerical Stability} } \]

than

\[ \boxed{ \text{Normal Equations} } \]

for moderately ill-conditioned problems.

4.3.17.16 Householder QR

One common implementation uses Householder transformations.

A Householder reflector has the form

\[ \boxed{ \mathbf H = \mathbf I - 2 \frac{ \mathbf v\mathbf v^T }{ \mathbf v^T\mathbf v }. } \]

It transforms a vector while preserving norm.

Householder QR is widely used because it is robust and numerically stable.

4.3.17.17 Givens Rotations

Another QR implementation uses Givens rotations.

For selected coordinates i,j,

\[ \boxed{ \mathbf G_{ij} = \begin{bmatrix} c & s\\ -s & c \end{bmatrix} } \]

with

\[ c^2+s^2=1. \]

Givens rotations are especially useful when measurements are processed incrementally or when sparse structure is important.

4.3.17.18 QR with Column Pivoting

If rank deficiency is possible, use column-pivoted QR:

\[ \boxed{ \mathbf G \mathbf P = \mathbf Q \mathbf R. } \]

The permutation matrix \mathbf P reorders columns so that the strongest independent directions are processed first.

This helps reveal numerical rank.

4.3.17.19 Rank-Revealing QR

In rank-revealing QR, the diagonal elements of \mathbf R provide information about independent versus weak columns.

If some diagonal values become extremely small, the system may be rank deficient or nearly so.

This is a useful alternative to SVD when computation is more constrained.

4.3.17.20 Cholesky Factorization

If the normal matrix

\[ \mathbf N = \mathbf G^T\mathbf W\mathbf G \]

is symmetric positive definite, then

\[ \boxed{ \mathbf N = \mathbf L \mathbf L^T } \]

or

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

This is the Cholesky factorization.

4.3.17.21 Cholesky Solution

Solve

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

first, then

\[ \mathbf L^T \delta\mathbf x = \mathbf y. \]

Cholesky is efficient because it exploits symmetry and positive definiteness.

4.3.17.22 Advantages of Cholesky

Cholesky provides

It is often useful for well-conditioned information or covariance matrices.

4.3.17.23 Limitation of Cholesky

Cholesky does not fix the conditioning problem created by forming the normal matrix.

If

\[ \mathbf N = \mathbf G^T\mathbf G, \]

the squared condition number remains.

Also, if \mathbf N is singular or not positive definite because of rank deficiency or numerical error, standard Cholesky can fail.

4.3.17.24 Cholesky Failure as a Diagnostic

If Cholesky unexpectedly fails, possible causes include

The failure should be investigated rather than merely bypassed.

4.3.17.25 Singular Value Decomposition

The singular value decomposition is

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

The singular values satisfy

\[ \boxed{ \sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_n \ge 0. } \]

SVD is one of the strongest tools for diagnosing numerical conditioning.

4.3.17.26 SVD Least-Squares Solution

For full rank,

\[ \boxed{ \delta\mathbf x = \mathbf V \boldsymbol\Sigma^{-1} \mathbf U^T \mathbf r. } \]

This makes the role of singular values explicit.

A measurement-space component associated with a small \sigma_i is amplified by

\[ \frac{1}{\sigma_i}. \]

4.3.17.27 Why Small Singular Values Are Dangerous

If

\[ \sigma_i \ll1, \]

then

\[ \frac{1}{\sigma_i} \gg1. \]

Measurement noise in that direction can therefore produce a very large parameter correction.

Thus,

\[ \boxed{ \text{Small Singular Value} \rightarrow \text{Noise Amplification}. } \]

4.3.17.28 SVD and Rank Deficiency

If

\[ \sigma_i=0, \]

the corresponding state direction is unobservable in the linearized system.

If

\[ \sigma_i \approx0, \]

it is weakly observable.

Therefore SVD simultaneously provides information about

4.3.17.29 Pseudoinverse

For a rank-deficient matrix,

\[ \mathbf G^+ \]

denotes the Moore-Penrose pseudoinverse.

Using the SVD,

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

where nonzero singular values are inverted.

Then

\[ \boxed{ \delta\mathbf x = \mathbf G^+ \mathbf r. } \]

4.3.17.30 Meaning of the Pseudoinverse Solution

For an underdetermined or rank-deficient system, the pseudoinverse gives the minimum-norm solution among possible solutions.

However,

\[ \boxed{ \text{Pseudoinverse} \neq \text{Recovered Lost Information}. } \]

It provides a mathematically defined solution but cannot create observability that the measurements do not contain.

4.3.17.31 Truncated SVD

One robust approach is to discard very small singular values.

If

\[ \sigma_i<\tau, \]

set its inverse contribution to zero.

Then

\[ \boxed{ \delta\mathbf x_{\rm TSVD} = \sum_{\sigma_i\ge\tau} \frac{ \mathbf v_i \mathbf u_i^T \mathbf r }{ \sigma_i }. } \]

This is truncated SVD.

4.3.17.32 Interpretation of Truncated SVD

TSVD says:

Do not attempt to estimate state combinations that the data constrain too weakly.

This reduces noise amplification but sacrifices resolution in weak directions.

Thus,

\[ \boxed{ \text{Robustness} \leftrightarrow \text{Information Retention}. } \]

4.3.17.33 Choosing the SVD Threshold

A numerical tolerance may depend on

A common numerical form is related to

\[ \boxed{ \tau \sim \epsilon_{\rm mach} \, \max(m,n) \, \sigma_{\max}. } \]

But an engineering threshold may need to be stricter than the purely numerical threshold.

4.3.17.34 Numerical Rank versus Physical Rank

A matrix may be mathematically full rank yet physically useless in one direction.

Therefore distinguish:

\[ \boxed{ \text{Mathematical Rank} } \]

from

\[ \boxed{ \text{Practical Numerical Rank}. } \]

Practical rank depends on how small a singular value can become before noise dominates its information.

4.3.17.35 Condition Number

For a nonsingular matrix,

\[ \boxed{ \kappa(\mathbf G) = \|\mathbf G\| \, \|\mathbf G^{-1}\|. } \]

Using the 2-norm,

\[ \boxed{ \kappa_2(\mathbf G) = \frac{ \sigma_{\max} }{ \sigma_{\min} }. } \]

4.3.17.36 Interpretation of the Condition Number

Approximately,

\[ \boxed{ \frac{ \|\delta\mathbf x_{\rm error}\| }{ \|\mathbf x\| } \lesssim \kappa(\mathbf G) \, \frac{ \|\delta\mathbf b\| }{ \|\mathbf b\| } } \]

for small perturbations in a linear solve.

Thus the condition number estimates how much input error may be amplified.

4.3.17.37 Condition Number and Significant Digits

If

\[ \kappa \sim10^p, \]

one may lose roughly p decimal digits of numerical accuracy in the worst case.

For example,

\[ \kappa\sim10^{10} \]

can be extremely problematic in double precision.

This is only a rough interpretation, but useful diagnostically.

4.3.17.38 Condition Number Is Unit Sensitive

If the state vector contains variables with very different units, the raw condition number may reflect scaling rather than physical observability.

For example:

Therefore scaling is essential before interpreting conditioning physically.

4.3.17.39 State Scaling

Define

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

where \mathbf S_x contains characteristic state scales.

Then

\[ \mathbf G \delta\mathbf x = \mathbf G \mathbf S_x \delta\tilde{\mathbf x}. \]

The scaled design matrix is

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

4.3.17.40 Example of State Scaling

Suppose

\[ r \sim10^7\ {\rm m}, \]

\[ v \sim10^3\ {\rm m/s}, \]

and a calibration parameter is

\[ \theta \sim10^{-4}. \]

Solving directly can produce very differently scaled columns.

Instead choose scales such as

\[ \boxed{ \mathbf S_x = \operatorname{diag} ( 10^7, 10^7, 10^7, 10^3, 10^3, 10^3, 10^{-4} ). } \]

4.3.17.41 Measurement Scaling

Measurements can also have very different units.

Examples include

Whitening by measurement covariance naturally provides statistically meaningful scaling.

Thus measurement whitening and state scaling solve different problems:

\[ \boxed{ \text{Whitening} \rightarrow \text{Measurement Scaling} } \]

\[ \boxed{ \text{State Scaling} \rightarrow \text{Unknown-Variable Scaling}. } \]

4.3.17.42 Column Normalization

A simple diagnostic is to normalize design-matrix columns:

\[ \boxed{ \tilde{\mathbf g}_j = \frac{ \mathbf g_j }{ \|\mathbf g_j\| }. } \]

This helps assess whether poor conditioning comes from different column magnitudes or from near-collinearity.

However, physical state scaling is generally more meaningful than arbitrary normalization for the final estimator.

4.3.17.43 Scaling Does Not Create Information

Scaling may greatly improve numerical conditioning.

But it does not change the actual measurement information.

If two columns are physically parallel,

\[ \mathbf g_1 = c\mathbf g_2, \]

scaling cannot make them independent.

Thus,

\[ \boxed{ \text{Scaling Fixes Numerical Problems} \quad\text{not}\quad \text{Observability Problems}. } \]

4.3.17.44 Square-Root Least Squares

Instead of working with the information matrix

\[ \mathbf G^T \mathbf W \mathbf G, \]

square-root least-squares methods work with a factor such as

\[ \boxed{ \mathbf R_s } \]

satisfying

\[ \boxed{ \mathbf R_s^T \mathbf R_s = \mathbf G^T \mathbf W \mathbf G. } \]

These methods avoid explicitly squaring conditioning where possible.

4.3.17.45 Square-Root Information Form

The information matrix is

\[ \boxed{ \boldsymbol\Lambda = \mathbf P^{-1}. } \]

Instead of storing \boldsymbol\Lambda, store a triangular factor

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

This can improve numerical stability and preserve positive definiteness.

4.3.17.46 Square-Root Covariance Form

Alternatively,

\[ \boxed{ \mathbf P = \mathbf S \mathbf S^T } \]

is maintained directly through its factor \mathbf S.

This is common in square-root Kalman filters.

4.3.17.47 Why Square-Root Methods Are Robust

Square-root methods generally provide better control over

They are especially useful for high-precision estimation and long recursive filtering.

4.3.17.48 Square-Root Filtering

In a conventional Kalman filter,

\[ \mathbf P \]

is propagated and updated directly.

In a square-root filter, a matrix factor is propagated instead.

Conceptually,

\[ \boxed{ \mathbf P = \mathbf S\mathbf S^T } \]

is maintained without repeatedly reconstructing and subtracting large nearly equal covariance terms.

4.3.17.49 Covariance Symmetry

A theoretical covariance satisfies

\[ \boxed{ \mathbf P = \mathbf P^T. } \]

Roundoff may cause

\[ \mathbf P \neq \mathbf P^T \]

numerically.

A common cleanup is

\[ \boxed{ \mathbf P \leftarrow \frac12 ( \mathbf P+\mathbf P^T ). } \]

However, repeated large asymmetry indicates a deeper numerical problem.

4.3.17.50 Positive Semidefinite Covariance

A covariance matrix should satisfy

\[ \boxed{ \mathbf x^T \mathbf P \mathbf x \ge0 } \]

for every \mathbf x.

If negative eigenvalues appear because of numerical error, the covariance is no longer physically meaningful.

Square-root and Joseph-form implementations help reduce this risk.

4.3.17.51 Joseph Covariance Update

The Joseph form is

\[ \boxed{ \mathbf P^+ = ( \mathbf I-\mathbf K\mathbf H ) \mathbf P^- ( \mathbf I-\mathbf K\mathbf H )^T + \mathbf K \mathbf R \mathbf K^T. } \]

This form is more numerically robust than the simplified update

\[ \mathbf P^+ = ( \mathbf I-\mathbf K\mathbf H ) \mathbf P^-. \]

4.3.17.52 Why the Joseph Form Helps

The Joseph form explicitly contains two positive-semidefinite contributions.

Therefore it better preserves covariance

under finite-precision arithmetic.

4.3.17.53 Rank Deficiency

A matrix is rank deficient when

\[ \boxed{ \operatorname{rank}(\mathbf G) < n. } \]

Then the least-squares solution is not unique.

Physically, some state or parameter combination is not determined by the data.

4.3.17.54 Exact Rank Deficiency

If a nonzero vector \mathbf q satisfies

\[ \boxed{ \mathbf G\mathbf q = \mathbf0, } \]

then moving the estimate along \mathbf q does not change the predicted measurements to first order.

No numerical method can uniquely recover that component from the measurements alone.

4.3.17.55 Near Rank Deficiency

More commonly,

\[ \boxed{ \mathbf G\mathbf q \approx \mathbf0. } \]

Then the state direction is only weakly constrained.

The solution may exist uniquely but be highly sensitive to noise.

This is often more dangerous than exact rank deficiency because a numerical solver may return a plausible-looking but unreliable number.

4.3.17.56 Detecting Rank Deficiency

Useful diagnostics include

Do not rely only on a solver returning successfully.

4.3.17.57 MATLAB Rank Diagnostics

Useful MATLAB commands include:

rank(G)
cond(G)
svd(G)
rcond(N)

and for column-pivoted QR:

[Q,R,E] = qr(G,0);

These should be interpreted after sensible scaling.

4.3.17.58 Reciprocal Condition Number

MATLAB's

rcond(A)

estimates approximately

\[ \boxed{ \frac{1}{\kappa(A)}. } \]

A very small rcond indicates a nearly singular matrix.

For example,

\[ \texttt{rcond} \rightarrow0 \]

means the solve may be unreliable.

4.3.17.59 Regularization

When the problem is weakly conditioned, one may introduce additional information or penalties.

A common regularized least-squares problem is

\[ \boxed{ \min_{\delta\mathbf x} \left[ \| \mathbf G\delta\mathbf x-\mathbf r \|^2 + \lambda \| \mathbf L\delta\mathbf x \|^2 \right]. } \]

This is Tikhonov regularization.

4.3.17.60 Tikhonov Normal Equations

The regularized system becomes

\[ \boxed{ ( \mathbf G^T\mathbf G + \lambda \mathbf L^T\mathbf L ) \delta\mathbf x = \mathbf G^T\mathbf r. } \]

For

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

this becomes ridge-type regularization.

4.3.17.61 Interpretation of Regularization

Regularization penalizes large or undesirable corrections.

It reduces noise amplification in weak directions.

But the solution becomes biased toward the imposed regularization model.

Thus,

\[ \boxed{ \text{Lower Variance} \leftrightarrow \text{Higher Bias}. } \]

4.3.17.62 Regularization Is Additional Information

Regularization should not be viewed as a purely numerical trick.

It imposes assumptions about the solution.

For example,

\[ \lambda \|\delta\mathbf x\|^2 \]

expresses a preference for small corrections.

Therefore its physical meaning should be understood and reported.

4.3.17.63 Prior Information as Regularization

Suppose

\[ \delta\mathbf x \sim \mathcal N( \delta\mathbf x_0, \mathbf P_0 ). \]

Then the prior contribution is

\[ \boxed{ ( \delta\mathbf x-\delta\mathbf x_0 )^T \mathbf P_0^{-1} ( \delta\mathbf x-\delta\mathbf x_0 ). } \]

This is statistically meaningful regularization.

4.3.17.64 Prior-Augmented Normal Matrix

With a prior,

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

The prior can stabilize weak measurement directions.

However, those directions remain weakly observed by the measurement data themselves.

4.3.17.65 Levenberg-Marquardt Regularization

For nonlinear least squares,

\[ \boxed{ ( \mathbf G^T \mathbf W \mathbf G + \lambda\mathbf I ) \delta\mathbf x = \mathbf G^T \mathbf W \mathbf r. } \]

This is closely related to Levenberg-Marquardt damping.

It improves robustness when the current linearization is poor.

4.3.17.66 Gauss-Newton versus Levenberg-Marquardt

Gauss-Newton uses

\[ \boxed{ \mathbf G^T \mathbf W \mathbf G. } \]

Levenberg-Marquardt adds

\[ \boxed{ \lambda\mathbf I. } \]

For small \lambda, behavior approaches Gauss-Newton.

For large \lambda, the step becomes more conservative.

4.3.17.67 Damping

Another simple strategy is

\[ \boxed{ \mathbf x_{k+1} = \mathbf x_k + \alpha \delta\mathbf x } \]

with

\[ 0<\alpha\le1. \]

Damping controls nonlinear iteration size but does not by itself fix fundamental rank deficiency.

4.3.17.68 Line Search

A line search chooses \alpha to reduce the objective function:

\[ \boxed{ J( \mathbf x+\alpha\delta\mathbf x ) < J(\mathbf x). } \]

This can improve convergence when a full Newton or Gauss-Newton step is too large.

4.3.17.69 Trust-Region Methods

Instead of accepting an arbitrarily large correction, trust-region methods constrain

\[ \boxed{ \| \delta\mathbf x \| \le \Delta. } \]

The local linear or quadratic model is trusted only within a neighborhood.

This can provide strong nonlinear robustness.

4.3.17.70 Scaling in Nonlinear Estimation

If parameter scales are very different, nonlinear solvers may take inappropriate steps.

Scaled variables

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

allow the optimization algorithm to treat comparable dimensionless corrections.

Scaling is often as important for nonlinear convergence as for matrix conditioning.

4.3.17.71 Residual Scaling

Likewise, residuals should be statistically scaled:

\[ \boxed{ \tilde{\mathbf r} = \mathbf R^{-1/2} \mathbf r. } \]

Then one unit of residual represents approximately one standard deviation regardless of measurement type.

4.3.17.72 Robust Least Squares

Classical least squares uses

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

This gives large residuals very strong influence.

Robust least squares replaces the quadratic penalty by

\[ \boxed{ J = \sum_i \rho(r_i). } \]

4.3.17.73 Huber Loss

A common robust function is the Huber loss:

\[ \boxed{ \rho(r) = \begin{cases} \frac12r^2, & |r|\le\delta\\ \delta( |r|-\frac12\delta ), & |r|>\delta. \end{cases} } \]

Small residuals remain quadratic.

Large residuals grow approximately linearly.

4.3.17.74 Why Robust Loss Helps

Robust losses reduce the influence of

They can improve solution stability when a small number of residuals are abnormal.

4.3.17.75 Robust Methods Do Not Fix Model Error

A systematic residual caused by

should not simply be downweighted.

Otherwise a robust solver may hide an important modeling problem.

Thus,

\[ \boxed{ \text{Robust Loss} \neq \text{Substitute for Correct Modeling}. } \]

4.3.17.76 Iteratively Reweighted Least Squares

Many robust estimators can be implemented through IRLS.

At iteration j,

\[ \boxed{ \mathbf W^{(j)} = \operatorname{diag} ( w_1^{(j)}, \ldots, w_m^{(j)} ) } \]

is constructed from current residuals.

Then solve a weighted least-squares problem.

The process repeats until weights and state stabilize.

4.3.17.77 Outlier Rejection versus Robust Weighting

Outlier rejection removes a measurement completely.

Robust weighting retains it but reduces its influence.

Both methods require care.

Aggressive rejection can bias the solution.

Weak weighting may allow bad data to contaminate it.

4.3.17.78 Sparse Systems

Large estimation problems often produce sparse matrices.

Examples include

Exploiting sparsity can dramatically reduce memory and computation.

4.3.17.79 Sparse QR

Sparse QR preserves sparsity better than dense factorization and avoids forming dense normal matrices.

It is useful for large least-squares systems with sparse measurement structure.

4.3.17.80 Sparse Cholesky

For sparse symmetric positive-definite information matrices,

\[ \mathbf N = \mathbf L\mathbf L^T \]

can be solved efficiently using sparse Cholesky.

Ordering strategies reduce fill-in.

4.3.17.81 Fill-In

During sparse factorization, originally zero matrix elements may become nonzero.

This is called fill-in.

Good variable ordering can substantially reduce fill-in and therefore computation.

4.3.17.82 Factor Graph Interpretation

Large estimation problems may be represented as a graph containing

The resulting least-squares system is often sparse.

Robust sparse linear algebra is therefore central to modern smoothing and factor-graph estimation.

4.3.17.83 Solver Choice by Problem Type

A practical guideline is:

No one solver is universally best.

4.3.17.84 Computational Cost

For an m\times n matrix with m\gg n,

normal-equation formation costs approximately

\[ O(mn^2), \]

while QR has similar asymptotic order but larger constants.

SVD is generally more expensive.

Therefore solver choice involves a tradeoff:

\[ \boxed{ \text{Robustness} \leftrightarrow \text{Computational Cost}. } \]

4.3.17.85 Do Not Optimize Too Early

A common mistake is choosing the fastest solver before understanding the conditioning.

During development, QR or SVD may provide valuable diagnostics.

After conditioning and rank are understood, a more efficient solver can be adopted if justified.

4.3.17.86 Floating-Point Precision

Double-precision arithmetic typically provides machine precision around

\[ \boxed{ \epsilon_{\rm mach} \approx 2.2\times10^{-16}. } \]

This does not mean solutions are accurate to 16 digits.

Conditioning determines how much precision is retained.

4.3.17.87 Single versus Double Precision

Single precision has much less numerical range and precision than double precision.

For poorly conditioned estimation problems, switching to single precision can cause significant degradation.

Onboard systems using reduced precision should therefore assess conditioning carefully.

4.3.17.88 Numerical Range

Large differences in variable magnitude can cause loss of numerical resolution even before formal singularity occurs.

Scaling helps keep variables within comparable numerical ranges.

4.3.17.89 Residual Cancellation

When predicted and measured quantities are both very large but their difference is small, subtractive cancellation can reduce numerical precision.

This can occur in some high-precision orbit or timing calculations.

Careful formulation and state centering can help.

4.3.17.90 Centered Variables

Instead of estimating absolute quantities of large magnitude, estimate perturbations around a reference:

\[ \boxed{ \mathbf x = \mathbf x_{\rm ref} + \delta\mathbf x. } \]

This often improves numerical precision and is standard in differential correction and error-state filtering.

4.3.17.91 Error-State Formulation

An error-state estimator propagates a nominal state separately from the small correction:

\[ \boxed{ \mathbf x_{\rm true} = \mathbf x_{\rm nominal} + \delta\mathbf x } \]

for additive states.

This keeps estimated corrections small and can improve conditioning.

For attitude, multiplicative error representations are typically used instead.

4.3.17.92 Regularization versus Parameter Removal

If a parameter is fundamentally unobservable, regularizing it may produce a finite number but not genuine information.

Sometimes the better approach is to remove or combine parameters.

Thus,

\[ \boxed{ \text{Unobservable Parameter} \rightarrow \text{Reparameterize or Remove} } \]

rather than automatically regularize.

4.3.17.93 Composite Parameters

If measurements depend only on

\[ \gamma = \theta_1\theta_2, \]

estimating \theta_1 and \theta_2 separately may be ill conditioned.

Estimating

\[ \boxed{ \gamma } \]

directly can greatly improve numerical and physical identifiability.

4.3.17.94 Physical Priors versus Artificial Regularization

A physical prior is supported by independent knowledge.

Artificial regularization is introduced mainly for numerical stability.

They should not be confused.

A well-documented estimator should state whether the result is

4.3.17.95 Covariance under Regularization

Once regularization or priors are added, covariance interpretation changes.

For a Bayesian prior,

\[ \boxed{ \mathbf P_{\rm post} = ( \mathbf G^T \mathbf R^{-1} \mathbf G + \mathbf P_0^{-1} )^{-1}. } \]

For arbitrary regularization, interpreting the inverse Hessian as a physical covariance may not always be valid.

4.3.17.96 Numerical Diagnostics to Report

A robust estimation analysis should consider reporting:

These diagnostics make numerical quality transparent.

4.3.17.97 Common Mistakes

Common numerical mistakes include:

  1. Explicitly computing matrix inverses.
  2. Using normal equations without checking conditioning.
  3. Ignoring that \mathbf G^T\mathbf G squares the condition number.
  4. Using Cholesky on a rank-deficient system.
  5. Ignoring scaling.
  6. Computing condition number before scaling physically different parameters.
  7. Treating solver success as proof of a trustworthy solution.
  8. Ignoring very small singular values.
  9. Regularizing an unobservable parameter without acknowledging the prior information introduced.
  10. Confusing numerical rank with physical estimability.
  11. Using single precision without evaluating conditioning.
  12. Allowing covariance to lose symmetry.
  13. Ignoring negative covariance eigenvalues.
  14. Using the simplified covariance update when numerical robustness is important.
  15. Downweighting systematic model error as if it were an outlier.
  16. Using pseudoinverse without interpreting the null space.
  17. Using SVD threshold blindly.
  18. Failing to distinguish scaling problems from geometry problems.
  19. Using excessive regularization to force apparent stability.
  20. Optimizing computational speed before establishing numerical correctness.

4.3.17.98 Practical Robust-Solution Workflow

A disciplined numerical workflow is

\[ \boxed{ \begin{array}{c} \text{Form Residuals and Jacobians}\\ \downarrow\\ \text{Apply Measurement Whitening}\\ \downarrow\\ \text{Scale States / Parameters}\\ \downarrow\\ \text{Inspect Rank and Singular Values}\\ \downarrow\\ \text{Compute Condition Number}\\ \downarrow\\ \text{Identify Weak Directions}\\ \downarrow\\ \text{Choose Solver}\\ \downarrow\\ \text{QR / Cholesky / SVD / Square-Root}\\ \downarrow\\ \text{Apply Regularization Only if Justified}\\ \downarrow\\ \text{Solve Without Explicit Inversion}\\ \downarrow\\ \text{Check Residual Reduction}\\ \downarrow\\ \text{Check Covariance Symmetry / PSD}\\ \downarrow\\ \text{Check Sensitivity to Scaling / Solver Choice}\\ \downarrow\\ \text{Validate Numerical and Physical Consistency} \end{array} } \]

4.3.17.99 Core Equations

Weighted least squares

\[ \boxed{ \min_{\delta\mathbf x} ( \mathbf r-\mathbf G\delta\mathbf x )^T \mathbf R^{-1} ( \mathbf r-\mathbf G\delta\mathbf x ) } \]

Normal equations

\[ \boxed{ \mathbf G^T \mathbf R^{-1} \mathbf G \, \delta\mathbf x = \mathbf G^T \mathbf R^{-1} \mathbf r } \]

QR factorization

\[ \boxed{ \tilde{\mathbf G} = \mathbf Q \mathbf R } \]

Cholesky factorization

\[ \boxed{ \mathbf N = \mathbf L \mathbf L^T } \]

Singular value decomposition

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

Condition number

\[ \boxed{ \kappa_2(\mathbf G) = \frac{ \sigma_{\max} }{ \sigma_{\min} } } \]

Tikhonov regularization

\[ \boxed{ ( \mathbf G^T \mathbf G + \lambda \mathbf L^T \mathbf L ) \delta\mathbf x = \mathbf G^T \mathbf r } \]

State scaling

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

Covariance factorization

\[ \boxed{ \mathbf P = \mathbf S \mathbf S^T } \]

Joseph covariance update

\[ \boxed{ \mathbf P^+ = ( \mathbf I-\mathbf K\mathbf H ) \mathbf P^- ( \mathbf I-\mathbf K\mathbf H )^T + \mathbf K \mathbf R \mathbf K^T } \]

4.3.17.100 Key Takeaways

Numerical conditioning determines whether an estimation problem can be solved reliably in finite-precision arithmetic.

A theoretically observable system may still be numerically fragile.

The classical normal equations

\[ \boxed{ \mathbf G^T \mathbf W \mathbf G \, \delta\mathbf x = \mathbf G^T \mathbf W \mathbf r } \]

are convenient but can worsen conditioning because

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

QR decomposition avoids explicitly forming the normal matrix and is usually more robust:

\[ \boxed{ \mathbf G = \mathbf Q\mathbf R. } \]

Cholesky is efficient when the system is symmetric positive definite and reasonably conditioned:

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

SVD provides the strongest diagnostic view:

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

Small singular values reveal weak or nearly unobservable estimation directions.

The condition number

\[ \boxed{ \kappa = \frac{ \sigma_{\max} }{ \sigma_{\min} } } \]

quantifies sensitivity to perturbations.

Scaling is essential because poorly chosen units can create apparent numerical problems:

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

But

\[ \boxed{ \text{Scaling Improves Numerical Representation} \neq \text{Creating New Information}. } \]

Regularization can stabilize weak problems:

\[ \boxed{ \mathbf G^T\mathbf G + \lambda\mathbf L^T\mathbf L } \]

but it introduces additional assumptions or prior information.

Square-root methods are especially valuable when numerical stability, covariance positivity, and long recursive operation matter.

The most important practical principle is:

\[ \boxed{ \text{Do Not Judge an Estimator Only by Whether the Solver Returns a Number.} } \]

A trustworthy solution should also be checked for:

The overall engineering strategy is therefore

\[ \boxed{ \text{Scale Properly} \rightarrow \text{Diagnose Rank and Conditioning} \rightarrow \text{Choose a Robust Solver} \rightarrow \text{Regularize Only When Justified} \rightarrow \text{Validate Numerically and Physically}. } \]

That is the foundation of robust numerical estimation.

Next: 4.3.18 Estimator Design Workflow and Algorithm Selection

Numerical conditioning determines whether a mathematically valid estimator can be solved reliably in finite precision. The next section brings the complete estimation framework together and explains how to choose among static, dynamic, batch, sequential, linear, nonlinear, probabilistic, smoothing, and parameter-estimation approaches.

4.3.18 Estimator Design Workflow and Algorithm Selection

Choosing among static, dynamic, batch, sequential, linear, nonlinear, probabilistic, smoothing, and parameter-estimation approaches.

Continue to Estimator Design Workflow and Algorithm Selection →