Face Recognition with Machine Learning in OpenCV_ Fusion of the results with the Localization Data of an Acoustic Camera for Speaker Identification

This contribution gives an overview of face recogni-tion algorithms, their implementation and practical uses. First, a training set of different persons' faces has to be collected and used to train a face recognizer. The resulting face model can be utilized to classify people in specific individuals or unknowns. After tracking the recognized face and estimating the acoustic sound source's position, both can be combined to give detailed information about possible speakers and if they are talking or not. This leads to a precise real-time description of the situation, which can be used for further applications, e.g. for multi-channel speech enhancement by adaptive beamformers.

I. INTRODUCTION
In recent times, the interest in and research on acoustic source localization and enhancement of certain sound sources has increased dramatically due to the growing desire for handsfree interaction with various devices [18]. Combining the ability to locate sound sources and to recognize possible speakers with a camera potentially enables machines to identify speakers. This makes a human-machine-interaction a lot easier, more adaptive and reliable. Comparable systems can be used for teleconferencing, smart rooms or ambient assisted living [5,8].
Combining a microphone array's ability to locate sound sources and the intuitive way of extracting information of a webcam's image, acoustic cameras have become quite popular for many industrial segments. They compute a color-coded sound map and thus, they visualize the sound pressure levels of a user-defined field of view. This way, acoustic cameras can locate sound sources quite accurately, which is why they are often used to identify unwanted noise sources [12,31].
Face recognition is a machine learning technique, which ideally allows detecting and identifying all faces seen in a picture or a video frame. It can be used for criminal detection, image processing, human computer interaction, etc. [33,34]. In the early development of face recognition systems, geometric facial features, e.g. eyes, nose and mouth, were explicitly used. Properties of these features and relations (e.g. positions, distances, angles) between them were used as descriptors for face recognition [15]. Today, holistic techniques, e.g. principal component analysis (see Eigenfaces) or linear discriminant analysis (see Fisherfaces), are used to identify individuals [3,13].

II. BASICS OF AN ACOUSTIC CAMERA
The implementation of an acoustic camera requires a suitable microphone array as well as beamforming algorithms to locate sound sources precisely. Given both, a color coded sound map of the measured sound pressure level can be computed and displayed as seen in Figure 3 [25,26].

A. A suitable microphone array
In [26] a appropriate microphone array (see Figure 1) for speaker and sound source localization has been developed and verified. It can be shown that double ring arrays with an odd number of microphones on each ring are desirable for locating speech sources [20]. In this project, the inner ring has a diameter of 0.2 m, while the outer ring is twice as large. An important part of acoustic cameras is a sound analysis and visualization software [31]. This software can for example be installed on a personal computer. The connection of microphones with any computer is achieved by using a microphone amplifier and a multi-channel sound card. In particular, two RME OCTAMIC XTC amplifiers and a RME MADIface USB are used. Both amplifiers digitize analog signals of up to eight channels completely synchronously. By interconnecting two RME OCTAMICs in series, up to 16 analog signals can be converted to digital values. Thus, in order to read out all signals synchronously, it is necessary to activate a delay compensation in each amplifier [27].
Similar to temporal undersampling, which causes temporal aliasing effects, spatial undersampling can lead to spatial aliasing. This effect can be observed in acoustic cameras' color maps as incorrectly detected sound sources [30]. In order to minimize spatial aliasing, there are multiple approaches possible. In [20] it is shown that an odd numbers of microphones on each ring of the array can reduce redundancy, which results in more robustness against spatial aliasing. Ring arrays in general decrease the redundancy of microphone arrays, because at a certain frequency only a few microphone pairs are affected by spatial aliasing, while others are not yet [9]. To build a sensor array, utilizing omnidirectional microphones, such as the selected condenser microphones AKG CK-92, has been shown advantageous [8,25,26]. Using more microphones and distribute them randomly on a plane are two possible improvements of an acoustic camera's microphone array.

B. Steered Response Power Beamforming Algorithm
Beamforming algorithms process signals in a way, so that desired directions are enhanced, while signals from all other directions are attenuated. This chosen direction can be called steering direction, with which a defined plane can be spatially sampled. The beamformer's output, when used in this way, is known as the steered response [12]. As seen in Figure 2, the Steered Response Power Algorithm with Phase Transform weighting (SRP-PHAT) calculates a color map by summing certain values of signal pairs' weighted cross correlation (GCC). It is broadly known that a correlation results in a signal's power and thus, the steered response outputs a power. This is why the described method is known as a steered response power beamformer [11,12].  [17] The mentioned GCC is similar to a regular cross correlation with the only difference, that weighted input signals are used. In order to get a PHAT weighting the Fourier transformed signals X1(ω) and the complex conjugate of X2(ω) are used as seen in (1). TDOA stands for time difference of arrival and it estimates the time difference of two sensors' signals. The TDOAs for a single microphone pair differ with the steering direction [11,25]. (1) In [25] it has been shown, that best sound source localization results can be achieved by combining the SRP-PHAT algorithm with a constantly weighted SRP beamformer. This is because the SRP-PHAT is unable to process narrowband signals, while being very robust against reverberations and sensor self noise. The SRP method is not as robust as the SRP-PHAT algorithm, but in contrast to that, it is able to locate narrowband signals such as sine waves or spoken vowels. A combination of both is implemented by utilizing a threshold for the signal's bandwidth. In this application, the bandwidth's threshold is set to 4 kHz, which is approximately an eighth of the chosen sampling rate. A typical output image of an acoustic camera can be seen in Figure 3. An ideal face detection system should be able to detect all faces shown in a picture or a video frame. For this task, it should neither matter in which position or orientation the faces are, nor which age, sex or ethnical origin the people to be classified belong to. Furthermore, an ideal face detection system should be insensitive to lighting changes or other external influences [16].
In OpenCV, a face library is implemented, which provides pre-trained face detectors as well as the possibility to train own classifiers [2,22]. Pre-trained classifiers for Haar-like and local binary pattern features support frontal face, facial landmark and whole person detection. If a self-trained classifier is used, several thousand pictures of non-faces and faces should be collected. A good training set considers faces with differences in age, sex, ethnical origin, facial hair, lighting and hairstyle [6,13,21]. Because of the complexity of the training process, only pre-trained classifiers for frontal faces are utilized in this contribution. While using the later described face detection algorithms utilizing Haar-like or local binary pattern features, often multiple faces result from one facial image. If these detected faces are located in a specific area close to each other, they are averaged in size and position to merge them into one detection result. This avoids multiple detections for a single person and reduces false positive rates [2,28].
Both, the Haar cascade and the local binary pattern classifier are implemented as cascaded classifiers to quickly reject non-faces but still keep a high accuracy for positive results (see Figure 4).

A. Haar Cascade Classifier
The Haar cascade classifier is a quite easy face detection method, and is therefore a very good basis for more complex algorithms. With a huge dataset, many different objects can be trained, e.g. faces, cars or whole persons. In order to classify images, Haar-like features (see Figure 5) are used and calculated extremely efficiently with integral images. Thus, regional knowledge can be considered. As very common, the face detector introduced in [35] can only handle grayscale images [6,15,35].  [35] In pictures with a resolution of 24x24 Pixel, more than 180,000 different Haar-like features can be found. Using a machine learning (ML) algorithm, the 6,061 most important features can be chosen and organized in a cascade structure. Training the chosen features fj results in a threshold θj and a parity pj. With features being all pixel values added in black blocks and subtracted from the sum of pixel values in the white blocks, the weak classifier can be described as [35] (2) In the first five stages of the cascade, 1, 10, 25, 25 and 50 of these weak classifiers are utilized to differ between face and non-face images. The number of weak classifiers per stage is determined by a defined false positive rate, which has to be achieved in each stage. Thus, it can be imagined that in the first few stages, only few features are necessary to get to this rate, but at the very last stage, very many are needed. Stages are added as long as a total false positive rate is met [35].
Using local binary patterns, more sophisticated Haar-like features or additional non-frontal face detectors, higher detection accuracies can be achieved [13,28,35].

B. Local Binary Pattern Classifier
Local binary patterns (LBPs) describe local relationships between neighboring pixels in a 3x3 environment. Starting in the top left corner and proceeding clockwise, the pixels' grayscale values are compared to the center pixels'. If the value of the center pixel at (n,m) is bigger than the neighbor's value, a 0 results, 1 else. These binary values can be put together and converted into a grayscale value of 0…255. Formally, this can be written as [16,28] (3) With s(x) being 0, if x<0, 1 else. It is quite clear, that the LBP is invariant to monotone grayscale transformations, e.g. changes in brightness or contrast. LBPs contain most information for maximum two changes between 0 and 1. Examples therefore can be seen in Figure 6 [28].  [23] Similar to the stage structure of a Haar cascade classifier, a weak classifier is formed by using a gray value histogram. With Hn(X) being the classifier for the n-th stage and h(n,m)(x) being the histogram value for listing x, the weak classifier can be described as in (4). The LBP classifier is not only faster than the Haar cascade classifier, but theoreticaly even more precise [3,28].

IV. FACE RECOGNITION
Even though face recognition is a much more challenging task than face detection is, today's face recognition systems are, at least under optimal conditions, very reliable [2,14]. Thus, many of today's applications use these methods to identify people in images, e.g. Facebook's Gallery or Apple iPhoto [32]. When classifying, problems mostly occur due to variations in light, perspective or facial expression [14,32]. Fur-thermore, similar looking individuals, e.g. father and son or twins, can cause uncertainties when differing them [15].
In order to get a higher face recognition accuracy, different approaches can be imagined. In general, it is recommended to use large datasets with many variations in pose, age and lighting conditions for training the model [3,13,32]. Another possibility to improve the recognition performance is to use infrared lighting to avoid shades or other disruptions. Furthermore, additional features, which only occur under invisible light, e.g. freckles and pigmentation, can be used to recognize faces [15].
In OpenCV, face recognizers using principal component analysis (Eigenfaces), linear discriminant analysis (Fisherfaces) and local binary pattern histograms are implemented [23]. Utilizing any of these methods, the face recognizer has to be trained with own face-images to differentiate between individuals. The classification is done by comparing the images' features in a high-dimensional feature space with a K-nearest neighbor algorithm [6].

A. Eigenface Classifier
The Eigenface classifier uses a principal component analysis (PCA) to reduce the dimensionality of the images. Utilizing a PCA, E eigenvectors with the highest eigenvalues can be selected to describe the given dataset. These eigenvectors span a quite low-dimensional face space, in which every image can be projected. Because the PCA's eigenvectors, after reshaping them into an image format, look very much like faces (see Figure 7), they are called Eigenfaces [4,34].  [23] Using the extracted Eigenfaces, unknown faces can be reconstructed (see Figure 8). This gives a good feeling about how many principal components are necessary to distinguish between individuals. Usually a number of 40 to 80 should be sufficient, but, depending on the dataset, sometimes up to 300 Eigenfaces should be used [13,23,36]. Figure 8 shows that the original face can be recognized starting at 20 Eigenfaces.
In order to calculate a dataset's (Γ1, Γ2, Γ3, …, ΓR) eigenvectors efficiently, its vectorized mean image Ψ and the differences Φi = Γi -Ψ are needed. Given a face matrix A = [Φ1 Φ2 … ΦR], the covariance matrix of all face-images can be calculated as [34] (5) This Matrix C inherits a dimension of N²xN² (for face images with a resolution of NxN), which means that N² eigenvektors uk have to be determined. This requires high computational resources and thus, it is unsuitible for real-time applications. For the common case that R<N², using a workaround by building a RxR dimensional Matrix L = A T A, only R eigenvectors υl are to be calculated. Utilizing them, the originally disired eigenvektors ul can be determined as [19,34] (6) Figure 8: Reconstruction of a face using Eigenfaces, according to [23].  (7). The resulting weights can be used to descibe a weight-vector Ω = (ω1, ω2, …, ωE) T [34].

B. Fisherface Classifer
Dimensionality reduction by linear discriminant analysis (LDA) can counter the PCA's disadvantage, not considering any class dependencies while projecting images in a feature space. Using Fisher's linear discriminant analysis (FLD), the classes stay linearly separable, which makes classification easier and more reliable, especially for changes in lighting condi-tions. Thus, E orthonormal vectors describe a matrix W, so that it maximizes between-class scatter (see (8)) but minimizes within-class scatter (see (9)). For both, matrices SB, SW have to be defined (8) (9) where C is the number of different classes, Ni the number of test images of class Xi and Ψi is its mean image. For face recogniton tasks, the LDA projection Wopt' can be written as in (10). Using for example a PCA, the dimension reduces to N-C, while the FLD reduces it further to C-1 [4].
The Fisherface methods provides better handling for background and lighting, than Eigenfaces do [13]. Furthermore, Fisherfaces are much more reliable when using a small training set or faces differing heavily from the training data, e.g. wearing glasses or facial expressions [4,23].

C. Local Binary Pattern Histogram Classifier
The classification using local binary pattern histograms (LBPH) is quite similar to the face detection with LBP. Differences are that in order to identify individuals, a K-nearest neighbor method is utilized and the LBP operator can be extended to get results that are more reliable. For this, multiple approaches are possible. Instead of considering eight direct neighbors, P neighboring pixels on a radius R can be used for the generalized LBP P,R operator [1,15,23]. Another option is to use a multi-block LBP to compare the avarage gray scale values of neighboring pixel blocks with the avarage of a centered region [15].
The LBPH classifier's main disadvantage is that it is quite slow and therefore unsuitable for fluent video playback in realtime situations [1]. Thus, as described in the following chapter, it cannot be used in the implemented application.

V. APPLICATION OF A FACE RECOGNITION SYSTEM
In order to compensate changes in lighting, face rotation, background and hairstyle, some preprocessing steps are taken before recognizing faces (see Figure 9). It can be assumed that this allows to apply the classifiers not only to constrained environments, but also to any [34].
As a suitable face detection algorithm, the Haar cascade classifier is chosen. Even though, OpenCV's LBP classifier shows in practice an approximately 61 % faster processing time, it is less accurate and is unable to detect faces reliably in an artificially lighted office room. In comparison to the LBP classifier, the Haar cascade classifier is a little slower, but still capable of consistently classify pictures into faces and nonfaces. In order to detect the left and right eye in an image, the corresponding Haar cascade eye detectors are utilized. To find the best face recognizer for the described application, several tests are considered. In particular, these tested the accuracy, training duration, model size and recognition speed over number of components (Eigenfaces/Fisherfaces) and classes in the training set as well as their total recognition accuracy over the number of training images (see Figure 10).
It can be seen, that, after applying the preprocessing steps, the Eigenface method constantly outperforms the Fisherface method, even though. For three and more faces, the LBPH is slightly more accurate than the Eigenface method. The test also shows that the training duration of Eigenfaces is for relatively small numbers of components slightly lower than the Fisherface method, while this changes for larger numbers of components. The training duration of LBPH is by far longer, especially when using larger radii and more neighbors. This also shows in comparing them at model sizes. LBPHs' model sizes are a lot larger than PCAs' and LDAs' models, which show the same magnitude, even though the LDA models are smaller. The Fisherfaces' recognition speeds are, starting at equally many classes and components, slightly higher than the Eigenfaces'. Before that point, both are almost identically high and about five times bigger than the LBPH's recognition speed. As seen in Figure 11, the Fisherfaces' recognition accuracy reaches its maximum approximately at the number of components being equal to the number of classes. This can be explained with the amount of class in the dataset. As described in section IV.B, the number of components is limited to C-1, which means adding additional components would not increase the recognition accuracy. Similar to this, the Eigenfaces' recognition accuracy reaches its maximum at R, the number of images in the training set (see section IV.A). Using five classes times ten images (minus two for testing), this limit is reached at approximately 40 components. Figure 11: Recognition accuracy vs. number of components These performance measures suggest using an Eigenface classifier because of its higher recognition accuracy and else similar properties. Further tests has shown, that using approximately 30 components, show best recognition results. This matches literatures suggestions [13,23,36]. [4] recommends not using the first three Eigenfaces to achieve even higher accuracies. Unfortunately this option is not supported by OpenCV [24].
In the final application, images for a face dataset have been collected. There, four individuals with approximately 1400 images are considered. The system implemented in OpenCV runs in real-time, providing approximately 15-18 frames per second.

VI. FUSION OF THE LOCALIZATION RESULTS FOR SPEAKER IDENTIFICATION
A fusion of the localization results of sound source localization and face recognition is able to enhance the reliability of a speaker detection system and enables it to identify the speaker. This can be used in smart rooms, improved speaker tracking for videoconferencing or applications for ambient assisted living [7,8,29]. Furthermore, an extension towards gesture recognition for human machine interaction is possible [7,10]. Therefore, a speaker identification algorithm is developed and introduced in this contribution.
In order to track and identify speakers, reliable sound source localization and face recognition are necessary. The sound source and therefore the potential speaker, is located by finding the color map's maximum. The face recognition provides a specific localization as well, but sometimes there are some uncertainties, which have to be eliminated. These could be falsely positive detected and recognized faces or wrongly classified individuals. To overcome these problems, three recognized faces are compared. If all of them are classified as the same person, the result is shown at the face's new position. Another possibility is that the recognized face matches one of the two previously identified individuals and is in approximately the same localization area as the currently detected. This results in a certain recognition at the currently classified face's position, too. If both options do not apply, no recognition result is being displayed and it is ignored like there never has been a face in the image. Using both, the speaker and face localization results, the overall localization and identification can be achieved as seen in Figure 12. The estimated outcome can differentiate between following: no result, identified speaker, face only, unknown speaker or loudspeaker, loudspeaker and known face at two different positions. Whenever possible, the localization position is chosen to the face's location, because it is well known that optical tracking algorithms have better spatial resolution than acoustic localization techniques [8].  Figure 13 shows a possible output image of the implemented speaker identification system. There, a speaker's face and sound source location are detected and merged for a more precise localization and identification of the sound source. The red circle marks the speaker's approximate position. Figure 13: Result of the speaker identification system VII. CONCLUSION This contribution briefly explains the basics of an acoustic camera and shows, why it makes sense to use a double ring array with an odd number of microphones. Additionally, it gives an overview of the implemented sound source localization methods. It can be shown, that a combination of SRP and SRP-PHAT algorithms is desirable for speech localization. Furthermore, this contribution gives an introduction to face detection and recognition methods. It is shown that Haar cascade classifiers outperform local binary pattern classifiers, when detecting faces. Similar, it is pointed out that for a face recognition system, Eigenfaces should be preferred to Fisherfaces and local binary pattern histograms. Finally, an algorithm for the fusion of localization results is introduced. This combines sound sources localization and face detection to identify speakers reliably.


Rapid Finger Motion Tracking on Low-Power Mobile Environments for Large Screen Interaction

Motion and gesture are garnering significant interest as the sizes of screens are getting larger. To provide lightweight finger motion tracking on low-power mobile environments, we propose an approach that breaks down the stereotypes of camera view points. By directing the camera view angle towards the ceiling, the proposed approach can reduce the problem complexity incurred by complicated background environments. Though this change incurs poor lighting conditions for image processing, by clustering and tracking the fragmented motion blobs from the motion image of the saturation channel, rapid finger motion can be tracked efficiently with low computational load. We successfully implemented and tested the proposed approach on a low-power mobile device with a 1.5 GHz mobile processor and a low specification camera with a capture rate of under 15 fps.

Introduction
As screens are getting larger, multi-modal interaction including gesture and voice is becoming a significant research topic due to their intuitiveness and convenience in order to handle more information. Since the introduction of the RGB-D sensor, hand gesture and motion tracking have become rapidly studied fields [3,6,7] among multi-modal interfaces.
However, the approaches adopting RGB-D sensors belong to an area that is struggling to secure popularity due to the requirement of specialized depth sensors for the motion interface. On the other hand, with the growth of the computer vision technology, recognizing hand posture using RGB sensors is also seeing enhanced accuracy [2,4,5]; although this approach still has problems dealing with complex backgrounds, and working in low-power environments. In addition, both RGB and RGB-D based approaches have restrictions on the distance between the human and the sensor if the sensor is installed on a screen.
In order to overcome such limitations, in this paper we propose an approach that tracks rapid finger motion in low-power mobile environments. By breaking down the stereotypes of camera view points, we successfully show the potential of using low-power mobile devices as remote motion capturing devices. In addition, by focusing on finger motion rather than gesture we can successfully implement the tracking system on low-power mobile devices, and configure it as a remote interface for large screen interaction, even in the case of blurry images.

Proposed Approach
Compared to traditional motion capturing systems on large screens, we separate the capturing sensor from the screen and hand its functionality to remote mobile devices that are more broadly distributed than specialized sensors like RGB-D. By directing the camera view angle towards the ceiling from around 20cm under the hand we can expect a relatively simple and static background compared to camera perspectives used in traditional motion capture systems. Challenging problems in the proposed system include minimizing lighting effects and detecting fast hand movement within a short distance between the hand and the mobile camera. In this paper we focus on the motion of two fingers to extract more robust finger position than using a single finger, and to open the possibility of extending the method to include recognition of gestures such as finger picking in the future. The overview of the proposed approach is illustrated in Fig. 1. We will describe each of the issues and how we resolved them in the following sub-sections.

Acquiring Image and Noise Removal
The reasons for most traditional motion capture systems used for large screens being directed towards the human body (along with potentially complex background information) include the camera's mount position, and the spatial relationship between the user and the camera. In the proposed system, by separating the capture function from the screen body to a remote mobile device, we can configure the spatial relationship between the user and the camera more flexibly. However, by directing the camera to the ceiling, the proposed system can not avoid the effect of lighting which is one of the most challenging problems in the field of computer vision. In order to minimize the effects of lighting, we capture an HSV (Hue, Saturation, Value) image rather than RGB image, and use only the saturation channel, which is known to be robust to changes in lighting. Fig. 2 (a) shows the captured saturation channel. From the saturation channel, we calculate the motion difference after the noise removal using a morphological operation.

Rapid Finger Motion Detection and Tracking
Though saturation is relatively robust to changes in lighting compared to RGB, it is not invariant; as such, there can occur fragments of motion blobs in the differential image. Many motion detection approaches adopted Motion History Image (MHI) [1] which accumulates image differentials over time to acquire robust motion blobs. However, in the case of low-power environments it is hard to use MHI because the frame rate of the device is too low to capture rapid finger motion. In the proposed approach, we emphasize differential image by horizontal morphological operation accepting fragmented blobs, as shown in Fig. 2 (b). In order to find the fingers in the fragmented motion image, we divide the entire image horizontally and compress each horizontal cell with the column-wise OR operation; Fig. 2 (c) shows these horizontal cells. By checking the number of black and white transitions from the top compressed cell to the bottom, we can find the first cell that includes the fingers. Based on the detected horizontal cell, we set up wider ROI (Region Of Interest) and extract the image blobs of each fragmented finger motion, as seen in Fig. 2 (d). In order to classify the motion blobs as left finger or right finger, we first find the upper-left and upper-right corners of the region spanned by the motion blobs (indicated by the light-blue dots in Fig. 2 (d)) by comparing the distance of the corners of each bounding box with upper-left and upper-right points of the ROI (indicated by red dots). Using the upper-left and upper-right corners of the region of motion, we can then classify each blob by the corner to which its center is nearest; illustrated in Fig. 2 (e). In order to track the motion robustly, we track each finger's motion over time and verify it by comparing size, aspect ratio and overlapping region with the previous frame.

Performance
We implemented the proposed approach in a low-power mobile environment with a 1.5 GHz mobile processor and a low specification camera with a capture rate of under 15 fps. Under these conditions, our method successfully detected and tracked rapid finger motion at 21 ms/frame in 320 × 240 resolution.

Conclusion
In this paper, we proposed a rapid finger motion tracking methodology on lowpower mobile environments for large screen interaction. In order to reduce the complexity of the problem caused by complicated and dynamic background conditions, we separated the capture function from the screen to a remote mobile device and changed the camera view, accepting fragmented motion caused by rapid finger motion and poor lighting conditions. By detecting finger motion efficiently from the fragmented motion blobs, we have successfully tracked rapid finger motion in a low-power environment. We implemented and verified the proposed approach on a low-power mobile device, which demonstrated real-time performance. In our future work, we will extend this approach to recognize finger motion gestures to provide more flexible interaction.


Viscous Fingering: A Topological Visual Analytic Approach

We present a methodology to analyze and visualize an ensemble of finite pointset method (FPM) simulations that model the viscous fingering process of salt solutions inside water. In course of the simulations the solutions form structures with increased salt concentration value, called viscous fingers. These structures are of primary interest to domain scientists since it is not deterministic when and where viscous fingers appear and how they evolve. To explore the aleatoric uncertainty embedded in the simulations we analyze an ensemble of simulation runs which differ due to stochastic effects. To detect and track the viscous fingers we derive a voxel volume for each simulation where fingers are identified as subvolumes that satisfy geometrical and topological constraints. Properties and the evolution of fingers are illustrated through tracking graphs that visualize when fingers form, dissolve, merge, and split. We provide multiple linked views to compare, browse, and analyze the ensemble in real-time.

: Particle system with 1,900,000 points rendered with a diverging color map filtered by salt concentration values: all particles (left), particles above concentration 100 (middle), and particles above concentration 200 (right).

Introduction
Viscous fingering is an instability phenomenon that occurs in porous media at the interface between two fluids of distinct viscosity. In particular, it appears when a less viscous fluid is injected within a more viscous one. This process is prominent in many fields of science and engineering [6,8]; including geology, hydrology, and chromatography. For instance, it plays a key role in the extraction process of oil reservoirs, and in the process of groundwater contamination. The phenomenon is characterized by the formation of a distinctive pattern called viscous fingers (Fig. 1 left). In particular, the geometrical evolution of these patterns provides good indications about the evolution of the penetration of the less viscous fluid. Thus, identifying, tracking, and analyzing the viscous fingers is of first importance for understanding the penetration process. Viscous fingering can be decomposed into three major regimes: 1. Launch: Initially, the interface between the two fluids is approximately planar. Then, the less viscous fluid starts to penetrate the more viscous one when sufficient injection force is applied to it. Viscous fingers of low and uniform amplitude start to appear in this phase.
2. Expansion: Once the reaction is launched, the difference of pressure between the two fluids tends to favor an acceleration of the penetration speed for the areas where the less viscous fluid penetrates the more viscous fluid the most.
3. Termination: Optionally, depending on the characteristics of the media and the fluids, the two fluids can eventually mix together in a termination state, where the finger pattern has completely disappeared, after finger structures merged and dissolved.
We analyze an ensemble of finite pointset method (FPM) simulations that model the viscous fingering process of salt solutions inside water. FPM simulations represent a medium as a numerical point cloud where each point stores properties of the medium such as salt concentration and velocity. The ensemble was provided for the 2016 scientific visualization contest [10] and consists of 50 FPM simulations that differ due to the aleatoric uncertainty of the mixing process. Each member simulates a cylinder that is filled with pure water and contains an infinite salt supply at its top. As soon as the salt mixes with the water, the resulting solutions sink down to the bottom as they have a higher density than the surrounding water. Fig. 2 shows the points of a FPM simulation where each particle has an associated salt concentration value. In course of the simulations the solutions form structures with increased salt concentration value called viscous fingers. However, it is not clear when and where viscous fingers appear and how they evolve. Furthermore, the accuracy of the mean solution rate and properties of the fingers also depend on the point cloud resolution. The ensemble consists of 50 simulations at three resolution levels, i.e., 250k, 650k, and 1900k particles. Each simulation run provides around 100 timesteps.   3 shows the outline of our approach that is an extension of the methodologies proposed by Aldrich et al. [2] and Favelier et al. [7]. First, we represent the computational domain as a voxel grid and estimate the concentration density distribution for each simulation based on the their particle clouds. Next, we identify viscous fingers as voxels with increased salt concentration density, i.e., superlevel sets in the voxel grid that exceed a certain concentration level (Fig. 3a). To identify individual fingers, noise, and the salt supply, we partition the superlevel sets according to the spatial Reeb graph (Fig. 3c). Each connected component of the partitioned, filtered superlevel set is a single viscous finger (Fig. 3d). To track fingers over time we test their respective volumes for spatial overlaps in time. Finally, we visualize the evolution of the fingers with an interactive tracking graph ( Fig. 1  right). Edges of the graph represent the lifespan of individual fingers, and nodes of the graph indicate the time when fingers appear, disappear, merge, and split. The tracking graph is linked to a 3D rendering of the particle data and a direct volume rendered image of the viscous fingers ( Fig. 1 left). Users can select fingers in either view to highlight parts of the voxel volume or branches of the tree to interactively explore the evolution of fingers.
We also store for each ensemble member computed metrics-such as size, integral, and bounding box-and direct volume rendered images of the viscous fingers in an image-based ensemble database that was built using the Cinema specification [1]. This database can be searched to compare specific fingers and their evolution. Furthermore, we can use the interactive tracking graph to query the database by selecting important timesteps or fingers in the graph. This enables users to effectively browse, compare, and examine the ensemble in real-time.

Related Work
There exists a huge body of work that demonstrates that tracking graphs are ideally suited to visualize the evolution and relation of components over time [9]. To utilize these methods we first compute a volumetric representation of the point clouds via an adaptation of Shepard's kernel method [13] for three-dimensional domains. Several approaches have been proposed to identify and track components of such time-varying scalar fields. Components can be superlevel sets within the domain that exceed a certain level [12,15,5,14,11], or subdomains that fulfill geometrical and topological constraints [3,4]. In our hybrid approach we compute superlevel set components and then filter components based on geometric and topological constraints. Subsequently, the extracted components can be tracked by testing for spatial overlaps in time [12,5,11,15], or applying methods of topological persistence [3,4,16]. Specifically, to track the three-dimensional finger structures we extend the tracking algorithm proposed by Lukasczyk et al. [12] for two-dimensional grid cells. We visualize the tracked fingers, their evolution, and the relationship between individual fingers via tracking graphs [3, 5, 12, 16]. Bremer et al. [5] use tracking graphs to visualize burning cells in large-scale combustion simulations where cells are defined as areas exceeding a fuel consumption rate threshold and are tracked via spatial overlaps. Building on Bremer et al. [5], Widanagamaachchi et al. [16] propose a dynamic version of tracking graphs capable of quick graph layout updates for varying thresholds.
To interactively browse and analyze the ensemble we utilize a Cinema [1] database. In general, such as database consists of pre-rendered images and pre-computed analysis results. The images and analysis results are derived for a range of parameter settings, e.g., camera angle, time, or voxel resolution. Later, a database-viewer can browse through the images and analysis results in real-time as if the results were computed at run time, i.e., emulating interactive data exploration.

Data Preprocessing
The ensemble consists of three sets of simulation runs (one set per resolution). Each run is represented as a time-varying particle system carrying salt concentration where the number, positions, and salt concentration values of particles vary over time. As proposed by Favelier et al. [7], we use an adaptation of Shepard's kernel method [13] to compute for each timestep a volumetric interpolation of the particles onto a 128 3 grid. However, our approach allows the use of other interpolation methods. Next, the regular grid is clipped to the original shape of the domain (a cylinder) to avoid boundary effects of the interpolation. The resulting tetrahedral meshes are then passed as input to our data analysis pipeline. The advantage of this preprocessing is threefold: 1. it brings continuity to the data, which enables the usage of a larger class of analysis algorithms; 2. the input for our analysis pipeline is independent of the particle cloud resolution, and; 3. it allows for interactive renderings of each timestep based on direct volume rendering, isosurfaces, or color mapped clipped views. Fig. 3 shows the outline of our detection algorithm. To identify fingers within the voxel data, we start by only considering voxels that exceed a user-specified salt concentration level (Fig. 3a). From our experience, we found that a salt concentration threshold of 10 gave consistent results along timesteps and across simulation runs. We denote with L the largest connected component of this superlevel set. To filter voxels which are considered noise or belong to the salt supply we derive the spatial Reeb graph (Fig. 3b) with the method proposed by Lukasczyk et al. [12]. This graph represents the evolution of the fingers in vertical direction (z-direction) for one single timestep, i.e., the graph is a topological skeleton of the fingers and as such captures when fingers form, dissolve, merge, and break into parts while sinking. We use this graph to partition the voxel volume in z-direction and subsequently check for each component of the graph whether its associated voxels have to be filtered out. Specifically, components are identified as part of the salt supply when their xy-bounding boxes are roughly the size of the can and are located at the top of the domain. Valid fingers also have to consists of more than 8 voxels and be at least two grid cells high to filter out small bubbles and bumps, respectively. As shown in Fig. 4, although the spatial Reeb graph detects bumps at t 0 and t 2 , they only get classified as fingers when they become significantly large enough in timesteps t 1 and t 4 , respectively. Those criteria, however, can easily be adjusted according to the interest of the domain scientists. After filtering noise and the salt supply (Fig. 3c), we determine the connected components within the remaining voxels and assign each voxel group a unique, positive integer ID across the entire simulation (Fig. 3d). Voxels not belonging to any group have the voxel group ID -1. Furthermore, for each voxel group (finger) we store its min/max value, bounding box, number of voxels, center of mass, and the salt concentration integral over its volume.

12
Physical Modeling for Virtual Manufacturing Systems and Processes Fig. 4: We use the spatial Reeb graph to filter out noise, the salt supply, and fingers that do not exceed a specified size threshold.
In order to implement this approach efficiently, the concentration values are represented as a 1D floating-point array and the component ID of each voxel is stored in an additional integer array of the same size called the component matrix (CM). Both arrays are passed as data textures to a WebGL shader that performs direct volume rendering. The shader renders for a given salt concentration the iso-surfaces and colors them according to the CM (Fig. 1 left). To enhance spatial perception we also render iso-lines in z-direction, an outline of the computational domain, and use screen-space ambient occlusion. It is also possible to view and superimpose the spatial Reeb graph and the underlying particle dataset. The rendering of the particle cloud is interactive, i.e., the analyst can filter particles by salt concentration, adjust the particle sizes, and choose the color map. Moreover, the analyst can examine the mapping between the particles and their associated fingers by color coding them according to the CM. By varying the concentration threshold for the finger detection, the analyst is able to examine the impact of said threshold in real time and choose a threshold that seems appropriate to partition the particles into fingers. The computation of the CM, the direct volume rendering (512 2 pixels), and the rendering of the particles (512 2 pixels for up to one million particles) is done at interactive framerates even on devices with integrated graphic cards.

Finger Tracking
In order to track fingers over time, we extended the method of Lukasczyk et al. [12] to derive tracking graphs for 3D structures. Hence, nodes in the graph represent 3D volumes (the fingers) and nodes of adjacent timesteps are connected by an edge iff their respective 3D volumes overlap. Critical nodes within this graph are nodes with either no previous nodes (birth nodes), no successor nodes (death nodes), multiple predecessor nodes (merge nodes), or multiple successor nodes (split nodes). Hence, a birth node represents the first time a finger is being detected, and a death node marks the last timestep before a finger dissolves. Merge nodes represent the event of multiple fingers joining, while split nodes indicate that fingers break apart. Fig. 5 shows the tracking graph for an entire simulation run where some timesteps are linked to a corresponding DVR image of the fingers. In this illustration of the tracking graph, the x-axis represents time, and the y-axis is used to minimize edge crossings to clearly show the evolution of fingers over time by reducing occlusion. Birth and death nodes are shown as discs, while merge and split nodes are shown as diamonds. Each finger is shown in a different color where colors are passed on from the largest predecessor to the largest successor. This makes it easy to follow the history of fingers and links them to the 3D rendering. The thickness of an edge encodes some metric of its corresponding finger, e.g., its average z-position, number of voxels, integral, or length. As demonstrated in Fig. 5, the tracking graph effectively illustrates the evolution of the fingers and provides a compact overview of an ensemble member. We also provide an alternative visualization that uses the y-axis of the tracking graph to show a metric of the fingers. For example, in Fig. 6 we use the y-axis to show the average z-position of a finger for the same simulation run as in Fig. 5. This visualization enables users to effectively understand the spatial relationship between fingers. However, this representation is highly occluded, especially when the number of tracked fingers increases. Therefore, it is possible to click on an edge to highlight all components that are connected to that edge, and fade out the rest. The user can easily switch between both representations and illustrated metrics. By clicking on nodes and edges the user can view information about the associated finger and highlight relevant volumes in the 3D rendering.  (top and bottom), the tracking graph (middle) effectively summarizes the evolution of the viscous fingers (shown in different colors), i.e., the time they appear (discs), disappear (discs), merge (diamonds), and split (diamonds). The x-axis represents time, the y-axis is used to minimize edge crossings, and the line thickness encodes finger volume.

Ensemble Summarization and Comparison
We compare and contrast the ensemble members by leveraging a visualization database framework constructed using the ParaView Cinema concept [1]. In this image-based approach, we create a database by analyzing and rendering the dataset with the methods described in the previous sections. Specifically, we store images of the particle dataset and the iso-surfaces for varying camera angles, concentration thresholds, and timesteps. To explore the image database, we implemented a web-based viewer that stores and displays the image database (Fig. 7). Each image in the viewer represents a single simulation; analyzed for a fixed concentration threshold. The camera of an image can be rotated by dragging it with the mouse, emulating real time volumetric interaction. The user can also use a slider to advance the current view over time. For comparison purposes, we also allow the user to sync both the view and the simulation time across all images. The user can browse through the massive image database by using an intuitive query interface. Each image in the database is linked to the parameters used to generate that image and statistical metrics we calculate based on the tracking graphs such as the average density, number, lifetime, or length of fingers. Once the user has discovered a dataset of interest, they can explore the simulation in the detail view (Fig. 1).
Using the analysis results stored in the Cinema database we can compute simulation summarizations and comparisons. For instance, to analyze the evolution of the descent of the dissolving salt we can query the minimum z-coordinates of the stored bounding boxes through time (Fig. 8). Three major observations come out of this plot. First, the descent seems linear with time after timestep 10. Second, all runs behave identically in the early timesteps and only little differences can be observed towards the end of the process. The differences that do occur increase linearly over time. Third, the 3 regimes (Launch, Expansion, Termination) are clearly visible: Launch: From timestep 0 to 10, the dissolving salt remains at the top of the domain; Expansion: From timestep 10 to 50-60, the dissolving salt travels down the domain; Termination: All runs hit the bottom of the domain in the same time-frame (timestep 50 to 60), after which a mixing of the two fluids can be observed. To confirm the bounds of the launch regime, we inspect the evolution of the average salt concentration of the dissolving salt L (Fig. 9). These curves exhibit a clear inflection pattern in the early stages of the simulation, with a clear local maximum at timestep 10 that marks the end of the launch phase.  To study the impact of the input data particle resolution on the characterization of the viscous fingering process, we visualize global summaries provided by global time-varying statistics for all runs in Fig. 10. This figure plots, for the dissolving salt L, the evolution over time of the minimum z-coordinate, the average salt concentration, the volume, and the average velocity of fingers. In particular, each of the three resolutions is represented with a distinct color. These three sets of curves exhibit similar global behaviors: a salt descent according to a linear slope (Fig. 10a), a linear decrease of salt concentration in the expansion regime (Fig. 10b), and a linear increase in both volume (Fig. 10c) and velocity (Fig. 10d). However, each resolution can be easily distinguished from the others as curves of the same color tend to cluster, with only minimal overlap with the other colors. This indicates clear distinctions between resolutions.

Conclusion and Future Work
We presented a topological data analysis framework adapted to the special case of analyzing finger structures produced across an ensemble of particle based fluid mixing simulations. Our analysis pipeline is based on a volumetric interpolation of the particle data, followed by finger identification and tracking based on geometric and topological constraints. Tracking these finger structures over time by spatial overlaps, we enable the evaluation of relevant time-varying statistics.
We applied our analysis framework on the data-sets provided in the scientific visualization contest 2016 [10]. From this analysis, the viscous fingering process can be characterized into three distinct regimes: launch, expansion, and termination. For these simulations, we were able to bound the occurrences for each of these phenomena. The interactive exploration features of our framework also enabled the quantitative validation of three characterizations within the expansion regime: larger finger structures tend to grow rapidly while smaller structures eventually dissolve. These observations provide insight to the consistent decrease of the numbers of fingers observed as the simulations progress. Additionally, we observed an unexpected behavior, that we named "detachment", where some fingers suddenly detach from the base of concentrated salt. These detachments typically occurred where the finger structures developed quickly, penetrated deeply, and consisted of a relatively low salt concentration. Our inter-resolution analysis revealed that simulation runs of distinct resolutions had globally similar behavior. However, later expansion starts, faster penetration rates, and fewer viscous fingers were observed as the resolution decreased (fewer particles represented the system as a whole). This comparison shows that the similarity to the highest resolution runs (according to the above criteria) decreases with the resolution, suggesting a convergence of the simulation code for increasing resolutions.
In the future we would like to validate our analysis framework with domain experts. Moreover, we would like to investigate the design of quantitative measures for a deeper evaluation of detachments, as we suspect they may have an impact on the efficiency of the penetration. Physical Modeling for Virtual Manufacturing Systems and Processes


s2fieldsofstudy:
Computer Science
Mathematics

text:
Automated Recurrence Analysis for Almost-Linear Expected-Runtime Bounds

We consider the problem of developing automated techniques for solving recurrence relations to aid the expected-runtime analysis of programs. Several classical textbook algorithms have quite efficient expected-runtime complexity, whereas the corresponding worst-case bounds are either inefficient (e.g., QUICK-SORT), or completely ineffective (e.g., COUPON-COLLECTOR). Since the main focus of expected-runtime analysis is to obtain efficient bounds, we consider bounds that are either logarithmic, linear, or almost-linear ($\mathcal{O}(\log n)$, $\mathcal{O}(n)$, $\mathcal{O}(n\cdot\log n)$, respectively, where n represents the input size). Our main contribution is an efficient (simple linear-time algorithm) sound approach for deriving such expected-runtime bounds for the analysis of recurrence relations induced by randomized algorithms. Our approach can infer the asymptotically optimal expected-runtime bounds for recurrences of classical randomized algorithms, including RANDOMIZED-SEARCH, QUICK-SORT, QUICK-SELECT, COUPONCOLLECTOR, where the worst-case bounds are either inefficient (such as linear as compared to logarithmic of expected-runtime, or quadratic as compared to linear or almost-linear of expected-runtime), or ineffective. We have implemented our approach, and the experimental results show that we obtain the bounds efficiently for the recurrences of various classical algorithms.

Introduction
Static analysis for quantitative bounds. Static analysis of programs aims to reason about programs without running them. The most basic properties for static analysis are qualitative properties, such as safety, termination, liveness, that for every trace of a program gives a Yes or No answer (such as assertion violation or not, termination or not). However, recent interest in analysis of resource-constrained systems, such as embedded systems, as well as for performance analysis, quantitative performance characteristics are necessary. For example, the qualitative problem of termination asks whether a given program always terminates, whereas the quantitative problem asks to obtain precise bounds on the number of steps, and is thus a more challenging problem. Hence the problem of automatically reasoning about resource bounds (such as time complexity bounds) of programs is both of significant theoretical as well as practical interest.
Worst-case bounds. The worst-case analysis of programs is the fundamental problem in computer science, which is the basis of algorithms and complexity theory. However, tion 2.4 for examples), our sound approach can obtain the asymptotically optimal expected-runtime bounds for the recurrences. In all the cases above, either the worst-case bounds (i) do not exist (e.g., COUPON-COLLECTOR), or (ii) are quadratic when the expected-runtime bounds are linear or almost-linear (e.g., QUICK-SELECT, QUICK-SORT); or (iii) are linear when the expected-runtime bounds are logarithmic (e.g., RANDOMIZED-SEARCH). Thus in cases where the worst-case bounds are either not applicable, or grossly overestimate the expectedruntime bounds, our technique is both efficient (linear-time) and can infer the optimal bounds. 3. Implementation. Finally, we have implemented our approach, and we present experimental results on the classical examples to show that we can efficiently achieve the automated expected-runtime analysis of randomized recurrence relations.
Novelty and technical contribution. The key novelty of our approach is an automated method to analyze recurrences arising from randomized recursive programs, which are not covered by Master theorem. Our approach is based on a guess-and-check technique. We show that by over-approximating terms in a recurrence relation through integral and Taylor's expansion, we can soundly infer logarithmic, linear and almost-linear bounds using simple comparison between leading terms of pseudo-polynomials.

Recurrence Relations
We present our mini specification language for recurrence relations for expectedruntime analysis. The language is designed to capture running time of recursive randomized algorithms which involve (i) only one function call whose expected-runtime complexity is to be determined, (ii) at most two integer parameters, and (iii) involve randomized-selection or divide-and-conquer techniques. We present our language separately for the univariate and bivariate cases. In the sequel, we denote by N, N 0 , Z, and R the sets of all positive integers, non-negative integers, integers, and real numbers, respectively.

Univariate Randomized Recurrences
Below we define the notion of univariate randomized recurrence relations. First, we introduce the notion of univariate recurrence expressions. Since we only consider single recursive function call, we use 'T' to represent the (only) function call. We also use 'n' to represent the only parameter in the function declaration.
Univariate recurrence expressions. The syntax of univariate recurrence expressions e is generated by the following grammar: e ::= c | n | ln n | n · ln n | 1 n | T (n − 1) | T n 2 | T n 2 | n−1 j=1 T(j) n | 1 n · n−1 j=⌈n/2⌉ T(j) + n−1 j=⌊n/2⌋ T(j) | c · e | e + e where c ∈ [1, ∞) and ln( ) represents the natural logarithm function with base e. Informally, T(n) is the (expected) running time of a recursive randomized program which involves only one recursive routine indicated by T and only one parameter indicated by n. Then each T( )-term in the grammar has a direct algorithmic meaning: -T (n − 1) may mean a recursion to a sub-array with length decremented by one; -T n 2 and T n 2 may mean a recursion related to a divide-and-conquer technique; finally, n−1 j=1 T(j) n and 1 n · n−1 j=⌈ n 2 ⌉ T(j) + n−1 j=⌊ n 2 ⌋ T(j) may mean a recursion related to a randomized selection of an array index.
Substitution. Consider a function h : N → R and univariate recurrence expression e. The substitution function, denoted by Subst(e, h), is the function from N into R such that the value for n is obtained by evaluation through substituting h for T and n for n in e, respectively. Moreover, if e does not involve the appearance of 'T', then we use the abbreviation Subst(e) i.e., omit h. For example, (i) if e = n + T(n − 1), and h : n → n · log n, then Subst(e, h) is the function n → n + (n − 1) · log(n − 1), and (ii) if e = 2 · n, then Subst(e) is n → 2n.
Univariate recurrence relation. A univariate recurrence relation G = (eq 1 , eq 2 ) is a pair of equalities as follows: where c ∈ (0, ∞) and e is a univariate recurrence expression. For a univariate recurrence relation G the evaluation sequence Eval(G) is as follows: Eval(G)(1) = c, and for n ≥ 2, given Eval(G)(i) for 1 ≤ i < n, for the value Eval(G)(n) we evaluate the expression Subst(e, Eval(G)), since in e the parameter n always decreases and is thus well-defined.
Finite vs infinite solution. Note that the above description gives a computational procedure to compute Eval(G) for any finite n, in linear time in n through dynamic programming. The interesting question is to algorithmically analyze the infinite behavior.
The function T G is unique and explicitly defined as follows: (1) Base Step. T G (1) := c; and (2) Recursive Step. T G (n) := Subst(e, T G )(n) for all n ≥ 2. The interesting algorithmic question is to reason about the asymptotic infinite behaviour of T G .

Motivating Classical Examples
In this section we present several classical examples of randomized programs whose recurrence relations belong to the class of univariate recurrence relations described in Section 2.1. We put details of pseudocode and how to derive the recurrence relations in this section in Appendix A. Moreover in all cases the base step is T(1) = 1, hence we discuss the recursive case.

Example 1 (RANDOMIZED-SEARCH
We note that the worst-case complexity for this algorithm is Θ(n).

⊓ ⊔
Example 2 (QUICK-SORT). Consider the QUICK-SORT algorithm [16, Chapter 7]. The recurrence relation for this example is: where T(n) represents the maximal expected execution time where n is the array length and the execution time of pivoting is represented by 2 · n. We note that the worst-case complexity for this algorithm is Θ(n 2 ). ⊓ ⊔

Example 3 (QUICK-SELECT).
Consider the QUICK-SELECT algorithm (cf. [16, Chapter 9]). The recurrence relation for this example is We note that the worst-case complexity for this algorithm is Θ(n 2 ). ⊓ ⊔

Example 4 (DIAMETER-COMPUTATION).
Consider the DIAMETER-COMPUTATION algorithm (cf. [40,Chapter 9]) to compute the diameter of an input finite set S of three-dimensional points. Depending on Eucledian or L 1 metric we obtain two different recurrence relations. For Eucledian we have the following relation: and for L 1 metric we have the following relation: We note that the worst-case complexity for this algorithm is as follows: for Euclidean metric it is Θ(n 2 · log n) and for the L 1 metric it is Θ(n 2 ). ⊓ ⊔ Example 5 (Sorting with QUICK-SELECT). Consider a sorting algorithm which selects the median through the QUICK-SELECT algorithm. The recurrence relation is directly obtained as follows: where T * ( ) is an upper bound on the expected running time of QUICK-SELECT (cf. Example 3). We note that the worst-case complexity for this algorithm is Θ(n 2 ). ⊓ ⊔

Separable Bivariate Randomized Recurrences
We consider a generalization of the univariate recurrence relations to a class of bivariate recurrence relations called separable bivariate recurrence relations. Similar to the univariate situation, we use 'T' to represent the (only) function call and 'n', 'm' to represent namely the two integer parameters.
Separable Bivariate Recurrence Expressions. The syntax of separable bivariate recurrence expressions is illustrated by e, h and b as follows: The differences are that (i) we have two independent parameters n, m, (ii) e now represents an expression composed of only T-terms, and (iii) h (resp. b) represents arithmetic expressions for n (resp. for m). This class of separable bivariate recurrence expressions (often for brevity bivariate recurrence expressions) stresses a dominant role on m and a minor role on n, and is intended to model randomized algorithms where some parameter (to be represented by n) does not change value.
Substitution. The notion of substitution is similar to the univariate case. Consider a function h : N × N → R, and a bivariate recurrence expression e. The substitution function, denoted by Subst(e, h), is the function from N × N into R such that Subst(e, h)(n, m) is the real number evaluated through substituting h, n, m for T, n, m, respectively. The substitution for h, b is defined in a similar way, with the difference that they both induce a univariate function.
Bivariate recurrence relations. We consider bivariate recurrence relations G = (eq 1 , eq 2 ), which consists of two equalities of the following form: where c ∈ (0, ∞) and e, h, b are from the grammar above.
Solution to bivariate recurrence relations. The evaluation of bivariate recurrence relation is similar to the univariate case. Similar to the univariate case, the unique solution T G : N × N → R to a recurrence relation G taking the form (8)

Motivating Classical Examples
In this section we present two classical examples of randomized algorithms where the randomized recurrence relations are bivariate. We put the detailed illustration for this two examples in Appendix B.
Example 6 (COUPON-COLLECTOR). Consider the COUPON-COLLECTOR problem [40,Chapter 3] with n different types of coupons (n ∈ N). The randomized process proceeds in rounds: at each round, a coupon is collected uniformly at random from the coupon types the rounds continue until all the n types of coupons are collected. We model the rounds as a recurrence relation with two variables n, m, where n represents the total number of coupon types and m represents the remaining number of uncollected coupon types. The recurrence relation is as follows: T(n, 1) = n · 1; T(n, m) = n/m + T(n, m − 1) where T(n, m) is the expected number of rounds. We note that the worst-case complexity for this process is ∞. ⊓ ⊔ Example 7 (CHANNEL-CONFLICT RESOLUTION). We consider two network scenarios in which n clients are trying to get access to a network channel. This problem is also called the RESOURCE-CONTENTION RESOLUTION [36,Chapter 13]. In this problem, if more than one client tries to access the channel, then no client can access it, and if exactly one client requests access to the channel, then the request is granted. In the distributed setting, the clients do not share any information. In this scenario, in each round, every client requests an access to the channel with probability 1 n . Then for this scenario, we obtain an over-approximating recurrence relation T(n, 1) = n · 1; T(n, m) = (n · e)/m + T(n, m − 1) for the expected rounds until which every client gets at least one access to the channel. In the concurrent setting, the clients share one variable, which is the number of clients which has not yet been granted access. Also in this scenario, once a client gets an access the client does not request for access again. For this scenario, we obtain an overapproximating recurrence relation T(n, 1) = 1 · 1; T(n, m) = 1 · e + T(n, m − 1) We also note that the worst-case complexity for both the scenarios is ∞. ⊓ ⊔

Expected-Runtime Analysis
We focus on synthesizing logarithmic, linear, and almost-linear asymptotic bounds for recurrence relations. Our goal is to decide and synthesize asymptotic bounds in the simple form: d · f + g, f ∈ {ln n, n, n · ln n} . Informally, f is the major term for time complexity, d is the coefficient of f to be synthesized, and g is the time complexity for the base case specified in (1) or (8).
Univariate Case: The algorithmic problem in univariate case is as follows: -Input: a univariate recurrence relation G taking the form (1) and an expression f ∈ {ln n, n, n · ln n}. -Output: Decision problem. Output "yes" if T G ∈ O(Subst(f)), and "fail" otherwise. -Output: Quantitative problem. A positive real number d such that for all n ≥ 1, or "fail" otherwise, where c is from (1).

Remark 1.
First note that while in the problem description we consider the form f part of input for simplicity, since there are only three possibilites we can simply enumerate them, and thus have only the recurrence relation as input. Second, in the algorithmic problem above, w.l.o.g, we consider that every e in (1) or (8) involves at least one T( )term and one non-T( )-term; this is natural since (i) for algorithms with recursion at least one T( )-term should be present for the recursive call and at least one non-T( )term for non-recursive base step. ⊓ ⊔ Bivariate Case: The bivariate-case problem is an extension of the univariate one, and hence the problem definitions are similar, and we present them succinctly below.
-Input: a bivariate recurrence relation G taking the form (8) and an expression f (similar to the univariate case). -Output: Decision problem. Output "yes" if T G ∈ O(Subst(f)), and "fail" otherwise; -Output: Quantitative problem. A positive real number d such that T G (n, m) ≤ d · Subst(f)(n, m) + c · Subst(h)(n) for all n, m ≥ 1, or "fail" otherwise, where c, h are from (8). Note that in the expression above the term b does not appear as it can be captured with f itself.
Recall that in the above algorithmic problems obtaining the finite behaviour of the recurrence relations is easy (through evaluation of the recurrences using dynamic programming), and the interesting aspect is to decide the asymptotic infinite behaviour.

The Synthesis Algorithm
In this section, we present our algorithms to synthesize asymptotic bounds for randomized recurrence relations.
Main ideas. The main idea is as follows. Consider as input a recurrence relation taking the form (1) and an univariate recurrence expression f ∈ {ln n, n, n · ln n} which specifies the desired asymptotic bound. We first define the standard notion of a guessand-check function which provides a sound approach for asymptotic bound. Based on the guess-and-check function, our algorithm executes the following steps for the univariate case.
1. First, the algorithm sets up a scalar variable d and then constructs the template h to be n → d · Subst(f)(n) + c for a univariate guess-and-check function.

Second, the algorithm computes an over-approximation OvAp(e, h) of Subst(e, h)
such that the over-approximation OvAp(e, h) will involve terms from n k , ln ℓ n (for k, ℓ ∈ N 0 ) only. Note that k, ℓ may be greater than 1, so the above expressions are not necessarily linear (they can be quadratic or cubic for example). 3. Finally, the algorithm synthesizes a value for d such that OvAp(e, h)(n) ≤ h(n) for all n ≥ 2 through truncation of [2, ∞) ∩ N into a finite range and a limit behaviour analysis (towards ∞).
Our algorithm for bivariate cases is a reduction to the univariate case.
Guess-and-Check functions. We follow the standard guess-and-check technique to solve simple recurrence relations. Below we first fix a univariate recurrence relation G taking the form (1). By an easy induction on n (starting from the N specified in Definition 1) we obtain Theorem 1.
Definition 1 (Univariate Guess-and-Check Functions). Let G be a univariate recurrence relation taking the form (1). A function h : N → R is a guess-and-check function Theorem 1 (Guess-and-Check, Univariate Case). If a function h : N → R is a guess-and-check function for a univariate recurrence relation G taking the form (1), We do not explicitly present the definition for guess-and-check functions in the bivariate case, since we will present a reduction of the analysis of separable bivariate recurrence relations to that of the univariate ones (cf. Section 4.2).
Overapproximations for Recurrence Expressions. We now develop tight overapproximations for logarithmic terms. In principle, we use Taylor's Theorem to approximate logarithmic terms such as ln (n − 1), ln ⌊ n 2 ⌋, and integral to approximate summations of logarithmic terms. All the results below are technical and depends on basic calculus (the detailed proofs are in the Appendix C).
Proposition 3. For all natural number n ≥ 2: Note that Proposition 3 is non-trivial since it approximates summation of reciprocal and logarithmic terms up to a constant deviation. For example, one may approximate n−1 j=1 ln j directly by n 1 ln x dx, but this approximation deviates up to a logarithmic term from Proposition 3. From Proposition 3, we establish a tight approximation for summation of logarithmic or reciprocal terms. Example 8. Consider the summation n−1 j=⌈ n 2 ⌉ ln j + n−1 j=⌊ n 2 ⌋ ln j (n ≥ 4). By Proposition 3, we can over-approximate it as x dx = n · ln n − n − ln n 2 + 1. By using Proposition 1, the above expression is roughly n · ln n − (1 − ln 2) · n + 1 2 · ln n + 0.6672 + 1 2·n (for details see Appendix C). ⊓ ⊔ Remark 2. Although we do approximation for terms related to only almost-linear bounds, Proposition 3 can be extended to logarithmic bounds with higher degree (e.g., n 3 ln n) since integration of such bounds can be obtained in closed forms. ⊓ ⊔

Algorithm for Univariate Recurrence Relations
We present our algorithm to synthesize a guess-and-check function in form (12) for univariate recurrence relations. We present our algorithm in two steps. First, we present the decision version, and then we present the quantitative version that synthesizes the associated constant. The two key aspects are over-approximation and use of pseudopolynomials, and we start with over-approximation. We relegate some technical details to Appendix D.
Definition 2 (Overapproximation). Let f ∈ {ln n, n, n · ln n}. Consider a univariate recurrence expression g, constants d and c, and the function h = d · Subst(f) + c. We define the over-approximation function, denoted OvAp(g, h), recursively as follows.
-Base Step A. If g is one of the following: c ′ , n, ln n, n · ln n, 1 n , then OvAp(g, h) := Subst(g).
-Base Step B. If g is a single term which involves T, then we define OvAp(g, h) from over-approximations Proposition 1-3. In details, OvAp(g, h) is obtained from Subst(g, h) by first over-approximating any summation through Proposition 3 (i.e., through those Γ ( ) functions defined below Proposition 3), then over-approximating any ln (n − 1), n 2 , n 2 , ln n 2 , ln n 2 by Proposition 1 and Proposition 2. The details of the important over-approximations are illustrated explicitly in Table 1.
Example 9. Consider the recurrence relation for Sherwood's RANDOMIZED-SEARCH (cf. (2)). Choose f = ln n and then the template h becomes n → d·ln n+1. From Example 8, we have that the over-approximation for 6 + 1 n · n−1 j=⌈ n ⊓ ⊔ Table 1. Illustration for Definition 2 where the notations are given in the top-left corner.
Over-approximation f, T-term Over-approximation n, e1 n − 1 n · ln n, e1 n · ln n − ln n − 1 + 1 n n, e2 n 2 n · ln n, e2 n · ln n, e5 Remark 3. Since integrations of the form x k ln l x dx can be calculated in closed forms (cf. Remark 2), Table 1 can be extended to logarithmic expressions with higher order, e.g., n 2 ln n.

⊓ ⊔
Pseudo-polynomials. Our next step is to define the notion of (univariate) pseudopolynomials which extends normal polynomials with logarithm. This notion is crucial to handle inductive arguments in the definition of (univariate) guess-and-check functions.

Definition 3 (Univariate Pseudo-polynomials).
A univariate pseudo-polynomial (w.r.t logarithm) is a function p : N → R such that there exist non-negative integers k, ℓ ∈ N 0 and real numbers a i , b i 's such that for all n ∈ N, W.l.o.g, we consider that in the form (13), it holds that (i) a 2 k +b 2 ℓ = 0, (ii) either a k = 0 or k = 0, and (iii) similarly either b ℓ = 0 or ℓ = 0.
Degree of pseudo-polynomials. Given a univariate pseudo-polynomial p in the form (13), we define the degree deg(p) of p by: deg(p) = k + 1 2 if k ≥ ℓ and a k = 0 and ℓ otherwise. Intuitively, if the term with highest degree involves logarithm, then we increase the degree by 1/2, else it is the power of the highest degree term.
Leading term p. The leading term p of a pseudo-polynomial p in the form (13) is a function p : N → R defined as follows: p(n) = a k · n k · ln n if k ≥ ℓ and a k = 0; and b ℓ · n ℓ otherwise ; for all n ∈ N. Furthermore, we define C p to be the (only) coefficient of p.
With the notion of pseudo-polynomials, the inductive argument of guess-and-check functions can be soundly transformed into an inequality between pseudo-polynomials.
Lemma 1. Let f ∈ {ln n, n, n · ln n} and c be a constant. For all univariate recurrence expressions g, there exists pseudo-polynomials p and q such that coefficients (i.e., a i , b i 's in (13)) of q are all non-negative, C q > 0 and the following assertion holds: for all d > 0 and for all n ≥ 2, Remark 4. In the above lemma, though we only refer to existence of pseudopolynomials p and q, they can actually be computed in linear time, because p and q are obtained by simple rearrangements of terms from OvAp(g, h) and h, respectively.
Example 10. Let us continue with Sherwood's RANDOMIZED-SEARCH. Again choose h = d · ln n + 1. From Example 9, we obtain that for every n ≥ 4, the inequality resulting from over-approximation and the inductive argument of guess-and-check holds for sufficiently large n. The following proposition provides a sufficient and necessary condition for checking whether d · p(n) ≥ q(n) holds for sufficiently large n.
Proposition 4. Let p, q be pseudo-polynomials such that C q > 0 and all coefficients of q are non-negative. Then there exists a real number d > 0 such that d · p(n) ≥ q(n) for sufficiently large n iff deg(p) ≥ deg(q) and C p > 0.
Note that by Definition 1 and the special form (12) for univariate guess-and-check functions, a function in form (12) needs only to satisfy the inductive argument in order to be a univariate guess-and-check function: once a value for d is synthesized for a sufficiently large N , one can scale the value so that the base condition is also satisfied. Thus from the sufficiency of Proposition 4, our decision algorithm that checks the existence of some guess-and-check function in form (12) is presented below. Below we fix an input univariate recurrence relation G taking the form (1) and an input expression f ∈ {ln n, n, n · ln n}.
Algorithm UniDec: Our algorithm, namely UniDec, for the decision problem of the univariate case, has the following steps.
4. Coefficient Checking. The algorithm examines cases on C p . If C p > 0 and deg(p) ≥ deg(q), then algorithm outputs "yes" meaning that "there exists a univariate guess-and-check function"; otherwise, the algorithm outputs "fail".
Theorem 2 (Soundness for UniDec). If UniDec outputs "yes", then there exists a univariate guess-and-check function in form (12) for the inputs G and f. The algorithm is a linear-time algorithm in the size of the input recurrence relation.
Example 11. Consider Sherwood's RANDOMIZED-SEARCH recurrence relation (cf. (2)) and f = ln n as the input. As illustrated in Example 9 and Example 10, the algorithm asserts that the asymptotic behaviour is O(ln n). ⊓ ⊔ Remark 5. From the tightness of our over-approximation (up to only constant deviation) and the sufficiency and necessity of Proposition 4, the UniDec algorithm can handle a large class of univariate recurrence relations. Moreover, the algorithm is quite simple and efficient (linear-time). However, we do not know whether our approach is complete. We suspect that there is certain intricate recurrence relations that will make our approach fail.
Analysis of examples of Section 2.2. Our algorithm can decide the following optimal bounds for the examples of Section 2.2.
1. For Example 1 we obtain an O(log n) bound (recall worst-case bound is Θ(n)).
In all cases above, our algorithm decides the asymptotically optimal bounds for the expected-runtime analysis, whereas the worst-case analysis grossly over-estimate the expected-runtime bounds.
Quantitative bounds. Above we have already established that our linear-time decision algorithm can establish the asymptotically optimal bounds for the recurrence relations of several classical algorithms. We now take the next step to obtain even explicit quantitative bounds, i.e., to synthesize the associated constants with the asymptotic complexity. To tackle these situations, we derive a following proposition which gives explicitly a threshold for "sufficiently large numbers". We first explicitly constructs a threshold for "sufficiently large numbers". Then we show in Proposition 5 that N ǫ,p,q is indeed what we need. Definition 4 (Threshold N ǫ,p,q for Sufficiently Large Numbers). Let p, q be two univariate pseudo-polynomials Then given any ǫ ∈ (0, 1), the number N ǫ,p,q is defined as the smallest natural number such that both x, y (defined below) is smaller than ǫ: where 1 deg(p)=deg(q) equals 1 when deg(p) = deg(q) and 0 otherwise.
With Proposition 5, we describe our algorithm UniSynth which outputs explicitly a value for d (in (12)) if UniDec outputs yes. Below we fix an input univariate recurrence relation G taking the form (1) and an input expression f ∈ {ln n, n, n · ln n}. Moreover, the algorithm takes ǫ > 0 as another input, which is basically a parameter to choose the threshold for finite behaviour. For example, smaller ǫ leads to large threshold, and vice-versa. Thus we provide a flexible algorithm as the threshold can be varied with the choice of ǫ.
Algorithm UniSynth: Our algorithm for the quantitative problem has the following steps: 1. Calling UniDec. The algorithm calls UniDec, and if it returns "fail", then return "fail", otherwise execute the following steps. Obtain the following inequality d · p(n) ≥ q(n) (n ∈ N) from the transformation step of UniDec.  Table 2). ⊓ ⊔

Algorithm for Bivariate Recurrence Relations
In this part, we present our results for the separable bivariate recurrence relations. The key idea is to use separability to reduce the problem to univariate recurrence relations. There are two key steps which we describe below.
Step 1. The first step is to reduce a separable bivariate recurrence relation to a univariate one.
Definition 5 (From G to Uni(G)). Let G be a separable bivariate recurrence relation taking the form (8). The univariate recurrence relation Uni(G) from G is defined by eliminating any occurrence of n and replacing any occurrence of h with 1.
Informally, Uni(G) is obtained from G by simply eliminating the roles of h, n. The following example illustrates the situation for COUPON-COLLECTOR example.
Example 13. Consider G to be the recurrence relation (9) for COUPON-COLLECTOR example. Then Uni(G) is as follows: T(n) = 1 n + T(n − 1) and T(1) = 1. ⊓ ⊔ Step 2. The second step is to establish the relationship between T G and T Uni(G) , which is handled by the following proposition, whose proof is an easy induction on m.

Proposition 6. For any separable bivariate recurrence relation G taking the form (8), the solution T G is equal to (n, m) → Subst(h)(n) · T Uni(G) (m).
Description of the Algorithm. With Proposition 6, the algorithm for separable bivariate recurrence relations is straightforward: simply compute Uni(G) for G and then call the algorithms for univariate case presented in Section 4.1.
Analysis of examples in Section 2.4. Our algorithm can decide the following optimal bounds for the examples of Section 2.4.
1. For Example 6 we obtain an O(n · log m) bound, whereas the worst-case bound is ∞.

For Example 7 we obtain an O(n · log m) bound for distributed setting and O(m)
bound for concurrent setting, whereas the worst-case bounds are both ∞.
Note that for all our examples, m ≤ n, and thus we obtain O(n · log n) and O(n) upper bounds for expected-runtime analysis, which are the asymptotically optimal bounds. In all cases above, the worst-case analysis is completely ineffective as the worstcase bounds are infinite. Moreover, consider Example 7, where the optimal number of rounds is n (i.e., one process every round, which centralized Round-Robin schemes can achieve). The randomized algorithm, with one shared variable, is a decentralized algorithm that achieves O(n) expected number of rounds (i.e., the optimal asymptotic expected-runtime complexity).
We consider the classical examples illustrated in Section 2.2 and Section 2.4. In Table 2 for experimental results we consider the following recurrence relations G: R.-SEAR. corresponds to the recurrence relation (2) for Example 1; Q.-SORT corresponds to the recurrence relation (3) for Example 2; Q.-SELECT corresponds to the recurrence relation (4) for Example 3; DIAM. A (resp. DIAM. B) corresponds to the recurrence relation (5) (resp. the recurrence relation (6)) for Example 4; SORT-SEL. corresponds to recurrence relation (7) for Example 5, where we use the result from setting ǫ = 0.01 in Q.-SELECT; COUPON corresponds to the recurrence relation (9) for Example 6; RES. A (resp. RES. B) corresponds to the recurrence relation (10) (resp. the recurrence relation (11)) for Example 7.
In the table, f specifies the input asymptotic bound, ǫ and Dec is the input which specifies either we use algorithm UniDec or the synthesis algorithm UniSynth with the given ǫ value, and d gives the value synthesized w.r.t the given ǫ ( for yes). We describe d 100 below. We need approximation for constants such as e and ln 2, and use the interval [2.7182, 2.7183] (resp., [0.6931, 0.6932]) for tight approximation of e (resp., ln 2).
The value d 100 . For our synthesis algorithm we obtain the value d. The optimal value of the associated constant with the asymptotic bound, denoted d * , is defined as follows. (1)). Then the sequence d z is increasing in z, and its limit is the optimal constant, i.e., d * = lim z→∞ d z . We consider d 100 as a lower bound on d * to compare against the value of d we synthesize. In other words, d 100 is the minimal value such that (12) holds for 1 ≤ n ≤ 100, whereas for d * it must hold for all n, and hence d * ≥ d 100 . Our experimental results show that the d values we synthesize for ǫ = 0.01 is quite close to the optimal value. We performed our experiments on Intel(R) Core(TM) i7-4510U CPU, 2.00GHz, 8GB RAM. All numbers in Table 2 are over-approximated up to 10 −3 , and the running time of all experiments are less than 0.02 seconds. From Table 2, we can see that optimal d are effectively over-approximated. For example, for QUICK-SORT (Eq. (3)) (i.e, Q.-SORT in the table), our algorithm detects d = 4.051 and the optimal one lies somewhere in [3.172, 4.051]. The experimental results show that we obtain the results extremely efficiently (less than 1/50-th of a second). For further details see Table 3 in Appendix E.

Related Work
Automated program analysis is a very important problem with a long tradition [46]. The following works consider various approaches for automated worstcase bounds [28,29,30,31,32,35,34,26,5,44] for amortized analysis, and the SPEED project [23,24,22] for non-linear bounds using abstract interpretation. All these works focus on the worst-case analysis, and do not consider expected-runtime analysis.
Our main contribution is automated analysis of recurrence relations. Approaches for recurrence relations have also been considered in the literature. Wegbreit [46] [2,3,4] transforms Java bytecode into recurrence relations and solves them through ranking functions. Moreover, The PURRS tool [6] addresses finite linear recurrences (with bounded summation), and some restricted linear infinite recurrence relations (with unbounded summation). Our approach is quite different because we consider analyzing recurrence relations arising from randomized algorithms and expected-runtime analysis through over-approximation of unbounded summations through integrals, whereas previous approaches either consider recurrence relations for worst-case bounds or combinatorial structures, or use generating functions or difference equations to solve the recurrence relations.
Proof rules have also been considered for recursive (probabilistic) programs in [25,33,41], but these methods cannot be automated and require manual proofs.

Conclusion
In this work we considered efficient algorithms for automated analysis of randomized recurrences for logarithmic, linear, and almost-linear bounds. Our work gives rise to a number of interesting questions. First, an interesting theoretical direction of future work would be to consider more general randomized recurrence relations (such as with more than two variables, or interaction between the variables). While the above problem is of theoretical interest, most interesting examples are already captured in our class of randomized recurrence relations as mentioned above. Another interesting practical direction would be automated techniques to derive recurrence relations from randomized recursive programs. Let T : N → N be the function such that for any n ∈ N, we have T (n) is the supremum of the expected execution times upon all inputs (ar, i, j) with j − i + 1 = n. We derive a recurrence relation for T as follows. Let n ∈ N and (ar, i, j), d be any input such that n = j − i + 1. We clarify two cases below:

ar[j] ≤ d or d < ar[i].
In both cases, we have T (1) = 1. In Case 1, we deduce from the pseudo-code in Fig. 1 that for all n ≥ 2, where the maximum ranges over all ℓ * := k * −i+1's. In Case 2, similarly we deduce that Thus a preliminary version G ′ of the recurrence relation is T(1) = 1 and for all n ≥ 2. Let T ′ : N → R be the unique solution to G ′ . Then from the fact that is attained at ℓ * = n 2 for all n ≥ 2. Then G ′ is transformed into our final recurrence relation as follows: We note that the worst-case complexity for this algorithm is Θ(n From the pseudo-code, the following recurrence relation is easily obtained: quicksort(ar, m + 1, j) ; end i f end i f } r e t u r n d 9 : e l s e 1 0 : r e t u r n diameter(S ′ ) end i f end i f } Fig. 4. DIAMETER-COMPUTATION for all n ≥ 2, the following recurrence relation is obtained: To fit our univariate recurrence expression, we use over-approximation, and the final recurrence relation for this example is We note that the worst-case complexity for this algorithm is Θ(n 2 ). ⊓ ⊔

Example 4.[DIAMETER-COMPUTATION]
Consider the DIAMETER-COMPUTATION algorithm (cf. [40,Chapter 9]) to compute the diameter of an input finite set S of threedimensional points. A pseudo-code to implement this is depicted in Fig. 4 Depending on Eucledian or L 1 metric we obtain two different recurrence relations. For Eucledian we have the following relation: T(n) = 2 + n + 2 · n · ln n + ( n−1 j=1 T(j))/n; with the execution time for lines 5-6 being taken to be 2 · n · ln n, and and for L 1 metric we have the following relation: T(j))/n with the execution time for lines 5-6 being taken to be 2 · n. We note that the worst-case complexity for this algorithm is as follows: for Euclidean metric it is Θ(n 2 · log n) and for the L 1 metric it is Θ(n 2 ).

:
i f ( m + 1 < j ) 6 : sortbyselect(ar, m + 1, j) ; end i f end i f }  N). The randomized process proceeds in rounds: at each round, a coupon is collected uniformly at random from the coupon types (i.e., each coupon type is collected with probability 1 n ); and the rounds continue until all the n types of coupons are collected. We model the rounds as a recurrence relation with two variables n, m, where n represents the total number of coupon types and m represents the remaining number of uncollected coupon types. The recurrence relation is as follows: T(n, 1) = n · 1; T(n, m) = n/m + T(n, m − 1) where T(n, m) is the expected number of rounds, n m represents the expected number of rounds to collect a new (i.e., not-yet-collected) coupon type when there are still m type of coupons to be collected, and n (for T(n, 1)) represents the expected number of rounds to collect a new coupon type when there is only one new coupon type to be collected. We note that the worst-case complexity for this process is ∞.

Example 7.[CHANNEL-CONFLICT RESOLUTION]
We consider two network scenarios in which n clients are trying to get access to a network channel. This problem is also called the RESOURCE-CONTENTION RESOLUTION [36,Chapter 13]. In this problem, if more than one client tries to access the channel, then no client can access it, and if exactly one client requests access to the channel, then the request is granted. While centralized deterministic algorithms exist (such as Round-Robin) for the problem, to be implemented in a distributed or concurrent setting, randomized algorithms are necessary.
Distributed setting. In the distributed setting, the clients do not share any information.
In this scenario, in each round, every client requests an access to the channel with probability 1 n . We are interested in the expected number of rounds until every client gets at least one access to the channel. At each round, let m be the number of clients who have not got any access. Then the probability that a new client (from the m clients) gets the access is m· 1 n ·(1− 1 n ) n−1 . Thus, the expected rounds that a new client gets the access is n m · 1 (1− 1 n ) n−1 . Since the sequence (1 − 1 n ) n−1 n∈N converges decreasingly to 1 e when n → ∞, this expected time is no greater than e · n m . Then for this scenario, we obtain an over-approximating recurrence relation T(n, 1) = n · 1; T(n, m) = (n · e)/m + T(n, m − 1) for the expected rounds until which every client gets at least one access to the channel. Note that in this setting no client has any information about any other client.
Concurrent setting. In the concurrent setting, the clients share one variable, which is the number of clients which has not yet been granted access. Also in this scenario, once a client gets an access the client does not request for access again. Moreover, the shared variable represents the number of clients m that have not yet got access. In this case, in reach round a client that has not access to the channel yet, requests access to the channel with probability 1 m . Then the probability that a new client gets the access becomes m · 1 m · (1 − 1 m ) m−1 . It follows that the expected time that a new client gets the access becomes 1 (1− 1 m ) m−1 which is smaller than e. Then for this scenario, we obtain an over-approximating recurrence relation T(n, 1) = 1 · 1; T(n, m) = 1 · e + T(n, m − 1) We also note that the worst-case complexity for both is ∞. ⊓ ⊔

C Details for Overapproximations
To prove results for overapproximations for recurrence expressions, we need the following well-known theorem.
We also recall that Moreover, we have the following result using integral-by-part technique and Newton-Leibniz Formula.
Lemma 2. For all a, b ∈ (0, ∞) such that a < b, the following assertions hold: Furthermore, we need the following simple lemmas. The following lemma provides a tight approximation for floored expressions, the proof of which is a simple case distinction between even and odd cases.
where ξ j,x is a real number in (j, j + x) obtained from Taylor's Theorem with Lagrange's Remainder. The first and fourth equalities come from the linear property of Riemann Integral; the second one follows from the variable substitution x ′ = x − j; the third one follows from Taylor's Theorem. Using the fact that ξ j,x ∈ (j, j + 1), one obtains that Then (14) follows from the facts that where in both situations we use the fact that 2 · j 2 ≤ 3 · j 3 for all j ∈ N.
Then we consider the difference where in both situations we use the fact that 12·j 2 ≤ 6·j 3 for all j ≥ 2. The inequalities above directly imply the inequalities in (15). Finally, we consider the difference n 1 x · ln x dx − n−1 j=m j · ln j .
Following similar approaches, we derive that for all natural numbers n ≥ 2, where ξ j,x ∈ (j, j + 1). Thus, one obtains that   16.000 0. 3 Table 3. Detailed experimental results where all running times (averaged over 5 runs) are less than 0.02 seconds (between 0.01 and 0.02 seconds).