Tiago Carneiro
email: tiago.carneiropessoa@uni.lu

Nouredine Melab
email: nouredine.melab@univ-lille.fr

Akihiro Hayashi
email: ahayashi@gatech.edu

Vivek Sarkar
email: vsarkar@gatech.edu

Towards Chapel-based Exascale Tree Search Algorithms: dealing with multiple GPU accelerators

Keywords: Backtracking, GPU, PGAS, Chapel

Tree-based search algorithms applied to combinatorial optimization problems are highly irregular and timeconsuming when it comes to solving big instances. Solving such instances efficiently requires the use of massively parallel distributed-memory supercomputers. According to recent Top500 trends, the degree of parallelism in these supercomputers continues to increase in size and complexity, with millions of heterogeneous (mainly CPU-GPU) cores. Harnessing this scale of computing resources raises at least three challenging issues which are described and addressed in this paper. Indeed, as a step towards exascale computing, we revisit the design and implementation of tree search algorithms dealing with multiple GPUs, in addition to scalability and productivity-awareness using Chapel. The proposed algorithm exploits Chapel's distributed iterators by combining a partial search strategy with pre-compiled CUDA kernels for more efficient exploitation of the intra-node parallelism. Extensive experimentation on big N-Queens problem instances using 24 GPUs shows that up to 90% of the linear speedup can be achieved.

I. INTRODUCTION

Tree search algorithms are strategies that implicitly enumerate a solution space, dynamically building a tree. Algorithms from this class, such as backtracking and branch-and-bound (B&B), are popular to efficiently solve permutation combinatorial optimization problems (PCOPs), especially when combined with parallel/distributed computing. As the decision version of PCOP is often NP-complete, the size of problems that can be solved to optimality is limited, even if largescale distributed computing is employed [START_REF] Crainic | Parallel branch-and-bound algorithms[END_REF]. In this sense, the use of GPUs became crucial because it enables solving to the optimality instances having prohibitive execution times on CPUs [START_REF] Mezmaz | A grid-enabled branch and bound algorithm for solving challenging combinatorial optimization problems[END_REF], [START_REF] Gmys | A computationally efficient branch-and-bound algorithm for the permutation flow-shop scheduling problem[END_REF]. Also in the context of large-scale distributed computing, the use of accelerators such as GPUs and FPGAs plays a special role, as the energy-efficiency of such devices helps to break the power barrier towards exascale [START_REF] Heldens | The landscape of exascale research: A datadriven literature analysis[END_REF].

It is important to point out that such large-scale heterogeneous systems are going to be complex to program, and efforts towards productivity are crucial for better exploiting the future generation of supercomputers [START_REF] Fiore | On the road to exascale: Advances in high performance computing and simulations-an overview and editorial[END_REF]. The DARPA High Productivity Computing System Program (HPCS) represents an effort for creating high-productivity languages for the next generation of supercomputers [START_REF] Lusk | Languages for high-productivity computing: the darpa hpcs language project[END_REF]. Among the HPCS highproductivity languages, Chapel stands out, as it is competitive to both C-OpenMP and MPI+X in terms of performance and scalability [START_REF] Chamberlain | Chapel comes of age: Making scalable programming productive[END_REF]. Existing work on distributed tree search [START_REF] Carneiro | Towards ultrascale branch-and-bound using a high-productivity language[END_REF] show that it is possible to achieve parallel efficiency and performance, but also high productivity, by using the Chapel language. However, due to Chapel's lack of official GPU support, related works do not employ GPUs, which is critical in the exact optimization field [START_REF] Gmys | A computationally efficient branch-and-bound algorithm for the permutation flow-shop scheduling problem[END_REF].

The objective of the present research is to deal with multiple GPU accelerators as a major step towards exascale using a high-productivity language. In this sense, we revisit the design and implementation of tree search algorithms taking into account GPU heterogeneity, in addition to scalability and productivity-awareness using Chapel. The proposed algorithm is a distributed GPU-based tree search algorithm for solving PCOPs. Chapel's lack of official GPU support is overcome by combining the language's features with pre-compiled CUDA kernels. Chapel's distributed iterators are used for work distribution and load balancing between computer nodes. For more efficient exploitation of the intra-node parallelism, each subproblem received by a computer node is processed locally on CPU and then offloaded to the local GPUs.

Extensive experimentation on big N-Queens problem instances using up to 24 GPUs shows that 90% of the linear speedup can be achieved compared to an optimized CUDA-C baseline. We also compare the proposed algorithm to another implementation using an existing approach from the literature that enables GPU support in Chapel [START_REF] Hayashi | GPUIterator: Bridging the gap between Chapel and GPU platforms[END_REF]. Results show that the proposed algorithm is up to 1.77× faster and almost 2× more efficient than its counterpart, mainly due to its lack of distributed load balancing. It is worth to mention that we use the N-Queens problem as a proof-of-concept that motivates further improvements in solving related combinatorial optimization problems.

The remainder of this document is structured as follows. Section II brings background information and the related work. Section III details the distributed GPU-based tree search algorithm in Chapel. In turn, Section IV presents the performance evaluation of the proposed algorithm. Next, Section V brings a discussion of the reported experimental results. Finally, conclusions are outlined in Section VI.

II. BACKGROUND AND RELATED WORK A. Tree-based Search Algorithms

Tree search algorithms are strategies that implicitly enumerate a solution space, dynamically building a tree [START_REF] Grama | A survey of parallel search algorithms for discrete optimization problems[END_REF]. Algorithms that belong to this class start with an initial (root) node, which represents the initial state of the problem to be solved. Nodes are branched during the search process, generating children nodes more constrained than their father node. The generated nodes are evaluated, and then, the valid and feasible ones are stored in a pool-like data structure called Active Set. The search generates and evaluates nodes until the data structure is empty or another termination criterion is satisfied.

During the search, if an undesirable state is reached, the algorithm discards this node and then chooses an unexplored (frontier) node in the active set. This action prunes some regions of the solution space, preventing the algorithm from unnecessary computations. However, the pruning of subproblems makes the shape of the tree irregular, which might result in severe load imbalance when parallel computing is used. In this sense, load balancing schemes are crucial for achieving parallel efficiency in tree-based search algorithms.

B. The Chapel Language

Chapel is an open-source parallel programming language designed to improve productivity in high-performance computing. In Chapel, the program is started with a single task, and parallelism is added through data or task-parallel features [START_REF] Chamberlain | Chapel comes of age: Making scalable programming productive[END_REF]. Furthermore, as Chapel belongs to the partitioned global address space languages (PGAS), the application has a global memory addressing space, and each segment of this space is assigned to a different locale [START_REF] Almasi | PGAS (partitioned global address space) languages[END_REF]. In Chapel, the term locale usually refers to a symmetric multiprocessing computer in a parallel system. Due to the PGAS, a task can refer to any variable lexically visible, whether this variable is placed on the same locale on which the task is running, or on the memory of another one. Moreover, indexes of data structures can be globally expressed, even in the case where the implementation of such data structures distributes the indexes across several locales.

Although Chapel is a parallel programming language, there is no official support for GPUs. There are solutions that translate the body of a forall loop into GPU code [START_REF]Performance portability with the chapel language[END_REF], [START_REF]GPGPU support in chapel with the radeon open compute platform (extended abstract)[END_REF]. However, these solutions may not always deliver the best possible performance on GPU, particularly when advanced GPU optimization techniques such as warp-shuffles and shared memory utilization are required.

C. Chapel's Parallel Iterators

Iterators in Chapel are similar to procedures that can be used to isolate iterations from the loop body. Each value yielded by the iterator corresponds to an iteration of the loop. Chapel provides different parallel and distributed iterators that implement load balancing between computer cores and locales. Related work on distributed tree search shows that the iterators provided by Chapel are the key feature to achieve a trade-off between productivity and parallel efficiency [START_REF] Carneiro | Towards ultrascale branch-and-bound using a high-productivity language[END_REF]. The complex communication pattern of a distributed master-worker algorithm is encapsulated by the iterators. As a consequence, there is no need for explicitly dealing with work distribution, load balancing, termination criteria, or metrics reduction. Due to Chapel's lack of official GPU support, related work do not deal with the CPU-GPU heterogeneity, which is critical in the exact optimization field.

As Chapel does not provide any iterator dedicated to GPUs, Hayashi et al. [START_REF] Hayashi | GPUIterator: Bridging the gap between Chapel and GPU platforms[END_REF] propose the GPUIterator module that allows Chapel programmers to perform forall loops on both CPU and multi-GPU with minor additions to the code. The GPUIterator automatically divides a given iteration space into CPU and GPU portions by looking at a userspecified parameter, namely CPUPercent. A forall loop that is responsible for the CPU portion is executed on CPUs. On the other hand, the user needs to prepare a callback function that is expected to invoke a GPU program that is responsible for the GPU portion. It is worth noting that the GPUiterator also facilitates distributed execution. Also, in the current implementation, the GPUIterator does not dynamically update the CPU-GPU percentage to mitigate loadimbalance during the execution, neither provides distributed load balancing.

III. THE PROPOSED ALGORITHM

This section presents a productivity-aware distributed CPU-GPU backtracking to solve permutation combinatorial problems. The algorithm follows the master-worker model and is implemented for enumerating all complete and feasible solutions to the N-Queens problem.

Backtracking is a fundamental tree-based search that dynamically enumerates a solution space in a depth-first fashion. Due to its low memory requirements and its ability to quickly find new solutions, depth-first search (DFS) is often preferred as a search strategy for branch-and-bound (B&B) search algorithms [START_REF] Zhang | Branch-and-bound search algorithms and their computational complexity[END_REF]. In turn, the N-Queens problem consists of placing N non-attacking queens on a N × N chessboard, and it is often used as a benchmark for novel tree-based search algorithms [START_REF] Bell | A survey of known results and research areas for n-queens[END_REF]. N-Queens is easily modeled as a permutation problem: position r of a permutation of size N designates the column in which a queen is placed in row r.

The concepts herein presented are similar to any PCOP and can be adapted for solving other problems with straightforward modifications [START_REF] Gmys | IVM-based parallel branch-and-bound using hierarchical work stealing on multi-GPU systems[END_REF], [START_REF] Carneiro Pessoa | GPU-accelerated backtracking using CUDA dynamic parallelism[END_REF]. In this paper, we use the N-Queens problem as a proof-of-concept that motivates further improvements in solving related combinatorial optimization problems.

A. Initial Premises

One can see in Section II-B that a couple of approaches try to mitigate the gap between Chapel and GPUs. Those solutions are not feasible for the implementation of distributed treebased search algorithms, as they require problem-specific data structures and several advanced programming features from the CUDA API [START_REF] Gmys | IVM-based parallel branch-and-bound using hierarchical work stealing on multi-GPU systems[END_REF], [START_REF] Carneiro Pessoa | GPU-accelerated backtracking using CUDA dynamic parallelism[END_REF]. In turn, the GPUIterator module fulfills this requirement since it is designed to facilitate the Task 0 use of pre-compiled optimized kernels. However, its lack of distributed load balancing makes it unfeasible for the purpose of this work, since the irregular loads produced by tree search algorithms demand load balancing for efficient use of the computational resources [START_REF] Carneiro | Towards ultrascale branch-and-bound using a high-productivity language[END_REF]. In turn, Chapel's iterators are not well suited for the GPU programming model, as the loop body is executed for each element e of an iteration space. However, one element is not enough to use efficiently a massively parallel device such as a modern GPU. In order to cope with the challenge, we proceed as follows.

(locale 0) 1-2 1 -2 -3 1 -2 -4 1 -2 -5 1-5 1 -5 -2 1 -5 -3 x 1 -5 -4 1 -2 -3 1 -2 -4 1 -2 -5 1 -5 -2 ... . . .
Firstly, the master locale generates initial load through a partial search on CPU, filling a pool of nodes. Chapel's distributed iterators are used for work distribution and load balancing among locales. Then, for more efficient exploitation of the intra-node parallelism, each subproblem received by a worker locale is processed locally on CPU, by another partial search, generating a task-local pool. Then, this set of subproblems is divided by the number of existing GPUs and then offloaded to the local devices. Details concerning the master locale and the initial search on CPU are presented in the next subsection.

Algorithm 1: The Master locale.

1 N ← get problem( ) 2 cutof f ← get cutof f depth( ) 3 second cutof f ← get scnd cutof f depth( ) 4 P ← {} N ode 5 metrics ← (0, 0) 6 metrics + = initial search(N, cutof f, P )

B. The Master Locale and the Initial Search

The algorithm starts with task 0 running on locale 0. As one can see in Algorithm 1, task 0 initially receives the size N of the problem, the first cutoff depth, and the second one (lines 1 -3). As illustrated in Figure 1, the initial pool of nodes P is generated though a partial search (line 6), called initial search. This procedure is performed by task 0 and implicitly enumerates all feasible and valid incomplete solutions containing cutof f elements of the permutation, keeping them into the pool P .

Lines 7 to 9 are responsible for defining the distributed pool of nodes P d . As one can see at line 7, initially the domain D (0..|P | -1) is mapped onto the locales according to a distribution. Then, the distributed pool P d is defined over the domain D. Finally, as depicted in Figure 2, P d is initialized via implicit bulk transfer, represented by the assignment operator at line 10, which performs L-1 transfers through the network, where L is the number of locales on which the application runs, instead of |P d | small ones.

The parallel distributed search takes place in line 11, adding parallelism by using the forall statement along with distributed iterators (DistributedIters), which are responsible for the assignment of nodes in P d to locales in a master-worker manner (distributed load balancing). Finally, there is no need for programming termination criteria or a reduction of the search metrics. The search finishes when the distributed active set P d is empty, and metrics are reduced by using the reduction intents provided by Chapel (+ reduce).

C. Exploiting Intra-node Parallelism

As commented in Section III-A, one node yielded by the distributed iterator used in Algorithm 1 (line 11) is not enough for the efficient use of all GPUs a locale has. For the efficient use of the GPUs the locale is equipped with, we proceed as outlined in Algorithm 2.

Before the intra-locale search starts, a task-local pool P l of type N ode is defined (line 2). Next, a second partial search is performed from depth cutof f until second cutof f , also storing into P l all feasible and valid incomplete solutions found at second cutof f (line 3). Finally, the multi-GPU search takes place (line 4), searching for valid and complete solutions, i.e., from depth second cutof f until the problem size (N ). The algorithm returns a tuple (metrics) containing the explored tree size and the number of valid and complete solutions found by Algorithm 3. The Intra-locale Multi-GPU Search: The GPUs of the locale are exploited as shown in Algorithm 3. Initially, the vectors used by the CUDA kernel to retrieve the metrics from the GPU are defined in lines 1 -2, and the variable γ receives the number of GPUs the computer node on which the task is running has (line 3). In turn, lines 5 to 14 are responsible for launching the searches on GPU. Initially, one task gpu id is created for each GPU. Next, the id of the GPU onto which the task gpu id is going to offload is set though function cuda set gpu(), which calls the the CUDA function cudaSetDevice() though the C-interoperability layer. At line 7, the number of subproblems each GPU processes is calculated. The function get_load() simply Algorithm 3: Exploiting multiple GPUs.

N , P , the second cutoff depth Output: A tuple containing the explored tree size and the number of complete and valid solutions found on GPU. verifies whether the division GP U load γ is exact. If it is not exact, the GPU with its id gpu id equals to γ -1 receives the remainder of the division. Next, it is required to specify which slices of P , sols h and tree h each one of the γ kernels are going to work on (lines 8-11). The task-local variable stride represents the beginning of the range belonging to each GPU and is calculated as shown in [START_REF] Crainic | Parallel branch-and-bound algorithms[END_REF].

1 tree h ← [0..|P | -1] int 2 sols h ← [0..|P | -1] int 3 γ ← cuda
gpu id * ( GP U load/γ ) (1)
Once the starting positions on the vectors are known, the task-local variables sols ptr, tree ptr, and pool pts are defined pointing to the beginning of the range the GPU gpu id is responsible for working on. Then, the function written in C responsible for launching the kernel is called (line 12). It is important to point out that host -to -device, device -to -host, error checking, and other CUDA-specific functions are performed in this function. Finally, a parallel reduction of the metrics is performed (lines 15-16), and then, the metrics are returned (line 17-18). One can see in Figure 3 an overview of the proposed algorithm.

D. Implementation Aspects

The master-worker behavior is obtained by setting the coordinated flag of the distributed iterator to true. This way, locale 0 does not search but coordinates the search process. To avoid losing the computational resources of the computer on which locale 0 resides, one locale more than the number of computer nodes available is launched by using the -nl launching option. Thus, the first computer of the reservation hosts two locales: locale 0 and locale 1.

Concerning the intra-node parallelism, the pointers arithmetic of Algorithm 3 are performed by using Chapel's c_ptrTo() function. The CUDA-C kernel is based on a serial and hand-optimized backtracking for solving permutation combinatorial problems [START_REF] Carneiro Pessoa | GPU-accelerated backtracking using CUDA dynamic parallelism[END_REF], and it is a non-recursive backtracking that does not use dynamic data structures, such as stacks.

The semantics of a stack is obtained by using a variable depth and by trying to increment the value of the vector board at position depth. If this increment results in a feasible and valid incomplete solution, the depth variable is then incremented, and the search proceeds to the next depth. After trying all configurations for a given depth, the search backtracks to the previous one.

It is important to point out that the Node data structure is similar to any permutation-based combinatorial/combinatorial optimization problem. It contains two integer vectors of size cutof f . The first one, identified by board, stores the feasible and valid incomplete solution. The second one keeps track of board lines by setting position n to 1 each time a queen is placed in the n-th line.

IV. PERFORMANCE EVALUATION

A. Experimental Protocol

In our experimental study, the following applications for enumerating all valid and complete configurations of the N-Queens problem are considered.

• Baseline (CUDA-C): multi-GPU implementation optimized for single-locale execution written in CUDA-C. It consists of an initial search on CPU until a cutoff depth, which generates a pool of nodes P . Then, this set is divided into γ subsets, where γ is the number of GPUs the computer node is equipped with. Finally, after the kernel execution, the CPU part of the code reduces the metrics (tree size and number of found solutions). • GPUIterator: distributed version of the baseline implementation written in Chapel. The main difference between the GPUIterator and the baseline is that this variant uses the GPUIterator module [START_REF] Hayashi | GPUIterator: Bridging the gap between Chapel and GPU platforms[END_REF], which has been introduced in Section II. Differently from the algorithm proposed in this paper, the GPUIterator implementation performs no distributed load balancing. In turn, it first divides P d into L subsets, where L is the number of locales on which the application runs. Then, on each locale, the iterator further divides each set L i , i ∈ {0, ..., L -1} into γ subsets to be offloaded to the GPUs. • ChplGPU: implementation of the productivity-aware distributed search detailed in the last section, which employs Chapel's distributed iterators for work distribution and load balancing. It is important to point out that all implementations introduced above employ the same kernel written in CUDA-C 1 . Finally, it is worth mentioning that it is not the objective of this work to compare Chapel vs. MPI+X for implementing distributed tree search in terms of productivity, performance, and scalability, as such a study has already been performed [START_REF] Carneiro | Towards ultrascale branch-and-bound using a high-productivity language[END_REF]. 1 All implementations can be found on the Chapel-based Optimization Project (ChOp) repository: https://github.com/tcarneirop/ChOp Instead, the objective of this evaluation is to investigate how the Chapel-based applications scale as the number of GPUs increases.

B. Parameters Settings

N-Queens problems of size (N ) ranging from 17 to 21 are considered. The experiments take from a few seconds (N = 17) to ten hours of parallel processing on two GPUs (N = 21). The number of computer nodes considered in the experiments ranges from 1 to 12. Computer nodes operate under Debian 4.9.0, 64 bits. They are equipped with two Intel Xeon E5-2650 v4 @ 2.00 GHz (a total of 24 cores/48 threads per node) and 128 GB RAM. Each computer node is equipped with two NVIDIA GeForce GTX 1080 Ti -Pascal generation (11 GB RAM and 3584 CUDA cores @ 1582 Mhz). Thus, the maximum number of GPUs used in the experiments is 24 (86, 016 CUDA cores). The computer nodes are interconnected through a 100 Gbps Intel Omni-Path network.

The ChplGPU and GPUIterator implementations are programmed with Chapel 1.22, and the default task layer (qthreads) is the one employed. The GPU kernel is programmed in CUDA and compiled with NVCC 10.1, whereas the function that calls the kernel is programmed in C and compiled with gcc 8.3. Both implementations written in Chapel are executed on top of GASNet, and appropriate environment variables should be set to fully utilize the capability of the target platform. Table I shows a summary of the runtime configurations for multi-locale execution. The OpenFabrics GASNet implementation is the one used for communication (CHPL_COMM_SUBSTRATE) along with SSH, which is responsible for getting the executables running on different locales (GASNET_PSM_SPAWNER).

Chapel provides two different distributed load balancing iterators: guided and dynamic, which are also similar to OpenMP's schedules of the same name. Experiments were also carried out to identify the best chunk for both load balancing strategies. As depicted in Figure 4, using the dynamic distributed iterator results in the best overall performance for ChplGPU. Therefore, it is the distributed iterator to be considered hereafter in the results presented for the ChplGPU implementation.

Experiments were also carried out to choose suitable cutoff depths for all implementations. This parameter directly influences the size of the pool of nodes, the granularity of the nodes, the GPU occupancy and the capacity of load balancing of the iterators [START_REF] Carneiro | Towards ultrascale branch-and-bound using a high-productivity language[END_REF], [START_REF] Carneiro Pessoa | GPU-accelerated backtracking using CUDA dynamic parallelism[END_REF], [START_REF] Jenkins | Lessons learned from exploring the backtracking paradigm on the GPU[END_REF]. One can see in Table II the best parameters experimentally found for all implementations introduced in Section IV-A

Parameter

ChplGPU GPUIterator CUDA-C cutof f

2 7 7 second cutof f 8 - - Distributed iterator dynamic - - Distributed chunk 1 - - C. Performance Results
1) Single-locale Execution: Figure 5 shows the singlelocale performance of both implementations written in Chapel compared to the CUDA-C baseline. Firstly, the GPUIterator implementation is equivalent to the baseline one in terms of performance. While it is 3% slower than the baseline when N = 17, the performance difference is not significant as N increases. This indicates that the use of the GPUIterator module does not incur significant overhead. However, this is not the case with ChplGPU. For example, ChplGPU is more than twice as slow as the baseline when N = 17.

The CUDA-C implementation is straightforward and optimized for single-locale execution. The initial search on CPU creates the active set and divides it into γ sets, one for each GPU, and only γ kernel calls are performed. The same is true for GPUIterator. In turn, ChplGPU has several sources of overhead for a single-locale execution. It consists of two partial searches on CPU before starting the search on GPU. The iterator distributes nodes to the worker locales, which perform another partial search on CPU for generating a tasklocal pool of nodes (P t ). In turn, local active sets are much smaller than the single one generated by CUDA-C, which results in less efficient use of the GPUs for the smallest sizes [START_REF] Jenkins | Lessons learned from exploring the backtracking paradigm on the GPU[END_REF].

It is also worth to point out that, unlike the CUDA-C baseline and the GPUIterator implementation that only launch γ kernels, ChplGPU launches |P d | * γ kernels. For instance, 0 0.5 2GPUs) to 12 computer nodes (24 GPUs). In the graph, the keys present the number of computer nodes followed by the number of GPUs in parenthesis -e.g., 12(24) means that the results are shown for 12 computer nodes and a total of 24 GPUs.

ChplGPU launches 240 * 2 kernels for N = 17. Additionally, each kernel launch requires host-to-device, device-to-host transfers, and metrics reductions on CPU, which is also a source of overhead. As the problem size grows, the overhead becomes less significant compared to the solution space to be evaluated. Therefore, ChplGPU, which is 2.15× slower than the baseline for N = 17, is only 15% and 10% slower for sizes N = 20 and 21, respectively.

2) Multi-locale Execution: Figure 6 presents the normalized execution times of ChplGPU and GPUIterator solving the N-Queens problems using from 4 to 24 GPUs. Firstly, consider the results for the smallest instance (N = 17). As solving this instance takes only a few seconds, GPUIterator is faster than ChplGPU for 4 to 20 GPUs. As more computer nodes are added, GPUIterator performs poorly due to the lack of distributed load balancing. In turn, the performance difference between ChplGPU and GPUIterator grows as more computer nodes are added. For N = 18 and 16 -24 GPUs, ChplGPU is from 1.32× -1.45× faster than its counterpart. As the problem size grows, GPUIterator becomes slower than ChplGPU even on 2 computer nodes -4 GPUs. For N = 20, ChplGPU is from 1.13× (4 GPUs) to 1.77× (24 GPUs) faster than the GPUIterator implementation.

Figure 7 shows the strong scaling of ChplGPU and GPUIterator on 2 GPUs (1 computer node) to 24 GPUs (12 computer nodes) over the CUDA-C baseline on 2 GPUs (1 computer node). Note that for each instance, the performance numbers are normalized to Baseline(2). In general, as the number of GPUs increases, the performance of both ChplGPU and GPUIterator improves except for the ChplGPU on 2 GPUs as with Figure 5. Also, ChplGPU achieves further performance improvements over GPUIterator as N , and the number of GPUs increases, due to the dynamic scheduling feature, which is not supported in GPUIterator.

To quantify these performance trends, we compare the performance numbers with the linear speedup -i.e., the CUDA-C baseline on nGP U s gives n× speedup. Initially, consider the results for ChplGPU (Figure 7a). When N ≥ 19, the speedups achieved for all < problem, #GP U S > configurations are at least 80% of the linear speedup. For N ≥ 20, results are even better, ranging from 87% to 91% of the linear speedup. In turn, due to its static load distribution scheme, the GPUIterator-based implementation reaches around 80% of the linear speedup for all sizes when only 4 GPUs (2 computer nodes) are used (refer to Figure 7b). For N ≥ 18, the speedups achieved by the GPUIterator implementation are around 50% of the linear speedup.

V. DISCUSSION

Initially, if only the single-locale programming is taken into account, the GPUIterator module is indeed an option in terms of time to a first implementation. Error-prone details such as division of the pool among locales and GPUs, pointers arithmetic, and dealing with the CUDA API are hidden to the programmer. As a consequence, the size of the single-locale version of the GPUIterator implementation has almost the same number of source lines of code (SLOCS) of the CUDA-C baseline implementation. Moreover, it also enables one to conceive a GPU-based distributed implementation by adding a few lines of code. For instance, an addition of 8 SLOCS is required to make GPUIterator a distributed application. However, the lack of load balancing of the GPUIterator module results in poor scalability, which justifies the more complex implementation of ChplGPU.

In terms of SLOCS, ChplGPU is 1.5× longer than the GPUIterator implementation. The partial searches necessary to exploit the distributed iterators, along with the Cinteroperability code related to the use of GPUs, amount for more than 50% of the SLOCS of ChplGPU. In the scope of this work, the higher programming effort of coding two levels of partial search and the functions responsible for offloading and getting data from the GPU pays off, as ChplGPU achieves parallel efficiency up to 2× higher and it is up to 1.77× faster than its Chapel-written counterpart. In this sense, a further necessary improvement to the GPUIterator module is implementing distributed GPU-oriented load balancing iterators. This effort would make it possible to program a distributed tree search that scales but having only a fraction of the size of ChplGPU.

Finally, it is important to mention that some of the challenges concerning employing GPUs for irregular tree search remain when using a high-productivity language. On the one hand, we can exploit high-level features provided by Chapel for several aspects of the search, such as work distribution and termination criteria. On the other hand, mixing programming models also brings GPU's challenges along with it. For instance, it is required to tune the chunk size of the distributed iterator taking into account that a subproblem yielded by the iterator must provide load enough for multiple GPUs. The performance of ChplGPU can be drastically decreased due to a bad parameter choice. This is also true for GPUIterator. Despite its straightforward implementation, it is still required to tune the depth parameters to generate a distributed pool that, after divided into L * γ subsets also can provide load enough to the GPUs, so they can be used efficiently.

VI. CONCLUSIONS AND FUTURE WORK

In this work, we revisited the design and implementation of distributed tree search algorithms dealing with multiple GPU accelerators, in addition to scalability and productivityawareness, using Chapel. A distributed backtracking for enumerating all valid solutions of the N-Queens problem was conceived. The proposed algorithm exploits Chapel's distributed iterators combined with pre-compiled CUDA kernels through Chapel's C-interoperability layer.

According to the performance results, the implementation of the proposed algorithm is slower than a hand-optimized CUDA-C implementation, when taking into account singlelocale execution. However, in distributed scenarios, the proposed algorithm scales, achieving more than 90% of the linear speedup in the biggest test-cases.

Permutation combinatorial optimization problems are commonly solved to optimality by using branch-and-bound search algorithms [START_REF] Crainic | Parallel branch-and-bound algorithms[END_REF]. Therefore, the first future research direction is to extend the proposed multi-locale backtracking into a distributed B&B. This way, it will be possible to solve challenging optimization problems, such as the Quadratic Assignment and the Flow-shop Scheduling Problem.

Going one step further towards exascale, the scalability should be increased considering more GPU-powered processing nodes. Also, heterogeneity can be improved to consider other accelerators such as FPGAs thanks to Chapel's interoperability capabilities. Finally, checkpointing-based fault tolerance should be addressed in addition to scalability, productivity-awareness, and CPU-GPU heterogeneity issues to achieve a holistic exascale-aware design and implementation.

xFig. 1 .

 1 Fig. 1. Schematic representation of an initial search on locale 0 -task 0 that generates the pool P for a problem size N = 4 and cutof f = 3. The figure depicts the branch that has the element 1 of the permutation as the root and generated 4 valid and feasible incomplete solutions at depth cutof f = 3.

7Fig. 2 .

 2 Fig.2. The distributed pool P d consists of several sets P i d , i ∈ {0, ..., L-1}, where L is the number of locales on which the application runs (based on[START_REF] Carneiro | Towards ultrascale branch-and-bound using a high-productivity language[END_REF]).

Algorithm 2 :

 2 Generating the task-local pool. Input: N , node, cutof f , scnd cutof f Output: A tuple containing the explored tree size and the number of complete and valid solutions found on GPU. 1 metrics ← (0, 0) 2 P l ← {} N ode // Task-local pool 3 metrics + = partial search(N, node, cutof f, scnd cutof f, P l ) 4 metrics + = Algorithm 3(N, P l , scnd cutof f ) 5 return metrics

Fig. 3 .

 3 Fig.3. The locale 0 (master) is responsible for generating the distribute pool P d and controlling the search. Each worker locale receives nodes from the master and generates a local pool (P l ) that is partitioned into γ subsets. As said in Section III-D, L locales are launched on L -1 computer nodes (Own representation adapted from[START_REF] Crainic | Parallel branch-and-bound algorithms[END_REF]).

Fig. 4 .

 4 Fig. 4. Performance comparison of the ChplGPU implementation using two different distributed iterator and static distribution of P d . Results are shown for 24 GPUs (12 computer nodes).

Fig. 5 .

 5 Fig. 5. Normalized execution time of the optimized CUDA-C implementation compared to its Chapel-based counterparts. Results are for one one computer node -two GPUs.

Fig. 6 .

 6 Fig. 6. Normalized execution time of GPUIterator and ChplGPU. Results are shown for 4 GPUs (2 computer nodes) to 24 GPUs (12 computer nodes).

Fig. 7 .

 7 Fig. 7. Speedup achieved by (a) ChplGPU and (b) GPUIterator implementations compared to the optimized CUDA-C baseline executed on one computer node (two GPUs). Results are shown for 1 computer node (2GPUs) to 12 computer nodes (24 GPUs). In the graph, the keys present the number of computer nodes followed by the number of GPUs in parenthesis -e.g.,12(24) means that the results are shown for 12 computer nodes and a total of 24 GPUs.

TABLE II BEST

 II PARAMETERS EXPERIMENTALLY FOUND FOR THE CHPLGPU, GPUITERATOR AND THE OPTIMIZED CUDA-C IMPLEMENTATION.

https://gitter.im/chapel-lang/chapel

ACKNOWLEDGMENTS

The experiments presented in this paper were carried out on the Grid'5000 testbed [19], hosted by INRIA and including several other organizations. We thank Bradford Chamberlain, Elliot Ronaghan, Vassily Litvinov, Michael Ferguson, and Lydia Duncan (Hewlett Packard Enterprise -HPE) for their support on Chapel's Gitter platform 2 .


M Aucejo

O De Smet

A generalized multiplicative regularization for input estimation

Keywords: Linear inverse problem, Force reconstruction, Multiplicative regularization

This paper implements a generalized multiplicative regularization for estimating the mechanical loads acting on a linear structure. The proposed strategy extends the ordinary multiplicative regularization, previously published by the authors, by introducing an extra tuning parameter, which is determined through an original iterative procedure. To assess the practical interest and the overall performances of the proposed approach, numerical and real-world applications are proposed. Obtained results illustrate the influence of the extra tuning parameter according to the measurement noise level and highlight the benefits brought by the generalized multiplicative regularization in terms of solution accuracy.

Introduction

Input estimation remains an important problem for the structural dynamics community as evidenced by the abundant literature dedicated to this topic in the recent years. Generally speaking, inverse methods can be classified into two groups. The first group includes methods that are specifically designed to solve the inverse problem in the time or frequency domains.

In the time domain, one can cite Kalman-like approaches [START_REF] Lourens | An augmented Kalman filter for force identification in structural dynamics[END_REF][START_REF] Maes | Joint input-state estimation in structural dynamics[END_REF][START_REF] Feng | Force localization and reconstruction based on a novel sparse Kalman filter[END_REF] or dynamic programming [START_REF] Busby | Optimal regularization of an inverse dynamics problem[END_REF][START_REF] Nordström | A dynamic programming algorithm for input estimationon linear time-variant systems[END_REF][START_REF] Gonzàlez | A general solution to the identification of moving vehicle forces on a bridge[END_REF] , while, in the frequency domain, methods based on the filtering of the equation of motion of structures, such as beams, cylindrical shells or plates, have been developed [START_REF] Zhang | Examples of using structural intensity and the force distribution to study vibrating plates[END_REF][START_REF] Pezerat | Force analysis technique: reconstruction of force distribution on plates[END_REF][START_REF] Djamaa | Mechanical Radial Force Identification of a Finite Cylindrical Shell by an Inverse Method[END_REF] The second group gathers more general approaches that can be used indistinctly in both domains. This is, for instance, the case of methods based on the virtual work principle [START_REF] Berry | Identification of dynamic loading on a bending plate using the Virtual Fields Method[END_REF][START_REF] Xu | Force identification of dynamic systems using virtual work principle[END_REF][START_REF] O'donoughue | Time-resolved identification of mechanical loadings on plates using the virtual fields method and deflectometry measurements[END_REF]. However, among all the methods belonging to this second group, regularization techniques are certainly the most widely used, because it allows including some prior information on the sources to identify. From a mathematical standpoint, the latter information is encoded in a regularization term generally incorporated in the formulation as an additive constraint to give rise to the so-called Tikhonov-like regularization [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF]. Although this form of additive regularization has proved its efficiency to solve force reconstruction problems [START_REF] Thite | The quantification of structure-borne transmission paths by inverse methods. Part 2 : Use of regularization techniques[END_REF][START_REF] Qiao | Sparse deconvolution for the large-scale ill-posed inverse problem of impact force reconstruction[END_REF][START_REF] Pan | Moving force identification based on redundant concatenated dictionary and weighted 1norm regularization[END_REF] , it requires the determination of the well-known regularization parameter. Unfortunately, this task reveals to be far from easy and automatic selection procedures have been developed for this purpose, such as the Generalized Cross-Validation [START_REF] Golub | Generalized Cross-Validation as a method for choosing a good ridge parameter[END_REF], the L-curve principle [START_REF] Hansen | Rank-deficient and discrete ill-posed problems: numerical aspects of linear inversion[END_REF] or the Bayesian indicator [START_REF] Pereira | Empirical bayesian regularization of the inverse acoustic problem[END_REF] to cite only only a few of them . However, it should be mentioned that all the above-cited selection techniques are generally computationally intensive and have been developed in the context of Tikhonov regularization (a.k.a. 2 -regularization). It results that they cannot be easily or directly applied to more general additive regularizations, such as 1 [START_REF] Efron | Least angle regression[END_REF],

q [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF] or mixed-norm regularizations [START_REF] Rezayat | Identification of dynamic forces using group-sparsity in frequency domain[END_REF]. To circumvent to this potential undesirable feature, we have recently introduced the ordinary multiplicative regularization in the context of structural source identification [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF][START_REF] Aucejo | A space-frequency multiplicative regularization for force reconstruction problems[END_REF][START_REF] Aucejo | Multi-parameter multiplicative regularization: An application to force reconstruction problems[END_REF]. In its more general form, the proposed multiplicative regularization is expressed as:

F = argmin F\{0} F(X, F) • R(F), (1)
where:

• X and F are the measured vibration field and the unknown excitation field respectively;

• F(X, F) is the data-fidelity term, which is a measure of the deviation of the estimated vibration field from the measured one. It actually encodes prior information related to the noise corrupting the data [START_REF] Le | A variational approach to reconstructing images corrupted by poisson noise[END_REF];

• R(F) is the regularization term that encodes prior information on the excitation field F [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF].

One of the remarkable properties of the previous formulation is that the corresponding solution is the point of intersection between the L-curve and a straight line with a slope equal to -1 [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF]. In this regard, it can be interesting to add an extra tuning parameter that brings a greater flexibility by modifying the slope of the straight line and so the identified solution.

Such an idea has been introduced by Regińska [START_REF] Regińska | A regularization parameter in discrete ill-posed problems[END_REF] and further developed by Viloche Bazán [START_REF] Bazán | Fixed-point iterations in determining the Tikhonov regularization parameter[END_REF] in the context of Tikhonov regularization to only compute the corresponding optimal regularization parameter α. In those papers, the regularization parameter is sought such that:

α = argmin α∈ R + * F(X, F α ) • R(F α ) µ , (2)
where F α is the solution of the Tikhonov regularization for a given regularization parameter α, while µ ∈ R + * is an extra tuning parameter. Practically, Eq. ( 2) is solved from a fixed-point algorithm that automatically adjusts the value of the extra tuning parameter [START_REF] Bazán | Fixed-point iterations in determining the Tikhonov regularization parameter[END_REF].

In the light of these works, this paper proposes to revisit the ordinary multiplicative regularization by introducing an extra tuning parameter in the framework we have previously developed to efficiently solve input estimation problems. This particular form of multiplicative regularization, referred to as Generalized Multiplicative Regularization (GMR) in the rest of the paper, differs from the works of Regińska and Viloche Bazán by its philosophy, since it considers more general regularization terms and aims at solving the formulation as a whole. In particular, this means that the identified excitation field, the related regularization parameter and the extra tuning parameter are computed within the same iterative resolution procedure. Hence, the basic motivation of this paper is to assess the applicability of the GMR for identifying mechanical loads acting on a structure. To this end, the mathematical formulation of the GMR and the related resolution algorithm are detailed in section 2. In section 3, a numerical experiment is carried out to thoroughly analyze the influence of the extra tuning parameter on the identified solution according to the measurement noise level. Based on these results, an original iterative procedure is proposed in section 4 to automatically compute an optimal value of the extra tuning parameter. Finally, a real-world experiment is carried out in section 5 to assess the ability of the whole resolution procedure in providing consistent reconstructions in operating conditions.

Generalized Multiplicative Regularization

This section introduces the mathematical formulation of the Generalized Multiplicative Regularization (GMR), as well as its resolution algorithm.

More specifically, this section highlights the main features of the GMR and focuses on the basic steps of the resolution process.

General formulation

As stated in the introduction, the GMR follows the idea proposed by

F = argmin F\{0} F(X, F) • R(F) µ . (3)
To explicitly define the data fidelity and regularization terms, F(X, F)

and R(F), three assumptions are made. The first one consists in considering the studied structures as linear and time-invariant systems. If we further assume that the mechanical problem is expressed in the frequency domain, then the dynamic behavior of the structure is completely determined by the transfer functions matrix H, relating the vibration field X to the excitation field F, such that:

X = HF + N, (4)
where N models the noise corrupting the measured data.

The second assumption is related to the statistical properties of the noise vector previously defined. Indeed, it has been recalled in the very beginning of this paper that the data-fidelity encodes any prior information available on the measurement noise. In this respect, the data-fidelity term has to reflect the main characteristics of the actual noise. Consequently, if the noise is supposed Gaussian and spatially white, the data-fidelity term is expressed as [START_REF] Antoni | A Bayesian approach to sound source reconstruction: optimal basis, regularization, and focusing[END_REF][START_REF] Aucejo | An optimal bayesian regularization for force reconstruction problems[END_REF]:

F(X, F) = X -HF 2 2 . (5)
The last assumption concerns the definition of the regularization term, which encodes one's prior knowledge of the spatial distribution of the excitation sources over the structure. Generally, forces of different nature can act simultaneously on a structure. As a result, it is supposed that the structure is excited in R different regions by local excitation fields F r of various types (localized or distributed, for instance). Formally, this naturally leads to introduced local regularization terms R(F r ), which are defined such that:

R(F r ) = F r qr qr , (6)
where • qr is the qr -norm (or quasi-norm), while q r is the related norm parameter defined in R + * . The latter parameter allows promoting either a sparse (or localized) excitation field when q r ≤ 1 or a smooth (or distributed) one when q r ≥ 21 .

Formally, the local regularization terms can be combined in an additive or a multiplicative manner [START_REF] Aucejo | Multi-parameter multiplicative regularization: An application to force reconstruction problems[END_REF]. The one retained in this paper is based on a Bayesian analysis and consists in defining the regularization term as the sum of the local regularization terms, that is:

R(F) = R r=1 F r qr qr . (7)
By bringing everything together, the general form of the GMR considered in this paper is finally written:

F = argmin F\{0} X -HF 2 2 • R r=1 F r qr qr µ . (8)
At this stage, two comments can be made regarding the global behavior of the GMR. First of all, whatever the particular form of the data-fidelity and regularization terms, the solution obtained from the GMR is the point at which the corresponding L-curve is tangent to a straight line having a slope equal to -1/µ. A general proof of this statement can be found in Ref. [START_REF] Regińska | A regularization parameter in discrete ill-posed problems[END_REF].

Then, the solution obtained from the GMR is directly conditioned to the value of the extra tuning parameter µ. Indeed, when µ = 0, the GMR reduces to an Ordinary Least Squares regression, whereas it reduces to the Ordinary Multiplicative Regularization (OMR) when µ = 1 [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF][START_REF] Aucejo | Multi-parameter multiplicative regularization: An application to force reconstruction problems[END_REF]. More generally, solutions minimizing the data-fidelity term are promoted when µ < 1, while solutions minimizing the regularization term are favored when µ > 1. To some extent, µ can be seen as a trade-off parameter, since it allows weighting the relative importance of the data-fidelity and regularization terms.

Resolution algorithm

The resolution of the GMR requires the implementation of an iterative procedure. From a very general standpoint, the resolution algorithm is divided into three main steps:

1. set k = 0 and initialize

F (0) 2. while convergence is not reached do a. Compute F (k+1)
b. Compute the convergence indicator end while

return F

In the present paper, the initialization step as well as convergence monitoring are not described and the interested reader can refer to Ref. [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF] for further information. Actually, the core of the proposed algorithm is the computation of the excitation field at iteration k + 1 (step 2.a), around which all the procedure is built. Practically, it is derived from the direct application of the first-order optimality condition to the functional:

J(F) = X -HF 2 2 • R r=1 F r qr qr µ . (9)
In doing so, the fixed-point iterate at iteration k + 1 is written:

F (k+1) = H H H + α (k+1) W (k+1) -1 H H X, (10)
where α (k+1) is the adaptive regularization parameter, while W (k+1) is a global weighting matrix. Formally, the adaptive regularization parameter α (k+1) is defined by the relation:

α (k+1) = µ X -H F (k) 2 2 R r=1 F (k) r qr qr . ( 11
)
On the other hand, the global weighting matrix W (k+1) is defined such that:

W (k+1) = diag W (k+1) 1 , . . . , W (k+1) r , . . . , W (k+1) R . ( 12
)
In the previous relation, each local weighting matrix W (k+1) r is a diagonal matrix given by:

W (k+1) r = diag w (k+1) r,1 , . . . , w (k+1) r,i
, . . . , w

(k+1) r,Nr (13)
with

w (k+1) r,i = q r 2 max r , F (k) ri qr-2 , ( 14
)
where N r is the number of identification point in the region r, F

ri is the i th component of the local excitation vector F (k) r , while r is a small positive number avoiding infinite weights when F (k) ri → 0 and q r < 2. The value of this parameter is chosen so that 5% of the values of F (0) are less than or equal to r [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF].

Numerical experiment

The numerical experiment proposed in this section aims at analyzing the behavior and the performances of the GMR on an academic application.

More specifically, the main objective of this application is twofold. First, it is intended to assess the influence of the extra tuning parameter µ on the solutions identified by the GMR. Second, it is sought to determine whether an overall trend emerges from the simulation results regarding the evolution of µ with respect to the measurement noise level.

Test case description

The present numerical experiment aims at identifying a unit point force Finally, it is necessary to define the reference force vector F ref , that could serve as a benchmark, as well as a set of qualitative indicators measuring the accuracy of the identified solutions. In the present case, the reference force vector must include the external point force as well as the reaction forces at boundaries, since their reconstruction is the primary objective of the proposed numerical experiment. This is for this particular reason that the reference force vector is computed as follows:

F ref = H -1 X ref . (15)
Fig. 1 presents the spatial distribution of the sources to identify. As expected from the problem description, the reference force vector exhibits smooth reaction forces at boundaries of the plate as well as a unit point force F 0 at (x 0 , y 0 ) = (0.42 m, 0.25 m).

GRE = F -F ref 1 F ref 1 . (16)
In a similar fashion, the relative error on the reaction forces is an indicator of the reconstruction quality of the reaction forces. Its mathematical definition is similar to the global relative error, since:

RERF = F (rf) -F (rf) ref 1 F (rf) ref 1 , (17)
where F (rf) is the force vector corresponding to the identified reaction forces at boundaries , while

F (rf)
ref is the related reference force vector. Contrary to the foregoing indicators, the peak error is a local indicator describing the reconstruction quality of the point force amplitude. Mathematically, it is defined such that:

PE = F p -F ref p F ref p , (18)
where F ref p is the point force amplitude associated to the reference force vector F ref , while F p is the point force amplitude associated to the identified solution

F at point (x 0 , y 0 ).
Finally, it worth mentioning that the quality of the identified excitation field is conditioned to the convergence of FE mesh used to build the the transfer functions matrix. Indeed, at a particular frequency, the FE model must be fine enough to properly describe the dynamic behavior of the structure, which is the case here. Another possibility, not used in this work, consists in computing the transfer functions matrix for a converged mesh and then applying selection matrices to fit with the measurement and reconstruction meshes.

Application

In this section, the behavior and the performances of the GMR with respect to the value of the extra tuning parameter µ and the measurement noise level are investigated. More precisely, by applying the GMR for various SNR values ranging from 30 dB to 0 dB (i.e. from weakly to extremely noisy data), one expects to see a trend emerge for the parameter µ with the measurement noise level. But before that, it remains to characterize the identification regions as well as the value of the related norm parameters q r . Following the description of the problem and the distribution of the sources presented in Fig. 1, two identification regions can be defined: (i) a central region associated to the norm parameter q 1 and including the external point force only and (ii) a region corresponding to the boundaries of the plate and associated to the norm parameter q 2 [see Fig. 2]. Practically, this division invites us to promote the sparsity of the excitation field in the central region and its continuity at boundaries. From the indications given in section 2.1 and our own experience of such problems [START_REF] Aucejo | An optimal bayesian regularization for force reconstruction problems[END_REF], choosing q 1 = 0.5 and q 2 = 2 fulfills this requirement.

To properly analyze the influence of µ on the reconstruction accuracy, the results obtained from the Ordinary Multiplicative Regularization (OMR), corresponding to the GMR for µ = 1, are first presented in Table 1 and Fig. reconstruction of the point force and that of the reaction forces. This is especially noticeable for extremely noisy data (i.e. SNR = 0 dB). In this case, the point force location is identified at ( x 0 , y 0 ) = (0.46 m, 0.25 m) instead of (x 0 , y 0 ) = (0.42 m, 0.25 m), corresponding to a location error less than 10%.

On the other hand, the PE computed from the point force amplitude identified at ( x 0 , y 0 ) is equal to 0.973 N instead of 1 N (see Table 2). Consequently, all things considered, this is a remarkable result for such an extreme configuration. Finally, it is worth noting that the optimal value of the extra tuning parameter decreases as the measurement noise level increases. More specifically, µ should be greater than 1 in case of weakly and moderately noisy data, while it should be less than or equal to 1 for highly and extremely noisy data. This observation perfectly reflects the subtle balance induced by the regularization process, i.e. the relative weighting of the data-fidelity and regularization terms controlled by the adaptive regularization parameter α. ( * ) PE computed from the amplitude at the identified point force location

Automatic computation of the extra tuning parameter

The results obtained in the previous section have shown that the optimal value of the extra tuning parameter decreases with the SNR. However, finding this value a priori is far from an easy task. This observation points out the need for implementing an automatic selection procedure. In this section, an original fixed-point algorithm is introduced for this purpose.

Basic principle

A careful analysis of the resolution algorithm, described in section 2.2,

shows that an optimal solution should be obtained provided that [see Eq. ( 11)]:

α R r=1 F r qr qr X -H F 2 2 ≈ µ, ( 19
)
where α is the value of the adaptive regularization parameter obtained at convergence of the resolution algorithm.

This observation suggests the definition of the following function:

ψ(α, µ) = α R r=1 F rα qr qr X -HF α 2 2 , ( 20
)
where F α is the solution of the additive counterpart of the GMR obtained for a constant regularization parameter α, that is 3 :

F α = argmin F X -HF 2 2 + α R r=1 F r q q µ , (21)
while F rα is the subset of F α restricted to the region r.

It results from Eqs. ( 19) and ( 20) that an optimal value of the extra regularization parameter µ should be found when:

µ = ψ ( α, µ) , (22)
which is equivalent to seek the solution such that F = F α .

This very basic observation suggests the implementation of the iterative algorithm presented in the next section and satisfying Eq. ( 22) within some tolerance when the iterative process has converged.

3 F α is computed from the algorithm presented in section 2.2. However, the actual

regularization parameter at iteration k + 1 is α (k+1) = µ α R r=1 F (k) rα qr qr µ-1
and not simply α.

Practical implementation

Before detailing the practical implementation of the proposed iterative procedure, it is important to fully understand its expected behavior. To this end, let us take another look at the numerical experiment presented in section 3 for an SNR of 30 dB. The starting point of our reflection is illustrated in Fig. 5 presenting the plot of ψ (α, µ 0 ) for µ 0 = 1 (µ 0 : initial choice of the extra tuning parameter). This figure shows that the point ( α 0 , µ 0 ), indicated by the marker ( * ), does not lie on the curve defined by ψ(α, µ 0 ), meaning that the corresponding initial solution ( F 0 , α 0 ) is not optimal in the sense of Eq. ( 22). For finding a value of µ adapted to the problem, it is chosen to apply the GMR once again using an updated value of µ, noted µ 1 and chosen such that µ 1 = ψ ( α 0 , µ 0 ) (see Fig. 6a). In doing so, the new point ( α 1 , µ 1 ), associated to the updated solution ( F 1 , α 1 ), is almost located on the curve defined by ψ(α, µ 1 ) (see Fig. 6b), meaning that F 1 ≈ F α 1 , i.e. µ 1 ≈ ψ( α 1 , µ 1 ). Following this first iteration, the process can be continued until the relative variation of µ between two successive iterations is less than or equal to some tolerance fixed by the user. The whole fixed-point algorithm used to determine the optimal triplet ( F, α, µ) is summarized by the pseudo-code presented in Algorithm 1.

-

Application

For the numerical experiment considered in section 3, the application of the algorithm previously described for a tolerance set to 10 -3 (i.e. 0.1%) allows obtaining the results gathered in Table 3 and Fig. 7. From the analysis of these results, some comments can be made. First of all, the value of the optimal extra tuning parameter found by the algorithm tends to decrease as the measurement noise level increases. This observation is actually in line with the conclusion drawn in section 3. Then, the value of µ, retained by Algorithm 1: Pseudo-code of the proposed iterative algorithm Inputs : X, H, q r , tol

Outputs: F, α, µ

Set j = 0, δ = 1, µ 0 = 1
Compute ( F 0 , α 0 ) from the GMR for µ 0 (see section 2.2)

while δ > tol do µ j+1 = ψ( α j , µ j ) Compute ( F j+1 , α j+1 ) from the GMR for µ j+1 δ = |µ j+1 -µ j | µ j j = j + 1 end F = F j , α = α j , µ = µ j
the algorithm, leads to a reasonable compromise between the reconstruction accuracy of the point force and that of the reaction forces. Obviously, better solutions can be obtained by manually tuning µ (see section 3), but this is a perilous undertaking without a precise knowledge of the target excitation field. Finally, it should be noted that, for the two lowest SNR values (i.e.

10 dB and 0 dB), the proposed algorithm either has difficulty converging or even diverges due to numerical issues during the process.

However, in the light of the previous results and those presented in Fig. 6, it is legitimate to question the interest of continuing the iterative process beyond the first iteration in terms of solution accuracy and computational efficiency. To address this issue, the results obtained after the first iteration of the algorithm are given in Table 4 and Fig. 8. It can be especially noted

Real-world experiment

In this section, a real-world application is carried out on a simple structure in order to assess the identification ability of the GMR in operating conditions.

Description of the experimental set-up

The structure under test is a steel parallelepiped box, excited on one of its faces by a shaker fed by a white noise signal and equipped with a force sensor (see Fig. 9). The dimensions and the thickness of the excited face are (L y , L z ) = (0.3 m, 0.35 m) and h = 5 mm respectively, while the excitation point is approximately located at (y 0 , z 0 ) = (0.10 m, 0.09 m).

As explained in section 3, the application of the GMR requires a measured vibration field X and a transfer function matrix H as input data. In the present experiment, the vibration field has been recorded by a scanning laser vibrometer on a grid of (N y × N z ) = (19 × 22) measurement points using the force sensor as phase reference. Furthermore, the measured velocity field is, in the following, normalized to the signal delivered by the force sensor to obtain the vibration field resulting from the application of a unit point force (i.e. F ref p = 1). Regarding now the definition of the transfer functions matrix, it has been chosen to derive it from the FE model of a free plate having the same dimensions and material as the excited face, assuming that the measured bending motions are the only available data. Accordingly, the FE mesh of the equivalent plate has been built from the measurement grid and is made up with 378 shell elements4 .

Application

The analysis of the experimental set-up suggests the definition of two identification regions similar to those defined in the numerical application.

More precisely, the first region, associated to the norm parameter q 1 , contains the external point force only, while the second one, associated to the norm parameter q 2 , corresponds to the boundaries of the excited face [see Fig. 10].

As explained in section 3.2, the values of the norm parameters must be chosen so as to reflect the spatial distribution of the forces to identify in the considered regions, i.e. the sparsity in the first region and the continuity of the boundary forces in the second region. Following the latter observation, the norm parameters q 1 and q 2 are respectively set to 0.5 and 2. In the present experiment, the GMR is applied at 625 Hz. This specific frequency has been chosen because it is not a natural frequency of the sys-tem 5 . In the following, because no reference excitation field is available, only the PE is evaluated from the results obtained for three values of the extra tuning parameter, corresponding respectively to the value used to initialize Algorithm 1 (µ 0 ), that computed after the first iteration (µ 1 ) and that reached at convergence ( µ) for a tolerance set to 10 -3 . Table 5 and Fig. 11 gather the results derived from the GMR. First of all, it is important to note that the least-squares solution is obtained by letting Algorithm 1 reach the convergence, which is obviously not the expected behavior. In the two other cases, the spatial distribution of the equivalent boundary forces are consistent with our expectations 6 and are somewhat equivalent whatever the value considered. On the contrary, the reconstruction of the excitation fields in the central region leads to noticeable differences. Indeed, Fig. 11 shows that the OMR (i.e. µ 0 = 1) allows reasonably well quantifying and locating the external point force despite the presence of small reconstruction artifacts, since it is identified at 5 At natural frequencies of a lightly damped structure, the mechanical response is driven by the modes of the system. This implies that the solution is generally not unique at these particular frequencies. 6 Because a force-only model has been adopted here, the identified boundary forces are actually the sum of pure forces and equivalent forces due to the potential existence of boundary moments.

( y 0 , z 0 ) = (0.103 m, 0.088 m) with a normalized amplitude F p of 0.98 (see Fig. 11a). Actually, the reconstruction artifacts, observed previously, are eliminated by applying the GMR for µ = µ 1 , while leading to a better quantification of the point force amplitude (see Table 5 and Fig. 11b). These findings are in line with the conclusion drawn in section 3.2 and suggesting that completing the iterative process is not necessarily beneficial in terms of solution accuracy. More specifically, stopping the process right after the first iteration allows obtaining consistent solutions while limiting the computational cost of the iterative scheme described in Algorithm 1.

Conclusion

A generalized multiplicative regularization has been proposed for solving linear input estimation problems. The proposed method extends the ordinary multiplicative regularization by introducing an extra tuning parameter.

In the course of this paper, the influence of this additional parameter according to the measurement noise level has been studied through a numerical experiment. It has been shown that the value of the extra tuning parameter tends to decrease as the noise level increases. However, because determining an optimal value of this parameter can be tricky, an original iterative procedure has been implemented. The application of this procedure both numerically and experimentally demonstrates the practical interest brought by the generalized multiplicative regularization in terms of solution accuracy and computational efficiency. More peculiarly, it has been shown that completing the iterative process is not necessarily beneficial in terms of solution accuracy and computational efficiency. Consequently, it is recommended to stop the process right after the first iteration. The application of the proposed approach to other forms of regularization terms, such as mixed-norm regularization terms, will be the matter of future works. ( * ) PE computed from the amplitude at the identified point force location more accurate than that obtained for µ = 1 for highly and moderately noisy data, which is in line with the results presented in sections 3.2 and 4.3. In the present case, the regularization term is favored, meaning that solutions based on prior knowledge of the sources to identify are promoted. This explains the disappointing results obtained for highly and extremely noisy data.

  Regińska and Viloche Bazán consisting in introducing an extra tuning parameter in the formulation of the functional to minimize. It results that the target excitation field F is sought as the solution of the following minimization problem:

  exciting, at 350 Hz, a simply-supported steel plate as well as the resulting reaction forces at boundaries. The dimensions of the plate are 0.6 m × 0.4 m × 0.005 m while the external point force is located at point (x 0 , y 0 ) = (0.42 m, 0.25 m), measured form the lower left corner of the plate. To carried out this numerical experiment, a synthesized vibration field and the transfer functions matrix of the structure must be provided as input data to the GMR. The generation of the experimental data is performed in two steps. First, the noiseless reference displacement field X ref is computed from a FE mesh of the plate made up with 187 shell elements 2 , assuming that only bending motions are measured. Then, a Gaussian white noise with a prescribed signal-to-noise Ratio (SNR) is added to the reference data to simulate the measured vibration field. As for the definition of the transfer function matrix H, it has been chosen to compute it from the FE model of the plate with free boundary conditions, assuming, here again, that only bending motions are measurable. Such a modeling enables the reconstruction of the point force and the reaction forces, insofar as they are then considered external to the system.

Figure 1 :

 1 Figure 1: Reference force vector F ref at 350 Hz

3

 3

Figure 2 : 2 (Figure 3 :

 223 Figure 2: Definition of the identification regions -(•) region 1 (Point force), (×) region 2 (Reaction forces) and (♦) location of the point force

Figure 4 :

 4 Figure 4: Reconstructed excitation field at 350 Hz from the GMR for µ = µ (manually tuned) -(a) SNR = 30 dB, (b) SNR = 20 dB, (c) SNR = 10 dB and (d) SNR = 0 dB

Figure 5 :

 5 Figure 5: Plot of ψ(α, µ 0 ) for µ 0 = 1 -(-) ψ (α, µ 0 ), (--) µ = µ 0 and ( * ) α 0 (the abscissa and the ordinate are scaled logarithmically using the common logarithm)

Figure 6 :

 6 Figure 6: Illustration of the first iteration of the process -(a) Selection of µ 1 -zoomed portion of ψ(α, µ 0 ) and (b) plot of ψ(α, µ 1 ) (the abscissa and the ordinate are scaled logarithmically using the common logarithm)

Figure 7 :Figure 8 :

 78 Figure 7: Reconstructed excitation field at 350 Hz from the GMR for µ = µ -(a) SNR = 30 dB, (b) SNR = 20 dB and (c) SNR = 10 dB

Figure 9 :

 9 Figure 9: Experimental set-up -A parallelepiped box excited on one of its faces by a shaker

Figure 10 :

 10 Figure 10: Definition of the identification regions -(•) region 1 (Point force), (×) region 2 (Boundary forces) and (♦) location of the point force

Figure A. 2 :

 2 Figure A.2: Reconstructed excitation field at 350 Hz from the GMR for µ = 2 -(a) SNR = 30 dB, (b) SNR = 20 dB, (c) SNR = 10 dB and (d) SNR = 0 dB

Table 1 :

 1 Performances of the OMR (GMR for µ = 1) for various measurement noise levels -N it : Number of iterations of the resolution algorithm

			SNR
	Indicators 30 dB 20 dB 10 dB 0 dB
	GRE (%)	33.61	30.27	30.74	52.88
	RERF (%) 35.79	31.81	31.10	45.52
	PE (%)	0.83	2.25	0.16	-99.45
	N it	27	13	15	16

Table 2 :

 2 Performances of the GMR for the optimal parameter µ (manually tuned) for various measurement noise levels -N it : Number of iterations of the algorithm

			SNR
	Indicators	30 dB	20 dB	10 dB	0 dB
	µ	5	2.35	1	0.4
	GRE (%)	15.10	20.04	30.74	71.39
	RERF (%)	15.44	20.63	31.10	51.97
	PE (%)	3.3×10 -3 4.2×10 -3	0.16	-2.70 (*)
	N it	12	13	15	51

Table 4 :

 4 Performances of the GMR for µ = µ 1 for various measurement noise levels

				SNR
	Indicators 30 dB 20 dB 10 dB	0 dB
	µ 1	1.53	1.19	0.97	0.50
	GRE (%)	26.37	27.42	31.16	63.69
	RERF (%) 28.15	28.26	31.52	46.17
	PE (%)	0.7	1.92	0.55	-15.40 (*)
	( * ) PE computed from the amplitude at the identified point force location

Table 5 :

 5 Peak error obtained at 625 Hz from the GMR Indicator µ 0 = 1 µ 1 = 3.66 µ = 2.76 × 10 -9

	PE (%)	-1.96	0.089	-7.76

  Table A.1: Performances of the GMR for µ = 0.5 for various measurement noise levels -N it : Number of iterations of the algorithm

				SNR
	Indicators 30 dB 20 dB 10 dB	0 dB
	GRE (%)	56.66	48.58	45.62	63.82
	RERF (%) 55.82	50.11	45.96	46.26
	PE (%)	1.98	3.19	6.91	-15.11 (*)
	N it	31	13	36	10

Table A .

 A 2: Performances of the GMR for µ = 2 for various measurement noise levels -N it : Number of iterations of the algorithm

			SNR
	Indicators 30 dB 20 dB 10 dB 0 dB
	GRE (%)	22.84	21.33	25.46	76.30
	RERF (%) 24.25	22.08	24.68	73.47
	PE (%)	0.59	0.59	-16.96 -99.99
	N it	12	13	16	24
			34

In the present paper, the subscript r (bold or italic) refers to a particular region of the structure. From a notation point of view, a bold character (r) is used for vectors and matrices, while an italic character (r) is used for scalar quantities.

The FE mesh has been defined using 8 elements per bending wavelength. Consequently, the current mesh is theoretically valid up to 590 Hz.

In the present experimental validation, the FE mesh has been defined using 8 elements per bending wavelength. Consequently, the reconstruction mesh is theoretically valid up to 3 kHz.

 4), while the reaction forces are reasonably well estimated. A closer look at Fig. 8d shows that the reconstructed point force is actually located at ( x 0 , y 0 ) = (0.46 m, 0.25 m)

with an amplitude of 0.84 N. Quantitatively, this corresponds to a PE equals to -15.4% for a location error less than 10%, which remains acceptable considering such a low SNR. In particular, it is worth noting that the reconstructed excitation field is less accurate than that obtained from the OMR (i.e. µ = 1) from weakly to highly noisy data. On the contrary, the opposite is observed for extremely noisy data, which is in agreement with the results presented in Figs. 4d and8d. Actually, these conclusions are quite logical, since the data-fidelity term is promoted when µ > 1. This is especially important for an SNR of 0 dB, since the information brought by the reconstruction model (see Eq. ( 4)) helps to obtain a physically consistent solution.