text:
\section{Introduction}
    \label{sec:Introduction}

In various applications, datasets are stored in a distributed manner among several sites or machines \citep[chap. 1.2]{fan2020statistical}.
Often, due to communication constraints as well as privacy restrictions, the raw data cannot be shared between the various machines.
Such settings have motivated the development of methods and supporting theory for distributed learning and inference. See, e.g., the reviews by \citet{huo2019aggregated}, \citet{gao2022review} and references therein.

In this paper we consider distributed learning of a  sparse linear regression model. Specifically,
we assume that the response $y\in\mathbb{R}$
and the vector $X\in\mathbb{R}^d$ of explanatory variables are linearly related via
\begin{align}
        y = X^{\top}\theta^* + w,
            \label{e:linear_model}
\end{align}
where $w\sim N(0,\sigma^2)$,   $\sigma>0$ is the noise level, and $\theta^*\in\re^d$ is an unknown vector of coefficients.
We further assume $X\in\re^{d}$ is random with mean zero and covariance matrix $\Sigma$.
We focus on a high-dimensional setting $d\gg 1$, and assume that $\theta^*$ is sparse with only $K\ll d$ nonzero coefficients.
The support set of $\theta^*\in\mathbb{R}^d$ is denoted
by $\mathcal{S}=\{i\in[d]\,|\, |\theta_i^*|>0  \}$.

Given $N$ samples $\{(X_i,y_i)\}_{i=1}^N$
from the model (\ref{e:linear_model}),
common tasks are to estimate
the vector $\theta^*$ and its support set $\cal S$.
Motivated by contemporary applications, we consider these tasks in a distributed setting where the data are randomly split among $M$ machines. Specifically,
we consider a star topology network, whereby the end machines
communicate only with a fusion center.

As reviewed in Section \ref{sec:background_bibreview},
estimating $\theta^*$ and its support in the above or similar distributed settings were studied by several authors, see for example  \citet{mateos2010distributed, chen2014split, lee2017communication, battey2018distributed,chen2020distributed, liu2021integrative, barghi2021distributed}
and references therein.
Most prior works on distributed regression
required communication of at least $O(d)$ bits per machine,  as in their schemes each machine
sends to the fusion center its full $d$-dimensional estimate of the unknown vector $\theta^*$. Some works in the literature denote this as communication efficient, in the sense that for a machine holding $n$ samples, an $O(d)$ communication is still significantly less than the size {$O(n\cdot d)$} of its data.

The design and analysis of communication efficient distributed schemes is important, as in various distributed settings the communication channel is the critical bottleneck.
Moreover, in some practical cases, such as mobile devices and sensor networks, the end machines may have very limited bandwidth.
Thus, in high dimensional settings with $d\gg 1$, it may not even be feasible for each machine to send messages of length $O(d)$.
In this work, we study such a restricted communication setting, assuming that each
machine is allowed to send to the fusion center only a limited number of bits, significantly lower than the dimension $d$.
Our goals are to develop low communication distributed schemes to estimate $\theta^*$ and its support and to theoretically analyze their performance.

We make the following contributions. On the methodology side, in Section \ref{sec:proposed_method} we present
several two round distributed schemes. The schemes vary slightly by the messages sent, but in all of them, the fusion center estimates the support set of $\theta^*$ in the first round and the vector $\theta^*$ in the second round.
In our schemes, each machine computes its own debiased lasso estimate. However, it sends to the center only the indices of its top few largest values, possibly along with their signs.
Hence, the communication per machine is significantly less than $d$ bits.
In the simplest variant, the fusion center estimates the support of $\theta^*$ by voting,
selecting the few indices that were sent by the largest number of machines.

Next, on the theoretical side,
we prove in Section \ref{sec:theoretical_results} that under suitable conditions,
with high probability
the first round of our scheme achieves exact support recovery with communication per machine sublinear in $d$. Specifically, we present support guarantees under two different parameter regimes.
Theorem \ref{teo:supprec_threshold_votes} considers a case with a relatively large
number of machines. Here, each machine sends a short message of $O(K \ln d)$ bits.
Next,
Theorem \ref{teo:supprec_threshold_votes_M_low}  considers a setting with  relatively few machines, $M=O(\ln d)$. Here, to achieve exact support recovery each machine sends a much longer message, of length $O(d^\alpha)$ for some suitable $\alpha<1$. This is still sublinear in $d$, and much less than the communication required if a machine were to send its full $d$-dimensional estimated vector.
The proofs of our theorems rely on recent results regarding the distribution of debiased lasso estimators, combined with sharp bounds on tails of binomial random variables.
Exact support recovery follows by showing that
with high probability,
all non-support indices receive fewer votes than support indices.

In Section \ref{sec:simulations} we  present simulations
comparing our schemes to previously proposed methods.
These illustrate that with our algorithms, the fusion center correctly detects the support of $\theta^*$ and consequently accurately estimates $\theta^*$, even at low signal to noise ratios where each machine is unable to do so.
Furthermore, this is achieved with very little communication per machine  compared to the dimension $d$.
One insight from both the simulations and our theoretical analysis is that for the fusion center to detect the correct support, it is not necessary to require $M/2$ votes as suggested in \citet{barghi2021distributed} and \cite{chen2014split}. Instead, {as few as $O(\ln d)$ votes suffice to distinguish support from non-support indices.}
Interestingly, under a broad range of parameter values, our schemes work as well as, and in some cases better than  more communication intensive approaches.
Our simulations also highlight the importance and advantages of a second round of communication.
Specifically, even though a single-round scheme based on averaging debiased lasso estimates, as proposed by
\citet{lee2017communication},
is minimax rate optimal and finds the correct support, it nonetheless may output an estimate with a larger mean squared error than that of our scheme.
We conclude with a summary and discussion in Section \ref{sec:summary}.
Proofs appear in the Appendix.

\paragraph{\bf Notation}
For an integer $k\geq 1$, we denote $[k]=\{1,2,\ldots,k\}$. The indicator function is denoted as $I(A)$, which equals one if condition $A$ holds and zero otherwise.
The $\ell_q$ norm of a vector $Y\in\re^n$ for {$q\geq 1$} is
$\|Y\|_q = (\sum_{i=1}^{n} |Y_i|^q)^{1/q}$,
whereas $\|Y\|_0 = \sum_{i=1}^{n} I(Y_i\neq 0)$ is its number of {nonzero} entries.
We denote by $|Y|$ the vector whose entries are $(|Y_1|,|Y_2|,\ldots,|Y_n|)$.
For a $d\times d$ matrix ${\bf A}=\{a_{ij}\}_{i,j=1}^d$, we denote $\| {\bf A}\|_{\infty} = \max_{1\leq i\leq d} \sum_{j=1}^d |a_{ij}|$.
We further denote by
$\sigma_{\min}({\bf A})$ and $\sigma_{\max}({\bf A})$ its  smallest and largest singular values, respectively.
For a subset $J\subset [d]$, ${\bf A}_J$ is the $d\times |J|$ matrix whose columns are those in the subset $J$. Similarly,
${\bf A}_{J,J}$ is the $|J|\times |J|$ submatrix whose rows and columns correspond to the indices in $J$.
The cumulative distribution function {(CDF)} of a standard Gaussian is denoted by $\Phi(\cdot)$ whereas $\Phi^c(\cdot)=1-\Phi(\cdot)$.
We write $a_n \gtrsim b_n$ for two sequences $\{a_n\}_{n\geq1}$ and $\{b_n\}_{n\geq1}$ if there are positive constants $C$ and $n_0$ such that $a_n \geq C b_n$ for all $n>n_0$.

\section{Previous works}
    \label{sec:previous_works}

\label{sec:background_bibreview}

Distributed linear regression schemes under various settings, not necessarily involving sparsity, have been proposed and theoretically studied in multiple fields, including sensor networks, statistics and machine learning, see for example  \citep{guestrin2004distributed,predd2006distributed, boyd2011distributed,zhang2013communication,heinze2014loco, rosenblatt2016optimality, jordan2019communication,chen2020distributed, dobriban2020wonder, zhu2021least, dobriban2021distributed}.

\cite{mateos2010distributed} were among the first to study distributed sparse linear regression in a general setting without a fusion center, where machines are connected and communicate with each other.
They devised a multi-round scheme whereby all the machines reach a consensus and jointly approximate the centralized solution, that would have been computed if all data were available at a single machine.
Several later works focused on the setting which we also consider in this paper, where machines are connected in a star topology to a fusion center, and only one or two communication rounds are made.
In a broader context of generalized sparse linear models,
\cite{chen2014split} proposed a divide-and-conquer approach where each machine
estimates $\theta^*$ by minimizing a penalized
objective with a sparsity inducing penalty, such as $\|\theta \|_1$.
Each machine sends its sparse estimate to the fusion center, which estimates the support by voting over the indices of the individual estimates of the $M$ machines.
Finally, the center estimates
$\theta^*$ by a weighted
average of these $M$ estimates.
For sparse linear regression, with each machine computing a lasso estimate of $\theta^*$, their method suffers from the well known  bias of the lasso, which is not reduced by averaging.

To overcome the bias of the lasso, in recent years several debiased lasso estimators were derived and theoretically studied, see \citet{zhang2014confidence,van2014asymptotically,javanmard2018debiasing}.
For distributed learning,
debiased estimators have been applied in various settings,
including hypothesis testing, quantile regression and more, see for example \citet{lee2017communication,battey2018distributed,liu2021integrative,lv2022debiased}.

In particular,
\citet{lee2017communication} proposed a single round scheme whereby each machine computes its own debiased lasso estimator, and sends it to the fusion center.
The center averages these debiased estimators and thresholds the result to estimate $\theta^*$ and recover
its support.
\cite{lee2017communication} proved that the resulting estimator achieves the same error rate as the centralized solution, and is minimax rate optimal.
However, their scheme requires a communication of $O(d)$ bits per machine and is thus not applicable in the restricted communication setting considered in this manuscript.
Moreover, as we demonstrate in the simulation section, unless the signal strength is very low, our two round scheme in fact achieves a smaller mean squared error, with a much lower communication. This highlights the potential sub-optimality of lasso and debiased lasso in sparse regression problems with sufficiently strong signals.

Most related to our paper is the recent work by \citet{barghi2021distributed}. In their method, each machine computes a debiased lasso estimator $\hat{\theta}$, but sends to the fusion center only the indices $i$ for which $|\hat{\theta}_i|$ is above a certain threshold.
The support set estimated by the fusion center consists of all indices that were sent by at least half of the machines, i.e., indices that received at least $M/2$ votes. Focusing on the consistency of feature selection,
\citet{barghi2021distributed} derive bounds on the type-I and type-II errors of the estimated support set.
Their results, however, are given as rates with unspecified multiplicative constants. As we show in this work, both theoretically and empirically, consistent support estimation is possible with a much lower voting threshold.
Furthermore, requiring at least $M/2$ votes implies that their scheme achieves exact support recovery only for much stronger signals.

We remark that voting is a natural approach for distributed support estimation under  communication constraints. \cite{amiraz2022distributed} analyzed voting-based distributed schemes in the context of a simpler problem of sparse Gaussian mean estimation.
They proved that even at low signal strengths, their schemes achieve exact support recovery with high probability using communication sublinear in the dimension.
Their setting can be viewed as a particular case of sparse linear regression but with a unitary design matrix. Their proofs, which rely on this property, do not extend to our setting.

\section{The lasso and debiased lasso estimators}
\label{sec:lasso_and_debiased_lasso}

For our paper to be self contained, we briefly review the lasso and debiased lasso and some of their theoretical properties.
The lasso \citep{tibshirani1996regression} is perhaps the most popular method to fit high-dimensional sparse linear models.
Given
a regularization parameter $\lambda>0$ and $n$ samples $(X_i, y_i)$, stacked in a design matrix ${\bf X}\in\re^{n\times d}$ and a response vector $Y\in\re^n$, the lasso estimator is
given by
\begin{align}
\tilde{\theta} = \tilde{\theta}({\bf X}, Y,\lambda)
= \arg\min_{\theta\in\re^{d}} \left\{ \frac{1}{2n}\|Y - {\bf X}\theta\|_2^2 + \lambda\|\theta\|_1 \right\}.
\label{e:lasso_estimator}
\end{align}
The lasso has two desirable properties.
First, computationally
Eq. (\ref{e:lasso_estimator}) is a convex problem for which there are fast solvers.
Second, from a theoretical standpoint,
it enjoys strong recovery guarantees,
assuming the data follows the model (\ref{e:linear_model}) with an exact or approximately sparse $\theta^*$, see for example
\citep{candes2005decoding,bunea2007sparsity,van2009conditions,hastie2015statistical}.
However, the lasso has two major drawbacks: it may output significantly biased estimates and it does not have a simple asymptotic distribution.
The latter is needed
for confidence intervals and hypothesis testing.
To overcome these limitations, and in particular derive confidence intervals for high-dimensional sparse linear models, several authors developed debiased lasso estimators \citep{zhang2014confidence, van2014asymptotically, javanmard2014hypothesis, javanmard2014confidence, javanmard2018debiasing}.

For random ${\bf X}$ with a known population
covariance matrix $\Sigma$, \cite{javanmard2014hypothesis} proposed $\frac{1}{n}\Sigma^{-1}{\bf X}^{\top}(Y - {\bf X}\tilde{\theta})$ as a debiasing term. .
As $\Sigma$ is often unknown, both \cite{van2014asymptotically} and \cite{javanmard2014confidence} developed methods to estimate its inverse $\Omega=\Sigma^{-1}$.
In our work, we estimate $\Omega$ using the approach of \cite{van2014asymptotically}, who assume that $\Omega$ is sparse.
In their method, presented in Algorithm \ref{alg:decorr_matrix},
$\hat{\Omega}$ is constructed by fitting a lasso regression with regularization $\lambda_{\Omega}>0$ to each column of ${\bf X}$ against all the other columns. {Hence,
it
requires solving $d$ separate lasso problems.}

Given the lasso estimate
$\tilde\theta$ of Eq. (\ref{e:lasso_estimator})
and the matrix
$\hat{\Omega}$, the debiased lasso is
\begin{align}
\hat{\theta} = \hat{\theta}(Y, {\bf X},\lambda,\lambda_{\Omega}) = \tilde{\theta} + \frac{1}{n}\hat{\Omega}{\bf X}^{\top}(Y - {\bf X} \tilde{\theta}).
\label{e:debiased_lasso}
\end{align}
An appealing property of $\hat{\theta}$ is that, under some conditions, it
is asymptotically unbiased with a Gaussian distribution.
For our analysis, we shall use the
following result \citep[Theorem 3.13]{javanmard2018debiasing}.

\begin{theorem}
    \label{teo:deblasso_asymp_dist}
Consider the linear model $Y = {\bf X}\theta^* + W$\!, where $W\sim N(0,\sigma^2{\bf I}_{n\times n})$ and ${\bf X}\in\re^{n\times d}$ has independent Gaussian rows with zero mean and covariance matrix $\Sigma\in\re^{d\times d}$. Suppose that $\Sigma$ satisfies the following conditions:
\begin{enumerate}[i.]
\item For all $i\in[d]$, $\Sigma_{ii}\leq 1$.
\item For some constants $C_{\max},C_{\min}>0$,
\begin{align}
    0<C_{\min}<\sigma_{\min}(\Sigma)\leq \sigma_{\max}(\Sigma)<C_{\max}.
    \label{e:singular_values_Sigma}
\end{align}
\item For $C_0=(32 C_{\max}/C_{\min}) + 1$,
and a constant $\rho > 0$,
$$\max_{J\subseteq [d],\, |J|\leq C_0 K}\| \Sigma_{J,J}^{-1}\|_{\infty} \leq \rho.$$
\end{enumerate}
Let $K_{\Omega}$ be the maximum row-wise sparsity of $\Omega=\Sigma^{-1}$, that is,
\begin{align*}
    K_{\Omega}=\max_{i\in [d]}\left| \{j\in[d]; \,\Omega_{ij}\neq 0, \, j\neq i\} \right|.
\end{align*}
Let $\tilde{\theta}$ be the lasso estimator computed using $\lambda = \kappa\sigma\sqrt{(\ln d)/n}$ for $\kappa\in [8, \kappa_{\max}]$, and let $\hat{\theta}$ be the debiased lasso estimator in Eq. \eqref{e:debiased_lasso} with $\hat{\Omega}$ computed by  Algorithm \ref{alg:decorr_matrix} with
$\lambda_{\Omega}=\kappa_{\Omega}\sqrt{(\ln d)/n}$ for some suitable large $\kappa_{\Omega}>0$. Let $\hat{\Sigma}={\bf X}^{\top}{\bf X}/n$ denote the empirical covariance matrix.
Then there exist constants $c$, $c_*$, $C$ depending solely on $C_{\min}$, $C_{\max}$, $\kappa_{\max}$ and $\kappa_{\Omega}$ such that, for $n \geq c \max\{K, K_{\Omega}\}\ln d$, the following holds:
\begin{align}
\sqrt{n}(\hat{\theta} - \theta^*) = Z + R, \quad Z|X \sim N(0, \sigma^2 \hat{\Omega}\hat{\Sigma}\hat{\Omega}^{\top}), \label{e:deblasso_estimator_representation_Z_R}
\end{align}
where $Z = n^{-1/2}\hat{\Omega}{\bf X}^{\top}W$ and $R = \sqrt{n}\left( \hat{\Omega}\hat{\Sigma} - {\bf I} \right)(\theta^* - \tilde{\theta})$, and with probability at least $1 - 2 d e^{-c_* n/K} - d e^{-c n} - 6 d^{-2}$,
\begin{equation}
\|R\|_{\infty} \leq
C\sigma\frac{\ln d}{\sqrt{n}}\left(
{\rho\sqrt{K}} + \min\{K, K_{\Omega}\}
\right).
\label{e:bound_bias_term_R}
\end{equation}
\end{theorem}

Assumptions $(i)$ and $(ii)$ in this theorem are common in the literature. Assumption $(iii)$ is satisfied, for example, by circulant matrices {$\Sigma_{ij}=\varsigma^{|i-j|}$, $\varsigma\in(0,1)$}.
The quantity $R$ in Eq. \eqref{e:deblasso_estimator_representation_Z_R} can be viewed as a bias term. By Theorem \ref{teo:deblasso_asymp_dist}, this bias is small if the sample size and dimension are suitably large, which in turn
implies that $\hat\theta_i$ is approximately Gaussian.
The following lemma, proven in the Appendix, bounds the error of this approximation. It will be used
in analyzing the probability of exact support recovery of our distributed scheme.

\begin{algorithm}[tp]
\caption{Computation of a precision matrix estimate $\hat{\Omega}$}\label{alg:decorr_matrix}
\begin{algorithmic}[1]
\vspace{0.5em}
\Require design matrix ${\bf X}\in\re^{n\times d}$, regularization parameter $\lambda_{\Omega}>0$.
\Ensure precision matrix estimate $\hat{\Omega}\in\re^{d\times d}$.
\Statex $\textbf{x}_i\in \re^{n}$ denotes the $i$-th column of ${\bf X}$.
\Statex ${\bf X}_{-i}\in \re^{n\times (d-1)}$ denotes the design matrix with the $i$-th column removed.
\For{$i=1,\ldots,d$}
	\State Fit a lasso with response
 $\textbf{x}_i$, design matrix ${\bf X}_{-i}$ and regularization parameter $\lambda_\Omega$.
    \State Let $\tilde{\gamma}_{i}=\{\tilde{\gamma}_{i,j}\}_{j=1,j\neq i}^d\in\re^{d-1}$ be the estimated regression coefficients of step 2.
	\State Compute  $\tilde{\tau}_i^2 = (2n)^{-1}\|\textbf{x}_i - {\bf X}_{-i}\tilde{\gamma}_{i}\|_2^2 + \lambda_{\Omega}\|\tilde{\gamma}_{i}\|_1$, $i\in[d]$.
\EndFor
\State Construct $d\times d$ matrix
\begin{align*}
\tilde{{\bf C}} = \left[\begin{array}{cccc}
1 & -\tilde{\gamma}_{1,2} & \cdots & -\tilde{\gamma}_{1,d} \\
-\tilde{\gamma}_{2,1} & 1 & \cdots & -\tilde{\gamma}_{2,d} \\
\vdots & \vdots & \ddots & \vdots \\
-\tilde{\gamma}_{d,1} & -\tilde{\gamma}_{d,2} & \cdots & 1 \\
\end{array}
\right].
\end{align*}
\State \Return $\hat{\Omega} = \mathrm{diag}\{\tilde{\tau}^{-2}_1,\ldots,\tilde{\tau}^{-2}_d\}\tilde{{\bf C}}$.
\end{algorithmic}
\end{algorithm}

\begin{lemma}
Under the assumptions of Theorem \ref{teo:deblasso_asymp_dist},
for any $\tau > 0$,
\begin{align}
\left| \Pr\left( \tfrac{\sqrt{n}(\hat{\theta}_i - \theta_i^*)}{\sigma\sqrt{c_{ii}}}  \leq \tau \right) - \Phi\left( \tau \right) \right|
\leq \frac{\delta_R}{\sigma\sqrt{c_{ii}}} \phi\left( \tau \right)
+  2de^{-c_* n/K} + de^{-c n} + \frac{6}{d^2},
\label{e:approx_cdf_dist_deblasso_estimator}
\end{align}
where $\phi(\cdot)$ denotes the Gaussian density function, $c_{ii} = (\hat{\Omega}\hat{\Sigma}\hat{\Omega}^{\top})_{ii}$, and $\delta_R$ is the upper bound on the bias term in Eq. \eqref{e:bound_bias_term_R}, namely
\begin{align}
\delta_R =
C\sigma\frac{\ln d}{\sqrt{n}}\left(
{\rho\sqrt{K}} + \min\{K, K_{\Omega}\}
\right).
\label{e:bound_bias_term_deltaR}
\end{align}
\label{lemma:approx_cdf_dist_deblasso_estimator}
\end{lemma}

\section{Distributed sparse regression with restricted communication}
\label{sec:proposed_method}

As described in Section \ref{sec:Introduction},  we consider a distributed setting with $M$ machines connected in a star topology to a fusion center.
For simplicity, we assume that each machine $m$ has a sample $({\bf X}^{m}, Y^{m})$ of $n=N/M$
i.i.d.
observations from the model (\ref{e:linear_model}), where $Y^{m}\in\re^n$ and ${\bf X}^{m}\in\re^{n\times d}$.
In describing our schemes, we further assume that the noise level $\sigma$ is known. If $\sigma$ is unknown,
it may be
consistently estimated, for example, by
the scaled lasso of \cite{sun2012scaled},  see also  \cite[Corollary 3.10]{javanmard2018debiasing}.

We present several two round distributed schemes to estimate the sparse vector $\theta^*$ of Eq. (\ref{e:linear_model}) under the constraint of limited communication between the $M$ machines and the fusion center.
Here we present the simplest scheme and discuss other variants in section \ref{subsec:variations_supp_estim}. In all variants,
the fusion center estimates the support of
$\theta^*$ in the first round, and
$\theta^*$ itself in the second round.

The first round of our scheme is described in Algorithm \ref{alg:count_votes}, whereas
the full
two round scheme is outlined in Algorithm \ref{alg:two_round_scheme}.
In the first round,
each machine $m\in[M]$ computes the following quantities
using its own data $({\bf X}^m, Y^m)$:
(i)~a lasso estimate $\tilde{\theta}^m$ by Eq. \eqref{e:lasso_estimator}; (ii) a  matrix $\hat{\Omega}^m$ by Algorithm \ref{alg:decorr_matrix}; and (iii) a debiased lasso $\hat{\theta}^m$ by Eq. \eqref{e:debiased_lasso}.
Up to this point, this is identical to  \citet{lee2017communication}.
The main difference is that in their scheme, each machine sends
to the center its
debiased lasso estimate $\hat{\theta}^m \in \mathbb{R}^d$, incurring $O(d)$ bits of communication.

In contrast, in our scheme each machine sends only a few indices.
Towards this end and
{in light of Eq. \eqref{e:approx_cdf_dist_deblasso_estimator} of Lemma \ref{lemma:approx_cdf_dist_deblasso_estimator},}
each machine computes a normalized vector $\hat{\xi}^m$ whose coordinates are given by
\begin{align}
    \hat{\xi}_k^{m} = \frac{\sqrt{n}\hat{\theta}^{m}_k}{\sigma(\hat{\Omega}^m\hat{\Sigma}^{m}(\hat{\Omega}^m)^{\top})_{kk}^{1/2}}
        ,
    \qquad  \forall k\in[d].
\label{e:standardized_debiased_lasso}
\end{align}
In the simplest variant, each machine  sends to the center only
indices $k$ such that $|\hat{\xi}_k^{m}|>\tau$ for some suitable threshold $\tau>0$.

\begin{algorithm}[t]
\caption{Distributed voting based scheme for support estimation}\label{alg:count_votes}
\begin{algorithmic}[1]
\vspace{1em}
\Require Data $({\bf X}^m, Y^m)\in\re^{n\times(d+1)}$, threshold $\tau$ and regularization parameters $\lambda_{\Omega}$ and $\lambda$.
\Ensure Support estimate $\hat{\mathcal{S}}$.
\smallskip
\Statex \textbf{At each local machine $m=1,\ldots,M$}
\State Compute a lasso estimator $\tilde{\theta}^m$ via Eq. (\ref{e:lasso_estimator}) with regularization parameter $\lambda$.
\State Compute a precision matrix estimate $\hat{\Omega}^m\in\re^{d\times d}$ by  Algorithm \ref{alg:decorr_matrix} with ${\bf X}^m$ and $\lambda_{\Omega}$.
\State Compute a debiased lasso estimate $\hat{\theta}^m\in\re^{d}$, Eq.  \eqref{e:debiased_lasso}, with data {$({\bf X}^m, Y^m)$}, $\lambda$ and $\hat{\Omega}^m$.
\State Calculate the empirical covariance matrix $\hat{\Sigma}^{m}=n^{-1}({\bf X}^m)^{\top} {\bf X}^m$.
\State Use $\hat{\Omega}^m$ and $\hat{\Sigma}^{m}$ to compute the standardized estimator $\hat{\xi}^{m}\in\re^d$,
Eq. \eqref{e:standardized_debiased_lasso}.
\State Set $\mathcal S^{m} = \{i; |\hat{\xi}^{m}_{i}| > \tau \}$ and send it to the fusion center.
\smallskip
\Statex \textbf{At the fusion center}
\State For each $i\in[d]$, compute $V_i = \sum_{m=1}^M I(i \in \mathcal S^m)$.
\State Sort $V_{j_1}\geq V_{j_2}\geq \cdots\geq V_{j_d}$.
\State \Return $\hat{\mathcal{S}}=\{j_1, \ldots, j_K \}$.
\end{algorithmic}
\end{algorithm}

Given the messages sent by the $M$ machines, the fusion center counts the number of votes received by each index.
If the sparsity level $K$ is known, its estimated support set $\hat{\mathcal{S}}$ consists of the $K$ indices with the largest number of votes.
Otherwise, as discussed
in Remark \ref{remark:unknown_sparsity}
below,
the center may estimate the support set by the indices whose number of votes exceed a suitable threshold.

Next, we describe the second round.
At its start, the fusion center sends the estimated support $\hat{\mathcal{S}}$ to all $M$ machines.
Next, each machine
computes the standard least squares regression solution, restricted to the set $\hat{\cal{S}}$, namely
\begin{equation}
    \hat\beta^{m} = \arg\min_{\beta}
        \| {\bf X}_{\hat{\mathcal{S}}}^m \beta - Y^m\|_2^2
            \label{e:beta_m}
    \end{equation}
where ${\bf X}_{\hat{\mathcal{S}}}^m \in \mathbb{R}^{n\times |\hat{\cal S}|}$ consists
of the columns of ${\bf X}^m$ corresponding to the indices in
$\hat{\cal S}$.
Each machine then sends its vector
$\hat \beta^m$
to the fusion center.
Finally,
the fusion center estimates $\theta^*$
by averaging these $M$ vectors,
\begin{equation}
    \hat \theta_i =
    \left\{
        \begin{array}{cc}
             \frac1M \sum_{m=1}^M \hat \beta^m_i &  i\in \hat{\cal S} \\
             0  & \mbox{otherwise}
        \end{array}
    \right.
            \label{e:theta_hat}
\end{equation}

In the next section we present several variants of this basic two round scheme.
Before that we make a few remarks and observations.

\begin{algorithm}[t]
\caption{Two round distributed scheme to estimate $\theta^*$}\label{alg:two_round_scheme}
\begin{algorithmic}[1]
\vspace{1em}
\Require Data $({\bf X}^m, Y^m)\in\re^{n\times(d+1)}$, sparsity $K$, threshold $\tau$ and regularizations $\lambda_{\Omega}$ and $\lambda$.
\Ensure A two-round estimate $\hat{\theta}\in\re^d$ of $\theta^*$.
\smallskip
\Statex \underline{First round}
\smallskip
\State The fusion center estimates $\hat{\mathcal{S}}$ with Algorithm \ref{alg:count_votes}.
\smallskip
\Statex \underline{Second round}
\smallskip
\State The fusion center sends $\hat{\mathcal{S}}$ to all $M$ machines.
\smallskip
\Statex \textbf{At each local machine $m=1,\ldots,M$}
\State Let ${\bf X}_{\hat{\mathcal{S}}}^m \in \re^{n\times K}$ be the $K$ columns of ${\bf X}^m$ corresponding to indices in $\hat{\mathcal{S}}$.
\State Compute $\hat{\beta}^m =
 \arg\min_\beta \| {\bf X}_{\hat{\mathcal{S}}}^m \beta - Y^m  \|_2^2$.
\State Send $\hat{\beta}^m$ to the fusion center.
\smallskip
\Statex \textbf{At the fusion center}
\State Given $\hat\beta^1,\ldots, \hat\beta^M$, compute the estimate
$\hat\theta$ according to Eq. \eqref{e:theta_hat}.
\State \Return $\hat{\theta}$.
\end{algorithmic}
\end{algorithm}

\begin{remark}
The communication of
the first round (Algorithm \ref{alg:count_votes}) depends on
the threshold $\tau$. A high threshold
leads to only few sent indices. However,
at low signal strengths,
the signal coordinates may not have the highest
values
$|\hat{\xi}^m_k|$ and thus may not be sent.
Hence,
for successful support recovery by the fusion center,
a lower threshold leading to many more sent coordinates
is required.
Since the maxima of $d$ standard Gaussian variables scales as $\sqrt{2\ln d}$, to comply with the communication constraints, the threshold $\tau$ should also scale as $O(\sqrt{\ln d})$.
In Section \ref{sec:theoretical_results}, we present suitable thresholds and
sufficient conditions on the
number of machines and on the signal strength, which guarantee support recovery by
Algorithm \ref{alg:count_votes}, with high probability and little communication per machine.
\end{remark}

\begin{remark}
With known $K$,
the communication per machine of the second round is $O(K\ln d)$ bits. For suitable choices of the threshold $\tau$ in the first round, this is negligible or at most comparable to the communication of the first round.
\end{remark}

\begin{remark}
A two round scheme, whereby given an estimated set $\hat{\cal S}$, the second round is identical to ours
was discussed by
\citet{battey2018distributed} in Eq. (A.2) of their supplementary. The difference is that in their first round, similar to \citet{lee2017communication}, each machine sends its full debiased lasso vector, with a communication of $O(d)$ bits.
\citet{battey2018distributed} showed that, under certain conditions, their two-round estimator attains an optimal rate.
In Section \ref{sec:theoretical_results}, {we prove} that for a sufficiently high SNR, our method achieves the same rate, but using much less communication.
\label{remark:two_round_estimator_Battey_etal_2018}
\end{remark}

\begin{remark}
    With a higher communication per machine in the second round, it is possible for the fusion center to compute the exact centralized least squares solution corresponding to the set $\hat{\cal S}$, denoted $
    \hat\theta^{\mbox{\tiny LS}}$.
    Specifically, suppose that each machine sends to the center
    both   the vector
    $ ({\bf X}_{\hat{\mathcal{S}}}^m)^\top Y^m$
    of length $|\hat{\cal S}|$, and the
    $|\hat{\cal S}| \times |\hat{\cal S}|$
    matrix $ ({\bf X}_{\hat{\mathcal{S}}}^m)^\top
    {\bf X}_{\hat{\mathcal{S}}}^m$.
    The center may then compute $\hat\theta^{\mbox{\tiny LS}}$
as follows
    \begin{equation}
        \hat\theta^{\mbox{\tiny LS}} = \left(
        \sum_{m=1}^M
        ({\bf X}_{\hat{\mathcal{S}}}^m)^\top
    {\bf X}_{\hat{\mathcal{S}}}^m
        \right)^{-1}
        \sum_{m=1}^M ({\bf X}_{\hat{\mathcal{S}}}^m)^\top Y^m .
        \label{e:centralized_ols_estimator}
    \end{equation}
With $K$ known and $|\hat{\cal S}|=K$, such a second round has a communication of $O(K^2)$ bits.
If the sparsity $K$ is non-negligible, this is much higher than the $O(K)$ bits of our original scheme.
In particular, if $K=O(d^{1/2})$, the
resulting communication
is
comparable to that of sending the full debiased lasso vector.
\end{remark}

\begin{remark}
In practice, the sparsity $K$ is often unknown. Instead of step 9 in Algorithm \ref{alg:count_votes}, one alternative is to estimate $\mathcal{S}$ by thresholding the number of votes. For some threshold $\tau_{votes}>0$, $\hat{\mathcal{S}}$ could be set as all indices $i$ such that $V_i > \tau_{votes}$.
{Lemma \ref{lemma:prob_max_nonsupp_votes_M_logaritmic} in Appendix \ref{sec:proofs} shows that, under suitable conditions, non-support indices have a small probability of receiving more than $2\ln d$ votes. Hence, $\tau_{votes} = 2\ln d$ is a reasonable choice for such a threshold value.}
\label{remark:unknown_sparsity}
\end{remark}

\subsection{Variations of Algorithm \ref{alg:count_votes}}
\label{subsec:variations_supp_estim}

Various adaptations of Algorithm \ref{alg:count_votes} are
possible and may offer better performance.
One example is a top $L$ algorithm where each machine sends to the center the indices of the $L$ largest
entries of $|\hat\xi^m|$, for some parameter $K\leq L \ll d$.
A similar approach was proposed in
\citet{amiraz2022distributed} for the simpler problem of sparse normal means estimation.
One advantage of this variant is that its communication per machine is fixed and known a priori $O(L \ln d)$. This is in contrast to the above
thresholding based
scheme, whose communication per machine is random.

A different variant is to use sums of signs
to estimate the support.
Here
machines send both the indices corresponding to the largest entries in $|\hat{\xi}^m|$ and their signs. Hence, in  step 5 of Algorithm \ref{alg:count_votes}
the message sent by machine $m$ is
\begin{align*}
    \mathcal{S}^m = \left\{ \left(i, \mathrm{sign}(\hat{\xi}_i^m) \right); |\hat{\xi}_i^m|>\tau \right\}.
\end{align*}
Next, the fusion center computes for each index $i\in[d]$ its corresponding sum of received signs, i.e.,
\begin{align}
    V_i^{sign} = \sum_{m=1}^M \mathrm{sign}(\hat{\xi}_i^m) I\left\{ \left(i, \mathrm{sign}(\hat{\xi}_i^m)\right) \in \mathcal{S}^m \right\}.
    \label{e:sum_signs}
\end{align}
For known $K$, the estimated support set are
the $K$ indices with largest values of $|V_i^{sign}|$.
This algorithm uses a few more bits than a voting scheme. However, sums of signs are expected to better distinguish between support and non-support coefficients when the number of machines
is large.
The reason is that
at non-support indices $j\not\in\cal S$,
the random variable
$V_j^{sign}$ has approximately zero mean, unlike sums of votes $V_j$,
whereas
at support indices
$|V_i^{sign}| \approx V_i$
since support indices are unlikely to be sent to the fusion center with the
opposite sign of $\theta^*_i$.
In the simulation section
we illustrate the improved performance of a sign-based over a votes-based distributed scheme.

\section{Theoretical results}
\label{sec:theoretical_results}

In this section, we present a theoretical analysis for one of our schemes.
Specifically, both Theorems
\ref{teo:supprec_threshold_votes} and \ref{teo:supprec_threshold_votes_M_low} show that under suitable conditions, with high probability Algorithm \ref{alg:count_votes}  achieves exact support recovery with little communication per machine.
In Theorem \ref{teo:supprec_threshold_votes}, the number of machines is relatively large, and the communication per machine is linear in the sparsity  $K$. In Theorem
\ref{teo:supprec_threshold_votes_M_low} the number of machines $M$ is logarithmic in $d$, in which case the communication per machine is much higher, though still sublinear in $d$.
Both theorems are based on the Gaussian approximation in Lemma \ref{lemma:approx_cdf_dist_deblasso_estimator} and on probability bounds for binomial random variables. Their proofs appear in Appendix \ref{sec:proofs}.

To put our theorems in context, let us briefly review previous results on exact support recovery in the simpler (non-distributed) sparse linear regression setting.
A key quantity characterizing the ability of exact support recovery is the signal strength, defined as  $\theta_{\min}=\min_{i\in\mathcal{S}}|\theta_i^*|$. As proven by \cite{wainwright2009sharp}, under suitable conditions on the design matrix,
the lasso estimator based on $n$ samples and an appropriately chosen regularization parameter $\lambda_n$,
achieves exact
support recovery with high probability, provided that $\theta_{\min} \gtrsim \sqrt{(\ln d)/n}$.
The same rate $\theta_{\min} \gtrsim \sqrt{(\ln d)/n}$ is also sufficient for support recovery using a debiased lasso estimator (see, e.g., section 2.2 of \cite{javanmard2014confidence}).
In a distributed setting, \cite{lee2017communication} proved that with high probability, their scheme achieves exact support recovery when $\theta_{\min} \gtrsim \sqrt{(\ln d)/(n M)}$.
While this result matches the centralized setting, their scheme
requires each machine to send to the center its $d$-dimensional debiased lasso estimate, incurring
$O(d)$ communication per machine.
Hence, an interesting range for the signal strength, for the study of support recovery under communication constraints, is $\sqrt{\frac{\ln d}{nM}} \lesssim \theta_{\min} \lesssim \sqrt{\frac{\ln d}{n}}$. In this range, individual machines may be unable to exactly recover the support using the lasso or debiased lasso estimators.

To derive support recovery guarantees, we assume the smallest {nonzero} coefficient of
$\theta^*$ is sufficiently large, namely
$|\theta_i^*| \geq \theta_{\min}$ for all $i\in \mathcal S$ and some suitable $\theta_{\min}>0$. For our analysis below,
conditional on the design matrices
${\bf X}^1,\ldots,{\bf X}^M$ at the $M$ machines,
it will be convenient to make the following change of variables
from $\theta_{\min}$ to the (data-dependent) SNR parameter $r$,
\begin{align}
    \theta_{\min} =  \theta_{\min}(d,\sigma,r,n,c_{\Omega}) =
    \sigma
    \sqrt{2\frac{c_{\Omega}}{n} r\ln d},
    \label{e:theta_min}
\end{align}
where $c_{\Omega}$ is defined as
\begin{align}
    c_{\Omega} = \max_{i\in[d],m\in[M]} \left( \hat{\Omega}^m\hat{\Sigma}^{m}(\hat{\Omega}^m)^{\top} \right)_{ii}.
    \label{e:sigma2_max}
\end{align}
Recall from Eq. \eqref{e:standardized_debiased_lasso} that by Theorem \ref{teo:deblasso_asymp_dist}, $\sigma^2\left( n^{-1}\hat{\Omega}^m\hat{\Sigma}^{m}(\hat{\Omega}^m)^{\top} \right)_{ii}$ is the asymptotic variance of $\hat\theta^m_i$.
Hence,
$\sigma^2 c_{\Omega}/n$ is the largest {variance} of all $d$ coordinates of the $M$ debiased estimators computed by the $M$ machines.
In terms of the SNR parameter $r$, the range of interest is thus $\frac1M < r < 1$.

Recall that our scheme is based on thresholding the normalized debiased lasso estimators $\hat{\xi}_k^m$ of Eq. \eqref{e:standardized_debiased_lasso}.
We denote the corresponding normalized signal by
\begin{align}
    \vartheta_k^m = \frac{\sqrt{n}\theta_k^*}{\sigma \left( \hat{\Omega}^m \hat{\Sigma}^m (\hat{\Omega}^m)^{\top} \right)^{1/2}_{kk}},
    \quad \forall k\in[d] .   \label{e:vartheta_min_standardized_deblasso}
\end{align}
Lemma \ref{lemma:approx_cdf_dist_deblasso_estimator} states that under suitable conditions, $\hat{\xi}_k^m - \vartheta_k^m$
has approximately a standard Gaussian distribution. This property plays an important role in our theoretical results.
Eq. \eqref{e:approx_cdf_dist_deblasso_estimator} of Lemma \ref{lemma:approx_cdf_dist_deblasso_estimator} provides a bound on the error
between the CDF of $\hat{\xi}_k^m - \vartheta_k^m$ and that of a standard Gaussian.
For a threshold $\tau$, let  $\epsilon(\tau)$ be the largest of these error bounds over all
$d$ coordinates in all $M$ machines,
\begin{align}
\epsilon(\tau) = \max_{k\in[d], m\in[M]} \left\{ \frac{\delta_R {\phi\left( \tau - \vartheta_k^m \right)} }{\sigma\left( \hat{\Omega}^m \hat{\Sigma}^m (\hat{\Omega}^m)^{\top} \right)^{1/2}_{kk}}
+  2de^{-c_* n/K} + de^{-c n} + \frac{6}{d^2} \right\}.
\label{e:largest_approx_error_cdf_xihat_Gaussian}
\end{align}
Recall that $\delta_R$, defined in
Eq. \eqref{e:bound_bias_term_deltaR},
is an upper bound on the bias $\hat{\theta}_k^m - \theta^*_k$.
By Lemma 5.4 of \cite{van2014asymptotically}, if the row sparsity of $\Omega$ satisfies $K_{\Omega} = o(n/\ln d)$ and $\hat{\Omega}^m$ is computed with regularization $\lambda_{\Omega} \propto \sqrt{(\ln d)/n}$, then $\left( \hat{\Omega}^m \hat{\Sigma}^m (\hat{\Omega}^m)^{\top} \right)_{kk}
\geq
\Omega_{kk} + o_P(1) \geq C_{\max}^{-1} + o_P(1)$ when $\frac{\ln d}{n} \to 0$.
Hence, when $n$ and $d$ are large,
all terms on the {right hand} side of
Eq. \eqref{e:largest_approx_error_cdf_xihat_Gaussian} are small, and the Gaussian approximation is accurate.

To prove that our scheme recovers $\mathcal{S}$ with high probability, we  assume that:
\begin{enumerate}[({C}1)]
    \item The $n$ samples in each of the $M$ machines are i.i.d. from the model \eqref{e:linear_model} and the conditions of Theorem \ref{teo:deblasso_asymp_dist} all hold. Additionally, all machines use the same regularization parameters $\lambda$ and $\lambda_{\Omega}$ to compute the lasso \eqref{e:lasso_estimator} and debiased lasso \eqref{e:debiased_lasso} estimators, respectively.  \label{cond:asymp_dist_deblasso}

    \item $|\theta_i^*| \geq \theta_{\min}(d,\sigma,r,n,c_{\Omega})$ for all $i\in \mathcal{S}$, where $\theta_{\min}$ and $c_{\Omega}$ are defined in Eqs. \eqref{e:theta_min} and \eqref{e:sigma2_max}, respectively. \label{cond:theta_larger_thetamin}
\end{enumerate}

The following theorem provides
a recovery guarantee for Algorithm \ref{alg:count_votes},
where the sparsity $K$ is assumed to be known to the fusion center.

\begin{theorem}
Suppose Algorithm \ref{alg:count_votes}
is run with threshold $\tau=\sqrt{2\ln d}$.
Assume that $d$ is sufficiently large and that the SNR in Eq. \eqref{e:theta_min} satisfies
$$
\frac14 \frac{\ln^2(48\sqrt{\pi}\ln^{3/2} d)}{\ln^2(d)} < r<1.
$$
Additionally, assume conditions C\ref{cond:asymp_dist_deblasso} and C\ref{cond:theta_larger_thetamin} hold and the approximation error in Eq. \eqref{e:largest_approx_error_cdf_xihat_Gaussian} satisfies $\epsilon(\tau) \leq 1/d$.
Then, if the number of machines satisfies
\begin{align}
    \frac{8\ln d}{\frac{\sqrt{2\left(1-\sqrt{r}\right)^2\ln d}}{\sqrt{2\pi} \left\{ 2\left(1-\sqrt{r}\right)^2\ln d + 1 \right\}} - \epsilon(\tau) d^{\left(1-\sqrt{r}\right)^2}} d^{\left(1-\sqrt{r}\right)^2} \leq M \leq \frac{d}{3},
    \label{e:number_machines_thresh_maxgauss}
\end{align}
with probability at least $1-\frac{K+1}{d}$, Algorithm \ref{alg:count_votes} achieves exact support recovery. Additionally, the expected communication per machine is $O\left( K \ln d \right)$ bits.
\label{teo:supprec_threshold_votes}
\end{theorem}

Let us make a few remarks regarding this theorem. The upper bound $M<d/3$ is rather artificial and stems from the fact that in our proof we assume $M<d$. It is possible to derive support guarantees also for the case $M>d$, though this setting
seems to be unlikely in practice.
The lower bound on the number of machines is required to guarantee that with high probability, all support indices receive more votes than any non-support coordinate.
The lower bound on the SNR $r$ ensures that the lower bound on the number of machines in Eq. (\ref{e:number_machines_thresh_maxgauss}) is indeed smaller than $d/3$, so the range of possible values for $M$ is not empty.
A similar lower bound on $r$ appeared in \cite{amiraz2022distributed} after their Theorem 1.B.

Another important remark is that the threshold $\tau=\sqrt{2\ln d}$ in Theorem
\ref{teo:supprec_threshold_votes} is relatively high, so each machine sends only few indices to the center.
However, to guarantee support recovery,
this requires a relatively large number of machines $M=polylog(d) \cdot d^{(1-\sqrt{r})^2}$.
In Theorem \ref{teo:supprec_threshold_votes_M_low}, we give sufficient conditions to still achieve a high probability of exact support recovery when the number of machines is much smaller, of order only logarithmic in $d$. The price to pay is a higher communication per machine, which nonetheless is still sub-linear in $d$, namely much lower than the communication required to send the whole debiased lasso vector. For the next theorem, we assume that a lower bound on the SNR is known to all machines, which set a threshold that depends on it.

\begin{theorem}
Suppose Algorithm \ref{alg:count_votes}
is run with threshold $\tau=\sqrt{2r\ln d}$.
Assume that $d$ is sufficiently large and that the SNR in Eq. \eqref{e:theta_min} satisfies
$$
\frac{\ln(16\ln d)}{\ln d} < r<1.
$$
Additionally, assume conditions C\ref{cond:asymp_dist_deblasso} and C\ref{cond:theta_larger_thetamin} hold and the approximation error in Eq. \eqref{e:largest_approx_error_cdf_xihat_Gaussian} satisfies $\epsilon(\tau) < 1/(4d^r)$.
If the number of machines satisfies
\begin{align}
    \frac{16\ln d}{1 - {2\epsilon(\tau)}} \leq M \leq d^r,
    \label{e:number_machines_thresh_M_low}
\end{align}
then with probability at least $1-\frac{K+1}{d}$, Algorithm \ref{alg:count_votes} achieves exact support recovery,  with  expected communication per machine
$O\left( d^{1-r} \ln d \right)$ bits.
\label{teo:supprec_threshold_votes_M_low}
\end{theorem}

Beyond support recovery, another quantity of interest is the accuracy of the distributed estimator $\hat{\theta}$ of Eq. \eqref{e:theta_hat}.
The following corollary, proven in Appendix \ref{sec:proofs},
shows that once $\mathcal{S}$ is precisely recovered, $\hat{\theta}$ is close to the oracle least squares estimator $\hat\theta^{\mbox{\tiny LS}}$ computed with all data
in a single machine and with knowledge of the true support.
Consequently, $\hat\theta$ is also close to the true vector $
\theta^*$.

\begin{corollary}
Assume the conditions of Theorem \ref{teo:supprec_threshold_votes} hold. Let $N=nM$ denote the total sample size in all $M$ machines. If $M = O\left( \frac{N K}{( \max\{K,\ln N\})^2} \right)$, then
\begin{align*}
\|\hat{\theta} - \hat\theta^{\mbox{\tiny LS}}\|_2 = O_{P}\left( \frac{\sqrt{M} \max\{K , \ln N\}}{N} \right)
\quad\text{and}\quad
\|\hat{\theta} - \theta^*\|_2 = O_{P}\left( \sqrt{\frac{K}{N}} \right),
\end{align*}
as $d, N \to \infty$ and $\frac{\ln d}{N/M} \to 0$,
where $\hat{\theta}$ is defined in Eq. \eqref{e:theta_hat} and
$\hat\theta^{\mbox{\tiny LS}}$ is the least squares solution using all $N$ samples and with a known $\cal S$, as in Eq. \eqref{e:centralized_ols_estimator}, appended by zeros at all coordinates $j\notin\mathcal S$.
\label{corol:supprec_threshold_votes}
\end{corollary}

Corollary \ref{corol:supprec_threshold_votes} shows that in
a high dimensional sparse setting, for a sufficiently strong signal,
Algorithm \ref{alg:two_round_scheme} with a threshold $\tau=\sqrt{2\ln d}$ achieves the same error rate as the oracle estimator.
Let us put this result in a broader context.
If the support $\cal S$ were known, then each machine could have computed its least squares solution restricted to $
\cal S$ and send it to the center for averaging.
As discussed in  \cite{rosenblatt2016optimality},
in a general setting of $M$-estimators,
if the number of machines is not too large, averaging
is optimal and to leading order coincides with the centralized solution.
Yet, while being rate optimal, we note that averaging does lead to a loss of accuracy and is not as efficient as the oracle estimator, see \cite{dobriban2021distributed}.

As mentioned in Remark \ref{remark:two_round_estimator_Battey_etal_2018}, \cite{battey2018distributed} also proposed a two-round estimator that attains the optimal rate in Corollary \ref{corol:supprec_threshold_votes}, but requires each machine to send at least $d$ values to the fusion center.
In contrast, $\hat{\theta}$ is computed using a much lower communication cost. Similar results can also be established for Algorithm \ref{alg:two_round_scheme} under the conditions of Theorem \ref{teo:supprec_threshold_votes_M_low}.

\subsection{Comparison to other works}
\label{subsec:comparison_other_works}

Theorems \ref{teo:supprec_threshold_votes} and \ref{teo:supprec_threshold_votes_M_low} can be viewed as analogous to Theorems 2.A and 2.B of \cite{amiraz2022distributed}, who studied distributed estimation of a sparse vector under communication constraints.
A key difference is that in their setting, the $d$ coordinates of the estimated vector at each machine are all independent and unbiased.
This allowed them to analyze a top-$L$ scheme since the probability of sending any non-support index was the same and could be easily bounded by $L/d$. As such, their proofs do not directly apply to our setting.
In our case, the debiased lasso $\hat\theta^m$ still has a small bias term and its $d$ coordinates are in general correlated. This implies that the probabilities of sending non-support indices are not all identical.
In our analysis, we bypass this issue by analyzing instead a thresholding approach at step 5 of Algorithm \ref{alg:count_votes}, where each  $\hat\theta^m_i$ for $i\in[d]$ is compared  separately to a fixed threshold. This way, we do not need to account for the complex dependence among different coordinates of the debiased lasso.

\cite{barghi2021distributed}
considered a similar distributed scheme to estimate the support of $\theta^*$. They did not normalize the debiased lasso estimate at each machine, and more importantly the estimated support set
consists only of those indices that received at least $M/2$ votes.
The authors performed a theoretical analysis of their scheme, though various quantities are described only up to unspecified multiplicative constants.
We remark that both theoretically as well as empirically, the SNR must be quite high for support indices to obtain at least $M/2$ votes. In our work, we present explicit expressions for the minimum SNR in Eq. \eqref{e:theta_min}, sufficient for exact support recovery, requiring far fewer votes at the fusion center. The simulations in the next section illustrate
the advantages of our scheme as compared to that of \cite{barghi2021distributed}.

\section{Simulations}
\label{sec:simulations}

We present simulations that illustrate the performance of our proposed methods in comparison to other distributed schemes.
We focus on methods based on debiased lasso estimates, and specifically
consider the following five
distributed schemes to estimate $\theta^*$
and its support.
\begin{itemize}
    \item \texttt{thresh-votes}: Algorithm \ref{alg:count_votes} with a threshold of $\tau = \sqrt{2\ln d}$.
    \item \texttt{top-L-votes}: The top $L$ algorithm presented in section \ref{subsec:variations_supp_estim}. Each machine sends the indices of its top $L$ values of $|\hat{\xi}_i^m|$ to the fusion center.
    \item \texttt{top-L-signs}: Each machine sends both the indices and signs of the top $L$ values of $|\hat{\xi}_i^m|$. The center forms $\hat{\mathcal{S}}$ using sums of signs as in Eq. (\ref{e:sum_signs}).
    \item \texttt{BNM21}: The algorithm proposed by \cite{barghi2021distributed} with a threshold $\tau = \sqrt{2\ln d}$. It is similar to Algorithm \ref{alg:count_votes}, the difference is that $\hat{\mathcal{S}}$ consists of the indices $i\in[d]$ with {$V_i \geq M/2$}.
    \item \texttt{AvgDebLasso}: Based on \cite{lee2017communication}, each machine sends its debiased lasso estimate $\hat{\theta}^m$. The center computes $\hat{\theta}^{avg} = \frac{1}{M} \sum_{m=1}^M \hat{\theta}^m$ and estimates the support as the indices with the $K$ largest values  $|\hat{\theta}^{avg}_i|$, $i\in[d]$.
    \end{itemize}
In all these algorithms, each machine computes a debiased lasso estimator using its own data. The methods differ by the content and length of the messages sent to the fusion center and hence by the manner in which the fusion center estimates $\mathcal S$ and $\theta^*$.
For a fair comparison, we run all methods with the same
regularization parameters $\lambda_{\Omega} = 2\sqrt{(\ln d)/n}$ and $\lambda = 8\sqrt{(\ln d)/n}$ in each machine to compute the precision matrix $\hat{\Omega}^m$ and the lasso estimator $\tilde{\theta}^m$, respectively.

We performed simulations with both known sparsity, where the methods were run as described above, as well as unknown sparsity.
In the latter case, for \texttt{thresh-votes}, \texttt{top-L-votes} and \texttt{top-L-signs}, the center computes $\hat{\mathcal{S}}$ as the indices $i$ such that $V_i$ or $|V_i^{sign}|$ is larger than $2\ln d$, see Remark \ref{remark:unknown_sparsity} in Section \ref{sec:proposed_method};
For
\texttt{AvgDebLasso},
we set $\hat{\mathcal{S}}$ as the indices $i$ such that $|\hat{\theta}_i^{avg}| > 11\frac{\ln d}{n} $. This scaling is motivated by Theorem 16 of \cite{lee2017communication}, whereby in our simulation setting this term is larger than the term $O\left(\sqrt{\ln d/(nM)}\right)$ in their bound for
$\|\hat \theta^{avg}-\theta^*\|_\infty$. The factor 11 was manually tuned for good results.
\texttt{BNM21} is unchanged, as it does not require knowledge of $K$.

We evaluate the accuracy of
an estimated support set $\hat S$ by the F-measure,
\begin{align*}
    \text{F-measure} = \frac{2\cdot \mathrm{precision} \cdot \mathrm{recall}}{\mathrm{precision} + \mathrm{recall}},
\end{align*}
where $\mathrm{precision}=|\mathcal{S}\cap \hat{\mathcal{S}}|/|\hat{\mathcal{S}}|$ and $\mathrm{recall}=|\mathcal{S}\cap \hat{\mathcal{S}}|/K$.
An F-measure equal to one indicates that exact support recovery was achieved.

Given a support estimate $\hat{\mathcal{S}}$, the vector $\theta^*$ is estimated as follows.
For all four methods excluding
\texttt{AvgDebLasso}, we perform a second round and compute the estimator
$\hat{\theta}$ given by Eq. \eqref{e:theta_hat}.
 \texttt{AvgDebLasso} is a single round scheme. Its estimate of $\theta^*$
 consists of $\hat{\theta}^{avg}$ restricted to the indices $i\in\hat{\mathcal{S}}$.
The error of an estimate $\hat\theta$ is measured by its $\ell_2$-norm $\|\hat{\theta} - \theta^*\|_2$.
As a benchmark for the achievable accuracy, we also computed the oracle centralized estimator $\hat\theta^{\mbox{\tiny LS}}$ that knows the support $\mathcal{S}$ and estimates $\theta^*$ by least squares on the whole data.

We generated data as follows:
The design matrix ${\bf X}^m\in \re^{n\times d}$ in machine $m$ has $n$ rows i.i.d. $N(0,\Sigma)$, with $\Sigma_{i,j}=0.5^{|i-j|}$. We then computed for each machine its matrix $\hat\Omega^m$, and the quantity
$c_\Omega$ in Eq. \eqref{e:sigma2_max}.
Next, we generated a $K$ sparse vector $\theta^*\in\re^d$,
whose {nonzero} indices are sampled uniformly at random from $[d]$.
Its nonzero coefficients have random $\pm 1$ signs
and their magnitudes are chosen from $K$ equally spaced values $\{\theta_{\min},\ldots, 2\theta_{\min}\}$, where
\begin{align*}
    \theta_{\min} = \sqrt{\frac{2 c_{\Omega} \ln d}{n}}.
\end{align*}
These matrices and the vector $\theta^*$ are then kept fixed.
For a simulation with SNR parameter $r$, we set
$\sigma = 1/\sqrt{r}$. Finally, in each realization
we generated the response $Y^m\in\re^n$ according to the model \eqref{e:linear_model}.

\begin{figure}[t]
\centering
\begin{subfigure}{0.49\linewidth}
      \centering
\includegraphics[scale=.45]{fmeasure_SNR_n250_p5000_s5_M100_circulant.png}
      \caption{}
\end{subfigure}
\begin{subfigure}{0.49\linewidth}
      \centering
\includegraphics[scale=.45]{ell2error_SNR_n250_p5000_s5_M100_circulant.png}
      \caption{}
\end{subfigure}
\caption{
Results for known sparsity, averaged over 500 realizations, as a function of SNR in the range $r\in[1/M,1]$. (a)
F-measure; (b) $\ell_2$ error, on a log scale.
}
\label{f:fmeas_risk_vary_SNR}
\end{figure}

Our first simulation compared the performance of various schemes as a function of the SNR, with a known sparsity $K=5$. We fixed the dimension $d=5000$, the sample size in each machine $n=250$, and the number of machines $M=100$.
The top $L$ method was run with $L=K$ (see Sec. \ref{subsec:variations_supp_estim}).
Figure \ref{f:fmeas_risk_vary_SNR}a displays the F-measure of each method, averaged over 500 realizations.
As expected, at low SNR values, \texttt{AvgDebLasso} achieved the best performance in terms of support recovery.
However, for stronger signals with $r>0.4$, both \texttt{top-K-votes} and \texttt{thresh-votes}
achieved an F-measure of one, in accordance with our theoretical results regarding exact recovery.
In particular, at sufficiently high SNR, our methods estimate the support as accurately as \texttt{AvgDebLasso},
but with 2-3 orders of magnitude less communication.
The scheme of \texttt{BNM21} achieves good performance only at higher SNR.

Figure \ref{f:fmeas_risk_vary_SNR}b shows the errors $\|\hat \theta - \theta^*\|_2$, averaged over 500 realizations. At low SNR, {$1/M < r<0.1$}, \texttt{AvgDebLasso} has the
smallest error.
However, for $r>0.4$, \texttt{thresh-votes} and {\texttt{top-K-votes}} yield more accurate estimates.
Thus, Figure \ref{f:fmeas_risk_vary_SNR}b shows the benefits of an accurate support estimate followed by a distributed least squares in a second round.
Indeed, at these SNR levels, our methods exactly recover the support.  Consequently, the second round reduces to a distributed ordinary least squares restricted to the correct support set $\cal S$. In accordance with Corollary \ref{corol:supprec_threshold_votes},
Algorithm \ref{alg:count_votes} then has the same error rate as the oracle.

\begin{figure}[t]
\centering
\begin{subfigure}{0.49\linewidth}
      \centering
\includegraphics[scale=.45]{unknownK_fmeasure_SNR_n250_p5000_s5_M100_circulant.png}
      \caption{}
\end{subfigure}
\begin{subfigure}{0.49\linewidth}
      \centering
\includegraphics[scale=.45]{unknownK_ell2error_SNR_n250_p5000_s5_M100_circulant.png}
      \caption{}
\end{subfigure}
\caption{
Results for unknown sparsity, averaged over 500 realizations, as a function of SNR in the range $r\in[1/M,1]$. (a)
 F-measure; (b) $\ell_2$ error on a log scale.
  }
\label{f:fmeas_circulant_unknown_K}
\end{figure}

Next we present simulation results for unknown sparsity, as a function of the SNR in the range $r\in[\frac1M,1]$.
As seen in Figure \ref{f:fmeas_circulant_unknown_K},
throughout this SNR range,
\texttt{AvgDebLasso} with a threshold of $11 (\ln d )/n$ achieves an F-measure close to one and $\ell_2$ errors close to those of the oracle.
In contrast, \texttt{thresh-votes}
achieves accurate estimates only for $r>0.6$.
These results illustrate that even when sparsity is unknown, our schemes can accurately estimate the vector $\theta^*$ and its support, albeit with a higher SNR as compared to the case of known sparsity.

\begin{figure}[t]
\centering
\begin{subfigure}{0.49\linewidth}
      \centering
\includegraphics[scale=.45]{ell2_error_vary_n_p5000_s5_M100_circulant.png}
      \caption{}
\end{subfigure}
\begin{subfigure}{0.49\linewidth}
      \centering
\includegraphics[scale=.45]{ell2_error_vary_M_p5000_s5_n250_circulant.png}
      \caption{}
\end{subfigure}
\caption{
$\ell_2$ error vs. sample size $n$ (a) and vs. number of machines $M$ (b), both on a log-log scale, at an SNR of $r=0.4$. Values are averaged over 1000 realizations.
}
\label{f:ell2_error_varying_n_varying_M}
\end{figure}

Finally, we present simulation results as a function of number of samples $n\in[n_{\min},n_{\max}] =[100,400]$ with $M=100$ machines,
and as a function of number of machines
$M\in[M_{\min},M_{\max}]=[40,160]$
with $n=250$ samples per machine.
Initially, for each machine $m\in[M_{\max}]$, we generated its full data {matrix ${\bf X}^{m}$} with $n_{\max}$ number of samples. We then computed the corresponding matrix $\hat\Omega^m$ and the quantity
$c_\Omega$ of Eq. \eqref{e:sigma2_max}. We then generated the sparse vector $\theta^*$ as described above.
To save on run-time, in simulations with a smaller number of samples $n<n_{\max}$ we nonetheless used the
decorrelation matrices $\Omega^m$ that correspond to $n_{\max}$ samples. This can be viewed as a semi-supervised setting, as also mentioned in \cite{javanmard2018debiasing}.
We fixed $r=0.4$ and $d=5000$, and the sparsity $K=5$ was known to the center.
Figure \ref{f:ell2_error_varying_n_varying_M} shows the resulting $\ell_2$ errors averaged over 1000 simulations.
In this simulation,  \texttt{top-K-votes} and \texttt{thresh-votes} are close to the centralized least squares oracle since their support estimates are accurate, as can be seen in Figure \ref{f:fmeas_risk_vary_SNR} for $r=0.4$.
Both {plots in Figure \ref{f:ell2_error_varying_n_varying_M}} show a linear dependence on a log-log scale with a slope of approximately $-1/2$, namely that the resulting errors decay as $1/\sqrt{n}$
and $1/\sqrt{M}$, respectively.
This is in agreement with our theoretical result in Corollary \ref{corol:supprec_threshold_votes}.

\begin{figure}[t]
\centering
\begin{subfigure}{0.49\linewidth}
      \centering
\includegraphics[scale=.45]{fmeasure_SNR_n250_p5000_s25_M100_circulant.png}
      \caption{}
\end{subfigure}
\begin{subfigure}{0.49\linewidth}
      \centering
\includegraphics[scale=.45]{ell2error_SNR_n250_p5000_s25_M100_circulant.png}
      \caption{}
\end{subfigure}
\caption{Results for schemes using sums of signs and sums of votes, averaged over 500 realizations, as a function of $r\!\in\![\frac1M, 0.1]$. (a) F-measure; (b) $\log_{10}$ of $\ell_2$ error.}
\label{f:fmeas_signs_vs_votes}
\end{figure}

\subsection{The advantages of sending signs}

We now illustrate the advantages of using sums of signs instead of sums of votes, in terms of both support recovery and parameter estimation.
In this simulation, we fixed $n=250$, $d=5000$, $M=100$ and $K=25$.
The results in Figure \ref{f:fmeas_signs_vs_votes} show that using sums of signs is more accurate than sums of votes for low SNR values in the range $r\in[1/M, 0.1]$.
Figure \ref{f:fmeas_signs_vs_votes} also shows that using $L=5K$ instead of $L=K$ significantly improves the accuracy of the support estimator, at the expense of increasing the communication. This
illustrates the potential trade-offs between the accuracy of support estimation and communication.

\section{Summary and Discussion}
    \label{sec:summary}

The development and analysis of
distributed statistical inference schemes having low communication
are important contemporary problems.
Given its simplicity and ubiquity, the sparse linear regression model has attracted significant
attention in the literature.
Most previous inference schemes for this model require communication per machine of at least $O(d)$ bits.
In this work we proved
theoretically and showed via simulations that, under suitable conditions,
accurate distributed inference for sparse linear regression is possible
with a much lower communication per machine.

Over the past years, several authors studied distributed statistical inference under communication constraints. Specifically, for sparse linear regression,
\cite{braverman2016communication} proved that without a lower bound on the SNR, to obtain a risk comparable to that of the minimax lower bound, a communication of at least $\Omega(M \min(n,d)/\log d)$ bits is required.
\citet{acharya2019distributed} proved that, under certain conditions, rate optimal estimates of a linear regression model can be computed using total communication sublinear in the dimension.
However, as they mention in their appendix B.3, a precise characterization of the ability to recover the support with sublinear communication in $d$ and its dependency on other parameters such as SNR and the number of machines is still an open problem.
In our theoretical results,
we presented explicit expressions for the minimal SNR at which our scheme is guaranteed to achieve exact recovery with high probability and with sublinear communication.
While we did not address the open problem of tight lower bounds,
our results highlight
the potential tradeoffs between SNR, communication and number of machines.

We believe that using more refined techniques, our theoretical analysis can be extended and improved.
For example, since the $d$ coordinates of a debiased lasso estimator are  correlated, sharp concentration bounds for dependent variables, like those of \cite{lopes2022sharp}, could improve our  analysis and extend it to other schemes such as top-$L$.
In our analysis, we focused on a setting where both the noise and covariates have Gaussian distribution. \citet{lee2017communication} and \citet{battey2018distributed}, for example, considered sub-Gaussian distributions for these terms. Our results can be adapted for this case, but a careful control of the various constants in probability bounds is needed to derive explicit expressions.

Finally, our low-communication schemes could also be applied to other problems,
such as sparse $M$-estimators, sparse covariance estimation and distributed estimation of jointly sparse signals. We leave these for future research.


text:
\subsection*{Exercises}}

\newcommand{\Rightarrow}{\Rightarrow}
\DeclareMathOperator{\Fin}{Fin}
\newenvironment{enalpha}{\begin{enumerate}
\renewcommand\labelenumi{(\alph{enumi})}}
{\end{enumerate}}

\newcounter{my_enumerate_counter}
\newcommand{\setcounter{my_enumerate_counter}{\value{enumi}}}{\setcounter{my_enumerate_counter}{\value{enumi}}}
\newcommand{\setcounter{enumi}{\value{my_enumerate_counter}}}{\setcounter{enumi}{\value{my_enumerate_counter}}}

\newenvironment{myenum}{\begin{enumerate}\setcounter{enumi}{\value{my_enumerate_counter}}}
{\setcounter{my_enumerate_counter}{\value{enumi}}\end{enumerate}}

\DeclareMathOperator{\Coll}{Coll}
\DeclareMathOperator{\dom}{dom}
\DeclareMathOperator{\ran}{ran}
\DeclareMathOperator{\Sym}{Sym}
\DeclareMathOperator{\cl}{cl}
\newcommand{\Co}[1]{\Coll(\aleph_1,<#1)}
\newcommand{\Coq}[1]{\Coll(#1,\delta)}
\newcommand{\Coqq}[2]{\Coll(#1,#2)}
\DeclareMathOperator{\crit}{cr}
\def\restriction{\restriction}
\newcommand{\mathfrak M}{\mathfrak M}
\newcommand{\mathcal M^*}{\mathcal M^*}
\newcommand{\FileName}[1]{\thanks{Filename: {\tt #1}}}
\newenvironment{enumerate_plus}
{\leftmargini=50pt
\begin{enumerate}}{\end{enumerate}}

\newcommand{\label}{\label}
\renewcommand{\land}{\wedge}
\newcommand{{\mathbf A}}{{\mathbf A}}
\newcommand{{\mathbf X}}{{\mathbf X}}
\newcommand{{\mathbf D}}{{\mathbf D}}
\newcommand{{\mathbf Y}}{{\mathbf Y}}
\newcommand{{\mathbf H}}{{\mathbf H}}
\newcommand{{\mathbf f}}{{\mathbf f}}
\newcommand{{\mathbf g}}{{\mathbf g}}
\newcommand{{\mathbf {\bar R}}}{{\mathbf {\bar R}}}
\newcommand{{\mathbf R}}{{\mathbf R}}
\newcommand{{\mathbf T}}{{\mathbf T}}
\newcommand{{\mathbf E}}{{\mathbf E}}
\newcommand{{\mathcal P}}{{\mathcal P}}
\newcommand{{\mathcal O}}{{\mathcal O}}
\newcommand{\sqsubseteq}{\sqsubseteq}
\newcommand{\underline\cup}{\underline\cup}

\newcommand{{\mathcal U}}{{\mathcal U}}
\newcommand{\mathcal V}{\mathcal V}
\newcommand{{}^\frown}{{}^\frown}
\newcommand{\mathbb Z}{\mathbb Z}
\newcommand{\bbZ/2\bbZ}{\mathbb Z/2\mathbb Z}

\newcommand{\mathbf a}{\mathbf a}
\newcommand{\mathbf b}{\mathbf b}
\newcommand{\mathbf c}{\mathbf c}
\newcommand{\mathbf P}{\mathbf P}
\newcommand{{\mathbf Q}}{{\mathbf Q}}
\newcommand{{\mathbf S}}{{\mathbf S}}
\newcommand{{\mathbf Q}}{{\mathbf Q}}
\newcommand{{\mathbf R}}{{\mathbf R}}
\newcommand{\text{even}}{\text{even}}
\newcommand{\text{odd}}{\text{odd}}
\newcommand{\mathbb B}{\mathbb B}
\newcommand{\mathbb M}{\mathbb M}
\newcommand{\bbN^{\uparrow\bbN}}{\mathbb N^{\uparrow\mathbb N}}
\newcommand{(\cU(1))^{\bbN}}{({\mathcal U}(1))^{\mathbb N}}
\newcommand{\mathcal I}{\mathcal I}
\newcommand{\mathbf K}{\mathbf K}
\newcommand{\mathfrak b}{\mathfrak b}
\newcommand{\mathfrak d}{\mathfrak d}
\newcommand{\mathbf U}{\mathbf U}
\newcommand{\mathbf Z}{\mathbf Z}
\DeclareMathOperator{\Ad}{Ad}
\newcommand{\left(\begin{matrix} 1 & 0 \\ 0 & -1 \end{matrix}\right)}{\left(\begin{matrix} 1 & 0 \\ 0 & -1 \end{matrix}\right)}

\newcommand{\mathcal F}{\mathcal F}
\newcommand{\mathcal G}{\mathcal G}
\DeclareMathOperator{\Aut}{Aut}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator{\graph}{graph}
\DeclareMathOperator{\proj}{proj}
\newcommand{\kappa}{\kappa}
\newcommand{\kappa'}{\kappa'}
\newcommand{\!\upharpoonright\!}{\!\upharpoonright\!}

\newcommand{\renewcommand{\labelenumi}{{\rm (\roman{enumi})}}\renewcommand{\itemsep}{0pt}}{\renewcommand{\labelenumi}{{\rm (\roman{enumi})}}\renewcommand{\itemsep}{0pt}}
\DeclareMathOperator{\Cu}{Cu}

\setlength{\topmargin}{0.0in}
\setlength{\textheight}{8in}

\setlength{\evensidemargin}{-0.3in}
\setlength{\oddsidemargin}{-0.3in}
\setlength{\textwidth}{7.2in}

\address{Department of Mathematics, University of Toronto, Toronto, Ontario M5S 2E4, Canada}
\email{elliott@math.toronto.edu}

\address{Department of Mathematics and Statistics, York University, 4700 Keele St., Toronto, Ontario, M3J 1P3, Canada}
\email{ifarah@mathstat.yorku.ca}

\address{Department of Mathematics, University of Houston, Houston, TX 77204-3476}
\email{vern@math.uh.edu}

\address{Department of Mathematics, Statistics, and Computer Science (M/C 249), University of Illinois at Chicago, 851 S. Morgan St., Chicago, IL 60607-7045, USA}
\email{rosendal.math@gmail.com}

\address{Department of Mathematics, Purdue University, 150 N. University St., West Lafayette IN 47906, USA}
\email{atoms@purdue.edu}

\address{Department of Mathematical Sciences, University of Copenhagen, Universitetsparken 5, 2100 Copenhagen, Denmark}
\email{asgert@math.ku.dk}

\title{The isomorphism relation for separable C*-algebras}

\author[Elliott et al.]{George A. Elliott, Ilijas Farah, Vern Paulsen, Christian Rosendal, Andrew S.  Toms and Asger T\"ornquist}

\date{\today}
\begin{document}

\begin{abstract}
We prove that the isomorphism relation for separable C$^*$-algebras, and also  the relations of complete and $n$-isometry for operator spaces and systems, are Borel reducible to the orbit equivalence relation of a Polish group action on a standard Borel space.
\end{abstract}

\maketitle

\section{Introduction}

The problem of classifying a collection of objects up to some notion of isomorphism can usually be couched as the study of an analytic equivalence relation on a standard Borel space parameterizing the said objects.  Such relations admit a notion of comparison, {\it Borel reducibility}, which allows one to assign a degree of complexity to the classification problem.  If $X$ and $Y$ are standard Borel spaces admitting equivalence relations $E$ and $F$ respectively then we say that $E$ is {\it Borel reducible} to $F$, written $E \leq_B F$, if there is a Borel map $\Theta:X \to Y$ such that
\[
xEy \iff \Theta(x)F\Theta(y).
\]
In other words, $\Theta$ carries equivalence classes to equivalence classes injectively.  We view $E$ as being ``less complicated" than $F$. There are some particularly prominent degrees of complexity in this theory which serve as landmarks for classification problems in general.  For instance, a relation $E$ is {\it classifiable by countable structures} (CCS) if it is Borel reducible to the isomorphism relation for countable graphs.  Classification problems in functional analysis (our interest here) tend not to be CCS, but may nevertheless be ``not too complicated" in that they are Borel reducible to the orbit equivalence relation of a Borel action of a Polish group on a standard Borel space;  this property is known as being {\it below a group action}.

The connections between Borel reducibility and operator algebras have received considerable attention lately.  Using Hjorth's theory of turbulence developed in \cite{Hj:Book}, it has been shown that several classes of operator algebras
are not CCS. This applies to von Neumann factors of every type \cite{sato09a},
ITPFI$_2$ factors \cite{SaTo:Turbulence},
separable simple unital nuclear C$^*$-algebras \cite{FaToTo:Turbulence},
and spaces of irreducible representations of non-type I C$^*$-algebras, \cite{KerLiPi:Turbulence,Fa:Dichotomy}. On the other hand, Elliott's K-theoretic classification of AF algebras together with the Borel computability of K-theory \cite{FaToTo:Descriptive} show that AF algebras are CCS.

We are interested in classifying separable C$^*$-algebras, operator spaces, and operator systems.  The classification problem for nuclear simple separable C$^*$-algebras was studied in \cite{FaToTo:Turbulence}, where the isomorphism relation was shown to be below a group action.  Establishing this upper bound was rather involved;  it required Borel versions of Kirchberg's $\mathcal{O}_2$ embedding and absorption theorems.  It of course invited the question of whether isomorphism of \emph{all} separable C$^*$-algebras is below a group action.  Here we give a surprisingly simple proof of the following result.

\begin{theorem}\label{T1}
Each of the following equivalence relations
is below a group action:
\begin{enumerate}
\item  The isomorphism relation for separable C$^*$\! -\! algebras.
\item Complete isometry, as well as $n$-isometry for any $n$,
 of separable operator spaces.
\item Complete isometry, as well as $n$-isometry for any $n$,
 of separable operator systems.
\item Unital, complete isometry of separable operator systems.
\item  Unital, complete order isomorphism of separable operator systems.
 \end{enumerate}
\end{theorem}
\noindent Part (1)  of course improves on one of the main results
of \cite{FaToTo:Turbulence}. Parts (2--5) give a partial answer to a question stated by Ed Effros during the January 2012 meeting ``Set theory and C*-algebras" at the American Institute of Mathematics. Theorem \ref{T1} follows from a more general result proved in \S\ref{s.PS}. We show that $\Iso(\mathbb U)$, the isometry group of the Urysohn
metric space $(\mathbb U,\delta)$, plays a role for the isometry relation  of what we term {\it Polish structures in a countable signature} (see \cite{BYBHU,FaHaSh:Model2}), analogous to the role played by the Polish group $S_\infty$ of all permutations of $\mathbb{N}$ for the isomorphism relation of countable structures.

\section{Polish structures}\label{s.PS}

The main ingredients for proving Theorem \ref{T1} are the notion of a \emph{Polish structure}, which we introduce now, and Theorem \ref{t.bga} below.

\begin{definition}\label{d.polstruct}
Let $\mathcal L=(l_1,\ldots)$ be a finite or infinite sequence in $\mathbb N$. A \emph{Polish $\mathcal L$-structure}\footnote{We consciously avoid the term ``metric structure'' since this is already used for a slightly different notion in continuous logic.
} is a triple
\[
\mathbb X=(X,d^{\mathbb X},(F_n^{\mathbb X})),
\]
 where $(X,d)$ is a separable complete metric space, called the \emph{domain}, and $F_n^{\mathbb X}\subseteq X^{l_n}$ is closed in the product topology. We think of $F_n^{\mathbb X}$ as a relation on $X$ and call $l_n$ the \emph{arity} of $F_n^{\mathbb X}$. The sequence $\mathcal L$ is called the \emph{signature} or \emph{arity sequence} of the structure $\mathbb X$. (Below we will suppress the superscript $\mathbb X$ whenever possible.)

Two Polish structures $\mathbb X$ and $\mathbb Y$ with the same signature $\mathcal L$ are said to be \emph{isometrically isomorphic} if there is an isometric bijection $h: X\to Y$ such that for all $n$ and $(x_1,\ldots,x_{l_n})\in X^{l_n}$ we have
$$
(x_1,\ldots,x_{l_n})\in F_n^{\mathbb X}\iff (h(x_1),\ldots, h(x_{l_n}))\in F_n^{\mathbb Y}.
$$
\end{definition}

Let $(Y,d)$ be a metric space, let $\mathcal L=(l_1,\ldots)$, and define $M(\mathcal L,Y,d)\subseteq (F(Y)\setminus\{\varnothing\})\times \prod_{n} F(Y^{l_n})$ by
$$
M(\mathcal L,Y,d)=\{(X,(F_n)): F_n\subseteq X^{l_n}\}.
$$
This set is Borel (by \cite[12.11]{Ke:Classical}), and to each $\mathbb X=(X,(F_n))\in M(\mathcal L,Y,d)$ we have the Polish $\mathcal L$-structure $(X,d\!\upharpoonright\! X,(F_n))$, which we also denote by $\mathbb X$. For $\mathbb X, \mathbb Y\in M(\mathcal L,Y,d)$, write $\mathbb X\simeq^{M(\mathcal L,Y,d)}\mathbb Y$ if $\mathbb X$ and $\mathbb Y$ are isomorphic Polish $\mathcal L$-structures.

When $(Y,d)$ is the Urysohn metric space $\mathbb U$, we will usually write $M(\mathcal L)$ rather than $M(\mathcal L,\mathbb U,\delta)$. This is motivated by the fact that every separable metric space can be isometrically embedded into $\mathbb U$, and so every Polish $\mathcal L$-structure is isomorphic to some $\mathbb X\in M(\mathcal L)$. Thus $M(\mathcal L)$ provides a parametrization of all Polish $\mathcal L$-structures as a standard Borel space. It is not hard to see that that this moreover is a good standard Borel parametrization in the sense of \cite{FaToTo:Turbulence}.

The group $\Iso(\mathbb U)$ acts diagonally on $\mathbb U^{l_i}$ for each $i$, and so it acts naturally on $M(\mathcal L)$ by $\sigma\cdot (X,(F_n))=(\sigma\cdot X,(\sigma\cdot F_n))$. This induces the orbit equivalence relation
$$
\mathcal X\equiv^{M(\mathcal L)}\mathcal Y\iff (\exists\sigma\in\Iso(\mathbb U))\,\sigma\cdot\mathcal X=\mathcal Y.
$$
The following theorem was inspired by   \cite[\S 2]{GaKe}.

\begin{theorem}\label{t.bga}
We have $\simeq^{M(\mathcal L)}\,\leq_B\,\equiv^{M(\mathcal L)}$, and so $\simeq^{M(\mathcal L)}$ is below a group action.
\end{theorem}

To prove this, we first need a uniformly Borel version of injective universality and homogeneity of $\mathbb U$. This will also be used several times later. For a countable set $A$, define
$$
\mathcal M_A=\{d\in \mathbb R^{A\times A}: d\text{ is a metric on }A\}.
$$
This is easily seen to form a $G_\delta$ subset of $\mathbb R^{A\times A}$, whence $\mathcal M_A$ is Polish in the subspace topology. When $A=\mathbb N$ we let $\mathcal M=\mathcal M_\mathbb N$.

\begin{lemma}\label{l.unifhom}
There are Borel functions $\vartheta_n:\mathcal M\to\mathbb U$, $n\in\mathbb Z$, such that for all $d,d'\in\mathcal M$:
\begin{enumerate}
\item the set $\{\vartheta_n(d):n\in\mathbb Z\}$ is dense in $\mathbb U$;
\item $\delta(\vartheta_n(d),\vartheta_m(d))=d(n,m)$ for all $n,m\in\mathbb N$;
\item any isometric bijection
$$
\sigma:\overline{\{\vartheta_n(d):n\in\mathbb N\}}\to\overline{\{\vartheta_n(d'):n\in\mathbb N\}}
$$ extends to an isometric automorphism of $\mathbb U$.
\end{enumerate}
\end{lemma}

\begin{proof}
We adapt the Katetov construction of the Urysohn metric space; see \cite{Kat} and \cite{GaKe}. Fix an enumeration $(A_n)_{n\in\mathbb N}$ of all finite non-empty subsets of $\mathbb Z$, and let $*$ be a point \emph{not} in $\mathbb Z$. For $d\in \mathcal M$, call a metric $\rho\in\mathcal M_\mathbb Z$ a \emph{Katetov extension} of $d$ if for all $n\in\mathbb N$, all $\varepsilon>0$, and every metric $\tilde d$ on $A_n\cup\{*\}$ satisfying
\begin{enumerate}
\item $\tilde d(x,y)=d(x,y)$ for all $x,y\in A_n$, and
\item $\tilde d(x,*)\in{\mathbb Q}_+$ for all $x\in A_n$
\end{enumerate}
there is is some $i\in\mathbb Z\setminus A_n$ such that $|\rho(i,x)-\tilde d(*,x)|<\varepsilon$ for all $x\in A_n$.
It follows directly from the definition that the set
$$
U=\{(d,\rho)\in\mathcal M\times \mathcal M_\mathbb Z:\rho\text{ is a Katetov extension of } d\}
$$
is $G_\delta$. The group $\Sym(\mathbb Z\setminus\mathbb N)$ of all permutations of $\mathbb Z\setminus\mathbb N$ acts naturally on $\mathcal M_\mathbb Z$, and the sections $U_d=\{\rho:(d,\rho)\in U\}$ are invariant under this action. The Katetov property guarantees that the $\Sym(\mathbb Z\setminus\mathbb N)$-orbits in $U_d$ are dense in $U_d$. It follows from Theorem \ref{t.hsp} (see appendix) that there is a Borel $\psi:\mathcal M\to \mathcal M_\mathbb Z$ such that $(d,\psi(d))\in U$ for all $d\in\mathcal M$. A standard approximate intertwining/back-and-forth argument (e.g., \cite{Ell:Towards}) shows that for any $(d,\rho)\in U$, the completion of $(\mathbb Z,\rho)$ is is isometrically isomorphic to $(\mathbb U,\delta)$, and a Borel coding of this argument analogous to \cite[Theorem 7.6]{FaToTo:Turbulence} shows that there are Borel maps $\vartheta_n:\mathcal M\to\mathbb U$, $n\in\mathbb Z$, satisfying (1) and (2) above. Finally, another approximate intertwining/back-and-forth argument can be used to establish (3).
\end{proof}

\begin{proof}[Proof of Theorem \ref{t.bga}]
We will define a Borel function $M(\mathcal L)\to M(\mathcal L):\mathbb X\mapsto\mathbb X'$ such that $\mathbb X\simeq^{M(\mathcal L)}\mathbb Y$ if and only if  $\mathbb X'\equiv^{M(\mathcal L)} \mathbb Y'$, and $\mathbb X\simeq^{M(\mathcal L)}\mathbb X'$. For $\mathbb X,\mathbb Y\in M(\mathcal L)$ with domains $X$ and $Y$ finite, the homogeneity of $\mathbb U$ immediately gives that $\mathbb X\simeq^{M(\mathcal L)} \mathbb Y$ if and only if  $\mathbb X\equiv^{M(\mathcal L)}\mathbb Y$, so here we may simply take $\mathbb X=\mathbb X'$.

For the case when $\mathbb X$ has infinite domain, recall that by the Kuratowski--Ryll-Nardzewski theorem we may find Borel functions $\psi_n^{k}: F(\mathbb U^k)\to\mathbb U^k$ such that $\{\psi_n^k(F):n\in\mathbb N\}$ is dense in $F$ whenever $F\neq\varnothing$. From the $\psi_n^k$ we can define a sequence of Borel functions $\tilde\psi_n:F(\mathbb U)\to\mathbb U$ such that whenever $\mathbb X=(X,(F_k))\in M(\mathcal L)$ and $X$ is infinite, the sequence $\tilde\psi_n(X)$ gives an injective enumeration of a dense subset of $X$ and for all $k$, $\psi^{l_k}_n(F_k)\subseteq \{\tilde\psi_n(X):n\in\omega\}^{l_k}$. Define an element of $d_X\in\mathcal M$ by $d_X(n,m)=\delta(\tilde\psi_n(X),\tilde\psi_m(x))$, and let $\gamma_n(X)=\vartheta_n(d_X)$, $n\in\mathbb Z$, where the $\vartheta_n$ are provided by Lemma \ref{l.unifhom}. Then the $\gamma_n$ are Borel, and if we define $\mathbb X'=(X',F_n')$ by letting
$$
X'=\cl(\{\gamma_n(X):n\in\mathbb N\})
$$
and
\begin{align*}
F_n'=&\cl(\{(\gamma_{i_1}(X),\ldots,\gamma_{i_{l_n}}(X)):\\
&(\exists (i_j))(\tilde\psi_{i_1}(X),\ldots,\tilde\psi_{i_{l_n}}(X))\in F_n\}),
\end{align*}
then $\mathbb X'\simeq^{M(\mathcal L)}\mathbb X$. That the map $\mathbb X\mapsto \mathbb X'$ is Borel follows from the Kuratowski--Ryll-Nardzewski theorem, and that $\mathbb X\simeq^{M(\mathcal L)}\mathbb Y$ precisely when $\mathbb X'\equiv^{M(\mathcal L)}\mathbb Y'$ is immediate from Lemma \ref{l.unifhom}.(3).
\end{proof}

\section{Proof of Theorem \ref{T1}}

\subsection{C*-algebras}\label{S.C*} A C*-algebra $A$ can be cast as a Polish structure with domain $A$ as follows: We have the relations $F_0=\{0\}\subseteq A$, $F_{+},F_{\cdot}\subseteq A^3$, which are the graphs of the functions $(a,b)\mapsto a+b$ and $(a,b)\mapsto ab$, respectively, $F_{*}\subseteq A^2$ which is the graph of the map $a\mapsto a^{*}$, and for each $q\in{\mathbb Q}(i)={\mathbb Q}+i{\mathbb Q}$ the relations
$$
F_q=\{(a,b)\in A^2: b=qa\}.
$$
The signature is $\mathcal L_{C^*}=(1,3,3,2,2,\ldots)$. It is easy to check directly that the set $M_{C^*}\subseteq M(\mathcal L_{C^*})$ which corresponds to C*-algebras is Borel. (See also Lemma \ref{L.2} below, where a more general statement is proven.) Let $\simeq^{M_{C^*}}$ denote the isomorphism relation in $M_{C^*}$. It can easily be shown that $M_{C^*}$ provides a standard Borel parametrization of the class of separable C*-algebras, in the sense of \cite[Definition~2.1]{FaToTo:Turbulence}). Clearly $\simeq^{M_{C^*}}$ is the restriction of $\simeq^{M(\mathcal L_{C^*})}$ to $M_{C^*}$, and so by Theorem \ref{t.bga} above $\simeq^{M_{C^*}}$ is below a group action.

To finish the proof of Theorem \ref{T1}.(1) we need only show that the parametrization $M_{C^*}$ is weakly equivalent to the parametrizations given in \cite{FaToTo:Turbulence}, in the sense of \cite[Definition 2.1]{FaToTo:Turbulence}. For this, recall from \cite[\S 2.4]{FaToTo:Turbulence} that $\hat \Xi$ is the space of countable normed ${\mathbb Q}(i)$-$*$-algebras with domain $\mathbb N$ which satisfy the C*-axiom. For each $\xi\in\hat\Xi$, let $C^*(\xi)$ be the C*-algebra obtained from completing $\xi$ and extending operations.  Let $\vartheta_n:\mathcal M\to\mathbb U$ be as in Lemma \ref{l.unifhom}. To each $\xi\in\hat\Xi$ we have the associated $d_\xi\in\mathcal M$ defined by $d_\xi(n,m)=\|n-m\|_\xi$, and the map $\xi\mapsto d_\xi$ is clearly Borel. It is then straightforward to define $\mathbb X^\xi=(X^\xi,F_0^\xi,F^{\xi}_+,F^\xi_{\cdot},F^\xi_{*},(F^\xi_{q})_{q\in{\mathbb Q}(i)})\in M_{C^*}$ directly from $d_\xi$ and the $\vartheta_n$ so that $\mathbb X^\xi$ is isomorphic to $C^*(\xi)$, letting $F_0=\{\vartheta_{0_\xi}(d_\xi)\}$,
$$
F^\xi_+=\{(\vartheta_n(d_\xi),\vartheta_m(d_\xi),\vartheta_k(d_\xi)): n+_{\xi}m=k\},
$$
and defining $F^\xi_{\cdot}$, $F^\xi_{*}$, $(F^\xi_q)_{q\in{\mathbb Q}(i)}$ analogously. The map $\xi\mapsto\mathbb X^\xi$ is then Borel by the Kuratowski--Ryll-Nardzewski theorem.

For the converse direction, recall from \cite[\S 2.4]{FaToTo:Turbulence} that $\Xi\subseteq\mathbb R^\mathbb N$ consists of all real sequences $\eta$ such that for some C*-algebra $A$ and some $y=(y_n)_{n\in\mathbb N}$ which is dense in $A$, we have that $\eta_n=\|\mathfrak p_n(y)\|_A$, where $(\mathfrak p_n)_{n\in\mathbb N}$ enumerates the non-commutative ${\mathbb Q}(i)$-$*$-polynomials. Let $f_n:F(\mathbb U)\setminus\{\varnothing\}\to\mathbb U$, $n\in \mathbb N$, be Borel functions provided by the Kuratowski--Ryll-Nardzewski theorem such that $\{f_n(X): n\in \mathbb N\}$ is a dense subset of $X$, for all $X\in F(\mathbb U)\setminus\{\varnothing\}$. For $\mathbb X\in M_{C^*}$, let $\eta^X_n=\delta(\mathfrak p_n((f_n(X)_{n\in\mathbb N}),0_\mathbb X)$ where $F^\mathbb X_0=\{0_\mathbb X\}$, and $\mathfrak p_n((f_n(X))$ is the evaluation of $\mathfrak p_n$ at $y=(f_n(X))$ in the C*-algebra coded by $\mathbb X$. It is easily seen that $X\mapsto \eta^X$ is Borel and $\mathbb X$ and $\eta^\mathbb X$ code isomorphic C*-algebras.

\subsection{Banach spaces and Banach algebras}
A separable (real or complex) Banach space $E$ can be cast as a Polish structure $(E,F_0,F_+,(F_q)_{q\in\mathbb K})$, where $\mathbb K={\mathbb Q}$ or $\mathbb K={\mathbb Q}(i)$, in the signature $\mathcal L_{B}=(1,3,2,2,\ldots)$ in analogy with the above parametrization of C*-algebras (omitting multiplication and involution). The subset $M_B(\mathbb U)$ of $M(\mathcal L_{B},\mathbb U)$ corresponding to Banach spaces is easily seen to be Borel, too. By an argument similar to that in \S\ref{S.C*}, one sees that this parameterization is equivalent to the standard parameterization of Banach spaces as closed subspaces of $C([0,1])$ (see e.g., \cite{FeLouRo} or \cite{Mell:Computing}).

In a similar vein, letting $\mathcal L_{BA}=(1,3,3,2,2,\ldots)$ (which happens to coincide with $\mathcal L_{C^*}$), we can parametrize separable Banach algebras, as well as involutive separable Banach algebras, by appropriate Borel subsets of $M(\mathcal L_{BA})$.

In all cases, Theorem \ref{t.bga} provides that the isometric
isomorphism relation is below a group action.

\subsection{Operator spaces and operator systems}\label{S.OS} To handle operator spaces and operator systems we will need a framework slightly more general than that of Polish structures. Formally, this could be handled by introducing \emph{multi-sorted} Polish structures, though here we make due with an ad hoc approach.

For metric spaces $(X,d_X)$ and $(Y,d_Y)$, let $UC(X,Y)$ denote the set of uniformly continuous functions $f:X\to Y$. Identifying each $f\in UC(X,Y)$ with $\graph(f)\subseteq X\times Y$, it is easily checked that $UC(X,Y)$ forms a Borel subset of $F(X\times Y)$.

Let $\mathbb U_n$, $n\in\mathbb N$, be a sequence of disjoint copies of the Urysohn metric space, with $\mathbb U_1=\mathbb U$. Identify $M_n(\mathbb U)$ with $\mathbb U^{n^2}$, and give this the supremum metric.
We define $\mathbb {OS}$ to be the set of sequences $\mathbb E=(E_n,f_n)_{n\in\mathbb N}$ such that:
\begin{enumerate}
\item $(E_n,f_n)\in M_B(\mathbb U_n)\times UC(M_n(\mathbb U),\mathbb U_n)$;
\item $f_n(M_n(E_1))= \dom(E_n)$.
\item $f_n\!\upharpoonright\! M_n(E_1)$ is linear;
\item for all $A,B\in M_n(\mathbb C)$ and all $X\in M_n(E_1)$ we have
$$
\|f_n(AXB)\|_{E_n}\leq\|A\|\|f_n(X)\|_{E_n}\|B\|;
$$
\item for all $m,n\in\mathbb N$ and all $X\in M_n(E_1)$, $Y\in M_m(E_1)$
$$
\|X\oplus Y\|_{E_{n+m}}=\max\{\|X\|_{E_n},\|Y\|_{E_m}\|\}.
$$
\end{enumerate}
Once again, $\mathbb {OS}$ forms a Borel set. Each $\mathbb E\in\mathbb {OS}$ codes a $L^\infty$-matrix-normed space, whence by Ruan's theorem (e.g. \cite[Theorem 13.4]{Paulsen}), $\mathbb {OS}$ parametrizes the class of separable operator spaces. Defining $\mathbb E\simeq^{\mathbb {OS}}\mathbb F$ if and only if $\mathbb E$ and $\mathbb F$ are completely isometrically isomorphic, it is easy to check that $\simeq^{\mathbb {OS}}$ is analytic, and so $\mathbb {OS}$ provides a good standard Borel parametrization of the separable operator spaces. Also, write $\mathbb E\simeq^{\mathbb {OS}}_n\mathbb F$ if and only $\mathbb E$ and $\mathbb F$ are $n$-isometric.

The group $\Iso(\mathbb U)\times\Iso(\mathbb U_n)$ acts in a Borel way on\break
 $UC(M_n(\mathbb U),\mathbb U_n)$ by
$$
((\sigma,\tau)\cdot f)(x_{ij})=\tau(f(\sigma^{-1}(x_{ij})).
$$
Thus we obtain a Borel action of $\prod_n \Iso(\mathbb U_n)$ on $\mathbb {OS}$ by
$$
(\sigma_n)\cdot(E_n,f_n)=(\sigma_n\cdot E_n,(\sigma_1,\sigma_n)\cdot f_n),
$$
and we write $\mathbb E\equiv^{\mathbb {OS}}\mathbb F$ if and only there is $(\sigma_n)$ such that $(\sigma_n)\cdot \mathbb E=\mathbb F$. Arguing as we did in the proof of Theorem \ref{t.bga}, we obtain that $\simeq^{\mathbb {OS}}\leq_B\equiv^{\mathbb {OS}}$, which proves that $\simeq^{\mathbb {OS}}$ is below a group action.

If we only consider the action of $\prod_{j\leq n} \Iso(\mathbb U_n)$ on $\mathbb {OS}$, and denote by $\equiv^{\mathbb {OS}}_n$ the induced orbit equivalence relation, the argument from Theorem \ref{t.bga} gives that $\simeq^{\mathbb {OS}}_n\,\leq_B\,\equiv^{\mathbb {OS}}_n$, thus showing $\simeq^{\mathbb {OS}}_n$ is below a group action. This establishes Theorem \ref{T1}.(2).

\medskip

A standard parameterization of operator systems is obtained by adding the adjoint operation to the structures parametrized by $\mathbb {OS}$. Then an analogous argument establishes Theorem \ref{T1}.(3).

\medskip

By adding  a constant for the multiplicative unit to the language one  sees that
the relation of unital complete isometry between operator systems is below a group action.
Finally, unital maps between operator systems are completely isometric if and only
if they are completely positive (see \cite[Proposition~3.6]{Paulsen}) and (4) and (5) of
Theorem~\ref{T1} follow.

\subsection{A remark about models of the logic for metric structures} \label{S.Met}
The isomorphism relation of countable structures, in the sense of \cite{Hj:Book}, is given by a continuous $S_\infty$-action.
The group $\Iso(\mathbb U)$ plays an analogous role for separable models of logic for metric structures. Such models consists of an underlying Polish space $X$, countably many functions $f_n\colon X^n\to X$, and countably many functions  $g_n\colon X^n\to \mathbb R$ (\emph{relations}). (Assuming  there are infinitely many functions and that the $n$'th function has $X^n$ as its domain is clearly not a loss of generality.) These functions are required to have a prescribed modulus of uniform continuity (see \cite{BYBHU}),
but we shall ignore this since it is not important for our present purposes.
To a fixed model $\mathcal X=(X,f_n, g_n: n\in \mathbb N)$, associate a Polish structure
$\mathbb X=(X,(F_n)_{n\in \mathbb N})$ where $F_n$, $n\in \mathbb N$, enumerate
the graphs of the $f_n$, as well as all the sets $\{\bar x\in X^n: g_n(x)\geq q_n\}$ where $(q_n)_{n\in\mathbb N}$ is a fixed enumeration of the rationals. This map is Borel (between the appropriate spaces) and $\mathcal X$ is isomorphic to $\mathcal Y$ if and only if $\mathbb X\simeq\mathbb Y$.
 Also, the map that sends an element of $\Met$ to a metric structure is  Borel.

The following is proved by a straightforward recursion analogous to the case of classical logic.

\begin{lemma} \label{L.2}
If $T$ is a theory in a countable language in the logic of metric structures then
the set of all $\mathbb X\in \Met$ that code a model of $T$ is
Borel. \qed
\end{lemma}

Combining this with Theorem \ref{t.bga}, we obtain:

\begin{theorem} If $T$ is a theory of the logic of metric structures in a separable
language, then the isometry relation of models of $T$
is Borel-reducible to an orbit equivalence relation of $\Iso(\mathbb U)$. \qed
\end{theorem}

\section{Concluding remarks}

(1) Aaron Tikuisis has pointed out that the above approach, using Polish structures, also can be used to provide a new proof that the isomorphism relation for von Neumann algebras with separable predual is below a group action. Previously, in \cite{sato09a}, it was shown that isomorphism of separably acting von Neumann algebras is below an action of the unitary group of $\ell_2$, using a completely different line of argument. It is at present not known if the latter provides a sharper upper bound on complexity than the former.

(2) In \cite{FaToTo:Turbulence} it was proved that separable unital AI-algebras
are not classifiable by countable structures, and that their isomorphism relation is below an action of $\Aut(\mathcal O_2)$. However, we do not at present know if the complexity
of the classification problem increases as one passes from nuclear to exact C*-algebras, or from exact to arbitrary C*-algebras. In particular, we do not know the answer to:

\begin{question}
Is the isomorphism relation for separable C*-algebras strictly more complicated, as measured by $\leq_B$, than that of nuclear separable C*-algebras?
\end{question}

Even if we restrict to nuclear \emph{simple} separable C*-algebras, we don't know the answer to this.

We may also ask if the upper bound on complexity provided by Theorem \ref{T1}.(1) is actually optimal. By \cite{GaKe}, the Borel actions of $\Iso(\mathbb U)$ realize the maximal complexity of equivalence relations induced by Polish group actions.

\begin{question}
Is the isomorphism relation for separable C*-algebras maximal among equivalence relations induced by a Polish group action?
\end{question}

We also don't know whether unital  $n$-order isomorphism of separable operator systems
is below a group action for~$n\in \mathbb N$.

\section{Acknowledgment}

The main result of this paper was proved in the aftermath of the BIRS
workshop on applications of descriptive set theory to functional analysis. We would
like to thank the staff at BIRS for providing an inspiring atmosphere.

G.A.E. and I.F.  partially supported by NSERC.
C.R. was partially supported by NSF grants 1201295 and 0901405.
A.T.  supported in part by grant no. 10-082689 from the Danish Council for Independent Research.

\section*{Appendix}

We prove the generalization of the homogeneous selection principle, \cite[Lemma 6.2 and 6.3]{FaToTo:Turbulence}, which is used to prove Lemma \ref{l.unifhom} above. For a subset $A\subseteq X\times Y$ of a Cartesian product and $x\in X$, we define  $A_x=\{y\in Y: (x,y)\in A\}$. A function $f:\proj_X(A)\to Y$ is a \emph{uniformization} of $A$ if $(x,f(x))\in A$ for all $x\in \proj_X(A)$.

\begin{appthm}\label{t.hsp}
Let $X, Y$ be Polish spaces, $d_Y$ a complete compatible metric on $Y$, and let $A\subseteq X\times Y$ be a $G_\delta$ set.

(1) If for some (any) sequence $(y_n)_{n\in\mathbb N}$ dense in $Y$ the set
$$
R=\{(x,n,\varepsilon)\in X\times\mathbb N\times{\mathbb Q}_+: (\exists y\in A_x) d_Y(y,y_n)<\varepsilon\}
$$
is Borel, then $\proj_X(A)$ is Borel and $A$ admits a Borel uniformization.

(2) If $G$ is a Polish group acting continuously on $Y$ by Borel automorphisms such that $A_x$ is $G$-invariant for all $x\in X$, and every $G$-orbit of a point $y\in A_x$ is dense in $A_x$, then $R$ defined in (1) is Borel, and so $A$ admits a Borel uniformization.
\end{appthm}
\begin{proof}
(1) We may assume that $d_Y$ is bounded by $\frac 1 2$. Clearly $\proj_X(A)=\{x\in X: (\exists n) (x,n,1)\in R\}$, so this set is Borel. To construct the uniformization, fix open sets $U_n\subseteq X\times Y$ such that $\bigcap_{n\in\mathbb N} U_n=A$, and fix an enumeration $(q_n)_{n\in\mathbb N}$ of ${\mathbb Q}_+$. Let $B(y,\varepsilon)$ denote the open $d_Y$-ball of radius $\varepsilon$ around $y\in Y$. We recursively define Borel maps $x\mapsto n_i(x)\in\mathbb N$ and $x\mapsto \varepsilon_i(x)\in{\mathbb Q}_+$, $i\in\mathbb N$, on $\proj_X(A)$, satisfying:
\begin{enumerate}
\item $n_i(x)=0$, $\varepsilon_0(x)=1$;
\item $\varepsilon_i(x)\leq \frac 1 {2^i}$;
\item $B(y_{n_i(x)},\varepsilon_i(x))\subseteq (U_i)_x$;
\item $(x,n_i(x),\varepsilon_i(x))\in R$;
\item $\overline{B(y_{n_{i+1}(x)},\varepsilon_{i+1})}\subseteq B(y_{n_i(x)},\varepsilon_i(x))$.
\end{enumerate}
If this can be done then $f(x)=\lim_{i\to\infty} y_{n_i(x)}$ is the desired Borel uniformization. Suppose $n_i(x)$ and $\varepsilon_i(x)$ have been defined for $i\leq k$. Let
$$
z\in B(y_{n_k(x)},\varepsilon_k(x))\cap A_x
$$
and let $0<\delta\leq 2^{k+1}$ be such that $2\delta+d_Y(z,y_{n_k(x)})<\varepsilon_k(x)$ and $B(z,2\delta)\subseteq U_{k+1}$. If $d(y_n,z)<\delta$, then $\overline{B(y_n,\delta)}\subseteq B(y_{n_k(x)},\varepsilon_k(x))$, $B(y_n,\delta)\cap A_x\neq\varnothing$, and $B(y_n,\delta)\subseteq U_{k+1}$. This shows that (1)--(5) above can be satisfied, and so we can define $n_{k+1}(x)$ and $\varepsilon_{k+1}(x)=q_n$, where $n_{k+1}(x)$ and $n$ are least possible satisfying (1)--(5) above. Since each requirement (1)--(5) are Borel, the maps $x\mapsto n_{k+1}(x)$ and $x\mapsto \varepsilon_{k+1}(x)$ are Borel, as required.

(2) It is clear from the definition that $R$ is analytic. Since
$$
R=\{(x,n,\varepsilon): (\forall y\in A_x)(\exists g\in G_0) d(g\cdot y,y_n)<\varepsilon\}
$$
for a countable dense $G_0\subseteq G$
it is also coanalytic, whence $R$ is Borel.
\end{proof}

\bibliographystyle{amsplain}

redpajama_set_name: RedPajamaArXiv