# An Efficient Approach Based on the Gradient Definition for Solving Conditional Nonlinear Optimal Perturbation

Authors: Bin Mu; Juhui Ren; Shijin Yuan
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Conditional nonlinear optimal perturbation (CNOP) has been widely applied to study the predictability of weather and climate. The classical method of solving CNOP is adjoint method, in which the gradient is obtained using the adjoint model. But some numerical models have no adjoint models implemented, and it is not realistic to develop from scratch because of the huge amount of work. The gradient can be obtained by the definition in mathematics; however, with the sharp growth of dimensions, its calculation efficiency will decrease dramatically. Therefore, the gradient is rarely obtained by the definition when solving CNOP. In this paper, an efficient approach based on the gradient definition is proposed to solve CNOP around the whole solution space and parallelized. Our approach is applied to solve CNOP in Zebiak-Cane (ZC) model, and, compared with adjoint method, which is the benchmark, our approach can obtain similar results in CNOP value and pattern aspects and higher efficiency in time consumption aspect, only 12.83 s, while adjoint method spends 15.04 s and consumes less time if more CPU cores are provided. All the experimental results show that it is feasible to solve CNOP with our approach based on the gradient definition around the whole solution space.


---

## Body

## 1. Introduction

In the study of weather and climate predictability, it is crucial to determine the fastest growing perturbation. To solve the fastest growing perturbation in a nonlinear system, Mu and Duan [1] proposed the concept of conditional nonlinear optimal perturbation (CNOP), which can represent the nonlinear initial perturbations that satisfy certain constraint conditions and result in the largest nonlinear evolution at the prediction time. Later, Mu et al. [2] extended the CNOP method to study the optimal parameter perturbation. The CNOP method has been widely applied to study the predictability of many phenomena and many research fields related to initial errors and model parameter errors, such as EI Niño-Southern Oscillation (ENSO) event [3–5], Kuroshio large meander, [6] and grassland ecosystem [7]; spring predictability barrier [8–11]; targeted observation of the atmosphere and ocean [12–15]; ensemble forecast [16–18]. It is obvious that the CNOP plays an important role in the study of weather and climate predictability.Solving CNOP is essentially an optimization problem of nonlinear objective function. In the current study, the approaches to solve CNOP can be classified into two types depending on whether the gradient is used. The gradient-based approaches solve CNOP by searching the optimum value along the direction of gradient descent, such as the spectral projected gradient (spg2) [19], the limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) [20], and the sequential quadratic programming (SQP). While the gradient-free approaches get the optimum value by searching randomly around the whole or partial solution space to solve CNOP, such as intelligent algorithm (IA). The random search method initially was applied to ideal numerical models with 3 to 20 dimensions [21, 22]. To apply this method to solve efficiently CNOP for practical models, some researchers reduced the high-dimensional solution space to a relative low one firstly and then employed IAs to solve CNOP in the reduced low-dimensional solution space [23–32]. They got good results when taking ZC model with 1080 dimensions as a case. But there inevitably exists the information loss because of dimension reduction. Therefore, the gradient descent method is widely used to calculate CNOP.Generally, the gradient is obtained using the adjoint model corresponding to the numerical model, which is referred to as the adjoint method. However, there are no corresponding adjoint models implemented for some numerical models [33–35], and it is not realistic to develop the corresponding adjoint model from scratch due to the huge amount of work. Wang and Tan [36, 37] attempted to obtain approximate gradient information based on ensemble technique, in which they employed the samples ensemble of initial perturbations and corresponding prediction increments to denote the approximate tangent linear matrix in the gradient formula of objective function. And a localization technique was introduced to ameliorate the spurious correlation between the two ensembles, in which the localization radius was achieved from artificial experience. This method calculated the gradient information only once during the whole optimization; therefore, it can obtain an approximate CNOP easier and more efficiently than adjoint-based method, but it depends on artificial experience. In fact, the gradient information can also be obtained using gradient definition, but the calculation efficiency will decrease dramatically with the sharp growth of the dimension. And in general, the dimensions of the numerical models for climate and weather are relatively high, which results in the fact that gradient definition method has rarely been applied in solving CNOP. At present, only Chen et al. [38] calculated the gradient in one way that is similar to the gradient definition, but this gradient is calculated in the feather space generated by dimension reduction. Firstly, they reduced the dimensions to the feather space using singular value decomposition (SVD) and represented the initial perturbations using the linear combination of base vectors. Consequently, the objective function was transformed into the function of the linear combination coefficients. The gradient was approximated by the differences, the linear combination coefficients, and prediction increments of the initial perturbations. In other words, the gradient calculated was formally the same as the definition of the gradient, but the small amount in the gradient definition equation is the increment of the coefficient, not the increment of the initial perturbations. This method can obtain an approximate CNOP, and time efficiency depends on the number of base vectors chosen.In this paper, an efficient approach based on the gradient definition around the whole solution space is proposed to solve CNOP, in which some parallel strategies are adopted to improve the calculation efficiency of gradient. In our approach, the gradient calculated is the gradient of objective function with respect to initial perturbation, and we solve the CNOP around the whole solution space, so the CNOP is more accurate. In addition, certain parallel strategies make our approach more efficient than adjoint method. Taking the ZC model as an example, which is a medium-complexity model to forecast ENSO event, our approach is applied to solve CNOP of ENSO event, and the experimental results show that our approach is feasible from CNOP value, CNOP pattern and time efficiency aspects.The remainder of this paper is organized as follows. The detailed introduction of ZC model and the concept of the CNOP are in Section2. Our efficient approach based on the gradient definition around the whole solution space accompanied by parallel strategies is described in Section 3. In Section 4, we employ ZC model as a case study and apply our approach to study the optimal precursor of an ENSO event; here we also show the results and compare the CNOP value, CNOP pattern, and time consumption with those from adjoint method. Finally, we summarize our conclusions and future works in Section 5.


## 2. Zebiak-Cane Model and CNOP

### 2.1. Zebiak-Cane (ZC) Model

ZC model is adopted as the case to verify the feasibility of our approach in solving CNOP. ZC model was developed to simulate and study EI Niño-Southern Oscillation (ENSO) phenomenon; it is a medium-complexity model. The model can calculate perturbations about a climatological mean state that is specified from observation [39]. It can also reproduce the warm events that possess a 3-4 years’ period of oscillation without anomalous external forcing, which is consistent with the real ENSO cycle. ZC model is a coupled atmosphere-ocean model, which has three components: the atmosphere, the ocean, and the coupling component.


#### 2.1.1. Atmosphere Component

The dynamics of atmosphere component follow the Gill model, which is described by the linear shallow-water equations on an equatorial beta plane. The circulation in atmosphere component is forced by a heating anomaly that depends on the sea surface temperature (SST) anomalies and moisture convergence. The atmospheric grid used in the atmosphere component lies in the region 101.25°E–73.125°W, 29°S–29°N.


#### 2.1.2. Ocean Component

The dynamics of the ocean component begin with the linear reduced-gravity model, which can successfully simulate the changes of thermocline depth anomalies and sea surface pressure during EI Niño events. In the ocean component, the surface intensification of wind-driven currents in the real ocean is simulated by the shallow frictional layer. The component can simulate the mean features of SST anomaly (SSTA) forced by ENSO composite wind anomalies. The oceanic gird used in the ocean component lies in the region 124°E–80°W, 28.75°S–28.75°N.


#### 2.1.3. Coupling Component

In the coupling component, the atmosphere component retains steady state and was run with certain monthly mean SSTA to simulate wind anomalies in advance. The ocean component is driven by the surface wind stress anomalies, which are produced by combining the background mean winds and surface wind anomalies generated by the atmosphere component. After coupling the ocean and atmosphere component, the region of the coupled model is shown in Figure1. The rectangle with black solid line represents the region of atmosphere component. The rectangle with black dash line represents the region of ocean component. The rectangle with red dash line represents the integration region of SSTA in the coupled model.Figure 1
The region of ZC model.


### 2.2. CNOP

CNOP can represent the initial perturbation subjected to a given physical constraint and result in the largest nonlinear evolution at the prediction time in the nonlinear weather and climate models. Suppose we have the following model:(1)∂X∂t+FX=0,Xt=0=X0, where X is an n-dimensional state vector of the model, while X0 is the n-dimensional initial state vectors at initial time (t = 0), and F is a nonlinear partial differential operator. The discrete form of (1) can be described as follows:(2)Xt=Mt0→tX0, where M is a nonlinear propagation operator, t0 and t are separately the initial optimization time and terminal time, Xt is the value of X at time t, and Xt represents the development of X0 from time t0 to t.The CNOP, represented usingx0δ, is the solution of the following optimization problem:(3)Jx0δ=maxx0≤δJx02,Jx0=Mt0→tX0+x0-Mt0→tX0, where x0 is the n-dimensional initial perturbation of X0 and δ is the constraint radius of the initial perturbation. Jx02 is the objective function.Obviously, solving (3) is to solve an optimization problem. So CNOP can be obtained by a nonlinear optimization algorithm. Generally, optimization algorithms, such as LBFGS, SQP, and SPG2, are designed to find the minimum value of the objective function. In this paper, the SPG2 algorithm is employed to solve CNOP. The SPG2 method is often applied to solve the problem of the following form:(4)minfxsubjecttox∈Ω, where Ω is a closed convex set in IRn. To use SPG2 algorithm to solve CNOP directly, we let Jox0δ=-J(x0δ); then the optimization problem in (3) is equivalent to the following optimization problem:(5)Jox0δ=minx0≤δ-Jx02.Now, the optimization problem has been converted into solving minimum value of the objective function, which is the same form with problem in (4); therefore, we can use SPG2 method directly.


## 2.1. Zebiak-Cane (ZC) Model

ZC model is adopted as the case to verify the feasibility of our approach in solving CNOP. ZC model was developed to simulate and study EI Niño-Southern Oscillation (ENSO) phenomenon; it is a medium-complexity model. The model can calculate perturbations about a climatological mean state that is specified from observation [39]. It can also reproduce the warm events that possess a 3-4 years’ period of oscillation without anomalous external forcing, which is consistent with the real ENSO cycle. ZC model is a coupled atmosphere-ocean model, which has three components: the atmosphere, the ocean, and the coupling component.


### 2.1.1. Atmosphere Component

The dynamics of atmosphere component follow the Gill model, which is described by the linear shallow-water equations on an equatorial beta plane. The circulation in atmosphere component is forced by a heating anomaly that depends on the sea surface temperature (SST) anomalies and moisture convergence. The atmospheric grid used in the atmosphere component lies in the region 101.25°E–73.125°W, 29°S–29°N.


### 2.1.2. Ocean Component

The dynamics of the ocean component begin with the linear reduced-gravity model, which can successfully simulate the changes of thermocline depth anomalies and sea surface pressure during EI Niño events. In the ocean component, the surface intensification of wind-driven currents in the real ocean is simulated by the shallow frictional layer. The component can simulate the mean features of SST anomaly (SSTA) forced by ENSO composite wind anomalies. The oceanic gird used in the ocean component lies in the region 124°E–80°W, 28.75°S–28.75°N.


### 2.1.3. Coupling Component

In the coupling component, the atmosphere component retains steady state and was run with certain monthly mean SSTA to simulate wind anomalies in advance. The ocean component is driven by the surface wind stress anomalies, which are produced by combining the background mean winds and surface wind anomalies generated by the atmosphere component. After coupling the ocean and atmosphere component, the region of the coupled model is shown in Figure1. The rectangle with black solid line represents the region of atmosphere component. The rectangle with black dash line represents the region of ocean component. The rectangle with red dash line represents the integration region of SSTA in the coupled model.Figure 1
The region of ZC model.


## 2.1.1. Atmosphere Component

The dynamics of atmosphere component follow the Gill model, which is described by the linear shallow-water equations on an equatorial beta plane. The circulation in atmosphere component is forced by a heating anomaly that depends on the sea surface temperature (SST) anomalies and moisture convergence. The atmospheric grid used in the atmosphere component lies in the region 101.25°E–73.125°W, 29°S–29°N.


## 2.1.2. Ocean Component

The dynamics of the ocean component begin with the linear reduced-gravity model, which can successfully simulate the changes of thermocline depth anomalies and sea surface pressure during EI Niño events. In the ocean component, the surface intensification of wind-driven currents in the real ocean is simulated by the shallow frictional layer. The component can simulate the mean features of SST anomaly (SSTA) forced by ENSO composite wind anomalies. The oceanic gird used in the ocean component lies in the region 124°E–80°W, 28.75°S–28.75°N.


## 2.1.3. Coupling Component

In the coupling component, the atmosphere component retains steady state and was run with certain monthly mean SSTA to simulate wind anomalies in advance. The ocean component is driven by the surface wind stress anomalies, which are produced by combining the background mean winds and surface wind anomalies generated by the atmosphere component. After coupling the ocean and atmosphere component, the region of the coupled model is shown in Figure1. The rectangle with black solid line represents the region of atmosphere component. The rectangle with black dash line represents the region of ocean component. The rectangle with red dash line represents the integration region of SSTA in the coupled model.Figure 1
The region of ZC model.


## 2.2. CNOP

CNOP can represent the initial perturbation subjected to a given physical constraint and result in the largest nonlinear evolution at the prediction time in the nonlinear weather and climate models. Suppose we have the following model:(1)∂X∂t+FX=0,Xt=0=X0, where X is an n-dimensional state vector of the model, while X0 is the n-dimensional initial state vectors at initial time (t = 0), and F is a nonlinear partial differential operator. The discrete form of (1) can be described as follows:(2)Xt=Mt0→tX0, where M is a nonlinear propagation operator, t0 and t are separately the initial optimization time and terminal time, Xt is the value of X at time t, and Xt represents the development of X0 from time t0 to t.The CNOP, represented usingx0δ, is the solution of the following optimization problem:(3)Jx0δ=maxx0≤δJx02,Jx0=Mt0→tX0+x0-Mt0→tX0, where x0 is the n-dimensional initial perturbation of X0 and δ is the constraint radius of the initial perturbation. Jx02 is the objective function.Obviously, solving (3) is to solve an optimization problem. So CNOP can be obtained by a nonlinear optimization algorithm. Generally, optimization algorithms, such as LBFGS, SQP, and SPG2, are designed to find the minimum value of the objective function. In this paper, the SPG2 algorithm is employed to solve CNOP. The SPG2 method is often applied to solve the problem of the following form:(4)minfxsubjecttox∈Ω, where Ω is a closed convex set in IRn. To use SPG2 algorithm to solve CNOP directly, we let Jox0δ=-J(x0δ); then the optimization problem in (3) is equivalent to the following optimization problem:(5)Jox0δ=minx0≤δ-Jx02.Now, the optimization problem has been converted into solving minimum value of the objective function, which is the same form with problem in (4); therefore, we can use SPG2 method directly.


## 3. The Efficient Approach Based on Gradient Definition

### 3.1. The Approach Based on Gradient Definition

The primary idea of our approach is to calculate the gradient of objective function with respect to the initial perturbationx0 using the gradient definition in mathematics firstly and then to apply spg2 method to solve CNOP based on the gradient information. In our case, the optimization problem is described in (5); obviously, it can be described as follows:(6)Fx0=-fx02,fx0=Jx0.In mathematics, the gradient is a generalization for the usual concept of derivative of a function in one dimension to a function in several dimensions. So the gradient of function in (6) is as follows:(7)gradF=-2fx0×f′x0, where f′x0 denotes the first-order partial derivative of function f(x0), namely, the gradient of function f(x0). Due to x0 being an n-dimensional perturbation, let x0=(x01,x02,x03,…,x0i,…,x0n), according to the definition of gradient in mathematics; then the gradient for function f(x0) in a rectangular coordinate systemic is(8)gradf=∂f∂x01e1+∂f∂x02e2+⋯+∂f∂x0iei+⋯+∂f∂x0nen, where gradf represents f′(x0), n is a positive nonzero integer, i=1,2,3,…,n, and ei are the orthogonal unit vectors pointing in the coordinate directions. Therefore, for a certain point a1,…,ai,…,an, the partial derivative of f at direction x0i is as follows:(9)∂f∂x0ia1,…,an=lim∇→0⁡fa1,…,ai+∇,…,an-fa1,…,ai,…,an∇, where ∇ is a real number which should approach 0 but never equals 0. We will provide detailed description on the setting the value of ∇ in Section 4. Using (9), once ∇ is determined, it becomes much easier to calculate the derivative of f at a certain direction for a certain point. In (8) and (9), if the dimension of the variable in a model is n-dimensional, the gradient vector is also n-dimensional.Algorithm1 shows the pseudocode of our approach. There are two main parts in the approach. First, we initialize the related parameters used in our approach; the meaning of parameters ∇,δ,n,x0 has been shown in the above-mentioned equations (3) and (9); the value of ∇ is determined in Section 4.1; δ is set as (1) in this paper. Then we use SPG2 algorithm to calculate CNOP, the maximum iteration steps are set as 20 for stopping criterion, thegradient(),values, andline_search represent related subroutines, thegradient subroutine calculates the gradient by implementing formulas (8) and (9), thevalues subroutine calculates the value of objective function in current position, andline_search subroutine searches the next position along the direction of gradient decent. Eventually, the program will output the CNOP as the result.Algorithm 1:The pseudocode of our approach.
Initialization:
(1) Set the parameters ∇,δ,n,x0
SPG2:
(2) Calculate the gradient of x0 with respect to the objective function using subroutine gradient(x0)
(3) Calculate the value of objective function in x0 using subroutine values(x0)
(4) While (the stopping criterion is not satisfied) do
(5) Calculate the new position x′ using subroutine line_search
(6) Calculate the gradient in x′ using subroutine gradient(x′)
(7) End while
Output: CNOP (thex′ when the value ofvalues(x′) is the minimum for all x′)


### 3.2. Parallel Strategies

As shown in Figure1, the outer rectangle with black solid line represents the region of atmospheric component with latitudinal resolution dx=5.625° and longitudinal resolution dy=2°. The middle rectangle with red dashed line denotes the region of the ocean component with the resolutions dx=2° and dy=0.5°, which forms a 30∗34 grid. After removing the unused marginal area, the inner rectangle with blue dashed line is the integration region of the SSTA with resolutions dx=5.625° and dy=2°, which forms a 20∗27 grid. When studying ENSO phenomenon, the two physical variables in ZC model involved in objective function are SSTA and thermocline height anomalies (THA). Thus, the dimension of ZC model is 1080 (20∗27∗2) after combining the two variables into one vector.Taking the ZC model with 1080 dimensions as an example, we implement serially our approach to solve CNOP descried above on TH-1A supercomputer system at National Supercomputer Center in Tianjin. The available resources for us are as follows: 20 available nodes, each node with two Intel Xeon X5670 processors at 2.93 GHz and 24 GB memory, total 240 CPU cores. We measure the time consumption of our serial approach with Intel VTune Amplifier, which is shown in Figure2; it costs 1482.069 s for a complete run., in which subroutinegradient occupies 99.9% of the entire time. We can conclude that the time consumption of subroutinegradient will dramatically increase with the increasing dimensions. Therefore, improving time efficiency of subroutinegradient is crucial and necessary. In this section, certain parallel strategies are designed.Figure 2
Time distribution in computation. The colored bar represents the distribution of the wait time according to the utilization levels (Idle, Poor, Ok, Ideal, and Over) defined by the VTune Amplifier XE. The longer the bar, the higher the value.Considering the dependence between the current iteration and next iteration, and the independency between every dimension of a certain gradient vector, we can perform parallel our approach when calculating one gradient vector in one iteration. To ensure the transportability and usability of the parallel strategy, we adopt MPI to realize parallelization on the cluster.To calculate in parallel, the gradient vector is divided into groups which are executed in parallel. The following is the way to decompose the gradient vector into groups: suppose we employm-processes to calculate one gradient vector concurrently; then we should divide the gradient vector into m-groups, and the size of each group for process ip is (10)g=1080m,ip∈1,2,3,…,m;whenr=0G2=1080m+1,ip≤r,ip≠01080m,ip>rorip=0,ip∈1,2,3,…,m;whenr≠0, where m represents the total number of processes used to compute the gradient, r represents the remainder from dividing 1080 by m, g1 and g2 represent the size of one group for different value of r, and ip represents the process number. The process ip calculates one group of one gradient vector; the dimensions for one gradient vector calculated by process ip can be described as follows:(11)g1∗ip-1~g1∗ip-1+1,ip≠0g1∗m-1+1~1080,ip=0;whenr=0g1∗ip-1~g1∗ip-1+1,ip≤r,ip≠0g1∗ip-1~g1∗ip-1+1+r,ip>r;whenr≠0g1∗m-1+1~1080,ip=0.The different parts of one gradient vectors calculated by different processes are collected as one whole gradient vector via the communication mechanism between processes which is implemented with MPI, specifically MPICH and the Intel compiler. The communication mechanism adopted is master-slave mode, process 0 as the master and others as the slaves. Supposing we usem processes, when calculating one gradient vector, process 1 to m-1 sends their part of gradient vector to process 0, respectively, and process 0 receives the messages from slaves and then combines all messages together into a complete gradient vector.


## 3.1. The Approach Based on Gradient Definition

The primary idea of our approach is to calculate the gradient of objective function with respect to the initial perturbationx0 using the gradient definition in mathematics firstly and then to apply spg2 method to solve CNOP based on the gradient information. In our case, the optimization problem is described in (5); obviously, it can be described as follows:(6)Fx0=-fx02,fx0=Jx0.In mathematics, the gradient is a generalization for the usual concept of derivative of a function in one dimension to a function in several dimensions. So the gradient of function in (6) is as follows:(7)gradF=-2fx0×f′x0, where f′x0 denotes the first-order partial derivative of function f(x0), namely, the gradient of function f(x0). Due to x0 being an n-dimensional perturbation, let x0=(x01,x02,x03,…,x0i,…,x0n), according to the definition of gradient in mathematics; then the gradient for function f(x0) in a rectangular coordinate systemic is(8)gradf=∂f∂x01e1+∂f∂x02e2+⋯+∂f∂x0iei+⋯+∂f∂x0nen, where gradf represents f′(x0), n is a positive nonzero integer, i=1,2,3,…,n, and ei are the orthogonal unit vectors pointing in the coordinate directions. Therefore, for a certain point a1,…,ai,…,an, the partial derivative of f at direction x0i is as follows:(9)∂f∂x0ia1,…,an=lim∇→0⁡fa1,…,ai+∇,…,an-fa1,…,ai,…,an∇, where ∇ is a real number which should approach 0 but never equals 0. We will provide detailed description on the setting the value of ∇ in Section 4. Using (9), once ∇ is determined, it becomes much easier to calculate the derivative of f at a certain direction for a certain point. In (8) and (9), if the dimension of the variable in a model is n-dimensional, the gradient vector is also n-dimensional.Algorithm1 shows the pseudocode of our approach. There are two main parts in the approach. First, we initialize the related parameters used in our approach; the meaning of parameters ∇,δ,n,x0 has been shown in the above-mentioned equations (3) and (9); the value of ∇ is determined in Section 4.1; δ is set as (1) in this paper. Then we use SPG2 algorithm to calculate CNOP, the maximum iteration steps are set as 20 for stopping criterion, thegradient(),values, andline_search represent related subroutines, thegradient subroutine calculates the gradient by implementing formulas (8) and (9), thevalues subroutine calculates the value of objective function in current position, andline_search subroutine searches the next position along the direction of gradient decent. Eventually, the program will output the CNOP as the result.Algorithm 1:The pseudocode of our approach.
Initialization:
(1) Set the parameters ∇,δ,n,x0
SPG2:
(2) Calculate the gradient of x0 with respect to the objective function using subroutine gradient(x0)
(3) Calculate the value of objective function in x0 using subroutine values(x0)
(4) While (the stopping criterion is not satisfied) do
(5) Calculate the new position x′ using subroutine line_search
(6) Calculate the gradient in x′ using subroutine gradient(x′)
(7) End while
Output: CNOP (thex′ when the value ofvalues(x′) is the minimum for all x′)


## 3.2. Parallel Strategies

As shown in Figure1, the outer rectangle with black solid line represents the region of atmospheric component with latitudinal resolution dx=5.625° and longitudinal resolution dy=2°. The middle rectangle with red dashed line denotes the region of the ocean component with the resolutions dx=2° and dy=0.5°, which forms a 30∗34 grid. After removing the unused marginal area, the inner rectangle with blue dashed line is the integration region of the SSTA with resolutions dx=5.625° and dy=2°, which forms a 20∗27 grid. When studying ENSO phenomenon, the two physical variables in ZC model involved in objective function are SSTA and thermocline height anomalies (THA). Thus, the dimension of ZC model is 1080 (20∗27∗2) after combining the two variables into one vector.Taking the ZC model with 1080 dimensions as an example, we implement serially our approach to solve CNOP descried above on TH-1A supercomputer system at National Supercomputer Center in Tianjin. The available resources for us are as follows: 20 available nodes, each node with two Intel Xeon X5670 processors at 2.93 GHz and 24 GB memory, total 240 CPU cores. We measure the time consumption of our serial approach with Intel VTune Amplifier, which is shown in Figure2; it costs 1482.069 s for a complete run., in which subroutinegradient occupies 99.9% of the entire time. We can conclude that the time consumption of subroutinegradient will dramatically increase with the increasing dimensions. Therefore, improving time efficiency of subroutinegradient is crucial and necessary. In this section, certain parallel strategies are designed.Figure 2
Time distribution in computation. The colored bar represents the distribution of the wait time according to the utilization levels (Idle, Poor, Ok, Ideal, and Over) defined by the VTune Amplifier XE. The longer the bar, the higher the value.Considering the dependence between the current iteration and next iteration, and the independency between every dimension of a certain gradient vector, we can perform parallel our approach when calculating one gradient vector in one iteration. To ensure the transportability and usability of the parallel strategy, we adopt MPI to realize parallelization on the cluster.To calculate in parallel, the gradient vector is divided into groups which are executed in parallel. The following is the way to decompose the gradient vector into groups: suppose we employm-processes to calculate one gradient vector concurrently; then we should divide the gradient vector into m-groups, and the size of each group for process ip is (10)g=1080m,ip∈1,2,3,…,m;whenr=0G2=1080m+1,ip≤r,ip≠01080m,ip>rorip=0,ip∈1,2,3,…,m;whenr≠0, where m represents the total number of processes used to compute the gradient, r represents the remainder from dividing 1080 by m, g1 and g2 represent the size of one group for different value of r, and ip represents the process number. The process ip calculates one group of one gradient vector; the dimensions for one gradient vector calculated by process ip can be described as follows:(11)g1∗ip-1~g1∗ip-1+1,ip≠0g1∗m-1+1~1080,ip=0;whenr=0g1∗ip-1~g1∗ip-1+1,ip≤r,ip≠0g1∗ip-1~g1∗ip-1+1+r,ip>r;whenr≠0g1∗m-1+1~1080,ip=0.The different parts of one gradient vectors calculated by different processes are collected as one whole gradient vector via the communication mechanism between processes which is implemented with MPI, specifically MPICH and the Intel compiler. The communication mechanism adopted is master-slave mode, process 0 as the master and others as the slaves. Supposing we usem processes, when calculating one gradient vector, process 1 to m-1 sends their part of gradient vector to process 0, respectively, and process 0 receives the messages from slaves and then combines all messages together into a complete gradient vector.


## 4. Experiments and Results Analysis

To demonstrate the effectiveness, validity, and time efficiency of our approach in solving CNOP, we employ ZC model as a case to study the optimal precursor of an ENSO event. Firstly, we calculate the gradient of objective function using our approach and then use the spg2 method to solve CNOP. The final solution of CNOP is the pattern of initial SSTA and THA that will cause the largest evolution at prediction time in the tropical Pacific, named as SSTA-CNOP and THA-CNOP, which are so-called optimal precursor. We optimize the ZC model for 9-month optimization period for different initial months (from January to December). For every initial month, there are corresponding SSTA-CNOP and THA-CNOP. We compare the results with those obtained from adjoint method which is referred to as the benchmark. Compared with adjoint method, our approach calculates the gradient using the gradient definition.When calculating gradient, the value of∇ is critical. Therefore, in Section 4.1, we conduct many experiments to decide the value of ∇. In the following Sections 4.2 and 4.3, compared with the adjoint method, we show the CNOP calculated by our approach from CNOP value and CNOP pattern aspects to verify its effectiveness and validity and then demonstrate the time consumption and speedup up to 240 CPU cores to verify the time efficiency.


### 4.1. Determination of∇

In Section3, we show the mathematical formula (see (9)) to calculate the gradient of objective function. In the equation, ∇ is a real constant which should approach 0 but never equals 0. For our case, the value of ∇ cannot be too small, because too small ∇ will lead no evolution for numerical model; that is, the limit value in (9) always equals or approaches 0; thus we cannot obtain the correct gradient direction. Therefore, what value ∇ should be in our case? We design the following two schemes to determine the value of ∇: (1) ∇ is constant value for every ai in (9); (2)  ∇ is changing with the value of ai in (9). For the CNOP calculated by different methods, the value of objective function (Jx0) is larger; then the CNOP is much better. So, in this paper, we will take the norm Jx0 to measure the magnitude of CNOP, and we take the value of Jx0 as the evaluating standard for CNOP value.For scheme (1), we conduct lots of experiments to solve CNOP, but it is found that the different initial month is corresponding to a different appropriate ∇ for CNOP. And it requires lots of experiments to determine the most appropriate ∇ for each initial month. Therefore, the conclusion is that the scheme (1) is not feasible in our case.For scheme (2), we let ∇=ai∗10-n(n=1,2,3). When n≥3, ∇ is too small, and CNOP value obtained is very small. When n<0, ∇ is too large to calculate the limit value. When n=1 and n=2, we compare the maximum value of objective function obtained by our approach with the results of adjoint method (shown in Figure 3). In Figure 3, the blue solid line represents the CNOP values calculated by adjoint method for different initial month, which is the baseline, while the red one (∇=ai∗10-1) and green one (∇=ai∗10-2) are the CNOP values using our approach.Figure 3
CNOP values calculated by our approach and adjoint method for different initial month. The blue line represents the values by adjoint method; the red one(∇=ai∗10-1) and green one (∇=ai∗10-2) are the values by our approach.In Figure3, the CNOP value obtained by our approach has similar tendency with the results of adjoint method and the CNOP value for every initial month is less than adjoint method, but the largest difference between them is 5.3 (Table 1), which is acceptable, so we can draw the conclusion that ∇=ai∗10-n(n=1,2) is appropriate for our approach.Table 1
The difference value between CNOP values calculated by adjoint method and our approach.
Initial month
1
2
3
4
5
6
7
8
9
10
11
12
a
i
∗
10
-
1
3.0
2.8
3.1
3.9
1.8
1.7
2.1
1.6
1.9
5.3
1.8
2.1
a
i
∗
10
-
2
3.1
3.1
2.5
3.6
3.1
1.3
1.3
1.1
3.3
2.6
2.9
2.7


### 4.2. Effectiveness and Validity

There are a corresponding CNOP value and CNOP pattern for every initial optimization month, so there are 12 CNOP values and 12 CNOP patterns for all 12 different initial months. In this section, we compare our approach with the adjoint method from CNOP value and CNOP pattern aspects to verify the effectiveness and validity.


#### 4.2.1. CNOP Value

In this section, we set∇=ai∗10-1 when initial month is 1, 2, 5, 9, 11, and others; we set ∇=ai∗10-2 to get higher CNOP values. Figure 4 depicts CNOP values from our approach and adjoint method for different initial month; x-axis represents the initial optimization time (from January to December) and the y-axis represents the CNOP values. We can see that the variation trend of CNOP values for the two methods is almost the same. In detail, red and blue lines show upward trends from January to March; from March to September, they go down; and from September to December, they go up again.Figure 4
CNOP values calculated from our approach and adjoint method for different initial month. The blue line represents the values from adjoint method; the red one is the values from our approach. The green one is the difference between them.


#### 4.2.2. CNOP Pattern

In this section, the spatial pattern of the optimal precursor (SSTA-CNOP and THA-CNOP) of ENSO phenomenon and corresponding SSTA evolution are compared to assess the validity of our approach. It is unnecessary to show all 12 CNOP patterns. We choose the patterns of the two initial optimization months which has the biggest (March) and smallest (September) CNOP value, respectively.Figure5 shows the patterns of SSTA-CNOP, THA-CNOP, and corresponding SSTA evolutions after 9 months obtained from our approach and the adjoint method while the initial month is March, while Figure 6 shows the patterns of September. (a, b) is the pattern of SSTA-CNOP, (c, d) is the pattern of THA-CNOP, and (e, f) is the pattern of the SSTA evolution; (a, c, e) is the pattern from our approach and (b, d, f) is the pattern from the adjoint method. These two optimal precursors obtained by the two methods both can evolve into an EI Niño event.Figure 5
Patterns of SSTA-CNOP (a, b), THA-CNOP (c, d), and corresponding SSTA evolutions (e, f) obtained from our approach (a, c, e) and adjoint method (b, d, f) while the initial month is March.
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
As in Figure5, the patterns corresponding to the initial month September.
(a)
(b)
(c)
(d)
(e)
(f)Figures5(a), 5(b), 6(a), and 6(b), the patterns of SSTA-CNOPs show almost the same spatial structure. The SST of western Pacific is abnormally high around the equatorial Pacific, while eastern Pacific is opposite. It is just the precursor of the EI Niño event. The difference is red and blue area of (b) is larger and darker. From Figures 5(c), 5(d), 6(c), and 6(d), the patterns of THA-CNOPs also show almost the same feature. The color deepened along the entire equatorial Pacific. The difference is that red area of (d) is larger. From Figures 5(e), 5(f), 6(e), and 6(f), evolutions of SSTA still show quite similar spatial feature. (f) of SSTA evolution is positive while (e) is negative. The difference is red area of (f) is larger.In a word, the CNOP pattern from our approach is quite similar to those from the adjoint method but is a little weaker. The result is in accordance with the results in Section4.2.1, which shows that the CNOP values from our approach are a bit smaller than those from the adjoint method. In conclusion, our approach can obtain the valid optimal precursor for ENSO phenomenon.


### 4.3. Time Efficiency

In this section, we demonstrate the time consumption and speedup up to 240 CPU cores to verify the efficiency of our approach. In this work, the average value of running the same program ten times is set as the final time consumption and speedup is the radio of serial execution time over the parallel execution time. We employ 12 CPU cores as a unit because each node in the cluster has 12 CPU cores. There are 20 nodes, 240 CPU cores totally.In Figure7, we show the time consumption diagram corresponding to the adjoint method, our serial approach, and our parallel approach with 240 CPU cores. When the CPU cores is 240, the time consumed is 12.83 s, which is less than the time spent by adjoint method and the speedup reaches 85.18.Figure 7
The time consumption diagram for adjoint method (15,4 s), our serial approach (1107.13 s), and our parallel approach with 240 CPU cores (12.83 s).To show the effectiveness of the parallel strategies designed in Section3.2, we show the time consumption and speedup with the number of CPU cores increasing from 12 to 12∗20 in Figure 8. The blue line stands for the time consumption and red line stands for speedup. With the number of CPU cores increasing from 12 to 12∗20, the time assumption is falling and the speedup grows almost linearly. From the trend of decreasing for the time assumption, we can expect less time consumption if more CPU cores are provided. And the speedup also has the trend of continually increasing if more CPU cores are provided. Of course, there exist bottlenecks for both time consumption and speedup with the increasing of CPU cores; we cannot find the bottleneck owing to the lack of computing resources.Figure 8
Time consumption (blue line) and speedup (red line) with the increasing CPU cores.


## 4.1. Determination of∇

In Section3, we show the mathematical formula (see (9)) to calculate the gradient of objective function. In the equation, ∇ is a real constant which should approach 0 but never equals 0. For our case, the value of ∇ cannot be too small, because too small ∇ will lead no evolution for numerical model; that is, the limit value in (9) always equals or approaches 0; thus we cannot obtain the correct gradient direction. Therefore, what value ∇ should be in our case? We design the following two schemes to determine the value of ∇: (1) ∇ is constant value for every ai in (9); (2)  ∇ is changing with the value of ai in (9). For the CNOP calculated by different methods, the value of objective function (Jx0) is larger; then the CNOP is much better. So, in this paper, we will take the norm Jx0 to measure the magnitude of CNOP, and we take the value of Jx0 as the evaluating standard for CNOP value.For scheme (1), we conduct lots of experiments to solve CNOP, but it is found that the different initial month is corresponding to a different appropriate ∇ for CNOP. And it requires lots of experiments to determine the most appropriate ∇ for each initial month. Therefore, the conclusion is that the scheme (1) is not feasible in our case.For scheme (2), we let ∇=ai∗10-n(n=1,2,3). When n≥3, ∇ is too small, and CNOP value obtained is very small. When n<0, ∇ is too large to calculate the limit value. When n=1 and n=2, we compare the maximum value of objective function obtained by our approach with the results of adjoint method (shown in Figure 3). In Figure 3, the blue solid line represents the CNOP values calculated by adjoint method for different initial month, which is the baseline, while the red one (∇=ai∗10-1) and green one (∇=ai∗10-2) are the CNOP values using our approach.Figure 3
CNOP values calculated by our approach and adjoint method for different initial month. The blue line represents the values by adjoint method; the red one(∇=ai∗10-1) and green one (∇=ai∗10-2) are the values by our approach.In Figure3, the CNOP value obtained by our approach has similar tendency with the results of adjoint method and the CNOP value for every initial month is less than adjoint method, but the largest difference between them is 5.3 (Table 1), which is acceptable, so we can draw the conclusion that ∇=ai∗10-n(n=1,2) is appropriate for our approach.Table 1
The difference value between CNOP values calculated by adjoint method and our approach.
Initial month
1
2
3
4
5
6
7
8
9
10
11
12
a
i
∗
10
-
1
3.0
2.8
3.1
3.9
1.8
1.7
2.1
1.6
1.9
5.3
1.8
2.1
a
i
∗
10
-
2
3.1
3.1
2.5
3.6
3.1
1.3
1.3
1.1
3.3
2.6
2.9
2.7


## 4.2. Effectiveness and Validity

There are a corresponding CNOP value and CNOP pattern for every initial optimization month, so there are 12 CNOP values and 12 CNOP patterns for all 12 different initial months. In this section, we compare our approach with the adjoint method from CNOP value and CNOP pattern aspects to verify the effectiveness and validity.


### 4.2.1. CNOP Value

In this section, we set∇=ai∗10-1 when initial month is 1, 2, 5, 9, 11, and others; we set ∇=ai∗10-2 to get higher CNOP values. Figure 4 depicts CNOP values from our approach and adjoint method for different initial month; x-axis represents the initial optimization time (from January to December) and the y-axis represents the CNOP values. We can see that the variation trend of CNOP values for the two methods is almost the same. In detail, red and blue lines show upward trends from January to March; from March to September, they go down; and from September to December, they go up again.Figure 4
CNOP values calculated from our approach and adjoint method for different initial month. The blue line represents the values from adjoint method; the red one is the values from our approach. The green one is the difference between them.


### 4.2.2. CNOP Pattern

In this section, the spatial pattern of the optimal precursor (SSTA-CNOP and THA-CNOP) of ENSO phenomenon and corresponding SSTA evolution are compared to assess the validity of our approach. It is unnecessary to show all 12 CNOP patterns. We choose the patterns of the two initial optimization months which has the biggest (March) and smallest (September) CNOP value, respectively.Figure5 shows the patterns of SSTA-CNOP, THA-CNOP, and corresponding SSTA evolutions after 9 months obtained from our approach and the adjoint method while the initial month is March, while Figure 6 shows the patterns of September. (a, b) is the pattern of SSTA-CNOP, (c, d) is the pattern of THA-CNOP, and (e, f) is the pattern of the SSTA evolution; (a, c, e) is the pattern from our approach and (b, d, f) is the pattern from the adjoint method. These two optimal precursors obtained by the two methods both can evolve into an EI Niño event.Figure 5
Patterns of SSTA-CNOP (a, b), THA-CNOP (c, d), and corresponding SSTA evolutions (e, f) obtained from our approach (a, c, e) and adjoint method (b, d, f) while the initial month is March.
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
As in Figure5, the patterns corresponding to the initial month September.
(a)
(b)
(c)
(d)
(e)
(f)Figures5(a), 5(b), 6(a), and 6(b), the patterns of SSTA-CNOPs show almost the same spatial structure. The SST of western Pacific is abnormally high around the equatorial Pacific, while eastern Pacific is opposite. It is just the precursor of the EI Niño event. The difference is red and blue area of (b) is larger and darker. From Figures 5(c), 5(d), 6(c), and 6(d), the patterns of THA-CNOPs also show almost the same feature. The color deepened along the entire equatorial Pacific. The difference is that red area of (d) is larger. From Figures 5(e), 5(f), 6(e), and 6(f), evolutions of SSTA still show quite similar spatial feature. (f) of SSTA evolution is positive while (e) is negative. The difference is red area of (f) is larger.In a word, the CNOP pattern from our approach is quite similar to those from the adjoint method but is a little weaker. The result is in accordance with the results in Section4.2.1, which shows that the CNOP values from our approach are a bit smaller than those from the adjoint method. In conclusion, our approach can obtain the valid optimal precursor for ENSO phenomenon.


## 4.2.1. CNOP Value

In this section, we set∇=ai∗10-1 when initial month is 1, 2, 5, 9, 11, and others; we set ∇=ai∗10-2 to get higher CNOP values. Figure 4 depicts CNOP values from our approach and adjoint method for different initial month; x-axis represents the initial optimization time (from January to December) and the y-axis represents the CNOP values. We can see that the variation trend of CNOP values for the two methods is almost the same. In detail, red and blue lines show upward trends from January to March; from March to September, they go down; and from September to December, they go up again.Figure 4
CNOP values calculated from our approach and adjoint method for different initial month. The blue line represents the values from adjoint method; the red one is the values from our approach. The green one is the difference between them.


## 4.2.2. CNOP Pattern

In this section, the spatial pattern of the optimal precursor (SSTA-CNOP and THA-CNOP) of ENSO phenomenon and corresponding SSTA evolution are compared to assess the validity of our approach. It is unnecessary to show all 12 CNOP patterns. We choose the patterns of the two initial optimization months which has the biggest (March) and smallest (September) CNOP value, respectively.Figure5 shows the patterns of SSTA-CNOP, THA-CNOP, and corresponding SSTA evolutions after 9 months obtained from our approach and the adjoint method while the initial month is March, while Figure 6 shows the patterns of September. (a, b) is the pattern of SSTA-CNOP, (c, d) is the pattern of THA-CNOP, and (e, f) is the pattern of the SSTA evolution; (a, c, e) is the pattern from our approach and (b, d, f) is the pattern from the adjoint method. These two optimal precursors obtained by the two methods both can evolve into an EI Niño event.Figure 5
Patterns of SSTA-CNOP (a, b), THA-CNOP (c, d), and corresponding SSTA evolutions (e, f) obtained from our approach (a, c, e) and adjoint method (b, d, f) while the initial month is March.
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
As in Figure5, the patterns corresponding to the initial month September.
(a)
(b)
(c)
(d)
(e)
(f)Figures5(a), 5(b), 6(a), and 6(b), the patterns of SSTA-CNOPs show almost the same spatial structure. The SST of western Pacific is abnormally high around the equatorial Pacific, while eastern Pacific is opposite. It is just the precursor of the EI Niño event. The difference is red and blue area of (b) is larger and darker. From Figures 5(c), 5(d), 6(c), and 6(d), the patterns of THA-CNOPs also show almost the same feature. The color deepened along the entire equatorial Pacific. The difference is that red area of (d) is larger. From Figures 5(e), 5(f), 6(e), and 6(f), evolutions of SSTA still show quite similar spatial feature. (f) of SSTA evolution is positive while (e) is negative. The difference is red area of (f) is larger.In a word, the CNOP pattern from our approach is quite similar to those from the adjoint method but is a little weaker. The result is in accordance with the results in Section4.2.1, which shows that the CNOP values from our approach are a bit smaller than those from the adjoint method. In conclusion, our approach can obtain the valid optimal precursor for ENSO phenomenon.


## 4.3. Time Efficiency

In this section, we demonstrate the time consumption and speedup up to 240 CPU cores to verify the efficiency of our approach. In this work, the average value of running the same program ten times is set as the final time consumption and speedup is the radio of serial execution time over the parallel execution time. We employ 12 CPU cores as a unit because each node in the cluster has 12 CPU cores. There are 20 nodes, 240 CPU cores totally.In Figure7, we show the time consumption diagram corresponding to the adjoint method, our serial approach, and our parallel approach with 240 CPU cores. When the CPU cores is 240, the time consumed is 12.83 s, which is less than the time spent by adjoint method and the speedup reaches 85.18.Figure 7
The time consumption diagram for adjoint method (15,4 s), our serial approach (1107.13 s), and our parallel approach with 240 CPU cores (12.83 s).To show the effectiveness of the parallel strategies designed in Section3.2, we show the time consumption and speedup with the number of CPU cores increasing from 12 to 12∗20 in Figure 8. The blue line stands for the time consumption and red line stands for speedup. With the number of CPU cores increasing from 12 to 12∗20, the time assumption is falling and the speedup grows almost linearly. From the trend of decreasing for the time assumption, we can expect less time consumption if more CPU cores are provided. And the speedup also has the trend of continually increasing if more CPU cores are provided. Of course, there exist bottlenecks for both time consumption and speedup with the increasing of CPU cores; we cannot find the bottleneck owing to the lack of computing resources.Figure 8
Time consumption (blue line) and speedup (red line) with the increasing CPU cores.


## 5. Correctness and Physical Meaning of the CNOP

To demonstrate the correctness of the CNOP calculated by the proposed approach, we calculate the change rate of the energy norm increment(Et-E0)/Eo from the CNOP over the integrating months according to [36], that is, the net growth rate of the energy (Figure 9). The energy norm E is defined as T, where T is the sea surface temperature and T is the 2-norm of T. Figure 9 shows that the energy from CNOP is increasing nonlinearly over the integrating months, and the energy increases around 35 times when integrating 12 months. Therefore, the CNOP calculated can show the fast nonlinear growth, which illustrates the physical definition of the CNOP. Furthermore, the CNOP patterns obtained from the proposed approach (Figures 5(a), 5(c), 6(a), and 6(c)) are similar to those from the adjoint method (Figures 5(b), 5(d), 6(b), and 6(d)), which also illustrates the correctness of the CNOP.Figure 9
The net growth rate of the energy from CNOP in 12 months.In physics, the CNOP can represent the optimal precursor that will induce the occurrence of certain physical events. As we know, when the El Niño event occurs, the sea surface temperature will present anomalously warm in the eastern and central tropical Pacific Ocean area. And the spatial patterns (Figures5(e) and 6(e)) of the CNOP evolution we calculated exactly show the same abnormity in the eastern and central tropical Pacific Ocean region, which means that the CNOP is just the optimal precursor of the El Niño event.


## 6. Conclusions and Future Works

In this paper, we proposed an efficient approach based on the gradient definition to solve CNOP around the whole solution space and some parallel strategies were designed to improve gradient calculation efficiency. It is the first time to solve CNOP using gradient definition around the whole solution space.To verify the effectiveness and validity of our approach, we applied our approach to solve CNOP to study the optimal precursor of an ENSO event in ZC model. The experiment results indicate that our approach can obtain good results, and the time consumed is less than the adjoint method, and the time consumption still has the trend of continually decreasing when providing more CPU cores.The cruciality of the proposed approach is to calculate the gradient of the objective function using the gradient definition. Zebiak-Cane model is of medium complexity (103-dimensional) and the objective function is differentiable. The proposed approach is applied to the complex models, whether the objective function is differentiable and the time efficiency must be taken into account. For nondifferentiable models, the approximate gradients information for those nondifferentiable points can be obtained by the proposed approach. Inevitably, the solving efficiency will go down dramatically along with the rapidly increasing dimensions. However, kinds of methods can be adopted to improve the time efficiency, such as the parallel of the numerical models based on the CPU/GPU, reducing the dimension of the original solution space using appropriate dimension reduction methods. At present, we are concentrating on applying the proposed approach to MM5 and WRF models, which are more than 105-dimensional; related papers will soon be published.



---

*Source: 3208431-2017-11-14.xml*