# On the convergence of the distributed Proximal point algorithm

Woocheol Choi<EMAIL_ADDRESS>
###### Abstract.

In this work, we establish convergence results for the distributed proximal
point algorithm (DPPA) for distributed optimization problems. We consider the
problem on the whole domain $\mathbb{R}^{d}$ and find a general condition on
the stepsize and cost functions such that the DPPA is stable. We prove that
the DPPA with stepsize $\eta>0$ exponentially converges to an
$O(\eta)$-neighborhood of the optimizer. Our result clearly explains the
advantage of the DPPA with respect to the convergence and stability in
comparison with the distributed gradient descent algorithm. We also provide
numerical tests supporting the theoretical results.

###### Key words and phrases:

Proximal point method, Distributed optimization, Distributed gradient
algorithm

###### 2010 Mathematics Subject Classification:

Primary 90C25, 68Q25

##  1\. Introduction

This works consider the distributed optimization

$\min_{x\in\mathbb{R}^{d}}f(x)=\sum_{i=1}^{n}f_{i}(x),$ (1.1)

where $n$ denotes the number of the agents in the network, and
$f_{i}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a differentiable local cost
only known agent $i$ for each $1\leq i\leq n$. In recent years, extensive
research has been conducted on distributed optimization due to its relevance
in various applications, including wireless sensor networks [1, 24], multi-
agent control [4, 22, 5], smart grids [9, 10], and machine learning [11, 21,
30, 3, 27].

Numerous studies have focused on solving distributed optimization problems on
networks. Relevant research works include [17, 18, 25, 26], and the references
therein. A key algorithm in this field is the distributed gradient descent
algorithm (DGD), introduced in [17]. Additionally, there exist various
versions of distributed optimization algorithms such as EXTRA [25] and
decentralized gradient tracking [18, 20, 28]. Lately, there has been growing
interest in designing communication-efficient algorithms for distributed
optimization [23, 8, 12, 2].

The distributed proximal point algorithm (DPPA) was proposed in [15, 16] as a
distributed counterpart of the proximal point method, analogous to the
relation between the distributed gradient descent algorithm and the gradient
descent method. The works [15, 16] established the asymptotic convergence of
the DPPA under the assumptions of a compact domain for each local cost $f_{i}$
and a decreasing step size. The work [14] designed the DPPA on a directed
graph and proved a convergence estimate with rate $O(1/\sqrt{t})$ when the
step size is set to $1/\sqrt{t}$ and each cost function has a compact domain.

It is a well-known fact that the proximal point method is more stable compared
to the gradient descent method for large choices of the stepsize. This fact
suggests that the DPPA may also exhibit be more stable than the DGD, as
mentioned in the previous works [13, 14, 16]. In this work, we provide
convergence results for the DPPA in the case of the entire domain and a
constant step size. Comparing the results with the convergence result [29] of
the DGD, we will find that the DPPA is more stable than the DGD when the
stepsize is large.

The DPPA for the problem (1.1) is described as follows:

$\begin{split}\hat{x}_{i}(t)&=\sum_{j=1}^{n}w_{ij}x_{j}(t)\\\
x_{i}(t+1)&=\textrm{argmin}_{x\in\mathbb{R}^{n}}\Big{(}f_{i}(x)+\frac{1}{2\eta}\|x-\hat{x}_{i}(t)\|^{2}\Big{)},\end{split}$
(1.2)

where $w_{ij}$ denotes the weight for the communication among agents in the
network desribed by an undirected graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$. Each node in $\mathcal{V}$ represents
an agent, and each edge $\\{i,j\\}\in\mathcal{E}$ means that $i$ can send
messages to $j$ and vice versa. We consider a graph $\mathcal{G}$ satisfying
the following assumption.

######  Assumption 1.

The communication graph $\mathcal{G}$ is undirected and connected, i.e., there
exists a path between any two agents.

We define the mixing matrix $W=\\{w_{ij}\\}_{1\leq i,j\leq n}$ as follows. The
nonnegative weight $w_{ij}$ is given for each communication link
$\\{i,j\\}\in\mathcal{E},$ where $w_{ij}\neq 0$ if $\\{i,j\\}\in\mathcal{E}$
and $w_{ij}=0$ if $\\{i,j\\}\notin\mathcal{E}$. In this paper, we make the
following assumption on the mixing matrix $W$.

######  Assumption 2.

The mixing matrix $W=\\{w_{ij}\\}_{1\leq i,j\leq n}$ is doubly stochastic,
i.e., $W\mathbf{1}=\mathbf{1}$ and $\mathbf{1}^{T}W=\mathbf{1}^{T}$. In
addition, $w_{ii}>0$ for all $i\in\mathcal{V}$.

In the following result, we establish that the sequence
$\\{x_{k}(t)\\}_{k=1}^{n}$ is stable, i.e., uniformly bounded for $t\geq 0$
under suitable conditions.

######  Theorem 1.1.

Suppose that the assumptions 1-2 hold, and assume that one of the following
conditions holds:

  1. (1)

The following function

$F_{\eta}(x):=\sum_{k=1}^{n}f_{k}(x_{k})+\frac{1}{2\eta}x^{T}(I-W)x$ (1.3)

is bounded below and has an optimal solution $(x_{1}^{*},\cdots,x_{n}^{*})$.

  2. (2)

Each function $f_{k}$ is $L$-smooth and $f$ is $\alpha$-stronlgy convex. In
addition, the stepsize $\eta>0$ satisfies the following inequality

$\eta^{2}(\alpha^{2}+\alpha
L)+\eta\Big{(}\alpha+L-\frac{(1-\rho_{W})\alpha^{2}}{L}\Big{)}<\frac{(1-\rho_{W})\alpha}{L},$
(1.4)

where $\rho_{W}$ is the spectral norm of the matrix
$W-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}$.

Then the sequence $\\{x_{k}(t)\\}_{k=1}^{n}$ of the DPPA with stepsize
$\eta>0$ is uniformly bounded for $t\geq 0$.

We now compare the stability result with that of the DGD described by

$x_{i}(t+1)=\sum_{j=1}^{N}w_{ij}x_{j}(t)-\eta\nabla f_{i}(x_{i}(t)).$

Let $\lambda_{n}(W)\in\mathbb{R}$ be the smallest eigenvalue of $W$. Then it
is known from [29] that the DGD is stable if the stepsize $\eta>0$ satisfies

$\eta\leq\frac{1+\lambda_{n}(W)}{L},$

provided that each cost function $f_{j}$ is convex and $L$-smooth.

| Condition on the costs | Stepsize | Paper
---|---|---|---
DGD |  Each $f_{j}$ is convex and $L$-smooth | $\eta\leq\frac{(1+\lambda_{n}(W))}{L}$ | [29]
DPPA |  $F_{\eta}$ is bounded below and has an optimizer | $\eta\in(0,\infty)$ | This work
Table 1. This table compares the condition on the stepsize for the stability
of the DGD and the DPPA.

Table 1 summarizes the condition on the stepsize $\eta>0$ for the stability of
the DGD and the DPPA. We observe that $F_{\eta}$ defined in (1.3) is bounded
below for any $\eta>0$ when each $f_{j}$ is convex. Therefore, the condition
of (1) in Theorem 1.1 is not so restrictive than the condition that $f_{j}$ is
convex which is required for the stability result [29] of the DGD. Hence the
result of Theorem 1.1 proves that the stepsize choice of $\eta>0$ is much
wider for the DPPA than that for the DGD.

For the convergence analysis of the DPPA, we assume the uniformly boundedness
property.

######  Assumption 3.

The sequence $\\{x_{k}(t)\\}_{k=1}^{n}$ is uniformly bounded for $t\geq 0$,
i.e., there is $R>0$ such that

$A_{t}\leq R\quad\textrm{and}\quad B_{t}\leq R\quad\forall~{}t\geq 0.$

Although the property is guaranteed for broad class of functions by the result
of Theorem 1.1, we formulate this assumption of the sake of simplicity in the
statement of the convergence result. We also consider the following
assumption.

######  Assumption 4.

The aggregate cost function $f$ is $\alpha$-strongly convex for some
$\alpha>0$ and each cost function $f_{j}$ is $L$-smooth for each $1\leq j\leq
n$, i.e.,

$\|\nabla f_{j}(x)-\nabla f_{j}(y)\|\leq L\|x-y\|$

for all $x,y\in\mathbb{R}^{d}$.

Under this assumption, there exists a unique optimizer
$x_{*}=\arg\min_{x\in\mathbb{R}^{d}}f(x)$. We let $D=\max_{1\leq i\leq
n}\|\nabla f_{i}(x_{*})\|$. Also we regard $x_{i}(t)$ as a row vector in
$\mathbb{R}^{1\times d}$, and define the variable
$\mathbf{x}(t)\in\mathbb{R}^{n\times d}$ and
$\bar{\mathbf{x}}(t)\in\mathbb{R}^{n\times d}$ by

$\mathbf{x}(t)=\left(x_{1}\left(t\right)^{T},\cdots,x_{m}\left(t\right)^{T}\right)^{T}\quad\textrm{and}\quad\bar{\mathbf{x}}(t)=\left(\bar{x}\left(t\right)^{T},\cdots,\bar{x}\left(t\right)^{T}\right)^{T},$
(1.5)

where $\bar{x}(t)=\frac{1}{n}\sum_{k=1}^{n}x_{k}(t)$. We let

$A_{t}=\|\bar{x}(t)-x_{*}\|\quad\textrm{and}\quad
B_{t}=\frac{1}{\sqrt{n}}\|\mathbf{\bar{x}}(t)-\mathbf{x}(t)\|=\Big{(}\frac{1}{n}\sum_{k=1}^{n}\|\bar{x}(t)-x_{k}(t)\|^{2}\Big{)}^{1/2}.$

We show that the DPPA converges exponentially to an $O(\eta)$-neighborhood of
the optimizer in the following result.

######  Theorem 1.2.

Suppose that the assumptions 1-3 hold and Then we have

$A_{t}\leq\frac{A_{0}}{(1+\eta\alpha)^{t}}+\frac{t\eta
LB_{0}\rho_{W}}{(1+\eta\alpha)}\max\Big{\\{}\rho_{W},\frac{1}{1+\eta\alpha}\Big{\\}}^{t-1}+\frac{\eta
L(2RL+D)}{\alpha(1-\rho_{W})}$ (1.6)

and

$B_{t}\leq(\rho_{W})^{t}B_{0}+\frac{\eta}{(1-\rho_{W})}(2RL+D).$

Upon knowing that the DGD is stable, the linear convergence result of the DGD
is proved when the stepsize $\eta$ satisfies $\eta\leq\frac{2}{\alpha+L}$. We
refer to [29, 7] for the detail. As for the DPPA, we note that there is no
restriction on the stepsize $\eta>0$ for the linear convergence as obtained in
Theorem 1.2. Table 2 compares the condition on the stepsize of the DGD and the
DPPA for the convergence of the algorithms.

| Condition on the costs | DGD | Estimate | Paper
---|---|---|---|---
DGD |  Each $f_{j}$ is $L$-smooth $f$ is $\alpha$-strongly convex | $\eta\leq\frac{2}{\alpha+L}$ | $O(e^{-ct})+O\Big{(}\frac{\eta}{1-\rho_{W}}\Big{)}$ | [29, 7]
DPPA |  Each $f_{j}$ is $L$-smooth $f$ is $\alpha$-strongly convex | $\eta\in(0,\infty)$ | $O(e^{-ct})+O\Big{(}\frac{\eta}{1-\rho_{W}}\Big{)}$ | This work
Table 2. This table compares the condition on the stepsize for the stability
of the DGD and the DPPA.

The rest of this paper is organized as follows. In Section 2, we establish two
sequential inequalities of $A_{t}$ and $B_{t}$. Section 3 is devoted to prove
the uniform boundedness result of Theorem 1.1. In Section 4, we prove the
convergence result of Theorem 1.2. Numerical results are presented in Section
5.

##  2\. Sequential estimates

In this section, we dervie two sequential inequalities of $A_{t}$ and $B_{t}$,
which are main ingredients for the stability and the convergence analysis of
the DPPA.

######  Proposition 2.1.

Assume that $f$ is $\alpha$-strongly convex. Then, for any stepsize $\eta>0$,
the sequence $\\{(A_{t},B_{t})\\}_{t\geq 0}$ satisfies the follwoing
inequality

$(1+\eta\alpha)A_{t+1}\leq A_{t}+\eta LB_{t+1}$ (2.1)

for all $t\geq 0$.

###### Proof.

From the minimality of (1.2), it follows that

$\nabla
f_{k}(x_{k}(t+1))+\frac{1}{\eta}\Big{(}x_{k}(t+1)-\sum_{j=1}^{n}w_{kj}x_{j}(t)\Big{)}=0.$
(2.2)

We reformulate this as

$x_{k}(t+1)+\eta\nabla f_{k}(x_{k}(t+1))=\sum_{j=1}^{n}w_{kj}x_{j}(t).$ (2.3)

By averaging this for $1\leq k\leq n$, we get

$\bar{x}(t+1)+\frac{\eta}{n}\sum_{k=1}^{n}\nabla
f_{k}(x_{k}(t+1))=\bar{x}(t).$ (2.4)

Using this and the fact that $\nabla f(x_{*})=0$, we find

$\begin{split}&\bar{x}(t+1)-x_{*}+\eta(\nabla f(\bar{x}(t+1))-\nabla
f(x_{*}))\\\ &=\bar{x}(t)-x_{*}+\eta\Big{(}\sum_{k=1}^{n}\nabla
f_{k}(\bar{x}(t+1))-\nabla f_{k}(x_{k}(t+1))\Big{)}.\end{split}$ (2.5)

Since $f$ is assumed to be $\alpha$-strongly convex, we have

$\|\nabla f(x)-\nabla f(y)\|^{2}\geq\alpha^{2}\|x-y\|^{2}$

and

$\langle x-y,~{}\nabla f(x)-\nabla f(y)\rangle\geq\alpha\|x-y\|^{2}.$

Combining these estimates, we get

$\begin{split}&\Big{\|}\bar{x}(t+1)-x_{*}+\eta\Big{(}\nabla
f(\bar{x}(t+1))-\nabla f(x_{*})\Big{)}\Big{\|}^{2}\\\
&=\|\bar{x}(t+1)-x_{*}\|^{2}+2\eta\langle\bar{x}(t+1)-x_{*},~{}\nabla
f(\bar{x}(t+1))-\nabla f(x_{*})\rangle\\\ &\qquad+\eta^{2}\|\nabla
f(\bar{x}(t+1))-\nabla f(x_{*})\|^{2}\\\
&\geq(1+2\eta\alpha+\alpha^{2})\|\bar{x}(t+1)-x_{*}\|^{2}.\end{split}$

Using this estimate in (2.5) and applying the triangle inequality, we get

$\begin{split}&(1+\eta\alpha)\|\bar{x}(t+1)-x_{*}\|\\\
&\leq\|\bar{x}(t)-x_{*}\|+\frac{\eta}{n}\Big{\|}\sum_{k=1}^{n}\nabla
f_{k}(\bar{x}(t+1))-\nabla f_{k}(x_{k}(t+1))\Big{\|}\\\
&\leq\|\bar{x}(t)-x_{*}\|+\frac{\eta
L}{n}\sum_{k=1}^{n}\|\bar{x}(t+1)-x_{k}(t+1)\|\\\
&\leq\|\bar{x}(t)-x_{*}\|+\frac{\eta
L}{\sqrt{n}}\|\bar{\mathbf{x}}(t+1)-\mathbf{x}(t+1)\|,\end{split}$

where we used the Cauchy-Schwartz inequality in the last inequality. This
gives the desired inequality. ∎

Next we will derive a bound of $B_{t+1}$ in terms of $A_{t}$ and $B_{t}$. For
this we will use the following result (see [19, Lemma 1]).

######  Lemma 2.2.

Suppose Assumptions 1 and 2 hold, and let $\rho_{W}$ be the spectral norm of
the matrix $W-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}.$ Then we have $\rho_{W}<1$
and

$\sum_{i=1}^{n}\Big{\|}\sum_{j=1}^{n}w_{ij}(x_{j}-\bar{x})\Big{\|}^{2}\leq(\rho_{W})^{2}\sum_{i=1}^{n}\|x_{i}-\bar{x}\|^{2},$

where $\bar{x}=\frac{1}{n}\sum_{k=1}^{n}x_{k}$ for any
$x_{i}\in\mathbb{R}^{d\times 1}$ and $1\leq i\leq n$.

In the following result, we find an estimate of $B_{t+1}$ in terms of $B_{t}$
and $A_{t+1}$.

######  Proposition 2.3.

Suppose that each $f_{j}$ is $L$-smooth. Then the sequence
$\\{(A_{t},B_{t})\\}_{t\geq 0}$ satisfies the following inequality

$B_{t+1}\leq\rho_{W}B_{t}+\eta LB_{t+1}+\eta LA_{t+1}+\eta D$ (2.6)

for all $t\geq 0$.

###### Proof.

We may write (2.3) and (2.4) in the following way

$\mathbf{x}(t+1)+\eta\nabla F(\mathbf{x}(t+1))=W\mathbf{x}(t)$

and

$\bar{\mathbf{x}}(t+1)+\eta\overline{\nabla
F}(\mathbf{x}(t+1))=\bar{\mathbf{x}}(t),$

where we have let

$\overline{\nabla F}(\mathbf{x}(t+1))=\Big{(}\nabla
f_{1}(x_{1}(t+1)),\cdots,\nabla f_{n}(x_{n}(t+1))\Big{)}^{T}.$

Combining the above equalities, we find

$\mathbf{x}(t+1)-\bar{\mathbf{x}}(t+1)=W(\mathbf{x}(t)-\bar{\mathbf{x}}(t))-\eta\Big{(}\nabla
F(\mathbf{x}(t+1))-\overline{\nabla F}(\mathbf{x}(t+1))\Big{)}.$

By applying the triangle inequality and Lemma 2.2, we deduce

$\begin{split}&\|\mathbf{x}(t+1)-\bar{\mathbf{x}}(t+1)\|\\\
&\leq\rho_{W}\|\mathbf{x}(t)-\bar{\mathbf{x}}(t)\|+\eta\|\nabla
F(\mathbf{x}(t+1))-\overline{\nabla F}(\mathbf{x}(t+1))\|.\end{split}$ (2.7)

Using the fact that the spectral radius of the matrix
$\Big{(}I_{n}-\frac{1}{n}1_{n}^{T}1_{n}\Big{)}$ is one, we obtain

$\begin{split}&\|\nabla F(\mathbf{x}(t+1))-\overline{\nabla
F}(\mathbf{x}(t+1))\|\\\ &\leq\|\nabla F(\mathbf{x}(t+1))\|\\\ &\leq\|\nabla
F(\mathbf{x}(t+1))-\nabla F(\bar{\mathbf{x}}(t+1))\|+\|\nabla
F(\bar{\mathbf{x}}(t+1))-\nabla F(x_{*})\|+\|\nabla F(x_{*})\|\\\ &\leq
L\|\mathbf{x}(t+1)-\bar{\mathbf{x}}(t+1)\|+\sqrt{n}L\|\bar{{x}}(t+1)-{x}_{*}\|+\sqrt{n}D.\end{split}$

Inserting this into (2.7) we obtain

$\begin{split}&\|\mathbf{x}(t+1)-\bar{\mathbf{x}}(t+1)\|\\\
&\leq\rho_{W}\|\mathbf{x}(t)-\bar{\mathbf{x}}(t)\|+\eta
L\|\mathbf{x}(t+1)-\bar{\mathbf{x}}(t+1)\|\\\
&\quad+\sqrt{n}L\eta\|\bar{x}(t+1)-x_{*}\|+\sqrt{n}D\eta,\end{split}$

which is the desired inequality. ∎

##  3\. Boundedness of the sequence

We prove the uniform boundedness result of Theorem 1.1 under the two
conditions separately in the below.

###### Proof of Theorem 1.1.

Assume the first condition of Theorem 1.1. We claim the following inequality

$\sum_{k=1}^{n}\|x_{k}(t+1)-x_{k}^{*}\|\leq\sum_{k=1}^{n}\|x_{k}(t)-x_{k}^{*}\|\quad\forall~{}t\geq
0.$ (3.1)

The optimizer $(x_{1}^{*},\cdots,x_{n}^{*})$ of (2.7) satisfies

$\nabla
f_{k}(x_{k}^{*})+\frac{1}{\eta}\Big{(}x_{k}^{*}-\sum_{j=1}^{n}w_{kj}x_{j}^{*}\Big{)}=0.$

Combining this with (2.2) gives

$\eta(\nabla f_{k}(x_{k}(t+1))-\nabla
f_{k}(x_{k}^{*}))+(x_{k}(t+1)-x_{k}^{*})=\sum_{j=1}^{n}w_{kj}(x_{j}(t)-x_{j}^{*}).$

From this we have

$\|x_{k}(t+1)-x_{k}^{*}\|\leq\Big{\|}\sum_{j=1}^{n}w_{kj}(x_{j}(t)-x_{j}^{*})\Big{\|}\leq\sum_{j=1}^{n}w_{kj}\|x_{j}(t)-x_{j}^{*}\|.$

Summing up this for $1\leq k\leq n$, we get

$\sum_{k=1}^{n}\|x_{k}(t+1)-x_{k}^{*}\|\leq\sum_{k=1}^{n}\sum_{j=1}^{n}w_{kj}\|x_{j}(t)-x_{j}^{*}\|=\sum_{j=1}^{n}\|x_{j}(t)-x_{j}^{*}\|,$

which proves the inequality (3.1). It gives the following bound

$\sum_{k=1}^{n}\|x_{k}(t)-x_{k}^{*}\|\leq\sum_{k=1}^{n}\sum_{j=1}^{n}w_{kj}\|x_{j}(t)-x_{j}^{*}\|=\sum_{j=1}^{n}\|x_{j}(0)-x_{j}^{*}\|\quad\forall~{}t\geq
0.$

Hence $A_{t}$ and $B_{t}$ are uniformly bounded.

Next we assume the second condition of Theorem 1.1. Then we claim that the
sequence $A_{t}$ and $B_{t}$ satisfies

$A_{t}\leq R\quad\textrm{and}\quad B_{t}\leq\frac{\alpha}{L}R,$

where $R>0$ is defined by

$R=\max\Bigg{\\{}A_{0},~{}\frac{L}{\alpha}B_{0},~{}\frac{L}{\alpha}B_{1},~{}\frac{\eta
D}{\frac{\alpha}{L}\Big{(}1-\eta
L-\frac{\eta^{2}L^{2}}{(1+\eta\alpha)}\Big{)}-\frac{\alpha}{L}\rho_{W}-\frac{\eta
L}{(1+\eta\alpha)}}\Bigg{\\}}.$ (3.2)

We argue by induction to prove this claim. First, we note that $A_{0}\leq R$,
$B_{0}\leq R$ and $B_{1}\leq R$ by the definition of $R$. Next we assume that

$A_{t}\leq R\quad\textrm{and}\quad B_{t+1}\leq cR$ (3.3)

for some $t\geq 0$ and $c=\frac{\alpha}{L}$. Then, we use these bounds in
(2.1) to find

$(1+\eta\alpha)A_{t+1}\leq(R+\eta LcR)=(1+\eta\alpha)R,$

which gives $A_{t+1}\leq R$. Next we recall the estimates (2.1) and (2.6) with
$t$ replaced by $(t+1)$ as

$B_{t+2}\leq\rho_{W}B_{t+1}+\eta LB_{t+2}+\eta LA_{t+2}+\eta D$

and

$\begin{split}(1+\eta\alpha)A_{t+2}&\leq A_{t+1}+\eta LB_{t+2}\\\ &\leq R+\eta
LB_{t+1},\end{split}$

where we used $A_{t+1}\leq R$ in the last inequality. Combining these
estimates with (3.3) gives

$\Big{(}1-\eta
L-\frac{\eta^{2}L^{2}}{(1+\eta\alpha)}\Big{)}B_{t+2}\leq(\rho_{W})(cR)+\frac{\eta
L}{1+\eta\alpha}R+\eta D.$

This gives $B_{t+2}\leq cR$ provided that

$\rho_{W}(cR)+\frac{\eta L}{(1+\eta\alpha)}R+\eta D\leq Rc\Big{(}1-\eta
L-\frac{\eta^{2}L^{2}}{(1+\eta\alpha)}\Big{)}.$

This holds true for $R>0$ defined in (3.2) and $\eta>0$ satisfying

$\frac{\alpha}{L}\Big{(}(1-\rho_{W})-\eta
L-\frac{\eta^{2}L^{2}}{(1+\eta\alpha)}\Big{)}>\frac{\eta L}{(1+\eta\alpha)},$

which is same with

$\eta^{2}(\alpha^{2}+\alpha
L)+\eta\Big{(}\alpha+L-\frac{(1-\rho_{W})\alpha^{2}}{L}\Big{)}<\frac{(1-\rho_{W})\alpha}{L}.$

The proof is done.

∎

##  4\. Convergence result

In this section, we prove the main convergence result of the decentralized
proximal point method.

###### Proof of Theorem 1.2.

By Proposition 2.1 and Proposition 2.3 we have the following inequalities

$(1+\eta\alpha)A_{t+1}\leq A_{t}+\eta LB_{t+1}$ (4.1)

and

$B_{t+1}\leq\rho_{W}B_{t}+\eta LB_{t+1}+\eta LA_{t+1}+\eta D.$

Using the assumption that $A_{t+1}\leq R$ and $B_{t+1}\leq R$, we have

$B_{t+1}\leq(\rho_{W})B_{t}+\eta(2RL+D)$

for all $t\geq 0$. Using this iteratively, we get

$\begin{split}B_{t}&\leq(\rho_{W})^{t}B_{0}+\eta(2RL+D)\Big{[}1+\rho_{W}+\cdots+(\rho_{W})^{t-1}\Big{]}\\\
&\leq(\rho_{W})^{t}B_{0}+\frac{\eta}{(1-\rho_{W})}(2RL+D).\end{split}$

Putting this inequality into (4.1) leads to

$\begin{split}&(1+\eta\alpha)A_{t+1}\\\ &\leq A_{t}+\eta
L\Big{[}(\rho_{W})^{t+1}B_{0}+\frac{\eta}{(1-\rho_{W})}(2RL+D)\Big{]}\\\
&=A_{t}+\eta
LB_{0}(\rho_{W})^{t+1}+\frac{\eta^{2}L}{(1-\rho_{W})}(2RL+D).\end{split}$
(4.2)

To analyize this sequential estimate, we consider two positive sequences
$\\{z_{t}\\}_{t\geq 0}$ and $\\{y_{t}\\}_{t\geq 0}$ satisfying

$(1+\eta\alpha)z_{t+1}=z_{t}+\eta LB_{0}(\rho_{W})^{t+1}$

and

$(1+\eta\alpha)y_{t+1}=y_{t}+\frac{\eta^{2}L}{(1-\rho_{W})}(2RL+D),$

where $z_{0}=A_{0}$ and $y_{0}=0$. It then follows from (4.2) that $A_{t}\leq
z_{t}+y_{t}$ for all $t\geq 0$.

We estimate $z_{t}$ as follows:

$\begin{split}z_{t}&=\frac{z_{0}}{(1-\eta\alpha)^{t}}+\frac{\eta
LB_{0}\rho_{W}}{(1+\eta\alpha)}\Big{[}\sum_{j=0}^{t-1}\frac{(\rho_{W})^{j}}{(1+\eta\alpha)^{t-1-j}}\Big{]}\\\
&\leq\frac{z_{0}}{(1-\eta\alpha)^{t}}+\frac{t\eta
LB_{0}\rho_{W}}{(1+\eta\alpha)}\max\Big{\\{}\rho_{W},\frac{1}{1+\eta\alpha}\Big{\\}}^{t-1}.\end{split}$

Next we estimate $y_{t}$ as

$\begin{split}y_{t}&=\frac{1}{(1+\eta\alpha)^{t}}y_{0}+\frac{\eta^{2}L}{(1+\eta\alpha)(1-\rho_{W})}(2RL+D)\sum_{j=0}^{t-1}\frac{1}{(1+\eta\alpha)^{j-1}}\\\
&\leq\frac{1}{(1+\eta\alpha)^{t}}y_{0}+\frac{\eta^{2}L}{(1+\eta\alpha)(1-\rho_{W})}(2RL+D)\sum_{j=0}^{\infty}\frac{1}{(1+\eta\alpha)^{j}}\\\
&=\frac{y_{0}}{(1+\eta\alpha)^{t}}+\frac{\eta
L(2RL+D)}{\alpha(1-\rho_{W})}.\end{split}$

Combining the above estimates gives

$\begin{split}A_{t}&\leq z_{t}+y_{t}\\\
&\leq\frac{A_{0}}{(1-\eta\alpha)^{t}}+\frac{t\eta
LB_{0}\rho_{W}}{(1+\eta\alpha)}\max\Big{\\{}\rho_{W},\frac{1}{1+\eta\alpha}\Big{\\}}^{t-1}+\frac{\eta
L(2RL+D)}{\alpha(1-\rho_{W})}.\end{split}$

This corresponds to the inequality (1.4) in the condition (2). Thus, it holds
that $A_{t+1}\leq R$ and $B_{t+2}\leq cR$. This completes the induction, and
so the claim is proved. Hence the uniform boundedness is proved. ∎

##  5\. Numerical tests

This section gives numerical experimental results of the DPPA. We consider the
cost function

$f(x)=\frac{1}{n}\sum_{k=1}^{n}\|A_{k}x-y_{k}\|^{2},$

where $n$ is the number of agents and for each $1\leq i\leq n$ and $A_{i}$ be
an $m\times d$ matrix whose element is chosen randomly following the normal
distribution $N(0,1)$. Also we set $y_{i}\in\mathbb{R}^{m}$ whose each element
is generated from the normal distribution $N(0,1)$. We choose $n=20$, $m=5$
and $d=10$.

For the communication matrix $W$, we link each two agents with probability
0.4, and define the weights $w_{ij}$ by

$w_{ij}=\left\\{\begin{array}[]{ll}1/\max\\{\textrm{deg}(i),\textrm{deg}(j)\\}&~{}\textrm{if}~{}i\in
N_{i},\\\ 1-\sum_{j\in N_{i}}w_{ij}&~{}\textrm{if}~{}i=j,\\\
0&~{}otherwise.\end{array}\right.$

We define the following value

$\eta_{c}=\frac{1+\lambda_{n}(W)}{L},$

where $\lambda_{n}(W)$ is the smallest eigenvalue of $W$ and $L$ is the
constant for the smoothness propery of $f_{j}$ for $1\leq j\leq n$. In our
experiment, the constants are computed as $L\simeq 29.7312$ and
$\lambda_{n}(W)\simeq-0.4009$. Therefore we find $\eta_{c}\simeq 0.0202$.

We perform the numerical simulation for the DPPA and the DGD with stepsize
chosen as $\eta=\eta_{c}+0.005$ and $\eta=\eta_{c}$. Figure 1 indicates the
graphs of the error $\log(\sum_{k=1}^{n}\|x_{k}(t)-x_{*}\|/n)$ with respect to
$t\geq 0$. The result shows that both the DPPA and the DGD are stable for
$\eta=\eta_{c}$ as the theoretical results of Table 1 guarantee. On the other
hand, if we choose $\eta=\eta_{c}+0.005$, then the DGD becomes unstable while
the DPPA is still stable. This supports the result of Theorem 1.1.

Figure 1. The graphs of the error $\log(\sum_{k=1}^{n}\|x_{k}(t)-x_{*}\|/n)$
for the DPPA and the DGD for $t\geq 0$ with stepsize $\eta=\eta_{c}+0.005$ and
$\eta=\eta_{c}$.

Next we perform the numerical test for the DPPA with various stepsize
$\eta\in\\{0.001,0.01,0.1,1,2\\}$. We measure the error
$\log(\sum_{k=1}^{n}\|x_{k}(t)-x_{*}\|/n)$ and the graphs are given in Figure
2. The result shows that the error decreases exponentially up to an $O(\eta)$
value, which supports the convergence result of Theorem 1.2.

Figure 2. The graphs of the error $\log(\sum_{k=1}^{n}\|x_{k}(t)-x_{*}\|/n)$
for the DPPA with respect to $t\geq 0$ and stepsize
$\eta\in\\{0.001,0.01,0.1,1,2\\}$.

## Acknowledgments

The author was supported by the National Research Foundation of Korea (NRF)
grants funded by the Korea government No. NRF-2016R1A5A1008055 and No.
NRF-2021R1F1A1059671.

## References

  * [1] J. A. Bazerque, G.B. Giannakis, Distributed spectrum sensing for cognitive radio networks by exploiting sparsity. IEEE Trans. Signal Process. 58(3), 1847–1862 (2010)
  * [2] A.S. Berahas, R. Bollapragada, N.S. Keskar, E. Wei, Balancing communication and computation in distributed optimization. IEEE Trans. Autom. Control 64(8), 3141–3155 (2018)
  * [3] L. Bottou, F. E. Curtis, and J. Nocedal, Optimization methods for large-scale machine learning, SIAM Review, vol. 60, no. 2, pp. 223–-311, 2018.
  * [4] F. Bullo, J. Cortes, S. Martinez, S. Distributed Control of Robotic Networks: A Mathematical Approach to Motion Coordination Algorithms, vol. 27. Princeton University Press, Princeton (2009)
  * [5] Y. Cao, W. Yu, W. Ren, and G. Chen, An overview of recent progress in the study of distributed multi-agent coordination, IEEE Trans Ind. Informat., 9 (2013), pp. 427–438
  * [6] W. Choi, D. Kim, S. Yun, Convergence results of a nested decentralized gradient method for non-strongly convex problems. J. Optim. Theory Appl. 195 (2022), no. 1, 172–204.
  * [7] W. Choi, J. Kim, On the convergence of decentralized gradient descent with diminishing stepsize, revisited, arXiv:2203.09079.
  * [8] A. I. Chen and A. Ozdaglar, A fast distributed proximal-gradient method, in Communication, Control, and Computing (Allerton), 2012 50th Annual Allerton Conference on. IEEE, 2012, pp. 601–608.
  * [9] G. B. Giannakis, V. Kekatos, N. Gatsis, S.-J. Kim, H. Zhu, and B. Wollenberg, Mon-itoring and optimization for power grids: A signal processing perspective, IEEE Signal Processing Mag., 30 (2013), pp. 107–128,
  * [10] V. Kekatos and G. B. Giannakis, Distributed robust power system state estimation, IEEE Trans. Power Syst., 28 (2013), pp. 1617–1626,
  * [11] P. A. Forero, A. Cano, and G. B. Giannakis, Consensus-based distributed support vector machines, Journal of Machine Learning Research, vol. 11, pp. 1663–-1707, 2010.
  * [12] D. Jakovetic, J. Xavier, and J. M. F. Moura, Fast Distributed Gradient Methods, IEEE Transactions on Automatic Control, vol. 59, no. 5, pp. 1131–1146, 2014.
  * [13] X. Li, G. Feng, L. Xie, Distributed proximal algorithms for multiagent optimization with coupled inequality constraints. IEEE Trans. Automat. Control 66 (2021), no. 3, 1223–1230.
  * [14] X. Li, G. Feng, L. Xie, Distributed Proximal Point Algorithm for Constrained Optimization over Unbalanced Graphs. In 2019 IEEE 15th International Conference on Control and Automation (ICCA) (pp. 824-829).
  * [15] K. Margellos, A. Falsone, S. Garatti, and M. Prandini, Proximal minimization based distributed convex optimization,” in 2016 American Control Conference (ACC), July 2016, pp. 2466–2471.
  * [16] K. Margellos, A. Falsone, S. Garatti, and M. Prandini, Distributed constrained optimization and consensus in uncertain networks via proximal minimization, IEEE Trans. Autom. Control, vol. 63, no. 5, pp. 1372–1387, May 2018.
  * [17] A. Nedić and A. Ozdaglar, Distributed subgradient methods for multi-agent optimization, IEEE Transactions on Automatic Control, vol. 54, no. 1, pp. 48, 2009.
  * [18] A. Nedić, A. Olshevsky, and W. Shi, Achieving geometric convergence for distributed optimization over time-varying graphs, SIAM Journal on Optimization, vol. 27, no. 4, pp. 2597–2633, 2017.
  * [19] S. Pu and A. Nedić, Distributed stochastic gradient tracking methods, Math. Program, pp. 1–49, 2018
  * [20] G. Qu and N. Li, Harnessing smoothness to accelerate distributed optimization, IEEE Transactions on Control of Network Systems, vol. 5, no. 3, pp. 1245–1260, 2018.
  * [21] H. Raja and W. U. Bajwa, Cloud K-SVD: A collaborative dictionary learning algorithm for big, distributed data, IEEE Transactions on Signal Processing, vol. 64, no. 1, pp. 173–-188, Jan. 2016.
  * [22] W. Ren, Consensus Based Formation Control Strategies for Multi-Vehicle Systems, in Proceedings of the Amer-ican Control Conference, 2006, pp. 4237–-4242.
  * [23] A. H. Sayed, Diffusion adaptation over networks. Academic Press Library in Signal Processing, 2013, vol. 3.
  * [24] I.D. Schizas, G.B. Giannakis, S.I. Roumeliotis, A. Ribeiro, Consensus in ad hoc WSNs with noisy links—part II: distributed estimation and smoothing of random signals. IEEE Trans. Signal Process. 56(4), 1650–1666 (2008)
  * [25] W. Shi, Q. Ling, G. Wu, W. Yin, Extra: an exact first-order algorithm for decentralized consensus optimization. SIAM J. Optim. 25(2), 944–966 (2015)
  * [26] W. Shi, Q. Ling, K. Yuan, G. Wu, W. Yin, On the linear convergence of the ADMM in decentralized consensus optimization. IEEE Trans. Signal Process. 62(7), 1750–1761 (2014)
  * [27] K. I. Tsianos, S. Lawlor, and M. G. Rabbat, Consensus-based distributed optimization: Practical issues and applications in large-scale machine learning, in Proceedings of the IEEE Allerton Conference on Communication, Control, and Computing, IEEE, New York, 2012, pp. 1543–1550,
  * [28] R. Xin and U. A. Khan, A linear algorithm for optimization over directed graphs with geometric convergence, IEEE Control Systems Letters, vol. 2, no. 3, pp. 315–320, 2018.
  * [29] K. Yuan, Q. Ling, and W. Yin, On the convergence of decentralized gradient descent, SIAM Journal on Optimization, vol. 26, no. 3, pp. 1835–1854, Sep. 2016.
  * [30] T. Yang, X. Yi, J. Wu, Y. Yuan, D. Wu, Z. Meng, Y. Hong, H. Wang, Z. Lin, and K. H. Johansson, A survey of distributed optimization, Annual Reviews in Control, vol. 47, pp. 278–305, 2019.


  *

# Ising model on the aperiodic Smith hat

Yutaka Okabe1, Komajiro Niizeki2111Professor emeritus and Yoshiaki Araki3
1Department of Physics, Tokyo Metropolitan University, Hachioji, Tokyo
192-0397, Japan
2Graduate School of Science, Tohoku University, Aoba-ku, Sendai 980-8578,
Japan
3Japan Tessellation Design Association, Musashino-shi, Tokyo 180-0003, Japan
<EMAIL_ADDRESS>
###### Abstract

Smith et al discovered an aperiodic monotile of 13-sided shape in 2023. It is
called the ‘Smith hat’ and consists of 8 kites. We deal with the statistical
physics of the lattice of the kites, which we call the ‘Smith-kite lattice’.
We studied the Ising model on the aperiodic Smith-kite lattice and the dual
Smith-kite lattice using Monte Carlo simulations. We combined the Swendsen-
Wang multi-cluster algorithm and the replica exchange method. We simulated
systems up to the total spin number $939201$. Using the finite-size scaling
analysis, we estimated the critical temperature on the Smith-kite lattice as
$T_{c}/J=2.405\pm 0.0005$ and that of the dual Smith-kite lattice as
$T^{*}_{c}/J=2.143\pm 0.0005$. Moreover, we confirmed the duality relation
between the critical temperatures on the dual pair of aperiodic lattices,
$\sinh(2J/T_{c})\sinh(2J/T^{*}_{c})=1.000\pm 0.001$. We also checked the
duality relation for the nearest-neighbor correlation at the critical
temperature, essentially the energy,
$\epsilon(T_{c})/\coth(2J/T_{c})+\epsilon(T^{*}_{c})/\coth(2J/T^{*}_{c})=1.000\pm
0.001$.

  * Keywords: Smith hat, aperiodic lattice, Ising model, Monte Carlo simulation, duality relation

##  1 Introduction

It is possible to fill a two-dimensional plane with a single shape using
regular polygons: equilateral triangles, squares, and regular hexagons. It has
long been a matter of interest to determine whether there are other shapes.
Penrose [2] showed that a plane can be filled with two shapes, a fat rhombus
and a thin rhombus, where the golden ratio, an irrational number, is involved.
The tilings obtained from these shapes do not have the usual periodic
structure but creates a quasiperiodic structure. It exhibits the 5-fold
(10-fold) symmetry. The quasicrystal [3] was studied by experimental studies,
and Shechtman et al discovered the 5-fold (10-fold) symmetry [4]. The Nobel
Prize in Chemistry 2011 was awarded to Shechtman “for the discovery of
quasicrystals”.

Smith, Myers, Kaplan and Goodman-Strauss discovered an aperiodic monotile in
2023 [5]. The shape is 13-sided and called the ‘hat’. It is able to fill an
infinite plane without overlaps or gaps in a pattern that not only never
repeats but also never can be repeated. Various studies have been conducted
following this discovery [6, 7, 8, 9]. In [7] and [8], the vertex set of the
monotile is shown to be a two-dimensional quasiperiodic lattice (2DQL) which
is obtained from a six-dimensional or a four-dimensional periodic lattice by
the cut-and-projection method. 2DQLs are divided into two groups according to
the nature of the projection window used in the cut-and-projection method. The
window of a 2DQL belonging to one of the two groups is a polygon but the one
belonging to the other group is a domain with a fractal boundary. A
representative 2DQL of the first group is Penrose tiling P1, whose major tile
is a regular pentagon. The monotile is shown, however, to belong to the second
group [7, 8].

The material science of quasilattices has been studied experimentally and
theoretically in various ways. Recent studies have focused on
superconductivity [10, 11], electric structure [12] and magnetic properies
[13]. There has been much interest in the phase transitions of spin systems on
quasicrystals compared to those of ordered and random systems. Okabe and
Niizeki [14] studied spin systems on a Penrose lattice mainly by means of
Monte Carlo simulation. They showed that the Ising model on the Penrose
lattice belongs to the same universality class as that on the two-dimensional
(2D) regular lattice. In addition, they studied the dual lattice [15] and
discussed the duality relation. A high-precision study was performed recently
up to $N=20,633,239$ using a GPU-based calculation [16]; A significant
improvement in accuracy was obtained.

In contrast, no one has investigated up to present spin systems on a 2DQL
belonging to one of the second group in the classification above. In this
paper, we study the Ising models on the Smith-kite lattice and its dual
lattice by means of the Monte Carlo simulation. We use the Swendsen-Wang (SW)
multi-cluster spin-flip algorithm [17] combined with the replica exchange
method of parallel tempering [18]. We organize the rest of the paper as
follows. In Sec. 2, we explain the model and numerical method. We give the
results in Sec. 3. We report the estimate of the critical temperature, and
discuss the duality relation. The final section is devoted to the summary and
discussions.

##  2 Model and Numerical Method

###  2.1 Lattice

Smith, Myers, Kaplan and Goodman-Strauss announced the discovery of an
aperiodic monotile in March 2023 [5]. This brand new 13-sided shape is called
‘the Smith hat’. We here study statistical physics of spin systems on the
Smith-hat tiling. As for the coordinate data, we use the code released by
Cheritat [19], which is a web application written in JavaScript. The
4-dimensional coordinate data of the hat vertices are given there, and by
projecting them into two dimensions, 2-dimensional coordinate data can be
obtained. Cheritat provides four sizes of Smith-hat tiling data. The number of
hats is 357, 2490, 17077 and 117051 for “small”, ”medium”, “big” and “bigger”
hat tiling, respectively. The Smith-hat tiling of the “small” size is
illustrated in figure 1. It is noteworthy to mention a family of aperiodic
monotiles [5]. The notation Tile($a,b$), where $a$ and $b$ are real numbers,
is used to specify a shape of hat tile. Tile($1,\sqrt{3}$), which is the most
popular, has 8 kites in the hat, whereas Tile($\sqrt{3},1$) has 10 kites in
the hat. Although we mainly treat Tile($1,\sqrt{3}$) of a 13-sided monotile
with 8 kites in this paper, we show both Tile($1,\sqrt{3}$) and
Tile($\sqrt{3},1$) in figure 1. In the site of Cheritat [19], the default
value of $\alpha_{1}$ is provided as 60 for Tile($1,\sqrt{3}$). When we treat
Tile($\sqrt{3},1$), we set $\alpha_{1}$ as 30.

Figure 1:  Smith-hat tiling of the “small” size ($N_{\rm hat}=357$). We use
the code released by Cheritat [19]. The left illusration is Tile($1,\sqrt{3}$)
of 8 kites, whereas the right illusration is Tile($\sqrt{3},1$) of 10 kites.

We consider the lattice based on the kites. A hat consists of 8 or 10 kites.
The vertices of the kites become the lattice points. There are extra 2 (4)
lattice points in an 8-kite (10-kite) hat. We need the information on the
connection of kite vertices; we make a table of connecting sites by checking
the vertices shared by neighboring hats. One vertex is shared by two, three or
four hats. In this way, we form a lattice made by the vertices of the kites,
and we call it the original Smith-kite lattice. We also consider the dual
lattice, which we call the dual Smith-kite lattice. The lattice points of the
dual lattice are set in a kite. We make the connection between the dual
lattice points by the two adjacent kites. These two lattices are dual, and we
will examine the duality relations for dual lattices discussed by Kramers and
Wannier [20]. The original and dual Smith-kite lattices are illustrated in
figure 2. The lattice points of the original lattice are denoted by green
circles, whereas those of the dual lattice are denoted by empty red circles.
The edges of the dual lattice are denoted by red dotted lines. We put spins on
the lattice points when we study the Ising model.

Figure 2:  Smith-kite lattices. The left panel is Tile($1,\sqrt{3}$) of 8
kites, whereas the right panel is Tile($\sqrt{3},1$) of 10 kites. The edges of
the Smith hat are denoted by blue lines. One hat is denoted by a shaded area.
A hat consists of 8 or 10 kites, and the edges of the kites are denoted by
thin green lines, which form the original Smith-kite lattice. The lattice
points of the original Smith-kite lattice are denoted by green circles,
whereas those of the dual Smith-kite lattice are denoted by empty red circles.
We can see that there are extra 2 (4) lattice points in one 8-kite (10-kite)
hat. The edges of the dual lattice are denoted by red dotted lines. In the
Ising model, spins are put on the lattice points.

The sizes of the lattices are tabulated in table 1. This is the data for
Tile($1,\sqrt{3}$) of 8 kites. The second column gives the number of hats for
four sizes. The number of center of kites is 8 times the number of hats,
whereas the number of vertex of kites is slightly larger because of the
boundary. When we discuss the averaged values for the spins in the simulation,
it is better to take an average over only inner sites, in order to reduce
surface effects. Thus, we also list the number of inner spins.

Table 1:  Size of Smith-hat tiling provided by Cheritat [19]. This is the data for Tile($1,\sqrt{3}$) of 8 kites. The original Smith-kite lattice is formed by the vertex of kites, whereas the dual Smith-kite lattice is formed by the center of kites. The number of the inner sites will be used when taking an average in the simulation, and it will be denoted by $N$.  | hats | vertex of kites | center of kites
---|---|---|---
|  |  total |  inner |  total |  inner
“small” |  357 |  3011 |  2807 |  2856 |  2655
“medium” |  2490 |  20388 |  19794 |  19920 |  19372
“big” |  17077 |  137686 |  136266 |  136616 |  135216
“bigger” |  117051 |  939201 |  935529 |  936048 |  932721

###  2.2 Model and simulation

We deal with the Hamiltonian of the Ising model given by

$\mathcal{H}=-J\sum_{\langle i,j\rangle}s_{i}s_{j},\quad s_{i}=\pm 1,$ (1)

where $J$ is the coupling and $s_{i}$ is the spin at lattice site $i$. The
lengths of the bonds in the kite are not all equal; since we are discussing
the network structure of the lattice, we assume that the interaction constants
all take constant $J$. The summation is taken over nearest-neighbor pairs
$\langle i,j\rangle$. Here, we consider the case of the ferromagnetic coupling
($J>0$).

In performing the Monte Carlo simulation, we employ the multi-cluster spin-
flip SW algorithm [17] to escape from the problem of the long autocorrelation
time near the critical temperature. We also combine the replica exchange
method of parallel tempering [18]. We used the 20,000 Monte Carlo steps (MCS)
after discarding the first 2,000 MCS for the “small” and “medium” sizes. We
used the 10,000 MCS after the first 1,000 MCS for the “big” size, and the
5,000 MCS after the 500 MCS for the “bigger” size. Since the cluster flip
algorithm [17] is used, the time required for equilibration is very short.
Even for a largest system size, 100 steps are sufficient for equilibration. We
discarded 1/10th of the measurement time steps. For temperature dependence,
precise measurement is possible by applying the replica exchange method [18].
As the system size increases, it becomes necessary to take a narrower
temperature difference to exchange temperatures. The detailed balance
conditions are fulfilled for both algorithms We made five independent runs
with different random number sequences for each size to estimate numerical
errors. In most plots the error bars are within the size of the marks. This
study dealt with a sample of one configuration for each size. Because we are
cutting finite sizes from an infinite system, we can take some other
configurations. However, since we are dealing with sufficiently large sizes,
we expect the dependence on configuration to be very small.

We make the Monte Carlo updates for the spins on all the sites. However, we
take the average of physical quantities over only the spins on inner sites to
reduce surface effects. Surface effects are one of the important topics in the
study of critical phenomena, and numerical investigation of surface effects is
often performed based on the difference between periodic and free boundary
conditions [21]. In this study, the bulk properties are the main subject, and
we will consider the quantity excluding the surface layer. The number of the
inner sites for both original and dual lattices are tabulated in table 1, and
in the following, we denote $N$ for the number of inner sites.

##  3 Results

###  3.1 Results of magnetization

We present the simulation results of Tile($1,\sqrt{3}$) of 8 kites. We first
show the results of the magnetization per spin, $m=(1/N)\sum_{i}M_{i}$. We
plot the temperature dependence of the squared magnetization, $\langle
m^{2}\rangle$, in figure 3. We show the data of the four sizes, “small”,
“medium”, “big”, and “bigger”, for both original and dual lattices. The
temperature is measured in units of $J$, that is, $J=1$. In figure 3, we
observe the spontaneous magnetization for lower temperature.

Figure 3:  Squared magnetization of the Ising model on the Smith-kite lattices
for both original and dual lattices.

Binder [22] pointed out that the moment ratio is effective for locating the
critical temperature based on the finite-size scaling (FSS) analysis. We plot
the temperature dependence of the moment ratio $\langle m^{4}\rangle/\langle
m^{2}\rangle^{2}$ in figure 4. We give the data of four sizes for both
original and dual lattices. This value approaches 1 for the ordered state
($T\to 0$), and approaches 3 for the disordered state ($T\to\infty$) because
of Gaussian fluctuations.

Figure 4:  Moment ratio of the Ising model on the Smith-kite lattices for both
original and dual lattices.

The correlation length $\xi$ divereges as the temperature approaches the
critical temperature as

$\xi(T)\propto|T-T_{c}|^{-\nu}$ (2)

with the correlation-length exponent $\nu$. The FSS [23, 24] postulates that
near the critical temperature the physical quantities of finite systems with
the linear size $L$ are scaled by $\xi/L$. The magnetization scales as

$\langle m^{2}\rangle=L^{-2\beta/\nu}f(tL^{1/\nu}),$ (3)

where $\beta$ is the spontaneous-magnetization exponent and $t=T-T_{c}$. Then,
the moment ratio becomes a single scaling variable such that

$\langle m^{4}\rangle/\langle m^{2}\rangle^{2}=\tilde{f}(tL^{1/\nu}).$ (4)

It means that the plots $\langle m^{4}\rangle/\langle m^{2}\rangle^{2}$ of
various sizes intersect at a single point, which leads to the estimate of the
critical temperature. It is to be noted that there are corrections to FSS for
small sizes. We see the crossing of data of different sizes in figure 4. We
here make a comment on the error bars. At very high temperatures, this ratio
becomes a 0/0 quantity, so the error bars may become large. However, just
above the critical point, the error is still small even for the largest system
size.

Figure 5:  FSS of moment ratio for the Ising model on the Smith-kite lattices;
$t=T-T_{c}$ and $\nu=1$. The left panel is the FSS plot for the original
lattice, where we assume $T_{c}=2.405$, whereas the right panel is that for
the dual lattice, where we assume $T^{*}_{c}=2.143$.

Assuming the critical temperature $T_{c}$, we try the FSS plots for both
original and dual lattices in figure 5. The linear size $L$ is $N^{1/2}$, and
we use $\nu=1$ of the 2D Ising value of the regular lattice. The universality
of the critical exponents for quasiperiodic spin systems has been confirmed,
and even if we analyze the critical exponent as a fitting parameter, we will
obtain almost the same value as the exact solution using the FSS analysis. The
left panel is the FSS plot for the original lattice, where we assume
$T_{c}=2.405$, whereas the right panel is that for the dual lattice, where we
assume $T^{*}_{c}=2.143$.

For comparison, in figure 6 we show the FSS trial when $T_{c}$ is assumed to
be 2.404 and 2.406 in the case of the original lattice. For clarity, we use
larger marks for the “bigger” lattice. If you compare figure 6 with the left
panel of figure 5, we can see that $T_{c}=2.405$ is a good choice. The trial
plots of $T_{c}=2.404$ and $T_{c}=2.406$ in figure 6 are choices of data to
highlight the differences. We can obtain an estimate of the error in the
narrow interval. Thus, our estimates of the critical temperatures are

${\rm original\ lattice}\quad T_{c}=2.405\pm 0.0005$ (5)

and

${\rm dual\ lattice}\quad T^{*}_{c}=2.143\pm 0.0005.$ (6)

Figure 6:  FSS trial of moment ratio for the Ising model on the original
Smith-kite lattices. Choose $T_{c}=2.404$ (left panel) and $T_{c}=2.406$
(right panel). For clarity, larger marks are used for the “bigger” lattice.

Next, we show the FSS of magnetization, (3), in figure 7. Here, the critical
exponents are those of the regular lattice, $\nu=1$ and $\beta=1/8$, and
$T_{c}$ are chosen as the values of (5) and (6). From figure 7, we can see
that the FSS works quite well.

Figure 7:  FSS of magnetization for the Ising model on the Smith-kite lattice;
$t=T-T_{c}$, $\nu=1$ and $\beta=1/8$. The left panel is the FSS plot for the
original lattice, where we assume $T_{c}=2.405$, whereas the right panel is
that for the dual lattice, where we assume $T^{*}_{c}=2.143$.

We now consider the duality relation. A simple relation is known between the
partition functions of the Ising models on a regular lattice and its dual
lattice [20, 25]. The duality relation between the critical temperatures,
$T_{c}$ and $T^{*}_{c}$, is given by

$\sinh(2J/T_{c})\sinh(2J/T^{*}_{c})=1.$ (7)

Substituting the estimates of the critical temperatures on the Smith-kite
lattice (5), and its dual lattice, (6), we obtain
$\sinh(2/T_{c})\sinh(2/T^{*}_{c})=1.000\pm 0.001$. The duality relation is
well satisfied for the aperiodic monotile lattice.

###  3.2 Results of Energy

Another check of the duality is the duality relation for the critical energy,

$\frac{\epsilon(T_{c})}{\coth(2J/T_{c})}+\frac{\epsilon(T^{*}_{c})}{\coth(2J/T^{*}_{c})}=1,$
(8)

where $\epsilon(T_{c})$ is the nearest neighbor spin correlation at the
critical temperature. This quantity is related to the energy per spin at the
critical temperature by $\epsilon(T_{c})=-E(T_{c})/2N$.

Figure 8:  Energy of the Ising model on the Smith-kite lattices for both
original and dual lattices.

The temperature dependence of the energy of the Ising model on the Smith-kite
lattice and its dual lattice is shown in figure 8. In order to estimate the
energy at the critical temperature, we use the FSS for the critical energy,
$(E(T_{c},L)-E(T_{c}))/N\propto tL^{-(1-\alpha)/\nu}$. Here, $\alpha$ is the
specific-heat exponent; it is 0 for the 2D Ising model. The leading anomaly of
the energy is $\alpha=0$, although the specific heat, the temperature
derivative of the energy, has logarithmic divergence. Using the critical
energies for each size, we plot the left-hand side of (8),
$\epsilon(T_{c})/\coth(2J/T_{c})+\epsilon(T^{*}_{c})/\coth(2J/T^{*}_{c})$, as
a function of $1/L=1/\sqrt{N}$ in figure 9. The result of polynomial fitting
(on the second order in powers of $1/\sqrt{N}$) is also shown by the dotted
curve. Then we obtain the estimate of the left-hand side of (8) as
$\epsilon(T_{c})/\coth(2J/T_{c})$ \+ $\epsilon(T^{*}_{c})/\coth(2J/T^{*}_{c})$
= $1.000\pm 0.001$. We also evaluated the critical energies by the same
fitting procedure as $E(T_{c})/N=-1.319$ for the original lattice and
$E(T^{*}_{c})/N=-1.505$ for the dual lattice.

Figure 9:  FSS of the critical energies of the Ising model on the Smith-kite
lattices. We plot the left-hand side of (8) as a function of $1/L=1/\sqrt{N}$.
We also give the result of polynomial fitting by the dotted curve.

The present result clearly indicates that the duality relations for the
critical temperature and the critical energy hold for an aperiodic monotile
lattice. In the case of the Penrose lattice, Komura and Okabe [16] presented
the analytical argument for the duality relation based on the discussion by
Kramers and Wannier [20]. They followed a detailed argument by Syozi [25]. The
point is that the duality relation comes from the local dual transformation of
the lattice, not from the regularity of the lattice. Thus, the same story can
be applied to the present Smith-kite aperiodic lattice.

###  3.3 Comparison with other lattices

Here, we compare the present results of the critical temperature with those of
other lattices. For the Ising models on the Penrose lattice and its dual
lattices, the estimates of the critical temperatures by a high-precision study
[16] are

$T_{c}^{\rm Penrose}=2.39781\pm 0.00005$ (9)

and

$T_{c}^{\rm dual\ Penrose}=2.14987\pm 0.00005.$ (10)

In the simulation of the high-precision study [16], the periodic Penrose
lattice with a large unit cell was used. The formation of the Penrose tiling
is associated with the golden ratio, which is the irrational number
$\tau=(1+\sqrt{5})/2$. The periodic Penrose lattice is constructed by the
continued-fraction expansion of the irrational number. The number of the
lattice sites in the unit cell depend on the truncated level of the continued-
fraction expansion.

It is also instructive to compare the present data with the data of regular
lattices, where the coordination number is 4 on average. The critical
temperatures of the Ising model on the kagome lattice and its dual lattice,
the diced lattice, were exactly obtained [26, 27]. See also the calculation of
the Fisher’s zero for the kagome and diced lattices [28]. The obtained exact
critical temperatures are

$T_{c}^{\rm kagome}=\frac{4}{\ln(3+2\sqrt{3})}=2.143319\cdots$ (11)

and

$T_{c}^{\rm diced}=\frac{2}{\cosh^{-1}\frac{1+\sqrt{3}}{2}}=2.405457\cdots.$
(12)

Of course, these critical temperatures satisfy the duality relation (7).

In table 2, we compare the $T_{c}$’s of the Ising models on the
(quasi-)lattices with the 4 coordination number. We also give the information
on the lattice structure in table 2. The square lattice is self-dual, and we
added it for reference;

$T_{c}^{\rm square}=\frac{2}{\ln(1+\sqrt{2})}=2.26919\cdots.$ (13)

The Smith-kite lattice, Penrose lattice, diced lattice are bipartite lattices,
and take several values for the coordination number; there is no frustration
for antiferromagnetic systems. On the other hand, the dual Smith-kite lattice,
dual Penrose lattice, kagome lattice are not bipartite, and the coordination
number is fixed as 4. Because there are “odd-numbered” rings, there appears
frustration for antiferromagnetic systems. The duality relations are satisfied
for all the dual pairs. The $T_{c}$’s of the Smith-kite lattice, Penrose
lattice, diced lattice are almost the same, which reflects upon the similarity
of the lattices.

Table 2:  Critical temperatures of Ising models on various (quasi-)lattices whose (average) coordination number is 4. We also give the information on the lattice structure in the table. The numbers in the parentheses are the average coordination number.  lattice | Smith-kite | d-Smith-kite | Penrose | d-Penrose
---|---|---|---|---
coordination number | 3,4,6 (4) | 4 | 3,4,5,6,7 (4) | 4
bipartite | yes | no | yes | no
frustration (AF) | no | yes | no | yes
$T_{c}/J$ | 2.405 | 2.143 | 2.39781 | 2.14987
lattice | diced | kagome | square |
coordination number | 3,6 (4) | 4 | 4 |
bipartite | yes | no | yes |
frustration (AF) | no | yes | no |
$T_{c}/J$ | 2.405457 $\cdots$ | 2.143319 $\cdots$ | 2.26919 $\cdots$ |

##  4 Summary and discussions

To summarize, we have studied the Ising models on the Smith-kite lattice and
the dual Smith-kite lattice by means of the SW multi-cluster spin-flip
algorithm. We have estimated the critical temperatures on these lattices, and
have shown that the finite-size scaling works quite well. The estimated
critical temperature on the Smith-kite lattice is $T_{c}/J=2.405\pm 0.0005$
and that on the dual Smith-kite lattice is $T^{*}_{c}/J=2.143\pm 0.0005$. We
have checked the duality relation. From the results of critical temperatures
on the Smith-kite lattice and its dual lattice, we have obtained
$\sinh(2J/T_{c})\sinh(2J/T^{*}_{c})=1.000\pm 0.001$. We have also shown the
duality relation for the energy, $\epsilon(T_{c})/\coth(2J/T_{c})$ \+
$\epsilon(T^{*}_{c})/\coth(2J/T^{*}_{c})$ = $1.000\pm 0.001$. These numerical
results clearly indicate that the duality relation is satisfied for an
aperiodic monotile lattice. We have compared the present results with those of
the Ising models on other lattices, including the Penrose lattice of
quasicrystal and the regular lattices (kagome and dice lattices). The average
coodination numbers of all these lattices are four.

We presented the results of Tile($1,\sqrt{3}$) with 8 kites. We also performed
the simulation for Tile($\sqrt{3},1$) with 10 kites. The results obtained are
in close numerical agreement. Although not in perfect agreement, the
difference is within 0.05% in the estimate of the critical temperatures. This
can be attributed to the fact that both are tiling families of the same
origin.

In the simulation of the Penrose lattice, the periodic lattice was used [14,
15, 16]. For the Smith-kite lattice, it will be helpful to use periodic
lattice for high-precision simulation. Similar procedure to approximate the
irrational number by a series of rational numbers could be applied to the case
of Smith-kite lattice.

From the viewpoint of statistical physics, we have considered the Smith-kite
lattice consisting of aperiodic monotiles. The study of electric systems, for
example, on this aperiodic lattice is desirable. Quite recently, the spectral
and transport properties of a tight-binding model defined on the Hat [29] and
the dimer model on the Hat [30] have been discussed.

Finally, we would like to comment on the study of tiling itself. The apeioric
Smith tiling still has much to be studied from a mathematical science
standpoint. The discovery of a new aperiodic monotile is also desirable.

## Data availability statement

All data that support the findings of this study are included within the
article.

## Acknowledgments

We thank Arnaud Cheritat for allowing us to use the program code he created.
This work was supported by JSPS KAKENHI Grant Number JP22K03472.

## ORCID iDs

Yutaka Okabe https://orcid.org/0000-0002-9434-5653

## References

## References

  * [1]
  * [2] Penrose R 1974 The role of aesthetics in pure and applied mathematical research Bull. Inst. Math. Appl. 10 266
  * [3] Levine D and Steinhardt P J 1984 Quasicrystals: A New Class of Ordered Structures Phys. Rev. Lett. 53 2477
  * [4] Shechtman D, Blech I, Gratias D and Cahn J W 1984 Metallic Phase with Long-Range Orientational Order and No Translational Symmetry Phys. Rev. Lett. 53 1951
  * [5] Smith D, Myers J S, Kaplan C S and Goodman-Strauss C 2023 An aperiodic monotile arXiv:2303.10798
  * [6] Smith D, Myers J S, Kaplan C S and Goodman-Strauss C 2023 A chiral aperiodic monotile arXiv:2305.17743
  * [7] Socolar J E S 2023 Quasicrystalline structure of the Smith monotile tilings arXiv:2305.01174
  * [8] Baake M, Gähler F and Sadun L 2023 Dynamics and topology of the hat family of tilings arXiv:2305.05639
  * [9] Akiyama S and Araki Y 2023 An alternative proof for an aperiodic monotile arXiv:2307.12322
  * [10] Deguchi K, Nakayama M, Matsukawa S, Imura K, Tanaka K, Ishimasa T and Sato N K 2015 Superconductivity of Au-Ge-Yb Approximants with Tsai-Type Clusters J. Phys. Soc. Japan 84 023705
  * [11] Kamiya K, Takeuchi T, Kabeya N, Wada N, Ishimasa T, Ochiai A, Deguchi K, Imura K and Sato N K 2018 Discovery of superconductivity in quasicrystal Nat. Commun. 9 154
  * [12] Fuchs J N and Vidal J 2016 Hofstadter buttery of a quasicrystal Phys. Rev. B 94 205437
  * [13] Nawa K, Avdeev M, Ishikawa A, Takakura H, Wang C W, Okuyama D, Murasaki R, Tamura R and Sato T J 2023 Magnetic properties of the quasicrystal approximant Au65Ga21Tb14 Phys. Rev. Mat. 7 054412
  * [14] Okabe Y and Niizeki K 1988 Monte Carlo Simulation of the Ising Model on the Penrose Lattice J. Phys. Soc. Japan 57 16
  * [15] Okabe Y and Niizeki K 1988 Duality in the Ising Model on the Quasicrystals J. Phys. Soc. Japan 57 1536
  * [16] Komura Y and Okabe Y 2016 High-Precision Monte Carlo Simulation of the Ising Models on the Penrose Lattice and the Dual Penrose Lattice J. Phys. Soc. Japan 85 044004\.
  * [17] Swendsen R H and Wang J S 1987 Nonuniversal critical dynamics in Monte Carlo simulations Phys. Rev. Lett. 58 86
  * [18] Hukushima K and Nemoto K 1996 Exchange Monte Carlo Method and Application to Spin Glass Simulations J. Phys. Soc. Japan 65 1604
  * [19] Cheritat A 2023 available at: https://www.math.univ-toulouse.fr/`~`cheritat/AppletsDivers/Monotile-4D-lift/3-outlines/
  * [20] Kramers H A and Wannier G H 1941 Statistics of the Two-Dimensional Ferromagnet. Part I Phys. Rev. 60 252
  * [21] Binder K 1983 Phase Transitions and Critical Phenomena vol 8, ed C Domb and J Lebowitz (New York: Academic Press)
  * [22] Binder K 1981 Finite size scaling analysis of ising model block distribution functions Z. Phys. B: Condens. Matter 43 119
  * [23] Barber M N 1983 Phase Transitions and Critical Phenomena vol 8, ed C Domb and J Lebowitz (New York: Academic Press)
  * [24] Cardy J L 1988 Finite Size Scaling ed J L Cardy (Amsterdam: North-Holland)
  * [25] Syozi I 1972 Phase Transitions and Critical Phenomena ed C Domb and M S Green (Academic Press, London) Vol. 1
  * [26] Syozi I 1950 Statistics of Kagome Lattice Prog. Theor. Phys. 6 306
  * [27] Kano K and Naya S 1953 Antiferromagnetism. The Kagome Ising Net Prog. Theor. Phys. 10 158
  * [28] Feldmann H, Shrock R and Tsai S H 1998 Complex-temperature partition function zeros of the Potts model on the honeycomb and kagome lattices Phys. Rev. E57 1335
  * [29] Schirmann J, Franca S, Flicker F and Grushin A G 2023 Physical properties of the Hat aperiodic monotile: Graphene-like features, chirality and zero-modes arXiv:2307.11054
  * [30] Singh S and Flicker F 2023 Exact Solution to the Quantum and Classical Dimer Models on the Spectre Aperiodic Monotiling arXiv:2309.14447