Scuba Half-Degree Extragalactic Survey: Data Reduction by Iterative Deconvolution by Kyle Quentin Lepage BASc, The University of British Columbia, 2000 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Department of Physics and Astronomy) We accept this thesis a$ conforming to the required /tjandard THE UNIVERSITY OF BRITISH COLUMBIA April 21, 2004 © Kyle Quentin Lepage, 2004 Library Authorization In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Name of Author (please print) Title of Thesis: ^>ceo Degree: KS-O \ ,.,L„ C Date (dd/mm/yyyy) yUlf - D ^ r ^ E * 4 r a ^.!« r-l,V r^> H^-Ur <v£ S p ' i ^ Department of \>u *<_s. & The University of British Columbia Vancouver, BC Canada y ^ V ^ c ^ - v Y e a r : ^&c>h Abstract n Abstract A generalization of an iterative algorithm used to reduce large data sets consisting of observations of the cosmic microwave background is presented. In particular, the algorithm is modified to reduce data sets obtained by arbitrary observations. The generalized algorithm is applied, in simulation, to a data set consisting of complicated observations taken during a Scuba HAlf Degree Extragalactic Survey (SHADES). The algorithm is found to be unstable. The generalized time-order algorithm is reformulated as a difference equation whose solution gives the N th iteration reduced data map in analytical form. Within the analytical framework, the convergence properties of the algorithm are determined, a stable, known work-around explained, and the converged map shown to be optimal in the least-squares sense. The known work-around, a modification of the generalized iterative algorithm, trades stability for ideal behaviour at the edge of the reduced data map. Using simulations of a sky populated by a realistic source count model, it is shown that the deleterious effect of the modified algorithm towards the edge of the map is negligible as compared to the statistical noise for most of the half-degree square map. Simulations of the modified algorithm applied to data taken with the SHADES measurement strategy with Gaussian noise added, reveal the presence of streaking and large scalefluctuationsin the converged map. The streaking which is demonstrated to be affected only by the measurement strategy, is indicative of a banddiagonal pixel-pixel covariance matrix, and is absent in noiseless simulations. The large scalefluctuationsare indicative of non-negligible, far off-diagonal terms in the pixel-pixel covariance matrix, and are also absent in noiseless simulations. The effect of different observation strategies on the streaking is investigated by simulation. The presence of the amplified large scale noise is explained by linear systems theory. Abstract iii The determination of the optimal observation strategy is formulated as the task of determining the linear combination of absolute power measurements which will yield a desired pixel-pixel covariance matrix. This problem is discussed. An un-deconvolved map of real SHADES data is produced. Sources are detected in the un-deconvolved map using Bayesian techniques. Contents iv Contents Abstract ii Contents iv List of Tables vii List of Figures ix Acknowledgements xviii I Thesis 1 1 The Scuba Half-Degree Extragalactic Survey 2 1.1 Introduction To The Submillimeter Galaxy Population 3 1.2 What Is The Cosmic History Of Massive Dust-Enshrouded Star Formation Activity? 4 1.3 Are SCUBA Sources The Progenitors Of Present-Day Massively Elliptical Galaxies? 5 1.4 What Fraction Of SCUBA Sources Harbour A Dust-Obscured Active Galactic Nucleus? 1.5 Legacy Value 5 6 1.5.1 Data Compression 6 1.5.2 Data Correlation: Integrated Sachs-Wolfe Effect 6 2 Observing The Scuba HAlf-Degree Extragalactic Survey 8 2.1 Submillimeter Common User Bolometer Array . 8 2.2 Data Acquisition 9 Contents v 2.3 SCUBA/JCMT Measurement Noise 11 2.3.1 Brief Overview Of Noise 11 2.3.2 13 SCUBA Noise 2.4 Data Preprocessing 2.4.1 Calibration 2.4.2 Nuisance Signals And Outlier Rejection 2.5 Data Reduction 14 14 18 21 2.5.1 Interpolative Reduction 23 2.5.2 27 Emerson Reconstruction 2.5.3 Direct Inversion 28 2.5.4 Iterative Inversion 29 3 Iterative Reduction of S H A D E S D a t a 3.1 Wright's Algorithm Generalized 3.1.1 Algorithm Checks 3.2 Generalized Algorithm Simulations 31 31 32 33 3.2.1 SHADES Simulation 34 3.2.2 Small Field Simulations 37 3.2.3 Minimum Comparable Size 43 3.3 Iterative Inversion As Newton-Raphson Root Finding 45 3.4 Formulation of Iterative Inversion as a Difference Equation 46 3.4.1 Formulation 46 3.4.2 Solution 48 3.4.3 Solution Check #1 49 3.4.4 Solution Check #2 50 3.4.5 Solution Check #3 50 3.4.6 Difference Equation Solution Simulation 52 3.4.7 Stability 54 3.4.8 Noise Properties 57 3.5 Algorithm Stabilization 59 3.5.1 Unstable Eigenvector Removal 59 Contents 3.5.2 Boundary Pixel Modification 61 3.5.3 Measurement And Pixel Update Decoupling 62 4 Central Update A l g o r i t h m 5 6 vi 71 4.1 Edge Effects 71 4.2 Gaussian Noise Simulation 81 4.3 More Realistic Noise Simulation 87 4.4 Chopping Effects On Pixel-Pixel Correlation 92 4.5 Chapter Results 92 S H A D E S Subaru X M M - D e e p Field 95 5.1 Data Reduction 95 5.2 Source Detection 98 5.3 Deconvolved Map Source Detection 106 Discussion 113 6.1 Main Results 113 6.2 Follow-Up Observation 114 6.3 Measurement Strategy 116 6.4 Optimal Measurement Strategy 118 6.5 Noise &; Correlation Removal 119 A Jordan Form 121 Bibliography 123 List of Tables vii List of Tables 3.1 Correlation M a t r i x : For a seven pixel, seven measurement, singledifference simulation 40 3.2 Correlation M a t r i x : Seven pixel, thirteen measurement, doubledifference simulation 42 3.3 Correlation M a t r i x : Seven pixel, fifteen measurement, doubledifference simulation 42 3.4 Eigenvalues of the characteristic matrix for the seven pixel, five double-difference measurement simulation 52 3.5 Eigenvectors of the characteristic matrix in the seven pixel, five double-difference measurement simulation. Each row represents a separate eigenvector. The first row is the eigenvector corresponding to the left-most eigenvalue in Table 3.4. For clarity, the eigenvector components have been rounded to three significantfigures.The eigenvectors are not linearly independent 52 3.6 Eigenvalues of the characteristic matrix for the seven pixel, 15 double-difference measurement simulation 53 3.7 Eigenvectors of the characteristic matrix for the seven pixel, 15 double-difference measurement simulation. Each row is a separate eigenvector; the eigenvector in the top row corresponds to the left most eigenvalue in Table 3.6 53 3.8 Eigenvalues of the characteristic matrix for the six pixel, four doubledifference measurement simulation with two fixed pixels. There is an absolute eigenvalue greater than one 56 List of Tables viii 3.9 Eigenvalues of the characteristic matrix in the simulation with six pixels, two double-difference measurements and two absolute measurements. There is an absolute eigenvalue greater than one 56 3.10 Eigenvalues of the characteristic matrixforthe six pixel,fourdoubledifference measurements and two absolute measurements simulation. There is an absolute eigenvalue greater than one 57 3.11 Eigenvalues of the characteristic matrixforthe six pixel,fourdoubledifference measurement and two absolute measurement. There are no absolute eigenvalues greater than one 58 3.12 Eigenvalues of the characteristic matrixforthe seven pixel, five double-difference measurement simulation 61 3.13 Characteristic eigenvalues for the edge-pixel enlarged pointing matrix 62 3.14 Characteristic matrix eigenvalues for the central update algorithm in the seven pixel, five double-difference measurement simulation. All of the absolute eigenvalues are less than or equal to one 63 5.1 Subaru X M M Deep Field Source Detection List - Part 1 110 5.2 Subaru X M M Deep Field Source Detection List - Part II .... I l l 5.3 Subaru X M M Deep Field Source Detection List - Part III ... . 112 List of Figures ix List of Figures 2.1 SCUBA 850 micrometer bolometer array. The plot is obtained from the positions of actual data on the sky. Note the spacing between equal distances along the two axes are different, the height and width of the array are equal 2.2 Tripos Coverage 2.3 Eight representative scuba bolometer timestreams and spectra for the 9 11 Subaru XMM Deep Field afterflatfielding,and extinction correction. The timestreams are non-stationary and strongly correlated by atmospheric emission. The timestream for bolometer 21 extends past the range of the y-axis. The standard deviation of the bolometer timestreams varies substantially from bolometer to bolometer. The magnitude of the Fourier transform of the bolometer timestreams are shown in the lower half of the plot. Note the flat spectra and the spikes on bolometers 17, 19, 20, 21, and possibly 23 and 24 16 2.4 Flux conversion factor as a function of observation file number. The observation files are not in time order. Note the considerable fluctuation of the FCF with observation number. Only calibrations with CRL618, Uranus and Mars are used. The chop throw of the calibrations is ignored 17 List of Figures x 2.5 Correspondence between sources detected at greater than 3a significance from calibrated data and from uncalibrated data. The calibrated data uses Method 3 in Section 2.4.1 to compute flux conversion factors. The source detection is described in Section 5.2, but differs in that the detected sources have not been iteratively cleaned from the map 19 2.6 Eight bolometer timestreams and their corresponding Fourier transform magnitudes for the Subaru XMM Deep Field after preprocessing. Note the lack of outliers in the preprocessed timestreams as well as the reduced rms when compared to the timestreams in Figure 2.3. The Fourier transform magnitudes are flatter than their non-preprocessed counter-parts. (Figure 2.3) 22 2.7 A measurement map of simulated measurements taken with the single pointing SHADES observation strategy of a fake sky with a single, centred source. Note the presence of negative sidelobes centred upon the six chop throws surrounding the source. Each square pixel in this map is three arcseconds square 25 3.1 Residual maps of a simulation of the generalized iterative algorithm. The measurements are performed on a statistically plausible fake sky using the full SHADES measurement strategy. The measurements are noiseless. High spatial frequency structure is removed from the residuals in the early iterations. In later iterations the unstable spatial frequencies begin to grow. The fake sky map saturates at pixel values of zero and ten. Pixels are three arcseconds square 35 3.2 Residual maps for iterations 1 through 6 and 25 in a simulation of the generalized iterative algorithm applied to measurements taken of a single source. The measurements comprise a single pointing in the SHADES measurement strategy. 36 List of Figures xi 3.3 Top: Fake sky map. Middle: Fake sky map estimates for iteration 1, 11 and 21. Bottom row: Residual maps for iteration 1, 11 and 21. The colour bar applies to the residual maps only 64 3.4 Residual mapsforiteration 1, 15 and 25. Note the reduction of highfrequency spatial structure between iteration 1 and 15 and the reappearance of high-frequency structure in iteration 25. The algorithm slowly diverges 65 3.5 Single-Difference: 7 Pixel Simulation, 6 Measurements 66 3.6 Single-Difference: 7 Pixel Simulation, 7 Measurements 67 3.7 Double-Difference: 7 Pixel Simulation, 5 Measurements 68 3.8 Double-Difference: 7 Pixel Simulation, 13 Measurements 69 3.9 Double-Difference: 7 Pixel Simulation, 15 Measurements 70 4.1 Standard deviation of the residual map pixels as a function of iteration for a central update algorithm data reduction simulation. The simulated data consist of measurements of a single, lOmJy, centred source, using the full SHADES measurement strategy. The computation neglects pixels within 150" from the map edge. The standard deviation of the residual map pixels is 0.0008 mJy at iteration 101. At iteration 1001 the standard deviation of the residual map pixels is less than 0.0001 mJy 72 4.2 Residual maps estimated on various iterations of a simulation of the central update algorithm using the full SHADES measurement strategy. The fake sky map is a single source, centred in the map, with a peak flux of ten. The pixel values are in mJy and the contrast has been stretched on the lower three panels. The first iteration map is a measurement map. Note the convergence and the absence of a visible effect caused by the edges of the map. The maximum absolute pixel value of the residual map on the 75 iteration is less than a tenth of th a percent of the source amplitude 73 List of Figures xii 4.3 Standard deviation of the residual map pixels as a function of iteration during a central update algorithm simulation. The simulated measurements are of an edge only, statistically plausible fake sky map and are noiseless. Pixels within 150" of the map edge are not used in the computation. At iteration 401, the residual map pixel standard deviation is within 0.02 mjy of the converged value of 0.07 mjy. The steps in the curve between iterations 401 and 801 are due to limited precision when reading and writing data to file. Note the maximum residual map pixel standard deviation of 0.185 mjy near iteration 60. At iteration one, the initial map estimate is the null map, and all pixels used in the computation of the standard deviation are equal to the fake sky map. Thus the standard deviation of the residual map pixels is zero, and rises with iteration as pixels inside the periphery are modified due to flux in the chop positions of the measurements centred upon pixels at the edge of the periphery. 75 4.4 An edge-only fake sky map and residual maps for a noiseless, full SHADES simulation of the central update algorithm. The fake sky map is populated with sources drawn from a realistic source count model. The edges are 150" wide. The pixel values are in mjy. Note the growth of the edge error. The maximum residual pixel standard deviation occurs around iteration 60. Note how the mean of the residual map pixels grows increasingly negative, reducing the standard deviation. At iteration 101, the absolute pixel values are less than two and a half percent of a 10 mjy source, 360" from the map edge. The white lines are spaced 200 pixels apart, each pixel is 3" wide. 76 4.5 The 801 * residual map for the statistically plausible, edge only, cens tral update algorithm simulation. The residual due to the edge effect is less than 0.15 mjy, 150" from the map edge 77 xiii List of Figures 4.6 Residual map standard deviation as a function of iteration during a simulation of data reduction using the central update algorithm. The data consist of measurements of a statistically plausible fake sky obtained with the full SHADES measurement strategy. Convergence to within 0.005 mJy of the converged standard deviation occurs by iteration 50. The converged standard deviation of 0.068 mJy is 0.002 mJy from the the converged standard deviation of the residual pixels in the edge-only simulation (Figure 4.3). The minimum of 0.0655 mJy in the standard deviation of the residual pixels occurs near iteration 40 79 4.7 Residual maps for a full SHADES simulation of the central update algorithm. The fake sky map is populated with sources drawn from a realistic source count model. Between iterations 50 and 90 the variance of the residual pixels reduces in the centre of the map and slightly increases toward the map edges. The pixel values are in mJy. Notice the high spatial frequency residuals, visible at iteration one, are largely gone by iteration 40 80 4.8 Residual map average and standard deviation as a function of iteration during a simulation of data reduction using the central update algorithm. Gaussian random noise is added to simulated measurements of a fake sky containing no flux. At iteration 100, the residual pixel standard deviation is 1.95 mJy higher than the expected ^measurement/ViV The central-update algorithm converges when re- ducing measurements with Gaussian additive noise 82 4.9 Measurements with Gaussian noise are iteratively reduced using the central update algorithm. The maps are, clockwise from upper left, the 1 iteration map, the 101 iteration map, the power spectrum st st of the 1 iteration map and the power spectrum of the 101 and 1 st st st iteration map. Note the structure that appears in the 101 iteration st maps 83 List of Figures xiv 4.10 No source streaking, other than the existing cross-hatch pattern, is apparent, in a simulated reduction of measurements of four large sources containing Gaussian additive noise using the central update algorithm. The pixel values are in mjy 85 4.11 Measurements of a statistically plausible fake sky with additive Gaussian noise are iteratively reduced using the central update algorithm. The maps are, clockwise from upper left, the 1 iteration map, the st 101 st iteration map, the power spectrum of the 1 iteration map st and the power spectrum of the 101 iteration map. The presence st of sources appears as low frequency power in the 1 iteration power st spectrum map 86 4.12 First four plots are of the randomly extended noise sequences for bolometers 1, 2, 9, and 35. The bolometer selection is arbitrary. Last four plots are the corresponding Fourier transform magnitudes. Note the lack of j noise in the spectrums. There is a bump at approx. 0.007 Hz in the bolometer 35 spectrum, possibly the Fourier transform spike described in chapter 2 88 4.13 Residual map standard deviation as a function of iteration in a central update algorithm simulation with a statistically plausible fake sky map and the full SHADES measurement strategy 89 4.14 Clockwise from upper left: Sky estimate maps for iteration one, 101, power spectrum map for iteration 101, power spectrum map for iteration one. Note the similarity between these maps and the maps produced in simulations with Gaussian additive noise. (Figure 4.12) The off-diagonal terms of the timestream noise-covariance matrix play a minor role in the pixel-pixel correlations present in the 101 iterast tion map. The pixel-pixel correlations are predominantly due to the off-diagonal terms in the P P matrix for the SHADES measurement T strategy. The P P matrix was introduced in chapter two T 91 List of Figures xv 4.15 Simulation of data reduced using the central update algorithm. The data are simulated measurements of a statistically plausible fake sky with blank field extended noise samples. One half of the chops are rotated 45 degrees from the east-west and north-south chop directions. Much of the cross-hatching visibile in the 101 iteration map st in Figure 4.14 has disappeared. The strong cross-shape in the 101 st power spectrum map in Figure 4.14 has been reduced. Large scale structure with half-wavelengths around 600" seem unaffected 94 5.1 The Subaru XMM-Deep Field measurement map. Units are in Jy, pixels are 3" square. The horizontal and vertical axes are specified in units of pixels 96 5.2 Subaru XMM-Deep Field noise map and hit map, top and bottom, respectively. The noise map units are in Jy 5.3 97 Mask used in source detection tofitthe amplitude of a point centred on each pixel. Note the negative sidelobes. The central Gaussian has a full-width half-maximum of 14.7". Each pixel is 3" square. The mask is 150" square. The mask mean pixel intensity is zero 100 5.4 Subaru XMM Deep Field: Clockwise from upper-left, most-probable amplitude map, amplitude uncertainty map, and ratio map. The source detection acts as a statistically weighted low-pass filter. The funny pattern along the edges of both the most probable amplitude map and the ratio map is due to the averaging of only a few noisy values in the assignment of the pixel values near the edges. The saturated, mask sized regions in the upper right and lower left of the most probable amplitude map and the ratio map are due to very large valued pixels along the edge of the Subaru measurement map. . . . 101 List of Figures xvi 5.5 Subaru XMM Deep Field most probable amplitude pixel histograms with noise model. Top, histograms produced from pixels with noise less than.4mJy. Bottom, histograms made from pixels with noise less than 3mJy. There is an excess of positive, and negativefluxin the uncleaned most probable amplitude map, due to the presence of sources. The apparent lack of false detections in the top panel is due to the greater pixel noise populating the top histogram's positive tail. 103 5.6 Subaru XMM Deep Field source detections. The black crosses mark the location of detected sources. The closest detections are around 70" from the map edge 104 5.7 Subaru XMM-Deep Field source list correspondence. There are 26 unpaired UBC detections, 25 unpaired detections from the SHADES community, and 35 matches 105 5.8 Sky estimate map pixel standard deviation as a function of iteration. The standard deviation is within 0.01 mjy of its converged value at iteration 100 106 5.9 401 sky estimate map. The pixel-pixel correlations are apparent. st Measurements centred upon pixels involved in the centre of measurements less than thirty times are discarded. Source detection on this map will have to deal with the large scale structure, and if it is to be optimal, the cross-hatching 107 5.10 Deconvoluted source detection correspondence with the sources detected on the measurement map (1 iteration map) and with the st source list provided by the rest of the SHADES community using the data acquired up to February 1 , 2004. There is no indication which st source detection method is more reliable 109 xvii List of Figures 6.1 Simulation of data reduced using the central update algorithm. The data are simulated measurements of a statistically plausible fake sky with blank field extended noise samples. An extra chop's worth of data has been added to measurements made with the full SHADES measurement strategy. The chop has a size of 44" and is rotated 45 degrees with respect to the original chop directions. Clockwise from upper left, first iteration sky estimate map, 101 iteration sky st estimate map, 101 iteration power spectrum map, first iteration st power specrum map. The power along the dominant cross in the 101 st power spectrum map in Figure 4.14 is somewhat reduced, though not as much as in Figure 6.1, and the cross-hatch pattern modified. . . . 115 6.2 Simulation of data reduced using the central update algorithm. The data are simulated measurements of a statistically plausible fake sky with blank field extended noise samples. Two extra chop's worth of data has been added to measurements made with the full SHADES measurement strategy. The chops have a size of 44" are rotated 45 degrees with respect to the original chop directions and are mutually orthogonal. Clockwise from upper left, first iteration sky estimate map, 101 iteration sky estimate map, 101 iteration power st st spectrum map, first iteration power specrum map. Power along the prominant cross in Figures 6.1 and Figure 4.14 is further reduced. At low frequencies the power has been concentrated to narrow lines along the direction of the prominant cross 117 Acknowledgements xviii Acknowledgements I would like to acknowledge the people whose suggestions and insight were invaluable: Mark Halpern, Colin Borys, Douglas Scott, and Guillaume Patanchon. Part I Thesis Chapter 1. The Scuba Half-Degree Extragalactic Survey 2 Chapter 1 The Scuba Half-Degree Extragalactic Survey The Scuba HAlf Degree Extragalactic Survey (SHADES) is a major extragalactic survey of 850/xm radiation obtained using the Submillimetre Common- User Bolometer Array (SCUBA) at the James Clerk Maxwell Telescope (JCMT). The survey is designed to cover half of a square degree to a 4cr detection limit of 8 mJy 1 sources. By constraining the evolution, and the energy source of the population of galaxies detected at 850 /jm, the survey is designed to maximize the impact of JCMT on cosmology over the next three years. To do this, three specific questions are addressed. [7] 1. What is the cosmic history of massive dust-enshrouded star-formation activity? 2. Are SCUBA sources the progenitors of present-day massive ellipticals? 3. What fraction of SCUBA sources harbour a dust-obscured active galactic nucleus (AGN)? In addition to these three main goals it is desirable to maximize the legacy value of the survey. ^mJy = 10- 2 6 Wm- Hz 2 _ 1 Chapter 1. The Scuba Half-Degree Extragalactic Survey 3 1.1 Introduction To The Submillimeter Galaxy Population Submillimetre radiation ranges from 200/xm to about one millimetre. It is advantageous to observe the early universe at these wavelengths because of an effect called positive k-correction. Positive k-correction is the compensation for inverse square distance attenuation of submillimetre wavelength detected radiation by the increase in rest frame luminosity of the emitted radiation detected in the submillimetre as a function of redshift. For observations at 850 /im, the k-correction largely compensates for the inverse square attenuation of electromagnetic radiation and the detected flux from an eight hundred fifty micron source is approximately constant between redshifts two and eight. Similarly, because of the negative k-correction at wavelengths other than the submillimetre band, high-redshift, submillimetre detected galaxies often lack counterpart detections. [4] High-redshift submillimetre bright sources are believed to be very luminous, dusty galaxies which account for a significant percentage of all energy produced by all galaxies during the age of the universe. [4] Understanding the clustering and evolution of the high-redshift population of submillimetre detected galaxies will constrain early universe structure formation models. To constrain the clustering of the submillimetre population requires a large survey of submillimetre flux and modest, unbiased redshift estimates of detected sources. [7] Obtaining redshifts for the high redshift submillimetre detected galaxy population is difficult, as the inferred submillimetre energy distribution of the population is a modified blackbody spectrum, lacking features and having few detectable 2 spectral lines. The detected spectral lines consist mostly of CO emission, detections emanating predominantly from strongly gravitationally lensed galaxies. [4] If a galaxy detected at submillimetre wavelengths can be detected at other wavelengths, and the spectral emission distribution is known, it is possible to determine 2 "grey-body" Chapter 1. The Scuba Half-Degree Extragalactic Survey 4 where T is the temperature of the radiating grey-body, and z is its redshift.[4] Redshift estimates require temperature constraints. Methods of constraining the spectral energy distributions and the temperature of the submillimetre galactic population involve studying: 1. Theoretical models of galactic spectral energy distributions. 2. Galacticly lensed high redshift galaxies. 3. Low redshift submillimetre galaxies. 4. Low-redshift analogous galaxies, such as ULIRGs. Large, local galaxies with high star formation rates are bright in far infrared radiation, though stars themselves radiate negligibly in the far infrared. The far infrared radiation arises through absorption and re-radiation of the stellar radiation by heavy elements in the interstellar medium. The heavy elements are called "dust". Studies of the above constrain high redshift submillimetre galaxy spectral energy distributions. The distribution of the temperatures of ULIRGs might constrain the expected temperature for high-redshift submillimetre galaxies, allowing for redshift estimation. Regardless, [10] [3] perform extensive Monte-Carlo simulations and claim the feasibility of determining the redshift of sources detected in the SHADES survey, to Sz = ±0.5, given follow-up observations from the Balloon Borne Large-Aperture Submillimeter Telescope (BLAST). BLAST will observe SHADES fields at three additional wavelengths, 250, 350, adn 500/mum. A submillimetre survey spanning a large fractin of a square degree is required to constrain the clustering of the submillimetre galaxy population. [7] 1.2 What Is The Cosmic History Of Massive Dust-Enshrouded Star Formation Activity? Star formation activity is determined by extrapolating observations of local star forming galaxies, detected in the infrared[15], to submillimetre observations. The Chapter 1. The Scuba Half-Degree Extragalactic Survey 5 suggested star-formation rate is SFR = 2.1xlO~ x iU ) (1.2) LeOfim is determined by assuming a spectral energy distribution for the submillimetre detected sources and extrapolating to sixty micrometer luminosity. Note that small errors in the submillimetre source temperature can lead to large errors in Leofim because of the T dependence of blackbody luminosity on temperature. Nonetheless, 4 when combined with weak redshift constraints (6z = ±0.5), an estimate of the star formation rate as a function of redshift should be possible.[16] 1.3 Are S C U B A Sources The Progenitors Of Present-Day Massively Elliptical Galaxies? Currently, observations of high-redshift submillimetre galaxies are consistent with luminous, massive, dusty galaxies, similar to the local universe ULIRG, and also consistent with observations of moderate sized galaxies with optical and ultra-violet counterparts experiencing a period of rapid star formation. [22] One tenth of a percent of all local galaxies are ULIRGs, and if theformertype of galaxy is the predominant type of SCUBA source, then it is likely that they are the progenitors of today's massive elliptical galaxies. Massive elliptical galaxies trace the regions of over-density in the universe and experience strong clustering. Thus, if SCUBA sources are the progenitors of the massive ellipticals, they also, are expected to be highly clustered. 1.4 What Fraction Of S C U B A Sources Harbour A Dust-Obscured Active Galactic Nucleus? The spectral energy distribution (SED) of sources detected in the submillimetre is largely thermal and lacks characteristic features. The absence of sharp features in the SED makes it difficult to determine the source of radiation which is absorbed and thermally re-radiated by the dust. The possibilities are dust-obscured stars and Chapter 1. The Scuba Half-Degree Extragalactic Survey 6 the accretion disks surrounding dust-obscured black holes with masses greater than one billion solar masses. To detect active galatic nuclei (AGN), rest-frame soft x-rayobservations of submillimetre sources are required. 1.5 Legacy Value The legacy value of the SHADES survey is maximized by placing the data in a form maximally conducive to useful analyses. It is difficult to predict future use; the data need to be reduced in a careful, yet practical, fashion. Reduced data are compressed from the original data to an extent which strikes an acceptable compromise between information content and data quantity. The reduction process should introduce minimal, well understood spatial correlations to the data. 3 1.5.1 D a t a Compression The raw data consist of many measurements of the SHADES region, and are too numerous to be easily used. Reduction is careful removal of the time component, and careful interpolation of the non-uniform data samples of the flux emanating from the SHADES field onto a uniformly spaced grid of points. Each point is associated with a small surrounding region of space and is assigned a value. The region of space, point, and value combination is a pixel. A pixelized image is a histogrammed approximation to a two dimensional function which is the flux radiated by the region observed as a function of spatial coordinates. 1.5.2 D a t a C o r r e l a t i o n : Integrated Sachs-Wolfe Effect The accuracy of the reduced data image, or map, should be adequate and its precision quantified. The degree to which spatial correlations present in the data are due to the data reduction algorithm should also be understood. If this is the case, then it is possible to detect physical effects which result in spatial correlations. An example of such a detection in maps of the Cosmic Microwave Background (CMB), is the detection of the integrated Sachs-Wolfe effect, an effect detectable in 3 This is a statement regarding the quantitative description of the data certainty. Chapter 1. The Scuba Half-Degree Extragalactic Survey 7 CMB images spanning much larger angular scales than one-half degree. The SachsWolfe effect is the net gravitational doppler shift experienced by a photon as it traverses a region of high density. If the gravitational well expands at a constant rate during the photon's traverse then the doppler shift of the photon is the redshift the photon would experience during motion through flat, expanding space. If, however, the growth of the gravitational well changes during the traverse of the photon, then there is an anomalous doppler shift in addition to the redshift experienced by the photon due to constant universal expansion. By correlating the redshift's of photons radiating from the sphere of last scattering which pass through overdense regions of space, it is possible to significantly detect an acceleration to the universe. This 4 detection would be difficult, if not impossible, if the images of the Cosmic Microwave Background (CMB) were significantly spatially correlated by the algorithm used to reduce the CMB data. [12] 4 Q A > 0 , 9 5 % C L , statistical errors only Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 8 Chapter 2 Observing The Scuba HAlf-Degree Extragalactic Survey The SHADES observation strategy, executed on the JCMT using SCUBA, is designed to provide sufficient and uniform coverage of a half-degree square patch of sky. The use of acquired data requires calibration prior to reduction. The subsequent sections describe SCUBA, and the acquisition, pre-processing, and reduction of the data acquired during SHADES. 1 2.1 Submillimeter Common User Bolometer Array SCUBA is a sparsely-sampled array of hexagonaly spaced, cryogenic bolometers, situated behind filters whose pass-band is centred on 450 and 850 micrometer wavelength radiation. Radiation at these wavelengths are relatively unattenuated by the atmosphere. The array is located at the JCMT focal plane. Figure 1 shows the spacing of the 850/xm SCUBA bolometers. 1 T y p i c a l l y i n t h e a s t r o n o m i c a l c o m m u n i t y d a t a p r e - p r o c e s s i n g is r e g a r d e d as p a r t o f d a t a r e - d u c t i o n . I n t h i s t h e s i s t h e d i s t i n c t i o n b e t w e e n d a t a r e d u c t i o n , as t h e a c t o f r e d u c i n g t h e n u m b e r of p a r a m e t e r s w h i c h d e s c r i b e t h e d a t a , a n d p r e - p r o c e s s i n g w h i c h is c a l i b r a t i o n , is m a d e f o r ease o f discourse. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 9 8 5 0 Micron B o l o m e t e r Array -100 8 -150 o -200 h -250 -300 -150 -100 -50 Right Ascension (arcseconds) 0 Figure 2.1: SCUBA 850 micrometer bolometer array. The plot is obtained from the positions of actual data on the sky. Note the spacing between equal distances along the two axes are different, the height and width of the array are equal. 2.2 Data Acquisition The SHADES survey consists of 180 nights of observation in grade two weather, where grade one is best and grade six is the worst. Data are collected at 1 Hz by averaging the readings from two locations on the sky, which are observed alternately at a frequency of 7.825 Hz. The sky locations are varied by chopping, a process which refers to rapidly tilting the secondary mirror, or chopper. Every 16 seconds , the primary mirror is moved such that the 2 position "chopped onto" is the position originally pointed to by the telescope. The data acquired in this new position are subtracted from the "single-differenced" data previously acquired, prior to the primary mirror movement. The result is "doubledifferenced" data, which removes any offsets, and signals of constant slope along the 2 T h e r e is a t h r e e s e c o n d l a g w h i l e t h e p r i m a r y m i r r o r is m o v e d . Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 10 direction of the chops, aiding in the removal of signal from atmospheric emission. Each sample from the double-differenced data is the normalized subtraction of two 3 single-difference samples acquired 19 seconds apart. The focal plane of the JCMT is undersampled by SCUBA both at the 450 pm and 850 fxm wavelengths. Thus, a small, 1 Hz motion of the secondary mirror, called "jiggling" is used to acquire data located between the bolometer detectors in the focal plane, and "fully sample the array" . In particular, the jiggle pattern is a 4 16 point pattern repeated every nineteen seconds designed to Nyquist sample the 850 detectors. A 64 point pattern, which repeats every four primary mirror movements, is used to nyquist sample the 450 //m detectors. The difference in the number of points in the jiggle patterns for the two wavelengths is due to the narrowing of the telescope point spread function from the 850 /zm wavelength to the 450 fim wavelength, ie. The spatial Nyquist frequency is increased for detection of 450 /jm radiation. The jiggle pattern is specified in elevation and azimuth, and also rotates with respect to the celestial sphere. Observation of each SHADES pointing involves a 64 point jiggle pattern with three chop "throws" specified in two directions on the celestial sphere. The chop 5 throws are 30", 44", and 68" in the east-west direction and 30", 44", and 68" in 6 the north-south direction. The chop throws are specified in declination and right ascension. These chop throws are chosen such that each chop throw measures a different spatial frequency harmonic and links many different pixels in different measurements. More on this topic will be discussed in section 6.3 on measurement strategy. The SHADES field is covered by many pointings of the type described above. A relatively uniform observation coverage of the field is obtained through "tripos" pointings, a scheme whereby three pointings, whose pointing centres are offset to form three overlapping circles of sky coverage, or a tripos, are panned over the SHADES region. Figure 2.2 is a cartoon sketch of the coverage of a single tripos 3 4 5 6 The double-difference is divided by two. Nyquist sample. Right Ascension/Declination (J2000 coordinates). 1" is one arcsecond. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 11 pointing. [17] Each circle in the tripos pointing is actually hexagonal and spans the 2 Figure 2.2: Tripos Coverage array size, roughly 150". As previously described, each pointing in the tripos consists of the six chop throws in the two directions, centred in the centre of one of the circles in Figure 2.2. 2.3 S C U B A / J C M T Measurement Noise The measurements acquired by the SCUBA detector are noisy. The precision of the reduced data depends upon the noise present in the unreduced measurements as well as how the reduction process affects the measurement noise. Much of the design of JCMT and SCUBA is motivated by a desire to mitigate the effects of noise. 2.3.1 Brief Overview O f Noise The measurement of a physical quantity is composed of the quantity of interest plus one or more contaminating or nuisance signals associated with the measurement process. In principle, the nuisance signals can be computed based upon classical physics and removed. In practice, a physical description of the nuisance signals 7 7 T h i s discussion assumes a l l measurements c a n be described b y classical physics. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 12 is very difficult and probablilistic descriptions of these signals are employed. A nuisance signal described by probabilistic properties is called noise. The spectrum of typical measurement noise is a combination of two processes 8 which are described by a particular function of frequency. The two processes are " j " noise and "white" noise, j exhibits a monotonic decrease with increasing frequency and white noise is constant over all frequencies, j noise dominates at low frequencies and white noise dominates at frequencies greater than the " j knee", or the frequency at which the white and j noise components to the noise spectrum are equal, j noise indicates power at low frequencies and corresponds to correlations in the noise. A typical assumption made of the two processes comprising the time stream spectrum are that they are Gaussian distributed. That is, the random variables comprising the process are described by Gaussian probability density functions. This approximation, when describing noise, is supported by the Central Limit Theorem. The Central Limit Theorem states that probability density functions of the sum of many random variables, regardless of their density functions, tends to a Gaussian probability density function as the number of random variables in the sum tends to infinity. [14] This claim is consistent with physical effects causing the noise; that is, the random nature of the noise is due to complexity. Alternative justification for the use of Gaussian probability density functions comes from a maximum entropy arguement, whose result is that a Gaussian probability density function assumes the least amount of knowledge about a process under certain, commonly encountered, conditions. [20] Samples from a random process can be dependent. Probability density func9 tions for the random variables comprising a process with dependent samples are multivariate and non-separable. Data sampled from such a random process results 10 in correlated samples. When detecting signals in a noisy data sequence containing correlations, it is necessary to account for the correlations. An optimal method of fitting to data with noise sampled from a Gaussian process is given in section 2.5 on 8 A random process is a collection of random variables[14]. When discussing a sequence of data points, the random process is the entire sequence and each data point is a random variable. Collection of random variables. A multivariate probability density function is called a joint-probability density function. 9 1 0 Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 13 data reduction. Gaussian probability density functions of a single variable are completely characterized by their mean value and their variance. With a multivariate Gaussian random variable, the process is completely described by its mean and its covariance matrix. The diagonal entries of the covariance matrix give the variances of the variables, and the off-diagonal terms give a measure of the amount with which each variable varies with the others. For a sequence of data values sampled from a random process, given by the vector n, the covariance matrix is N (2.1) = (nn ) T and the correlation matrix is N~ . Note that N is symmetric and has dimension l equal to the length of the vecotr n. Here, each data sample is considered a random variable in the Gaussian random process which is sampled to produce the vector of data samples, ft. When attempting to detect weak signals or large signals to high precision, one performs a low signal to noise ratio experiment. In such experiments, effort is made to reduce the noise correlations such that averaging TV measurements results in a reduction of the standard deviation of the averaged measurements by a factor of VN. 11 If all noise correlations have been removed from the signal, an arbitrar- ily large signal to noise ratio can be obtained given a sufficiently large number of measurements. A noise signal possessing correlations, such as a slowly varying drift, cannot be averaged to achieve an arbitrarily high signal to noise ratio because the average of the slow drift will eventually dominate the error. The drift appears as j noise in the power spectrum of the signal. 2.3.2 SCUBA Noise The dominant noise present at the SCUBA detector comes from atmospheric emission. The remaining sources of noise are comprised of confusion noise, photon detection statistical noise, Johnson noise, noise due to non-ideal electronic components n T h e variances of Gaussian random variables add as their inverse square. Thus, the variance of the sum of N measurements is a factor of y/N larger if all of the measurements are sampled from the same probability density function. Division by N to obtain the average value reduces the standard deviation of the average of the measurements by a factor i V ~ ' . 0 5 Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 14 (ie. physical components), quantization noise, etc. Combined, the effect of the noise is to give j noise and Gaussian white noise components to the signal. Atmospheric emission, to good approximation, is constant over the patch of sky subtended by the SCUBA bolometer array.[1] Astronomical observation of the SHADES region is a low signal to noise ratio experiment. The instantaneous signal to noise ratio is approximately 3% for a 10 mjy source, and considerable averaging is required to detect SHADES sources at the desired 4cr level, where u is the standard deviation of the beam convolved pixels, and is desired to be 2.5mJy.[17] 2.4 Data Preprocessing Prior to data reduction the data must be calibrated, outliers rejected, and nuisance signals understood, estimated, and if possible, removed. 2.4.1 Calibration SCUBA measurements are calibrated to account for unequal bolometer gains, atmospheric absorption and to acquire the correct conversion from detector signal in Volts to detected flux in milliJanskies. 12 The SHADES survey is granted time in grade two and upper grade three weather. Weather, in this context, refers to the current opacity of the atmosphere to 850 /um radiation. The weather is quoted in terms of a fitted parameter in a model of atmospheric absorption. The model is / = I e- TA 0 (2.2) where IQ is the intensity of 850/im radiation incident on the atmosphere, / the intensity of 850 //m radiation incident on the focal plane, r the 850 /an zenith opacity, and A the airmass. 13 Grade two to upper grade three weather corresponds to a zenith tau measured by the CSO lJansky = 1 0 12 1 3 - 2 6 Wm^Hz- 14 at 225 GHz between 0.05 and 0.10.[17] The 1 Airmass is inversely proportional to elevation and incorporates the increased amount of atmo- sphere along the line of site into equation 2.2. 14 Caltech Submillimeter Observatory Chapter 2. Observing The. Scuba HAlf-Degree Extragalactic Survey 15 225 GHz zenith tau is related to the 850 /im tau through relations available on the JCMT webpage. The relations are determined by correlating the CSO and have changed with time. The data are adjusted by the JCMT standard SURF procedures, flatfield and extinction, correcting for unequal bolometer gains and atmospheric absorption, respectively. The flux conversion factor, or FCF, in units of Janskies per Volt, is computed by reducing calibrators, fitting a two dimensional Gaussian to the resulting map of the calibrator, and comparing the integrated flux within the full-width half maximum of the fitted Gaussian to a known value. For a given night of observation there are several calibration observations, and the FCF computed from these observations is not constant. Figure 2.4 is a plot of the FCF multiplying the data as a function of data observation file. These FCFs are all computed using CRL618 using Method 3, described below. There is no discrimination between FCFs computed using calibration observations taken with different sized chop throws. The observation files contain measurements of the Subaru XMM-Deep Field from December 18, 2002 to June 22, 2003. There are 292 observation files each containing 320 double-difference measurements for each of the 37 bolometers in the SCUBA array. The observationfilesare not in time order. Deciding which FCF to apply to a given observation file of data is complicated. The following methods were tried: 1. Set the FCF to a nominal value of 225 Jy V . -1 2. Average the first and last computed FCF for a night of observation. Use a running average lOOJyV- . 1 of the FCFs if the FCF is above 400 Jy V - 1 or below 16 3. Average the FCFs computed from calibration observations performed most recently with respect to the time of the observation to be calibrated, and 1 5 A running average is an average of every FCF that is actually used on the data during the data reduction up to the current point in the program execution. The running average is re-computed with every computed F C F that will be used in the conversion of the data from V to Jy. 16 FCFs above 400 Jy V - 1 and below 100 Jy V - 1 are deemed bad and neglected from the com- putation. 400 and 100 are chosen as many standard deviations from the average. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 16 Dote? ^ Si~^!E ' Q l > $-s' 1 ' e e n £ rn J 3.3*11? j-o-ia 5, LSx^a* a.'B«iD* Dose Sequence: Bolometer 20 Bctio Sequence; Bolometer 1ft 3.a.fio* i.Qwin* i.a«io> i.sma* i.a*riffi^' s.c^cfl a.c*io* i.c^ic^ i.a-oa* 3.3*!ry ijQ*'H»* i.iiaa* Doioi Sequence; Bnlorn&tef 21 U H Dalo Sequence: Solomcicr 24 Qai« Sequence; Snipmeler 2.S I.O-IB* !.i»'ia* a.o*!D* 1,0**0* s.ax'.a 1 IJO-WIS* Sg^clryrn Mgqnij^rJB^ o.-oooi C-JMIQ aenco O.CDtn o.toan a.floai DiODIO g^ornftlf;r 19 0.01430 Spectrum Moqniiiide: Bdorr'Cis 20 r Spectrum Mcqniiudc: BoJorntjier 19 1 i'.e^ff j | p vmmt arena Sfj*u;tfum aomo a.oiflo iySaqnity«»i-; a.rjico c.iooc D-OS10 UDDQ 0.0 5 CO 0.! COEJ Sp^ryrr? MpqnlUi^p: gi^omsfor ?.Z EM^metsr 2'f C.DI'BO o.irxiC' 0.5030 Sp^-run-i v.rjcniiuiic: 3<srigrrBrtcr £3 i»-*E , . D.ocia . ' o.o'.co Q.IDOS to-4. . &OOIQ . 0.0130 . . G.IGOD ^3CO0 Figure 2.3: Eight representative scuba bolometer timestreams and spectra for the Subaru XMM Deep Field after flatfielding, and extinction correction. The timestreams are non-stationary and strongly correlated by atmospheric emission. The timestream for bolometer 21 extends past the range of the y-axis. The standard deviation of the bolometer timestreams varies substantially from bolometer to bolometer. The magnitude of the Fourier transform of the bolometer timestreams are shown in the lower half of the plot. Note the flat spectra and the spikes on bolometers 17, 19, 20, 21, and possibly 23 and 24. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey Flux C o n v e r s i o n F a c t o r vs Observation 17 Number 300 250 200 150 100 200 Observation Number 300 Figure 2.4: Flux conversion factor as a function of observation file number. The observationfilesare not in time order. Note the considerable fluctuation of the FCF with observation number. Only calibrations with CRL618, Uranus and Mars are used. The chop throw of the calibrations is ignored. the soonest to be performed after the observation to be calibrated has been completed. If there are only FCFs before or after the observation, use the FCF computed from a calibration observation that occurred closest in time to the observation to be calibrated. Use only FCFs which were computed after the most recent focus. Do not use any FCFs above 400 J y V - 1 or below 100 J y V . If under these conditions there are no FCFs, use the running -1 average. The locations and number of 3cr sources detected in a measurement map of the 17 Subaru XMM -Deep Field demonstrated a slight dependence upon the method of FCF computation. To quantify the effect, the sources detected in two maps of the Subaru XMM Deep Field are plotted in Figure 2.5. One of the maps is produced 17 A measurement map is a non-deconvolved map obtained by grouping all of the measurements centred upon a pixel, averaging them, and assigning the average to the pixel value. This is done for all pixels of the map. Measurement maps are described more fully in section 2.5.1 on Interpolative Reduction. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 18 using measurements calibrated with FCFs computed using the third method above, the other with measurements calibrated using a nominal FCF value of 225 JyV . -1 The source detection algorithm is described in section 5.2. The correspondence between sources detected from data reduced with Method 1 calibration and sources detected from data reduced with Method 3 calibration is 90%. Here correspondence referes only to location and not to the significance of the detection. For the remainder of this thesis, data reduction is performed with an FCF set to 225 Jy V . -1 Despite the correspondence demonstrated in 2.5, the usefulness of the calibrations depends on the degree of fluctuation of the FCF values, as well as on the intended analyses. For instance, the effect of calibration error on clustering computations is unknown and not constrained by this test. If it were desirable to estimate the flux to greater precision, and the statistical error were not dominant, then accounting for thefluctuationsin the FCF would be necessary. To estimate the effect of FCF fluctuation on the analysis intended, realistic simulations must be performed. Because of the quantity of calibration data taken during the SHADES observations, the statistical properties of the FCF can be determined and used in subsequent simulation. 18 2.4.2 Nuisance Signals And Outlier Rejection For SCUBA data taken by double differencing, the dominant nuisance signal is atmospheric emission. Atmospheric emission is approximately constant for all bolometers on any bolometer array sample. [1] Atmospheric emission and outliers present in the data are removed iteratively by removing 5cr outliers from the timestreams, then, for every array sample, subtracting the average of the bolometer samples, and repeating until convergence. [19] This process is further complicated by the presence of spikes at 1/16 Hz in the magnitude of some of the bolometer timestream Fourier transforms. 19 When subtracting the average of the array, if there is a coherent component to 18 19 T h i s idea was communicated by Tim Jenness to the SHADES collaboration through email. In this context, units of Hertz are not completely accurate, due to the delays in data acquisition between telescope nods. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey ~i 1 1 r - i — i — i — i i - r | i i i r 19 - + Calibrated O Uncalibrated 2h o c 73 Q O -2 _i -10 -5 i i_ 0 Right Ascension (arcmin) 5 Figure 2.5: Correspondence between sources detected at greater than 3a significance from calibrated data and from uncalibrated data. The calibrated data uses Method 3 in Section 2.4.1 to compute flux conversion factors. The source detection is described in Section 5.2, but differs in that the detected sources have not been iteratively cleaned from the map. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 20 the spike, then the value subtracted from the bolometer timestreams will contain the coherent average of the spiky bolometer values. To avoid this, when subtracting the average of the array at every time step from the individual bolometer readings, the array is divided into two parts. One part of the array is composed of bolometers which do not exhibit the spectrum spike, and the other part of the array is composed of bolometers which do. The spike is detected by banding the magnitude of the Fourier transform of each bolometer timestream into three equally sized bands in frequency. There are no gaps between the bands of frequencies, and the central band is centred upon ^ Hz. Within each band, the sum of the Fourier transform magnitude is performed over the band frequencies. The sums for the first frequency band and the last frequency band are averaged and compared to the sum for the central frequency band. If the central band sum is greater by a threshold, then the spike is considered detected. The average of the bolometers in the different axray parts is computed and subtracted from the bolometers which belong to the part of the array from which the subtracted average is computed. When computing the array averages, bolometers exhibiting a standard deviation three standard deviations away from the mean of the bolometer standard deviations for the observation are not included in the computation. I will refer to this computation as array average subtraction for the remainder of this thesis. The frequencies at which astrophysical sources appear in the data are at the harmonics of the 64 point jiggle pattern. That is, at l/64Hz, 1/32 Hz, 1/16 Hz, etc. The harmonics are formed from the repetition of the jiggle pattern every 64 seconds. Note that in the absence of the jiggle pattern, when pointed at a source, the chopping and nodding double-difference measurement strategy will produce a constant signal, which appears at zero frequency in the double-differenced data. Additional Fourier space structure will be caused as the bolometers move on and off source due to the jiggle pattern. This is expected to broaden the spectrum centred on the 1/64 Hz harmonics. Thus, bandstop filtering, with the stopband centred upon the 1/16 Hz spike is not an option. Lastly, the array average subtraction is performed iteratively along with 5a outlier rejection. That is, the array average subtraction is performed, then measure- Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 21 meritsfivestandard deviations from the mean of the measurements for any given observation file (320 double-difference measurements per observationfile)are removed. These two computations are repeated four times. Figure 2.5 is a plot of the clean timestreams and the magnitude of the Fourier transform for a representative sample of the bolometers. The preprocessed timestreams have a smaller standard deviation than their non-preprocessed counterparts and have no outliers. The magnitude of the Fourier transform of the preprocessed bolometer timestreams is flatter than the non-preprocessed bolometer timestream Fourier transform magnitudes. (Figure 2.3) 2.5 Data Reduction The final step in obtaining useable data is data reduction. The desired properties of the reduced data are stated in section 1.4. This section describes common methods of data reduction, their strengths and weaknesses, and relevant issues. It is beneficial when discussing the methods of data reduction to describe the process using linear algebra. To this end, beginning with the single-difference measurement, the telescope measurements are described as equations relating measurements to linear combinations of the observed sky pixels. When doing this, a useful approximation is to approximate the telescope point-spread function as a Diracdelta function. A double-difference measurement becomes a relationship between 20 three pixels, and a single-difference measurement is given by: M = - +pi lc Po (2.3) where po is the pixel in the chop location of the single difference, pi is the pixel at which the telescope points, and M\ is the left-chop single-difference measurement. c A second single-difference measurement is formed after the primary mirror is moved: M = - +p TC Pl 2 (2.4) Now the chop is looking at the position first pointed to by the telescope, and the telescope points at a new pixel, p2- The double-difference measurement is obtained 20 T h i s approximation can be relaxed by maintaining the Dirac-delta telescope point-spread func- tion, and convolving the true sky map with the actual telescope point-spread function. This is done in simulations later on. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey Quia Sei^yssnce: Botorne!;*? 1 Dfli« Sequence: BnlnmMer M i a.a-fcin' 0 J,Q*IB* j i.a*io* i.3*iD j 1 a.fl.«iD* u*!c J z.a*'*(P i.s*ia» s J !.Q"IC5 I.&K'.CJ* Z.&xifX- [Join _ Secuflncc: , Bulnrncic ,—. , Sflio Haoycnce; EJoipmeter *9 ;.L _ _ 22 S.C-'D' ;jO»lC ft ,— i.5*?C* , 2.rj«tC Dale Sflquflnc*;; Bolometer 22 Daio Sequence: Boiofnclef ? 1 <##» S.a-lD^ '.0**2* ^ft-ilO* 3.3*10* fi-DS'ioa to-*! - i.a».<c 3 1 1 ?.v*iiP LSMO* B.OO ! c a.sicra S^ctrM^^ *PgrtJfff? ?l O.CK!3 5.0* 1 C Spcctfurrt Msijoriuar:; B r i o ^ ' " i P L J CUB 1C c: {-C« 1 C * J .3*50* '2 Sir-1 CP _ SpcciHjm IVo^jiiJijde: B^yrrisigi JB 1 !}-*-£. o.scei L5«1tfi Da In 5eny«nco; Bolfifncis? 24 Dai-o Sentience: ti^'or'"*?^^ 23 :iL LO*tC* &. was fT^II!li!ll|p|(|||j 10" OJKJt C.OaiC Q.CIOD 5.WCO o.aooi oxflic o.aiao ciaco arafo cunon MOCC Hctlofneicf 10 0.0 < ' OB 3 a CO Stwc^rn Mcgniiijdi?; Bot&mdsr #1 1B~ O.OCQ! D.OJM 0 fl.'B f CM C, 10GO Spflctru^t MrjgfiHugrr: Bpiorrysiftr 7J- f>. 00.10 0.2103 Figure 2.6: Eight bolometer timestreams and their corresponding Fourier transform magnitudes for the Subaru XMM Deep Field after preprocessing. Note the lack of outliers in the preprocessed timestreams as well as the reduced rms when compared to the timestreams in Figure 2.3. The Fourier transform magnitudes are flatter than their non-preprocessed counterparts. (Figure 2.3) Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 23 by subtracting the two single- difference measurements and dividing by 2, obtaining: M Mdd i c = - M -po + 2pi - p r c g = 2 2 . ^ . ' Note that in actuality, po is the telescope point-spread function weighted sum of sky flux centred upon the centre of pixel po. The above approximation, using the Dirac-delta telescope point-spread function, is one which introduces a small amount of error into reduced data maps, and is used in all data reduction documented in this thesis. [19] Lastly, realistic measurements contain noise: (2.6) where n is the noise. Some methods of reducing observed data are: 1. Interpolative Reduction 2. Emerson Reconstruction^] 3. Direct Inversion 4. Iterative Inversion Methods one through four will be discussed in the following sections. 2.5.1 Interpolative Reduction Interpolative reduction is performed by fitting a function to measurements of the sky within some region of the pixel whose value is to be determined. The value of the fitted function at the pixel centre is assigned to the pixel. Fitting and pixel assignment are repeated for every pixel in the map. Commonly used fitting functions are two-dimensional Gaussians, planes, and step functions. Note that the act of fitting to a step function which is pixel shaped and non-zero only inside the pixel of interest is the same as averaging all measurements of the sky located within the pixel of interest. To see this consider x: 2 N (2.7) i=i Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey where y is the pixel value, or the height of the step function, di is the i 24 th data sample located within the pixel, and N is the number of data samples within the pixel. Taking the derivative of x and setting it to zero, obtain: 2 N N = (2-8) i=l i=l and (2-9) v = which gives the pixel value as the average of the N measurements lying within the pixel. Maps created using the above technique will be called "measurement maps" for the duration of the thesis. Figure 2.6 is a measurement map created by reducing simulated measurements of a source. The simulated measurements are taken of a fake sky consisting of a single source, using the SHADES observation strategy for a single SHADES pointing. There is an equivalent and informative, linear algebraic formulation of interpolative data reduction. The linear system, Px = M (2.10) where P is an N measurements by N pointings matrix, x is an iVp pixel vector of m p pixel values and M is an N vector of measurements. Equation 2.7 describes noise m free measurement of the sky estimate, x, to get the data, M. In the case of data reduction, M is known and the task is to determine x. This is the inverse problem. Since the linear system is overdetermined and the measurements noisy, an optimal map, x is desired. The inverse solution, optimal in a generalized least-squares sense, is x = (P N- P)~ P N- M T 1 1 T (2.11) 1 Where N is the noise covariance matrix defined in equation 2.1. The element of the N by N noise covariance matrix for the i th m m row and j th column is given by the covariance of nj with rij where nj is the noise added to the i measurement.[9] th With the correct choice of P, the least-squares solution, x, given by equation 2.6, is equivalent to determining the interpolated pixel values. Consider the case Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 25 Figure 2.7: A measurement map of simulated measurements taken with the single pointing SHADES observation strategy of a fake sky with a single, centred source. Note the presence of negative sidelobes centred upon the six chop throws surrounding the source. Each square pixel in this map is three arcseconds square. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 26 where the interpolating function is a two-dimensional step function, zero everywhere, except for the region of the pixel, where its value is one. Let the noise be uncorrelated and have unit variance. Then N is equal to the identity matrix. Further, let P be composed of rows that are all zero except for a single column in the each row, which is one. Now, consider the j h row of Px. The location of the one in the j h r o w t t of P is such that upon performing the multiplication, Px, the one in the row of columns multiplies the pixel upon which the measurement comprising the j h element of M t is centred. More concisely, ( )j PS mod N = 3 mod N = Mj (2.12) X p p Thus, P x forms a vector of length N whose i h entry is the sum of all measurements T p t centred upon the i h pixel, and (P P)~ T is an N by N diagonal matrix whose x p t p entries are -»r—, where -/Vbs is the number of times each pixel is centred upon in 0 a measurement. Thus, the assigned pixel values are averages of the measurements centred upon each pixel. If we generalize to uncorrelated noise with some variance, the noise covariance matrix becomes diagonal but not the identity, and the resulting pixels are the inverse-variance weighted sum of the measurements which were taken when the telescope was pointing at each pixel. It is now desirable to determine the uncertainty of the pixel values when the data are noisy. Let M = Ps + n (2.13) where s is the flux emanating from the sky sampled at the pixel centres, M is the pointing matrix of the above section, and n is the vector of noise samples of length N . m Every nf noise sample comes from a probability density, which is 1 in general multivariate, and non-Gaussian . Assuming that the noise is sampled from a Gaussian process with zero mean, then the noise probability distribution is completely specified by the noise covariance matrix, as previously mentioned in section 2.3.1. If N is diagonal (uncorrelated noise), and the noise distributed with a variance, a , then the pixel-pixel noise covariance matrix is given by 2 E = a (P P)~ 2 T 1 (2.14) Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 27 In general, the pixel-pixel covariance matrix is given by (xx ) = ss + £ T (2.15) T where the pixel-pixel noise covariance matrix, E, (2.16) £ = (M N~ M)~ T l l Here £ gives the uncertainty of the pixel parameters. 21 2.5.2 Emerson Reconstruction Emerson reconstruction^] is a method of removing the negative side-lobes around sources from a measurement map formed from differential measurements. The method uses Fourier techniques and requires the size and direction of the chop throws to be constant for all observations. The measurement map, x, can be considered as a two-dimensional convolution of a measurement function with the sky flux: x = c®s + n (2.17) where c is the convolution function and n is the noise. By inverse Fourier transforming the quotient of x and c one obtains the deconvolved map. p= s + T~ (^j 1 (2.18) where n is the temporal Fourier transform of the noise and c is the Fourier transform of the convolution function. In the pencil point spread function approximation for a one-dimensional map, c = -0.55 (x-x 0 Ax) + 6(x-x ) 0 - 0.56 (x - x + Ax) 0 (2.19) where x is the position pointed at by the telescope and Ax is the chop size. Thus, 0 c = exp ikx [exp ikx — 1 cos kAx\ 0 2 1 See, solution. (2.20) f o r e x a m p l e , [20] f o r t h e c o n n e c t i o n b e t w e e n l i n e a r r e g r e s s i o n a n d t h e m a x i m a l l y likely Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 28 and zeros in the magnitude of c cause spikes in j, introducing sinusoids into the deconvolved noise, correlating the pixels. This method of obtaining deconvolved reduced-data is not considered further in this thesis. For a full critique of this method see [11] [5] 2.5.3 D i r e c t Inversion Direct inversion is solving for x given M in equation 2.10 by brute force, where P, x, and M are given in the interpolative reduction example, but with a pointing matrix P, given by the actual measurement strategy. That is, the ,7 row of P has th a —0.5 in each column multiplying the chop pixels in the j measurement, and a 1 th measurement. in each column multiplying the pixel centred upon in the j t h There are two problems with this approach: 1. P P is singular. T 2. P P is huge. T The unmeasured modes of the sky manifest themselves as singularities in P P. T The inversion of P P can be computed, up to a linear combination of the singular T eigenvectors of P P by computing the pseudo-inverse of P P. [21] A pseudo-inverse T T is an inverse of the space spanned by the matrix, and yields the inverse of P up to a linear combination of the eigenvectors spanning the null space. The second problem is more serious. P P is N rows by N columns and T p p contains iV elements. For three arcsecond pixels, SHADES consists of 360000 2 pixels. 22 Thus, the computation involves the inversion of a 1.3 x 10 element 11 matrix. If each element is represented by eight bytes, the matrix is stored in 1000 Gigabytes of memory. This is 10 100 GByte hard drives. Standard matrix inversion techniques require order N operations to perform, in our case, this is roughly 10 3 33 floating point operations. On a GigaHertz machine, assuming 10 floatingpoint 9 operations per second, one completes the inversion in approximately 10 seconds, 24 or in 5.5 • 10 22 23 16 months. Exploitation of symmetries in the inversion may reduce 23 T h e survey is 1800 axcseconds by 1800 axcseconds in extent. Millions of Hubble times. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 29 the number of operations required by a factor, but will not reduce the O (N ) 3 dependence.[13] 24 Supercomputers are an option, but are not well suited to serial computations such as matrix inversion. They rely upon performing many parallel and independent computations simulateneously. The time required to perform the brute-force inversion is prohibitive, and due to the existence of a faster, iterative approach, it is not pursued. 2.5.4 Iterative Inversion In the mid-1990's Cosmic Microwave Background (CMB) experiments became large surveys of the sky. The number of pixel parameters to be determined in these surveys caused experimentalists to look for another way to solve equation 2.8. Ned Wright developed an iterative algorithm for determining the best fit pixel values which requires 0(N ) operations, and the idea was adopted and extended by others.[23] m The algorithm uses the data in timestream order, eliminating the need to store large matrices. For a single-difference experiment, the type of experiment in which the algorithm has been used, the iterative update scheme is given by ,j=l (2.21) where xf is the i pixel of map x at the n iteration, N is the number of measureih th m ments, j is the measurement index, m,j is the j th of the j t h measurement, <r is the variance 2 measurement, and A is the chop distance in pixels. [5] In words, when updating a pixel, assume the rest of the pixels comprising the map are correct. For each measurement involving the pixel to be updated, compute what the pixel must be if the other pixel in the measurement, and the measurement itself are correct. Thus, there will be a "correct" pixel update value for each measurement involving the pixel to be updated. Compute a weighted average of the "correct" pixel update values and assign the weighted average to the updated pixel. The weights are the variance of the j th 24 measurement. This is Wright's algorithm generalized to noisy N o t completely true; circulant matrices are at least one exception. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 30 measurements. A hallmark of this algorithm is that the next iteration pixel value 25 does not depend on the current pixel value. Note that on any one iteration there are conflicting assumptions. The update of a pixel requires that all of the other pixels are correct. Thus, when updating another pixel, the pixel just updated is assumed correct, and one of the pixels previously assumed correct, is now assumed incorrect. It is thought that this is why iteration is necessary, and that iteration will converge to the correct answer, a map optimal in the generalized least-squares sense. 26 Prior to this work, the algorithm described by equation 2.21 has been successfully employed in reducing the data from many CMB experiments. Most notably, the algorithm has been successfully used in the reduction of Wilkinson Microwave Anisotropy Probe (WMAP) data. [9] The success of the algorithm has inspired its use in the reduction of data acquired from the SHADES field using SCUBA. The application of this algorithm to SHADES data comprises the main topic of this thesis. 25 26 T h e algorithm in this form does not take into account correlations between the data samples. Provided the timestream noise covariance matrix is diagonal. Chapter 3. Iterative Reduction of SHADES Data 31 Chapter 3 Iterative Reduction of SHADES Data In this chapter the iterative reduction algorithm, equation 2.21, is generalized for use with an arbitrary measurement strategy. The generalized algorithm is studied theoretically and by simulation. 3.1 Wright's Algorithm Generalized Equation 2.21 is valid for single-difference data and must be modified for doubledifference data prior to its use in SHADES data reduction. The generalized algorithm retains the essence of the original algorithm. For each pixel involved in any one measurement, assume all of the remaining pixels in the measurement are correct. Compute the value of the pixel under consideration which will satisfy the equation between the remaining pixels , which are assumed correct, and the value of the measurement. Repeat this for all pixels involved in all measurements. For every pixel in the updated map, the above procedure will yield many estimates for each pixel. Averaging the estimates for each pixel gives the updated map. This is the same procedure used to determine equation 2.21, except now there are more pixels involved in a single measurement. The generalized version of equation 2.21 is: (3.1) Equation 3.1 gives the n th update of the i number of measurements, j is the j t h pixel in map x, where N m is the element of the measurement vector M, and is the gain or factor multiplying the i pixel in the computation (acquisition) of the th Chapter 3. Iterative Reduction of SHADES Data (AP) th 32 measurement. The denominator is a sum of all gains, in all measurements, multiplying the i t h pixel. In this form, the iterative deconvolution update can be used to reduce data taken with an arbitrary point-spread function, and an arbitrary measurement strategy. For the reduction of SHADES data, assuming Dirac-delta telescope point-spread function, the gains are - 0 . 5 for gains multiplying the pixels at the chopped positions and 1 for gains multiplying the pixel upon which the measurement is centred. This is the double-difference measurement described in Section 2.6: M = -0.5pieftchop +Pcentre ~ 0 . 5 p i t chop r 3.1.1 (3-2) gn A l g o r i t h m Checks There are two easy checks following from the considerations cited below. 1. The update of any pixel does not depend upon the previous value of the pixel being updated. This is a direct result of the assumptions made in the formation of equation 2.21 2. The limit of the pixel update should be zero, as the gain tends to zero for any single measurement. That is, when the pixel being updated is not involved in a measurement, the algorithm should be well-behaved. Check #1 Consider the update given by a single measurement, ie. N m update for the i th = 1, and consider the pixel. The update is given by s? = ^ 2 - " ' +x?- (3-3) X \9i) The gains in the numerator and denominator cancel when multiplying the x™~ 1 term in the numerator of the update, with the result that the trailing x™~ term is 1 cancelled by the in the first term of the update. Since this occurrs for every measurement, the update of any pixel does not depend on the previous estimate of itself. Chapter 3. Iterative Reduction of SHADES Data 33 Check #2 Consider the contribution to the update of the x\ pixel from a single M- in a h dataset with N measurements, where the j measurement does not involve the th m 7 pixel being updated. That is, consider the limit of equation 3.1 as gj tends to zero when g\ / 0, for all k ^ j: » 9 i m x n = & ^ A m * r -° + x > r ( 3 . 4 ) (si) Breaking the first term into two fractions, obtain, B ^° m ^ Ef5 (si) — ( Ef", 3 . 5 ) (si)' When k = i in the second term, the second term will cancel the right-most term. Explicitly performing the cancellation yields, lim *? = - A T , T ^ ~ 2 (3-6) which equals zero. 3.2 Generalized Algorithm Simulations To characterize the generalized iterative algorithm, simulations of data reduction using equation 3.1 are performed on data acquired during mock experiments. Each mock experiment yields measurements of a fake sky taken according to a measurement strategy. Fake skies measured vary from simple, seven bin, one dimensional 1 discrete skies with a single Kronecker-delta source, to two dimensional discrete skies populated with Gaussian convolved point sources sampled from a realistic source count model. Like the fake skies, simulated measurement strategies range from 2 simple, single-chop, single-difference measurement strategies, to the relatively complicated SHADES measurement strategy. 1 2 A fake sky is a map of the sky based upon some model of astrophysical flux. See [6] Chapter 3. Iterative Reduction of SHADES Data 34 In each simulation an initial map estimate is made, and equation 3.1 iteratively applied. Convergence is practically determined to be when the change in the sum of square pixel values is below a particular threshold. 3.2.1 SHADES Simulation The full SHADES measurement strategy is applied to a half-degree square fake sky populated with JCMT point-spread function convolved point sources sampled from a realistic submillimetre source count model. The original map estimate is taken 3 to be the null map. Figure 3.1 shows the fake sky map, and the residual maps 4 for various iterations. The colour bax at the bottom of the Figure applies to all of the maps except the fake sky map, which saturates at pixel values of zero and ten. A residual map is the difference between the estimate of the fake sky map and the actual fake sky map. There axe some noteworthy features in these maps. Small scale structure, present in the fake sky and in the residual maps, is gradually reduced in the residual maps during the first three iterations and is largely absent in the residual map at iteration 11. Between iterations 11 and 13, small-scale structure possessing a pattern different than the small-scale structure pattern at early iterations, begins to grow, and becomes the dominant residual pattern at iteration 25. Structure around the 150" scale remain constant during iteration. The edges of the maps also remain constant, because they are outside of the SHADES field and not measured. In this simulation the iterative algorithm diverges. The simulation described above is repeated. This time a single pointing of the SHADES measurement strategy is used to measure a fake sky populated by a single beam convolved point source with an amplitude of ten. Figure 3.2 depicts the residual maps for iterations 1 through 6 and iteration 25. As in the previous simulation, the residual flux oscillates and grows without bound. The initial map estimate might have an effect on algorithm convergence. To test this, the simulation used to produce Figure 3.1 is repeated, with an initial map 3 The J C M T PSF is approximated as a Gaussian with a full-width half-maxmum equal to 14.5" and an amplitude of one. A null map is a map whose pixel values are zero. 4 Chapter 3. Iterative Reduction of SHADES Data 35 Figure 3.1: Residual maps of a simulation of the generalized iterative algorithm. The measurements are performed on a statistically plausible fake sky using the full SHADES measurement strategy. The measurements are noiseless. High spatial frequency structure is removed from the residuals in the early iterations. In later iterations the unstable spatial frequencies begin to grow. The fake sky map saturates at pixel values of zero and ten. Pixels are three arcseconds square. Chapter 3. Iterative Reduction of SHADES Data 20 40 60 80 100 20 40 60 80 100 2 0 4 0 36 6 0 8 0 1 0 0 I t e r a t i o n 25 -40 -20 0 20 40 Figure 3.2: Residual maps for iterations 1 through 6 and 25 in a simulation of the generalized iterative algorithm applied to measurements taken of a single source. The measurements comprise a single pointing in the SHADES measurement strategy. Chapter 3. Iterative Reduction of SHADES Data 37 estimate equal to the fake sky map. Figure 3.3 shows the fake sky map, fake sky map estimates for iteration 1, iteration 11, iteration 21, and difference maps between the fake sky map estimate and the fake sky map for iteration 1, 11, and 21 for this simulation. The algorithm does not diverge from the perfect initial estimate. The initial map estimate is set to equal 1.1 times the fake sky map and the simulation is repeated. Figure 3.4 shows the residual maps for iteration 1, 15 and 25. By iteration 25 the unbounded mode in the final iterations of the simulation depicted in Figure 3.1 has re-appeared. This mode grows without bound and the algorithm diverges. It is interesting to note that the average pixel value in the residual maps of Figure 3.4 is around zero, whereas in Figure 3.1 the mean pixel value is between —0.9 and —0.4. The only change between the simulations producing Figure 3.1 and Figure 3.4 is the initial map estimate. Thus, to some extent, the initial map estimate affects the mean of the map pixels in this situation. At this stage, there is no noise, and all pixels involved in measurements are being updated. In this situation, assuming algorithm convergence, the edge pixels are well-defined, and no special treatment is required to determine the pixels at the map edges. In later chapters, not all pixels in each measurement will be updated and the measurements will be noisy. Under those circumstances, the edge behaviour is complicated and is addressed in section 4.1 on edge effects for a stable algorithm. The iterative algorithm instability might depend upon the size of the field measured, the measurement strategy, the conditions near the edge of the measured field, the initial map estimate, and the sky being measured. The single pointing simulation has around 10000 pixels and the full SHADES simulation half a million. To investigate the effect of boundary conditions and to compute simulation pixel-pixel correlation matrices, the subsequent simulations are reduced to seven pixels. 3.2.2 Small Field Simulations The algorithm is known to converge when reducing single-difference measurements. A maximally measured seven pixel, single-difference simulation is performed. The fake sky map is s T = [ 0 0 01 0 0 0 Chapter 3. Iterative Reduction of SHADES Data 38 where the s* element represents the flux from the i pixel. The vector of measureth h ments is obtained by, M = Ps where P is the pointing matrix and is given by: /-I 1 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 -1 V o 0 0 0 0 P = 1 (3.7) 0 -1 1 / Thus, M T = [ ° 0 1 -1 0 0 Note that there are six measurements and seven unknown pixel values, in the maximum number of measurements on a one-dimensional, non-periodic, seven pixel patch of sky. This is a manifestation of the non-invertibility of P P, and is due to the T inability of the single-difference measurement scheme to measure a constant sky offset. Figure 3.5 shows the measured fake sky map, the fake sky estimate maps for iterations 45 to 50 for the iterative reduction of the measurements M, and the standard deviation and mean of the map pixel values as a function of iteration. The initial map estimate is the zero map. There are two estimates of the fake sky map in the middle plot of Figure 3.5. The algorithm oscillates between these two estimates with each iteration. The third plot in Figure 3.5 shows the mean and standard deviation of the map pixels constantly oscillating about some value. This is due to the oscillation in the fake sky estimate maps and continues with iteration indefinitely. The algorithm is marginally stable. The simulation is repeated with a single extra, "wrap-around", measurement. Chapter 3. Iterative Reduction of SHADES Data 39 The new pointing matrix is 0 0 0 0 \ 0 0 0 0 0 0 - 1 1 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 -1 1 \ 1 0 0 0 0 0 -1/ / - l l 0 0 - 1 1 (3.9) and the transpose of the measurement vector is M 1 0 0 1-10 0 0 (3.10) Figure 3.6 once again depicts the fake sky map, the fake sky estimate maps and a plot of the mean and standard deviation of the fake sky estimate maps with iteration. This time there is a single estimate of the fake sky map, correct up to a a constant offset. After a small number of osciallations the mean and standard deviation of the estimated map settle to a single value and the algorithm converges. The addition of the extra measurement changes both the simulation boundary conditions and the simulation measurement strategy, and increases algorithm stability. From P in equation 3.9 we can see that the extra measurement constrains the first and last pixel. It is difficult to determine if the effect is due to boundary conditions, or due to a change in measurement strategy since the boundary conditions refer to the way the edges of the map are measured. Besides the pointing matrix, P, the pixel-pixel correlation matrix, P P, is very relevant in the solution to the linear system and T is listed in Table 3.1. It is given here, for the pointing matrix, P, in equation 3.9, 5 and referred to in later sections when comparing correlation matrices for different measurement strategies. The next step is to measure the same fake sky map with double-difference measurements. The first double-difference simulation uses the maximum number of measurements without wrap-around measurements and with a fixed single pixel 5 See section 2.5 Chapter 3. Iterative Reduction of SHADES Data Csd7 P i x 2 -1 0 0 0 0 -1 -1 2 -1 0 0 0 0 0 -1 2 -1 0 0 0 0 0 -1 2 -1 0 0 0 0 0 -1 2 -1 0 0 0 0 0 -1 2 -1 V -1 0 0 0 0 -1 2 7 Meas 40 Table 3.1: Correlation M a t r i x : For a seven pixel, seven measurement, singledifference simulation. chop size. The resulting pointing matrix is P = -.5 1 -.5 0 0 0 0 0 -.5 1 -.5 0 0 0 0 0 -.5 1 -.5 0 0 0 0 0 -.5 1 -.5 0 0 0 0 0 -.5 1 -.5 (3.11) The measurements are M = 0 -.5 1 -.5 0 (3.12) The double-difference measurements do not measure an offset or a constant slope, and thus there are onlyfivemeasurements. Figure 3.7 shows the map estimates for the last six iterations in a simulation with 25 iterations. The initial map estimate is the zero map. This time the algorithm is not marginally stable, as in the seven pixel, six measurment single-difference measurement simulation. Again on each iteration, the pixel values change sign and grow. The average also oscillates and grows, and the standard deviation exhibits something like exponential growth. The simulation is repeated with seven extra measurements added to the pointing matrix given by equation 3.11. The additional measurements include measurements with asymmetrical chop sizes, different chop sizes and "wrap-around" measurements. Figure 3.8 shows the simulation results. The estimate maps exhibit asymmetry, Chapter 3. Iterative Reduction of SHADES Data 41 and oscillate and grow with iteration. The asymmetry is due to the asymmetric measurements amongst the seven newly added measurements to the pointing matrix in equation 3.12. The pixel values grow without bound; however, the growth is slower than the growth depicted in Figure 3.7 for the double-difference, five measurement simulation. Adding two more measurements and repeating the simulation yields the results shown in Figure 3.9. The map is, once again, asymmetric and converged to within a constant of the correct fake sky map. The convergence is accurate to approximately 5% of the maximum pixel value in the fake sky map by iteration 100. The third plot in Figure 3.9 shows the map average converging within the first five iterations to a slightly positive value, and the map standard deviation decaying approximately exponentially. The stability of the iterative algorithm depends upon the measurement strategy. It is not clear what the dependence is and if other factors are contributing. It is also evident that there is a property of the single-difference pointing matrix, P, which is different from the double-difference pointing matrix and different between the 13 and 15 measurement, seven pixel double-difference simulations. A difference is the fewer unmeasured modes in the single-difference P P matrix, resulting T in more equations in the same number of unknowns. It is instructive to look at the pixel-pixel correlation matrices for the different simulations. 6 Both the pixel- pixel correlation matrix in Table 3.1 and the pixel-pixel correlation matrix in Table 3.3 correspond to converging iterative algorithm simulations. C d S many zero off-diagonal terms and dd 7Pix i5Meas agonal entries to off-diagonal entries than exhibits 7Pix 7Meas has generally a larger ratio of di- Cdd7Pixi3Meas- I the direct inversion of n the least-squares method of data reduction, the presence of large off-diagonal terms is indicative of correlated pixels in the data-reduced map, and is undesirable. Convergence depends upon measurement strategy, but the understanding is qualitative and it is not known whether it is the measurements at the edge of the map, as in the 6 T o a i d d e s c r i p t i o n , t h r o u g h o u t t h e t h e s i s p i x e l - p i x e l c o r r e l a t i o n m a t r i x refers t o P P. T This is a c c u r a t e , u p t o a f a c t o r , w h e n t h e t i m e l i n e n o i s e c o v a r i a n c e m a t r i x is d i a g o n a l , w i t h i d e n t i c a l entries. Chapter 3. Iterative Reduction of SHADES Data I 2.75 -1.5 0 -0.5 0.25 0 -1.5 2.75 -1 0.25 -0.5 -0.5 0.5 0.25 -0.5 -0.5 -1 0.25 -0.5 2.75 -1 -0.75 -1 2.75 -1 0 -1 2.75 -1 -0.5 0.25 -1 7 P i x 13 M e a s 42 2.5 0.25 -0.5 0.25 -1 0 -0.5 -0.5 0.25 0.5 V -1 -1 \ -0.5 -0.5 -0.75 -1 3.25 J Table 3.2: Correlation M a t r i x : Seven pixel, thirteen measurement, doubledifference simulation. Cdd 7 P i x 15Meas — 3 -1.5 0 -0.25 0.25 0 -1.5 -1.5 2.75 -1 0.25 -0.5 -0.5 0.5 0 -1 3 -1 0.25 -0.25 -1 -1 2.75 -1 0.25 -1 0.25 -1 2.75 -1 -0.75 0.25 -1 3 -1.5 -1 -0.75 -1.5 5.25 -0.25 0.25 0.25 0 -0.5 -0.5 -0.25 V -1.5 Table 3.3: Correlation 0.5 Matrix: difference simulation. -1 Seven pixel, fifteen measurement, double- 43 Chapter 3. Iterative Reduction of SHADES Data single-difference simulation in this section, or the initial map estimate, or the bulk of measurements at the centre of the field which determines the iterative algorithm convergence properties. Specifically, the following questions remain to be answered: 1. How does convergence depend on initial map estimate? 7 2. How does convergence depend on the unmeasured modes? 3. How does convergence depend on the boundary conditions? 8 4. Do the small field simulations sufficiently model conditions in the actual data reduction to draw valid conclusions? 5. How does convergence depend on noise? Note that the SHADES measurement strategy isfixedand cannot be manipulated to achieve convergence. To determine the answer to the above questions by simulation requires the fixing of one parameter, say the initial map estimate, and then varying the others through all parameter space, and then repeating for different values of the parameter under investigation. This is a difficult and exhaustive procedure, which the discovery of a reformulation of the iterative deconvolution algorithm made largely unnecessary. It is the subject of Section 3.4 onwards. 3.2.3 Minimum Comparable Size The pointing matrix in the 15 measurement double-difference simulation of the last section uses measurements with all physically possible chop throws, "wrap around" measurements, and measurements with asymmetrically sized chop throws of a onedimensional field. Despite this, the pixel-pixel correlation matrix in Table 3.3 is not very diagonal, and does not constitute an acceptable measurement strategy. Com9 pare this to measurement strategies such as that used by WMAP, which introduces very small correlations between the pixels when iteratively deconvolving, and whose 7 8 "How" here means "Does? And if so, in what way?" "Boundary conditions" refers to the tapering in the number of measurements of a region as the observation nears the edge. U p to 20% correlation between pixels. 9 Chapter 3. Iterative Reduction of SHADES Data 44 pixel-pixel correlation matrix is highly diagonal. WMAP does not perform every possible single difference measurement, yet achieves a relatively diagonal pixel-pixel covariance matrix. How is this accomplished? The number of unique measurements that can be made of a one-dimensional field is floor N= m (^) £ (3.13) i=l where floor (x) returns the truncated integer representation of x. For a seven pixel map measured using double-differences there are nine unique measurements, not involving wrap-around or asymmetric measurements. Thus, the number of possible measurements of a partial-sky field is expected to be O (-Np\). It is reasonable to x expect that observations of large fields can achieve negligible correlations without exhausting all measurement possibilities. The covariance matrix of the SHADES simulation is too big to be easily constructed and compared directly to the small simulation covariance matrices. Instead, a figure of merit is computed, which is the ratio of the central pixel's diagonal entry in the pixel-pixel correlation matrix to the sum of the absolute entries in the row. [2] This figure of merit has a maximal value of one and is applicable where the measurement strategy is uniform in the centre of the observed field. Thefigureof merit of the SHADES measurement strategy for a single pointing is 0.376, for the singledifference, seven pixel, seven measurement simulation 0.5, for the double-difference, seven pixel, thirteen measurement simulation 0.417 and for the double-difference, seven pixel,fifteenmeasurement simulation thefigureof merit is 0.423. Calculating thefigureof merit for the entire SHADES strategy, though requiring much less memory than the full SHADES correlation matrix, requires the full pointing matrix (minus the zero entries) which is still prohibitive. The iterative deconvolution algorithm applied to the full SHADES data set in simulation is slightly unstable and if thefigureof merit is a good indicator of algorithm stability then we expect a figure of merit for the full SHADES simulation between 0.417 and 0.423. Chapter 3. Iterative Reduction of SHADES Data 3.3 45 Iterative Inversion As Newton-Raphson Root Finding A method of determining the roots of a non-linear function of one variable is to make a guess as to the location of the root, compute the tangent of the function at this location, and then compute where the root would be located if the function were the tangent. The estimate of the root is the value of the function evaluated at the new estimated root location. The process is repeated until convergence. The rate of convergence in this algorithm depends upon the initial guess and the degree of non-linearity of the function. Divergence typically occurs when the slope of the tangent line to the curve at the estimated root location is zero. [13] Consider the equation Px = M (3.14) where P is the pointing matrix, x is the map of pixels, and M is the corresponding vector of measurements. Form an error vector E, such that E(x) = M-Px (3.15) The best estimate map will be the map corresponding to the roots of E. That is, the position in N dimensional space corresponding to an error vector of minimal p magnitude. Consider a small deviation, Sx, from the current map estimate, x, then NpiK BE Ei (x + Sx) = Ei (x) + V ——+ j=l where Ei is the 2 t h d X O {Sx ) 2 (3.16) j element of E. Let x* = x + Sx (3.17) We desire x' to be the root of E (x). Neglect the higher order terms in 3.16, and set E (a?) = 0. Solve for 8x and obtain 6x = -J~ E = -P l (3.18) where 13 dxj (3.19) Chapter 3. Iterative Reduction of SHADES Data 46 and P is the pointing matrix. The Jacobian, J, has dimensions N x iVpi and is m x not invertible. The least-squares optimal update is Sx = {P P)~ T 1 [P - Ps)] (M T (3.20) A couple of notes. First, P P is not invertible, is large, and its inversion is the T original problem, ie. 8x = x - (P P)~ T (3.21) PM 1 T Thus, the new map estimate is + SX= {P P)~ T £ n e w = 2?old 1 (3.22) PM T which is the least squares optimal solution to equation 3.14, is the topic of section 2.5.3, and is whose computation was hoped to be carried out by iteratively applying equation 3.1. Secondly, the error function is linear. That is, the higher order terms in equation 3.16 are zero, and convergence is expected in a single step, independent of the initial map estimate. The approach is abandoned, as it is simply a restatement of the original problem to which iterative deconvolution is addressed. 3.4 Formulation of Iterative Inversion as a Difference Equation The following formulation is motivated by a desire to mathematically describe the map update in a single equation. 3.4.1 Formulation Let the error between the data and the measurement of the map estimate be E = M - Px n (3.23) n where M is a vector of measurements, P is the pointing matrix and x is the n th n map estimate. Note that the error on the j th measurement is El = Mj - g- • xn (3.24) Chapter 3. Iterative Reduction of SHADES Data where g~j is a vector of gains corresponding to the j 47 measurement. Now the update t h for the 2 pixel, given by equation 3.1, is t h n-1 (3.25) +x 9i • 9i The numerator in the fraction of the above equation is an inner product, (3.26) ^ g i E ^ g - i - E and hence, 9i-E F + (3.27) n-l x g% • 9i Thus, each pixel update is an inner product between a gain vector for the 2 pixel t h and a vector of errors, scaled by the inverse square magnitude of the gain vector plus the previous pixel value. In a sense, to within a scale factor, the inner product between the gain vector and the error vector is the act of measuring the error contributed to the 2 pixel by the difference between the actual measurements and t h the measurements of the current map estimate. The last equation can be extended to give the update of the entire map. \ 1 90-90 / 91-91 go ^ 91 E -1 n \ 9N -9N p 9i Noting that \9N P P + ^n-l (3.28) 9N J P ) is equal to the tranpose of the pointing matrix P, and subJ stituting equation 3.4 for E -i, equation 3.9 can be simplified with the following n Chapter 3. Iterative Reduction of SHADES Data 48 algebra. go-go gi-gi P • ( M - P £ „ _ i ) + f„_i(3.29) T — \ §Np-gN J p — (I - DP P) = A • f„_i + (3.30) • x -i + £ > P • M T T n (3.31) 5 l go-go 9i-Si where D = 9Wp -9Af A = 7- / p DP P T and P = DP T •M . Equation 3.31 is a difference equation giving the n estimate of the map. th 3.4.2 Solution The solution is standard and is achieved by mathematical induction. From equation 3.31 the first three map estimates are (3.32) x\ = A • xo + B x = A • xi xz = A • x + B = A (A • x + AB + B) + B 2 + B = A • x + AB + B 2 0 2 2 = 0 (3.33) (3.34) A • x + A B + AB + B 3 2 0 By induction the solution is n-l x = A • xo + n A n r=0 T B (3.35) Chapter 3. Iterative Reduction of SHADES Data Equation 3.35 gives the estimate of the map produced by the n th 49 iteration of the generalized Wright algorithm. There are three easy tests of equation 3.35. They are: 1. Solution behaviour with a correct initial map estimate. 2. Solution behaviour with a diagonal pixel-pixel correlation matrix. 3. Solution behaviour as the number of iterations tends to infinity. The checks are presented in the following three sections. 3.4.3 Solution Check #1 If XQ is correct, then P-x 0 =M (3.36) The n map estimate should equal the initial map estimate if this is correct. The t h n th map estimate becomes: n-1 (3.37) and DP P T = (3.38) I-A then n-1 (3.39) n-1 (3.40) A • x + (J - A ) • x (3.41) x (3.42) n 0 0 as expected. 0 Chapter 3. Iterative Reduction of SHADES Data 3.4.4 50 Solution Check #2 In section 3.3.2, it is noted that the stability of the generalized iterative data reduction algorithm depends, at least in some way, upon the measurement strategy. A perfect measurement strategy results in a diagonal pixel-pixel correlation matrix. Furthermore, the diagonal elements of D in the generalized Wright algorithm scheme are equal to the inverse of the diagonal elements of the pixel-pixel correlation matrix. Thus, with a perfect measurement strategy, DP P is equal to the identity matrix, T and A equals zero. With A set to zero in equation 3.35, a single term survives: x = DP • M (3.43) T n Noting that DP P T = I: DP T =P (3.44) + Where P is the pseudo-inverse of P. Thus, + 10 x =P •M (3.45) + n When the sky is measured perfectly, the iterative algorithm gives the solution of P • x = M on any iteration up to a linear combination of the singular eigenvectors of P. 3.4.5 Solution Check #3 Consider the limit as n -> oo. It is desired that lim^ooXn = P + (3.46) •M Starting with the summation term in equation 3.35, note that n n+1 ^A -^A r r=0 r = /-A" +1 (3.47) r=l n A [I - A] = I - A r n+1 (3.48) r=0 n-1 and Y,A r r=0 °See, f o r i n s t a n c e , [21]. = [I- A ][I-A] n + (3.49) Chapter 3. Iterative Reduction of SHADES Data 51 Equation 3.35 becomes x = A • x + [I - A ] [I - A} n n n (3.50) + 0 Diagonalizing the N by N matrix A obtain p p A = vXv' (3.51) 1 and A = vXv^vXv' 2 = vX v~ 1 2 (3.52) 1 By induction A n = vX v~ n (3.53) l where v is a matrix composed of the eigenvectors of A, and A is a diagonal matrix whose diagonal entries are the eigenvalues of A. Note that the diagonalization is impossible if A does not have iVpi linearly independent eigenvectors. A is a x nearly symmetric, band-diagonal, invertible matrix, which has zero entries along the diagonal. A rigorous analysis can be performed using the Jordan Form. The Jordan form analysis yields the same results with a slight modification of the algorithm convergence when the absolute value of the eigenvalues approachs one. x = vX^' 1 n • x + [I- vX v~ ] [I - A} B n 0 l 11 Now (3.54) + and l i m ^ o o f n = vX°°v- -x + 1 [I- 0 vX^v' ] 1 [I - A} + B (3.55) If the absolute value of the eigenvalues of A are less than one, then we see that there is no dependence on the intial map estimate So and linw^ lim^oofn = [I-A] B = [DP P] = P = P -M (3.56) + T + DP + •M T (DP ) T + DP T •M + (3.57) (3.58) (3.59) which is the desired solution up to a linear combination of singular eigenvectors of DP P. If any of the absolute values of the eigenvalues of A are greater than one T then x will grow without bound with increasing n, and is unstable. n 1 1 See Appendix A. Chapter 3. Iterative Reduction of SHADES Data 52 -1.73 -1.33 -0.346 0.500 0.909 1.00 1.00 Table 3.4: Eigenvalues of the characteristic matrix for the seven pixel, five doubledifference measurement simulation. -0.522 0.298 -0.308 0.295 -0.308 0.298 -0.522 -0.606 0.303 -0.202 0.000 0.202 -0.303 0.606 0.617 -0.177 -0.140 0.369 -0.140 -0.177 0.617 -0.632 0.000 0.316 0.000 -0.316 0.000 0.632 0.627 0.215 -0.141 -0.285 -0.141 0.215 0.627 0.378 0.378 0.378 0.378 0.378 0.378 0.378 0.378 0.378 0.378 0.378 0.378 0.378 0.378 Table 3.5: Eigenvectors of the characteristic matrix in the seven pixel, five doubledifference measurement simulation. Each row represents a separate eigenvector. The first row is the eigenvector corresponding to the left-most eigenvalue in Table 3.4. For clarity, the eigenvector components have been rounded to three significantfigures.The eigenvectors are not linearly independent. 3.4.6 Difference E q u a t i o n S o l u t i o n S i m u l a t i o n The generalized Wright algorithm is simulated in a seven pixel simulation with five double-difference measurements. The 25 iteration map estimate is compared with th the map estimate computed using equation 3.35. The maps are identical up to machine precision. The eigenvalues and eigenvectors of the characteristic matrix A are given in Tables 3.4 and 3.5, respectively. There are two absolute eigenvalues greater than one, and two eigenvalues that are identically one. The eigenvectors corresponding to the eigenvalues of one are the same and correspond to the offset of the map to within floating point precision. The eigenvectors are not linearly independent and the characteristic matrix is not diagonalizable. Table 3.6 and 3.7 list the eigenvalues and eigenvectors for a seven pixel, 15 double-difference measurement simulation. All of the absolute eigenvalues are less Chapter 3. Iterative Reduction of SHADES Data 53 1.00 -0.987 -0.685 0.447 0.050 -0.060 0.235 Table 3.6: Eigenvalues of the characteristic matrix for the seven pixel, 15 doubledifference measurement simulation. 0.378 0.378 0.378 0.378 0.378 0.378 0.491 -0.573 0.401 -0.252 0.333 0.080 -0.298 0.174 -0.079 -0.195 0.553 -0.444 0.563 -0.325 0.509 -0.367 -0.415 -0.174 -0.159 0.606 0.105 0.438 -0.500 -0.293 -0.153 -0.454 0.403 0.064 -0.578 -0.003 0.159 0.354 -0.487 -0.550 0.006 0.378 0.493 0.428 -0.426 -0.100 0.329 0.571 -0.014 Table 3.7: Eigenvectors of the characteristic matrix for the seven pixel, 15 doubledifference measurement simulation. Each row is a separate eigenvector; the eigenvector in the top row corresponds to the left most eigenvalue in Table 3.6. Chapter 3. Iterative Reduction of SHADES Data 54 than or equal to one, and the simulation is stable (See Figure 3.8). The solution is correct up to a map offset which is the eigenvector corresponding to the eigenvalue of one. Examining the last plot in Figure 3.8, we see that the average of the map exhibits some dynamics at early iterations and settles to a positive constant value. Keeping this in mind we look at equation 3.35 in the context of the seven pixel, 15 double-difference measurement simulation. In this simulation, the initial map estimate is the zero map and equation 3.35 becomes n-1 s n = J2 ( - °) ATB 3 6 r=0 ie. x = P M - AP M + n n + (3.61) The first term in equation 3.61 is the desired solution and the second term we expect to die off with iteration. Unfortunately the eigenvectors in Table 3.7 are not orthogonal (as expected, since A is not symmetric) and are not linearly independent, and thus A is not diagonalizable. The behaviour of the map pixel average can be computed using the Jordan form, but it is preferable to analyze a simpler case. Nonetheless, qualitatively, with iteration the contributions to the map and to the mean of the map due to the eigenvectors with absolute eigenvalues less than one decrease, and we see the map mean fluctuate at early times and then stabilize. The reader is again referred to Appendix A for a discussion of the Jordan form. 3.4.7 Stability As previously mentioned, the absolute eigenvalues of the A matrix in equation 3.35 are dominant in determining the stability of the iterative algorithm. What type of measurement strategy results in absolute eigenvalues greater than one? The answer is measurement strategies that result in a wider band in the band diagonal P P T matrix. The A matrix, in equation 3.35 is zero along its diagonal. Heuristically, the larger off-diagonal bands in A force the eigenvectors of A to contain greater structure. The physical characteristics of measurement strategies resulting in narrow band-diagonal A matrices is addressed in the section on measurement strategy. In the next chapter the rate of convergence of simulations with the same measurement strategy will vary depending upon the fake sky that is measured. This is Chapter 3. Iterative Reduction of SHADES Data 55 expected; consider the difference between the x and x -i maps: n n (3.62) ^2 A DP M = r - ^2 T r=0 Ax„ (3.63) A DP M r T r=0 = (3.64) A ~ DP M n l T Let A be diagonalizable so that A = vXv (3.65) 1 then Ax n (3.66) v\ - v- DP M = n 1 1 T Convergence will vary depending on the degree of orthogonality between DP M T and the inverted eigenvectors corresponding to large absolute eigenvalues of A. Thus, one expects convergence to depend upon the measurements to a small extent. It might be possible that absolute eigenvalues greater than one correspond to attempting to invert a non-invertible matrix without constraining the unmeasured modes of P P. Fixing a pixel, in the context of the iterative algorithm, means to T never update a pixel involved in measurements. In this way the pixel is constant with iteration and stays fixed. A six pixel fake sky map is measured with the following pointing matrix. / -.5 P = \ 1 -.5 0 0 -.5 1 0 0 -.5 0 0 0 0 \ 0 -.5 0 0 -.5 0 (3.67) -.5 -.5 1 The pixel fixing can be performed by modifing the P matrix such that the first T and last rows are zero, i.e. 0 0 0 \ 1 -.5 0 0 -.5 1 -.5 0 0 -.5 1 -.5 0 0 -.5 1 0 0 / 0 P = 1 V o 0 / (3.68) Chapter 3. Iterative Reduction of SHADES Data The eigenvalues of the characteristic matrix A — I - DP P T 56 are listed in Table 3.8. There is an absolute eigenvalue greater than one. -1.3 0.64 1.0 -0.31 0.97 1.0 Table 3.8: Eigenvalues of the characteristic matrix for the six pixel, four doubledifference measurement simulation with two fixed pixels. There is an absolute eigenvalue greater than one. Perhaps setting, otherwise unmeasured pixels to a value will make the algorithm stable. In this case the pointing matrix, corresponding to two absolute measurements and two double-difference measurements, is f\ 0 0 -.5 P = 0 \0 P T 0 0 0 0 0 0 N 1 -.5 0 0 -.5 1 0 0 -.5 0 0 (3.69) 1 really equals the transpose of the pointing matrix and the eigenvalues of the characteristic A matrix are listed in Table 3.9. There are absolute eigenvalues greater 0 1 -1.8 -0.2 1 0 Table 3.9: Eigenvalues of the characteristic matrix in the simulation with six pixels, two double-difference measurements and two absolute measurements. There is an absolute eigenvalue greater than one. than one. Note that thefixedpixels in this case are not measured by the doubledifference measurements. Next add two more double-difference measurements linking the absolute mea- Chapter 3. Iterative Reduction of SHADES Data 57 surements. The resulting pointing matrix is: / 1 0 0 0 0 0 -.5 1 -.5 0 0 0 0 -.5 1 -.5 0 0 0 0 -.5 1 -.5 0 0 0 0 -.5 1 -.5 \ 0 0 0 0 0 (3.70) 1 ) Again, P really equals the transpose of the pointing matrix and the eigenvalues of T the characteristic A matrix are listed in Table 3.10. There is an absolute eigenvalue -.62 .72 .31 .03 0.98 -1.4 Table 3.10: Eigenvalues of the characteristic matrix for the six pixel, four doubledifference measurements and two absolute measurements simulation. There is an absolute eigenvalue greater than one. greater than one. As a prelude to the discussion of the central update algorithm in the section on algorithm stabilization and in chapter four, consider the following modified P T matrix: / 0 0 0 0 \ 10 P = l 0 0 0 10 0 (3.71) 0 0 10 0 0 0 1 \0 0 0 0/ The resulting eigenvalues are listed in Table 3.11. 3.4.8 Noise Properties If the prior probability distribrution is uniformly distributed and non-zero over a reasonable range of pixel values, and the measurement noise is Gaussian , the width Chapter 3. Iterative Reduction of SHADES Data -1.0 58 0.81 -0.81 0.31 -0.31 1.0 Table 3.11: Eigenvalues of the characteristic matrix for the six pixel, four doubledifference measurement and two absolute measurement. There are no absolute eigenvalues greater than one. of the likelihood function is given by the pixel-pixel noise covariance matrix. The pixel-pixel noise covariance matrix, C, is computed as follows: (3.72) C = {x xl) - (x ) {x ) T n n n To compute this, first the expected map estimate is determined using equation 3.33. Let the absolute eigenvalues of A be less than one, then the n th estimate is, expected map 12 <£„> = (x ) ((DP P) DP T + -(M + n)) T = (P {DP ) DP + n (3.73) 0 T + -(M + n)) T (3.74) 0 {£„) = (P+• M + P • n) (3.75) (x ) (3.76) + 0 = P+-M n 0 Substituting equation 3.75 for the second term in equation 3.72, C becomes C = ( P ( M + n)(M + n ) ( P ) ) - P M M^P + T 0 + + T +T 0 0 (3.77) since the expectation of n is zero, and the last term cancels the first term in the product, the pixel-pixel noise covariance matrix is C = {P rtn P ) + T = P NP +T + +T C = P NP + (3.79) +T C = (P N- P) T (3.78) 1 (3.80) + If P P is invertible, then this is equivalent to, T C = (P N- P)~ T 2 1 1 T h i s is c o r r e c t u p t o t h e n o n - i n v e r t e d s u b s p a c e o f t h e m a t r i c e s b e i n g p s e u d o - i n v e r t e d . (3.81) Chapter 3. Iterative Reduction of SHADES Data 59 where N is the timestream noise covariance matrix. Note that C is equivalent to the pixel-pixel noise covariance matrix obtained from the generalized least-squares formulation. [9] 3.5 Algorithm Stabilization The measurement strategy of SHADES is fixed, and if no modifications of the iterative algorithm (equation 3.35) are made, the algorithm does not converge. The goal is to modify the algorithm in a fashion that makes the algorithm stable while not appreciably affecting the quality of the solution. The following algorithm modifications are attempted: 1. Unstable eigenvector identification and removal. 2. Modifying the pixel size at the boundaries of the map. 3. Decoupling of measurements and pixel updates. 3.5.1 U n s t a b l e Eigenvector R e m o v a l The idea is to determine the eigenvectors and associated eigenvalues of A in equation 3.35 which have absolute eigenvalues greater than one, and subtract them, appropriately scaled, from the final map. The trick is to find the eigenvalues and eigenvectors without forming the matrix. Even if the matrix is formed, computation of the eigenvalues and eigenvectors is an 0(n ) operation.[13] A potential probe of 3 unstable eigenvectors is obtained by setting the measurement vector, M, in equation 3.35 to zero. Then, only the first term in equation 3.35 is computed. After many iterations, x =A • f = ^T n n v 0 v l •o (3-82) f i where Vi is a matrix whose columns are the eigenvectors with absolute eigenvalues greater than one. If SQ is guessed to be an eigenvector, then the n th estimate map is x = A£a?o n (3.83) Chapter 3. Iterative Reduction of SHADES Data 60 where A is the eigenvalue associated with an unstable eigenvector of A, namely e the guessed eigenvector xo- ls The ratio of the output maps from two successive iterations gives the associated eigenvalue. Assuming knowledge of the unstable eigenvectors and eigenvalues, the leastsquares optimal map can be recovered. Let u = [I - A]' DP M 1 (3.84) T Equation 3.35, with a zero-map initial guess becomes x = (l- A ) u n+l n (3.85) Consider x - x -i = (-A +A) u n+1 n n n (3.86) Diagonalize A in the limit of large n: A = v A2 v n+1 +1 (3.87) T e where v is a matrix whose columns are populated with the unbounded eigenvectors of A, A is a diagonal matrix populated with eigenvalues and the subscript e stands for explode. That is, all of the eigenvectors, corresponding to absolute eigenvalues less than one, are negligible in the decomposition. Let 14 c = v -u (3.88) i i be the the inner-product of u with the i e th x - unstable eigenvector. Then, = £ 4 (A* ) * (1 - XI) n n (3.89) i and the least-squares optimal map is n x = ^2 A DP M r T n r=0 +^ 4K i +lc (3-90) i Thus, the task is to determine the unstable eigenvectors, their associated eigenvalues and the inner-products of u with the unstable eigenvectors. If these are 1 3 Guessing an eigenvector is somewhat fortunate. 1 4 Again, all of this assumes linearly independent eigenvectors. The formulation can be generalized using the Jordan form. \ Chapter 3. Iterative Reduction of SHADES Data -1.73 61 -1.33 -0.346 0.500 0.909 1.00 1.00 Table 3.12: Eigenvalues of the characteristic matrix for the seven pixel, five double-difference measurement simulation. known, the least-squares optimal map is given by equation 3.90. The characteristic matrix, A, whose eigenvectors we wish to determine is iVp ixel x N j \. p xe Determining the eigenvectors is typically an O [n ] computation, with the number of operations 3 going to O [n ] per eigenvalue for a Hessenberg matrix. Thus, if there were 2 15 16 a way to determine which eigenvalues one would determine prior to applying the algorithm, it might be possible to determine the eigenvectors of A. [13] Future work might investigate non-standard techniques with sparse matrices; however, with the convergence of the central update algorithm given below, this method is abandoned. 3.5.2 Boundary Pixel Modification Changing the size of pixels at the edge of the map as a means of stabilizing the algorithm is investigated. Consider a seven pixel, one dimensional simulation of the iterative data reduction algorithm with the edge pixels enlarged to twice the size of the central pixels. The normal pointing matrix is / -.5 V 1 -.5 0 0 0 0 0 -.5 1 -.5 0 0 0 0 0 -.5 1 -.5 0 0 0 0 0 -.5 1 -.5 0 0 0 0 0 -.5 1 (3.91) -.5 / The associated eigenvalues of the characteristic matrix are given in Table 3.10, and reproduced here for convenience. With doubly enlarged edge pixels, the pointing 1 5 1 6 I f a l l e i g e n v a l u e s are d e s i r e d , t h i s is a g a i n , a n O [ra ] 3 computation. A H e s s e n b e r g m a t r i x c a n b e f o r m e d o u t o f a n a r b i t r a r y m a t r i x t h r o u g h a series o f H o u s e h o l d e r transformations. Chapter 3. Iterative Reduction of SHADES Data -1.57 62 -0.797 0.406 1.00 0.953 Table 3.13: Characteristic eigenvalues for the edge-pixel enlarged pointing matrix. matrix becomes .5 -.5 V 0 0 0 \ -.5 1 -.5 0 0 0 -.5 1 -.5 0 0 0 -.5 1 -.5 .5 / The eigenvalues of the characteristic matrix A are given in Table 3.10. By com0 0 0 -.5 paring the eigenvalues of the original characteristic matrix with the eigenvalues of the edge-pixel size increased characteristic matrix, it is apparent that the absolute value of the largest absolute eigenvalue in the two cases has decreased between the original characteristic matrix and the one with larger edge pixels, causing the iterative algorithm to become more stable. It is possible that algorithm stability might be accomplished by increasing the size of the edge pixels. In a sense, this reduces the number of unknowns while maintaining the same number of equations. Besides stability, enlarging the edge pixels might yield an improved estimate of the map by increasing the number of measurements involving an edge pixel. This is limited, however, since the resolution at the edges deteriorates as edge pixel size increases. That said, it is possible that a modest edge pixel enlargement would stabilize the full SHADES simulation of section 3.2.1. In light of the stabilizing method of the next section this method is not pursued further. 3.5.3 Measurement A n d P i x e l Update Decoupling A stable modification to Wright's algorithm was found by Colin BorysPersonalComm. The modification involves the D and P matrices in equation 3.35 All ocT curences of —0.5 in P are changed to zero, and any entries in D where, T = oo, are set to zero. This modification will henceforth be called the "central update algorithm". The pointing matrix is left intact. Thus, the pointing matrix and the pixel Chapter 3. Iterative Reduction of SHADES Data 63 -0.866 1.00 0.866 0.000 -0.500 0.500 1.00 Table 3.14: Characteristic matrix eigenvalues for the central update algorithm in the seven pixel, five double-difference measurement simulation. All of the absolute eigenvalues are less than or equal to one. update matrix are decoupled. The result is a band diagonal A with a thinner band and smaller absolute eigenvalues. The following eigenvalue computation of A for a seven pixel double-difference pointing matrix with modified transpose demonstrates the effect. This method works because the n 3.55 does not depend on the D or the P T th map estimate given by equation matrices. Thus, when considering only stability, they can be modified in any fashion, so long as the absolute eigenvalues of the resulting A matrix are less than one. The central update algorithm does not modify the pixel values at the edge of a map. An edge pixel that is not involved in the central part of a double-difference measurement is not updated. Thus, the initial map estimate at the edges is fixed, and any sky flux at the edges manifests itself as erroneous pixel values near the edge of the map. The effect decreases with distance from the edge and increases with the amount of flux at the edges of the field. In the next chapter performance of the central update algorithm when reducing SHADES data is characterized by simulation. Chapter 3. Iterative Reduction of SHADES Data 64 Figure 3.3: Top: Fake sky map. Middle: Fake sky map estimates for iteration 1, 11 and 21. Bottom row: Residual maps for iteration 1, 11 and 21. The colour bar applies to the residual maps only. Chapter 3. Iterative Reduction of SHADES Data Iteration Iteration 1 I l 1 Iteration 15 i \I 65 25 i \ ... 200 400 400 600 -i -0.1 -0.05 0 200 600 1 0.05 400 600 1 0.1 Figure 3.4: Residual maps for iteration 1, 15 and 25. Note the reduction of highfrequency spatial structure between iteration 1 and 15 and the reappearance of high-frequency structure in iteration 25. The algorithm slowly diverges. Chapter 3. Iterative Reduction of SHADES Data Figure 3.5: Single-Difference: 7 Pixel Simulation, 6 Measurements 66 Chapter 3. Iterative Reduction of SHADES Data Figure 3.6: Single-Difference: 7 Pixel Simulation, 7 Measurements 67 Chapter 3. Iterative Reduction of SHADES Data Figure 3.7: Double-Difference: 7 Pixel Simulation, 5 Measurements 68 Chapter 3. Iterative Reduction of SHADES Data Figure 3.8: Double-Difference: 7 Pixel Simulation, 13 Measurements 69 Chapter 3. Iterative Reduction of SHADES Data 70 Fake Sky Map 1.2 1.0 - 1 • • E- .... -z 0.8 r 0.8 - z~ _ 0.4 0.2 0.0 -0.2 _ + - + + + 0 1 2 + + + 4 3 6 -j - Map: Iteration 20 To Iterotion 25 1.5 - 1 1 1 1 1 1 1 1 1 1.0 =- 0.5 — 0.0 S -0.5 -1.0 + Iter. 2 0 - Iter 21 i 0 2 . Iter. 22 <>!ter. 2 3 .A. Map: Iteration 9 5 To iterotion 24 lter 4 3 PiK*i 5 Iter. 2 5 z 6 100 1.0 -E 0.5 0.0 • A— 9t — —-B-— • • 0.5 1.0 + 0 Iter. 9 5 1 • Iter 96 . Iter. 97 2 3 0 Iter. 98 .A 4 l t e r 99 5 Iter. 100 : 6 Map Mean and Mop Standard Deviotion vs. Iteration Number + L * Map S t a n d a r d Deviation Map Average + 0 . 2 0 0 0 0 0 H^*M^-H-H*H*H-H^i-H-H-H-H4i^^ H 1 1 H I H I H1111111111 B Iteration Number Figure 3.9: Double-Difference: 7 Pixel Simulation, 15 Measurements Chapter 4. Central Update Algorithm 71 Chapter 4 Central Update Algorithm The central update algorithm was introduced in section 3.5.3. In this chapter the edge effects, the algorithm behaviour when reducing measurements with additive Gaussian white noise, and the algorithm behaviour when reducing measurements with more realistic noise, are investigated by simulation with the SHADES measurement strategy. The application of the algorithm to real SHADES data, and subsequent source detection, are described in section 5.3. 4.1 Edge Effects To investigate the edge effects, a fake sky with a single, centred source is measured with the full SHADES measurement strategy and reduced using the central update algorithm. No noise has been added to the measurements. Figure 4.1 is a plot of the standard deviation of the residual map pixels as a function of iteration. Only pixels 150" from the map edge are included in the calculation. The standard deviation of the residual map pixels is 0.0008 mjy at iteration 101. At iteration 1001 the standard deviation of the residual map pixels is less than 0.0001 mjy Because there is no noise in this simulation, it is expected that the standard deviation of the residual pixels will asymptotically approach zero. Figure 4.2 shows residual maps for various, arbitrarily chosen, iterations. The central update algorithm converges and the edge effect, if there is one, is negligible by iteration 75. The simulation is repeated with a realistic fake sky map whose periphery is populated with sources drawn from a realistic source count model. The periphery is 150" wide (50 pixels). This is the same statistical source number count model used for various simulations of the generalized algorithm in Chapter 3. The purpose Chapter 4. Central Update Algorithm 72 Residual Map Standard Deviation 0.005 0.004 0.003 •5 0.002 0.001 0.000 400 600 Iteration Number 1000 Figure 4.1: Standard deviation of the residual map pixels as a function of iteration for a central update algorithm data reduction simulation. The simulated data consist of measurements of a single, 10 mjy, centred source, using the full SHADES measurement strategy. The computation neglects pixels within 150" from the map edge. The standard deviation of the residual map pixels is 0.0008 mjy at iteration 101. At iteration 1001 the standard deviation of the residual map pixels is less than 0.0001 mjy. Chapter 4. Central Update Algorithm Residual 600 -I 400 - J 200 • Residual 1 73 Residual 2 5 I B •11 •1 1 • 1 1 •1! M l 200 400 600 200 -0.4 Residual 11 600 J I L 0.2 0.4 0 Residual •ft 600 — -0.2 400 an i : n 1 200 — 200 400 600 T 200 400 400 Residual 21 * • 400 200 600 600 75 •I •1 i 200 i r 400 600 _i_ -0.01 -0.005 "1 1 0 0.005 0.01 Figure 4.2: Residual maps estimated on various iterations of a simulation of the central update algorithm using the full SHADES measurement strategy. The fake sky map is a single source, centred in the map, with a peak flux of ten. The pixel values are in mJy and the contrast has been stretched on the lower three panels. The first iteration map is a measurement map. Note the convergence and the absence of a visible effect caused by the edges of the map. The maximum absolute pixel value of the residual map on the 75 iteration is less than a tenth of a percent of the source th amplitude. Chapter 4. Central Update Algorithm 74 of this simulation is to determine the edge-effect produced by the central update algorithm in a realistic SHADES simulation. Inside the periphery the central portion of the fake sky map is set to the null map. Figure 4.3 is a plot of the standard deviation of the residual map pixels as a function of iteration. Only pixels inside the residual map periphery are used in the computation of the standard deviation of the residual map pixels. At iteration 401, the residual map pixel standard deviation is within 0.02 mJy of the converged value of 0.07mJy. The steps in the curve between iterations 401 and 801 are due to limited precision when reading and writing data to file. Note the maximum residual map pixel standard deviation of 0.185 mJy 1 near iteration 60. At iteration one, the initial map estimate is the null map, and all pixels used in the computation of the standard deviation are equal to the fake sky map. Thus the standard deviation of the residual map pixels is zero, and rises with iteration as pixels inside the periphery are modified due to flux in the chop positions of the measurements centred upon pixels at the edge of the periphery. Note that convergence occurs much faster with measurements of the single centred source than with measurements of the statistically plausible edge-only fake sky, even though the measurement strategy, and thus the eigenvectors and eigenvalues of the characteristic matrix, are identical. This is due to the coupling described in section 3.4.7 between the measurements and the characteristic matrix. Figure 4.4 shows the fake sky map and residual maps for several iterations. With iteration, edge-effect flux propagates into the centre of the residual map when reducing the SHADES data with the central update algorithm. Note that the residual map mean grows increasingly negative with iteration and reduces the residual map pixel standard deviation in the higher iterations with respect to the maximum of the residual pixel standard deviation near iteration 60. At iteration 101, near the peak of the residual map standard deviation, thefluxresidual due to the edge effect ranges from 0.25 mJy 360" from the map edge, to 0.5 mJy, 80" from the map edge. Figure 4.5 is the 801 * residual map. At iteration 801 the flux residual is s less than 0.15 mJy 150" from the map edge. The 50 pixels, 150", map edge is approximately the size of the outermost part of the map which is not measured in the 1 Author oversight. Chapter 4. Central Update Algorithm 1 0.201 Residuol Mop S t o n d o r d Deviotion Vs Iteration I ' I ' ' I ' 1 75 1 0.05 [- 0.001 - i 0 i i I 200 ) I i 400 Iteration Number i i I 600 i i i SCO Figure 4.3: Standard deviation of the residual map pixels as a function of iteration during a central update algorithm simulation. The simulated measurements are of an edge only, statistically plausible fake sky map and axe noiseless. Pixels within 150" of the map edge are not used in the computation. At iteration 401, the residual map pixel standard deviation is within 0.02 mJy of the converged value of 0.07 mJy. The steps in the curve between iterations 401 and 801 are due to limited precision when reading and writing data to file. Note the maximum residual map pixel standard deviation of 0.185 mJy near iteration 60. At iteration one, the initial map estimate is the null map, and all pixels used in the computation of the standard deviation are equal to the fake sky map. Thus the standard deviation of the residual map pixels is zero, and rises with iteration as pixels inside the periphery are modified due to flux in the chop positions of the measurements centred upon pixels at the edge of the periphery. Chapter 4. Central Update Algorithm 76 Fake Sky Map Residual 21 200 400 600 200 400 600 Residual 301 Figure 4.4: An edge-only fake sky map and residual maps for a noiseless, full SHADES simulation of the central update algorithm. The fake sky map is populated with sources drawn from a realistic source count model. The edges are 150" wide. The pixel values are in mJy. Note the growth of the edge error. The maximum residual pixel standard deviation occurs around iteration 60. Note how the mean of the residual map pixels grows increasingly negative, reducing the standard deviation. At iteration 101, the absolute pixel values are less than two and a half percent of a 10 mJy source, 360" from the map edge. The white lines are spaced 200 pixels apart, each pixel is 3" wide. Chapter 4. Central Update Algorithm Residual J -0.0007 -0.00065 -0.0006 Map 8 0 1 I -0.00055 77 L -0.0005 -0.00045 (Janskies) Figure 4.5: The 801 residual map for the statistically plausible, edge only, central st update algorithm simulation. The residual due to the edge effect is less than 0.15 mJy, 150" from the map edge. 78 Chapter 4. Central Update Algorithm SHADES measurement strategy. Each panel has 721 columns and 721 rows, which is 2163" wide. These results are encouraging. Provided noise is not greatly affected by the edge-effect, the edge-effect with a statistically plausible fake sky with the full SHADES measurement strategy is negligible compared to the expected statistical residuals. The simulation is repeated with the previous fake sky map except the null map centre is replaced with the centre of the statistically plausible fake sky map. The purpose of this simulation is to demonstrate the relatively fast disappearance of high spatial frequency structure with iteration, and the dominance of the edge effect residual in the converged residual map. (Figures 4.6 and 4.7) Plots of the residual map pixel standard deviation as a function of iteration are shown in Figure 4.6. Convergence to within 0.005 mJy of the converged standard deviation occurs by iteration 50. The converged standard deviation of 0.068 mJy is 0.002 mJy from the the converged standard deviation of the residual pixels in the edge-only simulation (Figure 4.3). The minimum of 0.0655 mJy in the standard deviation of the residual pixels occurs near iteration 40, and is consistent with the rising edge of the curve in Figure 4.3 at iterations less than 25, though the slope of the rise in Figure 4.3 must be less steep. This is purely heuristic, as the presence of sources in the fake map centre modifies the behaviour of the residual map mean with iteration. The 2 standard deviation in Figure 4.3 at iteration 30 is 0.18 mJy. The standard deviation in Figure 4.6 begins high, and then falls as the sky estimate map correctly estimates the sources outside of the map edge. With more iterations, the edge-effect residual propagates into the map and the standard deviation rises from its minimum to the converged value. Residual maps for arbitrary iterations between iteration one and 301 are depicted in Figure 4.7. 300" from the estimate map edges, the sky estimate is correct to within approximately 0.05 mJy, if the map average is estimated at —0.6 mJy, at iteration 90. The residual increases to approximately 0.1 mJy, 150" from the map edge and further increases with proximity to the map edge. Expected SHADES noise is around eight millijanskies per pixel and two millijanskies per beam. Thus, 2 In general, standard deviations of two different random variables do not add. 79 Chapter 4. Central Update Algorithm Residual Map S t a n d a r d Deviation 0.09 0.08 0.07 0.06 400 600 Iteration Number 1000 Figure 4.6: Residual map standard deviation as a function of iteration during a simulation of data reduction using the central update algorithm. The data consist of measurements of a statistically plausible fake sky obtained with the full SHADES measurement strategy. Convergence to within 0.005 mJy of the converged standard deviation occurs by iteration 50. The converged standard deviation of 0.068 mJy is 0.002 mJy from the the converged standard deviation of the residual pixels in the edge-only simulation (Figure 4.3). The minimum of 0.0655 mJy in the standard deviation of the residual pixels occurs near iteration 40. Chapter 4. Central Update Algorithm 80 » • ajgitfi • 1 : • -o.fl Residual 40 -o.6 Residual -o.a 50 mmm HUN 200 400 600 Residual 200 400 I Residual I 80 _1 200 600 -0.7 -0.65 -0.6 -0.55 90 — I 400 -0.5 Figure 4.7: Residual maps for a full SHADES simulation of the central update algorithm. The fake sky map is populated with sources drawn from a realistic source count model. Between iterations 50 and 90 the variance of the residual pixels reduces in the centre of the map and slightly increases toward the map edges. The pixel values are in mJy. Notice the high spatial frequency residuals, visible at iteration one, are largely gone by iteration 40. Chapter 4. Central Update Algorithm 81 in the SHADES survey, a 20% error on the measured flux emanating from a 10 mJy source is expected on 68% of the detected sources. At iteration 90 the algorithm has not converged and the dominant source of residual, over a 0.5 square degree area, is statistical. 4.2 Gaussian Noise Simulation The effect of additive Gaussian noise is determined by simulation. Random, white Gaussian noise is added to each measurement in a full SHADES simulation. The standard deviation of the noise, a, is chosen such that a/VN is around 8 mJy, where N is the number of times that a pixel is centred upon in a measurement. 8mJy is the pixel standard deviation due to noise expected from the SHADES measurement strategy. In this simulation a fake sky map with no sources is measured. Figure 4.8 shows a plot of the residual map standard deviation as a function of iteration. Pixels used in the computation of the residual map pixel standard deviation are located further than 150" from the map edge. The converged residual map pixel standard deviation is 9.95 mJy, 1.95 mJy higher than the expected value. The computation is repeated on the 1000 iteration residual map, using only pixels further than 600" th from the map edge. The resulting residual pixel standard deviation is 9.06 mJy. The reduction in standard deviation is due either to large scale structure in the map or due to noisier pixels at the map periphery. Figure 4.9 shows the first and 101 iteration maps and their power spectrum st maps. The power spectrum maps are invariant to 180° rotation. The lowest spatial frequency is centred in the power spectrum maps, the spectral resolution is 0.617 x 10~ (arcsec ) per pixel and the Nyquist frequency is | of a cycle per arcsecond. 3 -1 The 101 power spectrum map contains a bright cross centred in the map, as well st as a fainter, checkered pattern, corresponding to the chop directions. Along the right half of the middle row there are a number of bright peaks. The first four of these bright peaks, moving right from the centre of the 101 power spectrum map, st correspond to spatial wavelengths at 34.5", 22.8", 15.9", and 14.0", respectively. The first three of these wavelengths are roughly one-half of the chop throw sizes. The Chapter 4. Central Update Algorithm 82 Residual Map Standard Deviation 10.0 400 600 Iteration Number 800 Figure 4.8: Residual map average and standard deviation as a function of iteration during a simulation of data reduction using the central update algorithm. Gaussian random noise is added to simulated measurements of a fake sky containing no flux. At iteration 100, the residual pixel standard deviation is 1.95 mJy higher than the expected (7 asurement/v^/V The me central-update algorithm converges when reducing measurements with Gaussian additive noise. Chapter 4. Central Update Algorithm 200 200 400 -0.02 -0.01 1 L 1 1 0 0.01 PS I t e r a t i o n 1 83 •i 00 1 0.02 PS I t e r a t i o n 101 • m• i i 1 u• oilm si• • 11 HE ma • •• • • r»* • ji$jk j 4' * ! 100 200 300 400 500 100 I le-09 "1 2e-09 1 3e-09 200 300 400 500 I 1 4e-09 5e-09 Figure 4.9: Measurements with Gaussian noise are iteratively reduced using the central update algorithm. The maps are, clockwise from upper left, the 1 st iteration map, the 101 iteration map, the power spectrum of the 1 itst st eration map and the power spectrum of the 101 and 1 iteration map. st st Note the structure that appears in the 101 iteration maps. st Chapter 4. Central Update Algorithm 84 timestream noise covariance matrix is diagonal in this simulation, and the pixel-pixel correlations are due entirely to the measurement strategy. The power spectrum is the Fourier transform of the correlation function; structure in the power spectrum map is due to pixel-pixel correlations. Figure 4.4, by juxtaposing the first and 101 st iteration power spectrum maps, illustrates the growth of pixel-pixel correlations formed during the convergence of the iterative deconvolution algorithm. Examining the sky estimate maps we see a cross-hatch pattern and large scale fluctuations present in the 101 iteration map that are not present in thefirstitst eration sky estimate map. As mentioned in the previous paragraph, this structure appears in the power spectrum maps and is correlated pixel-pixel noise. An intuitive understanding of this structure is possible. The SCUBA double-differencing measurement strategy, for scales much larger than the chop distance, is approximately a double spatial derivative of the sky flux. A second derivative is a high-pass filter. Thus, deconvolving the chops requires amplifying the low spatial frequencies to undo the filter, resulting in large spatial scale noise amplification. The amplification of the large scale noise will increase the residual map pixel standard deviation. The simulation is repeated with a different fake sky, this time containing four, 30mJy sources and a single 50mJy source. Figure 4.10 shows the portion of the 1 and 101 sky estimate maps, and their corresponding power spectrum maps, st st which contain the large sources. Again, the 101 iteration map exhibits crossst hatching and large scale fluctuations. The regions surrounding the sources exhibit no correlations not exhibited by regions far from the sources in the same map. Though not quantified, the pixel-pixel noise correlations in the iterated map are not visibly affected by the presence of a 50 mJy source. The simulation is repeated a final time, this time using the statistical model of SCUBA sources to create the fake sky. The sky estimate maps and power spectrum maps for iteration one and 101 are shown in Figure 4.11. The iteration one power spectrum map contains new structure around the centre of the map. This structure is at wavelengths greater than 34.5" and is due to the presence of many sources in the map. The 101 power spectrum map appears very similar to the previous 101 st st power spectrum maps, the cross-hatching, and the bright peaks are all present, in Chapter 4. Central Update Algorithm Iteration 1 300 350 85 I t e r a t i o n 101 400 450 300 350 400 450 Figure 4.10: No source streaking, other than the existing cross-hatch pattern, is apparent, in a simulated reduction of measurements of four large sources containing Gaussian additive noise using the central update algorithm. The pixel values are in mJy. Chapter 4. Central Update Algorithm Iteration 1 200 I t e r a t i o n 101 400 600 K -0.02 200 i -0.01 0 l le-09 l 0.01 Power Spectrum 1 0 86 400 600 0.02 Power Spectrum 101 2e-09 3e-09 4e-09 5e-09 Figure 4.11: Measurements of a statistically plausible fake sky with additive Gaussian noise are iteratively reduced using the central update algorithm. The maps are, clockwise from upper left, the 1 iteration map, the st 101 iteration map, the power spectrum of the 1 iteration map and st st the power spectrum of the 101 iteration map. The presence of sources st appears as low frequency power in the 1 iteration power spectrum st map. Chapter 4. Central Update Algorithm 87 the same locations, at the same values. 4.3 More Realistic Noise Simulation To test the central update algorithm in a more realistic application, a simulation of the algorithm is performed with more realistic noise. Realistic noise is created by extending SCUBA data taken during three observations of a blank field to the 4.2 x 10 measurements required in the SHADES simulation. The extension is performed 7 as follows: the existing blank field measurements are divided into equally sized groups of contiguous measurement samples. The original blank field measurements are then extended by adding a randomly selected contiguous piece of data, drawn from the group of contiguous, existing data subsets. This process is continued, until a sufficient number of data samples is obtained. Figure 4.12 shows example timestream noise samples for four arbitrarily chosen bolometers and the magnitude of their Fourier transforms. Each bolometer timestream is pure noise and has 21,120 samples, requiring an extension factor of 22. The extension factor is arbitrarily chosen. The time streams and Fourier transform magnitudes have maximum values consistent with the timestreams and Fourier transform magnitudes of real, preprocessed SHADES data plotted in Figure 2.5. The randomly extended timestreams lack the unstationarity of the actual SHADES data and there are spikes in the 3 Fourier transform magnitudes not present in the real SHADES data. These spikes are likely due to thefixedsize of the timestream subsets used in the random extension of the timestreams. The randomly extended noise is added to simulations of the fake sky map with a statistically plausible fake sky in a full SHADES simulation of the central update algorithm. Plots of the residual map standard deviation are shown in Figure 4.13. At iteration 100, the standard deviation of the residual map pixels is approximately 0.35 mJy from the converged value of 11.55 mJy. Repeating the calculation for the 1000 residual map omitting the outer 600" results in a residual pixel standard th deviation of 8.89 mJy. This value is much nearer the expected value of 8mJy. The 3 The variance of each bolometer timestream in Figure 2.5 changes with time. That is, the data samples are drawn from time-dependent probability density functions. Chapter 4. Central Update Algorithm 88 Fourier Transform MoqnitudJ Fourier Tranitorm Mognituflc r™j«nej (tu) Figure 4.12: First four plots are of the randomly extended noise sequences for bolometers 1, 2, 9, and 35. The bolometer selection is arbitrary. Last four plots are the corresponding Fourier transform magnitudes. Note the lack of j noise in the spectrums. There is a bump at approx. 0.007 Hz in the bolometer 35 spectrum, possibly the Fourier transform spike described in chapter 2. Chapter 4. Central Update Algorithm Residual Map S t a n d a r d 200 89 Deviation 400 600 Iteration Number 800 1000 Figure 4.13: Residual map standard deviation as a function of iteration in a central update algorithm simulation with a statistically plausible fake sky map and the full SHADES measurement strategy. Chapter 4. Central Update Algorithm 90 larger decrease of the residual pixel standard deviation between this simulation and the simulated reduction of Gaussian distributed measurements of a blank sky is due either to greater contributions to the standard deviation from the map edges or due to a greater contribution to the pixel variance from large scale structure. By looking at the 101 iteration sky estimate map in Figure 4.14 and comparing st to the 101 iteration sky estimate map in Figure 4.9, we see that there is greater st large scalefluctuationin Figure 4.14, especially at the edges. This is consistent with the greater reduction in the residual pixel standard deviation with decreasing edge contributions to the calculation. Figure 4.14 depicts the first and 101 iteration maps with their corresponding st power spectrum maps. The maps in Figure 4.14 closely resemble those produced in the simulations of the central update algorithm with measurements containing Gaussian additive noise. In this simulation, it is unlikely that the timestream noise covariance matrix is diagonal, and one would expect to see extra pixel-pixel correlations. From Figure 4.15 this is not evident, and if present, has a relatively small effect. From iteration 100 to convergence, there is a 0.4 mJy change in the standard deviation of the residual map pixels (See Figure 4.13). Though this change is small, any outlier pixels are greatly suppressed in the computation due to the normalization by a large number of pixels. It is known that the algorithm will converge, neglecting edge-effects, to an optimal solution. In noiseless simulations, by iteration 90 residuals span hundreds of arcseconds and are below 0.1 mJy (See Figure 4.7). One wishes to cease iteration when improvements due to further convergence are negligible. Iteration amplifies the large scale noise and correlates the pixels, while removing the negative side-lobes present around sources in the first iteration sky estimate map. In principal, the converged map is optimal; however, subsequent optimal use of the converged map will require accounting for the pixelpixel correlations. The degree to which continued iterations between iteration 200 and iteration 1000 affect the science product will only be determined by simulation involving the method of accounting for the pixel-pixel correlations, and will depend on the objective of the analysis. The pixel-pixel correlations are complicated and it is not intuitive how their presence or enhancement will benefit, for instance, source Chapter 4. Central Update Algorithm Iteration 1 Iteration -0.02 -o.oi o 0.01 Power Spectrum 1 100 200 300 400 91 101 0.02 Power Spectrum 101 500 100 200 300 400 500 Figure 4.14: Clockwise from upper left: Sky estimate maps for iteration one, 101, power spectrum map for iteration 101, power spectrum map for iteration one. Note the similarity between these maps and the maps produced in simulations with Gaussian additive noise. (Figure 4.12) The off-diagonal terms of the timestream noise-covariance matrix play a minor role in the pixel-pixel correlations present in the 101 iterast tion map. The pixel-pixel correlations are predominantly due to the off-diagonal terms in the P P matrix for the SHADES measurement T strategy. The P P matrix was introduced in chapter two. T Chapter 4. Central Update Algorithm 92 detection. The next section demonstrates the effects of measurement strategy on pixel-pixel correlations by modifying the SHADES chopping scheme in simulation. 4.4 Chopping Effects On Pixel-Pixel Correlation In this section the effect on the pixel-pixel correlations is investigated through simulation by increasing the number of chops to 12. Specifically, one half of the measurements are performed along chops rotated 45 degrees from the orignal chop directions, adding six new chops to the measurement strategy. Figure 4.15 shows the first, and 101 sky estimate maps and their corresponding power spectrum maps. The total st number of measurements and the additive noise are the same as in the randomly extended noise simulations described in section 4.3. The fake sky is statistically plausible. The 101 sky estimate map in Figure 4.15 exhibits less cross-hatching st than the 101 sky estimate map in Figure 4.11. The large scale structure, with a st half-wavelength aroung 600", seems unaffected. Examining the corresponding 101 st power spectrum maps, we see a reduction in power along the very bright cross centred in Figure 4.11, and a change of the cross-hatch pattern to a star pattern. For both power spectra, the pattern corresponds to the direction of the chops. The power is distributed evenly in all directions in an inner square 100 pixels wide centred in the 101 power spectrum map in Figure 4.15. This distribution of correlations along st new directions is consistent with the arguement presented in Section 6.3 regarding the reduction of pixel-pixel correlations by involving every pixel in measurements involving many different pixels. 4.5 Chapter Results The central-update algorithm, as applied to SHADES, is evaluated has been evaluated by simulation. The residual due to the edge effect is found to be negligible compared to the statistical noise in the SHADES survey for most of the central half-square degree of the sky estimate map (Section 4.1). The presence of strong sources in the measured fake sky map, does not correlate the noise. (Section 4.2) For Chapter 4. Central Update Algorithm 93 a statistically plausible fake sky with randomly-extended realistic noise, the standard deviation of the residuals upon convergence is 11.56 mJy, 0.4 mJy higher than the standard deviation of the residual map pixels at iteration 101 (Section 4.3). In Section 4.4 the effect of changing the SHADES chopping pattern is investigated and the pixel-pixel correlations found to be significantly reduced. Besides further simulations, to use the iteratively deconvolved sky estimate map, either the pixel-pixel noise covariance matrix must be used with subsequent source detection, or the pixel correlations must be reduced. In the chapter on Discussion, the effect of follow-up observations aimed at reducing the pixel-pixel correlation is addressed. As well, the use of priors to constrain the correlation is briefly mentioned in the section on Noise & Correlation Removal. Chapter 4. Central Update Algorithm J Iteration I 1 I -0.02 I -0.01 0 0.01 94 Iteration 1 101 L 0.02 (Jansky) 100 200 300 400 500 100 200 300 400 500 Figure 4.15: Simulation of data reduced using the central update algorithm. The data are simulated measurements of a statistically plausible fake sky with blank field extended noise samples. One half of the chops are rotated 45 degrees from the east-west and north-south chop directions. Much of the cross-hatching visibile in the 101 iteration map in Figst ure 4.14 has disappeared. The strong cross-shape in the 101 power st spectrum map in Figure 4.14 has been reduced. Large scale structure with half-wavelengths around 600" seem unaffected. Chapter 5. SHADES Subaru XMM-Deep Field 95 Chapter 5 SHADES Subaru XMM-Deep Field A list of sources detected in a measurement map of the Subaru XMM-Deep Field is presented. The reduced data includes all SHADES observations up to the 1 of st Febuarary, 2004. 5.1 Data Reduction The data is pre-processed using the techniques described in Section 2.4, with the exception that the flux conversion factors are set to a nominal value of 225Jy V - 1 Figure 5.1 shows the measurement map of the Subaru XMM-Deep Field. Each pixel has a standard deviation of approximately 8 mJy, and is fully measured for most of the map. Figure 5.2 shows the hit map for the Subaru XMM-Deep Field. The value of each pixel in the hit map is the number of times each pixel is centred upon in a measurement and the noise maps are computed assuming a diagonal timestream noise covariance matrix. That is, the p h pixel in the noise map is it (5.1) where N is the number of measurements centred upon the i pixel and o is the th m n standard deviation of the set of measurements containing the n measurement. The th set of measurements is grouped by bolometer and observation file, and nominally contains 320 measurements. Chapter 5. SHADES Subaru XMM-Deep Field 96 Subaru XMM-Deep F i e l d Measurement Map 500 1 a» ii 400 300 • - ' 200 200 -0.03 400 300 -0.02 -0.01 I 0 500 I 0.01 0.02 0.03 (Jansky) Figure 5.1: The Subaru XMM-Deep Field measurement map. Units are in Jy, pixels are 3" square. The horizontal and vertical axes are specified in units of pixels. Chapter 5. SHADES Subaru XMM-Deep Field 97 Subaru XMM Deep F i e l d : Noise Map J ••fWlaiiiSiJi L L fa J r -h i n 200 0.006 1 1 300 400 0.008 0.01 r r 500 0.012 0.014 Subaru XMM Deep F i e l d : H i t Map * k mm iwiJI 200 0 300 50 400 100 500 150 Figure 5.2: Subaru XMM-Deep Field noise map and hit map, top and bottom, respectively. The noise map units are in Jy and the noise per pixel is between 7mJy and 11 mJy for the portions of the map fully measured by the SHADES strategy measurement strategy. The SHADES tripos pointing pattern is clearly visible. Chapter 5. SHADES Subaru XMM-Deep Field 5.2 98 Source Detection Source detection can be posed as the task of determining the model which will provide the best fit to the pixel values. Ideally one would determine the most likely model by computing the probabilities of different models by marginalizing over all model parameters. Proper computation of these probabilities will intrinisically incorporate an Occam's penalty, and will select the number of sources warranted by the data.[? , BayesTutorialjhe candidate models are no sources present, one source present, two sources present, etc. Because the practical solution to this problem is difficult, a more simple approach is taken.[5] [19] [18] The method is as follows: 1 for every pixel in the map, assume that there is a source centred upon the pixel under consideration, assume that the prior probabilities are uniformly distributed , the pixel noise is Gaussian and the source full-width half-maximum is given by the beam of the telescope. Compute the most probable amplitude for the source by finding the source amplitude which maximizes the likelihood. That is, given that there is a source centred upon the pixel under consideration, the likelihood is P(d| )ocexp(-^) (5.2) S where P (d\s) is the probability of the data, d, given the model, s, and fhzMt N Xi = i where j is the j t h (5 .3) 1 pixel under consideration (centred-upon), pi, is the ? pixel, /3 th is the source amplitude, fij, is the value of the source model centred upon pixel j, for the i th pixel, and o is related to the telescope beam and determines the beam full-width half-maximum. The most likely amplitude, J3, is obtained byfindingthe maximum of x and, for the j 2 th pixel, is given by A 1 The = ^^fr (5-4) c o m p u t a t i o n requires for every m o d e l considered a n integral of t h e l i k e l i h o o d over a l l source centres a n d a l l source a m p l i t u d e s allowed b y t h e priors. 99 Chapter 5. SHADES Subaru XMM-Deep Field in practice, will be negligible for i's significantly different than j, and the sums involving i over all of the map pixels can be limited to the contributing terms. Once the most probable amplitude is computed, the next thing to do is to compute the uncertainty in the amplitdue estimate. This is accomplished by taylor expanding the x about the most probable amplitude, and taking the second, and 2 final, term to be the uncertainty. Because the taylor expansion is about the maximum, the first order term is zero. The second order term gives an estimate of the width of the Gaussian likelihood. The uncertainty, is given by -0.5 (5.5) The source model, fij, used in the above computation is the measurement map formed from measuring a fake sky map containing a unitless source with amplitude one with the SHADES strategy. Figure 5.3 is a map of the mask. The next step in the algorithm is to compute the amplitude uncertainty using equation 5.4, and compute the ratio between the most probable amplitude and the amplitude uncertainty. This process is repeated for every pixel in the map. The most-probable amplitude map, uncertainty map, and the ratio map for the Subaru XMM-Deep Field are shown in Figures 5.4. Once the ratio map is created, the maximum pixel in the ratio map is found. A source, with an amplitude obtained from the corresponding location in the most probable amplitude map, is subtracted from the measurement map. The subtracted source is considered detected at the pixel value in the ratio map, and is recorded. The most probable amplitude, uncertainty, and ratio maps are re-computed, and the process repeated until there are no pixels greater than three sigma left in the ratio map. The list of recorded sources will be all sources greater than three standard deviations from the map mean of zero. These sources are considered detected at greater than three sigma. To check the source detection a histogram of the pixels in the uncleaned most probable amplitude map is compared with a histogram of the pixels in the source cleaned most probable amplitude map. The histograms are plotted in Figure 5.5, along with a scaled probability density function corresponding to the probability of the union of the pixel probabilities. The scaled probability density function Chapter 5. SHADES Subaru XMM-Deep Field 100 Figure 5.3: Mask used in source detection tofitthe amplitude of a point centred on each pixel. Note the negative sidelobes. The central Gaussian has a full-width half-maximum of 14.7". Each pixel is 3" square. The mask is 150" square. The mask mean pixel intensity is zero. Chapter 5. SHADES Subaru XMM-Deep Field Most P r o b a b l e A m p l i t u d e Map 101 A m p l i t u d e U n c e r t a i n t y Map >.002 0.003 n Figure 5.4: Subaru X M M Deep Field: Clockwise from upper-left, most-probable amplitude map, amplitude uncertainty map, and ratio map. The source detection acts as a statistically weighted low-pass filter. The funny pattern along the edges of both the most probable amplitude map and the ratio map is due to the averaging of only a few noisy values in the assignment of the pixel values near the edges. The saturated, mask sized regions in the upper right and lower left of the most probable amplitude map and the ratio map are due to very large valued pixels along the edge of the Subaru measurement map. Chapter 5. SHADES Subaru XMM-Deep Field 102 represents the expected histogram made by histogramming the pixels of a pure noise most probable amplitude map. The probability density function is constructed by summing the probability density functions for each of the pixels. Each pixel probability density function is given by a Gaussian whose width is specified by the corresponding pixel value in the amplitude uncertainty map. The sum of the varying width Gaussians gives the probability density function for the union of the most probable amplitude map pixels. The union probability density function is scaled by eye. Near the edge of the map the noise properties of the map have not been properly investigated. The source detection mask (Figure 5.3) is a measurement map of a single source and contains negative sidelobes, the largest of which is centred 68" from the source centre. Because of this, the mask is large and much data near the edges of the map would be lost if the convolution were performed by truncating the edge of the most probable amplitude map by half the mask size. Instead, the convolution is performed as an average, such that the mask does not have to lie completely within the map prior to getting a non-zero value. When the mask is only partly within the map, the computation is more susceptible to a few, very noisy pixels. This is clear in the edges of the most probable amplitude map depicted in Figure 5.5. If histograms of the noise were produced of the edge pixels and compared with a noise model constructed from the uncertainty map, then detections at map edge would be included. This test has not been done and the author chooses to remove sources with a chop substantially off the edge of the map from the list of detected sources. To check this removal, a map of the Subaru XMM-Deep Field with a cross drawn over the source positions is created. Figure 5.6 depicts the map. The lists of detected sources are given in Tables 5.1, 5.2 and 5.3, in three parts, for the Subaru XMM-Deep Field. The columns and rows of the sources corresponding to the locations in the maps shown in this chapter are included in the tables. The correspondence of the Subaru XMM-Deep Field source list with the list given by the SHADES community as of April 12, 2004 is demonstrated in Figure 5.7. Considerable non-correspondence exists between the sources detected by the two groups. The expected non-correspondence is small as the data prior to reduction Chapter 5. SHADES Subaru XMM-Deep Field 103 1OOOO F -0.02 -0.01 0.00 0.01 0.02 Pixel Value (mjy) Figure 5.5: Subaru XMM Deep Field most probable amplitude pixel histograms with noise model. Top, histograms produced from pixels with noise less than 4mJy. Bottom, histograms made from pixels with noise less than 3mJy. There is an excess of positive, and negative flux in the uncleaned most probable amplitude map, due to the presence of sources. The apparent lack of false detections in the top panel is due to the greater pixel noise populating the top histogram's positive tail. Chapter 5. SHADES Subaru XMM-Deep Field 104 srcOnRatioBlack.fits_0 __l I (SNR) Figure 5.6: Subaru XMM Deep Field source detections. The black crosses mark the location of detected sources. The closest detections are around 70" from the map edge. Chapter 5. SHADES Subaru XMM-Deep Field 105 S r c s Detected With SHADES Data Up To Feb 1, 2 0 0 4 1 , 1 | , , 1 , • 0 + - 5 + + «\ + + • T o o o o + + o o o + + + O + - o o 0 •r + + 1 + _ + - -10 O 0 -5 Rest of SHADES . Hi < * o o 0 < o o 1 - o * + o + o 1 0 0 + 1 - + i * " o 1 UBC - o + o - 1 + , 1 , 1 1 -5 o 1 0 1 1 , , , , 5 10 RA Offset (arcmins) Figure 5.7: Subaru XMM-Deep Field source list correspondence. There are 26 unpaired UBC detections, 25 unpaired detections from the SHADES community, and 35 matches. 106 Chapter 5. SHADES Subaru XMM-Deep Field ^ 4.4 in 4.2 4.0 0 100 200 Iteration 300 400 Figure 5.8: Sky estimate map pixel standard deviation as a function of iteration. The standard deviation is within 0.01 mJy of its converged value at iteration 100. are identical. That is, in principle, the source lists should only differ by differences in the reduction and source detection algorithms. The scatter may be due to improper calibration, or to an unknown observational problem. Different source detection algorithms are used between the groups, and mild disagreement on weakly detected sources is plausible. 5.3 Deconvolved Map Source Detection As a means of quantifying the usefulness of the deconvolved map, the central update algorithm is used to reduce the Subaru XMM-Deep Field data, and subsequent, nonoptimal source detection is performed. Plots of the standard deviation of the map pixels as a function of iteration are shown in Figure 5.12. The deconvolved map pixel standard deviation is within 0.01 mJy of its converged value at iteration 100. Figure 5.13 is the 401 iteration sky estimate map. st The large scale structure and cross-hatching seen in simulation in chapter four is evident. Measurements centred upon pixels involved in the centre of measurements less than thirty times Chapter 5. SHADES Subaru XMM-Deep Field 107 Figure 5.9: 401 sky estimate map. The pixel-pixel correlations are apparent. Meast surements centred upon pixels involved in the centre of measurements less than thirty times are discarded. Source detection on this map will have to deal with the large scale structure, and if it is to be optimal, the cross-hatching. Chapter 5. SHADES Subaru XMM-Deep Field 108 are discarded. The poorly measured edges in this map affect the stability of the centre update algorithm and increase the pixel-pixel correlations. Discarding poorly measured pixels from the data reduction, which means discarding the measurements of the poorly measured pixels, increases the stability of the algorithm and decreases the pixel-pixel correlations present in the resulting map. Source detection is performed on the 401 iteration map. The choice is arbitrary, st any map after iteration 100 would suffice. To remove some of the bias due to the large scale structure, the source detection is performed with the same mask used in the source detection performed on the 1 iteration map. This mask has a zero st mean and is evenly symmetrical. Odd functions and constant offsets will not be measured; however, fluctuations possessing even symmetry or no symmetry on the scale of the mask size will result in non-zero response and affect source detection. The source detection is performed using the algorithm described in section 5.2. Figure 5.14 shows the correspondence between the sources detected on the deconvolved map and the sources detected in section 5.2 for the Subaru XMM-Deep Field as well as the correspondence with the source list provided by the rest of the SHADES community. The source list resulting from the source detection on the deconvolved map has not been pruned to reject detections at the edge of the map. The correspondence between any of the two source lists is equal. There is no indication which list is more reliable. The optimal source detection, taking into account the full pixel-pixel correlation matrix, is expected to out-perform the source detection applied in section 5.2 to the measurement map. The correlation matrix will need to be large enough to account for the dominant correlations, and will need to be inverted. It is possible that the matrix will take a form amenable to fast inversion; however, the scrambling of the time series will likely result in a pixel- pixel correlation matrix lacking nice symmetries. Nonetheless, the investigation has not been done and the author is speculating. Chapter 5. SHADES Subaru XMM-Deep Field i -| 10 109 i i r + Iteration 1 O Iteration 4 0 ' o + O A Rest of Shadejs A A A A o * o O A UJ Q A -5h & A & A + + A -F A+ A A A -k> -101 o A _i i i i_ 10 RA ( a r c s e c s ) Figure 5.10: Deconvoluted source detection correspondence with the sources detected on the measurement map (1 iteration map) and with the source st list provided by the rest of the SHADES community using the data acquired up to February 1 , 2004. There is no indication which source st detection method is more reliable. Chapter 5. SHADES Subaru XMM-Deep Field 110 Subaru X M M Deep Field Source Detections - Part I RA (°) DEC (°) Significance Amplitude (mJy) Uncertainty (mJy) Col Row 34.3779 -4.9934 7.04 14.2 2.0 506 368 34.5151 -4.9242 5.66 11.9 2.1 342 451 34.4791 -4.8834 4.65 10.5 2.3 385 500 34.4130 -5.0892 4.54 9.1 2.0 464 253 34.4356 -4.9316 4.52 10.1 2.3 437 442 34.4255 -4.9417 4.56 9.6 2.1 449 430 34.5653 -4.9008 4.45 9.1 2.1 282 479 34.4983 -5.0842 4.16 8.5 2.0 362 259 34.4682 -5.0825 4.12 9.0 2.2 398 261 34.4289 -5.0742 4.08 9.1 2.2 445 271 34.5009 -5.1283 3.85 9.1 2.4 359 206 34.3762 -4.9975 3.76 6.9 1.8 508 363 34.4029 -5.0766 3.75 6.9 1.9 476 268 34.5787 -5.0467 3.69 9.2 2.5 266 304 34.4832 -5.1067 3.69 8.1 2.2 380 232 34.3762 -5.0167 3.67 6.7 1.8 508 340 34.4289 -5.0967 3.67 7.5 2.0 445 244 34.3695 -4.9850 3.61 7.4 2.1 516 378 34.4339 -5.0375 3.61 7.8 2.2 439 315 34.4113 -5.0608 3.60 6.8 1.9 466 287 34.3528 -4.9791 3.56 7.7 2.2 536 385 34.6054 -4.9541 3.52 9.2 2.6 234 415 34.5217 -4.9875 3.49 10.1 2.9 334 375 Table 5.1: Subaru X M M Deep Field Source Detection List - Part I Chapter 5. SHADES Subaru XMM-Deep Field 111 Subaru X M M Deep Field Source Detections - Part II RA (°) DEC (°) Significance Amplitude (mJy) Uncertainty (mJy) Col Row 34.5310 -4.9875 3.82 11.5 3.0 323 375 34.5636 -5.0308 3.45 7.8 2.3 284 323 34.4866 -4.8775 3.43 9.5 2.8 376 507 34.4022 -4.9325 3.42 7.7 2.3 477 441 34.4063 -5.0658 3.39 6.6 2.0 472 281 34.5326 -5.0675 3.39 7.6 2.3 321 279 34.3878 -5.0333 3.38 6.2 1.8 494 320 34.4189 -5.0217 3.36 7.7 2.3 457 334 34.6346 -4.9951 3.34 8.6 2.6 199 366 34.5267 -4.9258 3.33 6.7 2.0 328 449 34.5577 -4.9616 3.29 6.3 1.9 291 406 34.4415 -4.9650 3.28 6.5 2.0 430 402 34.4648 -5.0467 3.27 6.6 2.0 402 304 34.5686 -4.9200 3.25 7.2 2.2 278 456 34.3870 -4.9708 3.25 6.6 2.0 495 395 34.3879 -4.9525 3.25 7.0 2.2 494 417 34.4615 -5.1092 3.22 6.9 2.1 406 229 34.5151 -5.0708 3.22 7.3 2.3 342 275 34.4950 -4.9083 3.21 6.3 2.0 366 470 34.5041 -4.8867 3.19 7.2 2.3 355 496 34.3553 -4.9942 3.18 6.6 2.1 533 367 34.4515 -5.0008 3.15 8.1 2.6 418 359 34.5469 -5.0467 3.15 6.6 2.1 304 304 Table 5.2: Subaru X M M Deep Field Source Detection List - Part II Chapter 5. SHADES Subaru XMM-Deep Field 112 Subaru X M M Deep Field Source Detections - Part III RA (°) DEC (°) Significance Amplitude (mJy) Uncertainty (mJy) Col Row 34.3737 -4.9875 3.13 6.5 2.1 511 375 34.3921 -4.9825 3.15 5.6 1.8 489 381 34.5652 -4.9316 3.12 6.4 2.1 282 442 34.5845 -4.9458 3.11 6.9 2.2 259 425 34.5619 -4.9891 3.10 6.4 2.1 286 373 34.4239 -4.9758 3.10 6.2 2.0 451 389 34.4106 -4.8783 3.10 7.9 2.6 467 506 34.5452 -5.1292 3.09 10.5 3.4 306 205 34.4607 -4.9275 3.09 6.0 2.0 407 447 34.3712 -5.0458 3.09 7.3 2.4 514 305 34.4255 -4.9033 3.03 6.3 2.1 449 476 34.4214 -5.1341 3.03 9.1 3.0 454 199 34.4539 -5.0866 3.01 7.1 2.3 415 256 34.3946 -5.0775 3.00 6.7 2.2 486 267 34.4456 -5.0750 3.00 6.6 2.2 425 270 34.4673 -5.0783 3.04 7.0 2.3 399 266 Table 5.3: Subaru X M M Deep Field Source Detection List - Part III 113 Chapter 6. Discussion Chapter 6 Discussion This chapter summarizes the main results of the thesis, discusses follow-up observation aimed at reducing the pixel-pixel correlations in the deconvolved map, discusses general issues when choosing a measurement strategy, poses the problem of determining an optimal measurement strategy, and discusses future theoretical work regarding the reduction of pixel-pixel correlations in the deconvolved map. 6.1 Main Results Data calibration using a flux conversion factor computed from calibration observations of CRL618, Uranus and Mars is fournd to have a 10% difference in the correspondence of detected sources compared to the sources detected with the same data except using a flux conversion factor equalling 225 JyV . The chop sizes of -1 the differenct calibration observations is ignored (Section 2.4.1). The iterative deconvolution algorithm is generalized and formulated as the solution to a difference equation which gives the map for any iteration (Equation 3.35). The generalized deconvolution algorithm is shown to diverge in SHADES simulations. The convergence properties of the generalized iterative algorithm are determined by the eigenvalues of a characteristic matrix, with small dependence upon the measurements. The formulation is checked theoretically and by simulation. The iterative algorithm is theortically shown to converge to the least-squares optimal solution. (Section 3.59) A known work-around, the central-update algorithm, is explained by the theory (Section 3.5.3). The central-update algorithm, as applied to SHADES, is evaluated by simulation (Chapter 4). The residual due to the edge effect is found to be negligible compared to the statistical noise in the SHADES survey for most of the 6. Discussion Chapter 114 central half-square degree of the sky estimate map (Section 4.1). For a statistically plausible fake sky with randomly-extended realistic noise, the standard deviation of the residuals upon convergence is 11.56 mJy, 0.4 mJy higher than the standard deviation of the residual map pixels at iteration 101 (Section 4.3). In Section 4.4 the effect of changing the SHADES chopping pattern is investigated and the pixel-pixel correlations found to be significantly reduced. In chapter 5 a source list is presented for source detections in the Subaru XMMDeep Field using the first iteration, or measurement map. The source list is compared with a list of sources provided by the SHADES community. Source correspondence is around 50%. 6.2 Follow-Up Observation The effect of follow-up observation on the iterated map pixel-pixel correlations is investigated by modifying the SHADES chopping strategy in simulation. Two simulations are performed, the first with one extra chop, constituting an increase in the number of SHADES measurements by ^ the original number, and a second th simulation with two mutually orthogonal chop directions of the same size, increasing the number of SHADES measurements by | . th Figure 6.1 shows the first and 101 * sky estimate maps and their corresponding s power spectra, for a simulation with one extra chop per pointing. The chop is oriented at a 45 degree angle to the existing chop directions and has a size of 44". The extra chop per pointing constitutes an extra j j r of the SHADES observations. t h The fake sky is statistically plausible and the noise consist of randomly extended blank sky samples described in section 4.3. The power along the dominant cross in the 101 power spectrum map in Figure 4.14 is somewhat reduced, and the crossst hatch pattern modified. The maps are less noisy as there are more measurements. Despite this, the pixel-pixel correlations in the 101 power spectrum map are larger st than in the 12 chop simulation resulting in the 101 power spectrum map shown in st Figure 4.15. The effect of two extra chops is investigated. The simulation conditions are the Chapter 6. Discussion Iteration 1 115 I t e r a t i o n 101 -0.02 -0.01 0 0.01 0.02 (Jansky) Power Spectrum 1 100 200 300 400 Power Spectrum 101 500 100 200 300 400 500 Figure 6.1: Simulation of data reduced using the central update algorithm. The data are simulated measurements of a statistically plausible fake sky with blank field extended noise samples. An extra chop's worth of data has been added to measurements made with the full SHADES measurement strategy. The chop has a size of 44" and is rotated 45 degrees with respect to the original chop directions. Clockwise from upper left, first iteration sky estimate map, 101 iteration sky estimate map, 101 st st iteration power spectrum map, first iteration power specrum map. The power along the dominant cross in the 101 power spectrum map in st Figure 4.14 is somewhat reduced, though not as much as in Figure 6.1, and the cross-hatch pattern modified. 116 Chapter 6. Discussion same as in the above simulation except two extra chops for each pointing are added instead of one. The extra chops are rotated 90 degrees with respect to each other and are rotated 45 degrees from the east-west, north-south chop directions. Figure 6.2 shows the sky estimate maps and their corresponding power spectra for the first and 101 iteration maps. Power along the prominant cross in the 101 power spectrum st st maps in Figures 6.1 and Figure 4.14 is further reduced. At low frequencies the power has been concentrated to narrow lines along the chop directions. The cross-hatch pattern throughout the 101 power spectrum map is changed and reduced. st It is clear that follow-up observation with different chops reduces the pixel-pixel correlations present in the data. This effect should be two-fold; both a reduction in the pixel standard deviation by y/~N and the effect of the extra cross-linking described in the next section. However, the enhancement has a negligible effect on the large scale fluctuations in the deconvolved map (Compare the 101 sky st estimate maps in Figures 4.14, 6.1, and 6.2) and modestly improves the pixel-pixel correlations, while the required observation time is long. 6.3 Measurement Strategy The SHADES measurement strategy is not optimal when producing deconvolved maps. To produce a diagonal pixel-pixel noise covariance matrix, both the timeline noise covariance matrix and the P P matrix must be diagonal. Absolute power 1 measurements produce a diagonal P P matrix but often have a non-diagonal noise T covariance matrix. Many sources of noise correlation in the timelines are reduced by performing differential measurements; however, the reduction is traded for a diagonal P P. Experiments with good measurement strategies, such as WMAP, achieve T nearly diagonal pixel-pixel noise covariance matrices by taking differential measurements involving every pixel with as many other pixels in the map as possible, and by repeating the entire experiment many times. Each differential measurement correlates all of the pixels involved in the measurement, while increasing the appropriate diagonal term in P P by, for instance, one in the central update algorithm scheme T with noiseless measurements, and the off-diagonal terms by a small amount. By Chapter 6. Discussion I t e r a t i o n 101 Iteration 1 ~\ -0.02 117 -0.01 r 0 0.01 0.02 (Jansky) Power Spectrum 101 Power Spectrum 1 le-09 ~1 2<e-09 1 3e-09 1 4e-09 1 Se-09 Figure 6.2: Simulation of data reduced using the central update algorithm. The data are simulated measurements of a statistically plausible fake sky with blank field extended noise samples. Two extra chop's worth of data has been added to measurements made with the full SHADES measurement strategy. The chops have a size of 44" are rotated 45 degrees with respect to the original chop directions and are mutually orthogonal. Clockwise from upper left, first iteration sky estimate map, 101 iteration sky estimate map, 101 iteration power spectrum map, st st first iteration power specrum map. Power along the prominant cross in Figures 6.1 and Figure 4.14 is further reduced. At low frequencies the power has been concentrated to narrow lines along the direction of the prominant cross. Chapter 6. Discussion 118 involving many different pixels in many measurements of the same pixel, the diagonal terms grow with respect to the average off-diagonal term. Correlating every pixel with every other pixel, a map relatively free of pixel-pixel correlations can be constructed. An additional consideration, is that experiment repetition allows for the discrimination between unchanging astrophysical sources and j noise. The rate of experimental repitition will determine the frequencies at which this separation is possible, ie. The less the slow correlated noise changes between experimental repeats, the more of it that can be determined and removed. 6.4 Optimal Measurement Strategy The problem of determining the optimal measurement strategy can be stated as the task of determining the linear combination of absolute power measurements which will result in a maximally diagonal pixel-pixel covariance matrix. Consider the noise covariance matrix formed from the noise timestream of an absolute power experiment, N- 1 = (n-n ) (6.1) T Form differential data from the absolute power measurements by convolving with a measurement function, h: TV" = (^(h®nj (h®nj ^ 1 T (6.2) The new pixel-pixel noise covariance matrix, is C = P N~ P T 1 (6.3) where P is the new pointing matrix, containing the elements of h placed in the locations corresponding to the pixels pointed to in the absolute power measurements. —* The problem is to determine both h and where spatially the measurements that are weighted by the elements of h should be located. This is specified in the pointing matrix, P. The solution depends on N and, in general, seems difficult. 119 Chapter 6. Discussion 6.5 Noise & Correlation Removal The streaking, cross-hatching and large-scale structure present in the final iteration map are indicative of a non-diagonal pixel-pixel noise covariance matrix. The final iteration map is the least-squares optimal solution to (6.4) Px = M To improve upon the pixel correlation properties in the final iteration map, it is necessary to change the type of optimal solution that is desired. One such approach is through the use of priors in Bayes' equation. An example of a prior is a uniform probablility distribution for the pixels which is zero for negative pixel values. It is unclear which priors will provide an appropriate optimal map solution, or even how to determine them. Perhaps a more difficult alternative is to re-measure the measurements with a strategy which yields a more diagonal pixel-pixel noise covariance matrix. Consider m = Ps (6.5) + n where in is the existing vector of double-difference measurements, P is the unmodified pointing matrix, s the sky and ft is the noise. Let N N m m rhj = cijiPi • s + i=i where aji is the entry on the j th ji i a (6-6) n i=i row and the i th column of A, a re-measurement matrix. Pi is the i row of the pointing matrix, and rhj is the j th th linear combination of existing measurements. Let P = J2a Pi ] ji = (AP) j (6.7) i=l and AT, n i=i then M = Ps + fi (6.9) Chapter 6. Discussion 120 Forming the new pixel-pixel covariance matrix C = (6.10) P N- P T 1 Expanding equation 6.10 gives C = PA T T [ANA ] T 1 AP (6.11) Thus, the task is to determine A such that C is diagonal. The first conclusion to draw, is that if A is invertible, C = C. Solving for A looks tricky; however, we know of one solution. Consider the pointing matrix of an absolute power experiment. Every row of P contains a single one, and, for the sake of discussion, let N contain off diagonal terms consistent with slowly varying structure, or j noise. We perform a single-difference experiment by taking differences between the absolute power measurements. This removes the constant offset in the absolute-power measurements, and reduces the slowly varying structure, reducing the off-diagonal terms in N. At the same time, the new P P matrix becomes more band-diagonal. Note that in this example, A is T a timeline differencing matrix, and is singular. If the measurement strategy is very 1 good, then there are many measurements involving a single pixel with many other pixels and the extra diagonalization of N is accomplished with a nearly diagonal P P. T The determination of A in this example, is a restatement of the problem posed in the previous section. Thus, the task is to determine the singular matrix A that will form linear combinations of the existing SHADES data to produce a maximally diagonal C matrix. A is an JV x N 1 m m square matrix. Appendix A. Jordan Form 121 Appendix A Jordan Form For any square matrix, A, there exists a matrix M such that J = M~ AM (A.l) l where the columns of M consist of linearly independent "generalized eigenvectors" whose construction is given, for example, in [21]. The result is that the decomposition: A = MJM~ (A.2) l exists for any square matrix A and the matrix J consists of the eigenvalues of the A along the diagonal and ones and zeros just above the diagonal. The rest of the matrix is zero. Let A equal the characteristic matrix introduced in Section 3.4.1. We are interested in the behaviour of A n = (MJAT )" (A.3) = MJ M~ (A.4) 1 n l with large n. Consider the case where J in the decomposition of A in equation A.2 is / Ao 1 0 0 \ 0 Ao 0 0 0 0 Ai 1 V0 0 0 Ai (A.5) / 122 Appendix A. Jordan Form where An and Ai are the two eigenvalues of A. Now, ( J A \2 0 2A 0 0 \ \2 0 0 0 0 0 = A 0 0 A? 2Ai \ 0 0 0 A? / 3A, 0 2 / 0 \ 0 0 0 A? 3A 0 0 A? 0 Ag 0 V0 J (A.6) (A.7) 2 / By induction the n power of J is, th ( \o n A J = n nxr 1 0 0 0 0 0 0 0 A" nA^ \0 0 0 A" \ (A.8) 1 / For the elements of A to converge, the absolute eigenvalues of A must be less than one which is the same stability criterion developed in Chapter 3. Taylor expanding the function f(x) = (l-x) (A.9) n about x = 1, and neglecting higher order terms, obtain, f(x) = l n 2 (A.10) Thus, for absolute eigenvalues less than, but close to one, the non-zero off-diagonal terms in A.8 decay with increasing n in the same fashion as the diagonal terms. For absolute eigenvalues closer to zero, the decay of the off-diagonal terms with increasing n is slightly slower than the decay of the diagonal terms with increasing n. Bibliography 123 Bibliography [1] E. N. Archibald, T. Jenness, W. S. Holland, I. M. Coulson, N. E. Jessop, J. A. Stevens, E. I. Robson, R. P. J. Tilanus, W. D. Duncan, and J. F. Lightfoot. On the atmospheric limitations of ground-based submillimetre astronomy using array receivers. MNRAS, 336:1-13, October 2002. [2] R. Arendt, D. Fixen, and H. Moseley. Dithering strategiesforefficient selfcalibration of imaging arrays. ApJ, 2000. [3] I. Aretxaga, D. H. Hughes, E. L. Chapin, E. Gaztahaga, J. S. Dunlop, and R. J. Ivison. Breaking the 'redshift deadlock'- II. The redshift distribution for the submillimetre population of galaxies. MNRAS, 342:759-801, July 2003. [4] Andrew W. Blain, Ian Smail, R. J. Ivison, J.-P. Kneib, and David T. Prayer. Submillimeter galaxies, astro-ph, 2002. [5] C. Borys. Phd thesis: A sub-millimetre survey of dust enshrouded galaxies in the hubble deep field region. 2002. [6] C. Borys, S. Chapman, M. Halpern, and D. Scott. MNRAS, 2002. [7] James Dunlop. Shades proposal: The scuba/blast wide-field extragalactic survey. 2002. [8] D. T. Emerson, U. Klein, and C. G. T. Haslam. A multiple beam technique for overcoming atmospheric limitations to single-dish observations of extended radio sources. A&A, 76:92-105, June 1979. [9] G. Hinshaw, C. Barnes, C. L. Bennett, M. R. Greason, M. Halpern, R. S. Hill, N. Jarosik, A. Kogut, M. Limon, S. S. Meyer, N. Odegard, L. Page, D. N. Bibliography 124 Spergel, G. S. Tucker, J. L. Weiland, E. Wollack, and E. L. Wright. FirstYear Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Data Processing Methods and Systematic Error Limits. ApJS, 148:63-95, September 2003. [10] D. H. Hughes, I. Aretxaga, E. L. Chapin, E. Gaztanaga, J. S. Dunlop, M. J. Devlin, M. Halpern, J. Gundersen, J. Klein, C. B. Netterfield, L. Olmi, D. Scott, and G. Tucker. Breaking the 'redshift deadlock'- I. Constraining the star formation history of galaxies with submillimetre photometric redshifts. MNRAS, 335:871-882, October 2002. [11] D. Johnstone, C. D. Wilson, G. Moriarty-Schieven, J. GiannakopoulouCreighton, and E. Gregersen. Large-Area Mapping at 850 Microns. I. Optimum Image Reconstruction from Chop Measurements. ApJS, 131:505-518, December 2000. [12] M. R. Nolta, E. L. Wright, L. Page, C. L. Bennett, M. Halpern, G. Hinshaw, N. Jarosik, A. Kogut, M. Limon, S. S. Meyer, D. N. Spergel, G. S. Tucker, and E. Wollack. First Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Dark Energy Induced Correlation with Radio Sources. ApJ, May 2003. [13] Press, W. and Teukolsky, S. and Vetterling, W. and Flannery, S. Numerical Recipes in C. Cambridge University Press, second edition, 1992, cl988. [14] Ross, S. Introduction to Probability Models. Academic Press, seventh edition, 2000, cl972. [15] M. Rowan-Robinson. The Star Formation History of the Universe: An Infrared Perspective. ApJ, 549:745-758, March 2001. [16] M. Rowan-Robinson. Photometric redshifts in the Hubble Deep Fields: evolution of extinction and the star formation rate. MNRAS, 345:819-833, November 2003. Bibliography 125 [17] S. Scott. Observing the scuba half degree extragalactic survey (shades). July 2003. [18] S. E. Scott, M. J. Fox, J. S. Dunlop, S. Serjeant, J. A. Peacock, R. J. Ivison, S. Oliver, R. G. Mann, A. Lawrence, A. Efstathiou, M. Rowan-Robinson, D. H. Hughes, E. N. Archibald, A. Blain, and M. Longair. The SCUBA 8-mJy survey - I. Submillimetre maps, sources and number counts. MNRAS, 331:817-838, April 2002. [19] S. Serjeant, J. S. Dunlop, R. G. Mann, M. Rowan-Robinson, D. Hughes, A. Efstathiou, A. Blain, M. Fox, R. J. Ivison, T. Jenness, A. Lawrence, M. Longair, S. Oliver, and J. A. Peacock. Submillimetre observations of the Hubble Deep Field and Flanking Fields. MNRAS, 344:887-904, September 2003. [20] D. Sivia. Data Analysis A Bayesian Tutorial. Oxford University Press, 2003, cl996. [21] G. Strang. Linear Algebra And Its Applications. Harcourt Brace Jovanovich College Publishers, third edition, 1988, cl976. [22] T. Webb, S. Eales, S. J. Lilly, D. L. Clements, L. Dunne, W. Gear, K. Adelberger, M. Brodwin, H. Flores, S. Foucaud, O. Le Fevre, H. McCracken, A. Shapley, C. Steidel, and M. Yun. Dusty Star Forming Galaxies at HighRedshift: The Canada-UK Deep Submillimeter Survey. Bulletin of the American Astronomical Society, 34:571—h, December 2001. [23] N. Wright. ApJ, 1996.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Scuba half-degree extragalactic survey : data reduction...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Scuba half-degree extragalactic survey : data reduction by Iterative deconvolution Lepage, Kyle Quentin 2004
pdf
Page Metadata
Item Metadata
Title | Scuba half-degree extragalactic survey : data reduction by Iterative deconvolution |
Creator |
Lepage, Kyle Quentin |
Date Issued | 2004 |
Description | A generalization of an iterative algorithm used to reduce large data sets consisting of observations of the cosmic microwave background is presented. In particular, the algorithm is modified to reduce data sets obtained by arbitrary observations. The generalized algorithm is applied, in simulation, to a data set consisting of complicated observations taken during a Scuba HAlf Degree Extragalactic Survey (SHADES). The algorithm is found to be unstable. The generalized time-order algorithm is reformulated as a difference equation whose solution gives the Nth iteration reduced data map in analytical form. Within the analytical framework, the convergence properties of the algorithm are determined, a stable, known work-around explained, and the converged map shown to be optimal in the least-squares sense. The known work-around, a modification of the generalized iterative algorithm, trades stability for ideal behaviour at the edge of the reduced data map. Using simulations of a sky populated by a realistic source count model, it is shown that the deleterious effect of the modified algorithm towards the edge of the map is negligible as compared to the statistical noise for most of the half-degree square map. Simulations of the modified algorithm applied to data taken with the SHADES measurement strategy with Gaussian noise added, reveal the presence of streaking and large scale fluctuations in the converged map. The streaking which is demonstrated to be affected only by the measurement strategy, is indicative of a band-diagonal pixel-pixel covariance matrix, and is absent in noiseless simulations. The large scale fluctuations are indicative of non-negligible, far off-diagonal terms in the pixel-pixel covariance matrix, and are also absent in noiseless simulations. The effect of different observation strategies on the streaking is investigated by simulation. The presence of the amplified large scale noise is explained by linear systems theory. The determination of the optimal observation strategy is formulated as the task of determining the linear combination of absolute power measurements which will yield a desired pixel-pixel covariance matrix. This problem is discussed. An un-deconvolved map of real SHADES data is produced. Sources are detected in the un-deconvolved map using Bayesian techniques. |
Extent | 14086236 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-11-18 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0091468 |
URI | http://hdl.handle.net/2429/15287 |
Degree |
Master of Science - MSc |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2004-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2004-0253.pdf [ 13.43MB ]
- Metadata
- JSON: 831-1.0091468.json
- JSON-LD: 831-1.0091468-ld.json
- RDF/XML (Pretty): 831-1.0091468-rdf.xml
- RDF/JSON: 831-1.0091468-rdf.json
- Turtle: 831-1.0091468-turtle.txt
- N-Triples: 831-1.0091468-rdf-ntriples.txt
- Original Record: 831-1.0091468-source.json
- Full Text
- 831-1.0091468-fulltext.txt
- Citation
- 831-1.0091468.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0091468/manifest