Parameters optimization and real-time calibration of Measurement-Device-Independent Quantum Key Distribution Network based on Back Propagation Artificial Neural Network

TL;DR

This paper introduces a BPNN-based method for rapid parameter optimization and real-time calibration in large-scale MDI-QKD networks, reducing resource consumption compared to traditional optimization algorithms.

Contribution

It presents a novel application of Back Propagation Artificial Neural Networks for efficient parameter prediction and system calibration in MDI-QKD networks, enhancing scalability and resource efficiency.

Findings

BPNN predicts optimal parameters faster than traditional algorithms.

The method enables real-time system calibration using discarded data.

Resource savings in large-scale quantum key distribution networks.

Abstract

The parameters choosing (such as probabilities of choosing X-basis or Z-basis, intensity of signal state and decoy state, etc.) and system calibrating will be more challenging when the number of users of a measurement-device-independent quantum key distribution(MDI-QKD) network becomes larger. At present, people usually use optimization algorithms to search the best parameters. This method can find the optimized parameters accurately but may cost lots of time and hardware resources. It's a big problem in large scale MDI-QKD network. Here, we present a new method, using Back Propagation Artificial Neural Network(BPNN) to predict, rather than searching the optimized parameters. Compared with optimization algorithms, our BPNN is faster and more lightweight, it can save system resources. Another big problem brought by large scale MDI-QKD network is system recalibration. BPNN can support this work in real time, and it only needs to use some discarded data generated from communication process, rather than require us to add additional devices or scan the system

Tables

Table 1. Table 1: system conditions for numerical simulation Where d ​ c 𝑑 𝑐 dc denotes dark count, η d subscript 𝜂 𝑑 \eta_{d} denotes detector efficiency, e d subscript 𝑒 𝑑 e_{d} denotes misalignment, f 𝑓 f denotes error correcting efficiency and N 𝑁 N is data size

$dc$	${\eta }_d$	$e_d$	$f$	$N$
$6.02\ast {10}^{-8}$	$70\%$	$1.5\%$	$1.16$	${10}^{12}$

Table 2. Table 2: Charlie’s conditions

$dc$	${\eta }_d$	$f$	$N$
$6.02\times {10}^{-8}$	$70\%$	$1.16$	${10}^{12}$

Table 3. Table 3: system conditions

$L_a$	$L_b$	$dc$	${\eta }_d$	$f$	$N$
$10km$	$20km$	$6.02\times {10}^{-8}$	$70\%$	$1.16$	${10}^{12}$

Equations

$$\vec{v}=[{\mu }_{Za},{\nu }_{Za},{\mu }_{Xa},{\nu }_{Xa},P_{Za}^{\mu },P_{Za}^{\nu },P_{Xa}^{\mu },P_{Xa}^{\nu },$$

$${\mu }_{Zb},{\nu }_{Zb},{\mu }_{Xb},{\nu }_{Xb},P_{Zb}^{\mu },P_{Zb}^{\nu },P_{Xb}^{\mu },P_{Xb}^{\nu }]$$

$${\vec{v}}_{opt}=argma{x}_{\vec{v}\in V}[R(\vec{v})]$$

$${\vec{v}}_{opt}=f(\vec{s});{\vec{v}}_{opt}\in V,\vec{s}\in S$$

$${ϵ}_i=\sum _{m=1}^M{w}_{mi}\times x_m$$

$${\zeta }_i=f({ϵ}_i-{\theta }_i)$$

$${ϵ}_p=f(\sum _{j=1}^P{w}_{jp}\times {\zeta }_j-{\theta }_p)$$

$${\vec{v}}_{\pm }=[{\mu }_{Za}\pm \Delta s,{\nu }_{Za},{\mu }_{Xa},{\nu }_{Xa},P_{Za}^{\mu },P_{Za}^{\nu },P_{Xa}^{\mu },P_{Xa}^{\nu },$$

$${\mu }_{Zb},{\nu }_{Zb},{\mu }_{Xb},{\nu }_{Xb},P_{Zb}^{\mu },P_{Zb}^{\nu },P_{Xb}^{\mu },P_{Xb}^{\nu }]$$

$$r_{\pm }=S({\vec{v}}_{\pm })$$

Full text

††preprint: APS/123-QED

I Introduction

Quantum key distribution(QKD) allows two remote users, Alice and Bob, share random keys with information-theoretic security Shor and Preskill (2000); Lo and Chau (1999). However, there are still some loopholes since imperfect of devicesLütkenhaus (2000); Brassard et al. (2000); Fung et al. (2007); Xu et al. (2010); Zhao et al. (2008); Gerhardt et al. (2011); Makarov (2009). In present, measurement-device-independent quantum key distribution(MDI-QKD)Liu et al. (2013); Lo et al. (2012); Ma and Razavi (2012) may be the most promise and attractive QKD protocol, because all potential side channel attacks on measurement device are removedLo et al. (2012).

In recent years, people did many efforts to further improve the practicality of MDI-QKD. They present decoy state methodWang (2013); Yu et al. (2015); Zhou et al. (2016) and finite-key analysisCurty et al. (2014); Ma et al. (2012) to improve the practical security, Reference-Frame-Independent(RFI)Wang et al. (2015); Zhang et al. (2017) scheme and polarization randomizing schemeWang et al. (2017) to make the MDI-QKD system environmental robustness, etc.

It is quite natural to extend MDI-QKD structure to a ”star” networkTang et al. (2016); Fröhlich et al. (2013), in which all users connect to the untrusted relay Charlie. The problem is, when calculating MDI-QKD key rate and doing experiments, we always suppose two users have the same distance to Charlie (we could call it ’symmetric situation’, and call Charlie’s place ’symmetric site’. Otherwise, the ’asymmetric situation’ and ’asymmetric site’), which is difficult to reach in the real world, especially in making MDI-QKD network since there are too many users and is nearly impossible to place Charlie to symmetric site.

The secure key rate of MDI-QKD is heavily depend on symmetry level of two users’ channel losses since the protocol need Hong-Ou-Mandel(HOM) interference. When losses are asymmetric,the key rate decay fast. Recently, Wenyuan Wang, et al. present a ’7-intensity protocol’Wang et al. (2018) to deal with the problem. And we will introduce another asymmetric protocol in this paper.

The optimizations of parameters are crucial in QKD. It influences secure key rate significantlyXu et al. (2014); Wang et al. (2018). In MDI-QKD, each pair of users should keep a set of optimized parameters. Thus in large scale MDI-QKD network, parameter distribution could be a big challenge. For example, there are 10000 users in network, $10000\times 9999/2$, i.e about $5\times 10^{7}$ sets of parameters are needed. It may take unbelievably long time to generate them. If we need 1 second(actually may be far more than one second) to optimize a pair of parameters, $5\times 10^{7}$ seconds means more than one year! When a new user join the network. We should optimize 10000 pairs of new parameters for him. What’s worse, the condition on users’ side may change and we should adjust parameters to adapt this changing. This may cost plenty of resources too.

Previously, optimized parameters were generated by ’simulation and iterationXu et al. (2014); Wang et al. (2018). This method is mature and accurate —- if you give enough time and iteration times. However, simulation programs usually cost plenty of time, especially, when running optimization program, you have to iterate simulation program many times. Consider building a large scale MDI-QKD network, when there are a new user link to the network or some users need to re-optimize their parameters for some reasons, we have to optimize too many parameters since each pair of users should have their own optimized parameters.

Another problem brought by large scale MDI-QKD network is reference phase difference($\Delta{\phi}$) and misalignment($e_{d}$)Tang et al. (2014); Wang et al. (2015) on users’ side may change with time. Thus we must recalibrate users’ system in time. It’s really a big projection with traditional methods.

Previously, system parameters were re-calibrated by ’stop and scan’Wang et al. (2015) or ’add phase stabilization’Tang et al. (2014). The former needs too much time and the latter needs too many additional devices. In this paper, we present a scheme based on Back Propagation Artificial Neural Network(BPNN)Kohonen (1988); Rumelhart et al. (1986) to solve above problems. BPNN could not only calculate the optimized parameters directly, but also estimate $\Delta{\phi}$ and $e_{d}$. We could recalibrate the system according to the result given by BPNN rather than stopping and scanning the system or adding phase stabilizations. We validate this idea by simulation and find it is indeed useful.

The rest of this paper was organized as follows. In Sec.2, we briefly introduce an asymmetric protocol base on 3-intensity protocolWang (2013); Yu et al. (2015), our BPNN training data was generate by this protocol. In Sec.3, we introduce how to make a BPNN to calculate optimized parameters and compare it with traditional local search algorithm (LSA)Xu et al. (2014), and we will discuss the potential of BPNN. In Sec.4, we introduce methods to make and use BPNN to estimate the $\Delta{\phi}$ and $e_{d}$, we give simulation data to prove the validity of this way. Finally in the appendix A and B, we introduce the working principle of BPNN and process of LSA.

II Asymmetric Protocol Based on 3-intensity protocol

In asymmetric situations, asymmetric protocols perform much better than symmetric versionWang et al. (2018); Liu et al. (2018). In this section, we firstly introduce an asymmetric protocol derived from 3-intensity methodWang (2013); Yu et al. (2015). We will generate our BPNN dataset by this protocol. Original 3-intensity method considers signal and decoy state intensity in Z-basis($\mu_{Z}$ and $\nu_{Z}$) and X-basis($\mu_{X}$ and $\nu_{X}$), probability to choose signal and decoy state in Z-basis($P_{Z}^{\mu}$, $P_{Z}^{\nu}$) and in X-basis ($P_{X}^{\mu}$, $P_{X}^{\nu}$), i.e. 8 parameters in total. We don’t consider probability of vacuum state($P^{o}$) because $P^{o}=1-P_{X}^{\mu}-P_{X}^{\nu}-P_{Z}^{\mu}-P_{Z}^{\nu}$, and don’t distinguish vacuum in Z and X-basis since vacuum state don’t have any photons.

In our asymmetric protocol, Alice and Bob could choose different parameters. Thus we should optimize 16 parameters totally. We denote these parameters as vector $\vec{v}$ showed in formula(1). Subscripts a(b) denotes Alice(Bob).

[TABLE]

When all system conditions(including channel loss, misalignment, detector efficiency, data size and dark count rate)are fixed, secure key rate($r$) is function of $\vec{v}$, i.e. $r=R\left(\vec{v}\right)$. Parameters optimization could be regarded as searching for $\vec{v}_{opt}$ to maximize $r$ as showed in Eq.(2). Where $V$ denotes searching space of $\vec{v}$.

[TABLE]

Here we fix system conditions as showed in Tab.1 and compare key rate($r$) vs channel length($L_{a}$ , $L_{b}$) in different protocol. The result is generated from numerical simulationWang and Wang (2014) and showed in Fig.1, Fig.2 and Fig.3. Fig.1 is the result of our asymmetric protocol, we could find in asymmetric situations, it has higher key rate and reaches a longer distance than 3-intensity method showed in Fig.2. Then we limit $L_{a}+L_{b}=70km$ and compare two protocols. The result is showed in Fig.3 We could find that our protocol performs better than the 3-intensity protocol in strong asymmetric situations.

III predict optimized parameters by BPNN

We use ’predict’ in the title to emphasize our method is different from traditional methods. We calculate best parameters directly with the help of machine learning(ML) rather than use simulation and iteration to optimize them. ML is an algorithm that uses statistical techniques to give computer systems the ability to ’learn’. It is widely used in regression, classification, prediction, and clustering, etc.

Eq.(2), on the one hand, told us given system conditions, we could find a set of parameters $\vec{v}_{opt}$ to maximize $r$. It enlighten us there should be function relationship between system conditions($\vec{s}$) and $\vec{v}_{opt}$. We express it as Eq.(3). Where $V$ denotes searching space of $\vec{v}$, $S$ denotes space of $\vec{s}$.

[TABLE]

Function $f$ in Eq.(3) is quite difficult to find analytically. But we could approximate its performance by ML according to Universal Approximation TheoremCsáji (2001); Cybenko (1989), which says simple Artificial Neural Networks(ANN) can represent a wide variety of interesting functions when given sufficient data and appropriate hyperparameters.

We choose BPNN, a widely used ANN, to fitting the $f$ since its good ability of nonlinear mapping, self-adapting and generalization. As showed in Fig.4. BPNN include an input layer, several hidden layers and an output layer. Each layer consists of neurons and connects to adjacent layers by synapses.The data recieved by input layer, forward propagate and operated by synapses and nuerals, finally outputted by the output layer.

System conditions $dc$,$\eta_{d}$,$f$,$N$ are depended on operator’s devices. We call them Charlie’s conditions and regard them as fixed value. $L_{a}$,$L_{b}$ are depended on distance between users and Charlie. $e_{d}$ is mainly depended on users’ misalignment of quantum state preparation. We call them users’ conditions. Our purpose is to predict the best parameters according to users’ conditions. To improve the precision, each BPNN only predicts one element in $\vec{v}$. Thus we need totally 8 BPNNs, and we will introduce the BPNN which is used to predict $\mu_{Za}$ and $\mu_{Zb}$ as representative in the following.

Firstly we create an input layer with 3 neurons. They receive $L_{a}$,$L_{b}$ and $e_{d}$ as input respectively. Then we add a hidden layer consists of $20\times 20$ neurons to link with input layer. Finally we place an output layer with only one neurons, which will output the target parameter.

The next step is preparing training data for BPNN. This step uses traditional simulationWang and Wang (2014) , LSA and Particle Swarm Optimization Algorithm (PSO). The PSO can optimize the non-smooth function and non-convex function, but its speed is very slow. that is to say, PSO can find more accurate optimized parameters. The steps of dataset preparation is introduced as following

Step.1 Fixing Charlie’s conditions as Tab.2. Then choose users’ conditions as $\vec{u}^{0}=(L_{a}^{0},L_{b}^{0},e_{d}^{0})$. Using PSO to find the optimzied point $\vec{v}_{opt}^{0}$.

Step.2 Changing $\vec{u}^{0}$ to $\vec{u}^{1}$, but the difference should be very small. Choose $\vec{v}_{opt}^{0}$ as initial point and using LSA to optimize the parameters. The result is denoted as $\vec{v}_{opt}^{1}$. We trust the difference between $\vec{v}_{opt}^{0}$ and $\vec{v}_{opt}^{1}$ should be very small since the difference between $\vec{u}^{0}$ and $\vec{u}^{1}$ is small. LSA can be faster and more accurate when it has better initial point.

Step.3 repeat Step.2 to generate more pairs of $\vec{u}^{i}$ and $\vec{v}_{opt}^{i}$. When $\vec{u}^{i}$ is special point(such as $L_{a}=L_{a}^{0}$), using PSO algorithm rather than LSA.

Step.4 Record all $\vec{u}^{i}$ and $\vec{v}_{opt}^{i}$ as our dataset.

We only use PSO on some special user’s conditions. It’s a tradeoff between time and precision.

The final step is training. We introduce the BPNN which is used to predict $\mu_{Za}$ and $\mu_{Zb}$ as representative. We extract first element ${\mu_{Za}^{i}}$ from each $\vec{v}_{opt}^{i}$ and make a label dataset $V_{1}=\{\mu_{Za}^{i}|1\leq{i}\leq{M}\}$. Then put dataset $U$ into BPNN, and adjust weight of synapses with the back propagation algorithm(BPA) according to the difference between output of BPNN and label dataset $V_{1}$. We will introduce the working principle of BPNN in detail in Appendix. It is worth noting that, when predicting $\mu_{Zb}$, we should switch $L_{a}$ and $L_{b}$ first.

To verify the validity of our BPNN, we randomly generate 1500 sets of users’ conditions, then calculate optimized parameters by LSA and BPNN respectively in my PC(CPU:Intel Core i5 7500@ 3.40GHz; RAM: DDR4 8GBytes). We compare the running time and key rate.The result is showed in Fig.5. It indicates that, comparing with troditional LSA, the speed of our method is more than 100 times faster and key rate loss is less than 20% with the probability of 80%. The biggest superiority of BPNN is the average time running is only 0.13 second. It’s difficult for traditional methods to reach this speed in same hardware condition since the limitation of simulation program.

BPNN has another particular potential that their essence a is simple formula. It means that we could build BPNN by hardware(such as embedded systems) to further improve the speed(specially designed system usually runs faster). But it’s nearly impossible to design a hardware system for simulation program since it’s too complex.

IV Real-Time Calibration by BPNN

In MDI-QKD, reference phase difference between two users $\Delta{\phi}$ lead to errors in X-basis. Quantum states misalignment $e_{d}$ lead to errors in both Z and X-basis. To make the problem worse, these misalignments may change with time. Thus, MDI-QKD requires us to calibrate users’ system in time. As discussed in the introduction, it’s quite difficult to accomplish the mission in a large scale network. Consider there are $n$ users in the network. If we try to add phase stabilization, about $n\times(n-1)/2$ phase stabilizations are needed. If stopping and scanning the system, more than $n$ times scan are needed. Obviously, big $n$ may lead to system unable to work.

Consider a situation as follows. Two users in the network have fixed distance $L_{a}$,$L_{b}$. They optimized their parameters as $\vec{v}_{opt}$ under the conditions of $\Delta{\phi}=0$ and $e_{d}=0.002$. But $\Delta{\phi}$ and $e_{d}$ increase with time.

With the help of BPNN, we could estimate the $\Delta{\phi}$ and $e_{d}$ in real time, only need to statistic data of single photon error rate in X-basis ($e_{11}^{X}$), qubit error rate of signal state in Z-basis ($E_{\mu}^{Z}$) and secure key rate ($r$), instead of adding additional devices or stopping the system. It’s worth noting that these data are calculated for error correction and privacy amplification originally, i.e. they are recycled, not generate for BPNN specially.

An other superiority of our method is speed. BPNN could judge how much $\Delta{\phi}$ and $e_{d}$ are there in users’ system accurately and quickly. If we have known the value, we could recalibrate users’ system immediately.

Apart form these, there are two main superiorities compared with using traditional curve fitting to estimate $\Delta{\phi}$ and $e_{d}$. The first is the ability ’Non-linear mapping’. there is a theory called Universal approximation theorem Csáji (2001); Cybenko (1989). It states that a feed-forward network with a single hidden layer containing a finite number of neurons can approximate continuous functions on compact subsets, under mild assumptions on the activation function. The theorem thus states that simple neural networks can represent a wide variety of interesting functions when given appropriate parameters. When using traditional curve fitting, we should have some priori knowledge about the data we want to fitting and carefully choose a mode such as linear, logarithmic or polynomial etc. But, when using BPNN, we don’t care about the mode, we just prepare our dataset and push them to BPNN to do the training.

Another one is the ability to process high dimensional data. Artificial neural networks can easily process high dimensional data such as a photo with millions dimension. But curve fitting can’t be up to this job. Although in the problem of estimating $\phi$ and $e_{d}$, the dimension is not too high (input is three dimension and output is two dimension). But we want to give a demonstration and inspiration to others. Maybe in the future, we can use more input data to correct the polarization misalignment, time misalignment etc. And people may find more regular of QKD system by AI algorithm.

The BPNN for calibration need $e_{11}^{X}$, $E_{\mu}^{Z}$, and $r$ as its input. $e_{11}^{X}$ mainly implies the $\Delta{\phi}$, $E_{\mu}^{Z}$ implies the $e_{d}$, and $r$ is sensitive to all kinds of errors. Thus we choose them to be our input and use BPNN to find their relations.

Firstly we prepare our training dataset. The process is showed as following:

Step.1 We fix Charlie’s conditions and $L_{a}$ $L_{b}$ as showed in Tab.3. Choosing phase misalignment and misalignment error as $\Delta{\phi}^{0}$ and $e_{d}^{0}$. Then we simulate the MDI protocol, get the $e_{11}^{X0},E_{\mu}^{Z0},r^{0}$

Step.2 Changing phase misalignment and misalignment error to $\Delta{\phi}^{1}$ and $e_{d}^{1}$. Then we simulate the MDI protocol, get the $e_{11}^{X1},E_{\mu}^{Z1},r^{1}$

Step.3 Repeating the step.2 and record all $\Delta{\phi}^{i}$, $e_{d}^{i}$, $e_{11}^{Xi},E_{\mu}^{Zi},r^{i}$ as our training dataset. The $\Delta{\phi}^{i}$, $e_{d}^{i}$ are our output dataset and $e_{11}^{Xi},E_{\mu}^{Zi},r^{i}$ are our input dataset.

To verify the validity, we randomly generate 1000 pairs of $\Delta{\phi}^{i}$ and $e_{d}^{i}$, where $\Delta{\phi}$ is uniformly distribute in 0 to 0.5 and $e_{d}$ is uniformly distribute in 0.002 to 0.02. Then we run the simulation program to obtain related $e_{11}^{Xi}$, $E_{\mu}^{Zi}$ and $r^{i}$. We put $e_{11}^{X}$, $E_{\mu}^{Z}$ and $r$ into BPNN, then BPNN estimate these misaligments, denoted by $\Delta{\phi}_{n}et$ and $e_{d\_net}$. We compare the difference between true misaligments and estimated misaligments. The result is showed in Fig.6. It indicate that difference between estimated value and the true value is about 2 order of magnitude smaller than the true value, i.e. the estimation is quite accurate.

V Discussion

In summary, We firstly introduce a method to develop symmetric 3-intensity protocol to adapt the asymmetric situations and prove it indeed has higher key rate and reaches a longer distance by simulation. Our BPNN training dataset is generated by this method.

Then, We presented a totally new idea to get optimized parameters. With the help of BPNN, we don’t need to use optimization algorithm any more. We verified the performance of BPNN and found it could run far faster than traditional methods and analyzed the reason for its faster speed in principle. We finally discussed the potential of BPNN which optimization algorithm could never reach.

Also, we found BPNN could support to calibrate the system in time, and verified it by simulation. In simulation results, the estimations of BPNN are very accuracte. We trust this is a valuable idea since at present, calibration methods cost too much additional time or devices. Our method gives a possibility to improve the system performance only with discarded data.

VI ACKNOWLEDGMENTS

This work has been supported by the National Key Research and Development Program of China (2016YFA0302600), the National Natural Science Foundation of China (NSFC) (61775207, 61575183, 61771439, 61702469); Strategic Priority Research Program (B) of the Chinese Academy of Sciences (CAS) (XDB01030100, XDB01030300), and the Foundation of Science and Technology on Communication Security Laboratory (6142103040105).

VII Appendix A : Working principle of Back Propagation Artificial Neural Network

Artificial neural networks (ANN) consist of simple units with adaptability connect with each other. It is inspired by biological neural networks. The simple units come from neurons and connections come from synapses. A neuron could receive a signal, process it and transmit it to additional neurons connected to it. Error Backpropagation AlgorithmRumelhart et al. (1986) is used to calculate a gradient that is needed when training weights of synapses. The combination of them are called Back Propagation Artificial Neural Network(BPNN).

We will introduce BPNN by Fig.4. The red circles are input neurons, they receive $x_{m}$ as input signals, where $m=1,2,....,M$. Then the signals go through the synapses with certain weights and summed by hidden neurons’ input side as showed in Eq.(4).

[TABLE]

Then, the $i$’th hidden neurons, denoted by blue circles, receive $epsilon_{i}$ as its input signals and process them by active function $f$ as showed in Eq.(5). Where $\theta$ is threshold.

[TABLE]

The following steps are the same as above, $zeta_{i}$ go through the synapses with certain weights and summed as inputs of the next layer, processed buy active functions and so on, until the output layer gives results $y$ as Eq.(6).

[TABLE]

VIII Appendix B : Using Local Search algorithm to optimize parameters

We use LSA to generate optimized parameters as training data of BPNN. The process is described as follows.

Step.1 Setting up precision value $\tau$ to shut down the algorithm and step length $s$ to search better value.

Step.2 We choose a good initial point $\vec{v}$ as showed Eq.(1) and calculate the key rate by simulation program, i.e. $r=S(\vec{v})$ , where $S$ is simulation program which receive parameters $\vec{v}$ as input and output key rate $r$.

Step.3 The first element in $\vec{v}$ is $\mu_{Za}$, we move it to $\mu_{Za}+\Delta{s}$ and $\mu_{Za}-\Delta{s}$. Then we get two new points $\vec{v}_{+}$ and $\vec{v}_{-}$ as Eq.(7) and calculate key rate at these new points as Eq.(8).

[TABLE]

[TABLE]

Step.4 Compare $r$, $r_{+}$ and $r_{-}$. if $r_{+}$($r_{-}$) is the maximum one, move our searching point to $\vec{v}_{+}$($\vec{v}_{-}$), i.e. $\vec{v}=\vec{v}_{+}(\vec{v}_{-})$. Else if $r$ is the maximum, keep $\vec{v}$ unchanged. And record the maximum key rate as $r_{max}$.

Step.5 We compare $r$ and $r_{max}$, if the difference of them is less than $\tau$, shut down the algorithm and output $\vec{v}$ as optimized parameters. Else, let $r=r_{max}$ and repeat step.3 to step.5, but move the next element in $\vec{v}$ this time. If it has been the last element in this loop, return to the first element to move.

It is worth noting that, we should choose $s$ and $\tau$ carefully to balance the precision and speed.