AN AVERAGED PROJECTED ROBBINS-MONRO ALGORITHM FOR
ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL
DISTRIBUTION

arXiv:1606.04276v1 [math.ST] 14 Jun 2016

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

Abstract. The objective of this work is to propose a new algorithm to fit a sphere
on a noisy 3D point cloud distributed around a complete or a truncated sphere. More
precisely, we introduce a projected Robbins-Monro algorithm and its averaged version
for estimating the center and the radius of the sphere. We give asymptotic results such
as the almost sure convergence of these algorithms as well as the asymptotic normality
of the averaged algorithm. Furthermore, some non-asymptotic results will be given,
such as the rates of convergence in quadratic mean. Some numerical experiments
show the efficiency of the proposed algorithm on simulated data for small to moderate
sample sizes.

Keywords. Projected Robbins-Monro algorithm, Sphere fitting, Averaging, Asymptotic
properties.
1. Introduction
Primitive shape extraction from data is a recurrent problem in many research fields
such as archeology [Thom, 1955], medicine [Zhang et al., 2007], mobile robotics [Martins et al., 2008],
motion capture [Shafiq et al., 2001] and computer vision [Rabbani, 2006, Liu and Wu, 2014].
This process is of primary importance since it provides a high level information on the
data structure.
First works focused on the case of 2D shapes (lines, circles), but recent technologies
enable to work with three dimensional data. For instance, in computer vision, depth
sensors provide 3D point clouds representing the scene in addition to usual color images.
In this work, we are interested in the estimation of the center µ ∈ R3 and the radius
r > 0 of a sphere from a set of 3D noisy data. In practical applications, only a discrete
set of noisy measurements is available. Moreover, sample points are usually located only
near a portion of the spherical surface. Two kinds of problem can be distinguished:
shape detection and shape fitting.
Shape detection consists in finding a given shape in the whole data without any
prior knowledge on which observations belong to it. In that case, the data set may
represent several objects of different nature and may therefore contain a high number of outliers. Two main methods are used in practise to solve this problem. The
Hough transform [Abuzaina et al., 2013] performs a discretization of the parameter
space. Each observation is associated to a set of parameters corresponding to all
1

2

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

possible shapes that could explain the sample point. Then, a voting strategy is applied to select the parameter vectors of the detected shapes. The advantage of this
method is that several instances of the shape can be detected. However, a large
amount of memory is required to discretize the parameter space, especially in the
case of three dimensional models. The RANSAC (RANdom SAmple Consensus) paradigm [Fischler and Bolles, 1981, Schnabel et al., 2007] is a probabilistic method based
on random sampling. Observations are randomly selected among the whole data set and
candidate models are generated. Then, shapes can be detected thanks to an individual
scoring scheme. The success of the method depends on a given probability related to
the number of sampling and the fraction of points belonging to the shape.
The shape fitting problem assumes that all the data points belong to the shape.
For example, spherical fitting techniques have been used in several domains such as
industrial inspection [Jiang and Jiang, 1998], GPS localization [Beck and Pan, 2012],
robotics [Von Hundelshausen et al., 2005] and 3D modelling [Tran et al., 2015]. Geometric and algebric methods have been proposed [Landau, 1987, Rusu et al., 2003,
Al-Sharadqah, 2014] for parameters estimation. Moreover, let us note that fitting methods are generally applied for shape detection as a post-processing step in order to refine
the parameters of the detected shapes [Tran et al., 2015].
In a recent paper, Brazey and Portier [Brazey and Portier, 2014] introduced a new
spherical probability density function belonging to the family of elliptical distributions,
and designed to model points spread near a spherical surface. This probability density
function depends on three parameters, namely a center µ ∈ R3 , a radius r > 0 and a dispersion parameter σ > 0. In their paper, the model is formulated in a general form in Rd .
To estimate µ and r, a backfitting algorithm (see e.g. [Breiman and Friedman, 1985])
similar to the one used in [Landau, 1987] is employed. A convergence result is given
in the case of the complete sphere. However, no result is established in the case of a
truncated sphere while simulations showed the efficiency of the algorithm.
The objective of this work is to propose a new algorithm to fit a sphere on a noisy 3D
point cloud distributed around a complete or a truncated sphere. We shall assume that
the observations are independent realizations of a random vector X defined as

(1.1)

X = µ + r W UΩ ,

where W is a positive real random variable such that E [W ] = 1, UΩ is uniformly distributed on a measurable subset Ω of the unit sphere of R3 , W and UΩ are independent.
Parameters µ ∈ R3 and r > 0 are respectively the center and the radius of the sphere
we are trying to adjust to the point cloud. Random variable W allows to model the
fluctuations of points in the normal direction of the sphere. When Ω coincides with the
complete sphere, then the distribution of X is spherical (see e.g. [Muirhead, 2009]).
Indeed, if we set Y = (X − µ)/r, then the distribution of Y is rotationally invariant.

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

3

We are interested in estimating center µ and radius r. As kUΩ k = 1, we easily deduce
from (1.1) that


(X − µ)
(1.2)
µ = E X −r
kX − µk
(1.3)
r = E [kX − µk] .
It is clear that from these two equations, we cannot deduce explicit estimators of parameters µ and r using the method of moments since each parameter depends on
the other. To overcome this problem, we can use a backfitting type algorithm (as in
[Brazey and Portier, 2014]) or introduce a recursive stochastic algorithm, which seems
well-suited for this problem since equations (1.2) and (1.3) can also be derived from the
local minimization of the following quadratic criteria

1 
G(µ, r) := E (kX − µk − r)2 .
(1.4)
2
Stochastic algorithms, and more precisely Robbins-Monro algorithms, are effective and
fast methods (see e.g. [Duflo, 1997, Kushner and Yin, 2003, Robbins and Monro, 1951]).
They do not need too much computational efforts and can easily be updated, which make
of them good candidates to deal with big data for example. However, usual sufficient
conditions to prove the convergence of this kind of algorithm are sometimes not satisfied
and it is necessary to modify the basic algorithm. We can, for example, introduce a
projected version of the Robbins-Monro algorithm which consists in keeping the usual
estimators in a nice subspace with the help of a projection. Such an algorithm has been
recently considered in [Bercu and Fraysse, 2012] and [Lacoste-Julien et al., 2012].
In this paper, due to the non global convexity of function G, we estimate parameters
µ and r using a projected Robbins-Monro algorithm. We also propose an averaged
algorithm which consists in averaging the projected algorithm. In general, this averaged
algorithm allows to improve the rate of convergence of the basic estimators, or to reduce
the variance, or not to have to make a good choice of the step sequence, which can be
as exhaustive as to estimate the parameters. It is widely used when having to deal with
Robbins-Monro algorithms (see [Polyak and Juditsky, 1992] or [Pelletier, 1998] amoung
others).
This paper is organized as follows. In Section 2, we specify the framework and assumptions. After a short explanation on the non-convergence of the Robbins-Monro
algorithm, the projected algorithm and its averaged version are introduced in Section 3.
Section 4 is concerned with the convergence results. Some simulation experiments are
provided in Section 5, showing the efficiency of the algorithms. Proofs of the different
results are postponed in Appendix.
2. Framework and assumptions
We consider in this paper a more general framework than the one described in the
introduction. Let X be a random vector of Rd with d ≥ 2. Let F denotes the distribution
of X. We assume that X can be decomposed under the form
(2.1)

X = µ + r W UΩ .

4

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

where µ ∈ Rd , r > 0, W is a positive real continuous random variable (and with a
bounded density if d = 2), UΩ is uniformly distributed on a measurable subset Ω of the
unit sphere of Rd . Moreover, let us suppose that W and UΩ are independent.
Model (2.1) allows to model a point cloud of Rd spread around a complete or truncated
sphere of center µ ∈ Rd and radius r > 0. Random vector UΩ defines the position of
the points on the sphere and random variable W defines the fluctuations in the normal
direction of the sphere. As mentioned in the introduction, when Ω is the complete unit
sphere, then the distribution of X is spherical.
When W satisfies the condition E [W ] = 1, the radius r is identifiable and can be directly estimated. Indeed, since kUΩ k = 1, then kX − µk = rW and E [kX − µk] = r E [W ] = r.
However, this condition is sometimes not satisfied (as in [Brazey and Portier, 2014]) and
only r? := r E [W ] can be estimated. Therefore, in what follows, we are interested in
T
estimating θ := µT , r? , which will be denoted by (µ, r? ) for the sake of simplicity.
We suppose from now that the following assumptions are fulfilled:
• Assumption [A1]. The random vector X is not concentrated around µ:


E kX − µk−2 < ∞.

• Assumption [A2]. The random vector X admits a second moment:


E kX − µk2 < ∞.
These assumptions ensure that the values of X are concentrated around the sphere
and not around the center µ, without in addition too much dispersion. This framework
totally corresponds to the real situation that we want to model. Moreover, using (2.1),
Assumptions [A1] and [A2] reduce to assumptions on W . More precisely, [A1] reduces
to E [W −2 ] < ∞ and [A2] to E [W 2 ] < ∞.
Let us now introduce two examples of distribution allowing to model points spread
around a complete sphere and satisfying assumptions [A1] and [A2].
Example 2.1. Let us consider a random vector X taking values in Rd with a distribution
absolutely continuous with respect to the Lebesgue measure, with a probability density
function fδ defined for all δ > 0 by
(2.2)

fδ (x) =

Cd

1
,
kx − µkd−1 {kx − µk /r ∈ [1 − δ , 1 + δ]}

where Cd is the normalization constant. Then, we can rewrite X under the form (2.1)
with UΩ = U , W ∼ U([1 − δ, 1 + δ]) and E [W ] = 1 for any δ > 0.
Example 2.2. Let us consider the probability density function introduced in [Brazey and Portier, 2014].
It is defined for any x ∈ Rd by


2
1
(2.3)
f (x) = Kd exp − 2 kx − µk − r
,
2σ

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

5

where Kd is the normalization constant. Then, a random vector X with probability
density function f can be rewritten under the form (2.1), with E [W ] 6= 1, but E [W ] is
closed to 1 when the variance σ is negligible compared to the radius r.
To obtain points distributed around a truncated sphere, it is sufficient to modify the
previous densities by considering densities of the form fΩ (x) = CΩ f (x)1
(x − µ) ∈ Ω
d
where Ω is the set of points of R whose polar coordinates are given by (ρ, θ1 , . . . , θd−1 ∈
R∗+ × Θ) where Θ defines the convex part Ω of the surface of the unit sphere of Rd we
want to consider.
3. The algorithms
We present in this section two algorithms for estimating the unknown parameter θ
which can be seen as a local minimizer (under conditions) of a function. Indeed, let us
consider the function G : Rd × R −→ R defined for all y = (z, a) ∈ Rd × R∗+ by

1
1 
(3.1)
G(y) :=
E (kX − zk − a)2 = E [g (X, y)] ,
2
2
where we denote by g the function defined for any x ∈ Rd and y = (z, a) ∈ Rd × R∗+
by g(x, y) := (kx − zk − a)2 . The function G is Frechet-differentiable and we denote
by Φ its gradient, which is defined for all y = (z, a) ∈ Rd × R∗+ by



z−X
z − E [X] − a E
(3.2) Φ(y) := ∇G(y) = E [∇y g(X, y)] = 
kz − Xk 
a − E [kz − Xk] .
From (2.1) and definition of θ = (µ, r? ), we easily verify that ∇G(θ) = 0. Therefore,
since θ is a local minimizer of function G (under assumptions) or a zero of ∇G, an idea
could be to introduce a stochastic gradient algorithm for estimating θ.
3.1. The Robbins-Monro algorithm.
Let (Xn )n≥1 be a sequence of independent and identically distributed random vectors
of Rd following the same law as X and let (γn )n≥1 be a decreasing sequence of positive
real numbers satisfying the usual conditions
X
X
(3.3)
γn = ∞ and
γn2 < ∞.
n≥1

n≥1

When the functional G is convex or verifies nice properties, a usual way to estimate the
unknown parameter θ is to use the following recursive algorithm
(3.4)

θn+1 = θn − γn ∇y g(Xn+1 , θn ),

with θ1 chosen arbitrarily bounded. The term ∇y g (Xn+1 , θn ) can be seen as an estimate
of the gradient of G at θn , and the step sequence (γn ) controls the convergence of the
algorithm.
The convergence of such an algorithm is often established using the Robbins-Siegmund’s
theorem (see e.g. [Duflo, 1997]) and a sufficient condition to get it, is to verify that for
any y ∈ Rd × R∗+ , hΦ(y), y − θi > 0 where h., .i denotes the usual inner product and

6

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

k.k the associated norm. However, we can show that this condition is only satisfied for
y belonging to a subset of Rd × R∗+ to be specified. Thus, if at time (n + 1), the update
of θn (using (3.4)) leaves this subset, then it does not necessarily converge. Therefore,
we have to introduce a projected Robbins-Monro algorithm.
3.2. The Projected Robbins-Monro algorithm.
Let K be a compact and convex subset of Rd × R∗+ containing θ = (µ, r∗ ) and let
π : Rd × R∗+ −→ K be a projection satisfying

∀y, y 0 ∈ Rd × R∗+ , kπ(y) − π(y 0 )k ≤ ky − y 0 k
(3.5)
∀y ∈
/ K, π(y) ∈ ∂K
where ∂K is the frontier of K. An example will be given later.
Then, we estimate θ using the following Projected Robbins-Monro algorithm (PRM),
defined recursively by



(3.6)
θbn+1 = π θbn − γn ∇y g Xn+1 , θbn ,
where θb1 is arbitrarily chosen in K, and (γn ) is a decreasing sequence of positive real
numbers satisfying (3.3).
Of course the choice of subset K and projector π is crucial. It is clear that if K is
poorly chosen for a given projector, the convergence of the projected algorithm towards
θ will be slower, even if from a theoretical point of view, we shall see in the next
section dedicated to the theoretical results, that this algorithm is almost the same as
the traditional Robbins-Monro
algorithm since the updates of θbn , ie. the quantities

θbn − γn ∇y g Xn+1 , θbn , leave K only a finite number of times.
Let us now discuss the choice of K and π. The choice of K is directly related to the
following assumption that we introduce to ensure the existence of a compact subset on
which the scalar product hΦ(y), y − θi is positive.
• Assumption [A3]. There are two positive constants Rµ and Rr such that for
all y = (z, a) ∈ B(µ, Rµ ) × B(r∗ , Rr ),
(3.7)

 1 − kE [UΩ ]k2 /A
λmax Γ(z) <
,
r∗ + 32 Rr
z∈B(µ,Rµ )
sup

with A such that kE [UΩ ]k2 < A < 1, and λmax (M ) denotes the largest eigenvalue of matrix M , and



1
(X − z)(X − z)T
Γ(z) := E
Id −
.
kX − zk
kX − zk2
Remark 3.1. The less the sphere is troncated, the more kE [UΩ ]k is close to 0 and the
constraints on Rµ and Rr are relaxed. In particular, when the sphere is complete,
ie. UΩ = U where U denotes the random vector uniformly distributed on the whole

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

7

unit sphere of Rd , then E [UΩ ] = 0 and Assumption [A3] reduces to
sup


λmax Γ(z) <

z∈B(µ,Rµ )

r∗

1
.
+ 23 Rr

The main consequence of Assumption [A3] is the following proposition which is one
of the key point to establish the convergence of the PRM algorithm.
Proposition 3.1. Assume that [A1] to [A3] hold. Then, there is a positive constant c
such that for all y ∈ B(µ, Rµ ) × B(r∗ , Rr ),
hΦ(y), y − θi ≥ c ky − θk2 .
Proof. The proof is given in Appendix A.



Assumption [A3] is therefore crucial but only technical. It reflects the fact that the
sphere is not too much truncated and that the points are not too far away from the
sphere which corresponds to the real situations we want to model.
In a general framework, this technical assumption is difficult to verify since it requires to
specify the distribution of X. In the case of distribution of Example 2.1 with δ < 1/10,
we can easily exhibit constant Rµ and Rr . Indeed taking Rµ = Rr = r∗ /10, then
assumption [A3] holds. When the distribution of X is compactly supported with a
support included in [1 − δ, 1 + δ], it is fairly easy to find the constants provided that δ
is small enough. It is quite more difficult when dealing with distribution of Example 2.2.
Nevertheless, topological results can ensure that these constants exist.
From constants Rµ and Rr of Assumption [A3], it is then possible to simply define
a projector π which satisfies condition (3.5). Indeed, let us set K = Kµ × Kr with
Kµ = B(µ, Rµ ) and Kr = B(r∗ , Rr ), and define for any y = (z, a) ∈ Rd × R∗+ by
π(y) := (πµ (z), πr (a)), with

z
if z ∈ Kµ

(z − µ)
πµ (z) :=
otherwise
 µ + Rµ
kz − µk
and



a
if a ∈ Kr0
?
(a − r )
πr (a) :=
otherwise
 r + Rr
|a − r? |
Such projector satisfies the requested conditions. However, it is clear that this projector
can not be implemented since µ and r∗ are unknown. We shall see in the simulation
study how to overcome this problem.
We suppose from now that K is a compact and convex subset of B(µ, Rµ )×B(r∗ , Rr )
such that θ ∈ K, but θ ∈
/ ∂K, where ∂K is the frontier of K, i.e there is a positive
constant dmin such that B (θ, dmin ) ⊂ K.

8

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

3.3. The averaged algorithm.
Averaging is a usual method to improve the rate of convergence of Robbins-Monro
algorithms, or to reduce the variance, or finally not to have to make a good choice of the
step sequence (see [Polyak and Juditsky, 1992]), but for the projected algorithms, this
method is not widespread in the litterature. In this paper, we improve the estimation
of θ by adding an averaging step to the PRM algorithm. Starting from the sequence
(θbn )n≥1 given by (3.6), we introduce for any n ≥ 1,
n

1 Xb
θn =
θk ,
n k=1
which can also be recursively defined by

1 b
(3.8)
θn+1 = θn +
θn+1 − θn , and θ1 = θb1 .
n+1
We shall see in the following two sections, the gain provided by this algorithm.
4. Convergence properties
We now give asymptotic properties of the algorithms. All the proofs are postponed in
Appendix B. The following theorem gives the strong consistency of the PRM algorithm
as well as properties on the number of times we really use the projection.
Theorem 4.1. Let (Xn ) be a sequence of iid random vectors following the same law as
X. Assume that [A1] to [A3] hold, then
lim kθbn − θk = 0 a.s.

n→∞



Moreover, the number of times the random vectors θbn − γn ∇y g Xn+1 , θbn do not
belong to K is almost surely finite.
The following theorem gives the rate of convergence in quadratic mean and the
Lp rates of convergence of the PRM algorithm (under conditions)
 as wellas an upper bound of the probability that the random vector θbn − γn ∇y g Xn+1 , θbn does not
belong to K.
Theorem 4.2. Let (Xn ) be a sequence of iid random vectors following the same law
as X. Assume that [A1] to [A3] hold and consider a step sequence (γn ) of the form
γn = cγ n−α , with cγ > 0 and α ∈]1/2, 1[. Then, there is a positive constant C1 such
that for all n ≥ 1,


2
C1
b
E θn − θ
≤ α.
n


Moreover, for all positive integer p such that E kX − µk2p < ∞, there is a positive
constant Cp such that for all n ≥ 1,


2p
Cp
E θbn − θ
≤ pα ,
n

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

9

and for all n ≥ 1,
h


i
P θbn − γn ∇y g Xn+1 , θbn ∈
/K ≤

Cp
,
npα

d2p
min

where dmin := inf y∈∂K {ky − θk} and ∂K is the frontier of K.
We now focus on the asymptotic behavior of the averaged algorithm. First of all,
applying Theorem 4.1 and Toeplitz’s lemma for example, we easily obtain the strong
consistency of the averaged estimator θn . Introducing the following assumption, we can
specify its rate of convergence in quadratic mean as well as its asymptotic normality.
• Assumption [A4]. The Hessian of G at θ = (µ, r∗ ), denoted by Γθ and defined
by






X −µ
1
(X − µ) ⊗ (X − µ)
∗
E
Id − r E kX − µk Id −
kX − µk 
kX − µk2




Γθ := 
T

X −µ
E
1
kX − µk
is a positive definite matrix.
Note that thanks to topological results, this assumption also implies Proposition 3.1
but is not useful to obtain the constants Rµ and Rr . Nevertheless, this assumption is
crucial to establish the results of the two following theorems but it is satisfied as soon
as the sphere is not too much truncated and the dispersion around the sphere not too
important which corresponds to the real situations encountered. Using model (2.1), Γθ
rewrites under the form






(Id − β Id − E UΩ UΩT
E [UΩ ]
(4.1)
Γθ =
with β = E [W ] E W −1 .
T
E UΩ
1


When the sphere is complete, ie. UΩ = U , then E [UΩ ] = 0, E UΩ UΩT = (1/d)Id and
λmin (Γθ ) > 0 as soon as β < d/(d − 1). In the case of distribution of Example 2.1, we
have β = (log(1+δ)−log(1−δ))/(2δ) and [A4] is satisfied as soon as δ is small enough.
In the case of distribution of Example 2.2, [A4] is satisfied as soon as r >> σ. When
the sphere is not complete, we can easily show that a sufficient condition to ensure [A4]
is λmin (Var [UΩ ]) < 1 − 1/β, where Var [UΩ ] is the covariance matrix of the random
variable UΩ . In the case of the half sphere and d = 3, we have λmin (Var [UΩ ]) = 1/12
and Γθ is definite positive as soon as β < 12/11. This condition holds for distribution
of Example 2.1 with δ < 0.4 for instance, and distribution of Example 2.2 as soon as
r >> σ.
Theorem 4.3. Let (Xn ) be a sequence of iid random vectors following the same law
as X. Assume that [A1] to [A4] hold and consider a step sequence (γn ) of the form
γn = cγ n−α , with cγ > 0 and α ∈]1/2, 1[. Moreover, suppose that E[kX − µk12 ] < ∞.
Then there is a positive constant C such that for all n ≥ 1,
i C
h
2
≤ .
E θn − θ
n

10

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

With respect to results of Theorem 4.2, we clearly improve the rate of convergence
in quadratic mean. Note that the computed rate is the optimal one for such stochastic
algorithms. We finally give a central limit theorem which can be useful to build confidence
balls for the different parameters of the sphere.
Theorem 4.4. Let (Xn ) be a sequence of iid random vectors following the same law
as X and let us choose the step sequence (γn ) of the form γn = cγ n−α , with cγ > 0
12
and α ∈]1/2,
 1[. Assume that [A1] to [A4] hold and suppose that E[kX − µk ] < ∞.
Then θn satisfies
 L

√
−1
n θn − θ −→ N 0, Γ−1
(4.2)
Σ
Γ
θ
θ
n→∞

with

T 
(µ − X)
∗ (µ − X)
µ−X −r
µ−X −r
Σ := E 
kµ − Xk  
kµ − Xk   .
∗
∗
r − kµ − Xk
r − kµ − Xk


(4.3)

∗

From result (4.2) of Theorem 4.4, we easily derive that
 L
√ −1/2
(4.4)
nΣ
Γθ θn − θ −→ N (0, Id+1 ) .
n→∞

Therefore, in order to build confidence balls or statistical tests for the parameters of the
sphere, matrices Γθ and Σ must be estimated.
Let us decompose θn under the form (Z n , An ) where Z n ∈ Rd estimates the center µ
and An ∈ R∗+ the radius r∗ , and let us denote Un := (Xn − Z n )/ Xn − Z n . Then we
bn and Σ
b n iteratively as follows
can estimate Γθ and Σ by Γ


!
An
An
T
1
−
I
+
U
U
U
d
n
bn = (n − 1)Γ
bn−1 +
kXn −Z n k
kXn −Z n k n n
nΓ
,
T
1
Un


T
Z
+
A
U
Z
+
A
U
X
−
X
−
n
n
n
n
n
n
n
n
b n = (n − 1)Σ
b n−1 +
nΣ
,
An − kXn − Z n k
An − kXn − Z n k
b 1 = Id+1 and Γ
b1 = Id+1 to avoid usual problems of invertibility. It is not hard
where Σ
bn and Σ
b n respectively converge to Γθ and Σ and then deduce that
to show that Γ
 L
√ −1/2
b
bn θn − θ −→
(4.5)
Qn := n Σ
Γ
N (0, Id+1 ) .
n
n→∞

The simulation study of the next section will illustrate the good approximation of the
distribution of Qn by the standard gaussian for moderate sample sizes.
5. Some experiments on simulated data
We study in this section the behavior of the PRM and averaged algorithms on simulated data in the case d = 3, for small to moderate sample sizes. This section first
begins with the specification of the compact set involved in the definition of the PRM
algorithm which is of course a crucial point. We then study the performance of the
two algorithms in the case of the whole sphere with the distributions of Examples 2.1

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

11

and 2.2. Finally, we consider the case of the truncated sphere (a half-sphere) and we
compare our strategy with the one proposed by [Brazey and Portier, 2014].
In this simulation study, we shall always consider the same sphere defined by its center
µ = (0, 0, 0)T and its radius r = 50. In addition, to reduce sampling effects, our results
are based on 200 samples of size n. Finally, let us mention that simulations were carried
out using the statistical software R (see R Core Team, 2013).
5.1. Choice of the compact set and of the projection. We discuss here the crucial
point of the choice of the compact set K and of the projection π involved in the definition
of the PRM algorithm. The main problem is to find a compact set containing the
unknown parameter θ. We propose to build a preliminary estimation of θ, using a
geometric approach which consists in finding the center and the radius of a sphere of R3
from 4 non-coplanar distinct points. We denote by (µ0 , r0 ) this initial estimate of θ. From
this estimate, we define the compact set K by K := Kµ0 ×Kr0 with Kµ0 := B(µ0 , r0 /10)
and Kr0 := B(r0 , r0 /10), where the choice of the value r0 /10 for the radius of the balls
is justified by the discussion about Assumption [A3] in Section 3.2. We then define the
projector π as follows: for any y = (z, a) ∈ R3 × R∗+ , we set π(y) := (πµ0 (z), πr0 (a))
with

z
if z ∈ Kµ0

r0 (z − µ0 )
πµ0 (z) :=
otherwise
 µ0 +
10 kz − µ0 k
and

a
if a ∈ Kr0

r0 (a − r0 )
πr0 (a) :=
otherwise
 r0 +
10 |a − r0 |
With this strategy, we can raisonnably hope that if our initial estimate is not too poor,
then the true parameter belongs to K and the quadratic criteria G is convex on K.
We will see below that even if this preliminary estimation is rough, the true parameter
belongs to K and the PRM algorithm improves the estimation of θ.
Let us now describe our strategy to obtain a preliminary estimation of the parameter
θ = (µ, r? ). Since the data points are spread around the sphere, the estimation of the
parameters from only one quadruplet of points is not robust to random fluctuations. In
order to make the estimation more robust, we consider instead N quadruplets sampled
with replacement from the first K points of the sample X1 , . . . , Xn . For each quadruplet,
we calculate the center of the sphere which passes through these four points, which gives
a sequence of centers (b
µi )1≤i≤N . The initial estimate of the center, denoted by µ0 , is
then computed as the median point. Finally, we obtain an estimation of the radius by
calculating the empirical mean of the sequence (kXi − µ0 k)1≤i≤50 .
A simulation study carried out for various values of K and N in the case of the
whole and truncated sphere, shows that by taking K = 50 and N = 200, we obtain a
preliminary estimation of θ sufficiently good to ensure that the compact K contains θ.
To close this section, let us mention that although the initial estimate is quite accurate,
it is necessary to project the Robbins-Monro algorithm to ensure the convergence of the

12

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

estimator. Indeed, taking a step sequence of the form γn = cγ n−α , the results given
in Table 1 show that for some values of cγ and α, the parameter θ is poorly estimated
by the Robbins-Monro algorithm, while the PRM algorithm (Table 2) is less sensitive to
the step sequence choice.

cγ

α
0.51 0.6 0.66 0.75 0.99
1 0.27 0.14 0.09 0.06 0.23
5 108 106 105 104 105
10 1031 1018 1014 1010 106

Table 1. Robbins-Monro algorithm. Errors in quadratic mean of the
200 estimations of the center µ for samples of size n = 2000 in the case
of the distribution of Example 2.1.

cγ

α
0.51 0.6 0.66 0.75
1 0.28 0.15 0.09 0.05
5 1.55 0.76 0.48 0.24
10 3.22 1.35 0.94 0.43

0.99
0.24
0.05
0.08

Table 2. PRM algorithm. Errors in quadratic mean of the 200 estimations of the center µ for samples of size n = 2000 in the case of the
distribution of Example 2.1.

In the sequel of the simulation study, we take a step sequence of the form γn := n−2/3
(α = 2/3 is often considered as the optimal choice in the literature).

5.2. Case of the whole sphere. In what follows, we are interested in the behavior of
the PRM and averaged algorithms when samples are distributed on the whole sphere
according to the distribution of Example 2.1 with δ = 0.1.
Figure 1 shows that the accuracy of the estimations increases with the sample size. In
particular, as expected, the PRM algorithm significantly improves the initial estimations
of the center and the radius (see the first boxplots which correspond to the initial
estimations). Moreover, as expected in the case of the ”whole sphere”, we can see that
the three components of the center µ are estimated with the same accuracy.

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

13

Figure 1. Whole sphere with distribution of Example 2.1. From the
left to the right, boxplots of estimates of µx , µy , µz and r obtained with
the PRM algorithm for different sample sizes.

52

52

Let us now examine the gain provided by the use of the averaged algorithm. Figure 2
shows that for small sample sizes, the performances of the two algorithms are comparable,
but when n is greater than 500, the averaged algorithm is more accurate than the PRM
algorithm. We can even think that by forgetting the first estimates of the PRM algorithm,
we improve the behavior of the averaged algoritm when the sample size is small.

●

●

51

51

●

●

●
●

●

●

●
●

●
●
●

50

50

●

●
●

●
●

49

●

48

48

49

●
●

●
●

●

200

200

500

500

1000 1000 1500 1500 2000 2000

200

200

500

500

1000 1000 1500 1500 2000 2000

Figure 2. Whole sphere with distribution of Example 2.1. Boxplots of
estimates of µy (left) and r (right) obtained with the PRM algorithm (in
red) and with the averaged algorithm (in blue) for different sample sizes.
Finally, let us study the quality of the Gaussian approximation of the distribution of
Qn for a moderate sample size. This point is crucial for building confidence intervals or
statistical tests for the parameters of the sphere.
Figure 3 shows that this approximation is reasonable when n = 2000. Indeed, we can
see that the estimated density of each component of Qn is well superimposed with the

14

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

density of the N (0, 1). To validate these approximations, we perform a KolmogorovSmirnov test at level 5%. The test enables us to conclude that the normality is not
rejected for each component of Qn .

Figure 3. From the left to the right, estimated densities of each
components of Q2000 superimposed with the standard gaussian density.

5.3. Comparison with a backfitting-type algorithm in the case of a half-sphere.
In this section, we compare the performances of the averaged algorithm with the ones
of the backfitting algorithm introduced by [Brazey and Portier, 2014]. In what follows,
we consider samples coming from the distribution of Example 2.2, with σ = 1, in the
case of the half sphere defined by the set of points whose y-component is positive.
Results obtained with the two algorithms are presented in Figure 4. We focus on
parameter µy for the center since it is the more difficult to estimate. We can see that
even if the backfitting (BF for short) algorithm is better than the averaged algorithm, the
performances are globally good, which validates the use of our algorithm for estimating
the parameters of a sphere from 3D-points distributed around a truncated sphere. Recall
that convergence results are available for our algorithm in the case of the truncated
sphere, contrary to the backfitting algorithm for which no theoretical result is available
in that case.

15

1.0

51.0

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

●

●
●
●

50.5

0.5

●
●

●

●

●

●

●
●

●

●
●
●

0.0

50.0

●

●
●

49.5

−0.5

●
●
●

49.0

−1.0

●

●

500

500

1000

1000

1500

1500

2000

2000

500

500

1000

1000

1500

1500

2000

2000

Figure 4. Comparison of averaged and BF algorithms. Boxplots of the
estimates of µy (on the left) and r (on the right), obtained with the BF
algorithm (in blue) and with the averaged algorithm (in red) for the half
sphere in the case of Example 2.2.

6. Conclusion
We presented in this work a new stochastic algorithm for estimating the center and
the radius of a sphere from a sample of points spread around the sphere, the points
being distributed around the complete sphere or only around a part of the sphere.
We shown on simulated data that this algorithm is efficient, less accurate than the
backfitting algorithm proposed in [Brazey and Portier, 2014] but for which no convergence result is available for the case of the truncated sphere. Therefore, our main
contribution is to have proposed an algorithm for which we have given asymptotic results such as its strong consistency and its asymptotic normality which can be useful to
build confidence balls or statistical tests for example, as well as non asymptotic results
such as the rates of convergence in quadratic mean.
A possible extension of this work could be to extend the obtained results to the case of
the finite mixture model. This framework has been considered in [Brazey and Portier, 2014]
but no convergence result is established. Proposing a stochastic algorithm for estimating
the different parameters of the model and obtaining convergence results would be a nice
challenge.
Acknowledgements. The authors would like to thank Aurlien Vasseur for his contribution to the start of this work. We also would like to thank Peggy Cénac and Denis
Brazey for their constructive remarks and for their careful reading of the manuscript that
allowed to improve the presentation of this work. Finally, we would like to thank Nicolas
Godichon for his help in the creation of Figure 5.

16

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

Appendix A. Some convexity results and proof of proposition 3.1
The following lemma ensures that the Matrix in Assumption [A3] is well defined and
that the Hessian of G exists for all y ∈ Rd × R.
Lemma A.1. Assume [A1] holds. If d ≥ 3, there is a positive constant C such that for
all z ∈ Rd ,


1
E
≤ C.
kX − zk
Moreover, suppose that W admits a bounded density, then for all d ≥ 2, there is a
positive constant C such that for all z ∈ Rd ,



1
E
≤ C.
kX − zk
Note that for the sake of simplicity, we denote by the same way the two constants.
Proof of Lemma A.1. Step 1: d ≥ 3
By continuity and applying Assumption [A1], there are positive constants , C 0 such
that for all z ∈ B (µ, ),


1
≤ C 0.
E
kX − zk
Moreover, let z ∈ Rd such that kz − µk ≥ , we have
 Z +∞ 

1
1
P kX − zk ≤
E
=
dt
kX − zk
t
0
Z M
Z ∞




−1
=
P kX − zk ≤ t
dt +
P kX − zk ≤ t−1 dt
0
M
Z ∞


P kX − zk ≤ t−1 dt,
≤M+


M

with M positive and defined later. Moreover, let t ≥ M ,




P kX − zk ≤ t−1 = P kµ + rW UΩ − zk ≤ t−1



≤ P −t−1 + kz − µk ≤ rW ≤ t−1 + kz − µk , (µ + rW UΩ ) ∩ B z, t−1 6= ∅ ,
taking M = 2 . With previous condition on rW , calculating P [(µ + rW UΩ ) ∩ B (z, t−1 ) 6= ∅]
consists in measuring the intersection between a truncated sphere with radius bigger
than /2 with a ball of radius 1t , with 1t ≤ 2 . This is smaller than the surface of the
frontier of the ball (see the following figure).

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

17

Figure 5.
Thus, there is a positive constant k such that for all t ≥ M ,


k
(A.1)
P kX − zk ≤ t−1 ≤ d−1 .
t
Finally,


Z +∞
1
2
1
E
≤ +
k d−1 dt
2
kX − zk

t

2
d−2
+ k d−2
.

2 (d − 2)
n
o
d−2
0 2
We conclude the proof taking C = max C ,  + k 2d−2 (d−2) .
=

Step 2: d = 2 and W admits a bounded density
Let fmax be a bound of the density function of W . As in previous case, let z ∈ Rd such
that kz − µk ≥ ,





P kX − zk ≤ t−1 ≤ P −t−1 + kz − µk ≤ rW ≤ t−1 + kz − µk , (µ + rW UΩ ) ∩ B z, t−1 6= ∅
h
i

−1
−1
−1
= P (µ + rW UΩ ) ∩ B z, t
6= ∅ − t + kz − µk ≤ rW ≤ t + kz − µk

 −1
× P −t + kz − µk ≤ rW ≤ t−1 + kz − µk .

18

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

As in previous case, if t ≥ 2 , there is a positive constant k such that for all t ≥ 2 ,
i
h

P (µ + rW UΩ ) ∩ B z, t−1 6= ∅ − t−1 + kz − µk ≤ rW ≤ t−1 + kz − µk ≤ kt−1 .
Moreover, since fmax is a bound of the density function of W ,


1
1
2rfmax
P − + kz − µk ≤ rW ≤ + kz − µk ≤
t
t
t
Thus, for all t ≥ 2 ,

 2rfmax k
P kX − zk ≤ t−1 ≤
,
t2
and in a particular case,

2
1
≤ + krfmax ,
E
kX − zk



(A.2)

and one can conclude the proof taking C = max {C 0 , 2−1 + krfmax }.



Proof of Proposition 3.1. We want to show there is c > 0 such that for any y = (z, a) ∈ B(µ, εµ ) × B(r∗ , εr ),
P (y) := hy − θ , Φ(y)i ≥ c ky − θk. We have


 +
*

z−X
z − µ z − E [X] − a E
(A.3)
P (y) = P (z, a) =
,
kz − Xk  .
a − r∗
a − E [kX − zk]
For any z ∈ Rd , let us set F (z) := E [kX − zk] and f (z) := E [(z − X)/ kz − Xk].
Note that f is the gradient of F . Using (2.1), we deduce that F (µ) = r∗ , f (µ) = − E [UΩ ]
and E [X] = µ − r? f (µ). Then, (A.3) can be rewritten as
P (y) = kz − µk2 + r∗ hz − µ, f (µ)i − ahz − µ, f (z)i + (a − r∗ )2 − (a − r∗ )(F (z) − F (µ))
= kz − µk2 − (a − r? ) hz − µ, f (µ)i − a hz − µ, f (z) − f (µ)i + (a − r? )2 − (a − r? )(F (z) − F (µ)).
Moreover, using the following Taylor’s expansions,
1
F (z) = F (µ) + hz − µ, f (µ)i + (z − µ)T ∇f (c)(z − µ),
2
f (z) = f (µ) + h∇f (c0 ), z − µi,
with c, c0 ∈ [z, µ]. We get
P (y) = kz − µk2 − 2(a − r? ) hz − µ, f (µ)i − a(z − µ)T ∇f (c0 )(z − µ)
1
+ (a − r∗ )2 − (a − r? )(z − µ)T ∇f (c)(z − µ)
2
Now, remarking that for any positive constant A and real numbers x, y, we have
2xy ≤ A x2 + y 2 /A, we derive
1
P (y) ≥ kz − µk2 − A(a − r∗ )2 − kz − µk2 kf (µ)k2 − a k∇f (c)kop kz − µk2
A
1
+ (a − r∗ )2 − |a − r? | k∇f (c0 )kop kz − µk2 .
2
(A.4)

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

Let us denote by λM =

19

λmax ∇f (z) and choose A such that kf (µ)k2 = kE [UΩ ]k2 < A < 1.

sup
z∈B(µ,εµ )

Then, for any z ∈ B(µ, εµ ) and a ∈ B(r∗ , εr ), we have


1
3
2
∗
P (y) ≥ 1 − kf (µ)k − (r + εr ) λM kz − µk2 + (1 − A) (a − r∗ )2
A
2
Finally, using assumption [A3], we close the proof.

In order to linearize the gradient in the decompositions of the PRM algorithm and get
a nice decomposition of the averaged algorithm, we introduce the Hessian matrix of G,
denoted, for all y = (z, a) ∈ Rd × R, by Γy : Rd × R −→ Rd × R and defined by :
(A.5)






1
(X − z) ⊗ (X − z)
X −z
E
Id − aE kX − zk Id −
kX − zk2
kX − zk 


,
Γy = 



T


X −z
E
1
kX − zk
with, for all z, z 0 , z 00 ∈ Rd , z ⊗ z 0 (z 00 ) = hz, z 00 iz 0 . Applying Lemma A.1, the Hessian
matrix exists for all y ∈ Rd+1 .
Proposition A.1. Suppose [A1] to [A3] hold, there is a positive constant Cθ such that
for all y ∈ K,
kΦ(y) − Γθ (y − θ)k ≤ Cθ ky − θk2 .
Proof of Proposition A.1. Under Assumption [A1], by continuity, there are positive constants C 0 , 0 such that for z ∈ B (µ, 0 ),


1
≤ C 0.
E
kX − yk2
Moreover, note that for all y ∈ K,
Z

1

Γθ+t(y−θ) (y − θ)dt.

Φ(y) =
0

Thus, with analogous calculus to the ones in the proof of Lemma 5.1 in [Cardot et al., 2015],
one can check that there is a positive constant C 00 such that for all y ∈ B (θ, 0 ) ∩ K,
kΦ(y) − Γθ k ≤ C 00 ky − θk2 .
Moreover, for all y = (z, a) ∈ K and y 0 = (z 0 , a0 ) ∈ Rd × R,

h

i
h
i
hX−z,z 0 i(X−z)
1
X−z
0
0
z
−
yE
z
−
+
a
E
2
kX−zk
kX−zk
kX−zk
h
i
.
Γy (y 0 ) = 
hX−z,z 0 i
0
E kX−zk + a

20

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

Thus, applying Cauchy-Schwarz’s inequality,





hX − z, z 0 i (X − z)
X −z
1
0
0
0 2
0
z −
+aE
kΓy (y )k = z − aE
kX − zk
kX − zk
kX − zk2


2
hX − z, z 0 i
+ E
+ a0
kX − zk

2
1
2
2
2
2
0 2
0 2
≤ 3 kz k + 3 kak kz k E
+ 3 ka0 k + 2 kz 0 k + 2 ka0 k
kX − zk

2

Thus, applying Lemma A.1, there are positive constants A1 , A2 such that
kΓy (y 0 )k ≤ A1 ky 0 k + A2 kyk ky 0 k
Note that since K is compact and convex, there is a positive constant CK such that for
all y ∈ K and t ∈ [0, 1], kθ + t(y − θ)k ≤ CK , and in a particular case,
Γθ+(y−θ) (y − θ) ≤ (A1 + A2 CK ) ky − θk .
Thus, for all y ∈ K such that ky − θk ≥ 0 ,
Z
kΦ(y) − Γθ (y − θ)k ≤

1

Γθ+t(y−θ) (y − θ) dt

0

≤ (A1 + A2 CK ) ky − θk
1
≤ 0 (A1 + A2 CK ) ky − θk2 .


Thus, we conclude the proof taking Cθ = max C 00 , 10 (A1 + A2 CK ) .



Appendix B. Proof of Section 4
Proof of Theorem 4.1. Let us recall that there is a positive constant c such that for all
y ∈ K, hΦ(y), y − θi ≥ c ky − θk2 . The aim is to use previous inequality
and the fact


2

that the projection is 1-lipschitz in order to get an upper bound of E

θbn+1 − θ

|Fn

and apply Robbins-Siegmund theorem to get the almost sure convergence of the algorithm.
Almost sure convergence of the algorithm: Since π is 1-lipschitz,



2
2
θbn+1 − θ = π θbn − γn ∇y g Xn+1 , θbn − π (θ)


2
≤ θbn − γn ∇y g Xn+1 , θbn − θ
D


E


2
= θbn − θ − 2γn ∇y g Xn+1 , θbn , θbn − θ + γn2 ∇y g Xn+1 , θbn

2

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

21

Thus, since θbn is Fn -measurable,



D h




i
E
2
2
2
b
b
b
b
b
E θn+1 − θ |Fn ≤ θn − θ − 2γn E ∇y g Xn+1 , θn Fn , θn − θ + γn E ∇y g Xn+1 , θn



 2
D
E
2
2
b
b
b
b
= θn − θ − 2γn Φ(θn ), θn − θ + γn E ∇y g Xn+1 , θn
|Fn



 2
2
2
2
≤ θbn − θ − 2cγn θbn − θ + γn E ∇y g Xn+1 , θbn
|Fn
Moreover, let θbn := (Zn , An ) with Zn ∈ Rd and An ∈ R, we have
"



 2
Zn − Xn+1
E ∇y g Xn+1 , θbn
|Fn = E Zn − Xn+1 − An
kZn − Xn+1 k


≤ 4E kZn − Xn+1 k2 Fn + 4 (An )2



Fn + E |An − kZn − Xn+1 k|2 |Fn

2

2

− 2cγn θbn − θ

+ γn2 M.
2

Applying Robbins-Siegmund’s theorem (see [Duflo, 1997] for instance), θbn − θ
verges almost surely to a finite random variable, and in a particular case,
∞
X
2
γk θbk − θ < +∞.

con-

k=1

Thus, since

P

k≥1

γk = +∞,
2

(B.2)

lim

n→+∞

θbn − θ

= 0 a.s.

Number ofPtimes the projection is used
Let Nn := nk=1 1{θbk −γk ∇y (Xk+1 ,θbk )∈K
/ } . This sequence is non-decreasing, and suppose by contradiction that Nn goes to infinity. Thus, there
(nk ) such
 is a subsequence

that (Nnk ) is increasing, i.e for all k ≥ 1, θbnk − γn ∇y g Xnk +1 , θbnk ∈
/ K, and in a
particular case, θbnk +1 ∈ ∂K, where ∂K is the frontier of K. Let us recall that θ is in
the interior of K, i.e let dmin := inf y∈∂K kθ − yk, we have dmin > 0. Thus,
θbnk +1 − θ ≥ dmin

a.s,

and,
lim θbnk +1 − θ = 0 ≥ dmin > 0 a.s,

k→∞

which leads to a contradiction.

Fn

#

2



≤ 8kZn − µk2 + 8(r? )2 + 8 (An − r? )2 + 8E kµ − Xn+1 k2 |Fn
2
 
= 8 θbn − θ + 8(r? )2 + 8r2 E W 2 .
Let M := 8(r? )2 + 8r2 E [W 2 ], we have


2

b
(B.1)
E θn − θ |Fn ≤ 1 + 8γn2 θbn − θ



2



22

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

Proof of Theorem 4.2. Convergence in quadratic mean
The aim is to obtain an induction relation for the quadratic mean error. Let us recall
inequality (B.1),


2
2
2

b
E θn+1 − θ |Fn ≤ 1 + 8γn2 θbn − θ − 2cγn θbn − θ + γn2 M.
Then we have

E

2


≤ 1 − cγn +

θbn+1 − θ

8γn2




E

2



θbn − θ

+ M γn2 ,

and one can conclude the proof with the help of an induction (see [Godichon, 2015] for
instance) or applying a lemma of stabilization (see [Duflo, 1996]).
Lp rates of convergence
Let p ≥ 2, we now prove with the help of a strong induction that for all integer p0 ≤ p,
there is a positive constant Cp0 such that for all n ≥ 1,


2p0
Cp 0
b
≤ p0 α .
E θn − θ
n
This inequality is already checked for p0 = 1. Let p0 ≥ 2, we suppose from now that for
all integer k < p0 , there is a positive constant Ck such that for all n ≥ 1,


2k
Ck
≤ kα .
E θbn − θ
n
0

We now search to give an induction relation for the L2p -error. Let us recall that
D


E

 2
2
2
.
θbn+1 − θ ≤ θbn − θ − 2γn ∇y g Xn+1 , θbn , θbn − θ + γn2 ∇y g Xn+1 , θbn
 
We suppose from now that E [W 2p ] < +∞(and in aparticular case, E W k < +∞
for all integer k ≤ 2p) and let Un+1 := ∇y g Xn+1 , θbn . We have
2p0

θbn+1 − θ



2

θbn − θ

≤

+

γn2

kUn+1 k

2

p0

0

D

− 2p γn θbn − θ, Un+1

E

2

θbn − θ

+

γn2

(B.3)
0

+

p  0
X
p
k=2

k

γnk

D

θbn − θ, Un+1

E k

2

θbn − θ

+

γn2

kUn+1 k

2

p0 −k
.

The aim is to bound each term on the right-hand side of previous inequality. In this
purpose, we first need to introduce some technical inequalities.



 2

 2
2
b
b
kUn+1 k ≤ 2 ∇y g Xn+1 , θn
+ 2E ∇y g Xn+1 , θn
Fn


2
 2
2
2
?
2
.
≤ 16 2 θbn − θ + 2 (r ) + kµ − Xn+1 k + r E W

2

kUn+1 k

p0 −1

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

23

Thus, applying Lemma A.1 in [Godichon, 2015] for instance, for all integer k ≤ p0 ,


2k
 2 k
2k
2k
k
? 2k
2k
k b
k−1
k
.
+ 2 (r ) + kXn+1 − µk + r E W
kUn+1 k ≤ 4 16 2 θn − θ


In a particular case, since for all k ≤ p, E W 2k < +∞, there are positive constants
A1,k , A2,k such that for all n ≥ 1,


h
i
2k
 2 k
 2k 
2k
2k
k
? 2k
2k
k b
3k−1
+r E W
+ 2 (r ) + r E W
2 θn − θ
E kUn+1 k Fn ≤ 4
2k

≤ A1,k θbn − θ

(B.4)

+ A2,k .

We can now bound the expectation of the three terms on the right-hand side of inequality (B.3). First, since θbn is Fn - measurable, applying inequality (B.4), let

p0 
2
2
(∗) := E
θbn − θ + γn kUn+1 k
2p0



0



θbn − θ

=E

+

p  0
X
p
k=1

2p0


≤E

0



θbn − θ

+

p  0
X
p
k=1

Let B :=

Pp0

2k−2
A1,k ,
k=1 cγ

k

γn2k E


kUn+1 k2k θbn − θ

A1,k θbn − θ

2p0 −2k



2p0 −2k



2k

+ A2,k

θbn − θ

using previous inequality and by induction,
0


p  0
X
p
2k
(∗) ≤ 1 +
+
E
γn A2,k E θbn − θ
k
k=1
 X

p0  0 
0
2p

p 2k
Ck
≤ 1 + Bγn2 E θbn − θ
+
cγ A2,k (p0 +k)α
k
n
k=1


 0 
2p0

+ O γnp +1 .
≤ 1 + Bγn2 E θbn − θ
Bγn2

(B.5)

k

γn2k E



2p0





2p0 −2k



θbn − θ

In the same way, applying Cauchy-Schwarz’s inequality, let
"
p0 −1 #
D
E
2
2
0
2
(∗∗) := −2p γn E θbn − θ, Un+1
θbn − θ + γ kUn+1 k
n

0

≤ −2p γn E

D

θbn − θ, Un+1

E

"
+ 2p0 γn E

θbn − θ kUn+1 k

2p0 −2



θbn − θ

p0 −1  0
X
k=1


p − 1 2k
γn kUn+1 k2k θbn − θ
k

2p0 −2k

#
.



24

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

Moreover, since θbn is Fn -measurable, applying Proposition 3.1,


D
D
E
E
2p0 −2
2p0 −2
0
0
b
b
b
b
= −2p γn E θn − θ, E [Un+1 |Fn ] θn − θ
−2p γn E θn − θ, Un+1 θn − θ

D
E
2p0 −2
0
b
b
b
= −2p γn E θn − θ, Φ(θn ) θn − θ


2p0
0
.
≤ −2p cγn E θbn − θ
Moreover, since 2ab ≤ a2 + b2 , let
#
"
p0 −1  0
X p − 1
2p0 −2k
γn2k kUn+1 k2k θbn − θ
(∗∗0 ) := 2p0 γn E θbn − θ kUn+1 k
k
k=1
#
"
0 −1 
 pX

0
2p0 −2k
2
p
−
1
2k
2
θbn − θ
θbn − θ + kUn+1 k
γn2k kUn+1 k
≤ p0 γn E
k
k=1
0






p −1
X p0 − 1
2p0 −2k
2k+2 b
0
2k
≤ p γn
γn E kUn+1 k
θn − θ
+ E kUn+1 k2k θbn − θ
k
k=1

2p0 +2−2k


.

With analogous calculus to the ones for inequality (B.5), one can check that there is a
positive constant B 0 such that for all n ≥ 1,


 0

2p0
0
0 2
b
+ O γn(p +1)α .
(∗∗ ) ≤ B γn E θn − θ
Thus,
(B.6)
"
−2γn E γn

D

θbn − θ, Un+1

E

2

θbn − θ

+

γn2

2

kUn+1 k

p0 −1 #

0

≤ −2cp γn +

B 0 γn2



2p0


E



θbn − θ

 0

+ O γn(p +1)α .
Finally, applying Lemma A.1 in [Godichon, 2015] and since |ha, bi| ≤ 21 kak2 + 12 kbk2 ,
let
"
p0 −k #
p0  0 
D
E k
X
2
p
(∗ ∗ ∗) :=
γnk E
θbn − θ, Un+1
θbn − θ + γn2 kUn+1 k2
k
k=2
"
k 
p0 −k #
p0  0 
X
2
2
p
1
1
2
2
θbn − θ + kUn+1 k
θbn − θ + γn2 kUn+1 k
≤
γnk E
2
2
k
k=2



p0  0 
X
2k
2p0 −2k
0 −2k
p p0 −k−2 k
0 −2k
2k
2p
2p
kUn+1 k
≤
2
γn E
θbn − θ
+ kUn+1 k
θbn − θ
+ γn
k
k=2

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

25

Thus, with analogous calculus to the ones for inequality (B.5), one can check that there
is a positive constant B 00 such that for all n ≥ 1,
(B.7)
"
p0 −k #


p0  0 
D
E
X
2
k
p
2
00
2
k
2
≤ B γn E θbn − θ
γn E
θbn − θ, Un+1
θbn − θ + γn kUn+1 k
k
k=2
Finally, applying inequalities (B.5) to (B.7), there are positive constants B1 , B2 such
that for all n ≥ 1,




2p0
2p0

0
0
2
+ B2 γnp +1 .
≤ 1 − 2p cγn + B1 γn E θbn − θ
(B.8)
E θbn+1 − θ
Thus, with the help of an induction on n or applying a lemma of stabilization (see
[Duflo, 1996] for instance), one can check that there is a positive constant Cp0 such that
for all n ≥ 1,


2p0
Cp 0
≤ p0 α ,
E θbn − θ
n
which concludes the induction on p0 and the proof.
h


i
Bounding P θbn − γn ∇y g Xn+1 , θbn ∈
/K
Let us recall that dmin = inf y∈∂K ky − θk > 0 and that if W admits a2p-th moment,
2p
≤ nCpαp . Thus,
there is a positive constant Cp such that for all n ≥ 1, E θbn − θ
for all n ≥ 1,
Cp
≥E
(n + 1)pα



2p



θbn+1 − θ


2p



1{θbn −γn ∇y g(Xn+1 ,θbn )∈K
≥ E θbn+1 − θ
/ }
h


i
b
b /K .
≥ d2p
min P θn − γn ∇y g Xn+1 , θn ∈
Finally,
h


i
Cp
1
Cp 1
P θbn − γn ∇y g Xn+1 , θbn ∈
/ K ≤ 2p
≤ 2p pα .
pα
dmin (n + 1)
dmin n

Proof of Theorem 4.3. The aim is, in a first time, to exhibit a nice decomposition of
the averaged algorithm. In this purpose, let us introduce this new decomposition of the
PRM algorithm
 
(B.9)
θbn+1 − θ = θbn − θ − γn Φ θbn + γn ξn+1 + rn ,

2p0


+O



0
γnp +1



.

26

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

with


 
ξn+1 := −∇y g Xn+1 , θbn + Φ θbn ,





b
b
b
b
rn := π θn − γn ∇y g Xn+1 , θn − θn + γn ∇y g Xn+1 , θn .
Remark that (ξn ) is a sequence of martingale
adapted to the filtration (Fn )
 differences

b
b
and rn is equal to 0 when θn − γn ∇y g Xn+1 , θn ∈ K. Moreover, linearizing the
gradient, decomposition (B.9) can be written as


θbn+1 − θ = (IRd ×R − γn Γθ ) θbn − θ + γn ξn+1 − γn δn + rn ,

(B.10)

 


b
b
where δn := Φ θn − Γθ θn − θ is the remainder term in the Taylor’s expansion of
the gradient. This can also be decomposed as

 θb − θ θb − θ
rn
n+1
n
Γθ θbn − θ =
−
− δn +
+ ξn+1 .
γn
γn
γn
As in [Pelletier, 2000], summing these equalities, applying Abel’s transform and dividing
by n,
(B.11)
Γθ

1
θn − θ =
n



n
n
n 
 X
X
rk
θb1 − θ θbn+1 − θ X 1
1
b
θk − θ −
δk +
−
+
−
γ1
γn
γk γk−1
γk
k=1
k=1
k=2

!

We now give the rate of convergence in quadratic mean of each term using Theorem 4.2.
In this purpose, let us recall the following technical lemma.
Lemma B.1 ([Godichon, 2015]). Let Y1 , ..., Yn be random variables taking values in a
normed vector space such that for all positive constant q and for all k ≥ 1, E [kYk kq ] < ∞.
Thus, for all constants a1 , ..., an and for all integer p,
"
(B.12)

E

n
X
k=1

p#

ak Yk

≤

n
X

|ak | (E [kYk kp ])

k=1

The remainder terms: First, one can check that


2
 
b
1  θ1 − θ 
1
(B.13)
E
=o
.
2
n
γ1
n

1
p

!p

n

1X
+
ξk+1 .
n k=1

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

27

In the same way, applying Theorem 4.2,


2


b
2
1 1
1  θn+1 − θ 
b
= 2 2−2α E θn+1 − θ
E
n2
γn
cγ n
C1 1
c2γ n2−α
 
1
.
=o
n

≤
(B.14)

−1
α−1
Moreover, since γk−1 − γk−1
≤ 2αc−1
, applying Lemma B.1,
γ k

! 2
s 
n 
X
2
1
1
E θbk − θ
−
γ
γ
k
k−1
k=2
!2
n
X
C
4α2 c−2
1
1
γ
≤
n2
k 1−α/2
k=2


1
=O
n2−α
 
1
.
=o
n


2

n 

X
1 
1
1
≤ 1
−
θbk − θ
E
2
n
γk γk−1
n2
k=2


(B.15)

Thanks to Lemma A.1, there is a positive constant Cθ such that for all n ≥ 1,
2

kδn k ≤ Cθ θbn − θ

.

Thus, applying Lemma B.1 and Theorem 4.2, there is a positive constant C2 such that


!2
2
n
n q
X
X


1
1 
E kδk k2
E
δk  ≤ 2
n2
n
k=1
k=1
s 
! 2
n
4
Cθ2 X
≤ 2
E θbk − θ
n
k=1
!2
n
Cθ2 C2 X 1
≤
n2
kα
k=1


1
=O
n2α
 
1
=o
.
n

28

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER



Let Un+1 := ∇y g Xn+1 , θbn , note that if θbn − γn Un+1 ∈ K, then rn = 0. Thus,
applying Lemma B.1 and Cauchy-Schwarz’s inequality,


!2
2
n
n
q 
X

1 X 1
1 
rk 
≤ 2
E krk k2
E
n2
γ
n
γ
k
k
k=1
k=1
!
r h
n
i 2
1 X 1
2
= 2
E krk k 1θbk −γk Uk+1 ∈K
/
n
γk
k=1
!
n
h
i 14 2
1 


1 X 1
≤ 2
E krk k4 4 P θbk − γk Uk+1 ∈
.
/K
n
γ
k
k=1
Moreover, since π is 1-lipschitz,


4
4
b
krn k = π θn − γn Un+1 − θ + θ − θbn + γn Un+1
4
 

b
b
≤ π θn − γn Un+1 − π (θ) + θn − γn Un+1 − θ

4
b
≤ 2 θn − θ − γn Un+1
≤ 27 θbn − θ

4

+ 27 γn4 kUn+1 k2 .

Thus, applying inequality (B.4), there are positive constants A1 , A2 such that for all
n ≥ 1,
i
h
4
E krn k4 Fn ≤ A1 θbn − θ + A2 γn4 .
In a particular case, applying Theorem 4.2, there is a positive constant A3 such that for
all n ≥ 1,


A3
E krn k4 ≤ 2α .
n
Moreover, applying Theorem 4.2, there is a positive constant C6 such that for all n ≥ 1,
h
i
C6
P θbn − γn Un+1 ∈
/ K ≤ 12 6α .
dmin n
Then,


n
1  X rk
E
n2
γk
k=1

(B.16)

2



√

n
X
1
 ≤ C6 A 3
6
2
dmin cγ n
kα
k=1


1
=O
n2α
 
1
=o
.
n

!2

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

29

The martingale term: Since (ξn ) is a sequence of martingale differences adapted
to the filtration (Fn ),


2
n
n
n
n
X
1 
1 X 
2 X X
2

= 2
E [hξk+1 , ξk0 +1 i]
E
ξk+1
E kξk+1 k + 2
n2
n k=1
n k=1 k0 =k+1
k=1
=

n
n
n
1 X 
2 X X
2
E [hξk+1 , E [ξk0 +1 |Fk0 ]i]
E
kξ
k
+
k+1
n2 k=1
n2 k=1 k0 =k+1

=

n
1 X 
2
E
kξ
k
.
k+1
n2 k=1

Moreover,

hD
Ei
h
 2

2
2
b
b
E kξn+1 k = E kUn+1 k − 2E E [Un+1 |Fn ] , Φ(θn + E Φ(θn


2

2
= E kUn+1 k − E Φ(θbn )


≤ E kUn+1 k2 .
Finally, applying inequality (B.4) and Theorem 4.2, there is a positive constant M such
that


2

2
E kξn+1 k ≤ A1,1 E θbn − θ
+ A2,1
≤ M.
Then,

1 
E
n2

n
X

2



n
1 X
ξk+1  ≤ 2
M
n
k=1
k=1

=
which concludes the proof.

M
,
n


Proof of Theorem 4.4. Let us recall that the averaged algorithm can be written as follows
(B.17)
!

n
n
n
n 
 X
X

√
1
θb1 − θ θbn+1 − θ X 1
1
r
1 X
k
nΓθ θn − θ = √
−
+
−
θbk − θ −
δk +
+√
ξk+1 .
γ1
γn
γ
γ
γ
n
n
k
k−1
k
k=2
k=1
k=1
k=1
We now prove that the first terms on the right-hand side of previous equality converge
in probability to 0 and apply a Central Limit Theorem to the last one.

30

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

The remainder terms: Applying inequalities (B.13) to (B.16),
1 θb − θ P
√ 1
−−−→ 0,
n γ1 n→∞
1 θb − θ P
√ n+1
−−−→ 0,
n→∞
n γn


n


1 X 1
1
P
√
−
θbk − θ −−−→ 0,
n→∞
n k=2 γk γk−1
n

1 X rk P
√
−−−→ 0.
n k=1 γk n→∞

The martingale term Let θbn = (Zn , An ) ∈ Rd × R, then ξn+1 can be written as
0
ξn+1 = ξn+1
+ n+1 + 0n+1 , with
µ − Xn+1 − E [µ − Xn+1 |Fn ] − r?



µ−Xn+1
kµ−Xn+1 k

−E

h

µ−Xn+1
kµ−Xn+1 k

Fn

i!

0
ξn+1

:=

n+1

i
!
Fn
:= −

 ,
kZn − Xn+1 k − E [kZn − Xn+1 k Fn ] − kµ − Xn+1 k + E kµ − Xn+1 k Fn

,



r? − kµ − Xn+1 k − r? + E kµ − Xn+1 k Fn
(An − r? )

0n+1 := −

r

?



Zn −Xn+1
kXn+1 −Zn k

−



µ−Xn+1
kµ−Xn+1 k

Zn −Xn+1
kZn −Xn+1 k

−E
0

h

−E

Zn −Xn+1
kXn+1 −Zn k

h

Zn −Xn+1
kZn −Xn+1 k

−

µ−Xn+1
kµ−Xn+1 k

Fn

i!
.

Note that (ξn ) , (n ) , (0n ) are martingale differences sequences adapted to the filtration (Fn ). Thus,

1 
E
n

n
X
k=1

2

k+1



n
n
n
X

2X X
2
= 1
E kk+1 k +
E [hk+1 , k0 +1 i]
n k=1
n k=1 k0 =k+1

=

n
n
n
X
 2X
1X 
E kk+1 k2 +
E [hk+1 , E [k0 +1 |Fk0 ]i]
n k=1
n k=1 k0 =k+1

=


1X 
E kk+1 k2 .
n k=1

n

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

31

Moreover,
E kn+1 k2 |Fn = E




" 

−Xn+1
(An − r? ) kZZnn −X
n+1 k
kZn − Xn+1 k − kµ − Xn+1 k



2

#
Fn

!
h
i
n+1
(An − r? ) E kZZnn −X
F
n
−Xn+1 k
E [kZn − Xn+1 k − kµ − Xn+1 k |Fn ]

+

2

"*

!+
#
h
i

−Xn+1
Zn −Xn+1
?
(An − r? ) kZZnn −X
(A
−
r
)
E
F
n
n
kZn −Xn+1 k
n+1 k
,
Fn .
kZn − Xn+1 k − kµ − Xn+1 k
E [kZn − Xn+1 k − kµ − Xn+1 k |Fn ]

− 2E
Thus, one can check that



E kn+1 k2 |Fn ≤ E

" 

?

Zn −Xn+1
kZn −Xn+1 k

− (An − r )
kZn − Xn+1 k − kµ − Xn+1 k



2

#
Fn

≤ kAn − r? k2 + kZn − µk2
2

= θbn − θ

.

Thus, applying Theorem 4.2,


2
n
n
X

1 
1X 

k+1
E
=
E kk+1 k2
n
n k=1
k=1


n
2
1X
≤
E θbn − θ
n k=1
n

1 X C0
≤
n k=1 nα
 
1
=O
.
nα
As a particular case,
n

1 X
P
√
k+1 −−−→ 0.
n→∞
n k=1
Similarly, since

Zn −Xn+1
kZn −Xn+1 k

−

µ−Xn+1
kµ−Xn+1 k


h
E

0n+1

2

r?

i

|Fn ≤ E 



≤ 2,
Zn −Xn+1
kXn+1 −Zn k

0
? 2

≤ 2 (r ) E



−

µ−Xn+1
kµ−Xn+1 k

!



2

Fn 

Zn − Xn+1
µ − Xn+1
−
kXn+1 − Zn k kµ − Xn+1 k


Fn .

32

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

This last term is closely related to the gradient of the function we need to minimize
to get the geometric median (see [Kemperman, 1987] for example) and it is proved in
[Cardot et al., 2015] that since Lemma A.1 is verified, then


µ − Xn+1
Zn − Xn+1
−
E
Fn ≤ C kZn − µk ≤ C θbn − θ .
kXn+1 − Zn k kµ − Xn+1 k
Thus,


E kn+1 k2 |Fn ≤ 2C (r? )2 θbn − θ .
Finally, since (n ) is a sequence of martingale differences adapted to the filtration (Fn ),
applying Theorem 4.2 and Cauchy-Schwarz’s inequality,


2
n
n
i
X
1 
1X h 0
2
0

k+1
=
E k+1
E
n
n k=1
k=1
≤ 2C (r? )2

n
i
1X h b
E θn − θ
n k=1

n
p 1X
1
≤ 2C (r )
C1
α/2
n k=1 k


1
=O
.
nα/2
Note that with more assumptions on W , we could get a better rate but this one is
sufficient. Indeed, thanks to previous inequality,
n
1 X 0
P
√

−−−→ 0.
n k=1 k+1 n→∞
? 2

Finally, applying a Central Limit Theorem (see [Duflo, 1997] for example), we have the
convergence in law
n
1 X 0
L
√
(B.18)
ξ
−−−→ N (0, Σ) ,
n k=1 k+1 n→∞
with



µ−X 
µ−X 
µ − X − r? kµ−Xk
µ − X − r? kµ−Xk
Σ := E
⊗
,
r? − kµ − Xk
r? − kµ − Xk
which also can be written as
 ?
  ?

r UΩ − rW UΩ
r UΩ − rW UΩ
Σ=E
⊗
.
r? − rW
r? − rW
Thus, we have the convergence in law
 L
√
nΓθ θn − θ −−−→ N (0, Σ) ,
n→∞

and in a particular case,
√


 L
−1
n θn − θ −−−→ N 0, Γ−1
.
θ ΣΓθ
n→∞

ESTIMATING THE PARAMETERS OF A TRUNCATED SPHERICAL DISTRIBUTION

33



References
[Abuzaina et al., 2013] Abuzaina, A., Nixon, M. S., and Carter, J. N. (2013). Sphere detection in
kinect point clouds via the 3d hough transform. In Computer Analysis of Images and Patterns, pages
290–297. Springer.
[Al-Sharadqah, 2014] Al-Sharadqah, A. (2014). Further statistical analysis of circle fitting. Electronic
Journal of Statistics, 8(2):2741–2778.
[Beck and Pan, 2012] Beck, A. and Pan, D. (2012). On the solution of the gps localization and circle
fitting problems. SIAM Journal on Optimization, 22(1):108–134.
[Bercu and Fraysse, 2012] Bercu, B., Fraysse, P., et al. (2012). A robbins–monro procedure for estimation in semiparametric regression models. The Annals of Statistics, 40(2):666–693.
[Brazey and Portier, 2014] Brazey, D. and Portier, B. (2014). A new spherical mixture model for head
detection in depth images. SIAM Journal on Imaging Sciences, 7(4):2423–2447.
[Breiman and Friedman, 1985] Breiman, L. and Friedman, J. H. (1985). Estimating optimal transformations for multiple regression and correlation. Journal of the American statistical Association,
80(391):580–598.
[Cardot et al., 2015] Cardot, H., Cénac, P., and Godichon, A. (2015). Online estimation of the geometric median in hilbert spaces: non asymptotic confidence balls. arXiv preprint arXiv:1501.06930.
[Duflo, 1996] Duflo, M. (1996). Algorithmes stochastiques.
[Duflo, 1997] Duflo, M. (1997). Random iterative models, volume 22. Springer Berlin.
[Fischler and Bolles, 1981] Fischler, M. A. and Bolles, R. C. (1981). Random sample consensus: a
paradigm for model fitting with applications to image analysis and automated cartography. Communications of the ACM, 24(6):381–395.
[Gahbiche and Pelletier, 2000] Gahbiche, M. and Pelletier, M. (2000). On the estimation of the asymptotic covariance matrix for the averaged robbins–monro algorithm. Comptes Rendus de l’Académie des
Sciences-Series I-Mathematics, 331(3):255–260.
[Godichon, 2015] Godichon-Baggioni, A. (2015) Estimating the geometric median in Hilbert spaces
with stochastic gradient algorithms: Lp and almost sure rates of convergence. Journal of Multivariate
Analysis. Elsevier.
[Jiang and Jiang, 1998] Jiang, B. C. and Jiang, S. (1998). Machine vision based inspection of oil seals.
Journal of manufacturing systems, 17(3):159–166.
[Kemperman, 1987] Kemperman, J. (1987). The median of a finite measure on a banach space. Statistical data analysis based on the L1-norm and related methods (Neuchâtel, 1987), pages 217–230.
[Kushner and Yin, 2003] Kushner, H. J. and Yin, G. (2003). Stochastic approximation and recursive
algorithms and applications, volume 35. Springer Science & Business Media.
[Lacoste-Julien et al., 2012] Lacoste-Julien, S., Schmidt, M., and Bach, F. (2012). A simpler approach
to obtaining an o (1/t) convergence rate for the projected stochastic subgradient method. arXiv
preprint arXiv:1212.2002.
[Landau, 1987] Landau, U. (1987). Estimation of a circular arc center and its radius. Computer Vision,
Graphics, and Image Processing, 38(3):317–326.
[Liu and Wu, 2014] Liu, J. and Wu, Z.-k. (2014). An adaptive approach for primitive shape extraction
from point clouds. Optik-International Journal for Light and Electron Optics, 125(9):2000–2008.
[Martins et al., 2008] Martins, D. A., Neves, A. J., and Pinho, A. J. (2008). Real-time generic ball
recognition in robocup domain. In Proc. of the 3rd International Workshop on Intelligent Robotics,
IROBOT, pages 37–48.
[Muirhead, 2009] Muirhead, R. J. (2009). Aspects of multivariate statistical theory, volume 197. John
Wiley & Sons.
[Pelletier, 1998] Pelletier, M. (1998). On the almost sure asymptotic behaviour of stochastic algorithms.
Stochastic processes and their applications, 78(2):217–244.

34

ANTOINE GODICHON-BAGGIONI AND BRUNO PORTIER

[Pelletier, 2000] Pelletier, M. (2000). Asymptotic almost sure efficiency of averaged stochastic algorithms. SIAM Journal on Control and Optimization, 39(1):49–72.
[Polyak and Juditsky, 1992] Polyak, B. T. and Juditsky, A. B. (1992). Acceleration of stochastic approximation by averaging. SIAM Journal on Control and Optimization, 30(4):838–855.
[Rabbani, 2006] Rabbani, T. (2006). Automatic reconstruction of industrial installations using point
clouds and images. Technical report, NCG Nederlandse Commissie voor Geodesie Netherlands Geodetic
Commission.
[Robbins and Monro, 1951] Robbins, H. and Monro, S. (1951). A stochastic approximation method.
The annals of mathematical statistics, pages 400–407.
[Rusu et al., 2003] Rusu, C., Tico, M., Kuosmanen, P., and Delp, E. J. (2003). Classical geometrical
approach to circle fitting?review and new developments. Journal of Electronic Imaging, 12(1):179–193.
[Schnabel et al., 2007] Schnabel, R., Wahl, R., and Klein, R. (2007). Efficient ransac for point-cloud
shape detection. In Computer graphics forum, volume 26, pages 214–226. Wiley Online Library.
[Shafiq et al., 2001] Shafiq, M. S., Tümer, S. T., and Güler, H. C. (2001). Marker detection and trajectory generation algorithms for a multicamera based gait analysis system. Mechatronics, 11(4):409–437.
[Thom, 1955] Thom, A. (1955). A statistical examination of the megalithic sites in britain. Journal of
the Royal Statistical Society. Series A (General), pages 275–295.
[Tran et al., 2015] Tran, T.-T., Cao, V.-T., and Laurendeau, D. (2015). Extraction of cylinders and
estimation of their parameters from point clouds. Computers & Graphics, 46:345–357.
[Von Hundelshausen et al., 2005] Von Hundelshausen, F., Schreiber, M., and Rojas, R. (2005). A constructive feature detection approach for robotic vision. In RoboCup 2004: Robot Soccer World Cup
VIII, pages 72–83. Springer.
[Zhang et al., 2007] Zhang, X., Stockel, J., Wolf, M., Cathier, P., McLennan, G., Hoffman, E. A., and
Sonka, M. (2007). A new method for spherical object detection and its application to computer aided
detection of pulmonary nodules in ct images. In Medical Image Computing and Computer-Assisted
Intervention–MICCAI 2007, pages 842–849. Springer.
E-mail address: Antoine.Godichon@u-bourgogne.fr
E-mail address: Bruno.Portier@insa-rouen.fr
Université de Bourgogne, Institut de Mathématiques de Bourgogne, 9 Avenue
Alain Savary, 21000 Dijon, France
Normandie Université, INSA de Rouen, Laboratoire de Mathématiques, BP 08 Avenue de l’Université, 76131 Saint-Etienne-du-Rouvray, France

