SIAM J. NUMER. ANAL. c  2006 Society for Industrial and Applied Mathematics Vol. 44, No. 1, pp. 102–119 CONVERGENCE OF THE LLOYD ALGORITHM FOR COMPUTING CENTROIDAL VORONOI TESSELLATIONS ∗ QIANG DU † , MARIA EMELIANENKO ‡ , AND LILI JU § Abstract. Centroidal Voronoi tessellations (CVTs) are Voronoi tessellations of a bounded geometric domain such that the generating points of the tessellations are also the centroids (mass centers) of the corresponding Voronoi regions with respect to a given density function. Centroidal Voronoi tessellations may also be defined in more abstract and more general settings. Due to the natural optimization properties enjoyed by CVTs, they have many applications in diverse fields. The Lloyd algorithm is one of the most popular iterative schemes for computing the CVTs but its theoretical analysis is far from complete. In this paper, some new analytical results on the local and global convergence of the Lloyd algorithm are presented. These results are derived through careful utilization of the optimization properties shared by CVTs. Numerical experiments are also provided to substantiate the theoretical analysis. Key words. centroidal Voronoi tessellations, k-means, optimal vector quantizer, Lloyd algo- rithm, global convergence, convergence rate AMS subject classifications. 65D99, 65C20 DOI. 10.1137/040617364 1. Introduction. A centroidal Voronoi tessellation (CVT) is a special Voronoi tessellation of a given set such that the associated generating points are the centroids (centers of mass) of the corresponding Voronoi regions with respect to a predefined density function [7]. CVTs are indeed special as they enjoy very natural optimization properties which make them very popular in diverse scientific and engineering appli- cations that include art design, astronomy, clustering, geometric modeling, image and data analysis, resource optimization, quadrature design, sensor networks, and numer- ical solution of partial differential equations [1, 2, 3, 4, 7, 8, 9, 10, 11, 13, 14, 17, 15, 26, 29, 30, 31, 39, 44, 45]. In particular, CVTs have been widely used in the design of optimal vector quantizers in electrical engineering [25, 28, 40, 43]. They are also related to the so-called method of k-means [27] in clustering analysis. CVTs can also be defined in more general cases such as those constrained to a manifold [12, 11] or those corresponding to anisotropic metrics [16, 18], and other abstract settings [7, 9]. For modern applications of the CVT concept in large-scale scientific and engineer- ing problems, it is important to develop robust and efficient algorithms for construct- ing CVTs in various settings. Historically, a number of algorithms have been studied and widely used [7, 19, 25, 27, 38]. A seminal work is the algorithm first developed in the 1960s at Bell Laboratories by S. Lloyd which remains to this day one of the most popular methods due to its effectiveness and simplicity. The algorithm was later officially published in [35]. It is now commonly referred to as the Lloyd algorithm and is the main focus of this paper. ∗ Received by the editors October 20, 2004; accepted for publication (in revised form) August 15, 2005; published electronically February 8, 2006. This work was supported in part by NSF grants DMS-0409297, CCF-0430349, and ITR DMR-0205232. http://www.siam.org/journals/sinum/44-1/61736.html † Department of Mathematics, Pennsylvania State University, University Park, PA 16802 (qdu@ math.psu.edu). ‡ Department of Mathematical Sciences, Carnegie Mellon University, PA 15213 (masha@cmu.edu). § Department of Mathematics, University of South Carolina, Columbia, SC 29208 (ju@math.sc. edu). 102