{'fraction_non_alphanumeric': 0.0620671834625323, 'fraction_numerical': 0.018604651162790697, 'mean_word_length': 4.091566899092225, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Normalized Wolfe-Powell-type Local Minimax Method for Finding Multiple Unstable Solutions of Nonlinear Elliptic PDEs

The local minimax method (LMM) proposed in [Y. Li and J. Zhou, SIAM J. Sci. Comput. 23(3), 840--865 (2001)] and [Y. Li and J. Zhou, SIAM J. Sci. Comput. 24(3), 865--885 (2002)] is an efficient method to solve nonlinear elliptic partial differential equations (PDEs) with certain variational structures for multiple solutions. The steepest descent direction and the Armijo-type step-size search rules are adopted in the above work and playing a significant role in the performance and convergence analysis of traditional LMMs. In this paper, a new algorithm framework of the LMMs is established based on general descent directions and two normalized (strong) Wolfe-Powell-type step-size search rules. The corresponding algorithm named as the normalized Wolfe-Powell-type LMM (NWP-LMM) are introduced with its feasibility and global convergence rigorously justified for general descent directions. As a special case, the global convergence of the NWP-LMM algorithm combined with the preconditioned steepest descent (PSD) directions is also verified. Consequently, it extends the framework of traditional LMMs. In addition, conjugate gradient-type (CG-type) descent directions are utilized to speed up the LMM algorithms. Finally, extensive numerical results for several semilinear elliptic PDEs are reported to profile their multiple unstable solutions and compared for different algorithms in the LMM's family to indicate the effectiveness and robustness of our algorithms. In practice, the NWP-LMM combined with the CG-type direction indeed performs much better among its LMM companions.


Introduction
Various nonlinear problems in physics, chemistry, biology and materials sciences can be reduced to consider multiple solutions of the Euler-Lagrange equation associated with a continuously Fréchet-differentiable nonlinear functional E defined on a real Hilbert space H, i.e., E (u) = 0, u ∈ H, (1.1) where E is the Fréchet-derivative of E. Actually, solutions of the Euler-Lagrange equation (1.1) are called critical points of the functional E. The most well-known candidates for critical points are local extrema to which classical variational and optimization methods have contributed a lot. Nowadays, with the development of new experimental techniques, it has been possible to observe local unstable equilibria or transient excited states in numerous physical/chemical/biological systems. Consequently, their theoretical and numerical studies have attracted increasing attentions [2,4,15,21,27,30]. However, these local unstable equilibria or transient excited states are related to critical points that are not local extrema, and then called saddle points. Virtually, in terms of the instability analysis, for a critical point u * with its second-order Fréchet-derivative E (u * ) existing, its instability can be depicted by its Morse index (MI) [2], which is defined as the maximal dimension of subspaces of H on which the linear operator E (u * ) is negative-definite. In fact, for a nondegenerate critical point u * , i.e., E (u * ) is invertible, if its MI = 0, it is a strict local minimizer and then a stable critical point, while if its MI > 0, then it is an unstable critical point. Generally speaking, the higher MI means the more instability.
Owing to the instability and multiplicity of saddle points, the design and analysis of numerical methods for grasping saddle points in a stable way is much more challenging than that for stable critical points. In recent years, various numerical methods were developed to capture saddle points in a stable way, such as the (climbing) string method [9,22], the gentlest ascent dynamics [10], the (shrinking) dimer method [14,29] and the search extension method [3,26]. Nevertheless, most existing methods in the literature mainly focus on finding saddle points with MI = 1.
In this paper, we are mainly interested in stable and efficient numerical computations for multiple saddle points with high MIs. Motivated by classical minimax theorems in the critical point theory (see, e.g., [21] and references therein) and numerical researches of Choi-McKenna [5], Ding-Costa-Chen [8] and Chen-Zhou-Ni [4], Li and Zhou proposed a local minimax method (LMM) for various high-MI saddle points based on a local minimax characterization of them [16]. Then in [27], Xie et al. modified the LMM with a significant relaxation on the domain of the local peak selection, which is a crucial notion for the LMM and will be illustrated in details later. According to [16,27], the LMM grasps a saddle point with MI = n (n ∈ N + ) by dealing with a two-level local minimax problem as where S H = {v ∈ H : v = 1} is the unit sphere with · the norm in H, L ⊂ H is a given (n − 1)-dimensional closed subspace usually constructed based on some known or previously found critical points, and [L, v] = {tv + w L : t ≥ 0, w L ∈ L} denotes a closed half subspace. Actually, the inner local maximization is equivalent to an unconstrained optimization problem in Euclidean spaces, which can be solved efficiently by standard optimization algorithms. The outer constrained local minimization, which is generally infinite-dimensional and much more challenging in the numerical computation, is the major concern of the LMM. For the sake of handling this task, the LMM adopts a normalized iterative scheme (NIS) with gradient descent directions d SD k = −g k , i.e., v k+1 = v k (α k ) = v k − α k g k v k − α k g k , w k+1 = p(v k+1 ), k = 0, 1, 2, . . . , (1.3) where α k > 0 is a step-size and g k = ∇E(w k ) ∈ H is the gradient of E at w k = p(v k ) with p(v k ) representing a local maximizer of E on [L, v k ], known as the so-called local peak selection of E w.r.t. L at v k . Actually, the local peak selection p(v k ) can be expressed as p(v k ) = t k v k + w L k for some t k ≥ 0 and w L k ∈ L [16]. One of fundamental problems for the NIS (1.3) is how to choose a suitable stepsize α k . In earliest implementations of the LMM, a normalized exact step-size search rule was employed and aimed to find the step-size α k > 0 such that the functional E(p(v k (α))) attains its minimum [16], i.e., E(p(v k (α k ))) = min α>0 E(p(v k (α))), k = 0, 1, 2, . . . . (1.4) Nevertheless, such a step-size search rule is very expensive in practical computations and it is even hard to establish the global convergence of the corresponding LMM algorithm. To compensate for this shortage, several normalized inexact step-size search rules with the advantages of low computational cost and easy implementation have been introduced in the literature to choose the step-size α k > 0 such that the decrease amount E(p(v k )) − E(p(v k (α k ))) > 0 is acceptable. Note that a widely applied normalized inexact step-size search rule in traditional LMMs is the normalized Armijo-type step-size search rule, which was first introduced in [17] and further simplified in [28] as the following form In fact, the factor 1/4 can be replaced by any constant σ ∈ (0, 1) [18]. Thanks to this step-size search rule, global convergence results for the normalized Armijo-type LMM (NA-LMM) algorithm were established in [17,31]. However, the decrease condition (1.5) is satisfied for all sufficiently small step-sizes (see Fig. 1), hence some artificial safeguards are needed to prevent step-sizes from being too small and the algorithm from interminable backtracking [17,31]. Recently, a normalized Goldstein-type step-size search rule was proposed in [18] to guarantee the sufficient decrease of the functional and prevent step-sizes from being too small simultaneously. The feasibility and global convergence analysis of the normalized Goldstein-type LMM (NG-LMM) were also provided in [18]. Actually, the normalized Goldsteintype step-size search rule in [18] makes progress with two inequalities, which can be formulated as with constants σ and δ satisfying 0 < σ < δ < 1. Unfortunately, as shown in Fig. 1, the normalized Goldstein-type step-size search rule may exclude the minimizer α * of E(p(v k (α))) outside the acceptable interval [ᾱ 1 ,ᾱ 2 ]. Thus, more effective and reasonable step-size search rules are still demanded to devote to the LMM for credibly capturing saddle points of the functional E. Furthermore, the convergence rate is another fundamental problem of the NIS (1.3) for solving the outer minimization in the two-level local optimization problem (1.2). As shown in the NIS (1.3), the steepest descent direction d SD k = −g k = Figure 1: Illustration of normalized Armijo-and Goldstein-type step-size search rules in traditional LMMs: the acceptable interval of the Armijo-type step-size is (0,ᾱ 2 ], while the acceptable interval of the Goldstein-type step-size is [ᾱ 1 ,ᾱ 2 ]. −∇E(p(v k )) was chosen as a descent direction in all existing LMM algorithms. As a result, they have certain limitations in terms of convergence rate. It is well known that, in line search algorithms for unconstrained optimizations in Euclidean spaces, there are many iterative directions that have better performance than the steepest descent direction, such as the conjugate gradient (CG) direction and the quasi-Newton direction, which usually have rapid convergence rates. Therefore, designing an improved iterative scheme of the form as (1.3) by replacing the steepest descent direction d SD k = −g k by other more efficient descent directions d k ∈ H for the outer minimization process will be beneficial to improve the numerical performance and convergence rate of the LMMs' family.
Inspired by classical (strong) Wolfe-Powell line search rules and descent methods in the optimization theory in Euclidean spaces [20,24,25] and the work of traditional LMMs [16,17,18,27], this paper is aimed to develop a normalized (strong) Wolfe-Powell-type LMM (NWP-LMM) to capture multiple unstable solutions of the Euler-Lagrange equation (1.1). The key point is how to design feasible Wolfe-Powell-type step-size search rules in a Hilbert space. By employing some curvature properties, two types of normalized Wolfe-Powell-type step-size search rules will be introduced for the LMM with general descent directions. Their mathematical justifications and global convergence will be established rigorously by making full use of these curvature properties, which quitely distinguish from those of the NA-LMM [17,27] and NG-LMM [18]. Finally, two types of descent directions, i.e., the preconditioned steepest descent (PSD) direction and the CG-type descent direction, will be proposed and compared to be implemented in our NWP-LMM algorithm for computing multiple unstable solutions of several semilinear elliptic partial differential equations (PDEs), such as the nonlinear Schrödinger equation (NLSE), Hénon equation and Chandrasekhar equation. Indeed, it will be seen that the CG-type descent direction can greatly speed up our NWP-LMM algorithm.
The rest of this paper is organized as follows. Firstly, some preliminaries for the LMM are provided in section 2. Then, in section 3, the normalized (strong) Wolfe-Powell-type step-size search rules are introduced together with the NWP-LMM given out. Their feasibility and some related properties are also discussed in this section. Global convergence of the NWP-LMM algorithm with general descent directions is verified rigorously in section 4. In addition, two different types of descent directions, including the PSD and CG-type descent directions, are proposed and analyzed in section 5 to feasibly implement our NWP-LMM algorithms. Furthermore, section 6 reports the detailed numerical results in 2D including the numerical comparison of different LMM algorithms for above mentioned semilinear elliptic PDEs to illustrate the effectiveness and robustness of our approach. Finally, some conclusions are drawn in section 7.

Preliminaries
For the convenience of discussions later, we introduce some notations and basic lemmas for the LMM in this section.
Let (·, ·) and · be respectively the inner product and norm in the Hilbert space H. Denote 2 H as the set of all subsets of the Hilbert space H, S H = {v ∈ H : v = 1} as the unit sphere in H and X ⊥ as the orthogonal complement to a subspace X ⊂ H. Suppose that L, serving as a so-called support space later, is a given closed finite-dimensional subspace in H. Define the half subspace [L, v] := {tv + w L : t ≥ 0, w L ∈ L} for any v ∈ S H . Throughout this paper, we assume that the functional E has a local minimizer at 0 ∈ H and focus on finding nontrivial saddle points of E. The local peak selection, a crucial notion, is defined as follows (cf. [16,18,27]).
For a given v ∈ S H , we say that E has a local peak selection w.r.t. L at v if there is a neighborhood N (v) of v and a mapping p : Let ·, · denote the duality pairing between H and its dual space H * . The above definition implies that each local peak selection p(v), with v ∈ S H , belongs to the well-known Nehari manifold N E := {u ∈ H\{0} : E (u), u = 0}, which contains all nontrivial critical points of the functional E. In fact, the following orthogonality holds obviously. 18,27]). Assume that E ∈ C 1 (H, R) has a local peak selection p w.r.t.
The following property follows from direct computation and is frequently utilized in the feasibility and convergence discussions for the LMM-type algorithms.
Imitating similar lines of the proof for Theorem 2.1 in [16], the following local minimax principle can be obtained and is referred to Theorem 4.1 in [18]. Theorem 2.4. If E ∈ C 1 (H, R) has a local peak selection w.r.t. L atv ∈ S H \L, denoted by p(v) = tvv + w L v , satisfying (i) p is continuous atv; (ii) tv > 0; and (iii) v is a local minimizer of E(p(v)) on S H , then p(v) / ∈ L is a critical point of E.
Since the local peak selection p(v) is a local maximizer of E on the half subspace [L,v], Theorem 2.4 characterizes a saddle point of E as a solution to the local minimax problem (1.2), or equivalently, the local minimization problem of E on the solution submanifold When a saddle point u * of E (known as an unstable critical point in H) can be characterized as a local solution to the two-level optimization problem (1.2) of the form Virtually, suitable descent algorithms can work for the minimization process in the optimization problem (1.2) or (2.1). As noted above, in traditional LMMs, the steepest descent direction serves as a search direction to numerically solve the local minimax problem (1.2). While, in this paper, the NWP-LMM will be constructed and analyzed for general descent directions. It is worthwhile to point out that the classical NA-LMM [17,27] and NG-LMM [18] and their mathematical justifications can also be extended to the case of general descent directions by similar process.

Normalized Wolfe-Powell-type LMM
In this section, a NWP-LMM framework for general descent directions will be proposed to capture multiple saddle points of the functional E via solving the optimization problem (1.2). In order to hit this goal, two types of normalized Wolfe-Powelltype step-size search rules with general descent directions will be introduced and analyzed by involving curvature conditions. We adopt the same notations as those in the section 2 unless specified and begin with some essential properties. Let Further, we have the following property.
Set p to be a local peak selection of E w.r.t. L at v ∈ S H \L. The following lemma is crucial in constructing the normalized Wolfe-Powell-type step-size search rule.
Proof. The mean value theorem states that, when s ∈ R is close to zero, where ξ = p(v(α)) + θ(p(v(α + s)) − p(v(α))) for some θ = θ(α) ∈ (0, 1). By the continuity of E , p and v(α), we have The local Lipschitz continuity of p and Lemma 3.1 imply that the right-hand side of In other words, E(p(v(α))) is differentiable w.r.t. α and (3.3) holds. Finally, it is clear that the right-hand side of (3.3) is continuous w.r.t. α. The proof is completed.
Remark 3.3. If E possesses a higher regularity, say, E ∈ C 1 (H, R) and E : H → H * is locally γ 1 -Hölder continuous for some γ 1 ∈ (0, 1], then the regularity assumption of p in Lemma 3.2 can be relaxed to that p is locally γ 2 -Hölder continuous around v(α) for some γ 2 ∈ (1/(γ 1 + 1), 1]. Actually, the key step in the proof is to justify that the right-hand side of (3.4) is o(s d ). The local Hölder continuity of E and p implies that there exist two constants ) and θ = θ(α) ∈ (0, 1) are the same as in (3.4). Consequently, the right-hand side of (3.4) is o(s d ).

Normalized Wolfe-Powell-type step-size search rules
Now, we are ready to propose normalized Wolfe-Powell-type step-size search rules for general descent directions by constructing curvature conditions to not only prevent step-sizes from being too small, but also avoid excluding the local minimizer of the function E(p(v(α))) w.r.t. α. In other words, the optimal step-size is always contained in the feasible step-size interval.
Remark 3.5. It is noted that the above definition is slightly different from the usual definition of descent direction in optimization theory. In fact, the descent direction d in Definition 3.4 is not only required to decrease the functional E at p(v) (i.e., E (p(v)), d < 0) but also to satisfy the orthogonality condition d ⊥ [L, v]. We make some comments on the reasons for introducing the latter condition as follows.
Actually, this orthogonality plays very important role in both algorithm implementation and theoretical analysis for the classical LMMs. Based on this observation, preserving this orthogonality to the general descent direction d for the outerlevel minimization is a preferred choice.
For two given constants σ 1 and σ 2 with 0 < σ 1 < σ 2 < 1, we say that the step-size α > 0 satisfies • the normalized strong Wolfe-Powell-type step-size search rule at v, if there hold The condition (3.5a) or (3.6a) is referred to as the sufficient decrease condition, while conditions (3.5b) and (3.6b) are referred to as curvature conditions. It is pointed out that, if the steepest descent direction is employed, (3.5a) (or (3.6a)) is equivalent to the normalized Armijo-type condition used in traditional LMMs [17,27,28].
The feasibility of normalized (strong) Wolfe-Powell-type step-size search rules above is provided as follows.
Remark 3.8. Fig. 2 provides a geometric interpretation of the feasibility of the normalized (strong) Wolfe-Powell-type step-size search rule.

Normalized Wolfe-Powell-type local minimax algorithm
Following the idea of traditional LMMs [16,17,27], the vital steps of the NWP-LMM algorithm are described as below.
Step 2. Compute a descent direction E (w k ), d k < 0, which will be discussed in section 5.
Step 3. If the stopping criterion E (w k ) H * < ε is satisfied (or more criteria are satisfied if necessary), then output u n = w k and stop; otherwise, go to Step 4.
Step 4. Set v k (α) = v k +αd k v k +αd k and find a step-size α k > 0 satisfying the normalized Wolfe-Powell-type step-size search rule, i.e., or the normalized strong Wolfe-Powell-type step-size search rule, i.e., Step 5.
Then, update k := k + 1 and go to Step 2.
We remark here that, similar to the classical (strong) Wolfe-Powell line search algorithm in optimization theory (see, e.g., [11,19,23]), one can employ an interpolation approach to efficiently find a step-size α k > 0 satisfying the normalized Wolfe-Powell-type step-size search rule (3.11) or the normalized strong Wolfe-Powelltype step-size search rule (3.12) in Step 4 of Algorithm 3.1. The implementation details are skipped here for brevity.

Global convergence
In this section, we establish the global convergence of Algorithm 3.1. In order to hit this goal, the following concept of compactness is needed. We simply utilize the same notations as those in Algorithm 3.1 throughout this section. It is pointed out that the following lemma gives a significant behavior of the sequence generated by Algorithm 3.1, which does not depend on the choice of descent directions and step-sizes in the algorithm. The proof is similar to that of Lemma 2.3 in [27] and omitted here for brevity.
Lemma 4.2 states that once an initial ascent direction v 0 = v L 0 + v ⊥ 0 ∈ S\L is chosen in Algorithm 3.1, the closed subset contains all possible vectors v k that Algorithm 3.1 may generate. Thus, the domain of a peak selection can be limited to be V 0 instead of S H . We remark here that, in general, the local peak selection p(v) defined on V 0 is no longer a homeomorphism. The following weak version related to the homeomorphism property of the local peak selection p(v) plays a significant role for establishing the global convergence. It can be verified by an analogous argument to that of Theorem 2.1 in [27] and only the continuity of the peak selection p(v) on V 0 is sufficient. Consequently, we skip the proof for simplicity. Lemma 4.3. Suppose E ∈ C 1 (H, R) and let p be a peak selection of E w.r.t. L satisfying (i) p is continuous on V 0 and (ii) t k ≥ δ for some δ > 0, ∀ k = 0, 1, . . .. If the sequence {w k } generated by Algorithm 3.1 contains a subsequence {w k i } converging to some u * ∈ H, then the corresponding subsequence {v k i } converges to some v * ∈ V 0 with u * = p(v * ).
Before proving the global convergence of Algorithm 3.1, we give some assumptions on the general descent direction, which are quite reasonable and hold for many descent directions, especially for the steepest descent direction used in traditional LMMs.
(A3) if {v k } contains a subsequence {v k i } converging to somev with E (p(v)) = 0, then the corresponding descent direction subsequence {d k i } converges.
Here, the assumption (A1) serves as a strong descent condition, and the assumption (A2) requires the length of the descent direction d k to be controlled by that of the gradient. In addition, the assumption (A3) admits a certain weak version of the continuous dependency of the descent direction d k on v k . The global convergence of Algorithm 3.1 is as follows.
Further, if the functional E satisfies the (PS) condition, then (d) {w k } contains a subsequence converging to a critical point u * / ∈ L. In addition, if u * is isolated, then w k → u * as k → ∞.
Proof. The decreasing condition (3.11a) (or (3.12a)), condition (ii) and assumption (A1) say that, for k = 0, 1, . . ., Therefore, the sequence {E(w k )} is monotonically non-increasing. In addition, since the condition (iii) guarantees that E(w k ) is bounded from below for all k = 0, 1, . . ., the sequence {E(w k )} converges to some E ∞ := inf k≥0 E(w k ). Then, adding up (4.2), we can arrive at Hence, the left-hand side of (4.3) converges to E ∞ − E(w 0 ) which is finite. This immediately leads to the conclusion (a). To prove the conclusion (b), letū ∈ H be an accumulation point of the sequence {w k }. Then, there exists a subsequence {w k i } converging toū as i → ∞. Recalling the weak version of the homeomorphism property in Lemma 4.3, it leads to that the corresponding subsequence {v k i } converges to somev ∈ V 0 satisfyingū = p(v) = tvv + w L v . Lemma 4.2 and the condition (ii) yield that, for some δ > 0, and thereforeū / ∈ L. The following is to verify thatū is a critical point by taking full advantages of the curvature condition (3.11b) (or (3.12b)), which states that By the contradiction argument, suppose thatū is not a critical point, i.e., E (ū) = 0. Since the functional E ∈ C 1 (H, R), one can obtain Therefore, (4.6) and (4.7) immediately indicate that holds by the continuity of E and p. Moreover, reviewing Lemma 2.3 and the condition (ii), for some δ > 0, we have t k i → tv ≥ δ > 0 and t k i +1 → tv ≥ δ > 0, as i → ∞. Since E (ū) = E (p(v)) = 0, it follows from the assumption (A1), (4.5) and (4.7) that E (ū),d ≤ −c 1 E (ū) 2 H * < 0 for some c 1 > 0. Thus, we have tv E (ū),d < 0 and (4.10) contradicts the fact 0 < σ 2 < 1. Consequently, the accumulation point u = p(v) is a critical point. The conclusion (b) is obtained.
Next, we prove the conclusion (c) by the contradiction argument. Suppose that Then, for k large enough, E (w k ) H * ≥ δ 1 /2 > 0. Thus, the conclusion (a) admits In addition, Lemma 3.1 and the assumption (A2) yield that, for some Combining (4.11) and (4.12) results in a Cauchy sequence in the closed subset V 0 . Immediately, the completeness of the closed subset V 0 implies that there existsv ∈ V 0 s.t. v k →v as k → ∞. Further, by the continuity of p and E , we have This is a contradiction to the conclusion (b). Thus, the conclusion (c) is proved. The rest is to prove the conclusion (d). According to conclusions (a) and (c), there exists a subsequence {w k i } s.t. E(w k i ) → E ∞ and E (w k i ) → 0 as i → ∞. By the (PS) condition, {w k i } possesses a subsequence, still denoted by {w k i }, that converges to a critical point u * ∈ H. In addition, according to the conclusion (b), u * / ∈ L holds. Finally, under the assumption that u * is isolated and following the analogous lines in the proof of Theorem 2.4 in [31] for the global convergence, the proof is completed.

Descent directions
In this section, we propose two specific types of descent directions for implementing Algorithm 3.1 in details. One is the PSD direction and the other is the CG-type descent direction. We use the same notations as those in Algorithm 3.1 for subsequent discussions in this section, unless specified.

Preconditioned steepest descent direction
The gradient of E at w k = p(v k ), denoted by g k = ∇E(w k ), is defined by From the Riesz representation theorem, the gradient g k ∈ H exists uniquely and g k = E (p(v k )) H * . Consider the following PSD direction d k = −T k g k , k = 0, 1, . . . , (5.2) where T k = T (v k ), called a preconditioner at v k , is a positive-definite and selfadjoint bounded linear operator on H with an invariant subspace The following theorem provides the global convergence result of Algorithm 3.1 with d k taken as the PSD direction (5.2).
In order to prove the conclusion, it suffices to verify assumptions (A1)-(A3) of Theorem 4.4. In fact, firstly, the assumption (T2) states that which is (A1) (with c 4 = c 1 ). Further, by the assumption (T1), Thus, (A2) is verified by taking c 3 = c 2 . Finally, (A3) directly follows from the assumption (T3) and the continuity of E and p.
Taking T k in Theorem 5.1 simply as the identity operator on H yields the following corollary, which draws the global convergence of Algorithm 3.1 with the standard steepest descent direction utilized at each iterative step.

Conjugate gradient-type direction
Similar to the construction of the nonlinear CG method in the optimization theory (see, e.g., [6,13]), we consider the following CG-type direction for Algorithm 3.1: Here, β k ∈ R is a parameter to be determined, g k = ∇E(w k ) is the gradient of E at w k defined in (5.1), and Π k is the orthogonal projection from H onto [L, v k ] ⊥ . Since g k ∈ [L, v k ] ⊥ from Lemma 2.2, we have d k ∈ [L, v k ] ⊥ for all k ≥ 0. In addition, according to the definition of Π k , it holds that (g k , Π k d k−1 ) = (g k , d k−1 ), k ≥ 1. Thus, the CG-type direction (5.5) satisfies Remark 5.4. We remark here that Π k d k−1 can be explicitly expressed as Actually, by applying the facts that where Id is the identity operator on H, we can conclude that Taking the inner product with v ⊥ k and noting that Consequently, the expression (5.7) is true. The following lemma shows that, if the exact step-size search rule is applied in the previous iteration, the CG-type direction (5.5) with an arbitrary parameter β k is a descent direction.
Due to the fact that the exact step-size search rule is quite expensive in practical computations, we focus on how to ensure that the CG-type direction d k defined in (5.5) is a descent direction when a suitable inexact step-size search rule is used.
Inspired by the well-known Fletcher-Reeves CG method [12] in the optimization theory, we set It can be verified that the CG-type direction (5.5) with β k = β FR-like k is a descent direction if the step-size α k in each iteration satisfies the normalized strong Wolfe-Powell-type step-size search rule (3.12) with σ 2 ∈ (0, 1/2), as stated in the following lemma. The proof follows the lines of the proof for Theorem 1 in [1]. Lemma 5.6. For k = 0, 1, . . ., if g k = 0, d k in Algorithm 3.1 is defined in (5.5) with β k = β FR-like k , and α k is determined by the normalized strong Wolfe-Powell-type step-size search rule (3.12) with σ 2 ∈ (0, 1/2), then Consequently, (g k , d k ) < 0, k = 0, 1, . . . .
Remark 5.7. Under assumptions in Lemma 5.6, we have from (5.9) that (g k , d k ) ≤ −c 1 g k 2 with c 1 = (1 − 2σ 2 )/(1 − σ 2 ) > 0, i.e., the CG-type direction (5.5) with β k = β FR-like k satisfies the assumption (A1) in Theorem 4.4. Moreover, under the same assumptions, these conclusions in Lemma 5.6 can be extended directly to any choice of β k satisfying |β k | ≤ β FR-like k . Remark 5.8. It is currently unclear whether assumptions (A2) and (A3) in Theorem 4.4 hold for the NWP-LMM with the CG-type descent direction (5.5). Thus, the global convergence of it has not been verified yet and will be our future work. Indeed, the NWP-LMM with the CG-type descent direction (5.5) is very efficient for finding multiple solutions of semilinear elliptic PDEs, compared to traditional LMMs and the NWP-LMM with the steepest descent direction, which will be shown in the section 6. Actually, several different constructions of CG-type descent directions can also be designed based on similar ideas of various CG methods in the optimization theory, which can be found, e.g., in [6,13].

Numerical examples
In this section, we apply our NWP-LMM to find multiple unstable solutions of the semilinear elliptic boundary value problem (BVP) where Ω is a bounded domain in R N with a Lipschitz boundary ∂Ω, a ∈ L ∞ (Ω) and a(x) ≥ 0 (x ∈ Ω), and f :Ω × R → R satisfies the following standard hypotheses [21]: (f1) f (x, ξ) is locally Lipschitz onΩ × R; (f3) there are constants µ > 2 and R > 0 s.t., for |ξ| ≥ R, Define H = H 1 0 (Ω) with the a-dependent inner product and norm as Since a(x) is nonnegative and uniformly bounded, the norm · a is equivalent to the usual norm in H 1 0 (Ω), i.e., u = Ω |∇u(x)| 2 dx 1/2 . The variational energy functional associated to the BVP (6.1) is given as It is well known that, under hypotheses (f1)-(f4), E ∈ C 1 (H, R) and satisfies the (PS) condition. Each critical point of E is a weak solution and also a classical solution to the BVP (6.1) [21]. In addition, u ≡ 0 is a local minimizer of E. Moreover, in any finite-dimensional subspace of H, E(u) → −∞ uniformly as u a → ∞. Hence, for any finite-dimensional closed subspace L, the peak mapping P of E w.r.t. L is nonempty. According to [16], we assume, in addition to hypotheses (f1)-(f4), that For L = {0}, as shown in [16], under hypotheses (f1)-(f5), E has only one local maximizer in any direction, i.e., E has a unique peak selection p w.r.t. L = {0}. For any finite-dimensional closed subspace L, the uniqueness of the peak selection of E w.r.t. L implies its continuity. As a result, the unique peak selection It is clear that hypotheses (f1)-(f6) hold for all functions f in Cases 1-3. In addition, our approach is efficient for different domains such as a L-shaped domain, a ball or other complex domains in high dimensions. In our numerical experiments, the initial ascent direction v 0 is taken as v 0 = v 0 / ṽ 0 a withṽ 0 the solution to the Poisson problem Particularly, it is necessary to explain more about how to numerically compute the gradient of the energy functional and a peak selection in numerical experiments. According to the definition of the gradient g k = ∇E(w k ) of E at w k ∈ H in (5.1) and the definition of the inner product (·, ·) a in (6.2), the gradient g k can be expressed as g k = w k − φ k with φ k determined by the linear elliptic BVP which can be solved efficiently by a standard numerical method, such as the finite element method (FEM) or the finite difference method. In our numerical code, assempde, a FEM-based subroutine provided by the MATLAB PDE Toolbox, is implemented to accomplish this task. The square domain in Cases 1-2 and the dumbbell-shaped domain in Case 3 are, respectively, discretized with 32768 and 15552 triangular elements. In addition, the computation of a peak selection p(v k ) of E at v k ∈ S H w.r.t. a given (n − 1)-dimensional closed subspace L (n ≥ 1) is equivalent to solving an unconstrained optimization problem in n-dimensional Euclidean space. To do this, a MATLAB subroutine fminunc with the termination tolerance on the first-order optimality Tol = 10 −8 is called in our numerical code. Further, for the convenience of numerical comparisons, we introduce the following notations for three different algorithms: • SD-Armijo: the traditional LMM algorithm by utilizing the steepest descent direction d k = −g k and the normalized Armijo-type step-size search rule [17,27]. At each iterative step of this algorithm, the step-size α k is chosen as where constants σ, λ are fixed as σ = λ = 0.1.

Numerical results for the nonlinear Schrödinger equation
In this subsection, we report numerical results of Case 1, i.e., the NLSE in the focusing regime on the square domain Ω = (−1, 1) 2 as x ∈ ∂Ω. (6.6) Taking ω = 8, limited by paper length, only ten different solutions labeled by u 1 , u 2 , . . . , u 10 are shown in Fig. 3 for their profiles and features. The corresponding information on the support space L, subdomains Ω 1 and Ω 2 used in (6.3), and energy values of these solutions is listed in Table 1. In addition, numerical comparisons of the SD-StrongWolfe, CG-StrongWolfe and SD-Armijo in terms of CPU times for computing these solutions of the NLSE are presented in Fig. 4. From Figs. 3-4, Table 1 and additional results not shown here, we observe that three LMM algorithms considered are effective for finding multiple solutions of the NLSE in the focusing regime and the CG-StrongWolfe shows the best performance among its LMM companions. In addition, u 1 (Fig. 3(a)) is the only positive solution with the lowest energy value, i.e., it is the ground state solution.
Fix = 6, due to the space limitation, only twelve solutions we obtained, labeled by u 1 , u 2 , . . . , u 12 , are displayed in Fig. 5 for their profiles and features. The corresponding information on the support space L, subdomains Ω 1 and Ω 2 used in (6.3),  and energy values of these solutions is listed in Table 2. In addition, numerical comparisons of the SD-StrongWolfe, CG-StrongWolfe and SD-Armijo in terms of CPU times for computing these solutions of the Hénon equation are provided in Fig. 6. From Figs. 5-6, Table 2 and additional results not shown here, we observe that three LMM algorithms considered can effectively find multiple solutions of the Hénon equation. As expected, the CG-StrongWolfe also shows the best performance among its LMM companions. In addition, positive solutions are not unique in the case of = 6.

Numerical results for the Chandrasekhar equation
Here, we report numerical results of Case 3, i.e., the Chandrasekhar equation as −∆u(x) = (u 2 (x) + 2u(x)) 3/2 , x ∈ Ω, x ∈ ∂Ω. (6.8)  Fig. 5. Figure 5: Profiles of twelve different solutions of the Hénon equation: (a) the singlepeak positive (ground state) solution u 1 concentrated mainly on the corner; (b) a two-peak positive solution u 2 concentrated mainly on two adjacent corners; (c) a two-peak positive solution u 3 concentrated mainly on two diagonal corners; (d) a two-peak sign-changing solution u 4 concentrated mainly on two adjacent corners; (e) a two-peak sign-changing solution u 5 concentrated mainly on two diagonal corners; (f) a three-peak positive solution u 6 ; (g)-(h) two three-peak sign-changing solutions u 7 and u 8 ; (i) a four-peak positive solution u 9 ; (j)-(l) three four-peak sign-changing solutions u 10 , u 11 and u 12 .
We consider a dumbbell-shaped domain Ω as depicted in Fig. 7(a). It contains a smaller disk centered at (−1, 0) with radius 0.5 and a larger disk centered at (2, 0) with radius 1. A corridor of width 0.4, symmetric respect to the x 1 -axis, is constructed to link the two disks. In this case, we focus on finding multiple positive solutions. Limited by the length of the paper, we only present seven different positive   Fig. 7. The corresponding information on the support space L, subdomains Ω 1 and Ω 2 used in (6.3), and energy values of these solutions is listed in Table 3. In addition, numerical comparisons of the SD-StrongWolfe, CG-StrongWolfe and SD-Armijo in terms of CPU times for computing these solutions of the Chandrasekhar equation are provided in Fig. 8. From Figs. 7-8, Table 3 and additional results not shown here, we observe that three LMM algorithms considered can effectively find multiple positive solutions of the Chandrasekhar equation. Again, the CG-StrongWolfe as always shows the best performance in this case.
Above all, numerical experiments in this section indicate that the CG-type direction indeed speed up the LMM greatly.

Conclusions
In this paper, we introduced normalized Wolfe-Powell-type and strong Wolfe-Powelltype step-size search rules for general descent directions, and proposed a normalized Wolfe-Powell-type local minimax method (NWP-LMM) for finding multiple unstable solutions of semilinear elliptic problems. Under certain conditions on the local peak selection and general descent directions, the feasibility and global convergence of the NWP-LMM were rigorously verified. In addition, two feasible types of descent directions, i.e., preconditioned steepest descent directions and the conjugate gradient-type direction, were proposed and discussed. The global convergence of the NWP-LMM combined with the preconditioned steepest descent directions was also provided. Extensive numerical results for several semilinear elliptic equations, including the nonlinear Schrödinger equation, Hénon equation and Chandrasekhar equation in 2D, were reported with their multiple solutions displayed to illustrate the effectiveness and robustness of our approach. The superior numerical performance of the NWP-LMM combined with the conjugate gradient-type direction was observed in extensive numerical experiments, while the rigorous verification for its global convergence is ongoing. Furthermore, designing more efficient preconditioned steepest descent or preconditioned conjugate gradient-type directions within the framework of LMM by constructing appropriate preconditioners to further improve the efficiency of computing multiple solutions will be our future work.