On the numerical approximation of Blaschke-Santalo diagrams using Centroidal Voronoi Tessellations

. Identifying Blaschke-Santal´o diagrams is an important topic that essentially consists in determining the image Y = F ( X ) of a map F : X → R d , where the dimension of the source space X is much larger than the one of the target space. In some cases, that occur for instance in shape optimization problems, X can even be a subset of an infinite-dimensional space. The usual Monte Carlo method, consisting in randomly choosing a number N of points x 1 , . . . , x N in X and plotting them in the target space R d , produces in many cases areas in Y of very high and very low concentration leading to a rather rough numerical identification of the image set. On the contrary, our goal is to choose the points x i in an appropriate way that produces a uniform distribution in the target space. In this way we may obtain a good representation of the image set Y by a relatively small number N of samples which is very useful when the dimension of the source space X is large (or even infinite) and the evaluation of F ( x i ) is costly. Our method consists in a suitable use of Centroidal Voronoi Tessellations which provides efficient numerical results. Simulations for two and three dimensional examples are shown in the paper.


Introduction
Given a map  = ( 1 , . . .,   ) :  → R  , finding the image  () amounts to finding all optimal inequalities between the quantities   ,  = 1, . . ., .In particular, the image  () simultaneously characterizes solutions to all optimization problems of the form min rigidity, the eigenvalues of the Laplace operator −∆, and other similar geometrical or analytical quantities (see for instance [5,7,19]).
Even when the set  is a subset of a finite dimensional Euclidean space R  , the identification of the image  () through a numerical procedure may present deep and unexpected difficulties.The naive approach consisting in generating a random uniform sampling of  by means of a discrete set {  } ∈ , associated to the corresponding outputs { (  )} ∈ , may not produce a satisfactory approximation.Some parts of the image  () may be very rarely explored by a uniform random sampling of .In these cases one needs a very large number of samples in order to have a quite accurate description of the set  ().This phenomenon has a dramatic impact when the evaluation of  on the chosen sample requires the solution of one or more partial differential equations.In shape optimization, for instance, this procedure may be too expensive in terms of computational time.Such contexts require us to develop a new approach to get a precise description of the image  () using a relatively small number of sampling points.The choice of the sampling has to be adjusted carefully in order to comply with the complexity of the map  .
In this paper we develop a new method based on Voronoi tessellations which seems much more efficient than the standard random uniform approach.We describe in the following sections the numerical method, and we show some algebraic and geometric examples in which the efficiency of our procedure is clearly outlined.
The study of the range of scale invariant ratios between geometric quantities was initiated by Santaló in [23] and Blaschke in [2].If the geometric image of scale invariant ratios is completely characterized, then all possible inequalities between these quantities are known.In practice, often three geometric functionnals are used to generate at least two scale invariant ratios.
This approach has been investigated recently in a shape optimization context.We mention [7] concerning the diagram given by the area, the diameter and the inradius.In [14] the inequalities between volume, perimeter and the first Dirichlet-Laplace eigenvalue are investigated.The Cheeger inequality was investigated in [11] and inequalities involving the first Dirichlet eigenvalue, the torsion and the volume are studied in [12].
In most situations, a complete analytical understanding of the resulting Blaschke-Santaló diagrams is not available.This motivates the use of numerical tools.A first method is generating random shapes and computing quantities of interest.This method is illustrated in some works cited above.A more rigorous numerical approach is solving numerical optimization problems finding extreme points for vertical or horizontal slices of the diagram.This method is used in [13] using methods described in [1] and [3] to perform numerical shape optimization among convex sets.
This article proposes a completely new alternative approach which generates uniformly distributed samples in the Blaschke-Santaló diagram.Implicitly, our method also provides boundary points for the diagram.Compared to [13] where multiple constrained numerical optimizations are solved, we use a global iterative process, and we solve numerically a global optimization problem providing a geometrical description of the diagram.We illustrate the method for an algebraic example involving the trace and determinant of symmetric matrices with entries in [−1, 1].This simple example is already non-trivial starting from 3 × 3 matrices.In a second stage we investigate numerically the diagram associated to the area, perimeter and moment of inertia among convex shapes with two axes of symmetry.
Structure of the paper.Section 2 presents the approximation framework based on Centroidal Voronoi Tessellations.Theoretical, algorithmic and practical aspects are discussed.In Section 3 the algebraic example involving the (trace, determinant) diagram among symmetric matrices with bounded entries is illustrated.This simple example is already challenging starting from dimension three.Section 4 contains further algorithmic details regarding multi-grid strategies, extraction of the Blaschke-Santaló diagram from the Voronoi diagram and comparisons with the Monte-Carlo method.An example of a three dimensional diagram is also presented.Section 5 deals with shape optimization aspects regarding the approximation of the diagram associated to the perimeter, area and moment of inertia among convex sets with two axes of symmetry.A new parametrization of convex shapes is given in this context, involving only pointwise bounds on the parameters.

Optimal transport framework
Consider a continuous map  :  → R  , with  a compact metric space.In order to have a careful description of the image set  () we could randomly choose some points { 1 , . . .,   } in  and, for a large , the set { (  ) :  = 1, . . ., } would give an approximate description of the full image set  ().However, even if  is a subset of an Euclidean finite dimensional space, due to the nonlinearity of the map  , the number  that is necessary to have a rather accurate description of the image set  () could be extremely high.In other words, a uniform random choice of points { 1 , . . .,   } in  does not produce in general a well distributed sequence  ( 1 ), . . .,  (  ), and concentration/rarefaction effects very often occur.We should then make the random choice of the points { 1 , . . .,   } in  according to a probability measure that is not uniform and that depends on the function  , in order to obtain a well distributed sequence  ( 1 ), . . .,  (  ).If ℒ is the Lebesgue measure on  () the probability measure  governing the random choice of points on  should then be such that  #  = ℒ, where  # is the push-forward operator related to the function  , defined by for every measurable set  ⊂  ().We refer to [24] for all details related to optimal transport theory.
We prove the following elementary result which we did not find in the literature.
Then, for every probability measure  on  there exists a probability measure  on  such that Proof.Let  be a probability measure on  and let (  ) be a sequence of discrete probability measures on  with   ⇀  weakly*.Each   has the form where  , are suitable points in  .Since  =  () we may take  , ∈  such that  ( , ) =  , and define a discrete probability measure   on  by Then  #   =   and, possibly passing to a subsequence, we may assume that   ⇀  weakly* for some probability measure  on .Passing now to the limit as  → ∞, we obtain  #  = .
The usual uniform Monte Carlo method consists in taking  as the Lebesgue measure on the source space .In this case it often happens that the image points are unevenly distributed in  , making the numerical identification of the image set  of a rather poor quality.On the contrary, our goal is to construct (a discrete approximation of) a measure  in order to obtain a well distributed image measure , as close as possible to the Lebesgue measure on  .The proof of Theorem 2.1 is clearly non-constructive, since the image set  =  () is not a priori known.We then need a constructive method that provides the probability measure  in the theorem above through an approximation procedure.This is the goal of next sections.
To further motivate our approach, let us recall the classical optimal transport problem raised by Monge [20]: Given probability measures ,  on  and a cost function If the cost is simply given by the Euclidean distance squared ( 1 ,  2 ) = ‖ 1 −  2 ‖ 2 , the infimium in (2.1) is called the Wasserstein 2-distance and is denoted by  2 2 (, ).Our objective is to approximate the Lebesgue measure in the image  =  () by a discrete measure, given by a sum of Dirac masses supported on points "uniformly" distributed in  =  ().Suppose that the measure  does not give mass to ( − 1)-dimensional sets and the measure  is discrete, given by a sum of Dirac masses

Centroidal Voronoi Tessellations
Consider  a compact connected region in R  (typically a rectangular box).Given  points   ∈ , consider the associated Voronoi diagram consisting of a partition (  )  =1 of  such that In other words,   contains all points in  which are closer to   compared to the remaining images (  ) ̸ = .Note that in our convention, Voronoi cells are bounded and are subsets of .The Voronoi cells are, in general, different in volume and are not necessarily uniform, for a general distribution of points.See the example in Figure 1 where a Voronoi diagram corresponding to 10 random points in the square [−1, 1] 2 is shown.The Voronoi points are represented with red squares and the centroids of the Voronoi cells are represented with blue points.There exist, however, particular classes of Voronoi diagrams which have cells that are more uniform in size, called Centroidal Voronoi Tessellations (CVT).For such diagrams, the point   which determines the Voronoi cell   coincides with the centroid of the region   .An example of CVT is shown in Figure 1 where it can be observed that the Voronoi points overlap with the centroids of the associated cells.The Voronoi regions for such a configuration have uniform sizes.Relevant usages of CVTs involve optimal quantization, data compression, optimal quadrature and mesh generation.We refer to [8,9,17,18] for more details and references.There exist iterative algorithms that produce CVTs starting from general Voronoi diagrams, the most basic one being Lloyd's algorithm, described as follows.
Lloyd's Algorithm.For a prescribed number of iterations  ≥ 1, successively replace the Voronoi points   with the centroids   of their corresponding Voronoi cells   .The algorithm was introduced in [18] and has a straightforward implementation.Convergence properties for Lloyd's algorithm are investigated in [10,22].While being easy to implement, this algorithm is not the most efficient from a practical point of view for reasons expressed in the following.
In [9,17] it is shown that the Centroidal Voronoi Tessellation are critical points for the energy Indeed, an immediate computation which can be found, for example, in [9] shows that where   are the centroids of   , for 1 ≤  ≤  .For a more direct identification of the gradient one may consider the following alternative formula: Notice that the quadratic term is constant and can be ignored, leading to a simplified energy function.
In view of the gradient formula (2.5), Lloyd's algorithm can be viewed as a gradient descent algorithm for minimizing the energy (2.4) with no need for a step size control [9].Acceleration techniques for improving the convergence of Lloyd's algorithm are considered in [8], where a Newton algorithm is discussed.This proves efficient for a small number of Voronoi cells, but is computationally heavy for large  .
It is generally agreed in the numerical optimization practice that gradient descent algorithms have slow convergence for ill-conditioned problems [21].For large scale problems, quasi-Newton methods like the low memory BFGS algorithms (lbfgs) Chapter 7 of [21] have better convergence properties.In [17] such a quasi-Newton method was used for minimizing (2.4), providing a faster method for computing CVTs.It is well known that CVTs are not necessarily unique, but on the other hand CVTs obtained through a variational process have improved stability properties as already underlined in [17].Thus an alternative way of finding CVTs is the following.
Variational CVT.Minimize numerically the energy (2.4) using a quasi-Newton algorithm.Benchmarks presented in [17] show that such an algorithm is more efficient than Lloyd's Algorithm.

Approximation of Blaschke-Santaló diagrams using CVTs
The objective of this article is to provide numerical tools which allow the approximation of Blaschke-Santaló diagrams, i.e., the image of a mapping where  ⊂ R  is the set of admissible parameters, containing upper and lower bounds and other constraints.
Given  a positive integer, consider a random set of samples  1 , . . .,   ∈  ⊂ R  .As underlined in the introduction, the images  ( 1 ), . . .,  (  ) are not necessarily uniformly distributed in the image  () ⊂ R  .Our goal is to find a choice of the points  1 , . . .,   in such a way that their images  ( 1 ), . . .,  (  ) are uniformly distributed inside the image set  ().
In the following, we assume  () is bounded and consider  a bounding box containing  () strictly in its interior.Denote by   =  (  ) the images for the initial sampling.We obviously have   ∈ .Consider now the Voronoi diagram associated to the points   defined by (2.3).Since our goal is to obtain a more uniform distribution of the images, we search points   =  (  ) which produce a Voronoi diagram that is as close as possible to a CVT.Inspired from the results recalled in Section 2.2 we propose two algorithms for approximating Blaschke-Santaló diagrams.
Lloyd algorithm with projection.Lloyd's algorithm is simple to implement, for general CVTs, simply replacing the points with the corresponding centroids.When dealing with Blaschke-Santaló diagrams one would like, for each sample   ,  = 1, . . .,  to replace it with another admissible sample x such that  (x  ) is the centroid of the Voronoi region associated to  (  ).We are faced with two issues: (a) Given a point  in the image  (), find a sample realizing , i.e., find  ∈  such that  () = .(b) Given a general point  ∈ R  which may not be in the image  (), find  ∈  such that  () is the closest to  in a sense to be defined.
Both aspects enumerated above can be covered using a single optimization problem: In cases where  is a compact set and  is at least of class  1 problem (2.7) admits solutions and efficient approximations can be found using standard numerical optimization algorithms.A quick computation shows that the gradient of (2.7) (squaring the norm for simplicity) is   ()( () − ).The numerical optimization algorithm will produce an approximation  of the solution which is either a critical point (a local minimizer is possible) or it saturates some of the bound constraints for the parameter set.Thus, a numerical solution  will verify one of the following: -  ()( () − ) = 0, in which case, if   () is of full rank, we have  () = .If   () is singular, we may have  () ̸ = , see Remark 2.3; -if the bounds on  are active at the optimum, the minimum in (2.7) is strictly positive and it is highly probable that  is a boundary point of  ().We are thus lead to the following natural algorithm.In our implementation and in the following, the norm considered in problem (2.7) is the Euclidean one.However other choices are possible and may give different behavior for the algorithm.Supposing Algorithm 1 converges for a given threshold  we obtain a configuration of samples (  )  =1 such that: -Whenever   is such that the centroid   of   belongs to  () we have ‖ (  ) −   ‖ < .
The drawbacks of Algorithm 1 are similar to the ones of Lloyd's algorithm compared to the variational CVT.
One may interpret Algorithm 1 as a fixed point or gradient descent algorithm with projection.Such algorithms can be improved using quasi-Newton methods as described in the following.
Variational CVT for Blasche-Santaló diagrams.Similar to [17] we formulate a minimization problem.We propose to minimize the composition of (2.4) with the parametrization (2.6).For a given number of samples  we consider the functional  : with  defined in (2.4).In practice, we minimize  using quasi-Newton methods with bounds and linear constraints characterizing the parameter set .This type of problems can easily be handled using available implementations (fmincon in Matlab, Knitro see [6]).Assuming the function  defined in (2.6) is differentiable, the derivatives of  can be expressed with where  (  ) ∈ R × is the Jacobian of  evaluated at   and   is the centroid of the Voronoi region   ,  = 1, . . .,  .We arrive at the following algorithm.The minimization of the functional (2.8) is straightforward if the Voronoi diagram associated to a set of points can be computed.We use the routine compute RVD developed following the results in [17] from the library Geogram.https://github.com/BrunoLevy/geogram The optimization is performed in Matlab/Julia using fmincon or the Artelys Knitro software in Algorithm 2. Details regarding the optimization procedure and more specific aspects regarding the problem at hand are shown in the next section.

Application I: algebraic functions
We start with an algebraic example which is easy to state, but quickly becomes challenging.For  ≥ 2 consider the space Sym  ([−1, 1]) of symmetric  ×  matrices with real entries in the interval [−1, 1].We apply our algorithm to the study of the (tr, det) diagram (tr denotes the trace, det denotes the determinant).More precisely, consider the application  : Our goal is to identify the image  (Sym  ([−1, 1])) for some particular choice of .
While for  = 2 a complete analytical description of the diagram is possible, for  ≥ 3 the problem becomes challenging.On the other hand, the numerical method we propose is efficient and shows a clear description of the corresponding diagram.
Analysis of the two dimensional case.We have the following result.and fix since it is a second degree polynomial, concave in , having a maximum at  =  /2 ∈   and its minimum at one of the endpoints of the interval   .If  ≤ 0 the minimum is attained at  =  + 1 and if  ≥ 0 the minimum is attained at  =  − 1.Moreover, it is obvious that at  fixed, the image of the determinant is connected, since it is the continuous image of an interval.The conclusion follows, observing that  2 ∈ [0, 1].
Illustration of the numerical algorithms.In the following we apply the algorithms proposed in Section 2 to study the proposed diagram.The bounding boxes  are considered as follows: We have the following observations: Next we apply Algorithm 2 for  = 200 samples.The constrained numerical optimization problem is executed using the software Artelys Knitro with the active-set algorithm.The initial samples are randomly chosen with values in [−1, 1]  .The number of iterations is limited to 1000 and the optimality criterion tolerance is set to 10 −8 .All simulations reached the maximum number of iterations under these constraints.Results are shown in Figure 3. Similar observations can be underlined, recalling again that the diagrams are approximated remarkably well with 200 well distributed samples.
The function (2.4) is evaluated for the optimal configuration of both algorithms and shown in Table 1.It is clearly observed that Algorithm 2 provides a lower energy since this algorithm is focused on decreasing the objective value.Moreover, the use of quasi-Newton descent directions (compared to anti-gradient descent direction for Lloyd) is known to accelerate the convergence.On the other hand Lloyd's algorithm can make large changes in the optimization variables (each sample is replaced by another one closest to the corresponding centroid as possible) which is advantageous when applied to initial random configuration or to regions where (2.4) varies too slowly.
Concerning the computational cost, Algorithm 1 is more costly (in our direct implementation), possibly due to the large number of small size optimization problems (2.7) that need to be solved at each iteration.Algorithm 2 simply computes one Restricted Voronoi Diagram per iteration (global or line-search) and evaluates the associated cost function and its gradient.Table 1.Comparison of the energies (2.4) after 1000 iterations of Algorithms 1, 2 for the final configurations shown in Figures 2, 3. Applying the proposed algorithms directly for a large number of samples is rather inefficient, convergence being slow.In the following we propose a multi-grid strategy that accelerates convergence.

Practical aspects: re-centering, multi-grid
In many numerical applications, multi grid strategies are employed to accelerate simulations.An initial simulation is performed on a coarse grid.Then the discretization is refined and the current numerical optimizer is interpolated on the new grid and used as an initialization for a new optimization procedure.The refinement procedure is repeated until the desired precision is attained.
We intend to use a similar strategy described below.Given a point  =  () ∈  () we search for other points (  )  =1 in  such that  (  ) is close to  () for  = 1, . . ., .A natural idea is to find a small enough -dimensional ellipsoid around  which is approximately mapped onto a -dimensional sphere around  ().However, if  is a boundary point for  () it is not possible to find such an ellipsoid.Nevertheless, if the dimension  of the space of parameters is larger than the dimension  of the image, the Jacobian  will be most likely non-singular, thus allowing us to move the point in the interior of the constraint set  while preserving the same image.This idea is detailed below.

Re-centering procedure
We start by giving the following result.Suppose  > , which is always the case in the applications.
of class  2 and  0 ∈ , such that  ( 0 ) is of full rank.Then, if the affine hyperplane  0 + ker   ( 0 ) is not tangent to , it is always possible to find  in the interior of  such that  () =  ( 0 ).
In practice, we minimize the distance between  and the center of the box  when  () =  ( 0 ): Generic available software like fmincon in Matlab allows the implementation of the nonlinear constraint  () =  ( 0 ).The power  is chosen large enough ( = 10 in practice) such that the maximal difference between the coordinates of  and coordinates of the center of  becomes as small as possible.Problems (4.1) are computationally cheap for the algebraic application proposed previously.

Multi-grid procedure
Algorithms 1 and 2 proposed in Section 2 can converge slowly when the number of samples is large.This is especially problematic when random samples tend to concentrate mostly in some particular regions of the Blaschke-Santaló diagram.Thus, we are interested in ways of enriching the set of samples given by Algorithms 1 or 2 for a rather small initial number of samples.We propose two refinement procedures below.
(a) Spheres around current samples.Suppose  > .For a point  0 such that  ( 0 ) is not singular, we re-center it using (4.1).We start by computing the singular value decomposition in S −1 and consider points   =  0 +    , where  > 0 is a radius small enough.In practice  is equal to one third of the minimal distance among images  (  ).Among points   we select only those that belong to .If  0 is a boundary point for  () it is possible that only a few of the points   ,  = 1, . . .,  are admissible.Whenever a refinement is necessary, the procedure described above is repeated for every sample point (  )

𝑖=1
for which the singular values of  (  ) are above a certain threshold (10 −3 in our implementation).In Figure 4 an example of application of the methods described above is shown.We take a result given by Algorithm 2. At the end of the optimization process some samples may have images at the boundary of .Applying directly the multi-grid procedure, trying to add points on a circle around the current samples gives the result shown in the left picture in Figure 4.For some points in the interior of  () the algorithm fails to add the required number of points since the associated samples are on the boundary of .In the right picture, re-centering is performed before applying the refinement procedure.For all interior points of  () the algorithm manages to add the prescribed number of additional samples with images close to the previous ones.It can be observed that for the upper boundary, corresponding to a singular Jacobian in this case, no points are added, since the procedure described above cannot be applied.If the minimization (4.2) succeeds, the solution   is a sample for which  () =   .If   does not belong to  () then the minimization problem still produces a sample, as close as possible to   .There are two issues that motivate us to solve (4.2) only for particular edges in  .
-When  () is non-convex, the Delaunay triangulation  will contain triangles which are not contained in  ().Usually, such triangles have an obtuse angle, some sides being significantly larger than the others.-Some regions may contain a denser concentration of samples than others.Therefore, for edges with length too small compared to the average, we choose not to add the corresponding midpoints to the diagram.
In practice we compute the average length ℓ of sides of triangles in  and we solve (4.2) for midpoints of edges with length in [0.5ℓ, 1.5ℓ].Figure 5 shows the outcome of this refinement procedure for the same test case as the one showin in Figure 4.

Global algorithm and numerical examples
Taking the previous considerations into account leads us to the following practical approach.In particular, we combine the advantages of Algorithms 1, 2 and the refinement strategy.Algorithm 2 is run until convergence criteria are met or the maximal number of iterations is reached.Initialization: Choose  random samples in .Run Algorithm 1 for  1 iterations followed by Algorithm 2.
For  = 1 to  ref do: -Multigrid: do one of the following.
• Apply the re-centering procedure, solving (4.1) for each one of the resulting samples.Apply the first refinement procedure adding at most  add points around each sample.• Alternatively, use the midpoints of the triangles in the Delaunay triangulation to find new sample points.-Run Algorithm 1 for  1 iterations.
Running a few iterations of Algorithm 1 before running Algorithm 2 is motivated by the fact that Lloyd's algorithm can make large jumps when applied to a non-optimal initial configuration.
In the following we show how Algorithm 3 approximates Blascke-Santaló diagrams in practice.First we apply it for the (tr, det) diagram in the case  = 2.We start with a set of 30 samples and we perform three refinements.The simulations have 30, 104, 464 and 2423 samples, respectively.Initialization and the results of the successive stages of the algorithm are shown in Figure 6.
Similar simulations are made for  ∈ {3, 4}.The resulting final optimized configurations having 2684 and 2043 cells, respectively, are shown in Figure 7.
In all results presented up to this point, the centroids of Voronoi cells are shown in blue and the images of the samples are shown in red.In some situations, for  ∈ {3, 4} we may observe inner Voronoi points which do not coincide with the Voronoi cell's centroid.This happens especially close to the curves parametrized by

Extracting the Blaschke-Santaló diagram from the Voronoi diagram
The boundary points of the Blaschke-Santaló diagram can be recovered selecting only samples for which the image is far from the centroid of the associated Voronoi cell.However, extracting a polygon from these points is not straightforward.We use the following ideas to plot the Blaschke-Santaló diagram starting from the numerical results: -If the resulting diagram is convex, taking the convex hull of the images of the samples suffices.
-Since the optimized samples form a Centroidal Voronoi tessellation, except the boundary points, we exploit the associated Delaunay triangulation which covers the entire convex hull of the diagram.However, triangles which are outside our diagram are nearly flat (having a small or large angle).We eliminate from the Delaunay triangulation such triangles (with thresholds that are set case by case).
The resulting diagrams are shown in Figure 8.

Subset of a diagram
In some cases, more detail is needed regarding certain parts of a Blaschke-Santaló diagram.Rather than increasing the point density everywhere in the diagram, it is possible to focus only on the region of interest.Suppose we are interested in the region  () ∩ (, ).Then we can implement all algorithms presented previously adding the additional constraint The practical difficulty is that constraints (4.3) are nonlinear.General software like fmincon and Knitro allow the use of non-linear constraints.The behavior of the optimization algorithm is improved if the gradient of the non-linear constraints is computed explicitly, which is possible in our case.Nevertheless, adding the non-linear constraints (4.3) slows down the proposed algorithms.The speed loss is compensated by a lower number of samples, since we only focus on a subset of the desired diagram.
As an example, we show the (tr, det) diagram for  = 4 contained in the disk ((3.9,1), 0.2).We were interested in exploring extremal matrices in this region of the diagram, since the boundary parametrization seemed to change here.See Figure 9.

Details concerning 𝑑 ∈ {3, 4}
We already saw in Proposition 3.1 the case  = 2 admits an explicit characterization of the boundary of the (tr, det) diagram.For higher dimensions such a description is difficult to obtain.Nevertheless, identifying the extremal matrices for diagrams given in Figure 7 we are able to conjecture a precise parametrization of the corresponding boundaries.
We give a brief description of the extremal matrices and corresponding parametrizations of the boundary for  = 3: -Top boundary part 3: ) Similar observations can be made for  = 4. Due to symmetry reasons, we only detail the right-half of the boundary: -Right boundary: -Bottom right: The resulting parametrized boundaries are shown in Figure 10.

Higher-dimensions
The algorithm proposed generalizes in a straightforward way to three dimensional diagrams.All algorithmic aspects remain the same.Three dimensional Restricted Voronoi diagrams are computed again using the library Geogram [17].The example for the diagram (tr,  1  2 +  2  3 +  3  1 , det) for matrices in Sym 3 ([−1, 1]) is shown in Figure 11.As usual,  1 ,  2 ,  3 denote the eigenvalues of the 3 × 3 matrix.

Comparison with Monte Carlo method
The simplest approach to investigate the Blaschke-Santaló diagrams is generating random samples and computing the corresponding images.As underlined in Section 2.1, this choice does not necessarily produce images uniformly distributed in the desired diagram.In the following, we generate progressively, a fixed number of sample points and the corresponding images.We compare the quality of the result and the computational cost with the algorithms proposed in the previous sections.
We consider 10 4 and, respectively 10 6 random samples in Sym  ([−1, 1]), evaluate (3.1) and plot the corresponding points in R 2 .The corresponding results for  ∈ {2, 3, 4} are shown in Figure 12.The Blaschke-Santaló diagrams computed with the algorithms proposed in previous sections are represented as polygons while random samples are represented by points.The diagrams are rescaled to have the same width in the horizontal and vertical directions.
We notice that the Monte Carlo approach is inefficient, especially when the dimension increases.For  ∈ {3, 4} using one million random samples fails to give an accurate description of the diagram.
In comparison, we give an analysis of the computational cost for some of our simulations which give a high quality approximation of the Blaschke-Santaló diagrams.Simulations in Figures 2 and 3 use  = 200 samples with a limit of 1000 iterations.For Algorithm 2 at most 200 × 1000 function and gradient evaluations are preformed.The computation of the Voronoi diagrams using Geogram is very efficient.Algorithm 1 performs additional function evaluations when projecting the centroids on the space of samples, but remains of the same order of magnitude: ( × ), where  is the number of iterations.

Application II: example from convex geometry
We focus now on an application from convex geometry.Various other works investigate inequalities between geometric quantities using Blaschke-Santaló diagrams.Among these we mention [7,[11][12][13][14].In order to apply directly our computational framework we consider a particular case where functionals involved are smooth and the corresponding diagram is bounded.Consider the following three quantities: area (Ω), perimeter Per(Ω), momentum of inertia  (Ω) among two dimensional convex shapes Ω with two axes of symmetry.Since in this case the centroid is at the origin, the momentum of inertia is given by  = ∫︀ Ω || 2 d.We use this symmetry constraint to be able to compare our results with the theoretical results in [15].A second motivation is related to the simplification of the numerical framework.It is likely that removing the constraint will change the resulting diagram.
One can notice that both scale invariant ratios considered above are maximized by the disk.Theoretical details regarding the corresponding Blaschke-Santaló diagram are studied in [15].
We consider the mapping where the factor 100 is added so that the two quantities are comparable.Our objective is to approximate the image of the mapping ℱ defined above.Various methods were developed for parametrizing convex sets.We mention intersections of hyperplanes [16], the support function parametrized using truncated Fourier series in [1] or values on a discrete grid in [3].Methods proposed previously are generally based on linear inequality constraints on the set of parameters.In order to apply our framework directly, a more direct parametrization, using only bound constraints would be more appropriate.This leads us to propose an alternate, yet classical, discretization process.We focus on the class of convex sets with two orthogonal axes of symmetry.Since we are also working in a scale invariant setting, it is enough to parametrize concave and decreasing functions  : [0, 1] → R. Given a uniform discretization of [0, 1] using  + 1 points, observe that if  0 , . . .,   are samples of a concave decreasing function at   = /,  = 0, . . .,  then: -the first order differences are   =   −  +1 ,  = 0, . . .,  − 1 are increasing -the second order differences  +1 =  +1 −   ,  = 0, . . .,  − 2 are non-negative.Conversely, given non-negative values   , it is possible to construct samples of a concave decreasing function having   as second order differences.Therefore we take (  ) −1 =1 ,  0 :=  0 and   as variables in our parametrization.
We immediately obtain the following equalities which show that   , for  = 0, . . ., −1, can be expressed in terms of   and (  ) −2 =0 using the following expression: with A = (A  ) given by A  =  + 1 − max{, } for 1 ≤ ,  ≤ .The coordinates of the boundary points of the discrete convex set are given by (  ,   ),  = 0, . . .,  for the first quadrant.They are symmetrized to obtain the rest of the boundary.The area and the perimeter are computed in a straightforward way.For the momentum of inertia, we use the explicit formulas for polygons, found for example in [25], a direct consequence of Green's formulas.Since all computations are analytic in terms of the parameters, the partial derivatives of all quantities of interest are also computed analytically.
(a) Randomly generated shapes.Given the parametrization above and a number of parameters  ≥ 2 we can generate random convex shapes and plot the points given by (5.1).We generate 1000 random shapes

Computational costs
In this section we discuss the computational cost of the proposed algorithms.Following simulations shown in Figures 12,14 it is obvious that the naive implementation of the Monte Carlo method is inefficient.For the (tr, det) diagram, 10 6 random samples do not give a relevant description of the image for  ≥ 3.For the shape optimization example, different choices for the number of parameters give different descriptions: when considering two parameters, the upper boundary is well captured.Indeed, in [15] it is proved that the straight line appearing in the upper boundary correspond to rhombi.However, increasing the number of parameters does not improve the description of the diagram in this case, as the images of the samples tend to concentrate around the origin.
In order to be able to reproduce our results we publish a version of the codes at the following Github repository: https://github.com/bbogo/BlaschkeSantalo The computation times below are reported for the fmincon optimization software in Matlab.Using more specialized optimization software like KNITRO further decreases the computation time.All computation times reported in the following correspond to a machine with CPU Intel Core i7-8650 (quad core, max 4.2Ghz) with 32 GB of RAM.
The computational cost associated to results presented in Figures 2, 3 is reported in Table 2.Note that initial samples were chosen randomly, therefore multiple runs of the same algorithm may yield different results.Multiple conclusions can be drawn: -Applying Algorithms 1 and 2 directly for a large number of samples is not necessary the most efficient strategy.For  = 2 the computational costs are rather small, but for  ≥ 3 the cost increases.-Lloyd's algorithm generally uses more iterations and function evaluations than CVT, but the overall cost depends on the optimization algorithm used in CVT (fmincon's active-set in this case).-CVT manages to capture correctly the Blaschke-Santaló diagram using significantly fewer function evaluations than the Monte Carlo simulations shown in Figure 12.Compared to Lloyd's algorithm, CVT solves a global optimization problem, rather than a series of small ones.Thus, if the space of parameters has large dimensions and the simulation is done for many samples, CVT may not be feasible (at least not in the current implementation).
Let us now illustrate an example of computation cost for running Algorithm 3 using 20 iterations in Lloyd's algorithm before applying CVT.We start with 20 random initial samples and perform two refinements.The results are summarized in Table 3.It can be noted that in both cases  ∈ {2, 3} the outcome of Algorithm 3 using refinement strategies is more efficient than running a single algorithm (Lloyd or CVT) straight away for a large number of random initial samples.The more precise results in Figures 6, 7 (using thousands of samples)  took up to one hour of computational time using the KNITRO algorithms.However, one may argue, in view of the aspects enumerated above, that an accurate description of the Blaschke-Santaló diagrams can be obtained efficiently with using only a few hundred samples.
For the geometric application presented in Section 5 the software KNITRO was used to obtain results in Figure 15.The active-set algorithm in fmincon solves approximately the Karush-Kuhn-Tucker optimality conditions corresponding to the active constraints using a direct solver.For problems of small size, like in the algebraic application, this poses no problem.When dealing with the geometric application, each sample corresponds to 50 variables.Thus, the final configuration, uses 516 × 50 = 25800 variables.KNITRO uses an iterative solver with the active-set method, making this problem solvable.Simulations made for Figure 15 took less than one hour.
Nevertheless, Algorithm 1 (using fmincon for the numerical optimization) coupled with the refinement strategies can be applied to find the diagram in Figure 15.It only involves solving a sequence of numerical optimization problems having the same dimension as the space of parameters.Computational details regarding this example are presented in Table 4.The final result obtained is shown in Figure 17 and is similar to the one shown in Figure 15.This suggests that more complex functionals can be investigated with the proposed methods, taking into account the additional computational cost coming more costly objective functions.This is underlined in the next section.

Extensions and perspectives
In view of the presentation of the computational costs in the previous section, a natural continuation of this work should investigate efficient algorithms for minimizing numerically problems of the form (2.8).The implementations which produced results presented in this paper used existing software which did not exploit the particular structure of the problem: the variables  1 , . . .,   are almost independent in the cost function in the sense that the Hessian matrix for (2.8) has a sparse structure.More precisely, each   interacts at the second order only with the samples corresponding to the neighboring cells in the Voronoi diagram [17].
From the point of view of applications, other geometric Blaschke-Santaló diagrams could be investigated, involving the eigenvalues of the Laplace operator, for example (like in [14]).In the example presented in Section 5 the objective function considered (depending on area, perimeter and moment of inertia) was explicit in terms of the parameters, yielding its computation almost instantaneous.For a functional defined using a partial differential equation, the computational cost will involve an additional non-negligible number of function evaluations, which should be added to the computational cost presented in the previous section.Note that in the case where the objective function is shape differentiable, its gradient with respect to the parameters defining the convex shape, like in Section 5, can be computed (see [3] for example).Therefore, the proposed framework applies.Using very precise and efficient methods (spectral methods, for example) and limiting the number of iterations in Lloyd's algorithm (1) can make the computation of such diagrams feasible using the methods proposed in this article.

Conclusions
We propose efficient algorithms which approximate Blaschke-Santaló diagrams by generating samples having uniformly distributed images.The key ingredient is the search for images which produce Centroidal Voronoi Tessellations.The algorithms proposed, inspired from Lloyd's algorithm and the Variational method proposed in [17] are illustrated through examples coming from linear algebra and shape optimization.
We observe that using a reasonable computational cost, compared with the usual Monte Carlo methods which generate randomized samples, the algorithms proposed achieve a precise description of the Blaschke-Santaló diagrams.Using a multi-grid strategy, more samples can be considered, further improving the description of these diagrams.

Figure 1 .
Figure 1.General Voronoi diagram (left) vs. Centroidal Voronoi Diagram (right).The points (  ) are depicted by squares, the centroids of the Voronoi regions are depicted by dots.

Figure 2 .
Figure 2. Approximations of the (tr, det) diagram in Sym  ([−1, 1]) for cases  = 2 (left) and  = 3 (right) using Algorithm 1. Red dots are the images  (  ), blue dots are the centroids of the Voronoi cells   .Red points which do not coincide with blue points generally correspond to boundary points for the image  ().

Figure 4 .
Figure 4. Refinement procedure: application before re-centering fails to add the required number of points around all interior points in the image.(left) After re-centering, the multigrid procedure succeeds for all interior points.(right) Blue dots are the original images and red dots are the images of the new samples.
(b) Delaunay Triangulations.These triangulations are closely related to Voronoi diagrams and their computation is standard in computational geometry.The Delaunay triangulation is the dual graph of the Voronoi diagram, the circumcenters of the Delaunay triangles being the vertices of the Voronoi diagram.Given a set of samples (  )  =1 and the corresponding images   =  (  ),  = 1, . . .,  , start by computing the Delaunay triangulation  of (  )  =1 .Next, select all edges of triangles in  .For each such edge, take the midpoint   and solve the problem min ∈ ‖ () −   ‖. (4.2)

Figure 5 .
Figure 5. Refinement strategy using projections of midpoints of certain edges of the triangles in the Delaunay triangulation.Blue dots represent the images of the original samples.Red dots represent the images of the new samples.

Figure 8 .
Figure 8. Extracting the Blaschke-Santaló diagram from the optimal Voronoi tesellation by eliminating "flat" triangles from the associated Delaunay triangulation.

Figure 13 .
Figure13.Parametrization of concave decreasing function using second order differences.The first order differences   are decreasing.

Figure 16 .
Figure 16.Examples of shapes corresponding to points on the boundary of the (100   2 ,  2  ) diagram.

Figure 17 .
Figure 17.The diagram for the (100   2 ,  2  ) diagram using only Lloyd's algorithm and the Delaunay refinement strategy.
Proposition 2.2.Assume the parameter set  ⊂ R  is compact and  is  1 .Suppose  * 1 , . . .,  *  minimizes (2.8) on   and denote by    , i.e.,  ( *  ) is the centroid of the region   .In practice, if Proposition 2.3 does not apply, we may have the following situations.(a) If  (  ) is a boundary point for  () then  (  ) is not necessarily equal to the centroid   of the region   .(b) Interior points of  () for which the Jacobian  is not of full rank may behave like boundary points previously described in (a).In this case, a gradient based algorithm used for minimizing (2.8) may have stationary points for which some of the samples  In the practical examples, the sets where  () is singular correspond to one dimensional curves, denoted by   .The uniformity of the sampling in  () is not affected on either side of the singular curves   .
The matrices  and  are unitary and  is diagonal, containing the singular values on the diagonal.Assuming  ( 0 ) is not singular, the diagonal values in  are non-zero.Moreover, if (  ) 1≤≤ , (  ) 1≤≤ are columns if  and  , respectively and   , 1 ≤  ≤  are the singular values, we have the decomposition  ( 0 ) = Denote by  the matrix containing the first  columns of  as columns and  the  ×  diagonal matrix containing 1/  on its diagonal.Consider the vectors   , 1 ≤  ≤  which are columns of  •  •  .Then these vectors satisfy  ( 0 )  =   ,  ≤  ≤  where (  ) =1,..., is the canonical basis.Consider the matrix  =  •  •  containing   ,  = 1, . . .,  as columns.For  ≥ 3 consider  uniformly distributed points (  )  =1

Table 4 .
Computation times for Algorithm 1(Lloyd), with the Delaunay refinement strategy, applied for the geometric example in Section 5.
N. Samples N. Iterations N. func.eval Computation time