Franco Fusco

Olivier Kermorgant

Philippe Martinet

A Comparison of Visual Servoing from Features Velocity and Acceleration Interaction Models

Visual Servoing has been widely investigated in the last decades as it provides a powerful strategy for robot control. Thanks to the direct feed-back from a set of sensors, it allows to reduce the impact of some modeling errors and to perform tasks even in uncertain environments. The commonly exploited approach in this field is to use a model that expresses the rate of change of a set of features as a function of sensor twist. These schemes are commonly used to obtain a velocity command, which needs to be tracked by a low-level controller.

Another approach that can be exploited consists in going one step further and to consider an acceleration model for the features. This strategy allows also to obtain a natural and direct link with the dynamic model of the controlled system. This study aims at comparing the use of velocity and acceleration-based models in feed-back linearization for Visual Servoing. We consider the case of a redundant manipulator and discuss what this implies for both control techniques. By means of simulations, we show that controllers based on features acceleration give better results than those based on velocity in presence of noisy feedback signals.

I. INTRODUCTION

Visual Servoing is a control technique that allows precise relative positioning between a sensor and an observed object. It is robust to modeling errors and features a large convergence domain. The main research interests have been on modeling the interaction between descriptive features and sensor motion at the velocity level. Controllers are usually designed in two decoupled stages: a first high-level layer operates in the features space and produces a reference velocity [START_REF] Chaumette | Visual servo control. I. Basic approaches[END_REF], [START_REF]Visual servo control. II. Advanced approaches[END_REF] which is then forwarded to a low-level robot controller actuating the motors to track the desired velocity profile [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF], [START_REF] Fang | Adaptive 2.5D visual servoing of kinematically redundant robot manipulators[END_REF], [START_REF] Wang | Visual servoing considering sensing dynamics and robot dynamics[END_REF].

In recent years, the idea of directly linking the dynamics of robotics systems to sensory feedback has been considered. This had already been proposed over than two decades ago [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF], [START_REF] Zhang | Visual servoing with dynamics: control of an unmanned blimp[END_REF], with the appealing option of using a second-order model expressing the link between features acceleration and robot kinematics, altogether with the Dynamic model of the robot. This allows to define, e.g., Computed Torque Controllers directly in the feature space. These models have now been "rediscovered" and used in different contexts, such as vision-based aerial grasping [START_REF] Thomas | Toward image based visual servoing for aerial grasping and perching[END_REF], stero-vision systems [START_REF] Mohebbi | An acceleration command approach to robotic stereo image-based visual servoing[END_REF] and hybrid vision-force control for serial manipulators [START_REF] Vandernotte | Dynamic sensor-based control[END_REF].

In this work, we contribute to the analysis of second-order Visual Servoing control schemes by performing a systematic comparison with classical first-order controllers. In the recent literature, second order models have been mainly derived ad 1 Centrale Nantes, Laboratoire des Sciences du Numérique de Nantes LS2N, France. Mail: name.surname@ls2n.fr 2 Inria Sophia Antipolis, France. Mail: philippe.martinet@inria.fr hoc from the definition of considered features. We instead consider a more general approach, which does not rely on the selected set of descriptors. While this concept had been already described briefly in [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF], we contribute with more explicit results in the development.

In addition, we consider the case of a redundant manipulator featuring 7 degrees of freedom. Our objective in this sense is to focus the attention on the fact that second-order schemes must explicitly handle the redundancy even when no additional tasks are specified, whereas first-order controllers naturally cope with this added complexity.

The remainder of this paper is organized as follows: in Section II we recall basic modeling exploited in classical Visual Servoing and illustrate how acceleration models can be obtained. We then review how to link robot kinematics and dynamics to the feature space, and how to use these models to design controllers. We also discuss the implications of redundancy in the proposed control laws. In Section III we compare in a simulated environment control laws based on the first and second-order models. Reported results show that models based on features acceleration outperform those exploiting the cascade of Visual and low-level controllers. In particular, they exhibit good convergence properties and smoother control inputs, with less sensitivity to noise.

II. VISUAL SERVOING & DYNAMICS

A. First order Visual Servoing

We recall here the classical first-order modeling approach generally exploited in the context of Visual Servoing. To achieve the task of regulating the relative position between a sensor and an observed object, a set of descriptive features is firstly selected. Their numerical value is stored in the m-dimensional vector s. This set might include coordinates extracted from an image, the relative pose between the sensor frame and the object frame, distance measurements, etc. We will focus in the following on the case of Image Based Visual Servoing (IBVS), although the models presented below are generic and can be applied to a broader class of features and sensors.

Once the features have been chosen, it is necessary to relate them to the kinematics of the robotic system. This is achieved via the following differential relation:

ṡ = L s v + ∂s ∂t (1)
L s being the Interaction Matrix of the set of features and v T = v T ω T representing the kinematic screw, or twist, of the camera with respect to the object (expressed in the coordinate system of the sensor itself). The last term in (1) takes into account variations of the features that do not depend on the relative motion between the sensor and the object. It is generally possible to assume that additive component to be null. It is worth remarking that L s generally depends on both the observed features and a set of parameters. Among these, some can be assumed constant, such as intrinsic camera parameters. On the other side, there appear also quantities in L s whose value changes depending on the relative motion between the camera and the object. Examples of such parameters are the depth of image points, or plane parameters when dealing with planar image moments. In the following, the vector containing such parameters is denoted as z.

Starting from first-order models, it is shown in the next section how to obtain a relation that involves the second derivative of the features and the acceleration of the sensor.

B. Second order Visual Servoing

This class of models considers the second-order time derivative of the features and how it depends on the acceleration of sensors, in addition to their velocity. In related works the derivation of the second order model was achieved by considering the direct expression of the features and by evaluating explicitly s depending on its definition. We propose in the following a more general approach that we believe facilitates the derivation of these models. First of all, the derivative of (1) gives (assuming ∂s ∂t = 0):

s = L s d dt v + Ls v (2)
The derivative of the kinematic screw writes as follows:

d dt v = a + v × ω 0 (3)
which is obtained by considering the camera frame is not inertial with respect to the object. Regarding the differentiation of L s , it is convenient to consider this matrix as written in the form:

L s = 1 • • • m T (4)
i being the rows of L s , from which the derivative can be expressed as:

Ls =    ∂ 1 ∂s ṡ + ∂ 1 ∂z ż T • • • ∂ m ∂s ṡ + ∂ m ∂z ż T    (5)
For simplicity we assume in the following that an expression for ż is available, and that furthermore it takes the form ż = L z v. Under this assumption and by introducing

H i = ∂ i ∂s L s + ∂ i ∂z L z T : Ls =   v T H 1 • • • v T H m   (6)
As a result, (2) becomes:

s = L s a + L s v × ω 0 +   v T H 1 • • • v T H m   v (7)
The second term in the right-hand side can be rearranged in a nicer form that matches the structure of the last term.

For this purpose, one can consider the first three elements in each row of L s and place them in a set of vectors

(v)
i . After some simple manipulation, it is possible to rewrite [START_REF]Visual servo control. II. Advanced approaches[END_REF] in the form:

s = L s a +   v T (H 1 + L 1 ) v • • • v T (H m + L m ) v   = L s a + h s (8)
wherein matrices L i take the form:

L i = O 3 O 3 (v) i × O 3 (9)
having used the symbols O 3 and [•] × to denote respectively the null 3-by-3 matrix and the skew-symmetric matrix corresponding to the cross-product.

The main interest in the given development is that to get the second order model of any feature, the only required work is to evaluate ∂ i ∂s and ∂ i ∂z . We conclude this section by showing how a second order model is obtained from another set of features for which an acceleration model has already been obtained. Consider a set of features s a defined by some differentiable function of another set of features s b , i.e., s a = f (s b ). After development, one would get the following result:

sa = ∂f ∂s b L b a + ∂f ∂s b h b + d dt ∂f ∂s b L b v (10)
in which it is easy to identify the Interaction Matrix, multiplying the acceleration screw a, and the vector h a .

C. Robot Kinematics and Dynamics in Sensor Space

If we assume that the relative motion between the sensor and the object depends solely on the motion of the robot, the models presented in the previous section can be easily linked to the configuration space of the robotic system. Assuming an eye-in-hand setup, the velocity and acceleration of the sensor are given by:

v = TJ q (11a) a = TJq + T J q (11b)
wherein J is the 6-by-n Jacobian matrix of the robot, n being the number of degrees of freedom of the robot. The two vectors q, q ∈ R n collect respectively positions and velocities of the joints, while T is a twist-rotation matrix that projects the kinematic screw on the sensor frame. By injecting (11a) and (11b) into the first and second order models respectively, the following kinematic relations are obtained 1 :

ṡ = L s TJ q = J s q (12a) s = J s q + L s T J q + h s = J s q + h q (12b)
where the matrix J s = L s TJ is known as feature Jacobian.

We also recall here that the Inverse Dynamic Model (IDM) of a robotic system composed of a set of rigid bodies writes as:

τ = Mq + b (13)
wherein τ represents joint efforts while M is the generalized inertia matrix of the system, grouping inertial effects due to masses and inertiae of rigid bodies and motors. Finally, b groups Coriolis and centripetal forces, the efforts exerted on the system by gravity and the forces and moments acting on the actuators due to friction. Assuming that M is invertible, which is always true for a serial kinematic chain, the joint accelerations are given by q = M -1 (τ -b), thus leading to:

J s M -1 τ = s -h q + J s M -1 b (14)
This last equation provides a direct link between the commanded joint efforts and the acceleration of features.

D. Visual Servoing

Exploiting the illustrated kinematic and dynamic models, it is possible to build a variety of controllers that steer the observed features to a desired value s .

The objective of the control is to regulate the feature error e s . = s-s to zero. A proportional control with feed-forward can be applied when considering the first order model as:

q 1 = J + s (ṡ -λe s ) (15)
whereas a proportional and derivative one can be exploited if the second order model is used:

q 2 = J + s s -k d ės -k p e s -h q (16)
In both cases, J + s represents an estimation of the Moore-Penrose pseudo-inverse of J s . Using ( 16) also requires an estimation of h q and ṡ. Finally, the desired set of features is often constant, i.e., ṡ = 0 and s = 0, and therefore the controllers are slightly simplified.

The two reference commands correspond to the minimum norm-solutions of the respective minimization problems:

arg min q J s q -ṡ + λe s (17a) arg min q J s q -s + k d ės + k p e s + h q (17b)
They aim at enforcing an exponential decrease of the error according to the evolution of the autonomous linear systems ės + λe s = 0 and ës + k d ės + k p e s = 0, both asymptotically stable if positive gains are selected. Perfect linearization will not be achievable in general, but asymptotic stability is ensured for (15) provided that J s J + s > 0 [START_REF] Chaumette | Visual servo control. I. Basic approaches[END_REF]. Additional conditions need to be satisfied for the stability of ( 16) [START_REF] Keshmiri | Augmented image-based visual servoing of a manipulator using acceleration command[END_REF].

The output obtained from these kinematic control laws can be finally used in conjunction with the IDM to evaluate a set of effort commands to be sent to the actuators. For the first kind of regulation, the provided reference is a joint velocity command. In order to exploit the dynamics, a possible choice is to add another control layer providing the auxiliary acceleration signal

w 1 = qff + γ ( q 1 -q) (18)
given a positive proportional gain γ and the feed-forward term qff corresponding to the numerical derivative of [START_REF] Spica | Plane estimation by active vision from point features and image moments[END_REF]. This auxiliary command is finally inserted in the IDM, thus giving:

τ = Mw 1 + b (19)
Regarding ( 16), the produced output can be used directly to compute the efforts from the IDM.

E. Dealing with Redundancy

We address separately here the case in which the robot has more degrees of freedom than those required to perform the servoing task. Assuming that no singularities are encountered, we can state that the robot is redundant if n > rank (J s ). In this case, as the kernel of the matrix is not null, infinite solutions can be found that are compatible with (17a) and (17b). As outlined before, ( 15) and ( 16) give minimum norm solutions for the joint velocity and acceleration respectively.

Other controls which are compatible with the second-order Visual Servoing task can be found by using directly [START_REF] De Luca | On-line estimation of feature depth for image-based visual servoing schemes[END_REF] to compute the efforts from the desired task [START_REF] Vandernotte | Dynamic sensor-based control[END_REF], [START_REF] Thomas | Toward image based visual servoing for aerial grasping and perching[END_REF]. Two possible choices in this sense are to minimize τ or a weighted norm of the efforts τ T Dτ (with D > 0):

τ = J s M -1 + s -k d ės -k p e s -h q + J s M -1 b (20a) τ = (J s ) + D s -k d ės -k p e s -h q + J s M -1 b (20b) wherein (J s )
+ D represents a weighted pseudo-inverse. Regarding these two possibilities, it must be noted that the latter should generally be preferred, as the former might lead to high uncontrolled velocities [START_REF] De Luca | Efficient Dynamic Resolution of Robot Redundancy[END_REF].

Even though in this work we are interested only in a single task, it must be noted that, when dealing with a redundant manipulator, schemes based on the second-order model will not be enough to stabilize the robot to a given configuration. A classical first-order control, without supplementary tasks, will return the minimum norm joint velocity required to reach an equilibrium configuration. Regardless this being a global or local minimum, after the convergence the computed command will be zero, and therefore the robot will not move anymore (we do not consider perturbations here). The same is not true for the other class of controllers. Indeed, at equilibrium the command corresponds to a null acceleration, which is compatible with a non-null instantaneous velocity lying in the kernel of J s . It is thus necessary to ensure the given control law damps the velocity of the robot.

A possible solution is to add the following secondary task:

q = -k v q (21)
for a positive gain k v . The final control law would thus become:

q 2 = q 2 -ZZ T (k v q + q 2 ) (22)
with q 2 given as before by [START_REF] Mansard | Task sequencing for high-level sensorbased control[END_REF]. Z is an orthogonal basis of the kernel of J s , obtainable via, e.g., Complete Orthogonal Decomposition, so that the orthogonal projector is given in the efficient form ZZ T [START_REF] Escande | Hierarchical quadratic programming: Fast online humanoid-robot motion generation[END_REF]. The efforts can then be obtained by injecting into (13) q 2 instead of q 2 . As a concluding remark, it might be considered to inject the velocity damping task into (20a) or (20b), i.e., solving the problem directly at the effort level. However, the Servoing task altogether with the velocity damping objective are fully constraining the degrees of freedom of the robot, and, as direct consequence, the final result would not change. On the contrary, as these models require more operations, e.g., the inversion of M, in all our experiments they are not considered.

III. SIMULATION RESULTS

A comparison between classical first-order control schemes and second order ones has been performed in a simulated environment. We report in the following the details of our implementation and the results of this study.

Tests have been realized using Gazebo for physics simulation and ROS for the implementation of controllers and communication. ViSP [START_REF] Marchand | ViSP for Visual Servoing[END_REF] has been used for Visual Servoing related computations.

Our setup involves a positioning task from four coplanar points located at the vertices of a square. We consider a perspective camera mounted on the tip of a Kuka LWR4+ robot, a serial manipulator featuring 7 degrees of freedom. We assume a perfectly calibrated sensor with no distortion, and use directly the normalized coordinates of the observed points. We also do not consider uncertainties in the Geometric and Dynamic parameters of the robot.

The control period in simulation is set to 1 ms. We consider that the visual feed-back comes from a mixture of image processing and state estimation algorithm, and that the unknown parameters z can be estimated as well [START_REF] De Luca | On-line estimation of feature depth for image-based visual servoing schemes[END_REF], [START_REF] Spica | Plane estimation by active vision from point features and image moments[END_REF]. To take into account uncertainties due to these operations, we inject normally distributed noise in both s and z. In addition, the presence of encoders on each motor is simulated by quantization of the readings. More details on noise are presented in the next section.

Three control strategies are analyzed. The first one, namely First-Order (FO), corresponds to a velocity-based scheme without the feed-forward signal qff . The second strategy takes this term into account and will be referred to as First-Order with Feed-Forward (FOFF). Finally, a Second-Order (SO) scheme that evaluates joints acceleration via ( 16) is considered. Implementation and tuning of these controllers are discussed in Section III-B

A. Noise

Noise is introduced in the simulation to obtain a more realistic study-case. The first source of noise we consider is added to the current value of the features s. In a real setup this might come from, e.g., the image processing routine that detects the points of interest or estimators that produce a continuous feedback when an image is not available. We therefore add to each feature coordinate a random value, drawn from a truncated Gaussian distribution with zero mean. The truncation is performed symmetrically at ±3σ, σ being the standard deviation of the unbounded distribution.

In the experiments, we used σ = 2 • 10 -3 (corresponding to about 1 pixel). A similar strategy has been used to simulate uncertainties in the unknown depths of the observed points, by adding to z a noise with σ = 3 • 10 -3 .

Regarding the measurement of joint positions and velocity, we assume that only encoders are available. These devices are simulated by discretizing the values of joints angles with a given resolution ∆q. For this purpose we used ∆q = 10 µrad.

Noisy signals are then filtered using low-pass Butterworth 4 th -order filters and finally used in all the calculations performed by controllers.

B. Controllers Implementation and Parameters Tuning

1) First-Order (FO) Controller: One problem with the first order law [START_REF] Spica | Plane estimation by active vision from point features and image moments[END_REF] is that it presents a discontinuity in velocity at the very first iteration. In order to guarantee a "smooth start", we exploit the concepts presented in [START_REF] Mansard | Task sequencing for high-level sensorbased control[END_REF]. They showed that when a new task is activated, it is possible to achieve a smooth transition in the control input by adding a homogeneous term to the command. Applied to our case, the velocity reference q 1 is now evaluated as:

q 1 = J + s ṡ -λe s + e -µt e s (0) (23)
where e s (0) is the features error for t = 0, and we assume null desired features velocity at the beginning (ṡ (0) = 0). Due to this modification, the command at the first iteration is zero, thus enforcing a smooth transition of the desired velocity of the manipulator. The effect on the error dynamics is that it now behaves as an equivalent second order system -if perfect linearization was possible -according to the autonomous equation ës + (µ + λ) ės + µλe s = 0 (24)

One interesting consequence of this choice is that it is possible to enforce the same desired behavior of the features evolution with both a first and a second order control. Indeed, in a perfect world not affected by noise, the two control strategies are perfectly equivalent in terms of error evolution.

After computing the desired velocity (23), the effort command is computed as τ = γM ( q 1 -q) + b.

Regarding gains tuning, we decide to impose a critically damped behavior for the transient response, as this gives the fastest convergence without oscillations:

λ = µ = ω n (25)
wherein ω n > 0 defines the roots of the characteristic polynomial of the second order system. The autonomous response of the error would be given by e s (t) = e -ωnt (1 + ω n t) e s (0), which is strictly decreasing and converges to zero. We finally select the value of ω n so that after a time T the error is a fraction of its original value, i.e., e -ωnT (1 + ω n T ) = α, α ∈ (0, 1). Solving the equation numerically for α = 0.05 and T = 2 s gives ω n 2.37 s -1 . Concerning the proportional velocity gain γ in (18), the higher its value the better will be the tracking of the desired velocity profile. On the other hand, it will increase the control efforts τ , and therefore a trade-off has to be found. In our implementation, a value of γ = 8 s -1 has been chosen, which allows a velocity error to decrease of 80 % in 0.2 s.

2) First-Order with Feed-Forward (FOFF): The implementation of this controller is very similar to the one of FO. We use here as well the smooth velocity transition and the same set of gains. The feed-forward term qff is computed as the numerical derivative of the reference command q 1 . It is well-known that this highly amplifies noise, resulting in a generally unfeasible, if not just undesirable, control signal. To get a nicer behavior and reduce the impact of noise, (23) is thus filtered prior to differentiation. This considerably reduces the shakiness of the controller. The final efforts sent to the motors are evaluated using (18) and the IDM.

3) Second-Order (SO) Controller: This controller exploits the acceleration in the feature space. It computes joints accelerations using (22) to complete the visual task while damping the velocity of the manipulator, and subsequently exploits [START_REF] Marchand | ViSP for Visual Servoing[END_REF] to compute the actuation signal. To obtain an estimation of ės , (11a) was used. The control gains k p and k d are selected to obtain a critically damped evolution as in the FO case, with k p = ω 2 n , k d = 2ω n , and ω n = 2.37 s -1 . The last parameter to be tuned is k v , the gain appearing in the velocity damping task (21). As we assume the Servoing Task to be completed in about 2 s, this parameter is tuned so that the convergence of velocity to zero is achieved after a longer time. In practice, we used the value k v = 1 s -1 , which theoretically ensures that any non-null velocity will be reduced to about the 5 % of its initial value in 3 s.

C. Regulation to a Fixed Configuration

The first set of tests involves positioning the end-effector in front of a set of four points at a distance of 0.5 m, see Fig. 1. As the desired value of the features does not change during the tests, in the controllers the quantities ṡ and s are set to zero. The initial relative transformation from the camera to the object presents a rotation of about 90 • around the optical axis of the sensor, leading to the well-known retraction problem. Since this affects all control laws, we did not try to solve the problem to optimize the 3D behavior of the camera. We believe that a better choice of the descriptive features, such as almost decoupled image moments [START_REF] Chaumette | Image Moments: A General and Useful Set of Features for Visual Servoing[END_REF], would result in a trajectory of higher quality in all the cases.

We report results related to two cases, an ideal test run in absence of noise (Fig. 2) and a more realistic case with disturbances as discussed in Section III-A (Fig. 3). When signals are not influenced by any perturbation, all control strategies show nice convergence behavior and smooth control inputs. Compared to FOFF and SO, FO presents a deviation in convergence from the ideal evolution, justified by the intermediate velocity regulation not being able to perfectly track the desired command q 1 . FOFF and SO are instead able to better follow the desired behavior and the error evolves according to the expected second order exponential decay. The difference in these two controllers can be seen only in the control input, yet not significantly different. When noise is considered, more evident differences appear. FO's convergence becomes worse, with a little overshoot near to the equilibrium configuration, whereas FOFF and SO still feature a nice evolution of the error towards zero. In terms of control inputs, all strategies suffer from the added noise. In particular, the oscillations corresponding to FOFF are quite large, with peak-to-peak amplitude of about 30 N m even at equilibrium. Oscillations of FO and SO controllers are bounded instead to a reasonable amount, between 0.5 N m and 1 N m of amplitude in both cases. It is possible to enhance the convergence of FO by increasing the gain γ or selecting a faster behavior (higher ω n ), at the cost of higher oscillations in the commanded efforts. Similarly, to reduce the shakiness of the control signal in FOFF, the response has to be penalized with lower values of γ and/or ω n .

It is also worth mentioning that when noises are of higher magnitude the relative performances of these controllers are not significantly altered. FOFF and SO still give satisfactory results in terms of error convergence, and the oscillations in the control commands produced by FO and SO are quite smaller compared to those introduced by FOFF.

D. Tracking of a Visual Trajectory

In this section, we consider the case of varying desired features, in which case ṡ and s are non-zero. For each normalized point s

i = [x i , y i ]
T , the motion is defined as a circular trajectory starting at the instant t 0 and ending after a time period T :

x i (t) = x i (t 0 ) + R (cos ρ(t -t 0 ) -1) y i (t) = y i (t 0 ) + R sin ρ(t -t 0 ) (26)
In this definition, ρ(t) is a quintic polynomial whose value changes from 0 to 2π in T = 2 s, such that ρ(0) = ρ(T ) = 0 and ρ(0) = ρ(T ) = 0. We assume the desired initial position of each point s i (t 0 ) to be the one attained after the regulation task described in the previous section. The amplitude R is set to 0.2, and thus the radius of the physical trajectory is expected to be of 0.1 m (the sensor should be at 0.5 m from the four points, as explained before). In addition, we assume t 0 0, so that the homogeneous term in (23) does not affect the control, as the desired velocity is already null at t 0 .

Results are depicted in Fig. 4, which shows that FO is the worst controller in terms of tracking error. FOFF again presents high oscillations in the input command during the motion, but still with good tracking performances. Finally, SO gives the best results, with tracking errors of the same magnitude of FOFF and reduced effort oscillations.

IV. CONCLUSIONS

In this work, we reviewed and compared simple Visual Servoing control schemes based on both first and second order models. It has been shown that for a redundant manipulator it is desirable, if not necessary, to consider a secondary task in the case the Visual Task is formulated at the acceleration level in order to fully constrain the degrees of freedom and prevent uncontrolled motions spanning the kernel of the Feature Jacobian.

We showed via simulations that controllers can be designed both from velocity and acceleration commands to obtain a very similar evolution of visual features error. In particular, FOFF and SO were equally able to ensure good performances in terms of error regulation. However, results suggest that schemes based on second-order models are preferable in order to obtain smoother control efforts.

In this analysis, redundancy was not exploited to truly achieve any secondary objective. We believe that with further improvements it will be possible to enhance existing strategies by, e.g., minimizing the efforts along the trajectory while ensuring Visual Tasks to be properly executed. For this purpose, a possible future direction of work could consider the use of Model Predictive Controls.

We are working on a real implementation of the proposed strategies, in order to validate with experiments our simulated results. We are especially interested in performing Servoing tasks at higher speed, and to study more in depth the limits of FOFF and SO in such context.

Fig. 1 .

 1 Fig. 1. Initial (left) and desired (right) configuration of the observed object. The points used for the task are the centers of the four circles.

Fig. 2 .

 2 Fig. 2. Comparison of the visual error (left) and control input (right) using the different control strategies, in absence of noise. FOFF and SO are almost indistinguishable in terms of error evolution.

Fig. 3 .

 3 Fig. 3. Comparison of visual error (left) and control inputs (right) when noise is injected in simulations. FOFF and SO are still almost indistinguishable in terms of error evolution.

Fig. 4 .

 4 Fig. 4. Comparison of the visual error (left) and control input (right) during the execution of a trajectory. Time axes are relative to t 0 .

Obviously, (12b) can be obtained by differentiation of (12a) as well.

ACKNOWLEDGMENT

This work was carried out in the framework of the PROMPT project, a project funded by RFI Atlanstic 2020.


Christophe Sauret

Vanel Steve

Siyou Fotso

Joseph Bascou

Hélène Pillet

Engelbert Mephu-Nguifo

Pascale Fodé

Philippe Vaslin

V Siyou Fotso

E Mephu

N Guifo

Cluster analysis to investigate biomechanical changes during learning of manual wheelchair locomotion: a preliminary study

Keywords: manual wheelchair, learning, biomechanics, cluster analysis Cluster B

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Learning a new motor task is characterized by improvements in movement smoothness and kinematics repeatability [START_REF] Shmuelof | How is a motor skill learned? Change and invariance at the levels of task success and trajectory control[END_REF]. Hence, demonstrating biomechanical changes during the early phase of learning can be difficult due to high standard deviations that directly affect statistical analysis. Especially for curve analysis, timing variability results in wide corridor of normality. In addition, in the early phase of learning, each trial participates to the learning process, compromising the collection of redundant data. Besides, for further biomechanical analysis, curves selection/exclusion is often required, and is generally done manually based on a single parameter. The aim of this study was to investigate the interest of cluster analysis (i) to describe biomechanical changes during the early phase of learning manual wheelchair (MWC) locomotion and (ii) to define representative time courses of the locomotion parameters at different steps of the leaning process.

Materials & Methods

One young male able-bodied subject without any experience in MWC locomotion participated in this study. He was asked to train himself without indication at low intensity during 120 minutes distributed within 3 weeks. Two experimental sessions were performed allowing to analyze the 15 first and last minutes using an instrumented wheelchair (FRET-2) that allowed both the measurement of the mechanical actions applied on the handrims and the determination of the MWC velocity. The measurement session was repeated 6 months later without any practice. Each measurement session was composed of 3 blocks of 5 minutes, including about 10 start-ups each (84 start-ups in all). Data from the FRET-2 were used for cluster analysis based on k-mean algorithm [START_REF] Esling | Time-series data mining[END_REF], considering 6 parameters: total (P1) and push phase (P2) durations, total impulse of both (left and right) propelling torques (P3) and its distribution between the resultant of both propelling torques (P5), and linear velocity at the end of the push phase (P6). Because involved parameters did not have the same units and in order to attribute the same weight to all of them, data of each parameter were centred and reduced. Finally, the number of clusters was arbitrarily fixed to 8 and named from A to H.

Results

The 8 clusters gathered from 3 to 22 trials each. Main clusters where clusters B, D, E and H, gathering 14, 22, 11 and 12 trials, respectively. Time courses of the resultant of both (right & left) propelling torques for clusters B, D, E and H are depicted in Figure 1 and their features are reported in Table 1. The 1st block of the 1st session of practice gave the following sequence: C-B-B-B-B-C-B-F-F-F where each letter refers to the cluster to which belonged each startup: in this block, start-ups had mainly the features of clusters B and F. During the 3rd block of the 1st session, the sequence was B-D-B-F-D-F-E-D-G showing the apparition of start-ups of cluster D. During the first session (3 blocks), start-ups mainly belonged to cluster B (44%), whereas they mainly belonged to clusters E (35%) and H (31%) during the second session, and to cluster D (55%) during the session carried out after 6 months without any practice.

B D E H P 1 (s) 1.8 ± 0.1 1.7 ± 0.1 1.5 ± 0.1 1.6 ± 0.1 P 2 (s) 1.2 ± 0.1 1.1 ± 0.1 1.0 ± 0.1 1.2 ± 0.1 P 3 (Nm.s) -27 ± 3 -34 ± 2 -38 ± 2 -38 ± 2 P 4 (%)
40 ± 5 39 ± 3 37 ± 3 42 ± 3 P 5 (Nm) -40 ± 5 -52 ± 4 -64 ± 3 -59 ± 4 P 6 (m/s) 0.7 ± 0.1 0.9 ± 0.1 1.0 ± 0.1 1.0 ± 0.1

Table 1: Features of the four main clusters. Two distances between clusters have been computed: inter-centroids distance and minimal distance between the nearest points of two different clusters. These distances allow identifying what were the clusters that most resembled while considering all parameters. These distances showed that Cluster E and H were the most *Corresponding author. Email: christophe.sauret@ensam.eu resembling. Cluster B was closer to cluster D than to clusters E and H but cluster D resembled more to cluster E and H than to cluster B. So, cluster D can be seen as intermediate between cluster B on the one hand and clusters E and H on the other hand. Finally, using a more classical approach, standard deviations computed for each block of practice were generally from 1.5 to 2.0 times higher than within clusters, representing up to 30% of the nominal value (10% in average). In the same way, average time courses calculated for each block of practice could be significantly different from clusters means. For instance, curves obtained during the first and third blocks of the first session of practice were clearly different from those of the clusters (B and D) to which they mainly belonged (Figure 2).

Discussion

The symbolic representation method (Siyou Fotso et al., submitted) used in this study for classifying and analysing start-up pushes showed that even if large differences were observed between curves, no single trial was secluded in a cluster and thus considered as an outlier.

Cluster occurrence frequency allowed identifying the main clusters, which were depending on the expertise level. However, for the studied subject, features of cluster D, which was the most frequent cluster during the 3rd session (after 6 months), need further investigations to identify if the decrease of the dynamics properties (torque and impulse), compared to the 2nd session (after 120 minutes of training) where clusters E and H were the most frequent, was due to a lower selfselected intensity or to a loss of motor ability. The use of a symbolic representation (i.e. string) gives an easier readability of the strategies used by a subject to perform a task than the comparison of individual trials, especially when considering several parameters simultaneously (here: phase durations, torques, impulses and velocities).

The comparison with a more classical approach showed that cluster analysis yielded lower standard deviations that would allow a better definition of the different strategies and to obtain some significant statistical differences between motor patterns used during the learning process. However, one limitation of this method is that the number of strategies is defined by the number of clusters chosen as input in the algorithm. If some methods already exist to estimate the best number of clusters based on different criteria [START_REF] Sugar | Finding the number of clusters in a data set: An information theoretic approach[END_REF], their appropriateness with the present research question remains to be evaluated.

Conclusion

This preliminary study showed the potential of cluster analysis to study biomechanical changes during the learning of MWC locomotion. The interest of this processing method relies on both its direct use for locomotion analysis and its integration into the whole data processing workflow. It is, thus, able to provide representative data prior to further biomechanical analysis such as the use of musculoskeletal models, for instance. However, further work is still required before it could be used in an efficient and generalized way.

Figure 1 :

 1 Figure 1: Average time courses of the resultant of both propelling torques for clusters B, D, E and H.

Figure 2 :

 2 Figure 2: Average time courses of the resultant of both propelling torques for clusters B and D; and for the 1 st and the 3 rd blocks of practice.


Joseph Bascou

Christophe Sauret

Coralie Villa

François Lavaste

Hélène Pillet

Measurement of wheelchair adjustment effects on turning deceleration

Keywords: manual wheelchair, rotation, parameters

come   L'archive ouverte pluridisciplinaire

Introduction

Manual wheelchair (MWC) locomotion combines straightforward and turning motions, in everyday life as well as in sport practice. Many authors demonstrated the effects of various MWC properties, such as geometry or wheel type, for straightforward displacements [START_REF] Brubaker | Wheelchair Prescription: An Analysis of Factors That Affect Mobility and Performance[END_REF][START_REF] Medola | Users ' Perceptions on Mobility , Comfort and Usability of Manual Wheelchairs[END_REF], while only few studies have investigated their influence for turning motion [START_REF] Bascou | Turning Resistance of a Manual Wheelchair: A Theoretical Study[END_REF][START_REF] Caspall | Changes in Inertia and Effect on Turning Effort across Different Wheelchair Configurations[END_REF][START_REF] Kauzlarich | Wheelchair Caster Shimmy and Turning Resistance[END_REF]. In particular, the impact of wheelchair setup on its turning deceleration, which characterizes the MWC tendency to stop its turning motion, is unclear. This study aims at clarifying the effects of MWC adjustments on turning deceleration in the field, using a fractional factorial design.

Methods

Materials

An inertial measurement unit (MTi, X-sens, The Netherlands) was placed on the frame of a sport MWC whose front wheels were removed and replaced by a custom fork device (Figure 1). The MWC was loaded with 40 kg of additional mass, fixed on the seat. The fore-aft location of the additional mass and the custom fork allowed changing the fore-aft location of the total centre of mass (COM) with respect to the rear wheel axle (factor "A", ranging from 0.04m to 0.07m), the inclination of the fork hinge (factor "B, ranging from 0° to +3°), the fork trail distance (factor "C", ranging from 0.03m to 0.08m), the caster wheel diameter (factor "D", ranging from 0.06m to 0.08m) and the location of the fork hinge with respect to the rear wheel axle (factor "E", ranging from 0.50 m to 0.63m).

Fractional factorial design

8 sets of 5 turning deceleration tests were performed to define the influence of the MWC setup on its turning behaviour. In each trial an experimenter initiated the MWC rotation in clockwise direction and let the MWC turning freely during about 1 second. For each set, the 5 previously described MWC parameters (factors A to E) where changed according to a fractional factorial design [START_REF] Taguchi | System of Experimental Design: engineering methods to optimize quality and minimize costs[END_REF] with 2 levels (Table 1). The factorial design experiments provided a model allowing the assessment of the MWC angular deceleration according to its settings:

0 A A B B C C D D E E AB A B AD A D w a a l a l a l a l a l a l l a l l        
where w is the MWC angular deceleration, l A , l B , l C , l D , l E , are the level values (between +1 and -1) for the parameters A to E, and a A , a B , a C , a D , a E their respective effects. Second and third order interactions between factors were neglected, except for interactions AB and AD (with effects a AB , and a AD ). In order to evaluate the impact of possible ground inclination, set 1 was repeated (set 1b) in counterclockwise direction. Two additional sets (sets 10 and 11) were performed for model verification.

1 -1 -1 -1 -1 1 1 -1 2,8 (±0,4) 2 -1 -1 1 1 -1 -1 -1 1,8 (±0,1) 3 -1 1 -1 1 1 -1 1 2,1 (±0,1) 4 -1 1 1 -1 -1 1 1 3,0 (±0,2) 5 1 -1 -1 1 -1 1 -1 3,8 (±0,3) 6 1 -1 1 -1 1 -1 -1 4,2 (±0,3) 7 1 1 -1 -1 -1 -1 1 4,2 (±0,3) 8 1 1 1 1 1 1 1 4,0 (±0,1) 1b -1 -1 -1 -1 1 NA NA 2,5 (±0,2) 10 1 -1 1 1 1 NA NA 3,6 (±0,2) 11 1 1 1 1 1 NA NA 4,1 (±0,2)

Data treatment and hypothesis

The MWC angular velocity was computed assuming a planar movement: the rotation velocity was the resultant of the angular velocity measured by the Inertial Measurement Unit. The angular deceleration was obtained by time differentiating the MWC angular velocity.

Results and discussion

During the free turning phase, the MWC angular deceleration was constant (for every selected trials, the angular velocity followed a decreasing line with a correlation coefficient r² superior to 0.98), which supports the hypothesis that deceleration was not linked to MWC angular rotation velocity during the free turning phase and was directly linked to MWC settings.

Mean deceleration value a 0 was 3.2 rad/s 2 . Main effects were attributed to the fore-aft location of the total COM and the diameter of the front casters (Figure 2). The effect of fork angle and fore-aft position was significant, but represented less than one fourth of the caster diameter and on tenth of the total COM position. The fork trail had a very low effect in this movement as the rotation was already initiated, but should have an effect in the rotation initiation. The total COM position had a high effect as it changed both the total inertia and load distribution on the front and rear wheels, which is in accordance with previous theoretical results [START_REF] Bascou | Turning Resistance of a Manual Wheelchair: A Theoretical Study[END_REF]. The user can have a direct action on this parameter by modifying his posture, leaning forward to decrease (/ backward to increase) the rotation.

Figure 2 Effects of MWC settings on its deceleration

The wheel diameter also had a significant effect, possibly due to a variation of rolling (linked to the wheel radius) and swiveling resistances. This assumption is supported by the non negligible value of interaction between total COM position and the front wheel diameter (interaction AD).

Turning the MWC in counter-clockwise direction (trial 1b) resulted in a 8% difference with respect to the clockwise direction (trial 1), potentially due to a slight ground inclination.

Conclusions

This study allowed the experimental classification of the effects of various MWC settings on its angular deceleration in the field. The results underlie the importance of total COM position and the choice of front caster diameter. Further work should be conducted, particularly to assess the wheel swivelling resistance contribution to the MWC manoeuvrability and to cancel the effect of ground inclination.

Figure 1 :

 1 Figure 1: MWC, custom fork and additional masses

Table 1

 1 MWC settings during experimental plan and resulting angular deceleration

*Corresponding author. Email: joseph.bascou@invalides.

fr Trial/ factor A B C D E AD AB Angular deceleration rad/s²(±SD)

  Comparing the model results with the experimental ones (sets 10 and 11) resulted in 8% and 7% errors respectively, which is acceptable considering a fractional factorial design model. The hypothesis of negligible effects of AC and AE interactions was validated a posteriori considering the low effects that were observed for factors C and E.

Acknowledgements

The authors thank the INI/CERAH design and fabrication departments for the realization of the fork device and S.Houdet, elite MWC tennis player, for the loan of his sport MWC.


Van-Phu Ha
email: van-phu.ha@inria.fr

Tomofumi Yuki
email: tomofumi.yuki@inria.fr

Olivier Sentieys
email: olivier.sentieys@inria.fr

Noise Budgeting in Multiple-Kernel Word-Length Optimization

Keywords: Fixed-Point refinement, Word-Length Optimization, Multiple Kernel Optimization

Word-Length Optimization (WLO) is a key step in the implementation of Digital Signal Processing applications on hardware platforms. Existing approaches face scalability problems for applications with several kernels. In this paper, we present our work-in-progress to address the scalability problems when performing multi-kernel WLO. Our approach uses application-wide analysis to derive noise budgets for each kernel, followed by independent WLO. The main idea is to characterize the impact of approximating each kernel to the accuracy/cost through simulation and regression analysis. The constructed models can be used to determine the appropriate noise budget for each kernel. When applied to WLO for two kernels in Image Signal Processor, the noise budgets given by our analysis closely match those found empirically with global WLO.

I. INTRODUCTION AND MOTIVATION

Fixed-Point arithmetic is widely used for the implementation of Digital Signal Processing (DSP) systems on electronic devices. The Float-to-Fix conversion in hardware design flow always demands to optimize the word-length of variables in the system to find a good trade-off between the cost and quality requirement. The process of word-length determination is typically called Word-Length Optimization (WLO). WLO is known as an NP-Hard problem with the complexity growing exponentially when more variables are involved. In fact, this process may take up to 25-50% of design time [START_REF] Clark | Accelerating fixed-point design for mb-ofdm uwb systems[END_REF]. Therefore, WLO is still a problem of interest in reducing the time-tomarket in the electronics industry.

A common approach to WLO involves iterative search algorithms based on simulations [START_REF] Cantin | A comparison of automatic word length optimization procedures[END_REF]- [START_REF] Lee | Accuracy-guaranteed bit-width optimization[END_REF]. Others have proposed relaxing the problem with convex functions to directly solve for the optimal solution [8]- [START_REF] Parashar | A polynomial time algorithm for solving the word-length optimization problem[END_REF]. However, these methods cannot be directly extended to handle applications with several kernels. The accuracy evaluation by simulations is nearly impossible because the exponential growth of the problem leads to a huge number of simulations. Methods based on convex optimizations require the accuracy to be modeled as convex functions, which cannot be analytically constructed in general. In particular, existing analytical models focus on modeling the noise power, and cannot be used for other quality metrics such as Structural Similarity (SSIM). Recently, some techniques to address the scalability of WLO have been proposed [START_REF] Chung | Bit-width optimization by divide-and-conquer for fixed-point digital signal processing systems[END_REF]- [START_REF] Parashar | Accelerated performance evaluation of fixed-point systems with un-smooth operations[END_REF]. The work by Novo et al. [START_REF] Novo | Cracking the complexity of fixed-point refinement in complex wireless systems[END_REF] first performs local WLO with different accuracy constraints and later combine the local solutions with an iterative search. The others have explored hierarchical decomposition of the original system to reduce the problem size at each step in the exploration [START_REF] Chung | Bit-width optimization by divide-and-conquer for fixed-point digital signal processing systems[END_REF], [START_REF] Parashar | Accelerated performance evaluation of fixed-point systems with un-smooth operations[END_REF]. These work employ a large number of local WLO, which limits the scalability when considering applications with several kernels.

II. SPECIFIC PROBLEM AND SOLUTION

In this study, we discuss an approach to overcome scalability issues when applying WLO for applications with several kernels. We use Image Signal Processor, a post-processing pipeline for digital cameras illustrated in Figure 1, as our running example. The ISP consists of several kernels such as noise cancelling (Non-Local Means), linearization recovering, lightfall-off correcting (Vignetting Correction), color adjusting and quality improving by the remaining kernels. The challenges in performing WLO for applications with several kernels include:

• The number of operations is too high to perform global optimization across all kernels. • The kernels often contain diverse topologies, e.g., including non-linear structures and unsmooth operators, making it difficult to analytically model the impact on quality. • Analytically modeling quality metrics, other than noise power, is difficult. For instance, in ISP, the commonly used metric is SSIM, which is hard to model analytically. In this work, we view an application as a composition of kernels; a small enough unit of computation where local WLO is feasible. For our ISP example, the functional components shown in Figure 1 is directly used as kernels, but the granularity of a kernel can be arbitrarily changed. The loss of quality at the outputs of each kernel after global optimization indicates the amount of approximation that should be applied for each kernel to obtain the minimum cost. We call this amount of approximation to be applied for each kernel as noise budget. However, performing global optimization has scalability issues, and is not feasible for large application. We assume that performing local optimization to each kernel using the noise budget found in global solution gives a similar solution. Then, finding the noise budget without global optimization allows the global optimization problem to be simplified into local optimization of each kernel. For example, consider that we approximate two kernels, NLM and Unsharp, in the ISP with the target quality being SSIM ≥ 0.99. Applying global optimization, which is feasible for two kernels, gives a solution with the output quality being 0.99 SSIM. Applying the wordlengths in this solution independently to the two kernels gives output quality 0.9990 for NLM and 0.9923 for Unsharp. These are the noise budget given by the global optimization, which we would like to predict. Our approach consists of 3 steps described below:

1) The relation between cost and quality is found through exploration and regression. Each kernel is explored once with a heuristic search targeting high output quality to obtain cost and quality data samples. The cost is estimated from wordlength configurations based on the cost of each operator determined by synthesizing the operators targeting ASIC. These collected data points are used as inputs to non-linear regression with the function a × x b + c as a template. This gives us a continuous function that directly relates quality and cost without going through wordlengths. The use of power function in the template is motivated by the non-linear behavior of SSIM with respect to the wordlengths. 2) The noise budgets that gives minimal cost are found by performing constrained optimization over the functions.

3) The noise budgets found are used as constraints for independent WLO of each kernel. Then, the results are combined to form the final solution.

Our approach is scalable because (i) both the exploration and WLO are independent to each kernel, and hence has a linear complexity with respect to the number of kernels; and (ii) the functions of the constrained optimization are convex and hence convex optimization tools can be used.

Our preliminary result is summarized in Table I. We first perform local WLO using Tabu search [START_REF] Nguyen | Novel algorithms for word-length optimization[END_REF] to collect data points for regression. Figure 2 illustrates the constructed models and the data points from initial exploration. Note that only Pareto optimal points are used for regression. The predicted noise budgets are used to perform another local WLO for each kernel and the solutions are combined.

We are able to predict the noise budgets that are close to those inferred from global search results. The execution time is not significantly shorter than the global search for our experiments with two kernels, but we expect the difference to increase when more kernels are included in the exploration.

III. CONCLUSION

WLO in the context of large applications remains challenging. In this paper, we present our preliminary work towards improving the scalability of WLO in such contexts. The key idea in our approach is to analyze accuracy models constructed from sampling the design space to derive noise budgets for each kernel. Then, each kernel may be optimized independently, avoiding the combinatorial explosion of the search space as larger applications are considered.

The preliminary results indicate that the approach is feasible and effective for the ISP example. We are currently working on experimenting with approximating more than the two kernels in ISP, other examples, and different accuracy constraints to further evaluate our approach. we are also exploring how the model construction time can be further reduced.

Fig. 1 .

 1 Fig. 1. Image Signal Processor

TABLE I THE

 I GLOBAL SOLUTION VS. OUR SOLUTION BY CONSTRAINTS GIVEN BY THE NOISE BUDGETING PROCESS, AT REQUIRED SSIM = 0.99

	Kernel		SSIM		Cost (µm 2 )	Exec. Time (minutes)
		Ref.	Ours (Predicted)	Ours (Measured)	Ref.	Ours Ref	Ours
	NLM	0.9936	0.9961	0.9978	11570 12576 -	44
	Unsharp 0.9923	0.9939	0.9939	8251 8317 -	31
	Both	0.9900	-	0.9920	19821 20893 229 75 (+89)
		1
			NLM (Not used)
		0.9	NLM (Pareto-optimal)
	maximum value)	0.5 0.6 0.7 0.8	Unsharp (Not Used) Unsharp(Pareto-optimal) NLM Model Unsharp Model
	by
	(normalized	0.2 0.3 0.4
	Cost	0.1
		0
		0.985	0.99		0.995		1
				ISP SSIM
	Fig. 2. Model of Cost-SSIM interaction for NLM and Unsharp


Oumaima Matoussi

Yves Durand

Olivier Sentieys
email: olivier.sentieys@inria.fr

Anca Molnos
email: anca.molnos@cea.fr

Error Analysis of the Square Root Operation for the Purpose of Precision Tuning: a Case Study on K-means

Keywords: approximate computing, error analysis, roundoff error, algorithmic deviation, square root, Newton Raphson method, precision tuning, k-means, clustering, floating point

In this paper, we propose an analytical approach to study the impact of floating point (FLP) precision variation on the square root operation, in terms of computational accuracy and performance gain. We estimate the round-off error resulting from reduced precision. We also inspect the Newton Raphson algorithm used to approximate the square root in order to bound the error caused by algorithmic deviation. Consequently, the implementation of the square root can be optimized by fittingly adjusting its number of iterations with respect to any given FLP precision specification, without the need for long simulation times. We evaluate our error analysis of the square root operation as part of approximating a classic data clustering algorithm known as K-means, for the purpose of reducing its energy footprint. We compare the resulting inexact K-means to its exact counterpart, in the context of color quantization, in terms of energy gain and quality of the output. The experimental results show that energy savings could be achieved without penalizing the quality of the output (e.g., up to 41.87% of energy gain for an output quality, measured using structural similarity, within a range of [0.95,1]).

I. INTRODUCTION

The ever increasing volume, diversity and high dimensionality of data goes hand in hand with the rapid growth of energy consumption of computer systems. Thus, improving the energy efficiency of computer systems that try to keep pace with the constant growth of information is a critical concern.

In many application domains that deal with huge amounts of data such as multimedia processing (images, audio, video, etc.), data mining and machine learning, computations can be tolerant to some degree of error without critical degradation in the quality of the output. For example, a small image quality loss due to some modification of the color of a group of pixels can be hardly noticed by the user due to the limited capabilities of human perception. In such classes of applications, output accuracy could be traded for energy reduction.

Performance gains can be achieved at the application level using software techniques, thanks to a relatively new computing paradigm called approximate computing [START_REF] Xu | Approximate computing: A survey[END_REF]. Introducing inexactness in computations may lead to energy scaling with little to no loss in accuracy. A way to introduce inexactness in an application is by precision reduction of FLP variables and computations [START_REF] Ho | Efficient floating point precision tuning for approximate computing[END_REF], [START_REF] Barrois | Customizing Fixed-Point and Floating-Point Arithmetic -A Case Study in K-Means Clustering[END_REF]. Finding the sweet spot between the numerical accuracy and the performance of the application (e.g. energy consumption) by refining the bit-width of the FLP variables is known as precision tuning.

A loss in computational accuracy is inevitable due to reduced precision. However, to be able to determine the optimal precision (i.e. bit-width) of FLP variables that minimizes energy cost with the least impact on computational accuracy, two approaches can be considered: FLP simulation and analytical techniques to track round-off error (i.e. error introduced due to limited precision). FLP simulation consists in using libraries (e.g. MPFR) that allow the definition of FLP variables with adjustable bit-widths and provide the appropriate arithmetic operations. Nonetheless, the search space of the optimal bitwidth can be very large, in that every change in the format of a FLP variable requires a new simulation, which can be very time consuming.

Analytical approaches, on the other hand, try to determine a mathematical formula that models the impact of quantization error on the accuracy of the output [START_REF] Rocher | Analytical approach for numerical accuracy estimation of fixed-point systems based on smooth operations[END_REF], [START_REF] Parashar | Analytical approach for analyzing quantization noise effects on decision operators[END_REF]. Once this formula is established, it can be applied to any application with different FLP formats. Unfortunately, this only works with smooth operations (arithmetic operations like addition, subtraction, multiplication, division, etc.). So, analytical approaches help save the time spent on tuning FLP precision but fail to estimate the computational error in the presence of non-smooth operations. In this case, simulation becomes inescapable. Thus, we advocate through our work the combination of simulation and analytical techniques.

A plethora of work using simulation-based approaches for the purpose of FLP precision tuning was propounded [START_REF] Barrois | Customizing Fixed-Point and Floating-Point Arithmetic -A Case Study in K-Means Clustering[END_REF], [START_REF] Roy | Asac: Automatic sensitivity analysis for approximate computing[END_REF], [START_REF] Rubio-González | Precimonious: Tuning assistant for floating-point precision[END_REF], etc. As for analytical approaches, efforts were made in the context of fixed-point word-length optimization [START_REF] Rocher | Analytical approach for numerical accuracy estimation of fixed-point systems based on smooth operations[END_REF], [START_REF] Parashar | Accelerated performance evaluation of fixed-point systems with un-smooth operations[END_REF] and mainly focused on common arithmetic operations like addition and multiplication. Square root operation is significantly more costly, in terms of energy consumption, than addition or multiplication and requires several iterative cycles to complete. Moreover, unlike simple arithmetic operations, the square root is usually implemented using the Newton Raphson method, which is an approximation of the operation itself, and adds another type of error referred to as algorithmic deviation besides the round-off error. Consequently, error analysis of a square root operation is not a straightforward task and entails a twofold method.

The main contribution of this paper consists in conducting static error analysis of the square root operation and estimating a bound on the errors caused by limited precision and algorithmic deviation. Based on the formally derived error bound, the number of iterations of the square root operation is automatically adjusted for each precision. Therefore, the implementation of the square root (in terms of number of iterations) can be optimized with respect to any given FLP precision specification without the need for long simulation times.

In addition to the stand-alone analysis of the square root operation and for a more complete evaluation of our proposition, the second contribution of this paper consists in evaluating our analytical approach by studying the square root operation as part of a whole application. Our choice fell on a data clustering algorithm called K-means a.k.a. Lloyd's algorithm [START_REF] Hamerly | Accelerating Lloyd's Algorithm for k-Means Clustering[END_REF]. It is a type of unsupervized learning algorithm used to cluster a set of unlabeled data into k clusters based on data feature similarity [START_REF] Silva | Data stream clustering: A survey[END_REF]. The similarity is typically determined using the Euclidean distance measure, which involves square root operations. The distance function is the mainstay of the majority of clustering and classification algorithms, which will help us highlight the impact of square root optimization on the performance of the application as a whole.

It should be noted that the proposed analytical approach for the square root operation is independent of the application itself. Any clustering algorithm or any other application (e.g. digital signal and image processing, 3D graphics, spectrum analysis, wireless communications, etc.), for that matter, where square root computations can be found, would also benefit from our error analysis method. The importance of our contribution lies within the integration of our analytical approach, which is aimed for a smooth operation, namely the square root, in an application that also contains non-smooth operations. Hence, the combined effort of simulation and analytical results is needed for precision tuning. We quantify the efficiency of the proposed approach by measuring the quality degradation of K-means' output as well as its energy gain.

The rest of this paper is organized as follows. After overviewing contributions dealing with approximate computing in Section II, we detail our error analysis approach of the square root operation in Section III. We explain the process of sensitivity analysis of K-means, based on our square root analysis results and variable-precision FLP simulation in Section IV. We validate the efficiency of approximate Kmeans by measuring both SSIM and energy gain and we also discuss the experimental results in Section V, before concluding the paper in Section VI.

II. RELATED WORK

Precision tuning is a research direction that has received significant attention and is motivated by the fact that FLP operations contribute to the energy footprint of an application.

In [START_REF] Roy | Asac: Automatic sensitivity analysis for approximate computing[END_REF], a framework called ASAC using statistical methods was proposed. This framework helps discover approximable and non-approximable program parts through sensitivity analysis. The main idea is to perturb program variables and check the resultant output against a correct output, i.e. one fulfilling an acceptable QoS threshold.

Precimonious [START_REF] Rubio-González | Precimonious: Tuning assistant for floating-point precision[END_REF], a dynamic program analysis tool, was proposed to assist developers in choosing the lowest precision that satisfies accuracy and performance constraints. Given a set of FLP variables and their possible types, a search based on the delta-debugging algorithm is performed. The search outputs a type configuration that maps each variable to a type. However, Precimonious requires a representative set of program inputs provided by the user. If the same type configuration is applied on a much worse conditioned input then no guarantees can be made.

The method proposed in [START_REF] Ho | Efficient floating point precision tuning for approximate computing[END_REF] is not limited to determining the best mix of single or double precision in a program as in [START_REF] Rubio-González | Precimonious: Tuning assistant for floating-point precision[END_REF], but it computes the precision of FLP variables down to the bit level of the mantissa. A heuristic precision tuning algorithm based on a binary search is performed in order to find the smallest precision possible for each variable while keeping the output error within some user-given bound. The results show that the speedup gained by tuning some programs may be diminished or even eliminated because of the overhead due to data type conversions performed by the compiler.

In [START_REF] Barrois | Customizing Fixed-Point and Floating-Point Arithmetic -A Case Study in K-Means Clustering[END_REF], hardware implementations of fixed-point and FLP arithmetic operators are compared in terms of area, delay and energy. A custom FLP library called ct float was devised to vary the bit width of the variables. The authors concluded that FLP operators provide a better energy/accuracy tradeoff for small bitwidths but important area, delay and energy overhead for larger bitwidths, compared to fixed-point operators.

The common denominator in the majority of FLP precision tuning approaches is that they rely on varying the precision of the different FLP variables in an application (usually with the help of a multiple precision FLP library) and testing different combinations before deciding on the optimal precision. This process is prohibitively expensive, time wise. Furthermore, when big applications with a huge set of input is under optimization, the search space for optimal precision can be enormous and it might be impossible to cover fully.

Analytical approaches are faster than simulation because once a mathematical expression that models (or bounds) the round-off error is determined, it can be applied to different precisions. So, it is a one-time effort whose results can be applied to multiple FLP formats. Efforts have been made to find analytical models in the context of accuracy evaluation of fixedpoint systems. A series of propositions of analytical methods to measure the output noise of signal processing systems were presented in [START_REF] Rocher | Analytical approach for numerical accuracy estimation of fixed-point systems based on smooth operations[END_REF], which dealt with smooth operations, then in [START_REF] Parashar | Analytical approach for analyzing quantization noise effects on decision operators[END_REF], which catered for decision operators and recently in [START_REF] Parashar | Accelerated performance evaluation of fixed-point systems with un-smooth operations[END_REF], where a hybrid approach based on both simulation and analytical results was highlighted.

A parametric error analysis of a version of Goldschmidt's square root algorithm that uses directed rounding implemen-tations, was proposed in [START_REF] Seidel | A parametric error analysis of goldschmidt's square-root algorithm[END_REF]. The error analysis is based on relative errors of intermediate computations. The proposed error formulae were intended to help determine the optimal hardware implementation (i.e. multiplier dimensions). However, their analysis was not demonstrated on a real architecture.

In our work, we study analytically the square root operation, which is the backbone of many algorithms such as clustering, data mining and signal processing algorithms, with respect to FLP precision tuning. Moreover, we quantify the efficiency of the proposed error bound in the context of K-means by comparing the QoS of the inexact version of K-means to its exact counterpart. We also measure relative energy gains using energy costs obtained from [START_REF] Barrois | Customizing Fixed-Point and Floating-Point Arithmetic -A Case Study in K-Means Clustering[END_REF], [START_REF] Vasilakis | An instruction level energy characterization of arm processors[END_REF] and profiling information.

III. ANALYSIS OF THE SQUARE ROOT OPERATION

First, to eliminate any confusion, it should be noted that the terms accuracy and precision are not used interchangeably in this paper. Nearly all processors and programming languages support FLP numbers, which are defined in the IEEE-754 normalization. A FLP number f is represented by an exponent e, a mantissa m and a sign bit s: f = (-1) s × m × 2 e . We designate by precision the number of bits used to represent the mantissa of a FLP variable, whereas we refer by accuracy to the quality of the result, i.e., the degree of error in the output, compared to a reference or a golden output.

The square root operation is usually approximated with the Newton Raphson method. To compute y = √ a, a > 0, the Newton Raphson method starts with an initial guess (i.e. initial seed value) y 0 > 0. The initial guess is then refined by iterating over:

y n+1 = 1 2 (y n + a y n )
However, this conventional iteration for the square root computation is not frequently used because it entails a division at each step. Since division is generally much slower than multiplication, the square root reciprocal is usually advocated. The square root reciprocal converges to 1 √ a and iterates over:

x n+1 = xn 2 × (3 -a × x 2 n ).
In order to get √ a, the result is multiplied by a. To optimize the implementation of the Newton Raphson method, we aim at adjusting the number of its iterations according to the selected precision, without jeopardizing the accuracy of the result. To do so, we perform error analysis to statically determine at which iteration it is preferable to stop the computations, for a specific precision p, without causing further error. The two main causes of error that we investigated are:

• the round-off error caused by FLP representation and FLP operations, • and the systematic error (a.k.a. unavoidable error or algorithmic deviation), which results from the Newton Raphson approximation itself.

A. Bounding the Round-off Error

Let f l(x+y) = (x+y)(1+ add ) and f l(x×y) = (x×y)(1+ mul ), be approximations of the exact mathematical operations

+ and × respectively, ∀x, y ∈ R, such that | add,mul | ≤ m .
m , called machine epsilon or the unit roundoff, is defined as the smallest number such that 1 + m > 1 [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF]. m = β -p , where β is the base and p is the number of bits used for the magnitude of the mantissa of a FLP number represented as r 0 .r 1 r 2 ...r p-1 × β e . We assume that add = mul = and that multiplication by 1 2 does not incur round-off error. ∀x i ∈ R, let xi = f l(x i ) be the FLP representation of x i ; i = 1..n+1. Machine epsilon is an upper bound on the relative error in representing a FLP number: |xi-xi| |xi| ≤ m . We would like to bound the round-off error of the Newton Raphson reciprocal. To do so, we aim at expressing xn+1 as: xn+1 ≤ x n+1 × (1 + δ), where δ is the round-off error.

xn+1 = f l(x n+1 ) = f l xn 2 × f l 3 -f l a × f l(x 2 n ) xn+1 = xn 2 × 3 -a × x2 n (1 + ) 2 (1 + ) 2 (1) Let e(x n ) = 3 -a × x2 n (1 + ) 2 . |e(x n )| ≤ |3 -a × x2 n | + |2a × x2 n × | (2)
We have:

√ a × x n < 1 (3) Proof. x n+1 √ a -1 = -(1 + x n 1 2 √ a)(x n √ a -1) 2 < 0
Based on Equation 3we can conclude that

< (3 -a × x 2 n ) × 2 . (4)
Then, combining Equations 2 and 4 gives:

|e(x n )| < (3 -a × x 2 n ) × (1 + ), (5)
and combining Equations 1 and 5 gives

xn+1 < xn 2 × (3 -a × x2 n ) × (1 + ) 3 .
Disregarding O( n ) terms, n > 1, in the previous inequality [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF], we obtain

xn+1 < xn 2 × (3 -a × x2 n ) × (1 + 3 ), (6)
from which we can conclude that the round-off error is 3 per iteration. Note that xn+1 converges towards 1 √ a and that we want to bound the error of ŷn+1 , which converges towards √ a. So, only the last iteration should be multiplied by a:

x n+1 × a = y n+1 =⇒ f l(x n+1 × a) = f l(y n+1 ) =⇒ xn+1 × a × (1 + ) = ŷn+1 .
Consequently, the overall round-off error is δ ≤ n × 3 + .

B. Bounding the Systematic Error

The absolute systematic error in computing The relative systematic error in computing

1 √ a is |x i - 1 √ a |.
1 √ a is |x i -1 √ a | 1 √ a = | √ ax i -1|.
To find an upper bound for the systematic error, we start with a reasonable initial value defined as

1 2 × 1 √ a ≤ x 0 ≤ 3 2 × 1 √ a =⇒ 1 + 0.5 × √ a × x 0 ≤ 7 4 and √ a × x 0 -1 ≤ 1 2 . So, according to the proof of Equation 3: | √ a × x 1 -1| ≤ 4 7 × | 7 4 × ( √ a × x 0 -1)| 2 . Since 7 4 × ( √ a × x 0 -1) ≤ 7 8 then, | √ a × x 1 -1| ≤ 1 2 × 7 8 , It is safe to conclude that the systematic error in computing 1 √ a at iteration i = 0..n is | √ a × x i -1| ≤ 1 2 × ( 7 8
)

2 i -1 . (7)
The relative systematic error for 1 √ a is the same as the one for

√ a: |yi- √ a| √ a = |axi- √ a| √ a = | √ ax i -1|.
Both the systematic and round-off errors are plotted for different precision values (p) in Fig. 1. For space reason, only two of the resulting graphs are illustrated. The straight line corresponds to the round-off error, which increases proportionally to the number of iterations of the Newton Raphson method, whereas the curve designates the systematic error, which declines as the number of iterations grows.

The rise of the round-off error is remarkably rapid at lower precision values (e.g. Fig. 1(a)) but it slows down at higher precision (e.g. Fig. 1(b)).

The systematic error, on the other hand, is independent of the precision. It stabilizes at zero around the sixth iteration. This means that, disregarding the round-off error (round-off error equal to 0), six Newton Raphson iterations are sufficient in producing accurate (i.e. error=0) square root approximation. However, when the round-off error is in the picture, which is usually the case, the higher the number of iterations, the higher the round-off error is. Although at six sqrt iterations the systematic error is 0, the round-off error is 0.6 for p = 4 and 0.05 for p = 9.

The round-off error accumulates as the number of iterations rises, which counteracts the decline of the unavoidable error. That is why, the round-off error should not exceed the systematic error. Thus, the intersection between the two lines is indicative of the optimal number of iterations at which we have a good enough result.

IV. APPLICATION TO K-MEANS

So far, we studied the square root operation individually, in the context of error analysis for the purpose of precision tuning. It would be interesting to apply our study of the square root in a fullblown application that encompasses both smooth and non-smooth operations. We chose a clustering algorithm called K-means that uses the square root operation in the computation of the Euclidean distance function.

In this section, we start by describing K-means and its use in color quantization. Then, we present precision tuninig of K-means that involves both FLP simulation and our analytical results.

A. Color Quantization Using K-Means

The basic principle of the clustering problem is as follows. A data set X = (x 1 , x 2 , ..., x n ) is composed of data elements, where each data element x i is a d-dimensional vector. This data set is partitioned into k clusters of similar data points [START_REF] Jain | Data clustering: A review[END_REF].

In the particular case of color quantization, one of K-means' applications, an image represented as an array of size N × d where N is the number of pixels in the image and d is the color space (referred to as data dimension or data features), is represented with a smaller number of colors. Most commonly, the space is 3-dimensional (d = 3) and the coordinates encode the color. For example, RGB is a color space usually encoded as a 3-tuple of 8 bits each. The value of each dimension is within the range of [0, 255].

The flowchart in Fig. 2 gives an overview of the different steps of the clustering process in the context of color quantization with K-means.

The goal of the K-means clustering algorithm is to generate a compressed image out of the original image. To this end, a palette, i.e. a set of centroids C = (c 1 , c 2 , ..., c k ), is firstly chosen by selecting data elements (i.e. colors or pixels) that best represent the original image, for each cluster. In other words, the number of centroids is the number of colors that the palette is made of. Then, each data element is mapped to the closest [START_REF] Sinha | A perceptually based comparison of image similarity metrics[END_REF]. This results in a preliminary classification. Following this initial classification, k new centroids are re-computed as barycenters of the previously generated clusters and data elements are re-assigned to the new centroids based on the chosen distance measure. The algorithm is based on an iterative refinement technique and it iterates until centroids do not move anymore (i.e., no data elements change clusters), a maximum number of iterations is reached or the distance measure is minimized, as illustrated in Fig. 2. To quantify the difference between two pixels, a simpler distance function known as the Manhattan distance (sum of absolute differences) is usually used. However the optimal distance, in terms of minimizing within-class variance and producing higher-quality clusters, and more complicated, in terms of computational complexity, is the Euclidean distance (sum of the squared distances) [START_REF] Leeser | Design tradeoffs in a hardware implementation of the k-means clustering algorithm[END_REF]. Unlike other studies on approximate K-means, we use the optimal distance metric:

dist(p1, p2) = (R1 -R2) 2 + (G1 -G2) 2 + (B1 -B2) 2 .

B. Precision Tuning Using a Multiple-Precision Floating Point Library

We focus our precision analysis on the Euclidean distance function, which is the kernel of the K-means algorithm as it is in charge of computing the distance between a given data point and a cluster centroid. These distance values are pivotal in the correct assignment of pixels to their closest cluster and thus in providing a correct classification. We are going to include approximation in the FLP computations of the Euclidean distance function and leverage the results from our analysis in Section III to account for the impact of FLP precision variation on the square root operation.

Introducing inexactness by reducing the number of bits of the mantissa can sometimes provide results of the same accuracy but with better energy efficiency than the exact version of the program. To this aim, we studied the sensitivity of K-means by arbitrarily varying the number of bits of the mantissa (2-23 bits) and observing the repercussions on the output.

As was established in our analytical study of the square root operation (Section III), the number of bits of the mantissa has an impact on the number of square root iterations. So, to study the sensitivity of K-means, we re-wrote K-means' Euclidean distance function using the MPFR library, which is a smooth extension of the IEEE-754 standard where any FLP number can have its own precision [START_REF] Fousse | Mpfr: A multiple-precision binary floating-point library with correct rounding[END_REF], and we assigned to each precision its corresponding number of square root iterations according to the analysis in Section III. Our analytical results come in handy in the precision tuning of K-means in that, we do not have to vary the number of square root iterations for each precision simulation in order to determine the optimal pair (p, sqrt). The optimal number of square root iterations is statically determined for each precision (Section III). For example, for a precision p = 4, the FLP variables in K-means' Euclidean distance function are transformed into MPFR variables (mpf r t var), their precision is set to 4 (mpf r init2(var, 4)) and the number of iterations in the Newton Raphson method implementation is set directly to sqrt = 2.

In addition to the FLP variables, we tracked the operations that use these variables and changed them into MPFR operations. The Euclidean distance function encompasses multiplication, addition, subtraction and square root operations.

Listing 1: MPFR operations 1 mpfr sub ( r , r , r1 , MPFR RNDN ) ; / * r=r-r 1 ; * / 2 mpfr mul ( r , r , r , MPFR RNDN ) ; / * r=r * r ; * / 3 mpfr add ( r , r , r2 , MPFR RNDN ) ; / * r=r+r 2 ; * / 4 m p f r s q r t ( r tmp , r , MPFR RNDN ) ; / * r t m p= s q r t ( r ) * / Listing 1 showcases the transformation of different operations into MPFR operations. We implemented two inexact versions of the Euclidean distance function. We called these versions fused and unfused. In the first, the sub, mul and add operations are computed with full precision and the result is rounded once to N significant bits and then passed to the reduced-precision square root function. In the latter, all the operations are performed with reduced-precision, which means that the value is rounded four times before yielding the final result. We examine both fused and unfused versions in Section V to determine which version is more beneficial to implement.

The quality of the result as well as the energy consumption depend not only on the precision p of the FLP variables but also on other parameters such as the number of clusters k and the number of iterations n needed for the clustering process to converge. More clusters means that the compressed image will contain more colors, which is in favor of a better quality. With more iterations, there is a better chance for K-means to converge, which also contributes to the quality of the output. So, we investigated precision variation and its impact on the square root iterations over a range of different values of (k, n) pairs. Thus, we also varied the (k, n) parameters along with the precision and quantified the impact on the output in order to discern the best configuration that yields the desired energy-QoS tradeoff.

V. EXPERIMENTAL RESULTS

To determine whether a configuration (k, n) for a given precision p is favorable, we transform the original source code of K-means [START_REF] Yazdanbakhsh | Axbench : A benchmark suite for approximate computing across the system stack[END_REF] using the configuration in question, compile it to an ARM binary, execute it while using energy and instruction counters, and check for two criteria: energy consumption and QoS.

A. Experimental Setup

In order to quantify the impact of approximation on energy consumption, we make use of an energy model based on energy measurements of ARM Cortex-A7 instructions at a frequency of 500 MHz and without dependencies between instructions [START_REF] Vasilakis | An instruction level energy characterization of arm processors[END_REF]. We also make use of profile information, i.e., instruction types (e.g., add, mul, div) and the number of execution of each type. ARM's Cortex-A7 core is a dual issue in-order core that consists of one load/store, one multiply, one FLP and two integer units. A thorough characterization of the ARM Cortex-A7 instruction set with energy metrics for every instruction type could be found in [START_REF] Vasilakis | An instruction level energy characterization of arm processors[END_REF].

The energy values of integer and FLP instructions are normalized with reference to integer multiplication, i.e., integer multiplication is used as a unit of measure throughout the experiments. Normalized energy values of 32-bit integer and FLP instructions are listed in Table I. These normalized energy values are used to compute the energy consumption of exact K-means, which is considered as a touchstone. As for the energy values of variable precision (2-16 bits) operations, namely addition, multiplication and division, they are computed according to [START_REF] Barrois | Customizing Fixed-Point and Floating-Point Arithmetic -A Case Study in K-Means Clustering[END_REF] and shown in Fig. 3. It is noticeable from Fig. 3 that energy consumption is influenced by the number of bits of the mantissa; the higher the precision the higher the energy cost.

The overall energy consumption is computed by combining the profile information and the energy values

E total = #types i=1 op i × e i ,
where op i is the number of operations of type i and e i is the energy consumed per operation of type i (normalized to muli). We estimate the energy savings by comparing the energy consumption obtained by executing the exact version of K-means (i.e. with the highest precision p = 23) to the approximated version.

The quality of the result is measured using the SSIM index. SSIM is a perception-based metric that compares two images by incorporating a number of terms including contrast, luminance and structural information. The degradation of the image quality, due to compression with inexact K-means (reduced precision

B. Results

In this section, we study the impact of the precision and square root iterations (p, sqrt) on the quality of the output for various n (number of iterations needed for K-means to converge) and k (number of clusters) configurations. To do so, we start by jointly varying the precision bits p and the number of K-means iterations n, while fixing the number of clusters to k = 50, for both unfused and fused operations.

Then, we vary the number of clusters k and the precision p while fixing the number of K-means iterations to n = 10. The SSIM values are averaged over 10 different RGB images chosen from [START_REF]Uncompressed rgb images[END_REF]. An SSIM value is within the interval [0, 1], 1 designating the best image quality.

The first set of experiments that we conducted, i.e., SSIM(p, n) shows that the number of iterations n has little to no impact on SSIM for different precisions and for both the fused and unfused cases. This is true for other k values as well (10, 100, 500, 1000). Accordingly, the number of iterations n is fixed to 10 for the rest of the experiments.

However, by varying the precision with the number of clusters while fixing n to 10, a noticeable change in SSIM is observed and the results are reported in Fig. 4. For instance, in case of unfused operations (Fig. 4 SSIM values for the majority of the configurations. Taking into consideration the additional requirements that come with fused operations (i.e., dedicated hardware), which is not always supported by all architectures, as well as a compiler that allows the program to make use of the fused operations, we deem it too big an effort for an insignificant quality enhancement. Consequently, we advocate the use of unfused operations in the context of reduced-precision K-means.

It is also worth mentioning that we conducted the experiments from p = 2 to p = 23, but no noticeable change was detected starting from p = 16, as can be observed in Fig. 4. Thus, the rest of the graphs presented in this paper stop at p = 16.

To conclude whether approximating K-means is worthwhile, we also measured the energy gains. The percentages of energy savings are plotted against p for different (k, sqrt) configurations in Fig. 7. Each precision is attributed its optimal number of square root iterations according to our analytical study (Section III). For example, for p = 4, sqrt = 2. The graph indicates that there are energy gains for the different (p, k, sqrt) configurations but with different magnitudes. The energy gain of each configuration (p, k, sqrt) is computed with respect to the reference configuration (p = 23, k = 100, sqrt = 6). The smallest gain is observed at p = 16 and (k = 100, sqrt = 6) with a value of 14.09%. The biggest gain is at p = 2 and (k = 10, sqrt = 1) with a value of 96.57%. To be able to determine the best energy-QoS tradeoff Assuming that an output image with an SSIM value lying within the range of [0.95, 1] is considered a good quality image, then configuring K-means with a 6-bit mantissa, 100 clusters and 4 square root iterations is sufficient in producing such an image while saving 41.87% of energy.

VI. CONCLUSION

This paper focused on optimizing the implementation of the square root operation. An analytical examination of the Newton Raphson approximation of the square root showed that the number of square root iterations is influenced by the precision bits. Consequently, we associated to each precision its optimal number of Newton Raphson iterations. For a well-rounded evaluation of our proposition, we aimed at finding opportunities to reduce energy consumption of a classic clustering algorithm called K-means by varying the precision of its FLP variables and adjusting the number of square root iterations of its distance function.

The various approximated versions of K-means were compared to the exact version in terms of QoS, measured with SSIM, and relative energy gain. The obtained results can serve as a guideline in choosing the best configuration of precision bits, number of clusters and square root iterations (p, k, sqrt), i.e., one that yields the most energy gains, for a desired QoS (e.g., for an SSIM within [0.95, 1], an energy gain of 41.87% is achieved with a (6, 100, 4) configuration).

Fig. 1 :

 1 Fig. 1: Systematic and round-off errors for different precision values and iteration values of Newton Raphson

startFig. 2 :

 2 Fig. 2: Flowchart of the K-means algorithm

Fig. 3 :

 3 Fig. 3: Normalized energy values for different precisions

Fig. 4 :

 4 Fig. 4: Impact of precision variation and number of clusters on SSIM: for (a) unfused operations and (b) fused operations

Fig. 5 :Fig. 6 :

 56 Fig. 5: Executed instructions breakdown for exact K-means (p = 23) with different (k, sqrt) configurations

Fig. 7 :

 7 Fig. 7: Percentage of energy gain of different precisions p and with different (k, sqrt) configurations

TABLE I :

 I Energy values (normalized to muli) of a sample of integer and 32-bit FLP instructions

	Instr	addi	muli	si	li	addf mulf	sf	lf	divf
	Energy	1.02	1	2.41 1.79	1.16	1.16	2.34	1.92	7.80
	instri:integer instruction, instrf:FLP instruction