Dealing With Larger Networks

The Design Matrix
Dealing With Sparse Design Matrices
Solving The Normal Equations
Don't Invert N
Re-Order The Parameters In X
Phased And Sequential Solutions

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 =

  a1     b1 a2     b2 a3     b3

  , and Cm =

  s12     0     0 0     s12     0 0     0     s12

  (3 obs. solving 2 unk.)

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

  a12	+	a22	+	a32		a1b1	+	a2b2	+	a2b2
  s12		s22		s32		s12		s22		s32

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

  Bt C-1m m =	a12	+	a1b1	+	a22		a2b2	+	a32	+	a2b2
  	s12		s12		s22		s22		s32		s32
  	b1a1	+	b12	+	b2a2		b22	+	b2a2	+	b32
  	s12		s12		s22		s22		s32		s32

  =	1	a1 b1	[a1 b1] +	1	a2 b2	[a2 b2] +	1	a3 b3	[a3 b3]
  	s12			s22			s32

• 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:

  lm - lc =	fl	Dx2 +	fl	Dy2 +	fl	Dx4 +	fl	Dy4 +
  	fx2		fy2		fx4		fy4

  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	= -	x4 - x2	b =	fl	= -	y4 - y2
  	fx2		lc		fy2		lc
  c =	fl	= -	x4 - x2	d =	fl	= -	y4 - y2
  	fx4		lc		fy4		lc

  For this single observation:

  [ 0    0    0    a    b    0    0    0    0    c    d    0 ]	Dx1	= [ lm - lc ] + [v]
  	Dy1
  	Dz1
  	Dx2
  	Dy2
  	Dz2
  	Dx3
  	Dy3
  	Dz3
  	Dx4
  	Dy4
  	Dz4

• 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 s2	0	0	0	0	0	0	0	0	0	0	0	0	Dx1	=	1 s2	0
  		0	0	0	0	0	0	0	0	0	0	0	Dy1			0
  			0	0	0	0	0	0	0	0	0	0	Dz1			0
  				a2	ab	0	0	0	0	ac	ad	0	Dx2			a(lm - lc)
  					b2	0	0	0	0	bc	bd	0	Dy2			b(lm - lc)
  						0	0	0	0	0	0	0	Dz2			0
  							0	0	0	0	0	0	Dx3			0
  								0	0	0	0	0	Dy3			0
  			Sym						0	0	0	0	Dz3			0
  										c2	cd	0	Dx4			c(lm - lc)
  											d2	0	Dy4			d(lm - lc)
  												0	Dz4			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]

  Dx2 Dy2 Dx4 Dy4

  = [ lm - lc] + [v] (obs. Equation)

  1 s2	a2	ab	ac	ad	Dx2	=	1 s2	a(lm - lc)
  		b2	bc	bd	Dy2			b(lm - lc)
  			c2	cd	Dx4			c(lm - lc)
  				d2	Dy4			d(lm - lc)

• 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.