Rotated Eigenstructure Analysis for Source Localization without Energy-decay Models
Junting Chen Urbashi Mitra




Rotated Eigenstructure Analysis for 
Source Localization without
Energy-decay Models
    C. Garza
    
====================================================================================
















Ginwidth=0.5
	 
Ginwidth=1.0
	 

Herein, the problem of simultaneous localization of two sources given
a modest number of samples is examined. In particular, the strategy
does not require knowledge of the target signatures of the sources
a priori, nor does it exploit classical methods based on a
particular decay rate of the energy emitted from the sources as a
function of range. General structural properties of the signatures
such as unimodality are exploited. The algorithm localizes targets
based on the rotated eigenstructure of a reconstructed observation
matrix. In particular, the optimal rotation can be found by maximizing
the ratio of the dominant singular value of the observation matrix
over the nuclear norm of the optimally rotated observation matrix.
It is shown that this ratio has a unique local maximum leading to
computationally efficient search algorithms. Moreover, analytical
results are developed to show that the squared localization error
decreases at a rate n^-3 for a Gaussian field with a single source,
where n(log n)^2 scales proportionally to the number of samples
M.




§ INTRODUCTION




Underwater source detection and localization is an important but challenging
problem. Classical range-based or energy-based source localization
algorithms usually require energy-decay models and the knowledge of
the environment <cit.>.
However, critical environment parameters may not be available in many
underwater applications, in which case, classical model-dependent
methods may break down, even when the measurement snr is high. 

There have been some studies on source localization using nonparametric
machine learning techniques, such as kernel regressions and support
vector machines <cit.>.
However, these methods either require a large amount of sensor data,
or some implicit information of the environment, such as the choice
of kernel functions. For example, determining the best kernel parameters
(such as bandwidth) is very difficult given a small amount of data. 

This paper focuses on source detection and localization problems when
only some structural properties of the energy field generated by the
sources are available. Specifically, instead of requiring the knowledge
of how energy decays with distance to the source, the paper aims at
exploiting only the assumption that the closer to the source the higher
energy received, and moreover, the energy field of the source is spatially
invariant and decomposable. In fact, such a structural property is
generic in many underwater applications. The prior work <cit.>
studied the single source case, where an observation matrix is formed
from a few energy measurements of the field in the target area, and
the missing entries of the observation matrix are filled using matrix
completion methods. Knowing that the matrix would be rank-1 under
full and noise-free sampling of the whole area, svd is applied
to extract the dominant singular vectors, and the source location
is inferred from analyzing the peaks of the singular vectors. 

Herein, we propose to improve upon two shortcomings in <cit.>:
we make rigorous an estimation/localization bound (versus focusing
on the reduction of the search region) and we provide a method for
localizing two sources. In the two source case, we need to tackle
an additional difficulty that the svd of the observation matrix
does not correspond to the signature vectors of the sources. To resolve
this issue, a method of rotated eigenstructure analysis is proposed,
where the observation matrix is formed by rotating the coordinate
system such that the sources are aligned in a row or in a column of
the matrix. We develop algorithms to first localize the central axis
of the two sources, and then separate the sources on the central axis. 

To summarize, we derive algorithms to simultaneously localize up to
two sources based on only a few power measurements in the target area
without knowing any specific energy-decay model. The contributions
of this paper are as follows:


  * We derive the location estimators with analytical results to show
that the squared error decreases at a rate n^-3 for a Gaussian
field with a single source, where n(log n)^2 scales proportionally
to the number of samples M. 

  * We develop a localization algorithm for the double source case based
on a novel rotated eigenstructure analysis. We show that the two sources
can be separated even when their aggregate power field has a single
peak. 



The rest of the paper is organized as follows. Section <ref>
gives the system model and assumptions. Section <ref>
develops location estimator with performance analysis for single source
case. Section <ref> proposes rotated eigenstructure
analysis for double source case. Numerical results are given in Section
<ref> and Section <ref> concludes this
work. 



§ SYSTEM MODEL






Consider that there are K (K=1,2) sources with unknown locations
𝐬_k=(x_k^S,y_k^S)∈ℝ^2
located in a bounded area 𝒜. Suppose that the sensors
can only measure the aggregate power transmitted by the sources, and
is given by 

    h(x,y)=∑_kh_k(x,y)

for measurement location (x,y), where 

    h_k(x,y)=α u(x-x_k^S)u(y-y_k^S)

is the power density from source k, where α>0. The explicit
form of the density function h_k(x,y) is unknown to the system,
except that the characteristic function u(x) is known to
have the following properties


  a)  positive semi-definite, i.e., u(x)≥0 for all x∈ℝ

  b)  symmetric, i.e., u(x)=u(-x)

  c)  unimodal, i.e., u^'(x)<0 for x>0, 

  d)  smooth, i.e., |u^'(x)|<K_u for some K_u>0, and 

  e)  normalized, i.e., ∫_-∞^∞u(x)^2dx=1. 

Note that u(x) can be considered as the marginal power density
function. 

Consider that M power measurements {h^(l)} are taken over
distinct locations 𝐳^(l)=(x^(l),y^(l)), l=1,2,…,M,
uniformly at random in the target area 𝒜. The measurements
are assigned to a n_1× n_2 observation matrix 𝐇̂
as follows. First, partition the target area 𝒜 into n_1× n_2
disjoint cells 𝒢_ij, i=1,2,…,n_1 and j=1,2,…,n_2,
where n_1 and n_2 are to be determined. Second, assign the
power measurements h^(l) to the corresponding (i,j)th entry
of 𝐇̂ as 

    Ĥ_ij=s(𝒢_ij)h^(l)

if 𝐳^(l)∈𝒢_ij,where s(𝒢_ij)
measures the area of 𝒢_ij.[If multiple samples are close to each other and assigned to the same
entry of 𝐇̂, the value of that entry is the average
of the sample values.] Denote Ω as the set of observed entries of 𝐇̂,
i.e., (i,j)∈Ω if there exists 𝐳^(l)∈𝒢_ij
such that h^(l) is assigned to Ĥ_ij. 

For easy discussion, assume that 𝒜=[-L/2,L/2]×[-L/2,L/2],
n_1=n_2=n, and 𝒢_ij are rectangles centered
at (x_i,y_j), x_i=-L/2+L/2n+L/n(i-1),
y_j=-L/2+L/2n+L/n(j-1), and have identical
size with each other. Let 𝐇=α∑_k=1^K𝐮_k𝐯_k^T
be the matrix of ideal observation, where 

    𝐮_k    =L/N[u(x_1-x_k^S),u(x_2-x_k^S),…,u(x_n-x_k^S)]^T
    𝐯_k    =L/N[u(y_1-y_k^S),u(y_2-y_k^S),…,u(y_n-y_k^S)]^T

for k=1,2. Thus 𝐇 has rank at most K. For (i,j)∈Ω,
we have Ĥ_ij≈ H_ij, where the slight difference
is due to sampling away from the centers of the cells 𝒢_ij.
As a result, 𝐇̂ is a sparse and noisy observation
of the low rank matrix 𝐇. An application example is illustrated
in <ref>.

The goal of this paper is to find the approximate locations of the
sources using only the spatial invariant property (<ref>)
and the four generic properties of the characteristic function u(x).
Note that this problem is non-trivial. We insist on several features
of the algorithm to be developed: it should be robust to structural
knowledge of the signatures of the sources (as captured by g(x,y)
in (<ref>)). This disallows the use
of parametric regression or parameter estimation for source localization.
In addition, we wish to under-sample the target area using small M.
As such, maximum value entries may not represent the true locations
of the sources. While not a focus of the current work, we will use
matrix completion methods and the low rank property of 𝐇
as in <cit.> to cope with the under-sampled
observations. 







§ EIGENSTRUCTURE ANALYSIS FOR 
SINGLE SOURCE LOCALIZATION




To simplify the discussion, the following mild assumptions are made.[The two assumptions are mainly to avoid discussing the effects on
the boundary of 𝒜 and the high order noise term in the
sampling noise model (<ref>). Straight-forward
modifications can be made to handle the boundary effect in practical
algorithms.] 


  A1)  The observation area 𝒜 is large enough, such that there
is only negligible energy spreading outside the area 𝒜.

  A2)  The parameter n is not too small, such that u(x_i-x_k^S)^2δ^2≈∫_x_i^x_i+1u(x-x_k^S)^2dx
and u(y_i-y_k^S)^2δ^2≈∫_y_i^y_i+1u(y-y_k^S)^2dy
for all i=1,2,…,n. 

Mathematically, the above assumptions imply that the vectors 𝐮_k
and 𝐯_k have unit norm.



 §.§ Observation Matrix Construction


We first exploit the low rank property of 𝐇 to obtain
the full matrix 𝐇̂_c from the partially
observed matrix 𝐇̂. Let 𝒫_Ω(𝐗)
be a projection, such that the (i,j)th element of matrix 𝒫_Ω(𝐗)
is [𝒫_Ω(𝐗)]_ij=X_ij if (i,j)∈Ω,
and [𝒫_Ω(𝐗)]_ij=0 otherwise.
The completed matrix 𝐇̂_c can be found as
the unique solution to the following problem 

    𝐗minimize     𝐗_*
    subject to     𝒫_Ω(𝐗-𝐇̂)_F≤ϵ

where 𝐗_* denotes the nuclear norm of 𝐗
and ϵ is a small parameter to tolerate the discrepancy between
the two matrices. 

To choose a proper dimension n for the observation matrix 𝐇̂_c∈ℝ^n× n,
we consider the results in <cit.>. It has been shown that
under some mild conditions of 𝐇 (such as the strong incoherence
property and small rank property), the matrix 𝐇∈ℝ^n× n
can be exactly recovered with a high probability, if the dimension
n satisfies Cn(log n)^2≤ M and noise-free sampling, Ĥ_ij=H_ij
for (i,j)∈Ω, is performed. Here, C is a positive constant.
Given this, we propose to choose n=n_c as the largest
integer to satisfy n_c(log n_c)^2≤ M/C. 




 §.§ Location Estimator Exploiting Property of Symmetry


Consider the svd of the completed matrix 𝐇̂_c
as 𝐇̂_c=α_1𝐮̂_1𝐯̂_1^T+∑_i=2^n_cα_i𝐮̂_i𝐯̂_i^T.
We thus model the singular vectors of 𝐇̂_c as 𝐮̂_1=𝐮_1+𝐞_u
and 𝐯̂_1=𝐯_1+𝐞_v.

Note that the vectors 𝐮_1 and 𝐯_1 defined
in (<ref>) and (<ref>), respectively, contain the source
location information due to the unimodal property of u(x). However,
due to the noise vectors 𝐞_u and 𝐞_v,
the source location cannot be found by simply locating the peaks of
𝐮̂_1 and 𝐯̂_1. 

To resolve this difficulty, we exploit the symmetric property of u(x)
and develop a location estimator as follows.

Define a reflected correlation function as 

    R̂(t;𝐮̂_1)=∫_-∞^∞û(x)û(-x+t)dx

where û(x) is a (nonparametric) regression function from
vector 𝐮̂_1. For example, û(x) can be
obtained by û(x)=𝐮̂_1(i) if x=x_i, and
by linear interpolation between 𝐮̂_1(i) and 𝐮̂_1(i+1)
if x_i<x<x_i+1. Then the location estimator for x_1^S
is given by 

    x̂_1^S(𝐮̂_1)=1/2t∈ℝargmax R̂(t;𝐮̂_1).

The location estimator for y_1^S can be obtained in
a similar way.

The location estimator (<ref>) exploits the
fact that as 𝐮̂_1 is symmetric, the reflected correlation
(<ref>), which is the correlation between
𝐮̂_1 and a reflected and shifted version of 𝐮̂_1,
is maximized at the source location. Therefore, the estimator x̂_1^S(𝐮̂_1)
tries to the suppress the perturbation from the noise by correlating
over all the entries of 𝐮̂_1. 

We establish several properties for the estimator x̂_1^S(𝐮̂_1). 

Consider the autocorrelation for the characteristic function u(x)
as 

    τ(t)=∫_-∞^∞u(x)u(x-t)dx.

Then, the following property can be derived. 
  The autocorrelation function
τ(t) is non-negative and symmetric. In addition, τ(t)
is strictly decreasing in t>0.



Let the dominant singular vector of 𝐇̂_c
as the solution to (<ref>) be given by 𝐮̂_1=𝐮_1+𝐞_1,
where 𝐮_1 is the dominant singular of 𝐇.
Let 𝐞_1 be a vector with reverse elements
of 𝐞_1, i.e., the jth element of 𝐞_1
equals to the last but the jth element of 𝐞_1. Let
𝐞_1^-t be a vector obtained from the t-shift of
𝐞_1, i.e., for t>0, the first t
elements of 𝐞_1^-t are zeros and the remaining (n_c-t)
elements of 𝐞_1^-t are identical to the first (n_c-t)
elements of 𝐞_1; and for t<0, the
first (n_c-t) elements of 𝐞_1^-t are identical
to the last (n_c-t) elements of 𝐞_1
and the remaining t elements of 𝐞_1^-t are zeros.
With such a notion, we make the following assumption on the singular
vector 𝐮̂_1=𝐮_1+𝐞_1 of the
completed matrix 𝐇̂_c: 

    |𝐮_1^T𝐞_1^-t|≤ C_e|𝐮_1^T𝐞_1|

for any 0≤ t≤ n_c-1, where C_e<∞ is a positive
constant that only depends on the characteristic function u(x)
but not n_c or M. 

Such an approximation is motivated by two observations. First, the
entries of the vector 𝐞_1 may have roughly the same
chance to take positive values or negative values because both 𝐮_1
and 𝐮_1+𝐞_1 have unit norm. Second, the magnitude
of the elements in 𝐮_1 depends only on the characteristic
function u(x) but not n_c or M. Although it is difficult
to analytically validate the assumption (<ref>),
it can be roughly confirmed by massive simulation results. 

As a result, we have the following theorem to characterize the estimation
error of x̂_1^S.

Suppose that the sampling error of 𝐇̂ from the true
energy field matrix 𝐇 is bounded by 𝒫_Ω(𝐇̂-𝐇)_F≤ϵ̅
and the algorithm parameter ϵ in (<ref>)
is chosen as ϵ=ϵ̅. Then, with high probability,

    |x̂_1^S-x_1^S|≤1/2τ^-1(1-μ_uL^6n_c(M)^-3+o(n_c(M)^-3))

where τ^-1(r) is the inverse function of r=τ(t), μ_u=C_e128u(0)^2K_u^2,
and n_c(M) is the largest integer chosen such that M≥ Cn_c(log n_c)^2.



The specific performance from (<ref>) depends
on the characteristics of the energy field. Intuitively, if u(x)
has a sharp peak (large slope of the autocorrelation function τ(t)),
the localization error should be smaller. Consider a numerical example
where the energy field has a Gaussian characteristic function.

For a Gaussian characteristic function u(x)=(2γ/π)^1/4e^-γ x^2,
there exists a constant C_μ, which only depends on the characteristic
function u(x), such that with high probability, the squared estimation
error is upper bounded by 

    |x̂_1^S-x_1^S|^2+|ŷ_1^S-y_1^S|^2≤ C_μL^6n_c(M)^-3+o(n_c(M)^-3).




Theorem <ref> and Corollary <ref>
gives the asymptotic performance of the proposed localization algorithm
without knowing the energy-decay model. For large M, the worst
case squared error decays at a rate n_c(M)^-3. As a benchmark,
the squared error of a naive scheme, which estimates the source location
directly from the position of the measurement sample that observes
the highest power, decreases as M^-1, which is equivalent to
n_c(M)^-1(log n_c(M))^-2, much slower than that of the
proposed algorithm. This is because, the granularity of the original
observations is L/√(M). The results then confirm that by exploiting
the low rank property using matrix completion and the reflected correlation
technique, the proposed algorithm significantly improves the localization
resolution.



§ ROTATED EIGENSTRUCTURE ANALYSIS FOR 
DOUBLE SOURCE LOCALIZATION




The location estimator x̂_1^S in (<ref>)
is based on the intuition that the singular vectors of 𝐇
are just the vectors 𝐮_1 and 𝐯_1, which
contains the source location in their peaks. However, a similar technique
cannot be applied to the two source case, because 𝐮_k
and 𝐯_k may not be the singular vectors of 𝐇,
as the vectors {𝐮_k} may not be orthogonal.



 §.§ Optimal Rotation of the Observation Matrix


When there are two sources, the (ideal) observation matrix 𝐇
is not rank-1, expect for the special case where the two sources are
aligned on one of the axes of the coordinate system. 

wlog, assume that the sources are aligned with the x-axis,
where y_k^S=C for k=1,2. Consequently, we have 𝐯_1=𝐯_2,
and 𝐇=α(∑_k𝐮_k)𝐯_1^T,
which is rank-1. Hence, the right singular vector of 𝐇
is 𝐯_1 and, by analyzing the peak of 𝐯_1,
the central axis ŷ_k^S=C can be estimated. 

The above observations suggest that we rotate the coordinate system
such that the sources are aligned with one of the axes. Consider rotating
the coordinate system by θ. The entries of 𝐇̂_c
are rearranged into a new observation matrix 𝐇̂_θ,
where 

    [𝐇̂_θ]_(i,j)=[𝐇̂_c]_(p,q)

in which (p,q) is the index such that (x_p^',y_q^')
is the closest point in Euclidean distance to (x̅,y̅)
in the original coordinate system 𝒞_0, with x̅=dcos(β+θ)
and y̅=dsin(β+θ). Here β=∠(x_i,y_j)
is the angle of (x_i,y_j) to the x-axis of the rotated coordinate
system 𝒞_θ, and d=(x_i,y_j)_2. Note
that 1≤ i,j≤ n^', where n^'≤ n_c, since the rotation
of the axes induce truncation of some data samples. 

Let the orientation angle of the central axis of the sources wrt
the x-axis in the original coordinate system 𝒞_0
be θ_0, θ_0∈[0,π). Then the desired rotation
for coordinate system 𝒞_θ would be θ^*=θ_0
for θ_0<π/2, or θ^*=θ_0-π/2
for θ_0≥π/2. The desired rotation θ
can be obtained as 

    θ∈[0,π/2]maximize  ρ(θ)≜λ_1(𝐇̂_θ)/∑_kλ_k(𝐇̂_θ)

where λ_k(𝐀) is the kth largest singular value
of 𝐀. Note that ρ(θ)≤1 for all θ∈[0,π/2]
and ρ(θ^*)=1, where 𝐇̂_θ becomes
a rank-1 matrix when the sources are aligned with one of the axes.

The maximization problem (<ref>) is in general non-convex.
An exhaustive search for the solution θ^* is computationally
expensive, since for each θ, svd should be performed
to obtain the singular value profile of 𝐇̂_θ.
Therefore, we need to study the properties of the alignment metric
ρ(θ) in order to develop efficient algorithms for the source
detection.




 §.§ The Unimodal Property




We also show that the function ρ(θ) also has the unimodal
property defined as follows. 
 A function f(x) is called unimodal in a bounded
region (a,b), if there exists x_0∈[a,b], such that f^'(x)f^'(y)<0
for any a<x<x_0<y<b. 



The unimodality suggests that f(x) has a single peak in (a,b),
and hence f(x) has a unique local maximum (or minimum). 

The function ρ(θ) in (<ref>) is unimodal
in θ∈(θ^*-π/4,θ^*+π/4),
if 

    s·τ^'(t)>t·τ^'(s)

for all 0<s<t, where τ^'(t)≜d/dtτ(t).
In addition, ρ(θ) is strictly increasing over (θ^*-π/4,θ^*)
and strictly decreasing over (θ^*,θ^*+π/4). 



The result in Theorem <ref> is powerful,
since it confirms that the function ρ(θ) is unimodal within
a π/2-window, and there is a unique local maximum, when
the autocorrelation of the energy field characteristic function u(x)
agrees with the condition (<ref>). Note
that ρ(θ) is also symmetric wrt θ=θ^*.
As a result, a simple bisection search algorithm can efficiently find
the global optimal solution θ^* to (<ref>).
An example algorithm is given in Algorithm <ref>.



Note that condition (<ref>) can be satisfied
by a variety of energy fields. For example, for Laplacian field u(x)=√(γ)e^-γ|x|,
we have τ(t)=(1+γ t)e^-γ t, and τ^'(t)=-γ^2te^-γ t;
for Gaussian field u(x)=(2γ/π)^1/4e^-γ x^2,
we have τ(t)=e^-γ t^2/2, and τ^'(t)=-γ te^-γ t^2/2.
In both cases, condition (<ref>) is satisfied. 









 §.§ Source Detection






In the coordinate system 𝒞_θ under optimal rotation
θ=θ^* (assuming alignment on the x-axis), the left
and right singular vectors of 𝐇̂_θ can be
modeled as 𝐮̂_1=1/2(𝐮_1(θ^*)+𝐮_2(θ^*))+𝐞_u
and 𝐯̂_1=𝐯_1(θ^*)+𝐞_v,
respectively. Correspondingly, the y-coordinates of the sources
can be the found using estimator (<ref>) based
on reflected correlation 

    ŷ_1^S(𝐯̂_1;θ^*)=ŷ_2^S(𝐯̂_1;θ^*)=1/2t∈ℝargmax R̂(t;𝐯̂_1).


To find the x-coordinates, note that the function u_1(x)=1/2(u(x-x_1^S)+u(x-x_2^S))
is symmetric at x=1/2(x_1^S+x_2^S).
Therefore, the center of the two sources can be found by 

    ĉ=1/2t∈ℝargmax R̂(t;𝐮̂_1).


In addition, after estimating ŷ_1^S, the marginal
power density function u(x) can be obtained as û(y)=v̂_1(y-ŷ_1^S),
where v̂_1(y) is a regression function from 𝐯̂_1
(for example, by linear interpolation among y_1,y_2,…,y_n_c).
As a results, the x-coordinates of the two sources can be found
using similar techniques as spread spectrum early gate synchronization
<cit.>, and obtained as x̂_1^S(θ^*)=ĉ-d̂
and x̂_2^S(θ^*)=ĉ+d̂, where 

    d̂    =d≥0argmax   Q(d;𝐮̂_1,𝐯̂_1)

and 

    Q(d;𝐮̂_1,𝐯̂_1)≜1/2∫_-∞^∞û_1(x)(û(x-ĉ-d)+û(x-ĉ+d))dx.

It is straight-forward to show that Q(d;𝐮̂_1,𝐯̂_1)
is maximized at d^*=1/2|x_1^S-x_2^S|.

As a benchmark, consider a naive scheme that estimates x_1^S
and x_2^S by analyzing the peaks of 𝐮̂_1.
However, such naive strategy cannot work for small source separation,
because if d=1/2|x_1^S-x_2^S| is
too small, the aggregate power density function ũ_1(x)=u(x-x_1^S)+u(x-x_1^S-d)
would be unimodal and there is only one peak in 𝐮̂_1.
As a comparison, the proposed procedure estimator from procedure (<ref>)–(<ref>)
does not such a limitation. 







§ NUMERICAL RESULTS






In this section, we evaluate the performance of the proposed location
estimator in both single source and double source cases. Two sources
are placed in the area [-0.5,0.5]×[-0.5,0.5] uniformly and
independently at random, with the restriction that the distance between
the two sources is no more than 0.5.[When the two sources are far apart, the problem degenerates to two
single-source-localization problems.] The power field generated by each source in an underwater environment
is modeled as h_k(x,y)=e^-20(x-x_k^S)^2-20(x-y_k^S)^2,
k=1,2. There are M power measurements taken in the area 𝒜=[-1,1]×[-1,1]
uniformly at random. The parameter n_c of the proposed observation
matrix 𝐇̂∈ℝ^n_c× n_c is chosen
as the largest integer satisfying n_c(log n_c)^2≤ M/C,
for C=1. 

As a benchmark, the proposed location estimation is compared with
the naive scheme, which determines the source location directly form
the position of the measurement sample that observes the highest power.
In the two source case, the naive algorithm aims at detecting either
one of the sources, and the corresponding localization error is computed
as ℰ_ naive^2=min{𝐬̂_ naive-𝐬_1^2,𝐬̂_ naive-𝐬_2^2}.
As a comparison, the localization error of the proposed algorithm
is computed as ℰ^2=1/2(𝐬̂_1-𝐬_1^2+𝐬̂_2-𝐬_2^2).



Fig. <ref> depicts the mse of the source location
versus the number of samples M. In the single source case,the coefficient
of the worst case upper bound (<ref>)
is chosen as C_μ=1 to demonstrate the asymptotic decay rate
of the worst case squared error bound. The decay rate of the analytic
worst case error bound is roughly the same as the mse obtained
from the numerical experiment. It is expected that as M increases,
the two curves merge in an asymptotic way. As a benchmark, the proposed
scheme requires less than half of the samples to achieve similar performance
to that of the naive baseline even for small M (around 50).
More importantly, it demonstrates a higher mse decay rate, where
for medium M (around 200), the proposed scheme reduces the number
of samples to 1/10. In the double source case, there is an error
floor for the naive scheme, because the location that observes the
highest power may not be either one of the source locations. As a
comparison, there is no error floor in for proposed scheme as M
increases. 

Fig. <ref> shows an example on simultaneously localizing two
sources (red crosses). Although the aggregate power field has only
one peak, the algorithm (black circles) is able to separate the two
sources. 



§ CONCLUSIONS




This paper developed source localization algorithms from a few power
measurement samples, while no specific energy-decay model is assumed.
Instead, the proposed method only exploited the structural property
of the power field generated by the sources. Analytical results were
developed to demonstrate that the proposed algorithm decreases the
localization error at a higher rate than the baseline algorithm when
the number of samples increases. In addition, a rotated eigenstructure
analysis technique was derived for simultaneously localizing two sources.
Numerical results demonstrate the performance advantage in localizing
single or double sources.



§ ACKNOWLEDGMENTS


This research was supported, in part, by National Science Foundation
under Grant NSF CNS-1213128, CCF-1410009, CPS-1446901, Grant ONR N00014-15-1-2550,
and Grant AFOSR FA9550-12-1-0215.







§ APPENDIX




 §.§ Proof of Lemma <ref>



    τ^'(t)     =d/dt∫_-∞^∞u(x)u(x-t)dx
        =∫_-∞^∞-u(x)u^'(x-t)dx
        =-∫_-∞^0u(z+t)u^'(z)dz-∫_0^∞u(z+t)u^'(z)dz
        =-∫_-∞^0u(z+t)u^'(z)dz+∫_0^∞u(z+t)u^'(-z)dz
        =-∫_-∞^0u(z+t)u^'(z)dz+∫_-∞^0u(-w+t)u^'(w)dw
        =-∫_-∞^0[u(z+t)-u(-z+t)]u^'(z)dz
        =-∫_-∞^0[u(z+t)-u(z-t)]u^'(z)dz
        <0

where (<ref>) is due to the change of variable
z=x-t and u^'(z)=-u^'(-z), (<ref>) is
to change the variable z=-w, (<ref>) exploits
the fact that u(x)=u(-x), and the last inequality is due to u(z+t)-u(z-t)>0
and u^'(z)>0 for all z<0. 





 §.§ Proof of Theorem <ref>


To simplify the algebra, we only focus on the dominant terms wrt
n_c as n_c goes large. 



  §.§.§ Upper Bound of the Sampling Error


For notational convenience, define u_1(x)=u(x-x_1^S)
and v_1(y)=u(x-y_1^S). Consider the sampling position
(x,y)∈𝒢_ij. Using a Taylor expansion, we have 

    |h_1(x,y)-h_1(x_1,y_1)|
         =α|u_1(x)v_1(y)-u_1(x_1)v_1(y_1)|
         =α|(u_1(x_1)+u_1^'(x_1)(x-x_1))
               ×(v_1(y_1)+v_1^'(y_1)(y-y_1))-u_1(x_1)v_1(y_1)
                               +o(x-x_1)+o(y-y_1)|
        =α|u_1(x_1)v_1^'(y_1)(y-y_1)+v_1(y_1)u_1^'(x_1)(x-x_1)|
                               +o(x-x_1)+o(y-y_1)
       ≤α u(0)K_uL/n_c+o(L/n_c)

from the property u(x)≤ u(0) and |u^'(x)|≤ K_u. 

From (<ref>), we have 

    |Ĥ_ij-H_ij|     =(L/n_c)^2|h_1(x,y)-h_1(x_1,y_1)|
       ≤α u(0)K_uL^3/n_c^3+o(L^3/n_c^3).

As a result, 

    𝒫_Ω(𝐇̂_c-𝐇)_F^2    =∑_(i,j)∈Ω|Ĥ_ij-H_ij|^2
       ≤ M(α u(0)K_uL^3/n_c^3)^2≜ϵ̅^2.





  §.§.§ Matrix Completion with Noise and Singular Vector Perturbation


When there is sampling noise, the performance of matrix completion
can be evaluated by the following result. 

Consider that ϵ in (<ref>) is chosen
such that 𝒫_Ω(𝐇̂-𝐇)_F≤ϵ=ϵ̅.
Then, with high probability, 

    δ≜𝐇̂_c-𝐇_F≤4√((2+p)n_c/p)ϵ̅+2ϵ̅

where p=M/n_c^2.



As we focus on not too small n_c, which is chosen to be such
that M≈ Cn_c(log n_c)^2, the bound (<ref>)
can be simplified as 

    δ   ≤4√((2+Cn_c(log n_c)^2/n_c^2)n_c/Cn_c(log n_c)^2/n_c^2)ϵ+2ϵ
       ≈√(32/C)n_c/log n_cϵ.


Let 𝐮_1 and 𝐮̂_1=𝐮_1+𝐞_1
be the dominant left singular vectors of 𝐇 and 𝐇̂_c,
respectively. We exploit the following classical result from singular
vector perturbation analysis. 

Let σ_1 and σ_2 be the first and second dominant
singular values of 𝐇. Then, 

    sin∠(𝐮_1,𝐮̂_1)≤2𝐇̂_c-𝐇_F/σ_1-σ_2=2δ/σ_1-σ_2.




By exploiting Lemma <ref> for our
case, we have 

    sin∠(𝐮_1,𝐮̂_1)     =√(1-|𝐮_1^T(𝐮_1+𝐞_1)|^2)
        =√(-2𝐮_1^T𝐞_1+|𝐮_1^T𝐞_1|^2)
       ≈√(2|𝐮_1^T𝐞_1|)

where |·| denotes the absolute value operator, and we drop
the second order term |𝐮_1^T𝐞_1|^2,
since |𝐮_1^T𝐞_1| is small as we
focus on large n_c(M). We also note that 𝐮_1^T𝐞_1≤0. 

Consider that we have chosen M≈ Cn_c(log n_c)^2, and
moreover, 𝐇 is a rank-1 matrix with singular value σ_1=α.
As a result, 

    2|𝐮_1^T𝐞_1|≈sin^2∠(𝐮_1,𝐮̂_1)≤(2δ/α)^2   ≤4/α^232/C(n_c/log n_c)^2ϵ̅^2
        =128u(0)^2K_u^2L^6n_c^-3.





  §.§.§ Estimator based on Reflected Correlation


Let e(x)=û(x)-u_1(x).  Define a reflected correlation function
as 

    R(t;x_1^S)=∫_-∞^∞u(x-x_1^S)u(-x-x_1^S+t)dx.

Then, it follows that R(t;x_1^S)=τ(2x_1^S-t).
As a result, we have 

    R̂(t;𝐮̂_1)
         =∫_-∞^∞(u_1(x)+e(x))(u_1(-x+t)+e(-x+t))dx
         =∫_-∞^∞u_1(x)u_1(-x+t)dx+∫_-∞^∞u_1(x)e(-x+t)dx
               ∫_-∞^∞e(x)u_1(-x+t)dx+∫_-∞^∞e(x)e(-x+t)dx
         ≈ R(t;x_1^S)+∫_-∞^∞u_1(x)e(-x+t)dx
                   +∫_+∞^-∞e(-y+t)u_1(y)(-dy)
         =R(t;x_1^S)+2∫_-∞^∞u_1(x)e(-x+t)dx
         ≈ R(t;x_1^S)+2𝐮_1^T𝐞_1^-t

where the first approximation (<ref>) is by dropping
the second order term ∫_-∞^∞e(x)e(-x+t)dx, and
the second approximation (<ref>) is to use the inner
product 𝐮_1^T𝐞_1^-t to approximate
the integral based on assumptions A1 and A2 in Section <ref>.
As a result, we have R(t;x_1^S)-R̂(t;𝐮̂_1)≈-2𝐮_1^T𝐞_1^-t.

Recall that t̂=2x̂_1^S maximizes R̂(t̂;𝐮̂_1)
and t^*=2x_1^S maximizes R(t^*;x_1^S)=τ(2x_1^S-t^*).
We have 

    τ(0)-τ(2|x̂_1^S-x_1^S|)
           =R(t^*;x_1^S)-R̂(t̂;𝐮̂_1)
           ≈-2𝐮_1^T𝐞_1^-t
           ≤ C_e2|𝐮_1^T𝐞_1|
           ≤μ_uL^6n^-3+o(n_c^-3)

where μ_u=C_e128u(0)^2K_u^2 and o(n_c^-3) is
due to the fact that we keep omitting the higher order terms. Finally,
we obtain 

    τ(2|x̂_1^S-x_1^S|)=1-μ_uL^6n^-3+o(n_c^-3)

and hence,

    |x̂_1^S-x_1^S|≤1/2τ^-1(1-μ_uL^6n_c^-3+o(n_c^-3)).






 §.§ Proof of Theorem <ref>


We first study the singular vectors in double source case. 
 Let
𝐮_k(θ) and 𝐯_k(θ) be the vectors
defined following (<ref>) and (<ref>) in the rotated
coordinate system 𝒞_θ. The svd of 𝐇_θ
is given by 

    𝐇_θ=α_1𝐩_1𝐪_1^T+α_2𝐩_2𝐪_2^T

where α_1=α/2𝐮_1+𝐮_2𝐯_1+𝐯_2
and α_2=α/2𝐮_1-𝐮_2𝐯_1-𝐯_2
are the singular values, and

    𝐩_1=𝐮_1+𝐮_2/𝐮_1+𝐮_2,  𝐪_1=𝐯_1+𝐯_2/𝐯_1+𝐯_2


    𝐩_2=𝐮_1-𝐮_2/𝐮_1-𝐮_2,  𝐪_2=𝐯_1-𝐯_2/𝐯_1-𝐯_2

are the corresponding singular vectors. 



First, 

    𝐇_θ    =α̅(𝐮_1𝐯_1^T+𝐮_2𝐯_2^T)
        =α̅/2[(𝐮_1+𝐮_2)(𝐯_1+𝐯_2)^T+(𝐮_1-𝐮_2)(𝐯_1-𝐯_2)^T]
        =α_1𝐩_1𝐪_1^T+α_2𝐩_2𝐪_2^T

Hence, these four vectors form a decomposition of 𝐇_θ. 

Second, we have 

    𝐩_1^T𝐩_2    =c(𝐮_1+𝐮_2)^T(𝐮_1-𝐮_2)
        =c(𝐮_1^2-𝐮_2^2)
        =0

where c=1/(𝐮_1+𝐮_2𝐮_1-𝐮_2).
Similarly, 𝐪_1^T𝐪_2=0. In addition,
all the four vectors have unit norm. 

As a result, (<ref>) is the svd of 𝐇_θ.



Consider an arbitrary coordinate system. wlog (due to Assumption
1), assume that the first source is located at the origin, x_1^S=0
and y_1^S=0, and the second source is away from the
first source with distance D and angle θ to the x-axis,
x_2^S=Dcosθ and y_2^S=Dsinθ.
In addition, defining 

    u_c(x,θ)≜ u(x-Dcosθ),     u_s(x,θ)≜ u(x-Dsinθ)

we have 

    𝐮_1    =√(δ)[u(x_1),u(x_2),…,u(x_N)]^T
    𝐯_1    =√(δ)[u(y_1),u(y_2),…,u(y_M)]^T
    𝐮_2    =√(δ)[u_c(x_1,θ),u_c(x_2,θ),…,u_c(x_N,θ)]^T
    𝐯_2    =√(δ)[u_s(y_1,θ),u_s(y_2,θ),…,u_s(y_M,θ)]^T.


Based on assumption A1 and A2, we have 

    𝐮_k^2=(L/n)^2∑_i=1^Nu(x_i-x_k^S)^2   ≈∫_x_1^x_n-1u(x-x_k^S)^2dx
       ≈∫_-∞^∞u(x-x_k^S)^2dx=1

and similar integrals apply to 𝐯_k. 

As an equivalent statement to Theorem <ref>,
we need to show that ρ(θ) is a strictly increasing function
in θ∈(0,π/4). Equivalently, we should prove that
the function 

    λ_2(𝐇_θ)^2/λ_1(𝐇_θ)^2
         ≈∫_-∞^∞(u(x)-u_c(x,θ))^2dx/∫_-∞^∞(u(x)+u_c(x,θ))^2dx∫_-∞^∞(u(x)-u_s(x,θ))^2dx/∫_-∞^∞(u(x)+u_s(x,θ))^2dx
         ≜μ(θ)

is strictly increasing in θ∈(0,π/4), where the
approximated integrals are obtained from (<ref>).

To simplify the notation, define the integration operator ⟨·⟩
as 

    ⟨ f⟩≜∫_-∞^∞f(x,θ)dx

for a function f(x,θ). By definition, the integration operator
is linear and satisfies the additive property, i.e., ⟨ af⟩=a⟨ f⟩
and ⟨ f+g⟩=⟨ f⟩+⟨ g⟩, for a
constant a and a function g(x,θ). As a result, ⟨(u-u_c)^2⟩=⟨ u^2⟩+⟨ u_c^2⟩-2⟨ u· u_c⟩=2(1-⟨ u· u_c⟩),
and the function μ(θ) can be written as 

    μ(θ)=(1-⟨ u· u_c⟩)(1-⟨ u· u_s⟩)/(1+⟨ u· u_c⟩)(1+⟨ u· u_s⟩).


In addition, from the properties in calculus, if f(x,θ) and
∂/∂θf(x,θ) are continuous in θ,
then 

    d/dθ⟨ f⟩    =d/dθ∫_-∞^∞f(x,θ)dx
        =∫_-∞^∞∂/∂θf(x,θ)dx=⟨∂/∂θf⟩.

Therefore, defining 

    u_c^'(x,θ)    ≜d/dxu(x)|_x=x-Dcosθ
    
    u_s^'(x,θ)    ≜d/dxu(x)|_x=x-Dsinθ

we have

    d/dθ⟨ u· u_c⟩    =⟨ u·∂/∂θu_c(x,θ)⟩=⟨ u· u_c^'⟩ Dsinθ
    d/dθ⟨ u· u_s⟩    =⟨ u·∂/∂θu_s(x,θ)⟩=-⟨ u· u_s^'⟩ Dcosθ.


With some algebra, the derivative of μ(θ) can be obtained
as 

    d/dθμ(θ)     =η[Dcosθ⟨ u· u_s^'⟩(1-⟨ u· u_c⟩^2)
               -Dsinθ⟨ u· u_c^'⟩(1-⟨ u· u_s⟩^2)]
        =η[-t·τ^'(s)(1-τ(t)^2)+s·τ^'(t)(1-τ(s)^2)]

where η=2(1+⟨ u· u_c⟩)^-2(1+⟨ u· u_s⟩)^-2,
t=Dcosθ, and s=Dsinθ. 

Note that 0<s<t for 0<θ<π/4. Applying condition
(<ref>), we have 

    d/dθμ(θ)     >η· t·τ^'(s)[(1-τ(s)^2)-(1-τ(t)^2)]
        =η· t·τ^'(s)(τ(t)^2-τ(s)^2)
        >0

since τ^'(s)<0 and τ(t)<τ(s) for 0<s<t.

This confirms that μ(θ) is a strictly increasing function,
and hence ρ(θ) is a strictly increasing function in θ∈(0,π/4).
The results in Theorem <ref> is confirmed. 

IEEEtran
10

url@samestyle








BecStoLi:J08
A. Beck, P. Stoica, and J. Li, “Exact and approximate solutions of source
  localization problems,” IEEE Trans. Signal Process., vol. 56,
  no. 5, pp. 1770–1778, 2008.

QiXiuYua:J13
H.-D. Qi, N. Xiu, and X. Yuan, “A lagrangian dual approach to the
  single-source localization problem,” IEEE Trans. Signal Process.,
  vol. 61, no. 15, pp. 3815–3826, 2013.

SheHu:J05
X. Sheng and Y.-H. Hu, “Maximum likelihood multiple-source localization using
  acoustic energy measurements with wireless sensor networks,” IEEE
  Trans. Signal Process., vol. 53, no. 1, pp. 44–53, 2005.

MeeMitNar:J08
C. Meesookho, U. Mitra, and S. Narayanan, “On energy-based acoustic source
  localization for sensor networks,” IEEE Trans. Signal Process.,
  vol. 56, no. 1, pp. 365–377, 2008.

LiuHuPan:J12
Y. Liu, Y. H. Hu, and Q. Pan, “Distributed, robust acoustic source
  localization in a wireless sensor network,” IEEE Trans. Signal
  Process., vol. 60, no. 8, pp. 4350–4359, 2012.

ZisWax:C88
I. Ziskind and M. Wax, “Maximum likelihood localization of multiple sources by
  alternating projection,” Proc. IEEE Int. Conf. Acoustics, Speech, and
  Signal Processing, vol. 36, no. 10, pp. 1553–1560, 1988.

LefRea:J17
R. Lefort, G. Real, and A. Drémeau, “Direct regressions for underwater
  acoustic source localization in fluctuating oceans,” Applied
  Acoustics, vol. 116, pp. 303–310, 2017.

NguJorSin:J05
X. Nguyen, M. I. Jordan, and B. Sinopoli, “A kernel-based learning approach to
  ad hoc sensor network localization,” ACM Trans. on Sensor Networks,
  vol. 1, no. 1, pp. 134–152, 2005.

JinSohWon:J10
Y. Jin, W.-S. Soh, and W.-C. Wong, “Indoor localization with channel impulse
  response based fingerprint and nonparametric regression,” IEEE
  Trans. Wireless Commun., vol. 9, no. 3, pp. 1120–1127, 2010.

KimParYooKimPar:J13
W. Kim, J. Park, J. Yoo, H. J. Kim, and C. G. Park, “Target localization using
  ensemble support vector regression in wireless sensor networks,” IEEE
  Trans. on Cybernetics, vol. 43, no. 4, pp. 1189–1198, 2013.

ChoMit:C15
S. Choudhary and U. Mitra, “Analysis of target detection via matrix
  completion,” in Proc. IEEE Int. Conf. Acoustics, Speech, and Signal
  Processing, 2015, pp. 3771–3775.

ChoKumNarMit:J16
S. Choudhary, N. Kumar, S. Narayanan, and U. Mitra, “Active target
  localization using low-rank matrix completion and unimodal regression,”
  arXiv preprint arXiv:1601.07254, 2016.

CanPla:J10
E. J. Candes and Y. Plan, “Matrix completion with noise,” Proceedings
  of the IEEE, vol. 98, no. 6, pp. 925–936, 2010.

peterson:b95
R. Peterson, R. Ziemer, and D. Borth, Introduction to spread spectrum
  systems.1em plus 0.5em minus 0.4emEnglewood Cliffs, NJ:
  Prentice-Hall, 1995.

vu:J11singular
V. Vu, “Singular vectors under random perturbation,” Random Structures
  & Algorithms, vol. 39, no. 4, pp. 526–538, 2011.




