A gradient-based optimization method with functional principal component analysis for efficient structural topology optimization

Structural topology optimization (STO) is usually treated as a constrained minimization problem, which is iteratively addressed by solving the equilibrium equations for the problem under consideration. To reduce the computational effort, several reduced basis approaches that solve the equilibrium equations in a reduced space have been proposed. In this work, we apply functional principal component analysis (FPCA) to generate the reduced basis, and we couple FPCA with a gradient-based optimization method for the first time in the literature. The proposed algorithm has been tested on a large STO problem with 4.8 million degrees of freedom. Results show that the proposed algorithm achieves significant computational time savings with negligible loss of accuracy. Indeed, the density maps obtained with the proposed algorithm capture the larger features of maps obtained without reduced basis, but in significantly lower computational times, and are associated with similar values of the minimized compliance.


Introduction
The need for optimized solutions in structural applications has increased over the years and has become nowadays fundamental because of the limited availability of commodities, environmental impacts, increasingly stringent industrial time-to-market requests, and emerging manufacturing processes as additive manufacturing. In this framework, structural topology optimization (STO) aims at obtaining optimized performance from a structure while satisfying some functional constraints, e.g., bounds on total mass or stresses. In particular, we focus on STO, whose goal is to produce optimized structures by determining an optimal mass or volume distribution in a given design domain. Unlike other alternatives as shape and size optimization, which deal with predefined configurations, STO can virtually produce any mass distribution within the design domain.
In the literature, STO is treated as a constrained minimization problem, for which several approaches and numerical schemes have been proposed (Sigmund and Maute 2013;van Dijk et al. 2013;Munk et al. 2015). Almost all the proposed approaches are iterative, which means that the optimized mass distribution is determined by repeatedly performing structural finite element analyses (FEA) that involve the solution of the equilibrium equations for the problem under consideration. On the one hand, iterative approaches are very effective because they allow to obtain optimized solutions in a variety of situations without requiring additional assumptions. On the other hand, however, their bottleneck is the computational effort, as they require solving FEA a large number of times. For example, in a minimum compliance problem, up to 97% of the total computational time could be spent for numerically solving the equilibrium equations (Alaimo et al. 2018;Petersson 1999).
From these considerations, it becomes evident the need for reducing the computational time to solve the equilibrium equations in STO, which is a key requirement for large problems and three-dimensional applications. Reduced basis approaches represent a valid solution to such an issue, because they allow to reduce the dimensionality of the structural problem, i.e., the number of equilibrium equations to solve and, therefore, the related computational effort.
However, reduced basis approaches have been exploited in a very limited number of recent works. Gogu (Gogu 2015) constructed reduced basis on-the-fly by using the previously generated solutions of the equilibrium equations, and the reduced basis was adaptively enriched based on the STO convergence. Ferro et al. (Ferro et al. 2019) applied the principal component analysis (PCA) to the density map, obtaining an efficient numerical scheme for STO. Xiao et al. (Xiao et al. 2020) proposed a proper orthogonal decomposition (POD) scheme, related to PCA, to construct a reduced basis for the displacement field solution of FEA, to be used during the optimization; subsequently, the reduced basis is adaptively updated on-the-fly according to an error metric. Less recently, Wang et al. (Wang et al. 2007) proposed a method based on the Krylov subspaces, while Amir et al. (Amir et al. 2009) proposed the construction of a reduced order model using the combined approximations method. Finally, Yoon (Yoon 2010) used eigenmodes along with Ritz vectors for the construction of a reduced model in dynamic STO.
In this work, we propose and test a reduced basis method that relies on functional PCA (FPCA). FPCA for FEA is quite new in STO since it has only been exploited for the variable thickness sheet design problem in combination with a simulated annealing optimization heuristic (Alaimo et al. 2018); moreover, it was initially developed for uncertainty quantification (Bianchini et al. 2015). We extend the use of FPCA in STO by applying it to the class of gradientbased optimization methods, which are the most commonly used. We demonstrate that significant savings in STO computational times can be achieved with negligible loss of accuracy. In this sense, the novelty of our work is not to introduce a new numerical technique but to appropriately combine two consolidated STO methods in order to save computational time. Such methods are the gradient-based optimization, (see for example Sigmund and Maute (2013)), and the reduced basis method (FPCA), able to reduce the dimensionality of the problem.
More specifically, we employ FPCA along with a gradient-based approach, considering the code proposed in Sigmund (2001) and his gradient-based approach as a starting point for incorporating FPCA reduced basis. At each STO iteration, we convert the FEA problem into a smaller one by projecting it onto a reduced space by means of a data-driven FPCA reduced basis. The reduced problem is solved and, if error measures associated with the solution are below the given thresholds, the reduced solution is accepted; otherwise, the original FEA problem is solved and the complete solution is used to update the FPCA basis.
We decided to embed FPCA rather than PCA in the optimization because of the results reported in Bianchini et al. (2015), where FPCA has been employed for uncertainty quantification purposes. Indeed, the authors showed that FPCA clearly outperforms PCA, i.e., FPCA always turned out to be faster than PCA and well captured the behavior of the solution with less elements in the basis than PCA. As the description of the FEM solution in a reduced space is the same both for uncertainty quantification and our optimization purposes, we can consider that these results are also valid in our setting.
The remainder of the paper is structured as follows: the proposed gradient-based optimization method with FPCA is presented in Section 2; the test case used to validate our algorithm is described in Section 3, while the results are reported in Section 4; finally, discussions and conclusion of the work are drawn in Section 5.

Gradient-based optimization with FPCA
In this section, we present the addressed STO problem (Section 2.1), the gradient-based optimization of Sigmund (2001) (Section 2.2), and the novel contribution included in the optimization, i.e., the FPCA to generate the reduced basis (Section 2.3). Finally, we detail the overall algorithm in Section 2.4. Without loss of generality, we present the approach for two-dimensional cases; from a conceptual point of view, the approach is the same also for threedimensional cases while, from an operational point of view, the extension is very straightforward.

Structural topology problem
We address a STO problem in which a given mass is distributed over the domain Ω in order to minimize the structural compliance: where ε is the strain tensor, σ = Dε the stress tensor, and D the fourth-order linear elastic tensor. We assume that D depends on the mass density distribution as D = ρ pD , where ρ ∈ [ρ min , 1] is the mass density, p ≥ 1 a penalty coefficient,D the fourth order elastic tensor of the isotropic homogeneous bulk material, and ρ min a positive value close to 0.
Finally, ε is defined in terms of the displacement field u as ε = 1 2 ∇u + ∇u T , where superscript T denotes transposition.
The problem is discretized following a classical FEA approach, using four-node elements with bilinear displacement functions. Accordingly, the compliance c is expressed in the form: where N el is the number of elements in the discretized domain,û i the vector of the nodal displacements of element i, f i the vector of the nodal forces of element i, and K i the stiffness matrix of element i, whileû and f are the assembled (over all elements) displacement vector and the applied force vector, respectively. Assuming a uniform mass density ρ i over each element i, K i depends on ρ i in the form K i = ρ p iK i , whereK i is the stiffness matrix of element i analogous toD, which corresponds to a homogeneous and isotropic material.
The minimization of c, expressed as in (1), is subject to mass (volume) constraint expressed in the form: where ρ * is the imposed mean density that uniquely defines the overall mass, and V i and V are the volume of element i and the overall domain volume, respectively. Finally, the equilibrium equations are expressed in compact form as: where K is the assembled stiffness matrix.

Gradient-based optimization
The gradient-based minimization is an iterative procedure in which, at each iteration, the new solution in terms of ρ i in all elements i is obtained following an optimality criterion, as proposed in Bendsøe and Sigmund (1995) and applied in Sigmund (2001). The optimality criterion is based on the sensitivity of c with respect to the densities ρ i in the form: However, there is no explicit formula that links the displacements to the density distribution. For this reason, the computation of derivatives (3) is performed by adding an adjoint term to c, in the form: whereũ is an arbitrary term, and the difference in parentheses is equal to 0 thanks to (2). The gradient of c is then expressed, for each component, as: which is simplified if the arbitrary termũ is chosen equal to the actual displacement fieldû, turning into: Considering the expression of the stiffness matrix K and that ρ i is constant over each element i, the gradient of the compliance can be computed element-wise and expressed as: Let us remark that (7) is valid whenû represents the exact solution of the equilibrium equations. This is of particular importance when considering a solutionû rec reconstructed from the reduced space, which can be slightly different from the exact solutionû; in this case, the complete expression (5) must be considered, or the problem must be adapted in order to use (7). In this work, as detailed in Sections 2.3 and 2.4, we consider the second alternative and we solve the problem with a modified force f curr in such a way that Kû = f curr .
With (7), it is possible to iteratively update ρ i to minimize the compliance c step by step. In particular, the new density distribution can be computed with the following heuristic approach proposed by Bendsøe and Sigmund (Bendsøe and Sigmund 1995): where ρ N i denotes the new density of element i, m is a positive coefficient that limits density variation, η = 1 2 is a damping coefficient, and α i is given by the optimality condition as: where κ is a Lagrangian multiplier found through a bisectioning algorithm (Sigmund 2001).

FPCA for reduced model
The idea behind the use of FCPA in STO is to solve the equilibrium (2) in a reduced space for as many gradientbased iterations as possible. Such a reduced space, whose dimension is given by a limited number of principal displacement components, is generated by applying FPCA to a set of previously obtained solutions (Alaimo et al. 2018;Bianchini et al. 2015).
To this aim, let us consider M displacement solutionŝ u(ρ m ) of the equilibrium equations (with m = 1, . . . , M) corresponding to several realizations ρ m of ρ, where ρ denotes the vector including all ρ i .
The subspace spanned by these M solutions represents an approximation of the entire solution space and the covariance operator is: Based on that, each basis function solves the following eigenproblem: Operatively, we start from M displacement solutionŝ , where L is the number of nodes in the mesh and 2L the number of degrees of freedom.
The buffer matrix U is then split into two parts: the average part U av and the residual part U cov = U − U av that stores the fluctuations of the solution with respect to the average.
The orthonormal basis is obtained with the functional principal components of . We use the method of snapshots as in (Alaimo et al. 2018;Bianchini et al. 2015), which is valid if 2L > M as in our case. Accordingly, the basis of rank K is obtained as follows: where φ k are the finite element shape functions, and M is a 2L × 2L matrix that can also be intended as a mass matrix. 2. Solve the eigenvalue problem where the eigenvalues λ k are in descending order along k. 3. Choose the value of K by retaining all eigenfunctions that correspond to an eigenvalue λ k > Λ λ 1 , where λ 1 is the maximum eigenvalue and Λ assumes very small values (Dulong et al. 2007).
The projection of (2) onto the reduced space of dimension K is carried out using the 2L × K matrix [ξ 1 , . . . , ξ K ], where Let now consider the problem Kû = f curr , where f curr is the modified force that allows to exploit (7). With the above projection, we may solve the following reduced problem, instead of the complete one, in the subspace spanned by the principal components: whereû rs denotes the solution of the problem in the reduced space andū = 1 M M m=1û i is the average of the solutions stored in buffer U.
After solving (11) in the reduced space, the displacement solution is reconstructed as: where subscript rec indicates that the solution is reconstructed starting fromû rs . Althoughû rec is computationally cheaper to obtain, it is only an approximation of the actualû. Therefore, we must decide whether to acceptû rec or not based on the approximation error and the quality of the solution.
The approximation error is evaluated by computing the errors e 1 and e 2 between the reconstructed force f rec = Kû rec and f curr and between f rec and f, respectively: The global mass matrix M ρ is introduced to weigh less or exclude the elements with lower density ρ i from the computation of the errors. It is computed as the assembly of elementary mass matrices, obtained for each element as Finally, the compliance c generally decreases with the iterations in gradient-based approaches (Sigmund 2001). Therefore, to avoid too large deviations from this decreasing behavior, a reduced solution is also evaluated in terms of the difference δ c between the compliance c obtained in the current iteration and that at the previous one.

Overall integrated algorithm
The algorithm reflects the iterative structure of the gradientbased optimization in Sigmund (2001), enriched by the proposed reduced basis approach based on FPCA. The basic structure of each iteration g (with g = 1, . . . , G) consists of solving equilibrium equations for the current density solution ρ g , and computing the new solution ρ g+1 according to (8).
The buffer U contains a number M of solutions between M min and M max . This basic structure of the gradientbased optimization, which consists of solving the complete problem (2), is repeated for the first M min iterations while considering the nominal force f of the problem and storing the displacement solutions in U. At iteration g = M min , the matrix U rb is also built from U.
Starting from iteration g = M min + 1, the problem is projected onto the reduced space exploiting U rb , and the equilibrium equations are solved in the reduced space according to (11)(12)(13)(14).
As mentioned, the reconstructed solutionû rec is accepted if e 1 < τ 1 , e 2 < τ 2 , and δ c < τ 3 . In this case, the reconstructed force f rec is also kept as the force f curr to solve the problem at the next iteration.
Ifû rec is not accepted, the complete problem (2) is solved while considering the nominal force f. At the same time, the force f curr to solve the problem at the next iteration is reset equal to f. Moreover, the new complete solution is included in U; if the number of already stored solutions is equal to M max , the oldest solution is also removed. The new U cov is then computed with the updated U.
At the end of each iteration g, the new density map ρ g+1 is generated according to the updating scheme (8), either considering the solutionû of the complete problem (in the first M min iterations or when the reduced solution is rejected) or the reconstructed solutionû rec (when the reduced solution is accepted starting from iteration g = M min + 1).
A detailed pseudo-code of the algorithm is reported in Algorithm 1.
The dynamic management of the buffer dimension M starting from M min allows to calculate, and possibly accept, reduced solutions even in the very first iterations without waiting for the first M max complete solutions. In fact, accepting reduced solutions even in the very first iterations is a key requirement for increasing the computational efficiency. At the same time, the upper limit M max allows to replace the oldest solution with the new one, thus giving preference to more recent displacement solutions that are associated with closer density distributions than the older ones.
Let us remark that, due to the approximation induced by the reconstructed displacement solutionû rec , the quantity Kû rec − f is not exactly null and, therefore, expression (7) is not exact in relation to problem (2). In the literature, such an expression has been modified by adding additional terms that compensate the error, as in (Gogu 2015), but this approach requires additional computational effort at each iteration. In this work, we deal with the issue by replacing the nominal force f with f curr . In this way, expression (7) becomes correct for the problem Kû = f curr without the need to calculate additional terms. Obviously, the problem with f curr is representative of (2) if the distance between f curr and f is limited, as verified by the conditions e 1 < τ 1 and e 2 < τ 2 together.

Experimental validation
We compare our algorithm with the case in which the reduced basis approach is not exploited and all gradientbased iterations are performed with the complete problem, as in (Sigmund 2001).
Analyses have been conducted on the well-known Messerschmitt-Bolkow-Blohm (MBB) beam problem. It consists of a simply-supported rectangular domain of size 2n x ×n y elements, with a vertical downward force F applied in the midpoint of the lower side (force module equal to 1). Thanks to the symmetry, only half of the domain is studied, as represented in Fig. 1. We consider n x = 2000 and n y = 1200 elements in the analyses, which define a computationally demanding problem having more than 4.8 million degrees of freedom.
Moreover, the following parameter values are assumed: -penalty p = 3; In addition to these experiments, different M min and M max values were tested for ρ * = 0.5, r = 5, and τ 3 = 2.
Each setting of ρ * and r is also solved with no reduced basis to get a reference solution for comparison. In particular, a number of iterations G = 100 are assumed to get the reference compliancec of the setting without reduced basis. This value of G has shown to be large enough to guarantee convergence; indeed, we also tested greater numbers of iterations, up to G = 1000, and we observed very limited variations ofc, lower than 2%, between iteration g = 100 and iteration g = 1000. Table 1 shows the referencec obtained with no reduced basis with G = 100 iterations for each setting of ρ * and r, together with the corresponding computational time. Moreover, it shows the CPU time to reach the referencec plus 5% and 10% for each τ 3 , and the percentage reduction with respect to the CPU time to reach the same compliance without reduced basis.

Results
The results show that the proposed algorithm is almost always able to find similar values of compliance c with respect to the referencec obtained with no reduced basis by the original method of Sigmund (Sigmund 2001), for at least one value of τ 3 . Moreover, a significant CPU time saving is always obtained. A reduction up to 52% and equal to at least 25% in most cases is observed when considering the referencec plus 5%, and the reduction is even higher when consideringc plus 10%, being up to 68%. These higher reductions show that the largest gain is obtained in the first iterations, where the gradient-based optimization gives higher c reduction per iteration. Thus, the algorithm proved to be effective in these iterations, and this confirms the effectiveness of the dynamic buffer management to start using reduced problems early.
It is worth remarking that the reduced solutions are rejected because the condition δ c < τ 3 is not respected, while the conditions e 1 < τ 1 and e 2 < τ 2 are always respected even with the low imposed thresholds. Moreover, a deeper analysis of Table 1 shows that the computational gain is not strictly correlated to the value of τ 3 , due to the fact that the time gain is affected by competing factors. On the one hand, smaller τ 3 values prevent instability of the whole algorithm and the need of restarting from higher values of c; on the other hand, this greater stability is paid for with a higher number of complete FEA iterations, which are clearly more expensive in terms of CPU time. Figure 2a, b and c show the evolution of the best known compliance c as a function of the CPU time and the number of iterations, respectively, for ρ * = 0.5 and r = 5, under the different values of τ 3 together with the case with no reduced basis. We may observe a faster reduction of the objective function c over time when adopting our algorithm with FPCA reduced basis, which corresponds to the shorter CPU times reported in Table 1 to getc plus 5% and 10%. In fact, although the convergence is slower in terms of iterations, the lower computational time required to solve the reduced problem makes the overall computational time shorter when exploiting FPCA. The red points in Fig. 2a and c highlight that the best known c always decreases in correspondence of a complete solution, and that a complete solution is required on average in only 26% or 27% of the iterations. Anyway, even when the objective function increases due to a bad representation of the solution with few bases (and the best known c remains constant over a certain time period), our algorithm is still able to redirect to acceptable values of c, which are even better than the case with no reduced basis at the same CPU time. Figures 3, 4, and 5 show the density maps obtained with the proposed algorithm whenc plus 5% is reached and the corresponding maps from the case with no reduced basis after G = 100 iterations. The maps are reported for each setting of ρ * and r, and for τ 3 = 2. Whenc plus 5% is not reached within G = 100 iterations by our algorithm for τ 3 = 2 (i.e., with {ρ * , r} equal to {0.3, 15} and {0.5, 10}), the map obtained from the algorithm at iteration g = 100 is reported. The maps from our algorithm show to capture the larger features of the solution. Moreover, although topologically different for some details, compared with the maps obtained with no reduced basis, our maps correspond to similar values of compliance c.
To clean the maps obtained by the algorithm from zones with intermediate density, a post-processing of the solution is considered, which consists of executing 5 additional iterations in which the complete problem is solved. Results show the effectiveness of this post-processing. In fact, it does not significantly change the structure of the maps, but at the same time, it effectively cleans intermediate density regions and eliminates local irregularities. Table 2 also shows that the modulus of the compliance variation induced by post-processing is always less than 1.20% and Table 1 Results from the experiments for each combination of ρ i , r, and τ 3 with M min = 2 and M max = 20: referencec with G = 100 iterations without reduced basis and corresponding CPU time; time to reach the referencec plus 5% and 10% with the proposed algorithm and percentage reduction with respect to the computational time to reach the same value without reduced basis "-" denotes thatc plus both 5% or 10% was not reached by our algorithm in G = 100 iterations on average equal to 0.78%. Therefore, we may argue that the presence of intermediate material in the maps obtained with the proposed algorithm involves only a small variation on the compliance which is, on average, less than 0.78%. At the same time, the maps cleaned by the post-processing are different from those in the case with no reduced basis. This could highlight the presence of local optimal solutions, or that the gap requires other iterations to be eliminated. We highlight that the post-processing has a limited impact on the computational time. In fact, assuming that all complete solutions require the same time, the percent gain obtained by the algorithm is reduced by only the 5%.

Impact of buffer dimension M
Additional experiments are run to analyze the impact of the M min and M max . Figure 2b and  the case with ρ * = 0.5, r = 5, and τ 3 = 2. As longer tracts with a constant best known c correspond to solutions with unstable c over the iterations, we may notice that the solution is far more stable when using a higher buffer size in terms of both M min and M max . In particular, a higher M min value avoids the destabilization of the solution in the very first iterations of the algorithm (red line in the figure), but at the price of a higher computational cost during the initial phase when the buffer starts to collect the solutions to build up the first reduced basis.

Discussions and conclusion
In this work, we propose an algorithm that makes the gradient-based optimization proposed by Sigmund and Maute (2013) for STO more computationally efficient.
In particular, we couple this gradient-based optimization with the FPCA method to generate reduced basis, onto which the equilibrium equations are projected at several optimization iterations. Unlike few similar works available in the literature, we exploit FPCA rather than PCA or POD because FPCA appears to be generally faster and well captures the behavior of the solution with less elements in the basis (Alaimo et al. 2018). Moreover, to face the gap between the solution reconstructed from the reduced space and the exact solution of the complete problem, we have decided to modify the problem by means of f curr , which allows us to avoid the additional calculations required to compute corrective terms. These two innovations have made our approach alternative with respect to the others proposed in the literature and suitable for practical applications, where computational time is the bottleneck.
The proposed algorithm has been tested on a large problem with more than 4.8 million degrees of freedom, much larger than the problems considered to test similar approaches. Results confirm that the reduced FEA problem Fig. 3 Density maps for ρ * = 0.3 and r = 5 (left column), r = 10 (central column) or r = 15 (right column). Upper row shows the maps from the proposed algorithm without postprocessing with τ 3 = 2 whenc plus 5% is reached; central row shows the post-processed maps after 5 supplementary complete solutions; lower row shows the maps obtained with no reduced basis at 100 iterations. a r = 5, algorithm. b r = 10, algorithm. c r = 15, algorithm. d r = 5, post-processing. e r = 10, post-processing. f r=15, post-processing. g r = 5, complete. h r = 10, complete. i r = 15, complete is exploited in the 73-74% of the optimization iterations ( Fig. 2a and c), guaranteeing limited errors below τ 1 , τ 2 , and τ 3 . This allows to reduce the STO computational times by the 50.3% on average (with minimum 25.5% and maximum 68.8%) to reachc plus 10% and by the 31.4% on average (with minimum 7.2% and maximum 52.1%) to reachc plus 5% in the experiments. . Upper row shows the maps from the proposed algorithm without postprocessing with τ 3 = 2 whenc plus 5% is reached; central row shows the post-processed maps after 5 supplementary complete solutions; lower row shows the maps obtained with no reduced basis at 100 iterations. a r = 5, algorithm. b r = 10, algorithm. c r = 15, algorithm. d r = 5, post-processing. e r = 10, post-processing. f r = 15, post-processing. g r = 5, complete. h r = 10, complete. i r = 15, complete Density maps for ρ * = 0.7 and r = 5 (left column), r = 10 (central column), or r = 15 (right column). Upper row shows the maps from the proposed algorithm without postprocessing with τ 3 = 2 whenc plus 5% is reached; central row shows the post-processed maps after 5 supplementary complete solutions; lower row shows the maps obtained with no reduced basis at 100 iterations. a r = 5,algorithm. b r = 10,algorithm. c r = 15,algorithm. d r = 5, post-processing. e r = 10, post-processing. f r = 15, post-processing. g r = 5, complete. h r = 10, complete. i r = 15, complete We remark that, while the uniqueness of the STO optimal solution has been proved for p = 1 (Petersson 1999), illposedness of the problem for higher p values (e.g., p = 3 used in this work) is well known (Bendsøe and Sigmund 1995;Tovar and Khandelwal 2010). Consequently, rather than simply being an approximation, the maps obtained using FPCA could represent optimal alternative solutions. Table 2 Compliance c for different values of r, ρ * , and τ 3 = 2 before post-processing (from the proposed algorithm) and after postprocessing (with 5 supplementary complete solutions), with relative variation after post-processing In fact, by modifying the problem by means of f curr , we do not follow the same sequence of solutions followed without reduced basis, allowing the procedure to move to other regions where different minima can be found. Obviously, we stopped the iterations when we reached a predefined gap with respect to the reference compliancec to pursue time savings, but the evolution could lead to a different solution even in the subsequent iterations. We believe this is the reason why our density maps with 5% or 10% gap are practically equivalent to those with referencec and those at convergence when the gradientbased STO is performed without reduced basis, although they are only similar from a topological point of view. But, by exploiting reduced basis through FPCA, these solutions are obtained in a drastically lower CPU time.
As mentioned, our strategy to include f curr has avoided additional calculations to compute corrective terms to adapt relation (6), as carried out in Gogu (2015). However, in this way, the reduced problems generate solutions with an increase of c in a few iterations. But, at the same time, the compliance c is reduced so much in the next iteration with a complete FEA solution to provide significant gain in the overall computing time. Our future work will consist in defining the best thresholds τ 1 , τ 2 , and τ 3 in order to optimize this trade-off. Furthermore, we will compare the proposed strategy with f curr with different alternatives in which corrective terms are added in (6).
From an application point of view, we will apply the proposed algorithm to several practical cases. For example, we will consider the computational gain when performing a pre-screening of different layouts, such as orientations or alternative features, while limiting a finer analysis without reduced basis to the selected layout only. At the same time, we can exploit the computational gain to explore different initializations of the optimization procedure in a reduced time.
Funding Open access funding provided by Universitàdegli Studi di Pavia within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.

Conflict of Interests
The authors have no financial or proprietary interests in any material discussed in this article.

Replication of results
The full code written in MATLAB language is available at the following link: www.mi.imati.cnr.it/ettore/FPCA.