# One-Bit Quadratic Compressed Sensing:
From Sample Abundance to Linear Feasibility

Arian Eamaz⋆, Farhang Yeganegi⋆, Deanna Needell†, and Mojtaba Soltanalian⋆
The first two authors have contributed equally to this work. ⋆ECE Department,
University of Illinois Chicago, Chicago, USA
†Department of Mathematics, University of California Los Angeles, Los Angeles,
USA

###### Abstract

One-bit quantization with time-varying sampling thresholds has recently found
significant utilization potential in statistical signal processing
applications due to its relatively low power consumption and low
implementation cost. In addition to such advantages, an attractive feature of
one-bit analog-to-digital converters (ADCs) is their superior sampling rates
as compared to their conventional multi-bit counterparts. This characteristic
endows one-bit signal processing frameworks with what we refer to as _sample
abundance_. On the other hand, many signal recovery and optimization problems
are formulated as (possibly non-convex) quadratic programs with linear
feasibility constraints in the one-bit sampling regime. We demonstrate, with a
particular focus on quadratic compressed sensing, that the sample abundance
paradigm allows for the transformation of such quadratic problems to merely a
linear feasibility problem by forming a large-scale overdetermined linear
system; thus removing the need for costly optimization constraints and
objectives. To efficiently tackle the emerging overdetermined linear
feasibility problem, we further propose an enhanced randomized Kaczmarz
algorithm, called _Block SKM_. Several numerical results are presented to
illustrate the effectiveness of the proposed methodologies.

##  I Introduction

In the past two decades, sparsity-based processing methods have been
attracting a growing interest in statistical signal processing
applications[1]. Quadratic compressed sensing (QCS) is a widely used
formulation in sparse signal recovery; examples include when imaging a sparse
object using partially and spatially incoherent illumination[2], or phase
retrieval for sparse signals [3].

To approach the global optimum, the QCS problem was relaxed as a semidefinite
programming (SDP) problem, which involves minimizing the rank of a lifted
matrix while satisfying both the recovery constraints and the row sparsity
constraints on the signal[1, 4]. To retrieve the sparse solution, an iterative
thresholding algorithm was proposed that leverages a sequence of SDPs. This
approach is similar to the recent developments in the field of phase
retrieval, where similar semidefinite programming-based ideas have been
utilized[5, 4, 6, 7]. Unfortunately, these methods have a high complexity,
making them difficult to use for the QCS problem.

To overcome the computational challenges posed by convex optimization
techniques, non-convex methods have been introduced as an alternative
approach. These methods tackle the phase retrieval problem as a least-square
problem and aim to find a local optimum using various optimization
techniques[3, 8, 9]. In[3], they proposed the greedy sparse phase retrieval
(GESPAR), a fast local search method, to efficiently recover the signal from
measurements of magnitudes in the QCS problem which is more accurate than
existing local methods. However, the highly non-convex and non-unique nature
of the problem presents a challenge in finding an optimal local solution. To
enhance the performance of these local methods, various initialization
algorithms have been proposed to improve their outcomes[10, 11].

Sampling the signals of interest at high data rates with high-resolution ADCs
would dramatically increase the overall implementation cost and power
consumption of the sampling task. In multi-bit sampling scenarios, a very
large number of quantization levels is necessary in order to represent the
original continuous signal with high accuracy, which in practice, leads to a
considerable reduction in sampling rate [12, 13]. This attribute of multi-bit
sampling has served as a key motivator for the proliferation of
underdetermined signal processing tools [6, 14, 15]. An alternative solution
to such challenges is to deploy _one-bit quantization_ , which is an extreme
sampling scenario, where the signals are merely compared with given threshold
levels at the ADC, thus producing sign data ($\pm 1$). This enables signal
processing equipment to sample at a very high rate, with a considerably lower
cost and energy consumption, compared to their conventional counterparts that
employ multi-bit ADCs [16, 12, 17, 18].

The use of a fixed threshold in one-bit quantization can result in
difficulties in accurately estimating the signal amplitude. To address this
issue, recent studies have proposed the use of time-varying thresholds, which
have been shown to enhance signal recovery performance[19, 20, 21, 22].

In this paper, we consider the deployment of one-bit sampling with time-
varying thresholds on QCS, leading to an increased sample size and a _highly
overdetermined system_ as a result. Our proposed method can recover the
desired sparse signal from the _one-bit QCS_ by (i) generating abundant one-
bit measurements, in order to define a large scale overdetermined system where
a finite volume feasible set is created for QCS, and (ii) solving this
obtained linear feasibility problem by leveraging one of the efficient solver
families of overdetermined systems, namely the _Kaczmarz algorithms_.

The Kaczmarz method [23] is an iterative projection algorithm for solving
linear systems of equations and inequalities. It is usually applied to highly
overdetermined systems because of its simplicity. Many variants of this
iterative method and their convergence rates have been proposed and studied in
recent decades for both consistent and inconsistent systems including the
randomized Kaczmarz algorithm, the randomized block Kaczmarz algorithm and
most recently, the sampling Kaczmarz-Motzkin (SKM) method [24, 25, 26, 27].

To reconstruct the signal of interest from the one-bit sampled QCS, we employ
the novel variant of the Kaczmarz algorithm, _Block_ _S_ ampling _K_ aczmarz-
_M_ otzkin (Block SKM) whose theoretical guarantees will be discussed.

_Outline_ : Section II is dedicated to a review of QCS. In Section III, we
will briefly introduce the one-bit sampling via time-varying thresholds and
propose the _one-bit polyhedron_ for the QCS, which is a large-scale
overdetermined system. An accelerated Kaczmarz approach is proposed to find
the optimal point in the one-bit QCS polyhedron in Section IV. Also, the
convergence rate of proposed algorithm is investigated. Section V is devoted
to numerical results of the proposed Kaczmarz algorithm to show its recovery
performance in one-bit. Also, we compare the performance of proposed algorithm
to that of the well-known high-resolution method, GESPAR for the _phase
retrieval_ scenario, when the rank of middle matrix is one. Finally, Section
vi@ concludes the paper.

_Notation:_ We use bold lowercase letters for vectors and bold uppercase
letters for matrices. $\mathbb{C}$ and $\mathbb{R}$ represent the set of
complex and real numbers, respectively. $(\cdot)^{\top}$ and
$(\cdot)^{\mathrm{H}}$ denote the vector/matrix transpose, and the Hermitian
transpose, respectively. $\mathbf{I}_{N}\in\mathbb{R}^{N\times N}$ and
$\mathbf{0}_{N_{1}\times N_{2}}$ are the identity matrix of size $N$ and all-
zero matrix of size $N_{1}\times N_{2}$. $\operatorname{Tr}(.)$ denotes the
trace of the matrix argument. The Frobenius norm of a matrix $\mathbf{B}$ is
defined as
$\|\mathbf{B}\|_{\mathrm{F}}=\sqrt{\sum^{N_{1}}_{r=1}\sum^{N_{2}}_{s=1}\left|b_{rs}\right|^{2}}$
where $\\{b_{rs}\\}$ are elements of $\mathbf{B}$. The $\ell^{0}$-norm of a
vector counts the number of its non-zero elements. The Hadamard (element-wise)
product of two matrices $\mathbf{B}_{1}$ and $\mathbf{B}_{2}$ is denoted as
$\mathbf{B}_{1}\odot\mathbf{B}_{2}$. The vectorized form of a matrix
$\mathbf{B}$ is written as $\operatorname{vec}(\mathbf{B})$. $\mathbf{1}_{s}$
is the $s$-dimensional all-one vector. Given a scalar $x$, we define $(x)^{+}$
as $\max\left\\{x,0\right\\}$. The function $\operatorname{sgn}(\cdot)$ yields
the sign of its argument. The floor operation is denoted by $\lfloor\rfloor$.

##  II Quadratic Compressed Sensing

In QCS, a sparse high-dimensional signal is to be recovered from a quadratic
cost function [1, 3]:

$\displaystyle\min_{\mathbf{x}}$ $\displaystyle\left\|\mathbf{x}\right\|_{0}$
(1) s.t. $\displaystyle
y_{j}=\mathbf{x}^{\mathrm{H}}\mathbf{A}_{j}\mathbf{x},\quad
j\in\mathcal{J}=\left\\{1,\cdots,m\right\\},$

where $\mathbf{x}\in\mathbb{C}^{n}$ is the signal to be recovered,
$\left\\{y_{j}\right\\}$ are the measurements,
$\left\\{\mathbf{A}_{j}\right\\}\in\mathbb{R}^{n\times n}$ are the associated
sensing matrices, and $m$ is the number of measurements. The convex relaxation
of (1) is obtained by the matrix lifting procedure, given by

$\displaystyle
y_{j}=\mathbf{x}^{\mathrm{H}}\mathbf{A}_{j}\mathbf{x}=\operatorname{Tr}\left(\mathbf{A}_{j}\mathbf{xx}^{\mathrm{H}}\right)=\operatorname{Tr}\left(\mathbf{A}_{j}\mathbf{X}\right),$
(2)

where $\mathbf{X}=\mathbf{xx}^{\mathrm{H}}$, and
$\operatorname{Tr}\left(\mathbf{A}_{j}\mathbf{X}\right)=\operatorname{vec}\left(\mathbf{A}^{\top}_{j}\right)^{\top}\operatorname{vec}\left(\mathbf{X}\right)$.
The sparsity constraint on $\mathbf{x}$ can be dealt with by enforcing the
row-sparsity of $\mathbf{X}$. If $\mathbf{x}$ has $k$ non-zero elements, then
$\mathbf{X}$ has $k$ rows containing non-zero elements, and each of these rows
is also $k$-sparse. The row-sparsity of $\mathbf{X}$ may be promoted by adding
a quadratic constraint on $\mathbf{X}$, i.e.,
$\sum_{r}\left(\sum_{s}\left|\mathbf{X}_{rs}\right|^{2}\right)^{\frac{1}{2}}<\eta$,
where $\eta$ is a positive number[2]. Based on (2), the QCS problem can be
reformulated as:

find $\displaystyle\bm{X}$ (3) s.t.
$\displaystyle\operatorname{Tr}\left(\mathbf{A}_{j}\mathbf{X}\right)=y_{j},$
$\displaystyle\sum_{r}\left(\sum_{s}\left|\mathbf{X}_{rs}\right|^{2}\right)^{\frac{1}{2}}<\eta,$
$\displaystyle\operatorname{rank}\left(\mathbf{X}\right)=1,\leavevmode\nobreak\
\mathbf{X}\succeq 0.$

To have a convex program similar to [2], the problem (3) may be relaxed as,

$\displaystyle\min_{\mathbf{X}}$ $\displaystyle\operatorname{Tr}(\mathbf{X})$
(4) s.t.
$\displaystyle\operatorname{Tr}\left(\mathbf{A}_{j}\mathbf{X}\right)=y_{j},$
$\displaystyle\sum_{r}\left(\sum_{s}\left|\mathbf{X}_{rs}\right|^{2}\right)^{\frac{1}{2}}<\eta,\leavevmode\nobreak\
\mathbf{X}\succeq 0.$

The above problem is a semi-definite program (SDP). Similar SDP-based ideas
were recently utilized in the context of phase retrieval [6, 14]. However, the
SDP has a high computational complexity; in particular, the semi-definiteness
and row-sparsity constraint in the above problem render it computationally
demanding [28, 29, 3].

An interesting alternative to enforcing the _feasible set_ in problem (4),
denoted as $\mathcal{F}_{\mathbf{X}}$, emerges when one increases the number
of samples $m$, and solves the overdetermined linear system of equations with
$m\gg n$. In this sample abundance regimen, the linear constraint
$\operatorname{Tr}\left(\mathbf{A}_{j}\mathbf{X}\right)=y_{j}$ may actually
yield the optimum inside $\mathcal{F}_{\mathbf{X}}$. As a result of increasing
the number of samples, it is possible that the intersection of these
hyperplanes will achieve the optimal point without the need to consider other
costly constraints. However, this idea may face practical limitations in the
case of multi-bit quantization systems since ADCs capable of ultra-high rate
sampling are difficult and expensive to produce. Moreover, one cannot
necessarily expect these constraints to intersect with
$\mathcal{F}_{\mathbf{X}}$ in such a way to form a finite-volume space before
the optimum is obtained [15, 6].

In the next section, by deploying the idea of one-bit sampling with time-
varying thresholds, linear equality constraints are superseded by a massive
array of linear inequalities—thus forming a polyhedron that asymptotically
coincides with $\mathcal{F}_{\mathbf{X}}$.

##  III One-Bit QCS

In this section, we will briefly introduce the one-bit sampling with time-
varying scheme and the signal reconstruction problem in the one-bit
quantization scheme. We will demonstrate that the utilization of time-varying
thresholds in one-bit sampling results in a highly over-determined system,
represented as a polyhedron. Subsequently, by exploiting the _ample of
samples_ in the one-bit sampling approach, the one-bit sampled QCS problem
will be formulated as a _linear feasibility_ problem.

###  III-A One-Bit Sampling with Time-Varying Thresholds

Consider a bandlimited signal $y\in L^{2}$, which is to be represented by its
samples via the standard sampling formula[30],

$0<\mathrm{T}\leqslant\frac{\pi}{\Omega},\quad
y(t)=\sum_{k=-\infty}^{k=+\infty}y(k\mathrm{T})\operatorname{sinc}\left(\frac{t}{\mathrm{T}}-k\right),$
(5)

where $1/\mathrm{T}$ is the sampling rate and
$\operatorname{sinc}(t)=\frac{\sin(\pi t)}{(\pi t)}$ is an _ideal_ low-pass
filter. Suppose $y_{k}=y(k\mathrm{T})$ denotes the uniform samples of $y(t)$
with the sampling rate $1/\mathrm{T}$. Let $r_{k}$ denote the quantized
version of $y[k]$ with the formulation $r_{k}=Q(y_{k})$, where $Q$ denotes the
quantization effect. In one-bit quantization, compared to zero or constant
thresholds, time-varying sampling thresholds yield a better reconstruction
performance [31, 32]. These thresholds may be chosen from any distribution. In
this work, to be consistent with state-of-the-art [31, 33, 20], we consider a
Gaussian non-zero time-varying threshold vector
$\bm{\uptau}=\left[\tau_{k}\right]$ that follows the distribution
$\bm{\uptau}\sim\mathcal{N}\left(\mathbf{d}=\mathbf{1}d,\bm{\Sigma}\right)$.
In the case of one-bit quantization with such time-varying sampling
thresholds, the quantizer is simply written as
$r_{k}=\operatorname{sgn}\left(y_{k}-\tau_{k}\right)$. Let
$\mathbf{y}=[y_{k}]$ and $\mathbf{r}=[r_{k}]$. Then, the signal feasibility
based on the one-bit measurements takes the form

$\mathbf{r}\odot\left(\mathbf{y}-\uptau\right)\geq\mathbf{0},$ (6)

or equivalently

$\displaystyle\bm{\Omega}\mathbf{y}$
$\displaystyle\succeq\mathbf{r}\odot\uptau,$ (7)

where $\bm{\Omega}\triangleq\operatorname{diag}\left\\{\mathbf{r}\right\\}$.
Suppose $\mathbf{y},\uptau\in\mathbb{R}^{m}$, and that $\uptau^{(\ell)}$
denotes the time-varying sampling threshold in $\ell$-th experiment where
$\ell\in\mathcal{L}=\\{1,\cdots,m_{1}\\}$. According to (7), for the $\ell$-th
experiment we have

$\displaystyle\bm{\Omega}^{(\ell)}\mathbf{y}$
$\displaystyle\succeq\bm{r}^{(\ell)}\odot\uptau^{(\ell)},\quad\ell\in\mathcal{L},$
(8)

where
$\bm{\Omega}^{(\ell)}=\operatorname{diag}\left\\{\mathbf{r}^{(\ell)}\right\\}$.
In (8), we have $m_{1}$ linear system of inequalities which can be put
together and expressed as

$\tilde{\bm{\Omega}}\mathbf{y}\succeq\operatorname{vec}\left(\mathbf{R}\right)\odot\operatorname{vec}\left(\bm{\Gamma}\right),$
(9)

where $\mathbf{R}$ and $\bm{\Gamma}$ are matrices, with
$\left\\{\mathbf{r}^{(\ell)}\right\\}$ and $\left\\{\uptau^{(\ell)}\right\\}$
representing their columns, respectively, and $\tilde{\bm{\Omega}}$ is given
by

$\tilde{\bm{\Omega}}=\left[\begin{array}[]{c|c|c}\bm{\Omega}^{(1)}&\cdots&\bm{\Omega}^{(m)}\end{array}\right]^{\top},\quad\tilde{\bm{\Omega}}\in\mathbb{R}^{m_{1}m\times
m}.$ (10)

Utilizing the one-bit quantization technique with multiple time-varying
sampling threshold sequences allows for an increase in the number of samples
with little extra cost and serves as a gateway to the realm of few-bit
sampling. This can be especially beneficial in applications where measurement
limitations exist.

###  III-B One-Bit QCS as Linear Feasibility Problem

Hereafter, we will focus on (9) as an overdetermined linear system of
inequalities that is associated with the one-bit sampling scheme. If we apply
one-bit sampling to the QCS (1), referred to as one-bit QCS,

$r^{(\ell)}_{j}=\begin{cases}+1&\operatorname{Tr}\left(\mathbf{A}_{j}\mathbf{X}\right)>\uptau^{(\ell)}_{j},\\\
-1&\operatorname{Tr}\left(\mathbf{A}_{j}\mathbf{X}\right)<\uptau^{(\ell)}_{j}.\end{cases}$
(11)

As a result, by using the linear property of trace function
$\operatorname{Tr}\left(\mathbf{A}_{j}\mathbf{X}\right)=\operatorname{vec}\left(\mathbf{A}^{\top}_{j}\right)^{\top}\operatorname{vec}\left(\mathbf{X}\right)$,
the _one-bit QCS_ polyhedron can be written as

$\displaystyle\mathcal{P}=\left\\{\mathbf{X}\mid
r^{(\ell)}_{j}\operatorname{vec}\left(\mathbf{A}^{\top}_{j}\right)^{\top}\operatorname{vec}\left(\mathbf{X}\right)\geq
r^{(\ell)}_{j}\uptau^{(\ell)}_{j},\leavevmode\nobreak\
\ell\in\mathcal{L},\leavevmode\nobreak\ j\in\mathcal{J}\right\\},$ (12)

which is vectorized based on
$\mathbf{y}=\mathbf{V}\operatorname{vec}\left(\mathbf{X}\right)$, where
$\mathbf{V}$ is a matrix with
$\left\\{\operatorname{vec}\left(\mathbf{A}^{\top}_{j}\right)\right\\}$ as its
rows. The inequality (12) may be recast in the standard polyhedron form as

$\displaystyle\mathcal{P}=\left\\{\mathbf{X}\mid\mathbf{P}\operatorname{vec}\left(\mathbf{X}\right)\succeq\operatorname{vec}\left(\mathbf{R}\right)\odot\operatorname{vec}\left(\bm{\Gamma}\right)\right\\},$
(13)

where $\mathbf{P}=\tilde{\bm{\Omega}}\bm{V}$. By leveraging the sample
abundance in the one-bit sampling, the space constrained by (13), _shrinks_ to
become contained inside the _feasible region_. However, this shrinking space
always contains the globally optimal solution, with a volume that is
decreasing with an increasing number of one-bit samples. We will discuss our
approach to find the desired matrix $\mathbf{X}^{\star}$ below.

##  IV Proposed Algorithm

To recover the desired signal within the one-bit QCS polyhedron, we use an
accelerated variant of randomized Kaczmarz algorithm (RKA). Many variants of
this iterative method and their convergence rates have been proposed and
studied in recent decades for both consistent and inconsistent systems,
including the original randomized Kaczmarz algorithm, the randomized block
Kaczmarz algorithm and most recently, the sampling Kaczmarz-Motzkin (SKM)
method [24, 25, 27]. The block-structured nature of the one-bit QCS matrix has
motivated the development of the SKM method, designed specifically to handle
block-structured linear feasibility problems with efficiency. Further, the
proposed algorithm will be backed by theoretical guarantees.

###  IV-A SKM Method

The SKM is a _subconjugate gradient method_ to solve overdetermined linear
systems, i.e., $\mathbf{B}\mathbf{x}\preceq\mathbf{b}$, where $\mathbf{B}$ is
a $m_{1}m\times n$ matrix. The conjugate-gradient methods immediately turn
such an inequality to an equality of the following form:

$\left(\mathbf{B}\mathbf{x}-\mathbf{b}\right)^{+}=0,$ (14)

and then approach the solution by the same process as used for systems of
equations. Given a sample index set $\mathcal{J}$, without loss of generality,
rewrite (14) as the polyhedron

$\displaystyle\begin{cases}\mathbf{c}_{j}\mathbf{x}\leq
b_{j}&\left(j\in\mathcal{I}_{\leq}\right),\\\
\mathbf{c}_{j}\mathbf{x}=b_{j}&\left(j\in\mathcal{I}_{=}\right),\end{cases}$
(15)

where the disjoint index sets $\mathcal{I}_{\leq}$ and $\mathcal{I}_{=}$
partition $\mathcal{J}$ and $\\{\mathbf{c}_{j}\\}$ are the rows of
$\mathbf{B}$. The projection coefficient $\beta_{i}$ of the SKM at $i$
iteration is [25, 34, 35]

$\beta_{i}=\begin{cases}\left(\mathbf{c}_{j}\mathbf{x}_{i}-b_{j}\right)^{+}&\left(j\in\mathcal{I}_{\leq}\right),\\\
\mathbf{c}_{j}\mathbf{x}_{i}-b_{j}&\left(j\in\mathcal{I}_{=}\right).\end{cases}$
(16)

The central contribution of SKM lies in its innovative way of projection plane
selection. The hyperplane selection is done as follows. At iteration $i$ the
SKM algorithm selects a collection of $\gamma$ (denoted by the set
$\mathcal{T}_{i}$), uniformly at random out of $m_{1}m$ rows of the constraint
matrix $\mathbf{B}$. Then, out of these $\gamma$ rows, the row with maximum
positive residual is selected. Finally, the solution is updated as [27, 36]:
$\mathbf{x}_{i+1}=\mathbf{x}_{i}-\lambda_{i}\frac{\beta_{i}}{\|\mathbf{c}_{{}_{j^{\star}_{i}}}\|^{2}_{2}}\mathbf{c}^{\mathrm{H}}_{{}_{j^{\star}_{i}}}$,
where the index $j^{\star}_{i}$ is chosen as the _Motzkin sampling_ , i.e.,
$j^{\star}_{i}=\operatorname{argmax}\leavevmode\nobreak\
\left\\{\left(\mathbf{c}_{j}\mathbf{x}_{i}-b_{j}\right)^{+}\right\\},\leavevmode\nobreak\
j\in\mathcal{T}_{i}$ at iteration $i$, and $\lambda_{i}$ is a relaxation
parameter which for consistent systems must satisfy,
$0\leq\lim_{i\to\infty}\inf\lambda_{i}\leq\lim_{i\to\infty}\sup\lambda_{i}<2$,
to ensure convergence [24]. The convergence bound for SKM is given by

$\mathbb{E}\left\\{\left\|\mathbf{x}_{i}-\mathbf{x}_{\star}\right\|^{2}_{2}\right\\}\leq\left(1-\frac{2\lambda_{i}-\lambda^{2}_{i}}{\kappa^{2}\left(\mathbf{B}\right)}\right)^{i}\leavevmode\nobreak\
\left\|\mathbf{x}_{0}-\mathbf{x}_{\star}\right\|^{2}_{2},$ (17)

with
$\kappa\left(\mathbf{B}\right)=\|\mathbf{B}\|_{\mathrm{F}}\|\mathbf{B}^{\dagger}\|_{2}$
denoting the scaled condition number, and $\mathbf{x}_{\star}$ is the optimal
solution.

###  IV-B Block SKM Algorithm

The matrix $\mathbf{P}$ in (13) has a block structure with the following
formulation:

$\mathbf{P}=\left[\begin{array}[]{c|c|c}\mathbf{V}^{\top}\bm{\Omega}^{(1)}&\cdots&\mathbf{V}^{\top}\bm{\Omega}^{(m)}\end{array}\right]^{\top},\quad\mathbf{P}\in\mathbb{R}^{m_{1}m\times
n}.$ (18)

Therefore, it is useful to investigate the accelerated block-based RKA methods
to find the desired signal in (13) for further computational efficiency
enhancement. Our proposed algorithm, the _Block SKM_ , is described as
follows. Suppose we have a linear feasibility problem
$\mathbf{B}\mathbf{x}\preceq\mathbf{b}$ where
$\mathbf{B}=\left[\begin{array}[]{c|c|c}\mathbf{B}^{\top}_{1}&\cdots&\mathbf{B}^{\top}_{m_{1}}\end{array}\right]^{\top},$
$\mathbf{B}\in\mathbb{R}^{m_{1}m\times n}$, and
$\mathbf{b}=\left[\begin{array}[]{c|c|c}\mathbf{b}^{\top}_{1}&\cdots&\mathbf{b}^{\top}_{m_{1}}\end{array}\right]^{\top}$.
The proposed algorithm for sparse signal recovery, i.e., the Block SKM, may be
summarized as follows:

  1. 1.

Choose a block $\mathbf{B}_{j}$ uniformly at random with the probability
$\operatorname{Pr}\\{j=k\\}=\frac{\left\|\mathbf{B}_{k}\right\|^{2}_{\mathrm{F}}}{\|\mathbf{B}\|_{\mathrm{F}}^{2}}$.

  2. 2.

Compute $\mathbf{e}=\mathbf{B}_{j}\mathbf{x}-\mathbf{b}_{j}$.

  3. 3.

Let $\mathbf{e}^{\prime}$ denote the sorted version of $\mathbf{e}$ from
$e_{\text{max}}$ (the maximum element of $\mathbf{e}$) to $e_{\text{min}}$
(the minimum element of $\mathbf{e}$). This step is inspired by the idea of
the Motzkin sampling, presented in [27], to have an accelerated convergence.

  4. 4.

Select the first $k^{\prime}<n$ element of $\mathbf{e}^{\prime}$ and construct
the sub-problem
$\mathbf{B}_{j}^{\prime}\mathbf{x}\preceq\mathbf{b}_{j}^{\prime}$, where
$\mathbf{B}_{j}^{\prime}\in\mathbb{R}^{k^{\prime}\times n}$ and
$\mathbf{b}_{j}^{\prime}\in\mathbb{R}^{k^{\prime}\times 1}$. The reason behind
choosing $k^{\prime}<n$ is due to the computation of
$\left(\mathbf{B}_{j}^{\prime}\mathbf{B}_{j}^{\prime\top}\right)^{-1}$ in the
next step (Step $5$). For $k^{\prime}>n$, the matrix
$\mathbf{B}_{j}^{\prime}\mathbf{B}_{j}^{\prime\top}$ is rank-deficient and its
inverse is not available.

  5. 5.

Compute the Moore-Penrose of $\mathbf{B}_{j}^{\prime}$, i.e.,
$\mathbf{B}_{j}^{\prime\dagger}=\mathbf{B}_{j}^{\prime\top}\left(\mathbf{B}_{j}^{\prime}\mathbf{B}_{j}^{\prime\top}\right)^{-1}$.

  6. 6.

Update the solution
$\mathbf{x}_{i+1}=\mathbf{x}_{i}-\lambda_{i}\mathbf{B}_{j}^{\prime\dagger}\left(\mathbf{B}_{j}^{\prime}\mathbf{x}-\mathbf{b}_{j}^{\prime}\right)^{+}$.
This update process is inspired from the randomized block Kaczmarz method [37,
26] which takes advantage of the efficient matrix-vector multiplication, thus
giving the method a significant reduction in computational cost [34].

Particularly, in the case of the one-bit QCS polyhedron,
$\mathbf{B}=-\mathbf{P}$,
$\mathbf{x}=\operatorname{vec}\left(\mathbf{X}\right)$, and
$\mathbf{b}=-\operatorname{vec}\left(\mathbf{R}\right)\odot\operatorname{vec}\left(\bm{\Gamma}\right)$.

###  IV-C Convergence Analysis

It is worth pointing out that the Block SKM algorithm can be considered to be
a special case of the more general _sketch-and-project_ method, defined
as[38]:

$\mathbf{x}_{i+1}=\underset{\mathbf{x}}{\textrm{argmin}}\leavevmode\nobreak\
\left\|\mathbf{x}-\mathbf{x}_{i}\right\|^{2}_{2}\quad\textrm{subject
to}\quad\mathbf{S}^{\top}\mathbf{B}\mathbf{x}\preceq\mathbf{S}^{\top}\mathbf{b},$
(19)

where $\mathbf{S}\in\mathbb{R}^{m_{1}m\times k^{\prime}}$ is the sketch matrix
choosing a block uniformly at random from the main matrix as mentioned in step
$1$. The second step of the proposed algorithm follows the Motzkin sampling
where the index $j^{\star}_{i}$ is chosen in $i$-th iteration as follows:

$j^{\star}_{i}=\underset{j}{\textrm{argmax}}\left\\{\left(\left(\mathbf{S}^{\top}\mathbf{B}\right)_{j}\mathbf{x}_{i}-\left(\mathbf{S}^{\top}\mathbf{b}\right)_{j}\right)^{+}\right\\},$
(20)

with $(\cdot)_{i}$ denoting the $i$th row of the matrix argument.

In the Block SKM algorithm, the sketch matrix is given by

$\mathbf{S}=\left[\begin{array}[]{c|c|c}\mathbf{0}_{k^{\prime}\times
p}&\mathbf{I}_{k^{\prime}}&\mathbf{0}_{k^{\prime}\times(m_{1}m-k^{\prime}-p)}\end{array}\right]^{\top},\leavevmode\nobreak\
\mathbf{S}\in\mathbb{R}^{m_{1}m\times k^{\prime}},$ (21)

where $k^{\prime}$ is the block size and
$p=k^{\prime}\alpha,\leavevmode\nobreak\
\alpha\in\left\\{1,\cdots,\left\lfloor\frac{m_{1}m}{k\prime}\right\rfloor\right\\}$.
Note that the literature does not offer any theoretical guarantees for the
convergence of the Block SKM with the above sketch matrix [39]. To derive our
theoretical guarantees for the algorithm used to solve the one-bit QCS, we
change the sketch matrix to the _Gaussian_ sketch matrix as follows:

$\mathbf{S}=\left[\begin{array}[]{c|c|c}\mathbf{0}_{k^{\prime}\times
p}&\mathbf{G}&\mathbf{0}_{k^{\prime}\times(m_{1}m-k^{\prime}-p)}\end{array}\right]^{\top},\leavevmode\nobreak\
\mathbf{S}\in\mathbb{R}^{m_{1}m\times k^{\prime}},$ (22)

where $\mathbf{G}$ is a $k^{\prime}\times k^{\prime}$ Gaussian matrix, whose
entries are i.i.d. following the distribution $\mathcal{N}\left(0,1\right)$.
In this framework, we are able to provide some theoretical guarantees by
taking advantage of the favorable properties of Gaussian random variables.

Assume that $\mathcal{S}$ denotes a non-empty solution set of the polyhedron
(13). Owing to the fact that
$\mathbb{E}\left\\{\left\|\mathbf{x}_{i+1}-\mathcal{S}\right\|_{2}^{2}\right\\}\leq\mathbb{E}\left\\{\left\|\mathbf{x}_{i+1}-\mathbf{x}_{\star}\right\|_{2}^{2}\right\\}$
[25], then we proceed to prove the convergence rate by employing
$\mathbb{E}\left\\{\left\|\mathbf{x}_{i+1}-\mathbf{x}_{\star}\right\|_{2}^{2}\right\\}$.
Using the fact that $\mathbf{x}_{i+1}-\mathbf{x}_{\star}$ is orthogonal to
$(\mathbf{S}^{T}\mathbf{B})_{j^{\star}_{i}}$ [39], where $j^{\star}_{i}$ is
the index chosen based on the Motzkin sampling for the $i$-th iteration, we
have the following Pythagorean relation [39, 38]:

$\displaystyle\left\|\mathbf{x}_{i+1}-\mathbf{x}_{\star}\right\|_{2}^{2}=\left\|\mathbf{x}_{i}-\mathbf{x}_{\star}\right\|_{2}^{2}-\frac{\left\|\left(\left(\mathbf{S}^{\top}\mathbf{B}\right)_{j^{\star}_{i}}\mathbf{x}_{i}-\left(\mathbf{S}^{\top}\mathbf{b}\right)_{j^{\star}_{i}}\right)^{+}\right\|_{2}^{2}}{\left\|\left(\mathbf{S}^{\top}\mathbf{B}\right)_{i}\right\|_{2}^{2}}.$
(23)

In the linear inequality system, the Kaczmarz algorithms only updates the
solution when
$\mathbf{S}^{\top}\mathbf{B}\mathbf{x}_{i}\succeq\mathbf{S}^{\top}\mathbf{b}$
at $i$-th iteration. Therefore, one can readily rewrite (23) at iteration $i$
where the condition
$\mathbf{S}^{\top}\mathbf{B}\mathbf{x}_{i}\succeq\mathbf{S}^{\top}\mathbf{b}$
is met:

$\displaystyle\left\|\mathbf{x}_{i+1}-\mathbf{x}_{\star}\right\|_{2}^{2}=\left\|\mathbf{x}_{i}-\mathbf{x}_{\star}\right\|_{2}^{2}-\frac{\left\|\mathbf{S}^{\top}\mathbf{B}\mathbf{x}_{i}-\mathbf{S}^{\top}\mathbf{b}\right\|_{\infty}^{2}}{\left\|\left(\mathbf{S}^{\top}\mathbf{B}\right)_{j^{\star}_{i}}\right\|_{2}^{2}}.$
(24)

By taking the expectation over the error, we have

$\mathbb{E}_{\mathbf{S}}\left\\{\left\|\mathbf{x}_{i+1}-\mathbf{x}_{\star}\right\|_{2}^{2}\right\\}=\left\|\mathbf{x}_{i}-\mathbf{x}_{\star}\right\|_{2}^{2}-\mathbb{E}_{\mathbf{S}}\left\\{\frac{\left\|\mathbf{S}^{\top}\mathbf{B}\mathbf{x}_{i}-\mathbf{S}^{\top}\mathbf{b}\right\|_{\infty}^{2}}{\left\|\left(\mathbf{S}^{\top}\mathbf{B}\right)_{j^{\star}_{i}}\right\|_{2}^{2}}\right\\}.$
(25)

In addition, we have that

$\displaystyle\mathbb{E}_{\mathbf{S}}\left\\{\left\|\left(\mathbf{S}^{\top}\mathbf{B}\right)_{j^{\star}_{i}}\right\|_{2}^{2}\right\\}=\sum_{k=1}^{n}\mathbb{E}_{\mathbf{S}}\left\\{\left(\sum_{i_{1}=1}^{m_{1}m}\mathbf{S}_{ji_{1}}\mathbf{B}_{i_{1}i_{2}}\right)^{2}\right\\},$
(26)

or equivalently, in terms of $\mathbf{G}$ in (22),

$\displaystyle\sum_{i_{2}=1}^{n}\mathbb{E}_{\mathbf{G}}\left\\{\left(\sum_{i_{1}=1}^{k^{\prime}}\mathbf{G}^{\top}_{ji_{1}}\mathbf{B}_{i_{1}i_{2}}\right)^{2}\right\\}$
$\displaystyle=$ (27) $\displaystyle\sum_{i_{2}=1}^{n}$
$\displaystyle\sum_{i_{1}=1}^{k^{\prime}}\mathbb{E}_{\mathbf{G}}\left\\{\left(\mathbf{G}^{\top}_{ji_{1}}\right)^{2}\right\\}\mathbf{B}^{2}_{i_{1}i_{2}},$

with
$\mathbb{E}_{\mathbf{G}}\left\\{\left(\mathbf{G}^{\top}_{ji_{1}}\right)^{2}\right\\}=1$,
which helps to simplify (27) as

$\displaystyle\sum_{i_{2}=1}^{n}\sum_{i_{1}=1}^{k^{\prime}}\mathbf{B}^{2}_{i_{1}i_{2}}=\|\hat{\mathbf{B}}\|_{\mathrm{F}}^{2},$
(28)

where $\hat{\mathbf{B}}$ is the $k^{\prime}\times n$ submatrix of
$\mathbf{B}$. Due to the fact that the second term in the right-hand side of
(25) is an expectation over the convex function $f(x,y)=x^{2}/y$, we can apply
Jensen’s inequality as follows:

$\displaystyle\mathbb{E}_{\mathbf{S}}\left\\{\frac{\left\|\mathbf{S}^{\top}\mathbf{B}\mathbf{x}_{i}-\mathbf{S}^{\top}\mathbf{b}\right\|_{\infty}^{2}}{\left\|\left(\mathbf{S}^{\top}\mathbf{B}\right)_{j^{\star}_{i}}\right\|_{2}^{2}}\right\\}\geq\frac{\left(\mathbb{E}_{\mathbf{S}}\left\\{\left\|\mathbf{S}^{\top}\mathbf{B}\mathbf{x}_{i}-\mathbf{S}^{\top}\mathbf{b}\right\|_{\infty}\right\\}\right)^{2}}{\mathbb{E}_{\mathbf{S}}\left\\{\left\|\left(\mathbf{S}^{\top}\mathbf{B}\right)_{j^{\star}_{i}}\right\|_{2}^{2}\right\\}}.$
(29)

Since
$\mathbf{S}^{\top}\mathbf{B}\mathbf{x}_{\star}\preceq\mathbf{S}^{\top}\mathbf{b}$
and
$\mathbf{S}^{\top}\mathbf{B}\mathbf{x}_{i}\succeq\mathbf{S}^{\top}\mathbf{b}$,
one can conclude

$\left\|\mathbf{S}^{\top}\mathbf{B}\mathbf{x}_{i}-\mathbf{S}^{\top}\mathbf{b}\right\|_{\infty}\geq\left\|\mathbf{S}^{\top}\mathbf{B}\mathbf{x}_{i}-\mathbf{S}^{\top}\mathbf{B}\mathbf{x}_{\star}\right\|_{\infty}.$
(30)

It follows from the above that

$\displaystyle\mathbb{E}_{\mathbf{S}}\left\\{\frac{\left\|\mathbf{S}^{\top}\mathbf{B}\mathbf{x}_{i}-\mathbf{S}^{\top}\mathbf{b}\right\|_{\infty}^{2}}{\left\|\left(\mathbf{S}^{\top}\mathbf{B}\right)_{j^{\star}_{i}}\right\|_{2}^{2}}\right\\}\geq\frac{\left(\mathbb{E}_{\mathbf{S}}\left\\{\left\|\mathbf{S}^{\top}\mathbf{B}\left(\mathbf{x}_{i}-\mathbf{x}_{\star}\right)\right\|_{\infty}\right\\}\right)^{2}}{\|\hat{\mathbf{B}}\|_{\mathrm{F}}^{2}}.$
(31)

We can additionally take advantage of the estimate for the maximum of
independent normal random variables[39],

$\displaystyle\mathbb{E}_{\mathbf{S}}\left\\{\left\|\mathbf{S}^{\top}\mathbf{B}\left(\mathbf{x}_{i}-\mathbf{x}_{\star}\right)\right\|_{\infty}\right\\}$
$\displaystyle=\mathbb{E}_{\mathbf{S}}\left\\{\max_{t\in[k^{\prime}]}\left\langle\mathbf{s}_{t},\mathbf{B}\left(\mathbf{x}_{i}-\mathbf{x}_{\star}\right)\right\rangle\right\\}$
(32)
$\displaystyle=\mathbb{E}_{\mathbf{G}}\left\\{\max_{t\in[k^{\prime}]}\left\langle\mathbf{s}_{t},\hat{\mathbf{B}}\left(\mathbf{x}_{i}-\mathbf{x}_{\star}\right)\right\rangle\right\\}$
$\displaystyle\geq
c\|\hat{\mathbf{B}}\left(\mathbf{x}_{i}-\mathbf{x}_{\star}\right)\|_{2}\sqrt{\log
k^{\prime}},$

where $\mathbf{s}_{t}$ is the $t$-th column of $\mathbf{S}$,
$[k^{\prime}]=\left\\{1,2,\cdots,k^{\prime}\right\\}$, and $c$ is a positive
value. By plugging the inequality (32) into (25), and using the inequality,

$\left\|\hat{\mathbf{B}}\left(\mathbf{x}_{i}-\mathbf{x}_{\star}\right)\right\|_{2}^{2}\geq\sigma^{2}_{\textrm{min}}\left(\hat{\mathbf{B}}\right)\left\|\mathbf{x}_{i}-\mathbf{x}_{\star}\right\|^{2}_{2},$
(33)

where $\sigma^{2}_{\textrm{min}}$ is the minimum singular value. Thus, we
obtain

$\displaystyle\mathbb{E}\left\\{\left\|\mathbf{x}_{i+1}-\mathbf{x}_{\star}\right\|_{2}^{2}\right\\}$
$\displaystyle\leq\left\|\mathbf{x}_{i}-\mathbf{x}_{\star}\right\|_{2}^{2}-\frac{c\|\mathbf{B}\left(\mathbf{x}_{\star}-\mathbf{x}_{i}\right)\|_{2}^{2}\log
k^{\prime}}{\|\mathbf{B}\|_{\mathrm{F}}^{2}}$ (34)
$\displaystyle\leq\left\|\mathbf{x}_{i}-\mathbf{x}_{\star}\right\|_{2}^{2}-\frac{c\sigma_{\min}^{2}(\hat{\mathbf{B}})\log
k^{\prime}}{\|\hat{\mathbf{B}}\|_{\mathrm{F}}^{2}}\left\|\mathbf{x}_{i}-\mathbf{x}_{\star}\right\|_{2}^{2}$
$\displaystyle\leq\left(1-\frac{c\sigma_{\min}^{2}(\hat{\mathbf{B}})\log
k^{\prime}}{\|\hat{\mathbf{B}}\|_{\mathrm{F}}^{2}}\right)\left\|\mathbf{x}_{i}-\mathbf{x}_{\star}\right\|_{2}^{2},$

which can be recast as the following _convergence rate_ , after $K$ updates:

$\displaystyle\mathbb{E}\left\\{\left\|\mathbf{x}_{i+1}-\mathbf{x}_{\star}\right\|_{2}^{2}\right\\}\leq\left(1-\frac{c\sigma_{\min}^{2}(\hat{\mathbf{B}})\log
k^{\prime}}{\|\hat{\mathbf{B}}\|_{\mathrm{F}}^{2}}\right)^{K}\left\|\mathbf{x}_{0}-\mathbf{x}_{\star}\right\|_{2}^{2}.$
(35)

(a)  (b)

Figure 1: Average NMSE for the error between the desired $\mathbf{X}^{\star}$
and its recovered version $\bar{\mathbf{X}}$ for various numbers of time-
varying sampling threshold sequences $m_{1}\in\left\\{10,50,100,150\right\\}$
when the block SKM is utilized in with (a) $\|\mathbf{x}\|_{0}=5$, (b)
$\|\mathbf{x}\|_{0}=10$.

(a)  (b)

Figure 2: Average NMSE for the error between the desired $\mathbf{X}^{\star}$
and its recovered version $\bar{\mathbf{X}}$ for various numbers of time-
varying sampling threshold sequences $m_{1}\in\left\\{10,50,100,150\right\\}$
when the block SKM is utilized in with (a) $\|\mathbf{x}\|_{0}=5$, (b)
$\|\mathbf{x}\|_{0}=10$.

##  V Numerical Results

In this section, at first, we numerically scrutinize the capability of the
block SKM in the one-bit QCS problem by evaluating the squared Frobenius norm
of the error between the desired matrix $\mathbf{X}^{\star}$ and its estimate
$\bar{\mathbf{X}}$, normalized by the squared Frobenius norm of the desired
matrix:

$\mathrm{NMSE}\triangleq\frac{\left\|\mathbf{X}^{\star}-\bar{\mathbf{X}}\right\|^{2}_{\mathrm{F}}}{\left\|\mathbf{X}^{\star}\right\|^{2}_{\mathrm{F}}}.$
(36)

The input signal $\mathbf{x}\in\mathbb{R}^{64}$, is considered to be a sparse
signal with (i) $\|\mathbf{x}\|_{0}=5$, and (ii) $\|\mathbf{x}\|_{0}=10$. To
choose the time-varying sampling thresholds, we consider the framework
presented in [22], which relies on knowledge of the dynamic range of the
measurements $\mathbf{y}$. Assume
$\beta_{\mathbf{y}}=\left\|\mathbf{y}\right\|_{\infty}$ denotes the dynamic
range of the measurements. Then, herein we generate the time-varying sampling
thresholds as
$\left\\{\mathbf{\uptau}^{(\ell)}\sim\mathcal{N}\left(\mathbf{0},\frac{\beta_{\mathbf{y}}^{2}}{9}\mathbf{I}_{5000}\right)\right\\}_{\ell=1}^{m_{1}}$.
Each sensing matrix is generated based on
$\mathbf{A}_{j}=\mathbf{a}_{j}\mathbf{a}^{\mathrm{H}}_{j}$, where
$\mathbf{a}_{j}\sim\mathcal{N}\left(\mathbf{0},\mathbf{I}_{64}\right)$. We
solve the overdetermined one-bit QCS polyhedron in (13) via the Block SKM for
the number of time-varying sampling threshold sequences
$m_{1}\in\left\\{10,50,100,150\right\\}$. Fig. 1 appears to confirm the
possibility of recovering the desired matrix $\mathbf{X}^{\star}$ in the one-
bit QCS polyhedron (13) by applying Block SKM. As expected, the performance of
the recovery will be significantly enhanced as the number of time-varying
sampling threshold sequences grows large. The reason behind this observation
is the sample abundance condition which has been initially analyzed and proved
in [15, Theorem 1] and extended to another sampling scheme in [20]. Note that
the results in Fig. 1 are averaged over $15$ experiments. To examine the
performance of the proposed algorithm for the full rank $\mathbf{A}_{j}$
scenario, we generate a full rank $\mathbf{A}_{j}\in\mathbb{R}^{64\times 64}$
where its entries are i.i.d normal random variables. Similarly, we generate
time-varying sampling thresholds as
$\left\\{\mathbf{\uptau}^{(\ell)}\sim\mathcal{N}\left(\mathbf{0},\frac{\beta_{\mathbf{y}}^{2}}{9}\mathbf{I}_{5000}\right)\right\\}_{\ell=1}^{m_{1}}$.
Fig. 2 illustrates the recovery performance of the Block SKM in this case
while preserving the property of boosting the recovery error as the number of
time-varying sampling thresholds grows large. Each data point in Fig. 2 is
averaged over $15$ experiments. Note that the algorithm proposed employs only
low-resolution (one-bit) samples, but capitalizes on their abundance to
converge to the global solution with heightened precision as the quantity of
one-bit samples increases.

Figure 3: Comparing the recovery performance of the proposed Kaczmarz-based algorithm, namely the Block SKM, with that of SKM and RKA in terms of NMSE for a linear system of inequalities.  Table I: Comparing CPU times and $\operatorname{NMSE}$ of Block SKM and GESPAR. Algorithm | $m^{\star}$ | CPU time (s) | $\operatorname{NMSE}$
---|---|---|---
Block SKM (one-bit) | $5000$ | $0.0026$ | $3.1072e-7$
GESPAR | $128$ | $0.0041$ | $2.4382e-5$

Moreover, we numerically compare the RKA[25], SKM[27], and our proposed Block
SKM in linear systems of inequalities. We apply one-bit sampling to a system
of linear equalities $\mathbf{B}\mathbf{x}=\mathbf{y}$, resulting in the
creation of its corresponding system of linear inequalities as described in
(9). Herein, we consider $\mathbf{B}\in\mathbb{R}^{100\times 10}$,
$\mathbf{x}\in\mathbb{R}^{10}$, and $\mathbf{y}\in\mathbb{R}^{100}$. Each row
of $\mathbf{B}$ is generated as
$\mathbf{b}_{j}\sim\mathcal{N}\left(\mathbf{0},\mathbf{I}_{10}\right)$. Also,
the desired signal $\mathbf{x}$ is generated as
$\mathbf{x}\sim\mathcal{N}\left(\mathbf{0},\mathbf{I}_{10}\right)$.
Accordingly, we generate time-varying sampling thresholds as
$\left\\{\mathbf{\uptau}^{(\ell)}\sim\mathcal{N}\left(\mathbf{0},\frac{\beta_{\mathbf{y}}^{2}}{9}\mathbf{I}_{100}\right)\right\\}_{\ell=1}^{m_{1}}$
for $m_{1}=40$. The performance of the RKA, SKM, and Block SKM is illustrated
in Fig. 3. The results show that the Block SKM outperforms the other two
approaches, delivering a faster recovery and higher accuracy in the recovery
of the desired signal $\mathbf{x}$. The normalized mean square error for the
signal is defined as
$\operatorname{NMSE}\triangleq\frac{\left\|\mathbf{x}_{\star}-\bar{\mathbf{x}}\right\|_{2}^{2}}{\left\|\mathbf{x}_{\star}\right\|_{2}^{2}}$,
where $\mathbf{x}_{\star}$ and $\bar{\mathbf{x}}$ denote the true discretized
signal and its recovered version, respectively. The NMSE results in Fig. 3 are
averaged over $15$ experiments.

To further investigate the efficacy of the proposed algorithm in QCS, we
compare our proposed approach with the well-known GESPAR approach with the
initialization algorithm proposed in [3] in terms of NMSE and CPU time. As
presented in Table I, Block SKM outperforms GESPAR in terms of both NMSE and
CPU time. The results are obtained for $\mathbf{x}\in\mathbb{R}^{64}$ when the
optimal number of samples are utilized, and where $m^{\star}=2n=128$ (high-
resolution samples) and $m^{\star}=5000$ (one-bit samples) are considered for
the high-resolution method and one-bit QCS, respectively. Herein, the
optimality of sample sizes means that the number of samples utilized by
algorithms leads to their best performance (up to global phase), i.e.
satisfying the criterion
$\left\|\mathbf{x}_{\star}-\bar{\mathbf{x}}\right\|_{2}^{2}\leq 5\times
10^{-5}\left\|\mathbf{x}_{\star}\right\|_{2}^{2}$. By this comparison, we
remove the burden of the large number of samples from the GESPAR to fairly
compare their optimal shape deploying incomplete measurements with that of the
Block SKM. Note that the signal of interest is obtained from
$\bar{\mathbf{X}}=\bar{\mathbf{x}}\bar{\mathbf{x}}^{\mathrm{H}}$, where the
signal is the largest eigenvector of the recovered matrix.

##  VI Conclusion

We propose taking advantage of the abundant number of samples available in
one-bit sampling with time-varying thresholds to efficiently and globally
solve the quadratic compressed sensing problem. In particular, a state-of-the-
art randomized Kaczmarz algorithms is proposed to find the desired signal
inside the emerging confined feasible regions, named the one-bit polyhedron,
with an enhanced convergence rate. The numerical results showcased the
effectiveness of the proposed approaches for the quadratic compressed sensing
problem.

## References

  * [1] A. Beck and Y. Eldar, “Sparsity constrained nonlinear optimization: Optimality conditions and algorithms,” _SIAM Journal on Optimization_ , vol. 23, no. 3, pp. 1480–1509, 2013.
  * [2] Y. Shechtman, Y. Eldar, A. Szameit, and M. Segev, “Sparsity based sub-wavelength imaging with partially incoherent light via quadratic compressed sensing,” _Optics express_ , vol. 19, no. 16, pp. 14 807–14 822, 2011.
  * [3] Y. Shechtman, A. Beck, and Y. Eldar, “GESPAR: Efficient phase retrieval of sparse signals,” _IEEE Transactions on Signal Processing_ , vol. 62, no. 4, pp. 928–938, 2014.
  * [4] K. Jaganathan, Y. Eldar, and B. Hassibi, “Phase retrieval: An overview of recent developments,” _Optical Compressive Imaging_ , pp. 279–312, 2016\.
  * [5] K. Jaganathan, S. Oymak, and B. Hassibi, “Sparse phase retrieval: Convex algorithms and limitations,” in _2013 IEEE International Symposium on Information Theory_. IEEE, 2013, pp. 1022–1026.
  * [6] E. J. Candes, T. Strohmer, and V. Voroninski, “PhaseLift: Exact and stable signal recovery from magnitude measurements via convex programming,” _Communications on Pure and Applied Mathematics_ , vol. 66, no. 8, pp. 1241–1274, 2013.
  * [7] E. Candès and X. Li, “Solving quadratic equations via PhaseLift when there are about as many equations as unknowns,” _Foundations of Computational Mathematics_ , vol. 14, no. 5, pp. 1017–1026, 2014.
  * [8] G. Wang, G. Giannakis, and Y. Eldar, “Solving systems of random quadratic equations via truncated amplitude flow,” _IEEE Transactions on Information Theory_ , vol. 64, no. 2, pp. 773–794, 2017.
  * [9] T. Bendory, Y. Eldar, and N. Boumal, “Non-convex phase retrieval from STFT measurements,” _IEEE Transactions on Information Theory_ , vol. 64, no. 1, pp. 467–484, 2017.
  * [10] E. J. Candès, X. Li, and M. Soltanolkotabi, “Phase retrieval via Wirtinger flow: Theory and algorithms,” _IEEE Transactions on Information Theory_ , vol. 61, no. 4, pp. 1985–2007, Apr 2015.
  * [11] I. Waldspurger, “Phase retrieval for wavelet transforms,” _IEEE Transactions on Information Theory_ , vol. 63, no. 5, pp. 2993–3009, 2017.
  * [12] A. Eamaz, F. Yeganegi, and M. Soltanalian, “Modified arcsine law for one-bit sampled stationary signals with time-varying thresholds,” in _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021, pp. 5459–5463.
  * [13] P. Boufounos and R. Baraniuk, “1-bit compressive sensing,” in _2008 42nd Annual Conference on Information Sciences and Systems_. IEEE, 2008, pp. 16–21.
  * [14] E. J. Candes, Y. C. Eldar, T. Strohmer, and V. Voroninski, “Phase retrieval via matrix completion,” _SIAM review_ , vol. 57, no. 2, pp. 225–251, 2015\.
  * [15] A. Eamaz, F. Yeganegi, and M. Soltanalian, “One-bit phase retrieval: More samples means less complexity?” _IEEE Transactions on Signal Processing_ , vol. 70, pp. 4618–4632, 2022.
  * [16] A. Mezghani and A. L. Swindlehurst, “Blind estimation of sparse broadband massive MIMO channels with ideal and one-bit ADCs,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 11, pp. 2972–2983, 2018.
  * [17] L. Jacques, J. Laska, P. Boufounos, and R. Baraniuk, “Robust 1-bit compressive sensing via binary stable embeddings of sparse vectors,” _IEEE Transactions on Information Theory_ , vol. 59, no. 4, pp. 2082–2102, 2013.
  * [18] P. Boufounos, L. Jacques, F. Krahmer, and R. Saab, “Quantization and compressive sensing,” in _Compressed Sensing and its Applications: MATHEON Workshop 2013_. Springer, 2015, pp. 193–237.
  * [19] A. Eamaz, F. Yeganegi, and M. Soltanalian, “Covariance recovery for one-bit sampled stationary signals with time-varying sampling thresholds,” _Signal Processing_ , vol. 206, p. 108899, 2023.
  * [20] A. Eamaz, K. V. Mishra, F. Yeganegi, and M. Soltanalian, “UNO: Unlimited sampling meets one-bit quantization,” _arXiv preprint arXiv:2301.10155_ , 2022.
  * [21] C. Xu and L. Jacques, “Quantized compressive sensing with RIP matrices: The benefit of dithering,” _Information and Inference: A Journal of the IMA_ , vol. 9, no. 3, pp. 543–586, 2020.
  * [22] J. N. Laska, Z. Wen, W. Yin, and R. G. Baraniuk, “Trust, but verify: Fast and accurate signal recovery from 1-bit compressive measurements,” _IEEE Transactions on Signal Processing_ , vol. 59, no. 11, pp. 5289–5301, 2011.
  * [23] S. Kaczmarz, “Angenäherte auflösung von systemen linearer gleichungen (english translation by Jason Stockmann): Bulletin international de l’académie polonaise des sciences et des lettres,” 1937.
  * [24] T. Strohmer and R. Vershynin, “A randomized Kaczmarz algorithm with exponential convergence,” _Journal of Fourier Analysis and Applications_ , vol. 15, no. 2, pp. 262–278, 2009.
  * [25] D. Leventhal and A. S. Lewis, “Randomized methods for linear constraints: convergence rates and conditioning,” _Mathematics of Operations Research_ , vol. 35, no. 3, pp. 641–654, 2010.
  * [26] D. Needell and J. A. Tropp, “Paved with good intentions: Analysis of a randomized block Kaczmarz method,” _Linear Algebra and its Applications_ , vol. 441, pp. 199–221, 2014.
  * [27] J. De Loera, J. Haddock, and D. Needell, “A sampling Kaczmarz–Motzkin algorithm for linear feasibility,” _SIAM Journal on Scientific Computing_ , vol. 39, no. 5, pp. S66–S87, 2017.
  * [28] Z. Q. Luo, W. K. Ma, A. M. C. So, Y. Ye, and S. Zhang, “Semidefinite relaxation of quadratic optimization problems,” _IEEE Signal Processing Magazine_ , vol. 27, no. 3, pp. 20–34, 2010.
  * [29] M. M. Naghsh, M. Soltanalian, P. Stoica, M. Modarres-Hashemi, A. De Maio, and A. Aubry, “A doppler robust design of transmit sequence and receive filter in the presence of signal-dependent interference,” _IEEE Transactions on Signal Processing_ , vol. 62, no. 4, pp. 772–785, 2013.
  * [30] A. Bhandari, F. Krahmer, and R. Raskar, “On unlimited sampling and reconstruction,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 3827–3839, 2020.
  * [31] A. Ameri, J. Li, and M. Soltanalian, “One-bit radar processing and estimation with time-varying sampling thresholds,” in _2018 IEEE 10th Sensor Array and Multichannel Signal Processing Workshop (SAM)_. IEEE, 2018, pp. 208–212.
  * [32] A. Eamaz, F. Yeganegi, and M. Soltanalian, “Covariance recovery for one-bit sampled non-stationary signals with time-varying sampling thresholds,” _IEEE Transactions on Signal Processing_ , vol. 70, pp. 5222–5236, 2022.
  * [33] S. Khobahi and M. Soltanalian, “Signal recovery from 1-bit quantized noisy samples via adaptive thresholding,” in _2018 52nd Asilomar Conference on Signals, Systems, and Computers_. IEEE, 2018, pp. 1757–1761.
  * [34] J. Briskman and D. Needell, “Block Kaczmarz method with inequalities,” _Journal of Mathematical Imaging and Vision_ , vol. 52, no. 3, pp. 385–396, 2015.
  * [35] L. Dai, M. Soltanalian, and K. Pelckmans, “On the randomized Kaczmarz algorithm,” _IEEE Signal Processing Letters_ , vol. 21, no. 3, pp. 330–333, 2013.
  * [36] M. Sarowar Morshed and M. Saiful Islam, “Sampling Kaczmarz Motzkin method for linear feasibility problems: Generalization & acceleration,” _arXiv e-prints_ , pp. arXiv–2002, 2020.
  * [37] T. Elfving, “Block-iterative methods for consistent and inconsistent linear equations,” _Numerische Mathematik_ , vol. 35, no. 1, pp. 1–12, 1980.
  * [38] M. Dereziński and E. Rebrova, “Sharp analysis of sketch-and-project methods via a connection to randomized singular value decomposition,” _arXiv preprint arXiv:2208.09585_ , 2022.
  * [39] E. Rebrova and D. Needell, “On block gaussian sketching for the kaczmarz method,” _Numerical Algorithms_ , vol. 86, pp. 443–473, 2021.