\section{Introduction}

The non-invasivity of \gls{fmri} led it to dominate brain mapping research
since the early 1990s \citep{fmridomination}.
In the past three decades, tens of thousands of published studies acquired and
analysed \gls{fmri} signals, producing new understanding of the human brain and
the cognitive function of its different components.
Quickly, the idea of \emph{meta-analysing} this ever-growing amount of
neuroimaging studies flourished.
By gathering and synthesising the findings of a large corpora of neuroimaging
studies, can we derive new knowledge about the brain's mechanics?
Can we study consensus within the cognitive neuroscience community?
The lack of power of neuroimaging studies undermines the reproducibility of
their findings \citep{poldrack_scanning_2017, botvinik-nezer_variability_2020}.
Can we \emph{build} consensus by aggregating results from several underpowered
studies into more robust findings supported by past literature?
Neuroimaging studies traditionally report peak activation coordinates in a
standard stereotactic coordinate system.
This makes it possible to compare them from one study to another.
These coordinates result from statistical hypothesis testing and represent a
condensed synthesis of which regions of the brain are reported as activated by
a given study.
Directly meta-analysing whole-brain unthresholded statistical maps (i.e.\
\emph{image}-based meta-analysis) is known to yield statistically more powerful
results \citep{salimi-khorshidi_meta-analysis_2009}; and the field is moving in
that direction \citep{gorgolewski_neurovaultorg_2015}.
However, most meta-analyses have resorted to \gls{cbma}: the meta-analysis of
these peak activation coordinates.
In the past ten years, an ecosystem of \gls{cbma} databases and tools was
brought to life \citep[e.g.][]{brainmap, neurosynth}, becoming an integral part
of brain mapping research.
Automatic \gls{cbma} databases, like Neurosynth, extract both natural language
processing features and peak activation coordinates from neuroimaging studies.
These tools are used to derive activation patterns
\citep[e.g.][]{wager_fmri-based_2013, cole_global_2012} or reveal meaningful
cognitive processes through reverse inference
\citep[e.g.][]{smallwood_science_2015, seghier_angular_2013,
chang_decoding_2013, andrews-hanna_default_2014}.
With them, researchers can define more robust regions of interest supported by
past literature.

Nonetheless, currently available tools are limited in the complexity of
queries that they can express and solve on \gls{cbma} data.
Neurosynth \citep{neurosynth} presents brain map for single-term queries.
Although technically feasible, term-based \glspl{cq} lead to underpowered
meta-analyses due to the small number of studies matching the queries.
Methods for exploiting \gls{cbma} data have recently been proposed, but they
concern themselves with either developing new procedures for thresholding
statistical brain maps, or integrating spatial priors into probabilistic models
by correlating nearby voxel activations \citep{cbma}.
We look at neuroimaging meta-analysis from a different angle by improving upon
existing \gls{cbma} literature through the development of a \gls{dsl} that
leverages past research on probabilistic logic programming languages and
databases to \emph{formulate and solve more expressive \gls{cbma} queries}.
We believe that, with this approach, more could be wrung out of this type of
data.

Recently, NeuroQuery \citep{neuroquery} produced meta-analyses using
unstructured text-based queries.
By encoding the relationship between terms in a vocabulary using a regularised
linear model, NeuroQuery can produce brain maps for underrepresented terms (few
studies exactly match the term).
However, NeuroQuery's queries are distinct from and harder to interpret than
database queries, which have clear semantics.
Moreover, producing a brain map from studies related to some term $t_1$
\emph{and not related} to some other term $t_2$ is not possible because
NeuroQuery cannot express logic-based queries.
Finally, NeuroQuery is not a probabilistic model that can be plugged into a
richer hierarchical model combining meta-analyses with heterogeneous
modalities, such as neuroanatomical and ontological knowledge.

Since the 1970s, the computer science community has been working on extending
logic programming languages \citep{prolog, abiteboul_foundations_1995} with
probabilistic semantics to represent knowledge uncertainty inherent to
real-world data \citep[reviewed by][]{deraedt}.
A wide variety of efficient approaches to answering questions (queries) from
these programs were developed; alongside seminal theoretical understandings.
\Acrlongpl{dsl} are not new to the cognitive neuroscience community.
The White Matter Query Language \citep{wmql} was developed to help experts
formally describe white matter tracts in a near-to-English syntax.
To the best of our knowledge, applying these techniques to the formulation and
resolution of logic-based queries on probabilistic \gls{cbma} databases has yet
to be attempted.
This approach could make it possible to formulate elaborate hypotheses on the
brain's function and structure and test them against past cognitive
neuroscience literature.

Adopting a language-oriented programming approach, we use probabilistic logic
programming languages to formulate and solve logic-based queries on \gls{cbma}
databases.
This work fits into a broader project to design a \gls{dsl}, coined
\emph{NeuroLang}, for expressing and testing cognitive neuroscience hypotheses
that combine meta-analysis, neuroanatomical and ontological knowledge.
\emph{The work presented here focuses on the probabilistic semantics of
NeuroLang and its application to term-based \gls{cbma} queries.}

Contributions of this work are three-fold.
First, we investigate the feasibility and technicalities of applying
probabilistic logic programming to \gls{cbma}-based brain mapping.
We propose a way to encode a \gls{cbma} database as a probabilistic logic
program based on CP-Logic \citep{cplogic}, on which complex \gls{cbma} queries
can be solved.
We translate this program to an equivalent Bayesian network representation in
order to show that correct answers to probabilistic queries can be derived from
its factorised joint probability distribution.
Second, we explain how leveraging lifted query processing techniques
\citep{lifted_braz, dichotomytheorem} allows us to scale to the large size of
neuroimaging data at the voxel level.
Third, we propose a relaxed modeling of \gls{tfidf} features to better encode
the relationship between terms and studies and show that fewer samples are
needed to solve two-term \glspl{cq} than traditional approaches, on simulated
and real \gls{cbma} databases.

\section{Background}

\subsection{Term-based queries on \gls{cbma} databases}
\label{sec:term-based-queries}

An example of term-based query formulated in plain English is: ``for each
region of the brain, what is the probability that studies associated with both
terms \emph{insula} and \emph{speech} report its activation?''.
The result of term-based queries are used in forward inference to obtain a map
of the brain's activated regions reported by studies matching the query.

A \gls{cbma} database of $N$ studies with a fixed vocabulary of $M$ terms can
be represented as two matrices $\bm{X} \in \mathbb{R}^{N, M}$ and $\bm{Y} \in
\{0, 1\}^{N, K}$, where $\bm{X}_{ij}$ is a \gls{tfidf} feature measured for
term $j$ in study $i$ and $\bm{Y}_{ik} = 1$ if voxel $k$ is reported as
activated in study $i$.
In practice, $\bm{Y}$ is a sparse matrix because only a small proportion of
voxels are reported within a single study.

Forward inference brain maps are constructed from a probabilistic model where
binary random variables $A_k$ and $T_j$ respectively model the activation of
each voxel $v_k$ and the association of studies to each term $t_j$.
$\textbf{P}\left[A_k\middle|T_j\right]$ is the probability that voxel $k$
activates in studies conditioned on the studies being associated with term $j$
and $\textbf{P}\left[A_k\middle|T_\text{insula} \wedge T_\text{speech}\right]$
is the probability that voxel $k$ activates in studies conditioned on studies
being associated with both terms `insula' and `speech'.\footnote{We use
$\textbf{P}\left[A_k\middle|T_i, T_j\right]$ to denote $\textbf{P}\left[A_k =
1\middle|T_i = 1, T_j = 1\right]$.}

Neurosynth \citep{neurosynth} associates terms to studies by applying a
threshold $\tau$ to \gls{tfidf} features $\bm{X}$.
Forward inference maps are obtained by estimating, for each voxel $k$,
\begin{equation}
    \textbf{P}\left[A_k\middle|T_j\right]
    =
    \frac{%
        \sum_{i=1}^N
        \bm{Y}_{ik}
        \bm{1}[\bm{X}_{ij} > \tau]
    }{%
        \sum_{i=1}^N
        \bm{1}[\bm{X}_{ij} > \tau]
    }
    \label{eq:neurosynth-single-term}
\end{equation}
Solving a query with a $p$-term conjunction, $\varphi = T_1 \wedge \cdots
\wedge T_p$, is done by estimating, for each voxel $k$,
\begin{equation}
    \textbf{P}\left[A_k\middle|\varphi\right]
    =
    \frac{%
        \sum_{i=1}^N
        \bm{Y}_{ik}
        \bm{1}[\text{min}(\bm{X}_{i1}, \dots, \bm{X}_{ip}) > \tau]
    }{%
        \sum_{i=1}^N
        \bm{1}[\text{min}(\bm{X}_{i1}, \dots, \bm{X}_{ip}) > \tau]
    }
    \label{eq:neurosynth-two-term}
\end{equation}
As terms are added to this conjunction (and thus, complexity to the query), the
term $\bm{1}[\text{min}(\bm{X}_{i1}, \dots, \bm{X}_{ip}) > \tau]$ goes to
zero for an increasing number of studies.
Rapidly, obtaining a meaningful brain map becomes infeasible due to
statistically weak results.
A different model that relaxes the hard thresholding of \gls{tfidf} features is
proposed in \cref{sec:st}.
Note that, solving a \emph{disjunction} of two terms is done by replacing
$\text{min}$ with $\text{max}$, thereby requiring that only one of the
\gls{tfidf} features passes the threshold.
The more terms are added, the larger the number of studies that are included in
the estimation.
In that case, statistical power is thus not an issue.

\subsection{\Acrlong{plp}}

Before diving into how probabilistic logic programming can be used to encode
\gls{cbma} data, we give a brief introduction to those languages through the
example of CP-Logic.
We also define the syntactic restrictions of the subset of this language that
we use in our \gls{dsl}.

\subsubsection{CP-Logic}

We use CP-Logic \citep{cplogic} as an intermediate representation in the
compilation of our \gls{dsl}.
In CP-Logic, programs contain rules (also called \emph{CP-Events}) of the form
\begin{equation}
    (h_1 : p_1 \vee \cdots \vee h_n : p_n) \gets \varphi
\end{equation}
where $h_i$ are head predicates, $p_i$ are probabilities such that $\sum_i p_i
\leq 1$, and the implication rule's body (also called \emph{antecedent})
$\varphi$ is a first-order logic formula.
All variables occurring in the head (also called \emph{consequent}) of the rule
must also occur in $\varphi$.
Such rules are interpreted as `$\varphi$ being true causes one of the atoms
$h_i$ to be true'.
Which $h_i$ becomes true is drawn from the probability distribution defined by
probabilities $p_i$.
CP-Logic programs define a probability distribution over the set of
\emph{possible worlds} \citep{sato_statistical_1995} associated with possible
executions of the probabilistic program.

\subsubsection{Syntactic restrictions and probabilistic databases}

Only a \emph{subset} of CP-Logic's expressive syntax is necessary to encode a
\gls{cbma} database and formulate term-based queries on it.
In NeuroLang, two kinds of rules are allowed.
Deterministic rules $( h : 1 ) \gets \varphi$, where $\varphi$ is a conjunction
of predicates and $h$ is a single head predicate that is true with probability
$1$ whenever $\varphi$ is true.
Probabilistic rules whose body is $\top$ (always true).
If the head of the rule contains a single predicate, it is a probabilistic
fact.
If it contains more than one head predicate, it is a probabilistic choice.
Moreover, recursive rules such as $(\text{A}(x) : 0.3) \gets A(y) \wedge
\text{B}(x)$ are not permitted in the program.

With these syntactic restrictions, probabilistic rules define relations in a
probabilistic database.
If a rule has more than one head predicate, its tuples are mutually
exclusive and partition the space of possible worlds.
Queries with mutually exclusive predicates are rewritten to be compatible with
probabilistic tuple-independent databases.
Deterministic rules of the program define \emph{\glspl{ucq}} on these
relations.
A \gls{ucq} $Q(\bm{x})$ is defined by a disjunction $\text{CQ}_1(\bm{x}),
\dots, \text{CQ}_n(\bm{x})$, where $\text{CQ}_i(\bm{x})$ are \glspl{cq} which
conjunct logic literals.
One major theoretical result in the field of probabilistic databases is the
\emph{dichotomy theorem} \citep{dichotomytheorem}.
It classifies \glspl{ucq} based on their complexity: those that can be solved
in polynomial time and those that are \#P-hard, in the size of the database.
A set of rules analyses the syntax of a given \gls{ucq} $Q(\bm{x})$ to derive
an algebraic expression that solves $\textbf{P}[Q(\bm{x})]$: the probability of
$Q(\bm{x})$ being true over all possible groundings of the database (i.e.\
possible worlds).
This resolution strategy is called \emph{lifted query processing}.
Guarantees on the efficiency of query resolution is of particular interest in
the context of neuroimaging's high-dimensional space.

\section{Probabilistic \gls{cbma} databases}
\label{sec:pp}

We now describe how \gls{cbma} data and queries can be encoded as a CP-Logic
program.
We then show how this program can be translated to a Bayesian network.
We use its factorised joint probability distribution to analytically derive the
same solutions for term-based queries as in \cref{sec:term-based-queries}.
Finally, we describe our approach to solving queries on this program using
lifted query processing strategies.

\subsection{Encoding a \gls{cbma} database as a probabilistic logic program}

The program of \cref{fig:program} encodes a \gls{cbma} database.
The equiprobable choice on the SelectedStudy relation partitions the space of
possible worlds such that each one corresponds to a particular study.
VoxelReported and TermInStudy predicates encode matrices $\bm{Y}$ and $\bm{X}$.
\begin{figure*}[t]
    \centering
    \includestandalone[width=0.7\textwidth]{fig_program}
    \caption{%
        CP-Logic program encoding a probabilistic \gls{cbma} database.
        $\text{TermInStudy}(t, s)$ models the presence of term $t$ in study
        $s$.
        $\text{VoxelReported}(v, s)$ encodes whether voxel $v$ was reported in
        study $s$.
        The large $\text{SelectedStudy}$ equiprobable choice over studies makes
        each possible world correspond to a specific study.
        $\text{Activation}(v)$ and $\text{TermAssociation}(t)$ respectively
        model the activation of voxel $v$ and the association with term $t$.
        The SUCC query $\textbf{P} \left[ \text{Activation}(v) \right]$ gives
        the marginal probability of activation of voxels over all studies.
        The query $\textbf{P} \left[ \text{Activation}(t) \middle|
        \text{TermAssociation}(\text{insula}) \right]$ results in a forward
        inference map for the term \emph{insula}.
    }
    \label{fig:program}
\end{figure*}
We write the program such that solving the query $\textbf{P}\left[
    \text{Activation}(v) \middle| \varphi\right]$, where $\varphi$ conjuncts
and/or disjuncts $\text{TermAssociation}(t_j)$ atoms, produces the
probabilistic model of term-based \gls{cbma} queries described in
\cref{sec:term-based-queries}.
For instance, when defining
\begin{equation*}
    \varphi =
    \text{TermAssociation}(\text{insula})
    \wedge
    \text{TermAssociation}(\text{speech})
\end{equation*}
$\textbf{P}\left[\text{Activation}(v)\middle|\varphi\right]$ is
equivalent to the query $\textbf{P}\left[A_k\middle|T_\text{speech} \wedge
T_\text{insula}\right]$ described previously.
We show that in the next section.

\subsection{Equivalence between the program of \cref{fig:program} and the
\gls{cbma} approach of \cref{sec:term-based-queries}}

To justify the design of the program in \cref{fig:program}, we translate it to
an equivalent Bayesian network representation using the algorithm proposed by
\citet{cplogicbn}.
The resulting Bayesian network is depicted in \cref{fig:bn} using
plate-notation.
\begin{figure*}[t]
    \centering\includestandalone{fig_bn}
    \caption{%
        Plate-notation representation of the Bayesian network translated from
        the program described in \cref{sec:pp}.
        Each ground atom in the program (e.g. $\text{TermInStudy}(t_2,
        s_{21})$) becomes a binary random variable with a deterministic
        \gls{cpd}.
        Specific \textbf{AND} nodes encode the conjunctions in the antecedent
        of the rules of the program.
        Choice random variables $c^\text{SS}$, $c_{ji}^\text{TIS}$,
        $c_{ki}^\text{VR}$ represent probabilistic choices in the program.
    }
    \label{fig:bn}
\end{figure*}
To simplify the notation, we use $A_k$, $T_n$ and $T_m$ to denote random
variables $\text{Activation}(v_k)$, $\text{TermAssociation}(t_n)$, and
$\text{TermAssociation}(t_m)$.
From the joint probability distribution defined by the Bayesian network, it can
be derived that
\begin{align}
    &\textbf{P}[A_k, T_n, T_m] \\
    &=
    \sum_{i=1}^N
    \textbf{P}[c^\text{SS} = i]
    \textbf{P}[c_{ki}^\text{VR} = 1]
    \textbf{P}[c_{ni}^\text{TIS} = 1]
    \textbf{P}[c_{mi}^\text{TIS} = 1] \\
    &=
    \frac{1}{N}
    \sum_{i=1}^N
    \bm{Y}_{ik}
    \bm{1}[\bm{X}_{in} > \tau]
    \bm{1}[\bm{X}_{im} > \tau]
\end{align}
and, similarly, that
\begin{align}
    \textbf{P}[T_n, T_m]
    &=
    \sum_{i=1}^N
    \textbf{P}[c^\text{SS} = i]
    \textbf{P}[c_{ni}^\text{TIS} = 1]
    \textbf{P}[c_{mi}^\text{TIS} = 1] \\
    &=
    \frac{1}{N}
    \sum_{i=1}^N
    \bm{1}[\bm{X}_{in} > \tau]
    \bm{1}[\bm{X}_{im} > \tau]
\end{align}
From these two joint probability distributions, the solution of the conditional
query can be derived using that
$\textbf{P}[A_k|T_n, T_m] = \frac{\textbf{P}[A_k, T_n, T_m]}{\textbf{P}[T_n,
T_m]} $, which gives the formula of \cref{eq:neurosynth-two-term}, for $p = 2$.

The same can be shown for disjunctive queries $\textbf{P}\left[A_k\middle|T_n
\vee T_m\right]$ by summing the results of 3 two-term \glspl{cq} as follows
\begin{equation}
    \begin{split}
        \textbf{P}\left[A_k\middle|T_n \vee T_m\right]
        &=
        \textbf{P}\left[A_k\middle|T_n, T_m\right] + \\
        & \textbf{P}\left[A_k\middle|\neg T_n, T_m\right] +
        \textbf{P}\left[A_k\middle|T_n, \neg T_m\right]
    \end{split}
\end{equation}

This confirms that the probabilistic program of \cref{fig:program} is sound, as
solving queries on the program leads to the statistical estimation described in
the previous section.

\subsection{Solving queries on probabilistic \gls{cbma} databases}

We now explore query resolution techniques that scale to the size of large
probabilistic \gls{cbma} databases.
The estimation of a forward inference brain map for a two-term conjunction
corresponds to the query
\begin{equation*}
    \footnotesize
    \textbf{P}\left[\text{Activation}(v)\middle|\text{TermAssociation}(t_i),
    \text{TermAssociation}(t_j)\right]
\end{equation*}
We solve this task by defining two \glspl{cq}
\begin{align*}
    Q_1(v) & \gets
    \begin{aligned}
        & \text{Activation}(v), \text{TermAssociation}(t_i), \\
        & \text{TermAssociation}(t_j)
    \end{aligned} \\
    Q_2 & \gets
        \text{TermAssociation}(t_i), \text{TermAssociation}(t_j)
\end{align*}
such that $\frac{\textbf{P}[Q_1(v)]}{\textbf{P}[Q_2]}$ solves the initial
query.
The numerator corresponds the joint probability of a voxel activation and the
association to the two terms. The denominator corresponds to the joint
probability of the associations to the two terms.

\subsubsection{%
    \Gls{kc} approaches do not scale to the size of neuroimaging data
}

We implemented the program of \cref{fig:program} in ProbLog2 \citep{problog2}.
We observed that, when solving two-term \glspl{cq}, grounding and compiling the
program to \glspl{sdd} was impractical.
Solving a two-term \gls{cq} takes more than 30 minutes on a recent laptop.
This is due to the large number of voxels, terms and studies modeled in the
program, leading to a large number of ground literals.
To give perspective on the scale of \gls{cbma} and neuroimaging data, a brain
is typically partitioned into a grid of about \num[group-separator={,}]{230000}
$2\text{mm}^3$ voxels.
On average, studies in the Neurosynth database report
$\num[group-separator={,}]{3165}$ voxel activatons.
There are $\num[group-separator={,}]{14371}$ studies and
$\num[group-separator={,}]{3228}$ terms in the Neurosynth database.
We also tried compiling our program manually to \glspl{sdd} \citep{sdd}.
Despite our efforts, which did note include exploring recent tree-building
strategies~\citep{amarilli:hal-01336514}, the resolution of queries was still
too slow to be practical for real-world applications.
Currently available \gls{cbma} tools are capable of solving single-term queries
in seconds.
Resolution of more complex queries should have a similar time complexity.

\subsubsection{Lifted processing of \glspl{ucq} on probabilistic \gls{cbma}
databases}
\label{sec:lifted}

We leverage theoretical results which have identified classes of queries that
lifted inference can solve in polynomial time.
The dichotomy theorem \citep[proven in][]{dichotomytheorem} provides a
procedure for checking that \glspl{ucq} are \emph{liftable}.
This theorem is convenient because it guarantees that any query such that the
lifted processing rules apply is guaranteed to be solvable in PTIME.
If the query is not liftable, we resort to \gls{kc}-based resolution
techniques.
Because the language does not have probabilistic clauses and prevents
recursivity, we can use its deterministic rules to construct \glspl{ucq}
associated with a given probabilistic query $\textbf{P}[\psi(\bm{x})]$, where
$\psi(\bm{x})$ is a conjunction of intensional, extensional or probabilistic
literals.
This lifted approach makes it possible to solve \gls{cq} in a few seconds.
\emph{Extensional query plans} \citep[see 4.1 of][]{queryprocessing} are
obtained and evaluated to solve queries using a custom Python relational
algebra engine.

\section{Relating terms and studies probabilistically}
\label{sec:st}

The hard thresholding $\bm{1}[x > \tau]$ of \gls{tfidf} features $x$ presented
in \cref{sec:term-based-queries} misses studies that could be relevant to the
resolution of queries.
Because we are interested in solving more complex queries, in this section we
explore a relaxation by introducing the soft-thresholding function
\begin{equation}
    \omega(x ; \alpha, \tau)
    \coloneqq
    \sigma \left( \alpha ( x - \tau ) \right)
    \quad
    \in [0, 1]
    \label{eq:omega}
\end{equation}
where $\sigma$ is the logistic function and $\tau$ a threshold.
As $\alpha$ increases, $\omega( x ; \alpha, \tau )$ converges towards the
hard-thresholding function $\bm{1}[x > \tau]$.
With an appropriate $\alpha$, a larger proportion of studies is included in the
calculation of $\textbf{P}\left[A_k\middle|\varphi\right]$, giving better
estimates on small databases.
For example, results of two-term \gls{cq} $\textbf{P}\left[A_k\middle|T_1
\wedge T_2\right]$ and \gls{ucq} $\textbf{P}\left[A_k\middle|T_1 \vee
T_2\right]$ queries are estimated with
\begin{align}
    \textbf{P}\left[A_k\middle|T_1 \wedge T_2\right]
    &=
    \frac{%
        \sum\limits_{i = 1}^N
        \bm{Y}_{ik}
        \omega\left(
            \bm{X}_{i1} ; \alpha, \tau
        \right)
        \omega\left(
            \bm{X}_{i2} ; \alpha, \tau
        \right)
    }{%
        \sum\limits_{i = 1}^N
        \omega\left(
            \bm{X}_{i1} ; \alpha, \tau
        \right)
        \omega\left(
            \bm{X}_{i2} ; \alpha, \tau
        \right)
    }\label{eq:conjunction-st}
    \\
    \textbf{P}\left[A_k\middle|T_1 \vee T_2\right]
    &=
    \frac{%
        \sum\limits_{i = 1}^N
        \bm{Y}_{ik}
        \big(
            1 -
            \prod\limits_{j=1}^2
            \left( 1 -
                \omega\left(
                    \bm{X}_{ij} ; \alpha, \tau
                \right)
            \right)
        \big)
    }{%
        \sum\limits_{i = 1}^N
        \big(
            1 -
            \prod\limits_{j=1}^2
            \left( 1 -
                \omega\left(
                    \bm{X}_{ij} ; \alpha, \tau
                \right)
            \right)
        \big)
    }\label{eq:disjunction-st}
\end{align}
More generally, $\textbf{P}\left[A_k\middle|\varphi\right]$ can be estimated
for first-order logic formulas $\varphi$ that blend conjunctions and
disjunctions of Boolean random variables $T_j, j \in 1, \dots, M$.
For example, if $\varphi = (T_1 \vee T_2) \wedge (T_3 \vee T_4)$, we have
\begin{equation}
    \textbf{P}\left[A_k\middle|\varphi\right]
    =
    \frac{%
        \sum\limits_{i = 1}^N
        \bm{Y}_{ik}
        f(\bm{X}_{i1}, \bm{X}_{i2})
        f(\bm{X}_{i3}, \bm{X}_{i4})
    }{%
        \sum\limits_{i = 1}^N
        f(\bm{X}_{i1}, \bm{X}_{i2})
        f(\bm{X}_{i3}, \bm{X}_{i4})
    }
\end{equation}
where
$f(x_1, x_2)
=
1 -
(1 - \omega(x_1 ; \alpha, \tau))
(1 - \omega(x_2 ; \alpha, \tau))
$.

This modeling is implemented simply by integrating $\omega(\bm{X}_{ij} ;
\alpha, \tau)$ as the probabilities of probabilistic facts
$\text{TermInStudy}(t_j, s_i)$ in the program of \cref{fig:program}.

\section{Experiments and results}

We compare our method with Neurosynth's on simulated \gls{cbma} databases
sampled from a generative model and on the Neurosynth database.
Using both models, we solve 55 different two-term \glspl{cq}
$\textbf{P}\left[A_k\middle|T_i \wedge T_j\right]$.

\subsubsection{%
    Gain of statistical power when solving two-term \glspl{cq}
    $\textbf{P}\left[A_k\middle|T_i \wedge T_j\right]$ on smaller simulated
    \gls{cbma} databases
}
\label{sec:simulations}

We evaluate our method on simulated small \gls{cbma} databases obtained by
sampling from the generative model of \cref{fig:generative-model}.
This generative model provides the ground truth of which voxels activate in
studies matching a given query of interest.
This binary classification setting makes it possible to compare models by
measuring their ability to identify true voxel activations for multiple sample
sizes.
We experimented with multiple numbers of voxels ($K \in [100, 1000]$).
Preliminary results showed that varying the number of voxels in this range does
not alter the results.
We report results for $K = \num[group-separator={,}]{1000}$ voxels, of which
5\% are activated in studies matching the query.
Predicted voxel activations are obtained by thresholding $p$-values computed
from each model's estimation of $\textbf{P}\left[A_k\middle|\varphi\right]$
using a $G$-test of independence.
We use a $p$-value threshold of $0.01$ and a Bonferroni correction for multiple
comparisons.
\begin{figure}
    \centering
    \includestandalone{fig_generative_model}
    \caption{%
        Model for generating \gls{cbma} databases of size $N$.
        $Z_\text{TF}^{(i)}$ models term frequencies in study $i$ and
        follows a logistic-normal distribution.
        $Z_\text{IDF}$ computes inverse document frequencies from
        $\{Z_\text{TF}^{(i)}\}_{i \in N}$.
        $p_k$ is the probability of activation of voxel $v_k$.
        Vectors $\beta_k$ are obtained from a rejection sampling
        scheme that controls the proportion of voxels that activate when
        the query is verified.
        $Z_\text{IDF}$, $\bm{\mu}$ and $\bm{\Sigma}$ are estimated from
        4168 scrapped PubMed abstracts.
    }
    \label{fig:generative-model}
\end{figure}
Simulation results for two-term \glspl{cq} are presented in
\cref{fig:simulation}, where we compare our model's and Neurosynth's $F_1$
scores across 55 two-term \glspl{cq}.
These queries correspond to all two-term combinations out of 11 terms (depicted
on the $y$-axis of the bottom plot of \cref{fig:simulation}) associated with a
sufficiently large number of studies within the Neurosynth database to produce
meaningful forward inference map.
The $F_1$ score measures the performance of a binary classifier by combining
its precision and recall into a single metric.
We see the advantage of our approach over Neurosynth's for smaller generated
samples where activations related to the query can be identified more reliably
(higher $F_1$ scores).
Multiple values of $\alpha$ in the range $[100, 1000]$ were tried during our
experiments.
However when $\alpha$ is too small or too large, the model tends to include
either too many (and irrelevant) or too few studies in the estimation.
When $\alpha$ tends to $0$, $\omega$ becomes equivalent to Neurosynth's hard
thresholding.
We found that a sweet spot for $\alpha$ was around $300$ and report results for
that value.
Drawing the sigmoid curve for $\alpha = 300$ confirms that this transformation
of TFIDF features is adequate because it maintains Neurosynth's hard
thresholding's property of giving a $0$ or $1$ probability to the lowest and
highest TFIDF features (respectively).

The proposed approach did not show an advantage over Neurosynth for solving
two-term disjunctive queries.
This is expected, as such queries do not reduce the number of studies
incorporated in the estimation of $\textbf{P}\left[A_k\middle|T_i \vee
T_j\right]$, as explained in \cref{sec:term-based-queries}.

\begin{figure*}[t]
    \centering\includestandalone[width=0.8\textwidth]{fig_simulation_f1_score}
    \centering\includestandalone[width=0.8\textwidth]{fig_simulation_f1_mats}
    \caption{%
        Comparison of Neurosynth's and our method's $F_1$ scores across 55
        \glspl{cq} $\textbf{P}\left[A_k\middle|T_i \wedge T_j\right]$
        on simulated \gls{cbma} databases of varying sample sizes.
        For each sample size, 100 random sub-samples were used.
        \textbf{Above}, $F_1$ score distributions on all queries are compared
        across sample sizes.
        \textbf{Below}, $F_1$ score matrices (white is $0$, black is $1$) are
        compared across sample sizes.
        The upper triangular contains scores of our method and the lower
        triangular contains scores of Neurosynth.
        The threshold $\tau = 0.1$ is used in both models.
        The value $\alpha = 300$ was empirically chosen.
        Varying $\alpha$ near this value does not change the results
        noticeably.
        Sample sizes were taken on a logarithmic scale.
    }\label{fig:simulation}
\end{figure*}

\subsubsection{Gain of activation consistency on a real \gls{cbma} database}
\label{sec:real}

We evaluate our method on the Neurosynth \gls{cbma} database.
Because we don't have a ground truth of which voxels activate for a given
query, we resort to comparing models based on the \emph{consistency} of their
predicted activations over many random sub-samples of the Neurosynth database.

From predicted activation maps of $K$ voxels obtained from $M$ sub-samples of a
\gls{cbma} database, the \emph{consistency} for a two-term conjunctive formula
$\varphi$ is computed as
\begin{equation}
    C_\varphi
    \coloneqq
    \frac{1}{K}
    \sum_{k = 1}^K
    \left(
    1 - 2 \times \left|
    \frac{1}{M}
    \sum_{m = 1}^M
    \hat{y}_{mk}^\varphi
     - \frac{1}{2}
    \right|
    \right)
\end{equation}
where $\hat{y}_{mk}^\varphi = 1$ if voxel $k$ is predicted to be activated in
sub-sample $m$ when formula $\varphi$ is true.
The closer to one, the closer the average activation is to $0$ or $1$, which
indicates a higher consistency across sub-samples.
The closer to zero, the closer the average activation is to $0.5$ which
indicates that the predicted activations are highly variable across samples.

Results are reported in \cref{fig:real}, where the distribution of
consistencies across the same 55 \glspl{cq} as in
\cref{sec:simulations} are shown for multiple sample sizes.
For the largest sample sizes, consistency scores are closer to $1$ with our
method than with Neurosynth's.
For a sample size of 2395 (chosen on a logarithmic scale), the average
consistency of our method was 0.48 while Neurosynth's was 0.4 (+20\%) across
samples and queries. For a sample size of 3856, we notice a 10\% improvement.

\begin{figure}[t]
    \centering
    \includestandalone[width=0.45\textwidth]{fig_real_result}
    \caption{%
        Comparison of both models' distributions of voxel activation
        consistency across \num[group-separator={,}]{1000} sub-samples of the
        Neurosynth's database, for 55 two-term \glspl{cq} and for multiple
        sample sizes.
        As the sample size increases, our method finds more consistent
        activations than Neurosynth.
    }\label{fig:real}
\end{figure}

We did not experiment with larger sample sizes due to the computational cost of
running the experiment on many Neurosynth subsamples for all \glspl{cq}.
Also, we were mainly interested in whether our approach would be more
consistent for smaller sample sizes.
We observed that the consistency between Neurosynth and our approach was
similar when both models were estimated on the entire database.
This means that the proposed approach is more consistent on smaller sample
sizes but equivalently consistent on larger sample sizes.
The NeuroLang program implementing this model is available at
\url{https://github.com/NeuroLang/NeuroLang/tree/master/examples/plot_neurosynth_relaxed_tfidf.py}.

\section{Discussion}

This work fits into a broader approach to design a \acrlong{dsl} for expressing
logic-based cognitive neuroscience hypotheses that combine neuroimaging data,
neuroanatomical probabilistic maps, ontologies and meta-analysis databases to
produce fine-grained brain maps supported by past literature and heterogeneous
data.

The number of voxels ($K = \num[group-separator={,}]{1000}$) used in the
simulations of \cref{sec:simulations} is orders of magnitude lower than on the
typical whole-brain neuroimaging setting, where $K \simeq \mathcal{O}(10^5)$.
We chose to lower the dimension in the simulation setting for computational
practicality purposes.
However, we believe that maintaining the same proportion of reported
activations as in the Neurosynth database was enough to confirm our approach on
simulations before applying it to real data, as reported in \cref{sec:real}.

The flexible syntax of logic-based language allow to express all kinds of
queries.
$\textbf{P}\left[{\text{TermAssociation}(t)}\middle|{\varphi_\text{insula} \vee
\varphi_\text{fMRI}}\right]$ queries for terms that are most probably
associated with a given pattern of activation, where $\varphi_\text{insula}$ is
a conjunction of logic predicates $\text{Activation}(v_k)$ whose probabilities
come from a neuroanatomical probabilistic map of the insula, and where
$\varphi_\text{fMRI}$ is also a conjunction of predicates
$\text{Activation}(v_k)$ whose voxels $v_k$ are based on neuroimaging data
coming from a custom \gls{fmri} study.

The current version of NeuroLang is limited in what it can model, mainly due
to the syntactic limitations on programs and queries that we had to make in
order to use lifted processing strategies that scale to the size of \gls{cbma}
data.
In cognitive neuroscience, there is interest in using spatial priors to give
the incentive to nearby voxels to co-activate \citep{kong_spatial_2018}.
Spatial priors could be formulated as recursive probabilistic rules, such as
$\text{Activation}(v_1) : f(d) \gets \text{Activation}(v_2),
d = \texttt{euclidean}(v_1, v_2)$, where $d$ is the Euclidean distance measure
between two regions of the brain, and where $f$ maps $d$ to a proper
probability in $[0, 1]$.
The resolution of such queries remains a challenge both in terms of methodology
and tractability.
Future progress in the field of probabilistic programming languages could open
the door to other queries of interest to the cognitive neuroscience community.

\section{Conclusion}

This work is a step towards incorporating complex meta-analyses in brain
mapping models.
We encode a \gls{cbma} database in a probabilistic logic program on which
general logic-based queries can be solved.
Leveraging efficient query resolution strategies on probabilistic databases, we
are able to scale to the size of neuroimaging data.
We experimented with a new method for solving two-term \glspl{cq} using
\gls{tfidf} features more efficiently than the hard thresholding scheme used by
Neurosynth.
This is promising but further investigation should be conducted to know whether
this method extends to queries that conjunct more terms or queries that blend
conjunctions and disjunctions.
The proposed method requires the same computational power as Neurosynth.

\section*{Free software}

The source code of NeuroLang is freely available on GitHub%
\footnote{\url{https://github.com/NeuroLang/NeuroLang}}.

\section*{Acknowledgments}

This work was funded by the ERC-2017-STG NeuroLang grant.
We thank the AAAI-2021 reviewers and the senior program committee member for
the quality and rigor of their comments on this work.
We also thank the organisers of the conference for orchestrating the review of
such a large number of submissions.
We are grateful to some of our colleagues who shared insights on our work:
Antonia Machlouzarides-Shalit, Hicham Janati, and Louis Rouillard-Odera.


\section{Introduction}

Optimization problems have numerous applications among many fields such as imaging, data science or machine learning, to name a few. Many such problems in these areas are often formulated as
\begin{equation} \label{min}
    \hat{x} \in \arg\min_{x\in X} \sum_{i=1}^n f_i(A_ix) + g(x)
\end{equation}
where $f_i:Y_i\to\mathbb{R}\cup\{\infty\}$ and $g:X\to\mathbb{R}\cup\{\infty\}$ are convex functionals, and $A_i:X\to Y_i$ are linear operators between finite-dimensional real Hilbert spaces.

Examples of such problems are total variation regularized image reconstruction \cite{fercoqBianchi2019,ROF} such as image denoising \cite{chambolle2016introduction} or PET reconstruction \cite{ehrhardt2019faster}; regularized empirical risk minimization \cite{shalevZhang,zhangXiao} such as support vector machine (SVM) \cite{boserGuyonVapnik} or least  absolute  shrinkage  and  selection operator (LASSO) \cite{cevherBigData}; and optimization with large number of constraints \cite{fercoq_etal2019,patrascuNecoara2017}, among others.

While some classical approaches such as gradient descent are not applicable when the functionals $f_i$ or $g$ are not smooth \cite{chambolle2016introduction}, primal-dual methods are able to find solutions to (\ref{min}) without assuming differentiability. For convex, proper and lower-semicontinuous functionals $f_i,g$, a primal-dual formulation for (\ref{min}) reads
\begin{equation} \label{saddle}
    \hat{x},\hat{y} \in \arg\min_{x\in X}\max_{y\in Y} \sum_{i=1}^n \langle A_ix,y_i \rangle - f_i^*(y_i) + g(x)
\end{equation}
where $f^*$ is the \textit{convex conjugate} of $f$ and $Y=\Pi_{i=1}^n Y_i$. We refer to any solution $\hat{w}=(\hat{x},\hat{y})$ of (\ref{saddle}) as a \emph{saddle point}.

A well-known example of primal-dual methods that solve (\ref{saddle}) is the Primal-Dual Hybrid Gradient (PDHG) \cite{chambollePock,esserZhangChan,CremersChambollePock}, as presented by Chambolle \& Pock (2011).
It naturally breaks down the complexity of (\ref{min}) into separate optimization problems by doing separate updates for the primal and dual variables $x,y$, as shown in (\ref{det}).
PDHG is proven to converge to a solution of (\ref{saddle}), however its iterations become very costly for large-scale problems, e.g. when $n\gg1$ \cite{spdhg}.

More recently, Chambolle et al. proposed the Stochastic Primal-Dual Hybrid Gradient (SPDHG) \cite{spdhg} which reduces the per-iteration computational cost of PDHG by randomly sampling the dual variable: at each step, instead of the full dual variable $y$, only a random subset of its coordinates $y_i$ gets updated. This offers significantly better performance than the deterministic PDHG for large-scale problems \cite{spdhg}. Examples of similar random primal-dual algorithms are found in
\cite{fercoqBianchi2019,gao2019randomPD,latafat2019randomPD,zhangXiao}.

We are interested in the convergence of SPDHG. In \cite{spdhg}, it is shown that, for arbitrary convex functionals $f_i$ and $g$, SPDHG converges in the sense of \textit{Bregman distances}, which does not imply convergence in the norm.
In this paper we present a proof for the almost sure convergence of SPDHG for convex but not necessarily strongly convex functionals, using alternative arguments to the recently proposed proof by Alacaoglu et al. (\cite{alacaoglu}, Theorem 4.4). In contrast to~\cite{alacaoglu}, where they represent SPDHG by using projections of a single iterative operator, we present a more intuitive representation through a random sequence of operators (Lemma~\ref{fixpoints}).
A summary of the proof is laid out in Section~\ref{sec:ske} and the complete proof is detailed in Section~\ref{sec:conv}. A comparison with the proof of Alacaouglu et al. and other related work is discussed in Section~\ref{sec:rel}.

Finally, in Section~\ref{sec:numeric} we look into the application of parallel Magnetic Resonance Imaging (MRI) in order to compare the performance of SPDHG with that of the deterministic PDHG.

\section{Algorithm} In order to solve (\ref{saddle}), the deterministic PDHG with dual extrapolation~\cite{chambollePock} reads
\begin{equation} \label{det}
\begin{aligned}
    x^{k+1} & =\mathrm{prox}_{\tau g}(x^{k} - \tau A^T\bar{y}^{k}) \\
    y^{k+1} &  =
    \mathrm{prox}_{\sigma f^{*}}(y^{k}+\sigma Ax^{k+1})
\end{aligned}
\end{equation}
\vspace{0mm}
where $\bar{y}^k = 2y^k - y^{k-1}$ is an extrapolation on the previous iterates, $A$ and $f$ are given by $Ax=(A_1x,...,A_nx)$ and $f(y)=\sum_{i}f_i(y_i)$, and the \textit{proximity operator} of any functional $f$ is defined as $\mathrm{prox}_{\sigma f}(v) := \arg\min_{y\in Y}\frac{\|v-y\|^2}{2} + \sigma f(y)$.

SPDHG, in contrast, reduces the cost of iterations by only partially updating the dual variable $y=(y_i)_{i=1}^n$: at every iteration $k$, choose $j\in\{1,...,n\}$ at random with probability $p_i=\mathbb{P}(j=i)>0$, so that only the variable $y_{j}^{k+1}$ is updated, while the rest remain unchanged, i.e. $y_i^{k+1} = y_i^{k}$ for $i\neq j$.

\begin{algorithm}[SPDHG]
\label{alg_spdhg}
\begin{algorithmic}
\STATE{Choose $\tau,\sigma_i>0$ and $x^0\in X$. Set $y^0=\textbf{0}\in Y$ and $z^0=\bar{z}^0=\textbf{0}\in X.$ }
\STATE{For $k\geq0$ do}
\STATE{ \vspace{-6mm}
\begin{equation*}
\begin{aligned}
& \hspace{11mm} \text{select } j^{k} \in\{1,...\,,n\} \text{ at random} \hspace{64mm} \\
& \begin{aligned}
    \hspace{16mm} x^{k+1} & =\mathrm{prox}_{\tau g}(x^{k} - \tau\bar{z}^k) \\
    y_{i}^{k+1} &  = \left\{
    \begin{tabular}{ll}
    $\mathrm{prox}_{\sigma_i f_{i}^{*}}(y_{i}^{k}+\sigma_{i} A_{i}x^{k+1})$ \quad\quad &  if $i=j^{k}$ \\
    $y_{i}^{k}$ & else
    \end{tabular} \right. \\
    \delta^k  &= A^T_{j^{k}}(y_{j^{k}}^{k+1}-y_{j^{k}}^k) \\
           z^{k+1} &= z^k + \delta^k  \\
     \bar{z}^{k+1} &= z^{k+1} + {p^{-1}_{j^k}}\delta^k
\end{aligned}
\end{aligned}
\end{equation*} }
\end{algorithmic}
\end{algorithm}

\section{Main Result}
\label{sec:main}
We establish the almost sure convergence of SPDHG for any convex functionals, under the same step size conditions as in \cite{spdhg}:

\begin{assumption} \label{assu}
We assume the following to hold:
\begin{enumerate}
    \item The set of solutions to (\ref{saddle}) is nonempty.
    \item The functionals $g,f_i$ are convex, proper and lower-semicontinuous.
    \item The step sizes $\tau,\sigma_i>0$ satisfy
    \begin{equation} \label{assu_i}
        \tau\sigma_i\|A_i\|^2 < p_i \quad\text{for every } i.
    \end{equation}
\end{enumerate}
\end{assumption}

\begin{theorem}[\textbf{Convergence of SPDHG}]
\label{main}
Let $(w^k)_{k\in\mathbb{N}}=(x^k,y^k)_{k\in\mathbb{N}}$ be a random sequence generated by Algorithm \ref{alg_spdhg}. Under Assumption~\ref{assu},
the sequence $(w^k)_{k\in\mathbb{N}}$ converges almost surely to a solution of (\ref{saddle}).
\end{theorem}

\section{Sketch of the Proof} \label{sec:ske}

The following results lay out the proof of Theorem~\ref{main}. The complete proof is detailed in Section \ref{sec:conv}.
We use the notation $\|x\|^2_T=\langle Tx,x\rangle$, as well as the block diagonal operators $Q,S:Y\to Y$ given by $Q=\mathrm{diag}(p_1^{-1},...,p_n^{-1})$ and $S=\mathrm{diag}(\sigma_1,...,\sigma_n)$. The conditional expectation at time $k+1$ is denoted, for any functional $\varphi$, by
$ \mathbb{E}^{k+1}(\varphi(w^{k+1})) = \mathbb{E}(\varphi(w^{k+1}) | w^k) .$

The proof of Theorem \ref{main} uses the following important inequality from SPDHG, which is a consequence of (\cite{spdhg}, Lemma~4.4). This inequality is best summarized in (\cite{alacaoglu}, Lemma~4.1), which we have further simplified by using the fact that Bregman distances of convex functionals are nonnegative (\cite{spdhg}, Section~4).

\begin{lemma}[\cite{alacaoglu}, Lemma 4.1]
\label{ala}
Let $(w^k)_{k\in\mathbb{N}}$ be a random sequence generated by Algorithm~\ref{alg_spdhg} under Assumption~\ref{assu}. Then for every saddle point $\hat{w}$,
\begin{equation} \label{inequality0}
\begin{aligned}
    V^k(w^k-\hat{w}) \geq \mathbb{E}^{k+1}(V^{k+1}(w^{k+1}-\hat{w})) + V(x^{k+1}-x^k,y^k-y^{k-1})
\end{aligned}
\end{equation}
where
$V$ and $V^k$ are given by
$V(x,y) = \|x\|^{2}_{\tau^{-1}} + 2\langle QAx,y \rangle + \|y\|^{2}_{QS^{-1}}$ and
$$ V^k(x,y) = \|x\|^{2}_{\tau^{-1}} - 2\langle QAx,y^k-y^{k-1} \rangle + \|y^k-y^{k-1}\|^2_{QS^{-1}} + \|y\|^{2}_{QS^{-1}}. $$
\end{lemma}

The following result is the central argument of our proof. It makes use of inequality~(\ref{inequality0}) and a classical result from Robbins \& Siegmund (Lemma~\ref{lemma2.2}) to establish an important convergence result. Its proof is detailed in Section~\ref{sec:conv}.

\begin{proposition} \label{tozero1}
Let $(w^k)_{k\in\mathbb{N}}$ be a random sequence generated by Algorithm \ref{alg_spdhg} under Assumption \ref{assu} and let $\hat{w}$ be a saddle point. Then:
\begin{enumerate}[label=\roman*)]
    \item \label{i} The sequence $(w^k)_{k\in\mathbb{N}}$ is a.s.\@ bounded.
    \item \label{ii} The sequence $(V^k(w^k-\hat{w}))_{k\in\mathbb{N}}$ converges a.s.
    \item \label{iii} The sequence $(\|w^k-\hat{w}\|)_{k\in\mathbb{N}}$ converges a.s.
    \item \label{iv} If every cluster point of $(w^k)_{k\in\mathbb{N}}$ is a.s.\@ a saddle point, the sequence $(w^k)_{k\in\mathbb{N}}$ converges a.s.\@ to a saddle point.
\end{enumerate}
\end{proposition}

Lastly, we prove that every cluster point of $(w^k)_{k\in\mathbb{N}}$ is almost surely a solution to (\ref{saddle}). This is also explained in Section~\ref{sec:conv}.

\begin{proposition} \label{cluster_saddle}
Let $(w^k)_{k\in\mathbb{N}}$ be a random sequence generated by Algorithm~\ref{alg_spdhg} under Assumption~\ref{assu}. Then every cluster point of $(w^k)_{k\in\mathbb{N}}$ is almost surely a saddle point.
\end{proposition}

\begin{proof}[Proof of Theorem \ref{main}]
By Proposition \ref{cluster_saddle}, every cluster point of $(w^k)_{k\in\mathbb{N}}$ is almost surely a saddle point and, by Proposition~\ref{tozero1} \ref{iv}, the sequence $(w^k)_{k\in\mathbb{N}}$ converges almost surely to a saddle point.
\end{proof}

\section{Proof of Convergence} \label{sec:conv}

This section contains detailed proofs for our two main arguments, Propositions~\ref{tozero1} and \ref{cluster_saddle}. The proof of Proposition~\ref{tozero1} follows a similar strategy to that of Combettes \& Pesquet in (\cite{combettesPesquet}, Proposition 2.3), and we have divided it into three sections.

\subsection{Proof of Proposition \ref{tozero1} \ref{i}}
To show this first part, we borrow the following lemma from \cite{spdhg} (Lemma~4.2 with $c=1$, $v_i = \tau\sigma_i\|A_i\|^2$ and $\gamma^2=\max_i\frac{v_i}{p_i}$):

\begin{lemma}[\cite{spdhg}, Lemma 4.2]
\label{lemma4.2}
Let $p_i^{-1}\tau\sigma_i\|A_i\|^2 \leq \gamma^2<1$ for every $i$ and let $y^k$ be defined as in Algorithm \ref{alg_spdhg}. Then for every $x\in X$,
\begin{equation*}
    \mathbb{E}^k(V(x,y^k-y^{k-1})) \geq (1-\gamma)\mathbb{E}^k \big(\|x\|^{2}_{\tau^{-1}} + \|y^k-y^{k-1}\|^{2}_{QS^{-1}} \big) .
\end{equation*}
\end{lemma}

\begin{proof}[Proof of Proposition \ref{tozero1} \ref{i}]
By (\cite{alacaoglu}, Lemma~4.1), for any
saddle point $\hat{w}$ we have
\begin{equation} \label{inequality}
    \Delta^k \geq \mathbb{E}^{k+1}(\Delta^{k+1}) + V(x^{k+1}-x^k,y^k-y^{k-1})
\end{equation}
where $\Delta^k = V^k(w^k-\hat{w})$.
By Lemma \ref{lemma4.2},
\begin{equation*}
    \mathbb{E}^k(V(x^{k+1}-x^k,y^k-y^{k-1})) \geq  (1-\gamma)\mathbb{E}^k \Big\{\|x^{k+1}-x^k\|^{2}_{\tau^{-1}} + \|y^k-y^{k-1}\|^{2}_{QS^{-1}} \Big\}.
\end{equation*}
Hence, taking the full expectation in (\ref{inequality}) yields
\begin{equation*}
    \mathbb{E}(\Delta^k) \geq \mathbb{E}(\Delta^{k+1})
    + (1-\gamma)\mathbb{E}\big(\|x^{k+1}-x^{k}\|^{2}_{\tau^{-1}} + \|y^{k}-y^{k-1}\|^{2}_{QS^{-1}}\big) .
\end{equation*}
Taking the sum from $k=0$ to $k=N-1$ gives
\begin{equation} \label{sum}
    \Delta^0 \geq \mathbb{E}(\Delta^{N})
    + (1-\gamma)\mathbb{E}\Big\{\sum_{k=0}^{N-1} \|x^{k+1}-x^{k}\|^{2}_{\tau^{-1}} + \|y^{k}-y^{k-1}\|^{2}_{QS^{-1}}\Big\}
\end{equation}
where $y^{-1}=y^0$. This implies $\Delta^0 \geq \mathbb{E}(\Delta^{N})$ and, by Lemma~\ref{lemma4.2} we have
\begin{equation}\label{E(Delta^N)}
    \mathbb{E}(\Delta^N) \geq \mathbb{E}\big\{(1-\gamma)(\|x^N-\hat{x}\|^2_{\tau^{-1}} + \|y^N-y^{N-1}\|^2_{QS^{-1}}) + \|y^N-\hat{y}\|^2_{QS^{-1}}\big\}.
\end{equation}
It follows that $ \Delta^0 \geq (1-\gamma)\|x^N-\hat{x}\|^2_{\tau^{-1}} + \|y^N-\hat{y}\|^2_{QS^{-1}} $ a.s., from where it is clear that the sequence $(w^N)_{N\in\mathbb{N}}$ is bounded almost surely.
\end{proof}

\subsection{Proof of Proposition \ref{tozero1} \ref{ii}-\ref{iii}}
As in (\cite{combettesPesquet}, Proposition 2.3), we use a classical result from Robbins \& Siegmund:

\begin{lemma}[\cite{robbins1971}, Theorem 1]
\label{lemma2.2}
Let $\mathcal{F}_k$ be a sequence of sub-$\sigma$-algebras such that $\mathcal{F}_k\subset \mathcal{F}_{k+1}$ for every $k$, and let $\alpha_k$, $\eta_k$ be nonnegative $\mathcal{F}_k$-measurable random variables such that $\sum_{k=1}^\infty \eta_k <\infty$ almost surely and
$$ \mathbb{E}(\alpha_{k+1}\,|\,\mathcal{F}_k) \,\leq\, \alpha_k + \eta_k \;\text{ a.s.}$$
for every $k$. Then $\alpha_k$ converges almost surely to a random variable in $[0,\infty)$.
\end{lemma}

\begin{proof}[Proof of Proposition \ref{tozero1} \ref{ii}-\ref{iii}]
From (\ref{E(Delta^N)}) we have $\mathbb{E}(\Delta^N)\geq0$. Thus taking the limit as $N\to\infty$ in (\ref{sum}) yields
\begin{equation} \label{ntozero}
     \mathbb{E}\Big\{\sum_{k=0}^\infty \|x^{k+1}-x^{k}\|^{2}_{\tau^{-1}} + \|y^{k}-y^{k-1}\|^{2}_{QS^{-1}}\Big\} < \infty
\end{equation}
which implies
\begin{equation}\label{sum(y^n)}
    \sum_{k=0}^\infty \|x^{k+1}-x^{k}\|^{2}_{\tau^{-1}} + \|y^{k}-y^{k-1}\|^{2}_{QS^{-1}} < \infty \quad\text{a.s.}
\end{equation}
and, in particular,
\begin{equation}\label{yto0}
   \|y^k-y^{k-1}\|_{QS^{-1}}\to0 \quad\text{a.s.}
\end{equation}
Since $(w^k)_{k\in\mathbb{N}}$ is bounded a.s., so is $(x^k)_{k\in\mathbb{N}}$ and, since the operators $Q$, $A$ and $S$ are also bounded, there exists $M>0$ such that, for every $k$,
\begin{equation*} \label{Q}
    |\langle QA(x^k-\hat{x}),y^{k}-y^{k-1} \rangle| \leq \|QA\|\|x^k-\hat{x}\|\|y^k-y^{k-1}\| \leq M\|y^k-y^{k-1}\|_{QS^{-1}}
\end{equation*}
a.s. and therefore, by (\ref{yto0}),
\begin{equation} \label{qto0}
    \langle QA(x^k-\hat{x}),y^{k}-y^{k-1} \rangle \to 0 \quad \text{a.s.}
\end{equation}
The fact that $(w^k)_{k\in\mathbb{N}}$ is a.s.\@ bounded, together with (\ref{yto0}) and (\ref{qto0}) imply the sequence $(\Delta^k)_{k\in\mathbb{N}}$ is also a.s.\@ bounded. Thus there exists $\hat{M}\geq0$ such that $ \Delta^k+\hat{M} \geq 0 $ for every $k$. Let $\alpha_k = \Delta^k+\hat{M}$ and $\eta_k=2|\langle QA(x^{k+1}-x^{k}),y^{k}-y^{k-1} \rangle|$. From (\ref{inequality}) we deduce
\begin{equation} \label{supermartingale}
    \alpha_k + \eta_k \,\geq\, \mathbb{E}^{k+1}(\alpha_{k+1}) \quad\text{a.s. for every } k,
\end{equation}
where all the terms are nonnegative and,  for some $\tilde{M}>0$,
\begin{equation*}
\begin{aligned}
    \eta_k = 2|\langle QA(x^{k+1}-x^{k}),y^{k}-y^{k-1} \rangle|
    &\leq 2\|Q A\| \|x^{k+1}-x^k\| \|y^k-y^{k-1}\| \\
    &\leq 2\tilde{M}  \|x^{k+1}-x^k\|_{\tau^{-1}} \|y^k-y^{k-1}\|_{QS^{-1}} \\
    &\leq \tilde{M} \big( \|x^{k+1}-x^k\|_{\tau^{-1}}^2 + \|y^k-y^{k-1}\|^2_{QS^{-1}} \big)
\end{aligned}
\end{equation*}
which implies, by (\ref{sum(y^n)}), $\sum_{k=1}^\infty\eta_k<\infty$ a.s.. Thus (\ref{supermartingale}) satisfies all the assumptions of Lemma~\ref{lemma2.2} and it yields $\Delta^k \to \alpha$ a.s.
for some $\alpha\in[-\hat{M},\infty)$. Furthermore, from (\ref{yto0}) and (\ref{qto0}) we know some of the terms in $\Delta^k$ converge to 0 a.s., namely
$$ - 2\langle QA(x^{k}-\hat{x}),y^{k}-y^{k-1}\rangle
+ \|y^{k}-y^{k-1}\|^{2}_{QS^{-1}} \to 0 \quad\text{a.s.} $$
hence $ \|x^k-\hat{x}\|^2_{\tau^{-1}} + \|y^k-\hat{y}\|^2_{QS^{-1}} \to \alpha $ a.s.
Finally, the norm $\|w\|^2_R := \|x\|^2_{\tau^{-1}} + \|y\|^2_{QS^{-1}} $ is equivalent to the norm in $X\times Y$. Since the sequence $(\|w^k-\hat{w}\|_R)_{k\in\mathbb{N}}$ converges a.s., so does $(\|w^k-\hat{w}\|)_{k\in\mathbb{N}}$.
\end{proof}

\subsection{Proof of Proposition \ref{tozero1} \ref{iv}}
The two following lemmas are consequence of (\cite{combettesPesquet}, Proposition 2.3) and their proofs are not included here. We use the standard notation $(\mathbf{\Omega},\mathcal{F},P)$ for the probability space corresponding to the random iterations $w^k$.

\begin{lemma}[\cite{combettesPesquet}, Proposition 2.3 iii)]
\label{tozero2}
Let $\mathbf{F}$ be a closed subset of a separable Hilbert space and let $(w^k)_{k\in\mathbb{N}}$ be a sequence of random variables such that the sequence $(\|w^k-w\|)_{k\in\mathbb{N}}$ converges almost surely for every $w\in\mathbf{F}$.
Then there exists $\Omega\in\mathcal{F}$ such that $\mathbb{P}(\Omega)=1$ and the sequence $(\|w^k(\omega)-w\|)_{k\in\mathbb{N}}$ converges for all  $\omega\in\Omega$ and $w\in\mathbf{F}$.
\end{lemma}

\begin{lemma}[\cite{combettesPesquet}, Proposition 2.3 iv)]
\label{tozero3}
Let $\mathbf{G}(w^k)$ be the set of cluster points of a random sequence $(w^k)_{k\in\mathbb{N}}$. Assume there exists ${\Omega}\in\mathcal{F}$ such that $\mathbb{P}({\Omega})=1$ and for every $\omega\in{\Omega}$, $\mathbf{G}(w^k(\omega))$ is nonempty and the sequence $(\|w^k(\omega)-w\|)_{k\in\mathbb{N}}$ converges for all $w\in\mathbf{G}(w^k(\omega))$.
Then $(w^k)_{k\in\mathbb{N}}$ converges almost surely to an element of $\mathbf{G}(w^k)$.
\end{lemma}

\begin{proof}[Proof of Proposition \ref{tozero1} \ref{iv}] Let $\mathbf{F}$ be the set of solutions to the saddle point problem (\ref{saddle}).
By Proposition \ref{tozero1} iii) and Lemma~\ref{tozero2}, there exists $\Omega\in\mathcal{F}$ such that the sequence $(\|w^k(\omega)-w\|)_{k\in\mathbb{N}}$ converges for every $w\in\mathbf{F}$ and $\omega\in\Omega$. This implies, since $\mathbf{F}$ is nonempty, that $(w^k(\omega))_{k\in\mathbb{N}}$ is bounded and thus $\mathbf{G}(w^k(\omega))$ is nonempty for all $\omega\in\Omega$. By assumption, there exists $\tilde{\Omega}\in\mathcal{F}$ such that
$\textbf{G}(w^k(\omega))\subset\mathbf{F}$ for every $\omega\in\tilde{\Omega}$. Let $\omega\in\Omega\cap\tilde{\Omega}$, then  $(\|w^k(\omega)-w\|)_{k\in\mathbb{N}}$ converges for every $w\in\textbf{G}(w^k(\omega))\neq\emptyset$. By Lemma~\ref{tozero3}, we get the result.
\end{proof}

\subsection{Proof of Proposition \ref{cluster_saddle}}
The following lemma describes Algorithm~\ref{alg_spdhg} as a random sequence $(T_{j^k})_{k\in\mathbb{N}}$ of continuous operators on the primal-dual space $X\times Y$, such that the fixed points of the operators are saddle points of (\ref{saddle}):

\begin{lemma}
\label{fixpoints}
Denote $w=(w_i)_{i=0}^n=(x,y_1,...,y_n)$ and for every $j\in \{1,...,n\}$ let the operator $T_j:X\times Y \to X\times Y$ be defined by
\begin{equation*}
\begin{aligned}
\hspace{20mm}
(T_jw)_0 &= \mathrm{prox}_{\tau g}\big(x-\tau A^Ty -  \big(1+\frac{1}{p_j}\big)  \tau A_j^T((T_jw)_j-y_j)\big) \\
(T_jw)_i &= \left\{ \begin{tabular}{ll}
    $\mathrm{prox}_{\sigma_i f_i^*}(y_i+\sigma_iA_ix)$ \quad & if $i=j$ \\
    $y_i$ & else
\end{tabular} \right.
\hspace{9mm} \text{for } 1\leq i\leq n .
\end{aligned}
\end{equation*}
Then the iterations $w^k$ generated by Algorithm \ref{alg_spdhg} satisfy
\begin{equation} \label{itTs}
    T_{j^k}(x^{k+1},y^{k}) = (x^{k+2},y^{k+1}) .
\end{equation}
Furthermore, $\hat{w}$ is a solution to the saddle point problem (\ref{saddle}) if and only if it is a fixed point of $T_{j}$ for each $j\in \{1,...,n\}$.
\end{lemma}

\begin{proof}
By definition of the iterates in Algorithm \ref{alg_spdhg}, $(T_{j^k}(x^{k+1},y^{k}))_i = y_i^{k+1}$ for every $i\in\{1,...,n\}$. By induction it is easy to check that $z^k = A^Ty^k$, and thus
\begin{equation*}
\begin{aligned}
    \bar{z}^{k+1} &= z^k + (1+\frac{1}{p_{j^k}})\delta^k \\
    &= A^Ty^{k} + (1+\frac{1}{p_{j^k}}) A_{j^k}^T(y_{j^k}^{k+1}-y_{j^k}^{k}) \\
    &= A^Ty^{k} + (1+\frac{1}{p_{j^k}}) A_{j^k}^T ((T_{j^k}(x^{k+1},y^{k}))_{j^k}-y_{j^k}^{k}) .
\end{aligned}
\end{equation*}
Thus
$ (T_{j^k}(x^{k+1},y^{k}))_0 = \mathrm{prox}_{\tau g}(x^{k+1} - \tau \bar{z}^{k+1}) = x^{k+2} $
, which proves (\ref{itTs}). Now let $w$ be a fixed point of $T_j$ for every $j$. Then, for any $j$,
$$y_j = w_j= (T_j w)_j = \text{prox}_{\sigma_i f_j^*}(y_j+\sigma_j A_jx), $$
from where it follows that, for any $j$,
\begin{equation*}
\begin{aligned}
    x = w_0 = (T_j w)_0
    &= \text{prox}_{\tau g}(x-\tau A^Ty - (1+\frac{1}{p_j}) \tau A_j^T((T_j w)_j-y_j) \\
    &= \text{prox}_{\tau g}(x-\tau A^Ty).
\end{aligned}
\end{equation*}
These conditions on $x$ and $y$ define a saddle point (\cite{bredieslorenz}, 6.4.2). The converse result is direct.
\end{proof}

\begin{proof}[Proof of Proposition~\ref{cluster_saddle}]
Let $j^k$ be the sampling generated by the algorithm and let $z^k= (x^{k+1},y^{k})$. By Lemma \ref{fixpoints} we have $z^{k+1} = T_{j^k}z^k$ and, by~(\ref{ntozero}),
\begin{equation} \label{n-n+1}
    \mathbb{E}(\|z^{k} - z^{k-1}\|^2) =
    \mathbb{E}(\|x^{k+1} - x^k\|^2 + \|y^k-y^{k-1}\|^2)
    \to 0 .
\end{equation}
Furthermore, by the properties of the conditional expectation,
\begin{equation*}
\begin{aligned}
\mathbb{E}(\|z^{k+1} - z^k\|^2) = \mathbb{E}(\mathbb{E}^k(\|z^{k+1} - z^k\|^2))
&= \mathbb{E}\Big(\sum_{j=1}^n \mathbb{P}(j^k = j)\|T_j z^{k} - z^k\|^2\Big)\\
&= \sum_{j=1}^n \mathbb{P}(j^k = j)\mathbb{E}(\|T_j z^{k} - z^k\|^2) .
\end{aligned}
\end{equation*}
By assumption $p_j = \mathbb{P}(j^k=j) >0$, thus by (\ref{n-n+1}) we have
$\mathbb{E}(\|T_j z^k-z^k\|^2)\to0$ for every $j$
and therefore
\begin{equation} \label{Tstozero}
   T_j z^k-z^k \to 0 \quad\text{a.s. for every } j\in \{1,...,n\}.
\end{equation}
Assume now a convergent subsequence $w^{\ell_k}\to w^*$. From (\ref{n-n+1}), $ y^{k}-y^{k-1}\to0 $ a.s. and so $z^{\ell_k}$ also converges to $w^*$. By (\ref{Tstozero}) and the continuity of $T_j$, there holds
$$ w^* = \lim_{k\to\infty}z^{\ell_k} = \lim_{k\to\infty} T_{j}z^{\ell_k} = T_{j}\big(\lim_{k\to\infty}z^{\ell_k}\big) = T_jw^* \quad\text{a.s.}$$
for every $j$. Hence $w^*$ is almost surely a fixed point of $T_j$ for each $j$ and, by Lemma~\ref{fixpoints}, $w^*$ is a saddle point.
\end{proof}

\section{Relation to Other Work}\label{sec:rel}
\subsection{Chambolle et al. (2018)}
In the original paper for SPDHG \cite{spdhg}, it is shown that, under Assumption~\ref{assu}, the \textit{Bregman distance} to any solution $\hat{x},\hat{y}$ of (\ref{saddle}) converges to zero, i.e. the iterates $x^k,y^k$ of Algorithm \ref{alg_spdhg} satisfy
\begin{equation} \label{gap}
    D_g^{-A^T\hat{y}}(x^k,\hat{x}) + D_{f^*}^{A\hat{x}}(y^k,\hat{y}) \to 0 \quad\text{a.s.,}
\end{equation}
where the Bregman distance is defined by $ D_h^q(u,v) = h(u) - h(v) - \langle q,u-v \rangle $ for any functional $h$ and any point $q\in\partial h(v)$ in the subdifferential of $h$.
In \cite{spdhg}, it is shown that (\ref{gap}) implies $(x^k,y^k) \to (\hat{x},\hat{y})$ a.s.\@ if $f_i$ or $g$ are strongly convex.

\subsection{Combettes \& Pesquet (2014)}
In \cite{combettesPesquet}, Combettes \& Pesquet look into the convergence of random sequences $(w^k)_{k\in\mathbb{N}}$ of the form
\begin{equation} \label{T}
    w_i^{k+1} = \left\{ \begin{tabular}{ll}
        $(T w^k)_i$ & if $i\in\mathbb{S}^k$ \\
        $w_i^k$ & else
    \end{tabular} \right.
\end{equation}
where $\mathbb{S}^k\subset\{1,...,n\}$ is chosen at random. They use the Robbins-Siegmund lemma (Lemma~\ref{lemma2.2}) to prove that, for a nonexpansive operator $T$, the sequence $(w^k)_{k\in\mathbb{N}}$ converges to a fixed point of $T$. Later, Pesquet \& Repetti \cite{pesquetRepetti} used this to prove convergence for a wide class of random algorithms of the form (\ref{T}), where $T = (I+B)^{-1}$ is the resolvent operator of a monotone operator $B$.

\subsection{Alacaoglu et al. (2019)}
Recently Alacaoglu et al. also proposed a proof for the almost sure convergence of SPDHG (\cite{alacaoglu}, Theorem 4.4) using a strategy similar to ours. They describe SPDHG as a special case of a sequence $(\mathbf{w}^k)_{k\in\mathbb{N}}\subset\mathbb{R}^{dn+n^2}$ of the form (\ref{T}) for an operator $\textbf{T}:\mathbb{R}^{dn+n^2}\to\mathbb{R}^{dn+n^2}$, where $X=\mathbb{R}^d$ and $Y=\mathbb{R}^n$. They show that, under the step size condition (\ref{assu_i}),
the iterations $\mathbf{w}^k\in\mathbb{R}^{dn+n^2}$ converge to a fixed point of $\mathbf{T}$.
In contrast, we describe SPDHG using a random sequence $(T^k)_{k\in\mathbb{N}}$ of more intuitive operators $T^k:X\times Y\to X\times Y$, and then prove the sequence $(w^k)_{k\in\mathbb{N}}$ converges to a fixed point $w\in\bigcap_{k\in\mathbb{N}}\mathrm{Fix}T^k$.

\section{Numerical Examples}
\label{sec:numeric}

In this section we use MRI reconstruction as a case study to illustrate the performance of SPDHG in comparison to PDHG.
In parallel MRI reconstruction~\cite{mri}, a signal $x$ is reconstructed from multiple data samples $b_1,...,b_n$ of the form $b_i=A_ix+\eta_i$, where each $A_i:\mathbb{C}^d\to \mathbb{C}^m$ is an encoder operator from the signal space to the sample space, and $\eta_i\in \mathbb{C}^m$ represents noise added to the measurements.
A least-squares solution is given by
\begin{equation} \label{mri_model}
    \hat{x} \in \arg\min_x \sum_{i=1}^n \|A_ix-b_i\|^2 + g(x)
\end{equation}
where $g$ acts as a regularizer.
We recover our convex minimization template~(\ref{min}) by identifying $\mathbb{C}^d$ with $\mathbb{R}^{2d}$ and setting $X=\mathbb{R}^{2d}$, $Y_i=\mathbb{R}^{2m}$ and $ f_i(y) = \|y - b_i\|^2 $.

Here, we consider undersampled data with $n$ coils, i.e. $A_i = S\circ F\circ C_i$, where $S:\mathbb{C}^d\to\mathbb{C}^m$ is a subsampling operator, $F:\mathbb{C}^d\to\mathbb{C}^d$ represents the discrete Fourier transform and $C_ix = c_i\cdot x $ is the element-wise multiplication of $x$ and the $i$-th coil-sensitivity map $c_i\in \mathbb{C}^d$.
Data samples $b_i = A_ix^\dagger + \eta_i$ have been synthetically generated from an original reconstruction $x^\dagger$ from the BrainWave database \cite{MRIdata}.

In order to choose step-size parameters $\sigma_i$ and $\tau$ which comply with step-size condition (\ref{assu_i}), we consider $\sigma_i = \gamma(p_i/\|A_i\|)$ and $\tau = \gamma^{-1} (0.99/\max_i\|A_i\|)$, and then choose the value of $\gamma$ (among orders of magnitude $\gamma=10^{-5},10^{-4},...,10^{5}$) that gives the lowest objective $\Phi(x^k) = \sum_{i=1}^n f_i(Ax^k)+g(x^k)$ after 100 epochs.
Step sizes for PDHG have been optimized in the same way by setting $\sigma = \gamma/\|A\|$  and $\tau = \gamma^{-1}(0.99/\|A\|)$, which satisfies its step-size condition $\tau\sigma\|A\|^2<1$ \cite{chambollePock}.

\begin{figure}[t]
\centering
\subfloat[\scriptsize Ground truth  ]{\includegraphics[width=0.2\textwidth]{groundtruthA.png}} \quad
\subfloat[\scriptsize Target  ]{\includegraphics[width=0.2\textwidth]{targetA.png}} \quad
\subfloat[\scriptsize  SPDHG 4 coils \\ \centering{ 100 epochs}]{\includegraphics[width=0.2\textwidth]{m4_c3_0k_s.png}} \quad
\subfloat[\scriptsize  PDHG 4 coils \\ \centering{ 100 iterations}]{\includegraphics[width=0.2\textwidth]{m4_c3_0k_pd.png}} \\
\vspace{-3mm}
\subfloat[\scriptsize  SPDHG 8 coils \\ \centering{ 100 epochs}]{\includegraphics[width=0.2\textwidth]{m8_c3_0k_s.png}} \quad
\subfloat[\scriptsize  PDHG 8 coils \\ \centering{ 100 iterations}]{\includegraphics[width=0.2\textwidth]{m8_c3_0k_pd.png}} \quad
\subfloat[\scriptsize  SPDHG 8 coils \\ \centering{ 1,000 epochs}]{\includegraphics[width=0.2\textwidth]{m8_c3_1k_s.png}} \quad
\subfloat[\scriptsize  PDHG 8 coils \\ \centering{ 1,000 iterations}]{\includegraphics[width=0.2\textwidth]{m8_c3_1k_pd.png}}
\caption{Images reconstructed using different algorithms. Subfigure (a) shows the original ground truth $x^\dagger$ from which the noisy data samples have been synthetically generated, while subfigure (b) shows a reliable reconstruction $x^*$ to which the other figures should be compared. }
\label{fig: reconL2}
\end{figure}

Figure \ref{fig: reconL2} shows reconstructions obtained through SPDHG and PDHG.
Here, we considered model (\ref{mri_model}) with squared 2-norm regularizer $g(x) = 10^{-4}\|x\|^2$. The \emph{target} solution $x^*$ in (b) has been computed by running SPDHG for large number of epochs ($>10^4$). Comparing subfigures (e) and (f), the solution for SPDHG seems to be closer to the target (b) than that of PDHG at 100 epochs. At 1,000 epochs, solutions (g) and (h) appear visually similar to each other.

Figure \ref{fig: reconTV} shows reconstructions using model (\ref{mri_model}) with total-variation regularizer $g(x)=10^{-4}\|\nabla x\|_1$. As before, solutions by SPDHG appear closer to the target than those by PDHG. Notice as well how increasing the number of coils seems to yield more detailed although noisier reconstructions.

Figure \ref{fig: convergence plot} summarizes the performance of both algorithms for all of these reconstructions and some more.

\section{Conclusions}

In this paper we showed the almost sure convergence of SPDHG by using alternative arguments as provided by Alacaoglu et al \cite{alacaoglu}. We also investigated SPDHG in the context of parallel MRI where we observed a significant speed-up, thereby supporting the practical use of variational regularization methods for this application.

\begin{figure}[ht]
\centering
\subfloat[\scriptsize Ground truth ]{\includegraphics[width=0.2\textwidth]{groundtruthB.png}} \quad
\subfloat[\scriptsize Target ]{\includegraphics[width=0.2\textwidth]{TV8_c3_10k_s.png}} \quad
\subfloat[\scriptsize  SPDHG 4 coils \\ \centering{ 100 epochs}]{\includegraphics[width=0.2\textwidth]{TV4_c3_1k_s.png}} \quad
\subfloat[\scriptsize  PDHG 4 coils \\ \centering{ 100 iterations}]{\includegraphics[width=0.2\textwidth]{TV4_c3_1k_pd.png}} \\
\vspace{-3mm}
\subfloat[\scriptsize  SPDHG 8 coils \\ \centering{ 100 epochs}]{\includegraphics[width=0.2\textwidth]{TV8_c3_0k_s.png}} \quad
\subfloat[\scriptsize  PDHG 8 coils \\ \centering{ 100 iterations}]{\includegraphics[width=0.2\textwidth]{TV8_c3_0k_pd.png}} \quad
\subfloat[\scriptsize  SPDHG 8 coils \\ \centering{ 1,000 epochs}]{\includegraphics[width=0.2\textwidth]{TV8_c3_1k_s.png}} \quad
\subfloat[\scriptsize  PDHG 8 coils \\ \centering{ 1,000 iterations}]{\includegraphics[width=0.2\textwidth]{TV8_c3_1k_pd.png}}
\caption{Images reconstructed using noisy data samples synthetically generated from the ground truth~(a). Subfigure (b) shows the target reconstruction for comparison. }
\label{fig: reconTV}
\end{figure}

\begin{figure}[ht]
\centering
\subfloat[\scriptsize  $L^2$-regularizer]{\includegraphics[width=0.53\textwidth]{m_c3_10k.png}} \quad
\subfloat[\scriptsize  $TV$-regularizer]{\includegraphics[width=0.44\textwidth]{TV_c3_1k.png}}
\caption{ Relative objective $\Phi_r(x^k) = \frac{\Phi(x^k) - \Phi(x^*)}{\Phi(x^0) - \Phi(x^*)}$ at every epoch for image reconstruction, including the examples from Figures \ref{fig: reconL2} and \ref{fig: reconTV}.}
\label{fig: convergence plot}
\end{figure}

\bibliographystyle{splncs04}