Mathematical and NISTIR 6131

Computational

Sciences

Division

Information Technology Laboratory

TRUNCATING THE SINGULAR

VALUE DECOMPOSITION FOR

ILL-POSED PROBLEMS

Bert W. Rust

July 1998

U.S. DEPARTMENT OF COMMERCE

National Institute of Standards and Technology

Gaithersburg, MD 20899
TRUNCATING THE SINGULAR VALUE DECOMPOSITION

FOR ILL-POSED PROBLEMS

Bert W. Rust

Mathematical and Computational Sciences Division

Building 820, Room 365

National Institute of Standards and Technology

Gaithersburg, MD 20899

Abstract

Discretizing the rst-kind integral equations which model many physical measurement processes yields an ill-conditioned linear regression model

b = Ax + , where xis a vector representation of the function being measured, A is an instrument response matrix, b is a vector of measurements, and

 is a vector of unknown, random measuring errors. Least squares estimation

usually gives a sum of squared residuals much smaller than the expected value

and a wildly oscillating, physically implausible estimate of x. These symptoms suggest that the least squares estimate captures part of the variance that

properly belongs in the residuals. One strategy for shifting some of this variance to the residuals and simultaneously stabilizing the estimate is to truncate

the singular value decomposition A = UVTwhere U and V are orthogonal

matrices and  is a diagonal matrix of singular values. All of the singular values below some threshold value are reset to zero to give a new matrix tr , and

the estimated solution is calculated from the generalized inverse of the matrix

UtrVT. The most delicate part of this procedure is the determination of the

truncation threshold. Conventionally this has been regarded as a problem of

determining the \numerical rank" of A, but in most cases A is clearly not

rank-decient. This paper suggests an alternate strategy which uses the variances of the measuring errors to specify a truncation for the elements of the

rotated measurement vector UTb. The idea is to zero all of the components

that are dominated by the measurement errors and compute the estimate using

the full rank matrix. The problem of setting the truncation threshold becomes

one of deciding whether or not a measured value is signicantly dierent from

zero, a procedure familiar to most experimentalists. The paper also develops

some new diagnostics for the residuals which are useful not only for choosing

the truncation level for the (UTb)i , but also for assessing the quality of an

estimate obtained by any procedure.
1 First Kind Integral Equations with Uncertainties

First kind integral equations,

y(t) =

Z

b

a

K(t; )x() d ; (1.1)

where K(t; ) and y(t) are known functions, are routinely used to model instrument distortions in measuring an unknown function x(). In that context,

they are usually written as a system of equations

yi =

Z

b

a

Ki()x() d + i ; i = 1; 2; : : : ; m ; (1.2)

where the Ki() are known (previously measured or calculated) response functions of the instrument, the yi = y(ti) are measured values, corresponding to

a discrete mesh t1; t2; : : : ; tm, and the i are random, zero-mean measuring errors. In order to estimate x() it is necessary to further discretize the system,

in the process replacing it with a linear regression model

y = Kx+  ; (1.3)

where y is the m-vector of measurements, K is a known m n matrix, with

m  n, and xis an unknown n-vector whose components are either discrete

point estimates of x() on some mesh 1; 2; : : : ; n, or are the unknown coe-

cients in a truncated expansion of x() in terms of some set of basis functions

spanning the space of possible solutions. The vector  is an m-vector of random

measuring errors satisfying

E () = 0 ; E



 

T



= S2; (1.4)

where E is the expectation operator, 0 is the m-dimensional zero vector and

S2is the positive denite variance-covariance matrix for . In most problems

the measurement errors are statistically independent so

S2= diag(s

2

1; s

2

2; : : : ; s

2

m) ; (1.5)

where s1; s2; : : : ; sm are the standard deviations of the errors. If S2is not

diagonal, the model can be transformed into one with a diagonal variance

matrix by premultiplying (1.3) by the inverse of the lower triangular Cholesky

factor of S2.

1
There seems to be a general misapprehension in the numerical analysis

community that estimates of the si are seldom available, but in fact good

experimenters routinely provide them. Estimates of the measurement errors

are not considered to be something extra, but rather are an integral part of

the measurements, and published graphs of measured data will usually report

them as 1 error bars on the plotted points. An analyst who fails to use

this information implicitly assumes that S2 = s

2

Im where Im is the m-th order

identity matrix and s is an unknown scalar that can be, but usually is not, estimated from the sum of squared residuals for the least squares solution. Such

an assumption may be tenable if all the yi have roughly the same magnitude,

but in most cases, they span a range of values, and the magnitudes of the si

vary (usually nonlinearly) with the magnitudes of the yi.

In the following it will be assumed that S is a known diagonal matrix. It

will also be assumed, as is often the case, that the errors are samples from

a multivariate normal distribution, i.e., that   N( 0 ; S2). In x2, these

assumptions will be used to rescale the problem to have errors with a standard normal distribution and to derive a statistical diagnostic for the sum of

squared residuals for any estimate of the solution vector. In x3 a variant of the

well known Phillips problem [16] is subjected to this scaling and the ordinary

least squares estimate is calculated and found to be unsatisfactory. In x4 the

residuals for that estimate are subjected to diagnostic tests based on the expected value of their sum of squares and on their distribution when considered

as a time-series. The diagnostics that are developed there are useful for testing

any estimate, no matter how it may be obtained, and will be used throughout

the remainder of the paper.

The singular value decomposition is introduced in x5 and the conventional

method for stabilizing solution estimates by truncating the distribution of singular values is described in x6. This method is based on an assumption that

the matrix is rank decient and depends critically on the proper determination of its \numerical rank." But for most ill-posed problems, the matrix is

demonstrably not rank decient so there is no good reason for setting any of

the singular values to zero. An alternate approach to truncating the decomposition is given in x7. Rather than zeroing some of the singular values, one

instead zeroes small elements of the rotated right hand side vector formed by

premultiplying the vector of measurements by the transpose of the matrix of

left singular vectors, i.e., the leftmost factor in the singular value decomposition. The errors in this rotated measurement vector have a standard normal

distribution, so it is possible to establish a statistical criterion for determining whether or not a given element is signicantly dierent from zero. The

new method is applied to the modied Phillips problem to obtain satisfactory

2
estimates for 3 dierent truncation levels and the optimum level is taken to

be the one which best satises the diagnostic tests described in x4. In x8 the

method is successfully applied to real-world measurements of the energy spectrum of neutrons produced by a certain nuclear reaction. Finally, x9 gives a

brief discussion of algorithmic considerations and of how the method can be

extended when the knowledge of the measurement errors is not as complete as

might be desired.

2 Scaling the Problem

The linear regression model in the preceding section can be written

y = Kx+  ;   N( 0 ; S2) ; (2.1)

but it is advantageous to scale it with the matrix S1. Let

b  S1y ; A  S1K ;   S1 ; (2.2)

and note that by a standard theorem of multivariate statistics [2, Thm. 2.4.4],

  N( 0 ; S2) implies that   N( S10 ; S1S2[S1]

T

), so the scaled

model can be written

b = Ax+  ;   N( 0 ; Im ) ; (2.3)

or

b  N( Ax; Im ) : (2.4)

To see the advantage of this scaling, let x^ be an estimate of xand

^r = b  Ax^ ; (2.5)

be the corresponding residual vector. Since the regression model can also be

written

 = b  Ax; (2.6)

it is clear that an estimate x^ is acceptable only if ^r is a plausible sample from

the -distribution.

Since b  Ax  N( 0 ; Im ), it follows from another standard statistical

theorem [14, page 140] that

kb  Axk

2

 (b  Ax)

T

(b  Ax)  2(m) ; (2.7)

3
where 2

(m) denotes the Chi-squared distribution with m degrees of freedom,

and hence that

E fkb  Axk2g = m ; Var fkb  Axk2g = 2m : (2.8)

These two quantities provide rough bounds for the sum of squared residuals

that might be expected from a reasonable estimate of x. An estimate that

gives

m 

p

2m  kb  Ax^k

2

 m +

p

2m (2.9)

would be quite reasonable, but any x^ whose sum of squared residuals falls

outside the interval hm  2

p

2m; m + 2

p

2m

i

would be suspect. These rough

indicators can be sharpened considerably by using the cumulative distribution

function for 2(m). More details on this point are given in Appendix A.

3 Standard Linear Regression

The standard approach for the linear regression problem dened by [2.3] is to

assume that rank(A) = n and seek the minimum variance, unbiased estimator

for xby solving the least squares minimization problem

r

2

min = min

x2Rn

n

(b  Ax)

T

(b  Ax)

o

: (3.1)

The solution,

x^ =



ATA

1

ATb ; (3.2)

is called the best linear unbiased estimate of x, but it is well known [18,

Chapt. 1], [19, Chapt. 6], [21, Chapt. 2] that for regression models obtained

by discretizing rst kind integral equations, the elements of x^ are pathologically sensitive to small variations in the elements of b, so the presence of

the measuring errors leads to estimates that are totally unphysical, typically

oscillating wildly between extreme positive and negative values.

A useful test problem, with known solution, which shares many of the characteristics of real instrument correction problems can be obtained by discretizing a variant of the well known Phillips equation [16] and adding some random

errors to the discrete yi. The problem and a discretization with m = 300 and

n = 241 are described in detail in Appendices B and C. The functions y(t)

and x() are plotted on the right in Figure 1, and on the left are plotted the

functions K(t; j ) for the discrete values j = 3:0 ; 1:5 ; 0 ; 1:5 ; 3:0 . All of

the 241 K(t; k) have the same shape and subtend unit area. The 300 standard

deviations si range in value from 3:49 1011 to 7:78 106so the errors in

4
Figure 1: A variant of the Phillips problem

the yi are much smaller than the thickness of the curve in the plot of y(t).

But these small errors produce large oscillations in the least squares estimate

of the solution vector. This is shown in the plot in the upper left corner of

Figure 2 where the barely discernable dashed curve is the true solution x, and

the solid curve oscillating wildly about it is the least squares estimate x^. The

magnitudes of these oscillation are roughly 107times greater than the largest

random errors in the yi. Such a large amplication of the measuring errors is

the rule rather than the exception for ill posed problems.

4 Estimate Diagnostics

The most commonly used estimate diagnostic for linear regression problems

is the sum of squared residuals. For the problem in the preceding section the

least squares estimate gives

r

2

min = min

x2Rn

kb  Axk

2

= kb  Ax^k

2

= 43:01 (4.1)

which is much smaller than the value expected for the true solution which, by

Eq. (2.8), is

E

n

kb  Axk

2

o

= 300 : (4.2)

5
Figure 2: Least squares estimate and diagnostics for the modied Phillips

problem

From Eq. (2.8) it also follows that the standard deviation for this expected

value is p

600 = 24:49, so the least squares r

2

min is more than 10 standard

deviations smaller than the expected value. Thus, without even inspecting the

estimate itself, one could argue that x^ is unlikely to be a good approximation

to x.

Another good diagnostic for judging the acceptability of an estimate can be

gotten by considering the elements of the vectors  and ^r as time series, with

the element number taken as the time variable. Since the i are statistically

independent and identically normally distributed, they form a normally distributed white noise series, so the residuals for an acceptable estimate should

constitute a realization of such a series. The residuals for the least squares

estimate are plotted against element number at the upper right of Figure 2. It

is clear by inspection that they are not even approximately white noise, and

this is not surprising since it can be shown [10, Section 6.9.1] that

^r  N



0 ;



Im  A



ATA

1

AT



(4.3)

6
with the variance matrix becoming approximately diagonal only when m is

much larger than n.

A rigorous test for white noise is based on the periodogram [6, Chapt. 7]

which is an estimate of the power spectrum on the frequency interval 0 

f  1

2T

, where T is the sample spacing for the time variable. Here, the time

variable is the element number i, so T = 1. The graph at the lower left of

Figure 2 is a plot of a periodogram estimate at 4097 equally spaced frequency

points on the interval [0; 0:5]. The details on how it was computed are given

in Appendix D. For a white noise series, the variance would be distributed

uniformly over the frequency range, so the fact that there is much more power

on the interval [0; 0:25] than on [0:25; 0:5] suggests that the residuals should

not be regarded as white noise. A formal test of this hypothesis uses the

cumulative periodogram [6, Chapt. 7] (also described in Appendix D) which

is plotted as a solid curve at the lower right of the gure. The legend box

gives the length of this curve which can be compared with the value 1:11803

expected for a pure white noise spectrum. The two dashed lines enclose a 95%

condence band for white noise. The cumulative periodogram ordinates for

such a series should lie outside this band for at most 5% of the frequencies.

Since 2803 of the 4096 calculated ordinates lie outside the band, it is safe to

conclude that the least squares residuals are not white noise.

5 The Singular Value Decomposition

Insight into the failure of the least squares method to give a reasonable estimate

of xcan be gotten by substituting the singular value decompostion of A into

the model (2.3). Let

A = U

 



O

!

VT;  = diag(1; 2; : : : ; n) ; (5.1)

where 1  2      n, and

UTU = Im = UUT; VTV = In = VVT: (5.2)

Substituting (5.1) into (2.3) and premultiplying by UTgives

UTb =

 



O

!

VTx+ UT ; UT  N( 0 ; Im ) ; (5.3)

with the distribution of the UT vectors unchanged because premultiplication

by an orthogonal matrix simply rotates all the vectors in the distribution

7
through the same angle. Also, premultiplying a vector by an orthogonal matrix

leaves its length unchanged so Eq. (2.7) can be replaced by

kb  Axk

2

=
















UTb 

 



O

!

VTx
















2

 2(m) ; (5.4)

with the 2(m) distribution following from the fact that the (UT)i  n(0; 1),

idependently. Partitioning the m m matrix U,

U = (U1; U2) ; (5.5)

with U1 an m n and U2 an m (m  n) submatrix allows the least squares

problem to be written

min

x2Rn

kb  Axk

2

= min

x2Rn
















 

UT

1

UT

2

!

b 

 



O

!

VTx
















2

; (5.6)

from which it is quite evident that the least squares solution must satisfy

VTx^ = 1UT

1 b ; (5.7)

and that the minimum sum of squared residuals is

r

2

min = kb  Ax^k

2

= kUT

2 bk

2

: (5.8)

These last two equations can also be written



VTx^



i =



UTb



i

i

; i = 1; 2;    ; n ; (5.9)

and

r

2

min = kb  Ax^k

2

= Xm

i=n+1



UTb

2

i

; (5.10)

and it is easy to see that Eq. (4.3) can be written

^r  N



0 ; Im  U1UT

1



: (5.11)

For the modied Phillips problem the value of r

2

min was too small by 10 standard deviations which suggests that some of the (UT

b)i values in the sequence

(5.9) more properly belong in the sum in (5.10). This reasoning leads to the

idea of truncating the singular value decomposition.

8
6 Truncating the Singular Value Distribution

Historically, the idea of truncating the SVD has been treated as a problem

of determining the \numerical rank" of the matrix A. All of the computed

singular values smaller than some threshold value are treated as zeroes which

were corrupted by rounding errors into small non-zero quantities. Thus if p

is the smallest singular value greater than the truncation threshold, then one

replaces the matrix  by a truncated matrix

tr = diag(1; : : : ; p; 0; : : : ; 0) (6.1)

whose generalized inverse is given by

y

tr = diag 

1

1

; : : : ;

1

p

; 0; : : : ; 0

!

: (6.2)

The value p is said to be the \numerical rank" of A. If x~ is the estimate of

the solution to the truncated problem, then it is required to satisfy

VTx~ = y

trUT

1 b (6.3)

rather than Eq. (5.7), which means that Eqs. (5.9) and (5.10) are replaced by



VTx~



i =

8

>>><

>>>:



UTb



i

i

; i = 1; 2;    ; p ;

0 ; i = p + 1;    ; n ;

(6.4)

and

kb  Ax~k

2

= Xm

i=p+1



UTb

2

i

; (6.5)

and Eq. (5.11) is replaced by

~r  N

hIm  UpUT

p



Ax; Im  UpUT

p

i

; (6.6)

where Up is the submatrix formed by the rst p columns of U. Since p < n,

it is clear that the variance matrix for ~r is more nearly diagonal than that

for ^r, and that if p  m, then the r~i will be approximately independently

distributed.

The truncated singular value approach was rst suggested by Golub and

Kahan [8] who noted its similarity to the theoretical treatment given by

Smithies [17, Chapt. 8] for the singular functions and singular values of rst

9
kind integral equations. One of the rst to use the method was Richard Hanson [13] who suggested that the truncation threshold should be chosen to be

the smallest integer p such that

Xm

i=p+1



UTb

2

i

< m : (6.7)

In view of Eqs. (2.7) and (2.8), this seems a very sensible choice, but his suggestion was apparently ignored by most other workers in the eld who were

preoccupied with the problem of rigorously dening the notion of the \numerical rank" of the matrix. An especially in
uential paper in this line of

research was a technical report [9] by Golub, Klema and Stewart who used

techniques from perturbation analysis to dene a very complicated criterion

for rank determination which was norm-dependent and involved a triplet of

numbers (; ; p) with  and  being numbers satisfying p   >   p+1.

The idea was to nd a clear gap in the distribution of singular values and zero

all those on the low side of the gap. The trouble with this approach is that

matrices arising from discretizing rst kind integral equations almost never

display such a gap. This fact was clearly enunciated by P.C. Hansen [11] who

proposed \... a natural division of ill-conditioned matrices into two classes:

those with well-determined numerical rank and those with ill-determined numerical rank." By the latter he meant matrices whose singular values decayed

smoothly from 1 to n with no obvious gap between the larger and smaller

ones. More recently, in his new book [12, Chapt. 4], he has written \The main

feature of these problems is that all the singular values of the coecient matrix decay gradually to zero, with no gap anywhere in the spectrum. ... and

therefore the concept of `numerical rank' is not useful for these problems."

Nevertheless he devoted large segments of the book to problems of rank determination and rank revealing decompositions, and gave a good review of the

very sizable literature on these subjects.

The singular values for the modied Phillips problem are plotted as discrete squares in Figure 3 which also shows the elements of the rotated right

hand side vector UTb plotted as circles connected by straight line segments.

The decomposition was computed by subroutine SGESVD from the LAPACK

collection [1], using double precision arithmetic on a Sparc-20 workstation.

The largest and smallest singular values are

1(A) = 2:882 105; 241(A) = 2:622 102; (6.8)

which gives cond(A) = 1:099 107. The relative acccracy of the arithmetic

was mach = 2:22 1016, so the smallest singular value is many orders of

10
Figure 3: Singular values (squares) and rst n elements of jUT

bj (circles) for

the modied Phillips problem

magnitudes larger than the eective zero level. Thus there is no reason to

assume that rank(A) < n. Yet the problem obviously needs some sort of truncation to prevent the estimated solution from capturing variance that properly

belongs in the residuals. Fortunately, the desired result can be obtained by

leaving the singular values unchanged and truncating the rotated right hand

side vector UTb.

7 Truncating the Vector UTb

Even a cursory inspection of the jUTbj

i plotted in Figure 3 reveals a sharp

dichotomy in the distribution at i = 55. Before the break, the upper envelope

of the distribution is decreasing at a faster rate than that of the singular values.

After the break, the distribution is quite 
at with 64 of the 186 values above

the horizontal line jUTbj = 1 and 122 values below. That line is, by Equations

(5.3), just the one standard deviation level for the random errors (UT)i. If the

nal 186 elements of UTb were chosen randomly from an n(0; 1) distribution,

11
Figure 4: Singular values (squares) and rst n elements of jUT

bj (circles) for

the modied Phillips problem with no random errors

then one would expect roughly 62 of them to satisfy jUTbj

i > 1 and the other

124 to satisfy jUTbj

i  1. Since that is almost exactly what is observed, it is

quite reasonable to conclude that the 
at part of the distribution represents

elements of UTb that are dominated by the measuring errors. This conclusion

is reinforced by Figure 4 which gives plots of the i(A) and jUTbj

i for the

same problem with no random errors in the b vector. The singular values are

identical to the ones in Figure 3, but the distribution of the jUTbj

i satises the

discrete Picard condition over the whole range, and Eq. (5.9) yields estimates

of the elements of xthat are correct to 8 signicant digits.

The foregoing observations suggest a more logical way to truncate the

singular value decomposition. Rather than zeroing some of the singular values

one should instead zero those elements of jUTbj that are judged to be mostly

random error. The idea is to pick a truncation level  and require the estimated

12
solution x~ to satisfy

(VTx~)i =

8

>>><

>>>:



UTb



i

i

; if jUTbj

i >  ;

0 ; otherwise ;

i = 1; 2; : : : ; n : (7.1)

The sum of squared residuals is then given by

r

2

= kb  Ax~k

2

= X

i 2 I



UTb

2

i

; (7.2)

where the indexing set for the sum is

I =

n

i





 jUTbj

i   ; i = 1; 2; : : : ; n

o [

f n + 1; n + 2; : : : ; m g : (7.3)

If ui is the ith column of U, then Eq. (7.1) can be written

(VTx~)i =

8

>><

>>:



uT

i b



i

; if juT

i bj >  ;

0 ; otherwise ;

i = 1; 2; : : : ; n ; (7.4)

and if an m n matrix U~

1 is dened by

U~

1 = (u~1; u~2;    ; u~n) ; where u~i =

(

ui ; if juT

i bj >  ;

0 ; otherwise ;

(7.5)

then the estimated solution and residual vectors can be written

x~ = V1U~ T

1 b ; and ~r =



Im  U~

1U~ T1



b : (7.6)

It follows from a fundamental theorem in multivariate analysis [2, Thm. 2.4.5]

that

~r  N

hIm  U~

1U~ T1



Ax; Im  U~

1U~ T

1

i

: (7.7)

It is clear that as  increases, the number of non-zero columns in U~ decreases,

so the variance matrix for ~r more nearly approximates a diagonal matrix.

The success of the proposed method depends crucially on the choice of the

truncation level  . A safe and eective approach is to try several values of 

and use the diagnostics described in Section 4 and Appendix D to make the

nal choice. The results of this strategy for the modied Phillips problem are

summarized in Table 1. The labels for the rst, second and nal columns are

13
Table 1: Estimate diagnostics for the modied Phillips Problem



P

(b  Ax~)

2

i %(Ck 2 B95) LengthfCg jxjrms Figure

0.0 43.0 31.6 1.238 14.937 2

3.0 235.6 79.5 1.211 5.6865 5

4.0 267.7 100.0 1.205 0.0106 6

4.7 289.2 100.0 1.209 0.0122

5.0 311.4 100.0 1.207 0.0132

9.0 380.7 100.0 1.206 0.0148

self explanatory. The third column contains the percentage of the 4096 cumulative periodogram ordinates that lie inside the 95% condence band for white

noise. The fourth column contains the length of the cumulative periodogram

plot which can be compared with the length 1.11803 for a pure white noise

series. The fth column gives the root mean square average magnitude of the

errors in the elements of the estimated solution vector, i.e.,

jxjrms =

v

u

u

t

1

n

Xn

j=1

jx~j  x

j

j

2

: (7.8)

It can be included only because the true solution is known. The rst row in

the table gives the results for the least squares estimate shown in Figure 2.

A guideline for choosing a lower bound for the value of  is the fact that

most experimentalists would be reluctant to claim that a measured value is

signicantly dierent from zero if its magnitude does not exceed 3 standard

deviations for the error in the measurement. Since each jUTbj

i is scaled in

units of one standard deviation of its own random error, it suces to choose

 = 3:0. The results obtained for this choice are given in Figure 5 and row 2

of the table. The sum of squared residuals, 235.6, is less than the 1 percentile

value for the 2(300) distribution (see column 6 of Table 3, in Appendix A),

i.e.,

Pr nkb  Axk

2

 235:6

o

< Pr nkb  Axk

2

 245:3

o

= 0:01 : (7.9)

Furthermore, only 79:5% of the cumulative periodogram ordinates lie inside

the 95% band for white noise. The plot of the estimate reveals the persistence

of strong high frequency oscillations which, by inspection of Figure 3, can be

attributed to the component jUTbj205 = 3:547. Repeatedly measuring a standard normally distributed quantity will give a measurement with magnitude

14
Figure 5: Truncated solution and diagnostic plots for the modied Phillips

problem with  = 3

greater than 3 roughly once in every 385 measurements so the occurrence of

one such value in the nal 186 jUTbj

i is not an unlikely event. Its eect can

be transferred from the estimated solution to the residuals by increasing the

truncation threshold to  = 4. The results for that value are given in Figure 6

and row 3 of the table. The sum of squared residuals, 267.7, diers from the

expected value by about 1.3 standard deviations. While this is not a small

dierence, it is clear from column 6 of Table 3 that

0:05 = Pr nkb  Axk

2

 260:6

o

< Pr nkb  Axk

2

 267:7

o

; (7.10)

so 267.7 must be regarded as an acceptable value. This conclusion is reinforced

by the fact that none of the cumulative periodogram ordinates lie outside the

95% band. The plot of the estimated solution, at the upper left of Fig. 6, is

so close to that of the true solution that the two are almost indistinguishable.

Since there are only 241 discrete jUTbj

i, there are only 241 possible truncated estimates x~. Only 55 of them correspond to truncation levels greater

than  = 4. The smallest four of those are

jUTbj56 = 4:63; jUTbj51 = 4:71; jUTbj54 = 8:33; jUTbj52 = 9:49: (7.11)

15
Figure 6: Truncated solution and diagnostic plots for the modied Phillips

problem with  = 4

These values dene 3 ranges for  , in each of which the estimates are all

identical. For the lowest 2 ranges, the results for  = 4:7 and  = 5:0 are given

in the fourth and fth rows of Table 1. In both cases, the estimated sum of

squared residuals lie within one standard deviation of the expected value, and

all of the cumulative periodogram ordinates lie in the white noise band. The

highest range is represented by  = 9 with the results summarized in the row

6 of the table. Although all of the periodogram ordinates lie within the 95%

condence band for white noise, the sum of squared residuals, 380.7, is much

too large to be acceptable since, from column 6 of Table 3,

Pr fkb  Axk2 > 380g < Pr fkb  Axk2 > 359:1g

= 1  Pr fkb  Axk2  359:1g

= 0:01 :

(7.12)

Thus, the only estimates to give acceptable values for the sum of squared

residuals are the ones corresponding to truncation levels  = 4:0; 4:7; and 5:0.

Since the one for  = 4:0 gave the smallest LengthfCg, it was chosen as the

most acceptable, and this choice is validated by the fact that it also gave the

16
smallest jxjrms. But the other two estimates are quite good also, and are

qraphically almost indistinguishable from the  = 4:0 estimate.

8 A Real-World Example

The measurement of nuclear radiation spectra is one important source of rst

kind integral equations with uncertainties. Consider the energy spectrum of

the neutrons produced by the nuclear reaction

T ( d ; n )

4

H e ; (8.1)

i.e. tritium nuclei bombarded with deuterons to produce helium nuclei and

the neutrons being measured. If the bombarding particles are monoenergetic

then so are the neutrons produced, with energy in the range 12 - 22 MeV

(million electron volts), depending on that of the bombarding particles. Although the neutrons are monoenergetic, the measuring instrument both smears

and distorts the expected spectrum. The instrument in question, an NE-213

spectrometer, has been described by Verbinski, et. al. [20] and Burrus and

Verbinski [5]. Its response functions are plotted in Figure 7, where incident

energy corresponds to the variable  and pulse height to the variable t in Eqs.

(1.1). The incident neutrons are absorbed by a plastic scintillator which emits

a light pulse whose size is proportional to the incident energy. The light pulse

is in turn detected and amplied by a photomultiplier tube which transmits

a corresponding electrical pulse to an electronic circuit which then increments

a counter in one of m = 133 pulse-height bins designed to cover and subdivide the energy range of all possible incident neutrons. The integral equations

modelling the process are

ci =

Z

Eup

0

Ki(E)N(E)dE + "i ; i = 1; 2; : : : ; 113 ; (8.2)

where ci is the number of pulses counted in the ith bin, N(E) is the unknown

number of neutrons at energy E, and the Ki(E) are the instrument response

functions for the whole detector system. For a given value of E, the quantity

Ki(E)dE is the probability that a neutron with energy in the range E  1

2 dE

would produce a light pulse that would increment the count in the ith pulseheight bin. Hopefully, the bin number in which a neutron is counted would

increase as the energy increases, but an inspection of Figure 7 shows that this

does not always happen. Ideally the Ki(E) should be a series of narrow peaks

distributed along a diagonal line running from the upper left to the lower right

of the energy, pulse-height domain. The gure does exhibit a ridge along that

17
Figure 7: Instrument response functions for the Burrus neutron spectrum

problem. The quantity log10[1 + Ki(E)] is plotted in order to more clearly

show the structure for higher energies.

direction, but it attenuates to a barely discernible ripple for higher energies.

Even worse, the long energy tails for the lower pulse-height bins make it more

likely that a higher energy neutron will produce an increment in a lower rather

than a higher bin.

Figure 8 shows a measured spectrum obtained for the neutrons from the

reaction (8.1). Note that count rate rather than counts are plotted against

pulse-height. The experimental procedure is to accumulate counts for a time

sucient to allow good estimates to be made for the uncertainties in the count

totals and then to divide by the elapsed time to get rates. In Eq. (8.2) this

division replaces the ci with count rates yi and the number N(E) of neutron

18
Figure 8: Measured pulse-height spectrum for the Burrus problem. The two

piecewise linear curves are the  one standard deviation bounds for the measured values

with a neutron 
ux x(E). The result is

yi =

Z Eup

0

Ki(E)x(E)dE + i ; i = 1; 2; : : : ; 113 ; (8.3)

where the i are noise rates obtained from the "i. Estimates for the standard

deviations of the latter quantities were obtained from the usual assumption

about counting errors which is that the standard deviation for any bin is

approximately equal to the square root of the number of counts in that bin.

Although this procedure assumes a Poisson distribution for the error in any

bin, the number of counts accumulated was sucient to permit the use of the

normal approximation. Assuming that the measuring errors in any two bins

are statistically independent and dividing the estimated standard deviations

by the time of accumulation gives a diagonal variance matrix S2for the errors.

When the energy E is equated to the variable , the equations (8.3) become

identical to those in (1.2). They were discretized using simple rectangular

quadrature with an n = 77 point energy mesh, 0:2 M eV = E1 < E2 <    <

19
Figure 9: Singular values (squares) and rst n elements of jUT

bj (circles)

for the Burrus neutron spectrum problem. The horizontal line represents the

truncation level  = 3:5.

E77 = 18:91 M eV , to give the linear regression model (1.3). Multiplying

through by S1gave the normalized model (2.3). Although the two curves

plotted in Figure 8 indicate good statistics for the measurements, there is no

sign of the dominant peak expected for monoenergetic neutrons. The sharp

rise in the lower pulse-height bins contains most of the counts that should have

gone into that missing peak. The least squares estimate, which is overwhelmed

by the amplied noise, gives

r

2

min = kb  Ax^k

2

= 38:05 (8.4)

which is almost 5 standard deviations smaller than the expected value 113.

The singular values and rst n elements of jUTbj are shown in Figure 9.

The condition number of the matrix is 2:3503 105so it is clearly not rankdecient. The distribution of the jUT

bj

i exhibits a dichotomy at i = 47 where

it rst drops below the one sigma level for the error. Table 2 gives estimate

diagnostics for several values of  . The quantities tabulated are the same as

those in the rst 4 columns of Table 1 All of the non-zero  give estimates

20
Table 2: Estimate diagnostics for the Burrus neutron spectrum problem

 P(b  Ax)

2

i %(Ck 2 B95) LengthfC(f )g

0.0 38.0 37.4 1.255

3.0 93.8 100.0 1.211

3.5 125.7 100.0 1.212

4.0 150.7 100.0 1.211

5.0 213.9 63.6 1.216

Figure 10: Truncated solution and diagnostic plots for the Burrus neutron

spectrum problem with  = 3:5

21
with a strong peak at approximately 14 MeV, but only  = 3:0; 3:5; and 4:0

give residuals consistent with white noise. The sum of squared residuals for

 = 3:5 is closest to the expected value, so it was chosen as optimal even though

the length of the cumulative periodogram was slightly larger than those for

 = 3:0 and 4:0. The estimate and diagnostic plots for  = 3:5 are given in

Figure 10. The estimate is dominated by a single peak centered at 14 MeV

which is just what was expected for monoenergetic neutrons.

9 Discussion and Conclusions

The primary message of this paper is that for ill-posed problems in which the

dominant errors are the measurement uncertainties for the right hand side

vector, the instability in estimating the solution is attributable to components

of the measurements which are overwhelmed by the random error. In general,

these components will not correspond exactly to specic elements of the right

hand side vector. Scaling the problem by S1reduces it to one in which the

errors all have unit variance, and premultiplying the scaled measurement vector by UT

rotates it into a basis in which the noise dominated components

can be isolated by well understood statistical criteria. Zeroing the noise dominated components stabilizes the solution estimate without altering the kernel

matrix which is assumed to be known more accurately than the right hand

side vector. Of course this procedure does not guarantee that the estimated

solution will be close to the true solution since the noise may be large enough

to overwhelm important components of the signal. The only way to retrieve

such components would be to repeat the measurments with more accuracy.

The traditional method of truncating the singular value decomposition by

zeroing some of the singular values amounts to altering a well known quantity,

A, in order to compensate for uncertainties in an unknown quantity b = Ax.

In the process it reassigns the rank of the matrix, and since the truncation

level depends crucially on the size of the errors in the right hand side vector,

it follows that the altered rank depends on the size of those errors. As an

example, consider Figure 11 which is a plot of the i and jUTbj

i for a problem

obtained from the one illustrated in Figure 3 by increasing all the errors by

a factor of 10. Both gures display the 
at segment of the jUTbj curve,

which corresponds to error saturation, but for the problem with the larger

errors, this 
at segment intersects the singular value distribution at a point

higher and to the left of that for the problem with the smaller errors. Thus, if

the truncation required to stabilize the estimate were accomplished by zeroing

singular values, the problem with larger errors would assign a smaller numerical

22
Figure 11: Singular values (squares) and rst n elements of jUT

bj (circles)

obtained by increasing the errors in the modied Phillips problem by a factor

of 10

rank to the matrix than the one with smaller errors. Although it seems absurd

to have the rank of a matrix depend on the right hand side vector, precisely

this suggestion was made quite recently by P.C. Hansen [12, Chapt. 7] who

dened an \eective numerical rank" which is not even an integer. A better

approach might be to simply abandon the notion of \numerical rank", but the

huge vested interest in rank determination and rank revealing decompositions

will probably preclude such a course of action.

The truncation procedure described in this paper depends critically on

the availability of good estimates for the uncertainties in the elements of the

measured vector. A good experimenter will usually provide such estimates,

but it is not uncommon for them to be systematically too small or, more

rarely, too large. Often the form of the functional dependence of the errors

on the magnitudes of the measurements is correct but the estimated standard

deviations are all too small, or too large, by a constant factor. In such cases,

the diagnostics introduced in Section 4 can be used to pinpoint the discrepancy

and rectify it. The usual symptom is that the threshold  at which the estimate

23
begins to be stabilized is quite dierent from 3 or 4 and also dierent from the

value which gives residuals closest to white noise, as measured by the length of

the cumulative periodogram. In such cases it is a straightforward though often

tedious matter to adjust the scale of the estimated errors by a constant factor

which brings the two values of  into agreement somewhere in the interval

3    5.

The situation is more dicult when no estimates are given for the measurement errors. The usual approach is to implicitly assume that S2 = s

2

Im

where s is an unknown constant which can be estimated from the least squares

sum of squared residuals. Of course the latter quantity is usually much smaller

than the expected value m, but the size of s can then be adjusted upward to a

value such that some truncation level in the range 3:0    5:0 gives residuals

close to white noise and a sum of squared residuals close to the expected value.

This usually iterative technique works well if the standard deviations of the

measured errors are truly constant. If not, it may still be possible to assume a

functional form for the uncertainties which will determine the matrix S2up to

an unknown scaling constant which can then be adjusted for consistency in the

 values. The disadvantage in this approach is the amount of work involved

in validating a functional form and adjusting the scaling, but sometimes the

results obtained justify the eort.

A number of people have suggested that the truncation procedure advocated in this paper is incomplete because it has not been accompanied by an

algorithm for automatically determining the optimum  , using some criterion

like those used in choosing the smoothing parameter for regularization methods, e.g., generalized cross-validation, or the L-curve method. In practice, for

any given problem, only a few discrete truncations give estimates with residuals that have both a reasonable sum of squares and an acceptable spectral

distribution, and it is often the case that the two criteria are optimized by

dierent truncation levels. Thus some judgement is required to choose which

truncation level to use. Also there is little to be gained computationally by

conning the calculations to a single  . The amount of work required is dominated by the computation of the singular value decomposition, but once that

is accomplished it is relatively easy to isolate the few possibly acceptable truncation levels and, using a fast Fourier transform algorithm, to calculate all of

the diagnostics. Since the total eort involved is small in comparison to that

required to make the measurements, most users with real world problems will

not be unduly troubled by having to make the nal choice of the optimum

value. Of course the procedure can be automated by calculating the solutions

and diagnostics on a suitably chosen mesh for  . If the mesh spacing is ne

enough to nd all of the possibly acceptable solutions, it will produce many

24
duplicate solutions, but these will not seriously complicate the nal choice of

the parameter value. Users of regularization methods might also be well advised to consider the same two residual diagnostics in evaluating their favorite

criteria for choosing the regularization parameter.

Acknowledgements

I would like to thank Drs. W.R. Burrus and R.E. Funderlic for teaching me the

pitfalls of truncating the singular value distribution, Drs. D.P. O'Leary and

K.A. Remington for their reviews and suggestions for improving the numerical

analysis, and Dr. M.G. Vangel for many useful suggestions for improving the

statistical aspects of this manuscript.

Appendices

A Percentiles for the 2(m) Distribution

The cumulative distribution funtion for 2(m) is dened by

F (2) =

Z

2

0

1

2

m

2 



m

2

 z

m2

2 exp 



z

2



dz ; (A.1)

and for any probability , the -point of the distribution can be found by

solving F (2) = . If the solution is 2

, then

Pr nkb  Axk

2

 2



o

=  ; (A.2)

and

Pr nkb  Axk

2

> 2



o

= 1   : (A.3)

Standard tables (e.g., [3]) of the 2-distribution list 2

 for various values of

 and m, up to about m = 30. For m > 30, the quantity

p

22 

p

2m  1



is distributed approximately like the standard normal distribution, so 2

 is

well approximated by

2

 =

1

2

h

 +

p

2m  1

i

2

; (A.4)

where  is the -point of the standard normal distribution, i.e.,

Z



1

1

p

2

exp 

z

2

2

!

dz =  : (A.5)

25
Table 3: Selected percentiles for the standard normal and 2

distributions

  2



m = 100 m = 150 m = 200 m = 300

0.01 -2.326 69.4 112.0 155.7 245.3

0.05 -1.645 77.7 122.4 168.0 260.6

0.33 -0.440 93.4 142.0 190.8 288.8

0.50 0.0 99.5 149.5 199.5 299.5

0.67 0.440 105.8 157.2 208.4 310.4

0.95 1.645 124.1 179.3 233.7 341.1

0.99 2.326 135.0 192.4 248.7 359.1

Table 3 gives, for various values of , the corresponding values of  and, for

some representative values of m, the estimates of 2

 calculated by Eq. (A.4).

Good public domain software for computing 2

 can be found in the GAMS

collection [7, Class L5a1c].

B A Variant of the Phillips Problem

A useful test problem which shares many of the characteristics of real instrument correction problems is obtained by discretizing a variant of the well

known [16] Phillips equation. This modied Phillips problem can be written

y(t) =

Z

3

3

K(t; )x() d ; 6  t  6 ; (B.1)

with the kernel given by

K(t; ) =

8

<

:

1

6

n

1 + cos h



3

(  t)

io ; j  tj  3 ; jtj  6 ;

0 ; otherwise ;

(B.2)

and the exact solution by

x() = () + X

3

k=1

ck() ; (B.3)

where

() =

8

<

:

A0

h

1 + cos 



3



i

; jj  3 ;

0 ; jj > 3 ;

(B.4)

26
and

ck() =

8

<

:

Ak f1 + cos [2 (  k)]g ; j  kj  1

2

;

0 ; otherwise ;

(B.5)

with amplitude constants Ak and centering constants k chosen to be

A0 = 0:1 ; A1 = 0:5 ; A2 = 0:5 ; A3 = 1:0 ;

 1 = 1:5 ; 2 = 0:5 ; 3 = 1:5 :

(B.6)

The kernel diers from the Phillips original only in the inclusion of the normalizing factor 1

6 which was added to assure that for any ,

Z

3+

3+

K(t; ) dt = 1 : (B.7)

For a measuring instrument this condition assures that conservation laws are

not violated. Plots of K(t; ) for 5 representative values of  are given on the

left in Figure 1.

The exact solution to the original Phillips problem appears in a scaled down

form as the () term in the solution to the new problem. The scaling constant

A0 was chosen to reduce the original solution to the role of a background

function on which to superimpose the three discrete peaks represented by the

ck() terms. The three points  = k are the centers of these peaks and

the constants 2Ak are their heights above the background. The new solution

function is plotted as a dashed line on the right in Figure 1.

These changes in the Phillips problem were designed to make it more challenging and more reminiscent of real-world instrument correction problems.

Unfortunately they also make the representation of the function y(t) more

complicated. Substituting (B.3) into (B.1) gives

y(t) =

Z

3

3

K(t; )()d + X

3

k=1

Z

3

3

K(t; )ck()d ; 6  t  6 ; (B.8)

but care must be exercised in evaluating these integrals because K(t; ) = 0

on half of the rectangular domain f (t; ) j  6  t  6 ; 3    3 g and

each of the ck() is zero everywhere except on the interval k  1

2    k + 12

.

The last equation can also be written

y(t) = B(t) + X

3

k=1

Ck(t) ; (B.9)

where

B(t) 

Z

3

3

K(t; )()d ; (B.10)

27
and

Ck(t) 

Z

 k+ 1

2

 k 1

2

K(t; )ck()d ; k = 1; 2; 3 : (B.11)

Evaluating the integral for B(t) gives

B(t) =

1

6

A0



(6  jtj)



1 +

1

2

cos 



3

t

 +

9

2

sin 



3

jtj



; 6  t  6 ;

(B.12)

and the integrals for Ck(t) can be written

Ck(t) =

1

6

AkLk(t) ; k = 1; 2; 3 ; (B.13)

where

L1(t) =

8

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>:

0 ; 6  t  5;

t + 5 + 1

2

sin [(2t + 9)] + 3



sin h



3

(t + 2)i

+ 3

10

n

sin [(2t + 8)]  sin h



3

(t  1)io

+ 3

14

n

sin [(2t + 10)] + sin h



3

(t + 5)io ; 5  t  4;

1 + 3



n

 sin h



3

(1 + t)

i

+ sin h



3

(2 + t)

io

+ 3

10

n

sin h



3

(4 + t)

i

+ sin h



3

(1  t)

io

+ 3

14

n

sin h



3

(2  t)

i

+ sin h



3

(5 + t)

io ; 4  t  1;

2  t + 1

2

sin [(3  2t)]  3



sin h



3

(t + 1)i

+ 3

10

n

sin h



3

(t + 4)i+ sin [(2  2t)]o

+ 3

14

n

sin h



3

(2  t)

i

+ sin [(4  2t)]o; 1  t  2;

0 ; 2  t  6;

(B.14)

28
L2(t) =

8

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>:

0 ; 6  t  3;

t + 3 + 1

2

sin [(2t + 5)] + 3



sin h



3

t

i

+ 3

10

n

sin [(2t + 4)]  sin h



3

(t  3)io

+ 3

14

n

sin [(2t + 6)] + sin h



3

(t + 3)io ; 3  t  2;

1 + 3



n

sin h



3

(1  t)

i

+ sin h



3

t

io

+ 3

10

n

sin h



3

(2 + t)

i

+ sin h



3

(3  t)

io

+ 3

14

n

sin h



3

(4  t)

i

+ sin h



3

(3 + t)

io ; 2  t  3;

4  t + 1

2

sin [(7  2t)] + 3



sin h



3

(1  t)

i

+ 3

10

n

sin h



3

(t + 2)i sin [(2t  6)]o

+ 3

14

n

sin h



3

(4  t)

i

 sin [(2t  8)]o; 3  t  4;

0 ; 4  t  6;

(B.15)

and

L3(t) =

8

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>:

0 ; 6  t  2;

t + 2 + 1

2

sin [(2t + 3)] + 3



sin h



3

(t  1)i

+ 3

10

n

sin [(2t + 2)]  sin h



3

(t  4)io

+ 3

14

n

sin [(2t + 4)] + sin h



3

(t + 2)io ; 2  t  1;

1 + 3



n

sin h



3

(2  t)

i

 sin h



3

(1  t)

io

+ 3

10

n

sin h



3

(1 + t)

i

 sin h



3

(t  4)io

+ 3

14

n

sin h



3

(5  t)

i

+ sin h



3

(2 + t)

io ; 1  t  4;

5  t  1

2

sin [(2t  9)] + 3



sin h



3

(2  t)

i

+ 3

10

n

sin h



3

(t + 1)i sin [(2t  8)]o

+ 3

14

n

 sin h



3

(t  5)i sin [(2t  10)]o; 4  t  5;

0 ; 5  t  6;

(B.16)

The function y(t) is plotted as the solid curve on the right in Figure 1. All of

29
the details of the 3 peaks are so smeared together by the convolution with the

kernel function that that there is no hint of any structure in the underlying

x().

C Discretizing the Modied Phillips Problem

The modied Phillips problem described in Appendix B was discretized by

choosing m = 300 equally spaced mesh points on the interval 5:9625  t 

5:9625 to give

y(ti) =

Z

3

3

K(ti ; )x() d ; i = 1; 2; : : : ; 300 ; (C.1)

and by replacing each of the integrals by an n = 241 point trapezoidal quadrature sum, i.e.,

Z

3

3

K(ti; )x()d  X

241

j=1

!jK(ti ; j )x(j ); i = 1; : : : ; 300 ; (C.2)

where

( !1 ; !2 ; !3 ; : : : ; !240 ; !241 )

T

=

0:025

2

( 1 ; 2 ; 2 ; : : : ; 2 ; 1 )

T

: (C.3)

Dening the n-vector xand the m-vector yby

x

j = x(j ) ; j = 1 ; 2 ; : : : ; 241 ;

y

i = y(ti) ; i = 1 ; 2 ; : : : ; 300 ;

(C.4)

and the m n matrix K by

Ki;j = !jK(ti; j ) ; i = 1 ; 2 ; : : : ; 300 ;

; j = 1 ; 2 ; : : : ; 241 ;

(C.5)

gives

y= Kx+  ; (C.6)

where  is an m-vector of quadrature errors. A crucial assumption in replacing

the intregrals with quadrature sums is that the value of n is chosen large enough

so that the i are small relative to the random measuring errors i. To assure

that this assumption was satised for the test problem, the elements of the

vector ywere computed from the matrix-vector product

y= Kx(C.7)

30
rather than from Eqs. (B.9) - (B.16). More precisely, the matrix elements

Ki;j were computed from Eqs. (B.2), (C.3), (C.5), the vector elements x

j were

computed from Eqs. (B.3) - (B.6), and the vector ywas then computed from

Eq. (C.7). The \measured" vector y was then obtained by adding random

perturbations to the elements of this y. Each perturbation was chosen independently from a normal distribution with mean zero and standard deviation

si = (105)

p

y

i

, so the variance matrix was

S2= diag(s

2

1; s

2

2; : : : ; s

2

300) ; si = (105

)

q

y

i

: (C.8)

D The Cumulative Periodogram

For a given time series, a periodogram is an estimate of the power spectrum of

the series, i.e., an estimate of how the total variance in the series is distributed

in frequency. Such an estimate can be obtained at any desired number of frequencies fk; k = 1; 2;    ; N, in the interval 0  f  1

2T

, where T is the sample

spacing for the time variable, and N is chosen to be greater than or equal to

the number of sample points in the time series. If the time series is taken to

be the residuals ri for an estimated solution to an ill-posed problem, with the

element number i as the time variable, then T = 1, and the periodogram is

estimated on an N-point equally spaced mesh on the interval 0  f  0:5. It

is obtained by zero-padding the ri series to have N  m terms and computing

the discrete Fourier transform

Rk = T X

N

j=1

rj exp  i 2

kj

N

!

; k = 0; 1; 2; : : : ;

N

2

; (D.1)

where i =

p

1 and

rj =

(

bj  (Ax^)j ; j = 1; 2; : : : ; m ;

0 ; j = m + 1; m + 2; : : : ; N :

(D.2)

The zero-padding increases the density of the frequency mesh, but does not

change the value of the transform at any given frequency. Each Rk is associtated with the corresponding Fourier frequency fk = k

N T . The periodogram is

computed from the transform by

Pk = P(fk) =

1

N T jRkj

2

; k = 0; 1; 2; : : : ;

N

2

: (D.3)

Peaks in the plot of Pk versus fk indicate the presence of sinusoidal cycles in

the parent time series, and the height of each such peak gives an estimate of

the power of that cycle.

31
The cumulative peridogram is dened by

C0 = C(0) = 0 ; Ck = C(fk) =

1

S

X

k

j=1

Pj ; k = 1; 2; : : : ;

N

2

; (D.4)

where

S =

N=

X

2

j=1

Pj : (D.5)

Clearly C(f ) is a monotonic nondecreasing function of frequency ranging between the values C(0) = 0 and C(1=2T ) = 1. The relationship between the

periodogram and the cumulative periodogram is illustrated in Figures 12 and

13. In the rst case the time series is comprised of m = 128 independent

random samples xi from a standard normal distribution (mean = 0, variance

= 1), and in the second case the series is given by

xi = sin 

2i

10 

+ ni ; i = 1; 2; : : : ; 128 ; (D.6)

where the ni are independent samples from a zero-mean normal distribution

with variance = 0.1. Thus the rst series is a realization of a pure white noise

signal and the second is a single sinusoid, with unit amplitude and period

10, corrupted by zero-mean white noise with standard deviation 0.316. In

both cases the time series was rst detrended by subtracting out the mean

value which was approximately, but not exactly, zero. The periodogram for

the white noise case has the power distributed more or less uniformly on the

interval [0; 0:5], so the cumulative periodgram does not depart too much from

a diagonal line that would connect the points (0; 0) and (0:5; 1), i.e., the line

C(f ) = 2T f ; (D.7)

which is the theoretical distribution for pure white noise. One way to quantize

departures from white noise is to compare the length of that line,

q

(0:5)2 + (1:0)2 = 1:11803 ;

with that of the computed estimate,

Length(C) =

N=

X

2

k=1

q

(Ck  Ck1)

2

+ (fk  fk1)

2



; (D.8)

which is given in the legend box on the plot of the cumulative periodogram.

32
Figure 12: Periodogram and cumulative periodogram for a white noise time

series

Figure 13: Periodogram and cumulative periodogram for a single sinusoid plus

white noise time series

33
The periodogram for the sinusoid plus noise time series is dominated by

a single peak at frequency 0.1 (period = 10) which produces a wide deviation from the white noise diagonal at the same frequency in the cumulative

periodogram. One way to test whether of not such a deviation represents

a statistically signicant departure from white noise is to construct the two

parallel o-diagonal lines dened by

C(f ) =  + 2T f ; (D.9)

where  is the 5% point of the Kolmogorov-Smirnov statistic for a sample of

size m=2. These two lines, which are plotted in both Figures 12 and 13, enclose

a 95% condence band which can be used to test the hypothesis that the time

series is a realization of white noise. The cumulative periodogram ordinates

for such a series should lie outside this band for at most 5% of the frequencies.

None of the ordinates lie outside the band for the white noise case, but more

than 50% of them do for the sinusoid plus noise case.

References

[1] Anderson, E., Bai, Z., Bischof, C., Demmel, J., Dongarra, J., Du Croz,

J., Greenbaum, A., Hammarling, S., McKenney, A., Ostrouchov, S. and

Sorensen, D. (1992), LAPACK Users' Guide, SIAM, Philadelphia.

[2] Anderson, T.W. (1958), An Introduction to Multivariate Statistical Analysis, John Wiley & Sons, New York, Chapt. 2.

[3] Beyer, W.H. (1968), Handbook of Tables for Probability and Statistics,

CRC Press, Boca Raton, pp. 293-294.

[4] Boisvert, R.F., Howe, S.E. and Kahaner, D.K. (1984), The guide to available mathematical software (GAMS), PB 84-171305, National Technical

Information Service, Springeld, VA.

[5] Burrus, W. and Verbinski, V. (1969), Fast neutron spectroscopy with thick

organic scintil lators, Nuclear Instruments and Methods, 67, pp.181-196.

[6] Fuller, W.A. (1976), Introduction to Statistical Time Series, John Wiley

& Sons, New York, Chapt. 7.

[7] GAMS (1998), Guide to Available Mathematical Software,

http://gams.nist.gov/.

34
[8] Golub, G. and Kahan, W. (1965), Calculating the singular values and

pseudo-inverse of a matrix, J. SIAM Numer. Anal. Ser. B, 2, pp. 205-224.

[9] Golub, G., Klema, V. and Stewart, G.W. (1976), Rank Degeneracy and

Least Squares Problems, Technical Report TR-456, Department of Computer Science, University of Maryland.

[10] Graybill, F.A. (1976), Theory and Application of the Linear Model,

Duxbury Press, North Scituate, Chapt. 6.

[11] Hansen, P.C. (1987), The truncated SVD as a method for regularization,

BIT, 27, pp. 534-553.

[12] Hansen, P.C. (1998), Rank-Decient and Discrete Il l-Posed Problems,

SIAM, Philadelphia, Chapt. 4.

[13] Hanson, R.J. (1971), A numerical method for solving Fredholm integral

equations of the rst kind using singular values, SIAM J. Numer. Anal.,

8, pp. 616-622.

[14] Hogg, R.V. and Craig, A.T. (1965), Introduction to Mathematical Statistics, Macmillan, New York, Sect. 4.6.

[15] Kahaner, D., Moler, C., and Nash, S. (1989), Numerical Methods and

Software, Prentice Hall, Englewood Clis, Chapt. 3.

[16] Phillips, D.L. (1962), A technique for the numerical solution of certain

integral equations of the rst kind, J. Assoc. Comput. Mach., 9, pp. 84-

97.

[17] Smithies, F. (1958), Integral Equations, Cambridge University Press,

Cambridge, Chapt. 8.

[18] Tarantola, A. (1987), Inverse Problem Theory, Elsevier, Amsterdam,

Chapt. 1.

[19] Twomey, S. (1977), Introduction to the Mathematics of Inversion in Remote Sensing and Indirect Measurements, Elsevier, Amsterdam, Chapt.

6.

[20] Verbinski, V., Burrus, W., Love, T., Zobel, W., Hill, N., and Textor, R.

(1968), Calibration of an organic scintil lator for neutron spectrometry,

Nuclear Instruments and Methods, 65, pp. 8-25.

35
[21] Wing, G.M. (1991), A Primer on Integral Equations of the First Kind,

SIAM, Philadelphia, Chapt. 2.

36