# Insight of the Green’s function as a defect state in a boundary value
problem

Jose D. H. Rivero  Li Ge  College of Staten Island, CUNY, Staten Island, New
York 10314, USA  The Graduate Center, CUNY, New York, NY 10016, USA

###### Abstract

A new perspective of the Green’s function in a boundary value problem as the
only eigenstate in an auxiliary formulation is introduced. In this treatment,
the Green’s function can be perceived as a defect state in the presence of a
$\delta$-function potential, the height of which depends on the Green’s
function itself. This approach is illustrated in one-dimensional and two-
dimensional Helmholtz equation problems, with an emphasis on systems that are
open and have a non-Hermitian potential. We then draw an analogy between the
Green’s function obtained this way and a chiral edge state circumventing a
defect in a topological lattice, which shines light on the local minimum of
the Green’s function at the source position.

##  I Introduction

Initially distributed to just 51 private members of a subscription library in
1828 [1], the essay of George Green on the application of mathematical
analysis to the theories of electricity and magnetism has inspired generations
of physicists and mathematicians. By studying the three-dimensional Laplace
equation governing the static electric potential, Green recognized the unique
utility of a function inversely proportional to the distance between two
charged bodies, which now bears his name and preceded Dirac’s introduction of
the $\delta$-function by more than a century. The importance of the Green’s
function, together with the Green’s theorem, was recognized by masterminds of
modern physics such as Lord Kelvin and Julian Schwinger.

Besides its effectiveness as a theoretical and numerical tool to solve linear
differential equations [2, 3], the Green’s function provides a fundamental
connection between the source and the field [4]. When extended to address
dynamical equations and their spectral representations, the Green’s function
becomes the propagator that is fundamental to field theories [5]. As such, it
is a powerful and essential tool to study a variety of transport and
scattering problems, ranging from condensed matter physics [6, 7], optics and
photonics [8, 9, 10] to high energy physics [11].

In optical physics, one area of interest in this article, the Green’s function
has been routinely employed to analyze rich optical phenomena including light-
matter interaction in lasers and other light sources [12, 13, 14], as well as
to devise applications for optical computing, information processing, and
telecommunications. What we offer here is a novel perspective of the Green’s
function, i.e., treating it as the single eigenstate in an auxiliary boundary
value problem. In addition to further enrichment and shaping of our physical
intuition through the Green’s function, we find exceptional parallels between
the Green’s function and defect states due to a local potential, including a
chiral edge state circumventing a defect on its path in a topological lattice.

Below we introduce formally the Green’s function problem and lay out the
fundamentals of our approach. The Green’s function of an operator
$\mathcal{L}$ in a variety of physics problems can be defined by

$[z-\mathcal{L}]G(\bm{r},\bm{r^{\prime}};z)=\delta(\bm{r}-\bm{r^{\prime}})$
(1)

with a proper boundary condition. Here $z$ is a parameter (e.g., the energy)
and $\bm{r},\bm{r^{\prime}}$ are the coordinates of the field and the source,
respectively. A new perspective of the Green’s function can be obtained using
the following auxiliary eigenvalue problem

$[z-\mathcal{L}]\psi_{m}(\bm{r},\bm{r^{\prime}})=\lambda_{m}(\bm{r^{\prime}})f(\bm{r},\bm{r^{\prime}})\psi_{m}(\bm{r},\bm{r^{\prime}}).$
(2)

Below we suppress the $\bm{r^{\prime}}$-dependence of $\psi_{m}$ and
$\lambda_{m}$ for conciseness. As we will show, this auxiliary problem has a
single eigenvalue $\lambda_{0}$ and eigenstate $\psi_{0}(\bm{r})$, when
$f(\bm{r},\bm{r^{\prime}})$ is chosen to be $\delta(\bm{r}-\bm{r^{\prime}})$.
The Green’s function is then uniquely determined by

$G(\bm{r},\bm{r^{\prime}})=\frac{\psi_{0}(\bm{r})}{\lambda_{0}\psi_{0}(\bm{r^{\prime}})},$
(3)

and $\lambda_{0}^{-1}$ gives the value of the Green’s function at the source
$\bm{r}=\bm{r^{\prime}}$. As a bonus, we obtain directly the local density of
states (LDOS) that is proportional to the imaginary part of
$G(\bm{r},\bm{r};z)$ [7], i.e., $\text{Im}[\lambda_{0}^{-1}]$. The reciprocity
of the Green’s function, though implicit in Equation (3), can be easily
verified in the absence of an effective magnetic field, as shown in Appendix
A].

Before we apply this approach to various Hermitian and non-Hermitian problems,
we first prove Equation (3) and discuss the general properties of
$\psi_{m}(\bm{r})$. With the choice of $f(\bm{r},\bm{r^{\prime}})$ mentioned
above, the right hand side of Equation (2) for $\psi_{0}(\bm{r})$ is simply

$\lambda_{0}\delta(\bm{r}-\bm{r^{\prime}})\psi_{0}(\bm{r})=\lambda_{0}\delta(\bm{r}-\bm{r^{\prime}})\psi_{0}(\bm{r^{\prime}}),$
(4)

from which Equation (3) follows directly by comparing with Equation (1), with
the requirement $\psi_{0}(\bm{r^{\prime}})\neq 0$. If there were another
eigenstate $\psi_{m\neq 0}(\bm{r})$, then by repeating the same procedure we
would find that the Green’s function is proportional to $\psi_{m\neq
0}(\bm{r})$ as well, which contradicts the uniqueness of the Green’s function
with a properly imposed boundary condition. When implemented numerically,
there do exist spurious eigenvectors $\psi_{m\neq 0}(\bm{r})$, which however
can be easily discarded due to their ill-behaved $\lambda_{m}$’s as we will
discuss in Section II.1.

The auxiliary eigenvalue approach equips us with a conceptually different way
to treat the Green’s function, i.e., as a defect state [15] emerging due to
the $\delta$-function potential:

$\left[\mathcal{L}+\lambda_{0}\delta(\bm{r}-\bm{r^{\prime}})\right]\psi_{0}(x)=z\psi_{0}(x).$
(5)

As we will show, this point of view is particularly interesting and helpful in
a topological system with chiral edge states [16, 17, 18, 19, 20]. For
example, if a point source is placed at the edge of a two-dimensional (2D)
topological insulator, the auxiliary eigenvalue approach indicates that an
analogy exists between the Green’s function and a chiral edge state
circumventing a defect at the same location. This interpretation provides an
intuitive understanding of the local minimum of the Green’s function at the
source position, which we will illustrate using a 2D square lattice with a
$\pi/2$ Landau gauge.

The rest of the paper is organized into two main parts, where we validate our
method and discuss the new insight it provides, respectively. In Section II,
we first validate our method in a 1D Hermitian (closed) system where the
analytical form of the Green’s function is available. We then extend the
validation to two non-Hermitian 1D systems with parity-time ($\cal PT$)
symmetry [21, 22], focusing on the Green’s function at an exceptional point
(EP). An EP is a unique degeneracy found only in non-Hermitian systems, where
two or more eigenstates of the system coalesce. It has led to a plethora of
intriguing phenomena, and in particular, it was shown recently that the
Green’s function can be fully decoupled from the coalesced eigenstate in a
photonic system, which is instead given by the Jordan vector or the “missing
dimension” of the Hilbert space [23]. We show that our method based on
Equation (2) captures this extraordinary behavior nicely in a ring cavity,
besides describing correctly the Green’s function in a $\cal PT$-symmetric
photonic molecule. We further validate our method in quasi-1D waveguides,
which are frequently employed to study disordered mesoscopic and optical
systems [24, 25, 26]. In Section III, we first discuss how the defect state
corresponding to the Green’s function is conceptually different from previous
studies of defect states in a 1D continuous Hermitian system and a 1D non-
Hermitian lattice. We then highlight the intriguing manifestation of the
linkage between the Green’s function and a defect state in the aforementioned
2D topological lattice.

##  II Validation

###  II.1 1D Hermitian case

Figure 1: The Green’s function (a) and DOS (b) in a 1D dielectric cavity with
perfect mirrors. The green dashed lines are obtained using the analytical
expression (8), while the black solid lines are from our auxiliary eigenvalue
approach (3). $n=3$, $k=30/L$ and $x^{\prime}=0.45L$ (marked by the arrow) are
used in (a), and 2000 grid points are used for the finite difference
implementation of Equation (2).

We start with an exmaple where the analytical form of the Green’s function is
available. Consider the scalar Helmholtz equation in 1D with a uniform
refractive index $n\in\mathbb{R}$:

$\mathcal{L}=-\frac{1}{n^{2}}\partial_{x}^{2},\quad z=k^{2}.$ (6)

Here $k$ is the real-valued wave vector in free space. Below we take the speed
of light in vacuum to be 1 and do not distinguish $k$ from the (circular)
frequency. We impose the Dirichlet boundary conditions $G(x,x^{\prime})=0$,
$\psi_{m}(x)=0$ at $x=0,L$, which renders the system Hermitian. Consequently,
it can be shown that the Green’s function at the source position (i.e.,
$\lambda_{0}^{-1}$) is real:

$\lambda_{0}=\frac{k}{n}\left[\frac{1}{\tan[nk(x^{\prime}-L)]}-\frac{1}{\tan(nkx^{\prime})}\right].$
(7)

When solved using a finite difference scheme [27], the numerical value of
$\lambda_{0}$ given by Equation (3) shows a good agreement with the analytical
result given by Equation (7). The corresponding Green’s function, which we
shown in Figure 1(a), also agrees nicely with its analytical expression

$G(x,x^{\prime})=\begin{cases}\frac{\sin[nk(x-L)]}{\lambda_{0}\sin[nk(x^{\prime}-L)]},&(x>x^{\prime})\\\
\frac{\sin(nkx)}{\lambda_{0}\sin(nkx^{\prime})}.&(x\leq
x^{\prime})\end{cases}$ (8)

The numerical implementation of the $\delta$-function is usually taken to be
the limit of a sharp analytical distribution, such as a Gaussian with the
standard deviation $\sigma\rightarrow 0$. At the same time, another small
quantity that requires attention in the finite difference implementation of
Equation (2) is the lattice spacing $\Delta x$. If we choose $\sigma>\Delta
x$, we find spurious eigenstates $\psi_{m\neq 0}(x)$ with almost identical
spatial dependence away from the source position (not shown), but they have
structures that reside in the finite extension of the approximated
$\delta$-function (e.g., fast oscillations). This problem is remedied by
letting $\sigma<\Delta x$, which practically leads to an approximation of the
$\delta$-function that has a single non-zero element at the position of the
source, the value of which is given by $1/\Delta x$. This choice warrants that
the integration of the $\delta$-function is 1 over any range enclosing the
source and 0 otherwise. With this choice, we find that all spurious
eigenvalues of Equation (2) approaches infinity (e.g., $|\lambda_{m\neq
0}L|>10^{17}$ in the case shown in Figure 1), and the sole, physical one
$\lambda_{0}$ is easily obtained by setting the numerical routine to search
for the eigenvalue with the smallest modulus.

As mentioned in the introduction, our auxiliary eigenvalue approach also
produces the LDOS directly. For a Hermitian system, a limiting procedure is
needed to regulate the singularity at real-valued resonant frequencies:

$\text{LDOS}(x;k)=\lim_{s\to 0^{+}}-\frac{2k}{\pi}\text{Im}[G(x,x;k+is)],$ (9)

i.e., a small positive imaginary part is added to the frequency here,
constructing in this way the retarded Green’s function [See Appendix A].
Equation (9) also applies to non-Hermitian cases, as long as the complex
resonant frequencies are on or below the real axis. The integration of
LDOS$(x;k)$ over the whole system then gives the density of states (DOS) as a
function of the frequency.

To calculate DOS using our approach based on Equation (2), we choose a small
$s=0.03/L$ and calculate $-\frac{2k}{\pi}\text{Im}[\lambda_{0}^{-1}]$
numerically. The result agrees well with the analytical result [see Figure
1(b)], where $\lambda_{0}$ given by Equation (7) is used. The latter leads to

$\text{DOS}(k)=\sum_{m}\delta\left(k-\frac{m\pi}{nL}\right)$ (10)

as expected once $s\to 0$, where $k_{m}=m\pi/nL\,(m=1,2\ldots)$ are the real-
valued resonant frequencies.

###  II.2 1D non-Hermitian cases

For our next validation, we study the Green’s function at an EP in a photonic
molecule [28]. Such a case presents a serious challenge to the standard
approach based on the eigenvalues of $\cal L$, i.e., the bilinear expansion

$\displaystyle G(\bm{r},\bm{r^{\prime}};z)=\sum_{m}$
$\displaystyle\frac{\overline{\phi}_{m}(\bm{r})\phi_{m}(\bm{r^{\prime}})}{(z-z_{m})(m,m)},$
(11) $\displaystyle\mathcal{L}\phi_{m}(\bm{r})$
$\displaystyle=z_{m}\phi_{m}(\bm{r}).$ (12)

Depending on the symmetry of $\mathcal{L}$, the partner function
$\overline{\phi}_{m}(\bm{r})$ may or may not be the same as
$\phi_{m}^{*}(\bm{r})$, and $(m,n)$ is the resulting inner product of
$\overline{\phi}_{m}(\bm{r})$ and $\phi_{n}(\bm{r})$. At an EP, the inner
product $(m,m)$ becomes zero due to the coalescence of two or more eigenstates
of the system [29, 30, 22]. Although this divergence can be eliminated using a
non-Hermitian perturbation theory [10, 23], it requires a priori knowledge of
the EPs and the generalized eigenvectors, which adds to the complexity of the
problem. Our auxiliary eigenvalue formulation, on the other hand, does not
suffer from this drawback.

Below we exemplify our method and compare it to the result of the perturbation
theory, without the need to invoke the Jordan vector that completes the
Hilbert space at the EP [23]. Our photonic molecule is composed of two half-
wavelength cavities coupled by a distributed Bragg reflector (DBR) placed in
air [See Figure 2(a); inset]. If an unbalanced loss distribution is introduced
in the two half-wavelength cavities, an emerging effective $\cal PT$ symmetry
governs the system [31].

Here we consider the outgoing boundary condition at a real-valued frequency,
which corresponds to the laser and is different from quasi-bound states or
resonances with a complex frequency throughout the entire space [12]. The
resulting complex eigenvalues of $\mathcal{L}$ given in Equation (6) are found
in the lower half of the complex frequency plane, as a result of the non-
Hermiticity caused by the cavity openness. These complex eigenvalues are known
as continuous flux (CF) states [12, 27], which have also been used to study
nuclear decays [32, 33]. We note that the EPs of CF states have not been
studied before, unlike their counterparts of quasi-bound states or resonances.

Figure 2: Green’s function of a $\mathcal{PT}$ photonic molecule at an EP. The
source is placed at the center of the right cavity in (a) and in between the
two cavities in (b). Inset in (a): Schematic of the photonic molecule and the
imposed boundary condition, and labels denoting the refractive indices of each
component of the heterostructure. Green dashed lines correspond to the
perturbative expression in Equation (16) and black solid lines correspond to
our auxiliary eigenvalue approach. The locations of the sources are marked
with arrows. (c) Difference between the bilinear expansion with the
perturbative correction and the auxiliary method in (b) at the source
position.

In the vicinity of an EP of frequency $k_{0}$, the CF states in our system can
be expressed in terms of waves confined in the left- and right-cavity, i.e.,
$\psi(x)\approx a_{l}\psi_{l}(x)+a_{r}\psi_{r}(x)$, where the amplitudes
$a_{l,r}$ are determined by the Helmholtz equation. Without the
$\mathcal{PT}$-symmetric perturbation, our photonic molecule is symmetric and
hence $a_{l}=\pm a_{r}$ in the symmetric and anti-symmetric modes with CF
frequencies $\tilde{k}_{S},\tilde{k}_{A}$. The introduction of a weak
$\mathcal{PT}$-symmetric perturbation in the dielectric function couples the
amplitudes $a_{l}$ and $a_{r}$, determined by their spatial overlap $C$ with
the non-Hermitian perturbation, which represents the strength of gain and
loss.

The eigenfrequencies of the perturbed system are then found to be
$q_{\pm}^{2}=k_{0}^{2}\left(1\pm\Delta\varsigma\right)$. Here $k_{0}$ is the
CF frequency of a single half-wavelength cavity sandwiched by two DBRs and
$\Delta$ is a dimensionless detuning defined by
$(\tilde{k}_{S}^{2}-\tilde{k}_{A}^{2})/2k_{0}^{2}$. We have also defined
$\varsigma\equiv\sqrt{1-\beta^{2}}$, where $\beta=C/\Delta$. One EP is reached
when $\varsigma=0$, resulting in $q^{2}_{\pm}=q_{0}^{2}$. The corresponding
eigenstates $\psi_{\pm}(x)$ are given by

$\psi_{\pm}(x)=\psi_{l}(x)+\beta_{\pm}\psi_{r}(x),$ (13)

where $\beta_{\pm}\equiv i\beta\pm\varsigma$. The inner product of the
eigenstates is defined as

$(i,j)\equiv\int\epsilon_{0}(x)\psi_{i}(x)\psi_{j}(x)dx,\quad(i,j=\pm)$ (14)

where $\epsilon_{0}(x)$ is the dielectric function before the $\cal
PT$-symmetric perturbation is introduced. In other words, the partner
functions in Equation (11) are chosen as
$\overline{\psi}_{\pm}(x)=\epsilon_{0}(x)\psi_{\pm}(x)$, similar to the
Hermitian case discussed earlier. This definition of the inner product
warrants the biorthogonality $(+,-)=(-,+)=0$, and we find
$(+,+)=2\beta_{+}\varsigma\to 0$, $(-,-)=-2\beta_{-}\varsigma\to 0$ as the
system approaches the EP.

As seen from Equation (11), the vanishing inner products $(+,+),\,(-,-)$ here
at the EP cause a catastrophe in the calculation of the Green’s function,
because the two corresponding terms diverge independent of the frequency:

$G(x,x^{\prime};k)\approx\frac{\overline{\psi}_{+}(x^{\prime})\psi_{+}(x)}{(q^{2}-q_{+}^{2})(+,+)}+\frac{\overline{\psi}_{-}(x^{\prime})\psi_{-}(x)}{(q^{2}-q_{-}^{2})(-,-)}.$
(15)

However, it can be shown that the diverging behaviors in the two terms cancel
each other precisely [10, 23], leading to

$\displaystyle G(x,x^{\prime};k)\approx$
$\displaystyle\;\epsilon_{0}(x^{\prime})\frac{\psi_{l}(x^{\prime})\psi_{l}(x)+\psi_{r}(x^{\prime})\psi_{r}(x)}{k^{2}-k_{0}^{2}}$
$\displaystyle-i\Delta
k_{0}^{2}\epsilon_{0}(x^{\prime})\frac{\psi_{EP}(x^{\prime})\psi_{EP}(x)}{(k^{2}-k_{0}^{2})^{2}}$
(16)

to the leading order of the small perturbation parameter $\varsigma$. Here
$\psi_{EP}=\psi_{\pm}|_{\beta=1}$ is the coalesced eigenstate at the EP. As
expected [9, 34], a second-order pole appears in the second term due to this
coalescence. Details on the perturbative analysis can be found in Appendix B.

To verify the robustness of our method based on Equation (3) in the vicinity
of the EP, we choose the heterostructure of length $L$ shown in Figure 2. It
consists of DBRs of refractive indexes $n_{1}=2+0.001i$ and $n_{2}=3+0.001i$,
and each layer accommodates a quarter of their respective wavelengths at
$ka=1.570$, where $a$ is the lattice constant. The two half-wavelength
cavities with loss are fine-tuned to achieve an EP at
$k_{\text{EP}}a=1.570-0.006913i$ using
$n_{l}=2+0.001i$,$n_{r}=1.9999+0.0029i$, and $k=\text{Re}[k_{\text{EP}}]$. We
note that the loss in the right cavity is stronger than that in the left
cavity, creating an effective $\mathcal{PT}$-symmetric system with its complex
eigenvalue spectrum shifted along the imaginary axis [35].

Figure 2 shows the Green’s function of this system when the source is placed
at different locations. First, we place the source at the center of the right
cavity, where we expect the approximation (16) using just $\psi_{\pm}(x)$ [and
$\psi_{l,r}(x)$] to hold. As Figure 2(a) shows, it is nearly identical to the
result of our auxiliary eigenvalue approach (3), and the inclusion of more CF
states away from the EP barely changes the Green’s function (not shown). If
instead, we place the source at the center of the heterostructure, Equation
(16) alone is insufficient to capture the features of the Green’s function
[see Figure 2(b)]. However, once a large set of additional eigenfunctions of
$\mathcal{L}$ are also included, a good agreement to our auxiliary eigenvalue
approach is again observed [see Figure 2(c)]. These results show that Equation
(3) provides a reliable and convenient method to calculate the Green’s
function, even in the presence of EPs.

Figure 3: Chirality-reversal Green’s function in a $\cal PT$-symmetric ring
cavity. (a) Phase of the Green’s function (solid line) and the coalesced
eigenstate (dotted line). Dashed lines point to the phase of the Green’s
function at the source (marked by the arrow). Shaded regions show the half
periods with higher loss. (b) False color plots showing the constant
amplitudes of the coalesced eigenstate (top) and the Green’s function
(bottom). A finite width is imposed on the ring for visual clarity. Here
$n_{0}=3$, $\delta n=0.003$, $l=1$, $k_{\text{EP}}L=2.0944-0.0021i$,
$kL=2.0944$, and $x^{\prime}=L/8$.

A similar but more striking behavior of the Green’s function at an EP was
recently reported in an effective $\cal PT$-symmetric ring cavity with
refractive index [23]

$n(\theta)=(n_{0}+i\delta n)+\delta n(\cos 2l\theta+i\sin 2l\theta).$ (17)

Here $\theta$ is the azimuthal angle, and below we will use the arc length
$x=R\theta\in[0,L]$ as the coordinate, where $R$ is the radius of the ring and
$L=2\pi R$ is the circumference. $n_{0}$ and $\delta n$ are the real and
imaginary parts of the background index, and the latter includes both
absorption and radiation losses and is positive. $l$ is a positive integer,
and the complex index grating, proportional to $e^{2il\theta}$, scatters the
clockwise (CW) wave of angular momentum $-l$ to the counterclockwise (CCW)
wave of angular momentum $l$ but not vice versa. Note that the chirality and
the sign of $\text{Im}[n]$ here are defined with respect to the temporal
dependence $e^{-i\omega t}$.

Consequently, an EP appears at $k_{\text{EP}}R=l/(n_{0}+i\delta n)$ with the
coalesced CCW eigenstate $\psi(x)=e^{il\theta}$ [36]. The Green’s function, on
the other hand, can be fully decoupled from this mode even on resonance, if
the source is placed at $\theta=(m+1/4)\pi/l$ where $m$ is a non-negative
integer [23], i.e., at one of the most lossy spots in the passive ring cavity.
The Green’s function is given by the corresponding Jordan vector $J(x)\propto
e^{-il\theta}$ instead, i.e., the “missing dimension” of the Hilbert space at
the EP in the CW direction.

This extraordinary behavior is captured nicely using our auxiliary eigenvalue
approach (see Figure 3). In addition, a perturbative approach shows that the
Green’s function at the source is given by [23]

$G(x^{\prime},x^{\prime};k)\approx\frac{1}{2(k-k_{\text{EP}})k_{\text{EP}}L},$
(18)

which is almost imaginary on resonance (i.e., $k=\text{Re}[k_{\text{EP}}]$)
for a high-Q resonance with
$|\text{Im}[k_{\text{EP}}]|\ll\text{Re}[k_{\text{EP}}]$. In the case shown in
Figure 3, this value is $(-113.99i+0.228)R$ and nearly identical to that given
by our auxiliary eigenvalue approach, i.e., $(-113.99i+0.224)R$.

###  II.3 Quasi-1D waveguides

Quasi-1D waveguides are frequently used in the study of disordered mesoscopic
and optical systems [24, 25, 26], and the Green’s function plays a crucial
role to construct the scattering and transfer matrices [37, 38]. Here again we
consider the scalar Helmholtz equation in a waveguide with background
refractive index $n(\bm{r})$:

$\mathcal{L}=-\frac{1}{n^{2}(\bm{r})}(\partial_{x}^{2}+\partial_{y}^{2}),\quad
z=k^{2}.$ (19)

A finite width $L_{y}$ in the transverse direction and the Dirichlet boundary
conditions at $y=0,L_{y}$ lead to a set of transverse modes (“channels”)
$f_{m}(y)=\sin[nk_{m}^{(y)}y]$, where $nk_{m}^{(y)}=m\pi/L_{y}$ is the
transverse wave number and $m$ is a positive integer. At a given frequency
$k$, the longitudinal wave number in the $m$th channel is given by
$nk_{m}^{(x)}=n\\{k^{2}-[k_{m}^{(y)}]^{2}\\}^{1/2}$, and this channel is
propagating (evanescent) if $k_{m}^{(x)}$ is real (imaginary).

Figure 4: Green’s function in a uniform waveguide with $L_{y}=2\,\mu$m and
$n=1.5$ everywhere. (a) and (b) show the slices along $x$ and $y$ at the
source position $x^{\prime}=0,\,y^{\prime}=3/13\,\mu$m. The free-space
wavelength is $1550$ nm.

To validate our auxiliary eigenvalue approach in quasi-1D waveguides, we first
consider a uniform waveguide with the outgoing boundary condition. In this
case an analytical expression exists for the Green’s function, which can be
written as the following infinite sum using the channel functions:

$G(\bm{r},\bm{r^{\prime}};k)=\sum_{m}n\frac{\sin[nk_{m}^{(y)}y^{\prime}]\sin[nk_{m}^{(y)}y]}{ik_{m}^{(x)}L_{y}}e^{\pm
ink_{m}^{(x)}x}.$ (20)

Here the source point $\bm{r^{\prime}}=(x^{\prime},y^{\prime})$ is placed at
$x^{\prime}=0$, $y^{\prime}\neq 0,L_{y}$, and the “$+\,(-)$” sign in the
exponent applies to a positive (negative) $x$. Clearly, only the propagating
channels of a finite number affect the farfield behavior of the Green’s
function, while the logarithmic divergence of the Green’s function at the
source, a generic property in 2D (including quasi-1D), is reflected by the
infinite number of evanescent channels in the summation. Numerically, this
divergence is truncated either by the inclusion of only a finite number of
evanescent channels or the finite resolution of the spatial discretization. We
also note that the Green’s function is dimensionless in quasi-1D.

Figure 4 shows one example where our auxiliary eigenvalue approach,
implemented using the finite difference method [39, 27] (solid line), is
compared with the analytical result given by Equation (20). With just the
three propagating channels available in this case (dotted line), Equation (20)
describes the Green’s function well far from the source ($>2\,\mu{m}$); with
more channels included (e.g., 100; dashed line), a good agreement between
Equation (20) and our auxiliary eigenvalue approach is observed, where the
grid spacings $\Delta x=1/30\,\mu$m, $\Delta y=3/130\,\mu$m used in the finite
difference method are comparable to the shortest evanescent tail in the
summation (i.e., $1/k_{m=100}^{(x)}\approx 0.01\,\mu$m.)

For structured or disordered quasi-1D waveguides connected to two semi-
infinite regions (“leads”), the Helmholtz equation is no longer separable in
$x$ and $y$, and an analytical expression for the Green’s function in the form
of Equation (20) does not exist. Nevertheless, one can still compare to the
bilinear expansion (11) as we show in Figs. 5(a) and (b) for a $\cal
PT$-symmetric waveguide, despite its slow convergence [see Figure 5(c)].

Figure 5: (a,b) Same as Figure 4 but for the Green’s function in a $\cal
PT$-symmetric waveguide. Here $L_{x}=9\,\mu$m, $L_{y}=3\,\mu$m, and the source
is at $x^{\prime}=-1/60\,\mu$m, $y^{\prime}=1.5\,\mu$m. Inset in (a):
Schematic of the waveguide. $\text{Im}[n(x)]=\pm 0.05$ for $x\in[-3,0]\mu$m
and $[0,3]\mu$m respectively, where $\text{Re}[n(x)]=2$ in the central half of
the width. $\text{Re}[n(x)]=1.2$ elsewhere inside the waveguide and $1$
outside. (c) Convergence of the bilinear expansion at the source (upper) and
at $x=4.5\,\mu$m, $y=1.5\,\mu$m (lower).

Here we also briefly review a standard technique used to calculate the Green’s
function in quasi-1D waveguides, i.e., the recursive Green’s function method
[40, 41]. The comparison of the Green’s functions obtained in our proposed
approach and by the recursive Green’s function method will be present in
Section III.

The heart of the recursive Green’s function method lies in the celebrated
Dyson’s equation:

$\bm{G}=\bm{G}^{0}+\bm{G}^{0}\bm{VG}.$ (21)

$\bm{G}$ is the Green’s matrix expressed in the spatial basis here, i.e., its
element $\bm{G}_{\underline{i,j}}$ gives the value of the Green’s function at
point $i$ when the source is placed at point $j$. In a quasi-1D system with
$N$ segments, $\bm{V}$ represents the couplings between the $(N-1)$th and
$N$th segments due to $\mathcal{L}$, and $\bm{G}^{0}$ is the value of $\bm{G}$
when $\bm{V}$ is taken to be zero. In other words, $\bm{G}^{0}$ contains the
Green’s matrices of two separate systems, one for the first $N-1$ segments on
the left as a whole and one for the last segment on the right. This recursive
procedure then starts with a single segment on the left and proceeds with more
segments added from the right, one at a time.

When implemented using the finite difference method that is equivalent to a
tight-binding lattice [41], the recursive Green’s function gives the identical
result as our auxiliary eigenvalue approach, because they solve exactly the
same discretized equation. Each segment of the system is simply a single
column along $y$. We note, however, these two methods serve different
purposes. The auxiliary eigenvalue approach gives the values of the Green’s
function inside the entire waveguide for a given position of the source, which
can be both in the interior of the waveguide or on its boundaries. The
recursive Green’s function method, on the other hand, is particularly suitable
for solving the Green’s function connecting the two boundary layers, with the
source placed in either one of them [See Appendix C]. Extra steps are needed
to solve the Green’s function inside the waveguide, which is more intense
numerically.

##  III New insights

Having validated the proposed auxiliary eigenvalue approach to calculate the
Green’s function in various systems, we now turn to the new insight this
method provide, i.e., viewing the Green’s function as a defect state as
manifested by Equation (5).

We start by revisiting the 1D Hermitian system described by the scalar
Helmholtz equation in Figure 1. The eigenstate $\psi_{0}(x)$ in our auxiliary
eigenvalue equation (2) (and the Green’s function) for a given energy
$z=k^{2}$ then corresponds to the wave function given by the Schrödinger
equation inside a closed box with a $\delta$-function potential:

$\left[-\frac{1}{n^{2}}\partial_{x}^{2}+\lambda_{0}\delta(x-x^{\prime})\right]\psi_{0}(x)=z\psi_{0}(x),$
(22)

and the sign of $\lambda_{0}$ determines whether the potential is attractive
or repulsive, which can be found analytically using Equation (7).

This property of $\lambda_{0}$ is a crucial difference between our approach
and previous studies of defect states, where the defect strength (i.e.,
$\lambda_{0}$) is treated as a free parameter to find all possible $z$’s. In a
nutshell, a given defect strength leads to an infinite set of wave functions
with different energies $\\{z_{i}\\}$ [42], while for a given $z$ there is a
single wave function of Equation (22) [and more generally, Equation (5)],
i.e., the Green’s function, the value of which at the source gives
$\lambda_{0}^{-1}$. Such a one-to-one mapping from $z$ to $\lambda_{0}$ as an
eigenvalue problem has not been recognized previously to the best of our
knowledge, which has prevented the association of a Green’s function as a
defect state in the past.

To illustrate this point, we show in Figure 6(a) several energy eigenvalues
$z_{i}$’s of Equation (22) as a function of $\lambda_{0}$. It is important to
note that the “bandwidths” of the eigenstates do not overlap [shaded areas in
in Figure 6(a)], even when the range of $\lambda_{0}$ is extended from
$-\infty$ to $\infty$. This feature is determined by the uniqueness of the
Green’s function: for a given energy $z$ and source position $x^{\prime}$, the
Green’s function is uniquely determined by the boundary condition, and so is
its inverse height at $x^{\prime}$, which gives the only allowed value of
$\lambda_{0}$ for this pair of $z$ and $x^{\prime}$. Each of the “bandgaps” in
Figure 6(a) shrinks to a single point when $\lambda_{0}$ is extended from
$-\infty$ to $\infty$, and they correspond to the energy eigenvalues of the
two subsystems separated by the $\delta$-potential.

Figure 6: (a) Eigenvalues of the Schrödinger equation (22) as a function of
$\lambda_{0}$ for $n=3$ and $x^{\prime}=0.45L$. Only four of them are shown
near $z=900/L^{2}$. (b) Schematics showing the potential as a function of the
position. Shaded areas indicate the Dirichlet boundary condition and the
$\delta$-function potential with $\lambda_{0}=38.6/L$. Two unnormalized wave
functions $|\psi(x)|^{2}$ are also shown.

A particularly interesting type of defect states has a spatially localized
wave function. It takes place when the energy of the defect state is inside a
bandgap, causing the defect state to be largely decoupled from the rest of the
system. Here we show that the Green’s function in a 1D non-Hermitian lattice
can have the same property, which further stresses the new perspective our
method provides. This 1D lattice features a non-Hermitian flatband [15],
enabled by imposing non-Hermitian particle-hole (NHPH) symmetry [43, 44, 45,
46] via staggered gain and loss along a lattice of optical resonators. In the
tight-binding Hamiltonian description, it is given by

$H\phi_{n}=[z_{0}+(-1)^{n-1}i\gamma]\phi_{n}+g(\phi_{n+1}+\phi_{n-1}),$ (23)

where $n$ labels the lattice sites from 1 to $N$, $g\in\mathbb{R}$ is the
nearest neighbor coupling, and the alternate positive and negative imaginary
on-site detuning $\gamma$ models gain and loss. Below we shift the energy
reference such that the on-site energy $z_{0}$ of a resonator mode is zero.
NHPH symmetry then warrants a symmetric energy spectrum about the origin of
the complex plane, i.e., $z_{i}=-z_{j}^{*}$. At a critical value of the gain
and loss strength $\gamma=2g$, the real part of the energy bands collapse to
0, with all the eigenvalues $\\{z_{i}\\}$ completely imaginary.

Now let us consider the Green’s function of such a system, i.e.,

$(z-H)G=\delta_{nn^{\prime}}.$ (24)

In the auxiliary eigenvalue problem, it can be rewritten as

$(H+\lambda_{0}\delta_{nn^{\prime}})\psi_{n}=z\psi_{n},\quad
G=\frac{\psi_{n}}{\lambda_{0}\psi_{n^{\prime}}}.$ (25)

In other words, the Green’s function corresponds to a single defect state,
where the on-site energy at the position of the source $n^{\prime}$ is shifted
by $\lambda_{0}$.

We emphasize that this shift $\lambda_{0}$ is complex in general and hence
different from the Hermitian case shown in Figure 6 and other typical
considerations of a defect state, including that in Reference [15]. One
example is shown in Figure 7, where we place the source at the left edge
($n^{\prime}=1$) of a lattice with 60 lattice sites. $\lambda_{0}$ stays in
the lower complex plane when $z$ increases from $-g$ to $g$, and the stronger
defect strength (and lower LDOS) at $z=\pm g$ gives the Green’s function a
clearly localized profile, as reflected by the half participation ratio (HPR)
defined by $\sum_{n}\,(|\psi_{n}|^{2})^{2}/(2\sum_{n}\,|\psi_{n}|^{4})$ in the
units of the lattice constant [see Figure 7(b)]. We note that the HPR of an
exponentially localized intensity
$|\psi_{n}|^{2}=e^{-n/\xi}\,(n=1,2,\ldots,N)$ is approximately equal to its
localization length $\xi$ [47] when $\xi,N\gg 1$. The localization of the
Green’s function is weakened as the defect strength $\lambda_{0}$ becomes
smaller near $z=0$, with its HPR now greater than 15 lattice sites.

Figure 7: (a) Complex defect strength $\lambda_{0}$ when $z$ increases from
$-g$ to $g$. The lattice has 60 sites and $\gamma=2g$. (b) HPR as a function
of $z$, showing a strong localization (HPR $\sim 1$) at $|z|=t$ to a weaker
localization (HPR $>15$) at $z=0$.

Finally, we reveal the most intriguing connection between the Green’s function
and a defect state, i.e., in a 2D topological lattice with a chiral edge
state. This system breaks Lorentz reciprocity [48, 49], which can be achieved
in an electronic system by imposing a magnetic field. An analogy can be
introduced to photonic systems by imposing an artificial gauge field, achieved
experimentally by asymmetric couplings between two neighboring lattice sites
on a tight-binding lattice [16, 17, 18, 19]: while the couplings in both the
$x$ and $y$ directions still have the same amplitude, their phases are now
different. Here we consider a square lattice (Figure 8) with uniform vertical
coupling and horizontal asymmetric couplings of the same amplitude $g$,
realizing a Landau gauge with a $\pi/2$-flux through the smallest plaquette
[16, 17]. Its bottom right corner is pierced by the opposite flux, and an on-
site potential shift of $-2g$ is also introduced to decouple it from the rest
of the system.

Figure 8: Analogy between the Green’s function and a chiral edge state
circumventing a defect in a topological square lattice. (a) The Green’s
function calculated with $z=1.85g$ and the source at $x=a,y=10a$. Solid,
dashed, and dotted lines show the results of our auxiliary eigenvalue
approach, the recursive Green’s function, and the bilinear expansion,
respectively. (b,c) The values of the Green’s function along $y$ at $x=40a$
and $a$, respectively.

Due to its sublattice symmetry [50], the energy spectrum of the main region is
symmetric about $z=0$, set at the value of the identical on-site potential. It
has an edge band with CCW chiral edge states in $z/g\in[1.08,2.61]$, with
their CW counterparts in $z/g\in[-2.61,-1.08]$. We probe the response of the
system by placing a point source at $x=a,y=10a$, where $a=1$ is the lattice
constant. As Figure 8(a) shows, the Green’s function with $z$ at the middle of
the CCW chiral edge band displays the same characteristic as the eigenstates
of $\mathcal{L}$, i.e., localized near the edge of the topological region. Our
auxiliary eigenvalue approach gives identical result as the recursive Green’s
function method [solid and dashed lines in Figs. 8(b) and (c)], and the
Green’s function displays two noticeable features when compared with the
bilinear expansion: while only a small number ($\sim$20) of edge states are
needed to capture the Green’s function on the opposite (right) edge of the
system (dotted line), they cannot describe the behavior of the Green’s
function near the source, where it displays a local minimum.

One way to understand this behavior is offered by our auxiliary eigenvalue
approach, where the Green’s function is effectively a defect state perturbed
by the $\delta$-function potential [see Equation (5)]. In this point of view,
the local minimum of the Green’s function is the consequence of the
topological protection of the chiral edge states: it circumvents this defect
as can be clearly seen in Figure 8(a). This change, of course, creates a local
perturbation so strong that can only be captured by a large number of the
unperturbed eigenstates of $\mathcal{L}$ in the bilinear expansion, including
even the chiral edge states in the opposite direction that are much different
in energy (not shown). The far-field, on the other hand, is not affected by
this local perturbation, as we have seen in Figure 8(b). We note that the
Green’s function does not display a local minimum at all energies: if $z$ is
almost on resonance with an unperturbed eigenstate in the edge band, the
Green’s function is given essentially by this single chiral edge mode, which
goes through the defect created by the $\delta$-function potential with little
scattering.

##  IV Discussion and conclusions

In summary, we have introduced in this work a new perspective of the time-
independent Green’s function, which is obtained as a single eigenstate of an
auxiliary eigenvalue formulation that embodies a defect state created by the
$\delta$-function potential. The height of the $\delta$-function potential is
determined by the inverse value of the Green’s function at the source
position, which is given directly in the form of a generalized eigenvalue
problem given by Equation (2). It is the only well behaved and finite
eigenvalue, easy to identify numerically. Therefore, our approach differs both
conceptually and computationally from previous investigations of eigenstates
of a $\delta$-function potential, which were not related to the Green’s
function.

The uniqueness of the eigenstate that gives the Green’s function in our
approach should be distinguished from the single bound state in an attractive
1D $\delta$-function potential [51]. In our case the $\delta$-function
potential can be repulsive or attractive in a Hermitian system [see Figure
6(a)], and it becomes complex in general in non-Hermitian systems.
Furthermore, the eigenstate in our approach does not depend on the original
potential included in the operator $\cal L$, while the aforementioned bound
state assumes a vacuum background. Finally, our eigenstate exists in higher
dimensions as well, with clearly different boundary conditions and spatial
profiles from a bound state [see Figure 8(a), for example].

We have also verified the Green’s function obtained in our method by comparing
both to analytical results when available and to two other frequently used
numerical methods, i.e., the recursive Green’s function method and the
bilinear expansion in the basis of the system’s eigenstates. At an EP, where a
perturbative treatment of the bilinear expansion becomes necessary, our method
is still robust as seen from the examples in two $\mathcal{PT}$-systems. Our
method also gives identical result to the recursive Green’s function method,
implemented by finite difference on the same tight-binding lattice.

Our defect-state approach can also be applied to more numerically demanding
cases, e.g., in the study of diffusive transport and wave localization [24,
25]. We have investigated disordered quasi-1D waveguides over 100 wavelengths
long and with over 60 transverse channels (not shown). In these cases, the
memory storage for the recursive Green’s function may become an issue, when
the values of the Green’s function in the interior of the waveguides are also
computed for light-matter interaction or laser emissions. This is because the
Green’s matrix obtained from the recursive procedure mingles the values of the
Green’s function generated by sources placed across the entire system. The
defect-state approach excels in this regard, because each source is treated
independently. We should also mention that this advantage becomes less
noticeable or even a disadvantage if the Green’s function needs to be
evaluated with the source at many different locations.

More importantly, our approach offers a previously unexplored physical insight
that both the recursive Green’s function method and the bilinear expansion
lack, i.e., the linkage between the Green’s function and a defect state. We
have exemplified an intriguing manifestation of this linkage using a
topological chiral edge state, where the local minimum of the Green’s function
is analogous to a chiral edge state circumventing a boundary defect.
Therefore, even though our discussions have focused on the Helmholtz equation
for scalar optical waves, we expect this new perspective and the physical
insight it offers to have important applications in other related fields as
well, including condensed matter systems, acoustics, electronic circuits and
so on.

Acknowledgements

We thank Douglas Stone, Yidong Chong and Azriel Genack for helpful
discussions. This project is supported by the NSF under Grant No. PHY-1847240.

##  Appendix A Reciprocity in the auxiliary eigenvalue approach

We derive some relations between the Green’s functions
$G(\bm{r},\bm{r^{\prime}})$ and $G(\bm{r^{\prime}},\bm{r})$, based on the
symmetries of $\mathcal{L}$. We start from a matrix description in which the
Green’s function operator $G(z)$ is defined from
$G(\bm{r},\bm{r^{\prime}};z)=\langle\bm{r}|G(z)|\bm{r^{\prime}}\rangle$. We
use the same symbol $\mathcal{L}$ for the matrix representation of the system.
The Green’s function satisfies the following equation:

$\displaystyle[z\bm{1}-\mathcal{L}]G(z)=\bm{1},\hskip
14.22636ptG(z)[z\bm{1}-\mathcal{L}]=\bm{1},$ (26)

where $\bm{1}$ is the identity matrix. Here, we identify the Green’s function
as being both the right and left matrix inverses of $[z\bm{1}-\mathcal{L}]$.
Considering that in non-Hermitian physics, there is a distinction between left
and right sets of eigenvectors [52], one may attempt to distinguish between
the left and right Green’s functions. However, it can be shown that they are
identical, since $[z\bm{1}-\mathcal{L}]$ is a matrix of full rank. It is
possible to obtain relations between $G(z)$ and $[G(z)]^{T}$, or
$[G(z)]^{\dagger}$, when performing matrix transposition or hermitian
conjugation on Equation (26).

We start by analyzing the case $\mathcal{L}^{T}=\mathcal{L}$: matrix
transposition of Equation (26) leads to the relation $G(z)^{T}=G(z)$, or more
explicitly $G(\bm{r^{\prime}},\bm{r};z)=G(\bm{r},\bm{r^{\prime}};z)$, which is
the usual Lorentz reciprocity condition. This relation holds when
$\mathcal{L}$ is a symmetric matrix, describing a system with symmetric
couplings and hence without an effective magnetic field. It is valid for non-
Hermitian systems with this property as well, such as the case with a non-
Hermitian flatband discussed in Figure 7 of the main text.

In a general Hermitian system, $\mathcal{L}^{\dagger}=\mathcal{L}$, and
uniqueness of the matrix inverse demands that $[G(z)]^{T}=[G(z^{*})]^{*}$, or
$G(\bm{r^{\prime}},\bm{r};z)=[G(\bm{r},\bm{r^{\prime}};z^{*})]^{*}$. It is
customary to state this relation in terms of Green’s functions $G^{(\pm)}$
with outgoing and incoming boundary conditions, respectively. This
nomenclature is specially useful if $z=k+i\eta$, with $\eta>0$ small compared
to $k$:

$G^{(\pm)}(\bm{r},\bm{r^{\prime}};k)\equiv\lim\limits_{\eta\to
0^{+}}G(\bm{r},\bm{r^{\prime}};k\pm i\eta).$ (27)

Therefore, we can write the reciprocity condition as
$G^{(\pm)}(\bm{r^{\prime}},\bm{r};k)=[G^{(\mp)}(\bm{r},\bm{r^{\prime}};k)]^{*}$.
The solutions $G^{(\pm)}$ are oftentimes referred to as the retarded and
advanced Green’s functions, inspired from the association to sources
generating outgoing or incoming waves. Systems that are non-Hermitian and do
not have symmetric couplings lead to more complicated reciprocity relations,
involving the Green’s functions of two different systems. This is expected,
since the scattering properties derived from the Green’s functions of systems
related by a time-reversal transformation are nontrivial and yet to be
understood. Reference [53] provides interesting scattering relations between
two “time-reversal partners” in 1D and quasi 1D systems.

Our Green’s function approach for time-reversal or non-Hermitian symmetric
systems lead to a relation between the single eigenstates:

$\frac{\psi_{0}(\bm{r},\bm{r^{\prime}})}{\lambda_{0}(\bm{r^{\prime}})\psi_{0}(\bm{r^{\prime}},\bm{r^{\prime}})}=\frac{\psi_{0}(\bm{r^{\prime}},\bm{r})}{\lambda_{0}(\bm{r})\psi_{0}(\bm{r},\bm{r})},$
(28)

or

$\psi_{0}(\bm{r^{\prime}},\bm{r})=\psi_{0}(\bm{r},\bm{r^{\prime}})\left[\frac{\lambda_{0}(\bm{r})\psi_{0}(\bm{r},\bm{r})}{\lambda_{0}(\bm{r^{\prime}})\psi_{0}(\bm{r^{\prime}},\bm{r^{\prime}})}\right],$
(29)

whereas for the general non-Hermitian problem, we find

$\bar{\psi}_{0}(\bm{r^{\prime}},\bm{r})=[\psi_{0}(\bm{r},\bm{r^{\prime}})]^{*}\left[\frac{\bar{\lambda}_{0}(\bm{r})\bar{\psi}_{0}(\bm{r},\bm{r})}{\lambda_{0}(\bm{r^{\prime}})\psi_{0}(\bm{r^{\prime}},\bm{r^{\prime}})}\right],$
(30)

where we used the notation
$\bar{\psi}_{0}(\bm{r},\bm{r^{\prime}},z)\equiv\psi_{0}(\bm{r},\bm{r^{\prime}},z^{*})$.
The same result is obtained if we identify $\bar{\psi}_{0}$ with the advanced
solution and $\psi_{0}$ with the retarded solution, or vice-versa.

##  Appendix B Perturbation theory at an EP

Below we review a generic procedure to approach effective
$\mathcal{PT}$-symmetric systems with the outgoing boundary condition based on
the coupled-mode theory. In particular, we describe a photonic molecule that
consists on identical components $\mathbb{L}$ and $\mathbb{R}$ that support
localized states $\psi_{l}(x)$ and $\psi_{r}(x)$ in each component. The system
can be, for instance, two identical whispering gallery resonators that are
brought to proximity, such that their modes overlap and produce a pair of
bonding and antibonding states. Alternatively, one can think of a pair of
half-wavelength cavities coupled by a DBR, which is the system discussed in
Section 2.2 of the main text. Note that due to the outgoing boundary
condition, such systems are not truly $\cal PT$-symmetric, because time
reversal changes the outgoing boundary condition to the incoming boundary
condition. However, if the states are highly localized as in our discussion
below, the minute flux in or out of the system can be omitted in the
discussion of the coupled mode theory.

The left component satisfies the Helmholtz equation:

$\displaystyle\left[\partial_{x}^{2}+\epsilon_{l}(x)k_{0}^{2}\right]\psi_{l}(x)$
$\displaystyle=0,\quad(x\in\mathbb{L})$ (31)
$\displaystyle\left[\partial_{x}^{2}+k^{2}\right]\psi_{l}(x)$
$\displaystyle=0.\quad(x\in\text{elsewhere})$ (32)

Here $k$ is the real-valued external frequency at which the Green’s function
is going to be evaluated, and $k_{0}\sim k$ is one complex-valued CF
frequency. The right component is defined similarly, and the localized states
are normalized by
$\int_{x\in\mathbb{L},\mathbb{R}}\epsilon_{l,r}\psi_{l,r}^{2}dx=1$,
respectively. The composite system is defined by
$\epsilon_{0}(x)=\epsilon_{l}(x)+\epsilon_{r}(x)$ and satisfies

$\displaystyle\left[\partial_{x}^{2}+\epsilon_{0}(x)\tilde{k}^{2}\right]\tilde{\psi}(x)$
$\displaystyle=0,\quad(x\in\mathbb{L}\cup\mathbb{R})$ (33)
$\displaystyle\left[\partial_{x}^{2}+k^{2}\right]\tilde{\psi}(x)$
$\displaystyle=0,\quad(x\in\text{elsewhere})$ (34)

where $\tilde{k}$ is the CF frequency of a “supermode”

$\tilde{\psi}(x)=a_{l}\psi_{l}(x)+a_{r}\psi_{r}(x)$ (35)

to be determined. We note that $\tilde{\psi}(x)$ automatically satisfies the
outgoing boundary condition at the two ends of the composite system due to the
same property of both $\psi_{l,r}(x)$.

We make the following assumption about the system: the spatial overlapping of
the localized states in one cavity is very weak, i.e.,
$\int_{x\in\mathbb{L},\mathbb{R}}\epsilon_{l,r}\psi_{l}\psi_{r}dx\equiv
g\;(|g|\ll 1)$. In addition, we note that the frequency shift due to
$|\int_{x\in\mathbb{R}}\epsilon_{r}\psi_{l}^{2}dx|=|\int_{x\in\mathbb{L}}\epsilon_{l}\psi_{r}^{2}dx|$
is even smaller, which is then neglected in our analysis below.

To derive the coupled-mode theory, we first insert Equation (35) into Eqs.
(33), (34) and simplify the result using Eqs. (31), (32) for $\psi_{l}(x)$ and
their counterparts for $\psi_{r}(x)$. Then, we multiply the resulting equation
by either $\psi_{l}(x)$ or $\psi_{r}(x)$ and perform an integration in the
left and right component of the system respectively, obtaining in this way two
independent equations:

$\begin{bmatrix}\tilde{k}^{2}-k_{0}^{2}&g\tilde{k}^{2}-tk^{2}\\\
g\tilde{k}^{2}-tk^{2}&\tilde{k}^{2}-k_{0}^{2}\end{bmatrix}\begin{bmatrix}a_{l}\\\
a_{r}\end{bmatrix}=0.$ (36)

It is a system of equations for the unknown amplitudes $a_{l,r}$ and frequency
$\tilde{k}$, where $t$ is defined by
$\int_{x\in\mathbb{L}}\psi_{r}^{2}dx=\int_{x\in\mathbb{R}}\psi_{l}^{2}dx$.
Since the matrix operator is symmetric and has identical diagonal elements, a
zero-determinant condition implies that $a_{r}=\pm a_{l}$, giving rise to the
symmetric and anti-symmetric modes $\tilde{\psi}_{S}(x)$ and
$\tilde{\psi}_{A}(x)$ respectively, with frequencies:

$\displaystyle\tilde{k}_{S,A}^{2}=\frac{k_{0}^{2}\pm tk^{2}}{1\pm g}.$ (37)

With the aforementioned approximations $|g|\ll 1$ and $k_{0}\sim k$, these
frequencies reduce to

$\displaystyle\tilde{k}_{S}^{2}\approx(1+t-g)k_{0}^{2},\hskip
14.22636pt\tilde{k}_{A}^{2}\approx\left(1-t+g\right)k_{0}^{2},$ (38)

and their splitting equal to
$\tilde{k}_{S}^{2}-\tilde{k}_{A}^{2}=2(t-g)k_{0}^{2}\equiv 2\Delta k_{0}^{2}$.

Now, we introduce a $\mathcal{PT}$-symmetric perturbation $i\epsilon_{1}(x)$,
which is an odd function of $x$. The effective $\cal PT$-symmetric system
satisfies

$[\partial_{x}^{2}+(\epsilon_{0}(x)+i\epsilon_{1}(x))q^{2}]\psi(x)=0,$ (39)

where $q$ is the perturbed CF frequency to be determined together with the
amplitudes $A$ and $S$ in the expansion of the perturbed eigenstate
$\psi(x)=A\tilde{\psi}_{A}(x)+S\tilde{\psi}_{S}(x)$. Following an analogous
procedure to the one described above, we obtain again a system of equations
that can be represented by the following matrix equation for $A$, $S$:

$\begin{bmatrix}q^{2}-\tilde{k}_{S}^{2}&iCq^{2}\\\
iCq^{2}&q^{2}-\tilde{k}_{A}^{2}\end{bmatrix}\begin{bmatrix}S\\\
A\end{bmatrix}=0,$ (40)

where $C\equiv\int\epsilon_{1}(x)\tilde{\psi}_{A}(x)\tilde{\psi}_{S}(x)dx$ is
integrated over the whole system. We note that this equation and the following
perturbative results apply when $\epsilon_{0}(x)$ is complex, i.e., the system
does not need to have balanced gain and loss to be effectively $\cal{PT}$
symmetric [35].

To rewrite it in a more familiar form, we note that since the $\cal
PT$-symmetric perturbation is assumed to be weak, we approximate $Cq^{2}$ by
$Ck_{0}^{2}$ that leads to

$\tilde{H}\begin{bmatrix}S\\\
A\end{bmatrix}\equiv\begin{bmatrix}k_{S}^{2}&-iCk_{0}^{2}\\\
-iCk_{0}^{2}&k_{A}^{2}\end{bmatrix}\begin{bmatrix}S\\\
A\end{bmatrix}=q^{2}\begin{bmatrix}S\\\ A\end{bmatrix}.$ (41)

The perturbed frequency $q$ is then found to be:

$q^{2}_{\pm}\approx k_{0}^{2}\left(1\pm\sqrt{\Delta^{2}-C^{2}}\right),$ (42)

where we have used the forms of $\tilde{k}_{S,A}^{2}$ from Equation (38). The
EP then emerges when the $C^{2}=\Delta^{2}$, where $q_{\pm}\equiv
k_{\text{EP}}\approx k_{0}$. In the main text we have introduced
$\varsigma\equiv\sqrt{1-\beta^{2}}$ and $\beta\equiv C/\Delta$.

This system of equations can be transformed to the basis of left- and right-
components, in terms of the amplitudes $a_{l}$, $a_{r}$:

$H\begin{bmatrix}a_{l}\\\ a_{r}\end{bmatrix}\equiv
k_{0}^{2}\begin{bmatrix}1-iC&\Delta\\\
\Delta&1+iC\end{bmatrix}\begin{bmatrix}a_{l}\\\
a_{r}\end{bmatrix}=q^{2}\begin{bmatrix}a_{l}\\\ a_{r}\end{bmatrix}.$ (43)

It has the most familiar form of a $\mathcal{PT}$-symmetric Hamiltonian. This
is the model that we have used to describe the photonic molecule in Section
2.2.

The perturbative result of the Green’s function given by Equation (16) in the
main text is derived from the following expression:

$\displaystyle
G(x,x^{\prime};k)=\frac{\epsilon_{0}(x^{\prime})}{(k^{2}-k_{0}^{2})^{2}-\Delta^{2}\varsigma^{2}}\\{(k^{2}$
$\displaystyle-
k_{0}^{2})[\psi_{l}(x^{\prime})\psi_{l}(x)+\psi_{r}(x^{\prime})\psi_{r}(x)]\,$
$\displaystyle-\Delta
k_{0}^{2}[\,i\beta[\psi_{l}(x^{\prime})\psi_{l}(x)-\psi_{r}(x^{\prime})\psi_{r}(x)]-[\psi_{l}(x^{\prime})\psi_{r}(x)+\psi_{r}(x^{\prime})\psi_{l}(x)]\,]\\}$
(44)

by dropping terms proportional to $\varsigma^{2}$, including that in
$\beta=\sqrt{1-\varsigma^{2}}\approx 1$. We note that Equation (44) includes
terms to all orders of $\varsigma$. Here it is unnecessary to invoke the
Jordan vector $J(x)$ that completes the Hilbert space at the EP [23], but to
compare with the result given in Reference [10], we note that it can be chosen
as

$J(x)=\frac{i\psi_{l}(x)+\psi_{r}(x)}{2\sqrt{\Delta}k_{0}}e^{-i\pi/4},$ (45)

which satisfies

$[H-k_{\text{EP}}^{2}\mathbf{1}]J(x)=\psi_{\text{EP}}(x)$ (46)

with the coalesced eigenstate

$\psi_{\text{EP}}(x)=2\sqrt{\Delta}k_{0}e^{-i\pi/4}[\psi_{l}(x)+i\psi_{r}(x)]$
(47)

at the EP. Here $\mathbf{1}$ is the identity matrix, and $\psi_{\text{EP}}(x)$
is normalized differently from the main text to satisfy
$(\psi_{\text{EP}},J)=1$ [see the definition of the inner product given by
Equation (14) in the main text]. We further note that the $\cal PT$-symmetric
Hamiltonian given by Equation (43) above is symmetric, and hence it is
unnecessary to distinguish left and right eigenstates (and Jordan vectors)
because they are identical. Using $J(x)$ and $\psi_{\text{EP}}(x)$ specified
above, we find that the Green’s function derived in Reference [10] at the EP
agrees with our result given by Equation (16) in the main text, once the
additional factor $\epsilon_{0}(x^{\prime})$ is accounted for that originates
from rewriting the Helmholtz equation as Equation (6) in the main text.

##  Appendix C The recursive Green’s function

Below we briefly review the recursive Green’s function method and our
implementation of it using finite difference. Assuming that the Green’s
function defined by Equation (1) has been obtained at every grid point for a
quasi-1D waveguide with $N-1$ columns and $N_{y}$ rows, we can write down a
Green’s matrix $\bm{G}^{(N-1)}$ that satisfies

$\left[z\bm{1}^{(N-1)}-\bm{L}^{(N-1)}\right]\bm{G}^{(N-1)}=\bm{1}^{(N-1)}.$
(48)

Its individual element $\bm{G}^{(N-1)}_{\underline{i,j}}$ gives the value of
the Green’s function at point $i$ when the source is placed at point $j$. We
reserve the notation $\bm{G}^{(N-1)}_{i,j}$ without the underscore below the
subscripts for the $N_{y}\times N_{y}$ block of the Green’s matrix evaluated
in the $i$th column with the source in the $j$th column, and we apply the same
notation to other matrices as well. The $(N-1)N_{y}$ grid points are labeled
from the first to the last point in the first column, continued in the same
way onto the second and remaining columns. $\bm{1}^{(N-1)}$ is the identity
matrix with $(N-1)N_{y}$ rows and columns, and $\bm{L}^{(N-1)}$ is the finite
difference implementation of the operator $\mathcal{L}$ in Equation (1), again
with $(N-1)N_{y}$ rows and columns.

To calculate the Green’s matrix with an additional column added from the
right, we separate $\bm{L}^{(N)}$ into $\bm{L}^{0}+\bm{V}$ and construct an
ancillary matrix $\bm{G}^{0}$, where

$\bm{L}^{0}\equiv\begin{pmatrix}\bm{L}^{(N-1)}&\\\
&\bm{L}^{(N)}_{N,N}\end{pmatrix},\quad\bm{G}^{0}\equiv\begin{pmatrix}\bm{G}^{(N-1)}&\\\
&\bm{G}^{0}_{N,N}\end{pmatrix}.$ (49)

$\bm{V}$ only contains the couplings between the $(N-1)$th and $N$th columns
in $\bm{L}^{(N)}$, and

$\left[z\bm{1}^{(1)}-\bm{L}^{(N)}_{N,N}\right]\bm{G}^{0}_{N,N}=\bm{1}^{(1)}$
(50)

defines the Green’s matrix $\bm{G}^{0}_{N,N}$ of the isolated $N$th column. It
is straightforward to show that

$\left[z\bm{1}^{(N)}-\bm{L}^{0}\right]\bm{G}^{0}=\bm{1}^{(N)}$ (51)

and we derive the well known Dyson’s equation

$\bm{G}^{(N)}=\bm{G}^{0}+\bm{G}^{(N)}\bm{V}\bm{G}^{0},$ (52)

where $\bm{G}^{(N)}$ is defined by

$\left[z\bm{1}^{(N)}-\bm{L}^{(N)}\right]\bm{G}^{(N)}=\bm{1}^{(N)}$ (53)

similar to Equation (48). Equation (52) can be verified directly by
multiplying $[z\bm{1}^{(N)}-\bm{L}^{(N)}]$ from the left.

Below we drop the superscript of $\bm{G}^{(N)}$. Equation (52) differs from
its more familiar form in mesoscopic physics [41], i.e.,

$\bm{\tilde{G}}=\bm{\tilde{G}}^{0}+\bm{\tilde{G}}^{0}\bm{V}\bm{\tilde{G}},$
(54)

because there the Green’s matrix is defined by
$\bm{\tilde{G}}[z\bm{1}^{(N)}-\bm{L}^{(N)}]=\bm{1}^{(N)}$ instead. We have
already proven the equivalence between the left and right matrix inverses of
$[z\bm{1}-\mathcal{L}]$ in Appendix A, thus $\tilde{\bm{G}}$ and $\bm{G}$ are
identical.

Using Equation (52), the recursive relations between the blocks of the Green’s
matrices $\bm{G}_{i,j}$ connecting the two boundary layers (i.e., the 1st and
$N$th columns) can be derived:

$\displaystyle\bm{G}_{NN}=\;A^{-1},$ (55)
$\displaystyle\bm{G}_{1N}\;=\;\bm{G}^{0}_{1,N-1}\bm{V}_{N-1,N}A^{-1},$ (56)
$\displaystyle\bm{G}_{N1}\;=\;\bm{G}^{0}_{N,N}\bm{V}_{N,N-1}\bm{G}^{0}_{N-1,1},$
(57)
$\displaystyle\bm{G}_{11}\;\;=\;\bm{G}^{0}_{1,1}+\bm{G}_{1,N}\bm{V}_{N,N-1}\bm{G}^{0}_{N-1,1}.$
(58)

A single matrix inversion for

$A=z\bm{1}^{(1)}-\bm{L}^{(N)}_{N,N}-V_{N,N-1}\bm{G}^{0}_{N-1,N-1}V_{N-1,N}$
(59)

is needed when adding one layer from the right, and its last term is often
referred to as the self-energy [40]. Note that besides the operator $L^{(N)}$,
all the pieces of information needed for this recursive scheme are the same
four matrix for the system with $(N-1)$ layers. As a result, the recursive
Green’s function method is highly efficient to treat transport problems, which
does not require the knowledge of the Green’s function inside the system. If
we need the full Green’s matrix to investigate, for example, light-matter
interaction in the bulk of the system, more steps are required in the
recursive Green’s function, which scale as $N^{2}$ and are more intense
numerically.

## References

  * Green [1871] G. Green,  _Mathematical papers of the late George Green_ (Macmillan, London, 1871).
  * Stakgold and Holst [2011] I. Stakgold and M. Holst,  _Green’s functions and boundary value problems_, 3rd ed., Pure and Applied Mathematics (John Wiley & Sons, Hoboken, 2011).
  * Melnikov and Borodin [2017] Y. A. Melnikov and V. N. Borodin,  _Green’s Functions_ , Developments in Mathematics No. 48 (Springer, Cham, 2017).
  * Morse and Feshbach [1953] P. M. Morse and H. Feshbach,  _Methods of Theoretical Physics_ , Vol. 1 (McGraw-Hill, New York, 1953).
  * Schwinger [1998] J. Schwinger,  _Particles, Sources, and Fields_ (Perseus, Reading, Mass, 1998).
  * Rickayzen [1980] G. Rickayzen,  _Green’s functions and Condensed Matter_ (Academic Press, New York, 1980).
  * Economou [2007] E. Economou,  _Green’s functions in Quantum Physics_, 3rd ed. (Springer, 2007).
  * Davy  _et al._ [2015] M. Davy, Z. Shi, J. Wang, X. Cheng, and A. Z. Genack, Transmission Eigenchannels and the Densities of States of Random Media, Phys. Rev. Lett. 114, 033901 (2015).
  * Lin  _et al._ [2016] Z. Lin, A. Pick, M. Lončar, and A. W. Rodriguez, Enhanced Spontaneous Emission at Third-Order Dirac Exceptional Points in Inverse-Designed Photonic Crystals, Phys. Rev. Lett. 117, 107402 (2016).
  * Pick  _et al._ [2017] A. Pick, B. Zhen, O. D. Miller, C. W. Hsu, F. Hernandez, A. W. Rodriguez, M. Soljačić, and S. G. Johnson, General theory of spontaneous emission near exceptional points, Opt. Express 25, 12325 (2017).
  * Newton [1982] R. G. Newton,  _Scattering Theory of Waves and Particles_, 2nd ed. (Springer, Berlin, 1982).
  * Türeci  _et al._ [2006] H. E. Türeci, A. D. Stone, and B. Collier, Self-consistent multimode lasing theory for complex or random lasing media, Phys. Rev. A 74, 043822 (2006).
  * Türeci  _et al._ [2008] H. E. Türeci, L. Ge, S. Rotter, and A. D. Stone, Strong interactions in multimode random lasers, Science 320, 643 (2008).
  * Ge  _et al._ [2010] L. Ge, Y. D. Chong, and A. D. Stone, Steady-state ab initio laser theory: Generalizations and analytic results, Phys. Rev. A 82, 063824 (2010).
  * Qi  _et al._ [2018] B. Qi, L. Zhang, and L. Ge, Defect States Emerging from a Non-Hermitian Flatband of Photonic Zero Modes, Phys. Rev. Lett. 120, 093901 (2018).
  * Hafezi  _et al._ [2011] M. Hafezi, E. A. Demler, M. D. Lukin, and J. M. Taylor, Robust optical delay lines with topological protection, Nat. Phys. 7, 907 (2011).
  * Hafezi  _et al._ [2013] M. Hafezi, S. Mittal, J. Fan, A. Migdall, and J. M. Taylor, Imaging topological edge states in silicon photonics, Nat. Photonics 7, 1001 (2013).
  * Bandres  _et al._ [2018] M. A. Bandres, S. Wittek, G. Harari, M. Parto, J. Ren, M. Segev, D. N. Christodoulides, and M. Khajavikhan, Topological insulator laser: Experiments, Science 359, eaar4005 (2018).
  * Zhao  _et al._ [2019] H. Zhao, X. Qiao, T. Wu, B. Midya, S. Longhi, and L. Feng, Non-Hermitian topological light steering, Science 365, 1163 (2019).
  * Barik  _et al._ [2018] S. Barik, A. Karasahin, C. Flower, T. Cai, H. Miyake, W. DeGottardi, M. Hafezi, and E. Waks, A topological quantum optics interface, Science 359, 666 (2018).
  * Bender and Boettcher [1998] C. M. Bender and S. Boettcher, Real spectra in non-Hermitian Hamiltonians having PT symmetry, Phys. Rev. Lett. 80, 5243 (1998).
  * Feng  _et al._ [2017] L. Feng, R. El-Ganainy, and L. Ge, Non-Hermitian photonics based on parity–time symmetry, Nat. Photonics 11, 752 (2017).
  * Chen  _et al._ [2020] H. Z. Chen, T. Liu, H. Y. Luan, R. J. Liu, X. Y. Wang, X. F. Zhu, Y. B. Li, Z. M. Gu, S. J. Liang, H. Gao, L. Lu, L. Ge, S. Zhang, J. Zhu, and R. M. Ma, Revealing the missing dimension at an exceptional point, Nat. Phys. 16, 571 (2020).
  * Beenakker [1997] C. W. J. Beenakker, Random-matrix theory of quantum transport, Rev. Mod. Phys. 69, 731 (1997).
  * Imry [2008] Y. Imry,  _Introduction to Mesoscopic Physics_ , 2nd ed. (Oxford University Press, Oxford, 2008).
  * Mello and Kumar [2004] P. A. Mello and N. Kumar,  _Quantum Transport in Mesoscopic Systems_ (Oxford University Press, 2004).
  * Ge [2017a] L. Ge, Constructing the scattering matrix for optical microcavities as a nonlocal boundary value problem, Photonics Res. 5, B20 (2017a).
  * Liertzer  _et al._ [2012] M. Liertzer, L. Ge, A. Cerjan, A. D. Stone, H. E. Türeci, and S. Rotter, Pump-induced exceptional points in lasers, Phys. Rev. Lett. 108, 173901 (2012).
  * Heiss [2004] W. D. Heiss, Exceptional points of non-Hermitian operators, J. Phys. A 37, 2455 (2004).
  * Berry [2004] M. V. Berry, Physics of nonhermitian degeneracies, Czech. J. Phys. 54, 1039 (2004).
  * Ge and El-Ganainy [2016] L. Ge and R. El-Ganainy, Nonlinear modal interactions in parity-time (PT) symmetric lasers, Sci. Rep. 6, 24889 (2016).
  * Kapur and Peierls [1937] P. L. Kapur and R. Peierls, The dispersion formula for nuclear reactions, P. Roy Soc. A-Math. Phy. 166, 277 (1937).
  * Goldberger and Watson [1964] M. L. Goldberger and K. M. Watson, Lifetime and decay of unstable particles in S-matrix theory, Phys. Rev. 136, B1472 (1964).
  * Heiss [2015] W. D. Heiss, Green’s Functions at Exceptional Points, Int. J. Theor. Phys. 54, 3954 (2015).
  * Guo  _et al._ [2009] A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti, M. Volatier-Ravat, V. Aimez, G. A. Siviloglou, and D. N. Christodoulides, Observation of PT-symmetry breaking in complex optical potentials, Phys. Rev. Lett. 103, 93902 (2009).
  * Miao  _et al._ [2016] P. Miao, Z. Zhang, J. Sun, W. Walasik, S. Longhi, N. M. Litchinitser, and L. Feng, Orbital angular momentum microlaser, Science 353, 464 (2016).
  * Fisher and Lee [1981] D. S. Fisher and P. A. Lee, Relation between conductivity and transmission matrix, Phys. Rev. B 23, 6851 (1981).
  * Stone and Szafer [1988] A. D. Stone and A. Szafer, What is measured when you measure a resistance? The landauer formula revisited, IBM J. Res. Dev. 32, 384 (1988).
  * Ge  _et al._ [2015] L. Ge, K. G. Makris, D. N. Christodoulides, and L. Feng, Scattering in PT- and RT-symmetric multimode waveguides: Generalized conservation laws and spontaneous symmetry breaking beyond one dimension, Phys. Rev. A 92, 062135 (2015).
  * Thouless and Kirkpatrick [1981] D. J. Thouless and S. Kirkpatrick, Conductivity of the disordered linear chain, J. Phys. C 14, 235 (1981).
  * Sols  _et al._ [1989] F. Sols, M. Macucci, U. Ravaioli, and K. Hess, Theory for a quantum modulated transistor, J. Appl. Phys. 66, 3892 (1989).
  * Patil [2006] S. H. Patil, Harmonic oscillator with a $\updelta$-function potential, Eur. J. Phys. 27, 899 (2006).
  * Malzard  _et al._ [2015] S. Malzard, C. Poli, and H. Schomerus, Topologically Protected Defect States in Open Photonic Systems with Non-Hermitian Charge-Conjugation and Parity-Time Symmetry, Phys. Rev. Lett. 115, 200402 (2015).
  * Ge [2017b] L. Ge, Symmetry-protected zero-mode laser with a tunable spatial profile, Phys. Rev. A 95, 023812 (2017b).
  * Rivero and Ge [2021] J. D. H. Rivero and L. Ge, Chiral symmetry in non-hermitian systems: Product rule and clifford algebra, Phys. Rev. B 103, 014111 (2021).
  * Kawabata  _et al._ [2019] K. Kawabata, S. Higashikawa, Z. Gong, Y. Ashida, and M. Ueda, Topological unification of time-reversal and particle-hole symmetries in non-Hermitian physics, Nature Commun. 10, 297 (2019).
  * Baboux  _et al._ [2016] F. Baboux, L. Ge, T. Jacqmin, M. Biondi, E. Galopin, A. Lemaitre, L. Le Gratiet, I. Sagnes, S. Schmidt, H. E. Tureci, A. Amo, and J. Bloch, Bosonic Condensation and Disorder-Induced Localization in a Flat Band, Phys. Rev. Lett. 116, 066402 (2016).
  * Landau and Lifshitz [1960] L. D. Landau and E. M. Lifshitz,  _Electrodynamics of Continuous Media_ (Pergamon Press, Oxford, 1960).
  * Ge and Feng [2016] L. Ge and L. Feng, Optical-reciprocity-induced symmetry in photonic heterostructures and its manifestation in scattering PT-symmetry breaking, Phys. Rev. A 94, 043836 (2016).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * Gottfried and Yan [2003] K. Gottfried and T.-M. Yan,  _Quantum Mechanics: Fundamentals_ , 2nd ed. (Springer, New York, 2003).
  * Rivero and Ge [2020] J. D. H. Rivero and L. Ge, Pseudochirality: A Manifestation of Noether’s Theorem in Non-Hermitian Systems, Phys. Rev. Lett. 125, 083902 (2020).
  * Rivero and Ge [2019] J. D. H. Rivero and L. Ge, Time-reversal-invariant scaling of light propagation in one-dimensional non-hermitian systems, Phys. Rev. A 100, 023819 (2019).


# Lifelong Matrix Completion with Sparsity-Number

Ilqar Ramazanli
Carnegie Mellon University
Pittsburgh, USA
<EMAIL_ADDRESS>
###### Abstract

Matrix completion problem has been previously studied under various adaptive
and passive settings. Previously, researchers have proposed passive, two-phase
and single-phase algorithms using coherence parameter, and multi phase
algorithm using sparsity-number. It has been shown that the method using
sparsity-number reaching to theoretical lower bounds in many conditions.
However, the aforementioned method is running in many phases through the
matrix completion process, therefore it makes much more informative decision
at each stage. Hence, it is natural that the method outperforms previous
algorithms. In this paper, we are using the idea of sparsity-number and
propose and single-phase column space recovery algorithm which can be extended
to two-phase exact matrix completion algorithm. Moreover, we show that these
methods are as efficient as multi-phase matrix recovery algorithm. We provide
experimental evidence to illustrate the performance of our algorithm.

##  1 Introduction

Lifelong Adaptive Matrix Completion problem has been the center of the
attention of research for many recent years. Previously matrix completion has
been mainly studied in passive sampling setting by uniformly sampling and
observing some subset of entries of the matrix in the beginning of the
algorithm and applying nuclear norm minimization. In a series of work in Gross
(2011), Candes & Recht (2009) and Recht (2011) it has been shown that
uniformly sampling and observing
$\Omega((m+n)r\mathrm{max}(\mu_{0}^{2},\mu_{1}^{2})\log^{2}{n_{2}})$ many
entries are enough to recover an $m\times n$ matrix of rank $r$ with $\mu_{0}$
and $\mu_{1}$ column and row space coherence parameters using this idea.
Moreover, Candés & Tao (2010) showed this passive uniform sampling setting we
need at least $\Omega(mr\mu_{0}\log{n})$ observations to recover the matrix
exactly. Therefore, this result was concluding almost optimality of the
nuclear norm minimization where sampling method is passive.
However, in modern data analysis it has been show that adaptive sensing and
sampling methods are outperforming passive sampling methods Ramazanli (2022);
Haupt et al. (2011); Warmuth & Kuzmin (2008); Ramazanli et al. (2020). In a
follow up work the paper Krishnamurthy & Singh (2013) has proposed an adaptive
matrix completion algorithm that can actually recover underlying matrix using
at most $\mathcal{O}(\mu_{0}r^{1.5}\log{\frac{r}{\epsilon}})$ many
observation. Later authors further improved their existing results in
Krishnamurthy & Singh (2014) to show that the adaptive matrix completion
algorithm can actually successfully perform by just observing
$\mathcal{O}(\mu_{0}r\log^{2}{\frac{r}{\epsilon}})$ many observation.
The latest result has been even further optimized in the paper Balcan & Zhang
(2016) to present an adaptive algorithm which performs only with observation
complexity of just $\mathcal{O}(\mu_{0}r\log{\frac{r}{\epsilon}})$. Recently
Ramazanli & Poczos (2020) proposed a new parametrization method called
sparsity-number and a multi-phase algorithm which shows that using this new
method, the observation complexity can go as low as $\mathcal{O}((m+n-r)r)$.
The idea of this method to use the sparsest vector of column and row spaces
instead of coherence in the recovery process. Studying the sparsest vector of
the column and row space has been focus in many research problems (McCormick
(1983), Qu et al. (2015), Coleman & Pothen (1986)), and Ramazanli & Poczos
(2020) showed how to use this parameter efficiently for adaptive matrix
completion. In this paper, we propose a matrix completion method which merges
ideas of single phase algorithms such as Balcan & Zhang (2016) with the idea
of sparsity number and multi-phase method Ramazanli & Poczos (2020) which is
resulting in a single phased algorithm which is as efficient as multi-phased
one.
Low-rank matrix completion plays a significant role in many real-world
applications, including camera motion inferring, multi-class learning,
positioning of sensors, and gene expression analysis Balcan & Zhang (2016);
Krishnamurthy & Singh (2013). In gene expression analysis, the target matrix
represents expression levels for various genes across several conditions.
Measuring gene expression, however, is expensive, and we would like to
estimate the target matrix with a few observations as possible. In this paper,
we provide an algorithm that can be used for matrix completion from limited
data. Roughly speaking, to find each unknown expression level, we are supposed
to do multiple measurements. Each of the additional measurements has its extra
cost. Naturally, we aim to solve the entire problem using the least possible
measurement cost.

##  2 Preliminaries

In this section we start by providing notations and definitions those are used
throughout the paper. Then we will provide a single phase Balcan & Zhang
(2016) and multi-phase Poczos et al. (2020) matrix recovery algorithms.
Throughout the paper, we denote by $\mathbf{M}$ the target underlying $m\times
n$ sized rank-$r$ matrix that we want to recover. $\|x\|_{p}$ denote the
$L_{p}$ norm of a vector $x\in\mathbb{R}^{n}$. We call $x_{i}$ the $i$’th
coordinate of $x$. For any, $\Omega\subset[n]$ let $x_{\Omega}$ denote the
induced subvector of $x$ from coordinates $\Omega$. For any
$\mathbf{R}\subset[m]$, $\mathbf{M}_{\mathbf{R}:}$ stands for an
$|\mathbf{R}|\times n$ sized submatrix of $\mathbf{M}$ that rows are
restricted by $\mathbf{R}$. We define $\mathbf{M}_{:\mathbf{C}}$ in a similar
way for restriction with respect to columns. Intuitively,
$\mathbf{M}_{\mathbf{R}:\mathbf{C}}$ defined for
$|\mathbf{R}|\times|\mathbf{C}|$ sized submatrix of $\mathbf{M}$ with rows
restricted to $\mathbf{R}$ and columns restriced to $\mathbf{C}$. Moreover,
for the special case $\mathbf{M}_{i:}$ stands for $i$-th row and
$\mathbf{M}_{:j}$ stands for the $j$’th column. Similarly,
$\mathbf{M}_{i:\mathbf{C}}$ will represent the restriction of the row $i$ by
$C$ and $\mathbf{M}_{\mathbf{R}:j}$ represents restriction of the column $j$
by $\mathbf{R}$.
The following definitions are key parameters in this paper, first we define
coherence which is an important parameter in the matrix completion research
domain. Then, we visit the definition of the sparsity-number.

######  Definition 1.

Coherence parameter of a matrix $\mathbf{M}$ with column space $\mathbb{U}$ is
defined as following where $\mathcal{P}_{\mathbb{U}}$ represents the
orthogonal projection operator onto the subspace $\mathbb{U}$.

$\displaystyle\mu(\mathbb{U})=\frac{n}{r}\underset{1\leq j\leq
n}{\max}||\mathcal{P}_{\mathbb{U}}e_{j}||^{2},$

One can observe that if $e_{j}\in\mathbb{U}$ for some $j\in[n]$, then the
coherence will attain its maximum value: $\mu(\mathbb{U})=\frac{n}{r}$.

######  Definition 2.

Nonsparsity-number of a vector $x\in\mathbb{R}^{m}$ represented by $\psi(x)$
and defined as $\psi(x)=\|x\|_{0}$. Moreover, nonsparsity-number of an
$m\times n$ sized matrix $\mathbf{M}$ and a subspace
$\mathbb{U}\subseteq\mathbb{R}^{m}$ also represented by $\psi$ and defined as:

$\displaystyle\psi(\mathbf{M})=\mathrm{min}\\{\psi(x)|x=\mathbf{M}z\text{ and
}z\notin\mathrm{null}(\mathbf{M})\\}\hskip
42.67912pt\psi(\mathbb{U})=\mathrm{min}\\{\psi(x)|x\in\mathbb{U}\text{ and
}x\neq 0\\}$

Sparsity-number is just completion of nonsparsity-number and is donated by
$\overline{\psi}$:

$\displaystyle\overline{\psi}(x)=m-\psi(x)\hskip
56.9055pt\overline{\psi}(\mathbf{M})=m-\psi(\mathbf{M})\hskip
56.9055pt\overline{\psi}(\mathbb{U})=m-\psi(\mathbb{U})$

###  2.1 Single Phase Matrix Completion

Here we provide the generic framework for single-phase matrix completion
algorithm that has been used by Krishnamurthy & Singh (2013, 2014); Balcan &
Zhang (2016). The algorithm uniformly sample $d$-many entry in each column and
observe them. Based on the observation, if it has been decided the column is
linearly independent with previous columns then algorithm decides to observe
entire column. Initially, it has been shown that choosing
$d=\mathcal{O}(n\mu_{0}r^{1.5}\log{r})$ would be enough this algorithm to
perform successfully in Krishnamurthy & Singh (2013). Then, authors showed
setting $d=\mathcal{O}(n\mu_{0}r\log^{2}{r})$ would also be enough
Krishnamurthy & Singh (2014). In a following work Balcan & Zhang (2016)
optimized the complexity by reducing one of the $\log{r}$ factors from $d$ and
proving $d=\mathcal{O}(n\mu_{0}r\log{r})$ still works successfully with
probability $1-\epsilon$.

HN2016: Exact recovery Balcan & Zhang (2016).

Input: $d=\mathcal{O}(\mu_{0}r\log{\frac{r}{\epsilon}})$
Initialize: $k=0,\widehat{\mathbf{U}}^{0}=\emptyset$

1: Draw uniformly random entries $\Omega\subset[m]$ of size $d$

2: for $i$ from $1$ to $n$ do

3: if
$\|\mathbf{M}_{\Omega:i}-{\mathcal{P}_{\mathbf{\widehat{U}_{\Omega}^{k}}}}\mathbf{M}_{\Omega:i}\|>0$

4: Fully observe $\mathbf{M}_{:i}$

5:
$\widehat{\mathbf{U}}^{k+1}\leftarrow\widehat{\mathbf{U}}^{k}\cup\mathbf{M_{:i}}$,
Orthogonalize $\widehat{\mathbf{U}}^{k+1}$, $k=k+1$

6: otherwise:
$\widehat{\mathbf{M}}_{:i}=\widehat{\mathbf{U}}^{k}{\widehat{\mathbf{U}}^{k^{+}}_{\Omega}}\widehat{\mathbf{M}}_{\Omega:i}$

Output: $\widehat{\mathbf{M}}$

###  2.2 Multi Phase Matrix Completion

Here we provide multi-phase matrix completion algorithm from Ramazanli &
Poczos (2020). The are two main differences between this algorithm and
previous ones. First, previous algorithms goes through the columns of the
matrix just once, however, here we go through each column many times. Second,
there is a difference in terms of available pre-information. One can see that,
in order to perform previous algorithms, we need to have at least estimated
information of rank $r$ and column space coherence $\mu_{0}$. However, in
order to perform $\mathbf{ERRE}$ all we need to know is either column or row
space is not highly coherent.

ERRE: Exact recovery while rank estimation. Ramazanli & Poczos (2020)

Input: $T$-delay parameter at the end of algorithm
Initialize: $R=\emptyset,C=\emptyset,\widehat{r}=0,delay=0$

1: while $delay<T$ do

2: $delay=delay+1$

3: for $j$ from $1$ to $n$ do

4: Uniformly pick an unobserved entry $i$ from $\mathbf{M}_{:j}$

5: $\widehat{R}={R}\cup\\{i\\},\widehat{C}={C}\cup\\{j\\}$

6: If $\mathbf{M}_{\widehat{R}:\widehat{C}}$ is nonsingular :

7: Fully observe $\mathbf{M}_{:j}$ and $\mathbf{M}_{i:}$ and set
$R=\widehat{R}$ , $C=\widehat{C}$ , $\widehat{r}=\widehat{r}+1,delay=0$

8: Orthogonalize column vectors in $C$ and assign to $\mathbf{U}$

9: for each column $j\in[n]\setminus C$ do

10:
$\widehat{\mathbf{M}}_{:j}=\widehat{\mathbf{U}}{\widehat{\mathbf{U}}_{R:}}^{+}\widehat{\mathbf{M}}_{R:j}$

Output: $\widehat{\mathbf{M}},\widehat{r}$

The authors proved that the observation complexity could be bounded by the
following expression which is more optimal than previous methods and many
times it matches $\mathcal{O}((m+n-r)r)$ which is theoretical lower bound for
any kind of rank-$r$ matrix recovery algorithm.

######  Theorem 1.

Let $r$ be the rank of underlying $m\times n$ sized matrix $\mathbf{M}$ with
column space $\mathbb{U}$ and row space $\mathbb{V}$. Then, $\mathbf{ERRE}$
exactly recovers the underlying matrix $\mathbf{M}$ while estimating rank with
probability at least
$1-(\epsilon+e^{-T\frac{\psi(\mathbb{U})\psi(\mathbb{V})}{m}})$ using number
of observations at most:

$\displaystyle(m+n-r)r+Tn+\mathrm{min}\Big{(}2\frac{mn}{\psi(\mathbb{U})}\log{(\frac{r}{\epsilon})},\frac{\frac{2m}{\psi(\mathbb{U})}(r+2+\log{\frac{1}{\epsilon}})}{\psi(\mathbb{V})}n)\Big{)}$

##  3 Main Results

In this section, we provide single-phase column space recovery algorithm using
sparsity-number which can be easily extended to two-phase matrix recovery
algorithm. The algorithm works in a similar flavor to Balcan & Zhang (2016)
that it requires a preinformation regarding properties of the matrix in order
to execute.
However, there is one clear difference in terms of preinformation that, here
we require estimated information for sparsity number rather than coherence,
and also we use the information regarding row-space as well rather than simply
columns space. The reason behind this is, in deep mathemtical analysis of
aforementioned methods shows us all these methods work efficiently under the
condition row spaces are highly coherent and in-depth analysis of these
methods tells us that these algorithms designed for the highest value of row
space coherence, and they also work perfectly well for remaining cases.
In the following algorithm, we extend these algorithms to a method that can
enjoy the properties of row space; meanwhile, performing similarly good for
highly coherent row space matrices. Note taht, the algorithm below designed
and analyzed under the condition that estimated r, $\psi(\mathbb{U})$, and
$\psi(\mathbb{V})$ provided. However, under the condition just estimated $r$
and $\mu_{0}$ provided, we can set $\psi(\mathbb{V})$ to be equal to 1 and use
the inequality $\mu_{0}r>\frac{m}{k}$ to transfer the information of $\mu_{0}$
to $\psi(\mathbb{U})$. Then the algorithm below will perform as good as Balcan
& Zhang (2016) with observation complexity of
$mr+\mathcal{O}(nr\mu_{0}\log{\frac{r}{\epsilon}})$.

EREI: Exact recovery with estimated information

Input: $r,\psi(\mathbb{U}),\psi(\mathbb{V})$
Initialize:
$R=\emptyset,k=0,\widehat{\mathbf{U}}^{0}=\emptyset,d=\mathrm{min}\Big{(}2\frac{m}{\psi(\mathbb{U})}\log{(\frac{r}{\epsilon})},\frac{\frac{2m}{\psi(\mathbb{U}}(r+2+\log{\frac{1}{\epsilon}})}{\psi(\mathbb{V})}\Big{)}$

1: Draw uniformly random entries $\Omega\subset[m]$ of size $d$

2: for $i$ from $1$ to $n$ do

3: if
$\|\mathbf{M}_{\Omega:i}-{\mathcal{P}_{\mathbf{\widehat{U}_{\Omega}^{k}}}}\mathbf{M}_{\Omega:i}\|>0$

4: Fully observe $\mathbf{M}_{:i}$

5:
$\widehat{\mathbf{U}}^{k+1}\leftarrow\widehat{\mathbf{U}}^{k}\cup\mathbf{M_{:i}}$,
Orthogonalize $\widehat{\mathbf{U}}^{k+1}$, $k=k+1$

6: Select a row $a\in\Omega\setminus R$ that,
$\mathbf{\widehat{U}_{R\cup\\{a\\}}^{k+1}}$ is rank $k+1$ then $R\leftarrow
R\cup\\{a\\}$

7: Draw uniformly random entries $\Delta\subset[m]\setminus R$ of size $d$ and
$\Omega=\Delta\cup R$

8: Observe entire $\mathbf{M}_{R:}$

9: for $i$ from $1$ to $n$ do

10: if $\mathbf{M}_{:i}$ not fully observed then :
$\widehat{\mathbf{M}}_{:i}=\widehat{\mathbf{U}}^{k}{\widehat{\mathbf{U}}^{k^{+}}_{R:}}\widehat{\mathbf{M}}_{R:i}$

Output: $\widehat{\mathbf{M}}$x

Below we show the execution of the algorithm for $d=2$ and $r=1$. $a$ stands
for entries that is observed randomly, $a$ stands for entries that is observed
deterministically and $a$ stands for recovered entries without observing. Note
that $\Omega_{1}=\\{1,5\\},\Omega_{2}=\\{2,5\\},\Omega_{3}=\\{1,3\\}$ and
after 3-rd iteration $R$ becomes $\\{3\\}$ therefore, in fourth column
$\mathbf{M}_{3:4}$ observed deterministically, besides together random
observations $\Omega_{4}=\\{5,6\\}$. After all of the iterations completed, we
observe the entire $\mathbf{M}_{3:}$, and in the next step, we recover
remaining entries.

$\left[\begin{array}[]{cccc}\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{white}1}&\hbox{\pagecolor{white}3}&\hbox{\pagecolor{white}2}&\hbox{\pagecolor{white}3}\\\
\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{white}2}&\hbox{\pagecolor{white}6}&\hbox{\pagecolor{white}4}&\hbox{\pagecolor{white}6}\\\
\end{array}\right]\left[\begin{array}[]{cccc}\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{white}0}&\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{white}1}&\hbox{\pagecolor{white}3}&\hbox{\pagecolor{white}2}&\hbox{\pagecolor{white}3}\\\
\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{white}2}&\hbox{\pagecolor{white}6}&\hbox{\pagecolor{white}4}&\hbox{\pagecolor{white}6}\\\
\end{array}\right]\left[\begin{array}[]{cccc}\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{white}0}&\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{white}1}&\hbox{\pagecolor{white}3}&\hbox{\pagecolor{gray}2}&\hbox{\pagecolor{white}3}\\\
\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{white}2}&\hbox{\pagecolor{white}6}&\hbox{\pagecolor{white}4}&\hbox{\pagecolor{white}6}\\\
\end{array}\right]\left[\begin{array}[]{cccc}\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{darkgray}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{white}0}&\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{darkgray}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{white}1}&\hbox{\pagecolor{white}3}&\hbox{\pagecolor{darkgray}2}&\hbox{\pagecolor{darkgray}3}\\\
\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{darkgray}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{darkgray}0}&\hbox{\pagecolor{gray}0}\\\
\hbox{\pagecolor{white}2}&\hbox{\pagecolor{white}6}&\hbox{\pagecolor{darkgray}4}&\hbox{\pagecolor{gray}6}\\\
\end{array}\right]\left[\begin{array}[]{cccc}\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{darkgray}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{white}0}&\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{darkgray}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{darkgray}1}&\hbox{\pagecolor{darkgray}3}&\hbox{\pagecolor{darkgray}2}&\hbox{\pagecolor{darkgray}3}\\\
\hbox{\pagecolor{white}0}&\hbox{\pagecolor{white}0}&\hbox{\pagecolor{darkgray}0}&\hbox{\pagecolor{white}0}\\\
\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{darkgray}0}&\hbox{\pagecolor{gray}0}\\\
\hbox{\pagecolor{white}2}&\hbox{\pagecolor{white}6}&\hbox{\pagecolor{darkgray}4}&\hbox{\pagecolor{gray}6}\\\
\end{array}\right]\left[\begin{array}[]{cccc}\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{lightgray}0}&\hbox{\pagecolor{darkgray}0}&\hbox{\pagecolor{lightgray}0}\\\
\hbox{\pagecolor{lightgray}0}&\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{darkgray}0}&\hbox{\pagecolor{lightgray}0}\\\
\hbox{\pagecolor{darkgray}1}&\hbox{\pagecolor{darkgray}3}&\hbox{\pagecolor{darkgray}2}&\hbox{\pagecolor{darkgray}3}\\\
\hbox{\pagecolor{lightgray}0}&\hbox{\pagecolor{lightgray}0}&\hbox{\pagecolor{darkgray}0}&\hbox{\pagecolor{lightgray}0}\\\
\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{gray}0}&\hbox{\pagecolor{darkgray}0}&\hbox{\pagecolor{gray}0}\\\
\hbox{\pagecolor{lightgray}2}&\hbox{\pagecolor{lightgray}6}&\hbox{\pagecolor{darkgray}4}&\hbox{\pagecolor{gray}6}\\\
\end{array}\right]$

######  Theorem 2.

Let $r$ be the rank of underlying $m\times n$ sized matrix $\mathbf{M}$ with
column space $\mathbb{U}$ and row space $\mathbb{V}$. Then, with probability
$1-\epsilon$ the algorithm $\mathbf{EREI}$ exactly recovers the underlying
matrix $\mathbf{M}$ using number of observations at most

$\displaystyle(m+n-r)r+\mathrm{min}\Big{(}2\frac{mn}{\psi(\mathbb{U})}\log{(\frac{r}{\epsilon})},\frac{\frac{2m}{\psi(\mathbb{U})}(r+2+\log{\frac{1}{\epsilon}})}{\psi(\mathbb{V})}n\Big{)}$

###### Proof.

We show that the observation complexity is upper bounded both by

$(m+n-r)r+2\frac{m}{\psi(\mathbb{U})}\log{(\frac{r}{\epsilon})}$

and also

$(m+n-r)r+\frac{\frac{2m}{\psi(\mathbb{U})}(r+2+\log{\frac{1}{\epsilon}})}{\psi(\mathbb{V})}$

Throught the proof we will refer to the proof Theorem 4 - complexity bound
regarding $\mathbf{ERRE}$ of the paper Ramazanli & Poczos (2020) also we are
denoting $\psi(\mathbb{U})$ by $k$ to have the consistency with that proof.

###  3.1 $(m+n-r)r+2\frac{m}{\psi(\mathbb{U})}\log{(\frac{r}{\epsilon})}$

We first start with the case that if minimum of these two quantities is
$2\frac{m}{k}\log{(\frac{r}{\epsilon})}$. As the matrix has rank $r$, there
exists at least one set of linearly independent columns with $r$ columns. We
select the set of linear independent columns$-C$ that has lexicographically
smallest indices. We show sampling $2\frac{m}{k}\log{(\frac{1}{\epsilon})}$
entries from each column will give us the probability of at least
$1-r\epsilon$ correctly recovery.
From the step 5 of proof of the theorem $\mathbf{ERRE}$ we can see sampling
$2\frac{m}{k}\log{(\frac{1}{\epsilon})}$ entries from an active column, would
give guarantee of probability of at least $1-\epsilon$ detection of
independence. As $C$ is lexicographically smallest, each column is active on
the time entries sampled from it, and each of $r$ columns will succeed with
probability at least $1-\epsilon$. Therefore, using union bound, with
probability $1-r\epsilon$ all of the columns in C will succeed, which
guarantees the exact recovery.
Replacing $\epsilon$ by $\frac{\epsilon}{r}$ will conclude the result that
sampling $2\frac{m}{k}\log{(\frac{r}{\epsilon})}$ from each column will
guarantees the correctness of the algorithm with probability at least
$1-\epsilon$.

###  3.2
$(m+n-r)r+\frac{\frac{2m}{\psi(\mathbb{U})}(r+2+\log{\frac{1}{\epsilon}})}{\psi(\mathbb{V})}$.

From the follow up of lemma 13 in Ramazanli & Poczos (2020), we conclude that
in a process of $\frac{k}{m}$ probability success and $1-\frac{k}{m}$
probability of failure, having
$\frac{2m}{k}\big{(}r+2+\log{\frac{1}{\epsilon}}\big{)}$ trial is enough to
guarantee getting $r$ many success with at least probability $1-\epsilon$.
Failure probability of the algorithm $\mathbf{EREI}$ is equal to failing
finding $r$ linearly independent columns. Consider following equation:

$\displaystyle P\big{(}\mathbf{EREI}fails\big{)}$
$\displaystyle=P\Big{(}\mathbf{EREI}fails\text{ and
TNAO}\geq\frac{2m}{k}\big{(}r+2+\log{\frac{1}{\epsilon}}\big{)}\Big{)}$
$\displaystyle+P\Big{(}\mathbf{EREI}fails\text{ and
TNAO}<\frac{2m}{k}\big{(}r+2+\log{\frac{1}{\epsilon}}\big{)}\Big{)}$

where we denote TNAO as total number of active observations. Recall that we
call an observation active, if it is active in the execution time (the column
is still not contained in the current column space). Intuitively we represent
NAO by number of active observations executed by the algorithm $\mathbf{EREI}$
in the given specific time.
From lemma 18, if there exists an active column, then there exists at least
$t$ many active columns. Therefore, failure of the algorithm is equivalent to
existence of an active column when algorithm terminates. Moreover each of our
observations in those columns were active observations and considering the
fact that we observed
$\frac{\frac{2m}{k}\big{(}r+2+\log{\frac{1}{\epsilon}}\big{)}}{t}$ many
entries in each of them, total number of active observations is at least

$\displaystyle
t\frac{\frac{2m}{k}\big{(}r+2+\log{\frac{1}{\epsilon}}\big{)}}{t}=\frac{2m}{k}\Big{(}r+2+\log{\frac{1}{\epsilon}}\Big{)}$

Therefore
$P\Big{(}\text{NAO}<\frac{2m}{k}\big{(}r+2+\log{\frac{1}{\epsilon}}\big{)}|\mathbf{EREI}fails\Big{)}=0$
and using Bayesian rule we conclude

$\displaystyle P\Big{(}\mathbf{EREI}fails\text{ and
TNAO}\leq\frac{2m}{k}\big{(}r+2+\log{\frac{1}{\epsilon}}\big{)}\Big{)}=0$

Then, following equation simplly satisfied:

$\displaystyle
P\big{(}\mathbf{EREI}fails\big{)}=P\Big{(}\mathbf{EREI}fails\text{ and
TNAO}\geq\frac{2m}{k}\big{(}r+2+\log{\frac{1}{\epsilon}}\big{)}\Big{)}$

We can observe the following inequality as $\mathbf{EREI}$ may tamporarily
fail at the point that the number of active observations is
$\frac{2m}{k}\big{(}r+2+\log{\frac{1}{\epsilon}}\big{)}$ but it can succeed
finding remaining independent columns later during the execution:

$\displaystyle P\Big{(}\mathbf{EREI}\text{ currently fail when NAO =}$
$\displaystyle\frac{2m}{k}\big{(}r+2+\log{\frac{1}{\epsilon}}\big{)}\Big{)}\geq$
$\displaystyle P\Big{(}\mathbf{EREI}fails\text{ and
TNAO}\geq\frac{2m}{k}\big{(}r+2+\log{\frac{1}{\epsilon}}\big{)}\Big{)}$

Therefore we conclude that:

$\displaystyle P\big{(}\mathbf{EREI}fails\big{)}\leq
P\Big{(}\mathbf{EREI}\text{ currently fail when NAO =}$
$\displaystyle\frac{2m}{k}\big{(}r+2+\log{\frac{1}{\epsilon}}\big{)}\Big{)}$

Remember the fact that at each active observation probability of
$\mathbf{EREI}$ detecting linear independence of is larger or equal than
$\frac{k}{m}$. From the previous discussion if the probability is exactly
equal to $\frac{k}{m}$ then still not finding $r$ linearly independent column
at $\frac{2m}{k}\big{(}r+2+\log{\frac{1}{\epsilon}}\big{)}$ observations is
less than $\epsilon$. Therefore, $\mathbf{EREI}$ not detecting r linearly
independent column after
$\frac{2m}{k}\big{(}r+2+\log{\frac{1}{\epsilon}}\big{)}$ observations is
smaller than $\epsilon$, which is equivalent to

$\displaystyle P\big{(}\mathbf{EREI}fails\big{)}\leq\epsilon$

as desired. ∎

######  Corollary 3.

Lets assume that we have estimated values of rank $r$, column space coherence
number $\mu_{0}$ and estimated row space sparsity-number is
$\psi(\mathbb{V})$. Then, if $\psi(\mathbb{V})$ is $\mathcal{O}(1)$ then
observation complexity is buounded by
$(m+n-r)r+\mathcal{O}(nr\mu_{0}\log{\frac{r}{\epsilon}})$ and if
$\psi(\mathbb{V})$ is $\mathcal{O}(n)$ then with probabiliy
$\frac{1}{2^{\mathcal{O}(r)}}$ the bound is $\mathcal{O}((m+n-r)r)$ which is
theoretical lower bound for exact completion problem.

## References

  * Balcan & Zhang (2016) M.F Balcan and H. Zhang.  Noise-tolerant life-long matrix completion via adaptive sampling.  _Advances in Neural Information Processing Systems_ , 2016.
  * Candes & Recht (2009) E. J. Candes and B. Recht.  Exact matrix completion via convex optimization.  _Foundations of Computational Mathematics_ , 2009.
  * Candés & Tao (2010) E. J. Candés and T. Tao.  The power of convex relaxation: Near-optimal matrix completion.  _IEEE Transactions on Information Theory_ , 2010.
  * Coleman & Pothen (1986) T. F. Coleman and A. Pothen.  The null space problem i. complexity.  _SIAM Journal on Algebraic Discrete Methods_ , 1986.
  * Gross (2011) D. Gross.  Recovering low-rank matrices from few coefficients in any basis.  _Information Theory, IEEE Transactions on._ , 2011.
  * Haupt et al. (2011) J. Haupt, R. Castro, and R. Nowak.  Adaptive sampling for sparse detection and estimation.  _IEEE Transactions on Information Theory_ , 2011.
  * Krishnamurthy & Singh (2013) A. Krishnamurthy and A. Singh.  Low-rank matrix and tensor completion via adaptive sampling.  In _Advances in Neural Information Processing Systems_ , 2013.
  * Krishnamurthy & Singh (2014) A. Krishnamurthy and A. Singh.  On the power of adaptivity in matrix completion and approximation.  _arXiv preprint_ , 2014.
  * McCormick (1983) S. T. McCormick.  A combinatorial approach to some sparse matrix problems.  _DTIC Document_ , 1983.
  * Poczos et al. (2020) Barnabas Poczos et al.  Optimal adaptive matrix completion.  _arXiv preprint arXiv:2002.02431_ , 2020.
  * Qu et al. (2015) Qing Qu, Ju Sun, and John Wright.  Finding a sparse vector in a subspace: linear sparsity using alternating directions.  _IEEE TransactionRANSACTION onON information theoryINFORMATION THEORY_ , 2015.
  * Ramazanli (2022) Ilqar Ramazanli.  Performance of distribution regression with doubling measure under the seek of closest point.  _arXiv preprint arXiv:2203.00155_ , 2022.
  * Ramazanli & Poczos (2020) Ilqar Ramazanli and Barnabas Poczos.  Optimal exact matrix completion under new parametrization.  _arXiv preprint arXiv:2002.02431_ , 2020.
  * Ramazanli et al. (2020) Ilqar Ramazanli, Han Nguyen, Hai Pham, Sashank J Reddi, and Barnabas Poczos.  Adaptive sampling distributed stochastic variance reduced gradient for heterogeneous distributed datasets.  _arXiv preprint arXiv:2002.08528_ , 2020.
  * Recht (2011) B. Recht.  A simpler approach to matrix completion.  _The Journal of Machine Learning Research_ , 2011.
  * Warmuth & Kuzmin (2008) M. K. Warmuth and D. Kuzmin.  Randomized online pca algorithms with regret bounds that are logarithmic in the dimension.  _Journal of Machine Learning Research_ , 2008.