# ORKA: Accelerated Kaczmarz Algorithms for Signal Recovery from One-Bit
Samples

Arian Eamaz, , Farhang Yeganegi, Deanna Needell, , and
Mojtaba Soltanalian This work was supported in part by National Science
Foundation Grant CCF-1704401. The first two authors contributed equally to
this work.A. Eamaz, F. Yeganegi and M. Soltanalian are with the Department of
Electrical and Computer Engineering, University of Illinois Chicago, Chicago,
IL 60607, USA (_Corresponding author: Arian Eamaz_). D. Needell is with the
Department of Mathematics, University of California Los Angeles, Los Angeles,
CA 90095 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 the nuclear norm minimization, 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 make this achievable, we propose enhanced randomized Kaczmarz
algorithms to tackle these highly overdetermined feasibility problems. Several
numerical results are presented to illustrate the effectiveness of the
proposed methodologies.

###### Index Terms:

Convex-relaxed problems, nuclear norm minimization, one-bit quantization, one-
bit ADCs, randomized Kaczmarz algorithm, statistical signal processing, time-
varying sampling thresholds.

##  I Introduction

We consider an optimization problem of the form

$\displaystyle\min_{\mathbf{X}}$ $\displaystyle f(\mathbf{X})$ (1) s.t.
$\displaystyle\mathcal{A}\left(\mathbf{X}\right)=\mathbf{y},$
$\displaystyle\mathbf{X}\in\Omega_{c},$

where $f(.)$ is a cost function, $\mathbf{X}\in\mathbb{C}^{n_{1}\times n_{2}}$
is the matrix of unknowns, $\mathbf{y}\in\mathbb{R}^{n}$ is the measurement
vector, and $\mathcal{A}$ is a linear transformation mapping
$\mathbb{C}^{n_{1}\times n_{2}}$ into $\mathbb{R}^{n}$.

This problem has been used as a relaxed version of some well-known NP-hard
problems, and emerging in wide variety of statistical signal processing
applications. Although many problems can be expressed in the form in (1), the
applications we will focus on in this paper include some specific problems of
interest in statistical signal processing, which can take advantage of low-
resolution (and particularly one bit) sampling and processing:

  * •

Low-rank matrix recovery: The task of recovering a low-rank matrix from its
linear measurements plays a central role in computational science. The problem
occurs in many areas of applied mathematics such as signal processing [1, 2,
3, 4, 5, 6, 7], machine learning [8, 9, 10, 11, 12, 13], and computer
vision[14]. In this scenario, the cost function of (1), $f(.)$, is typically
to be the _nuclear norm_ or the Frobenius norm, and the constraint set
$\Omega_{c}$ would be a amplitude restriction on the elements of matrix
$\mathbf{X}$; see [15, 2].

  * •

Phase retrieval: Phase retrieval has received a great deal of interest as it
aims to recover an unknown signal solely from phaseless measurements that
depend on the signal through a linear observation, commanding numerous
applications in applied physics and statistical signal processing communities
over the past decades [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28]. To
have a convex formulation, the phase retrieval has been relaxed into semi-
definite programs where the problem boils down to a trace minimization while
considering the positive semi-definite constraint [16, 1].

  * •

Compressed sensing: Compressed sensing (CS) offers a framework for
simultaneous sensing and compression of finite dimensional vectors, that
relies on linear dimensionality reduction. Through a CS formulation, sparse
signals may be recovered from highly incomplete measurements[29]. The problem
(1) can be adopted in the CS content when
$f\left(\mathbf{X}\right)=\left\|\operatorname{vec}\left(\mathbf{X}\right)\right\|_{1}$.

  * •

Magnetic resonance imaging: Reconstructing magnetic resonance images commonly
involves collecting a series of frames of data in which a radio frequency
excitation produces new transverse magnetization, which is then sampled along
a particular trajectory in $k$-sparse representation. Due to meet various
physical and physiological constraints, most MRI methods utilize a sequence of
acquisitions, each of which partially samples the representation. Let the
acquired sequence of measurements be represented by $y_{i}$, where $i$ is the
sequence index, and $\left\\{\mathcal{A}_{i}\left(\mathbf{X}\right)\right\\}$
denote a linear transformation, chosen in a manner that promotes sparsity in
the range space. In this example, the cost function can be considered to be
the $\ell_{1}$-norm, and the sequence of acquisitions are used as linear
constraints in (1)[30, 31].

Sampling the signals of interest at high data rates with high-resolution ADCs
would dramatically increase the overall manufacturing cost and power
consumption of such ADCs. In multi-bit sampling scenarios, a very large number
of quantization levels is necessary in order to represent the original
continuous signal in with high accuracy, which in turn leads to a considerable
reduction in sampling rate [32]. This attribute of multi-bit sampling is the
key reason for the general emergence of underdetermined systems
$n_{1}n_{2}\geq n$ in (1) [27, 1, 16]. 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 ADCs, producing sign data ($\pm 1$). This enables signal processing
equipments to sample at a very high rate, with a considerably lower cost and
energy consumption, compared to their counterparts which employ multi-bit ADCs
[33, 34, 32, 35].

In traditional one-bit sampling schemes, the signal recovery is accomplished
by comparing the signal with a fixed threshold, usually zero. This creates
some difficulties in estimating signal parameters. In contrast, recent works
have employed time-varying sampling thresholds, which exhibit enhanced
recovery performance for the signal parameters [32, 36, 37, 16, 38, 39, 40].

In this paper, we consider the deployment of one-bit sampling with time-
varying thresholds,leading to an increased sample size and a _highly
overdetermined system_ as a result. The proposed _O_ ne-bit aided _R_
andomized _K_ aczmarz _Algorithm_ , which we refer to as ORKA, can find the
desired signal $\mathbf{X}^{\star}$ in (1) 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 (1), and (ii) solving this obtained
linear feasibility problem by leveraging one of the efficient solver families
of overdetermined systems, _Kaczmarz algorithms_. The Kaczmarz method [41] 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. Each iteration projects onto the solution space
corresponding to one row in the linear system, in a sequential regimen. The
method has been applied to various applications in image reconstruction,
digital signal processing, and computer tomography [42, 43, 16]. 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 method [44, 45, 46,
47, 48].

###  I-A Contributions of the Paper

In [16], we showed that the sheer number of measurements acquired in one bit
sampling facilitates recovering the signal of interest in a less costly manner
by making costly constraints such as semidefiniteness and rank redundant.
Then, a simple randomized Kaczmarz algorithm (RKA) was utilized to solve the
obtained linear feasibility problem. This idea is generalized in this paper to
(1) where we generate the abundant samples and eventually introduce a one-bit
linear feasibility region named the _one-bit polyhedron_. In other words, by
using this technique, we make (1) a large-scale overdetermined system which is
the desired application setting for Kaczmarz algorithms.

To solve our highly overdetermined system, we propose two novel variants of
RKA which will be compared with the existing RKA variants. Furthermore, an
algorithm is proposed based on our model to adaptively evaluate the time-
varying sampling thresholds. The convergence rate of the proposed algorithm is
investigated based on the moments generating function of recovery errors and
the scaled condition number of the constraint matrix. Finally, the performance
of the proposed method is examined in nuclear norm minimization-based
problems.

###  I-B Organization of the Paper

Section ii@ is dedicated to a review of proximal methods which have been
utilized to tackle (1) by projecting the final solution on the desired
feasible set. In Section iii@, we will introduce our algorithm to solve (1),
ORKA, which tackles the problem as a large-scale overdetermined system and
finds the optimal point in the one-bit polyhedron by an accelerated Kaczmarz
approach. Moreover, two new variants of the Kaczmarz algorithms are proposed
that enhance the convergence rate and the computational complexity of these
solvers. To investigate the convergence rate of ORKA, at first, we will
introduce a penalty function in Section iv@ based on the Chernoff bound.
Section v@ discusses an iterative algorithm to achieve optimized time-varying
sampling threshold sequences which benefit the signal recovery process with
enhanced accuracy. As a representative application, in Section vi@, ORKA and
other proposed algorithms will be applied in the context of low-rank matrix
recovery in the form of a nuclear norm minimization problem. Finally, Section
ix@ 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}$ is the
identity matrix of size $N$. $\operatorname{Tr}(.)$ denotes the trace of the
matrix argument.
$\left\langle\mathbf{B}_{1},\mathbf{B}_{2}\right\rangle=\operatorname{Tr}(\mathbf{B}_{1}^{\mathrm{H}}\mathbf{B}_{2})$
is the standard inner product between two matrices. The nuclear norm of a
matrix $\mathbf{B}\in\mathbb{C}^{N_{1}\times N_{2}}$ is denoted
$\left\|\mathbf{B}\right\|_{\star}=\sum^{M}_{i=1}\sigma_{i}$ where $M$ and
$\left\\{\sigma_{i}\right\\}$ are the rank and singular values of
$\mathbf{B}$, respectively. 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^{k}$-norm of a
vector $\mathbf{b}$ is defined as
$\|\mathbf{b}\|^{k}_{k}=\sum_{i}|b|^{k}_{i}$. 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}$. Additionally, the Kronecker product is
denoted as $\mathbf{B}_{1}\otimes\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\\}$. $f\asymp g$ means $f$ and $g$
are asymptotically equal. $\operatorname{Diag}\left\\{\mathbf{b}\right\\}$
denotes a diagonal matrix with $\\{b_{i}\\}$ as its diagonal elements.

##  II Projections On Convex Sets: Dealing With Costly Constraints

To tackle (1), many non-convex and local optimization algorithms have been
developed over the years. Nevertheless, in recent decades, convex programming
formulations via _relaxation_ have come to the fore to approximate _global_
solutions. In the convex framework, various iterative methods have been
proposed to tackle the problem with a Lagrangian formulation such as _Uzawa’s
algorithm_ and the proximal forward-backward splitting method (PFBS) [49, 50,
2]. Moreover, to keep the problem solution inside the constraint set
$\Omega_{c}$, the orthogonal projection $\mathcal{P}_{\Omega_{c}}$ is applied
to solutions in each iteration. This process is briefly explained below.

The Lagrangian for (1) is written as [2],

$\mathcal{L}\left(\mathbf{X},\bm{\lambda}\right)=f\left(\mathbf{X}\right)+\left\langle\bm{\lambda},\mathbf{y}-\mathcal{A}\left(\mathbf{X}\right)\right\rangle,$
(2)

where $\bm{\lambda}\in\mathbb{R}^{n}$. Uzawa’s algorithm aims to find a saddle
point $\left(\mathbf{X}^{\star},\bm{\lambda}^{\star}\right)$, where
$\sup_{\bm{\lambda}}\inf_{\mathbf{X}}\mathcal{L}\left(\mathbf{X},\bm{\lambda}\right)=\inf_{\mathbf{X}}\sup_{\bm{\lambda}}\mathcal{L}\left(\mathbf{X},\bm{\lambda}\right)$,
with the iterative procedure:

$\displaystyle\begin{cases}\mathcal{L}\left(\mathbf{X}^{k},\bm{\lambda}^{k-1}\right)=\min_{\mathbf{X}}\mathcal{L}\left(\mathbf{X},\bm{\lambda}^{k-1}\right),\\\
\bm{\lambda}^{k}=\mathcal{P}_{\Omega_{c}}\left(\bm{\lambda}^{k-1}+\alpha_{k}\left(\mathbf{y}-\mathcal{A}\left(\mathbf{X}^{k}\right)\right)\right),\end{cases}$
(3)

where $\alpha_{k}$ is the step size. This iterative steps can be rewritten as

$\displaystyle\begin{cases}\mathbf{X}^{k}=\operatorname{Prox}_{f}\left(\mathcal{A}^{\star}\left(\bm{\lambda}^{k-1}\right)\right),\\\
\bm{\lambda}^{k}=\mathcal{P}_{\Omega_{c}}\left(\bm{\lambda}^{k-1}+\alpha_{k}\left(\mathbf{y}-\mathcal{A}\left(\mathbf{X}^{k}\right)\right)\right),\end{cases}$
(4)

where $\operatorname{Prox}_{f}$ is the proximal operator minimizing the
Lagrangian function, and $\mathcal{A}^{\star}$ is the adjoint of
$\mathcal{A}$.

Since every linear equation can be reformulated in standard form, we recast
$\mathcal{A}\left(\mathbf{X}\right)=\mathbf{y}$ as
$\mathbf{A}\mathbf{x}=\mathbf{y}$, where $\mathbf{A}\in\mathbb{C}^{n\times
n_{1}n_{2}}$ is a matrix version of the operator $\mathcal{A}$, and
$\mathbf{x}=\operatorname{vec}\left(\mathbf{X}\right)$[15]. The optimization
problem (1) is equivalently given by [2, 1]

$\displaystyle\min_{\mathbf{X}}$ $\displaystyle
g\left(\mathbf{X}\right)=\frac{1}{2}\left\|\mathbf{y}-\mathbf{A}\operatorname{vec}\left(\mathbf{X}\right)\right\|^{2}_{2}+\lambda
f(\mathbf{X})$ (5) s.t. $\displaystyle\mathbf{X}\in\Omega_{c}.$

To solve this problem, instead of using proximal methods, a projected gradient
method such as Nesterov iterative approach may be utilized, i.e.,
$\mathbf{X}^{k}=\mathcal{P}_{\Omega_{c}}\left(\mathbf{X}^{k-1}-\alpha_{k}\nabla
g\left(\mathbf{X}^{k-1}\right)\right)$.

Famous examples for $\operatorname{Prox}_{f}$ and $\mathcal{P}_{\Omega_{c}}$,
are the singular value thresholding operator (SVT) and the semi-definite
orthogonal projector, respectively. SVT is useful when
$f\left(\mathbf{X}\right)=\left\|\mathbf{X}\right\|_{\star}$, mathematically
defined as [2]:

$\displaystyle\mathcal{D}_{\delta}=\mathbf{U}\operatorname{Diag}\left\\{\left(\sigma_{k}-\delta\right)^{+}\right\\}\mathbf{V}^{\top},$
(6)

where $\mathbf{U}$ and $\mathbf{V}$ are unitary matrices from singular value
decomposition (SVD), and $\\{\sigma_{k}\\}$ are the singular values.
Furthermore, the semi-definite projector emerges in semi-definite programming
where the convex constraint set is a positive semi-definite (PSD) matrix. It
compares eigenvalues of the solution in each iteration with zero or a fixed
threshold [1], i.e.,

$\displaystyle\mathcal{P}_{\Omega_{c}}=\mathbf{U}\operatorname{Diag}\left\\{\left(\lambda_{k}-\delta\right)^{+}\right\\}\mathbf{U}^{\top},$
(7)

where $U$ is the unitary matrix coming from the Schur decomposition. In the
case of both operators, the approximate solution should be projected onto a
feasible convex set at each iteration via recovering all singular values and
eigenvalues and comparing their smaller elements with a threshold, which is
quite expensive [1].

An interesting alternative to enforcing the feasible set
$\mathcal{F}_{\mathbf{X}}=\left\\{\operatorname{Prox_{f}\cap\Omega_{c}}\right\\}$
in (1) emerges when one increases the number of samples $n$, and solves the
overdetermined linear system of equations with $n\geq n_{1}n_{2}$. In this
sample abundance regimen, the linear constraint
$\mathcal{A}\left(\mathbf{X}\right)=\mathbf{y}$ 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 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 [16, 27].

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 in forming the feasible polyhedron. Therefore, by
increasing the number of samples, a finite-volume space may be created inside
$\mathcal{F}_{\mathbf{X}}$ with shrinking size; making projections on
$\Omega_{c}$ redundant. _From a practical point of view, one-bit sampling is
done efficiently at a very high rate with a significantly lower cost compared
to its high-resolution counterpart. It has been examined in[16] that even
though only partial information is made available to one-bit signal processing
algorithms, they can achieve acceptable recovery performance with less
complexity compared to the high-resolution scenario. Thus, it is both
practical and necessary to study the ground-breaking opportunities that emerge
in the context of the wide array of problems formulated as (1) due to the
availability of a large number of one-bit samples._

##  III Proposed Algorithm

In this section, at first we begin by presenting a summarized review of
randomized Kaczmarz algorithms. Then, we propose a novel Kaczmarz method
variant formulated based on the _sampling Kaczmarz-Motzkin algorithm_ (SKM)
and a _preconditioning approach_. One-bit sampling via time-varying thresholds
will be combined with the proposed randomized Kaczmarz method to create highly
overdetermined linear inequalities. This paves the way for the recovery of the
desired signal $\mathbf{X}^{\star}$ in (1) without solving the original
optimization problem; merely by tacking accounts of its linear constraints. We
name our algorithm _O_ ne-bit aided _R_ andomized _K_ aczmarz _A_ lgorithm
(ORKA). Due to the block structure of the linear feasibility in ORKA, we will
propose a block-based Kaczmarz algorithm accordingly.

###  III-A Randomized Kaczmarz Algorithm (RKA)

The randomized Kaczmarz algorithm (RKA) is a _sub-conjugate gradient method_
to solve a linear feasibility problem, i.e,
$\mathbf{C}\mathbf{x}\preceq\mathbf{b}$ where $\mathbf{C}$ is a ${m\times n}$
matrix with $m>n$ [45, 44]. Conjugate-gradient methods immediately turn the
mentioned inequality to an equality in the following form:

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

and then, approach the solution by the same process as used for systems of
equations. Without any loss of generality, consider (8) to be a polyhedron:

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

where the disjoint index sets $I_{\leq}$ and $I_{=}$ partition our sample
index set $\mathcal{J}$, and $\\{\mathbf{c}_{j}\\}$ denote the rows of
$\mathbf{C}$. Based on this problem, the projection coefficient $\beta_{i}$ of
the RKA is defined as [45, 47, 51]:

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

Also, the unknown column vector $\mathbf{x}$ is iteratively updated as

$\mathbf{x}_{i+1}=\mathbf{x}_{i}-\frac{\beta_{i}}{\left\|\mathbf{c}_{j}\right\|^{2}_{2}}\mathbf{c}^{\mathrm{H}}_{j},$
(11)

where, at each iteration $i$, the index $j$ is chosen independently at random
from the set $\mathcal{J}$, following the distribution

$P\\{j=k\\}=\frac{\left\|\mathbf{c}_{k}\right\|^{2}_{2}}{\|\mathbf{C}\|_{\mathrm{F}}^{2}}.$
(12)

If the system (9) is consistent and its feasible region is nonempty, RKA
converges linearly in expectation [44, 45]:

$\mathbb{E}\left\\{\hbar\left(\mathbf{x}_{i},\mathbf{x}^{\star}\right)\right\\}\leq
q^{i}\leavevmode\nobreak\
\hbar\left(\mathbf{x}_{0},\mathbf{x}^{\star}\right),$ (13)

where
$\hbar\left(\mathbf{x}_{i},\mathbf{x}^{\star}\right)=\left\|\mathbf{x}_{i}-\mathbf{x}^{\star}\right\|_{2}^{2}$
is the distance function between two points in the space, $\mathbf{x}^{\star}$
is a desired point, $i$ is the number of required iterations for RKA, and
$q\in\left(0,1\right)$ is given as

$q=1-\frac{1}{\kappa^{2}\left(\mathbf{C}\right)},$ (14)

with
$\kappa\left(\mathbf{C}\right)=\|\mathbf{C}\|_{\mathrm{F}}\|\mathbf{C}^{\dagger}\|_{2}$
denoting the scaled condition number.

###  III-B Sampling Kaczmarz-Motzkin Algorithm (SKM)

The SKM combines the ideas of both the RKA and the Motzkin method. Its
generalized convergence theorem, and a validation of feasibility, which has
been formulated based on the convergence analysis of RKA and sampling Motzkin
method for solving linear feasibility problem have been fully explored in
[48].

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
$\tau_{i}$), uniformly at random out of $m$ rows of the constraint matrix
$\mathbf{C}$. Then, out of these $\gamma$ rows, the row with maximum positive
residual is selected. Finally, the solution is updated as [48, 52]:

$\mathbf{x}_{i+1}=\mathbf{x}_{i}-\lambda_{i}\frac{\beta_{i}}{\left\|\mathbf{c}_{j^{\star}}\right\|^{2}_{2}}\mathbf{c}^{\mathrm{H}}_{j^{\star}},$
(15)

where $j^{\star}=\operatorname{argmax}\leavevmode\nobreak\
\left\\{\left(\mathbf{c}_{j}\mathbf{x}_{i}-b_{j}\right)^{+}\right\\},\leavevmode\nobreak\
j\in\tau_{i}$, and $\lambda_{i}$ is a relaxation parameter which for
consistent systems must satisfy [44],

$0\leq\lim_{i\to\infty}\inf\lambda_{i}\leq\lim_{i\to\infty}\sup\lambda_{i}<2,$
(16)

to ensure convergence. The convergence bound for SKM is given by

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

In the case where the constraint matrix is normalized, i.e.
$\|\mathbf{c}_{j}\|^{2}_{2}=1$, $s_{i}$ is the number of satisfied constraints
after iteration $i$, and $L_{i}=\max\left\\{m-s_{i},m-\gamma\right\\}$, for
the $(i+1)$th iteration we have [48],

$\mathbb{E}\left\\{\hbar\left(\mathbf{x}_{i},\mathbf{x}^{\star}\right)\right\\}\leq\left(1-\frac{\sigma^{2}_{min}\left(2\lambda_{i}-\lambda^{2}_{i}\right)}{V_{i}}\right)^{i}\leavevmode\nobreak\
\hbar\left(\mathbf{x}_{0},\mathbf{x}^{\star}\right).$ (18)

This recovery error bound is tighter than (17).

###  III-C Our Contribution: Preconditioned SKM (PrSKM)

According to the convergence rate formula of RKA, if we can reduce the value
of the scaled condition number, the convergence is accelerated, and the upper
bound of the recovery error
$\mathbb{E}\left\\{\left\|\mathbf{x}_{i}-\mathbf{x}^{\star}\right\|_{2}^{2}\right\\}$
decreases. Moreover, by having a lower value of $q$, a lower number of
iterations is required to achieve a specific recovery error bound, usually
considered to be the algorithm’s termination criterion. Consequently, let $I$
is the number of iterations, the computational cost of RKA which behaves as
$\mathcal{O}\left(In\right)$, is diminished as well. To make this happen, one
can start from reducing $q=1-\frac{1}{\kappa^{2}\left(\mathbf{C}\right)}$
which occurs when the scaled condition number $\kappa\left(\mathbf{C}\right)$
is diminished; a condition that can be satisfied by considering the following
theorem.

######  Theorem 1.

The infimum scaled condition number of a matrix
$\mathbf{C}\in\mathbb{R}^{m\times n}$ is given by

$\inf_{\mathbf{C}}\kappa\left(\mathbf{C}\right)=\sqrt{n},$ (19)

which is achieved if and only if $\mathbf{C}$ is of the form
$\mathbf{C}=\alpha\mathbf{U}$, where $\mathbf{U}$ is an orthonormal matrix and
$\alpha\in\mathbb{R}$ is a scalar.

###### Proof.

The condition number of the matrix $\mathbf{C}$ is defined as
$\varrho\left(\mathbf{C}\right)=\frac{\sigma_{max}}{\sigma_{min}}$, where
$\sigma_{max}$ and $\sigma_{min}$ are its minimum and maximum singular values,
respectively [53]. The scaled condition number can be written as
$\kappa\left(\mathbf{C}\right)=\frac{\|\mathbf{C}\|_{\mathrm{F}}}{\sigma_{min}}$.
Therefore, the scaled condition number has the following relation with
$\varrho\left(\mathbf{C}\right)$:

$\kappa\left(\mathbf{C}\right)=\frac{\|\mathbf{C}\|_{\mathrm{F}}}{\sigma_{max}}\varrho\left(\mathbf{C}\right).$
(20)

Furthermore, the condition number can be considered to be an upper bound for
the scaled condition number as well based on the readily-known inequality
relation between norm-$2$ and the Frobenius norm [53]:

$\displaystyle\|\mathbf{C}\|_{\mathrm{F}}$
$\displaystyle\leq\sqrt{n}\|\mathbf{C}\|_{\mathrm{2}},$ (21)
$\displaystyle\frac{\|\mathbf{C}\|_{\mathrm{F}}}{\sigma_{\min}}$
$\displaystyle\leq\sqrt{n}\frac{\|\mathbf{C}\|_{\mathrm{2}}}{\sigma_{\min}},$

or equivalently,

$\displaystyle\kappa\left(\mathbf{C}\right)\leq\sqrt{n}\varrho\left(\mathbf{C}\right).$
(22)

Thus, lowering $\varrho\left(\mathbf{C}\right)$ generally leads to a
decreasing scaled condition number. Additionally, the lowest possible value
for $\varrho$ is $1$ which is achieved for scaled unitary matrices
$\mathbf{U}$ as if we let $\mathbf{S}=\alpha\mathbf{U}$, and
$\mathbf{O}=\mathbf{S}^{\top}\mathbf{S}=\alpha^{2}\mathbf{I}_{n}$, then,
$\sigma_{min}=\sigma_{max}=\alpha$, and $\varrho=1$. Vice versa, if
$\varrho=1$, it means $\sigma_{min}=\sigma_{max}$ which leads to a diagonal
matrix $\mathbf{O}=\alpha^{2}\mathbf{I}_{n}$. It is straightforward to verify
that the decomposition of $\mathbf{O}$ results in an $\mathbf{S}$ that is a
scaled-version of an orthonormal matrix. As a result, the lowest achievable
upper bound for the scaled condition number is obtained as
$\kappa\left(\mathbf{C}\right)\leq\sqrt{n}$, and according to (20),
$\kappa\left(\mathbf{C}\right)=\frac{\alpha\|\mathbf{U}\|_{\mathrm{F}}}{\alpha}=\sqrt{n}$.
∎

Accordingly, it would be enough to make our matrix $\mathbf{C}$ unitary by a
process which is referred to as the _preconditioning method_. In
preconditioning, the linear feasibility is rewritten as

$\mathbf{C}\mathbf{M}\mathbf{z}\preceq\mathbf{b},$ (23)

where $\mathbf{M}$ is the preconditioner and
$\mathbf{x}=\mathbf{M}\mathbf{z}$. The straightforward way to approach this
task is to use QR decomposition where the constraint matrix is decomposed as
$\mathbf{C}=\mathbf{Q}_{c}\mathbf{R}_{c}$, with unitary
$\mathbf{Q}_{c}\in\mathbb{R}^{m\times n}$, and
$\mathbf{R}_{c}\in\mathbb{R}^{n\times n}$ is an upper triangular matrix,
leading to

$\mathbf{Q}_{c}=\mathbf{C}\mathbf{R}^{-1}.$ (24)

Thus, the good choice for the preconditioner is
$\mathbf{M}=\mathbf{R}^{-1}_{c}$. To find the desired point
$\mathbf{z}^{\star}$, the SKM is selected in order to apply to the linear
feasibility (23), then the desired signal $\mathbf{x}^{\star}$ is obtained
from $\mathbf{x}^{\star}=\mathbf{R}^{-1}_{c}\mathbf{z}^{\star}$. We refer to
this method _Pr_ econditioned _SKM_ (PrSKM).

######  Proposition 1 (PrSKM).

The proposed algorithm, PrSKM, can be summarized as follows:

  1. 1.

Calculate the QR decomposition of the constraint matrix $\mathbf{C}$ to obtain
the preconditioner $\mathbf{M}$.

  2. 2.

Using the change of variables, $\mathbf{x}=\mathbf{M}\mathbf{z}$, obtain
$\mathbf{H}\mathbf{z}\preceq\mathbf{b}$, where
$\mathbf{H}=\mathbf{C}\mathbf{M}$ and $\mathbf{M}=\mathbf{R}^{-1}_{c}$.

  3. 3.

Choose a sample set of $\gamma$ constraints (denoted as $\tau_{i}$) uniformly
at random from the rows of $\mathbf{H}$.

  4. 4.

From these $\gamma$ constraints, choose
$j^{\star}=\operatorname{argmax}\leavevmode\nobreak\
\left\\{\left(\mathbf{h}_{j}\mathbf{z}_{i}-b_{j}\right)^{+}\right\\},\leavevmode\nobreak\
j\in\tau_{i}$ where $\mathbf{h}_{j}$ is the $j$th row of $\mathbf{H}$.

  5. 5.

Update the solution via the iterations
$\mathbf{z}_{i+1}=\mathbf{z}_{i}-\lambda_{i}\frac{\left(\mathbf{h}_{j^{\star}}\mathbf{z}_{i}-b_{j^{\star}}\right)^{+}}{\left\|\mathbf{h}_{j^{\star}}\right\|^{2}_{2}}\mathbf{h}^{\mathrm{H}}_{j^{\star}}$
until convergence.

  6. 6.

Recover the desired signal from the final solution of SKM as
$\mathbf{x}^{\star}=\mathbf{M}\mathbf{z}^{\star}$.

The scaled condition number of PrSKM is obtained as
$\kappa\left(\mathbf{H}\right)=\sqrt{n}$, which implies $q=\frac{n-1}{n}$.

###  III-D One-Bit Polyhedron

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

$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),$
(25)

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 $\mathrm{T}$. Let $r_{k}$ denote the quantized version
of $y[k]$ with the formulation

$r_{k}=Q(y_{k}),$ (26)

where $Q$ denotes the quantization effect. Consider a non-zero time-varying
Gaussian threshold $\bm{\uptau}=\left[\uptau_{k}\right]$ with 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, (26) is simply written as

$r_{k}=\operatorname{sgn}\left(y_{k}-\tau_{k}\right),$ (27)

where $\operatorname{sgn}(.)$ is the sign function.

Algorithm 1 Block SKM Algorithm

1:_Input:_ Matrix $\mathbf{B}\in\mathbb{R}^{mn\times d}$ where
$\mathbf{B}=\left[\begin{array}[]{c|c|c}\mathbf{B}^{\top}_{1}&\cdots&\mathbf{B}^{\top}_{m}\end{array}\right]^{\top}$,
right-hand side $\mathbf{b}$ with dimension $mn$, initial value of
$\mathbf{x}_{0}$ with dimension $n$, convergence tolerance $\epsilon>0$, and
relaxation parameter $\lambda_{i}$.

2:_Output:_ An estimate $\mathbf{x}$ for the solution to the linear
feasibility problem $\mathbf{B}\mathbf{x}\preceq\mathbf{b}$.

3:Initiate the following loop by setting $i=0$.

4:while
$\left\|\left(\mathbf{B}\mathbf{x}_{i}-\mathbf{b}\right)^{+}\right\|_{2}\leq\epsilon$
do

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

6: 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}$).

7: Select the first $k^{\prime}<d$ element of $\mathbf{e}^{\prime}$ and
construct the problem
$\mathbf{B}_{j}^{\prime}\mathbf{x}\preceq\mathbf{b}_{j}^{\prime}$, where
$\mathbf{B}_{j}^{\prime}\in\mathbb{R}^{k^{\prime}\times d}$ and
$\mathbf{b}_{j}^{\prime}\in\mathbb{R}^{k^{\prime}\times 1}$.

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

9: Update the solution $\mathbf{x}_{i+1}$ as:
$\mathbf{x}_{i+1}\leftarrow\mathbf{x}_{i}-\lambda_{i}\mathbf{B}_{j}^{\prime\dagger}\left(\mathbf{B}_{j}^{\prime}\mathbf{x}-\mathbf{b}_{j}^{\prime}\right)^{+}.$

10: Increase $i$ by one.

11:end while

The information gathered through the one-bit sampling with time-varying
thresholds presented here may be formulated in terms of an overdetermined
linear system of inequalities. In Eq. (27),

$r_{k}=\begin{cases}+1&y_{k}>\uptau_{k},\\\ -1&y_{k}<\uptau_{k}.\end{cases}$
(28)

Therefore, one can formulate the geometric location of the signal as

$r_{k}\left(y_{k}-\uptau_{k}\right)\geq 0.$ (29)

Let $\mathbf{y}=[y_{k}]$ and $\mathbf{r}=[r_{k}]$. Then, the vectorized
representation of (29) is

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

or equivalently

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

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

$\displaystyle\begin{cases}\bm{\Omega}^{(1)}\mathbf{y}\succeq\mathbf{r}^{(1)}\odot\bm{\uptau}^{(1)}&\quad\mathbf{r}^{(1)}=\operatorname{sgn}\left(\mathbf{y}-\bm{\uptau}^{(1)}\right),\\\
\bm{\Omega}^{(2)}\mathbf{y}\succeq\mathbf{r}^{(2)}\odot\bm{\uptau}^{(2)}&\quad\mathbf{r}^{(2)}=\operatorname{sgn}\left(\mathbf{y}-\bm{\uptau}^{(2)}\right),\\\
\vdots&\quad\vdots\\\
\bm{\Omega}^{(m)}\mathbf{y}\succeq\mathbf{r}^{(m)}\odot\bm{\uptau}^{(m)}&\quad\mathbf{r}^{(m)}=\operatorname{sgn}\left(\mathbf{y}-\bm{\uptau}^{(m)}\right),\end{cases}$
(32)

where
$\bm{\Omega}^{(\ell)}=\operatorname{diag}\left\\{\mathbf{r}^{(\ell)}\right\\}$.
In Eq. (32), we have $m$ linear system of inequalities which can also be
merged in one inequality as

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

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

$\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}^{mn\times
n}.$ (34)

Assuming a large number of samples which is a common situation in one-bit
sampling scenarios, hereafter, we consider (33) as an overdetermined linear
system of inequalities associated with the sampling scheme presented in (27).
The inequality (33) can be recast by a polyhedron as

$\displaystyle\mathcal{P}_{\mathbf{y}}=\left\\{\mathbf{y}\mid\tilde{\bm{\Omega}}\mathbf{y}\succeq\operatorname{vec}\left(\mathbf{R}\right)\odot\operatorname{vec}\left(\bm{\Gamma}\right)\right\\},$
(35)

which we refer to as the _one-bit polyhedron_.

###  III-E ORKA: Towards Circumventing Costly Constraints

If one applies one-bit sampling with time-varying sampling thresholds to the
measurement $\mathbf{y}\in\mathbb{R}^{n}$ from (1) following the process
defined in Subsection III-D, the arising inequality system is simply given by

$\tilde{\bm{\Omega}}\mathbf{A}\mathbf{x}\succeq\operatorname{vec}\left(\mathbf{R}\right)\odot\operatorname{vec}\left(\bm{\Gamma}\right),$
(36)

where $\mathbf{A}\mathbf{x}=\mathbf{y}$, and
$\mathbf{x}=\operatorname{vec}\left(\mathbf{X}\right)$. Consequently, the one-
bit polyhedron for this problem is obtained as

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

where $\mathbf{P}=\tilde{\bm{\Omega}}\mathbf{A}$.

By taking advantage of one-bit sampling, in the asymptotic scenario of with
sample abundance, the space restricted by the one-bit polyhedron
$\mathcal{P}_{\mathbf{x}}$, _shrinks_ to become contained inside the feasible
set $\mathcal{F}_{\mathbf{X}}$. Note that this shrinking space always contains
the global minima, with a volume that is diminished with an increased sample
size. As a result, instead of using proximal operators and orthogonal
projectors, it is enough to find the desired signal $\mathbf{x}^{\star}$ in
(37). To do so, one can use the PrSKM algorithm proposed in Section III-C. It
is worth noting that the PrSKM is a _row-based_ algorithm, where at each
iteration, the row index is chosen independently at random. However, the
matrix $\mathbf{P}$ in (37) has a block structure with the following
formulation

$\mathbf{P}=\left[\begin{array}[]{c|c|c}\mathbf{A}^{\top}\bm{\Omega}^{(1)}&\cdots&\mathbf{A}^{\top}\bm{\Omega}^{(m)}\end{array}\right]^{\top},\quad\mathbf{P}\in\mathbb{R}^{mn\times
d},$ (38)

where $d=n_{1}n_{2}$. Therefore, it is useful to investigate the block-based
RKA methods to find the desired signal in $\mathcal{P}_{\mathbf{x}}$ for
further efficiency enhancement. Our proposed algorithm, _Block SKM_ , is
described as follows.

Algorithm 2 Architecture of ORKA.

1:_Input:_ The measurement vector $\mathbf{y}$ obtained as
$\mathcal{A}\left(\mathbf{X}\right)=\mathbf{y}$ from (1), $m$ sequences of
time-varying sampling thresholds generated as
$\left\\{\bm{\uptau}^{(\ell)}\sim\mathcal{N}\left(\mathbf{0},\mathbf{I}\right);\ell\in\mathcal{L}\right\\}$.

2:_Output:_ Recovered optimal signal
$\mathbf{x}^{\star}=\operatorname{vec}\left(\mathbf{X}^{\star}\right)$.

3:Apply one-bit sampling on $\mathbf{y}$ and generate sequences of one-bit
measurements from:
$\mathbf{r}^{(\ell)}\leftarrow\operatorname{sgn}\left(\mathbf{x}-\bm{\uptau}^{(\ell)}\right),\quad\ell\in\mathcal{L}.$

4:Construct a linear feasibility region from the one-bit sampled data as:
$\bm{\Omega}^{(\ell)}\mathbf{y}\succeq\mathbf{r}^{(\ell)}\odot\bm{\uptau}^{(\ell)},\quad\ell\in\mathcal{L}.$

5:Define a highly-overdetermined system, the one-bit polyhedron, based on
obtained inequalities:
$\mathcal{P}_{\mathbf{x}}=\left\\{\mathbf{x}\mid\mathbf{P}\mathbf{x}\succeq\operatorname{vec}\left(\mathbf{R}\right)\odot\operatorname{vec}\left(\bm{\Gamma}\right)\right\\},$
where $\mathbf{R}$ and $\bm{\Gamma}$ are matrices with
$\left\\{\mathbf{r}^{(\ell)}\right\\}$ and
$\left\\{\bm{\uptau}^{(\ell)}\right\\}$ representing their columns,
respectively.

6:Employ RKA variants (PrSKM or Block SKM) to recover $\mathbf{X}^{\star}$
within the one-bit polyhedron.

######  Proposition 2 (Block SKM).

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}\end{array}\right]^{\top}$,
and
$\mathbf{b}=\left[\begin{array}[]{c|c|c}\mathbf{b}^{\top}_{1}&\cdots&\mathbf{b}^{\top}_{m}\end{array}\right]^{\top}$.
The proposed algorithm for feasible signal recovery, Block SKM, can be
summarized as follows:

  1. 1.

Choose a block $\mathbf{B}_{j}$ uniformly at random with the probability
$P\\{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 Section III-B, to have an accelerated
convergence.

  4. 4.

Select the first $k^{\prime}<d$ 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 d}$ and
$\mathbf{b}_{j}^{\prime}\in\mathbb{R}^{k^{\prime}\times 1}$. The reason behind
choosing $k^{\prime}<d$ 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}>d$, 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}$,
$\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 [55,
46] which takes advantage of the efficient matrix-vector multiplication, thus
giving the method a significant reduction in computational cost [47].

The steps of the block SKM and ORKA are summarized in Algorithm 1 and
Algorithm 2, respectively. To examine the performance of the block SKM, we
compare it to the PrSKM, SKM and RKA.

(a)

(b)

Figure 1: Comparing the NMSE recovery performance of the two proposed Kaczmarz
algorithms, namely the PrSKM and the block SKM, with that of SKM and RKA for:
(a) a linear equation system, (b) a linear inequality system.

###  III-F Comparing RKA, SKM, PrSKM and Block SKM

In this section, we numerically compare the RKA, SKM, PrSKM, and Block SKM in
linear systems of equalities as well as those formed by inequalities.

Linear feasibility of equalities: Herein, we consider a block linear system of
equalities $\mathbf{A}\mathbf{x}=\mathbf{b}$, where
$\mathbf{A}=\left[\begin{array}[]{c|c|c}\mathbf{A}_{1}^{\top}&\cdots&\mathbf{A}_{100}^{\top}\end{array}\right]^{\top}$,
$\mathbf{A}_{i}\in\mathbb{R}^{10\times 10}$, $\mathbf{x}\in\mathbb{R}^{10}$,
and $\mathbf{b}\in\mathbb{R}^{1000}$. Each row of $\mathbf{A}_{i}$ is
generated as
$\mathbf{a}_{j}^{i}\sim\mathcal{N}\left(\mathbf{0},\mathbf{I}_{10}\right)$.
Also, the unknown signal $\mathbf{x}$ is generated as
$\mathbf{x}\sim\mathcal{N}\left(\mathbf{0},\mathbf{I}_{10}\right)$. The
normalized mean square error (NMSE) is defined as

$\operatorname{NMSE}\triangleq\frac{\left\|\mathbf{x}^{\star}-\bar{\mathbf{x}}\right\|_{2}^{2}}{\left\|\mathbf{x}^{\star}\right\|_{2}^{2}},$
(39)

where $\mathbf{x}^{\star}$ and $\bar{\mathbf{x}}$ denote the true discretized
signal and its recovered version, respectively.

Fig. 1 illustrates the performance of RKA, SKM, PrSKM, and Block SKM in the
recovery of $\mathbf{x}$ from the system $\mathbf{A}\mathbf{x}=\mathbf{b}$
with NMSE results. As can be observed, the Block SKM outperforms the other
three approaches in the recovery task. Also, it can be seen that the PrSKM has
a better recovery performance compared to that of the RKA and the SKM.

Linear feasibility of inequalities: We utilize ORKA to make a linear equation
$\mathbf{B}\mathbf{x}=\mathbf{y}$ linear inequalities system, where the number
of time-varying sampling threshold sequences is $m=40$,
$\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}\in\sim\mathcal{N}\left(\mathbf{0},\mathbf{I}_{10}\right)$. Each
time-varying sampling threshold sequence $\bm{\uptau}^{(\ell)}$ is considered
to have the distribution
$\bm{\uptau}^{(\ell)}\sim\mathcal{N}\left(\mathbf{0},\mathbf{I}_{10}\right)$.
The performance of the RKA, SKM, PrSKM, and Block SKM is illustrated in Fig.
1. Similar to the linear feasibility of equalities, it can be seen that the
Block SKM has a better accuracy in the recovery of the desired signal
$\mathbf{x}$ in the one-bit polyhedron (37) compared to the other three
approaches. The NMSE results in Fig. 1 are averaged over $15$ experiments.

##  IV Probabilistic Effect of Sample Abundance In ORKA

An integral part of our proposed recovery algorithm is RKA, whose recovery
error was readily given in (17). As shown in [44], the convergence rate of RKA
does not depend on the number of equations in the system. We will show that
the convergence rate of ORKA for linear feasibility is the same as RKA.
Nevertheless, when we face a non-linear constraint in our problem, as is
generally the case in (1), it is desirable are made redundant by using the
opportunity of having a large number of samples; as typically provided via
one-bit sampling. In such a case, The offered convergence rate appears to be
insufficient since we must have enough number of samples to fulfill costly
constraints. So, an extra term as a _penalty_ must be considered to present
the importance of sample size in our algorithm [16].

By adding more inequality constraints in (37) as a result of extra one-bit
samples, the shrinkage of the said polyhedron will put a downward pressure on
the distance between the desired signal $\mathbf{x}^{\star}$ and its
surrounding hyperplanes, each presenting an informative measurement that will
shrink the feasibility space. We will show that by judicious sampling, the
average of these distances will be bounded, which may be considered to be a
finite-volume space created around $\mathbf{x}^{\star}$. Moreover, as a result
of using an overdetermined linear system of inequalities, the convergence of
the RKA is guaranteed [16, 45, 46, 48].

###  IV-A Recovery Error Upper Bound for ORKA

As the scaled condition number is the central parameter governing the recovery
error of Kaczmarz algorithms, we will evaluate it for ORKA-created matrix
$\mathbf{P}$ in the following, starting with $\sigma_{min}$. The singular
values of $\mathbf{P}$ may be determined based on the following theorem, which
thus unveils the value of $\sigma_{min}$.

######  Theorem 2.

In ORKA, the rank of $\mathbf{P}$ is equal to that of the constraint matrix
$\mathbf{A}$, and its singular values are given by

$\left\\{\sigma_{i}\right\\}=\sqrt{m}\left\\{\sigma_{i\mathbf{A}}\right\\},$
(40)

where $\left\\{\sigma_{i\mathbf{A}}\right\\}$ are singular values of
$\mathbf{A}$, and $m$ is the number of time-varying sampling threshold
sequences. Moreover, the scaled condition number of the ORKA-created matrix
$\mathbf{P}$ is equal to that of the constraint matrix $\mathbf{A}$:

$\kappa\left(\mathbf{P}\right)=\kappa\left(\mathbf{A}\right).$ (41)

###### Proof.

To obtain the singular values of $\mathbf{P}$, the matrix
$\mathbf{W}=\mathbf{P}^{\top}\mathbf{P}$ is computed as

$\displaystyle\mathbf{W}$ $\displaystyle=\mathbf{P}^{\top}\mathbf{P},$ (42)
$\displaystyle=\left[\mathbf{A}^{\top}\bm{\Omega}^{(1)}\vdots\mathbf{A}^{\top}\bm{\Omega}^{(2)}\vdots\cdots\vdots\mathbf{A}^{\top}\bm{\Omega}^{(m)}\right]\left[\begin{array}[]{c}\bm{\Omega}^{(1)}\mathbf{A}\\\
\cdots\\\ \bm{\Omega}^{(2)}\mathbf{A}\\\ \vdots\\\ \cdots\\\
\bm{\Omega}^{(m)}\mathbf{A}\end{array}\right],$
$\displaystyle=\mathbf{A}^{\top}\bm{\Omega}^{(1)}\bm{\Omega}^{(1)}\mathbf{A}+\cdots+\mathbf{A}^{\top}\bm{\Omega}^{(m)}\bm{\Omega}^{(m)}\mathbf{A},$
$\displaystyle=m\mathbf{A}^{\top}\mathbf{I}\mathbf{A}=m\mathbf{A}^{\top}\mathbf{A},$

which means the singular values of $\mathbf{P}$ are
$\left\\{\sigma_{i}\right\\}=\sqrt{m}\left\\{\sigma_{i\mathbf{A}}\right\\}$.

Also, the Frobenius norm of $\mathbf{P}$ is obtained as

$\displaystyle\|\mathbf{P}\|^{2}_{\mathrm{F}}$
$\displaystyle=\operatorname{Tr}\left(\mathbf{P}^{\top}\mathbf{P}\right),$
(43)
$\displaystyle=\operatorname{Tr}\left(m\mathbf{A}^{\top}\mathbf{A}\right)=m\|\mathbf{A}\|^{2}_{\mathrm{F}}.$

Consequently, the scaled condition number is independent of the number of
time-varying thresholds sequences. It follows that
$\kappa\left(\mathbf{P}\right)=\kappa\left(\mathbf{A}\right)$. ∎

######  Corollary 1.

For $\mathbf{A}=\mathbf{I}$, corresponding to the one-bit sampled signal
recovery problem, the scaled condition number of
$\bm{\Omega}=\left[\bm{\Omega}^{(1)}\vdots\bm{\Omega}^{(2)}\vdots\cdots\vdots\bm{\Omega}^{(m)}\right]$
is $\kappa\left(\Omega\right)=\sqrt{n}$, which is the infimum of the scaled
condition number as shown in Theorem 1.

Note that the convergence bound (17) for ORKA is independent of the number of
time-varying sampling threshold sequences $m$, and it cannot take into account
the effect of an increasing number of time-varying threshold sequences.
Inspired by [16], we augment (17) with a sample size-dependent penalty
function to make it useful in a sample abundance scenario:

######  Proposition 3 (Convergence rate of ORKA).

In the proposed recovery approach, it is deemed necessary to have a sufficient
number of samples (inequalities) in order to guarantee a finite-volume
feasible region and a bounded recovery error. Once our search area is located
inside $\mathcal{F}_{\mathbf{X}}$, we may effectively deploy (17) for the
convergence rate. The convergence rate of the Kaczmarz variants is useful when
we have a linear feasibility problem. On the other hand, in (1), the main
constraints are non-linear and they may be considered to be redundant by
deploying enough samples [16]. A sample size-aware convergence rate for ORKA
may be formulated as:

$\displaystyle\mathbb{E}\left\\{\hbar\left(\mathbf{x}_{i},\mathbf{x}^{\star}\right)\right\\}\leq\left(1-\frac{2\lambda_{i}-\lambda^{2}_{i}}{\kappa^{2}\left(\mathbf{A}\right)}\right)^{i}\hbar\left(\mathbf{x}_{0},\mathbf{x}^{\star}\right)+\Upsilon\left(m\right),$
(44)

where $\Upsilon(.)$ is an asymptotically decreasing function, such that if the
number of time-varying threshold sequences is enough for the one-bit
polyhedron to fit inside $\mathcal{F}_{\mathbf{X}}$, the sample size-dependent
penalty function $\Upsilon\left(m\right)$ approaches zero.

To propose an appropriate sample size-dependent penalty function, we will
utilize the first theorem in [16], which studies the possibility of creating a
finite-volume space around the optimal signal.

###  IV-B Sample Size-Dependent Penalty Function via Moment Generating
Functions

We investigate the convergence of ORKA through a probabilistic lens. To do so,
define the distance between the optimal point $\mathbf{x}^{\star}$ and the
$j$-th hyperplane presented in (37) as

$\displaystyle d_{j}\left(\mathbf{x}^{\star},\bm{\uptau}^{(\ell)}\right)$
$\displaystyle=\left|r_{j}\odot\left(\mathbf{a}_{j}\mathbf{x}^{\star}-\tau^{(\ell)}_{j}\right)\right|^{2},\quad
j\in\left\\{1,\cdots,m^{\prime}\right\\},$ (45)

where $r_{j}\odot\mathbf{a}_{j}$ is the $j$-th row of $\mathbf{P}$,
$\mathbf{a}_{j}$ is the $j$-th row of $\mathbf{A}$, and
$\mathbf{a}_{j}=\mathbf{a}_{j+n}$. It is easy to observe that by generally
reducing the distances between $\mathbf{x}^{\star}$ and the constraint-
associated hyperplanes, the possibility of _capturing_ the optimal point is
increased. For a specific sample size $m^{\prime}=mn$, when the volume of the
finite space around the optimal point is reduced, the average of
$\left\\{d_{j}\left(\tilde{\mathbf{x}}^{\star},\bm{\uptau}^{(\ell)}\right)\right\\}$
is diminished as well. This average of distance can be written as [45]:

$\displaystyle
T_{\text{ave}}=\frac{1}{m^{\prime}}\sum^{m^{\prime}}_{j=1}d_{j}\left(\tilde{\mathbf{x}}^{\star},\bm{\uptau}^{(\ell)}\right).$
(46)

The possibility of creating a finite-volume, and the importance of the number
of samples in the recovery performance of RKA, can be captured by the
_Chernoff bound_ as illustrated below.

######  Theorem 3 (See[16]).

Assume the distances
$\left\\{d_{j}\left(\tilde{\mathbf{x}}^{\star},\bm{\uptau}^{(\ell)}\right)\right\\}$
between the desired point $\tilde{\mathbf{x}}^{\star}$ and the hyperplanes of
the polyhedron defined in (12) are i.i.d. random variables. Then:

  * •

The Chernoff bound of $T_{\text{ave}}$ is given by

$\operatorname{Pr}\left(\frac{1}{m^{\prime}}\sum^{m^{\prime}}_{j=1}d_{j}\left(\tilde{\mathbf{x}}^{\star},\bm{\uptau}^{(\ell)}\right)\leq
a\right)\geq 1-\inf_{t\geq 0}\frac{\Psi_{T}}{e^{ta}},$ (47)

where $a$ is an average distance point in space at which the finite-volume
space around the desired signal is created, and $\Psi_{T}$ is the moment
generating function (MGF) of the error recovery, given as

$\Psi_{T}=\left(1+t\frac{\mu^{(1)}_{d_{j}}}{m^{\prime}}+\cdots+t^{\kappa}\frac{\mu^{(\kappa)}_{d_{j}}}{\kappa!m^{\prime\kappa}}+\mathcal{R}\left(m^{\prime}\right)\right)^{m^{\prime}},$
(48)

with $\mu^{(\kappa)}_{d_{j}}=\mathbb{E}\left\\{d^{\kappa}_{j}\right\\}$, and
$\mathcal{R}$ denoting a bounded remainder associated with truncating the
Taylor series expansion of $\Psi_{T}$.

  * •

$\Psi_{T}$ is decreasing with an increasing sample size in the sample
abundance scenario, leading to an increasing lower bound in (47).

The MGF $\Psi_{T}$ contains two parts. The first part has an increasing trend
until a specific sample size $m^{\star}$, which indicates the existence of an
abundant number of samples. After that, the function has a decreasing
behavior. Therefore, $\Psi_{T}$ with $m\geq m^{\star}$ can be a good choice
for a sample size-dependent penalty function. Particularly, the penalty
function can be chosen as $\Psi_{T}-\Psi_{\infty}$, where
$\Psi_{\infty}=\lim_{m\rightarrow\infty}\Psi_{T}$, to ensure
$\Upsilon(m)\rightarrow 0$ as $m\rightarrow\infty$.

Since, we do not have access to the probability density function of
$\left\\{d_{j}\right\\}$, thus, the MGF must be evaluated by the truncated
Taylor series expansion, which may be accurately approximated by a rational
function such as a _Padé approximation_ (PA). The decreasing part of
$\Psi_{T}$ in $m>m^{\star}$ with PA is modeled as follows[32, 38]:

$\displaystyle\Upsilon(m)$
$\displaystyle\asymp\left(1+\cdots+t^{\kappa}\frac{\mu^{(\kappa)}_{d_{j}}}{\kappa!m^{\prime\kappa}}\right)^{m^{\prime}}-\Psi_{T},$
(49)
$\displaystyle=\frac{a_{0}+\frac{a_{1}}{m}}{b_{0}+\frac{b_{1}}{m}}-\frac{a_{0}}{b_{0}},$

where $\left\\{a_{0},a_{1},b_{0},b_{1}\right\\}$ are the PA coefficients as
given by

$\displaystyle a_{0}$ $\displaystyle=e^{u}\left(12u^{2}-24v\right),$ (50)
$\displaystyle a_{1}$
$\displaystyle=e^{u}\left(-3u^{4}+8u^{3}+12u^{2}v-24uv-12v^{2}\right),$
$\displaystyle b_{0}$ $\displaystyle=12u^{2}-24v,$ $\displaystyle b_{1}$
$\displaystyle=3u^{4}+8u^{3}-12u^{2}v-24uv+12v^{2},$

where $u=\mu_{d_{j}}^{(1)}t$ and $v=\frac{\mu_{d_{j}}^{(2)}t^{2}}{2}$.

##  V Judicious Sampling With Adaptive Thresholding for ORKA

By the spirit of using the iterative RKA, a suitable time-varying threshold
can be selected in order to enhance the recovery performance. In ORKA, we face
a highly overdetermined linear feasibility problem creating a finite-volume
space. To capture the desired signal $\mathbf{x}^{\star}$ more efficiently,
the right-hand side of the inequalities in (37), i.e.
$\operatorname{vec}\left(\mathbf{R}\right)\odot\operatorname{vec}\left(\bm{\Gamma}\right)$,
must be determined in a way that each associated hyperplane passes through the
desired feasible region within $\mathcal{F}_{\mathbf{X}}$. Therefore, an
algorithm is proposed to ensure that this occurs in practice.

We propose an iterative algorithm generating an adaptive sampling threshold to
accurately obtain the desired solution. To have the smaller area of the
finite-volume space around the desired signal $\mathbf{x}^{\star}$, one can
somehow choose thresholds to reduce distances between them and the desired
point. To do so, we update the time-varying threshold for
$\ell\in\left\\{1,\cdots,m\right\\}$ as

$\mathbf{A}\mathbf{x}_{k}-\mathbf{r}^{(\ell)}_{k}\odot\bm{\epsilon}^{(\ell)}_{k}=\bm{\uptau}^{(\ell)}_{k+1},$
(51)

where $\bm{\epsilon}^{(\ell)}_{k}$ are positive vectors in the $k$-th
iteration of the algorithm. This updating process is based on

$r_{j}=\begin{cases}+1&\mathbf{a}_{j}\mathbf{x}>\uptau_{j},\\\
-1&\mathbf{a}_{j}\mathbf{x}<\uptau_{j},\end{cases}$ (52)

where $\\{\mathbf{a}_{j}\\}$ are the rows of $\mathbf{A}$. The one-bit
measurements $\left\\{\mathbf{r}^{(\ell)}_{k}\right\\}$ are updated in the way
to satisfy (37) in iteration $k$, i.e. the inequalities
$\Omega^{(\ell)}_{k}\mathbf{A}\mathbf{x}^{\star}\geq\mathbf{r}^{(\ell)}_{k}\odot\bm{\uptau}^{(\ell)}_{k}$.
The reason behind this updating is to ensure that each halfspace associated
with a threshold in iteration $k$ is getting closer to the optimal point in
the correct direction which means the main side of the halfspace is forced to
cover the optimal solution.

######  Proposition 4 (ORKA with Adaptive Thresholding).

Consider applying ORKA to a linear feasibility problem
$\mathbf{A}\mathbf{x}=\mathbf{y}$ as part of the linear constraints of (1).
Suppose the initial time-varying threshold sequences are
$\left\\{\bm{\uptau}^{(\ell)}_{0}\right\\}\sim\mathcal{N}\left(0,1\right)$
(with the same length as $\mathbf{r}^{(\ell)}$), and
$\left\\{\delta^{(\ell)}\right\\}$ are positive numbers. Also,
$\mathbf{x}_{k}$, $\bm{\uptau}^{(\ell)}_{k}$, $\mathbf{r}^{(\ell)}_{k}$ and
$\bm{\epsilon}_{k}$ denote their associated values at iteration $k$. The
proposed sampling algorithm is summarized as follows:

  1. 1.

Find a point inside the following polyhedron with proposed accelerated
Kaczmarz algorithms, i.e. the PrSKM or the block SKM for
$\bm{\uptau}=\left\\{\bm{\uptau}^{(\ell)}_{k}\right\\}$:

$\mathcal{P}_{k}=\left\\{\mathbf{x}_{k}\mid\tilde{\bm{\Omega}}_{k}\mathbf{A}\mathbf{x}_{k}\succeq\mathbf{b}_{k}\right\\},$
(53)

where
$\mathbf{b}_{k}=\operatorname{vec}\left(\mathbf{R}_{k}\right)\odot\operatorname{vec}\left(\bm{\Gamma}_{k}\right)$,
$\mathbf{R}_{k}$ and $\bm{\Gamma}_{k}$ are matrices with
$\left\\{\mathbf{r}^{(\ell)}_{k}\right\\}$ and
$\left\\{\bm{\uptau}^{(\ell)}_{k}\right\\}$ representing their columns,
respectively.

  2. 2.

Update $\bm{\Gamma}_{k+1}$ as:

$\operatorname{vec}\left(\mathbf{R}_{k}\right)\odot\operatorname{vec}\left(\bm{\Gamma}_{k+1}\right)=\tilde{\bm{\Omega}}_{k}\mathbf{A}\mathbf{x}_{k}-\frac{\bm{\epsilon}_{k}}{2}.$
(54)

  3. 3.

Compute $\bm{\epsilon}_{k}$, a block vector containing
$\left\\{\bm{\epsilon}^{(\ell)}_{k}\right\\}$, as:

$\bm{\epsilon}_{k}=\tilde{\bm{\Omega}}_{k}\mathbf{A}\mathbf{x}_{k}-\operatorname{vec}\left(\mathbf{R}_{k}\right)\odot\operatorname{vec}\left(\bm{\Gamma}_{k}\right).$
(55)

  4. 4.

Update $\mathbf{R}_{k+1}$ based on (52):

$\mathbf{r}^{(\ell)}_{k+1}=\operatorname{sgn}\left(\mathbf{y}-\bm{\uptau}^{(\ell)}_{k+1}\right),\quad\ell\in\left\\{1,\cdots,m\right\\}.$
(56)

  5. 5.

Increase $k$ by one.

  6. 6.

Stop when
$\left\|\bm{\uptau}^{(\ell)}_{k+1}-\bm{\uptau}^{(\ell)}_{k}\right\|_{2}\leq\delta^{(\ell)}$.

One can observe that by deploying this adaptive thresholding algorithm,
smaller values of $\\{d_{j}\\}$ will emerge which leads to their moments
$\left\\{\mu^{(\kappa)}_{d_{j}}\right\\}$ to further diminish. Therefore,
$\Psi_{T}$ is smaller in this scenario and a smaller number of time-varying
sampling threshold sequences can be utilized in ORKA with similar recovery
performance. Additionally, non-informative sampling thresholds, which appear
as extra inequality constraints in the random time-varying sampling thresholds
scenario, may be efficiently removed by choosing the adaptive thresholds with
closer hyperplanes to the desired point.

(a) $m=2$

(b) $m=6$

(c) $m=60$

(d) $m=2$

(e) $m=6$

(f) $m=60$

Figure 2: Shrinkage of the one-bit polyhedron (37) in blue, ultimately placed
within the unit ball of the nuclear norm
$\left\|\mathbf{X}\right\|_{\star}\leq 1$ shown with black cylindrical region
and its red contours, when the number constraints (samples) grows large. The
arrows point to the half-space associated with each inequality constraint. The
evolution of the feasible regime is depicted with increasing samples in three
cases: (a) and (d) small sample-size regime, constraints not forming a finite-
value polyhedron; (b) and (e) medium sample-size regime, constraints forming a
finite-volume polyhedron, parts of which are outside the cylinder; (c) and (f)
large sample-size regime, constraints forming a finite-volume polyhedron
inside the nuclear norm cylinder, making its constraint redundant. The optimal
point representing the signal to be recovered is shown by yellow.

##  VI Low-Rank Matrix Recovery Via ORKA

As mentioned earlier, low-rank matrix recovery is an excellent example for
problems that assume the form in (1), and that can be tackled using our
methodology. In this section, at first, we will briefly introduce the _nuclear
norm minimization_ form of the problem. Accordingly, we will apply ORKA to
this problem without considering the associated costly constraints. At the
end, the recovery performance of ORKA will be numerically evaluated
considering different matrix ranks and sample sizes to investigate the
existence of a sample abundance scenario.

###  VI-A Problem Formulation

The problem of the low-rank matrix recovery can be formulated as:

find $\displaystyle\mathbf{X}$ (57) s.t.
$\displaystyle\mathcal{A}\left(\mathbf{X}\right)=\mathbf{y},$
$\displaystyle\operatorname{rank}\left(\mathbf{X}\right)\leq M,$
$\displaystyle\mathbf{X}\in\Omega_{c},$

where $\mathbf{X}\in\mathbb{C}^{n_{1}\times n_{2}}$ is the matrix of unknowns,
$\mathbf{y}\in\mathbb{R}^{n}$ is the measurement vector, and $\mathcal{A}$ is
a linear transformation mapping $n_{1}\times n_{2}$ into $\mathbb{R}^{n}$. In
general, $\Omega_{c}$ can be chosen such as the set of semi-definite matrices,
symmetric matrices, upper or lower triangle matrices, Hessenberg matrices and
a specific constraint on the matrix elements
$\left\|\mathbf{X}\right\|_{\infty}\leq\alpha$ or on its eigenvalues, i.e.,
$\lambda_{i}\leq\epsilon$ where $\left\\{\lambda_{i}\right\\}$ are eigenvalues
of $\mathbf{X}$ [4, 1, 53].

The problem (57) can be rewritten as an optimization problem:

$\displaystyle\min_{\mathbf{X}}$
$\displaystyle\operatorname{rank}\left(\mathbf{X}\right)$ (58) s.t.
$\displaystyle\mathcal{A}\left(\mathbf{X}\right)=\mathbf{y},$
$\displaystyle\mathbf{X}\in\Omega_{c}.$

This problem is known to be NP-hard, whose solution is difficult to
approximate [56, 15]. Recall that the rank of $\mathbf{X}$ is equal to the
number of nonzero singular values. In the case when the singular values are
all equal to one, the sum of the singular values is equal to the rank. When
the singular values are less than or equal to one, the sum of the singular
values is a convex function that is strictly less than the rank. Therefore, it
is been popular for this problem to replace the rank function with the sum of
the singular values of $\mathbf{X}$; i.e., its nuclear norm. The nuclear norm
minimization alternative of the problem is given by [2, 57, 15]:

$\displaystyle\min_{\mathbf{X}}$
$\displaystyle\left\|\mathbf{X}\right\|_{\star}$ (59) s.t.
$\displaystyle\mathcal{A}\left(\mathbf{X}\right)=\mathbf{y},$
$\displaystyle\mathbf{X}\in\Omega_{c}.$

In this problem, the feasible set $\mathcal{F}_{\mathbf{X}}$ is obtained as

$\mathcal{F}_{\mathbf{X}}=\left\\{\mathcal{P}^{\star}\cap\Omega_{c}\right\\},$
(60)

where $\mathcal{P}^{\star}$ is defined as follows

$\mathcal{P}^{\star}=\left\\{\mathbf{X}\mid\left\|\mathbf{X}\right\|_{\star}\leq\tau\right\\},\quad\tau\in\mathbb{R}^{+}.$
(61)

Next, we will apply ORKA to (58) to make its costly constraints redundant by
using abundant number of time-varying sampling thresholds $m$.

A numerical investigation of (37) when it is achieved for the nuclear norm
minimization, reveals that by increasing the number of time-varying sampling
threshold sequences $m$, the space formed by the intersection of half-spaces
(inequality constraints) can fully shrink to the desired signal
$\mathbf{X}^{\star}$ inside the feasible region of (61) which is shown by the
cylindrical space [15]—see Fig. 2 for an illustrative example of this
phenomenon. As can be seen in this figure, the blue lines displaying the
linear feasibility form a finite-volume space around the optimal point
displayed by the yellow circle inside the cylinder (the elliptical region) by
growing the number of threshold sequences or one-bit samples. In (a)/(d),
constraints are not enough to create a finite-volume space, whereas in (b)/(e)
such constraints can create the desired finite-volume polyhedron space which,
however, is not fully inside the cylinder. Lastly, in (c)/(f), the created
finite-volume space shrinks to be fully inside the cylinder.

(a)

(b)

Figure 3: Average NMSE for the Frobenius norm of error for the recovery of the
matrix $\mathbf{X}$ associated with different time-varying sampling threshold
sequences sizes when the PrSKM and the block SKM are utilized in ORKA: (a)
$\operatorname{rank}\left(\mathbf{X}\right)=1$, (b)
$\operatorname{rank}\left(\mathbf{X}\right)=4$.

(a)

(b)

Figure 4: Comparing the average NMSE for the Frobenius norm of error for the
recovery of the matrix $\mathbf{X}$ using ORKA when (i) a random threshold and
(ii) the adaptive sampling threshold are adopted when the PrSKM and the block
SKM are utilized in ORKA: (a) $\operatorname{rank}\left(\mathbf{X}\right)=1$,
(b) $\operatorname{rank}\left(\mathbf{X}\right)=4$.

###  VI-B Numerical Illustrations

In this section, we numerically scrutinize the capability of the ORKA in the
nuclear norm minimization problem (59) instead of (59) by the squared
Frobenius norm of the error normalized by the squared Frobenius norm of the
desired matrix $\mathbf{X}^{\star}$, defined as

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

We solve the overdetermined one-bit polyhedron in (37) via the PrSKM and the
Block SKM. To make this happen, we obtain the one-bit polyhedron from a linear
feasibility problem $\mathbf{A}\mathbf{x}=\mathbf{y}$, where
$\mathbf{A}\in\mathbb{R}^{200\times 25}$, $\mathbf{x}\in\mathbb{R}^{25}$
($\mathbf{x}=\operatorname{vec}\left(\mathbf{X}\right)$ where
$\mathbf{X}\in\mathbb{R}^{5\times 5}$), and $\mathbf{y}\in\mathbb{R}^{200}$.
We consider the number of time-varying sampling threshold sequences to be
$m\in\left\\{10,20,30,40,50,60\right\\}$. Each row of $\mathbf{A}$ is
generated as
$\mathbf{a}_{j}\sim\mathcal{N}\left(\mathbf{0},\mathbf{I}_{25}\right)$. For
the desired matrix $\mathbf{X}$, we generate
$\mathbf{X}=\mathbf{K}\mathbf{K}^{\top}$ where (i)
$\mathbf{K}\in\mathbb{R}^{5\times 4}$ is the Gaussian matrix, and (ii)
$\mathbf{K}\in\mathbb{R}^{5\times 1}$ is the Gaussian vector. Also, each time-
varying sampling threshold $\bm{\uptau}^{(\ell)}$ is considered to have the
distribution
$\bm{\uptau}^{(\ell)}\sim\mathcal{N}\left(\mathbf{0},\mathbf{I}_{200}\right)$.
Fig. 3 appears to confirm the possibility of recovering the desired matrix
$\mathbf{X}^{\star}$ in the one-bit polyhedron (37) by ORKA. As expected, the
performance of the recovery will be significantly enhanced as the number of
time-varying sampling threshold sequences grows large. Also, similar to
before, it can be seen that the Block SKM outperforms the PrSKM in the low-
rank matrix recovery problem.

To improve the recovery performance, we proposed the adaptive time-varying
sampling threshold in Section V. Fig. 4 illustrates the performance of the
Block SKM in the low rank matrix recovery in the one-bit polyhedron (37) when
we have the high-dimensional input signal $\mathbf{x}\in\mathbb{R}^{128}$ and
$\mathbf{A}\in\mathbb{R}^{20000\times 128}$, with (i) a random threshold, and
(ii) an adaptive time-varying threshold. As can be seen, the recovery
performance is significantly enhanced when the Block SKM is equipped with the
adaptive time-varying threshold.

(a)

(b)

Figure 5: Average NMSE for the error between the desired signal
$\mathbf{x}^{\star}$ and its recovered version $\bar{\mathbf{x}}$ for
different time-varying sampling threshold sequences sizes when the PrSKM and
the block SKM are utilized in ORKA with (a) $k=2$, (b) $k=4$.

##  VII One-Bit Compressed Sensing:
From Optimization to Linear Feasibility

Compressed sensing (CS) is an interesting and rapidly growing area of research
that has attracted considerable attention in electrical engineering, applied
mathematics, statistics, and computer science [29, 4]. In CS, a sparse high-
dimensional signal is to be recovered by incomplete measurements such a
recovery may be formulated as [29]:

$\displaystyle\min_{\mathbf{x}}$ $\displaystyle\left\|\mathbf{x}\right\|_{1}$
(63) s.t. $\displaystyle\mathbf{A}\mathbf{x}=\mathbf{y},$

where $\mathbf{A}\in\mathbb{R}^{m\times n}$, and $m\ll n$. One of the
important applications of CS emerges in the signal recovery from a sequence of
acquisitions $\\{y_{i}\\}$ obtained from a sparse linear transformation
(wavelet transformations are known for such a property, for instance) in the
magnetic resonance imaging (MRI). The reconstruction problem of the desired
signal $\mathbf{x}^{\star}$ is given by

$\displaystyle\min_{\mathbf{x}}$ $\displaystyle\left\|\mathbf{x}\right\|_{1}$
(64) s.t. $\displaystyle\mathcal{A}_{i}\left(\mathbf{x}\right)=y_{i},\quad
i\in\left\\{1,\cdots,n\right\\}.$

In this section, we first formulate the optimization problem of the _one-bit
compressed sensing_. Then, by taking advantage of one-bit sampling, we
increase sample size in (64) and create an associated one-bit polyhedron.

###  VII-A Problem Formulation

Let $\bm{\uptau}$ denotes the time-varying threshold vector. The one-bit
samples are generated as

$\displaystyle r_{i}$
$\displaystyle=\begin{cases}+1&\mathbf{a}^{\top}_{i}\mathbf{x}\geq\uptau_{i},\\\
-1&\mathbf{a}^{\top}_{i}\mathbf{x}<\uptau_{i},\end{cases}$ (65)

where
$\mathcal{A}_{i}\left(\mathbf{x}\right)=\mathbf{a}^{\top}_{i}\mathbf{x}$. The
occurrence probability vector $\mathbf{p}$ for the one-bit measurement
$\mathbf{r}$ is given as [16],

$\displaystyle p_{i}$
$\displaystyle=\begin{cases}\Phi\left(\mathbf{a}^{\top}_{i}\mathbf{x}\right)&\text{for}\quad\\{r_{i}=+1\\},\\\
1-\Phi\left(\mathbf{a}^{\top}_{i}\mathbf{x}\right)&\text{for}\quad\\{r_{i}=-1\\},\end{cases}$
(66)

where $\Phi(.)$ is the CDF of $\bm{\uptau}$. The log-likelihood function of
the sign data $\mathbf{r}$ is given by

$\displaystyle\mathcal{L}_{\mathbf{r}}(\bm{\mu},\mathbf{x})$
$\displaystyle=\sum^{m}_{i=1}\left\\{\mathbb{I}_{(r_{i}=+1)}\log\left(\Phi(\mathbf{a}^{\top}_{i}\mathbf{x})\right)\right.$
(67)
$\displaystyle\left.+\mathbb{I}_{(r_{i}=-1)}\log\left(1-\Phi(\mathbf{a}^{\top}_{i}\mathbf{x})\right)\right\\}.$

Therefore, the maximum likelihood estimation (MLE) for the one-bit compressed
sensing can be written as

$\min_{\mathbf{x}}\mathcal{L}_{\mathbf{r}}(\bm{\mu},\mathbf{x})+\lambda\|\mathbf{x}\|_{1}.$
(68)

The alternative formulations for one-bit compressed sensing can be found in
[58].

Nevertheless, as discussed earlier, by deploying one-bit sampling, the
opportunity exists to increase the number of samples in (64). The one-bit
compressed sensing is thus solely accomplished by creating a highly-
constrained one-bit polyhedron. In other words, instead of solving an
optimization problem with costly constraints, the problem may be tackled by
the proposed accelerated Kaczmarz algorithms; namely, PrSKM and the block SKM.

(a)

(b)

Figure 6: Comparing the average NMSE for the desired signal
$\mathbf{x}^{\star}$ and its recovered signal using ORKA when (i) a random
threshold and (ii) the adaptive sampling threshold are adopted with (a)
$k=20$, (b) $k=40$

###  VII-B Numerical results

To examine the performance of ORKA in CS and to validate the theoretical
results described in this paper, we consider signal recovery with different
number of time-varying sampling threshold sequences
$m\in\left\\{10,20,30,40,50,60\right\\}$. Input signals
$\mathbf{x}^{\star}\in\mathbb{R}^{10}$ are generated with sparsity orders
$k=2$ and $k=4$, respectively. The sparsity order $k$ is defined as the number
of non-zero elements in a vector. Time-varying sampling thresholds and the
constraint matrix $\mathbf{A}$ are generated as in Subsection VI-B. To compare
two proposed algorithms, the NMSE defined in (39) is utilized and the results
are averaged over $15$ experiments.

As can be seen in Fig. 5, by increasing the number of time-varying sampling
threshold sequences, the performance of our method is improved. Beside the
possibility of increasing the number of measurements $n$, the higher number of
samples are available in ORKA by comparing the measurements with multiple
threshold sequences $\ell\in\left\\{1,\cdots,m\right\\}$. In other words, we
have the opportunity to increase $n$ and $m$ simultaneously, when the number
of samples is $m^{\prime}=mn$.

Same as Subsection VI-B, the adaptive thresholding algorithm is applied to
ORKA for the high-dimensional input signal $\mathbf{x}\in\mathbb{R}^{128}$ in
order to enhance its recovery performance, whose outcome is presented in Fig.
6. The NMSE results are reported with sparsity orders $k=20$ and $k=40$.

Table I: Comparing CPU times and $\operatorname{NMSE}$ of ORKA and $\ell^{1}$-minimization. Algorithm | $m^{\star}$ | CPU time (s) | $\operatorname{NMSE}$
---|---|---|---
ORKA | $500$ | $3.1240e-04$ | $3.2052e-12$
$\ell^{1}$ | $100$ | $0.0071$ | $2.4572e-11$

To further investigate the efficacy of ORKA in CS, we compare our proposed
approach with the well-known $\ell^{1}$-minimization approach formulated in
(64) in terms of NMSE and CPU time. As presented in Table I, ORKA outperforms
$\ell^{1}$-minimization in terms of both NMSE and CPU time. The results are
obtained for $\mathbf{x}\in\mathbb{R}^{128}$ when the optimal number of
samples are utilized, and where $m^{\star}=4k\log(n/k)$ and $m^{\star}=500$
are considered for the high-resolution method and ORKA, respectively. Herein,
optimality of sample sizes means that the number of samples utilized by
algorithms leads to their best performance, i.e. satisfying the criterion
$\left\|\mathbf{x}_{i}-\mathbf{x}^{\star}\right\|_{2}^{2}\leq 5\times
10^{-11}\left\|\mathbf{x}^{\star}\right\|_{2}^{2}$. By this comparison, we
remove the burden of the large number of samples from the
$\ell^{1}$-minimization to fairly compare their optimal shape deploying
incomplete measurements with that of ORKA.

It is worth pointing out that for a $64$-bit ADC, $m=100$ corresponds to
$6400$ bits of information while ORKA solely employs $500$ bits. Therefore, it
appears from Table I that ORKA achieves a better accuracy in terms of NMSE
with not only fewer information bits but also a smaller computational cost.

##  VIII Conclusion

We proposed a novel algorithm, ORKA, that takes advantage of the abundant
number of samples available in one-bit sampling with time-varying thresholds
to efficiently and globally solve some well-studied problems in the form of
(1); including low-rank matrix recovery and compressed sensing. Moreover, two
state-of-the-art randomized Kaczmarz algorithms are proposed to use in ORKA 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 low-
rank matrix recovery and compressed sensing problems.

## References

  * [1] 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.
  * [2] J. Cai, E.J. Candès, and Z. Shen,  “A singular value thresholding algorithm for matrix completion,”  SIAM Journal on optimization, vol. 20, no. 4, pp. 1956–1982, 2010\.
  * [3] M. A. Davenport, Y. Plan, E. Van Den Berg, and M. Wootters,  “1-bit matrix completion,”  Information and Inference: A Journal of the IMA, vol. 3, no. 3, pp. 189–223, 2014.
  * [4] M.A. Davenport and J. Romberg,  “An overview of low-rank matrix recovery from incomplete observations,”  IEEE Journal of Selected Topics in Signal Processing, vol. 10, no. 4, pp. 608–622, 2016.
  * [5] E.J. Candes and Y. Plan,  “Tight oracle inequalities for low-rank matrix recovery from a minimal number of noisy random measurements,”  IEEE Transactions on Information Theory, vol. 57, no. 4, pp. 2342–2359, 2011.
  * [6] Y. Chen, A. Jalali, S. Sanghavi, and C. Caramanis,  “Low-rank matrix recovery from errors and erasures,”  IEEE Transactions on Information Theory, vol. 59, no. 7, pp. 4324–4337, 2013.
  * [7] X. Li, L. Huang, H. So, and B. Zhao,  “A survey on matrix completion: Perspective of signal processing,”  arXiv preprint arXiv:1901.10885, 2019.
  * [8] B. Haeffele, E. Young, and R. Vidal,  “Structured low-rank matrix factorization: Optimality, algorithm, and applications to image processing,”  in International conference on machine learning. PMLR, 2014, pp. 2007–2015.
  * [9] F. Nie, H. Huang, and C. Ding,  “Low-rank matrix recovery via efficient schatten p-norm minimization,”  in Twenty-sixth AAAI conference on artificial intelligence, 2012\.
  * [10] David M Blei,  “Probabilistic topic models,”  Communications of the ACM, vol. 55, no. 4, pp. 77–84, 2012.
  * [11] G. Obozinski, B. Taskar, and M. Jordan,  “Joint covariate selection and joint subspace selection for multiple classification problems,”  Statistics and Computing, vol. 20, no. 2, pp. 231–252, 2010.
  * [12] M. Pontile, T. Evgeniou, and A. Argyriou,  “Convex multi-task feature learning,”  Journal of Machine Learning, vol. 10, pp. 243–272, 2007.
  * [13] A. Argyriou, T. Evgeniou, and M. Pontil,  “Multi-task feature learning,”  Advances in neural information processing systems, vol. 19, 2006\.
  * [14] C. Tomasi and T. Kanade,  “Shape and motion from image streams under orthography: a factorization method,”  International journal of computer vision, vol. 9, no. 2, pp. 137–154, 1992.
  * [15] B. Recht, W. Xu, and B. Hassibi,  “Null space conditions and thresholds for rank minimization,”  Mathematical programming, vol. 127, no. 1, pp. 175–202, 2011.
  * [16] 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.
  * [17] R. P. Millane,  “Phase retrieval in crystallography and optics,”  J. Opt. Soc. Am. A, vol. 7, no. 3, pp. 394–411, Mar 1990.
  * [18] W. Kim and M. H. Hayes,  “The phase retrieval problem in x-ray crystallography,”  in Proceedings of the Acoustics, Speech, and Signal Processing, 1991\. ICASSP-91., 1991 International Conference, USA, 1991, ICASSP ’91, p. 1765–1768, IEEE Computer Society.
  * [19] J. R. Fienup, J. C. Marron, T. J. Schulz, and J. H. Seldin,  “Hubble space telescope characterized by using phase-retrieval algorithms,”  Appl. Opt., vol. 32, no. 10, pp. 1747–1767, Apr 1993.
  * [20] J. E. Krist and C. J. Burrows,  “Phase-retrieval analysis of pre- and post-repair hubble space telescope images,”  Appl. Opt., vol. 34, no. 22, pp. 4951–4964, Aug 1995.
  * [21] A. M. Sarnik,  “Phase retrieval: A practical application for the space telescope,”  in Inverse Optics II, R.H.T. Bates and A. J. Devaney, Eds. International Society for Optics and Photonics, 1985, vol. 0558, pp. 85 – 94, SPIE.
  * [22] R. W. Gerchberg and W. O. Saxton,  “Phase determination for image and diffraction plane pictures in the electron microscope,”  Optik, vol. 34, pp. 275–284, 1971.
  * [23] R. W. Gerchberg and W. O. Saxton,  “A practical algorithm for the determination of phase from image and diffraction plane pictures,”  Optik, vol. 35, pp. 227–246, 1972.
  * [24] J. R. Fienup,  “Comments on ”the reconstruction of a multidimensional sequence from the phase or magnitude of its Fourier transform”,”  IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 31, no. 3, pp. 738–739, Jun 1983.
  * [25] J. R. Fienup,  “Phase retrieval algorithms: a comparison,”  Appl. Opt., vol. 21, no. 15, pp. 2758–2769, Aug 1982.
  * [26] J. R. Fienup,  “Reconstruction of an object from the modulus of its Fourier transform,”  Opt. Lett., vol. 3, no. 1, pp. 27–29, Jul 1978.
  * [27] 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.
  * [28] E.J. 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.
  * [29] Y.C. Eldar and G. Kutyniok,  Compressed sensing: theory and applications,  Cambridge university press, 2012.
  * [30] M. Lustig, D. Donoho, J.M. Santos, and J.M. Pauly,  “Compressed sensing MRI,”  IEEE signal processing magazine, vol. 25, no. 2, pp. 72–82, 2008\.
  * [31] J.A. Fessler,  “Model-based image reconstruction for MRI,”  IEEE signal processing magazine, vol. 27, no. 4, pp. 81–89, 2010\.
  * [32] 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.
  * [33] Texas Instruments,  “Analog embedded processing,”  Active Filter Design Techniques, Literature Number SLOA088 Texas Instruments, excerpted from OP-Amps for Everyone, Literature Number SLOD006A.
  * [34] 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.
  * [35] S. Sedighi, B. Shankar, M. Soltanalian, and B. Ottersten,  “One-bit DoA estimation via sparse linear arrays,”  in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2020, pp. 9135–9139.
  * [36] C. Qian and J. Li,  “ADMM for harmonic retrieval from one-bit sampling with time-varying thresholds,”  in 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2017, pp. 3699–3703.
  * [37] C. Gianelli, Luzhou Xu, Jian Li, and Petre Stoica,  “One-bit compressive sampling with time-varying thresholds for sparse parameter estimation,”  in 2016 IEEE Sensor Array and Multichannel Signal Processing Workshop (SAM). IEEE, 2016, pp. 1–5.
  * [38] Arian Eamaz, Farhang Yeganegi, and Mojtaba 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.
  * [39] P. Wang, J. Li, M. Pajovic, P. T. Boufounos, and Philip V Orlik,  “On angular-domain channel estimation for one-bit massive MIMO systems with fixed and time-varying thresholds,”  in 51st Asilomar Conference on Signals, Systems, and Computers. IEEE, 2017, pp. 1056–1060.
  * [40] F. Xi, Y. Xiang, S. Chen, and A. Nehorai,  “Gridless parameter estimation for one-bit MIMO radar with time-varying thresholds,”  IEEE Transactions on Signal Processing, vol. 68, pp. 1048–1063, 2020.
  * [41] 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\.
  * [42] H.G. Feichtinger and T. Strohmer,  “A Kaczmarz-based approach to nonperiodic sampling on unions of rectangular lattices,”  in SampTA’95: 1995 Workshop on Sampling Theory and Applications, 1995, pp. 32–37.
  * [43] M. Sezan and H. Stark,  “Applications of convex projection theory to image recovery in tomography and related areas,”  Image Recovery: Theory and Application, pp. 155–270, 1987.
  * [44] 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.
  * [45] 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.
  * [46] 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\.
  * [47] J. Briskman and D. Needell,  “Block Kaczmarz method with inequalities,”  Journal of Mathematical Imaging and Vision, vol. 52, no. 3, pp. 385–396, 2015.
  * [48] 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.
  * [49] I. Daubechies, M. Defrise, and C. De Mol,  “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,”  Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, vol. 57, no. 11, pp. 1413–1457, 2004.
  * [50] P.L. Combettes and V. Wajs,  “Signal recovery by proximal forward-backward splitting,”  Multiscale modeling & simulation, vol. 4, no. 4, pp. 1168–1200, 2005.
  * [51] L. Dai, M. Soltanalian, and K. Pelckmans,  “On the randomized Kaczmarz algorithm,”  IEEE Signal Processing Letters, vol. 21, no. 3, pp. 330–333, 2013\.
  * [52] M. Sarowar Morshed and M. Saiful Islam,  “Sampling Kaczmarz Motzkin method for linear feasibility problems: Generalization & acceleration,”  arXiv e-prints, pp. arXiv–2002, 2020.
  * [53] C. F. Van Loan and G. Golub,  “Matrix computations (Johns Hopkins studies in mathematical sciences),”  1996\.
  * [54] A. Bhandari, F. Krahmer, and R. Raskar,  “On unlimited sampling and reconstruction,”  IEEE Transactions on Signal Processing, vol. 69, pp. 3827–3839, 2020.
  * [55] T. Elfving,  “Block-iterative methods for consistent and inconsistent linear equations,”  Numerische Mathematik, vol. 35, no. 1, pp. 1–12, 1980.
  * [56] R. Meka, P. Jain, C. Caramanis, and I.S. Dhillon,  “Rank minimization via online learning,”  in Proceedings of the 25th International Conference on Machine learning, 2008, pp. 656–663.
  * [57] B. Recht, M. Fazel, and P. A. Parrilo,  “Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization,”  SIAM review, vol. 52, no. 3, pp. 471–501, 2010.
  * [58] S. Khobahi and M. Soltanalian,  “Model-based deep learning for one-bit compressive sensing,”  IEEE Transactions on Signal Processing, vol. 68, pp. 5292–5307, 2020.