# An Improved PTAS for Covering Targets with Mobile Sensors

Nonthaphat Wongwattanakij E-mail<EMAIL_ADDRESS>Department of Computer
Engineering, Kasetsart University, Bangkok, Thailand. Supported by Graduate
Student Research Grant, Faculty of Engineering, Kasetsart University. Nattawut
Phetmak E-mail<EMAIL_ADDRESS>Department of Computer Engineering,
Kasetsart University, Bangkok, Thailand. Chaiporn Jaikaeo E-mail:
<EMAIL_ADDRESS>Department of Computer Engineering, Kasetsart University,
Bangkok, Thailand. Jittat Fakcharoenphol E-mail<EMAIL_ADDRESS>Department
of Computer Engineering, Kasetsart University, Bangkok, Thailand. Supported by
Thailand Reseach Fund, Grant RSA-6180074

###### Abstract

This paper considers a movement minimization problem for mobile sensors. Given
a set of $n$ point targets, the $k$-Sink Minimum Movement Target Coverage
Problem is to schedule mobile sensors, initially located at $k$ base stations,
to cover all targets minimizing the total moving distance of the sensors. We
present a polynomial-time approximation scheme for finding a $(1+\epsilon)$
approximate solution running in time $n^{O(1/\epsilon)}$ for this problem when
$k$, the number of base stations, is constant. Our algorithm improves the
running time exponentially from the previous work that runs in time
$n^{O(1/\epsilon^{2})}$, without any target distribution assumption. To devise
a faster algorithm, we prove a stronger bound on the number of sensors in any
unit area in the optimal solution and employ a more refined dynamic
programming algorithm whose complexity depends only on the width of the
problem.

##  1 Introduction

Wireless sensor networks (WSNs) are used in various types of applications such
as military, industry, and disaster responses [ASSC02, YMG08]. To effectively
monitor a region of interest, sensor nodes are placed at specific locations so
that they can collectively provide overall sensing coverage for the entire
area [Wan11, WZLD20]. WSNs may consist of static sensor nodes, mobile sensor
nodes, or both. With presence of mobile sensor nodes, monitor areas may be
dynamically adjusted by having these mobile nodes navigate to appropriate
positions. This capability helps reducing the overall number of sensor nodes
and/or provides better sensing coverage over regions of interest.

In this paper we study the $k$-Sink Minimum Movement Target Coverage Problem,
where we would like to schedule sensors, with sensing range $r$, from $k$
stations to cover all targets, while minimizing the total distance for every
sensor. Gao, Chen, Wu, and Chen [GCWC17] proposed a polynomial-time
approximation scheme for this problem that runs in time
$n^{O(1/\epsilon^{2})}$ and returns a solution whose cost is at most
$(1+\epsilon)$ of the optimal cost when the targets are distributed uniformly
and randomly in the surveillance region. This result first appeared as
[CGWC16] in INFOCOM’16.

Our contribution is as follows. We observe that, as in many PTAS’s for other
geometric problems, their general framework, based on exhaustive search and
shifting, does not depend on the target distribution; thus it works without
the distribution assumption. More over, we notice that the running time bound
of their algorithm depends on an erroneous claim of the bound on the number of
sensors needed in the optimal solution (in Theorem 3 in [CGWC16] and in
Theorem 5 in [GCWC17]) for a particular region (see Appendix A) and provide a
correction. Our proof actually provides a stronger guarantee, allowing us to
devise a more refined dynamic programming subroutine, the main contribution of
this work, used in the algorithm that runs in time $n^{O(1/\epsilon)}$, which
is an exponential improvement on the running time. While we consider our
result a fairly theoretical one, we hope that the dynamic programming approach
presented here might find more applications in other networking problems.

We describe the network model, problem statements and a review of Gao et al.
[GCWC17, CGWC16] algorithm, called Energy Effective Movement Algorithm (EEMA),
in Section 3. Section 4 describes the improved algorithm and proves the
running time and its performance guarantee. Our running time analysis depends
on the bound on the sensors needed in the optimal solution, proved in Section
5.

##  2 Related works

Coverage problems has been widely studied in computational geometry [O’R87].
These problems are important to wireless sensor networks (see recent survey in
[Wan11] and the book by Wu et al. [WZLD20]). There are numerous heuristic
algorithms for the placement and deployment problem, e.g., the algorithms for
full coverage based on packing has been proposed by Wang and Tseng [WT08] for
areas without obstacles and by Wang, Hu, and Tseng [WHT08] for areas with
obstacles, for coverage enhancement using virtual forces [ZC04], for
strategies that deals with coverage holes using Voronoi polygons [MMA+14], and
for sensor relocation using Delaunay triangulation [KJ17]. This paper as well
as that of Gao et al. [GCWC17] studies movement scheduling problems. For
dynamic area coverage, Liu, Dousse, Nain, and Towsley [LDNT13] considered
coverage using mobile sensor networks. Liu et al [LDWS08] and He et al.
[HCL+12] considered scheduling problem when the goal is to protect an area
using a barrier of sensors. Ammari [Amm12] tried to minimize sensor movement
cost while ensuring area coverage. For strategies for improving coverage,
Wang, Lim, and Ma [WLM09] provided a survey. Guo and Jafarkhani [GJ19]
presented centralized and decentralized heuristics for optimizing sensor
movement under network lifetime constraints. Elhoseny et al. [ETYH18] employed
the genetic algorithm as a heuristic to deal with coverage requirements.

When theoretical guarantees are needed, approximation algorithms for many
versions of coverage problems have been studied. Section 3.2 describes related
work on polynomial-time approximation schemes on geometric graphs closely
related to this work. In this section, we list a few notable works on
approximation algorithms for coverage problems. For target coverage, if the
goal is to schedule sensor activation groups to maximize the life time, Ding
et al. [DWW+12] devised a $(4+\xi)$-approximation algorithm for any $\xi>0$.
For restricted form of $k$-coverage, Xu and Song [XS14] gave a
$3$-approximation algorithm. When a threshold for movement of each sensor is
given and the objective is to maximize target coverage, a recent result by
Liang, Shen, and Chen [LSC21] gave an algorithm with an approximation ratio of
$(1-1/e)$. Another important class of coverage problem is to find dominating
sets in unit disk graphs. Marathe et al. [MBH+95] described a 5-approximation
algorithm for finding dominating sets and a 10-approximation algorithm for
finding connected dominating sets. For weighed dominating sets, Erlebach and
Mihalák [EM10] and Zou et al. [ZWX+11] independently devised a
$(4+\epsilon)$-approximation algorithms. Zou et al. [ZWX+11] also present a
$(5+\epsilon)$-approximation algorithm for the connected version of the
problem. A recent book by Wu et al. [WZLD20] is an excellent source for these
developments.

##  3 Preliminary and Problem Statement

###  3.1 Model and Problem Statement

Following Gao et al. [GCWC17, CGWC16], the homogeneous network considered in
this paper can be modeled as $G=(T,P,r,S)$, where
$T=\\{t_{1},t_{2},\ldots,t_{n}\\}$ is the target set,
$P=\\{p_{1},p_{2},\ldots,p_{k}\\}$ is the station set, $r$ is the sensing
range, and $S=\\{s_{1},s_{2},\ldots,s_{m}\\}$ is the mobile sensor set. The
targets are static points on the plane that we want to cover using sensors
with sensing range $r$. We say that a sensor set $S$ covers a point $t_{i}$ if
there exists a sensor $s_{j}\in S$ such that $d(s_{j},t_{i})\leq r$, where
$d(a,b)$ denotes the Euclidean distance between two points $a$ and $b$ on the
plane. The coverage requirement ensures that every target $t_{i}\in T$ must be
covered. Each sensor $s_{j}\in S$ must be scheduled to move to its location
from some station $p_{\ell}\in P$; the movement distance is
$d(s_{j},p_{\ell})$. Since the station set $P$ is static we can easily find
$p_{\ell}\in P$ that minimizes $d(s_{j},p_{\ell})$. We would later refer to
the distance just as $d(s_{j})$ defined to be $min_{p_{\ell}\in
P}d(s_{j},p_{\ell})$.

Gao et al. [GCWC17, CGWC16] defined the $k$-Sink Minimum Movement Target
Coverage Problem ($k$-MMTC) as follows. Given a set of targets $T$, station
set $P$, and sensing radius $r$, find a set of sensors $S$ covering $T$ such
that the total movement distance is minimized. They proposed an algorithm
called Energy Effective Movement Algorithm (EEMA) that is claimed to be a PTAS
that solves this problem (reviewed in Subsection 3.3).

We would like to note that Gao et al. [GCWC17, CGWC16] also assumed that the
targets are distributed uniformly and randomly in the surveillance region. Our
work does not need this requirement. In fact, we observe that their algorithm
and proofs work even without this distribution assumption as well.

###  3.2 Polynomial-Time Approximation Scheme

A polynomial-time approximation scheme (PTAS) is an algorithm that, given a
constant $\epsilon>0$, returns a solution to the problem of cost within
$(1+\epsilon)$ of the optimal in time polynomial in $n$, the size of the
problem, assuming that $\epsilon$ is a constant. Note that the running time
may depend exponentially on $1/\epsilon$, i.e., the running time might
increase very rapidly as the solution quality increases (as
$1/\epsilon\rightarrow\infty$ when $\epsilon\rightarrow 0$).

A general framework for developing PTAS for geometric problems introduced by
Hunt et al. [HMR+98] follows the approach used in planar graphs by Baker
[Bak94]. First the problem is decomposed into much smaller subproblems which
can be solved separately exhaustively either by dynamic programming or brute-
force search. Then, the “locally-optimal” solutions from these subproblems are
combined to obtain the final solution. Clearly, the combined solution might
not be globally optimal because the dependency between each subproblem is
neglected. However, one can typically ensure that the neglected cost is small,
i.e., at most $1/\epsilon$ fractions of the optimal cost, when the size of
each subproblem is large enough using an averaging argument. For geometric
problems mostly on unit disk graphs, Hunt et al. [HMR+98] consider an area of
size $O(1/\epsilon)\times O(1/\epsilon)$ and use the fact that there is a
feasible solution of size $O(1/\epsilon^{2})$ to devise many PTAS’s that run
in time $n^{O(1/\epsilon^{2})}$. For dominating sets, which is another form of
coverage problems, if the nodes are unweighted, there are several PTAS for the
standard problem [HMR+98, NH05] and for the case when the dominating set is
required to be connected [ZGWD09]. (See [WZLD20] for a comprehensive treatment
of this topics.)

The dynamic programming approach for covering problems is widely used. Erik
Jan van Leeuwen [vL05, vL06, vL09] has employed roughly the same speed-up idea
to the dominating set problem in unit disk graphs.

###  3.3 Review of EEMA [GCWC17, CGWC16]

EEMA proposed by Gao, Chen, Wu, and Chen [GCWC17, CGWC16] works in two phases.
The first phase divides the surveillance region into a set of subareas. This
changes the nature of the problem from continuous optimization problem to a
combinatorial one. The second phase approximates the solution to the
combinatorial problem constructed in the first phase. We shall describe both
phases and our contribution is essentially on the second phase of their
proposed EEMA.

####  3.3.1 First phase of EEMA

For each target $t_{i}$, the detection cycle $D(t_{i})$ is a circle of radius
$r$ centered at $t_{i}$. Note that a sensor in $D(t_{i})$ covers $t_{i}$.
Detection cycles may intersect; thus, a single sensor in the intersection can
detect all targets associated with these cycles. The first phase of EEMA
identifies all minimal subareas defined by intersections of detection cycles
(see Figure 1). Each subarea is called a curved-boundary. Gao et al. proved
the following lemma (implicitly in Theorem 2 in [GCWC17]).

$t_{1}$$t_{2}$$t_{3}$$t_{4}$ Figure 1: Detection cycles for
$t_{1},t_{2},t_{3},$ and $t_{4}$. There are 9 minimal subareas to be
considered.

######  Lemma 3.1 ([GCWC17]).

The number of minimal subareas defined by intersections of all $n$ detection
cycles is $O(n^{2})$

They also gave an $O(n^{2})$-time algorithm (called the Graph Conversion
algorithm) for finding all minimal subareas.

####  3.3.2 Second phase of EEMA

During the second phase, Gao et al. noted that placing a sensor anywhere in a
particular minimal subarea covers the same set of target, defined as covered-
set of that subarea. Therefore, to minimize the total distance, if one has to
schedule a sensor to that subarea, the sensor must travel from the closest
station. For the subarea $\sigma_{i}$, they gave an algorithm for finding
$w(\sigma_{i})$, the minimum distance from any point in $\sigma_{i}$ to any
station in $P$. Given the representation of all subareas, computing
$w(\sigma_{i})$ for all subareas $\sigma_{i}$ takes only $O(n^{2}k)$.

After all these preprocessing steps, we are left with a combinatorial problem
of choosing a subset of subareas to cover all targets.

Let $Q$ be the bounding box containing all targets. As in many geometric PTAS,
Gao et al. then employed the algorithm, called the Partition algorithm, that
divides $Q$ into cells $cell(Q)$, each of size $2mr\times 2mr$, solves the
selection optimally for each cell $e\in cell(Q)$, and returns the union of
selected subareas.

The essential subproblems for the Partition algorithm deals with each cell
$e$. To solve each subproblem, Gao et al. used a brute-force search that runs
in time $n^{O(m^{2})}$. To bound the running time they proved the following
key lemma (implicitly in Theorem 5).

######  Lemma 3.2 ([GCWC17, CGWC16]).

The number of sensors (and subareas selected) in the optimal solution in each
cell $e$ is at most $O(m^{2})$.

However, this lemma is incorrect (see Appendix A). We provide a fix to this
lemma (in Section 5) that implies that the number of sensors is $O(m^{2}+k)$,
yielding that their algorithm for each subproblem actually runs in time
$n^{O(m^{2}+k)}$. More generally, if we have an upperbound of $L$ on the
number of sensors, a brute-force algorithm would run in time $n^{O(L)}$. We
also provide a stronger bound that allows us to improve the running time to
$n^{O(m+k)}$, improving the running time exponentially.

After solving each subproblem, each solution is combined straight-forwardly.
Gao et al. used shifting techniques to show that if the targets are
distributed uniformly in $Q$ (as described in the model in Section III-A), the
total cost is $(1+3/m)$ times the optimal cost, providing the guarantee of
$1+\epsilon$ of the PTAS. We note that their proof works even without the
assumption on the distribution of the targets, for completeness, we also
provide the proof of this performance guarantee in Section 4.2.

##  4 The Algorithm

We focus on the second phase of EEMA where we work entirely on the
combinatorial version of the problem (after the transformation using the Graph
Conversion algorithm of [GCWC17, CGWC16]). Recall that we are given a set of
$h$ minimal subareas $R=\\{\sigma_{1},\sigma_{2},\ldots,\sigma_{h}\\}$ defined
by intersections of detection cycles. Recall also that $h=O(n^{2})$. Also for
each subarea $\sigma_{i}$, we also know the minimum distance $w(\sigma_{i})$
needed if we schedule a sensor into this subarea. The goal is to select a
subset of subareas $R^{\prime}\subseteq R$ that minimizes the total distance
while ensuring that these subareas cover all targets.

Our main subroutine works with a square region of size $2mr\times 2mr$,
referred to as a cell. In Subsection 4.1, we present our key contribution, an
algorithm for solving the problem optimally when all targets are inside a cell
that runs in time $n_{e}^{O(m+k)}$, where $n_{e}$ is the number of subareas
that can cover targets in $e$. Since later on we set $m=\Theta(1/\epsilon)$,
this implies the total running time of $n^{O(1/\epsilon)}$, an improvement
over the brute-force algorithm of [GCWC17] that runs in time
$n^{O(1/\epsilon^{2})}$.

To solve the whole problem, the surveillance region is partitioned into many
cells, where the dynamic programming can be applied to. Let $Q$ be the
bounding box covering the surveillance region. More specifically, assume that
the lower-left corner of $Q$ is at $(0,0)$ and its upper-right corner is at
$(M,M)$. To facilitate the shifting, we further assume that targets are inside
the square whose corners are $(2mr,2mr)$ and $(M,M)$, i.e., there are no
targets whose $x$ or $y$ co-ordinates are less than $2mr$. We divide $Q$
uniformly into square cells of size $2mr\times 2mr$. There will be $\lceil
M/2mr\rceil\cdot\lceil M/2mr\rceil$ cells. Each cell is a square whose lower-
left and upper-right corners are of the form $(i\cdot 2mr,j\cdot 2mr)$ and
$((i+1)\cdot 2mr,(j+1)\cdot 2mr)$. Let $cell(Q)$ be the set of all cells. See
Figure 2 (a).

(a)

(b)

Figure 2: (a) the input region with targets (shown in shaded area), the
bounding box $Q$, and the set of all cells of size $2mr\times 2mr$. (b) The
shifted bounding box $Q_{f}$ and its cells $cell(Q_{f})$.

The algorithm then applies the dynamic programming subroutine to every cell
$e\in cell(Q)$, and returns the union of the selected subareas. Gao et al.
[GCWC17] proved that the solution returned is $4$-approximate.

Using standard technique for PTAS design, as also in [GCWC17], we shift the
bounding box to “average out” the cost of the optimal to obtain the
$(1+\epsilon)$-approximate solution. We discuss this step in Subsection 4.2
and present the proof of the performance guarantee.

Finally, using the bounds to be proved in Section 5, we give the analysis of
the running time in Subsection 4.3.

###  4.1 Dynamic Programming

In this section, we describe a dynamic programming that, for a given cell $e$
of size $2mr\times 2mr$, finds the optimal schedule for sensors to cover all
the targets in $e$. The cell $e$ is partitioned into vertical strips each of
size $r\times 2mr$ (defined formally later). If we know that the number of
sensors in the optimal solution for each strip of $e$ is at most $L$, the
dynamic programming algorithm in this section runs in time $n_{e}^{O(L)}$,
where $n_{e}$ is the number of subareas containing any targets inside cell
$e$. We give a proof in Lemma 5.1 that $L=O(m+k)$; since we assume that $k$ is
a constant, the running time of the dynamic programming is $n_{e}^{O(m)}$.

Because we use dynamic programming, we shall define other types of
subproblems. To avoid confusion later in this section, we shall refer to this
main subproblem in a cell as the cell problem. We further divide the cell into
$2m$ vertical strips, each of size $r\times 2mr$. We refer to these strips as
strips $H_{1},H_{2},\ldots,H_{2m}$. See Figure 3.

$H_{1}$$H_{2}$$H_{i}$$H_{2m}$$\cdots$$\cdots$ Figure 3: Strips
$H_{1},\ldots,H_{2m}$ and subproblem $A_{i}$ (shown in shade).

We define the following subproblems $\mathcal{A}_{i}$ for each $1\leq i\leq
2m$ as follows.

Subproblem $\mathcal{A}_{i}$: find subset of subareas $R^{\prime}$ that covers
all targets in $H_{1}\cup H_{2}\cup\cdots\cup H_{i}$ with the minimum total
distance.

Let $a_{i}$ denote the optimal total distance for $\mathcal{A}_{i}$. Recall
that with this definition $a_{2m}$ is the optimal cost for the cell problem.
Here we only consider the algorithm that finds the minimum cost; to also find
the actual solution, one only need to keep the best solutions together with
its costs in the dynamic programming states.

For strip $H_{i}$, let $R_{i}\subseteq R$ be subset of subareas that cover
some target in $H_{i}$. For a subset $U\subseteq R_{i}$ of subareas of size at
most $L$ covering some targets in $H_{i}$, we also define subproblems
$\mathcal{B}_{i}(U)$ for each $1\leq i\leq 2m$ as follows.

Subproblem $\mathcal{B}_{i}(U)$: find subset of subareas $R^{\prime}$ that
covers all targets in $H_{1}\cup H_{2}\cup\cdots\cup H_{i}$ such that
$R^{\prime}\cap R_{i}=U$ with the minimum total distance.

Note that $\mathcal{B}_{i}(U)$ considers the same problem as $\mathcal{A}_{i}$
but with additional conditions that to cover targets in $H_{i}$ one must use
subareas in $U$. Let $b_{i,U}$ denote the optimal total distance for
$\mathcal{B}_{i}(U)$. If it is impossible to do so, we let $b_{i,U}=\infty$.
Therefore, we know that

$a_{i}=\min_{U\subseteq R_{i}:|U|\leq L}b_{i,U}.$

Hence, it is enough to solve $\mathcal{B}_{i}(U)$.

As the base case, to solve $\mathcal{B}_{1}(U)$ for each $U\subseteq R_{1}$,
we can easily check if $U$ covers all targets and keep its cost or keep the
value of $\infty$ otherwise. Because $|U|\leq L$, there are at most
$O(n_{e}^{L})$ subproblems to consider. Since it takes $O(n^{3})$ for each
subproblem, the total running time is $n_{e}^{O(L)}$.

To solve $\mathcal{B}_{i}(U)$ for each $U\subseteq R_{i}$, first note that
since the strip width is $r$, any subareas $\sigma^{\prime}$ covering any
targets in $H_{i-2}$ cannot be in $R_{i}$, i.e., they cannot help covering
targets in $H_{i}$. On the other hand, subareas in $R_{i}$ cannot cover any
targets in $H_{1},\ldots,H_{i-2}$. Therefore, when considering solutions, we
only need to worry about dependencies between subareas in $U\in R_{i}$ and
$U^{\prime}\in R_{i-1}$. Thus to find the solution for $\mathcal{B}_{i}(U)$,
we enumerate all $U^{\prime}\in R_{i-1}$ and checks

  * •

if $U$ and $U^{\prime}$ are compatible (i.e., for subareas $\sigma$ in both
$R_{i-1}$ and $R_{i}$, $U$ and $U^{\prime}$ make the same decision, i.e.,
$\sigma\in U$ iff $\sigma\in U^{\prime}$), and

  * •

if $U\cup U^{\prime}$ covers all targets in $H_{i}$.

If both conditions are met, the set of subareas for
$\mathcal{B}_{i-1}(U^{\prime})$ and $U$ form a feasible solution for
$\mathcal{B}_{i}(U)$. Finally, we choose $U^{\prime}$ that minimizes the cost.

For each $U$, there are $O(n_{e}^{L})$ subsets $U^{\prime}$ to consider.
Therefore to solve all $\mathcal{B}_{i}(U)$ for all $U\in R_{i}$, the
algorithm considers $O(n_{e}^{L})\cdot O(n_{e}^{L})=n_{e}^{O(L)}$, resulting
in the running time of $n_{e}^{O(L)}$ as well. The total running time for
solving all $\mathcal{A}_{i}$ is thus $m\cdot n_{e}^{O(L)}=m\cdot
n_{e}^{O(m+k)}$, since $L=O(m+k)$.

###  4.2 Shifting and the Performance Guarantee

The shifting procedure proceeds in $m$ rounds. For each round $f$, for $0\leq
f\leq m-1$, the bounding box $Q$ is translated so that its lower-left corner
is at $(2f,2f)$. Let $Q_{f}$ denote the translated bounding box and
$cell(Q_{f})$ be the cells in $Q_{f}$. Let $S_{f}$ be the set of subareas
returned from round $f$. The shifting procedure selects the solution with the
minimum cost over returned solutions $S_{f}$’s. See Figure 2 (b).

We shall prove that there exists round $f$ such that the cost of $S_{f}$ is
$1+4/m$ of the optimal cost. Note that Gao et al. also presented a similar
proof in Theorem 8 of [GCWC17] (appeared as Theorem 6 in [CGWC16]) but they
did not realize that this proof works for general case without the target
distribution assumption. We include the proof here for completeness.

######  Theorem 4.1.

The shifting procedure returns a solution whose cost is at most $1+4/m$ of the
optimal cost.

###### Proof.

For any set of subareas $S^{\prime}$, let $w(S^{\prime})$ be its cost, i.e,

$w(S^{\prime})=\sum_{\sigma\in S^{\prime}}w(\sigma).$

Let $S^{*}$ be the optimal set of subareas and $w(S^{*})$ be the optimal total
movement cost. Consider round $f$ and the returned solution $S_{f}$. For each
$e\in cell(Q_{f})$, for a set of subareas $S^{\prime}$, let $S^{\prime}(e)$ be
the set of subareas in $S^{\prime}$ that cover some target in $e$. From the
optimality of the dynamic programming algorithm, we have that for any $e$,

$\sum_{\sigma\in S_{f}(e)}w(\sigma)\leq\sum_{\sigma\in S^{*}(e)}w(\sigma).$

For a cell $e$ and a set of subareas $S^{\prime}$, let $S^{\prime}(\delta(e))$
be the set of subareas in $S^{\prime}$ that cover some target in $e$ and also
intersect the boundary of $e$. We can split the cost of $S^{*}(e)$ into two
parts, i.e.,

$\sum_{\sigma\in S^{*}(e)}w(\sigma)=\sum_{\sigma\in S^{*}(e)\setminus
S^{*}(\delta(e))}w(\sigma)+\sum_{\sigma\in S^{*}(\delta(e))}w(\sigma).$

Note that for different $e,e^{\prime}\in cell(Q_{f})$, $S^{*}(e)\setminus
S^{*}(\delta(e))$ and $S^{*}(e^{\prime})\setminus S^{*}(\delta(e^{\prime}))$
are disjoint. Thus summing the cost for $S_{f}$, we have that $w(S_{f})$
equals

$\displaystyle\sum_{e}\sum_{\sigma\in S_{f}(e)}w(\sigma)$
$\displaystyle\leq\sum_{e}\sum_{\sigma\in S^{*}(e)}w(\sigma)$
$\displaystyle=\sum_{e}\sum_{\sigma\in S^{*}(e)\setminus
S^{*}(\delta(e))}w(\sigma)\;+$ $\displaystyle\;\;\;\;\sum_{e}\sum_{\sigma\in
S^{*}(\delta(e))}w(\sigma),$ $\displaystyle\leq
w(S^{*})+\sum_{e}\sum_{\sigma\in S^{*}(\delta(e))}w(\sigma),$

where $e$ is summed over $e\in cell(Q_{f})$, and the last inequality follows
from the disjointness of subareas strictly inside a cell $e$.

Taking the average over $w(S_{f})$ for all $f$, we get

$\displaystyle\frac{1}{m}\sum_{f=0}^{m-1}w(S_{f})$
$\displaystyle\leq\frac{1}{m}\sum_{f=0}^{m-1}\left(w(S^{*})+\sum_{e\in
cell(Q_{f})}\sum_{\sigma\in S^{*}(\delta(e))}w(\sigma)\right)$
$\displaystyle=w(S^{*})+\frac{1}{m}\sum_{f=0}^{m-1}\left(\sum_{e\in
cell(Q_{f})}\sum_{\sigma\in S^{*}(\delta(e))}w(\sigma)\right).$

Consider a particular subarea $\sigma\in S^{*}$. Note that there are at most
$4$ cells $e$ over all shifted instances such that $\sigma\in
S^{*}(\delta(e))$, because $\sigma$ may intersect at most one vertical cell
boundary twice and another one horizontal cell boundary twice. This implies
that the second term is at most

$\sum_{f=0}^{m-1}\left(\sum_{e\in cell(Q_{f})}\sum_{\sigma\in
S^{*}(\delta(e))}w(\sigma)\right)\leq 4w(S^{*}).$

Plugging this inequality into the previous bound, we have that

$\frac{1}{m}\sum_{f=0}^{m-1}w(S_{f})\leq
w(S^{*})+\frac{4}{m}w(S^{*})=(1+4/m)w(S^{*}).$

Since all $w(S_{f})$’s are non-negative, at least one $f$ exists such that
$w(S_{f})\leq(1+4/m)w(S^{*})$, as required. ∎

###  4.3 Running time

From Theorem 4.1, we need to set $m=4/\epsilon$ to guarantee that the
algorithm gives a $(1+\epsilon)$-approximate solution. For each shifted
bounding box $Q_{f}$ and for each cell $e\in cell(Q_{f})$, the dynamic
programming algorithm from Section 4.1 runs in time $m\cdot n_{e}^{O(L)}$
where $n_{e}$ is the number of subareas covering some target in $e$ and $L$ is
the upperbound on the number of sensors needed.

In Section 5, we show that $L=O(m+k)=O(m)$ since we assume that $k$ is a
constant. Therefore, the dynamic programming algorithm runs in time
$n_{e}^{O(m)}$. Combining the running time for the dynamic programming for all
cells in all shifted bounding boxes, we get that the algorithm runs in time

$\sum_{f}\sum_{e\in cell(Q_{f})}m\cdot
n_{e}^{O(m)}=m^{2}h^{O(m)}=m^{2}n^{O(m)}=(1/\epsilon^{2})n^{O(1/\epsilon)},$

as claimed.

##  5 Bounds on the number of sensors

Our key result needed for bounding the running time of the dynamic programming
in Section 4.1 is the following lemma.

######  Lemma 5.1.

For each cell of size $2mr\times 2mr$, the number of sensors in a square of
size $\sqrt{r/2}\times\sqrt{r/2}$ in an optimal solution is at most
$O(1)+O(k)=O(k)$. Moreover, for any vertical strip of size $r\times 2mr$ in a
cell, there exists an optimal solution that employs at most $O(m+k)$ sensors
in this strip.

The following two lemmas are sufficient to establish Lemma 5.1. Lemma 5.2
deals with sensors which are far from any stations, i.e., sensors whose moving
distances are at least $r/2$. Lemma 5.3 considers sensors which stay close to
a station. We remark that we make no attempt to optimize the constants in the
proofs.

######  Lemma 5.2.

For each cell of size $2mr\times 2mr$, the number of sensors in a square $q$
of size $\sqrt{r/2}\times\sqrt{r/2}$ in the optimal solution whose moving
distances are at least $r/2$ is at most a constant.

###### Proof.

Note that sensors in $q$ can cover targets in adjacent squares (all of size
$\sqrt{r/2}\times\sqrt{r/2}$). However since the sensing radius is at most
$r$, they cannot cover targets that much further. We draw additional 20
adjacent squares as shown in Figure 4. Note that moving a sensor inside $q$ to
any of these squares requires moving it for additional distance of $r$; thus
increasing the cost by this much. However, since every sensor considered in
this lemma has already moved for the distance of $r/2$, removing one sensor
reduces the cost by $r/2$.

$q$ Figure 4: Adjacent squares inside which targets can be covered by sensors
in $q$

We shall prove the lemma by contradiction. Let $m_{q}$ be the number of
sensors inside $q$. Assume that there are at least $63$ sensors in $q$, i.e.,
$m_{q}\geq 63$. Consider another solution where we remove $m_{q}-21$ sensors,
arbitrarily, from $q$, and move another $20$ sensors to all $20$ adjacent
squares. Note that the new schedule covers all previously covered targets. On
the one hand, removing sensors decreases the cost by at least $(m_{q}-21)\cdot
r/2\geq(63-21)\cdot r/2=20r$. On the other hand, moving sensors to adjacent
squares increases the cost by at most $20r$. Thus, we obtain a cheaper
solution and this contradicts the assumption. ∎

It remains to analyze sensors that are close to the station. The following
Lemma 5.3 shows that for a station $p\in P$, there exists an optimal solution
that schedules at most $O(1)$ sensors within the distance of $r/2$ from $p$.
We again remind that we do not try to optimize the constant in the lemma.

######  Lemma 5.3.

In an optimal solution, the number of sensors whose moving distances are at
most $r/2$ from a particular station $p_{\ell}$ is at most a constant.

Consider station $p\in P$. Since it is free to place a sensor at $p$, we
assume that all targets within the radius $r$ of $p$ are covered. We start
with an important observation.

######  Lemma 5.4.

For any distance $0<a\leq r/2$, the number of sensors originating from $p$
whose distance from $p$ is in the range $[a/2,a]$ is at most a constant.

###### Proof.

Sensors in this range may cover targets $t$ whose distance $d(p,t)$ satisfies
$r\leq d(p,t)\leq r+a$. Call the area in that distance a ring $Ring(r,r+a)$.
Note that we can place $6$ sensors around a circle of radius $2a\leq r$ to
cover entirely $Ring(r,r+a)$, paying the cost of $12a$. Since each sensor
whose distance to station $p$ falls in the range $[a/2,a]$ must travel a
distance of at least $a/2$, having more than $24$ sensors in this range costs
more than this simple feasible solution, leading to a contradiction. ∎

While the previous lemma is a good step in proving the bound, it can be
applied infinitely (as distances can be made smaller and smaller) resulting in
an unbounded number of sensors. To resolve that we need a notion of progress.
We first give an overview for the proof of Lemma 5.3.

To prove this lemma, we shall apply Lemma 5.4 repeatedly to the optimal
solution. If we can show that the number of times we apply the lemma is
constant, we know that the total number of sensors is also constant. Later in
this section, we mainly state the required preliminary lemmas; their proofs
are deferred to Appendix B.

Consider the furthest sensor from station $p$ in the optimal solution whose
distance from $p$ fall in the range $(0,r/2]$. Let $a_{1}$ be the distance and
denote that sensor by $s_{1}$. Lemma 5.4 ensures that the number of sensors
whose distance from $p$ are in the range of $[a_{1}/2,a_{1}]$ is constant.

To apply Lemma 5.4 again, we look for the furthest sensor which is closer to
$p$ than $a_{1}/2$. Let $a_{2}$ be its distance. We repeat this process,
yielding a sequence of distances $a_{1},a_{2},\ldots$, where

$a_{i+1}\leq a_{i}/2.$

These distances define a sequence of levels of sensors, where $L_{i}$ is the
set of sensors whose distances from $p$ are in the range $(a_{i}/2,a_{i}]$.

To show that the number of levels is at most a constant, we define the
following notion of covered angles to capture the progress. From the station
$p$ located at $(x_{p},y_{p})$, fix a reference line $\ell$ to be a (half)
straight line originating from $(x_{p},y_{p})$ to $(\infty,y_{p})$. For any
target $t$, given $p$ and line $\ell$, we can define the angle $\theta_{t}$ to
be the angle determined by segment $\overline{pt}$ and line $\ell$. We say
that $t$ is at angle $\theta_{t}$ and we also refer to $\theta_{t}$ as the
angle of $t$. We say that the interval $[\theta_{1},\theta_{2}]$ of angles is
covered at radius $a^{\prime}$ if every target at angle
$\theta\in[\theta_{1},\theta_{2}]$ at distance between $r$ and $a^{\prime}$ is
covered.

To be precise, we have to describe how to account for covered angles. When
dealing with sensors at level $L_{i}$, we consider a single sensor $s$ in that
level furthest from $p$. Let $a$ be the distance from $s$ to station $p$. At
this level we consider the interval of angles covered by $p$ at radius $a/2$.
When we consider the next level $L_{i+1}$, we again consider a single sensor
$s^{\prime}\in L_{i+1}$ at furthest distance $a^{\prime}$ from $p$ and try to
account for covered angles at radius $a^{\prime}/2$. Since previously covered
intervals at radius $a/2$ remain covered at radius $a^{\prime}/2<a/2$; we can
take all covered intervals from the previous level to this level as well. We
note that while sensors at lower levels $s$ may cover more angles at radius
$a^{\prime}/2$ (for higher level analysis), we keeps the covered intervals
fixed when we analyze the additional covered angles; thus the intervals of
covered angles act as lowerbounds on the actual angles covered.

Thus, to account for the progress, we maintain a set $S$ of (circular)
disjoint intervals of angles and the current guarantee radius $\alpha$, with
the condition that any target $t$ with distance at most $\alpha$ from station
$p$ at angle $\theta_{t}$ such that $\theta_{t}\in I$ for some interval $I\in
S$ is covered by the current set of sensors. We remark that if this condition
holds for some set $S$ with radius $\alpha$, it also holds for any radius
$\alpha^{\prime}<\alpha$.

The following lemma (whose proof is in Appendix B.1) shows that if we place a
sensor at distance $a$ from $p$, it would cover an interval of angles at
radius $a^{\prime}$ of width at least $2\cdot\arccos(13/20)$, where
$a^{\prime}\leq a/2$. (Also, see Figure 8.)

######  Lemma 5.5.

Assume that station $p$ is at $(0,0)$. If we place sensor $s$ at distance $a$
from $p$, w.l.o.g. at $(a,0)$, where $a\leq r/2$, the interval

$[-\arccos(13/20),\arccos(13/20)]$

is covered at radius $a^{\prime}$, where $a^{\prime}\leq a/2$.

Lemma 5.5 essentially allows us to make progress from one level to the next.
However, we need more structures to ensure that when we have circles from many
levels. Consider two sensors $s$ and $s^{\prime}$ at two consecutive levels
(i.e., $s\in L_{i}$ and $s^{\prime}\in L_{i+1}$). Assume that $s$ is at
distance $a$ from the station $p$ where $a\leq r/2$; this also implies that
$s^{\prime}$ is at distance $a^{\prime}\leq a/2$ from $p$. We consider the
regions covered additional by $s^{\prime}$. Formally, let circles
$C_{p},C_{s},$ and $C_{s^{\prime}}$ denoted circles of radius $r$ centered at
$p,s,$ and $s^{\prime}$ respectively. The regions covered additionally by
$s^{\prime}$ is $C_{s^{\prime}}\setminus(C_{p}\cup C_{s})$, shown in Figure 5.
The following lemma, regarding these regions, is crucial.

$p$$s$$s^{\prime}$$R_{1}$$R_{2}$ Figure 5: Two regions $R_{1}$ and $R_{2}$ in
$C_{s^{\prime}}\setminus(C_{p}\cup C_{s})$. Blue dashed circle shows the
distance $r+a^{\prime}/2$.

######  Lemma 5.6.

Let $s$ and $s^{\prime}$ be sensors at levels $L_{i}$ and $L_{i+1}$ with
distances $a$ and $a^{\prime}$ from station $p$, respectively, where
$a^{\prime}\leq a/2\leq r/4$. Consider $C_{s^{\prime}}\setminus(C_{p}\cup
C_{s})$. If the set contains two disjoint regions $R_{1}$ and $R_{2}$, the
further points in both regions is at distance at most $r+a^{\prime}/2$ from
station $p$.

The proof of Lemma 5.6 appears in Appendix B.2. Given Lemma 5.6, we can argue
that if $C_{s^{\prime}}\setminus(C_{p}\cup C_{s})$ contains two regions, it is
cheaper to place two sensor at distance strictly less than $a^{\prime}/2$ from
$p$; as the cost would be $2\cdot a^{\prime}/2<a^{\prime}$; contradicting the
optimality of the optimal solution. Therefore, the additional covering area of
$C_{s^{\prime}}$ contains only one region. In that case, we can definitely
rotate $C_{s^{\prime}}$ to maximize the additional covering angle at distance
$a^{\prime}/2$ to be at least $\arccos(13/20)$, half of the angle guaranteed
by Lemma 5.5.

Figure 6 shows additional angles covered by $s^{\prime}$.

$p$$s$$s^{\prime}$$v$$w$ Figure 6:  How sensor $s^{\prime}$ (level $L_{i+1}$)
at distance $a^{\prime}=a/2$ from station $p$ covers additional angle interval
$[2\theta,\theta]$ at radius $a^{\prime}/2$ (shown in green), where
$\theta=\arccos(13/20)$.

Let us state the formal statement of the discussion above regarding two
sensors at two consecutive levels. Its proof appears in Appendix B.3.

######  Lemma 5.7.

Consider two sensor $s$ and $s^{\prime}$ at consecutive levels $L_{i}$ and
$L_{i+1}$. Let $a$ and $a^{\prime}$ be the distances from $s$ and $s^{\prime}$
to $p$. Let the interval of angles covered by $s$ at radius $a/2$ be $I$. The
additional angle (w.r.t. $I$) covered by $s^{\prime}$ at radius $a^{\prime}/2$
is at least $\arccos(13/20)$.

We are ready to prove our key Lemma 5.3

###### Proof.

(of Lemma 5.3)

We show that the number of levels is constant. In our analysis we maintain a
set $S$ of disjoint covered angle intervals and the radius $b$.

Recall that $a_{1},a_{2},\ldots,$ are the furthest distances of sensors
$s_{1},s_{2},\ldots,$ at levels $L_{1},L_{2},\ldots$ from station $p$.
Initially, starting at level $L_{1}$, considering the furthest sensor
$s_{1}\in L_{1}$, from Lemma 5.5, we have that at radius $b=a_{1}/2$, $s_{1}$
covers the interval of width $2\cdot\arccos(13/20)$. Thus, we have that $S$
contains a single interval of width $2\cdot\arccos(13/20)$ of covered angles
at radius $b=a_{1}/2$.

As we consider sensors at higher levels, we make progress by either increase
the total width of the intervals or decrease the number of intervals.
Intuitively, we show that while iterating over these levels, one of these
events occurs:

  * (1)

when considering a sensor $s_{i}$ in level $L_{i}$ with distance $a_{i}$, the
sum of the width of the covered intervals at radius $a_{i}/2$ increases by at
least a constant (more precisely at least $\arccos(13/20)\approx 0.863$), or

  * (2)

if $s_{i}$ does not satisfy the previous property, it means that $s_{i}$ is
“between” two other sensors from the lower levels which are well-separated; in
this case the show that the number of covered intervals decreases.

Consider sensor $s_{i}$ at distance $a_{i}$ at level $L_{i}$. We consider the
intervals of covered angles at radius $b=a_{i}/2$. As discussed before every
target whose angle is in any interval in $S$ is covered at the new lower
radius $b$. Let interval $J_{i}$ be the interval of angles covered by $s_{i}$
at radius $a_{i}/2$. Again from Lemma 5.5, we have that $J_{i}$ is of width
$2\cdot\arccos(13/20)$.

We would like to add $J_{i}$ to $S$. If $J_{i}$ only intersects at most one
interval $J^{\prime}$ already in $S$, we know that either the additional angle
is at least $\arccos(13/20)$, or from Lemma 5.7, we can move $s_{i}$ along the
circle of radius $a_{i}$ so that the additional angle after combining $J_{i}$
with $J^{\prime}$ is at least $\arccos(13/20)$; resulting in event (1). The
only case this is not true is when after moving $s_{i}$, $J_{i}$ intersects
two intervals; thus, adding $J_{i}$ decreases the number of intervals in $S$
and event (2) occurs.

Since the number of intervals cannot be less than 1, and each time event (1)
occurs we increase the total width by a constant and introduce at most one new
interval, we know that both events occurs for a constant number of times as
well, because the total width is at most $2\pi$. (To be precise, the number of
levels is at most $\lceil 2\cdot 2\pi/\arccos(13/20)\rceil=15$.) This implies
that the number of levels we encounter is at most constant, yielding the
lemma. ∎

##  6 Conclusions

We present a polynomial-time approximation scheme for finding a
$(1+\epsilon)$-approximate solution running in time $n^{O(1/\epsilon)}$ for
this problem when $k$, the number of base stations, is constant. This
improves, exponentially, over the PTAS presented by Gao et al. [GCWC17] that
runs in time $n^{O(1/\epsilon^{2})}$. We prove a stronger bound on the number
of sensors in any unit area in the optimal solution and employ a more refined
dynamic programming algorithm whose complexity depends only on the width of
the problem. We hope that the dynamic programming approach presented here
might find more applications in other networking problems.

## References

  * [Amm12] Habib M. Ammari.  On the problem of k-coverage in mission-oriented mobile wireless sensor networks.  Computer Networks, 56(7):1935–1950, 2012.
  * [ASSC02] I.F. Akyildiz, W. Su, Y. Sankarasubramaniam, and E. Cayirci.  Wireless sensor networks: a survey.  Computer Networks, 38(4):393 – 422, 2002.
  * [Bak94] Brenda S. Baker.  Approximation algorithms for np-complete problems on planar graphs.  J. ACM, 41(1):153–180, January 1994.
  * [CGWC16] Z. Chen, X. Gao, F. Wu, and G. Chen.  A ptas to minimize mobile sensor movement for target coverage problem.  In IEEE INFOCOM 2016 - The 35th Annual IEEE International Conference on Computer Communications, pages 1–9. IEEE INFOCOM 2016 - The 35th Annual IEEE International Conference on Computer Communications, 2016.
  * [DWW+12] Ling Ding, Weili Wu, James Willson, Lidong Wu, Zaixin Lu, and Wonjun Lee.  Constant-approximation for target coverage problem in wireless sensor networks.  In 2012 Proceedings IEEE INFOCOM, pages 1584–1592, 2012.
  * [EM10] Thomas Erlebach and Matúš Mihalák.  A (4 + $\epsilon$)-approximation for the minimum-weight dominating set problem in unit disk graphs.  In Evripidis Bampis and Klaus Jansen, editors, Approximation and Online Algorithms, pages 135–146, Berlin, Heidelberg, 2010. Springer Berlin Heidelberg.
  * [ETYH18] Mohamed Elhoseny, Alaa Tharwat, Xiaohui Yuan, and Aboul Ella Hassanien.  Optimizing k-coverage of mobile wsns.  Expert Systems with Applications, 92:142–153, 2018.
  * [GCWC17] Xiaofeng Gao, Zhiyin Chen, Fan Wu, and Guihai Chen.  Energy efficient algorithms for $k$-sink minimum movement target coverage problem in mobile sensor network.  IEEE/ACM Trans. Netw., 25(6):3616–3627, dec 2017.
  * [GJ19] Jun Guo and Hamid Jafarkhani.  Movement-efficient sensor deployment in wireless sensor networks with limited communication range.  IEEE Transactions on Wireless Communications, 18(7):3469–3484, 2019\.
  * [HCL+12] Shibo He, Jiming Chen, Xu Li, Xuemin Shen, and Youxian Sun.  Cost-effective barrier coverage by mobile sensor networks.  In 2012 Proceedings IEEE INFOCOM, pages 819–827, 2012.
  * [HMR+98] Harry B Hunt, Madhav V Marathe, Venkatesh Radhakrishnan, S.S Ravi, Daniel J Rosenkrantz, and Richard E Stearns.  Nc-approximation schemes for np- and pspace-hard problems for geometric graphs.  Journal of Algorithms, 26(2):238 – 274, 1998.
  * [KJ17] Thamarak Khampeerpat and Chaiporn Jaikaeo.  Mobile sensor relocation for nonuniform and dynamic coverage requirements.  IEICE Trans. Inf. Syst., 100-D(3):520–530, 2017.
  * [LDNT13] Benyuan Liu, Olivier Dousse, Philippe Nain, and Don Towsley.  Dynamic coverage of mobile sensor networks.  IEEE Transactions on Parallel and Distributed Systems, 24(2):301–311, 2013.
  * [LDWS08] Benyuan Liu, Olivier Dousse, Jie Wang, and Anwar Saipulla.  Strong barrier coverage of wireless sensor networks.  MobiHoc ’08, page 411–420, New York, NY, USA, 2008. Association for Computing Machinery.
  * [LSC21] Dieyan Liang, Hong Shen, and Lin Chen.  Maximum target coverage problem in mobile wireless sensor networks.  Sensors, 21(1), 2021.
  * [MBH+95] M. V. Marathe, H. Breu, H. B. Hunt, S. S. Ravi, and D. J. Rosenkrantz.  Simple heuristics for unit disk graphs.  Networks, 25(2):59–68, 1995.
  * [MMA+14] Hamid Mahboubi, Kaveh Moezzi, Amir G. Aghdam, Kamran Sayrafian-Pour, and Vladimir Marbukh.  Distributed deployment algorithms for improved coverage in a network of wireless mobile sensors.  IEEE Transactions on Industrial Informatics, 10(1):163–174, 2014\.
  * [NH05] Tim Nieberg and Johann Hurink.  A ptas for the minimum dominating set problem in unit disk graphs.  In WAOA, pages 296–306. Springer, 2005.
  * [O’R87] Joseph O’Rourke.  Art Gallery Theorems and Algorithms.  Oxford University Press, Inc., USA, 1987.
  * [vL05] Erik Jan van Leeuwen.  Approximation algorithms for unit disk graphs.  In Dieter Kratsch, editor, Graph-Theoretic Concepts in Computer Science, 31st International Workshop, WG 2005, Metz, France, June 23-25, 2005, Revised Selected Papers, volume 3787 of Lecture Notes in Computer Science, pages 351–361. Springer, 2005.
  * [vL06] Erik Jan van Leeuwen.  Better approximation schemes for disk graphs.  In Lars Arge and Rusins Freivalds, editors, Algorithm Theory – SWAT 2006, pages 316–327, Berlin, Heidelberg, 2006. Springer Berlin Heidelberg.
  * [vL09] Erik Jan van Leeuwen.  Optimization and approximation on systems of geometric objects.  PhD thesis, University of Amsterdam, 2009.
  * [Wan11] Bang Wang.  Coverage problems in sensor networks: A survey.  ACM Comput. Surv., 43(4):32:1–32:53, October 2011.
  * [WHT08] You-Chiun Wang, Chun-Chi Hu, and Yu-Chee Tseng.  Efficient placement and dispatch of sensors in a wireless sensor network.  IEEE Transactions on Mobile Computing, 7(2):262–274, 2008.
  * [WLM09] Bang Wang, Hock Beng Lim, and Di Ma.  A survey of movement strategies for improving network coverage in wireless sensor networks.  Computer Communications, 32(13):1427–1436, 2009.
  * [WT08] You-Chiun Wang and Yu-Chee Tseng.  Distributed deployment schemes for mobile wireless sensor networks to ensure multilevel coverage.  IEEE Transactions on Parallel and Distributed Systems, 19(9):1280–1294, 2008.
  * [WZLD20] Weili Wu, Zhao Zhang, Wonjun Lee, and Ding-Zhu Du.  Optimal Coverage in Wireless Sensor Networks.  Number 978-3-030-52824-9 in Springer Optimization and Its Applications. Springer, December 2020.
  * [XS14] Xiaohua Xu and Min Song.  Restricted coverage in wireless networks.  In IEEE INFOCOM 2014 - IEEE Conference on Computer Communications, pages 558–564, 2014.
  * [YMG08] Jennifer Yick, Biswanath Mukherjee, and Dipak Ghosal.  Wireless sensor network survey.  Computer Networks, 52(12):2292 – 2330, 2008.
  * [ZC04] Yi Zou and Krishnendu Chakrabarty.  Sensor deployment and target localization in distributed sensor networks.  ACM Trans. Embed. Comput. Syst., 3(1):61–91, feb 2004.
  * [ZGWD09] Zhao Zhang, Xiaofeng Gao, Weili Wu, and Ding-Zhu Du.  A ptas for minimum connected dominating set in 3-dimensional wireless sensor networks.  Journal of Global Optimization, 45(3):451–458, Nov 2009.
  * [ZWX+11] Feng Zou, Yuexuan Wang, Xiao-Hua Xu, Xianyue Li, Hongwei Du, Pengjun Wan, and Weili Wu.  New approximations for minimum-weighted dominating sets and minimum-weighted connected dominating sets on unit disk graphs.  Theor. Comput. Sci., 412(3):198–208, January 2011.

##  Appendix A A Counter Example to the Sensor Bound in Gao et al.[GCWC17]

The algorithm proposed by Gao et al. [GCWC17] (first appeared as [CGWC16])
relies on the fact that the number of sensors in the optimal solution required
for each cell $e$ of size $2mr\times 2mr$ is at most $O(m^{2})$. Note that
this bound is independent of $k$, the number of stations. The proof presented
in the paper uses the fact that there exists a feasible solution with that
many sensors. However, the actual cost of $k$-MMTC is the total distance, not
the number of sensors. In this section, we provide a counter example to the
sensor bound in Theorem 5 in [GCWC17] that requires $k$ sensors.

$\alpha$$\alpha{+}\beta$$r$$c$$t_{1}$$t_{2}$$t_{k}$$p_{1}$$p_{2}$$p_{k}$
Figure 7: A counter example that requires $k$ sensors.

We assume that $m$ is large enough, which is true when the error requirement
$\epsilon$ is required to be small enough as $m=\Omega(1/\epsilon)$. Let $c$
be the point at the center of cell $e$. For any $\alpha>0$ and very small
$\beta>0$, place $k$ targets uniformly at distance $\alpha$ from $c$ and for
each target $t_{i}$ place a station $p_{i}$ at distance $r+\alpha+\beta$ from
$c$ in the same direction from $c$ as $t_{i}$. Figure 7 illustrates this
example. Consider a solution that for each target $t_{i}$ we move a sensor
from $p_{i}$ with distance $\beta$ to cover it, resulting in the cost of
$k\beta$ with $k$ sensors. If one wants to use smaller number of sensors, some
sensor is required to moved further with distance that depends on $\alpha$ and
$k$; thus we can always choose $\beta$ to be small enough so that that extra
distance is larger than $k\beta$. This ensures that the optimal solution needs
at least $k$ sensors as claimed.

##  Appendix B Proofs of lemmas

###  B.1 Proof of Lemma 5.5

See 5.5

###### Proof.

Let $s$ the sensor at distance $a$ from $p$. Consider two circles $C_{p}$ and
$C_{s}$ centered at $p$ and $s$ with radius $r$. Let $C^{\prime}$ be another
circle centered at $p$ of radius $r+a^{\prime}\leq r+a/2$. Since targets
within distance $r+a^{\prime}$ from $p$ in $C_{s}$ and $C_{p}$ are covered, an
uncovered target must be in $C^{\prime}\setminus(C_{s}\cup C_{p})$. The worst
case would attain when $C^{\prime}$ is largest, so we assume that
$a^{\prime}=a/2$. Also, note that the covered interval would be smallest when
$a$ is largest; thus, we also assume that $a=r/2$. This case is shown in
Figure 8.

$r+a^{\prime}$$r$$a$$q$$p$$s$$v$ Figure 8:  How a single sensor at $(a,0)$
covers targets in angle interval $[-\theta,\theta]$. The black circles are
$C_{p}$ and $C_{s}$; the dashed circle is $C^{\prime}$. Covered area
considered in the lemma is shaded.

We calculate the angle $\theta$ of $vps$ which is half the angle covered. From
the assumptions, we know that the length of $\overline{pv}$ is
$r+a^{\prime}=r+r/4=5r/4$, the length of $\overline{sv}$ is $r$, and the
length of $\overline{ps}=a=r/2$. Let $q$ be the point on $\ell$ such that
angle $pqv$ is the right angle. Elementary calculation shows that the length
of $\overline{sq}$ = $5r/16$. Thus, the angle $\theta$ equals

$\arccos\left(\frac{r/2+5r/16}{5r/4}\right)=\arccos(13/20).$

Since the covered interval is
$[-\theta,\theta]=[-\arccos(13/20),\arccos(13/20)]$, the lemma follows. ∎

###  B.2 Proof of Lemma 5.6

See 5.6

###### Proof.

To prove this lemma, we assume that $C_{s^{\prime}}\setminus(C_{p}\cup C_{s})$
contains two regions $R_{1}$ and $R_{2}$ and there is some point $t_{1}$ in
the regions at distance at least $r+a^{\prime}/2$ from station $p$. Without
loss of generality, we assume that region $R_{1}$ is the region of points with
positive $y$ co-ordinates and $R_{2}$ is the region of points withs negative
$y$ co-ordinates. We further assume that $t_{1}\in R_{1}$.

Let $\delta=a^{\prime}/2\leq a/4$. Given this configuration, we note that we
can certainly move $s^{\prime}$ such that (1) $t_{1}$ is at the distance
exactly at $r+a^{\prime}/2=r+\delta$ from $p$ and also lies on the boundary of
$C_{s}$, and (2) region $R_{2}$ contains exactly one point $t_{2}$ at the
intersection of the boundary of $C_{p}$ and $C_{s}$. Furthermore, this can be
done without increasing the distance from $s^{\prime}$ to $p$. Figure 9
illustrates this assumption where $\delta=a^{\prime}/2$. We note that
$s^{\prime}$ is not shown.

$r{+}\delta$$r$$r$$r$$a$$p$$s$$t_{1}$$t_{2}$ Figure 9: Points $t_{1}\in R_{1}$
and $t_{2}\in R_{2}$ in our proofs.

Our goal is to show that for any point $s^{\prime}$ satisfying the assumption,
$s^{\prime}$ must be too far from $p$; more precisely,

$d(p,s^{\prime})>2\delta=a^{\prime};$

thus, we have the contradiction as required.

Throughout this section, we consider a simple quadrilateral $\square
pt_{2}st_{1}$ (see Figure 9). In our analysis, we let $p$ be at location
$(0,0)$ and $s$ be at $(a,0)$. If we fix $t_{1}$ and $t_{2}$, the point
$s^{\prime}$ must be in $C_{t_{1}}\cup C_{t_{2}}$, to be refered to later on
as the feasible region $R_{f}=C_{t_{1}}\cup C_{t_{2}}$. Without loss of
generality, we assume that $s^{\prime}$ be the point in $R_{f}$ closest to
$p$. (See Figure 10.) The following lemma states the location of $s^{\prime}$.

######  Lemma B.1.

Point $s^{\prime}$ is at

$\left(\frac{\delta^{2}+2r\delta+a^{2}}{2a}-\frac{a}{2},\frac{\sqrt{((2r+\delta)^{2}-a^{2})(a^{2}-\delta^{2})}}{2a}-\frac{\sqrt{(2r)^{2}-a^{2}}}{2}\right).$

###### Proof.

The location for $t_{2}$ is at

$\left(\frac{a}{2},-\frac{\sqrt{(2r)^{2}-a^{2}}}{2}\right),$

and the location for $t_{1}$ is at

$\left(\frac{\delta^{2}+2r\delta+a^{2}}{2a},\frac{\sqrt{((2r+\delta)^{2}-a^{2})(a^{2}-\delta^{2})}}{2a}\right).$

Let $m$ be a midpoint between $t_{2}$ and $t_{1}$, the coordinate for $m$ is
the average of coordinates for $t_{1}$ and $t_{2}$. We remark that $s\in
R_{f}$. Furthermore, $s$ is the furthest point in $R_{f}$ from $p$. Thus, the
point $s^{\prime}$ is the reflection of $s$ on the point $m$; its coordinate
can be calculated to be as stated given that $s$ is at $(a,0)$. ∎

$\theta$$R_{f}$$r{+}\delta$$r$$r$$r$$a$$p$$s$$t_{1}$$t_{2}$$s^{\prime}$ Figure
10: The feasible region $R_{f}$ and $s^{\prime}$.

With Lemma B.1, considering only the $x$ coordinates of $s^{\prime}$, we have
that

$\displaystyle\frac{\delta^{2}+2r\delta+a^{2}}{2a}-\frac{a}{2}$
$\displaystyle=\frac{\delta^{2}+2r\delta+a^{2}-a^{2}}{2a}$
$\displaystyle=\frac{\delta^{2}+2r\delta}{2a}$
$\displaystyle\geq\frac{\delta^{2}+2r\delta}{r}$
$\displaystyle=\frac{\delta^{2}}{r}+2\delta$ $\displaystyle>2\delta,$

where the first inequality follows from the fact that $a\leq r/2$ and the
second strict inequality is true because $\delta>0$. Thus,
$d(p,s^{\prime})>2\delta$ as required.

∎

###  B.3 Proof of Lemma 5.7

See 5.7

###### Proof.

Let $C_{s}$ and $C_{s^{\prime}}$ be two coverage circles for sensor $s$ and
$s^{\prime}$, and $C_{p}$ be the circle of radius $r$ centered at $p$. Lemma
5.6 implies that if the area $C_{s^{\prime}}\setminus(C_{s}\cup C_{p})$
contains two regions, they can be covered by two circles located at distance
at most $a^{\prime}/2$ from $p$, with the total cost of less than
$a^{\prime}$; but this leads to a contradiction.

Thus, $C_{s^{\prime}}\setminus(C_{s}\cup C_{p})$ contains only on region. Let
$v$ be the intersection of $C_{s^{\prime}}$ and $C_{s}$ outside $C_{p}$. (See
Figure 6.) We claim that the additional angle covered by $C_{s^{\prime}}$ is
at least $\arccos(13/20)$. To see this, assume otherwise. In this case, the
distance from $v$ to $p$ is less than $r+a^{\prime}$; thus we can move
$s^{\prime}$ along the circle of radius $a^{\prime}$ centered at $p$ with the
same movement cost while increasing the coverage area outside $C_{s}\cup
C_{p}$. ∎


# Optimal Structured Matrix Approximation for Robustness to Incomplete
Biosequence Data

Chris Salahub, Jeffrey Uhlmann Chris Salahub is a PhD student in Statistics at
the University of Waterloo, Jeffrey Uhlmann is an Associate Professor in
Electrical Engineering and Computer Science at the University of Missouri.

###### Abstract

We propose a general method for optimally approximating an arbitrary matrix
$\mathbf{M}$ by a structured matrix $\mathbf{T}$ (circulant, Toeplitz/Hankel,
etc.) and examine its use for estimating the spectra of genomic linkage
disequilibrium matrices. This application is prototypical of a variety of
genomic and proteomic problems that demand robustness to incomplete
biosequence information. We perform a simulation study and corroborative test
of our method using real genomic data from the Mouse Genome Database [1]. The
results confirm the predicted utility of the method and provide strong
evidence of its potential value to a wide range of bioinformatics
applications. Our optimal general matrix approximation method is expected to
be of independent interest to an even broader range of applications in applied
mathematics and engineering.

###### Index Terms:

Structured matrix, circulant matrix, Toeplitz matrix, Frobenius norm, optimal
approximation, genetics, linkage diseqbuilibrium.

##  I Introduction

Large matrices are ubiquitous to Big Data applications, especially in
bioinformatics, to encode structural relationships among components of complex
systems. These relationships may be physical in the case of connectivities
among molecular components of genomic or proteomic sequences [2, 3, 4], or
they may represent statistical relationships encoded in covariance matrices
[5, 6, 7, 8]. Matrix encodings permeate virtually all areas of science and
engineering because they admit mathematical operations that can reveal subtle
but important properties of a system that otherwise would be practically
impossible to discern. Unfortunately, the size of matrices that arise in many
bioinformatics applications are extremely large, ranging from $n\times n$
matrices with $n$ in the thousands up to tens of millions [9]. This introduces
enormous challenges both because virtually all nontrivial matrix operations
have $O(n^{3})$ complexity111Theoretical algorithms [10, 11] exist that reduce
the exponent from 3 to as low as 2.36 [12], but the improved theoretical
complexity does not translate to improved practical efficiency even for the
largest matrices encountered in current applications [13]. and uncertainty in
data collection frequently makes these matrices incomplete.

In many cases, the computational and statistical costs associated with these
large matrices are prohibitive unless practical approximations can be found.
The most natural and common approach is to simply perform analysis on a small
subset of the data, e.g., by global downsampling [14]. This approach may be
effective in applications involving statistics that will be preserved, but at
a cost of coarser resolution in the form of larger error variances. However,
subsampling cannot be performed for data points corresponding to critical
features that thus cannot be entirely discarded. In this case, the only option
may be to transform the dataset to one of lower complexity that retains
approximate, though degraded, information about all features [15]. An example
would be reduction of a large transportation network by removing roads for
which alternative routes are known to exist so that effective route planning
may be performed, though with potentially suboptimal results.

In this paper we describe a general mechanism for approximating an arbitrary
$n\times n$ matrix, i.e., having $O(n^{2})$ parametric dimensionality, with a
structured $n\times n$ matrix which is defined by as few as $O(n)$ parameters.
In certain cases, this comes with enormous computational and memory savings
for most matrix operations, e.g., spectral decompositions are reduced from
$O(n^{3})$ to $\tilde{O}(n)$, where $\tilde{O}$ hides logarithmic and
polylogarithmic factors. For example, all common matrix operations performed
on a circulant matrix have $O(n\log(n))=\tilde{O}(n)$ time complexity [16], as
opposed to $O(n^{3})$ for a general matrix. The first key question is whether
optimal structured approximations for arbitrary matrices can be constructed
efficiently. The second key question is whether the resulting approximations
retain properties relevant to practical applications. The first question is
answered in the following section, and the second question is subsequently
addressed in Section III with a study involving genomic analysis using data
from the Mouse Genome Database [1].

##  II Structured Matrix Approximation

Suppose we would like to approximate the $n\times n$ matrix

$\mathbf{M}=\begin{bmatrix}m_{00}&m_{01}&\dots&m_{0,n-1}\\\
m_{10}&m_{11}&\dots&m_{1,n-1}\\\ \vdots&\vdots&\ddots&\vdots\\\
m_{n-1,0}&m_{n-1,1}&\dots&m_{n-1,n-1}\end{bmatrix}$ (1)

$m_{ij}\in\mathbb{C}$ with a structured matrix
$\mathbf{T}\in\mathbb{C}^{n\times n}$ for computational or analytical reasons,
e.g., circulant $\mathbf{T}$ for preconditioning [17, 18] and Toeplitz-Hankel
$\mathbf{T}$ for physical modelling [19]. For any application, we prefer the
approximating matrix $\mathbf{T}$ to be optimal by some measure, commonly the
Frobenius norm of the difference $\mathbf{T}-\mathbf{M}$ defined as

$||{\mathbf{T}-\mathbf{M}}||_{F}=\sqrt{\operatorname{trace}\left({(\mathbf{T}-\mathbf{M})}^{\ast}(\mathbf{T}-\mathbf{M})\right)},$
(2)

where ${\mathbf{A}}^{\ast}$ is the conjugate matrix of
$\mathbf{A}\in\mathbb{C}^{n\times n}$. Minimizing this quantity for a
circulant $\mathbf{T}$ was the express goal of [17] and was noted as a
positive feature of the approximation in [18]. Here we present a general
result that can construct an optimal structured matrix approximation to an
arbitrary given general matrix.

###  II-A Structured matrices

We say $\mathbf{T}\in\mathbb{C}^{n\times n}$ is a structured matrix if it has
entries $t_{ij}$, $0\leq i,j\leq n-1$, following a pattern in $i$ and $j$,
that is if

$t_{ij}=t_{f(i,j)}$ (3)

where $f:\\{0,1,\dots,n-1\\}^{2}\mapsto\\{0,1,2,\dots,K\\}$ is the index
function defining the membership of the index pair $i,j$ to an index set with
a constant value. This implies a $k^{\text{th}}$ index set
$\mathcal{T}_{k}=\\{(i,j)|f(i,j)=k\\}$ with cardinality
$\lvert\mathcal{T}_{k}\rvert=n_{k}>0$.

The index function $f(\cdot,\cdot)$ defines the structure of $\mathbf{T}$ by
defining the $\mathcal{T}_{k}$. Common structures and corresponding index
functions are shown in Table I, though these functions are not unique for a
given structure. Many candidate functions define identical index sets, for
example Hankel matrices can take either $f(i,j)=j+i$ or $f(i,j)=2(n-1)-j-i$.

TABLE I: Some common examples of structured index functions. Structure | $f(i,j)$
---|---
Circulant | $(i-j)\bmod n$
Toeplitz | $j-i+n$
Hankel | $i+j$

###  II-B Optimizing the Frobenius norm

Using the notation defined above, we obtain the following theorem.

######  Theorem 1 (Means minimize $||{\mathbf{T}-\mathbf{M}}||_{F}$).

The Frobenius optimal approximating structured matrix $\mathbf{T}$ with index
function $f(i,j)$ for $\mathbf{M}$ is given by $\mathbf{T}_{M}$ with

$t_{ij}=t_{f(i,j)}=\overline{m}_{f(i,j)}$ (4)

where

$\overline{m}_{k}:=\frac{1}{n_{k}}\sum_{\mathcal{T}_{k}}m_{ij}.$ (5)

is the mean of entries in $\mathbf{M}$ over the corresponding index set.
Furthermore, $\frac{1}{\sqrt{n}}||{\mathbf{T}_{M}-\mathbf{M}}||_{F}$ is the
total within-group standard deviation of entries in $\mathbf{M}$ over all
index sets.

###### Proof.

Take $\overline{m}_{k}$ to be the mean of entries in $\mathbf{M}$ for the
$k^{\text{th}}$ index set as in Equation 5, define the vector of all such
means

$\overline{\mathbf{m}}={(\overline{m}_{0},\overline{m}_{1},\dots,\overline{m}_{K})}^{\mkern-1.5mu\mathsf{T}}.$

Further, denote the vector of unique $t_{k}$ as

$\mathbf{t}={(t_{0},t_{1},\dots,t_{K})}^{\mkern-1.5mu\mathsf{T}}$

and the diagonal matrix of $n_{k}$ as

$\mathbf{N}=\operatorname{diag}(n_{0},n_{1},\dots,n_{K}).$

As Equation 2 is always positive, any $\mathbf{T}$ which minimizes
$||{\mathbf{T}-\mathbf{M}}||_{F}$ will also minimize
$||{\mathbf{T}-\mathbf{M}}||_{F}^{2}$. Expanding gives

$\displaystyle\operatorname{trace}\left({(\mathbf{T}-\mathbf{M})}^{\ast}(\mathbf{T}-\mathbf{M})\right)=$
$\displaystyle\operatorname{trace}{\mathbf{M}}^{\ast}\mathbf{M}$
$\displaystyle-\operatorname{trace}{\mathbf{M}}^{\ast}\mathbf{T}$
$\displaystyle-\operatorname{trace}{\mathbf{T}}^{\ast}\mathbf{M}$
$\displaystyle+\operatorname{trace}{\mathbf{T}}^{\ast}\mathbf{T}.$ (6)

${\mathbf{M}}^{\ast}\mathbf{M}$ is constant in $\mathbf{T}$, so can be
ignored. The latter three terms can be considered individually to give
$\operatorname{trace}{\mathbf{T}}^{\ast}\mathbf{T}=\sum_{k=0}^{K}n_{k}{t}^{\ast}_{k}t_{k}$,
$\operatorname{trace}{\mathbf{M}}^{\ast}\mathbf{T}=\sum_{k=0}^{K}n_{k}t_{k}{\overline{m}}^{\ast}_{k}$,
and
$\operatorname{trace}{\mathbf{T}}^{\ast}\mathbf{M}=\sum_{k=0}^{K}n_{k}{t}^{\ast}_{k}\overline{m}_{k}.$
So we seek to minimize

$F(\mathbf{t})=\sum_{k=0}^{K}n_{k}{t}^{\ast}_{k}t_{k}-\sum_{k=0}^{K}n_{k}{t}^{\ast}_{k}\overline{m}_{k}-\sum_{k=0}^{K}n_{k}t_{k}{\overline{m}}^{\ast}_{k},$

which we can write in matrix form as

$\displaystyle F(\mathbf{t})$
$\displaystyle={\mathbf{t}}^{\ast}\mathbf{N}\mathbf{t}-{\mathbf{t}}^{\ast}\mathbf{N}\overline{\mathbf{m}}-{\overline{\mathbf{m}}}^{\ast}\mathbf{N}\mathbf{t}$
$\displaystyle={\left(\mathbf{t}-\overline{\mathbf{m}}\right)}^{\ast}\mathbf{N}\left(\mathbf{t}-\overline{\mathbf{m}}\right)-{\overline{\mathbf{m}}}^{\ast}\mathbf{N}\overline{\mathbf{m}}.$

As $n_{k}>0$ for all $k=0,1,\dots,K$, $\mathbf{N}$ is positive definite, and
so the quadratic form ${\mathbf{x}}^{\ast}\mathbf{N}\mathbf{x}$ has a minimum
of zero when $\mathbf{x}=\mathbf{0}$. Therefore $F(\mathbf{t})$ is minimized
for $\mathbf{t}=\overline{\mathbf{m}}$ and has a minimum of

$F(\overline{\mathbf{m}})=-{\overline{\mathbf{m}}}^{\ast}\mathbf{N}\overline{\mathbf{m}}=-\sum_{k=0}^{K}n_{k}||{\overline{m}_{k}}||^{2}.$
(7)

So $\mathbf{T}_{M}$ is the Frobenius-optimal structured matrix $\mathbf{T}$
approximating $\mathbf{M}$. The residual $\mathbf{T}_{M}-\mathbf{M}$ has a
squared Frobenius norm of

$\displaystyle||{\mathbf{T}_{M}-\mathbf{M}}||_{F}^{2}$
$\displaystyle=\sum_{k=0}^{K}n_{k}\left(\sum_{\mathcal{T}_{k}}\frac{||{m_{ij}}||^{2}}{n_{k}}-||{\overline{m}_{k}}||^{2}\right)$
$\displaystyle=\sum_{k=0}^{K}n_{k}\sigma_{k}^{2}$ (8)

where
$\sigma_{k}^{2}=\frac{1}{n_{k}}\sum_{\mathcal{T}_{k}}\left(m_{ij}-\overline{m}_{k}\right)^{2}$
is the variance of the $m_{ij}$ in the index set $\mathcal{T}_{k}$. Therefore
we have

$\frac{1}{\sqrt{n}}||{\mathbf{T}_{M}-\mathbf{M}}||_{F}=\sqrt{\sum_{k=0}^{K}\frac{n_{k}}{n}\sigma_{k}^{2}},$

which is the total within-group standard deviation from the structured means.
∎

The above proof applies to any matrix with structured equalities, and it
implies a reduction in its space complexity from $O(n^{2})$ to
$O(\lvert{\mathbf{t}}\rvert)$. In the particular case of a circulant matrix,
this reduction achieves an $O(n)$ representation, and the optimal structured
matrix can be computed using the discrete Fourier transform.

###  II-C Circulant matrices

Circulant matrices, denoted $\mathbf{C}$, have an index function

$f(i,j)=(i-j)\bmod n$ (9)

and so contain $n$ unique values denoted
$t_{0},t_{1},\dots,t_{n-1}\in\mathbb{C}$. They see widespread use in signal
processing, computation, and physical modelling both due to their close
relationship with the Fourier transform and their known eigensystem [17, 16,
19]. $\mathbf{C}$ has eigenvalues $\lambda_{k}$ for $k=0,1,\dots,n-1$ given by
$\lambda_{k}=t_{0}+\sum_{l=1}^{n-1}t_{l}\omega^{lk}$ and a corresponding
$k^{\text{th}}$ eigenvector given by
$\mathbf{x}_{k}={(1,\omega^{k},\omega^{2k},\dots,\omega^{(M-1)k})}^{\mkern-1.5mu\mathsf{T}}$
where $\omega=\exp(\frac{2\pi i}{n})$ is the complex $n^{\text{th}}$ root of
unity and $i=\sqrt{-1}$.

Much of the utility of circulant matrices arises from this eigensystem. The
$n\times n$ matrix of eigenvectors of $\mathbf{C}$ scaled to be unitary,

$\mathbf{F}=\frac{1}{\sqrt{n}}[\mathbf{x}_{0}|\mathbf{x}_{1}|\mathbf{x}_{2}|\dots|\mathbf{x}_{n-1}]={\mathbf{F}}^{\mkern-1.5mu\mathsf{T}},$
(10)

is simply the discrete Fourier transform (DFT), which provides an alternate
route to compute $\mathbf{T}_{M}$ in the circulant case. Consider the simple
approximation algorithm in 1.

Algorithm 1 Optimal circulant approximation.

OPTCIRC$(\mathbf{M})$

$\mathbf{D}\leftarrow\operatorname{diag}\left(\mathbf{F}\mathbf{M}{\mathbf{F}}^{\ast}\right)$

$\mathbf{C}_{D}={\mathbf{F}}^{\ast}\mathbf{D}\mathbf{F}$

return $\mathbf{C}_{D}$

As $\mathbf{D}$ is diagonal and $\mathbf{F}$ is the matrix of circulant
eigenvectors, $\mathbf{C}_{D}$ is circulant with eigenvalues $d_{jj}$, the
diagonal values of $\mathbf{D}$. To determine the elements
$(\mathbf{C}_{D})_{ij}$ in terms of $\omega$, $\mathbf{x}$, and $\mathbf{M}$,
first note

$d_{jj}=\frac{1}{n}\sum_{k=1}^{n}\sum_{l=1}^{n}\omega^{(k-l)(j-1)}m_{lk}.$
(11)

Analogously, taking ${\mathbf{F}}^{\ast}\mathbf{D}\mathbf{F}$ gives an $i,j$
element

$\displaystyle(\mathbf{C}_{D})_{ij}$ $\displaystyle=$
$\displaystyle{\mathbf{F}}^{\ast}\mathbf{D}\mathbf{F}$ (12) $\displaystyle=$
$\displaystyle\frac{1}{n^{2}}\sum_{l=1}^{n}\sum_{k=1}^{n}m_{lk}{\mathbf{x}}^{\ast}_{k-l}\mathbf{x}_{i-j}$
$\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{l=1}^{n}\sum_{k=1}^{n}m_{lk}I\big{(}i-j\equiv
k-l\bmod n\big{)}$

where $I(A)$ is the indicator function which returns $1$ if $A$ is true and
$0$ if $A$ is false. As $n_{k}=n$ for all index sets in the circulant case,
Equation 12 indicates that $\mathbf{C}_{D}$ is generated by replacing the
values of $M$ along each circulant diagonal by the corresponding diagonal
mean. Therefore $\mathbf{C}_{D}=\mathbf{T}_{M}$ when $\mathbf{T}$ is
restricted to be circulant and so by Theorem 1 $\mathbf{C}_{D}$ is Frobenius
optimal.

Though related results have been noted for the case where $\mathbf{M}$ is
Toeplitz by [17] and more generally in [18], both of these are particular
examples of the more general result of Theorem 1. Specifically, Algorithm 1
makes no assumptions about $\mathbf{M}$ and can be applied to any matrix.

##  III Application to Genetic Linkage Disequilibrium

In modern genome-wide association studies (GWAS), assessing the impact of a
region of the genome on a measured trait requires an adjustment for linkage
disequilibrium (LD) [20, 21]. Briefly, LD is the correlation between the
inheritance of genetic regions, or markers, over successive generations. For
regions $i$ and $j$ on chromosomes $c_{i}$ and $c_{j}$, $LD(i,j)$ is
frequently computed from an additive measure of genetic distance $d(i,j)$
(often centiMorgans or cM) using the equation

$LD(i,j)=I(c_{i}=c_{j})\gamma e^{-\frac{d(i,j)}{50}}$ (13)

where $\gamma$ is a constant determined by the population characteristics
[22]. As $LD(i,j)$ is constant in $d(i,j)$, this implies the LD matrix is a
structured matrix for $c_{i}=c_{j}$, and zero otherwise.

In practice, LD matrices are used in multiple testing adjustment to account
for correlations between individual tests of markers [23, 24]. Occasionally,
missing data complicate direct use of the LD matrix by forcing pairwise
computation of the LD matrix, which is therefore not guaranteed to be positive
definite and can lead to negative eigenvalues if used directly. Certain
multiple testing methods assume strictly positive eigenvalues and so require
some adjustment to address these negative eigenvalues, such as setting them to
zero [25].

We propose a preprocessing step to improve the robustness of pairwise LD
matrices to missing data. Instead of using the pairwise LD matrix, we can
leverage the known structure of linkage disequilibrium given by Equation 13
and Theorem 1 to compute the optimal theoretical matrix bytaking appropriate
means of the pairwise LD. This not only replaces the observed matrix with one
which follows the prescribed theoretical structure, but as the operation is
equivalent to taking structured means, it should reduce the impact of
individual missing entries.

###  III-A Simulated data

To test this proposal, a simple simulation study was performed. 100 synthetic
populations of 100 individuals were generated measured at 20 markers with
$d(i,j)=15$ cM and $c_{i}=c_{j}=1$ for each pair, a setting corresponding to a
Toeplitz LD matrix. For each population, the complete data was used to compute
the observed LD matrix and determine its eigenvalues. An increasing proportion
of observations were then removed completely at random from the data to
simulate different levels of data completeness, and at each proportion the
pairwise LD matrix was computed and used to generate the nearest Toeplitz
matrix. The minimum eigenvalue and sum of squared errors in the ordered
eigenvalues from the complete data LD matrix were computed and recorded for
both the pairwise LD matrix and the nearest Toeplitz matrix. Figure 1 shows
the results.

---
(a)
(b)
Figure 1: Paired boxplots of the (a) minimum eigenvalues and (b) sum of
squared errors in the ordered eigenvalues for the pairwise LD matrix and the
nearest Toeplitz matrix by the proportion of data missing. The nearest
Toeplitz, displayed to the right of the line for each pair of boxplots, is
more robust to missing data than the pairwise LD matrix, displayed to the left
of the line for each pair, but is biased when the data are complete.

Figure 1(a) shows that the nearest Toeplitz is more robust to negative
eigenvalues than the pairwise LD matrix alone. A vast majority of the 100
simulated populations have nearest Toeplitz matrices with no negative
eigenvalues until more than a third of the data is missing, and even then only
about 25% produce negative eigenvalues. In contrast, the pairwise LD matrix
has negative eigenvalues more than 75% of the time when as little as 15% of
the data is missing. Figure 1(b) indicates this robustness also extends to the
sum of squared differences between ordered eigenvalues, which does not depend
greatly on the completeness of the data for the nearest Toeplitz but does for
the pairwise LD. The cost of this robustness and the reduction in space
complexity from $O(n^{2})$ to $O(2n)$ is potential bias when the data is
(mostly) complete.

###  III-B Real data

Following the simulated study of the previous section, we now corroborate our
findings by replicating the experiment with real data. The Mouse Genome
Database (MGD) is a public repository of information on mouse genomics,
including references which describe hundreds of thousands of markers and their
measured values in several experiments on live mice [1]. One of these data
sets records markers measured on the backcross of a mutated strain of _Mus
musculus_ and a strain of _Mus spretus_ in [26]. Partial measurement of 5951
markers across all chromosomes on 94 animals is reported, with complete
observations present on 2624 of the markers used to create a 2624$\times$2624
observed correlation matrix. In contrast to the previous simulated example,
these markers are not positioned uniformly across the chromosomes; the cM
distances between adjacent markers differ greatly.

Instead of focusing on all of the markers, however, we make the real data
comparable to the simulated case by considering only those markers on
chromosome 1. This leaves 199 markers measured across the chromosome with
adjacent distances ranging from less than 0.01 cM to 4.37 cM. The cM distances
were used to generate a theoretical correlation matrix according to Equation
13, these theoretical correlations were rounded to two decimal places, and
these rounded correlations were treated as the level sets of the theoretical
structured matrix. Note that unlike the examples illustrated above, these
level sets do not correspond with a named structured matrix. Nonetheless,
Theorem 1 dictates that the optimal structured approximation in the Frobenius
norm is given by means computed over the level sets. An unfortunate
consequence of the lack of correspondence with a named matrix is that there is
no guarantee that the result will be a valid correlation matrix.

After these steps, we have two correlation matrices: the observed matrix on
the full data and the theoretical correlation matrix based on cM distances. To
examine the impact of missing data on the computed eigenvalues, marker
measurements were removed completely at random 100 times for each of a range
of proportions of missingness and the pairwise LD matrix and the optimal
structured approximation based on the theoretical level sets were generated.
For both the pairwise LD matrix and the optimal structured matrix, the minimum
eigenvalue and sum of squared errors from the matrix of complete observations
was computed. Figure 2 displays the results.

---
(a)
(b)
Figure 2: Paired boxplots of the (a) minimum eigenvalues and (b) sum of
squared errors in the ordered eigenvalues for the pairwise LD matrix (to the
left of the corresponding line) and the optimal structured matrix (to the
right of the corresponding line) by the proportion of data missing. The bias
from approximation is more serious in this case than the simulated example:
negative eigenvalues are produced for the complete data.

First, note that the larger correlation matrix for the real data example has
resulted in less relative variability for our summaries of eigenvalues in both
the pairwise LD and optimal structured matrices, making the difference in
perfomance even starker. Just as in the simulated case, the optimal structured
matrix proves far more robust to missing data. Both the minimum eigenvalue and
the sum of squared errors barely change on average as the proportion of
missing data increases. In contrast, the pairwise LD matrix has an error that
grows in the proportion of missing data and a minimum eigenvalue which
continues to decrease as data is removed. When as little as 15% of the data is
missing, the optimal structured matrix has a better distribution in both
metrics.

The bias observed in the simulated case has more serious consequences in the
real data, however. For mostly complete data, the optimal structured matrix
produces negative eigenvalues, and so the optimal structured matrix clearly
does not correspond with a correlation matrix. This limits its use for cases
with nearly complete data, where an optimal structured approximation should
only be used if the structure is Toeplitz, circulant, or of another class with
known qualities to its eigenspectrum. Future work will examine a more
sophisticated approximation in the form of a structured matrix plus a diagonal
matrix. This form retains a reduced space complexity, but the diagonal
components provide variables that can potentially be optimized to ensure
desired properties, e.g., postive semidefiniteness for correlation matrices.

##  IV Conclusion

In this paper, we have presented a general method for obtaining structured
approximations to arbitrary matrices. Specifically, we have proved that the
Frobenius norm of $\mathbf{M}$ to a matrix $\mathbf{T}$ with structural
equality of certain entries is minimized by replacing the entries in
$\mathbf{M}$ by their means for each index set of $\mathbf{T}$. The value of
this minimum corresponds to the weighted standard deviation of entries within
each index set, giving an intuitive way to measure the expected error for this
approximation. In the circulant case, this is equivalent to taking the
circulant matrix with eigenvalues
$\operatorname{diag}\left(\mathbf{F}\mathbf{M}{\mathbf{F}}^{\ast}\right)$,
where $\mathbf{F}$ is the DFT matrix.

We believe that our general structured approximation method is of significant
independent theoretical and practical interest. Theoretically, it promises a
simple and computationally-efficint algorithm to compute the Frobenius-optimal
matrix with structured equalities to an arbitrary matrix. Evidence for its
practical value comes from our example involving the creation of a structured
approximation of the pairwise genetic linkage disequilibrium matrix that is
highly robust to missing data. In both simulated and real data examples, the
eigenvalues changed little for the optimal structured matrix when as much as
40% of the data was missing, though at the cost of some bias when the data was
nearly complete. Despite this remarkable result, it must be emphasized that
the new structured approximation method is potentially applicable to any
problem involving large-scale matrices for which standard linear algebraic
methods cannot practically be applied.

##  V References

## References

  * [1] C. J. Bult, J. A. Blake, C. L. Smith, J. A. Kadin, J. E. Richardson, and t. M. G. Database Group, “Mouse genome database (MGD) 2019,” _Nucleic Acids Research_ , vol. 47, no. D1, pp. D801–D806, 2019.
  * [2] Z. Bai and J. W. Silverstein, _Spectral analysis of large dimensional random matrices_. Springer, 2010, vol. 20.
  * [3] K. Chandramouli and P. Qian, “Proteomics: challenges, techniques and possibilities to overcome biological sample complexity,” _Hum Genomics Proteomics_ , 2009.
  * [4] R. Poulos, P. Hains, R. Shah, and et al., “Strategies to enable large-scale proteomics for reproducible research,” _Nature Communications_ , vol. 11, no. 3793, 2020.
  * [5] Y. Zhou, “Set-based differential covariance testing for genomics,” _Stat (Int Stat Inst)_ , vol. 8, no. 1, pp. 251–280, 2019.
  * [6] P. Danaher and P. Wang, “Covariance-based analyses of biological pathways,” _Biometrika_ , vol. 102, no. 3, pp. 533–544, 2015.
  * [7] J. Baik and J. W. Silverstein, “Eigenvalues of large sample covariance matrices of spiked population models,” _Journal of multivariate analysis_ , vol. 97, no. 6, pp. 1382–1408, 2006.
  * [8] J. Li and S. X. Chen, “Two sample tests for high-dimensional covariance matrices,” 2012.
  * [9] J. Liu, Y. Liang, and N. Ansari, “Spark-based large-scale matrix inversion for big data processing,” _IEEE Access_ , vol. 4, pp. 2166–2176, 2016.
  * [10] V. Strassen, “Gaussian elimination is not optimal,” _Num. Math_ , vol. 13, p. 354–356, 1969.
  * [11] D. Coppersmith and S. Winograd, “Matrix multiplication via arithmetic progressions,” _J. Symb. Comput._ , vol. 9, no. 3, pp. 251–280, 1990.
  * [12] J. Alman and V. V. Williams, “A refined laser method and faster matrix multiplication,” in _Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms, SODA 2021, Virtual Conference, January 10 - 13, 2021_ , D. Marx, Ed. SIAM, 2021, pp. 522–539.
  * [13] S. Robinson, “Toward an optimal algorithm for matrix multiplication,” _SIAM News_ , vol. 38, no. 9, 2005.
  * [14] A. U. Jörn Lötsch, Sebastian Malkusch, “Optimal distribution-preserving downsampling of large biomedical data sets (opdisdownsampling),” _PLOS One_ , 2021.
  * [15] C. Jain, “Coverage-preserving sparsification of overlap graphs for long-read assembly,” _Bioinformatics_.
  * [16] R. M. Gray, _Toeplitz and circulant matrices: A review_. now Publishers Inc., 2006.
  * [17] T. F. Chan, “An optimal circulant preconditioner for Toeplitz systems,” _SIAM Journal on Scientific and Statistical Computing_ , vol. 9, no. 4, pp. 766–771, 1988.
  * [18] M. Hariprasd and M. Venkatapathi, “Circulant decomposition of a matrix and the eigenvalues of toeplitz type matrices,” _arXiv preprint arXiv:2105.14805_ , 2022.
  * [19] O. Narayan and B. S. Shastry, “Generalized Toeplitz–Hankel matrices and their application to a layered electron gas,” _Journal of Physics A: Mathematical and Theoretical_ , vol. 54, no. 17, p. 175201, 2021.
  * [20] E. S. Lander and D. Botstein, “Mapping mendelian factors underlying quantitative traits using RFLP linkage maps,” _Genetics_ , vol. 121, no. 1, pp. 185–199, 1989.
  * [21] E. Uffelmann, Q. Q. Huang, N. S. Munung, J. de Vries, Y. Okada, A. R. Martin, H. C. Martin, T. Lappalainen, and D. Posthuma, “Genome-wide association studies,” _Nature Reviews Methods Primers_ , vol. 1, no. 1, pp. 1–21, 2021\.
  * [22] C. Salahub, “A structural model of genome-wide association studies,” _arXiv preprint arXiv:2205.10391_ , 2022.
  * [23] K. N. Conneely and M. Boehnke, “So many correlated tests, so little time! Rapid adjustment of P values for multiple correlated tests,” _The American Journal of Human Genetics_ , vol. 81, no. 6, pp. 1158–1168, 2007.
  * [24] B. Han, H. M. Kang, and E. Eskin, “Rapid and accurate multiple testing correction and power estimation for millions of correlated markers,” _PLoS genetics_ , vol. 5, no. 4, p. e1000456, 2009.
  * [25] J. Li and L. Ji, “Adjusting multiple testing in multilocus analyses using the eigenvalues of a correlation matrix,” _Heredity_ , vol. 95, no. 3, pp. 221–227, 2005.
  * [26] L. Rowe, J. Nadeau, R. Turner, W. Frankel, V. Letts, J. Eppig, M. Ko, S. Thurston, and E. Birkenmeier, “Maps from two interspecific backcross DNA panels available as a community genetic mapping resource,” _Mammalian Genome_ , vol. 5, no. 5, pp. 253–274, 1994.