# PROMPT: Parallel Iterative Algorithm for $\ell_{p}$ norm linear regression
via Majorization Minimization with an application to semi-supervised graph
learning

R. Jyothi and P. Babu

###### Abstract

In this paper, we consider the problem of $\ell_{p}$ norm linear regression,
which has several applications such as in sparse recovery, data clustering,
and semi-supervised learning. The problem, even though convex, does not enjoy
a closed-form solution. The state-of-the-art algorithms are iterative but
suffer from convergence issues, i.e., they either diverge for $p>3$ or the
convergence to the optimal solution is sensitive to the initialization of the
algorithm. Also, these algorithms are not generalizable to every possible
value of $p$. In this paper, we propose an iterative algorithm : Parallel
IteRative AlgOrithM for $\ell_{\textbf{P}}$ norm regression via MajorizaTion
Minimization (PROMPT) based on the principle of Majorization Minimization and
prove that the proposed algorithm is monotonic and converges to the optimal
solution of the problem for any value of $p$. The proposed algorithm can also
parallelly update each element of the regression variable, which helps to
handle large scale data efficiently, a common scenario in this era of data
explosion. Subsequently, we show that the proposed algorithm can also be
applied for the graph based semi-supervised learning problem. We show through
numerical simulations that the proposed algorithm converges to the optimal
solution for any random initialization and also performs better than the
state-of-the-art algorithms in terms of speed of convergence. We also evaluate
the performance of the proposed algorithm using simulated and real data for
the graph based semi-supervised learning problem.

###### Index Terms:

$\ell_{p}$ norm linear regression, Majorization Minimization, Parallel
algorithm, Graph based semi-superivsed learning

##  I Introduction and Problem Formulation

Linear regression is a statistical model which is used to estimate the
relationship between a dependent variable and one or more independent
variables [1]. This technique is commonly used in data forecasting [2], time
series analysis [3] and for risk assessment [4]. In linear regression, the
dependent variable $y$ is related to the unknown independent variables
${\boldsymbol{x}}=[x_{1},x_{2},\ldots,x_{n}]^{T}\in\mathbb{R}^{n}$ using a
linear model as follows:

$\begin{array}[]{ll}y=a_{1}x_{1}+a_{2}x_{2}+\cdots+a_{n}x_{n}+e={\boldsymbol{a}}^{T}{\boldsymbol{x}}+e\end{array}$
(1)

where $a_{i}$’s are the known regression parameters and $e$ is the unknown
noise usually modeled as a Gaussian random variable. Suppose we have $m$ such
observations, then the model in (1) can be re-written compactly as:

$\begin{array}[]{ll}{\mathbf{y}}={\mathbf{A}}{\boldsymbol{x}}+{\boldsymbol{e}}\end{array}$
(2)

where ${\mathbf{y}}\in\mathbb{R}^{m}$ is obtained by stacking the observations
$\\{y_{i}\\}_{i=1}^{m}$,
${\mathbf{A}}=[{\boldsymbol{a}}_{1},{\boldsymbol{a}}_{2},\ldots,{\boldsymbol{a}}_{m}]^{T}\in\mathbb{R}^{m\times
n}$ and ${\boldsymbol{e}}=[e_{1},e_{2},\ldots,e_{m}]^{T}\in\mathbb{R}^{m}$ is
the unknown error vector. Usually, $m>n$ i.e., there are more equations than
unknowns and hence the linear system of equations in (2) does not have any
solution. In such a case, a classical approach adopted to solve the set of
equations in (2) is to find an ${\boldsymbol{x}}$ such that the sum of the
squares of the error i.e., the $\ell_{2}$ norm of the error vector is minimum:

$\begin{array}[]{ll}\underset{{\boldsymbol{x}}}{\rm
arg\>min}\>\\{f({\boldsymbol{x}})=\|{\mathbf{y}}-{\mathbf{A}}{\boldsymbol{x}}\|^{2}_{2}\\}\end{array}$
(3)

where $\|\cdot\|_{2}$ is the Euclidean norm. The estimate ${\boldsymbol{x}}$
obtained by solving the above problem is known as the least-squares estimator
and is given by
${\boldsymbol{x}}=\left({\mathbf{A}}^{T}{\mathbf{A}}\right)^{-1}\left({\mathbf{A}}^{T}{\mathbf{y}}\right)$.
However, it has been observed that the performance of the least-squares
estimator is degraded when the underlying noise is not Gaussian distributed
[5]. To overcome this issue, the authors in [6] and [7] show that the estimate
${\boldsymbol{x}}$ obtained by minimizing the $\ell_{p}$ norm of the error
vector, for $p\neq 2$, is less sensitive to the noise distribution. In
practice, the noise distribution may be non-Gaussian and hence in this paper
we focus on solving the $\ell_{p}$ norm linear regression problem which is
defined as:

$\begin{array}[]{ll}\textrm{LR:}\quad\underset{{\boldsymbol{x}}}{\rm
arg\>min}\>\\{f_{{}_{LR}}({\boldsymbol{x}})=\|{\mathbf{y}}-{\mathbf{A}}{\boldsymbol{x}}\|_{p}\\}\end{array}$
(4)

Even though the problem in (4) is convex and differentiable, one cannot obtain
a closed-form solution for it using the KKT conditions. Hence, the state-of-
the art methods used to solve the problem in (4) are iterative in nature. The
benchmarking algorithm used to solve the $\ell_{p}$ norm problem for
$p\in[1,\infty)$ is the Iterative Reweighted Least Squares (IRLS) algorithm
([8, 9, 10, 11]). It is a simple iterative algorithm which replaces the
objective function in (4) with a weighted $\ell_{2}$ norm function at every
iteration, which when minimized admits a closed-form solution. However, it was
found that this method diverges for $p>3$ ([10],[12]). Since then, many
algorithms have been proposed to overcome the issue of convergence faced by
IRLS algorithm ([13, 14, 15]). We review these methods briefly in the next
section. Apart from being a robust regressor, the $\ell_{p}$ norm linear
regression for $p>2$ has applications in semi-supervised learning commonly
used in classification problems ([16, 17, 18]), which we discuss next.
In today’s world, classification problems are becoming pervasive wherein the
goal is to categorize the given data into its corresponding class using a
mapping function [1]. In the supervised learning setting, the training data
set (which consists of data points and its corresponding class labels) are
used to learn the mapping function which can be used later to predict the
class label of an unknown test data point. However, one setback of this
approach is that it requires large amount of labeled data points to learn the
mapping function. Data labeling is a laborious task (e.g. writing a transcript
for speech recognition) and also requires expert input (deciding whether a
brain scan is healthy or not). However, one can easily obtain large amounts of
unlabeled data. Semi-supervised learning utilizes both the labeled and
unlabeled data points to learn the mapping function. By doing so, semi-
supervised learning learns a better mapping function when compared to the
supervised learning which uses only the labeled data points [12].
Under the semi-supervised learning approach, graph based semi-supervised
learning is commonly used, wherein the labels of the unknown data points are
learned by propagating the known label information using graphs. Suppose we
are given $n$ data points ${\boldsymbol{v}}_{i}\in\mathbb{R}^{d}$ out of which
$u$ are unlabeled data points and $l$ are labeled points with labels
$h_{i}=f({\boldsymbol{v}}_{i})$ (where $f$ is an unknown real-valued
function). Generally, the number of labeled points are much smaller than the
number of unlabeled points i.e., $l<<u$. An undirected graph
${\boldsymbol{G}}=(V,E)$ with vertex set $V=\\{v_{1},v_{2},\ldots v_{n}\\}$
are used to represent the $n$ data points. The subset ${L}\subset V$ of the
vertex set are used to represent the labeled points and the remaining vertices
$V\setminus L$ represents the unlabeled points. The edges $E$ of the graph are
equipped with non-negative edge weights ${W}=\\{w_{xy}\\}_{x,y\in{E}}$, which
are chosen such that $w_{xy}\approx 1$ when $x$ is similar to $y$ and
$w_{xy}\approx 0$ when $x$ and $y$ are dissimilar. The task of graph based
semi-supervised learning is to assign labels to the remaining vertices without
changing the labels at the vertex subset $L$. Since, this problem has
infinitely many solutions, one usually makes the _semi-supervised smoothness
assumption_ which states that if two data points ${\boldsymbol{v}}_{1}$,
${\boldsymbol{v}}_{2}$ in a high density region are close then so should be
their corresponding outputs $h_{1}$, $h_{2}$; hence the function $f$ is
required to be smooth in the high density regions [19]. A common approach used
to label the data points at remaining vertices such that the mapping function
does not change rapidly in the high density regions and also agrees with the
labels at the vertices $L$ is by using the $\ell_{2}$-based Laplacian
regularizer, which is given as:

$\begin{array}[]{ll}\hat{f}(v)=\underset{f}{\rm
arg\>min}\>\displaystyle\sum_{i,j\in
E}w_{ij}\left(f(v_{i})-f(v_{j})\right)^{2}\\\ \hskip
28.45274pt{\rm{subject\>to}}\>f(v_{i})=h_{i},i\in{L}\end{array}$ (5)

The above model was first introduced by [20]. From an intuitive perspective,
minimizing the objective function in (5) will ensure that the similar points
(for which $w_{ij}$ is larger) will have similar labels, making
$\left(f(v_{i})-f(v_{j})\right)^{2}$ smaller, while for dissimilar points
($w_{ij}$ is small) and hence are allowed to have dissimilar labels. However,
it was found that when the number of unlabeled data points $u$ is much greater
than the number of labeled data points $l$, the learned function $\hat{f}(v)$
becomes constant everywhere, with sharp spikes near the labeled data. We
illustrate the same in Fig.1.a. for the number of labeled points $l$ equal to
$2$ and number of unlabeled points $u$ equal to $105$. In Fig. 1, the $X$ and
the $Y$-axis represent the values taken by the data points and the $Z$-axis
represent the value of the function surface at the corresponding data points.
From Fig.1.a it can be seen that for $p=2$ the function $\hat{f}(v)$ is
constant everywhere with sharp spikes near the labeled data point. We also
illustrate the function surface for $p=2.5$ and $p=\infty$ in Fig.1.b and Fig.
1.c, respectively.

---
Figure 1: Illustration of the function surface $\hat{f}(v)$ for different
values of $p$. For $p=2$ the surface is constant with peaks near the labeled
points and for $p>2$ the surface becomes smoother. This example was generated
for two labeled and 105 unlabeled data points [12].

From Fig. 1.b and Fig. 1.c it can be seen that as the value of $p$ increases
the function becomes smooth. This observation led to the formulation of
$\ell_{p}$-based Laplacian regularizer problem [21]:

$\begin{array}[]{ll}\underset{f}{\rm arg\>min}\>\displaystyle\sum_{i,j\in
E}w_{ij}\left|f(v_{i})-f(v_{j})\right|^{p}\\\ \hskip
8.53581pt{\rm{subject\>to}}\>f(v_{i})=h_{i},i\in L\end{array}$ (6)

The above problem can be cast as an $\ell_{p}$ norm linear regression problem
(which we will show in sec. IV). Hence, the benchmarking algorithm used to
solve the problem in (6) is the modified IRLS algorithm [15]. The paper is
organized as follows. In the next section, we discuss the methods used to
solve the $\ell_{p}$ norm regression problem and also list the contributions
made in this manuscript. In sec. III, we give an overview of Majorization
Minimization (MM) procedure and in sec. IV we propose a parallel algorithm to
solve the problem in (4) using the MM principle for any value of $p$. At the
end of the same section we discuss the application of the proposed algorithm
for the graph based semi-supervised problem. In sec. V, we compare the
proposed algorithm with the state-of-the art methods via computer simulations
and finally conclude the paper in sec. VI.

##  II Related Work and Contributions

The conventional IRLS algorithm used to solve the $\ell_{p}$ norm problem was
first developed by Karlovitz [8] and also independently by several other
researchers ([9, 10, 11]). A brief history of IRLS algorithm can be found in
[22]. To compute the value of the next iterate ${\boldsymbol{x}}^{k+1}$, the
IRLS algorithm solves the following weighted least-squares problem:

$\begin{array}[]{ll}{\boldsymbol{x}}^{k+1}=\quad\underset{{\boldsymbol{x}}}{\rm
arg\>min}\>({\mathbf{A}}{\boldsymbol{x}}-{\mathbf{y}})^{T}{\boldsymbol{W}}^{k}({\mathbf{A}}{\boldsymbol{x}}-{\mathbf{y}})\end{array}$
(7)

where ${\boldsymbol{W}}$ is a diagonal matrix with diagonal elements
$\\{w_{1},w_{2}\cdots w_{m}\\}$ and ${\boldsymbol{W}}^{k}$ is the value taken
by ${\boldsymbol{W}}$ at the $k^{th}$ iteration. The IRLS algorithm starts
with ${\boldsymbol{W}}^{0}={\mathbf{I}}$ and solves the problem in (7) whose
solution is the least-squares estimator
${\boldsymbol{x}}^{0}=\left({\mathbf{A}}^{T}{\mathbf{A}}\right)^{-1}\left({\mathbf{A}}^{T}{\mathbf{y}}\right)$.
It then updates the weights matrix
${\boldsymbol{W}}^{k}=\textrm{diag}(|{\mathbf{A}}{\boldsymbol{x}}^{k}-{\mathbf{y}}|^{p-2})$
and solves the weighted least-squares problem in (7) whose solution is
${\boldsymbol{x}}^{k}=({\mathbf{A}}^{T}{\boldsymbol{W}}^{k}{\mathbf{A}})^{-1}{\mathbf{A}}^{T}{\boldsymbol{W}}^{k}{\mathbf{y}}$
and this is repeated until convergence. The pseudocode of IRLS algorithm is
shown in Table 1.

Table 1: Pseudocode of IRLS algorithm
---
Input: Noisy observed data ${\mathbf{y}}$, data matrix ${\mathbf{A}}$ and $p$
Initialize: Set _k_ = 0. Initialize
${{\boldsymbol{x}}^{0}}=\left({\mathbf{A}}^{T}{\mathbf{A}}\right)^{-1}\left({\mathbf{A}}^{T}{\mathbf{y}}\right)$.
Repeat:
1) Update
${\boldsymbol{W}}^{k}=\textrm{diag}(|{\mathbf{A}}{\boldsymbol{x}}^{k}-{\mathbf{y}}|^{p-2})$
2)
${\boldsymbol{x}}^{k+1}=({\mathbf{A}}^{T}{\boldsymbol{W}}^{k}{\mathbf{A}})^{-1}{\mathbf{A}}^{T}{\boldsymbol{W}}^{k}{\mathbf{y}}$
$k\leftarrow k+1$
until
$\left|\dfrac{{f_{{}_{\textrm{LR}}}({{{{\boldsymbol{x}}}^{k}}})}-{f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}^{k-1}})}}{{f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}^{k-1}})}}\right|\leq
10^{-3}$

The IRLS algorithm has several drawbacks. Firstly, the IRLS algorithm involves
computing the inverse of a square matrix of size $n$ at every iteration -
making the algorithm computationally expensive for large $n$. Also, at any
iteration, if ${\mathbf{y}}$ becomes equal to
${\mathbf{A}}{\boldsymbol{x}}^{k}$, the matrix
$({\mathbf{A}}^{T}{\boldsymbol{W}}^{k}{\mathbf{A}})$ becomes ill-conditioned
and the algorithm will run into numerical issues. Another drawback of IRLS
algorithm is that it fails to converge for $p>3$ ([10],[12]). To overcome the
issue of convergence, the authors in [8] calculated the value of the next
iterate by the following weighted combination of the previous iterate
${\boldsymbol{x}}^{k-1}$ and the IRLS iterate $\tilde{{\boldsymbol{x}}}^{k}$
with the updated weights:

$\begin{array}[]{ll}{\boldsymbol{x}}^{k+1}=q\tilde{{\boldsymbol{x}}}^{k}+(1-q){\boldsymbol{x}}^{k-1}\end{array}$
(8)

However, this method is very slow because one has to choose an optimal value
for $q$ at every iteration. The authors in [15] proposed a modified IRLS
algorithm wherein to avoid ill-conditioned matrix, they systematically added
non-zero values to the diagonal elements of the matrix ${\boldsymbol{W}}^{k}$.
Also, to improve the stability of the IRLS algorithm, the next iterate
${\boldsymbol{x}}^{k+1}$ is calculated by doing a line search along the line
joining the previous iterate ${\boldsymbol{x}}^{k-1}$ and the standard IRLS
iterate with the modified weights. Homotopy based approaches have also been
proposed to solve the $\ell_{p}$ norm regression problem [13], wherein first a
simpler optimization problem is solved i.e., the $\ell_{p}$ norm problem for
$p=2$ is solved and then $p$ is slowly increased until the desired minimum of
the original objective function is reached. The authors in [23] used interior
point methods to solve the $\ell_{p}$ norm problem. However, this method
converges very slowly for large dimension of input data [15].

In this paper, we propose a parallel iterative algorithm to solve the
$\ell_{p}$ norm regression problem for any $p$ based on the Majorization
Minimization (MM) procedure (which will be introduced briefly in the next
section) and show through numerical simulations that the proposed algorithm
converges to the stationary point of the problem for any random initialization
and enjoys faster speed of convergence. The major contributions of this paper
are:

  1. 1.

A MM based parallel algorithm - Parallel IteRative AlgOrithM for
$\ell_{\textbf{P}}$ norm regression via MajorizaTion Minimization (PROMPT) is
proposed to solve the $\ell_{p}$ norm regression problem for $p\in[1,\infty]$.
The proposed algorithm can update each element of ${\boldsymbol{x}}$ parallely
- which is useful for solving large dimensional problem.

  2. 2.

We prove that the proposed algorithm converges to the stationary point of the
problem.

  3. 3.

We show that the proposed algorithm can be applied for the graph based semi-
supervised learning problem.

  4. 4.

Numerical simulations are conducted to compare the proposed algorithm with the
state-of-the-art algorithms.

##  III Majorization Minimization

In this section, we give an overview of MM procedure which has been applied to
develop PROMPT. To explain the MM framework, we consider the following
optimization problem:

$\begin{array}[]{ll}\underset{{\boldsymbol{x}}\in\chi}{\rm
minimize}\>f({\boldsymbol{x}})\end{array}$ (9)

where $f({\boldsymbol{x}})$ could be either a convex or non-convex function
and $\chi$ is the constraint set. Under the MM framework, at every iteration,
a surrogate function $g({\boldsymbol{x}}|{\boldsymbol{x}}^{k})$ which
majorizes the function $f({\boldsymbol{x}})$ is constructed and minimized i.e.

${\boldsymbol{x}}^{k+1}\in\underset{{\boldsymbol{x}}\in\chi}{\rm
arg\>min}\quad g\left({\boldsymbol{x}}|{\boldsymbol{x}}^{k}\right)$ (10)

where ${\boldsymbol{x}}^{k+1}$ is the value taken by ${\boldsymbol{x}}$ at the
$(k+1)^{th}$ iteration. A surrogate function
$g({\boldsymbol{x}}|{\boldsymbol{x}}^{k})$ is said to majorize a function
$f({\boldsymbol{x}})$ if it satisfies the following properties:

$g\left({\boldsymbol{x}}^{k}|{\boldsymbol{x}}^{k}\right)=f\left({\boldsymbol{x}}^{k}\right)$
(11) $g\left({\boldsymbol{x}}|{\boldsymbol{x}}^{k}\right)\geq
f\left({\boldsymbol{x}}\right)$ (12)

Using (10), (11) and (12) it can be shown that the sequence of points
$\\{{\boldsymbol{x}}^{k}\\}$ generated by the MM procedure monotonically
decrease the objective function:

$\begin{array}[]{ll}f({\boldsymbol{x}}^{k+1})\leq
g({\boldsymbol{x}}^{k+1}|{\boldsymbol{x}}^{k})\leq
g({\boldsymbol{x}}^{k}|{\boldsymbol{x}}^{k})=f({\boldsymbol{x}}^{k})\end{array}$
(13)

where the first inequality and the last equality are obtained by using (11)
and (12). The second inequality is by (10). The computational complexity and
the convergence speed of the algorithm depends on the choice of the surrogate
function. For instance, in the case of a multivariate optimization problem, a
surrogate function could make the parameters separable and hence each of them
could be updated parallely - which is particularly useful for a large scale
problem. An overview of the various surrogate functions used can be found in
([24], [25]).

##  IV Proposed Algorithm for $\ell_{p}$ norm linear regression problem

In this section, we first propose an iterative algorithm PROMPT based on the
MM principle to solve the $\ell_{p}$ norm linear regression problem for
$p\in(1,\infty)$. The proposed algorithm can parallely update each element of
${\boldsymbol{x}}$ which is useful to handle large scale data efficiently. We
then discuss the extensions of the proposed algorithm to solve the $\ell_{1}$
and the $\ell_{\infty}$ norm linear regression problems. Next, we show that
the proposed algorithm converges to the stationary point of the $\ell_{p}$
norm linear regression problem. Finally, at the end of the section, we discuss
the application of the proposed algorithm for the graph based semi-supervised
learning problem.

###  IV-A Parallel IteRative AlgOrithM for $\ell_{\textbf{P}}$ norm regression
via MajorizaTion Minimization (PROMPT)

The objective function $f_{{}_{LR}}({\boldsymbol{x}})$ in (4) is not separable
in each element of ${\boldsymbol{x}}$, which can be observed by first
rewriting $f_{{}_{LR}}({\boldsymbol{x}})$ as:

$\begin{array}[]{ll}f_{{}_{LR}}({\boldsymbol{x}})=\|{\mathbf{y}}-{\mathbf{A}}{\boldsymbol{x}}\|_{p}=\displaystyle\sum_{i=1}^{m}|y_{i}-{\boldsymbol{a}}_{i}^{T}{\boldsymbol{x}}|^{p}\end{array}$
(14)

where ${\boldsymbol{a}}_{i}\in\mathbf{R}^{n}$ is the $i^{th}$ row of
${\mathbf{A}}$ matrix and $y_{i}$ is the $i^{th}$ element of ${\mathbf{y}}$.
Expanding $|y_{i}-{\boldsymbol{a}}_{i}^{T}{\boldsymbol{x}}|^{p}$, one gets
terms coupled in the elements of ${\boldsymbol{x}}$ which makes the parallel
minimization of $f_{{}_{LR}}({\boldsymbol{x}})$ challenging. In this section,
we develop a parallel algorithm using the MM principle in which we form a
surrogate function $g(x_{j}|{\boldsymbol{x}}^{k})$ which majorizes
$f_{{}_{LR}}({\boldsymbol{x}})$. The surrogate function is separable in each
element of ${\boldsymbol{x}}$ and hence each element of ${\boldsymbol{x}}$ can
be updated parallely. To develop the surrogate function we make use of the
following lemma:

###### Lemma IV.1

Given any $\tilde{{\boldsymbol{x}}}={\tilde{{\boldsymbol{x}}}^{k}}$, the
function $\left|{{\boldsymbol{c}}^{T}}\tilde{{\boldsymbol{x}}}\right|^{p}$
($p\geq 1$), where
${\boldsymbol{c}}\>\rm{and}\>\tilde{{\boldsymbol{x}}}\in\mathbf{R}^{n+1}$, can
be upper bounded as:

$\begin{array}[]{ll}\left|\left({{\boldsymbol{c}}^{T}}\tilde{{\boldsymbol{x}}}\right)\right|^{p}\leq\displaystyle\sum_{j=1}^{n+1}\dfrac{1}{n+1}\left|\left(n+1\right)c_{j}\left(\tilde{x}_{j}-\tilde{x}_{j}^{k}\right)+{\boldsymbol{c}}^{T}{\tilde{{\boldsymbol{x}}}^{k}}\right|^{p}\end{array}$
(15)

###### Proof:

We replicate the proof from [24] for the sake of clarity. Note that the
function $\left|\left(\cdot\right)\right|^{p}$ ($p\geq 1$) is convex and hence
by using the Jensen’s inequality we get:

$\begin{array}[]{ll}\left|\left(\displaystyle\sum_{j=1}^{n+1}\dfrac{s_{j}}{n+1}\right)\right|^{p}\leq\displaystyle\sum_{j=1}^{n+1}\dfrac{|\left(s_{j}\right)|^{p}}{n+1}\end{array}$
(16)

Letting
$s_{j}=(n+1)c_{j}\left(\tilde{x}_{j}-\tilde{x}_{j}^{k}\right)+{\boldsymbol{c}}^{T}(\tilde{{\boldsymbol{x}}}^{k})$
and substituting it in (16), the inequality in (15) is achieved. ∎

Let ${\boldsymbol{c}}_{i}=[c_{i,1}\cdots
c_{i,n+1}]^{T}=[-{\boldsymbol{a}}_{i}\>\>y_{i}]^{T}$ and let
$\tilde{{\boldsymbol{x}}}=[{\boldsymbol{x}},1]^{T}$. Then the objective
function $f_{{}_{LR}}({\boldsymbol{x}})$ in (14) can be rewritten as:

$\begin{array}[]{ll}f_{{}_{LR}}({\boldsymbol{x}})=\displaystyle\sum_{i=1}^{m}|y_{i}-{\boldsymbol{a}}_{i}^{T}{\boldsymbol{x}}|^{p}=\displaystyle\sum_{i=1}^{m}|{\boldsymbol{c}}_{i}^{T}\tilde{{\boldsymbol{x}}}|^{p}\end{array}$
(17)

Using lemma IV.1, we majorize the term
$|{\boldsymbol{c}}_{i}^{T}\tilde{{\boldsymbol{x}}}|^{p}$ and hence majorize
the objective function $f_{{}_{LR}}({\boldsymbol{x}})$. At any given
${\boldsymbol{x}}^{k}$, the majorization function
$g(x_{j}|{\boldsymbol{x}}^{k})$ is given by:

$\begin{array}[]{ll}g(x_{j}|{\boldsymbol{x}}^{k})=\displaystyle\sum_{i=1}^{m}\displaystyle\sum_{j=1}^{n+1}\dfrac{1}{n+1}\left|\left(n+1\right)c_{ij}\left(\tilde{x}_{j}-\tilde{x}_{j}^{k}\right)+{\boldsymbol{c}}_{i}^{T}{\tilde{{\boldsymbol{x}}}^{k}}\right|^{p}\end{array}$
(18)

where $c_{ij}$ represents the $j^{th}$ element of the ${i}^{th}$ vector
${\boldsymbol{c}}_{i}$ . The above surrogate function can be rewritten as:

$\begin{array}[]{ll}g(x_{j}|{\boldsymbol{x}}^{k})=\displaystyle\sum_{i=1}^{m}\displaystyle\sum_{j=1}^{n}\left|c_{ij}\tilde{x}_{j}+d^{k}_{ij}\right|^{p}+\\\
\displaystyle\sum_{i=1}^{m}\left|c_{i,n+1}\tilde{x}_{n+1}+d^{k}_{i,n+1}\right|^{p}=\displaystyle\sum_{i=1}^{m}\displaystyle\sum_{j=1}^{n}\left|-a_{ij}{x}_{j}+d^{k}_{ij}\right|^{p}+\\\
\hskip
42.67912pt\displaystyle\sum_{i=1}^{m}\left|y_{i}+d^{k}_{i,n+1}\right|^{p}\end{array}$
(19)

where
$d^{k}_{ij}=-c_{ij}\tilde{x}^{k}_{j}+\dfrac{1}{n+1}{\boldsymbol{c}}_{i}^{T}{\boldsymbol{x}}^{k}$.
From above, it can be observed that if $a_{ij}$ is zero then the product
$a_{ij}x_{j}$ becomes equal to zero, hence the first term of
$g(x_{j}|{\boldsymbol{x}}^{k})$ has to be computed only over the non-zero
elements of $a_{ij}$ i.e.,

$\begin{array}[]{ll}g(x_{j}|{\boldsymbol{x}}^{k})=\displaystyle\mathop{\sum_{i=1}^{m}\sum_{j=1}^{n}}_{a_{ij}\neq
0}\left|-a_{ij}{x}_{j}+d^{k}_{ij}\right|^{p}+\displaystyle\sum_{i=1}^{m}\left|y_{i}+d^{k}_{i,n+1}\right|^{p}\end{array}$
(20)

Ignoring the constant terms in $g(x_{j}|{\boldsymbol{x}}^{k})$, the surrogate
minimization problem at any iteration, given ${\boldsymbol{x}}^{k}$ becomes:

$\begin{array}[]{ll}\underset{{\boldsymbol{x}}}{\rm
arg\>min}\displaystyle\mathop{\sum_{i=1}^{m}\sum_{j=1}^{n}}_{a_{ij}\neq
0}w_{ij}\left|{x}_{j}-\dfrac{d^{k}_{ij}}{a_{ij}}\right|^{p}\end{array}$ (21)

where $w_{ij}=|a_{ij}|^{p}$. The above problem is separable in each element of
${\boldsymbol{x}}$. Hence, at every iteration we solve the following surrogate
minimization problem parallely over each element of ${\boldsymbol{x}}$:

$\begin{array}[]{ll}\underset{x_{j}}{\rm
arg\>min}\displaystyle\mathop{\sum_{i=1}^{m}}_{a_{ij}\neq
0}w_{ij}\left|{x}_{j}-\dfrac{d^{k}_{ij}}{a_{ij}}\right|^{p}\end{array}$ (22)

Therefore, each element of ${\boldsymbol{x}}$ can be updated parallely. The
problem in (22) does not have a closed-form solution. However, we can employ a
parameter free bisection method to solve the above problem. The initial
interval $[a,b]$ of the bisection method for the problem in (22) is found
using the following lemma:

###### Lemma IV.2

Consider the following problem wherein given the positive weights
$(w_{1},w_{2}\cdots w_{N})$ and the data points $(a_{1},a_{2}\cdots a_{N})$,
the problem is to estimate $x$ such that it minimizes the weighted sum of
absolute deviations raised to the $p^{th}$ power, i.e.

$\begin{array}[]{ll}\underset{x}{\rm
arg\>min}\>f(x)=\displaystyle\sum_{j=1}^{N}w_{j}|x-a_{j}|^{p}\end{array}$ (23)

The solution $x$ of the above problem will always lie in-between
$\rm{min}({\boldsymbol{a}})$ and $\rm{max}({\boldsymbol{a}})$ where
${\boldsymbol{a}}=[a_{1},a_{2}\cdots a_{N}]^{T}$.

###### Proof:

The gradient of the problem in (23) is given as:

$\begin{array}[]{ll}f^{\prime}(x)=\displaystyle\sum_{j=1}^{N}p\>w_{j}|x-a_{j}|^{p-2}(x-a_{j})=\displaystyle\sum_{j=1}^{N}v_{j}(x-a_{j})\end{array}$
(24)

where $v_{j}=p\>w_{j}|x-a_{j}|^{p-2}$ and is always positive. Note that for
$x=\rm{min}({\boldsymbol{a}})$ and $x=\rm{max}({\boldsymbol{a}})$, the
gradient $f^{\prime}(x)$ will always be lesser than and greater than zero,
respectively. Then, according to the Intermediate Value Theorem [26], the
gradient will be equal to zero for $x$ in the interval
[$\rm{min}({\boldsymbol{a}})$, $\rm{max}({\boldsymbol{a}})$]. ∎

Using lemma IV.2, the initial interval $[a,b]$ for the bisection search
approach to solve the problem in (22) is
$[\rm{min}({\boldsymbol{s}}_{j})\>\>\rm{max}({\boldsymbol{s}}_{j})]^{T}$ where
the $i^{th}$ element of ${\boldsymbol{s}}_{j}$ is given as
$s_{ij}=\dfrac{d^{k}_{ij}}{a_{ij}}$. Further, by the principle of MM, it is
sufficient that the bisection search approach finds a value
${\boldsymbol{x}}^{k}$ at the $k^{th}$ iteration such that
$f_{{}_{\textrm{LR}}}({{\boldsymbol{x}}}^{k})<f_{{}_{\textrm{LR}}}({{\boldsymbol{x}}}^{k-1})$
and not find the exact minimum of the surrogate minimization problem in (22).
Pseudocode of the proposed algorithm is shown in Table 2. We would like to
point out here that unlike the IRLS algorithm, the proposed algorithm does not
have the problem of getting stuck - since the minimum in (22) is found over
the non-zero elements of $a_{ij}$’s. Also, when compared to IRLS and the
modified IRLS algorithm, the proposed algorithm does not involve computing
matrix inverse.

Table 2: Pseudocode of the proposed algorithm for $\ell_{p}$ norm linear
regression problem for $p\in(1,\infty)$.
---
Input: Noisy observed data ${\mathbf{y}}$, data matrix ${\mathbf{A}}$ and $p$
Initialize: Set $k=0$. Initialize ${{{\boldsymbol{x}}}^{0}}$
Precompute: $w_{ij}=|a_{ij}|^{p}$, $i\in(1,\cdots m)$ and $j\in(1,\cdots n)$
Repeat:
Compute the following parallely over all the elements of ${\boldsymbol{x}}$:
1) Compute
$s_{ij}=\dfrac{d^{k}_{ij}}{a_{ij}}=x^{k}_{j}+\dfrac{y_{i}-{\boldsymbol{a}}_{i}^{T}{\boldsymbol{x}}^{k}}{(n+1)a_{ij}}$,
$i\in(1,2,\cdots m)$.
2) $x_{j}^{k+1}$ is obtained by solving (22) using bisection method with
$a=\rm{min}\left({\boldsymbol{s}}_{j}\right)$ and
$b=\rm{max}\left({\boldsymbol{s}}_{j}\right)$.
$k\leftarrow k+1$, until
$\left|\dfrac{{f_{{}_{\textrm{LR}}}({{{{\boldsymbol{x}}}^{k}}})}-{f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}^{k-1}})}}{{f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}^{k-1}})}}\right|\leq
10^{-3}$

Before ending this subsection, we will here discuss the computational
complexity of the proposed algorithm. First, we would like to point out that
the weights $w_{ij}$ which are required to compute $x_{j}^{k+1}$ are
independent of the optimization variable and hence can be precomputed with a
complexity of $\mathcal{O}(p)$. Next, at every iteration, the complexity of
the algorithm is only dictated by the computation of the inner product
${\boldsymbol{a}}_{i}^{T}{\boldsymbol{x}}^{k}$ and the bisection search step.
The complexity of the inner product
${\boldsymbol{a}}_{i}^{T}{\boldsymbol{x}}^{k}$ is $\mathcal{O}(n)$ and the
complexity of the bisection search algorithm is
$\mathcal{O}({\textrm{log}}\>l)$, where $l$ is the number of sub-intervals.
Hence, the total complexity of the proposed algorithm is
$\mathcal{O}\left(p+k\left(n+{\textrm{log}}\>l\right)\right)$. Also as
mentioned earlier, the complexity of the state-of-the-art IRLS algorithm is
dictated by the computation of the inverse of the matrix
${\mathbf{A}}^{T}{\boldsymbol{W}}^{k}{\mathbf{A}}$ and hence the complexity of
the IRLS algorithm is $\mathcal{O}(n^{3})$. Hence, the proposed algorithm has
the least computational complexity when compared to the state-of-the-art
algorithm.

###  IV-B Special Cases

When compared to the $\ell_{p}$ norm regression problem for $p\in(1,\infty)$,
the $\ell_{1}$ and $\ell_{\infty}$ norm regression problems are non-
differentiable - thereby one cannot use the bisection search algorithm
developed in the previous subsection to solve the $\ell_{1}$ and
$\ell_{\infty}$ norm regression problem. Hence, in this subsection, we discuss
an alternate approach to solve the surrogate minimization problem and extend
PROMPT algorithm to solve the $\ell_{1}$ and $\ell_{\infty}$ norm regression
problems.
1\. $\ell_{1}$ norm linear regression - Consider the following $\ell_{1}$ norm
linear regression problem:

$\begin{array}[]{ll}f_{{}_{LR}}({\boldsymbol{x}})=\|{\mathbf{y}}-{\mathbf{A}}{\boldsymbol{x}}\|_{1}=\displaystyle\sum_{i=1}^{m}|y_{i}-{\boldsymbol{a}}_{i}^{T}{\boldsymbol{x}}|\end{array}$
(25)

Similar to the $\ell_{p}$ norm regression problem ($p>2$) in (14), we rewrite
the above problem as:

$\begin{array}[]{ll}f_{{}_{LR}}({\boldsymbol{x}})=\displaystyle\sum_{i=1}^{m}|{\boldsymbol{c}}_{i}^{T}\tilde{{\boldsymbol{x}}}|\end{array}$
(26)

where ${\boldsymbol{c}}_{i}=[c_{i,1}\cdots
c_{i,n+1}]^{T}=[-{\boldsymbol{a}}_{i}\>\>y_{i}]^{T}$ and let
$\tilde{{\boldsymbol{x}}}=[{\boldsymbol{x}},1]^{T}$. Then by using lemma IV.1,
we majorize the term $|{\boldsymbol{c}}_{i}^{T}\tilde{{\boldsymbol{x}}}|$ and
obtain the following surrogate minimization problem:

$\begin{array}[]{ll}\underset{x_{j}}{\rm
arg\>min}\displaystyle\mathop{\sum_{i=1}^{m}\sum_{j=1}^{n}}w_{ij}\left|{x}_{j}-s_{ij}\right|\end{array}$
(27)

where $w_{ij}=|a_{ij}|$, $s_{ij}=\dfrac{d^{k}_{ij}}{a_{ij}}$ and
$d^{k}_{ij}=-c_{ij}\tilde{x}^{k}_{j}+\dfrac{1}{n+1}{\boldsymbol{c}}_{i}^{T}{\boldsymbol{x}}^{k}$.
Similar to the surrogate mimization problem in (22), the above problem is
separable in each element of ${\boldsymbol{x}}$ and hence can be updated
parallely. In addition, unlike the surrogate minimization problem in (22), the
above surrogate minimization problem has a closed-form solution and is given
by the weighted median of the points $s_{ij}$ [27]. Given the points
$s_{1j},s_{2j},\cdots,s_{mj}$ and its associated weights $w_{ij}$, the
weighted median $x_{j}$ is calculated as
a. Sort the points $s_{ij}$ in ascending order.
b. Normalize the weights $w_{ij}$ associated with the data points $s_{ij}$
such that the sum of the weights is equal to one.
c. The weighted median $x_{j}$ is one of the points $s_{qj}$ satisfying
$\displaystyle\sum_{i=1}^{q-1}w_{ij}\leq\dfrac{1}{2}$ and
$\displaystyle\sum_{i=q+1}^{m}w_{ij}\leq\dfrac{1}{2}$.
The pseudocode of the proposed algorithm for $\ell_{1}$ norm regression
problem is shown in Table 3.

Table 3: Pseudocode of the proposed algorithm for $\ell_{1}$ norm linear
regression problem
---
Input: Noisy observed data ${\mathbf{y}}$, data matrix ${\mathbf{A}}$
Initialize: Set $k=0$. Initialize ${{{\boldsymbol{x}}}^{0}}$
Precompute: $w_{ij}=|a_{ij}|$, $i\in(1,\cdots m)$ and $j\in(1,\cdots n)$
Repeat:
Compute the following parallely over all the elements of ${\boldsymbol{x}}$:
1) Compute
$s_{ij}=\dfrac{d^{k}_{ij}}{a_{ij}}=x^{k}_{j}+\dfrac{y_{i}-{\boldsymbol{a}}_{i}^{T}{\boldsymbol{x}}^{k}}{(n+1)a_{ij}}$,
$i\in(1,2,\cdots m)$.
Table 3: Pseudocode of the proposed algorithm for $\ell_{1}$ norm linear
regression problem
---
2) $x_{j}^{k+1}$ is obtained by finding the weighted median of the points
$s_{ij}$ as discussed in Subsection IV.B.
$k\leftarrow k+1$, until
$\left|\dfrac{{f_{{}_{\textrm{LR}}}({{{{\boldsymbol{x}}}^{k}}})}-{f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}^{k-1}})}}{{f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}^{k-1}})}}\right|\leq
10^{-3}$

2\. $\ell_{\infty}$ norm linear regression - The $\ell_{\infty}$ norm linear
regression problem is given as:

$\begin{array}[]{ll}\underset{{\boldsymbol{x}}}{\rm
arg\>min}\\{f_{{}_{LR}}({\boldsymbol{x}})=\|{\mathbf{y}}-{\mathbf{A}}{\boldsymbol{x}}\|_{\infty}=\\\
\hskip
71.13188pt\displaystyle\max_{1,2,\cdots,m}|y_{i}-{\boldsymbol{a}}_{i}^{T}{\boldsymbol{x}}|\\}\end{array}$
(28)

Similar to the development of the PROMPT algorithm, we majorize the objective
function in (28) using lemma IV.1 and obtain the following surrogate
minimization problem:

$\begin{array}[]{ll}\underset{x_{j}}{\rm
min}\displaystyle\max_{1,2,\cdots,m}w_{ij}\left|{x}_{j}-s_{ij}\right|\end{array}$
(29)

Note that the above surrogate minimization problem is separable in each
element of ${\boldsymbol{x}}$ and therefore we can solve the above problem
parallely over each element of ${\boldsymbol{x}}$: To solve the above problem,
we first rewrite the problem in the epigraph form as:

$\begin{array}[]{ll}\underset{x_{j},z}{\rm min}\>\>z\\\ \textrm{subject
to}\>\>w_{ij}^{2}\left|{x}_{j}-s_{ij}\right|^{2}\leq
z\>\>i=1,2,\cdots,m\end{array}$ (30)

The Lagrange of the above problem is:

$\begin{array}[]{ll}L(x_{j},z,\lambda)=z+\displaystyle\sum_{i=1}^{m}\lambda_{i}\left(w_{ij}^{2}\left|{x}_{j}-s_{ij}\right|^{2}-z\right)\end{array}$
(31)

where $\lambda_{i}$’s are the Lagrange multiplier associated with each
inequality. Next, we derive the KKT conditions for the problem in (30).

  1. 1.

Minimizing the Lagrange function with respect to $z$ we get:

$\begin{array}[]{ll}\displaystyle\sum_{i=1}^{m}\lambda_{i}=1\end{array}$ (32)

  2. 2.

Minimizing the Lagrange function with respect to $x_{j}$ we get:

$\begin{array}[]{ll}x_{j}=\dfrac{\displaystyle\sum_{i=1}^{m}\lambda_{i}w_{ij}^{2}s_{ij}}{\displaystyle\sum_{i=1}^{m}\lambda_{i}w_{ij}^{2}}\end{array}$
(33)

  3. 3.

By the complementary slackness property we have:

$\begin{array}[]{ll}\lambda_{i}\left(z-w_{ij}^{2}|x_{j}-s_{ij}|^{2}\right)=0\end{array}$
(34)

  4. 4.

Each dual variable $\lambda_{i}$ must be non-negative. Since the problem in
(30) is convex, any pair of $(x_{j},\lambda_{i})$ that satisfies the KKT
conditions would be the primal and dual optimal solutions. Hence, we use the
above KKT conditions to solve the problem in (30) iteratively [28]. The pseudo
code of the proposed algorithm to solve the $\ell_{\infty}$ norm linear
regression problem is shown in Table. 3.

Table 3: Pseudocode of the proposed algorithm for $\ell_{\infty}$ norm linear
regression problem
---
Input: Noisy observed data ${\mathbf{y}}$, data matrix ${\mathbf{A}}$
Initialize: Set $k=0$. Initialize ${{{\boldsymbol{x}}}^{0}}$ and
${{\lambda}_{i}^{0}}$.
Precompute: $w_{ij}=|a_{ij}|$, $i\in(1,\cdots m)$ and $j\in(1,\cdots n)$
Repeat:
Compute the following parallely over all the elements of ${\boldsymbol{x}}$:
1) Compute
$s_{ij}=\dfrac{d^{k}_{ij}}{a_{ij}}=x^{k}_{j}+\dfrac{y_{i}-{\boldsymbol{a}}_{i}^{T}{\boldsymbol{x}}^{k}}{(n+1)a_{ij}}$,
$i\in(1,2,\cdots m)$.
2) Using the KKT conditions solve the problem in (30) iteratively to obtain
$x_{j}^{k+1}$
$k\leftarrow k+1$, until
$\left|\dfrac{{f_{{}_{\textrm{LR}}}({{{{\boldsymbol{x}}}^{k}}})}-{f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}^{k-1}})}}{{f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}^{k-1}})}}\right|\leq
10^{-3}$

###  IV-C Proof of convergence

Given that the proposed algorithm is based on MM procedure, it is guaranteed
that the sequence of points $\\{{{\boldsymbol{x}}}^{k}\\}$ generated by MM
algorithm will monotonically decrease the problem in (4). Moreover, since
$f_{{}_{\textrm{LR}}}({{\boldsymbol{x}}})$ in (4) is bounded below by zero, it
is ensured that the sequence $f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}})$ will
converge to a finite value.
We now show that the sequence $\\{{{\boldsymbol{x}}}^{k}\\}$ converges to the
stationary point of the problem in (4). From the monotonic property of MM we
have:

$\begin{array}[]{ll}f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}^{0}})\geq
f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}^{1}})\geq
f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}^{2}})\end{array}$ (35)

Assume that there is a subsequence ${{\boldsymbol{x}}}^{r_{j}}$ converging to
a limit point $\tilde{{\boldsymbol{q}}}$. Then from (11), (12) and (35) we
get:

$\begin{array}[]{ll}g({{\boldsymbol{x}}}^{r_{j+1}}|{{\boldsymbol{x}}}^{r_{j+1}})=f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}^{r_{j+1}}})\leq
f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}^{r_{j}+1}})\leq
g({{\boldsymbol{x}}}^{r_{j}+1}|{{\boldsymbol{x}}}^{r_{j}})\\\ \hskip
71.13188pt\leq g({{\boldsymbol{x}}}|{{\boldsymbol{x}}}^{r_{j}})\end{array}$
(36)

where $g(\cdot)$ is the surrogate function as defined in (20). Then, letting
$j\rightarrow\infty$, we get:

$\begin{array}[]{ll}g(\tilde{{\boldsymbol{q}}}|\tilde{{\boldsymbol{q}}})\leq
g({{\boldsymbol{x}}}|\tilde{{\boldsymbol{q}}})\end{array}$ (37)

which implies
$g^{\prime}(\tilde{{\boldsymbol{q}}}|\tilde{{\boldsymbol{q}}})\geq 0$. Since
the first order behavior of surrogate function is same as function
$f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}})$ ([12]),
$g^{\prime}(\tilde{{\boldsymbol{q}}}|\tilde{{\boldsymbol{q}}})\geq 0$ implies
$f^{\prime}_{{}_{\textrm{LR}}}({\tilde{{\boldsymbol{q}}}})\geq 0$. Hence,
$\tilde{{\boldsymbol{q}}}$ is the stationary point of
$f_{{}_{\textrm{LR}}}({{{\boldsymbol{x}}}})$ and therefore the proposed
algorithm converges to the stationary point of the problem in (4).

###  IV-D Graph Based semi-supervised learning

In this subsection we discuss the application of proposed algorithm for the
graph based semi-supervised learning problem. Graph based semi-supervised
learning consists of two steps:

  1. 1.

Construction of graph $G=(V,E)$ from the given data.

  2. 2.

Using the constructed graph estimating the labels of the unlabeled data points
using an appropriate algorithm.

We now discuss each step in detail. In the case of first step, as suggested in
[12], one can construct a symmetric $K$-nearest neighbour (K-NN) graph from
the given data. In K-NN graph, each vertex are used to represent a data point.
An edge between the data points ${\boldsymbol{v}}_{i}$ and
${\boldsymbol{v}}_{j}$ exists if and only if
${\boldsymbol{v}}_{i}\in{\textrm{KNN}}({\boldsymbol{v}}_{j})$ or
${\boldsymbol{v}}_{j}\in{\textrm{KNN}}({\boldsymbol{v}}_{i})$, where
${\textrm{KNN}}({\boldsymbol{v}}_{j})$ stands for the $K$-nearest neighbours
of the ${\boldsymbol{v}}_{j}$ among the data points
${\boldsymbol{v}}_{1},\cdots,{\boldsymbol{v}}_{j-1},{\boldsymbol{v}}_{j+1},\cdots,{\boldsymbol{v}}_{n}$.
The $K$-nearest neighbours of ${\boldsymbol{v}}_{j}$ are found by measuring
the Euclidean distance between ${\boldsymbol{v}}_{j}$ and each of the
remaining data points
${\boldsymbol{v}}_{1},\cdots,{\boldsymbol{v}}_{j-1},{\boldsymbol{v}}_{j+1},\cdots,{\boldsymbol{v}}_{n}$
and choosing the data points corresponding to the $K$ smallest distance. Each
edge is associated with a non-negative edge weight which is defined as:

$\displaystyle w_{ij}$ $\displaystyle=$
$\displaystyle\begin{cases}\textrm{exp}\left(-\dfrac{\|{\boldsymbol{v}}_{i}-{\boldsymbol{v}}_{j}\|_{2}^{2}}{\sigma^{2}}\right)&;{\boldsymbol{v}}_{i}\in{\textrm{KNN}}({\boldsymbol{v}}_{j})\\\
\textrm{exp}\left(-\dfrac{\|{\boldsymbol{v}}_{j}-{\boldsymbol{v}}_{i}\|_{2}^{2}}{\sigma^{2}}\right)&;{\boldsymbol{v}}_{j}\in{\textrm{KNN}}({\boldsymbol{v}}_{i})\\\
0&;\textrm{otherwise}\end{cases}$ (38)

where
$\sigma=\dfrac{1}{2}{\rm{max}}\\{\|{\boldsymbol{v}}_{i}-{\boldsymbol{v}}_{j}\|_{2}^{2}:w_{ij}>0\\}$.
From the constructed graph $G$, we now estimate the labels of the unlabeled
points by solving the problem in (6). For the ease of readability, we re-write
the problem once again here:

$\begin{array}[]{ll}\underset{f}{\rm arg\>min}\>\displaystyle\sum_{i,j\in
E}w_{ij}\left|f(v_{i})-f(v_{j})\right|^{p}\\\ \hskip
8.53581pt{\rm{subject\>to}}\>f(v_{i})=h_{i},i\in L\end{array}$ (39)

Let $f_{i}=f(v_{i})$ and ${\boldsymbol{f}}=[f_{1},f_{2}\cdots
f_{n}]\in\mathbf{R}^{n}$ and
${\boldsymbol{h}}=[h_{1},h_{2},\cdots,h_{l}]\in\mathbf{R}^{l}$. Then the
constrained problem in (39) can be re-written as:

$\begin{array}[]{ll}\underset{{\boldsymbol{f}}}{\rm
arg\>min}\>\displaystyle\sum_{i=1}^{u}\displaystyle\sum_{j=1}^{u}w_{ij}\left|f_{i}-f_{j}\right|^{p}+\displaystyle\sum_{i=1}^{u}\displaystyle\sum_{j=1}^{l}w_{i,j+u}|f_{i}-h_{j}|^{p}\end{array}$
(40)

We now convert the above problem to $\ell_{p}$ norm regression problem. To do
so, we first define the edge vertex incidence matrix ${\boldsymbol{B}}$. Each
row of the edge incidence matrix ${\boldsymbol{B}}$ corresponds to an edge and
each column corresponds to a vertex and its entries are given as:

$\displaystyle{\boldsymbol{B}}_{ev}$ $\displaystyle=$
$\displaystyle\begin{cases}1&;{\textrm{if v is the head of e}}\\\
-1&;{\textrm{if v is the tail of e}}\\\ 0&;\textrm{otherwise}\end{cases}$ (41)

Let ${\boldsymbol{G}}={\boldsymbol{W}}^{1/p}{\boldsymbol{B}}$, where
${\boldsymbol{W}}$ is a diagonal matrix with the weights of the edges and
${\boldsymbol{z}}=-[{\boldsymbol{B}}]_{:,u+1:u+l}{\boldsymbol{h}}$, where the
notation $[{\boldsymbol{B}}]_{:,u+1:u+l}$ is used to represent all the rows of
$u+1$ to $u+l$ columns of the matrix ${\boldsymbol{B}}$. Then the problem in
(40) can be re-written as:

$\begin{array}[]{ll}\underset{{\boldsymbol{f}}}{\rm
arg\>min}\|{\boldsymbol{G}}{\boldsymbol{f}}-{\boldsymbol{z}}\|_{p}\end{array}$
(42)

which is an $\ell_{p}$-norm regression problem and the proposed algorithm can
be used to solve it.

##  V Numerical simulation

In this section, we compare the performance of the proposed algorithm with the
state-of-the art algorithms used to solve the $\ell_{p}$ norm linear
regression problem. In particular, we compare the proposed algorithm with the
IRLS algorithm ([8, 9, 10, 11]) and the modified IRLS algorithm [15]. Also, we
compare the performance of the proposed algorithm with the modified IRLS
algorithm for the Graph based semi-supervised learning problem. The algorithms
were implemented in MATLAB. All the simulations were carried out on a PC with
$2.40$ GHz Intel Xeon Processor with $64$ GB RAM.

1\. In this simulation, we fix $p$ equal to $5$, the dimension of the
observation vector ${\mathbf{y}}\in\mathbb{R}^{m}$ as $50$, the dimension of
${\boldsymbol{x}}\in\mathbb{R}^{n}$ as $20$ and compare the convergence of the
proposed algorithm with the state-of-the art algorithms, IRLS algorithm ([8,
9, 10, 11]) and the modified IRLS algorithm [15]. The elements of the matrix
${\mathbf{A}}$ and the observed data ${\mathbf{y}}$ were randomly generated
from Normal distribution with zero mean and unit variance. All the algorithms
were initialized at
${{\boldsymbol{x}}^{0}}=\left({\mathbf{A}}^{T}{\mathbf{A}}\right)^{-1}\left({\mathbf{A}}^{T}{\mathbf{y}}\right)$,
i.e., at the optimal solution of $\ell_{2}$ norm regression problem. Fig. 2.a.
shows the objective value vs run time of the proposed algorithm and the
modified IRLS algorithm and Fig. 2.b. shows the objective value vs run time of
the IRLS algorithm. The ground truth shown in Fig. 2.a. is obtained using CVX
[29]. From Fig. 2 it can be seen that the proposed algorithm and the modified
IRLS algorithm converges to the true solution, while the IRLS algorithm fails
to converge. Also, from Fig. 2.a. it can be seen that the proposed algorithm
converges faster than the modified IRLS algorithm which could be because the
proposed algorithm can parallely update each elements of ${\boldsymbol{x}}$
when compared to the non-parallel modified IRLS algorithm. Although we show
the convergence of the algorithms for single run in Fig. 2, we observed
similar behavior for many runs.

(a) Objective value vs run time of the proposed algorithm and the modified
IRLS algorithm

(b) Objective value vs run time of the IRLS algorithm

Figure 2: Comparison of convergence of the proposed algorithm with the
modified IRLS and the IRLS algorithm with ${\boldsymbol{x}}^{0}$ equal to the
optimal point of the $\ell_{2}$ norm regression problem.

2\. In this simulation, we compare the convergence of the algorithms with the
elements of ${\boldsymbol{x}}^{0}$ generated randomly from Normal distribution
with mean zero and unit variance. Similar to the previous experiment, we fix
$p$ equal to $5$, the dimension of the observation vector ${\mathbf{y}}$ as
$50$ and the dimension of ${\boldsymbol{x}}$ equal to $20$ and generate the
elements of the data matrix ${\mathbf{A}}$ and the observed data
${\mathbf{y}}$ randomly from Normal distribution with zero mean and unit
variance. Fig. 3.a. shows the objective value vs run time of the proposed
algorithm and the modified IRLS algorithm for $10$ different random
initializations ${\boldsymbol{x}}^{0}$. Fig. 3.b. shows the objective value vs
run time of the IRLS algorithm for one such random initialization. From Fig. 3
it can be seen that the proposed algorithm converges to the ground truth while
the IRLS and the modified IRLS algorithms fails to converge. Hence, the
modified IRLS algorithm converges to the optimal solution of the $\ell_{p}$
norm regression problem only when initialized with the solution of $\ell_{2}$
norm regression problem. This initialization scheme has two issues. Firstly,
it requires the matrix ${\mathbf{A}}^{T}{\mathbf{A}}$ to be invertible and
secondly, the computational complexity associated with the inversion of
${\mathbf{A}}^{T}{\mathbf{A}}$ is $\mathcal{O}(n^{3})$ \- which make this
initialization expensive for large value of $n$.

(a) Objective value vs run time of the proposed and the algorithm in [15]

(b) Objective value vs run time of the IRLS algorithm

Figure 3: Comparison of convergence of the proposed algorithm with the
modified IRLS and the IRLS algorithm with ${\boldsymbol{x}}^{0}$ generated
randomly from Normal distribution

3\. In this simulation, we show that the proposed algorithm converges to the
optimal solution irrespective of the initialization. To show the same we
generate elements of ${\boldsymbol{x}}^{0}$ from three different distributions
- normal distribution with zero mean and unit variance, uniform distribution
from $[0,1]^{n}$ and exponential distribution with mean equal to $0.1$.
Similar to the previous experiments, the dimension of the observation vector
and the dimension of ${\boldsymbol{x}}$ were set equal to $50$ and $20$,
respectively and their elements were randomly generated from Normal
distribution with zero mean and unit variance. The algorithm was made to run
until the following condition was met:

$\begin{array}[]{ll}f_{{}_{LR}}({\boldsymbol{x}}^{k})-f_{{}_{LR}}({\boldsymbol{x}}^{*})\leq
10^{-3}\end{array}$ (43)

where $f_{{}_{LR}}({\boldsymbol{x}}^{*})$ is the optimal value of the problem
in (4) and is obtained by solving the problem using CVX. Fig. 4 shows the
objective value vs run time of the proposed algorithm with the elements of
${\boldsymbol{x}}^{0}$ generated randomly from different distributions and for
different values of $p$. From Fig. 4 it can be seen that the proposed
algorithm converges to the optimal solution irrespective of the initialization
scheme. Hence, unlike the modified IRLS algorithm, the proposed algorithm is
not sensitive to the initialization.

(a) Objective value vs run time of the proposed algorithm for $p=10$

(b) Objective value vs run time of the proposed algorithm for $p=30$

(a) Objective value vs run time of the proposed algorithm for $p=80$ Figure 5:
Objective value vs run time of the proposed algorithm for different random
initialization.

4\. In this simulation, we vary the dimension of the observation vector
${\mathbf{y}}$ and compare the convergence speed of the proposed algorithm and
the modified IRLS algorithm. The dimension of the observation vector was
varied i.e. $m$ was varied from $10000$ to $20000$ in steps of $1000$, $n$ was
fixed at $1000$ and $p$ was kept equal to $10$. The elements of the data
matrix ${\mathbf{A}}$ and the observed data ${\mathbf{y}}$ were generated from
a Normal distribution with mean zero and unit variance. Both the algorithms
were initialized at the optimal point of $\ell_{2}$ norm regression problem.
The run time was averaged over 50 trials. Fig. 6 shows the performance of the
algorithms for the varying dimension of the observation vector ${\mathbf{y}}$
and as can be seen from the figure the proposed algorithm converges about five
times faster than the modified IRLS algorithm.

---
Figure 6: Comparison of average run time of the proposed algorithm with the
modified IRLS algorithm for varying size of the observation vector
${\mathbf{y}}$

5\. In this simulation, we show that the proposed algorithm converges to the
optimal solution for the two special cases discussed in section. IV i.e., for
the $\ell_{1}$ and the $\ell_{\infty}$ norm linear regression problems. In the
case of $\ell_{1}$ norm linear regression, we also compare the convergence
speed of the proposed algorithm with the IRLS algorithm. Since, neither the
IRLS algorithm nor the modified IRLS algorithm can be extended to solve the
$\ell_{\infty}$ norm regression problem, we do not compare the convergence
speed of the proposed algorithm with them. For both the special cases, the
elements of the matrix ${\mathbf{A}}$ and the observed data ${\mathbf{y}}$
were randomly generated from Normal distribution with zero mean and unit
variance. Fig. 7 and Fig. 8 shows the objective value vs. run time of the
proposed algorithm for the $\ell_{1}$ and $\ell_{\infty}$ norm regression
problems for different sizes of the data matrix $A$ and for different random
initialization ${\boldsymbol{x}}^{0}$. From Fig. 7 and Fig. 8 it can be seen
that the proposed algorithm converges to the true solution for all the
initialization. Also, from Fig. 7 it can be seen that the proposed algorithm
converges faster than the IRLS algorithm.

(a) Objective value vs run time for problem dimension $m=800$ and $n=3$

(b) Objective value vs run time for problem dimension $m=1000$ and $n=5$.

Figure 7: Objective value vs run time of the proposed and the IRLS algorithm
for the $\ell_{1}$ norm regression problem

(a) Objective value vs run time for problem dimension $m=500$ and $n=3$

(b) Objective value vs run time for problem dimension $m=1000$ and $n=5$.

Figure 8: Objective value vs run time of the proposed algorithm for the
$\ell_{\infty}$ norm regression problem

6\. In this simulation we compare the performance of the proposed algorithm
and the modified IRLS algorithm for the Graph based semi-supervised problem.
To do so, we simulated the labeled and unlabeled data points randomly from a
uniform distribution from $[0,1]^{d}$ for $d=10$. The labels were also
generated randomly from a uniform distribution from $[0,1]$. A K-NN graph was
constructed with the value of $K=10$ using the code in [30]. The non-negative
edge weights were generated using (38). We compared the convergence of the
algorithms by varying the number of unlabeled data points $u$ from $50$ to
$500$ in steps of $50$ and for fixed number of labeled data points equal to
$10$. Fig. 9 compares the average run time of the algorithms for varying
number of unlabeled points. The run time was average over $50$ monte carlo
trials. From Fig. 9 it can be seen as the number of unlabeled points increases
the run time of the modified IRLS algorithm increases exponentially while that
that of the proposed algorithm increases linearly.

Figure 9: Comparison of average run time of the proposed algorithm with the
modified IRLS algorithm for simulated data.

Next, we compare the performance of the algorithms for the graph based semi-
supervised problem on two different real data sets obtained from [31]. We
describe each data set here:

  1. 1.

IRIS data set - This data set contains fifty samples of three different Iris
flower - Iris Setosa, Iris Versicolour and Iris Virginica. Four features were
measured from each sample - sepal length, sepal width, petal length and the
petal width, all in centimeter. Hence, the dimension $d$ of the data points is
equal to four. The task is to learn a classifier which can can classify the
data points into its corresponding species.

  2. 2.

SPECT data set - It is a multivariate data set containing $267$ samples. Each
sample is described by $23$ features obtained from the Single Proton Emission
Computed Tomography (SPECT) images. Hence, the dimension $d$ of each data
point is equal to $23$. The task is to classify the given data point as normal
or abnormal perfusion.

Similar to the previous experiment we compare the performance of the
algorithms by varying the number of unlabeled data points. The unlabeled data
points were obtained by taking $l_{c}$ data points from each class such that
the total number of labeled data points is equal to $l$ and treating the
remaining data as unlabeled data. Table. I and Table. II compares the run time
of the algorithms for varying number of labeled points per class and from the
Table it can be seen that the proposed algorithm converges faster than the
modified IRLS algorithm for both the data sets.

Table I: Comparison of run time of algorithms in seconds for IRIS data set No. of labeled data points per class |  Proposed algorithm |  Modified IRLS [15]
---|---|---
$2$ |  $0.18$ |  $0.29$
$4$ |  $0.14$ |  $0.28$
$6$ |  $0.11$ |  $0.25$
$8$ |  $0.08$ |  $0.24$
$10$ |  $0.07$ |  $0.21$
Table II: Comparison of run time of algorithms in seconds for SPECT data set No. of labeled data points per class |  Proposed algorithm |  Modified IRLS [15]
---|---|---
$2$ |  $0.19$ |  $1.01$
$4$ |  $0.16$ |  $0.97$
$6$ |  $0.11$ |  $0.91$
$8$ |  $0.09$ |  $0.86$
$10$ |  $0.08$ |  $0.78$

##  VI Conclusion

In this paper, we proposed an iterative algorithm PROMPT based on the MM
procedure to solve the $\ell_{p}$ norm regression problem for any value of
$p$. PROMPT can parallely update each element of ${\boldsymbol{x}}$ which is
helpful to handle the large scale data efficiently. Also, unlike the state-of-
the-art algorithm, the proposed algorithm converges to to the optimal point
irrespective of the initialization scheme. We show through computer
simulations that the proposed algorithm has faster speed of convergence when
compared to the state-of-the-art algorithms and in the end we also evaluate
the performance of the proposed algorithm for the graph based semi-supervised
problem using real and simulated data.

## References

  * [1] C. M. Bishop, _Pattern recognition and machine learning_. springer, 2006.
  * [2] R. J. Hyndman and G. Athanasopoulos, _Forecasting: principles and practice_. OTexts, 2018.
  * [3] M. S. Paolella, _Linear Models and Time-Series Analysis: Regression, ANOVA, ARMA and GARCH_. John Wiley & Sons, 2018.
  * [4] R. T. Rockafellar, S. Uryasev, and M. Zabarankin, “Risk tuning with generalized linear regression,” _Mathematics of Operations Research_ , vol. 33, no. 3, pp. 712–729, 2008.
  * [5] G. Agro, “Maximum likelihood and lp-norm estimators,” _Statistica Applicata_ , vol. 4, no. 2, pp. 171–182, 1992.
  * [6] A. B. Forsythe, “Robust estimation of straight line regression coefficients by minimizing pth power deviations,” _Technometrics_ , vol. 14, no. 1, pp. 159–166, 1972.
  * [7] H. Leon Harter, “Nonuniqueness of least absolute values regression,” _Communications in Statistics-Theory and Methods_ , vol. 6, no. 9, pp. 829–838, 1977.
  * [8] L. Karlovitz, “Construction of nearest points in the $\ell_{p}$, $p$ even, and $\ell_{\infty}$ norms.” _Journal of Approximation Theory_ , vol. 3, no. 2, pp. 123–127, 1970.
  * [9] C. L. Lawson, “Contribution to the theory of linear least maximum approximation,” _Ph. D. dissertation, Univ. Calif._ , 1961.
  * [10] M. R. Osborne, _Finite algorithms in optimization and data analysis_. John Wiley & Sons, Inc., 1985\.
  * [11] I. F. Gorodnitsky and B. D. Rao, “Sparse signal reconstruction from limited data using focuss: A re-weighted minimum norm algorithm,” _IEEE Transactions on signal processing_ , vol. 45, no. 3, pp. 600–616, 1997.
  * [12] M. F. Rios, J. Calder, and G. Lerman, “Algorithms for $\ell_{p}$-based semi-supervised learning on graphs,” _arXiv preprint arXiv:1901.05031_ , 2019\.
  * [13] S. Kahng, “Best $\ell_{p}$ approximation,” _Mathematics of Computation_ , vol. 26, no. 118, pp. 505–508, 1972.
  * [14] H. Ekblom, “Calculation of linear best $\ell_{p}$ approximations,” _BIT Numerical Mathematics_ , vol. 13, no. 3, pp. 292–300, 1973.
  * [15] D. Adil, R. Peng, and S. Sachdeva, “Fast, provably convergent irls algorithm for p-norm linear regression,” in _Advances in Neural Information Processing Systems_ , 2019, pp. 14 166–14 177.
  * [16] A. Subramanya and P. P. Talukdar, “Graph-based semi-supervised learning,” _Synthesis Lectures on Artificial Intelligence and Machine Learning_ , vol. 8, no. 4, pp. 1–125, 2014.
  * [17] M. Guillaumin, J. Verbeek, and C. Schmid, “Multimodal semi-supervised learning for image classification,” in _2010 IEEE Computer society conference on computer vision and pattern recognition_. IEEE, 2010, pp. 902–909.
  * [18] A. Elmoataz, X. Desquesnes, and M. Toutain, “On the game p-laplacian on weighted graphs with applications in image processing and data clustering,” _European Journal of Applied Mathematics_ , vol. 28, no. 6, pp. 922–948, 2017\.
  * [19] O. Chapelle, B. Schölkopf, A. Zien _et al._ , “Semi-supervised learning, vol. 2,” _Cambridge: MIT Press. Cortes, C., & Mohri, M.(2014). Domain adaptation and sample bias correction theory and algorithm for regression. Theoretical Computer Science_, vol. 519, p. 103126, 2006.
  * [20] X. Zhu, Z. Ghahramani, and J. D. Lafferty, “Semi-supervised learning using gaussian fields and harmonic functions,” in _Proceedings of the 20th International conference on Machine learning (ICML-03)_ , 2003, pp. 912–919.
  * [21] A. El Alaoui, X. Cheng, A. Ramdas, M. J. Wainwright, and M. I. Jordan, “Asymptotic behavior of $\backslash$ell_p-based laplacian regularization in semi-supervised learning,” in _Conference on Learning Theory_ , 2016, pp. 879–906.
  * [22] C. S. Burrus, “Iterative reweighted least squares,” _OpenStax CNX. Available online: http://cnx. org/contents/92b90377-2b34-49e4-b26f-7fe572db78a1_ , vol. 12, 2012.
  * [23] Y. Nesterov and A. Nemirovskii, _Interior-point polynomial algorithms in convex programming_. Siam, 1994, vol. 13.
  * [24] Y. Sun, P. Babu, and D. P. Palomar, “Majorization-minimization algorithms in signal processing, communications, and machine learning,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 3, pp. 794–816, 2016.
  * [25] D. R. Hunter and K. Lange, “A tutorial on MM algorithms,” _The American Statistician_ , vol. 58, no. 1, pp. 30–37, 2004.
  * [26] E. W. Weisstein, “Bolzano’s theorem,” 2008.
  * [27] Y. Vardi and C.-H. Zhang, “The multivariate l1-median and associated data depth,” _Proceedings of the National Academy of Sciences_ , vol. 97, no. 4, pp. 1423–1426, 2000.
  * [28] H. Juel, “Minimax location.” [Online]. Available: http://www2.imm.dtu.dk/courses/02721/kap6.pdf
  * [29] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” http://cvxr.com/cvx, Mar. 2014.
  * [30] M. F. Rios, “Laplacian $\ell_{p}$ graph ssl,” 2019. [Online]. Available: https://github.com/mauriciofloresML/Laplacian_Lp_Graph_SSL
  * [31] D. Dua and C. Graff, “UCI machine learning repository,” 2017. [Online]. Available: http://archive.ics.uci.edu/ml


M.M.: E-mail<EMAIL_ADDRESS>
# A practical compact source of heralded single photons for simple detection
LIDAR

Mateusz P. Mrozowski  Department of Physics, University of Strathclyde, John
Anderson Building,
107 Rottenrow East, Glasgow G4 0NG, United Kingdom  John Jeffers  Department
of Physics, University of Strathclyde, John Anderson Building,
107 Rottenrow East, Glasgow G4 0NG, United Kingdom  Jonathan D. Pritchard
Department of Physics, University of Strathclyde, John Anderson Building,
107 Rottenrow East, Glasgow G4 0NG, United Kingdom

###### Abstract

Optical quantum technologies such as quantum sensing, quantum cryptography and
quantum computation all utilize properties of non-classical light, such as
precise photon-number and entangled photon-pair states, to surpass
technologies based on the classical light. A common route for obtaining
heralded single photons is spontaneous four-wave mixing in optical fibers,
allowing for a well-defined spatial mode, for high efficiency integration into
optical fiber networks. These fibers are typically pumped using large,
commercial, pulsed lasers requiring high-power ($\sim$10 W) pump lasers and
are limited to $\sim$MHz repetition rate. Here we propose a cost-efficient,
compact and mobile alternative. Photon pairs at 660 nm and 960 nm will be
created using four-wave mixing in commercial birefringent optical fiber,
pumped using transform limited picosecond pulses with GHz repetition rates
derived from a 785 nm CW laser diode using cavity-enhanced optical frequency
comb generation. The pulses are predicted to have average power of 275 mW, a
peak power of $>$40 W, and predicted photon yield of $>$2000 pairs detected
per second. This design will be later utilized to implement a quantum
illumination scheme based on a coincidence count between idler and signal
photons – instead of joint measurement between signal and idler. This will
allow for quantum advantage over classic LIDAR without the requirement for
maintaining an interferometric stability in free space.

###### keywords:

Quantum information, Quantum optics, Nonlinear optics, Fiber optics, Ultrafast
optics

##  1 INTRODUCTION

Photon pairs are an essential resource, either as factorisable states of pure
photons for applications such as quantum computation with linear optics [1],
and quantum simulations using photons [2], or as entangled pairs used to
enhance timing measurements [3], improve clock synchronization [4], cancel
dispersion in interferometers [5] and for quantum illumination [6, 7, 8, 9,
10].

Whilst entangled light sources enable quantum illumination to achieve a
theoretical detection enhancement over classical illumination[7], developing a
detection scheme able to utilise this regime in practice - without sensitivity
to decoherence or need for an interferometric stability on a length scale of
an optical wavelength, was so far possible only within an optical fiber [11].
Recent proposals have shown that quantum illumination using heralded rather
than entangled photons can also be exploited in LIDAR schemes[12, 13], in
which simple detection based on coincidence measurements using single photon
detector modules has been demonstrated to enable rejection of thermal
background and robustness against spoofing[7].

Two common routes to obtain heralded photons are Spontaneous Parametric Down
Conversion (SPDC) in non-linear crystals [14], and Spontaneous Four-Wave
Mixing (SFWM) in birefringent optical fibers [15]. The latter offers
advantages in the form of a well defined spatial mode for high efficiency
coupling into optical fiber networks developed for quantum information
protocols [16, 17, 18], as well as the lower material cost of optical fibers,
compared to the cost of non-linear photonic crystals. Recent studies have
shown that using SFWM in commercial grade polarisation-maintaining single-mode
optical fibers provides a high purity and wavelength tuneable source of
heralded single photons, adaptable to both visible[15] and telecoms
wavelengths[19].

In typical experiments, photon pairs are obtained by pumping birefringent
fibers using bulky, high power ($\sim 10$ W) pulsed laser systems, which come
with significant price-tag ($\sim\pounds 100$k), and a fixed repetition rate
($\sim 1-100$ MHz). Here we propose an alternative low cost and compact,
heralded photon source, offering a high repetition rate ($\sim$ GHz), and a
tuneable pulse width ($1-10$ ps) that will allow for temporal control of the
generated photons. This source uses an optical frequency comb generator (OFCG)
based on an intra-cavity electro-optical modulator (EOM) to convert a CW pump
into a mode-locked pulse train [20, 21, 22]. A dual-cavity design is adopted
to achieve a high efficiency output [23], overcoming the major limitation of a
single cavity OFCG sources. We estimate a photon pair detection rate of
$>2000$ photon pairs/s following simulation of our pump source, currently
limited by the damage threshold of the intra-cavity EOM.

##  2 Optical Frequency Comb Generation

###  2.1 Single cavity OFCG

A simple OFCG source involves use of an EOM inside an optical cavity driven
resonantly with the cavity free spectral range (FSR) to generate a comb of
phase-coherent sidebands that convert a CW pump laser into a train of short
optical pulses with a repetition rate equal to twice the modulation frequency.
The output electric field can be calculated by modelling the round-trip
propagation as illustrated in Fig. 1(a)

Figure 1: (a) Model of single cavity OFCG. Here we are using an approximation,
in which light passes through the EOM only once per round trip, while the
crystal is assumed to be twice as long and modulation depth is doubled. (b)
Spectral output of the single cavity OFCG. Efficiency of the system
$\eta=0.3\%$. (c) Temporal output of the single cavity OFCG with a FWHM pulse
$\Delta\tau=1.3$ ps as indicated by the red shaded area. The asymmetry is an
artifact of the crystal dispersion.

using the mathematical model presented in Sec. A.1, resulting in an output
field $\vec{E}_{\mathrm{out}}$ equal to

$\vec{E}_{\mathrm{out}}=t_{\mathrm{c}}\hat{\phi}\left(\frac{\phi_{k}}{2}\right)\times
M(\delta)\times
t_{2}\left(\left[\mathbbm{1}-r_{1}r_{2}t_{\mathrm{c}}^{2}\hat{\phi}(\phi_{k})\times
M(2\delta)\right]^{-1}\times t_{1}\vec{E}_{\mathrm{in}}\right),$ (1)

where $\hat{\phi}(\phi_{k})$ and $M(2\delta)$ are the full round trip phase
change operator and full round trip coupling matrix respectively, $\delta$ is
the modulation depth in radians, $\phi_{k}$ is the round-trip phase change of
$k^{\mathrm{th}}$ sideband, $t_{i}$ and $r_{i}$ are transmittance and
reflectivity of the mirror $M_{i}$, $t_{\mathrm{c}}$ is transmission through
the EOM crystal, and $\vec{E}_{\mathrm{in}}$ is the single frequency optical
field represented as a vector with central mode ($k=0$) having an intensity
equal to 1. This provides the electric field represented in the frequency
domain from which the temporal response can be obtained using a Fourier
transform. The frequency spectrum of the output exhibits an exponential decay
in the electric field as a function of sideband order[21], therefore the pulse
shape will be Lorentzian of the form [22]

$E(t)\propto\left[\frac{1}{1+(4f_{m}\delta Ft)^{2}}\right]E_{\mathrm{in}},$
(2)

where $f_{m}$ is the modulation frequency, and $F$ is the cavity finesse. The
pulse width $\Delta\tau_{p}$ at FWHM intensity can then be approximated as[22]

$\Delta\tau_{p}\approx\frac{\sqrt{\sqrt{2}-1}}{2f_{m}\delta F}.$ (3)

In the system consisting of an EOM placed in the cavity, the dispersion comes
from both the dispersive effect of the electro-optic crystal, as well as the
dispersion of the cavity mirrors. The latter can be neglected because it is
small compared to the dispersive effect of the crystal. The full round trip
phase change per sideband is given by [24]

$\phi_{k}=\frac{k^{2}}{2}\mathrm{GVD}\left(2\pi f_{m}\right)^{2}L_{c},$ (4)

where $k$ represents the sideband order, $f_{m}$ is the frequency with which
the EOM is modulated, $\mathrm{GVD}$ is the group velocity dispersion of the
material, and $L_{c}$ is the crystal’s length. The phase difference between
the carrier and the sidebands generated by the EOM increases quadratically
with sideband order. Eventually that difference in phase will be large enough
for the sideband to no longer be resonant with the optical cavity, because of
the destructive interference that light mode will be suppressed, and in effect
the frequency comb will be truncated. This dispersion limit produces a sharp
cut-off in the optical frequency comb spectrum, limiting it’s span to [24]

$\Delta f=\frac{\left[2\left(-\beta+\delta\right)G\right]^{1/2}}{\pi},$ (5)

where $\beta$ is the normalized detuning between the input laser frequency
$f_{o}$ and the nearest cavity resonance $f_{r}$, defined as
$\beta=\left[\pi\left(f_{o}-f_{r}\right)\right]/\mathrm{FSR}$, $\delta$ is
modulation depth, and $G$ is related to the EOM material dispersion by $G=2\pi
c/D\lambda_{o}^{2}L_{c}$, where $L_{c}$ is the crystal length, $\lambda_{o}$
is the input laser wavelength in vacuum and
$D=\left.\frac{\lambda_{o}}{c}\frac{\partial^{2}n}{\partial\lambda^{2}}\right|_{\lambda=\lambda_{o}}$
is the material dispersion which can be related to the group velocity
dispersion by use of [20] $\mathrm{GVD}\cong D\lambda_{o}^{2}/2\pi c$.

Using this model we extract typical performance for a single cavity system
using the following parameters. Mirrors 1 and 2 were set to equal reflectivity
$R_{1}=R_{2}=0.99$, with the EOM modulation depth $\delta=\pi/2$, and
$f_{m}=2.39$ GHz. Due to the light passing twice through the EOM, it
experiences a modulation of $2\delta$ meaning for $\delta>\pi/2$ the light
becomes resonant with the next cavity mode leading to an unstable output rate
and uneven intervals between output pulses[20, 25].

The EOM is modelled as an $L_{c}=20$ mm magnesium-doped lithium tantalate
crystal (MG:LiTaO3) chosen for its linear group velocity dispersion
($\mathrm{GVD}=308$ fs2/mm[26]) at the central wavelength of $\lambda_{0}=785$
nm and assumed to have an anti-reflective (AR) coating with $R=0.01\%$ on each
facet. The output electric field generated in the single cavity can be seen in
Fig. 1(b). As expected, in the frequency domain we can observe an exponential
relationship between sidebands leading to a Lorentzian temporal output pulse,
as shown in Fig. 1(c). The resulting FWHM pulse width $\Delta\tau_{p}=1.3$ ps
and the efficiency of the system $\eta=0.3\%$ at a repetition rate of
$2f_{m}=4.78$ GHz.

###  2.2 A Coupled Cavity OFCG

A major limitation of the single cavity OFCG is the trade-off between
increasing the modulation depth and the cavity finesse to achieve shorter
pulses, and the increased loss from the central cavity mode resulting in an
impedance mis-match that inhibits efficient mode-matching into the cavity and
hence a low output power [20]. For the parameters presented in previous
section the efficiency of the device was predicted to be only $\eta=0.3\%$,
with $<1\%$ typical for experiments [20, 27].
This limitation can be overcome using a coupled cavity configuration as
demonstrated in [23] and shown schematically in Fig. 2(a). Here a planar
mirror is placed between the input mirror M1 and the EOM crystal surface - to
split the OFCG into two separate cavities - an empty coupling cavity (M1 +
M2), and a second cavity (M2 + M3) containing the EOM, referred to as the EOM
cavity. Here the coupling cavity is designed to be almost lossless, with
impedance matched mirrors $R_{1}=R_{2}$ to allow for nearly $100\%$ coupling
into the EOM cavity.
Using the mathematical model derived in Sec. A.2 we obtain the steady state of
the system by solving the coupled equations

$\displaystyle\vec{E}_{1}=\left[\mathbbm{1}-\left(r^{\prime}_{1}t^{\prime}_{2}r_{3}t_{2}t_{\mathrm{c}}^{2}\hat{\phi}_{1}\times\hat{\phi}_{1}\times\hat{\phi}\left(\phi_{k}\right)\times
M(2\delta)\times\left\\{\mathbbm{1}-r^{\prime}_{2}r_{3}t_{\mathrm{c}}^{2}\hat{\phi}\left(\phi_{k}\right)\times
M(2\delta)\right\\}^{-1}+r^{\prime}_{1}r_{2}\hat{\phi}_{1}\hat{\phi}_{1}\right)\right]^{-1}\times
t_{1}\hat{\phi}_{1}\vec{E}_{\mathrm{in}},$ (6a)
$\displaystyle\vec{E}_{\mathrm{out}}=t_{3}t_{\mathrm{c}}\hat{\phi}\left(\frac{\phi_{k}}{2}\right)\times
M(\delta)\left[\mathbbm{1}-r^{\prime}_{2}r_{3}t_{\mathrm{c}}^{2}\hat{\phi}\left(\phi_{k}\right)\times
M(2\delta)\right]^{-1}\times t_{2}\vec{E}_{1},$ (6b)

where $t_{i}$ and $r_{i}$ are the transmitance and reflectivity of the mirror
$M_{i}$, $\hat{\phi_{1}}$ is a phase change operator representing half round
trip phase change between the sidebands in coupling cavity,
$\hat{\phi}\left(\phi_{k}\right)$ is a phase operator representing full round
trip phase change between the sidebands in the EOM cavity.

Figure 2: (a) Model of a coupled cavity OFCG. Again we are simplifying the
maths by assuming the crystal is twice as long and light passes it only once
per round-trip. (b) Spectral output of the coupled cavity OFCG. Incidental
light beam is assumed to have a normalised intensity of 1. Parameters used to
obtain this plot were $R_{1}=R_{2}=0.99$, $R_{3}=0.96$, free spectral range of
the coupling cavity was set to $21.13$ GHz, while the second cavity was set to
match the 5th subharmonic of frequency of modulation $f_{m}=2.39$ GHz.
Modulation depth of the crystal $\delta=\pi/2$, resulting in efficiency of the
system of $90.4\%$. (c) Temporal output from the coupled cavity OFCG with a
FWHM pulse $\Delta\tau_{p}=1.33$ ps as indicated by the red area. Asymmetry is
an artifact of the crystal dispersion and finite coupling matrix.

For the two cavity setup, an additional restriction on the frequency width of
the generated comb arises when the $i^{\mathrm{th}}$ sideband becomes resonant
with the $j^{\mathrm{th}}$ mode of the input coupling cavity. To maximise the
comb-width, the FSR of the input cavity is chosen with a large minimum common
denominator with the FSR of the EOM cavity, and a large finesse to reduce the
cavity linewidth. The first sideband that is completely resonant with the
coupling cavity is called the characteristic mode, which corresponds to the
maxiumum sideband order of the comb.

To enhance photon pair generation using the source, the coupled cavity design
was optimised to maximise the peak pulse energy whilst remaining below the
damage threshold of the EOM. For the EOM crystal detailed in Sec. 2.1 to
achieve $\delta=\pi/2$ at a modulation frequency $f_{m}=2.39$ GHz the crystal
aperture is limited to $2\times 2$ mm, with a damage threshold of 20 W/mm2 set
by the AR coating. To mitigate this restriction the EOM cavity is extended to
set $\mathrm{FSR_{2}}=f_{m}/5$, equal to the 5th subharmonic of the modulation
frequency. The source performance is optimised through adjusting mirror
reflectivities and radius of curvature to define beam waist in the EOM and the
FSR of the input cavity. Possible values of $\mathrm{FSR_{1}}$ were dictated
by the curvatures of the mirrors allowing for desirable beam waist on the
crystal facet. Final value of $\mathrm{FSR_{1}}=21.13$ GHz was chosen due to
the characteristic mode appearing farthest away from the central mode in the
spectrum. The peak power per pulse was calculated from

$P_{\mathrm{pk}}=\frac{\eta P_{\mathrm{in}}}{\Delta\tau_{p}2f_{m}},$ (7)

where $\eta$ is the efficiency of the EOM cycle, $P_{\mathrm{in}}$ is the
maximal average pump power of the system (up to 300 mW) - limited by the
damage threshold of the EOM crystal. Fig. 3(a) shows results of optimisation
of the finesse of the EOM cavity. Fig. 3(b) shows the variation of the
modulation depth in the system running with optimised parameters.

Figure 3: (a) Optimisation plot for the Coupled Cavity OFCG. Finesse of the
EOM cavity was varied whilst $\delta=\pi/2$, $f_{m}=$ 2.39 GHz,
$\mathrm{FSR_{1}}=$ 21.13 GHz, $\mathrm{FSR_{2}}=$ 478 MHz, and
$R_{1}=R_{2}=0.99$ \- the optimal value of the coupling cavity’s finesse ($F$
= 300). Additionally input power, was scaled so the intra-cavity intensity was
below the damage threshold of the EOM crystal. On this plot the relation
between peak power, pulse width and finesse can be observed. Increasing the
finesse allows for shorter pulses but quickly ramps up the intra-cavity
intensity, which places a restriction on $P_{\mathrm{in}}$ and has negative
impact on the peak power. (b) Variation of modulation depth while the system
is operating with $f_{m}=2.39$ GHz and for optimised reflectivities of the
mirrors ($R_{1}=R_{2}=$ 0.99 & $R_{3}=$ 0.96). As we decrease the modulation
depth below $\pi/2$, the FWHM pulses $\Delta\tau_{p}$ become longer,
respectively 1.3, 2.2 and 3.6 ps, which has negative effect on the peak power.

Following the optimisation the coupling cavity is chosen to have
$\mathrm{FSR_{1}}=21.13$ GHz and $R_{1}=R_{2}=0.99$ yielding a finesse of 300.
The EOM cavity has $\mathrm{FSR_{2}}=478$ MHz leading to the characteristic
mode of 1880 and finesse of 120. The simulated output field spectrum and
temporal profile are shown in Fig. 2(b-c) at the optimised peak power, using
system parameters $R_{1}=R_{2}=0.99$, $R_{3}=0.96$, $\delta=\pi/2$,
$f_{m}=2.39$ GHz, FSR of the coupling cavity was 21.13 GHz, while the EOM
cavity was set to match the $5^{\mathrm{th}}$ subharmonic of frequency of
modulation. For these parameters FWHM of the pulse width $\Delta\tau_{p}=$
1.33 ps, $\eta=90.4\%$, the repetition rate equal to 4.78 GHz, peak power
$42.4$ W, and beam waist $\omega_{0}=$ 144 $\mathrm{\mu m}$.

##  3 Photon pair generation

Figure 4: (a) Theoretical birefringent phase-matching contour as a function of
the pump central wavelength $\lambda_{p}$. Signal wavelengths marked in green,
idler wavelengths marked in red and blue dotted line marking the central
frequency of the pump used. Values for the Fibercore HB800G optical fiber,
operating in normal regime $\Delta n=4.3\times 10^{-4}$ [15]. (b) Resulting
joint spectral intensity of the generated signal & idler pair, with purity of
71%, fiber is assumed to be $L=10$ cm long.

For conversion of the output of the coupled cavity OFCG into spectrally pure
and factorisable photon pairs for quantum illumination we propose use of FWM
in a birefringent fiber, where photon pair frequencies are determined by
energy conservation and by phase matching contour

$\kappa=k_{s}+k_{i}-2k_{p}+\frac{2}{3}\gamma P_{p},$ (8)

where, $\kappa$ is the phase-mismatch between the propagation constants of the
signal ($s$), idler ($i$) and two pump waves ($p$) - caused by the chromatic
dispersion in the fiber, with $k_{j}=\frac{n_{j}\omega_{j}}{c}$
$\left(j=s,i,p\right)$ for light wave with angular frequency $\omega_{j}$
propagating in the medium with refractive index $n_{j}$.

For a HB800G fibre, a 785 nm pump wavelength is converted to a $\sim 660$ nm
and $\sim 960$ nm signal and idler pair - shown in Fig. 4(a-b), ideally suited
for demonstration of quantum LIDAR due to working with conventional GaAs SPAD
detectors. The spectral purity of the created pair is determined by the joint
spectral amplitude (JSA), $f(\omega_{s},\omega_{i})$[19]. In the ideal case of
pure state generation the JSA of the output signal and idler photon pair
should be separable - $f(\omega_{s},\omega_{i})=S(\omega_{s})I(\omega_{i})$.
Here functions $S(\omega_{s})$ and $I(\omega_{i})$ represent the spectral
amplitudes of the signal and idler fields.

To quantify the degree of inseparability (the degree of correlation between
signal and idler modes) Schmidt decomposition of the JSA can be utilised[28].
The resulting Schmidt number $K$ represents the number of modes excited in FWM
($K=1$) for a pure state, which is related to purity via $\mathcal{P}=1/K$
[28]. For our system we consider a fiber length 10 cm corresponding to an
estimated purity of 71 % - shown on Fig. 4(b), which can be increased using
shorter fibers at the expense of reduced pair generation rates.

To estimate our expected pair detection rate we scale experimental data from
[15], featuring a commercial pulsed Ti:Sapphire laser generating bandwidth
limited pulses with $\Delta\lambda=3$ nm centred at 704 nm, with a repetition
rate $R_{r}=80$ MHz. For an $L_{c}=$ 10 cm fiber, pair rates of
$N_{\mathrm{exp}}=$ 23000 pairs/s are obtained from an average power of 15 mW.
Assuming $\mathrm{sech^{2}}$ pulse shape, this corresponds to pulse width of
0.17 ps with a peak pulse power of $P_{\mathrm{pk,ref}}=$ 1.1 kW. The pair
production rate is linearly proportional to the fiber length L and
quadratically proportional to peak power squared, $N\propto
LP_{\mathrm{pk}}^{2}$[29]. This allows us to estimate the probability of
creating a photon pair per pulse with $\alpha=N_{\mathrm{exp}}/R_{r},$ which
then can be re-scaled using the data from Sec. 2.2 to obtain the predicted
pair detection rate

$N=\left(\frac{P_{\mathrm{pk}}}{P_{\mathrm{pk,ref}}}\right)^{2}\times\alpha\times
2f_{m},$ (9)

where $P_{\mathrm{pk}}$ is peak power per pulse estimated from our model (42.4
W), and $2f_{m}$ is the repetition rate of our system. We use the above to
predict $>2000$ photon pairs detected per second.

##  4 Conclusion

We have presented a highly efficient coupled cavity OFCG optimised for use in
generating spectrally pure photon pairs for use in quantum illumination
experiments. The source provides a compact, low power and low cost approach to
creating picosecond pulses with high (GHz) repetition rate. We predict using
photon pair generation in a commercial birefriengent single mode optical fiber
to yield $>2000$ photon pairs/s, with photon pair purity of 71%.

Unlike commercial pump sources, this design provides a dynamically tunable
pulse width through controlling the EOM modulation depth. This method can be
adapted to work at other wavelengths and in future extended to telecoms
wavelengths. Exploration of the EOM materials offering higher damage
thresholds will enable higher peak pulse energy and shorter pulse duration to
boost the attainable pair generation rate. The purity of the photon pair
generated can be increased by shortening the fibre to $L=6.81$ cm to achieve
highest pair purity attainable for this source of 80% at the cost of a
reduction in photon pair production. Alternatively, a longer fiber can be
chosen sacrificing purity for the increased pair generation rates.

###### Acknowledgements.

This project is funded by the UK Ministry of Defence.

##  Appendix A Mathematical models

###  A.1 Single Cavity OFCG Mathematical model

The output electric field can be calculated by considering the situation in
Fig. 1(a). We are able to represent all of those electric fields as follows

$\displaystyle\vec{E}_{1}$
$\displaystyle=t_{1}\vec{E}_{\mathrm{in}}+r_{1}\vec{E}_{2},$ (10a)
$\displaystyle\vec{E}_{2}$
$\displaystyle=t_{\mathrm{c}}^{2}\hat{\phi}(\phi_{k})M(2\delta)\times\vec{E}_{1},$
(10b)

where $r_{i}$ and $t_{i}$ stand for reflectivity and transmittance of mirror
Mi and $t_{\mathrm{c}}$ is transmittance trough the cavity. Operator
$\hat{\phi}(\phi_{k})$ is the full round-trip phase change due to crystal
dispersion and can be expressed as

$\hat{\phi}(\phi_{k})=\begin{bmatrix}e^{-i\phi_{k}}&0&\ldots&\ldots&0\\\
0&\ddots&\ddots&&\vdots\\\ \vdots&\ddots&e^{-i\phi_{0}}&\ddots&\vdots\\\
\vdots&&\ddots&\ddots&0\\\ 0&\ldots&\ldots&0&e^{-i\phi_{-k}}\end{bmatrix},$
(11)

where $k$ is the sideband order and $\phi_{k}$ is the phase change due to
crystal dispersion per sideband, given by Eq. (4). $M(2\delta)$ is the
transformation matrix of the EOM. This matrix can be represented as

$M(2\delta)=\begin{bmatrix}J_{0}(2\delta)&-J_{1}(2\delta)&J_{2}(2\delta)&\ldots&J_{2k}(2\delta)\\\
J_{1}(2\delta)&J_{0}(2\delta)&-J_{1}(2\delta)&\ldots&-J_{2k-1}(2\delta)\\\
J_{2}(2\delta)&J_{1}(2\delta)&J_{0}(2\delta)&\ldots&J_{2k-2}(2\delta)\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
J_{2k}(2\delta)&J_{2k-1}(2\delta)&J_{2k-2}(2\delta)&\ldots&J_{0}(2\delta)\end{bmatrix},$
(12)

where $J_{k}(\delta)$ is Bessel function identity of order $i$, and an origin
$\delta$, and where $\delta$ is the modulation depth.

Using Eq. (10a) - (10b) we can solve the field state in $\vec{E}_{1}$
depending on incident field $\vec{E}_{\mathrm{in}}$

$\vec{E}_{1}=\left[\mathbbm{1}-r_{1}r_{2}t_{\mathrm{c}}^{2}\hat{\phi}(\phi_{k})\times
M(2\delta)\right]^{-1}\times t_{1}\vec{E}_{\mathrm{in}},$ (13)

using this solution, we can calculate output electric field of the OFCG -
$\vec{E}_{\mathrm{out}}$, to be

$\vec{E}_{\mathrm{out}}=t_{\mathrm{c}}\hat{\phi}\left(\frac{\phi_{k}}{2}\right)\times
M(\delta)\times(t_{2}\vec{E}_{1}).$ (14)

###  A.2 Coupled Cavity OFCG Mathematical Model

To solve the coupled cavity OFCG we take situation in Fig. 2(a)

$\displaystyle\vec{E}_{1}$
$\displaystyle=\hat{\phi}_{1}\times\left(t_{1}\vec{E}_{\mathrm{in}}+r^{\prime}_{1}\vec{E}_{4}\right),$
(15a) $\displaystyle\vec{E}_{2}$
$\displaystyle=t_{2}\vec{E}_{1}+r^{\prime}_{2}\vec{E}_{3},$ (15b)
$\displaystyle\vec{E}_{3}$
$\displaystyle=t_{\mathrm{c}}^{2}\hat{\phi}_{2}\times M(2\delta)\times
r_{3}\vec{E}_{2},$ (15c) $\displaystyle\vec{E}_{4}$
$\displaystyle=\hat{\phi}_{1}\times\left(t^{\prime}_{2}\vec{E}_{3}+r_{2}\vec{E}_{1}\right),$
(15d)

where $\hat{\phi}_{1}$ is a phase operator, representing half round trip phase
change between sidebands, in first cavity, $\hat{\phi}_{2}$ is a phase
operator, representing full round trip phase change between sidebands in the
EOM cavity, M(2$\delta$) is full round trip transformation matrix, $r(t)_{i}$
represent mirror Mi reflectivity (transmitance) on left side of mirrors
respectively, while $r^{\prime}(t^{\prime})_{i}$ represent mirror Mi
reflectivity (transmitance) on the right side of the mirrors. relation between
those can be represented as

$\displaystyle r^{\prime}$ $\displaystyle=r,$ (16a) $\displaystyle t^{\prime}$
$\displaystyle=-t,$ (16b)

this ensures the transmission phase is preserved. By substituting Eq. (15c)
into Eq. (15b) we obtain

$\displaystyle\vec{E}_{2}$
$\displaystyle=t_{2}\vec{E}_{1}+r^{\prime}_{2}r_{3}t_{\mathrm{c}}^{2}\left(\hat{\phi}_{2}\times
M(2\delta)\times\vec{E}_{2}\right),$ (17a) $\displaystyle\vec{E}_{2}$
$\displaystyle=\left[\mathbbm{1}-r^{\prime}_{2}r_{3}t_{\mathrm{c}}^{2}\hat{\phi}_{2}\times
M(2\delta)\right]^{-1}\times t_{2}\vec{E}_{1},$ (17b)

after that we substitute Eq. (15c) into Eq. (15d) to obtain

$\vec{E}_{4}=\hat{\phi}_{1}\left(t^{\prime}_{2}r_{3}t_{2}t_{\mathrm{c}}^{2}\hat{\phi}_{2}\times
M(2\delta)\left[\mathbbm{1}-r^{\prime}_{2}r_{3}t_{\mathrm{c}}^{2}\hat{\phi}_{2}M(2\delta)\right]^{-1}\times\vec{E}_{1}+r_{2}\vec{E}_{1}\right),$
(18)

finally we substitute Eq. (17b) and Eq. (18) into (15a), to solve the steady
state of the system

$\vec{E}_{1}=\left[\mathbbm{1}-\left(r^{\prime}_{1}t^{\prime}_{2}r_{3}t_{2}t_{\mathrm{c}}^{2}\hat{\phi}_{1}\times\hat{\phi}_{1}\times\hat{\phi}_{2}\times
M(2\delta)\times\left\\{\mathbbm{1}-r^{\prime}_{2}r_{3}t_{\mathrm{c}}^{2}\hat{\phi}_{2}\times
M(2\delta)\right\\}^{-1}+r^{\prime}_{1}r_{2}\hat{\phi}_{1}\hat{\phi}_{1}\right)\right]^{-1}\times
t_{1}\hat{\phi}_{1}\vec{E}_{\mathrm{in}},$ (19)

now to obtain the output electric field we substitute Eq. (19) into Eq. (17b),
and multiply the result by transmittance of the third mirror, half round trip
phase operator in the EOM cavity, and half round trip transformation matrix.

$\vec{E}_{\mathrm{out}}=t_{3}t_{\mathrm{c}}\hat{\phi}\left(\frac{\phi_{k}}{2}\right)\times
M(\delta)\vec{E_{2}}.$ (20)

## References

  * [1] Kok, P., Munro, W. J., Nemoto, K., Ralph, T. C., Dowling, J. P., and Milburn, G. J., “Linear optical quantum computing with photonic qubits”, Rev. Mod. Phys. 79, 135–174 (2007).
  * [2] Loredo, J. C., Broome, M. A., Hilaire, P., Gazzano, O., Sagnes, I., Lemaitre, A., Almeida, M. P., Senellart, P., and White, A. G., “Boson sampling with single-photon fock states from a bright solid-state source”, Phys. Rev. Lett. 118, 130503 (2017).
  * [3] Harris, S. E., “Chirp and compress: Toward single-cycle biphotons”, Phys. Rev. Lett. 98 (2007).
  * [4] V. Giovannetti, S. L. and Maccone, L., “Quantum-enchanced positioning and clock synchronization”, Nature 412, 417–419 (2001).
  * [5] Franson, J. D., “Nonclassical nature of dispersion cancellation and nonlocal interferometry”, Phys. Rev. A 80 (2009).
  * [6] England, D. G., Balaji, B., and Sussman, B. J., “Quantum-enhanced standoff detection using correlated photon pairs”, Phys. Rev. A 99 (Feb 2019).
  * [7] Lloyd, S., “Enhanced sensitivity of photodetection via quantum illumination”, Science 321, 1463–1465 (2008).
  * [8] Tan, S.-H., Erkmen, B. I., Giovannetti, V., Guha, S., Lloyd, S., Maccone, L., Pirandola, S., and Shapiro, J. H., “Quantum illumination with gaussian states”, Phys. Rev. Lett. 101, 253601 (2008).
  * [9] Shapiro, J. H. and Lloyd, S., “Quantum illumination versus coherent-state target detection”, New J. Phys. 11, 063045 (2009).
  * [10] Guha, S. and Erkmen, B. I., “Gaussian-state quantum-illumination receivers for target detection”, Phys. Rev. A 80, 052310 (2009).
  * [11] Zhang, Z., Mouradian, S., Wong, F. N., and Shapiro, J. H., “Entanglement-enhanced sensing in a lossy and noisy environment”, Phys. Rev. Lett. 114 (2015).
  * [12] Frick, S., McMillan, A., and Rarity, J., “Quantum rangefinding”, Opt. Express 28, 37118 (2020).
  * [13] Yang, H., Roga, W., Pritchard, J. D., and Jeffers, J., “Gaussian state-based quantum illumination with simple photodetection”, Opt. Express 29, 8199 (2021).
  * [14] Kaneda, F., Garay-Palmett, K., U’Ren, A. B., and Kwiat, P. G., “Heralded single-photon source utilizing highly nondegenerate, spectrally factorable spontaneous parametric downconversion”, Opt. Express 24, 10733–10747 (2016).
  * [15] Smith, B. J., Mahou, P., Cohen, O., Lundeen, J. S., and Walmsley, I. A., “Photon pair generation in birefringent optical fibers”, Opt. Express 17, 23589–23602 (2009).
  * [16] Clark, A. S., Fulconis, J., Rarity, J. G., Wadsworth, W. J., and O’Brien, J. L., “All-optical-fiber polarization-based quantum logic gate”, Phys. Rev. A 79, 030303 (Mar 2009).
  * [17] Matthews, J. C. F., Politi, A., Stefanov, A., and O’Brien, J. L., “Manipulating multi-photon entanglement in waveguide quantum circuits”, Nat. Photonics 3 (2009).
  * [18] Smith, B. J., Kundys, D., Thomas-Peter, N., Smith, P. G. R., and Walmsley, I. A., “Phase-controlled integrated photonic quantum circuits”, Opt. Express 17 (2009).
  * [19] Lugani, J., Francis-Jones, R. J. A., Boutari, J., and Walmsley, I. A., “Spectrally pure single photons at telecommunications wavelengths using commercial birefringent optical fiber”, Opt. Express 28, 5147 (2020).
  * [20] Kovacich, R. P., Sterr, U., and Telle, H. R., “Short-pulse properties of optical frequency comb generators”, Appl. Opt. 39, 4372–4376 (2000).
  * [21] Kourogi, M., Nakagawa, K., and Ohtsu, M., “Wide-span optical frequency comb generator for accurate optical frequency difference measurement”, IEEE J. Quantum Electron. 29, 2693–2701 (1993).
  * [22] Kobayashi, T., Sueta, T., Cho, Y., and Matsuo, Y., “High‐repetition‐rate optical pulse generator using a fabry‐perot electro‐optic modulator”, Appl. Phys. Lett. 21, 341–343 (1972).
  * [23] Macfarlane, G. M., Bell, A. S., Riis, E., and Ferguson, A. I., “Optical comb generator as an efficient short-pulse source”, Opt. Lett. 21, 534–536 (1996).
  * [24] Kourogi, M., Widiyatomoko, B., Takeuchi, Y., and Ohtsu, M., “Limit of optical-frequency comb generation due to material dispersion”, IEEE J. Quantum Electron. 31, 2120–2126 (1995).
  * [25] Xiao, S., Hollberg, L., Newbury, N. R., and Diddams, S. A., “Toward a low-jitter 10 ghz pulsed source with an optical frequency comb generator”, Opt. Express 16 (2008).
  * [26] Moutzouris, K., Hloupis, G., Stavrakas, I., Triantis, D., and Chou, M.-H., “Temperature-dependent visible to near-infrared optical properties of 8 mol% mg-doped lithium tantalate”, Opt. Mater. Express 1, 458–465 (2011).
  * [27] Xiao, S., Hollberg, L., and Diddams, S. A., “Generation of a 20 ghz train of subpicosecond pulses with a stabilized optical-frequency-comb generator”, Opt. Lett. 34, 85–87 (2009).
  * [28] Zielnicki, K., Garay-Palmett, K., Cruz-Delgado, D., Cruz-Ramirez, H., O’Boyle, M. F., Fang, B., Lorenz, V. O., U’Ren, A. B., and Kwiat, P. G., “Joint spectral characterization of photon-pair sources”, J. Mod. Opt. 65, 1141–1160 (2018).
  * [29] Lin, Q., Yaman, F., and Agrawal, G. P., “Photon-pair generation in optical fibers through four-wave mixing: Role of raman scattering and pump polarization”, Phys. Rev. A 75, 023803 (2007).