Dealing With Larger Networks

The Design Matrix

Consider the observation equations model and solution:

Bx = m + v with C_m

B^t C^-1_mBx = B^t C^-1_m m or Nx = l (normal equations)

= N^-1 l
For survey networks, the design matrix B is:
- Sparse - the horizontal angle and slope distance observations have 6 differential coefficients, most other observations only have 4 (horizontal distance, bearing), 2 (level difference), or 1 (direct observation of a coordinate).
- Can be very large in dimension. If there are 100 unknowns and 500 observations B will have dimension 500 * 100, where the normal matrix will only be 100 *100 ! This may become an issue with a computer implementation of network adjustment.

Fortunately the full design matrix does not need to be formed. Observations can be added directly to the normal equations one at a time if they are not correlated with other observations. Correlated observations must be added to the normal equations as a group.

Let B =
a₁     b₁

a₂     b₂

a₃     b₃
, and C_{m =}
s₁²     0     0

0     s₁²     0

0     0     s₁²
(3 obs. solving 2 unk.)

Forming the normal matrix B^t C^-1_m m gives:

a₁² + a₂² + a₃² a₁b₁ + a₂b₂ + a₂b₂

s₁² s₂² s₃² s₁² s₂² s₃²

Notice how the normal matrix is a sum of components unique to each observation, ie:

B^t C^-1_m m = a₁² + a₁b₁ + a₂² a₂b₂ + a₃² + a₂b₂

s₁² s₁² s₂² s₂² s₃² s₃²

b₁a₁ + b₁² + b₂a₂ b₂² + b₂a₂ + b₃²

s₁² s₁² s₂² s₂² s₃² s₃²

= 1
a₁

b₁
[a₁ b₁] + 1
a₂

b₂
[a₂ b₂] + 1
a₃

b₃
[a₃ b₃]

s₁² s₂² s₃²

Thus the normal matrix can be formed one observation at a time, where the observations are not correlated with each other (that is C_m is diagonal). If observations are correlated they can be added to the normal equations as a group, for example the 3 components of a GPS baseline would be added to the normal equations together.

The right hand side of the normal equation (B^t C^-1_m m ) can be formed in exactly the same manner, observation by observation. (Again with the correlation free stipulation).

Dealing With Sparse Design Matrices

Example: Consider a network with 4 unknown stations in 3D, and one observation of slope distance l with precision s between stations 2 and 4.

The observation equation is: l = (x₄ - x₂)² + (y₄ - y₂)²

Linearised:

l_m - l_c = fl Dx₂ + fl Dy₂ + fl Dx₄ + fl Dy₄ +

fx₂ fy₂ fx₄ fy₄

where:

l_m is the measured value

l_c is the length computed with assumed coordinates

The derivatives (with reference to previous notes on planar observation equations) are:

a = fl = - x₄ - x₂ b = fl = - y₄ - y₂

fx₂ l_c fy₂ l_c

c = fl = - x₄ - x₂ d = fl = - y₄ - y₂

fx₄ l_c fy₄ l_c

For this single observation:

[ 0 0 0 a b 0 0 0 0 c d 0 ] Dx₁ = [ l_m - l_c ] + [v]

Dy₁

Dz₁

Dx₂

Dy₂

Dz₂

Dx₃

Dy₃

Dz₃

Dx₄

Dy₄

Dz₄

Note the number of zeroes in the design matrix, and imagine the case if there were 100 unknowns rather than 12.

When the normal equations are formed there will be no entries for coefficient not involved in the observation - stations 1 and 3, and the heights of 2 and 4:

1

s²
0 0 0 0 0 0 0 0 0 0 0 0 Dx₁ =
1

s²
0

0 0 0 0 0 0 0 0 0 0 0 Dy₁ 0

0 0 0 0 0 0 0 0 0 0 Dz₁ 0

a² ab 0 0 0 0 ac ad 0 Dx₂ a(l_m - l_c)

b² 0 0 0 0 bc bd 0 Dy₂ b(l_m - l_c)

0 0 0 0 0 0 0 Dz₂ 0

0 0 0 0 0 0 Dx₃ 0

0 0 0 0 0 Dy₃ 0

Sym 0 0 0 0 Dz₃ 0

c² cd 0 Dx₄ c(l_m - l_c)

d² 0 Dy₄ d(l_m - l_c)

0 Dz₄ 0

The pattern of where entries in the normal equations will appear, established in the above matrix, is usually used by adjustment software to reduce the size of the design matrix. All that really needs to be computed is:

[a b c d]
Dx₂

Dy₂

Dx₄

Dy₄
= [ l_m - l_c] + [v] (obs. Equation)

1

s²
a² ab ac ad Dx₂ =
1

s²
a(l_m - l_c)

b² bc bd Dy₂ b(l_m - l_c)

c² cd Dx₄ c(l_m - l_c)

d² Dy₄ d(l_m - l_c)

Then the elements of this minimal representation of the normal equations are added (as preivously demonstrated) to the appropriate place in the full normal equations.
Note that the design matrix is also required to compute the covariance matrix of the adjusted measurements C and the residuals C (see previous notes). The above method of forming only the required part of the design matrix can be also be applied in this situation.

Solving The Normal Equations

Solving the normal equations is the most computationally expensive part of the least squares adjustment.
The set of linear equations Nx = l must be solved - either by inverting N, or solving a set of n linear equations where n is the size (rows and columns) of N.
The time (computer instructions required) taken to invert a matrix is a fuction of the size of the matrix cubed. For example if the size of the normal matrix is doubled, the time taken to solve the normal equations is increased by a factor of 8.
For large networks the solution time for the normal equations can become excessive. There are a few techniques that can be applied:

Don't Invert N
Re-order the parameters in x, to reduce the bandwidth of the matrix
Phased solutions

Don't Invert N

Nx = l can be solved without finding N^-1 (like solving a set of linear equations by repeatedly eliminating unknowns then backsubstituting , e.g. Gaussian elimination). Solving a set of linear equations in such a manner is approximately 3 times faster than solving the equations by computing a matrix inverse.
The problem with not inverting N is that you do not get the covariance matrix for the adjusted parameters (and residuals and adjusted measurements) - and the whole basis of our approach has been to use these covariance matrices to statistcally evaluate the network !
One approach is to solve Nx = l without inversion until the network converges, and only invert N on the final iteration of the solution

Re-Order The Parameters In X, To Reduce The Bandwidth Of The Matrix

Depending on network geometry and the measurement configuration, N can be sparse (have a lot of zero elements).
There are certain algorithms for inverting a matrix (for example the Cholesky decomposition) that allow a large number of the zero terms to be ignored in the computation of the matrix inverse.
Generally the order of the parameters in x is re-organised. This re-organisation does not increase or decrease the number of non-zero terms in N, but moves the zero terms into a favourable position for a more efficient matrix inversion. The aim is to minimise bandwidth (The maximum number of non-zero elements from the leading diagonal, until the rest of the elements are zero).
This technique is not always worth applying. The matrix may not be sparse enough to be worth the computation effort of re-organising the order of the parameters - there is some cost in optimizing the order of the parameters in x.
In some network situations - such as linear features, N is sparse and structured in a minimal bandwidth format naturally.

Table Of Contents

The Design Matrix

Dealing With Sparse Design Matrices

Solving The Normal Equations

Don't Invert N

Re-order the parameters in x, to reduce the bandwidth of the matrix

Phased solutions

Don't Invert N

Re-Order The Parameters In X, To Reduce The Bandwidth Of The Matrix