# Parallelization Strategies for the Randomized Kaczmarz Algorithm on Large-
Scale Dense Systems

Inês Ferreira<EMAIL_ADDRESS>Juan A. Acebrón
<EMAIL_ADDRESS>José Monteiro<EMAIL_ADDRESS>
###### Abstract

The Kaczmarz algorithm is an iterative technique designed to solve consistent
linear systems of equations. It falls within the category of row-action
methods, focusing on handling one equation per iteration. This characteristic
makes it especially useful in solving very large systems. The recent
introduction of a randomized version, the Randomized Kaczmarz method, renewed
interest in the algorithm, leading to the development of numerous variations.
Subsequently, parallel implementations for both the original and Randomized
Kaczmarz method have since then been proposed. However, previous work has
addressed sparse linear systems, whereas we focus on solving dense systems. In
this paper, we explore in detail approaches to parallelizing the Kaczmarz
method for both shared and distributed memory for large dense systems. In
particular, we implemented the Randomized Kaczmarz with Averaging (RKA) method
that, for inconsistent systems, unlike the standard Randomized Kaczmarz
algorithm, reduces the final error of the solution. While efficient
parallelization of this algorithm is not achievable, we introduce a block
version of the averaging method that can outperform the RKA method.

###### keywords:

Dense matrices , Linear systems , Iterative algorithms , Kaczmarz algorithm ,
Parallel and distributed computing , Least-Squares problem , Randomized
Algorithms

###### MSC:

15A06, 15A52, 65F10, 65F20, 68W20, 65Y05, 68W10, 68W15

††journal: Journal of Parallel and Distributed Computing

[inst1]organization=INESC-ID,addressline=Instituto Superior Técnico,
city=Universidade de Lisboa, country=Portugal

[inst2]organization=Department of Mathematics,addressline=Carlos III
University of Madrid, country=Spain

##  1 Introduction

Solving linear systems of equations is a common problem in many areas of
engineering such as image reconstruction and signal processing. Given a matrix
of real elements, $A\in\mathbb{R}^{m\times n}$, and a vector of constants,
$b\in\mathbb{R}^{m}$, solving a linear system of equations can be reduced to
finding the solution, $x\in\mathbb{R}^{n}$, that satisfies

$Ax=b\>.$ (1)

In this paper, we consider overdetermined systems ($m\geq n$) since many real-
world problems fall into this category. More specifically, in problems reliant
on data measurements, it is common for the matrix’s row count to scale
proportionally with the volume of collected data. Consequently, when data is
plentiful, the resulting matrix tends to be rectangular and the system is
overdetermined. Two examples are the following. First, in the context of
camera calibration, it is possible to characterize the relationship between a
set of 3D points in the real world and their projection onto the image plane
of a camera using a transformation. The parameters of this transformation can
be computed by solving a linear system that is overdetermined if there are
more than 4 data entries. For more information see Section 3.2 of [1]. A
second example is during a Computed Tomography (CT scan). The problem of
reconstructing an image of a scanned body in a CT scan can be reduced to
solving a linear system. The size of this linear system is dependent on
factors such as the image’s pixel count, the number of measurement angles
employed, and the detector array’s size [2]. Depending on the relationship
between these parameters, overdetermined systems are not an uncommon
occurrence. If a solution for (1) exists we say that the system is consistent.
For full-rank consistent systems, there exists a single solution, which we
denote by $x^{*}$. When there is no solution, the system is said to be
inconsistent, and we are usually interested in finding the least-squares
solution, that is

$x_{LS}=arg\min_{x}\|Ax-b\|^{2}\>,$ (2)

where $\|\>.\>\|$ is the Euclidean $L^{2}$ norm. The direct computation of the
least-squares solution can be performed using
$x_{LS}=(A^{T}A)^{-1}A^{T}b\>=A^{\dagger}b$, where $A^{\dagger}$ represents
the Moore–Penrose inverse or pseudoinverse. Real-world linear systems are
usually constructed using data obtained from measuring physical quantities. As
a result of inherent measurement errors, it is more common to have
inconsistent than consistent systems.

There are two classes of numerical methods to solve linear systems of
equations: direct and iterative. When matrices are large, iterative methods
are usually used to the detriment of direct methods, since they do not require
the manipulation of the input matrix. A particular class of iterative methods
are row-action methods [3] that use a single equation per iteration. These
have two advantages over other iterative methods: they can be used in real-
time while data is still being collected and they can be used to solve very
large systems that do not fit in a single machine. The Kaczmarz method [4] is
one of such methods.

Improving the performance of algorithms can be achieved through
parallelization. Parallelizing an algorithm can not only decrease the
execution time but, in the case of distributed memory parallelizations, it can
also enable the processing of large amounts of data, something that could not
be attained if we were to use one machine only. However, parallelization is
not a straightforward technique since we must consider data dependencies
between tasks, communication overhead between processing units, and
synchronization points, among others.

In this paper, we discuss several parallelization strategies for the Kaczmarz
algorithm targeting both shared and distributed memory systems. More
specifically, we implemented the Randomized Kaczmarz with Averaging [5], and
we argue that it is not possible to achieve an efficient parallelization of
this method. We introduce a new block version of this method and show that
using, the same relaxation parameters, the parallelization of this algorithm
is more efficient than the parallelization of the averaging method without
blocks.

The organization of this document is as follows. In Section 2, we present the
original and randomized versions of the Kaczmarz algorithm. Furthermore, we
introduce some previous parallelization strategies for the algorithm. In
Section 3, we present the implementation details together with experimental
results. Finally, in Section 4, we conclude this paper and discuss future
work.

##  2 Background

###  2.1 The Kaczmarz Method

The Kaczmarz method was proposed in [4] as an iterative algorithm to find the
solution of a consistent linear system of equations. Each iteration of the
algorithm can be computed using

$x^{(k+1)}=x^{(k)}+\alpha_{i}\>\frac{b_{i}\>-\langle
A^{(i)},x^{(k)}\rangle}{\|A^{(i)}\|^{2}}\>{A^{(i)}}^{T}\>,\quad\text{with}\quad
i=k\text{ mod }m\>,$ (3)

where $A^{(i)}$ is the $i$-th row of $A$, $b_{i}$ is the $i$-th coordinate of
$b$, $k$ is the iteration number (that starts at $0$), and
$\alpha_{i}\in(0,2)$ is a relaxation parameter. The initial guess $x^{(0)}$ is
typically set to zero, and $\langle,\rangle$ represents the dot product of two
vectors. Note that the rows of matrix $A$ are used cyclically, which makes
this algorithm also known as the Cyclic Kaczmarz (CK) method. The Kaczmarz
method has a simple geometric interpretation if $\alpha_{i}=1$: each iteration
$x^{(k+1)}$ is the projection of $x^{(k)}$ onto a hyperplane defined by a
different row of the system, that is, $H_{i}=\\{x:\langle
A^{(i)},x\rangle=b_{i}\\}$. In each iteration, the estimate of the solution
will satisfy a different constraint until a point is reached where all the
constraints are satisfied. In the case of underdetermined systems, the
Kaczmarz method converges to the least Euclidean norm solution [6]. For
inconsistent systems, the Kaczmarz method does not converge to the least-
squares solution since it remains always a certain distance from $x_{LS}$.

###  2.2 Randomized Kaczmarz Method

More recently, it has been observed that choosing the rows of $A$ randomly
instead of cyclically can accelerate the algorithm. One example of this
behavior is for systems with highly coherent matrices [7], that is, matrices
for which the angle between consecutive rows is small. In this case, the
convergence of the cyclical Kaczmarz is quite slow since the change in the
solution estimate in each iteration is minimal, as we can see in Figure 1a.
However, using the rows of the system randomly can lead to a faster
convergence, as Figure 1b shows.

(a) Cyclical selection of rows.

(b) Random selection of rows.

Figure 1: Convergence of the Kaczmarz method for a consistent system in 2
dimensions using two different row selection criteria.

Kaczmarz [4] proved the convergence of the method to $x^{*}$ for consistent
linear systems. However, it would be helpful to quantify the rate of
convergence of the algorithm in terms of the condition number, $k(A)$, so that
we can compare it to other iterative algorithms. This has proved to be a
difficult task since this algorithm depends on the order on which the
different rows of the matrix are processed.

Following these two statements, a randomized version of the algorithm was
introduced by Strohmer and Vershynin in [8] where rows are chosen with
probability proportional to their rows. In each iteration, we use the row with
index $i$, chosen at random from the probability distribution

$P\\{i=l\\}=\frac{\|A^{(l)}\|^{2}}{\|A\|_{F}^{2}}\quad(l=0,1,2,...,m-1)\>,$
(4)

where $\|A\|_{F}=\sqrt{\sum_{i=1}^{m}\|A^{(i)}\|^{2}}$ is the Frobenius norm
of a matrix. This new method called the Randomized Kaczmarz (RK) method,
converges linearly to the solution $x^{*}$. Similarly to the original Kaczmarz
method, the RK method does not converge to the least-squares solution of an
inconsistent system. However, Needell [9] proved that RK reaches a solution
estimate within a fixed distance from the solution, called the convergence
horizon. After Strohmer and Vershynin presented the RK method, many more
versions of the Kaczmarz method were developed, along with randomized versions
of other iterative methods. In [10], many of these variations are
experimentally compared by analyzing their performance in terms of the number
of iterations and execution time for overdetermined large dense systems.

###  2.3 Parallel Implementations of the Kaczmarz Method

It is not a trivial task to parallelize the Kaczmarz algorithm, as it is an
iterative algorithm where each iteration depends on the previous one. There
are, nonetheless, two main strategies for the parallelization of iterative
algorithms, both of which divide the equations into blocks: block-sequential
(also called block-iterative) and block-parallel. When using a block-
sequential approach, the blocks of equations are processed sequentially but
the work inside each block is parallelized. For a block-parallel approach,
blocks of equations are distributed among the processors and computed in
parallel, after which the results are combined. A block-sequential approach is
only efficient if there is a substantial amount of work in each iteration to
make up for the overhead of parallelization. This is not the case for the
Kaczmarz algorithm, where the work in each iteration consists of computing the
internal product $\langle A^{(i)},x^{(k)}\rangle$ and updating each entry of
the solution, a relatively small task compared to methods that require, for
example, matrix-vector multiplication. Thus, as the amount of work in each
iteration of the Kaczmarz algorithm is proportional to the number of columns
in matrix $A$, $n$, a block-sequential approach will be inadequate for systems
with small $n$.

For this reason, previously proposed parallel implementations focus on a
block-parallel strategy, and are reviewed next. It is also important to note
that these implementations were developed for sparse matrices, whereas our
implementation considers dense systems.

####  2.3.1 Randomized Kaczmarz with Averaging Method

The Randomized Kaczmarz with Averaging (RKA) method [5] was introduced by
Moorman, Tu, Molitor, and Needell with two goals: the first one was to speed
up the convergence of RK (while taking advantage of parallel computation); the
second was to decrease the convergence horizon when dealing with inconsistent
systems since the RK method does not converge to the least squares solution in
this case. In each iteration, multiple updates corresponding to different rows
are computed and then gathered and averaged before advancing to the next
iteration. If $q$ is the number of threads to use, then we need at each
iteration a set $\tau_{k}$ of $q$ rows randomly sampled. Then an iteration of
the RKA method can be written as

$x^{(k+1)}=x^{(k)}+\frac{1}{q}\sum_{i\in\tau_{k}}w_{i}\frac{b_{i}-\langle
A^{(i)},x^{(k)}\rangle}{\|A^{(i)}\|^{2}}{A^{(i)}}^{T}$ (5)

with $w_{i}$ representing the row weights, similar to the relaxation parameter
$\alpha_{i}$ in (3). The rows in set $\tau_{k}$ are chosen using the
probability distribution in (4). In [5] it is shown that RKA maintains the
linear convergence of RK. Moreover, when using multiple threads, RKA can
decrease the convergence horizon for inconsistent systems.

In the case of consistent systems and uniform weights, that is, using
$w_{i}=\alpha$, the authors computed the optimal value for $\alpha$:

$\alpha^{*}=\begin{cases}{\frac{q}{1+(q-1)s_{min}}},&s_{max}-s_{min}\leq{\frac{1}{q-1}}\\\
{\frac{2q}{1+(q-1)(s_{min}+s_{max})}},&s_{max}-s_{min}>{\frac{1}{q-1}}\end{cases}$
(6)

where $s_{min}=\sigma_{min}^{2}(A)/\|A\|_{F}^{2}$ and
$s_{max}=\sigma_{max}^{2}(A)/\|A\|_{F}^{2}$.

The introduction of parallel computing was made by noting that the projections
corresponding to each row in set $\tau_{k}$, that is, the calculations inside
the sum in (5), can and should be done in parallel. However, the authors do
not discuss an actual parallel implementation of the algorithm, meaning that
no results regarding speedups are available.

####  2.3.2 Component-Averaged Row Projections

The Component-Averaged Row Projections (CARP) method was introduced to solve
partial differential equations (PDEs) by Gordon and Gordon in [11]. Previous
methods used a block-sequential approach and are not adequate to solve PDEs:
some methods require computation and storage of data related to the inverses
of submatrices; others are very robust but do not have a consistent way to
partition the matrix into submatrices. CARP is a block-parallel method of the
Kaczmarz method where the blocks of equations are assigned to processors and
the results are then merged to create the next iteration. CARP is more memory
efficient than previous block-based projection methods since the only memory
requirement for each processor is the submatrix of the local equations
coefficients and respective entries of the $b$ vector and solution vector. In
terms of parallel behavior, CARP was developed in both shared and distributed
memory and it exhibits an almost linear speedup ratio. It is also important to
note that the optimal performance of CARP depends on the choice of the
relaxation parameter and the number of inner iterations in each block, which
are not known. Furthermore, this method was studied for systems derived from
PDEs which are sparse.

####  2.3.3 Asynchronous Parallel Randomized Kaczmarz Algorithm

Liu, Wright, and Sridhar [12] developed the Asynchronous Parallel Randomized
Kaczmarz (AsyRK) algorithm for shared memory systems. The asynchronous
parallel technique used in AsyRK was developed by Niu et al [13] for HOGWILD!,
a parallelization scheme for the Stochastic Gradient Descent (SGD) method
developed without locking. HOGWILD! was built specifically for sparse problems
to minimize the time spent in synchronization. Since SGD is an iterative
algorithm, similar to CARP or RKA, we can calculate several iterations in
parallel and then update the solution vector $x$. Normally we would think that
$x$ should be locked because, otherwise, the processors would override each
other’s updates. But, if the data is sparse, memory overwrites are minimal
since each processor only modifies a small part of $x$. This method reduces
significantly the overhead of synchronization; even if memory overwrites do
occur, the error introduced by them is very small and has minimal impact.
There is a small restriction regarding the update of $x$: there is a maximum
number of updates, $\rho$, that can occur between the time at which any
processor reads the current $x$ and the time at which it makes its update. The
asynchronous parallel variant of the Stochastic Gradient Descent method
exhibits near-linear speedup with the number of processors when using sparse
data. The application of this no-lock technique to the asynchronous
parallelization of the Stochastic Coordinate Descent (SCD) method [14], called
AsySCD [15], has also shown improvements in the rate of convergence when
compared to sequential versions.

The AsyRK algorithm was built by applying the HOGWILD! technique to RK. Liu,
Wright, and Sridhar show that the AsyRK method can outperform AsySCD in
running time by an order of magnitude. AsyRK includes the following steps:

  * •

The matrix $A$ is partitioned into blocks of equal size, such that a subset of
equations is assigned to each thread;

  * •

Each thread randomly orders its rows such that the rows in each iteration are
sampled without replacement, which was experimentally shown to have better
performance than sampling with replacement. After a full scan of a thread’s
submatrix, the thread reshuffles the rows.

  * •

Then, each thread repeats the following cycle: a row with index $i$ is
sampled; the quantity ${A^{(i)}}^{T}(b_{i}-\langle A^{(i)},x\rangle)$ is
computed; and the non-zero entries of $x$ are updated by a multiple of the
previously computed variable.

Although AsyRK has good results when compared to AsySCD, it is important to
remember that this method was developed for sparse matrices only and that
linear speedup is only observed if the number of processors fulfills certain
requirements related with the maximum eigenvalue of matrix $A^{T}A$.

##  3 Implementation and Results for the Parallelizations of the Kaczmarz
Method

In this section, we present several shared and distributed memory approaches
to the parallelization of RK using both block-sequential and block-parallel
strategies.

The organization of the section is as follows. In Section 3.1, we discuss the
experimental setup for our simulations. In Section 3.2, we discuss the
implementation and results of parallelization using a block-sequential
approach. The remaining parallelizations are block-parallel approaches. In
Section 3.3 we parallelize the RKA method presented in Section 2.3.1. In
Section 3.4, we propose a new block version of RKA. For these three sections,
we will use consistent systems. Finally, in Section 3.5 we discuss the
application of some of these parallelization approaches to solving
inconsistent systems.

###  3.1 Experimental Setup

Simulations were developed using the C++ programming language (version 11);
its source code and corresponding documentation are publicly available 111Code
available here: https://github.com/inesalfe/Review-Par-Kaczmarz.git.

The Kaczmarz method converges for the unique solution of consistent systems
but not for the least squares solution of inconsistent systems. However, the
RKA method (see Section 2.3.1) can decrease the convergence horizon of
inconsistent systems. For these reasons, we generated two data sets with
overdetermined systems: a consistent and inconsistent data set.

  * •

For the consistent data set, matrix entries were sampled from normal
distributions with average $\mu\in[-5,5]$ and standard deviation
$\sigma\in[1,20]$. The matrix with the largest dimension was generated by
using a different $\mu$ and $\sigma$ for each row and smaller-dimension
matrices were obtained by “cropping” the largest matrix. This keeps some
similarities between matrices of different dimensions for comparison purposes.
The solution $x$ is also sampled from the same probability distribution used
for matrix elements and the vector of constants $b$ is calculated as the
product of $A$ and $x$. Note that since matrix elements are random, the
systems will be full rank and the solution unique. The chosen matrix
dimensions were the following: $2000$, $4000$, $20000$, $40000$, $80000$, or
$160000$ for the number of rows and $50$, $100$, $200$, $500$, $750$, $1000$,
$2000$, $4000$, $10000$ or $20000$ for the number of columns.

  * •

For the inconsistent data set, we added an error term to the consistent
systems from the previous data set. Let $b$ and $b_{LS}$ be the vectors of
constants of the consistent and inconsistent systems. The latter was defined
such that $b_{LS}=b+\xi$, where $\xi$ is a random vector of the same size as
$b$, with entries given by a random number sampled from a gaussian
distribution $N(0,1)$. The least-squares solution, $x_{LS}$, was obtained
using the Conjugate Gradient method for Least-Squares (CGLS).

Since the Kaczmarz method is an iterative method, we need to define a stopping
criterion. Our goal here is to assess the performance of the iterative part of
the parallelized algorithms, ignoring the part of the algorithm corresponding
to the evaluation of the stopping criterion. To accomplish this we used the
following procedure: first, we determine the number of iterations, $k$, that
parallel implementations take to achieve a given error, that is,
$\|x^{(k)}-x^{*}\|^{2}<\varepsilon$; then we measure the runtime using that
previously calculated value as the maximum number of iterations. The value of
$\varepsilon$ controls the desired accuracy of the solution and was defined as
$\varepsilon=10^{-8}$ during our simulations.

Note that, in the Randomized Kaczmarz method, the sequence of sampled rows,
chosen according to the probability distribution in (4), will be different
according to the random number generator used. In our simulations, we used the
C++ class discrete_distribution
222https://en.cppreference.com/w/cpp/numeric/random/discrete_distribution as a
probability distribution and mt19937
333https://cplusplus.com/reference/random/mt19937/ as the random number
generator. The former requires a seed for initialization. Consequently, even
if we use the same initial estimate guess for the solution ($x^{(0)}=0$), runs
of the algorithms using generators with different seeds will require different
numbers of iterations to achieve the same stopping criterion. For this reason,
to get a robust estimate for the average number of iterations, for each input,
the algorithm is run 10 times with different seeds. As previously explained,
the execution time is then measured as the total time of 10 runs of the
algorithm using the average number of iterations as the maximum number of
iterations. The number 10 was chosen since it is enough to achieve a
percentual standard deviation of $1\%$ in the execution time.

The shared memory implementations of the several methods were accomplished
using the OpenMP444https://www.openmp.org/ C++ API (version 5.0). All
experiments were carried out on a computer with an AMD EPYC 9554P 64-Core
processor, running at 3.1GHz, and 256GB memory.

The distributed memory implementations were developed using the
MPI555https://www.open-mpi.org C++ API (version 4.1.1). All experiments were
carried out on the Navigator Cluster 666https://www.uc.pt/lca/computing-
resources/navigator-cluster/. We had access to a partition with 43 nodes, each
with two 2.7 GHz central processing units (Intel Xeon E5-2697 v2 CPU) with 12
cores and a 96GB memory.

###  3.2 Parallelization of Each Iteration

In a first attempt at parallelization, we show that using a block-sequential
approach that parallelizes the work inside each iteration does not always
exhibit speedup and, when it does, it is far from ideal. In this section, we
will only discuss the results of the shared memory implementation. Our
simulations use 1, 2, 4, 8, 16, and 64 threads.

As previously mentioned, the main computations in each iteration are the
calculation of the internal product, $\langle A^{(i)},x^{(k)}\rangle$, and the
update of the solution. The former can be effortlessly parallelized using the
OpenMP reduce command with the sum operation. The latter is easily handled by
distributing the entries of the solution by the available threads using the
OpenMP for command since the update of the entries of $x$ can be done
independently.

(a) Speedup for systems with $n=1000$.

(b) Speedup for systems with $n=10000$.

Figure 2: Speedups for the parallel implementation of RK using a block-
sequential approach using a fixed number of columns.

Figures 2a and 2b present the speedup defined as the quotient of the
sequential execution time and the parallel time using different numbers of
threads. Figure 2a shows that, regardless of the number of threads, there is
no speedup. Moreover, computational time increases when using more threads. In
Figure 2b, it is shown the speedup obtained for systems with a larger number
of columns. As expected the results are now better than the results for
systems with a smaller number of columns, since we are increasing the work
inside each iteration. However, the speedups are still far from the ideal
value. Note that this corresponds to a speedup equal to the number of threads.
Furthermore, we see that the parallelization is slower for 64 than it is for
16 threads, which means that the overhead of parallelization depends on the
number of threads.

Since it is clear that it is not possible to have an efficient implementation
of RK by parallelizing the computation of $\langle A^{(i)},x^{(k)}\rangle$, we
decided not to explore distributed memory implementations for this approach.

###  3.3 Randomized Kaczmarz with Averaging

We now move on to block-parallel implementations. In this section, we
implement the RKA algorithm (Section 2.3.1) in both shared and distributed
memory.

####  3.3.1 Implementation and Results for Shared Memory

In each iteration of RKA, each thread samples a row of the matrix, computes an
updated version of the estimate of the solution, and then the results for all
threads are averaged. In this implementation, we decided to use uniform row
weights ($w_{i}=\alpha$), meaning that we can rewrite (5) such that

$x^{(k+1)}=x^{(k)}+\frac{\alpha}{q}\sum_{i\in\tau_{k}}\frac{b_{i}-\langle
A^{(i)},x^{(k)}\rangle}{\|A^{(i)}\|^{2}}{A^{(i)}}^{T}\>,$ (7)

where $q$ is the number of threads. Note that, if $q=1$, we recover the RK
method. Using the previous equation, we implemented a sequential version of
RKA so that we could validate the results obtained by the authors of RKA. We
will now go through the details of the computations involved in each iteration
of the parallel implementation of RKA, presented in Algorithm 1. Note that
this represents a parallel section where all threads are running this code
concurrently.

Algorithm 1 Pseudocode for an iteration of the parallel implementation of RKA.
$\mathcal{D}$ is a probability distribution of a random variable taking the
row indices as values and with probabilities proportional to their norms
described by (4).

1:$\mathit{it}\leftarrow\mathit{it}+1$

2:omp barrier

3:omp for $i=0,...,n-1$ do

4: $x^{(prev)}_{i}\leftarrow x_{i}$

5:$\mathit{row}\leftarrow$ sampled from $\mathcal{D}$

6:$\mathit{scale}\leftarrow\alpha\times{\frac{b_{row}-\langle
A^{(row)},x^{(prev)}\rangle}{q\>\|A^{(row)}\|^{2}}}$

7:omp critical

8: for $i=0,...,n-1$ do

9: $x_{i}\leftarrow x_{i}+\mathit{scale}\times A^{(row)}_{i}$

In lines 3 and 4 of Algorithm 1, we store the estimate of the solution from
the previous iteration, $x^{(prev)}$. This is necessary since, if the scale
factor in line 6 was computed using the current estimate of the solution, $x$,
one thread could be computing the scale factor while another thread was
updating $x$ in the critical section ahead and the scale factor would not have
the correct value.

In line 5 of Algorithm 1, a row is sampled according to a probability
distribution with probabilities proportional to the norms of the rows of
matrix $A$. However, it does not make sense to have threads sampling the same
sequence of rows since we would be averaging identical results. This can be
easily avoided by giving each thread a different seed for the random number
generator that samples rows. Nonetheless, it is not only possible that
different threads sample the same row, but also likely if rows have very
different norms.

$0$$1$$2$$3$$4$$5$$6$$7$$8$$0$$1$$2$$3$$4$$5$$6$$7$$8$$0$$1$$2$$3$$4$$5$$6$$7$$8$$x^{(k+1)}_{\text{Thread
2}}$$x^{(k+1)}_{\text{Thread 1}}$$x^{(k+1)}_{\text{Thread 0}}$$x^{(k+1)}$
Computed by Thread 0 Computed by Thread 1 Computed by Thread 2

Figure 3: Averaging of the solution estimate using a matrix. Example for 3
threads and $n=9$.

In lines 7 to 9 of Algorithm 1, the results are combined. To ensure that all
threads update the estimate of the solution and that no two threads are
updating $x$ simultaneously, a critical section must be created, meaning that
the gathering of results is done sequentially. To try to introduce some
parallelism in this step, three other approaches were tested:

  * •

In the first, different threads start updating $x$ in a different entry, and
the atomic command is used. Note that the solution estimate is a shared
variable and data near the boundary between entries computed by two different
threads can fall in the same cache block. This causes many cache block
invalidations, consequently significantly increasing the average memory access
time, making this approach slower than with the critical section;

  * •

In the second we use the reduce command on the solution vector with the sum
operation, but this requires previously setting $x$ to zero, which can explain
why this approach is once more slower than the critical section;

  * •

In the third approach, instead of using the $x^{(prev)}$ array, we use a
matrix $x^{(k+1)}_{\text{Thread}}$, also a shared variable, that has $q$ rows
and $n$ columns. Each row of the matrix contains the solution estimate
computed by each thread. Figure 3 contains a scheme of the computation of the
solution estimate in a given iteration for 3 threads and a system with 9
columns. After each thread computes the solution estimate relative to the row
it sampled, the averaging process can be done in parallel. This is true since
threads can average the results for different entries simultaneously. However,
this approach requires an additional synchronization point before the
averaging. Moreover, similarly to the first approach, the update of the
solution estimate during the averaging process can be quite slow since there
are cache blocks that have entries that “belong” to different threads. For
these reasons, this approach didn’t prove to be faster than using the critical
section.

We can already identify two problems in the parallel implementation of this
algorithm: not only do we have a synchronization point in each iteration, but
also the averaging of results is done sequentially.

We now discuss the results for the RKA method, obtained with the same number
of threads as in the previous section (1, 2, 4, 8, 16, and 64). Two choices
for the uniform weight parameter $\alpha$ in (7) were made: $\alpha=1$ and
$\alpha=\alpha^{*}$, where $\alpha^{*}$ is the optimal parameter for RKA for
consistent systems, given by (6). The sequential and parallel versions of RKA
were tested for several numbers of threads (one thread is the RK method). We
start by showing how RKA can be used to accelerate the convergence of RK.

(a) Iterations.

(b) Speedup.

Figure 4: Results for RKA using 2, 4, 8, 16, and 64 threads and row weights
$\alpha=1$ for several overdetermined systems with $n=4000$ with a varying
number of rows.

Figure 4 shows the results for $\alpha=1$. It is clear, from Figure 4a that,
as expected, the number of iterations of RKA is inferior to that of RK.
Furthermore, when the number of threads is increased, fewer iterations are
needed for RKA to converge. However, the difference between the number of
iterations for two different values of $q$ gets smaller when $q$ increases
(for example, the decrease in iterations from 2 to 4 threads is larger than
for 4 to 8 threads). Figure 4b represents the speedup computed as the quotient
between the execution time of RK and the execution time of the parallel
implementation of RKA. The results show that, regardless of the number of
threads, RKA is a slower method than RK using $\alpha=1$. Note that, for the
RKA method, the work done by each thread during the computation of the results
corresponding to one row of matrix $A$ is the same as RK (lines 5 and 6 of
Algorithm 1). This means that the low speedup values for RKA can only be due
to the averaging of results. Although RKA requires fewer iterations to
converge than RK, this reduction is not enough to overcome the time spent in
updating $x^{(k)}$, which has to be done sequentially. Furthermore, the time
spent in updating the solution is proportional to the number of threads, which
explains why speedups decrease when more threads are used.

(a) Iterations.

(b) Speedup.

Figure 5: Results for RKA using 2, 4, 8, 16, and 64 threads and row weights
$\alpha=\alpha^{*}$ for several overdetermined systems with $n=4000$ with a
varying number of rows.

Figure 5 shows the results for $\alpha=\alpha^{*}$. Regarding the number of
iterations, we can see from Figure 5a that the reduction in the number of
iterations is much more significant when using the optimal values,
$\alpha=\alpha^{*}$, than when using unitary weights (Figure 4a). Furthermore,
that decrease in iterations is proportional to the number of threads used
(with the exception of 64 threads). This is translated into execution time,
since Figure 5b shows that speedups increase between 2 and 16 threads, unlike
the results for $\alpha=1$ (Figure 4b). However, for 64 threads the speedup is
smaller than for 16 threads. This is because the communication/synchronization
is more expensive for 64 threads than for 16 and the number of iterations is
not small enough to make up for averaging the results for 64 threads. In
summary, the speedups for the optimal row weights are better than the ones for
unitary row weights due to the significant decrease in iterations but the
values are still far from being ideal. Furthermore, it is important to
remember that we are ignoring the time needed to compute $\alpha^{*}$, which
is considerably high since the optimal values of the row weights require
computing the singular values of the matrix.

As previously mentioned, several threads may sample the same row in a given
iteration since all threads have access to the full matrix of the system. To
force threads to sample different rows, we developed a variation of the RKA
algorithm where matrix $A$ is distributed among threads. Therefore, the
parallel implementation of this algorithm can still be described by Algorithm
1 with the exception of line 5 where rows are sampled. To partition the system
between threads we used a distributed approach. For a given thread with
identifier $\mathit{t_{id}}$, its block of rows spans from index
$\mathit{low}=\lfloor t_{id}\frac{m}{q}\rfloor$ to index
$\mathit{high}=\lfloor(t_{id}+1)\frac{m}{q}\rfloor-1$, inclusively.

In the following, we propose another variation in the implementation of the
RKA algorithm: instead of using the optimal value $\alpha^{*}$, which requires
computing the singular values of the full matrix $A$, each thread can compute
a particular value for $\alpha$ resorting exclusively to the partial data of
the matrix they own. The submatrix $A$ used to compute $\alpha$ in each thread
is assigned by using the same distributed approach as before. This variation
was implemented since it is faster to have threads simultaneously computing
the optimal $\alpha$ of a submatrix than to have a single thread computing
$\alpha$ using the full matrix.

Table 1: Average number of iterations until convergence for a system
$40000\times 10000$. The numbers in parentheses correspond to the differences
in the number of iterations compared to the second column.

Threads | Full Matrix $\alpha$ | Partial Matrix $\alpha$
---|---|---
Full Matrix Access | Distributed Approach | Full Matrix Access | Distributed Approach
2 | $418718$ | $416164\>(-2554)$ | $418715\>(-3)$ | $416159\>(-2559)$
4 | $210381$ | $209061\>(-1320)$ | $210373\>(-8)$ | $209053\>(-1328)$
8 | $105499$ | $106083\>(584)$ | $105490\>(-9)$ | $106076\>(577)$
16 | $53576$ | $53742\>(166)$ | $53660\>(84)$ | $53825\>(249)$

In summary, we have two ways of sampling rows: threads sample rows from the
entire matrix or threads sample rows from nonoverlaping submatrices (Full
Matrix Access or Distributed Approach in Table 1); and we have two ways of
computing parameter $\alpha$: using the optimal value $\alpha^{*}$ for all
threads computed with the entire matrix or having a different $\alpha$ value
for each thread, computed using different submatrices (Full Matrix $\alpha$ or
Partial Matrix $\alpha$ in Table 1). Table 1 shows the average number of
iterations until convergence for a system $40000\times 10000$ for the four
possible scenarios. Note that the second column corresponds to the scheme used
to obtain the results in Figure 5. It is clear that, whether threads can
sample rows from the entire matrix or not, there is little difference in the
number of iterations using the optimal $\alpha$ computed with the total matrix
or $\alpha$ computed with submatrices (columns 2 and 4 are very similar and so
are columns 3 and 5). However, that difference increases when we increase the
number of threads, meaning that the singular values in the submatrices get
further away from the singular values of the full matrix for larger numbers of
threads. Regarding the difference between threads accessing the entire matrix
or not when sampling rows note that (comparison between columns 2 and 3 or 4
and 5), independently of how $\alpha$ is calculated, for smaller numbers of
threads, it is better to have threads choose rows from disjoints parts of the
matrix. However, this difference is small (the percentual difference in
iterations is less than 1%). Furthermore, for a larger number of threads, the
effect is the opposite, and it is better to have threads sample rows from the
entire matrix. Note that the submatrices owned by each thread get smaller when
the number of threads increases. Therefore, it is reasonable that, when using
more threads, the Distributed Approach requires more iterations than the Full
Matrix Access option since there is a relationship between the number of rows
available for sampling and information. For a more in-depth analysis of the
relationship between rows and information see Section 5.1.1 of [10].

From this analysis, we can conclude that, if time is an important factor when
computing the optimal value for $\alpha$, we can use a different $\alpha$ for
each thread computed using a submatrix of the original problem without having
a noticeable difference in the number of iterations. Moreover, having threads
sample rows from disjoint parts of the matrix can slightly improve or worsen
the convergence of the algorithm, depending on the number of threads used.
This analysis will be relevant for the distributed memory implementation of
RKA.

####  3.3.2 Implementation and Results for Distributed Memory

Since one of the advantages of using distributed memory is to be able to
process data sets that are not able to be stored in a single machine, for the
distributed memory implementation of RKA we divide matrix $A$ and vector $b$
between the several available machines. Recall that, from the analysis of the
shared memory approach, we concluded that there is no major difference in
terms of the number of iterations between partitioning the system or not among
processes. The details in the implementation of a single iteration of the
parallel version of RKA for distributed memory are shown in Algorithm 2.

Algorithm 2 Pseudocode for an iteration of the parallel distributed
implementation of RKA. $\mathcal{D}$ is the aforementioned probability
distribution of a random variable taking the row indices as values and with
probabilities proportional to their norms described by (4).

1:$\mathit{it}\leftarrow\mathit{it}+1$

2:$\mathit{row}\leftarrow$ sampled from $\mathcal{D}$

3:$\mathit{scale}\leftarrow\alpha\times{\frac{b_{row}-\langle
A^{(row)},x^{(prev)}\rangle}{\|A^{(row)}\|^{2}}}$

4:for $i=0,...,n-1$ do

5: $x_{i}\leftarrow{\frac{x_{i}+\mathit{scale}\times A^{(row)}_{i}}{np}}$

6:mpi Allreduce ($x,+$)

Note that Algorithm 2 is much simpler than the implementation of RKA for
OpenMP (Algorithm 1). This is due to the fact that in distributed memory we do
not need to worry about conflicts between processes reading and writing in the
same memory position. This eliminates the need for saving the estimated
solution from the previous iteration ($x^{(prev)}$ in Algorithm 1). The
averaging of the results was accomplished with the Allreduce command with the
sum operation. In OpenMP we had to update the solution sequentially using a
critical section. The dependency of the communication time on the number of
processes in MPI is different than the relationship between the communication
time and the number of threads in OpenMP. In OpenMP the averaging of the
results takes $O(q)$ time, where $q$ is the number of threads, while here it
takes $O(log(np))$, where $np$ is the number of processes, since the mpi
Allreduce operation is implemented with a hypercube topology. However,
communication time has other dependencies. In MPI, the time required to send
messages between processors depends on the bandwidth, that is, the amount of
data that can be transmitted over the network in one unit of time. In OpenMP,
communication time depends on how fast threads can access memory. This means
that, although communication in MPI has a smaller dependency on the number of
processes, communication between processes is usually slower than when using
OpenMP.

Since a single node of the cluster has 2 central processing units with 12
cores each, we experimented with 2 configurations for distributed memory. In
the first, we use all 24 cores in each node and have one process per core. In
the second option, we have two mpi processes per node. Our simulations were
run for 1, 2, 4, 8, 12, 24, and 48 processes, meaning that we need at most 2
nodes to run the first configuration and 24 nodes to run the second
configuration.

(a) Speedup for overdetermined systems with $n=1000$.

(b) Speedup for overdetermined systems with $n=10000$.

Figure 6: Results for RKA using 2, 4, 8, 12, 24, and 48 processes and row
weights $\alpha=\alpha^{*}$ using two distributed memory configurations.

Figure 6 shows the results for the two discussed process/node configurations.
Note that, for smaller systems, shown in Figure 6a, maximizing the number of
processes per node is a faster option than distributing processes across
several nodes. This is due to the cost of communication being higher for
processes on different nodes than for processes on the same node.

However, the same conclusion cannot be drawn for larger systems, shown in
Figure 6b. From 2 to 12 processes, there is little difference between the two
configurations. For 24 processes, two processes per node is a faster option.
Note that the submatrix owned by each process does not fit in the L1 cache,
meaning that processes on the same node contend for entries in the L3 cache.
This effect is more noticeable for larger systems for which there is a lower
probability that the sampled rows are stored in the cache. Therefore, although
communication inside the same node is faster than between nodes, the memory
access time has a larger effect than communication time for larger systems,
This effect is more noticeable when more processes are being used, explaining
the results for 24 processes. For 48 processes, even when utilizing all CPUs
in a node, we need to use two nodes, which explains why the configuration with
2 processes per node is slower than the other configuration option.
Furthermore, note that the speedups for 48 processes in Figures 6a and 6b is
smaller than the speedups for 24 processes. Similarly to what happened for the
shared memory implementation, the number of iterations for 48 processes is
smaller than for 24 but it is not small enough to make up for the increased
cost in communication when using more processes. In summary, depending on the
dimension of the problem, one might benefit from using several nodes even if
we are not using all the CPUs in that node. Regardless of the configuration
used, similar to the shared memory implementation, the speedup results are
poor.

Our conclusion is that it is not possible to efficiently parallelize RKA due
to the large overhead in synchronization/communication between
threads/processes, especially when averaging the results. For efficient
parallelization, we need to find a way to average the results in parallel, or
to decrease the number of times that the results are averaged and, therefore,
decrease the impact of communication. We explore this last approach next.

###  3.4 Randomized Kaczmarz with Averaging with Blocks

To decrease the impact of communication we developed a new variant of the RKA
method called the Randomized Kaczmarz with Averaging with Blocks (RKAB)
method.

####  3.4.1 Description of the Method

In a single iteration of RKA, each task only makes use of one row of the
matrix before the results are averaged. In RKAB, instead of using only one row
per iteration, each process computes the results corresponding to a block of
several rows, meaning that the results from several processes are only
gathered every once in a while. We can already identify one advantage and one
disadvantage of RKAB in regards to RKA: on one hand, communication will happen
much more sporadically; on the other hand, since results are only shared
between processes after processing a block of rows, there is less shared
information, which might lead to a slower convergence in regards to the number
of iterations.

We now describe how to update the solution in a given iteration $k$ of the
RKAB method. We start by defining $v_{\gamma}$ as the solution estimate in
thread $\gamma$. In each iteration, this variable is set to the solution
estimate given in the previous iteration, meaning that,
$v_{\gamma}^{(0)}=x^{(k)}$. Then, each task samples a given number of rows,
defined by parameter $\mathit{block\>size}$, and updates its solution
estimate, such that

$v_{\gamma}^{(j+1)}=v_{\gamma}^{(j)}+\alpha_{i}\>\frac{b_{i}\>-\langle
A^{(i)},v^{(j)}\rangle}{\|A^{(i)}\|^{2}}\>{A^{(i)}}^{T}\>,$ (8)

with $j$ starting at $0$ and ending at iteration $\mathit{block\>size}-1$,
which we denote by $bs$. Similarly to RK, rows are sampled using the
probability distribution (4). Finally, we average the results obtained in each
thread with

$x^{(k+1)}=\frac{1}{q}\sum_{\gamma=1}^{q}v_{\gamma}^{(\mathit{bs+1})}\>.$ (9)

The RKAB is distinct from the CARP method (Section 2.3.2) for several reasons:
first, the CARP method is a variation of the Kaczmarz method and the RKAB is a
variation of the Randomized Kaczmarz method, meaning that rows are sampled
according to different criteria; second, in the RKAB method, we average the
results from all threads but, in the CARP method, since it is used for sparse
matrices, this is not necessary since threads only change parts of the
solution.

####  3.4.2 Implementation and Results for Shared Memory

Algorithm 3 Pseudocode for an iteration of the parallel implementation of
RKAB. $\mathcal{D}$ is the aforementioned probability distribution of a random
variable taking the row indices as values and with probabilities proportional
to their norms described by (4).

1:$\mathit{it}\leftarrow\mathit{it}+1$

2:omp barrier

3:$\mathit{row}\leftarrow$ sampled from $\mathcal{D}$

4:$\mathit{scale}\leftarrow\alpha\times{\frac{b_{row}-\langle
A^{(row)},x\rangle}{\|A^{(row)}\|^{2}}}$

5:for $i=0,...,n-1$ do

6: $v_{i}\leftarrow x_{i}+\mathit{scale}\times A^{(row)}_{i}$

7:for $b=0,...,\mathit{block\>size}-1$ do

8: $\mathit{row}\leftarrow$ sampled from $\mathcal{D}$

9: $\mathit{scale}\leftarrow\alpha\times{\frac{b_{row}-\langle
A^{(row)},v\rangle}{\|A^{(row)}\|^{2}}}$

10: for $i=0,...,n-1$ do

11: $v_{i}\leftarrow v_{i}+\mathit{scale}\times A^{(row)}_{i}$

12:for $i=0,...,n-1$ do

13: $v_{i}\leftarrow v_{i}-x_{i}$

14:omp barrier

15:omp critical

16: for $i=0,...,n-1$ do

17: $x_{i}\leftarrow x_{i}+{\frac{v_{i}}{q}}$

A detailed explanation of the work inside a single iteration of the parallel
shared memory implementation of RKAB is presented in Algorithm 3. Here,
instead of saving the previous iteration (Algorithm 1), every thread has a
private variable $v$ that stores the current results for that thread. Again,
all threads execute this same code in parallel.

When starting a new block, threads use the estimate of the solution from the
previous iteration (lines 3 to 6 of Algorithm 3). For the remaining rows of
the block, threads use their local estimative of the solution, $v$, to compute
the scale factor (lines 7 to 11 of Algorithm 3). Note that the size of the
block, $\mathit{block\>size}$, has to be determined by the user, and using the
RKAB method with $\mathit{block\>size}=1$ is equivalent to the RKA method.
Later on, we will discuss how to choose $\mathit{block\>size}$.

The process of averaging the results is slightly different from RKA. In each
thread, after the results of the entire block are computed, we subtract $x$ to
$v$ so that we can update $x$ by summing the differences. Another option would
be setting $x$ to $0$ and summing the full results but the computational
effort is the same. Note that we need the barrier in line 14 of Algorithm 3 to
avoid having one thread updating $x$ in the critical region while another
thread is in the process of computing $v$ in line 13.

During the implementation of RKA, we used uniform row weights, $\alpha$, and
analyzed the results using $\alpha=1$ and $\alpha=\alpha^{*}$, where
$\alpha^{*}$ are the optimal parameters for consistent systems. However, RKAB
is a different algorithm from RKA and there is no a priori estimated optimal
value for $\alpha$. For that reason, we will first analyze the dependency of
RKAB on other parameters such as the $\mathit{block\>size}$ using $\alpha=1$.
After this analysis, we will evaluate how several choices of row weights,
$\alpha$, can impact the performance of RKAB.

(a) Iterations for a system $80000\times 1000$.

(b) Number of lines used for a system $80000\times 1000$.

(c) Total computational time for a system $80000\times 1000$.

Figure 7: Results for the parallel implementation of RKAB using 2, 4, 8, 16,
and 64 threads for an overdetermined system $80000\times 1000$ for several
values of $\mathit{block\>size}$.

We will start by analyzing how the execution time of RKAB depends on the size
of the blocks by setting
$\mathit{block\>size}=\\{5,10,100,500,1000,2000,4000,10000\\}$. Similarly to
the previous parallelization attempts, the parallel implementation of RKAB was
tested for 1, 2, 4, 8, 16, and 64 threads. Figure 7 shows the results for RKAB
for a system with dimensions $80000\times 1000$. Figure 7a shows that, for all
threads, when we increase $\mathit{block\>size}$, the number of iterations
decreases. This was expected since, by processing a larger number of rows in
each iteration, the estimate of the solution will converge faster.
Furthermore, for fixed $\mathit{block\>size}$, larger numbers of threads
require fewer iterations. However, as in the case of the RKA algorithm with
$\alpha=1$ (Figure 4a), this reduction is small. As a result, the total number
of used rows, shown in Figure 7b, increases accordingly to the number of
threads involved (the total number of used rows is given by the product of the
number of iterations, threads, and $\mathit{block\>size}$). Figure 7b also
shows that, for a given number of threads, the total number of used rows
remains stable until reaching $\mathit{block\>size}=n=1000$, after which it
increases. The results in Figure 7c can be easily explained by using the
number of iterations and the total number of rows. First, for blocks of a
fixed size, the computational time is similar for different numbers of threads
except for 64 threads. This is because the overhead in
communication/synchronization, due to the larger number of threads, overcomes
the observed reduction in the number of iterations. Second, increasing
$\mathit{block\>size}$ reduces the computational time. This is because of the
smaller number of times that the threads have to communicate to average the
results, despite the total amount of work being similar between
$\mathit{block\>size}$s (see Figure 7b). However, note that for larger blocks
($\mathit{block\>size}>n$) the computational time increases. For a full rank
matrix, having the same number of rows as columns is enough information to
solve the system and find its unique solution. From this observation, we can
conclude that it is useless to use $\mathit{block\>size}$ larger than $n$,
since the estimated solution obtained from each thread is already close to the
true solution making it pointless to average similar results.

(a) Total computational time for a system $80000\times 4000$.

(b) Total computational time for a system $80000\times 10000$.

Figure 8: Results for the sequential version of RK and the parallel
implementation of RKAB using 2, 4, 8, 16, and 64 threads for several
overdetermined consistent systems.

To further validate the conclusions regarding $\mathit{block\>size}$, Figure 8
shows the results for a system with other numbers of columns and it also shows
the sequential time of RK. Once more, for a system $80000\times 4000$ (Figure
8a), time slightly increases when we go from $\mathit{block\>size}=4000$ to
$\mathit{block\>size}=10000$, especially for lower number of threads. However,
for the system with $80000\times 10000$ (Figure 8b), time does not increase
for the larger $\mathit{block\>size}$. In summary, we can use the number of
columns as a rule of thumb to select the $\mathit{block\>size}$. Still, just
like RKA, depending on the number of threads and, in this case, also on
$\mathit{block\>size}$, the parallel implementation might not be faster than
RK and, when it is, the difference in execution time is small.

In Section 3.3, we analyzed the difference in the number of iterations of the
RKA method between threads sampling rows from any part of matrix $A$ (Full
Matrix Access) or distributing the matrix such that threads sample rows from
disjoint parts of the matrix (Distributed Approach). Here we will make the
same analysis for the RKAB method. Figure 9 shows the results for the two
sampling schemes for a system with dimension $40000\times 10000$. Figure 9a
shows a noticeable difference in iterations between the two methods for larger
$\mathit{block\>sizes}$, especially for 64 threads. That difference in
iterations is even more visible when we consider the total number of used
rows, shown in Figure 9b. Consequently, this translates into the execution
time (Figure 9c) since the time for the Distributed Approach starts to
increase for smaller values of $\mathit{block\>size}$ while the time using
Full Matrix Access is stable. These results show that we cannot use the same
rule of thumb for the optimal $\mathit{block\>size}$ when using a Distributed
Approach. In the next section, where we discuss the distributed memory
approach for this method, we include a more detailed discussion of this topic.

(a) Average number of iterations.

(b) Total number of rows used.

(c) Execution time.

Figure 9: Comparison between two different sampling schemes for a system
$40000\times 10000$.

We now discuss how different row weights, given by parameter $\alpha$, can
affect the performance of RKAB. For each number of threads, several values of
$\alpha$ were chosen between $1$ and the optimal value of $\alpha$ for RKA,
$\alpha^{*}$, described by (7). Since the values of $\alpha^{*}$ for a system
$80000\times 1000$ are $1.999$ and $3.992$ for 2 and 4 threads, we choose
$\alpha=\\{1.0,1.2,1.3,1.5,1.8,1.999\\}$ for 2 threads and
$\alpha=\\{1.0,1.5,2.0,2.5,3.0,3.991\\}$ for 4 threads, so that we have more
or less evenly spaced test values between $1$ and $\alpha^{*}$. Figure 10
shows the results for the system with dimensions $80000\times 1000$. Note
that, for 4 threads, there are a few values of $\alpha$ for which the number
of iterations is not shown (for example, using $\mathit{block\>size}=500$ and
$\alpha=3.0$). The results for these values are not shown because the RKAB
does not converge. This conclusion was drawn by realizing that the error
$\|x^{(k)}-x^{*}\|$ grows with the number of iterations, meaning that the
estimated solution diverges from $x^{*}$. Figure 10a shows that, for 2
threads, $\alpha^{*}$ is not the optimal value of $\alpha$ for RKAB. Although,
regardless of $\mathit{block\>size}$, RKAB converges for $\alpha^{*}$, there
are other values of $\alpha$ for which the number of iterations is lower. It
also appears that the optimal value of $\alpha$ for RKAB decreases with the
increase of $\mathit{block\>size}$. For 4 threads (Figure 10b), not only is
$\alpha^{*}$ not the optimal value for $\alpha$, but also, depending on the
$\mathit{block\>size}$, RKAB might not even converge to this value. Just like
for 2 threads, the optimal $\alpha$ seems to decrease when the
$\mathit{block\>size}$ is increased.

(a) Number of iterations as a function of $\alpha$ using 2 threads for
different $\mathit{block\>size}$s.

(b) Number of iterations as a function of $\alpha$ using 4 threads for
different $\mathit{block\>size}$s.

Figure 10: Number of iterations for the RKAB method as a function of parameter
$\alpha$ for a system with dimensions $80000\times 1000$.

In summary, RKA may converge for a given $\alpha$ while RKAB does not; the
optimal values of $\alpha$ for RKA are not the same as the ones for RKAB; the
optimal value of $\alpha$ for RKAB depends on $\mathit{block\>size}$.

We now compare the performance of the RKAB method with the RKA method. In
Table 2 we show the execution times using different numbers of threads for
RKAB and RKA. Note that the RKAB method is always faster than the RKA using
$\alpha=1$. The same cannot be said when comparing with the RKA method with
the optimal $\alpha$ parameter. However, if we take into account the computing
time of $\alpha^{*}$, using the RKAB method with $\alpha=1$ is a faster
option. In conclusion, if $\alpha^{*}$ is not available, the RKAB is a faster
method than the RKA method. Nonetheless, neither of these methods can
consistently beat the sequential RK algorithm.

Table 2: Execution times of the several algorithms in seconds for a system
$80000\times 10000$. The last column corresponds to the computation time of
the optimal $\alpha$ parameter for the RKA algorithm. For the RKAB algorithm,
we used $\mathit{block\>size}=n$. The execution time of the sequential
version, the RK method, is $50$ seconds.

Threads | RKAB ($\alpha=1$) | RKA ($\alpha=1$) | RKA ($\alpha=\alpha^{*}$) | Computing $\alpha^{*}$
---|---|---|---|---
2 | $47$ | $157$ | $89$ | $2500$
4 | $46$ | $235$ | $66$ | $2514$
8 | $46$ | $322$ | $46$ | $2513$
16 | $47$ | $389$ | $31$ | $2506$
64 | $81$ | $883$ | $45$ | $2503$

####  3.4.3 Implementation and Results for Distributed Memory

Similarly to the distributed memory implementation of RKA (Section 3.3.2), the
system will be partitioned between the several processes. Regarding the
parallel implementation of RKAB for distributed memory, the details of a
single iteration of RKAB for distributed memory are shown in Algorithm 4.

Algorithm 4 Pseudocode for an iteration of the parallel distributed
implementation of RKAB. $\mathcal{D}$ is the aforementioned probability
distribution of a random variable taking the row indices as values and with
probabilities proportional to their norms described by (4).

1:$\mathit{it}\leftarrow\mathit{it}+1$

2:for $b=0,...,\mathit{block\>size}-1$ do

3: $\mathit{row}\leftarrow$ sampled from $\mathcal{D}$

4: $\mathit{scale}\leftarrow\alpha\times{\frac{b_{row}-\langle
A^{(row)},x\rangle}{\|A^{(row)}\|^{2}}}$

5: for $i=0,...,n-1$ do

6: $x_{i}\leftarrow x_{i}+\mathit{scale}\times A^{(row)}_{i}$

7:$\mathit{row}\leftarrow$ sampled from $\mathcal{D}$

8:$\mathit{scale}\leftarrow\alpha\times{\frac{b_{row}-\langle
A^{(row)},x\rangle}{\|A^{(row)}\|^{2}}}$

9:for $i=0,...,n-1$ do

10: $x_{i}\leftarrow{\frac{x_{i}+\mathit{scale}\times A^{(row)}_{i}}{np}}$

11:mpi Allreduce ($x,+$)

Algorithm 4 is much simpler than the parallel implementation for RKAB using
shared memory (Algorithm 3) since we do not need the extra variable $v$.

Similarly to the analysis of the distributed memory implementation of the RKA
method, for the RKAB method we also tested two process/node configurations.

(a) System with dimension $80000\times 1000$.

(b) System with dimension $80000\times 10000$.

Figure 11: Execution time for the distributed memory implementation of the
RKAB method for two different systems using two configurations.

Figure 11 shows the results for two systems. We will first analyze the
difference between the two configurations and then the relationship between
time and $\mathit{block\>size}$.

For the smaller system (Figure 11a) and smaller $\mathit{block\>size}$s, the
configuration using two processes per node is the slower option. Similarly to
RKA, this can be explained since the cost of communication is cheaper for
processes in the same nodes than between nodes. However, two processes per
node are the faster option for larger $\mathit{block\>size}$s. Note that, when
the $\mathit{block\>size}$ is increased, the number of iterations decreases
and, consequently, so does the amount of times that the processes communicate
between themselves. Therefore, the communication part of the runtime of the
algorithm is negligible compared with the computation part, and the memory
access time will have more impact than the communication time between
processes. When more processes are used, there is more competition to access
memory since these systems cannot be stored in the cache. This effect can also
be seen in Figure 11b where, regardless of the $\mathit{block\>size}$, it is
always faster to use two processes per node.

Regarding the optimal choice for $\mathit{block\>size}$, just like we
discussed in the previous section, using the number of columns of the system
for the $\mathit{block\>size}$ is not a good choice when distributing the
matrix by the processes. To explain why this is the case we will use an
example. Suppose that we are studying these two systems using 40 processes.
The size of the submatrix in each process for the smaller and larger systems
will be $2000\times 1000$ and $2000\times 10000$, respectively. Note that, for
the smaller system, we have an overdetermined system and, for the second
system, we will have an underdetermined system. For the first system, using
$\mathit{block\>size}$ to be equal to the number of columns should be enough
for us to get a solution close to the solution of the entire system in each
process. However, the same is not necessarily true for the second system since
we have fewer rows than columns. Furthermore, if we use
$\mathit{block\>size}=n=10000$, we will have to use the same row more than
once and reuse information. In summary, the choice of the
$\mathit{block\>size}$ for the distributed memory approach will depend not
only on the size of the system in each process but also on the relationship
between the number of rows and columns of that subsystem. Further
investigation into this topic is necessary to find a systematic way to choose
$\mathit{block\>size}$.

###  3.5 Application of RKA and RKAB to Inconsistent Systems

In Sections 3.3.1 and 3.4.2 we discussed the parallel implementations of RKA
and RKAB. We concluded that, in general, neither RKA nor RKAB can consistently
beat the sequential RK in terms of execution time. However, RKA can decrease
the convergence horizon for inconsistent systems when more than one thread is
used, something that is not possible for RK (Section 2.3.1). In this section,
we show that RKAB is also able to achieve this, which is relevant since the
RKAB method can be faster than the RKA method.

To analyze how the convergence horizon can change for several numbers of
threads, we will show the evolution of the norm of the error,
$\|x^{(k)}-x_{LS}\|$, and the norm of the residual $\|Ax^{(k)}-b\|$. To obtain
these variables, we ran the RKA and RKAB algorithms with a maximum number of
iterations and stored the error and residual norm every few iterations. That
interval of iterations is given by variable $\mathit{step}$. The inconsistent
system used in this section has dimensions $80000\times 1000$.

(a) Error norm, given by $\|x^{(k)}-x_{LS}\|$.

(b) Residual norm, given by $\|Ax^{(k)}-b\|$.

Figure 12: Results for RKA (with $\alpha=1$) for a system $80000\times 1000$.
Here we show the first 30000 iterations and the error and residual were stored
every $\mathit{step}=100$ iterations.

Figure 12 shows the error and residual evolution for the RKA algorithm using a
single row weight, $\alpha=1$. Figure 12a shows that using a higher number of
threads, $q$, decreases the error value upon which the error stabilizes.
Figure 12b shows that the residual for $q=20$ and $q=50$ stabilizes around the
residual value for the least-squares solution. However, this does not mean
that the solution obtained using these numbers of threads is the least-squares
solution since the error is not zero.

(a) Error norm, given by $\|x^{(k)}-x_{LS}\|$.

(b) Residual norm, given by $\|Ax^{(k)}-b\|$.

Figure 13: Results for RKA (with $\alpha=\alpha^{*}$) for a system
$80000\times 1000$. Here we show the first 30000 iterations and the error and
residual were stored every $\mathit{step}=100$ iterations.

Figure 13 shows the error and residual evolution for the RKA algorithm using
the optimal row weights for consistent systems, $\alpha=\alpha^{*}$. Note
that, although it is not expected for $\alpha^{*}$ to be optimal for these
systems since they are inconsistent, we would like to analyze the impact of
different $\alpha$ values in the convergence horizon. Note that increasing $q$
only decreases the final error for $q=50$. Moreover, for $\alpha=1$ (Figure
12a), it was guaranteed that increasing $q$ decreases the minimum error value,
something that is not true for $\alpha=\alpha^{*}$. Figure 13b also shows that
the residual only decreases when using $q=50$, but that decrease is smaller
than the one for $\alpha=1$ (Figure 12b). Although using $\alpha=\alpha^{*}$
may or may not be helpful in decreasing the convergence horizon (it depends on
the value of $q$), Figures 13a and 13b show that the error and residual
stabilize in less iterations than for $\alpha=1$.

(a) Error norm, given by $\|x^{(k)}-x_{LS}\|$.

(b) Residual norm, given by $\|Ax^{(k)}-b\|$.

Figure 14: Results for RKAB (with $\alpha=1$) for a system $80000\times 1000$.
Here we show the first 30 iterations and the error and residual were stored
every $\mathit{step}=1$ iterations.

Figure 14 shows the error and residual evolution for the RKAB algorithm using
block sizes equal to the number of columns of the system
($\mathit{block\>size}=1000$) and unitary row weights ($\alpha=1$). Figure 14
is very similar to Figure 12 in terms of the relationship between the error
and residual for increasing numbers of threads, which shows that the RKAB
method, like RKA, can decrease the convergence horizon. RKAB requires fewer
iterations to converge since the work per iteration is much more significant
than the work in one iteration of RKA.

In conclusion, for inconsistent systems, the RKA and RKAB methods can be used
to decrease the convergence horizon. For RKA, unitary row weights ($\alpha=1$)
can decrease the convergence horizon more than using the optimal row weights
for consistent systems ($\alpha=\alpha^{*}$). However, the optimal values for
consistent systems $\alpha=\alpha^{*}$ can accelerate the method’s
convergence, that is, the error and residual stabilize in fewer iterations
than when using $\alpha=1$. The RKAB method (with $\mathit{block\>size}=n$)
has the same effect in decreasing the convergence horizon as the RKA method
when both methods use unitary row weights ($\alpha=1$).

##  4 Conclusion

In this paper, we explored multiple approaches to parallelize the Randomized
Kaczmarz method using large-scale dense systems. We conclude that, either
using shared or distributed memory, is not possible to parallelize this method
efficiently for these types of systems. More specifically, we implemented the
Randomized Kaczmarz with Averaging (RKA) method. We showed that the parallel
algorithm is not efficient at all. This is because the high overhead of
synchronization/communication outweighs the clear advantage of effectively
reducing the number of iterations when compared with RK. To overcome these
overheads, we propose here a new block version of the averaging method called
the Randomized Kaczmarz with Averaging with Blocks (RKAB) method. This method
is capable, as it was the former RKA, of reducing the convergence horizon for
inconsistent systems. Although the parallel implementation of this method
cannot consistently surpass the execution times of the sequential RK, it is
faster than RKA when unitary row weights are used (that is, setting the
relaxation parameter to 1). For the shared memory implementation of RKAB, we
showed that a good choice for the block size parameter is to use the number of
columns of the matrix of the system.

In summary, if we are dealing with an inconsistent system and the goal is not
to find the minimum error solution but rather to regularize the solution
filtering out the underlying noise, the RKAB method is a better option than
the RKA method since we must consider the time used to compute the optimal
$\alpha$ parameter for RKA. This is a recurrent problem in many practical
applications that use real-world measurements, for example when trying to
reconstruct images of scanned bodies during a Computed Tomography.

## Acknowledgements

This work was supported by national funds through FCT, Fundação para a Ciência
e a Tecnologia, under projects URA-HPC PTDC/08838/2022 and UIDB/50021/2020
(DOI:10.54499/UIDB/50021/2020). JA was funded by Ministerio de Universidades
and specifically the requalification program of the Spanish University System
2021-2023 at the Carlos III University.

## References

  * [1] K. D. Gremban, C. E. Thorpe, T. Kanade, Geometric camera calibration using systems of linear equations, in: Proceedings. 1988 IEEE International Conference on Robotics and Automation, IEEE, 1988, pp. 562–567.
  * [2] P. C. Hansen, J. Jørgensen, W. R. Lionheart, Computed tomography: algorithms, insight, and just enough theory, SIAM, 2021.
  * [3] X. Chen, The kaczmarz algorithm, row action methods, and statistical learning algorithms, Frames and harmonic analysis 706 (2018) 115–127.
  * [4] 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).
  * [5] J. D. Moorman, T. K. Tu, D. Molitor, D. Needell, Randomized kaczmarz with averaging, BIT Numerical Mathematics 61 (1) (2021) 337–359.
  * [6] A. Ma, D. Needell, A. Ramdas, Convergence properties of the randomized extended gauss–seidel and kaczmarz methods, SIAM Journal on Matrix Analysis and Applications 36 (4) (2015) 1590–1604.
  * [7] T. Wallace, A. Sekmen, Deterministic versus randomized kaczmarz iterative projection, arXiv preprint arXiv:1407.5593 (2014).
  * [8] T. Strohmer, R. Vershynin, A randomized kaczmarz algorithm with exponential convergence, Journal of Fourier Analysis and Applications 15 (2007) 262–278.
  * [9] D. Needell, Randomized kaczmarz solver for noisy linear systems, BIT Numerical Mathematics 50 (2) (2010) 395–403.
  * [10] I. A. Ferreira, J. A. Acebrón, J. Monteiro, Survey of a class of iterative row-action methods: The kaczmarz method, arXiv preprint arXiv:2401.02842 (2024).
  * [11] D. Gordon, R. Gordon, Component-averaged row projections: A robust, block-parallel scheme for sparse linear systems, SIAM Journal on Scientific Computing 27 (3) (2005) 1092–1117.
  * [12] J. Liu, S. J. Wright, S. Sridhar, An asynchronous parallel randomized kaczmarz algorithm, arXiv preprint arXiv:1401.4780 (2014).
  * [13] B. Recht, C. Re, S. Wright, F. Niu, Hogwild!: A lock-free approach to parallelizing stochastic gradient descent, Advances in neural information processing systems 24 (2011).
  * [14] Y. Nesterov, Efficiency of coordinate descent methods on huge-scale optimization problems, SIAM Journal on Optimization 22 (2) (2012) 341–362.
  * [15] J. Liu, S. Wright, C. Ré, V. Bittorf, S. Sridhar, An asynchronous parallel stochastic coordinate descent algorithm, in: International Conference on Machine Learning, PMLR, 2014, pp. 469–477.


# The partition function modulo 4

Ken Ono  Department of Mathematics, University of Virginia, Charlottesville,
VA. 22904<EMAIL_ADDRESS>
###### Abstract.

It is widely believed that the parity of the partition function $p(n)$ is
“random.” Contrary to this expectation, in this note we prove the existence of
infinitely many congruence relations modulo 4 among its values. For each
square-free integer $1<D\equiv 23\pmod{24},$ we construct a weight 2
meromorphic modular form that is congruent modulo 4 to a certain twisted
generating function for the numbers
$p\big{(}\frac{Dm^{2}+1}{24}\big{)}\pmod{4}$. We prove the existence of
infinitely many linear dependence congruences modulo 4 among suitable sets of
holomorphic normalizations of these series. These results rely on the theory
of class numbers and Hilbert class polynomials, and generalized twisted
Borcherds products developed by Bruinier and the author.

###### Key words and phrases:

Borcherds products, Heegner divisors, harmonic Maass forms, partition function

###### 2020 Mathematics Subject Classification:

(Primary) 11P83; (Secondary) 05A17

K.O. thanks the Thomas Jefferson Fund and the NSF (DMS-2002265 and
DMS-2055118) for their support.

##  1\. Introduction and statement of results

A partition of a non-negative integer $n$ is a non-increasing sequence of
positive integers that sums to $n.$ The partition function $p(n)$ counts the
number of partitions of $n.$ Ramanujan’s celebrated congruences [39, 40]
assert, for every non-negative integer $n,$ that

$\begin{split}p(5n+4)&\equiv 0\pmod{5},\\\ p(7n+5)&\equiv 0\pmod{7},\\\
p(11n+6)&\equiv 0\pmod{11}.\end{split}$

These congruences have inspired many works (for example, see [2, 3, 4, 5, 6,
7, 13, 20, 24, 26, 27, 34, 43] to name a few). Atkin [5] and Watson [43]
notably proved infinite families of congruences, where the moduli are
arbitrary powers of 5, 7 and 11.

In the 1960s, Atkin [6] discovered further congruences, such as

$p(17303n+237)\equiv 0\pmod{13},$

where the moduli involve primes other than 5, 7, and 11. In 2000, the author
revisited these examples, and he employed the Deligne-Serre theory of modular
$\ell$-adic Galois representations and Shimura’s theory of half-integral
weight modular forms [34] to prove the existence of infinitely many such
congruences for every prime modulus $\ell\geq 5$. Ahlgren and the author [2,
3] later extended this result to include all moduli coprime to 6.

In contrast to this rich theory of congruences, much less is known about
$p(n)$ modulo 2 and 3. Apart from a theorem of Radu [38], which proves that
there are no congruences of the form

$p(an+b)\equiv 0\pmod{\ell},$

where $\ell\in\\{2,3\\},$ very little is known. The parity of $p(n)$ is
conjectured to be “random”. Indeed, a well-known conjecture of Parkin and
Shanks [37] asserts that

$\lim_{X\rightarrow+\infty}\frac{\\#\\{n\leq X\ :\ p(n)\ {\text{\rm is even
(resp. odd)}}\\}}{X}=\frac{1}{2}.$

Although there have been many works (for example, see [1, 10, 16, 18, 19, 21,
23, 25, 28, 29, 30, 31, 32, 33, 37, 38, 42], to name a few) on the parity of
$p(n)$, this conjecture seems far out of reach. Indeed, the strongest results
do not preclude the possibility that

$\lim_{X\rightarrow+\infty}\frac{\\#\\{n\leq X\ :\ p(n)\ {\text{\rm is even
(resp. odd)}}\\}}{X^{\frac{1}{2}+\varepsilon}}=0.$

Even less is known about $p(n)$ modulo 3. Indeed, it isn’t known that $p(n)$
is infinitely often a multiple of 3, or takes values in any fixed residue
class modulo 3 infinitely often.

In view of these difficult problems, in 2010 the author initiated [36]
research in a different direction. Instead of concentrating on the values of
the partition function in arithmetic progressions, he identified special
quadratic polynomial arguments. For each square-free integer $1<D\equiv
23\pmod{24},$ he constructed modular forms that capture the parity of the
values $p\big{(}\frac{Dm^{2}+1}{24}\big{)}.$ Using their modularity, he proved
general theorems on the infinitude of the number of even (resp. odd) numbers
among such values. We refine this work to the modulus 4.

Throughout, we suppose that $1<D\equiv 23\pmod{24}$ is square-free. We define
the formal twisted generating function111We note that $p(\alpha)\coloneqq 0$
for non-integral $\alpha.$

(1.1) $P(D;q)\coloneqq\sum_{m,n\geq
1}\chi_{-D}(n)\chi_{12}(m)p\left(\frac{Dm^{2}+1}{24}\right)q^{mn},$

where $\chi_{d}(\cdot)$ is the discriminant $d$ Kronecker character.

######  Remark.

Apart from $p(0)=1$, each partition number occurs in a single $P(D;q),$ as
$24n-1=Dm^{2}$ determines the square-free $1<D\equiv 23\mod 24$. Furthermore,
modulo 2 we note that

$P(D;q)\coloneqq\sum_{\begin{subarray}{c}m,n\geq 1\\\
\gcd(n,D)=1\end{subarray}}p\left(\frac{Dm^{2}+1}{24}\right)q^{mn}\pmod{2}.$

Our first result establishes that each $P(D;q)$ is congruent modulo 4 to a
specific weight 2 meromorphic modular form that arises as the logarithmic
derivative of a twisted Borcherds product. To make this precise, we require
the notion of a Heegner point. Let $\mathcal{Q}_{N,-D}$ be the set of positive
definite integral binary quadratic forms $Q(x,y)=ax^{2}+bxy+cy^{2},$ where
$N\mid a,$ with discriminant $b^{2}-4ac=-D.$ The congruence subgroup
$\Gamma_{0}(N)$ acts on $\mathcal{Q}_{N,-D}$ with finitely many orbits. If
$Q\in\mathcal{Q}_{N,-D},$ then its Heegner point is

(1.2)
$\tau_{Q}\coloneqq\frac{-b+\sqrt{-D}}{2a}\in\mathbb{H}\coloneqq\\{\mbox{Im}(\tau)>0\\}.$

Since $\Gamma_{0}(N)$ acts on $\mathcal{Q}_{N,-D}$ with finitely many orbits,
there are finitely many such points in any fundamental domain for the action
of $\Gamma_{0}(N)$ on $\mathbb{H}.$

Every $P(D;q)$ is congruent modulo $4$ to the Fourier expansion at $i\infty$
(i.e. $q\coloneqq e^{2\pi i\tau}$) of a weight 2 meromorphic modular form with
simple poles supported at $\Gamma_{0}(6)$ discriminant $-D$ Heegner points.

######  Theorem 1.1.

If $1<D\equiv 23\pmod{24}$ is square-free, then there is a weight 2
meromorphic modular form $\mathcal{L}_{D}(\tau)$ on $\Gamma_{0}(6)$ for which

$P(D;q)\equiv\mathcal{L}_{D}(\tau)\pmod{4}.$

The poles of $\mathcal{L}_{D}(\tau)$ are simple and are at $\Gamma_{0}(6)$
discriminant $-D$ Heegner points.

Theorem 1.1 implies the existence of infinitely many linear dependence
congruences modulo 4 among $q$-series assembled from the $P(D;q).$ Although
these congruence relations generally require many $q$-series, they are
ubiquitous. Therefore, the partition function is not truly random modulo 4, as
its values satisfy a complicated infinite web of congruences.

To make this precise, let $S$ be a finite set of square-free integers
$1<D\equiv 23\pmod{24},$ and let $h_{S}:=\max\left\\{h(-D)\ :\ D\in
S\right\\},$ the maximum of the class numbers $h(-D).$ For each $D\in S,$ we
define the holomorphic normalizations

(1.3) $\widehat{P}_{S}(D;q)\coloneqq P(D;q)\cdot\Delta(\tau)^{h_{S}}\cdot
H_{-D}(1/\Delta(\tau)).$

For fundamental discriminants $-D<0,$ we recall the Hilbert class polynomial

(1.4)
$H_{-D}(X)\coloneqq\prod_{Q\in\mathcal{Q}_{-D}/\operatorname{PSL}_{2}(\mathbb{Z})}(X-j(\tau_{Q}))\in\mathbb{Z}[X],$

where $\mathcal{Q}_{-D}$ denotes the discriminant $-D$ positive definite
integral binary quadratic forms, and

$j(\tau)\coloneqq\frac{\left(1+240\sum_{n=1}^{\infty}\sum_{d\mid
n}d^{3}q^{n}\right)^{3}}{q\prod_{n=1}^{\infty}(1-q^{n})^{24}}=\frac{E_{4}(\tau)^{3}}{\Delta(\tau)}=q^{-1}+744+196884q+\dots$

is a hauptmodul for the field of modular functions on $\text{\rm
SL}_{2}(\mathbb{Z}).$ The class number $h(-D)$ is the degree of $H_{-D}(X).$

######  Theorem 1.2.

Assuming the notation and hypotheses above, if $\\#S>12h_{S}+2,$ then the
$q$-series $\left\\{\widehat{P}_{S}(D;q)\ :\ D\in S\right\\}$ are linearly
dependent modulo $4.$

######  Remarks.

(1) Let $S_{t}=\\{D_{1},D_{2},\dots,D_{t}\\}$ be the set of the first $t$
square-free integers $1\leq D\equiv 23\pmod{24}.$ In earlier work, Griffin,
Tsai and the author (see Lemma 2.2 of [17]) proved for fundamental
discriminants $-D<0$ that

$h(-D)\leq\frac{\sqrt{D}(\log D+2)}{\pi}.$

Since a positive proportion of the integers $n\equiv 23\pmod{24}$ are square-
free, this implies the existence of a linear dependence relation modulo 4
among the

$\\{\widehat{P}_{S_{t}}(D_{1};q),\widehat{P}_{S_{t}}(D_{2};q),\dots,\widehat{P}_{S_{t}}(D_{t};q)\\}$

for every sufficiently large $t$. In fact, there is a linear dependence for
any subset with more than $12h_{S_{t}}+2$ many $q$-series.

(2) If there are $t>12m+2$ distinct square-free integers $D_{i}\equiv
23\pmod{24}$ for which $h(-D_{i})=m,$ then the proof of Theorem 1.2 shows that
the $q$-series

$\\{P(D_{1};q)H_{-D_{1}}(1/\Delta(\tau)),P(D_{2};q)H_{-D_{2}}(1/\Delta(\tau)),\dots,P(D_{t};q)H_{-D_{t}}(1/\Delta(\tau))\\}$

are linearly dependent modulo 4. This follows from the fact that each
$\widehat{P}_{S}(D_{i};q)$ has the same factor $\Delta(\tau)^{m}$ that can be
factored out.

This note is organized as follows. In Section 2 we recall some twisted
Borcherds products that arise from Ramanujan’s 3rd order mock theta functions.
In Section 3 we study the logarithmic derivatives of these products, and we
study their reductions modulo 4. Then we prove Theorem 1.1. Finally, we prove
Theorem 1.2 in Section 4.

## Acknowledgements

The author thanks Andreas Mono and Badri Pandey for useful conversations
related to this paper.

##  2\. Borcherds products for Ramanujan’s 3rd order mock theta functions

In their work on derivatives of modular $L$-functions, the author and Bruinier
[15] produced generalized twisted Borcherds Products arising from weight 1/2
harmonic Maass forms. These constructions are extensions and generalizations
of previous automorphic infinite products obtained by Borcherds [8, 9] and
later Zagier [44]. The general results (see Theorems 6.1 and 6.2 of [15]) give
modular forms with twisted Heegner divisor whose infinite product expansions
arise from weight 1/2 harmonic Maass forms (see [11, 14] for background on
harmonic Maass forms).

Ramanujan’s mock theta functions, which are special examples of holomorphic
parts of harmonic Maass forms of weight 1/2 (for example, see [12, 13, 35, 45,
46, 47]), can be used to construct such products. We recall one example which
involves a vector-valued form assembled from the third order mock theta
functions $f(q)$ and $\omega(q),$ which are defined by

(2.1)
$f(q):=1+\sum_{n=1}^{\infty}\frac{q^{n^{2}}}{(1+q)^{2}(1+q^{2})^{2}\cdots(1+q^{n})^{2}},$

and

(2.2)
$\begin{split}\omega(q):=\sum_{n=0}^{\infty}\frac{q^{2n^{2}+2n}}{(q;q^{2})_{n+1}^{2}}=\frac{1}{(1-q)^{2}}+\frac{q^{4}}{(1-q)^{2}(1-q^{3})^{2}}+\frac{q^{12}}{(1-q)^{2}(1-q^{3})^{2}(1-q^{5})^{2}}+\cdots.\end{split}$

It is important to note that $f(q)$ and $\omega(q)$ have integer coefficients.

Using these mock theta functions, we define the vector-valued function

$\displaystyle R(\tau)\coloneqq(R_{0}(\tau),R_{1}(\tau),\dots,R_{11}(\tau)),$

with components

(2.3) $R_{j}(\tau)\coloneqq\begin{cases}0\ \ \ \ \ &{\text{\rm if}}\
j=0,3,6,9,\\\ \chi_{-12}(j)q^{-1}f(q^{24})\ \ \ \ &{\text{\rm if}}\
j=1,5,7,11\\\ 2q^{8}\left(-\omega(q^{12})+\omega(-q^{12})\right)\ \ \ \
&{\text{\rm if}}\ j=2,\\\ -2q^{8}\left(\omega(q^{12})+\omega(-q^{12})\right)\
\ \ \ &{\text{\rm if}}\ j=4,\\\
2q^{8}\left(\omega(q^{12})+\omega(-q^{12})\right)\ \ \ \ &{\text{\rm if}}\
j=8,\\\ 2q^{8}\left(\omega(q^{12})-\omega(-q^{12})\right)\ \ \ \ &{\text{\rm
if}}\ j=10.\end{cases}$

For convenience, we denote the coefficients of the components by

$R_{j}(\tau)\eqqcolon\sum_{n\geq n_{j}}C_{R}(j;n)q^{n}$

For each square-free $1<D\equiv 23\pmod{24}$, we have the rational function

(2.4) $P_{D}(X):=\prod_{b\\!\\!\\!\\!\mod
D}(1-e(-b/D)X)^{\left(\frac{-D}{b}\right)},$

where $e(\alpha):=e^{2\pi i\alpha}$ and $\left(\frac{-D}{b}\right)$ is the
Kronecker character for the negative fundamental discriminant $-D$. We define
the generalized Borcherds product $\Psi_{D}(\tau)$ by

(2.5)
$\Psi_{D}(\tau):=\prod_{m=1}^{\infty}P_{D}(q^{m})^{C_{R}(\overline{m};Dm^{2})}.$

Here $\overline{n}$ denotes the canonical residue class of $n$ modulo 12.

######  Theorem 2.1.

[§8.2 of [15]] Assuming the notation and hypotheses above, the function
$\Psi_{D}(z)$ is a weight 0 meromorphic modular form on $\Gamma_{0}(6)$ with a
discriminant -$D$ twisted Heegner divisor (see §5 of [15] for the explicit
divisor).

##  3\. Proof of Theorem 1.1

We recall the holomorphic differential operator

$\displaystyle\mathbb{D}\coloneqq q\frac{d}{dq}=\frac{1}{2\pi
i}\cdot\frac{\mathrm{d}}{\mathrm{d}\tau}.$

Let $\Psi_{D}(\tau)$ be the generalized Borcherds product in (2.5). Then, we
define the function

$\displaystyle\mathcal{L}_{D}(\tau)\coloneqq\frac{1}{-\sqrt{-D}}\frac{\left(\mathbb{D}\Psi_{D}\right)(\tau)}{\Psi_{D}(\tau)}.$

We prove the following properties about $\mathcal{L}_{D}(\tau)$.

######  Lemma 3.1.

Assuming the notation and hypotheses of Theorem 2.1, the following are true.
(1) The function $\mathcal{L}_{D}(\tau)$ is a weight 2 meromorphic modular
form on $\Gamma_{0}(6)$ with $q$-expansion

$\displaystyle\mathcal{L}_{D}(\tau)=\sum_{m=1}^{\infty}C_{R}(\overline{m};Dm^{2})m\sum_{\begin{subarray}{c}n=1\\\
\gcd(n,D)=1\end{subarray}}^{\infty}\left(\frac{-D}{n}\right)q^{mn}.$

(2) The poles of $\mathcal{L}_{D}(\tau)$ in the fundamental domain of
$\Gamma_{0}(6)$ are simple and are located at $\Gamma_{0}(6)$ Heegner points
with discriminant $-D.$

###### Proof.

Since $\mathbb{D}$ agrees with the Ramanujan-Serre derivative operator in
weight $0$, the function $\left(\mathbb{D}\Psi_{D}\right)(g;\tau)$ has weight
$2$. Together with Theorem 2.1, this proves the first assertion of part (1).
To verify the claimed Fourier expansion, we calculate that

$\displaystyle\frac{\left(\mathbb{D}P_{D}\right)(q^{m})}{P_{D}(q^{m})}$
$\displaystyle=q\frac{d}{dq}\sum_{b\text{ mod
}D}\left(\frac{-D}{b}\right)\log\left(1-e^{-2\pi
i\frac{b}{D}}q^{m}\right)=-m\sum_{b\text{ mod
}D}\left(\frac{-D}{b}\right)\sum_{n=1}^{\infty}\left(e^{-2\pi
i\frac{b}{D}}q^{m}\right)^{n}.$

By [22, (3.12)], we have

$\displaystyle\left(\frac{-D}{n}\right)\sum_{b\text{ mod
}D}\left(\frac{-D}{b}\right)e^{-\frac{2\pi ib}{D}}=\sum_{b\text{ mod
}D}\left(\frac{-D}{b}\right)e^{-\frac{2\pi ibn}{D}}\qquad\text{if
}\gcd(n,D)=1,$

and both sides of this identity vanish if $\gcd(n,D)>1$. We use [22, Theorem
3.3] to evaluate the left hand side of this identity assuming $\gcd(n,D)=1$,
and note that $-D\equiv 1\pmod{4}$. Therefore, we find that

$\displaystyle\mathcal{L}_{D}(\tau)$
$\displaystyle=\frac{1}{-\sqrt{-D}}\sum_{m=1}^{\infty}C_{R}(\overline{m};Dm^{2})\cdot
q\frac{d}{dq}\log\left(P_{D}(q^{m})\right)$
$\displaystyle=\sum_{m=1}^{\infty}C_{R}(\overline{m};Dm^{2})m\sum_{\begin{subarray}{c}n=1\\\
\gcd(n,D)=1\end{subarray}}^{\infty}\left(\frac{-D}{n}\right)q^{mn},$

which proves the claimed Fourier expansion of $\mathcal{L}_{D}(\tau)$. The
second claim follows from Theorem 2.1 as every zero or pole of
$\Psi_{D}(\tau)$ becomes a simple pole of its logarithmic derivative. ∎

Euler observed that

$P(q)\coloneqq\sum_{n=0}^{\infty}p(n)q^{n}=\prod_{n=1}^{\infty}\frac{1}{1-q^{n}}=\frac{q^{1/24}}{\eta(\tau)}.$

We have the following alternate identity for $P(q)$, which in turn provides a
crucial congruence between $P(q)$ and Ramanujan’s third order mock theta
function $f(q).$

######  Lemma 3.2.

As formal power series, we have that $P(q)\equiv f(q)\pmod{4}$.

###### Proof.

For every positive integer $m$, the $q$-series

$\frac{1}{(1-q)(1-q^{2})\cdots(1-q^{m})}=\sum_{n=0}^{\infty}a_{m}(n)q^{n}$

is the generating function for $a_{m}(n)$, the number of partitions of $n$
whose summands do not exceed $m$. Therefore, the $q$-series

$\frac{q^{m^{2}}}{(1-q)^{2}(1-q^{2})^{2}\cdots(1-q^{m})^{2}}=\sum_{n=0}^{\infty}b_{m}(n)q^{n}$

is the generating function for $b_{m}(n)$, the number of partitions of $n$
with a Durfee square of size $m^{2}$. The identity

$P(q)=1+\sum_{m=1}^{\infty}\frac{q^{m^{2}}}{(1-q)^{2}(1-q^{2})^{2}\cdots(1-q^{m})^{2}}$

follows by summing in $m$. The claimed congruence follows trivially from the
definition of $f(q)$ as

$\frac{q^{m^{2}}}{(1-q)^{2}(1-q^{2})^{2}\cdots(1-q^{m})^{2}}\equiv\frac{q^{m^{2}}}{(1+q)^{2}(1+q^{2})^{2}\cdots(1+q^{m})^{2}}\pmod{4}.$

∎

###### Proof of Theorem 1.1.

We note that $2(\omega(\pm q^{12})+\omega(\mp q^{12})\equiv 0\pmod{4}$, which
then implies that

$R_{j}(\tau)=\sum_{n\geq
n_{j}}C_{R}(j;n)q^{n}\equiv\begin{cases}\chi_{-12}(j)q^{-1}f(q^{24})\equiv\chi_{-12}(j)q^{-1}P(q^{24})\pmod{4}\
\ \ &{\text{\rm if}}\ j=1,5,7,11,\\\ 0\pmod{4}\ \ \ \ \ &{\text{\rm
otherwise.}}\end{cases}$

The theorem follows as the definition of $P(D;q)$ is congruent to
$\mathcal{L}_{D}(\tau)$ modulo 4 thanks to Lemmas 3.1 and 3.2, and the fact
that $\chi_{-12}(m)m\equiv\chi_{12}(m)\pmod{4}.$ ∎

##  4\. Proof of Theorem 1.2

Here we prove Theorem 1.2, which follows easily from a well-known theorem of
Sturm [41] on congruences between holomorphic modular forms.

###### Proof of Theorem 1.2.

For each square-free $1\leq D\equiv 23\pmod{24},$ we have that
$\mathcal{L}_{D}(\tau)$ is a weight 2 meromorphic modular form on
$\Gamma_{0}(6)$ with simple poles at discriminant $-D$ Heegner points.
Therefore, in terms of the modular function $j(\tau)=q^{-1}+744+196884q+\dots$
and the weight 12 cusp form $\Delta(\tau)=q-24q^{2}+\dots,$ we find that

$\mathcal{L}_{D}(\tau)\cdot\Delta(\tau)^{h(-D)}\cdot H_{-D}(j(\tau))$

is a weight $12h(-D)+2$ holomorphic modular form on $\Gamma_{0}(6).$ Moreover,
since $j(\tau)=E_{4}(\tau)^{3}/\Delta(\tau)\equiv 1/\Delta(\tau)\pmod{4},$ we
have that

$\mathcal{L}_{D}(\tau)\cdot\Delta(\tau)^{h(-D)}\cdot
H_{-D}(j(\tau))\equiv\mathcal{L}_{D}(\tau)\cdot\Delta(\tau)^{h(-D)}\cdot
H_{-D}(1/\Delta(\tau))\pmod{4}.$

Given the set $S$, we then find that each $\widehat{P}_{S}(D;q)$ is congruent
to a weight $12h_{S}+2$ holomorphic modular form on $\Gamma_{0}(6).$ However,
by a theorem of Sturm, two weight $k$ holomorphic modular forms are congruent
if their Fourier expansions are congruent for more than the first

$[\text{\rm SL}_{2}(\mathbb{Z}):\Gamma_{0}(6)]\cdot(12h_{S}+2)/12=12h_{S}+2$

many terms. ∎

## References

  * [1] S. Ahlgren, _Distribution of parity of the partition function in arithmetic progressions_ , Indagationes Math. (N.S.) 10 (1999), pages 173-181.
  * [2] S. Ahlgren, _The partition function modulo composite integers $M$_, Math. Ann. 318 (2000), pages 793-803.
  * [3] S. Ahlgren and K. Ono, _Congruence properties for the partition function_ , Proc. Natl. Acad. Sci., USA 98 (2001), pages 12882-12884.
  * [4] G. E. Andrews, _The theory of partitions_ , Cambridge Univ. Press, Cambridge, 1984.
  * [5] A. O. L. Atkin, _Proof of a conjecture of Ramanujan_ , Glasgow Math. J. 8 (1967), pages 14-32.
  * [6] A. O. L. Atkin, _Multiplicative congruence properties and density problems for $p(n)$_, Proc. London Math. Soc. (3) 18 (1968), pages 563-576.
  * [7] A. O. L. Atkin and J. N. O’Brien, _Some properties of $p(n)$ and $c(n)$ modulo powers of 13_, Trans. Amer. Math. Soc. 126 (1967), pages 442-459.
  * [8] R. E. Borcherds, _Automorphic forms on $O_{s+2,2}(\mathbb{R})$ and infinite products_, Invent. Math. 120 (1995), pages 161-213.
  * [9] R. E. Borcherds, _Automorphic forms with singularities on Grassmannians_ , Invent. Math. 132 (1998), pages 491–562.
  * [10] M. Boylan and K. Ono, _Parity of the partition function in arithmetic progressions II_ , Bull. London Math. Soc. 33 (2001), pages 558-564.
  * [11] K. Bringmann, A. Folsom, K. Ono, and L. Rolen, _Harmonic Maass forms and mock modular forms: Theory and applications_ , Vol. 64, Amer. Math. Soc., Providence, 2017.
  * [12] K. Bringmann and K. Ono, _The $f(q)$ mock theta function conjecture and partition ranks_, Invent. Math. 165 (2006), 243-266.
  * [13] K. Bringmann and K. Ono, _Dyson’s ranks and Maass forms_ , Ann. of Math., 171 (2010), pages 419-449.
  * [14] J. H. Bruinier and J. Funke, _On two geometric theta lifts_ , Duke Math. J. 125 (2004), 45-90.
  * [15] J. H. Bruinier and K. Ono, _Heegner divisors, $L$-functions and harmonic weak Maass forms_, Ann. of Math. 172 (2010), 2135-2181.
  * [16] F. Garvan and D. Stanton, _Sieved partition functions and $q$-binomial coefficients_, Math. Comp. 55 (1990), pages 299-311.
  * [17] M. Griffin, K. Ono, and W.-L. Tsai, _Heights of points on elliptic curves over $\mathbb{Q}$_, Proc. Amer. Math. Soc. 149 (2021), 5093-5100.
  * [18] M. Hirschhorn, _On the residue mod 2 and 4 of $p(n)$_, Acta Arith. 38 (1980), pages 105-109.
  * [19] M. Hirschhorn, _On the parity of $p(n)$ II_, J. Comb. Th. (A) 62 (1993), pages 128-138.
  * [20] M. Hirschhorn and D.C. Hunt, _A simple proof of the Ramanujan conjecture for powers of $5$_, J. Reine Angew. Math. 326 (1981), pages 1–17.
  * [21] M. Hirschhorn and M. Subbarao, _On the parity of $p(n)$_, Acta Arith. 50 (1988), pages 355-356.
  * [22] H. Iwaniec and E. Kowalski, Analytic number theory, American Mathematical Society Colloquium Publications, 53, American Mathematical Society, Providence, RI, 2004.
  * [23] O. Kolberg, _Note on the parity of the partition function_ , Math. Scand. 7 (1959), pages 377-378.
  * [24] J. Lovejoy and K. Ono, _Extension of Ramanujan’s congruences for the partition function modulo powers of 5_ , J. Reine Ange. Math. 542 (2002), pages 123-132.
  * [25] L. Mirsky, _The distribution of values of the partition function in residue classes_ , J. Math. Anal. Appl. 93 (1983), pages 593–598.
  * [26] M. Newman, _Periodicity modulo $m$ and divisibility properties of the partition function_, Trans. Amer. Math. Soc. 97 (1960), pages 225-236.
  * [27] M. Newman, _Congruences for the coefficients of modular forms and some new congruences for the partition function_ , Canadian J. Math. 9 (1957), pages 549-552.
  * [28] J.-L. Nicolas, _On the parity of generalized partition functions. II_ , Period. Math. Hungar. 43 (2001), pages 177–189.
  * [29] J.-L. Nicolas, _Valeurs impaires de la fonction de partition $p(n)$_, Int. J. Number Theory 2 (2006), pages 469–487.
  * [30] J.-L. Nicolas, _Parité des valeurs de la fonction de partition $p(n)$ et anatomie des entiers (French. English summary) [Parity of the values of the partition function $p(n)$ and anatomy of the integers_, Anatomy of integers, pages 97–113, CRM Proc. Lecture Notes, 46, Amer. Math. Soc., Providence, RI, 2008.
  * [31] J.-L. Nicolas, I. Z. Ruzsa, and A. Sárközy, _On the parity of additive representation functions (appendix by J.-P. Serre)_ , J. Number Th. 73 (1998), pages 292-317.
  * [32] K. Ono, _Parity of the partition function in arithmetic progressions_ , J. Reine Angew. Math. 472 (1996), pages 1-15.
  * [33] K. Ono, _Odd values of the partition function_ , Discrete Math. 169 (1997) pages 263-268.
  * [34] K. Ono, _The partition function modulo $m$_, Ann. of Math. 151 (2000), pages 293-307.
  * [35] K. Ono, _Unearthing the visions of a master: harmonic Maass forms and number theory_ , Proceedings of the 2008 Harvard-MIT Current Developments in Mathematics Conference, Somerville, Ma.,2009, pages 347-454.
  * [36] K. Ono, _Parity of the partition function_ , Adv. Math. 225 (2010), 349-366.
  * [37] T. R. Parkin and D. Shanks, _On the distribution of parity in the partition function_ , Math. Comp. 21 (1967), pages 466-480.
  * [38] C.-S. Radu, _A proof of Subbarao’s Conjecture_ , J. Reine Angew. Math. Math. 672 (2012), 161-175.
  * [39] S. Ramanujan, _Congruence properties of partitions_ , Proc. London Math. Soc. (2) 19, pages 207-210.
  * [40] S. Ramanujan, _Ramanujan’s unpublished manuscript on the partition and tau-functions_ (Commentary by B. C. Berndt and K. Ono) The Andrews Festschrift (Ed. D. Foata and G.-N. Han), Springer-Verlag, Berlin, 2001, pp. 39-110.
  * [41] J. Sturm, _On the congruence of modular forms_ , Springer Lect. Notes in Math. 1240, Springer-Verlag, Berlin (1984), pages 275-280.
  * [42] M. Subbarao, _Some remarks on the partition function_ , Amer. Math. Monthly 73 (1966), pages 851-854.
  * [43] G. N. Watson, _Ramanujan’s vermutung über zerfällungsanzahlen_ , J. Reine Angew. Math. 179 (1938), pages 97-128.
  * [44] D. Zagier, _Traces of singular moduli_ , in Motives, polylogarithms and Hodge theory, Part I (Irvine, CA, 1998), 211–244, Int. Press Lect. Ser., 3, I, Int. Press, Somerville, MA.
  * [45] D. Zagier, _Ramanujan’s mock theta functions and their applications [d’aprés Zwegers and Bringmann-Ono]_ (2006), Séminaire Bourbaki, no. 986.
  * [46] S. P. Zwegers, _Mock $\vartheta$-functions and real analytic modular forms_, $q$-series with applications to combinatorics, number theory, and physics (Ed. B. C. Berndt and K. Ono), Contemp. Math. 291, Amer. Math. Soc., (2001), pages 269-277.
  * [47] S. Zwegers, _Mock theta functions_ , Ph.D. Thesis (Advisor: D. Zagier), Universiteit Utrecht, 2002.