Yohann Uguen
email: yohann.uguen@insa-lyon.fr

Florent De

Steven Derrien
email: steven.derrien@univ-rennes1.fr

A high-level synthesis approach optimizing accumulations in floating-point programs using custom formats and operators

High-level synthesis (HLS) is a big step forward in terms of design productivity. However, it restricts data-types and operators to those available in the C language supported by the compiler. The present work lifts this restriction: it is a case study of enhancing an HLS design flow with non-standard operators, which can then be automatically optimized for their application context. The focus here is on widely used summation-reduction patterns. A source-to-source compiler rewrites, inside critical loop nests of the input C code, selected floating-point additions into sequences of simpler operator using non-standard arithmetic formats. This enables hoisting floating-point management out the loop. What remains inside the loop is a sequence of fixedpoint additions whose size is computed to enforce a user-specified, application-specific accuracy constraint on the result. Evaluation of this method demonstrates significant improvements in the speed/resource usage/accuracy trade-off.

I. INTRODUCTION

Many case studies have demonstrated the potential of Field-Programmable Gate Arrays (FPGAs) as accelerators for a wide range of applications, from scientific or financial computing to signal processing and cryptography. FPGAs offer massive parallelism and programmability at the bit level. This enables programmers to exploit a range of techniques that avoid many bottlenecks of classical von Neumann computing: dataflow operation without the need of instruction decoding; massive register and memory bandwidth, without contention on a register file and single memory bus; operators and storage elements tailored to the application in nature, number and size.

However, to unleash this potential, development costs for FPGAs are orders of magnitude higher than classical programming. High performance and high design costs are the two faces of the same coin.

Hardware design flow and High-level synthesis: To address this, languages such as C or Java are increasingly being considered as hardware description languages. This has many advantages. The language itself is more widely known than any HDL. The sequential execution model makes designing and debugging much easier. One can even use software execution on a processor for simulation. All this drastically reduces development time.

The process of compiling a software program into hardware is called High-Level Synthesis (HLS), with tools such as Vi-vadoHLS [START_REF]Vivado Design Suite User Guide: High-Level Synthesis[END_REF] or Catapult C1 among others [START_REF] Nane | A survey and evaluation of fpga high-level synthesis tools[END_REF]. These tools are in charge of turning a C description into a circuit. This task requires to extract parallelism from sequential programs constructs (e.g. loops) and expose this parallelism in the target design. Today's HLS tools are reasonably efficient at this task, and can automatically synthesize highly efficient pipelined dataflow architectures.

They however miss one important feature: they are not able to tailor operators to the application in size, and even less in nature. This comes from the C language itself: its high-level datatypes and operators are limited to a small number (more or less matching the hardware operators present in mainstream processors). Any high-level C description must therefore live with this constraint. The broader objective of this work is to address this limitation.

Arithmetic in HLS: To better exploit the freedom offered by hardware and FPGAs, HLS vendors have enriched the C language with integer and fixed-point types of arbitrary size 2 . However the operations on these types remain limited to the basic arithmetic and logic ones. Exotic or complex operators (for instance for finite-field or floating-point arithmetic) may be encapsulated in a C function that is called to instantiate the operator. This is actually what happens when ones processes floating-point code through HLS.

Such low-level descriptions are currently taken from a parameterized library or operator generators [START_REF] De | High-Performance Computing Using FPGAs[END_REF]. The latter enables many opportunities of compiler optimizations to the operators (not all of which are currently exploited):

• Constant propagation is well understood in compilers, and allows to use hardware multipliers and dividers specialized for a constant argument. • Other algebraic transformations may be detected using classical instruction selection techniques, and exploited. For instance, in hardware, a squarer or a cuber may cost much less than a multiplier [START_REF] De | High-Performance Computing Using FPGAs[END_REF]. • Operators may be shared and fused [START_REF] Langhammer | Floating point datapath synthesis for FPGAs[END_REF]. In these examples, the compiler directs the construction of the arithmetic, but the opposite can also be true. For instance, the pipeline levels necessary in large floating-point operators can be exploited as tiling memory [START_REF] Alias | FPGA-specific synthesis of loop-nests with pipelined computational cores[END_REF]. In this case, it is the arithmetic datapath inside the loop that defines the parameters of a loop transformation.

The case study in this work is a program transformation that applies to floating-point additions on a loop's critical path. It decomposes them into elementary steps, resizes the corresponding subcomponents to guarantee some userspecified accuracy, and merges and reorders these components to improve performance. The result of this complex sequence of optimizations could not be obtained from an operator generator, since it involves global loop information.

For this purpose, we envision a compilation flow involving one or several source-to-source transformations, as illustrated by Figure 1. Before detailing it, we must digress a little on the subtleties of the management of floating-point arithmetics by compilers.

HLS faithful to the floats: Most recent compilers, including the HLS ones [START_REF] Hrica | Floating-point design with vivado HLS[END_REF], attempt to follow established standards, in particular C11 and, for floating-point arithmetic, IEEE-754. This brings the huge advantage of almost bit-exact reproducibility -the hardware will compute exactly the same results as the software. However, it also greatly reduces the freedom of optimization by the compiler. For instance, as floating point addition is not associative, C11 mandates that code written a+b+c+d should be executed as ((a+b)+c)+d, although (a+b)+(c+d) would have shorter latency. This also prevents the parallelization of loops implementing reductions. A reduction is an associative computation which reduces a set of input values into a reduction location. Listing 1 provides the simplest example of reduction, where acc is the reduction location.

The first column of sequential semantics is also possible, but very expensive [START_REF] Kapre | Optimistic parallelization of floatingpoint accumulation[END_REF].

Towards HLS faithful to the reals: Another point of view, chosen in this work, is to assume that the floatingpoint C program is intended to describe a computation on real numbers. In other words, the floats are interpreted as real numbers in the initial C, thus recovering the freedom of associativity (among other). Indeed, most programmers will perform the kind of non-bit-exact optimizations illustrated by Listing 2 (sometimes assisted by source-to-source compilers or "unsafe" compiler optimizations). In a hardware context, we may also assume they wish they could tailor the precision (hence the cost) to the accuracy requirements of the application -a classical concern in HLS [START_REF] Gort | Range and bitmask analysis for hardware optimization in high-level synthesis[END_REF], [START_REF] Caffarena | High-level synthesis of multiple word-length DSP algorithms using heterogeneous-resource FPGAs[END_REF]. In this case, a pragma should specify the accuracy of the computation with respect to the exact result. A high-level compiler is then in charge of determining the best way to ensure the prescribed accuracy.

The proposed approach uses number formats that are larger or smaller than the standard ones. These, and the corresponding operators, are presented in Section II. Then Section III presents and evaluates the compiler side of the proposed technique, and section IV evaluates it on the FPMark benchmark suite.

II. THE ARITHMETIC SIDE: AN APPLICATION-SPECIFIC

ACCUMULATOR IN VIVADOHLS Kulisch advocated a very large floating-point accumulator [START_REF] Kulisch | The exact dot product as basic tool for long interval arithmetic[END_REF] whose 4288 bits would cover the entire range of double precision floating-point. Such an accumulator would remove rounding errors from all the possible floating-point additions and sums of products, with the added bonus that addition would become associative. So far, Kulisch's full accumulator has proven too costly to appear in mainstream processors. However, in the context of application acceleration with FPGAs, it can be tailored to the accuracy requirements of applications. Its cost then becomes comparable to classical floating point operators, although it vastly improves accuracy [START_REF] Florent De Dinechin | An FPGA-specific approach to floating-point accumulation and sumof-products[END_REF]. This operator can be found in the FloPoCo [START_REF] De | High-Performance Computing Using FPGAs[END_REF] generator and in Altera DSP Builder Advanced. Its core idea, illustrated on Figure 2, is to use a large fixed-point register into which the mantissas of incoming floating-point summands are shifted (top) then accumulated (middle). A third component (bottom) converts the content of the accumulator back to the floating-point format. The sub-blocks visible on this figure (shifter, adder, and leading zero counter) are essentially the building blocks of a classical floating-point adder.

The accumulator described in this section presents two improvements over the one offered in FloPoCo [START_REF] Florent De Dinechin | An FPGA-specific approach to floating-point accumulation and sumof-products[END_REF]:

• In FloPoCo, Float-to-Fix and Accumulator form a single component, which restricts its application to simple accumulations similar to Listing 1. The two components of Figure 2 enable a generalization to arbitrary summations within a loop, as Section III will show. • Our implementation supports subnormal numbers.

A. The parameters of a large accumulator

The main feature of this approach is that the internal fixed-point representation is configurable in order to control accuracy. It has two parameters:

• MSB A is the weight of the most significant bit of the accumulator. For example, if MSB A = 20, the accumulator can accommodate values up to a magnitude of 2 20 ≈ 10 6 . • LSB A is the weight of the least significant bit of the accumulator. For example, if LSB A = -50, the accumulator can hold data accurate to 2 -50 ≈ 10 -15 .

The accumulator width w a is then computed as MSB A -LSB A +1, 71 bits in the previous example. 71 bits represents a wide range and high accuracy, and still additions on this format will have one-cycle latency for practical frequencies on recent FPGAs. If this is not enough the frequency can be improved thanks to partial carry save [START_REF] Florent De Dinechin | An FPGA-specific approach to floating-point accumulation and sumof-products[END_REF]

B. Implementation details

This accumulator has been implemented in C, using arbitrary-precision fixed point types (ap_int). The addition is then written +, the shift is written using the C operator <<. The leading zero count, bit range selections and other operations are implemented using VivadoHLS built-in functions. A few VivadoHLS compiler directives are also used to get the best performance from the design. Altogether these components are written as two C functions of 28 lines of code for the FloatToFix, 22 lines for the FixToFloat.

C. Validation

To evaluate and refine this implementation, we used Listing 3, which we compared to Listings 1 and 2. In the latter, the loop was unrolled by a factor 7, as it is the latency of a floating-point adder on our target FPGA (Kintex-7).

For test data, we use as in Muller et al. [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF] the input values c[i]=(float)cos (i), where i is the input array's index. Therefore the accumulation computes

i c[i].
The parameters chosen for the accumulator are:

• MSBA = 17. Indeed, as we are adding cos(i) 100K times, an upper bound is 100K, which can be encoded in 17 bits. • MAXMSBx = 1 as the maximum input value is 1.

• LSBA = -50: the accumulator itself will be accurate to the 50th fractional bit. Note that a float input will see its mantissa rounded by FloatToFix only if its exponent is smaller than 2 -25 , which is very rare. In other words, this accumulator is much more accurate than the data that is thrown to it. The results are reported in Table I for simple and double precision. The Accuracy line of the table reports the number of correct bits of each implementation, after the result has been rounded to a float. All the data in this table was obtained by generating VHDL from C synthesis using VivadoHLS followed by place and route from Vivado v2015.4, build 1412921. This table also reports synthesis results for the corresponding FloPoCo-generated VHDL, which doesn't include the array management. VivadoHLS uses DSPs to implement the shifts in its floating-point adders. Even if the shifts were implemented in LUTs, the first column would remain well below 500 LUTs: it has the best resource usage. However the latency of one iteration is 7 cycles, hence 100K iterations takes 700K cycles. When unrolling the loop, VivadoHLS is using almost 4 times more LUTs for floats, and 3 times more for doubles. The unrolled versions improves latency over naive versions. Nevertheless, our approach gets even betters latencies for a reasonable LUT usage. Also, we achieve maximum accuracy for the float format which caps at 24 bits (the internal representations of the double, unrolled double and our approach have a higher accuracy than 24 bits, but are then casted to the 24 bits of the float format). Finally, our results are very close to FloPoCo ones, both in terms of LUTs usage, DPSs and latency.

Using this implementation method, we also created an exact floating-point multiplier with the final rounding removed as in [START_REF] Florent De Dinechin | An FPGA-specific approach to floating-point accumulation and sumof-products[END_REF]. This function is called ExactProductFloatToFix. Due to lack of space we do not present it in detail. As the output of this multiplier is not standard, we also created an adapted Float-to-fix block called ExactProductFloatToFix. These functions represent 44 lines of code for ExactProduct and 21 lines of code for ExactProductFloatToFix.

III. THE COMPILER SIDE: GECOS SOURCE-TO-SOURCE

TRANSFORMATIONS

Now that we showed that VivadoHLS is able to generate specialized operators of similar quality to FloPoCo's ones, we want to provide a tool that not only transforms Listing 1 into Listing 3, but also generalizes this transformation to many more situations.

We chose to develop this work in GeCoS [START_REF] Floc | GeCoS: A framework for prototyping custom hardware design flows[END_REF], an opensource, extensible source-to-source compiler framework built upon model-driven engineering. Indeed, all the present work is already freely available as a GeCoS plugin.

This work focuses on two computational patterns, namely the accumulation and the sum of product. Both are specific instances of the reduction pattern, which can be optimized by many compilers or parallel run-time environments. Such patterns are therefore exposed to the compiler/runtime either by the user through directives, or automatically identified using static analysis techniques [START_REF] Redon | Detection of scans in the polytope model[END_REF], [START_REF] Doerfert | Polly's polyhedral scheduling in the presence of reductions[END_REF].

Since our focus is not on the detection of reductions, our tool follows a simple approach based on a combination of user directive and (simple) program analysis. More specifically, the user has to identify target accumulation variable through a pragma, and provide additional information such as the dynamic range of the accumulated data along with the target accuracy.

This local approach is easier, more general and less invasive than approaches that attempt to convert a whole floating-point program into fixed-point [START_REF] Sentieys | Automatic Fixed-Point Conversion: a Gateway to High-Level Power Optimization[END_REF].

A. Compiler directive

In imperative languages such as C, reduction are implemented using for or while constructs. Our compiler directive must therefore appear inside such a construct. Listing 4 illustrates its usage on the code of Listing 1.

The pragma must contain the following information:

• The keyword FPacc, which triggers the transformations.

• The name of the variable in which the accumulation is performed, preceded with the keyword VAR. In the example, the accumulation variable is acc.

• The maximum value that can be reached by the accumulator through the use of the MaxAcc keyword. This value is used to determine the weight MSB A . • The desired accuracy of the accumulator using the epsilon keyword. This value is used to determine the weight LSB A . • Optional: The maximum value of the inputs of the accumulator in the MaxInput field. This value is used to determine the weight MaxMSB X . If this information is not provided, then MaxMSB X is set to MSB A .

In the case when no size parameters are given, a full Kulisch accumulator is produced. Also note that the user can quietly overestimate the maximum value of it's accumulator without major impact on area. For instance, overestimating MaxAcc by a factor 10 only adds 3 bits to the accumulator width.

B. Code transformation

The proposed transformation consists in an algorithm that navigates and modifies the intermediate representation (IR) of the program. To illustrate this algorithm, consider the sample program shown in Listing 5. It performs a reduction into the variable sum, involving both sums and sums of product. The IR associated to the loop body is the directed acyclic graph (DAG) of Figure 3a. The keywords FPMul and FPAdd represent the use of floating-point multipliers and adders respectively. The dotted arrows represent the data dependency between two consecutive iterations of the loop. With a floatingpoint adder needing 7 cycles, the next iteration will be stalled during this operation. As the proposed algorithm pushes the floating-point normalization out of the loop, the adder only have a 1-cycle delay. Therefore we are able to have an iteration iterval of 1 as shown in Figure 4a.

; i++){ sum+=in1[i] * in2[i-1]; sum+=in1[i]; sum+=in2[i+1]; } return sum; }
The code transformation is applied to every block within the for block annotated with our pragma. It is a bottom-up walk of the program DAG, starting from a FPAdd node that writes to the accumulation variable. During this walk, the following actions are performed depending on the visited nodes:

• A node with the summation variable is ignored.

• A FPAdd node is transformed to an accurate fixed-point adder. The analysis is then recursively launched on that node. • A FPMul node is replaced with a call to the ExactProduct function followed by a call to ExactProdFloatToFix. • Any other node has a call to FloatToFix inserted. This algorithm rewrites the DAG from Figure 3a into the new DAG shown on Figure 4a.

in2[i-1] in1[i] float mul float add in2[i+1] float add float add sum sum (a) DAG in2[i-1] in1[i] in2[i+1]
Then the corresponding C code is generated. In addintion, we provide a statement on the accumualation loop that tries to reduce the iteration interval to 1. The synthesized codes from before and after the transformations result in the architectures from Figure 3b and Figure 4b respectively. Also, at the end of the transformed loop, a call to FixToFloat is inserted in order for the programmer to retrieve his value as a floatingpoint.

C. Evaluation of the toy example of Listing 5

The proposed transformations work on non-trivial examples such as Listing 5. depends on epsilon, all the other parameters being those given in the pragma of Listing 5. All these versions where synthesised for 100MHz, but note that they all can run at 400MHZ. Our transformed code makes VivadoHLS use more LUTs for less DSPs. Again this is due to Vivado's use of DSP for the shifters. In all cases, on this example, the transformed code has its latency reduced by a factor 20.

IV. EVALUATION ON FPMARK BENCHMARKS

In order to evaluate the relevance of the proposed transformations on real-life programs, we used the EEMBC FPMark benchmark suite [START_REF]EEMBC -the embedded microprocessor benchmark consortium[END_REF]. Table II: Comparison between the naive code from Listing 5 and its transformed equivalent. All these versions run at 100MHz.

This suite consists of 10 programs. A first result is that half of these programs contain visible accumulations:

• Enhanced Livermore Loops (1/16 kernels contains one accumulation) • LU Decomposition (multiple accumulations)

• Neural Net (multiple accumulations)

• Fourier Coefficients (one accumulation)

• Black Scholes (one accumulation)

The following focuses on these, and ignores the other half (Fast Fourier Transform, Horner's method, Linpack, ArcTan, Ray Tracer).

Most benchmarks come in single-precision and doubleprecision versions, and we focus here on the single-precision ones.

A. Benchmarks and accuracy: methodology

Each benchmark comes with a golden reference against which the computed results are compared. As the proposed transformations are controlled by the accuracy, it may happen that the transformed benchmark is less accurate than the original. In this case, it will not pass the benchmark verification test, and rightly so.

A problem is that the transformed code will also fail the test if it is more accurate than the original. Indeed, the golden reference is the result of a certain combination of rounding errors using the standard FP formats, which we do not attempt to replicate.

To work around this problem, each benchmark was first transformed into a high-precision version where the accumulation variable is a 10,000-bit floating-point numbers using the MPFR library. We used the result of this highly-accurate version as a "platinum" reference, against which we could measure the accuracy of the benchmark's golden reference. This allowed us to choose our epsilon parameter such that the transformed code would be at least as accurate as the golden reference.

B. Benchmarks improved by the proposed transformation

Enhanced Livermore Loops: This program contains 16 kernels of loops that compute physics equations. Among these kernels, there is one that performs a sum-of-product (banded linear equations). This kernel computes 20000 sumsof-products. The values accumulated are taken from an array. This is a perfect candidate for the proposed transformations.

For this benchmark, the optimal accumulation parameters were found as: MaxAcc=50000 epsilon=1e-5 MaxInput=22000

Synthesis results of both codes (before and after transformation) are given in Table III Both of these programs accumulate values from registers or memory that are already computed. It makes these programs good candidates for the proposed transformations.

Due to lack of time, and also because VivadoHLS is unable to predict a latency for these implemented designs due to their variable loop sizes, we do not present complete results for these two benchmarks. Still, in order to show that the approach works on these examples, the LU inner loops were transformed and synthesized. Table IV shows the results obtained for the smallest (8 terms) and the largest (45 terms) sums-of-products: the latency is vastly improved even for the smallest one. The Fourier coefficients program contains an accumulation which is performed in simple precision. This program comes in three different configurations: small, medium and big. Each of them computes the same algorithm but with a different amount of iterations. The "big" version is supposed to compute the most accurate answer. We get similar results for the three versions of this program, so we only present the "big" version here. In this version, 2K terms are accumulated multiple times during the computation. The accumulator is reset at every call.

The parameters determined for this benchmark were the following: MaxAcc=6000 epsilon=1e-7 MaxInput=10

This results in an accumulator using 14 bits for the integer part and 24 bits for the fractional part. The synthesis results obtained for the original and transformed codes are given in Table V. Here, area is again comparable, accuracy is improved by one order of magnitude, but latency is not improved. To understand why, Listing 6 provides an example illustrating the problem. The variable x from statement 1 (S1) is used as a function parameter in S2. This introduces an intra-iteration dependency between these statements. As we transform S2 to be an application-specific accumulator with a 1 cycle latency, we still have to wait for x to be ready at each iteration. This makes the overall latency of one iteration to be limited by the speed at which x can be computed.

As x is itself an accumulation, an idea is to apply the transformation to x as well. The problem is that S2 needs its value as a float at each iteration. If x is computed as a non-standard fixed point number, it will need to be converted to a float at each iteration, which takes several cycles.

In principle there is a solution, which is to decouple the two accumulation loops and insert a fix-to-float pipeline between them. However this transformation is well beyond the scope of this article and is left for future work. Black Scholes: This program contains an accumulation that sums 200 terms. The result of this computation is divided by a contant (that could be optimized by using transformations based on [START_REF] De | Table-Based Division by Small Integer Constants[END_REF]). This process is performed 5000 times.

Here the optimal accumulator parameters are the following: MaxAcc=245000 epsilon=1e-4 MaxInput=278 This gives us an accumulator that uses 19 bits for the integer part and 10 bits for the fractional part. The result of the synthesis are provided in Table VI Again, for comparable area, accuracy is vastly improved by latency is not improved. This time it is for a different reason, illustrated by Listing 7. Let us consider the latency of the floating-point accumulation to be N , and the latency of the function f to be M . The latency of one iteration is then going to be max(N, M ). In the case where M < N our transformations will improve the overall latency by reducing it to M . Listing 7: Illustration of a computation latency preventing a latency reduction float acc=0; for(int i=0; i<N; i++){ acc += f(...); } However, in the context of Black Sholes, M is greater than N . Therefore, reducing the latency of the accumulation has no impact on performance, as we have to wait for the result of the computation of f at each iteration. Again, there is an architectural solution to this, but VivadoHLS is currently unable to find it.

V. CONCLUSION

The main result of this work is that HLS tools have the potential to generate efficient designs for handling floatingpoint computations in a completely non-standard way. The use of application-specific intermediate formats can provide both performance and accuracy at a competitive cost. For this, we have to sacrifice the strict respect of the IEEE-754 and C11 standards. It is replaced by the strict respect of a high-level accuracy specification.

Classically, designers have to face a trade-off between performance and cost. This approach adds computation accuracy to this trade-off. Some designers may not like this. To convince them, consider that established performance benchmarks compute results which are accurate only to a few bits. If only a few bits are important, do we really need to instanciate 32bit or 64-bit floating-point operators to compute them ? Isn't this accuracy information worth investigating and exploiting?

This work also provides a practical tool that improves a given C program. The input to the tool is application-specific information representing high-level domain knowledge such as the range and desired accuracy of a variable. The resulting code is compatible with VivadoHLS.

The proposed transformation already works very well on 3 of the 10 FPMarks where it improves both latency and accuracy by an order of magnitude for comparable area. For 2 more benchmarks, the latency is not improved (but not degraded either) due to current limitations of HLS tools. This defines short-term future work, which could be addressed either within the back-end HLS tool, or as more source-tosource transformations.

In the longer term, we believe there is much more to come. The arithmetic optimizations that a classical compiler can do are very limited by the fixed hardware of classical processors. With compilers of high-level software to hardware, there is much more freedom, hence many more opportunities to build application-specific arithmetic operators. Future work will attempt to explore this new realm, starting with operator specialisation, operator fusion, and compile-time generation of application-specific cores, and building upon compiler progresses in program analysis.

Figure 1 :

 1 Figure 1: The proposed compilation flow

Figure 2 :

 2 Figure 2: The conversion from float to fixed-point (top), the fixed-point accumulation (middle) and the conversion from the fixed-point format to a float (bottom).

Listing 3 :

 3 Sum of floats using the large fixed-point accumulator #define N 100000 float acc = 0; ap_int<68> long_accumulator = 0; for(int i = 0; i < N; i++) { long_accumulator += FloatToFix(in[i]); } acc = FixToFloat(long_accumulator);

Listing 4 :

 4 Illustration of the use of a pragma for the naive accumulation #define N 100000 float accumulation(float in[N]){ float acc = 0; #pragma FPacc VAR=acc MaxAcc=100000 epsilon=1E-15 MaxInput=1 for(int i=0; i<N; i++){ acc+=in[i]; } return acc; }

Listing 5 :

 5 Simple reduction with multiple accumulation statements #define N 100000 float computeSum(float in1[N], float in2[N]){ float sum = 0; #pragma FPacc VAR=sum MaxAcc=300000 epsilon=1e-15 MaxInput=3 for (int i=1; i<N-1

Figure 3 :

 3 Figure 3: DAG of the loop body from Listing 5 (top) and it's corresponding architecture (bottom)

Figure 4 :

 4 Figure 4: DAG of the loop body from Listing 5 (top) and it's corresponding architecture (bottom) after transformations

Listing 6 :

 6 Illustration of a loop intra-iteration dependency preventing a latency reduction float x=0, acc=0; for(int i=0; i<N; i++){ S1: x += dx; S2: acc += f(...,x,...); }

Table I

 I

	Listing 1: Naive reduction
	#define N 100000
	float acc = 0;
	for(int i=0; i<N; i++){
	acc+=in[i];
	}
	Listing 2: Unrolled reduction
	#define N 100000
	float acc = 0, tmp1=0, ... , tmp10=0;
	for(int i=0; i<N; i+=10){
	tmp1+=in[i];
	...
	tmp10+=in[i+9];
	}
	acc=tmp1+...+tmp10;
	shows how VivadoHLS synthe-
	sizes Listing 1 on Kintex7. The floating-point addition takes
	7 cycles, and the adder is only active one cycle out of 7 due
	to the loop-carried dependency. Listing 2 shows a manually
	unrolled version of Listing 1. VivadoHLS will not transform
	Listing 1 into Listing 2, because they are not semantically
	equivalent (the floating-point additions are reordered as if
	they were associative). However Listing 2 expresses more
	parallelism, which VivadoHLS is able to exploit (second
	column of Table I). The main adder is now active at each cycle
	on a different sub-sum. Note that a parallel execution with the

Table I :

 I Different ways of implementing a simple accumulation.

		Listing 1	Listing 2	Listing 1	Listing 2	Listing 3	FloPoCo VHDL
		(float)	(float)	(double)	(double)	(71 bits)	(71 bits)
	LUTs	266	907	801	2193	736	719
	DSPs	2	4	3	6	0	0
	Latency	700K	142K	700K	142K	100K	100K
	Accuracy	17 bits	17 bits	24 bits	24 bits	24 bits	24 bits

  but this was not useful in the present work. For comparison, for the same frequency, a floating-point adder has a latency of 7 to 10 cycles, depending on the target.

	Exponent Mantissa	Sign
	MaxMSB X	w e	w f
	-
	Shifter Negate MaxMSB X -LSB A + 1 Shift value	FloatToFix
	w A	+ Registers		Accumulator
		Negate
	LZC + Shifter w A w e w f		FixToFloat
	Exponent Mantissa	Sign

Table II

 II

	in2[i-1]	in1[i]
	Exact Product
	ExactProduct
	FloatToFix
			FloatToFix
	in2[i+1]
			Fixed Acc
	FloatToFix
	long accumulator	FixedAcc
	Fixed Acc
	long accumulator
		(a) DAG
	in2[i+1]	in2[i-1]	in1[i]
	Exact Product
	ExactProduct
	FloatToFix
			FloatToFix
	FloatToFix
		Fixed Add
		Fixed Add
	Fixed Add
	long accumulator
	(b) Architecture
				shows how resource consumption

Table III :

 III . As in the previous toy examples, latency and accuracy are vastly improved for comparable area. Synthesis results of one Livermore Loops kernel before and after transformations.LU Decomposition and Neural Net: Both the LU decomposition and the neural net programs contain multiple nested small accumulations. In the LU decomposition program, an inner loop accumulates between 8 and 45 values. Such accumulations are performed more than 7M times. In the neural net program, inner loops accumulate between 8 and 35 values, and such accumulations are performed more than 5K times.

		Latency LUTs DSPs	Accuracy
	Original Code	80K	384	5	11 bits
	Transformed Code	20K	576	2	13 bits

Table IV :

 IV Synthesis results of the loops from LU benchmark.

			Latency LUTs	DSPs
	8 terms	Original Code	82	809	5
		Transformed Code	17	1007	2
	45 terms	Original Code	452	819	5
		Transformed Code	54	1034	2
	C. Benchmarks that suggest future work
	Fourier Coefficients:

Table V :

 V Synthesis results of the Fourier coefficients program before and after transformations.

		Latency	LUTs	DSPs	Accuracy
	Original Code	8K	34596	64	6 bits
	Transformed Code	8K	34681	59	11 bits

  .

		Latency	LUTs	DSPs	Accuracy
	Original Code	3M	15640	175	19 bits
	Transformed Code	3M	15923	175	23 bits
	Table VI: Synthesis results of the Black Scholes program
	before and after transformations. The last column provides
	the results for a hand tuned version where we split the loop
	in two.

Catapult C Synthesis, Mentor Graphics,

2011, http://calypto.com/en/ products/catapult/overview/ 2 Arbitrary-size floating-point should follow some day, it is well supported by mature libraries and tools


Rui Zengh
email: zengrui@seu.edu.cn.jswu@seu.edu.cn.lotfi.senhadji@univ-rennesi.fr.shu.list@seu.edu.cn

Jiasong Wu

Senhadjitt- Lotfi

Huazhong Shuh

TENSOR OBJECT CLASSIFICATION VIA MULTILINEAR DISCRIMINANT ANALYSIS NETWORK

Keywords: Deep learning, MLDANet, PCANet, L DANet, tensor object classification Recently, Chan et al. [5]

This paper proposes an multilinear discriminant analysis net work (M LDANet) for the recognition of multidimensional objects, knows as tensor objects. The MLDANet is a variation of linear discriminant analysis network (LDANet) and prin cipal component analysis network (PCANet), both of which are the recently proposed deep learning algorithms. The ML DANet consists of three parts: 1) The encoder learned by MLDA from tensor data. 2) Features maps obtained from de coder.

3) The use of binary hashing and histogram for feature pooling. A learning algorithm for MLDANet is described. Evaluations on UCFll database indicate that the proposed MLDANet outperforms the PCANet, LDANet, MPCA+LDA, and MLDA in terms of classification for tensor objects.

INTRODUCTION

One key ingredient for the success of deep learning in visual content classification is the utilization of convolution archi tectures [START_REF] Bengio | Learnin .,.g deep architectures for AI[END_REF][START_REF] Lecun | Gradient-based learning applied to document recogni tion[END_REF][START_REF] Krizhevsky | Image net classification with deep convolutional neural networks[END_REF], which are inspired by the structure of human visual system [START_REF] Lecun | Convolutional networks for images, speech, and time series[END_REF]. A convolution neural network (CNN) [START_REF] Lecun | Gradient-based learning applied to document recogni tion[END_REF] consists of multiple trainable stages stacked on the top of each other, following a supervised classifier. Each stage of CNN is organized in two layers: convolution layer and pooling layer. the dictionary in the convolution layer and the pooling lay er is composed of the simplest binary hashing and histogram. The PCANet leads to some pleasant and thought-provoking surprises: such a basic network has achieved the state-of-the art performance in many visual content datasets. Meanwhile, Chan et al. [START_REF] Chan | P canet: A simple deep learning baseline for image classi fication?[END_REF] proposed linear discriminant analysis network (LDANet) as a variation of PCANet.

However, PCANet and LDANet are deteriorated when dealing with visual content, which is naturally represent ed as tensor objects. That is because when using PCANet or LDANet, the multidimensional patches, taken from visual content, are simply converted to vector to learn the dictionary. It is well known that vector representation of patches breaks the natural structure and correlation in the original visual content. Moreover, it may also, suffers from the so-called curse of dimensionality [START_REF] Shakhnarovich | Face recogni tion in subspaces[END_REF].

Recently, there is growing interest in the tensorial exten sion of deep learning algorithms. Yu et al. [START_REF] Yu | The deep tensor neu ral network with applications to large vocabulary speech recognition[END_REF] proposed deep tensor neural network, which can be seen as a tensorial ex tension of deep neural network (DNN), outperforms DNN in large vocabulary speech recognition. Hutchinson et al. [START_REF] Hutchinson | Tensor deep s tacking networks[END_REF] p resented the tensorial extension of deep stack neural network, which has been successfully used in MNIST handwriting im age recognition, phone classification, etc. However, the simi lar tensorial extension research has not been reported for deep learning algorithms with convolutional architecture.

In this paper, we propose a simple deep learning algorith m for tensor object classification, that is, multilinear discrim inant analysis network, which is a tensorial extension of P CANet and LDANet. The simulation on UCF 11 database [START_REF]UCF11[END_REF] demonstrates that the MLDANet outperforms PCANet and L DANet in terms of classification accuracy for tensor objects.

REVIEW OF MLDA

In this section, we briefly review MLDA [10], which is a mul tilinear extension of LDA. The MLDA obtains discriminative features through maximizing the Fisher's discrimination cri terion, which is described as follows.

An Nth tensor object is denoted as X E R / l X I 2 X ... X I N .

It is addressed by N indices i n , n = 1,2, ... ,N, and each i n addresses the n-mode of X. The n-mode tensor product of X by a matrix U E n J " xI " is defined as: ( X x."U)(i" .. ,i"_l,j,,,i,,+l, .. ,iN) = L X (i1, .. ,iN rU(j,,,i,,)• (1) The projection from tensor X E n I t xI 2 x•••xI N to a s calar Y can be described as follows:

(2)

where { u ( n ) T }t:=1 is a set of unit projection vectors. This ten sor to scalar projection is called elementary multilinear pro jection (EMP), which consists of one projection vector in each mode. An EMP of a tensor X E R /,XI 2 XI 3 is illustrated in Fig. 1.

The tensor-to-vector projection (TVP) from a tensor X E n I t xI 2 x ... xI N to a vector Y E R P is to find a vector set { u � n ) 1', n = 1, ... , N}:=l' which are able to do P times EM P. The process can be described as:

-X N { ( n ) 1' -N} P Y -X n =l u p ,n -1, ... , p=l , ) 1' { u � n ,n = 1, ... ,N } = arg max F;

(5)

THE ARCHITECTURE OF MLDANET

Fig. 3 shows the architecture of MLDANet for third-order tensor objects classification. It contains two convolutional layers and one pooling layer. The filter bank in each convolu tional layer is learned independently. We use binary hashing and histogram as pooling operation for the features extracted from the first two convolutional layers.

1. The first stage of MLDANet

For the given M third-order tensor objects {Xm}�=l E R / l X I t X I 3 , which contains C classes, we take tensor patch es around the I-mode and 2-mode by taking all 3-mode elements of mth tensor object, i.e., the tensor patch size is k1 x k2 X 1 3 , we collect all (overlapping) It x 12 ten sor patches from Xm. The tensor patches have the same class with Xm. We put these tensor patches into a set tm = { tm ,q E n k , x k 2 X I s } ��� I 2 . By repeating the above process for every tensor objects, we can get all tensor patches t = { t m}�=l for learning filter bank in the first stage. Let L1 be the number of filters in the first stage. We apply

1 L { (I) ( 2 ) ( 3 ) } L l MLDA to t to
F ml = mat(t x3 _ { u ( n ) T }) E n J,xI2 m,q n -1 I , q = 1, ... ,h X 12, (6)
where mat( v ) is a function that maps v E n 1 , h on a matrix FEnJ ,xI 2.

For each tensor object, we can obtain L1 feature map s of size h x 12. We denote these feature maps of mth tensor object in the first stage with {F ml E nJ ,xI 2,1 = 1, ... , L1};;; = 1. The feature maps of each tensor object cap ture the main variation of original data.

The second stage of MLDANet

Through the first stage, the tensor object is already mapped into low-dimensional tensor feature, the dimension of the 3mode is much lower than that of the I-mode and 2-mode. That is to say, the redundancy of 3-mode has been greatly reduced. Therefore, for the simplicity of computation and the conve nience of building network, we use the conventional LDA in the second stage to learn the filter bank. The number of filters in the second stage is L2.

Around each pixel, we take a k1 x k2 patch, and col lect all (overlapping) patches of all the feature maps F ml, .

{ } J,Xh nk ,x k 2 h h d t th

I .e., rml,q q =l E l'v we re eac rml,q eno es e qth vectorized patch in the lth feature map of mth tensor ob ject. We then subtract the patch mean from each patch, and construct the matrix R ml = [rml,l, rml,2,• .. , rml,I , XI 2 ] for them, where R ml belongs to the same class with the mth ten sor object. Let S e is the set of matrix R ml in class c. We then compute the class mean f e and the mean of class f as follows:

So, the within-class scatter matrix and the between-class scatter matrix are defined respectively as follows:

C m 1 h = R ml * v h,h = 1, ... ,L 2,1 = 1, ... ,L 1, (11)
where * denotes 2D convolution [START_REF] Lecun | Gradient-based learning applied to document recogni tion[END_REF], and the boundary of R ml is zero-padded before convolving with V h so as to make Cm 1 h have the same size as R ml. The number of output feature maps of the second stage is L1L2. One or more additional stages can be built if a deep architecture is found to be bene ficial.

The pooling layer in MLDANet

First, we binarize each feature map by using Heaviside step function H ( . ), whose value is one for positive entries and zero otherwise. The binarized feature maps are denoted by C m 1 h E { O , 1 }J, X 12. Owing to every feature map capture d ifferent variations by Vh. Cmlh should be weighted to convert into a single integer-valued feature map:

L2 W ml = L 2 h -1 cmlh . (12)
h =l note that each entry of feature map w ml are integers in the range [0,2L2 -1].

Next, we partition w ml into B blocks, and then compute the histogram (with 2L2 bins) of the decimal values in each block. All the B histograms are concatenated into one vector as the lth feature vector Bhis t (w ml ) of tensor object X m. The final feature of input tensor object Xm is then defined as the set of feature vector, i.e., fm=[Bhis t (w m1 ), Bhis t (W m2 ), ... , Bhis t (W m L, )]. (13) Note that the local block can be either overlapping or non overlapping depending on applications [START_REF] Chan | P canet: A simple deep learning baseline for image classi fication?[END_REF].

EXPERIMENTAL RESULTS

We evaluate the performance ofMLDANet on UeFll dataset for tensor object classification. UeFll is a sport action video dataset which contains 11 action categories: basketball shoot ing, biking, diving, golf swinging, horseback riding, soccer juggling, swinging, tennis swing, trampoline jumping, vol leyball spiking, and walking with a dog. All videos in UeFll are manually collected from YouTube and their sizes are all 240 x 320 pixels. For each category, the videos are grouped into 25 groups with more than 4 action clips in it. The video clips in the same group have COlmnon scenario. This video dataset is very challenging in classification due to large vari ations in camera motion, object appearance and pose, object scale, viewpoint, cluttered background, illumination, etc.

In this experiment, we only choose the first ten groups in each category. The total number of experimental videos is For each group, half videos are randomly selected for training and others for testing. Every video is resized to 48 x 64 in order to reduce the computational complexity. Almost every videos have variation in frames. For those frame larger than 20, we only choose the first twenty frames. For a few videos, whose frames are less than 20, we just copy the last frame to fill them.

We then compare the proposed MLDANet with PCANet [START_REF] Chan | P canet: A simple deep learning baseline for image classi fication?[END_REF], LDANet [START_REF] Chan | P canet: A simple deep learning baseline for image classi fication?[END_REF], MPCA+LDA [START_REF] Lu | M pca: Multilinear principal component analysis of tensor objects[END_REF], and MLDA [10]. The model parameters of MLDANet, PCANet, and LDANet all include the patch size k1 x k2, the number of filters in each stage L1, L2, the number of stages, overlapping ratio of block, and the block size. Chan et al. [START_REF] Chan | P canet: A simple deep learning baseline for image classi fication?[END_REF] have shown that the ap propriate number of filters L1, L2 in PCANet and LDANet is L1 = L2 = 8. By considering that MLDA is the ten sorial extension of conventional LDA. Thus, we always set L1 = L2 = 8 for all networks. The patch size k1 x k2 are changed from 3 x 3 to 7 x 7 and three block sizes 6 x 8, 12 x 16, 24 x 32 are considered here. The overlapping ratio is set to 50%. Unless stated otherwise, we use linear SVM classifier. The recognition rates of above networks averaged over 5 different random splits are shown in Fig. 4 (a)-(c).

For conventional tensor object classification by using MP CA+LDA, we change the dimensions of input feature of LDA from 10 to 100. The dimensions of feature vector extract- ed from MLDA vary from 10 to 100. We draw the recog nition accuracy of MPCA+LDA and MLDA in Fig. 4(d).

The best performance of MLDANet, PCANet, LDANet, MP CA+LDA, and MLDA are listed in Table 1.

We see that all one-stage networks outperform two con ventional tensor object classification algorithms, that is, M PCA+LDA and MLDA. The reason is that the convolution al architecture imitate the brain networks, which can provide more robust feature than other methods for visual content [START_REF] Lecun | Convolutional networks for images, speech, and time series[END_REF]. LDANet-l achieves the best performance in the one-stage networks, but the improvement from LDANet-l to LDANet-2 is not larger as that ofMLDANet. PCANet-l performs worse than those based on LDA algorithm networks like LDANet and MLDANet. It is because that LDA type algorithm pro vides the features which have the best classification perfor mance, however, PCA maximize the directional variation in the features. MLDANet-l is not as good as LDANet-l be cause the feature extracted from MLDANet-l is not appropri ate as the direct input of linear SVM. For two-stage network s, MLDANet-2 achieves the best performance. Surprisingly, the performance of PCANet-2 increase not more than 18.24% compared to that of PCANet-l, but it is better than that of LDANet-2.

CONCLUSION

In this paper, we have proposed and implemented a novel deep learning architecture, that is, MLDANet, which takes full advantage of the structure information in tensor objects by convolutional architecture. MLDANet is composed of t wo convolutional layers, which use MLDA and LDA to learn filter banks respectively and one pooling layer. We have e valuated the performance of the MLDANet on UCFll and show that our model performs well in tensor object classifica tion. This work provides the inspiration for other convolution al deep architectures in tensor object classification. As future works, we will focus on the tensorial extension of CNN.

Fig. ! .

 ! Fig. !. The process of elementary multilinear projection (EM P).

Fig. 2 .

 2 Fig. 2. Tensor-to-vector projection (TVP).

( 3 )

 3 whose pth component is obtained from the pth EMP as y ( p ) = ( l ) T ( 2 ) T ( N ) T . . X X 1 u p X 2U p ... X N Up • Flg. 2 shows the schematlc plot for TVP. Suppose that we are given M input tensor objects, {Xm}�=l E R / I X ... xI N , which contains C classes. The pth projected scalar of {Xm}�=l are defined as { y m p }�=l ' where Y m p = Xm Xt:=l { u � n ) T }t:=1' So the between-class s catter matrix and the within-class scatter matrix for pth scalar

Fig. 3 .

 3 Fig. 3. The architecture of two-stage MLDANet.

  earn the 1 vector sets u p , u p , u p p=l ' Each tensor patch is converted into Ll scalars by using {U ( l ) U ( 2 ) U ( 3 ) }L, Thus the lth feature map of tensor p , p , P p =l• , object Xm in the first stage is defined as:

  We then get w* E n k , k 2 X L2 by maximizing the Fisher's discriminant criterion as follows: * Tr( wTS Bw) w = arg max (T S ) ' s.t. wT w = h 2 • (10) w Tr w WW By using mat function, each column of w* is converted into matrix { v hE n k , x k 2 } t� 1 . These matrices are treated as the filter bank in the second stage. Let the hth output of the lth feature map for the mth tensor object in the second stage be

  PCANeH 20 '-=:!C: :..!P :l<, CA fl!: N !!i et i; ;o -2 '--____ _

Fig. 4 .

 4 Fig. 4. (a)-( c): Recognition rate of networks in different patch size. (d) are the performance of MPCA+LDA and MLDA

Table 1 .

 1 The best performance of MLDANet, LDANet, P CANet, MPCA+LDA and MLDA.

	Methods	Accuracy
	MLDANet-l	64.55
	MLDANet-2	78.93
	LDANet-l	73.58
	LDANet-2	76.59
	PCANet-l	58.68
	PCANet-2	76.92
	MPCA+LDA	45.15
	MLDA	38.46

This work was supported by the National Basic Research Program of China under Grant 201lCB707904, by the NSFC under Grants 61201344, 61271312, 11301074, by the SRFDP under Grants 20110092110023 and 20120092120036, the Project-ponsored by SRF for ROCS, SEM, by Nat ural Science Foundation of Jiangsu Province under Grant BK2012329 and by Qing Lan Project. This work is supported by INSERM postdoctoral fel lowship.


Majd Saleh

Ahmad Karfoul

Amar Kachenoura

Lotfi Senhadji

Laurent Albera

may come

Low cost and efficient kurtosis-based deflationary ICA

method: application to MRS sources separation problem.

Low Cost and Efficient Kurtosis-Based Deflationary ICA Method:

Application to MRS Sources Separation Problem M. Saleh * , A. Karfoul † ‡ , A. Kachenoura † ‡ , L. Senhadji † ‡ and L. Albera † ‡

I. INTRODUCTION

Independent Component Analysis (ICA) consists of estimating statistically independent source signals from a mixture of them. This mixture can be linearly modeled as follows:

x[m] = As[m] + ν[m] (1)
where m is a natural number and x[m] is a realization of an N -dimensional random vector process {x[m]} of the observations, A (N × P ) denotes a tall mixing matrix, s[m] stands for a realization of a P -dimensional random process of sources {s[m]} with components assumed to be statistically independent and ν[m] denotes a realization of an N -dimensional process of noise assumed to be independent from the source. Note that all processes considered in this paper are real-valued. ICA algorithms estimate sources by maximizing a measure of nongaussianity such as the kurtosis (i.e. Fourth order cumulants) of the data [START_REF] Cardoso | Jacobi angles for simultaneous diagonalization[END_REF], [START_REF] Comon | Independent component analysis -a new concept[END_REF], [START_REF] Delfosse | Adaptive blind separation of independent sources: a deflation approach[END_REF], [START_REF] Hyvärinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF], [START_REF] Oja | The fastica algorithm revisited: convergence analysis[END_REF], [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF] to cite a few. Some ICA approaches [START_REF] Cardoso | Jacobi angles for simultaneous diagonalization[END_REF], [START_REF] Oja | The fastica algorithm revisited: convergence analysis[END_REF] recover independent sources in only one single-step, these are the so-called symmetric approaches. Others, referred to as deflationary approaches extract one source at a time [START_REF] Comon | Independent component analysis -a new concept[END_REF], [START_REF] Delfosse | Adaptive blind separation of independent sources: a deflation approach[END_REF], [START_REF] Hyvärinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF], [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF]. The key advantage of the deflationary scheme is that contrasts for single-source extraction can easily be proven to be free from spurious local solutions, so that all their local optima are associated with valid solutions [START_REF] Comon | Handbook of Blind Source Separation[END_REF]- [START_REF] Hyvärinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF], [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF]. Besides, deflationary scheme can be stopped early once the sources of interest has already been extracted. Hence, a considerable reduction in processing time is guaranteed [START_REF] Sardouie | An efficient jacobi-like deflationary ICA algorithm: Application to EEG denoising[END_REF].

RobustICA algorithm is a common deflationary kurtosisbased ICA algorithm. It performs kurtosis maximization in an efficient way using a gradient-ascent algorithm with global plane search scheme [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF]. However both the execution time and the numerical complexity of the RobustICA could be further improved while maintaining a good estimation quality especially for situations where both the previous issues are crucial. For example, in the case of high dimensional problems often encountered in practice such as in High Resolution-ElectroEncephaloGraphic (HR-EEG) data [START_REF] Oostenveld | The five percent electrode system for high-resolution EEG and ERP measurements[END_REF], in the context of Magnetic Resonance Spectroscopy (MRS) [START_REF] Jansen | 1h mr spectroscopy of the brain: absolute quantification of metabolites[END_REF], etc.

The purpose of this paper is to investigate the possible ways of improving both the numerical complexity and the execution time of the RobustICA approach by i) comparing the relative speed of different optimization schemes used in the RobustICA method and ii) investigating the strategies to accelerate the optimization process (by assessing the validity of these strategies). With these aims in mind, a natural extension of the original RobustICA is first given by resorting to the first order optimization scheme, namely the nonlinear Conjugate Gradient (CG) method. Next, beyond this first order implementation, two second order-based approaches are proposed. The Newton method wherein the Hessian is exactly computed and the quasi-Newton using the BFGS (Named after its inventors, Broyden, Fletcher, Goldfarb, and Shanno) scheme where an approximation of the Hessian is considered instead. The efficiency of the proposed approaches is evaluated in the context of the MRS sources separation and compared to both the conventional RobustICA and the deflation-based FastICA [START_REF] Hyvärinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF] methods.

II. ALGORITHMS

Both the first order and the second order-based RobustICA approaches are considered in this section. Regarding the Newton-based variant, an exact computation of the Hessian is also provided.

A. First order-based optimization scheme (beyond the classical RobustICA implementation)

Basically, the RobustICA algorithm solves the real-valued ICA problem (1) by maximizing the following contrast function [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF]:

F (w (p) ) = |C 4 (y p )| E[y 2 p ] 2 = | E[y 4 p ] E[y 2 p ] 2 -3| (2)
where C 4 (y p ) denotes the fourth-order cumulant (kurtosis) of the p-th estimated source, y p , such that y p = w (p) T x where w (p) is the associated p-the separation vector. Besides, E[.] stands for the mathematical expectation and |.| denotes the absolute value of its argument. Although data pre-whitening is not mandatory for RobustICA, employing it would improve the convergence especially with high-dimensional data [START_REF] Hyvarinen | Independent Component Analysis[END_REF]. Therefore, the observations x are considered to be pre-whitened in the sequel. The considered optimization problem is then defined as:

max w (p) F (w (p) ) = |E[y 4 p ] -3| s. t. w (p) = 1 (3)
Originally, a first-order optimization scheme is adopted to maximize F (2) with the following update rule [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF]:

w (p) k+1 ← w (p) k + µ (k) opt g (p) k (4) where g (p)
k is typically the gradient direction at the current iteration k associated to the p-th source and µ (k) opt is the optimal step size associated to the considered direction g (p)

k . Note that from now on and for the sake of clarity, the subscript p will be omitted in the rest of this paper. It is noteworthy that an optimal computation of µ showed i) a very low computational complexity and ii) a robustness to the presence of spurious local extrema and saddle points in the contrast F (2) [START_REF] Comon | Handbook of Blind Source Separation[END_REF], [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF]. However, the gradient-ascent algorithm often finds itself taking steps in the same direction as the earlier steps. Hence more iterations are needed for convergence and consequently a considerable numerical complexity. A natural way to improve the convergence speed of the gradient-ascent based RobustICA is to use the nonlinear Conjugate Gradient (CG) instead. Indeed, computing nonlinear CG directions has the same complexity as the gradient descent/ascent directions, but in general, they are more effective [START_REF] Nocedal | Numerical Optimization[END_REF]. The CG-based RobustICA algorithm, denoted by RobustICA CG , based on the Polak and Ribière scheme, is summarized as follows:

1) Initialization: random initialization of w 0 , k ← 0 and g 0 ← ∇ w0 F , where ∇ w0 F is the gradient of F with respect to w evaluated at w = w 0 . 2) Repeat until convergence or a maximum number of iterations is reached a) Calculate the optimal step size µ opt in the search direction g k as in [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF]. b) Update w k+1 ← w k + µ opt g k . c) Normalize w k+1 ← w k+1 / w k+1 . d) Calculate ∇ w k+1 F as given by [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF].

e) Calculate β P R k+1 = ∇ w k+1 F T (∇ w k+1 F -∇ w k F ) ∇ w k F 2 . f) Update the search direction: g k+1 ← ∇ w k+1 F + β P R k+1 g k . g) k ← k + 1.
3) Deflate the observations using the Gram Schmidt approach [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF] and go back tp step [START_REF] Boyd | Convex optimization[END_REF].

Note that deflation step can be also performed using a linear regression strategy [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF]. For the sake of clarity, only the Gram Schmidt approach is considered in this paper. CG method typically converges in few iterations thanks to the conjugate direction strategy and shows very low numerical complexity and hence small memory allocation compared to its gradient ascent counterpart. However, the convergence speed of the RobustICA is expected to be further improved when second order information about the objective function is incorporated. In this respect, a second-order optimization scheme such as the Newton [START_REF] Boyd | Convex optimization[END_REF] could be adopted giving rise to the Newton-based RobustICA.

B. Second order optimization scheme

Newton method exploits further information about the surface of the objective function such as the second order derivative (i.e. the Hessian) which is a good way to accelerate the local convergence [START_REF] Madsen | Methods for non-linear least squares problems[END_REF]. Thanks to the global plane search combined with the Hessian computation, the proposed second order implementation of the RobustICA using the Newton method would provide a faster convergence.

Using the Lagrangian of F (3), denoted by L(w) we obtain:

L(w) = |E[y 4 ] -3| + β 2 (1 -w 2 ) (5)
First, the gradient of ( 5) is obtained as follows:

∇ w L(w) = 4 E[y 3 x] -βw (6)
where = sign(E[y 4 ] -3). The Hessian is then obtained as follows:

H(w) = 12 E[y 2 xx T ] -βI (7)
Subject to w = 1, the constant Lagrange multiplier β is calculated at the extrema where ∇ w L(w) = 0 as follows:

4 E[y 3 x] -βw = 0 ⇒ w T wβ = 4 w T E[y 3 x] ⇒ w 2 β = 4 w T E[xy 3 ] ⇒ β = 4 E[w T xy 3 ] ⇒ β = 4 E[y 4 ] (8)
By substituting (8) in ( 6) and ( 7), the gradient and the exact Hessian expressions of L (5) are respectively obtained as follows:

∇ w L(w) = 4 (E[y 3 x] -E[y 4 ]w) (9)
H(w) = 4 (3E[y 2 xx T ] -E[y 4 ]I) (10)
The following pseudo-code summarizes the Newton-based RobustICA approach, denoted by RobustICA N :

1) Initialization: random initialization of w 0 , k ← 0, Hessian initialization H 0 = I (where I is the identity matrix) and calculate ∇ w0 L. 2) Repeat until convergence or a maximum number of iterations is reached a) Compute the search direction g k = -H -1 k ∇ w k L where H k is the Hessian of L (5). b) Compute the optimal step size µ opt in the search direction g k as in [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF] c) Update w k+1 ← w k + µ opt g k . d) Normalize w k+1 ← w k+1 / w k+1 . e) Compute ∇ w k+1 L using (9). f) Compute H k+1 using (10). g) k ← k + 1. 3) Deflate the observations using the Gram Schmidt approach [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF] and go back to step (2). Despite its high convergence speed (few iterations), the Newton direction scheme requires an exact computation of the Hessian which is time consuming and computationally demanding per iteration. Therefore, quasi-Newton directions have been proposed wherein the Hessian is rather well approximated using the gradient of the considered objective function. The BGFS approach is one of the most popular quasi-Newton methods. It enjoys the Hessian's autocorrection property. That is to say, the approximated Hessian tends to correct itself in few iterations once an adequate line search defined by Wolfe conditions [START_REF] Nocedal | Numerical Optimization[END_REF] is adopted. It is worth noted that thanks to the global plane search used in the RobustICA, those conditions are intrinsically verified when implementing the BFGS-based RobustICA approach, denoted by RobustICA BFGS . The following pseudocode summarizes the proposed RobustICA BFGS approach:

1) Initialization: random initialization of w 0 . Set k ← 0,

H -1 0 = I and calculate ∇ ω0 F . 2) Repeat until convergence or a maximum number of iterations is reached a) Compute the search direction g k = -H -1 k ∇ w k F . b)
Compute the optimal step size µ opt in the search direction g k as in [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF]. c) Update w k+1 ← w k + µ opt g k . d) Normalize w k+1 ← w k+1 / w k+1 . e) Calculate ∇ w k+1 F as given by [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF].

f) Define s k = w k+1 -w k and y k = ∇ w k+1 F - ∇ ω k F . g) calculate ρ k = 1/(y k T s k ). if k = 0, H -1 k = y k T s k y k T y k I; end if. H -1 k+1 ← (I -ρ k s k y k T )H -1 k (I -ρ k y k s k T ) + ρ k s k s k T h) k ← k + 1.
3) Deflate the observations using the Gram Schmidt approach [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF] and go back to step (2).

III. EXPERIMENTAL ANALYSIS

This section is devoted to compare the performance of the proposed implementations of the RobustICA, namely, RobustICA CG and the other proposed Newtonbased schemes, the RobustICA N and RobustICA BGFS to the RobustICA method in its typical gradient-based implementation denoted by RobustICA G [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF] and the FastICA algorithm in its deflationary version [START_REF] Hyvärinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF]. The FastICA algorithm indeed stands for the most popular method for ICA, due to its simplicity, convergence speed, and satisfactory results in numerous applications. Experiments are carried out on synthetic in vivo MRS data acquired using an array of 32 sensors. Particularly, a realistic MRS of two metabolites (Creatine and Myo-inositol) is generated using Lorentzian and Gaussian functions [START_REF] Moussaoui | Non-negative source separation using the maximum likelihood approach[END_REF]. The performance of the considered implementations is evaluated in terms of the Normalized Mean Square Error (NMSE) [START_REF] Albera | ICA-based EEG demonising: a comparative analysis of fifteen methods[END_REF], the execution time, the global numerical complexity in floating point operations (flops) and the total number of iterations required to extract the P sources. Table I shows the numerical complexity per iteration of the considered methods. The aforementioned performance criteria are evaluated for different SNR values, i.e. SNR ∈ {0, 5, 10, 15, 20, 25, 30, 35, 40} dB. Results are averaged over 40 Monte Carlo realizations wherein the noise and the mixing matrix were generated randomly.

Figure 1 shows a mixture of the two above metabolites for SNR equal to 5 dB (Figure 1 (Top)) and the estimated sources using the considered approaches compared to the original sources (red line) (Figure 1 (bottom)). For the sake of clarity i) signals in this figure are slightly vertically shifted with respect to each other and ii) only the extracted signals using the RobustICA G , the RobustICA CG , the RobustICA BFGS are shown where similar behavior for the rest of the considered methods was observed. As shown in figure 1 (bottom), all the proposed RobustICA's variants show quasi similar behavior with comparable estimation accuracy. This claim is justified by the NMSE values shown in figure 2 for all SNR ones. Furthermore, the latter figure shows the benefit behind extending the gradient ascent scheme to the conjugate gradient one in terms of numerical complexity and execution time. Figure 2 shows also the efficiency of the proposed RobustICA BFGS . Indeed, it provides a competitive estimation quality, a lower execution time and numerical complexity (reduction around 50%) compared to the original RobustICA G . Generally the new RobustICA CG implementation and the proposed RobustICA BFGS provide a reduction of both the execution time and the numerical complexity over the conventional RobustICA G and the conventional FastICA algorithms. Note that, although the FastICA has the lowest numerical complexity per iteration, the entier numerical complexity of both the RobustICA CG and the RobustICA BFGS stands for the lowest as they converge in much less iterations. Regarding the RobustICA N , it shows the smallest number of iterations for convergence but the highest numerical complexity among the all considered methods. This high numerical complexity is due to the necessary exact Hessian computation step.

IV. CONCLUSION

In this paper, two Newton-based deflationary ICA algorithms have been investigated as a variant of the exact line search kurtosis-based algorithm, namely RobustICA. They are either based on an exact Hessian computation, RobustICA N or on an approximation of this latter using the BFGS scheme giving rise to the RobustICA BFGS approach. Both the second order optimization-based approaches were also compared to the deflation-based FastICA and the original RobustICA methods. Besides, the nonlinear CG-based implementation, the RobustICA CG , was considered in this study as a natural extension of the gradient ascent-based implementation initially adopted in the RobustICA method. All of the considered variants inherit the exact line search property from the original RobustICA algorithm. Experimental results on MRS data showed around 50% reduction in both of execution time and numerical complexity using the proposed quasi-Netwton (i.e. the BGFS method) and the CG-based schemes compared to the original gradient ascent procedure for the RobustICA method. This is while maintaining a good estimation accuracy. Regarding the Newton-based implementations, it showed, despite its fast convergence (i.e. very small number of iterations), a higher execution time and numerical complexity compared to the conventional RobustICA method and the other proposed implementations notably its quasi-Newton counterpart, the RobustICA BFGS . This behavior was due to the heavy numerical operations required to compute the exact Hessian.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Top: 2 observations of the Mixture of two metabolites (Creatine and Myo-inositol) for SNR = 5 dB. Bottom: Original sources (Red) vs. estimated sources using RobustICA G (Blue), RobustICA CG (Green) and RobustICA BFGS (Black).

TABLE I THE

 I NUMERICAL COMPLEXITY PER ITERATION IN FLOPS FOR THE FIVE CONSIDERED KURTOSIS-BASED METHODS. N IS THE NUMBER OF OBSERVATIONS, T IS THE SAMPLE SIZE AND Kp IS THE TOTAL NUMBER OF ITERATIONS NEEDED TO EXTRACT THE P-TH SOURCE.

		Numerical complexity per
		iteration [flops]
	RobustICA G	(5N + 12)T
	RobustICA CG	(5N + 12)T
	RobustICA BFGS	(5N + 12)T
	RobustICA N	( N 3 2Kp + N 2 2 + 5N + 12)T
	FastICA	(2N + 2)T


Emmanuel Lhuillier

Sandrine Ithurria

Armel Descamps-Mandine

Thierry Douillard

Remi Castaing

Zhen Xiang

Xu

Pierre-Louis Taberna

Patrice Simon

Herve Aubin

Benoit Dubertret

Investigating the n-and p-Type Electrolytic Charging of Colloidal Nanoplatelets

We investigate the ion gel gating of 2D colloidal nanoplatelets. We propose a simple, versatile, and air-operable strategy to build electrolyte-gated transistors. We provide evidence that the charges are injected in the quantum states of the nanocrystals. The gating is made possible by the presence of large voids into the NPL films and is sensitive to the availability of the nanocrystals surface.

INTRODUCTION

Nanoplatelets 1 (NPLs) are 2D colloidal materials, and the interest for this new type of colloidal nanoparticles has been mostly driven by their original optical properties. In particular, the atomic-scale control of the thickness [START_REF] Ithurria | Colloidal nanoplatelets with twodimensional electronic structure[END_REF][START_REF] Ithurria | Quasi 2D Colloidal CdSe Platelets with Thicknesses Controlled at the Atomic Level[END_REF] enables the synthesis of samples in which no inhomogeneous broadening is observed, leading to particularly narrow optical features. Thanks to this behavior, NPL appears a good candidate for optoelectronic devices with well controlled optical features. In addition, compared with 0D colloidal quantum dots (CQDs), NPL can also offer some potential improvement in charge transport. The 2D shape of NPL can lead to a larger overlap of the electron or hole wave function with neighboring NPL, while reducing the number of hopping steps to reach the electrodes thanks to their larger lateral size. The integration of NPL into transistors [START_REF] Lhuillier | Electrolyte Gated Field Effect Transistor to Probe the Surface Defects and Morphology in Films of Thick CdSe Colloidal Nanoplatelets[END_REF] and phototransistors [START_REF] Lhuillier | Electrolyte gated colloidal nanoplatelets based phototransistor and its use for bicolor detection[END_REF] has recently been reported. The photoconductivity of NPL has even been shown to be enhanced by 8 decades when they are integrated with nanometer scale electrodes [START_REF] Lhuillier | Nanoplatelets bridging a nanotrench: a new architecture for photodetectors with increased sensitivity[END_REF] compared with the same material operated in the hopping regime. For all of these devices, electrolyte gating was achieved through the use of ion gel electrolytes. It has been observed that this allows a fine-tuning of the conductance and photoconductance of the system. The use of ion gel gating instead of liquid electrolyte gating [START_REF] Wang | Electrochromic nanocrystal quantum dots[END_REF][START_REF] Vanmaekelbergh | Electrochemical gating: a method to tune and monitor the (opto) electronic properties of functional materials[END_REF] was first motivated by the need for an all-solid system that allows a size reduction and thus easier integration in final devices. While the interest in electrolytic gating of NPL for optoelectronic purpose is high, the gating mechanism in films of 2D NPL still requires a deeper understanding.

Electrolyte gating was first introduced in liquid form for films of quantum dots by the Guyot-Sionnest 7,9 and Vanmaekelbergh groups. [START_REF] Vanmaekelbergh | Electrochemical gating: a method to tune and monitor the (opto) electronic properties of functional materials[END_REF][START_REF] Germeau | Optical Transitions in Artificial Few-Electron Atoms Strongly Confined inside ZnO Nanocrystals[END_REF] Compared with solid-state back gating, electrolyte gating achieves a higher gate capacitance. This results from the fact that in back gating the capacitance is driven by the ratio of the dielectric contant of the gate material over the thickness of the insulating layer, which is at least a few multiples of the lattice parameter (greater than a few nanometers in practice). On the contrary, in the case of electrolytes, this is the size of the ion (<0.1 nm), rather than the electrolyte thickness, which drives the capacitance value and leads to a higher capacitance. As a consequence, while back gating remains limited to sheet densities <10 13 cm -2 , capacities >10 14 cm -2 have been obtained using electrolytes. [START_REF] Kang | High carrier densities achieved at low voltages in ambipolar PbSe nanocrystal thinfilm transistor[END_REF] The latter consequently allows injection of several charges per nanoparticle, which is generally not possible via back gating. Another key advantage for electrolyte gating is the ability to charge thick films, thanks to the percolation of ions within the film of nanoparticles. More recently, polymer and ion gel gating has been adapted to gate films of CQD, [START_REF] Kang | Sizedependent electrical transport in CdSe nanocrystal thin films[END_REF][START_REF] Kang | Influence of Silver Doping on Electron Transport in Thin Films of PbSe Nanocrystals[END_REF][START_REF] Kang | Size-and temperature-dependent charge transport in PbSe nanocrystal thin films[END_REF][START_REF] Bisri | Low driving voltage and high mobility ambipolar field-effect transistors with PbS colloidal nanocrystals[END_REF] in particular, as a probe of their doping.

We investigate the ion gel gating of NPL films. We demonstrate that this method allows fabrication of air-operable transistors with n-or p-type doping, depending on the cadmium chalcogenide compound used. The charging mechanism is further characterized by spectro-electrochemistry and impedance spectroscopy measurements. Finally, we study the effects of the physical and chemical properties (surface chemistry, film morphology, and temperature) of the films on the charging mechanism.

EXPERIMENTAL SECTION

2.1. Material Synthesis. CdSe Nanoplatelets. In a first step, cadmium(II) myristate (Cd(Myr) 2 ) is prepared following ref 16. In a typical synthesis, 240 mg of Cd(Myr) 2 and 25 mg Se powder are mixed in 30 mL of octadecene (ODE). The mixture is degassed under vacuum for 20 min at room temperature. Then, the atmosphere is switched to Ar, and the temperature is set to 240 °C. At 204 °C, 160 mg of cadmium(II) acetate (Cd(OAc) 2 ) is quickly added. The mixture is held for 12 min at 240 °C. One mL of oleic acid is quickly injected to quench the reaction and the solution is cooled. The NPLs are precipitated by adding ethanol. After centrifugation, the obtained solid is redispersed in hexane. The cleaning procedure is repeated three times.

CdSe/CdS Core-Shell Nanoplatelets. The procedure follows the one described by Ithurria et al. in ref 19. NPLs in hexane are transferred to N-methylformamide (NMFA) by the addition of sodium sulfide in NMFA ([Na 2 S] = 2g•L -1 ). It is equivalent to the growth of a half layer. The CdSe-S 2-NPLs are then precipitated with acetonitrile to remove the excess sulfide and redispersed in NMFA. The second half layer is grown by the addition of Cd(OAc) 2 in NMFA. (The reaction is quasi-instantaneous.) Excess precursors are removed by precipitating the nanocrystals with a mixture of 5:1 toluene/ acetonitrile. To finish the CdSe-S 2--Cd(OAc) 2 , NPLs are resuspended in NMFA. It is possible to grow several layers using Na 2 S and Cd(OAc) 2 as precursors of sulfide and cadmium. Here the sample has 3.5 layers of CdS.

CdTe Nanoplatelets. The procedure follows the one described by Pedetti et al. in ref 17. Cadmium(II) propanoate (Cd(Prop) 2 ) is first prepared by mixing 1 g of CdO in 10 mL of propionic acid under Ar for 1 h. Then, the flask is then opened to air and the temperature is raised to 140 °C to concentrate the solution by a factor of 2. The resulting white mixture is precipitated by the addition of acetone. After centrifugation the solid is dried under vacuum for 24 h. In the glovebox, 1 M Te in tri-n-octylphosphine (TOPTe) is prepared by stirring 2.55 g of Te pellets in 20 mL of TOP for 4 days at room temperature. In a three-necked flask 0.13 g of Cd(Prop) 2 , 160 μL of oleic acid, and 10 mL of ODE are degassed for 90 min at 95 °C. Then, the atmosphere is switched to Ar and the temperature is raised to 210 °C. 0.2 mL of 1 M TOPTe is quickly injected into the flask. After 20 min, the reaction is quenched by adding 1 mL of oleic acid and cooling to room temperature. Ethanol is added to precipitate the CdTe NPL. The solid obtained after centrifugation is redispersed in hexane. This procedure is repeated three times.

Additional material characterizations and details of the synthesis of other materials are given in the Supporting Information.

2.2. Device Fabrication. Electrolyte Preparation. The electrolyte is a mixture of polyethylene glycol (PEG), with a given molar weight (from 6 × 10 3 to 88 × 10 3 g mol -1 ), and ions (typically LiClO 4 ). For a typical electrolyte, 50 mg of LiClO 4 and 230 mg of PEG (MW = 6000 g mol -1 ) are heated together at 150 °C on a hot plate in the glovebox for 2 h.

Electrode Fabrication. Electrodes are fabricated using a regular lift-off procedure on a Si/SiO 2 wafer with an oxide thickness of 500 nm. Thermally evaporated Cr (2.5 nm) and Au (17.5 nm) are used for the contact. Electrodes are composed of 25 interdigitated pairs, which are 2.5 mm long. Two sets of electrodes have been designed with a spacing of 10 and 20 μm. The W/L ratio is 6250 for the 10 μm spaced electrodes and 3125 for the 20 μm spaced electrodes.

Ligand Exchange Procedure. The ligand exchange is done directly on the nanoparticles in solution. For a typical ligand exchange, 15 mg of Na 2 S is dissolved in 0.5 mL of NMFA by sonication. Then, 1 mL of the solution of nanoparticles in a nonpolar solvent is added. The mixture is then stirred until the nanoparticles have transferred to the polar phase. Then, 2 mL of hexane is added, and after sonication, the nonpolar phase is removed. This step is repeated a second time. Finally, the nanoparticles are precipitated by the addition of acetonitrile, and after centrifugation the solid is redispersed in NMFA. For narrow band gap materials (PbS, PbSe, HgTe), the nanoparticles are dispersed in a mixture of hexane and octane (9:1 volume ratio). The electrodes are warmed at 125 °C on a hot plate for 2 min and then cooled to room temperature. The solution is drop-cast onto the interdigitated electrodes. The typical thickness for a film is 30 nm. Then, the film is dipped in a solution of Na 2 S in ethanol (1-5 g L -1 ) for 1 min for PbS and PbSe and 30s for HgTe. The film is then rinsed in pure ethanol. Finally, the film is annealed at 80 °C for 1 min.

Field Effect Transistor Preparation. To prepare the electrolyte transistor, we softened the electrolyte by warming it on a hot plate at 80 °C. The electrolyte is brushed on the film of nanoparticles processed with the S 2-ligand exchange. Typical thickness is ∼1 mm. The device is finally cooled down. Finally, a TEM grid is deposited on the electrolyte as a gate electrode. The sample can be brought up to the measurement setup in air. Figure 1 summarizes the different steps of the preparation of the electrolytic transistor.

2.3. Characterization. XRD measurements are made using a Philips X'pert X-ray diffractometer operating at 40 kV and 40 mA with an X-ray wavelength of 1.54 Å. UV-visible spectra are acquired using a Cary 5E spectrometer. Thin films are prepared for the SEM and TEM imaging in Figure 3 using a Ga 2+ Focus Ion Beam (NVision 40 from ZEISS). SEM images were made using the electron beam in the FIB instrument. Transmission electron microscopy images are acquired on a JEOL 2010F TEM.

Electrical measurements are performed either using a probe station operating under air or in an Oxford CFV cryostat. Bias and current measurements are made using two Keithley 2400 source meters. For the transfer curve (I DS vs V GS ) the bias step is 1 mV and the sweep rate is 2 or 3 mV s -1 . Cyclic voltammetry measurements are made with a film of nanoparticles with a known mass on a ∼1 cm 2 Cu substrate. A liquid electrolyte is used made of a mixture of 1 M LiPF 6 in a mixture of 1:1 ethylene carbonate and dimethyl carbonate. A Li counter electrode is used and the cell is connected to a multipotentiostat (VMP3 from Bio-Logic) and the bias sweep rate is 100 μV•s -1 . For electrochemical impedance spectroscopy data, the same VMP3 setup is used. In this case a 2 V DC bias was applied and modulated by a 1 mV AC signal.

RESULTS AND DISCUSSION

3.1. Ion-Gel-Gated Nanoplatelet Films. To build an iongel-gated transistor, we first start from the NPL solution (CdSe, CdSe/CdS or CdTe) and perform a ligand exchange step [START_REF] Nag | Metal-f ree inorganic ligands for colloidal nanocrystals: S 2[END_REF] to S 2-. Early motivation for this ligand exchange was to boost the carrier mobility, but it will also impact the carrier density, as discussed later in the text. The obtained solution of NPL is then deposited on Cr/Au electrodes. Finally, a gel electrolyte is brushed onto the NPL film and a Cu metallic grid is used as top gate electrode. The procedure is described in Figure 1.

In the case of CdSe/CdS NPL, we observe a strong n-type gating. The on-off ratio can be as high as 10 8 while operating the transistor under low drain (0.5 V) and gate bias (<2 V); see Figure 2a. This current modulation is one of the highest reported for colloidal nanocrystal films and is 1 order of magnitude higher than that for pure CdSe NPL. [START_REF] Lhuillier | Electrolyte Gated Field Effect Transistor to Probe the Surface Defects and Morphology in Films of Thick CdSe Colloidal Nanoplatelets[END_REF] In the off state, the system is limited by the leakage current through the electrolyte. The sub threshold slope can be as low a 78 mV/ decade, which is close to the thermodynamic limit (60 mV/ decade at room temperature). The difference is likely a result of filling of defect states introduced by the C-ALD procedure [START_REF] Boehme | Density of Trap States and Auger-mediated Electron Trapping in CdTe Quantum-Dot Solids[END_REF] used to grow the CdS shell.

The gating of CdTe nanocrystals has been far less investigated and was only recently observed. [START_REF] Boehme | Density of Trap States and Auger-mediated Electron Trapping in CdTe Quantum-Dot Solids[END_REF] CdTe NPL shows p-type charging; see Figure 2b. The performance is reduced compared with the CdSe/CdS film; nevertheless, the on-off ratio of 10 3 is the highest value reported for CdTe colloidal system, to the best of our knowledge.

3.2. Electrochemical Properties. One of the key advantages of our method to build the electrolytic transistor is their air operability. All of the measurements have been conducted in air, and cycling over a week (>150 cycles) barely affects the transistor performances. The thick ion gel layer probably protects the NPL film. This is a significant improvement compared with liquid electrolytic gating, which required both air-free preparation and use of the device.

The method is not limited to NPL and can be operated with a large range of nanocrystals: CdSe QD, CdS nanorods, PbS QD and cubes, PbSe QD, and HgTe QD; see Figure S10. In the following, we choose to focus on a film of CdSe/CdS NPL. We want to address how the 2D aspect of the NPL combined with the ion gel nature of the electrolyte affects the charging mechanism.

We start by performing spectroelectrochemistry on the NPL film. We follow, as a function of time and charging, the optical absorption of the NPL; see Figure 3a. The CdSe/CdS NPL absorption spectra typically present two peaks associated with the heavy-hole to conduction band (λ ≈ 630 nm) and the lighthole to conduction band transition. As a function of charging, we observe a bleach of the optical features resulting from the filling of the conduction band by the injected electron. We can conclude that the nanoparticle quantum states are filled by the electrolyte gating.

The viscosity of the electrolyte combined with the 2D shape of the NPL may slow the percolation process compared with liquid electrolytes. We observe that in response to a step of gate bias the gate current typically has two time constants: one around 1 s and one around 1 min. To further characterize the dynamics of the charging, we use electrochemical impedance spectroscopy. We typically apply a 1 mV AC bias over a 2 V DC gate bias, which is above the turn on voltage of the electrolyte transistor, and measure the current, allowing the evaluation of the associated impedance and phase. At high frequency (f between 50 Hz and 1 kHz), the phase is close to 90°, which is the expected value for a pure capacitor; see Figure 3c. This behavior corresponds to the flattened half circle in the impedance plot; see Figure 3b. We can relate this behavior to the formation of the electrostatic double layer. [START_REF] Yuan | Electrostatic and Electrochemical Nature of Liquid-Gated Electric-Double-Layer Transistors Based on Oxide Semiconductors[END_REF] By further slowing down the dynamics, the phase drops and gets close to 0°, as expected for resistive behavior. In the impedance plot, the low frequency part of the diagram corresponds to the tail at high Re(Z). The resistive behavior is usually associated with diffusion-limited process. To model the electrochemical properties of the NPL film coupled to the ion gel electrolyte, we use the equivalent circuit [START_REF] Lu | In Supercapacitors: Materials, Systems and Applications[END_REF] proposed in Figure 3b. The circuit is composed of an electrolyte resistor (R el ≈ 500Ω) in series with a RC circuit made of a capacitor C DL modeling the double-layer gating and R CT which accounts for the chargetransfer resistance. To take into account that the half circle of the impedance plot is flattened we model the impedance with the expression C DL = (1/(Y O (jw) a )) where j 2 = -1, w is the pulsation and, Y O is the magnitude. We obtain a good fit with α = 0.8, while α = 1 will have described an ideal RC circuit. Finally, we add a diffusion term W (Warburg diffusion impedance) to model the tail occurring at low frequency. In general, W is expressed as W = (A/√w) + (A/(j√w)), where A is an amplitude that in our case is 6 ± 2 kΩ•s -1/2 . This term leads to a 45°slope in the impedance curve and accounts for the diffusion of the ion into the electrolyte. However, we typically observe a much smaller slope, which comes as a deviation compared with the equivalent circuit modeling; see Figure 3b. We attribute this deviation to a second source of diffusion-limited process, which is the percolation of the Li + ions within the nanoparticle film. This diffusion is known to be sensitive to the presence of voids, which leave space for the ions to diffuse. [START_REF] Boehme | Electrochemical Charging of CdSe Quantum Dot Films: Dependence on Void Size and Counterion Proximity[END_REF] In a previous study we demonstrated the fact that the gating was quite sensitive to the packing of the NPL. [START_REF] Lhuillier | Electrolyte Gated Field Effect Transistor to Probe the Surface Defects and Morphology in Films of Thick CdSe Colloidal Nanoplatelets[END_REF] In particular, the NPL core tends to stack, which prevents the Li + percolation and makes the gating ineffective. We conducted a study of the CdSe/CdS NPL film morphology by extracting a layer of the film using focused ion beam processing. The film is then further characterized using SEM (Figure 4a) and TEM (Figure 4b,c). We observe that despite the cubic aspect of the NPL the film is full of voids with a broad size distribution between 1 and 50 nm. These voids are responsible for the gating, even for thick films (up to micrometers). The less effective gating obtained with the CdTe NPL compared with the CdSe/CdS may not result from the material but instead from the difference in lateral size. The CdTe NPLs typically have a lateral extension >100 nm, and this makes the percolation path for the cation more complex.

3.3. Structural Properties of the Film. A striking feature related to this gating process is its dependence on the surface chemistry. As previously stated, the ligand exchange was first motivated by the enhancement of the interparticle coupling; however, cyclic voltammetry (CV, see Figure S13) showed significant variations with different capping ligands; see Figure 4d.

The integrated charge extracted from the CV curve (see Figure 4d) is significantly higher for atomically short (S 2-) passivation by a factor 5 to 10. This latter property remains true even after several cycles. Over the first cycles, the reduction of undesired organic species in the solution leads to a large asymmetry between the oxidation and reduction current. The overall increase in the injected charge into the film results from the higher availability of nanocrystal surface, which leads as a response to an easier tunnel injection of the electrons into the quantum states by the electrodes. This is consistent with the use of similar nanoparticle surface chemistries in the field of batteries. [START_REF] Oszajca | Precisely Engineered Colloidal Nanoparticles and Nanocrystals for Li-Ion and Na-Ion Batteries: Model Systems or Practical Solutions?[END_REF] In the last section we discuss the temperature dependence of this electrolyte gating and its stability. Electrolyte gating generally may lead to some charge fluctuation even under constant potential. To limit this fluctuation, it is useful to freeze the electrolyte by cooling the sample below the melting point of the matrix. In the case of ion gel gating the higher viscosity of the matrix limits this fluctuation. While at room temperature (20 °C) we observe an efficient gating, the charge injection is completely stopped already at 0 °C. In the PEG matrix, the Li + diffuses under a solvated form where an ether crown is formed, and its rearrangement allows the ionic displacement. At 0 °C the PEG matrix freezes, as observed on the XRD diffractogram (see Figure 4e), while the latter presents no fine structure at room temperature. The formation of crystallines zones prevents the rearrangement of the ether crown and the ionic displacement is suppressed. Having a polymer electrolyte with a freezing temperature close to room temperature consequently reduces the charge fluctuation compared with electrolyte with a lower freezing point.

Figure 1 .

 1 Figure 1. (a) TEM images of CdSe/CdS nanoplatelets. (b) TEM images of CdTe NPL. (c) Schematic of the ligand exchange procedure. (d) Schematic of the final electrolytic transistor.

Figure 2 .

 2 Figure 2. Transfer curve (drain current as a function of gate bias) for CdSe/CdS NPL (a) and CdTe (b) film. The gate current is also plotted to estimate the leakage current.

Figure 3 .

 3 Figure 3. (a) Top: absorption spectrum of the CdSe/CdS film. Bottom: the change in the absorption at different times under a 2 V gate bias. (b) Nyquist plot (imaginary part as a function of the real part) of the impedance of a CdSe/CdS film. The data are fitted using the proposed equivalent circuit given in inset. (d) Phase plot associated with the impedance plot. For panels b and c, a 2 V DC gate bias is modulated by a 1 mV AC bias.

Figure 4 .

 4 Figure 4. (a,a′) SEM images of a slice of CdSe/CdS film cut using FIB. (b,c) TEM images of the same film as panel a. (d) Apecific charge for film of CdSe/CdS NPL capped with butylamine and S 2-extracted from the cyclic voltammetry. (e) XRD for as-prepared electrolyte and after 1 h at 0 °C.

CONCLUSIONS

We propose a versatile method for electrolyte gating of NPL and nanocrystal films. The obtained devices, almost solid thin solids, are air operable. We provided evidence that the charges are indeed injected in the quantum states of the nanocrystals. This efficient gating is made possible by the large porosity of the NPL film with large voids letting the ions percolate even through thick films of NPL.

■ ASSOCIATED CONTENT

* S Supporting Information

The Supporting Information is available free of charge on the ACS Publications Web site at DOI: The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.5b05296.

Additional data giving details of the material synthesis and complementary information related to the transistor characterization. (PDF)

■

Notes

The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS

B.D. thanks ANR (grant SNAP) for funding. We are grateful to the METSA network and the CLYM laboratory for performing imaging of NPL film. We thank Sean Keuleyan for the corrections of the manuscript.