\section{Introduction}
\label{sec:intro}
Graph statistics is useful for finding meaningful connection patterns in network data, and 
subgraph counting is known as a fundamental task in graph analysis. 
For example, a triangle is a cycle of size three, 
and a $k$-star consists of a central node connected to $k$ other nodes.  
These subgraphs 
can be used to calculate 
a 
clustering coefficient ($=\frac{3 \times \text{\#triangles}}{\text{\#2-stars}}$). 
In a social graph, the clustering coefficient measures 
the tendency of nodes (users) to form a cluster with each other. 
It also represents the average probability that a friend's friend is also a friend \cite{Newman_PRL09}. 
Therefore, the clustering coefficient is useful for analyzing the effectiveness of friend suggestions. 
Another example of the subgraph is a 4-cycle, 
a cycle of size four. 
The 4-cycle count is useful for measuring the clustering ability in bipartite graphs (e.g., 
online dating networks, mentor-student networks \cite{Kutty_WWW14}) 
where a triangle never appears \cite{Lind_PRE05,Robins_CMOT04,Sanei-Mehri_CIKM19}. 
Figure~\ref{fig:subgraphs} shows examples of triangles, 2-stars, and 4-cycles. 
Although these subgraphs are important for analyzing the connection patterns or clustering tendencies, their exact numbers can leak sensitive edges (friendships) \cite{Imola_USENIX21}. 

\begin{figure}[t]
  \centering
  \includegraphics[width=0.65\linewidth]{fig/subgraphs.pdf}
  \vspace{-2mm}
  \caption{Examples of subgraph counts.
  }
  \label{fig:subgraphs}
\end{figure}

DP (Differential Privacy) \cite{Dwork_ICALP06,DP} -- the gold standard of privacy notions -- has been widely used to strongly protect edges in graph data \cite{Day_SIGMOD16,Ding_TKDE21,Hay_ICDM09,Imola_USENIX21,Imola_USENIX22,Karwa_PVLDB11,Kasiviswanathan_TCC13,Qin_CCS17,Sun_CCS19,Ye_ICDE20,Ye_TKDE21}. 
In particular, recent studies 
\cite{Imola_USENIX21,Imola_USENIX22,Qin_CCS17,Ye_ICDE20,Ye_TKDE21} 
have applied LDP (Local DP) \cite{Kasiviswanathan_FOCS08} to 
graph data. 
In the graph LDP model, each user obfuscates her neighbor list (friends list) by herself and sends the obfuscated neighbor list to a data collector. 
Then, the data collector estimates graph statistics, such as subgraph counts. 
Compared to central DP where a central server has personal data of all users (i.e., the entire graph), LDP does not have a risk that all personal data are leaked from the server by 
cyberattacks 
\cite{Henriquez_breach2021} or insider attacks \cite{Kohen_insider_threats}. 
Moreover, LDP can be applied to decentralized social networks \cite{Paul_CN14,Salve_CSR18} (e.g., diaspora* \cite{Diaspora}, Mastodon \cite{Mastodon}) where no server can access the entire graph; e.g., 
the entire graph is distributed across many servers, or no server has any original edges. 
It is reported in \cite{Imola_USENIX21} that $k$-star counts can be accurately estimated in this model. 

However, it is much more challenging to accurately count more complicated subgraphs such as triangles and 4-cycles under LDP. 
The root cause of this is its local property -- a user cannot see edges between others. 
For example, user $v_1$ cannot count triangles or 4-cycles including $v_1$, as she cannot see edges between others, e.g., $(v_2,v_3)$, $(v_2,v_4)$, and $(v_3,v_4)$. 
Therefore, the existing algorithms \cite{Imola_USENIX21,Imola_USENIX22,Ye_ICDE20,Ye_TKDE21} 
obfuscate 
each bit of the neighbor list 
rather than the subgraph count by 
the RR (Randomized Response) \cite{Warner_JASA65}, which randomly flips 0/1. 
As a result, their algorithms suffer from extremely large estimation errors because it makes all edges noisy. 
Some studies \cite{Imola_USENIX21,Imola_USENIX22} 
significantly improve the accuracy by introducing 
an additional round of interaction between users and the data collector. 
However, multi-round interaction may be impractical in many applications, as it requires a lot of user effort and synchronization; 
in  \cite{Imola_USENIX21,Imola_USENIX22}, every user must respond twice, and the data collector must wait for responses from all users 
in each round. 

In this work, we focus on a \textbf{one-round} of interaction between users and the data collector and propose accurate subgraph counting algorithms by introducing a recently studied privacy model: the \textit{shuffle model} \cite{Erlingsson_SODA19,Feldman_FOCS21}. 
In the shuffle model, each user sends her (encrypted) obfuscated data to an intermediate server called the shuffler. 
Then, the shuffler randomly shuffles the obfuscated data of all users and sends the shuffled data to the data collector (who decrypts them). 
The shuffling amplifies DP guarantees of the obfuscated data under the assumption that the shuffler and the data collector do not collude with each other. 
Specifically, it is known that DP strongly protects user privacy when a parameter (a.k.a. privacy budget) $\epsilon$ is small, e.g., $\epsilon \leq 1$ \cite{DP_Li}. 
The shuffling significantly reduces $\epsilon$ and therefore significantly improves utility at the same value of $\epsilon$. 
To date, the shuffle model has been successfully applied to tabular data \cite{Meehan_ICLR22,Wang_PVLDB20} and 
gradients \cite{Girgis_AISTATS21,Liu_AAAI21} in federated learning. 
We 
apply the shuffle model to graph data to accurately count subgraphs within one round. 

The main challenge in subgraph counting in the shuffle model is that each user's neighbor list is \textit{high-dimensional data}, i.e., $n$-dim binary string where $n$ is the number of users. 
Consequently, applying the RR to each bit of the neighbor list, as in the existing work \cite{Imola_USENIX21,Imola_USENIX22,Ye_ICDE20,Ye_TKDE21}, results in an extremely large privacy budget $\epsilon$ even after applying the shuffling (see Section~\ref{sub:technical} for more details). 

We address this issue by introducing a new, basic technique called \textit{wedge shuffling}. 
Our wedge shuffle technique obfuscates a wedge (2-hop path) between a specific user-pair. 
For example, in Figure~\ref{fig:subgraphs}, users $v_1$ and $v_4$ have a wedge between $v_2$ and $v_3$, whereas $v_5, \ldots, v_8$ do not. 
We obfuscate such wedge information by the RR and shuffle them. 
Because the wedge information is \textit{one-dimensional binary data}, it can be sent with small noise and small $\epsilon$. 
In addition, the wedge is the main component of several subgraphs, such as triangles, 4-cycles, and 3-hop paths \cite{Sun_CCS19}. 
Since the wedge has little noise, we can accurately count these subgraphs based on wedge shuffling. 

We apply wedge shuffling to triangle and 4-cycle counting tasks with several additional techniques. 
For triangles, we first 
propose an algorithm that counts triangles involving the user-pair at the endpoints of the wedges by locally sending an edge between the user-pair to the data collector. 
Then we propose an algorithm to count triangles in the entire graph by sampling independent user-pairs, which share no common users. 
We also propose a technique to reduce the variance of the estimate by ignoring sparse user-pairs, where either of the two users has a very small degree. 
For 4-cycles, we propose an algorithm to calculate an unbiased estimate of the 4-cycle count from that of the wedge count via bias correction. 

We provide upper bounds on the estimation error for our triangle and 4-cycles counting algorithms. 
Through comprehensive evaluation, we show that 
our algorithms 
accurately estimate these subgraph counts within one round under the shuffle model. 



\smallskip
\noindent{\textbf{Our Contributions.}}~~Our contributions are as follows: 
\begin{itemize}
    \item We propose a wedge shuffle technique to enable privacy amplification 
   
    of graph data. 
    To our knowledge, we are the first to apply the shuffle model to 
   
    graph data. 
    \item We propose one-round triangle and 4-cycle counting algorithms 
   
    based on our wedge shuffle technique. 
   
    For triangles, we propose three additional techniques: sending local edges, sampling independent user-pairs, and variance reduction by ignoring sparse user-pairs. 
    For 4-cycles, we propose a bias correction technique. 
    We show upper bounds on the estimation error for each algorithm. 
    \item We evaluate our algorithms using two real graph datasets. 
    Our experimental results show that our one-round shuffle algorithms 
    significantly outperform one-round local algorithms in terms of accuracy 
   
    and achieve a small estimation error (relative error $\ll 1$) with a reasonable privacy budget, e.g., smaller than $1$ in edge DP. 
\end{itemize}
In Appendix~\ref{sec:cluster}, we show that our triangle algorithm is also useful for accurately estimating the clustering coefficient within one round. 
We can use our algorithms to analyze the clustering tendency or the effectiveness of friend suggestions in decentralized social networks by introducing a shuffler. 
We implemented our algorithms in C/C++. 
Our code is available on GitHub \cite{SubgraphShuffle}. 
The proofs of all statements in the main body are given in Appendices~\ref{sec:triangle_proof} and \ref{sec:4cycle_proof}. 

\section{Related Work}
\label{sec:related}
\noindent{\textbf{Non-private Subgraph Counting.}}~~Subgraph counting has been extensively studied in a non-private setting (see \cite{Ribeiro_CS21} for a recent survey). 
Examples of subgraphs include triangles \cite{Bera_PODS20,Eden_FOCS15,Kolountzakis_IM12,Wu_TKDE16},  4-cycles \cite{Bera_STACS17,Kallaugher_PODS19,Manjunath_ESA11,McGregor_PODS20}, $k$-stars \cite{Aliakbarpour_Alg18,Gonen_DM11}, and 
$k$-hop paths \cite{Bjorklund_ICALP19,Kartun-Giles_TWC18}. 

Here, the main challenge is to reduce the computational time of counting these subgraphs in large-scale graph data. 
One of the simplest approaches is edge sampling \cite{Bera_PODS20,Eden_FOCS15,Wu_TKDE16}, which randomly 
samples edges in a graph. 
Edge sampling outperforms other sampling methods 
(e.g., node sampling, triangle sampling) \cite{Wu_TKDE16} and is also adopted in \cite{Imola_USENIX22} 
for private triangle counting. 

Although our triangle algorithm also samples user-pairs, ours is different from edge sampling in two ways. 
First, our algorithm does not sample an edge but samples a pair of users who may or may not be a friend. 
Second, our algorithm samples user-pairs that share no common users to avoid the increase of the privacy budget $\epsilon$ as well as to reduce the time complexity 
(see Section~\ref{sec:triangle} for details). 

\smallskip
\noindent{\textbf{Private Subgraph Counting.}}~~Differentially private subgraph counting has been widely studied, and the previous work assumes either the central \cite{Ding_TKDE21,Karwa_PVLDB11,Kasiviswanathan_TCC13} or local \cite{Imola_USENIX21,Imola_USENIX22,Sun_CCS19,Ye_ICDE20,Ye_TKDE21} models. 
The central model assumes a centralized social network and has a data breach issue, as explained in Section~\ref{sec:intro}. 

Subgraph counting in the local model has recently attracted attention. 
Sun \textit{et al.} \cite{Sun_CCS19} propose subgraph counting algorithms 
assuming 
that each user knows all friends' friends. 
However, this assumption does not hold in many social networks; e.g., Facebook users can change their settings so that anyone cannot see their friend lists. 
Therefore, we make a minimal assumption -- each user knows only her friends. 


In this setting, recent studies propose triangle \cite{Imola_USENIX21,Imola_USENIX22,Ye_ICDE20,Ye_TKDE21} and $k$-star \cite{Imola_USENIX21} counting algorithms. 
For $k$-stars, Imola \textit{et al.} \cite{Imola_USENIX21} propose a one-round algorithm that is order optimal and show that it provides a very small estimation error. 
For triangles, they
propose a one-round algorithm that applies the RR to 
each bit of the neighbor list 
and then calculates an unbiased estimate of triangles from the noisy graph. 
We call this algorithm \textsf{RR}$_{\triangle}${}. 
Imola \textit{et al.} \cite{Imola_USENIX22} show that \textsf{RR}$_{\triangle}${} provides a much smaller estimation error than the one-round triangle algorithms in \cite{Ye_ICDE20,Ye_TKDE21}. 
In \cite{Imola_USENIX22}, they also reduce the time complexity of \textsf{RR}$_{\triangle}${} 
by using the ARR (Asymmetric RR), which samples each 1 (edge) after applying the RR. 
We call this algorithm \textsf{ARR}$_{\triangle}${}. 
In this paper, we use \textsf{RR}$_{\triangle}${} and \textsf{ARR}$_{\triangle}${} as baselines in triangle counting. 
For 4-cycles, there is no existing algorithm under LDP, to our knowledge. 
Thus, we compare our shuffle algorithm with its local version, which does not shuffle the obfuscated data. 

For triangles, Imola \textit{et al.} also propose a two-round local algorithm in \cite{Imola_USENIX21} and significantly reduce its download cost in \cite{Imola_USENIX22}. 
Although we focus on one-round algorithms, 
we show in 
Appendix~\ref{sec:two-round} 
that our one-round algorithm is comparable to the two-round algorithm in \cite{Imola_USENIX22}, which requires a lot of user effort and synchronization, in terms of accuracy. 



\smallskip
\noindent{\textbf{Shuffle Model.}}~~The privacy amplification by shuffling has been recently studied in \cite{Balle_CRYPTO19,Cheu_EUROCRYPT19,Erlingsson_SODA19,Feldman_FOCS21}. 
Among them, the privacy amplification bound by Feldman \textit{et al.} \cite{Feldman_FOCS21} is the state-of-the-art -- it provides a smaller $\epsilon$ than other bounds, such as \cite{Balle_CRYPTO19,Cheu_EUROCRYPT19,Erlingsson_SODA19}. 
Girgis \textit{et al.} \cite{Girgis_CCS21} show a better bound than the bound in \cite{Feldman_FOCS21} when used with composition. 
However, the bound in \cite{Feldman_FOCS21} outperforms the bound in \cite{Girgis_CCS21} when used without composition (which is the case with this work). 
Therefore, we use the bound in \cite{Feldman_FOCS21}. 

The shuffle model has been applied to tabular data \cite{Meehan_ICLR22,Wang_PVLDB20} and gradients \cite{Girgis_AISTATS21,Liu_AAAI21} in federated learning. 
Meehan \textit{et al.} \cite{Meehan_ICLR22} construct a graph from public auxiliary information and determine a permutation of obfuscated data using the graph to reduce re-identification risks. 
Liew \textit{et al.} \cite{Liew_SIGMOD22} propose network shuffling, which shuffles obfuscated data via random walks on a graph. 
Note that both \cite{Meehan_ICLR22} and \cite{Liew_SIGMOD22} use graph data to shuffle another type of data. 
To our knowledge, our work is the first to shuffle graph data itself. 

\section{Preliminaries}
\label{sec:preliminaries}
In this section, we describe some preliminaries for our work. 
Section~\ref{sub:notations} defines the basic notations used in this paper. 
Sections~\ref{sub:privacy} and \ref{sub:shuffle} introduce DP on graphs and the shuffle model, respectively. 
Section~\ref{sub:utility} explains utility metrics. 

\subsection{Notations}
\label{sub:notations}
Let $\mathbb{R}$, $\mathbb{R}_{\ge0}$, $\mathbb{N}$, and $\mathbb{Z}_{\ge0}$ be the sets of real numbers, non-negative real numbers, natural numbers, and non-negative integers, respectively. 
For $a\in\mathbb{N}$, let $[a]$ be the set of natural numbers that do not exceed $a$, i.e., $[a] = \{1, 2, \ldots, a\}$. 

We consider an undirected social graph $G=(V,E)$, where $V$ represents a set of nodes (users) and $E \subseteq V \times V$ represents a set of edges (friendships). 
Let $n\in\mathbb{N}$ be the number of nodes in $V$, and $v_i \in V$ be the $i$-th node, i.e., $V=\{v_1,\ldots,v_n\}$. 
Let $I_{-(i,j)}$ be the set of indices of users other than $v_i$ and $v_j$, i.e., $I_{-(i,j)} = [n]\setminus\{i,j\}$. 
Let $d_i \in \mathbb{Z}_{\ge0}$ be a degree of $v_i$, $d_{avg} \in \mathbb{R}_{\ge0}$ be the average degree of $G$, and $d_{max} \in \mathbb{N}$ be the maximum degree of $G$. 
In most real graphs, $d_{avg} \ll d_{max} \ll n$ holds. 
We denote a set of graphs with $n$ nodes by $\mathcal{G}$. 
Let $f^\triangle: \mathcal{G} \rightarrow \mathbb{Z}_{\ge0}$ and $f^\square: \mathcal{G} \rightarrow \mathbb{Z}_{\ge0}$ be triangle and 4-cycle count functions, respectively. 
The triangle count function takes $G \in \mathcal{G}$ as input and outputs the number $f^\triangle(G)$ of triangles in $G$, 
whereas the 4-cycle count function takes $G$ as input and outputs the number $f^\square(G)$ of 4-cycles. 

Let $\mathbf{A}=(a_{i,j}) \in \{0,1\}^{n \times n}$ be 
an adjacency matrix corresponding to $G$. 
If $(v_i,v_j) \in E$, then $a_{i,j} = 1$; otherwise, $a_{i,j} = 0$. 
We call $a_{i,j}$ an \textit{edge indicator}. 
Let $\mathbf{a}_i \in \{0,1\}^n$ be a neighbor list of user $v_i$, i.e., the $i$-th row of $\mathbf{A}$. 
We summarize the basic notations in Table~\ref{tab:notations} of Appendix~\ref{sec:notation_table}. 

\subsection{Differential Privacy}
\label{sub:privacy}
\noindent{\textbf{DP and LDP.}}~~We use differential privacy, and more specifically $(\epsilon,\delta)$-DP \cite{DP}, as a privacy metric: 

\begin{definition} [$(\epsilon,\delta)$-DP \cite{DP}] \label{def:DP} 
Let $n \in \mathbb{N}$ be the number of users. 
Let $\epsilon \in \mathbb{R}_{\ge0}$ and $\delta \in [0,1]$. 
Let $\mathcal{X}$ be the set of input data for each user. 
A randomized algorithm $\mathcal{M}$ with domain $\mathcal{X}^n$ 
provides \emph{$(\epsilon,\delta)$-DP} if for any neighboring databases $D,D' \in \mathcal{X}^n$ that differ in a single user's data and any 
$S \subseteq \mathrm{Range}(\mathcal{M})$, 
\begin{align*}
\Pr[\mathcal{M}(D) \in S] \leq e^\epsilon \Pr[\mathcal{M}(D') \in S] + \delta.
\end{align*}
\end{definition}
$(\epsilon,\delta)$-DP guarantees that two neighboring datasets $D$ and $D'$ are almost equally likely when $\epsilon$ and $\delta$ are close to $0$. 
The parameter $\epsilon$ is called the privacy budget. 
It is well known that $\epsilon \leq 1$ is acceptable and $\epsilon \geq 5$ is unsuitable in many practical scenarios \cite{DP_Li}. 
In addition, the parameter $\delta$ needs to be much smaller than $\frac{1}{n}$ \cite{Barber_arXiv14,DP}. 

LDP \cite{Kasiviswanathan_FOCS08} is a special case of DP where $n=1$. 
In this case, a randomized algorithm is called a \textit{local randomizer}. 
We denote the local randomizer by $\mathcal{R}$ to distinguish it from the randomized algorithm $\mathcal{M}$ in the central model. 
Formally, LDP is defined as follows: 
\begin{definition} [$\epsilon$-LDP \cite{Kasiviswanathan_FOCS08}] \label{def:LDP} 
Let $\epsilon \in \mathbb{R}_{\ge0}$. 
Let $\mathcal{X}$ be the set of input data for each user. 
A local randomizer $\mathcal{R}$ with domain $\mathcal{X}$ 
provides \emph{$\epsilon$-LDP} if for any $x,x' \in \mathcal{X}$ and any $S \subseteq \mathrm{Range}(\mathcal{R})$, 
\begin{align}
\Pr[\mathcal{R}(x) \in S] \leq e^\epsilon \Pr[\mathcal{R}(x') \in S].
\label{eq:LDP}
\end{align}
\end{definition}


\smallskip
\noindent{\textbf{Randomized Response.}}~~We use Warner's RR (Randomized Response) \cite{Warner_JASA65} to provide 
LDP. 
Given $\epsilon \in \mathbb{R}_{\ge0}$, Warner's RR $\mathcal{R}_{\epsilon}^W: \{0,1\} \rightarrow \{0,1\}$ maps $x \in \{0,1\}$ to $y \in \{0,1\}$ with the probability: 
\begin{align*}
    \Pr[\mathcal{R}_{\epsilon}^W(x) = y] = 
    \begin{cases}
    \frac{e^\epsilon}{e^\epsilon + 1}   &   \text{(if $x=y$)} \\
    \frac{1}{e^\epsilon + 1}   &   \text{(otherwise)}.
    \end{cases}
\end{align*}
$\mathcal{R}_{\epsilon}^W$ 
provides $\epsilon$-LDP in Definition~\ref{def:LDP}, where $\mathcal{X} = \{0,1\}$. 
We refer to Warner's RR $\mathcal{R}_{\epsilon}^W$ with parameter $\epsilon$ as \textit{$\epsilon$-RR}. 


\smallskip
\noindent{\textbf{DP on Graphs.}}~~For graphs, we can consider two types of DP: 
\textit{edge DP} and \textit{node DP} \cite{Hay_ICDM09,Raskhodnikova_Encyclopedia16}. 
Edge DP hides the existence of one edge, whereas node DP hides the existence of one node along with its adjacent edges. 
In this paper, we focus on edge DP because existing one-round local triangle counting algorithms \cite{Imola_USENIX21,Imola_USENIX22,Ye_ICDE20,Ye_TKDE21} use edge DP. 
In other words, we are interested in 
how much the estimation error is reduced at the same value of $\epsilon$ in edge DP by shuffling. 
Although node DP is much stronger than edge DP, it is much harder to attain and often results in a much larger $\epsilon$ \cite{Chen_SIGMOD13,Sajadmanesh_arXiv22}. 
Thus, we leave an algorithm for shuffle node DP with small $\epsilon$ (e.g., $\epsilon \leq 1$) for future work. 
Another interesting avenue of future work is establishing a lower bound on the estimation error for node DP. 

Edge DP assumes that anyone (except for user $v_i$) can be an adversary who infers edges of user $v_i$ and that the adversary can obtain all edges except for edges of $v_i$ as background knowledge. 
Note that the central and local models have different definitions of neighboring data in edge DP. 
Specifically, edge DP in the central model \cite{Raskhodnikova_Encyclopedia16} considers two graphs that differ in one edge. 
In contrast, edge LDP 
\cite{Qin_CCS17} considers two neighbor lists that differ in one bit: 

\begin{definition} [$(\epsilon,\delta)$-edge DP \cite{Raskhodnikova_Encyclopedia16}] \label{def:edge_DP} 
Let $n \in \mathbb{N}$, $\epsilon \in \mathbb{R}_{\ge0}$, and $\delta \in [0,1]$. 
A randomized algorithm $\mathcal{M}$ with domain $\mathcal{G}$ provides \emph{$(\epsilon, \delta)$-edge DP} 
if for any two neighboring graphs $G, G' \in \mathcal{G}$ that differ in \textbf{one edge} and any $S \subseteq \mathrm{Range}(\mathcal{M})$, 
\begin{align*}
\Pr[\mathcal{M}(G) \in S] \leq e^\epsilon \Pr[\mathcal{M}(G') \in S] + \delta.
\end{align*}
\end{definition}

\begin{definition} [$\epsilon$-edge LDP~\cite{Qin_CCS17}] \label{def:edge_LDP} 
Let $\epsilon \in \mathbb{R}_{\ge0}$. 
A local randomizer $\mathcal{R}$ with domain $\{0,1\}$ provides \emph{$\epsilon$-edge LDP} if for any two neighbor lists $\mathbf{a}_i, \mathbf{a}'_i \in \{0,1\}^n$ that differ in \textbf{one bit} and any $S \subseteq \mathrm{Range}(\mathcal{R})$, 
\begin{align*}
\Pr[\mathcal{R}(\mathbf{a}_i) \in S] \leq e^\epsilon \Pr[\mathcal{R}(\mathbf{a}'_i) \in S].
\end{align*}
\end{definition}

As with edge LDP,  
we define \textit{element DP}, which considers two adjacency matrices that differ in one bit, in the central model:

\begin{definition} [$(\epsilon,\delta)$-element DP] \label{def:element_DP} 
Let $n \in \mathbb{N}$, $\epsilon \in \mathbb{R}_{\ge0}$, and $\delta \in [0,1]$. 
A randomized algorithm $\mathcal{M}$ with domain $\mathcal{G}$ provides \emph{$(\epsilon, \delta)$-element DP} 
if for any two neighboring graphs $G, G' \in \mathcal{G}$ that differ in \textbf{one bit} in the corresponding adjacency matrices $\mathbf{A}, \mathbf{A}' \in \{0,1\}^{n \times n}$
and any $S \subseteq \mathrm{Range}(\mathcal{M})$, 
\begin{align*}
\Pr[\mathcal{M}(G) \in S] \leq e^\epsilon \Pr[\mathcal{M}(G') \in S] + \delta.
\end{align*}
\end{definition}



Although element DP and edge DP have different definitions of neighboring data, they are closely related to each other:

\begin{proposition}\label{prop:element_edge_DP}
If a randomized algorithm $\mathcal{M}$ provides $(\epsilon, \delta)$-element DP, it also provides $(2\epsilon, 2\delta)$-edge DP. 
\end{proposition}
\begin{proof}
Adding or removing one edge affects two bits in an adjacency matrix. 
Thus, by group privacy \cite{DP}, any 
$(\epsilon, \delta)$-element DP 
algorithm $\mathcal{M}$ 
provides $(2\epsilon, 2\delta)$-edge DP. 
\end{proof}
Similarly, if 
a randomized algorithm $\mathcal{M}$ in the central model 
applies a local randomizer $\mathcal{R}$ providing $\epsilon$-edge LDP to each neighbor list $\mathbf{a}_i$ ($1 \leq i \leq n$), it provides $2\epsilon$-edge DP \cite{Imola_USENIX21}. 

In this work, we use the shuffling technique to provide $(\epsilon, \delta)$-element DP and then Proposition~\ref{prop:element_edge_DP} to provide $(2\epsilon, 2\delta)$-edge DP. 
We also compare our shuffle algorithms providing $(\epsilon, \delta)$-element DP and $(2\epsilon, 2\delta)$-edge DP with local algorithms providing $\epsilon$-edge LDP and $2\epsilon$-edge DP to see how much the estimation error is reduced by introducing the shuffle model and a very small $\delta$ ($\ll \frac{1}{n}$). 

\subsection{Shuffle Model}
\label{sub:shuffle}
We consider the following shuffle model. 
Each user $v_i \in V$ obfuscates her personal data 
using a local randomizer $\mathcal{R}$ providing $\epsilon_L$-LDP for $\epsilon_L \in \mathbb{R}_{\ge0}$. 
Note that $\mathcal{R}$ is common to all users. 
User $v_i$ encrypts the obfuscated data and sends it to a shuffler. 
Then, the shuffler randomly shuffles the encrypted data and sends the results to a data collector. 
Finally, the data collector decrypts them. 
The common assumption in the shuffle model is that the shuffler and the data collector do not collude with each other. 
Under this assumption, the shuffler cannot access the obfuscated data, and the data collector cannot link the obfuscated data to the users. 
Hereinafter, we omit the encryption/decryption process because it is clear from the context. 

We use the privacy amplification result by Feldman \textit{et al.} \cite{Feldman_FOCS21}:
\begin{theorem} [Privacy amplification by shuffling \cite{Feldman_FOCS21}] \label{thm:shuffle}
Let $n \in \mathbb{N}$ and $\epsilon_L \in \mathbb{R}_{\ge0}$. 
Let $\mathcal{X}$ be the set of input data for each user. 
Let $x_i \in \mathcal{X}$ be input data of the $i$-th user, and 
$x_{1:n} = (x_1, \cdots, x_n) \in \mathcal{X}^n$. 
Let $\mathcal{R}: \mathcal{X} \rightarrow \mathcal{Y}$ be a local randomizer providing $\epsilon_L$-LDP. 
Let $\mathcal{M}_S: \mathcal{X}^n \rightarrow \mathcal{Y}^n$ be an algorithm that given a dataset $x_{1:n}$, computes $y_i = \mathcal{R}(x_i)$ for $i \in [n]$, samples a uniform random permutation $\pi$ over $[n]$, and outputs $y_{\pi(1)}, \ldots, y_{\pi(n)}$. 
Then for any $\delta \in [0,1]$ such that $\epsilon_L \leq \log (\frac{n}{16 \log (2/\delta)})$, $\mathcal{M}_S$ provides $(\epsilon, \delta)$-DP, where
\begin{align}
\epsilon = f(n, \epsilon_L, \delta)
\label{eq:shuffle_epsilon_f}
\end{align}
and 
\begin{align}
f(n, \epsilon_L, \delta) = \log \left( 1 + \frac{e^{\epsilon_L}-1}{e^{\epsilon_L}+1} \left( \frac{8\sqrt{e^{\epsilon_L} \log(4/\delta)}}{\sqrt{n}} + \frac{8 e^{\epsilon_L}}{n} \right) \right).
\label{eq:shuffle_epsilon}
\end{align}
\end{theorem}
Thanks to the shuffling, the shuffled data $y_{\pi(1)}, \ldots, y_{\pi(n)}$ available to the data collector provides $(\epsilon, \delta)$-DP, where $\epsilon \ll \epsilon_L$. 

Feldman \textit{et al.} \cite{Feldman_FOCS21} also propose an efficient method to numerically compute a tighter upper bound than the closed-form upper bound in Theorem~\ref{thm:shuffle}. 
We use both the closed-form and numerical upper bounds in our experiments. 
Specifically, we use the numerical upper bounds in Section~\ref{sec:experiments} and compare the numerical bound with the closed-form bound in Appendix~\ref{sec:numerical_closed}. 

By (\ref{eq:shuffle_epsilon}), when we treat $\epsilon$ and $\delta$ as constants, $\epsilon_L$ can be expressed as 
$\epsilon_L = \log(n) + O(1)$. 
Similarly, the privacy amplification bound in \cite{Cheu_EUROCRYPT19} can also be expressed as $\epsilon_L = \log(n) + O(1)$. 
We use the bound in \cite{Feldman_FOCS21} because it is the state-of-the-art, as described in Section~\ref{sec:related}. 

Note that the privacy amplification result in \cite{Feldman_FOCS21} assumes that the data elements are shuffled \textit{before} applying the local randomizers, i.e., shuffle-then-randomize. 
However, 
shuffling outputs of the same local randomizers (randomize-then-shuffle) 
is equivalent to first shuffling input data and then applying the local randomizers (shuffle-then-randomize) as described in \cite{Erlingsson_SODA19}. 
Thus, if all users adopt the same local randomizer, the result in \cite{Feldman_FOCS21} can be applied to the randomize-then-shuffle model, which gives Theorem~\ref{thm:shuffle}. 

\subsection{Utility Metrics}
\label{sub:utility}
Following the existing work, we use the expectation of the $l_2$ loss \cite{Kairouz_ICML16,Murakami_USENIX19,Wang_USENIX17} in our theoretical analysis and the relative error \cite{Bindschaedler_SP16,Chen_CCS12,Xiao_SIGMOD11} in our experiments. 
The $l_2$ loss is a squared error between a true value and its estimate. 
Let $f: \mathcal{G} \rightarrow \mathbb{Z}_{\ge0}$ be a subgraph count function that can be instantiated by $f^\triangle$ or $f^\square$. 
Let $\hat{f}: \mathcal{G} \rightarrow \mathbb{R}$ be the corresponding estimator. 
Let 
$l_2^2: \mathbb{Z}_{\ge0} \times \mathbb{R} \rightarrow \mathbb{R}_{\ge0}$ 
be the expected $l_2$ loss function, which maps 
the true count $f(G)$ and the estimate $\hat{f}(G)$ to the expected $l_2$ loss. 
Then it can be expressed as $l_2^2(f(G), \hat{f}(G)) = \operatorname{\mathbb{E}}[(f(G) - \hat{f}(G))^2]$, where the expectation is taken over the randomness in the estimator $\hat{f}$. 
By the bias-variance decomposition \cite{mlpp}, the expected $l_2$ loss can be expressed as a summation of the squared bias $(\operatorname{\mathbb{E}}[\hat{f}(G)] - f(G))^2$ and the variance $\operatorname{\mathbb{V}}[\hat{f}(G)] = \operatorname{\mathbb{E}}[(\hat{f}(G) - \operatorname{\mathbb{E}}[\hat{f}(G)])]^2$. 
Thus, for an unbiased estimator $\hat{f}$ satisfying $\operatorname{\mathbb{E}}[\hat{f}(G)] = f(G)$, the expected $l_2$ loss is equal to the variance, i.e., $l_2^2(f(G), \hat{f}(G)) = \operatorname{\mathbb{V}}[\hat{f}(G)]$. 

Although the $l_2$ loss is suitable for theoretical analysis, it tends to be large when the number $n$ of users is large. 
This is because the true triangle and 4-cycle counts are very large when $n$ is large -- $f^\triangle(G) = O(n d_{max}^2)$ and $f^\square(G) = O(n d_{max}^3)$. 
Therefore, we use the relative error in our experiments. 
The relative error is an absolute error divided by the true value and is given by $\frac{|f^\triangle(G) - \hat{f}^\triangle(G)|}{\min\{f^\triangle(G), \eta\}}$, where $\eta \in \mathbb{R}_{\ge0}$ is a small positive value. 
Following the convention \cite{Bindschaedler_SP16,Chen_CCS12,Xiao_SIGMOD11}, we set $\eta = \frac{n}{1000}$. 
When the relative error is well below $1$, the estimate is accurate. 

\section{Shuffle Model for Graphs}
\label{sec:shuffle}
In this work, we apply the shuffle model to graph data to accurately estimate subgraph counts, such as triangles and 4-cycles. 
Section~\ref{sub:technical} explains our technical motivation. 
In particular, we explain why it is challenging to apply the shuffle model to graph data. 
Section~\ref{sub:wedge_shuffle} proposes a wedge shuffle technique to overcome the technical challenge. 

\subsection{Our Technical Motivation}
\label{sub:technical}

\begin{figure}[t]
  \centering
  \includegraphics[width=0.99\linewidth]{fig/shuffle.pdf}
  \vspace{-2mm}
  \caption{Shuffle model for graphs. 
 
 
  }
  \label{fig:shuffle_model}
\end{figure}

The shuffle model has been introduced to dramatically reduce the privacy budget $\epsilon$ (hence the estimation error at the same $\epsilon$) in tabular data \cite{Meehan_ICLR22,Wang_PVLDB20} or 
gradients 
\cite{Girgis_AISTATS21,Liu_AAAI21}. 
However, it is very challenging to apply the shuffle model to graph data, 
as explained below. 

Figure~\ref{fig:shuffle_model} shows the shuffle model for graph data, where each user $v_i$ has her neighbor list $\mathbf{a}_i \in \{0,1\}^n$. 
The main challenge here 
is that 
the shuffle model uses a \textit{standard} definition of LDP for the local randomizer and that a neighbor list is \textit{high-dimensional data}, i.e., $n$-dim binary string. 
Specifically, LDP in Definition~\ref{def:LDP} requires any pair of inputs $x$ and $x'$ to be indistinguishable; i.e., 
the inequality (\ref{eq:LDP}) must hold for all pairs of possible inputs. 
Thus, if we use the entire neighbor list 
as input data (i.e., $\mathbf{a}_i = x_i$ in Theorem~\ref{thm:shuffle}), 
either privacy or utility is destroyed for large $n$. 

To illustrate this, consider the following example. 
Assume that $n=10^5$ and $\delta=10^{-8}$. 
Each user $v_i$ applies 
$\epsilon_0$-RR 
with $\epsilon_0=1$ to each bit of her neighbor list $\mathbf{a}_i$. 
This mechanism is called the randomized neighbor list \cite{Qin_CCS17} and provides $\epsilon_0$-edge LDP. 
However, the privacy budget $\epsilon_L$ in the standard LDP (Definition~\ref{def:LDP}) 
is extremely large -- by group privacy \cite{DP}, $\epsilon_L = n \epsilon_0 = 10^5$. 
Because 
$\epsilon_L$ 
is much larger than $\log (\frac{n}{16 \log (2/\delta)}) = 8.09$, we cannot use the privacy amplification result in Theorem~\ref{thm:shuffle}. 
This is evident from the fact that the shuffled data $y_{\pi(1)}, \ldots, y_{\pi(n)}$ are easily re-identified when $n$ is large. 
If we use 
$\epsilon_0$-RR 
with $\epsilon_0 = \frac{1}{n}$, we can use the amplification result (as $\epsilon_L = n \epsilon_0 = 1$). 
However, it makes obfuscated data almost a random string and destroys the utility because $\epsilon_0$ is too small. 

In this work, we address this issue by introducing a basic technique, which we call \textit{wedge shuffling}. 

\subsection{Our Approach: Wedge Shuffling}
\label{sub:wedge_shuffle}


Figure ~\ref{fig:wedge_shuffle} shows the overview of our wedge shuffle technique. 
This technique calculates the number of wedges (2-hop paths) between 
a specific pair of users $v_i$ and $v_j$. 

Algorithm~\ref{alg:WShuffle} shows our wedge shuffle algorithm, which we call \textsf{WS}{}. 
Given users $v_i$ and $v_j$, 
each of the remaining users $v_k$ ($k \ne i, j$) 
calculates a \textit{wedge indicator} $w_{i-k-j} = a_{k,i} a_{k,j}$, 
which 
takes $1$ if 
a wedge $v_i$-$v_k$-$v_j$ exists and $0$ otherwise (line 2). 
Then, $v_k$ obfuscates $w_{i-k-j}$ using $\epsilon_L$-RR and sends it to the shuffler (line 3). 
The shuffler randomly shuffles the noisy wedges using a random permutation $\pi$ over 
$I_{-(i,j)}$ ($=[n]\setminus\{i,j\}$) 
to provide $(\epsilon, \delta)$-DP with $\epsilon \ll \epsilon_L$ (line 5). 
Finally, the shuffler sends the shuffled wedges to the data collector (line 6). 
The only information available to the data collector is  the number of wedges from $v_i$ to $v_j$, i.e., 
common friends of $v_i$ and $v_j$. 

\begin{figure}[t]
  \centering
  \includegraphics[width=0.99\linewidth]{fig/wedge_shuffle.pdf}
  \vspace{-2mm}
  \caption{Overview of wedge shuffling with inputs $v_i$ and $v_j$. 
 
  }
  \label{fig:wedge_shuffle}
\end{figure}

\setlength{\algomargin}{5mm}
\begin{algorithm}[t]
  \SetAlgoLined
  \KwData{Adjacency matrix $\mathbf{A} \in \{0,1\}^{n \times n}$, $\epsilon_L \in \mathbb{R}_{\ge0}$, user-pair $(v_i,v_j)$.}
  \KwResult{Shuffled wedges $\{y_{\pi(k)} | k \in I_{-(i,j)}\}$.}
  \ForEach{$k \in I_{-(i,j)}$}{
    [$v_k$] $w_{i-k-j} \leftarrow a_{k,i} a_{k,j}$\;
    [$v_k$] $y_k \leftarrow \mathcal{R}_{\epsilon_L}^W(x)(w_{i-k-j})$; Send $y_k$ to the shuffler\;
  }
  [s] Sample a random permutation $\pi$ over $I_{-(i,j)}$\;
  [s] Send $\{y_{\pi(k)} | k \in I_{-(i,j)}\}$ to the data collector\;
  [d] \KwRet{$\{y_{\pi(k)} | k \in I_{-(i,j)}\}$}
  \caption{Our wedge shuffle algorithm \textsf{WS}{}. 
  [$v_k$], [s], and [d] represent that the process is run by user $v_i$, the shuffler, and the data collector, respectively. 
  }\label{alg:WShuffle}
\end{algorithm}

Our wedge shuffling has two main features. 
First, 
the wedge indicator $w_{i-k-j}$ is \textit{one-dimensional binary data}. 
Therefore, it can be sent with small noise and small $\epsilon$, unlike the $n$-dimensional neighbor list. 
For example, when $n=10^5$, $\delta=10^{-8}$, and $\epsilon=1$, the value of $\epsilon_L$ in (\ref{eq:shuffle_epsilon_f}) and (\ref{eq:shuffle_epsilon}) is $\epsilon_L = 5.44$. 
In this case, $\epsilon_L$-RR rarely flips $w_{i-k-j}$ -- the flip probability is $0.0043$. 
In other words, the shuffled wedges are almost free of noise. 

Second, the wedge is the main component of many subgraphs such as triangles, 
$k$-triangles \cite{Karwa_PVLDB11}, 
3-hop paths \cite{Sun_CCS19}, 
and 4-cycles. 
For example, a triangle consists of one wedge and one edge, e.g., $v_i-v_{\pi(3)}-v_j$ and $(v_i, v_j)$ in Figure ~\ref{fig:wedge_shuffle}. 
More generally, a $k$-triangle consists of $k$ triangles sharing one edge. 
Thus, it can be decomposed into $k$ wedges and one edge. 
A 3-hop path consists of one wedge and one edge. 
A 4-cycle consists of two wedges, e.g., $v_i-v_{\pi(1)}-v_j$ and $v_i-v_{\pi(3)}-v_j$ in Figure ~\ref{fig:wedge_shuffle}. 
Because the shuffled wedges have little noise, we can accurately count these subgraphs based on wedge shuffling, compared to local algorithms in which all edges are noisy. 

In this work, we focus on triangles and 4-cycles and 
present algorithms with upper bounds on the estimation error based on our wedge shuffle technique. 




\section{Triangle Counting Based on Wedge Shuffling}
\label{sec:triangle}
Based on our wedge shuffle technique, we first propose a
one-round triangle counting algorithm.
Section~\ref{sub:triangle_overview} describes the overview of our algorithms.
Section~\ref{sub:wedge} proposes an algorithm for counting triangles involving a specific user-pair
as a building block of our triangle counting algorithm.
Section~\ref{sub:triangle} proposes
our triangle counting algorithm.
Section~\ref{sub:var_red} proposes a technique to significantly reduce the variance in our triangle counting algorithm.


\subsection{Overview}
\label{sub:triangle_overview}
Our wedge shuffle technique tells the data collector the number of common friends of $v_i$ and $v_j$.
However, this information is not sufficient to count triangles in the entire graph.
Therefore, we introduce
three additional techniques:
\textit{(i) sending local edges},
\textit{(ii) sampling independent user-pairs},
and
\textit{(iii) variance reduction by ignoring sparse user-pairs}.
Below, we briefly explain each technique.

\begin{figure}[t]
  \centering
  \includegraphics[width=0.99\linewidth]{fig/local_edges.pdf}
  \vspace{-2mm}
  \caption{Overview of our \textsf{WSLE}{} (Wedge Shuffling with Local Edges) algorithm with inputs $v_i$ and $v_j$.
 
 
  }
  \label{fig:local_edges}
\vspace{2mm}
  \centering
  \includegraphics[width=0.7\linewidth]{fig/sampling_pairs.pdf}
  \vspace{-2mm}
  \caption{Overview of our triangle counting algorithm.
  We use our
 
  \textsf{WSLE}{}
  algorithm with each user-pair.
 
  }
  \label{fig:triangle_count}
\end{figure}

\smallskip
\noindent{\textbf{Sending Local Edges.}}~~First, we consider the problem of counting triangles involving a specific user-pair $(v_i, v_j)$ and propose an algorithm to send
\textit{local edges} between $v_i$ and $v_j$, along with shuffled wedges, to the data collector.
We call this the \textsf{WSLE}{} (Wedge Shuffling with Local Edges) algorithm.

Figure~\ref{fig:local_edges} shows the overview of \textsf{WSLE}{}.
In this algorithm, users $v_i$ and $v_j$ obfuscate edge indicators $a_{i,j}$ and $a_{j,i}$, respectively, using $\epsilon$-RR and send them to the data collector directly (or through the shuffler without shuffling).
Then, the data collector calculates an unbiased estimate of the triangle count
from the shuffled wedges and the noisy edges.
Because $\epsilon$ is small, a large amount of noise is added to the edge indicators.
However, \textit{only one edge} is noisy (the other two have little noise) in any triangle the data collector sees.
This brings us an advantage over the one-round local algorithms in which all three edges are noisy.

\smallskip
\noindent{\textbf{Sampling Independent User-Pairs.}}~~Next, we consider the problem of counting triangles in the entire graph $G$.
A naive solution to this problem is to use our
\textsf{WSLE}{} algorithm
with all $\binom{n}{2}$ user-pairs as input.
However, it results in very large $\epsilon$ and $\delta$ because it uses each element of the adjacency matrix $\mathbf{A}$ many times.
To address this issue,
we propose a triangle counting algorithm that samples
independent user-pairs, which share no common users.

Figure~\ref{fig:triangle_count} shows the overview of our triangle algorithm.
The data collector sends the sampled user-pairs to users.
Then, users apply
\textsf{WSLE}{}
with each user-pair
and send the results to the data collector.
Finally, the data collector calculates an unbiased estimate of the triangle count from the results.
Because our triangle algorithm uses each element of
$\mathbf{A}$ \textit{at most once}, it provides $(\epsilon,\delta)$-element DP hence $(2\epsilon,2\delta)$-edge DP.
In addition, our triangle algorithm reduces the time complexity from $O(n^3)$ to $O(n^2)$ by sampling user-pairs rather than using all user-pairs.

We prove that the expected $l_2$ loss of our triangle counting algorithm is
$O(n^3)$ when we ignore the factor of $d_{max}$.
When we do not shuffle wedges, the
$l_2$ loss is
$O(n^4)$.
In addition, the $l_2$ loss of the existing one-round local algorithm \cite{Imola_USENIX22} with the same time complexity is $O(n^6)$, as shown in
Appendix~\ref{sec:upper}.
Thus, our
algorithm
provides a dramatic improvement over the local algorithms.

\smallskip
\noindent{\textbf{Variance Reduction.}}~~Although our
algorithm
dramatically
improves
the $l_2$ loss, the factor of $n^3$ may still be large.
Therefore, we propose a variance reduction technique that ignores sparse user-pairs, where either of the two users has a very small degree.
Our basic idea is that the number of triangles involving such a user-pair is very small
and can be approximated by $0$.
By ignoring the sparse user-pairs, we can significantly reduce the variance at the cost of introducing a small bias.
We prove that our variance reduction technique reduces the expected $l_2$ loss from $O(n^3)$ to $O(n^\gamma)$ where $\gamma\in[2,3)$ and makes one-round triangle counting more accurate.


\subsection{WSLE (Wedge Shuffling with Local Edges)}
\label{sub:wedge}

\noindent{\textbf{Algorithm.}}~~We first
propose
the \textsf{WSLE}{} algorithm
as a building block of our triangle counting algorithm.
\textsf{WSLE}{} counts
triangles involving a specific user-pair $(v_i,v_j)$.

Algorithm~\ref{alg:WSLE} shows \textsf{WSLE}{}.
Let $f_{i,j}^\triangle: \mathcal{G} \rightarrow \mathbb{Z}_{\ge0}$ be a function that takes $G \in \mathcal{G}$ as input and outputs the number $f_{i,j}^\triangle(G)$ of triangles involving $(v_i,v_j)$ in $G$.
Let $\hat{f}_{i,j}^\triangle(G) \in \mathbb{R}$ be an estimate of $f_{i,j}^\triangle(G)$.

\setlength{\algomargin}{5mm}
\begin{algorithm}[t]
  \SetAlgoLined
  \KwData{Adjacency matrix $\mathbf{A} \in \{0,1\}^{n \times n}$,
   
   
    $\epsilon \in \mathbb{R}_{\ge0}$, $\delta \in [0,1]$,
    user-pair $(v_i,v_j)$.
  }
  \KwResult{Estimate $\hat{f}_{i,j}^\triangle(G)$ of the number $f_{i,j}^\triangle(G)$ of triangles involving $(v_i,v_j)$.}
 
 
  $\epsilon_L \leftarrow \texttt{LocalPrivacyBudget}(n,\epsilon,\delta)$\;
  \tcc{Wedge shuffling}
  $\{y_{\pi(k)} | k \in I_{-(i,j)}\} \leftarrow$ \textsf{WS}{}$(\mathbf{A}, \epsilon_L, (v_{\sigma(i)}, v_{\sigma(i+1)}))$\;
  \tcc{Send local edges}
  [$v_i$] $z_i \leftarrow \mathcal{R}_{\epsilon}^W(x)(a_{i,j})$; Send $z_i$ to the data collector\;
  [$v_j$] $z_j \leftarrow \mathcal{R}_{\epsilon}^W(x)(a_{j,i})$; Send $z_j$ to the data collector\;
  \tcc{Calculate an unbiased estimate}
  [d] $q_L \leftarrow \frac{1}{e^{\epsilon_L}+1}$; $q \leftarrow
  \frac{1}{e^\epsilon+1}$\;
  [d] $\hat{f}_{i,j}^\triangle(G) \leftarrow \frac{(z_i+z_j-2q)\sum_{k \in
  I_{-(i,j)}} (y_{k} - q_L)}{2(1-2q)(1-2q_L)}$\;
  [d] \KwRet{$\hat{f}_{i,j}^\triangle(G)$}
  \caption{\textsf{WSLE}{}
 
  (Wedge Shuffling with Local Edges).
  \textsf{WS}{} is shown in Algorithm~\ref{alg:WShuffle}.
 
 
  }\label{alg:WSLE}
\end{algorithm}

We first call the function \texttt{LocalPrivacyBudget}, which calculates a local privacy budget $\epsilon_L$ from $n$, $\epsilon$, and $\delta$ (line 1).
Specifically, this function calculates $\epsilon_L$
such that $\epsilon$ is a closed-form upper bound (i.e., $\epsilon = f(n-2, \epsilon_L, \delta)$ in (\ref{eq:shuffle_epsilon_f})) or numerical upper bound in the shuffle model with $n-2$ users.
Given $\epsilon_L$, we can easily calculate the closed-form or numerical upper bound $\epsilon$ by (\ref{eq:shuffle_epsilon}) and the open source code in \cite{Feldman_FOCS21}\footnote{\url{https://github.com/apple/ml-shuffling-amplification}.}, respectively.
Thus, we can also easily calculate $\epsilon_L$ from $\epsilon$ by calculating a lookup table for pairs $(\epsilon, \epsilon_L)$ in advance.

Then, we run our wedge shuffle algorithm \textsf{WS}{} in Algorithm~\ref{alg:WShuffle} (line 2); i.e., each user $v_k \in I_{-(i,j)}$ sends her obfuscated wedge indicator
$y_k = \mathcal{R}_{\epsilon_L}^W(w_{i-k-j})$ to the shuffler, and the shuffler sends
shuffled wedge indicators $\{y_{\pi(k)} | k \in I_{-(i,j)}\}$ to the data collector.
Meanwhile, user $v_i$ obfuscates her edge indicator $a_{i,j}$ using $\epsilon$-RR $\mathcal{R}_{\epsilon}^W$ and sends the result $z_i = \mathcal{R}_{\epsilon}^W(a_{i,j})$ to the data collector
(line 3).
Similarly, $v_j$ sends $z_j = \mathcal{R}_{\epsilon}^W(a_{j,i})$ to the data collector (line 4).


Finally, the data collector estimates $f_{i,j}^\triangle(G)$ from $\{y_{\pi(k)} | k \in I_{-(i,j)}\}$, $z_i$, and $z_j$.
Specifically, the data collector calculates the estimate $\hat{f}_{i,j}^\triangle(G)$ as follows:
\begin{align}
    \textstyle{\hat{f}_{i,j}^\triangle(G) = \frac{(z_i+z_j-2q)\sum_{k \in I_{-(i,j)}} (y_{k} -
    q_L)}{2(1-2q)(1-2q_L)},}
    \label{eq:hfij_triangle}
\end{align}
where $q_L = \frac{1}{e^{\epsilon_L}+1}$ and $q = \frac{1}{e^\epsilon+1}$ (lines 5-6).
Note that this estimate involves simply summing over the set $\{y_{\pi(k)}\}$ and does not require knowing the value of $\pi$. This is
consistent with the shuffle model.
As we prove later, $\hat{f}_{i,j}^\triangle(G)$ in (\ref{eq:hfij_triangle}) is an unbiased estimate of $f_{i,j}^\triangle(G)$.

\smallskip
\noindent{\textbf{Theoretical Properties.}}~~Below, we show some theoretical properties of \textsf{WSLE}{}.
First, we prove that
the estimate $\hat{f}_{i,j}^\triangle(G)$
is unbiased:
\begin{theorem}
\label{thm:unbiased_I}
  For any indices $i,j \in [n]$, the estimate produced by
  $\textsf{WSLE}{}$ satisfies $\operatorname{\mathbb{E}}[\hat{f}_{i,j}^\triangle(G)] = f_{i,j}^\triangle(G)$.
\end{theorem}

Next,
we show the expected $l_2$ loss ($=$ variance). 
Recall that in the shuffle model, $\epsilon_L = \log n + O(1)$ when $\epsilon$ and $\delta$ are constants.
We show
the $l_2$ loss for a general case and for the shuffle model:
\begin{theorem}
\label{thm:l2-loss_I}
  For any indices $i,j \in [n]$, the estimate produced by
  \textsf{WSLE}{} provides the following utility guarantee:
\begin{align}
& l_2^2(f_{i,j}^\triangle(G), \hat{f}_{i,j}^\triangle(G)) = \operatorname{\mathbb{V}}[\hat{f}_{i,j}^\triangle(G)] \nonumber\\
  & \leq \frac{n q_L + q(1-2q_L)^2 d_{max}^2}{(1-2q)^2(1-2q_L)^2} \triangleq
  err_{\textsf{WSLE}}(n,d_{max},q,q_L).
\label{eq:l2_I_general}
\end{align}
  When $\epsilon$ and $\delta$ are constants and $\epsilon_L = \log n + O(1)$, we
  have
\begin{align}
  & err_{\textsf{WSLE}{}}(n,d_{max}, q, q_L) = O(d_{max}^2).
\label{eq:l2_I_shuffle}
\end{align}
\end{theorem}
Because \textsf{WSLE}{} is a building block for our triangle counting algorithms, we
introduce the notation $err_{\textsf{WSLE}{}}(n,d_{max}, q, q_L)$ for our upper bound
in~\eqref{eq:l2_I_general}. Observing
\eqref{eq:l2_I_general}, if we do not use the shuffling technique (i.e.,
$\epsilon_L = \epsilon$), then $err_{\textsf{WSLE}{}}(n,d_{max}, q, q_L) = O(n + d_{max}^2)$ when we
treat $\epsilon$ and $\delta$ as constants.
In contrast, in the shuffle model where we have $\epsilon_L = \log n + O(1)$,
then $err_{\textsf{WSLE}{}}(n,d_{max}, q, q_L) = O(d_{max}^2)$.
This means that wedge shuffling reduces the $l_2$ loss from $O(n + d_{max}^2)$ to $O(d_{max}^2)$, which is significant when $d_{max} \ll n$.

\subsection{Triangle Counting}
\label{sub:triangle}
\noindent{\textbf{Algorithm.}}~~Based on
\textsf{WSLE}{},
we propose an algorithm that counts triangles in the entire graph $G$.
We denote this algorithm by \textsf{WShuffle}$_{\triangle}${}, as it applies wedge shuffling to triangle counting.

Algorithm~\ref{alg:wshuffle_triangle} shows \textsf{WShuffle}$_{\triangle}${}.
First, the data collector samples independent user-pairs, which share no common users.
Specifically, it calls the function \texttt{RandomPermutation}, which samples a uniform random permutation $\sigma$ over $[n]$ (line 1).
Then, it samples independent user-pairs as
$(v_{\sigma(1)}, v_{\sigma(2)}), (v_{\sigma(3)}, v_{\sigma(4)}), \ldots, (v_{\sigma(2t-1)}, \allowbreak v_{\sigma(2t)})$, where $t \in [\lfloor \frac{n}{2} \rfloor]$.
The parameter $t$ represents the number of user-pairs and controls the trade-off between the $l_2$ loss and the time complexity;
when $t = \lfloor \frac{n}{2} \rfloor$, the $l_2$ loss is minimized and the time complexity is maximized.
The data collector sends the sampled user-pairs to users (line 2).

\setlength{\algomargin}{5mm}
\begin{algorithm}[t]
  \SetAlgoLined
  \KwData{Adjacency matrix $\mathbf{A} \in \{0,1\}^{n \times n}$, $\epsilon \in \mathbb{R}_{\ge0}$, $\delta \in [0,1]$, $t \in [\lfloor \frac{n}{2} \rfloor]$.
  }
  \KwResult{Estimate $\hat{f}^\triangle(G)$ of $f^\triangle(G)$.}
  \tcc{Sample independent user-pairs}
  [d] $\sigma \leftarrow$\texttt{RandomPermutation}$(n)$\;
  [d] Send $(v_{\sigma(1)}, v_{\sigma(2)}), \ldots, (v_{\sigma(2t-1)}, v_{\sigma(2t)})$ to users\;
  \ForEach{$i \in \{1, 3, \ldots, 2t-1\}$}{
    $\hat{f}_{\sigma(i), \sigma(i+1)}^\triangle(G) \leftarrow$ \textsf{WSLE}{}$(\mathbf{A}, \epsilon, \delta, (v_{\sigma(i)}, v_{\sigma(i+1)}))$\;
  }
  \tcc{Calculate an unbiased estimate}
  [d] $\hat{f}^\triangle(G) \leftarrow \frac{n(n-1)}{6t} \sum_{i=1, 3, \ldots, 2t-1} \hat{f}_{\sigma(i),\sigma(i+1)}^\triangle(G)$\;
  [d] \KwRet{$\hat{f}^\triangle(G)$}
  \caption{Our triangle counting algorithm \textsf{WShuffle}$_{\triangle}${}.
 
 
  \textsf{WSLE}{} is shown in Algorithm~\ref{alg:WSLE}.
  }\label{alg:wshuffle_triangle}
\end{algorithm}

Then, we run our wedge algorithm \textsf{WSLE}{} in Algorithm~\ref{alg:WSLE} with each sampled user-pair as input (lines 3-5).
Finally, the data collector estimates the triangle count $f^\triangle(G)$ as follows:
\begin{align}
    \textstyle{\hat{f}^\triangle(G) = \frac{n(n-1)}{6t} \sum_{i=1, 3, \ldots, 2t-1} \hat{f}_{\sigma(i),\sigma(i+1)}^\triangle(G)}
   \label{eq:hf_triangle_II}
\end{align}
(line 6). Later, we prove that $\hat{f}^\triangle(G)$ in (\ref{eq:hf_triangle_II}) is unbiased.

\smallskip
\noindent{\textbf{Theoretical Properties.}}~~We prove that
\textsf{WShuffle}$_{\triangle}${} provides DP:
\begin{theorem}
\label{thm:DP_II}
\textsf{WShuffle}$_{\triangle}${} provides $(\epsilon, \delta)$-element DP and $(2\epsilon, 2\delta)$-edge DP.
\end{theorem}
Theorem~\ref{thm:DP_II} comes from the fact that
\textsf{WSLE}{} with a user-pair $(v_i,v_j)$ provides $(\epsilon,\delta)$-DP for each element in the $i$-th and $j$-th columns of the adjacency matrix $\mathbf{A}$ and that \textsf{WShuffle}$_{\triangle}${} samples independent user-pairs, i.e., it uses each element of $\mathbf{A}$ at most once.

Note that running \textsf{WSLE}{} with all $\binom{n}{2}$ user-pairs requires a very large privacy budget --
it provides $((n-1) \epsilon, (n-1) \delta)$-DP, as it uses each element of $\mathbf{A}$ at most $n-1$ times.
The privacy budget is very large, even using the advanced composition \cite{DP,Kairouz_ICML15}.
We avoid this issue by sampling user-pairs that share no common users.

We also prove that
\textsf{WShuffle}$_{\triangle}${} provides an unbiased estimate:
\begin{theorem}
\label{thm:unbiased_II}
The estimate produced by \textsf{WShuffle}$_{\triangle}${} satisfies $\allowbreak \operatorname{\mathbb{E}}[\hat{f}^\triangle(G)] = f^\triangle(G)$.
\end{theorem}

Next, we analyze the expected $l_2$ loss ($=$ variance) of \textsf{WShuffle}$_{\triangle}${}.
This analysis is 
non-trivial 
because \textsf{WShuffle}$_{\triangle}${} samples each user-pair \textit{without replacement}.
In this case, the sampled user-pairs are not independent.
However, we can prove that
$t$ estimates in (\ref{eq:hf_triangle_II}) are negatively correlated with each other (Lemma~\ref{lem:sampling_replacement_var}
in Appendix~\ref{sub:l2-loss_II_proof}).
Thus, the variance of the sum of $t$ estimates
in (\ref{eq:hf_triangle_II}) is upper bounded by the sum of their variances, each of which is given by Theorem~\ref{thm:l2-loss_I}.
This brings us to the following result:

\begin{theorem}
\label{thm:l2-loss_II}
The estimate produced by \textsf{WShuffle}$_{\triangle}${} provides the following utility guarantee:
\begin{align}
& l_2^2(f^\triangle(G), \hat{f}^\triangle(G)) = \operatorname{\mathbb{V}}[\hat{f}^\triangle(G)] \nonumber\\
  &\leq
  \frac{n^4}{36t}err_{\textsf{WSLE}{}}(n,d_{max},q,q_L) +
  \frac{n^3}{36t}d_{max}^3, \label{eq:thm:l2_II}
\end{align}
where $err_{\textsf{WSLE}{}}(n,d_{max},q,q_L)$ is given by (\ref{eq:l2_I_general}).
  When $\epsilon$ and $\delta$ are constants,  $\epsilon_L = \log(n) + O(1)$, and
  $t = \lfloor\frac{n}{2}\rfloor$, we have
\begin{align}
l_2^2(f^\triangle(G), \hat{f}^\triangle(G))
  &\leq
  O(n^3d_{max}^2). \label{eq:thm:l2_II_simp}
\end{align}
\end{theorem}

The first and second terms in~(\ref{eq:thm:l2_II}) are caused by
Warner's RR
and the sampling of independent user-pairs, respectively.
In other words, the $l_2$ loss of \textsf{WShuffle}$_{\triangle}${} can be decomposed into two factors: the RR
and user-pair sampling.

For example, assume that $t = \lfloor \frac{n}{2} \rfloor$.
When we do not shuffle wedges (i.e., $\epsilon_L = \epsilon$), then
$err_{\textsf{WSLE}{}}(n,d_{max},q,q_L) = O(n + d_{max}^2)$, and
$l_2$ loss in (\ref{eq:thm:l2_II}) is $O(n^4 + n^3 d_{max}^2)$.
When we shuffle wedges, the $l_2$ loss is $O(n^3 d_{max}^2)$.
Thus, when we ignore the factor of $d_{max}$, our wedge shuffle technique reduces the $l_2$ loss from $O(n^4)$ to $O(n^3)$ in triangle counting.
The factor of $n^3$ is caused by the RR for local edges.
This is intuitive because a large amount of noise is added to the local edges.

Finally, we analyze the time complexity of \textsf{WShuffle}$_{\triangle}${}.
The time complexity of running \textsf{WSLE}{} with all $\binom{n}{2}$ user-pairs is $O(n^3)$, as there are $O(n^2)$ user-pairs in total and \textsf{WSLE}{} requires the time complexity of $O(n)$.
In contrast, the time complexity of \textsf{WShuffle}$_{\triangle}${} with $t = \lfloor \frac{n}{2} \rfloor$ is $O(n^2)$ because it samples $O(n)$ user-pairs.
Thus, \textsf{WShuffle}$_{\triangle}${} reduces the time complexity from $O(n^3)$ to $O(n^2)$ by user-pair sampling.
We can further reduce the time complexity at the cost of increasing the $l_2$ loss
by setting $t$ small,
i.e., $t \ll \lfloor \frac{n}{2} \rfloor$.

\subsection{Variance Reduction}
\label{sub:var_red}
\noindent{\textbf{Algorithm.}}~~\textsf{WShuffle}$_{\triangle}${}
achieves the $l_2$ loss of $O(n^3)$ when we ignore the factor of $d_{max}$.
To provide a smaller estimation error,
we propose a variance reduction technique that ignores sparse user-pairs.
We denote our triangle counting algorithm with the variance reduction technique by \textsf{WShuffle}$_{\triangle}^*${}.

As explained in Section~\ref{sub:triangle}, the factor of $n^3$ is caused by the RR for local edges.
However, most user-pairs $v_i$ and $v_j$ have a very small minimum degree
$\min\{d_i, d_j\} \ll d_{max}$,
and there is no edge $(v_i, v_j)$ between them in almost all cases.
In addition, even if there is an edge $(v_i, v_j)$, the number of triangles involving the sparse user-pair is very small
(at most $\min\{d_i, d_j\}$)
and can be approximated by $0$.
By ignoring such sparse user-pairs, we can dramatically reduce the variance of the RR for local edges at the cost of a small bias.
This is an intuition behind our variance reduction technique.

Algorithm~\ref{alg:wshuffle_triangle_vr} shows \textsf{WShuffle}$_{\triangle}^*${}.
This algorithm detects sparse user-pairs based on the degree information.
However, user $v_i$'s degree $d_i$ can leak the information about edges of $v_i$.
Thus, \textsf{WShuffle}$_{\triangle}^*${} calculates a differentially private estimate of $d_i$ within one round.

\setlength{\algomargin}{5mm}
\begin{algorithm}[t]
  \SetAlgoLined
  \KwData{Adjacency matrix $\mathbf{A} \in \{0,1\}^{n \times n}$, $\epsilon_1, \epsilon_2 \in \mathbb{R}_{\ge0}$, $\delta \in [0,1]$, $t \in [\lfloor \frac{n}{2} \rfloor]$, $c \in \mathbb{R}_{\ge0}$.
  }
  \KwResult{Estimate $\hat{f}^\triangle(G)$ of $f^\triangle(G)$.}
  \tcc{Sample independent user-pairs}
  [d] $\sigma \leftarrow$\texttt{RandomPermutation}$(n)$\;
  [d] Send $(v_{\sigma(1)}, v_{\sigma(2)}), \ldots, (v_{\sigma(2t-1)}, v_{\sigma(2t)})$ to users\;
  \ForEach{$i \in \{1, 3, \ldots, 2t-1\}$}{
    $\hat{f}_{\sigma(i), \sigma(i+1)}^\triangle(G) \leftarrow$ \textsf{WSLE}{}$(\mathbf{A}, \epsilon_2, \delta, (v_{\sigma(i)}, v_{\sigma(i+1)}))$\;
  }
  \tcc{Send noisy degrees}
  \For{$i=1$ \KwTo $n$}{
    [$v_i$] $\tilde{d}_i \leftarrow d_i + \textrm{Lap}(\frac{1}{\epsilon_1})$; Send $\tilde{d}_i$ to the data collector\;
  }
  \tcc{Calculate a variance-reduced estimate}
  [d] $\tilde{d}_{avg} \leftarrow \frac{1}{n} \sum_{i=1}^n \tilde{d}_i$; $d_{th} \leftarrow c \tilde{d}_{avg}$\;
  [d] $D \leftarrow \{i | i = 1, 3, \ldots, 2t-1,
  \min\{\tilde{d}_{\sigma(i)}, \tilde{d}_{\sigma(i+1)}\} > d_{th}\}$\;
  [d] $\hat{f}^\triangle(G) \leftarrow \frac{n(n-1)}{6t} \sum_{i \in D} \hat{f}_{\sigma(i),\sigma(i+1)}^\triangle(G)$\;
  [d] \KwRet{$\hat{f}^\triangle(G)$}
  \caption{Our triangle counting algorithm with variance reduction \textsf{WShuffle}$_{\triangle}^*${}.
 
  \textsf{WSLE}{} is shown in Algorithm~\ref{alg:WSLE}.
  }\label{alg:wshuffle_triangle_vr}
\end{algorithm}

Specifically, \textsf{WShuffle}$_{\triangle}^*${} uses two privacy budgets: $\epsilon_1, \epsilon_2 \in \mathbb{Z}_{\ge0}$.
The first budget $\epsilon_1$ is for privately estimating $d_i$, whereas the second budget $\epsilon_2$ is for
\textsf{WSLE}{}.
Lines 1 to 5 in Algorithm~\ref{alg:wshuffle_triangle_vr} are the same as those in Algorithm~\ref{alg:wshuffle_triangle}, except that Algorithm~\ref{alg:wshuffle_triangle_vr} uses $\epsilon_2$
to provide $(\epsilon_2, \delta)$-element DP.
After these processes, each user $v_i$ adds the Laplacian noise $\textrm{Lap}(\frac{1}{\epsilon_1})$ with mean $0$ and scale $\frac{1}{\epsilon_1}$ to her degree $d_i$ and sends the noisy degree $\tilde{d}_i$ ($= d_i + \textrm{Lap}(\frac{1}{\epsilon_1})$) to the data collector (lines 6-8).
Because the sensitivity \cite{DP} of $d_i$ (the maximum distance of $d_i$ between two neighbor lists that differ in one bit) is $1$, adding $\textrm{Lap}(\frac{1}{\epsilon_1})$ to $d_i$ provides $\epsilon_1$-element DP.

Then, the data collector estimates the average degree $d_{avg}$ as $\tilde{d}_{avg} = \frac{1}{n}\sum_{i=1}^n \tilde{d}_i$ and sets a threshold $d_{th}$ of the minimum degree to $d_{th} = c \tilde{d}_{avg}$, where $c \in \mathbb{R}_{\ge0}$ is a small positive number, e.g.,
$c \in [1,10]$ (line 9).
Finally, the data collector estimates $f^\triangle(G)$ as
\begin{align}
    \textstyle{\hat{f}^\triangle(G) = \frac{n(n-1)}{6t} \sum_{i \in D} \hat{f}_{\sigma(i),\sigma(i+1)}^\triangle(G),}
   \label{eq:hf_triangle_II_ast}
\end{align}
where
\begin{align*}
D = \{i | i = 1, 3, \ldots, 2t-1,
  \min\{\tilde{d}_{\sigma(i)}, \tilde{d}_{\sigma(i+1)}\} > d_{th}\}
\end{align*}
(lines 10-11).
The difference between (\ref{eq:hf_triangle_II}) and (\ref{eq:hf_triangle_II_ast}) is that (\ref{eq:hf_triangle_II_ast}) ignores sparse user-pairs $v_{\sigma(i)}$ and $v_{\sigma(i+1)}$ such that $\min\{\tilde{d}_{\sigma(i)}, \tilde{d}_{\sigma(i+1)}\} \leq d_{th}$.
Since
$d_{avg} \ll d_{max}$ in practice,
$d_{th} \ll d_{max}$
holds for small $c$.

The parameter $c$ controls the trade-off between the bias and variance of the estimate $\hat{f}^\triangle(G)$.
The larger $c$ is, the more user-pairs are ignored.
Thus, as $c$ increases, the bias is increased, and the variance is reduced.
In practice,
a small $c$ not less than $1$ 
results in a small $l_2$ error
because most real graphs are scale-free networks that have a power-law degree distribution \cite{NetworkScience}.
In the scale-free networks, most users' degrees are smaller than the average degree $d_{avg}$.
For example, in the BA (Barab\'{a}si-Albert) graph model \cite{NetworkScience,Hagberg_SciPy08},
most users' degrees are $\frac{d_{avg}}{2}$.
Thus, if we set
$c \in [1,10]$, for example,
then most user-pairs are ignored (i.e., $|D| \ll t$), which leads to a significant reduction of the variance at the cost of a small bias.

\smallskip
\noindent{\textbf{Theoretical Properties.}}~~As with \textsf{WShuffle}$_{\triangle}${},
\textsf{WShuffle}$_{\triangle}^*${} provides the following privacy guarantee:

\begin{theorem}
\label{thm:DP_II_ast}
\textsf{WShuffle}$_{\triangle}^*${} provides $(\epsilon_1 + \epsilon_2, \delta)$-element DP and $(2(\epsilon_1 + \epsilon_2), 2\delta)$-edge DP.
\end{theorem}

Next, we analyze the bias of \textsf{WShuffle}$_{\triangle}^*${}. Here, we assume most users have a small degree 
using parameters $\lambda \in \mathbb{R}_{\ge0}$ and $\alpha \in [0,1)$: 

\begin{theorem}
\label{thm:bias_II_ast}
  Suppose that in $G$, there exist $\lambda \in \mathbb{R}_{\ge0}$ and $\alpha \in [0,1)$
  such that at most $n^\alpha$ users have a degree larger than
  $\lambda d_{avg}$. Suppose \textsf{WShuffle}$_{\triangle}^*${} is run with $c \geq \lambda$.
  Then, the estimator produced by
  \textsf{WShuffle}$_{\triangle}^*${} provides the following bias guarantee:
\begin{align}
    Bias[\hat{f}^\triangle(G)]
    = |\operatorname{\mathbb{E}}[\hat{f}^\triangle(G)] - f^\triangle(G)|
    \leq
    \frac{n c^2 d_{avg}^2}{3} +
   
    \frac{4 n^\alpha}{3 \epsilon_1^2}.
   
    \label{eq:thm:bias_II_ast}
\end{align}
\end{theorem}
The values of $\lambda$ and $\alpha$ depend on the original graph $G$.
In the scale-free networks, $\alpha$ is small for a moderate value of $\lambda$.
For example, in the BA graph with $n=107614$ and $d_{avg}=200$ used in Appendix~\ref{sec:BA-graph}, $\alpha = 0.5$, $0.6$, $0.8$, and $0.9$ when $\lambda = 10.1$, $5.4$, $1.6$, and $0.9$, respectively.
When $c$ and $\epsilon_1$ are constants,
the bias can be expressed as $O(n d_{avg}^2)$.

Finally, we show the variance of \textsf{WShuffle}$_{\triangle}^*${}. 
This result assumes that $c$ is bigger $(= \frac{(1-\alpha) \log n}{\epsilon_1 d_{avg}})$ than $\lambda$. 
We assume this 
because otherwise, many sparse users (with $d_i \leq \lambda d_{avg}$) have a noisy degree $\tilde{d}_i \geq c \tilde{d}_{avg}$, causing the set $D$ to
be noisy. In practice, the gap between $c$ and $\lambda$ is small because 
$\log n$ is much smaller than $d_{avg}$. 


\begin{theorem}
\label{thm:var_II_ast}
  Suppose that in $G$, there exist $\lambda \in \mathbb{R}_{\ge0}$ and $\alpha \in [0,1)$
  such that at most $n^\alpha$ users have a degree larger than
  $\lambda d_{avg}$. Suppose \textsf{WShuffle}$_{\triangle}^*${} is
  run with $c \geq \lambda + \frac{(1-\alpha) \log n}{\epsilon_1 d_{avg}}$.
  Then, the estimator produced by
  \textsf{WShuffle}$_{\triangle}^*${} provides the following variance guarantee:
\begin{align}
  &\operatorname{\mathbb{V}}[\hat{f}_*^\triangle(G)] \leq \nonumber\\
  &\frac{n^2 d_{max}^4}{9} + \frac{2n^{2+2\alpha}}{9t}err_{\textsf{WSLE}}(n,d_{max},q,q_L) + \frac{n^{2+\alpha}
  d_{max}^3 }{36t}.
  \label{eq:thm:l2_II_ast}
\end{align}
When $\epsilon_1$, $\epsilon_2$, and $\delta$ are constants,
$\epsilon_L = \log n + O(1)$, and $t = \lfloor\frac{n}{2}\rfloor$,
\begin{align}
  \operatorname{\mathbb{V}}[\hat{f}_*^\triangle(G)] \leq O(n^2 d_{max}^4 + n^{1+2\alpha} d_{max}^2).
  \label{eq:thm:l2_II_ast_simp}
\end{align}
\end{theorem}


The first\footnote{The first term in (\ref{eq:thm:l2_II_ast}) is actually $\frac{(\sum_{i=1}^n d_i^2)^2}{9}$ and is much smaller than $\frac{n^2 d_{max}^4}{9}$. 
We express it as $O(n^2 d_{max}^4)$ in (\ref{eq:thm:l2_II_ast_simp}) for simplicity. 
See Appendix~\ref{sub:var_II_ast_proof} for details.}, second, and third terms in (\ref{eq:thm:l2_II_ast}) are caused by the randomness in the choice of $D$, the RR, and user-pair sampling, respectively. 
By (\ref{eq:thm:l2_II_ast_simp}), our variance reduction technique reduces the variance from $O(n^3)$ to $O(n^\gamma)$ where $\gamma\in[2,3)$ when we ignore the factor of $d_{max}$.
Because the expected $l_2$ loss is the sum of the squared bias and the variance, it is also $O(n^\gamma)$.

The value of $\gamma$ in our bound $O(n^\gamma)$ depends on the parameter $c$ in \textsf{WShuffle}$_{\triangle}^*${}.
For example, in the BA graph ($n=107614$, $d_{avg} = 200$), $\gamma=2$, $2.2$, $2.6$, and $2.8$ ($\alpha=0.5$, $0.6$, $0.8$, and $0.9$) when $c=10.4$, $5.6$, $1.7$, and $1.0$, respectively, and $\epsilon_1=0.1$.
Thus, the variance decreases with increase in $c$.
However, by (\ref{eq:thm:bias_II_ast}), a larger $c$ results in a larger bias. 
In our experiments, we show that \textsf{WShuffle}$_{\triangle}^*${} provides a small estimation error when $c=1$ to $4$.
When $c=1$, \textsf{WShuffle}$_{\triangle}^*${} empirically works well despite a large $\gamma$ because most users' degrees are smaller than $d_{avg}$ in practice, as explained above.
This indicates that our upper bound in (\ref{eq:thm:l2_II_ast_simp}) might not be tight when $c$ is around $1$.
Improving the bound is left for future work.



\subsection{Summary}
\label{sub:summary}

Table~\ref{tab:upper_bounds_triangle} summarizes the
performance guarantees
of one-round triangle algorithms providing edge DP.
Here,
we consider a variant of \textsf{WShuffle}$_{\triangle}${} that does not shuffle wedges (i.e., $\epsilon_L = \epsilon$) as a one-round local algorithm.
We call this variant \textsf{WLocal}$_{\triangle}${} (Wedge Local).
We also show the variance of \textsf{ARR}$_{\triangle}${} \cite{Imola_USENIX22} and \textsf{RR}$_{\triangle}${} \cite{Imola_USENIX21}.
The time complexity of \textsf{RR}$_{\triangle}${} is $O(n^3)$\footnote{Technically speaking, the algorithms of \textsf{RR}$_{\triangle}${} and
the one-round local algorithms in \cite{Ye_ICDE20,Ye_TKDE21} involve counting the number of triangles in a dense
graph. This can be done in time $O(n^\omega)$, where $\omega \in [2,3)$ and $O(n^\omega)$ is the time required for matrix multiplication. However, these algorithms are of theoretical interest,
and they do not outperform naive matrix multiplication except
for very large matrices~\cite{Alman_2021}. Thus, we assume implementations that use naive matrix multiplication in
$O(n^3)$ time.}, and that
of \textsf{ARR}$_{\triangle}${} is $O(n^2)$ when we set the sampling probability $p_0 \in [0,1]$ of the ARR to $p_0=O(n^{-1/3})$.
We prove the variance of \textsf{ARR}$_{\triangle}${} in this case and \textsf{RR}$_{\triangle}${}
in Appendix~\ref{sec:upper}.
We do not show the other one-round local algorithms \cite{Ye_ICDE20,Ye_TKDE21} in Table~\ref{tab:upper_bounds_triangle} for two reasons: (i) they have the time complexity of $O(n^3)$ and suffer from a larger estimation error than \textsf{RR}$_{\triangle}${} \cite{Imola_USENIX22};
(ii) their upper-bounds on the variance and bias are unclear.

\begin{table}[t]
  \caption{Performance guarantees
 
  of one-round triangle counting algorithms providing edge DP.
 
  $\alpha \in [0,1)$. 
  See also footnote 2 for the variance of \textsf{WShuffle}$_{\triangle}^*${}.
 
  }
  \vspace{-4mm}
  \centering
  \begin{tabular}{|l|c|c|c|c|}
    \hline
    Algorithm & \hspace{-0.8mm}Model\hspace{-0.8mm} & Variance & Bias & Time \\ \hline
   
   
   
    \textsf{WShuffle}$_{\triangle}^*${} & \hspace{-0.8mm}shuffle\hspace{-0.8mm} & \hspace{-1mm}$O(n^2
    d_{max}^4 \hspace{-0.5mm}+\hspace{-0.5mm} n^{1+2\alpha} d_{max}^2)$\hspace{-1mm} & \hspace{-1.5mm}$O(n d_{avg}^2)$\hspace{-1.5mm} & $O(n^2)$\\
    \hline
   
    \textsf{WShuffle}$_{\triangle}${} & \hspace{-0.8mm}shuffle\hspace{-0.8mm} & $O(n^3 d_{max}^2)$ & $0$ & $O(n^2)$ \\ \hline
    \textsf{WLocal}$_{\triangle}${} & local & $O(n^4 + n^3 d_{max}^2)$ & $0$ & $O(n^2)$ \\ \hline
    \textsf{ARR}$_{\triangle}${} \cite{Imola_USENIX22} & local & $O(n^6)$ & $0$ & $O(n^2)$ \\ \hline
    \textsf{RR}$_{\triangle}${} \cite{Imola_USENIX21} & local & $O(n^4)$ & $0$ & $O(n^3)$ \\ \hline
  \end{tabular}
  \label{tab:upper_bounds_triangle}
\end{table}

Table~\ref{tab:upper_bounds_triangle} shows that our \textsf{WShuffle}$_{\triangle}^*${} dramatically outperforms the three local algorithms -- when we ignore $d_{max}$, the $l_2$ loss of \textsf{WShuffle}$_{\triangle}^*${} is 
$O(n^\gamma)$ where $\gamma\in[2,3)$, 
whereas that of the local algorithms is $O(n^4)$ or $O(n^6)$.
We also show this through experiments.





\section{4-Cycle Counting Based on Wedge Shuffling}
\label{sec:4cycle}

Next, we propose a one-round 4-cycle counting algorithm in the shuffle model.
Section~\ref{sub:4cycle_overview} explains its overview.
Section~\ref{sub:4cycle_details} proposes our 4-cycle counting algorithm and shows its theoretical properties.

\subsection{Overview}
\label{sub:4cycle_overview}
We apply our wedge shuffling technique to 4-cycle counting with two additional techniques: \textit{(i) bias correction} and \textit{(ii) sampling independent user-pairs}. 
Below, we briefly explain each of them. 

\smallskip
\noindent{\textbf{Bias Correction.}}~~As with triangles, we begin with the problem of counting 4-cycles involving specific users $v_i$ and $v_j$. 
We can leverage the noisy wedges output by our wedge shuffle algorithm \textsf{WS}{} to estimate such a 4-cycle count. 
Specifically, let $f^\square_{i,j}: \mathcal{G} \rightarrow \mathbb{Z}_{\ge0}$ be a function that, given $G \in \mathcal{G}$, outputs the number $f^\square_{i,j}(G)$ of $4$-cycles for which users $v_i$ and $v_j$ are opposite nodes, i.e. the number of unordered pairs $(k,k')$ such that $v_i-v_k-v_j-v_{k'}-v_i$ is a path in $G$.
Each pair $(k,k')$ satisfies the above requirement if and only if $v_i-v_k-v_j$ and $v_i-v_{k'}-v_j$ are wedges in $G$.
Thus, we have $f^\square_{i,j}(G) = \binom{f^\wedge_{i,j}}{2}$, where $f^\wedge_{i,j}$ is the number of wedges between $v_i$ and $v_j$. 
Based on this, we calculate an unbiased estimate $\hat{f}^\wedge_{i,j}$ of the wedge count using \textsf{WS}{}. 
Then, we calculate an estimate of the 4-cycle count as $\binom{\hat{f}^\wedge_{i,j}}{2}$. 
Here, it should be noted that the estimate $\binom{\hat{f}^\wedge_{i,j}}{2}$ 
is \textit{biased}, as proved later. 
Therefore, 
we perform bias correction -- we subtract a positive value from the estimate to obtain an unbiased estimate $\hat{f}_{i,j}^\square(G)$ of the 4-cycle count. 

Note that unlike \textsf{WSLE}{}, no edge between $(v_i,v_j)$ needs to be sent. 
In addition, thanks to the privacy amplification by shuffling, all wedges can be sent with 
small noise.

\smallskip
\noindent{\textbf{Sampling Independent User-Pairs.}}~~Having an estimate $\hat{f}_{i,j}^\square(G)$, we turn our attention to estimating 4-cycle count $f^\square(G)$ in the entire graph $G$.
As with triangles, a naive solution using estimates $\hat{f}_{i,j}^\square(G)$ for all $\binom{n}{2}$ user-pairs $(v_i,v_j)$ results in very large $\epsilon$ and $\delta$. 
To avoid this, we sample independent user-pairs and obtain an unbiased estimate of $f^\square(G)$ from them. 

\subsection{4-Cycle Counting}
\label{sub:4cycle_details}

\setlength{\algomargin}{5mm}
\begin{algorithm}[t]
  \SetAlgoLined
  \KwData{Adjacency matrix $\mathbf{A} \in \{0,1\}^{n \times n}$, $\epsilon \in \mathbb{R}_{\ge0}$, $\delta \in [0,1]$, $t \in [\lfloor \frac{n}{2} \rfloor]$.
  }
  \KwResult{Estimate $\hat{f}^\square(G)$ of $f^\square(G)$.}
  $\epsilon_L \leftarrow \texttt{LocalPrivacyBudget}(n,\epsilon,\delta)$\;
  [d] $q_L \leftarrow \frac{1}{e^{\epsilon_L}+1}$\;
  \tcc{Sample independent user-pairs}
  [d] $\sigma \leftarrow$\texttt{RandomPermutation}$(n)$\;
  [d] Send $(v_{\sigma(1)}, v_{\sigma(2)}), \ldots, (v_{\sigma(2t-1)}, v_{\sigma(2t)})$ to users\;
  \ForEach{$i \in \{1, 3, \ldots, 2t-1\}$}{
    $\{y_{\pi_i(k)} | k \hspace{-1mm}\in I_{-(\sigma(i),\sigma(i+1))}\} \hspace{-1mm}\leftarrow\hspace{-1mm}$ \textsf{WS}{}$(\mathbf{A}, \epsilon_L, (v_{\sigma(i)}, v_{\sigma(i+1)}))$\;
    [d] $\hat{f}_{\sigma(i),\sigma(i+1)}^\wedge(G) \leftarrow \sum_{k \in I_{-(\sigma(i),\sigma(i+1))}} \frac{y_{k} - q_L}{1-2q_L}$\;
    [d] $\hat{f}_{\sigma(i),\sigma(i+1)}^\square(G) \leftarrow  \frac{\hat{f}_{\sigma(i),\sigma(i+1)}^\wedge(G)(\hat{f}_{\sigma(i),\sigma(i+1)}^\wedge(G)-1)}{2} - \frac{n-2}{2}\frac{q_L(1-q_L)}{(1-2q_L)^2}$\;
  }
  \tcc{Calculate an unbiased estimate}
  [d] $\hat{f}^\square(G) \leftarrow \frac{n(n-1)}{4t}\sum_{i=1,3,\ldots,2t-1} \hat{f}_{\sigma(i),\sigma(i+1)}^\square(G)$\;
 
  [d] \KwRet{$\hat{f}^\square(G)$}
  \caption{Our 4-cycle counting algorithm \textsf{WShuffle}$_{\square}${}.
 
  \textsf{WS}{} is shown in Algorithm~\ref{alg:WShuffle}.
  }\label{alg:wshuffle_cycle}
\end{algorithm}

\noindent{\textbf{Algorithm.}}~~Algorithm \ref{alg:wshuffle_cycle} shows our 4-cycle counting algorithm. 
We denote it by \textsf{WShuffle}$_{\square}${}. 
First, we set a local privacy budget $\epsilon_L$ from $n$, $\epsilon$, and $\delta$ in the same way as \textsf{WSLE}{} (line 1). 
Then, we sample $t$ disjoint pairs of users using the permutation $\sigma$ (lines 3-4). 
Each pair is given by $(v_{\sigma(i), \sigma(i+1)})$ for $i \in \{1, 3, \ldots, 2t-1\}$.

For each $i \in \{1, 3, \ldots, 2t-1\}$, we compute an unbiased estimate $\hat{f}_{\sigma(i),
\sigma(i+1)}^\square(G)$ of the 4-cycle count involving $v_{\sigma(i)}$ and $v_{\sigma(i+1)}$ (lines 5-9). 
To do this, 
we call \textsf{WS}{} on
$(v_{\sigma(i)} v_{\sigma(i+1)})$ to obtain an unbiased estimate $\hat{f}_{\sigma(i),
\sigma(i+1)}^\wedge(G)$ of the wedge count (lines 6-7). 
We calculate an estimate $\hat{f}_{i,j}^\wedge(G)$ of the number $f_{i,j}^\wedge(G)$ of wedges between $v_i$ and $v_j$ in $G$ as follows:
\begin{equation}\label{eq:wedge-estimate}
  \textstyle{\hat{f}^\wedge_{i,j}(G) = \sum_{k \in I_{-(i,j)}} \frac{y_k - q_L}{1-2q_L}.}
\end{equation}
Later, we will prove that $\hat{f}^\wedge_{i,j}(G)$ is an unbiased estimator. 
As with (\ref{eq:hfij_triangle}), this estimate involves the sum over the set $\{y_{\pi(k)}\}$ and does not require knowing the permutation $\pi$ produced by the shuffler. 
Then, we obtain an unbiased estimator of $f_{\sigma(i), \sigma(i+1)}^\square$ as follows:
\begin{equation}\label{eq:4-cycle-ij-estimate}
  \textstyle{\hat{f}_{i,j}^\square(G) = \frac{\hat{f}^\wedge_{i, j}(G)(\hat{f}^\wedge_{i, j}(G)-1)}{2} - \frac{n-2}{2}\frac{q_L(1-q_L)}{(1-2q_L)^2}}
\end{equation}
(line 8). 
Note that 
there is a quadratic relationship between $f_{i, j}^\square(G)$ and $f_{i, j}^\wedge(G)$, i.e., $f_{i, j}^\square(G) = \binom{f_{i, j}^\wedge(G)}{2}$. 
Thus, 
even though $\hat{f}_{i, j}^\wedge(G)$ is unbiased, 
we must subtract a term from $\binom{\hat{f}_{i, j}^\wedge(G)}{2}$ (i.e., bias correction) to obtain an unbiased estimator $\hat{f}_{i, j}^\square(G)$. 
This forms the righthand side of
~\eqref{eq:4-cycle-ij-estimate} and ensures that $\hat{f}_{i,j}^\square(G)$ is unbiased. 

Finally, 
we sum and scale 
$\hat{f}_{\sigma(i), \sigma(i+1)}^\square(G)$ for each $i$ to obtain an estimate $\hat{f}^\square(G)$ of the 4-cycle count $f^\square(G)$ in the entire graph $G$: 

\begin{equation}\label{eq:4-cycle-estimate}
  \textstyle{\hat{f}^\square(G) = \frac{n(n-1)}{4t}\sum_{i=1,3,\ldots,2t-1} \hat{f}_{\sigma(i),\sigma(i+1)}^\square(G)}
\end{equation}
(line 10). As we prove later, $\hat{f}^\square(G)$ is an unbiased estimate of $f^\square(G)$. 

\smallskip
\noindent{\textbf{Theoretical Properties.}}~~First, \textsf{WShuffle}$_{\square}${} guarantees DP:
\begin{theorem}\label{thm:DP_IV}
\textsf{WShuffle}$_{\square}${} provides $(\epsilon, \delta)$-element DP and $(2\epsilon, 2\delta)$-edge DP.
\end{theorem}

In addition, thanks to the design of~\eqref{eq:4-cycle-ij-estimate}, we can show that \textsf{WShuffle}$_{\square}${} produces an unbiased estimate of $f^\square(G)$:
\begin{theorem}\label{thm:unbiased_IV}
  The estimate produced by \textsf{WShuffle}$_{\square}${} satisfies $\allowbreak \operatorname{\mathbb{E}}[\hat{f}^\square(G)] =
  f^\square(G)$.
\end{theorem}

Finally, we show the expected $l_2$ loss ($=$ variance) of $f^\square(G)$:

\begin{theorem}\label{thm:l2-loss_IV}
  The estimate produced by \textsf{WShuffle}$_{\square}${} satisfies 
  \begin{align}
    & l_2^2(f^\square(G), \hat{f}^\square(G)) = \operatorname{\mathbb{V}}[\hat{f}^\square(G)] \nonumber \\
    &\leq \frac{9n^5 q_L(d_{max} +
    nq_L)^2}{16t(1-2q_L)^4} + \frac{n^3d_{max}^6}{64t}. \label{eq:4-cycle-var}
  \end{align}
  When $\epsilon$ and $\delta$ are constants, $\epsilon_L = \log n + O(1)$, and $t
  = \lfloor \frac{n}{2} \rfloor$, we have
  \begin{equation}\label{eq:4-cycle-var-simp}
    l_2^2(f^\square(G), \hat{f}^\square(G)) = \operatorname{\mathbb{V}}[\hat{f}^\square(G)] 
    = O\left( n^3 d_{max}^2 + n^2d_{max}^6 \right).
  \end{equation}
\end{theorem}
The first and second terms in (\ref{eq:4-cycle-var}) are caused by the RR and the sampling of independent user-pairs, respectively. 

\begin{table}[t]
  \caption{Performance guarantees of one-round 4-cycle counting algorithms providing edge DP.
 
  }
  \vspace{-4mm}
  \centering
  \begin{tabular}{|l|c|c|c|c|}
    \hline
    Algorithm & Model & Variance & Bias & Time \\ \hline
    \textsf{WShuffle}$_{\square}${} & shuffle & $O(n^3 d_{max}^2 + n^2 d_{max}^6)$ & $0$ & $O(n^2)$ \\ \hline
    \textsf{WLocal}$_{\square}${} & local & $O(n^6 + n^2 d_{max}^6)$ & $0$ & $O(n^2)$ \\ \hline
  \end{tabular}
  \label{tab:upper_bounds_4cycle}
\end{table}

\smallskip
\noindent{\textbf{Summary.}}~~Table~\ref{tab:upper_bounds_4cycle} summarizes the performance guarantees of the 4-cycle counting algorithms. 
As a one-round local algorithm, we consider a local model version of \textsf{WShuffle}$_{\square}${} that does not shuffle wedges (i.e., $\epsilon_L = \epsilon$). 
We denote it by \textsf{WLocal}$_{\square}${}. 
To our knowledge, \textsf{WLocal}$_{\square}${} is the first local 4-cycle counting algorithm. 

By (\ref{eq:4-cycle-var}), when $t = \lfloor \frac{n}{2} \rfloor$, the $l_2$ loss of \textsf{WLocal}$_{\square}${} can be expressed as $O(n^6 + n^2 d_{max}^6)$. 
Thus, our wedge shuffle technique dramatically reduces the $l_2$ loss from $O(n^6 + n^2 d_{max}^6)$ to $O(n^3 d_{max}^2 + n^2 d_{max}^6)$. 
Note that the square of the true count $f^\square(G)$ is $O(n^2 d_{max}^6)$. 
This indicates that our \textsf{WShuffle}$_{\square}${} may not work well in an extremely sparse graph where $d_{max} < n^{\frac{1}{4}}$. 
However, $d_{max} \gg n^{\frac{1}{4}}$ holds in most social graphs; e.g., the maximum number $d_{max}$ of friends is much larger than $100$ when $n=10^8$. 
In this case, \textsf{WShuffle}$_{\square}${} can accurately estimate the 4-cycle count, as shown in our experiments. 



\section{Experimental Evaluation}
\label{sec:experiments}
Based on the performance guarantees summarized in Tables~\ref{tab:upper_bounds_triangle} and \ref{tab:upper_bounds_4cycle}, we pose the following research questions: 

\begin{description}[leftmargin=8.75mm]
    \item[RQ1.] 
   
    How much do our entire algorithms (\textsf{WShuffle}$_{\triangle}^*${} and \textsf{WShuffle}$_{\square}${}) outperform the local algorithms?
   
    \item[RQ2.] 
   
    For triangles, how much does our variance reduction technique decrease the relative error?
    \item[RQ3.] 
   
    How small relative errors do our entire algorithms 
   
    achieve with a small privacy budget?
\end{description}
We designed experiments to answer these questions. 

\subsection{Experimental Set-up}
\label{sub:set-up}
We used the following two real graph datasets: 
\begin{itemize}
    \item \textbf{Gplus}: The first dataset is the Google+ dataset \cite{McAuley_NIPS12} denoted by \textsf{Gplus}{}. 
   
   
    This dataset includes a social graph $G=(V,E)$ with $n=107614$ users and $12238285$ edges, where an edge $(v_i, v_j) \in E$ represents that a user $v_i$ follows or is followed by $v_j$. 
    The average and maximum degrees are $d_{avg} = 227.4$ and $d_{max} = 20127$, respectively. 
    \item \textbf{IMDB}: The second dataset is the IMDB (Internet Movie Database)~\cite{IMDB_GD05} denoted by \textsf{IMDB}{}. 
    This dataset includes a bipartite graph between $896308$ actors and $428440$ movies. 
    From this, we extracted a graph $G=(V,E)$ with $n=896308$ actors and $57064358$ edges, where an edge represents that two actors have played in the same movie. 
    The average and maximum degrees are $d_{avg} = 127.3$ and $d_{max} = 15451$, respectively; i.e., \textsf{IMDB}{} is more sparse than \textsf{Gplus}{}. 
\end{itemize}
In Appendix~\ref{sec:BA-graph}, we also evaluate our algorithms using the Barab\'{a}si-Albert graphs \cite{NetworkScience,Hagberg_SciPy08}, which have a power-law degree distribution. 
Moreover, in Appendix~\ref{sec:bipartite}, we evaluate our 4-cycle algorithms using bipartite graphs generated from \textsf{Gplus}{} and \textsf{IMDB}{}. 

For triangle counting, we evaluated the following four one-round algorithms: \textsf{WShuffle}$_{\triangle}^*${}, \textsf{WShuffle}$_{\triangle}${}, \textsf{WLocal}$_{\triangle}${}, and \textsf{ARR}$_{\triangle}${} \cite{Imola_USENIX22}. 
We did not evaluate \textsf{RR}$_{\triangle}${} \cite{Imola_USENIX21}, because it was too inefficient -- it was reported in \cite{Imola_USENIX21} that when $n=10^6$, \textsf{RR}$_{\triangle}${} would require over $30$ years even on a supercomputer. 
The same applies to the one-round local algorithms in \cite{Ye_ICDE20,Ye_TKDE21} with the same time complexity ($=O(n^3)$). 

For 4-cycle counting, we compared \textsf{WShuffle}$_{\square}${} with \textsf{WLocal}$_{\square}${}. 
Because \textsf{WLocal}$_{\square}${} is the first local 4-cycle counting algorithm (to our knowledge), we did not evaluate other algorithms. 

In our shuffle algorithms \textsf{WShuffle}$_{\triangle}^*${}, \textsf{WShuffle}$_{\triangle}${}, and \textsf{WShuffle}$_{\square}${}, we set $\delta = 10^{-8}$ ($\ll \frac{1}{n}$) and $t=\frac{n}{2}$. 
We used the numerical upper bound in \cite{Feldman_FOCS21} for calculating $\epsilon$ in the shuffle model. 
In Appendix~\ref{sec:numerical_closed}, we compare the numerical bound with the closed-form bound in Theorem~\ref{thm:shuffle}. 
In \textsf{WShuffle}$_{\triangle}^*${}, we set 
$c\in[0.1,4]$ 
and 
divided the total privacy budget $\epsilon$ as $\epsilon_1 = \frac{\epsilon}{10}$ and $\epsilon_2 = \frac{9\epsilon}{10}$. 
Here, we assigned a small budget to $\epsilon_1$ because a degree $d_i$ has a very small sensitivity ($=1$) and $\textrm{Lap}(\frac{1}{\epsilon_1})$ is very small. 
In \textsf{ARR}$_{\triangle}${}, we set the sampling probability $p_0$ to $p_0 = n^{-1/3}$ or $0.1n^{-1/3}$ so that the time complexity is $O(n^2)$. 

We ran each algorithm $20$ times and evaluated the average relative error over the $20$ runs. 

\subsection{Experimental Results}
\label{sub:results}

\noindent{\textbf{Relative Error vs. $\epsilon$.}}~~We first evaluated the relation between the relative error and 
$\epsilon$ in element DP or edge LDP, i.e., $2\epsilon$ in edge DP. 
We also measured the time to estimate the triangle/4-cycle count from the adjacency matrix $\mathbf{A}$ using a supercomputer \cite{ABCI} with two Intel Xeon Gold 6148 processors (2.40 GHz, 20 Cores) and 412 GB main memory. 

\begin{figure}[t]
  \centering
  \includegraphics[width=0.99\linewidth]{fig/res1_eps.pdf}
  \vspace{-4mm}
  \caption{Relative error vs. $\epsilon$ 
 
  ($n=107614$ in \textsf{Gplus}{}, $n=896308$ in \textsf{IMDB}{}, $c=1$). 
  $p_0$ is the sampling probability in the ARR. 
 
  }
  \label{fig:res1_eps}
\end{figure}


\begin{table}[t]
  \caption{Relative error (RE) when $\epsilon=0.5$ or $1$ and computational time ($n=107614$ in \textsf{Gplus}{}, $n=896308$ in \textsf{IMDB}{}, $c=1$). 
  The lowest relative error is highlighted in bold.
  }
  \vspace{-4mm}
  \centering
  (a) \textsf{Gplus}{}\\
  \begin{tabular}{|c|c|c|c|}
    \hline
    & RE ($\epsilon=0.5$) & RE ($\epsilon=1$) & Time (sec)\\ \hline
    \textsf{WShuffle}$_{\triangle}^*${} & $\textbf{0.298}$ & $\textbf{0.277}$ & $36.0$ \\ \hline
    \textsf{WShuffle}$_{\triangle}${} & $0.312$ & $0.279$ & $36.2$ \\ \hline
    \textsf{WLocal}$_{\triangle}${} & $6.10$ & $1.14$ & $58.3$ \\ \hline
    \textsf{ARR}$_{\triangle}${} ($p_0=n^{-1/3}$) & $49.0$ & $4.93$ & $715$ \\ \hline
    \textsf{ARR}$_{\triangle}${} ($p_0=0.1n^{-1/3}$) & $1.88 \times 10^3$ & $197$ & $34.8$ \\ \hline \hline
    \textsf{WShuffle}$_{\square}${} & $\textbf{0.145}$ & $\textbf{0.147}$ & $34.7$ \\ \hline
    \textsf{WLocal}$_{\square}${} & $2.08$ & $0.596$ & $57.0$ \\ \hline
  \end{tabular}\\
  (b) \textsf{IMDB}{}\\
  \begin{tabular}{|c|c|c|c|}
    \hline
    & RE ($\epsilon=0.5$) & RE ($\epsilon=1$) & Time (sec)\\ \hline
    \textsf{WShuffle}$_{\triangle}^*${} & $\textbf{0.488}$ & $\textbf{0.308}$ & $2.39 \times 10^3$ \\ \hline
    \textsf{WShuffle}$_{\triangle}${} & $1.41$ & $0.522$ & $2.40 \times 10^3$ \\ \hline
    \textsf{WLocal}$_{\triangle}${} & $74.6$ & $26.3$ & $3.96 \times 10^3$ \\ \hline
    \textsf{ARR}$_{\triangle}${} ($p_0=0.1n^{-1/3}$) & $2.98 \times 10^4$ & $3.27 \times 10^3$ & $2.81 \times 10^3$ \\ \hline \hline
    \textsf{WShuffle}$_{\square}${} & $\textbf{0.303}$ & $\textbf{0.308}$ & $2.29 \times 10^3$ \\ \hline
    \textsf{WLocal}$_{\square}${} & $282$ & $59.1$ & $3.91 \times 10^3$ \\ \hline
  \end{tabular}
  \label{tab:res1_eps_tri_time}
\end{table}

Figure~\ref{fig:res1_eps} shows the relative error ($c=1$). 
Here, we show the performance of \textsf{WShuffle}$_{\triangle}${} when we do not add the Laplacian noise (denoted by \textsf{WShuffle}$_{\triangle}${} (w/o Lap)). 
In \textsf{IMDB}{}, we do not show \textsf{ARR}$_{\triangle}${} with $p_0 = n^{-1/3}$, because it takes too much time (longer than one day). 
Table~\ref{tab:res1_eps_tri_time} highlights the relative error when $\epsilon=0.5$ or $1$. 
It also shows the running time of counting triangles or 4-cycles when $\epsilon=1$ (we verified that the running time had little dependence on $\epsilon$). 

Figure~\ref{fig:res1_eps} and Table~\ref{tab:res1_eps_tri_time} show that our shuffle algorithms dramatically improve the local algorithms. 
In triangle counting, 
\textsf{WShuffle}$_{\triangle}^*${} outperforms \textsf{WLocal}$_{\triangle}${} by one or two orders of magnitude and \textsf{ARR}$_{\triangle}${} by even more\footnote{Note that \textsf{ARR}$_{\triangle}${} uses only the lower-triangular part of the adjacency matrix $\mathbf{A}$ and therefore provides 
$\epsilon$-edge DP (rather than $2\epsilon$-edge DP); i.e., it does not suffer from the doubling issue explained in Section~\ref{sub:privacy}. However, Figure~\ref{fig:res1_eps} shows that \textsf{WShuffle}$_{\triangle}^*${} significantly outperforms \textsf{ARR}$_{\triangle}${} 
even if we double $\epsilon$ for only \textsf{WShuffle}$_{\triangle}^*${}.}. 
\textsf{WShuffle}$_{\triangle}^*${} also requires less running time than \textsf{ARR}$_{\triangle}${} with $p_0 = n^{-1/3}$. 
Although the running time of \textsf{ARR}$_{\triangle}${} can be improved by using a smaller $p_0$, it results in a higher relative error. 
In 4-cycle counting, 
\textsf{WShuffle}$_{\square}${} significantly outperforms \textsf{WLocal}$_{\square}${}. 
The difference between our shuffle algorithms and the local algorithms is larger in \textsf{IMDB}{} because it is more sparse; i.e., the difference between $d_{max}$ and $n$ is larger in \textsf{IMDB}{}. 
This is consistent with our theoretical results in Tables~\ref{tab:upper_bounds_triangle} and \ref{tab:upper_bounds_4cycle}. 

Figure~\ref{fig:res1_eps} and Table~\ref{tab:res1_eps_tri_time} also show that \textsf{WShuffle}$_{\triangle}^*${} outperforms \textsf{WShuffle}$_{\triangle}${}, especially when $\epsilon$ is small. 
This is because the variance is large when $\epsilon$ is small. 
In addition, \textsf{WShuffle}$_{\triangle}^*${} significantly outperforms \textsf{WShuffle}$_{\triangle}${} in \textsf{IMDB}{} because \textsf{WShuffle}$_{\triangle}^*${} significantly reduces the variance when $d_{max} \ll n$, as shown in Table~\ref{tab:upper_bounds_triangle}. 
In other words, this is also consistent with our theoretical results. 
For example, when $\epsilon=0.5$, our variance reduction technique reduces the relative error from $1.41$ to $0.488$ (about one-third) in \textsf{IMDB}{}. 

Furthermore, Figure~\ref{fig:res1_eps} shows that the relative error of \textsf{WShuffle}$_{\triangle}^*${} is hardly changed by adding the Laplacian noise. 
This is because the sensitivity of each user's degree $d_i$ is very small ($=1$). 
In this case, the Laplacian noise is also very small. 

\begin{figure}[t]
  \centering
  \includegraphics[width=0.99\linewidth]{fig/res2_n.pdf}
  \vspace{-4mm}
  \caption{Relative error vs. $n$ ($\epsilon=1$, $c=1$).
  }
  \label{fig:res2_n}
\end{figure}

Our \textsf{WShuffle}$_{\triangle}^*${} achieves a relative error of $0.3$ ($\ll 1$) 
when the privacy budget is $\epsilon = 0.5$ or $1$ in element DP ($2\epsilon = 1$ or $2$ in edge DP). 
\textsf{WShuffle}$_{\square}${} achieve a relative error of $0.15$ to $0.3$ with a smaller privacy budget (e.g., $\epsilon = 0.2$) because it  does not send local edges -- the error of \textsf{WShuffle}$_{\square}${} is mainly caused by user-pair sampling that is independent of $\epsilon$. 

In summary, our \textsf{WShuffle}$_{\triangle}^*${} and \textsf{WShuffle}$_{\square}${} significantly outperform the local algorithms and achieve a relative error much smaller than $1$ with a reasonable privacy budget, i.e., $\epsilon \leq 1$. 

\smallskip
\noindent{\textbf{Relative Error vs. $n$.}}~~Next, we evaluated the relation between the relative error and $n$. 
Specifically, we randomly selected $n$ users from all users and extracted a graph with $n$ users. 
Then we set $\epsilon = 1$ and changed $n$ to various values starting from $2000$. 

Figure~\ref{fig:res2_n} shows the results ($c=1$). 
When $n=2000$, \textsf{WShuffle}$_{\triangle}${} and \textsf{WShuffle}$_{\square}${} provide 
relative errors close to \textsf{WLocal}$_{\triangle}${} and \textsf{WLocal}$_{\square}${}, respectively. 
This is because the privacy amplification effect is limited when $n$ is small. 
For example, when $n=2000$ and $\epsilon=1$, 
the numerical bound is $\epsilon_L=1.88$. 
The value of $\epsilon_L$ increases with increase in $n$; e.g., when $n=107614$ and $896308$, the numerical bound is $\epsilon_L= 5.86$ and $7.98$, respectively. 
This explains the reason that our shuffle algorithms significantly outperform the local algorithms when $n$ is large in Figure~\ref{fig:res2_n}. 




\smallskip
\noindent{\textbf{Parameter $c$ in \textsf{WShuffle}$_{\triangle}^*${}.}}~~Finally, we evaluated our \textsf{WShuffle}$_{\triangle}^*${} while changing the parameter $c$ that controls the bias and variance. 
Recall that as $c$ increases, the bias is increased, and the variance is reduced. 
We set $\epsilon=0.1$ or $1$ and changed $c$ from $0.1$ to $4$. 


\begin{figure}[t]
  \centering
  \includegraphics[width=0.99\linewidth]{fig/res4_thr.pdf}
  \vspace{-4mm}
  \caption{Relative error vs. parameter $c$ in \textsf{WShuffle}$_{\triangle}^*${} ($n=107614$ in \textsf{Gplus}{}, $n=896308$ in \textsf{IMDB}{}).
  }
  \label{fig:res4_thr}
\end{figure}

Figure~\ref{fig:res4_thr} shows the results. 
Here, we also show the relative error of \textsf{WShuffle}$_{\triangle}${}. 
We observe that the optimal $c$ is different for $\epsilon=0.1$ and $\epsilon=1$. The optimal $c$ is around $3$ to $4$ for $\epsilon=0.1$, whereas the optimal $c$ is around $0.5$ to $1$ for $\epsilon=1$. 
This is because the variance of \textsf{WShuffle}$_{\triangle}${} is large (resp.~small) when $\epsilon$ is small (resp.~large). 
For a small $\epsilon$, a large $c$ is effective in significantly reducing the variance. 
For a large $\epsilon$, a small $c$ is effective in keeping a small bias. 

We also observe that \textsf{WShuffle}$_{\triangle}^*${} is always better than (or almost the same as) \textsf{WShuffle}$_{\triangle}${} when $c=1$ or $2$. 
This is because most users' degrees are smaller than the average degree $d_{avg}$, as described in Section~\ref{sub:var_red}. 
When $c=1$ or $2$, most user-pairs are ignored. 
Therefore, we can significantly reduce the variance at the cost of a small bias. 

\smallskip
\noindent{\textbf{Summary.}}~~In summary, our answers to the three questions at the beginning of Section~\ref{sec:experiments} are as follows. 
RQ1: Our \textsf{WShuffle}$_{\triangle}^*${} and \textsf{WShuffle}$_{\square}${} outperform the one-round local algorithms by one or two orders of magnitude (or even more). 
RQ2: Our variance reduction technique significantly reduces the relative error (e.g., by about one-third) 
for a small $\epsilon$ in a sparse dataset. 
RQ3: 
\textsf{WShuffle}$_{\triangle}^*${} achieves a relative error of $0.3$ ($\ll 1$) when $\epsilon=0.5$ or $1$ in element DP ($2\epsilon=1$ or $2$ in edge DP). 
\textsf{WShuffle}$_{\square}${} achieves a relative error of $0.15$ to $0.3$ with a smaller privacy budget: $\epsilon=0.2$. 

\section{Conclusion}
\label{sec:conclusion}
In this paper, we made the first attempt (to our knowledge) to apply the shuffle model to graph data. 
We proposed wedge shuffling as a basic technique and then applied it to 
one-round triangle and 4-cycle counting with several additional techniques. 
We showed upper bounds on the expected $l_2$ loss for each algorithm. 
We also showed through comprehensive experiments that our one-round shuffle algorithms significantly outperform the one-round local algorithms and achieve a small relative error with a reasonable privacy budget, e.g., smaller than $1$ in edge DP. 

For future work, we would like to apply wedge shuffling to other subgraphs such as 3-hop paths \cite{Sun_CCS19} and $k$-triangles \cite{Karwa_PVLDB11}. 


\bibliographystyle{ACM-Reference-Format}
