Mon. Not. R. Astron. Soc. 000, 1–18 (2007)

Printed 7 January 2013

(MN LATEX style file v2.2)

Clustered Calibration: An Improvement to Radio Interferometric
Direction Dependent Self-Calibration
S. Kazemi1⋆, S. Yatawatta2 , S. Zaroubi1

1 Kapteyn

Astronomical Institute, University of Groningen, P.O. Box 800, 9700 AV Groningen, the Netherlands
Postbus 2, 7990 AA Dwingeloo, the Netherlands

arXiv:1301.0633v1 [astro-ph.IM] 3 Jan 2013

2 ASTRON,

7 January 2013

ABSTRACT

The new generation of radio synthesis arrays, such as LOFAR and SKA, have been designed to
surpass existing arrays in terms of sensitivity, angular resolution and frequency coverage. This
evolution has led to the development of advanced calibration techniques that ensure the delivery of accurate results at the lowest possible computational cost. However, the performance
of such calibration techniques is still limited by the compact, bright sources in the sky, used
as calibrators. It is important to have a bright enough source that is well distinguished from
the background noise level in order to achieve satisfactory results in calibration. This paper
presents “clustered calibration” as a modification to traditional radio interferometric calibration, in order to accommodate faint sources that are almost below the background noise level
into the calibration process. The main idea is to employ the information of the bright sources’
measured signals as an aid to calibrate fainter sources that are nearby the bright sources. In the
case where we do not have bright enough sources, a source cluster could act as a bright source
that can be distinguished from background noise. For this purpose, we construct a number
of source clusters assuming that the signals of the sources belonging to a single cluster are
corrupted by almost the same errors. Under this assumption, each cluster is calibrated as a
single source, using the combined coherencies of its sources simultaneously. This upgrades
the power of an individual faint source by the effective power of its cluster. The solutions thus
obtained for every cluster are assigned to each individual source in the cluster. We give performance analysis of clustered calibration to show the superiority of this approach compared
to the traditional un-clustered calibration. We also provide analytical criteria to choose the
optimum number of clusters for a given observation in an efficient manner.
Key words: methods: statistical, methods: numerical, techniques: interferometric

1 INTRODUCTION
Low frequency radio astronomy is undergoing a revolution as a
new generation of radio interferometers with increased sensitivity and resolution, such as the LOw Frequency ARray (LOFAR)1 ,
the Murchison Wide-field Array (MWA) 2 and the Square Kilometer Array (SKA)3 are being devised and some are already becoming operational. These arrays form a large effective aperture by the
combination of a large number of antennas using aperture synthesis (Thompson et al. 2001). In order to achieve the full potential
of such an interferometer, it is essential that the observed data is
properly calibrated before any imaging is done. Calibration of radio interferometers refers to the estimation and reduction of errors
introduced by the atmosphere and also by the receiver hardware,

⋆
1
2
3

E-mail: kazemi@astro.rug.nl
http://www.lofar.org
http://www.mwatelescope.org
http://www.skatelescope.org

c 2007 RAS

before imaging. We also consider the removal of bright sources
from the observed data part of calibration, that enable imaging the
weak background sources. For low frequency observations with a
wide field of view, calibration needs to be done along multiple directions in the sky. Proper calibration across the full field of view is
required to achieve the interferometer’s desired precision and sensitivity giving us high dynamic range images.
Early radio astronomy used external (classical) calibration
which estimates the instrumental unknown parameters using a radio source with known properties. This method has a relatively
low computational cost and a fast execution time, but it generates
results strongly affected by the accuracy with which the source
properties are known. In addition, existence of an isolated bright
source, which is the most important requirement of external calibration, in a very wide field of view is almost impractical, and
even when it does, external calibration would only give information along the direction of the calibrator. The external calibration
was then improved by self-calibration (Pearson & Readhead 1984).
Self-calibration has the advantage of estimating both the source

2

Kazemi et al.

and instrumental unknowns, without the need of having a prior
knowledge of the sky, only utilizing the instrument’s observed data.
Moreover, it can achieve a superior precision by iterating between
the sky and the instrumental parameters.
The accuracy of any calibration scheme, regardless of the
used technique, is determined by the level of Signal to Noise Ratio (SNR). This limits the feasibility of any calibration scheme
to only bright sources that have a high enough SNR to be distinguished from the background noise (Bernardi et al. 2010; Liu et al.
2009; Pindor et al. 2010). Note that interferometric images are produced using the data observed during the total observation (integration) time. However, calibration is done at shorter time intervals of that total duration. This increases the noise level of the data
for which calibration is executed compared to the one in the images. In other words, there are some faint sources that appear well
above the noise in the images while they are well below the noise
level during calibration. It has been a great challenge to calibrate
such faint sources having a very low SNR. In this paper we present
clustered self-calibration, introduced in Kazemi et al. (2011a), and
emphasize its better performance compared to un-clustered calibration below the calibration noise level. Existing un-clustered calibration approaches (Intema et al. 2009; Smirnov 2011) can only handle a handful of directions where strong enough sources are present
and we improve this, in particular for subtraction of thousands of
sources over hundreds of directions in the sky, as in the case for
LOFAR (Bregman 2012).
The implementation of clustered calibration is performed by
clustering the sources in the sky, assuming that all the sources in a
single cluster have the same corruptions, and calibrating each cluster as a single source. At the end, the obtained calibration solutions
for every source cluster is assigned to all the sources in that cluster.
This procedure improves the information used by calibration by incorporating the total of signals observed at each cluster instead of
each individual source’s signal. Therefore, when SNR is very low, it
provides a considerably better result compared to un-clustered calibration. Moreover, calibrating for a set of source clusters instead
of for all the individual sources in the sky reduces the number of
directions along which calibration has to be performed, thus considerably cutting down the computational cost. However, there is
one drawback in clustered calibration: The corruptions of each individual source in one given cluster will almost surely be slightly
different from the corruptions estimated for the whole cluster by
calibration. We call this additional error as “clustering error” and
in order to get the best performance in clustered calibration, we
should find the right balance between the improvement in SNR as
opposed to degradation by clustering error.
Recently, clustering methods have gained a lot of popularity in
dealing with large data sets (Kaufman & Rousseeuw 1990). There
are various clustering approaches that could be applied to calibration based on their specific characteristics. An overview of different
clustering methods is given in Xu & Wunsch (2008). Clustering of
radio sources should take into account (i) their physical distance
from each other and (ii) their individual intensity. The smaller the
angular separations of sources are, the higher the likelihood that
they share the same corruptions in their radiated signals. Moreover, in order to get the best accuracy in the calibration results,
there should be a balance between effective intensities of different
clusters. Thus, in the clustering procedure, every source should be
weighted suitably according to its brightness intensity.
The brightness distribution of radio source in the sky is a
power law and the spatial distribution is Poisson. Therefore, clustering the sources via probabilistic clustering approaches is com-

putationally complex. Hierarchical clustering (Johnson 1967) is a
well-known clustering approach with a straight forward implementation suitable for our case. But, its computational cost grows exponentially with the size of its target data set which can be a disadvantage when the number of sources is huge. Weighted K-means
clustering (Kerdprasop et al. 2005; MacQueen 1967) is also one of
the most used of clustering schemes applicable in clustered calibration. The advantage of this clustering technique is its low computational cost, which is proportional to the number of clusters and the
size of the target data set. Nonetheless, we emphasize that the computational time taken by any of the aforementioned clustering algorithms is negligible compared with the computational time taken
by the actual calibration. Therefore, we pursue all clustering approaches in this paper. However, the use of Fuzzy C-means clustering (Bezdek 1981) in clustered calibration requires major changes
in the calibration data model and will be explored in future work.
This paper is organized as follows: First, in section 2 we
present the general data model used in clustered calibration. In section 3, we present modified weighted K-means and divisive hierarchical clustering for clustering sources in the sky. Next, in section 4, we focus on analyzing the performance of clustered calibration, with an a priori clustered sky model, and compare the results with un-clustered calibration. In clustered calibration, there is
contention between the improvement of SNR by clustering sources
and the errors introduced by the clustering of sources itself. Thus,
we relate the clustered calibration’s performance to the effective
Signal to Interference plus Noise Ratio (SINR) obtained for each
cluster. For this purpose, we use statistical estimation theory and
the Cramer-Rao Lower Bounds (CRLB) (Kay 1993). In section 5,
we derive criteria for finding the optimum number of clusters for a
given sky. We use the SINR analysis and adopt Akaike’s Information Criterion (AIC) (Akaike 1973) and the Likelihood Ratio Test
(LRT) (Graves 1978) to estimate the optimum number of clusters.
We present simulation results in section 6 to show the superiority
of clustered calibration to un-clustered calibration and the performance of the presented criteria in detecting the optimum number
of clusters. Finally, we draw our conclusions in section 7. Through
this paper, calibration is executed by the Space Alternating Generalized Expectation maximization (SAGE) (Fessler & Hero 1994;
Yatawatta et al. 2009; Kazemi et al. 2011b) algorithm. Moreover,
in our simulations, radio sources are considered to be uniformly
distributed in the sky and their flux intensities follow Raleigh distribution, which is the worst case scenario. In real sky models, there
usually exist only a few (two or three) number of bright sources
which dominate the emission. In the presence of these sources and
the background noise, it is impractical to solve for the other faint
sources in the field of view individually. Therefore, obtaining a better calibration performance via the clustered calibration approach,
compared to the un-clustered one, is guaranteed. Therefore, in section 6 we illustrate this using simulations in which the brightness
distribution of sources is a power law with a very steep slope. On
top of that, Kazemi et al. (2012); Yatawatta et al. (submitted 2012)
also present the performance of clustered calibration on real observations using LOFAR.
The following notations are used in this paper: Bold, lowercase letters refer to column vectors, e.g., y. Upper case bold letters
refer to matrices, e.g., C. All parameters are complex numbers, unless stated otherwise. The inverse, transpose, Hermitian transpose,
and conjugation of a matrix are presented by (.)−1 , (.)T , (.)H , and
(.)∗ , respectively. The statistical expectation operator is referred to
as E{.}. The matrix Kronecker product and the proper (strict) subset are denoted by ⊗ and (, respectively. In is the n × n identity
c 2007 RAS, MNRAS 000, 1–18

Clustered Radio Interferometric Calibration
matrix and ∅ is the empty set. The Kronecker delta function is presented by δij . R and C are the sets of Real and Complex numbers,
respectively. The Frobenius norm is shown by ||.||. Estimated pac All logarithmic calculations are
rameters are denoted by a hat, (.).
to the base e. The multivariate Gaussian and Complex Gaussian
distributions are denoted by N and CN , respectively.

2 CLUSTERED SELF-CALIBRATION DATA MODEL
In this section, we present the measurement equation of a polarimetric clustered calibration in detail (Hamaker et al. 1996;
Hamaker 2006). Suppose we have a radio interferometer consisting of N polarimetric antennas where each antenna is composed
of two orthogonal dual-polarization feeds that observe K compact
sources in the sky. Every i-th source signal, i ∈ {1, 2, . . . , K},
e pi = [vXpi vY pi ]T at X and Y
causes an induced voltage of v
dipoles of every p-th antenna, p ∈ {1, 2, . . . , N }. In practice,
T

epi = Jpi ei ,
v

(1)

where ei = [eXi eY i ] is the source’s electric field vector and
Jpi represents the 2 × 2 Jones matrix (Hamaker et al. 1996) corresponding to the direction-dependent gain corruptions in the radiated signal. These corruptions are originated from the instrumental
(the beam shape, system frequency response, etc.) and the propagation (tropospheric and ionospheric distortions, etc.) properties
which later on, in this section, will be explained in more detail.
The signal vp obtained at every antenna p is a linear suepi where i ∈
perposition of the K sources corrupted signals, v
{1, 2, . . . , K}, plus the antenna’s thermal noise. The multitude of
ignored fainter sources also contributes to this additive noise.
The voltages collected at the instrument antennas get corrected
for geometric delays, based on the location of their antennas, and
some instrumental effects, like the antenna clock phases and electronic gains. Then, they are correlated in the array’s correlator to
generate visibilities (Hamaker et al. 1996). The visibility matrix of
the baseline p − q, E{vp ⊗ vqH }, is given by
!
K
X
H
Vpq = Gp
Jpi (θθ )Ci{pq} Jqi (θθ ) GH
(2)
q + Npq .
i=1

In Eq. (2), θ ∈ CP , P = 4KN , is the unknown instrumental
and sky parameter vector, Npq is the additive 2 × 2 noise matrix
of the baseline p − q, and Ci{pq} is the Fourier transform of the
i-th source coherency matrix Ci = E{ei ⊗ eH
i } (Born & Wolf
1999; Hamaker et al. 1996). If the i-th source radiation intensity
is Ii , then Ci = I2i I2 . Considering this source to have equatorial
coordinate, (Right Ascension α, Declination δ), equal to (αi , δi ),
and the geometric components of baseline p − q to be (u, v, w),
then
√
2
2
(3)
Ci{pq} = e−2π(ul+vm+w( 1−l −m −1)) Ci ,
where 2 = −1 and,
l = sin(αi − α0 )cos(δi ),

m = cos(δ0 )sin(δi ) − cos(αi − α0 )cos(δi )sin(δ0 ),
are the source direction components corresponding to the observation phase reference of (α0 , δ0 ) (Thompson et al. 2001). The errors
common to all directions, such as the receiver delay and amplitude
errors, are given by Gp and Gq . Initial calibration at a finer time
c 2007 RAS, MNRAS 000, 1–18

3

and frequency resolution is performed to estimate and correct for
Gp -s and the corrected visibilities are obtained as
−H
e pq = G−1
V
.
p Vpq Gq

(4)

The remaining errors are unique to a given direction, but residual
errors in Gp -s are also absorbed into these errors, which are denoted by Jpi in the usual notation.
Vectorizing Eq. (4), the final visibility vector of the baseline
p − q is given by
vpq =

K
X
i=1

J∗qi (θθ ) ⊗ Jpi (θθ )vec(Ci{pq} ) + npq .

(5)

Stacking up all the cross correlations (measured visibilities) and
noise vectors in y and n, respectively, the un-clustered selfcalibration measurement equation is given by
y=

K
X

si (θθ ) + n.

(6)

i=1

In Eq. (6), y, n ∈ CM , M = 2N (N − 1), the noise vector is considered to have a zero mean Gaussian distribution with covariance
Π , n ∼ N (0, Π M ×M ), and the nonlinear function si (θθ ) is defined
as


J∗2i (θθ ) ⊗ J1i (θθ )vec(Ci{pq} )


J∗3i (θθ ) ⊗ J1i (θθ )vec(Ci{pq} )


(7)
si (θθ ) ≡ 
.
..


.
J∗Ni (θθ ) ⊗ J(N−1)i (θθ )vec(Ci{pq} )

Calibration is essentially the Maximum Likelihood (ML) estimation of the unknown parameters θ (P complex values or 2P real
values), or of the Jones matrices J(θθ ), from Eq. (6) and removal of
the K sources. Note that calibration methods could also be applied
to the uncorrected visibilities of (4) to estimate Gp and Gq errors
as well.
The Jones matrix Jpi , for every i-th direction and at every p-th
antenna, is given as
Jpi ≡ Epi Zpi Fpi .

(8)

In Eq. (8), Epi , Zpi , and Fpi are the antenna’s voltage pattern,
ionospheric phase fluctuation, and Faraday Rotation Jones matrices, respectively. In practice, the E, Z, and F Jones matrices obtained for nearby directions and for a given antenna are almost the
same. Thus, for every antenna p, if the i-th and j-th sources have a
small angular separation from each other, we may consider
Jpi ∼
= Jpj .

(9)

This is the underlying assumption for clustered calibration.
Clustered calibration first assigns source clusters, Li for i ∈
{1, 2, . . . , Q} where Q ≪ K, on which the sky variation is considepi which
ered to be uniform. Then, it assumes there exists a unique J
is shared by all the sources of the i-th cluster Li , i ∈ {1, 2, . . . , Q},
at receiver p, p ∈ {1, 2, . . . , N }. Based on that, the visibility at every baseline p − q, given by Eq. (4), is reformulated as
e pq =
V

Q
X
i=1

epi (θe){
J

X

l∈Li

e
eH
e
Cl{pq} }J
qi (θ ) + Npq .

(10)

e pq is the clustered calibration’s effective noise at
In Eq. (10), N
baseline p − q which will be explicitly discussed at section 4. Note
that clustered calibration estimates the new unknown parameter

4

Kazemi et al.

e
θe ∈ CP where Pe = 4QN . Denoting the effective signal of every i-th cluster at baseline p − q by
X
e i{pq} ≡
C
Cl{pq} ,
(11)
l∈Li

the clustered calibration visibility vector at this baseline (vectorized
form of Eq. (10)) is
vpq =

Q
X
i=1

epi (θe)vec(C
e i{pq} ) + n
e∗qi (θe) ⊗ J
e pq .
J

(12)

Finally, stacking up the visibilities of all the instrument’s baselines
in vector y, clustered calibration’s general measurement equation
is resulted as
y=

Q
X
i=1

e
e.
si (θe) + n

(13)

In Eq. (13) e
si is defined similar to si in Eq. (7) where J and C are
e and C,
e respectively.
replaced by J
Because of the similarity between the clustered and the unclustered calibration’s measurement equations presented by Eq.
(13) and Eq. (6), respectively, they could utilize the same calibration techniques. Thus, the only difference between these two
types of calibration is that clustered calibration solves for clusters
of sources instead of for the individual ones. That upgrades the signals which should be calibrated for from Eq. (3) to Eq. (11).

3 CLUSTERING ALGORITHMS
Clustering is grouping a set of data so that the members of the same
group (cluster) have some similarities (Kaufman & Rousseeuw
1990). This similarity is defined based on the application of the
clustering method.
We need to define clustering schemes in which two radio
sources merge to a single cluster based on the similarity in their
direction dependent gain errors ( Eq. (8) ). Radiations of sources
that are close enough to each other in the sky are assumed to be
affected by almost the same corruptions ( Eq. (9) ). Based on that,
we aim to design source clusters with small angular diameters. On
the other hand, every cluster’s intensity is the sum of the intensities
of its members ( Eq. (11) ). In order to keep a balance between different clusters’ intensities, we intend to apply weighted clustering
techniques in which the sources are weighted proportional to their
intensities.
Suppose that the K sources, x1 , x2 , . . . , xK have
(α1 , δ1 ), (α2 , δ2 ), . . . , (αK , δK ) equatorial coordinates, respectively. TheP
aim is to provide Q clusters so that the objective
Q
function f =
q=1 D(Lq ) is minimized. D(Lq ) is the angular
diameter of cluster Lq , for q ∈ {1, 2, . . . , Q}, defined as
D(Lq ) ≡ max {d(xi , xj )|xi , xj ∈ Lq },

For defining the centroids, we associate a weight to every source
xi , as
wi = w(xi ) ≡

Ii
,
I∗

for i ∈ {1, 2, . . . , K},

(16)

where Ii is the source’s intensity and I ∗ = min {I1 , I2 , . . . , IK }.
Applying the weights to the clustering procedure, the centroids of
the clusters lean mostly towards the brightest sources. That causes
a tendency in faint sources to gather with brighter sources close
to them into one cluster. Thus, their weak signals are promoted
being added up with some brighter sources’ signals. Moreover, very
strong sources will be isolated such that their signals are calibrated
individually, without being affected by the other faint sources.
We cluster radio sources using weighted K-means
(Kerdprasop et al. 2005) and divisive hierarchical clustering
(Johnson 1967) algorithms. Since the source clustering for calibration is performed offline, its computational complexity is
negligible compared with the calibration procedure itself.
3.1 Weighted K-means clustering
Step 1. Select the Q brightest sources, x1∗ , x2∗ , . . . , xQ∗ , and initialize the centroids of Q clusters by their locations as
for q ∈ {1, 2, . . . , Q}, q ∗ ∈ {1∗ , 2∗ , . . . , Q∗ }.
(17)
Step 2. Assign each source to the cluster with the closest centroid,
defining the membership function
(
1, if d(xi , cq ) = min{d(xi , cj )|j = 1, . . . , Q}
mLq (xi ) =
0, Otherwise
cq ≡ [αq∗ , δq∗ ],

Step 3. Update the centroids by
PK
mLq (xi ) wi xi
cq = Pi=1
,
K
i=1 mLq (xi ) wi

for q ∈ {1, 2, . . . , Q}.

(18)

Repeat steps 2 and 3 until there are no reassignments of sources to
clusters.
3.2 Divisive hierarchical clustering
Step 1. Initialize the cluster counter Q′ to 1, assign all the K
sources to a single cluster L1 and ∅ to a set of null clusters A.
Step 2. Choose cluster Lq∗ , for q ∗ ∈ {1, 2, . . . , Q′ } − A, with the
largest angular diameter
D(Lq∗ ) = max{D(Lq )|q ∈ {1, 2, . . . , Q′ } − A}.

(19)

Step 3. Apply the presented weighted K-means clustering technique to split Lq∗ into two clusters, L′q∗ and L′′q∗ .
Step 4. If D(L′q∗ ) + D(L′′q∗ ) < D(Lq∗ ), then set Q′ = Q′ + 1,
Lq∗ ≡ L′q∗ , LQ′ ≡ L′′q∗ , and A = ∅, otherwise set A = A∪{q ∗ }.
Repeat steps 2, 3, and 4 until Q′ = Q.

(14)

and d(., .) is the angular separation between any two points on the
celestial sphere. Having two radio sources a and b with equatorial
coordinates (αa , δa ) and (αb , δb ), respectively, the angular separation d(a, b), in radians, is obtained by
p
cos2 δb sin2 ∆α + [cosδa sinδb − sinδa cosδb cos∆α]2
−1
tan
,
sinδa sinδb + cosδa cosδb cos∆α
(15)
where ∆α = αb − αa .

3.3 Comparison of clustering methods
Hierarchical clustering method tends to design clusters with almost the same angular diameters, whereas, the K-means clustering
method tends to keep the same level of intensity at all its clusters.
In practice, since hierarchical clustering method makes less errors
in dedicating the same solutions to sources in small clusters, it performs better than Weighted K-means clustering in a clustered calibration procedure. But, when the number of sources (and clusters)
c 2007 RAS, MNRAS 000, 1–18

42

42

41

41

40

40
δ [deg.]

δ [deg.]

Clustered Radio Interferometric Calibration

39
38

39
38

37

37

36

36

35

278

276

274
α [deg.]

272

35

270

Figure 1. A simulated 8 by 8 degrees sky of fifty point sources with intensities below 3 Jy. The source positions and their brightness are following
uniform and Rayleigh distributions, respectively. The marker sizes are proportional to sources intensities.

In this section, we explain the reasons for clustered calibration’s
better performance, compared to un-clustered calibration, at a low
c 2007 RAS, MNRAS 000, 1–18

276

274
α [deg.]

272

270

42
41

3.3.1 Example 1: Weighted K-means and hierarchical clustering

4 PERFORMANCE ANALYSIS

278

Figure 2. Fifty point sources are clustered into ten source clusters by
Weighted K-means clustering technique. There is not a good balance between different clusters angular diameters.

is very large (Q > 100), its prohibitive computational costs makes
the fast K-means clustering method preferable.

40
δ [deg.]

We simulate an 8 by 8 degrees sky with fifty point sources with intensities below 3 Jansky (Jy). The source positions and their brightness follow uniform and Rayleigh distributions, respectively. The
result is shown by Fig. 1 in which the symbol sizes are proportional
to intensities of sources. Weighted K-means and divisive hierarchical clustering methods are applied to cluster the fifty sources into
ten source clusters. The results are presented in Fig. 2 and Fig. 3,
respectively. Fig. 2 shows that the Weighted K-means clustering
could design source clusters with considerably large angular diameters. Assigning the same calibration solutions to the sources of
these large clusters could cause significant errors. However, as Fig.
3 shows, this is not the case for the hierarchical clustering and it
constructs clusters with almost the same angular diameters.
Since the number of sources in this simulation is not that large
(K = 50), the difference between execution time of the two clustering methods is not significant. Hence in such a case, the use of hierarchical clustering method, rather than the Weighted K-means, is
advised. However, this is not the case when we have a large number
of sources, and subsequently a large number of source clusters, in
the sky. To demonstrate this, we use the two clustering techniques
for clustering thousand of sources (K = 1000) into Q source clusters, Q ∈ {3, 4, . . . , 100}. The methods’ computational times versus the number of clusters are plotted in Fig. 4. As Fig. 4 shows,
for large Qs, the computational cost of Weighted K-means is much
cheaper than the one of the hierarchical clustering. That can make
the Weighted K-means clustering method more suitable than the
hierarchical clustering for such a case.

5

39
38
37
36
35

278

276

274
α [deg.]

272

270

Figure 3. Fifty point sources are clustered into ten source clusters via hierarchical clustering method. Different clusters have almost the same angular
diameters.

SNR (Kazemi et al. 2012). This superiority is in the sense of
achieving an unprecedented precision in solutions with a considerably low computational complexity, given the optimum clustering
scheme. In the next section, we present different criteria for finding
the optimum number of clusters at which the clustered calibration
performs the best.

4.1 Cramer-Rao Lower Bounds
The most fundamental assumption in clustered calibration is that
the sources at the same cluster have exactly the same corruptions
in their radiated signals. This assumption is of course incorrect,
nonetheless, it provides us with a stronger signal, the sum of the
signals in the whole cluster. We present an analytic comparison of
clustered and un-clustered calibration where we use the Cramer-

6

Kazemi et al.

Computational time (s)

800

In other words, the variance of any unbiased estimator of the unknown parameter vector θ is bounded from below by the diagonal
elements of [I(θθ )]−1 .
Using (23) and (26), the Fisher information matrix of the visibility vector y is obtained as

Weighted K−means
Divisive Hierarchical

600

I(θθ ) = 2σ −2 Re(JsH Js ),
400

where Js is the Jacobian matrix of s with respect to θ
Js (θθ ) =

200
0
0

(27)

20

40
60
Number of clusters

80

100

Figure 4. Weighted K-means and divisive hierarchical clustering methods
computational costs. For small number of source clusters, there is no difference between execution times of the two clustering methods. But, when
the number of source clusters is large, the computational cost of Weighted
K-means becomes much cheaper than the one of the hierarchical clustering.

Rao Lower Bound (CRLB) (Kay 1993) as a tool to measure the
performance of the calibration.
4.1.1 Estimations of CRLB for two sources at a single cluster
For simplicity, first consider observing two point sources at a single
baseline, lets say baseline p − q. Based on Eq. (4), the visibilities
are given by
H
e pq = Jp1 C1{pq} JH
V
q1 + Jp2 C2{pq} Jq2 + Npq ,

(20)

e pq , the visiin the un-clustered calibration strategy. Vectorizing V
bility vector is
y = J∗q1 ⊗ Jp1 vec(C1{pq} ) + J∗q2 ⊗ Jp2 vec(C2{pq} ) + npq .

Assuming npq ∼ CN (0, σ 2 I4 ), we have

y ∼ CN (s(θθ ), σ 2 I4 ),

(21)

where
s(θθ ) ≡

X

i=1,2

J∗qi (θθ ) ⊗ Jpi (θθ )vec(Ci{pq} ).

(22)

Using Eq. (21), the log-likelihood function of the visibility vector
y is given by
π
L(θθ |y) = −4ln{ 2 } − σ −2 (y − s(θθ ))H (y − s(θθ )).
(23)
σ
CRLB is a tight lower bound on the error variance of any unbiased parameter estimators (Kay 1993). Based on its definition, if
the log-likelihood function of the random vector y, L(θθ |y), satisfies the “regularity” condition


∂
Ey
L(θθ |y) = 0, for all θ ,
(24)
∂θθ
then for any unbiased estimator of θ , θb,
var(θbi ) > [I −1 (θθ )]ii ,

for i ∈ {1, . . . , M },

where I(θθ ) is the Fisher information matrix defined as
 2

∂ L(θθ |y)
I(θθ ) = − Ey
.
∂θθ ∂θθ T

(25)

(26)

2
X
∂
{J∗ ⊗ Jpi }[I4 ⊗ vec(Ci{pq} )].
θ qi
∂θ
i=1

(28)

Thus, variations of any unbiased estimator of parameter vector θ ,
lets say θb, is bounded from below by the CRLB as
Var(θb) > [2σ −2 Re(JsH Js )]−1 .

(29)

Lets try to bound the error variations of the clustered calibration parameters assuming that the two sources construct a single
cluster, called cluster number 1. We reform Eq. (20) as
e p1 (C1{pq} +C2{pq} )J
eH
Vpq = J
q1 +Γ1{pq} +Γ2{pq} +Npq , (30)

where Γi{pq} , referred to as the “clustering error” matrices, are
given by
e
eH
Γi{pq} = Jpi Ci{pq} JH
qi − Jp1 Ci{pq} Jq1 ,

(31)

e pq ≡ Γ1{pq} + Γ2{pq} + Npq .
N

(32)

e∗q1 ⊗ J
ep1 vec(C1{pq} + C2{pq} ) + n
e pq ,
y=J

(33)

ep1 (θe) is the clustered calibration solution at receiver p.
and J
e pq in
Eq. (30) implies that what is considered as the noise matrix N
the clustered calibration data model, Eq. (10), is in fact
Vectorizing Eq. (30), the clustered calibration visibility vector is
obtained by

e pq ).
e pq = vec(N
where n
We point out that depending on the observation as well as the
positions of the two sources on the sky, the clustering error Γi{pq}
will have different properties. However, in order to study the performance of the clustered calibration in a statistical sense, and to
simplify our analysis, we make the following assumptions.
(i) Consider statistical expectation over different observations
and over different sky realizations where the sources are randomly
e → E{J}
distributed on the sky. In that case, almost surely E{J}
and consequently
E{Γi{pq} } → 0.

(34)

In other words, we assume the clustering error to have zero mean
over many observations of different parts of the sky.
(ii) We assume that the closer the sources are together in the sky,
the smaller the errors introduced by clustering would be. Therefore,
given a set of sources, the clustering error will reduce as the number of clusters increase. In fact this error introduced by clustering
vanishes when the number of clusters is equal to the number of
sources (each cluster contains only one source). Therefore, given a
set of sources, the variance of Γi{pq} will decrease as the number
of clusters increase.
Using Eq. (34) and bearing in mind that E{Npq } = 0, we
eH
e pq ∼ CN (0, σ
e 2 I4 .
can consider n
e 2 I4 ) where E{e
npq n
pq } = σ
Therefore,
y ∼ CN (e
s, σ
e 2 I4 ),

c 2007 RAS, MNRAS 000, 1–18

Clustered Radio Interferometric Calibration

The CRLBs obtained for the un-clustered and clustered calibrations in Eq. (29) and Eq. (35), respectively, are both almost
equal to the inverse of the Signal to Interference plus Noise Ratio
(SINR), SINR−1 . In un-clustered calibration, the effective signal
for the faintest source is C2{pq} where the noise is Npq . Therefore, SINR for this source is

2000
Clustered calibration
Un−clustered calibration

CRLB

1500

1000

||C2{pq} ||2
.
||Npq ||2

SINR2 =

(36)

But, in clustered calibration, the effective signal and noise are
e {pq} ≡ C1{pq} + C2{pq} and N
e pq , respectively. Thus, SINR
C
for the cluster is

500

0
0

7

1000
2000
3000
4000
Clustered calibration effective noise power

5000

Figure 5. Clustered and un-clustered calibrations CRLB. When the effece 2 , is small enough, then
tive noise power of the clustered calibration, ||N||
its CRLB is lower than of the un-clustered calibration’s and it reveals a
superior performance.

e∗q1 ⊗ J
e p1 vec(C1{pq} + C2{pq} ),
e
s≡J

and similar to Eq. (29), we have

b
σ −2 Re(JesH Jes )]−1 .
Var(θe) > [2e

(35)

We use numerical simulations to compare the un-clustered and
clustered calibrations performances via their CRLBs which are
given by Eq. (29) and Eq. (35), respectively.

4.1.2 Example 2: CRLB for two sources and one cluster
We simulate a twelve hour observation of two point sources with
intensities I1 = 11.25 and I2 = 2.01 Jy at sky coordinates (l, m)
equal to (−0.014, −0.005) and (−0.011, −0.010) radians, respectively. We use the uv-coverage of Westerbork Synthesis Radio Telescope (WSRT) with 14 receivers in this simulation.
We Consider the J Jones matrices in Eq. (20) to be diagonal.
Their amplitude and phase elements follow U(0.75, 0.95) and
U(0.003, 0.004) distributions, respectively. The background noise
is N ∼ CN (0, 10I). Jones matrices of the clustered calibration,
e p1 for p = 1, 2, are obtained as Jp1 + U(0.02, 0.40)ejU (0.5,2) .
J
e matrices, we calculated CRLB of the unFor 20 realizations of J
clustered and clustered calibrations using Eq. (29) and Eq. (35),
respectively. The results are presented in Fig. 5. As shown in this
figure, for small enough errors matrices Γ of Eq. (31), the clustered calibration’s CRLB stands below the un-clustered calibration’s CRLB. On the other hand, with increasing power of error mae the un-clustered calibratrices, or the power of effective noise N,
tion’s CRLB becomes lower than the clustered calibration’s CRLB.
4.1.3 Analysis of CRLB
Generally, if source 1 is considerably brighter than source 2,
||C1{pq} || ≫ ||C2{pq} ||, and if the weak source power is much
lower than the noise level, ||C2{pq} || ≪ ||Npq ||, then clustered
calibration’s performance is better than un-clustered calibration.
Note that the worst performance of both calibrations is at the
faintest source and we are more concerned to compare the CRLBs
for this source.
c 2007 RAS, MNRAS 000, 1–18

SINRc =

e {pq} ||2
||C
.
e pq ||2
||N

(37)

Clustered calibration has an improved performance when
SINRc ≫ SINR2 .

(38)

Consider the two possible extremes in a clustered calibration
procedure:
(i) Clustering many sources in a large field of view to a very
small number of clusters. In this case, the angular diameter of a
cluster is probably too large for the assumption of uniform corruptions to apply. Subsequently, dedicating a single solution to all the
sources of every cluster by clustered calibration introduces clustering error matrices Γ with a large variance (see Eq. (31)). Having
high interference power, the clustered calibration effective noise
e of Eq. (32) becomes very large. Therefore, clustered calibration
N
SINR will be very low and it does not produce high quality results.
(ii) Clustering sources in a small field of view to a very large
number of clusters. In this case, the variance of Γ matrices are almost zero while the signal powers of source clusters are almost as
low as of the individual sources. Therefore, the SINR of clustered
calibration is almost equal to the un-clustered calibration SINR and
the calibration performance is expected to be almost the same as
well.
Thus, the best efficiency of clustered calibration is obtained at the
smallest number of clusters for which Eq. (38) is satisfied. We use
the SINR of Eq. (38) as an efficient criterion for detecting the optimum number of clusters.

4.1.4 Generalization to many sources and many clusters
For the visibility vector y of un-clustered calibration’s general data
model, presented by Eq. (6), we have
y ∼ CN (

K
X

si (θθ ), Π ).

(39)

i=1

Therefore, the CRLB of un-clustered calibration is

X
−1
K
K
X
var(θθ ) > 2 Re (
,
Jsi (θθ ))H Π −1 (
Jsi (θθ ))
i=1

i=1

where Jsi is the Jacobian matrix of si with respect to θ .
Computing the exact CRLB is more complicated when we
have clustered calibration. In the clustered calibration measurement
equation, given by Eq. (13), we have
e≡
n

K
X
i=1

Γ i + n,

(40)

8

Kazemi et al.
−10

where n is the un-clustered calibration’s noise vector,
(41)

and Γ i{pq} is given by Eq. (31). Due to the existence of the nuisance
parameters Γ i in the clustered calibration data model, calculation
of its conventional CRLB is impractical. This leads us to the use of
Cramer-Rao like bounds devised in the presence of the nuisance parameters (Gini & Reggiannini 2000). We apply the Modified CRLB
(MCRLB) (Gini et al. 1998) to the performance of clustered calibration.
b̃
The MCRLB for estimating the errors of θ in the presence of
the nuisance parameters Γ (clustering error) is defined as



−1
∂ ∂
b̃
Γ
θ
var(θ ) > Ey,ΓΓ − Ey|ΓΓ
ln{P
(y|Γ
;
θ̃
)}
,
T
∂θ̃θ ∂θ̃θ
(42)
Γ; θ̃θ ) is the Probability Density Function (PDF) of the
where P (y|Γ
visibility vector y assuming that the Γ matrices of Eq. (41) are a
priori known. Since n ∼ CN (0, Π ), from Eq. (13) we have
Q
K
X
X
Γ ∼ CN ([
Γ i ], Π ),
y|Γ
s̃i +

2
1
0



∂
∂
T
θ ∂θ̃
∂θ̃
θ


Γ
θ
ln{P (y|Γ ; θ̃ )} , which

20
30
Number of clusters

40

50

−6

10

Q
Q
K
K
X
X
X
X
JΓ i (θ̃θ )]}.
JΓ i (θ̃θ )]H Π −1 [
Js̃i (θ̃θ ) +
2Re{[
Js̃i (θ̃θ ) +
i=1

Clustered calibration
Un−clustered calibration

i=1

Note that Ey,ΓΓ in Eq. (42) could be estimated by Monte-Carlo
method.
As a rule of thumb, reducing the heavy computational cost of
MCRLB, one can interpret the SINR test of Eq. (38) as follows: If
in average the effective SINR of clustered calibration, SINRc , gets
higher than the effective SINR of un-clustered calibration obtained
for the weakest observed signal, SINRw ,
E{SINRc } ≫ E{SINRw },

x 10

9
8
SINR

i=1

10

Figure 6. Clustered calibration’s MCRLB. For very small Q, where the
effect of interference is large, MCRLB is high. By increasing the number
of clusters, MCRLB decreases and reaches its minimum where the best
performance of the clustered calibration is expected. After that, due to the
dominant effect of the background noise, MCRLB starts to increase until it
reaches the un-clustered calibration CRLB.

is called the modified Fisher information, is equal to

i=1

3

i=1

i=1

and therefore in Eq. (42), − Ey|ΓΓ

(43)

Clustered calibration
Un−clustered calibration

4
MCRLB

Γi{12} )T . . . vec(Γ
Γi{(N−1)N} )T ]T ,
Γ i = [vec(Γ

5 x 10

7
6
5
4

(44)

then clustered calibration can achieve a better results. In Eq. (44),
the expectation is taken with respect to the noise N, error matrices
Γ , and all the baselines.
4.1.5 Example 3: MCRLB and SINR estimations
We simulate WSRT including N = 14 receivers which observe
fifty sources with intensities below 15 Jy. The source positions
and their brightness follow uniform and Rayleigh distributions, respectively. The background noise is N ∼ CN (0, 15IM ), where
M = 2N (N − 1) = 364. We cluster sources using divisive
hierarchical clustering, into Q number of clusters where Q ∈
e are gener{3, 4, . . . , 50}. Clustered calibration’s Jones matrices, J,
jU (0,0.2)
ated as U(0.9, 1.1)e
. Since for smaller number of clusters,
we expect larger interference (errors)
soP in clustered calibration’s
150
lutions, for every Q, we consider 50
i=1 Γi ∼ CN (0, Q IM ). The
choice of the complex Gaussian distribution for the error matrices
Γ is due to the central limit theorem and the assumptions made in
section 4.1.1.
We proceed to calculate the clustered calibration’s MCRLB,
given by Eq. (42), and E{SINRc }, utilizing the Monte-Carlo
method. Jacobian matrices for MCRLB are calculated numerically
and in computation of E{SINRc }, signal power of every cluster is
obtained only using the cluster’s brightest and faintest sources. The

3
0

10

20
30
Number of clusters

40

50

Figure 7. Clustered calibration’s SINR. SINR is low for small Q, when
the interference is large. By increasing the number of clusters the SINR
increases and gets its highest level for which the best performance of the
calibration is expected. After that, it decreases due to the dominant effect of
the background noise, and converges to the un-clustered calibration SINR.

estimated results of MCRLB and E{SINRc } are presented by Fig.
6 and Fig. 7, respectively. As we can see in Fig. 6, for very small
Q, where the effect of interference is large, MCRLB is high. By increasing the number of clusters, MCRLB decreases and reaches its
minimum where the best performance of the clustered calibration is
expected. After that, due to the dominant effect of the background
noise, MCRLB starts to increase until it reaches the CRLB of unclustered calibration. The same result is derived from E{SINRc }
plot of Fig. 7. As Fig. 7 shows, E{SINRc } is low for very small Q,
when the interference (i.e., the error due to clustering) is large. By
increasing the number of clusters, E{SINRc } increases and gets its
peak for which the clustered calibration performs the best. After
that, it decreases and converges to the E{SINR} of un-clustered
calibration.
c 2007 RAS, MNRAS 000, 1–18

Clustered Radio Interferometric Calibration
4.2 Computational cost
In the measurement equation of un-clustered calibration, presented
in Eq. (6), we have M = 2N (N − 1) constraints given by the
visibility vector y, and need to solve for P = 4KN unknown
parameters θ . If P > M , then Eq. (6) will be an under-determined
non-linear system. This clarifies the need of having a small enough
N (number of antennas) and a large enough K (number of sources)
for estimating θ . However, clustered calibration, Eq. (13), has the
advantage of decreasing the number of directions, K, relative to the
number of source clusters, Q ≪ K. This considerably cuts down
the number of unknown parameters P that needs to be calibrated,
thus reducing the computational cost of calibration.

5 SELECTION OF NUMBER OF CLUSTERS
Consider a clustered calibration procedure with a predefined clustering scheme. There is no guarantee that the calibration results
for Q number of clusters, where Q ∈ {1, 2, . . . , K} is randomly
chosen, is the most accurate. Thus, we seek the optimum number
of clusters at which the clustered calibration performs the best. In
this section, we describe the use of: (i) Akaike’s Information Criterion (AIC) (Akaike 1973), as well as (ii) Likelihood-Ratio Test
(LRT) (Graves 1978), in finding this optimum Q for a given observation. Some other alternative criteria could also be found in
Wax & Kailath (1985). Note that for different clustering schemes
the optimum Q is not necessarily the same.

We utilize Akaike’s Information Criterion (AIC) to find the optimum Q for clustered calibration.
e ∼ CN (0, σ
Consider having n
e 2 IM ) in the general data model of
clustered calibration, Eq. (13). Then, the log-likelihood of the visibility vector y is given by
=
−

−M log π − M log σ
e2

Q
Q
X
X
1
e))H (y −
e
e
(y
−
s
(
θ
si (θe)).
i
σ
e2
i=1
i=1

1
c
(y −
σ
e2 =
M

i=1

e
si (θe))H (y −

Q
X
i=1

=

−M log π − M

−

M log{

e
si (θe)).

(46)

Q
Q
X
X
1
e
e
(y −
si (θe))H (y −
si (θe))}. (47)
M
i=1
i=1

Using Eq. (47), the AIC is given by

b
AIC(Q) = −2L(θe|y) + 2(2Pe).

(48)

The optimum Q is selected as the one that minimizes AIC(Q).
5.2 Likelihood-Ratio Test (LRT)
Errors in clustered calibration originate from the system (sky and
instrumental) noise, “clustering errors” introduced in section 4.1.1,
c 2007 RAS, MNRAS 000, 1–18

θepi = [Re(Je11,p ) Im(Je11,p ) . . . Re(Je22,p ) Im(Je22,p )]Ti . (50)
T
T
zpkl = [θepk θepl ]T .

(51)

In fact, we are concatenating the solutions of the same antenna for
two different directions (clusters) together.
We define the null H0 model as

(45)

Substituting Eq. (46) in Eq. (45), we arrive at the maximum likelihood estimation of θe,
b
L(θe|y)

i

Then, the parameter vector θepi (corresponding to the i-th direction
and p-th antenna) is obtained by

H0 : zpkl ∼ N (m, Σ0 ).

(52)

m = [m̄(θepk )T m̄(θepl )T ]T ,

(53)

where

The maximum likelihood estimation of the noise variance σ
e 2 is
Q
X

and “solver noise” which is referred to as errors produced by the
calibration algorithm itself. We assume that the true Jones matrices
along different directions (clusters) at the same antenna are statistically uncorrelated. Note that this assumption is only made for
the LRT to produce reasonable results and this assumption is not
needed for clustered calibration to work. Therefore, if such correlations exist, they are caused by the aforementioned errors. Consequently, the more accurate the clustered calibration solutions are,
the smaller their statistical similarities would be. Based on this general statement, the best number of clusters in a clustered calibration
procedure is the one which provides us with the minimum correlations in the calibrated solutions. Note that for a fixed measurement,
the correlation due to the system noise is fixed. Therefore, differences in the statistical similarities of solutions obtained by different
clustering schemes are only due to “clustering errors” and “solver
noise”.
To investigate the statistical interaction between the gain solutions we apply the Likelihood-Ratio Test (LRT).
e pi (θe) for directions i,
Consider the clustered calibration solution J
i ∈ {1, 2, . . . , Q}, at antennas p, where p ∈ {1, 2, . . . , N },
"
#
Je11,p Je12,p
e
.
(49)
Jpi =
Je21,p Je22,p

Let us define for each antenna p and each pair of directions k and
l, where k and l are belong to {1, 2, ..., Q}, a vector zpkl as

5.1 Akaike’s Information Criterion (AIC)

L(θe|y)

9

and

Σ0 =

"

s2 (θepk )
0

0
s2 (θepl )

#

.

(54)

In Eq. (53) and Eq. (54), m̄ and s2 are denoting sample mean and
sample variance, respectively. Note that having a large number of
samples in hand, the assumption of having a Gaussian distribution
for solutions is justified according to the Central Limit theorem.
The structure of the variance matrix Σ0 tells us that the statistical
correlation between the components of the random vector zqkl , or
between the solutions θepk and θepl , is zero. This is the ideal case in
which there are no estimation errors.
To investigate the validity of the null model compared with the
case in which there exists some correlation between the solutions,
we define the alternative H1 model as
H1 : zpkl ∼ N (m, Σ1 ),
where the variance matrix Σ1 is given by


s2 (θepk )
Cov(θepk , θepl )
.
Σ1 = 
Cov(θepk , θepl )T
s2 (θepl )

(55)

(56)

10

Kazemi et al.

Cov(θepk , θepl ) in Eq. (56) is the 8 × 8 sample covariance matrix.
Using the above models, the Likelihood-Ratio is defined as


Likelihood for null model
Λ = −2ln
,
(57)
Likelihood for alternative model
which has a χ2 distribution with 64 degrees of freedom. As Λ becomes smaller, the null model, in which the statistical correlation
between the solutions is zero, becomes more acceptable rather than
the alternative model. Therefore, the smaller the Λ is, the less the
clustered calibration’s errors are, and vice-versa.

6 SIMULATION STUDIES
We use simulations to compare the performance of un-clustered
and clustered calibration. Working with simulations has the advantage of having the true solutions available, which is not the case in
real observations. That makes the comparison much more objective. Nevertheless, the better performance of the clustered calibration in comparison with the un-clustered ones in calibrating for real
observations of LOFAR is also shown by Kazemi et al. (2011a);
Yatawatta et al. (submitted 2012).
We simulate an East-West radio synthesis array including 14
antennas (similar to WSRT) and an 8 by 8 degrees sky with fifty
sources with low intensities, below 3 Jy. The source positions and
their brightness follow uniform and Rayleigh distributions, respectively. The single channel simulated observation at 355 MHz is
shown in Fig. 8.
We proceed to add gain errors, multiplying source coherencies
by the Jones matrices, as it is shown in Eq. (2), to our simulation.
The amplitude and phase of the Jones matrices’ elements are generated using linear combination of sine and cosine functions. We aim
at simulating a sky with almost uniform variations on small angular
scales. In other words, we provide very similar Jones matrices for
sources with small angular separations. To accomplish this goal, for
every antenna, we first choose a single direction as a reference and
simulate its Jones matrix as it is explained before. Then, for the
remaining forty nine sources, at that antenna, the Jones matrices
(amplitude and phase terms) are that initial Jones matrix multiplied
by the inverse of their corresponding angular distances from that
reference direction. The result of adding such gain errors to our
simulation is shown in Fig. 9.

6.1 Performance comparison of the Clustered and
un-clustered calibrations at SNR=2
We add noise n ∼ CN (0, σ 2 I) with σ 2 = 28, as it is shown in Eq.
(6), to our simulation. The result has a SN R = 2 and is presented
in Fig. 10. We have chosen to present the case of SN R = 2 since
for this particular simulated observation both the divisive hierarchical and the weighted K-means clustered calibrations achieve their
best performances at the same number of clusters, as will be shown
later in this section.
We apply un-clustered and clustered calibration on the simulation to compare their efficiencies. The fifty sources are grouped
into Q ∈ {3, 4, . . . , 49} number of clusters, using the proposed
divisive hierarchical and weighted K-means clustering algorithms.
Self-calibration is implemented via Space Alternating Generalized Expectation Maximization (SAGE) algorithm (Fessler & Hero

1994; Yatawatta et al. 2009; Kazemi et al. 2011b) with nine iterations. Plots of the averaged Frobenius distance between the simulated (true) Jones matrices and the obtained solutions is shown in
Fig. 11. As we can see in Fig. 11, for both clustering schemes, increasing the number of clusters decreases this distance and the minimum is reached at approximately thirty three clusters (Q = 33).
Beyond this number of clusters, it increases until the fifty individual
sources become individual clusters. This shows that the best performance of both the divisive hierarchical and the weighted K-means
clustered calibrations is approximately at thirty three clusters and
is superior to that of the un-clustered calibration.
The Frobenius distance curves in Fig. 11, the MCRLB curve in
Fig. 6, and the SINR curve in Fig. 7 illustrate that clustered calibration with an extremely low number of clusters does not necessarily
perform better than the un-clustered calibration. The reason is that
when there are only a small number of clusters, the interference, or
the so-called “clustering errors” introduced in section 4.1.1, is relatively large. Therefore, the effect of this interference dominates the
clustering of signals. On the other hand, we reach the theoretical
performance limit approximately after twenty five number of clusters and therefore increasing the number beyond this point gives
highly variable results, mainly because we are limited by the number of constraints as opposed to the number of parameters that we
need to solve for. But, this is not the case for the plots in Fig. 6 and
Fig. 7. The reason is that the MCRLB results of Fig. 6 as well as the
SINR results of Fig. 7 are obtained by Monte-Carlo method with
iterations over fifty different sky and noise realizations. However,
Fig. 11 is limited to the presented specific simulation with only one
sky and one noise realization.
The residual images of the un-clustered calibration as well as
the divisive hierarchical and weighted K-means clustered calibrations for Q = 33 are shown by Fig. 12 and Fig. 13, respectively.
As it is shown by Fig. 12, in the result of un-clustered calibration,
the sources are almost not subtracted at all and there is a significant residual error remaining. The residuals have asymmetric Gaussian distribution with variance σ 2 = 82.29 which is much larger
than the simulated (true) noise variance σ 2 = 10.85. On the other
hand, the sources have been perfectly subtracted in the case of clustered calibration, Fig. 13, and the residuals converge to the simulated background noise distribution. The residuals of the divisive
hierarchical and Weighted K-means clustered calibrations follow
symmetric zero mean Gaussian distributions with σ 2 = 20.17 and
σ 2 = 18.76, respectively. These variances are closer to the simulated one σ 2 = 10.85 and this indicates the promising performance
of clustered calibration. As we can see, hierarchical clustered calibration provides a slightly better result compared to the K-means
one. This is due to the fact that hierarchical clustering constructs
clusters of smaller angular diameters and thus it assigns the same
calibration solutions to the sources with smaller angular separations.
We also calculate the Root Mean Squared Error of Prediction
(RMSEP) to assess the performance of clustered and un-clustered
calibrations’ non-linear regressions. The results of log(RMSEP),
presented by Fig. 14, also justify that the best efficiencies of the hierarchical and K-means clustered calibrations are obtained at thirty
three number of clusters. But, note that there is a difference between the behavior of log(RMSEP) plot of Fig. 14 and the plots
of MCRLB, SINR, and Frobenius distance between the simulated
Jones matrices and solutions in Fig. 6, Fig. 7, and Fig. 11, respectively. In Fig. 14, log(RMSEP) of clustered calibration is less than
that of un-clustered calibration, even for extremely low number of
clusters. This means that even with those low number of clusters,
c 2007 RAS, MNRAS 000, 1–18

Clustered Radio Interferometric Calibration

11

Jy
1.4
1.2
1
0.8
0.6
0.4
0.2
0
−0.2

Figure 8. Single channel simulated observation of fifty sources, with intensities below 3 Jy. The source positions and their brightness are following uniform
and Rayleigh distributions, respectively. The image size is 8 by 8 degrees at 355 MHz. There are no gain errors and noise in the simulation.

Jy
0.6
0.5
0.4
0.3
0.2
0.1
0
−0.1

Figure 9. Simulated image with added gain errors. The errors, the complex 2 × 2 Jones matrices, are generated as linear combinations of sin and cos functions.
The variation of the sky is almost uniform on small angular scales.

c 2007 RAS, MNRAS 000, 1–18

12

Kazemi et al.
Jy
0.6
0.5
0.4
0.3
0.2
0.1
0
−0.1
−0.2

Figure 10. Simulated image of sky, corrupted by gain errors and by noise, as in Eq. (6). The simulated noise vector n has zero mean complex Gaussian
distribution and the SNR is equal to 2.

Jy

0.1
0.08
0.06
0.04
0.02
0
−0.02
−0.04
−0.06
−0.08
−0.1
Figure 12. Residual image of the un-clustered SAGE calibration for fifty sources. The sources are almost not subtracted at all and there are significant residual
errors around them. The residuals have asymmetric Gaussian distribution with variance σ2 = 82.29 which is much larger than the true noise variance
σ2 = 10.85.

c 2007 RAS, MNRAS 000, 1–18

Clustered Radio Interferometric Calibration

13

Jy

0.1
0.08
0.06
0.04
0.02
0
−0.02
−0.04
−0.06
−0.08
−0.1

||True Jones−Calibration solutions||

Figure 13. The residual images of the clustered calibration using hierarchical (right) and Weighted K-means (left) clustering methods with thirty three
source clusters. Calibration is implemented by SAGE algorithm. The sources are subtracted perfectly and the residuals converge to the simulated background
noise distribution. The hierarchical clustering and Weighted K-means residuals follow symmetric zero mean Gaussian distributions with σ2 = 20.17 and
σ2 = 18.76, respectively, where the simulated noise distribution is, CN (0, 10.85I).

we have

35
Weighted K−means
Divisive Hierarchical
un−clustered

30

y=

i=1

25
20

e
si (θe) +

K
X

Γ i + n,

(58)

i=1

After executing a calibration for the above data model, there is a
distance between the target parameters θe and the estimated solub
tions θe. This is the so-called “solver noise”, mentioned in section

15

5.2. Thus, the residuals are given as

10
5
0
0

Q
X

10

20
30
Number of Clusters

40

50

Figure 11. The average Frobenius distance between the simulated (true)
Jones matrices and the solutions of clustered and un-clustered calibrations.
The two curves represent clustered calibration via divisive hierarchical and
weighted K-means clustering algorithms. By increasing the number of clusters, for both clustering methods, this distance is decreased and gets its minimum approximately at thirty three clusters. After that, it is increased till the
fifty individual sources. That shows that the best performance of the clustered calibration is at around thirty three clusters.

clustered calibration performs better than the un-clustered calibration. This is somewhat in disagreement with the scenarios shown in
Fig. 6, Fig. 7, and Fig. 11. For a better understanding of the reason
behind this contrast, first lets see how the residual errors are originated.
Based on Eq. (13) and Eq. (40), in the clustered calibration strategy,
c 2007 RAS, MNRAS 000, 1–18

y−

Q
X
i=1

Q

K

i=1

i=1

X
X
b
b
e
Γ i + n.
{e
si (θe) − e
si (θe)} +
si (θe) =

(59)

From Eq. (59), we immediately see that the background noise n is
fixed
Pand the “clustering errors” are calculated for all the sources
as K
i=1 Γ i . However, since we solve only for Q directions and
not for all the K sources individually, the “solver noise” part,
PQ
b
{e
si (θe) − e
si (θe)}, is also calculated only for Q clusters and
i=1

not for all the K sources. It is clear that for a very small Q, this
term could be much less than for Q ≃ K. Therefore, in Fig. 11, the
result of RMSEP at a very low number of clusters is still less than
the ones at Q = K. When Q is not at the two extremes of being
very small or very large (almost equal to the number of individual
sources), then the result of RMSEP is promising. Moreover, applying more accurate calibration methods or increasing the number of
iterations, the “solver noise” will decrease and subsequently, we
expect the RMSEP curve to have the same behavior as the curves
of Fig. 7 and Fig. 11.

14

Kazemi et al.
6

−0.5 x 10

3
Weighted K−means
Divisive Hierarchical
un−clustered

Weighted K−means
Divisive Hierarchical
un−clustered

−1

1
−1.5
AIC

log(RMSEP)

2

0

−2
−1
−2.5

−2
−3
0

10

20
30
Number of clusters

40

−3
0

50

Figure 14. The RMSEP for clustered and un-clustered calibrations. The results are obtained using a base ten logarithmic scale. The two curves are corresponding to clustered calibration via divisive hierarchical and weighted
K-means clustering algorithms. By increasing the number of clusters, the
results are decreased and the minimum result is obtained at around thirty
three clusters. After that, the results are increased till the fifty individual
sources. That shows the superior performance of the clustered calibration
compared to the un-clustered one. The best performance of the clustered
calibration for both of the applied clustering methods is at around thirty
three number of clusters.

40

50

1800

We utilize AIC and LRT to select the optimum number of clusters for which clustered calibration achieves its best performance.
The methods are applied to our simulation for the case of SNR=2.
The AIC and LRT results are shown in Fig. 15 and Fig. 16, respectively. They both agree on Q = 33 as the optimum number of
clusters for the divisive hierarchical and Weighted K-means clustered calibrations. Likelihood-Ratio plot of Fig. 16 has almost the
same behavior as the plot of Frobenius distance between the simulated Jones matrices and the obtained solutions presented by Fig.
11. The reason is that the results of the both plots are obtained using
the solutions themselves as the input data. However, since AIC results are computed using the residual errors as inputs, which is also
the case for obtaining the RMSEP curves of Fig. 14, AIC curves of
Fig. 15 are slightly steeper than the Frobenius distance between the
simulated Jones matrices and solutions and the Likelihood-Ratio
curves of Fig. 11 and Fig. 16, respectively.

1400

Likelihood−Ratio

1600

We start changing the noise in our simulation to see how it effects the clustered calibration’s efficiency. We simulate the cases
for which SNR∈ {1, 2, . . . , 15} and apply clustered calibration on
them. Since the sky model does not change, the clusters obtained by
divisive hierarchical and weighted K-means methods for the case of
SNR= 2 remain the same. Fig. 17 shows the optimum number of
clusters, on which the best performances of clustered calibrations
are obtained for those different SNRs. As we can see in Fig. 17,
for low SNRs, the optimum Q is small. By increasing the SNR, the
optimum Q is increased till it becomes equal to the number of all
the individual sources that we have in the sky, i.e., K. This means
that when the SNR is very low, the benefit of improving signals
by clustering sources is much higher than the payoff of introducing
“clustering errors” in a clustered calibration procedure. Therefore,
clustered calibration for Q ≪ K has a superior performance com-

20
30
Number of clusters

Figure 15. AIC plot for clustered and un-clustered calibrations. Both the
weighted K-means and divisive hierarchical clustered calibrations get their
minimum AIC at about thirty three clusters. This illustrates that their best
performances are obtained at this number of clusters. Also, their AIC results at thirty three clusters is lower than the un-clustered calibration’s AIC,
which shows their better performances compared to the un-clustered calibration.

6.2 Optimum number of clusters for SNR=2

6.3 Clustered calibration’s efficiency at different SNRs

10

Weighted K−means
Divisive Hierarchical
un−clustered

1200
1000
800
600
400
200
0

10

20
30
Number of clusters

40

50

Figure 16. Likelihood-ratio of the gain solutions obtained by clustered and
un-clustered calibrations. In the both cases of weighted K-means and divisive hierarchical clustered calibrations, the minimum Likelihood-ratio values belong to approximately thirty three number of clusters. These minimums are also lower than the un-clustered calibration’s Likelihood-ratio
result. Therefore, clustered calibration via both the clustering methods performs better than the un-clustered calibration and it achieves the best accuracy in its solutions at thirty three clusters.

pared to the un-clustered calibration. While for a high enough SNR,
the situation becomes the opposite. In this case, the un-clustered
calibration achieves better results compared to clustered calibration
having the disadvantage of introducing “clustering errors”.

6.3.1 Empirical estimation of SINR
Having the results of Fig. 17 in hand, we could find an empirical
model for estimating the optimum number of clusters for various
SNRs. Note that by changing the observation and the instrument
characteristics, this model will also be changed.
As it is explained in section 4.1.3, the best performance of
clustered calibration is obtained when the SINR is at its highest
c 2007 RAS, MNRAS 000, 1–18

Divisive Hierarchical optimum Q

55
50
45
40
35
30
25

Weighted K−means
Divisive Hierarchical
un−clustered

20
15
0

5

SNR

10

Figure 17. The optimum number of clusters, on which the best performance
of clustered calibration is obtained, at different SNRs. For low SNRs, the
efficiency of the clustered calibration is superior to the un-clustered calibration. As the SNR gets higher, clustered calibrations achieve their best
solutions utilizing a higher number of clusters. Finally, when the SNR is
high enough, the performance of un-clustered calibration becomes better
than the clustered one.

level. Fig. 17 shows the number of clusters on which the maximum
SINR was obtained, where the signal (sky) and the noise powers are
known a priori. Thus, the only unknown for estimating the SINR
is the interference, or the “clustering errors”, for which we need
to have a prediction model. After estimating the SINR using this
model, finding the optimum Q will be straightforward.
Consider the definition of “clustering errors” given by Eq.
(31). It is logical that for every source, the difference between its
e is
true Jones matrix and the clustered calibration solution, ||J − J||,
a function of the angular distance between the source and the centroid of the cluster that it belongs to. Based on this and using Eq.
(31) and Eq. (34), for the interference of the i-th cluster at baseline
p − q we assume that
X
e i{pq} ||2 {D(Li )}ν I2 ),
Γl{pq} ∼ CN (0, η||C
(60)
l∈Li

where η and ν are unknowns. Eq. (60), in fact, considers an intere i{pq} ||2 {D(Li )}ν for every i-th
ference power (variance) of η||C
cluster, i ∈ {1, 2, . . . , Q}, at baseline p − q. Assuming the interferences of different clusters to be statistically independent from each
other, and bearing in mind that the baseline’s additive noise Npq
has also a complex Gaussian distribution independent from those
interferences’, then the noise plus interference power for the i-th
cluster at baseline p − q is obtained by
e i{pq} ||2 {D(Li )}ν + ||Npq ||2 .
η||C

(61)

Fitting suitable η and ν to Eq. (61), the SINR for the i-th cluster at
e i{pq} ||2 ,
baseline p − q is equal to the cluster’s signal power, ||C
divided by the result of Eq. (61). Subsequently, estimation of
E{SINRc } will be straight forward where the expectation is calculated with respect to all the source clusters and all the baselines. Note that simulation provides us with the true noise power,
||Npq ||2 . In the case of having a real observation, this power could
be estimated by Eq. (46).
Fig. 18 shows the number of clusters on which divisive hierarchical and weighted K-means clustered calibrations achieve their
maximum estimated E{SINRc }. The results are calculated for
SNR∈ {1, 2, . . . , 15}. For the hierarchical clustering η = 1550
c 2007 RAS, MNRAS 000, 1–18

15

50

40

30

20
True value
SINR estimation
10
0

15

Weighted K−means optimum Q

Optimum number of clusters (Q)

Clustered Radio Interferometric Calibration

5

SNR

10

15

50

40

30

20
True value
SINR estimation
10
0

5

SNR

10

15

Figure 18. Optimum number of clusters at which the divisive hierarchical (top) and Weighted K-means (bottom) clustered calibrations perform
the best. For both of the clustering methods, the results obtained by SINR
estimations mostly match the true optimum number of clusters.

and ν = 0.3, and for the K-means clustering η = 2500 and
ν = 0.003. As we can see, for both the clustering methods, these
maximum E{SINRc }s are mostly obtained at the true optimum
number of clusters for which clustered calibration performed the
best. Introducing more refined models compared to Eq. (60) could
even improve the current result.

6.4 Realistic sky models
So far, we have limited our studies to sky models in which the
brightness and position of the radio sources follow Rayleigh and
uniform distributions, respectively. These characteristics provide us
with a smooth and uniform variation of flux intensities in our simulated skies. In such a case, the effects of the background noise on
the faintest and the strongest signals are almost the same. Therefore, if clustered calibration performs better than the un-clustered
calibration that would be only based on upgrading the signals
against the noise. Although, in nature, we mostly deal with the sky
models in which the distribution of the flux intensities is a power
law, with a steep slope, and the spatial distribution is Poisson.
Hence, there exist a few number of very bright sources, whose signals are considerably stronger than the others, and they are sparse in
the field of view. The corruptions of the background noise plus the
interferences of the strong signals of those few bright sources make

Kazemi et al.

the calibration of the other faint point sources impractical. Thus,
there is the need for utilizing the clustered calibration which applies
the solutions of the bright sources to their closed by fainter ones or
solves for upgraded signals obtained by adding up a group of faint
signals together. This has been shown by Kazemi et al. (2011a);
Yatawatta et al. (submitted 2012), when comparing the efficiency
of the clustered and un-clustered calibrations on LOFAR real observations. In this section, using simulations, we also reveal the
superiority of clustered calibration compared to the un-clustered
calibration for such sky models.
We simulate a sky of 52 radio point sources which are obtained
by modified Jelić et al. (2008) foreground model. The brightness
distribution of the point sources follows the source count function
obtained at 151 MHz (Willott et al. 2001), while the angular clustering of the sources are characterized by a typical two-point correlation function,
ρ(d) = Ad−0.8 .

(62)

In (62), ρ is the two point correlation function, d is the angular separation, and A is the normalization amplitude of ρ. The flux cut off
is 0.1 Jy.
We corrupt the signals with gain errors which are linear combinations of sin and cos functions, as in our previous simulations. At
the end, a zero mean Gaussian thermal noise with a standard deviation of 3 mJy is added to the simulated data. The result is shown
in Fig. 19. In Fig. 19, all the bright sources are concentrated on the
right side of the image, rather than being uniformly distributed in
the field of view, and the rest of the sources are so faint that are
almost invisible.
We apply the clustered and un-clustered calibrations on Q ∈
{3, 4, . . . 51} number of clusters and K = 52 number of individual sources, respectively. The clustering method used is the divisive hierarchical and the calibrations are executed via SAGE algorithm with nine number of iterations. The residual noise variances obtained are demonstrated in Fig. 20. As Fig. 20 shows, the
level of the residual noise obtained by the clustered calibration for
Q ∈ {15, 16, . . . 45} number of clusters is always below the result
of the un-clustered calibration. This proves the better performance
of the clustered calibration. The best result of the clustered calibration, with the minimum noise level, is achieved for Q = 27 number
of clusters.
The residual images of the clustered calibration with Q = 27
number of source clusters, and the un-clustered calibration for
K = 52 individual sources are shown by Fig. 21. In the right side
of the residual image of the un-clustered calibration there exist artificial stripes caused by over and under estimating the brightest
sources of the field of view. That shows the problematic performance of the un-clustered calibration. However, clustered calibration has generated much less artificial effects after subtracting these
sources. On top of that, the zoomed in window in the left side of
the images of Fig. 21 show that the faint sources are not removed
by the un-clustered calibration at all, while being almost perfectly
subtracted by the clustered calibration. Moreover, the residual noise
of the clustered calibration follows a symmetric zero mean Gaussian distributions with a standard deviation of 4.2 mJy, while the
one from the un-clustered calibration has an asymmetric Gaussian
distribution with mean and standard deviation equal to -1.2 and 5.3
mJy, respectively. Taking to account that the simulated noise distribution is a zero mean Gaussian distribution with a variance of
3 mJy, the superior performance of the clustered calibration compared to the un-clustered one is evident.
As the final conclusion of this simulation, calibrating below

8
Divisive Hierarchical
Un−clustered
Noise variance

16

7

6

5

4
0

10

20
30
Number of clusters

40

50

Figure 20. The noise variances of the residual images, obtained by clustered and un-clustered calibrations, in mJy. The level of the residual noise
obtained by the clustered calibration for Q ∈ {15, 16, . . . 45} number of
clusters stands below the result of the un-clustered calibration. That reveals
the superior performance of the clustered calibration in comparison with the
un-clustered one. The best result of the clustered calibration which achieves
the minimum noise level is at Q = 27 number of clusters.

the noise level, clustered calibration always performs better than
the un-clustered calibration. This is regardless of the sky model and
is only based on the fact that solving for individual sources with
very poor signals is impractical. Nevertheless, when some sources
are very close to each other, the sky corruption on their signals
would be exactly the same and there is no point in solving for every
of them individually.

7 CONCLUSIONS
In this paper, we demonstrate the superior performance of “clustered calibration” compared to un-clustered calibration especially
in calibrating sources that are below the calibration noise level. The
superiority is in the sense of having more accurate results by the
enhancement of SNR as well as by the improvement of computational efficiency by reducing the number of directions along which
calibration has to be performed.
In a “clustered calibration” procedure, sky sources are
grouped into some clusters and every cluster is calibrated as a single source. That replaces the coherencies of individual sources by
the total coherency of the cluster. Clustered calibration is applied to
these new coherencies that carry a higher level of information compared with the individual ones. Thus, for the calibration of sources
below the noise level it has a considerably better performance compared to un-clustered calibration. Another way of looking at clustering is to consider the distribution of source flux densities, i.e.
the number of sources vs. the source flux density curve. Regardless
of the intrinsic flux density distribution, clustering makes the number of clusters vs. the cluster flux density curve more uniform, thus
yielding superior performance. An analytical proof of this superiority, for an arbitrary sky model, is presented using MCRLB and
SINR analysis.
KLD and LRT are utilized to detect the optimum number of
clusters, for which the clustered calibration accomplishes its best
performance. A model for estimating SINR of clustered calibration
is also presented by which we could find the optimum number of
clusters at low computational cost.
c 2007 RAS, MNRAS 000, 1–18

Clustered Radio Interferometric Calibration

17

Jy

0.8

0.6

0.4

0.2

0

−0.2
Figure 19. Simulated observation of fifty two sources which are obtained by modified Jelić et al. (2008) foreground model. The corrupting gain errors are
generated as linear combinations of sin and cos functions. The image size is 8 by 8 degrees and the additive thermal noise is a zero mean Gaussian noise with
a standard deviation of 3 mJy.

Divisive hierarchical as well as Weighted K-means clustering
methods are used to exploit the spatial proximity of the sources.
Simulation studies reveal clustered calibration’s improved performance at a low SNR, utilizing these clustering algorithms. Both the
clustering methods are hard clustering techniques which divide data
to distinct clusters. However, we expect more accurate results using
fuzzy (soft) clustering, which constructs overlapping clusters with
uncertain boundaries. Application and performance of this type of
clustering for clustered calibration will be explored in future work.

ACKNOWLEDGMENT
The first author would like to gratefully acknowledge NWO grant
436040 and to thank N. Monshizadeh for the useful discussions. We
also thank the referee for a careful review and valuable comments
that helped enhance this paper.

REFERENCES
Akaike H., 1973, in Second International Symposium on Information Theory (Tsahkadsor, 1971), Akadémiai Kiadó, Budapest,
pp. 267–281
Bernardi G., Mitchell D. A., Ord S. M., Greenhill L. J., Pindor B.,
Wayth R. B., Wyithe J. S. B., 2010, ArXiv e-prints
Bezdek J. C., 1981, Pattern Recognition with Fuzzy Objective
Function Algorithms. Kluwer Academic Publishers, Norwell,
MA, USA
Born M., Wolf E., 1999, Principles of Optics. Cambridge University Press
Bregman J. D., 2012, PhD Thesis, University of Groningen
c 2007 RAS, MNRAS 000, 1–18

Fessler J. A., Hero A. O., 1994, IEEE Transactions on Signal Processing, 42, 2664
Gini F., Reggiannini R., 2000, IEEE Transactions on Communications, 48, 2120
Gini F., Reggiannini R., Mengali U., 1998, IEEE, 46, 52
Graves S., 1978, British J. Philos. Sci., 29, 1
Hamaker J. P., 2006, A&A, 456, 395
Hamaker J. P., Bregman J. D., Sault R. J., 1996, A&AS, 117, 137
Intema H. T., van der Tol S., Cotton W. D., Cohen A. S., van
Bemmel I. M., Röttgering H. J. A., 2009, A&A, 501, 1185
Jelić V., Zaroubi S., Labropoulos P., Thomas R. M., Bernardi G.,
Brentjens M. A., de Bruyn A. G., Ciardi B., Harker G., Koopmans L. V. E., Pandey V. N., Schaye J., Yatawatta S., 2008, MNRAS, 389, 1319
Johnson S. C., 1967, Psychometrika, 32, 241
Kaufman L., Rousseeuw P., 1990, Finding Groups in Data An Introduction to Cluster Analysis. Wiley Interscience, New York
Kay S., 1993, Fundamentals of Statistical Signal Processing: Estimation theory, Prentice Hall signal processing series. PrenticeHall PTR
Kazemi S., Yatawatta S., Zaroubi, 2012, IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP),
2533
Kazemi S., Yatawatta S., Zaroubi S., 2011a, IEEE Workshop on
StatisticalSignal Processing, SSP2011, 414, 597
Kazemi S., Yatawatta S., Zaroubi S., de Bruyn A. G., Koopmans
L. V. E., Noordam J., 2011b, MNRAS, 414, 1656
Kerdprasop K., Kerdprasop N., Sattayatham P., 2005, in In
DaWaK, pp. 488–497
Liu A., Tegmark M., Zaldarriaga M., 2009, MNRAS, 394, 1575
MacQueen J. B., 1967, in Proc. of the fifth Berkeley Symposium
on Mathematical Statistics and Probability, Cam L. M. L., Ney-

18

Kazemi et al.

Jy
0.03
0.025
0.02
0.015
0.01
0.005
0
−0.005
−0.01
−0.015
Figure 21. The residual images of the clustered calibration for Q = 27 number of clusters (right) and the un-clustered calibration for K = 52 (left).
Calibration is implemented by SAGE algorithm with nine number of iterations and the clustering method applied is divisive hierarchical. In the right side of
the residual image of the un-clustered calibration we see stripes of over and under estimations due to problematic performance of the un-clustered calibration
in subtracting the brightest sources. The zoomed in window in the left side of the image also shows that the faint sources are not removed at all. However,
clustered calibration could remove all the faint sources almost perfectly and has generated much less artefacts after subtracting for the brightest sources in the
right side of the field of view. Moreover, the residual noise of the clustered calibration follows a symmetric zero mean Gaussian distributions with a standard
deviation of 4.2 mJy, while the one from the un-clustered calibration has an asymmetric Gaussian distribution with mean and standard deviation equal to -1.2
and 5.3 mJy, respectively. Taking to account that the simulated noise distribution is a zero mean Gaussian distribution with a standard deviation of 3 mJy, the
better performance of the clustered calibration compared to the un-clustered one is evident.

man J., eds., Vol. 1, University of California Press, pp. 281–297
Pearson T. J., Readhead A. C. S., 1984, ARA&A, 22, 97
Pindor B., Wyithe J. S. B., Mitchell D. A., Ord S. M., Wayth R. B.,
Greenhill L. J., 2010, PASA, 28, 46
Smirnov O. M., 2011, A&A, 527, A108
Thompson A. R., Moran J. M., Swenson G. W., 2001, Interferometry and Synthesis in Radio Astronomy, 2nd edn. Wiley-VCH
Wax M., Kailath T., 1985, Acoustics, Speech and Signal Processing, IEEE Transactions on, 33, 387
Willott C. J., Rawlings S., Blundell K. M., Lacy M., Eales S. A.,
2001, MNRAS, 322, 536
Xu R., Wunsch D., 2008, Clustering. John Wiley and Sons
Yatawatta S., de Bruyn A. G., LOFAR Collaboration, submitted
2012, Astronomy and Astrophysics
Yatawatta S., Zaroubi S., de Bruyn G., Koopmans L., Noordam J.,
2009, Digital Signal Processing Workshop and 5th IEEE Signal
Processing Education Workshop, 2009. DSP/SPE 2009. IEEE
13th, 150
This paper has been typeset from a TEX/ LATEX file prepared by the
author.

c 2007 RAS, MNRAS 000, 1–18

