The Monge-Ampere equation: various forms and numerical methods

TL;DR

This paper introduces three new forms of the Monge-Ampere equation, emphasizing a Fourier integral form with proven positivity and bounds, leading to a novel numerical method for solving space-periodic problems with applications in cosmology.

Contribution

The paper develops a new Fourier integral form of the Monge-Ampere equation and a corresponding numerical method with proven properties, advancing computational approaches in this area.

Findings

Established positivity and bounds for the Fourier integral kernels.

Developed a convergent numerical method for space-periodic Monge-Ampere problems.

Demonstrated method effectiveness on a cosmological test problem.

Abstract

We present three novel forms of the Monge-Ampere equation, which is used, e.g., in image processing and in reconstruction of mass transportation in the primordial Universe. The central role in this paper is played by our Fourier integral form, for which we establish positivity and sharp bound properties of the kernels. This is the basis for the development of a new method for solving numerically the space-periodic Monge-Ampere problem in an odd-dimensional space. Convergence is illustrated for a test problem of cosmological type, in which a Gaussian distribution of matter is assumed in each localised object, and the right-hand side of the Monge-Ampere equation is a sum of such distributions.

Tables

Table 1. Table 1: Number of evaluations (NE) of the determinant of the Hessian, performed by AICDM in the course of computation of approximate solutions ( 7 ) of varying accuracy to the test MAE with the r.h.s. defined by ( 22 ), ( 23 ) and ( 36 ), discrepancy d ∞ ​ ( η ) subscript 𝑑 𝜂 d_{\infty}(\eta) for these solutions, and wallclock duration of runs (DR, seconds) for computation of approximate solutions to the test MAE by EAICDM.

$d\left(\eta \right)$	${10}^{-3}$	${10}^{-4}$	${10}^{-5}$	${10}^{-6}$	${10}^{-7}$	${10}^{-8}$	${10}^{-9}$	${10}^{-10}$
$d_{\infty }\left(\eta \right)$	$0.64\times {10}^{-2}$	$0.86\times {10}^{-3}$	$1.02\times {10}^{-3}$	$0.23\times {10}^{-4}$	$0.28\times {10}^{-5}$	$0.39\times {10}^{-6}$	$0.36\times {10}^{-7}$	$0.41\times {10}^{-8}$
NE	64	126	251	439	1 018	1 515	2 030	2 653
DR	3.38	4.69	11.17	42.7	81	203	435	774

Equations

$$\det \Vert u_{x_i{x}_j}\Vert =f(\mathbf{x})$$

$$u={\nabla }^{-2}\sqrt{u_{x_1{x}_1}^2+u_{x_2{x}_2}^2+2u_{x_1{x}_2}^2+2f}$$

$$u=\int_{R^{N}}\widetilde{u}(\mbox{\boldmath$\omega$\unboldmath})e^{i\mbox{\boldmath$\omega$\unboldmath}\cdot\bf x}\,d\mbox{\boldmath$\omega$\unboldmath}$$

$$\det \Vert a_{ij}\Vert =\frac{1}{N!}\sum _{i_1,...,i_N,j_1,...,j_N}{\epsilon }_{i_1...i_N}{\epsilon }_{j_1...j_N}\prod _{n=1}^N{a}_{i_n{j}_n}$$

$$\det\|u_{x_{i}x_{j}}\|={(-1)^{N}\over N!}\,\sum_{i_{1},...,i_{N},j_{1},...,j_{N}}\varepsilon_{i_{1}...i_{N}}\varepsilon_{j_{1}...j_{N}}\prod_{n=1}^{N}\int_{R^{N}}\widetilde{u}(\mbox{\boldmath$\omega$\unboldmath})\omega_{i_{n}}\omega_{j_{n}}e^{i\mbox{\boldmath$\omega$\unboldmath}\cdot\bf x}\,d\mbox{\boldmath$\omega$\unboldmath}$$

$$={(-1)^{N}\over N!}\,\int_{R^{N}}...\int_{R^{N}}\left(\,\sum_{i_{1},...,i_{N}}\varepsilon_{i_{1}...i_{N}}\prod_{n=1}^{N}\omega^{n}_{i_{n}}\right)\left(\,\sum_{j_{1},...,j_{N}}\varepsilon_{j_{1}...j_{N}}\prod_{n=1}^{N}\omega^{n}_{j_{n}}\right)\left(\,\prod_{n=1}^{N}\widetilde{u}(\mbox{\boldmath$\omega$\unboldmath}^{n})\right)$$

$$\times\exp\left(i\sum_{n=1}^{N}\mbox{\boldmath$\omega$\unboldmath}^{n}\cdot{\bf x}\right)\,d\mbox{\boldmath$\omega$\unboldmath}^{1}...\,d\mbox{\boldmath$\omega$\unboldmath}^{N}$$

$$={(-1)^{N}\over N!}\,\int_{R^{N}}...\int_{R^{N}}{\det}^{2}\left\|\mbox{\boldmath$\omega$\unboldmath}^{1},...,\mbox{\boldmath$\omega$\unboldmath}^{N-1},\mbox{\boldmath$\omega$\unboldmath}-\sum_{n=1}^{N-1}\mbox{\boldmath$\omega$\unboldmath}^{n}\right\|$$

$$\times\left(\,\prod_{n=1}^{N-1}\widetilde{u}(\mbox{\boldmath$\omega$\unboldmath}^{n})\right)\widetilde{u}\!\left(\mbox{\boldmath$\omega$\unboldmath}-\sum_{n=1}^{N-1}\mbox{\boldmath$\omega$\unboldmath}^{n}\right)e^{i\mbox{\boldmath$\omega$\unboldmath}\cdot\bf x}\,d\mbox{\boldmath$\omega$\unboldmath}^{1}...\,d\mbox{\boldmath$\omega$\unboldmath}^{N-1}d\mbox{\boldmath$\omega$\unboldmath}$$

$$\!={(-1)^{N}\over N!}\!\int_{R^{N}}\!...\!\int_{R^{N}}\!{\det}^{2}\!\left\|\mbox{\boldmath$\omega$\unboldmath}^{1}\!\!,...,\mbox{\boldmath$\omega$\unboldmath}^{N-1}\!\!,\mbox{\boldmath$\omega$\unboldmath}\right\|\left(\,\prod_{n=1}^{N-1}\widetilde{u}(\mbox{\boldmath$\omega$\unboldmath}^{n})\!\right)\widetilde{u}\!\left(\!\mbox{\boldmath$\omega$\unboldmath}-\!\!\sum_{n=1}^{N-1}\mbox{\boldmath$\omega$\unboldmath}^{n}\!\right)\!e^{i\mbox{\boldmath$\omega$\unboldmath}\cdot\bf x}d\mbox{\boldmath$\omega$\unboldmath}^{1}...\,d\mbox{\boldmath$\omega$\unboldmath}^{N-1}d\mbox{\boldmath$\omega$\unboldmath}.$$

$$\frac{1}{N!}\sum _{i_1,...,i_N,j_1,...,j_N}{\epsilon }_{i_1...i_N}{\epsilon }_{j_1...j_N}{\left(u_{x_{i_1}{x}_{j_1}}...u_{x_{i_{N-1}}{x}_{j_{N-1}}}u\right)}_{x_{i_N}{x}_{j_N}}=f.$$

$$\frac{1}{N!}\sum _{i_1,...,i_N,j_1,...,j_N}{\epsilon }_{i_1...i_N}{\epsilon }_{j_1...j_N}{\int }_{R^N}{u}_{x_{i_1}{x}_{j_1}}...u_{x_{i_{N-1}}{x}_{j_{N-1}}}u{\phi }_{x_{i_N}{x}_{j_N}}d\mathbf{x}={\int }_{R^3}f\phi d\mathbf{x}.$$

$$⟨f⟩\equiv \underset{R\to \infty }{\lim }\frac{1}{|B_R|}{\int }_{B_R}f(\mathbf{x})d\mathbf{x},$$

$${\int }_{T^N}fd\mathbf{x}=0,$$

$$u=c\left(\frac{|\mathbf{x}{|}^2}{2}+u^{\prime }\right),$$

$$c^N=⟨f⟩.$$

$$\nabla^{2}u^{\prime}=\int_{R^{N}}\widetilde{\eta}(\mbox{\boldmath$\omega$\unboldmath})e^{i\mbox{\boldmath$\omega$\unboldmath}\cdot\bf x}\,d\mbox{\boldmath$\omega$\unboldmath},\quad u^{\prime}=\int_{R^{N}}\widetilde{u}^{\prime}(\mbox{\boldmath$\omega$\unboldmath})e^{i\mbox{\boldmath$\omega$\unboldmath}\cdot\bf x}\,d\mbox{\boldmath$\omega$\unboldmath},\quad\hbox{where \ }\widetilde{u}^{\prime}(\mbox{\boldmath$\omega$\unboldmath})=-\widetilde{\eta}(\mbox{\boldmath$\omega$\unboldmath})/|\mbox{\boldmath$\omega$\unboldmath}|^{2},$$

$$f/c^{N}=\int_{R^{N}}\widetilde{f}(\mbox{\boldmath$\omega$\unboldmath})e^{i\mbox{\boldmath$\omega$\unboldmath}\cdot\bf x}\,d\mbox{\boldmath$\omega$\unboldmath}.$$

$$\det\|u^{\prime}_{x_{i}x_{j}}\|={1\over N!}\,\int_{R^{N}}...\int_{R^{N}}{\det}^{2}\left\|{\bf i}_{\mbox{\boldmath$\omega$\unboldmath}^{1}},...,{\bf i}_{\mbox{\boldmath$\omega$\unboldmath}^{N-1}},{\bf i}_{\mbox{\boldmath$\omega$\unboldmath}-\sum_{n=1}^{N-1}\mbox{\boldmath$\omega$\unboldmath}^{n}}\right\|$$

$$\times\left(\,\prod_{n=1}^{N-1}\widetilde{\eta}(\mbox{\boldmath$\omega$\unboldmath}^{n})\right)\widetilde{\eta}\!\left(\mbox{\boldmath$\omega$\unboldmath}-\sum_{n=1}^{N-1}\mbox{\boldmath$\omega$\unboldmath}^{n}\right)e^{i\mbox{\boldmath$\omega$\unboldmath}\cdot\bf x}\,d\mbox{\boldmath$\omega$\unboldmath}^{1}...\,d\mbox{\boldmath$\omega$\unboldmath}^{N-1}d\mbox{\boldmath$\omega$\unboldmath},$$

$${(-1)^{m}\over N!}\sum_{i_{1},...,i_{N},j_{1},...,j_{N}}\!\!\varepsilon_{i_{1}...i_{N}}\varepsilon_{j_{1}...j_{N}}\sum_{|\sigma|=m}\!\!\left(\,\prod_{n:~{}i_{n},j_{n}\in\sigma}\int_{R^{N}}\widetilde{u}^{\prime}(\mbox{\boldmath$\omega$\unboldmath})\omega_{i_{n}}\omega_{j_{n}}e^{i\mbox{\boldmath$\omega$\unboldmath}\cdot\bf x}\,d\mbox{\boldmath$\omega$\unboldmath}\right)\prod_{n:~{}i_{n}{~{}\rm or~{}}j_{n}\in\!\!\!/\sigma}\!\!\delta_{i_{n}j_{n}}$$

$$={(-1)^{m}\over m!}\,\sum_{1\leq p_{1}<...<p_{N-m}\leq N}\,\int_{R^{N}}...\int_{R^{N}}\!\left(\,\sum_{j_{1},...,j_{m}}\varepsilon_{j_{1}...j_{m}p_{1}...p_{N-m}}\prod_{n=1}^{m}\mbox{\boldmath$\omega$\unboldmath}^{n}_{j_{n}}\right)^{\!\!2}$$

$$\times\left(\,\prod_{n=1}^{m}\widetilde{u}^{\prime}(\mbox{\boldmath$\omega$\unboldmath}^{n})\right)\exp\left(i\sum_{n=1}^{m}\mbox{\boldmath$\omega$\unboldmath}^{n}\cdot{\bf x}\right)\,d\mbox{\boldmath$\omega$\unboldmath}^{1}...\,d\mbox{\boldmath$\omega$\unboldmath}^{m}$$

$$=\int_{R^{N}}...\int_{R^{N}}A_{m}\left({\bf i}_{\mbox{\boldmath$\omega$\unboldmath}^{1}},...,{\bf i}_{\mbox{\boldmath$\omega$\unboldmath}^{m-1}},{\bf i}_{\mbox{\boldmath$\omega$\unboldmath}-\sum_{n=1}^{m-1}\mbox{\boldmath$\omega$\unboldmath}^{n}}\right)$$

$$\times\left(\,\prod_{n=1}^{m-1}\widetilde{\eta}(\mbox{\boldmath$\omega$\unboldmath}^{n})\right)\widetilde{\eta}\!\left(\mbox{\boldmath$\omega$\unboldmath}-\sum_{n=1}^{m-1}\mbox{\boldmath$\omega$\unboldmath}^{n}\right)e^{i\mbox{\boldmath$\omega$\unboldmath}\cdot{\bf x}}\,d\mbox{\boldmath$\omega$\unboldmath}^{1}...\,d\mbox{\boldmath$\omega$\unboldmath}^{m-1}d\mbox{\boldmath$\omega$\unboldmath},$$

$$A_m({\mathbf{i}}^1,...,{\mathbf{i}}^m)\equiv \frac{1}{m!}\sum _{1\le p_1<...<p_{N-m}\le N}{M}_{p_1...p_{N-m}}^2({\mathbf{i}}^1,...,{\mathbf{i}}^m)$$

$$M_{p_1...p_{N-m}}({\mathbf{i}}^1,...,{\mathbf{i}}^m)\equiv \sum _{j_1,...,j_m}{\epsilon }_{j_1...j_m{p}_1...p_{N-m}}({\mathbf{i}}^1{)}_{j_1}...({\mathbf{i}}^m{)}_{j_m},$$

$${\mathcal{M}}_m\equiv \Vert {\mathbf{i}}^1,...,{\mathbf{i}}^m\Vert ,$$

$$\widetilde{\eta}(\mbox{\boldmath$\omega$\unboldmath})+\sum_{m=2}^{N}\int_{R^{N}}...\int_{R^{N}}A_{m}\left({\bf i}_{\mbox{\boldmath$\omega$\unboldmath}^{1}},...,{\bf i}_{\mbox{\boldmath$\omega$\unboldmath}^{m-1}},{\bf i}_{\mbox{\boldmath$\omega$\unboldmath}-\sum_{n=1}^{m-1}\mbox{\boldmath$\omega$\unboldmath}^{n}}\right)$$

$$\times\left(\,\prod_{n=1}^{m-1}\widetilde{\eta}(\mbox{\boldmath$\omega$\unboldmath}^{n})\right)\widetilde{\eta}\!\left(\mbox{\boldmath$\omega$\unboldmath}-\sum_{n=1}^{m-1}\mbox{\boldmath$\omega$\unboldmath}^{n}\right)d\mbox{\boldmath$\omega$\unboldmath}^{1}...\,d\mbox{\boldmath$\omega$\unboldmath}^{m-1}=\widetilde{f}(\mbox{\boldmath$\omega$\unboldmath}),$$

Full text

1 Introduction

The Monge–Ampère equation (MAE)

[TABLE]

is encountered in many areas of numerical analysis and physics, ranging from image processing [1, 2, 3] to cosmology [4, 5, 6]. Here the subscripts $x_{i}$ denote derivatives in the respective spatial variables; $\|u_{x_{i},x_{j}}\|$ is the Hessian of $u$, i.e., the matrix comprised of second derivatives of $u$. Existence and regularity of its solutions was considered in [7, 8, 9]. Various strategies were proposed for its numerical solution. A provably convergent method for solving the Dirichlet problem for the two-dimensional MAE in bounded convex domains was presented in [10]. It is directly linked to the geometric interpretation of the equation that had given an opportunity to demonstrate existence of weak solutions [13]. However, the actual numerical examples in [10] have just up to 25 grid points, which is clearly insufficient to establish the practical feasibility of the method.111More recently, a similar algorithm was used [11] to design a beam–shaping reflector; computations on a $55\times 55$ grid required about 15 minutes of a Pentium 4 processor. This problem of geometric optics can be recast as the Monge–Kantorovich mass transfer problem, reducing to a PDE of the MAE type, see [12]. The MAE can be recast as a minimisation problem; application of algorithms for saddle-point optimisation to a two-dimensional MAE was considered in [14, 15], and least-squares minimisation was advocated in [16]. Efficient methods for solution of the discrete optimal transportation problem were developed for application to cosmological problems [4, 5, 6]. It may be perceived that discrete methods correspond well to the physics of mass transportation in the Universe — after all, galaxies, as observed by astronomers, are clearly well-localised, discrete objects! This argument, however, becomes less persuading, when one recalls that visible matter constitutes only a very small fraction of different, “dark” kinds of matter, whose density distribution is supposed to be continuous (as opposed to a nearly discrete distribution in clusters of the visible matter).

In another group of methods the MAE is treated as a generic nonlinear partial differential equation. Application of Galerkin’s and finite element methods to a regularised MAE was considered in [17, 18]. A pseudospectral Newton’s algorithm was reported to perform well for a MAE in $R^{2}$ with a smooth r.h.s. [19]. However, when we implemented this algorithm for a three-dimensional MAE, we found that it failed to converge, unless the r.h.s. was a smooth slowly varying function. This is consistent with the observation [10] that Newton’s method, applied to equations equivalent to the MAE, is useful to improve an approximate solution only if the approximation is accurate enough, and easily fails otherwise. Two methods for numerical solution of a Dirichlet problem for an elliptic MAE in a two-dimensional convex region $\Omega$ are examined in [20]. The first one employs a finite-difference discretisation of the equation; it is advocated for application to the MAE with a possibly singular solution. The second one is an iterative method for the MAE in the form of a fixed point problem

[TABLE]

(here $\nabla^{-2}$ denotes the inverse Laplacian); it is claimed to perform better, when the solution is regular (i.e., belongs to the Sobolev space $W^{2}_{2}(\Omega)$ ).

A so-called “inexact” iterative Newton–Krylov solver with preconditioning was applied for two- [21] and three-dimensional [22] grid generation with the properties of equidistribution and minimum distortion. The nature of the problem solved in [21, 22] required only a modest accuracy of solutions, discrepancies the order of $10^{-3}-10^{-4}$ being acceptable. Examples of computations of such moderately accurate solutions illustrating the scalability of the algorithm with respect to the employed spatial resolution were provided ibid.

In Sections 2-4 we derive three alternative forms of (1): the “second-order divergence” form, the Fourier integral form and the “convolution” form, which, to the best of our knowledge, were never presented in the literature before. The second form suggests two related methods for computation of space-periodic solutions to (1). The methods and results of their test applications are presented in Section 6; the Monge-Ampère problem of the cosmological type, which we use to test our algorithms, is presented in Section 5, after a short statement of the cosmological reconstruction problem. In the Concluding remarks we briefly discuss some open questions.

We note that the MAE’s of various kinds are considered in the literature. For the MAE arising in the mass transportation problem, typically the r.h.s. of (1) depends on the gradient of the unknown function (this form of the equation can be found in the original memoir by Ampère [24]). In the theory of PDE’s, usually one considers the MAE with a known r.h.s. (see [9, 23]). The MAE of this kind arises in the reconstruction of the early Universe (see discussion at the beginning of Section 5). Since this cosmological problem is the main original motivation for our work, we restrict our discussion to the case of the “mathematician’s” MAE. However, a straightforward reformulation of our algorithms makes them applicable for the case of the r.h.s. depending on the unknown function.

2 The “second-order divergence” form of the MAE

The divergence form of the MAE (1), in which the l.h.s. of the equation is represented as a sum of the first derivatives of certain quantities, is well known (see, e.g., [17]). In this section we derive a representation of the l.h.s. of the MAE as a sum of second derivatives, which we call the second-order divergence form of the MAE.

Consider a Fourier integral solution

[TABLE]

to the MAE in $R^{N}$. Substituting the integral into (1) and using the identity for $N\times N$ matrices

[TABLE]

where each of the indices $i_{1},...,i_{N},j_{1},...,j_{N}$ takes the values $1,...,N$ and $\varepsilon_{p_{1}...p_{N}}$ denotes the unit antisymmetric tensor of rank $N$, we find

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here \left\|\mbox{\boldmath\omega\unboldmath}^{1},...,\mbox{\boldmath\omega\unboldmath}^{N}\right\| denotes a $N\times N$ matrix comprised of $N$ columnar vectors \mbox{\boldmath\omega\unboldmath}^{1},...,\mbox{\boldmath\omega\unboldmath}^{N}. The first factor in the integrand of (3) is a quadratic function of $\omega$, hinting that the l.h.s. of (1) can be transformed into a sum of second derivatives. Indeed, “reverse engineering” of (3) reveals an alternative, “second-order divergence” form of (1) in $R^{N}$:

[TABLE]

Equivalence of (1) and (4) is now easily established directly: Using the standard rule for differentiation of the products in the l.h.s. of (4), we render it as a sum of products of derivatives of $u$. Third-order derivatives enter such a product only in pairs of the form $u_{x_{i_{p}}x_{i_{N}}x_{j_{p}}}u_{x_{i_{q}}x_{j_{N}}x_{i_{q}}}$, and hence all products involving third-order derivatives cancel out due to antisymmetry of the tensors $\varepsilon_{i_{1}...\,i_{N}}$. Similarly, all terms involving fourth-order derivatives cancel out. Consequently, each term in the l.h.s. of (4) gives rise to a single term $\displaystyle\varepsilon_{i_{1}...\,i_{N}}\varepsilon_{j_{1}...\,j_{N}}u_{x_{i_{1}}x_{j_{1}}}...\,u_{x_{i_{N}}x_{j_{N}}}$, and hence the sum is the determinant of the Hessian of $u$, as required.

This form has an interesting consequence. Suppose, for the sake of simplicity, that a space-periodic solution $u\in T^{N}$ is sought. Let $\varphi$ be a smooth function with a finite support (in this context such $\varphi$ are called test functions). Multiplying the MAE by $\varphi$ and twice integrating by parts each term in the sum, we obtain by virtue of (4)

[TABLE]

As usual in the theory of partial differential equations, a weak solution to (1), $u$, can be defined as a function satisfying the integral identity (5) for any test function $\varphi$. (Other definitions of weak solutions to the MAE are natural: generalised solutions obtained by geometric constructions [13] and the so-called viscosity solutions [9]; the two definitions are equivalent [23].) Commonly (e.g., see [17, 25]), only one integration by parts is performed in this integral identity. Our form is advantageous in that a lesser regularity of the weak solution is required for (5) to be well-defined. The following argument illustrates this for $N>2$: For any $u$ from the Sobolev space $W^{2}_{N-1}(T^{N})$ the integrals in the l.h.s. of (5) are well-defined (because by the Sobolev embedding theorem this implies $\nabla u\in L_{2(N-1)}(T^{N})$ and hence $u\in L_{\infty}(T^{N})$ ). By contrast, integrals in the similar identity obtained by just a single integration by parts are not well-defined for $u\in W^{2}_{N-1}(T^{N})$. (Note that a viscosity solution for fully nonlinear second-order equations is only required to be continuous [9].)

3 The Fourier integral form of the MAE; positivity and bounds of kernels

3.1 Derivation

In the Monge–Ampère–Kantorovich approach to the cosmological reconstruction problems [4, 5], the MAE arises for the potential of the inverse Lagrangian map, in which the function $f$ is a ratio between matter density at the current epoch and at the much earlier epoch of matter-radiation decoupling (at the earlier epoch, the distribution is very close to uniform); thus, $f>0$.

Here, we consider a problem in which the spatial mean of $f$ is positive. Note that for odd $N$, if $\langle f\rangle<0$, the change of the unknown function $u\to-u$ reverses the sign of $f$ and $\langle f\rangle$ admits the desirable sign (here $\langle\cdot\rangle$ denotes space averaging:

[TABLE]

and $|B_{R}|$ is the volume of the ball $B_{R}\subset R^{N}$ of radius $R$). Suppose that $u$ and $f$ have the same spatial periodicity; integration of (4) over a periodicity cell $T^{N}$ then yields

[TABLE]

which is incompatible with $\langle f\rangle\neq 0$. Thus we assume henceforth

[TABLE]

where $u^{\prime}$ has the periodicity of $f$ and a zero spatial mean (this being just a normalisation), and $c$ is a constant. Substitution of (7) into (1) and integration over a periodicity cell yields

[TABLE]

To derive this, note that $u_{x_{i}x_{j}}=c(\delta_{ij}+u^{\prime}_{x_{i}x_{j}})$, where $\delta_{ij}$ is the Kronecker symbol, and hence the derivatives of $u^{\prime}$ in the l.h.s. of (1) are present only in $\det\|u^{\prime}_{x_{i}x_{j}}\|$ and a linear combination of minors of smaller sizes of the Hessian of $u^{\prime}$. By the algebraic transformation presented in the previous section each such minor can be converted into a sum of second derivatives of products of $u^{\prime}$ with its second derivatives. Thus the spatial mean of the l.h.s. of (1) over any periodicity cell, and hence the mean defined by (6), does not involve $u^{\prime}$ and is equal to $c^{N}$.

In a more general formulation, we seek a solution (7) satisfying (8) and $\langle u^{\prime}\rangle=0$, assuming that $f$ and $u^{\prime}$ can be represented as Fourier integrals:

[TABLE]

[TABLE]

(If the problem, defined by (1), (7) and (8), is considered for space-periodic $f$ and $u^{\prime}$, integrals in wave vectors in what follows are replaced by the respective Fourier sums.)

Eq. (2) is equivalent to

[TABLE]

[TABLE]

where $\bf i_{a}$ denotes a unit vector in the direction of $\bf a$. Our immediate goal is to derive a similar expression for the terms in (1) of lower orders in $u^{\prime}$. The term of order $m$ is

[TABLE]

(here the sum $\sum_{|\sigma|=m}$ is over all subsets $\sigma\subset\{1,...,N\}$ of cardinality $m$)

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

is the sum of squares of all minors of rank $m$,

[TABLE]

obtained by crossing out rows of numbers $p_{1}<...<p_{N-m}$ from the $N\times m$ matrix

[TABLE]

comprised of $m$ columnar vectors ${\bf i}^{1},...,{\bf i}^{m}$.

Therefore, the problem (1) is reduced after the substitution (7) to the system of integral equations, which we call the Fourier integral form of the problem (1) and (7):

[TABLE]

[TABLE]

which is now stated in terms of the Fourier coefficients \widetilde{\eta}(\mbox{\boldmath\omega\unboldmath}) of the Laplacian of the unknown function $u^{\prime}$. Equations (11) are valid for all \mbox{\boldmath\omega\unboldmath}\neq 0; the respective equation for \mbox{\boldmath\omega\unboldmath}=0 is (8).

If the r.h.s. of (1) is zero-mean ($\langle f\rangle=0$), the MAE admits the Fourier integral form (11) (where \widetilde{\eta}(\mbox{\boldmath\omega\unboldmath}) are now the Fourier coefficients of $\nabla^{2}u$), where the l.h.s. is reduced to a single term for $m=N$ in the sum $\sum_{m=2}^{N}$.

3.2 Bounds for the kernels $A_{m}$ in the Fourier integral form

In this subsection we establish bounds

[TABLE]

provided all vectors ${\bf i}^{s}$ have a unit norm. These bounds will play a crucial rôle for our numerical algorithm.

Addition to a column ${\bf i}^{s}$ of any linear combination of columns ${\bf i}^{s^{\prime}}$ for $s^{\prime}<s$ does not change the value of any minor $M_{p_{1}...p_{N-m}}({\bf i}^{1},...,{\bf i}^{m})$. Using the Gram–Schmidt orthogonalisation process, we change all ${\bf i}^{s}$ in ${\cal M}_{m}$ to ${\bf j}^{s}$ such that $(i)$ ${\bf j}^{1}={\bf i}^{1}$, $(ii)$ for any $s$, ${\bf j}^{s}$ differs from ${\bf i}^{s}$ by a linear combination of vectors ${\bf i}^{s^{\prime}}$ for $s^{\prime}<s$ and thus $A_{m}$ remains unaltered, $(iii)$ for any $s$, ${\bf j}^{s}$ is orthogonal to all ${\bf j}^{s^{\prime}}$ for $s^{\prime}<s$. Hence

[TABLE]

where $\theta_{s}$ is the angle between ${\bf i}^{s}$ and the subspace spanned by $\{{\bf i}^{s^{\prime}}|~{}s^{\prime}<s\}$. Consequently,

[TABLE]

We denote ${\cal M}^{\prime}_{m}\equiv\left\|{\bf j}^{1},...,{\bf j}^{m}\right\|$ and ${}^{t}{\cal M}^{\prime}_{m}$ the transpose of ${\cal M}^{\prime}_{m}$. The identity [26]

[TABLE]

(which does not require orthogonality of $\{{\bf j}^{1},...,{\bf j}^{m}\}$) can be easily proved using the formula

[TABLE]

We enlarge the set $\{{\bf j}^{1},...,{\bf j}^{m}\}$ by vectors ${\bf j}^{s}$ for $s>m$ to a complete orthonormal basis in $R^{N}$. Let $\cal U$ be an orthogonal matrix comprised of the $N$ columnar vectors ${\bf j}^{1},...,{\bf j}^{N}$, and $\cal E$ be the $N\times m$ matrix, whose all entries are 0 except for ${\cal E}_{ss}=1$ for all $1\leq s\leq m$. Then ${\cal M}^{\prime}_{m}=\cal UE$ and hence

[TABLE]

where ${\cal I}_{m}$ is the identity matrix of size $m$. Consequently,

[TABLE]

This demonstrates that (12) are sharp bounds.

3.3 Solution of the MAE for a weakly fluctuating r.h.s.

If the fluctuating part of the r.h.s. in (1), $f-\langle f\rangle$, is small relative the mean $\langle f\rangle$, (11) suggests that \widetilde{\eta}(\mbox{\boldmath\omega\unboldmath})\approx\widetilde{f}(\mbox{\boldmath\omega\unboldmath}) and the nonlinear terms in (11) are small. Then the system can be solved by iteration:

[TABLE]

[TABLE]

Theorem

$1^{\circ}$. Suppose positive constants $C_{0}$ and $C_{1}$ satisfy the inequality

[TABLE]

[TABLE]

and $\widetilde{\eta}_{0}$ satisfies the inequality

[TABLE]

for $K=0$. Then iterates (14) are globally bounded: (17) holds true for all $K>0$.

$2^{\circ}$. Under the same conditions, the following inequalities are satisfied for all $K>0$:

[TABLE]

[TABLE]

where it is denoted

[TABLE]

Furthermore, if

[TABLE]

then iterations (14) converge to a solution to (11), which is unique in the ball (17).

The proof is elementary: Inequalities (17)–(19) stem from the bounds (12). By the contraction mapping principle, inequalities (18)–(20) imply convergence of iterations (14) in the norms of $L_{1}(R^{N})$ and $C(R^{N})$ in the space of Fourier coefficients, yielding a solution satisfying (17).

Evidently, equations (15) and (20), where equality is assumed instead of the inequalities, have positive solutions for any $N$. (For instance, $C_{0}=\sqrt{3}-4/3,\ C_{1}=1/3$ for $N=3$.) If for the chosen values of $C_{0}$ and $C_{1}$ (15) holds true, $C_{2}=1$ and the inequality (16) is strict, then all conditions of the Theorem become satisfied for a slightly smaller value of $C_{0}$.

Note that in view of (20) $C_{0}<1$ and hence (16) implies $f>0$ everywhere. Consequently, our Theorem is mostly of interest as a statement about convergence of iterations (14). Existence of weak solutions was proved by geometric methods in [13] for the MAE with an arbitrary positive r.h.s. in a compact convex domain. Although our Theorem is significantly weaker because of the strong restrictions on the r.h.s., we have proved it for a non-compact domain — the entire space.

4 The “convolution” form of the MAE

Following the same algebraic ideas, the MAE can be partially “integrated”. We again use the Fourier transform of $u^{\prime}$:

[TABLE]

Let

[TABLE]

where

[TABLE]

[TABLE]

That these conditions can be enforced, is elementary for \mbox{\boldmath\omega\unboldmath}\neq 0; the case \mbox{\boldmath\omega\unboldmath}=0 is not problematic, because $\widetilde{\xi}(0)=0$. Then

[TABLE]

is a real-valued function. (For a given $u^{\prime}$ it is not uniquely defined.)

Let us render the term of order $m$ in $u^{\prime}$ in the l.h.s. of (1) in the terms of $\xi({\bf x})$ employing the expression (9), (10), (13) and the “identity”

[TABLE]

(as usual understood in the sense of generalised functions):

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In particular, the linear term ($m=1$) is $\nabla^{2}u^{\prime}=\int_{R^{N}}\nabla\xi({\bf x}^{1})\cdot\nabla\xi({\bf x-x}^{1})\,d{\bf x}^{1}$.

Therefore, the problem (1) is equivalent, after the substitution (7) and (8), to

[TABLE]

We call this integral equation the “convolution” form of the MAE. Since it does not involve second-order derivatives, it may prove useful for development of an iterative algorithm for numerical solution of the problem defined by (1), (7) and (8), which involves suitable transformations of the spatial variable.

If the r.h.s. of (1) is zero-mean, the MAE also admits the convolution form (11), where $c=1$, the l.h.s. is reduced to a single term for $m=N$ in the sum $\sum_{m=1}^{N}$, and \widetilde{u}^{\prime}(\mbox{\boldmath\omega\unboldmath}) in the definition of \widetilde{\xi}(\mbox{\boldmath\omega\unboldmath}) denotes the Fourier transform of $u$.

5 A test problem with a cosmological flavour

One important area of application of the MAE is the reconstruction of the dynamical history of the Universe from present observations of the spatial distribution of masses (galaxies, clusters, including their dark-matter components). Let us briefly recall the background. Peebles was the first to propose that, from the sole knowledge of the current positions of galaxies (from the Local Group which includes our own galaxy) without knowledge of their (proper) velocities, reconstruction of the full dynamical history is a meaningful goal [27]. Indeed, the very strong constraint on the distribution of masses at the epoch of decoupling — which has to be almost uniform — makes reconstruction a possibly well-posed two-point boundary problem that can then be solved by variational techniques. In fact, using convexity techniques, a theorem of uniqueness of reconstruction was proved in [5], using the Euler–Poisson equations, which describes the dynamics of matter on sufficiently large scales (of the order of a few million light years). In practice, reconstruction on such scales was done so far via the Monge–Ampère–Kantorovich (MAK) method [4] (see also [5, 6]) which assumes that the Lagrangian map from initial to current mass locations is the gradient of a convex potential. This holds exactly for the Zel’dovich approximation [28] (in which the equation for the velocity of matter reduces to the inviscid Burgers equation with a zero r.h.s., and hence the peculiar velocities remain constant along trajectories of point masses) and also for its refinement, the first-order Lagrangian perturbation approximation [29]. It is then easy, using mass conservation, to derive a MAE for the potential of the inverse Lagrangian map (the latter is the Legendre–Fenchel transform of the potential of the direct map). By a theorem of Brenier [30] the MAE becomes a Monge–Kantorovich mass transportation problem with quadratic cost which, after discretisation, can be solved by optimisation techniques (see [5] for details). The r.h.s. of the MAE is equal to the ratio of densities at the present and initial positions. Thus, if the initial density was not uniformly distributed, we would have to solve the MAE with the r.h.s. depending on the gradient of the unknown potential. However, at the epoch of decoupling (about 380 thousand years after the Big Bang) the density of matter was close to uniform, which implies a significant simplification of the MAE to be solved: its r.h.s. is known.

Here we explore for the first time the possibility of directly solving the three-dimensional MAE without discretisation on a toy model with a cosmological flavour.

We assume that the dimension of space is $N=3$, and the r.h.s. of the MAE has the following structure:

[TABLE]

The function (22) describes mass distribution for $G$ “objects”, which are “localised”, if the value of the parameter $\delta$ is small compared to the distance between objects. In the cosmological context galaxies or clusters of galaxies can be regarded as such objects; then $f^{(g)}$ describes the total distribution of mass of both visible matter and dark matter, in which the galaxies are embedded. Although observations attest, that the distribution of visible matter in galaxies is clearly discontinuous, it is astrophysically sound to expect that the distribution of all types of matter is smooth due to the prevailing smoothness of the dark matter. We assume that objects have density distributions with a Gaussian shape:

[TABLE]

The $g$-th object of mass $m^{(g)}>0$ is located at ${\bf r}^{(g)}$. All $m^{(g)}$, $\sigma^{(g)}$ and ${\bf r}^{(g)}$ are independent of the small parameter $\delta>0$.

We will be seeking a solution (7) with a space-periodic $u^{\prime}$ to the Monge–Ampère problem (1), where a space-periodic r.h.s. $\widehat{f}$ is the sum of “clones” of (22) over all periodicity cells. Without loss of generality we can assume that the total mass is normalised:

[TABLE]

thus (8) implies $c=1$ in (7).

It is instructive to consider particular solutions to this Monge–Ampère problem.

5.1 An exact solution to the MAE for a spherically symmetric mass distribution

In spherical coordinates centred at ${\bf r}^{(g)}$, (1) becomes, for spherically symmetric $u$ and $f$,

[TABLE]

where $\rho=|{\bf r-r}^{(g)}|$. This equation has an obvious solution

[TABLE]

whose first derivatives are uniformly bounded:

[TABLE]

the second derivatives are $O(\delta^{-1})$:

[TABLE]

This estimate follows from the following inequalities:

$\bullet$ $|x_{i}x_{m}/\rho^{2}|\leq 1$;

$\bullet$ for $\rho<\delta$, $f(\rho/\delta)$ is uniformly bounded and

[TABLE]

for some positive constant $\underline{c}$ and $\overline{c}$;

$\bullet$ for $\rho\geq\delta$, $(\rho/\delta)^{2}f(\rho/\delta)$ is uniformly bounded and

[TABLE]

for some constant $\underline{c}^{\prime}$ and $\overline{c}^{\prime}$.

Moreover, if $\rho$ is larger than a positive constant, $f(\rho/\delta)=o(\delta^{3})$, and hence all second derivatives of $u$ are uniformly bounded outside any sphere of a fixed radius.

The properties of the solution discussed in this subsection are implied just by a fast decay of $f$ at infinity and the spherical symmetry of the r.h.s.; the assumption that the profile is Gaussian is unnecessary.

5.2 A one-cell solution for $G>1$ spherically symmetric localised objects

A solution for $G>1$ objects can be expected to admit a power series expansion

[TABLE]

where (by analogy with (24)) each $u_{n}({\bf r},\delta)$ and its first derivatives are $O(\delta^{0})$, and the second derivatives are $O(\delta^{-1})$.

Naively one could expect interaction of localised objects to be asymptotically unimportant, and hence to obtain in the leading order the sum of individual one-object solutions:

[TABLE]

Here $u^{(g)}(|{\bf r-r}^{(g)}|,\delta)$ is the one-object solution (24) for the $g$-th object. Let us inspect, to what extent such a conjecture might be true.

Consider a neighbourhood of an object $\gamma$. In the leading order, the l.h.s. of (1) is $\det\|{u_{0}\,}_{x_{i}x_{j}}\|$. Substitution of (25) into the determinant yields a sum of triple products of second derivatives of various $u^{(g)}$. As shown in the previous subsection, for $g\neq\gamma$ each second derivative of $u^{(g)}$ is $O(\delta^{0})$ in this neighbourhood (the distance between the objects $g$ and $\gamma$ is a fixed positive constant). Thus, only triple products of second derivatives of $u^{(\gamma)}$ contribute terms of the leading order, $O(\delta^{-3})$, the one of the l.h.s. These products constitute $\det\|u^{(\gamma)}_{x_{i}x_{j}}\|$ and thus match the term $\delta^{-3}f^{(\gamma)}$ in the r.h.s. Nonlinear interaction of pairs of solutions for individual objects is at a lower order, $O(\delta^{-2})$, and that of triplets of one-object solutions - still weaker, $O(\delta^{-1})$. Hence, after the finite sum (24) is substituted into the MAE, the highest order terms do cancel out, and in this respect the naive conjecture is confirmed.

Nevertheless, the one-cell solution (25) does not represent the leading order term of the solution to the problem at hand. To see this, note that at large distances (25) exhibits a linear growth in $|\bf r|$, and not a quadratic one, as the form of the solution (7) requires. Furthermore, after subtraction of the quadratic profile $c|{\bf r}|^{2}$, the remaining part $u^{\prime}$ must be space-periodic, which is clearly not the case for (25). We cannot use (25) to construct a global solution conformant with this periodicity requirement by the procedure, used to obtain the periodic r.h.s. $\widehat{f}$ from the individual density profiles (22), because the sum of periodically distributed “clones” of one-cell solutions (25) is infinite. Moreover, for such hypothetical space-periodic sum there would be infinitely many pairwise interactions between objects yielding products of the order $O(\delta^{-2})$; hence, even if the sum of such products is finite, it will not necessarily remain $O(\delta^{-2})$.

We observe that the pairwise interaction does not become weaker when the distance between the interacting objects grows. Thus, even in the high-contrast limit the solution that we are seeking does not reduce to individual interactions between objects; their interaction is essentially collective, which makes it hard to predict the structure of the solution.

6 Computation of a solution to the MAE with a space-periodic r.h.s.

In this section we present two iterative algorithms for numerical solution of the space-periodic Monge–Ampère problem. One version (AICDM), employing numerical improvement of convexity and discrepancy minimisation stabilising the iterative process, is suitable for computation of solutions for everywhere positive right-hand sides $f$ (see Subsection 6.4). Another one (ACPDM), involving continuation in parameter and discrepancy minimisation, does not require this condition to be satisfied (see Subsection 6.2). They both rely on the basic algorithm for iterative solution of a fixed point problem for the MAE (see Subsection 6.1). Test applications of the two algorithms are considered in Subsections 6.3 and 6.4 . We assume $\langle f\rangle\neq 0$, however, reformulation of ACPDM for the case of a zero-mean space-periodic r.h.s. is straightforward.

6.1 A basic solver (FPAR)

If the amplitude of fluctuation of $f$ is small compared to the mean, — more precisely, if (16) is satisfied, — the Theorem (see Subsection 3.3) establishes existence of a solution to (1), that is a perturbation of $c|{\bf x}|^{2}/2$. If, in addition, (20) holds true, the Theorem guarantees that iterations defined by (14) converge to the Laplacian of the perturbation; this offers an algorithm for numerical solution of the MAE. If either of the conditions (16) or (20) are violated (which is the practically interesting case), iterations (14) do not necessarily converge. Different algorithms are necessary, taking this into account.

Suppose the kernels $A_{m}$ in (11) are frozen and take constant values $a_{m}$, respectively. Then the new system of equations (11) can be easily solved, since it is the Fourier form of the polynomial equation

[TABLE]

in the physical space. This observation suggests the following algorithm. We express (11) in the terms of $\eta=\nabla^{2}u^{\prime}$ as

[TABLE]

where

[TABLE]

and implement iterations

[TABLE]

in the physical space.

The applicability of this algorithm depends crucially on the choice of the coefficients $a_{m}$. In view of the positivity and the bounds for the kernels (12), we impose

[TABLE]

Furthermore, it seems practical to set

[TABLE]

because for $a_{1}=1$ the linear in $u^{\prime}$ term is treated exactly,

[TABLE]

and the median values (28) of $A_{m}$ for $m>1$ minimise the ranges of the kernels

[TABLE]

in (29). Since for the coefficients (28)

[TABLE]

the l.h.s. of (27) “captures” the nonlinear behaviour of the l.h.s. of (11), the algorithm has chances to converge. However, for the extreme values of $A_{m}$ the values of the kernels in (29) are as large as the respective medians, and therefore convergence of iterations (27), (28) is not guaranteed.

At least two strategies can be proposed for setting the value of $a_{0}$ (which is a free parameter in the sense that (11) is not required to be satisfied for \mbox{\boldmath\omega\unboldmath}=0):

$1^{\circ}$. At each iteration $a_{0}$ is tuned, so that $\langle\eta\rangle=0$ (which holds true for $\eta=\nabla^{2}u^{\prime}$).

$2^{\circ}$. We set $a_{0}=0$.

Note that although in the subspace of space-periodic functions the inverse Laplacian is defined for zero-mean scalar fields only, $\langle\eta_{K}\rangle$ is not required to vanish, because the mean is removed by differentiation in the Hessian before the inverse Laplacian is evaluated.

For any odd $N$, the equation (27) in $\eta_{K}$ with the coefficients (28) has a unique root for any r.h.s. To check this, it is enough to establish that the derivative $D_{N}(\eta)$ of the l.h.s. of (27) is positive for any $\eta$. For the choice of coefficients (28),

[TABLE]

At a minimum $D_{N}^{\prime}(\eta)=0$ and hence (30) implies that at the minimum $D_{N}>0$, proving monotonicity of the l.h.s. of (27). It can be shown similarly, that for any even $N$ the number of roots is 0 or 2. Thus, the algorithm is guaranteed to be applicable for odd $N$ only.

It can be proved, that the Theorem also applies for iterations (27) (with the same constant $C_{1}$ bounding the same norm of solutions).

In an application to a test problem inspired by cosmology (see Subsection 6.3), this algorithm, which we call FPAR (Fixed Point Algorithm for the Regular part of the MAE) produces a sequence of iterations, initially converging, but subsequently blowing up: linear instability sets in, and hence the respective unstable mode must be removed.

6.2 A more advanced solver (ACPDM)

The behaviour of the basic algorithm FPAR suggests that it should be embedded as an engine within a more advanced algorithm. In this subsection we present such an advanced solver, ACPDM (Algorithm with Continuation in a Parameter and Discrepancy Minimisation).

Consider a generalisation of (26):

[TABLE]

where

[TABLE]

Here the coefficients (28) are assumed, and the new parameter $p$ is confined to the interval $[0,1]$. For $p=0$, (31) is just a set of polynomial equations of degree $N$; for $p=1$ it reduces to (26) which is equivalent to (1). Continuation in the parameter $p$ is implemented: (31) is solved for a set of values $p_{j}$, increasing from 0 to 1, and an initial approximation of the solution for $p=p_{j}$ is obtained by polynomial extrapolation of solutions for all $p_{j^{\prime}}<p_{j}$. (Numerical extrapolation requires performing quadruple precision computations, if the number of nodes exceeds roughly a dozen.) For any $p$, solution of (31) involves iterations similar to (27):

[TABLE]

The r.m.s. discrepancy

[TABLE]

is used in the termination condition. Note that $\langle Q(\eta;1)\rangle=0$ and hence

[TABLE]

for $p=1$; in particular, $d(\eta)=0$ for $p=1$, if and only if for $u^{\prime}=\nabla^{-2}\eta$ and the normalisation (8) the field (7) is a solution to the MAE.

Let $(\cdot,\cdot)$ denote a scalar product and $\|\cdot\|$ the induced norm of a scalar field: $\|{\bf v}\|\equiv\sqrt{({\bf v,v})}$. In the test runs reported in the next subsection, the scalar product of the functional Lebesgue space $L^{2}(T^{3})$ has been assumed:

[TABLE]

Other products (with different weight functions introduced in the above integral) have been also considered (see Subsection 6.4).

If $\eta^{\prime}$ and $\eta^{\prime\prime}$ are approximate solutions to the generalised MAE (31), then

[TABLE]

where $A$ is the linearisation of (31) around the solution. Consequently, the concept of minimisation of the residual in Krylov spaces and approaches for its realisation can be borrowed from solvers for linear problems (such as the Generalised Conjugate Gradients Method [31]). ACPDM involves sequences of stabilised iterations described below, which exploit this concept. To make such an iteration, the following data is required: an approximate solution $\eta_{K}$, and two sets of $S$ scalar fields, $v_{s}({\bf x})$ and $w_{s}({\bf x})$, $0\leq s\leq S$, where $0\leq S\leq S_{\max}$ (as well as some other quantities computed at the previous stabilised iteration), and $S_{\max}$ is a parameter of the algorithm. It is supposed that all $v_{s}$ are mutually orthogonal with respect to the scalar product $(\cdot,\cdot)$, and

[TABLE]

for small $w_{s}$.

A sequence of stabilised iterations of ACPDM is initialised using the current approximation $\eta_{0}$ by computing $\eta_{1}=\eta^{\prime}_{1}$ as a solution to (32) for $K=1$, and setting $S=0$ (i.e. the sets $w_{s}({\bf x})$ and $v_{s}({\bf x})$ are empty). For $K>1$, a stabilised iteration of ACPDM consists of the following steps:

$i$. At each grid point in the physical space solve equation (32):

[TABLE]

$ii$. Compute $Q(\eta^{\prime}_{K})$.

$iii$. Orthogonalise $v^{\prime}\equiv Q(\eta^{\prime}_{K})-Q(\eta^{\prime}_{K-1})$ to all $v_{s}$, $1\leq s\leq S$, with respect to the scalar product $(\cdot,\cdot)$, and set

[TABLE]

$iv$. Compute

[TABLE]

$v$. If $S<S_{\max}$, increase $S$ by 1; otherwise (i.e., if $S=S_{\max}$), discard $v_{1}$ and $w_{1}$, and decrease by 1 the indices $s$ of the remaining $S_{\max}$ fields $v_{s}$ and $w_{s}$.

$vi$. Compute the r.h.s. of (32) for $\eta_{K+1}$ substituted in place of $\eta_{K-1}$, the field $Q(\eta_{K+1})$ and $d(\eta_{K+1})$. If $d(\eta_{K+1})$ is less than a given small threshold, then $\eta_{K+1}$ is the desired approximate solution and computation for the present $p$ is finished. If $\|Q(\eta_{K+1})\|>\|Q(\eta_{K})\|$, then the current sequence is terminated.

A few comments are in order. Clearly, $v_{S+1}$ and $w_{S+1}$ obtained in step $iii$, possess the required properties: $v_{S+1}$ is $(\cdot,\cdot)-$orthogonal to all $v_{s}$ for $s\leq S$, and (35) holds for $v_{S+1}$ and $w_{S+1}$ by virtue of (34). The coefficients in the sum computed in step $iv$ minimise the discrepancy $\|Q(\eta^{\prime}_{K}-\sum_{s=1}^{S}q_{s}w_{s})\|$, assuming all $q_{s}w_{s}$ are small and hence their quadratic contributions are negligible. The minimisation plays a dual rôle: on the one hand, it stabilises basic iterations (32), removing the instability modes as soon as they become substantial; on the other, it significantly increases the efficiency of the algorithm (up to a factor 20 compared to other algorithms relying on iterations (32)). The assumption that nonlinear terms are small can be incorrect; also, as a result of accumulation of the neglected nonlinear terms after a number of steps, at some stage $v_{s}$ can cease to approximate $Aw_{s}$ accurately enough. If the inequality $\|Q(\eta_{K+1})\|>\|Q(\eta_{K})\|$ is found to hold true in step $vi$, we interpret this as an indication that the adverse effect of nonlinearity has become significant, and then the algorithm breaks the current sequence. In order to reduce the adverse influence of nonlinearity, a small number $S_{\max}$ may be chosen; in our runs, $S_{\max}=5$.

Now we can assemble the algorithm from the building blocks, discussed above. ACPDM performs continuation in the parameter $p$. For a given $p$, it starts by carrying out basic iterations (32). As soon as convergence slows down (we have used the condition $d^{2}(\eta_{K})>d^{2}(\eta_{K-1})/2$ ), the algorithm switches to perform a sequence of stabilised iterations. If the sequence terminates in step $vi$ because nonlinear effects became significant, the algorithm proceeds by performing basic iterations (32). Usually it takes a small number of them for the instability to set in, and as soon as ACPDM detects that the inequality $d^{2}(\eta_{K})>\alpha d^{2}(\eta_{K-1})$ holds true, it starts a new sequence. Here $\alpha$ is a parameter, which can slightly exceed 1 (in order to allow the transients to die off and the dominant instability modes to set in, so that the latter could be efficiently removed in subsequent stabilised iterations); in our test runs we have used $\alpha=1.05$ .

6.3 Application to a test problem inspired by cosmology

We have chosen to employ a problem of the kind discussed in Section 5 as a test-bed for our algorithm, because cosmological applications of the MAE are probably the most important ones. The positiveness of the r.h.s. (implying convexity of solutions [13]) is not required for application of the algorithm.

We are seeking a solution (7) with a space-periodic $u^{\prime}$ to the Monge–Ampère problem (1), where a space-periodic r.h.s. $\widehat{f}$ is a sum of “clones” of (22) over all periodicity cells. Our poor man’s Universe involves $G=3$ objects in the periodicity cell $T^{3}([0,1]^{3})$, centred at the origin: $-1/2\leq x_{i}\leq 1/2$. They are described by the following parameter values in (23):

[TABLE]

[TABLE]

[TABLE]

Noting that $\widehat{f}({\bf x})$ achieves its maximum at ${\bf r}^{(3)}$ and its minimum at $\displaystyle{1\over 4}(1,1,1)$, we can estimate the contrast number of the problem as

[TABLE]

(the factor 4 in the denominator accounts for the replicas of (23), located in the neighbour periodicity boxes at a distance $1/\sqrt{2}$).

Iterations defined by (32) have been performed on a uniform grid comprised of $64^{3}$ points in the periodicity cell $T^{3}([0,1]^{3})$. We have evaluated $\det\|u^{\prime}_{x_{i}x_{j}}+\delta_{ij}\|$ by the pseudospectral method with dealiasing (for $N=3$ this requires computation of the derivatives $u^{\prime}_{x_{i}x_{j}}$ in the physical space on the twice finer grid involving $128^{3}$ points). Apparently dealiasing does not play any important rôle in these computations; this is consistent with the fast decay of the energy spectrum of $\eta=\nabla^{2}u^{\prime}$ by 11 orders of magnitude.

We have tested both strategies for choosing $a_{0}$ (see Subsection 6.1), as well as a hybrid strategy, in which $a_{0}=0$, but after solving (27) we reset $\eta_{K}:=\eta_{K}-\langle\eta_{K}\rangle$. It turns out that the hybrid strategy and $1^{\circ}$ are very close in the number of iterations necessary to obtain a solution to the MAE to the same accuracy $d(\eta_{K})<10^{-10}$. However, this implies that the hybrid strategy is several times more efficient in the terms of CPU time, since it requires just one evaluation of solutions to (27), while evaluation of $a_{0}$ following strategy $1^{\circ}$ requires several such evaluations. The strategy $2^{\circ}$ has proved to be the most efficient one, both in the terms of CPU time and the number of required iterations. (The same nodes $p_{j}$ were used in all the test runs.) In the remaining part of the subsection we discuss convergence of the advanced method with $a_{0}=0$.

ACPDM performs a combination of basic iterations (32) and stabilised iterations, which have different computational costs. The most time-consuming operation, determining the cost, is computation of the determinant of the Hessian in $F(\eta)$. There are two such operations in a stabilised iteration, and a single one in a basic iteration. Consequently, the former is approximately twice longer than the latter. To enable an accurate comparison of performance of various versions of our code, we report the number of evaluations of determinant in each run (a run is terminated when the obtained approximate solution $\eta$ satisfies the accuracy requirement $d(\eta)=10^{-10}$). Note however, that the main bulk of computations is performed by doing stabilised iterations; the number of basic iterations is usually below 1% of the number of stabilised ones.

In all runs, the initial iteration is a solution to (31) for $p=0$. If ACPDM is applied to this field for $p=1$, iterations quickly become chaotic and cease to converge (this has necessitated to include into the algorithm continuation in the parameter $p$). Initially, ACPDM has been applied for uniformly distributed nodes $p_{j}=j/J$, $j=0,...,J-1$ for $J=20$, for which the algorithm has shown a remarkably fast convergence (see Fig. 1). A polynomial 20-node extrapolation yields an approximate solution to our test problem to the accuracy $d(\eta)=0.47\times 10^{-2}$ and $d_{\infty}(\eta)=0.02$, where

[TABLE]

When it is used as an initial approximation for a run for $p=1$, convergence is slow and the pattern of convergence is erratic. It takes 38 685 evaluations of the determinant of the Hessian for the two discrepancy norms to decrease to $0.99\times 10^{-9}$ and $0.39\times 10^{-7}$, respectively. At this stage convergence stalls. We have therefore got into a local minimum of $d(\eta)$, out of which no exit can be found; we will refer to it as a spurious minimum solution.

A better approximation to the solution (7) to our test problem is obtained by adding more nodes near the right endpoint $p=1$. We have chosen to add $J^{\prime}=13$ nodes, the complete $p$-mesh being

[TABLE]

(information on convergence at the new nodes $p_{J+j}^{\phantom{2}}$ is also included in Fig. 1). Initial approximations at each new $p_{J+j}^{\phantom{2}}$ become more and more accurate in the interval $0<p\leq 0.55$ ($j=11$), then the discrepancy $d(\eta)$ of the initial approximations starts growing and admits the maximum $0.52\times 10^{-3}$ for $p=0.9875$ ($j=21$), and subsequently decreases again. The rate of convergence progressively falls down as $p_{J+j}^{\phantom{2}}$ approaches 1: the number of evaluations of the determinant of the Hessian yielding solutions of the desired accuracy $d(\eta(p_{J+j}^{\phantom{2}}))<10^{-10}$, markedly increases. However, ACPDM does not stall any more. A polynomial extrapolation involving the 33 nodes delivers an approximate solution $u^{\prime}$ for $p=1$ to the accuracy $d(\eta)=0.78\times 10^{-7}$ and $d_{\infty}(\eta)=0.20\times 10^{-5}$, and after further 1 885 evaluations of the determinant ACPDM yields an approximate solution with the two discrepancy norms down to $10^{-10}$ and $0.26\times 10^{-8}$, respectively. In total, 9 814 evaluations of the determinant of the Hessian are involved in computations on the mesh (37) providing the solution to the MAE.

The geometry of the obtained solution — isosurfaces of $u^{\prime}$ and $\nabla^{2}u^{\prime}$ — is shown in Figs. 2 and 3, respectively. The structures disclosed by the isosurfaces of $\nabla^{2}u^{\prime}$ (Fig. 3) are clearly associated with the three objects incorporated into the r.h.s. (22) of the test problem. The figures also reveal a subtle interaction of the objects along the lines connecting their centres (see in Fig. 3 the tube-like structures, which connect the regions of higher values of $\nabla^{2}u^{\prime}$ encompassing the centres of objects).

Figures 4 and 5 display the structure of the field $\nabla u^{\prime}$, determining the displacement of mass in the test Universe. It turns out that $u^{\prime}$ possesses hidden symmetries, namely, mirror reflection symmetries about any plane that is parallel to a coordinate plane and contains a pair of objects. (Deliberately, our computations are not sped up by exploiting these symmetries.) Because of the symmetries, the component of $\nabla u^{\prime}$ normal to such a plane vanishes on this plane. To illustrate the behaviour of the components of the gradient not shown in Fig. 4, we present in Fig. 5 their isolines on coordinate planes, which pass through the centre of the periodicity cube $T^{3}([0,1]^{3})$ and hence are displaced from the planes shown in Fig. 4 by a quarter of period.

6.4 A solver for MAE with an everywhere positive r.h.s. (AICDM)

Clearly, it is desirable to avoid the refinement of solutions for nodes $0.95<p_{j}<1$, which has proved necessary in application of ACPDM to our test problem. Comparison of the spurious minimum solution, obtained with the 20 nodes $p_{j}$, and the more accurate one is instructive. The maximum discrepancy between the two approximations of $\eta=\nabla^{2}u^{\prime}$ for these solutions, equal to 0.132, is attained at the minimum of the r.h.s., $(1/4,1/4,1/4)$, and all points, where the discrepancy is larger than 0.02 are located within the distance 1/16 from this point. The Hessian of (7) computed for the spurious minimum solution is negative at this point, i.e. this solution (7) is not a convex function, as it has to be [13]. This has suggested to develop a modification of ACPDM, presented here, which we call AICDM (Algorithm with Improvement of Convexity and Discrepancy Minimisation).

The algorithm proceeds as an extension of the ACPDM operating with a single node $p=1$; computation of a good initial approximation by extrapolation in $p$ is unnecessary. AICDM involves an additional procedure: improvement of convexity of an approximate solution, which is performed once the condition $d^{2}(\eta_{K})<\beta d^{2}(\eta_{K^{\prime}})$ is satisfied. Here $K$ is the number of the current iteration, $K^{\prime}$ is the number of the iteration, at which the last previous improvement of convexity was performed, and $\beta<1$ is a constant factor (we have chosen $\beta=0.01$ in the test computation reported in this subsection).

Improvement of convexity in AICDM could be performed by computation of a convex hull of the iterate, but we prefer to apply a numerically simpler procedure. At each grid point in the physical space, we compute the eigenvalues $\lambda$ of the Hessian of $u^{\prime}=\nabla^{-2}(\eta_{K}-\langle\eta_{K}\rangle)$ (they are all real, since the Hessian is a symmetric matrix). If they are all larger than $-1$, then the approximate solution given by (7) for the current iterate $\eta_{K}$ is locally convex at the point under consideration (in this discussion we assume that, after normalisation, $c=1$). Accordingly, if the minimum eigenvalue $\lambda_{\min}$ exceeds $-1$, no action is taken. Suppose now the contrary, i.e. $\lambda_{\min}<-1$. The minimum eigenvalue would become $-1$, if at this point each second derivative $u^{\prime}_{x_{i}x_{i}}$ is increased by $-1-\lambda_{\min}$, and hence the Laplacian $\eta_{K}=\nabla^{2}u^{\prime}$ is increased by $3(-1-\lambda_{\min})$. Following this observation, at the points where $\lambda_{\min}<-1$ we increase $\eta_{K}$ by $6(-1-\lambda_{\min})$ (by choosing the factor 6 instead of 3 we are “overimproving” $\eta_{K}$). It is not guaranteed, of course, that $|{\bf r}|^{2}/2+\nabla^{-2}(\eta^{\prime}_{K}-\langle\eta^{\prime}_{K}\rangle)$ is convex for the resultant $\eta^{\prime}_{K}$, since this procedure changes all mixed derivatives at each point and does not guarantee that each second derivative $u^{\prime}_{x_{i}x_{i}}$ at the points of local non-convexity is increased by the same amount, or that non-convexity does not appear at new grid points. Hence we proceed, repeating the procedure till each eigenvalue at each grid point becomes larger than $-1-d(\eta_{K})/2$. In our experience, only a small number of such iterations is necessary to enforce this condition. The quantities $6(-1-\lambda_{\min})$ are small, being at most comparable with the global discrepancy $d(\eta_{K})$; although the discrepancy $d(\eta^{\prime}_{K})$ for the “improved” approximation $\eta^{\prime}_{K}$ can exceed the discrepancy $d(\eta_{K})$ for the original approximation $\eta_{K}$, it thus turns out that the increase (compared to $d(\eta_{K})$) is usually modest.

An approximate solution to the test problem (formulated in Subsection 6.3), satisfying $d(\eta)=10^{-10}$ and $d_{\infty}(\eta)=0.41\times 10^{-8}$, has been obtained by the AICDM. This has required 2 653 evaluations of the determinant of the Hessian. Table 1 illustrates the deterioration of convergence in this run, as better accuracy numerical solutions are successively found. Note that computational cost of each iteration in improvement of convexity of an approximate solution slightly exceeds that of computation of the determinant of the Hessian, and each iteration is included into evaluation counts (one evaluation per a convexity improving iteration) presented in Tables 1 and 2.

In this run, as in runs reported in Subsection 6.3, the Lebesgue space scalar product (33) has been assumed. We have also inspected convergence of AICDM employing scalar products

[TABLE]

with weight functions

[TABLE]

(see Table 2). For $q>0$, more prominence is given to discrepancy in the regions, where the r.h.s. of the MAE, $\widehat{f}$, admits relatively high values; for $q<0$, the opposite happens, i.e. discrepancy in the regions, where the values of $\widehat{f}$ are relatively low, is given more weight. Duration of the shortest run, for $q=-1/2$, with a sequential code on a 3.16 GHz Intel Core Duo processor is 25 min. 19 sec. The computations do not reveal any clear dependence of the efficiency of AICDM on the power $q$ — in a large interval of $q$ the variation of the number of iterations is just several per cent. The low sensitivity to the value of $q>0$ is probably linked to the relative smallness of the region where $\widehat{f}>\langle\widehat{f}\rangle$.

6.5 An enhanced version of AICDM

The efficiency of AICDM can be further improved without employing new mathematical ideas. For instance, we have explored the possibility of gradual refinement of approximate solutions with increasing spatial resolution. The test problem, considered in the two previous subsections, was solved with the resolution of $(16M)^{3}$ Fourier harmonics, where the integer $M$ is successively increased from 1 to 4. When seeking a solution with the resolution $64^{3}$ harmonics, which satisfies the accuracy requirement $d(\eta)<10^{-\kappa}$, the intermediate-stage computations with the resolution of $(16M)^{3}$ harmonics are terminated as soon as AICDM finds a solution to the accuracy

[TABLE]

We call this algorithm enhanced AICDM, or EAICDM. Durations of runs for integer $\kappa\geq 3$ are shown in Table 2 (we do not present total numbers of evaluations of the determinant of the Hessian in runs by EAICDM, since the times of their computations vary significantly with the resolution).

A further acceleration is likely to be achieved by avoiding dealiasing, but we did not explore this possibility.

7 Concluding remarks

We have presented new forms of the Monge–Ampère equation in $R^{N}$: the second-order divergence (4), Fourier integral (11) and convolution (21) forms. They have been derived under the assumption that the MAE (1) has a r.h.s. with a non-vanishing spatial mean, and hence (1) admits solutions (7). The first form gives an opportunity to relax the regularity requirements for weak solutions to the MAE, such that the integral in the l.h.s. of the identity (5) is well-defined. The third form is an integro-differential equation of the first order, which might be useful for development of an algorithm for computation of a solution (7) using transformations of spatial variables. This paper is mostly concerned with the Fourier integral form of the MAE used to prove existence of a small-amplitude solution of the form (7) for a weakly varying r.h.s. in (1) and to formulate an algorithm for computation of a solution (7) to (1) in an odd-dimensional space. In a test application to a three-dimensional MAE with the r.h.s. reminiscent of mass transportation problems considered in cosmology, we have demonstrated that a solution to the MAE with a smooth positive r.h.s. can be efficiently obtained by two versions of this algorithm, ACPDM (the algorithm with continuation in a parameter and discrepancy minimisation) and AICDM (the algorithm with improvement of convexity and discrepancy minimisation). While the latter is suitable for computation of a convex solution for an everywhere positive r.h.s. $f$, the former requires only $\langle f\rangle\neq 0$. However, modification of ACPDM for the case of a zero-mean space-periodic r.h.s. is straightforward.

There are some analogies between our method and the inexact Newton–Krylov solver with preconditioning, used in [21]. In our algorithms, solving (27) or (32), where the easily invertible polynomial part of the MAE is separated out, can be regarded as preconditioning; we believe that it is optimal. Indeed, its design is based on the algebraic nature of the MAE and primarily on the positivity and bounds of the kernels in the Fourier integral form, a very special — and so far not reported — property of the MAE. Furthermore, our algorithms are of the Krylov type. However, we do not rely on Newton iterations: instead of solving a succession of Newton problems, each with a complexity comparable to that of the MAE, we tackle the MAE directly.

In this paper, space periodicity and Fourier decompositions play an important rôle in two respects: First, we rely on the Fourier methods to find the optimal coefficients (28) in the equations (27) or (32). Note that space periodicity is not required for the derivations — we consider Fourier integrals, and not Fourier sums, and the algorithms do not require computation of these integrals. Second, our algorithms are realised in the spectral form, because we seek space-periodic solutions. In a periodicity domain the required computations of the inverse Laplacian, as well as numerical differentiation become trivial, when spectral methods are applied. Recall that Fourier methods have also the advantage of being more accurate: when the resolution is increased, Fourier series converge to a solution together with all the derivatives that the solution possesses, whilst finite differences do not approximate derivatives beyond their fixed order. We would also use the spectral approach, for instance, to solve numerically Dirichlet or Neumann boundary value problems for $u^{\prime}$ (see (7)) in a rectangular box.

However, the geometry of the domain, where the solution is sought, may prohibit the use of spectral methods. What happens when the MAE must be solved for more complicated boundary conditions in regions whose geometry is not rectangular (parallelepiped-shaped)? Our algorithms do not inherently rely on spectral methods and will remain applicable. Since we are solving the MAE in terms of the Laplacian of the unknown function, we would have to invert the Laplacian with suitable regular boundary conditions. Finite differences can be applied in conjunction with our method to carry out this task and for computation of derivatives involved in the function $F$ appearing in the r.h.s. of (32) and (32). Key equations for our algorithms, (27) and (32), are formulated in the physical space and can be used in conjunction with our method for a variety of boundary conditions. Their convergence properties remain of course to be investigated — for instance, whether the theoretically predicted values (28) will still be the optimal choice for the coefficients of the polynomial in the l.h.s. of (32).

An attractive feature of our algorithm is its simplicity: our Fortran-95 source code realising EAICDM is below 18 KB (672 lines long, not including the source for the Fast Fourier Transform).

Hereafter we list some open questions. Under which conditions does the generalised MAE

[TABLE]

with the coefficients (28), which is solved for a set of $p$ by ACPDM, possess zero-mean space-periodic solutions for all $0<p<1$? Does its discretisation always have a solution, as long as the generalised MAE itself does? Does the solution of (39) depend analytically on the parameter $p$, and hence polynomial extrapolation for $p=1$, that we employ, is mathematically sensible, or should another asymptotics near $p=1$ be assumed? What is the optimal choice of the sequence of values of $p$ to be used by ACPDM? Which scalar product is optimal for acceleration of convergence of stabilised iterations? Can Chebyshev techniques [32] be used to improve efficiency of the iterative processes?

Additional questions emerging outside the main topic of this paper — numerical methods for solution of the MAE — also cannot be avoided, since any information concerning the structure of solutions to Monge–Ampère problems of the cosmological type with a high-contrast r.h.s., formulated in Section 5, can be incorporated into specialised solvers in order to improve their performance (like it has proved possible to accelerate computations about 4 times just by taking into account in AICDM convexity of solutions to the MAE with a positive r.h.s.). The questions are: What is the asymptotics of solutions in the small parameter $\delta$ determining the width of the localised objects? Figures 4 and 5 show that the space is divided into regions of dominant influence of each object. Also, notable are almost axisymmetric structures in the Laplacian of the solution around the lines connecting centres of objects, seen on Figure 3. Can these geometric features be identified by performing the asymptotic analysis of the problem? How does the contrast number measure numerical complexity of the MAE and, in particular, the condition number of the linearisation near the solution (for a given spatial discretisation)?

Acknowledgments. We are grateful to S. Colombi and A. Sobolevskii for discussions. Computations have been carried out on the computer “Mésocentre SIGAMM (Simulations Interactives et Visualisation en Géophysique, Astronomie, Mathématique et Mécanique)” hosted by Observatoire de la Côte d’Azur, France. Research visits of VZ and OP to the Observatoire de la Côte d’Azur in the autumns of 2007 and 2008 were supported by the French Ministry of Education. The authors were partially financed by the grant ANR-07-BLAN-0235 OTARIE from Agence nationale de la recherche (France). VZ and OP were also partially supported by the grant 07-01-92217-CNRSL_a from the Russian foundation for basic research.