A nonmonotone spectral projected gradient method for large-scale topology optimization problems

An efficient gradient-based method to solve the volume constrained topology optimization problems is presented. Each iterate of this algorithm is obtained by the projection of a Barzilai-Borwein step onto the feasible set consisting of box and one linear constraints (volume constraint). To ensure the global convergence, an adaptive nonmonotone line search is performed along the direction that is given by the current and projection point. The adaptive cyclic reuse of the Barzilai-Borwein step is applied as the initial stepsize. The minimum memory requirement, the guaranteed convergence property, and almost only one function and gradient evaluations per iteration make this new method very attractive within common alternative methods to solve large-scale optimal design problems. Efficiency and feasibility of the presented method are supported by numerical experiments.


Introduction
The goal of topology optimization is to find optimal material distribution in the given design domain subject to some constraints governed by certain physical properties and/or some other practical constraints during the design. In the past two decades, advances in the theory of homogenization, optimization, numerical analysis as well as newly developed engineering approaches make topology optimization techniques to become a standard tool of engineering design, in particular in the field of structural mechanics. For more detailed literature review, one may refer [1,3] and the references therein. In topology optimization, the design parameters are often material properties (e.g., conductivity or stiffness tensor) and the objective is often to minimize an integral functional defined on the spatial domain with state variables satisfying a partial differential equation (PDE) corresponding to certain physical law. In addition, some bound constraints on the design variables and usually a global material resource constraint are also often enforced during the design.
One typical large class of topology optimization problems have the structure that a nonlinear non-convex objective functional is minimized over a feasible region defined by a second order elliptic PDE together with bilateral bound and a single equality constraints. In this paper, we focus on solving this class of optimal design problems which have a wide range of applications in engineering design, such as compliance mechanics, fluid dynamics, heat transfer, functionally graded and composite materials [1,3].
In the nonlinear optimization literature, there are many well-known optimization methods for finding solutions of general optimization problem. Among them, the sequential quadratic programming (SQP) [14] method is widely used and generally considered to be an efficient method for smooth nonlinear optimization problems with constraints. Trust region methods [7] are another class of well-studied methods which have strong global and local convergence properties. More recently, interior-point methods receives much more attention because of their polynomial complexity. However, most of the above mentioned methods require the evaluation or a certain type of approximations of the hessian at each iteration. In addition, these methods also often require to solve a linear system of equations (usually indefinite, dense and is difficult to solve by iterative methods) to a certain level at each iteration. Therefore, when the problem size is very large, which often occurs after the discretization of the topology optimization problem, obtaining the Hessian information as well as solving very large linear system of equations at each iteration could be very expensive. Hence, although those second order methods enjoy fast local convergence properties, they are generally not efficient to solve very largescale problems, especially when a solution with very high accuracy is not strictly required.
On the other hand, some classes of optimization algorithms are developed within engineering community to solve large-scale engineering design problems. The first method of this kind is known as CONLIN or convex linearization [12]. More advanced version of CONLIN, called method of moving asymptotes (MMA), was introduced by Svanberg [19] . Within each iteration of this method the optimization problem is approximated by a convex separable sub-problem, for which efficient solvers are available. The globalization of the method is performed by either linesearch [21] procedures or conservative sub-iterations [20].
Despite the improvements of these methods for general structural optimization problems, a key issue to design an efficient optimization method is to exploit the specific structure of the desired problem. The optimality criteria (OC) method [3] is the first method developed particularly to solve the resource constrained topology optimization problems. Recently, OC becomes very popular and is the most widely used method in the engineering community for such problems. However, OC is not globally convergent and its application is limited to some special problems like thermal/structural compliance minimization [?, see:]ch.5]allaire2002soh. Moreover, the convergence rate of OC is not very promising.
By applying the recent techniques developed in the field of nonlinear optimization, in this paper we would like to design an optimization algorithm for solving very large-scale topology optimization problems. Each iterate of this algorithm is obtained by performing an adaptive nonmonotone line search along the line segment connected by the current point and the projection point of a cyclic Barzilai-Borwein (CBB) step onto the feasible set. Hence, the method can be called the projected cyclic Barzilai-Borwein (PCBB) method. Because of the special structure of the feasible set of the problem, which consists of box constants and a single linear constraint, it is possible to do the projection on the feasible set very efficiently with linear time complexity. The main attractive features of the presented algorithm are: producing strictly feasible iterations; using almost one objective function and gradient evaluations per iteration; O(n) memory consumption (6n working memory); easily to be implemented; only first order (gradient) information being required.

Barzilai-Borwein methods
Consider to solve the following finite-dimensional unconstrained optimization problem, where f : R n → R is continuously differentiable, and R n denotes the Euclidean space with dimension n. Suppose that x 0 is the starting point, x k is the current point, and g k is the gradient of f at x k , i.e., g k = ∇f (x k ). Then gradient methods generate the next iterative point by where the stepsize α k is computed by some line search techniques. Two classical ways of selecting initial stepsize in the line searches are given by the so called Steepest Descent (SD) and Minimal Gradient (MG) methods, which minimize f (x k −αg k ) and g(x k − αg k ) along the search direction g k , respectively: where · denotes the Euclidean norm of a vector. However, it is well-known that SD and MG methods can be very slow when the Hessian of f is singular or nearly singular at the local minimum. In this case the iterates could approach the minimum very slowly in a zigzag fashion [13]. The basic idea of Barzilai-Borwein (BB) [2] method is to use the matrix D(α k ) = 1 α k I, where I denotes the identity matrix, to approximate of the Hessian ∇ 2 f (x k ) by imposing a quasi-Newton condition on D(α k ): , and k 2. By straightforward calculation, BB stepsize obtained from (2.5) is By symmetry, another alternative BB stepsize could be computed by: which gives In contrast to the SD or MG methods in which the non-trivial and expensive optimization problem (2.3) or (2.4) need to be solved to obtain the initial stepsize, the BB stepsize is readily available during the iterations by formula (2.6) or (2.8). In practice, to keep the stability of the numerical procedure, it is often to project the BB stepsize onto a safeguard interval, that is to set The following Lemma shows the spectral property of Barzilai-Borwein methods and it is due to this property that these methods are usually called spectral gradient methods.
Lemma 2.1. The Barzilai-Borwein stepsize, α BB k , is the inverse of the Rayleigh quotient, related to vector s k−1 , of the averaged Hessian of the objective function between two consecutive iterations k − 1 and k.
Proof. By the Mean-Value Theorem and straightforward computations, one has Therefore, which complete the proof.
By Lemma 2.1, one has Λ min (α BB k ) −1 Λ max , where Λ max and Λ min are the maximum and minimum eigenvalues of the averaged Hessian matrix 1 0 ∇ 2 f (x k + ts k ) dt respectively. Therefore, α BB k I can be considered as an approximation of the inverse of the averaged Hessian of the objective function. This shows that the Barzilai-Borwein methods could incorporate certain useful second order information with extremely little additional expense of computational and memory cost compared with SD or MG methods. Due to its easy implementation, efficiency and low storage requirement, BB-type methods have been widely used in many applications. Exceptional good performances of BB-type methods have been observed for solving large-scale problems in particular when only approximate (not very accurate) solutions are desired, cf. [11].
It has been shown that if the exact steepest descent step (2.3) is reused in a cyclic fashion, the convergence speed of gradient methods can be greatly accelerated. However, it is often very expensive or impractical to find the exact stepsize along the steepest descent searching direction for large dimensional problems, unless the objective function is a quadratic function. Hence, for non-quadratic objective functions, it is unrealistic to apply some cyclic fashion of the steepest descent method. On the contrary, the BB stepsize (2.6) or (2.8) can still be easily calculated. Hence, analogous to the cyclic steepest descent method, the cyclic BB (CBB) has been introduced in [9] for general nonlinear optimization. The numerical results given in [9] show the CBB methods have excellent numerical performances compared with SD methods and even be competitive to some well-known nonlinear conjugate gradient methods. Given an integer m 1, which is the cycle length, CBB stepsize can be expressed as The R-linear convergence of CBB method for a strongly convex quadratic objective function has been proved in [8], while the local R-linear convergence for the CBB method at a local minimizer for general nonlinear objective function has been established in [9].

Projected Barzilai-Borwein methods
Consider the following constrained counterpart of problem (2.1) where D ⊆ R n is a closed non-empty convex set. Because of the constraints in (3.1), the iterations generated by (2.2) may lie outside of the feasible set D. Therefore, (2.2) need to be modified in order to maintain the feasibility of iterations. Gradient projection methods (see: [18]) keep the feasibility of iterates by frequently projecting trial steps generated from (2.2) onto the feasible set. Considering x as the trial step, the projection point y of x onto D, denoted by P D [x], can be computed by solving the following minimization problem Since the feasible set is convex, problem (3.2) always has an unique solution. However, in general (3.2) is still a convex constrained quadratic programming problem, which could be as difficult as the original problem. And for general convex constrained large-scale problems, this projection step at each iteration could be very time consuming and is normally the most expensive part of gradient projection methods. Hence, there is little interests on applying the gradient projection methods for large-scale problem unless the gradient projection step can be performed very efficiently. However, in some special cases when efficient algorithms for calculating the projection (3.2) exist, for example there are only box or single ball constraints, these gradient projection methods will be attractive. Fortunately, as we will see in the next section, the projection step can be performed very efficiently for the volume constrained topology optimization problems. This makes it possible for us to design gradient projection type algorithms for volume constrained topology optimization. Combining the Barzilai-Borwein stepsize rule and the gradient projection method, the projected Barzilai-Borwein (PBB) method was first introduced in [5]. In this method the iteration updating formula (2.2) is modified to is computed by connecting the current point to the projection point of the trial iterate (2.2) based on the BB stepsize, that is In [5], d α k was called spectral projected gradient. It is not difficult to show that d α k is a descent direction (see Lemma 3.1). This together with the convexity of the feasible set D would imply that for a sufficiently small β k , an iterate of (3.3) will reduce the objective function value while simultaneously preserve the feasibility of the iterates.
(ii) d α k (x) = 0 if and only if x is a stationary point for (3.1). Proof. see the Proposition 2.1 in [16].

Globalization by nonmonotone line search
Using descent directions and simply applying the SD, MG or BB stepsizes in gradient projection methods is generally not sufficient to ensure global convergence of the iterates starting from an arbitrary initial point. Hence, to deal with general nonlinear objective function, a globalization strategy is required to guarantee the global convergence of the algorithm.
Monotonic gradient-based methods usually generate a sequence of iterates for which a sufficient decrease in the objective function (or the related merit function) is enforced at every iteration. In many cases, the globalization strategy accepts the stepsize in the search direction, if it satisfy the well-known Wolfe or Armijo type conditions (cf. [?, ]ch. 3]nocedal2006no). This can be accomplished using either of monotonic line searches or trust region methods.
Since the search direction is parallel to the negative direction of the gradient (projected gradient) in BB (PBB) methods, the globalization of BB (PBB) methods by monotonic function value reduction often reduces them to the classic SD (projected SD) method. Hence, these monotonic methods will often distroy all of the advantages of BB (PBB) methods in contrast to SD (projected SD) method. To maintain the inherit spirits of BB-type methods, it is essential to accept the initial BB-type stepsize as frequently as possible while simultaneously ensure the global convergence. Hence, some nonmonotone line search techniques need to be developed to globalize the BB-type methods. The first nonmonotone line search technique was developed in [15] in which the main goal was to accept the full Newton step as much as possible. Combing the type of nonmonotone line search in [15], the first globalized version of BB nethod (GBB) was introduced in [17]. Following [17], the globalized PBB method (GPBB) was suggested in [5]. In these methods the following (weaker) objective function value decrease condition is enforced during each iteration where δ ∈ (0, 1), d k is the search direction and m k is a nonnegative nondecreasing integer, bounded by some fixed integer M . More precisely Based on the same motivations, the more efficient and adaptive nonmonotone line searches were particularly designed for BB-type methods in [16]. The globalized projected cyclic Barzilai-Borwein (PCBB) algorithm based on these new nonmonotone line search techniques can be described as follows: Algorithm 4.1. Parameters: • δ ∈ (0, 1), descent parameter used in Armijo line search.
The variable f r k in Algorithm 4.1 denotes the so called "reference" function value in nonmonotone line search. It can be seen that the traditional monotone line search simply corresponds to the choice of setting f r k = f (x k ) at each iteration. And the nonmonotone line search developed in [15] corresponds to the choice of setting f r k = f max k . In our present study, f r k is chosen based on Algorithm 4.2 adapted from [16]. Let f k denote f (x k ). In the algorithm 4.2, the integer a counts the number of consecutive iterations for which β k = 1 in Algorithm 4.1 is accepted and the Armijo line search in step 5 is skipped. The integer l counts the number of iterations since the function value is strictly decreased by an amount ∆ > 0.
R2. Set f R as follows in step 4 of Algorithm 4.1: ∆ for given decrease parameter ∆ > 0 and integer L > 0. Now lets give more details about computation ofᾱ k in step 1 of Algorithm 4.1. This parameter is computed based on the safeguarded CBB scheme as follows. Let j as an integer that counts the number of times in which the current BB step has been reused and let m as the CBB memory in (2.10), i.e., the maximum number of times the BB step will be reused. S0. If k = 0 chooseᾱ 0 ∈ [α min , α max ] and a parameter θ < 1 near 1; set j = 0 and flag=1. If k > 0 set flag = 0. S1. 0 < |d ki | <ᾱ k |g ki | for some i (component of vector) then set flag = 1. S2. If β k = 1 in Algorithm 4.1 then set j=j+1; Else (β k < 1) set flag =1. S3. If j m or flag=1 or s T k y k / s k y k θ then: S3.1. If s T k y k 0 then S3.1.1. If j > 1.5m then set t = min{ x k ∞ , 1}/ d 1 (x k ) ∞ , α k+1 = min{α max , max{α min , β k }} and j = 0; 2 Else setᾱ k+1 = min{α max , max{α min , α BB k }} and j = 0. In Algorithm 4.3, the former BB stepsize is reused for the current iterate unless one of the following conditions happens in which the new BB stepsize is computed (see S3.2) : (I) the previous BB stepsize is truncated by the projection step, i.e., when the trial point using BB stepsize lies outside of the feasible domain and the gradient projection was performed (see S1); (II) the previous BB stepsize is truncated by the line search step (see S2 where β k < 1); (III) the number of times the BB stepsize was reused reaches to its bound, i.e., j m (see S3); (IV) s T k y k / s k y k is close to 1 (see [?, ]section 4]dai2006cbb for details about the justification for this decision). The condition s T k y k < 0 (see S3.1) is equivalent to detection of the negative curvature in the searching direction. Assuming that the objective function can be well approximated by a quadratic function in the vicinity of the current iterate, a relatively large stepsize should be used in the next iteration (see S3.1.1) to reduce the function as much as possible once a negative curvature is detected. This strategy is similar to the original SPG algorithm (see [?, ]section 2]birgin2001ass). Now lets briefly review the convergence theory of Algorithm 4.1.
Theorem 4.4. Let L be the level set defined by We assume the following conditions hold: G1. f is bounded from below in L and d max = sup k d k < ∞. G2. IfL is the collection of x ∈ D whose distance to L is at most d max , then ∇f is Lipschitz continuous onL. Then either algorithm 4.1 with ǫ = 0 terminates in a finite number of iterations at a stationary point, or we have lim inf Proof. see the Theorem 2.2 in [16].
When f is a strongly convex function, (3.1) has a unique minimizer x * and the conclusion of the global convergence Theorem 4.4 can be strengthened as follows.

Projection onto the feasible set
The introduced algorithm in the previous section is general and can be applied to any type of convex constrained optimization problems which satisfy the requirements of Theorem 4.4. However, the performance of the method is significantly affected by the efficiency of the projection step. In this section, we explore the special structure of the feasible set in the volume constrained topology optimization and introduce a very efficient algorithm to do the projection step.
After discretization, the feasible set in the volume constrained topology optimization problems has the following form (see Figure 1) where a ∈ R n , a i ∈ R + , b ∈ R + , l, u ∈ R n , 0 < l i u i < ∞. Henceforth, in this section we denote by D the feasible domain defined in (5.1). Considering (3.2) together with (5.1), for a given trial vector x, P D [x] is the unique minimizer of the following box constrained Lagrangian: where λ ∈ R is a proper Lagrange multiplier related to the volume constraint and the box constrained set B is defined as B = {z ∈ R n : l z u}. For any fixed value of λ, (5.2) is a convex separable minimization problem with respect to z, which has the following explicit solution: where the max and min operators are understood as componentwise. Since solutions resulted from (5.3) satisfy the bound constraints, to solve the projection problem (3.2) the remaining task then turns out to find the Lagrange multiplier λ * such that (5.3) satisfies the volume constraint. Hence, an efficient algorithm is needed to find the (unique) root λ * of the following nonlinear non-smooth one-dimensional equation: Considering (5.3) together with (5.4), it can be seen that the graph of g(λ) has 2n breakpoints at: . . , n. Since l i u i and a i > 0 for all i, we have λ u i λ l i . So, each z i (λ) in (5.3) can be expressed in the following form: if λ λ l i . From (5.5), it is obvious that for each i, z i (λ) is a continuous piecewise linear and non-increasing function of λ (see Figure 2). Therefore, by a i > 0 for all i, g(λ) is a continuous piecewise linear and non-increasing function of λ. Then, with some straitforward calculations, we can see that the root λ * of g(λ) is unique and λ * ∈ [λ min , λ max ] with λ min 0 λ max , where λ min = min{λ u i } and λ max = max{λ l i }. So, starting from interval [λ min , λ max ] and using some classical one dimensional root finding method, it is possible to find λ * up to the machine precision with a priori known computational complexity. However, it is possible to use more advanced root founding methods to improve the performance of this step.
In our approach, the Brent's root finding method [6] is employed to solve (5.4). The Brent's method does not assume the function differentiability and is enable to manage the limited precision of the computed arithmetic very well. In the worst condition, the convergence of this method is never slower that of the bisection method. The Brent's method has been proved to be very efficient and robust in practice and it is currently accepted as a standard method for one dimensional root finding problem (cf. [?, ]chapter 9]press2007nre). The specific implementation details of Brent's method are available in the Numerical Recipe (see [?, ]chapter 9]press2007nre)..

Numerical solution of topology optimization problems
To show efficiency of the proposed method for the volume constrained topology optimization [4], in this section we solve the standard model problems presented in [10]. The description of the problem is as follows: considering two isotropic conducting materials, with thermal conductivities k α and k β , 0 < k α < k β , in a simply connected design domain Ω ⊂ R d (d = 2, 3), the goal is to mix these materials with a fixed ratio in Ω, such that the total temperature gradient in Ω is minimized under the thermal load f ∈ L 2 (Ω). More specifically, the problem can be formulated as 2 Ω ∇θ · ∇θ dx, subject to: where θ ∈ H 1 (Ω) is the state variable, p 1 is a penalization factor and w ∈ L 2 (Ω) is the control parameter (topology indicator field). By some standard derivations, the first order optimality condition of problem P can be written as the following: in Ω, in Ω, in Ω, where η ∈ H 1 0 (Ω) is the adjoint state, G is the L 2 gradient of the objective functional with respect to w and P D (u) denotes the L 2 projection of function u onto the admissible set D, By discretization of Ω into an n control volumes, we have the finite dimensional counterpart of problem (P). We also assume that the state variable and the design parameter are defined at the center of each control volume. Under these assumptions, the admissible design domain D forms a simplex in R n which is identical to the continuous knapsack constraints in (5.1). At each iterate of the control parameter w, we solve the associated state PDE. Then the discretized optimization problem would have the general format of problem (3.1), which has a nonlinear objective function with convex continuous knapsack constraints.
In our numerical experiments, the problem (P) was solved in two and three dimensions for Ω = [0, 1] 2 and Ω = [0, 1] 3 . The spatial domain is divided into 127 2 and 31 3 control volumes in two and three dimensions respectively. In all of experiments, we set k α = 1, f (x) = 1 and R = 0.4. And in these experiments two conductivity ratios 2 and 100 were tested, which is equivalent to setting k β = 2, 100 respectively. The penalization factor p is taken to be 1 and 10 for conductivity ratio 2 and 100 respectively. The governing PDE are solved by cell centered finite volume method using central difference scheme and the related system of linear equations are solved by a preconditioned conjugate gradient method with relative convergence threshold 10 −20 . The optimization is performed for 15 iterations in these numerical experiments. Notice that using finite volume method we do not observe any topological instability phenomena (the checkerboard pattern), which often occurs by using finite element method for topology optimization problems.
To evaluate the efficiency of the presented method, we compare our results with the results obtained by the method of moving asymptotic (MMA) [19]. The implementation of MMA available in SCPIP code [22] 1 is used in our experiments. All default parameters are used in SCPIP code, except the parameter for the constraint violation, for which the threshold 10 −7 is used in this study. That SCPIP code used in MMA algorithm has two globalization strategies. The first one is identical to that of the original MMA by Svanberg [20], while the second one is globalization by monotone line search method (see: [21]). The first strategy is employed in our experiments, since our results with second strategy was significantly worse in terms of the computational cost. Note that in our procedure the PDE constraint is solved upto the accuracy of the finite volume method we have applied for solving this PDE, and the constraints (5.1) in PCBB algorithm are satisfied upto the machine precision.
The results of our numerical experiments including the variation of the objective functional during the optimization process and the final resulted topology (w-field) are shown in figures 3, 4, 5 and 6. The plots in figures 3, 4 show the success of the presented method to solve these topology optimization problems. Roughly speaking, both methods behave very competitively in terms of computational cost and final results. The main differences in the results are related to the conductivity ratio 100, in which the PCBB performs superior and its final objective function values are considerably lower than that of MMA. The sign of the differences is clear in the figures 5 and 6 related to the final topologies. But for the conductivity ratio equal to 2, it seems MMA behaves slightly better than PCBB. However, the differences in terms of the objective function values as well as the final topologies are almost negligible .
In all our numerical experiments, the final total number of function and gradient evaluations was 15 which is equal to the number of optimization cycles. This result shows that both methods used only one function and gradient evaluations per iteration in practice. We believe that this is a key property for the success of MMA and makes the method well accepted in the engineering design community. In fact, MMA behaves very conservatively and uses small steps to proceed toward a local minimum. More clearly, it does not use (expensive) line search globalization strategy (unlike alternative methods), but uses reasonably small steps such that hopefully the merit function decreases sufficiently during each iteration. Of course, whenever the merit function increases (or sufficient decrease in merit function value violates), which rarely happens in practice, it performs sub-cycles to ensure the desired monotonic behavior. On the other hand, the nonmonotone PCBB methods enjoys such property using an alternative strategy. In practice, by applying the nonmonotonic line search, PCBB often uses only one function evaluation per optimization cycle. As our results clearly show, this property makes the nonmonotone PCBB method as a very competitive alternative to MMA for these class of problems.
It is important to note that in all of our experiments presented here, the projection step is used actively in PCBB method. Therefore, the cyclic reuse of step size never employed in our testing problems (cf. algorithm 4.3). Moreover, for the conductivity ratio 100, PCBB method explored directions of negative curvature 1 The SCPIP code (in Fortran) is freely available through personal request from its original author (Christian Zillober: christian.zillober(at)uni-wuerzburg.de).
after a few iterations, and so used large step-sizes which considerably accelerate the convergence. This property plays an important role for the superior results of PCBB compared with MMA in these cases (cf . table 1).
Besides the presented numerical results here, we have also applied the presented method successfully to many families of topology optimization problems with very satisfied results. But we omit the details of these experiments here due to the limited space for this paper.    ). α(2) and α(100) denote values of α k for the conductivity ratio 2 and 100 respectively.

Closing remarks
Recent studies on the spectral projected gradient methods show that these class of methods are very promising for solving the large-scale convex constrained optimization problems when the projection on the feasible set can be performed efficiently. In this paper, we presented a particular spectral projected gradient method called PCBB (Projected Cyclic Barzilai-Borwein) for solving volume constrained topology optimization problems. This method applies the cyclic Barzilai-Borwein stepsizes and uses the most recent adaptive nonmonotone line search techniques, which greatly improves the efficiency of the method as well as ensures its global convergence. By exploring the structure of the admissible set of the volume constrained topology optimization problem, the projection step can be performed very efficiently. In addition to high efficiency, our presented method also enjoys the following features: easy implementation, minimum memory requirement, no need for second order (Hessian) information, not sensitive to data noise, feasibility being strictly maintained during optimization process. All these features are essential for a successful numerical method to solve large-scale topology optimization problems.
Our numerical results indicate our presented methods are very promising and well-suited for the class of topology optimization problems considered in this paper. Comparing our results with those of MMA, a well accepted method in shape and topology optimization community, it seems the presented methods could be a very competitive alternative choice for solving the class of topology design problems. Finally, our results suggest that including spectral step size and a nonmonotone globalization strategy in MMA, the MMA algorithm could be further significantly improved. This would be a topic for our continuing research.