F Rousseau

A Non-local Approach for Image Super-Resolution using Intermodality Priors

, or generating a high resolution brain image from an input low-resolution image, with the help of another high resolution brain image. We propose an approach for image super-resolution by using anatomical intermodality priors from a reference image. Contrary to interpolation techniques, in order to be able to recover fine details in images, the reconstruction process is based on a physical model of image acquisition. Another contribution to this inverse problem is a new regularization approach that uses an example-based framework integrating non-local similarity constraints to handle in a better way repetitive structures and texture. The effectiveness of our approach is demonstrated by experiments on realistic Brainweb Magnetic Resonance images and on clinical images from ADNI, generating automatically high-quality brain images from low-resolution input.

Introduction

Image processing algorithm performances are often limited by the image resolution. Resolution is for instance a key point in brain segmentation in Magnetic Resonance (MR) imaging for which partial volume effect is a limiting factor for fine image analysis. Then, it clearly appears that improving image resolution is still one of the main challenges in medical image processing. In medical imaging, a so-called low resolution (LR) 3D image is usually a stack of 2D thick slices. As a result, 3D data are generally not isotropic. Figure 1 shows for instance a high resolution (HR) T1-weighted image and a LR T2-weighted image of the same patient (in this case, these images come from the Brainweb database [START_REF] Cocosco | Online Interface to a 3D MRI Simulated Brain Database Proceedings of 3-rd International Conference on Functional Mapping of the[END_REF]). In many cases, it may be necessary to put HR and LR images into a same coordinate system. In medical imaging, this is usually done by applying interpolation techniques [15]. Image interpolation is a very common image processing technique in medical imaging pipelines and it may have a strong impact on other processing steps such as segmentation or registration. Interpolation methods can be divided into two groups: scene-based and object-based methods [START_REF] Grevera | An objective comparison of 3-D image interpolation methods[END_REF]. Scene-based approaches use only image intensities to determine the interpolated intensity (for instance: nearest neighbor, linear interpolation, spline-based interpolation). Such scene-based methods produce perceptually unsatisfactory results with blurred edges and textures. Many edge-preserving interpolation techniques have been reported to handle this problem. However, these techniques rely on accurate edge information that is not obtainable from coarse data. In order to guide the interpolation process, object-based methods make use of additional information extracted from images. As an example, non-rigid registration techniques have been used to drive interpolation methods between adjacent slices [START_REF] Frakes | A New Method for Registration-Based Medical Image Interpolation[END_REF], [START_REF] Penney | Registration-Based Interpolation[END_REF]. However, scene-based and object-based techniques do not take advantage of a model of the imaging process.

Super-resolution (SR) is another model-based technique which relies on modeling the imaging processes and using regularization methods describing a priori constraints. The principle of this approach is usually combining LR images to produce an image that has a higher spatial resolution than the original images [START_REF] Bose | Special Issue: High Resolution Image Reconstruction Int[END_REF]. In medical imaging, several SR methods have been proposed to combine LR images to reconstruct one HR image [START_REF] Peeters | The use of super-resolution techniques to reduce slice thickness in functional[END_REF][START_REF] Carmi | [END_REF][START_REF] Rousseau | Registration-Based Approach for Reconstruction of High-Resolution in Utero Fetal[END_REF] (usually by modifying the acquisition protocol). SR is a large research field encompassing many applications. The work we present in this paper is related to the single-frame SR framework [START_REF] Van Ouwerkerk | Image super-resolution survey Image and[END_REF], meaning that only one LR image is used to generate an HR image. More specifically, we focus on studies involving MR imaging for which an anatomical HR image and several other LR images are acquired to keep acquisition Figure 1: Example of MR Brainweb data [START_REF] Cocosco | Online Interface to a 3D MRI Simulated Brain Database Proceedings of 3-rd International Conference on Functional Mapping of the[END_REF]. First row: high resolution T1-weighted image (1×1mm 2 in-plane resolution, 1mm slice thickness); second row: low resolution T2-weighted image (1×1mm 2 in-plane resolution, 3mm slice thickness). The purpose of the proposed approach is to reconstruct a high-resolution T2-weighted image using information from the T1-weighted image.

time at an acceptable level for the patient (see Figure 1). This is the case for routinely performed clinical MR acquisitions such as follow-up of neurodegenerative diseases, brain tumor evolution or diffusion tensor imaging. Typically, one isotropic HR T1-weighted image and several anisotropic LR images (such as T2-weighted, FLAIR or proton density images) are acquired. One can note that this issue of relative difference in image resolution may also occur in other clinical settings. For instance, patient follow-up may require an analysis of scans taken at different magnetic field strength, potentially at different clinical sites.

In such context, we propose a new approach for image SR by using information from an HR MR image to drive the image reconstruction of the LR MR image. The general framework developed in this paper can be applied to other domains where image resolution is an important issue, such as remote sensing. This paper is built on previously published work [START_REF] Rousseau | Brain Hallucination European Conference on Computer Vision, Part I[END_REF]. Many experiments have been added in order to provide further insight to the proposed approach. In Section 2, we present the image SR problem using a model-based framework and some recently proposed example-based approaches. Section 3 details our non-local approach for image SR using an HR reference MR image. In Section 4, results obtained on realistic brainweb images are presented and discussed.

Image super-resolution

In this section, we present the model-based framework for image SR which leads to an ill-posed inverse problem. Then, recently proposed regularization approaches are described.

Model-based framework

Contrary to interpolation approaches, model-based approaches use a generic observation model such as:

y = DBW x + n (1)
where y denotes the LR image, x is the HR image, n represents observation noise, D is the sub-sampling matrix, B a blur matrix, W is the geometric transformation.

The three operators can be combined into a single matrix H: H = DBW . The matrix H thus incorporates motion compensation, degradation effects, and sub-sampling for the LR image y. In this paper, we assume that the degradation operator H and the noise characteristics are known and we focus only on the reconstruction step.

Based on this observation model, the SR image can be estimated by minimizing a least-square cost function such as:

x = arg min x y -Hx 2 .

(2)

For such an inverse problem, some form of regularization plays a crucial role and must be included in the cost function to stabilize the problem or constrain the space of solutions [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF]. Thus, the HR image is computed by considering the following equation:

x = arg min x L(x, y, H) + λR(x). (3)
L(x, y, H) is a data fidelity term related to the physical model that penalizes inconsistency between the estimated HR image x and the observed LR image y. This term depends on the observation model used. In our case, the data fidelity term is simply defined as following:

L(x, y, H) = y -Hx 2 . (4)
R(x) is a regularization term. A common approach for regularization is to take explicitly into account the image geometry and to introduce a global weight λ that balances the contribution of prior smoothness terms and a fidelity term.

Image regularization

Examples of popular pixel-based regularizers are:

• Tikhonov regularization: R(x) = P p=0 Ω c p |x (p) | 2 where c p are positive coefficients, Ω is the domain of integration (all the voxels) and x (p) is p th -order derivative of x.

• Variational approach: R(x) = P p=0 Ω c p |x (p) | q where q is different from 2 (q = 2 is the Tikhonov regularization).

• Markov random field a priori image model: R(x) = α s φ(x s ) + β rs φ(x rx s ) where φ can be a L 2 L 1 function for instance.

However, the use of such a regularization term assumes a specific image model: for instance, variational approaches can be based on the assumption that images are made of smooth regions separated by sharp edges. Such prior regularization approach introduces then strong constraints into the reconstruction process.

Explicit models for the many regularities and geometries seen in local patterns are needed to develop better image reconstruction algorithms. In contrast to the pixel-based regularization approach discussed above, example-based methods consists of modeling non-local pairwise interactions from training data or a library of image patches. The principle of example-based SR methods is to add a similarity constraint between the voxels of HR image and the nearest examples present in the learning database [START_REF] Baker | Limits on Super-Resolution and How to Break Them[END_REF], [START_REF] Datsenko | Example-based single document image super-resolution: a global MAP approach with outlier rejection Multidim[END_REF]. Such an approach suggests that the reconstructed image should locally look like examples existing in the learning database. The regularization term can then be defined in a general way as follows:

R(x, E) = v,k∈Ω(v) w v,k f (x(v)) -E v,k 2 (5)
where f (x(v)) is an operator on the HR image x at the voxel v, Ω(v) is a neighborhood of v, E is the learning database, E v,k is a element of E related to v and w v,k is a local weight. It is important to note that in examplebased methods the regularization term R depends on the learning database E. For instance (see [START_REF] Baker | Limits on Super-Resolution and How to Break Them[END_REF]), f (x(v)) can stand for the gradient value of the SR image x at the voxel v and E v,k would be the gradient value of the best matching pixel with respect to v in the learning database.

Although the possibility to introduce new image priors makes the example-based approach very attractive, the key point (which can be a major drawback depending on the application) is the need of a relevant learning database.

In this work, we proposed to take advantage of the existence of an HR reference image in the MR imaging context and to investigate the use of a patch-based approach without any learning database by assuming that there exist related patterns in the LR anisotropic image and an HR isotropic image of the same patient. From the point of view of example-based techniques, the HR isotropic image can be seen as a specific learning image database since information present in the HR image is used to constrain the ill-posed reconstruction problem. 3 Proposed Approach

Overview

The main idea of this work is to reconstruct an HR image using one LR image and intermodality priors from another HR image. In this context, we propose to use a non-local patch-based approach to define the regularization term in order to take into account complex spatial interactions within images. Moreover, in contrast to example-based approaches for image modeling, the proposed method is unsupervised and thus uses no image patch learning database and no computationally intensive training algorithms.

The three key points of the proposed approach are:

• a non-local regularization functional which can handle repetitive patterns in the image domain,

• the assumption that an HR reference image E x of a patient may contain relevant examples which should be used to reconstruct an HR image x from a LR image y of the same patient,

• no use of external learning database.

The principle of non-local techniques is briefly described in Section 3.2. Section 3.3 presents how the patterns of the reference image E x are used to drive the reconstruction process of x. An overview of the algorithm is provided in Section 3.4.

A non-local approach

Recently, Buades et al. [START_REF] Buades | A review of image denoising algorithms, with a new one Multiscale Modeling & Simulation[END_REF] have proposed a very efficient denoising algorithm relying on a non-local framework. Since then, this non-local strategy has been studied and applied in several image processing applications such as non-local regularization functionals in the context of inverse problems [START_REF] Rousseau | Brain Hallucination European Conference on Computer Vision, Part I[END_REF][START_REF] Mignotte | A non-local regularization strategy for image deconvolution[END_REF][START_REF] Kinderman | Deblurring and Denoising of Images by Nonlocal Functionals Multiscale Model[END_REF][START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF][START_REF] Peyré | Non-local Regularization of Inverse Problems European Conference on Computer Vision, Part III[END_REF].

Let us consider a weighted graph w that links together voxels v, k over the image domain with a weight w(v, k). This weighted graph w is a representation of non-local interactions between image elements. In [START_REF] Buades | A review of image denoising algorithms, with a new one Multiscale Modeling & Simulation[END_REF], the graph w is used for denoising purpose using a non-local neighborhood averaging strategy (called non-local means (NLM)):

∀v, d N LM (x(v)) = k∈Ω(v) w(v, k)x(k) (6)
where

d N LM (x(v)) is a denoised version of x(v).
The computation of the weighted graph is described in the Appendix.

Buades et al. have shown that, for 2D natural images, the NLM filter outperforms state-of-the-art denoising methods such as the Rudin-Osher-Fatemi Total Variation minimization scheme or the Perona-Malik Anisotropic diffusion. The NLM method tries to take advantage of the high degree of redundancy of any natural image and appears to be an unsupervised example-based denoising method.

In [START_REF] Rousseau | Brain Hallucination European Conference on Computer Vision, Part I[END_REF][START_REF] Mignotte | A non-local regularization strategy for image deconvolution[END_REF], a non-local regularization functional for inverse problems has been proposed relying on a NLM denoised version of the reconstructed image:

R(x) = v x(v) -d N LM (x(v)) 2 . ( 7
)
Such a regularization functional has the advantage over popular variational regularization functionals to handle in a better way repetitive structures and texture without incorporating strong priors on the reconstructed image.

A non-local regularization functional using intermodality priors

We propose to introduce into the regularization term a similarity constraint between local patterns of the reconstructed HR image x and the HR reference image E x . Let w Ex the weighted graph estimated on the HR reference image E x and w x the weighted graph estimated on the HR reconstructed image x.

A naive approach to use intermodality priors for SR would be to directly compute the weighted graph w Ex and to estimate the SR image using the regularization functional defined in Equation 7. In this case, the denoised version of the reconstructed HR image x would be computed using the weighted graph w = w Ex in Equation 6:

∀v, d N LM (x(v), E x ) = k∈Ω(v) w Ex (v, k)x(k) (8)
This simple strategy relies on the assumption that local patterns in the HR reference image E x could be directly used as examples to regularize the reconstructed HR image x1 . However, there are some cases where this assumption may not hold. In the context of MR imaging, images E x and y are not acquired with the same MR sequence. Multimodal MR data do not reveal the same tissue specificity; for instance, Multiple Sclerosis (MS) lesions are clearly visible in T2-weighted images but not in T1-weighted images (see Figure 3). In this case, the HR T1-weighted image may not be a relevant candidate to guide the reconstruction process. In order to handle the case of possible outliers, we propose to modify the regularization term by locally analyzing the correlation between the two graphs w x and w Ex .

For each voxel v, we define a scalar α as the correlation between the two weighted graphs w x and w Ex . To handle outliers (such as MS lesions), we propose an adaptive scheme by modifying the way to compute the denoised version of the HR reconstructed image x:

d N LM (x(v), E x ) = α(v) k∈Ω(v) w Ex (v, k)x(k) + (1 -α(v)) k∈Ω(v) w x (v, k)x(k). ( 9
)
d N LM (x(v), E x
) is now a weighted average of two denoised versions of x(v). The first term is the denoised version of x(v) computed by using the HR reference image E x and the second term is a denoised version of x(v) obtained with the current estimate of the reconstructed HR image x. If the two weighted graphs are correlated, the HR reference image E x is likely to be a relevant candidate to guide the reconstruction process and thus, α is close to 1. If the weighted graphs are uncorrelated, the presence of outliers is detected and α is close to 0. Defining d N LM (x(v), E x ) in such way allows us to choose the best examples to regularize the reconstructed HR image x.

We define then a regularization functional R(x, E x ) taking advantage of intermodality priors as follows:

R(x, E x ) = v x(v) -d N LM (x(v), E x ) 2 (10)
where d N LM (x(v), E x ) is a denoised version of x(v) using the HR reference image E x (see Equation 9). Compute the weighted graph w x 4:

Update α 5:

Estimate d N LM (x, E x ) 6:
Update x by minimizing L(x, y, H) + λR(x, E x ) 7: until convergence

Results

To explore the ability to reconstruct high resolution images of realistic typical anatomical brain structures, we applied the algorithm on MRI images of Brainweb [START_REF] Cocosco | Online Interface to a 3D MRI Simulated Brain Database Proceedings of 3-rd International Conference on Functional Mapping of the[END_REF] and real images (ADNI data).

Brainweb images

Brainweb is a simulated brain database which is often used as a gold standard for the analysis of in vivo MR data. Using the HR image provided by Brainweb, we have generated LR images with an in-plane resolution of 1mm x 1mm and a slice thickness of 3mm, using the observation model described by Equation 1. As in the Brainweb simulator, the point spread function is modeled with a boxcar function. Intensities of HR Brainweb images have been normalized between 0 and 255.

The database contains simulated brain MRI data based on two anatomical models: normal and multiple sclerosis (MS lesions have been extracted from real MRI data). We also have investigated the influence of Gaussian noise (in this case, both LR images and HR reference images are corrupted with Gaussian noise). The proposed methodology has been implemented on 4 cases: To investigate the influence of the parameters, we use the PSNR in decibels (dB) to quantify the quality of reconstructed images: In all experiments, a gradient descent method is used to optimize the cost function. The parameter λ is set to 0.01 (this value has been estimated in order to maximize the PSNR of the recovered image on Brainweb database).

P SN R(x) = 10 log 10 d 2 M SE(x) , M SE(x) = 1 |Ω| v∈Ω (x(v) -x(v)) 2 (11)
By comparing the ground truth image with the reconstructed HR brain images, the Brainweb dataset allows us to evaluate the influence of each parameter of the proposed approach.

Influence of the smoothing parameter β

Figure 4 shows the influence of the smoothing parameter β for the computation of non-local weights w N LM . For denoising purpose, Buades et al. [START_REF] Buades | A review of image denoising algorithms, with a new one Multiscale Modeling & Simulation[END_REF] show that this parameter should be set to 1. Using the Brainweb image ground truth for normal and pathological cases, we investigate different values of β. Based on this experiment, it appears that a small value of β (such as 0.1 or 0.25) may significantly decrease the quality of the reconstructed image. However, for higher values of β, the quality of the reconstructed image does not seem to be strongly affected which leads to the assumption that the proposed non-local approach is not very sensitive to this setting.

Influence of the size of the search volume and the neighborhood

Figures 5 and6 show the influence of the size of the search volume and the local neighborhood. Increasing the number of voxels in the search volume does not seem to affect the PSNR when the half-size of the search volume is greater than 3. Practically, using a small value for search volume size prevents useless computations. Concerning the patch size, it clearly appears that a patch of 3 × 3 × 3 is of optimal size with respect to the PSNR and the computation time.

Influence of the image correlation

Figure 7 shows the correlation map computed between the HR T1-weighted image and the HR T2-weighted image for the multiple sclerosis Brainweb case. This figure shows that for regions containing outliers (i.e. MS lesions), α is low. For these regions, the final image is only reconstructed using information of LR T2-weighted image. Moreover, we also report results obtained with fixed values of α in Table 1 (α = 0 or 1). When α = 0 (which corresponds to the case of SR single-frame with non-local regularization), the quality of the reconstructed image is better than using interpolation methods. However, better results are obtained using an HR reference image during the reconstruction process (adaptive α using graph correlation or α = 1). We have reported PSNR results obtained with the different methods in Table 1. Results obtained with the proposed approach compare favorably with B-spline interpolation. Figures 8 and9 show the results for non pathological and pathological MR Brainweb images for axial, coronal and sagittal views. Figure 10 is a zoom in on reconstruction results and allows a better visual comparison between interpolation method (trilinear and cubic B-spline) and the proposed SR method. It can be seen that details have been recovered and contrast between structures has been improved. Typical interpolation artefacts (see the skull for instance) are less visible on images reconstructed using the proposed method.

Quantitative comparison with interpolation methods on Brainweb images

Figure 11 shows the difference between the T2-weighted ground truth image and the reconstructed images (non pathological case). It highlights the fact that the proposed framework provides a better reconstructed image, specifically at boundaries of brain structures. As pointed out in Section 3.3, the presence of brain lesions might introduce bias in reconstructed images using intermodality priors. Figure 12 shows results in MS lesion areas. Again, it appears that the proposed method using intermodality priors compares favorably with interpolation methods. More specifically, no spurious artefacts have been detected in reconstructed images. These results tend to show that our reconstruction method is robust to outliers such as MS lesions.

In MS, lesions are usually less visible in T1-weighted images than in T2-weighted images. To verify that our algorithm does not add any reconstruction artefacts, an HR T1-weighted image has been reconstructed using an HR T2-weighted image as reference image in a pathological case. Figure 12 shows details of reconstructed images in a MS lesion area. It seems that the proposed method does not add "hallucinated" (or false) lesion. This particular point is discussed in Section 5.

In addition to PSNR values and visual assessment, an error analysis has been done in MS lesion areas. Figures 13 and14 show the histograms of error reconstruction for a reconstructed T2-weighted image (using an HR T1-weighted image as reference) and a reconstructed T1-weighted image (using an HR T2-weighted image as reference). In both cases, the standard deviation of the error reconstruction distribution is smaller for our method than for interpolation method. This tends to confirm that the proposed method allows a better reconstruction (specifically at boundaries) and that it is robust to outliers such as MS lesions.

Results on ADNI images

Data used in the preparation of this article were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (www.loni.ucla.edu/ADNI). The ADNI was launched in 2003 by the National Institute on Aging (NIA), the National Institute of Biomedical Imaging and Bioengineering (NIBIB), the Food and Drug Administration (FDA), private pharmaceutical companies and non-profit organizations, as a $60 million, 5year public-private partnership. The primary goal of ADNI has been to test whether serial magnetic resonance imaging (MRI), positron emission tomography (PET), other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of mild cognitive impairment (MCI) and Figure 15: Reconstruction results on ADNI data (the LR T2-weighted image is firstly registered on the HR T1weighted image). First row : nearest-neighbor interpolation, Second row : third order B-spline interpolation, Third row : the proposed approach. early Alzheimer's disease (AD). Determination of sensitive and specific markers of very early AD progression is intended to aid researchers and clinicians to develop new treatments and monitor their effectiveness, as well as lessen the time and cost of clinical trials. The Principle Investigator of this initiative is Michael W. Weiner, M.D., VA Medical Center and University of California -San Francisco. ADNI is the result of efforts of many co-investigators from a broad range of academic institutions and private corporations, and subjects have been recruited from over 50 sites across the U.S. and Canada. The initial goal of ADNI was to recruit 800 adults, ages 55 to 90, to participate in the research -approximately 200 cognitively normal older individuals to be followed for 3 years, 400 people with MCI to be followed for 3 years, and 200 people with early AD to be followed for 2 years." For up-to-date information see www.adni-info.org.

The image size and resolution are respectively 160 × 192 × 192 voxels and 1.2mm × 1.250mm × 1.250mm for the T1-weighted image, and 228 × 256 × 55 voxels and 0.938mm × 0.938mm × 3.000mm for the LR axial T2-weighted image. The T1-weighted image is first registered to the LR T2-weighted image using a rigid transform and then the proposed SR technique is applied. Figures 15 and16 show typical results obtained on the ADNI data using a T1-weighted image as the reference image to drive the reconstruction process of an axial T2-weighted image (the resolution of the reconstructed T2-weighted image is then : 0.938mm × 0.938mm × 0.938mm). We can observe on these Figures that for the reconstructed image using the proposed approach, the contrast of brain structures is higher than using interpolation techniques. Typical interpolation artefacts (like staircase artefact) are less visible on the reconstructed image. Moreover, it can be noticed that the non-local regularization framework allows to jointly denoise and reconstruct the T2-weighted image.

To understand how the HR T1-weighted is used to drive the reconstruction process, the correlation map is displayed in Figure 17. It appears that the parameter α has a high value near to the edges which means that the HR reference image is mainly used for brain region with high contrast (interface between white matter and grey matter for instance). Conversely, uniform and/or textured regions such as caudate nucleus or putamen are only reconstructed using information existing in the LR T2-weighted image.

Looking at Figure 15, one can argue that the proposed reconstructed framework leads to an over smooth image which may not be realistic. As Figure 18 shows, the reconstructed T2-weighted image has indeed a very similar aspect to the HR T1-weighted image denoised using the Non-Local Mean algorithm [START_REF] Buades | A review of image denoising algorithms, with a new one Multiscale Modeling & Simulation[END_REF]. Therefore, this property of the proposed algorithm appears to be beneficial since the proposed processing leads to a denoised HR T2-weighted MR image.

Discussion

The first contribution of this work concerns a simple idea for image SR which is the use of an HR reference image to improve the resolution of the LR image. As previously stated, such approach is directly related to the medical context we are interested in. However, we believe that such new image SR approach may have a substantial impact in the image processing research field. Moreover, possible further work can focus on the unification of this framework with other previously proposed approaches such as object-based interpolation and registration-based interpolation.

The second contribution is about the modeling of the SR problem introducing a non-local regularization term. We have shown that example-based approaches such as non-local means can be embedded into the reconstruction process to enhance the performance of super resolution techniques. If an HR brain image of the patient is available, LR images of the same patient can be enhanced by exploiting non-local pairwise interactions. We believe that this work is a new original way to tackle the problem of image enhancement in MR imaging by exploiting the multimodality aspect of MR data. The key point of this work is the use of a nonlocal approach to define a new regularization term in the reconstruction process. The proposed methodology is related to example-based SR methods such as the one developed by Baker et al. [START_REF] Baker | Hallucinating Faces Fourth Int. Conf. on Automatic Face and Gesture Recognition[END_REF] but in our case the learning database is reduced to one reference image. Moreover, in this work we exploit non-local interactions between voxels by defining non-local weights w (also called non-local graphs) which take into account non-local interactions in the reconstructed image and the reference image.

Experimental results show that the developed algorithm compares favorably with interpolation approaches. The two key points of the proposed approach, with respect to interpolation methods, are the use of an observation model (as in SR approaches) and the use of a reference HR image which drives the reconstruction process. Experiments on Brainweb images show that even with the presence of lesions visible only in the LR T2-weighted image, the reconstruction method we propose is able to recover such "outliers" without introducing artefacts which may come from the HR T1-weighted image (where lesions are not clearly visible). In this particular case, we would have thought that the lesions would disappear. Although our experiments tend to show that our approach is robust to such outliers, this is a crucial point that needs to be further investigated (small lesions, tumors, etc.). In general, this point raises the following question: what is the influence of an interpolation or image SR algorithm on image analysis (segmentation, detection, etc.) ?

This points out that such image enhancement technique needs to be validated using physical phantoms. Moreover, an evaluation on large database is required to prove that SR based methods can have a significant impact on medical image processing pipelines (specifically for registration, segmentation or change detection). Such validation can be application dependent and thus requires a specific evaluation framework. As it has been shown in Section 4, in addition to PSNR and visual assessment, there is a clear need to develop standard criteria to measure image quality reconstruction of SR technique. This is a key point for image enhancement technique particularly in medical imaging where image artefacts have to be avoided.

On one hand, this work shows that SR based techniques provides higher MR image quality than standard interpolation algorithms. On the other hand, interpolation techniques do not require any rigid registration step and are less time consuming (SR based methods may require several hours per reconstruction). Based on our experiments (not shown in the paper), registration does not introduce artefacts in the reconstructed images using the proposed approach. If images are not correctly registered, the correlation coefficient α tends to zero which means that the reconstruction process is only driven by the T2-weighted image.

High resolution imaging is a key point for MR brain image analysis in order to study anatomical details. Since few brain image analysis algorithms take into account an observation model such as in Eq. 1, the proposed reconstruction approach may have a substantial impact on segmentation or registration. Therefore, such model based HR image reconstruction algorithm may represent an important step towards multimodal brain analysis at fine scale.

While in this work we are only focused on brain MR image reconstruction, the proposed approach relying intermodality priors might have potential applications to other medical image modalities. Future work would involve studying a similar method for multimodal images (computerized tomography, diffusion MRI, ultrasound etc.). Depending on the intended medical application, such approach could allow to save time of data acquisition by enabling image quality improvement as a post-acquisition step. The corner stone would concern then the criterion for linking the different modalities (How to calculate the local correlation map ? What type of information can be used between the imaging modalities ?).

A Computation of weighted graphs w x and w E x wx and wE x are computed in the same way. For the sake of clarity, we describe only the computation of wx.

The weight wx(v, k) between two voxels v and k is computed as follows [START_REF] Coupé | An Optimized Blockwise Non Local Means Denoising Filter for 3D Magnetic Resonance Images[END_REF]:

wx(v, k) = 1
Zv e -P (x(v))-P (x(k)) 2 2 2N i β σ2 [START_REF] Grevera | An objective comparison of 3-D image interpolation methods[END_REF] where P (x(k)) is a 3D patch of x centered in voxel k and P (x(v)) is the 3D patch of x centered in voxel v; Zv is a constant of normalization; Ni is the number of voxels of a 3D patch. The distance between the 3D patches is the sum over voxels of patches of intensity differences using the L2 norm. The specific case concerning the central point of the patch is computed as in [START_REF] Kinderman | Deblurring and Denoising of Images by Nonlocal Functionals Multiscale Model[END_REF]. With the assumption of Gaussian noise in images, β is set to 1 (see [START_REF] Buades | A review of image denoising algorithms, with a new one Multiscale Modeling & Simulation[END_REF] for theoretical justifications) and the standard deviation of noise is estimated via pseudo-residuals ǫv as defined in [START_REF] Gasser | Residual variance and residual pattern in nonlinear regression[END_REF]. For each voxel v of x, let us define:

ǫv = r 6 7 0 @ x(v) - 1 6
X k∈Ω(v)

x(v)

1 A (13)
where Ω(v) is the 6-neighborhood at voxel v. The standard deviation of noise is computed as the least square estimator:

σ2 = 1 n X v ǫ 2 v ( 14
)
where n is the number of voxels of x.

Moreover, as suggested by Mahmoudi and Sapiro in [START_REF] Mahmoudi | Fast image and video denoising via nonlocal means of similar neighborhoods[END_REF], voxel preselection can avoid useless computation and also improve the result of denoising. The preselection of relevant voxels is based on the mean and variance of patches [START_REF] Coupé | An Optimized Blockwise Non Local Means Denoising Filter for 3D Magnetic Resonance Images[END_REF]. wx(v, k) is set to 0 if one of these conditions is not fulfilled: µ < P (x(v))

P (x(k)) < 1 µ (15)
σ 2 < var(P (x(v))) var(P (x(k)))

< 1 σ 2 (16)
with µ = 0.95 and σ 2 = 0.5.

Figure 2 :

 2 Figure 2: Comparison of single image SR framework and the proposed SR approach. Single image SR techniques aim at estimating an HR image using only one LR image as input. In the proposed strategy, the SR estimation is done using jointly the LR image and a reference HR image.

Figure 3 :

 3 Figure 3: Illustration of possible outliers. The MS lesion (in the circle) is less visible in the HR T1-weighted image than in the LR T2-weighted image. In this case, the assumption that local patterns in the HR image E x could be used as examples to regularize the reconstructed HR image x does not hold.

3. 4

 4 Overview of the algorithm Require: LR image y, starting estimate of HR image x, HR reference image E x 1: Compute the weighted graph w Ex 2: repeat 3:

1 .

 1 normal anatomical model without noise, 2. multiple sclerosis anatomical model without noise, 3. normal anatomical model with Gaussian noise (σ = 10), 4. multiple sclerosis anatomical model with Gaussian noise (σ = 10),

Figure 4 :

 4 Figure 4: Influence of the smoothing parameter β on the PSNR computed with the reconstructed image compare the Brainweb ground truth. •: non-pathological noise free Brainweb image, ♦: Multiple Sclerosis noise free Brainweb image, ×: non-pathological noisy Brainweb image, ⋆: Multiple Sclerosis noisy Brainweb image.

Figure 5 :

 5 Figure 5: Influence of the size of the search volume on the PSNR computed with the reconstructed image compare the Brainweb ground truth. •: non-pathological noise free Brainweb image, ♦: Multiple Sclerosis noise free Brainweb image, ×: non-pathological noisy Brainweb image, ⋆: Multiple Sclerosis noisy Brainweb image.

Figure 6 :

 6 Figure 6: Influence of the size of the neighborhood on the PSNR computed with the reconstructed image compare the Brainweb ground truth. •: non-pathological noise free Brainweb image, ♦: Multiple Sclerosis noise free Brainweb image, ×: non-pathological noisy Brainweb image, ⋆: Multiple Sclerosis noisy Brainweb image.

Figure 7 :

 7 Figure 7: Estimated correlation map for the multiple sclerosis Brainweb image. First row : HR T2-weighted image, Second row: HR T1-weighted image, Third row: Correlation map (values of the parameter α. A linear intensity scale between 0 and 1 is used to display α). Three arrows show MS lesion locations (axial view).

Figure 8 :

 8 Figure 8: Reconstruction results (non pathological case). From left to right: Ground truth, Trilinear interpolation, third order B-spline interpolation, the proposed approach.

Figure 9 :

 9 Figure 9: Reconstruction results in presence of MS lesions. From left to right: Ground truth, Trilinear interpolation, third order B-spline interpolation, the proposed approach.

Figure 10 :

 10 Figure 10: Zoom on reconstruction results. From left to right: Ground truth, Trilinear interpolation, third order B-spline interpolation, the proposed approach.

Figure 11 :

 11 Figure 11: Difference between ground truth image and reconstruction results (non pathological case, axial view). From left to right: Trilinear interpolation, third order B-spline interpolation, the proposed approach.The same linear intensity scale is used for the display of image differences.

Figure 12 :

 12 Figure 12: Difference between ground truth image and reconstruction results (pathological case). First row : reconstruction of a T2-weighted image, second row : reconstruction of a T1-weighted image. From left to right: Trilinear interpolation, third order B-spline interpolation, the proposed approach, MS lesion mask. The same linear intensity scale is used for the display of image differences.

Figure 13 :Figure 14 :

 1314 Figure13: Histogram of error reconstruction (difference between ground truth T2-weighted image and reconstruction results) in MS lesion mask (pathological case). Mean and standard deviation of these distributions are: trilinear(0.32 ; 97.62), cubic B-spline(2.10 ; 71.17), proposed method(0.74 ; 67.72).

Figure 16 :

 16 Figure 16: Zoom on reconstruction results on ADNI data (coronal view). A) Nearest Neighbor Interpolation, B) Trilinear Interpolation, C) third order B-spline interpolation, D) Proposed approach.

Figure 17 :

 17 Figure 17: Estimated correlation map for ADNI images. First row : HR T1-weighted image, Second row: reconstructed T2-weighted image, Third row: Correlation map (values of the parameter α. A linear intensity scale between 0 and 1 is used to display α).

Figure 18 :

 18 Figure 18: Zoom on reconstruction results on ADNI data (coronal view). A) Nearest Neighbor Interpolation, B) Proposed approach, C) T1-weighted image used as reference image, D) denoised T1-weighted image.

Table 1 :

 1 Performances of image reconstruction methods (PSNR in dB) using noise free Brainweb images. The proposed method is compared with several interpolation methods. The influence of the correlation map (value of the parameter α) on the reconstruction results is also reported.

	Method	Non pathological Image Multiple Sclerosis Image
	Nearest Neighbor Interpolation	23.74	20.89
	Trilinear Interpolation	24.65	21.61
	3rd order B-spline Interpolation	24.76	21.69
	5th order B-spline Interpolation	24.73	21.67
	Proposed method	32.20	29.48
	Proposed method with α = 0	28.30	26.05
	Proposed method with α = 1	32.09	29.16

It is important to note that in example-based SR methods, the elements of the learning database used for reconstruction are image intensities. In the proposed approach, the elements of the learning database are weighted graphs, which describe non-local interactions between image voxels.

Acknowledgment

The research leading to these results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 207667).

The author would like to thank Sylvain Faisan and Vincent Noblet from LSIIT (CNRS/University of Strasbourg) for fruitful discussions. Data collection and sharing for this project was funded by the Alzheimer's Disease Neuroimaging Initiative (ADNI; Principal Investigator: Michael Weiner; NIH grant U01 AG024904). ADNI is funded by the National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering (NIBIB), and through generous contributions from the following: Pfizer Inc., Wyeth Research, Bristol-Myers Squibb, Eli Lilly and Company, GlaxoSmithKline, Merck & Co. Inc., AstraZeneca AB, Novartis Pharmaceuticals Corporation, Alzheimer's Association, Eisai Global Clinical Development, Elan Corporation plc, Forest Laboratories, and the Institute for the Study of Aging, with participation from the U.S. Food and Drug Administration. Industry partnerships are coordinated through the Foundation for the National Institutes of Health. The grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer's Disease Cooperative Study at the University of California, San Diego. ADNI data are disseminated by the Laboratory of Neuro Imaging at the University of California, Los Angeles.

The research leading to these results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 207667). Data used in the preparation of this article were obtained from the Alzheimer[Pleaseinsert"PrerenderUnicode-˝intopreamble]s Disease Neuroimaging Initiative (ADNI) database (www.loni.ucla.edu/ADNI). As such, the investigators within the ADNI contributed to the design and implementation of ADNI and/or provided data but did not participate in analysis or writing of this report. ADNI investigators include (complete listing available at www.loni.ucla.edu/ADNI/Collaboration/ADNI Manuscript Citations.pdf)


Markus Weimer
email: weimer@acm.org

Alexandros Karatzoglou

Alex Smola
email: alexis@ci.tuwien.ac.at

Improving Maximum Margin Matrix Factorization

Collaborative filtering is a popular method for personalizing product recommendations. Maximum Margin Matrix Factorization (MMMF) has been proposed as one successful learning approach to this task and has been recently extended to structured ranking losses. In this paper we discuss a number of extensions to MMMF by introducing offset terms, item dependent regularization and a graph kernel on the recommender graph. We show equivalence between graph kernels and the recent MMMF extensions by Mnih and Salakhutdinov [1]. Experimental evaluation of the introduced extensions show improved performance over the original MMMF formulation.

Introduction

Collaborative filtering has gained much attention in the machine learning community due to its applications in electronic commerce sites such as those of Amazon, Apple and Netflix. Such sites typically offer personalized recommendations to their customers. The quality of these suggestions is crucial to the overall success, since good recommendations will increase the propensity of a purchase.

However, suggesting the right items is a highly nontrivial task: (1) There are many items to choose from. [START_REF] Hoffman | Latent semantic models for collaborative filtering[END_REF] Customers are willing to consider only a small number of recommendations (typically in the order of ten). Collaborative filtering addresses this problem by learning the suggestion function for a user from ratings provided by this and other users on items.

The known data can be thought of a sparse n×m matrix Y of rating/purchase information, where n denotes the number of users and m is the number of items. In this context, Y ij indicates the rating of item j given by user i. Typically, the rating is given on a five star scale and thus Y ∈ {0, . . . , 5} n×m , where the value 0 indicates that a user did not rate an item. In this sense, 0 is special since it does not indicate that a user dislikes an item but rather that data is missing.

Related Work A common approach to collaborative filtering is to fit a factor model to the data. For example by extracting a feature vector for each user and item in the data set such that the inner product of these features minimizes an explicit or implicit loss functional (e.g. [START_REF] Hoffman | Latent semantic models for collaborative filtering[END_REF] following a probabilistic approach).

The underlying idea behind these methods is that both user preferences and item properties can be modeled by a number of factors.

The basic idea of matrix factorization approaches is to fit the original Y matrix with a low rank approximation F . More specifically, the goal is to find such an approximation that minimizes the sum of the squared distances between the known entries in Y and their predictions in F . One possibility of doing so is by using a Singular Value Decomposition of Y and by using only a small number of the vectors obtained by this procedure. In the information retrieval community this numerical operation is commonly referred to as Latent Semantic Indexing.

Note, however, that this method does not do justice to the way Y was formed. An entry Y ij = 0 indicates that we did not observe a (user,object) pair. It does, however, not indicate that user i disliked object j. In [START_REF] Srebro | Weighted low-rank approximations[END_REF], an alternative approach is suggested which is the basis of the method described in this paper. We aim to find two matrices U and M where U ∈ R n×d and M ∈ R d×m such that F = U M with the goal to approximate the observed entries in Y rather than approximating all entries at the same time.

In general, finding a globally optimal solution of the low rank approximation problem is unrealistic: in particular the aproach proposed by [START_REF] Srebro | Weighted low-rank approximations[END_REF] for computing the a weighted factorization, which is relavant in collaborative filtering settings, requires semidefinite programming which is feasible only for hundreds, at most, thousands of terms. Departing from the goal of minimizing the rank, Maximum Margin Matrix Factorization (MMMF) aims at minimizing the Froebenius norms of U and M , resulting in a set of convex problems when taken in isolation and thus tractable by current optimization techniques. It was shown in [START_REF] Srebro | Maximum-margin matrix factorization[END_REF][START_REF] Srebro | Rank, trace-norm and max-norm[END_REF] that optimizing the Froebenius norm is a good proxy for optimizing the rank in its application to model complexity control. Similar ideas based on matrix factorization have been also proposed in [START_REF] Salakhutdinov | Probabilistic matrix factorization[END_REF][START_REF] Takács | Major components of the gravity recommendation system[END_REF].

Recently [START_REF] Weimer | Cofi rank -maximum margin matrix factorization for collaborative ranking[END_REF] proposed to extend the general MMMF framework in order to minimize structured (ranking) losses instead of the sum of squared errors on the known ratings. Key in the reasoning is that collaborative filtering is often at the heart of recommender systems. For those, only the ranking of unrated items in terms of the user preferences matter. To enable effective optimization of the structured ranking loss, a novel optimization technique [START_REF] Smola | Bundle methods for machine learning[END_REF] was used to minimize the loss in terms of the Normalized Discounted Cumulative Gain (NDCG).

Our Contribution Building upon the results outlined above, we introduce a number of extensions to the these MMMF models.

-An efficient means of computing the gradient in multiclass ordinal regression.

-A bias for movies and users to deal with heterogeneity of movies and users.

-An automatic adaptive regularization scheme which can deal with the varying number of movies per user and likewise users per movie. -A graph kernel that captures similarities between users and items in the recommender graph in the spirit of [START_REF] Salakhutdinov | Probabilistic matrix factorization[END_REF][START_REF] Basilico | Unifying collaborative and content-based filtering[END_REF]. We prove that both methods are essentially equivalent and equal to a graph kernel on the recommender graph.

For each of these extensions, we show how they can be integrated into the MMMF framework such that structured estimation as proposed in [START_REF] Weimer | Cofi rank -maximum margin matrix factorization for collaborative ranking[END_REF] is still feasible. The paper is organized as follows: Section 2 describes the general MMMF model, its generalization to structured estimation and the use of state-of-the-art optimization methods to train the model. Section 3 describes our extensions to that model. In section 4, we discuss our experimental evaluation and section 5 concludes the paper with remarks on future work.

Maximum Margin Matrix Factorization

Optimization Problem

MMMF computes a dense approximation F of the sparse matrix Y which forms the training data. The approximation is based on the modeling assumption that any particular rating of item j by user i is a linear combination of item and user features. Thus, the approximation can be written as F = U M . Here, U i * represents the feature vector for user i and M * j is the feature vector for item j. The predicted rating of item j by user i is then the inner product between these feature vectors:

F ij = U i * , M * j
Finding the appropriate matrices U and M is achieved by minimizing the regularized loss functional where the Froebenius norm ( U 2 F = tr U U ⊤ ) of the U and M matrices is used for capacity control and thus overfitting prevention. The Froebenius norm has been introduced to the MMMF framework and shown to be a proper norm on F [START_REF] Srebro | Rank, trace-norm and max-norm[END_REF]. This leads us to the following optimization problem: minimize

U,M L(F, Y ) + λ m 2 M 2 F + λ u 2 U 2 F ( 1
)
Here λ m , λ u are the regularization parameters for the M and U matrix respectively and F = U M . Moreover, L(F, Y ) is a loss measuring the discrepancy between Y and F . The optimization problem can be solved exactly by using a semi definite reformulation [START_REF] Srebro | Maximum-margin matrix factorization[END_REF]. However, this dramatically limits the size of the problem to several thousand users / movies. Instead, we exploit the fact that the problem is convex in U and M respectively when the other variables are fixed to perform subspace descent [START_REF] Rennie | Fast maximum margin matrix factoriazation for collaborative prediction[END_REF].

The Loss

Squared Loss In the original MMMF formulation, L(F, Y ) was chosen to be the sum of the squared errors [START_REF] Srebro | Maximum-margin matrix factorization[END_REF]:

L(F, Y ) = 1 2 n i=0 m j=0 S ij (F ij -Y ij ) 2 where S ij = 1 if user i rated item j 0 otherwise (2)
This loss decomposes for the non zero elements of Y and consequently it is amenable to efficient minimization by repeatedly solving a linear system of equations for each row / column of U and M separately (i.e. in parallel) -the objective function in (1) is convex quadratic in U and M respectively whenever the other term is fixed. The gradient of L(F, Y ) with respect to F can be computed efficiently, since

∂ Fij L(F, Y ) = S ij F ij -Y ij .
This means that we have

∂ F L(F, Y ) = S. * (F -Y ) (3)
where . * implies element-wise multiplication of S with F -Y . In other words, the gradient of the loss is a sparse matrix.

Non separable Loss This decomposition into losses, depending on Y ij and F ij alone, fails when dealing with structured losses that take an entire row of predictions, i.e. all predictions for a given user into account. Such losses are closer to what is needed in recommender systems, since users typically want to get good recommendations about which movies they are interested in. A fairly accurate description of which movies they hate is probably less desirable. The recent paper [START_REF] Weimer | Cofi rank -maximum margin matrix factorization for collaborative ranking[END_REF] describe an optimization procedure which is capable of dealing with such problems. In general, a non separable loss takes on the following form:

L(F, Y ) := n i=1 l(F i * , Y i * ) (4)
Gradients of L(F, Y ) decompose immediately into ∂ Fi * l(F i * , Y i * ). This allows for efficient gradient computation.

Ordinal Loss We now discuss an extension of a common ranking loss, namely the ordinal regression score, as suggested in [START_REF] Herbrich | Large margin rank boundaries for ordinal regression[END_REF]. For simplicity of notation we only study a row-wise loss l(f, y), where we assume that f := F i * and y := Y i * have already been compressed to contain only nonzero entries in Y i * with the corresponding entries of F i * having been selected accordingly. Assume that y is of length m containing m j movies of score j, that is j m j = m. For a given pair of movies (u, v) we consider them to be ranked correctly whenever y u > y v implies that also

f u > f v . A loss of 1 is incurred whenever this implication does not hold. That is, we count yu>yv C(y u , y v ) {f u ≤ f v } . ( 5
)
Here C(y u , y v ) denotes the cost of confusing a movie with score y u with one of score y v . Since there are n = 1 2 m 2j m 2 j terms in the sum we need to renormalize the error by n in order to render losses among different users comparable. Moreover, we need to impose a soft-margin loss on the comparator Algorithm 1.1. (l, g) = l(f, y, C) input Vectors f and y, score matrix C output Loss l and gradient g initialize l = 0 (loss) and g = 0 (gradient) for i = 1 to n do bi = 0 (lower counter) and ui = mi (upper counter)

end for Let c = [f -1 2 , f + 1 2 ] ∈ R 2m Rescale C ← 2C/(m 2 -u 2
2 ) index = Argsort(c) (find overlaps between pairs) for i = 1 to 2m do j = index(i) mod m and z = yj if index(i) ≤ m (from the first half) then for k = 1 to z -1 (we should be better than those) do l ← l -C(k, z)u k cj and gj ← gj -C(k, z)u k end for bz ← bz + 1 (there are now bz + 1 elements below us) else for (they should be better than us

) k = z + 1 to n do r ← r + C(z, k)b k cj+m and gj ← gj + C(z, k)b k end for uz ← uz -1 (
we've just seen one more term from above) end if end for {f u ≤ f v } to obtain a convex differentiable loss. This yields the loss

l(f, y) = 2 m 2 - j m 2 j -1 yu>yv C(y u , y v ) max(0, 1 -f u + f v ). ( 6
)
The gradient ∂ f l(f, y) can be computed in a straightforward fashion via

∂ f l(f, y) = -2 yu>yv C(y u , y v ) m 2 -j m 2 j (7)
In general, computing losses using preferences such as ( 6) is an O(m 2 ) operation. However, we may extend the reasoning of [START_REF] Joachims | Training linear SVMs in linear time[END_REF] to more than binary scores to obtain an O(m log m) algorithm instead. Algorithm 1.1 relies on sorting f before taking sums. It uses the decomposition of the soft margin loss via max(0, 1 -

f u + f v ) = max(0, (f v + 0.5) -(f u -0.5)) = max(0, c v+m -c v )
where c = [f -0.5, f + 0.5] to disentangle upper and lower bounds. It then traverses the sorted list of c to check for how many terms an upper or lower bound is violated by means of auxiliary counters b and u.

Fig. 1. A convex function (solid) is bounded from below by Taylor approximations of first order (dashed). Adding more terms improves the bound.

Optimization

Although (1) is not jointly convex in U and M , it is still convex in U if M is kept fixed and convex in M if U is kept fixed. We thus resort to alternating subspace descent as proposed by [START_REF] Rennie | Fast maximum margin matrix factoriazation for collaborative prediction[END_REF] by keeping U fixed and minimizing over M and repeating the process for M with U fixed. We have the following procedure: repeat For fixed M minimize (1) with respect to U . For fixed U minimize (1) with respect to M . until no more progress is made or a maximum iteration count is reached. While this approach does not ensure that a global minimum is reached, it proves to be rather efficient and scalable for problems of up to 10 8 nonzero entries in Y (Netflix) [START_REF] Weimer | Cofi rank -maximum margin matrix factorization for collaborative ranking[END_REF].

Our implementation uses a bundle method solver for the two optimization steps over U and M . Recently, bundle methods have been introduced with promising results for optimizing regularized risk functions in supervised machine learning [START_REF] Smola | Bundle methods for machine learning[END_REF]. Bundle methods are especially suited to structured output estimation problems as they require relatively few function and gradient computations which can be costly in the structured estimation setting [START_REF] Tsochantaridis | Large margin methods for structured and interdependent output variables[END_REF], while they can be also used to minimize non-smooth loss functions where only subdifferentials are available.

The key idea behind bundle methods is to compute successively improving linear lower bounds of an objective function through first order Taylor approximations as shown in Figure 1. Several lower bounds from previous iterations are bundled in order to gain more information on the global behavior of the function. The minimum of these lower bounds is then used as a new location where to compute the next approximation, which leads to increasingly tighter bounds and convergence.

The main computational cost in using the bundle method solver is the computation of the gradients with respect to M and U . Using the chain rule yields This computation is easily parallelizable, since terms in ∂ F L(F, Y ) with respect to each user can be computed separately.

∂ M L(F, Y ) = U ⊤ ∂ F L(F, Y ) and ∂ U L(F, Y ) = [∂ F L(F, Y )] ⊤ M. (8)
We assume that the loss decomposes per user. Thus, we can optimize U by optimizing each row of U on its own as shown in algorithm 1.2. Note that we effectively construct and solve a regularized risk model for each user for a dense data matrix X and parameters w. On the other hand, when minimizing with respect to M we need to deal with the entire loss jointly. The main issue to solve is to compute the loss and its gradient with respect to M . Algorithm 1.3 shows an efficient way to compute the gradient which decomposes again for all users besides the final multiplication.

Extensions

So far we described a number of practical implementation details and extensions in the context of the loss function. We now discuss extensions of the model to take prior knowledge about the function class into account.

Offset Individual users may have different standards when it comes to rating movies. For instance, some users may rarely award a 5 while others are quite generous with it. Likewise, movies have an inherent quality bias. For instance, 'Plan 9 from Outer Space' will probably not garner high ratings with any movie buff while other movies may prove universally popular. This can be taken into account by means of an offset per movie. This can be incorporated via

F ij = U i * , M j * + u i + m j . (9)
Here u and m are bias vectors for movies and users alike. In practice, we simply extend the dimensionalities of U and M by one for each bias while pinning the corresponding coordinate of the other matrix to assume the value of 1. In this form no algorithmic modification for the U and M optimization is needed. The computational cost of this extension is near zero, as the feature vectors of the convex optimization problem are extended by only one or two dimensions.

Please note that this offset is different from a simple normalization of the input data: It is learned for each user and for each movie, while a normalization cannot cater for both appropriately at the same time.

Adaptive regularization Fixing a single regularization parameter for movies and users respectively is not a very effective choice. For instance, estimation for a user who has seen very few movies will likely suffer from overfitting unless heavily regularized. On the other hand, for users who have rated lots of movies, we should be more gentle in terms of regularization. Likewise, the number of ratings per movie varies widely and the regularization should take this into account.

Those issues can be dealt with by sample-size adaptive regularization for both movies and users. Denote by D u and D m diagonal matrices corresponding to movies and users. Setting

D u ii = n -α i and D M jj = m -α j
where n i denotes the number of movies user i has rated and m j denotes the number of users which have rated movie j, we obtain a sample size dependent regularizer as follows: minimize

U,M L(U M, Y ) + λ m 2 tr M ⊤ D m M + λ u 2 tr U ⊤ D u U.
In our experiments we found that α = 0.5 provides best generalization performance. This is equivalent to the regularization scales provided in a maximum a posteriori setting where the log-prior is fixed whereas the evidence scales linearly with the number of observations. As the computation for this scaling can be done in advance, the computational cost of the adaptive regularizer is not important when compared to the overall runtime. The needed statistics can even be pre-computed and reused in many experiments.

Graph kernels So far we ignored a crucial piece of information, namely the fact that the ratings themselves are not random. For instance, knowing that a user rated 'Die Hard', 'Die Hard 2', and 'Top Gun' makes it likely that this user is interested in action movies. This information can be gained without even looking at the score matrix Y . One would expect that we should be able to take advantage of this structural information in addition to the actual scores.

One possibility, proposed by [START_REF] Basilico | Unifying collaborative and content-based filtering[END_REF] is to use the inner product between the movies two users rent as a kernel for comparing two different users. Denote by S ij = {C ij > 0}. In this case they define the kernel between users i and i ′ to be S i * , S i ′ * . It is well known that such a model is equivalent to using a linear model with user-features given by S. We can use this to improve the user matrix U to U + SA for a suitably chosen feature matrix A.

Independently, [START_REF] Salakhutdinov | Probabilistic matrix factorization[END_REF] recently developed a related line of thought by assuming that the user matrix is given by U + SA, where U and A are normally distributed and S is a row-normalized version of S, that is Si * = S i * -1 1 S i * . While their optimization strategy is very different (they use Markov Chain Monte Carlo sampling), it should already be clear at this point that the outcome is very similar to that of [START_REF] Basilico | Unifying collaborative and content-based filtering[END_REF].

We now show that both approaches, which are approximately equivalent (barring the normalization of S to S), are also equivalent to the use of graph kernels on the bipartite ranking graph defined between users and movies. For this purpose we require the following lemma: Lemma 1. Denote by f : R n → R some function and let A ∈ R n×d . Moreover, let U ∈ R n×d and S ∈ R n×m . Then the following problems are equivalent:

minimize U,A f (V ) + U 2 + A 2 (10) minimize V f (V ) + U ⊤ (1 + SS ⊤ ) -1 U (11)
Proof. Denote by (U * , A * ) the optimal solution of [START_REF] Rennie | Fast maximum margin matrix factoriazation for collaborative prediction[END_REF]. Clearly in this case for

V := U * + SA * the optimization problem minimize A f (V ) + V -SA 2 + A 2 (12)
needs to have A * as its solution. What remains is to express A as a function of V and to show that in this case ( 10) and ( 11) are equivalent. Taking derivatives of [START_REF] Herbrich | Large margin rank boundaries for ordinal regression[END_REF] with respect to A yields

∂ A V -SA 2 + A 2 = 2S ⊤ (SA -V ) + 2A.
Hence the gradient of the objective function vanishes for A * = (1+S ⊤ S) -1 S ⊤ V . Plugging this back into (11) yields the objective function

f (V ) + 1 -S(1 + S ⊤ S) -1 S ⊤ V 2 + (1 + S ⊤ S) -1 S ⊤ V 2 = f (V ) + (1 + SS ⊤ ) -1 V 2 + S ⊤ (1 + SS ⊤ ) -1 V 2 = f (V ) + V ⊤ (1 + SS ⊤ ) -1 V
Here we used the Sherman-Morrison-Woodbury identity to transform the second term in the second line. The third term follows from the fact that left and right singular vectors associated with S constitute the eigenvectors of (1 + SS ⊤ ) and (1 + S ⊤ S) respectively. Hence we may "push" S to the left in the third term. The last equality follows by direct calculation. *

This lemma shows that the parameterization U + SA (or U + SA respectively) is equivalent to using a kernel (1 + SS ⊤ ) -1 as regularization. The latter is well known as the inverse Laplacian kernel, since SS ⊤ encodes the undirected graph obtained by connecting all users which watched the same movie. The connection strength in SS ⊤ denotes the number of movies both users shared.

The net result of this reparameterization is that ( 10) is a computationally more efficient way of dealing with such symmetries rather than computing the inverse of (1 + SS ⊤ ). Since we may have millions of users the latter would be computationally infeasible.

Note also the connection to the spectral theory of graphs [START_REF] Smola | Kernels and regularization on graphs[END_REF]: the eigenvalues and eigenvectors of 1 + S S⊤ are close to those of the bipartite graph Laplacian, which can be used for clustering between movies and users respectively. This means that for similar users and movies we end up using similar parameters respectively.

Experiments

To evaluate the the extension described in section 3 we performed experiments on the eachmovie and movielens data sets. The data sets are summarized in table 1.

Data set Users Movies Ratings EachMovie 61265 1623 2811717 MovieLens 983 1682 100000

Table 1. Data set statistics

In order to compare the results to those in the literature, we evaluate our system using the Normalized Discounted Cumulative Gain (NDCG) measure:

DCG(Y i * , π)@k = m j=0 2 Y iπ[j] -1 log 2 (j + 1) , N DCG(Y i * , π)@k = DCG(Y i * , π)@k DCG(Y i * , π s )@k (13)
The permutation π is computed as the argsort of the predicted values: π = argsort(F i * ). The perfect permutation π s is the argsort of the true ratings given by the user: π s = argsort(Y i * ). A NDCG of 1.0 indicates that the model sorts the movies in the same order as the user. The parameter k is a cut-off beyond which the actual ranking does no longer matter. This follows the intuition that typical recommender systems can only present a limited amount of items to the user. In all our experiments, we evaluated using NDCG@10.

Following [START_REF] Yu | Collaborative ordinal regression[END_REF] and [START_REF] Weimer | Cofi rank -maximum margin matrix factorization for collaborative ranking[END_REF] , we distinguish two different evaluation scenarios: strong and weak generalization.

Weak generalization:

The train set is built by sampling 10,20 or 50 movies randomly from the seen movies of each user. The ranking performance is than evaluated using the NDCG@10 score on the remaining movies.

Strong generalization: The model is evaluated on users that were not present at training time. We follow the procedure described in [START_REF] Yu | Collaborative ordinal regression[END_REF]: Movies with less than 50 ratings are discarded. The 100 users with the most rated movies are selected as the test set and the methods are trained on the remaining users. In evaluation, 10, 20 or 50 ratings from those of the 100 test users are selected. For those ratings, the user training procedure is applied to optimize U . M is kept fixed in this process to the values obtained during training. The remaining ratings are tested using the same procedure as for the weak generalization.

It is important to note that in both cases, the system is evaluated on the vast majority of movies. This mimics the true situation of a recommender system which always has more unrated items to recommend than already rated items at its disposal. However, these evaluation schemes cannot show the usefulness of the per user offset, as the number of items in the train set is fixed for all users. To evaluate the influence of this extension, we performed experiments following the weak evaluation scheme with exchanged train and test sets: The system is tested on 10, 20 or 50 items per user and trained on the remaining items.

In all experiments, the regularization parameters were fixed and not formally tuned. The dimension of U and M is fixed to d = 10 in all experiments. We did not observe significant performance drop when compared to a value of d = 100 in [START_REF] Weimer | Cofi rank -maximum margin matrix factorization for collaborative ranking[END_REF] or d = 30 in [START_REF] Salakhutdinov | Probabilistic matrix factorization[END_REF]. Additionally, we present results only for the least squares regression loss, as the emphasis is on the evaluation of the extensions to the original model. All experiments were performed ten times with different random draws of the train and test set from the data set. In total, we report results on 960 experiments.

Results

Weak generalization Table 2 shows the results for the weak generalization setting. Overall we can see that adding the offset term does improve the performance while enabling just the Graph kernel does not seem to yield significant gains in performance. This can be attributed to the fact that the impact of the graph kernel will be most profound in cases where users have rated relatively few movies compared to the average user. An additional factor to take into account is that due to computational time constraints we did not adjust the regularization parameter for each configuration. This is particularly significant in the graph kernel phase where we optimize over an additional set of parameters. Nevertheless we observe that in some cases the graph kernel seems to bring significant performance particularly when combined with the offset term. This leads us to believe that with proper parameter tuning, the graph kernel can add significant performance increases.

Inverted weak generalization

The results of this set of experiments table 3, again confirm our observations in the previous setting. The offset term along often Method N=10 N=20 N=50 EachMovie Plain 0.625 ± 0.000 0.639 ± 0.000 0.641 ± 0.000 Offset 0.646 ± 0.000 0.653 ± 0.000 0.647 ± 0.000 GraphKernel 0.583 ± 0.000 0.585 ± 0.000 0.590 ± 0.001 OffsetGK 0.576 ± 0.000 0.597 ± 0.000 0.580 ± 0.001 MovieLens Plain 0.657 ± 0.000 0.658 ± 0.000 0.686 ± 0.000 Offset 0.678 ± 0.000 0.680 ± 0.000 0.701 ± 0.000 GraphKernel 0.624 ± 0.001 0.644 ± 0.000 0.682 ± 0.000 OffsetGK 0.670 ± 0.001 0.681 ± 0.000 0.682 ± 0.000 Table 2. The NGDC@10 accuracy over ten runs and the standard deviation for the weak generalization evaluation.

combined with the graph kernel or the adaptive regularization brings significant performance gains. Again we have to note that the overall regularization parameters were not tuned for the individual configurations.

Method N=10 N=20 N=50 EachMovie Plain 0.859 ± 0.000 0.731 ± 0.000 0.627 ± 0.000 Offset 0.859 ± 0.000 0.734 ± 0.000 0.631 ± 0.000 GraphKernel 0.837 ± 0.000 0.693 ± 0.000 0.553 ± 0.001 AdaReg 0.858 ± 0.000 0.729 ± 0.000 0.635 ± 0.000 AdaRegGK 0.832 ± 0.000 0.692 ± 0.000 0.578 ± 0.000 OffsetGK 0.832 ± 0.000 0.689 ± 0.000 0.587 ± 0.001 OffsetAdareg 0.859 ± 0.000 0.728 ± 0.000 0.637 ± 0.000 All 0.836 ± 0.000 0.702 ± 0.000 0.585 ± 0.000 MovieLens Plain 0.875 ± 0.000 0.750 ± 0.000 0.673 ± 0.000 Offset 0.886 ± 0.000 0.764 ± 0.000 0.703 ± 0.001 GraphKernel 0.845 ± 0.000 0.720 ± 0.001 0.667 ± 0.000 AdaReg 0.873 ± 0.000 0.736 ± 0.000 0.652 ± 0.001 AdaRegGK 0.835 ± 0.000 0.694 ± 0.001 0.645 ± 0.001 OffsetGK 0.882 ± 0.000 0.773 ± 0.000 0.703 ± 0.000 OffsetAdareg 0.874 ± 0.000 0.750 ± 0.000 0.681 ± 0.000 All 0.869 ± 0.000 0.730 ± 0.002 0.645 ± 0.005 Table 3. The NGDC@10 accuracy over ten runs and the standard deviation for the inverted weak generalization evaluation.

Strong generalization For the strong generalization setting, we compare our results to Gaussian Process Ordinal Regression (GPOR) [START_REF] Chu | Gaussian processes for ordinal regression[END_REF] Gaussian Process Regression (GPR), the collaborative extensions (CPR, CGPOR) as well as the original MMMF implementation [START_REF] Yu | Collaborative ordinal regression[END_REF], [START_REF] Weimer | Cofi rank -maximum margin matrix factorization for collaborative ranking[END_REF]. Table 4 shows our results compared to the ones from [START_REF] Yu | Collaborative ordinal regression[END_REF].

For the Movielens data, our system with the offset and graph kernel extensions outperforms the other systems we compare to. Additionally, the system with both extensions performs consistently better than the ones with only one extension. On the Eachmovie data, our system performs the best with the offset parameters enabled. it appears that that the graph kernel in the Eachmovie data set does not improve the performance but again this can be attributed to a poor choice of the regularization parameters for this data set.

On both data sets, the results are very good in the strong generalization setting. This is not as surprising as it may seem at first: In the strong generalization phase, our system solves a convex problem for each test user in order to learn the right U i * for that user. If the system learned reasonable features for the items in the weak generalization phase of this evaluation, good results for the strong evaluation phase are to be expected. We believe that this is an important benefit of our method in many applications, as it allows for fast accurate predictions for new users without the need to retrain the whole system. 4. The NGDC@10 accuracy over ten runs and the standard deviation for the strong generalization evaluation.

The variance over the ten runs on different data in all experiments is surprisingly low, especially given the fact that we are optimizing a non convex function. The same is true for the variance on the objective function. The low variance may mean that we always reach the same local minimum or that this minimum is indeed a global one.

Overall, the experiments show that the performance of the basic model can be significantly increased by the extensions proposed in this paper. Please note that we did not optimize the regularization parameters, which might improve the performance even further. Additional improvements are to be expected when applying other loss functions like ordinal regression or the ranking losses as described in [START_REF] Weimer | Cofi rank -maximum margin matrix factorization for collaborative ranking[END_REF]. Both have shown to yield better ranking performance, yet we could not evaluate them here as each loss function would have added another 960 experiments.

Conclusion

In this paper, we have shown several extensions to the original MMMF model that add up to recent advances in the optimization procedure and the losses. We introduced offset terms, item dependent regularization and a graph kernel on the recommender graph. We also showed that recent extensions to MMMF [START_REF] Salakhutdinov | Probabilistic matrix factorization[END_REF] as well as well known approaches [START_REF] Basilico | Unifying collaborative and content-based filtering[END_REF] are both instances of our graph kernel formulation.

The extensions have been introduced in a way that preserves recent extensions of MMMF to structured loss [START_REF] Weimer | Cofi rank -maximum margin matrix factorization for collaborative ranking[END_REF]. Additionally, this still allows us to use state of the art optimizers based on bundle methods which have recently been proposed for the regularized risk minimization problem [START_REF] Smola | Bundle methods for machine learning[END_REF].

On all evaluated data sets in all evaluation settings, one combination of the proposed extensions yielded significantly improved results, even though we did not tune the parameters of the model. Thus, even better results are to be expected in real world applications of this method.

The software developed to evaluate the methods described in this paper will be available on http://www.cofirank.org.

Algorithm 1 . 2 . 3 .

 123 Optimization over U input Matrix U and M , data Y output Matrix U for i = 1 to n do Select idx as the nonzero set of Yi * Initialize w = Ui * Ui,idx = argmin w l(wMidx, * , Yi,idx) + λu 2 Computation of ∂M L input Matrix U and M , data Y output ∂M L = D ⊤ M for i = 1 to n do Update w ← Ui * Find index ind where Yi * = 0 X ← M [ind, :] Update Di * ← ∂F L(wX, Yi * [ind]) end for∂M L = D ⊤ M

Acknowledgements

Markus Weimer has been funded under Grant 1223 by the German Science Foundation (DFG). Alexandros Karatzoglou was supported by a grant of the ANR -CADI project. We gratefully acknowledge support by the Frankfurt Center for Scientific Computing in running our experiments.