Join-Up-To(m): improved hyperscalable load balancing

Abstract

Various load balancing policies are known to achieve vanishing waiting times in the large-scale limit, that is, when the number of servers tends to infinity. These policies either require a communication overhead of one message per job or require job size information. Load balancing policies with an overhead below one message per job are called hyperscalable policies. While these policies often have bounded queue length in the large-scale limit and work well when the overhead is somewhat below one, they show poor performance when the communication overhead becomes small, that is, the mean response time tends to infinity when the overhead tends to zero even at low loads. In this paper, we introduce a hyperscalable load balancing policy, called Join-Up-To(m), that remains effective even when the communication overhead tends to zero. To study its performance under general job size distributions, we make use of the “queue at the cavity" approach. We provide explicit results for the first two moments of the response time, the generating function of the queue length distribution and the Laplace transform of the response time. These results show that the mean response time only depends on the first two moments of the job size distribution.

Appendices

A On the convergence to the queue at the cavity

In this appendix, we make a number of observations related to the convergence to the queue at the cavity for exponential job sizes and finite buffers of size B. The results in this section can also be generalized to phase-type distributed job sizes (which complicates notations). In (1*), we show that the stochastic system consisting of N servers is a density-dependent population process as defined by Kurtz [17]. In (2*), we present an expression for the drift function, which is not continuous everywhere, and define a differential inclusions based on the drift function. Leveraging the framework in Gast and Gaujal [18] allows us to show that the sample paths of the stochastic systems converge to the set of solutions of the differential inclusion over finite time scales as N tends to infinity. If the differential inclusion has multiple solutions, the system may converge to any solution of the DI, depending on its random innovations. In (3*), we argue that there exists a solution of the differential inclusion that makes a so-called sliding motion in a certain region of the state space, and this region contains a fixed point that corresponds to our queue at the cavity.

Given the above three results, weak convergence of the steady-state measures to the Dirac measure of the fixed point in (3*) follows due to Gast and Gaujal [18, Section 4.2] provided that we can show that the trajectory of any solution of the differential inclusion converges to this fixed point. In other words, it suffices to show that the fixed point mentioned in (3*) is a global attractor of any solution to the differential inclusion in order to exchange the limits of t and N. This proof of global attraction is still missing. We comment on a possible approach at the end of this section.

(1*) We first show that the stochastic system consisting of N servers with exponential service times and finite buffers of size B is a density-dependent population process. Define the variables $Y^{(N)}_{i,j}(t)$, for $0\le j\le i \le B$ as the fraction of the N servers that have queue length j and for which the dispatcher has an estimated queue length equal to $i \ge j$ at time t. Clearly, due to the exponential job sizes the variables $Y^{(N)}_{i,j}(t)$ form a continuous time Markov chain on the state space $\mathcal {S}^{(N)} =$ $\{y_{i,j} \mid 0\le j\le i \le B, \sum _{i,j} y_{i,j}=1, Ny_{i,j} \in \mathbb {N}\}$ $\subseteq \mathbb {Z}^{(B+1)(B+2)/2}/N$. This Markov chain is a density-dependent population process if there exists a finite set $\mathcal {L} \subset \mathbb {Z}^{(B+1)(B+2)/2}$ (with $0 \not \in \mathcal {L}$), such that for each $\ell \in \mathcal {L}$ and $y \in S^{(N)}$, the rate of transition from y to $y + \ell /N$ is of the form $N \beta _\ell (y) \ge 0$, where $\beta _{\ell }(\cdot )$ does not depend on N. Let $e_{(i,j)} \in S^{(N)}$ be the vector with $y_{i,j}=1$ (and zeros elsewhere). For the JUT(m) system, we have three types of transitions. (1) We can have an arrival that is assigned to a queue with length j and estimated length i. These transitions are denoted as $\ell _{a(i,j)} = -e_{(i,j)}+e_{(i+1,j+1)}$ (for $i < B$) as they change the queue state from (i, j) to $(i+1,j+1)$. Let $\kappa (y)$ be the minimum of m and the smallest estimated queue length when the system is in state y, that is, $\kappa (y) = \min (m, \min \{i \mid \exists j: y_{i,j} > 0 \})$. As the job arrivals occur at rate $\lambda N$ and a job is assigned to a queue with the smallest estimated queue length if $\kappa (y) < m$ and at random otherwise, we have

$$\begin{aligned} \beta _{\ell _{a(i,j)}}(y) = \left\{ \begin{array}{ll} 0 &{} i \not = \kappa (y)< m\\ \lambda y_{i,j} / \sum _{s} y_{i,s} &{} i = \kappa (y) < m\\ \lambda y_{i,j} &{} \kappa (y) = m \end{array} \right. \end{aligned}$$

(2) A service completion can occur in a server with length j and estimated queue length i. We denote these transitions as $\ell _{s(i,j)} = -e_{(i,j)} + e_{(i,j-1)}$ for $i \ge j > 0$. As service completions do not depend on other queues, we have $\beta _{\ell _{s(i,j)}}(y) = y_{i,j}$ due to the exponential service times with mean 1. (3) The last type of transition that can occur is an update from an idle server, which changes its state from (i, 0) to (0, 0) for $i > 0$. We denote these as $\ell _{u(i,0)} = -e_{(i,0)}+e_{(0,0)}$. As such updates occur at rate $\delta _0$ in any idle queue, we have $\beta _{\ell _{u(i,0)}}(y)=\delta _0 y_{i,0}$. The functions $\beta _\ell (\cdot )$ do not depend on N, therefore the Markov chain is a density-dependent population process.

(2*) The drift function f(y), with components $f_{(i,j)}(y)$ in our case, of a density-dependent population process are defined as $f(y) = \sum _{\ell \in \mathcal {L}} \beta _{\ell }(y) \ell$. Let $u_i = \sum _{j=0}^i y_{i,j}$. Given the above discussion on the transitions in $\mathcal {L}$, we have

$$\begin{aligned} f_{(i,j)}&(y) = -1[j>0] y_{i,j}+1[i> j]y_{i,j+1} \nonumber \\&- 1[0=j<i] \delta _0 y_{i,0} +\delta _0 1[i=j=0] \sum _{s>0}y_{s,0} \nonumber \\ {}&- 1[i< m] \frac{\lambda }{u_i} y_{i,j} 1[\kappa (y)=i] + 1[0 < j \le i \le m] \frac{\lambda }{u_{i-1}}y_{i-1,j-1} 1[\kappa (y)=i-1] \nonumber \\&-1[i \ge m] \lambda y_{i,j} 1[\kappa (y) = m] +1[i> m, j > 0] \lambda y_{i-1,j-1} 1[\kappa (y) = m], \end{aligned}$$

(19)

where $1[A]=1$ if A is true and $1[A]=0$ otherwise. The first two terms are due to the service completions, the next two due to the updates and the remaining ones are a result of the arrivals. Note that the $1[i \ge m]$ and $1[i > m]$ conditions on the last two terms can be dropped as $y_{i,j} = 0$ for $i < m$ when $\kappa (y) = m$. Further, for ease of presentation, the changes needed due to having a finite B are omitted.

When the drift function f(y) is Lipschitz continuous Kurtz showed that the sample paths of the stochastic system converge to the solution of the set of ODEs given by $dy(t)/dt = f(y(t))$ over any finite time interval [0, T]. In our case, the drift function f is clearly not continuous due to the presence of the $\kappa (\cdot )$ function. The result of Kurtz was however generalized in Gast and Gaujal [18, Theorem 5] to systems with drifts that contain discontinuities. More specifically, define the differential inclusion (DI) $dy(t)/dt \in F(y(t))$ with $y(0)=y_0$ where F(y) is the convex closure of the set of all f(y) values that can be obtained as $f(y) = \lim _n f(y_n)$ with $\lim _n y_n = y$. Let $\mathcal {G}_T(y_0)$ be the set of solutions to the DI on [0, T] with $y(0)=y_0$, where a solution is an absolutely continuous function y such that $df(y)/dt \in F(y(t))$ almost everywhere. Gast and Gaujal [18, Theorem 5] then implies that

$$\begin{aligned} \inf _{y \in \mathcal {G}_T(y_0)} \sup _{t \in [0,T]} \Vert Y^{(N)}(t) - y(t) \Vert \rightarrow 0,\end{aligned}$$

in probability provided that $\sup _y \sum _{\ell \in \mathcal {L}} \beta _\ell (y) < \infty$ and $\sum _{\ell \in \mathcal {L}} \Vert \ell \Vert \sup _y \beta _\ell (y) < \infty$. Both conditions hold in our case as $\mathcal {L}$ is finite and $\sup _y \beta _\ell (y) \le \max (1,\delta _0)$.

To define the set-valued function F(y), we introduce the vectors $w^k(y)$ for $k=0,\ldots , \kappa (y)-1$ with (i, j)-th component given by:

$$\begin{aligned} w^k_{(i,j)}(y)&= -1[j>0] y_{i,j} + 1[i>j]y_{i,j+1} \nonumber \\&- 1[i>j=0] \delta _0 y_{i,0} +\delta _0 1[i=j=0] \sum _{s>0}y_{s,0} \nonumber \\&- \lambda 1[k=i] + \lambda 1[j>0, k=i-1]. \end{aligned}$$

(20)

Looking at (19) one finds that the set F(y) is defined as the convex closure of the set $\{w^0(y),\ldots ,w^{\kappa (y)-1}(y),f(y)\}$. When $\kappa (y)=m$ this means that F(y) contains all functions ${\tilde{f}}(y)$, with components ${\tilde{f}}_{(i,j)}(y)$ of the form

$$\begin{aligned} {\tilde{f}}_{(i,j)}(y) =&-1[j>0] y_{i,j}+1[i> j]y_{i,j+1} \nonumber \\&- 1[i>j=0] \delta _0 y_{i,0} +\delta _0 1[i=j=0] \sum _{s>0}y_{s,0} \nonumber \\ {}&- 1[i< m] \lambda \alpha _i + 1[0 < j \le i \le m] \lambda \alpha _{i-1} \nonumber \\&- \lambda y_{i,j} \alpha _m +1[j > 0]\lambda y_{i-1,j-1} \alpha _m, \end{aligned}$$

(21)

with $\alpha _{i} \in [0,1]$ and $\sum _{i=0}^{m} \alpha _{i} =1$.

(3*) Suppose now that we are in a region of the state space where $\kappa (y)= m$ and $q_0(y)=\sum _{i\ge 0} y_{i,0} \le 1-\lambda$. In order to remain in this part of the state space by making a so-called sliding motion, the drift of $y_{0,0}$ should be zero, such that $y_{0,0}$ remains zero. By (21) and the fact that $y_{i,j} = 0$ for $i <m$ when $\kappa (y)=m$ shows that $\alpha _0 = \delta _0 q_0(y)/\lambda$. Further, if we demand that $y_{i,j}$ remains zero for $0< i < m$, then (21) indicates that $\alpha _{i}=\alpha _{i-1}$. As the sum of all $\alpha$’s equals one, we have $\alpha _m = 1-\delta _0 q_0(y) m/\lambda \le 1$ when $q_0(y)\le 1-\lambda$ due to our assumption throughout the paper that $\lambda > \delta m$ and the fact that $\delta _0 = \delta /(1-\lambda )$. If we now focus on the region with $q_0(y)=1-\lambda$ during this sliding motion, we find that $\alpha _{i} = \delta /\lambda$ for $0< i < m$ and $\alpha _m = 1-\delta m/\lambda = {\tilde{\lambda }}/\lambda$. When we plug in the above $\alpha$ values in (21), we find

$$\begin{aligned} {\tilde{f}}_{(i,j)}(y) =&-1[j>0] y_{i,j}+1[i> j]y_{i,j+1} \nonumber \\&\ \ - 1[i>j=0] \delta _0 y_{i,0} +\delta _0 1[i=j=0] \sum _{s>0}y_{s,0} \nonumber \\ {}&\ \ - 1[i=j< m] \delta + 1[0 < j = i \le m] \delta \nonumber \\&\ \ - {\tilde{\lambda }} y_{i,j} +1[j> 0]{\tilde{\lambda }} y_{i-1,j-1}, \nonumber \\ =&-1[j>0] y_{i,j} +1[i> j] y_{i,j+1} \nonumber \\&\ \ - 1[i>j=0] \delta _0 y_{i,0} +1[i=j=m] \delta _0 (1-\lambda ) \nonumber \\&\ \ - {\tilde{\lambda }} y_{i,j}+1[j > 0]{\tilde{\lambda }} y_{i-1,j-1}, \end{aligned}$$

(22)

where the second equality is due to $\sum _{s >0}y_{s,0} = \sum _{s \ge m}y_{s,0} = 1-\lambda = \delta /\delta _0$ when $\kappa (y)= m$.

If we sum these drifts over i and use the fact that $\sum _{i>0} y_{i,0}=1-\lambda$, we find

$$\begin{aligned} \sum _{i \ge j} {\tilde{f}}_{(i,j)}(y) =&-1[j>0] y_{j} + y_{j+1} - 1[j=0]\delta _0 (1-\lambda ) \nonumber \\&\ \ +1[j=m] \delta _0 (1-\lambda ) - {\tilde{\lambda }} y_{j}+1[j > 0]{\tilde{\lambda }} y_{j-1}, \end{aligned}$$

(23)

where $y_j = \sum _{i\ge j} y_{i,j}$. Recall now that the queue at the cavity for the JUT(m) policy with exponential job sizes is defined as an M/M/1 queue with arrival rate ${\tilde{\lambda }} = \lambda - \delta m$, except that when the queue is empty, there are also batch arrivals of size m that occur at rate $\delta _0 = \delta /(1-\lambda )$ such that the probability that the queue is idle is given by $1-\lambda$. By demanding that $\sum _{i\ge j} {\tilde{f}}_{(i,j)}(y) = 0$ and by replacing $1-\lambda$ by $y_0$, we obtain the balance equations of such an M/M/1 queue given by

$$\begin{aligned} y_0 (\delta _0 + {\tilde{\lambda }})&= y_1, \\ y_m (1+{\tilde{\lambda }})&= {\tilde{\lambda }} y_{m-1} + y_{m+1} + y_0 \delta _0 \\ y_j (1+{\tilde{\lambda }})&= y_{j+1} + {\tilde{\lambda }} y_{j-1}, \\ \end{aligned}$$

for $m \not = j > 0$.

This completes items (1*) to (3*). To prove convergence of the stationary measures, we must show that the fixed point of (3*) is a global attractor for any solution of the differential inclusion. This could be done by first showing that there is a unique solution and subsequently showing that any trajectory of this solution (for any starting point $y_0$, including all points with $\kappa (y_0) < m$) converge to this fixed point. A sufficient condition such that we have at most one solution is that the set valued function F(y) is one-sided Lipschitz. This means that for any $y, y' \in \mathbb {R}^{(B+1)(B+2)/2}$ and any $z \in F(y), z' \in F(y')$, we have

$$\begin{aligned}\langle y-y',z-z' \rangle \le L \Vert y-y' \Vert ^2, \end{aligned}$$

for some constant L, where $\langle x,y \rangle$ is the inner product. The following example indicates that the set-valued function F(y) that characterizes our differential inclusion is not one-sided Lipschitz. Let $m=1$ and let y be such that $y_{0,0}=\epsilon , y_{1,1}=1-\epsilon$ which implies that the only nonzero f(y) components are $f_{0,0}(y)=-\lambda , f_{1,0}(y)=1-\epsilon$ and $f_{1,1}(y)=\lambda +\epsilon -1$. Let $y'$ be such that $y'_{1,1}=1-2\epsilon$ and $y'_{2,1}=2\epsilon$, then the nonzero components of $f(y')$ are given by $f_{1,0}(y')=1-2\epsilon , f_{1,1}(y')=-1+2\epsilon -\lambda$, $f_{2,2}(y')=\lambda$, $f_{(2,0)}(y')=2\epsilon$ and $f_{(2,1)}(y')=-2\epsilon$. As $f(y) \in F(y)$ and $\Vert y-y' \Vert ^2 = O(\epsilon ^2)$, we must have that $\langle y-y',f(y)-f(y') \rangle = O(\epsilon ^2)$. However,

$$\begin{aligned} \langle y-y',f(y)-f(y') \rangle&= \epsilon (f_{(0,0)}(y)-f_{(0,0)}(y')) + \epsilon (f_{(1,1)}(y)-f_{(1,1)}(y'))\\&\ \ \ \ \ -2\epsilon (f_{(2,1)}(y)-f_{(2,1)}(y')) \\&= -\epsilon \lambda + \epsilon (2 \lambda - \epsilon ) -2\epsilon (-2\epsilon )= \lambda \epsilon + O(\epsilon ^2). \end{aligned}$$

Hence, F(y) is not one-sided Lipschitz and the uniqueness of the solution of the differential inclusion must be proven in some other manner. One possible approach could be to find a change of variables such that the set-valued drift does become one-sided Lipschitz. Once uniqueness of the solution is established, one can try to use monotonicity arguments to prove global attraction of the fixed point in (3*).

B Calculation of second (raw) moment of the response time

In this appendix, we derive a formula for the second (raw) moment $E[R^2]$ of the response time, which combined with E[R] yields a formula for the variance Var[R]. We have $E[R^2] = R^{*\prime \prime }(0)$. By using (8) together with $G^{*\prime }(0) = -E[G] = -1$, $-{Y^*}'(0)=E[Y]=E[G^2]/2$ and $G^{*\prime \prime }(0) = E[G^2]$, we obtain

$$\begin{aligned} E[R^2]&= {\tilde{\lambda }} Y^{*\prime \prime }(0) \\&\ \ \ +\frac{{\tilde{\lambda }}}{\lambda }\left( \pi ^{\prime \prime }(1) +2\pi ^\prime (1)E[G^2]+(1-\lambda )E[G^2] \right) \\ {}&\ \ \ + \frac{\delta }{\lambda } \left( E[G^2]\frac{m(m+1)}{2} +\frac{(m-1)m(m+1)}{3} \right) . \end{aligned}$$

One readily checks that $Y^{*\prime \prime }(0)= E[G^3]/3$ and we already know that $\pi ^\prime (1) = \lambda E[R]$. From (4), we have

$$\begin{aligned} \pi ^{\prime \prime }(1)&= \frac{\lambda }{{\tilde{\lambda }}}\left( \beta ^{\prime \prime }(1)+2\beta ^\prime (1)\xi ^\prime (1) + \xi ^{\prime \prime }(1) \right) \\ {}&\ \ \ -\frac{\delta (m-2)(m-1)m}{3{\tilde{\lambda }}}. \end{aligned}$$

Making use of (5), one finds

$$\begin{aligned} \beta ^{\prime \prime }(1)= \frac{\delta }{\lambda (1-\tilde{\lambda })}\sum _{i=1}^{m-1} i(i-1) = \frac{\delta (m-2)(m-1)m}{3\lambda (1-{\tilde{\lambda }})}. \end{aligned}$$

We still need to find $\xi ^{\prime \prime }(1)$. Denote respectively by ${\tilde{R}}$ and W the response and waiting time of an ordinary M/G/1 queue with arrival rate ${\tilde{\lambda }}$. By using [13, (5.30)], we get $\xi ^{\prime \prime }(1) = {\tilde{\lambda }}^2 E[{\tilde{R}}^2]$. We have $E[{\tilde{R}}^2]=E[(W+G)^2] = E[W^2]+2E[W]E[G]+ E[G^2]$. As $E[G] = 1$ and $E[W] = \frac{\tilde{\lambda }E[G^2]}{2(1-{\tilde{\lambda }})}$, we obtain $2E[W]E[G]+ E[G^2] = E[G^2]/(1-{\tilde{\lambda }})$. $E[W^2]$ is given by [19, p.256]:

$$\begin{aligned} E[W^2] = \frac{{\tilde{\lambda }}^2E[G^2]^2}{2(1-{\tilde{\lambda }})^2} + \frac{{\tilde{\lambda }} E[G^3]}{3(1-{\tilde{\lambda }})} . \end{aligned}$$

It follows that

$$\begin{aligned} \xi ^{\prime \prime }(1) = \frac{{\tilde{\lambda }}^4 E[G^2]^2}{2(1-{\tilde{\lambda }})^2} + \frac{{\tilde{\lambda }}^3 E[G^3]}{3(1-{\tilde{\lambda }})} + \frac{{\tilde{\lambda }}^2E[G^2]}{1-\tilde{\lambda }}. \end{aligned}$$

Putting everything together, we get

$$\begin{aligned} E[R^2]&= \frac{{\tilde{\lambda }} E[G^3]}{3} + \frac{\delta \tilde{\lambda }(m-2)(m-1)m}{3\lambda (1-{\tilde{\lambda }})} \\&+{\tilde{\lambda }} \frac{\delta m(m-1)}{\lambda (1-{\tilde{\lambda }})}\left( 1 + \frac{{\tilde{\lambda }} E[G^2]}{2(1-{\tilde{\lambda }})}\right) \\&+ \frac{{\tilde{\lambda }}^4E[G^2]^2}{2(1-{\tilde{\lambda }})^2}+ \frac{{\tilde{\lambda }}^3E[G^3]}{3(1-{\tilde{\lambda }})} + \frac{{\tilde{\lambda }}^2E[G^2]}{1-{\tilde{\lambda }}} \\&+ 2{\tilde{\lambda }} E[G^2]\left( 1 + \frac{{\tilde{\lambda }} E[G^2]}{2(1-{\tilde{\lambda }})} + \frac{\delta }{\lambda } \frac{m(m-1)}{2(1-{\tilde{\lambda }})}\right) \\&+\frac{{\tilde{\lambda }}(1-\lambda )}{\lambda }E[G^2] \\&+ \frac{\delta }{\lambda } \left( E[G^2]\frac{m(m+1)}{2} +\frac{(m-1)m(m+1)}{3} \right) . \end{aligned}$$