A Gradually Reinforced Sample-Average-Approximation Differentiable Homotopy Method for a System of Stochastic Equations

This paper is concerned with the problem of finding an $x\in\mathbb{R}^{n}$ satisfying

where $\xi$ is a stochastic parameter in a probability space $\Omega$, $f(\cdot,\xi):\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is a continuously differentiable mapping and $\mathbb{E}[\cdot]$ denotes the expected value over $\Omega$. As an effective paradigm for addressing risks and uncertainty, the SSE (2) can be regarded as a natural extension of a system of equations (SE) to a stochastic environment and has been extensively studied and applied in a number of areas including economics, game theory, management science and engineering. An excellent review on these applications can be found in [43]. Furthermore, the SSE is closely associated with several other stochastic problems such as stochastic variational inequalities (SVI) [6, 7, 17, 29, 42], stochastic nonlinear complementarity problems (SNCP) [5, 18, 31] and stochastic optimization (SO) [12, 27, 30, 39].

To solve the deterministic SE, several well-known methods have been proposed in the literature, which include Newton methods, homotopy or path-following methods and their variants. However, the convergence of Newton methods very much depends on the starting points. The homotopy methods are a class of powerful methods for solving SEs and possess the desired property of global convergence, which play an inevitable and lasting role in various fields; see, for instances, [11, 16, 28, 53]. The homotopy methods can be classified as simplicial homotopy methods and differentiable homotopy methods. The simplicial homotopy methods, dating back to the seminal paper by [14, 44], are powerful mechanisms for solving equations with highly nonlinearities or non-smooth structures, but they cannot benefit much from the differentiability of problems and can be extremely time-consuming especially when the problem sizes are large. The differentiable homotopy methods, introduced in [33], perfectly overcome the deficiency of the simplicial homotopy methods and are capable of following a smooth path.

Over the past few decades, the differentiable homotopy methods have been substantively investigated in the literature and applied to various areas such as general equilibrium theory and game theory for solving the problems that can be reformulated as systems of differentiable equations. Two direct proofs were given in [20] to show the feasibility of linear tracing procedure, which was made differentiable for computing Nash equilibria in [22]. A stochastic version of linear tracing procedure was established in [24] for the computation and selection of stationary equilibria or stationary Markov perfect equilibria in stochastic games. Later, a more efficient homotopy method based on the idea of interior-point method was developed in [11]. An all-solution differentiable homotopy method was proposed in [32] to find all equilibria for static and dynamic games with continuous strategies. A convex quadratic penalty differentiable homotopy method was developed in [8] to compute perfect equilibria for large scale noncooperative games. A differentiable homotopy method was described in [45] to solve general equilibrium models with incomplete asset markets and its reliability was evidenced by numerous numerical experiments. A generically convergent differentiable homotopy method was constituted in [21] to compute an equilibrium for a finance version of the general equilibrium problem considered in [45]. A proximal block coordinate homotopy framework was presented in [50] to solve large scale sparse least squares problems through numerically following a piecewise smooth path. A polyhedral homotopy-baed method was designed in [36] to find solutions to generalized Nash equilibrium problems. More differentiable homotopy methods and their applications can be found in the literature such as [23, 35, 37, 51, 54] and the references therein.

Unfortunately, all the existing methods fail to directly solve a stochastic system since the underlying mapping $\mathbb{E}_{\xi}[f(x,\xi)]$ of an SSE is in a form of an expected value with respect to a certain stochastic parameter vector $\xi$ and needs to be evaluated first. As demonstrated in [30], the evaluation of $\mathbb{E}_{\xi}[f(x,\xi)]$ is generally a tough job, because the distribution of the stochastic parameter is usually unknown and can only be simulated with historical data. Even if the distribution is given, computing multidimensional integrals is very costly. To deal with these difficulties, two classes of methods have been proposed in the literature: stochastic approximation (SA) methods and simulation based approaches. It was proved in [41] that the SA method is almost surly convergent to a solution of an SSE, only when $f$ satisfies certain conditions and the samples and stepsize are suitably chosen. Notwithstanding, this method is very sensitive to the choice of the stepsize at each iteration and sometimes performs poorly in practice. A modified SA method with a better performance was proposed in [40]. This better performance, however, occurs just for a special class of convex stochastic optimization and saddle point problems.

The simulation based approaches are fashionable tools to address the uncertainty. Among them, the sample-average-approximation (SAA) is one of the most popular representatives. The basic idea of the SAA is rather simple: randomly generating samples for the stochastic parameter $\xi$ and approximating the “true” underlying mapping by the average of several deterministic mappings corresponding to these samples. As demonstrated in [26, 34, 46], an appropriate incorporation of the SAA into a numerical algorithm can lead to a reasonable performance in solving general stochastic problems, and an exponential convergence rate of the SAA for the SVI and SNCP has been established under some mild conditions. Inspired by this success, this paper intends to incorporate the SAA into a differentiable homotopy method for a solution to the SSE (2). To apply the SAA in the existing methods, one needs to make a tradeoff between a cheap coarse estimate and an expensive finer estimate. A common practice is to use a gradually reinforced SAA scheme, that is, initially selecting a small sample size and then gradually increasing the sample size to approximate the target expected value. This allows a rapid progress at early stages and reduces the overall compuational cost for finding a desired solution. However, this practice may fail to converge due to the changes in sample size [4].

It is well known that a differentiable homotopy method attains global convergence through a continuous deformation process. This naturally raises the question: can we devise a differentiable homotopy method, which is able to progressively increase the sample size of SAA during the deformation process? The question seems quite challenging and were not considered in the existing work due to the intuitive collision between the discreteness of samples and the continuous deformation of the homotopy methods. To overcome this hurdle, we introduce a sequence of continuously differentiable functions of the homotopy parameter $t$ ranging between zero and one and incorporate with these functions a gradually reinforced sample-average-approximation (GRSAA) scheme into a differentiable homotopy method. As a result of this incorporation, we establish a differentiable homotopy system and reap a GRSAA differentiable homotopy method. As $t$ descends from one to zero, the homotopy system gradually increases the sample size and eventually reaches the desired SAA (4) at $t=0$. The solution set of the homotopy system contains an everywhere smooth path, which starts from an arbitrary point and ends at a solution to the desired SAA or a satisfactory approximate solution to the SSE (2). The GRSAA differentiable homotopy method serves as a bridge to link a gradually reinforced SAA scheme and a differentiable homotopy method and retains the inherent property of global convergence of a standard differentiable homotopy method. To evince the benefit of differentiability, we also present a GRSAA simplicial homotopy method. To make numerical comparisons with a standard differentiable homotopy method and the GRSAA simplicial homotopy method, we have implemented the GRSAA differentiable homotopy method to solve several important applications of the SSE such as the SVI and market equilibrium problems. Numerical results further verify that two main features of the GRSAA differentiable homotopy method, the gradual reinforcement in sample size and differentiability, can significantly enhance the effectiveness and efficiency.

The rest of the paper is organized as follows. In Section 2, we first give a brief review about the gradually reinforced SAA scheme and differentiable homotopy methods and then develop a gradually reinforced sample-average-approximation (GRSAA) differentiable homotopy method to solve the SAA system for the SSE (2). The convergence properties of the proposed method are discussed in Section 3. For numerical comparisons, Section 4 describes a GRSAA simplicial homotopy method, which belongs to a type of non-differentiable homotopy method. In Section 5, we employ the GRSAA differentiable homotopy method to solve some numerical examples and compare the performance of the GRSAA differentiable homotopy method with that of the GRSAA simplicial homotopy method and a standard differentiable homotopy method. We also apply the GRSAA differentiable homotopy method to solve several large-scale SSEs. All the numerical results are reported in Section 5. The paper is concluded in Section 6.

To further evince the advantage of the GRSAA differentiable homotopy method, we describe in this section a gradually reinforced SAA (GRSAA) simplicial homotopy method. As an underlying triangulation of the method, we make use of the $D_{2}$-triangulation of $(0,1]\times\mathbb{R}^{n}$ with continuous refinement of grid size, whose definition and pivot rules are given in [10]. Given $\tau_{0}>0$, let $t_{\ell}=1/(1+\tau_{0}\ell)$, $\ell=0,1,\ldots$. Then, $t_{\ell}$, $\ell=0,1,\ldots$, form a descent sequence with $t_{0}=1$ and $\lim\limits_{\ell\rightarrow\infty}t_{\ell}=0$. Given $\tau_{1}>0$, we define $q_{\ell}=\tau_{1}\ell$, $\ell=1,2,\ldots$. Therefore, as $\ell\rightarrow\infty$, $q_{\ell}\rightarrow\infty$. Given the initial grid size $\varpi_{0}>0$ and $r_{\ell}\in\{1/j\;|\;j=1,2,\ldots\}$, let $\varpi_{\ell+1}=r_{\ell}\varpi_{\ell}$, $\ell=0,1,\ldots$, which is the gird size of a simplex of the $D_{2}$-triangulation in $\{t_{\ell}\}\times\mathbb{R}^{n}$. Let $x^{0}\in\mathbb{X}$ be a given point. We still utilize the continuously differentiable function $\theta_{\ell}(t)$ defined in Section 2 to form a simplicial homotopy system,

where $\bar{d}(t,x)=\theta(t)\dfrac{1}{q_{\ell}}\sum\limits_{i=1}^{q_{\ell}}f(x,\xi_{i})+(1-\theta(t))\dfrac{1}{q_{\ell+1}}\sum\limits_{i=1}^{q_{\ell+1}}f(x,\xi_{i})$ for $t_{\ell+1}\leq t\leq t_{\ell}$. Clearly, when $t=1$, the system (11) has a unique solution $x^{0}$. As follows, we exploit the homotopy mapping $g(t,x)$ to attain a GRSAA simplicial homotopy method.

Let $v^{c}=(0,x^{0}-\frac{\varpi_{0}}{n+1}e)$ with $e=(1,1,\ldots,1)^{\top}\in\mathbb{R}^{n}$. Let $D_{2}$ be the collection of all simplices of the $D_{2}$-triangulation of $(0,1]\times\mathbb{R}^{n}$ after the translation of $v^{c}$, where the translation ensures that $(1,x^{0})$ is in the interior of a unique simplex in $\{1\}\times\mathbb{R}^{n}$. We denote a $q$-dimensional simplex in $D_{2}$ by $\sigma=\langle v^{0},v^{1},\ldots,v^{q}\rangle$, where $v^{i}$ is a vertex of $\sigma$ for $i=0,1,\ldots,q$.

Following a similar argument to that in [8], we know that there is a unique complete simplex contained in $\{1\}\times\mathbb{R}^{n}$. Given these notations, a simplicial homotopy method for computing an approximate solution to the SSE can be stated as follows.

Step 0:

Let $\delta_{0}$ be a given sufficiently small positive number. Let $\tau_{0}=\langle v^{0},v^{1},\ldots,v^{n}\rangle$ be the unique complete simplex in $\{1\}\times\mathbb{R}^{n}$ with $v^{0}=(1,x^{0})$ and $\sigma_{0}$ be the unique $(n+1)$-dimensional simplex in $[t_{1},1]\times\mathbb{R}^{n}$ with $\tau_{0}$ as its facet. Let $v^{+}$ be the vertex of $\sigma_{0}$ opposite to $\tau_{0}$. Let $\ell=0$ and go to Step 1.

Step 1:

Perform a linear programming step to bring $\left(\begin{array}[]{c}1\\ g(v^{+})\end{array}\right)$ into the system, $\sum\limits_{l=0}^{n}\zeta_{l}\left(\begin{array}[]{c}1\\ g(v^{l})\end{array}\right)=\left(\begin{array}[]{c}1\\ 0\end{array}\right).$ Suppose that $\left(\begin{array}[]{c}1\\ g(v^{l})\end{array}\right)$ leaves the system. Let $\tau_{\ell+1}$ be the facet of $\sigma_{\ell}$ opposite to $v^{l}$ and $v^{l}=v^{+}$. Go to Step 2.

Step 2:

If $\tau_{\ell+1}\subset\{t_{\ell+1}\}\times\mathbb{R}^{n}$ and $t_{\ell+1}\leq\delta_{0}$, the method terminates. Otherwise, let $\sigma_{\ell+1}$ be the unique simplex that shares together with $\sigma_{\ell}$ a common facet $\tau_{\ell+1}$. Let $v^{+}$ be the vertex of $\sigma_{\ell+1}$ opposite to $\tau_{\ell+1}$ and $\ell=\ell+1$. Go to Step 1.

As follows, we show that the simplicial path generated by the GRSAA simplicial homotopy method leads to a solution to the SSE (2). Let $g^{-1}(0)=\{(t,x)\in[0,1]\times\mathbb{R}^{n}\;|\;g(t,x)=0\}.$ Then it follows from Assumption 1 that $g^{-1}(0)$ is a compact set. As a result of Definition 1 and the continuity of $f$, it is easy to verify that all the complete simplices are contained in a bounded set of $[0,1]\times\mathbb{R}^{n}$. Under the nondegeneracy assumption, following the same argument as that in Todd [48], one can derive that all the simplices generated by the GRSAA simplicial homotopy method are different, that is, no cycle will occur. Let $\tau_{\ell}$, $\ell=1,2,\ldots$, be the complete simplices generated by the method. Let $\bar{t}_{\ell}$ be the smallest value in $(0,1]$ such that $\tau_{\ell}\subset(0,\bar{t}_{\ell}]\times\mathbb{R}^{n}$. Since all $\tau_{\ell}$, $\ell=1,2,\ldots$, are contained in a bounded set, we must have $\bar{t}_{\ell}\rightarrow 0$ as $\ell\rightarrow\infty$.

Let $\sigma=\{t\}\times\langle x^{0},x^{1},\ldots,x^{n}\rangle$ be an $n$-dimensional simplex in $\{t\}\times\mathbb{R}^{n}$. For $y=\sum\limits_{l=0}^{n}\zeta_{l}(y^{l})$ with $\zeta_{l}\geq 0$ and $\sum\limits_{l=0}^{n}\zeta_{l}=1$, a piecewise linear approximation of $g(t,y)$ is given by $\bar{g}(t,y)=\sum\limits_{l=0}^{n}\zeta_{l}g(t,y^{l}).$

Since $g(t,y)$ is uniformly continuous on $[0,1]\times\mathbb{R}^{n}$, there is a constant $\rho_{0}>0$ such that $\|g(t,y)-g(\hat{t},\hat{y})\|\leq\rho_{0}\|(t,y)-(\hat{t},\hat{y})\|$ for any $(t,y)$ and $(\hat{t},\hat{y})$ in $[0,1]\times\mathbb{R}^{n}$.