The complete knowledge of a quantum state allows the prediction of the probability of all possible measurement outcomes, a crucial step in quantum mechanics. It can be provided by tomographic methods [START_REF] Vogel | Determination of quasiprobability distributions in terms of probability distributions for the rotated quadrature phase[END_REF] which have been applied to atomic [START_REF] Leibfried | Experimental determination of the motional quantum state of a trapped atom[END_REF][START_REF] Kurtsiefer | Measurement of the Wigner function of an ensemble of helium atoms[END_REF] , molecular [START_REF] Dunn | Experimental determination of the quantum-mechanical state of a molecular vibrational mode using fluorescence tomography[END_REF] , spin [START_REF] Shulman | Demonstration of entanglement of electrostatically coupled singlet-triplet qubits[END_REF][START_REF] Medford | Self-consistent measurement and state tomography of an exchange-only spin qubit[END_REF] and photonic states. For optical [START_REF] Smithey | Measurement of the Wigner distribution and the density matrix of a light mode using optical homodyne tomography: application to squeezed states and the vacuum[END_REF][START_REF] Breitenbach | Measurement of the quantum states of squeezed light[END_REF][START_REF] Lvovsky | Quantum state reconstruction of the single-photon Fock state[END_REF] or microwave [START_REF] Bertet | Direct measurement of the Wigner function of a one-photon Fock state in a cavity[END_REF][START_REF] Houck | Generating single microwave photons in a circuit[END_REF][START_REF] Hofheinz | Synthesizing arbitrary quantum states in a superconducting resonator[END_REF][START_REF] Eichler | Experimental tomographic state reconstruction of itinerant microwave photons[END_REF] photons, standard tomography is obtained by mixing the unknown state with a large-amplitude coherent photon field. However, for fermions such as electrons in condensed matter, this approach is not applicable because fermionic fields are limited to small amplitudes (at most one particle per state), and so far no determination of an electron wavefunction has been made. Recent proposals involving quantum conductors suggest that the wavefunction can be obtained by measuring the time-dependent current of electronic wave interferometers [START_REF] Haack | Glauber coherence of single-electron sources[END_REF] or the current noise of electronic Hanbury-Brown/Twiss interferometers [START_REF] Samuelsson | Quantum sate tomography with quantum shot noise[END_REF][START_REF] Grenier | Single-electron quantum tomography in quantum Hall edge channels[END_REF][START_REF] Ferraro | Wigner function approach to single electron coherence in quantum Hall edge channels[END_REF] . Here we show that such measurements are possible despite the extreme noise sensitivity required, and present the reconstructed wavefunction quasiprobability, or Wigner distribution function [START_REF] Ferraro | Wigner function approach to single electron coherence in quantum Hall edge channels[END_REF] , of single electrons injected into a ballistic conductor. Many identical electrons are prepared in well-controlled quantum states called levitons [START_REF] Dubois | Minimal-excitation states for electron quantum optics using levitons[END_REF] by repeatedly applying Lorentzian voltage pulses to a contact on the conductor [START_REF] Levitov | Electron counting statistics and coherent states of electric current[END_REF][START_REF] Lebedev | Generating spin-entangled electron pairs in normal conductors using voltage pulses[END_REF][START_REF] Keeling | Minimal excitation states of electrons in one-dimensional wires[END_REF] . After passing through an electron beam splitter, the levitons are mixed with a weak-amplitude fermionic field formed by a coherent superposition of electron-hole pairs generated by a small alternating current with a frequency that is a multiple of the voltage pulse frequency [START_REF] Grenier | Single-electron quantum tomography in quantum Hall edge channels[END_REF] . Antibunching of the electrons and holes with the levitons at the beam splitter changes the leviton partition statistics, and the noise variations provide the energy density matrix elements of the levitons. This demonstration of quantum tomography makes the developing field of electron quantum optics with ballistic conductors a new test-bed for quantum information with fermions [START_REF] Lebedev | Generating spin-entangled electron pairs in normal conductors using voltage pulses[END_REF][START_REF] Beenakker | Proposal for production and detection of entangled electron-hole pairs in a degenerate electron gas[END_REF][START_REF] Beenakker | Optimal spin-entangled electron-hole pair pump[END_REF][START_REF] Sherkunov | Optimal pumping of orbital entanglement with single-particle emitters[END_REF] . These results may find direct application in probing the entanglement of electron flying quantum bits [START_REF] Yamamoto | Electrical control of a solid-state flying qubit[END_REF] , electron decoherence [START_REF] Ferraro | Wigner function approach to single electron coherence in quantum Hall edge channels[END_REF] and electron interactions. They could also be applied to cold fermionic (or spin-1/2) atoms [START_REF] Brantut | Conduction of ultracold fermions through a mesoscopic channel[END_REF] .

A quantum state y contains all the information about a particle or a system. Disregarding spin for simplicity, a practical representation is given by the product of the wavefunction Q which contains the spatial information and an occupation state in a Fock space representation. Although the determination of the latter requires similar experimental resources for a fermion and for a single boson, the complete determination of the wavefunction by tomographic methods is fundamentally more demanding for fermions than for bosons. To emphasize this let us consider a single particle propagating in a single spatial mode (an electron in an effectively one-dimensional ballistic conductor or a photon in an optical medium). In the space and time domain representation, the electron quantum state at position x and time t is Q(t 2 x/v F )jf ae, where v F is the velocity, Q is the wavefunction and jf ae is the fermionic occupation state, j1ae or j0ae. Similarly, a photonic state is characterized by E 0 u(t 2 x/c)jbae, where E 0 is the electric field amplitude of a single photon, u is the spatial mode which plays the role of Q, jbae is, for example, a number state or a Glauber coherent state, and c is the speed of light. In quantum optics, the determination of the mode u can be done by mixing with a coherent field (local oscillator) with amplitude ffiffiffiffiffiffiffiffi

N LO p E 0 u t{x=c ð Þ,
where N LO is the mean photon number. Then a classical measurement of u can be made because the fundamental quantum measurement uncertainty , 1

. ffiffiffiffiffiffiffiffi N LO p vanishes for large ffiffiffiffiffiffiffiffi N LO p . Only the determination of the bosonic state jbae requires purely quantum tomography. For fermions, however, the state occupation being limited to 1, no classical amplitude level can be reached and a quantum tomography is necessary for both the wavefunction Q and the fermion state jf ae. In the present work, we prepare many electrons, each in the known fermionic state jf ae 5 j1ae. We show that, although extremely demanding, a quantum tomography of the wavefunction of a time-resolved single electron is possible using shot noise, and we determine the first-order coherence QÃ e 0 ð ÞQ e ð Þ in the energy representation (where Q e ð Þ is the Fourier transform of Q(t) and the asterisk denotes complex conjugation).

We first explain how we prepare the single-electron state in a quantum conductor and then how the quantum tomography is performed. The quantum conductor is a quantum point contact (QPC) placed at low temperature (electron temperature, T e < 35 mK). It is formed in a high-mobility two-dimensional electron gas using split gates deposited on top of a GaAs/GaAlAs heterojunction (Fig. 1a and Supplementary Information). Applying negative voltage V G to the gates creates a constriction transmitting a limited number of electronic modes. Here we select a single mode whose transmission D(V G ) is obtained by measuring the conductance G~(2e 2 =h)D, where e is the electron charge and h is Planck's constant. The on-demand injection of charges into the conductor is obtained by applying a voltage pulse V L (t) on one of the ohmic contacts adjacent to the QPC, say the left. If the voltage flux satisfies e Ð z? {? V L t ð Þ dt~h then the charge is unity. Moreover, a Lorentzian pulse V t ð Þ~2Bw=e(t 2 zw 2 ), where B~h=2p and 2w is the pulse width at mid height, creates a clean single-electron state [START_REF] Levitov | Electron counting statistics and coherent states of electric current[END_REF] called a leviton [START_REF] Dubois | Minimal-excitation states for electron quantum optics using levitons[END_REF] with an energy just greater than the Fermi energy of the conductor. Remarkably, although the leviton originates from a collective displacement of the Fermi sea towards positive energy, it can be considered a single-particle state [START_REF] Lebedev | Generating spin-entangled electron pairs in normal conductors using voltage pulses[END_REF][START_REF] Keeling | Minimal excitation states of electrons in one-dimensional wires[END_REF] disentangled from the Fermi sea [START_REF] Keeling | Minimal excitation states of electrons in one-dimensional wires[END_REF] . A wavefunction representation of the leviton in the energy domain is Q e ð Þ~ffi ffiffiffiffiffiffiffiffiffi ffi 2w=B p h e ð Þe {ew=B , where h is the Heaviside step function. Experimentally, periodic trains of levitons have been generated [START_REF] Dubois | Minimal-excitation states for electron quantum optics using levitons[END_REF] . When partitioned by a QPC [START_REF] Dubois | Minimal-excitation states for electron quantum optics using levitons[END_REF] they generate a current shot noise corresponding to electrons prepared in the occupation state j1ae. Using shot-noise spectroscopy and electron Hong-Ou-Mandel noise correlation techniques, their wavefunction probability in the energy and time domain was found to be consistent with expectations. Thus, levitons provide a clean way of testing wavefunction quantum tomography in a quantum conductor.

We now discuss how to probe the first-order coherence from which the full wavefunction can be extracted. A first approach could be wave interferometry. In quantum optics, this gives G 1 (t 2 t9), which is the firstorder coherence AEy 1 (t9)y(t)ae averaged over the mean time t~tzt 0 ð Þ=2 (here y 1 is the Hermitian conjugate of y). The information is incomplete, however, and other approaches must be used [START_REF] Polycarpou | Adaptive detection of arbitrarily shaped ultrashort quantum light states[END_REF] . Wave interferometry is more promising with electrons: measuring the time-dependent current of an electronic Mach-Zehnder interferometer has been proposed [START_REF] Haack | Glauber coherence of single-electron sources[END_REF] as an elegant way to determine the full coherence AEy 1 (t9)y(t)ae, but so far no such experiments have been done. Another approach is performing shot-noise measurements. This is the approach followed in the present experimental work. Photon shot noise has been exploited using adaptive homodyne tomography [START_REF] Polycarpou | Adaptive detection of arbitrarily shaped ultrashort quantum light states[END_REF] to infer the unknown spectrotemporal properties of a single photon mode u. A large local oscillator field is mixed with the single photon in a beam splitter and its temporal shape is tuned to maximize the output noise, providing a classical copy of the unknown mode (Fig. 1b). For electrons, a large fermion field cannot be used but the principle of homodyne mixing can be kept. A recent proposal is to reverse the amplitude hierarchy [START_REF] Grenier | Single-electron quantum tomography in quantum Hall edge channels[END_REF] . Because it is the single-electron source corresponding to the maximum fermionic amplitude (as jf ae 5 j1ae), this is the local oscillator fermionic probe source whose amplitude has to be weak so that simple information can be provided by the noise variation at the beam-splitter output when mixing particles coming from the two sources. For simplicity, let us consider zero temperature. When a flux-quantized Lorentzian voltage pulse is applied to the left contact, a leviton is sent towards the QPC beam splitter. Repeating the experiment and recording the charge at the right contact gives the mean charge AEQae 5 De and the charge partition noise AEDQ 2 ae 5 D(1 2 D)e 2 . Now let us apply a d.c. voltage V R to the right contact, rising the electrochemical potential by eV R . The part of the leviton partition noise proportional to

Ð eVR 0 Q e ð Þ Ã Q e ð Þde, corresponding to the energy range e g [0, eV R ], is
replaced by the partition noise of the right electrons reduced by the factor 1{

Ð eVR 0 Q e ð Þ Ã Q e ð
Þde because of antibunching. Differentiating the noise with respect to V R provides the diagonal part of the energy density matrix of the leviton

Q e~eV R ð Þ j j 2 .
Remarkably, very little modification is necessary to access the non-diagonal part AEy 1 (e9 ? e)y(e)ae.

Let us superimpose on the right d.c. voltage the a.c. voltage V LO (t) 5 g LO (hn/e)cos(2pn(t 2 t)) with small amplitude g LO = 1, frequency n and time delay t referred to the leviton emission time. This generates at the QPC a weak fermionic field made of electrons in a superposition of states with their initial energies shifted by 6hn and with probability amplitudes g LO e i2pnt /2 and 2g LO e 2i2pnt /2, respectively, the probability of keeping the initial energy being 1 2 g LO 2 < 1 (Fig. 1c). Their mixing with the levitons at the beam splitter gives the antibunching probability amplitude

a g LO ,e ð Þ~Q Ã e ð Þ 1{g LO 2 À Á z 1 2 QÃ ezhn ð Þg LO e i2pnt { 1 2 QÃ e{hn ð Þg LO e {i2pnt
Differentiating the noise with respect to eV R , and keeping only the a.c. contribution, gives a noise variation proportional to

a V G V G V(t) V LO (t) b Photon source LO Correlation Detector Detector c Leviton ϕ(ε + hν 0 ) ~VLO (t) V R 0 μ R Energy Energy ϕ(ε) ~ϕ(ε -hν 0 ) φ ~ƒR (ε) ε V(t)
g LO Q e ð ÞQ ezkhn 0 ð Þ {Q e{khn 0 ð Þ Q e ð Þ ð Þ ,
g LO Re½Q Ã eV R ð ÞQ eV R zhn ð Þ e {i2pnt {Q Ã eV R ð ÞQ eV R {hn ð Þ e i2pnt
Fully quantum tomography can then be done by varying the voltage frequency and time delay t, providing the real and imaginary parts of the first-order coherence [START_REF] Grenier | Single-electron quantum tomography in quantum Hall edge channels[END_REF] . Because the leviton wavefunction in the energy domain is real, we get

d DQ 2 h i dV R !g LO cos 2pnt ð Þ½Q Ã eV R ð ÞQ eV R zhn ð Þ {Q Ã eV R ð ÞQ eV R {hn ð Þ ð1Þ
Because only positive energy is concerned, equation (1) gives

Q eV R ð ÞQ eV R zhn ð Þfor 0 # eV R , hn, Q eV R ð ÞQ eV R zhn ð Þ {Q eV R {hn ð Þ Q eV R ð Þ for hn # eV R ,
2hn and so on, providing a full determination of the energy density matrix.

In our experiment, single-electron levitons are injected at a frequency n 0 by applying periodic Lorentzian voltage pulses from the left contact. The full-width at mid-height is 2w 5 30 ps and the period is T 5 n 0 21 5 166 ps. To extract the diagonal and off-diagonal parts of the energy density matrix, we follow the experimental approach discussed above with current noise measured instead of charge fluctuations. With V R 5 0 and in the absence of a.c. voltage, the partitioning of levitons generates the low-frequency current noise spectral density S I 0 (0) 5 4n 0 AEDQ 2 ae 5 4n 0 e 2 D(1 2 D) where double counting of the mode from spin degeneracy is included. Because of the periodicity,

y z e 0 ð Þy e ð Þ Pk d e 0 {e{khn 0 ð Þ QÃ ezkhn 0 ð Þ Q e ð Þ
, where d is the Dirac delta function and k is an integer. Only energies separated by multiple of hn 0 contribute to off-diagonal terms, and the frequency of V LO (t) must be n 5 6kn 0 . Because high harmonics give too small a signal, measurements are restricted to k 5 0, 1 and 2.

The diagonal part (k 5 0) is obtained by measuring the difference between the partition noise of single-electron levitons injected at frequency n 0 , using the periodic Lorentzian voltage V L (t), and that of continuously injected electrons provided by a constant voltage equal to the mean value of the voltage pulse V L (t)~hn 0 =e. The zero-temperature noise difference is [START_REF] Dubois | Minimal-excitation states for electron quantum optics using levitons[END_REF][START_REF] Dubois | Integer and fractional charge Lorentzian voltage pulses analyzed in the framework of photon-assisted shot noise[END_REF]

DS k~0 I m R ð Þ~S 0 I hn 0 ð m R 0 1{Q e ð Þ 2 À Á de 0 @ { ð ? m R Q e ð Þ 2 de{ m R {hn 0 j j 1 C A ð2Þ
On the right-hand side, the first term represents the noise of right electrons emitted at energy below m R 5 eV R , whose shot noise is reduced by antibunching with the leviton. The second term is the partition noise of levitons for energy above m R . It is followed by the subtracted d.c. shot noise. Because Q(e) varies by steps in hn 0 intervals, equation ( 2) displays linear variation by parts. Introducing the notation Q e ð Þ 2 ~rl,l for lhn 0 , e , (l 1 1)hn 0 , this yields 0 for m R , 0, (2(1 2 r 0,0 )m R 1 hn 0 2 m R )/ hn 0 for 0 , m R , hn 0 , 2(1 2 r 0,0 ) 2 2r 1,1 (m R 2 hn 0 )/hn 0 for hn 0 , m R , 2hn 0 , and so on. The finite electron temperature T e < 35 mK < 0.12hn 0 /k B smoothens the singularities separating linear variations. A finite-temperature expression, obtained in a Floquet scattering approach [START_REF] Pedersen | Scattering theory of photon-assisted electron transport[END_REF][START_REF] Dasenbrook | Floquet theory of electron waiting times in quantum-coherent conductors[END_REF] , is given in Supplementary Information. Figure 2a shows measurements of DS I k 5 0 versus m R for levitons (open blue circles). From the data we extract r l,l directly from a fit including the small-temperature effects. The solid red curve corresponds to the best four-parameter fit, from which the four first r l,l (that is, l 5 0 to 3) are extracted (Supplementary Information). Self-heating including a small temperature increase (35 mK for V R < 0 to 47 mK for V R < 3hn 0 /e) has been included. The r l,l are plotted in red in Fig. 2b. Comparison with theoretical values based on a Floquet scattering modelling (in black) is good within the experimental uncertainties.

For the measurement of the off-diagonal part (k ? 0), we measure the noise difference when switching on and off the amplitude of V LO (t) 5 g LO (hn/e)cos(2pn(t 2 t)) superimposed on the right d.c. voltage V R and synchronized with the Lorentzian pulses. The zerotemperature shot noise difference is

DS k I m R ð Þ~S 0 I hn 0 2kg LO cos 2pkn 0 t ð Þ | ð m R 0 Q e ð ÞQ ezkhn 0 ð Þ {Q e ð ÞQ e{khn 0 ð Þ ð Þ de ð3Þ
This expression also varies linearly by parts because Q e ð ÞQ ezkhn 0 ð Þ ~rl,lzk for lhn 0 , e , (l 1 1)hn 0 is a step function of the energy. For example, for k 5 1 the integral in equation ( 3) is 0 for m R , 0, r 0,1 m R / hn 0 for 0 , m R , hn 0 , r 0,1 2 r 1,2 (m R 2 hn 0 )/hn 0 for hn 0 , m R , 2hn 0 , and so on. Finite-temperature expressions are given in Supplementary Information.

Figure 3a, b shows measurements of DS I k versus t for k 5 1, g LO 5 0.095 and k 5 2, g LO 5 0.0425 at fixed m R < hn 0 . Clear oscillations are observed. The delay period in Fig. 3b is twice that observed in Fig. 3a. Before going further, we emphasize that observation of these oscillations demonstrates off-diagonal coherence in an electronic orbital state, here a leviton. No experiment has given direct evidence of off-diagonal coherence of electrons in a conductor yet. We now compare the offdiagonal terms extracted from the measurements with those predicted for a leviton. Figure 3c,d shows the experimental variation of DS I k with m R for k 5 1 and 2 (open blue circles), fixing t 5 0 to maximize the signal. In both plots, the solid black curves represent the theoretical signal expected for a leviton. The red dashed curves are the best fits, from which the off-diagonal terms r l,l1k can be extracted. For k 5 2 a four-parameter fit is used and l 5 0 to 3, and for k 5 1 a five-parameter fit is used because a finite, but very small, hole excitation contribution

RESEARCH LETTER

(l 5 21) was found necessary to improve the fit. They are shown in Fig. 3e, f. The errors bars are obtained from the standard error of the fit. They are much larger than those in Fig. 2b because the measured noise is about ten times smaller. The noise oscillations with phase in Fig. 3a, b and the noise variations with voltage in Fig. 3c,d are well reproduced by a model including only the 35-50 mK thermal rounding. Supplementary Fig. 2 provides additional tomographic noise measurements performed on two-electron sine-wave pulses. Here again good agreement is found between data and theory. Deviation arising from decoherence effects due to electronelectron interaction is not expected here because the two-dimensional Landau quasi-particle lifetime is longer than the (50 mK) thermal time for the typical energies probed here (Supplementary Information).

From the experimental values of Q e ð ÞQ ezkhn 0 ð Þobtained above we can now reconstruct part of the zero-temperature Wigner distribution function [START_REF] Ferraro | Wigner function approach to single electron coherence in quantum Hall edge channels[END_REF] (WDF) of the levitons: W t,e ð Þ~ð z?

{?

y z ezd=2 ð Þ y e{ ð d=2Þie {id t=B dd. The WDF gives the full information on the complex wavefunction. It has been calculated for a single Lorentzian pulse [START_REF] Haack | Glauber coherence of single-electron sources[END_REF] and for periodic pulses [START_REF] Ferraro | Wigner function approach to single electron coherence in quantum Hall edge channels[END_REF] . In the latter case, the WDF is also periodic. Its expression in terms of the quantities r l,l1k is given in Supplementary Information. Figure 4a and Fig. 4b respectively give three-dimensional plots of the experimental and theoretical WDFs. Because the experimental WDF is limited to the first two harmonics and the energy range 2hn 0 # e , 3.5hn 0 , the theoretical WDF is truncated accordingly for comparison. To better show the agreement between experiment and theory, we provide in Fig. 4c a cut of the WDF at energy e 5 1.25hn 0 (red, experiment; black, theory). Using the experimental errors of r l,l1k (l 5 0) in Fig. 3e, f, we have generated the two red dashed curves in Fig. 4c between which the experimental WDF is expected to be located. For reduced times t/T around ,0.3 and 0.7 we clearly find robust regions of negative WDF, a hallmark of a quantum state. Finally, as a particular application of the WDF, the integration over energy e gives the probability jQ(t)j 2 . Its experimental (red solid curve) and theoretical (black dashed curve) variations with time are given in Fig. 4d. The excellent agreement validates our experimental determination of the WDF from the tomographic noise procedure of ref. 16. This provides better information on the temporal shape of the wavefunction than the convolution jAEQ(t)jQ(t 1 t)aej 2 given by the electronic Hong-Ou-Mandel method used for levitons in ref. 18 or for single electrons in ref. 31.

Using the known leviton state, we have demonstrated that the quantum wave tomography of itinerant electrons can be done experimentally using shot noise. Probing the off-diagonal energy density matrix opens the possibility of studying more complex electron quantum states such as entangled electron quantum bits in ballistic conductors or complex electron quantum states built by interaction. Examples of the latter include Kondo electronic states, fractional quantum Hall quasi-particles or electron Andreev pairs built by superconducting proximity effects.

LETTER doi: 10

 10 .1038/nature13821 Quantum tomography of an electron T. Jullien 1 *, P. Roulleau 1 *, B. Roche 1 , A. Cavanna 2 , Y. Jin 2 & D. C. Glattli 1

Figure 1 |

 1 Figure 1 | Schematics of quantum wave tomography. a, Experiment principle. Periodic Lorentzian voltage pulses V(t) applied to the left contact of a two-dimensional electron gas inject unit charge pulses called levitons which are partitioned by an electronic beam splitter called a quantum point contact (QPC). The split-gate voltage V G controls the transmission D of the onedimensional electronic mode formed at the QPC. A d.c. voltage V R (not shown) and a weak a.c. voltage V LO (t) 5 (kg LO hn 0 /e)cos(2pkn 0 (t 2 t)) are applied to the right contact The latter generates a small flux of electrons and holes (red-blue wavy line) which interferes with the incoming levitons (blue bumps) and is analogous to the local oscillator used in quantum optics. The red and blue arrows indicate the directions of electron-hole and leviton excitation propagation, respectively. Bicolour symbols indicate the outgoing states made of levitons mixed with the small flux of electron and hole excitations. Measuring the low-frequency current shot noise while varying V R for different harmonics (k 5 0, 1 and 2) provides a tomographic measurement of the leviton energy density matrix. b, Optical analogue. The local oscillator (LO) emits a monochromatic photon field, here of high amplitude, which interferes in the beam splitter with photons emitted by the source to be analysed. The photons detected provide a tomographic measurement of the photonic state. c, Experiment principle in energy representation. All quantities are step functions of the energy with energy scale hn 0 . Left: the energy dependence of the wavefunction of periodic levitons emitted from the left contact. Right: the energy distribution of the right reservoir when the voltage V R 1 V LO (t) is applied. For small g LO , an electron emitted with energy e arrives at the QPC in a superposition of states of energy e 6 khn 0 . Its interference with the levitons changes the leviton partition noise by a quantity proportional to g LO Q e ð ÞQ ezkhn 0ð Þ {Q e{khn 0 ð Þ Q e ð Þ ðÞ , giving information on the

  giving information on the non-diagonal energy density matrix. f R 5 | Q(e) | 2 . d, The 6 GHz periodic Lorentzian voltage pulses are obtained from a frequency synthesizer by combining four harmonics with appropriate phases and amplitudes and sent to the left ohmic contact of the sample via 40 GHz d.c. transmission lines. To generate V LO harmonics 1 and 2 are derived from the synthesizer, attenuated and sent to the right contact. A time delay is added using a phase shifter. Bias tees separate the high-frequency components from the d.c. voltage bias and the detected low-frequency current noise. The latter is converted into voltage fluctuations. A computer performs fast Fourier transform cross-correlations after cryogenic amplification and fast acquisition. Appropriate d.c. voltage biases are applied to the left and right contacts.

Figure 2 |Figure 3 |Figure 4 |

 234 Figure 2 | Measurement of the diagonal part of the energy density matrix. a, The data (blue circles) display the shot noise measured when applying a d.c. voltage hn 0 /e, minus the shot noise measured when sending single-charge levitons at frequency n 0 . The data are plotted versus the d.c. voltage V R applied on the opposite contact. The red curve represents the best fit to equation (2), including finite temperature, from which the parameters r l,l are extracted. The temperature T e is ,35 mK at V R 5 0 and ,48 mK at V R 5 3hn 0 /e. The transmission is D 5 0.19 and S I 0 5 9.48 3 10 229 A 2 Hz 21 . Error bars, s.e.m. calculated from the cross-correlation noise spectrum standard deviation in a 800 kHz bandwidth and after n < 100,000 acquisitions. b, Diagonal part of the energy density matrix versus energy using the r l,l extracted from the fit. Comparison with theory for levitons is in black. Error bars, s.e. calculated from the best fit.

Acknowledgements We acknowledge the ERC Advanced Grant 228273 MeQuaNo and thank P. Jacques for technical help, P. Pari, P. Forget and M. de Combarieu for cryogenic support, and P. Degiovanni and C. Grenier for discussions improving the manuscript.

Supplementary Information is available in the online version of the paper.

Author Contributions D.C.G. designed the project. T.J. and P.R. made the measurements and did the data analysis. B.R. contributed to the data analysis. P.R., T.J., B.R. and D.C.G. wrote the article. The sample was provided by Y.J. on wafer from A.C.

Author Information Reprints and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interests. Readers are welcome to comment on the online version of the paper. Correspondence and requests for materials should be addressed to D.C.G. (christian.glattli@cea.fr).

LETTER RESEARCH


Ahmad Abdelfattah

Marc Baboulin

Veselin Dobrev

Jack J Dongarra

Christopher Earl

Joël Falcou

Azzam Haidar

Ian Karlin

Tzanio V Kolev

Ian Masliah

S Tomov

High-performance Tensor Contractions for GPUs

Keywords: Tensor contractions, Tensor HPC, GPU, Batched linear algebra, FEM, Applications

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The development of high-performance tensor algebra is important due to tensors' frequent use in physics and engineering, where tensors provide a foundational mathematical tool for brief, yet comprehensive, formulations and solutions of problems in areas such as elasticity, fluid mechanics, multi-physics, quantum chemistry, general relativity, and many others [START_REF] Kolda | Tensor decompositions and applications[END_REF]. Advances in microprocessors and storage technologies have made it feasible to target higher dimension and accuracy computational approaches that model mutilinear relations, e.g., in recent areas of high interest such as various data analysis applications and machine learning; that can also be formulated through tensors. At the same time, to enable these applications to efficiently use tensor computations on current hardware, and in particular GPUs, a number of research challenges must be addressed, including advances in the development of scalable, tensor-based algorithms, autotuning, and code generation techniques, that are targets of this paper, towards setting the foundations for a high-performance tensor algebra library for accelerators [START_REF] Baboulin | Towards a high-performance tensor algebra package for accelerators[END_REF].

Tensors are multi-dimensional arrays that can be used to describe physical properties featuring multilinear relations. Well known mathematical objects like scalars, vectors, and matrices can be generalized to tensors that are of order zero, one, and two, respectively. Also, tensor transformations like flattening of a tensor to matrices or reshaping of matrices into tensors, can be used to link tensor computations to the developments in high-performance numerical linear algebra (LA). Therefore, similar to many applications, tensor computations can also significantly benefit from representing their computations in terms of BLAS, as well as from other dense LA algorithms and techniques for multicore and GPU architectures. While a comprehensive LAPACK-style initiative to tensors may be still far away [START_REF]Future directions in tensor-based computation and modeling[END_REF], this work concentrates on the development of tensor contractions -a building block for tensor computations -through leveraging the current LA developments for GPU and multicore architectures in libraries like BLAS and MAGMA [START_REF] Tomov | Towards dense linear algebra for hybrid GPU accelerated manycore systems[END_REF], and more specifically the MAGMA Batched computational framework [START_REF] Dong | LU Factorization of Small Matrices: Accelerating Batched DGETRF on the GPU[END_REF][START_REF] Haidar | Batched matrix computations on hardware accelerators based on GPUs[END_REF][START_REF] Haidar | A Framework for Batched and GPU-Resident Factorization Algorithms Applied to Block Householder Transformations[END_REF].

Microprocessor and storage technology advances have significantly influenced the design of high-performance numerical algorithms and libraries over the years. While humans perceive well algorithms expressed in terms of scalar computations (tensors of order zero), advances towards vector machines in the 70's lead to the development of the LINPACK library to use vector operations (or 1 st -order tensors), which was redesigned for performance into LAPACK in the 80's to better use cache-based machines through matrix-matrix operations (2 nd -order tensors). Operations on higher-dimensional data were added in the 90's for distributed-memory systems, where ScaLAPACK was designed for 2D-block cyclic matrix distributions. For the emerging in the 00's multicore architectures, the PLASMA library introduced tiled algorithms and tiled data layouts (4 th -order tensors; see Figure 1). In the 2010's, the MAGMA libraries were designed for heterogeneous architectures, including current investigations on new data layouts, functionalities, batched computations, and possibly generalizations to tensors, in order to provide applications new functionalities to deal efficiently with multi-dimensional data. Figure 1 Left illustrates the notion of tensor and tensor contraction in DLA, as well as one of our tensor contraction design using a DSEL . Figure 1 Right illustrates the need of tensor contractions in machine learning. The computational characteristics in this case are common to many applications: the operations of interest are in general small, they must be batched for efficiency, and various transformations may have to be explored to transform the batched small computations to regular and therefore efficient to implement operations, e.g., GEMM.

Filters F F n Output O n n,k O n,k O = k,i D i ∑ n,i
Plans to modernize numerical libraries for up-coming supercomputers must consider the projected changes in system architectures that will define the landscape on which HPC software will execute. For example, one of the next generation of supercomputers, exceeding 150 petaflop peak performance, is aimed to be deployed by 2018 and to include hybrid nodes, coupling powerful latency-optimized processors (IBM Power9) with highly parallel throughputoptimized accelerators (NVIDIA Volta GPUs) through NVIDIA NVLink interconnect. The platforms will also benefit from 3D-stacked memories. The challenges are severe to prepare libraries for systems like this. The efforts must be multidisciplinary, incorporating LA, languages, code generations and optimizations, domain science and application-specific numerical algorithms. Of particular interest is BLAST, a high-order FEM hydrodynamics research code currently developed at LLNL. Compared to low-order methods, BLAST improves the accuracy of Arbitrary Lagrangian-Eulerian (ALE) hydrodynamic simulations and provides a viable path to extreme parallel computing and exascale architectures. A core requirement to enable this path however is the development of advanced methods and software for tensor contractions, as introduced, on GPU and multicore architectures.

Tensors in numerical libraries

The use of tensors in physics and engineering applications has motivated an extensive research in both algorithms for tensor algebra and computational aspects for tensor-based simulations. See the survey [START_REF] Kolda | Tensor decompositions and applications[END_REF], and the references cited there, for an overview of the linear algebra aspects of tensors research, including higher-order tensor decompositions, their applications, and available software. A general overview with current and future direction in tensor research is presented in the Future directions in tensor-based computation and modeling workshop report [START_REF]Future directions in tensor-based computation and modeling[END_REF].

There has been a number of software packages developed for tensors. Some are designed to be used through numerical computing mathematical environments like the Tensor Toolbox1 for Matlab, GRTensor II2 for Maple, or Ricci3 for Mathematica, and therefore do not target accelerators, and high-performance computing in general. A few are specialized for particular applications, most notably for quantum chemical computations. For example, [START_REF] Ma | Acceleration of Streamed Tensor Contraction Expressions on GPGPU-Based Clusters[END_REF] presents early work on using accelerators to accelerate tensor contractions on GPUs. The approach uses code generation techniques and is incorporated in the NW Chem computational chemistry suite. More recent work [START_REF] Nelson | Generating efficient tensor contractions for gpus[END_REF] also uses code generation techniques and autotuning (with a system called Baracuda, based on CUDA-CHiLL and a high-level module Optimizing Compiler with Tensor OPeration Intelligence (OCTOPI)) to report significant acceleration compared to NW Chem on particular tensor contractions.

While we use code generation and autotuning, our approach also relies on context knowledge, and in particular that tensor reshaping techniques can often transform tensor contractions from many applications into many GEMMs, similar to the example on Figure 1, Right from machine learning. This transformation can not easily be detect automatically, and moreover, even if the transformation is somehow indicated, it is well known that general GEMM optimizations using only compiler techniques can not match in performance best-tailored solutions. Similar to our approach but for CPU systems only, [START_REF] Solomonik | A massively parallel tensor contraction framework for coupled-cluster computations[END_REF] executes contractions via GEMM on a properly ordered and structured tensor. As in the other approaches in quantum chemistry, large tensor contractions are targeted. In contrast, we target many but small contractions, that are often sub-warp/vector in size (see next section), resulting in large-scale computations. Tensor reshuffle operations are done to cast the contractions to GEMMs, when possible, and a batched approach with custom-built BLAS kernels is used to efficiently execute them.

Additional closely related efforts include [START_REF] Markidis | Openacc acceleration of the nek5000 spectral element code[END_REF][START_REF] Shin | Speeding up nek5000 with autotuning and specialization[END_REF] and [START_REF] Stock | Model-driven simd code generation for a multi-resolution tensor kernel[END_REF].

3 Tensor formulation for high-order FEM

We derive tensor formulations for the high-order FE kernels of the Lagrangian phase of BLAST, which solves the Euler equations of compressible hydrodynamics in a moving Lagrangian frame [START_REF] Dobrev | High-order curvilinear finite element methods for lagrangian hydrodynamics[END_REF][START_REF] Dong | A step towards energy efficient computing: Redesigning a hydrodynamic application on CPU-GPU[END_REF]. On a semi-discrete level, the conservation laws of Lagrangian hydrodynamics can be written as:

Momentum Conservation: M V dv dt = -F • 1, (1)
Energy Conservation:

de dt = M -1 E F T • v , (2)
Equation of Motion:

dx dt = v, where (3)
v, e, and x are the unknown velocity, specific internal energy, and grid position, respectively; M V and M E are independent of time velocity and energy mass matrices; and F is the generalized corner force matrix depending on (v, e, x) that needs to be evaluated at every time step.

To illustrate the tensor formulation of this problem, consider the finite element (FE) mass matrix for an element (zone) E with a weight ρ, as a 2-dimensional tensor: are the FE basis functions on the reference element, |J E | is the determinant of the element transformation, and {q k } nq k=1 and {α k } nq k=1 are the points and weights of the quadrature rule. Tensors can be introduced by taking the nq × nd matrix B, B ki = ϕ i (q k ), and the nq × nq

(M E ) ij = nq k=1 α k ρ(q k ) ϕ i (q k ) ϕ j (q k ) |J E (q k )| , i, j = 1, . . . ,
diagonal matrix D E , (D E ) kk = α k ρ(q k ) |J E (q k )|. Then, (M E ) ij = nq k=1 B ki (D E ) kk B kj , i.
M i1,i2,j1,j2 = k1,k2 (B 1d k1,i1 B 1d k1,j1 )(B 1d k2,i2 B 1d k2,j2 )D k1,k2 , (4)
where

B 1d is a dense O(p) × O(p) matrix and D is a dense O(p) × O(p) matrix. In 3D, M i1,i2,i3,j1,j2,j3 = k1,k2,k3 (B 1d k1,i1 B 1d k1,j1 )(B 1d k2,i2 B 1d k2,j2 )(B 1d k3,i3 B 1d k3,j3 )D k1,k2,k3
, where (5)

D is a dense O(p) × O(p) × O(p) tensor, and i = (i 1 , • • • , i d ), j = (j 1 , • • • , j d ), k = (k 1 , • • • , k d )
are the decompositions of the dof and quadrature point indices into the tensor components along logical coordinate axes. Note that if we store the tensors as column-wise 1D arrays, then

M nd1×nd2×nd1×nd2 i1,i2,j1,j2 = M nd×nd i,j
= M nd 2 i+nd j = M nd 2 i1+nd1i2+nd(j1+nd1j2) , i.e. we can reinterpret M as a nd × nd matrix, or a fourth order tensor of dimensions nd 1 × nd 2 × nd 1 × nd 2 , or a vector of dimension nd 2 , without changing the underlying storage.

The action of the operator, U = M V , can be written in the tensor product case as

U i1,i2 = k1,k2,j1,j2 B 1d k1,i1 B 1d k2,i2 D k1,k2 B 1d k1,j1 B 1d k2,j2 V j1,j2 , and (6)
U i1,i2,i3 = k1,k2,k3,j1,j2,j3 B 1d k1,i1 B 1d k2,i2 B 1d k3,i3 D k1,k2,k3 B 1d k1,j1 B 1d k2,j2 B 1d k3 ,j3 V j1,j2,j3 . (7)
Given B 1d , D, and V , tensors M and U must be computed based on the generalized contractions ( 4)- [START_REF] Dong | A step towards energy efficient computing: Redesigning a hydrodynamic application on CPU-GPU[END_REF]. The matrix sizes are relatively small, e.g., with p = 2..8 and d = 2 or 3.

Tensor contraction interfaces and code generation

To develop high-quality HPC software for tensor contractions, we impose the following three main requirements on our interface design:

Convenience of use: Interfaces of HPC libraries are extremely important for the libraries' adoption by the community. Interfaces must provide convenience and practicality of use, including ease of interoperability between libraries. Ideally, a tensor data type standard must be defined. The standard for a dense matrix is a pointer, sizes, leading dimension, and assumption for column-major data layout, e.g., as in BLAS and LAPACK. For tensors, motivated by reviewing interfaces in available libraries and the success of BLAS, we propose to represent a tensor by its dimensions and a data layout abstraction. The abstraction is to provide a choice of predefined layouts, and ways to switch it, or to easily add user-specified layouts, without changing the underlying algorithms. In general, a specific layout provides the formula of mapping the multi-dimensional tensor to the memory, e.g., a second order tensor can be stored as a column-major matrix A with leading dimension lda, in which case the abstraction maps A i,j to A[i + j * lda] (see Listing 1).

Readability: Numerical software must be understandable, which is needed for ease of maintenance, as well as code optimizations, and testing. While we can easily implement any interface, e.g., even expressing the interface and tensor APIs in a DSEL if needed (plus code generation techniques to translate the DSEL to a standard language; see the example in Einstein tensor notations on Figure 1 Left), we determined that it is better to provide implementations relying on a standard language and the code generation features provided within the language. The C++14 standard for example is expressive enough to allow us to implement readable and easy to use interfaces.

Performance: While we expect to obtain high performance mostly through algorithmic design and autotuning (see Section 5), we do not want to compromise on optimization opportunities like removing parameters checking and loop unrolling to eliminate jumps and loop count decrements. These optimizations are essential especially for the small computations that we target. Therefore, in our design we consider the use of compiler features related to code generation (e.g., templates, etc.), as further discussed below.

Related to performance, a lost optimization opportunity is if static checking and compile time information is not provided as part of the scientific libraries used. As an example, LA-PACK routines start by checking entry parameters dynamically, which results in extra time for small size computations. The type of algorithms that we intend to use, e.g., as the MAGMA Batched algorithms, also require specific tuning [START_REF] Haidar | Batched matrix computations on hardware accelerators based on GPUs[END_REF] as performance will greatly vary depending on numerous parameters. To avoid these performance drawbacks, and also benefit from optimization techniques like compiler loop unrolling for static dimension problems, we use features of the new C++14 standard. In particular, the constexpr specifier enables to evaluate the value of a function or variable at compile time. We use this feature with integral constants and template specialization to design our tensor dimensions layout: We can then freely extract the dimensions (Listing 1) and use them to specify our CPU and GPU kernels at compile time. Our tensor model is based on our layout, the data type of the tensor and its locality. The memory buffer is based on vector from the Standard Template Library for CPU and thrust [START_REF] Bell | Thrust: A 2 6. GPU Computing Gems Jade Edition[END_REF] for GPU. To generate a tensor, we need to pass a data type and locality as template parameter and the size to the constructor (Listing 2).

Transfers between CPU and GPU tensors can be expressed through the function copy (Listing 2). This function will use the stream 0 by default but a stream can be passed for asynchronous copy. We have designed two models for batched computing (Listing 3). The first one is based on allocating a single memory block for all tensors to improve data transfers and locality while the other is a group of same size tensors. The generalized tensor contractions ( 4)-( 7) can be summarized to a few kernels [START_REF] Baboulin | Towards a high-performance tensor algebra package for accelerators[END_REF]. One natural way to implement them is as a sequence of pairwise contractions, i.e. by evaluating the sums one at a time. While this is an efficient approach, especially when coupled with a batched approach as in MAGMA [START_REF] Haidar | Batched matrix computations on hardware accelerators based on GPUs[END_REF][START_REF] Haidar | A Framework for Batched and GPU-Resident Factorization Algorithms Applied to Block Householder Transformations[END_REF], there is enough complexity here that maybe one can come up with something better. Indeed, a number of the kernels needed can be reduced further to a tensor index reordering (generalized transpose) plus the dgemm A, B → A T B. This is because contraction by the first index, for example nd2nd3) . If the contraction is not by the first index, reducing it to A T B requires data reordering. However, we note that there is a way not to do this explicitly since the matrices are very small; we organize our algorithms (next) to have a phase of reading A and B to shared memory, followed by a computational phase. Thus, the reading can be from consecutive or not data, and in either case to benefit from data reuse (of the small A and B in shared memory) using the same computational phase.

C i1,i2,i3 = k1 A k1,i1 B k1,i2,i3 , can be written as Reshape(C) nd1×(nd2nd3) = A T Reshape(B) nq1×(
In summary, all of the ( 6)-( 7) operations can be implemented through reshaping and the kernels A i,j,k,l = s B i,s,j C k,s,l and A i,k,l,j = s B i,s,j C k,s,l . The matrix assembly contractions (4)-( 5) can also be reduced to a common kernel (plus reshaping): A k,i,l,j = s B s,i B s,j C k,s,l .

Performance model

To evaluate the efficiency of our algorithms we derive theoretical bounds for the maximum achievable performance P max = F/T min , where F is the flops needed and T min is the fastest time to solution. For simplicity, consider C = αAB + βC on square matrices of size n. We have F ≈ 2n 3 and T min = min T (T Read(A,B,C) + T Compute(C) + T W rite(C) ). Note that we have to read/write 4n 2 elements, or 32n 2 Bytes for double precision (DP) calculations. Thus, if the maximum achievable bandwidth is B (in Bytes/second), we take T min = 4n 2 /B in DP. Note that this time is theoretically achievable if the computation totally overlaps the data transfers and does not disrupt the maximum rate B of read/write to the GPU

P max = 2n 3 B 32n 2 = nB 16 in DP.
The achievable bandwidth can be obtained by benchmarks, e.g., using NVIDIA's bandwidthTest.

For example, on a K40 GPU with ECC on the peak is 180 GB/s, so in that case P max is 11.25 n GFlop/s (denoted by the dashed line on Figure 3). Thus, when n = 16 for example, we expect a theoretical maximum performance of 180 GFlop/s in DP.

Algorithms for tensor contractions through batched GEMMs

As described, we transform the tensor contractions to batched GEMM. To achieve HP on batched GEMM for small matrices that are sub-warp in size, we apply the following techniques:

• Using templates and constexpr specifiers we manage to avoid parameters checking and get compiler-unrolled loops in our kernels;

• We avoid passing arrays of pointers to the batched matrices, which are quite large, by passing them through formula/function that is part of the tensor's structure definition;

• The kernels are organized as follows: (1) Read A and B into shared memory; (2) Compute AB in registers; (3) Update C. Reading A, B, and C is through functions in the tensor's structure definition, which allows us to work with matrices that are not stored in the standard matrix format. Thus, we do not need explicit data reordering for some contractions, in order to efficiently benefit from the GEMM computation (step 2 above);

• We versions based on: how A and B are read; prefetching or not C by intermixing reading C with the computation in (2); number of matrices handled by a thread block, combined with prefetching variations for A and B; number of threads per thread block; and combinations of the above.

In all cases, to maximize data reuse, a single GEMM is done on a single thread block. Batching of the computation is done as in [START_REF] Haidar | A Framework for Batched and GPU-Resident Factorization Algorithms Applied to Block Householder Transformations[END_REF]. A GPU GEMMs technique, used for large matrices since the Fermi architecture [START_REF] Nath | An Improved Magma Gemm For Fermi Graphics Processing Units[END_REF], is to apply hierarchical communications/blocking on all memory levels, including a new register blocking. Register blocking can be added by making a single thread compute more than one element of a resulting matrix (in the same row or column, so that when an element from A or B is loaded into a register, it gets used more than once). Current results show that register blocking may not be needed in this case.

Autotuning

The complexity of tuning our algorithms is handled through autotuning [START_REF] Li | A note on auto-tuning GEMM for GPUs[END_REF][START_REF] Abdelfattah | Performance, Design, and Autotuning of Batched GEMM for GPUs[END_REF]. We note that although tuning is important, the algorithmic design (as in Section 5.2) is the more critical part for obtaining high-performance, as the overall success is limited by the quality of the algorithms used. With that in mind, there are generally two kinds of approaches for doing autotuningmodel-driven optimization and empirical optimization. We use a combination. The modeldriven part comprises compiler code generation and optimizations (as in Section 4). For the empirical optimization, a large number of parametrized code variants are generated and run on a given platform to discover the one that gives the best performance. The effectiveness of empirical optimization depends on the chosen parameters to optimize, and the search heuristic used. The algorithm designs and implementations using templates are very convenient in setting up our autotuning framework. The process is illustrated on Figure 2. Further details on the autotuning framework can be found in [START_REF] Abdelfattah | Performance, Design, and Autotuning of Batched GEMM for GPUs[END_REF], while details on the kernels, including their extension to multicore CPUs are available in our technical report [START_REF] Masliah | High-performance matrix-matrix multiplications of very small matrices[END_REF]. Our performance is within 90% of the theoretical maximum, as derived in Section 5.1, and denoted by a dashed line. Slight improvement may be possible through register blocking and tuning.

Experiments on tensor computations

Conclusions and future directions

We presented an approach based on tensors for high-order FEM and we described it as a multidisciplinary effort that includes a tensor formulation of the numerical algorithms of a domain science problem, the definition of tensor interfaces and code generation, and the design of algorithms design and tuning for performance and portability. We developed a proof-ofconcept software for it, showing that it has the potential to get optimal performance. We currently achieve within 90% of a theoretically derived peak on GPUs, using on-the-fly tensor reshaping to cast tensor contractions to small but many GEMMs, executed using a batched computing approach, custom-built GEMM kernels for small matrices, and autotuning. As future work, there are many options that need to be coded and tuned. The resulting HP tensor contractions package will be released as open source through the MAGMA library. It remains to integrate the developments in the BLAST code, and complete the autotuning in anticipation for easy portability to future supercomputers like Sierra, and new GPU architectures like the Volta with its 3D-stacked memory, expected to further favor tensor computations.

∑

  Matrix A in tiled data-layout as a 4 th -order tensor: // Declare a 4 th -order Tensor A on the GPU Tensor<64, 64, 9, 8, gpu_t> A;  // DSEL design using Einstein notation: repeated // index k means a summation/contraction.  // Range of the other indices is full/range as // given through the left assignment operand A(i, j, m:1..8, n:1..7) -= A(i, k,m,0) * A(k, j,0,n);  A rank-64 update as tensor contraction on index k (for i = 0..63 for j = 0..63 for m = 1..8 for n = 1..7):Tensor contractions in machine learning (Convolutional Neural Networks in computer vision)

Figure 1 :

 1 Figure 1: Left: Example of a 4 th -order tensor resulting from tile matrix layout used in dense LA, a tensor contraction, and a possible tensor contractions design using Einstein summation notation and a Domain Specific Embedded Language (or DSEL ) . Right: Illustration of tensor contractions needed and viable approaches to solve them in machine learning.

  e., M = B T DB (omitting the element index E) . Using FEs of order p with a quadrature rule of order p, we have nd = O(p d ) and nq = O(p d ), so B is a dense O(p d ) × O(p d ) matrix. If the FE basis and the quadrature rule have tensor product structure, then in 2D,

/

  / Template s p e c i a l i z a t i o n c o n s t e x p r a u t o l a y o u t = o f s i z e <5 ,3 >() ; // U s i n g I n t e g r a l c o n s t a n t c o n s t e x p r a u t o l a y o u t 1 = o f s i z e ( 5 c , 3 c ) ; // U s i n g dynamic d i m e n s i o n s c o n s t e x p r a u t o l a y o u t 2 = o f s i z e ( 5 , 3 ) ; // A c c e s s D i m e n s i o n s a t c o m p i l e t i m e c o n s t e x p r a u t o dim1 = l a y o u t ( 1 ) ; Listing 1: Dimensions for Tensors // C r e a t e a r a n k 2 t e n s o r o f s i z e 5 , 3 on GPU c o n s t e x p r t e n s o r <f l o a t , gp u > d t s ( o f s i z e <5 ,3 >() ) ; // C r e a t e a r a n k 2 t e n s o r o f s i z e 5 , 3 on CPU c o n s t e x p r t e n s o r <f l o a t > t s ( o f s i z e <5 ,3 >() ) ; // Use t h r u s t t o f i l l d t s w i t h 9 t h r u s t : : f i l l ( d t s . b e g i n ( ) , d t s . end ( ) , 9 ) ; // Copy d t s f r o m GPU t o t s on CPU c o p y ( d t s , t s ) ; Listing 2: Create Tensor and copy As seen in Listing 1, we propose 3 ways to specify dimensions using the function of size which returns a layout containing the static or dynamic dimension. Each operator inside the layout uses the constexpr keyword which enables us to return sizes at compile time if possible.

/

  / C r e a t e a b a t c h t h a t w i l l c o n t a i n 15 t e n s o r s o f s i z e 5 , 3 , 6 c o n s t e x p r a u t o b a t c h<f l o a t , gpu > b = m a k e b a t c h ( o f s i z e ( 5 c , 3 c , 6 c ) , 1 5 ) ; // A c c e s s i n g a t e n s o r f r o m t h e b a t c h r e t u r n s a v i e w on i t c o n s t e x p r a u t o v i e w b = b ( 0 ) ; // C r e a t e a g r o u p i n g o f t e n s o r s o f same s i z e t e n s o r s c o n s t e x p r a u t o group<f l o a t , gp u > g ( o f s i z e ( 5 c , 3 c ) ) ; // Add a t e n s o r t o t h e g r o u p c o n s t e x p r a u t o t e n s o r <f l o a t , gp u > d t s ( o f s i z e ( 5 c , 3 c ) ) ; g . p u s h b a c k ( d t s ) ; Listing 3: Batched tensors Once we have defined these functions we can call the kernel to compute a batched dgemm on tensors of dimension 2.

  c o n s t e x p r a u t o b a t c h<f l o a t , gpu > b = m a k e b a t c h ( o f s i z e ( 5 c , 3 c ) , 1 5 ) ; c o n s t e x p r a u t o b a t c h<f l o a t , gpu > b1 = m a k e b a t c h ( o f s i z e ( 5 c , 3 c ) , 1 5 ) ; // P r o d u c t o f two t e n s o r b a t c h e d o f d i m e n s i o n 2 f o r m a t r i x p r o d u c t u s i n g C++ o p e r a t o r c o n s t e x p r a u t o c = b * b1 ; // P r o d u c t u s i n g t h e b a t c h dgemm f u n c t i o n t h a t c a n be s p e c i a l i z e d d e p e n d i n g on p a r a m e t e r s c o n s t e x p r a u t o c = batch gemm ( b , b1 ) ; Listing 4: Tensor Operations

Figure 2 :

 2 Figure 2: Autotuning framework for high-performance tensor contractions in MAGMA.

Figure 3 Figure 3 :

 33 Figure3shows the DP performance of four of our autotuned kernels from Section 5.2. Each of the kernels gives best results for certain dimensions. Shown is also the performance on 16

1) Kernel variants: performance parameters are exposed through a templated kernel interface

	template< typename T, int DIM_X, int DIM_Y,
	int BLK_M, int BLK_N, int BLK_K,
	int DIM_XA, int DIM_YA, int DIM_XB, int DIM_YB,
	int THR_M, int THR_N, int CONJA, int CONJB >
	static __device__ void tensor_template_device_gemm_nn( int M, int N, int K, …
	2)

CPU interfaces that call the GPU kernels as a Batched computation

		3) Python scripts that generate the search space for the
		parameters DIM_X, DIM_Y …
		index, DIM_X, DIM_Y, …
		#define NN_V_0	4,	8, 8, 24, 8, 1, 4, 8, 4, 8
		#define NN_V_1	4,	8, 8, 32, 8, 1, 4, 8, 4, 8
		#define NN_V_2	4,	8, 8, 40, 8, 1, 4, 8, 4, 8
	template<typename T, int DIM_X, int DIM_Y, … >	…
	void tensor_template_batched_gemm_nn( int m, int n, int k, … ) { …	4)
	tensor_template_device_gemm _nn<T, DIM_X, DIM_Y, … ><<<dimGrid, dimBlock, 0, queue>>>(m, n, k,…);
	}

Scripts that run all versions in the search space, analyze the results, and return the best combination of parameters, which is stored in the library for subsequent use.

http://www.sandia.gov/∼tgkolda/TensorToolbox/

http://grtensor.phy.queensu.ca/

http://www.math.washington.edu/∼lee/Ricci/

Acknowledgments

This material is based upon work supported by the National Science Foundation under Grants No. CSR 1514286 and ACI-1339822, NVIDIA, the Department of Energy, and in part by the Russian Scientific Foundation, Agreement N14-11-00190. This work was further performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. LLNL release number LLNL-CONF-681073-DRAFT.