text:
# First results from the VSA – III. The CMB power spectrum

## 1 Introduction

Anisotropies in the Cosmic Microwave Background (CMB) have now been detected by many experiments (most recently Netterfield et al. (2001), Lee et al. (2001), Halverson et al. (2002), Padin et al. (2001)). At present the most successful model for explaining the origin of these fluctuations postulates that they are seeded in the very early universe by quantum perturbations with random phase, followed by a period of inflationary expansion. The fluctuations in the CMB are predicted to be Gaussian in nature, and hence can be completely characterised through their power spectrum. A further prediction is that the power spectrum will show acoustic peaks due to plasma oscillations on scales smaller than the sound horizon at the surface of last scattering.

In this paper we present the power spectrum of the CMB fluctuations detected by the Very Small Array (VSA) on spherical harmonic modes $`\mathrm{}`$150–900. We outline how the fully calibrated time-stream data are converted into a power spectrum and the various data checks that we have performed to confirm the validity of our analysis. This paper is the third in a series of four papers which report the results of the first season of observations made using the VSA in its compact configuration. Paper I (Watson et al., 2002) describes the design of the VSA and our experimental method; the observational strategy, foreground removal and reduction techniques for the data analysed in this paper are described in Paper II (Taylor et al., 2002); finally, the cosmological implications of the VSA power spectrum are discussed in Rubiño-Martin et al. (2002) (Paper IV).

## 2 Observations and Initial Data Processing

### 2.1 The observations

The Very Small Array (VSA) is a 14-element interferometer array for Cosmic Microwave Background (CMB) observations; it operates at a frequency between 26 and 36 GHz with a receiver bandwidth of 1.5 GHz. In its compact configuration, used here, the instrument is most sensitive to angular structure in the range $`\mathrm{}`$150–900. As well as the 14-element array, there is also a single baseline interferometer used for radio source flux measurements. A fuller description of the instrument is given in Paper I.

The present observations were made at a frequency of 34 GHz in the period 2000 September to 2001 September and were centred on three separate areas of sky. Each VSA observation maps a region with FWHM of $`4.6^{}`$. The process of field selection is discussed in Paper II; selection was based primarily on low Galactic (synchrotron, free-free and dust) emission and an absence of known bright foreground sources. Overlapping fields were observed in each area in order to reduce the sample variance, to increase the resolution in $`\mathrm{}`$-space and to allow for direct assessment of data reliability and the detection of any residual instrumental effects (Section 4.2). The array configuration used was designed to provide an approximately uniform spread of interferometer baselines while retaining a reasonable aperture filling factor.

The amplitude and phase calibrations of the individual interferometer baselines were checked both by frequent short measurements of a number of secondary calibration sources and also by regular longer observations of a small number of primary calibration sources. The overall calibration procedure is described fully in Paper I; the overall accuracy of the calibration in flux density and temperature is $`3.5\%`$ .

An important aspect of the VSA is the inclusion of a separate interferometer, comprising two 3.7-m dishes on a 9-m north-south baseline, for determining the flux densities of foreground sources (radiogalaxies and quasars) falling within the observed fields. The positions of all sources which might affect the VSA observations were obtained from survey observations using the Ryle Telescope (Waldram et al., 2002); these positions were subsequently observed concurrently with the main VSA observations in a series of regular pointed observations and the measured flux densities subtracted from the visibilities observed by the main array (Paper II).

### 2.2 Initial Data Processing

The data are calibrated and processed as described in Papers I and II. Early tests of the telescope revealed an unwanted local spurious signal, particularly evident on the shorter baselines. The procedures for removing this signal from the data are described fully in Paper II; the tests presented in Section 4.2 below indicate that the removal of this signal is complete. Similar filtering processes were used to eliminate signals arising from the Sun and the Moon.

## 3 Derivation of the Power Spectrum

Since the three sets of mosaiced fields are widely separated from one another on the sky, the data corresponding to each field are analysed individually. The derivation of the CMB power spectrum for each mosaic is performed using the maximum-likelihood method presented in Hobson & Maisinger (2002), which we summarise here. The results from each mosaic are then combined, as outlined below, to produce the final estimate of the power spectrum.

Since the number of individual visibility measurements for each mosaic is very large ($``$ 800,000), it is first necessary to compress these data in some way. For each separate field in the mosaic, the visibilities are binned into cells in the $`uv`$-plane; for the case in which the instrumental noise covariance matrix is diagonal, this corresponds to the maximum-likelihood solution for the value of the binned signal visibility in each cell. This is analogous to the ‘map-making’ step in the analysis of single-dish CMB observations, in which time-ordered data are binned into pixels on the sky (see, e.g. Borrill (1999)). Since we are not interested here in making accurate CMB maps from the binned visibility data, the $`uv`$-plane is simply divided into equal-area cells of size $`\mathrm{\Delta }u=3`$ wavelengths. As the aperture function of the compact VSA is well-modelled by a Gaussian with a FWHM of 12 wavelengths, this cell size ensures that the $`uv`$-plane is comfortably over-sampled, while reducing the total number of visibility measurements for each field significantly (to $``$ 2500).

The binned visibilities are thus the basic input to the likelihood analysis for the CMB power spectrum. The compact VSA is sensitive to the multipole range $`\mathrm{}150900`$, and the effective aperture function (after mosaicing) has a FWHM of $`\mathrm{\Delta }\mathrm{}=83`$. We thus divide the total $`\mathrm{}`$-range into 10 spectral bands each of width $`\mathrm{\Delta }\mathrm{}`$, in order to limit the correlation between adjacent bands. In each band, we assume $`\mathrm{}(\mathrm{}+1)C_{\mathrm{}}`$ to have a constant value $`\overline{𝒞}_B`$ $`(B=1,2,\mathrm{},10)`$. These flat power levels in each band are then the parameters to be determined in our likelihood analysis of the data for each mosaic. We denote these band powers collectively by the parameter vector $`𝒂`$.

Assuming the CMB emission and the instrumental noise to be Gaussian random fields, the log-likelihood of obtaining the binned visibility data vector $`𝒗`$ (which contains the real and imaginary parts separately of each binned visibility), given some set of flat band powers $`𝒂`$, is given by

$$\mathrm{ln}=\text{constant}\frac{1}{2}\left[\mathrm{ln}|𝑪(𝒂)|+𝒗^\mathrm{t}𝑪^1(𝒂)𝒗\right],$$

where $`𝑪(𝒂)=𝑺(𝒂)+𝑵`$ is the sum of the predicted signal covariance matrix and the noise covariance matrix. The maximum-likelihood CMB power spectrum $`\widehat{𝒂}`$ is calculated using a simple iterative numerical maximisation algorithm. Starting from an initial guess $`𝒂_0`$ (which is unimportant), independent line maximisations are performed for each band power $`𝒞_B`$ in turn, while keeping the others fixed. The whole solution vector $`𝒂`$ is then updated and the process repeated until convergence is obtained. This typically requires around 5 iterations.

The well-defined correlation structure of visibility data in the $`uv`$-plane allows each line-maximisation to be performed using only the subset of visibilities that are sensitive to the band power being varied, thereby speeding-up the evaluation of the likelihood function, while keeping the calculation exact. For a single VSA mosaic, the maximum-likelihood solution can be obtained in around one hour on an 8-node Beowulf Cluster with 1.8-GHz AMD Athlon processors.

The uncertainties in the derived maximum-likelihood CMB power spectrum are estimated in two complementary ways. Assuming the likelihood function in the parameter space $`𝒂`$ to be well-approximated by a multivariate Gaussian near its peak, the covariance matrix of the parameter uncertainties is given simply by (minus) the inverse of the curvature (or Hessian) matrix at the peak $`\widehat{𝒂}`$. This matrix is easily evaluated numerically in a few hours of CPU time. The square-root of the diagonal elements of the resulting covariance matrix give the standard error $`\mathrm{\Delta }\overline{𝒞}_\mathrm{B}`$ on each band power, whereas the off-diagonal elements provide a measure of the correlation between the band power estimates in different spectral bins. We find, typically, that adjacent bins are anti-correlated at around the $`520\%`$ level, and for more widely separated bins the correlation is negligible.

Using the covariance matrix clearly produces symmetric error bars on each band power $`\overline{𝒞}_\mathrm{B}`$, which may be a poor representation of the uncertainty, especially for poorly constrained band powers. An alternative approach is to make use of the fact that the band power estimates in different spectral bins are only weakly correlated, so that the off-diagonal elements of the covariance matrix are small compared to those on the diagonal. In this case, a better representation of the uncertainty in each estimate may be obtained by directly evaluating the likelihood function through the peak $`\widehat{𝒂}`$, along each parameter direction in turn. In the ideal case, where the band power estimates are independent, the resulting curves would be the marginal distributions of each band power. Since the likelihood function can be evaluated very quickly along each direction in parameter space, the resulting ‘marginal’ distributions can be calculated for each VSA mosaic in around 2 hours, and provide a useful cross-check of the standard errors obtained from the covariance matrix at the peak. We find that for spectral bins in which the flat band power is tightly-constrained, the shape of the likelihood function is very close to Gaussian. However, for bins in which the power level less well constrained, the likelihood function is better described by an offset log-normal distribution; this is discussed further in Paper IV.

Once the likelihood functions for the flat band power $`\overline{𝒞}_\mathrm{B}`$ in each spectral bin have been obtained separately for each VSA mosaic, the mosaics are combined simply by multiplying together the respective likelihood functions in each bin. This assumes that each mosaic provides an independent measurement of the CMB power spectrum in each spectral bin, which is valid given that the three sets of mosaiced fields are widely separated on the sky. The correlation between the resulting band power estimates in different bins is easily obtained by calculating the covariance matrix, as described above, at the new joint optimum $`\widehat{𝒂}`$.

All the above functions are implemented using the Madcow analysis package (Hobson & Maisinger, 2002).

### 3.1 Window functions

Because of various instrumental effects (e.g. non-uniform $`uv`$ coverage and the finite size of the primary beam of the telescope), a bin samples the underlying power spectrum $`C_l`$ through a window function $`W(l)`$. For the Bth bin, the measured power corresponds to

$$p_B=\underset{\mathrm{}}{}\frac{W_B(\mathrm{})}{\mathrm{}}𝒞_{\mathrm{}}$$

where $`𝒞_{\mathrm{}}=\mathrm{}(\mathrm{}+1)C_{\mathrm{}}/(2\pi )`$.

The window function for a given bin is determined as follows. The range of Fourier modes to which a given visibility is sensitive are given by the Fourier transform of the primary beam for a single field observation (see e.g. Hobson & Maisinger (2002)). For a Gaussian primary beam,

$$\stackrel{~}{A}(\text{u})=2\pi \sigma ^2\mathrm{exp}(2\pi ^2\sigma ^2|\text{u}|^2),$$

where u is the vector in the visibility plane (measured in wavelengths). For mosaiced observations, the effective beam is a superposition of displaced primary beams. We can think of this superposition as a convolution of a centred primary beam with a sum of delta functions at the beam centres. In the Fourier domain this corresponds to a multiplication of the aperture function with the ‘Fraunhofer diffraction pattern’ of the beam centres

$$\stackrel{~}{A}_{\text{eff}}(\text{u})=\mathrm{\Sigma }_j\stackrel{~}{A}(\text{u})e^{2\pi i\text{u}\text{x}_j},$$

where $`\text{x}_j`$ are the positions of the beam centres from the origin, measured in radians. Note that in general this aperture function is complex.

For a given bin, the weighted complex sensitivity map is

$$S(\text{u})=\underset{k}{}w_k\stackrel{~}{A}_{\text{eff}}(\text{u}\text{v}_k),$$

where index $`k`$ runs over visibilities in a given bin. These visibilities have instrumental weights $`w_k`$ and positions on the visibility plane $`\text{v}_k`$.

The un-normalised window function is then given by:

$$W(\mathrm{})=𝒩_{\varphi =0}^{\varphi =2\pi }|S(\mathrm{},\varphi )|^2d\varphi ,$$

where $`\mathrm{}=2\pi |\text{u}|`$. The normalisation constant $`𝒩`$ can be trivially found using $`_{\mathrm{}}(2\mathrm{}+1)W(\mathrm{})/(2\mathrm{}(\mathrm{}+1))=1`$. The window functions for the combined VSA data set are plotted in Figure 1, showing the small degree of correlation between adjacent bins. These functions are used in Paper IV in the estimation of cosmological parameters.

## 4 Data Checks

The complete process of editing and filtering the data and also the subsequent stages of data reduction were carried out independently by the Cambridge group and, jointly, by the combined IAC and JBO teams. A comparison of the two sets of results showed good agreement, the effect of any differences being small compared to the intrinsic uncertainties on the final power spectra.

### 4.1 Test of data reduction procedure

Aspects of our data reduction procedure, such as filtering and calibration, could potentially have introduced systematic errors into the VSA data. In order to test this, we produced a realisation of the CMB sky and simulated a mock VSA observation including such instrumental effects as visibility quadrature errors, phase steps due to path compensation and thermal noise. We analysed these simulated data using our standard reduction procedure and then produced a power spectrum as described in section 3. We found this to be entirely consistent with our input model and so concluded that our method of data reduction did not introduce any significant systematic errors. In particular, it demonstrates that any correlations in the visibilities introduced by the filtering process has not observable effect on the derived power spectrum.

### 4.2 The Local Spurious Signal

The local spurious signal, fully discussed in Paper I, was found to depend only on the antenna tracking angle, and not on the table elevation. Therefore, we would expect the spurious signal to be identical for different fields with the same declination observed over the same hour angle range. This is easily confirmed by combining the (unfiltered) baseline timeseries of two such fields by both addition and subtraction. We find that adding the two fields enhances the spurious signal, whilst subtraction entirely removes it. Whilst this technique of addition and subtraction is adequate for detecting the presence of spurious signal in unfiltered data, it is insufficiently sensitive to detect possible low level residual signal once filtering has been applied.

The increase in detection sensitivity that we require in order to test the filtered data for residual spurious signal can, however, be obtained using a modified version of a MEM algorithm used for extraction of CMB signal from VSA data (Maisinger et al., 1997).

We add an extra term to the MEM reconstruction which is the signal that is identical in the CMB datasets. We then consider the case of two fields at identical declination. As the CMB signal in the two fields will not be identical, and the noise is random in each case separately, any common component to the two fields will be spurious signal.

We tested this algorithm by applying it to pairs of simulated observations of CMB fields assuming CDM primordial fluctuations, to which we added an identical component with an rms level such that it was not the dominant signal. We used a variety of common components, including a scaled down version of the unfiltered spurious signal. We found that the MEM algorithm was able to reconstruct the common signals well, recovering the structure excellently, and recovering the amplitude of the signal to within 10%. As we are primarily interested in whether or not the spurious signal is still present, as opposed to any accurate quantification, this is perfectly adequate. Note that, even for entirely independent fields, the reconstructed shared signal power spectra are not zero (Figure 2), but rather show a value increasing as $`\mathrm{}^2`$, consistent with the correlation between two white noise signals.

The two VSA2 fields had identical declination to allow testing of our filtering procedures. The MEM algorithm was applied to these fields, and the results compared with CDM realisations with identical uv-coverage and thermal noise in which no shared signal is present. The results (Figure 2) show that the common signal found from the real data fields and the simulations agree within the errors, the difference between the two sets of points being less than 4% of the measured CMB for values of $`\mathrm{}<750`$. In the same manner, we find no evidence for residual spurious signal in the pairs of VSA1,3 fields at similar declinations.

Non-Gaussianity analysis of the binned visibilities allows us to locate and remove the few remaining visibilities contaminated with spurious signals down to a low level. The removal of these points has a negligible effect on the final power spectrum, giving us confidence that we are subtracting the spurious signal to a level well below that which could affect our results. The full details of this analysis will be published in Savage et al. in prep.

### 4.3 Foreground Source Subtraction

Radio galaxies and quasars are a significant contaminant of the CMB at microwave frequencies and in the higher $`\mathrm{}`$-ranges will dominate the CMB signal, making it essential to remove their contribution. Tests with simulated fields have been used to assess the potential contribution, before subtraction, of these sources to the final CMB spectrum. We generated ten realisations of the CMB sky using a particular CDM model and added to these the sources that we have observed to be present in the VSA1 field. We have then compared the power spectra recovered from simulated VSA observations with and without sources. Two of these simulations, representing the range of results obtained from the ten simulations, are shown in Figure 3. It is apparent that, with the fairly small number of sources in any one VSA field (typically 12 sources), the impact on the power spectrum is unpredictable. Although the main contribution tends to be in the higher $`\mathrm{}`$-bins (with errors potentially reaching $`\pm `$100% at $`\mathrm{}=900`$), changes of up to $`\pm `$10% can occur in the lowest $`\mathrm{}`$-bin; these are due to the chance superposition of interference fringes.

To determine the effect of residual sources on the VSA results, two further simulations without a CMB contribution have been carried out. We base these simulations on the 15-GHz source counts from Taylor et al. (2001) and extrapolate up to 34 GHz using a mean spectral index of $`\alpha =0.55`$. The first (Figure 4 upper plot) comprises the contributions of the known point sources and a statistical distribution of weaker sources; the second (lower plot) includes only the statistical contribution of weaker sources and gives an indication of the possible residual contribution to the VSA power spectrum after subtraction of known sources. It is clear that, with no source subtraction, the CMB data are significantly compromised for $`\mathrm{}`$-values $`600`$. After subtraction of known sources, the contribution is reasonably small for $`\mathrm{}`$-values up to about 1000. As demonstrated by the results of Figure 3, the contribution of sources to the observed spectra can not generally be predicted from a simple combination of simulated power spectra, although such an approach becomes feasible in the limiting case of many weak sources per synthesised beam.

We can also estimate the residual source contribution to the power spectrum using the preliminary 34 GHz source count derived from the source subtractor observations in Paper II. Integrating the count from zero flux to our complete source subtraction limit of 80 mJy and converting to units of $`\mathrm{\Delta }T/T_0`$, the source power is given by $`C_{\mathrm{}}=7.7\times 10^{16}`$, corresponding to a power spectrum value of $`\frac{\mathrm{}(\mathrm{}+1)}{2\pi }C_{\mathrm{}}=580\mu `$K<sup>2</sup> at $`\mathrm{}=800`$. For the case with effectively no source subtraction, taking the upper flux limit to be $`0.5`$ Jy, we find $`C_{\mathrm{}}=3.7\times 10^{15}`$, equivalent to $`2800\mu `$K<sup>2</sup> at $`\mathrm{}=1000`$; these results are in good agreement with the extrapolations from the 15 GHz counts.

The precise residual contribution due to sources is subject to uncertainties in our knowledge of the true weak-source distribution; the values shown in Figure 4 are likely to be over-estimates since some sources with flux densities less than 80 mJy were actually subtracted. In order to assess the possible impact of these residual sources on derived parameters (Paper 4) the cosmic parameter analysis has been repeated after subtracting the residual contributions shown Figure 4 from the observed power spectrum. The changes in the fitted parameters are small, the largest effect being a reduction of 0.05 in $`n_s`$, representing a change of about 0.5 s.d.

### 4.4 Galactic Foregrounds

The diffuse Galactic foregrounds are discussed in Paper II, where it is suggested that the total contribution from the three components of the Galactic foreground (synchrotron, free-free and possibly spinning dust) amount to no more than about $`\mathrm{\Delta }T=5\mu `$K at an angular scale of $`1^{}`$ ($`\mathrm{}200`$). Here we discuss their effect on the results presented in this paper. The free-free and synchrotron components are relatively well-known and are each expected to contribute about $`\mathrm{\Delta }T=12\mu `$K. The spinning dust component is more uncertain, but may be the dominant component, perhaps contributing up to $`\mathrm{\Delta }T5\mu `$K. This is based on a dust-correlated component with a correlation coefficient of $`10\mu `$K/(MJy sr<sup>-1</sup>) at $`100\mu `$m. This coefficient is still to be clearly demonstrated observationally. The power spectra of the diffuse foregrounds falls with increasing $`\mathrm{}`$ (for example, Giardino et al. (2001)) and is likely to be half these values at $`\mathrm{}500`$ .

Any Galactic contribution adds in quadrature with the CMB signal and hence at the position of the first CMB peak ($`\mathrm{}200`$) which has $`\mathrm{\Delta }T_{\mathrm{CMB}}75\mu `$K, $`5\mu `$K of foreground signal will increase the observed signal by $`0.17\mu `$K. Similarly, at $`\mathrm{}500`$ where $`\mathrm{\Delta }T_{\mathrm{CMB}}45\mu `$K, the increase will be $`0.28\mu `$K at the most. We see therefore, for the VSA fields, and at a frequency of 34 GHz, the contibution from Galactic foregrounds is likely to be negligible ($`<1`$ percent). A similar conclusion, albeit based on observations of different regions of sky, was reached by Halverson et al. (2002) who found that Galactic emission made a negligible contribution to the observed CMB power spectrum. A more complete cross-correlation analysis to investigate the contribution from dust-correlated emission is in progress (Dickinson et al. in prep.).

### 4.5 Noise Estimation

The likelihood analysis used for power-spectrum estimation requires an accurate estimate of the rms noise level. Since on individual baselines the contribution of the CMB to the individual data samples is very small, the noise level can be obtained directly from the standard deviation of the data, after the filtering and flagging processes have been completed. The noise level associated with each of the binned visibilities (Section 3) is obtained from an appropriately weighted combination of the noise levels of the data in each bin. As an additional check, the scatter in the data points contributing to each bin has also been used to provide a noise estimate. Consistent estimates of the noise level were obtained by the two methods. The overall noise estimate is accurate to 2.5%.

The sensitivity of the likelihood analysis to errors in the noise level has been tested by analysing the same dataset with different assumed levels of noise; the corresponding changes in estimated CMB power produced by the above uncertainty in estimated noise level are less than 1% for $`\mathrm{}`$-values up to $`750`$.

### 4.6 $`\chi ^2`$ Tests

As an additional check on the consistency of the data, we computed the $`\chi ^2`$ statistic for a variety of splits on the visibility data. The data from each VSA field were split into two and the difference vector, $`𝚫=(𝒗_1𝒗_2)`$, formed where $`𝒗_1`$ and $`𝒗_2`$ are the visibility vectors for each half of the data. The $`\chi ^2`$ statistic was then calculated using $`\chi ^2=𝚫^\mathrm{t}𝑪_\mathrm{N}^1𝚫`$, where $`𝑪_\mathrm{N}`$ is the noise covariance matrix given by $`𝑪_\mathrm{N}=(𝑪_{\mathrm{N1}}+𝑪_{\mathrm{N2}})`$. Since for each field the sky signal measured in each of the two halves of the data should be the same, the $`\chi ^2`$ statistic can be used to test for the presence of systematic errors in the data. The data for each VSA field were split in two according to observing epoch and the $`\chi ^2`$ values and associated significances are given in Table1. The data split on the VSA3 field also corresponds to a split between day and night observations. This was the only field to be observed in this way, and its $`\chi ^2`$-value confirms the consistency of daytime versus nightime observation.

The consistency of the power spectra derived from each of the 3 VSA mosaiced fields (VSA1M,VSA2M and VSA3M) was also compared by forming the $`\chi ^2`$ statistic on pairs of power spectra. In this case the $`\chi ^2`$ value is given by

$$\chi ^2=(𝒂_1𝒂_2)^t(𝑯_1^1+𝑯_2^1)^1(𝒂_1𝒂_2),$$

where $`𝒂_1`$ and $`𝒂_2`$ are the two sets of bandpowers and $`𝑯_1^1`$ and $`𝑯_2^1`$ are the corresponding inverse Hessian matrices. Since $`\chi ^2`$ statistic assumes that the likelihoods are Gaussian, we use the Hessians calculated using an offset log-normal approximation (see Paper IV for further details). The $`\chi ^2`$ values (and significances) for the VSA1M/VSA2M, VSA1M/VSA3M and VSA2M/VSA3M power spectra comparisons are 7.6 (0.67), 9.82 (0.46) and 5.03 (0.89) respectively. In each case there are ten degrees of freedom in the power spectrum analysis.

## 5 The Power Spectrum

The filtered and source-subtracted data for each of the VSA fields have been analysed using the Madcow software package (Hobson & Maisinger, 2002) as described in Section 3. The combination of the three mosaiced spectra is shown in Figure 5. The bin width used is $`\mathrm{\Delta }\mathrm{}=83`$, which gives weakly-correlated errors in each bin. The actual correlations between bins are given by the correlation matrix (Table 3). To reduce the bias in assessing features in the power spectrum caused by the settings of the bin centres, we have also calculated the power spectrum with bin centres shifted by one half a bin width to the right of the original bin centres. These results are shown in Figure 5 with dashed error bars. Adjacent ‘double-binned’ points are highly correlated but do sensibly sample the power spectrum of our data. Numerical values for both binnings are given in Table 2.

The plotted error bars contain the contributions from both thermal noise and sample variance, but not calibration errors, which introduce a completely correlated uncertainty in all the points of $`\pm 7`$ percent. Errors from pointing and primary beam uncertainties are negligible. Since the temperature sensitivity of the VSA compact configuration falls off dramatically after $`\mathrm{}=800`$, all data above this have been binned together.

## 6 Discussion

### 6.1 The VSA power spectrum

The power spectrum shown in Figure 5 shows a clear detection of the first peak at $`\mathrm{}220`$, and power at the level of about $`2000\mu `$K<sup>2</sup> between $`\mathrm{}=300`$ and $`\mathrm{}=900`$. We have attempted to quantify the detection of a second peak at $`\mathrm{}550`$, as this is the region of the power spectrum with the largest anti-correlations between adjacent bins (see Table 3); bin 6 centred at $`\mathrm{}=556`$ is anti-correlated with its neighbours at the $`20`$ percent level. We made Monte Carlo simulations of the five $`C_{\mathrm{}}`$ points between $`\mathrm{}=390`$ and $`\mathrm{}=722`$. Points were draw from the 5-dimensional Gaussian distribution described by the actual correlation matrix of these points, but with mean values equal to the weighted mean of the five actual points. Following Hobson & Magueijo (1996), who define a normalised convexity about a power spectrum point, we calculated the change in normalised convexity of the power spectrum about the inner 3 points, defined by

$$\mathrm{\Sigma }=\frac{𝒞}{\sigma },$$

where $`𝒞=(C_4+C_83(C_5+C_7))/2+2C_6`$ and $`\sigma `$ is the overall error in $`𝒞`$ given from the errors in the individual points $`\sigma _i`$ by $`\sigma ^2=(\sigma _4^2+\sigma _8^2+9(\sigma _5^2+\sigma _7^2))/4+4\sigma _6^2`$. This effectively compares the hypothesis that the power spectrum in this region is flat, to the one that it is described by ‘trough–peak–trough’.

In 1000 realisations, we found only 27 instances of $`\mathrm{\Sigma }`$ being larger than the value observed in the real data of $`\mathrm{\Sigma }=2.2`$. We therefore conclude that the observed second peak is detected at $`97`$ percent confidence. The power spectrum is completely consistent with the adiabatic inflationary models, fits to which are discussed in Paper IV.

### 6.2 Comparison with other experiments

In Figure 6 we compare the new VSA power spectrum plotted with those from BOOMERANG (Netterfield et al., 2001), DASI (Halverson et al., 2002), MAXIMA (Lee et al., 2001). Only single-binned (weakly-correlated) points are shown. We have attempted to compare the random and correlated errors on the various experiments in a consistent way, difficult though this is on a single plot. Two sets of error bars are shown for each plot; the smaller bars indicate 68 percent confidence limits from the random (thermal and sample variance) errors, while the larger error bars represent systematic (calibration and beam) errors as reported – all the points from a single experiment are able to move up or down within the larger error bars. It is significant that, although the overall errors of the different experiments are comparable, the relative contribution of systematic errors (including beam uncertainty) is much smaller for the interferometric data.

The agreement between the experiments on the existence, heights and positions of two peaks and of power at higher $`\mathrm{}`$, is evident. This is particularly significant given the very different experimental techniques involved and the different foregrounds and systematic errors faced by the different experiments. The points in Figure 6 have been obtained over a frequency range of $`26`$$`150`$ GHz, and by ground-based interferometers and balloon-borne scanned total-power telescopes. They are all from different regions of the sky, the calibrations are all independent and based on different absolute calibration sources, and for the two low-frequency experiments foreground sources have been subtracted in different ways, yet the agreement of the power spectra is striking.

A detailed comparison between the experiments is difficult to do from the data points alone because of the correlated errors between points for each experiment. To make a meaningful comparison it is necessary to fit the underlying power spectrum to each data set, taking into account the correlations, and to compare the parameters that describe the power spectrum. In principle many parameterisations of the power spectrum would suffice for this comparison, but in practice it is obviously sensible to use the standard adiabatic cold dark matter power spectrum models, and fits to these models for the VSA and other experiments are considered in detail in Paper IV.

## 7 Conclusions

We have derived the power spectrum of the CMB anisotropies from the first year’s VSA observations, made using its compact array configuration. We measure the flat band power in 10 weakly-correlated bins of width $`\mathrm{\Delta }\mathrm{}=83`$ between $`\mathrm{}150`$ and $`\mathrm{}900`$. The results are subject to a calibration uncertainty of $`\pm 7`$ percent in power, with negligible beam uncertainty. The contribution to the power spectrum from diffuse Galactic emission and residual radio sources is also negligible. Our results are in excellent agreement with other recent measurements as regards the amplitude and position of two peaks in the power spectrum; power is also detected out to the resolution limit of the experiment at $`\mathrm{}900`$.

Band powers, correlation matrices and window functions are available from

`http://mrao.cam.ac.uk/telescopes/vsa/results.html`.

## ACKNOWLEDGEMENTS

We thank the staff of the Mullard Radio Astronomy Observatory, Jodrell Bank Observatory and the Teide Observatory for invaluable assistance in the commissioning and operation of the VSA. The VSA is supported by PPARC and the IAC. Partial financial support was provided by Spanish Ministry of Science and Technology project AYA2001-1657. A. Taylor, R. Savage, B. Rusholme, C. Dickinson acknowledge support by PPARC studentships. K. Cleary and J. A. Rubiño-Martin acknowledge Marie Curie Fellowships of the European Community programme EARASTARGAL, “The Evolution of Stars and Galaxies”, under contract HPMT-CT-2000-00132. K. Maisinger acknowledges support from an EU Marie Curie Fellowship. A. Slosar acknowledges the support of St. Johns College, Cambridge. We thank Professor Jasper Wall for assistance and advice throughout the project.


format: text

text:
# A Primal-Dual Parallel Approximation Technique Applied to Weighted Set and Vertex Cover

## 1 Introduction

The linear-programming primal-dual method for obtaining sequential algorithms for exact optimization problems is well studied \[Ch83\]. Primal-dual techniques have also been used to obtain sequential approximation algorithms (e.g., for NP-hard problems \[Ch79, Ho82, etc.\] and for on-line problems \[You\]). In this paper, we apply primal-dual techniques to obtain a deterministic parallel approximation algorithm for the minimum-weight vertex- and set-cover problems and their duals, maximum-weight edge and element packing.

The result demonstrates that linear-programming primal-dual techniques can lead to fast, efficient parallel algorithms. The algorithm is natural, yet fits no existing paradigm for such algorithms, being unique in that it uses only “local” information at each step, yet is deterministic. Generally, such algorithms have required randomization (e.g., \[II86, ABI86, Lu86\]).

Given an $`n`$-vertex, $`m`$-edge graph $`G=(V,E)`$ with vertex weights and an $`ϵ>0`$, our algorithm returns a vertex cover of weight at most $`2/(1ϵ)`$ times the minimum. It uses $`O(\mathrm{ln}^2m\mathrm{ln}\frac{1}{ϵ})`$ time and $`m/\mathrm{ln}^2m`$ processors (i.e., $`O(m\mathrm{ln}\frac{1}{ϵ})`$ operations) on an EREW-PRAM. More generally, the algorithm finds a set cover of weight at most $`r/(1ϵ)`$ times the minimum, using $`O(r\mathrm{ln}^2m\mathrm{ln}\frac{1}{ϵ})`$ time and $`M/\mathrm{ln}^2m`$ processors (i.e., $`O(rM\mathrm{ln}\frac{1}{ϵ})`$ operations). Here $`m`$ is the number of elements, $`M`$ is the sum of the set sizes, and $`r`$ is the maximum number of sets in which any element occurs. (For vertex cover, $`r=2`$.) In each case, the algorithm also implicitly finds a near-maximal dual solution (an edge or element packing) that is also within the corresponding factor of optimal.

The algorithm can be implemented using only integer arithmetic (see §4.2). If the weights are integers and $`1/ϵ`$ is less than the sum of the vertex (resp. set) weights, then the weight of the cover is at most 2 (resp. $`r`$) times the minimum.

### 1.1 Related Work

The first $`r`$-approximation algorithm for weighted vertex/set cover was due to Hochbaum \[Ho82\]. She considered the relaxation of the natural integer linear program for the problem. The dual of this program is maximum edge packing. The so-called complimentary-slackness conditions are that a (fractional) cover and a packing are optimal provided (i) every vertex in the cover has its constraint met in the packing and (ii) every edge with non-zero packing weight has exactly one vertex in the cover. Hochbaum observed that an optimal packing was necessarily maximal, that for any maximal packing, the vertex set formed by the vertices whose packing constraints are met with equality form a cover, that such a packing and cover satisfy (i), and that (i) is sufficient to guarantee $`r`$-approximation because (ii) is approximately satisfied in that every edge has at most $`r`$ vertices in the cover. Since an optimal dual solution can be found in polynomial time by solving the linear program, Hochbaum obtained a polynomial-time algorithm. Bar-Yehuda and Even \[BE81\] observed that sequentially raising the edge-packing weights as much as possible yields a maximal edge packing, thus obtaining a linear-time algorithm. For our algorithm, we relax (i) further, insisting only that every vertex in the cover nearly have its constraint met, and we show how to simultaneously raise many edge-packing weights so that the packing quickly becomes nearly maximal and the weight of the cover formed by the vertices that nearly have their packing constraints met with equality is within $`r/(1ϵ)`$ of optimal.

In \[Cl83\], Clarkson showed that in a restricted class of graphs, approximation ratios better than 2 could be obtained for vertex cover. Clarkson gave the first parallel approximation algorithm — a relatively complicated randomized algorithm \[Cl91\]. According to Motwani’s lecture notes on approximation algorithms \[Mot92\], which contain a survey of results on vertex cover, the best approximation ratio known is $`2\frac{\mathrm{log}\mathrm{log}n}{2\mathrm{log}n}`$, due to Bar-Yehuda and Even \[BE85\] and to Monien and Speckenmeyer \[MS85\]. In \[Ho83\], Hochbaum gives a $`(22/k)`$-approximation algorithm, where $`k`$ is the maximum vertex degree, and she conjectures that there is no polynomial-time $`c`$-approximation algorithm for any $`c<2`$ unless P$`=`$NP.

Chvátal’s weighted-set-cover algorithm guarantees a set cover of weight at most $`\mathrm{ln}\mathrm{\Delta }`$ times the minimum, where $`\mathrm{\Delta }`$ is the maximum set size \[Ch79, Lo75, Jo74\]. Berger, Rompel, and Shor \[BRS89\] give a parallel algorithm that guarantees a factor of $`(1+ϵ)\mathrm{ln}\mathrm{\Delta }`$. Their algorithm uses a linear number of processors and runs in polylogarithmic time with some restrictions on the weights.

The intuition behind our complexity analysis relies on a lemma of general interest for parallel graph algorithms (Lemma 6). The lemma has previously found application in the analyses of randomized parallel graph algorithms: Israeli and Itai’s maximum-matching algorithm \[II86\] and Alon, Babai and Itai’s maximal-independent-set algorithm \[ABI86\].

In concurrent independent work, Cohen gives a parallel approximation algorithm for maximum flow in shallow networks \[Co92\]. If network flows are viewed as packings of source-to-sink paths, then maximal packings correspond to blocking flows. Cohen gives an $`ϵ`$-blocking flow algorithm that is similar in spirit to our algorithm, although a number of different issues arise. In a more recent work, Luby and Nisan give a parallel primal-dual approximation algorithm for positive linear programming \[LN93\]. Hochbaum’s original algorithm can be parallelized by employing Luby and Nisan’s algorithm; the resulting algorithm would obtain an approximation ratio comparable to ours and have an incomparable running time (growing linearly with $`1/ϵ`$, but not with $`r`$). Previously, Goldberg et al. \[GPST92\] gave a parallel primal-dual algorithm to find (exactly) maximum-weight bipartite matchings. Their algorithm appears to be the first parallel algorithm to use primal-dual techniques, but it requires polynomial time.

### 1.2 Problem Definitions

Let $`G=(V,E2^V)`$ be a given hypergraph with vertex weights $`w:V\mathrm{}^+`$. Let $`E(v)`$ denote the set of edges incident to vertex $`v`$. Let $`G`$ have $`m`$ edges. Let $`r`$, the rank of $`G`$, be the maximum size of any edge. (For an ordinary graph, $`r=2`$.) Let $`M`$, the size of $`G`$, be the sum of the edge sizes. For any real-valued function $`f`$ and a subset $`S`$ of its domain, let $`f(S)`$ denote $`_{xS}f(x)`$.

#### Vertex Cover.

A vertex cover for $`G`$ is a subset $`CV`$ of the vertices such that for each edge $`eE`$, some vertex in $`e`$ is in $`C`$. The (minimum-weight) vertex-cover problem is to find a vertex cover with minimum total weight $`w(C)`$.

#### Edge Packing.

An edge packing is an assignment $`p:E\mathrm{}^+`$ of non-negative weights to the edges of the hypergraph such that the total weight $`p(E(v))`$ assigned to the edges incident to any vertex $`v`$ is at most $`w(v)`$. The (maximum-weight) edge-packing problem is to find an edge packing maximizing $`p(E)`$, the weight of $`p`$. The fractional relaxations of the vertex cover and edge packing problems are linear programming duals.

### 1.3 Related Problems

Let $`𝒞`$ be a family of sets with weights $`w:𝒞\mathrm{}^+`$. Let $`U`$ denote $`_{S𝒞}S`$.

#### Set Cover.

A set cover is a subfamily $`𝒞^{}𝒞`$ such that $`_{S𝒞^{}}S=U`$ — in words, every element of $`U`$ is in some set in the cover. The (minimum-weight) set-cover problem is to find a set cover of minimum total weight $`w(𝒞^{})`$.

#### Element Packing.

An element packing is an assignment of non-negative weights to the elements such that the total weight assigned to the elements of any set $`S`$ is at most $`w(S)`$. The (maximum-weight) element-packing problem is to find an element packing maximizing the net weight assigned to elements.

### 1.4 Equivalences

The vertex cover problem in hypergraphs is equivalent to the minimum-weight set-cover problem as follows. For each $`S𝒞`$ we have a vertex $`v_S`$ in the hypergraph. For each element $`xU`$, we have an edge that contains $`v_S`$ if and only if $`xS`$. The number of edges $`m`$ is the number of elements. The rank $`r`$ is the maximum number of sets in which any element occurs. The size $`M`$ is the sum of the set sizes. The dual problems are also equivalent.

## 2 Reduction of Vertex Cover to $`ϵ`$-Maximal Packing

We first reduce our problem to the problem of finding what we call an $`ϵ`$-maximal packing. This reduction generalizes \[Ho82, BE81\], who considered $`ϵ=0`$.

###### Lemma 1 (Duality)

Let $`C`$ be an arbitrary vertex cover and $`p`$ an arbitrary edge packing. Then $`p(E)w(C)`$.

Proof:

$$p(E)=\underset{eE}{}p(e)\underset{eE}{}|eC|p(e)=\underset{vC}{}p(E(v))\underset{vC}{}w(v)=w(C).$$

###### Lemma 2 (Approximate Complimentary Slackness)

Let $`C`$ be a vertex cover and $`p`$ be a packing such that $`p(E(v))(1ϵ)w(v)`$ for every $`vC`$. Then $`(1ϵ)w(C)rp(E)`$. By duality, the weights of $`C`$ and $`p`$ are within a factor of $`r/(1ϵ)`$ from their respective optima.

Proof: Since $`(1ϵ)w(v)p(E(v))`$ for $`vC`$,

$$(1ϵ)w(C)=(1ϵ)\underset{vC}{}w(v)\underset{vC}{}p(E(v))=\underset{eE}{}|eC|p(e)rp(E).$$

We tighten Lemma 2 slightly when the weights are integers.

###### Lemma 3

In Lemma 2, if the weights are integers and $`ϵ<1/w(V)`$, then the weight of $`C`$ is at most $`r`$ times the minimum.

Proof: Let $`C^{}`$ be a minimum-weight cover. From Lemma 2, $`(1ϵ)w(C)rw(C^{})`$, so $`w(C)rw(C^{})+ϵw(C)=rw(C^{})`$.

Given a packing $`p`$, define $`C_p=\{vV:p(E(v))(1ϵ)w(v)\}`$. If $`C_p`$ is a vertex cover, then we say $`p`$ is $`ϵ`$-maximal. Note that $`p`$ is $`0`$-maximal if and only if $`p`$ is maximal. By Lemma 2, if $`p`$ is $`ϵ`$-maximal, then $`C_p`$ and $`p`$ are within a factor of $`r/(1ϵ)`$ from their respective optima.

## 3 The Algorithm

We have reduced the problem to finding an $`ϵ`$-maximal packing. The algorithm maintains a packing $`p`$ and the partial cover $`C_p=\{vV:p(E(v))(1ϵ)w(v)\}`$. The algorithm increases the individual $`p(e)`$’s until $`p`$ is $`ϵ`$-maximal and $`C_p`$ is a cover. When a vertex $`v`$ enters $`C_p`$, $`v`$ and the edges containing $`v`$ are deleted from the hypergraph. Let $`E_p`$ denote the set of remaining edges, let $`E_p(v)`$ denote the remaining edges incident to vertex $`v`$, and let $`d_p(v)`$ be the degree of $`v`$ in $`G_p=(V,E_p)`$. Define the residual weight $`w_p(v)`$ of vertex $`v`$ to be $`w(v)p(E(v))`$.

In a single round of the algorithm, for each remaining edge $`e`$, $`p(e)`$ is raised. To ensure that $`p`$ remains a packing, each vertex $`v`$ limits the increase in each $`p(e)`$ for $`ev`$ to at most $`w_p(v)/d_p(v)`$. Each $`p(e)`$ is then increased as much as possible subject to the limits imposed by all the $`ve`$. That is, each $`p(e)`$ is increased by $`\mathrm{min}_{ve}w_p(v)/d_p(v)`$. The algorithm repeats this basic round until $`p`$ converges to an $`ϵ`$-maximal packing. It then returns $`C_p`$. To implement the algorithm we maintain $`w_p`$ instead of $`p`$:

| Cover$`(G=(E,V),w,ϵ)`$ —  Returns a vertex cover of hypergraph $`G`$ of weight at most $`r/(1ϵ)`$ times the minimum. |
| --- |
| 1 | for | $`vV`$ par-do $`w_p(v)w(v)`$; $`E_p(v)E(v)`$; $`d_p(v)|E(v)|`$ |
| 2 | while | edges remain do |
| 3 |  | for each remaining edge $`e`$ par-do $`\delta (e)\mathrm{min}_{ve}w_p(v)/d_p(v)`$ |
| 4 |  | for | each remaining vertex $`v`$ par-do |
| 5 |  |  | $`w_p(v)w_p(v)_{eE_p(v)}\delta (e)`$ |
| 6 |  |  | if | $`w_p(v)ϵw(v)`$ then |
| 7 |  |  |  |  | delete $`v`$ and incident edges, updating $`E_p()`$ and $`d_p()`$ |
| 8 | return the set of deleted vertices |

As noted above, the limit on the increase in each $`p(e)`$ ensures that $`p`$ remains a packing. Consequently, the correctness and the approximation ratio of the algorithm are established by Lemmas 2 and 3. Using standard techniques \[Ja92\], each iteration of the while loop beginning with $`q`$ remaining edges can be done in $`O(\mathrm{ln}q)`$ time and $`O(rq)`$ operations on an EREW-PRAM.

## 4 Complexity Analysis

In this section, we prove our main theorem:

Main Theorem The algorithm requires $`O(r\mathrm{ln}^2m\mathrm{ln}\frac{1}{ϵ})`$ time and $`M/\mathrm{ln}^2m`$ processors, i.e., $`O(rM\mathrm{ln}\frac{1}{ϵ})`$ operations.

We use a potential function argument. Given a packing $`p`$, define

$$\varphi _p=\underset{vV}{}d_p(v)\mathrm{ln}\frac{w_p(v)}{ϵw(v)}.$$

The next lemma shows that during an iteration of the while loop $`\varphi _p`$ decreases by at least the number of edges remaining at the end of the loop. This is how we show progress.

###### Lemma 4

Let $`p`$ and $`p^{}`$, respectively, be the packing before and after an iteration of the while loop. Then $`\varphi _p\varphi _p^{}|E_p^{}|`$.

Proof: During the iteration, we say that a vertex $`v`$ limits an incident edge $`eE_p`$ if $`v`$ determines the minimum in the computation of $`\mathrm{min}_{ve}w_p(v)/d_p(v)`$. For each vertex $`v`$, let $`v`$ limit $`L(v)`$ edges, so that $`w_p^{}(v)w_p(v)(1L(v)/d_p(v)).`$ Let $`V^{}`$ denote the set of vertices that remain after the iteration. Then

$`\varphi _p\varphi _p^{}`$ $`=`$ $`{\displaystyle \underset{vV}{}}\left(d_p(v)\mathrm{ln}{\displaystyle \frac{w_p(v)}{ϵw(v)}}d_p^{}(v)\mathrm{ln}{\displaystyle \frac{w_p^{}(v)}{ϵw(v)}}\right)`$
$``$ $`{\displaystyle \underset{vV^{}}{}}d_p(v)\mathrm{ln}{\displaystyle \frac{w_p(v)}{w_p^{}(v)}}`$
$``$ $`{\displaystyle \underset{vV^{}}{}}d_p(v)\mathrm{ln}(1L(v)/d_p(v))`$
$``$ $`{\displaystyle \underset{vV^{}}{}}L(v)`$
$``$ $`|E_p^{}|`$
The second-to-last step follows because $`\mathrm{ln}(1x)x`$. The last step follows because each of the edges that remains is limited by some vertex in $`V^{}`$.

###### Lemma 5

There are at most $`(1+r\mathrm{ln}\frac{1}{ϵ})(1+\mathrm{ln}m)`$ iterations.

Proof: Let $`p`$ and $`p^{}`$, respectively, be the packing before and after any iteration. Let $`a=r\mathrm{ln}\frac{1}{ϵ}`$.

Clearly $`\varphi _p^{}|E_p^{}|a`$. By Lemma 4, $`\varphi _p^{}\varphi _p|E_p^{}|`$. Thus, $`\varphi _p^{}\varphi _p(11/(a+1))`$. Before the first iteration, $`\varphi _pma`$. Inductively, before the $`i`$th iteration,

$$\varphi _pma(11/(a+1))^{i1}ma\mathrm{exp}((i1)/(a+1)).$$

The last inequality follows from $`e^x1+x`$ for all $`x`$. Fixing $`i=1+(a+1)\mathrm{ln}m`$, we have $`\mathrm{exp}((i1)/(a+1))\mathrm{exp}(\mathrm{ln}m)=1/m`$, so before the $`i`$th iteration, $`\varphi _pa`$.

During each subsequent iteration, at least one edge remains, so $`\varphi _p`$ decreases by at least $`1`$. Thus, $`\varphi _p0`$ before an $`i+a`$th iteration can occur.

#### Time.

As each iteration requires $`O(\mathrm{ln}m)`$ time, the above lemma implies that the total time is $`O(r\mathrm{ln}^2m\mathrm{ln}\frac{1}{ϵ})`$.

#### Operations.

Recall that an iteration with $`q`$ edges requires $`O(rq)`$ operations. Consequently, the total number of operations is bounded by an amount proportional to $`r`$ times the sum, over all iterations, of the number of edges at the beginning of that iteration.

By Lemma 4, in a given iteration, $`\varphi _p`$ decreases by at least the number of edges remaining at the end of the iteration. Thus, the sum over all iterations of the number of edges during the iteration is at most $`m+\varphi _p`$ for the initial $`p`$. This is $`m+M\mathrm{ln}\frac{1}{ϵ}`$. Hence there are $`O(rM\mathrm{ln}\frac{1}{ϵ})`$ operations.

#### Processors.

Using standard techniques, the operations can be efficiently scheduled without increasing the time or the operations by more than a constant, so by the Work-Time Scheduling Principle \[Ja92\], the number of processors required is $`M/\mathrm{ln}^2m`$ — the work divided by the time. This establishes the Main Theorem.

### 4.1 The Intuition for Ordinary Graphs

The potential function analysis, while easy to verify, hides an interesting combinatorial principle that gives a good intuitive understanding of the algorithm for ordinary graphs. Recall that, during a single iteration of the while loop, a vertex $`v`$ limits an edge $`e`$ if $`v`$ determines the minimum in the calculation of $`\mathrm{min}_{ve}w_p(v)/d_p(v),`$ and that $`L(v)`$ denotes the number of edges limited by $`v`$. If a vertex $`v`$ limits at least a third of its incident edges, then $`w_p(v)`$ decreases by at least one third its value. Call such a vertex good. (After $`O(\mathrm{ln}\frac{1}{ϵ})`$ iterations of the while loop in which $`v`$ is good, $`v`$ will enter $`C_p`$.) In a given iteration, few vertices might be good. However, at least half of the remaining edges touch good vertices. This is a consequence of the following lemma:

###### Lemma 6 (\[II86, ABI86\])

Consider a directed graph. Call a vertex good if more than one-third of its incident edges are directed into it. Then at least half of the edges are directed into good vertices.

Proof: If a vertex is not good, call it bad. The in-degree of any bad vertex is at most half its out-degree, so the number of edges directed into bad vertices is at most half the number of edges directed out of bad vertices. Thus, the number of edges directed into bad vertices is at most half the number of edges. Thus, at least half the edges are directed into good vertices.

To see why the lemma applies, imagine directing each remaining edge into a vertex that limits it. Then the lemma shows that at least half of the remaining edges touch vertices that are good. Thus, in a given iteration, at least half the edges touch vertices whose residual weights decrease by more than a factor of $`1/3`$. This, intuitively, is why the algorithm makes progress.

This lemma is of independent interest: it drives the analyses of the running times of Israeli and Itai’s randomized maximal matching algorithm \[II86\] and of Alon, Babai, and Itai’s randomized maximal independent set algorithm \[ABI86\]. Interestingly, the natural generalization of the lemma to hypergraphs is not strong enough to give an analysis as tight as our potential function analysis.

### 4.2 Using Integer Arithmetic

If arithmetic precision is an issue, we can uniformly scale the original (integer) vertex weights so that the smallest weight is at least $`m/ϵ`$, and then use integer division (taking the floor) when computing the $`w_p(v)/d_p(v)`$’s. Essentially the same analysis carries through. (If $`w(v)m/ϵ`$, then $`w_p(v)m`$ while $`v`$ remains, so $`w_p(v)/d_p(v)1`$, hence $`w_p(v)/d_p(v)(w_p(v)/d_p(v))/2`$, and the net reduction in a $`w_p(v)`$ during an iteration is at least half what it would have been without taking the floor. Thus, the analysis will go through by doubling the potential function.)

Assuming without loss of generality that $`ϵ1/(2w(V))`$, if the original weights are $`k`$-bit integers, then the largest weight after scaling is bounded by

$$m\frac{1}{ϵ}2^{k+1}2^{k+2}mw(V)2^{2k+3}m|V|.$$

Hence the scaled weights are $`(2k+3+\mathrm{log}_2m+\mathrm{log}_2|V|)`$-bit integers. Subsequently all operations involve only integer arithmetic on smaller, non-negative integers.

Acknowledgments: We would like to thank Michael Luby and an anonymous referee for pointing out connections to the randomized algorithms of \[ABI86, Lu86\] and to the work by Hochbaum \[Ho82\], respectively.