Nicolas Desmars
email: nicolas.desmars@ec-nantes.fr

Yves Perignon
email: yves.perignon@ec-nantes.fr

Guillaume Ducrozet
email: guillaume.ducrozet@ec-nantes.fr

Charles-Antoine Guérin
email: charles-antoine.guerin@univ-tln.fr

Stephan T Grilli
email: grilli@uri.edu

Pierre Ferrant
email: pierre.ferrant@ec-nantes.fr

Yves P Érignon

Charles-Antoine Gu

Phase-Resolved Reconstruction Algorithm and Deterministic Prediction of Nonlinear Ocean Waves From Spatio-Temporal Optical Measurements

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.

ABSTRACT

We investigate a nonlinear phase-resolved reconstruction algorithm and models for the deterministic prediction of ocean waves based on a large number of spatio-temporal optical measurements of surface elevations. We consider a single sensor (e.g., LIDAR, stereo-video, etc.) mounted on a fixed offshore structure and remotely measuring fields of free surface elevations. Assuming a uniform distribution of measurement points over the sensor aperture angles, the density of free surface observation points geometrically decreases with the distance from the sensor. Additionally, wave shadowing effects occur, which become more important at small viewing angles (i.e., grazing incidence on the surface). These effects result in observations of surface elevation that are sparsely distributed. Here, based on earlier work by [START_REF] Nouguier | Nonlinear ocean wave reconstruction algorithms based on spatiotemporal data acquired by a flash lidar camera[END_REF], we present and discuss the characteristics of an algorithm, aimed at assimilating such sparse data and able to deterministically reconstruct and propagate ocean surface elevations for their prediction in time and space. This algorithm could assist in the automatic steering and control of a variety of surface vehicles. Specifically, we compare prediction results using linear wave theory and the weakly nonlinear Choppy Wave Model [START_REF] Nouguier | choppy wave" model for nonlinear gravity waves[END_REF][START_REF] Nouguier | Scattering from nonlinear gravity waves: the "choppy wave" model[END_REF], extended here to an "improved" second order formulation. The latter model is based on an efficient Lagrangian formulation of the free surface and was shown to be able to model wave properties that are important to the proper representation of nonlinear free surfaces, namely wave shape and celerity. Synthetic datasets from highly nonlinear High Order Spectral simulations are used as reference oceanic surfaces. Predicted results are analyzed over an area that evolves in time, using the theoretical amount of information assimilated during the reconstruction of the wave field. For typical horizons of prediction, we discuss the capabilities of our assimilation process for each wave model considered.

INTRODUCTION

The availability of real-time phase-resolved wave fields is important to many offshore applications, such as the optimal maneuvering and operations of surface vessels or ocean renewable energy harvesting systems. However, many standard ocean wave models [START_REF] Booij | The swan ẅave model for shallow water[END_REF][START_REF] Chawla | A multigrid wave forecasting model: A new paradigm in operational wave forecasting[END_REF] only predict phase-averaged quantities based on spectral representations (e.g., directional wave energy spectrum), which do not provide a detailed description of deterministic wave properties (e.g., free surface elevation and fluid pressure fields). Nevertheless, under the linear wave superposition assumption, such spectral representations can be used to generate phase-resolved wave fields, by combining frequency wave component amplitudes with a set of (typically random) phases. Based on this data, the initial wave field can be computed through an inverse Fourier transform (FT) and its propagation in space and time performed using a linear or nonlinear propagation model (e.g., High Order Spectral (HOS) model [START_REF] Ducrozet | Hos-ocean: Open-source solver for nonlinear waves in open ocean based on high-order spectral method[END_REF]).

Here, we instead propose a method for directly reconstructing a phase-resolved nonlinear wave field, based on spatiotemporal optical measurements of the ocean surface. Future sea state conditions are then predicted by propagating the reconstructed waves in time and space using a weakly nonlinear model. This problem was studied in some earlier work on the basis of free surface elevation time series measured at fixed wave probes [START_REF] Blondel | Deterministic non-linear wave prediction using probe data[END_REF][START_REF] Takagi | Prediction of wave time-history using multipoint measurements[END_REF][START_REF] Simanesew | Surface wave predictions in weakly nonlinear directional seas[END_REF]. However, in situ measurements, typically made from a moving vessel or vehicle, are more challenging since wave reconstruction must rely on data acquired at constantly updated locations surrounding the path of the vehicle, which leads to practical limitations. A successful solution to this problem was proposed based on X-band radar measurements made from an onboard sensor, combined with a 3D-FT to reconstruct a large patch of free surface elevation surrounding the sensor [START_REF] Dankert | Wind-and wave-field measurements using marine x-band radar-image sequences[END_REF][START_REF] Xiao | Large-scale deterministic predictions of nonlinear ocean wave-fields[END_REF][START_REF] Hilmer | Deterministic wave predictions from the wamos ii[END_REF][START_REF] Qi | Phase-resolved wave field simulation calibration of sea surface reconstruction using noncoherent marine radar[END_REF]. This initial estimate was then used in a direct numerical simulation of future sea states.

Our work is based on the new nonlinear ocean surface reconstruction algorithm proposed by Nouguier et al. [START_REF] Nouguier | Nonlinear ocean wave reconstruction algorithms based on spatiotemporal data acquired by a flash lidar camera[END_REF][START_REF] Grilli | Ocean wave reconstruction algorithms based on spatio-temporal data acquired by a flash lidar camera[END_REF], on the basis of simulated free surface elevation datasets, such as would be acquired by an optical sensor that can remotely and simultaneously measure surface elevations at many geo-referenced locations at a high frequency (e.g., a generic LIDAR camera, stereovideo acquisition system). Their algorithm was based on applying the weakly nonlinear Choppy Wave Model [START_REF] Nouguier | choppy wave" model for nonlinear gravity waves[END_REF][START_REF] Nouguier | Scattering from nonlinear gravity waves: the "choppy wave" model[END_REF], which is Lagrangian and thus more efficient for propagating wave fields in the time domain than Eulerian models based on higher-order Stokes expansions. In this paper, we both validate and extend this approach by using simulation results from a HOS model as reference data. Linear wave propagation models can be used to provide short-term forecasts in calm sea states. However, as time between the prediction horizon and the initial sea state reconstruction increases and/or the sea state becomes more severe, nonlinear wave effects play an increasingly large role in the proper representation of wave fields, their propagation (e.g., due to am-plitude dispersion effects on wave celerity), and on the accuracy of the forecast. Using a nonlinear wave reconstruction and propagation algorithm thus becomes crucial to an accurate sea state prediction.

In this study, the performance of the weakly nonlinear Choppy Wave Model (CWM) and its improved second-order formulation developed hereby, are compared with predictions of the higher-order, but much more computationally demanding, "HOS-ocean" model [START_REF] Ducrozet | Hos-ocean: Open-source solver for nonlinear waves in open ocean based on high-order spectral method[END_REF]. Earlier work on ocean surface reconstruction and forecasting algorithms [START_REF] Wu | Direct simulation and deterministic prediction of large-scale nonlinear ocean wave-field[END_REF][START_REF] Blondel | Reconstruction et prévision déterministe de houle à partir de données mesurées[END_REF][START_REF] Qi | Predictable zone for phase-resolved reconstruction and forecast of irregular waves[END_REF] has shown that the size of the dataset assimilated during the reconstruction phase yields a well-defined time evolving spatial area, over which a sea-state prediction is theoretically possible. Accordingly, in this paper, we also evaluate the efficiency of our prediction method with respect to this theoretical prediction zone.

This paper is organized as follows. We first describe the wave models used for the free surface reconstruction and propagation. Second, for each wave model, we present the inverse problem that must be solved to perform the initial free surface reconstruction. Third, we detail the method for generating synthetic spatio-temporal free surface elevation datasets, such as would be acquired by a the depicted optical sensor. Fourth, for one-dimensional measurements used for simplicity, we analyze the effect on the forecast accuracy of reconstruction parameters, such as high cut-off frequency of the reconstructed spectrum, heterogeneity of the observation grid, and number of spatial observation times. For a two-dimensional surface reconstruction, we finally investigate the effect of directional wave spreading.

WAVE PROPAGATION MODELS

Two types of wave models are applied in this work to reconstruct ocean surfaces, based on a measured spatio-temporal dataset. These are based on: (i) linear wave theory (LWT); and (ii) the first-and second-order Lagrangian CWM, which include some intrinsic nonlinear wave properties. Note that the following developments are done under deep water assumption, but the extension to finite depth is straightforward.

Linear Wave Model

We consider a Cartesian coordinate system (x, y, z), with x and y axes located on the mean water surface and the z axis being vertical and positive upward. A linear ocean surface representation is derived by superposing N individual wave components, propagating in the horizontal plane r r r = (x, y), of amplitude A n , wavelength λ n and direction θ n with respect to the x-axis (n = 1, ..., N), yielding,

η (r r r,t) = N ∑ n A n cos (k k k n • r r r -ω n t + ϕ n ) , (1)
where t is time, ϕ n are random phases uniformly distributed in [0, 2π], and

k k k n = k n k k k n = (k n cos θ n , k n sin θ n )
where k n = 2π/λ n are wavenumbers related to angular frequencies ω n through the deep water dispersion relationship, ω 2 n = k n g, where g denotes the acceleration of gravity. To simplify the mathematical developments, we will use the equivalent formulation,

η (r r r,t) = N ∑ n (a n cos ψ n + b n sin ψ n ) , (2)
where

ψ n = k k k n • r r r -ω n t

Choppy Wave Model

Optical measurements of the ocean surface provide, at each time, the surface elevation measured for a set of spatial points at irregular but defined Eulerian locations in a reference coordinate system. Hence, the wave model used in the reconstruction algorithm must also be able to provide and use comparable information. On the other hand, Lagrangian representations of the free surface of a given order can model nonlinear wave properties that are not included in Eulerian developments of the same order [START_REF] Nouguier | choppy wave" model for nonlinear gravity waves[END_REF][START_REF] Nouguier | Scattering from nonlinear gravity waves: the "choppy wave" model[END_REF][START_REF] Gjosund | A lagrangian model for irregular waves and wave kinematics[END_REF]. Here, the CWM following our developments meets both of these requirements, as it is derived as a solution of Lagrangian dynamical equations in an Eulerian system.

Let R R R L and Z L denote the horizontal and vertical displacements of a particle with initial horizontal position r r r 0 and vertical position z 0 = 0 (i.e., on the free surface), respectively. A Lagrangian perturbation expansion in wave steepness of these parameters yields [START_REF] Pierson | Models of random seas based on the lagrangian equations of motion[END_REF][START_REF] Pierson | Perturbation analysis of the navierstokes equations in lagrangian form with selected linear solutions[END_REF][START_REF] Nouguier | Second-order lagrangian description of tri-dimensional gravity wave interactions[END_REF],

R R R L (r r r 0 ,t) = r r r 0 + ∑ i D D D i (r r r 0 ,t) , Z L (r r r 0 ,t) = ∑ i Z i (r r r 0 ,t) , (3)
where D D D i and Z i are the i th -order terms in the Lagrangian expansion. Expressing the particle coordinates from Eqn. (3) in an Eulerian framework requires finding an explicit relationship between R R R L and Z L , which can be achieved through the implicit relationship [START_REF] Grilli | Ocean wave reconstruction algorithms based on spatio-temporal data acquired by a flash lidar camera[END_REF],

Z L R R R L -∑ i D D D i (r r r 0 ,t) ,t ≃ Z L R R R L -∑ i D D D i (R R R L ,t) ,t . (4)
The CWM surface elevation based on this expansion are then reformulated as,

η (r r r,t) = ∑ i Z i r r r -∑ i D D D i (r r r,t) ,t . (5)
The zeroth-order Lagrangian solution is the particle position at rest, i.e., D D D 0 = (0, 0) and Z 0 = 0. The first-order solution yields [START_REF] Pierson | Models of random seas based on the lagrangian equations of motion[END_REF][START_REF] Pierson | Perturbation analysis of the navierstokes equations in lagrangian form with selected linear solutions[END_REF][START_REF] Nouguier | Second-order lagrangian description of tri-dimensional gravity wave interactions[END_REF],

D D D 1 (r r r,t) = ∑ N n k k k n (-a n sin ψ n + b n cos ψ n ) , Z 1 (r r r,t) = ∑ N n a n cos ψ n + b n sin ψ n . (6)
Assuming a single wave component, this first-order solution has been shown to be consistent with a third-order Stokes expansion for a surface frozen at t = 0 (see Eqn. (63) in [START_REF] Nouguier | choppy wave" model for nonlinear gravity waves[END_REF], with a misprint in the sign of the cos (2Kx) factor). Apart from the mean surface level this makes the first-order CWM efficient at representing nonlinear geometrical properties of ocean surface waves (i.e., sharper crests and flatter troughs). However, at this order, water particles located on the free surface move around circles around their initial position (r r r, z = 0) with a uniform angular velocity and, hence, their mean position is at z = 0. Thus, the asymmetry observed in nonlinear wave shape is not properly accounted, as it would require a non-zero mean surface level [START_REF] Nouguier | choppy wave" model for nonlinear gravity waves[END_REF][START_REF] Gjosund | A lagrangian model for irregular waves and wave kinematics[END_REF][START_REF] Longuet-Higgins | Eulerian and lagrangian aspects of surface waves[END_REF].

Increasing the order of the Lagrangian expansion, we find that the second-order term Z 2 includes a vertical shift of the mean water level [START_REF] Pierson | Models of random seas based on the lagrangian equations of motion[END_REF][START_REF] Nouguier | Second-order lagrangian description of tri-dimensional gravity wave interactions[END_REF], which causes particles to oscillate around a higher position. This important feature was added to the reconstruction process by reformulating the first-order CWM as,

η (r r r,t) = N ∑ n (a n cos ψn + b n sin ψn ) + 1 2 k n a 2 n + b 2 n , ψn = k k k n • (r r r -D D D 1 (r r r,t)) -ω n t. ( 7
)
which, throughout the paper, will be referred to as CWM1.

In addition, we also developed and used an improved formulation for a second-order CWM, denoted ICWM,

η (r r r,t) = N ∑ n a n cos Ψn + b n sin Ψn + 1 2 k n a 2 n + b 2 n , Ψn = k k k n • [r r r -(D D D I (r r r,t) + U s t)] -Ω n t, (8)
with,

D D D I (r r r,t) = N ∑ n k k k n (-a n sin Ψ n + b n cos Ψ n ) , Ψ n = k k k n • r r r -Ω n t, ( 9
)
where

Ω n = ω n -1/2k k k n • U s and U s = ∑ N n ω n k k k n a 2 n + b 2
n is the Stokes drift vector. Note that, for a monochromatic wave, this introduces corrections of the phase and group velocities consistent with third-order Stokes theory. In fact, here, we only consider self-interaction terms at second order, thus neglect the effect of interactions of different components. The detailed derivation of this formulation will be presented in future work [START_REF] Guérin | An improved "choppy wave" model for the time evolution of nonlinear surfaces[END_REF].

METHODS

We present a theoretical, yet realistic, setup for an optical system able to instantaneously remotely measure a large dataset of ocean surface elevations. This dataset is used in combination with various wave models (e.g., LWT, CWM1, ICWM) to perform a nowcast (i.e., initial) ocean surface reconstruction, and based on this to issue short-term forecasts of ocean surface elevations over a specified area. To this effect, we discuss the suitable choice of cutoff frequencies and directions of the reconstructed wave field, and based on this define the area over which the accuracy of the sea state forecast will be evaluated, as compared to a prediction made with a higher-order HOS model. We detail the free surface reconstruction process for linear and nonlinear wave models, and the regularization procedure that is required for solving the resulting inverse problem.

Setup Description

Let us consider an optical system (denoted by OS throughout this paper) simply mounted on a fixed offshore structure at location (x c , y c , z c ) (z c = 30 m is used in the present applications) and facing the ocean surface with viewing angles (α, β ) = (76 • , 0 • ) and aperture angles (α a , β a ) (here α a = 20 • ), and remotely measuring free surface elevations using J rays (here J = 64 × 64 = 4096), uniformly distributed over the aperture angles. This yields the observation zone shown in Fig. 1. During an assimilation time T a , the OS acquires K sets of spatial observations, at a constant sampling frequency f s (here 1 Hz), leading to J × K = L spatio-temporal observations.

We assume a random incident wave field, with its main di-rection along the x-axis, described by a directional JONSWAP (JS) spectrum, S (k, θ ) = F (k) × G (θ ) [START_REF] Hasselmann | Mea-surements of wind-wave growth and swell decay during the joint north sea wave project (jonswap)[END_REF], of significant wave height H s , peak period T p , and peakedness parameter γ. The directional spreading function is standard and defined as,

G (θ ) = 1 σ cos 2 π(θ -θ dir ) 2σ , if |θ -θ dir | ≤ π 2 , 0, if |θ -θ dir | > π 2 , (10)
where θ dir = 0 is the mean direction of propagation of the wave field, and σ is a directional parameter. The combination of wave field parameters selected in the present simulations is:

H s = 3 m, T p = 10 s, γ = 3.3, σ = π/4,

Cutoff Frequencies and Directions

In practical applications, the reconstructed wave field components should be limited in frequency and direction to the bandwidths that are meaningful to the dynamic response of the structure or ocean vehicle of interest. As the study of a specific structural response is beyond the scope of this paper, we first aim at accurately reconstructing wave fields without specifying any response bandwidth. Hence, the cutoff limits k min,max , θ min,max are calculated with respect to the reference ocean wave spectrum.

Additionally, both the size of the observation area and the spatial sampling resolution are related to the minimum and maximum wave frequencies that can be reconstructed, respectively. For the sake of simplicity, here, we only consider spatial observations (i.e., K = 1) to estimate cutoff frequencies. In this case, the low cutoff wavenumber (or frequency since ω = √ gk) is defined by the largest distance L obs between two observation points as, k min ≥ 2π/L obs . In the following, we set k min to its minimal value by determining L obs at the first observation time t obs . When reconstructing a signal over a regular observation grid (i.e., with constant spatial sampling frequency), the maximum high cutoff frequency must satisfy Shannon's condition, i.e., k max ≤ 2π/ (2ℓ obs ), where ℓ obs is the distance between two observation points. Since the observation grid is highly irregular, we will determine k max by evaluating its influence on the quality of the reconstruction.

Cutoff limits in wave directions are set by considering the fraction of total wave energy beyond which the remaining energy can be considered as negligible, i.e.,

θ min -∞ G (θ ) dθ = +∞ θ max G (θ ) dθ = µ +∞ -∞ G (θ ) dθ , (11)
with µ the fraction of total energy considered as negligible. In the following µ = 1%, which corresponds to θ min,max ≃ ±35 • for the spreading function defined above.

Prediction Zone Definition

Only a limited amount of information is assimilated during the reconstruction process, due to the physical spatio-temporal limitations of optical free surface observations (e.g., limited size of observation zone) and the finite bandwidths of the reconstructed wave field. Hence, only a similarly limited spatiotemporal region, called prediction zone, is accessible for sea state forecast, as a result of propagating the assimilated information [START_REF] Blondel | Reconstruction et prévision déterministe de houle à partir de données mesurées[END_REF][START_REF] Qi | Predictable zone for phase-resolved reconstruction and forecast of irregular waves[END_REF]. For a one-dimensional (1D) wave field, this region is bounded by the time evolution of the amplitude and phase of the slowest and fastest wave components at the end (x max obs ) and beginning (x min obs ) of the observation zone, respectively. It can be shown that reconstructed information attached to a specific wave component (i.e., wave amplitude and phase) is traveling at the wave group velocity C g = ω/ (2k) [START_REF] Wu | Direct simulation and deterministic prediction of large-scale nonlinear ocean wave-field[END_REF]. Thus, in 1D, as shown in Fig. 2a, a point (x,t ≥ t rec ) is included in the prediction zone if,

x min obs +C max g (t -t rec ) ≤ x ≤ x max obs +C min g (t -t obs ) , (12)
where t rec and t obs are the reconstruction and observation times, respectively. Both the LWT and CWM1 models rely on the linear deep water dispersion relationship, while ICWM includes a second order correction term which affects both phase and group velocity. The two constrains described by Eqn. ( 12) eventually cross each other (at point P max = (x max ,t max ); Fig. 2a), which means that the measured information is obsolete and can no longer be used in a prediction. Note that increasing the assimilation time T a leads to increasing the size of the prediction zone.

In two-dimensions (2D), the wave field is decomposed in individual wave components propagating in direction θ ∈ θ min , θ max . The corresponding prediction zone can be calculated for each direction using Eqn. [START_REF] Hilmer | Deterministic wave predictions from the wamos ii[END_REF], as if it were a 1D case, by replacing x obs by d obs = r r r obs • k k k = x obs cos θ + y obs sin θ , the distance along the considered direction. The intersection of each zone forms the 2D prediction zone. To save computational time, the prediction zone area can be approximated based on the two extreme directions θ min and θ max , as shown in Fig. 2b. The farthest reachable prediction zone P max corresponds to a segment (Fig. 2b).

Data Assimilation Process

Linear Reconstruction For a linear wave field represented by Eqn. (1), the reconstruction process is based on the minimization of a quadratic cost function representing the mean square difference between observations and predictions of the wave model. Here, assuming there are no measuring and model errors, the cost function is,

C = 1 L L ∑ ℓ=1 ( η (r r r ℓ ,t ℓ ) -η ℓ ) 2 , ( 13
)
where η (r r r ℓ ,t ℓ ) and η ℓ are the predicted and measured surface elevations at observation points r r r ℓ and at time t ℓ (ℓ = 1, ..., L), respectively. This problem is equivalent to a least-square minimization of a cost function, which can be achieved by specifying (for n, m = 1, ..., N; note that index summation is implied for repeated indices),

∂C ∂ a m = 0, ∂C ∂ b m = 0 ⇐⇒ A mn X n = B m , (14)
where,

X n = a n , X N+n = b n , (15)
is the control vector of 2N unknown wave parameters and,

B m = L ∑ ℓ=1 η ℓ cos ψ mℓ , B N+m = L ∑ ℓ=1 η ℓ sin ψ mℓ , (16)
is a vector containing the observation points information, with

ψ mℓ = k k k m • r r r ℓ -ω m t ℓ ,

and,

A mn = L ∑ ℓ=1 cos ψ mℓ cos ψ nℓ , A m,N+n = L ∑ ℓ=1 cos ψ mℓ sin ψ nℓ , A N+m,n = L ∑ ℓ=1 sin ψ mℓ cos ψ nℓ , A N+m,N+n = L ∑ ℓ=1 sin ψ mℓ sin ψ nℓ , (17)
is the wave model matrix. The linear system of Eqns. ( 14) to [START_REF] Qi | Predictable zone for phase-resolved reconstruction and forecast of irregular waves[END_REF] can then be solved for the wave model parameters (a n , b n ).

Nonlinear Reconstruction For a nonlinear sea state reconstruction the cost function of Eqn. ( 13) is also minimized, using one of the nonlinear wave models, CWM1 or ICWM, to represent surface elevations (Eqns. [START_REF] Blondel | Deterministic non-linear wave prediction using probe data[END_REF] or [START_REF] Takagi | Prediction of wave time-history using multipoint measurements[END_REF]). Here, we detail the system construction for CWM1. Following the same methodology as for the linear model, we use the cost function minimization Eqn. ( 14) with,

B m = L ∑ ℓ=1 η ℓ P mℓ , B N+m = L ∑ ℓ=1 η ℓ Q mℓ , (18)
with,

P mℓ = (1 + b m k m sin ψ mℓ ) cos ψmℓ -a m k m sin ψ mℓ sin ψmℓ + k m a m , Q mℓ = (1 + a m k m cos ψ mℓ ) sin ψmℓ -b m k m cos ψ mℓ cos ψmℓ + k m b m . (19)
Similarly, we have,

A mn = L ∑ ℓ=1 P mℓ cos ψnℓ + 1 2 k n a n , A m,N+n = L ∑ ℓ=1 P mℓ sin ψnℓ + 1 2 k n b n , A N+m,n = L ∑ ℓ=1 Q mℓ cos ψnℓ + 1 2 k n a n , A N+m,N+n = L ∑ ℓ=1 Q mℓ sin ψnℓ + 1 2 k n b n . (20)
Since both A mn and B m now depend on wave parameters (a n , b n ), the system of Eqns. ( 14) becomes nonlinear and is solved iteratively. To do so, equations are linearized by computing A 16) and [START_REF] Qi | Predictable zone for phase-resolved reconstruction and forecast of irregular waves[END_REF]. Convergence is typically achieved within 5 to 10 iterations.

A similar approach can be used for ICWM, but the inverse problem becomes more complex due to the higher-order of expansion in the formulation. The complete derivation of the solution for this model and its reconstruction system will be presented in future work.

Regularization

In applications, the inverse reconstruction problem defined above (Eqn. ( 14)) can become illconditioned due to practical constraints, such as the heterogeneous distribution of spatial observation points, the limited ocean area observed by the OS, and the the frequency and direction bandwidth cutoffs in the reconstructed wave field. To be able to find consistent results independently of the conditioning of the matrix to invert (i.e., mn ), a Tikhonov regularization method is applied, which yields,

min ||A mn X n -B m || 2 -r 2 ||X m || 2 , ( 21
)
where r is the regularization parameter and ||.|| denotes the Euclidean norm. The optimal regularization parameter is found using the "L-curve" method, which consists of finding the r value corresponding to the point of maximum curvature (i.e., corner) of the (log ||A mn X n -B m || , log ||X m ||) function, hence, providing an optimal compromise between minimizing the residual error and ensuring that the norm of the solution does not become too large. The L-curve corner can be determined analytically through solving a singular value decomposition problem [START_REF] Calvetti | L-curve and curvature bounds for tikhonov regularization[END_REF][START_REF] Hansen | Computational inverse problems in electrocardiology[END_REF].

AND DISCUSSION

In the following applications, we illustrate and discuss key of our proposed ocean wave reconstruction algorithms, for optically measured surface elevation data. In the absence of field data, the algorithms are applied to reference synthetic surface elevation datasets, generated first in 1D using a higher-order HOS model. The influence of reconstruction parameters on the accuracy of the solution with respect to the reference data is studied, in particular, the number of wave components, the high cutoff frequency, the heterogeneity of the observation grid, and the number of observation times, while highlighting differences between linear and nonlinear prediction results. We then present a 2D case and show how the horizontal aperture angle can help in the reconstruction.

Synthetic Dataset and Error Definition

The open source HOS-ocean model [START_REF] Ducrozet | Hos-ocean: Open-source solver for nonlinear waves in open ocean based on high-order spectral method[END_REF] is used to generate synthetic reference ocean surfaces from which datasets of observation points are extracted. Since its first developments [START_REF] Dommermuth | A highorder spectral method for the study of nonlinear gravity waves[END_REF][START_REF] West | A new numerical method for surface hydrodynamics[END_REF], this model, which simulates the nonlinear propagation of directional wave fields over large spatio-temporal domains, has undergone numerous validations and been used in various applications, such as for simulating freak waves [START_REF] Ducrozet | 3-d hos simulations of extreme waves in open seas[END_REF], developing a numerical wave tank [START_REF] Bonnefoy | Advances in Numerical Simulation of Nonlinear Water Waves[END_REF], or in wave prediction [START_REF] Blondel-Couprie | Experimental validation of non-linear deterministic prediction schemes for long-crested waves[END_REF][START_REF] Qi | Nonlinear phase-resolved reconstruction of irregular water waves[END_REF]. Based on an initialization with a JS spectrum and a set of random phases, HOSocean computes the time evolution of the complex amplitude of wave components. These are then used to construct snapshots of the reference ocean surface at selected observation times. In the presented case, the order of the HOS expansion is set to 1 or 5, to generate a linear or a nonlinear reference surface, respectively. In each case, the HOS numerical grid and parameters are selected to ensure converged results.

Based on the computed reference surfaces, synthetic datasets of surface elevations data such as would be acquired by the depicted OS are generated by computing the intersection between J optical rays with the reference ocean surface (see, e.g., Fig. 3 for a 1D surface).

The following relative standard error (RSE) is used as an indicator of the prediction accuracy,

ε RSE (t) = ∑ N r q=1 (η (r r r q ,t) -η q ) 2 ∑ N r q=1 η 2 q , ( 22
)
where r r r defines the footprint of the spatial prediction zone at time t and N r is the number of reconstructed points, which is selected large enough to achieved a converged error estimate. RSE quantifies the scaled mean squared difference between the reconstructed and reference surfaces over the observation points in the prediction zone. Since the quality of the prediction depends on the local surface geometry at the observation time (e.g., through shadowing effects), an overall RSE is computed by averaging RSE values for N s = 50 surface reconstructions of ocean surface generated at the same time with HOS-ocean for different initial random phases,

εRSE (t) = 1 N s N s ∑ q=1 ε RSE q (t) , ( 23
)
using the same initial wave spectrum.

One-Dimensional Reconstruction

Using a JS spectrum with parameters as defined above, 1D linear surfaces are first generated and then reconstructed based on data acquired at a single time (i.e., K = 1), with L = 64 observation points uniformly distributed over the observation area. Figure 4 shows the RSE at reconstruction time t rec = t obs = 0 (i.e., for a nowcast), as a function of the number of wave components N used in the reconstructed wave field, for different high cutoff frequencies. Note that Shannon's criterion yields, k max ≤ 13.9k

p , but the HOS surface was generated using a cutoff frequency k max ≃ 16k p . As expected, results show that the reconstruction error decreases as the high cutoff frequency grows, and eventually converges to a minimum, here approximately 10%, for k max = 10k

p . In all cases, the number of wave components (N = 15 to 60) appears to have a very limited influence on the accuracy of the reconstruction. However, this number must be

FIGURE 5. PREDICTION ERROR OF A NONLINEAR OCEAN SURFACE AS A FUNCTION OF TIME, FOR DIFFERENT WAVE RECONSTRUCTION MODELS, OVER A REGULAR GRID FOR

K = 1.
large enough to allow for a proper inversion of the reconstruction system. We find: N ≥ 30 for k max ≥ 8k p and N ≥ 40 for k max ≥ 12k p . In the following we use k max = 10k p and N = 50. Next, nonlinear reference surfaces are generated using the same parameters and the reconstruction error over the prediction zone is computed as a function of time t = t rec = 0 to 3T p ≃ 0.8t max , for the same L = 64 uniformly distributed observation points, using either the LWT, CWM1, or ICWM model in the reconstruction. For each model, the ocean surface is first reconstructed at t = t rec = 0 and then propagated to a later time where its RSE error with respect to the reference HOS solution is evaluated; for comparison, the case of Fig. 4 propagated in time is also marked on the figure (i.e., reconstruction using LWT of a linear reference surface). Figure 5 shows that at t = t rec = 0 all models give the same results, which is expected since at initial time the different reconstruction algorithms all solve equivalent inverse problems. Then, as time increases, each wave model behaves according to its own properties. While the CWM1 model provides a better prediction than the LWT model, the ICWM further improves the prediction accuracy as time increases, due to the additional nonlinear phase speed corrections. Here, the improvement of the nonlinear over the linear prediction is on the order of a few percent. The linear prediction of a linear reference surface, shown on the same figure, is an indication of the best reconstruction accuracy that could be achieved for this problem and dataset, if all nonlinear effects were properly accounted for in the nonlinear wave models. Hence, the difference between this linear and nonlinear forecasts using ICWM increases with time due to the incomplete representation of time-dependent nonlinear effects in the latter model, to approximately reach 5% at t = 3T p . Ocean observations made by an OS such as a LIDAR camera yield highly irregular grids, since the density of observation points geometrically decreases with the distance from the OS, which is made worse by wave shadowing effects becoming more important at grazing incidence angles (i.e., for the most distant observation points; see Fig. 3). The effect of using an irregular grid on the accuracy of the ocean surface reconstruction is assessed by recalculating the case of Fig. 5, using the ICWM model, for a dataset of L = 64 irregularly spaced observation points. First, considering a single observation time (i.e., K = 1), Fig. 6 shows that the reconstruction error greatly increases (by more than a factor of 2) as compared to a regular dataset. As suggested in earlier work [START_REF] Nouguier | Nonlinear ocean wave reconstruction algorithms based on spatiotemporal data acquired by a flash lidar camera[END_REF], however, this deterioration of the ocean reconstruction and forecasting can be significantly reduced by increasing the number of observation times. Indeed, this has the effect of making some portions of the free surface, invisible to the OS at a given time due to shadowing effects, visible and hence measurable at a later time. Thus, for K = 2, Fig. 6 shows that the error already drops by 10% as compared to K = 1, and then further drops when increasing K to reach a maximum of ∼ 20% for K = 10. The latter, however, is still 5% larger than the error achieved when using a regular grid, for the same number of observation times (also shown in the figure). Comparing results in Figs. 5 and6 we also see that, when using a regular grid of observation points, increasing K has no significant effect on the prediction error; for both the LWT and ICWM models, maximum errors indeed remain on the order of 15 to 20% for K = 1 and 10, even if ICWM leads to a consistently better prediction. Figure 6 also shows that the reconstruction error decreases when K increases form 1 to 10, but then gradually increases for larger K values. This might result from increasing numerical errors in the inverse problem solution, as the observed dataset becomes larger. Moreover, since reconstructed waves are only modeled up to second-order in nonlinearity in ICWM, higherorder time-dependent nonlinear effects included in the HOSocean solution, which affect the reconstruction process when using temporal data, are not well represented in the wave model as time increases. Therefore, the model becomes unable to properly represent the ocean surface when the assimilation time T a is too long (e.g., phase shifts resulting from amplitude dispersion effects). Varying the temporal sampling frequency (set here to 1 Hz) would allow finding out whether it is the assimilation time or the total number of observations L that is more likely to improve, and then limit, the quality of the prediction.

K = 1 K = 2 K = 3 K = 5 K = 10 K = 20 K = 30 K =
In general, increasing the assimilation time T a results both in a decrease of the prediction error and in an increase of the size of the prediction area. Both effects are illustrated in Fig. 7, which shows the reference and reconstructed surface elevations at reconstruction time t rec , for K = 1 and 30. Observation points are shown at the first observation time t obs . With K = 30, the reconstructed surface is much closer to the reference data than for K = 1, and the prediction area is larger. The figure shows that, for K = 30, the reconstructed surface remains quite close to the reference data even beyond the prediction zone (x/λ p ≥ 3.8). This suggests that the estimate for the upper limit of the prediction zone, which is calculated based on the distance traveled by the slowest wave component during the assimilation time, is conservative (Eqn. ( 12)). Since longer waves travel faster, the prediction of low frequency components can still be accurate beyond the upper limit. A similar effect affects the high frequency components at distances below the lower limit, when time increases (t > t rec ). This suggests the presence of a "gray area" surrounding the prediction zone and growing with time, over which the amount of measured information could still be enough for an accurate prediction, but with decreasing accuracy as distance increases from the prediction zone boundary.

In Fig. 7, when the density of observation points decreases (K = 1; 1.2 ≤ x/λ p ≤ 2.3), the reconstructed surface tends not to fit well the reference surface elevation. As Shannon's condition is not met over this part of the observation zone, an aliasing phenomenon leads to overestimating high frequency component amplitudes, which causes the observed high frequency surface oscillations. However, since this poorly reconstructed part of the observation zone is located close to the lower limit of the prediction zone, and the inaccuracy only affects high frequency components, the growing "gray area" effect tends to remove that part of the surface from the prediction zone, while the other part remains within it. This explains the decreasing prediction error as a function of time for small K values observed in Fig. 6. With several observation times, this phenomenon disappears due to an artificial refinement of the observation grid: as indicated before, parts of the surface that were initially not illuminated by the OS rays become visible at a different observation time.

Two-Dimensional Reconstruction

We now generate a 2D nonlinear wave field with HOSocean, for the same input JS spectrum as before, and reconstruct it using an observation grid of 64 × 64 spatial points (L = 4096) with k max = 10k p , K = 10, and 50 wave frequencies distributed over 30 directions, yielding N = 1500 wave components. Since the number of spatial observations is significantly larger than before, calculating the exact intersection of L optical rays with time-varying reference surfaces becomes numerically demanding. Thus, we estimate the horizontal location of observation points, r r r ℓ , as if the reference surface were plane, neglecting shadowing effects. Results are comparable since the major contribution to the spatial spreading of observation points, i.e., the increasingly grazing incidence of OS rays at large distances, is accounted for. Also, the reconstruction error is now calculated for one surface and not averaged over 50 surfaces as in 1D, and the study presented here is limited to the linear wave reconstruction, even if the reference oceanic surface is nonlinear.

In the following, we estimate the effect of the directional spreading of the observation grid on the ocean reconstruction, which is function of the horizontal aperture angle β a of the OS (Fig. 1). Figure 8 shows the reconstruction error calculated as a function of time over the 2D prediction zone. With β a = 10 • , the error is ∼ 70% at reconstruction time and almost monotonically decreases as time increases up to t = 3T p . Similar to 1D cases, this decrease in error results from the growing "gray area" surrounding the prediction zone boundary. Poorly reconstructed parts of the ocean surface close to the initial prediction zone boundary are gradually removed as time increases. However, with the decreasing size of the prediction zone, the error estimate becomes more subject to slight variations, depending on which waves are crossing the prediction zone at the considered time; this explains the increasing error at t = 2.5T p . Such effects did not occur in the 1D cases since errors were averaged over 50 different surfaces. Figure 8 shows that increasing β a results in a significantly reduced error, which converges to a minimum on the order of 15% for β a = 110 • ; this error is similar to 1D cases. This large error reduction is related to the area of the ocean surface covered by the observation zone, as compared to the area of the prediction zone, at reconstruction time. As shown in Fig. 2b, the two-dimensional prediction zone depends on the range of directions of propagation of wave compo-nents in the reconstructed wave field. Here, as θ min,max ≃ ±35 • and the OS is facing the main wave direction, the optimal aperture angle that leads to a maximal prediction zone coverage is 180 -|θ max | + θ min ≃ 110 • .

Figure 9 shows the horizontal distribution of the local error |η -η HOS | /H s , where η is the reconstructed linear surface and η HOS the nonlinear reference surface, at reconstruction time and after two peak periods of propagation, using an aperture angle β a = 110 • . The prediction error is clearly much smaller within the prediction zone than outside of it. At both selected times, the "gray area" with reduced local error, which surrounds the prediction zone boundary, can clearly be identified. At reconstruction time, the "gray area" is close to the the upper limit of the prediction zone (right side of the boundary in Fig. 9 top) because the reconstructed longer waves propagate faster, and thus reach beyond this limit during assimilation time. After 2T p , the "gray area" has expanded and is now completely surrounding the prediction zone. This confirms that the method used here for estimating the prediction zone boundary is conservative.

CONCLUDING REMARKS

We showed that both 1D and 2D ocean wave fields (or sea states) can be accurately reconstructed on the basis of irregular datasets of realistic, but synthetic (and thus error free), optical measurements. Reconstruction algorithms based on linear and various nonlinear wave models were applied and validated. Despite its simple analytical formulation, the Choppy Wave Model (CWM) is able to simulate relevant nonlinear effects in ocean wave fields. The new improved Choppy models that better predict wave shape (CWM1) and nonlinear phase velocity (ICWM), allow improving the accuracy of the wave reconstruction and prediction.

We quantified effects on the reconstruction of the main characteristic of optical measurements, which is to yield highly sparsely distributed datasets, making for a particularly challenging spatial reconstruction. However, we showed that using spatio-temporal information allows overcoming this issue.

We quantified the influence of several reconstruction parameters, which are shown to be of significant importance to the prediction error. In particular, adequately selecting the cutoff limits of the reconstructed spectrum (in both frequency and direction) and the optical system viewing angles, allows accurately predicting the wave field up to several peak periods of propagation. The influence on the prediction results of a growing "gray area" surrounding the prediction zone boundary was highlighted.

The numerical algorithms used in the present work were not optimized for computational efficiency nor parallelized. On the basis of theoretical and practical estimations, it is expected that optimized and efficiently parallelized algorithms (e.g., on multicore or GPU processors) will easily cope with the increased complexity and computational demand resulting from including non-linear effects in a 2D ocean reconstruction. At present, however, the assessment of such computational performance of the method is beyond the scope of this paper and will be left out for future work.

FIGURE 1 .

 1 FIGURE 1. GEOMETRICAL PARAMETERS LEADING TO THE OBSERVATION ZONE.

FIGURE 2 .

 2 FIGURE 2. (a) SPATIO-TEMPORAL EVOLUTION OF THE PRE-DICTION ZONE FOR A 1D SET OF OBSERVATIONS, AND (b) LO-CATION OF THE PREDICTION ZONE FOR A 2D SET OF OBSER-VATIONS FOR A SINGLE TIME RECONSTRUCTION.

m

  based on wave parameters obtained at the previous iteration p, when solving for X (p+1) n at iteration p + 1. The solution is initialized at p = 1, using A (0) mn and B (0) m based on the linear reconstruction problem and thus on Eqns. (

FIGURE 3 .

 3 FIGURE 3. 1D SPATIAL SAMPLING OF A SYNTHETIC NONLINEAR FREE SURFACE. STRAIGHT LINES ARE OS RAY TRAJECTORIES.

FIGURE 4 . 1

 41 FIGURE 4. LINEAR RECONSTRUCTION ERROR OF A LINEAR SURFACE AT TIME t rec = t obs = 0 (NOWCAST), AS A FUNCTION OF THE NUMBER OF WAVE COMPONENTS, FOR DIFFERENT HIGH CUTOFF FREQUENCIES, OVER A REGULAR GRID FOR K = 1.

FIGURE 6 .

 6 FIGURE 6. PREDICTION ERROR OF A NONLINEAR OCEAN SURFACE AS A FUNCTION OF TIME, FOR DIFFERENT NUM-BERS OF OBSERVATION TIMES, USING ICWM OVER AN IR-REGULAR GRID. SOME REGULAR GRID RESULTS ARE PLOT-TED FOR REFERENCE.

FIGURE 7 .

 7 FIGURE 7. ELEVATION AT RECONSTRUCTION TIME OF THE REFERENCE (HOS) AND RECONSTRUCTED (ICWM) SURFACES, FOR DIFFERENT NUMBERS OF OBSERVATION TIMES. VERTICAL LINES DEFINE THE SPATIAL LIMITS OF THE INITIAL PREDICTION ZONE.

FIGURE 8 .

 8 FIGURE 8. PREDICTION ERROR OF A NONLINEAR OCEAN SURFACE AS A FUNCTION OF TIME, USING ICWM FOR DIF-FERENT HORIZONTAL APERTURE ANGLES.

FIGURE 9 .

 9 FIGURE 9. LOCAL 2D LINEAR PREDICTION ERROR |η -η HOS | /H s , AT RECONSTRUCTION TIME (TOP) AND AFTER 2T p OF PROPAGATION (BOTTOM). LINES MARK THE BOUNDARY OF THE PREDICTION ZONE.

  are spatio-temporal phases, and (a n , b n ) = (A n cos ϕ n , A n sin ϕ n ) are wave parameters describing the ocean surface.

  which corresponds to a realistic moderately steep sea state, frequently observable in open oceans, with a corresponding, deep water characteristic steepness, H s /λ p ≃ 2%.

ACKNOWLEDGMENT

This work is supported by Naval Group and Centrale Nantes through their Joint Laboratory for Marine Technology, and by the PREDEMO-NAV project, grant ANR-15-ASTR-0006. The authors would like to acknowledge the efforts of reviewers for their constructive comments.


Filip Pawłowski
email: filip.pawlowski@ens-lyon.fr

Rob H Bisseling
email: r.h.bisseling@uu.nl

Bora Uc ¸ar ¶

A N Yzelman

Combinatorial Tiling for Sparse Neural Networks

Sparse deep neural networks (DNNs) emerged as the result of search for networks with less storage and lower computational complexity. The sparse DNN inference is the task of using such trained DNN networks to classify a batch of input data. We propose an efficient, hybrid model-and data-parallel DNN inference using hypergraph models and partitioners. We exploit tiling and weak synchronization to increase cache reuse, hide load imbalance, and hide synchronization costs. Finally, a blocking approach allows application of this new hybrid inference procedure for deep neural networks. We initially experiment using the hybrid tiled inference approach only, using the first five layers of networks from the IEEE HPEC 2019 Graph Challenge, and attain up to 2× speedup versus a data-parallel baseline.

I. INTRODUCTION

The Graph Challenge encourages new developments in graph analytics, which, in its 2019 edition [START_REF] Kepner | Sparse Deep Neural Network Graph Challenge[END_REF] focused on fast inference for sparse deep neural networks. It described the problem using the GraphBLAS C API [START_REF] Buluc | The GraphBLAS C API specification[END_REF], which advocates the use of matrix-based algorithms. Sparsely-connected neural networks exhibit lower computational complexity and lower memory requirements compared to their dense counterparts. They may originate by pruning a dense network as in the Banded Sparse Neural Networks [START_REF] Van Der Lugt | Banded Sparse Neural Networks and their parallel computation[END_REF], or result from training a fixed sparse topology as in the RadiX-Net [START_REF] Alford | Pruned and structurally sparse neural networks[END_REF]. The input data matrix may also be sparse, due to feature extraction techniques generating sparse representations (from, e.g., image, video, or signal data), or because input may be naturally sparse (e.g., graph inputs). One dataset with sparse representations is the MNIST database of handwritten digits [START_REF] Lecun | MNIST handwritten digit database[END_REF].

We summarise the sparse DNN challenge in terms of input and neural networks (NNs) used. A sparse neural network consists of d ∈ N sparse weight matrices W (k) ∈ R n k ×n k+1 called layers, where k, n k , n k+1 ∈ N for all 0 ≤ k < d. A feature matrix X (0) ∈ R n×n0 consists of n ∈ N sparse feature vectors, one for each data instance to be classified. Sparse inference refers to the computation of the final classification matrix X (d) ∈ R n×c from the input feature matrix X (0) :

X (d) = f (• • • f (f (X (0) W (0) + e n b (0) T
)W (1) + e n b (1)

T ) • • • W (d-1) + e n b (d-1) T ) , (1)
where c = n d is the number of classes, f : R → R is an activation function to be applied element-wise, b (k) ∈ R n k+1 ×1 is a vector of bias at layer k, and e n = (1, . . . , 1) T .

The challenge's input feature matrix consists of 60 000 images from the MNIST dataset with each image flattened to a single vector and thresholded such that the values are either 0 or 1. Images are interpolated to the number of neurons in the neural networks: 1024, 4096, 16384, and 65536. Several deep sparse neural networks are generated using RadiX-Net [START_REF] Alford | Pruned and structurally sparse neural networks[END_REF], with the number of neurons, layers, and bytes in Table I. This size, in bytes, assumes Compressed Row Storage (CRS) using four-byte values and indices. Inference employs the rectified linear unit (ReLU): for Y = f (X) and X, Y ∈ R m×n ,

Y ij = f (X ij ) = 0 if X ij < 0 X ij , otherwise , ∀0 ≤ i < m, 0 ≤ j < n.
The bias b in Equation ( 1) is -0.3,-0.35,-0.4,-0.45 for each of the 1k, 4k, 16k, and 64k neuron NNs, respectively. This correction can be integrated with the application of f . Without loss of generality we hence omit b from the remainder text. Sparse inference may be viewed as a repeated sparse matrix-sparse matrix multiplication (SpGEMM) C = AB with sparse matrices A, B, and C where nonzero output elements are filtered using the function f . Thus, our goal is to compute X (k+1) = f (X (k) W (k) ) for each layer. Parallelization strategies typically work by partitioning the input feature matrix (data-parallel), partitioning the neural network (modelparallel), and/or by pipelining [START_REF] Ben-Nun | Demystifying parallel and distributed deep learning: An in-depth concurrency analysis[END_REF]. The previous sparse DNN challenge submissions exploit data-parallelism and propose methods to maintain load balance since the number of nonzeros in the rows of the input matrix differ arbitrarily. There were six submissions [START_REF] Bisson | A GPU implementation of the sparse deep neural network graph challenge[END_REF]- [START_REF] Mofrad | Multithreaded layer-wise training of sparse deep neural networks using compressed sparse column[END_REF] providing performance data for analysis [START_REF] Kepner | Graphchallenge.org sparse deep neural network performance[END_REF]. Most employ high-performance frameworks, such as SuiteSparse:GraphBLAS, GraphBLAST, or Kokkos, to take advantage of highly-optimized SpGEMM kernels and achieve shorter development time. Since such frameworks are designed for (trans)portability, their use typically disallows optimizations across the individual steps of the inference.

Our approach assumes a shared-memory architecture consisting of p s connected processors. Each supports p t threads, for a total of p = p s p t threads. We pin threads to specific cores, so to prevent it moving from one core to another during the computation. To analyse the shared-memory algorithms, we distinguish between two types of data movement, intra-socket and inter-socket. These are quantified via the intra-socket cost

(k+1) = X (k) W (k) . WE ASSUME f HAS BEEN APPLIED AFTER ALL LAYERS, EXCEPT FOR X (0) .
g and the inter-socket cost h, both in seconds per byte. We write L (in seconds) for the time to complete a barrier.

The main contribution of this paper is twofold. First, we introduce a hypergraph model of the sparse neural network inference process. Second, we propose an efficient sharedmemory hybrid algorithm which combines the data-and model-parallel methods to achieve tiling through multiple layers of a neural network. The latter is enabled by the use of existing hypergraph partitioners, yielding an effective partitioning that leverages the density of weight matrices. We evaluate our techniques on the data and neural networks defined by the sparse DNN challenge.

A. The challenge dataset

We extract some statistics of the 1k-neuron, 120-layer network to demonstrate important properties of the weight matrices. Let the density of an m × n matrix A be nz (A) mn , and the compressed density be nz (A) m n , where nz(•) returns the number of nonzeros in a matrix and m is the number of nonempty rows. Table II shows the density and the compressed density of the feature matrices X (k) for a few of the first fifteen layers of the network, as well as the size of the memory touched by the SpGEMM operation processing that layer. The density of feature matrices continuously decreases from the second layer onwards and remains constant at 3.02% after the 14th layer. On the other hand, the compressed density increases and reaches 100%, and remains so after layer 14. We believe that compressed densities tend to hundred percent in most neural networks beyond RadiX-Net since classification tends to result in nonzero probabilities for all classes.

B. Sequential SpGEMM

There are a number of efficient SpGEMM algorithms that differ in the way the matrices are stored, the nonzeros are visited, and the way scalar multiples are accumulated [START_REF] Nagasaka | High-performance sparse matrix-matrix products on Intel KNL and multicore architectures[END_REF]. Our own sequential codes, alike Davis et al.'s [START_REF] Davis | Write quick, run fast: Sparse deep neural network in 20 minutes of development time in SuiteSparse:GraphBLAS[END_REF], use a variant of Gustavson's algorithm with modifications to fuse the filtering step (ReLU) and thresholding of the nonzeros. We replace the index list typically used in a sparse accumulator SPA with a dense array SPAC, which flags columns belonging to the active row i. Contrary to an index list, this array need not be reset between iterations and is more efficient when rows of the output matrix hold relatively many nonzeros-which, as argued in Section I-A, holds for NN inference. Finally, since our inference algorithms require two consecutive SpGEMMs with a common input matrix, we implement a fused variant to prevent accessing that matrix twice.

The work complexity of X (k+1) = X (k) W (k) using Gustavson's approach is Θ(

n-1 i=0 nz(X (k) {:,i} ) nz(W (k)
{i,:} )). The CRS matrices of X (k) and X (k+1) are streamed to the CPU only once, which incurs a data movement of Θ(nz (X (k) ) + nz (X (k+1) ) + n). The CRS of W (k) as well as the sparse accumulator structures (SPA and SPAC of size Θ(max k n k+1 )) are accessed Θ(nz(X (k) )) times; the whole W (k) is streamed at most once per row of X (k) , and only once at minimum.

C. Sequential inference

The work complexity of the sparse inference (1) derived by consecutively applying the preceding SpGEMM analysis is

Θ d-1 k=0 n k -1 i=0 nz(X (k) {:,i} ) nz(W (k)
{i,:} ) . Recall that while input and output feature matrices are touched only once, the weight matrix W (k) is touched multiple times. If this matrix can be cached, the data movement cost is

Θ d k=0 nz(X (k) ) + n + d-1 k=0 nz(W (k) ) + n k g . If W (k) cannot be cached this increases with O( d-1 k=0 n k -1 i=0 nz(X (k) {:,i} ) nz(W (k)
{i,:} )g). We use two buffers to store input and output feature matrices throughout inference, as each SpGEMM may reuse the preceding SpGEMM's input buffer to store its output while re-using the preceding output buffer as input. The total storage requirement thus becomes

Θ max k nz(X (k) ) + d-1 k=0 nz(W (k) ) + n + d k=0 n k .

II. APPROACH

We present the data-parallel inference, which is the defacto standard in deep learning in Section II-B. We then introduce the model-parallel approach in II-C which we extend to include tiling in Section II-D. Finally, we combine the tiling model-parallel inference with the data-parallel inference and describe how to apply the resulting method for deep neural networks in Section II-E. Our model-parallel, tiling, and hybrid tiling inference methods all rely on hypergraph modeling of the inference process. We introduce this first.

         W (0) T W (1) 0 • • • 0 0 W (2) T W (3) • • • 0 . . . . . . . . . 0 • • • W (d-4) T W (d-3) 0 0 • • • 0 W (d-2) T W (d-1)          Figure 1.
The staircase matrix of a neural network consisting of d layers.

(Here, d is even.)

A. Hypergraph models for inference

A hypergraph H = (V, N ) consists of a set of vertices V and a set of hyperedges N . A hyperedge (a net) h ∈ N is a subset of vertices, h ⊆ V. The hypergraph partitioning problem [START_REF] Lengauer | Combinatorial Algorithms for Integrated Circuit Layout[END_REF] is the task of dividing the vertices of a hypergraph into K parts of roughly equal sizes, while minimizing an objective (cost) function defined in terms of nets; the most common variants of this problem are NP-complete. That is one seeks a partition Π V = {V 0 , . . . , V K-1 } of the vertex set of H. This simultaneously induces a (K + 1)-way partition on the nets, Π N = {N 0 , . . . , N K-1 ; N S }, where nets in N i for 0 ≤ i < K have vertices only in V i and are internal, while nets in N S have vertices in more than one part and are external.

The connectivity λ h of a net h is equal to the number of parts it connects. A partition Π V is said to be balanced if for each part V k , w(V k ) ≤ w(V/K)(1 + ε),
where ε is a constraint to the partitioner, and w is the amount of work involved (usually measured as nonzeros in a corresponding matrix). C ¸atalyürek and Aykanat [START_REF] Aykanat | A fine-grain hypergraph model for 2d decomposition of sparse matrices[END_REF] proposed a fine-grain hypergraph model for sparse matrices A, under which each nonzero of A becomes a unique vertex while each column and row of A corresponds to a net consisting of the nonzeros it contains.

We contribute a hypergraph model of the whole neural network inference process which can be interpreted as applying the fine-grain model on the staircase matrix defined in Figure 1. This matrix connects the consecutive layers W (k) into a staircase, thus preserving the relationship between the individual layers. A partitioning of this hypergraph results in a distribution of weight matrices, thus defining which thread processes which part of the neural network. We do this under two constraints: 1) a load balance criterion on weights in each individual layer, which ensures threads have a roughly equal amount of work when collaboratively processing a layer, and 2) minimization of the number of external rows and columns in N S (which is known as the cutnet metric [START_REF] Aykanat | PaToH: A Multilevel Hypergraph Partitioning Tool[END_REF]). We use PaToH [START_REF] Aykanat | PaToH: A Multilevel Hypergraph Partitioning Tool[END_REF], a multi-constraint hypergraph partitioner, to obtain such a balanced K-partitioning of the staircase matrix.

It may be that a column of W (i) is internal while its corresponding row of W (i+1) is external, or vice versa. This creates dependencies between internal and external parts of the resulting partitioning of layers. Handling these efficiently during inference is possible using known techniques (see, e.g., Yzelman and Roose [START_REF] Yzelman | High-level strategies for parallel sharedmemory sparse matrix-vector multiplication[END_REF]); we save that exercise for future work, however. Instead and for simplicity, we here artificially move the internal rows or columns to the external set N S so that such conflicts will not arise in the subsequent presentation.

B. Data-parallel inference

Data-parallel inference partitions the feature matrix X (0) ∈ R n×n0 into p parts using a row-wise partitioning. We map row i to thread i/ n/p to achieve this, so that each thread performs sequential inference on different sets of (approximately) n/p input data elements in an embarrassingly parallel fashion.

Each thread q allocates two buffers to store the local feature matrices X (k) {sq,:} and X (k+1) {sq,:} when processing W (k) , where s q is the set of data indices assigned to thread q, while the weight matrices are allocated once and shared. This approach has no work overhead compared to a sequential method. There is, however, a minimum data movement overhead of Ω(

d-1 k=0 (nz (W (k) ) + n k )(p t g + (p -p t )h
)) as all threads need to refer to the shared weights; replication of the weights on each socket would reduce this to Ω(

d-1 k=0 (nz (W (k) ) + n k )p t g).
This bound increases for each weight matrix whose accesses to, depending on the density pattern of X (k) , do not fit in cache. The storage overhead is Θ(p) since every thread maintains a CRS matrix.

C. Model-parallel inference

Model-parallel inference partitions the neural network, which in our case corresponds to a partitioning of the staircase matrix in Figure 1 mapping each individual weight to a unique thread. Let W (k) q be the matrix consisting of the set of weights of W (k) distributed to thread q, and R (k) = {r

(k) 0 , . . . , r (k) p-1 ; r (k)
S } be a partitioning of row indices of W (k) such that all weights on rows in R (ki) are in W (ki) (internal rows) while R (k S ) are rows shared by multiple parts of W (k) (external rows). Let Z (k) = {z refer to the connectivities of row i and column j, respectively. When the layer number k can be inferred from context we shall omit the superscripts from the r

(k) q , z (k)
q , and λ (k) {i,j} . The matrix W (k) q is stored locally by thread q, as are the corresponding input and output feature matrices X (k) {:,rq∪r S } and X (k+1) {:,zq∪z S } Since the column sets r (k) q , z (k) q are disjoint, each thread q stores a disjoint set of features. The parts corresponding to separator columns r

(k) S , z (k)
S , however, are replicated across all threads. Since a thread q may refer to remote data such as the separator features X (k) {sq,r S } r for 0 ≤ r < p, r = q, we retain the subscript r to emphasize which data is allocated to what thread.

The nonzeros of W (k) q can be logically subdivided into four parts from W . We assume the number of rows n is larger than or equal to p and employ twophase reduction to collaboratively compute p-1 r=0 X (k+1) {:,z S } r . This subdivides the total number of rows n into p parts and makes each thread responsible for reducing approximately n/p rows only; thread q hence computes X (k+1)

{sq,z S } = p-1 r=0 X (k+1) {sq,z S } r
and stores that result locally in X (k+1) {sq,z S } q . Algorithm 1 summarizes the resulting approach.

Algorithm 1 The model-parallel layer-k inference at thread q.

Input: Local input X (k) {:,rq } q , separator part X (k) {sq ,r S } q , and weight matrix W (k) {rq ∪r S ,zq ∪z S } q . Output: Local output X (k+1) {:,zq } q and separator X (k+1) {sq ,z S } q . 1: for r ← 0 to p -1 do 2: X (k+1) {sr ,zq } q ← f (X (k) {sr ,rq } q W (k) {rq ,zq } q + 3: X (k) {sr ,r S } r W (k) {r S ,zq } q ) 4: X (k+1) {sr ,z S } q ← X (k) {sr ,rq } q W (k) {rq ,z S } q + 5: X (k) {sr ,r S } r W (k) {r S ,z S } q Without ReLU 6: exec barrier Only for k = 0 7: X (k+1) {sq ,z S } q ← f ( p-1 r=0 X (k+1) {sq ,z S } r ) 8: exec barrier
There is no work overhead in flops as the summations at two-phase reduction correspond to summations of a sequential SpGEMM, now delayed due to having to reduce them over multiple threads. Ignoring the symbolic phase and assuming the W (k) q can be cached, memory movement overhead is solely due to reading the input separator parts X (k) {sr,r S } r at Lines 3 and 5 as well as both writing and reading of partial results at Lines 5 and 7, respectively. This yields a data movement overhead of

Θ     d-1 k=0 nz(X (k) {:,r S } ) + max q p-1 r=0, r =q nz(X (k+1) {sq,z S } r ) + n ph     .
Comparing to the movement overhead of the data-parallel method, Θ(p d-1 k=0 nz(W (k) )h), the model-parallel variant is preferable whenever the combined size of all X (k) {:,r S } and X (k+1) {:,z S } is smaller than that of all W (k) . The former combined size is proportional to the number of external columns, which is exactly the objective function (the cutnet metric) of the partitioning process. If the W (k) q do not fit in cache, an additional Θ(

d-1 k=0 max q [nz(W (k) q ) + |r q ∪ r S |]g) data movement occurs, resulting in an additional overhead of Θ(((p d-1 k=0 max q |r q | + |r S | -n k /p) + ε k nz(W (k) ))g),
where ε is the attained load imbalance of the partitioning.

The parallel storage overhead is bound by the maximum number of extra nonzeros stored in the partial results X (k+1) {:,z S } r plus the extra CRS indexing arrays for pd local weight matrices, p buffer matrices, and p 2 two-phase reduction buffers:

O max 0≤k<d p-1 q=0 nz(X (k+1) {:,z S } q ) +p(n + d k=1 |r s |) .
The synchronization overhead is Θ(dpL).

D. Tiling model-parallel inference

Tiling refers to a technique which increases data reuse by prematurely pausing a loop's iteration such that its output may be used by the subsequent iteration immediately and while they are still cached. Tiling techniques may incur overheads such as recomputation of factors needed to continue a paused iteration. We achieve a tiling in the model-parallel algorithm just presented by processing a batch of data elements through all layers of the network, before moving on to another batch. This allows the intermediate results X (k) , 0 < k < d, to remain in cache, which is especially effective when also the local parts of the neural network fit in cache.

To compute the tiling (model-parallel) inference, each thread q invokes Algorithm 2 to run a model-parallel inference for each of the t = n/b size batches of b size rows of X (k) . To reduce the barrier overhead of Θ(dpL), each thread q overlaps the computation of local results X (k+1) {:,zq} q with the barrier, implemented as weak point-to-point synchronizations, before processing the contributions of remote threads X S | by t. Hence, there exists a trade-off between choosing higher block sizes that tile for smaller caches versus how well underlying neural networks can be partitioned. Since there is no explicit barrier, synchronization overheads lie far below O(tdpL) in practice.

E. Hybrid tiling and deep inference

We may combine the data-parallel and tiling inference methods. Assume p = p 0 p 1 threads, then hybrid inference splits n into p 0 parts using a block distribution and processes each part using the tiling inference with p 1 threads. Thus, p 0 thread groups use the same p 1 -way partitioned weights. We select the minimum p 1 for which weights fit in the combined local caches. Higher p 1 would only be sensible if the separator size would decrease, which normally is not the case. Choosing p 0 ≥ p/p 1 then enables use of (almost) all available cores.

Deep neural networks where weights do not fit local cache sizes can still benefit from tiling. First, we cut the network's layers into successive blocks, and then apply our proposed tiling algorithm for each block. For example, if d = 120 and we create blocks of 5 layers, inference should perform 24 calls to the tiling algorithm. This requires streaming from main memory not only X (0) , but all of X (5k) , 0 ≤ k < 24. An automated method for selecting an appropriate cut could proceed greedily, 'growing' blocks of layers for which local parts remain cacheable and separator sizes remain below a Algorithm 2 Tiling model-parallel inference at thread q.

1: for m ← 0 to n/bsize -1 do

2:

Let bm = {mbsize , . . . , min{(m + 1)bsize , n} -1}.

3:

Let {s0, . . . , sp-1} be a block partitioning of bm into p parts.

4:

Let Cq be a bm × n k+1 in-cache matrix 5:

Let Dq be a bm × n k in-cache matrix resetAvail (C {:,r S } q ), resetAvail (D {sq ,r S } q ) 8:

for k = 0 to d do if checkAvail (D {:,r S } r ) then 13:

D {sq ,r S } q ← D {sq ,r S } q + D {sq ,r S } r 14:

done ← done + 1 15:

r ← (r + 1) mod p 16:

if r equals q then 17:

r = (q + 1) mod p 18:

D {sq ,r S } q ← f (D {sq ,r S } q ) Delayed ReLU 19: resetAvail (C {sq ,r S } q ), signalAvail (D {sq ,r S } q ) 20:
if k = d then Intra-and inter-thread data movement 21: r = q, done = 0 22:

while done < p do 23: if checkAvail (D {sr ,r S } r ) then 24: C {sr ,zq } q ← f (D {sr ,rq } q W (k) {rq ,zq } q + 25: D {sr ,r S } r W (k) {r S ,zq } q ) 26: C {sr ,z S } q ← D {sr ,rq } q W (k) {rq ,z S } q + 27: D {sr ,r S } r W (k) {r S ,z S } q Without ReLU 28: done ← done + 1 29:
r ← (r + 1) mod p 30:

resetAvail (D {:,r S } q ), signalAvail (C {:,r S } q )

31:

Swap C and D 32:

X (d)
{bm,zq } ← D {:,zq }q Write back to memory 33:

X (d)
{sq ,z S } ← D {sq ,z S }q certain threshold, using the hybrid method to deal with case where this method ends up with relatively low p 1 . More elaborate schemes may be envisioned; whatever these secondary details may be, a blocked tiling approach as presented here should lie at the core of a competitive sparse neural network inference method. Assuming an appropriate partitioning, a blocked tiling inference identifies two parameters: a block size for which the combined consecutive layers fit cache, and a tile size for which the intermediate results can be cached as well. We proceed with a series of experiments demonstrating this general approach works for the Graph Challenge dataset.

III. EXPERIMENTS

We use two machines for experiments: an Ivy Bridge node consisting of two sockets each equipped with ten cores, and a Cascade Lake node with two sockets each with 22 cores. The L2 cache size on the Ivy Bridge node is 256 kB while on the Cascade Lake it is 1 MB per core. Their L3 cache sizes are 2.5 MB and 1.25 MB per core, respectively. Both have 32 KB of L1 data cache size per core.

We use single-precision floats for values and integers for indices, following earlier challenge submissions [START_REF] Davis | Write quick, run fast: Sparse deep neural network in 20 minutes of development time in SuiteSparse:GraphBLAS[END_REF]. Our implementation builds on an internal C++ GraphBLAS code base used in sequential mode ensuring that threads allocate data using a local allocation policy enforced by the libnuma library. OpenMP is used for parallelization together with a custom ANSI C module that implements the (almost) synchronizationfree mechanism of our tiling method. Code compiles using GCC 9.2.0. Both machines run Linux with kernel version 3.10.0 on Ivy Bridge and 5.4.0 on Cascade Lake.

The performance of the tiling inference exceeds that of the data-parallel algorithm only if the feature matrix X (0) and weights W (k) do not fit cache; otherwise, the same beneficial cache effects of tiling naturally take place for the data-parallel algorithm. For the tiling methods, we take the maximum tile size for which all thread-local tiles of X (k) remain in L3 cache, or both the thread-local tiles of X (k) and W (k) q for all k remain in L3 cache if the combined network does not fit L3 cache.

The density of C and D are unknown and change as inference proceeds. We assume they are dense and correct the resulting storage requirements using the density of the input matrix X (0) for the density parameter α ≤ 1. This accounts for the density of X (k) assuming it does not grow with k. We ensure that b size > 0 is divisible by p 1 .

We experiment on the first five layers of the 120-layer HPEC networks: 1) to confirm the beneficial effects of caching the intermediate results, and 2) to confirm our understanding of the performance characteristics of the proposed tiling and hybrid tiling inference methods. As the X (k) fit in cache, we expect no benefit from the tiling method versus our data-parallel baseline implementation on the Ivy Bridge machine for both the 1k-neuron and 4k-neuron HPEC neural networks (where additionally the X (k) have decreasing densities as k increases, see Table I). Indeed, the results in Table III confirm slowdowns of the tiling method versus the data-parallel baseline, which grow with p due to the parallel overhead increasing with the separator size. We omit the results for the 4k-neuron network as the same effect is observed.

For the 16kand 64k-neuron networks, the X (0) do not fit in cache. Due to decreasing density, however, it may be that later X (k) could fit; this would be beneficial to data-parallel inference. Regardless, we certainly expect better performance from our tiling variant for the first five layers. For increasing p, since the parallel overhead grows with the number of external rows and columns, we expect this benefit to decrease. Table III confirms both expected behaviors.

Having done pure tiling for various p, we are able to select suitable p 1 equal to the best performing p. We use the remaining threads for the hybrid inference. Table IV shows the results for the 5-layer 16kand 64k-neuron network. These confirm that the hybrid tiling inference scales, maintaining similar speedups versus the data-parallel baseline as for the pure tiling method with p = p 1 threads. From Table IV, we additionally conclude that using hyperthreads benefits both the data-parallel and tiling methods, in the latter case presumably made possible by the almost synchronization-free method employed between (hyper)threads.

The fastest inference on the first five HPEC layers using the full Ivy Bridge machine are obtained by our hybrid tiling method, which are 41% faster versus the data-parallel method for the 16k-neuron network, and 50% faster on the 64k one.

Finally, we repeat both experiments for the Cascade Lake machine in Tables V and VI. For the hybrid tiling experiments, we took various combinations of p 0 p 1 some of which result in using less than the 22 available cores per socket due to divisibility. The results in both tables confirm the same expected behavior, and thus show that the method applies to different architectures. On Cascade Lake, the fastest inference is 101% and 94% faster compared to the data-parallel method for the 16kand 64k-neuron experiment, respectively, using all 22 cores with hyperthreads in a 22 • 2 configuration.

IV. CONCLUSIONS

We propose a hypergraph model of sparsely connected neural networks by applying the fine-grain model on the staircase matrix. This matrix defines a block structure of all weight matrices that accurately captures dependencies and communication between layers. By applying multi-constraint load balance criteria we use this model to partition the layers while minimizing overhead during model-parallel inference.

We then propose an efficient, hybrid model-and dataparallel sparse neural network inference method that performs tiling through layers to increase cache reuse, and exploits weak synchronization to hide load imbalance and inter-thread synchronization costs. The hybrid tiling algorithm offers best results for inference on the first five layers of the deep RadiX-Net NNs, provided that the input feature matrix is sufficiently large as confirmed using the 16kand 64k-neuron networks. We demonstrated the approach works across architectures.

Optimizations such as the use of algorithms for sparse matrix-dense matrix [START_REF] Acer | Improving performance of sparse matrix dense matrix multiplication on large-scale parallel systems[END_REF] multiplication may be used when the compressed density reaches 100% during inference. Another potential improvement relates to the use of the cutnet metric: it should be possible to achieve overheads proportional to the λ -1-metric instead, since similar bounds were achieved in earlier work on sparse matrix-vector multiplication [START_REF] Yzelman | High-level strategies for parallel sharedmemory sparse matrix-vector multiplication[END_REF]. Allowing permutations between layers, which can also be modeled via hypergraphs [START_REF] Uc | Partitioning sparse matrices for parallel preconditioned iterative methods[END_REF], may allow for further decreased separator size. Work is in progress to efficiently apply the proposed hybrid tiling method, potentially with the mentioned optimizations, on batches of input layers successively. This will confirm the applicability and benefits of the proposed approach on a broader set of neural network inference tasks.

  {rq∪r S ,zq∪z S } q , each with different data dependencies: (i) processing W (k) {rq,zq} q only depends on the locally computed dataX (k) {:,rq} q ; (ii) processing W (k) {r S ,zq} q requires collaboratively computed data X (k) {sq,r S } r for all r; (iii) processing W (k) {rq,z S } q requires collaborative computation of the output p-1 r=0 X (k+1) {:,z S } r; and (iv) processing W (k) {r S ,z S } q requires both parallel read access on X (k) {sq,r S } r z S } r

  (k) {sr,r S } r =q . The costs of tiling inference are similar to that of modelparallel, with overheads multiplied by t while substituting b size for n. This retains overheads proportional to n, and magnifies overheads proportional to |r (k) S | and |z (k)

THE

  THE SIZES (IN MB) OF THE INPUT FEATURE MATRIX X (0) AND THE HPEC NEURAL NETWORKS WITH DIFFERENT NUMBERS OF LAYERS AND NEURONS. COMPRESSED DENSITY AND THE DENSITY OF FEATURE MATRICES X (k) . THE LAST COLUMN REPORTS THE TOTAL MEMORY (IN MB) OCCUPIED BY THE THREE CRS MATRICES IN X

	Number of	Input	Number of layers d
	neurons	matrix	120		480	1920
	1k	48.86	30.47	121.88	487.51
	4k	191.11	121.88	487.50	1950.01
	16k	754.46	487.50 1950.00	7800.01
	64k 2992.42 1950.00 7800.00 31200.01
			Table I
	Layer Compressed			Memory
	k	density Density	(in MB)
	0		10.38	10.38	136.70
	1		31.08	28.33	149.15
	2		20.65	13.94	80.99
	4		34.28		6.86	45.98
	6		62.56		4.26	29.07
	8		97.23		3.45	23.63
	10		99.93		3.12	21.93
	15	100.00		3.02	21.46
			Table II

Table III THE

 III TILING INFERENCE RESULTS OVER THE FIRST 5 LAYERS OF 1K-, 16K-, AND 64K-NEURON NETWORKS ON IVY BRIDGE. THE SEPARATOR IS GIVEN IN PERCENTAGE OF THE TOTAL INPUT AND OUTPUT COLUMN SIZE. SPEEDUPS ARE RELATIVE TO THE DATA-PARALLEL BASELINE.

				1k-neuron					16k-neuron				64k-neuron
		Separator	Tile			Time (in s)	Separator Tile	Time (in s)	Separator	Tile	Time (in s)
	p	(in %)	size Baseline Tiling (speedup)	(in %)	size	Baseline Tiling (speedup)	(in %)	size	Baseline Tiling (speedup)
	2	11.79 11392		2.01	2.22 ( 0.91 )	0.80	158	33.62	19.86 ( 1.69 )	2.08	220	119.90	67.64 ( 1.77 )
	3	20.61	8067		1.32	1.63 ( 0.81 )	5.37	117	22.07	14.77 ( 1.49 )	2.41	210	77.99	46.05 ( 1.69 )
	4	16.25	6628		1.02	1.32 ( 0.77 )	4.89	116	16.92	10.35 ( 1.63 )	4.55	200	61.57	34.57 ( 1.78 )
	5	34.79	4575		0.81	1.55 ( 0.52 )	6.02	95	13.44	10.01 ( 1.34 )	1.69	200	47.22	27.65 ( 1.71 )
	6	26.67	4278		0.68	1.09 ( 0.62 )	4.24	114	11.60	7.82 ( 1.48 )	2.53	180	39.24	23.91 ( 1.64 )
	7	40.70	3304		0.58	1.15 ( 0.50 )	7.34	63	9.55	6.80 ( 1.40 )	7.79	140	33.53	24.75 ( 1.35 )
	8	21.25	3160		0.53	0.79 ( 0.67 )	7.97	72	9.78	5.87 ( 1.67 )	3.68	160	32.64	18.53 ( 1.76 )
	9	41.81	2601		0.46	0.99 ( 0.46 )	6.77	81	7.56	6.81 ( 1.11 )	3.26	180	26.28	17.58 ( 1.49 )
	10	42.71	2500		0.42	0.83 ( 0.51 )	8.27	90	7.00	5.37 ( 1.30 )	5.23	100	24.45	17.93 ( 1.36 )
	20	46.16	1340		0.22	0.53 ( 0.53 )	16.74	180	3.66	3.87 ( 0.95 )	6.17	200	12.58	10.06 ( 1.25 )
							Time (in s)
		Network	p p 0	p 1	Baseline Tiling (speedup)
		16k 20	5	4	3.66	2.20 ( 1.66 )
		16k 40	20	2	2.56	1.81 ( 1.41 )
		64k 20	4	5	12.58	7.62 ( 1.65 )
		64k 40	8	5	9.72	6.49 ( 1.50 )
					Table IV

THE

  HYBRID TILING INFERENCE RESULTS OVER THE FIRST 5 LAYERS OF 16K-, AND 64K-NEURON NETWORKS ON IVY BRIDGE.

						Time (in s)
	Network	p p 0	p 1	Baseline Tiling (speedup)
	16k 20	5	4	4.65	2.61 ( 1.78 )
	16k 22	11	2	4.49	2.23 ( 2.01 )
	16k 24	6	4	3.96	2.19 ( 1.81 )
	16k 44	22	2	2.36	1.18 ( 2.00 )
	16k 44	4	11	2.36	1.63 ( 1.45 )
	64k 20	5	4	16.99	9.06 ( 1.88 )
	64k 22	11	2	15.85	8.49 ( 1.87 )
	64k 40	5	8	8.91	4.75 ( 1.88 )
	64k 44	22	2	8.68	4.47 ( 1.94 )
	64k 44	4	11	8.68	5.10 ( 1.70 )
				Table V
	THE HYBRID TILING INFERENCE RESULTS OVER THE FIRST 5 LAYERS OF
	16K-,

AND 64K-NEURON NETWORKS ON CASCADE LAKE.

Table VI THE

 VI TILING INFERENCE RESULTS OVER THE FIRST 5 LAYERS OF 16K-, AND 64K-NEURON NETWORKS ON CASCADE LAKE.

			16k-neuron		64k-neuron
		Tile		Time (in s)	Tile		Time (in s)
	p	size Baseline Tiling (speedup)	size Baseline Tiling (speedup)
	2	158	43.53	23.05 ( 1.89 )	220	166.34	82.58 ( 2.01 )
	3	117	28.80	17.28 ( 1.67 )	210	103.39	56.03 ( 1.85 )
	4	116	21.84	12.27 ( 1.78 )	200	79.19	42.03 ( 1.88 )
	5	95	17.94	11.93 ( 1.50 )	200	62.33	33.29 ( 1.87 )
	6	114	15.10	9.33 ( 1.62 )	180	52.40	28.63 ( 1.83 )
	7	63	12.80	8.29 ( 1.54 )	140	45.48	31.42 ( 1.45 )
	8	72	12.17	7.32 ( 1.66 )	160	42.87	22.54 ( 1.90 )
	9	81	10.47	8.07 ( 1.30 )	180	34.88	21.04 ( 1.66 )
	10	90	9.16	6.52 ( 1.40 )	100	32.36	21.86 ( 1.48 )
	11	99	9.10	6.28 ( 1.45 )	110	32.24	19.90 ( 1.62 )
	16	144	6.73	5.26 ( 1.28 )	160	23.55	12.50 ( 1.88 )
	22	198	4.73	5.20 ( 0.91 )	220	15.85	13.41 ( 1.18 )