Stephane Redon
email: redon@cs.unc.edu

Ming C Lin
email: lin@cs.unc.edu

An Efficient, Error-Bounded Approximation Algorithm for Simulating Quasi-Statics of Complex Linkages

Keywords: Design and Analysis, Simulation, Virtual Prototyping

Design and analysis of articulated mechanical structures, commonly referred to as linkages, is an integral part of any CAD/CAM system. The most common approaches formulate the problem as purely geometric in nature, though dynamics or quasi-statics of linkages should also be considered. Existing optimal algorithms that compute forward dynamics or quasi-statics of linkages have a linear runtime dependence on the number of joints in the linkage. When forces are applied to a linkage, these techniques need to compute the accelerations of all the joints and can become impractical for rapid prototyping of highly complex linkages with a large number of joints.

We introduce a novel algorithm that enables adaptive refinement of the forward quasi-statics simulation of complex linkages. This algorithm can cull away joints whose contribution to the overall linkage motion is below a given user-defined threshold, thus limiting the computation of the joint accelerations and forces to those that contribute most to the overall motion. It also allows a natural trade-off between the precision of the resulting simulation and the time required to compute it. We have implemented our algorithm and tested its performance on complex benchmarks consisting of up to 50,000 joints. We demonstrate that in some cases our algorithm is able to achieve up to two orders of magnitude of performance improvement, while providing a high-precision, error-bounded approximation of the quasi-statics of the simulated linkage.

Introduction

Modeling and design of articulated mechanical structures, commonly referred to as "linkages", is an integral part of any CAD/CAM system. The existing approaches typically focus on sizing mechanical constraints that guide the movement of the endeffector of the system within a workspace. The function of the device is usually prescribed as a set of positions reachable by the endeffector, as well as the mechanical constraints formed by joints that limit relative movement. The objective is to compute all the linkage configurations that can achieve a specific task. Formulated in this way the design problem is purely geometric in nature, often with an emphasis on systems consisting of lower numbers of degrees of freedom to reduce the complexity of the problem.

Complex linkages, such as articulated systems of many links, molecular chains, or articulated-body representations of deformable bodies, often require consideration of additional physical constraints and present new computational challenges. Forward simulation of articulated mechanical systems, i.e. computing the joint accelerations and forces depending on external forces and applied torques, is a fundamental issue that should be taken into account in designing highly complex linkages. Especially, a quasi-static analysis can provide insights into mechanisms at or close to equilibrium, or which have small dynamic effects (slow evolution or small-scale motions). However, such computations are frequently expensive to perform and sometimes ignored in rapid prototyping of highly complex mechanical linkages. Existing optimal algorithms that compute forward dynamics or quasi-statics of linkages have a linear runtime dependence on the number of joints in the linkage. When forces are applied to the linkage, these techniques need to compute the accelerations of all the joints, which can become impractical for rapid prototyping of complex linkages with a large number of joints.

In this paper, we present a novel algorithm for progressive refinement of the quasi-static simulation of complex linkages, which offers three key advantages:

• Error-bounded Approximation: Given any user-defined error threshold, our algorithm is able to determine an approximation of the linkage acceleration within the required error bound.

• Progressive Refinement: Our algorithm operates by progressively increasing the quality of the approximate linkage acceleration, thus providing a trade-off between the precision of the approximation and the time required to compute it.

• Efficiency: Our algorithm determines the approximate acceleration by computing only some of the joints' accelerations and associated coefficients, resulting in a significant performance improvement when the ratio on the number of processed joints to the total number of joints is small, i.e. when the required precision is low, or when the forces applied to the linkage have a local effect only on the linkage motion.

To the best of our knowledge, this algorithm is the first technique that offers a method to progressively refine the computation of the quasi-statics of a linkage while guaranteeing the precision of the calculations. We have implemented our algorithm in C++ and tested it on some complex mechanical systems consisting of up to 50,000 links. In some cases, we were able to achieve up to two orders of magnitude of performance improvement over an existing linear-time algorithm, while providing a high-precision, errorbounded approximation of the quasi-statics of the simulated linkage (see an example in Figure 1).

The rest of the paper is organized as follows. Section 2 presents a brief survey of related work on rigid body systems. Section 3 gives an overview of the basic mathematical notation and representations that are fundamental to our approach. Section 4 demonstrates how we obtain an approximate acceleration of a linkage. Section 5 shows how, using the hierarchical reference frames representation we introduce, the linkage state update can be restricted to some coefficients only. Section 6 presents some results obtained with our new approach, and Section 7 concludes by suggesting some potential research directions.

Related Work

The computation of the kinematics, quasi-statics, or dynamics of articulated bodies has been a subject of study for many years (cf [START_REF] Hooker | The dynamical attitude equations for an n-body satellite[END_REF] for example). We only provide a brief overview here and refer the readers to [START_REF] Featherstone | Robot dynamics: Equations and algorithms[END_REF] for a recent state-of-the-art survey.

Essentially, an important step in the development of simulation algorithms has been the introduction of methods which have a linear complexity in the number of bodies in the multibody system. Some of these linear methods [START_REF] Bae | A recursive formulation for constrained mechanical systems dynamics: Part i. open-loop systems[END_REF][START_REF] Brandl | A very efficient algorithm for the simulation of robots and similar multibody systems without inversion of the mass matrix[END_REF][START_REF] Featherstone | Robot Dynamics Algorithms[END_REF][START_REF] Hollerbach | A recursive lagrangian formulation of manipulator dynamics and a comparative study of dynamics formulation complexity[END_REF][START_REF] Mcmillan | Efficient computation of articulated-body inertias using successive axial screws[END_REF], which rely on a recursive formulation of the motion equations, have now been known for more than two decades. In parallel, several authors have contributed to simplify the motion equations by developing new notations and formulations. Some of these include the spatial notation [START_REF] Featherstone | Robot Dynamics Algorithms[END_REF]], and its new version [Featherstone 1999a;Featherstone 1999b], the spatial operator algebra [START_REF] Rodriguez | Spatial operator algebra for manipulator modelling and control[END_REF], as well as Lie-Group formulations [START_REF] Mueller | A lie-group formulation of kinematics and dynamics of constrained mbs and its application to analytical mechanics[END_REF]. More recently, several researchers have proposed parallel algorithms to compute the forward dynamics of articulated bodies using multiple processors [START_REF] Fijany | Parallel o(log n) algorithms for computation of manipulator forward dynamics[END_REF]Featherstone 1999a;Featherstone 1999b].

To our knowledge, no algorithm has been proposed for obtaining an error-bounded approximation of the quasi-statics of a multibody system. All existing techniques compute the accelerations of all joints in the multibody system. They cannot detect nor take advantage of the situations where the interaction between the environment and the linkage is localized, i.e. the situations where computing only some of the joint accelerations is enough to obtain a useful, error-bounded approximation of the linkage motion.

Preliminaries

In this section, we first provide an overview of the divide-andconquer algorithm (DCA) introduced by Featherstone [Featherstone 1999a;Featherstone 1999b], upon which our algorithm is built. We restrict the presentation to the quasi-statics case, where the body velocities are zero at each instant. Then, we give an overview of our approach.

Overview of the DCA

Consider a system of rigid bodies connected by joints. Assume the initial body positions and the acceleration-independent forces are known. The acceleration-independent forces may include gravity, applied external forces, and active joint forces produced by actuators. The forward quasi-statics problem is to determine the body accelerations and the kinematic constraint forces based on the initial body positions and the acceleration-independent forces, assuming the velocities of the rigid bodies are zero at each instant. To solve this problem (and the more general forward dynamics problem, where the body velocities are non zero), Featherstone [START_REF] Featherstone | Robot Dynamics Algorithms[END_REF]] defines an articulated body as a system of rigid bodies with a set of handles, i.e. locations attached to some of the rigid bodies where external forces may be applied. Figure 2(a) shows one rigid body with three handles which, although abstract in nature, have been materialized for clarity. Note that the number of handles m in the articulated body is independent of the number of bodies or joints in the articulated body. That is, the handles are a set of arbitrary locations on the articulated body that are used to simplify the equations of motion by focusing on some of the rigid bodies only. More specifically, [START_REF] Featherstone | Robot Dynamics Algorithms[END_REF] shows that the handles accelerations can be expressed as affine functions of the forces applied to the handles according to the following articulated-body equation:

    a 1 a 2 . . . a m     =     Φ 1 Φ 12 • • • Φ 1m Φ 21 Φ 2 • • • Φ 2m . . . . . . . . . . . . Φ m1 Φ m2 • • • Φ m         f 1 f 2 . . . f m     +     b 1 b 2 . . . b m     (1)
where a i is the 6 × 1 spatial acceleration of handle i, f i is the 6 × 1 spatial force applied to handle i, b i is the 6 × 1 bias acceleration of handle i (the acceleration handle i would have if all the external forces were zero), Φ i is the 6 × 6 inverse articulated-body inertia of handle i, and Φ i j is the 6 × 6 cross-coupling inverse inertia between handles i and j.

In the DCA [Featherstone 1999a;Featherstone 1999b], an articulated body is recursively defined by assembling two (floating) articulated bodies together. Thus, the complete articulated body is described by a binary assembly tree, which specifies how articulated bodies (subassemblies) are repeatedly connected using joints. The leaf nodes in the assembly tree are rigid bodies, and the root node is a special node with only one child, which describes how the complete articulated body relates to the world coordinate system, potentially connecting the articulated body to a fixed base. The root node is the only node where an articulated body can be attached to a fixed base.

The DCA involves four passes on the assembly tree, i.e. four com-plete traversals of the assembly tree. Two preliminary passes compute the new body positions from the current values of the joint positions and accelerations, and the coordinate transformation matrices needed by the algorithm, while the next two passes perform the actual dynamics computations:

• Main pass: The main pass of the DCA determines the coefficients in the articulated-body equations (1) for each node (i.e. articulated body) in the assembly tree, from the leaf nodes to the root node.

• Back-substitution pass: Once the articulated-body equations are computed, the back-substitution pass computes the joint accelerations and the kinematic constraint forces, from the root node to the leaf nodes. The articulated-body equation of A is:

Main pass

     a A 1 1 . . . a A 1m a A 2      =      Φ A 1 11 • • • Φ A 1 1m Φ A 1 1 2 . . . . . . . . . . . . Φ A 1 m1 • • • Φ A 1mm Φ A 1m2 Φ A 21 1 • • • Φ A 21m Φ A 2           f A 1 1 . . . f A 1m f A 2      +      b A 1 1 . . . b A 1m b A 2      .
This can be rewritten such as to single out the handle

H A 2 : a A 1 a A 2 = Φ A 1 Φ A 12 Φ A 21 Φ A 2 f A 1 f A 2 + b A 1 b A 2 , where Φ A 1 = mat(Φ A 1 i j
) is a 6m × 6m composite matrix of crosscoupling and inverse articulated-body inertias involving handles

H A 1 1 , . . . , H A 1m , Φ A 12 = (Φ A 21 ) T = col(Φ A 1 i 2 )
a B 1 a B 2 = Φ B 1 Φ B 12 Φ B 21 Φ B 2 f B 1 f B 2 + b B 1 b B 2 ,
where handle H B 1 is singled out, and the quantities involving the other handles (and potentially the handle H B 1 ) are composite vectors or matrices.

Finally, the articulated-body equation of C can be written:

a A 1 a B 2 = Φ C 1 Φ C 12 Φ C 21 Φ C 2 f A 1 f B 2 + b C 1 b C 2 ,
where all quantities are composite vectors or matrices involving the remaining m + n handles H A 1 1 , . . . , H A 1m and H B 2 1 , . . . , H B 2n . Note that the handles H A 2 and H B 1 are consumed by the assembly operation, and the forces f A 2 and f B 1 , that were external forces with respect to the articulated bodies A and B, are now internal kinematic constraint forces for the articulated body C.

Assume S denotes the 6 × d J joint's motion subspace, where d J is the number of degrees of freedom of the joint. A typical value for a revolute joint in link coordinates is S = (0, 0, 1, 0, 0, 0) T . Assume Q denotes the d J × 1 vector of active joint forces. Setting

V = (Φ A 2 + Φ B 1 ) -1 W = V -VS(S T VS) -1 S T V β = b A 2 -b B 1 γ = Wβ + VS(S T VS) -1 Q,
it can be shown (cf [Featherstone 1999a;Featherstone 1999b]) that the articulated-body equation of C can be obtained from those of A and B:

Φ C 1 = Φ A 1 -Φ A 12 WΦ A 21 Φ C 2 = Φ B 2 -Φ B 21 WΦ B 12 Φ C 21 = Φ B 21 WΦ A 21 = (Φ C 12 ) T b C 1 = b A 1 -Φ A 12 γ b C 2 = b B 2 + Φ B 21 γ (2)
For leaf nodes, which represent rigid bodies, the articulated-body equation coefficients in the case of quasi-statics are:

Φ i = Φ i j = I -1 b i = I -1 f k , ( 3
)
where I is the rigid body spatial inertia, and f k is an accelerationindependent external force applied to the rigid body (any number of handles can be attached to a rigid body). Since we assume the evolution of the rigid-body system is quasi-static, the velocitydependent term in the bias acceleration is zero.

Back-substitution pass

Once the coefficients of the articulated-body equations have been updated in the main pass, the back-substitution pass determines the joint accelerations and the kinematic constraint forces, from the root node to the leaf nodes.

Based on the values of the forces f A 1 and f B 2 computed when processing the parent nodes, the d J × 1 joint acceleration vector is:

q = (S T VS) -1 (Q -S T V(Φ A 21 f A 1 -Φ B 12 f B 2 + β)), (4)
and the 6 × 1 kinematic constraint forces are:

f B 1 = WΦ A 21 f A 1 -WΦ B 12 f B 2 + γ f A 2 = -f B 1 (5)
The root node calculations take two possible forms, depending on the existence of a connection between the articulated body and a fixed base. Since the root node describes the placement of the complete articulated body in the world frame, it has only one articulated-body child B. The external force f B 2 is assumed to be known (usually f B 2 = 0, since all acceleration-independent forces have already been accounted for in the bias accelerations).

If the articulated body has a floating base, then

f B 1 = 0 and a B 1 = Φ B 12 f B 2 + b B 1 .
However, if the articulated body is attached to a fixed base through a joint, the joint acceleration is:

q = (S T (Φ B 1 ) -1 S) -1 (Q + S T (Φ B 1 ) -1 (Φ B 12 f B 2 + b B 1 )), (6)
and the kinematic constraint force is:

f B 1 = (Φ B 1 ) -1 (Sq -Φ B 12 f B 2 -b B 1 ). (7)
Again, the velocity-dependent terms have been removed. Note that it is possible to simulate gravity by substituting b B 1 -a g for b B 1 , which suppresses the need to account for gravity through the action of an external force upon all leaf nodes.

Overview of Our Algorithm

Assume the initial acceleration-independent forces, the initial body positions and velocities, the initial coordinates transformation matrices, and the initial articulated-body equations are known. We can view the DCA as a two-step process which is iterated to simulate the quasi-statics or dynamics of a linkage:

• Joint accelerations and forces update: The joint accelerations and kinematic constraint forces are computed based on the current articulated-body equations.

• Linkage state update: The body positions (and velocities in the dynamics case), the coordinates transformation matrices, and the articulated-body equations are updated based on the joint accelerations computed in the first step and the new acceleration-independent forces.

When run on a single processor, the DCA has a linear complexity in the number of joints in the complete articulated body, since each pass in the algorithm needs to traverse the complete assembly tree. There are no mechanisms in the DCA which allow to predict the joints with largest accelerations, i.e. to determine them without having to compute all joint accelerations.

In order to avoid this linear complexity, our algorithm extends the original DCA by modifying each of its two steps in the following way:

• Approximate linkage acceleration: We introduce a class of approximate linkage accelerations. Given any user-defined error threshold, we are able to determine an approximation of the linkage acceleration, while guaranteeing that the error caused by the approximation is below the given user-defined threshold. We introduce and calculate some total acceleration equations, which help us order the computations in the back-substitution pass and cull away some joints whose contribution to the linkage acceleration is found to be insignificant. This allows us to compute an approximate acceleration which is provably error-bounded, without having to complete the back-substitution pass and compute all the joint accelerations.

• Simplified linkage state update: We introduce a hierarchical reference frames representation in which we express the state of the linkage and the total acceleration equations. We show that, using this representation, only some of the quantities involved in the second step have to be updated, thereby avoiding a linear-time update.

Acceleration Approximation

In this section, we introduce a method which allows us to adaptively refine the computation of the joint accelerations and (kinematic constraint) forces, and obtain an approximation of the linkage acceleration which is provably error-bounded.

The key to our approach is the demonstration that it is possible to compute the total acceleration of an articulated body, i.e. the sum of the (squared) norms of the joint accelerations of all the joints in the articulated body, without computing each joint acceleration separately. We use the total acceleration during the back-substitution pass of the DCA to progressively refine the computation of the quasi-statics of the linkage, and stop the refinement when a sufficient precision has been obtained.

As can be seen in Section 3.1, the DCA relies on a bottom-up pass (the main pass), which calculates the articulated-body equations of all the nodes in the assembly tree, to subsequently compute the joint accelerations and forces during a top-down pass (the backsubstitution pass). We demonstrate that, similarly, it is possible to calculate some total acceleration equations in a bottom-up pass, which allow us to compute the total accelerations in a top-down fashion.

Definitions

Let C be a linkage, and let q C denote the composite vector of the positions of the N C joints in C:

q C =    q 1 . . . q N C    .
Each joint position q i is a d J i × 1 vector, where d J i is the number of degrees of freedom for joint i in C, and q C is a d J C × 1 vector, where

d J C = ∑ N C
i=1 d J i is the total number of degrees of freedom in C. Similarly, the composite vector qC contains the accelerations of the N C joints in C.

In this paper, we propose to define an approximation of the linkage acceleration as a partial computation of the joint accelerations. Specifically, we define an approximate acceleration in the following way.

Definition 1 (Approximate Acceleration). Consider an articulated body C and its actual acceleration qC . An approximate acceleration s C of the exact linkage acceleration qC is a composite vector in which some block components are equal to the actual joint accelerations in C, and the other block components are zero. All blocks in s C have the same dimension as their counterpart in qC , and there exists a partition of {1, . . . , N C } into two subsets E and Z such as:

s C =    s 1 . . . s N C    and s C i = qi , i ∈ E s C j = 0, j ∈ Z .
The subset E is the set of indices of the joints whose acceleration has been computed exactly, while the subset Z is the set of indices of the joints whose acceleration has been assumed to be zero, resulting in a potential approximation of the actual linkage acceleration qC .

In order to measure the quality of an approximate acceleration, we use several classical error measures.

Definition 2 (Acceleration Error Measures). Consider an articulated body C, its actual acceleration qC , and an approximation s C of its acceleration. The acceleration error ε(C) will alternatively denote:

• The absolute linkage acceleration error:

ε(C) = ||s C -qC || 2 = ∑ j∈Z qT j qj
• The relative linkage acceleration error:

ε(C) = ||s C -qC || 2 ||q C || 2 = ∑ j∈Z qT j qj ∑ N C j=1 qT j qj
• The absolute joint acceleration error:

ε(C) = max j∈Z ||q j || 2
• The relative joint acceleration error:

ε(C) = max j∈Z ||q j || 2 max j=1...N C ||q j || 2
In the rest of this section we show that, given any user-defined maximum acceleration error ε max , and whichever error measure ε(C) we choose, we are able to find an approximate acceleration s C which is sufficiently close to the actual acceleration of the linkage, i.e. such as ε(C) ε max , without having to compute each joint acceleration in C.

We obtain an approximate linkage acceleration by performing a partial back-substitution pass, thereby computing the exact value of some joint accelerations, and implicitly assuming that the other joint accelerations are zero. We show that it is possible to not complete the back-substitution pass and still determine an approximate linkage acceleration within the required error bound.

The key to our approach is the use of the total acceleration of a linkage, defined as follows.

Definition 3 (Total Acceleration). Consider an articulated body C with actual acceleration qC . The total acceleration A (C) of C is the sum of the squared L 2 -norms of the accelerations of the N C joints in C:

A (C) = ||q C || 2 2 = N C ∑ i=1 qT i qi .
Alternatively, we will use an index-free notation:

A (C) = ∑ C qT q,
where a summation over a linkage C is a shorthand for the summation over all the joints in C.

Total Acceleration Equations

We demonstrate that the total acceleration of an articulated body is a quadratic function of the external forces acting upon it, as stated in the following lemma.

Lemma 1 (Total Acceleration Equation). The total acceleration A (C) of an articulated body C, formed by assembling an articulated body A containing m + 1 handles and an articulated body B containing n + 1 handles, is a quadratic function of the external forces acting upon C:

A (C) = f A 1 f B 2 T Ψ C 1 Ψ C 12 Ψ C 21 Ψ C 2 f A 1 f B 2 + f A 1 f B 2 T p C 1 p C 2 + η C , (8)
where

Ψ C 1 is a 6m × 6m matrix, Ψ C 2 is a 6n × 6n matrix, Ψ C 12 = (Ψ C 21 ) T is a 6m × 6n matrix, p C 1 is a 6m × 1 vector, p C
2 is a 6n × 1 vector, and η C is in IR. We call equation ( 8) the total acceleration equation of C.

Proof. We demonstrate the lemma by induction on the height of the subassembly tree h. For subassembly tree C of height zero representing a single (floating) rigid body, the total acceleration is zero, since there is no joint in a rigid body, and can thus be expressed as a quadratic function with

Ψ C 1 = 0, Ψ C 2 = 0, Ψ C 21 = (Ψ C 12 ) T = 0, p C 1 = 0, p C 2 = 0
, and η C = 0. Now assume that the result holds for any subassembly tree of height strictly smaller than h, and consider the assembly tree of an articulated body C with height h. Assume that C is formed by assembling two articulated bodies A and B. We will treat the case of the root node, with only one child B, afterwards.

The heights of the assembly trees of A and B are strictly smaller than h, and the total accelerations of A and B can thus be expressed as quadratic functions of the forces respectively acting upon them. Our goal is to show that the total acceleration of C is a quadratic function of the forces acting upon C.

The total acceleration A (C) of the articulated body C can be computed from those of A and B, respectively denoted by A (A) and A (B), and the acceleration q• of the joint connecting A and B:

A (C) = ∑ C qT q = qT • q• + ∑ A qT q + ∑ B qT q = qT • q• + A (A) + A (B).
Hence, it suffices to express each of the three terms in the last sum as a quadratic function of the forces acting upon C to prove that the result holds for C as well.

By the inductive hypothesis, the total acceleration A (A) can be expressed as a quadratic function of the forces acting upon A:

A (A) = f A 1 f A 2 T Ψ A 1 Ψ A 12 Ψ A 21 Ψ A 2 f A 1 f A 2 + f A 1 f A 2 T p A 1 p A 2 + η A .
However, equation ( 5) shows that the force f A 2 is an affine function of f A 1 and f B 2 , and thus all forces acting upon A are affine functions of the forces acting upon C:

f A 1 f A 2 = I d 0 -WΦ A 21 WΦ B 12 f A 1 f B 2 + 0 -γ ,
where I d is the identity matrix. Consequently, A (A) can be expressed as a quadratic function of the forces acting upon C.

Similarly, by the inductive hypothesis, the total acceleration A (B) can be expressed as a quadratic function of the forces acting upon B:

A (B) = f B 1 f B 2 T Ψ B 1 Ψ B 12 Ψ B 21 Ψ B 2 f B 1 f B 2 + f B 1 f B 2 T p B 1 p B 2 + η B But, since f B 1 = -f A 2 , the force f B
1 is an affine function of f A 1 and f B 2 , and thus the forces acting upon B are affine functions of the forces acting upon C:

f B 1 f B 2 = WΦ A 21 -WΦ B 12 0 Id f A 1 f B 2 + γ 0 ,
which shows that the total acceleration A (B) is a quadratic function of the forces acting upon C.

All that remains is to express qT • q• as a quadratic function of the forces acting upon C. In the case of a generic internal node, q• is computed using equation (4). To simplify the expression of qT

• q• in this case, let us set

U = (S T VS) -1 S T V N = U T U R = (S T VS) -1 (Q -S T Vβ).
We then have:

qT • q• = f A 1 f B 2 T Φ A 12 NΦ A 21 -Φ A 12 NΦ B 12 -Φ B 21 NΦ A 21 Φ B 21 NΦ B 12 f A 1 f B 2 +2 f A 1 f B 2 T -Φ A 12 U T R Φ B 21 U T R + R T R,
which is a quadratic function of the forces acting upon C.

The root node is a special case for which A is non-existent, and a supplementary joint acceleration is potentially defined when the articulated body is attached to a fixed base.

If the articulated body is floating, then q• = 0 and f

B 1 = 0. Thus A (C) = A (B) is a quadratic function of f B
2 only, i.e. the forces acting upon C.

If the articulated body is attached to a fixed base, however, q• is computed using equation ( 6), and f B 1 is computed using equation ( 7). Both quantities are affine functions of f B 2 , and thus A (C) = qT

• q• + A (B) is again a quadratic function of f B 2 , i.e. the forces acting upon C.

It should be clear from the proof that the coefficients in A (C) can be computed from the coefficients in A (A) and A (B). Setting

Y = N + W(Ψ A 2 + Ψ B 1 )W Z = 2U T R + W(p A 2 -p B 1 ) -2W(Ψ A 2 + Ψ B 1 )γ,
the coefficients of the total acceleration of a generic internal node can be written:

Ψ C 1 = Ψ A 1 + Φ A 12 YΦ A 21 -(Φ A 12 WΨ A 21 + Ψ A 12 WΦ A 21 ) Ψ C 2 = Ψ B 2 + Φ B 21 YΦ B 12 -(Φ B 21 WΨ B 12 + Ψ B 21 WΦ B 12 ) Ψ C 21 = -Φ B 21 YΦ A 21 + Φ B 21 WΨ A 21 + Ψ B 21 WΦ A 21 = (Ψ C 12 ) T p C 1 = p A 1 -Φ A 12 Z -2Ψ A 12 γ p C 2 = p B 2 + Φ B 21 Z + 2Ψ B 21 γ η C = η A + η B + R T R + γ T (Ψ A 2 + Ψ B 1 )γ -γ T (p A 2 -p B 1 ) (9)
The details of the calculations are omitted here. For the root node, the coefficients in A (C) are obtained by zeroing the coefficients relative to A in equation ( 9).

Consequently, the total acceleration equations can be computed in a bottom-up pass: the pass starts from the leaf nodes with all coefficients set to zero, and recursively applies the formulas (9) until the root of the assembly tree is reached. Assuming these total acceleration equations have been calculated, it is now possible to compute the total acceleration of a subassembly during the backsubstitution pass, without having to traverse all the descendents of the subassembly to determine the joint acceleration and forces of all the joints. Therefore, when a node C is processed during the back-substitution pass of the DCA, the forces acting upon this node are available, and equation ( 8) can be used to obtain A (C) without processing C's descendents.

In summary, the forces acting upon an articulated body can be used during the back-substitution pass to compute not only the principal joint's acceleration and forces, but also the total acceleration of the articulated body, before (or without) processing its descendents.

Approximate linkage acceleration

Consider an articulated body C, and let ε max 0 denote a given user-defined error threshold. Assume for now that the error measure we use is the absolute linkage acceleration error. The goal is to determine an approximate acceleration s C for which ε(C) ε max , by performing a partial back-substitution pass: some joint accelerations will be computed exactly, and some others will be implicitly set to zero.

Assume the state of the linkage and the acceleration-independent forces are known. Assume the articulated-body equations and the total acceleration equations are all up-to-date, and we are about to perform the back-substitution pass of the DCA. Since, when the back-substitution pass begins, the forces f A 1 and f B 2 acting on the complete articulated body C are assumed to be known, the total acceleration A (C) = ||q C || 2 2 can be obtained from equation (8). Now suppose A (C) = 0. This assumption implies that all the joint accelerations of all the joints in the articulated body are zero, since A (C) is the sum of the (positive) squared norms of all the joint accelerations. In this case, we have determined all the joint accelerations without descending the assembly tree, and we can stop Figure 4: Partial computation of the joint accelerations. Our algorithm computes only some joint accelerations in the articulated body, while guaranteeing that the error that results from this partial computation is smaller than a user-defined threshold.

the back-substitution pass at the root node. In this case, the approximate acceleration s C is actually equal to the exact one qC , and ε(C) = 0. Suppose now that A (C) is not zero, but still smaller than (ε max ) 2 . This means that all joint accelerations are sufficiently close to zero. If we stop the back-substitution pass right at the root node and set s C = 0, we still obtain an approximate acceleration which is close enough to the exact one:

ε(C) = A (C) ε max .
Assume finally that A (C) is strictly larger than (ε max ) 2 . We do not stop the traversal of the assembly tree right at the root node. We compute the joint acceleration and forces of C's principal joint, and recursively traverse the descendents of C using a priority queue and maintaining the value of the current acceleration error ε c , i.e. the error that would result from stopping the recursion.

The priority queue is used to order the traversal of the subassemblies. A node A has a higher priority than a node B if A (A) > A (B). During the traversal, we thus descend first in the portions of the linkages which have the highest total acceleration. We simultaneously maintain the current acceleration error ε c . At the beginning of the back-substitution pass, we set ε c = A (C), where C is the root node of the assembly tree, since this is the acceleration error that would result from stopping the back-substitution before having processed any node. During the back-substitution pass, whenever the acceleration q• of the principal joint of a node C is computed, we subtract its contribution qT

• q• to the current acceleration error:

ε c ← (ε c ) 2 -qT • q• .
Then, if C has children A and B, we compute their total acceleration A (A) and A (B) and push them into the priority queue. Since qT

• q• = A (C) -A (A) -A (B), the current acceleration error ε c is always equal to (the square root of) the sum of the total accelerations of the nodes contained in the priority queue:

ε c = ∑ C∈q A (C).

Since the norm qT

• q• of the principal joint acceleration of any node C is positive, the current acceleration error can only decrease as more nodes are processed. This allows us to stop the backsubstitution pass as soon as the current acceleration error becomes The case of the relative acceleration error measure can be treated similarly. Since we are able to compute the current acceleration error and, in the beginning of the back-substitution pass, the norm ||q C || 2 = A (C) of the linkage acceleration, we are able to control the relative linkage acceleration error as well.

Finally, since max C ||q|| 2 A (C), we can also control the absolute and relative joint acceleration errors.

In summary, we have shown that, provided the linkage state and the total acceleration equations are up-to-date, we are able to determine an approximate acceleration s C up to any level of precision. The time required to compute the approximate acceleration is O(N E ), where N E is the number of computed joint accelerations, that is, the size of the traversed sub-tree of the assembly tree (cf Figure 4). This algorithm can be simply written. Let q denote the priority queue. Assume that popping q returns the articulated body with the highest priority, and assume C denotes the node currently processed by the back-substitution pass. Let A and B denote the potential child nodes of C. The algorithm is:

ALGORITHM (ACCELERATION APPROXIMATION)
C ← root node of the assembly tree

ε c ← A (C) q ← push(C) while (ε c > ε max ) { C ← pop(q) compute joint acceleration q• and forces of C's principal joint ε c ← (ε c ) 2 -qT • q• for any articulated child A of C { compute A (A) q ← push(A) } }
When this algorithm terminates (which is guaranteed, provided the user-defined threshold ε max is greater than or equal to zero), an approximate linkage acceleration has been determined such that the total acceleration error is smaller than ε max .

Simplified Linkage State Update

As the linkage configuration and the forces applied to a linkage evolve, the coefficients describing the linkage state and the total acceleration equations need to be updated before performing the next back-substitution pass. In this section, we introduce a hierarchical reference frames representation in which we express the state of the linkage and the total acceleration equations. We show that, using this representation, only some of the quantities involved in the second step need to be updated, thereby avoiding a linear-time update of the equation coefficients. Our hierarchical reference frames representation resembles the hierarchy of transforms in the Chain-Tree data structure introduced by Lotan et al. [START_REF] Lotan | Efficient maintenance and self-collision testing for kinematic chains[END_REF] to describe the kinematics of protein chains. However, our representation is able to handle branches in the linkages.

Hierarchical reference frames

Consider a linkage C formed by assembling two articulated bodies A and B. We call the handles used to connect A and B the principal handles of A and B. The root node has no principal handle, since it is not connected to any other assembly.

As in Featherstone [Featherstone 1999a], each handle has a reference frame in which all the quantities related to that handle are expressed. During the simulation, we do not require maintaining the positions of the bodies in the world frame, but only three types of 4 × 4 homogeneous coordinates transformation matrices that we define as follows:

• Principal joint transforms: consider an articulated body C with a left child node A and a right child node B. By convention, we call principal joint transform the coordinates transformation matrix T C

• from the reference frame of B's principal handle to the reference frame of A's principal handle. This transform is a function of the position of C's principal joint only. The principal joint transform of the root node is the coordinates transformation matrix from the reference frame of B's principal handle to the world frame.

• Handles cross-transforms: consider an articulated body C

with handles H C 1 , . . . , H C h , where h is the number of handles of C. We maintain the coordinates transformation matrices T C i j from the reference frame of each handle H C j to the reference frame of other handles H C i , for 1 i, j h.

• Principal handles transforms: in each node, we maintain the coordinates transformation matrix from the reference frame of the node's principal handle to the reference frame of its parent node's principal handle. For example, if B is a node with parent C, we maintain the coordinates transformation matrix T C B from the reference frame of B's principal handle to the reference frame of C's principal handle. The root node has no such transform, and the principal handle transform of the only child of the root node is equal to the root node's principal joint transform, i.e. a transform from the reference frame of B's principal handle to the world frame. These 4 × 4 homogeneous coordinate transformation matrices are converted at will into 6 × 6 spatial motion or spatial force transformation matrices when necessary [Featherstone 1999a].

Note that these three types of coordinates transforms are the only ones we need to maintain. Especially, the position of a rigid body is never expressed in the world frame, except when an external force is applied to this rigid body. We note that these world transforms are not necessary in order to display the linkage, when a graphics library which adopts other hierarchical representations is used (eg. OpenGL).

These transforms can be computed recursively from the leaf nodes to the root node of the assembly tree. Figure 3 shows these transforms in the case of the articulated body assembled in Figure 2. Each handle has a reference frame rigidly attached to it, and the transform associated to a pair of handles are represented by arrows. We detail the computations in the rest of this section.

Handle's cross-transforms

The handles cross-transforms are computed from the leaf nodes to the root node of the assembly tree. If C is a leaf node of the assembly tree, its cross-transforms are constant and do not have to be ever updated, since the handles are fixed in the rigid body's reference frame. If C has some children A and B, however, its crosstransforms are computed from those of A and B, and the position of C's principal joint. Since the handles of C are either handles of A or B, we only need to compute the cross-transforms which involve one handle in A and one handle in B. Indeed, the cross-transforms from handles in A to handles in A have already been computed in A, as have the cross-transforms from handles in B to handles in B already been computed in B. The cross-transforms T C i j from a handle H C j in B to a handle H C i in A are computed using the principal joint transform T C

• . Let H A l denote A's principal handle, H B u denote B's principal handle, and assume

H C i = H A k and H C j = H B v .
Then:

T C i j = T A kl T C • T B uv , (10)
where T A kl is the cross-transform from H A l to H A k , and T B uv is the cross-transform from H B v to H B u . The cross-transform T C ji is computed by inverting T C i j .

Principal handles transforms

Similarly, the principal handles transforms are computed from the leaf nodes to the child nodes of the root node of the assembly tree. Assume C is a linkage with children A and B. Without loss of generality, assume we want to compute the transform T C B from the principal handle of B to the principal handle of C. If the principal handle of C is a handle of B, then T C

B is equal to one of B's crosstransforms. However, if the principal handle of C is a handle of A, we need to compute T C B using the principal joint transform T C • . Assume H C i = H A k is the principal handle of C, and let H A l denote the principal handle of A. Then:

T C B = T A kl T C • . (11)
The principal handle transform T C A can be computed in a similar way.

Simplified Coefficients Update

During the simulation we need to maintain, for each node in the assembly tree, the node's transforms, the node's inverse inertias and bias accelerations, and the node's total acceleration coefficients.

For any node C, most coefficients only depend on their counterparts in A and B and on the current principal joint transform of C. This is the case of the cross-transforms, the principal handles transforms, the inverse inertias, and the quadratic coefficients of the total acceleration equations.

These quantities depend on the joints positions and the accelerationindependent forces, as can be seen from their recursive expressions (2) and ( 9). Precisely:

• The principal joint transform T C

• of a node's principal joint depends on the principal joint's position q • .

• The cross-transforms T C i j of an internal node depend on the position of the node's principal joint, and the cross-transforms of its child nodes.

• The principal handle transform T C B of a node B with parent C depends on the node's cross-transforms, and potentially the principal joint transform of the node's parent if C's principal handle is not a handle of B.

• The inverse inertias Φ C 1 , Φ C 2 , Φ C 12 and Φ C 21 and the quadratic coefficients Ψ C 1 , Ψ C 2 , Ψ C 12 and Ψ C 21 of a node C depend on the position of the node's principal joint, and the inverse inertias and quadratic coefficients of the child nodes. The inverse inertias and the quadratic coefficients of a leaf node are constant.

• The bias accelerations b C

1 and b C 2 , the linear coefficients p C 1 and p C 2 , and the constant coefficient η C of a node have a similar dependency on the joint positions, the bias accelerations, the linear and the constant coefficients of their child nodes, but they also depend on the acceleration-independent forces. If an active joint force Q or an external force f k applied to a leaf node change, these quantities and the corresponding ones in all the parent nodes have to be updated.

Assuming an explicit integration scheme, the joints positions q • only depend on their acceleration q• . Since, during the (partial) back-substitution pass, an approximate linkage acceleration s C has been obtained by computing exactly some joint accelerations, and approximating the accelerations of the other joints as zero, we only need to update the positions q • , and the principal joint transforms T C

• , of the joints whose acceleration has been computed. If N E joints accelerations have been computed, the complexity of this update is O(N E ).

Similarly, since the handles cross-transforms, the principal handles transforms, the inverse inertias and the quadratic coefficients of a node only depend on the principal joint's position and their counterparts in the node's children, we only need to update these coefficients in the nodes whose principal joint's acceleration has been computed, and the complexity of this update is O(N E ).

Because the remaining coefficients also depend on the principal joints positions, we need to update them at least where the joints accelerations have been computed. Thus, the complexity of this update is at least O(N E ). However, we also need to take into account changes in the acceleration-independent forces. We do so by first retrieving the appropriate 4 × 4 transforms and then converting them to 6 × 6 spatial force transforms.

Assume the leaf node is C l and that the nodes on the path from the root to C l are C 0 (the root node), C 1 , . . . , C l-1 and C l . The transform T 0 l from the principal handle of C l to the world frame is obtained by accumulating the principal handle transforms on the path to C l :

T 0 l = T 0 1 . . . T l-1 l ,
and the transform T l 0 from the world frame to the leaf node's principal handle is obtained by inverting T 0 l . The transforms from the world to all the leaf node's handles are then obtained by multiplying T 0 l by the appropriate handles cross-transforms of C l , and are converted to 6 × 6 spatial force transforms to express f k in the reference frames of the leaf node's handles. The complexity of an update in an applied external force f k is thus also O(log(N J )).

Although each separate update would have a O(log(N J )) complexity, it is possible to perform the updates level by level as in [START_REF] Lotan | Efficient maintenance and self-collision testing for kinematic chains[END_REF], from the leaf nodes to the root node, and avoid processing a same node twice. Assuming there are N U nodes where either an active joint force Q or an applied external force f k is updated, the complexity of the coefficients update is then O(N U log(N J /N U )). This complexity grows to O(N J ) when N J nodes require an update in an acceleration-independent force Q or f k .

Thanks to this hierarchical state representation, it is thus possible to limit the update of the coefficients in the linkage state and the total acceleration equations to some nodes only. The number of these required nodes is at least O(N E ), where N E is the number of joint accelerations which have been computed during the backsubstitution pass to obtain the approximate acceleration, and depend on the number of acceleration-independent forces that have to be updated.

Before the beginning of the first time-step, the coefficients in the linkage state and the total acceleration equations are computed for all nodes in the assembly tree. Then, given an user-defined error threshold ε max , our algorithm repeatedly computes an approximate acceleration, with an acceleration error smaller than ε max , and performs the simplified coefficient update based on the approximate acceleration and the changes in the acceleration-independent forces.

Results and Discussion

In this section, we present several preliminary benchmarks used to test our algorithm. The algorithm has been implemented in C++, and the benchmarks have been performed on a 2.8 GHz Pentium PC with 2GBytes of RAM. In all benchmarks, the error metric we used was the relative joint acceleration error, which is independent of both the linkage size and the magnitude of the linkage acceleration.

Bicycle Chain

We first provide an informal, "visual" description of the effect of the acceleration error threshold on the quasi-static simulation of a linkage. For this first test, we used a bicycle chain composed of 30 links, and applied a constant force on the last link of the chain for

Millipede Robot

Our second benchmark is a millipede-like robot. The robot has 1,000 legs with 10 links each and a spine of 3,000 links, for a total of 13,000 links and 12,999 joints. Figure 5 In this benchmark, we have applied a single force to a leg of the millipede during 4,000 time steps. We have performed the simulation both with the DCA and our algorithm for various precisions (ε max = 2 -i , 0 i 19), and we have measured the average running times for each simulation, per time step. We have also measured the average number of joint accelerations that our algorithm had to compute in order to guarantee that the relative joint acceleration error was below ε max . The results are plotted in Figure 6. As expected, as more and more joint accelerations need to be computed when ε max decreases, the average running time of our algorithm increases as well. When too many nodes have to be updated (i.e. when the error threshold becomes smaller than 10 -4 ), our algorithm becomes less efficient than the DCA due to the higher cost of a node update. Below this threshold, however, our algorithm can quickly obtain an approximate simulation. Figure 6: Average running time of our algorithm when a force is applied to the millipede robot. As the maximum error decreases, more and more joint accelerations have to be computed and the average running time of our algorithm increases. When too many nodes have to be updated, our algorithm becomes less efficient than the DCA. Below this threshold, however, our algorithm can quickly obtain an approximate simulation.

Molecule model

Our third benchmark involves a randomly generated molecule-like linkage with 50,000 rotational degrees of freedom (cf Figure 7).

In this benchmark, we have varied the number of forces N f applied to the linkage during the simulation (1, 10, 100, and 1,000), as well as the desired precision (ε max = 2 -i , 0 i 19). For a given number of forces N f , we have selected a random initial linkage configuration, have performed 4,000 simulation steps with the DCA, and have determined the average time required to perform one time step with the DCA. This average running time was approximately independent of the number of applied forces and found to be about 875 milliseconds per time step. We have then performed 4,000 simulation steps with our algorithm for each precision ε max , using the same initial configuration each time the simulation was restarted. For each pair of parameters (N f , ε max ), we have determined the average running time of our algorithm per time step, and the average number of nodes for which the joint acceleration has been computed during the back-substitution pass. The results are plotted in the graphs shown in Figure 8, and some of the data are also given in Table 1.

Discussion

Several remarks can be made:

• Performance vs precision: The average running time and the average number of computed joint accelerations increase with the required precision. It can be noted that the running time and the number of nodes processed by our algorithm do not seem to increase much after the required precision is smaller than a certain threshold. This observation suggests that it might be possible to define some motion complexity measures, which would characterize the behavior of complex linkages and could provide some insight when designing linkages. For example, comparing Figures 6 and 8 suggests that the impact of a force on a molecule-like linkage is more local than in the case of the millipede robot: for a given precision, the ratio of the number of joints affected by an applied force to the total number of joints is smaller in the case of the molecule-like linkage.

• Performance vs number of applied forces: There is a clear correlation between the number of acceleration-independent forces that have to be updated at each time step, and the performance of our algorithm. As the number of applied forces increases, more and more joint accelerations have to be computed for a given precision, and more and more nodes have to be processed in the lazy linkage update step.

• Performance improvement: Even when the required precision is high (about 10 -6 ), our algorithm offers a significant performance improvement over the DCA. In the case of the molecule-like linkage, when 1,000 forces are applied at each time step, the average running time of our algorithm when ε max = 2 -19 10 -6 is about 50 milliseconds, resulting in about 17.5 times performance improvement. When only one force is applied, the average running time of our algorithm with the same precision is about 2.5 milliseconds, providing a factor of 350 times speed-up.

Relative Displacement Error. Because an error in a single joint acceleration may induce a large error in many body positions expressed in world coordinates, especially with large linkages and rotational degrees of freedom, we need to make sure that, although we compute the joint accelerations only up to some precision, the approximation still allows us to obtain a good approximation of the linkage configuration in the world frame. In this regard, the molecule-like linkage is the most challenging of the three benchmarks presented here. In order to quantify the quality of the simulation provided by our algorithm in this case, we have thus also measured the relative displacement error of the simulated linkage. Even with high precision requirements, however, our algorithm provides up to a factor of 350 times speed-up over the DCA.

Error Threshold DCA The table shows the average running time of our algorithm (in milliseconds), the number of joint accelerations computed, and the relative displacement error when the maximum error threshold ε max and the number of forces N f applied to the linkage vary. By rapidly determining the joints that contribute most to the linkage acceleration, our algorithm is able to obtain a performance gain up to two orders of magnitude over the DCA, with an error threshold of 2 -19 10 -6 .

N f 2 -3 2 -7 2 -11 2 -15 2 -19
Whenever a simulation with parameters (N f , ε max ) completed, we determined the displacement d i (N f , ε max ) incurred by the center of gravity of each rigid body i in the linkage during the simulation, in world coordinates. By comparing with the displacement d i (N f ) obtained with the DCA, we determined the relative displacement error e(N f , ε max ) as follows:

e(N f , ε max ) = max i ||d i (N f , ε max )|| max i ||d i (N f )||

Figure 9 shows how the average relative displacement error per time step evolves when the precision on the joints accelerations increases. As expected, increasing the precision of the computation of the joint accelerations increases the precision on the configuration of the molecule expressed in world coordinates. Specifically, enforcing a precision of about 10 -6 on the joint accelerations allows us to obtain a final linkage configuration with a relative displacement error smaller than 10 -6 . Considering the previous performance benchmarks, we conclude that an application for which this level of precision is sufficient can benefit substantially from our algorithm.

Conclusion

We have introduced a novel algorithm for efficiently computing an approximation of the forward quasi-statics of an acyclic linkage. Given any user-defined error threshold and acceleration error measure, our algorithm is able, at each time step, to rapidly determine the joints which contribute most to the overall linkage motion, and determine an approximate linkage acceleration within the required error bound. Our algorithm allows a natural trade-off between the precision of the resulting simulation and the time required to compute it. We have implemented our algorithm and tested its performance on complex benchmarks consisting of up to 50,000 links. We demonstrate that in some cases our algorithm is able to achieve up to two orders of magnitude of performance improvement, while providing a high-precision, error-bounded approximation of the quasi-statics of the simulated linkage.

We would like to extend this work to mechanisms with loops. In our current framework, indeed, loops can only be handled through corrective penalty forces, i.e. forces applied to correct drifts where the linkage should be closed. Because the general version of the DCA can handle loops, however (at the price of an increased computational complexity [Featherstone 1999b]), we believe that it might be possible to extend our algorithm to handle loops analytically, and not only through penalty forces.

An important extension of this work is the generalization of the proposed algorithm to simulate the full dynamics of complex linkages. This extension is non-trivial, because of the presence of velocitydependent coefficients in the articulated-body equations and the total acceleration equations. Independent of the reference frames in which they are expressed, the body velocities relatively to the world frame have to be available at each time, which, a priori, would require updating all the equations coefficients.

We have recently shown that a simplified update mechanism is still possible in the dynamics case, and have proposed an adaptive articulated-body dynamics algorithm [Redon et al. 2005a]. Nonetheless, we believe progressive quasi-statics can be useful in several areas, and the present work suggests several future research directions:

• Approximate collision detection and response: We would like to complement this work with new approximate methods for collision detection and response. For example, we would like to investigate partial update mechanisms for boundingvolume hierarchies typically used in collision detection [START_REF] Redon | Fast continuous collision detection between rigid bodies[END_REF]Redon et al. 2005b].

• Progressive quasi-statics with force feedback: We believe the progressive quasi-statics framework introduced in this paper might be useful for fast approximate force-feedback simulation of complex linkages, which typically require update frequencies greater than 1000Hz (as opposed to 25Hz for visual simulation).

• Perception studies: The precision-performance trade-off allowed by our framework seems to make it possible to study human perception of articulated body quasi-static motion, through precise evaluations of the perceptibility of errors in the motion of an articulated body. This would allow for effective and pertinent simplification of visual and haptic (forcefeedback) articulated-body simulation.

Finally, beside CAD/CAM, we would also like to investigate the application of our approach in drug design, digital actors, robotics and automation, computational biology, and design of medical instrumentation.

Figure 1 :

 1 Figure 1: Error-bounded approximation of a complex linkage's quasi-statics. (a) A close-up view of the model. (b) The initial and final positions of the linkage. (c)-(e) The simulation results of our algorithm vs. the DCA for varying degrees of precision, with both the view of the chain and a close-up of the end position.

Figure 2 :

 2 Figure 2: An articulated body C is formed by connecting two articulated bodies A and B with a joint. (a) A rigid body with three handles. (b) The subassembly A has four handles H A 1 1 , H A 1 2 , H A 1 3 and H A 2 , while the subassembly B has three handles H B 1 , H B 2 1 and H B 2 2 . The handles H A 2 and H B 1 are consumed during the assembly operation, and only five handles remain in C. (c) The assembly tree describes the hierarchy of articulated bodies.

  Figure 2(b)-(c) shows an example of such an assembly operation. The goal of the main pass is to obtain the articulated-body equation of C from the articulated-body equations of A and B.

Figure 3 :

 3 Figure 3: The three types of transforms maintained during the simulation. Left: The cross-transforms T A i j and T B i j of A and B and the principal joint transform T C • of C, from which the cross-transforms of C can be computed. Right: The principal handles transforms of A and B, which relate the principal handles of A and B to the principal handle of C (assumed to be H C 3 in this case), can be computed from the cross-transforms of A and B and the principal joint transform T C• .

Figure 5 :

 5 Figure 5: A millipede-like robot with 13,000 links. (a) The building blocks of the millipede robot. (b) A force f k is applied on a leg of the millipede during the simulation.

  (a) shows a close-up on the building blocks of the millipede, and Figure5(b)shows a snapshot of the simulation.

Figure 8 :

 8 Figure8: Performance of our algorithm on a molecule-like model with 50,000 rotational degrees of freedom, when the number of forces applied to the linkage and the required precision vary. As the precision enforced on the computation of the joint accelerations increases, the average running time of our algorithm increases as well (left), because more joint accelerations have to be computed (right). Even with high precision requirements, however, our algorithm provides up to a factor of 350 times speed-up over the DCA.

Figure 7 :

 7 Figure 7: A randomly generated molecule-like linkage with 50,000 rotational degrees of freedom. (a) The complete model. (b) A close-up view.

Figure 9 :

 9 Figure 9: Influence of the relative joint acceleration error on the relative displacement error. As the precision enforced on the computation of the joint acceleration increases, the precision in the final configuration of the linkage expressed in the world frame increases as well.

Table 1 :

 1 Performance of our algorithm with the molecule model.

			Time (milliseconds)
	1	0.6	0.9	1.4	1.9	2.4	875.1
	10	0.9	1.3	1.8	2.6	3.1	873.0
	100	2.7	4.1	6.4	7.8	8.4	875.0
	1,000	17.8	21.4	42.6	49.2	49.9	880.0
		Number of computed joint accelerations
	1	25	37	47	59	73
	10	32	45	56	75	85	50,000
	100	39	80	145	189	206
	1,000	41	149	830	1,060	1,086
			Relative displacement error
	1	6E-05 2E-06 1E-07 4E-09 2E-09
	10	1E-04 1E-04 4E-07 4E-08 1E-08	0.0
	100	8E-04 6E-04 7E-05 3E-06 1E-07
	1,000 6E-03 2E-03 7E-05 2E-06 6E-08

Acknowledgements

The authors would like to acknowledge Dr. Roy Featherstone for his help with the DCA. This project is supported in part by the Army Research Office, Intel Corporation, the National Science Foundation, and the Office of Naval Research.


Li-Thiao-Té Sébastien
email: lithiaote@cmla.ens-cachan.fr

Semiparametric Estimation of a Noise Model with Quantization Errors

Keywords: parametric estimation, quantization effects, life sciences signal processing Mathematics Subject Classification (2000) MSC 62F10 Point estimation (Parametric inference), 65G30 Interval

The detectors in mass spectrometers are precise enough to count ion events, but in practice large quantization errors affect the observations. To study the statistics of low intensity chemical noise, we model the detector signal as X = ⌊τ N⌋ and estimate both τ and N in a semi-parametric approach where the integer valued random variable N represents the number of ions and τ represents the gain parameter of the detector. When τ ≤ 1, we explain why the gain parameter cannot be recovered without a priori information on N. When τ > 1 however, N can be deduced from X and a sufficiently precise estimate of τ. To perform parametric estimation of τ, we first study simple estimators which provide useful upper bounds. We then introduce the concept of compatible lattices and we derive an optimal estimator that is independent of the law of N.

Introduction

Ion Detectors

Mass spectrometers are instruments that ionize the compounds of a sample, separate the ions, then quantify the ions at each mass to charge ratio. The resulting signal is a histogram that represents ion intensity as a function of the mass to charge ratio of the ions. With sufficient precision in the separation and the mass to charge measurement, the components of the sample can be identified and quantified. Mass spectrometers are widely used for analysing very diverse mixtures, e.g. detecting explosives for airport security or analysing oil products. See [START_REF] Aebersold | Mass spectrometry-based proteomics[END_REF][START_REF] Lane | Mass spectrometry-based proteomics in the life sciences[END_REF][START_REF] Siuzdak | Mass spectrometry for biotechnology[END_REF] for an introduction to mass spectrometry in the life sciences.

We consider detectors similar to microchannel plate detectors that are used in most mass spectrometers [START_REF] Wiza | Microchannel Plate Detectors[END_REF]. When an ion hits the detector plate, it produces an analog signal that is amplified, quantized, then reported to the computer. The level of quantization is quite high as there may be only 2 11 = 2048 levels1 in some instruments, and small signals as well as chemical noise are strongly affected by quantization effects. Specific difficulties have appeared with high-throughput analyses of biological material. In particular, biological samples may contain trace amounts of molecules of interest. These are difficult to distinguish from chemical noise which produces patterns similar to real signals [START_REF] Andreev | A universal denoising and peak picking algorithm for LC-MS based on matched filtration in the chromatographic time domain[END_REF][START_REF] Krutchinsky | On the nature of the chemical noise in MALDI mass spectra[END_REF][START_REF] Windig | A noise and background reduction method for component detection in liquid chromatography/mass spectrometry[END_REF]. In 2004, [START_REF] Anderle | Quantifying reproducibility for differential proteomics: noise analysis for protein liquid chromatographymass spectrometry of human serum[END_REF] suggested Poisson-like behaviour for the ion intensity based on a linear relationship between the mean and variance of the noise. This linear relationship suggests that the amplification factor of the detector may be unaccounted for in the data set.

To study chemical noise in the experimental data, we interpret the amplification factor as an overdispersion parameter in a semi-parametric approach. To study chemical noise in the experimental data, we estimate the amplification factor and an unknown distribution for the chemical noise in a semi-parametric approach. Let N denote the number of chemical noise ions that reach the detector, we consider the following observation model:

X = ⌊τN + ε⌋
where the noisy signal τN + ε is truncated before observation. τ represents the amplification factor of the detector and ε represents electronic noise. In this paper, we make the assumption that ε = 0 or equivalently that there are only quantization errors in the measurements. The observation model is associated with the statistical structure (N, B, P θ ), with θ = (τ, N) where τ is a positive real number and N is a probability distribution on N.

We believe a priori that N is Poisson distributed as it models rare events (ion counts). Consequently, τ can be interpreted as an overdispersion parameter affecting Poisson distributed observations. This has been tackled in the framework of double exponential families, as presented by Efron in [START_REF] Efron | Double exponential families and their use in generalized linear regression[END_REF]. In [START_REF] Antoniadis | Estimation semi-paramétrique dans les familles doublement poissonniennes et application aux spectres de diffraction[END_REF], Antoniadis et al use double exponential families in a regression model to analyse diffraction spectra. This corresponds to estimating the regression function µ i in the model X ′ i = τN i where N i is Poisson distributed with varying parameter µ i . However, this framework does not explicitly take into account quantization errors and thus provides poor parameter estimates as we will show in Section 4.2. Moreover, we wish to confirm the Poisson hypothesis using non parametric estimation.

Our approach is to first estimate τ given a set of observations of X, then deduce the distribution of N from the estimate. We show that the estimate is precise enough to allow complete disambiguation of the observations.

Estimation of the Ion Statistics

With negligible quantization error, the observation model becomes X = τN. Estimation of τ is trivial; all that is required is to observe the event {N = 1} i.e. {X = τ × 1}, or the two events {x 1 = τi} and {x 2 = τ(i + 1)} and compute the difference x 2x 1 .

To recover N, it then suffices to consider X/τ. The quantization error may be neglected when τ ≫ 1 in the observation model X = ⌊τ N⌋ and the previous estimates provide τ with a precision on the order of the quantization error.

In the general case, we can recover the samples of N from the samples of X when the mapping x → ⌊τ x⌋ is injective. The inverse mapping is y → ⌈y/τ⌉. We call this situation the distinguishible case. It occurs if and only if τ ≥ 1 (see proof in the Appendix, Prop 6). In this situation, the semi-parametric approach can be separated into parametric estimation of the gain parameter τ then non parametric estimation of the distribution of N from iid samples.

Example 1 > data = floor(1.32 * n) % Distinguishible case > n
[1] 1 2 3 4 5 6 7 8 9 10 > data

[1] 1 2 3 5 6 7 9 10 11 13

When τ is smaller than 1, the truncation error merges adjacent values of N. In the following example, the events {N = 3} and {N = 4} cannot be distinguished in the data set. This is because the corresponding observation is {X = 2} in both cases. In the distinguishible case, it is natural to sort and index the observed values in order to determine the mapping x → ⌊τ x⌋. This is not sufficient in practice because of missing values or outliers which can modify the indexes.

Observation Set

The gain parameter and the law of N have separate effects on X. In the distinguishible case, the distribution function of X is a transformation of the distribution function of N by the mapping x → ⌊τ x⌋. The gain parameter and the truncation error only distort the position of each peak, whereas the relative frequencies are unchanged.

Consequently, the support S of the empirical distribution is sufficient information for estimating τ whereas the empirical frequencies are sufficient information for the distribution of N.

In the non distinguishible case, the set of observed integers is always N for large samples (see Section 3). As a consequence, τ cannot be estimated based on that set alone. A semi-parametric approach is not feasible either. For instance, we cannot estimate the mean E[N] but only E[X] = τE[N]. To separate τ and N, we have to provide prior assumptions on the distribution of N like a Poisson parametric family.

In the following, we study properties of the set S of observed integers. This set can be constructed from the dataset in O(n log(n)) time using sorting for example. The algorithmic complexity of the following algorithms is governed by the size of S , and in particular, the maximum integer in S .

We focus on the distinguishible case, and perform parametric estimation of the gain parameter from a random set of integers. As the support of the empirical distribution is a sufficient statistic for τ, we use the statistical structure Ω = 2 N , T, P τ , τ ∈ ]1, +∞[ where Ω is the power set of N and T is the exhaustive σ -algebra on Ω . (P τ , τ ∈ ]1, +∞[ ) is a parametric family of distributions on Ω that is implicitly generated in the following way. For a fixed integer n and fixed but unknown integer-valued random variable N, P τ is the distribution of the random variable S which is the set of observed integers in an independent identically distributed sample (X 1 , . . . , X n ) of X = ⌊τN⌋.

Organization of the paper

To estimate τ in the distinguishible case, we first provide simple estimators for τ in Section 2. These are later used as a starting point for improved estimators and to restrict the search space for τ.

In Section 3, we define the notion of compatible values and provide a few properties of the set of compatible values. In particular, the true parameter τ is a compatible value and is close to the highest compatible value. This leads to an optimal estimator that is described in 3.3.

We show the results of some simulations in section 4 and compare with the Maximum Likelihood Estimator obtained from the Double Poisson Family, an estimator based on linear regression and another one based on Fourier transform.

Estimators and Upper Bounds for τ

The results in this section are based on the following idea. Two points in S are separated by at least ⌊τ⌋. Consequently, when τ is large, then S is a sparse set, whereas S is dense when τ is near 1. For instance, there are consecutive points in S if and only if τ ≤ 2 (see Proposition 7 in the Appendix).

A better estimate can be obtained by combining more than 2 consecutive points. Let x, y denote the set of integers between x and y. If x, y is a subset of S , then

τ < 1 + 1 y -x
. Consequently, τ can be estimated by 1 + 1 yx with a precision on the order of the inverse of the length of the interval 1 yx . However, this estimator is strongly affected by missing values in S .

Instead of considering all the segments in S , we propose to use the overall density of the set, which is easier to compute algorithmically. Let x = ⌊τ n⌋ denote the largest integer in S . Then τ <

x + 1 n . When n is unknown (because of potential missing values), let n denote the number of non zero observed integers i.e. the number of elements in S . Then τ <

x + 1 n ≤

x + 1 n .

Consequently,

x + 1 n is an estimate of τ with precision on the order of 1/n (without missing values). As it uses the whole data, it is usually more precise than the previous bound. We will use this in the rest of the paper to restrict the search space for τ.

Let us compare the previous bounds on an example. Suppose that τ = 1.32 and S = {1, 2, 3, 5, 6, 7, 9, 10, 11, 13}. As there are consecutive integers in S we obtain τ < 2 . Using the interval 5, 7 , we obtain τ < 1 + 1/2. Using the interval 9, 11 , we obtain τ < 1 + 1/2 as well. The density upper bound is τ < 14/10.

Remark

We only provided upper bounds in this section because lower bounds can only be deduced from the integers that cannot be generated in the model. These are difficult to distinguish from missing values, which are integers that can be generated in the model, but do not appear in the set S of observed integers.

Compatible Values

The upper bounds that we proposed in the previous section are easy to compute but rather poor because they only take into account the proportion of observed integers. In this section, we describe an algorithm with higher computational load but which can leverage the information in the location of each observed integer in the data set.

Lattices of Integers

In the observation model X = ⌊τN⌋ where N is integer valued, only specific integers can be generated. Given a strictly positive real number t, let us define the set of possible values for x as the lattice associated to t, i.e. the infinite set of integers S t = {x = ⌊t k⌋, k ∈ N}. The set of observed integers S is also called the empirical lattice.

With infinitely many observations, the parameter τ is completely characterized by the empirical lattice as the following proposition shows. This justifies that S is sufficient information for estimating τ.

n⌋ = ⌊t 1 n 1 ⌋ with n 1 > n.
In the latter case, distinguishibility implies that there are strictly more elements in S t 1 ∩ A than in S t 2 ∩ A where A denotes the set of integers 0, ⌊t 2 n⌋ .

The Set of Compatible Values

Proposition 1 is not sufficient for estimating τ because in practice we only observe a finite set S S τ . Consequently we define the notion of compatible lattices and equivalently compatible values. For any positive real t, we say that t is compatible with the data if S ⊂ S t . Likewise, for any two sets A and B, A is compatible with B if B ⊂ A. Being compatible with the data set is a necessary condition for a valid estimator of τ.

The set of values that are compatible with the infinite lattice S τ is adequate for estimating τ because of the following proposition.

Proposition 2 τ is the largest real number in C(S τ ).

Proof τ is a compatible value, we only have to show that it is the largest. Let u denote a real number greater than τ, and let α denote a positive real number such that τ < τ + α < u. We will prove that u is not compatible with τ by constructing an element in S τ that cannot be in S u . Let a denote a positive integer such that a > 1 α , and n = ⌊τa⌋. Suppose that S τ ⊂ S u then n belongs to S u , and there exists a positive integer a ′ such that n = ⌊ua ′ ⌋. a ′ ≥ a because in the distinguishible case, a and a ′ correspond to their indices in the sets S τ and S u and S τ ⊂ S u . Moreover, as τ ≤ u we have ⌊τa⌋ ≤ ⌊ua⌋ ≤ ⌊ua ′ ⌋. For all three terms to be equal to n in the distinguishible case requires that a = a ′ . Consequently, both τ and u lie in the interval

[ n a , n+1 a [. As a result, |τ -u| ≤ 1 a which contradicts a > 1 α .
The set C(S τ ) has an intricate structure. It contains the positive real numbers smaller than 1 and the harmonics τ k , k ∈ N * , but these are not the only values. For example, 4/3 is compatible with 2 because every even integer can be written as ⌊k × 4/3⌋, k ∈ N. Indeed, let k be an even integer. Either k is a multiple of 4, in which case

k = 4i = ⌊ 4 3 × 3i⌋, or k = 4i + 2 = ⌊ 4 3 × (3i + 2)⌋.

Estimation with a Finite Lattice

In practice, the empirical lattice is finite and can contain missing values and outliers. We say that an integer is missing from S when it is in the theoretical lattice S τ , smaller than x = max S , but not in S . The set of compatible values with S is a finite union of intervals and compatible values are never isolated. As the data contains less information, the true parameter τ is not the supremum of C(S ), but it is still maximal in the following sense. Proof Let t > 1 denote a compatible value. For each observed value x ∈ S , there exists an integer n such than x = ⌊tn⌋. Consequently,

t verifies t ∈ [ x n , x+1 n [. The inter- section of the constraints t ∈ [ x n , x+1 n [ for all x ∈ S is an interval t ∈ [ x 1 n 1 , x 2 n 2 [. All values t ∈ [ x 1 n 1 , x 2 n 2
[ verify all of the constraints and are thus compatible. The length of this interval is x 2 n 2 -x 1 n 1 which is at least 1/(n 1 n 2 ). In the distinguishible case, n 1 < max S and n 2 < max S , which implies that the length is at least 1/(max S ) 2 .

Lemma 2 Let t be a positive real number that is compatible with the empirical lattice. Then t <

x + 1 n .

Proof This follows directly from the upper bounds in Section 2. See Proposition 9 in the Appendix.

To complete the proof, it suffices to show that τ belongs to the largest interval.

Proof

′ < t. For all n ∈ N , ⌊nτ⌋ ≤ ⌊nt⌋. As t / ∈ [a ′ , b ′ [, t
breaks at least one of the constraints, that is to say, there is an integer x in S such that x = ⌊nτ⌋ < ⌊nt⌋. x is skipped in S t and thus t is not compatible.

The previous proposition suggests that it suffices to find the largest compatible interval to estimate τ, and this is our proposed estimator τ. More precisely, the set C(S ) = ∪ J j=1 [a j , b j [ is a union of J intervals, with (a j ) and (b j ) increasing sequences, then τ = a J + b J 2 .

We use the following algorithm to compute τ. This also computes the mapping x = ⌊τn⌋ → n and the precision.

-compute the set of observed values by sorting the data set and removing multiple occurences -compute the upper bound τ < B =

x + 1 n where n = cardS -find an approximation of the largest compatible value t by testing the compatibility of the real numbers t k = B -k x2 -deduce the indexes from t, that is to say for all x ∈ S , find i such that x = ⌊ti⌋

-compute the interval [a, b[ as the intersection of the constraints t ∈ x i , x + 1 i , for
all x in S -return a + b 2 as an estimator for τ

Properties of the Estimator

According to the previous results, the estimator performs well when there are no missing values or outliers. Its precision is (ba)/2 and can be computed inside the algorithm. The precision is at least 1/ n, but depending on the value of τ it can reach a precision on the order of 1/ n2 . In all cases, the precision is better than the density bound, and there is a lower bound.

If there are missing values or outliers, the algorithm may find an interval of compatible values that does not contain τ. For example, if the dataset is {0, 2, 4, 6, 8}, a reasonable estimator would answer 2 and not τ = 4/3 with missing values 1 and 5. In practice, such cases are rare, and are related to arithmetic properties of the set S . However, the largest compatible value is never an erroneous answer to the problem. It is a parcimonious answer in the sense that it is the smallest lattice which may explain the dataset.

The estimator is optimal in the sense that the algorithm finds an interval of positive real numbers that are all plausible. Given a dataset (x 1 = ⌊τ i 1 ⌋, . . . , x n = ⌊τ i n ⌋) of size n, there is an interval of compatible values that can generate (x 1 , . . . , x n ) from the same realization (i 1 , . . . , i n ) of N. Let [a J , b J [ with b J = sup C(S ), the following proposition holds.

Proposition 4 Given a realization

(i 1 , . . . , i n ) of N, all values in [a J , b J [ generate the same data set (x 1 , . . . , x n ), i.e. ∀t ∈ [a J , b J [, ∀ j ∈ 1, n , x j = ⌊τ i j ⌋ = ⌊t i j ⌋ Proof As in the proof of Proposition 3, [a J , b J [ is the intersection of the constraints x j = ⌊t i j ⌋.
The data set does not contain enough information to distinguish the values in [a J , b J [. In particular, even if the realization (i 1 , . . . , i n ) is given, then the values are not distinguishible. Note that if x 0 is known not to be in S τ , then for all integers i, τ ≥ x 0 +1 i or τ < x 0 i . These inequalities are not informative because they are already contained in x = ⌊ti⌋, ∀x ∈ S .

The program is quite fast. First because is relies only on the set S which is much smaller than the dataset when τ is near 1 and N is independent identically distributed, because repeats of N are discarded. As the following proposition shows, with few missing values, the density bound is precise and the algorithm is quicker. All compatible values can be retrieved by testing B x2 numbers.

k ncardS + 1 n( x+1)
which can be approximated by C ≃ τ 2 x(k + 1 τ ).

Testing for the compatibility of a real t is linear in the size of S , so the whole procedure is at most quadratic. The full set of compatible values can be obtained in cubic time. The vertical axis represents values of τ. The set of compatible values is composed of several intervals and represented on the left. For each interval, we select one compatible value t and represent the lattice S t . All reals in the same interval generate the same lattice, up to max(S ).

Results and Discussion

Compatible Values Estimator

For comparison, Figure 2 displays S t for several values that are not compatible with the data. For example, 5 and 11 are in the dataset but not in S 1.2 .

Two sources of variation affect the estimate τ. First, the estimator is not perfect because the dataset is finite. Second, the data set S is random. Figure 3 shows the performance of the estimator with a fixed dataset (N ∈ 1, 10 ) for several values of τ. The intervals shown correspond to the intervals in C(S ) that contain the largest compatible value.

We can see that the precision of the estimator varies with τ. Only the range [START_REF] Aebersold | Mass spectrometry-based proteomics[END_REF][START_REF] Anderle | Quantifying reproducibility for differential proteomics: noise analysis for protein liquid chromatographymass spectrometry of human serum[END_REF] is shown because the precision only depends on the rest τ -⌊τ⌋ modulo 1. Consequently, the absolute precision is roughly constant, whereas the relative precision is

O 1 τ
. With small quantization error (τ ≫ 1) the estimation problem is easier. Figure 4 shows the distribution of τ when the dataset S is the result of 15 samples of X = ⌊τN⌋ where τ = 1.32 and N is distributed according to a Poisson random variable with mean 5.5. The distribution of the estimator value is obtained from the 200 repeats shown in the bottom of the plot thanks to a kernel estimate, even if the distribution is a sum of Dirac point masses. The interval shows the interval obtained with the (complete) dataset {⌊1.32 * n⌋, n ∈ 1, 13 }.

The Double Poisson Family

In this section, we briefly recall the results from [START_REF] Efron | Double exponential families and their use in generalized linear regression[END_REF], and deduce an estimator for our model. Let g µ (y) = e -µ µ y /y! denote the distribution function of a Poisson ran- dom variable with mean µ. The double Poisson distribution with parameters θ , µ is defined as: where c is a normalization constant.

f θ ,µ (y) = c(θ , µ)θ 1/2 g µ (y) θ g y (y) 1-θ = c(θ , µ) θ 1/2 e -θ µ e -
Maximum Likelihood Estimation leads to the following estimators. Let (Y 1 , . . . ,Y n ) be independent identically distributed random variables with distribution f θ ,µ , then

μ = 1 n n ∑ i=1 Y i θ = n 2 ∑ n i=1 I(Y j , µ)
where

I(µ 1 , µ 2 ) = µ 1 (log(µ 1 ) -log(µ 2 )) -(µ 1 -µ 2 ).
Let Y θ ,µ be a random variable with distribution function f θ ,µ , then according to [START_REF] Efron | Double exponential families and their use in generalized linear regression[END_REF] Y θ ,µ has approximately the same distribution as X/θ where X is Poisson distributed with mean µθ . With Poisson distribution for the ion counts, our observation model becomes Y = ⌊τN⌋ where N is Poisson distributed with mean λ . Consequently, estimates for τ and λ can be deduced from θ and μ with the following relations:

τ = 1 θ λ = μ θ
The double Poisson distribution is a correct approximation of the distribution of X/θ for large µ, and in that case, τ and λ are unbiased estimates of τ and λ . The standard deviation of τ is τ √ 2 √ n . Figure 5 shows the distribution of τ with flooring noise, i.e. in the model Y = ⌊τN⌋ (solid line) and without flooring noise in the model Y = τN (dotted line). The plot was generated with 2000 repeats with data sets of size 500. We observe a large standard deviation compared with the compatible values estimator, even on a much larger data set. With large values of µ, truncation has limited effect on the estimate. For modeling rare ion count events, we need to study the estimators with small values of λ and τ. In that case, τ is strongly biased for both models as shown on Figure 6. This implies that the approximation is not suited to this range of parameters, and that the flooring noise makes a significant difference there. Figure 6 was generated using 2000 repeats with data sets of size 500. For comparison, we show the optimal interval obtained by the compatible values estimator on the data set {⌊1.32 * n⌋, n ∈ 1, 13 }.

Fourier Estimator

From the set τN we can construct the signal f : t → ∑ k∈N δ (x -τk) where δ denotes the Dirac function, that is to say a periodic series of pulses. The period τ may thus be estimated using Fourier transform. Likewise, we define the estimator 1/τ F as the maximum of the Fourier transform of the quasi-periodic signal f : t → ∑ k∈N δ (x -⌊τk⌋).

As τ can be seen as a quasi-period, our estimation problem is closely linked to the "harmonic retrieval problem". Many approaches have been proposed in that domain and the main focus is on the estimation of the Power Spectral Density [START_REF] Hayes | Statistical Digital Signal Processing and Modeling[END_REF]. However, the signal is usually perturbed by additive noise whereas in this paper we consider a distortion of the time axis.

We use the following algorithm:

-sample the signal f at the points x i = i for the integers i in 0, max(S ) -compute the Discrete Fourier Transform -compute an upper bound using Proposition 9 : τ < B =

x + 1 n -find the frequency with highest absolute Fourier coefficient -return the corresponding period (inverse of the frequency) This estimator has a precision that corresponds to the sampling rate in time space around the true value. In the Fourier space, the sampling rate is uniform with steps of length 1/ max(S ) which is equivalent to 1/(τ × x). In the time space, as P = 1/ f , then ∆ P = -∆ f / f 2 and the sampling rate is non uniform. For f = 1/τ we obtain the precision of the Fourier estimator as τ/ x. This suggests that the precision decreases with τ. However, the signal frequency 1/τ is near x/ max(S ) which is one of the sampling points. As a result, in practice, the absolute precision is on the order of 1/ x and independent of τ. Remark When oversampling by a factor k, i.e. sampling at the points x i = i k for the integers i in 0, max(S ) × k , the harmonics of 1 Hz increase in magnitude. Therefore it is necessary to weed out the frequencies above 1 Hz in the distinguishible case. Moreover, oversampling increases the maximum frequency that can be represented in the Fourier space and does not improve the precision of the estimator.

On a random dataset, the Fourier estimator suffers greatly from missing values. Figure 8 shows the distribution of τ F with 200 simulations and a dataset of size 15 where N is distributed according to a Poisson random variable with mean 5.5. The precision of the estimator is much worse than the compatible values estimator (see the plotted interval). The Fourier estimate τ F is compatible with the dataset in only about 1% of the simulations. Fig. 9 The linear regression estimator on a dataset without missing values.

Linear Regression Estimator

The observation model X = ⌊τN⌋ may be written X = τN + ε where ε is an error term. Even if ε is not Gaussian, linear regression can yield a reasonable estimate of the regression coefficient τ as Figure 9 shows.

We use the following algorithm:

-compute the empirical lattice {x i } by sorting and removing duplicates in the dataset -compute the indexes {n i } according to the sorting index -fit a regression line of the form x i = an i + 0.5 -return a Figure 9 shows the linear regression estimator on the dataset S = ⌊τ 1, 10 ⌋ (no missing values). For each element in the dataset, if the regression line intersects the length 1 interval then the estimate is compatible with the data point. Note that the truncation error is not centered. Consequently, we compute the regression coefficient in the the model X + 0.5 = τN + ε. For the same reason, the regressors are below the regression line.

The main difficulty in the linear regression is that the values of the regressor variable N are unknown. In the distinguishible case, it is possible to reconstruct them when there are no missing values, i.e. S τ ∩ 0, n = S where n = max S . Otherwise, the regressors will be shifted and that affects strongly the estimate. Figure 10 shows such a case. The regressors inferred in the linear regression estimator and the regression line are shown in solid line. For comparison, the true regressors are displayed in dotted line. The compatible values estimator finds the true regressors and its regression line is shown in dotted line.

Conclusion

In the observation model X = ⌊τN⌋, the parameter τ can be reliably estimated independently from N. This allows the full recovery of the statistics of N prior to modeling. The structure of N may then be studied at length afterwards.

The estimator based on compatible values is optimal and reasonably quick. It is resistant to missing values in practice, and in the worst case returns an acceptable (parcimonious) answer without hypotheses on the law of N.

Unfortunately, this estimator only takes into account truncation noise, and yields poor results on real data. We are currently pursuing an extension of the model that mixes electronic noise and truncation effects.

Compared to the other three estimators, the compatible values estimator performs much better but also more slowly. The Double Poisson Family is simply not a suitable model in our range of parameters, but there is room for improvement for the other estimators. For example, the main difficulty in the linear regression is computing the indexes. With some knowledge about the law of N, quantile regression could be applied.

The Fourier estimator suggests a strong relationship with the harmonic retrieval problem, although the signal is not periodic. Although the truncation error considered in this paper is very different from Gaussian errors usually considered in harmonic retrieval, some algorithms from that field may make a better compromise between speed and precision for the current problem.

Proposition 3 Lemma 1

 31 The set of compatible values C(S ) contains exactly ]0, 1] and intervals of length at least 1/ x2 where x = max S . In particular, if there are no outliers or missing values in S then τ belongs to the interval [a, b[ such that b = sup C(S ). The proof is based on the following two lemmas. The set of compatible values contains exactly ]0, 1] and a finite number of intervals of the form [a, b[ of length at least 1/ x2 where x = max S .

  There are only finitely many intervals of length at least 1/ x2 in [0, x+1 n max ], so there exists such an interval [a, b[. Let N denote the set N = {n|⌊nτ⌋ ∈ S }, i.e. the set of values for N that generate S . τ belongs to a certain interval [a ′ , b ′ [ which is the intersection of the constraints τ ∈ [ x n , x+1 n [, for all x = ⌊nτ⌋ in S . We show that b ′ = b, i.e. no positive real is both greater than b ′ and compatible. Let t such that b

Proposition 5

 5 If there are no missing values, the largest compatible value is found after at most x2 n ≃ τ x steps. With a small number of missing values k ≪ n, the number of steps is on the order of τ 2 x k + 1 τ where k = n -cardS is the number of missing values. Proof The procedure begins at B = x+1 n , ends before x n because τ ≥ a ≥ x n , and proceeds in steps of length 1/ x2 . Consequently, there are at most C = x2 ( x+1 n -x n ) steps. Let k = n -cardS denote the number of missing values. We make the following three approximations: k ≪ n, 1 ≪ x and τ ≃ x n . Then C = x2 ( x + 1)

Figure 1

 1 Figure 1 illustrates the compatible values estimator on a simulated dataset. The dataset {6, 6, 11, 5, 3, 5, 2, 6, 5, 13, 2, 7, 7, 7, 6} is obtained from the observation model X = ⌊1.32 * N⌋ where N is distributed according to a Poisson random variable with mean 5.5. It is first reduced to the lattice S = {2, 3, 5, 6, 7, 11, 13} and is shown at the bottom.

Fig. 1

 1 Fig. 1 Comparison of a few compatible lattices and the dataset.

2 Fig. 2

 22 Fig.2Comparison of the dataset and a few lattices that are not compatible.

Fig. 3

 3 Fig. 3 Length of the maximal interval for several values of τ with N ∈ 1,10 .

Fig. 4

 4 Fig. 4 Kernel density estimate of the distribution of the compatible values estimator on a random dataset.

Fig. 5

 5 Fig. 5 Kernel density estimate of the distribution of the double poisson estimate of τ with flooring noise (solid line) and without (dotted line).

Fig. 6

 6 Fig. 6 Kernel density estimate of the distribution of the double poisson estimate of τ with flooring noise (solid line) and without (dotted line).

Fig. 7

 7 Fig.[START_REF] Krutchinsky | On the nature of the chemical noise in MALDI mass spectra[END_REF] Fourier transform of the quasi-periodic signal, in Fourier space (left) and period space (right). The vertical line shows the upper bound from Proposition 9.

Figure 7

 7 Figure 7 shows in the frequency and period space the Fourier transform of the quasiperiodic signal obtained from the dataset N ∈ 1, 10 . The vertical line corresponds to the upper bound from Proposition 9.

Fig. 8

 8 Fig. 8 Kernel density estimate of the distribution of the Fourier estimator.

Fig. 10

 10 Fig.[START_REF] Windig | A noise and background reduction method for component detection in liquid chromatography/mass spectrometry[END_REF] The linear regression estimator in the case of missing values. The compatible values estimator is shown in dotted line.

Proposition 1 (Equivalence between lattices and numbers)

  In the distinguishible case, let t 1 and t 2 denote two real numbers such that t 1 ≥ 1 and t 2 ≥ 1. Then S t 1 = S t 2 if and only if t 1 = t 2 . Proof Obviously, if t 1 = t 2 then S t 1 = S t 2 . Let us prove the converse, i.e. S t 1 = S t 2 implies t 1 = t 2 or equivalently if t 1 = t 2 then S t 1 = S t 2 . Suppose that t 1 < t 2 . There exists n ∈ N such that ⌊t 1 n⌋ < ⌊t 2 n⌋. Either ⌊t 2 n⌋ / ∈ S t 1 , in which case S t 1 = S t 2 , or ⌊t 2

Single precision floating point numbers have

24 ∼ 16.10 6 levels of precision plus sign and exponent.

Appendix

Proposition [START_REF] Hayes | Statistical Digital Signal Processing and Modeling[END_REF] The mapping x → ⌊τ x⌋ is injective if and only if τ ≥ 1.

Proof If τ = 1, the mapping is the identity function. Suppose τ > 1 and let n 1 and n 2 denote two (positive) integers such that n 1 < n 2 . Then τn 2 -τn 1 > τ > 1 and ⌊τn 2 ⌋ > ⌊τn 1 ⌋. When τ < 1, the mapping is not injective because ⌊τ × 1⌋ = ⌊τ × 0⌋ = 0.

Upper Bounds on τ Proposition 7 (Any two observations) Let x and y be two distinct elements of the set S of observed values. Then τ < 1 + |x -y|.

Proof Let i and j be the values of N corresponding to x and y i.e. x = ⌊τi⌋ and y = ⌊τ j⌋. Then we have the inequalities: x ≤ τi < x + 1, y ≤ τ j < y + 1, and thus τ( ji) < yx + 1. Assuming x < y, we obtain τ < y-x+1 j-i < yx + 1.

Proposition 8 (Observed intervals) Let x, y denote the set of integers between x and y. If x, y is a subset of S , then τ < 1 + 1 yx .

Proof Using the same notations as in the proof of Proposition 7, τ < y-x+1 j-i < y-x j-i +

Proposition 9 (Density Upper Bound) Let x = ⌊τ n⌋ denote the largest integer in S . Then τ <

x + 1 n . When n is unknown (because of potential missing values), let n denote the number of non zero observed integers. Then τ <

x + 1 n ≤

x + 1 n .