killcontents

# Resilient Distributed Optimization Algorithms for Resource Allocation

César A. Uribe†, Hoi-To Wai†, Mahnoosh Alizadeh CAU and HTW have contributed
equally. CAU is with LIDS, MIT, Cambridge, MA, USA. HTW is with Dept. of SEEM,
CUHK, Shatin, Hong Kong. MA is with Dept. of ECE, UCSB, Santa Barbara, CA,
USA. This work is partially supported by UCOP Grant LFR-18-548175 and CUHK
Direct Grant #4055113. E-mails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Distributed algorithms provide flexibility  over centralized algorithms for
resource allocation problems, e.g., cyber-physical systems. However, the
distributed nature of these algorithms often makes the systems susceptible to
man-in-the-middle attacks, especially when messages are transmitted between
price-taking agents and a central coordinator. We propose a resilient strategy
for distributed algorithms under the framework of primal-dual distributed
optimization. We formulate a robust optimization model that accounts for
Byzantine attacks on the communication channels between agents and
coordinator. We propose a resilient primal-dual algorithm using state-of-the-
art robust statistics methods. The proposed algorithm is shown to converge to
a neighborhood of the robust optimization model, where the neighborhood’s
radius is proportional to the fraction of attacked channels.

##  1 Introduction

Consider the following multi-agent optimization problem involving the average
of parameters in the constraints:

$\begin{array}[]{rl}\displaystyle\min_{\bm{\theta}_{i}\in\mathbb{R}^{d},\forall
i}&U(\bm{\theta})\mathrel{\mathop{:}}=\frac{1}{N}\sum_{i=1}^{N}U_{i}(\bm{\theta}_{i})\\\
{\rm s.t.}&g_{t}\left(\frac{1}{N}\sum_{i=1}^{N}\bm{\theta}_{i}\right)\leq
0,~{}t=1,...,T,\vspace{.1cm}\\\
&\bm{\theta}_{i}\in{\mathcal{}C}_{i},~{}i=1,...,N,\end{array}$ (1)

where both $U_{i}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ and
$g_{t}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ are continuously differentiable,
convex functions, and ${\mathcal{}C}_{i}$ is a compact convex set in
$\mathbb{R}^{d}$. We let ${\bm{0}}\in{\mathcal{}C}_{i}$ and

$\max_{\bm{\theta},\bm{\theta}^{\prime}\in{\mathcal{}C}_{i}}\|\bm{\theta}-\bm{\theta}^{\prime}\|\leq
R,~{}i=1,...,N,$ (2)

such that $R$ is an upper bound on the diameters of ${\mathcal{}C}_{i}$.

Problem (1) arises in many resource allocation problems with a set of
potentially _nonlinear_ constraints on the amount of allowable resources, see
Section 1.1 for a detailed exploration.

We consider a system where there exists a central coordinator and $N$ agents.
In this context, the function $U_{i}(\bm{\theta}_{i})$ and parameter
$\bm{\theta}_{i}$ are the utility of the $i$th agent and the resource
controlled by agent $i$, respectively. As the agents work independently, it is
desirable to design algorithms that allow the $N$ agents to solve (1)
cooperatively through communication with the central coordinator. Among
others, the primal-dual optimization methods [1] have been advocated as they
naturally give rise to decomposable algorithms that favor distributed
implementation [2]. In addition to their practical success, these methods are
supported by strong theoretical guarantees where fast convergence to an
optimal solution of (1) is well established. However, the distributed nature
of these methods also exposes the system to vulnerabilities not faced by
traditional centralized systems. Precisely, existing algorithms assume the
agents, and the communication links between central server and agents, to be
_completely trustworthy_. However, an attacker can take over a sub-system
operated by the agents, and deliberately edit the messages in these
communication links, i.e., a Byzantine attack. This might result in an
unstable system with possible damages to hardware and the system overall.

In this paper, we propose strategies for securing primal-dual distributed
algorithms, e.g., in [1], tailored to solving a relaxed version of the
resource allocation problem (1). A key observation is that the existing
algorithms depend on reliably computing the average of a set of parameter
vectors, $\\{\bm{\theta}_{i}\\}_{i=1}^{N}$, transmitted by the agents. As a
remedy, we apply robust statistics techniques as a subroutine, therefore
proposing a _resilient_ distributed algorithm that is proven to converge to a
neighborhood of the optimal solution of a robust version of (1).

Vulnerabilities of various types of distributed algorithms have been
identified and addressed in a number of recent studies. Relevant examples are
[3, 4, 5, 6, 7] which study secure decentralized algorithms on a general
network topology but consider consensus-based optimization models. Moreover,
[8, 9, 10] consider a similar optimization architecture as this paper, yet
they focus on securing distributed algorithms for machine learning tasks which
assumes i.i.d. functions, a fundamentally different setting from the current
paper. Our work is also related to the literature on robust statistics [11,
12], and particularly, with the recently rekindled research efforts on high
dimensional robust statistics [13, 14, 15]. These works will be the working
horse for our attack resilient algorithm.

Our contributions and organization are as follows. First, we derive a formal
model for attack resilient resource allocation via a conservative
approximation for the robust optimization problem [cf. Section 3]. Second, we
apply and derive new robust estimation results to secure distributed resource
allocation algorithms [cf. Section 4]. Third, we provide a non-asymptotic
convergence guarantee of the proposed attack resilient algorithm [cf. Section
4.1]. In particular, our algorithm is shown to converge to a
${\mathcal{}O}(\alpha^{2})$ neighborhood to the optimal solution of (1), where
$\alpha\in[0,\frac{1}{2})$ is the fraction of attacked links.

Notations. Unless otherwise specified, $\|\cdot\|$ denotes the standard
Euclidean norm. For any $N\in\mathbb{N}$, $[N]$ denotes the finite set
$\\{1,...,N\\}$.

###  1.1 Motivating Examples

Our set-up here can be employed in a wide range of optimization problems for
resource allocation and networked control in multi-agent systems, e.g., in the
pioneering example of congestion control in data networks [16, 17]; in
determining the optimal price of electricity and enabling more efficient
demand supply balancing (a.k.a. demand response) in smart power distribution
systems [18, 19]; in managing user transmit powers and data rates in wireless
cellular networks [20]; in determining optimal caching policies by content
delivery networks [21]; in optimizing power consumption in wireless sensor
networks with energy-restricted batteries [22, 23]; and in designing
congestion control systems in urban traffic networks [24]. These examples
would have different utility functions and constraint sets that can be handled
through our general formulation in (1). For example, in the power/rate control
problem in data networks, the cost functions are usually logarithmic functions
associated with rate $\theta_{i}$, e.g.,
$U_{i}(\theta_{i})=-\beta_{i}\log(\theta_{i})$. In demand response
applications in power distribution systems, the utilities capture the users’
benefits from operating their electric appliances under different settings.
For example, we can capture the cost function of temperature $\theta_{i}$
controlled by a price-responsive air conditioner as
$U_{i}(\theta_{i})=b_{i}(\theta_{i}-\theta_{\rm{comf}})^{2}-c_{i}$ [19]. In
terms of constraints, our general nonlinear constraint formulation can not
only capture common linear resource constraints such as link capacity in data
networks [16, 17], but can also handle important non-linear constraints
arising in many different applications. For example, in radial power
distribution systems, nonlinear convexified power flow constraints can be
included for distributed demand response optimization (to see a description of
distribution system power flow constraints, see, e.g., [25, 26]). This can
enable our algorithm to perform demand supply balancing in power disribution
systems in a distributed and resilient fashion.

##  2 Primal-dual Algorithm for Resource Allocation

This section reviews the basic primal-dual algorithm for resource allocation.
Let $\bm{\lambda}\in\mathbb{R}_{+}^{T}$ be the dual variable. We consider the
Lagrangian function of (1):

$\begin{split}&{\mathcal{}L}(\\{\bm{\theta}_{i}\\}_{i=1}^{N};\bm{\lambda})\mathrel{\mathop{:}}=\frac{1}{N}\sum_{i=1}^{N}U_{i}(\bm{\theta}_{i})+\sum_{t=1}^{T}\lambda_{t}\\!~{}g_{t}\Big{(}\frac{1}{N}\sum_{i=1}^{N}\bm{\theta}_{i}\Big{)}.\end{split}\vspace{-.1cm}$

Assuming strong duality holds (e.g., under the Slater’s condition), solving
problem (1) is equivalent to solving its dual problem:

$\max_{\bm{\lambda}\in\mathbb{R}_{+}^{T}}~{}\min_{\bm{\theta}_{i}\in{\mathcal{}C}_{i},\forall
i}~{}{\mathcal{}L}(\\{\bm{\theta}_{i}\\}_{i=1}^{N};\bm{\lambda}).$ (P)

For a given $\bm{\lambda}$, the inner minimization of (P) is known as the
Lagrangian relaxation of (1), which can be interpreted as a _penalized_
resource allocation problem [19].

In a distributed setting, the goal is to solve (1) where the agents only
observe a _pricing signal_ received from the central coordinator, and this
pricing signal is to be updated iteratively at the central coordinator. As
suggested in [1], we apply the primal-dual algorithm (PDA) to a regularized
version of (P). Let us define

$\begin{split}&{\mathcal{}L}_{\upsilon}(\\{\bm{\theta}_{i}\\}_{i=1}^{N};\bm{\lambda})\mathrel{\mathop{:}}=\\\
&\textstyle{\mathcal{}L}(\\{\bm{\theta}_{i}\\}_{i=1}^{N};\bm{\lambda})+\frac{\upsilon}{2N}\sum_{i=1}^{N}\|\bm{\theta}_{i}\|^{2}-\frac{\upsilon}{2}\|\bm{\lambda}\|^{2},\end{split}$
(3)

such that ${\mathcal{}L}_{\upsilon}(\cdot)$ is $\upsilon$-strongly convex and
$\upsilon$-strongly concave in $\\{\bm{\theta}_{i}\\}_{i=1}^{N}$ and
$\bm{\lambda}$, respectively. Let $k\in\mathbb{Z}_{+}$ be the iteration index,
$\gamma>0$ be the step sizes, the PDA recursion is described by:

$\displaystyle\bm{\theta}_{i}^{(k+1)}=$ (4a)
$\displaystyle~{}~{}~{}~{}{\mathcal{}P}_{{\mathcal{}C}_{i}}\big{(}\bm{\theta}_{i}^{(k)}-\gamma\\!~{}{\nabla}_{\bm{\theta}_{i}}{\mathcal{}L}_{\upsilon}(\\{\bm{\theta}_{i}^{(k)}\\}_{i=1}^{N};\bm{\lambda}^{(k)})\big{)},\forall~{}i\in[N]$
$\displaystyle\bm{\lambda}^{(k+1)}=\big{[}\bm{\lambda}^{(k)}+\gamma\\!~{}{\nabla}_{\bm{\lambda}}{\mathcal{}L}_{\upsilon}(\\{\bm{\theta}_{i}^{(k)}\\}_{i=1}^{N};\bm{\lambda}^{(k)})\big{]}_{+}$
(4b)

where ${\mathcal{}P}_{{\mathcal{}C}_{i}}(\cdot)$ is the Euclidean projection
operator, $[\cdot]_{+}$ denotes $\max\\{0,\cdot\\}$, and the gradients are:

$\begin{split}&{\nabla}_{\bm{\theta}_{i}}{\mathcal{}L}_{\upsilon}(\\{\bm{\theta}_{i}^{(k)}\\}_{i=1}^{N};\bm{\lambda}^{(k)})=\textstyle\frac{1}{N}\Big{(}{\nabla}_{\bm{\theta}_{i}}U_{i}(\bm{\theta}_{i}^{(k)})+\upsilon\\!~{}\bm{\theta}_{i}^{(k)}\\\
&\hskip
45.52458pt\textstyle+\sum_{t=1}^{T}\lambda_{t}^{(k)}{\nabla}_{\bm{\theta}}g_{t}(\bm{\theta})\Big{|}_{\bm{\theta}=\frac{1}{N}\sum_{i=1}^{N}\bm{\theta}_{i}^{(k)}}\Big{)},\\\\[-17.07182pt]
\end{split}$ (5)
$\begin{split}&\big{[}{\nabla}_{\bm{\lambda}}{\mathcal{}L}_{\upsilon}(\\{\bm{\theta}_{i}^{(k)}\\}_{i=1}^{N};\bm{\lambda}^{(k)})\big{]}_{t}=g_{t}\Big{(}{\textstyle\frac{1}{N}\sum_{i=1}^{N}}\bm{\theta}_{i}^{(k)}\Big{)}-\upsilon\\!~{}\lambda_{t}^{(k)},\end{split}$
(6)

for all $i$, $t$. We denoted $[{\bm{x}}]_{t}$ as the $t$th element of
${\bm{x}}\in\mathbb{R}^{T}$. In particular, observe that (4) performs a
projected gradient descent/ascent on the primal/dual variables.

From the above, both gradients with respect to (w.r.t.) $\bm{\theta}_{i}$ and
$\lambda_{t}$ depend only on the average parameter
$\overline{\bm{\theta}}^{(k)}\mathrel{\mathop{:}}=\frac{1}{N}\sum_{i=1}^{N}\bm{\theta}_{i}^{(k)}$.
We summarize the primal dual distributed resource allocation (PD-DRA)
procedure in Algorithm 1. In addition to solving the general problem (1),
Algorithm 1 also serves as a general solution method to popular resource
allocation problems [19].

Algorithm 1 PD-DRA Procedure.

1: for $k=1,2,...$ do

2: _(Message exchanges stage)_ :

  1. (a)

Central coordinator receives $\\{\bm{\theta}_{i}^{(k)}\\}_{i=1}^{N}$ from
agents and computes $\overline{\bm{\theta}}^{(k)}$,
$\\{{\nabla}_{\bm{\theta}}g_{t}(\overline{\bm{\theta}}^{(k)})\\}_{t=1}^{T}$.

  2. (b)

Central coordinator broadcasts the vectors $\overline{\bm{\theta}}^{(k)}$,
$\overline{\bm{g}}^{(k)}\mathrel{\mathop{:}}=\sum_{t=1}^{T}\lambda_{t}^{(k)}{\nabla}_{\bm{\theta}}g_{t}(\overline{\bm{\theta}}^{(k)})$
to agents.

3: _(Computation stage)_ :

  1. (a)

Agent $i$ computes the update for $\bm{\theta}_{i}^{(k+1)}$ according to (4a)
using the received $\overline{\bm{\theta}}^{(k)}$.

  2. (b)

The central coordinator computes the update for $\bm{\lambda}^{(k+1)}$
according to (4b).

4: end for

As the regularized primal-dual problem is strongly convex/concave in
primal/dual variables, Algorithm 1 converges linearly to an optimal solution
[1]. To study this, let us denote
${\bm{z}}^{(k)}=(\\{\bm{\theta}_{i}^{(k)}\\}_{i=1}^{N},\bm{\lambda}^{(k)})$ as
the primal-dual variable at the $k$th iteration,

$\bm{\Phi}({\bm{z}}^{(k)})\mathrel{\mathop{:}}=\left(\begin{array}[]{c}{\nabla}_{\bm{\theta}}{\mathcal{}L}_{\upsilon}(\\{\bm{\theta}_{i}^{(k)}\\}_{i=1}^{N},\bm{\lambda}^{(k)})\\\
{\nabla}_{\bm{\lambda}}{\mathcal{}L}_{\upsilon}(\\{\bm{\theta}_{i}^{(k)}\\}_{i=1}^{N},\bm{\lambda}^{(k)})\end{array}\right).\vspace{-.1cm}$
(7)

######  Fact 1.

[1, Theorem 3.5] Assume that the map $\bm{\Phi}({\bm{z}}^{(k)})$ is $L_{\Phi}$
Lipschitz continuous. For all $k\geq 1$, we have

$\|{\bm{z}}^{(k+1)}-{\bm{z}}^{\star}\|^{2}\leq(1-2\gamma\upsilon+\gamma^{2}L_{\Phi}^{2})\\!~{}\|{\bm{z}}^{(k)}-{\bm{z}}^{\star}\|^{2}\;,$
(8)

where ${\bm{z}}^{\star}$ is a saddle point to the regularized version of (P).
Setting $\gamma=\upsilon/L_{\Phi}^{2}$ gives
$\|{\bm{z}}^{(k+1)}-{\bm{z}}^{\star}\|^{2}\leq\big{(}1-\upsilon^{2}/L_{\Phi}^{2}\big{)}\|{\bm{z}}^{(k)}-{\bm{z}}^{\star}\|^{2}$,
$\forall~{}k\geq 1$.

##  3 Problem Formulation

Despite the simplicity and the strong theoretical guarantee, the PD-DRA method
is susceptible to attacks on the channels between the central coordinator and
the agents, as described below.

…Central CoordinatorAgent $i$Agent $j$get $\overline{\bm{g}}^{(k)}$send
$\bm{\theta}_{i}^{(k)}$Attacked!

Figure 1: Illustrating the PD-DRA algorithm under attack. The uplink for agent
$j$ is compromised such that the correct $\bm{\theta}_{j}^{(k)}$ is not
transmitted to the central node. The up/downlink for agent $i$ are operating
properly.

Attack Model. We consider a situation when _uplink_ channels between agents
and the central coordinator are compromised [see Fig. 1]. Let
${\mathcal{}A}\subset[N]$ be the set of _compromised uplink channels_ , whose
identities are unknown to the central coordinator. We define
${\mathcal{}H}\mathrel{\mathop{:}}=[N]\setminus{\mathcal{}A}$ as the set of
trustworthy channels. At iteration $k$, instead of receiving
$\bm{\theta}_{i}^{(k)}$ from each agent $i\in[N]$ [cf. Step 2(a)], the central
coordinator receives the following messages:

${\bm{r}}_{i}^{(k)}=\begin{cases}\bm{\theta}_{i}^{(k)},&\text{if}~{}i\in{\mathcal{}H},\\\
{\bm{b}}_{i}^{(k)},&\text{if}~{}i\in{\mathcal{}A}.\end{cases}\vspace{-.2cm}$
(9)

We focus on a Byzantine attack scenario such that the messages,
${\bm{b}}_{i}^{(k)}$, communicated on the attacked channels can be arbitrary.
Under such scenario, if the central coordinator forms the naive average
$\widehat{\bm{\theta}}^{(k)}=1/N\sum_{i=1}^{N}{\bm{r}}_{i}^{(k)}$ and computes
the gradients ${\nabla}g_{t}(\widehat{\bm{\theta}}^{(k)})$ accordingly, this
may result in uncontrollable error since the deviation
$\widehat{\bm{\theta}}^{(k)}-(1/N)\sum_{i=1}^{N}\bm{\theta}_{i}^{(k)}$ can be
arbitrarily large. It is anticipated that the PD-DRA method would not provide
a solution to the regularized version of (P).

Robust Optimization Model. In light of the Byzantine attack, it is impossible
to optimize the original problem (P) since the contribution from
$U_{i}(\cdot):i\in{\mathcal{}A}$ becomes unknown to the central coordinator.
As a compromise, we focus on optimizing the cost function of agents with
trustworthy uplinks and the following robust optimization problem as our
target model:

$\displaystyle\min_{\begin{subarray}{c}\bm{\theta}_{i}\in{\mathcal{}C}_{i},i\in{\mathcal{}H}\end{subarray}}$
$\displaystyle~{}{\textstyle\frac{1}{|{\mathcal{}H}|}\sum_{i\in{\mathcal{}H}}U_{i}(\bm{\theta}_{i})}$
(10a) $\displaystyle~{}{\rm s.t.}$
$\displaystyle\displaystyle\max_{\bm{\theta}_{j}\in{\mathcal{}C}_{j},j\in{\mathcal{}A}}~{}g_{t}\Big{(}{\textstyle\frac{1}{N}\sum_{i=1}^{N}\bm{\theta}_{i}}\Big{)}\leq
0,~{}\forall~{}t,\vspace{.1cm}$ (10b)

note that $\\{\bm{\theta}_{j}\\}_{j\in{\mathcal{}A}}$ is taken away from the
decision variables and we have included (10b) to account for the _worst case_
scenario for the resource usage of the agents with compromised uplinks. This
is to ensure that the physical operation limit of the system will not be
violated under attack. Consider the following assumption which will be assumed
throughout the paper:

######  H 1.

For all $\bm{\theta}\in\mathbb{R}^{d}$, the gradient of $g_{t}$ is bounded
with $\|{\nabla}g_{t}(\bm{\theta})\|\leq B$ and is $L$-Lipschitz continuous.

We define

$\overline{g}_{t}(\bm{\theta})\mathrel{\mathop{:}}=g_{t}(\bm{\theta})+{\textstyle\frac{|{\mathcal{}A}|}{N}}\big{(}RB+{\textstyle\frac{1}{2}}LR^{2}\big{)},$
(11)

######  Lemma 1.

Under H1. The following problem yields a _conservative_ approximation of (10),
i.e., its feasible set is a subset of the feasible set of (10):

$\begin{array}[]{rl}\displaystyle\min_{\bm{\theta}_{i}\in{\mathcal{}C}_{i},i\in{\mathcal{}H}}&\frac{1}{|{\mathcal{}H}|}\sum_{i\in{\mathcal{}H}}U_{i}(\bm{\theta}_{i})\vspace{.1cm}\\\
{\rm
s.t.}&\displaystyle\overline{g}_{t}\left({\textstyle\frac{1}{N}\sum_{i\in{\mathcal{}H}}\bm{\theta}_{i}}\right)\leq
0,~{}\forall~{}t\in[T],\end{array}$ (12)

Similar to PD-DRA, we define the regularized Lagrangian function of (12) as:

$\begin{split}&\overline{\mathcal{}L}_{\upsilon}(\\{\bm{\theta}_{i}\\}_{i\in{\mathcal{}H}};\bm{\lambda};{\mathcal{}H})\\\
&\mathrel{\mathop{:}}={\textstyle\frac{1}{|{\mathcal{}H}|}\sum_{i\in{\mathcal{}H}}}U_{i}(\bm{\theta}_{i})+{\textstyle\sum_{t=1}^{T}}\lambda_{t}\\!~{}\overline{g}_{t}\left({\textstyle\frac{1}{N}\sum_{i\in{\mathcal{}H}}\bm{\theta}_{i}}\right)\\\
&\hskip
14.22636pt\textstyle+\frac{\upsilon}{2|{\mathcal{}H}|}\sum_{i\in{\mathcal{}H}}\|\bm{\theta}_{i}\|^{2}-\frac{\upsilon}{2}\|\bm{\lambda}\|^{2}.\end{split}$
(13)

Again, the regularized Lagrangian function is $\upsilon$-strongly convex and
concave in $\bm{\theta}$ and $\bm{\lambda}$, respectively.

Our main task is to tackle the following modified problem of (P) under
Byzantine attack on (some of) the uplinks:

$\max_{\bm{\lambda}\in\mathbb{R}_{+}^{T}}\min_{\bm{\theta}_{i}\in{\mathcal{}C}_{i},\forall
i\in{\mathcal{}H}}~{}\overline{\mathcal{}L}_{\upsilon}(\\{\bm{\theta}_{i}\\}_{i\in{\mathcal{}H}};\bm{\lambda};{\mathcal{}H}),$
(P’)

and we let
$\widehat{\bm{z}}^{\star}=(\widehat{\bm{\theta}}^{\star},\widehat{\bm{\lambda}}^{\star})$
be the optimal solution to (P’). Notice that (P’) bears a similar form as (P)
and thus one may apply the PD-DRA method to the former naturally. However,
such application requires the central coordinator to compute the sample
average

$\textstyle\overline{\bm{\theta}}^{(k)}_{\mathcal{}H}\mathrel{\mathop{:}}=\frac{1}{|{\mathcal{}H}|}\sum_{i\in{\mathcal{}H}}\bm{\theta}_{i}^{(k)},$
(14)

at each iteration. However, the above might not be computationally feasible
under the attack model, since the central coordinator is oblivious to the
identity of ${\mathcal{}H}$. This is the main objective in the design of our
scheme.

##  4 Robust Distributed Resource Allocation

In this section, we describe two estimators for approximating
$\overline{\bm{\theta}}^{(k)}_{\mathcal{}H}$ [cf. (14)] from the received
messages (9) without knowing the identity of links in ${\mathcal{}H}$. To
simplify notations, we define $\alpha\geq|{\mathcal{}A}|/N$ as a known upper
bound to the fraction of compromised channels and assume $\alpha<1/2$ where
less than half of the channels are compromised.

As discussed after (14), the problem at hand is _robust mean estimation_ ,
whose applications to robust distributed optimization has been considered in
the machine learning literature [9, 10, 14] under the assumption that the
‘trustworthy’ signals are drawn i.i.d. from a Gaussian distribution. Our
setting is different since the signals $\bm{\theta}_{i}^{(k)}$,
$i\in{\mathcal{}H}$ are variables from the previous iteration whose
distribution is non-Gaussian in general. Our analysis will be developed
without such assumption on the distribution.

We first consider a simple median-based estimator applied to each coordinate
$j=1,...,d$. First, define the coordinate-wise median as:

$\big{[}{\bm{\theta}}_{\sf med}^{(k)}\big{]}_{j}={\sf
med}\big{(}\\{[{\bm{r}}_{i}^{(k)}]_{j}\\}_{i=1}^{N}\big{)},$ (15)

where ${\sf med}(\cdot)$ computes the median of the operand. Then, our
estimator is computed as the mean of the nearest $(1-\alpha)N$ neighbors of
$\big{[}{\bm{\theta}}_{\sf med}^{(k)}\big{]}_{j}$. To formally describe this,
let us define:

${\mathcal{}N}_{j}^{(k)}=\\{i\in[N]:\big{|}\big{[}{\bm{r}}_{i}^{(k)}-{\bm{\theta}}_{\sf
med}^{(k)}\big{]}_{j}\big{|}\leq r_{j}^{(k)}\\},$ (16)

where $r_{j}^{(k)}$ is chosen as $|{\mathcal{}N}_{j}^{(k)}|=(1-\alpha)N$. Our
estimator is:

$\textstyle[\widehat{\bm{\theta}}^{(k)}_{{\mathcal{}H}}]_{j}=\frac{1}{(1-\alpha)N}\sum_{i\in{\mathcal{}N}_{j}^{(k)}}[{\bm{r}}_{i}^{(k)}]_{j}.$
(17)

The following bounds the performance of (17).

######  Proposition 1.

Suppose that
$\max_{i\in{\mathcal{}H}}\big{\|}\bm{\theta}_{i}^{(k)}-\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}\big{\|}_{\infty}\leq
r$, then for any $\alpha\in(0,\frac{1}{2})$, it holds that

$\big{\|}\widehat{\bm{\theta}}^{(k)}_{\mathcal{}H}-\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}\big{\|}\leq\frac{\alpha}{1-\alpha}\Big{(}2+\sqrt{\frac{(1-\alpha)^{2}}{1-2\alpha}}\Big{)}r\sqrt{d}\;.\vspace{-.1cm}$
(18)

Under mild assumptions, the condition
$\max_{i\in{\mathcal{}H}}\big{\|}\bm{\theta}_{i}^{(k)}-\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}\big{\|}_{\infty}\leq
r$ can be satisfied with $r=\Theta(R)$, as implied by the compactness of
${\mathcal{}C}_{i}$ [cf. (2)]. Moreover, for sufficiently small $\alpha$, the
right hand side on (18) can be approximated by ${\mathcal{}O}(\alpha
R\sqrt{d})$. However, this median-based estimator may perform poorly for large
$\alpha$ (especially when
$\alpha\rightarrow{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{1}/{2}}$)
or dimension $d$. For these situations, a more sophisticated estimator is
required, as detailed next.

Algorithm 2 Recovering the mean of a set [15].

1: Input: $\alpha$, $\bm{\theta}_{i}^{(k)}$, $c_{i}=1$ for all $i=1,\ldots,N$,
and $\mathcal{B}=\\{1,\ldots,N\\}.$

2: Set $X_{\mathcal{B}}=[\cdots\bm{\theta}_{j}^{(k)}\cdots]^{\top}$ for
$j\in\mathcal{B}$ as the concatenated data matrix.

3: Let $Y\in\mathbb{R}^{d\times d}$ and
$W\in\mathbb{R}^{\mathcal{B}\times\mathcal{B}}$ be the maximizer/minimizer of
the saddle point problem  $\displaystyle\max_{\begin{subarray}{c}Y\succeq
0,\\\ \text{tr}(Y)\leq
1\end{subarray}}\min_{\begin{subarray}{c}0\leq{W_{ij}},\\\
W_{ij}\leq\frac{4{-}\alpha}{\alpha(2{+}\alpha){\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}n}},\\\
\sum_{j}W_{ji}=1\end{subarray}}\sum_{i\in\mathcal{B}}c_{i}(\bm{\theta}_{i}^{(k)}{-}X_{\mathcal{B}}w_{i})^{\top}Y(\bm{\theta}_{i}^{(k)}{-}X_{\mathcal{B}}w_{i})$

4: Let
$\tau_{i}^{*}=(\bm{\theta}_{i}^{(k)}-X_{\mathcal{B}}w_{i})^{\top}Y(\bm{\theta}_{i}^{(k)}-X_{\mathcal{B}}w_{i})$.

5: if $\sum_{i\in\mathcal{B}}c_{i}\tau_{i}^{*}>4n\sigma^{2}$ then

6: For $i\in\mathcal{B}$, replace $c_{i}$ with
$\big{(}1-\frac{\tau^{*}_{i}}{\max_{j\in\mathcal{B}}\tau^{*}_{j}}\big{)}c_{i}$.

7: For all $i$ with $c_{i}<\frac{1}{2}$, remove from $\mathcal{B}$.

8: Go back to Line 3.

9: end if

10: Set $W_{1}$ as the result of zeroing out all singular values of $W$ that
are greater than $0.9$.

11: Set $Z=X_{\mathcal{B}}W_{0}$ where
$W_{0}=(W{-}W_{1})(I{-}W_{1})^{\text{-}1}$.

12: if $\text{rank}(Z)=1$ then

13: Output: $\widehat{\bm{\theta}}_{\mathcal{}H}^{(k)}$ as average of the
columns of $X_{\mathcal{B}}$.

14: else

15: Output: $\widehat{\bm{\theta}}_{\mathcal{}H}^{(k)}$ as a column of $Z$ at
random.

16: end if

To derive the second estimator, we apply an auxiliary result from [15] which
provides an algorithm for estimating
$\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}$, as summarized in Algorithm 2. We
observe:

######  Proposition 2.

[15, Proposition 16] Suppose that
$\lambda_{\max}(\frac{1}{|{\mathcal{}H}|}\sum_{i\in{\mathcal{}H}}(\bm{\theta}_{i}^{(k)}-\overline{\bm{\theta}}_{\mathcal{}H}^{(k)})(\bm{\theta}_{i}^{(k)}-\overline{\bm{\theta}}_{\mathcal{}H}^{(k)})^{\top})\leq\sigma^{2}$.
For any $\alpha\in[0,\frac{1}{4})$, Algorithm 2 produces an output such that
$\|\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}-\widehat{\bm{\theta}}_{\mathcal{}H}^{(k)}\|={\mathcal{}O}(\sigma\sqrt{\alpha})$.

Again, similar to Proposition 1, the required condition above can be satisfied
with $\sigma=\Theta(R)$ under mild conditions. Thus, Proposition 2 states that
Algorithm 2 recovers $\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}$ up to an
error of ${\mathcal{}O}(\sqrt{\alpha}R)$. Note that this bound is dimension
free unlike the median estimator analyzed in Proposition 1.

The idea behind Algorithm 2 is to sequentially identify and remove the subset
of points that cannot be re-constructed from the mean of the data points. The
solution of the optimization problem in Line 3 measures how well can we
recover the data points as an average of the other $|{\mathcal{}H}|$ points.
The bounded sample variance assumption guarantees that one can re-construct
any element in the set ${\mathcal{}H}$ from its mean, thus, all such points
that introduce a large error, as quantified by $c_{i}$ can be safely removed.
Line 5 quantifies the magnitude of the optimal point of Line 3, and if such
value is large, such points that introduce a large error are down-weighted.
The process is repeated until the optimal solution of Line 3 is small enough
and a low rank approximation of the optimal $W$ can be used to return the
sample mean estimate.

Algorithm 3 Resilient PD-DRA

1: Input: Each agent has initial state $\bm{\theta}_{i}^{(0)}$.

2: for $k=1,2,...$ do

3: _(At the Central Coordinator)_ :

  1. (a)

Receives $\\{{\bm{r}}_{i}^{(k)}\\}_{i=1}^{N}$, see (9), from agents.

  2. (b)

Computes robust mean $\widehat{\bm{\theta}}_{\mathcal{}H}^{(k)}$ using the
estimator (17) or Algorithm 2.

  3. (c)

Broadcasts the vectors $\widehat{\bm{\theta}}_{\mathcal{}H}^{(k)}$ and
$\widehat{\bm{g}}^{(k)}_{\mathcal{}H}\mathrel{\mathop{:}}=\sum_{t=1}^{T}\lambda_{t}^{(k)}{\nabla}_{\bm{\theta}}\overline{g_{t}}(\widehat{\bm{\theta}}_{\mathcal{}H}^{(k)})$
to agents.

  4. (d)

Computes the update for $\bm{\lambda}^{(k+1)}$ with (20).

4: _(At each agent $i$)_:

  1. (a)

Agent receives $\widehat{\bm{\theta}}_{\mathcal{}H}^{(k)}$ and
$\widehat{\bm{g}}^{(k)}_{\mathcal{}H}$.

  2. (b)

Agent computes update for $\bm{\theta}_{i}^{(k+1)}$ with (19).

5: end for

Attack Resilient PD-DRA method. The above section provides the enabling tool
for developing the resilient PD-DRA method, which we summarize in Algorithm 3.
The algorithm behaves similarly as Algorithm 1 applied to (P’), with the
exception that the central coordinator is oblivious to ${\mathcal{}H}$, and it
uses a robust mean estimator to find an approximate average for the signals
sent through the trustworthy links. This approximate value is used to compute
the new price signals, and sent back to agents. In particular, the primal-dual
updates are described by

$\bm{\theta}_{i}^{(k+1)}={\mathcal{}P}_{{\mathcal{}C}_{i}}\big{(}\bm{\theta}_{i}^{(k)}-{\textstyle\frac{\gamma}{N}}\big{(}\widehat{\bm{g}}^{(k)}_{\mathcal{}H}+{\nabla}U_{i}(\bm{\theta}_{i}^{(k)})+\upsilon\bm{\theta}_{i}^{(k)}\big{)}\big{)},$
(19)
$\lambda_{t}^{(k+1)}=\big{[}\lambda_{t}^{(k)}+\gamma\big{(}\overline{g_{t}}({\textstyle\frac{|{\mathcal{}H}|}{N}}\widehat{\bm{\theta}}_{\mathcal{}H}^{(k)})-\upsilon\lambda_{t}^{(k)}\big{)}\big{]}_{+}.$
(20)

######  Lemma 2.

Algorithm 3 is a primal-dual algorithm [1] for (P’) with perturbed gradients:

$\displaystyle\widehat{\bm{g}}_{\bm{\theta}}^{(k)}$
$\displaystyle={\nabla}_{\bm{\theta}}\overline{\mathcal{}L}_{\upsilon}(\bm{\theta}^{(k)};\bm{\lambda}^{(k)};{\mathcal{}H})+{\bm{e}}_{\bm{\theta}}^{(k)},$
(21a) $\displaystyle\widehat{\bm{g}}_{\bm{\lambda}}^{(k)}$
$\displaystyle={\nabla}_{\bm{\lambda}}\overline{\mathcal{}L}_{\upsilon}(\bm{\theta}_{i}^{(k)};\bm{\lambda}^{(k)};{\mathcal{}H})+{\bm{e}}_{\bm{\lambda}}^{(k)},$
(21b)

where we have used concatenated variable as
$\bm{\theta}=(\bm{\theta}_{1},...,\bm{\theta}_{N})$ and
$\bm{\lambda}=(\lambda_{1},...,\lambda_{T})$. Under H1 and assuming that
$\lambda_{t}^{(k)}\leq\overline{\lambda}$ for all $k$, we have:

$\|{\bm{e}}_{\bm{\theta}}^{(k)}\|\leq\overline{\lambda}LT\|\widehat{\bm{\theta}}_{\mathcal{}H}^{(k)}-\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}\|,$
(22) $\|{\bm{e}}_{\bm{\lambda}}^{(k)}\|\leq
BT\|\widehat{\bm{\theta}}_{\mathcal{}H}^{(k)}-\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}\|.\vspace{-.1cm}$
(23)

The assumption $\lambda_{t}^{(k)}\leq\overline{\lambda}$ can be guaranteed
since
$\overline{g_{t}}({\textstyle\frac{|{\mathcal{}H}|}{N}}\widehat{\bm{\theta}}_{\mathcal{}H}^{(k)})$
is bounded.

###  4.1 Convergence Analysis

Finally, based on Lemma 2, we can analyze the convergence of Algorithm 3. Let
$\widehat{\bm{z}}^{\star}=(\widehat{\bm{\theta}}^{\star},\widehat{\bm{\lambda}}^{\star})$
be a saddle point of (P’) and define

$\overline{\bm{\Phi}}({\bm{z}}^{(k)})\mathrel{\mathop{:}}=\left(\begin{array}[]{c}{\nabla}_{\bm{\theta}}\overline{\mathcal{}L}_{\upsilon}(\\{\bm{\theta}_{i}^{(k)}\\}_{i\in{\mathcal{}H}},\bm{\lambda}^{(k)};{\mathcal{}H})\\\
-{\nabla}_{\bm{\lambda}}\overline{\mathcal{}L}_{\upsilon}(\\{\bm{\theta}_{i}^{(k)}\\}_{i\in{\mathcal{}H}},\bm{\lambda}^{(k)};{\mathcal{}H})\end{array}\right),$
(24)

######  Theorem 1.

Assume the map $\overline{\bm{\Phi}}({\bm{z}}^{(k)})$ is $L_{\Phi}$-Lipschitz
continuous. For Algorithm 3, for all $k\geq 0$ it holds

$\begin{split}&\|{\bm{z}}^{(k+1)}-\widehat{\bm{z}}^{\star}\|^{2}\leq\big{(}1-\gamma\upsilon+2\gamma^{2}L_{\Phi}^{2}\big{)}\|{\bm{z}}^{(k)}-\widehat{\bm{z}}^{\star}\|^{2}\\\
&+\big{(}\frac{4\gamma}{\upsilon}+2\gamma^{2}\big{)}E_{k}.\end{split}$ (25)

where
$E_{k}\mathrel{\mathop{:}}=\|{\bm{e}}_{\bm{\theta}}^{(k)}\|^{2}+\|{\bm{e}}_{\bm{\lambda}}^{(k)}\|^{2}$
is the total perturbation at iteration $k$. Moreover, if we choose
$\gamma<{\upsilon}/{2L_{\Phi}^{2}}$ and $E_{k}$ is upper bounded by
$\overline{E}$ for all $k$, then

$\limsup_{k\rightarrow\infty}\|{\bm{z}}^{(k)}-\widehat{\bm{z}}^{\star}\|^{2}\leq\frac{\frac{4}{\upsilon}+2\gamma}{\upsilon-2\gamma
L_{\Phi}^{2}}\overline{E}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>(26)

Combining the results from the last subsection, the theorem shows the desired
result that the resilient PD-DRA method converges to a
${\mathcal{}O}(\alpha^{2}R^{2}d)$ neighborhood of the saddle point of (P’), if
the median-based estimator (17) is used [or ${\mathcal{}O}(\alpha R^{2})$ if
Algorithm 2 is used], where $\alpha$ is the fraction of attacked uplink
channels. Moreover, it shows that the convergence rate to the neighborhood is
linear, which is similar to the classical PDA analysis [1].

Interestingly, Theorem 1 illustrates a trade-off in the choice of the step
size $\gamma$ between convergence speed and accuracy. In specific, (25) shows
that the rate of convergence factor $1-\gamma\upsilon+2\gamma^{2}L_{\Phi}^{2}$
can be minimized by setting $\gamma=\upsilon/(4L_{\Phi}^{2})$. However, in the
meantime, the asymptotic upper bound in (26) is increasing with $\gamma$ and
it can be minimized by setting $\gamma\rightarrow 0$. This will be a design
criterion to be explored in practical implementations.

##  5 Conclusions

In this paper, we studied the strategies for securing a primal-dual algorithm
for distributed resource allocation. Particularly, we propose a resilient
distributed algorithm based on primal-dual optimization and robust statistics.
We derive bounds for the performance of the studied algorithm and show that it
converges to a neighborhood of a robustified resource allocation problem when
the number of attacked channels is small.

## Acknowledgement

The authors would like to thank the anonymous reviewers for feedback, and Mr.
Berkay Turan (UCSB) for pointing out typos in the original submission of this
paper.

##  Appendix A Proof of Lemma 1

Since $g_{t}$ is $L_{t}$-smooth, the following holds

$\begin{split}&\textstyle
g_{t}\big{(}\frac{1}{N}\sum_{i=1}^{N}\bm{\theta}_{i}\big{)}\leq
g_{t}\big{(}\frac{1}{N}\sum_{i\in{\mathcal{}H}}\bm{\theta}_{i}\big{)}\\\
&\textstyle+\frac{1}{N}\sum_{j\in{\mathcal{}A}}\Big{\langle}\bm{\theta}_{j},{\nabla}g_{t}\big{(}\frac{1}{N}\sum_{i\in{\mathcal{}H}}\bm{\theta}_{i}\big{)}\Big{\rangle}+\frac{L_{t}}{2N^{2}}\big{\|}\sum_{j\in{\mathcal{}A}}\bm{\theta}_{j}\big{\|}^{2}\end{split}$
(27)

Furthermore, observe that the gradient of $g_{t}$ is uniformly bounded by $B$
and the diameter of ${\mathcal{}C}_{j}$ is $R$, then the right hand side of
(27) can be upper bounded by

$g_{t}\big{(}{\textstyle\frac{1}{N}\sum_{i\in{\mathcal{}H}}\bm{\theta}_{i}}\big{)}+{\textstyle\frac{1}{N}}\sum_{j\in{\mathcal{}A}}\big{(}RB+{\textstyle\frac{1}{2}}LR^{2}\big{)}$
(28)

As such, defining
$c_{t}\mathrel{\mathop{:}}=\frac{|{\mathcal{}A}|}{N}\big{(}RB+\frac{1}{2}LR^{2}\big{)}$,
it can be seen that

$g_{t}\big{(}{\textstyle\frac{1}{N}\sum_{i\in{\mathcal{}H}}\bm{\theta}_{i}}\big{)}+c_{t}\leq
0,~{}t=1,...,T$ (29)

implies the desired constraint in (10).

##  Appendix B Proof of Proposition 1

Fix any $j\in[d]$. The assumption implies that for all $i\in{\mathcal{}H}$,
one has

$\big{|}[\bm{\theta}_{i}^{(k)}-\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}]_{j}\big{|}\leq
r.$ (30)

We observe that $|{\mathcal{}H}|\geq(1-\alpha)N$. Applying [8, Lemma 1] shows
that the median estimator111At each coordinate, the median is the geometric
median estimator of one dimension in [8]. satisfies

$\big{|}\big{[}{\bm{\theta}}_{\sf
med}^{(k)}-\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}\big{]}_{j}\big{|}\leq(1-\alpha)\sqrt{\frac{1}{1-2\alpha}}r.$
(31)

The above implies that for all $i\in{\mathcal{}H}$, we have

$\big{|}[\bm{\theta}_{i}^{(k)}-{\bm{\theta}}_{\sf
med}^{(k)}]_{j}\big{|}\leq\Big{(}1+\sqrt{\frac{(1-\alpha)^{2}}{1-2\alpha}}\Big{)}r.$
(32)

This implies that
$r_{j}^{(k)}\leq\Big{(}1+\sqrt{\frac{(1-\alpha)^{2}}{1-2\alpha}}\Big{)}r$
since $|{\mathcal{}H}|\geq(1-\alpha)N$. We then bound the performance of
$\widehat{\bm{\theta}}^{(k)}$:

$\begin{split}&(1-\alpha)N[\widehat{\bm{\theta}}^{(k)}]_{j}=\sum_{i\in{\mathcal{}N}_{j}^{(k)}}[{\bm{r}}_{i}^{(k)}]_{j}\\\
&=\sum_{i\in{\mathcal{}H}}[{\bm{r}}_{i}^{(k)}]_{j}-\sum_{i\in{\mathcal{}H}\setminus{\mathcal{}N}_{j}^{(k)}}[{\bm{r}}_{i}^{(k)}]_{j}+\sum_{i\in{\mathcal{}A}\cap{\mathcal{}N}_{j}^{(k)}}[{\bm{r}}_{i}^{(k)}]_{j},\end{split}$
(33)

thus

$\begin{split}&(1-\alpha)N[\widehat{\bm{\theta}}^{(k)}-\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}]_{j}\\\
&=-\sum_{i\in{\mathcal{}H}\setminus{\mathcal{}N}_{j}^{(k)}}[{\bm{r}}_{i}^{(k)}-\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}]_{j}+\sum_{i\in{\mathcal{}A}\cap{\mathcal{}N}_{j}^{(k)}}[{\bm{r}}_{i}^{(k)}-\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}]_{j}.\end{split}$

Notice that $|{\mathcal{}A}\cap{\mathcal{}N}_{j}^{(k)}|\leq\alpha N$ and thus
$|{\mathcal{}H}\setminus{\mathcal{}N}_{j}^{(k)}|\leq\alpha N$. Gathering terms
shows

$|[\widehat{\bm{\theta}}^{(k)}-\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}]_{j}|\leq\frac{\alpha
N}{(1-\alpha)N}\Big{(}2+\sqrt{\frac{(1-\alpha)^{2}}{1-2\alpha}}\Big{)}r.$ (34)

The above holds for all $j\in[d]$. Applying the norm equivalence shows the
desired bound.

##  Appendix C Proof of Lemma 2

Comparing the equations in (21) with (19), (20), we identify that

$\big{[}{\bm{e}}_{\bm{\theta}}^{(k)}\big{]}_{i}=\frac{1}{N}{\sum_{t=1}^{T}}\lambda_{t}^{(k)}\big{(}{\nabla}_{\bm{\theta}}\overline{g_{t}}({\textstyle\frac{|{\mathcal{}H}|}{N}}\widehat{\bm{\theta}}_{\mathcal{}H}^{(k)})-{\nabla}_{\bm{\theta}}\overline{g_{t}}({\textstyle\frac{|{\mathcal{}H}|}{N}}\overline{\bm{\theta}}_{\mathcal{}H}^{(k)})\big{)}$
(35)
$\big{[}{\bm{e}}_{\bm{\lambda}}^{(k)}\big{]}_{t}=\overline{g}_{t}({\textstyle\frac{|{\mathcal{}H}|}{N}}\widehat{\bm{\theta}}_{\mathcal{}H}^{(k)})-\overline{g}_{t}({\textstyle\frac{|{\mathcal{}H}|}{N}}\overline{\bm{\theta}}_{\mathcal{}H}^{(k)})$
(36)

where $\big{[}{\bm{e}}_{\bm{\theta}}^{(k)}\big{]}_{i}$ denotes the $i$th block
of ${\bm{e}}_{\bm{\theta}}^{(k)}$. Using H1 and the said assumptions, we
immediately see that

$\|\big{[}{\bm{e}}_{\bm{\theta}}^{(k)}\big{]}_{i}\|\leq\overline{\lambda}\frac{LT}{N}\|\widehat{\bm{\theta}}_{\mathcal{}H}^{(k)}-\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}\|$
(37)

which then implies (22). H1 implies that $\overline{g}_{t}$ is $B$-Lipschitz
continuous, therefore

$|\big{[}{\bm{e}}_{\bm{\lambda}}^{(k)}\big{]}_{t}|\leq
B\|\widehat{\bm{\theta}}_{\mathcal{}H}^{(k)}-\overline{\bm{\theta}}_{\mathcal{}H}^{(k)}\|,$
(38)

which implies (23).

##  Appendix D Proof of Theorem 1

Based on Proposition 2, our idea is to perform a perturbation analysis on the
PDA algorithm. Without loss of generality, we assume $N=1$ and denote
$\bm{\theta}=\bm{\theta}_{1}$. To simplify notations, we also drop the
subscript, denote the modified and regularized Lagrangian function as
${\mathcal{}L}=\overline{\mathcal{}L}_{\upsilon}$. Furthermore, we denote the
saddle point to (P’) as
${\bm{z}}^{\star}=(\bm{\theta}^{\star},\bm{\lambda}^{\star})$.

Using the fact that
$\bm{\theta}^{\star}={\mathcal{}P}_{{\mathcal{}C}}(\bm{\theta}^{\star})={\mathcal{}P}_{{\mathcal{}C}}\big{(}\bm{\theta}^{\star}-\gamma{\nabla}_{\bm{\theta}}{\mathcal{}L}(\bm{\theta}^{\star},\bm{\lambda}^{\star})\big{)}$,
we observe that in the primal update:

$\displaystyle\|\bm{\theta}^{(k+1)}-\bm{\theta}^{\star}\|^{2}$
$\displaystyle\overset{(a)}{\leq}\|\bm{\theta}^{(k)}-\bm{\theta}^{\star}\|^{2}-2\gamma\langle\widehat{\bm{g}}_{\bm{\theta}}^{(k)}-{\nabla}_{\bm{\theta}}{\mathcal{}L}(\bm{\theta}^{\star},\bm{\lambda}^{\star}),\bm{\theta}^{(k)}-\bm{\theta}^{\star}\rangle$
$\displaystyle\hskip
22.76228pt+\gamma^{2}\|\widehat{\bm{g}}_{\bm{\theta}}^{(k)}-{\nabla}_{\bm{\theta}}{\mathcal{}L}(\bm{\theta}^{\star},\bm{\lambda}^{\star})\|^{2}$

where (a) is due to the projection inequality
$\|{\mathcal{}P}_{{\mathcal{}C}}({\bm{x}}-{\bm{y}})\|\leq\|{\bm{x}}-{\bm{y}}\|$.
Furthermore, using the Young’s inequality, for any $c_{0},c_{1}>0$, we have

$\displaystyle\|\bm{\theta}^{(k+1)}-\bm{\theta}^{\star}\|^{2}$
$\displaystyle\leq\|\bm{\theta}^{(k)}-\bm{\theta}^{\star}\|^{2}$
$\displaystyle\hskip
14.22636pt-2\gamma\langle{\nabla}_{\bm{\theta}}{\mathcal{}L}(\bm{\theta}^{(k)},\bm{\lambda}^{(k)})-{\nabla}_{\bm{\theta}}{\mathcal{}L}(\bm{\theta}^{\star},\bm{\lambda}^{\star}),\bm{\theta}^{(k)}-\bm{\theta}^{\star}\rangle$
$\displaystyle\hskip
14.22636pt{+\gamma^{2}(1+c_{0})\|{\nabla}_{\bm{\theta}}{\mathcal{}L}(\bm{\theta}^{(k)},\bm{\lambda}^{(k)})-{\nabla}_{\bm{\theta}}{\mathcal{}L}(\bm{\theta}^{\star},\bm{\lambda}^{\star})\|^{2}}$
$\displaystyle\hskip
14.22636pt{-2\gamma\langle{\bm{e}}_{\bm{\theta}}^{(k)},\bm{\theta}^{(k)}-\bm{\theta}^{\star}\rangle+\gamma^{2}\big{(}1+\frac{1}{c_{0}}\big{)}\|{\bm{e}}_{\bm{\theta}}^{(k)}\|^{2}}$
$\displaystyle\leq(1+2c_{1}\gamma)\\!~{}\|\bm{\theta}^{(k)}-\bm{\theta}^{\star}\|^{2}$
$\displaystyle\hskip
14.22636pt-2\gamma\langle{\nabla}_{\bm{\theta}}{\mathcal{}L}(\bm{\theta}^{(k)},\bm{\lambda}^{(k)})-{\nabla}_{\bm{\theta}}{\mathcal{}L}(\bm{\theta}^{\star},\bm{\lambda}^{\star}),\bm{\theta}^{(k)}-\bm{\theta}^{\star}\rangle$
$\displaystyle\hskip
14.22636pt+\gamma^{2}(1+c_{0})\|{\nabla}_{\bm{\theta}}{\mathcal{}L}(\bm{\theta}^{(k)},\bm{\lambda}^{(k)})-{\nabla}_{\bm{\theta}}{\mathcal{}L}(\bm{\theta}^{\star},\bm{\lambda}^{\star})\|^{2}$
$\displaystyle\hskip
14.22636pt+\Big{(}\frac{2\gamma}{c_{1}}+\gamma^{2}+\frac{\gamma^{2}}{c_{0}}\Big{)}\|{\bm{e}}_{\bm{\theta}}^{(k)}\|^{2}.$

Similarly, in the dual update we get,

$\displaystyle\|\bm{\lambda}^{(k+1)}-\bm{\lambda}^{\star}\|^{2}$
$\displaystyle\leq\|\bm{\lambda}^{(k)}-\bm{\lambda}^{\star}\|^{2}+\gamma^{2}\|\widehat{\bm{g}}_{\bm{\lambda}}^{(k)}-{\nabla}_{\bm{\lambda}}{\mathcal{}L}(\bm{\theta}^{\star},\bm{\lambda}^{\star})\|^{2}$
$\displaystyle\hskip
14.22636pt+2\gamma\langle\widehat{\bm{g}}_{\bm{\lambda}}^{(k)}-{\nabla}_{\bm{\lambda}}{\mathcal{}L}(\bm{\theta}^{\star},\bm{\lambda}^{\star}),\bm{\lambda}^{(k)}-\bm{\lambda}^{\star}\rangle$
$\displaystyle\leq(1+2c_{1}\gamma)\\!~{}\|\bm{\lambda}^{(k)}-\bm{\lambda}^{\star}\|^{2}$
$\displaystyle\hskip
14.22636pt+2\gamma\langle{\nabla}_{\bm{\lambda}}{\mathcal{}L}(\bm{\theta}^{(k)},\bm{\lambda}^{(k)})-{\nabla}_{\bm{\lambda}}{\mathcal{}L}(\bm{\theta}^{\star},\bm{\lambda}^{\star}),\bm{\lambda}^{(k)}-\bm{\lambda}^{\star}\rangle$
$\displaystyle\hskip
14.22636pt+\gamma^{2}(1+c_{0})\|{\nabla}_{\bm{\lambda}}{\mathcal{}L}(\bm{\theta}^{(k)},\bm{\lambda}^{(k)})-{\nabla}_{\bm{\lambda}}{\mathcal{}L}(\bm{\theta}^{\star},\bm{\lambda}^{\star})\|^{2}$
$\displaystyle\hskip
14.22636pt+\Big{(}\frac{2\gamma}{c_{1}}+\gamma^{2}+\frac{\gamma^{2}}{c_{0}}\Big{)}\|{\bm{e}}_{\bm{\lambda}}^{(k)}\|^{2}.$

Summing up the two inequalities gives:

$\begin{split}&\|{\bm{z}}^{(k+1)}-{\bm{z}}^{\star}\|^{2}\\\
&\leq(1+2c_{1}\gamma)\\!~{}\|{\bm{z}}^{(k)}-{\bm{z}}^{\star}\|^{2}+\Big{(}\frac{2\gamma}{c_{1}}+\gamma^{2}+\frac{\gamma^{2}}{c_{0}}\Big{)}E_{k}\\\
&\hskip
7.11317pt-2\gamma\langle\bm{\Phi}({\bm{z}}^{(k)})-\bm{\Phi}({\bm{z}}^{\star}),{\bm{z}}^{(k)}-{\bm{z}}^{\star}\rangle\\\
&\hskip
14.22636pt+\gamma^{2}(1+c_{0})\|\bm{\Phi}({\bm{z}}^{(k)})-\bm{\Phi}({\bm{z}}^{\star})\|^{2}\\\
&\overset{(a)}{\leq}\Big{(}1+2\gamma(c_{1}-\upsilon)+\gamma^{2}(1+c_{0})L_{\Phi}^{2}\Big{)}\|{\bm{z}}^{(k)}-{\bm{z}}^{\star}\|^{2}\\\
&\hskip
14.22636pt+\Big{(}\frac{2\gamma}{c_{1}}+\gamma^{2}+\frac{\gamma^{2}}{c_{0}}\Big{)}E_{k},\end{split}$

where (a) uses the strong monotonicity and smoothness of the map $\bm{\Phi}$.
Setting $c_{1}=\upsilon/2$ yields

$\begin{split}&\|{\bm{z}}^{(k+1)}-{\bm{z}}^{\star}\|^{2}\\\
&\leq\Big{(}1-\gamma\upsilon+\gamma^{2}(1+c_{0})L_{\Phi}^{2}\Big{)}\|{\bm{z}}^{(k)}-{\bm{z}}^{\star}\|^{2}\\\
&\hskip
14.22636pt+\Big{(}\frac{4\gamma}{\upsilon}+\gamma^{2}+\frac{\gamma^{2}}{c_{0}}\Big{)}E_{k}.\end{split}$
(39)

Observe that we can choose $\gamma$ such that
$1-\gamma\upsilon+\gamma^{2}(1+c_{0})L_{\Phi}^{2}<1$. Moreover, the above
inequality implies that $\|{\bm{z}}^{(k)}-{\bm{z}}^{\star}\|^{2}$ evaluates to

$\begin{split}&\|{\bm{z}}^{(k+1)}-{\bm{z}}^{\star}\|^{2}\\\
&\leq(1-\gamma\upsilon+\gamma^{2}(1+c_{0})L_{\Phi}^{2})^{k}\|{\bm{z}}^{(0)}-{\bm{z}}^{\star}\|^{2}+\\\
&\sum_{\ell=1}^{k}(1-\gamma\upsilon+\gamma^{2}(1+c_{0})L_{\Phi}^{2})^{k-\ell}\Big{(}\frac{4\gamma}{\upsilon}+\gamma^{2}+\frac{\gamma^{2}}{c_{0}}\Big{)}E_{\ell}\end{split}$

If $E_{k}\leq\overline{E}$ for all $k$, then ${\bm{z}}^{(k)}$ converges to a
neighborhood of ${\bm{z}}^{\star}$ of radius

$\limsup_{k\rightarrow\infty}\|{\bm{z}}^{(k)}-{\bm{z}}^{\star}\|^{2}\leq\frac{\frac{4\gamma}{\upsilon}+\gamma^{2}+\frac{\gamma^{2}}{c_{0}}}{\gamma\upsilon-\gamma^{2}(1+c_{0})L_{\Phi}^{2}}\overline{E}$
(40)

Setting $c_{0}=1$ concludes the proof.

## References

  * [1] J. Koshal, A. Nedić, and U. V. Shanbhag, “Multiuser optimization: Distributed algorithms and error analysis,” SIAM Journal on Optimization, vol. 21, no. 3, pp. 1046–1081, 2011.
  * [2] D. P. Palomar and M. Chiang, “A tutorial on decomposition methods for network utility maximization,” IEEE JSAC, vol. 24, no. 8, pp. 1439–1451, 2006\.
  * [3] S. Sundaram and C. N. Hadjicostis, “Distributed function calculation via linear iterative strategies in the presence of malicious agents,” IEEE TAC, vol. 56, no. 7, pp. 1495–1508, 2011.
  * [4] F. Pasqualetti, A. Bicchi, and F. Bullo, “Consensus computation in unreliable networks: A system theoretic approach,” IEEE TAC, vol. 57, no. 1, pp. 90–104, 2012.
  * [5] R. Gentz, S. X. Wu, H.-T. Wai, A. Scaglione, and A. Leshem, “Data injection attacks in randomized gossiping,” IEEE TSIPN, vol. 2, no. 4, pp. 523–538, 2016.
  * [6] S. Sundaram and B. Gharesifard, “Distributed optimization under adversarial nodes,” IEEE TAC, 2018.
  * [7] Y. Chen, S. Kar, and J. Moura, “Resilient distributed estimation: Sensor attacks,” IEEE Transactions on Automatic Control, 2018.
  * [8] J. Feng, H. Xu, and S. Mannor, “Distributed robust learning,” arXiv preprint arXiv:1409.5937, 2014.
  * [9] D. Yin, Y. Chen, K. Ramchandran, and P. Bartlett, “Byzantine-robust distributed learning: Towards optimal statistical rates,” arXiv preprint arXiv:1803.01498, 2018.
  * [10] D. Alistarh, Z. Allen-Zhu, and J. Li, “Byzantine stochastic gradient descent,” in NeurIPS, pp. 4618–4628, 2018.
  * [11] D. L. Donoho and P. J. Huber, “The notion of breakdown point,” A festschrift for Erich L. Lehmann, vol. 157184, 1983.
  * [12] P. J. Huber, Robust statistics.  Springer, 2011.
  * [13] S. Minsker et al., “Geometric median and robust estimation in banach spaces,” Bernoulli, vol. 21, no. 4, pp. 2308–2335, 2015.
  * [14] I. Diakonikolas, G. Kamath, D. M. Kane, J. Li, A. Moitra, and A. Stewart, “Robust estimators in high dimensions without the computational intractability,” in IEEE FOCS, pp. 655–664, 2016.
  * [15] J. Steinhardt, M. Charikar, and G. Valiant, “Resilience: A criterion for learning in the presence of arbitrary outliers,” arXiv preprint arXiv:1703.04940, 2017.
  * [16] F. P. Kelly, A. K. Maulloo, and D. K. H. Tan, “Rate control for communication networks: shadow prices, proportional fairness and stability,” Journal of the Operational Research Society, vol. 49, pp. 237–252, Mar 1998.
  * [17] S. H. Low and D. E. Lapsley, “Optimization flow control. i. basic algorithm and convergence,” IEEE/ACM Transactions on Networking, vol. 7, pp. 861–874, Dec 1999.
  * [18] P. Samadi, A. Mohsenian-Rad, R. Schober, V. W. S. Wong, and J. Jatskevich, “Optimal real-time pricing algorithm based on utility maximization for smart grid,” in IEEE SmartGridComm, pp. 415–420, Oct 2010\.
  * [19] N. Li, L. Chen, and S. H. Low, “Optimal demand response based on utility maximization in power networks,” in 2011 IEEE Power and Energy Society General Meeting, pp. 1–8, July 2011.
  * [20] M. Chiang and J. Bell, “Balancing supply and demand of bandwidth in wireless cellular networks: utility maximization over powers and rates,” in IEEE INFOCOM 2004, vol. 4, pp. 2800–2811 vol.4, March 2004.
  * [21] M. Dehghan, L. Massoulie, D. Towsley, D. Menasche, and Y. C. Tay, “A utility optimization approach to network cache design,” in IEEE INFOCOM, pp. 1–9, April 2016.
  * [22] M. Zhao, J. Li, and Y. Yang, “Joint mobile energy replenishment and data gathering in wireless rechargeable sensor networks,” in Proceedings of the 23rd International Teletraffic Congress, ITC ’11, pp. 238–245, International Teletraffic Congress, 2011.
  * [23] R. Deng, Y. Zhang, S. He, J. Chen, and X. Shen, “Maximizing network utility of rechargeable sensor networks with spatiotemporally coupled constraints,” IEEE JSAC, vol. 34, pp. 1307–1319, May 2016.
  * [24] N. Mehr, J. Lioris, R. Horowitz, and R. Pedarsani, “Joint perimeter and signal control of urban traffic via network utility maximization,” in IEEE ITSC, pp. 1–6, Oct 2017.
  * [25] J. Lavaei, D. Tse, and B. Zhang, “Geometry of power flows and optimization in distribution networks,” IEEE Transactions on Power Systems, vol. 29, no. 2, pp. 572–583, 2014.
  * [26] J. A. Taylor, Convex optimization of power systems.  Cambridge University Press, 2015.


# Benefits of non-adiabatic quantum control in quantum computation through
spin qubit systems

Nirupam Dutta<EMAIL_ADDRESS>Independent Researcher, 25B Shanti Ghosh
Street, Kolkata, West Bengal, India

(August 28, 2024)

###### Abstract

This is evident that the controllable quantum systems can be the reliable
building blocks for Quantum computation. In reality we are witnessing the
progress towards making the idea tractable enough, though optimistic but the
threshold is not very near to us. The dawn of quantum computation has begun.
In the future, we hope to see a full fledged operationally stable quantum
computer which can solve the problems beyond the scope of classical digital
computers. We may call it quantum supremacy. Nevertheless, we should not
forget that there are problems which demand classical computers to be in the
game for a better performance in comparison to the same through quantum
devices. In the current stage of computing technology, the most beneficial
area is nothing but an hybrid approach and that is for no doubt will reign the
market for the next five to ten years. This hybrid aspect has several
directions such as simulating quantum computation on a classical computer.
Keeping both the aspect, computation through real physical devices and
simulation on a classical computer by accessing available quantum computers
for cloud computing, some advantages have been discussed in this article which
will be elaborated as well in future articles. These advantages are inherent
if we can achieve proper non-adiabatic control over the spin system in the
laboratory. Otherwise these aspects can always be simulated by using quantum
algorithms to see whether they can be useful in comparison to a purely
classical computing machine. This is no doubt a new window for progress in the
direction of quantum computation.

We need controllable quantum systems in order to employ their quantum
mechanical properties to perform computational tasks zero ; one ; two  that
are either beyond the scope of any classical computer or quicker in processing
information in comparison to the same performed on a classical computing
machine three ; four . This agenda is ruling the current state of the art in
the domain of computation. A large number of theorists and experimentalists at
various organisations are participating in this race for achieving the so
called quantum supremacy five ; six ; seven ; eight . Physical systems at the
atomic and subatomic scale are quite difficult to control in order to
manipulate them according to our computational needs. For example a two level
quantum system which is the most basic text book topic of quantum mechanics
can not be prepared easily in the laboratory. It has already been four decades
that scientists are advocating in favour of quantum computation as in many
ways it can revolutionise the subject of computation and perhaps can push the
current technological horizon. Theoretical understanding has already made us
able to realise that the quantum superposition, interference and entanglement
are three important pillars nine  for making this dream a reality. In today’s
scenario, many have achieved creating qubits in various different systems
under sophisticated laboratory conditions.

Quantum computation, we all know is different from classical computation in
the following way. Instead of classical beats which are either on or off like
a switch in the classical circuit, here the basic component is a superposed
state of a two level system which we name as qubit. But unless there is
entanglement among the qubits, they are no where more efficient than classical
bits as it somehow creates the technical essence of parallel computing. With
the help of the qubits, like we do for a classical computing device, we can
create gates for performing basic operations through qubits and at the end a
measurement of this physical process though destroys the superposition but
strategically designing the circuit for constructive interference we assure
the desired results. This is the general process of gate based quantum
computation but this is not the only mode of quantum computation. There are
quantum algorithms written for collection of qubits which do not require
anything like a gate structure for computing complex problems like
optimisation problems or factorisation problems. Quantum annealing ten ;
nishi1 ; nishi2  is one of the examples which deals with adiabatic quantum
computation to solve optimisation problems. Almost decade later the D-WAVE
device king  is successfully handling various tasks and has already been
updated to a controlled device dealing with more than thousand qubits.

There are few types of quantum mechanical systems which have been successfully
employed so far to create qubits in the laboratory. These are for example
superconducting qubits, trapped ions, laser controlled photons and controlled
spin qubits. In order to achieve operationally functional computing devices
which can run through different quantum algorithms, hundreds of thousands of
stable qubits are needed. In this context it is worth mentioning that this
number is not actually close tho the number of qubits needed for achieving
quantum supremacy of the computing device. On the other hand, these large
numbers of qubits have to be entangled with each other in order to become
clusters of logical qubits to be able to perform computation. It is evident
that to perform different types of tasks based on different algorithms, the
number of active quantum qubits are different but a sustainable correlation
among them can not be assured trivially. The reason is simple and that is
nothing but the measurement sensitivity of quantum systems.

The most important thing that we still should keep in mind in this context is
that these different types of protocols and schemes that I have discussed
mostly employ the adiabatic control over the quantum systems. The basics of
quantum annealing are strictly based on the adiabatic control over the system.
The Hamiltonian of the system is strictly maintained within the slow rate of
change such that it should not break the adiabatic condition. On the other
hand, gate based quantum computation has also been developed through adiabatic
control.

In this article I will talk about some benefits of nonadiabatic controls over
quantum systems. In my demonstration I will first talk about physical spin
systems in which formation of physical qubits can be tailored by non adiabatic
evolution and in the same context, I will additionally discuss how
nonadiabatic effects can give rise to the creation of several qubits almost in
no time from parent qubits. These are important in real physical systems when
they can be controlled within a confined region of space but this challenge
obviously can be tackled by the experimentalist as far as the viability of the
process is concerned.

Furthermore, a full fledged computer is still a matter of waiting for the
sufficient number of controllable qubits. Considering the state of the art,
the recipe we present here for non adiabatic control can be useful in
simulating various schemes of quantum computation through a hybrid approach
that combines the cloud computing facilities of available quantum computers
and simulation of various methods of solving complex problems through specific
quantum algorithms.

##  I A. Qubit Multiplier

Figure 1: This is a schematic diagram of a Stern-Gerlach device which shows
that the inhomogeneous magnetic field actually separates two spin states in
the up and down direction.

Let me start with a two level system for example a physical spin qubit. A spin
qubit can always be prepared through any available mechanism. Now if we allow
the qubit $a_{0}|0\rangle+a_{1}|1\rangle$ to pass through an inhomogeneous
time dependent or independent magnetic field, the evolved state should be
determined through the nature of the evolution. We are using $0$, $1$ to
denote the up and down spin along the $z$ axis of the specified coordinate
system 1. If the system evolves through an adiabatic process the final state
remains one qubit but the state will be entangled with the wave functions in
the up and down directions. This phenomena is known as path spin entanglement
home . The final state $|\Psi(t)\rangle$ can be expressed as,

$|\Psi(t)\rangle=a_{0}e^{-i\omega_{0}t}|0\rangle\otimes\phi_{u}+a_{1}e^{i\omega_{0}t}|1\rangle\otimes\phi_{d}$
(1)

By path spin entanglement, I mean the spatial wave functions $\phi_{u}$ and
$\phi_{d}$ corresponding to upward and downward direction along the $z$ axis
are coupled to the up and down spin states. This is evident from the above
equation that spatial and spin part of the state are not product separable.
For an operationally and formally ideal situation of the Stern-Gerlach set up
1, the up and down part of the wave function can be well separated. In that
case, if an ensemble of $N$ number of spin half particles undergo this
magnetic field configuration, the up and down spins will be accumulated within
a small region of the space in the upward and downward directions
respectively. The probability of having the up spin in the spatially upward
direction is given by,

$P_{0,u}={a_{0}}^{2}\int{\phi_{u}}^{2}dX.$ (2)

Where $X$ represents spatial variables.

So the total number of accumulated particles $N_{0,u}$ in the state
$|0\rangle$ in the upward direction is $N$ times the probability $P_{0,u}$.
Similarly total number of accumulated particles $N_{0,d}$ having the down spin
$|1\rangle$ in the downward direction is $NP_{1,d}$ as the probability of
having down spin in the downward direction is given by,

$P_{1,d}={|{a_{1}}|}^{2}\int{\phi_{d}}^{2}dX$ (3)

Actually each of the $N_{0,u}$ number of particles are in the up state and at
the same time the cluster itself can be used as one logical bit $|0\rangle$.
The same argument applies to the cluster of particles which are confined in
the downward direction.

###  I.1 A. 1. Converter Oracle

This adiabatic evolution and path spin entanglement has an interesting feature
in the context of computations through the qubits. Suppose we start with a
qubit state and with one step boosting through regular ideal S-G set up, we
can create well separated two streams of two different states. Now, among the
two streams, if we allow only the upward one to be released, when considered
as a final output, the manipulation results in conversion of a superposed
state to a pure state which is equivalent to a classical bit.

Figure 2: This is a schematic diagram of the converter which can transform a
superposed state toa pure state with the help of ideal S-G apparatus.

The entire process in this case serves like an oracle which changes an input
qubit to an output which serves as a classical bit. A reverse process can
always restore to the initial state and hence can keep the unitarity
unviolated through the process. This oracle can be used as a gate that
transforms quantum to classical beats and can be useful in various algorithms
as a component of an integrated circuit to solve specific computational
problems.

###  I.2 A. 2. The non-adiabatic multiplier

A very different situation arises when the system evolves through a non-
adiabatic evolution. For this purpose, the Hamiltonian of the spin system
needs to be changing quite fast with respect to time. How fast it should be is
determined with respect to the a characteristic time scale $\tau$ of the
system which is the inverse of the energy gap between the Zeeman levels, I
mean the energy states of the two level system. If the characteristic time
$\tau$ is much less than the time needed for the unit amount change of the
Hamiltonian, the system evolves adiabatically otherwise it evolves non-
adiabatically. The precise criteria for non-adiabatic evolution is as follows
aharonov ; amin ; tong .

$\Big{|}\frac{\langle
0(t)|\dot{H}(t)|1(t)\rangle}{(E_{0}-E_{1})^{2}}\Big{|}<<1.$ (4)

$|0(t)\rangle$ and $|1(t)\rangle$ are the instantaneous eigenstates of $H(t)$
and $E_{0}$ and $E_{1}$ are the corresponding energy eigenvalues.

There can be many different possible magnetic field configurations, which can
violates this condition and serve our purpose. For a simple demonstration, let
me consider a very well known configuration which has been understood
previously in the context of spin states in other areas of physics. Here we
will be enjoying a few advantages for the purpose of quantum computing. The
considered field configuration is a rotating magnetic field and in order to
have the desired outcome, I will present a specific case when a rotating field
component will be added to a formally and operationally ideal Stern-Gerlach
setup. The rotating field is confined in the X-Y plane only and it rotates
about the Z axis. Another method of violating adiabatic condition was
presented in one of our previous articles in the context of evolution of spin
states when there is azimuthal inhomogeneity inside the S-G set up dutta . In
both the cases the mathematical construction of the problem is same. This
field strength has to be decided in such a way that the time scale of spin
states dynamics is much smaller than that of the spatial wave function. Then,
We can represent the Hamiltonian in the following way.

$\displaystyle H(t)=\gamma B(\sigma_{x}\sin\theta\cos\phi(t)$
$\displaystyle+\sigma_{y}\sin\theta\sin\phi(t)+\sigma_{z}\cos\theta)$
$\displaystyle=\omega_{0}(\sigma_{x}\sin\theta\cos\omega t$ (5)
$\displaystyle+\sigma_{y}\sin\theta\sin\omega t+\sigma_{z}\cos\theta).$

Where $\gamma$ is the gyromagnetic ratio and $\sigma_{x}$, $\sigma_{y}$,
$\sigma_{z}$ are three different Pauli spin matrices in the chosen orientation
of the total spin $S$. The rotating components makes the field rotating in the
x-y plane with a frequency $\omega$. $\theta$ is the angle between the
z-component of the magnetic moment ($\mu=\gamma S$) and the magnetic field
$\vec{B}$ and $\phi(t)$ is the azimuthal angle that magnetic field makes in
the $X-Y$plane considered here. The matrix form of the Hamiltonian is given
below.

$\displaystyle\begin{aligned}
H(t)=\frac{\omega_{0}}{2}\left(\begin{array}[]{cc}\cos\theta&e^{-i\omega
t}\sin\theta\\\ e^{i\omega t}\sin\theta&-\cos\theta\\\
\end{array}\right)\end{aligned}$

The quantity $\omega_{0}$ is the Larmor frequency of the magnetic moment of
the spin$1/2$ particle.

The instantaneous Eigenstates of the time dependent Hamiltonian can be
expressed as,

$|0(t)\rangle=\begin{pmatrix}e^{-i\omega t/2}\sin\theta/2\\\ e^{i\omega
t/2}\cos\theta/2\end{pmatrix}$

and

$|1(t)\rangle=\begin{pmatrix}e^{i\omega t/2}\sin\theta/2\\\ -e^{-i\omega
t/2}\cos\theta/2\end{pmatrix}$

.

Plugging these two states in the equation 4, we arrive to the condition for
adiabatic evolution as,

$\frac{\omega}{2\omega_{0}}\sin\theta<<1.$ (6)

This equation tells us that whether the evolution is adiabatic or not depends
on three different parameters: $\omega_{0}$, $\omega$ and $\theta$. One can
employ this freedom to break the adiabatic condition and in that case the
evolution of spin states will be determined by solving the Schrödinger
equation.

Figure 3: This is a schematic diagram of the oracle ”qubit multiplier” which
is based on the non-adiabatic control inside a formally and operationally
ideal Stern-Gerlach device.

Now let me come back to the story that we have started in this section.
Suppose the initial qubit is a superposition of two spin states $|0\rangle$
and $|1\rangle$. The moment they enter into the magnetic field, the up and
down states will be coupled to the spatial up and down component of the wave
function and will evolve according to the time dependent Hamiltonian. The
interaction Hamiltonian $H$ does not contain any spatial part hence the time
dependence will affect the spin part only. The spatial inhomogeneity of the
basic ideal S-G setup will help separating the wave functions in the upward
and downward direction along the Z axis. So we just need to study the
evolution of up and down spin states under the influence of the additional
rotating component.

We have to calculate the evolution of the up state which will be deflected in
the upward direction and for the down state which will be deflected in the
downward direction. For both purposes we have to solve the Schrödinger
equation for the time dependent Hamiltonian with two different initial
conditions with the help of the following instantaneous Eigenstates. Let’s
write down the evolved state in the upward direction in the following way,

$\Psi_{u,z}=|s_{zu}(t)\rangle\otimes\Phi_{u}.$ (7)

Where,

$|s_{zu}(t)\rangle=\alpha_{0}(t)|0(t)\rangle+\alpha_{1}(t)|1(t)\rangle$

Now, We see that the $N_{u}$ number of particles forms the cluster of qubits
which will be accumulated in the upward direction. On the other hand the
$N_{d}$ number of particles will form a cluster of qubits in the downward
direction. The state can be written as,

$\Psi_{d,z}=|s_{zd}(t)\rangle\otimes\Phi_{d}.$ (8)

Where,

$|s_{zd}(t)\rangle=\beta_{0}(t)|0(t)\rangle+\beta_{1}(t)|1(t)\rangle$

Now plugging this in the Schrödinger equation with the initial conditions
$a_{0}(t=0)=1$ for the upward stream and $a_{1}(t=0)=1$ for the down stream,
we are left with the following solutions for $\alpha_{0}$, $\alpha_{1}$,
$\beta_{0}$ and $\beta_{1}$ as follows.

$\displaystyle\begin{aligned}
\alpha_{0}(t)&=\cos\frac{\bar{\omega}t}{2}-i\frac{\omega_{0}-\omega\cos\theta}{\bar{\omega}}\sin\frac{\bar{\omega}t}{2}\\\
\alpha_{1}(t)&=i\frac{\omega\sin\theta}{\bar{\omega}}\sin\frac{\bar{\omega}t}{2}\\\
\beta_{0}(t)&=i\frac{\omega\sin\theta}{\bar{\omega}}\sin\frac{\bar{\omega}t}{2}\\\
\beta_{1}(t)&=\cos\frac{\bar{\omega}t}{2}+i\frac{\omega_{0}-\omega\cos\theta}{\bar{\omega}}\sin\frac{\bar{\omega}t}{2}\end{aligned}$
(9)

Where,
$\bar{\omega}=\sqrt{\omega_{0}^{2}+\omega^{2}-2\omega_{0}\omega\cos\theta}$

The first advantage of this nonadiabatic process is that we are able to create
two different logical qubits through two different clusters of physical
qubits. This is like an oracle which doubles the number of logical qubits just
by giving one order of non-adiabatic boosting through the time dependent
magnetic field inside the ideal S-G setup. With such $n$ number of non
adiabatic booster the logical qubits can be scaled by a factor of $2^{n}$.

The second advantage is that we can create many basic one qubit logic gates
through the non-adiabatic control. Let us discuss just a few among these gates
and the associated values of the control parameters. There are four different
parameters such as $\omega_{0}$, $\omega$, $\theta$ and $t$ which can be
adjusted to achieve the desired output. The scope of exploration is almost
unlimited but let me just show one particular choice of these free parameters
and how we can achieve three different one qubit quantum logic gates through
this.

Lets start with a particular choice which is
$\cos\theta=\frac{\omega_{0}}{\omega}$. The ratio of the two frequencies have
to be chosen in such a way so that it does not satisfy the adiabatic condition
6. It is very much possible as there is a big domain for which the ratio
$\frac{\omega_{0}}{\omega}$ can be fixed accordingly.

Case 1. Quantum Not Gate: Suppose the incoming state to the S-G setup is
$a_{0}|0\rangle+a_{1}|1\rangle$ The evolved state in the upward direction
becomes $a_{1}|0\rangle+a_{0}|1\rangle$ at some instant $t=\tau_{not}$ for
which

$\cos\frac{\bar{\omega}\tau_{not}}{2}=a_{1}\hskip 8.5359pt\text{and}\hskip
8.5359pti\sin\frac{\bar{\omega}\tau_{not}}{2}=a_{0}$ (10)

For this purpose the time dependent Hamiltonian should change up to
$H(\tau_{not})$ from $t=0$ to $t=\tau_{not}$ and should remain fixed
afterwards. The qubit at the final instant can be treated as the output of a
quantum NOT gate and can be utilised accordingly. Under this circumstance, the
qubits in the downward direction is the same as the input qubit. This
operation is equivalent to an Identity operation on the input state

Case 2. The Z Gate: With the same value of $\theta$, if we consider the final
state in the upward direction at some instant $\tau_{z}$ so that it satisfies
the following condition.

$\cos\frac{\bar{\omega}\tau_{z}}{2}=a_{0}\hskip 8.5359pt\text{and}\hskip
8.5359pti\sin\frac{\bar{\omega}\tau_{z}}{2}=-a_{1}$ (11)

In this case the final outcome in the upward direction is $a_{0}|0\rangle-
a_{1}|1\rangle$. Obviously, the Hamiltonian in this case shall remain
unchanged from $t=\tau_{z}$ with the final value $H(\tau_{z})$.

Case 3. The Hadamard Gate: For the same value of $\theta$, if we impose the
following condition,

$\cos\frac{\bar{\omega}\tau_{H}}{2}=\frac{a_{0}+a_{1}}{\sqrt{2}}\hskip
8.5359pt\text{and}\hskip
8.5359pti\sin\frac{\bar{\omega}\tau_{H}}{2}=\frac{a_{0}-a_{1}}{\sqrt{2}}$ (12)

The final state in the upward direction is
$a_{0}\frac{|0\rangle+|1\rangle}{\sqrt{2}}+a_{0}\frac{|0\rangle-|1\rangle}{\sqrt{2}}$.
This is the output of a single qubit Hadamard gate.

##  II B. Non-adiabatic control over multi qu-bit entangled states

So far, we have talked regarding the non-adiabatic control through the
rotating time dependent magnetic field inside a formally and operationally
ideal S-G apparatus. In this section we will discuss the control we can have
over multi qubit entangled states which is one of the important ingredients
for quantum computation. I won’t be making this section prolonged as in an
upcoming article I will show preparation of various many qubit quantum gates
and specific controls for introducing algorithms to solve complex problems.
Let me give an introductory description along with a few important hints that
can create some ambition for future study in the domain of non-adiabatic
quantum control.

Instead of S-G set up, here we consider Mach-zehnder setup in conjunction with
an S-G device for the purpose of creating entanglement among different multi-
particle qubits pan . This is actually done through the mechanism of swapping
of path-spin entanglement to spin spin entanglement. The issue has been well
studied previously. I am not going into the details of the mechanism as we are
trying to emphasise on the improvement that we are bringing in through the
application of an additional time dependent rotating magnetic field. Through a
specific process first the mechanism creates two particle (particle $1$ and
$2$) states which are coupled to a pseudo path like spatial variable and later
a third particle 3 is introduced which previously has not interacted with the
two particles. The swapping then creates various entanglement between 2 and 3.
Intuitively, the process is the following. Particles ‘2’ and ‘3’ are kept
separated without allowing them to interact throughout the process. Now the
scheme allows the particle ‘1’ to interact independently with the particles
‘2’ and ‘3’, without the need of having ‘2’ and ‘3’ in the near vicinity of
each other.

Now the qubit undergoes through two different S-G apparatus and these devices
provide an unitary transformation /cite. Upon measurement mentioned on the
state of the particle $1$ we can create two different kinds of entangled
states which are given below.

$1)$. $a|00\rangle_{23}+b|11\rangle_{23}$

$2)$. $a|00\rangle_{23}-b|11\rangle_{23}$

But if we use the non adiabatic control over the outgoing state of the Mach-
Zehender state through a time dependent rotating magnetic field, we can create
four different types of entangled states in the two particle tensor product
space of particle $2$ and $3$. Furthermore, by adjusting the parameter values,
we can create various two qubits gates in this context. In this article we are
not going into the details of those gates and quantum circuits that can be
employed to create various quantum algorithms for solving complex issues. This
issue I leave for a different article in near future.

##  III Discussion

My purpose in this article is to introduce the key idea of non adiabatic
control and the various benefits that it can offer for computing. Now
implementation of the control in the form of a physical device needs necessary
engineering and that might suffer from several difficulties in the laboratory.
Nevertheless the procedure opens up new possibilities for integrated devices
which can be useful in solving various complicated problems with the help of
proper algorithms. On the other hand, a mathematical simulation of non-
adiabatic control can be a good approach if it is implemented through quantum
algorithms running on basic quantum computers that are available for cloud
computing. In a different article we will show a few such algorithms and how
they can be useful in solving various computations in the simulation of many
body physics and chemistry as well as in the domain optimization problems.

## References

  * (1) R. P. Feynman, Simulating physics with computers, Int. J. Theor. Physics 21, 467-488 (1982).
  * (2) Deutsch, D., “Quantum theory, the Church-Turing principle and the universal quantum computer”, Proceedings of the Royal Society, London, vol. A400, 1985, pp. 97 – 117.
  * (3) Deutsch, D., “Quantum computational networks”, Proceedings of the Royal Society, London, vol. A425, 1989, pp. 73 – 90.
  * (4) Berthiaume, A. and Brassard, G., “The quantum challenge to structural complexity theory”, Proceedings of 7th IEEE Conference on Structure in Complexity Theory, 1992, pp. 132 – 137.
  * (5) Berthiaume, A. and Brassard, G.,  “Oracle quantum computing”, Journal of Modern Optics, vol. 41, no. 12, December 1994, pp. 2521 – 2535.
  * (6) Boixo, S. et al. “Characterizing quantum supremacy in near-term devices”, Nat. Phys. 14, 595 (2018). arXiv:1608.00263.
  * (7) Harrow, A. W. & Montanaro, A.  “Quantum computational supremacy”, Nature 549, 203 (2017).
  * (8) Aaronson, S. & Chen, L.  “Complexity-theoretic foundations of quantum supremacy experiments”, arXiv:1612.05903 (2016).
  * (9) Neill, C. et al.  “A blueprint for demonstrating quantum supremacy with superconducting qubits”, arXiv:1709.06678 (2017).
  * (10) Preskill, J.  “Quantum computing and the entanglement frontier (2012)”, 25th Solvay Conf.
  * (11) G. E. Santoro and E. Tosatti.  “Optimization using quantum mechanics: quantum annealing through adiabatic evolution”, J. Phys. A: Math. Gen., 39, R393 (2006).
  * (12) Kadowaki, Tadashi and Nishimori, Hidetoshi,  “Quantum annealing in the transverse Ising model”, PhysRevE.58. 5355, (1998).
  * (13) S. Morita and H. Nishimori.  “Mathematical foundation of quantum annealing”, J. Math. Phys., 49, 125210 (2008).
  * (14) Andrew D. King et. al.  “Coherent quantum annealing in a programmable 2000-qubit Ising chain”, Nature Physics Vol. 18, 1324-1328 (2022).
  * (15) Home D, Pan A K, Ali M M and Majumdar A S  “Aspects of nonideal Stern–Gerlach experiment and testable ramifications” Journal of Physics A: Mathematical and Theoretical 40 13975 (2007).
  * (16) Y. Aharonov, J. Anandan Phase change during a cyclic quantum evolution Phys. Rev. Lett. 58, 1593 (1987)
  * (17) M. H. S. Amin Consistency of the Adiabatic Theorem Phys. Rev. Lett. 102, 220401 (2009)
  * (18) D. M. Tong Quantitative Condition is Necessary in Guaranteeing the Validity of the Adiabatic Approximation Phys. Rev. Lett 104, 120401 (2010)
  * (19) N. Dutta and A. Dey Blurred path-spin entanglement in Stern-Gerlach apparatus: interplay between magnetic inhomogeneity and Larmor precession arXiv:1410.4396 [nucl-th] (2014)
  * (20) A. K. Pan et. al, ”Swapping path-spin intraparticle entanglement onto spin-spin interparticle entanglement”, Euro. Phys. Lett. 89, 10005(2010).


Laser-induced alignment of nanoparticles and macromolecules for single-particle-imaging
Muhamed Amin
Jean-Michel Hartmann
Amit K. Samanta
Jochen Küpper

Laser-induced alignment of particles and molecules was long envisioned to support
three-dimensional structure determination using single-particle imaging with x-ray free-electron
lasers [PRL 92, 198102 (2004)]. However, geometric alignment of isolated macromolecules has not
yet been demonstrated. Using molecular modeling, we analyzed and demonstrated how the alignment
of large nanorods and proteins is possible with standard laser technology, and performed a
comprehensive analysis on the dependence of the degree of alignment on molecular properties and
experimental details. Calculations of the polarizability anisotropy of about 150,000 proteins
yielded a skew-normal distribution with a location of 1.2, which reveals that most of these
proteins can be aligned using appropriate, realistic experimental parameters. Moreover, we
explored the dependence of the degree of alignment on experimental parameters such as particle
temperature and laser-pulse energy.

§ INTRODUCTION

X-ray free-electron lasers (XFELs) promise the diffractive imaging of single molecules and
nanoparticles at atomic resolution <cit.>. Ultra-short, high intensity x-ray
pulses interact with individual molecules and using the “diffraction-before-destruction”
approach <cit.> large series of diffraction patterns are
collected. Two-dimensional diffraction patterns from randomly oriented samples would then be
computationally assembled to a diffraction volume to retrieve the three-dimensional
structure <cit.>.

In the standard single-particle imaging (SPI) method, the camera images do not contain a
priori information of the molecules' orientation and, so far, the uncertainty is attacked
in silico <cit.>. Significant efforts were made in improving the
reconstruction process and the achievable resolution, but it is still a highly challenging task,
especially for weakly scattering particles like individual proteins, where diffraction signals from
single molecules are generally not strong enough to allow for the averaging and sorting
approaches <cit.>. This is one of the major bottlenecks for
atomic-spatial-resolution SPI.

Imaging samples with controlled alignment or orientation <cit.> allows to
significantly mitigate this problem by allowing to sum the diffraction signals from many
identically aligned molecules to provide a much stronger signal, thus improving the image
reconstruction step and paving the way toward atomic-resolution SPI <cit.>. Careful analysis of simulated diffraction patterns of
laser-aligned proteins demonstrated that it is possible to observe the secondary structure of
proteins with only reasonably-strong degrees of alignment
$\costhreeD\geq0.9$ <cit.>.

The alignment of small molecules using external electric fields from moderately intense, nonresonant
light pulses was studied extensively using experimental and theoretical
methods <cit.>. Strong alignment was achieved for linear, symmetric
top, and asymmetric top molecules in the adiabatic <cit.>,
intermediate <cit.>, and impulsive
regimes <cit.>. Considerable efforts were made
for laser-induced alignment of large molecules <cit.> and even for some complex and floppy polyatomic
molecules <cit.> as well as for weakly bound molecular
complexes <cit.>. Such aligned-molecule samples were also studied by
electron <cit.> and x-ray diffractive imaging
experiments <cit.>. Possibilities to laser-align large biomolecules without
deterioration of the secondary structure were proposed <cit.>,
but, so far, no alignment for such systems has been reported.

For macromolecules, there are several challenges for achieving the required alignment. Theoretical
predictions supported by in silico analysis can be a very important step to guide the
experiments. However, atomistic molecular dynamics (MD) simulations are computationally very expensive
for large particles, especially when the interaction time extends to hundreds of
nanoseconds <cit.>. In addition, to accurately predict the “ensemble averaged”
single-particle diffraction images of macromolecules, simulation of a large distribution of
particles is essential rather than a single particle. Ensemble computations are also crucial for
studying the temperature effects, as this can be an important factor in preserving the secondary
structures of macromolecules, like proteins, and this step makes MD simulations even more
computationally expensive.

In this work, we predicted and analyzed the laser-induced alignment of
(bio-)nanoparticles. Nanoparticles and proteins were treated as rigid bodies, supported by previous
molecular-dynamics simulations of structural changes in strong electric
fields <cit.>. The key parameters for such simulations are the overall
polarizability, shape and their anisotropy, the temperature of the macromolecules, and the alignment
laser field. We disentangled how these parameters can be tuned for maximum alignment of gold-nanorod
model systems and how this can be exploited for the strong laser alignment of biological
macromolecules, , proteins.

§ COMPUTATIONAL METHODS

The response of the particles to a nonresonant external electric field was calculated based on their
polarizability tensors, which yielded the time-dependent induced dipole moments. For metallic
nanorods, the polarizabilities were directly obtained by solving Laplace’s equation with Dirichlet
boundary conditions and using Monte Carlo path integral methods <cit.>. For
proteins, the polarizabilities were derived from the same calculations using regression-based
scaling <cit.>.

Molecular ensembles were set up with random initial phase-space distributions, , with initial
angular velocities according to a Boltzmann distribution at the given temperature and with random
initial angular positions. Each run included 20,000 particles. The angular positions of particles
were stored in quaternions throughout the calculation. The inertial tensors of the artificial
nanorods were calculated for cylinders. However, for proteins the tensors were calculated based on
the atomic masses and coordinates as available from the protein data bank
(PDB) <cit.>.

Electric fields of the alignment-laser pulses were represented by Gaussian functions with variable
peak intensities. For simplicity, we used a temporal full-width at half maximum (FWHM) of 8 ns,
corresponding to the typical pulse duration of standard Q-switched Nd:YAG lasers, often used in
alignment experiments. We also assumed typical linearly polarized laser pulses with intensities of
$10^{10}\ldots10^{12}$ .

The particles' angular phase-space positions were propagated in time by integrating Euler's
equations of angular motion using a new, Python-based, openly available software package
CMIclassirot <cit.>, which was based on and checked against previous
classical-alignment computations <cit.>. All simulations propagated the
particles for 50 ns, except for the alignment of gold nanoparticles in Fig. <ref>, for
which the simulation time was extended to 200 ns. The effect of resonances was ignored in the
simulation, , the laser-field frequency was assumed to be far off resonance from the molecules
and nanoparticles absorption.

§ RESULTS

The computed time-dependent degrees of laser-induced alignment of nanorods of three different sizes
are shown in Fig. <ref> for 298 K and 4 K and alignment-laser peak intensities of
$10^{11}\ldots10^{12}~\Wpcmcm$. Throughout this manuscript, nanorod sizes are represented as
($height$/nm, $diameter$/nm), with the calculations performed for 102, 5010,
10020, respectively.

Degree of alignment of gold nanorods for different temperatures ($T$) and laser-field
intensities (I). a) $T=298$ K and $I=10^{11}~\Wpcmcm$. b) $T=298$ K and $I= 10^{12}~\Wpcmcm$.
and c) $T=4$ K and $I=10^{11}~\Wpcmcm$. The degree of alignment of 102 rods is
shown in red, 5010 in blue, and 10020 in black. The temporal laser
profile is indicated by the shaded gray area and the corresponding field intensities are
specified on the secondary axes. Please see the text for details and definitions of the

During the pulse and with particles initially at room temperature, the degree of alignment for all three
particles' sizes, Fig.<ref> a, b, follows the temporal laser profile on the rising
edge in a quasi-adiabatic fashion <cit.>. The smallest particle
102 exhibits the largest alignment of $\costhreeD=0.96$.

However, the laser turn-off dynamics show considerable differences: The two larger nanorods exhibit
permanent alignment after the laser pulse is turned off, whereas the small nanorod
quasi-adiabatically follows the temporal laser profile to an isotropic field-free angular
distribution, , $\costhreeD=1/3$. This can be rationalized by comparing the rotation periods of
the nanorods, , the average time needed for each nanorod to rotate around its center of mass by
360, to the alignment-laser-pulse duration. The rotation periods depend on the temperature,
because the initial phase space is assigned based on a corresponding Boltzmann distribution. For
large nanorods at room temperature, the rotational periods are within 10 (, $\omega\approx6.3\cdot10^{5}$ rad/s), which is three orders of
magnitude larger than the pulse duration. Thus, these particles exhibit non-adiabatic alignment,
resulting in a permanent alignment after the pulse is off, Fig. <ref> a, b. Already
at low intensity, $I=10^{11}~\Wpcmcm$, the particles are confined even after the laser is off to
fully rotate in a plane containing the laser polarization vector, corresponding to $\costhreeD=0.5$,
Fig. <ref> a. Thus, no increase in the permanent alignment is possible at
$I=10^{12}~\Wpcmcm$, Fig. <ref> b. For small nanorod, the rotational period reduces to
20 ns, which is comparable to the alignment-laser pulse duration. Thus, these particles exhibit
quasi-adiabatic alignment <cit.> and no permanent
alignment is observed after the pulse, Fig. <ref> a, b.

However, at 4 K, which is experimentally achievable <cit.>, the rotational
period increases to 400 ns for the 102 rods, resulting in the transition to the
non-adiabatic regime and thus a field-free permanent alignment is observed even for this smallest
rod after the pulse is off, see Fig. <ref> c.

For the larger nanorods (blue, black) the degree of permanent alignment after the pulse shows a
strong oscillation that decays with time, which is a result of the slow dephasing between the
rotations of the confined nanorods, which start to rotate with different angular velocities, but
significantly slower than for the small particles. This variation in the angular velocities causes the
decay of the pattern with time. Furthermore, to ensure that these oscillations were not a result of
undersampling, we performed a convergence study in which we increased the number of particles for
the large nanorods to 200,000 particles and we obtained the same oscillations.

Although all our simulations take into account the full polarizability tensor, it is instructive to
study the degree of alignment as a function of the nanorod shape described by the ratio of its
principal moments of polarizability, , the ratio of $\alpha_{\parallel}$ to $\alpha_{\perp}$,
which correlates with the polarizability anisotropy:

\begin{align}
\label{eqn:pol-anisotropy}
\alpha_\parallel &= \max\left(\alpha_{11},\alpha_{22},\alpha_{33}\right) \notag \\
\alpha_\perp &= \left(3\alpha-\alpha_\parallel\right)/2 \\
\alpha_r &= \alpha_\parallel / \alpha_\perp \notag
\end{align}

with the static-polarizability components in the principal axes of polarizability frame
$\alpha_{ii}$ and $\alpha=\sum\alpha_{ii}$. For $\alpha_r$ the scaling factor applied to the
elements of the polarizability tensors to account for the dielectric medium of proteins cancels out,
and does not have to be included in this discussion of shape.

a) Maximum degree of alignment, max, as a function of the polarizability ratio
$\alpha_r$ at different temperatures and intensities for a nonorod of $10\pi~\text{nm}^3$. At
low temperature a high degree of alignment is achieved even at very low intensity. b) The
$\alpha_r$ distribution of about 150,000 proteins in the protein data bank (PDB) approximately
follows a skew-normal distribution with a location of $1.2$, scale of $0.7$ and skewness
parameter of $6.5$.

Fig. <ref> a shows the dependence of the maximum degree of alignment , maxon
$\alpha_r$ for two different laser peak intensities of $I= 10^{10}~\Wpcmcm$ and $I= 10^{11}~\Wpcmcm$.
There is a quick rise of maxfor $\alpha_r$ in the range
$1.2\ldots2.5\cdot10^{10}~\Wpcmcm$, depending on temperature, confirming that with increasing
particle anisotropy there is a significant increase in the maximum degree of alignment, as could be
expected. Further increasing $\alpha_r$ led to a further slow increase in the maximum alignment
toward an asymptotic maximum. Furthermore, for higher intensity one observes stronger alignment,
especially for small values of $\alpha_r$. The same holds for lower temperatures, which enable
significantly increased alignment even at lower laser intensities <cit.>. Specifically, for $10^{10}~\Wpcmcm$ and room-temperature particles it
requires $\alpha_r>8$ to obtain a reasonable degree of alignment of $\costhreeD>0.8$, whereas at low
temperature (4 K) <cit.>, strong alignment of $\costhreeD=0.9$ is achieved
even at $\alpha_r=1.5$ .

Proteins and similar biological macromolecules are the principal target for single particle imaging.
Thus, we calculated the polarizability ratio $\alpha_r$ of about 150,000 proteins using the database
we built previously <cit.>, see Fig. <ref> b. These data were
obtained by calculating the polarizability tensors of proteins based on their PDB structures using
ZENO <cit.>. These tensors were diagonalized to put all
molecules in the polarizability frame, their $\alpha_r$ values were computed according to
(<ref>), and summarized in the histogram in Fig. <ref> b. The $\alpha_r$
distribution follows a skew-normal distribution with a location of $\ordsim1.2$. As $96~\%$ of the
data have an $\alpha_r$ larger than the location value, this distribution directly indicates that a
corresponding, significant fraction of these proteins have sufficiently anisotropic polarizabilities
that make them directly amenable to strong laser alignment.

To provide further understanding of the effect of shape of the object on the maximum achieved
alignment, we studied the effect of laser intensity and particle temperature for a gold nanorod
with $\alpha_r=1.2$ and size 3.63.2, see Fig. <ref>.

a) The effect of laser pulse intensity on the degree of alignment of a gold nanorod with
$\alpha_r=1.2$ and size 3.63.2 at 298 K (blue) and 4 K (orange). A good degree of
alignment $\max\costhreeD=0.9$ is achieved at 4 K with intensities as low as of
$10^{10}~\Wpcmcm$ (filled circles). However, a significantly stronger intensities (more than
$5\cdot10^{11}~\Wpcmcm$) are required to obtain a good degree of alignment at 298 K. b) The
effect of temperature on the degree of alignment of the same gold nanorod (polarizability ratio
of 1.2) for two different laser intensities, $10^{10}~\Wpcmcm$ (blue) and $10^{11}~\Wpcmcm$
(orange). Strong alignment is achieved for both intensities at very low temperatures. However,
a quick decay in the maxis observed for the low intensity with increasing the

Fig. <ref> a shows the effect of different laser intensities on the degree of alignment at
298 K and 4 K. The degree of alignment quickly increases at low intensities, depending on
temperature, and saturates toward the asymptotic limit (vide supra). Again, the achievable
degree of alignment at 4 K is significantly higher than at room temperature. In the latter case,
strong alignment of $\max\costhreeD=0.9$ is achieved for the highest intensities of
$10^{12}~\Wpcmcm$, whereas a 4 K sample enables very strong alignment of $\max\costhreeD=0.95$ for
intensities as low as $5\cdot10^{10}~\Wpcmcm$.

Similarly, Fig. <ref> b shows the degree of alignment for the same particle versus the
temperature at two different laser intensities, demonstrating the decrease of alignment with
increasing temperature, especially so for weaker laser intensities.

Thus, based on our model system, which is a nanorod with an anisotropy similar to that of most of
the proteins,  Fig. <ref>, a good degree of alignment is achievable for most proteins at 298 K
using 10 ns alignment-laser pulses with peak intensities around $10^{11}~\Wpcmcm$. This alignment
can be improved significantly exploiting cryogenically cooled (4 K) <cit.>
proteins even with a much weaker laser pulse. Furthermore, the cooling will also reduce the chance
of structural damage of proteins by intense laser fields <cit.>.

Degree of alignment of the green fluorescent protein, depicted in the inset, at 298 K
(red) and 4 K (blue) using $10^{12}~\Wpcmcm$ pulse intensity (shaded area in gray). In
consistency with Fig. <ref>, a good degree of alignment of $\costhreeD=0.84$ is
achieved at 4 K even at very low intensities, , in the raising flank of the laser pulse at
$5\cdot10^{10}~\Wpcmcm$. However, a stronger intensity of $10^{12}~\Wpcmcm$ is required to
achieve similar alignment at room temperature.

To simulate the alignment of an actual protein, we have applied our method to simulate laser-induced
alignment of the prototypical green fluorescent protein (GFP), which has a cylindrical shape
($\alpha_r=1.5$) and thus strong alignment could be expected, see Fig. <ref>. We calculated the
inertia and polarizability tensors based on the PDB structure with PDBID
1GFL <cit.>. The elements of the polarizability tensors were scaled by a factor
of $0.4$ to account for the dielectric medium of proteins
($\epsilon_r=3.2$ <cit.>). Then, we applied a Gaussian laser pulse with a
FWHM of 8 ns. As shown in Fig. <ref>, moderate alignment of $\max\costhreeD=0.84$ was
obtained at room temperature for a peak intensity of $10^{12}~\Wpcmcm$. However, at 4 K
$\max\costhreeD=0.85$ was already achieved for $2\cdot10^{10}~\Wpcmcm$ and for $10^{12}~\Wpcmcm$ we
obtained $\max\costhreeD=0.94$.

In summary, using simulations based on the classical dynamics of rigid bodies we showed that
significant laser-induced alignment of nanorods and biological macromolecules, with typical
polarizability ratios $\alpha_{\parallel}/\alpha_{\perp}$, can be achieved. The dependence of the
degree of alignment on the alignment-laser intensity, sample temperature, and molecular size and
polarizability were analyzed. We showed that a very high degree of alignment can be achieved for
cryogenically-cooled <cit.> proteins at a moderate laser power of
$10^{10}~\Wpcmcm$, which should not cause radiation damage <cit.>. This high
degree of control of $\costhreeD\geq0.94$ paves the way for future atomic
resolution <cit.> single particle x-ray and electron diffractive imaging
experiments. Furthermore, the achievable atomic spatial and femtosecond temporal resolution provide
the prerequisites for future time resolved studies of ultrafast biochemical dynamics.

Our approach provides clear insight into the optical control of macromolecules and enables better
modeling of the experimental parameters for successful laser-induced alignment experiments, which are
currently in progress in our group. Envisioned future experiments plan to make use of our cryogenic
nanoparticle cooling setup <cit.> together with efficient laser control to
achieve a very high degree of alignment for shock-frozen proteins and, in turn, sub-nanometer
resolution in single particle x-ray imaging. We believe that this framework will prove useful for
furthering the field of single-particle x-ray imaging and would allow us to observe atomically
resolved snapshots of ultrafast chemical dynamics.

The achievable strong laser alignment of nanoscopic objects could have further applications in
nanoscience <cit.> as well as in nanoscale quantum optics or quantum
sensing <cit.>.

§ ACKNOWLEDGMENT

We acknowledge financial support by Deutsches Elektronen-Synchrotron DESY, a member of the Helmholtz
Association (HGF) and the use of the Maxwell computational resources operated at DESY. This work has
been supported by the European Research Council under the European Union’s Seventh Framework Program
(FP7/2007-2013) through the Consolidator Grant COMOTION (614507) and the Cluster of Excellence
“Advanced Imaging of Matter” (AIM, EXC 2056, ID 390715994) of the Deutsche Forschungsgemeinschaft

§ CODE AVAILABILITY

Classical rotational dynamics simulations were performed using CMIclassrot, available at