Juan Carlos Arceo

Jimmy Lauber

Emilie Simoneau

Sylvain Cremoux

Nonlinear Convex Control for Robotic Assistive Therapy via LMIs

A device for personalized robotic assisted therapy is introduced, the process for system identification and control design for a class of nonlinear discrete-time systems with output disturbance attenuation is presented. The results obtained are in terms of linear matrix inequalities. Simulation and real-time results are shown as well.

I. INTRODUCTION

According to World Health Organization there are 15 million stroke cases each year, 5 million of these people will be disabled. The economic cost for this illness is equivalent to 53.6 billion dollars in USA, more than 4% of the total health care cost in UK and 3% in Netherlands [START_REF] Mackay | The atlas of heart disease and stroke[END_REF].

One of the goals in robotic rehabiliation is restore the motor function of stroke patients and this can be achieved by constraint-induced movement. Therapy must adapt to user's requirements [START_REF] Riener | Humancentered rehabilitation robotics[END_REF] in aspects such as the degree of disability and its economic cost [START_REF] Díaz | Lower-limb robotic rehabilitation: literature review and challenges[END_REF], given this context and inspired by [START_REF] Xu | A closed-loop brain-computer interface triggering an active ankle-foot orthosis for inducing cortical neural plasticity[END_REF] the motoBOTTE shown in Fig. 1 is a device designed for robotic assisted therapy with a low production cost.

This system will be coupled with a hybrid brain-computer interface and a force sensor to offer personalized therapy and a first step is track a reference with the ankle (denoted by φ). Therefore, this report deals with modeling and control of the motoBOTTE. Convex systems (also named Takagi-Sugeno) [START_REF] Taniguchi | Model construction, rule reduction and robust compensation for generalized form of Takagi-Sugeno fuzzy systems[END_REF] combined with Lyapunov's direct method can be used for nonlinear controller design; a common approach is parallel compensated distribution (PDC) [START_REF] Tanaka | Fuzzy Control Systems Design and Analysis: A linear matrix inequality approach[END_REF], which mimics the nonlinearities in the system by means of the interpolation functions, this approach yields to conditions in terms of linear matrix inequalities (LMIs) that can be solved via optimization algorithms in polynomial time with commercially available software [START_REF] Lofberg | YALMIP : a toolbox for modeling and optimization in matlab[END_REF]; performance parameters like disturbance attenuation, input constraint and decay rate can easily be included in the design [START_REF] Tanaka | Fuzzy Control Systems Design and Analysis: A linear matrix inequality approach[END_REF].

In Section II system identification and nonlinear controller design processes are shown. In Section III simulation and real-time results are exhibited. In Section IV a conclusion based on the results is given.

II. PRELIMINARIES

The system is actuated by a linear piston [START_REF]CAHB-21 Linear Actuator: Installation, operation and maintenance manual[END_REF], the angle is limited mechanically φ ∈ [1.273, 2.5212] in radians; the sampling time (denoted by T s ) is 0.01 seconds, there is an *This work was supported by the University of Valenciennes. 1 The authors are with the Univ. Valenciennes, LAMIH UMR CNRS 8201, F-59313 Valenciennes, France. (corresponding author:

Jimmy.Lauber@univ-valenciennes.fr) input constraint of ±10 Volts and if the input equals zero the system will not move. Based on the given information we can describe the dynamics of the motoBOTTE with: Defining the filtered reference as r f (t + 1) = Ar f (t) + (1 -A)r(t), with |A| < 1, r(t) is the reference, the error e(t + 1) = φ(t + 1) -r f (t + 1), its integral e i (t + 1) = e i (t)+T s e(t) and its time derivative e d (t+1) = e(t+1)

φ(t+1) = φ(t)+(b 0 +b 1 φ(t)) u(t), (1)
Ts -e(t) Ts ; using ∆ 1 = r r f and ∆ 2 = r f φ to gather the expressions as Ex(t + 1) = A(∆ 1 )x(t) + B(φ)u(t) + P (∆ 1 , ∆ 2 , φ) with real estimated Fig. 2. Estimation process.

x(t) = e(t) e i (t) e d (t)

T . If r f = r, then ∆ 1 = 1 and P (•) = 0, we will design a controller for this. As there exists E -1 , our error system is x(t+1) = Āx(t)+ B(φ)u(t), with

E =   1 0 0 0 1 0 -1 0 T s   , A =   1 0 0 T s 1 0 -1 0 0   , B =   b 0 +b 1 φ(t) 0 0   , Ā=E -1 A and B(φ) = E -1 B.
The nonconstant term φ(t) can be rewritten by means of the nonlinear sector approach [START_REF] Taniguchi | Model construction, rule reduction and robust compensation for generalized form of Takagi-Sugeno fuzzy systems[END_REF], it is bounded by φ ∈ [1.273, 2.5212]. Define the interpolation functions as h 1 = φmax-φ(t) φmax-φmin and h 2 = φ(t)-φmin φmax-φmin , notice that h 1 +h 2 = 1 and 0 ≤ h i ≤ 1 for i ∈ {1, 2}, the convex system is x(t+1) = For control design with disturbance attenuation consider a disturbance w(t) measured in the output y(t) = Cx(t) + Dw(t). The origin of the system is asymptotically stable with a PDC control law [START_REF] Tanaka | Fuzzy Control Systems Design and Analysis: A linear matrix inequality approach[END_REF]

u(t) = 2 j=1 F j x(t) = F h x(t),
decay rate α 2 and γ is minimized subject to y(t) w(t) ≤ γ if there are matrices X > 0, γ > 0, M j for (i, j) ∈ {1, 2}, where F i = M i X -1 and 2Υ ii +Υ ij +Υ ji < 0 [START_REF] Tuan | Parameterized linear matrix inequality techniques in fuzzy control system design[END_REF] holds with :

Υ ij =     -α 2 X XC T D XA T i +M T j B T i XC T D T CX D T D-γ 2 I 0 0 A i X+B i M j 0 -X 0 CX 0 0 -I     .

III. RESULTS

The output and disturbance matrices are C = I 3 and D = 0.1 × 1 0.01 100 T , we want a decay rate of α 2 = 0.935. A solution was found for the LMIs with γ = 10.0241, F 1 = 146.7459 0.1757 and F 2 = 285.8307 1127.1404 0.3419 . The closed-loop response is shown in Fig. 3, where the dashed line is the filtered reference, the solid black line is the output angle in simulation, the blue line is the real-time output angle and the red line is the real-time output angle where the disturbance is a leg pushing in an opposite direction to r f . The control signal used for each case is illustrated in Fig. 4.

IV. CONCLUSIONS

The mathematical model of a device used for robotic assisted therapy has been obtained, then it was used for nonlinear controller design with output disturbance attenuation guaranteed. Conditions in terms of LMI have been provided and results in simulation and real-time have been shown. The performance of the closed-loop system might be improved by combining a feedforward with a feedback control approach considering that the trajectory its known a priori.

Fig. 1 .

 1 Fig. 1. The motoBOTTE.

where b 0

 0 and b 1 are unknown parameters to be estimated via input and output data and the Parameter Estimation Toolbox in MATLAB. The real-time signals (in black) and the estimated (in blue) are shown in Fig. 2, a squared error cost function was optimized to a value of ss cf = 1.0475 via the Levenberg-Marquardt algorithm [4] the parameters values obtained were b 0 = -0.00203 and b 1 = 0.00054.

2 i=1

 2 h i Āx(t)+ Bi u(t) , with Bi = B(φ)| hi=1 constant terms B1 =

Fig. 3 .Fig. 4 .

 34 Fig. 3. Simulation and real-time response.


Pablo J Pavan
email: pjpavan@inf.ufrgs.br

Matheus S Serpa
email: msserpa@inf.ufrgs.br

Emmanuell D Carreño
email: edcarreno@inf.ufpr.br

Víctor Martínez
email: victor.martinez@inf.ufrgs.br

Edson L Padoin
email: padoin@unijui.edu.br

Philippe O A Navaux
email: navaux@inf.ufrgs.br

Jairo Panetta
email: jairo.panetta@gmail.com

Jean-François Mehaut
email: jean-francois.mehaut@imag.fr

Improving Performance and Energy Efficiency of Geophysics Applications on GPU Architectures

Keywords: Geophysics applications, Manycore systems, Energy efficiency, GPU

Energy and performance of parallel systems are an increasing concern for new large-scale systems. Research has been developed in response to this challenge aiming the manufacture of more energy efficient systems. In this context, this paper proposes optimization methods to accelerate performance and increase energy efficiency of geophysics applications used in conjunction to algorithm and GPU memory characteristics. The optimizations we developed applied to Graphics Processing Units (GPU) algorithms for stencil applications achieve a performance improvement of up to 44.65% compared with the read-only version. The computational results have shown that the combination of use read-only memory, the Z-axis internalization and reuse of specific architecture registers allow increase the energy efficiency of up to 54.11% when shared memory was used and increase of up to 44.53% when read-only was used.

Introduction

Several applications in areas, such as physics simulation, weather forecast, oil exploration, climate modeling and atomic simulation require high processing power and efficient models. Some of these scientific applications make use of stencil computations that include both implicit and explicit Partial Differential Equations (PDE) solvers [START_REF] Datta | Stencil computation optimization and auto-tuning on state-of-the-art multicore architectures[END_REF]. Besides the scientific importance of stencils, they are interesting as an architectural evaluation benchmark because they have abundant parallelism and low computational intensity, offering opportunities for on-chip parallelism and challenges for associated memory systems [START_REF] Datta | Stencil computation optimization and auto-tuning on state-of-the-art multicore architectures[END_REF]. Today, PetaFlops systems allow reaching increasingly accurate results these scientific applications. To respond to the high processing demand of stencil applications, High Performance Computing (HPC) systems gather the processing power of several computational resources to solve these problems.

Scientific simulations may consume weeks of supercomputer time and most of this time is spent in stencil computations [START_REF] De La Cruz | Towards a multi-level cache performance model for 3d stencil computation[END_REF]. Continuous changes in the fabrication process of the microprocessors industry have increased the performance of its products and influenced state-of-the-art HPC systems. However, this exponential increase in computational performance also leads to an exponential growth in power demand [START_REF] Dong | Power measurements and analyses of massive object storage system[END_REF][START_REF] Laros | Topics on measuring real power usage on high performance computing platforms[END_REF][START_REF] Padoin | ARMbased cluster: Performance, Scalability and Energy Efficiency[END_REF]. Reductions in the total execution time of applications are also relevant for energy consumption, energy is saved when hardware resources are used for a shorter time.

However, it is possible to achieve even greater energy savings if the application is able to exploit the different memory levels available. Today, the combined use of Graphics Processing Units (GPUs) and CPUs in HPC systems has become a popular choice among the top ranked and yet to come platforms. Stencil computing is typically memory-bound, memory performance is particularly important for most stencil kernels. GPUs have several processing elements inside a single die and different memory levels. For this reason, one of the most important strategies for optimizing the performance of stencil computing is the optimization of memory access. Besides, stencil computing can be ported to GPUs with significantly improved performance when compared to implementations performed on CPUs [START_REF] Maruyama | Optimizing stencil computations for nvidia kepler gpus[END_REF].

The performance of a stencil for a given architecture can be estimated through the roofline model [START_REF] Williams | Roofline: An insightful visual performance model for multicore architectures[END_REF]. This model relates the maximum performance of code to its computational intensity, considering the speed of memory access and the processing capacity of the machine [START_REF] Nasciutti | Impacto da arquitetura de memria de gpgpus na velocidade de computao de estnceis[END_REF]. In this paper, we improved the performance and achieved increase energy efficiency of stencil applications by improving methods and optimizations of GPU code. Those are used in conjunction with specific GPU memory characteristics. We focus on analyzing the impact of the stencil size and usage of different memory hierarchies and registers of the GPU to improve performance, power demand, energy consumption and energy efficiency.

The remaining sections of this paper are organized as follows. Section 2 discusses some of the related works on energy consumption In Section 3, we present the stencils application and details of our versions and optimization developed.

In Section 4, we present the evaluation methodology used in the conducted experiments and the stencils and their implementation details. In addition, in Section 5, we address the results obtained from the experiments. Finally, the Section 6 emphasizes the scientific contribution of the work and notes several challenges that we can address in the future.

Related Work

Several studies have evaluated performance of stencils to improve their energy efficiency in CPUs and GPUs systems. Despite that, the processors and accelerators remain as the component with the highest power demand of the systems [6]. GPUs are made aiming at massively parallel processing, to achieve this they use hundreds of processing units working together. These characteristics lead to its superior energy efficiency if compared with CPUs systems [START_REF] Padoin | Evaluating application performance and energy consumption on hybrid cpu+ gpu architecture[END_REF].

Micikevicius et al. [START_REF] Micikevicius | 3d finite difference computation on gpus using cuda[END_REF] compared the performance of a stencil ported from CPU to GPU. Their version of the stencil running in a GPU achieved an order of magnitude higher than running in a contemporary CPU. They conclude that it is possible to improve their results by the usage of shared memory to reduce communication overhead.

Bauer et al. [START_REF] Bauer | Cudadma: Optimizing gpu memory bandwidth via warp specialization[END_REF] showed that the main bottleneck in GPU applications are related to the memory system. To reduce its impact, they used DMA warps to improve memory transfer between on-chip and off-chip memories. They achieved a speedup up to 3.2 times on several kernels form scientific applications.

Schäfer and Fey [START_REF] Schafer | High performance stencil code algorithms for gpgpus[END_REF] evaluate a set of algorithms on Fermi GPUs. They evaluate micro benchmarks using shared memory and found that using only L1 cache creates a problem for its limited throughput. Also, the L2 cache is not a good option because of cache blocking. They conclude that a new alternative to use shared memory was needed to overcome communication bottleneck.

Falch and Elster [START_REF] Falch | Register caching for stencil computations on gpus[END_REF] proposed the usage of a manually managed cache to combine the memory from multiple threads. Using their technique, they achieved a speedup of up to 2.04 in a synthetic stencil. They concluded that manual caching is an effective approach to improve memory access and that applications with regular access patterns are suitable to implement their technique.

Zhou et al. [START_REF] Zhou | A novel gpuaccelerated strategy for contingency screening of static security analysis[END_REF] points that the use of GPUs enables considerable gains in performance compared to using CPU. They have applied GPUs successfully in many computations and memory intensive realms due to its superior performances in float-pointing calculation, memory bandwidth, and power consumption. The results obtained show a speedup of up to 50 times using GPU algorithm rather than CPU algorithm. In similar works, Zhou et al. [START_REF] Zhou | Hands-on performance tuning of 3d finite difference earthquake simulation on gpu fermi chipset[END_REF] obtained a speedup between 10 and 15 times using a GPU rather than CPU.

Xue et al. [START_REF] Xue | An efficient gpu implementation for locating micro-seismic sources using 3d elastic wave time-reversal imaging[END_REF] also make comparisons between GPU and CPU implementation. They obtained a speedup up to 18 times in the GPU-based implementation of a time-reversal imaging micro-seismic event location. Also, Nikitin et al. [START_REF] Nikitin | Parallel algorithm of 3d wave-packet decomposition of seismic data: Implementation and optimization for gpu[END_REF] obtained average speedup up to 46 times using GPU for compared to CPU for processing a synthetic seismic data set (data compression, de-noising, and interpolation).

Maruyama and Aoky [START_REF] Maruyama | Optimizing stencil computations for nvidia kepler gpus[END_REF] presents a method for stencil computations on the NVIDIA Kepler architecture that uses shared memory for better data locality combined with warp specialization for higher instruction throughput, their method achieves approximately 80% of the value from roof line model estimation.

Hamilton et al. [START_REF] Hamilton | Large stencil operations for gpubased 3-d acoustics simulations[END_REF] investigate the computational performance of GPU-based stencil operations using stencils of varying shape and size (ranging from seven to more than 450 points in size). They found that using an NVIDIA K20 GPU, data movement, rather than computing, was the bottleneck, and as such, the performance obtained can be attributed to the effects of the L2 and texture caches on the card.

Compact stencils are more efficient using the texture cache and require fewer reads from global memory. The leggy stencils schemes required a significant portion of global memory bandwidth in order to achieve similar performance as compact stencils of similar size in points.

Nasciutti and Panetta [START_REF] Nasciutti | Impacto da arquitetura de memria de gpgpus na velocidade de computao de estnceis[END_REF] did a performance analysis of 3D stencils on GPUs focusing on the proper use of the memory hierarchy. They conclude that the preferred code is the combination of read only cache reuse, inserting the Z loop into the kernel and register reuse.

Different to other approaches that allocate workload on CPU and GPU architectures, or works that use GPUs to achieve considerable performance gains when compared to traditional CPU architecture, our goal aims to increase the performance and energy efficiency of stencil application applying methods and optimization to use different memory levels of the GPUs.

Geophysical Model Optmizations

The model simulates the collection of data in a seismic wave propagation. At intervals of, equipment coupled to the ship emits waves that reflect and refract on changes of the medium in the subsoil. Eventually, these waves return to the surface of the sea, being collected by specific microphones (geophones) coupled to cables towed by the ship. The set of signals received by each geophone over time constitutes a seismic trace. For each wave emission, the seismic traces of all cable geophones are recorded. The ship continues to sailing and emits signals over time.

Acoustic wave propagation approximation is the current backbone for seismic imaging tools. It has been extensively applied to imaging potential oil and gas reservoirs beneath salt domes. We consider the model formulated by the isotropic acoustic wave propagation under Dirichlet boundary conditions over a finite 3D rectangular domain, prescribing to all boundaries, and the isotropic acoustic wave propagation. Propagation speed depends on variable density, the acoustic pressure, and the media density. These applications are modeling and solved using stencil computations. In this context, the computational performance of GPU-based stencils have a great scientific importance as it is used in many areas of scientific computing. Regarding the capatiblities of current GPU architectures, the NVIDIA Kepler provide memories with different characteristics compared with CPUs. One of the main differences between GPUs and CPUs is the way their memory subsystem work. In a CPU, access to memory is done by obtaining their data from caches. Usually looking on L1, L2, L3, and DRAM in that order. On the other hand, in a GPU the L1 memory cache, is used specifically for accesses to the stack and register spill, i.e., when too many local variables do not fit in the register file, and thus some of it has to be cached. L2 memory is used for global accesses requested by stream processors.

The current GPU have also registers files, a shared memory. They are a texture memory and a global memory with different characteristics such as size, speed, read-only memory and in the way that is possible to use them. These registers were not available in Nvidia GPUs before Kepler architecture. In Figure 1 is shown an overview of the Kepler GPU architecture, which have different SP (Stream Processor) in each SMX (Streaming Multiprocessors).

To exploit the use of different memory levels available on current GPU, we develop three versions of a stencil kernel using each one of the GPU memories. Each stencil version, give us a different insight of the performance and capabilities of the GPU memory subsystem.

-The first version called naive take no advantage of any of the GPU highspeed memories and access data only from global memory. -The second version called shared stores one part of the stencil data in the shared memory scratchpad. The shared-memory version also uses the GPU resources that the naive version uses, the main difference is that this version also uses the shared-memory available on each SMX (Streaming Multiprocessors). Each one of the SMX have one internal shared-memory to store data as shown in Figure 1. In this version, data is manually allocated by the programmer through the use of the shared directive, indicating such data will be shared among all the GPU threads. The compiler automatically configures the space division between the L1 cache memory and the shared cache memory, choosing one of three options: 16 KB for the L1 cache and 48 KB for the shared cache, 32 KB for each, or 48 KB for the L1 cache and 16 KB for the shared cache. -The third version called read-only stores most read data in a read-only texture memory which is faster than shared memory but works with read-only data. This version takes advantage of the read-only cache, this cache is the SMX memory bank that stores only read data, it is also called texture memory. Originally it was used only for textures, but starting with the Kepler architecture any data can be stored in this cache by using the C-99 directive const restrict. The programmer may also explicitly use this cache through the intrinsic lgd().

We developed two optimizations for each of the versions to evaluate improvements in performance and energy efficiency by reusing the Z direction data. Reusing Z direction data is named internalization.

-The int.z version takes advantage of data locality by storing stencil data for direction Z. This optimization consists of the internalization of the Zaxis into the threads. Doing the internalization ensures that neighbouring Z-blocks execute sequentially, increasing the reuse of L2 cache data. Direction Z data is used to calculate subsequents points in the X-Y direction. -The int.z.reg version consists of combining the int.z with the usage of registers to store the Z direction points. For example, to calculate the point Z3 in a 13 points stencil, the neighbouring points in X and Y, as well as points Z1, Z2, Z3, Z4 and Z5 are required. In order to calculate the points in Z4, points Z2, Z3, Z4 and Z5 would be availed, and it is necessary to request the global memory only points Z6, as well as the neighbours in X and Y.

Experimental Methodology

Our experiments were developed in a NVIDIA K20m GPU card. The K20m is a Kepler architecture GPU with 2496 CUDA cores. Each Streaming Multiprocessor has a configurable on-chip memory that can be configured as 48/32/16 KB shared memory with 16/32/48 KB of L1 cache. They also have a faster 48 KB read-only cache and a 1280 KB shared L2 cache. Table 1 describes in detail the environment we used. We used NVIDIA Management Library (NVML) to measure the power usage. Regarding the energy efficiency measurement, we used the metric of performance achieved divided by average power. Each experiment was executed 30 times, we show average values, as well a 95% confidence interval calculated with Student's t-distribution.

Results

This section shows the optimizations techniques we used to improve the performance and energy efficiency of a stencil application. The stencil we used In the following subsections, we describe each optimization and analyze how they address the performance and energy efficiency improvements. We also show the results obtained by using the three different memories and the results of the optimizations applied in each of them, on a NVIDIA Kepler architecture.

Performance and Energy Efficiency improvements over Naive version

This subsection shows the improvements obtained by using two optimization techniques over the naive version of the stencil computation. Figure 2 shows the performance and energy efficiency of the naive version and the optimizations. The first optimization technique used was the int.z which stores data from direction Z in local variables aiming to take advantage of the data locality by reusing these data in the subsequent iterations. The performance and energy efficient were improved by up to 4.65% and up to 4.55% using the int.z technique over a naive version. This improvement occurs due to the reuse of L2 cache data made by this optimization. The number of access in global memory is reduced by increasing L2 cache hits.

We propose a second optimization called int.z.reg which consists of the int.z optimization along with the use of the register file to store the Z points. In int.z Z points were stored only in local variables. Using this optimization performance overtakes the previous versions with an improvement of up to 34.31% compared with the naive version. The energy efficiency was also improved by up to 34.30%. The results show that the usage of registers, which are faster than local variables, allow us to obtain more performance with a better energy efficency.

Performance and Energy Efficiency improvements over Shared Memory

In the previous subsection, we showed the optimizations applied in the naive version. Although the performance and energy efficiency was improved by both optimizations techniques, the naive version does not take advantage of fast GPU memories as shared memory. Thus, we improved the naive version by using the shared memory scratch pad to store a slice of data that is reused by the threads of the same block. The data was manually allocated using the shared directive, indicating a piece of data shared among all threads. We also applied the int.z and int.z.reg optimizations aiming to improve the performance of the memory operations.

The performance and energy efficiency results are showed in Figure 3. Using this optimization, performance was improved by up to 2.25% and 54.46% in the int.z and int.z.reg optimizations compared with the shared memory version. It occurs due to the data stored in scratch pad is reused by the threads in the following iterations. The energy efficiency was improved by up to 2.02% and 54.11% using these optimizations.

Performance and Energy Efficiency improvements over Read-only Memory

In this subsection, we are showing the improvements obtained when we use both optimizations and the read-only memory. Since the data we store in the shared memory was not update we may take advantage of the read-only memory. The read-only memory is faster than shared memory but exclusively used for readonly operations. We can explicitly define that global memory reads be stored in the read-only memory using the lgd() intrinsic. The int.z optimization over the read.only version achieve a performance improvement of up to 34.30%. Implementing the int.z.reg that also uses the register file the performance was improved by up to 44.65%. The energy efficiency was also improved by these optimizations. The int.z version improved the energy efficiency by up to 34.20% while the int.z.reg improved the efficiency by up to 44.53%.

Conclusion

Several scientific applications make use of stencil computations to their model simulations. Stencils have both implicit and explicit PDE being so also interesting as an architectural evaluation benchmark. The computing present in these applications are low intensity, once that they are typically memory-bound. In this form, memory optimizations are important for to use the fastest memories available in GPUs and increase their the energy efficiency.

In this paper, aim to achieve energy savings, we introduce methods and optimization to stencil application that exploit the different memory levels available. Our developed methods, which are used in conjunction with specific GPU memory characteristics, allow to use the read-only cache and also the shared memory. Also, our developed optimization allows to combine the Z-axis internalization of stencil application with the reuse of registers of GPU architecture.

The main contribution of this paper is performance and energy efficiency increases when applied GPU-algorithms and optimization over stencil application. Our developed GPU-optimized algorithms for stencil applications achieve performance improvement of up to 54.11% and 44.53% when were used shared memory and read-only cache respectively over the naive version. This increase in computational performance also improves the energy efficiency in an equivalent value, once that our methods and optimization do not increase the power demand.

Changes in the GPU architecture, as in the case of the introduction of the read-only cache in the Kepler architecture, can generate changes in the results presented in this work. In the future, we plan to investigate methods and optimization to achieve gains in stencil applications over new NVIDIA architecture and Intel Xeon Phi.

Fig. 1 :

 1 Fig. 1: Sample of the memory subsystem on a NVIDIA Kepler architecture

Fig. 2 :

 2 Fig. 2: Improvements over Naive version which uses Global Memory.

Fig. 3 :

 3 Fig. 3: Improvements over Shared Memory.

  read.only read.only.int.z read.only.int.z.

Fig. 4 :

 4 Fig. 4: Improvements over Read-only Memory.

Table 1 :

 1 Configuration of GPU system. simulates the propagation of a single wavelet over time. To create the simulation, it solves the isotropic acoustic wave propagation with constant density under Dirichlet boundary conditions over a 3D domain. The stencil is a 13-arm with the following input sizes: (1024 × 256 × 256), (2048 × 256 × 256), (4096 × 256 × 256), and (7168 × 256 × 256).

Acknowledgments

This research has received funding from the EU H2020 Programme and from MCTI/RNP-Brazil under the HPC4E Project, grant agreement n. o 689772. It was also supported by Intel under the Modern Code project, and the PETRO-BRAS oil company under Ref. 2016/00133-9. We also thank to RICAP, partially funded by the Ibero-American Program of Science and Technology for Development (CYTED), Ref. 517RT0529.


Julien Mandon

Evgeni Sorokin

Irina T Sorokina

Guy Guelachvili

Dr Nathalie Picqué
email: nathalie.picque@u-psud.fr

Supercontinua for high

Keywords: optics, fibers, 060.2390 Fiber optics, infrared, 140.7090 Ultrafast lasers, 230.6080 Optical devices, Sources, 300.6390 Spectroscopy, molecular, 140.3070 Infrared and far-infrared lasers

come

Bright, coherent and broadband radiation brings a number of advantages for high resolution spectroscopy, so that a number of several facilities are, for instance, presently using [1] -or planning to use -synchrotron radiation for gas phase spectroscopy. As a much simpler tabletop low-cost alternative, femtosecond (fs) mode-locked lasers represent the first versatile coherent sources providing broad bandwidth and high brightness. Several spectroscopic methods using fs mode-locked lasers are currently being implemented. Most of them make use of multichannel grating spectrometers [START_REF] Diddams | Molecular fingerprinting with the resolved modes of a femtosecond laser frequency comb[END_REF]3], which limit the resolution and applicable spectral domains. Multiplex spectroscopy does not meet these restrictions and various approaches have been successfully demonstrated [4][5][6], in particular in the infrared where detector arrays for multichannel recording are expensive or inexistent. These approaches are rather easy and cheap to be put in practice with commercial and laboratory instruments widely used in chemistry and physics. However, their full potential is currently limited by the narrow spectral range of available fs sources in the infrared. The supercontinuum (SC) generation may be a convenient and inexpensive way to overcome this restriction.

SC can be conveniently generated from femtosecond lasers or any high energy pulsed lasers when launched into nonlinear optical fibers. In particular, Photonic Crystal Fibers (PCF) have permitted more than one octave spectral broadening in the near-infrared and visible [7]. Currently, SC is boosting several fields like frequency metrology, pulse compression and optical coherence tomography. Its application to optical gas monitoring is presently restricted to wavelength agile methods [8][9][10]. Multiplex Fourier Transform (FT) spectroscopy is an interesting opportunity to benefit from these broad sources for high accuracy and high sensitivity gas phase spectroscopy. The aim of this paper is to show that SC may be fully taken advantage of in multiplex spectroscopy and to illustrate the SC potential for high resolution FT spectroscopy on the example of a SC generated in a highly nonlinear fiber coupled to a 1.5 µm fs laser.

The experimental setup is shown in Fig 1 . A chirped-mirror-controlled Cr 4+ :YAG oscillator [11], pumped by a Nd:YVO 4 laser, is modelocked by a commercial SESAM (BATOP Mod. SAM-1550-1). It is operated in the chirped-pulse positive dispersion regime [12]. In this regime, details of which are described elsewhere [13], the laser emits positively chirped 1.4 ps pulses centered at 1525 nm with a repetition rate of 142 MHz and 100 mW average power. These pulses are compressed down to 120 fs in a 3.2 m piece of a standard SMF28 single mode silica fiber, and launched into the nonlinear fiber with about 550 pJ input energy. We have tested several kinds of nonlinear fibers, including dispersion-shifted germanium-doped silica fibers (DSF) and PCFs. The most effective broadening has been obtained with a 40-cm long PCF made of SF6 glass [14] in which the light, launched with an efficiency of about 32%, is guided through a 4.5 µm-diameter solid core with nonlinear index n 2 = 2.2×10 -15 cm 2 .W -1 , surrounded by a periodic arrangement of air holes. After recollimation, the SC beam has about 25 mW average power and passes through a single-pass 70-cm-long absorption sample cell. The beam is attenuated and analyzed by a commercial rapid-scan FT spectrometer (Bruker IFS66) just as a steady-state source.

The achieved spectral broadening is shown on Fig. 2 at low resolution (50 cm -1 -ca. 10 nm). The spectral domain of the laser emission itself (without broadening) extends over 200 nm at -30 dB level. In the high resolution spectra reported below, this 30 dB attenuation corresponds to the limit beyond which no spectroscopic data are measurable. At the PCF output, the spectrum extends in the infrared over about 1100 nm, from 1200 (8330 cm -1 ) to 2200 nm (4545 cm -1 ) at -30 dB level, and in the visible-very near infrared over about 275 nm, from 675 (14815 cm -1 ) to 950 nm (10525 cm -1 ) at -30 dB level. This broadening is consistent with previous experimental data [ [14][15][16][17] on the same type of fiber with different core sizes and various seeding sources.

Each high resolution (3.6 GHz -0.12 cm -1 ) spectrum of the PCF output, is recorded in in just 13.1 s. It is made of 67290 independent useful spectral elements covering 1375 nm (8075 cm -1 ), with a maximum signal to noise ratio (SNR) about 2400. At this resolution, fringes scrambling spectral information are observed. Their origin is not completely clear: they do not seem to be due to birefringence from the PCF as rotating the input light polarization has no effect. They are however specific to the employed PCF, because the SC generated in a single mode DSF is fringe-free. We believe that these fringes are due to the mode beating in the multimode PCF that could not be observed previously due to insufficient resolution [16]. The overall short-term stability made possible the efficient suppression of the fringes, by dividing the spectra with and without gas inside the absorption cell. Since the broadening was much stronger in the multimode PCF, we decided to use it to test the SC potential for high resolution spectroscopy, trading the slightly degraded SNR against the larger bandwidth.

As a first demonstration, C 2 H 2 and NH 3 spectra are recorded. Figure 3 shows a portion of the acetylene absorption spectrum at two different pressures: 28 and 160 hPa, between 1511 and 1548 nm (6620-6460 cm -1 ) at 3.6 GHz (0.12 cm -1 ) resolution. All absorption lines displayed in Fig. 3 are due to acetylene. The most intense band is the ν 1 +ν 3 band of 12 C 2 H 2 .

The rotational lines from the P and R branches of the ν 1 +ν 3 +ν 4 1 -ν 4 1 and ν 1 +ν 3 +ν 5 1 -ν 5 1 hot bands of 12 C 2 H 2 and the ν 1 +ν 3 cold band of 12 C 13 CH 2 are also observed. The ammonia spectrum is shown on the lower trace of Fig. 4 at 15 GHz (0.50 cm -1 ) resolution. The region of the ν 2 +ν 3 +ν 4 band of NH 3 is observed. The spectrum is complex and crowded and to our knowledge, no rotational assignment is available in the literature. For comparison, the upper trace provides a portion of spectrum recorded with a tungsten lamp, which is the most widely used white source in near-infrared absorption FTS. The SNR of this reference spectrum is about 1.7 times better than the SC spectrum, but the recording time is 200 times longer. With the same recording time as with the tungsten lamp, the SNR of the SC spectrum would have been 33 000. One advantage of SC sources is their high spectral radiance, which for the present 25 mW SC source equals 2.8×10 7 W.m -2 .sr -1 .nm -1 . It is 1×10 5 times stronger than a tungsten lamp comparable to a 3000 K black body source. A SNR enhancement of about 350 could have then been achieved. This was actually not possible: the detector was saturated by the 25 mW laser beam and a power attenuation by more than one order of magnitude before entering in the interferometer was necessary. On the other hand, tungsten lamps have a flat spectrum, which enables to benefit from the whole spectrum. Further optimization of SC generation providing flat spectra without strong intensity variations would improve their applicability for spectroscopic measurements.

It is also informative to compare the potential of SC sources with the performance of synchrotron radiation. The expected brightness [18] of the infrared beamline at the French facility Soleil is 10 19 photons/s/1% of BW/mm 2 .steradian for the near infrared (1 µm), which converts into a spectral radiance equal to 2×10 5 W.m -2 .sr -1 .nm -1 . This is two orders of magnitude weaker than the present SC, which can still be improved in terms of output power and spectral coverage. Moreover, since such SC source is a table-top device made of commercially available parts, it could even be incorporated into industrial instruments. Presently the wavelength agile techniques are the fastest gas monitoring methods using the SC sources at comparable resolution [8][9][10]. They, however, feature frequency-dependent spectral resolution and require long fibers as dispersers, that currently limit their applicability to the telecom bands and have wavelength sweeps that may be too fast. A FT spectrometer requires longer acquisition times, seconds instead of milliseconds, but is much more appropriate for accurate measurements of spectroscopic interest. The FT spectral resolution may be improved easily up to 30 MHz, the accuracy in the wavenumber scale may reach 10 -9 , and several millions of spectral elements may be measured in a single experiment, resulting in spectra spanning more than one octave. In addition, SC can be implemented with the frequency comb spectrometric approach that we reported recently with a mode-locked laser [6], which enhances sensitivity and enables simultaneous recording of broadband absorption and dispersion spectra.

In this article, supercontinuum femtosecond laser systems are demonstrated to be powerful high-brightness sources for high resolution spectroscopy. The corresponding accurate analysis may be well performed by FT spectrometers over the whole emission range. SC sources become commercially available in the visible and near-infrared ranges and they may be easily implemented on the basis of existing spectroscopic instruments. The development of a tabletop ultrabroadband coherent source in the mid/far-infrared should significantly impact the practical applications of spectroscopy, and SC generated by chalcogenide or fluoride fibers [19] could likely be an efficient approach.

O. Okhotnikov (Tampere University of Technology) is warmly thanked for providing the samples of nonlinear Ge-doped fiber, P. Jacquet (Laboratoire de Photophysique Moléculaire, Orsay) for his participation and M. Hanna (Laboratoire Charles Fabry de l'Institut d'Optique) for useful advices. Efficient technical support was provided by J. Caillaud. I.T Sorokina acknowledges an Invited Professor position with Université Paris-Sud. The Programme Pluri-Formation de l'Université Paris-Sud "Détection de traces de gaz" and the Austrian Fonds zur Förderung der wissenschaftlichen Forschung (project P17973) have funded this work.

Figure captions

Fig. 1 .

 1 Fig. 1. Experimental setup. The chirped pulses generated by the Cr 4+ :YAG oscillator are compressed by a single mode fiber and spectrally broadened in the PCF. The SC light probes a cell and is analyzed by a FT spectrometer. CM: chirped mirrors, OC: output coupler.

Fig. 2 .

 2 Fig.2. Low resolution spectra on a logarithmic vertical scale of the laser output (blue line, plot filled in white) and the 4.5 µm-core-diameter SF6 PCF output (red line, plot filled in magenta).

Fig. 3 .Fig. 4 .

 34 Fig.3. SC high resolution spectra on a linear vertical scale of the acetylene molecule in the 1.5 µm region at 28 and 160 hPa. The strongest spectral features are the P and R branches of the ν 1 +ν 3 band of12 C 2 H 2 .

Fig. 1 .Fig. 2 .

 12 Fig. 1.Experimental setup. The chirped pulses generated by the Cr 4+ :YAG oscillator are compressed by a single mode fiber and spectrally broadened in the PCF. The SC light probes a cell and is analyzed by a FT spectrometer. CM: chirped mirrors, OC: output coupler.

Fig. 3 . 3 Fig. 4 .

 334 Fig.3.SC high resolution spectra on a linear vertical scale of the acetylene molecule in the 1.5 µm region at 28 and 160 hPa. The strongest spectral features are the P and R branches of the ν 1 +ν 3 band of12 C 2 H 2 .


Bastien Mussard
email: bastien.mussard@colorado.edu

Random Phase Approximation in Projected Oscillator Orbitals

The projected oscillator orbitals (pOOs) are localized virtual orbitals constructed by multiplying localized occupied orbitals by harmonics. Following a recent paper by Mussard and Ángyán 1 , further developments of projected oscillator orbitals are shown, notably the equations for pOOs of general order as well as their overlaps are derived. The performance of these localized virtual orbitals is demonstrated up to third order. It is found that a good fraction of the aug-cc-pVQZ RPA correlation energy is recovered by use of a smaller number of pOOs. This is especially true where considering only the long-range correlation energy, which is important for the description of London dispersion forces.

INTRODUCTION

In the context of density functional approximations, the functionals from the so-called fifth rung of the Jacob's ladder not only use the density, its gradient, and a set of occupied orbitals but require also a priori the knowledge of the set of virtual orbitals. Such methods have the obvious drawback that the size of the virtual orbital space can be very large.

A solution often used to keep the size of matrices in reasonable limits employs an auxiliary basis set to expand the occupied-virtual products. Such approaches are known in quantum chemistry as resolution of identity 2 or density-fitting 3,4 methods and similar advantages can be found by use of Cholesky decomposition 5 techniques.

Further gain can be achieved by local correlation methods, which take advantage of the locality of purposefully made orbitals to reduce the number of significant terms to calculate in a given method. Hence, both the localization of orbitals [6][7][8][9][10][11][12][13] and the derivation of local correlation methods [14][15][16][17][18][19] have lately been an active area of research.

In the present work it is shown that accurate random phase approximation (RPA) energies can be evaluated with approximations which avoid any explicit reference to virtual orbitals and involve quantities that are computable from occupied orbitals alone, much like what was shown by Surján for MP2 20 . The strategy followed here uses localized virtual orbitals called the projected oscillator orbitals (pOOs) 1 , in the spirit of the projected atomic orbitals techniques [21][22][23] .

In a fairly recent paper 1 , we showed that the pOOs, an original idea by Foster and Boys [24][25][START_REF] Boys | Quantum Theory of Atoms, Molecules, and the Solid State, A Tribute to[END_REF] , could be used in the context of the RPA equations. The core idea behind the pOOs is to construct a set of virtual orbitals directly from a set of occupied localized orbitals (LMOs) by multiplying them with solid spherical harmonics. The orthogonality of these oscillator orbitals with the occupied space is ensured by projection to obtain the projected oscillator orbitals. The pOOs are non-orthogonal among each other which is a source of some complications, although only to the same extent and in the same manner as in the case of the projected atomic orbitals, and the work that has been done over the years in the context of lo-cal correlation methods using projected atomic orbitals 27 can be used. As mentioned in another paper 1 , the pOOs show similarities with other works [START_REF] Kirkwood | [END_REF][29][30][31][32][33][34] .

The main interest of this paper is not to reproduce the full correlation energy with high numerical precision, but only a well-defined part of it, namely the long-range dynamical correlation energy which is usually responsible for the London dispersion forces. It is indeed now well-documented that most of the conventional density functional calculations in the Kohn-Sham framework are unable to grasp the physics of these long-range forces, unless special corrections are added to the total energy [35][36][37][38][39][40][41] . On the other hand it has been demonstrated in earlier works [42][43][44][45] that the essential physical ingredients of London dispersion forces are contained in the rangeseparated hybrid RPA method, where the short-range correlation effects are described within a density functional approximation and the long-range exchange and correlation are handled at the long-range Hartree-Fock and long-range RPA levels. The range-separated hybrid method has the additional advantage of showing a more favorable convergance with basis set size 46 . Note that the use of localized orbitals for dispersion energy calculations has already been proposed since the early works on local correlation methods [47][48][49][50][51] .

In the following, new derivations of the projected oscillator orbitals framework are presented as well as a very brief recall of the local formulation of the random phase approximation equations before showing the performance on full-range and long-range correlation energies. Detailed derivations are provided in the Appendix when needed.

PROJECTED OSCILLATOR ORBITALS

A set of oscillator orbitals (OOs) is produced by multiplying a set of localized molecular orbitals (LMOs) {|i } by solid spherical harmonics centered on the barycenter D i of the LMOs (here for the first order solid spherical harmonic) :

|OO = (r α -D i α )|i (1)
where rα is the α = x, y, z component of the position operator. The set of localized virtual orbitals {|i α } that are the projected oscillator orbitals (pOOs) is constructed by projecting the oscillator orbitals out of the occupied space :

|pOO = P (r α -D i α )|i . = |i α (2)
where P = ( 1 -|k k|) and where the composite index "i α " refers to a pOO generated from the i-th LMO by using the rα harmonic. The pOOs are not orthogonal among each other, and have an overlap (here again between the first order pOOs) :

S iαj β = i|r α P rβ |j = i|r α rβ |j - N LMO k i|r α |k k|r β |j (3)
where N LMO is the number of (localized) occupied orbitals.

In Eq. ( 2), the tail emerging from the projection term |k k|r α |i obviously damages the locality of the pOOs. The Foster-Boys localization criterion for the LMOs 24,[START_REF] Boys | Quantum Theory of Atoms, Molecules, and the Solid State, A Tribute to[END_REF] precisely ensures that the sum of the offdiagonal elements of the rα operators taken between the occupied orbitals is minimized 52 . Hence the pOOs are optimally localized when used in conjunction with LMOs obtained by use of the Foster-Boys criterion.

The virtual pOOs have by design a nodal surface intersecting the occupied orbital centroid and coinciding with the region of the highest electron density of the occupied orbital : this ensures an optimal description of the correlation. Figure 1. Number of virtual orbitals in terms of multiples of the number of occupied orbitals for the atoms and molecules considered in this paper. In filled symbols are the cc-pVXZ basis sets, denoted "vdz", "vtz" and "vqz", and in shallow symbols are the aug-cc-pVXZ, denoted "avdz", "avtz" and "avqz". The sets of pOOs up to a given order each yield a constant number of virtual orbitals in terms of multiples of the number of occupied orbitals, shown by the horizontal lines.

H B r A r C l P S A l C l 2 S i C S 2 M g N a F H C l N e C O S S O 2 H 2 S N O C O 2 N 2 O C O N 2 B C H F S iH 4 H 2 O C 2 H 2 C H 3 C H O C H 3 O H B e C H 3 O C H 3 L i C 2 H 4 C 3 H 6 C H 3 N H 2 H H e N H 3 C 2 H 6 C

Higher order pOOs

Higher order pOOs that have further nodes can be obtained in a similar way, using higher order polynomials. A pOO of general order n hence reads :

|i α...η = P (r α -D i α ) . . . (r η -D i η )|i (4)
It can be shown (see Appendix) that pOOs can be written as a function of pOOs of lower order :

|i α...η = P rα . . . rη |i - n-1 k=1 D i . . . D i n-k |K pOO of order k (5)
In this equation, every term of order k in the sum is actually replicated n k times, each time assigning a different set of k indices to the pOO |K and the remaining indices to the vector component D i . This explains why the terms in the sum are not written with their indices, and might be made clearer with an example, found in the Appendix. As a result of Eq. ( 5), the overlaps involving any order of pOOs are functions of lower order overlaps :

N |M = i| r . . . r n P r . . . r m |j - n k=1 m l=1 k+l =n+m D i . . . D i n-k K|L D j . . . D j m-l (6)
where again the terms of the sum are replicated n k m l times with a different set of k indices assigned to the pOO K|, a different set of l indices assigned to the pOO |L , and the remaining indices to the vectors D i and D j . The pOOs of a given order n and all equations involving those pOOs hence only need ingredients (overlap, etc. . .) already constructed for lower order pOOs plus at most the matrix elements of the multipole operator of order 2n.

Size of the virtual orbital space

The number of first order pOOs generated is 3.N LMO (the x, y and z pOOs constructed for each LMO). In terms of multiples of the number of LMOs, the number of pOOs generated at each order n of spherical harmonics used is given by a figurate number (linear number for first order, triangular number for second order, tetrahedral number for third order, etc. . .), which simply corresponds to the number of unique n-tuplets that can be made with x, y, z

P (n) = 1 2 (n + 1)(n + 2) = n + 2 n (7)
This is summed-up in Table I, which also shows the cummulated number of pOOs generated up to order n for the first 5 orders. As comparison, Figure 1 shows the number of virtual orbital generated by the widely used Dunning basis sets cc-pVDZ, cc-pVTZ, cc-pVQZ, aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ for the atoms and molecules on which the correlation energy are calculated in this paper (those are taken from the set of the small systems compiled by Tkatchenko and Scheffler 39,53 ). This shows that for example the set of pOOs up to third order yields a size of virtual space greater than the one generated by the cc-pVTZ and on par with the one generated by aug-cc-pVTZ. The random phase approximation can be formulated in a number of different ways 44,45,[54][55][56] , but in a local formulation of the "ring coupled-cluster double" (rCCD), the RPA correlation energy can be given by a sum over (localized) pair contributions :

E c = 1 2 Tr K ij T ji (8)
where the amplitudes T are found by solving the Riccati equations R ij = 0, where :

R ij = B ij + A im T mj S + ST im A mj + ST im B mn T nj S (9)
where the matrices are written in the LMOs/pOOs basis :

A ij kαl β = ij|k α l β -k α i|l β j + ij kαl β B ij kαl β = ij|k α l β -ij|l β k α ij kαl β = δ ij f kαl β -S kαl β f ij ( 10
)
where f is the Fock matrix in LMOs/pOOs basis and the two-body integrals between LMOs and pOOs are in physicist's notation. The Riccati equations are solved iteratively in a pseudo-canonical basis ; see our original paper 1 for more detailed derivations and a detailed workflow in the appendix.

CORRELATION ENERGIES

The performance of the pOOs is investigated for the calculation of the RPA correlation energies of a set of Each dot represents a result for one of those molecule, obtained with several different basis. The gray boxes gather the gray dots corresponding to reference results that are obtained with the unmodified cc-pVDZ, aug-cc-pVDZ, cc-pVTZ, aug-cc-pVTZ, cc-pVQZ and aug-cc-pVQZ basis sets (denoted as vdz, avdz, etc. . .). The color dots are the pOO results, with pOOs of first (second and third) order being in blue (red and green). Those pOOs can be obtained in the context of the different reference Dunning basis, hence the nuances of blue, red and green. The black line in both plots roughly connects the centers of the gray boxes.

small molecules 39 already used in our original paper 1 . Both post-HF calculations and range-separated hybrid calculations using the RSHPBE functional with a parameter µ = 0.5 are performed. The occupied orbitals obtained are localized, and the pOOs are generated before calculating either the full-range or the long-range RPA correlation energies.

All calculations were performed with a development version of the quantum chemical program package Molpro 57 , where the construction of the pOOs up to third order has been implemented for the purpose of this paper, together with the corresponding overlaps.

Results of calculations using the cc-pVDZ, cc-pVTZ, cc-pVQZ as well as the aug-cc-pVDZ, aug-cc-pVTZ and aug-cc-pVQZ basis sets are shown. In practice, because of linear dependencies built in the pOOs (see briefly the workflow in the Appendix), the RPA equations are carried with an effective number of independant pseudocanonical basis functions 1 which can be lower than the multiple of the number of occupied orbitals mentioned in Table I and in Figure 1. The pertinent quantities to look at are hence the RPA correlation energies and the number of pseudo-canonical orbitals involved in the calculations. Figure 3 shows this in terms of the fractions of the aug-cc-pVQZ RPA correlation energies recovered as a function of the fractions of the aug-cc-pVQZ number of virtual orbitals involved in the calculations. On the left are shown post-HF results and on the right long-range RPA correlation energy results.

The RPA correlation energies obtained in each basis sets without the use of localized orbitals are shown in gray dots and gathered for clarity within gray boxes. The rest of the data (the RPA correlation energies obtained when using LMOs and pOOs) is shown in colored sym-bols. The results from pOOs obtained within a (aug-)cc-pVXZ basis are shown in diamonds (squares). The pOOs of first order are shown in shades of blue, of second order in shades of red and of third order in shades of green.

A first interesting result is that the performances of a given order of pOOs is not significantly changed by the basis used to construct them : pOOs will perform a certain way independently of the core basis set used (this is even more striking in the case of the long-range results, right of Figure 3). Furthermore, a positive finding is that the calculations involving pOOs recover a good portion of the targeted RPA correlation energies (full-or long-range) for the size of the virtual space involved. As a marker, the black line in the figure is drawn roughly between the center of the gray boxes holding the non-localized RPA energies, and shows the fraction of the aug-cc-pVQZ correlation energy recovered by smaller basis sets. The pOO results lie virtually systematically above that line, meaning the pOOs recover a larger portion of the targeted RPA energies than the basis set having similar numbers of virtual orbitals. Finally, as expected, the fraction of long-range RPA correlation energy recovered is significantly higher than the recovered portion of full-range correlation energy. This is a nice result, as the framework of pOOs was largely thought as a way to build approximations to model the London dispersion forces, as can be seen in our original paper 1 .

CONCLUSION

It was shown that localized virtual orbitals can be constructed easily and systematically by multiplying localized occupied orbitals by harmonic functions of higher than first order. The size of the virtual space obtained is comparable to the size of that generated when using the Dunning basis sets, and the RPA correlation energies calculated with the pOO localized virtual orbitals are close to the energies without local approximations. Because of an interesting iterative feature emerging in the construction of the projected oscillator orbitals, it is easy to generate higher order pOO from generated lower order pOOs. A great number of additional approximation, developed in our original paper 1 in the context of first order projected oscillator orbitals, such that multipolar expansion to approximate the electron repulsion, spherical approximation, etc. could be applied in the future to the higher order pOOs derived here.

ACKNOWLEDGEMENTS

This research was conducted while in the University of Colorado at Boulder and hence was supported through the startup package of Sandeep Sharma. The author would like to thank the late János G. Ángyán for the initial insightful discussions on this subject and on many others. János was a great mentor to me and he is missed very much.

APPENDIX Relations between pOOs

The natural expression for a pOO of order n is given by the binomial theorem :

|i α...η = P (r α -D i α ) . . . (r η -D i η )|i = n k=1 A n k ( 11
)
where A n k is defined for convenience :

A n k = (-1) n-k (D i . . . D i ) n-k P r . . . r k |i (12)
and actually contains n k repetitions with different assignations of the indices α . . . η. For each of those repetitions, a different set of k indices among α . . . η is assigned to the k position operators and the remaining indices are given to the vector components D i . For example, consider a pOO |i αβγ of order 3. The term A 3 1 will be :

A 3 1 = (D i α D i β ) P rγ |i A 3 1 = (D i α D i γ ) P rβ |i A 3 1 = (D i β D i γ ) P rα |i
i.e. it contains 3 1 = 3 repetitions with different assignations of the indices α, β and γ.

We seek to prove Eq. ( 5), which amounts with the current notations to :

|i α...η = A n n - n-1 k=1 D i . . . D i n-k |K pOO of order k (13)
where the same repetitions with different assignations of the indices α . . . η are found now for each terms of the sum. We begin by inserting in Eq. ( 13) the binomial expression of Eq. ( 11) for the pOO |K :

|i α...η = A n n - n-1 k=1 (D i . . . D i ) n-k k m=1 A k m ( 14
)
where we can rewrite the term in the sums using :

(D i . . . D i ) n-k A k m = (-1) k-m (D i . . . D i ) n-k (D i . . . D i ) k-m P r . . . r m |i (15)
= (-1) k+n A n m This yields :

|i α...η = A n n + n-1 k=1 k m=1 (-1) k+(n-1) A n m = A n n + n-1 m=1 n-1 k=m (-1) k+(n-1) A n m ( 16
)
where in the last line the double summation was simply rearranged (think of summations over lines versus over columns of a matrix whose columns are composed of A n m ). In terms of assignations of the indices α . . . η, in Eq. ( 15) there is n k ways to assign indices to the first string of n -k "D i " and k m ways to assign indices to the second string of k -m "D i ". This hence becomes a combinatorial problem, and recovering Eq. ( 11) from Eq. ( 16) boils down to proving that :

n-1 k=m (-1) k+(n-1) n k k m = n m (17)
which should be the number of repeated occurences of the A n m term with different indices in Eq. ( 16). This amounts to prove that :

n-1 k=m (-1) k+(n-1) n k k m n m = 1 (18)
where, manipulating the fraction, the expression to study is :

n-1 k=m (-1) k+(n-1) n -m k -m = n k=m (-1) k+(n-1) n -m k -m -(-1) -1 n -m n -m = n-m k =0 (-1) k n -m k (-1) n-m+1 + 1 = 1 (19)
where the first step consists simply in adding and substracting the n-th element of the sum, the second step is a change of the dummy index of the sum and the final step is a realization that the expression in curvy brackets is the binomial theorem for (1 + (-1)) n-m = 0. This proves Eq. ( 5).

Workflow

Although in our original paper 1 relations to construct all ingredients for the pOOs using solely data from the occupied space are presented (this is the purpose of the projected schemes), for the clarity of the proof-of-concept calculations in this paper, the N pOO × N vir matrix V to transform virtual orbitals to pOOs is introduced :

We assume here that the pOOs can be exanded in the virtual basis, i.e. that the pOO space is a subspace of the virtual space. This should be more rigourously checked in later work. For example, the V matrix corresponding to the first order pOO reads :

V (1) iαa = i|r α |a (21)
and, given Eq. ( 5), repeated in the Appendix as Eq. ( 13), the V matrix corresponding to a N -th order pOO simply reads : ...a matrix for order k pOO (22) The overlap between pOOs is S = VV † and the pseudocanonical basis is obtained by solving the generalized eigenvalue equation

V (N )
fX = SX (23)
Because of linear dependencies built in the construction of the pOOs, the overlap matrix S is in general only positive semidefinite instead of positive definite. Hence the pseudo-canonical basis of effective size N eff ≤ N pOO can alternatively by found with the N pOO × N eff matrix X that similtaneously diagonalizes f and S.

From these relations, one can then construct from the N AO ×N vir tranformation matrix C the following sequential transformation matrices :

C AO→vir → CV † AO→pOO → CV † X AO→pseudo-cano (24)
to easily either try out the pOOs, draw them, or calculate energies.

Figure 2 .

 2 Figure 2. Examples of pOOs generated from the oxygen lone pair of the formaldehyde (left panels) by applying first-, secondand third order harmonics (second, third and last panels). The first row shows contour taken on the blue plane in the sketch molecule on the left, the second row shows contour taken on the red plane.

Figure 3 .

 3 Figure 3. Results of calculations of the RPA correlation energies (left : full-range, right : long-range) of a small set of molecules.Each dot represents a result for one of those molecule, obtained with several different basis. The gray boxes gather the gray dots corresponding to reference results that are obtained with the unmodified cc-pVDZ, aug-cc-pVDZ, cc-pVTZ, aug-cc-pVTZ, cc-pVQZ and aug-cc-pVQZ basis sets (denoted as vdz, avdz, etc. . .). The color dots are the pOO results, with pOOs of first (second and third) order being in blue (red and green). Those pOOs can be obtained in the context of the different reference Dunning basis, hence the nuances of blue, red and green. The black line in both plots roughly connects the centers of the gray boxes.

  |i α...η = Nvir a V iα...ηa |a

  iα...ηa = i|r α . . . rη |a -

Table I .

 I Number of virtual orbitals generated at and up to a given order of harmonics used.

				LOCAL RPA EQUATIONS
	order	figurate	number of	cumulant
	of the	number	pOOs	number of
	pOOs		generated	pOOs
	1	3 1	3 × NLMO	3 × NLMO
	2	4 2	6 × NLMO	9 × NLMO
	3	5 3	10 × NLMO	19 × NLMO
	4	6 4	15 × NLMO	34 × NLMO
	5	7 5	21 × NLMO	51 × NLMO


Davide Giraudo
email: davide.giraudo1@univ-rouen.fr

HÖLDERIAN WEAK INVARIANCE PRINCIPLE UNDER THE MAXWELL AND WOODROOFE CONDITION

Keywords: Mathematics Subject Classification. 60F05; 60F17 Invariance principle, martingales, Hölder spaces, strictly stationary process

come

Hölderian weak invariance principle under Maxwell and

Woodroofe condition

Introduction and main results

Let (Ω, F , µ) be a probability space and let T : Ω → Ω be a measure-preserving bijective and bi-measurable map. Let M be a sub-σ-algebra of F such that T M ⊂ M. If f : Ω → R a measurable function, we denote S n (T, f ) := n-1 j=0 f • T j and W (n, f, T, t) := S [nt] (T, f ) + (nt -[nt])f • T [nt] .

(1.1)

We shall write S n (f ) and W (n, f, t) for simplicity, except when T is replaced by T 2 . An important problem in probability theory is the understanding of the asymptotic behavior of the process (n -1/2 W (n, f, t), t ∈ [0, 1]) n 1 . Conditions on the quantities E[S n (f ) | T M] and S n (f )-E[S n (f ) | T -n M] have been investigated. The first result in this direction was obtained by Maxwell and Woodroofe [MW00]: if f is M-measurable and

+∞ n=1 E [S n (f ) | M] 2 n 3/2 < +∞, (1.2)
then (n -1/2 S n (f )) n 1 converges in distribution to η 2 N , where N is normally distributed and independent of η. Then Volný [Vol06] proposed a method to treat the nonadapted case. Peligrad and Utev [PU05] proved the weak invariance principle under condition (1.2). The nonadapted case was addressed in [Vol07]. Peligrad and Utev also showed that condition (1.2) is optimal among conditions on the growth of the sequence (

E [S n (f ) | M] 2 ) n 1 : if +∞ n=1 a n E [S n (f ) | M] 2 n 3/2 < ∞ (1.3)
for some sequence (a n ) n 1 converging to 0, the sequence (n -1/2 S n (f )) n 1 is not necessarily stochastically bounded (Theorem 1.2. of [PU05]). Volný constructed [START_REF]Martingale approximation and optimality of some conditions for the central limit theorem[END_REF] an example satisfying (1.3) and such that the sequence S n (f ) -1 2 S n (f ) n 1 admits two subsequences which converge weakly to two different distributions.

Let us denote by H α the space of Hölder continuous functions, that is, the functions x : [0, 1] → R such that x Hα := sup 0 s<t 1 |x(t) -x(s)| /(t -s) α + |x(0)| is finite. Since the paths of Brownian motion belong almost surely to H α for each α ∈ (0, 1/2) as well as W (n, f, •), we can investigate the weak convergence of the sequence (n -1/2 W (n, f, •)) n 1 in the the space H α , for 0 < α < 1/2. The case of i.i.d. sequences and stationary martingale difference sequences have been addressed respectively by Račkauskas and Suquet (Theorem 1 of [RS03]) and Giraudo (Theorem 2.2 of [START_REF] Giraudo | Holderian weak invariance principle under a Hannan type condition[END_REF]). In this note, we focus on conditions on the sequences

(E[S n (f ) | M]) n 1 and (S n (f ) -E[S n (f ) | T -n M]) n 1 . Theorem 1.1. Let p > 2 and f ∈ L p . If +∞ k=1 E[S k (f ) | M p k 3/2 < +∞, +∞ k=1 S k (f ) -E[S k (f ) | T -k M p k 3/2 < +∞, (1.4)
then the sequence n -1/2 W (n, f ) n 1 converges weakly to the process √ ηW in H 1/2-1/p , where W is the Brownian motion and the random variable η is independent of W and is given by η = lim n→+∞ E S n (f ) 2 | I /n (where I is the σ-algebra of invariant sets and the limit is in the L 1 sense).

Of course, if f is M-measurable, all the terms of the second series vanish and we only have to check the convergence of the first series.

Remark 1.2. If the sequence (f • T j ) j 0 is a martingale difference sequence with respect to the filtration (T -i M), then condition (1.4) is satisfied if and only if the function f belongs to L p , hence we recover the result of [START_REF] Giraudo | Holderian weak invariance principle under a Hannan type condition[END_REF]. However, if the sequence (f • T j ) j 0 is independent, (1.4) is stronger than the sufficient condition t p µ {|f | > t} → 0. This can be explained by the fact that the key maximal inequality (2.9) does not include the quadratic variance term which appears in the martingale inequality. In Remark 1 (after the proof of Theorem 1) in [START_REF] Peligrad | A maximal Lp-inequality for stationary sequences and its applications[END_REF], a version of (2.9) with this term is obtained. In our context it seems that it does not follow from an adaptation of the proof.

Remark 1.3. In [START_REF] Giraudo | Holderian weak invariance principle under a Hannan type condition[END_REF], the conclusion of Theorem 1.1 was obtained for an M-measurable f under the condition

∞ i=1 E f | T i M -E f | T i+1 M p < ∞, (1.5)
which holds as soon as

+∞ k=1 E f • T k | M p k 1/p < +∞, (1.6)
while (1.4) holds as soon as

+∞ k=1 E f • T k | M p √ k < +∞. (1.7)
Therefore, (1.7) gives a better sufficient condition than (1.6) if we seek for conditions relying only on

E f • T k | M p k 1 .
However, (1.5) gives the existence of a martingale approximation in the following sense: there exists a martingale difference m ∈ L p (M) such that

W (n, f ) -W (n, m) H 1/2-1/p p,∞ = o( √ n).
(1.8) Indeed, define for an integrable function h and a non-negative integer i, P i (h

) := E h | T i M - E h | T i+1 M . If f satisfies (1.5), then we set m := i 0 P 0 U i f . Then for any K 1, the equality f -m = K i=0 P i (f ) -P 0 U i f + +∞ i=K+1 P i (f ) -P 0 U i f holds. Since K i=0 P i (f ) -P 0 U i f
may be written as (I -U )g K , where g K is such that t p µ {|g K | > t} → 0 as t goes to infinity, we get, by inequalities (2.4) and (2.5) of [START_REF] Giraudo | Holderian weak invariance principle under a Hannan type condition[END_REF] that

lim sup n→+∞ 1 √ n W (n, f ) -W (n, m) H 1/2-1/p p,∞ i K+1 lim sup n→+∞ 1 √ n W (n, P i (f ))) H 1/2-1/p p,∞ + W n, P 0 U i (f ) H 1/2-1/p p,∞ .
We conclude by Proposition 2.3 of [START_REF] Giraudo | Holderian weak invariance principle under a Hannan type condition[END_REF].

The following condition (in the spirit of Maxwell and Woodroofe's one) is sufficient for a martingale approximation in the sense of (1.8):

+∞ k=1 E[S k (f ) | M p k 1+1/p < +∞.
(1.9) Indeed, Theorem 2.3 of [CM14] gives a martingale differences sequence m • T i i 0 such that lim n→+∞ n -1/p S n (f -m) p = 0. Using Serfling arguments (see [Ser70]), we get that (1.9) implies

lim n→+∞ n -1/p max 1 i n |S i (f -m)| p = 0.
(1.10)

Note that for a function h, by Lemma A.2 of [START_REF] Markevičiūtė | Functional central limit theorems for sums of nearly nonstationary processes[END_REF],

n -1/2 W (n, h) H 1/2-1/p p,∞
2n -1/p max 1 j n |S j (f )| p,∞ , hence by (1.10), the martingale approximation (1.8) holds. Furthermore, using the construction given in [DV08,Dur09], in any ergodic dynamical system of positive entropy one can construct a function satisfying condition (1.4) but not (1.5) and vice versa.

Remark 1.4. For the ρ-mixing coefficient defined by

ρ(n) = sup Cov(X, Y )/( X 2 Y 2 ), X ∈ L 2 (σ(f • T i , i 0), Y ∈ L 2 (σ(f • T i , i n)) ,
Lemma 1 of [START_REF] Peligrad | A maximal Lp-inequality for stationary sequences and its applications[END_REF] shows that for an adapted process, condition (1.4) is satisfied if the series

∞ n=1 ρ 2/p (2 n ) converges. However, the conclusion of Theorem 1.1 holds if t p µ {|f | > t} → 0 and ∞ n=1 ρ(2 n ) converges (see Theorem 2.3, [Gir16a]), which is less restrictive.
It turns out that even in the adapted case, condition (1.4) is sharp among conditions on E[S k (f ) | M p in the following sense.

Theorem 1.5. For each sequence (a n ) n 1 converging to 0 and each real number p > 2, there exists a strictly stationary sequence (f

• T j ) j 0 and a sub-σ-algebra M such that T M ⊂ M, ∞ n=1 a n n 3/2 E[S n (f ) | M] p < ∞, (1.11) but the sequence n -1/2 W (n, f, t) n 1 is not tight in H 1/2-1/p .
Remark 1.6. Using the inequalities in [START_REF] Peligrad | A maximal Lp-inequality for stationary sequences and its applications[END_REF] in order to bound E [S n (f ) | T M] 2 , we can see that the constructed f in the proof of Theorem 1.5 satisfies the classical Maxwell and Woodroofe condition (1.2) (the fact that p is strictly greater than 2 is crucial), hence the weak invariance principle in the space of continuous functions takes place. However, it remains an open question whether condition (1.11) implies the central limit theorem or the weak invariance principle (in the space of continuous functions).

Proofs

We may observe that condition (1.4) implies by Theorem 1 of [START_REF] Peligrad | A maximal Lp-inequality for stationary sequences and its applications[END_REF] that the sequence (S n (f )/ √ n) n 1 is bounded in L p ; nevertheless the counter-example given in Theorem 2.6 of [Gir16a] shows that we cannot deduce the weak invariance principle from this.

We shall rather work with a tighness criterion. The analogue of the continuity modulus in

C[0, 1] is ω α , defined by ω α (x, δ) = sup 0<|t-s|<δ |x(t) -x(s)| |t -s| α , x : [0, 1] → R, δ ∈ (0, ]. Define H o α [0, 1] := {x ∈ H α [0, 1], lim δ→0 ω α (x, δ) = 0}.
We shall essentially work with the space H o α [0, 1] which, endowed with

• α : x → ω α (x, 1) + |x(0)|, is a separable Banach space (while H α [0, 1] is not). Since the canonical embedding ι : H o α [0, 1] → H α [0, 1] is continuous, each convergence in distribution in H o α [0, 1] also takes place in H α [0, 1].
Let us state the tighness criterion we shall use (Theorem 13 of [Suq99]).

Proposition 2.1. Let α ∈ (0, 1). A sequence of processes (ξ n ) n 1 with paths in H o α [0, 1] and such that ξ n (0) = 0 for each n is tight in H o α [0, 1] if and only if ∀ε > 0, lim δ→0 sup n→+∞ µ {ω α (ξ n , δ) > ε} = 0. (2.1)
In order to prove the weak convergence in H o α [0, 1], it suffices to prove the convergence of the finite dimensional distributions and establish tighness in this space.

A maximal inequality. For p > 2, we define

h p,∞ := sup A∈F µ(A)>0 1 µ(A) 1-1/p E[|h| 1 A ]. (2.2)
This norm is linked to the tail function of h by the following inequalities (see Exercice 1.1.12 p. 13 in [Gra14]):

sup t>0 t p µ {|h| > t} 1/p h p,∞ p p -1 sup t>0 t p µ {|h| > t} 1/p . (2.3)
As a consequence, if N is an integer and h 1 , . . . , h n are functions, then max

1 j N |h j | p,∞ p p -1 N 1/p max 1 j N |h j | p,∞ .
(2.4)

For a positive n 1, a function f : Ω → R and a measure-preserving map T , we define

M (n, f, T ) := max 0 i<j n |S j (T, f ) -S i (T, f )| (j -i) 1/2-1/p . (2.5) By Lemma A.2 of [MSR12]
, the Hölderian norm of a polygonal line is reached at two vertices, hence

M (n, f, T ) = n 1/2-1/p W (n, f, T, •) H 1/2-1/p (2.6) Applying Proposition 2.3 of [Gir16b],
we can find for each p > 2 a constant C p depending only on p such that if (m • T i ) i 1 is a martingale difference sequence, then for each n,

1 √ n W (n, m, T, •) H 1/2-1/p p,∞ C p m p .
(2.7)

In the sequel, fix such a constant C p that we shall choose greater than 6 • 2 1/p p/(p -1). We denote by U the Koopman operator associated with T , that is, for each f : Ω → R and each

ω ∈ Ω, (U f )(ω) = f (T ω).
Definition 2.2. Let H be a closed subspace of L p . Let P be a linear operator from H to itself.

We say that (H, P ) satisfies condition (C) if (1) the inclusion U -1 H ⊂ H holds (respectively the inclusion U H ⊂ H holds);

(2) P is power bounded on H, that is, for each h ∈ H,

K(P ) := sup n 1 sup h∈H\{0} P n h p h p < +∞ ;
(2.8)

(3) if h ∈ H is such that P h = 0, then the sequence (h • T i ) i 0 is a martingale difference sequence with respect to the filtration T -i M i 0 (respectively

T -i-1 M i 0 ); (4) P U -1 f = f for each f ∈ H (respectively P U f = f for each f ∈ H).
Let us give two examples of subspace H and operator P satisfying condition (C).

(1) Let H be the subspace of L p which consists of M-measurable functions and

P h := E [U h | M].
Then (H, P ) satisfies condition (C).

(2) Let H be the subspace of L p which consists of functions h such that E [h | M] = 0 and

P h := U -1 h -E U -1 h | M . Then (H, P ) satisfies condition (C).
The goal of this subsection is to establish the following maximal inequality.

Proposition 2.3. Let T : Ω → Ω be a bijective and bi-measurable measure-preserving map.

Let H be a closed subspace of L p . Let r be a positive integer. For each , operator P from H to itself such that (H, P ) satisfies condition (C), each f ∈ H and each integer n satisfying

2 r-1 n < 2 r , M (n, f, T ) p,∞ C p n 1/p   (1 + K (P )) f p + K p r-1 j=0 2 -j/2 2 j -1 i=0 P i f p   , (2.9)
where

K p = 2 1/p-1/2 + 2 1/2 (1 + K(P )).
If H is a closed subspace of L p and P : H → H an operator such that (H, P ) satisfies condition (C), we define for f ∈ H the quantity (2.11)

f MW(p,P ) := +∞ j=0 2 -j/2 2 j -1 i=0 P i f p (2.
Note that MW(p, P ) endowed with • MW(p,P ) is a Banach space. Combining Proposition 2.3 and (2.6), we derive the following bound for the Hölderian norm of the partial sum process.

Corollary 2.4. Let H be a closed subspace of L p and let P be an operator from H to itself such that (H, P ) satisfies the condition (C). Then there exists a constant C = C(p, P ) such that for each n, and each h ∈ H,

1 √ n W (n, h) H 1/2-1/p p,∞ C h MW(p,P ) (2.12)
The proof of Proposition 2.3 is in the same spirit as the proof of Theorem 1 of [START_REF] Peligrad | A maximal Lp-inequality for stationary sequences and its applications[END_REF], which is done by dyadic induction. To do so, we start from the following lemma: Lemma 2.5. For each positive integer n, each function h : Ω → R and each measure-preserving map T : Ω → Ω, the following inequality holds:

M (n, h, T ) 6 max 0 k n h • T k + 1 2 1/2-1/p M n 2 , h + h • T, T 2 .
(2.13)

Proof. First, notice that if 1 j n, then j = 2 j 2 or j = 2 j 2 + 1, hence

S j (h) -S 2[ j 2 ] (h) max 0 k n h • T k . (2.14) Similarly, we have S i (h) -S 2[ i+2 2 ] (h) 2 max 0 k n h • T k . (2.15) It thus follows that M (n, h, T ) 4 max 0 k n h • T k + max 0 i<j n S 2[ j 2 ] (h) -S 2[ i+2 2 ] (h) (j -i) 1/2-1/p . (2.16)
Notice that if j i + 4, then

1 j 2 - i + 2 2 j -i 2 ,
(2.17) and we derive the bound max

0 i<j n S 2[ j 2 ] (h) -S 2[ i+2 2 ] (h) (j -i) 1/2-1/p 1 2 1/2-1/p max 0 u<v [ n 2 ] S v (T 2 , h + h • T ) -S u (T 2 , h + h • T ) (v -u) 1/2-1/p + + max 0 i<j n j i+4 S 2[ j 2 ] (h) -S 2[ i+2 2 ] (h) .
Since for j i + 4, the number of terms of the form h

• T q involved in S 2[ j 2 ] (h) -S 2[ i+2 2 ] (h) is at most 2, we conclude that max 0 i<j n S 2[ j 2 ] (h) -S 2[ i+2 2 ] (h) (j -i) 1/2-1/p 1 2 1/2-1/p M n 2 , h + h • T, T 2 + + 2 max 0 k n h • T k .
Combining this inequality with (2.16), we obtain (2.13), which concludes the proof of Lemma 2.5. Now, we establish inequality (2.9) by induction on r.

Proof of Proposition 2.3. We first assume that P U -1 = Id and U -1 H ⊂ H. We check the case r = 1. Then necessarily n = 1 and the expression M (n, f, t) reduces to f . Since C p and K p are greater than 1, the result is a simple consequence of the triangle inequality applied to f -U -1 P f and U -1 P f . Now, assume that Proposition 2.3 holds for some r and let us show that it takes place for r + 1. We thus consider an integer n such that 2 r n < 2 r+1 , a function f ∈ H, a measurepreserving map T : Ω → Ω bijective and bi-measurable, and a sub-σ-algebra M satisfying T M ⊂ M, a closed subspace H of L 2 such that U -1 H ⊂ H and an operator P : H → H such that (H, P ) satisfies condition (C) with P U -1 = Id and we have to show that (2.9) holds with r + 1 instead of r. First, using inequality M (n, f ) M (n, f -U -1 P f ) + M (n, U -1 P f ) and Lemma 2.5 with h := U -1 P f , we derive

M (n, f, T ) M n, f -U -1 P f, T + 6 max 0 k n U -1 P f • T k + + 1 2 1/2-1/p M n 2 , (I + U )U -1 P f, T 2 , (2.18)
hence taking the norm • p,∞ , we obtain by (2.4) that

M (n, f, T ) p,∞ M (n, f -U -1 P f, T ) p,∞ + 6(n + 1) 1/p p p -1 U -1 P f p + + 1 2 1/2-1/p M n 2 , (I + U )U -1 P f, T 2 p,∞ . (2.19)
By inequality (2.7) and accounting the fact that 6 • (n + 1) 1/p p/(p -1) C p n 1/p , we obtain

M (n, f, T ) p,∞ C p n 1/p f -U -1 P f p + C p n 1/p P f p + + 1 2 1/2-1/p M n 2 , (I + U )U -1 P f, T 2 p,∞ . (2.20) Since 2 r-1
[n/2] < 2 r , we may apply the induction hypothesis to the integer [n/2], the function h := (I + U -1 )P f , T 2 instead of T and P 2 instead of P . This gives

n 2 -1/p M n 2 , h, T 2 p,∞ C p 1 + K P 2 h p + + C p K p r-1 j=0 2 -j/2 2 j -1 i=0 P 2i I + U -1 P f p , (2.21)
where K p = 2 1/p-1/2 + 2 1/2 1 + K(P 2 ) . Notice that h p 2 P f p , and by item 4 of Definition 2.2, it follows that

2 j -1 i=0 P 2i I + U -1 P f = 2 j -1 i=0 P 2i+1 f + P 2i f = 2 j+1 -1 i=0 P i f.
(2.22)

Accounting the inequality K P 2 K (P ) and K p K p , we have

n 2 -1/p M n 2 , h, T 2 p,∞ 2 (1 + K (P )) C p P f p + C p K p r-1 j=0 2 -j/2 2 j+1 -1 i=0 P i f p = 2 (1 + K (P )) C p P f p + 2 1/2 C p K p r j=1 2 -j/2 2 j -1 i=0 P i f p and we infer M n 2 , h, T 2 p,∞ n 2 1/p 2 (1 + K (P )) -K p √ 2) C p P f p + n 1/p 2 1/2-1/p C p K p r j=0 2 -j/2 2 j -1 i=0 P i f p . (2.23)
Pluggling this into (2.20), we derive

M (n, f, T ) p,∞ C p n 1/p (1 + K (P )) f p + n 1/p C p K p r j=0 2 -j/2 2 j -1 i=0 P i f p + + C p n 1/p 1 + 2 1-1/p (1 + K(P )) -2 1/2-1/p K p P f p . (2.24)
The definition of K p implies that 2 1/p-1/2 -√ 2(1 + K(P )) -K p = 0, hence (2.9) is established. This concludes the proof of Proposition 2.3 in the case P U -1 = Id.

When P U = Id and U H ⊂ H we do the same proof, but replacing each occurrence of U -1 by U . This ends the proof of Proposition 2.3. 2.2. Proof of Theorem 1.1. Since the convergence of the finite dimensional distributions is contained in the main result of [Vol07], the only difficulty in proving Theorem 1.1 is to establish tightness. To this aim, we shall proceed as in the proof of Theorem 5.3 in [START_REF] Cuny | Limit theorems under the Maxwell-Woodroofe condition in Banach spaces[END_REF].

Proposition 2.6. Let T be a measure preserving map, H a closed subspace of L p (p > 2) and let P be an operator from H to itself such that (H, P ) satisfies condition (C). Assume that h is an element of H such that h MW(p,P ) < +∞ Then the sequence (n

-1/2 W (n, h)) n 1 is tight in H 1/2-1/p . Proof. Let us define V n := n-1 i=0 P i . Using V n V k p K(P ) min k V n p , n V k p , we
derive that for each f ∈ MW(p, P ),

V 2 n f MW(p,P ) 2 n K(P )    V 2 n f p 2 n/2 + k n+1 V 2 k f p 2 k/2    (2.25)
which goes to 0 as n goes to infinity. If m 1 is an integer and if n is such that

2 n m < 2 n+1 , then V m f MW(p,P ) m K(P ) m n k=0 V 2 k f MW(p,P ) K(P ) m n k=0 2 k ε k , (2.26)
where (ε k ) k 1 is a sequence converging to 0. This entails that the operator P is mean-ergodic on MW(p, P ). Furthermore, since P has no non trivial fixed points on the Banach space MW(p, P ), • MW(p,P ) , we derive by Theorem 1.3 p.73 of [Kre85] that the subspace (I -P )MW(p, P ) is dense in MW(p, P ) for the topology induced by the norm • MW(p,P ) .

Let h ∈ H be such that h MW(p,P ) < +∞ and x > 0. We can find f ∈ (I -P )MW(p, P ) such that h -f MW(p,P ) < x. Consequently, using Corollary 2.4, we derive that for each positive ε and δ,

µ ω 1/2-1/p 1 √ n W (n, h), δ > 2ε ε -p x + µ ω 1/2-1/p 1 √ n W (n, f ), δ > ε . (2.27)
Now, since the function f belongs to (I -P )MW(p, P ), we can find

f ′ ∈ MW(p, P ) such that f = f ′ -P f ′ . If P U -1 = Id, then we write f = f ′ -U -1 P f ′ + (U -1 -I)f ′ and if P U = Id, then f = f ′ -U P f ′ + (U -I)f ′ .
In other words, f admits a martingale-coboundary decomposition in L p (since f ′ belongs to L p ). Consequently, by Corollary 2.5 of [START_REF] Giraudo | Holderian weak invariance principle under a Hannan type condition[END_REF], the sequence (n -1/2 W (n, f )) n 1 is tight in H 1/2-1/p . By Proposition 2.1 and (2.27), we derive that for each positive ε and x,

lim δ→0 lim sup n→+∞ µ ω 1/2-1/p 1 √ n W (n, h), δ > 2ε ε -p x.
(2.28)

Since x is arbitrary we conclude the proof of (2.6) by using again Proposition 2.1.

Proof of Theorem 1.1.

Writing f = E [f | M] + f -E [f | M],
the proof reduces (as mentioned in the begining of the section) to establish tightness in

H o 1/2-1/p [0, 1] of the sequences (W n ) n 1 := n -1/2 W (n, E [f | M]) n 1 and (W ′ n ) n 1 := n -1/2 W (n, f -E [f | M]) n 1 . • Tightness of (W n ) n 1 . We define P (f ) := E [U f | M] and H := {f ∈ L p , f is M-measurable} .
(2.29)

Then (H, P ) satisfies condition (C). Since

n-1 i=0 P i (E [f | M]) = E [S n (f ) | M] ,
(2.30) the convergence of the first series in (1.4) is equivalent to f ∈ MW(p, P ) (by Lemma 2.7 of [PU05]). By Proposition 2.6, we derive that the sequence (

W n ) n 1 is tight in H o 1/2-1/p [0, 1]. • Tightness of (W ′ n ) n 1 . We define P (f ) := U -1 f -E U -1 f | M and H := {f ∈ L p , E [f | M] = 0} .
(2.31)

Since for each f ∈ H and each k 1, P k f p 2 f p , (H, P ) satisfies condition (C) (see the proof of Proposition 2 in [Vol07] for the other conditions). Since

P (E [f | M]) = 0, we have n i=1 P i (f -E [f | M]) = n i=1 P i f = U -n S n (f ) -E S n (f ) | T -n M , (2.32)
hence the convergence of the second series in (1.4) implies that f belongs to MW(p, P ) (by Lemma 37 of [START_REF] Merlevède | Rosenthal-type inequalities for the maximum of partial sums of stationary processes and examples[END_REF]). By Proposition 2.6, we derive that the sequence (W

′ n ) n 1 is tight in H o 1/2-1/p [0, 1]
. This ends the proof of Theorem 1.1. 2.3. Proof of Theorem 1.5. We take a similar construction as in the proof of Proposition 1 of [START_REF] Peligrad | A maximal Lp-inequality for stationary sequences and its applications[END_REF]. We consider a non-negative sequence (a n ) n 1 , and a sequence (u k ) k 1 of real numbers such that

u 1 = 1, u 2 = 2, u p/2+1 k + 1 < u k+1 for k 3 and a t k -2 for t u k .
(2.33)

Notice that since p > 2, the conditions (2.33) are more restrictive than that of the proof of Proposition 1 of [START_REF] Peligrad | A maximal Lp-inequality for stationary sequences and its applications[END_REF]. If i = u j for some j 1, then we define p i := cj/u 1+p/2 j and p i = 0 otherwise. Let (Y k ) k 0 be a discrete time Markov chain with the state space Z + and transition matrix given by p k,k-1 = 1 for k 1 and p 0,j-1 := p j , j 1. We shall also consider a random variable τ which takes its values among non-negative integers, and whose distribution is given by µ(τ = j) = p j . Then the stationary distribution exists and is given by

π j = π 0 ∞ i=j+1 p i , j 1, where π 0 = 1/E[τ ].
(2.34)

We start from the stationary distribution (π j ) j 0 and we take g(x) := 1 x=0 -π 0 , where π 0 = µ {Y 0 = 0}. We then define f • T j = X j := g(Y j ).

It is already checked in [PUW07] that the sequence (X j ) j 0 satisfies (1.11), where M = σ(X k , k j) and S n = n j=1 X j . To conclude the proof, it remains to check that the sequence n -1/2 W (n, f, T ) n 1 is not tight in H o 1/2-1/p , which will be done by disproving (2.1) for a particular choice of ε. To this aim, we define

T 0 = 0, T k = min {t > T k-1 | Y t = 0} , τ k = T k -T k-1 , k 1.
(2.35)

Then (τ k ) k 1 is an independent sequence and each τ k is distributed as τ and

S T k = k j=1 (1 -π 0 τ j ) = k -π 0 T k . (2.36)
Let us fix some integer K greater than E[τ ]. Let δ ∈ (0, 1) be fixed and n an integer such that 1/n < δ. Then the inequality

1 (nK) 1/p max 0 i<j nK j-i nδ |S j -S i | (j -i) 1/2-1/p 1 (nK) 1/p 1 {T n Kn} × × max 1 k n S T k -S T k-1 (T k -T k-1 ) 1/2-1/p 1 {|T k -T k-1 | nδ} (2.37)
takes place. By (2.35) and (2.36), this can be rewritten as Since τ 1 almost surely, the following inclusions take place for n > (2/π 0 ) p :

1 (nK) 1/p max 0 i<j nK j-i nδ |S j -S i | (j -i) 1/2-1/p 1 (nK) 1/p 1 {T n Kn} × × max 1 k n |1 -π 0 τ k | τ 1/2-1/p k 1 {τ k nδ} .
A n (π 0 /(2K 1/p )) ⊃ π 0 τ 1/2+1/p -τ -1/2+1/p π 0 /(2K 1/p )(Kn) 1/p ∩ {τ nδ} ⊃ τ 1/2+1/p 1 + π 0 n 1/p /2 π 0 ∩ {τ nδ} ⊃ τ 1/2+1/p n 1/p ∩ {τ nδ} = n 2/(p+2) τ nδ .

Consequently, for j large enough, we deduce that the last term of (2.44) is equal to 1. Since 1 (nK) 1/p max 0 i<j nK j-i nδ

µ
|S j -S i | (j -i) 1/2-1/p ω 1/2-1/p 1 √ nK W (nK, f ), δ , (2.46)
we derive that (2.1) does not hold with ε = π 0 /(2K 1/p ). This finishes the proof of Theorem 1.5.

  10) and the vector space MW(p, P ) := f ∈ H | f MW(p,P ) < +∞ .

1 -

 1 (2.38) Defining for a fixed C the eventA n (C) := |1 -π 0 τ | τ 1/2-1/p C(Kn) 1/p ∩ {τ nδ} , (2.39)we obtain by independence of (τ k ) 1/p max0 i<j nK j-i nδ |S j -S i | (j -i) 1/2-1/p C    1 -(1 -µ(A n (C))) n -µ {T n > Kn} . (2.40)By the law of large numbers, we obtain, accountingK > E[τ ](1 -µ(A n (C))) n . (2.41)We choose C := π 0 /(2K 1/p ). Considering the integers n of the form u 1/p max 0 i<j nK j-i nδ |S j -S i | (j -i) 1/2-1/p

  |S j -S i | (j -i) 1/2-1/p 1 -µ {τ = u j })

  A u (p+2)/2Since τ takes only integer values among u l 's and u

								(p+2)/2 j	δ < u j+1 (by (2.33) and the fact
	that δ < 1), we obtain in view of (2.42), that
	lim sup n→∞	µ	  	1 (nK) 1/p max 0 i<j nK j-i nδ
				j	π 0 2K 1/p	µ u	(p+2)/2 j	2/(p+2)	τ	u	(p+2)/2 j	δ .	(2.43)

Acknowledgements

The author would like to thank an anonymous referee for many valuable comments which improved the presentation of the paper and led to a shorter proof of Theorem 1.1.