Florentin Millour
email: fmillour@oca.eu

Olivier Chesneau

Anthony Meilland

Nicolas Nardetto

Optical interferometry and adaptive optics of bright transients

Introduction

Bright optical transients (i.e. transients typically visible with the naked eye) are populated mainly by novae eruptions plus a few supernovae (among which the SN1987a event). Indeed, usually one bright nova happen every two year, either in the North or the South hemisphere (see Fig. 1). It occurs that current interferometers have matching sensitivities, with typically visible or infrared limiting magnitudes in the range 5-7. The temporal development of the fireball, followed by a dust formation phase or the appearance of many coronal lines can be studied with the VLTI. The detailed geometry of the first phases of novae in outburst remains virtually unexplored. This paper summarizes the work which has been done to date using mainly the Very Large Telescope Interferometer.

We invite the reader to have a look to the extensive review on the topic by Chesneau & Banerjee (2012) for a complete description of the science on transients that can be achieved with optical/infrared long-baseline interferometers. We give a short summary of the content of this paper in the next section.

2 Why observing novae with optical interferometers?

Optical interferometers represent a breakthrough in terms of spatial resolution, that can provide crucial information related to the nova phenomenon. All targets in the 3-5 kpc range can be potentially resolved by current interferometers (CHARA, VLTI, NPOI).

The VLT Interferometer can provide measurements of the angular diameters of the nova ejecta in continuum and lines, from the near-IR to the mid-IR in the very first moments of the outburst. The primary outcome of these observations is a direct 1 estimate of the expansion parallax, thus the distance to novae. Of importance is the possibility to spatially and spectrally resolve different near-IR emission lines to estimate the physical conditions throughout the wind and the ejecta. For those objectives, the medium spectral resolution of the near-IR instrument AMBER is an asset. If the nova appears to form dust (CO novae ), an in-depth study of the dust forming regions can be carried out with the MIDI instrument. Using a set of flexible observing runs, we shall follow the outburst from the first days up to several months.

Nova as a spherical fireball

Up to now, the program of observations has focused on novae with magnitudes reachable by the VLTI (South declination & K magnitude ≥ 7). Past observations and theoretical work on the nova phenomenon have provided a substantial knowledge about the physical nature of these binary systems and the outburst. However, these investigations are naturally limited by the difficulty of estimating the distance, which is usually inferred indirectly and with large errors. Spherical symmetry is a basic tenant adopted in the derivation of relationships that link the non-spatially resolved photometric and spectroscopic observations to the physical parameters of the system (Gehrz et al. 1998, Table 3). Spherical symmetry is implicitly assumed when the uv coverage is not sufficient to perform a better analysis. This was the case for the Nova V1280 Sco (Chesneau et al. 2008) which was also observed exclusively using 2 telescope recombination. Two years after, observations with the AO system NACO mounted at the UT4 telescope revealed an impressive dusty bipolar nebula (Chesneau et al. 2012, see Fig. 2).

Figure 2: 2010 NACO K band image after a PSF subtraction revealing the impressive bipolar nebula. Mid-IR images also show that there is no dust emission in the equatorial plane.

A Bipolar fireball from the first blink

An interferometer is mainly sensitive to the angular size of a nova in its early stages. Measuring the size of the fireball in different orientations on-sky allows us to infer the axis ratio and orientation of an individual nova shell. This is relatively easy to obtain for a 3-6-telescopes interferometer. This was the case for the outburst of the recurrent nova RS Oph (Chesneau et al. 2007).

The highly collimated outflow from the RS Ophiuchi has been imaged by the HST (Bode et al. 2007) and in the radio (O'Brien et al. 2006). The AMBER observations showed that the jet was already existing 5.5 days after the discovery, and provided a unique view of radial-velocities which could afterwards complement the expansion rates derived by the HST and radio images (Ribeiro et al. 2009).

The signature of a bipolar jet in interferometric data is now well identified provided that the emission lines are spectrally resolved (R∼1500). The famous nova T Pyx exhibited a spherical appearance in broadband PIONIER data, but the signature of a bipolar kinematics was clearly detected in our spectrally resolved AMBER data (Chesneau et al. 2011). The numerous peculiarities of the T Pyx eruptions can be explained in the frame of recurrent nearly face-on eruptions that launch fast material in the line-of-sight and slow material perpendicular to it, building up the slow expansion shell imaged by the HST.

Intermediate Luminosity Optical Transients

Intermediate-Luminosity Optical Transients (ILOTs), are eruptive stars with peak luminosity between those of novae and supernovae that have been also called Red Novae or red Transient. The powering processes and whether they are due to binary interaction or are formed through single star evolution are debated. High angular resolution techniques can play a role by tracking bipolarity and the formation of disks. Furthermore, one can also study the remaining central star when a merger is highly suspected, for instance by detecting the deformation due to very high rotational rate.

Sakurai's object

In 1996, Sakurai's object (V4334 Sgr) suddenly brightened in the center of a faint Planetary Nebula (PN). This very rare event was interpreted as being the reignition of a hot white dwarf that caused a rapid evolution back to the cool giant phase. From 1998 on, a copious amount of dust has formed continuously, screening out the star that remained embedded in this expanding high optical-depth envelope. Mid-IR interferometry performed in 2008 with the MIDI/VLTI instrument discovered a unexpectedly compact (30 40 milli-arc-second, 105 140 AU assuming a distance of 3.5 kpc), highly inclined, dust disk (Chesneau et al. 2009). The major axis of the disk is aligned with an asymmetry seen in the old PN. This implies that the mechanism responsible for shaping the dust envelope surrounding Sakurai's object was already at work when the old PN formed, a strong argument for binary interaction.

V 838 Mon

V838 Monocerotis erupted in 2002, brightening by 9 magnitudes in a series of outbursts, and eventually developing a spectacular light echo. A very red star emerged surrounded by copious amount of new dust that condensed from the expanding ejecta of the outbursts. V838 Mon is the close-by archetype of the ILOT sources which are triggering very active research currently. MIDI/VLTI observations obtained over the last few months showed that the dust resides in the form of a flattened structure (∼15-50 mas from 8 to 13 µm), i.e. a 90x300 AU flattened structure for a distance of 6.2 kpc (Sparks et al. 2008). The modelling of this extended structure is in progress but it is incomplete without a much better knowledge of the central source which is seen as a very cool M-L type super-giant. AMBER observations were also obtained in 2013 to measure the size of the central source, its shape (since it has potential to be a fast rotator) and study the cool photosphere / dusty disk transition. The AMBER data were essentially acquired with small baselines (≤ 30 m) and are quite noisy (σV 2 ≈ 0.05), leading to an angular size of 3 mas, but with a large uncertainty. We can basically only set an upper limit to the diameter of the HK-bands object, of 4.7 mas. This would make the HK-bands object smaller than 30 AU. This size (and shape) difference between the HK-bands and N band is striking and reminiscent of super-giant stars dusty disks (Millour et al. 2011;Wheelwright et al. 2012). Further AMBER observations would enable us to pinpoint more precise properties of this intriguing object.

Conclusion and prospectives

We presented here a few results from the VLTI campaigns on Novae and ILOTs. These campaigns present a challenge in terms of scheduling and observatory response but provide unique insights on the early processes at stake when a nova explode. With the upcoming infrared instruments like MATISSE [START_REF] Lopez | Society of Photo-Optical Instrumentation Engineers[END_REF] or GRAVITY (Eisenhauer et al. 2005), getting a finer idea of the geometry will be much faster than with current instruments. The development of visible interferometric instruments at the CHARA array or at the VLTI would also be an asset to get the sharpest multiwavelength picture of these objects a few days after outburst before the advent of the ELT in the mid-2030s which will enable direct snapshot pictures of the fireball at 10's-milli-arc-seconds resolution for a much larger number of novae.

Figure 1 :

 1 Figure 1: Typical frequency of bright novae from 1975 to 2013. About one nova every two years occurs, either in the North or South hemisphere 2 . V1663 Aql & V458 Vul do not appear in this plot as too faint in V (but bright in K). The average rate per year is plotted as a dashed line. The arrows show Novae observed with optical interferometers.

data from http://www.cbat.eps.harvard.edu/nova_list.html, complemented by data from http://asd.gsfc.nasa.gov/Koji.Mukai/novae/novae.html . Some peak magnitudes were corrected using AAVSO data.

2 data from http://www.cbat.eps.harvard.edu/nova_list.html, complemented by data from http://asd.gsfc.nasa.gov/Koji.Mukai/novae/novae.html . Some peak magnitudes were corrected using AAVSO data.


Adrien Todeschini
email: adrien.todeschini@inria.fr

Franc ¸ois Caron

Marie Chavent
email: marie.chavent@u-bordeaux2.fr

Probabilistic Low-Rank Matrix Completion with Adaptive Spectral Regularization Algorithms

We propose a novel class of algorithms for low rank matrix completion. Our approach builds on novel penalty functions on the singular values of the low rank matrix. By exploiting a mixture model representation of this penalty, we show that a suitably chosen set of latent variables enables to derive an Expectation-Maximization algorithm to obtain a Maximum A Posteriori estimate of the completed low rank matrix. The resulting algorithm is an iterative soft-thresholded algorithm which iteratively adapts the shrinkage coefficients associated to the singular values. The algorithm is simple to implement and can scale to large matrices. We provide numerical comparisons between our approach and recent alternatives showing the interest of the proposed approach for low rank matrix completion.

Introduction

Matrix completion has attracted a lot of attention over the past few years. The objective is to "complete" a matrix of potentially large dimension based on a small (and potentially noisy) subset of its entries [START_REF] Srebro | Maximum-Margin Matrix Factorization[END_REF][START_REF] Candès | Exact matrix completion via convex optimization[END_REF][START_REF] Candès | Matrix completion with noise[END_REF]. One popular application is to build automatic recommender systems, where the rows correspond to users, the columns to items and entries may be ratings or binary (like/dislike). The objective is then to predict user preferences from a subset of the entries.

In many cases, it is reasonable to assume that the unknown m × n matrix Z can be approximated by a matrix of low rank Z ≃ AB T where A and B are respectively of size m × k and n × k, with k ≪ min(m, n). In the recommender system application, the low rank assumption is sensible as it is commonly believed that only a few factors contribute to user's preferences. The low rank structure thus implies some sort of collaboration between the different users/items [START_REF] Mazumder | Spectral regularization algorithms for learning large incomplete matrices[END_REF].

We typically observe a noisy version X ij of some entries (i, j) ∈ Ω where Ω ⊂ {1, . . . , m} × {1, . . . , n}. For (i, j) ∈ Ω

X ij = Z ij + ε ij , ε ij iid ∼ N (0, σ 2 ) (1)
where σ 2 > 0 and N (µ, σ 2 ) is the normal distribution of mean µ and variance σ 2 . Low rank matrix completion can be adressed by solving the following optimization problem minimize

Z 1 2σ 2 (i,j)∈Ω (X ij -Z ij ) 2 + λ rank(Z) (2)
where λ > 0 is some regularization parameter. For general subsets Ω, the optimization problem ( 2) is computationally hard and many authors have advocated the use of a convex relaxation of ( 2) [START_REF] Fazel | Matrix rank minimization with applications[END_REF][START_REF] Candès | Enhancing sparsity by reweighted l1 minimization[END_REF][START_REF] Mazumder | Spectral regularization algorithms for learning large incomplete matrices[END_REF], yielding the following convex optimization problem minimize

Z 1 2σ 2 (i,j)∈Ω (X ij -Z ij ) 2 + λ Z * (3)
where Z * is the nuclear norm of Z, or the sum of the singular values of Z. [START_REF] Mazumder | Spectral regularization algorithms for learning large incomplete matrices[END_REF] proposed an iterative algorithm, called Soft-Impute, for solving the nuclear norm regularized minimization [START_REF] Candès | Matrix completion with noise[END_REF].

In this paper, we show that the solution to the objective function (3) can be interpreted as a Maximum A Posteriori (MAP) estimate when assuming that the singular values of Z are independently and identically drawn (iid) from an exponential distribution with rate λ. Using this Bayesian interpretation, we propose alternative concave penalties to the nuclear norm, obtained by considering that the singular values are iid from a mixture of exponential distributions. We show that this class of penalties bridges the gap between the nuclear norm and the rank penalty, and that a simple Expectation-Maximization (EM) algorithm can be derived to obtain MAP estimates. The resulting algorithm iteratively adapts the shrinkage coefficients associated to the singular values. It can be seen as the equivalent for matrices of reweighted ℓ 1 algorithms [START_REF] Candès | Enhancing sparsity by reweighted l1 minimization[END_REF] for multivariate linear regression. Interestingly, we show that the Soft-Impute algorithm of [START_REF] Mazumder | Spectral regularization algorithms for learning large incomplete matrices[END_REF] is obtained as a particular case. We also discuss the extension of our algorithms to binary matrices, building on the same seed of ideas, in the supplementary material. Finally, we provide some empirical evidence of the interest of the proposed approach on simulated and real data.

Complete matrix X

Consider first that we observe the complete matrix X of size m × n. Let r = min(m, n). We consider the following convex optimization problem minimize

Z 1 2σ 2 X -Z 2 F + λ Z * (4)
where • F is the Frobenius norm. The solution to Eq. ( 4) in the complete case is a soft-thresholded singular value decomposition (SVD) of X [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF][START_REF] Mazumder | Spectral regularization algorithms for learning large incomplete matrices[END_REF], i.e. Z = S λσ 2 (X) where S λ (X) = U D λ V T with D λ = diag(( d 1λ) + , . . . , ( d rλ) + ) and t + = max(t, 0). X = U D V T is the singular value decomposition of X with D = diag( d 1 , . . . , d r ).

The solution Z to the optimization problem (4) can be interpreted as the Maximum A Posteriori estimate under the likelihood (1) and prior p(Z) ∝ exp (-λ Z * ) Assuming Z = U DV T , with D = diag(d 1 , d 2 , . . . , d r ) this can be further decomposed as p(Z) = p(U )p(V )p(D) where we assume a uniform Haar prior distribution on the unitary matrices U and V , and exponential priors on the singular values d i , hence

p(d 1 , . . . , d r ) = r i=1 Exp (d i ; λ) (5)
where Exp(x; λ) = λ exp(-λx) is the probability density function (pdf) of the exponential distribution of parameter λ evaluated at x. The exponential distribution has a mode at 0, hence favoring sparse solutions.

We propose here alternative penalty/prior distributions, that bridge the gap between the rank and the nuclear norm penalties. Our penalties are based on hierarchical Bayes constructions and the related optimization problems to obtain MAP estimates can be solved by using an EM algorithm.

Hierarchical adaptive spectral penalty

We consider the following hierarchical prior for the low rank matrix Z. We still assume that Z = U DV T , where the unitary matrices U and V are assigned uniform priors and D = diag(d 1 , . . . , d r ).

We now assume that each singular value d i has its own regularization parameter γ i . We further assume that the regularization parameters are themselves iid from a gamma distribution

p(γ 1 , . . . , γ r ) = r i=1 p(γ i ) = r i=1 Gamma(γ i ; a, b)
where Gamma(γ i ; a, b) is the pdf of the gamma distribution of parameters a > 0 and b > 0 evaluated at γ i . The marginal distribution over d i is thus a continuous mixture of exponential distributions

p(d i ) = ∞ 0 Exp(d i ; γ i ) Gamma(γ i ; a, b)dγ i = ab a (d i + b) a+1 (6)
It is a Pareto distribution which has heavier tails than the exponential distribution. Figure 1 shows the marginal distribution p(d i ) for a = b = β. The lower β, the heavier the tails of the distribution. When β → ∞, one recovers the exponential distribution of unit rate parameter. Let

pen(Z) = -log p(Z) = - r i=1 log(p(d i )) = C 1 + r i=1 (a + 1) log(b + d i ) (7)
be the penalty induced by the prior p(Z). We call the penalty [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF] the Hierarchical Adaptive Spectral Penalty (HASP). On Figure 3 (top) are represented the balls of constant penalties for a symmetric 2 × 2 matrix, for the HASP, nuclear norm and rank penalties. When the matrix is assumed to be diagonal, one recovers respectively the lasso, hierarchical adaptive lasso (HAL) [START_REF] Candès | Enhancing sparsity by reweighted l1 minimization[END_REF][START_REF] Lee | A hierarchical Bayesian framework for constructing sparsity-inducing priors[END_REF] and ℓ 0 penalties, as shown on Figure 3 (bottom).

The penalty (7) admits as special cases the nuclear norm penalty λ||Z|| * when a = λb and b → ∞.

Another closely related penalty is the log-det heuristic [START_REF] Fazel | Matrix rank minimization with applications[END_REF][START_REF] Fazel | Log-det heuristic for matrix rank minimization with applications to hankel and euclidean distance matrices[END_REF] penalty, defined for a square matrice Z by log det(Z + δI) where δ is some small regularization constant. Both penalties agree on square matrices when a = b = 0 and δ = 0.

EM algorithm for MAP estimation

Using the exponential mixture representation [START_REF] Candès | Enhancing sparsity by reweighted l1 minimization[END_REF], we now show how to derive an EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] to obtain a MAP estimate

Z = arg max Z [log p(X|Z) + log p(Z)]
i.e. to minimize We use the parameters γ = (γ 1 , . . . , γ r ) as latent variables in the EM algorithm. The E step is obtained by

L(Z) = 1 2σ 2 X -Z 2 F + r i=1 (a + 1) log(b + d i ) (8)
Q(Z, Z * ) = E [log(p(X, Z, γ))|Z * , X] = C 2 - 1 2σ 2 X -Z 2 F - r i=1 E[γ i |d * i ]d i
Hence at each iteration of the EM algorithm, the M step consists in solving the optimization problem minimize

Z 1 2σ 2 X -Z 2 F + r i=1 ω i d i (9)
where

ω i = E[γ i |d * i ] = ∂ ∂d * i [-log p(d * i )] = a+1 b+d * i . ( 9
) is an adaptive nuclear norm regularized optimization problem, with weights ω i . Without loss of generality, assume that

d * 1 ≥ d * 2 ≥ . . . ≥ d * r . It implies that 0 ≤ ω 1 ≤ ω 2 ≤ . . . ≤ ω r (10)
The above weights will therefore penalize less heavily higher singular values, hence reducing bias. As shown by [START_REF] Gaïffas | Weighted algorithms for compressed sensing and matrix completion[END_REF][START_REF] Chen | Reduced rank regression via adaptive nuclear norm penalization[END_REF], a global optimal solution to Eq. ( 9) under the order constraint [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is given by a weighted soft-thresholded SVD

Z = S σ 2 ω (X) (11)
where

S ω (X) = U D ω V T with D ω = diag(( d 1 -ω 1 ) + , . . . , ( d r -ω r ) + ). X = U D V T is the SVD of X with D = diag( d 1 , . . . , d r ) and d 1 ≥ d 2 . . . ≥ d r .
Algorithm 1 summarizes the Hierarchical Adaptive Soft Thresholded (HAST) procedure to converge to a local minimum of the objective [START_REF] Lee | A hierarchical Bayesian framework for constructing sparsity-inducing priors[END_REF]. This algorithm admits the soft-thresholded SVD operator as a special case when a = bλ and b = β → ∞. Figure 2 shows the thresholding rule applied to the singular values of X for the HAST algorithm (a = b = β, with β = 2 and β = 0.1) and the soft-thresholded SVD for λ = 1. The bias term, which is equal to λ for the nuclear norm, goes to zero as d i goes to infinity.

Setting of the hyperparameters and initialization of the EM algorithm

In the experiments, we have set b = β and a = λβ where λ and β are tuning parameters that can be chosen by crossvalidation. As λ is the mean value of the regularization parameter γ i , we initialize the algorithm with the soft thresholded SVD with parameter σ 2 λ. It is possible to estimate the hyperparameter σ within the EM algorithm as described in the supplementary material. In our experiments, we have found the results not very sensitive to the setting of σ, and set it to 1.

Algorithm 1 Hierarchical Adaptive Soft Thresholded (HAST) algorithm for low rank estimation of complete matrices

Initialize Z (0) . At iteration t ≥ 1

• For i = 1, . . . , r, compute the weights ω

(t) i = a+1 b+d (t-1) i • Set Z (t) = S σ 2 ω (t) (X) • If L(Z (t-1) )-L(Z (t) )
L(Z (t-1) )

< ε then return Z = Z (t)

Matrix completion

We now show how the EM algorithm derived in the previous section can be adapted to the case where only a subset of the entries is observed. It relies on imputing missing values, similarly to the EM algorithm for SVD with missing data, see e.g. [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Srebro | Weighted low-rank approximations[END_REF].

Consider that only a subset Ω ⊂ {1, . . . , m} × {1, . . . , n} of the entries of the matrix X is observed.

Similarly to [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF], we introduce the operator P Ω (X) and its complementary P ⊥ Ω (X)

P Ω (X)(i, j) = X ij if (i, j) ∈ Ω 0 otherwise and P ⊥ Ω (X)(i, j) = 0 if (i, j) ∈ Ω X ij otherwise
Assuming the same prior (6), the MAP estimate is obtained by minimizing

L(Z) = 1 2σ 2 P Ω (X) -P Ω (Z) 2 F + (a + 1) r i=1 log(b + d i ) (12)
We will now derive the EM algorithm, by using latent variables γ and P ⊥ Ω (X). The E step is given by (details in supplementary material)

Q(Z, Z * ) = E log(p(P Ω (X), P ⊥ Ω (X), Z, γ))|Z * , P Ω (X) = C 4 - 1 2σ 2 P Ω (X) + P ⊥ Ω (Z * ) -Z 2 F - r i=1 E[γ i |d * i ]d i
Hence at each iteration of the algorithm, one needs to minimize

1 2σ 2 X * -Z 2 F + r i=1 ω i d i (13)
where ω i = E[γ i |d * i ] and X * = P Ω (X) + P ⊥ Ω (Z * ) is the observed matrix, completed with entries in Z * . We now have a complete matrix problem. As mentioned in the previous section, the minimum of ( 13) is obtained with a weighted soft-thresholded SVD. Algorithm 2 provides the resulting iterative procedure for matrix completion with the hierarchical adaptive spectral penalty.

Algorithm 2 Hierarchical Adaptive Soft Impute (HASI) algorithm for matrix completion

Initialize Z (0) . At iteration t ≥ 1

• For i = 1, . . . , r, compute the weights ω

(t) i = a+1 b+d (t-1) i • Set Z (t) = S σ 2 ω (t) P Ω (X) + P ⊥ Ω (Z (t-1) ) • If L(Z (t-1) )-L(Z (t) ) L(Z (t-1) ) < ε then return Z = Z (t)
Related algorithms Algorithm 2 admits the Soft-Impute algorithm of [START_REF] Mazumder | Spectral regularization algorithms for learning large incomplete matrices[END_REF] as a special case when a = λb and b = β → ∞. In this case, one obtains at each iteration ω (t) i = λ for all i. On the contrary, when β < ∞, our algorithm adaptively updates the weights so that to penalize less heavily higher singular values. Some authors have proposed related one-step adaptive spectral penalty algorithms [START_REF] Bach | Consistency of trace norm minimization[END_REF][START_REF] Gaïffas | Weighted algorithms for compressed sensing and matrix completion[END_REF][START_REF] Chen | Reduced rank regression via adaptive nuclear norm penalization[END_REF]. However, in these procedures, the weights have to be chosen by some procedure whereas in our case they are iteratively adapted.

Initialization The objective function [START_REF] Chen | Reduced rank regression via adaptive nuclear norm penalization[END_REF] is in general not convex and different initializations may lead to different modes. As in the complete case, we suggest to set a = λb and b = β and to initialize the algorithm with the Soft-Impute algorithm with regularization parameter σ 2 λ.

Scaling Similarly to the Soft-Impute algorithm, the computationally demanding part of Algorithm 2 is S σ 2 ω (t) P Ω (X) + P ⊥ Ω (Z (t-1) ) which requires calculating a low rank truncated SVD. For large matrices, one can resort to the PROPACK algorithm [START_REF] Larsen | Lanczos bidiagonalization with partial reorthogonalization[END_REF][START_REF] Larsen | Propack-software for large and sparse svd calculations[END_REF] as described in [START_REF] Mazumder | Spectral regularization algorithms for learning large incomplete matrices[END_REF]. This sophisticated linear algebra algorithm can efficiently compute the truncated SVD of the "sparse + low rank" matrix P Ω (X) + P ⊥ Ω (Z (t-1) ) = P Ω (X) -P Ω (Z (t-1) )

sparse

+ Z (t-1)
low rank and can thus handle large matrices, as shown in [START_REF] Mazumder | Spectral regularization algorithms for learning large incomplete matrices[END_REF].

Experiments

Simulated data

We first evaluate the performance of the proposed approach on simulated data. Our simulation setting is similar to that of [START_REF] Mazumder | Spectral regularization algorithms for learning large incomplete matrices[END_REF]. We generate Gaussian matrices A and B respectively of size m × q and n × q, q ≤ r so that the matrix Z = AB T is of low rank q. A Gaussian noise of variance σ 2 is then added to the entries of Z to obtain the matrix X. The signal to noise ratio is defined as

SNR = var(Z) σ 2 .
We set m = n = 100 and σ = 1. We run all the algorithms with a precision ǫ = 10 -9 and a maximum number of t max = 200 iterations (initialization included for HASI). We compute err, the relative error between the estimated matrix Z and the true matrix Z in the complete case, and err Ω ⊥ in the incomplete case, where

err = || Z -Z|| 2 F ||Z|| 2 F and err Ω ⊥ = || P ⊥ Ω ( Z) -P ⊥ Ω (Z)|| 2 F ||P ⊥ Ω (Z)|| 2 F
For the HASP penalty, we set a = λβ and b = β. We compute the solutions over a grid of 50 values of the regularization parameter λ linearly spaced from λ 0 to 0, where λ 0 = ||P Ω (X)|| 2 is the largest singular value of the input matrix X, padded with zeros. This is done for three different values β = 1, 10, 100. We use the same grid to obtain the regularization path for the other algorithms.

Complete case

We first consider that the observed matrix is complete, with SNR = 1 and q = 10.

The HAST algorithm 1 is compared to the soft thresholded (ST) and hard thresholded (HT) SVD. Results are reported in Figure 4(a). The HASP penalty provides a bridge/tradeoff between the nuclear norm and the rank penalty. For example, value of β = 10 show a minimum at the true rank q = 10 as HT, but with a lower error when the rank is overestimated. Incomplete case Then we consider the matrix completion problem, and remove uniformly a given percentage of the entries in X. We compare the HASI algorithm to the Soft-Impute, Soft-Impute+ and Hard-Impute algorithms of [START_REF] Mazumder | Spectral regularization algorithms for learning large incomplete matrices[END_REF] and to the MMMF algorithm of [START_REF] Rennie | Fast maximum margin matrix factorization for collaborative prediction[END_REF]. Results, averaged over 50 replications, are reported in Figures 4(b-c) for a true rank q = 5, (b) 50% of missing data and SNR = 1 and (c) 80% of missing data and SNR = 10. Similar behavior is observed, with the HASI algorithm attaining a minimum at the true rank q = 5. We then conduct the same experiments, but remove 20% of the observed entries as a validation set to estimate the regularization parameters (λ, β) for HASI, and λ for the other methods. We estimate Z on the whole observed matrix, and use the unobserved entries as a test set. Results on the test error and estimated ranks over 50 replications are reported in Figure 5. For 50% missing data, HASI is shown to outperform the other methods. For 80% missing data, HASI and Hard Impute provide the best performances. In both cases, it is able to recover very accurately the true rank of the matrix.

Collaborative filtering examples

We now compare the different methods on several benchmark datasets. We first consider the Jester datasets [START_REF] Goldberg | Eigentaste: A constant time collaborative filtering algorithm[END_REF]. The three datasets1 contain one hundred jokes, with user ratings between -10 and +10. We randomly select two ratings per user as a test set, and two other ratings per user as a validation set to select the parameters λ and β. The results are computed over four values β = 1000, 100, 10, 1.

We compare the results of the different methods with the Normalized Mean Absolute Error (NMAE)

NMAE = 1 card(Ωtest) (i,j)∈Ωtest |X ij -Z ij | max(X) -min(X)
where Ω test is the test set. The mean number of iterations for Soft-Impute, Hard-Impute and HASI (initialization included) algorithms are respectively 9, 76 and 76. Computations for the HASI algorithm take approximately 5 hours on a standard computer. The results, averaged over 10 replications (with almost no variability observed), are presented in Table 1. The HASI algorithm provides very good performance on the different Jester datasets, with lower NMAE than the other methods. Second, we conducted the same comparison on two MovieLens datasets2 , which contain ratings of movies by users. We randomly select 20% of the entries as a test set, and the remaining entries are split between a training set (80%) and a validation set (20%). For all the methods, we stop the regularization path as soon as the estimated rank exceeds r max = 100. This is a practical consideration: given that the computations for high ranks demand more time and memory, we are interested in restricting ourselves to low rank solutions. Table 1 presents the results, averaged over 5 replications. For the MovieLens 100k dataset, HASI provides better NMAE than the other methods with a low rank solution. For the larger MovieLens 1M dataset, the precision, maximum number of iterations and maximum rank are decreased to ǫ = 10 -6 , t max = 100 and r max = 30. On this dataset, MMMF provides the best NMAE at maximum rank. HASI provides the second best performances with a slightly lower rank.

Conclusion

The proposed class of methods has shown to provide good results compared to several alternative low rank matrix completion methods. It provides a bridge between nuclear norm and rank regularization algorithms. Although the related optimization problem is not convex, experiments show that initializing the algorithm with the Soft-Impute algorithm of [START_REF] Mazumder | Spectral regularization algorithms for learning large incomplete matrices[END_REF] provides very satisfactory results.

In this paper, we have focused on a gamma mixture of exponentials, as it leads to a simple and interpretable expression for the weights. It is however possible to generalize the results presented here by using a three parameter generalized inverse Gaussian prior distribution (see e.g. [START_REF] Zhang | EP-GIG priors and applications in Bayesian sparse learning[END_REF]) for the regularization parameters γ i , thus offering an additional degree of freedom. Derivations of the weights are provided in the supplementary material. Additionally, it is possible to derive an EM algorithm for low rank matrix completion for binary matrices. Details are also provided in supplementary material.

While we focus on point estimation in this paper, it would be of interest to investigate a fully Bayesian approach and derive a Gibbs sampler or variational algorithm to approximate the posterior distribution, and compare to other full Bayesian approaches to matrix completion [START_REF] Seeger | Fast variational Bayesian inference for non-conjugate matrix factorization models[END_REF][START_REF] Nakajima | Global analytic solution of fullyobserved variational Bayesian matrix factorization[END_REF].

Figure 1 :

 1 Figure 1: Marginal distribution p(d i ) with a = b = β for different values of the parameter β. The distribution becomes more concentrated around zero with heavier tails as β decreases. The case β → ∞ corresponds to an exponential distribution with unit rate.

Figure 2 :

 2 Figure 2: Thresholding rules on the singular values d i of X for the soft thresholding rule (λ = 1), and hierarchical adaptive soft thresholding algorithm with a = b = β.

  Figure 3: Top: Manifold of constant penalty, for a symmetric 2 × 2 matrix Z = [x, y; y, z] for (a) the nuclear norm, (b-c) hierarchical adaptive spectral penalty with a = b = β (b) β = 1 and (c) β = 0.1, and (d) the rank penalty. Bottom: contour of constant penalty for a diagonal matrix [x, 0; 0, z], where one recovers the classical (e) lasso, (f-g) hierarchical lasso and (h) ℓ 0 penalties.

Figure 4 :

 4 Figure 4: Test error w.r.t. the rank obtained by varying the value of the regularization parameter λ. Results on simulated data are given for (a) complete matrix with SNR=1 (b) 50% missing and SNR=1 and (c) 80% missing and SNR=10.

Figure 5 :

 5 Figure 5: Boxplots of the test error and ranks obtained over 50 replications on simulated data.

Figure 6

 6 Figure6shows the NMAE in function of the rank. Low values of β exhibit a bimodal behavior with two modes at low rank and full rank. High value β = 1000 is unimodal and outperforms Soft-Impute at any particular rank.

3 Figure 6 :

 36 Figure 6: NMAE on the test set of the (a) Jester 1 and (b) Jester 3 datasets w.r.t. the rank obtained by varying the value of the regularization parameter λ. The curves obtained on the Jester 2 dataset are hardly distinguishable from (a) and hence are not displayed due to space limitation.

Table 1 :

 1 Results on the Jester and MovieLens datasets

		Jester 1	Jester 2	Jester 3	MovieLens 100k MovieLens 1M
		24983 × 100	23500 × 100	24938 × 100	943 × 1682	6040 × 3952
		27.5% miss.	27.3% miss.	75.3% miss.	93.7% miss.	95.8% miss.
	Method	NMAE Rank NMAE Rank NMAE Rank	NMAE Rank	NMAE Rank
	MMMF	0.161	95	0.162	96	0.183	58	0.195	50	0.169	30
	Soft Imp	0.161	100	0.162	100	0.184	78	0.197	156	0.176	30
	Soft Imp+ 0.169	14	0.171	11	0.184	33	0.197	108	0.189	30
	Hard Imp	0.158	7	0.159	6	0.181	4	0.190	7	0.175	8
	HASI	0.153	100	0.153	100	0.174	30	0.187	35	0.172	27

Jester datasets can be downloaded from the url http://goldberg.berkeley.edu/jester-data/

MovieLens datasets can be downloaded from the url http://www.grouplens.org/node/73.

Acknowledgments

F.C. acknowledges the support of the European Commission under the Marie Curie Intra-European Fellowship Programme. The contents reflect only the authors views and not the views of the European Commission.