An Inquiry-Based Introduction to Linear Algebra and Applications

Subsection Application: Global Positioning System

You are probably familiar with the Global Positioning System (GPS). The system allows anyone with the appropriate software to accurately determine their location at any time. The applications are almost endless, including getting real-time driving directions while in your car, guiding missiles, and providing distances on golf courses.

The GPS is a worldwide radio-navigation system owned by the US government and operated by the US Air Force. GPS is one of four global navigation satellite systems. At least twenty four GPS satellites orbit the Earth at an altitude of approximately 11,000 nautical miles. The satellites are placed so that at any time at least four of them can be accessed by a GPS receiver. Each satellite carries an atomic clock to relay a time stamp along with its position in space. There are five ground stations to coordinate and ensure that the system is working properly.

The system works by triangulation, but there is also error involved in the measurements that go into determining position. Later in this section we will see how the method of least squares can be used to determine the receiver's position.

Subsection Introduction

A singular value decomposition has many applications, and in this section we discuss how a singular value decomposition can be used in image compression, to determine how sensitive a matrix can be to rounding errors in the process of row reduction, and to solve least squares problems.

Subsection Image Compression

The digital age has brought many new opportunities for the collection, analysis, and dissemination of information. Along with these opportunities come new difficulties as well. All of this digital information must be stored in some way and be retrievable in an efficient manner. A singular value decomposition of digitally stored information can be used to compress the information or clean up corrupted information. In this section we will see how a singular value decomposition can be used in image compression. While a singular value decomposition is normally used with very large matrices, we will restrict ourselves to small examples so that we can more clearly see how a singular value decomposition is applied.

The first step in compressing an image is to digitize the image. There are many ways to do this and we will consider one of the simplest ways and only work with gray-scale images, with the scale from 0 (black) to 255 (white). A digital image can be created by taking a small grid of squares (called pixels) and coloring each pixel with some shade of gray. The resolution of this grid is a measure of how many pixels are used per square inch. As an example, consider the 16 by 16 pixel picture of a flower shown in Figure 30.1.

To store this image pixel by pixel would require \(16 \times 16 = 256\) units of storage space (1 for each pixel). If we let \(M\) be the matrix whose \(i,j\)th entry is the scale of the \(i,j\)th pixel, then \(M\) is the matrix

Recall that if \(U \Sigma V^{\tr}\) is a singular value decomposition for \(M\text{,}\) then we can also write \(M\) in the form

given in (29.2). For this \(M\text{,}\) the singular values are approximately

Notice that some of these singular values are very small compared to others. As in Preview Activity 30.1, the terms with the largest singular values contain most of the information about the matrix. Thus, we shouldn't lose much information if we eliminate the small singular values. In this particular example, the last 4 singular values are significantly smaller than the rest. If we let

then we should expect the image determined by \(M_{12}\) to be close to the image made by \(M\text{.}\) The two images are presented side by side in Figure 30.2.

This small example illustrates the general idea. Suppose we had a satellite image that was \(1000 \times 1000\) pixels and we let \(M\) represent this image. If we have a singular value decomposition of this image \(M\text{,}\) say

if the rank of \(M\) is large, it is likely that many of the singular values will be very small. If we only keep \(s\) of the singular values, we can approximate \(M\) by

and store the image with only the vectors \(\sigma_1 \vu_1\text{,}\) \(\sigma_2 \vu_2\text{,}\) \(\ldots\text{,}\) \(\sigma_{s}\vu_s\text{,}\) \(\vv_1\text{,}\) \(\vv_1\text{,}\) \(\ldots\text{,}\) \(\vv_s\text{.}\) For example, if we only need 10 of the singular values of a satellite image (\(s = 10\)), then we can store the satellite image with only 20 vectors in \(\R^{1000}\) or with \(20 \times 1000 = 20,000\) numbers instead of \(1000 \times 1000 = 1,000,000\) numbers.

A similar process can be used to denoise data.⁵²

Subsection Calculating the Error in Approximating an Image

In the context where a matrix represents an image, the operator aspect of the matrix is irrelevant — we are only interested in the matrix as a holder of information. In this situation, we think of an \(m \times n\) matrix as simply a long vector in \(\R^{mn}\) where we have folded the data into a rectangular array. If we are interested in determining the error in approximating an image by a compressed image, it makes sense to use the standard norm in \(\R^{mn}\) to determine length and distance. This leads to what is called the Frobenius norm of a matrix. The Frobenius norm \(||M||_F\) of an \(m \times n\) matrix \(M = [m_{ij}]\) is

There is a natural corresponding inner product on the set of \(m \times n\) matrices (called the Frobenius product) defined by

where \(A = [a_{ij}]\) and \(B = [b_{ij}]\) are \(m \times n\) matrices. Note that

If an \(m \times n\) matrix \(M\) of rank \(r\) has a singular value decomposition \(M = U \Sigma V^{\tr}\text{,}\) we have seen that we can write \(M\) as an outer product

where the \(\vu_i\) are the columns of \(U\) and the \(\vv_j\) the columns of \(V\text{.}\) Each of the products \(\vu_i \vv_i^{\tr}\) is an \(m \times n\) matrix. Since the columns of \(\vu_i \vv_i^{\tr}\) are all scalar multiples of \(\vu_i\text{,}\) the matrix \(\vu_i \vv_i^{\tr}\) is a rank 1 matrix. So (30.3) expresses \(M\) as a sum of rank 1 matrices. Moreover, if we let \(\vx\) and \(\vw\) be \(m \times 1\) vectors and let \(\vy\) and \(\vz\) be \(n \times 1\) vectors with \(\vy = [y_1 \ y_2 \ \ldots \ y_n]^{\tr}\) and \(\vz = [z_1 \ z_2 \ \ldots \ z_n]^{\tr}\text{,}\) then

Using the vectors from the singular value decomposition of \(M\) as in (30.3) we see that

It follows that

When we used the singular value decomposition to approximate the image defined by \(M\text{,}\) we replaced \(M\) with a matrix of the form

We call \(M_k\) the rank \(k\) approximation to \(M\text{.}\) Notice that the outer product expansion in (30.5) is in fact a singular value decomposition for \(M_k\text{.}\) The error \(E_k\) in approximating \(M\) with \(M_k\) is

Once again, notice that (30.6) is a singular value decomposition for \(E_k\text{.}\) We define the relative error in approximating \(M\) with \(M_k\) as

Now (30.4) shows that

In applications, we often want to retain a certain degree of accuracy in our approximations and this error term can help us accomplish that.

In our flower example, the singular values of \(M\) are given in (30.2). The relative error in approximating \(M\) with \(M_{12}\) is

Errors (rounded to 4 decimal places) for approximating \(M\) with some of the \(M_k\) are shown in Table 30.3

Subsection The Condition Number of a Matrix

A singular value decomposition for a matrix \(A\) can tell us a lot about how difficult it is to accurately solve a system \(A \vx = \vb\text{.}\) Solutions to systems of linear equations can be very sensitive to rounding as the next exercise demonstrates.

Notice that a simple rounding in the \((2,2)\) entry of the coefficient matrix led to a significantly different solution. If there are rounding errors at any stage of the Gaussian elimination process, they can be compounded by further row operations. This is an important problem since computers can only approximate irrational numbers with rational numbers and so rounding can be critical. Finding ways of dealing with these kinds of errors is an area of on-going research in numerical linear algebra. This problem is given a name.

A matrix that is not ill-conditioned is said to be well-conditioned. Since small changes in entries of ill-conditioned matrices can lead to large errors in computations, it is an important problem in linear algebra to have a way to measure how ill-conditioned a matrix is. This idea will ultimately lead us to the condition number of a matrix.

Suppose we want to solve the system \(A \vx = \vb\text{,}\) where \(A\) is an invertible matrix. Activity 30.4 illustrates that if \(A\) is really close to being singular, then small changes in the entries of \(A\) can have significant effects on the solution to the system. So the system can be very hard to solve accurately if \(A\) is close to singular. It is important to have a sense of how “good” we can expect any calculated solution to be. Suppose we think we solve the system \(A \vx = \vb\) but, through rounding error in our calculation of \(A\text{,}\) get a solution \(\vx'\) so that \(A \vx' = \vb'\text{,}\) where \(\vb'\) is not exactly \(\vb\text{.}\) Let \(\Delta \vx\) be the error in our calculated solution and \(\Delta \vb\) the difference between \(\vb'\) and \(\vb\text{.}\) We would like to know how large the error \(||\Delta \vx||\) can be. But this isn't exactly the right question. We could scale everything to make \(||\Delta \vx||\) as large as we want. What we really need is a measure of the relative error \(\frac{||\Delta \vx||}{||\vx||}\text{,}\) or how big the error is compared to \(||\vx||\) itself. More specifically, we want to know how large the relative error in \(\Delta \vx\) is compared to the relative error in \(\Delta \vb\text{.}\) In other words, we want to know how good the relative error in \(\Delta \vb\) is as a predictor of the relative error in \(\Delta \vx\) (we may have some control over the relative error in \(\Delta \vb\text{,}\) perhaps by keeping more significant digits). So we want know if there is a best constant \(C\) such that

This best constant \(C\) is the condition number — a measure of how well the relative error in \(\Delta \vb\) predicts the relative error in \(\Delta \vx\text{.}\) How can we find \(C\text{?}\)

Since \(A \vx' = \vb'\) we have

Distributing on the left and using the fact that \(A\vx = \vb\) gives us

We return for a moment to the operator norm of a matrix. This is an appropriate norm to use here since we are considering \(A\) to be a transformation. Recall that if \(A\) is an \(m \times n\) matrix, we defined the operator norm of \(A\) to be

One important property that the norm has is that if the product \(AB\) is defined, then

To see why, notice that

Now \(\frac{||A(B\vx)||}{||B\vx||} \leq ||A||\) and \(\frac{||B\vx||}{||\vx||} \leq ||B||\) by the definition of the norm, so we conclude that

for every \(\vx\text{.}\) Thus,

Now we can find the condition number. From \(A \Delta \vx = \Delta \vb\) we have

so

Similarly, \(\vb = A\vx\) implies that \(||\vb|| \leq ||A|| \ ||\vx||\) or

Combining (30.8) and (30.9) gives

This constant \(||A^{-1}|| \ ||A||\) is the best bound and so is called the condition number of \(A\text{.}\)

How does a singular value decomposition tell us about the condition number of a matrix? Recall that the maximum value of \(||A\vx||\) for \(\vx\) on the unit \(n\)-sphere is \(\sigma_1\text{.}\) So \(||A|| = \sigma_1\text{.}\) If \(A\) is an invertible matrix and \(A = U \Sigma V^{\tr}\) is a singular value decomposition for \(A\text{,}\) then

where

Now \(V \Sigma^{-1} U^{\tr}\) is a singular value decomposition for \(A^{-1}\) with the diagonal entries in reverse order, so

Therefore, the condition number of \(A\) is

Subsection Pseudoinverses

Not every matrix is invertible, so we cannot always solve a matrix equation \(A \vx = \vb\text{.}\) However, every matrix has a pseudoinverse \(A^+\) that acts something like an inverse. Even when we can't solve a matrix equation \(A \vx = \vb\) because \(\vb\) isn't in \(\Col A\text{,}\) we can use the pseudoinverse of \(A\) to “solve” the equation \(A \vx = \vb\) with the “solution” \(A^+ \vb\text{.}\) While not an exact solution, \(A^+ \vb\) turns out to be the best approximation to a solution in the least squares sense. We will use the singular value decomposition to find the pseudoinverse of a matrix.

Only some square matrices have inverses. However, every matrix has a pseudoinverse. A pseudoinverse \(A^{+}\) of a matrix \(A\) provides something like an inverse when a matrix doesn't have an inverse. Pseudoinverses are useful to approximate solutions to linear systems. If \(A\) is invertible, then the equation \(A \vx = \vb\) has the solution \(\vx = A^{-1}\vb\text{,}\) but when \(A\) is not invertible and \(\vb\) is not in \(\Col A\text{,}\) then the equation \(A \vx = \vb\) has no solution. In the invertible case of an \(n \times n\) matrix \(A\text{,}\) there is a matrix \(B\) so that \(AB = BA = I_n\text{.}\) This also implies that \(BAB = B\) and \(ABA = A\text{.}\) To mimic this situation when \(A\) is not invertible, we search for a matrix \(A^+\) (a pseudoinverse of \(A\)) so that \(AA^+A = A\) and \(A^+AA^+ = A^+\text{,}\) as we saw in Preview Activity 30.7. Then it turns out that \(A^+\) acts something like an inverse for \(A\text{.}\) In this case, we approximate the solution to \(A \vx = \vb\) by \(\vx^* = A^+\vb\text{,}\) and we will see that the vector \(A\vx^* = AA^+\vb\) turns out to be the vector in \(\Col A\) that is closest to \(\vb\) in the least squares sense.

A reasonable question to ask is how we can find a pseudoinverse of a matrix \(A\text{.}\) A singular value decomposition provides an answer to this question. If \(A\) is an invertible \(n \times n\) matrix, then 0 is not an eigenvalue of \(A\text{.}\) As a result, in the singular value decomposition \(U \Sigma V^{\tr}\) of \(A\text{,}\) the matrix \(\Sigma\) is an invertible matrix (note that \(U\text{,}\) \(\Sigma\text{,}\) and \(V\) are all \(n \times n\) matrices in this case). So

where

In this case, \(V \Sigma^{-1} U^{\tr}\) is a singular value decomposition for \(A^{-1}\text{.}\)

To understand in general how a pseudoinverse is found, let \(A\) be an \(m \times n\) matrix with \(m \neq n\text{,}\) or an \(n \times n\) with rank less than \(n\text{.}\) In these cases \(A\) does not have an inverse. But as in Preview Activity 30.7, a singular value decomposition provides a pseudoinverse \(A^+\) for \(A\text{.}\) Let \(U \Sigma V^{\tr}\) be a singular value decomposition of an \(m \times n\) matrix \(A\) of rank \(r\text{,}\) with

The matrices \(U\) and \(V\) are invertible, but the matrix \(\Sigma\) is not if \(A\) is not invertible. If we let \(\Sigma^+\) be the \(n \times m\) matrix defined by

then \(\Sigma^{+}\) will act much like an inverse of \(\Sigma\) might. In fact, it is not difficult to see that

where \(\Sigma\Sigma^{+}\) is an \(m \times m\) matrix and \(\Sigma^{+}\Sigma\) is an \(n \times n\) matrix.

The matrix

is a pseudoinverse of \(A\text{.}\)

Pseudoinverses satisfy several properties that are similar to those of inverses. For example, we had an example in Preview Activity 30.7 where \(AA^{+}A = A\) and \(A^+AA^+ = A^+\text{.}\) That \(A^+\) always satisfies these properties is the subject of the next activity.

Activity 30.9 shows that \(A^{+}\) satisfies properties that are similar to those of an inverse of \(A\text{.}\) In fact, \(A^{+}\) satisfies several other properties (that together can be used as defining properties) as stated in the next theorem. The conditions of Theorem 30.6 are called the Penrose or Moore-Penrose conditions.⁵³ Verification of the remaining parts of this theorem are left for the exercises.

Also, there is a unique matrix \(A^+\) that satisfies these properties. The verification of this property is left to the exercises.

Subsection Least Squares Approximations

The pseudoinverse of a matrix is also connected to least squares solutions of linear systems as we encountered in Section 24. Recall from Section 24, that if the columns of \(A\) are linearly independent, then the least squares solution to \(A\vx = \vb\) is \(\vx = \left(A^{\tr}A\right)^{-1}A^{\tr} \vb\text{.}\) In this section we will see how to use a pseudoinverse to solve a least squares problem, and verify that if the columns of \(A\) are linearly dependent, then \(\left(A^{\tr}A\right)^{-1}A^{\tr}\) is in fact the pseudoinverse of \(A\text{.}\)

Let \(U \Sigma V^{\tr}\) be a singular value decomposition for an \(m \times n\) matrix \(A\) of rank \(r\text{.}\) Then the columns of

form an orthonormal basis for \(\R^m\) and \(\{\vu_1, \vu_2, \ldots, \vu_r\}\) is a basis for \(\Col A\text{.}\) Remember from Section 25 that if \(\vb\) is any vector in \(\R^m\text{,}\) then

is the least squares approximation of the vector \(\vb\) by a vector in \(\Col A\text{.}\) We can extend this sum to all of columns of \(U\) as

It follows that

where

Now, if \(\vz = A^{+} \vb\text{,}\) then

and hence the vector \(A \vz = AA^+ \vb\) is the vector \(A\vx\) in \(\Col A\) that minimizes \(||A \vx - \vb||\text{.}\) Thus, \(A\vz\) is in actuality the least squares approximation to \(\vb\text{.}\) So a singular value decomposition allows us to construct the pseudoinverse of a matrix \(A\) and then directly solve the least squares problem.

If the columns of \(A\) are linearly independent, then we do not need to use an SVD to find the pseudoinverse, as the next activity illustrates.

Therefore, to calculate \(A^{+}\) and solve a least squares problem, Theorem 30.7 shows that as long as the columns of \(A\) are linearly independent, we can avoid using a singular value decomposition of \(A\) in finding \(A^{+}\text{.}\)

Subsection Examples

What follows are worked examples that use the concepts from this section.

Subsection Summary

• The condition number of an \(m \times n\) matrix \(A\) is the number \(||A^{-1}|| \ ||A||\text{.}\) The condition number provides a measure of how well the relative error in a calculated value \(\Delta \vb\) predicts the relative error in \(\Delta \vx\) when we are trying to solve a system \(A \vx = \vb\text{.}\)

• A pseudoinverse \(A^{+}\) of a matrix \(A\) can be found through a singular value decomposition. Let \(U \Sigma V^{\tr}\) be a singular value decomposition of an \(m \times n\) matrix \(A\) of rank \(r\text{,}\) with

  \begin{equation*} \Sigma = \left[ \begin{array}{ccccc|c} \sigma_1\amp \amp \amp \amp \amp \\ \amp \sigma_2\amp \amp 0\amp \amp \\ \amp \amp \sigma_3\amp \amp \amp 0 \\ \amp 0 \amp \amp \ddots \amp \amp \\ \amp \amp \amp \amp \sigma_r \\ \hline \amp \amp 0\amp \amp \amp 0 \end{array} \right] \end{equation*}

  If \(\Sigma^+\) is the \(n \times m\) matrix defined by

  \begin{equation*} \Sigma^{+} = \left[ \begin{array}{ccccc|c} \frac{1}{\sigma_1}\amp \amp \amp \amp \amp \\ \amp \frac{1}{\sigma_2}\amp \amp 0\amp \amp \\ \amp \amp \frac{1}{\sigma_3}\amp \amp \amp 0 \\ \amp 0 \amp \amp \ddots \amp \amp \\ \amp \amp \amp \amp \frac{1}{\sigma_r} \\ \hline \amp \amp 0\amp \amp \amp 0 \end{array} \right]\text{,} \end{equation*}

  then \(A^{+} = V\Sigma^{+}U^{\tr}\text{.}\)

• A pseudoinverse \(A^{+}\) of a matrix \(A\) acts like an inverse for \(A\text{.}\) So if we can't solve a matrix equation \(A \vx = \vb\) because \(\vb\) isn't in \(\Col A\text{,}\) we can use the pseudoinverse of \(A\) to “solve” the equation \(A \vx = \vb\) with the “solution” \(A^+ \vb\text{.}\) While not an exact solution, \(A^+ \vb\) turns out to be the best approximation to a solution in the least squares sense.

Subsection Project: GPS and Least Squares

In this project we discuss some of the details about how the GPS works. The idea is based on intersections of spheres. To build a basic understanding of the system, we begin with a 2-dimensional example.

Project Activity 30.11.

Suppose that there are three base stations \(A\text{,}\) \(B\text{,}\) and \(C\) in \(\R^2\) that can send and receive signals from your mobile phone. Assume that \(A\) is located at point \((-1,-2)\text{,}\) \(B\) at point \((36,5)\text{,}\) and \(C\) at point \((16,35)\text{.}\) Also assume that your mobile phone location is point \((x,y)\text{.}\) Based on the time that it takes to receive the signals from the three base stations, it can be determined that your distance to base station \(A\) is \(28\) km, your distance to base station \(B\) is \(26\) km, and your distance to base station \(C\) is \(14\) km using a coordinate system with measurements in kilometers based on a reference point chosen to be \((0,0)\text{.}\) Due to limitations on the measurement equipment, these measurements all contain some unknown error which we will denote as \(z\text{.}\) The goal is to determine your location in \(\R^2\) based on this information.

If the distance readings were accurate, then the point \((x,y)\) would lie on the circle centered at \(A\) of radius \(29\text{.}\) The distance from \((x,y)\) to base station \(A\) can be represented in two different ways: \(28\) km and \(\sqrt{(x+1)^2 + (y+2)^2}\text{.}\) However, there is some error in the measurements (due to the receiver clock and satellite clocks not being snychronized), so we really have

\begin{equation*} \sqrt{(x+1)^2 + (y+2)^2} + z = 28\text{,} \end{equation*}

where \(z\) is the error. Similarly, \((x,y)\) must also satisfy

\begin{equation*} \sqrt{(x-36)^2 + (y-5)^2} + z = 26 \end{equation*}

and

\begin{equation*} \sqrt{(x-16)^2 + (y-35)^2} + z = 14\text{.} \end{equation*}

(a)

Explain how these three equations can be written in the equivalent form

\begin{align} (x+1)^2+(y+2)^2\amp =(28-z)^2\tag{30.11}\\ (x-36)^2 + (y-5)^2 \amp = (26-z)^2\tag{30.12}\\ (x-16)^2 + (y-35)^2 \amp = (14-z)^2\text{.}\tag{30.13} \end{align}

(b)

If all measurements were accurate, your position would be at the intersection of the circles centered at \(A\) with radius \(28\) km, centered at \(B\) with radius \(26\) km, and centered at \(C\) with radius \(14\) km as shown in Figure 30.12. Even though the figure might seem to imply it, because of the error in the measurements the three circles do not intersect in one point. So instead, we want to find the best estimate of a point of intersection that we can. The system of equations (30.11), (30.12), and (30.13) is non-linear and can be difficult to solve, if it even has a solution. To approximate a solution, we can linearize the system. To do this, show that if we subtract corresponding sides of equation (30.11) from (30.12) and expand both sides, we can obtain the linear equation

\begin{equation*} 37x + 7y + 2z = 712 \end{equation*}

in the unknowns \(x\text{,}\) \(y\text{,}\) and \(z\text{.}\)

(c)

Repeat the process in part (b), subtracting (30.11) from (30.13) and show that we can obtain the linear equation

\begin{equation*} 17x + 37y + 14z = 1032 \end{equation*}

in \(x\text{,}\) \(y\text{,}\) and \(z\text{.}\)

(d)

We have reduced our system of three non-linear equations to the system

\begin{alignat*}{4} {37}x \amp {}+{} \amp {7}y \amp {}+{} \amp {2}z \amp = \amp {} \amp 712\\ {17}x \amp {}+{} \amp {37}y \amp {}+{} \amp {14}z \amp = \amp {} \amp 1032 \end{alignat*}

of two linear equations in the unknowns \(x\text{,}\) \(y\text{,}\) and \(z\text{.}\) Use technology to find a pseudoinverse of the coefficient matrix of this system. Use the pseudoinverse to find the least squares solution to this system. Does your solution correspond to an approximate point of intersection of the three circles?

Project Activity 30.11 provides the basic idea behind GPS. Suppose you receive a signal from a GPS satellite. The transmission from satellite \(i\) provides four pieces of information — a location \((x_i,y_i,z_i)\) and a time stamp \(t_i\) according to the satellite's atomic clock. The time stamp allows the calculation of the distance between you and the \(i\)th satellite. The transmission travel time is calculated by subtracting the current time on the GPS receiver from the satellite's time stamp. Distance is then found by multiplying the transmission travel time by the rate, which is the speed of light \(c=299792.458\) km/s.⁵⁵ So distance is found as \(c(t_i-d)\text{,}\) where \(d\) is the time at the receiver. This signal places your location within in a sphere of that radius from the center of the satellite. If you receive a signal at the same time from two satellites, then your position is at the intersection of two spheres. As can be seen at left in Figure 30.13, that intersection is a circle. So your position has been narrowed quite a bit. Now if you receive simultaneous signals from three spheres, your position is narrowed to the intersection of three spheres, or two points as shown at right in Figure 30.13. So if we could receive perfect information from three satellites, then your location would be exactly determined.

There is a problem with the above analysis — calculating the distances. These distances are determined by the time it takes for the signal to travel from the satellite to the GPS receiver. The times are measured by the clocks in the satellites and the clocks in the receivers. Since the GPS receiver clock is unlikely to be perfectly synchronized with the satellite clock, the distance calculations are not perfect. In addition, the rate at which the signal travels can change as the signal moves through the ionosphere and the troposphere. As a result, the calculated distance measurements are not exact, and are referred to as pseudoranges. In our calculations we need to factor in the error related to the time discrepancy and other factors. We will incorporate these errors into our measure of \(d\) and treat \(d\) as an unknown. (Of course, this is all more complicated that is presented here, but this provides the general idea.)

To ensure accuracy, the GPS uses signals from four satellites. Assume a satellite is positioned at point \((x_1,y_1,z_1)\) at a distance \(d_1\) from the GPS receiver located at point \((x,y,z)\text{.}\) The distance can also be measured in two ways: as

and as \(c(t_1-d)\text{.}\) So

Again, we are treating \(d\) as an unknown, so this equation has the four unknowns \(x\text{,}\) \(y\text{,}\) \(z\text{,}\) and \(d\text{.}\) Using signals from four satellites produces the system of equations

We conclude this project with a final note. At times a GPS receiver may only be able to receive signals from three satellites. In these situations, the receiver can substitute the surface of the Earth as a fourth sphere and continue the computation.