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

This file is documentation for the SIAM \LaTeX\ style, including how
to typeset the main document, the {\scshape Bib}\TeX\xspace\ file, and any supplementary
material. More information
about SIAM's editorial style can be found in the style manual, available
at \url{https://www.siam.org/journals/pdf/stylemanual.pdf}.
The major changes in the SIAM standard class are summarized in \cref{sec:changes}.
The SIAM \LaTeX\@ files can be found at
\url{https://www.siam.org/journals/auth-info.php}. The files that
are distributed for the standard macros are given below. 
\begin{itemize}
\item \texttt{siamart220329.cls} (required): Main SIAM standard \LaTeX\ class file.
\item \texttt{siamplain.bst} (required): Bibliographic style file for
  {\scshape Bib}\TeX\xspace.
\item \texttt{docsiamart.tex}: Produces this documentation.
\item \texttt{references.bib}: {\scshape Bib}\TeX\xspace\ database for this
  documentation and examples.
\item \texttt{ex\_article.tex}: Template for article.
\item \texttt{ex\_supplement.tex}: Template for supplement.
\item \texttt{ex\_shared.tex}: Template for shared information for
  article and supplement.
\end{itemize}
To use these files, put \texttt{siamart220329.cls} and
\texttt{siamplain.bst} in the directory with your
paper or, alternatively, into your \LaTeX\@ and {\scshape Bib}\TeX\xspace\@ paths,
respectively. 
The outline of a SIAM \LaTeX\ article is shown in
\cref{ex:outline}. Templates are provided and discussed in more detail
in \cref{sec:template}.

\begin{example}[label={ex:outline},listing only,%
  listing options={style=siamlatex,{morekeywords=[1]{maketitle},
      morekeywords=[2]{siamart220329}},}]%
  {Document outline}
\documentclass{siamart220329}
\begin{document}
\maketitle
\end{document}
\end{example}

\section{Class options}
\label{sec:class-options}

Class options can be included in the bracketed argument of the
command, separated by commas. The possible class options are:
\begin{itemize}
\item \code{review} --- Recommended for submitting your manuscript to
  a SIAM journal.  Adds line numbers as well as the statement ``This
  manuscript is for review purposes only'' to the bottom of each page.
\item \code{final} --- Turns off the black boxes that help authors
  identify lines that are too long. The final published version will
  have this option on.
\item \code{supplement} --- Specifies that the file is a supplement
  and not the main document, causing changes in the appearance of the
  title and numbering; see \cref{sec:supplement} for details.
\item \code{hidelinks} --- Turns off colors on hyperlinks; see
  \cref{sec:cr+hyp}. The hyperlinks still exist, but there is no color
  to differentiate them.
  The final published version will have this option on.
\end{itemize}


\section{Front matter}
\label{sec:front}
The title and author parts are formatted using the standard
\code{\title}, \code{\author}, and \code{\maketitle} commands as
described in Lamport \cite{La86}. The title and author should be
declared in the preamble. The title and author names are automatically
converted to uppercase in the document.
If there is more than one author, each additional author should be preceded by
the \code{\and} command.  
The addresses and support acknowledgments are added via
\code{\thanks}. Each author's thanks should specify their address.
The support acknowledgment should be put in the title thanks,
unless specific support needs to be specified for individual authors,
in which case it should follow the author address.  
The header for this file was produced by the code in \cref{ex:header}, including an
example of a shared footnote. Each thanks produces a footnote, so the
footnote of the second author is \#3.
The command
\code{\headers{title}{authors}} command, with
the title (possibly shortened to fit) and the authors' names, creates
the page headers, automatically converted to uppercase.

\examplefile[label={ex:header},listing only,%
listing options={style=siamlatex,%
deletetexcs={and,thanks,title,author},%
{moretexcs=[2]{and,thanks,title,author,maketitle,headers,email}}}
]{Title and authors in preamble}{tmp_\jobname_header.tex}

\newpage
Following the author and title is the abstract, key words listing, and MSC codes, designated using the \code{abstract}, \code{keywords}, and \code{MSCcodes}
environments. 
Authors are responsible for providing MSC codes which can be found
on the AMS web site \cite{AMSMSC2010}. The abstract, keywords, and MSC codes
for this document are specified in \cref{ex:abstract}.

\examplefile[label={ex:abstract},%
before upper={\preamble{\bs newcommand\{\bs BibTeX\}\{\{\bs scshape Bib\}\bs TeX\bs xspace\}}},
listing only,%
listing options={style=siamlatex,%
{morekeywords=[2]{abstract,keywords,MSCcodes}}}
]{Abstract, keywords, and MSC codes}{tmp_\jobname_abstract.tex}

A more complete example, including a PDF supplement, that uses the
included files \texttt{ex\_article.tex}, \texttt{ex\_supplement.tex},
and \texttt{ex\_shared.tex} is discussed in
\cref{sec:template}. The example files can be used as a starting point
for producing a document.

\section{Cross references and hyperlinks}
\label{sec:cr+hyp}

SIAM now supports cross references and hyperlinks via the
\texttt{cleveref} and \texttt{hyperef} packages, which are loaded by
the class file.

\subsection{Cleveref}
\label{sec:cleveref}

SIAM strongly recommends using the commands provided by 
the \texttt{cleveref} package for cross referencing.
The package is automatically loaded and already customized
to adhere to SIAM's style guidelines. 
To create a cross reference, use the command \code{\cref} (inside
sentence) or \code{\Cref} (beginning of a sentence) in place of the
object name and \code{\ref}. 
The \texttt{cleveref}
package enhances \LaTeX's cross-referencing features, allowing
the format of cross references to be determined automatically
according to the ``type" of cross reference (equation, section, etc.)
and the context in which the cross reference is used.
So, the package
\emph{automatically} inserts the object name as well as the
appropriate hyperlink; see \cref{ex:cref}.
It may require two \LaTeX\@ compilations for the references to show up
correctly.
Additional examples are shown in the sections below
for equations, tables, figures, sections, etc.

\begin{example}[label=ex:cref,bicolor,listing options={style=siamlatex,%
{morekeywords=[2]{cref,ref}}}]{Advantage of using cleveref}
The normal way to get a cross reference with a hyperlink requires a
lot of typing: \hyperref[thm:mvt]{Theorem~\ref*{thm:mvt}}.
The \texttt{cleveref} package gets both the name and hyperlink
automatically using a single macro: \cref{thm:mvt}.
It also handles multiple references with the same macro, such as
\cref{thm:mvt,fig:pgfplots,fig:testfig}.
\end{example}


\subsection{Hyperef}
\label{sec:hyperef}

Hyperlinks are created with the \code{\href} and \code{\url} commands,
as shown in \cref{ex:href}.
SIAM has also defined the \code{\email} command, as shown in
\cref{ex:header}.
You can hide links (i.e., turn off link colors) with the \code{hidelinks} option.

\begin{example}[label={ex:href},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{href,url}}}]{Creating hyperlinks}
The \href{https://www.siam.org}{SIAM homepage} has general information.
Note that the colored text will \emph{not} appear in the print version
nor will the hyperlink be active, so the writer may want to specify
the location explicitly instead by using \url{https://www.siam.org}.
\end{example}

Note that homepage links via \code{\url} in the \code{\thanks} environment require
special formatting for the tilde (\string~) character. The formatting
is used in the template and shown in \cref{ex:shared}.


\section{Math and equations}
\label{sec:math}
Here we show some example equations, with numbering, and examples of
referencing the equations. SIAM now includes the package 
\texttt{amsmath} by default, and we include some of its features as
well, although the reader should consult the package user manual for
further guidance \cite{amsmath,shortmath}.
 Several of the example are adapted 
from Mittlebach and Goossen's guide to \LaTeX~\cite{MiGo04}.

\Cref{ex:textmath} is a straightforward example of inline mathematics
equations that does not use any special packages or features.

\begin{example}[label={ex:textmath},bicolor]{Inline math}
The following shows an example of math in text: 
Let $S=[s_{ij}]$ ($1\leq i,j\leq n$) be a $(0,1,-1)$-matrix of order $n$.
\end{example}

In \cref{ex:bbm}, we show the recommended method for getting
blackboard fonts using the \texttt{amsfonts} package. This is not
loaded by default and must be included in the preamble.

\begin{example}[label={ex:bbm},bicolor,before upper={\preamble{\bs
      usepackage\{amsfonts\}}},%
listing options={style=siamlatex,%
{morekeywords=[2]{mathbb}}}]{Blackboard math}
Blackboard bold characters, such as $\mathbb{C}$ and $\mathbb{R}$,
should be created with the \texttt{amsfonts} package, although this
is not included by default.
\end{example}

\Cref{ex:smallmatrix} shows the \code{smallmatrix} environment for an
inline matrix from the \texttt{amsmath} package, which is included by
default.

\begin{example}[label={ex:smallmatrix},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{smallmatrix}}}]{Inline matrix}
Matrices of no more than two rows appearing in text can be created
as shown in the next example:
$B = \bigl[ \begin{smallmatrix} B_{11} & B_{12} \\ 
  B_{21} & B_{22} \end{smallmatrix} \bigr]$.
\end{example}

Bigger matrices can be rendered with environments from
the \texttt{amsmath} package, such as \code{bmatrix} and
\code{pmatrix} used in
\cref{ex:matrices}.

\begin{example}[label={ex:matrices},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{bmatrix,pmatrix}}}]{Creating matrices}
Display matrices can be rendered using environments from \texttt{amsmath}:
\begin{equation}\label{eq:matrices}
  S=\begin{bmatrix}1&0\\0&0\end{bmatrix}
  \quad\text{and}\quad
  C=\begin{pmatrix}1&1&0\\1&1&0\\0&0&0\end{pmatrix}.
\end{equation}
\Cref{eq:matrices} shows some example matrices.
\end{example}

\newpage
\Cref{ex:dmo} shows how to use the \code{\DeclareMathOperator} command
from the \texttt{amsopn} package to declare the \code{\Range} macro. 
(This example also uses the \texttt{braket} package for the
\code{\set} macro, but this is not necessarily recommended by SIAM.)
\begin{example}[label={ex:dmo},%
before upper={\preamble{\bs usepackage\{braket,amsfonts,amsopn\}}\\
\noindent\preamble{\bs DeclareMathOperator\{\bs Range\}\{Range\}}},%
bicolor,%
listing options={style=siamlatex,%
{moretexcs=[2]{Range}}}
]{Declaring math operators}
An example of a math operator:
\begin{equation}\label{eq:range}
  \Range(A) = \set{ y \in \mathbb{R}^n | y = Ax }.
\end{equation}
\end{example}

\Cref{ex:foo} shows how to use the \code{align} environment from
\texttt{amsmath} to easily align multiple equations.

\begin{example}[label={ex:foo},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{align}}}]{Aligned equations}
\Cref{eq:a,eq:b,eq:c} show three aligned equations.
\begin{align} 
  f &= g, \label{eq:a} \\
  f' &= g', \quad\text{and} \label{eq:b} \\
  \mathcal{L}f &= \mathcal{L}g \label{eq:c}.
\end{align}
\end{example}

Another way to number a set of equations is the
\code{subequations} environment from \texttt{amsmath}, as shown in \cref{ex:aligned}.

\begin{example}[label={ex:aligned},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{subequations}}}]{Subequations}
We calculate the Fr\'{e}chet derivative of $F$ as follows:
\begin{subequations}
\begin{align}
  F'(U,V)(H,K) 
  &= \langle R(U,V),H\Sigma V^{T} + U\Sigma K^{T} -
  P(H\Sigma V^{T} + U\Sigma K^{T})\rangle \label{eq:aa} \\
  &= \langle R(U,V),H\Sigma V^{T} + U\Sigma K^{T}\rangle 
  \nonumber \\
  &= \langle R(U,V)V\Sigma^{T},H\rangle + 
  \langle \Sigma^{T}U^{T}R(U,V),K^{T}\rangle. \label{eq:bb}
\end{align}
\end{subequations}
\Cref{eq:aa} is the first line, and \cref{eq:bb} is the last line.
\end{example}

~

For an equation split over multiple lines, \cref{ex:ml} shows the
usage of the \code{multline} environment provided by \texttt{amsmath}.

~

\begin{example}[label={ex:ml},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{multline}}}]{Equation split across lines}
We claim that the projection $g(U,V)$ is given by the pair of matrices:
\begin{multline} \label{eq:ml}
  g(U,V) = \biggl( \frac{R(U,V)V\Sigma^{T}U^{T} 
    - U\Sigma V^{T}R(U,V)^{T}}{2}U,\\
  \frac{R(U,V)^{T}U\Sigma V^{T}-V \Sigma^{T}U^{T}R(U,V)}{2}V \biggr).
\end{multline}
\end{example}

\section{Theorem-like environments}
\label{sec:thm}
SIAM loads \texttt{ntheorem} package and uses it to define the following
theorem-like environments:
\code{theorem}, 
\code{lemma}, 
\code{corollary},
\code{definition}, and
\code{proposition}.
SIAM also defines a \code{proof} environment that automatically
inserts the symbol ``$\,\proofbox\,$'' at the end of any proof, even if it ends in an
equation environment. \emph{Note that the document may need to be
  compiled twice for the mark to appear.}
Some of the calculus examples were adapted from \cite{CalcI}. 
\Cref{ex:theorem} shows usage of the \code{theorem} environment. 
An optional argument can be used to name the theorem.
\Cref{ex:cor} illustrates a corollary, without a name, and the
proof environment.

~

\begin{example}[label=ex:theorem,bicolor,parbox=false,%
listing options={style=siamlatex,%
{morekeywords=[2]{theorem}}}]{Theorem}
\begin{theorem}[Mean Value Theorem]\label{thm:mvt}
  Suppose $f$ is a function that is continuous on the closed interval
  $[a,b]$.  and differentiable on the open interval $(a,b)$.
  Then there exists a number $c$ such that $a < c < b$ and
  \begin{displaymath}
    f'(c) = \frac{f(b)-f(a)}{b-a}.
  \end{displaymath}
  In other words, $f(b)-f(a) = f'(c)(b-a)$.
\end{theorem}
\end{example}

\begin{example}[label=ex:cor,bicolor,parbox=false,%
listing options={style=siamlatex,%
{morekeywords=[2]{corollary,proof}}}]%
{Corollary and proof}
\begin{corollary}
  Let $f(x)$ be continuous and differentiable everywhere. If $f(x)$
  has at least two roots, then $f'(x)$ must have at least one root.
\end{corollary}
\begin{proof}
  Let $a$ and $b$ be two distinct roots of $f$.
  By \cref{thm:mvt}, there exists a number $c$ such that
  \begin{displaymath}
    f'(c) = \frac{f(b)-f(a)}{b-a} = \frac{0-0}{b-a} = 0.
  \end{displaymath}
\end{proof}
\end{example}

SIAM also defines commands to create your own theorem- and
remark-like environments:
\begin{itemize}
\item \code{newsiamthm} --- Small caps header, italized body. 
\item \code{newsiamremark} --- Italics header, roman body.
\end{itemize}
Each command takes two arguments. The first is the environment name,
and the second is the name to show in the document. These commands should
be used instead of \code{\newtheorem}.
\Cref{ex:claim,ex:ref} shows how to use the commands above, including how to
specify the plural version for \texttt{cleveref} if it is unusual.

\begin{example}[label=ex:claim,bicolor,%
before upper={\preamble{\bs newsiamthm\{claim\}\{Claim\}}\\
\noindent\preamble{\bs newsiamremark\{hypothesis\}\{Hypothesis\}}\\
\noindent\preamble{\bs crefname\{hypothesis\}\{Hypothesis\}\{Hypotheses\}}},%
parbox=false,%
listing options={style=siamlatex,%
{morekeywords=[2]{claim,proof,hypothesis}}}]{New theorem-like environment}
\begin{claim}\label{cl:constant}
  If $f'(x) = 0$ for all $x \in (a,b)$ then $f(x)$ is constant on $(a,b)$.
\end{claim}
\begin{hypothesis}\label{hyp1}
The function $f$ is continuously differentiable.
\end{hypothesis}
\begin{hypothesis}\label{hyp2}
The random variable is normally distributed.
\end{hypothesis}
\end{example}


\begin{example}[label=ex:ref,bicolor,listing options={style=siamlatex,%
{morekeywords=[2]{cref}}}]{References}
We can reference multiple types of objects with a single reference:
\cref{cl:constant,thm:mvt,hyp1,hyp2}.
\end{example}


\section{Tables}
\label{sec:tab}
Table captions should go above the tables.
\Cref{ex:simpletable} shows the code to generate a
\cref{tab:simpletable}.
A more complicated example is shown in 
\cref{ex:table}, which generates \cref{tab:KoMa14}. This
example uses subfloats via the \texttt{subfig} package, as well as
special column options from the \texttt{array} package.

\begin{tcbverbatimwrite}{tmp_\jobname_simpletable.tex}
\begin{table}[tbhp]
\footnotesize
  \caption{Example table}\label{tab:simpletable}
\begin{center}
  \begin{tabular}{|c|c|c|} \hline
   Species & \bf Mean & \bf Std.~Dev. \\ \hline
    1 & 3.4 & 1.2 \\
    2 & 5.4 & 0.6 \\ \hline
  \end{tabular}
\end{center}
\end{table}
\end{tcbverbatimwrite}

\examplefile[label={ex:simpletable},%
listing only, listing options={style=siamlatex}]%
{Example table.}{tmp_\jobname_simpletable.tex}

\input{tmp_\jobname_simpletable.tex}

\begin{tcbverbatimwrite}{tmp_\jobname_table.tex}
\newcolumntype{R}{>{$}r<{$}}  %
\newcolumntype{V}[1]{>{[\;}*{#1}{R@{\;\;}}R<{\;]}} %
\begin{table}[tbhp]
\footnotesize
 \captionsetup{position=top}
  \caption{Example table adapted from Kolda and Mayo \rm{\cite{KoMa14}}.}\label{tab:KoMa14}
\begin{center}
  \subfloat[$\beta=1$]{
    \begin{tabular}{|r|R|V{3}|c|r@{\,$\pm$\,}l|} \hline
occ. & \multicolumn{1}{c|}{$\lambda$} & \multicolumn{4}{c|}{$\mathbf{x}$} & 
fevals &  \multicolumn{2}{c|}{time (sec.)}\\ \hline
 718 & 11.3476 &  0.5544 &  0.3155 &  1.2018 &  0.0977 &   45 & 0.17 & 0.06 \\ \hline 
 134 &  3.7394 &  0.2642 & -1.1056 &  0.2657 & -0.3160 &   31 & 0.12 & 0.05 \\ \hline 
   4 & \multicolumn{6}{c|}{\emph{--- Failed to converge ---}} & 0.21 & 0.10 \\ \hline 
 \end{tabular}}

  \subfloat[$\beta=-1$]{
    \begin{tabular}{|r|R|V{3}|c|r@{\,$\pm$\,}l|} \hline
occ. & \multicolumn{1}{c|}{$\lambda$} & \multicolumn{4}{c|}{$\mathbf{x}$} & 
fevals &  \multicolumn{2}{c|}{time (sec.)}\\ \hline
  72 & -1.1507 &  0.2291 &  0.6444 &  0.3540 & -0.8990 &   34 & 0.14 & 0.06 \\ \hline 
 624 & -6.3985 &  0.1003 &  0.1840 &  0.5305 &  1.2438 &   48 & 0.19 & 0.08 \\ \hline 
   2 & \multicolumn{6}{c|}{\emph{--- Failed to converge ---}} & 0.23 & 0.02 \\ \hline 
 \end{tabular}}
\end{center}
\end{table}
\end{tcbverbatimwrite}

\examplefile[label={ex:table},%
before upper={\preamble[\scriptsize]{\bs usepackage\{array\}}\\[-0.4em]
\noindent\preamble[\scriptsize]{\bs usepackage[caption=false]\{subfig\}}},%
listing only, listing options={%
  style=siamlatex,basicstyle=\ttfamily\scriptsize}]%
{Example table with subtables.}{tmp_\jobname_table.tex}

\input{tmp_\jobname_table.tex}

\section{Figures}
\label{sec:fig}

It is recommended that all figures be generated in high resolution.
In the past, SIAM has required encapsulated
postscript (EPS) format for final production. This is still an
acceptable format, but SIAM also now allows high-resolution 
PDF, JPEG, and PNG figures. 
If working with EPS images and using \texttt{pdflatex}, we
recommend the package \texttt{epstopdf} to automatically convert EPS
images to PDF for inclusion in PDF documents created by
\texttt{pdflatex}.  Please note that \texttt{epstopdf} requires \texttt{texlive-font-utils} which may not be part of a standard \LaTeX\ installation.
\Cref{ex:fig} shows the code to generate \cref{fig:testfig}.  This
example uses the \texttt{graphicx} package for the
\code{\includegraphics} command. 


\begin{tcbverbatimwrite}{tmp_\jobname_fig.tex}
\begin{figure}[tbhp]
  \centering
  \subfloat[$\epsilon_{\max}=5$]{\label{fig:a}\includegraphics{lexample_fig1}}
  \subfloat[$\epsilon_{\max}=0.5$]{\label{fig:b}\includegraphics{lexample_fig2}}
  \caption{Example figure using external image files.}
  \label{fig:testfig}
\end{figure}
\end{tcbverbatimwrite}

\examplefile[label={ex:fig},%
before upper={\preamble[\scriptsize]{\bs usepackage\{graphicx,epstopdf\}}\\[-0.4em]
\noindent\preamble[\scriptsize]{\bs usepackage[caption=false]\{subfig\}}},%
listing only, listing options={%
  style=siamlatex,basicstyle=\ttfamily\scriptsize}]%
{Example figure with subfigures and external files}{tmp_\jobname_fig.tex}

\input{tmp_\jobname_fig.tex}

Another option for figures is a graphics-generator that is platform- and
format-independent.
PGF is a TeX macro package for generating  such graphics and works together with the most important TeX backend drivers, including pdftex and dvips. 
The user-friedly syntax layer called TikZ.
Here we show an example using \texttt{PGFPLOTS}, useful for drawing high-quality 
plots directly in \LaTeX.
\Cref{ex:data} and \cref{ex:pgfplots} shows the data and code, respectively, to generate \cref{fig:pgfplots}, adapted from \cite{pgfplots}.

\examplefile[label={ex:data},listing only,
listing options={style=siamlatex,basicstyle=\ttfamily\scriptsize}]%
{Example data file (data.dat)}{data.dat}

\begin{tcbverbatimwrite}{tmp_\jobname_tikz.tex}
\begin{figure}[tbhp]
  \centering
  \begin{tikzpicture}
    \begin{loglogaxis}[height=2.75in, grid=major,
      xlabel={Degrees of Freedom}, ylabel={$L_2$ Error},
      legend entries={$d=2$,$d=3$}]
      \addplot table [x=d2_dof,y=d2_l2_err] {data.dat};
      \addplot table [x=d3_dof,y=d3_l2_err] {data.dat};
    \end{loglogaxis}
  \end{tikzpicture}
  \caption{Example \texttt{PGFPLOTS} figure.}
  \label{fig:pgfplots}
\end{figure}
\end{tcbverbatimwrite}

\examplefile[label={ex:pgfplots},%
before upper={\preamble[\scriptsize]{\bs usepackage\{pgfplots\}}},%
listing only, listing options={%
  style=siamlatex}]%
{Example TikZ/PGF for platform-independent graphics.}{tmp_\jobname_tikz.tex}

\input{tmp_\jobname_tikz.tex}


\section{Algorithms}
\label{sec:algs}

SIAM automatically includes the \texttt{algorithm} package in the
class definition. This provides the float environment.
Users have the choice of \texttt{algpseudocode},
\texttt{algorithmic}, and other packages for actually formatting the algorithm.
For example, \cref{alg:buildtree} is produced by the code in \cref{ex:alg}.
In order to reference lines within the algorithm, we need to tell the
\texttt{cleveref} package how to do the referencing, which is the
second line of \cref{ex:alg}. Then we can use the code
\code{\cref{line3}} to produce \cref{line3}.

\begin{tcbverbatimwrite}{tmp_\jobname_alg.tex}
\begin{algorithm}
\caption{Build tree}
\label{alg:buildtree}
\begin{algorithmic}[1]
\STATE{Define $P:=T:=\{ \{1\},\ldots,\{d\}$\}}
\WHILE{$\#P > 1$}
\STATE\label{line3}{Choose $C^\prime\in\mathcal{C}_p(P)$ with $C^\prime := \operatorname{argmin}_{C\in\mathcal{C}_p(P)} \varrho(C)$}
\STATE{Find an optimal partition tree $T_{C^\prime}$ }
\STATE{Update $P := (P{\setminus} C^\prime) \cup \{ \bigcup_{t\in C^\prime} t \}$}
\STATE{Update $T := T \cup \{ \bigcup_{t\in\tau} t : \tau\in T_{C^\prime}{\setminus} \mathcal{L}(T_{C^\prime})\}$}
\ENDWHILE
\RETURN $T$
\end{algorithmic}
\end{algorithm}
\end{tcbverbatimwrite}

\examplefile[float=htpb,label={ex:alg},%
before upper={\preamble[\scriptsize]{\bs usepackage\{algorithmic\}}\\[-0.4em]
\preamble[\scriptsize]{\bs Crefname\{ALC@unique\}\{Line\}\{Lines\}}},%
listing only, listing options={%
  style=siamlatex,basicstyle=\ttfamily\scriptsize}]%
{Example algorithm}{tmp_\jobname_alg.tex}

\input{tmp_\jobname_alg.tex}

\section{Sections}
\label{sec:sec}
Sections are denoted using standard \LaTeX\ section commands, i.e.,
\code{\section}, \code{\subsection}, etc. 
If you wish to end the section title with something other that a period (the default), you have to add the command \code{\nopunct} at the end of the title.

Appendices are created with the normal sectioning commands, following
the command \code{
\section{Introduction} \label{sec:introduction}
We consider distributed optimization (DO) where multiple agents cooperate to minimize a global objective function, expressed as a sum of local objective functions, subject to some constraints. 
In this paper we focus on the following DO model:
\begin{subequations}
\label{model:dist}
\begin{align}
\min_{w \in \mathbb{R}^n} \ &  \sum_{p=1}^P f_p(w) \\
\mbox{s.t.} \ 
& w \in \mathcal{W} := \cap_{p=1}^P \mathcal{W}_p, 
\end{align}
\end{subequations}
where $P$ is the number of agents, $w \in \mathbb{R}^n$ is a global variable vector, $\mathcal{W}_p$ is a local feasible region assumed to be a compact convex set for each agent $p$, and $f_p$ is a local convex objective function. 
The model \eqref{model:dist} appears in many applications, including distributed control in power systems \cite{molzahn2017survey, patari2021distributed} and distributed machine learning (ML) as well as federated learning (FL) \cite{li2020federated, 9599369}.
Numerous DO algorithms \cite{yang2019survey} have been developed, but most of the algorithms are susceptible to data leakage. 
Specifically, each agent iteratively solves a local optimization model constructed by its own data and communicates some information (e.g., local solutions or local gradients) with its neighbors in a communication network until a global solution $w^*$ is obtained.
Even though locally stored data are not shared with other agents, it is still possible to reconstruct the data from the information communicated among agents, as pointed out in \cite{shokri2017membership, zhu2019deep}.

To address the data privacy concern, we propose a privacy-preserving DO algorithm for solving \eqref{model:dist}.
The proposed algorithm provides a statistical guarantee of data privacy, known as differential privacy (DP), and a sequence of iterates that converges to an optimal solution of \eqref{model:dist} in expectation. 
Under the existence of convex constraints (i.e., $\mathcal{W}_p \subset \mathbb{R}^n$), we show that the proposed algorithm provides more accurate solutions compared with the existing DP algorithm, while both algorithms provide the same level of data privacy.
As a result, our algorithm can mitigate a trade-off between data privacy and solution accuracy (i.e., learning performance in the context of ML), which is one of the main challenges in developing DP algorithms, as described in \cite{dwork2014algorithmic}.

\subsection{Related literature}
DP is a privacy-preserving technique that randomizes an output of a query on data such that any single data sample cannot be inferred by the randomized output \cite{dwork2006calibrating,dwork2014algorithmic}.
Typically, the randomized output is constructed by adding either Laplacian noise or Gaussian noise to the true output of the query, known as Laplace and Gaussian mechanisms, respectively.
Because the true output is perturbed by adding the noise, this DP technique is often referred to as an output perturbation method.
We will discuss details of DP in Section \ref{sec:prelim_dp} since the algorithm proposed in this paper guarantees data privacy in an DP manner.

DP has been integrated in various learning algorithms that train ML models by using data centralized at a single server.
Depending on where to inject noise to randomize an output of the ML models, DP can be categorized by output and objective perturbation methods.
For ease of exposition, we denote the output and objective perturbation by \texttt{OutPert} (i.e., adds a noise vector to the model parameters) and \texttt{ObjPert} (i.e., adds a noisy linear function to the objective function before training), respectively.
These two methods have been applied to convex empirical risk minimization models in \cite{chaudhuri2011differentially, kifer2012private}, where the authors show that both theoretically and empirically, \texttt{ObjPert} is superior to \texttt{OutPert} in managing the inherent trade-off between privacy and learning performance. 

Within the context of distributed ML as well as FL, DP has been introduced in numerous algorithms \cite{agarwal2018cpsgd, naseri2020local, wei2020federated, zhang2016dynamic, huang2019dp, huang2020differentially} and frameworks \cite{ryu2022appfl,he2020fedml,ziller2021pysyft} mainly because of its computational efficiency \cite{wang2020hybrid}.
Specifically, local ML models solved for every iteration of the algorithms are randomized to communicate the randomized outputs (i.e., local model parameters or local gradients depending on the algorithms), which cannot be utilized for reconstructing any single data sample in the local dataset.
Most DP algorithms in the literature rely on \texttt{OutPert} to ensure DP \cite{naseri2020local,wei2020federated,huang2019dp,huang2020differentially}, which adds noise to the intermediate outputs of local ML models.
A few exceptions include \cite{zhang2016dynamic}, where \texttt{ObjPert} is applied to the alternating direction method of multipliers (ADMM) algorithm to randomize the outputs of local models to ensure DP.
For the unconstrained convex ML models, the authors show that ADMM with \texttt{ObjPert} outperforms that with \texttt{OutPert}.
However, the ADMM algorithm requires each agent to compute an exact solution of the convex subproblem for every iteration, which can be computationally expensive.
This computational challenge has been circumvented by employing a linearized ADMM \cite{huang2019dp}, where the convex loss function is linearly approximated, resulting in solving a convex quadratic subproblem that has a closed-form solution.
This work has been extended by introducing multiple local updates to reduce communication costs in \cite{huang2020differentially} and shown to outperform the other existing DP algorithms.

The aforementioned DP algorithms have been developed for unconstrained models; but, in general, constraints are necessary, for example, in most optimal control problems including distributed control of power flow. 
Also, imposing hard constraints on ML models (instead of penalizing constraints in the objective function as in \cite{shen2021agnostic,du2021fairness}) is increasingly considered for purposes such as improving accuracy, explaining decisions suggested by ML models, promoting fairness, and observing some physical laws (e.g., \cite{gallego2022controlled, goh2016satisfying, marquez2017imposing, zafar2019fairness, zhao2019physics}).
This calls for developing DP algorithms suitable for the general constrained optimization model, like the form \eqref{model:dist}.
Indeed the state-of-the-art DP algorithm \cite{huang2020differentially} developed for unconstrained models can be directly applied to solve \eqref{model:dist}, but \texttt{OutPert} used in the algorithm could make each agent  communicate \emph{infeasible} randomized outputs of the constrained local optimization models for every iteration of the algorithm, which in turn could negatively affect the overall convergence of the algorithm and the quality of the final solution.
To address this infeasibility issue arising from the existing DP algorithms applied to the general constrained optimization models, 
we propose a new DP algorithm based on \texttt{ObjPert} for solving \eqref{model:dist}, which allows each agent to communicate \emph{feasible} randomized outputs for every iteration of the algorithm while ensuring data privacy.
Specifically, the proposed algorithm generalizes the existing linearized ADMM by introducing a multiple local updates technique to reduce communication costs and incorporating \texttt{ObjPert} into local constrained optimization models to communicate feasible randomized ouputs that ensure DP. 
To the best of our knowledge, no previous work studies the privacy and convergence analyses of \texttt{ObjPert} in the linearized ADMM with multiple local updates for constrained convex distributed optimization models, which are included in this paper.




 
\subsection{Contributions} 
Our major contribution is to provide  theoretical and numerical analyses of \texttt{ObjPert} in the linearized ADMM with multiple local updates for the constrained model of the form \eqref{model:dist}.
First, under the existence of convex constraints $\mathcal{W}_p$ defined for every agent $p$, we propose how to perturb the local objective function $f_p$ such that the resulting outputs of the local model cannot be utilized for reconstructing data, thus preserving data privacy in an DP manner. 
Specifically, we propose Laplace and Gaussian mechanisms within the context of the linearized ADMM with \texttt{ObjPert}.
It turns out that the proposed mechanisms are different from the existing Laplace and Gaussian mechanisms used in \texttt{OutPert} with respect to the sensitivity computation.
Second, we show that the sequence of iterates produced by our algorithm converges to an optimal solution of \eqref{model:dist} in expectation and that its associate convergence rate is sublinear. 
Specifically, the rate is $O(\frac{1}{\bar{\epsilon} E \sqrt{T}})$ for a smooth convex function, $O(\frac{1}{\bar{\epsilon}^2 E  \sqrt{T}})$ for a nonsmooth convex function, and $O(\frac{1}{\bar{\epsilon}^2 E T})$ for a strongly convex function, where $T$ is the number of ADMM rounds, $E$ is the number of local updates, and $\bar{\epsilon}$ is a privacy budget parameter resulting in stronger data privacy as $\bar{\epsilon}$ decreases.
Third, we numerically demonstrate the effectiveness of our DP algorithm in two different application areas:  distributed control of power flow and federated learning.
 
\subsection{Notation and organization of the paper}
We use $\langle \cdot, \cdot \rangle$ and $\|\cdot \|$ to denote the scalar product and the Euclidean norm, respectively.
We define $[a]:= \{1,\ldots,a\}$.
We use $f'$ and $\nabla f$ to denote a subgradient and a gradient of $f$, respectively.
We use \textbf{pdf}, \textbf{det}, and \textbf{relint} to stand for probability density function, determinant, and relative interior, respectively.

We present some preliminaries in Section \ref{sec:preliminaries}, the proposed algorithm in Section \ref{sec:dpadmm}, privacy analysis in Section \ref{sec:privacy}, convergence analysis in Section \ref{sec:convergence}, numerical demonstration of the proposed algorithm in Section \ref{sec:experiments}, and conclusions in Section \ref{sec:conclusion}. 

\section{Preliminaries} \label{sec:preliminaries}
In this section we present some preliminaries on ADMM and DP for the completeness of this paper. 

\subsection{ADMM} \label{sec:ADMM}
Numerous distributed optimization algorithms \cite{liu2017convergence, nedic2009distributed, johansson2008subgradient, boyd2011distributed, shi2014linear} have been developed for solving \eqref{model:dist}. 
Among them, ADMM and its variants have been widely used for solving distributed optimization. 
Since the proposed algorithm in this paper is built on ADMM, we briefly describe a classical ADMM.
To this end, we first reformulate \eqref{model:dist} by introducing a local copy $z_p \in \mathbb{R}^{n}$ for every agent $p \in [P]$:
\begin{subequations}
\label{model:dist_1}  
\begin{align}
  \min_{w \in \mathbb{R}^n} \ &  \sum_{p=1}^P f_p(z_p)  \\
  \mbox{s.t.} \ 
  & z_p \in \mathcal{W}_p, \ \forall p \in [P], \\
  & w  = z_p, \ \forall p \in [P], \label{dist_1_consensus}
\end{align}
\end{subequations}
where $w$ represents global variables and $z_p$ represents local variables for agent $p$.

A typical iteration $t$ of ADMM for solving \eqref{model:dist_1} is given by
\begin{subequations}
\label{ADMM}
\begin{align}  
w^{t+1} & \leftarrow \arg \min_{w \in \mathbb{R}^{n}} \ \sum_{p=1}^P \{ \langle \lambda^t_p, w \rangle + \frac{\rho^t}{2} \|w-z^t_p\|^2 \}, \label{ADMM-1} \\
z^{t+1}_p & \leftarrow \arg \min_{z_p \in \mathcal{W}_p} \ f_p(z_p) - \langle \lambda^t_p, z_p \rangle + \frac{\rho^t}{2} \|w^{t+1}-z_p\|^2, \ \forall p \in [P], \label{ADMM-2} \\
\lambda^{t+1}_p &  \leftarrow  \lambda^{t}_p + \rho^t (w^{t+1}-z^{t+1}_p), \ \forall p \in [P], \label{ADMM-3}
\end{align}
\end{subequations}
where $\lambda_p$ is the Lagrange multiplier associated with the consensus constraints \eqref{dist_1_consensus} and $\rho^t > 0$ is a penalty parameter.


By noting that \eqref{ADMM-2} can be a computational bottleneck, solving \eqref{ADMM-2} inexactly has been proposed in the literature without affecting overall convergence. For example, linearized ADMM \cite{gao2018information, lin2017extragradient} has been proposed whose convergence rate is proven to be $\mathcal{O}(1/T)$ when the objective function to be linearized is smooth.
In this paper we approximate \eqref{ADMM-2} as follows:
\begin{align}
  z^{t+1}_p \gets \arg \min_{z_p  \in \mathcal{W}_p } \ \langle f'_p(z^t_p), z_p \rangle + \frac{1}{2\eta^{t}}\|z_p - z^t_p\|^2 +  \frac{\rho^t}{2} \|w^{t+1}-z_p + \frac{1}{\rho^t} \lambda^t_p  \|^2,  \label{ADMM-2-Prox} 
\end{align}
which is derived from \eqref{ADMM-2} by (i) replacing the convex function $f_p$ in \eqref{ADMM-2} with its lower approximation $\widehat{f}_p(z_p) := f_p(z^t_p) + \langle f'_p(z^t_p), \ z_p - z_p^t \rangle$, where $f'_p(z^t_p)$ is a subgradient of $f_p$ at $z^t_p$, and (ii) adding a proximal term $\frac{1}{2\eta^{t}}\|z_p - z^t_p\|^2$ with a parameter $\eta^{t} > 0$ that controls the proximity of a new solution $z^{t+1}_p$ from $z^t_p$ computed from the previous iteration. 
We note that when $\mathcal{W}_p := \mathbb{R}^n$ or $\mathcal{W}_p := \{ z_p \in \mathbb{R}^n: l \leq z_p \leq u \}$ for some $l$ and $u$, a closed-form solution expression can be derived from the optimality condition. 


\subsection{Differential privacy} \label{sec:prelim_dp}
In this section we describe DP, a privacy-preserving technique that provides statistical guarantee on data privacy.  
\begin{definition}[Differential privacy \cite{dwork2014algorithmic}] \label{def:differential_privacy_1}
A randomized function $\mathcal{A}$ satisfies $(\bar{\epsilon},\bar{\delta})$-DP if, for all neighboring datasets $\mathcal{D}$ and $\mathcal{D}'$ that differ in a single entry and for all possible events $\mathcal{S}$, the following inequality holds:
\begin{align}
 \mathbb{P}( \mathcal{A}(\mathcal{D}) \in \mathcal{S} ) \leq  e^{\bar{\epsilon}} \mathbb{P}( \mathcal{A}(\mathcal{D}') \in \mathcal{S} ) + \bar{\delta} ,  \label{def_differential_privacy}
\end{align}
where $\bar{\epsilon} \geq 0$, $\bar{\delta} \geq 0$, and $\mathcal{A}(\mathcal{D})$ (resp. $\mathcal{A}(\mathcal{D}')$) represent randomized outputs of $\mathcal{A}$ on input $\mathcal{D}$ (resp. $\mathcal{D}'$). 
Note that in this paper we refer to 
$\bar{\epsilon}$ as the  privacy budget,
$\ln ( \frac{ \mathbb{P}( \mathcal{A}(\mathcal{D}) \in \mathcal{S} ) }{ \mathbb{P}( \mathcal{A}(\mathcal{D}') \in \mathcal{S} )} )$ as the  privacy loss, and a function $\mathcal{A}$ satisfying DP (i.e., \eqref{def_differential_privacy}) as a  mechanism.
\end{definition}

If $\bar{\delta}=0$, we say that $\mathcal{A}$ is $(\bar{\epsilon},0)$-DP, namely, pure DP. 
One can easily see from \eqref{def_differential_privacy} that the absolute value of privacy loss, $|\ln ( \frac{ \mathbb{P}( \mathcal{A}(\mathcal{D}) \in \mathcal{S} ) }{ \mathbb{P}( \mathcal{A}(\mathcal{D}') \in \mathcal{S} )} )|$, is bounded by $\bar{\epsilon}$.
Therefore, as the privacy budget $\bar{\epsilon}$ decreases, it becomes harder to distinguish any neighboring datasets $\mathcal{D}$ and $\mathcal{D}'$ by analyzing the randomized outputs, thus providing stronger data privacy on any single entry in the dataset.
If $\bar{\delta}>0$, we say that $\mathcal{A}$ is $(\bar{\epsilon},\bar{\delta})$-DP, namely, an approximate DP that provides weaker privacy guarantee compared with the pure DP.
According to \cite[Lemma 3.17]{dwork2014algorithmic}, $(\bar{\epsilon},\bar{\delta})$-DP is equivalent to saying that the absolute value of privacy loss is bounded by $\bar{\epsilon}$ with probability $1-\bar{\delta}$. Note that $\bar{\delta} < 1$.

In what follows, we present the existing Laplace and Gaussian mechanisms.
\begin{lemma}[Laplace mechanism \cite{dwork2014algorithmic}]\label{lem:existing-laplace}
Given any function $\mathcal{R}$ that takes data $\mathcal{D}$ as input and outputs a vector in $\mathbb{R}^n$, the Laplace mechanism is defined as
\begin{subequations}
\label{LaplaceMechanism}
\begin{align}
\mathcal{A}(\mathcal{D}) := \mathcal{R}(\mathcal{D}) + \tilde{\xi}, \label{A}
\end{align}
where $\tilde{\xi} \in \mathbb{R}^n$ whose elements are independent and identically distributed (IID) random variables drawn from a Laplace distribution with zero mean and a variance $2 (\bar{\Delta}_1 / \bar{\epsilon})^2$. 
Here, $\bar{\epsilon}$ is a privacy budget from \eqref{def_differential_privacy}, and $\bar{\Delta}_1$ is an $L_1$-sensitivity defined as
\begin{align}
& \bar{\Delta}_1 := \max_{\forall (\mathcal{D},\mathcal{D}')} \| \mathcal{R}(\mathcal{D}) - \mathcal{R}(\mathcal{D}') \|_1,  \label{L1_sensitivity}
\end{align}
\end{subequations}
where $(\mathcal{D},\mathcal{D}')$ are two neighboring datasets.
The Laplace mechanism $\mathcal{A}$ in \eqref{LaplaceMechanism} satisfies $(\bar{\epsilon},0)$-DP.
\end{lemma}


\begin{lemma}[Gaussian mechanism \cite{dwork2014algorithmic}]\label{lem:existing-gaussian}
The Gaussian mechanism takes the form of \eqref{A} where $\tilde{\xi} \in \mathbb{R}^n$ is a noise vector whose elements are IID random variables drawn from a normal distribution with zero mean and a variance $2 \ln(1.25 / \bar{\delta}) (\bar{\Delta}_2 / \bar{\epsilon})^2$, where
\begin{align}
\bar{\Delta}_2:= \max_{\forall (\mathcal{D},\mathcal{D}')} \| \mathcal{R}(\mathcal{D}) - \mathcal{R}(\mathcal{D}') \|_2  \label{L2_sensitivity}
\end{align}
represents $L_2$-sensitivity. The Gaussian mechanism $\mathcal{A}$ satisfies $(\bar{\epsilon},\bar{\delta})$-DP, where $\bar{\delta} > 0$.
\end{lemma}
Even though the Gaussian mechanism that ensures  $(\bar{\epsilon},\bar{\delta})$-DP with $\bar{\delta}>0$ provides a weaker privacy guarantee compared with the Laplace mechanism that ensures $(\bar{\epsilon},0)$-DP, the Gaussian mechanism could lead to less noise since its variance relies on $L_2$ sensitivity, which is smaller than $L_1$ sensitivity used in the Laplace mechanism.
For this reason, the Gaussian mechanism has been widely adopted in numerous learning algorithms since the amount of noise introduced affects the learning performance.
Moreover, the Gaussian mechanism is favorable in FL algorithms where DP is applied for every iteration of the algorithms. 
The reason is that a mechanism that permits $T$ adaptive interactions with the Gaussian mechanisms ensures $( \bar{\epsilon} \sqrt{2T\ln(1/\bar{\delta}')} + 2k\epsilon^2, T \bar{\delta} + \bar{\delta}' ) $-DP, which is stronger than a mechanism that permits $T$ adaptive interactions with the Laplace mechanisms and ensures $(T\bar{\epsilon}, T\bar{\delta})$-DP, as shown in \cite{dwork2010boosting}.
Moreover, Abadi et al.~\cite{abadi2016deep} proposed a stronger composition theorem for the Gaussian mechanisms.

\begin{lemma}[Stronger composition theorem \cite{abadi2016deep}] \label{lem:composition}
An algorithm $\mathcal{A}$ that consists of a sequence of $T$ adaptive Gaussian mechanisms, namely, $\mathcal{A}^1, \ldots, \mathcal{A}^T$, each of which is proven to be $(\bar{\epsilon},\bar{\delta})$-DP with $\bar{\delta}>0$, is $(\epsilon', \bar{\delta})$-DP, where 
\begin{subequations}
\begin{align}
\ln(\bar{\delta}) = \ &  \min_{\tau \in \mathbb{Z}_+} \sum_{t=1}^T \alpha^t (\tau) - \tau \epsilon',  \label{lem:composition-1} \\
\alpha^{t} (\tau) := & \max_{\forall (\mathcal{D}, \mathcal{D}')} \ln \Bigg( \mathbb{E}_{\psi} \Bigg[ \Big( \frac{\textbf{pdf} \big( \mathcal{A}^{t} (\mathcal{D}) = \psi \big)}{\textbf{pdf} \big(  \mathcal{A}^{t} (\mathcal{D}') = \psi \big)} \Big)^{\tau} \Bigg] \Bigg), \ \forall t \in [T]. \label{lem:composition-2}
\end{align}
\end{subequations} 
Here $\alpha^{t} (\tau)$ represents all possible log of the moment-generating function at $\tau \in \mathbb{Z}_+$.
\end{lemma} 

In this paper we propose Laplace and Gaussian mechanisms in the context of the linearized ADMM with \texttt{ObjPert}, which are different from the Laplace and Gaussian mechanisms in Lemmas \ref{lem:existing-laplace} and \ref{lem:existing-gaussian}, respectively, with respect to the sensitivity computations as in \eqref{L1_sensitivity} and \eqref{L2_sensitivity}, respectively.
We also show that the  proposed Gaussian mechanism allows the stronger composition theorem to be applied.
  
 

 
 
\section{Differentially private linearized ADMM with multiple local updates} \label{sec:dpadmm}
In this section we propose a differentially private linearized ADMM algorithm for solving the distributed convex optimization model \eqref{model:dist}.
The proposed algorithm is equipped with (i) \textit{multiple local updates} that can reduce costs of communication between agents and (ii) \textit{objective perturbation} (i.e., \texttt{ObjPert}), which ensures data privacy in a DP manner without losing feasibility.
The proposed algorithm is especially effective when $\mathcal{W}_p \subset \mathbb{R}^n$, mainly because each agent communicates randomized (to ensure data privacy) but feasible local solutions that could affect the overall convergence of the algorithm.
In the rest of this section we describe the proposed algorithm (Section \ref{sec:proposed_algo}) and present a detailed comparison with the existing DP algorithm based on \texttt{OutPert} (Section \ref{sec:baseline}).
The privacy and convergence analyses are presented in Sections \ref{sec:privacy} and \ref{sec:convergence}, respectively.

\subsection{The proposed algorithm} \label{sec:proposed_algo}
The proposed algorithm is built on the linearized ADMM, that is, $\{ \eqref{ADMM-1} \rightarrow \eqref{ADMM-2-Prox} \rightarrow \eqref{ADMM-3} \}_{t=1}^T$.
Note that each agent solves \eqref{ADMM-2-Prox} using its local data and machines, which can be considered as a single local update.
In our algorithm, each agent conducts multiple local updates.
Specifically, for every $e \in [E]$, where $E$ is the number of local updates, each agent solves the following problem:
\begin{align}
z_p^{t,e+1} \leftarrow  \arg \min_{z_p  \in \mathcal{W}_p } \ &  \langle f'_p(z^{t,e}_p), z_p \rangle + \frac{1}{2\eta^{t}}\|z_p - z^{t,e}_p\|^2 + \frac{\rho^t}{2} \|w^{t+1}-z_p + \frac{1}{\rho^{t}} \lambda^t_p  \|^2.  \label{ADMM-2-MLU} 
\end{align}
The multiple local updates technique has been utilized in the federated averaging algorithm \cite{mcmahan2017communication} to reduce the communication costs. 
Inspired by that work, we introduce this technique in the linearized ADMM.

We note that the objective function of \eqref{ADMM-2-MLU} is strongly convex and $\mathcal{W}_p$ is a closed convex set; thus it has a unique optimal solution.
Therefore, the mapping from $f'_p(z_p^{t,e})$ to $z_p^{t,e+1}$ is injective (i.e., every $z_p^{t,e+1}$ has at most one $f'_p(z_p^{t,e})$). Therefore, the release of $z_p^{t,e+1}$ can lead to the leakage of $f'_p(z_p^{t,e})$, which contains the data information. 
It has been shown that the data can be reconstructed from the gradients (e.g., \cite{zhu2019deep, zhao2020idlg, jeon2021gradient}).

To protect data, one can simply add some calibrated noise $\tilde{\xi}^{t,e}_p$ to the output $z^{t,e+1}_p$, which is known as the output perturbation (i.e., \texttt{OutPert}) method. Note that the Laplace and Gaussian mechanisms described in Section \ref{sec:prelim_dp} belong to \texttt{OutPert} since they add Laplacian or Gaussian noise to the output directly to obtain randomized outputs to release. 
By noting that the resulting randomized outputs (that ensure data privacy) may be infeasible to $\mathcal{W}_p$, which may affect overall convergence of the algorithm, we alternatively add noise to the objective function of the constrained subproblem \eqref{ADMM-2-MLU} so that the resulting randomized outputs to release are always feasible to $\mathcal{W}_p$. 
Specifically, we add an affine function $\frac{1}{2\rho^t} \| \tilde{\xi}^{t,e}_p \|^2 - \langle w^{t+1} - z_p + \frac{1}{\rho^t} \lambda^t_p, \tilde{\xi}^{t,e}_p \rangle$ to \eqref{ADMM-2-MLU}, resulting in
\begin{align}
  z^{t,e+1}_p \leftarrow & \arg \min_{z_p \in \mathcal{W}_p} \ \langle f'_p(z^{t,e}_p), z_p \rangle + \frac{1}{2\eta^{t}} \| z_p - z_p^{t,e}\|^2 + \frac{\rho^t}{2} \|w^{t+1}-z_p + \frac{1}{\rho^t}(\lambda^t_p - \tilde{\xi}^{t,e}_p) \|^2,  \label{DPADMM-2-Prox}
\end{align}  
where $\tilde{\xi}^{t,e}_p \in \mathbb{R}^{n}$ is a noise vector that is guaranteed to preserve data privacy in a DP manner. In Section \ref{sec:privacy} we discuss how to generate such noise.
 

The proposed algorithm is presented in Algorithm \ref{algo:DP-IADMM-Prox}.
The computation at the central server is described in lines 1--9, while the local computation for each agent $p$ is described in lines 11--23.
In lines 2--3, the initial points are sent from the server to all agents.
In line 5, the central server computes a global solution $w^{t+1}$ via the equation derived from the optimality condition of \eqref{ADMM-1}.
In line 6, the server broadcasts $w^{t+1}$ to all local agents.
In lines 15--22, each agent receives $w^{t+1}$ from the server, conducts local updates for $E$ times, and sends the resulting local solution $z_p^{t+1}$ to the server.
Note that $z_p^{t+1}$ is an average of $E$ numbers of randomized outputs, each of which ensures data privacy.
In line 8 and line 23, the dual updates \eqref{ADMM-3} are performed at the server and the local agents individually. 
Note that these dual updates are identical since the initial points at the server and the local agents are the same.  

\begin{algorithm}[!ht]
\caption{Proposed DP linearized ADMM with multiple local updates}
\label{algo:DP-IADMM-Prox}
\begin{algorithmic}[1]
\STATE \textbf{(Server):}
\STATE Initialize $\lambda^1_1, \ldots, \lambda^1_P, z^1_1, \ldots, z^1_P$.
\STATE \texttt{Send} $\lambda^1_p, z^1_p$ to all agent $p \in [P]$ (\textbf{to line 12}).
\FOR{$t \in [T]$}
\STATE $w^{t+1} \leftarrow \frac{1}{P} \sum_{p=1}^P ( z_p^t - \frac{1}{\rho^t} \lambda^t_p) $.
\STATE \texttt{Send} $w^{t+1}$ to all agents (\textbf{to line 15}).
\STATE \texttt{Receive} $z^{t+1}_p$ from all agents (\textbf{from line 22}).
\STATE $\lambda^{t+1}_p \leftarrow \lambda^t_p + \rho^t(w^{t+1}-z^{t+1}_p)$ for all $p \in [P]$.
\ENDFOR
\STATE  
\STATE \textbf{(Agent $p \in [P]$):}  
\STATE \texttt{Receive} $\lambda^1_p, z^1_p$ from the server (\textbf{from line 3}).
\STATE Initialize $z^{0,E+1}_p = z^1_p$.
\FOR{$t \in [T]$}
\STATE \texttt{Receive} $w^{t+1}$ from the server (\textbf{from line 6}).    
\STATE Set $z_{p}^{t,1} \leftarrow z_p^{t-1,E+1}$
\FOR{$e \in [E]$}
\STATE Compute $z_p^{t,e+1}$ via \eqref{DPADMM-2-Prox}.
\ENDFOR        
\STATE $z^{t+1}_p \leftarrow \frac{1}{E} \sum_{e=1}^E z_p^{t,e+1}$.    
\STATE \texttt{Send} $z^{t+1}_p$ to the server (\textbf{to line 7}).
\STATE $\lambda^{t+1}_p \leftarrow \lambda^t_p + \rho^t(w^{t+1}-z^{t+1}_p)$.
\ENDFOR
\end{algorithmic}
\end{algorithm}
    
The benefits of Algorithm \ref{algo:DP-IADMM-Prox} are as follows.
First, 
solving \eqref{DPADMM-2-Prox} can be computationally cheaper than its exact version. Furthermore, for some special cases when $\mathcal{W}_p$ is a box constraint, \eqref{DPADMM-2-Prox} admits a closed-form solution expression.
Second, the quality of the solution can be improved via the multiple local updates that could result in saving communication rounds.
Third, the amount of communication is reduced by excluding the communication of the dual information.
This is different from the existing work \cite{zhou2021communication} that considers both multiple local primal and dual updates, namely, solving \eqref{ADMM-2-Prox} and \eqref{ADMM-3} multiple times, resulting in communicating not only local solution $z^{t+1}_p$ but also dual information $\lambda^{t+1}_p$.
Fourth and most important, data privacy is guaranteed in a DP manner for any communication rounds.


\subsection{Comparison with the existing DP algorithm} \label{sec:baseline}
We note that the authors in \cite{huang2019dp, huang2020differentially} proposed a differentially private linearized ADMM with multiple local updates when $\mathcal{W}_p=\mathbb{R}^n$.
Specifically, they employed a Gaussian mechanism (described in Section \ref{sec:prelim_dp}) that adds Gaussian noise $\tilde{\xi}^{t,e}_p$ to the output $z_p^{t,e+1}$ in \eqref{ADMM-2-MLU}.
We note that when $\mathcal{W}_p=\mathbb{R}^n$, adding Gaussian noise to the output (i.e., \texttt{OutPert}) is exactly the same as adding $- ( 1/\eta^t + \rho^t) \langle \tilde{\xi}^{t,e}_p, z_p \rangle$ to the objective function of \eqref{ADMM-2-MLU} (i.e., \texttt{ObjPert}).
We formalize the argument in the following remark.
\begin{remark}
Suppose that $\mathcal{W}_p=\mathbb{R}^n$ in \eqref{ADMM-2-MLU}.
The output perturbation is given by 
\begin{align}
z_p^{t,e+1} = \Big(\frac{1}{\eta^t} +\rho^t \Big)^{-1} \Big( - f'_p(z_p^{t,e}) + \frac{1}{\eta^t} z_p^{t,e} +\rho^t w^{t+1} + \lambda^t_p \Big) + \tilde{\xi}^{t,e}_p, \label{remark_opt_cond}
\end{align}
where $\tilde{\xi}^{t,e}_p$ is the Gaussian noise and the first term is derived by the optimality condition of \eqref{ADMM-2-MLU} with $\mathcal{W}_p=\mathbb{R}^n$.
Now consider the following problem constructed by the objective perturbation method:
\begin{align}
z_p^{t,e+1} \leftarrow  \arg \min_{z_p  \in \mathbb{R}^n } \  & \langle f'_p(z^{t,e}_p), z_p \rangle + \frac{1}{2\eta^{t}}\|z_p - z^{t,e}_p\|^2 + \frac{\rho^t}{2} \|w^{t+1}-z_p + \frac{1}{\rho^{t}} \lambda^t_p  \|^2 \nonumber \\
& - ( 1/\eta^t + \rho^t) \langle \tilde{\xi}^{t,e}_p, z_p \rangle. \nonumber
\end{align}
From its optimality condition, one can derive \eqref{remark_opt_cond}.
\end{remark}

We also note that when $\mathcal{W}_p = \mathbb{R}^n$, the multiple local update used in \cite{hager2020convergence} is exactly the same as solving \eqref{ADMM-2} by a subgradient descent algorithm up to $E$ iterations given by  
\begin{align}
z^{t, e+1}_p = z^{t,e}_p - s \Big( f_p'(z_p^{t,e})  - \rho^t (w^{t+1}-z_p^{t,e} + \frac{1}{\rho^t} \lambda^t )  \Big), \ \forall e \in [E], \nonumber
\end{align}
where $s = \eta^t/(\eta^t \rho^t +1)$ is a constant step size. 
To see this, one can derive the optimality condition of \eqref{ADMM-2-MLU} when $\mathcal{W}_p = \mathbb{R}^n$:
\begin{align*}
0 & = f'_p(z_p^{t,e}) + \frac{1}{\eta^t} (z_p^{t,e+1} - z_p^{t,e}) - \rho^t(w^{t+1} - z_p^{t,e+1}  + \frac{1}{\rho^t} \lambda^t_p ) \\  
& = f'_p(z_p^{t,e}) + \frac{1}{\eta^t} (z_p^{t,e+1} - z_p^{t,e}) - \rho^t(w^{t+1} - z_p^{t,e}  + \frac{1}{\rho^t} \lambda^t_p ) - \rho^t (z_p^{t,e} - z_p^{t,e+1}) \\
& = \eta^t \Big( f'_p(z_p^{t,e}) - \rho^t(w^{t+1} - z_p^{t,e}  + \frac{1}{\rho^t} \lambda^t_p )  \Big) - (\eta^t \rho^t+1) (z_p^{t,e} - z_p^{t,e+1}) \\
& = -\frac{\eta^t}{\eta^t \rho^t+1} \Big( f'_p(z_p^{t,e})  + \rho^t(w^{t+1} - z_p^{t,e}  + \frac{1}{\rho^t} \lambda^t_p )  \Big)  +  (z_p^{t,e} - z_p^{t,e+1}).
\end{align*}
To summarize, when $\mathcal{W}_p = \mathbb{R}^n$, there is no significant difference between \texttt{OutPert} and \texttt{ObjPert}.
Also, the multiple local updates simply involve solving the exact problem \eqref{ADMM-2} by the subgradient descent algorithm for a finite number of iterations.

Different from the literature,  this paper focuses on a situation when $\mathcal{W}_p \subset \mathbb{R}^n$, which appears in many application areas. 
In this case, we aim to show that \texttt{ObjPert} outperforms \texttt{OutPert}.
To do so, we construct a baseline algorithm based on \texttt{OutPert}, which is a form of Algorithm \ref{algo:DP-IADMM-Prox} where \eqref{DPADMM-2-Prox} in line 18 is replaced with
\begin{align}
z_p^{t,e+1} \gets z_p^{t,e+1} \text{ from \eqref{ADMM-2-MLU} } + \tilde{\xi}^{t,e}_p, \label{baseline}
\end{align}
where $\tilde{\xi}^{t,e}_p$ should be sampled to ensure DP.
However, the randomized output $z_p^{t+1} = \frac{1}{E} \sum_{e=1}^E z_p^{t,e+1}$, where $z_p^{t,e+1}$ is computed by \eqref{baseline}, may not be feasible to $\mathcal{W}_p$, whereas the randomized output $z_p^{t+1}$ in line 20 of Algorithm \ref{algo:DP-IADMM-Prox} is always feasible (see Figure \ref{fig:feasible}).

\begin{figure}[!ht]  
\centering
\includegraphics[width=0.45\textwidth]{./Figures/obj_perturb.pdf}  
\includegraphics[width=0.45\textwidth]{./Figures/out_perturb.pdf}  
\caption{Comparison of \eqref{DPADMM-2-Prox} and \eqref{baseline} for fixed $t$, $p$, and $E=3$. 
Every randomized outputs produced by \eqref{DPADMM-2-Prox} is feasible, whereas those produced by \eqref{baseline} may not be.}
\label{fig:feasible}
\end{figure}
While the privacy and convergence analyses of the baseline algorithm can be immediately borrowed from \cite{huang2020differentially}, it is not the case for the proposed algorithm (i.e., Algorithm \ref{algo:DP-IADMM-Prox}) because \texttt{ObjPert} is different from \texttt{OutPert} when $\mathcal{W}_p \subset \mathbb{R}^n$. 
Notably, noise generated for ensuring DP in Algorithm \ref{algo:DP-IADMM-Prox} is based on the sensitivity of the subgradients, which is different from the sensitivity of the true outputs (i.e., \eqref{L1_sensitivity} or \eqref{L2_sensitivity}) in the literature. 
Also, \texttt{ObjPert} under the existence of $\mathcal{W}_p \subset \mathbb{R}^n$ requires  construction of a set of mechanisms, each of which is $(\bar{\epsilon}, \bar{\delta})$-DP, that provides a sequence of the randomized outputs that converges to a feasible randomized output, which is not needed in the privacy analysis for the case when $\mathcal{W}_p = \mathbb{R}^n$ as in \cite{huang2020differentially}. 
More details on our privacy analysis will be discussed in Section \ref{sec:privacy}.
Furthermore, we show that the rate of convergence of Algorithm \ref{algo:DP-IADMM-Prox} is sublinear without any error bound, which is different from the convergence analysis in \cite{huang2020differentially}, where the rate of convergence is $\mathcal{O}(1/\sqrt{T}) + \mathcal{O}(1/\bar{\epsilon}^2)$, where the error bound $\mathcal{O}(1/\bar{\epsilon}^2)$ increases as $\bar{\epsilon}$ decreases for stronger data privacy. 
More details on our convergence analysis  will be discussed in Section \ref{sec:convergence}.







 

\section{Privacy analysis} \label{sec:privacy}
In this section we provide a privacy analysis of the proposed algorithm.
Specifically, we will use the following lemma to show that the constrained subproblem \eqref{DPADMM-2-Prox} for fixed $t$, $e$, and $p$ provides a randomized output that guarantees data privacy in a DP manner. 
 
\begin{lemma}[Theorem 1 of \cite{kifer2012private}] \label{lemma:kifer}
Let $\mathcal{A}$ be a randomized algorithm induced by the random noise $\tilde{\xi}$ that provides output $\phi(\mathcal{D},\tilde{\xi})$.
Let $\{\mathcal{A}_{\ell}\}_{\ell=1}^{\infty}$ be a sequence of randomized algorithms, each of which is induced by $\tilde{\xi}$ and provides output $\phi^{\ell}(\mathcal{D},\tilde{\xi})$.
If $\mathcal{A}_{\ell}$ is $(\bar{\epsilon},\bar{\delta})$-DP for all $\ell$ and satisfies a pointwise convergence condition, namely, $\lim_{{\ell} \rightarrow \infty} \phi^{\ell}(\mathcal{D},\tilde{\xi}) = \phi(\mathcal{D},\tilde{\xi})$, then $\mathcal{A}$ is also $(\bar{\epsilon},\bar{\delta})$-DP, where $\bar{\delta} \geq 0$.
\end{lemma} 

In what follows, we construct a sequence $\{\mathcal{A}_{\ell}\}_{\ell=1}^{\infty}$ of randomized algorithms that satisfy the pointwise convergence condition in Section \ref{sec:randomized_algo}, and we propose two different mechanisms $\mathcal{A}_{\ell}$ for any $\ell$ as in Section \ref{sec:mechanisms}.



\subsection{A sequence of randomized algorithms} \label{sec:randomized_algo}
For the rest of this section we fix $t \in [T]$, $e \in [E]$, and $p \in [P]$. 
Also, we denote the objective function of \eqref{DPADMM-2-Prox}, which is strongly convex, by  
\begin{align}
  & G^{t,e}_p(z_p; \tilde{\xi}^{t,e}_p) :=  \langle f'_p(z^{t,e}_p), z_p \rangle   +   \frac{1}{2\eta^{t}} \| z_p - z_p^{t,e}\|^2 + \frac{\rho^t}{2} \|w^{t+1}-z_p + \frac{1}{\rho^t}(\lambda^t_p - \tilde{\xi}^{t,e}_p) \|^2 \label{fn_G} 
\end{align}
and the feasible region of \eqref{DPADMM-2-Prox} by
\begin{align*}
     \mathcal{W}_p := \{ z_p \in \mathbb{R}^{n} : h_m (z_p) \leq 0, \ \forall m \in [M_p] \}, 
\end{align*}
where $h_m$ is convex and twice continuously differentiable and $M_p$ is the total number of inequalities.

By utilizing an indicator function $\mathcal{I}_{\mathcal{W}_p}(z_p)$ that outputs zero if $z_p \in \mathcal{W}_p$ and $\infty$ otherwise, \eqref{DPADMM-2-Prox} can be expressed by the following problem:
\begin{align*}
\min_{z_p \in \mathbb{R}^{n}} \ G^{t,e}_p (z_p; \tilde{\xi}^{t,e}_p) + \mathcal{I}_{\mathcal{W}_p}(z_p).
\end{align*}
We observe that the indicator function can be approximated by the following penalty function:
\begin{align}
& g_{\ell} (z_p) :=  \sum_{m=1}^{M_p} \ln ( 1+ e^{\ell h_m(z_p)}), \label{logfn-prox}  
\end{align}
where $\ell > 0$. 
Increasing $\ell$ enforces the feasibility, namely, $h_m(z_p) \leq 0$, resulting in $g_{\ell}(z_p) \rightarrow 0$.
It is similar to the logarithmic barrier function (LBF) $-(1/\ell) \sum_{m=1}^M \ln (-h_m(z_p))$, in that the approximation becomes closer to the indicator function as $\ell \rightarrow \infty$.
However, the penalty function $g_{\ell}$ is different from LBF in that $z_p$ is not restricted.
See Figure \ref{fig:penalty} for an example of the penalty function compared with LBF.
\begin{figure}[!ht]  
\centering
\includegraphics[width=0.45\textwidth]{./Figures/penalty.jpg}  
\includegraphics[width=0.45\textwidth]{./Figures/barrier.jpg}  
\caption{Example of the penalty function in \eqref{logfn-prox} with $M_p=1$ and $h_1(z_p) = z_p$ compared with the logarithmic barrier function.}
\label{fig:penalty}
\end{figure}
By replacing the indicator function with the penalty function in \eqref{logfn-prox}, we construct the following \textit{unconstrained} problem (which is $\mathcal{A}_{\ell}$ in the context of Lemma \ref{lemma:kifer}):
\begin{align}
z^{t,e+1}_{p\ell} \leftarrow  \arg \min_{z_p \in \mathbb{R}^{n}} \ & G^{t,e}_p (z_p; \tilde{\xi}^{t,e}_p)  + g_{\ell}(z_p), \label{ADMM-2-Prox-log}
\end{align}
where the objective function is strongly convex because $g_{\ell}$ is convex over all domains and $G^{t,e}_p$ is strongly convex. Therefore, $z^{t,e+1}_{p\ell}$ is the unique optimal solution.

 
\begin{theorem}\label{thm:pointwise_convergence}
  For fixed $t$, $e$, and $p$,  
  \begin{align}
  \lim_{\ell \rightarrow \infty} z^{t,e+1}_{p \ell}  = z_p^{t,e+1},  
  \end{align}
  where $z_p^{t,e+1}$ and $z^{t,e+1}_{p \ell}$ are from \eqref{DPADMM-2-Prox} and \eqref{ADMM-2-Prox-log}, respectively.
\end{theorem}
\begin{proof}
 
  We develop the proof by contradiction.
  To this end, we suppose that $z^{t,e+1}_{p\ell}$ converges to $\widehat{z} \neq z^{t,e+1}_p$ as $\ell$ increases.
  Consider $\zeta := \| \widehat{z} - z^{t,e+1}_p \| / 2$.
  Since $z^{t,e+1}_{p\ell}$  converges to $\widehat{z}$, there exists $\ell' > 0$ such that $\| \widehat{z} - z^{t,e+1}_{p\ell} \| < \zeta$ for all $\ell \geq \ell'$.
  By the triangle inequality, we have
  \begin{subequations}
  \begin{align}
  \| z^{t,e+1}_{p\ell} - z^{t,e+1}_p \| \geq \| \widehat{z} - z^{t,e+1}_p \| - \| \widehat{z} - z^{t,e+1}_{p\ell} \| > 2\zeta - \zeta = \zeta, \ \forall \ell \geq \ell'. \label{triangle}
  \end{align}
  Since $G^{t,e}_p$ is strongly convex with a constant $\mu > 0$, we have
  \begin{align}
  G^{t,e}_p(z^{t,e+1}_{p\ell}) - G^{t,e}_p(z^{t,e+1}_p) \geq \textstyle \frac{\mu}{2} \| z^{t,e+1}_{p\ell} - z^{t,e+1}_p \|^2 >  \frac{\mu \zeta^2}{2}, \ \forall \ell \geq \ell', \label{strong_triangle}
  \end{align}
  where the last inequality holds by \eqref{triangle}.
  By adding $g_{\ell} (z^{t,e+1}_{p\ell}) \geq 0$ to the left-hand side of \eqref{strong_triangle}, we derive the following inequality:
  \begin{align}
  \big\{ G^{t,e}_p(z^{t,e+1}_{p\ell}) + g_{\ell} (z^{t,e+1}_{p\ell})  \big\} - G^{t,e}_p(z^{t,e+1}_p)  > \textstyle \frac{\mu \zeta^2}{2}, \ \forall \ell \geq \ell'. \label{diff_0}
  \end{align}  
  The continuity of $G^{t,e}_p: \mathcal{W}_p \mapsto \mathbb{R}$ at $z^{t,e+1}_p$ implies that for every $\epsilon > 0$, there exists a $\delta > 0$ such that for all $z \in \mathcal{W}_p$:
  \begin{align}
  z \in \mathcal{B}_{\delta}(z^{t,e+1}_p) := \{ z \in \mathbb{R}^{n} : \| z - z^{t,e+1}_p \| < \delta \} \ \Rightarrow \  G^{t,e}_p(z) - G^{t,e}_p(z^{t,e+1}_p)  < \epsilon. \label{continuity}
  \end{align}
  Consider $\tilde{z} \in \mathcal{B}_{\delta}(z^{t,e+1}_p) \cap \textbf{relint}(\mathcal{W}_p)$, where \textbf{relint} indicates the relative interior.
  Since $h_m(\tilde{z}) < 0$ for all $m \in [M_p]$, $g_{\ell}(\tilde{z})$ goes to zero as $\ell$ increases.
  Hence, there exists $\ell'' > 0$ such that
  \begin{align}
    g_{\ell}(\tilde{z}) =  \sum_{m=1}^{M_p} \ln (1 + e^{\ell h_m(\tilde{z})}) < \epsilon, \ \forall \ell \geq \ell''. \label{ineq_z_tilde}
  \end{align}
  For all $\ell \geq \ell''$, we derive the following inequalities:
  \begin{align}
    G^{t,e}_p(z^{t,e+1}_{p\ell}) + g_{\ell}(z^{t,e+1}_{p\ell}) \leq G^{t,e}_p(\tilde{z}) + g_{\ell}(\tilde{z}) < G^{t,e}_p(\tilde{z}) + \epsilon < G^{t,e}_p(z^{t,e+1}_p) + 2\epsilon, \label{sandwich_ineq}
  \end{align}
  where
  the first inequality holds because $z^{t,e+1}_{p\ell}$ is the optimal solution of \eqref{ADMM-2-Prox-log},
  the second inequality holds by \eqref{ineq_z_tilde}, and
  the last inequality holds by \eqref{continuity}.
  To see the contradiction, consider $\epsilon \in (0, \frac{\mu \zeta^2}{4})$.
  Then we have
  \begin{align}
  \big\{ G^{t,e}_p(z^{t,e+1}_{p\ell}) + g_{\ell} (z^{t,e+1}_{p\ell})  \big\} - G^{t,e}_p(z^{t,e+1}_p) \underbrace{<}_{\text{from } \eqref{sandwich_ineq}} 2 \epsilon < \frac{\mu \zeta^2}{2}, \ \forall \ell \geq \ell''. \label{diff_1}
  \end{align}  
  Therefore, for all $\ell \geq \max \{ \ell', \ell'' \}$, \eqref{diff_0} and \eqref{diff_1} contradict.
\end{subequations}
\end{proof}


  
\subsection{Two mechanisms} \label{sec:mechanisms}
In this section we propose two mechanisms, namely, variants of the existing Laplace and Gaussian mechanisms in Lemmas \ref{lem:existing-laplace} and \ref{lem:existing-gaussian}, respectively, based on the proposed \eqref{ADMM-2-Prox-log} for fixed $t$, $e$, $p$, and $\ell$.

First, we propose to sample $\tilde{\xi}^{t,e}_p$ in \eqref{ADMM-2-Prox-log} from a Laplace distribution with zero mean and a variance $2 (\bar{\Delta}^{t,e}_{p,1} / \bar{\epsilon})^2$, where $\bar{\Delta}^{t,e}_{p,1}$ is defined as follows:
\begin{align}
  \bar{\Delta}^{t,e}_{p,1} := \max_{ \forall (\mathcal{D}_p, \mathcal{D}_p')} \| f'_p(z^{t,e}_p;\mathcal{D}_p) - f'_p(z^{t,e}_p;\mathcal{D}'_p)\|_1, \label{redefined_L1sensitivity}
\end{align}
which is different from \eqref{L1_sensitivity} in that the proposed $L_1$ sensitivity \eqref{redefined_L1sensitivity} is based on the difference of subgradients (i.e., inputs of \eqref{ADMM-2-Prox-log}) while the existing $L_1$ sensitivity \eqref{L1_sensitivity} is based on the difference of true outputs.
In the following theorem, we show that the proposed Laplace mechanism, namely, \eqref{ADMM-2-Prox-log} with the Laplacian noise, can achieve DP.
 
\begin{theorem} \label{thm:laplace}  
  For fixed $t$, $e$, $p$, and $\ell$, the proposed \eqref{ADMM-2-Prox-log} with the aforementioned Laplacian noise is $(\bar{\epsilon},0)$-DP: 
  \begin{align}
    \label{DP_1}
    \mathbb{P} ( z^{t,e+1}_{p \ell} (\mathcal{D}_p )\in \mathcal{S} ) \leq  e^{\bar{\epsilon}} \mathbb{P} ( z^{t,e+1}_{p \ell} (\mathcal{D}'_p )\in \mathcal{S} ), 
  \end{align}  
  for all possible events $\mathcal{S} \subset \mathbb{R}^{n}$ and neighboring datasets $(\mathcal{D}_p, \mathcal{D}'_p)$.
\end{theorem}
\begin{proof}
  It suffices to show that the following is true:
  \begin{subequations}
  \begin{align}
   
   
   
    \textbf{pdf} \big( z^{t,e+1}_{p\ell} (\mathcal{D}_p) = \psi \big)
    \leq e^{\bar{\epsilon}} \ \textbf{pdf} \big( z^{t,e+1}_{p\ell} (\mathcal{D}'_p) = \psi \big), \ \forall \psi \in \mathbb{R}^{n}, \label{DP_PDF}    
  \end{align} 
  where \textbf{pdf} stands for probability density function.
  This is because integrating $\psi$ in \eqref{DP_PDF} over $\mathcal{S}$ yields \eqref{DP_1}.

  Consider $\psi \in \mathbb{R}^{n}$.
  If we have $z^{t,e+1}_{p \ell} (\mathcal{D}_p)=\psi$, then $\psi$ is the unique minimizer of \eqref{ADMM-2-Prox-log} because the objective function in \eqref{ADMM-2-Prox-log} is strongly convex.
  From the optimality condition of \eqref{ADMM-2-Prox-log}, we derive
  \begin{align}
    \tilde{\xi}^{t,e}_p (\psi; \mathcal{D}_p) = & - f'_p(z^{t,e}_p;\mathcal{D}_p) + \rho^t(w^{t+1}-\psi) + \lambda^t_p - \nabla g_{\ell}(\psi)  -  \frac{1}{\eta^{t}} \big( \psi - z_p^{t,e} \big), \label{correspondence}
  \end{align}
  where $\nabla g_{\ell}(\psi) = \sum_{m=1}^{M_p} \frac{\ell e ^{\ell h_m(\psi)}}{1+e^{\ell h_m(\psi)}} \nabla h_m(\psi)$, and we use $\tilde{\xi}^{t,e}_p (\psi; \mathcal{D}_p)$ and $\tilde{\xi}^{t,e}_p$ interchangeably.
  Note that the mapping from $\psi$ to $\tilde{\xi}^{t,e}_p$ via \eqref{correspondence} is injective. 
  Also the mapping is surjective because for all $\tilde{\xi}^{t,e}_p$, there exists $\psi$ (i.e., the unique minimizer of \eqref{ADMM-2-Prox-log}) such that \eqref{correspondence} holds.
  Therefore, the relation between $\psi$ and $\tilde{\xi}^{t,e}_p$ is bijective, which allows us to utilize the inverse function theorem \cite[Theorem 17.2]{billingsley}, namely,
    \begin{align}
      \textbf{pdf} \big( z^{t,e+1}_{p\ell} (\mathcal{D}_p) = \psi \big) \cdot \big|\textbf{det}[\nabla \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}_p) ] \big| = \text{Lap} \big( \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}_p); 0, \bar{\Delta}_{p,1}^{t,e} /\bar{\epsilon} \big), \label{inverse}
    \end{align}
    where 
    $\text{Lap}(\cdot; 0,\bar{\Delta}^{t,e}_{p,1} / \bar{\epsilon} ) $ represents a probability density function of Laplace distribution with zero mean and a variance $2(\bar{\Delta}^{t,e}_{p,1} / \bar{\epsilon})^2$, 
    \textbf{det} represents a determinant of a matrix, and
    $\nabla \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}_p)$ represents a Jacobian matrix of the mapping from $\psi$ to $\tilde{\xi}^{t,e}_p$ in \eqref{correspondence}, namely,
    \begin{align}
    \nabla \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}_p) =  (-\rho^t - 1/\eta^{t}) \mathbb{I}_{n} - \nabla \Big( \sum_{m=1}^{M_p} \frac{\ell e ^{\ell h_m(\psi)}}{1+e^{\ell h_m(\psi)}} \nabla h_m(\psi) \Big) , \label{Jacobian}
    \end{align}
    where $\mathbb{I}_{n}$ is an identity matrix of $n \times n$ dimensions.
    Since the Jacobian matrix is not affected by the dataset, we have
    \begin{align}
    \nabla \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}_p) = \nabla \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}'_p).  \label{jacobian}
    \end{align}
    Based on \eqref{inverse} and \eqref{jacobian}, we derive the following inequalities:
    \begin{align}
      & \ln \Bigg( \frac{\textbf{pdf} \big( z^{t,e+1}_{p \ell} (\mathcal{D}_p) = \psi \big)}{ \textbf{pdf}  \big( z^{t,e+1}_{p \ell} (\mathcal{D}'_p ) = \psi \big)}  \Bigg) = \ln \Bigg( \frac{\text{Lap} \big( \tilde{\xi}^{t,e}_p (\psi; \mathcal{D}_p); 0, \bar{\Delta}_{p,1}^{t,e}/\bar{\epsilon} \big) }{\text{Lap} \big( \tilde{\xi}^{t,e}_p (\psi; \mathcal{D}'_p); 0, \bar{\Delta}_{p,1}^{t,e} /\bar{\epsilon} \big) } \Bigg) \label{laplacian-1} \\
      = \ &    (\bar{\epsilon}/\bar{\Delta}_{p,1}^{t,e})(\| \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}'_p) \|_1 - \| \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}_p) \|_1)  \nonumber \\
      \leq \ &   (\bar{\epsilon}/\bar{\Delta}_{p,1}^{t,e})(\| \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}'_p) - \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}_p) \|_1)  \label{triangle_inequality} \\
      = \ &   (\bar{\epsilon}/\bar{\Delta}_{p,1}^{t,e})(\| f'_p(z^{t,e}_p;\mathcal{D}_p) - f'_p(z^{t,e}_p;\mathcal{D}'_p)\|_1) \leq \bar{\epsilon} , \label{4.8b}  
    \end{align}    
    \end{subequations}
    where
    the equality in \eqref{laplacian-1} holds because of \eqref{inverse} and \eqref{jacobian},
    the inequality \eqref{triangle_inequality} holds because of the triangle inequality,
    the equality in \eqref{4.8b} holds because of \eqref{correspondence}, and
    the inequality in \eqref{4.8b} holds by the definition of $\bar{\Delta}_{p,1}^{t,e}$ in \eqref{redefined_L1sensitivity}.    
\end{proof}

Next, we propose to sample $\tilde{\xi}^{t,e}_p$ in \eqref{ADMM-2-Prox-log} from a Gaussian distribution with zero mean and a variance $(\sigma^{t,e}_p)^2 := 2 \ln(1.25/\bar{\delta}) (\bar{\Delta}^{t,e}_{p,2} / \bar{\epsilon})^2$, where $\bar{\Delta}^{t,e}_{p,2}$ is defined as follows:
\begin{align}
  \bar{\Delta}^{t,e}_{p,2} := \max_{ \forall (\mathcal{D}_p, \mathcal{D}_p')} \| f'_p(z^{t,e}_p;\mathcal{D}_p) - f'_p(z^{t,e}_p;\mathcal{D}'_p)\|_2, \label{redefined_L2sensitivity}
\end{align}
which is different from \eqref{L2_sensitivity} in that the proposed $L_2$ sensitivity \eqref{redefined_L2sensitivity} is based on the difference of subgradients wheres the existing $L_2$ sensitivity is based on the difference of true outputs.
In the following theorem, we show that the proposed Gaussian mechanism, namely, \eqref{ADMM-2-Prox-log} with the Gaussian noise, can achieve DP.

\begin{theorem} \label{thm:gaussian}  
  For fixed $t$, $e$, $p$, and $\ell$, the proposed \eqref{ADMM-2-Prox-log} with the aforementioned Gaussian noise is $(\bar{\epsilon},\bar{\delta})$-DP: 
  \begin{align}
    \label{DP_2}
    \mathbb{P} ( z^{t,e+1}_{p \ell} (\mathcal{D}_p )\in \mathcal{S} ) \leq  e^{\bar{\epsilon}} \mathbb{P} ( z^{t,e+1}_{p \ell} (\mathcal{D}'_p )\in \mathcal{S} ) + \bar{\delta}, 
  \end{align}  
  for all possible events $\mathcal{S} \subset \mathbb{R}^{n}$ and neighboring datasets $(\mathcal{D}_p, \mathcal{D}'_p)$.
\end{theorem}

\begin{proof}  
  We develop the proof based on the proof in Theorem \ref{thm:laplace}. 
  Specifically, \eqref{laplacian-1} can be rewritten as 
  \begin{subequations}
  \begin{align}
    & \ln \Bigg( \frac{\textbf{pdf} \big( z^{t,e+1}_{p \ell} (\mathcal{D}_p) = \psi \big)}{\textbf{pdf} \big( z^{t,e+1}_{p \ell} (\mathcal{D}'_p ) = \psi \big)} \Bigg) = \ln \Bigg( \frac{\mathcal{N} \big( \tilde{\xi}^{t,e}_p (\psi; \mathcal{D}_p); 0, (\sigma^{t,e}_p)^2 \big) }{ \mathcal{N} \big( \tilde{\xi}^{t,e}_p (\psi; \mathcal{D}'_p); 0, (\sigma^{t,e}_p)^2 \big) } \Bigg) \label{gaussian-1}  \\
    = & \frac{1}{2 (\sigma^{t,e}_p)^2} \Big(\| \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}'_p) \|^2 - \| \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}_p) \|^2 \Big) \nonumber   \\
    = &   \frac{1}{2 (\sigma^{t,e}_p)^2} \Big(\| \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}_p) + f'_p(z_p^{t,e};\mathcal{D}_p) - f'_p(z_p^{t,e};\mathcal{D}'_p)  \|^2 - \| \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}_p) \|^2 \Big)  \label{gaussian-2} \\ 
    = &   \frac{1}{2 (\sigma^{t,e}_p)^2} \Big( 2 \langle \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}_p), f'_p(z_p^{t,e};\mathcal{D}_p) - f'_p(z_p^{t,e};\mathcal{D}_p') \rangle \nonumber \\
    & + \|f'_p(z_p^{t,e};\mathcal{D}_p) - f'_p(z_p^{t,e};\mathcal{D}_p') \|^2 \Big), \nonumber 
  \end{align}
  where $\mathcal{N} \big( \cdot ; 0, (\sigma^{t,e}_p)^2 \big)$ in \eqref{gaussian-1} represents a probability density function of Gaussian distribution with zero mean and a variance $(\sigma^{t,e}_p)^2 := 2 \ln(1.25/\bar{\delta}) (\bar{\Delta}^{t,e}_{p,2} / \bar{\epsilon})^2$, where $\bar{\Delta}^{t,e}_{p,2}$ is from \eqref{redefined_L2sensitivity}, and 
  the equality in \eqref{gaussian-2} holds because of \eqref{correspondence}. 

  Let $x := \langle \tilde{\xi}^{t,e}_p (\psi;\mathcal{D}_p), f'_p(z_p^{t,e};\mathcal{D}_p) - f'_p(z_p^{t,e};\mathcal{D}_p') \rangle$. 
  Then $x \sim \mathcal{N}(0, (\sigma^{t,e}_p)^2 \|f'_p(z_p^{t,e};\mathcal{D}_p) - f'_p(z_p^{t,e};\mathcal{D}_p') \|^2 )$ since each element of $\tilde{\xi}^{t,e}_p$ is drawn from $\mathcal{N}(0, (\sigma^{t,e}_p)^2)$.    
  Let $y := x/(\sigma^{t,e}_p)^2 + \|f'_p(z_p^{t,e};\mathcal{D}_p) - f'_p(z_p^{t,e};\mathcal{D}_p') \|^2 / (2 (\sigma^{t,e}_p)^2)$ be the privacy loss random variable defined in \eqref{gaussian-1}. Then the mean and variance of $y$ are given by
  \begin{align}
    & \frac{\|f'_p(z_p^{t,e};\mathcal{D}_p) - f'_p(z_p^{t,e};\mathcal{D}_p') \|^2}{2 (\sigma^{t,e}_p)^2} \text { and } \frac{\|f'_p(z_p^{t,e};\mathcal{D}_p) - f'_p(z_p^{t,e};\mathcal{D}_p') \|^2}{(\sigma^{t,e}_p)^2}, \label{privacy_loss_rv_mean_var}
  \end{align}
  respectively.  
  Let $z \sim \mathcal{N}(0,1)$. Then the privacy loss random variable can be rewritten as 
  \begin{align}
    y := \frac{\|f'_p(z_p^{t,e};\mathcal{D}_p) - f'_p(z_p^{t,e};\mathcal{D}_p') \|}{\sigma^{t,e}_p} z + \frac{\|f'_p(z_p^{t,e};\mathcal{D}_p) - f'_p(z_p^{t,e};\mathcal{D}_p') \|^2}{2 (\sigma^{t,e}_p)^2}.
  \end{align}
  It suffices to prove the following:
  \begin{align}
    \mathbb{P}\big( |y| > \bar{\epsilon} \big) \leq \bar{\delta}, \label{eps_delta_dp}
  \end{align}
  which is equivalent to the definition of $(\bar{\epsilon},\bar{\delta})$-DP according to \cite[Lemma 3.17]{dwork2014algorithmic}.
  To this end, we derive an upper bound on the left-hand side of \eqref{eps_delta_dp} as 
  \begin{align}
    & \mathbb{P}\big( |y| > \bar{\epsilon} \big)  \nonumber \\
    = \ & \mathbb{P}\Big( |z| >\frac{\bar{\epsilon} \sigma^{t,e}_p}{\|f'_p(z_p^{t,e};\mathcal{D}_p) - f'_p(z_p^{t,e};\mathcal{D}_p') \|}   - \frac{\|f'_p(z_p^{t,e};\mathcal{D}_p) - f'_p(z_p^{t,e};\mathcal{D}_p') \|}{2 \sigma^{t,e}_p}\Big) \nonumber  \\
    \leq \ & \mathbb{P}\Big( |z| >  \sqrt{2 \ln(1.25/\bar{\delta})}   - \frac{\bar{\epsilon} }{2 \sqrt{2 \ln(1.25/\bar{\delta})}  }\Big)  \label{def_delta} \\
    \leq \ & \frac{2}{\sqrt{2 \pi}} \exp \Big(\frac{-1}{2} \Big(\sqrt{2 \ln(1.25/\bar{\delta})}   - \frac{\bar{\epsilon} }{2 \sqrt{2 \ln(1.25/\bar{\delta})}  } \Big)^2 \Big)  \label{tail_bound}  \\
    \leq \ & \bar{\delta},
  \end{align}
  where the inequality \eqref{def_delta} holds because of the definition of $\sigma^{t,e}_p$;
  the inequality \eqref{tail_bound} holds because of the tail bound of $\mathcal{N}(0,1)$ is given by $\mathbb{P}(|z| > r)  \leq \frac{2}{\sqrt{2\pi}} \exp(\frac{-r^2}{2})$; and
  the inequality holds because if $t=\sqrt{2 \ln(1.25/\bar{\delta})}$, then
  $(t - \bar{\epsilon}/t)^2 \geq 2 \ln (\frac{2}{\bar{\delta}} \frac{1}{\sqrt{2\pi}})$ for $\bar{\epsilon} \leq 1$ (see \cite[Appendix A]{dwork2014algorithmic} for more details). 
  \end{subequations}
\end{proof}


\subsection{Total privacy leakage}

The results from Sections \ref{sec:randomized_algo} and \ref{sec:mechanisms} allow us to use Lemma \ref{lemma:kifer} to prove that \eqref{DPADMM-2-Prox} preserves data privacy in a DP manner.

\begin{theorem} \label{thm:privacy}
  For fixed $t \in [T]$ and $e \in [E]$, each agent $p \in [P]$ can achieve (i) $(\bar{\epsilon}, 0)$-DP by solving \eqref{DPADMM-2-Prox}, where $\tilde{\xi}^{t,e}_p$ is drawn from a Laplace distribution with zero mean and a variance $2(\bar{\Delta}^{t,e}_{p,1}/\bar{\epsilon})^2$ and $\bar{\Delta}^{t,e}_{p,1}$ is from \eqref{redefined_L1sensitivity} and (ii) $(\bar{\epsilon}, \bar{\delta})$-DP with $\bar{\delta}>0$ by solving \eqref{DPADMM-2-Prox}, where $\tilde{\xi}^{t,e}_p$ is drawn from a Gaussian distribution with zero mean and a variance $2\ln(1.25/\bar{\delta})(\bar{\Delta}^{t,e}_{p,2}/\bar{\epsilon})^2$ and $\bar{\Delta}^{t,e}_{p,2}$ is from \eqref{redefined_L2sensitivity}.
\end{theorem} 

We note that the average randomized output $z^{t}_p$ in line 20 of Algorithm \ref{algo:DP-IADMM-Prox} is the communicated information among agents while the locally updated outputs $\{z_p^{t,1}, \ldots, z_p^{t,E} \}$ are not.
Nevertheless, we randomize all locally updated outputs to prevent data leakage from any one of them, which considers an extreme case when an adversary knows the average randomized output $z^{t}_p$ as well as any $E-1$ locally updated outputs (e.g.,  $\{z_p^{t,1}, \ldots, z_p^{t,E-1} \}$) that can be utilized to compute the  last unknown output (e.g., $z_p^{t,E}$).

Now we study a total privacy leakage of Algorithm \ref{algo:DP-IADMM-Prox}, which can be considered as a $TE$-fold adaptive algorithm.
This involves an extreme situation in which every communication round leaks the output.
According to the existing composition theorem \cite[Theorem 3.16]{dwork2014algorithmic}, Algorithm \ref{algo:DP-IADMM-Prox} guarantees $(TE \bar{\epsilon}, TE \bar{\delta})$-DP because the proposed mechanism \eqref{ADMM-2-Prox-log} guarantees $(\bar{\epsilon}, \bar{\delta})$-DP where $\bar{\delta} \geq 0$. 
This composition theorem implies that data privacy becomes weaker as $T$ and $E$ increase.
We note that this theorem holds for both Laplace and Gaussian mechanisms.

For the Gaussian mechanism, on the other hand, Abadi et al.~\cite{abadi2016deep} proposed a stronger composition theorem based on the moments accountant method, which states that a $T$-adaptive algorithm that guarantees $(\bar{\epsilon}, \bar{\delta})$-DP with $\bar{\delta} > 0$ for every iteration by the existing Gaussian mechanism in Lemma \ref{lem:existing-gaussian} is $(\mathcal{O}(\sqrt{T}) \bar{\epsilon}, \bar{\delta})$-DP, which is a tighter bound than $(T \bar{\epsilon}, T\bar{\delta})$-DP from the composition theorem \cite[Theorem 3.16]{dwork2014algorithmic}.
In the following theorem, we utilize the moments accountant method in \cite{abadi2016deep} to show that our algorithm with the Gaussian mechanism proposed in Theorem \ref{thm:gaussian} provides a tighter bound. 
\begin{theorem}
Consider Algorithm \ref{algo:DP-IADMM-Prox} where \eqref{DPADMM-2-Prox} in line 18 is replaced with \eqref{ADMM-2-Prox-log} where $\ell$ is sufficiently large such that $z^{t,e+1}_{p \ell} \approx z^{t,e+1}_p $ as shown in Theorem \ref{thm:pointwise_convergence}.
This algorithm guarantees $(\sqrt{\frac{TE \ln (1/\bar{\delta})}{\ln(1.25/\bar{\delta})}} \bar{\epsilon}, \bar{\delta} )$-DP if \eqref{ADMM-2-Prox-log} is $(\bar{\epsilon}, \bar{\delta})$-DP, as shown in Theorem \ref{thm:gaussian}.
\end{theorem}
\begin{proof}
\begin{subequations}
We develop the proof based on \cite[Theorem 2]{abadi2016deep}, which is restated in Lemma \ref{lem:composition} in this paper. 
Since we consider an algorithm composed of a sequence of $TE$ adaptive Gaussian mechanisms, we rewrite \eqref{lem:composition-2} within our context, for fixed $t$, $e$, $p$, and sufficiently large $\ell$, as
\begin{align}
  \alpha^{t,e}_p (\tau) := & \max_{\forall (\mathcal{D}_p, \mathcal{D}_p')} \ln \big( \mathbb{E} [ x^{\tau} ] \big),  \ \ \ \
  x := \frac{\textbf{pdf} \big( z^{t,e+1}_{p \ell} (\mathcal{D}_p) = \psi \big)}{\textbf{pdf} \big( z^{t,e+1}_{p \ell} (\mathcal{D}'_p ) = \psi \big)},
\end{align} 
where 
$(\mathcal{D}_p, \mathcal{D}_p')$ are two neighboring datasets, 
$z^{t,e+1}_{p \ell}$ is the optimal solution of \eqref{ADMM-2-Prox-log},
and
$x$ is a random variable that follows a log-normal distribution because in Theorem \ref{thm:gaussian} we show that $\ln(x)$ is a random variable that follows a Gaussian distribution whose mean $\mu$ and variance $\sigma^2$ are presented in \eqref{privacy_loss_rv_mean_var}.
Noting that the $\tau$th moment of the random variable $x$ is given by $\mathbb{E}[x^{\tau}] = \exp( \tau \mu + (\tau^2/2) \sigma^2 )$, we derive
\begin{align}
  \alpha^{t,e}_p(\tau) = & \max_{\forall (\mathcal{D}_p, \mathcal{D}_p')} \ln  ( \mathbb{E}[ x^{\tau}] )  =  \max_{\forall (\mathcal{D}_p, \mathcal{D}_p')}    (\tau+\tau^2) \Big( \frac{\|f'_p(z_p^{t,e};\mathcal{D}_p) - f'_p(z_p^{t,e};\mathcal{D}_p') \|^2}{2 (\sigma^{t,e}_p)^2}  \Big) \nonumber \\
  = &\frac{\tau (\tau+1) (\bar{\Delta}^{t,e}_{p,2})^2}{2 (\sigma^{t,e}_{p})^2} = \frac{\tau(\tau+1) \bar{\epsilon}^2}{4 \ln (1.25/\bar{\delta})}, \label{def_moment}
\end{align}
where the derivation in \eqref{def_moment} holds because of \eqref{redefined_L2sensitivity} and a variance $(\sigma^{t,e}_p)^2 := 2 \ln(1.25/\bar{\delta}) (\bar{\Delta}^{t,e}_{p,2} / \bar{\epsilon})^2$.

Next we restate \eqref{lem:composition-1} within our context as follows:
\begin{align}
  \ln(\bar{\delta}) = \ &  \min_{\tau \in \mathbb{Z}_+} \sum_{t=1}^T \sum_{e=1}^E \alpha^{t,e}_p(\tau)  - \tau \epsilon' = \min_{\tau \in \mathbb{Z}_+} \tau ( a \tau + (a-\epsilon') ), \label{tau_problem} \\
  \text{ where } a := \ &  \frac{TE \bar{\epsilon}^2}{4 \ln (1.25/\bar{\delta})}.
\end{align}
We note that the objective function $\tau ( a \tau + (a-\epsilon') )$ in \eqref{tau_problem} is strictly convex because $a > 0$, whose value is zero at $\tau = 0$ and $\tau = -(a-\epsilon')/a$. 
Since $\ln (\bar{\delta}) < 0$ as $\bar{\delta} \in (0,1)$, the optimal value should be negative; and this is possible when $\tau = 1$ is a feasible solution to the optimization problem in \eqref{tau_problem}.
To this end, we should have
\begin{align}
\frac{\epsilon'-a}{a} > 1  \Leftrightarrow  \epsilon' > \frac{TE \bar{\epsilon}^2}{2 \ln (1.25/\bar{\delta})}. \label{thm:composition-0}
\end{align} 
Now consider
\begin{align}
  \ln(\bar{\delta}) = \ &  \min_{\tau \in \mathbb{Z}_+} \sum_{t=1}^T \sum_{e=1}^E \alpha^{t,e}_p(\tau)  - \tau \epsilon'  \nonumber \\
  \geq \ & \min_{\tau \in \mathbb{R}} \sum_{t=1}^T \sum_{e=1}^E \alpha^{t,e}_p(\tau) - \tau \epsilon' = - \frac{(\epsilon')^2 \ln(1.25/\bar{\delta})}{TE \bar{\epsilon}^2} - \frac{TE \bar{\epsilon}^2}{16 \ln(1.25/\bar{\delta})} + \frac{\epsilon'}{2} \label{thm:composition-1} \\
  \geq \ & - \frac{(\epsilon')^2 \ln(1.25/\bar{\delta})}{TE \bar{\epsilon}^2} - \frac{\epsilon'}{8}  + \frac{\epsilon'}{2}  \geq  - \frac{(\epsilon')^2 \ln(1.25/\bar{\delta})}{TE \bar{\epsilon}^2},\label{thm:composition-2}
\end{align}
where 
the inequality in \eqref{thm:composition-1} holds as the integrality restriction is relaxed;
the equality in \eqref{thm:composition-1} holds because $\tau^* = (\epsilon'-b) / (2b)$, where $b = \bar{\epsilon}^2 / (4 \ln (1.25/\bar{\delta}))$, by the optimality condition; and
the first inequality in \eqref{thm:composition-2} holds because of \eqref{thm:composition-0}.
Therefore we have
\begin{align}
  \bar{\epsilon} \sqrt{ \frac{TE \ln (1/ \bar{\delta})}{\ln(1.25/\bar{\delta})} } \leq \epsilon'.
\end{align}
\end{subequations}
\end{proof}


\section{Convergence analysis} \label{sec:convergence}
In this section we show that a sequence of iterates generated by Algorithm \ref{algo:DP-IADMM-Prox} converges to an optimal solution of \eqref{model:dist} in expectation and that its associated convergence rate is sublinear without any error bound, unlike \cite{huang2020differentially} that results in a nontrivial positive error bound as stated in Section \ref{sec:baseline}.
Moreover, our results characterize the impact of privacy budget parameter $\bar{\epsilon}$ and the number of local updates $E$ as well as the number of ADMM rounds $T$ on the convergence rates.

\subsection{Assumptions}
We make the following assumptions:
\begin{assumption}\label{assump:convergence}  
  \noindent
  \begin{enumerate}
    \item[(i)] The ADMM penalty parameter $\rho^t$ is nondecreasing and bounded; that is, $\rho^1 \leq \rho^2 \leq \ldots \leq \rho^T \leq \rho^{\text{max}}$.
    \item[(ii)] There exists $\gamma > 0$ such that $\gamma \geq 2\| \lambda^*\|$, where $\lambda^*$ is a dual optimal. 
    \item[(iii)] The local objective function $f_p$ is $H$-Lipschitz over a set $\mathcal{W}_p$ with respect to the Euclidean norm.    
  \end{enumerate}  
\end{assumption}
Assumption \ref{assump:convergence} is typically used for the convergence analysis of ADMM (see Chapter 15 of \cite{beck2017first}). 
We adopt the assumptions because the proposed algorithm is a variant of ADMM.
Specifically, Algorithm \ref{algo:DP-IADMM-Prox} with $E=1$ and $\tilde{\xi}^{t,1}_p=0$ for all $t$ is equivalent to the linearized ADMM.
Under Assumption \ref{assump:convergence}, we show that the rate of convergence in expectation is $\mathcal{O}( 1 / (\bar{\epsilon}^2 E \sqrt{T}))$ (resp., $\mathcal{O}( 1 / (\bar{\epsilon} E \sqrt{T}))$) when the local objective function $f_p$ is smooth (resp. nonsmooth).

In addition, we show that the rate becomes $\mathcal{O}( 1 / (\bar{\epsilon}^2 E T))$ when the local objective function $f_p$ is strongly convex. 
This result requires additional assumptions:
\begin{assumption}\label{assump:convergence_stronglyconvex}
\noindent
\begin{enumerate}
\item[(i)] There exists $\gamma > 0$ such that $\| \lambda^t \| \leq \gamma, \ \forall t$.
\item[(ii)] $\rho^t \leq \frac{t}{t-1} \rho^{t-1}, \ \forall t$.
\end{enumerate}
\end{assumption}
Assumption \ref{assump:convergence_stronglyconvex} (i) can be strict in practice. As indicated in \cite{azadi2014towards}, however, it can be considered as a price that we have to pay for faster convergence (see \cite[Assumption 3]{azadi2014towards}).
Assumption \ref{assump:convergence_stronglyconvex} (ii) is not strict since it is satisfied by a constant penalty (i.e., $\rho^t = \rho, \ \forall t$) that is commonly considered in the literature (e.g.,~\cite{beck2017first}).

\subsection{Main results}
In this section we present the convergence analysis results in Theorems \ref{thm:smooth}, \ref{thm:nonsmooth}, and \ref{thm:strong} when the local objective function $f_p$ is smooth, nonsmooth, and strongly convex, respectively; we relegate all the detailed proofs to the Appendix.

\begin{theorem} \label{thm:smooth}
  Suppose that Assumption \ref{assump:convergence} holds, the local objective function $f_p$ is $L$-smooth convex, and the parameter $\eta^t$ of \eqref{DPADMM-2-Prox} is updated as
  \begin{subequations}    
  \begin{align}
    \eta^t  = \frac{1}{L + \sqrt{t}/\bar{\epsilon}}, \ \forall t, \label{smooth_proximity}
  \end{align}
  where $\bar{\epsilon} > 0$ is a privacy budget parameter.
  Then we have  
  \begin{align}
    & \mathbb{E} \Big[ F(z^{(T)}) - F(z^*) + \gamma \| Aw^{(T)} - z^{(T)} \| \Big]   \label{smooth_inequality} \\
    & \leq \frac{nPU_3 + U_2^2/(2E) }{ \bar{\epsilon} \sqrt{T} } + \frac{U_2^2  (\rho^{\text{max}} +L/E) +  (\gamma + \|\lambda^1\|)^2 / \rho^1 }{2T}, \nonumber
  \end{align}    
  where  
  \begin{align}    
    & z := [z_1^{\top}, z_2^{\top}, \ldots, z_P^{\top}]^{\top}, \ \ z^{(T)} :=  \frac{1}{TE} \sum_{t=1}^T \sum_{e=1}^E z^{t,e+1}, \ \ F(z) :=  \sum_{p=1}^P f_p(z_p), \label{Common_Notation} \\ 
    & A^{\top} := \begin{bmatrix}
        \mathbb{I}_J \ \cdots \  \mathbb{I}_J        
      \end{bmatrix}_{J \times PJ}, \ \ w^{(T)} :=  \frac{1}{T} \sum_{t=1}^T w^{t+1}, \nonumber \\
    & U_1 := \max_{u \in \mathcal{W}_p, p \in [P]} \| f'_p(u) \|, \ \ U_2 := \max_{u, v \in \mathcal{W}_p, p \in [P] } \|u-v\|, \nonumber \\
    & U_3 := 
    \begin{cases}
    \max_{p \in [P]}  2 ( \bar{\Delta}^{t,e}_{p,1} )^2  & \text{ for the Laplace mechanism} \\
    \max_{p \in [P]}  2 \ln(1.25/\bar{\delta}) ( \bar{\Delta}^{t,e}_{p,2}  )^2  &   \text{ for the Gaussian mechanism} 
    \end{cases}    
    \nonumber
  \end{align}  
  and $z^*$ is an optimal solution of \eqref{model:dist_1}.
  \end{subequations}  
\end{theorem}
\begin{proof}
See Appendix \ref{apx-thm:smooth}.
\end{proof}
The inequality \eqref{smooth_inequality} implies that the rate of convergence in expectation is 
$\mathcal{O}(1/ (\bar{\epsilon} \sqrt{T}))$ for $\bar{\epsilon} \in (0, \infty)$, 
while in a nonprivate setting it is $\mathcal{O}(1/T)$ because the first term in \eqref{smooth_inequality} is zero when $\bar{\epsilon}=\infty$.
  

\begin{theorem} \label{thm:nonsmooth}
Suppose that Assumption \ref{assump:convergence} holds, the local objective function $f_p$ is a nonsmooth convex function, and the parameter $\eta^t$ of \eqref{DPADMM-2-Prox} is updated as
\begin{subequations}
\begin{align}
\eta^{t} = \frac{1}{\sqrt{t}}, \ \forall t.  \label{nonsmooth_proximity}
\end{align}
Then we have
\begin{align}
& \mathbb{E} \Big[ F(z^{(T)}) - F(z^*) + \gamma \| Aw^{(T)} - z^{(T)} \| \Big] \leq  \label{nonsmooth_inequality} \\
&  \frac{ nPU_3 /\bar{\epsilon}^2 + PU_1^2 + U_2^2 / (2E) }{ \sqrt{T} } + \frac{U_2^2 \rho^{\text{max}} +  (\gamma + \|\lambda^1\|)^2 / \rho^1 + 2\gamma U_2}{2T}, \nonumber
\end{align}  
where $z$, $F(z)$, $A$, $w^{(T)}$, $U_1$, $U_2$, $U_3$ are from \eqref{Common_Notation} while we define
\begin{align}
z^{(T)} :=  \frac{1}{TE} \sum_{t=1}^T \sum_{e=1}^E z^{t,e}. \label{nonsmooth_z}
\end{align}
\end{subequations}
\end{theorem}
\begin{proof}
  See Appendix \ref{apx-thm:nonsmooth}.
\end{proof}
The inequality \eqref{nonsmooth_inequality} implies that the rate of convergence in expectation is 
$\mathcal{O}(1/ (\bar{\epsilon}^2 \sqrt{T}))$ for $\bar{\epsilon} \in (0, \infty)$, 
while in a nonprivate setting it is $\mathcal{O}(1/\sqrt{T})$ because the first term in \eqref{nonsmooth_inequality} is $(PU_1^2 + U_2^2/(2E))/\sqrt{T}$ when $\bar{\epsilon}=\infty$. 
 
\begin{theorem} \label{thm:strong}
Suppose that Assumptions \ref{assump:convergence} and \ref{assump:convergence_stronglyconvex} hold, the local objective function $f_p$ is $\alpha$-strongly convex, and the parameter $\eta^t$ of \eqref{DPADMM-2-Prox} is updated as
\begin{subequations}
\begin{align}
\eta^{t} = \frac{2}{\alpha(t+2)}, \ \forall t. 
\end{align}
Then we have
\begin{align}
  & \mathbb{E} \Big[ F(z^{(T)}) - F(z^*) + \gamma \| Aw^{(T)} - z^{(T)} \| \Big] \leq  \label{strong_inequality} \\
  &\frac{1}{T+1} \Big\{2 U_2 \gamma + U_2^2 \rho^{\text{max}} + 4 \gamma^2 / \rho^1 + \alpha U_2^2/(2E)  + 2P(U_1^2 + nU_3/\bar{\epsilon}^2)/\alpha \Big\}, \nonumber
\end{align}
where $z$, $F(z)$, $A$, $w^{(T)}$, $U_1$, $U_2$, $U_3$ are from \eqref{Common_Notation} while we define
\begin{align}
  & z^{(T)} := \frac{2}{T(T+1)} \sum_{t=1}^T t  (\frac{1}{E} \sum_{e=1}^E z^{t,e}), \ \ w^{(T)} := \frac{2}{T(T+1)} \sum_{t=1}^T t w^{t+1}. \label{strong_z}   
\end{align}  
\end{subequations}
\end{theorem}
\begin{proof}
  See Appendix \ref{apx-thm:strong}. 
\end{proof}
The inequality \eqref{strong_inequality} derived under Assumptions \ref{assump:convergence} and \ref{assump:convergence_stronglyconvex} implies that the rate of convergence in expectation is 
$\mathcal{O}(1/ (\bar{\epsilon}^2 T))$ for $\bar{\epsilon} \in (0, \infty)$, 
while in a nonprivate setting it is $\mathcal{O}(1/T)$ when $\bar{\epsilon}=\infty$. 
  
From the convergence analysis we note that increasing the number $E$ of local updates decreases the values on the right-hand side of \eqref{smooth_inequality}, \eqref{nonsmooth_inequality}, and \eqref{strong_inequality}.
This implies that the gap between $F(z^{(T)})$ and $F(z^*)$ can become smaller by increasing $E$ for fixed $T$.
This may result in better performance by introducing the multiple local updates, thus reducing communication costs, as will be numerically demonstrated in Section \ref{sec:experiments}.

  
 



\section{Numerical demonstration} \label{sec:experiments}
In this section we  numerically demonstrate the effectiveness of the proposed algorithm based on the objective perturbation and multiple local updates for distributed optimization models with constraints.
Specifically, we compare our approach with the popular output perturbation method, which has been widely used in the literature \cite{dvorkin2020differentially, huang2019dp, huang2020differentially}.
For ease of exposition, we denote by $\texttt{ObjL}$ and $\texttt{ObjG}$ the proposed algorithm with Laplace and Gaussian mechanisms, respectively, which will be compared with $\texttt{OutL}$ and $\texttt{OutG}$, the output perturbation methods.
For the comparison of the algorithms, we consider two applications: distributed control of power flow in Section \ref{sec:numerical_power_flow} and federated learning in Section \ref{sec:numerical_fl}.

\subsection{Privacy-preserving distributed control of power flow} \label{sec:numerical_power_flow}
In power systems, an increasing penetration of distributed energy resources in the network has motivated the distributed control of power flow.
Such a problem is a form of \eqref{model:dist}.
We specifically consider that the power network can be decomposed into several zones, each of which is controlled by the agent $p \in [P]$. 
These agents cooperate to find a global power flow solution $w^*$ by iteratively providing a local power flow solution resulting from the model composed of local objective function $f_p$ and local constraints $\mathcal{W}_p$.

In particular we consider an optimization problem that determines a distributed control of power flow such that the deviation of power balance is minimized as considered in \cite{holzer2021grid, dandurand2021bilevel, kim2016data}.
Such a problem can be formulated with the following local objective function:
\begin{align}
f_p (z_p) := \sum_{i \in \mathcal{N}_p} (a_{pi}^{\top} z_{pi} + d_i)^2, \ \forall p \in [P] , \label{localobj}
\end{align}
where $\mathcal{N}_p$ is a set of buses controlled by agent $p$, $d_i$ is a given demand data, $z_{pi}$ is a vector composed of power flow and generation variables connected to a node $i \in \mathcal{N}_p$, and $a_{pi}$ is a given coefficient vector whose element has a value in $[-1,1]$.
Note that the  objective function is derived from the power balance equations $a_{pi}^{\top} z_{pi} = d_i$ for all $i \in \mathcal{N}_p$. 
The local constraint $\mathcal{W}_p$ is composed of convex relaxation constraints for observing some physical laws and operational constraints.
Detailed formulation is presented in Appendix \ref{apx:distpf}.

In this setting, agents are not required to share electric power loads consumed at their zones, but it is still possible to reverse-engineer local solutions communicated between agents to estimate the sensitive power load data \cite{ryu2021privacy}.
To protect data, one can utilize DP algorithms that admit an inevitable trade-off between data privacy and solution quality.
We aim to show that our algorithms \texttt{ObjG} and \texttt{ObjL} provide better solution quality than do the existing algorithms \texttt{OutG} and \texttt{OutL} while providing the same level of data privacy.

\subsubsection{Experimental settings} \label{sec:1_exp_settings}
For the power network instances, we consider case 14 and case 118, each of which is decomposed into three zones as in \cite{ryu2021privacy}.
For the ADMM parameters, we set $\rho^t = 100$ for all $t$ and $\eta^t = 1/\sqrt{t}$.
We consider various $\bar{\epsilon} \in \{0.05, 0.1, 0.5, 1\}$ where smaller $\bar{\epsilon}$ ensures stronger data privacy, as described in Definition \ref{def:differential_privacy_1}. 
For the Gaussian mechanism, we fix $\bar{\delta} = 10^{-2}$ because $\bar{\delta} < n^{-2}$ is typically chosen where $n$ is the number of data points.
For the algorithms based on output perturbation (i.e., \texttt{OutG} and \texttt{OutL}), we compute the sensitivity $\bar{\Delta}_1$ from \eqref{L1_sensitivity} and $\bar{\Delta}_2$ from \eqref{L2_sensitivity} by considering the worst-case scenario as described in \cite[Proposition 1]{dvorkin2020differentially}, resulting in $\bar{\Delta}_1 = \bar{\Delta}_2 = \beta$, where $\beta$ measures the adjacency of two neighboring datasets, namely, $\mathcal{D}_p := \{d_i\}_{i \in \mathcal{N}_z}$ and $\mathcal{D}_p':= \{d_i'\}_{i \in \mathcal{N}_z}$, which are $\beta$-adjacent datasets if  $\exists i \text{ s.t. } |d_i - d_i'| \leq \beta, \ d_j = d_j', \forall j \neq i$.
For a fair comparison, we also compute the worst-case sensitivity for our algorithms \texttt{ObjG} and \texttt{ObjL}, resulting in $\bar{\Delta}_1 = \bar{\Delta}_2 = 2\beta$, which is twice  the sensitivity used in \texttt{OutG} and \texttt{OutL} based on \eqref{redefined_L1sensitivity}, \eqref{redefined_L2sensitivity}, and \eqref{localobj}.
In this experiment we set $\beta = 1 \%$ as in \cite{dvorkin2020differentially}.
For all experiments, we solve optimization models by Ipopt \cite{wachter2006implementation} via Julia 1.8.0 on Bebop, a 1024-node computing cluster at Argonne National Laboratory.
Each computing node has 36 cores with Intel Xeon E5-2695v4 processors and 128 GB DDR4 of memory.



\begin{figure*}[!ht]  
\centering
\begin{subfigure}[b]{0.23\textwidth}
\centering        
\includegraphics[width=\textwidth]{./Figures_1/case14_noise_eps_0.05}  
\includegraphics[width=\textwidth]{./Figures_1/case14_objval_mlu_1_eps_0.05}   
\includegraphics[width=\textwidth]{./Figures_1/case14_objval_mlu_5_eps_0.05}   
\includegraphics[width=\textwidth]{./Figures_1/case14_prim_res_object_gaussian_eps_0.05_mlu}  
\caption{$\bar{\epsilon}=0.05$}
\end{subfigure}
\begin{subfigure}[b]{0.23\textwidth}
\centering        
\includegraphics[width=\textwidth]{./Figures_1/case14_noise_eps_0.1}  
\includegraphics[width=\textwidth]{./Figures_1/case14_objval_mlu_1_eps_0.1}   
\includegraphics[width=\textwidth]{./Figures_1/case14_objval_mlu_5_eps_0.1}   
\includegraphics[width=\textwidth]{./Figures_1/case14_prim_res_object_gaussian_eps_0.1_mlu}  
\caption{$\bar{\epsilon}=0.1$}
\end{subfigure}
\begin{subfigure}[b]{0.23\textwidth}
\centering      
\includegraphics[width=\textwidth]{./Figures_1/case14_noise_eps_0.5}  
\includegraphics[width=\textwidth]{./Figures_1/case14_objval_mlu_1_eps_0.5}   
\includegraphics[width=\textwidth]{./Figures_1/case14_objval_mlu_5_eps_0.5}   
\includegraphics[width=\textwidth]{./Figures_1/case14_prim_res_object_gaussian_eps_0.5_mlu}  
\caption{$\bar{\epsilon}=0.5$}
\end{subfigure}     
\begin{subfigure}[b]{0.23\textwidth}
\centering      
\includegraphics[width=\textwidth]{./Figures_1/case14_noise_eps_1.0}  
\includegraphics[width=\textwidth]{./Figures_1/case14_objval_mlu_1_eps_1.0}   
\includegraphics[width=\textwidth]{./Figures_1/case14_objval_mlu_5_eps_1.0}   
\includegraphics[width=\textwidth]{./Figures_1/case14_prim_res_object_gaussian_eps_1.0_mlu}  
\caption{$\bar{\epsilon}=1$}
\end{subfigure}       
\caption{Comparison of \texttt{ObjG}, \texttt{ObjL}, \texttt{OutG}, \texttt{OutL} in terms of noise magnitude ($1^{\text{st}}$ row), objective value ($2^{\text{nd}}$ and $3^{\text{rd}}$ rows), and feasibility ($4^{\text{th}}$ row) by using case 14.}
\label{fig:case14}
\end{figure*}  




\begin{figure*}[!ht]  
\centering
\begin{subfigure}[b]{0.23\textwidth}
\centering        
\includegraphics[width=\textwidth]{./Figures_1/case118_noise_eps_0.05}  
\includegraphics[width=\textwidth]{./Figures_1/case118_objval_mlu_1_eps_0.05}   
\includegraphics[width=\textwidth]{./Figures_1/case118_objval_mlu_5_eps_0.05}   
\includegraphics[width=\textwidth]{./Figures_1/case118_prim_res_object_gaussian_eps_0.05_mlu}  
\caption{$\bar{\epsilon}=0.05$}
\end{subfigure}
\begin{subfigure}[b]{0.23\textwidth}
\centering        
\includegraphics[width=\textwidth]{./Figures_1/case118_noise_eps_0.1}  
\includegraphics[width=\textwidth]{./Figures_1/case118_objval_mlu_1_eps_0.1}   
\includegraphics[width=\textwidth]{./Figures_1/case118_objval_mlu_5_eps_0.1}   
\includegraphics[width=\textwidth]{./Figures_1/case118_prim_res_object_gaussian_eps_0.1_mlu}  
\caption{$\bar{\epsilon}=0.1$}
\end{subfigure}
\begin{subfigure}[b]{0.23\textwidth}
\centering      
\includegraphics[width=\textwidth]{./Figures_1/case118_noise_eps_0.5}  
\includegraphics[width=\textwidth]{./Figures_1/case118_objval_mlu_1_eps_0.5}   
\includegraphics[width=\textwidth]{./Figures_1/case118_objval_mlu_5_eps_0.5}   
\includegraphics[width=\textwidth]{./Figures_1/case118_prim_res_object_gaussian_eps_0.5_mlu}  
\caption{$\bar{\epsilon}=0.5$}
\end{subfigure}     
\begin{subfigure}[b]{0.23\textwidth}
\centering      
\includegraphics[width=\textwidth]{./Figures_1/case118_noise_eps_1.0}  
\includegraphics[width=\textwidth]{./Figures_1/case118_objval_mlu_1_eps_1.0}   
\includegraphics[width=\textwidth]{./Figures_1/case118_objval_mlu_5_eps_1.0}   
\includegraphics[width=\textwidth]{./Figures_1/case118_prim_res_object_gaussian_eps_1.0_mlu}  
\caption{$\bar{\epsilon}=1$}
\end{subfigure}       
\caption{Comparison of \texttt{ObjG}, \texttt{ObjL}, \texttt{OutG}, \texttt{OutL} in terms of noise magnitude ($1^{\text{st}}$ row), objective value ($2^{\text{nd}}$ and $3^{\text{rd}}$ rows), and feasibility ($4^{\text{th}}$ row) by using case 118.}
\label{fig:case118}
\end{figure*}  

\subsubsection{Performance comparison}
We report the performance of \texttt{ObjG}, \texttt{ObjL}, \texttt{OutG}, and \texttt{OutL} for case 14 and case 118 instances in Figure \ref{fig:case14} and \ref{fig:case118}, respectively. 

First, we compare the four algorithms with respect to the magnitude of noise introduced during $T=100$ ADMM rounds, which is defined as
\begin{align}
  \sum_{p=1}^P  \sum_{e=1}^E  \sum_{i=1}^n | \tilde{\xi}^{t,e}_{p i} |, \ \forall t \in [T], \label{form:noise_mag}
\end{align}
where $\tilde{\xi}^{t,e}_{pi} \in \mathbb{R}$ is sampled from either Laplace or Gaussian distribution.
The noise magnitude \eqref{form:noise_mag} under various $\bar{\epsilon}$ is reported in figures from the first row of Figure \ref{fig:case14} and \ref{fig:case118}.
From these figures, we make the following observations:
\begin{enumerate}
  \item The noise magnitude increases as $\bar{\epsilon}$ decreases.
  \item For fixed $\bar{\epsilon}$, the proposed algorithm \texttt{ObjG} (resp., \texttt{ObjL}) requires more noise compared with \texttt{OutG} (resp., \texttt{OutL}).
  \item For fixed $\bar{\epsilon}$, algorithms with the Gaussian mechanism (i.e., \texttt{ObjG} and \texttt{OutG}) require more noise compared with algorithms with the Laplace mechanism (i.e., \texttt{ObjL} and \texttt{OutL}).
\end{enumerate} 
To explain these observations, we report in Table \ref{table:var} the variance of distributions used for sampling each element of $\tilde{\xi}^{t,e}_p$.
For all algorithms considered in the table, decreasing $\bar{\epsilon}$ increases the variance, which can explain the first observation.
From the table we note that \texttt{ObjG} has the largest variance, followed by \texttt{OutG}, \texttt{ObjL}, and \texttt{OutL}. 
These explain the second and third observations.

\begin{table}[H]
\centering
\caption{Variance of distributions used for sampling noise in the four algorithms.}
\begin{tabular}{ c c c c}
\hline
\texttt{ObjG} & \texttt{ObjL}  & \texttt{OutG}  & \texttt{OutL} \\ \hline    
$2 \ln ( \frac{1.25}{\bar{\delta}}) (\frac{2\beta}{\bar{\epsilon}})^2$ 
& $2  (\frac{2\beta}{\bar{\epsilon}})^2$ 
& $2 \ln ( \frac{1.25}{\bar{\delta}}) (\frac{\beta}{\bar{\epsilon}})^2$ 
& $2  (\frac{\beta}{\bar{\epsilon}})^2$   \\
\hline
\end{tabular}
\label{table:var}
\end{table}
  


Second, we report the objective function values (which should be zero at optimality in a non-private setting) produced by the four algorithms under various $\bar{\epsilon}$ in the figures located at the second row of Figures \ref{fig:case14} and \ref{fig:case118}, when $E=1$, namely, single local update.
We notice that the proposed algorithms \texttt{ObjG} and \texttt{ObjL} produce lower objective function values compared with \texttt{OutG} and \texttt{OutL} even though our algorithm requires more noise to be introduced for guaranteeing the same level of data privacy. 
This shows the effectiveness of our approach, which ensures feasibility of the intermediate solutions randomized for DP.
We also note that \texttt{ObjL} performs better than \texttt{ObjG}. The reason is  that \texttt{ObjL} requires less noise in this experiment.

Third, we set the number of local updates to $E=5$ and report the objective values in the figures located at the third row of Figures \ref{fig:case14} and \ref{fig:case118}. For larger $\bar{\epsilon}$, \texttt{OutG} and \texttt{OutL} sometimes produce better objective function values. Clearly, however,  our algorithms provide better objective values for smaller $\bar{\epsilon}$ (i.e., stronger data privacy).
The reason is mainly that there is no $\bar{\epsilon}$-dependent error bound in our convergence rate presented in Section \ref{sec:convergence} while the error bound exists for the existing DP algorithms \cite{huang2020differentially} based on the output perturbation, which increases as $\bar{\epsilon}$ decreases. 

Fourth, we report feasibility, namely, $\|Aw^t+Bz^t\|_2$, where $w$ is a global solution, $z$ is a local solution, and $Aw+Bz$ represents the consensus constraints (which should be zero at optimality in a non-private setting), produced by \texttt{ObjG} for $E \in \{1,5\}$ in figures located at the last row of Figures \ref{fig:case14} and \ref{fig:case118}. 
We observe that increasing the number of local updates provides a solution with better quality. 
Similar behavior is observed in other algorithms \texttt{ObjL}, \texttt{OutG}, and \texttt{OutL},  numerically demonstrating the effectiveness of introducing the multiple local update.


\subsection{Federated learning} \label{sec:numerical_fl}
In ML, transferring data into a central server for training can be limited because of, for example, low bandwidth or data privacy issues. Such limitations have motivated the use of distributed optimization, also known as federated learning. 
In most FL literature, a distributed empirical risk minimization model is utilized that is a form of \eqref{model:dist}.
Specifically, there are multiple $P$ agents that have their own data and machine for local training. 
These agents cooperate to find a global model parameter $w^*$ by iteratively training using the local objective function $f_p (w) := (\frac{I_p}{I}) \frac{1}{I_p}\sum_{i \in \mathcal{I}_p}  l (w; x_i, y_i)$, where $\mathcal{I}_p$ is an index set of local data samples, $I_p := |\mathcal{I}_p|$ is the number of local data samples, $I:=\sum_{p=1}^P I_p$ is the total number of data samples, 
$l$ represents a loss function, 
and $x_i$ and $y_i$ represent data features and data labels, respectively.
We note that most FL literature assumes $\mathcal{W} := \mathbb{R}^n$, but more recently the importance of imposing hard constraints instead of soft constraints (i.e., penalizing in the objective function) has been discussed in the ML community \cite{gallego2022controlled, marquez2017imposing}. 
 
In particular, we consider a multiclass logistic regression model for classifying image data.
Such problems can be formulated as the form of \eqref{model:dist_1} with the following local objective function:
\begin{align}
f_p (z_p) :=  - \frac{1}{I} \sum_{i=1}^{I_p} \sum_{k=1}^K y_{pik} \ln \Big(\frac{\exp(\sum_{j=1}^J  x_{pij} z_{pjk}) }{\sum_{k'=1}^K \exp(\sum_{j=1}^J  x_{pij} z_{pjk'}) } \Big), \ \forall p \in [P], \label{FL_local_obj}
\end{align}
where $z_p \in \mathbb{R}^{J \times K}$ represents local model parameters and $x_p \in \mathbb{R}^{I_p \times J}$ and $y_p \in \mathbb{R}^{I_p \times K}$ are data features and labels, respectively. Note that for each local data sample $i \in [I_p]$, there are $J$ data features and $K$ data labels.
In this experiment we consider a simple local constraint $\mathcal{W}_p :=  [-u, u]$ for some $u \in \mathbb{R}^{J \times K}_{++}$. 

In this setting agents are not required to share their local image data with other agents. 
 One can, however,  reconstruct the law data if local model parameters are exposed.
To protect data, one can utilize a DP algorithms if one can admit the inevitable trade-off between data privacy and learning performance. 
In this experiment we aim to show that our algorithm \texttt{ObjG} outperforms the existing algorithm \texttt{OutG} while ensuring the same level of data privacy under the existence of convex constraints. 
We note that \texttt{OutG} proposed in \cite{huang2019dp, huang2020differentially} demonstrated better learning performance compared with other existing DP algorithms, such as DP-SGD in \cite{abadi2016deep} and DP-ADMM in \cite{zhang2016dynamic}, under the same level of data privacy, when there are no constraints to satisfy.

\subsubsection{Experimental settings} \label{sec:numerical_fl_exp}
We consider two publicly available datasets for image classification: MNIST \cite{lecun1998mnist} and FEMNIST \cite{caldas2018leaf}.
For the MNIST dataset, we split the 60,000 training data points over $P=10$ agents, each of which is assigned to have the same number of IID datasets.
For the FEMNIST dataset, we follow the preprocess procedure\footnote{https://github.com/TalwalkarLab/leaf/tree/master/data/femnist} to sample 5\% of the entire 805,263 data points in a non-IID manner, resulting in 36,708 training samples distributed over $P=195$ agents.
For local constraints $\mathcal{W}_p :=  [-u, u]$, we set $u=0.1$ for MNIST and $u=1.0$ for FEMNIST.
We set $\bar{\delta}=10^{-6}$ and consider various $\bar{\epsilon} \in \{0.05, 0.1, 0.5, 1\}$ where smaller $\bar{\epsilon}$ ensures stronger data privacy.
For \texttt{OutG}, we compute the $L_2$ sensitivity of the form \eqref{L2_sensitivity} as discussed in \cite{huang2019dp}.
We note that the $L_2$ sensitivity in \cite{huang2019dp} is derived from the optimality condition when $\mathcal{W}_p:=\mathbb{R}^n$.
It can be utilized for the simple box constraint considered in this experiment, but it cannot be utilized straightforwardly for a general convex $\mathcal{W}_p$. 
For \texttt{ObjG}, we compute the $L_2$ sensitivity of the form \eqref{redefined_L2sensitivity}.
Specifically, we add a data sample $(x_{pi^*}, y_{pi^*})$ to the existing dataset $\mathcal{D}_p$ to construct a neighboring dataset $\mathcal{D}_p'$. 
Then we can compute the $L_2$ sensitivity \eqref{redefined_L2sensitivity} based on \eqref{FL_local_obj} as follows:
\begin{align}
\bar{\Delta}^{t,e}_{p,2} = \max_{ \forall (\mathcal{D}_p, \mathcal{D}_p')}  \Big[ \sum_{j=1}^J \sum_{k=1}^K \Big\{ \frac{1}{I+1}   x_{pi^* j} \Big(-y_{pi^*k} + \frac{  \exp(\sum_{j'=1}^J  x_{pi^*j} z^{t,e}_{p j' k})}{\sum_{k'=1}^K \exp(\sum_{j'=1}^J  x_{pi^*j'} z^{t,e}_{pjk'})}  \Big) \Big\}^2 \Big]^{0.5}, \label{fl_sensitivity}
\end{align} 
which is computed by finding $(x_{pi^*}, y_{pi^*})$ from the universe of datasets that maximizes \eqref{fl_sensitivity} for given $z^{t,e}_p$.

A parameter $\rho^t$ is chosen after conducting some tuning process discussed in Appendix \ref{apx:hyper}. We note that the chosen parameter $\hat{\rho}^t$ is nondecreasing and bounded above, thus satisfying Assumption \ref{assump:convergence} (i).

We implemented the algorithms in Python and ran the experiments on Swing, a 6-node GPU computing cluster at Argonne National Laboratory. Each node of Swing has 8 NVIDIA A100 40 GB GPUs, as well as 128 CPU cores.

 


\begin{figure*}[!ht]  
\centering
\begin{subfigure}[b]{0.23\textwidth}
\centering        
\includegraphics[width=\textwidth]{./Figures_2/MNIST_bound_0.1_eps_0.05_noise}
\includegraphics[width=\textwidth]{./Figures_2/MNIST_bound_0.1_eps_0.05_training}  
\includegraphics[width=\textwidth]{./Figures_2/MNIST_bound_0.1_eps_0.05_testing}  
\caption{$\bar{\epsilon}=0.05$}
\end{subfigure}
\begin{subfigure}[b]{0.23\textwidth}
\centering        
\includegraphics[width=\textwidth]{./Figures_2/MNIST_bound_0.1_eps_0.1_noise}
\includegraphics[width=\textwidth]{./Figures_2/MNIST_bound_0.1_eps_0.1_training}  
\includegraphics[width=\textwidth]{./Figures_2/MNIST_bound_0.1_eps_0.1_testing}  
\caption{$\bar{\epsilon}=0.1$}
\end{subfigure}
\begin{subfigure}[b]{0.23\textwidth}
\centering      
\includegraphics[width=\textwidth]{./Figures_2/MNIST_bound_0.1_eps_0.5_noise}
\includegraphics[width=\textwidth]{./Figures_2/MNIST_bound_0.1_eps_0.5_training}  
\includegraphics[width=\textwidth]{./Figures_2/MNIST_bound_0.1_eps_0.5_testing}  
\caption{$\bar{\epsilon}=0.5$}
\end{subfigure}     
\begin{subfigure}[b]{0.23\textwidth}
\centering      
\includegraphics[width=\textwidth]{./Figures_2/MNIST_bound_0.1_eps_1.0_noise}
\includegraphics[width=\textwidth]{./Figures_2/MNIST_bound_0.1_eps_1.0_training}  
\includegraphics[width=\textwidth]{./Figures_2/MNIST_bound_0.1_eps_1.0_testing}  
\caption{$\bar{\epsilon}=1$}
\end{subfigure}           
\caption{Comparison of \texttt{ObjG} and \texttt{OutG} with respect to the noise magnitude ($1^{\text{st}}$ row), objective value ($2^{\text{nd}}$ row), and testing error ($3^{\text{rd}}$ row) by using MNIST.}
\label{fig:MNIST}
\end{figure*}   


\begin{figure*}[!ht]  
  \centering
  \begin{subfigure}[b]{0.23\textwidth}
  \centering        
  \includegraphics[width=\textwidth]{./Figures_2/FEMNIST_bound_1.0_eps_0.05_noise}  
  \includegraphics[width=\textwidth]{./Figures_2/FEMNIST_bound_1.0_eps_0.05_training}  
  \includegraphics[width=\textwidth]{./Figures_2/FEMNIST_bound_1.0_eps_0.05_testing}  
  \caption{$\bar{\epsilon}=0.05$}
  \end{subfigure}
  \begin{subfigure}[b]{0.23\textwidth}
  \centering        
  \includegraphics[width=\textwidth]{./Figures_2/FEMNIST_bound_1.0_eps_0.1_noise}  
  \includegraphics[width=\textwidth]{./Figures_2/FEMNIST_bound_1.0_eps_0.1_training}  
  \includegraphics[width=\textwidth]{./Figures_2/FEMNIST_bound_1.0_eps_0.1_testing}    
  \caption{$\bar{\epsilon}=0.1$}
  \end{subfigure}
  \begin{subfigure}[b]{0.23\textwidth}
  \centering      
  \includegraphics[width=\textwidth]{./Figures_2/FEMNIST_bound_1.0_eps_0.5_noise}  
  \includegraphics[width=\textwidth]{./Figures_2/FEMNIST_bound_1.0_eps_0.5_training}  
  \includegraphics[width=\textwidth]{./Figures_2/FEMNIST_bound_1.0_eps_0.5_testing}  
  \caption{$\bar{\epsilon}=0.5$}
  \end{subfigure}     
  \begin{subfigure}[b]{0.23\textwidth}
  \centering      
  \includegraphics[width=\textwidth]{./Figures_2/FEMNIST_bound_1.0_eps_1.0_noise}  
  \includegraphics[width=\textwidth]{./Figures_2/FEMNIST_bound_1.0_eps_1.0_training}  
  \includegraphics[width=\textwidth]{./Figures_2/FEMNIST_bound_1.0_eps_1.0_testing}  
  \caption{$\bar{\epsilon}=1$}
  \end{subfigure}           
  \caption{Comparison of \texttt{ObjG} and \texttt{OutG} with respect to the noise magnitude ($1^{\text{st}}$ row), objective value ($2^{\text{nd}}$ row), and testing error ($3^{\text{rd}}$ row) by using FEMNIST.}
  \label{fig:FEMNIST}
  \end{figure*}  


\subsubsection{Performance comparison} 
We report the performance of the proposed algorithm \texttt{ObjG} and the existing algorithm \texttt{OutG} for MNIST and FEMNIST image classification in Figure \ref{fig:MNIST} and \ref{fig:FEMNIST}, respectively.

First, we report in the first row of Figures \ref{fig:MNIST} and \ref{fig:FEMNIST} the maginitude of noise introduced during $T=5000$ ADMM rounds, which is defined as
\begin{align}
  \sum_{p=1}^P \sum_{e=1}^E \Big(  \frac{1}{n}  \sum_{i=1}^n | \tilde{\xi}^{t,e}_{p i} | \Big), \ \forall t \in [T],
\end{align}
where $\tilde{\xi}^{t,e}_{pi} \in \mathbb{R}$ is sampled from a Gaussian distribution with zero mean and variance $\sigma^2 = 2 \ln (1.25/\bar{\delta}) ( \bar{\Delta}^{t,e}_{p,2} / \bar{\epsilon})^2$, where $\bar{\Delta}^{t,e}_{p,2}$ is the $L_2$ sensitivity computed as described in Section \ref{sec:numerical_fl_exp}.
We observe that the noise magnitude for \texttt{ObjG} is greater than that for \texttt{OutG} mainly because the sensitivity of \texttt{ObjG} is greater.

Second, we report the objective values (i.e., training costs) produced by \texttt{ObjG} and \texttt{OutG} in the second row of Figures \ref{fig:MNIST} and \ref{fig:FEMNIST}, respectively.
We observe that the objective value increases as $\bar{\epsilon}$ decreases, which indicates the trade-off between data privacy and learning performance, well known in the literature on DP algorithms \cite{dwork2014algorithmic}.
Also, the objective values of \texttt{ObjG} are lower than those of \texttt{OutG}.
Furthermore, increasing $E$ lowers the objective values, showing the efficacy of the multiple local updates.
  
Third, we report the testing errors produced by \texttt{ObjG} and \texttt{OutG} in the third row of Figures \ref{fig:MNIST} and \ref{fig:FEMNIST}, respectively.
We observe that the testing errors of \texttt{ObjP} with multiple local updates significantly outperforms the other algorithms.
This result implies that our algorithm can mitigate the trade-off between data privacy and learning performance.   

\section{Conclusion} \label{sec:conclusion}
In this paper we proposed a differentially private distributed optimization algorithm for general convex optimization problems. 
Our algorithm addresses a data privacy concern in distributed optimization by communicating randomized outputs of the local optimization models, which prevents reconstructing data stored at the local machine.
Different from the existing DP algorithms in the literature, most of which have been developed for unconstrained machine learning models, our algorithm is developed for constrained convex optimization models, which can be beneficial because constraints are necessary in most optimization control problems and which become more important in some machine learning applications.

Specifically, our algorithm generalizes the existing linearized ADMM by introducing multiple local updates to reduce communication costs and by incorporating the objective perturbation method to the local optimization models so that the resulting randomized outputs are always feasible (i.e., satisfying the constraints).
This is different from the popular output perturbation method used in the existing DP algorithms where the resulting randomized outputs constructed by adding noise to the true output, which could be infeasible. 
We presented two different mechanisms, namely, Laplace and Gaussian mechanisms, that ensure differential privacy for every iteration of the proposed algorithm.
These are similar to the Laplace and Gaussian mechanisms used in the output perturbation method except for the sensitivity computation. 
We also provided a convergence analysis that showed that the rate of convergence in expectation is sublinear without any error bound.  
We discussed two different application areas, namely, distributed control of power flow and federated learning, where the proposed algorithms can be utilized.
From the experiments, we numerically demonstrate the outperformance of the proposed algorithm over the existing DP algorithms with the output perturbation method, when there are constraints to satisfy.

We note that the privacy and convergence analyses in this paper are limited to the convex setting and thus not directly applicable to the nonconvex setting. 


































\newpage




\section{A detailed example}

Here we include some equations and theorem-like environments to show
how these are labeled in a supplement and can be referenced from the
main text.
Consider the following equation:
\begin{equation}
  \label{eq:suppa}
  a^2 + b^2 = c^2.
\end{equation}
You can also reference equations such as \cref{eq:matrices,eq:bb} 
from the main article in this supplement.

\lipsum[100-101]

\begin{theorem}
An example theorem.
\end{theorem}

\lipsum[102]
 
\begin{lemma}
An example lemma.
\end{lemma}

\lipsum[103-105]

Here is an example citation: \cite{KoMa14}.

\section[Proof of Thm]{Proof of \cref{thm:bigthm}}
\label{sec:proof}

\lipsum[106-112]

\section{Additional experimental results}
\Cref{tab:smfoo} shows additional
supporting evidence. 

\begin{table}[htbp]
\footnotesize
  \caption{Example table.}\label{tab:smfoo}
\begin{center}
  \begin{tabular}{|c|c|c|} \hline
   Species & \bf Mean & \bf Std.~Dev. \\ \hline
    1 & 3.4 & 1.2 \\
    2 & 5.4 & 0.6 \\ \hline
  \end{tabular}
\end{center}
\end{table}


\bibliographystyle{siamplain}
