Y Benezeth

S Bobbia

K Nakamura

R Gomez

J Dubois

Probabilistic Signal Quality Metric for Reduced Complexity Unsupervised Remote Photoplethysmography

Keywords: Biomedical monitoring, heart rate measurements, signal quality metric

Remote photoplethysmography (rPPG) is a recent technique for estimating heart rate by analyzing the pulsatility of skin hue using regular cameras. To determine the quality of the measurement, many existing methods are based on the signal-tonoise ratio (SNR) calculated in the frequency domain. However, the Fast Fourier Transform (FFT) operation is performed with a minimal complexity of O(n log n). Therefore, the use of this quality metric in an unsupervised rPPG framework in which this metric is estimated a large number of times will tend to greatly increase the complexity of the solution. In this paper, we propose a new probabilistic formulation of a cardiac signal quality index, with lower complexity, based on the Bayesian information criterion (BIC) that encapsulates the characteristic shape of the rPPG signal. The results of this study, obtained on a public database, have demonstrated that the proposed probabilistic metric outperforms the regular SNR metric with a lower computation complexity.

I. INTRODUCTION

Since 2010 and pioneering work of Poh et al. [START_REF] Poh | Non-contact automated cardiac pulse measurements using video imaging and blind source separation[END_REF], remote photoplethysmography has become a very active field of research. In recent years, several methods have been proposed to estimate the physiological signal with fewer and fewer environmental constraints. Thus, an expensive configuration is no longer necessary [START_REF] Verkruysse | Remote plethysmographicimaging using ambient light[END_REF], since a consumer grade webcam can be used to make the measurements. The basic principle of remote photoplethysmography stems from reflective photoplethysmography where the light reaching a camera is modulated by the blood pulsations of skin. The rhythmic beating of the heart results in the pulsating blood volume which in turn results in minute changes in the color of the skin which can be quantified using different signal processing techniques to generate a cardiac signal.

The problem of rPPG measurement has been attacked from different directions. In the initial research related to rPPG measurement, the blood volume pulse (BVP) signal was usually extracted from the RGB temporal traces using blind source separation (BSS) such as ICA [START_REF] Poh | Non-contact automated cardiac pulse measurements using video imaging and blind source separation[END_REF] or PCA [START_REF] Lewandowska | Measuring pulse rate with a webcam -a noncontact method for evaluating cardiac activity[END_REF] or more recently PVM [START_REF] Macwan | Periodic Variance Maximization using Generalized Eigenvalue Decomposition applied to Remote Photoplethysmography estimation[END_REF]. Another class of methods focuses on exploiting the physical characteristics of skin tissue and its distinctive interaction with light. Methods such as CHROM [START_REF] De Haan | Robust pulse rate from chrominance-based rppg[END_REF], POS [START_REF] Wang | Algorithmic principles of remote-ppg[END_REF] and PBV [START_REF] De Haan | Improved motion robustness of remote-ppg by using the blood volume pulse signature[END_REF] fall into this category and are based on a detailed Supported by the Conseil Régional de Bourgogne Franche Comté through the plan d'actions régional pour l'innovation (PARI) and the European Union through the PO FEDER-FSE Bourgogne 2014/2020 programs. formulation of models based on physiological properties of the skin. The principal advantages of these methods lies in their computational simplicity, owing to their analytic formulation.

As discussed in different review papers (e.g. [START_REF] Sun | Photoplethysmography revisited: from contact to noncontact, from point to imaging[END_REF], [START_REF] Mcduff | A Survey of Remote Optical Photoplethysmographic Imaging Methods[END_REF]), many rPPG methods implement a general pipeline-based framework: regions of interest (ROI) are detected and tracked over frames, RGB channels are then combined to estimate the pulse signal, which is filtered and analyzed to extract physiological parameters such as heart rate or respiration rate. This pipeline-based framework emphasizes the importance of the common first step of ROI segmentation. Usually, ROI segmentation is based on the result of classical face detection [START_REF] Poh | Non-contact automated cardiac pulse measurements using video imaging and blind source separation[END_REF] and tracking algorithms and possibly refined with skin pixel classification [START_REF] Wang | Unsupervised subject detection via remote ppg[END_REF], [START_REF] Macwan | Heart rate estimation using remote photoplethysmography with multi-objective optimization[END_REF]. As an alternative to this pipeline approach, data-driven methods exploit the pulse-signal as a feature to segment the ROI in an unsupervised manner using voxels [START_REF] Wang | Unsupervised subject detection via remote ppg[END_REF] or temporal superpixels [START_REF] Bobbia | Real-Time Temporal Superpixels for Unsupervised Remote Photoplethysmography[END_REF], [START_REF] Bobbia | Remote Photoplethysmography Based on Implicit Living Skin Tissue Segmentation[END_REF] video segmentation. They are called unsupervised rPPG methods to emphasize the difference with methods that require a trained classifier to determine the ROI.

For these so-called unsupervised techniques, the quality of the different ROIs is estimated in order to identify skin areas [START_REF] Wang | Unsupervised subject detection via remote ppg[END_REF] or to favor areas where information is predominant [START_REF] Bobbia | Remote Photoplethysmography Based on Implicit Living Skin Tissue Segmentation[END_REF], [START_REF] Bobbia | Unsupervised skin tissue segmentation for remote photoplethysmography[END_REF]. The signal-to-noise ratio (SNR), defined as the ratio of the power of the main pulsatile component and the power of background noise, is generally used for this purpose. However, because of the video decomposition in voxels or temporal superpixels, a large number of FFT which has a complexity of O(n log n), has to be calculated. Moreover, with this frequency domain index, the shape of the temporal signal is not considered. However, in the case of rPPG signals, the temporal form of the signal is characteristic and is an indication of signal quality.

From these observations, we propose a probabilistic framework to model the rPPG signals and derive a signal quality metric in the temporal domain. Interestingly, the model encapsulates the difference in duration between the growth and decay phases of the different stages of the cardiac cycle and has a complexity of O(n).

The rest of the paper is organized as follows. The unsupervised rPPG framework is described in section II. Then the probabilistic rPPG signal quality metric is described in section III. Finally, experiments and results are presented in section IV while conclusion is presented in section V.

II. UNSUPERVISED RPPG FRAMEWORK

The overview of the proposed method is shown in Fig. 1. The algorithm can be decomposed into four main steps. First, the input video frames are decomposed into several temporal superpixels using the IBIS method [START_REF] Bobbia | Real-Time Temporal Superpixels for Unsupervised Remote Photoplethysmography[END_REF]. The segmentation step is performed by implicitly identifying the superpixel boundaries. Hence, only a fraction of the image is used to perform the segmentation which reduces greatly the computational burden of the process. The segmentation quality remains comparable to state of the art methods while computational time is significantly reduced. Second, a tentative rPPG signal is estimated from each temporal superpixel. The term tentative means that the estimate can also be done on background superpixels. It is important to note here that this approach can be used with any existing rPPG method. For the sake of simplicity, the CHROM method [START_REF] De Haan | Robust pulse rate from chrominance-based rppg[END_REF] has been selected for this work. Indeed, this method consists of a simple linear combination of the RGB channels and is therefore very fast to calculate. In our case, this computational simplicity is a very interesting property because this operation is performed for all temporal superpixels. Finally, a signal quality measure, i.e. a pulsatility measure, is estimated for each tentative pulse signal to determine contributive signals and a weighted average of all the signals is computed where the weights are given by the pulsatility measure. Since the most pulsatile signals provide high quality information, i.e. higher weight, this method implicitly selects regions of interest that represents living skin tissue and favors areas where the pulse trace is more predominant. This method is employed using a sliding window and consecutive resultant signals are then overlapped.

III. HIDDEN MARKOV MODELS FOR RPPG SIGNAL

QUALITY MEASUREMENT

A. rPPG signal modeling

In order to measure the quality of rPPG signals and the weights in the fusion process, it is possible to estimate signalto-noise ratio (SNR) defined as the ratio of the power of the main pulsatile component and the power of background noise, computed in dB due to the wide dynamic range of the signals. This approach has been used successfully in previous works [START_REF] Bobbia | Remote Photoplethysmography Based on Implicit Living Skin Tissue Segmentation[END_REF], [START_REF] Bobbia | Unsupervised skin tissue segmentation for remote photoplethysmography[END_REF]. However, the estimation of the rPPG signal quality of all the tentative signals with a frequency-domain measurement, imposes the computation of a large number of FFT which has a complexity of O(n log n). Moreover, with the SNR, the shape of the temporal signal is not considered. However, in the case of rPPG signals, the temporal form of the signal is characteristic and is an indication of the quality of the information. Fig. 2 presents a cardiac pulse signal that highlights the systole and diastole phenomena as well as the dicrotic notch. However, it is often not possible to observe the dicrotic notch in many remote PPG scenarios. In this case, the measured cardiac signal has two distinctive phases: rapid growth corresponding to systole followed by a slow decrease corresponding to the diastole. This difference in duration between the growth and decay phases is not captured by signal quality metrics in the frequency domain (such as SNR).

In this paper, we use a probabilistic framework to model the rPPG signal and derive a signal quality metric in the temporal domain. This approach has been used successfully in BayesHeart [START_REF] Xiangmin | Bayesheart: A probabilistic approach for robust, low-latency heart rate monitoring on camera phones[END_REF] where hidden Markov models are used to model temporal patterns of the different stages of the cardiac cycle. The signal is modeled by an alternation of two states

B. Cardiac signal quality measure

The quality metric of the rPPG signal is then determined using the Bayesian Information Criterion B [START_REF] Schwarz | Estimating the dimension of a model[END_REF] which provides the likelihood of the observations O given the model parameters λ:

B λ = -2.ln(P (O|λ)) + k.(ln(n) -ln(2π)). ( 1
)
where n is the number of points in the signal. Fig. 4 shows the model likelihood, as defined in (1) and a regular SNR metric estimated on all tentative rPPG signals. Thus, it is possible to observe that in both cases the face is correctly selected and that the areas that contribute the most (cheekbones and forehead) have a better score. Finally, in order to determine the quality of the tentative pulse signals, we use two two-state models. The first model is trained on a set of rPPG signals and is denoted λ s . The second model is an equiprobable model and is denoted λ n which models a random signal. The quality index v is eventually defined by:

v = 1 - B λs B λn . ( 2
)
We have empirically observed that it is often better to keep only the x% best contributions. We therefore define ṽ which represents the set of v such that:

ṽ = {v, v > max i (v i ). x 100 }, (3)
with i the i th temporal superpixel and max i (v i ) the most important contribution over all temporal superpixels. Weights are then obtained normalizing ṽi with:

w i = ṽi i∈K ṽi , (4)
with K the number of temporal superpixels. Finally, the fusion is carried out using the hidden states H i (t) determined by the Viterbi method to obtain, for each temporal superpixel, the most probable sequence of states that led to the provided observations. This method acts as a filtering and helps to mitigate disturbances. The final pulse signal S(t) is constructed with the weighted average of the state sequences of the different temporal superpixels weighted by the weights w i :

S(t) = i∈K w i H i (t). (5)
IV. EXPERIMENTS AND RESULTS The in-house UBFC-RPPG database [START_REF] Bobbia | Unsupervised skin tissue segmentation for remote photoplethysmography[END_REF] is used to evaluate the proposed cardiac signal quality metric in the unsupervised rPPG framework. The REALISTIC subset of this database, used in this work, comprises of 46 videos where the subjects were required to play a time sensitive mathematical game in order to vary the heart rate and simultaneously emulate the scenario of the typical activity of using a computer. All the videos were taken under ambient light with limited illumination variations. The experimental setup with sample images is depicted in Fig. 5. The UBFC-RPPG database is made publicly available along with the ground truth data from the pulse oximeter for rPPG measurement analysis 1 .

The video frames were obtained with a custom C++ application using a Logitech C920 web camera placed at about 1m from the subject with a resolution of 640x480 in 8-bit uncompressed RGB format at approximately 30 frames per second. A CMS50E transmissive pulse oximeter was used to obtain the ground truth PPG data. The estimation of the model parameters is performed with the Baum-Welch method [START_REF] Baum | An Inequality and Associated Maximization Technique in Statistical Estimation for Probabilistic Functions of Markov Processes[END_REF] on the SIMPLE subset of the UBFC-RPPG dataset. This subset is composed of 7 simple videos in which the volunteers were asked to sit still. The transition matrix A and the emission matrix B are initialized with:

A = 0.8 0.2 0.2 0.8 , B = 0.6 0.1 0.3 0.1 0.3 0.1 0.6 0.1 . (6)
1 https://sites.google.com/view/ybenezeth/ubfcrppg Thus, the probability of transition from one state to the other is set to 20% while the probability of observation of the symbols are function of the shape of the signal. Since the most represented symbols with a rPPG signal are growth and decay, the emission probabilities are set to 60% for the observation of a growth in state 1 (i.e. growth) and 30% for a continuous decay in the state 1. These probabilities are reversed in state 2 (i.e. decay). The probabilities of observing a concave or convex curve variation are necessarily lower and are thus sett to 10%.

We give below the matrices obtained after the training:

A = 0.
Table I shows the heart rate estimation accuracy comparison between the proposed signal quality metric based on the probabilistic framework and the regular SNR metric for the fusion process. The evaluation metrics used are Mean Absolute Error (M AE), Root Mean Square Error (RM SE), Pearsons correlation coefficient (r) and P recision between heart rate calculated using the rPPG signal and the heart rate calculated using the ground truth PPG waveform. The P recision metric represents the percentage of estimations where the absolute error is under a threshold (2.5 or 5 bpm). In order to obtain these results, the CHROM method [2] was selected for its simplicity, we use a resolution of 150 temporal superpixels and x = 5%. Values in bold indicate the best result for each metric. It can first be noted that the average error M AE and the squared error RM SE are both lower with BIC. That means that the number of errors is lower and that the wrong esimates have values which are closer from the ground truth with BIC. However, the results are very close, and the heart rate estimation is very good in both cases using this dataset. Fig. 6 shows the robustness to motion of the two signal quality indices. The background is generally very efficiently eliminated by the SNR metric and the difference between the weights is larger. However, the movements and deformation of superpixels tend to bring out false signal information induced by the variations of superpixel boundaries. On the contrary, BIC quality metric can be more resilient to these disturbances.

V. CONCLUSION

In the present study, we have described and evaluated a new cardiac signal metric in an unsupervised rPPG framework in which the regions of interest segmentation is performed by implicitly selecting living skin tissue via their distinct pulsatility feature. Photoplethysmogram signals are estimated with the weighted fusion of several tentative rPPG signals computed on a set of temporal superpixels. The cardiac signal quality metric is based on a probabilistic framework where the characteristic shape of the temporal signal is considered in the hidden Markov model. It has been shown that the fusion of several rPPG signals estimated on several regions of interest generally improves the quality of the resulting cardiac signal (e.g. [START_REF] Bobbia | Unsupervised skin tissue segmentation for remote photoplethysmography[END_REF]). However, this approach can be very intensive in calculations because the operations performed on each area are repeated several times. Therefore, it is interesting to note that the complexity of the proposed metric is O(n) whereas the complexity of the FFT (necessary for the SNR) is O(n log n).

Based on a publicly available dataset of 46 videos, the results of this study have demonstrated that the probabilistic metric outperforms the regular SNR metric with a lower computation complexity.

With regard to future developments, it would be interesting to quantify the saving of computing time more precisely and to study the possibility of using a four-states model by considering the dicrotic notch for scenarios where the signal is of good quality (for example in controlled environments with professional cameras). Next, we also plan to examine the use of this measure to analyze other biomedical signals.

Fig. 1 .

 1 Fig. 1. Flowchart of the temporal superpixel-based heart rate estimation framework: (1) Input video stream is decomposed into temporal superpixels. (2) Tentative rPPG signal is extracted from each temporal superpixel. (3) A signal quality index is estimated for each ROI. (4) A weighted average of all the tentative rPPG signals is finally computed.

Fig. 2 .

 2 Fig. 2. Characteristic waveform of PPG with systolic and diastolic phases and dicrotic notches.

Fig. 3 .

 3 Fig. 3. Possible discrete observations.

Fig. 4 .

 4 Fig. 4. Comparison of SNR (left) and BIC (right) as a cardiac signal quality measure. Blue means low pulsatility measures and yellow/orange means high. K is the number of temporal superpixels.

Fig. 5 .

 5 Fig. 5. Experimental Setup (top) and sample images from the UBFC-RPPG database (bottom).

Fig. 6 .

 6 Fig. 6. Example of results obtained with SNR (left) and BIC (right).

TABLE I RESULTS

 I FOR HEART RATE ESTIMATION WITH THE PROPOSED BIC SIGNAL QUALITY METRIC COMPARED TO THE REGULAR SNR METRIC.

	Metrics	SNR	BIC
	Precision 5 bpm	0.940	0.975
	Precision 2.5bpm	0.899	0.951
	RMSE	4.00	2.41
	MAE	2.05	1.21
	r	0.83	0.82


Peixi Li
email: peixili@etu.u-bourgogne.fr

Yannick Benezeth

Keisuke Nakamura

Randy Gomez

Fan Yang

Model-based Region of Interest Segmentation for Remote Photoplethysmography

Keywords: Remote Photoplethysmography (rPPG), Heart Rate (HR), Region of Interest Segmentation

. We have shown that this modification in how the spatial averaging of the ROI pixels is calculated can significantly increase the final performance of heart rate estimate.

INTRODUCTION

The photoplethysmography (PPG) is an optical technique to non invasively detect the blood volume pulse (BVP). PPG sensors have two basic components, a LED light source and a photodetector. The light source illuminates the human tissue and the photodetector detects the light variation. Since the amount of light variation is proportional to blood intensity variation, the BVP can be obtained and the HR is calculated based on BVP.

The PPG is a low-cost and non-invasive technique, however it cannot work when motion and longterm monitoring are required. Therefore, remote PPG (rPPG) has been recently developed to address these drawbacks. The principle of rPPG is very similar with contact PPG. However, instead of using a LED light source and a specified photodetector, the rPPG can simply makes use of a web camera and ambient light. By detecting the variation of the light reflected by the face, the BVP and HR can be measured.

According to many previous studies [START_REF] Sun | Photoplethysmography revisited: from contact to noncontact, from point to imaging[END_REF][START_REF] Mcduff | A survey of remote optical photoplethysmographic imaging methods[END_REF], the rPPG methods for HR measurement share a common pipeline-based framework: ROI selection for each frame, RGB sig-nals selection and combination to get the pulse signal, and finally HR extraction from the filtered pulse signal. The ROI selection is the first step of the whole framework. The state of the art methods for ROI selection includes manual selection [START_REF] Verkruysse | Remote plethysmographic imaging using ambient light[END_REF], skin pixel classification [START_REF] Wang | A novel algorithm for remote photoplethysmography: Spatial subspace rotation[END_REF] and conventional face detection and tracking [START_REF] Poh | Advancements in noncontact, multiparameter physiological measurements using a webcam[END_REF].

Lots of research work has indicated that the ROI detection has significant effect on the HR measurement [START_REF] Bousefsaf | Continuous wavelet filtering on webcam photoplethysmographic signals to remotely assess the instantaneous heart rate[END_REF]. This is because the numerical error of the signals will increase if the number of the effective pixels, i.e. skin pixels, is too small. Conversely, the number of non-skin pixels also has a significant negative influence on signal quality. Moreover, it has been shown that rPPG signal is not distributed homogeneously on the face. Some skin regions contain more rPPG signal than others. For example, it has been shown that the cheeks and forehead contain much more information than other areas of the face [START_REF] Tur | Basal perfusion of the cutaneous microcirculation: measurements as a function of anatomic position[END_REF]. From this premise, some studies have suggested to simply select the cheeks and/or forehead as ROI [START_REF] Scalise | Heart rate measurement in neonatal patients using a webcamera[END_REF][START_REF] Lewandowska | Measuring pulse rate with a webcam a non-contact method for evaluating cardiac activity[END_REF]. This selection has the undeniable advantage of simplicity but it is also possible to weight the skin pixels by their relative importance. In doing so, we explicitly favor areas of the face where there is more information.

In this paper, we propose a model-based ROI segmentation that explicitly favors the most important facial regions. The model describing the spatial distribution of rPPG information on the face was trained with ten videos. This model is then used to weight the pixels during spatial averaging. This approach has been validated using our in-house publicly available video dataset [START_REF] Bobbia | Unsupervised skin tissue segmentation for remote photoplethysmography[END_REF], called UBFC-RPPG specifically geared towards research on rPPG techniques. We have shown that this modification in how the spatial averaging of the ROI pixels is calculated can significantly increase the final performance of heart rate estimate compared with other state of the art methods such as face detection, skin classification and landmarks detection.

In section 2, the state of the art methods are briefly introduced. The proposed ROI segmentation method is explained in details in section 3. The experiment is described in section 4. The conclusion is presented in section 5.

STATE OF THE ART OF ROI SEGMENTATION METHODS

Most ROI segmentation techniques are based on the result of classical face detection and tracking algorithms. The ROI is then possibly refined with skin pixel classification or more precise ROI definition based on a set of landmarks. In this section we present several state of the art ROI segmentation techniques. Since all of the video dataset are the recordings of heads and faces, the most straightforward method to detect the ROI is to use face detector and tracker (later called face). For this algorithm, the classical Viola-Jones face detector [START_REF] Viola | Rapid object detection using a boosted cascade of simple features[END_REF] and Kanade-Lucas-Tomasi tracking (Lucas et al., ) algorithm (cf. Fig 1(a)) can be used for implementation.

Since the rPPG information is only present on skin pixels, skin/non-skin classification (later called skin) is a popular improvement over the classical face detection and tracking. For instance, some researchers [START_REF] Macwan | Remote photoplethysmography measurement using constrained ica[END_REF] used this algorithm [START_REF] Conaire | Detector adaptation by maximising agreement between independent data sources[END_REF] in their rPPG research work. The skin detection algorithm is achieved by thresholding of a non-parametric histogram which is trained by manually classified skin/non-skin pixels. The significant advantage of this algorithm is that it works very fast since it is based on a Look-Up-Table (LUT). Fig. 1(b) is an example of this method.

The ROI detection can also be segmented by defining a facial contour with a set of landmarks (later called landmarks). For this algorithm the method proposed by Kazemi can be used for implementation [START_REF] Kazemi | One millisecond face alignment with an ensemble of regression trees[END_REF]. One example is shown in Fig. 1(c).

The rPPG signal is not distributed homogeneously on skin. Some skin regions contain more rPPG signal than others. For example, it has been shown that SNR of rPPG signals extracted from forehead or cheekbones are significantly higher than other face regions. This assertion has already been used by different ROI segmentation techniques. For example, in some work only the cheeks and forehead were selected [START_REF] Lewandowska | Measuring pulse rate with a webcam a non-contact method for evaluating cardiac activity[END_REF]. In a previous study, ROI segmentation, based on temporal superpixels, implicitly favors regions of interest where the pulse trace is more prominent [START_REF] Bobbia | Unsupervised skin tissue segmentation for remote photoplethysmography[END_REF]. However, this datadriven method is very sensitive to motion and errors in superpixel tracking induce incorrect segmentation.

MODEL BASED ROI SEGMENTATION

In this paper, we propose an effective technique for explicitly favoring certain areas of the face during the spatial averaging step of RGB pixels. The model that encapsulates the spatial distribution of rPPG information was trained using an in-house database of 10 videos recorded under very favorable conditions. For this experiment, we used a EO-23121C camera recording 1024 × 768 uncompressed images at 30 fps.

The average length of each video is about one minute. Subjects sat on an chair with back support. To make sure that the face is fixed in a specified position, we used a shelf and asked the volunteers to put the heads onto the shelf. Fig. 2 shows two sample images from the dataset. Then, the face sequence is aligned based on the location of the eyes. The video frames were filtered by a 25×25 averaging filter to decrease quantization noise. The rPPG signal is extracted using the chrominancebased method [START_REF] De Haan | Robust pulse rate from chrominance-based rppg[END_REF]. This method is very fast and the computational complexity is very low. It linearly combined the RGB channels by projecting them onto two orthogonal chrominance vectors:

X(t) = 3y R (t) -2y G (t), Y (t) = 1.5y R (t) + y G (t) -1.5y B (t).
(1)

Where y c (t) is the RGB signal after filtering, c ∈ {R, G, B} are the color channels, and X and Y are two orthogonal chrominance vectors. The pulse signal S is obtained by S(t) = X(t) -αY (t) where α = σ(X)/σ(Y ). This σ is the standard deviation calculated over time.

The SNR for each pixel is then estimated as the ratio of the area under the power spectrum curve in a region surrounding the maximum peak in the frequency spectrum, divided by the area under the curve in the rest of the frequency spectrum:

SNR ( x, y) = 10 × log 10 ( P signal P noise ) (2)
Fig. 3 presents the spatial distribution of the SNR values averaged for all the videos of this dataset. As expected, we can observe that cheeks and forehead have in average, higher SNR (about 6 dB) than other face locations (e.g. 0 dB for the chin).

This spatial map is then used during the spatial averaging of the pixels of the ROI. Instead of using equal weights for all pixels in the ROI, we use a weighted average where weights are defined based on our model.

The weights are calculated as:

ω (x,y) = a SNR (x,y) ∑ (x,y) a SNR (x,y) (3)
where a is a constant. With higher a, the weights of the region of higher SNR are higher.

The optimal selection of a is done empirically. Experiments and details are given in Section 4.

Eventually, RGB triplet is obtained for each frame with the following weighted average: Where I t (x, y) is the RGB value at time t of a pixel at location (x, y) and ω (x,y) is the corresponding weight of the pixel.

S(t) = ∑ x,y I t (x, y) × ω (x,y) (4)
It is interesting to note that this technique can be advantageously combined with all the rPPG methods that perform a spatial average of the ROI pixels.

EXPERIMENTS

The proposed model-based ROI segmentation (later called model) is experimentally compared with the state of the art methods, namely face detection (face), skin detection (skin) and landmarks detection (landmarks). This section presents the details of the experiments. First, we describe the video dataset and experimental setup. Then, we present the system framework. Evaluation metrics and experimental results are finally presented.

Video Dataset

All of the ROI selection methods were assessed by the performance of remote HR measurement. For the experiments, we used our in-house video datasets UBFC-rPPG [START_REF] Bobbia | Unsupervised skin tissue segmentation for remote photoplethysmography[END_REF] which have 46 videos. The volunteers were asked to sit at the distance of one meter from the web camera and were asked to play a time sensitive video game. Ambient light was used in the experiment to create diffuse reflections. The experimental set up is shown in Fig. 4. A C++ program was used to record the videos and synchronize the videos with the signals of the contact PPG sensor which was used as the ground truth. The web camera is a Logitech C920. The resolution of the video frame is 640x480 and the frame rate is 30/second. The format is 8-bit uncompressed RGB. The contact PPG sensor is a CMS50E transmissive pulse oximeter. The experimental dataset with the ground truth can be downloaded from our project webpage 1 . Some sample images are shown in Fig. 5.

System Framework

As explained previously, most methods share a common pipeline-based framework where ROI are first detected and tracked over frames, RGB channels are then combined to estimate the pulse signal, which is filtered and analyzed to extract heart rates. In this study, we used the same procedure to compare our ROI segmentation method with other segmentation techniques. First, RGB signals are normalized, detrended by smoothness priors approach and filtered by band-pass Butterworth filter. Then the RGB signal is selected and combined by chrominance based method (De Haan and Jeanne, 2013).

1 https://sites.google.com/view/ybenezeth/ubfcrppg With this pulse signal, we used Welch's method to obtain the periodogram with a sliding window of 20 seconds, and step size is one second. Identically, this process was used on the PPG signal of the contact sensor, which was utilized as the ground truth. The system framework is briefly shown in Fig. 6.

The Evaluation Metrics

Seven metrics were used to assess the performance of the ROI selection algorithms in the framework of HR measurements:

• Pearson Correlation Coefficient (R) is used to evaluate the correlation between the RPPG measurement and contact PPG measurement (ground truth).

• Mean Absolute Error (MAE) in beats per minute (bpm) is calculated as the absolute differences between HR measured from RPPG and HR measured from contact PPG signals which can be represented as |HR rPPG -HR PPG |.

• MAE5 is the MAE which discards all the outliers with an error larger than 5 bpm.

• Precision5 and Precision2.5 represent the percentage of estimations where the absolute errors is under a threshold (2.5 or 5 bpm). • Root mean square error (RMSE) is the square root of the average of squared differences between HR measured from RPPG and ground truth.

• Signal-to-Noise Ratio (SNR) is calculated as the ratio of the power of the main pulsatile component and the power of background noise, computed in dB due to the wide dynamic range.

Results

First, we present an experiment to empirically select the optimal scalar a of equation (3). To do this, several values for a were tested in the system framework. The result is shown in table 1. It can be observed that the ranking is different according to the different metrics. a = 2 provides a higher SNR while a larger value of a provides a better HR estimation (given by Correlation R, MAE and RMSE). For the rest of the experiments, we chose to favor the quality of the signal (given by the SNR) and will use a = 2. Second, the proposed model-based ROI segmentation method is compare with other state of the art methods. The average results are shown in table 2. It is obvious that the ROI segmentation algorithms have significant effects on the results. All the metrics are very different with different ROI selection methods. It is also very clear that the proposed algorithm performs much better than all the other algorithms. It has the lowest MAE5 and RMSE. It has the highest correlation R, which means the RPPG method matches the ground truth better than other methods. SNR metric is also the highest which indicates that the method offers the best quality of the signal. The skin method is actually very good, although it is slightly worse than the proposed algorithm.

CONCLUSION

A good ROI is very important for rPPG algorithms to work properly. The ROI must contain as much critically useful information as possible. Most of the state of the art methods focus on the improvement of the face and skin detection, and their objectives are mostly to get as many skin pixels as possible and to discard as many non-skin pixels as possible. However, the rPPG signals are not distributed homogeneously on the human face, so it is reasonable to consider putting weights on the pixels before other processing. In this paper, we presented a supervised model based on SNR weights to improve the ROI segmentation. This algorithm is compared with other algorithms by evaluating the performance of HR measurement. The experiments were done with a low-cost web camera and a contact PPG sensor as the ground truth. 56 videos were used in the experiments and 10 of them were used to create the SNR weight map and the other 46 were tested with this map. The results showed that the ROI segmentation affects the HR measurement significantly and our new algorithm performed better than all the other state of the art methods.

Since the algorithm is based on a supervised model, it is important to generate a reliable SNR weight map. The model may be different in different environmental conditions. It may fail if the registration of the model is imprecise due to the low accuracy of the landmarks detection. The ROI may have different ways to be resized or reshaped to fit the model. For instance, an elliptical model can be generated with the contour provided by the landmarks detection. These offer good opportunities to improve this algorithm.

Figure 1 :

 1 Figure 1: The State of Art ROI segmentation examples: (a) face (b) skin and (c) landmarks.

Figure 2 :

 2 Figure 2: The subjects for the model.

Figure 3 :

 3 Figure 3: Different SNR in different locations.

Figure 4 :

 4 Figure 4: Experimental Setup.

Figure 5 :

 5 Figure 5: Experimental sample images from the UBFC-RPPG database.

Figure 6 :

 6 Figure 6: System Framework.

Table 1 :

 1 The evaluation values for different weight calculation methods.

		Correlation R Precision 2.5 Precision 5 MAE MAE5 RMSE SNR
	a = 10 0.737	0.757	0.874	3.60	1.27	4.84	2.82
	a = 5	0.732	0.755	0.873	3.68	1.27	4.93	3.08
	a = e	0.728	0.755	0.876	3.68	1.29	4.92	3.31
	a = 2	0.695	0.750	0.872	3.99	1.29	5.55	3.35
	a =1.5 0.676	0.740	0.862	4.54	1.27	6.23	3.14
	a = 1	0.674	0.737	0.861	4.56	1.28	6.25	3.09

Table 2 :

 2 The average evaluation values for ROI detection.

		Correlation R	Precision 2.5 Precision 5 MAE MAE5 RMSE SNR
	face	0.531	0.612	0.766	8.72	1.50	12.4	0.570
	skin	0.683	0.716	0.851	5.10	1.36	6.78	3.21
	landmarks 0.649	0.714	0.842	5.24	1.32	7.08	2.83
	model	0.695	0.750	0.872	3.99	1.29	5.55	3.35


Artur Je Ż

Alexander Okhotin
email: alexander.okhotin@utu.fi

COMPLEXITY OF SOLUTIONS OF EQUATIONS OVER SETS OF NATURAL NUMBERS

Systems of equations over sets of natural numbers (or, equivalently, language equations over a one-letter alphabet) of the form Xi = ϕi(X1, . . . , Xn) (1 i n) are considered. Expressions ϕi may contain the operations of union, intersection and pairwise sum A+B = {x+y | x ∈ A, y ∈ B}. A system with an EXPTIME-complete least solution is constructed, and it is established that least solutions of all such systems are in EXPTIME. The general membership problem for these equations is proved to be EXPTIME-complete.

Introduction

The study of expressions over sets of numbers and of the computational complexity of their properties began in the paper by Stockmeyer and Meyer [START_REF] Stockmeyer | Word problems requiring exponential time[END_REF], who considered subsets of N 0 = {0, 1, 2, . . .} as formal languages over a one-letter alphabet. In this case, concatenation of languages turns into a pairwise addition of elements of sets: X+Y = {x+y|x ∈ X, y ∈ Y }. Stockmeyer and Meyer established that the membership problem for expressions with union, intersection and addition is NP-complete. Some extensions of this result were obtained by Yang [START_REF] Yang | Integer circuit evaluation is PSPACE-complete[END_REF], who considered integer circuits (that is, expressions in which subexpressions may be shared) with one more operation of pairwise multiplication, and established similar complexity results. A systematic study of complexity of expressions and circuits with different sets of operations was carried out by McKenzie and Wagner [START_REF] Mckenzie | The complexity of membership problems for circuits over sets of natural numbers[END_REF][START_REF] Mckenzie | The complexity of membership problems for circuits over sets of natural numbers[END_REF].

In this paper we consider equations over sets of natural numbers, which are a more general device than expressions and circuits, and study the computational complexity of their least solutions, as well as of their membership problem. These equations naturally correspond to language equations over a one-letter alphabet. Language equations have recently become an active area of research, see a recent survey by Kunc [START_REF] Kunc | What do we know about language equations?[END_REF]. In particular, unexpected hardness results on language equations have been obtained by Kunc [START_REF] Kunc | The power of commuting with finite sets of words[END_REF] and by Okhotin [START_REF] Okhotin | Decision problems for language equations with Boolean operations[END_REF][START_REF] Okhotin | Unresolved systems of language equations: expressive power and decision problems[END_REF], and this connection gives another motivation for our study. Recent results by Jeż [START_REF] Jeż | Conjunctive grammars can generate non-regular unary languages[END_REF] on the expressive power of conjunctive grammars provide a technical foundation for our results.

We consider equations in the resolved form      X 1 = ϕ 1 (X 1 , . . . , X n ) . . .

X n = ϕ n (X 1 , . . . , X n ) (*)
in which every variable X i assumes value of a set of nonnegative integers. The right-hand side ϕ i of each equation may contain the operations of union, intersection and +, as well as singleton constants. Every such system has a least solution with respect to componentwise inclusion, which can be obtained by fixpoint iteration. Our result, established in Section 3, is a construction of a system (*), such that testing the membership of numbers in its least solution is an EXPTIME-hard problem (with the numbers given in binary notation). The result is obtained by a new kind of arithmetization of an alternating linear-space Turing machine. It is also shown that for every system (*) the membership of numbers in its least solution can be tested in exponential time, which makes the constructed set the hardest.

Let us compare our result to the existing results on expressions and circuits on sets of numbers. Previous research was concerned with the complexity of the general membership problem, where it was sufficient to encode an instance of some hard problem for numbers in an expression or a circuit. In our case, the task is to construct a system that represents a class of problems, while instances of that problem are to be encoded as numbers.

As compared to the research on language equations, our present approach studies a similar problem of constructing a representation of a hard set (cf. Kunc [START_REF] Kunc | The power of commuting with finite sets of words[END_REF], Okhotin [START_REF] Okhotin | Decision problems for language equations with Boolean operations[END_REF][START_REF] Okhotin | The hardest linear conjunctive language[END_REF], Jeż [START_REF] Jeż | Conjunctive grammars can generate non-regular unary languages[END_REF]). However, while encoding a computation of a Turing machine as a string over {a, b} is an ordinary task, in our case we have to encode similar objects as numbers, that is, as strings over a one-letter alphabet. These strings have no apparent structure, and hence the proposed arithmetization is quite unobvious.

This result allows us to establish the complexity of the general membership problem for equations with {∪, ∩, +}, which is stated as follows: "Given a system and a number n 0 in binary notation, determine whether n is in the first component of the least solution of the system". For integer expressions and integer circuits with the operations {∪, ∩, +}, it is known from Stockmeyer and Meyer [START_REF] Stockmeyer | Word problems requiring exponential time[END_REF] and from McKenzie and Wagner [START_REF] Mckenzie | The complexity of membership problems for circuits over sets of natural numbers[END_REF][START_REF] Mckenzie | The complexity of membership problems for circuits over sets of natural numbers[END_REF] that a similar problem is PSPACE-complete. Another weaker model are equations with {∪, +}, that is, without intersection, for which the corresponding problem is NP-complete due to the result of Huynh [START_REF] Huynh | Commutative grammars: the complexity of uniform word problems[END_REF] on the commutative case of the context-free grammars. In our case of equations with {∪, ∩, +}, the general membership problem is EXPTIME-complete, which is established in Section 4. An exponential algorithm for solving this problem is given by a parsing algorithm on conjunctive grammars [START_REF] Okhotin | A recognition and parsing algorithm for arbitrary conjunctive grammars[END_REF].

Language equations and conjunctive grammars

While our results are on the complexity of equations in sets of numbers, our methods are derived from the domain of formal language theory, in particular, from some recent results on language equations.

In language equations, the unknowns are formal languages over an alphabet Σ. If |Σ| = 1, they coincide with equations over sets of numbers, while for larger alphabets they constitute a more general notion. The main object of this study are equations of the resolved form (*), in which variables assume values of sets of non-negative integers, and the right-hand sides may contain the operations of union, intersection and addition of sets. These equations obviously correspond to language equations over a one-letter alphabet with the operations of union, intersection and concatenation, and the recent results on language equations of this kind provide a theoretical foundation, as well as a second motivation, for the present research.

The first type of language equations to be studied were equations of the same form (*) containing union and concatenation, but no intersection: Ginsburg and Rice [START_REF] Ginsburg | Two families of languages related to ALGOL[END_REF] established that these equations provide a natural semantics for the context-free grammars. Equations with added intersection therefore constitute a generalization of the context-free grammars. Definition 2.1 (Okhotin [START_REF] Okhotin | Conjunctive grammars[END_REF]). A conjunctive grammar is a quadruple G = (Σ, N, P, S), in which Σ and N are disjoint finite non-empty sets of terminal and nonterminal symbols respectively; P is a finite set of grammar rules, each of the form

A → α 1 & . . . &α n (where A ∈ N , n 1 and α 1 , . . . , α n ∈ (Σ ∪ N ) * )
while S ∈ N is a nonterminal designated as the start symbol.

The semantics of conjunctive grammars is defined by the least solution of the following system of language equations:

A = A→α 1 &...&αm∈P m i=1 α i (for all A ∈ N ) (2.1)
The component corresponding to each A ∈ N is then denoted by L G (A), and L(G) is defined as L G (S).

The operations used in the right-hand sides of systems (2.1) are union, intersection and concatenation. Since they are monotone and continuous, a least solution always exists and can be obtained by fixpoint iteration as

i 0 ϕ i (∅, . . . , ∅), (2.2)
where ϕ is the right-hand side of (2.1) as a vector operator on |N |-tuples of languages, while denotes pairwise union of vectors of sets. An equivalent definition of conjunctive grammars can be given using term rewriting [START_REF] Okhotin | Conjunctive grammars[END_REF], which generalizes Chomsky's word rewriting. The importance of these grammars lies with the fact that their expressive power is substantially greater than that of the contextfree grammars, while the generated languages can still be parsed in time O(n 3 ), and the practical context-free parsing algorithms, such as recursive descent and generalized LR, admit generalization to conjunctive grammars without an increase in their complexity.

The question of whether conjunctive grammars can generate any non-regular unary language has been an open problem for some years, until recently solved by Jeż [START_REF] Jeż | Conjunctive grammars can generate non-regular unary languages[END_REF], who constructed a grammar for the language {a 4 n | n 0}. Let us reformulate this grammar as the following resolved system of four equations over sets of numbers:

Example 2.2 (Jeż [5]). The system        X 1 = (X 2 + X 2 ) ∩ (X 1 + X 3 ) ∪ {1} X 2 = (X 12 + X 2 ) ∩ (X 1 + X 1 ) ∪ {2} X 3 = (X 12 + X 12 ) ∩ (X 1 + X 2 ) ∪ {3} X 12 = (X 3 + X 3 ) ∩ (X 1 + X 2 )
has least solution

X i = { | base-4 notation of is i0 . . . 0}, for i = 1, 2, 3, 12.
Sets of this kind can be conveniently specified by regular expressions for the corresponding sets of base-k notations of numbers, which in this case are 10 * , 20 * , 30 * and 120 * , respectively. In the following we shall omit some parentheses in the right-hand sides of equations, and assume the following default precedence of operations: addition has the highest precedence, followed by intersection, and then by union with the least precedence.

Using the same technique in a more elaborate construction, a general theorem on the expressive power of unary conjunctive grammars was established. It can be reformulated for equations over sets of numbers as follows:

Theorem 2.3 (Jeż [START_REF] Jeż | Conjunctive grammars can generate non-regular unary languages[END_REF]). For every k 2 and for every finite automaton M over the alphabet {0, . . . , k -1} there exists a system of resolved language equations over N 0 using ∪, ∩, +, such that its least solution is

(S 1 , S 2 , . . . , S n ), where S i ⊆ N 0 and S 1 = { | k-ary notation of is in L(M )}.
Let us note in passing a recent paper by Jeż and Okhotin [START_REF] Jeż | Conjunctive grammars over a unary alphabet: undecidability and unbounded growth[END_REF] establishing a generalization of this result to a larger family of automata recognizing positional notations.

Though representing sets of numbers with a regular positional notation using this type of formal grammars was an unexpected and strong result in terms of language theory, it has no implications on computational complexity, as all these sets are computationally easy. More general representation theorem of Jeż and Okhotin [START_REF] Jeż | Conjunctive grammars over a unary alphabet: undecidability and unbounded growth[END_REF] also does not imply any better complexity results than P-completeness, which, as the present paper shows, is much below the actual complexity of these equations.

Therefore, a new method of constructing such equations is needed to understand their complexity. This step is made in the next section, which introduces an arithmetization technique based upon addition of sets of numbers.

Representing an EXPTIME-complete language

In this section it will be shown that languages defined by least solutions of resolved language equations using +, ∪ and ∩ can be EXPTIME-complete, and this is the hardest language in this family. Denote this family by EQ(∪, ∩, +).

Theorem 3.1. The family EQ(∪, ∩, +) is contained in EXPTIME and contains an EXPTIME-complete language.

The proof is by constructing such a system of equations. The given system encodes a computation of a linear-bounded alternating Turing machine (ATM) It is known that such machines can recognize some EXPTIME-complete languages [START_REF] Chandra | Alternation[END_REF].

In our case we shall consider ATMs operating on a circular tape and moving to the right at every step. Its tape originally contains the input word, and the squares containing it constitute all space available to the machine. Obviously, such machines are as powerful as linear-bounded ATMs of the general form.

Formally, such a machine is defined as M = (Ω, Γ, Q E , Q A , δ, q 0 , q fin ), where Ω is the input alphabet, Γ = {a 0 , a 1 , . . . , a max } ⊃ Ω is the tape alphabet, Q E and Q A are disjoint sets of existential and universal states, respectively, Q = Q E ∪Q A and q 0 , q fin ∈ Q. Given an input w ∈ Ω + , M starts in state q 0 with the head over the first symbol of w. The transition function is δ : Q × Γ → 2 Q×Γ , and the head is moved one symbol to the right at every step. Once the head moves beyond the right-most symbol, it is moved back over the first symbol of w, maintaining its current state; this implements a circular tape. For technical reasons, assume that (q, a ) / ∈ δ(q, a) for all q ∈ Q and a, a ∈ Σ, (that is, the machine never stays in the same state), and that δ(q, a) = ∅ for all q ∈ Q A and a ∈ Σ.

Our construction of a system of equations over sets of numbers simulating a computation is based upon representing instantaneous descriptions of the ATM as numbers. We shall think of these numbers as written in base-(8+|Q|+max(|Q|+7, |Γ|)) positional notation, and the entire argument is based upon mapping the symbols used by the machine to digits, and then using addition to manipulate individual digits in the positional notation of numbers. It must be noted that this positional notation is only a tool for our understanding of the constructions, while the actual equations deals with numbers as they are.

Let Σ = {0, 1, . . . , 7 + |Q| + max(|Q| + 7, |Γ|)} be the alphabet of digits, and define the mapping of symbols to digits, • : Q ∪ Γ → Σ, as follows:

q i = 7 + i (for q i ∈ Q) a i = 7 + |Q| + i (for a i ∈ Γ)
Furthermore, let Q = { q | q ∈ Q} and Γ = { a | a ∈ Γ}. Now the tape of the ATM containing symbols a i 1 . . . a in , with the head over the j-th symbol and the machine in state q, is represented as the following string of digits: 0 a i 1 . . . 0 a i j-1 q a i j 0 a i j+1 . . . 0 a in 0 ∈ Σ * For technical reasons, configurations in which the head has just moved over the last symbol but has not yet jumped to the first position are considered separately, and will be represented as strings of the form 0 a i 1 . . . 0 a in q , where q is the current state. Note that digits denoting letters are written only in even positions, while odd positions are reserved for the states of the Turing machine. The set of all strings of digits representing valid encodings of tapes is specified by the following regular expression over Σ: Tape = (0 Γ ) * Q ( Γ 0) * \ Q The set Tape should be a considered as a formal language over Σ, which will be used later as a part of representations of some sets of numbers. Subsets of this set representing tapes with different states will be denoted as follows:

Tape u = {w | w ∈ Tape, u is a substring of w} Tape u = {w | w ∈ Tape, u is a prefix of w}
Besides the contents of the tape, the encoding for Turing machine configurations uses a counter of rotations of the circular tape. This counter specifies the number of passes through the tape the machine is still allowed to make before it must halt. It is represented in binary notation using digits {0, 1}, and the set of valid counter representations is Counter = 1{0, 1} * Normally the counter uses only digits {0, 1}, but in order to implement the incrementation of the counter we shall use strings with one digit 2 representing zero with carry. The set of valid representations of counters with a carry is

Counter = 1{0, 1} * 2{0, 1} * ∪ 2{0, 1} *
For every string c k-1 . . . c 0 ∈ Counter ∪ Counter , define its value as

Value(c k-1 . . . c 0 ) = k-1 j=0 c j • 2 j .
Now define the mapping from configurations of the Turing machine to numbers. A configuration with the tape contents, head position and current state given by a string of digits w ∈ Tape, and with the counter value given by x ∈ Counter is represented by a string of digits x55w, where two marker digits 55 separate the values. This string of digits in base-|Σ| positional notation specifies a certain number, which accordingly represents the configuration.

The key property of this encoding is that every transition of the ATM reduces the numerical value of its configuration. Indeed, if the head is moved to the right, then a digit q is replaced with 0 and all other modifications are done on less signigicant digits. If the head jumps from the end to the beginning, then the counter is decremented, and since the counter occupies more significant positions in the number than the tape, this transition decreases the value of the configuration as well. This monotonicity allows us to encode dependence of configurations on each other by using addition of nonnegative numbers only.

The construction of equations representing the computation of the ATM begins with some expressions that will be used in the right-hand sides of equations. These expressions contain some constant sets of numbers given as regular languages over the alphabet Σ. Every such language represents the set of all numbers with |Σ|-ary notation of the given form. According to Theorem 2.3, every such set can be represented by a separate system of equations using only singleton constants. All these subsystems are assumed to be included in the constructed system, and each of the regular expressions in the system can be formally regarded as a reference to one of the auxiliary variables.

Definitions of a few of these regular languages incorporate positional notations of numbers obtained by subtracting one number from another. For convenience, these values are given in the form u v, with u, v ∈ Σ * being positional notations of two numbers (the former shall be greater or equal to the latter). One can write, e.g., (u v)0 * for the set of all numbers with their |Σ|-ary notation beginning with fixed digits determined by the given difference, followed with any number of zeroes.

Step(X) = q∈Q E a∈Γ (q ,a )∈δ(q,a)

Move q ,a ,q,a (X) ∪ q∈Q A a∈Γ (q ,a )∈δ(q,a)

Move q ,a ,q,a (X)

Move q,a,q ,a (X) = X ∩ Counter 55 Tape a q + q a 0 a q (00) * ∩ Counter 55 Tape q a Jump(X) = q X ∩ Counter 55 Tape q + (1000 q )(00 The construction uses two variables, X and Y . Either variable represents the set of proper configurations of the machine, starting from which the machine accepts. The variable X represents configurations belonging to the set Counter 55 Tape, while Y represents configurations from (Counter ∪ Counter )55 Tape, in which the counter may contain one carry digit 2 that needs to be propagated to higher positions. The equations, using the above auxiliary functions, are as follows:

) + + q ∩ (Counter ∪ Counter )55 Tape q Carry(Y ) = Y ∩ {0,
X = Final ∪ Step(X) ∪ Y ∩ Counter 55 Tape (3.1) Y = Jump(X) ∪ Carry(Y ) (3.2)
In order to determine the least solution of this system, let us first establish some properties of the auxiliary functions.

The first quite elementary property is their distributivity over infinite union, which allows us to study these operations as operations on individual numbers, and then infer their action on sets of numbers.

Lemma 1 (Distributivity). Each function f ∈ {Move q,a,q ,a , Jump, Carry} is distributive over infinite union, in the sense that f (S) = n∈S f ({n}) for every S ⊆ N 0 . This follows from the fact that each of these expressions consists of intersections with constant sets, sums with constant sets and unions. On the other hand, note that if an expression contains intersections or sums of multiple expressions involving X, then it is not necessarily distributive over infinite union; in particular, Step need not be distributive.

One of the main technical devices used in these functions is addition of a constant set of numbers with |Σ|-ary notation u0 * (that is, a set m • |Σ| i i 0 ) with one, two or three non-zero digits in u. The following lemma establishes that this addition can never rewrite the double markers 55, that is, every sum in which these markers are altered does not represent a valid tape contents. This means that every such addition manipulates the counter and the tape separately, and the changes do not mix.

Lemma 2 (Marker preservation). For every x, x ∈ {0, 1, 2, 3} * \ 0Σ * and w, w ∈ Tape, if

x 55w ∈ x55w + (Σ 3 ∪ Σ 2 ∪ Σ)0 * , then |w| = |w |.
The next statement describes the operation of Carry: applied to a configuration with the counter having a single carry digit 2, Carry changes this digit to 0 and increments the next digit, making it 1 or 2. Note that all operations are in |Σ|-ary notation. The tape contents is not altered. The next lemma states the functionality of Jump, which can be described as follows. If Jump is applied to a configuration in which the head scans over the first symbol, then the result of the operation is the previous configuration, in which the head is at the right-most position beyond the end of the string, while the value of the counter x is greater by 1.

Lemma 5. Let x = xc ∈ Counter with c ∈ {0, 1} and w = q w0 ∈ Tape with q ∈ Q, that is, w encodes a configuration with the head over the first symbol. Then Jump(x55w) = { x(c + 1)550 w q }.

For any string α ∈ Σ * of a different form, Jump(α) = ∅.

It follows from Lemma 5 that Jump is a reversible function, that is, the previous configuration given by Jump(x55w) corresponds to x55w only. This is stated as follows: Lemma 6. Let x 55w ∈ Jump(x55w). Then w = 0 w q and w = q w0 for some state q, and Value(x ) = Value(x) + 1.

Let us now proceed with specifying the action of Move, which represents symbol manipulation, head movement and state change of a Turing machine according to the membership of states and symbols specified in δ. Generally, when Move q,a,q ,a is applied to a valid configuration, it computes the preceding configuration of the machine. This configuration is unique because of the restriction built in Move q,a,q ,a in its subscripts. The symbols and states used as the subscript restrict its applicability to the following case: in the current configuration the machine is in state q and the symbol to the left rewritten at the previous step is a, while in the previous configuration the machine was in state q and scanned the symbol a . For all other configurations and in all other cases, the function produces the empty set. Lemma 7. Let q, q ∈ Q and a, a ∈ Γ. Let x ∈ Counter and w = w0 a q w ∈ Tape for some w ∈ (0 Γ ) * and w ∈ ( Γ 0) * . Then Move q,a,q ,a (x55w) = x55 w q a 0 w.

For every string α ∈ Σ * of a different form, Move q,a,q ,a (α) = ∅.

Similarly to Lemma 6, reversibility of Move q,a,q ,a directly follows from Lemma 7.

Lemma 8. Let x55w ∈ Move q ,a ,q,a (x 55w ). Then w = w q a 0 w and w = w0 a q w for some w ∈ (0 Γ ) * and w ∈ ( Γ 0) * , and x = x .

Having established a solution complexity theorem for equations over sets of numbers, let us discuss its implications on conjunctive grammars over a one-letter alphabet.

Every conjunctive language is in P [START_REF] Okhotin | Conjunctive grammars[END_REF], and some conjunctive languages over a multipleletter alphabet are known to be P-complete [START_REF] Okhotin | The hardest linear conjunctive language[END_REF]. The case of a unary alphabet is special, as it is known that no sparse language, in particular no unary language, can be P-complete unless DLOGSPACE = P [START_REF] Ogihara | Sparse hard sets for P yield space-efficient algorithms[END_REF][START_REF] Cai | Sparse hard sets for P: resolution of a conjecture of Hartmanis[END_REF], that is, unless the notion of P-completeness is trivial. However, from Theorem 3.1 one can infer the following result slightly weaker than P-completeness: Corollary 3.2. There exists a EXPTIME-complete set of numbers S ⊆ N, such that the language L = {a n |n ∈ S} of unary notations of numbers from S is generated by a conjunctive grammar.

Note that for every unary language generated by a conjunctive grammar, the corresponding set of numbers is in EXPTIME. The set constructed in Corollary 3.2 can thus be regarded as the computationally hardest among unary conjunctive languages.

A simple consequence of Corollary 3.2 refers to the complexity of parsing for conjunctive grammars.

Corollary 3.3. Unless PSPACE = EXPTIME, there is no logarithmic-space parsing algorithm for conjunctive languages over a unary alphabet.

The membership problem

Consider the general membership problem for our equations, stated as follows: "Given a system X i = ϕ i (X 1 , . . . , X m ) and a number n in binary notation, determine whether n is in the first component of the least solution of the given system". Its complexity is now easy to establish. Theorem 4.1. The membership problem for resolved systems of equations over sets of numbers with operations {∪, ∩, +} is EXPTIME-complete.

Proof. Membership in EXPTIME. The algorithm begins with representing the given system as a conjunctive grammar over a unary alphabet, with a linearly bounded blow-up. The given number n is represented as a string a n with an exponential blow-up. Then it is sufficient to apply the known polynomial-time algorithm for solving the membership problem for conjunctive grammars [START_REF] Okhotin | A recognition and parsing algorithm for arbitrary conjunctive grammars[END_REF].

The EXPTIME-hardness of the general membership problem immediately follows from Theorem 3.1 by fixing the system of equations.

Let us conclude by comparing the complexity of the membership problem for expressions, circuits and equations, as well as the families of sets representable by their solutions. All known results are given in Table 1.

The new complexity results for the equations over sets of numbers naturally fit into the framework of the existing research. On the other hand, the new results on the expressive power of equations come in a sharp contrast with the previous work: these equations can represent non-trivial sets of numbers, which are computationally as hard as the general membership problem for this class.

It remains an open question, what is the exact family of sets of natural numbers defined by these equations. For instance, is it possible to represent the set of all primes?

  1} * 2{0, 1} * 55 Tape + 10 * ∩ {0, 1} * 3{0, 1} * 55 Tape + 10 3 0 * ∩ {0, 1} + ∪ {0, 1} * 2{0, 1} * 55 Tape In addition, define the set of final configurations of the machine: Final = Counter 55 Tape q fin

Lemma 3 (

 3 Carry propagation). For every x ∈ Counter and for every w ∈ Tape, Carry {x55w} = {x 55w}, where x ∈ Counter ∪ Counter and Value(x ) = Value(x). If x ∈ Counter , then the position of 2 in x is greater than the position of 2 in x. According to Lemma 3, Carry moves the carry by one position higher. The next lemma shows that sufficiently many iterations of Carry always eliminate the carry digit: given a counter with the notation x = x01 k-1 2, Carry k transforms it to x = x10 k-1 0. Lemma 4 (Termination of carry propagation). For every x ∈ Counter ∪ Counter and w ∈ Tape there exists x ∈ Counter and k 0, such that Carry k (x55w) = x 55w and Value(x) = Value(x ).

Table 1 :

 1 Comparison of formalisms over sets of integers.

		Representable sets	Membership problem
	expressions with {∪, +}	Finite	NP-complete [17]
	circuits with {∪, +}	Finite	NP-complete [4, 9, 10]
	equations with {∪, +}	Ultimately periodic	NP-complete [4]
	expressions with {∪, ∩, +}	Finite	PSPACE-complete [17]
	circuits with {∪, ∩, +}	Finite	PSPACE-complete [9, 10]
	equations with {∪, ∩, +}	EXPTIME, contains EXPTIME-complete set	EXPTIME-complete

This work is licensed under the Creative Commons Attribution-NoDerivs License. To view a copy of this license, visit http://creativecommons.org/licenses/by-nd/3.0/.

(A. Jeż) Supported by MNiSW grant number N206 024 31/3826, 2006-2008, and by a short visit grant from the European Science Foundation under project AutoMathA, reference number 1763. (A. Okhotin) Supported by the Academy of Finland under grant 118540. c A. Jeż and A. Okhotin

The flow control of the alternating Turing machine includes existential and universal nondeterminism in the corresponding states, and a single step is in fact a disjunction or conjunction of several transitions as specified in Move. This logic is transcribed in the expression Step(X), which computes the set of all previous configurations, from which machines in a universal state make all their transitions to configurations in X and machines in an existential state make at least one of their transitions to some configuration in X. This implements one step of the computation of the machine, backwards. Lemma 9. Let x ∈ Counter and w ∈ Tape, let q ∈ Q be the state encoded in w. Then x55w ∈ Step(X) if and only if

• the configuration w has the head not in the position beyond the right-most symbol, that is, w = w q w0 for some w, w ∈ Σ * . • if q ∈ Q E , then for some string w encoding next configuration of the ATM there holds x55w ∈ X. • if q ∈ Q A , then for every string w encoding next configuration of the ATM there holds x55w ∈ X.

Having established the formal meaning of the auxiliary operations, let us return to the equations. The equation for X states that a configuration leads to acceptance if and only if it is accepting itself (Final), or one can directly proceed from it to a configuration leading to acceptance (Step(X)), or that it is a configuration obtained in Y . The equation for Y specifies circular rotation of the tape by Jump(X) and implements iterated carry propagation as in Lemma 4 by a self-reference Carry(Y ). Altogether, the least solution of these equations corresponds to the computation of the machine as follows: It remains to observe that the number of steps of the machine is exponentially bounded, hence the acceptance of a word by the machine is represented by the following number in the least solution of the constructed system: Main Lemma. ATM M accepts a string a 1 . . . a n ∈ Ω + if and only if

Proof of Theorem 3.1. The system of equations constructed above has an EXPTIMEcomplete least solution.

To see that the least solution of every system is in EXPTIME, it is sufficient to represent it as a conjunctive grammar over a unary alphabet. Then, given a number n, its membership in the least solution can be tested by supplying the string a n to a known cubic-time parsing algorithm for conjunctive grammars [START_REF] Okhotin | Conjunctive grammars[END_REF]. Its time is cubic in n, hence exponential in the length of the binary notation of n.