Lucio Compagno
email: lcompagno@dii.unict.it

Natalia Trapani
email: ntrapani@dii.unict.it

Reliability Driven Standardization of Mechanical Seals for Petrochemical Applications

Keywords: Mechanical seal, reliability data, standardization, maintenance, life table 1

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

Today's challenge in Oil & Gas (O&G) industry is to become more "sustainable", that is to produce in a more efficient, economical and ecological way. To obtain this result, production systems require quality products and reliable equipments, able to guarantee high availability and service efficacy and efficiency.

As it is known, a significant proportion of industrial emissions occurs through spurious leaks in process system (fugitive emissions), and particularly from flanges, valves and rotating shafts. Fugitive emissions from European refineries were quantified up to 10,000 tons of dangerous substances per year [START_REF]Sealing Technology -BAT guidance note[END_REF][START_REF] Onat | A review of fugitive emissions[END_REF]. The primary purpose of a sealing system is to avoid a fluid loss to protect the environment from pollution due to hazardous substances and mechanical seals are the most diffuse solution for rotating shaft sealing because of their high safety performances and low maintenance costs.

A careful selection of mechanical seals, which requires an analysis of fitness for service (i.e. treated fluid, operating conditions), regular inspection and maintenance are the minimum requirements to obtain high performances [START_REF] Huebner | Advancements In Mechanical Sealing -API 682 Fourth Edition[END_REF], i.e. high reliability and low level of VOC (Volatile Organic Compound) emissions [START_REF]Sealing Technology -BAT guidance note[END_REF]. The concept of Best Available Techniques (BAT) was introduced as a key principle in the European IPPC (integrated Pollution Prevention and Control) Directive [START_REF]EC of the European Parliament and of the Council of 15 January 2008 concerning integrated pollution prevention and control (Codified version)[END_REF]. In this Directive the word 'Techniques' includes both the technology used and the way in which the installation is designed, built, managed, maintained, operated and decommissioned according to a lifecycle management approach.

The selection criteria is often based on technical standards edited by API (American Petroleum Institute), such as API 682 [START_REF]API Standard 682[END_REF] and ISO 21049 [START_REF]2011(E): Pumps -Shaft Sealing Centrifugal and Rotary Pumps[END_REF], EPA (Environmental Protection Agency), such as Method 21 [START_REF]Determination of volatile organic compound leak[END_REF], or ESA (European Sealing Association) [START_REF]Sealing Technology -BAT guidance note[END_REF].

In order to satisfy customer needs, manufacturers have developed very large catalogues with a wide mix of products for each use and also customized solutions, and they usually offers to their clients also the maintenance service.

For manufacturers which want or have (i.e. in a global service contract) to guarantee high maintenance service performances (as it is known [START_REF]Sealing Technology -BAT guidance note[END_REF], some solutions generically have a lower Total Emission Level/Lifecycle cost ratio), this could result in high inventory costs but also in higher production costs and maintenance time. In literature there are a lot of papers that define methods to estimate assets reliability and useful life left, in order to support condition based maintenance and prognostics. Most interesting is the GAMM method proposed by Barbera et al. [START_REF] Barberá | A case study of GAMM (Graphical Analysis for Maintenance Management) applied to water pumps in a sewage treatment plant, Chile[END_REF], which uses a graphical approach based on a nonparametric estimation of reliability function using historical data; Abiri-Jaromi et al. [START_REF] Abiri-Jahromi | A Two-Stage Framework for Power Transformer Asset Maintenance Management -Part I: Models and Formulations[END_REF] use an operational research approach to optimize maintenance expected costs; other methods to assess remaining useful life use knowledgebased approach (expert and fuzzy), life expectancy models (stochastic and statistical), Artificial Neural Networks and Physical, showing that the research on how estimate reliability of in service components is still popular [START_REF] Si | Remaining useful life estimation -A review on the statistical data driven approaches[END_REF][START_REF] Sikorska | Prognostic modelling options for remaining useful life estimation by industry[END_REF] and a unique solution method is not available.

In order to give a more efficient maintenance service and to reduce inventory and production costs for a mechanical seal manufacturer, a statistical data analysis approach was defined to obtain a standardization of installed base driven by operational conditions and reliability performances.

Case study

The mechanical seal manufacturer has an installed base of 3,846 items in 8 plants of different clients. The seals can be classified into two macro categories (single or double) which generates 9 families and 21 different typologies (see Fig. 1).

Fig. 1. Classification of seals into families and typologies

The installed base operate in different operational conditions that are caused by the potential combinations of the factors specified in Tab. 1.

Table 1. Operational conditions for the installed base

Methodology

The methodology for the research is composed of three phases:

1. Installed mechanical seals census, to verify the correct classification according to Fig. 2 and the completeness of operational data according to Table 1; 2. Failure data analysis of seals typology, using Life Table technique [START_REF] Cox | Regression models and Life-Table[END_REF] for assessing a "MTTF Family" (Mean Time To Failure) inside the present database; 3. Identify, within Families with similar application field, the most reliable and verify if it can substitute the others. Having verified the inability to perform statistical analysis on the failures suffered by each item, due to the very low number of single item failure data that would make non-significant the results of the analysis, the Company was required to group seals in uniform classes. This consolidation has led to the definition of 35 families containing seals with similar structural characteristics and service conditions.

Seals

The operational conditions database was analyzed in order to verify the completeness of information. Unfortunately the database had complete information only for 1,968 items of the 3,846 installed items (51.2%, 26 families) with a significant reduction of the useful data. For the operational conditions a codification was defined based on operational data defined in Table 1, thus obtaining 21 codes for dual seals and 23 ones for single seals, e.g. DCRND3T1P1 that is a Dual (D), Cartridge (CR), process fluid not hydrocarbon (N), fluid specific mass γ > 0.8 kg/dm³ (D3), process temperature T 200 °C (T1), process pressure P 50 bar (P1).

Failure Data Analysis

The Company provided the records of the dates of failure for each item, surveyed in 64 months, stating the causes of failure identified through a Root Cause Analysis (RCA).

The reliability analysis is based on the evaluation of failure times of items that are part of each family using Life Table theory to exploit as statistically significant also the reliability data of seals that have never failed in the observation period.

The data provided by the Company have the following characteristics:

1. the installation period is known for each component, but not the exact date of installation, then there are doubts on the "service start date"; 4. seals are installed at different times and at each failure they were replaced with different seals (new or reconditioned); assuming that the reconditioning is able to return the seal "as new" it is as if, for any fault, a new component was installed; 5. dates of failure are note from 1 January 2006 to 30 April 2011 (64 months); 6. many components didn't fault in the period of observation.

Therefore, the data cannot be considered to be complete and adequate to run a "classical" statistical analysis. This led to the decision of executing an analysis using the Life Table methodology which allows to use data relating to seals that have never failed and data relating to the last period of service (from the last failure to April 30, 2011) by treating them as "censored data".

The Life Table is a table of extended frequency that, at discrete time intervals, counts the number of items still running at the beginning of the interval (Nr), of those failed in a period (Nf), and of those no longer controllable (e.g. lost to control or censored, Nc) to assess the elements at risk Ne, by the formula: Ne = Nr -Nc/2. Using this technique it is possible to reach reliable estimates of the four main functions for the study of the life of a component (survival function, probability den-sity function, failure rate, MTTF). The minimum size required for the analysis is 30 elements, including fault data and censored data.

In table 2 the life table of a dual seal family (DOTMS) is reported. The table shows all the variables that can be calculated (the calculated standard error often show very low values). The input data for the construction of the life table in this work are the times between repair and the service interval of items lost to control (censored data), according to the following assumptions:

1. installation date was considered to overlap the last failure date available in 2005, or the date of 01/01/2006 for all those items whose installation was before 2006; 2. if the item has no faults until 30 April 2011, it was assumed a service period overlapping the entire period of observation; 3. for items that have unknown installation period, the first failure date was considered as installation date.

A goodness-of-fit test for the PDF of each family was then executed by a statistical software application to determine whether it follows a known theoretical distribution (exponential, normal,...). Most of them are normally distributed, some of them are exponentially distributed, as shown in Fig. 2.

Fig. 2. Goodness of fit test of probability density function for two sets of data

In any case it is possible to calculate the MTTF using the formula (2):

n 1 k k k k Δt (t) f t MTTF ( 2
)
where t k is the mean value of the Δt k interval and f k (t) is the interval estimation of the PDF f(t). The MTTF value will be the base for standardization.

Standardization

Crossing the MTTF data with the seals operational conditions it is possible to establish which families are more reliable for each service (operational conditions as specified in Table 1 and paragraph 3.1), as shown in Fig. 3.

Fig. 3. Families MTTF vs construction and operational conditions

Discussion

Analysis of Figure 3 allows to establish that for the same range of operational conditions the manufacturer proposes from five to eight different seal families solutions but from the MTTF analysis it is clear that only a few of these can be considered Code reliable in the specific range, e.g. considering the code DCRID3T1P1 four of the installed families have a MTTF less than 800 days and four over 1200 days. The substitution of less reliable seals with the more reliable ones could generate benefits both for the manufacturer who could significantly reduce his catalogue and spare parts inventory with related costs, and for the customer who obtains reduced downtimes, reduced production costs, reduced potential environmental costs due to the better control of fugitive emissions, reduced seals maintenance and repair costs.

The cumulative results of standardization is synthesized in table 3. The overall results show that from 26 seal families it was suggested to reduce to only 16 families; the seal families replacement, that will be executed gradually (at the first maintenance intervention of each item), is referred to 646 item of 13 families which could be replaced by seal of only the three most reliable families specified in Table 3. In some cases, e.g. DOTA replacement of DSP or DSPAP generate a MTTF more than double, with obvious advantages for the reliability of the equipment.

Conclusions

Standardization effects are immediately visible both for the producer and customers. For the former in terms of stock costs reduction thanks to an increased interchangeability of spare parts which allows to reduce the size of the order quantity maintaining the same level of risk; for the customers in terms of increased reliability and availability of rotating shafts on which they are installed, with a direct effect on people and environment safety.

The manufacturer also experimented other less tangible benefits which consist of a greater maintenance specialization, a production mix reduced with benefits on maintenance lead time, customer service level, production lead time and ultimately lower costs of production.

Further studies on the long range effects of the suggested standardization on asset management costs will be probably available after a new census of the current installed base.

6.

Table 2 .

 2 Life table for DOTMS family

	Int. Start	Nr. Items	Nr.	Nr.	Nr.
	[days]	running	censored	exposed	Dead
	0	90	9	85.5	21
	250	60	3	58.5	13
	500	44	3	42.5	9
	750	32	0	32.0	4
	1000	28	4	26.0	6
	1250	18	1	17.5	3
	1500	14	1	13.5	2
	1750	11	11	5.5	0
	Int. Start	Survived	Cum. Survived	PDF	Failure
	[days]	fraction	fraction		rate
	0	0.75439	1.00000	0.00098	0.0011
	250	0.77778	0.75439	0.00067	0.0010
	500	0.78824	0.58674	0.00050	0.0009
	750	0.87500	0.46249	0.00023	0.0005
	1000	0.76923	0.40468	0.00037	0.0010
	1250	0.82857	0.31129	0.00021	0.0008
	1500	0.85185	0.25793	0.00015	0.0006
	1750	0.90909	0.21972	-	-

Table 3 .

 3 Seal families replacement as effect of standardization

	Old Family MTTF old [d]	ITEM Nr. Total Item New Family MTTF new [d]
	DOTA	1318	161
	DOTD	1310	97
	DSTL	1246	9
	DODAP	1173	61	450	DOAP	1512
	DSTD	741	26
	DOP	689	48
	DST	645	48
	SD	1210	11
	SOGP	1068	8	115	SSF	1498
	SSL	986	96
	DSP DSPAP	575 560	56 25	81	DOTA	1318
	Total suggested replacement	646	646


Philippe Thomas
email: philippe.thomas@univ-lorraine.fr

Marie-Christine Suhner
email: marie-christine.suhner@univ-lorraine.fr

André Thomas
email: andre.thomas@univ-lorraine.fr

CART for supply chain simulation models reduction

Keywords: Decision tree, CART, reduced model, simulation, supply chain

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Planning or scheduling scenario evaluation is an important issue for internal/external Supply Chains (SC) control. Simulation is useful to perform this evaluation and allows to highlight evolution of the resources states, work in process and queues. This information allows to perform a "Predictive scheduling" [START_REF] Lopez | Ordonnancement de la production[END_REF] which concerns MPS which is initially established with the Manufacturing Planning and Control System. This first approach must not be confused with "Reactive scheduling" which gives new MPS solution, established after significant events occur on the shop-floor. The real time systems performing production reporting lead to obtain information into the management system very quickly [START_REF] Khouja | An aggregate production planning framework for the evaluation of volume flexibility[END_REF]. However, it is difficult to use this huge of information to take a decision [START_REF] Pritsker | Simulation for planning and scheduling[END_REF][START_REF] Roder | Visibility is the key to scheduling success[END_REF]. At this level of planning, load/capacity balancing is obtained via the "management of critical resource capacity" function or Rough-Cut Capacity Planning which essentially concerns bottlenecks [START_REF] Vollmann | Manufacturing, Planning and Systems Control[END_REF]. Goldratt and Cox [4] propose to manage all the supply chain by bottlenecks control and call it Theory of Constraints (TOC). For this purpose, the use of dynamic discrete event simulation of material flow is helpful [START_REF] Thomas | Reducing simulation models for scheduling manufacturing facilities[END_REF]. Simulation models of real industrial cases are very com-

Code de champ modifié

plex and can lead to problems of scale [START_REF] Page | An aggregate production planning framework for the evaluation of volume flexibility[END_REF] that is why it is interesting to use simplest (reduced/aggregated) models of simulation [START_REF] Brooks | Simplification in the simulation of manufacturing systems[END_REF][START_REF] Chwif | Discret event simulation model reduction: A causal approach[END_REF][START_REF] Ward | Argument for constructively simple models[END_REF]. To design these reduced model, different approaches may be used as using continuous flow model to approximate discrete manufacturing environments [START_REF] Ho | Performance evaluation and perturbation analysis of discrete event dynamics systems[END_REF]15], or using metamodels (linear regression, splines, Kriging…) in order to perform simulation model [START_REF] Kleijnen | A methodology for fitting and validating metamodels in simulation[END_REF]. Neural networks are also been used to perform this task [START_REF] Thomas | Multilayer Perceptron for Simulation Models Reduction: Application to a Sawmill Workshop[END_REF][START_REF] Thomas | A neural network for the reduction of a Product Driven System emulation model[END_REF]. Another tool, regression tree, is able to find a suitable model useful for the model reduction. Classification and regression trees are machine learning methods able to fit a model to data. Models are obtained by partitioning data space and thanks to the association of a simple prediction model to each subspace [START_REF] Breiman | Classification and regression trees[END_REF]. The main goal of this paper is to evaluate the ability of regression trees to be used in order to reduce simulation model. In the second part, the proposed approach of reduction model and the regression trees are presented. The third part will be devoted to the presentation of an industrial application which is a sawmill flow shop case. In order to evaluate the pertinence of the proposition, the reduced simulation model of sawmill internal SC case is presented and, before to conclude, is compared to the same work performed by using neural network in the last part.

2

The reduced model

The algorithm

Model complexity is relative to the number of elements, connections and model calculations [START_REF] Zeigler | Theory of modelling and simulation[END_REF]. Different approaches have been proposed to reduce model. To see a brief overview, see [START_REF] Thomas | A neural network for the reduction of a Product Driven System emulation model[END_REF]. The proposed reduction algorithm is a modification of those presented by Thomas and Thomas [START_REF] Thomas | Multilayer Perceptron for Simulation Models Reduction: Application to a Sawmill Workshop[END_REF]. The proposed approach is based on the association of discrete event models and continuous models (regression tree) in order to design a simulation model. Our objective is to maximize the bottleneck utilization rate and, at the same time, simplify simulation model design for modelers. Its main steps are recalled and explained below:

1. Identify the structural bottleneck (work center (WC) which for several years has been mainly constrained in capacity). 2. Identify the conjunctural bottleneck for the bundle of Manufacturing Orders (MO) of the MPS under consideration. 3. Among the WC not listed in 1 and 2, identify the one (synchronisation WC) satisfying these two conditions: (a) present at least in one of the MO using a bottleneck, (b) widely used considering the whole MO. 4. If all MO have been considered go to 5 if not go to 3.

Use regression tree for modelling the intervals between WC that have been found

during preceding steps (figure 1).

Work Centers (WC) remaining in the model are either conjunctural or structural bottlenecks or WC which are vital to the synchronization of the MO. Other WC are incorporated in "aggregated blocks" upstream or downstream of the bottlenecks [START_REF] Thomas | A neural network for the reduction of a Product Driven System emulation model[END_REF].

Regression tree

Decision trees are sequential models, which logically combine a sequence of simple tests. A regression tree looks like a classification tree, except that the output variable is a continuous one and a regression model is fitted to each node to give the predicted value of the output [START_REF] Loh | Classification and regression trees[END_REF]. The first regression tree (AID) has been proposed by Morgan and Sonquist [START_REF] Morgan | Problems in the analysis of survey data, and a proposal[END_REF]. The principle of this regression tree is to begin to the root node, and then, for each node, find the test on the input patterns which minimizes the sum of the node impurities in its two child nodes. Classification and Regression Tree (CART) algorithm works on the same principle [START_REF] Breiman | Classification and regression trees[END_REF]. The main difference between AID and CART relies with the pruning and estimation process. CART is based on the generalization of the binomial variance called Gini index. The growing procedure of the tree is recursive and the process is stopped when [START_REF] Lewis | An introduction to classification and regression tree (CART) analysis[END_REF]:

─ there is only one observation in each of the child nodes, ─ all observations within each child node have the identical distribution of predictor variables, making splitting impossible, ─ an external limit on the number of levels in the tree or on impurity decrease level threshold is reached.

The resulting tree presents generally overfiting problem and needs a pruning phase, based on the impurity term used during the growing phase (Gini) associated to a penalty term corresponding to the number of terminal nodes. The final tree selection is performed by using a cross-validation procedure [START_REF] Breiman | Classification and regression trees[END_REF].

Overview of the sawmill

The proposed approach is applied to build a simulation model of a sawmill workshop. In this actual case, managers need a tool to help them in their weekly decision-making Master Production schedule (MPS) process. They want:

─ to evaluate the effectiveness of its MPS, ─ to maximize its load rate, and so, its global productivity, ─ to explain some unexplained congestion phenomena of the trimmer WC.

A first work [START_REF] Thomas | Reducing simulation models for scheduling manufacturing facilities[END_REF] with a complete model has allowed to represent the congestion phenomena and to use this representation in order to improve the load rate. This model, has showed that a bottleneck load rate too high (higher than 60%) degrades productivity of bottleneck, and so, sawmill productivity. The difficulty is that bottleneck is the last sawmill WC but all influent factors on bottleneck productivity depend of the first WC. However, this complete model is unusable on a real case for the MPS dynamic evaluation because of time needed to modify it. Sawmill can be described from a process point of view. It can be represented by two linear parallel flows for main and secondary products. This fact, associated to variation of log dimensions lead process to be non-linear. Therefore, the physical industrial production system can be divided into three main parts. To explain the process, the trip of a log from its admission into the process to its exit in planks form will be described. The first process part is the Canter line presented figure 2. Dashed arrows indicate the products flow. The considered log is taken into the process by using conveyors RQM1, RQM2, and RQM3. According to its characteristics (determined by scanner MS), the log is driven to conveyors RQM4 or RQM5, which are used as input inventory for the Canter line. After that, logs go on the first canter's machine and later on the CSMK saw, which transforms logs into square-shaped parallelepipeds. This step gives the two first sides of the parallelepipeds (the added value on the log) and produces two planks (secondary products), which are taken out of the Canter line by the BT4 and BT5 conveyors. The log is then driven on the RQM6 conveyors, rotated 90°, and stored in RQM7 awaiting its second passage on the CSMK saw. After the second passage, the squared is completed, and two other secondary products are taken out of the Canter line by BT4 and BT5 toward the second part of the process, the Kockums line. The squared log is cut on the MKV saw into three planks (called main products). These main products are driven to the third part of the process, the trimmer line. Figure 3 shows the second part of the process, where the main machine is the Kockums saw. Only secondary products are driven on this part. As previously said, the secondary products are taken in the line by the BT4 and BT5 conveyors. They are cut by the QM11 saw, after which they reach the Kockums saw, which optimizes the planks according to the products needed. The alignment table is used as the input inventory of the Kockums saw. The secondary products are finally sent to the third part of the process by the exit conveyor.

The third part of the process is the trimmer line, which is presented in figure 4. This line performs the final operation of cross cutting. This operation consists in cutting up products to length. The inputs of the line are from collectors 1 and 2, which collect the secondary and main products from Kockums and Canter lines respectively. Saw 1 is used to perform default bleeding and Saw 2 cuts up products to length. A previous work [START_REF] Thomas | Reducing simulation models for scheduling manufacturing facilities[END_REF] has shown that this last machine, the trimmer saw, is the bottleneck of the entire process, and, as said previously, the productivity of the trimmer depends to the decisions taken on the canter work center. So, the impact of bad decisions are seen too much late to be corrected. So in order to evaluate decisions, managers need a simulation tool.

The simulation models

As said previously, the bottleneck of this line is the trimmer. Consequently, modeling the function of the inventories RQM4, RQM5, and RQM7, and of the Canter and Kockums lines is unnecessary and gave no direct and useful information for the evaluation of the MPS. In fact, only the arrival times of the products in the trimmer queue are useful for simulating the load of this bottleneck. So, in the complete model presented figure 5, all the parts surrounded with dashed line is unnecessary and may be replaced by a continuous model, here a regression tree.

According to the modeling process recalled in part 2, the specific sawmill model could be designed. To build this model, we need to identify the input variables. Thomas et al. [START_REF] Thomas | A neural network for the reduction of a Product Driven System emulation model[END_REF] collected the available input data which can be classified into three categories: Data related to the products (here the logs), data related to the process and data related to the bill of material or routing (here the cutting plan). The data related to the products are mainly dimensional ones as length (lg) and three values for timber diameters (diaPB, diaGB, and diaMOY). The data related to the process are the process variables collected at the time of log arrival. In particular, we require the input stock and the utilization rate of the bottleneck, here the trimmer (Q_trim, and U_trim, respectively). The number of logs present in the process between the inputs and the exit of the Canter line (Q_RQM) and in the input inventories (Q_RQM4, Q_RQM5 and Q_RQM7) are needed. The data related to the routing correspond here to the information related to the cutting plan of the logs which must be cut into main and secondary products. This information is given by the type of part (T_piece). Consequently, the model input variables are: Lg, diaGB, diaMOY, diaPB, T_piece, Q_trim, U_trim, Q_RQM, Q_RQM4, Q_RQM5 and Q_RQM7. In our application, 12775 products are simulated with the complete model. These data are used to fit the behavior of the reduced model to the complete one which serve as reference model.

The next step is to identify the output variable. Our objective is to estimate the delay ( T) corresponding to the duration of the throughput time for the 12775 products.

T is measured between the process input time and the trimmer queue input time. For regression problems, the goal is to reduce the distance between the model and the data. Hence, the RMSE criterion is suitable, and is used in this work: where N denotes the number of data, y n is the n th actual data point, and ˆn y is the pre- dicted value. Figure 6 presents the evolution of the RMSE in function of the number of nodes in the tree. This figure shows that classically, even if RMSE always decreases on the learning data set, it is not the same on the validation data set. It begins by decreasing before to increase. This fact illustrates the overfitting problem and the selection of the best model must be performed by using a cross validation strategy. In our case, the selected model comprises 296 nodes. This best model must be compared with those obtained in preceding work which use neural network model [START_REF] Thomas | A neural network for the reduction of a Product Driven System emulation model[END_REF]. These results are presented table 1. They show that regression tree outperforms neural model and allows to improve results of 20%.

Conclusion

The use of regression tree in order to build a reduced simulation model is investigated here. The regression tree model is used to model the functioning of a part of the process that is not constrained in capacity. This approach has been applied to the modeling of a sawmill workshop. The results show that the complete and reduced models gave similar results even if the log arrival rule is changed. Moreover, the results are compared with results obtained with a neural network model and the comparison shows that regression tree outperforms neural network on this application. This means that it seems efficient to use a regression tree to model a part of a process instead of constructing the complete model. Assuming that the construction of a regression tree is a quasi-automated task, in which the modeler only collects and selects the input data set. It is faster and easier to design this kind of reduced model. This approach allows the modeler to focus on the management of bottlenecks. Our intentions for future works are the validation of this approach on different applications, particularly

Fig. 11 .

 11 Fig. 11. Reduction model algorithm

Fig. 22 .

 22 Fig. 22. The Canter line

Fig. 44 .

 44 Fig. 44. The trimmer line

Fig. 55 .

 55 Fig. 55. Complete model

Fig. 66 .

 66 Fig. 66. Evolution of RMSE in function of nodes numberidentification data set in blue dashedvalidation data set in red

Table 11 .

 11 Comparison of RMSE for Neural network and regression tree models

  Police :Gras, Anglais (États Unis), Ne pas vérifier l'orthographe ou la grammaire on several external supply chains, such that at least one particular enterprise belongs to different supply chains.

	RMSE
	350
	300
	250
	200
	150	0	500	1000	1500	2000	2500	Nodes number
				RMSE Identification	Validation
				Neural model	408.45	413.93
				Tree model	263.47	324.10

Mis en forme : Police :Gras, Ne pas vérifier l'orthographe ou la grammaire Mis en forme :

Mis en forme : Police :Gras, Ne pas vérifier l'orthographe ou la grammaire

Mis en forme : Police :Gras, Ne pas vérifier l'orthographe ou la grammaire Mis en forme : Police :Gras, Ne pas vérifier l'orthographe ou la grammaire


Julien Mairal
email: julien.mairal@inria.fr

End-to-End Kernel Learning with Supervised Convolutional Kernel Networks

In this paper, we introduce a new image representation based on a multilayer kernel machine. Unlike traditional kernel methods where data representation is decoupled from the prediction task, we learn how to shape the kernel with supervision. We proceed by first proposing improvements of the recently-introduced convolutional kernel networks (CKNs) in the context of unsupervised learning; then, we derive backpropagation rules to take advantage of labeled training data. The resulting model is a new type of convolutional neural network, where optimizing the filters at each layer is equivalent to learning a linear subspace in a reproducing kernel Hilbert space (RKHS). We show that our method achieves reasonably competitive performance for image classification on some standard "deep learning" datasets such as CIFAR-10 and SVHN, and also for image super-resolution, demonstrating the applicability of our approach to a large variety of image-related tasks.

Introduction

In the past years, deep neural networks such as convolutional or recurrent ones have become highly popular for solving various prediction problems, notably in computer vision and natural language processing. Conceptually close to approaches that were developed several decades ago (see, [START_REF] Le Cun | Efficient backprop[END_REF]), they greatly benefit from the large amounts of labeled data that have been made available recently, allowing to learn huge numbers of model parameters without worrying too much about overfitting. Among other reasons explaining their success, the engineering effort of the deep learning community and various methodological improvements have made it possible to learn in a day on a GPU complex models that would have required weeks of computations on a traditional CPU (see, e.g., [START_REF] He | Deep residual learning for image recognition[END_REF][START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF][START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF]).

Before the resurgence of neural networks, non-parametric models based on positive definite kernels were one of the most dominant topics in machine learning [START_REF] Schölkopf | Learning with kernels: support vector machines, regularization, optimization, and beyond[END_REF]. These approaches are still widely used today because of several attractive features. Kernel methods are indeed versatile; as long as a positive definite kernel is specified for the type of data considered-e.g., vectors, sequences, graphs, or sets-a large class of machine learning algorithms originally defined for linear models may be used. This family include supervised formulations such as support vector machines and unsupervised ones such as principal or canonical component analysis, or K-means and spectral clustering. The problem of data representation is thus decoupled from that of learning theory and algorithms. Kernel methods also admit natural mechanisms to control the learning capacity and reduce overfitting [START_REF] Schölkopf | Learning with kernels: support vector machines, regularization, optimization, and beyond[END_REF].

On the other hand, traditional kernel methods suffer from several drawbacks. The first one is their computational complexity, which grows quadratically with the sample size due to the computation of the Gram matrix. Fortunately, significant progress has been achieved to solve the scalability issue, either by exploiting low-rank approximations of the kernel matrix [START_REF] Williams | Using the Nyström method to speed up kernel machines[END_REF][START_REF] Zhang | Improved Nyström low-rank approximation and error analysis[END_REF], or with random sampling techniques for shift-invariant kernels [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF]. The second disadvantage is more critical; by decoupling learning and data representation, kernel methods seem by nature incompatible with end-to-end learning-that is, the representation of data adapted to the task at hand, which is the cornerstone of deep neural networks and one of the main reason of their success. The main objective of this paper is precisely to tackle this issue in the context of image modeling. Specifically, our approach is based on convolutional kernel networks, which have been recently introduced in [START_REF] Mairal | Convolutional kernel networks[END_REF]. Similar to hierarchical kernel descriptors [START_REF] Bo | Object recognition with hierarchical kernel descriptors[END_REF], local image neighborhoods are mapped to points in a reproducing kernel Hilbert space via the kernel trick. Then, hierarchical representations are built via kernel compositions, producing a sequence of "feature maps" akin to convolutional neural networks, but of infinite dimension. To make the image model computationally tractable, convolutional kernel networks provide an approximation scheme that can be interpreted as a particular type of convolutional neural network learned without supervision.

To perform end-to-end learning given labeled data, we use a simple but effective principle consisting of learning discriminative subspaces in RKHSs, where we project data. We implement this idea in the context of convolutional kernel networks, where linear subspaces, one per layer, are jointly optimized by minimizing a supervised loss function. The formulation turns out to be a new type of convolutional neural network with a non-standard parametrization. The network also admits simple principles to learn without supervision: learning the subspaces may be indeed achieved efficiently with classical kernel approximation techniques [START_REF] Williams | Using the Nyström method to speed up kernel machines[END_REF][START_REF] Zhang | Improved Nyström low-rank approximation and error analysis[END_REF].

To demonstrate the effectiveness of our approach in various contexts, we consider image classification benchmarks such as CIFAR-10 [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] and SVHN [START_REF] Netzer | Reading digits in natural images with unsupervised feature learning[END_REF], which are often used to evaluate deep neural networks; then, we adapt our model to perform image super-resolution, which is a challenging inverse problem. On the SVHN and CIFAR-10 datasets, we obtain a competitive accuracy, with about 2% and 10% error rates, respectively, without model averaging or data augmentation. For image up-scaling, we outperform recent approaches based on classical convolutional neural networks [START_REF] Dong | Learning a deep convolutional network for image super-resolution[END_REF][START_REF] Dong | Image super-resolution using deep convolutional networks[END_REF].

We believe that these results are highly promising. Our image model achieves competitive performance in two different contexts, paving the way to many other applications. Moreover, our results are also subject to improvements. In particular, we did not use GPUs yet, which has limited our ability to exhaustively explore model hyper-parameters and evaluate the accuracy of large networks. We also did not investigate classical regularization/optimization techniques such as Dropout [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], batch normalization [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF], or recent advances allowing to train very deep networks [START_REF] He | Deep residual learning for image recognition[END_REF][START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF]. To gain more scalability and start exploring these directions, we are currently working on a GPU implementation, which we plan to publicly release along with our current CPU implementation.

Related Deep and Shallow Kernel Machines. One of our goals is to make a bridge between kernel methods and deep networks, and ideally reach the best of both worlds. Given the potentially attractive features of such a combination, several attempts have been made in the past to unify these two schools of thought. A first proof of concept was introduced in [START_REF] Cho | Kernel methods for deep learning[END_REF] with the arc-cosine kernel, which admits an integral representation that can be interpreted as a one-layer neural network with random weights and infinite number of rectified linear units. Besides, a multilayer kernel may be obtained by kernel compositions [START_REF] Cho | Kernel methods for deep learning[END_REF]. Then, hierarchical kernel descriptors [START_REF] Bo | Object recognition with hierarchical kernel descriptors[END_REF] and convolutional kernel networks [START_REF] Mairal | Convolutional kernel networks[END_REF] extend a similar idea in the context of images leading to unsupervised representations [START_REF] Mairal | Convolutional kernel networks[END_REF].

Multiple kernel learning [START_REF] Sonnenburg | Large scale multiple kernel learning[END_REF] is also related to our work since is it is a notable attempt to introduce supervision in the kernel design. It provides techniques to select a combination of kernels from a predefined collection, and typically requires to have already "good" kernels in the collection to perform well. More related to our work, the backpropagation algorithm for the Fisher kernel introduced in [START_REF] Sydorov | Deep Fisher kernels -end to end learning of the Fisher kernel GMM parameters[END_REF] learns the parameters of a Gaussian mixture model with supervision. In comparison, our approach does not require a probabilistic model and learns parameters at several layers. Finally, we note that a concurrent effort to ours is conducted in the Bayesian community with deep Gaussian processes [START_REF] Damianou | Deep Gaussian processes[END_REF], complementing the Frequentist approach that we follow in our paper.

Learning Hierarchies of Subspaces with Convolutional Kernel Networks

In this section, we present the principles of convolutional kernel networks and a few generalizations and improvements of the original approach of [START_REF] Mairal | Convolutional kernel networks[END_REF]. Essentially, the model builds upon four ideas that are detailed below and that are illustrated in Figure 1 for a model with a single layer.

Idea 1: use the kernel trick to represent local image neighborhoods in a RKHS. Given a set X , a positive definite kernel K : X × X → R implicitly defines a Hilbert space H, called reproducing kernel Hilbert space (RKHS), along with a mapping ϕ : X → H. This embedding is such that the kernel value K(x, x ) corresponds to the inner product ϕ(x), ϕ(x ) H . Called "kernel trick", this approach can be used to obtain nonlinear representations of local image patches [START_REF] Bo | Object recognition with hierarchical kernel descriptors[END_REF][START_REF] Mairal | Convolutional kernel networks[END_REF].

More precisely, consider an image I 0 : Ω 0 → R p0 , where p 0 is the number of channels, e.g., p 0 = 3 for RGB, and Ω 0 ⊂ [0, 1]2 is a set of pixel coordinates, typically a two-dimensional grid. Given two image patches x, x of size e 0 × e 0 , represented as vectors in R p0e 2 0 , we define a kernel K 1 as

K 1 (x, x ) = x x κ 1 x x , x x if x, x = 0 and 0 otherwise, (1)
where . and ., . denote the usual Euclidean norm and inner-product, respectively, and κ 1 ( ., . ) is a dot-product kernel on the sphere. Specifically, κ 1 should be smooth and its Taylor expansion have non-negative coefficients to ensure positive definiteness [START_REF] Schölkopf | Learning with kernels: support vector machines, regularization, optimization, and beyond[END_REF]. For example, the arc-cosine [START_REF] Cho | Kernel methods for deep learning[END_REF] or the Gaussian (RBF) kernels may be used: given two vectors y, y with unit 2 -norm, choose for instance

κ 1 ( y, y ) = e α1( y,y -1) = e -α 1 2 y-y 2 2 . (2)
Then, we have implicitly defined the RKHS H 1 associated to K 1 and a mapping ϕ 1 : R p0e 2 0 → H 1 .

Idea 2: project onto a finite-dimensional subspace of the RKHS with convolution layers.

The representation of patches in a RKHS requires finite-dimensional approximations to be computationally manageable. The original model of [START_REF] Mairal | Convolutional kernel networks[END_REF] does that by exploiting an integral form of the RBF kernel. Specifically, given two patches x and x , convolutional kernel networks provide two vectors ψ 1 (x), ψ 1 (x ) in R p1 such that the kernel value ϕ 1 (x), ϕ 1 (x ) H1 is close to the Euclidean inner product ψ 1 (x), ψ 1 (x ) . After applying this transformation to all overlapping patches of the input image I 0 , a spatial map M 1 : Ω 0 → R p1 may be obtained such that for all z in Ω 0 , M 1 (z) = ψ 1 (x z ), where x z is the e 0 × e 0 patch from I 0 centered at pixel location z. 2 With the approximation scheme of [START_REF] Mairal | Convolutional kernel networks[END_REF], M 1 can be interpreted as the output feature map of a one-layer convolutional neural network.

A conceptual drawback of [START_REF] Mairal | Convolutional kernel networks[END_REF] is that data points ϕ 1 (x 1 ), ϕ 1 (x 2 ), . . . are approximated by vectors that do not live in the RKHS H 1 . This issue can be solved by using variants of the Nyström method [START_REF] Williams | Using the Nyström method to speed up kernel machines[END_REF], which consists of projecting data onto a subspace of H 1 with finite dimension p 1 .

For this task, we have adapted the approach of [START_REF] Zhang | Improved Nyström low-rank approximation and error analysis[END_REF]: we build a database of n patches x 1 , . . . , x n randomly extracted from various images and normalized to have unit 2 -norm, and perform a spherical K-means algorithm to obtain p 1 centroids z 1 , . . . , z p1 with unit 2 -norm. Then, a new patch x is approximated by its projection onto the p 1 -dimensional subspace F 1 = Span(ϕ(z 1 ), . . . , ϕ(z p1 )).

The projection of ϕ 1 (x) onto F 1 admits a natural parametrization ψ 1 (x) in R p1 . The explicit formula is classical (see [START_REF] Williams | Using the Nyström method to speed up kernel machines[END_REF][START_REF] Zhang | Improved Nyström low-rank approximation and error analysis[END_REF] and Appendix A), leading to

ψ 1 (x) := x κ 1 (Z Z) -1/2 κ 1 Z x x if x = 0 and 0 otherwise, (3)
where we have introduced the matrix Z = [z 1 , . . . , z p1 ], and, by an abuse of notation, the function κ 1 is applied pointwise to its arguments. Then, the spatial map M 1 : Ω 0 → R p1 introduced above can be obtained by (i) computing the quantities Z x for all patches x of the image I (spatial convolution after mirroring the filters z j ); (ii) contrast-normalization involving the norm x ; (iii) applying the pointwise non-linear function κ 1 ; (iv) applying the linear transform κ 1 (Z Z) -1/2 at every pixel location (which may be seen as 1×1 spatial convolution); (v) multiplying by the norm x making ψ 1 homogeneous. In other words, we obtain a particular convolutional neural network, with non-standard parametrization. Note that learning requires only performing a K-means algorithm and computing the inverse square-root matrix κ 1 (Z Z) -1/2 ; therefore, the training procedure is very fast.

Then, it is worth noting that the encoding function ψ 1 with kernel (2) is reminiscent of radial basis function networks (RBFNs) [START_REF] Broomhead | Radial basis functions, multi-variable functional interpolation and adaptive networks[END_REF], whose hidden layer resembles (3) without the matrix κ 1 (Z Z) -1/2 and with no normalization. The difference between RBFNs and our model is nevertheless significant. The RKHS mapping, which is absent from RBFNs, is indeed a key to the multilayer construction that will be presented shortly: a network layer takes points from the RKHS's previous layer as input and use the corresponding RKHS inner-product. To the best of our knowledge, there is no similar multilayer and/or convolutional construction in the radial basis function network literature.

I 0 x x kernel trick projection on F 1 M 1 ψ 1 (x) ψ 1 (x ) I 1 linear pooling Hilbert space H 1 F 1 ϕ 1 (x) ϕ 1 (x )
Figure 1: Our variant of convolutional kernel networks, illustrated between layers 0 and 1. Local patches (receptive fields) are mapped to the RKHS H 1 via the kernel trick and then projected to the finite-dimensional subspace F 1 = Span(ϕ(z 1 ), . . . , ϕ(z p1 )). The small blue crosses on the right represent the points ϕ(z 1 ), . . . , ϕ(z p1 ). With no supervision, optimizing F 1 consists of minimizing projection residuals. With supervision, the subspace is optimized via back-propagation. Going from layer k to layer k + 1 is achieved by stacking the model described here and shifting indices.

Idea 3: linear pooling in F 1 is equivalent to linear pooling on the finite-dimensional map M 1 .

The previous steps transform an image

I 0 : Ω 0 → R p0 into a map M 1 : Ω 0 → R p1
, where each vector M 1 (z) in R p1 encodes a point in F 1 representing information of a local image neighborhood centered at location z. Then, convolutional kernel networks involve a pooling step to gain invariance to small shifts, leading to another finite-dimensional map I 1 : Ω 1 → R p1 with smaller resolution:

I 1 (z) = z ∈Ω0 M 1 (z )e -β1 z -z 2 2 . ( 4
)
The Gaussian weights act as an anti-aliasing filter for downsampling the map M 1 and β 1 is set according to the desired subsampling factor (see [START_REF] Mairal | Convolutional kernel networks[END_REF]), which does not need to be integer. Then, every point I 1 (z) in R p1 may be interpreted as a linear combination of points in F 1 , which is itself in F 1 since F 1 is a linear subspace. Note that the linear pooling step was originally motivated in [START_REF] Mairal | Convolutional kernel networks[END_REF] as an approximation scheme for a match kernel, but this point of view is not critically important here.

Idea 4: build a multilayer image representation by stacking and composing kernels. By following the first three principles described above, the input image

I 0 : Ω 0 → R p0 is transformed into another one I 1 : Ω 1 → R p1 .
It is then straightforward to apply again the same procedure to obtain another map

I 2 : Ω 2 → R p2 , then I 3 : Ω 3 → R p3 , etc.
By going up in the hierarchy, the vectors I k (z) in R p k represent larger and larger image neighborhoods (aka. receptive fields) with more invariance gained by the pooling layers, akin to classical convolutional neural networks.

The multilayer scheme produces a sequence of maps (I k ) k≥0 , where each vector I k (z) encodes a point-say f k (z)-in the linear subspace F k of H k . Thus, we implicitly represent an image at layer k as a spatial map

f k : Ω k → H k such that I k (z), I k (z ) = f k (z), f k (z ) H k for all z, z .
As mentioned previously, the mapping to the RKHS is a key to the multilayer construction. Given I k , larger image neighborhoods are represented by patches of size e k × e k that can be mapped to a point in the Cartesian product space H e k ×e k k endowed with its natural inner-product; finally, the kernel K k+1 defined on these patches can be seen as a kernel on larger image neighborhoods than K k .

3 End-to-End Kernel Learning with Supervised CKNs

In the previous section, we have described a variant of convolutional kernel networks where linear subspaces are learned at every layer. This is achieved without supervision by a K-means algorithm leading to small projection residuals. It is thus natural to introduce also a discriminative approach.

Backpropagation Rules for Convolutional Kernel Networks

We now consider a prediction task, where we are given a training set of images I 1 0 , I 2 0 , . . . , I n 0 with respective scalar labels y 1 , . . . , y n living either in {-1; +1} for binary classification and R for regression. For simplicity, we only present these two settings here, but extensions to multiclass classification and multivariate regression are straightforward. We also assume that we are given a smooth convex loss function L : R × R → R that measures the fit of a prediction to the true label y.

Given a positive definite kernel K on images, the classical empirical risk minimization formulation consists of finding a prediction function in the RKHS H associated to K by minimizing the objective

min f ∈H 1 n n i=1 L(y i , f (I i 0 )) + λ 2 f 2 H , (5)
where the parameter λ controls the smoothness of the prediction function f with respect to the geometry induced by the kernel, hence regularizing and reducing overfitting [START_REF] Schölkopf | Learning with kernels: support vector machines, regularization, optimization, and beyond[END_REF]. After training a convolutional kernel network with k layers, such a positive definite kernel may be defined as

K Z (I 0 , I 0 ) = z∈Ω k f k (z), f k (z) H k = z∈Ω k I k (z), I k (z) , (6)
where I k , I k are the k-th finite-dimensional feature maps of I 0 , I 0 , respectively, and f k , f k the corresponding maps in Ω k → H k , which have been defined in the previous section. The kernel is also indexed by Z, which represents the network parameters-that is, the subspaces F 1 , . . . , F k , or equivalently the set of filters Z 1 , . . . , Z k from Eq. ( 3). Then, formulation (5) becomes equivalent to

min W∈R p k ×|Ω k | 1 n n i=1 L(y i , W, I i k ) + λ 2 W 2 F , (7)
where . F is the Frobenius norm that extends the Euclidean norm to matrices, and, with an abuse of notation, the maps I i k are seen as matrices in R p k ×|Ω k | . Then, the supervised convolutional kernel network formulation consists of jointly minimizing [START_REF] Dong | Learning a deep convolutional network for image super-resolution[END_REF] with respect to W in R p k ×|Ω k | and with respect to the set of filters Z 1 , . . . , Z k , whose columns are constrained to be on the Euclidean sphere.

Computing the derivative with respect to the filters Z 1 , . . . , Z k . Since we consider a smooth loss function L, e.g., logistic, squared hinge, or square loss, optimizing [START_REF] Dong | Learning a deep convolutional network for image super-resolution[END_REF] with respect to W can be achieved with any gradient-based method. Moreover, when L is convex, we may also use fast dedicated solvers, (see, e.g., [START_REF] Lin | A universal catalyst for first-order optimization[END_REF], and references therein). Optimizing with respect to the filters Z j , j = 1, . . . , k is more involved because of the lack of convexity. Yet, the objective function is differentiable, and there is hope to find a "good" stationary point by using classical stochastic optimization techniques that have been successful for training deep networks.

For that, we need to compute the gradient by using the chain rule-also called "backpropagation" [START_REF] Le Cun | Efficient backprop[END_REF]. We instantiate this rule in the next lemma, which we have found useful to simplify the calculation.

Lemma 1 (Perturbation view of backpropagration.) Consider an image I 0 represented here as a matrix in R p0×|Ω0| , associated to a label y in R and call I Z k the k-th feature map obtained by encoding I 0 with the network parameters Z. Then, consider a perturbation E = {ε 1 , . . . , ε k } of the set of filters Z. Assume that we have for all j ≥ 0,

I Z+E j = I Z j + ∆I Z,E j + o( E ), ( 8
)
where E is equal to k l=1 ε l F , and ∆I Z,E j is a matrix in R pj ×|Ωj | such that for all matrices U of the same size,

∆I Z,E j , U = ε j , g j (U) + ∆I Z,E j-1 , h j (U) , (9)
where the inner-product is the Frobenius's one and g j , h j are linear functions. Then,

∇ Zj L(y, W, I Z k ) = L (y, W, I Z k ) g j (h j+1 (. . . h k (W)), (10)
where L denote the derivative of the smooth function L with respect to its second argument.

The proof of this lemma is straightforward and follows from the definition of the Fréchet derivative. Nevertheless, it is useful to derive the closed form of the gradient in the next proposition.

Proposition 1 (Gradient of the loss with respect to the the filters Z 1 , . . . , Z k .) Consider the quantities introduced in Lemma 1, but denote I Z j by I j for simplicity. By construction, we have for all j ≥ 1,

I j = A j κ j (Z j E j (I j-1 )S -1 j )S j P j , (11)
where I j is seen as a matrix in R pj ×|Ωj | ; E j is the linear operator that extracts all overlapping e j-1 × e j-1 patches from a map such that E j (I j-1 ) is a matrix of size p j-1 e 2 j-1 × |Ω j-1 |; S j is a diagonal matrix whose diagonal entries carry the 2 -norm of the columns of E j (I j-1 ); A j is short for κ j (Z j Z j ) -1/2 ; and P j is a matrix of size |Ω j-1 | × |Ω j | performing the linear pooling operation. Then, the gradient of the loss with respect to the filters Z j , j = 1, . . . , k is given by [START_REF] He | Deep residual learning for image recognition[END_REF] with

g j (U) = E j (I j-1 )B j - 1 2 Z j κ j (Z j Z j ) (C j + C j ) h j (U) = E j Z j B j + E j (I j-1 ) S -2 j M j UP j -E j (I j-1 ) Z j B j , ( 12
)
where U is any matrix of the same size as I j , M j = A j κ j (Z j E j (I j-1 )S -1 j )S j is the j-th feature map before the pooling step, is the Hadamart (elementwise) product, E j is the adjoint of E j , and

B j = κ j Z j E j (I j-1 )S -1 j A j UP j and C j = A 1/2 j I j U A 3/2 j . ( 13
)
The proof is presented in Appendix B. Most quantities that appear above admit physical interpretations: multiplication by P j performs downsampling; multiplication by P j performs upsampling; multiplication of E j (I j-1 ) on the right by S -1 j performs 2 -normalization of the columns; Z j E j (I j-1 ) can be seen as a spatial convolution of the map I j-1 by the filters Z j ; finally, E j "combines" a set of patches into a spatial map by adding to each pixel location the respective patch contributions.

Computing the gradient requires a forward pass to obtain the maps I j through [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF] and a backward pass that composes the functions g j , h j as in [START_REF] He | Deep residual learning for image recognition[END_REF]. The complexity of the forward step is dominated by the convolutions Z j E j (I j-1 ), as in convolutional neural networks. The cost of the backward pass is the same as the forward one up to a constant factor. Assuming p j ≤ |Ω j-1 |, which is typical for lower layers that require more computation than upper ones, the most expensive cost is due to E j (I j-1 )B j and Z j B j which is the same as Z j E j (I j-1 ). We also pre-compute A 1/2 j and A 3/2 j by eigenvalue decompositions, whose cost is reasonable when performed only once per minibatch. Off-diagonal elements of M j UP j -E j (I j-1 ) Z j B j are also not computed since they are set to zero after elementwise multiplication with a diagonal matrix. In practice, we also replace A j by (κ j (Z j Z j ) + εI) -1/2 with ε = 0.001, which corresponds to performing a regularized projection onto F j (see Appendix A). Finally, a small offset of 0.00001 is added to the diagonal entries of S j .

Optimizing hyper-parameters for RBF kernels. When using the kernel (2), the objective is differentiable with respect to the hyper-parameters α j . When large amounts of training data are available and overfitting is not a issue, optimizing the training loss by taking gradient steps with respect to these parameters seems appropriate instead of using a canonical parameter value. Otherwise, more involved techniques may be needed; we plan to investigate other strategies in future work.

Optimization and Practical Heuristics

The backpropagation rules of the previous section have set up the stage for using a stochastic gradient descent method (SGD). We now present a few strategies to accelerate it in our context.

Hybrid convex/non-convex optimization. Recently, many incremental optimization techniques have been proposed for solving convex optimization problems of the form [START_REF] Dong | Learning a deep convolutional network for image super-resolution[END_REF] when n is large but finite (see [START_REF] Lin | A universal catalyst for first-order optimization[END_REF] and references therein). These methods usually provide a great speed-up over the stochastic gradient descent algorithm without suffering from the burden of choosing a learning rate. The price to pay is that they rely on convexity, and they require storing into memory the full training set. For solving [START_REF] Dong | Learning a deep convolutional network for image super-resolution[END_REF] with fixed network parameters Z, it means storing the n maps I i k , which is often reasonable if we do not use data augmentation. To partially leverage these fast algorithms for our non-convex problem, we have adopted a minimization scheme that alternates between two steps: (i) fix Z, then make a forward pass on the data to compute the n maps I i k and minimize the convex problem [START_REF] Dong | Learning a deep convolutional network for image super-resolution[END_REF] with respect to W using the accelerated MISO algorithm [START_REF] Lin | A universal catalyst for first-order optimization[END_REF]; (ii) fix W, then make one pass of a projected stochastic gradient algorithm to update the k set of filters Z j . The set of network parameters Z is initialized with the unsupervised learning method described in Section 2.

Preconditioning on the sphere. The kernels κ j are defined on the sphere; therefore, it is natural to constrain the filters-that is, the columns of the matrices Z j -to have unit 2 -norm. As a result, a classical stochastic gradient descent algorithm updates at iteration t each filter z as follows z ← Proj . 2 =1 [z-η t ∇ z L t ], where ∇ z L t is an estimate of the gradient computed on a minibatch and η t is a learning rate. In practice, we found that convergence could be accelerated by preconditioning, which consists of optimizing after a change of variable to reduce the correlation of gradient entries. For unconstrained optimization, this heuristic involves choosing a symmetric positive definite matrix Q and replacing the update direction ∇ z L t by Q∇ z L t , or, equivalently, performing the change of variable z = Q 1/2 z and optimizing over z . When constraints are present, the case is not as simple since Q∇ z L t may not be a descent direction. Fortunately, it is possible to exploit the manifold structure of the constraint set (here, the sphere) to perform an appropriate update [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. Concretely, (i) we choose a matrix Q per layer that is equal to the inverse covariance matrix of the patches from the same layer computed after the initialization of the network parameters. (ii) We perform stochastic gradient descent steps on the sphere manifold after the change of variable

z = Q 1/2 z , leading to the update z ← Proj . 2=1 [z -η t (I -(1/z Qz)Qzz )Q∇ z L t ].
Because this heuristic is not a critical component, but simply an improvement of SGD, we relegate mathematical details in Appendix C.

Automatic learning rate tuning. Choosing the right learning rate in stochastic optimization is still an important issue despite the large amount of work existing on the topic, see, e.g., [START_REF] Le Cun | Efficient backprop[END_REF] and references therein. In our paper, we use the following basic heuristic: the initial learning rate η t is chosen "large enough"; then, the training loss is evaluated after each update of the weights W. When the training loss increases between two epochs, we simply divide the learning rate by two, and perform "back-tracking" by replacing the current network parameters by the previous ones.

Active-set heuristic. For classification tasks, "easy" samples have often negligible contribution to the gradient (see, e.g., [START_REF] Le Cun | Efficient backprop[END_REF]). For instance, for the squared hinge loss L(y, ŷ) = max(0, 1y ŷ) 2 , the gradient vanishes when the margin y ŷ is greater than one. This motivates the following heuristic: we consider a set of active samples, initially all of them, and remove a sample from the active set as soon as we obtain zero when computing its gradient. In the subsequent optimization steps, only active samples are considered, and after each epoch, we randomly reactivate 10% of the inactive ones.

Experiments

We now present experiments on image classification and super-resolution. All experiments were conducted on 8-core and 10-core 2.4GHz Intel CPUs using C++ and Matlab.

Image Classification on "Deep Learning" Benchmarks

We consider the datasets CIFAR-10 [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] and SVHN [START_REF] Netzer | Reading digits in natural images with unsupervised feature learning[END_REF], which contain 32 × 32 images from 10 classes. CIFAR-10 is medium-sized with 50 000 training samples and 10 000 test ones. SVHN is larger with 604 388 training examples and 26 032 test ones. We evaluate the performance of a 9-layer network, designed with few hyper-parameters: for each layer, we learn 512 filters and choose the RBF kernels κ j defined in (2) with initial parameters α j = 1/(0.5 2 ). Layers 1, 3, 5, 7, 9 use 3×3 patches and a subsampling pooling factor of √ 2 except for layer 9 where the factor is 3; Layers 2, 4, 6, 8 use simply 1 × 1 patches and no subsampling. For CIFAR-10, the parameters α j are kept fixed during training, and for SVHN, they are updated in the same way as the filters. We use the squared hinge loss in a one vs all setting to perform multi-class classification (with shared filters Z between classes). The input of the network is pre-processed with the local whitening procedure described in [START_REF] Paulin | Local convolutional features with unsupervised training for image retrieval[END_REF]. We use the optimization heuristics from the previous section, notably the automatic learning rate scheme, and a gradient momentum with parameter 0.9, following [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF]. The regularization parameter λ and the number of epochs are set by first running the algorithm on a 80/20 validation split of the training set. λ is chosen near the canonical parameter λ = 1/n, in the range 2 i /n, with i = -4, . . . , 4, and the number of epochs is at most 100. The initial learning rate is 10 with a minibatch size of 128.

We present our results in Table 1 along with the performance achieved by a few recent methods without data augmentation or model voting/averaging. In this context, the best published results are obtained by the generalized pooling scheme of [START_REF] Lee | Generalizing pooling functions in convolutional neural networks: Mixed, gated, and tree[END_REF]. We achieve about 2% test error on SVHN and about 10% on CIFAR-10, which positions our method as a reasonably "competitive" one, in the same ballpark as the deeply supervised nets of [START_REF] Lee | Deeply-supervised nets[END_REF] or network in network of [START_REF] Lin | Network in network[END_REF]. Due to lack of space, the results reported here only include a single supervised model. Preliminary experiments with no supervision show also that one may obtain competitive accuracy with wide shallow architectures. For instance, a two-layer network with (1024-16384) filters achieves 14.2% error on CIFAR-10. Note also that our unsupervised model outperforms original CKNs [START_REF] Mairal | Convolutional kernel networks[END_REF]. The best single model from [START_REF] Mairal | Convolutional kernel networks[END_REF] gives indeed 21.7%. Training the same architecture with our approach is two orders of magnitude faster and gives 19.3%. Another aspect we did not study is model complexity.

Here as well, preliminary experiments are encouraging. Reducing the number of filters to 128 per layer yields indeed 11.95% error on CIFAR-10 and 2.15% on SVHN. A more precise comparison with no supervision and with various network complexities will be presented in another venue.

Image Super-Resolution from a Single Image

The blurry estimates are processed by a 9-layer network, with 3 × 3 patches and 128 filters at every layer without linear pooling and zero-padding. Pixel values are predicted with a linear model applied to the 128-dimensional vectors present at every pixel location of the last layer, and we use the square loss to measure the fit. The optimization procedure and the kernels κ j are identical to the ones used for processing the SVHN dataset in the classification task. The pipeline also includes a pre-processing step, where we remove from input images a local mean component obtained by convolving the images with a 5 × 5 averaging box filter; the mean component is added back after up-scaling.

For the evaluation, we consider three datasets: Set5 and Set14 are standard for super-resolution; Kodim is the Kodak Image database, available at http://r0k.us/graphics/kodak/, which contains high-quality images with no compression or demoisaicing artefacts. The evaluation procedure follows [START_REF] Dong | Learning a deep convolutional network for image super-resolution[END_REF][START_REF] Dong | Image super-resolution using deep convolutional networks[END_REF][START_REF] Timofte | Anchored neighborhood regression for fast example-based superresolution[END_REF][START_REF] Wang | Deep networks for image super-resolution with sparse prior[END_REF] by using the code from the author's web page. We present quantitative results in Table 2. For x3 upscaling, we simply used twice our model learned for x2 upscaling, followed by a 3/4 downsampling. This is clearly suboptimal since our model is not trained to up-scale by a factor 3, but this naive approach still outperforms other baselines [START_REF] Dong | Learning a deep convolutional network for image super-resolution[END_REF][START_REF] Dong | Image super-resolution using deep convolutional networks[END_REF][START_REF] Wang | Deep networks for image super-resolution with sparse prior[END_REF] that are trained end-to-end. Note that [START_REF] Wang | Deep networks for image super-resolution with sparse prior[END_REF] also proposes a data augmentation scheme at test time that slightly improves their results. In Appendix D, we also present a visual comparison between our approach and [START_REF] Dong | Image super-resolution using deep convolutional networks[END_REF], whose pipeline is the closest to ours, up to the use of a supervised kernel machine instead of CNNs. A Orthogonal Projection on the Finite-Dimensional Subspace F 1

First, we remark that the kernel K 1 is homogeneous such that for every patch x and scalar γ > 0,

ϕ 1 (γx) = γϕ 1 (x).
Thus, we may assume x to have unit 2 -norm without loss of generality and perform the projection on F 1 = Span(ϕ(z 1 ), . . . , ϕ(z p1 )) of the normalized patch, before applying the inverse rescaling.

Then, let us denote by f x the orthogonal projection of a patch x with unit 2 -norm defined as

f x := arg min f ∈F1 ϕ 1 (x) -f 2 H1 ,
which is equivalent to

f x := p1 j=1 α j ϕ 1 (z j ) with α ∈ arg min α∈R p 1 ϕ 1 (x) - p1 j=1 α j ϕ 1 (z j ) 2 H1
.

After short calculation, we obtain

f x = p1 j=1 α j ϕ 1 (z j ) with α ∈ arg min α∈R p 1 1 -2α κ 1 (z x) + α κ 1 (Z Z)α ,
since the vectors z j provided by the spherical K-means algorithm have unit 2 -norm. Assuming κ 1 (Z Z) to be invertible, we have α = κ 1 (Z Z) -1 κ 1 (Z x). After projection, normalized patches x, x may be parametrized by α = κ 1 (Z Z) -1 κ 1 (Z x) and α = κ 1 (Z Z) -1 κ 1 (Z x ), respectively. Then, we have

f x , f x H1 = α κ 1 (Z Z)α = ψ 1 (x), ψ 1 (x ) ,
which is the desired result.

When κ 1 (Z Z) is not invertible or simply badly conditioned, it is also common to use instead

ψ 1 (x) = κ 1 (Z Z) + εI -1/2 κ 1 (Z x),
where ε > 0 is a small regularization that improves the condition number of κ 1 (Z Z). Such a modification can be interpreted as performing a slightly regularized projection onto the finitedimensional subspace F 1 .

B Computation of the Gradient with Respect to the Filters

To compute the gradient of the loss function, we use Lemma 1 and start by analyzing the effect of perturbing every quantity involved in [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF] such that we may obtain the desired relations ( 8) and [START_REF] Goodfellow | Maxout networks[END_REF]. Before proceeding, we recall the definition of the set Z + E = {Z 1 + ε 1 , . . . , Z k + ε k } and the precise definition of the Landau notation o( E ), which we use in [START_REF] Dong | Image super-resolution using deep convolutional networks[END_REF]. Here, it simply means a quantity that is negligible in front of the norm

E = k j=1 ε j F -that is, lim E→0 I Z+E j -I Z j -∆I Z,E j F E = 0.
Then, we start by initializing a recursion: I Z 0 is unaffected by the perturbation and thus ∆I Z,E 0 = 0. Consider now an index j > 0 and assume that (8) holds for j -1 with ∆I Z,E j-1 = O( E ). First, we remark that E j (I Z+E j-1 ) = E j (I Z j-1 ) + E j (∆I Z,E j-1 ) + o( E ). Then, the diagonal matrix S j becomes after perturbation

S j + S -1 j E j (I Z j-1 ) E j (∆I Z,E j-1 ) ∆Sj +o( E ).
The inverse diagonal matrix S -1 j becomes

S -1 j -S -3 j E j (I Z j-1 ) E j (∆I Z,E j-1 ) ∆S -1 j +o( E ),
and the matrix A j becomes κ j (Z j +ε j ) (Z j +ε j )

-1 2 = κ j Z j Z j +ε j Z j +Z j ε j + o( ε j F )) -1 2 = κ j Z j Z j +κ j Z j Z j Z j ε j +ε j Z j +o( E ) -1 2 = A 1 2 j I+A j κ j Z j Z j Z j ε j +ε j Z j A j +o( E ) -1/2 A 1 2 j = A j - 1 2 A 3 2 j κ j Z j Z j Z j ε j +ε j Z j A 3 2 j ∆Aj +o( E ),
where we have used the relation

(I + Q) -1/2 = I -1 2 Q + o( Q F ).
Note that the quantities ∆A j , ∆S j , ∆S -1 j that we have introduced are all O( E ). Then, by replacing the quantities A j , S j , S -1 j , I j-1 by their perturbed versions in the definition of I j given in [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF], we obtain that

I Z+E j is equal to (A j +∆A j )κ j (Z j +ε j ) E j (I Z j-1 )+E j (∆I Z,E j-1 ) (S -1 j +∆S -1 j ) (S j +∆S j )P j +o( E ).
Then, after short calculation, we obtain the desired relation

I Z+E j = I Z+E j + ∆I Z,E j + o( E ) with ∆I Z,E j =∆A j κ j (Z j E j (I Z j-1 )S -1 j )S j P j + A j κ j (Z j E j (I Z j-1 )S -1 j ) (ε j E j (I Z j- 1
)) P j + A j κ j (Z j E j (I Z j-1 )S -1 j ) (Z j E j (∆I Z,E j-1 )) P j + A j κ j (Z j E j (I Z j-1 )S -1 j ) (Z j E j (I Z j-1 )∆S -1 j S j ) P j + A j κ j (Z j E j (I Z j-1 )S -1 j )∆S j P j .

First, we remark that ∆I Z,E j = O( E ), which is one of the induction hypothesis we need. Then, after plugging in the values of ∆A j , ∆S j , ∆S -1 j , and with further simplification, we obtain

∆I Z,E j = - 1 2 A 3 2 j κ j Z j Z j Z j ε j + ε j Z j A 1 2 j I Z j + A j κ j (Z j E j (I Z j-1 )S -1 j ) ε j E j (I Z j-1 ) P j + A j κ j (Z j E j (I Z j-1 )S -1 j ) Z j E j (∆I Z,E j-1 ) P j -A j κ j (Z j E j (I Z j-1 )S -1 j ) Z j E j (I Z j-1 ) S -2 j E j (I Z j-1 ) E j (∆I Z,E j-1 ) P j + M Z j S -2 j E j (I Z j-1 ) E j (∆I Z,E j-1 ) P j ,
where M Z j is the j-th feature map of I 0 before the j-th linear pooling step-that is, I Z j = M Z j P j . We now see that ∆I Z,E j is linear in ε j and ∆I Z,E j-1 , which guarantees that there exist two linear functions g j , h j that satisfy [START_REF] Goodfellow | Maxout networks[END_REF]. More precisely, we want for all matrix U of the same size as ∆I Z,E j ε j , g j (U) = -

1 2 A 3 2 j κ j Z j Z j Z j ε j + ε j Z j A 1 2 j I Z j , U + A j κ j (Z j E j (I Z j-1 )S -1 j ) ε j E j (I Z j-1 ) P j , U , and
∆I Z,E j-1 , h j (U) = A j κ j (Z j E j (I Z j-1 )S -1 j ) Z j E j (∆I Z,E j-1 ) P j , U -A j κ j (Z j E j (I Z j-1 )S -1 j ) Z j E j (I Z j-1 ) S -2 j E j (I Z j-1 ) E j (∆I Z,E j-1 ) P j , U + M Z j S -2 j E j (I Z j-1 ) E j (∆I Z,E j-1 ) P j , U .
Then, it is easy to obtain the form of g j , h j given in [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], by using in the right order the following elementary calculus rules: (i) UV, W = U, WV = V, U W , (ii) U, V = U , V , (iii) U V, W = U, V W for any matrices U, V, W of appropriate sizes, and also (iv) E j (U), V = U, E j (V) , by definition of the adjoint operator. We conclude by induction.

C Preconditioning Heuristic on the Sphere

In this section, we present a preconditioning heuristic for optimizing over the sphere S p-1 , inspired by second-order (Newton) optimization techniques on smooth manifolds [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. Following [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF], we will consider gradient descent steps on the manifold. A fundamental operation is thus the projection operator P z onto the tangent space at a point z. This operator is defined for the sphere by

P z [u] = (I -zz )u,
for any vector u in R p . Another important operator is the Euclidean projection on S p-1 , which was denoted by Proj . 2=1 in previous parts of the paper.

Gradient descent on the sphere S p-1 is equivalent to the projected gradient descent in R p . When optimizing on a manifold, the natural descent direction is the projected gradient P z ∇L(z). In the case of the sphere, a gradient step on the manifold is equivalent to a classical projected gradient descent step in R p with particular step size:

Proj . 2=1 [z -ηP z [∇L(z)]] = Proj . 2=1 z -η I -zz ∇L(z) = Proj . 2=1 1 + ηz ∇L(z) z -η∇L(z) = Proj . 2=1 z - η 1 + ηz ∇L(z)
∇L(z) .

In R p with no constraint, pre-conditioning is equivalent to performing a change of variable.

For unconstrained optimization in R p , faster convergence is usually achieved when one has access to an estimate of the inverse of the Hessian ∇ 2 L(z)-assuming twice differentiability-and using the descent direction (∇ 2 L(z)) -1 ∇L(z) instead of ∇L(z); then, we obtain a Newton method. When the exact Hessian is not available, or too costly to compute and/or invert, it is however common to use instead a constant estimate of the inverse Hessian, denoted here by Q, which we call pre-conditioning matrix. Finding an appropriate matrix Q is difficult in general, but for learning linear models, a typical choice is to use the inverse covariance matrix of the data (or one approximation). In that case, the preconditioned gradient descent step consists of the update z -ηQ∇L(z). Such a matrix Q is defined similarly in the context of convolutional kernel networks, as explained in the main part of the paper. A useful interpretation of preconditioning is to see it as optimizing after a change of variable. Define indeed the objective L(w) = L(Q 1/2 w).

Then, minimizing L is equivalent to minimizing L with respect to z, with the relation z = Q 1/2 w. Moreover, when there is no constraint on z and w, the regular gradient descent algorithm on L is equivalent to the preconditioned gradient descent on L:

w ← w -η∇ L(w) ⇐⇒ w ← w -ηQ 1/2 ∇L(Q 1/2 w)

⇐⇒ z ← z -ηQ∇L(z) with z = Q 1/2 w.

We remark that the Hessian ∇ 2 L(w) is equal to Q 1/2 ∇ 2 L(Q 1/2 w)Q 1/2 , which is equal to identity when Q coincides with the inverse Hessian of L. In general, this is of course not the case, but the hope is to obtain a Hessian ∇ 2 L that is better conditioned than ∇ 2 L, thus resulting in faster convergence.

Bicubic

Sparse coding [START_REF] Zeyde | On single image scale-up using sparse-representations[END_REF] CNN2 [START_REF] Dong | Image super-resolution using deep convolutional networks[END_REF] SCKN (Ours)

Figure 2: Visual comparison for x3 image up-scaling. Each column corresponds to a different method (see bottom row). RGB images are converted to the YCbCr color space and the up-scaling method is applied to the luminance channel only. Color channels are up-scaled using bicubic interpolation for visualization purposes. CNN2 and SCKN perform similarly in textured areas, but SCKN provides significantly sharper artefact-free edges (see in particular the butterfly image). Best seen by zooming on a computer screen with an appropriate PDF viewer that does not smooth the image content.

Bicubic

Sparse coding [START_REF] Zeyde | On single image scale-up using sparse-representations[END_REF] CNN2 [START_REF] Dong | Image super-resolution using deep convolutional networks[END_REF] SCKN (Ours)

Figure 3 :

 3 Figure 3: Another visual comparison for x3 image up-scaling. See caption of Figure 2. Best seen by zooming on a computer screen with an appropriate PDF viewer that does not smooth the image content.

Table 1 :

 1 Test error in percents reported by a few recent publications on the CIFAR-10 and SVHN datasets without data augmentation or model voting/averaging.

		Stoch P. [29] MaxOut [9]	NiN [17]	DSN [15]	Gen P. [14]	SCKN (Ours)
	CIFAR-10	15.13	11.68	10.41	9.69	7.62	10.20
	SVHN	2.80	2.47	2.35	1.92	1.69	2.04

Table 2 :

 2 Reconstruction accuracy for super-resolution in PSNR (the higher, the better). All CNN approaches are without data augmentation at test time. See Appendix D for the SSIM quality measure.

	Fact. Dataset Bicubic SC [30] ANR [26] A+[26] CNN1 [7] CNN2 [8] CSCN [27] SCKN
		Set5	33.66	35.78	35.83	36.54	36.34	36.66	36.93	37.07
	x2	Set14	30.23	31.80	31.79	32.28	32.18	32.45	32.56	32.76
		Kodim	30.84	32.19	32.23	32.71	32.62	32.80	32.94	33.21
		Set5	30.39	31.90	31.92	32.58	32.39	32.75	33.10	33.08
	x3	Set14	27.54	28.67	28.65	29.13	29.00	29.29	29.41	29.50
		Kodim	28.43	29.21	29.21	29.57	29.42	29.64	29.76	29.88

To simplify, we use zero-padding when patches are close to the image boundaries, but this is optional.

Acknowledgments

This work was supported by ANR (MACARON project ANR-14-CE23-0003-01).

Preconditioning on a smooth manifold requires some care. Unfortunately, using second-order information (or simply a pre-conditioning matrix) when optimizing over a constraint set or over a smooth manifold is not as simple as optimizing in R p since the quantities Q∇L(z), P z [Q∇L(z)], QP z [∇L(z)] may not be feasible descent directions. However, the point of view that sees pre-conditioning as a change of variable will give us the right direction to follow.

Optimizing L on S p-1 is in fact equivalent to optimizing L on the smooth manifold

which represents an ellipsoid. The tangent plane at a point w of Sp-1 being defined by the normal vector Qw/ Qw 2 , it is then possible to introduce the projection operator Pw on the tangent space:

Then, we may define the gradient descent step rule on Sp-1 as

With the change of variable z = Q 1/2 w, this is equivalent to

This is exactly the update rule we have chosen in our paper, as a heuristic in a stochastic setting.

D Additional Results for Image Super-Resolution

We present a quantitative comparison in Table 3 using the structural similarity index measure (SSIM), which is known to better reflect the quality perceived by humans than the PSNR; it is commonly used to evaluate the quality of super-resolution methods, see [START_REF] Dong | Image super-resolution using deep convolutional networks[END_REF][START_REF] Timofte | Anchored neighborhood regression for fast example-based superresolution[END_REF][START_REF] Wang | Deep networks for image super-resolution with sparse prior[END_REF]. Then, we present a visual comparison between several approaches in Figures 2,3, and 4. We focus notably on the classical convolutional neural network of [START_REF] Dong | Image super-resolution using deep convolutional networks[END_REF] since our pipeline essentially differs in the use of our supervised kernel machine instead of convolutional neural networks. After subjective evaluation, we observe that both methods perform equally well in textured areas. However, our approach recovers better thin high-frequency details, such as the eyelash of the baby in the first image. By zooming on various parts, it is easy to notice similar differences in other images. We also observed a few ghosting artefacts near object boundaries with the method of [START_REF] Dong | Image super-resolution using deep convolutional networks[END_REF], which is not the case with our approach.

Table 3: Reconstruction accuracy of various super-resolution approaches. The numbers represent the structural similarity index (SSIM), the higher, the better.

Fact. Dataset Bicubic SC [START_REF] Zeyde | On single image scale-up using sparse-representations[END_REF] ANR [START_REF] Timofte | Anchored neighborhood regression for fast example-based superresolution[END_REF] A+ [START_REF] Timofte | Anchored neighborhood regression for fast example-based superresolution[END_REF]