# How to avoid X’es around point sources in maximum likelihood CMB maps

Sigurd K. Næss  Center for Computational Astrophysics, Flatiron Institute

###### Abstract

The maximum likelihood estimator for CMB map-making is optimal and unbiased as
long as the data model is correct, but in practice it rarely is, with model
errors including sub-pixel structure and instrumental problems like time-
variable gain and pointing errors. In the presence of such errors, the
solution is biased, with the local error in each pixel leaking outwards along
the scanning pattern by a noise correlation length. The most important sources
of such leakage are strong point sources, and for common scanning patterns the
leakage manifests as an X around each such source. I discuss why this happens,
and present several old and new methods for mitigating and/or eliminating this
leakage, along with a small stand-alone TOD simulator and map-maker in Python
that implements them.

##  1 Introduction

Cosmic Microwave Background (CMB) telescopes map the sky by scanning an array
of detectors repeatedly across the sky, resulting in a time-series of samples
$d$. This process is usually modeled via the linear equation (Tegmark, 1997)

$\displaystyle d$ $\displaystyle=Pm+n$ (1)

where $d$ is the vector of samples, $m$ is a vector representing the pixels of
the sky we are trying to reconstruct, and $P$ is a response matrix that
encodes where the telescope was pointing on the sky when each sample was taken
and how it responded to that. $n$ represents the noise contribution to each
sample, and is usually taken to be normal distributed with some covariance
$N$. The maximum likelihood solution to this equation is

$\displaystyle\hat{m}$ $\displaystyle=(P^{T}N^{-1}P)^{-1}P^{T}N^{-1}d\equiv
M_{PN}d$ (2)

where $M_{N}$ is the map-making operator using the response matrix $P$ and the
noise model $N$. It’s simple to see that this solution is unbiased by
inserting our model for $d$:

$\displaystyle\langle\hat{m}\rangle$
$\displaystyle=(P^{T}N^{-1}P)^{-1}P^{T}N^{-1}(Pm+\overbrace{\langle
n\rangle}^{0})=m$ (3)

Given this result it might be surprising to hear that most maximum-likelihood
maps made using eq. 2 are biased to some extent, with the archetypical example
being linear artifacts extending away from bright point sources in the map,
often in the shape of an X (but this depends on the telescope’s scanning
pattern). The source of these artifacts is _model error_ : the failure of eq.
1 to accurately describe the real data; and in this paper I will describe the
main types of model error in maximum-likelihood mapmaking, simulate their
effect and investigate several old and new methods for mitigating or
eliminating them.

##  2 Sources of model error

###  2.1 Sub-pixel structure

The most obvious problem with eq. 1 is that it models the sky as a finite
vector of pixels, which will leave out any signal on scales smaller than the
pixel spacing. This is exacerbated by the canonical choice of a _nearest-
neighbor_ response model in $P$. Every CMB map-maker I am aware of uses
nearest-neighbor, including destripers like MADAM (Keihänen et al., 2010) and
SRoll (Planck Collaboration, 2018), filter+bin map-makers like those used in
SPT (Schaffer et al., 2011) and BICEP (BICEP2 Collaboration, 2014) or the
maximum likelihood map-makers used in ACT (ACT Collaboration, 2017) and QUIET
(Ruud et al., 2015). In this model, the value of each sample is simply given
by the value of the closest pixel to it, regardless of how far away from the
pixel center it is. This results in a response matrix with a very simple
structure: For each row (representing a sample) a single element will be 1
(representing the pixel hit by that sample), and all others will be zero.

Multiplication by a matrix with this structure is simple and efficient. Given
a map m[ny,nx] and arrays y[nsamp], x[nsamp] containing the x and y pixel
coordinates of each sample, the forward operation $d=Pm$ can be implemented as

⬇

for i in range(nsamp):

d[i] = m[round(y[i]),round(x[i])]

and the transpose operation $m=P^{T}d$ can be implemented as

⬇

for i in range(nsamp):

m[round(y[i]),round(x[i])] += d[i]

However, this simplicity comes at a cost: The signal model implied by this
response matrix has a constant value inside each pixel and a discontinuous
jump at the edge of each pixel. This model and the resulting error is
illustrated in figure 1.

uncorrelated | correlated
---|---
|
Figure 1: A simulated data vector consisting of a noiseless scan through a
smooth, Gaussian profile, compared to its nearest-neighbor-pixellated maximum
likelihood model. Left: When an uncorrelated noise model is used ($N\propto
I$) the residual is large but stays in the pixel where they were sourced. This
error can be characterized as a simple pixel window. Right: When a correlated
noise model is used the residuals are still large, but are now spread out over
a larger area than the signal that sourced them.

A step-function model like this clearly cannot match the smooth behavior of
the beam-convolved sky, and the only term in eq. 1 that can absorb such a
mismatch between the data and the model is the noise term $n$. How this model
error manifests is therefore determined by our noise model $N$.

If $N$ is assumed to be uncorrelated ($N$ diagonal in time domain), then a
high value of the residual (honorary noise as far as our data model is
concerned) in one sample does not tell us anything about what the noise is
doing in any other samples, and so the residual just stays in the pixel that
birthed it. The value in each pixel is therefore just the mean of the samples
that hit it (left panel of figure 1). This is approximately the mean value of
the signal inside the area covered by the pixel, which is a straightforward
and useful thing to measure, and is not the sort of bias we are worried about
here.

If, on the other hand, $N$ is assumed to be correlated, then the presence of
high “noise” (actually residual) in some samples will make the maximum
likelihood prection for other samples within a noise correlation length have a
similar value for the noise. For long correlation lengths this could affect
samples many pixels away from the source of the residual. This is illustrated
in the right panel of figure 1. Here large sub-pixel residuals in the center
result in a non-zero expectation value for the noise at large distances. Given
that the actual data has negligibly small values there, the best fit model
must cancel the expected noise, resulting in a non-zero model extending far
beyond the area with significant signal. Effectively the correlated noise
model is causing local model errors to leak roughly a correlation length into
the surrounding pixels.

###  2.2 Signal or instrumental variability

Real telescopes don’t have 100% accurate pointing and gain, and even if they
did, the sky is not static and unchanging. These effects mean the telescope
effectively sees the sky jitter around slightly while varying in brightness.
There is no room for this in eq. 1 aside from the noise term $n$, so time-
variable instrumental errors or variable point sources will lead to exactly
the same sort of artifacts as sub-pixel errors do.

##  3 Mitigation methods

The artifacts are ultimately caused by a combination of two factors:

  1. 1.

Multiple samples are mapped onto the same pixel, but don’t agree on what the
signal in that pixel should be.

  2. 2.

The correlated noise model causes a non-local response to such errors.

This points to two approaches for eliminating them: Improve the model so it
matches the data more accurately, or modify the noise model to make it more
local.

###  3.1 Higher-order mapmaking

An obvious way to improve the model is to replace nearest-neibhor
interpolation in the response matrix $P$ with smoother interpolation functions
such as bilinear or bicubic spline interpolation. These are popular in image
processing, but are as far as I am aware unused in CMB map making.

Bilinear interpolation is considerably slower and more complicated than
nearest neighbor, and each sample now interacts with the four closest pixels
instead of just one. The forward operation $d=Pm$ is now (based on the
implementaiton of scipy.ndimage.map_coordinates):

⬇

for i in range(nsamp):

py1, px1 = floor(y[i]), floor(x[i])

py2, px1 = py1+2, px1+1

ry, rx = y[i]-py1, x[i]-px1

vx1 = m[py1,px1] * (1-rx) + m[py1,px2] * rx

vx2 = m[py2,px1] * (1-rx) + m[py2,px2] * rx

d[i] = vx1 * (1-ry) + vx2 * ry

The transpose operation $m=P^{T}d$ is:

⬇

for i in range(nsamp):

py1, px1 = floor(y[i]), floor(x[i])

py2, px1 = py1+2, px1+1

ry, rx = y[i]-py1, x[i]-px1

vx1 = d[i] * (1-ry)

vx2 = d[i] * ry

m[py1,px1] += vx1 * (1-rx)

m[py1,px2] += vx1 * rx

m[py2,px1] += vx2 * (1-rx)

m[py2,px2] += vx2 * rx

Bicubic spline interpolation is another large step up in complexity. See the
appendix for details.

Figure 2 shows how these higher-order interpolation schemes perform for the
same 1D Gaussian simulation we considered in figure 1

Figure 2: Comparison of the maximum-likelihood models for nearest neighbor,
linear and cubic interpolation in the response matrix for uncorrelated and
correlated noise models, for the same data as in figure 1. Linear and cubic
interpolation fit the central part of the profile much better than nearest
neighbor, leading to lower discrepancy between the data and model at high
radius in the correlated case (solid curves). But this situation is reversed
for the uncorrelated noise model (dashed curves). This happens because the
linear and cubic models are inherently non-local, and are willing to sacrifice
accuracy in areas where the signal is very weak in order to get a better match
in areas with strong signal.

Both linear and cubic interpolation fit the peak of the Gaussian much better
than nearest neighbor, resulting in roughly 4 orders of magnitude lower
spurious signal at high radius for the case with a correlated noise model.
However, they only reach this level at large radii. Before this, the spurious
signal is dominated by a new effect that was not present for nearest neighbor
interpolation: Higher order interpolation is inherently non-local, which means
that each pixel is correlated with its neighbors (and neighbors-neighbors for
bicubic). This is responsible for the exponentially decaying ringing pattern
we see in the figure for these interpolation types, even for uncorrelated
noise. Effectively the fit is sacrificing accuracy at large radius in order to
fit the strong, central peak slightly better.

So we see that high interpolation order has smaller total residuals and hence
lower non-local artifacts from the noise correlations, but low-order
interpolation has shorter correlation length for the interpolation functions
themselves, leading to a faster decay of the ringing. For the example in
figure 2, the sweet spot is linear interpolation. Even better locality (at
even higher cost) should be possible with an interpolation function based on
the actual point spread function of the instrument.

Aside from the high cost and non-locality inherent in these interpolators,
they are also limited by only addressing sub-pixel errors. Higher-order map-
making does nothing to remove artifacts from variable instrumental errors or
intrinsic sky variability.

####  3.1.1 Iterative approximation

Since the response matrix often is the bottleneck in maximum-likelihood
mapmaking, it may be too expensive to make this step several times slower by
using higher-order mapmaking. A simple approximation that avoids most of this
slowdown is to first make a normal map using the fast, nearest-neighbor
response matrix, and then subtract this estimate of the sky from the time-
ordered data using a bilinear or bicubic spline response matrix. This process
can be repeated if necessary, though in my tests there was little point in
doing more than 1-3 iterations.

$\displaystyle\hat{m}_{0}$ $\displaystyle=0$ $\displaystyle\hat{m}_{i+1}$
$\displaystyle=\hat{m}_{i}+M_{PN}(d-\bar{P}\hat{m}_{i})$ (4)

Here $P$ is the standard nearest-neighbor response matrix and $\bar{P}$ is the
higher-order one we are approximating. Note that this iteration converges on
an _approximation_ to the higher-order map, not the exact solution we would
have gotten by using $\bar{P}$ in the equation 2. However, as we shall see in
the results section it does a pretty good job.

If one is already doing iterative mapmaking to avoid bias from estimating the
noise model $N$ from the full data $d$ rather than just the noise $n$, then
this iterative approximation to higher-order mapmaking can simply be done as a
part of that, and will incur almost no extra cost.

Alternatively, if one is solving for each map using an iterative method like
preconditioned conjugate gradients (CG), then the small scales that are the
main source of sub-pixel error often converge much faster than the larger
ones. In this case only a few CG steps are needed for all but the last
iteration of equation 4, again resulting in almost no overhead compared to
standard mapmaking.

###  3.2 Point source subtraction

In figure 1 we saw that non-local sub-pixel artifacts could leak far away from
the signal that sourced them, but it also bears noting that their amplitude
was about $10^{4}$ times lower than the source. For a smooth, low-contrast
signal like the cosmic microwave background, each pixel is contaminated by
leakage from other pixels in the vicinity, but since those other pixels have
signal of similar amplitude but random phase, one can expect each pixel to
only be contaminated at the $10^{-4}$ level, which is negligible. Non-local
model errors are therefore only a concern when a strong, compact signal is
close to a region of weaker signal. In CMB map-making the typical case would
be a strong point source like a quasar.

If one knows the location and amplitude of the bright point sources in the
sky, one can simply subtract those point sources in the data before making the
map. This results in a map free of both the point sources (which can be useful
in itself) and the artifacts they would have sourced. If needed, the source
model can then be projected onto the sky using an uncorrelated noise model $W$
(for example $W=I$) and added back to the map.

$\displaystyle\hat{m}_{\textrm{srcfree}}$
$\displaystyle=M_{PN}(d-\sum_{i}A_{i}B_{i})$
$\displaystyle\hat{m}_{\textrm{tot}}$
$\displaystyle=\hat{m}_{\textrm{srcfree}}+M_{PW}\sum_{i}A_{i}B_{i}$ (5)

Here $A_{i}$ is the amplitude of source $i$, and $B_{i}$ is a vector of each
sample’s response to an instrument beam centered on the position of that
source.

This simple approach was used in e.g. Dünner et al. (2013), and succeeds in
getting rid of the artifacts sourced by sub-pixel structre in the point
sources. However, unless one has a point source model that takes into account
time-variable pointing, gain or beams errors, or intrinsic amplitude
variations, then this method does nothing to reduce artifacts from
variability. It is also limiting that it requires some properties of the sky
to already be known, and that it can only be applied to signals for which an
accurate, non-pixellated model is available.

###  3.3 White source mapping

Model errors do not produce non-local artifacts when using an uncorrelated
nosie model, but of course, simply mapping using such a noise model would
produce a horribly stripy map in the presence of realistic correlated noise.
However, we don’t need everything to be uncorrelated to get rid of leakage, we
just need to decouple the bright point sources from those correlations. This
suggests the following approach:

  1. 1.

Make a point source free set of samples $d^{\prime}$ from $d$ by smoothly in-
painting the samples in $d$ that hit the bright sources, ideally using maximum
likelihood in-painting, but the quality of the in-painting only matters for
the optimality of the method - bad in-painting will not cause a bias111 This
method splits the data into two parts,
$d_{\textrm{tot}}=d_{\textrm{inpaint}}+(d_{\textrm{tot}}-d_{\textrm{inpaint}})$,
and maps the first term with a correlated noise model and the last term with
an uncorelated noise model. As long as both map-making methods are unbiased,
the result will be unbiased regardless of the exact value of
$d_{\textrm{inpaint}}$. In practice map-making with a correlated noise model
can be biased due to model error artifacts, as we have seen. But for
reasonable, smooth in-painting methods this bias will be completely
negligible..

  2. 2.

Make a maximum-likelihood map $\hat{m}^{\prime}$ of the in-painted TOD using
the full noise model.

$\displaystyle\hat{m}^{\prime}$ $\displaystyle=M_{PN}d^{\prime}$ (6)

  3. 3.

If the in-painted regions were small compared to the noise correlation length,
then the missing signal $d-d^{\prime}$ from the previous map will have
approximately white noise, so we can map these bothersome samples using a
simple white noise model:

$\displaystyle\Delta m$ $\displaystyle=M_{PW}(d-d^{\prime})$ (7)

  4. 4.

Add these maps to get the total, unbiased sky map:
$\hat{m}=\hat{m}^{\prime}+\Delta m$.

###  3.4 Source sub-sampling

We can completely eliminate model error by having at least one degree of
freedom per sample, but if we did that to every sample, there would be no
averaging down of the noise, and we wouldn’t really get a map, just a
reorganized version of the input data. However, there is nothing wrong with
oversampling like this just around strong point sources and other problematic
regions, as long as these regions are small compared to the noise correlation
length. This results in the model

$\displaystyle d$ $\displaystyle=Pm+Gs+n$
$\displaystyle=\begin{bmatrix}P&G\end{bmatrix}\begin{bmatrix}m&s\end{bmatrix}^{T}+n$
(8)

Here $s$ are the new per-sample degrees of freedom for the samples that hit
bright sources, and $G$ the matrix that maps them 1-to-1 only the
corresponding samples. To be explicit, $a=Gs$ could be implemented as a[mask]
= s, and $b=G^{T}d$ as b = d[mask], where mask is True for samples that hit
bright sources and False otherwise.

One can then solve jointly for the maximum-likelihood values $\hat{m}$ and
$\hat{s}$. This can still be done using equation 2, since equation 8 is still
of the same form as equation 1. $m$ and $s$ are in isolation not well defined:
Some samples are both covered by pixels from $m$ and by special per-sample
degrees of freedom from $m$, so signal can move freely from one to the other
without changing the residuals. This can be avoided by modifying $P$ to have
it skip those samples, but it is not actually necessary, for the total map

$\displaystyle\hat{m}_{\textrm{tot}}$ $\displaystyle=\hat{m}+M_{PW}\hat{s}$
(9)

_is_ well-defined in either case, and is the artifact-free map of the sky we
are after.

An advantage of this method is that it is very similar to a common way to
handle per-sample cuts in maximum-likelihood map-making. There one also solves
for an extra degree of freedom per cut sample. The only difference is that we
here don’t throw that part of the solution away afterwards, but instead
project them back into the final map using a white noise model (see section
3.5).

###  3.5 Source cutting

Perhaps the conceptually simplest way to avoid model errors from small
problematic regions of the sky is to simply exclude samples hitting those
regions from the mapping process. In practice, however, most correlated noise
models are best represented in Fourier space, and fast Fourier transforms
cannot easily handle missing samples. Hence, simply dropping those samples
isn’t straightforward. A trick that effectively achieves the same thing is to
modify the pointing matrix to give the bad samples a dummy degree of freedom
each instead of pointing them at the sky (Dünner et al., 2013),

$\displaystyle d$ $\displaystyle=\bar{P}m+Gs+n$ (10)

Here $G$ has the same definition as in eq. 8, selecting only the sources
hitting the bright objects, and $\bar{P}$ is identical to $P$ except that
those samples are skipped. The solution to this is equivalent to
$\hat{m}_{\textrm{tot}}$ from eq. 9 with the cut region masked. To see this,
notice that replacing $\bar{P}$ with $P$ will not change the solution outside
the cut area because $\bar{P}$ already has enough degrees of freedom to
perfectly match the data there, resulting in zero local residual, and adding
even more degrees of freedom in this area will not change that.

###  3.6 PGLS

Postprocessed Generalized Least Squares (PGLS) is a method developed for
mapmaking with high-S/N, high dynamic range data from the infrared telescope
Herschel, where artifacts are much more pronounced than in CMB maps (Piazzo et
al., 2012). Unlike the other methods discussed here, this method does not try
to prevent artifacts from forming, but instead attempts to measure and
subtract them in a second pass of map-making. We can model a standard maximum-
likelihood map built using eq. 2 as in practice being $\hat{m}=m+a+q$ where
$m$ is the map we would have gotten for a noise-free dataset with an
uncorrelated noise model, $a$ are the unwanted artifacts and $q$ is the pixel-
space noise. The goal of PGLS is build an estimator $\hat{a}$ for the
artifacts in the map, and then subtract this from $\hat{m}$ to produce a clean
map.

We can exploit the fact that $a$ is a locally poor fit to the data to isolate
it from the signal. Projecting back into time domain we get

$\displaystyle r$ $\displaystyle=P\hat{m}-d=Pm-d+Pa+Pq$ (11)

Ideally the term $d-Pm$ would equal our time-domain noise $n$, but as we have
seen it also contains sub-pixel errors $e$ (which is why the artifacts arise
in the first place). Hence, we have

$\displaystyle r$ $\displaystyle=Pa+Pq-n-e$ (12)

The sub-pixel errors average to zero when averaged over a pixel using an
uncorrelated noise model222 This follows from our definition of $m$ above.
$M_{PW}e=M_{PW}(s-PM_{PW}s)=0$, where $s$ is the noise-free sky signal., so if
we could get rid of the correlated noise in the term $Pq-n$ we could recover
$\hat{a}$ by mapping $r$ with an uncorrelated noise model. PGLS does this by
applying a running median filter $f$ to $r$, resulting in the following
algorithm:

$\displaystyle\hat{m}_{0}$ $\displaystyle=M_{PN}d$ $\displaystyle\hat{m}_{i}$
$\displaystyle=\hat{m}_{i-1}-M_{PW}f(P\hat{m}_{i-1}-d)$ (13)

The median filter length is a compromise between the need to remove as much
correlated noise as possible and to retain the artifacts themselves, and must
be tuned to the instrument. The balance between artifacts and noise can also
be adjusted by iterating the scheme, as indicated by the index $i$. Each
iteration recovers more of the artifacts at the cost of letting through more
noise. For the toy example I test in this paper I found a median filter width
of 21 samples and 5 iterations to work well.

###  3.7 XGLS

The reason why local errors in areas of high signal end up leaking into the
surroundings is ultimately that the noise model does not expect these areas to
be more prone to high “noise” than other regions. We can remedy this if we
have an estimate of sub-pixel variance by addin this as an extra term $X$ in
the noise model such that $N_{\textrm{tot}}=N+X$, where
$X_{ij}=\delta_{ij}\sigma_{i}^{2}$ and where $\sigma_{i}$ is the standard
deviation of the sub-pixel signal variability in the pixel hit by sample $i$.
Given this noise model one can produce an optimal, artifact-free map as
$\hat{m}=M_{PN_{\textrm{tot}}}d$. This is the Pixel Noise Generalized Least
Squares (XGLS) estimator (Piazzo, 2017), and like PGLS it was developed for
high S/N, high-dynamic range infrared data from the Herschel space telescope.

In many ways this is the proper way to solve the problem - instead of more or
less ad-hoc workarounds one tackes the problem at its root cause by telling
the noise model about this extra source of variance. This elegance comes at a
cost, though:

  1. 1.

One needs an estimate of the sub-pixel signal variance. If one knew the noise-
free signal $s$, then this could be estimated as
$X=\textrm{diag}(PM_{PW}(s-PM_{PW}s)^{2})$. Here $\textrm{diag}(x)$ indicates
the diagonal matrix with $x$ on the diagoanal, and the squaring operation is
performed sample-wise. For simplicity this is what I do in the simulations in
this paper, but for real data where $s$ is unknown approximations are needed -
see Piazzo (2017) for a discussion.

  2. 2.

While the sample-diagonal $X$ and the (usually) Fourier-diagonal $N$ are both
easy to invert on their own, their sum is not diagonal in any simple basis and
is therefore much harder to invert, typically requiring an iterative method
like conjugate gradients (CG) Eq. 2 itself is also usually solved with CG,
resulting in a solution process where a full sequence of CG steps for the
noise matrix inversion must be completed for every single CG step for the map.
This makes XGLS dramatically slower than any of the other methods I tested.
This is confounded by the outer CG iteration’s sensitivity to inaccuracy in
the inversion in the inner CG, necessitating a very strict convergence
criterion there, further slowing things down. Piazzo (2017) report that they
made this reasonably performant with some tuning, but still found the starting
point of the iteration to have large effect on the solution. This is an
indicator of convergence problems.

###  3.8 Summary of the methods

Table 1 shows a summary of the methods we will investigate. The ideal method
would be blind (no prior knowledge about which parts of the sky have strong
signal), and would handle both sub-pixel errors and inconsistency at low cost,
but none of the methods tick all those boxes, and in particular none of the
blind methods can handle inconsistency.

| std | it.lin | it.cubic | bilin | bicubic | pgls | xgls | srcsub | srcmask | srcwhite | srcsamp
---|---|---|---|---|---|---|---|---|---|---|---
blind | yes | yes | yes | yes | yes | yes | yes | no | no | no | no
subpix | no | yes | yes | yes | yes | yes | yes | yes | yes | yes | yes
incon. | no | no | no | no | no | no | yes | no | yes | yes | yes
holes | no | no | no | no | no | no | no | no | yes | no | no
cost | 1 | 1* | 1* | 4 | 16 | 1 | 500 | 1 | 1 | 1 | 1
compl. | v.low | high | v.high | high | v.high | low | high | low | low | low | low
pixwin | std | small | small | small | small | std | std | std | std | std | std
new | no | yes | yes | yes | yes | no | no | no | no | yes | yes
Table 1: Summary of the methods discussed in this paper. The rows are: blind:
Whether the method can be used on the whole map with no prior knowledge of
which regions have large model errors; subpix: Whether the method targets sub-
pixel errors; incon.: Whether the method targets inconsistency errors like
gain and pointing errors; cost: A rough indication of time cost of the method,
relative to standard map-making. Iterative linear/cubic mapmaking naively have
a cost factor of $N_{\textrm{iter}}$, but in practice this can usually be
avoided. compl.: A rough indication of how complex the method is to implement
(this is somewhat subjective). pixwin: The pixel window the method produces.
Higher-order interpolation methods results in a “smaller” (closer to unity)
pixel window.new: Whether this is a new method or not. Those marked “yes” are,
to my knowledge, not previously published.

##  4 Simulations

I wrote a minimal maximum-likelihood mapmaking library toy.py in about 300
lines of Python to demonstrate the performance of these methods, along with a
driver script examples.py that generates the plots in this paper. They can be
found at https://github.com/amaurea/model_error. The library emphasizes
simplicity and is neither general nor optimized for speed, and only depends on
numpy and scipy (with the exception of bilinear and bicubic map-making which
also depends on pixell. That said, the techniques described here have also
been implemented in the full Atacama Cosmology Telescope (ACT) mapmaker, and
despite the simplicity of these test cases the results are representative of
running a real map-maker on real data.

###  4.1 2D simulations

The main simulations use two sets of 243 x 243 equi-spaced samples in a square
grid, the first ordered vertically and the second horizontally (the ordering
matters for the noise correlation structure). These are mapped onto an 81 x 81
pixel grid so that there are 3 x 3 samples per pixel per set. Using this setup
I simulated two signals:

  1. 1.

A Gaussian point source with a beam standard deviation of 1 pixel and an
amplitude of $2\cdot 10^{4}$ (500 times higher than the color scale used in
the plots in the results section). This will be representative of e.g. a
bright quasar as seen by a high-resolution CMB telesope like ACT.

  2. 2.

A CMB-like field built by simulating a $1/l$ spectrum with an angular
resolution of 1 pixel and an RMS of 1. This will test whether model error is a
problem in the absence of bright sources.

I simulated two types of noise model:

  1. 1.

An uncorrelated noise model $N_{\textrm{uncorr}}=I$, representing the baseline
case where no leakage is expected.

  2. 2.

A correlated noise model
$N_{\textrm{corr}}(f_{1},f_{2})=\sigma^{2}(1+(f_{1}/f_{\textrm{knee}})^{\alpha})\delta_{f_{1},f_{2}}$.
This is diagonal in Fourier space. $\alpha=-4$ represents an atmosphere-like
steeply falling spectrum, meeting a noise floor with RMS $\sigma=1$ at the
knee frequency $f_{\textrm{knee}}=0.02$ in pixel units.

Becuase it is the non-local weighting implied by a correlated noise model that
causes leakage, not the noise itself, one does not actually have to include
the noise in the simulations. To maximize S/N in the figures most of the
simulations are noise-free, but I also include some noisy simulations to show
how effectively each method suppresses noise.

I simulated 3 types of model errors:

  1. 1.

Sub-pixel signal was simulated generating the signals at higher resolution
than the output maps, as described above. This was included in every
simulation.

  2. 2.

Gain errors were simulated by decreasing/increasing the signal amplitude by
0.5% for the vertical/horizontal data set, representing a total gain mismatch
of 1%. This was included in the simulations labeled “gain error”.

  3. 3.

Pointing errors were simulated by offsetting the horizontal dataset by 0.01
pixels diagonally relative to the vertical dataset. This was includes in the
simulations labeled “pt. error”.

###  4.2 1D simulations

I also produced as smaller number of 1-dimensional simulations for figure 1
and figure 2. These used the same noise models and point source signal as the
2-dimensional simulations, but consisted of just a single scan of length 1100
samples which was mapped onto 100 pixels.

##  5 Results

|  standard |  standard |  it.lin |  it.cubic |  bilin |  bicubic |  pgls |  xgls |  srcsub |  srccut |  srcwhite |  srcsamp
---|---|---|---|---|---|---|---|---|---|---|---|---
|  uncorr |  corr |  corr |  corr |  corr |  corr |  corr |  corr |  corr |  corr |  corr |  corr
plain |  |  |  |  |  |  |  |  |  |  |  |
gain error |  |  |  |  |  |  |  |  |  |  |  |
pt. error |  |  |  |  |  |  |  |  |  |  |  |
plain + noise |  |  |  |  |  |  |  |  |  |  |  |
gain error + noise |  |  |  |  |  |  |  |  |  |  |  |
pt. error + noise |  |  |  |  |  |  |  |  |  |  |  |
plain + cmb + noise |  |  |  |  |  |  |  |  |  |  |  |
gain error + cmb + noise |  |  |  |  |  |  |  |  |  |  |  |
pt. error + cmb + noise |  |  |  |  |  |  |  |  |  |  |  |
Figure 3: Model errors around a bright point source for the different methods
discussed in this paper, for a simulation with 3x3 samples uniformly
distributed inside each pixel for each of a vertical and horizontal scanning
pattern. The top third is for a noise free simulation, the second third
includes strongly correlated noise, and the bottom third adds a CMB-like
smooth signal. The rows are: plain: no special conditions; gain error: the
vertical scanning pattern sees a 1% higher signal than the horizontal scanning
pattern; and pt. error: the horizontal scanning pattern sees the signal offset
by 0.01 pixel diagonally. The column headers refer to the different mapmaking
algorithms described in the text. The sub-headers indicate whether a
correlated or uncorrelated noise model was used. The color range is $\pm
10^{-4}$ of the peak source amplitude.

| standard | it.lin | it.cubic | bilin | bicubic | pgls | xgls | srcsub | srccut | srcwhite | srcsamp
---|---|---|---|---|---|---|---|---|---|---|---
near | 0.868 | 0.143 | 0.270 | 0.172 | 1.777 | 0.295 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
far | 3.518 | 1.146 | 0.392 | 0.889 | 0.178 | 0.201 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
near + pt. | 2.544 | 2.016 | 1.965 | 2.058 | 3.614 | 0.355 | 0.000 | 1.909 | 0.000 | 0.000 | 0.000
far + pt. | 3.603 | 1.241 | 1.248 | 1.212 | 1.239 | 0.257 | 0.000 | 1.235 | 0.000 | 0.000 | 0.000
near + gain | 5.668 | 4.593 | 4.928 | 4.605 | 6.425 | 2.065 | 0.000 | 4.799 | 0.000 | 0.000 | 0.000
far + gain | 3.823 | 2.896 | 3.026 | 2.877 | 3.034 | 1.488 | 0.000 | 3.116 | 0.000 | 0.000 | 0.000

Table 2: The maximum absolute value of the leakage from a point source in
units of $10^{-4}$ of the peak amplitude, measured in two distance bins: near:
5-10 pixels from from the center (just outside the point where the beam
becomes negligible), and far: beyond 10 pixels from the center. The leakage is
measured as the difference between a map using the given model and a
correlated noise model and a map using standard nearest neighbor with a white
noise model. For the beam, pixel size and noise model used here, sub-pixel
errors result in O($10^{-4}$) leakage. Higher-order interpolation reduces this
by a factor of 0.5-20 depending on the distance and method (yes, in the worst
case it makes it 2x worse near the source), while the targeted methods
completely eliminate sub-pixel errors. Once gain or pointing errors are added,
only source cutting, white source mapping and source sub-sampling provide any
meaningful improvement (though this will depend on the size of the errors
chosen).

Figure 3 shows maximum likelihood maps for each method for a strong point
source. In the noise free case (top 3rd of the figure) standard nearest-
neighbor mapmaking with an uncorrelated noise model (first column) results in
no leakage as expected, whether in the presence of normal sub-pixel signal,
gain errors or pointing errors. However, in the presence of realistic
correlated noise (mid 3rd of the figure) this method is extremely suboptimal,
resulting in a map completely dominated by large-scale noise.

Standard map-making with a correlated noise model (second column) down-weights
the noise correctly, but causes an X-shaped pattern of leakage around the
source, regardless of whether the noise is actually present. This only gets
worse in the presence of gain and pointing errors. The size of the leakage
relative to the peak source amplitude is shown in table 2, and is typically of
order $3\cdot 10^{-4}$ for this choice of beam, pixel size and noise model.
While this is small relative to the source itself, it might be much brighter
than any surrounding CMB for a sufficiently bright source.

Full and iterative linear and bicubic interpolation all behave similarly. In
the absence of gain and pointing errors they provide a moderate leakage
reduction, typically by a factor of 3-10, but varying by method and distance
from the source. The computationally heaviest method, bicubic interpolation,
provides the best reduction (about a factor of 20) of these at large distance,
but as we saw in figure 2, this comes at the cost of making things worse near
the source. The iterative approximation to bicubic interpolation might be the
best compromise between low-range and short-range leakage, reducing both to
$<0.4\cdot 10^{-4}$. However, since these methods only target sub-pixel
errors, it comes as no surprise that they do not help at all with data
inconsistencies like gain and pointing errors, where the X returns in full
force.

Source subtraction completely eliminates the artifacts in the ideal case, but
performs no better than the blind interpolation methods in the presence of
gain and pointing errors. The other targeted methods fare better: Source
cutting avoids artifacts even for inconsistent data, at the cost of a hole in
the map, while white source mapping and source sub-sampling avoid both these
problems, eliminating artifacts even in the presence of inconsistent data
while still providing a map of the source.

|  standard |  standard |  it.lin |  it.cubic |  bilin |  bicubic |  pgls |  xgls |  srcsub |  srccut |  srcwhite |  srcsamp
---|---|---|---|---|---|---|---|---|---|---|---|---
|  uncorr |  corr |  corr |  corr |  corr |  corr |  corr |  corr |  corr |  corr |  corr |  corr
plain |  |  |  |  |  |  |  |  |  |  |  |
gain error |  |  |  |  |  |  |  |  |  |  |  |
pt. error |  |  |  |  |  |  |  |  |  |  |  |
plain resid |  |  |  |  |  |  |  |  |  |  |  |
gain error resid |  |  |  |  |  |  |  |  |  |  |  |
pt. error resid |  |  |  |  |  |  |  |  |  |  |  |
plain wdiff |  |  |  |  |  |  |  |  |  |  |  |
gain error wdiff |  |  |  |  |  |  |  |  |  |  |  |
pt. error wdiff |  |  |  |  |  |  |  |  |  |  |  |
Figure 4: Like figure 3, but for a noiseless CMB-like signal instead of a
strong point source. Top: the resulting maps. Middle: the root-mean-square of
the residual (data-model) per pixel. Bottom: The difference between each map
and its uncorrelated version. The color range is $\pm 4$ in arbitrary units at
the top and $\pm 0.1$ in the same units in the middle and bottom. The small
gain and pointing errors we simulate here have almost no effect for low
dynamic range fields like the CMB. The source-specific methods assume a zero-
amplitude source in the patch center, even though none is actually present.

| standard | it.lin | it.cubic | bilin | bicubic | pgls | xgls | srcsub | srccut | srcwhite | srcsamp
---|---|---|---|---|---|---|---|---|---|---|---
plain | 36.58 | 21.48 | 12.70 | 5.02 | 0.88 | 150.84 | 36.38 | 36.58 | 79.67 | 274.90 | 36.82
pt. error | 38.68 | 33.55 | 26.29 | 14.35 | 13.78 | 153.88 | 38.47 | 38.68 | 80.69 | 275.80 | 38.97
gain error | 48.87 | 40.07 | 36.52 | 34.24 | 34.25 | 152.00 | 48.70 | 48.87 | 80.69 | 273.01 | 48.97

Table 3: The standard deviation of the difference between a map made with a
correlated and uncorrelated noise model, in units of $10^{-4}$ of the standard
deviation of the map itself, for each mapmaking method, and for a smooth,
noiseless CMB-like field. The fractional leakage is typically O($4\cdot
10^{-3}$). This is larger than in table 2 because every point is a source of
leakage, not just a single source, but since it is relative to the CMB instead
of a bright point source, the amplitude is negligible in absolute terms. The
white noise mapping method is an outlier, with relative error up to almost 3%
of the signal itself, but since this is a source-targeted method this would
only happen near the targeted sources. Unlike table 2 the same mapmaking
method is used for the maps being subtracted, to avoid penalizing methods with
different pixel windows.

Figure 4 shows maximum-likelihood maps, residual RMS and deviation from the
uncorrelated solution for each method for a smooth, CMB-like field. It doesn’t
really make sense to apply the source-targeted methods when there is no source
present, but I still include them with a dummy point source with zero
amplitude in the patch center to see how these methods interact with the
surrounding CMB.

With a smooth, low-dynamic range field like this, the leakage is small enough
that the maps are practically identical. From table 3 and the bottom 3rd of
figure 4 we see that the leakage is typically more than 200 times smaller than
the signal itself, even in the presence of gain and pointing errors.

The exception is white source mapping, which has X-shaped leakage at 3% of the
signal strength. This happens because the maximum likelihood in-painting of
the samples that hit the source is done independently for the horizontal and
vertical data set, which therefore end up disagreeing with each other in this
region. As we have seen before, such disagreement is interpreted as noise, and
spreads out by a noise correlation length. Like the X around strong point
sources, this too is a form of model error bias, but since it’s sourced by the
CMB itself rather than a potentially very bright point source, it is never
more than a small fraction of the CMB. A 3%-level contamination in a small
fraction of the sky is too small to matter in most circumstances, but to be
safe one can always use source sub-sampling instead, which doesn’t have this
issue.

##  6 A note on polarization

The toy example used for illustration in this article does not include
polarization in order to keep the code as simple and easy to understand as
possible. That said, the algorithms described here are general and
straightforward to extend to polarization (amounting mainly to a change in the
implementation of the response matrix $P$), and with the exception of PGLS and
XGLS have been implemented successfully in the Atacama Cosmology Telescope
mapmaking framework enlib.

That said, it can be useful to consider the what the polarized properties of
model error artifacts are. The maximum-likelihood map $\hat{m}=M_{PN}d$ has a
noise covariance matrix $C=(P^{T}N^{-1}P)^{-1}$, and sub-pixel variance in
some pixel $i$ source artifacts $a_{j}\sim C^{\frac{1}{2}}_{ij}$. That is: the
artifacts leaking from a pixel are realizations of the noise covariance
structure around that pixel. We can use this to predict whether artifacts will
introduce significant I-to-P leakage. For a polarized map we can slice the
covariance matrix $C$ into Stokes and pixel dimensions: $C_{\alpha\beta ij}$,
where greek letters index Stokes parameters and latin one index pixels.
Artifacts will have significant I-to-P leakage if $C_{IP}$ is not very small
relative to the pure polarization components $C_{PP}$. Here $P$ is a stand-in
for Q or U.

How large should we expect these terms to be a priori? In telescopes that
measure polarization via simultanous differencing, where pairs (or more) of
detectors with different polarization orientation measure the same spot on the
sky simultanously, atmospheric noise can be very effectively suppressed in
polarization, leading to low $C_{IP}$. The same is the case for telescopes
with (rapidly) rotating half-wave plates or other polarization modulators.
This group of telescopes would have negligible I-to-P leakage in artifacts.
P-to-P artifacts - that is artifacts in polarization maps thate are sourced by
strong, local sources in polarization, like polarized point sources - would
also be greatly reduced because less correlated noise in polarization results
in a short correlation length and hence much more compact artifacts.

On the other hand, telescopes that rely on combining non-simultanous
measurements to separate the Stokes components would not suppress atmospheric
noise in polarization. An example of this would be a telescope with non-
polarization-sensitive detectors that uses a stepping half-wave plate for
polarization separation. Here $C_{IP}$ would be large, and artifacts sourced
by unpolarized sources could end up being almost 100% polarized.

##  7 Summary

The model error mitigation strategies I have looked at can be grouped into
four main classes:

  1. 1.

Higher-order interpolation map-making, which attempts to reduce sub-pixel
errors by making the model interpolate smoothly between pixel centers instead
of using a step function. These are blind, in that no prior knowledge of the
sky is needed to apply them. They deal moderately well with sub-pixel errors,
typically reducing them by a factor of 3-20 in these tests (except very close
to the source). They do not help at all with variability errors, however. As
far as I am aware, these are new to the literature.

  2. 2.

Source targeted methods, which target a pre-defined set of regions of the sky,
typically based on a point source database. With the exception of source
subtraction, these can completely eliminate leakage from the targeted objects
both from sub-pixel errors and from variability errors. Of these, source white
mapping and source sub-sampling are new.

  3. 3.

Artifact subtraction, of which PGLS is the only example. It has the advantage
of being a fast, blind method that handles both sub-pixel and variability
errors moderately well, at the cost of some bias and increased correlated
noise in the maps.

  4. 4.

Signal-aware noise models, of which XGLS is the only example. It is by far the
highest quality blind method, at the cost of high complexity and extremely
high computational cost.

Given the high implementational complexity and limited improvement from
higher-order map-making, I do not think these methods are worth it, especially
given the low impact of model errors outside the neighborhood of strong point
sources or similar high-contrast objects. XGLS performs much better, tieing
with the best source-targeted methods without needing any targeting itself.
But its cost is probably prohibitively high for realistic CMB datasets.
Instead, I recommend source targeted methods, and in particular the source
sub-sampling method. Of the four source-targeted methods I considered, it has
the best performance, eliminating leakage even for variability errors without
leaving a hole in the map (as source cutting does) or introducing any
secondary leakage (as white source mapping does). It is also quite
straightforward to apply, sharing most of its implementation with an existing
technique for handling per-sample data cuts. PGLS is also worth considering if
a blind method is needed, due to its simplicity and speed.

## Acknowledgments

The Atacama Cosmology Telescope data set provided the inspiration for this
paper. I would like to thank Jon Sievers and Reijo Keskitalo for useful
discussion about prior art. Flatiron Institute is supported by the Simons
Foundation.

## References

  * ACT Collaboration (2017) ACT Collaboration. 2017, arXiv:1610.02360, J. Cosmology Astropart. Phys, 6, 031, The Atacama Cosmology Telescope: two-season ACTPol spectra and parameters
  * BICEP2 Collaboration (2014) BICEP2 Collaboration. 2014, arXiv:1403.3985, BICEP2 I: Detection Of B-mode Polarization at Degree Angular Scales
  * Dünner et al. (2013) Dünner, R., et al. 2013, arXiv:1208.0050, ApJ, 762, 10, The Atacama Cosmology Telescope: Data Characterization and Mapmaking
  * Keihänen et al. (2010) Keihänen, E., Keskitalo, R., Kurki-Suonio, H., Poutanen, T., & Sirviö, A. S. 2010, arXiv:0907.0367, A&A, 510, A57, Making cosmic microwave background temperature and polarization maps with MADAM
  * Piazzo (2017) Piazzo, L. 2017, IEEE Transactions on Image Processing, 26, 5232, Least Squares Image Estimation for Large Data in the Presence of Noise and Irregular Sampling
  * Piazzo et al. (2012) Piazzo, L., Ikhenaode, D., Natoli, P., Pestalozzi, M., Piacentini, F., & Traficante, A. 2012, IEEE Transactions on Image Processing, 21, 3687, Artifact Removal for GLS Map Makers by Means of Post-Processing
  * Planck Collaboration (2018) Planck Collaboration. 2018, arXiv:1807.06207, arXiv e-prints, arXiv:1807.06207, Planck 2018 results. III. High Frequency Instrument data processing and frequency maps
  * Ruud et al. (2015) Ruud, T. M., et al. 2015, arXiv:1508.02778, ApJ, 811, 89, The Q/U Imaging Experiment: Polarization Measurements of the Galactic Plane at 43 and 95 GHz
  * Schaffer et al. (2011) Schaffer, K. K., et al. 2011, arXiv:1111.7245, ApJ, 743, 90, The First Public Release of South Pole Telescope Data: Maps of a 95 deg2 Field from 2008 Observations
  * Tegmark (1997) Tegmark, M. 1997, astro-ph/9611130, ApJ, 480, L87, How to Make Maps from Cosmic Microwave Background Data without Losing Information

##  Appendix A Bicubic spline interpolation

The forward operation $d=Pm$ can be done in two parts. The first is a
convolution of the whole map by a spline pre-filter, followed by the actual
per-sample interpolation. This is implemented in
scipy.ndimage.map_coordinates, and when specializing to 2 dimensions and
ignoring boundary conditions and optimization, it can be written in pseudo-
Python like this:

⬇

def interpol(m, y, x):

m = pre_filter(m)

d = eval(m, y, x)

return d

def pre_filter(m):

m = m.copy()

ny, nx = m.shape

for py in range(ny): filter_1d(m[py,:])

for px in range(nx): filter_1d(m[:,px])

return m

def filter_1d(a):

n = len(a)

p = sqrt(3)-2

for i in range(1,n):

a[i] += p*a[i-1]

a[n-1] *= -p/(1-p**2)

for i in range(n-2,-1,-1):

a[i] = p*(a[i+1]-a[i])

def eval(m, y, x):

nsamp = len(y)

for i in range(nsamp):

py, px = floor(y[i]), floor(x[i])

# offset of point from floored position

ry, rx = y[i]-py, x[i]-px

wy, wx = get_weights(ry), get_weights(rx)

for jy in range(4):

for jx in range(4):

d[i] += wy[jy]*wx[jx]*m[py+jy-1,px+jx-1]

return d

def get_weights(r):

for j in range(4):

# distance of each reference pixel from point

d = abs(r+j-1)

w[j] = (d*d*(d-2)*3+4)/6 if d < 1 else (2-d)**3/6

return w

The transpose operation $m=P^{T}d$ is obtained by reversing the data flow:

⬇

def interpol_trans(d, y, x):

m = eval_trans(d, y, x)

m = pre_filter_trans(m)

return m

def pre_filter_trans(m):

m = m.copy()

ny, nx = m.shape

for px in range(nx): filter_1d_trans(m[:,px])

for py in range(ny): filter_1d_trans(m[py,:])

return m

def filter_1d_trans(a):

n = len(a)

p = sqrt(3)-2

a[0] *= -p

for i in range(1, n-1):

a[i] = p*(a[i-1]-a[i])

a[n-1] -= a[n-2]

a[n-1] *= -p/(1-p**2)

for i in range(n-2,0,-1):

a[i] += p*a[i+1]

def eval_trans(d, y, x):

nsamp = len(y)

m = zeros([ny,nx])

for i in range(nsamp):

py, px = floor(y[i]), floor(x[i])

# offset of point from floored position

ry, rx = y[i]-py, x[i]-px

wy, wx = get_weights(ry), get_weights(rx)

for jy in range(4):

for jx in range(4):

m[py+jy-1,px+jx-1] += wy[jy]*wx[jx]*d[i]

return m

This transposed bicubic spline operation is missing from sicpy, but is
available in pixell.interpol.map_coordinates. See its source code for the full
details (https://github.com/simonsobs/pixell).


1.25in1.25in0.8in1.5in

# Bounds via spectral radius-preserving row sum expansions

Joseph P. Stover111Department of Mathematics, Gonzaga University, 502 E Boone
Ave, Spokane, WA 99258. E-mail<EMAIL_ADDRESS>
###### Abstract

We show a simple method for constructing larger matrices but preserving the
spectral radius. This yields a sufficient criteria for two square matrices of
arbitrary dimension have the same spectral radius, a way to compare spectral
radii of two matrices, and a way to derive new upper and lower bounds on
spectral radius which give the standard row sum bounds as a special case.

††footnotetext: Keywords: block matrix; spectral radius; nonnegative
matrix††footnotetext: AMS 2020 Subject Classification: 15A18, 15A42

##  1 Introduction

For non-negative square matrices it is well known that the spectral radius is
between the minimum and maximum row sums. If the row rums are constant, then
that sum is the spectral radius. It is also know that equality between the
spectral radius and either the minimal or maximal row sum for an irreducible
matrix implies its row sums are constant. See Theorem 1.1 in [2]. Here we will
take a square matrix and increase its dimension by expanding one of its
diagonal elements (or a square diagonal block) into a larger square block
while still preserving the spectral radius of the original matrix. Then we use
this to derive new bounds on the spectral radius and show the standard row sum
bounds are a special case.

Let $M\geq 0$ be a nonnegative square matrix with spectral radius $\rho(M)$.
Consider diagonal element $M_{ii}=d$. Via a permutation we can simply assume
$d$ is the lower right most diagonal component

$M=\left[\begin{matrix}A&\mathbf{b}\\\ \mathbf{c}^{\top}&d\end{matrix}\right]$

where $A$ is an $(n-1)\times(n-1)$ block and $\mathbf{b}$ and $\mathbf{c}$ are
$(n-1)\times 1$ vectors.

We now create a new matrix

$\tilde{M}=\left[\begin{matrix}A&B\\\ C&D\end{matrix}\right].$

where $D$ is a square matrix of arbitrary dimension $s\times s$ with constant
row sums equal to $d$, $B$ is a matrix whose $i^{th}$ row sum equals the
$i^{th}$ element of $\mathbf{b}$ (but arbitrary entries otherwise), and $C$ is
a matrix with each row identical to $\mathbf{c}^{\top}$. We refer to
$\tilde{M}$ as a row sum expansion of $M$ on $d$. The overal dimensions of
$\tilde{M}$ are now $(n-1+s)\times(n-1+s)$.

It is in fact the case that $\rho(M)=\rho(\tilde{M})$ (which we prove in
Theorem 2.1 below). Similarly, one can perform column sum expansions instead
(with the rules on $B$ and $C$ symmetrically adjusted). One can also expand
multiple diagonal elements simultaneously, some using row sums and some using
column sums, as long as the adjacent blocks are created appropriately. In
fact, one can perform any sequence of permutations, transposes, and various
mixes of column and row sum expansions (always transforming diagonal elements
into square blocks) while still preserving the spectral radius. The resulting
matrix (with larger dimension) always has the same spectral radius.

We first give some examples of row and column sum expansions to familiarize
the reader with the process. Then in Section 2 we prove that this process
preserves the spectral radius. Finally in Section 3 we give bounds on the
spectral radius of a matrix based on this idea and give a way to compare
spectral radii of two matrices with different dimensions.

### Examples of matrix expansions

Let $A$ and $B$ be nonnegative square matrices, we write $A\sim B$ if and only
if $\rho(A)=\rho(B)$.

When we perform an expansion of a diagonal element into a block, we sometimes
add additional parentheses to illustrate the new block when was previously a
single component. Here is an introductory example:

$[5]\sim\left[\begin{matrix}2&3\\\
4&1\end{matrix}\right]\sim\left[\begin{matrix}2&4\\\
3&1\end{matrix}\right]\sim\left[\\!\\!\begin{array}[]{ccc}\left(\begin{array}[]{cc}1&1\\\
2&0\end{array}\right)&\left(\begin{array}[]{c}4\\\ 4\end{array}\right)\\\
\left(\begin{array}[]{cc}1&2\end{array}\right)&1\end{array}\\!\\!\right].$

where the first step we simply created a $2\times 2$ matrix with constant row
sums of 5 and then transposed the result in the next step. The final step, we
row-sum expanded the first diagonal element. This requires us to repeat the
horizontal off-diagonal component 4 (for a row sum expansion we always
replicate exactly any rows to the left or right) and we can expand the
vertical off-diagonal element 3 into an arbitrary $1\times 2$ block as long as
its row sum is 3 (for a row sum expansion we are allowed to replicate any
above or below matrix elements into arbitrary rows that each sum to the
original component they are expanded from).

Here is an example with a simultaneous row sum expansion on both diagonal
elements into different sized diagonal blocks. In this example, we simply
require the row sums in each block to be constant and equal to the original
matrix component being expanded.

$\left[\begin{matrix}5&7\\\
2&4\end{matrix}\right]\sim\left[\\!\\!\begin{array}[]{ccc}\left(\begin{array}[]{cc}3&2\\\
0&5\end{array}\right)&\left(\begin{array}[]{ccc}1&4&2\\\
3&1.5&2.5\end{array}\right)\\\ \left(\begin{array}[]{cc}1&1\\\ 0&2\\\
1.5&0.5\end{array}\right)&\left(\begin{array}[]{ccc}3&0.5&0.5\\\ 1&1&2\\\
0&4&0\end{array}\right)\end{array}\\!\\!\right].$

Here is an example where we perform a column expansion on the first diagonal
element and a row expansion on the second. In this case, the lower left block
must be the lower left matrix component repeated in all positions, but the
upper right block just needs to sum to the original upper right component:

$\left[\begin{matrix}5&8\\\
7&3\end{matrix}\right]\sim\left[\\!\\!\begin{array}[]{ccc}\left(\begin{array}[]{cc}4&2\\\
1&3\end{array}\right)&\\!\\!\\!\left(\begin{array}[]{ccc}0&1\\\
2&5\end{array}\right)\\\ \left(\begin{array}[]{cc}7&7\\\
7&7\end{array}\right)&\\!\\!\\!\left(\begin{array}[]{ccc}3&0\\\
2&1\end{array}\right)\end{array}\\!\\!\right].$

Again, we could go on concocting complicated mixes of transposes,
permutations, and row and column sum expansions. This allows the overall mix
of row and column sums to become almost any collection of nonnegative numbers,
thus the standard bounds on spectral radius do not make it obvious that the
spectral radius is preserved.

##  2 Row sum expansions preserve spectral radius

We now prove that row sum expansions preserve the spectral radius. Here it is
important that we are expanding diagonal elements into square blocks. In
general, a row expansion that does not expand diagonal elements into square
blocks does not preserve the spectral radius.

######  Theorem 2.1.

Let $M$ be a non-negative matrix of dimension $n\times n$, and let $\tilde{M}$
be a row sum expansion of $M$ on any one of its diagonal elements. Then
$\rho(M)=\rho(\tilde{M})$.

###### Proof.

Without loss of generality, we prove the theorem for the last diagonal
component of $M$ by permutational similarity. Let

$M=\left[\begin{matrix}A&\mathbf{b}\\\ \mathbf{c}^{\top}&d\end{matrix}\right]$

where $A$ is a $(n-1)\times(n-1)$ block, $\mathbf{b}$ and $\mathbf{c}$ is a
$(n-1)\times 1$ vectors, and $d$ is a single diagonal component.

Since $M\geq 0$, we know that its spectral radius $\rho\geq 0$ is an
eigenvalue with nonnegative, nonzero left eigenvector $(\mathbf{x}^{\top},y)$
where $\mathbf{x}$ has $n-1$ components and $y$ is a single component (see
Theorem 8.3.1 in [1]). Assume the eigenvector sums to one
$y+\sum_{i=1}^{n-1}x_{i}=1$. Note that we do not require any uniqueness, so
there may be several such eigenvectors.

We have that

$\displaystyle(\mathbf{x}^{\top},y)\left[\begin{matrix}A&\mathbf{b}\\\
\mathbf{c}^{\top}&d\end{matrix}\right]$
$\displaystyle=(\mathbf{x}^{\top}A+y\mathbf{c}^{\top},\mathbf{x}^{\top}\mathbf{b}+yd)$
$\displaystyle=(\rho\mathbf{x}^{\top},\rho y)$

thus giving

$\mathbf{x}^{\top}\mathbf{b}=y(\rho-d)$ (2.1) $\mathbf{x}^{\top}(\rho
I_{n-1}-A)=y\mathbf{c}^{\top}$ (2.2)

It is already known that $\rho\geq d$. If $\rho=d$, then $\mathbf{x}$ is
orthogonal to $\mathbf{b}$ and thus has zeros wherever $\mathbf{b}$ has
positive components (recall that everything is nonnegative). If $\rho$ is also
an eigenvalue of principle submatrix $A$, then we can find infinitely many
$\mathbf{x}$, but, again, we do not require uniqueness here. We just require
some nonzero, nonnegative eigenvector to exist.

We now create a new $(n-1+s)\times(n-1+s)$ matrix $\tilde{M}\geq 0$ by
expanding $d$ into a $s\times s$ matrix $D$ with constant row sums
$\sum_{j=1}^{s}D_{ij}=d$ for all $i=1,2,\ldots,s$.

$\tilde{M}=\left[\begin{matrix}A&B\\\ C&D\end{matrix}\right]$

where $B$ is an $(n-1)\times s$ block whose row sums equal the elements of
vector $\mathbf{b}$: $\sum_{j=1}^{s}B_{ij}=b_{i}$ for $i=1,2,\ldots,n-1$, and
$C$ is an $s\times(n-1)$ block where each row is identical
$C_{i-}=\mathbf{c}^{\top}$ for each $i=1,2,\ldots,s$.

Similarly, this expanded matrix has a nonnegative spectral radius
$\tilde{\rho}\geq 0$ which is an eigenvalue with nonnegative, nonzero
eigenvector $(\tilde{\mathbf{x}}^{\top},\mathbf{y}^{\top})$. Again we assume
that $\sum_{i=1}^{s}y_{i}+\sum_{i=1}^{n-1}x_{i}=1$. We will show that
$\tilde{\mathbf{x}}=\mathbf{x}$ and $\sum y_{i}=y$ and $\tilde{\rho}=\rho$.

First, we show that $\rho,({\mathbf{x}}^{\top},\mathbf{y}^{\top})$ is an
eigenpair for $\tilde{M}$ with $\mathbf{y}$ being a vector that solves

${\mathbf{x}}^{\top}B=\mathbf{y}^{\top}(\rho I_{s}-D).$ (2.3)

There may be several such vectors $\mathbf{y}$ (if $\rho=d$ the constant row
sum for matrix $D$). We choose a $\mathbf{y}$ so that this eigenvector sums to
one requiring $y=\sum_{i=1}^{s}y_{i}=1-\sum_{i=1}^{n-1}x_{i}$. This is
possible since the row sums of $\mathbf{x}^{\top}B$ are all equal to
$\mathbf{x}^{\top}\mathbf{b}$ and the row sums of $\mathbf{y}^{\top}(\rho
I_{s}-D)$ are all equal to
$\left(\sum_{i=1}^{s}y_{i}\right)(\rho-d)=y(\rho-d)$ which are equivalent by
(2.1).

Since block $C$ simply has vector $\mathbf{c}^{\top}$ repeated $s$ times for
its rows, we have

$\mathbf{y}^{\top}C=\left(\sum_{i=1}^{s}y_{i}\right)\mathbf{c}^{\top}=y\mathbf{c}^{\top}.$
(2.4)

Using (2.2) and (2.4) we see that

${\mathbf{x}}^{\top}A+\mathbf{y}^{\top}C={\mathbf{x}}^{\top}A+y\mathbf{c}^{\top}=\rho{\mathbf{x}}^{\top}.$
(2.5)

Using (2.3), we see that

${\mathbf{x}}^{\top}B+\mathbf{y}^{\top}D=\rho\mathbf{y}^{\top}.$ (2.6)

Now we have by using (2.5) and (2.6) that

$\displaystyle({\mathbf{x}}^{\top},\mathbf{y}^{\top})\left[\begin{matrix}A&B\\\
C&D\end{matrix}\right]$
$\displaystyle=({\mathbf{x}}^{\top}A+\mathbf{y}^{\top}C,{\mathbf{x}}^{\top}B+\mathbf{y}^{\top}D).$
$\displaystyle=(\rho{\mathbf{x}}^{\top},{\rho}\mathbf{y}^{\top})$

This shows that $\rho,({\mathbf{x}}^{\top},\mathbf{y}^{\top})$ is an eigenpair
for $\tilde{M}$. This establishes that $\tilde{\rho}\geq\rho$.

Now we show that $\tilde{\rho}=\rho$ by contradiction. Assume that
$\tilde{\rho}>\rho$ and that we have eigenpair
$\tilde{\rho},(\tilde{\mathbf{x}}^{\top},\tilde{\mathbf{y}}^{\top})$ for
$\tilde{M}$. Then we have

$\tilde{\mathbf{x}}^{\top}B=\tilde{\mathbf{y}}^{\top}(\tilde{\rho}I_{s}-D)$

which, after checking the row sums on each side, yields

$\tilde{\mathbf{x}}^{\top}\mathbf{b}=\left(\sum_{i=1}^{s}\tilde{y}_{i}\right)(\tilde{\rho}-d).$
(2.7)

And of course, we also have that

$\tilde{\mathbf{x}}^{\top}(\tilde{\rho}I_{n-1}-A)=\tilde{\mathbf{y}}^{\top}C=\left(\sum_{i=1}^{s}\tilde{y}_{i}\right)\mathbf{c}^{\top}$
(2.8)

Together, (2.7) and (2.8) and using (2.1) and (2.2) with
$\tilde{\rho},\tilde{\mathbf{x}},\tilde{y}=\sum_{i=1}^{s}\tilde{y}_{i}$
instead of $\rho,\mathbf{x},y$ show that
$\tilde{\rho},(\tilde{\mathbf{x}}^{\top},\tilde{y})$ is an eigenpair for $M$
which contradicts the fact that $\rho$ is its spectral radius. Hence
$\tilde{\rho}=\rho$.

∎

This theorem applies to column sum expansions as well (since this is
equivalent to row sum expansions on $A^{\top}$), and hence it applies be any
sequence of row and column sum expansions, permutations, and transposes since
each step preserves the spectral radius.

##  3 Spectral radius bounds via contractions

If we were to perform the reverse procedure, contracting our matrix to lower
its dimension, then unless the row or column sums in all blocks are constant
(maybe this is achievable with some permutation or transpose), we would have
to break the preservation of spectral radius. However, we use this to derive
bounds on the spectral radius.

Again, letting $M$ be an $n\times n$ nonnegative square matrix, we define
upward and downward row and column sum contractions.

######  Definition 3.1.

Let $M$ be a nonnegative square matrix. A _downward row sum contraction_ of
$M$ is denoted $M^{\downarrow}$ and is a matrix created as follows. First we
partition $M$ into a block matrix with $k^{2}$ blocks denoted $M=(A_{ij})$ for
$i,j$ in $\\{1,2,...,k\\}$ requiring that all diagonal blocks $A_{ii}$ are
square matrices. In each block, we take the minimum row sum: if block $A_{ij}$
has size $s\times t$ then $a_{ij}=\min_{u}\\{\sum_{v=1}^{t}(A_{ij})_{uv}\\}$.
Then we have $M^{\downarrow}_{ij}=a_{ij}$.

Similarly we define _upward row sum contraction_ of $M$ denoted by
$M^{\uparrow}$ which uses the maximum row sum in each block:
$M^{\uparrow}=(b_{ij})$ where
$b_{ij}=\max_{u}\\{\sum_{v=1}^{t}(B_{ij})_{uv}\\}$. Naturally, we also define
upward and downward column sum contractions (which are identical to applying
row sum contractions to $M^{\top}$). We could also perform any sequence of
permutations, transposes, and row and column sum contractions to arrive at a
new matrix of lower dimension.

### An example

Consider the matrix $M$ partitioned as follows:

$M=\left[\\!\begin{array}[]{c|ccc|cc}1&2&0&3&7&1\\\ \hline\cr 4&1&0&5&1&0\\\
0&2&3&0&0&5\\\ 5&2&0&1&3&1\\\ \hline\cr 0&0&3&5&0&2\\\
4&1&3&0&1&3\end{array}\\!\right].$

Then taking the minimum and maximum row sum in each block gives

$M^{\downarrow}=\left[\\!\begin{array}[]{ccc}1&5&8\\\ 0&3&1\\\
0&4&2\end{array}\\!\right],\quad
M^{\uparrow}=\left[\\!\begin{array}[]{ccc}1&5&8\\\ 5&6&5\\\
4&8&4\end{array}\\!\right].$

Notice that before reducing the dimensions of the matrix, we can increase or
decrease components of $M$ in order to create an ordering. For this example,
here is one possibility with the changed components in bold:

$\tilde{M}^{\downarrow}=\left[\\!\begin{array}[]{c|ccc|cc}1&2&0&3&7&1\\\
\hline\cr\mathbf{0}&1&0&\mathbf{2}&1&0\\\ 0&2&\mathbf{1}&0&0&\mathbf{1}\\\
\mathbf{0}&2&0&1&\mathbf{0}&1\\\ \hline\cr 0&0&3&\mathbf{1}&0&2\\\
\mathbf{0}&1&3&0&1&\mathbf{1}\end{array}\\!\right]\leq\left[\\!\begin{array}[]{c|ccc|cc}1&2&0&3&7&1\\\
\hline\cr 4&1&0&5&1&0\\\ 0&2&3&0&0&5\\\ 5&2&0&1&3&1\\\ \hline\cr
0&0&3&5&0&2\\\
4&1&3&0&1&3\end{array}\\!\right]\leq\left[\\!\begin{array}[]{c|ccc|cc}1&2&0&3&7&1\\\
\hline\cr\mathbf{5}&1&0&5&\mathbf{5}&0\\\ \mathbf{5}&2&3&\mathbf{1}&0&5\\\
5&2&\mathbf{3}&1&3&\mathbf{2}\\\ \hline\cr\mathbf{4}&0&3&5&0&\mathbf{4}\\\
4&1&3&\mathbf{4}&1&3\end{array}\\!\right]=\tilde{M}^{\uparrow}$

It should be intuitively clear that this will always be possible, to
arbitrarily decrease or increase some components to make all row sums in each
block constant and equal to the minimum or maximum row sum. It is important
that we keep the diagonal blocks as square principle submatrices in order to
use Theorem 2.1. As mentioned above already, a row expansion that does not
expand diagonal elements into square blocks does not generally preserve the
spectral radius. Hence if our diagonal block of constant sum are now square,
then reducing the dimension can change the spectral radius, which we wish to
avoid here.

It is a standard result for ordered square nonnegative matrices that $A\leq B$
(component-wise: $A_{ij}\leq B_{ij}$ for all $i,j$) that the spectral radii
are also ordered $\rho(A)\leq\rho(B)$ (see Corollary 8.1.20 in [1]). Thus we
get

$\rho(\tilde{M}^{\downarrow})\leq\rho(M)\leq\rho(\tilde{M}^{\uparrow})$

for any arbitrary downward and upward adjustment of $M$ to make constant row
sums in each block. We formalize this as a theorem in the general case.

Let $C_{s}^{\downarrow}(M)$ be the set of all $s\times s$ downward
contractions of matrix $M$ and $C_{s}^{\uparrow}(M)$ be the set of all
$s\times s$ upward contractions. Each $M^{\downarrow}\in
C_{s}^{\downarrow}(M)$ is an $s\times s$ matrix whose components are equal to
the minimum row sum of each block for some partition of $M$, and similarly for
any $M^{\uparrow}\in C_{s}^{\uparrow}(M)$. Note that these are distinct from
$\tilde{M}^{\downarrow}$ and $\tilde{M}^{\uparrow}$ which have the same
dimensions as $M$ but with some components adjusted downward or upward,
respectively, but resulting in the ordering: $\tilde{M}^{\downarrow}\leq
M\leq\tilde{M}^{\uparrow}$. By Theorem 2.1 we have that
$\rho(M^{\downarrow})=\rho(\tilde{M}^{\downarrow})$ even though they have
different sizes and also that $\rho(M^{\uparrow})=\rho(\tilde{M}^{\uparrow})$.

Note that trivially, $[\min_{i}\\{\sum_{j=1}^{n}M_{ij}\\}]\in
C_{1}^{\downarrow}(M)$ and $[\max_{i}\\{\sum_{j=1}^{n}M_{ij}\\}]\in
C_{1}^{\uparrow}(M)$ for the minimal and maximal row sums of $M$ viewed as
$1\times 1$ matrices and that $M\in C_{n}^{\downarrow}(M)\cap
C_{n}^{\uparrow}(M)$.

Let $C^{\downarrow}(M)=\cup_{j=1}^{n-1}C_{j}^{\downarrow}(M)$ and
$C^{\uparrow}(M)=\cup_{j=1}^{n-1}C_{j}^{\uparrow}(M)$ be the sets of all
possible downward and upward contractions of $M$ (with dimensions less than
$M$). It should be clear that each of these sets is finite. The process of
creating a downward or upward contraction involves some finite sequence of
permutations, transposes, and dimension decreases. Eventually, including more
permutations or transposes will create no new matrices, and obviously the
dimension can only be decreased to $1$ at the lowest. Transpose and
permutations also preserve the spectral radius, so it is only the dimension
reductions that can make the spectral radius smaller or larger. Note that a
matrix in $C^{\downarrow}(M)$ may include several dimension decrease steps in
its creation, which could mean the spectral radius undergoes several decreases
relative to $M$, and similar for $C^{\uparrow}(M)$.

######  Theorem 3.2.

Let $C^{\downarrow}(M)$ and $C^{\uparrow}(M)$ be the set of all downward and
upward contractions of nonnegative square matrix $M$. Then

$\max_{M^{\downarrow}\in
C^{\downarrow}(M)}\\{\rho(M^{\downarrow})\\}\leq\rho(M)\leq\min_{M^{\uparrow}\in
C^{\uparrow}(M)}\\{\rho(M^{\uparrow})\\}.$

###### Proof.

For any $M^{\downarrow}\in C^{\downarrow}(M)$, we can create a matrix
$\tilde{M}^{\downarrow}$ by adjusting individual components of $M$ downward
and have $\rho(M^{\downarrow})=\rho(M^{\downarrow})$ by Theorem 2.1. We also
have that $\tilde{M}^{\downarrow}\leq M$, and similarly for any
$M^{\uparrow}\in C^{\uparrow}(M)$. So that we always have
$\rho(M^{\downarrow})=\rho(\tilde{M}^{\downarrow})\leq\rho(M)\leq\rho(\tilde{M}^{\uparrow})=\rho(M^{\uparrow})$.
The result follows by taking the maximum on the left and minimum on the right.
∎

The standard bounds on $\rho(M)$ using its row sums is a special case of
Theorem 2.1.

######  Corollary 3.3.

We have that $[\min_{i}\\{\sum_{j=1}^{n}M_{ij}\\}]\in C_{1}^{\downarrow}(M)$
and $[\max_{i}\\{\sum_{j=1}^{n}M_{ij}\\}]\in C_{1}^{\uparrow}(M)$ which gives
$\min_{i}\\{\sum_{j=1}^{n}M_{ij}\\}\leq\rho(M)\leq\max_{i}\\{\sum_{j=1}^{n}M_{ij}\\}$.

There is no guarantee that Theorem 3.2 gives better bounds on $\rho(M)$ than
the standard row sum bounds, but it does improve the bounds in at least some
cases (some examples are given below).

Here is another direct corollary which just considers particular $2\times 2$
row sum contractions, for which calculating the spectral radius is
straightforward.

######  Corollary 3.4.

Let $A=(a_{ij})$ be a $2\times 2$ downward row sum contraction of $M$ and
$B=(b_{ij})$ be a $2\times 2$ upward row sum contraction of $M$. Then

$\frac{\operatorname{tr}(A)}{2}+\sqrt{\left(\frac{\operatorname{tr}(A)}{2}\right)^{2}-\det(A)}\leq\rho(M)\leq\frac{\operatorname{tr}(B)}{2}+\sqrt{\left(\frac{\operatorname{tr}(B)}{2}\right)^{2}-\det(B)}$

Again we call attention to the fact that we have required the diagonal blocks
to be square submatrices as it is not true in general when the diagonal blocks
are allowed to have arbitrary dimensions. In the gneral case, taking the
minimum row sums of each block will not necessarily give a lower spectral
radius. This is because we are relying on Theorem 2.1 to preserve the spectral
radius when contracting the (downward and upward component-adjusted) matrices
to a smaller size.

Note that rather than performing contractions, one could perform expansions
and then adjust elements up or down appropriately to get different bounds on
the spectral radius, but generally calculating spectral radius for a larger
matrix is more difficult or computationally intensive, so we don’t state that
result separately.

Here is an example application of Corollary 3.4. We give matrix $M$ and the
bounds on its spectral radius from row sums alone.

$M=\left[\\!\begin{array}[]{ccc}1&3&2\\\ 5&1&1\\\
2&4&3\end{array}\\!\right]\quad 6\leq\rho(M)\leq 9$

Now we row-sum contract the upper left $2\times 2$ block of $M$ to produce
downward and upward contractions with spectral radii given below.

$\displaystyle M^{\downarrow}=\left[\\!\begin{array}[]{ccc}4&1\\\
6&3\end{array}\\!\right],M^{\uparrow}=\left[\\!\begin{array}[]{ccc}6&2\\\
6&3\end{array}\\!\right]\Rightarrow\rho(M^{\downarrow})=6\leq\rho(M)\leq
8.3\approx\rho(M^{\uparrow})$

Thus our method gives a slight refinement over considering the row sums of $M$
alone. If we considered all possible $2\times 2$ contractions of this matrix,
the best bounds are approximately $6.3\approx
2+\sqrt{19}\leq\rho(M)\leq(9+\sqrt{57})/2\approx 8.3$. By considering column
sums of $M$ we get $6\leq\rho(M)\leq 8$ and by considering all $2\times 2$
column sum contractions we get approximately
$6.5\approx(5+\sqrt{65})/2\leq\rho(M)\leq 8$, still a slight refinement over
the basic method.

Here is an example where there is a slight refinement of the lower bound but a
more significant refinement on the upper bound (relative to row or column sums
alone).

$M=\left[\\!\begin{array}[]{ccccc}\framebox{1}&3&\framebox{2}&1&\framebox{2}\\\
7&\underline{1}&1&\underline{3}&3\\\
\framebox{2}&4&\framebox{3}&1&\framebox{0}\\\
1&\underline{1}&5&\underline{2}&2\\\
\framebox{4}&3&\framebox{0}&2&\framebox{1}\end{array}\\!\right]\quad
8\leq\rho(M)\leq 15$

The best bounds are found with contracting the principle submatrices on
indices $\\{1,3,5\\}$ (boxed above) and on $\\{2,4\\}$ (underlined above) each
into a single component (i.e. using permutation $(1,2,3,4,5)\to(1,3,5,2,4)$
and partitioning into a $3\times 3$ upper left diagonal block and $2\times 2$
lower right diagonal block) which results in the downward and upward row sum
contraction matrices

$M^{\downarrow}=\left[\\!\begin{array}[]{ccc}5&4\\\
8&3\end{array}\\!\right],M^{\uparrow}=\left[\\!\begin{array}[]{ccc}5&5\\\
11&4\end{array}\\!\right]$

giving the bounds

$9.74\approx
4+\sqrt{33}=\rho(M^{\downarrow})\leq\rho(M)\leq\rho(M^{\uparrow})=\frac{9+\sqrt{153}}{2}\approx
11.93$

with the actual value being approximately $\rho(M)\approx 10.995$.

Here is another straightforward corollary which discusses comparing the
spectral radii of two matrices of arbitrarily dimensions.

######  Corollary 3.5.

Let $A,B$ be two square nonnegative matrices or arbitrary and possibly
distinct dimension. If we can apply some sequence of permutations, transposes,
and upward row or column sum contractions or expansions to $A$ and also some
sequence of these operations to $B$ but with downward row or column sum
contractions or expansions to create matrices satisfying $A^{\uparrow}\leq
B^{\downarrow}$ (with both being of the same dimension) or more generally
$\rho(A^{\uparrow})\leq\rho(B^{\downarrow})$ (with neither
$A^{\uparrow}\not\leq B^{\downarrow}$ nor their being of the same dimension
required), then $\rho(A)\leq\rho(B)$.

For example, consider square matrices of possibly distinct arbitrary
dimensions written in block form with blocks of arbitrary dimensions (as long
as the diagonal blocks are square submatrices)
$A=\left[\begin{matrix}A_{11}&A_{12}\\\ A_{21}&A_{22}\end{matrix}\right]$ and
$B=\left[\begin{matrix}B_{11}&B_{12}\\\ B_{21}&B_{22}\end{matrix}\right]$. Let
$a_{ij}$ represent the maximum row sum of block $A_{ij}$,
$A^{\uparrow}=(a_{ij})$, and $b_{ij}$ the minimum row sum of block $B_{ij}$,
$B^{\downarrow}=(b_{ij})$. If we have that

$\frac{\operatorname{tr}(A^{\uparrow})}{2}+\sqrt{\left(\frac{\operatorname{tr}(A^{\uparrow})}{2}\right)^{2}-\det(A^{\uparrow})}\leq\frac{\operatorname{tr}(B^{\downarrow})}{2}+\sqrt{\left(\frac{\operatorname{tr}(B^{\downarrow})}{2}\right)^{2}-\det(B^{\downarrow})}$

then $\rho(A)\leq\rho(B)$.

Here is an example of comparing two matrices.

$A=\left[\\!\begin{array}[]{c|cc|c}\framebox{2}&1&1&\framebox{2}\\\ \hline\cr
1&1&3&0\\\ 0&0&2&1\\\ \hline\cr\framebox{1}&2&0&\framebox{4}\\\
\end{array}\\!\right]\to\tilde{A}=\left[\\!\begin{array}[]{cc|cc}1&3&1&0\\\
0&2&0&1\\\ \hline\cr 1&1&\framebox{2}&\framebox{2}\\\
2&0&\framebox{1}&\framebox{4}\\\ \end{array}\\!\right]\to
A^{\uparrow}=\left[\\!\begin{array}[]{cc}4&1\\\ 2&5\\\ \end{array}\\!\right]$

Note that we have applied a permutation to $A$ indicated by the boxed entries:
$(1,2,3,4)\to(3,1,2,4)$.

$B=\left[\\!\begin{array}[]{cc|c}1&2&2\\\ 3&1&3\\\ \hline\cr 1&1&5\\\
\end{array}\\!\right]\to B^{\downarrow}=\left[\\!\begin{array}[]{cc}3&2\\\
2&5\\\ \end{array}\\!\right]$

In this case it is clear that $\rho(A^{\uparrow})=6$ due to constant column
sums and that $\rho(B^{\downarrow})=4+\sqrt{16-11}>6$ implying that
$\rho(A)<\rho(B)$ even though there is no obvious relationship between their
spectral radii by simply looking as row or column sums or comparing matrix
components. Note that simply partitioning $A$ without any permutation does not
give the same spectral radius comparison and thus isn’t useful.

## References

  * [1] R. A. Horn and C. R. Johnson.  Matrix Analysis.  Cambridge, 2 edition, 2012.
  * [2] H. Minc.  Nonnegative Matrices.  Wiley, 1 edition, 1988.