Pluripotential Numerics

TL;DR

This paper develops numerical methods to approximate key quantities in Pluripotential Theory, including extremal functions, transfinite diameter, and equilibrium measures, using polynomial meshes and orthonormalization.

Contribution

It introduces new algorithms for approximating pluripotential quantities with proven convergence and demonstrates their effectiveness through numerical tests.

Findings

Convergence of transfinite diameter approximation.

Uniform convergence of extremal function approximation.

Successful numerical tests in simple cases.

Abstract

We introduce numerical methods for the approximation of the main (global) quantities in Pluripotential Theory as the \emph{extremal plurisubharmonic function} $V_E^*$ of a compact $\mathcal L$-regular set $E\subset \C^n$, its \emph{transfinite diameter} $\delta(E),$ and the \emph{pluripotential equilibrium measure} $\mu_E:=\ddcn{V_E^*}.$ The methods rely on the computation of a \emph{polynomial mesh} for $E$ and numerical orthonormalization of a suitable basis of polynomials. We prove the convergence of the approximation of $\delta(E)$ and the uniform convergence of our approximation to $V_E^*$ on all $\C^n;$ the convergence of the proposed approximation to $\mu_E$ follows. Our algorithms are based on the properties of polynomial meshes and Bernstein Markov measures. Numerical tests are presented for some simple cases with $E\subset \R^2$ to illustrate the performances of the proposed methods.

Figures

Equations

$${\Lambda }_k(z_0,\dots ,z_k):=\underset{z\in E}{\sup }\underset{f\in \mathcal{C}(E),f\ne 0}{\sup }\frac{|I_k[f](z)|}{\Vert f{\Vert }_E}\le (k+1)$$

$$\mathrm{cap}(E):=\exp \left(-\underset{\mu \in {\mathcal{M}}^1(E)}{\min }{\int }_E{\int }_E\log \frac{1}{|z-\zeta |}d\mu (z)d\mu (\zeta )\right)$$

$$\delta (E):=\underset{k}{\lim }{\left|\mathrm{Vdm}(z_0,\dots ,z_k)\right|}^{\frac{1}{\dim {\mathcal{P}}^k}},\text{ for }(z_0,\dots ,z_k)\text{ Fekete array}.$$

$$g_E(z,\infty ):=\underset{\zeta \to z}{\mathrm{limsup}}\left(\sup \left\{u(\zeta ):u\in \mathcal{L}(\mathbb{C}),u{|}_E\le 0\right\}\right).$$

$$g_E(z,\infty )=\underset{\zeta \to z}{\mathrm{limsup}}\left(\sup \left\{\frac{1}{\mathrm{deg}p}\log |p(\zeta )|,p\in \mathcal{P},\Vert p{\Vert }_E\le 1\right\}\right).$$

$$|p(z)|\le \Vert p{\Vert }_E\exp (\mathrm{deg}p{g}_E(z,\infty )),\forall p\in \mathcal{P}$$

$$\Delta g_E(z,\infty )={\mu }_E$$

$$V_E(\zeta ):=$$

$$V_E^{\ast }(z):=$$

$${\tilde{V}}_E(\zeta ):=$$

$${\tilde{V}}_E^{\ast }(z):=$$

$$\alpha \prec \beta \text{ if }$$

$$\bar{j}:=\min \{j\in \{1,2,\dots ,n\}:{\alpha }_j\ne {\beta }_j\}.$$

$$\alpha :\mathbb{N}⟶{\mathbb{N}}^n.$$

$$e:{\mathbb{N}}^n\to \mathcal{P}({\mathbb{C}}^n)$$

$$\alpha \mapsto e_{\alpha }(z):=z_1^{{\alpha }_1}\cdot z_2^{{\alpha }_2}\cdot \cdots \cdot z_n^{{\alpha }_n}$$

$${\mathcal{P}}^k({\mathbb{C}}^n)=\mathrm{span}\{e_i(z),1\le i\le N_k\}=:\mathrm{span}{\mathcal{M}}_n^k,$$

$$N_k:=\dim {\mathcal{P}}^k({\mathbb{C}}^n)=\left(\genfrac{}{}{0ex}{}{n+k}{n}\right).$$

$${\mathrm{Vdm}}_k({\mathbf{z}}_1,\dots {\mathbf{z}}_{N_k}):=\det [e_i({\mathbf{z}}_j){]}_{i,j=1,\dots ,N_k}.$$

$$\left|{\mathrm{Vdm}}_k({\mathbf{z}}_1,\dots {\mathbf{z}}_{N_k})\right|=\underset{{\mathbf{\zeta }}_1,\dots {\mathbf{\zeta }}_{N_k}\in E}{\max }\left|{\mathrm{Vdm}}_k({\mathbf{\zeta }}_1,\dots {\mathbf{\zeta }}_{N_k})\right|$$

$${\delta }_k(E):={\left(\underset{{\mathbf{\zeta }}_1,\dots {\mathbf{\zeta }}_{N_k}\in E}{\max }\left|{\mathrm{Vdm}}_k({\mathbf{\zeta }}_1,\dots {\mathbf{\zeta }}_{N_k})\right|\right)}^{\frac{n+1}{nk{N}_k}}$$

$$\delta (E):=\underset{k}{\lim }{\delta }_k(E).$$

$$\Vert p{\Vert }_E\le C_k\Vert p{\Vert }_{A_k},$$

$$\mathrm{Card}{A}_k\ge N_k:=\dim {\mathcal{P}}^k({\mathbb{C}}^n)=\left(\genfrac{}{}{0ex}{}{k+n}{k}\right)=\mathcal{O}(k^n)$$

$${\left|{\mathrm{Vdm}}_k(z_1^{(k)},\dots ,z_{N_k}^{(k)})\right|}^{\frac{n+1}{nk{N}_k}}\to \delta (E)$$

$$\underset{k}{\mathrm{limsup}}{\left(\underset{p\in {\mathcal{P}}^k\setminus \{0\}}{\sup }\frac{\Vert p{\Vert }_E}{\Vert p{\Vert }_{L_{\mu }^2}}\right)}^{1/k}\le 1.$$

$$B_k^{\mu }(z):=K_k^{\mu }(z,z)=\sum _{j=1}^{N_k}|q_j(z,\mu ){|}^2.$$

$$\frac{B_k^{\mu }}{N_k}\mu {\rightharpoonup }^{\ast }{\mu }_E,$$

$$\underset{k}{\mathrm{limsup}}\Vert B_k^{\mu }{\Vert }_E^{1/(2k)}\le 1.$$

$${\mu }_k:=\frac{1}{\mathrm{Card}{A}_k}\sum _{j=1}^{\mathrm{Card}{A}_k}{\delta }_{z_j^{(k)}},$$

$$B_k(\hat{z})=\sum _{j=1}^{N_k}{c}_j{q}_j(\hat{z},\mu ):=p(\hat{z}),c_j:=\overline{q_j}(\hat{z},\mu ).$$

$$\Vert p{\Vert }_E\le {\left(\sum _{j=1}^{N_k}|c_j{|}^2\right)}^{1/2}\underset{z\in E}{\max }{\left(\sum _{j=1}^{N_k}|q_j(z,\mu ){|}^2\right)}^{1/2}=B_k(\hat{z})=\Vert B_k{\Vert }_E.$$

Full text

The methods rely on the computation of a polynomial mesh for $E$ and numerical orthonormalization of a suitable basis of polynomials. We prove the convergence of the approximation of $\delta(E)$ and the uniform convergence of our approximation to $V_{E}^{*}$ on all $\mathbb{C}^{n};$ the convergence of the proposed approximation to $\mu_{E}$ follows. Our algorithms are based on the properties of polynomial meshes and Bernstein Markov measures.

Numerical tests are presented for some simple cases with $E\subset\mathbb{R}^{2}$ to illustrate the performances of the proposed methods.

1. 1 Introduction

2. 2 Preliminaries

3. 2.1 Pluripotential theory: some definitions

4. 2.2 Admissible meshes and Bernstein Markov measures

5. 3 Approximating the extremal function

6. 3.1 Theoretical results

7. 3.2 The SZEF and SZEF-BW algorithms

8. 3.2.1 SZEF Algorithm Implementation

9. 3.2.2 SZEF-BW Algorithm Implementation

10. 3.3 Numerical Tests of SZEF and SZEF-BW

11. 3.3.1 Test case 1

12. 3.3.2 Test case 2

13. 3.3.3 Test case 3

14. 4 Approximating the transfinite diameter

15. 4.1 Theoretical result

16. 4.2 Implementation of the TD-GD algorithm

17. 4.3 Numerical test of the TD-GD algorithm

18. 4.3.1 TD-GD test case 1: the unit ball in $\mathbb{R}^{2}$

19. 4.3.2 TD-GD test case 2: the unit simplex in $\mathbb{R}^{2}$

20. 5 Approximating the equilibrium measure

1. Introduction

Let $E\subset\mathbb{C}$ be a compact infinite set. Polynomial interpolation of holomorphic functions on $E$ and its asymptotic are intimately related with logarithmic potential theory, i.e., the study of subharmonic function of logarithmic growth on $\mathbb{C}$. This is a well established classical topic whose study goes back to Bernstein, Fekete, Leja, Szegö, Walsh and many others; we refer the reader to [37], [47] and [41] for an extensive treatment of the subject.

To study polynomial interpolation on a given compact set $E$ one introduces the Vandermonde determinant (usually with respect to the monomial basis) and, for any degree $k\in\mathbb{N},$ tries to maximize its modulus among the tuples of $k+1$ distinct points. Any such array of points is termed Fekte array of order $k$ for $E.$ A primary interest on Fekete arrays is that one immediately has the bound

[TABLE]

for the Lebesgue constant $\Lambda_{k}$ (i.e., the norm of the polynomial interpolation operator $I_{k}$) for any Fekete array of order $k$ for $E.$

On the other hand Fekete arrays provide the link of polynomial interpolation with potential theory. Indeed, the logarithmic energy $\operatorname*{cap}(E)$ of a unit charge distribution on $E$ at equilibrium

[TABLE]

turns out to coincide with certain asymptotic of the modulus of the Vandermonde determinant computed at any sequence of Fekete points and with respect to the monomial basis. Since the considered asymptotic is a geometric mean of mutual distances of Fekete points, it is termed transfinite diameter of $E$ and denoted by $\delta(E),$

[TABLE]

The Fundamental Theorem of Logarithmic Potential Theory asserts that $\delta(E)=\operatorname*{cap}(E)=\tau(E)$, where $\tau(E)$ is the Chebyshev constant of $E$ and is defined by means of asymptotic of certain normalized monic polynomials. Moreover, provided $\delta(E)\neq 0$, for any sequence of arrays having the same asymptotic of the Vandermonde determinants as Fekete points the sequence of uniform probability measures supported on such arrays converges weak star to the unique minimizer of the logarithmic energy minimization problem, that is the equilibrium measure $\mu_{E}$ of $E$.

The other fundamental object in this theory is the Green function with pole at infinity $g_{E}(\cdot,\infty)$ for the domain $\mathbb{C}\setminus\hat{E}$, where $\hat{E}$ is the polynomial hull of $E.$

[TABLE]

Here $\mathcal{L}(\mathbb{C})$ is the Lelong class of subharmonic function of logarithmic growth, i.e., $u(z)-\log|z|$ is bounded near infinity.

It turns out that, provided $\delta(E)\neq 0$, one has

[TABLE]

It follows that the Green function of $E$ is intimately related to polynomial inequalities and polynomial interpolation: for instance, one has the Bernstein Walsh Inequality

[TABLE]

and the Bernstein Walsh Theorem [47], that relates the rate of decrease of the error of best polynomial uniform approximation on $E$ of a function $f\in hol(\operatorname{int}E)\cap\mathscr{C}(E)$ to the possibility of extending $f$ holomorphically to a domain of the form $\{g_{E}(z,\infty)<c\}.$

Lastly, it is worth to recall that, using the fact that $\log|\cdot|$ is the fundamental solution of the Laplace operator in $\mathbb{C}$, it is possible to show that

[TABLE]

in the sense of distributions.

When we move from the complex plane to the case of $E\subset\mathbb{C}^{n}$, $n>1$, the situation becomes much more complicated. Indeed, one can still define Fekete points and look for asymptotic of Vandermonde determinants with respect to the graded lexicographically ordered monomial basis computed at these points, but this is no more related to a linear convex functional on the space of probability measures (the logarithmic energy) neither to a linear partial differential operator (the Laplacian) as in the case of $\mathbb{C}.$

During the last two decades (see for instance [25], [27]), a non linear potential theory in $\mathbb{C}^{n}$ has been developed: Pluripotential Theory is the study of plurisubharmonic functions, i.e., functions that are upper semicontinuous and subharmonic along each complex line. Plurisubharmonic functions in this setting enjoy the role of subharmonic functions in $\mathbb{C}$, while maximal plurisubharmonic functions replace harmonic ones; the geometric-analytical relation among the classes of functions being the same. It turned out that this theory is related to several branches of Complex Analysis and Approximation Theory, exactly as happens for Logarithmic Potential Theory in $\mathbb{C}$.

It was first conjectured by Leja that Vandermonde determinants computed at Fekete points should still have a limit and that the associated probability measures sequences should converge to some unique limit measure, even in the case $n>1.$ The existence of the asymptotic of Vandermonde determinants and its equivalence with a multidimensional analogue of the Chebyshev constant was proved by Zaharjuta [48], [49]. Finally the relation of Fekete points asymptotic with the equilibrium measure and the transfinite diameter in Pluripotential Theory (even in a much more general setting) has been explained by Berman Boucksom and Nymström very recently in a series of papers; [6], [5]. Indeed, the situation in the several complex variables setting is very close to the one of logarithmic potential theory, provided a suitable translation of the definitions, though the proof and the theory itself is much much more complicated.

Since Logarithmic Potential Theory has plenty of applications in Analysis, Approximation Theory and Physics, many numerical methods for computing approximations to Greens function, transfinite diameters and equilibrium measure has been developed following different approaches as Riemann Hilbert problem [31], numerical conformal mapping [24], linear or quadratic optimization [39, 38] and iterated function systems [28].

On the other hand, to the best authors’ knowledge, there are no algorithms for approximating the corresponding quantities in Pluripotential Theory; the aim of this paper is to start such a study. This is motivated by the growing interest that Pluripotential Theory is achieving in applications during the last years. We mention, among the others, the quest for nearly optimal sampling in least squares regression [30, 42, 22], random polynomials [14, 50] and estimation of approximation numbers of a given function [46].

Our approach, first presented in the doctoral dissertation [32, Part II Ch. 6], is based on certain sequences, first introduced by Calvi and Levenberg [21], of finite subsets of a given compact set termed admissible polynomial meshes having slow increasing cardinality and for which a sampling inequality for polynomials holds true. The core idea of the present paper is inspired by the analogy of sequences of uniform probability measures supported on an admissible mesh with the class of Bernstein Markov measures (see for instance [32], [15] and [8]).

Indeed, we use $L^{2}$ methods with respect to these sequences of measures: we can prove rigorously that our $L^{2}$ maximixation procedure leads to the same asymptotic as the $L^{\infty}$ maximization that appears in the definitions of the transfinite diameter (or other objects in Pluripotential Theory), this is due to the sampling property of admissible meshes. On the other hand the slow increasing cardinality of the admissible meshes guarantees that the complexity of the computations is not growing too fast.

We warn the reader that, though all proposed examples and tests are for real sets $E\subset\mathbb{R}^{n}\subset\mathbb{C}^{n}$, our results hold in the general case $E\subset\mathbb{C}^{n}.$ This choice has been made essentially for two reasons: first, the main examples for which we have analytical expression to compare our computation with are real, second, the case of $E\subset\mathbb{R}^{2}$ is both computationally less expensive and easier from the point of view of representing the obtained results.

The methods we are introducing in the present work are suitable to be extended in at least three directions. First, one may consider weighted pluripotential theory (see for instance [12]) instead of the classical one: the theoretical results we prove here can be recovered in such a more general setting by some modifications. It is worth to mention that a relevant part of the proofs of our results rest upon this weighted theory even if is not presented in such a framework, since we extensively use the results of the seminal paper [6]. However, in order to produce the same algorithms in the weighted framework, one needs to work with weighted polynomials, i.e., functions of the type $p(z)w^{\deg p}(z)$ for a given weight function $w$ and typically $E$ needs no more to be compact in this theory, these changes cause some theoretical difficulties in constructing suitable admissible meshes and may carry non trivial numerical issues as well.

Second, we recall that pluripotential theory has been developed on certain ”lower dimensional sets” of $\mathbb{C}^{n}$ as sub-manifolds and affine algebraic varieties. If $E$ is a compact subset of an algebraic subset $\mathcal{A}$ of $\mathbb{C}^{n}$, then one can extend the pluripotential theory of the set $\mathcal{A}_{reg}$ of regular points of $\mathcal{A}$ to the whole variety and use traces of global polynomials in $\mathbb{C}^{n}$ to recover the extremal plurisubharmonic function $V_{E}^{*}(z,\mathcal{A});$ see [40]. This last direction is probably even more attractive, due to the recent development of the theory itself especially when $E$ lies in the set of real points of $\mathcal{A},$ see for instance [29] and [7].

Lastly, we mention an application of our methods that is ready at hand. Very recently polynomial spaces with non-standard degree ordering (e.g., not total degree nor tensor degree) start to attract a certain attention in the framework of random sampling [22], Approximation Theory [46], and Pluripotential Theory [18]. For instance, one can consider spaces of polynomials of the form $\mathscr{P}^{k}_{q}:=\operatorname*{span}\{z^{\alpha},\alpha_{i}\in\mathbb{N}^{n},q(\alpha)<k\}$, where $q$ is any norm or even, more generally, $\mathscr{P}^{k}_{P}:=\operatorname*{span}\{z^{\alpha},\alpha_{i}\in\mathbb{N}^{n},\alpha\in kP\}$ for any $P\subset\mathbb{R}_{+}^{n}$ closed and star-shaped with respect to [math] such that $\cup_{k\in\mathbb{N}}kP=\mathbb{R}_{+}^{n}.$ Since this spaces are being used only very recently, many theoretical questions, from the pluripotential theory point of view, arise. Our methods can be used to investigate conjectures in this framework by a very minor modification of our algorithms.

The paper is structured as follows. In Section 2 we introduce admissible meshes and all the definitions and the tools we need from Pluripotential Theory.

Then we present our algorithms of approximation: for each of them we prove the convergence and we illustrate their implementation and their performances by some numerical tests; we stress that, despite the fact that we will consider only cases of $E\subset\mathbb{R}^{2}$ for relevance and simplicity, our techniques are fine for general $E\subset\mathbb{C}^{n}.$ We consider the extremal function $V_{E}^{*}$ (Pluripotential Theory counterpart of the Green function, see (3) below) in Section 3, the transfinite diameter $\delta(E)$ in Section 4 and the pluripotential equilibrium measure $\mu_{E}:=\left(\operatorname*{dd^{c}}V_{E}^{*}\right)^{n}$ in Section 5.

All the experiment are performed with the MATLAB software PPN package, see

http://www.math.unipd.it/$\sim$fpiazzon/software.

2. Preliminaries

2.1. Pluripotential theory: some definitions

Let $\Omega\subset\mathbb{C}$ be any domain and $u:\Omega\rightarrow\mathbb{R}\cup\{-\infty\}$, $u$ is said to be subharmonic if $u$ is upper semicontinuous and $u(z)\leq\frac{1}{2\pi r}\int_{|z-\zeta|=r}u(\zeta)ds(\zeta)$ for any $r>0$ and $z\in\Omega$ such that $B(z,r)\subseteq\Omega.$

A function $u:\Omega\rightarrow\mathbb{R}\cup\{-\infty\}$, where $\Omega\subset\mathbb{C}^{n}$, is said to be plurisubharmonic if $u$ is upper semicontinuous and is subharmonic along each complex line (i.e., each complex one dimensional affine variety); the class of such functions is usually denoted by $\mathcal{PSH}(\Omega).$ It is worth to stress that the class of plurisubharmonic functions is strictly smaller than the class of subharmonic function on $\Omega$ as a domain in $\mathbb{R}^{2n}.$

The Lelong class of plurisubharmonic function with logarithmic growth at infinity is denoted by $\mathcal{L}(\mathbb{C}^{n})$ and $u\in\mathcal{L}(\mathbb{C}^{n})$ iff $u\in\mathcal{PSH}(\mathbb{C}^{n})$ is a locally bounded function such that $u(z)-\log|z|$ is bounded near infinity.

Let $E\subset\mathbb{C}^{n}$ be a compact set. The extremal function $V_{E}^{*}$ (also termed pluricomplex Green function) of $E$ is defined mimicking one of the possible definitions in $\mathbb{C}$ of the Green function with pole at infinity; see (1).

[TABLE]

It turns out that the extremal function enjoys the same relation with polynomials of the Green function; precisely Siciak introduced [44]

[TABLE]

and shown (the general statement has been proved by Zaharjuta) that $\tilde{V}_{E}^{*}\equiv V_{E}^{*}$ and $\tilde{V}_{E}=_{\text{q.e.}}V_{E}.$ Here q.e. stands for quasi everywhere and means for each $z\in\mathbb{C}^{n}\setminus P$ where $P$ is a pluripolar set, i.e., a subset of the $\{-\infty\}$ level set of a plurisubharmonic function not identically $-\infty$. In the case that $V_{E}^{*}$ is also continuous, the set $E$ is termed $\mathcal{L}$-regular.

A remarkable consequence of this equivalence is that the Bernstein Walsh Inequality (3) and Theorem (see [44]) hold even in the several complex variable setting simply replacing $g_{E}(z,\infty)$ by $V_{E}^{*}(z)$; see for instance [27].

In Pluripotential Theory the role of the Laplace operator is played by the complex Monge Ampere operator $\left(\operatorname*{dd^{c}}\right)^{n}.$ Here $\operatorname*{dd^{c}}=2i\bar{\partial}\partial$, where $\partial u(z):=\sum_{j=1}^{n}\frac{\partial u(z)}{\partial z_{j}}dz_{j}$, $\bar{\partial}u(z):=\sum_{j=1}^{n}\frac{\partial u(z)}{\partial\bar{z}_{j}}d\bar{z}_{j}$ and $\left(\operatorname*{dd^{c}}u\right)^{n}=\operatorname*{dd^{c}}u\wedge\operatorname*{dd^{c}}u\dots\wedge\operatorname*{dd^{c}}u.$ For $u\in\mathscr{C}^{2}(\mathbb{C}^{n},\mathbb{R})$ one has $\left(\operatorname*{dd^{c}}u\right)^{n}=c_{n}\det[\partial^{2}u/\partial z_{i}\partial\bar{z}_{j}]d\text{Vol}_{\mathbb{C}^{n}}.$ The Monge Ampere operator extends to locally bounded plurisubharmonic functions as shown by Bedford and Taylor [4], [3], $\left(\operatorname*{dd^{c}}u\right)^{n}$ being a positive Borel measure.

The equation $\left(\operatorname*{dd^{c}}u\right)^{n}=0$ on a open set $\Omega$ (in the sense of currents) characterize the maximal plurisubharmonic functions; recall that $u$ is maximal if for any open set $\Omega^{\prime}\subset\subset\Omega$ and any $v\in\mathcal{PSH}(\Omega)$ such that $v|_{\partial\Omega^{\prime}}\leq u|_{\partial\Omega^{\prime}}$ we have $v(z)\leq u(z)$ for any $z\in\Omega^{\prime}.$

Given a compact set $E\subset\mathbb{C}^{n},$ two situations may occur: either $V_{E}^{*}\equiv+\infty$, or $V_{E}^{*}$ is a locally bounded plurisubharmonic function. The first case is when $E$ is pluripolar, it is too small for pluripotential theory. In the latter case $\left(\operatorname*{dd^{c}}V_{E}^{*}\right)^{n}=0$ in $\mathbb{C}^{n}\setminus E$ (i.e., $V_{E}^{*}$ is maximal on such a set), in other words the positive measure $\left(\operatorname*{dd^{c}}V_{E}^{*}\right)^{n}$ is supported on $E$. Such a measure is usually denoted by $\mu_{E}$ and termed the (pluripotential) equilibrium measure of $E$ by analogy with the one dimensional case.

Let us introduce the graded lexicographical strict order $\prec$ on $\mathbb{N}^{n}.$ For any $\alpha,\beta\in\mathbb{N}^{n}$ we have

[TABLE]

This is clearly a total (strict) well-order on $\mathbb{N}^{n}$ and thus it induces a bijective map

[TABLE]

On the other hand the map

[TABLE]

is a isomorphism having the property that $\deg e_{\alpha}(z)=|\alpha|.$

Thus, if we denote by $e_{i}(z)$ the $i$-th monomial function $e_{\alpha(i)}(z)$, we have

[TABLE]

where

[TABLE]

From now on we will refer to $\mathcal{M}_{n}^{k}$ as the graded lexicographically ordered monomial basis of degree $k.$

For any array of points $\{\boldsymbol{z}_{1},\dots\boldsymbol{z}_{N_{k}}\}\subset E^{N_{k}}$ we introduce the Vandermonde determinant of order $k$ as

[TABLE]

If a set $\{\boldsymbol{z}_{1},\dots,\boldsymbol{z}_{N_{k}}\}$ of points of $E$ satisfies

[TABLE]

it is said to be a array of Fekete points for $E.$ Clearly Fekete points do not need to be unique. Zaharjuta proved in his seminal work [48] that the sequence

[TABLE]

does have limit and defined, by analogy with the case $n=1$, the transfinite diameter of $E$ as

[TABLE]

It turns out that the condition $\delta(E)=0$ characterize pluripolar subsets of $\mathbb{C}^{n}$ as it characterize polar subset of $\mathbb{C}.$

Let $\{\boldsymbol{z}^{(k)}\}_{k\in\mathbb{N}}:=\{(\boldsymbol{z}_{1}^{(k)},\dots,\boldsymbol{z}_{N_{k}}^{(k)})\}_{k\in\mathbb{N}}$ be a sequence of Fekete points for the compact set $E$, we consider the canonically associated sequence of uniform probability measures $\nu_{k}:=\sum_{j=1}^{N_{k}}\frac{1}{N_{k}}\delta_{\boldsymbol{z}_{j}^{(k)}}.$

Berman and Boucksom [5] showed111The results of [5] and [6] hold indeed in the much more general setting of high powers of a line bundle on a complex manifold. that the sequence $\nu_{k}$ converges weak star to the pluripotential equilibrium measure $\mu_{E}$ as it happens in the case $n=1.$ We will use both this result (in Section 5) and a remarkable intermediate step of its proof (in Section 3 and Section 4) termed Bergman Asymptotic, see (21) below.

2.2. Admissible meshes and Bernstein Markov measures

We recall that a compact set $E\subset\mathbb{R}^{n}$ (or $\mathbb{C}^{n}$) is said to be polynomial determining if any polynomial vanishing on $E$ is necessarily the null polynomial.

Let us consider a polynomial determining compact set $E\subset\mathbb{R}^{n}$ (or $\mathbb{C}^{n}$) and let $A_{k}$ be a subset of $E$. If there exists a positive constant $C_{k}$ such that for any polynomial $p\in\mathscr{P}^{k}(\mathbb{C}^{n})$ the following inequality holds

[TABLE]

then $A_{k}$ is said to be a norming set for $\mathscr{P}^{k}(\mathbb{C}^{n}).$

Let $\{A_{k}\}$ be a sequence of norming sets for $\mathscr{P}^{k}(\mathbb{C}^{n})$ with constants $\{C_{k}\}$, suppose that both $C_{k}$ and $\operatorname*{Card}(A_{k})$ grow at most polynomially with $k$ (i.e., $\max\{C_{k},\operatorname*{Card}(A_{k})\}$ $=\mathcal{O}(k^{s})$ for a suitable $s\in\mathbb{N}$), then $\{A_{k}\}$ is said to be a weakly admissible mesh (WAM) for $E$; see222The original definition in [21] is actually slightly weaker (sub-exponential growth instead of polynomial growth is allowed), here we prefer to use the present one which is now the most common in the literature. [21]. Observe that necessarily

[TABLE]

since a (W)AM $A_{k}$ is $\mathscr{P}^{k}(\mathbb{C}^{n})$-determining by definition.

If $C_{k}\leq C$ $\forall k$, then $\{A_{k}\}_{\mathbb{N}}$ is said an admissible mesh (AM) for $E$; in the sequel, with a little abuse of notation, we term (weakly) admissible mesh not only the whole sequence but also its $k$-th element $A_{k}$. When $\operatorname*{Card}(A_{k})=O(k^{n})$, following Kroó [26], we refer to $\{A_{k}\}$ as an optimal admissible mesh, since this grow rate for the cardinality is the minimal one in view of equation (14).

Weakly admissible meshes enjoy some nice properties that can be also used together with classical polynomial inequalities to construct such sets. For instance, WAMs are stable under affine mappings, unions and cartesian products and well behave under polynomial mappings. Moreover any sequence of interpolation nodes whose Lebesgue constant is growing sub-exponentially is a WAM.

It is worth to recall other nice properties of (weakly) admissible meshes. Namely, they enjoy a stability property under smooth mapping and small perturbations both of $E$ and $A_{k}$ itself; [34]. For a survey on WAMs we refer the reader to [16].

Weakly admissible meshes are related to Fekete points. For instance assume a Fekete triangular array $\{\boldsymbol{z}^{(k)}\}=\{(z_{1}^{(k)},\dots,z_{N_{k}}^{(k)})\}$ for $E$ is known, then setting $A_{k}:=\boldsymbol{z}^{(k\log k)}$ for all $k\in\mathbb{N}$ we obtain an admissible mesh for $E$; [11].

Conversely, if we start with an admissible mesh $\{A_{k}\}$ for $E$ it has been proved in [19] that it is possible to extract (by numerical linear algebra) a set $\boldsymbol{z}^{(k)}:=\{z_{1}^{(k)},\dots,z_{N_{k}}^{(k)}\}\subset A_{k}$ from each $A_{k}$ such that the sequence $\{\boldsymbol{z}^{(k)}\}$ is an asymptotically Fekete sequence of arrays, i.e.,

[TABLE]

as happens for Fekete points. By the deep result of Berman and Boucksom [6] (and some further refining, see [10]) it follows that the sequence of uniform probability measures $\{\nu_{k}\}$ canonically associated to $\boldsymbol{z}^{(k)}$ converges weak star to the pluripotential equilibrium measure.

In [6] authors pointed out, among other deep facts, the relevance of a class of measures for which a strong comparability of uniform and $L^{2}$ norms of polynomials holds, they termed such measures Bernstein Markov measures. Precisely, a Borel finite measure $\mu$ with support $S_{\mu}\subseteq E$ is said to be a Bernstein Markov measure for $E$ if we have

[TABLE]

Let us denote by $\{q_{j}(z,\mu)\}$, $j=1,\dots,N_{k}$ the orthonormal basis (obtained by Gram-Schmidt orthonormalizaion starting by $\mathcal{M}_{n}^{k}$) of the space $\mathscr{P}^{k}$ endowed by the scalar product of $L^{2}_{\mu}.$ The reproducing kernel of such a space is $K_{k}^{\mu}(z,\bar{z}):=\sum_{j=1}^{N_{k}}\bar{q_{j}}(z,\mu)q_{j}(z,\mu)$, we consider the related Bergman function

[TABLE]

As a side product of the proof of the asymptotic of Fekete points Berman Boucksom and Nyström deduces the so called Bergman Asymptotic

[TABLE]

for any positive Borel measure $\mu$ with support on $E$ and satisfying the Bernstein Markov property.

Note that, by Parseval Identity, the property (16) above can be rewritten as

[TABLE]

Bernstein Markov measures are very close to the Reg class defined (in the complex plane) by Stahl and Totik and studied in $\mathbb{C}^{n}$ by Bloom [8]. Recently Bernstein Markov measures have been studied by different authors, we refer the reader to [15] for a survey on their properties and applications.

Our methods in the next sections rely on the fact admissible meshes are good discrete models of Bernstein Markov measures; let us illustrate this. From now on

• •

we assume $E\subset\mathbb{C}^{n}$ to be a compact $\mathcal{L}$-regular set and hence polynomial determining,

• •

we denote by $\mu_{k}$ the uniform probability measure supported on $A_{k}=\{z_{1}^{(k)},\dots,z_{M_{k}}^{(k)}\}$, i.e.,

[TABLE]

• •

we denote by $B_{k}(z)$ the function $B_{k}^{\mu_{k}}(z)$ and by $K_{k}(z,\zeta)$ the function $K_{k}^{\mu_{k}}(z,\zeta).$

Assume that an admissible mesh $\{A_{k}\}$ of constant $C$ for the compact polynomial determining set $E\subset\mathbb{C}^{n}$ is given. Now pick $\hat{z}\in E$ such that $B_{k}(\hat{z})=\max_{E}B_{k},$ we note that

[TABLE]

By Parseval inequality we have

[TABLE]

Therefore we an write

[TABLE]

thus

[TABLE]

On the other hand for any polynomial $p\in\mathscr{P}^{k}$ we have

[TABLE]

Recall that it follows by the definition of (weakly) admissible meshes that $(C^{2}\operatorname*{Card}A_{k})^{1/(2k)}\to 1.$

Thus, the sequence of probability measures associated to the mesh has the property

[TABLE]

which closely resembles (18).

Conversely, assume $\{\mu_{k}\}$ to be a sequence of probabilities on $E$ with $\operatorname*{Card}\operatorname*{supp}\mu_{k}=\mathcal{O}(k^{s})$ for some $s$, then we have

[TABLE]

Therefore, if $\|B_{k}\|_{E}=\mathcal{O}(k^{t})$ for some $t$, the sequence of sets $\{\operatorname*{supp}\mu_{k}\}$ is a weakly admissible mesh for $E$.

3. Approximating the extremal function

3.1. Theoretical results

In this section we introduce certain sequences of functions, namely $u_{k},v_{k},\tilde{u}_{k}$ and $\tilde{v}_{k}$, that can be constructed starting by a weakly admissible mesh, all of them having the property of local uniform convergence to $V_{E}^{*},$ provided $E$ is $\mathcal{L}$-regular.

Theorem 3.1**.**

Let $E\subset\mathbb{C}^{n}$ be a compact $\mathcal{L}$-regular set and $\{A_{k}\}$ a weakly admissible mesh for $E$, then, uniformly in $\mathbb{C}^{n}$, we have

[TABLE]

Proof.

We first prove (21), for we introduce

[TABLE]

The sequence of function $\Phi_{E}^{(k)}$ has been defined by Siciak and has been shown to converge to $\exp\tilde{V}_{E}^{*}$ (see equation (7)) for $E$ $\mathcal{L}$-regular, moreover we have $V_{E}^{*}\equiv\tilde{V}_{E}^{*};$ [44], see also [43].

Let us denote by $\mathcal{F}^{(k)}_{2}$ the family $\{p\in\mathscr{P}^{k}:\|p\|_{L^{2}_{\mu_{k}}}\leq 1\}$ we notice that, due to the Parseval Identity, we have

[TABLE]

Let us pick $p\in\mathcal{F}^{(k)}_{2}$, we have $\|p\|_{E}\leq\sqrt{\|B_{k}\|_{E}}\|p\|_{L^{2}_{\mu_{k}}}$ for the reason above, thus $q:=p\|B_{k}\|_{E}^{-1/2}\in\mathcal{F}_{E}^{(k)}.$

Hence

[TABLE]

It follows that

[TABLE]

On the other hand, since $\mu_{k}$ is a probability measure, we have $\|p\|_{E}\geq\|p\|_{L^{2}_{\mu_{k}}}$ for any polynomial. Hence if $p\in\mathcal{F}_{E}^{(k)}$ it follows that $p\in\mathcal{F}_{2}^{(k)}.$ Thus $v_{k}(z)\geq\log\Phi_{E}^{(k)}(z).$ Therefore we have

[TABLE]

Note that we have $\limsup_{k}\|B_{k}\|_{E}^{1/2k}\leq 1$ since $\{A_{k}\}$ is weakly admissible (see equation (20)), hence we can conclude that locally uniformly we have

[TABLE]

This concludes the proof of (21), let us prove (22).

It follows by Cauchy-Schwarz and Holder Inequalities and by $\int B_{k}d\mu_{k}=N_{k}$ that

[TABLE]

Thus it follows that

[TABLE]

On the other hand, for any $p\in\mathscr{P}^{k}$ we have

[TABLE]

hence, using the definition of Siciak function,

[TABLE]

Finally, using (23), we have

[TABLE]

uniformly in $\mathbb{C}^{n},$ this concludes the proof of (22) since $N_{k}^{1/(2k)}\to 1$ and both $v_{k}$ and $\log\Phi_{E}^{(k)}$ converge to $V_{E}^{*}$ uniformly in $\mathbb{C}^{n}.$ ∎

It is worth to notice that both $u_{k}$ and $v_{k}$ are defined in terms of orthonormal polynomials with respect to $\mu_{k}$, hence they can be computed with a finite number of operations at any point $z\in\mathbb{C}^{n}$, indeed we have

[TABLE]

Also we note that Theorem 3.1 can be understood as a generalization of the original Siciak statement [44, Th. 4.12]. Indeed, if we take $A_{k}=\{z_{1}^{(k)},\dots,z_{N_{k}}^{(k)}\}$ a set of Fekete points of order $k$ for $E$ we get $q_{j}(z,\mu_{k})=\sqrt{N_{k}}\ell_{j,k}(z)$, where $\ell_{j,k}(z)$ is the $j$-th Lagrange polynomial, hence we have

[TABLE]

Here $\Lambda_{A_{k}}(z)$ is the Lebesgue function of the interpolation points $A_{k}.$ Therefore, for $A_{k}$ being a Fekete array of order $k$, we have $u_{k}(z)=\log\left(\Lambda_{A_{k}}(z)\right)^{1/k}$, this is precisely $\exp\left(k\Phi_{k}^{(2)}(z)\right)$ in the Siciak notation.

In Section 5 we will deal with measures of the form $\nu_{k}:=\frac{B_{k}^{\mu}}{N_{k}}\mu$ for a Bernstein Markov measure $\mu$ for $E$, or, more generally $\nu_{k}:=\frac{B_{k}^{\mu_{k}}}{N_{k}}\mu_{k}$, where the sequence $\{\mu_{k}\}$ has the property $\limsup_{k}\|B_{k}^{\mu_{k}}\|_{E}^{1/2k}=1$; we refer to such a sequence $\{\mu_{k}\}$ as a Bernstein Markov sequence of measures. Due to a modification of the Berman Boucksom and Nymstrom result, $\nu_{k}$ converges weak star to $\mu_{E}$ (see Proposition 5.1 below). Note that $\nu_{k}$ is still a probability measure since $B_{k}(z)\geq 0$ for any $z\in\mathbb{C}^{n}$ and

[TABLE]

Here we point out another (easier, but very useful to our aims) property of the ”weighted” sequence $\nu_{k}$: actually they are a Bernstein Markov sequence of measure, more precisely the following theorem holds.

Theorem 3.2**.**

Let $E\subset\mathbb{C}^{n}$ be a compact $\mathcal{L}$-regular set and $\{A_{k}\}$ a weakly admissible mesh for $E.$ Let us set $\tilde{\mu}_{k}:=\frac{B_{k}^{\mu_{k}}}{N}\mu_{k},$ where $\mu_{k}$ is the uniform probability measure on $A_{k}.$ Then the following holds.

1. i)

For any $k\in\mathbb{N}$ and any $z\in\mathbb{C}^{n}$

[TABLE]

Thus $\limsup_{k}\|B_{k}^{\tilde{\mu}_{k}}\|_{E}^{1/2k}=1.$ 2. ii)

We have

[TABLE]

uniformly in $\mathbb{C}^{n}.$

From now on we use the notations

[TABLE]

where $B_{k}$ is as above $B_{k}^{\mu_{k}}$ and $\mu_{k}$ will be clarified by the context.

Proof.

We prove (25), then $\limsup_{k}\|\tilde{B}_{k}\|_{E}^{1/2k}=1$ follows immediately by $\limsup_{k}\|B_{k}\|_{E}^{1/2k}=1$ and $\lim_{k}N_{k}^{1/k}=1.$ The proof of (26) and (27) is identical to the ones of Theorem 3.1 so we do not repeat them.

We simply notice that, for any sequence of polynomials $\{p_{k}\}$ with $\deg p_{k}\leq k$, we have

[TABLE]

Now, for any $z\in E$, we pick a sequence $\{p_{k}\}$ such that it maximizes (for any $k$) the ratio $(|q(z)|\,\|q\|_{L^{2}_{\tilde{\mu}_{k}}}^{-1})^{1/k}$ among $q\in\mathscr{P}^{k}$ and we get

[TABLE]

Here the last inequality follows by the definition of $B_{k}$. Note in particular that $B_{k}(z)=1+|q_{2}(z)|^{2}+\dots$, thus $\frac{N}{\min_{E}B_{k}}\leq N=\mathcal{O}(k^{n})$ and $\|\tilde{B}_{k}\|_{E}^{1/2k}\leq N^{1/2k}\|B_{k}\|_{E}^{1/2k}\sim\|B_{k}\|_{E}^{1/2k}$ as $k\to\infty.$ ∎

Remark 3.3**.**

We stress that the upper bound (25) is in many cases quite rough, though sufficient to prove the convergence result (26). Indeed, since $\int B_{k}d\mu_{k}=1$ for any $k$, it follows that $\frac{N}{\min_{E}B_{k}}$ is always larger than $1$, but we warn the reader that $\|\tilde{B}_{k}\|_{E}$ does not need to be larger than $\|B_{k}\|_{E}$ in general. Hence the measure $\tilde{\mu}_{k}$ may be more suitable than $\mu_{k}$ for our approximation purposes.

3.2. The SZEF and SZEF-BW algorithms

The function $V_{E}^{*}$, at least for a regular set $E$, can be characterized as the unique continuous solution of the following problem

[TABLE]

It is rather clear that writing a pseudo-spectral or a finite differences scheme for such a problem is a highly non trivial task, as one needs to deal with a unbounded computational domain $\mathbb{C}^{n}\setminus E$, with a positivity constraint on $\left(\operatorname*{dd^{c}}u\right)^{n}$, and with a prescribed growth rate at infinity (both encoded by $u\in\mathcal{L}(\mathbb{C}^{n})$).

Here we present the SZEF and SZEF-BW algorithms (which stands for for Siciak Zaharjuta Extremal Function and Siciak Zaharjuta Extremal Function by Bergman weight) to compute the values of the functions $u_{k}$ and $v_{k}$ (see Theorem 3.1) and the functions $\tilde{u}_{k}$ and $\tilde{v}_{k}$ (see Theorem 3.2) respectively at a given set of points. In our methods both the growth rate and the plurisubharmonicity are encoded in the particular structure of the approximated solutions $u_{k}$ or $v_{k}$, while the unboundedness of $\mathbb{C}^{n}\setminus E$ does not carry any issue since all the sampling points used to build the solutions lie on $E.$ Indeed, once the approximated solution is computed on a set of points and the necessary matrices are stored, it is possible to compute $u_{k}$ or $v_{k}$ on another set of points by few very fast matrix operations; this will be more clear in a while.

To implement our algorithms we make the following assumptions.

• •

Let $E\subset\mathbb{C}^{n}$ be a compact regular set, for simplicity let us assume $E$ to be a real body (i.e., the closure of a bounded domain), but notice that this assumption is not restrictive neither from the theoretical nor from the computational point of view.

• •

We further assume that we are able to compute a weakly admissible mesh $\{A_{k}\}=\{z_{1},\dots,z_{M_{k}}\}$ for $E$ with constants $C_{k},$ $k\in\mathbb{N}.$ Note that an algorithmic construction of an admissible mesh is available in the literature for several classes of sets [36, 35, 33], since the study of (weakly) admissible meshes is attracting certain interest during last years.

• •

We assume $n=2$ to hold the algorithm complexity growth.

The implementation is based on the following choices.

• •

Let us fix a computational grid

[TABLE]

with finite cardinality $L:=L_{1}\cdot L_{2}$ and let us denote by $\Omega_{E}$ the (possibly empty) set $\Omega\cap E$, while by $\Omega_{0}$ the set $\Omega\setminus\Omega_{E}.$ We will reconstruct the values $u_{k}(\zeta),v_{k}(\zeta)$ $\forall\zeta\in\Omega$, however we will test the performance of our algorithm (see Subsection 3.3 below) only in term of error on $\Omega_{0}.$ This choice is motivated as follows. First, we have to mention the fact that the point-wise error of $u_{k}$ and $v_{k}$ exhibits two rather different behaviours depending whether the considered point $\zeta$ lies on $E$ or not: the convergence for $\zeta\in E$ is much slower. In contrast, for any regular set $E,$ the function $V_{E}^{*}$ identically vanishes on $E$, hence there is no point in trying to approximate it on $E.$ Note that the function $V_{E}^{*}(\zeta)\neq 0$ for any $\zeta\in\mathbb{C}^{n}\setminus E$, thus we can measure the error of $u_{k}$ of $v_{k}$ on $\Omega_{0}\subset\mathbb{C}^{2}\setminus E$ both in the absolute and in the relative sense.

• •

Let $L:\mathbb{C}^{2}\rightarrow\mathbb{C}^{2}$ be the invertible affine map, mapping $A_{k}$ in the square $[-1,1]^{2}$ and defined by $P_{i}(z):=\frac{2}{b_{i}-a_{i}}\left(z_{i}-\frac{a_{i}+b_{i}}{2}\right),$ $a_{i}:=\min_{z\in A_{k}}z_{i},$ $b_{i}:=\max_{z\in A_{k}}z_{i}.$ We denote by $T_{j}(z)$ the classical $j$-th Chebyshev polynomial and we set

[TABLE]

where $\alpha:\mathbb{N}\to\mathbb{N}^{2}$ is the one defined in (11). The set $\mathcal{T}^{k}_{P}:=\{\phi_{j}(z),\,1\leq j\leq(k+1)(k+2)/2\}$ is a (adapted Chebyshev) basis of $\mathscr{P}^{k}(\mathbb{C}^{2}).$ We will use the basis $\mathcal{T}^{k}_{P}$ for computing and orthogonalizing the Vandermonde matrix of degree $k$ computed at $A_{k}$. This choice has been already fruitfully used, for instance in [19, 16], when stable computations with Vandermonde matrices are needed and is on the background of the widely used matlab package ChebFun2 [45].

3.2.1. SZEF Algorithm Implementation

The first step of SZEF is the computation of the Vandermonde matrices with respect to the basis $\mathcal{T}^{k}_{P}$

[TABLE]

Here $VT$ is computed simply by the formula $T_{h}(x)=\cos(h\arccos(x))$ while for $WT$ we prefer to use the recursion algorithm to improve the stability of the computation since the points $P(\zeta_{i}),$ $\zeta_{i}\in\Omega,$ may in general lie outside of $[-1,1]^{2}$ or even be complex.

The second step of the algorithm is the most delicate: we perform the orthonormalization of $VT$ and store the triangular matrix defining the change of basis. More precisely the orthonormalization procedure is performed by applying the QR algorithm twice (following the so called twice is enough principle). First we apply the QR algorithm to $VT$ and we store the obtained $R_{1}$, then we apply the QR algorithm to $VT\setminus R_{1},$ we obtain $Q$ and $R_{2}$ that we store. Here $\setminus$ is the matlab backslash operator implementing the backward substitution; this is much more stable than the direct inversion of the matrix $R_{1}$, which may be ill conditioned. Note that $Q_{i,j}=p_{j}(z_{i}),$ where, since $Q$ is an orthogonal matrix,

[TABLE]

Therefore $\sqrt{M_{k}}Q_{i,j}=q_{j}(z_{i}).$

Step three. We compute the orthonormal polynomials evaluated at the points $\{\zeta_{1},\dots\zeta_{L}\}=\Omega.$ Again we prefer the backslash operator rather than the matrix inversion to cope with the possible ill-conditioning of $R_{1}$ and $R_{2}$: we compute

[TABLE]

Note that $W_{i,j}=p_{j}(\zeta_{i})$ and thus $q_{j}(\zeta_{i})=W_{i,j}\sqrt{M_{k}}.$ We also compute the matrix $K=Q\cdot W^{T}$, here we have

[TABLE]

Step 4. Finally we have

[TABLE]

and

[TABLE]

3.2.2. SZEF-BW Algorithm Implementation

First the step 1 and step 2 of SZEF algorithm are performed.

Step 3. The Bergman weight $\sigma=(B_{k}/N(z_{1}),\dots,B_{k}(z_{M})/N)$ is computed by

[TABLE]

Then the weighted Vandermonde matrix $V^{(w)}_{i,j}:=\sqrt{\sigma_{i}}Q_{i,j}$ is computed by a matrix product.

Step 4. We compute the orthonormal polynomials by another orthonormalization by the QR algorithm: $V^{(w)}=Q^{(w)}\cdot R^{(w)}$. We get

[TABLE]

We also compute the evaluation of $\tilde{q}_{j}$s at $\Omega$ by

[TABLE]

Step 5. We perform the step 4 of the SZEF algorithm, where $Q$ and $W$ are replaced by $\tilde{Q}$ and $\tilde{W}.$

3.3. Numerical Tests of SZEF and SZEF-BW

The extremal function $V_{E}^{*}$ can be computed analytically for very few instances as real convex bodies and sub-level sets of complex norms. For the unit real ball $B$ Lundin (see for instance [25]) proved that

[TABLE]

where $h(\zeta):=\zeta+\sqrt{1-\zeta^{2}}$ denotes the inverse of the Joukowski function and maps conformally the set $\mathbb{C}\setminus[-1,1]$ onto $\mathbb{C}\setminus\{z\in\mathbb{C}:|z|\leq 1\}.$

Let $E$ be any real convex set, one can define the convex dual set $E^{*}$ as

[TABLE]

Baran [1, 2] proved, that if $E$ is a convex compact set symmetric with respect to [math] and containing a neighbourhood of [math], the following formula holds.

[TABLE]

Here $extrF$ is the set of points of $F$ that are not interior points of any segment laying in $F.$

In order to be able to compare our approximate solution with the true extremal function, we build our test cases as particular instances of the Baran’s Formula:

• •

test case 1, $E_{4}:=[-1,1]^{2},$

• •

test case 2, $E_{m}$ the real $m$-agon centred at $0,$

• •

test case 3, $E_{\infty}:=B$, the real unit disk.

We measure the error with respect to the true solution $V_{E}$ in terms of

[TABLE]

which should be intended as a quadrature approximation for the $L^{1}$ norm of the error.

Also we consider, as an approximation of the relative $L^{1}$ error, the quantity

[TABLE]

and the following approximation of the error in the uniform norm

[TABLE]

In all the tests we performed the rate of convergence is experimentally sub-linear, i.e., $\|s_{k+1}\|/\|s_{k}\|\to 1$ as $k\to\infty$, where $s_{k}:=f_{k+1}-f_{k}$ for $f_{k}$ being one of $u_{k},\tilde{u}_{k},v_{k},\tilde{v}_{k}$ and $\|\cdot\|$ being one of the pseudo-norms we used above for defining $e_{1}$ and $e_{\infty}.$ This slow convergence may be overcame by extrapolation at infinity with the vector rho algorithm (see [20]). We present below experiments regarding both the original and the accelerated algorithm.

3.3.1. Test case 1

In this case equation (33) reads as

[TABLE]

To build an admissible mesh for $[-1,1]^{2}$ we can use the Cartesian product of an admissible mesh $X_{k}$ for one dimensional polynomials up to degree $k;$ see [16]. We can pick $X_{k}:=-\cos(\theta_{j}),$ $\theta_{j}:=j\pi/\lceil mk\rceil,$ with $j=0,1,\dots,\lceil mk\rceil$.

First we compare the performance of our four approximations in terms of $L^{1}$ error behaviour as $k$ grows large. Also, in order to better understand the rate of convergence, we compute the ratios

[TABLE]

for $f_{k}$ in $\{u_{k},v_{k}\tilde{u}_{k},\tilde{v}_{k}\}$ and $\Omega=[-2,2]^{2}.$ We report the results in Figure 1 and 2 respectively. The profile of convergence is slow but monotone, indeed the asymptotic constants $s_{k}$ become rather close to $1.$ This suggest a sub-linear convergence rate.

It is worth to say that in all the tests we made the point-wise error is much smaller far from $E$ than for points that lie near to $E$: in the experiment reported in Figure 1 and Figure 2 $\Omega$ is an equispaced real grid in $[0,2]^{2}$, but the results are the same for larger bounds and becomes much better if $\Omega\cap E=\emptyset$ or when $\Omega$ is a purely imaginary set, i.e., $\Omega\cap\mathbb{R}^{2}=\emptyset.$

Fortunately, even if the convergence rate shown by the experiment of Figure 1, the vector rho algorithm, see for instance [20], works effectively on our approximation sequences and actually allows to produce much better approximations than the original sequences $\{u_{k}\},\{\tilde{u}_{k}\},\{v_{k}\},\{\tilde{v}_{k}\}.$ We report in Figure 3 a comparison among the errors, on $[-2,2]^{2}$ and $[0,20]^{2}$ respectively, of the original sequence $\{u_{k}\}$ and the accelerated sequence produced by the vector rho algorithm, together with a linear (with respect to log-log scale) fitting of the original errors. In contrast with the sub-linear convergence enlightened above, the accelerated sequence exhibits (in the considered interval for $k$) a super-quadratic behaviour.

3.3.2. Test case 2

To construct a weakly admissible mesh on the real regular $m$-agon

[TABLE]

we use the algorithm proposed in [17]. First the convex polygon is subdivided in non overlapping triangles, then an admissible mesh $\{\hat{A}_{k,l}\}_{k}$ for each triangle $l,$ $l=1,2,\dots,m$ is created by the Duffy transformation; finally, the union $\{\tilde{A}_{k}\}:=\{\cup_{l=1}^{m}\hat{A}_{k,l}\}$ of all meshes (with no point repeated) is an admissible mesh for the polygon by definition.

The resulting mesh $\tilde{A}_{k}$ has good approximation properties: for instance it can be constructed in such a way to have a very small constant, however $\tilde{A}_{k}$ is not tailored to our problem. Points of $\tilde{A}_{k}$ cluster, as one could expect, at any corner of the regular polygon, but also they cluster near the edges of each triangle in the subdivision. This ”spurious” clustering is not coming from the geometry of the problem, instead it is an effect of the method we used to solve it. More importantly, this issue tends to deteriorate the convergence of SZEF algorithm. Let us briefly give an insight on why this phenomena occurs in the following remark.

Remark 3.4**.**

In Section 5 we will prove (see Proposition 5.1) that the sequence of measures

[TABLE]

where $\rightharpoonup^{*}$ denotes the convergence in the weak∗ topology on Borel measures. Assume for simplicity that $\mu_{E}$ is absolutely continuous with respect to the Lebesgue measure $m$ restricted to $E$ and consider an admissible mesh that tends to cluster points on a ball $B(z_{0},r)\subset E$ for which $\frac{\mu_{E}}{dm}(z),z\in B(z_{0},r),$ is very small. The asymptotic above implies in particular that $B_{k}^{\mu_{k}}(z)$ will be very small for $z\in\operatorname*{supp}\mu_{k}\cap B(z_{0},r).$

On the other hand we have $N_{k}^{-1}\int B_{k}^{\mu}d\mu=1$ for any measure and thus there will be some points $\hat{z}\in\operatorname*{supp}\mu_{k}\setminus B(z_{0},r)$ for which $B_{k}^{\mu_{k}}(\hat{z})$ is very large. Recall that the uniform convergence of Theorem 3.1 relies on the fact that

[TABLE]

thus we should aim to get an admissible mesh $\{A_{k}\}$ having small constant, whose cardinality is slowly increasing, and whose Bergman function is ”not too large”, e.g., is ”not too oscillating”. From the observation above a good heuristic to achieve such an aim is to mimic the density of the equilibrium measure.

In order to overcome this issue we can apply two strategies.

The first one is to get rid of these extra points. To this aim we can use the AFP algorithm [19] to extract a set of discrete Fekete points $A_{k}$ of order $2k$ from $\tilde{A}_{2k}.$

Let us motivate this choice. First $A_{k}$ has been shown (numerically) to be a weakly admissible mesh in many test cases, see [19] and reference therein. Also, even more importantly, the resulting Bergman function is (experimentally) much less oscillating than the one that can be computed starting by $\tilde{A}_{k}$ and this turns in a smaller uniform norm on $E_{m}$ of such a function. Lastly, heuristically speaking in view of Remark 3.4, we would like to use a mesh which is mimicking the distribution of the equilibrium measure. We invite the reader to compare this last discussion with [10] where the definition of optimal measures is introduced.

We can also perform a different choice which rests upon Theorem 3.2. Indeed the SZEF-BW algorithm uses a Bergman weighting of $\mu_{k}$ to prevent the phenomena we discussed in Remark 3.4. We tested our algorithm SZEF-BW for different values of $m$ and several choices of $\Omega,$ here we consider the case $m=6.$

Again, the convergence is quite slow but (numerically) monotone and, apart from the region of $\Omega_{0}$ close to $\partial_{\mathbb{R}^{n}}E,$ it is not affected by the particular choice of $\Omega;$ that is, we can appreciate numerically the global uniform convergence proved in Theorem 3.2. Moreover, extrapolation at infinity is successful even in this test case. We report profiles of convergence relative to two possible choices of $\Omega$ in Figure 4.

3.3.3. Test case 3

Lastly we consider $E=E_{\infty}:=\{(z_{1},z_{2}):\Im(z_{i})=0,\Re(z_{1})^{2}+\Re(z_{2})^{2}\leq 1\};$ in such a case the true solution is computed by the Lundin Formula (32). We can easily construct an admissible mesh of degree $k$ for the real unit disk following [16], for, it is sufficient to take a set of Chebyshev-Lobatto points $\{\eta_{0},\eta_{1},\dots,\eta_{s}\},$ $s>k$ and set

[TABLE]

We report the behaviour of the error and the convergence profile in Figure 5, again the sub-linear rate of convergence is rather evident.

Also, we compare the profile of the error function $e_{k}(z):=|\tilde{v}_{k}(z)-V_{E_{\infty}}^{*}(z)|$ on different two real dimensional squares in Figure 6. The global uniform convergence of $\tilde{v}_{k}$ to $V_{E_{\infty}}^{*}$ theoretically proven in Theorem 3.2 reflects on our experiments: $e_{k}$ is small (approximately $10^{-2}$ for $k\geq 30$) and very flat away from $E_{\infty}$ while it attains its maximum on $E_{\infty}$ with a fast oscillation near $\partial_{\mathbb{R}^{2}}E_{\infty}.$ Again, the extrapolation at infinity improves the quality of our approximation.

4. Approximating the transfinite diameter

In this section we present a method for approximating the transfinite diameter of a real compact set.

4.1. Theoretical result

Given any basis $\mathcal{B}_{k}=\{b_{1},\dots,b_{N_{k}}\}$ of the space $\mathscr{P}^{k}$ and a measure $\mu$ inducing a norm on $\mathscr{P}^{k}$, we denote by $G_{k}(\mu,\mathcal{B}_{k})$ its Gram matrix, namely we set

[TABLE]

The hermitian matrix $G_{k}(\mu,\mathcal{B}_{k})$ has square root, indeed introducing the generalized Vandermonde matrix

[TABLE]

we have $V_{k}(\mu,\mathcal{B}_{k})^{\textsc{H}}V_{k}(\mu,\mathcal{B}_{k})=G_{k}(\mu,\mathcal{B}_{k}).$ Note that, for $\mu$ being the uniform probability measure supported on an array of unisolvent points of degree $k$, $V_{k}(\mu,\mathcal{B}_{k})$ is the standard Vandermonde matrix for the basis $\mathcal{B}_{k}$ and divided by $\sqrt{N_{k}}.$

We recall, see [13], the following relation between Gram determinants and $L^{2}$ norms of Vandermonde determinants.

[TABLE]

Where $\mathcal{M}_{k}$ denotes the (graded lexicographically ordered) monomial basis.

Here is the main result this section.

Theorem 4.1**.**

Let $E\subset\mathbb{C}^{n}$ be a compact $\mathcal{L}$-regular set and $\{A_{k}\}$ a (weakly) admissible mesh for $E$ then, denoting by $\mu_{k}$ the uniform probability measure on $A_{k}$, we have

[TABLE]

Proof.

We use equation (34). Since $\mu_{k}$ is a probability measure, it follows that

[TABLE]

hence

[TABLE]

On the other hand, by the sampling property of admissible meshes it follows that, for any polynomial $p$ of degree at most $k$ we have $\|p\|_{E}\leq C\|p\|_{A_{k}}\leq C\sqrt{\operatorname*{Card}A_{k}}\|p\|_{L^{2}_{\mu_{k}}}.$ Thus, since $\det\operatorname{{Vdm}}_{k}(\zeta_{1},\dots,\zeta_{N_{k}})$ is a polynomial in each variable $\zeta_{i}\in\mathbb{C}^{n}$ of degree not larger than $k$, we get

[TABLE]

Since $(\operatorname*{Card}A_{k})^{1/k}\to 1$, being $\{A_{k}\}$ weakly admissible, it follows that

[TABLE]

∎

In principle Theorem 4.1 provides an approximation procedure for $\delta(E)$, given $\{A_{k}\}$, however the straightforward computation of the left hand side of (35) leads to stability issues. In the next subsection we present our implementation of an algorithm based on Theorem 4.1 and we discuss a possible way to overcome such difficulties.

4.2. Implementation of the TD-GD algorithm

We recall that we denote by $T_{j}(z)$ the classical $j$-th Chebyshev polynomial and we set $\mathcal{T}_{k}:=\{\phi_{1},\dots,\phi_{N_{k}}\},$ where

[TABLE]

where $\alpha:\mathbb{N}\to\mathbb{N}^{2}$ is the one defined in (11) and $\pi_{h}$ is the $h$-coordinate projection.

We denote by $V_{k}=V_{k}(A_{k},\mathcal{T}_{k})$ the Vandermonde matrix of degree $k$ with respect the mesh $A_{k}:=\{(x_{1},y_{1}),\dots,(x_{M_{k}},y_{M_{k}})\}$ and the basis $\mathcal{T}_{k}$, that is

[TABLE]

similarly we define $W_{k}:=V_{k}(A_{k},\mathcal{M}_{k})$ where the chosen reference basis is the lexicographically ordered monomial one.

Now we notice that, setting $M_{k}:=\operatorname*{Card}A_{k},$

[TABLE]

thus we have

[TABLE]

The direct application of this procedure leads to a unstable computation that actually does not converge.

On the other hand, the computation of the Gram determinant in the Chebyshev basis,

[TABLE]

is more stable and we have

[TABLE]

Here the matrix $P_{k}$ is the matrix of the change of basis. Again the numerical computation of $\det P_{k}$ becomes severely ill-conditioned as $k$ grows large.

Instead, our approach is based on noticing that $P_{k}$ does not depend on the particular choice of $E$, thus we can compute the term $\left(\det(P_{k})\right)^{\frac{n+1}{nkN_{k}}}$ once we know $(\det G_{k}(\hat{\mu}_{k}))^{\frac{n+1}{2nkN_{k}}}$ and $(\det\tilde{G}_{k}(\hat{\mu}_{k}))^{\frac{n+1}{2nkN_{k}}}$ for a particular $\hat{\mu}_{k}$ which is a Bernstein Markov measure for $\hat{E}\subseteq[-1,1]^{2}$ as

[TABLE]

Also we can introduce a further approximation, since $\det G_{k}(\hat{\mu}_{k})^{\frac{n+1}{2nkN_{k}}}\to\delta(\hat{E}),$ we replace in the above formula $\det G_{k}(\hat{\mu}_{k})^{\frac{n+1}{2nkN_{k}}}$ by $\delta(\hat{E}).$ Finally, we pick $\hat{E}:=[-1,1]^{2}$ and $\hat{\mu}_{k}$ uniform probability measure on an admissible mesh for the square, for instance the Chebyshev Lobatto grid with $(2k+1)^{2}$ points, thus our approximation formula becomes

[TABLE]

where we used $\delta([-1,1]^{2})=1/2;$ [9].

Finally, to compute the determinants of the Gram matrices on the right hand side of equation (38) we use the SVD factorization of the square root of the Gram matrices, note that for instance

[TABLE]

where $V_{k}$ and $M_{k}$ has been defined above and $S_{k}=diag(\sigma_{1},\dots,\sigma_{N_{k}})$ is the diagonal matrix with the singular values the matrix $V_{k}/\sqrt{M_{k}}.$

4.3. Numerical test of the TD-GD algorithm

In order to illustrate how our algorithm works in practice, we perform two numerical tests for real compact sets whose transfinite diameters have been computed analytically in [9]. Namely, we consider the case of the unit disk $B^{2}=\{x\in\mathbb{R}^{2}:|x|\leq 1\}$ and the unit simplex $S^{2}:=\{x\in\mathbb{R}_{+}^{2}:x_{1}+x_{2}\leq 1\}$ For such sets Bos and Levenberg computed formulas that in the specific case of dimension $n=2$ read as

[TABLE]

4.3.1. TD-GD test case 1: the unit ball in $\mathbb{R}^{2}$

To compute the approximation of $(G_{k}(\hat{\mu}_{k}))^{\frac{n+1}{2nkN_{k}}}$ we pick $\hat{\mu}_{k}$ to be the probability measure on the point set

[TABLE]

This is a well known admissible mesh of constant $2$ of the square $[-1,1]^{2}$ for the space of $k$ tensor degree polynomials that in particular include the space $\mathscr{P}^{k}$; [23, 16].

For $\mu_{k}$ we use an admissible mesh $A_{k}$ of degree $k$ built as in [16] using the radial symmetry of the unit disk. Here

[TABLE]

The admissible meshes $A_{k}$ and $\hat{A}_{k}$ are displayed in Figure 7.

We compute the right hand side of equation (38) for the sequence of values $k=4,6,\dots,28$ and we report both the absolute and the relative errors in Figure 8 (continuous line without and with diamonds respectively). On one hand we notice that the convergence rate is very slow, but on the other hand the sequence of approximations is monotone and the error structure is good for the application of the extrapolation. Indeed we report the absolute and relative error (dashed line without and with diamonds respectively) of the sequence obtained by the diagonal of the rho table (rho algorithm) in the same figure. Notice that the absolute error of the accelerated sequence at degree $28$ is $4.9105\cdot 10^{-6}$, that is, * six digits of $\delta(B^{2})$ are computed exactly.*

4.3.2. TD-GD test case 2: the unit simplex in $\mathbb{R}^{2}$

The first part of the algorithm for the computation of $\delta(S^{2})$, i.e., the computation of the factor in (38) coming from (37), is identical to the one we performed for $\delta(B^{2}).$

Then we pick an admissible mesh on $S^{2}$ following [16]. Our mesh $A_{k}$ at the $k$-th stage is the image under the Duffy transformation of a Chebyshev grid on $[-1,1]^{2}$ formed by $(4k+1)^{2}$ points, see Figure 9. We recall that the Duffy transformation (with a suitable choice of parameters) maps the unit square onto the simplex and any degree $k$ polynomial on the simplex is pulled back by the Duffy transformation onto the square to a polynomial of degree not larger than $2k.$ It follows that $\{A_{k}\}$ is an admissible mesh of constant $2$ for the simplex; [16].

Once the mesh has been defined the numerical computations to get the right hand side of (38) are performed as above. The obtained results, both in terms of absolute and relative errors, are displayed in Figure 10.

Again, the defined algorithm is very slowly converging, nevertheless using extrapolation at infinity by the rho algorithm we get a sequence rather fast converging. Indeed, more than six exact digits of $\delta(S^{2})$ can be computed in less than 10 seconds even on a rather outdated laptop, e.g., Intel CORE i3-3110M CPU, 4 Gb RAM.

5. Approximating the equilibrium measure

Fekete points are (at least theoretically) the first tool to investigate how the equilibrium measure looks like for a given regular compact set $E\subset\mathbb{C}^{n}$. Indeed, the main result of [6] asserts that the sequence of uniform probability measures supported at the $k$-th stage on a Fekete array of order $k$ is converging weak∗ to the equilibrium measure of the considered set. Unfortunately Fekete arrays are known analytically for very few instances and they are characterized in general as solutions of an extremely hard optimization problem, hence, though its strong theoretical motivation, this method is not of practical interest.

However, the results in [6] are in fact more general (as shown also in [10]): one can take asymptotically Fekete arrays (see equation (15)) and obtain the same result. This is indeed the approach of [16, Th. 1], where the asymptotically Fekete arrays are produced by a discretizing the optimization problem using an admissible mesh as optimization domain. A $50$-th stage of an asymptotic Fekete array for a regular hexagon is reported in Figure 11.

Another strategy to get a sequence of (weighted) point masses approaching the equilibrium measure in the weak∗ sense is based on the Bergman Asymptotic (21) and the use of admissible meshes. We summarize this in the following proposition, which is a consequence of the results of [6] and [10].

Proposition 5.1**.**

Let $E\subset\mathbb{C}^{n}$ a regular compact set. Let $\{A_{k}\}$ be a weakly admissible mesh for $E$ and let $\mu_{k}$ be the uniform probability measure supported on $A_{k}.$ Denote by $\tilde{\mu}_{k}$ the measure $N_{k}^{-1}B_{k}^{\mu_{k}}.$ Then we $\tilde{\mu}_{k}$ converges weak∗ to $\mu_{E}.$

Sketch of the proof.

First we note that the sequence of measures $\{\mu_{k}\}$ leads to the transfinite diameter, i.e. one has the asymptotic result (35) as proven above. This is the starting point for applying the machinery of [10].

Indeed for such a sequence of measures one still has Lemma 2.7 and Lemma 2.8 of [10], thus using the derivative of the Aubin-Mabuchi energy functional and Lemma 3.1 of [6], one gets

[TABLE]

∎

Approximations of the equilibrium measure built by means of point masses may have certain interest when one aims to perform approximated computations with equilibrium measure, for instance computing orthogonal series, since the approximation is given in terms of a quadrature rule. On the other hand, such a kind of approximation can not be easily represented to get an insight on how the equilibrium measure looks like for a given $E$; this property becomes relevant if one aims to test or argue conjectures.

In the rest of the section we will introduce an approximation scheme for $\mu_{E}$ based on absolutely continuous measures with respect to the standard Lebesgue measure.

Our method is based on the following lemma.

Lemma 5.1**.**

Let $\mu$ be a positive Borel measure. Let us set, for any $k\in\mathbb{N},$

[TABLE]

then we have

[TABLE]

Here $\operatorname*{adj}$ denotes the adjugate of a matrix.

Proof.

First notice that $B_{k}^{\mu}(z)=\sum_{h=1}^{N_{k}}q_{h}(z,\mu)\overline{q_{h}}(z,\mu)$ is a smooth function never vanishing in $\mathbb{C}^{n}$, hence we can use classical differentiation with no problems.

We have

[TABLE]

Also, using the linearity of differentiation and the tensor structure of $B_{k}^{\mu}=b^{H}b$ we get

[TABLE]

So we can write

[TABLE]

Lastly we use the Matrix Determinant Lemma, i.e., $\det(A+uv^{T})=\det A+v^{T}\operatorname*{adj}(A)u,$ and the fact that $\det(\lambda A)=\lambda^{n}\det A$ to get equation (39). ∎

We already shown, see Theorem 3.1, that, for the sequence $\{\mu_{k}\}$ of uniform probability measures supported on a weakly admissible mesh for $E$, one has the asymptotic

[TABLE]

locally uniformly. We recall also that the Monge Ampere operator is continuous under the local uniform limit (see for instance [25]), thus

[TABLE]

where the limit is to be intended in the sense of the weak∗ topology of Borel measures. Therefore we have the following.

Theorem 5.1**.**

Let $E\subset\mathbb{C}^{n}$ a regular compact set. Let $\{A_{k}\}$ be a weakly admissible mesh for $E$ and denote by $\mu_{k}$ the uniform probability measure supported on $A_{k}.$ Let us denote by $\eta_{k}$ the sequence of functions

[TABLE]

The sequence $(2i)^{n}\eta_{k}d\operatorname*{Vol}_{\mathbb{C}^{n}}$ converges weak∗ to $\mu_{E}.$ In particular, when $D$ has full rank, we have

[TABLE]

where $S=diag(\sigma_{1},\dots,\sigma_{n})$ is the diagonal matrix $R^{H}R$ and $D=QR$ is the standard QR factorization of $D.$

Proof.

The only thing that remains to prove is equation (40). It is sufficient to simply notice that if $A$ is any invertible matrix, then $\operatorname*{adj}(A)=(\det(A))A^{-1}.$ In our specific case, in which $S=R^{H}R$ for a triangular matrix $R$, $\det S=\prod_{l=1}^{n}\sigma_{l}$ factors out and (40) follows ∎

Remark 5.2**.**

Note that the measures $\eta_{k}$ are not a priori supported on $E$, however it follows trivially by the above theorem that also the sequence of measures having density $\eta_{k}\chi_{E}$ (i.e., the restriction to $E$ of $\eta_{k}$) has the same weak∗ limit.

Acknowledgements

The findings of this work are essentially a part of the doctoral dissertation [32]. Consequently, much of what we present here have been deeply influenced by the discussions with the Advisor Prof. N. Levenberg (Indiana University).

All the software used in the numerical tests we performed above has been developed in collaboration with Prof. M. Vianello (University of Padova). The author deeply thanks him both for the scientific collaboration and the support.