Robin Chatelin

Dominique Anne-Archard

Marlène Murris-Espin

Marc Thiriet

Philippe Poncet

Numerical and experimental investigation of mucociliary clearance breakdown in cystic fibrosis

Keywords: Computational biology, Cystic fibrosis, Mucociliary clearance, Modeling, Rheology, Generalized Stokes equations, Penalization, Grid-particle methods

The human tracheobronchial tree surface is covered with mucus. A healthy mucus is an heterogeneous material flowing toward the esophagus and a major defense actor against local pathogen proliferation and pollutant deposition. An alteration of mucus or its environment such as in cystic fibrosis dramatically impacts the mucociliary clearance. In the present study, we investigate the mechanical organization and the physics of such mucus in human lungs by means of a joint experimental and numerical work. In particular, we focus on the influence of the shear-thinning mucus mobilized by a ciliated epithelium for mucociliary clearance. The proposed robust numerical method is able to manage variations of more than 5 orders of magnitude in the shear rate and viscosity. It leads to a cartography that allows to discuss major issues on defective mucociliary clearance in cystic fibrosis. Furthermore, the computational rheological analysis based on measurements shows that cystic fibrosis shear-thinning mucus tends to aggregate in regions of lower clearance. Yet, a rarefaction of periciliary fluid has a greater impact than the mucus shear-thinning effects.

Introduction

The human respiratory tract (nose, pharynx, larynx, trachea, bronchi, bronchioles, and alveoli) and its related cavities and conduits (sinuses and Eustachian tubes) has a wetted surface covered with mucus flowing toward the esophagus. Human conducting airways are lined with a pseudostratified, secretory, and ciliated epithelium (with ciliated, secretory, and basal cells) coating submucosal glands as well as cartilaginous structures in some segments of the respiratory tract such as partial or complete rings in the proximal part of the tracheo-bronchial tree.

An efficient mucociliary transport requires an appropriate mucus composition for optimal clearance, adequate production of periciliary fluid and number of functioning ciliated cells, and coordinated motion of cilia for mucus propulsion, as well as respiratory epithelium integrity over long distances. The mucociliary clearance protects the airway epithelium against inhaled pathogens and noxious entrapped micro-and nanoparticles. The mucociliary clearance also influences drug delivery to airways, drug absorption at a given station of the respiratory tract depending on the local rate of mucus transport. 1 Corresponding author: philippe.poncet@univ-pau.fr

The present study focuses on impaired mucociliary transport in cystic fibrosis. Cystic fibrosis (CF) is one of the most common life-shortening genetic diseases. It is observed on 1 in every 2500 infants in Europe and North America. It is most common among Caucasian population (1 in 22 individuals is a heterozygous carrier). The heterozygote frequency reaches 4% of the population of Western world. The prevalence of cystic fibrosis per 10000 inhabitants in the United States of America is 0.797 and is on average 0.737 in the 27 European Union countries (with a remarkable 2.98 in the Republic of Ireland) [START_REF] Farrell | The prevalence of cystic fibrosis in the European Union[END_REF].

Defective mucus composition alters mucociliary transport. Impaired chloride transport in cystic fibrosis prevents normal hydration of mucus, leading to highly viscous accumulations. Cystic fibrosis airways are characterized by airway surface liquid volume depletion and mucus accumulation, thereby hampering airway clearance of inhaled pollutants and pathogens [START_REF] Derichs | Hyperviscous airway periciliary and mucous liquid layers in cystic fibrosis measured by confocal fluorescence photobleaching[END_REF].

Our numerical and rheological investigations aim at providing an independent and self-consistent analysis of mucociliary clearance altered by cystic fibrosis allowing full reproducibility. Questions needed to be answered encompass the role on mucus viscosity of the disease itself, associated infection, and influence of the periciliary fluid (PCL). Current models assume two Newtonian layers of different viscosities or a single complex fluid [START_REF] Li | An immersed boundary-lattice Boltzmann method for single-and multicomponent fluid flows[END_REF], or two layers with a constant visco-elastic parameter superposed to a spatially varying diffusion [START_REF] Guo | Mucociliary Transport in Healthy and Diseased Environments[END_REF]. In the present article, a rheological hydration/maturation law incorporates the mucin transfer through a Newtonian periciliary fluid and mucin meshwork formation in a non-Newtonian mucus layer (ML) with smooth and sharp transitions.

Many rheological laws have been identified or used in the literature, from shear-thinning to viscoelastic models but nonshear-thinning (Maxwell [START_REF] Mitran | Metachronal wave formation in a model of pulmonary cilia[END_REF][START_REF] Smith | Modelling mucociliary clearance[END_REF], Oldroyd B model [START_REF] Norton | Model of ciliary clearance and the role of mucus rheology[END_REF][START_REF] Nawroth | Beyond the mucus escalator: Complex ciliary hydrodynamics in disease and function[END_REF][START_REF] Guo | Mucociliary Transport in Healthy and Diseased Environments[END_REF]), viscoplastic (Bingham [START_REF] Zamankhan | Steady motion of Bingham liquid plugs in two-dimensional channels[END_REF], investigation of 2D mucus without cilia, Carreau [START_REF] Montenegro-Johnson | Physics of rheologically enhanced propulsion: Different strokes in generalized Stokes[END_REF] with spatially constant parameters, and Herschel-Bulkley [START_REF] Kocevar-Nared | Comparative rheological investigation of crude gastric mucin and natural gastric mucus[END_REF] constitutive laws). Mucus rheology depends on the context (mucus conservation, illness severity, treatment effect, mucus source, sampling mode, etc.). As mentioned by Besseris and Yeates [START_REF] Besseris | Rotating magnetic particle microrheometry in biopolymer fluid dynamics: Mucus microrheology[END_REF]: "no constitutive law ever was proposed for this material because it was never clearly classified either as a chemically or physically linked gel".

The major issue of this study is to quantify the mucociliary clearance and to investigate the influence of the highly shearthinning nature of mucus. To do that, a state-of-the-art solver is used to compute the variable-viscosity Stokes problem in a mobile domain due to beating cilia [START_REF] Chatelin | A hybrid grid-particle method for moving bodies in 3D Stokes flow with variable viscosity[END_REF].

This model deals with noncartilaginous ciliated walls of the tracheobronchial tree, where the ciliated surface is assumed to be periodically flat.

This article is structured as follows. Firstly, the experimental method is described from samples collection to measurement validation. Secondly, we introduce the mathematical model for flow and maturation which enables the transition from the Newtonian PCL to non-Newtonian ML. The numerical method solving efficiently this heterogeneous fluid dynamics around beating cilia is then described. It involves only standard elliptic problems and avoid matrix assembly, thereby leading to a low storage method [START_REF] Chatelin | Hybrid grid-particle methods and penalization: A Sherman-Morrison-Woodbury approach to compute 3D viscous flows using FFT[END_REF]. This is then applied to a wide range of rheological parameters, which gives a cartography of mucociliary efficiency. This allows comparison of our data with literature data.

In addition to the agreement between experimental observations and numerical parametric studies, the situation of healthy and pathological mucus is achieved on a map. A collection of these maps is finally provided in order to quantify the impact of PCL reduction on mucociliary clearance.

Rheology measurements and experimental setup

The main and mandatory property for numerical simulations of shear-thinning fluids is the viscosity function η(γ) with respect to the shear-rate γ, the full contraction of the symmetric part of the gradient of the fluid velocity u:

γ = (2D(u):D(u)) 1/2
where D(u) = (∇u + ∇u T )/2. Two types of equipments are used to assess this parameter: classical rheometers with a cone and plate or parallel plates geometry [START_REF] Moores | The effect of highfrequency ventilation on non-Newtonian properties of bronchial mucus[END_REF][START_REF] Serisier | Macrorheology of cystic fibrosis, chronic obstructive pulmonary disease and normal sputum[END_REF][START_REF] Suk | The penetration of fresh undiluted sputum expectorated by cystic fibrosis patients by non-adhesive polymer nanoparticles[END_REF][START_REF] Tomaiuolo | A new method to improve the clinical evaluation of cystic fibrosis patients by mucus viscoelastic properties[END_REF], and magnetic microrheometers which only allow small amplitude oscillatory motions [START_REF] Jeanneret-Grosjean | Sampling technique and rheology of human tracheobronchial mucus[END_REF][START_REF] Rubin | Collection and analysis of respiratory mucus from subjects without lung disease[END_REF][START_REF] Zayas | Tracheal mucus rheology in patients undergoing diagnostic bronchoscopy. Interrelations with smoking and cancer[END_REF].

Experiments were carried out on a Mars III rheometer (Thermo Sc.) equipped with a cone and plate geometry (angle 1 o , diameter 35 mm) which requires a small volume of fluid (0.2 ml). Measurements were conducted at the body temperature (37 o C). The lower plate Peltier system for temperature regulation was completed by a home-made sample hood with a specific temperature regulation and wet sponges to saturate the measuring chamber and to prevent samples evaporation.

The rheological characterization consists of dynamic measurements (small amplitude oscillations) with determination of the linear domain followed by measurement of the storage G and loss modulus G as functions of the oscillation frequency ω. G (ω) and G (ω) are respectively the elastic (in-phase) and viscous (out-of-phase) components of the complex modulus G * = G + iG which is the ratio between shear stress and strain during oscillation measurements. The viscosity η(γ) was then determined by a steady shear flow test for shear rates ranging from 10 -3 s -1 to 10 3 s -1 . A special attention was paid to reach a steady state regime for low shear rate values.

Methodology validation

Mucus was sampled during consultation (in stabilized CF patients as well as in a patient with bronchiectasis without CF) or hospitalization (both decompensated or stabilized posttherapeutic states). Samples were frozen before transport and processing. Comparison between fresh and frozen mucus (freezing duration 21 and 30 days) were carried out on three samples and confirmed that the freezing does not alter the rheological propoerties as it was already observed by Sanders et al. [START_REF] Sanders | Cystic fibrosis sputum: a barrier to the transport of nanospheres[END_REF] and Rubin et al. [START_REF] Rubin | Collection and analysis of respiratory mucus from subjects without lung disease[END_REF].

Thixotropy is sometimes mentioned in mucus studies and related measurements [START_REF] Puchelle | Elasto-thixotropic properties of bronchial mucus and polymer analogs. I: Experimental results[END_REF]. For a given patient, when the mucus sample volume was sufficient, shear viscosity was successively measured on the same and on different samples. The dispersion observed using these two series of measures was similar, that is, thixotropy was comparable to the natural dispersion arising from heterogeneity of the sample. The weak thixotropy was already observed in viscosity measurements for increasing and decreasing shear rate [START_REF] Yeates | Physicochemical properties of mucus and its propulsion[END_REF]. In addition, no phase separation was observed during defrosting and during measurement.

Fifteen mucus samples were analyzed and among them two were rejected. They were visually too heterogeneous and gave rise to a large dispersion in viscosity and modulus values.

Rheological results

A shear-thinning behavior is observed in the entire shear rate range. The storage modulus G (ω) and the viscous modulus G (ω) clearly display the viscoelastic character of mucus. These curves cross each other at a frequency ω c . The characteristic time defined from ω c is of the order 10 2 s.

A power law model η = K γN-1 was fitted on the shear rate range [0.01 -100 s -1 ], which is representative of shear rate involved in the studied flow. Values for the power-law index N and the consistency K are obtained by mean-square log-fitting of measures (see Fig. 1 1: Rheological parameters (K, N) from collected samples (CFM1-CFM11 standing for Cystic Fibrosis, and BM12 standing for Bronchiectasis) and the literature (CFL1-CFL5 and HL6-HL9 from the complex modulus G*). Aliases CFL and HL stand for "Cystic Fibrosis from Literature" and "Healthy cases from Literature". three measurements, each on different mucus samples are displayed in Fig. 1. These curves illustrate both the usual intraand inter-subject variability of any biological material.

Few literature data on human airway mucus are available for steady state viscosity, whereas numerous data exist for complex viscosity η * = G * /ω, giving

| η * | = | G * |/ω = G (ω) 2 + G (ω) 2 / ω (1)
The | η * | modulus is considered as an admissible estimation of viscosity when shear viscosity is not obtainable (Cox-Merz rule [START_REF] Bird | Dynamics of Polymeric Liquids[END_REF]), and this was largely used. Nevertheless, all authors recognize a strong shear-thinning behavior. This is illustrated by the power law parameters extracted from literature data on G and G and displayed on Table 1, which constitutes a comparison basis for the shear viscosity measured in this study. Whereas some authors [START_REF] Moores | The effect of highfrequency ventilation on non-Newtonian properties of bronchial mucus[END_REF][START_REF] Dawson | Enhanced viscoelasticity of human cystic fibrotic sputum correlates with increasing microheterogeneity in particle transport[END_REF][START_REF] Suk | The penetration of fresh undiluted sputum expectorated by cystic fibrosis patients by non-adhesive polymer nanoparticles[END_REF] observed a small Newtonian plateau at low shear rates (0.1 -1.0 s -1 ), others [START_REF] Yeates | Physicochemical properties of mucus and its propulsion[END_REF][START_REF] Puchelle | Elasto-thixotropic properties of bronchial mucus and polymer analogs. I: Experimental results[END_REF][START_REF] Puchelle | Rheological properties controlling mucociliary frequency and respiratory mucus transport[END_REF] did not monitor any trend to a plateau although exploring shear rate from 0.01 s -1 . Our own experience does not reveal any plateau even for shear rate values of 10 -3 s -1 . This behavior is generally linked to viscoplasticity, but the τ(γ) curves obtained on the tested samples do not display any tendency to a yield stress. So, from our measurements, it appears that, if a yield stress exists, it is likely very low (typically < 0.01 Pa).

Figure 1: Curves from CFM2 and CFM7 measurements from Table 1 with their best-fitting power-law on γ ∈ [10 -2 , 10 2 ].

The power-law index N values reported in Table 1, which derive from rheograms η(γ) ranges from 0.2 to 0.4; they are similar to literature data [START_REF] Puchelle | Rheological properties controlling mucociliary frequency and respiratory mucus transport[END_REF][START_REF] Suk | The penetration of fresh undiluted sputum expectorated by cystic fibrosis patients by non-adhesive polymer nanoparticles[END_REF]. Values of consistency K and viscosity are also similar. All viscosity curves are limited by two envelope curves defined by (K min = 0.45, N min = 0.39) and (K max = 5.8, N max = 0.20). The large shear rate range explored in this study gives an appropriate description of viscosity for the mucus film displacement.

Modeling of mucus maturation

This section first presents the airway surface liquid structure, composed of mucus and periciliary fluid surrounding the lung epithelial cells. The two next subsections present how the fluid velocity u and the fluid viscosity η G are linked, and finally how the viscosity η G is modeled for such a mixture, by a Carreau law with spatially varying coefficients.

Airway surface liquid

The airway surface liquid (ASL) covering bronchial walls is a stratified medium composed of a mucus layer (ML) overlying a periciliary fluid layer (PCL), which behaves as a Newtonian fluid with a viscosity close to that of water. Both PCL and ML are chiefly composed of water (>98%). Several works in the literature [START_REF] Fulford | Muco-ciliary transport in the lung[END_REF][START_REF] Craster | Surfactant transport on mucus films[END_REF][START_REF] Smith | A viscoelastic traction layer model of muco-ciliary transport[END_REF][START_REF] Jayathilake | Threedimensional numerical simulations of human pulmonary cilia in the periciliary liquid layer by the immersed boundary method[END_REF][START_REF] Li | An immersed boundary-lattice Boltzmann method for single-and multicomponent fluid flows[END_REF] (an exhaustive reference list is given in [START_REF] Smith | Modelling mucociliary clearance[END_REF]) have considered two layers of fluid with constant fluid properties.

The periciliary fluid wetting the bronchial wall has a height generally estimated to be at least 2/3 of the total ASL thickness. Once they are released by goblet cells, polymerized mucins are further hydrated, expand, and travel freely in the PCL away from the wetted surface [START_REF] Davis | Regulated Airway Goblet Cell Mucin Secretion[END_REF][START_REF] Davis | Coupling of airway ciliary activity and mucin secretion to mechanical stresses by purinergic signaling[END_REF][START_REF] Ehre | Cystic fibrosis, An inherited disease affecting mucin-producing organs[END_REF][START_REF] Bansil | Mucin Biophysics[END_REF].

The innermost mucus layer results from progressive maturation and structuration of mucins in the core of the fluid. Yet, the transition between the PCL and ML layers is supposed to be continuous, assuming a progressive mucin maturation between the epithelial surface and mucus, characterized by a progressive meshwork loop size [START_REF] Button | Periciliary brush promotes the lung health by separating the mucus layer from airway epithelia[END_REF]. 13) coupled to the 3D Stokes equation [START_REF] Button | Periciliary brush promotes the lung health by separating the mucus layer from airway epithelia[END_REF]. Bottom picture: static mixing given by Eq. ( 12).

The present model considers this continuous heterogeneity. The ASL rheological parameters depends on the mucin maturation parameter α, defined as the volumic fraction of matured mucins at a given ASL height, that is, on the cross-linking magnitude of mucin polymers [START_REF] Georgiades | Particle tracking microrheology of purified gastrointestinal mucins[END_REF][START_REF] Raynal | Calcium-dependent protein interactions in MUC5B provide reversible cross-links in salivary mucus[END_REF][START_REF] Verdugo | Molecular mechanism of product storage and release in mucin secretion. II. The role of extracellular Ca++[END_REF]. The maturation α varies from 0 at the PCL-epithelium interface to 1 in the ML with a sigmoid or tan -1 profile.

Flow modeling and domain geometry

The computational domain Ω contains a full ciliated cell (a segment of the respiratory epithelium, cf. Fig. 2). The ciliated cell is the time-dependent domain B(t) and the ASL domain

F (t) = Ω \ B(t).
The rectangular domain is bounded by the mucus-air interface at its top and the wetted epithelial surface at its bottom, with periodical conditions in the streamwise direction.

The Reynolds, Womersley, and Froude numbers are sufficiently small to suppose a mucus regime following the stationary 3D Stokes equations in a time-dependent fluid domain F (t), that is to say the flow is always at its quasi-equilibrium [START_REF] Enault | Mucus dynamics subject to air and wall motion[END_REF]:

-div(2η G D(u)) + ∇p = f in F (t), ( 2
)
where u is the velocity field satisfying divu = 0 in F (t), p the pressure, D(u) = (∇u + ∇u T )/2 the strain rate tensor, η G the dynamic viscosity, and u = ū on boundary ∂F (t method [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF] is then used to take into account the influence of the beating cilia on the fluid. It is a robust fictitious domain method compatible with any spatial discretization, in particular regular grids which permits to use fast solvers, as detailed below.

The cilia velocity ū is computed from the beating model presented in [START_REF] Chatelin | A hybrid grid-particle method for moving bodies in 3D Stokes flow with variable viscosity[END_REF], based on damped Young modulus propagation, and compatible with the cilia beating cycle [START_REF] Fauci | Biofluidmechanics of Reproduction[END_REF][START_REF] Gheber | Extraction of cilium beat parameters by the combined application of photoelectric measurements and computer simulation[END_REF][START_REF] Sanderson | Ciliary activity of cultured rabbit tracheal epithelium, Beat pattern and metachrony[END_REF][START_REF] Mitran | Metachronal wave formation in a model of pulmonary cilia[END_REF]. Since the cilia velocities are prescribed, this is a one-way fluidstructure interaction. This assumption is valid only for early stage of the sickness or for stabilized patients, that is to say for moderate deviation from the healthy case, and does not concern critical care. In particular, the modification and adaptation of the beating patterns and frequency with respect to the modification of the mucus behavior (viscosity and non-Newtonian effect) is not taken into account in this numerical model. This had been investigated experimentally [START_REF] Hill | Force Generation and Dynamics of Individual Cilia under External Loading[END_REF], and the development of a coherent numerical model remains, to our knowledge, an open problem.

The set of cilia is beating, involving a sufficiently large number of cilia, conforming in-vivo [START_REF] Sanderson | Ciliary activity of cultured rabbit tracheal epithelium, Beat pattern and metachrony[END_REF] and in-vitro [START_REF] Hussong | Experimental investigation of the flow induced by artificial cilia[END_REF] cilia behavior, leading to the flow displayed on Fig. 2. The cell we consider involves 25 cilia, yielding to a sufficiently dense cilia array to get a velocity no longer depending on the number of cilia. Their configuration and simulation parameters are displayed on Table 2. Cilia beat in a coordinated manner and the resulting beating exhibits an antipleptic wave. This parameter remains constant in all the computations, its influence on the mucuciliary clearance was already investigated in the literature [START_REF] Hussong | Experimental investigation of the flow induced by artificial cilia[END_REF][START_REF] Mitran | Metachronal wave formation in a model of pulmonary cilia[END_REF][START_REF] Smith | Modelling mucociliary clearance[END_REF][START_REF] Guo | Cilia beating patterns are not hydrodynamically optimal[END_REF][START_REF] Ding | Mixing and transport by ciliary carpets: a numerical study[END_REF].

Rheology and maturation modeling

As mentioned in introduction and in section 2.2, the mucus layer exhibits both shear-thinning, visco-elastic and/or viscoplastic features. While several studies have focused on the visco-elastic aspects [START_REF] Guo | Mucociliary Transport in Healthy and Diseased Environments[END_REF][START_REF] Guo | Cilia beating patterns are not hydrodynamically optimal[END_REF][START_REF] Ding | Mixing and transport by ciliary carpets: a numerical study[END_REF][START_REF] Smith | A viscoelastic traction layer model of muco-ciliary transport[END_REF][START_REF] Smith | Modelling mucociliary clearance[END_REF][START_REF] Sedaghat | On the effect of mucus rheology on the mucociliary transport[END_REF], visco-plastic effects are less investigated (see [START_REF] Montenegro-Johnson | Physics of rheologically enhanced propulsion: Different strokes in generalized Stokes[END_REF] for spatially constant Carreau law in 2D films and [START_REF] Craster | Surfactant transport on mucus films[END_REF] for analytic computations). The present study focus on the shear-thinning behavior, and bring a new light on this kind of non-Newtonian feature. In order to take in account the transition between the Newtonian PCL and the shear-thinning mucus, an exponential growth is considered between the two fluids of viscosities η PCL and η ML (γ) = K ML γN ML -1 . The parameters of the power-law model depend on α:

K(α) = η PCL K ML η PCL α and N(α) = α • N ML + (1 -α) (3)
where η PCL = 10 -3 Pa • s is the Newtonian solvent viscosity assumed to the water viscosity.

The power-law model needs to be regularized in the vicinity of zero as viscosity tends to infinity. Otherwise, the effective viscosity has a lower bound at high shear rates, the solvent viscosity η PCL . The Carreau model matches these two features:

η G (α, u) = η PCL + (η 0 (α) -η PCL ) 1 + (βγ) 2 N(α)-1 2 (4)
where the regularization parameter β is the inverse of the cutshear rate (low shear rates; β = 4 × 10 3 s). The parameter η 0 (α) is the cut-viscosity and is determined by the choice of the power-law coefficients and cut-shear rate. This leads to:

η 0 (α) = η PCL K ML η PCL α β α(1-N ML ) (5)
This class of functions is able to generate variable concentration features [START_REF] Soby | Viscoelastic Properties of Solutions of Ovine Submaxillary Mucin[END_REF][START_REF] Delgado-Reyes | Flow and Dynamic Viscoelastic Characterization of Non-Purified and Purified Mucin Dispersions[END_REF]. Equation 4 corresponds to a Newtonian fluid η G (γ, α) ≡ η PCL when α = 0 and evolves to the Carreau law when α increase to 1 in the ML:

η ML (γ) = η PCL + η ML -η PCL 1 + (βγ) 2 N ML -1 2 (6)

Numerical method

Solving numerically the equation ( 2) in a mobile and complex geometry, could have been done using various approaches, such as accelerated kernel methods [START_REF] Tornberg | A Fast Multipole Method for the Three Dimensional Stokes Equations[END_REF], immersed boundary methods designed for heterogeneous flows (See [START_REF] Fai | Immersed boundary method for variable viscosity and variable density problems using fast constant-coefficient linear solvers I: Numerical method and results[END_REF][START_REF] Poncet | Analysis of an immersed boundary method for three-dimensional flows in vorticity formulation[END_REF] for its time-dependent formulation), or improved penalization techniques [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF][START_REF] Chatelin | A hybrid grid-particle method for moving bodies in 3D Stokes flow with variable viscosity[END_REF][START_REF] Chatelin | Hybrid grid-particle methods and penalization: A Sherman-Morrison-Woodbury approach to compute 3D viscous flows using FFT[END_REF]]. Among the different numerical methods applied to a fluid bathing a set of cilia, one can find Stokeslets methods [START_REF] Guo | Cilia beating patterns are not hydrodynamically optimal[END_REF][START_REF] Fauci | Biofluidmechanics of Reproduction[END_REF], finite element methods [START_REF] Montenegro-Johnson | Physics of rheologically enhanced propulsion: Different strokes in generalized Stokes[END_REF][START_REF] Smith | A viscoelastic traction layer model of muco-ciliary transport[END_REF][START_REF] Smith | Modelling mucociliary clearance[END_REF] or immersed boundary coupled to lattice Boltzmann methods [START_REF] Jayathilake | Threedimensional numerical simulations of human pulmonary cilia in the periciliary liquid layer by the immersed boundary method[END_REF][START_REF] Li | An immersed boundary-lattice Boltzmann method for single-and multicomponent fluid flows[END_REF], this being not exhaustive.

In this study, we have used the equation ( 2) penalized by ε -1 (uū) over the cilia domain B(t), in order to satisfy the continuity between the cilia and the fluid. This is a nonlinear Stokes problem, analyzed mathematically in [START_REF] Chatelin | Analysis of the penalized 3D variable viscosity Stokes equations coupled to diffusion and transport[END_REF], which is robust and well-fitted to deformable geometry:

-div(2η G (α, u)D(u)) + χ(t) ε (u -ū) = f in Ω, (7)
in the whole computational domain Ω, u(t) is the cilium velocity deduced from the beating model, χ(t) the characteristic function of B(t), and ε 1 the penalization parameter.

In [START_REF] Chatelin | A hybrid grid-particle method for moving bodies in 3D Stokes flow with variable viscosity[END_REF], an iterative projection algorithm was developed to compute separately the pressure and the velocity field of a Newtonian variable viscosity Stokes problem. This numerical method efficiently computes the velocity field close to the boundaries (an inherent problem of classical projection methods [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]). This original methodology allows (in the same iteration) the explicit integration of the nonlinearity of the power law so that the generalized algorithm becomes a multi-criterion fixed point.

First, we introduce a projection on divergence-free fields denoted Π(u * ) = u * -∇ζ where ζ is the solution to the following linear elliptic partial differential equation:

         -∆ζ = -divu * in the domain Ω ∂ζ ∂n = u * • n on the boundary ∂Ω ( 8
)
where n is the outer unit vector normal to ∂Ω. This numerical algorithm relies on fast solvers [START_REF] Sweet | A Parallel and Vector Variant of the Cyclic Reduction Algorithm[END_REF], based on FFT decompositions, guaranteeing a quasi linear computational cost with respect to the number of grid points, which suits very well to large 3D computations. A globally second order numerical solution of Eq. 7 is the limit of the sequence u k = Π(u * k ) where u * k is defined by:

                 -η G (α, u k )∆u * k+1 + χ(t) ε (u * k+1 + Π(u * k ) -u * k -ū) = f + 2D(u k ) + (divu * k )Id ∇η G (α, u k ) in Ω u * k+1 = g -Π(u * k ) + u * k on ∂Ω (9)
The solution obtained by this method satisfies the boundary conditions u = g, which is usually not the case when such a Dirichlet conditions is set together with projection methods [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]. The equation ( 9) can be written

-∆u * k+1 + χ(t) [ε η G ] -1 u * k+1 = RHS (u * k , ū) (10)
where

RHS (u * k , ū) =η -1 G f + χ(t) [ε η G ] -1 (ū + u * k -Π(u * k )) + 2D(u k ) + (divu * k )Id ∇ log η G (α, u k ) (11)
The successive iterations of equation ( 10) have been firstly solved by a multigrid efficient solver [START_REF] Adams | Modified multigrid for 3D elliptic equations with cross derivatives[END_REF][START_REF] Chatelin | A hybrid grid-particle method for moving bodies in 3D Stokes flow with variable viscosity[END_REF]. However, this had been improved by the use of matrix perturbation theory [START_REF] Chatelin | Hybrid grid-particle methods and penalization: A Sherman-Morrison-Woodbury approach to compute 3D viscous flows using FFT[END_REF]: Indeed, the penalization operator χ(t) [ε η G ] -1 Id in (10) exhibits large jumps, but localized, and can be seen as a large perturbation of the -∆ operator. By means of the Sherman-Morrison-Woodbury formula, this elliptic PDE is reduced to a linear system whose size is the number of grid point in the cilia body B(t), and whose any evaluation requires a call to a FFT solver (see [START_REF] Chatelin | Hybrid grid-particle methods and penalization: A Sherman-Morrison-Woodbury approach to compute 3D viscous flows using FFT[END_REF] for more details).

The mucin maturation parameter α can be only a function of distance from the airway wall (called static mixing), such as the following sigmoid expression:

α(x, y, z) = tan -1 (λ(z/H -δ)) -tan -1 (-λδ) tan -1 (λ(1 -δ)) -tan -1 (-λδ) (12)
where is the ASL height is H, PCL height is δH, and the dimensionless transition coefficient is λ. It can also satisfy a transport equation (called convected mixing):

∂ t α + u•∇α -σ∆α = 0 (13)
where Eq. ( 12) then stands for the initial condition, where σ is a numerical smoothing diffusion chosen as small as possible in order to avoid spurious jumps occurring in pure transport (10 -11 m 2 .s -1 for instance). When the convected mixing is used, convection is dominant and a Lagrangian method is used to integrate in time this partial differential equation [START_REF] Chatelin | A hybrid grid-particle method for moving bodies in 3D Stokes flow with variable viscosity[END_REF]. Once it is coupled with the Eulerian 3D Stokes solver, it yields to a robust numerical tool, allowing large time steps without significant loss on accuracy, up to 34-fold the time step satisfying the usual transport stability condition u max ∆t < ∆x.

In the present cases, choosing static or convected mixing, whose resulting α-isovalues are displayed on Fig. 2, leads only to a 3% difference in the mucus mean velocity. Furthermore, the case of Newtonian fluids had been investigated in [START_REF] Chatelin | A parametric study of mucociliary transport by numerical simulations of 3D non-homogeneous mucus[END_REF], with linear Stokes equation. The present study focuses specifically on non-Newtonian effects on the mucociliary clearance, involving a non-linear Stokes model. The resulting non-Newtonian effects and their impact on mucociliary clearance are discussed thereafter.

Phase cartography

Measurements on mucus samples and numerical simulations enable to plot phase cartography. Resulting maps yield positions of healthy and CF mucus from the present study and literature. The effect of PCL height, which markedly impairs the mucociliary clearance, can be investigated.

The balance between PCL and ML quantified by the scalar field α following Eq. ( 12), with a PCL height is chosen as δH where δ = 0.6, for a total ASL height H = 12 µm, representing the proximal compartment of the tracheobronchial tree, and a transition coefficient λ = 5. In the second part of this study, the PCL height varies, δ ranging from 0.3 to 0.7.

Moreover, in order to quantify the efficiency of the mucociliary clearance, the velocity is averaged both in space, over the full ML, and in time, over 5 beating cycles:

U = 1 5T meas(Ω) 5T 0 Ω u x (x, y, z, t) dV dt ( 14
)
where x denotes the streamwise direction given by the cilia beating.

Cystic fibrotic mucus positioning in phase cartography

Because the cilium beating is asymmetrical, this velocity is positive and a net displacement of the ML is observed in the direction of the esophagus. A snapshot of a simulation is presented on Fig. 2. The mean mucus velocity is computed for different values of the {N, K} couple. Simulations were carried out for 169 {N, K} couples to obtain the phase cartography of Fig. 3.

Viscosity values are not fully meaningful for shear-thinning fluids. The consistency index K takes generally lower values 1).

than that in healthy subject samples. The power index varies with a greater extent in CF samples compared to values observed in healthy subject samples (except for HL9).

As depicted in Fig. 3, the mucus transport velocity does not evolve monotonously with neither N nor K, a moderate decrease in the K index leading to a lower mucus transport with a loss of velocity in the range 6-14% for our mucus samples (Table 1). An isolated drop of 6% in mucus mean velocity may suffice to be pathological. Although the HL9 datum [START_REF] Serisier | Macrorheology of cystic fibrosis, chronic obstructive pulmonary disease and normal sputum[END_REF] is very different from others [START_REF] Jeanneret-Grosjean | Sampling technique and rheology of human tracheobronchial mucus[END_REF][START_REF] Zayas | Tracheal mucus rheology in patients undergoing diagnostic bronchoscopy. Interrelations with smoking and cancer[END_REF] (but this might be linked to a nebulization realized prior to the sampling), the mucus velocity magnitude varies less than 5% in these healthy mucus cases. The shear-thinning behavior has a strong effect, as the mucus transport velocity decays to a large extent when N decreases to 0.2. This result is quite counterintuitive, but might result from a larger backflow of the mucus layer when viscosity decreases. This phenomenon depends highly on the respective ML and PCL height that affect the cilium impingement on mucus (strength and duration).

Clinically, two classes of CF mucus can be considered: stabilized and decompensated. An acute secondary infection exacerbates chronic infection with ciliated cell death, further disrupting the mucus propulsion machinery. The expectoration rheology of the CFM3 patient was measured both in stabilized and decompensated states: post-treatment sputum, after 17 days of treatment [CFM3b, stabilized] and 2 days after hospitalization [CFM3a, acute] displayed on Fig. 3 and Table 1. According to measures and computations, the mucociliary clearance improves, as expected using a combined antibiotic and steroid treatment. Figure 4: Loss of mucus velocity with respect to rheological parameters (K, N) for several PCL/ML transition heights δ, given in the bottom right plots. The reference velocity (0%) is the healthy cases from [START_REF] Jeanneret-Grosjean | Sampling technique and rheology of human tracheobronchial mucus[END_REF][START_REF] Zayas | Tracheal mucus rheology in patients undergoing diagnostic bronchoscopy. Interrelations with smoking and cancer[END_REF], also displayed in Fig. 3 and 1 and Table 1.

Impact of PCL thickness

Altered mucus velocity is shown in {K, N} diagrams (Fig. 4) for several heights of periciliary fluid (i.e., several δ values), representing the effect of variable CFTR deficiency [START_REF] Thornton | Structure and Function of the Polymeric Mucins in Airways Mucus[END_REF][START_REF] Ehre | Cystic fibrosis, An inherited disease affecting mucin-producing organs[END_REF]. This alteration is given in percent of the mean healthy data (i.e., HL6 and HL8 in Table 1 [START_REF] Jeanneret-Grosjean | Sampling technique and rheology of human tracheobronchial mucus[END_REF][START_REF] Zayas | Tracheal mucus rheology in patients undergoing diagnostic bronchoscopy. Interrelations with smoking and cancer[END_REF]).

The PCL/ML thickness ratio impacts markedly the mucus transport. The mucus velocity is reduced from -14% to -21% when the PCL thickness decreases (δ falls from 0.6 to 0.3). Hence, for a mucus with a given set of values of the parame-ter couple {N, K}, a weak decrease in δ significantly reduces the mucus velocity. A limitation of this work is linked to the absence of interference between the mucus rheology and motion amplitude and frequency of respiratory cilia for a given state of ciliated cells. However, pathophysiological data are missing to handle this interaction, which can yet be explored using the present model.

Conclusion

The present work had two major objectives. Firstly, it was aimed at using an efficient solver that enables the analysis of all influent shear-thinning parameters which control the mucus motion in the respiratory tract, using a mucus maturation model closer to the biological process than a double-deck fluid domain with sharp PCF-ML transition. This mucus baths beating cilia, through which travel immersed mucin polymers that hydrate and cross-link. Secondly, it was aimed at providing a set of experimental parameters both in order to confront the numerical results (by means of healthy/pathologic classification) and to provide input rheological parameters to the numerical simulations. Despite the complexity and richness of our model, many data that are mandatory to perfectly handle the process are still missing (effect of CF on the respective thickness of the PCL and ML layers, effect of mucus viscosity on cilium beating, etc.). Nevertheless, we have shown that shear-thinning effects can of the same order as the visco-elastic effects in such environments [START_REF] Sedaghat | On the effect of mucus rheology on the mucociliary transport[END_REF].

The highly shear-thinning character of the mucus leads to very high ratio of the dynamic viscosity in the flow, around 10 5 when the shear rate γ evolves in the range [10 -2 -10 3 ]. Therefore the numerical method based on projective iteration [START_REF] Chatelin | A hybrid grid-particle method for moving bodies in 3D Stokes flow with variable viscosity[END_REF][START_REF] Chatelin | Hybrid grid-particle methods and penalization: A Sherman-Morrison-Woodbury approach to compute 3D viscous flows using FFT[END_REF] extended herein to nonlinear viscosity proves to be extremely robust.

Cartographies show a strong influence of the power law index N that must be properly incorporated in numerical simulations. In stabilized cystic fibrosis, the mucus viscosity observed for a given shear rate does not significantly differ from values observed in healthy people. Nevertheless, as the mucus velocity does not evolve monotonously with the rheological indices N and K, mucus transport in CF patients can be less efficient. A vicious circle is then created, as a drop in mucuciliary clearance allow local pathogen proliferation and further deterioration of the respiratory epithelium function.

Figure 2 :

 2 Figure 2: Evolution of isovalues of mucin maturation α ∈ [0.4, 0.9] (in red). Top picture: solution of transport equation (13) coupled to the 3D Stokes equation[START_REF] Button | Periciliary brush promotes the lung health by separating the mucus layer from airway epithelia[END_REF]. Bottom picture: static mixing given by Eq. (12).

Figure 3 :

 3 Figure 3: Phase cartography of the mucus transport velocity, with respect to rheological parameters {N, K}. Collected samples are displayed in red, while the results collected from literature are displayed in blue (Table1).

Table

  and Table 1). Three examples of steadystate shear viscosity curves for different mucus from two to

	Alias	N	K Source
	CFL1 0.408 0.544 [54]
	CFL2 0.182 1.397 [61]
	CFL3 0.166 2.245 [61]
	CFL4 0.163 2.423 [61]
	CFL5 0.17 7.34 [15]
	HL6 0.19 14.1 [34]
	HL7 0.195 15.49 [34]
	HL8 0.215 13.18 [67]
	HL9 0.482 0.355 [54]
	CFM1 0.20 5.80
	CFM2 0.39 0.45
	CFM3 0.32 1.24
	CFM4 0.39 1.60 CFM5 0.30 2.37 CFM6 0.23 4.22 CFM7 0.21 3.10 CFM8 0.23 3.53 CFM9 0.37 0.51	Present study
	CFM10 0.36 2.17
	CFM11 0.22 4.75
	BM12 0.40 1.43

Table 2 :

 2 Parameters used for the numerical simulations, gathered from[START_REF] Sanderson | Ciliary activity of cultured rabbit tracheal epithelium, Beat pattern and metachrony[END_REF][START_REF] Smith | Modelling mucociliary clearance[END_REF][START_REF] Fahy | Airway mucus function and dysfunction[END_REF][START_REF] Thiriet | Tissue Functioning and Remodeling in the Circulatory and Ventilatory Systems[END_REF] or from our own measures.

	Name	Value
	Cilia length	8 µm
	Cilia diameter	0.3 µm
	Number of Cilia	25
	PCL viscosity (η PCL )	0.001 Pa.s
	Viscosity transition stiffness (λ)	5
	Transition length (δ)	0.6
	ASL height (H)	12 µm
	Beating frequency	10 Hz
	Inverse low-frequency cut-off (β)	4 10 3 s
	Distance between adjacent cilia	1.5 µm
	Metachronal wavelength	87.5 µm

). The penalty

Acknowledgments

This work was partially supported by the ANR Grant BioFiReaDy, under the contract number ANR-2010-JCJC-0113-01. This work was also granted access to the HPC resources of P2CHPD (from Fédération Lyonnaise de Simulation et de Modélisation Numériques). We also acknowledge the use of the rheometer from the cooperative structure of joint research FERMaT (FR3089).

Conflicts of interest declaration

We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome. We confirm that the article has been read and approved by all named authors. We further confirm that the order of authors listed in the paper has been approved by all of us.


Florian Yger
email: florian.yger@dauphine.fr

Sylvain Chevallier
email: sylvain.chevallier@uvsq.fr

Quentin Barthélemy
email: q.barthelemy@foxstream.fr

Sinno Jialin Pan

Masashi Sugiyama

S Sra

Geodesically-convex optimization for averaging partially observed covariance matrices

Keywords: SPD matrices, average, missing data, data imputation

Symmetric positive definite (SPD) matrices permeates numerous scientific disciplines, including machine learning, optimization, and signal processing. Equipped with a Riemannian geometry, the space of SPD matrices benefits from compelling properties and its derived Riemannian mean is now the gold standard in some applications, e.g. brain-computer interfaces (BCI). This paper addresses the problem of averaging covariance matrices with missing variables. This situation often occurs with inexpensive or unreliable sensors, or when artifact-suppression techniques remove corrupted sensors leading to rank deficient matrices, hindering the use of the Riemannian geometry in covariance-based approaches.

An alternate but questionable method consists in removing the matrices with missing variables, thus reducing the training set size. We address those limitations and propose a new formulation grounded in geodesic convexity. Our approach is evaluated on generated datasets with a controlled number of missing variables and a known baseline, demonstrating the robustness of the proposed estimator. The practical interest of this approach is assessed on real BCI datasets. Our results show that the proposed average is more robust and better suited for classification than classical data imputation methods.

Introduction

With from the pioneering work of Candès and Recht (2009), a new light was shed on the matrix completion problem and the imputation of missing data. It demonstrated that data imputation could be formulated as an optimization problem benefiting from the underlying structure of the data. Such method has been successfully applied to recommendation systems (Rennie and Srebro, 2005), multitask learning (Obozinski et al., 2010) or manifold learning (Weinberger and Saul, 2006). The bound of the number of entries for an optimal reconstruction of a low-rank matrix are provided in Candès and Recht (2009) and improved in Recht (2011), with only minimal assumptions on the coherence of the matrix to recover.

Several problems connected to low-rank matrix completion have been considered, including recovering spectrally sparse objects from missing time-domain samples (Chen and Chi, 2013), or in graph data (Narang et al., 2013). Positive definite and positive semi-definite matrices are indeed good candidates for matrix completion approaches as those matrices lives in a low-dimensional space, and have been the subject of early works in matrix completion (Johnson, 1990;Johnson and Tarazaga, 1995). Symmetric positive-definite (SPD) matrices are commonly found in covariance-based approaches, such as financial time series processing (Bingham, 2014), radar detection [START_REF] Arnaudon | Riemannian medians and means with applications to radar signal processing[END_REF], computer vision (Harandi et al., 2017) or medical imaging (Pennec et al., 2006).

In this contribution, we focus on yet another application of positive-definite matrices: the brain-computer interfaces (BCI). Since the seminal work of [START_REF] Barachant | Multiclass brain-computer interface classification by Riemannian geometry[END_REF] that introduced a new approach for BCI, the covariance matrices are directly handled on the manifold of positive-definite matrices. This implied to reformulate common machine learning algorithms and processing pipelines for this specific geometry. These approaches became the new gold standard for BCI (Yger et al., 2017;Congedo et al., 2017) with high accuracy on various kinds of brain signals, such as motor imagery (Jayaram and Barachant, 2018), evoked potentials (Korczowski et al., 2015) or steady-state visually-evoked potential (Chevallier et al., 2020), and won several Kaggle competitions.

In all these applications, from BCI to financial time series, machine learning models rely on average centers of SPD matrices. The estimation of average (or center of gravity) is the key to obtain reliable results for the target application, and must be robust to missing variables. However, in BCI for example, unreliable sensors or artifact-suppression techniques removing corrupted sensors leading to rank deficient matrices, hindering the use of the usual Riemannian averaging. To our knowledge, there is no systematic review of the processing of SPD matrices with missing data in the literature. This is the main objective of this paper, by confronting the existing approaches and by introducing a new framework.

Our contributions are the following: i) a sound framework for SPD matrices averaging with missing variables, and a theoretical analysis of the geodesic convexity of this optimization,

ii) an evaluation of the convergence of the proposed approach on a synthetic dataset, iii) a comparison with Euclidean mean and common imputation methods in machine learning, iv) a validation on a real dataset for BCI application.

Review of related works

There are several types of missing data that could be encompassed in the literature. In an attempt to make a summary of the studies on the subject, we could describe the different cases that are illustrated on Fig. 1. Our approach comes from practical concerns occurring when considering multivariate data. We call X ∈ R c×n a matrix composed of c input variables (sensors, or components, shown in rows) and of n samples (observations, visible in columns). We call Σ ∈ R c×c its covariance matrix, which is SPD1 . As discussed hereafter, the data can be corrupted and missing data can occur with particular patterns.

The first case one occurs when some multivariate samples have not been acquired by the system, resulting in missing columns in the matrix X. Missing values in this case is often interpolated or imputed from average or median values. In the second case, some sensors have a trouble during acquisition. Consequently, some variables/rows are missing in X, and these rows are thus missing in Σ, but also their corresponding columns due to the symmetry.

In Fukuda et al. (2001), the authors rely on matrix completion techniques to propose a solver in semi-definite programming that is able to cope with partial information in the constraint space of the primal-dual interior-point methods. Tsuda et al. (2003) propose to complete missing variables for SVM kernel estimation in genetic applications, relying on an auxiliary kernel matrix that integrates another source of information.

Some works consider the case where some elements of the matrix are missing (Bishop and Byron, 2014;Choi et al., 2019), but without a particular link between these elements, excepted symmetry. The work of Johnson and Tarazaga (1995) is on matrix completion for graph data with positive semi-definite constraint, for missing entries outside the diagonal. Another approach by Laurent (2001), using the same setup, try to estimate missing information on the diagonal.

Several problems can be stated when processing incomplete matrices, i.e. matrices with missing data. The first one concerns the estimation of covariance from an incomplete matrix X, with the objective to build a robust estimator (Little, 1998;Schneider, 2001;Lounici, 2014). The second category of problems is the covariance completion. From an incomplete covariance matrix Σ, the goal is to find the best way to complete its missing values (Johnson and Tarazaga, 1995;Bishop and Byron, 2014).

A two-step approach is presented in (Rodrigues et al., 2019): firstly, covariance matrices are imputed on the manifold (performing on the matrix a whitening, a scaling, completing it by adding 1s on the diagonal where variables are missing, and then de-whitening), secondly a classical covariance averaging is applied on completed matrices. From several incomplete covariance matrices, the goal is to estimate the best average covariance matrices. In the work of Choi et al. ( 2019), the mean is limited to 2 covariance matrices, having identical missing elements. To sum up this state-of-the-art, all these works are classified in Table 1 along two criteria: the type of missing data (horizontal) and the type of problems (vertical). The contribution of this paper belongs to covariance averaging, providing a way to average several SPD matrices with missing variables. It is important to note that missing variables can change for each matrix to average. The resulting matrix is defined on the union (and not the intersection) of variables of input matrices.

Estimating Riemannian mean with missing variables

After a brief recall of the Riemannian geometry to process SPD matrices, this section presents the problem of estimating an average matrix from a set of SPD matrices with missing variables, introducing masks to tackle missing variables. The cost function is shown to be geodesically convex under some restricted cases.

Preliminaries

We consider the set P c of symmetric positive definite (SPD) of size c × c. Formally, this set is defined as:

P c = Σ ∈ R c×c | Σ = Σ , ∀x ∈ R c x Σx > 0 , (1)
and it is endowed with the Löwner partial order, defined between Σ 1 , Σ 2 ∈ P c as:

Σ 1 Σ 2 ⇔ Σ 1 -Σ 2 ∈ P c .
The most straightforward way to handle P c is to consider it as a subspace of S c , the space of symmetric matrices. As such, it is a Euclidean space and the tools of conic geometry and semi-definite programming could apply. However, this Euclidean geometry comes at the cost of dramatic drawbacks, like the swelling effect 2 (Fletcher et al., 2004; Horev et al., 2015). Non-Euclidean alternatives have been successfully applied in some applications (Yger et al., 2017;Congedo et al., 2017;Harandi et al., 2017).

The most common estimator of the covariance matrix Σ is the sample covariance matrix, defined as:

Σ = 1 n -1 XX .
(2)

More robust estimators have been proposed (Daniels and Kass, 2001) and have been evaluated in the context of EEG-based BCI (Chevallier et al., 2018). The most common are the shrinkage estimator that combines the sample covariance matrix with a target covariance matrix, that could be similar to the identity matrix, to ensure that the estimated matrices avoid any ill-conditioning problem.

Due to the curvature of the space of SPD matrices, an adequate metric could rely on Σ 1 t Σ 2 , the geodesic between Σ 1 , Σ 2 ∈ P c and defined as:

Σ 1 t Σ 2 = Σ 1 2 1 Σ -1 2 1 Σ 2 Σ -1 2 1 t Σ 1 2 1 = Σ 1 2 1 Exp t Log Σ -1 2 1 Σ 2 Σ -1 2 1 Σ 1 2 1 , t ∈ [0, 1] . (3)
For clarity reason, when t is omitted Σ 1 Σ 2 = Σ 1 0.5 Σ 2 , the geometric center of Σ 1 and Σ 2 . We will consider hereafter the affine-invariant Riemannian (AIR) distance, that derived from the geodesic and that is widely used for its geometric properties.

Definition 1 Affine-invariant Riemannian distance (Förstner and Moonen, 1999). For Σ 1 , Σ 2 ∈ P c , two SPD matrices, the Riemannian distance is defined as

δ R (Σ 1 , Σ 2 ) = Log(Σ -1 2 1 Σ 2 Σ -1 2 1 ) F = c k=1 log 2 λ k 1 2 , ( 4
)
where • F is the Frobenius norm, and

λ k , k = 1, . . . , c, are the eigenvalues of Σ -1 2 1 Σ 2 Σ -1 2 1 .
Among the numerous properties of the affine-invariant Riemannian distance [START_REF] Bhatia | Positive definite matrices[END_REF], the congruence-invariance3 states that, for any invertible matrix

W ∈ R c×c : δ R (Σ 1 , Σ 2 ) = δ R (W Σ 1 W, W Σ 2 W ) . (5)
We now want to estimate the mean of a set of m SPD matrices {Σ i } m i=1 . It may actually be helpful to estimate a weighted mean, where the weights {w i } m i=1 are non-negative and sum up to a constant (in practice, 1 or m). These weights can be fixed to encode prior knowledge, or tuned by cross-validation to optimize some given criterion.

Definition 2 Riemannian mean (Moakher, 2005;Fletcher et al., 2004). The Riemannian mean Σ ∈ P c of a set of m SPD matrices {Σ i } m i=1 , associated to non-negative weights {w i } m i=1 , is defined as

Σ = arg min Σ∈Pc f R (Σ) = arg min Σ∈Pc 1 2 m i=1 w i δ 2 R (Σ i , Σ) . (6)
This mean, also known as geometric mean, has no closed form and therefore has to be computed iteratively, for example through a gradient descent algorithm. The gradient of the cost function f R is defined as (Fletcher et al., 2004):

∇f R (Σ) = - m i=1 w i Σ 1 2 Log(Σ -1 2 Σ i Σ -1 2 )Σ 1 2 = - m i=1 w i Log Σ (Σ i ) , (7)
where Log Σ (Σ i ) is the logarithmic map, projecting the matrix Σ i from manifold to tangent space at point Σ. The geometric mean is the unique matrix Σ of the manifold P c nullifying the sum of the m tangent vectors.

As defined in Zadeh et al. ( 2016) and explained in Boumal (2020, Chap 11), the geometric mean is a geodesically convex function. Even if it is not convex in the Euclidean sense, geodesic convexity ensures that functions are convex along the manifold and enjoying similar properties to convex functions.

Definition 3 Geodesically convex function (Sra and Hosseini, 2015). Let P be a Riemannian manifold and Q ⊂ P a geodesically convex set. A function f : Q → R is said geodesically convex, if for all points Σ 1 and Σ 2 of this set, it satisfies

f (Σ 1 t Σ 2 ) ≤ tf (Σ 1 ) + (1 -t)f (Σ 2 ) , t ∈ [0, 1] .
(8)

Masked Riemannian mean

A mask M ∈ R c×(c-p) is defined as the identity matrix, but with the columns of the p missing variables removed, with 0 ≤ p < c. In our problem, the mask is not constrained to be identical for all input matrices, ie. each input matrix can be associated to a particular mask.

If an input X ∈ R c×n is incomplete, with p missing variables, the available (non-missing) variables can be extracted thanks to the product by a mask M ∈ R c×(c-p) . The resulting complete submatrix is X = M X ∈ R (c-p)×n . Thus, its covariance matrix Σ ∈ P c-p estimated with Eq. ( 2) is:

Σ = 1 n -1 X X = 1 n -1 M XX M = 1 n -1 M ΣM , (9)
where Σ = XX is an incomplete covariance matrix. Incomplete matrices X and Σ are illustrated in Fig. 1(b). Consequently, a mask is a full column-rank matrix such that the submatrix M ΣM is SPD, of size (c -p) × (c -p) 4 . Now, we consider the set {Σ i } m i=1 of incomplete SPD matrices, ie. each matrix Σ i has p i missing variables, with 0 ≤ p i < c. Consequently, their associated masks {M i } m i=1 are potentially all different. The masked Riemannian mean Σ ∈ P c is defined as:

Σ = arg min Σ∈Pc f R M (Σ) = arg min Σ∈Pc 1 2 m i=1 w i δ 2 R (M i Σ i M i , M i ΣM i ) . ( 10
)
Note that, each matrix M i Σ i M i being SPD, the distance δ R can be applied and the cost function is well defined. As explained in [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF] and Boumal (2020, Chap 4), the Riemannian gradient can be deduced from the Euclidean gradient and then plugged into a Riemannian gradient descent. Combining the chain rule of derivatives with the gradient formula given in Eq. ( 7), the Euclidean gradient of the cost function f R M (Σ) can be written as:

∇f R M (Σ) = - m i=1 w i M i Log M i ΣM i M i Σ i M i M i . (11)
The cost being computed on subparts of the matrix Σ, the gradient is null for the elements not involved on the compressed matrix M i ΣM i . Hence, the Euclidean gradient of δ 2 R , although computed on a subpart of Σ can be broadcast to the whole space using M i (•)M i operator. Remark that Eq. ( 11) falls back to Eq. ( 7) when there is no missing variables, i.e. when all masks M i = I c .

The illustration of the averaging of two masked matrices is displayed in Fig. 2.

Link with the NaN-mean

If we replace the Riemannian distance δ R by the Euclidean distance δ E (Σ 1 , Σ 2 ) = Σ 1 -Σ 2 F into Eq. ( 10), the gradient of the masked Euclidean mean would have been:

∇f E M (Σ) = - m i=1 w i M i M i Σ i M i -M i ΣM i M i = - m i=1 w i M i M i (Σ i -Σ) M i M i , (12)
where the operator M i M i (•)M i M i extracts only the non-missing values of Σ i .

Defining I u,v as the set containing the indices of matrices without missing values for entry (u, v), the gradient can be expressed element-wise:

∇f E M (Σ)(u, v) = - i∈Iu,v w i (Σ i (u, v) -Σ(u, v)) . ( 13
)
This gradient is nullified by the matrix:

Σ(u, v) = i∈Iu,v w i Σ i (u, v) i∈Iu,v w i . ( 14
)
When all weights are equal, this masked Euclidean mean is commonly called NaN-mean and is written as:

Σ(u, v) = 1 |I u,v | i∈Iu,v Σ i (u, v) , (15)
where |I u,v | denotes the cardinality of I u,v .

Σold

Σnew Σ 2 ) containing p 1 (resp. p 2 ) missing variables encoded by mask M 1 (resp. M 2 ).

P c P c-p 1 M 1 (•)M 1 M 1 (•)M 1 M 1 Σold M 1 M 1 Σ 1 M 1 Log M 1 Σold M 1 (M 1 Σ 1 M 1 ) P c-p 2 M 2 (•)M 2 M 2 (•)M 2

Geodesic convexity

In practice, we noticed that a simple Riemannian gradient descent always converged to the same minimum (independently of its initialization) on our experiments and we tested it numerically in Pymanopt (Townsend et al., 2016). This lead us to formulate a conjecture on the geodesic convexity of the problem. Unfortunately, we could generate some counterexamples to this conjecture in simplistic cases, but we also found some restricted cases where the geodesic convexity holds.

Applying the congruence-invariance of Eq. ( 5) with

W = (M i Σ i M i ) -1 2 ∈ R (c-p i )×(c-p i )
, we can express the cost function f R M , defined in Eq. ( 10), as follows:

f R M (Σ) = 1 2 m i=1 w i δ 2 R I (c-p i ) , ((M i Σ i M i ) -1 2 ) M i ΣM i (M i Σ i M i ) -1 2 = 1 2 m i=1 w i δ 2 R I (c-p i ) , M i Σ Mi , (16)
where

Mi = M i (M i Σ i M i ) -1 2 ∈ R c×(c-p i )
is a full column-rank matrix.

In the case where we have M M I (c-p) , then, for any Σ 1 , Σ 2 I c :

I (c-p) M M M (Σ 1 Σ 2 ) M ( M Σ 1 M ) ( M Σ 2 M ) , (17)
and because Log(•) 2 F = δ 2 R (I, •) is monotonically increasing above I, we have:

δ 2 R I (c-p) , M (Σ 1 Σ 2 ) M ≤ δ 2 R I (c-p) , ( M Σ 1 M ) ( M Σ 2 M ) . (18)
From this, g-convexity follows upon using g-convexity of δ R (Sra and Hosseini, 2015) in the usual way for the very last occurrence of δ R . This proof is detailed in the supplementary material of this article.

In the case where we have M i Mi I (c-p i ) , ∀i, then, each element of the sum in f R M (Σ) is g-convex on Σ on the set of SPD matrices such that Σ i I c , where matrices Σ i ∈ P c ideally exist but are only partially observed.

Experimental analysis

The first experiment evaluates the proposed approach on an artificial dataset and compares it with an Euclidean approach. The second experiment demonstrates the advantage of averaging incomplete matrices over the classical methods of data imputation, that are the matrix deletion and channel (or variable) deletion. The third experiment, provided in the supplementary material of this article, evaluates the convergence rate of the new algorithm. The last experiment is conducted on a real dataset, acquired for a motor imagery experiment in BCI. For all experiments, we used uniform weights w i for the matrices (hence leading to un-weighted averages). Experiments are made using pyRiemann 0.2.6 and Pymanopt 0.2.5 (Townsend et al., 2016).

Evaluation of convergence on synthetic dataset

For the first experiment, the dataset is defined such that the training samples are generated according to a distribution parametrized by dispersion σ and reference G. To this end, the reference matrix is generated as G = U DU , where the diagonal matrix D ∈ R c×c has strictly positive values drawn from a triangular distribution and where the orthogonal matrix U is obtained as the eigenvectors of a random matrix A ∈ R c×2c . The sample matrices of the dataset {Σ i } m i=1 are built as Σ i = U (D + i )U , where i is an additive white noise drawn from a Gaussian distribution N (0, σ 2 ). This perturbation is introduced to control the dispersion of matrices on the manifold during generation. Indeed, if lim m→∞ Ḡ = G, where Ḡ is the geometric mean of {Σ i } m i=1 , in practice it is equal to Ḡ = U (D + ¯ )U , with ¯ = 1 m i i . Then, the masks {M i } m i=1 are generated with p missing variables. In these experiments, m = 1000, c = 10, and p = 1, 2, 3. The masked Riemannian mean RM of incomplete matrices is computed by solving the optimization problem Eq. ( 10). The masked Euclidean mean ĒM is computed using Eq. ( 14). Visualisations of this dataset are displayed in the supplementary material of this article.

To evaluate of convergence of the masked Riemannian mean RM , we compute its distance to the groundtruth mean Ḡ. This distance is compared to the distance between the masked Euclidean mean ĒM and the groundtruth. To have a good overview, the comparison is repeated 10 times, computing the mean and standard deviations of distances. Moreover, increasing the standard deviation σ of additive noises j allows to control the dispersion of matrices. Different number of missing variables p are also tested in this experiment. Figure 3 compares the distances δ 2 R ( RM , Ḡ) and δ 2 R ( ĒM , Ḡ) when varying parameters σ and p. We observe that all distances increase with dispersion of matrices. Moreover, the masked Riemannian mean is always closer to the groundtruth than the Euclidean one. This result is extremely consistent across all tested parameters, showing that the masked Riemannian mean captures more information allowing to converge to the groundtruth faster than masked Euclidean one. This difference tends to decrease when p increases. We can conclude that the masked mean is more robust to missing data than its Euclidean counterpart, and that the convergence is effective even in case where only a fraction of the information is available.

Comparison with other strategies on synthetic dataset

When some variables are missing or have been removed after a trimming step (outlier rejection of very noisy variables), several methods of machine learning exist for data imputation. Notably, the matrix deletion consists of removing a matrix if any single variable is missing [START_REF] Allison | Missing data[END_REF]. Another common approach is to remove each variable that is missing in some of the matrices, called variable deletion.

Indeed, when the number of matrices with missing variables is small, matrix deletion is a good choice, that could even act as some sort of regularization. When the proportion of matrices with missing variables increases, the variable deletion may be a good strategy as it ensures all matrices are considered (even if it is with less information in each). The objective of this experiment is to compare the Riemannian masked mean RM to the Riemannian mean applied after a matrix deletion strategy, denoted Rmd , and after a variable deletion strategy, denoted Rvd .

We use the same data generation as previous experiment, using m = 1000, c = 10 and σ = 0.5. However, in this experiment only a given proportion of matrices (from 10 to 60 %) have missing variables and are associated with a mask. Distances between the groundtruth and the different means are computed, and the comparison is repeated 10 times. As previously, different values of missing variables are tested. Figure 4 compares the distances δ 2 R ( RM , Ḡ), δ 2 R ( Rmd , Ḡ) and δ 2 R ( Rvd , Ḡ) when varying the number of missing variables p and the proportion of masked matrices. This experiment shows that the Riemannian masked mean is close to the groundtruth and better than the other strategies, matrix deletion and variable deletion. This indicates that masked mean captures more information than other strategies, and that this information helps to converge to the groundtruth. As observed before, this difference tends to decrease when p increases.

We can conclude that masked mean is a good alternative to the deletion strategies, allowing to keep more statistical power in analysis (Olinsky et al., 2003) and avoiding to introduce a bias like the matrix deletion when missing variables are not random.

Validation on BCI dataset

This BCI experiment is conducted on a motor imagery dataset, acquired by Yi et al. (2014), where 10 subjects are performing a left hand/right hand movement and that should be detected using only their brain waves recorded with EEG. There are m = 160 trials, 80 for left hand movement and 80 for right hand, recorded for each subject on c = 60 channels/sensors.

Following the previous experiments, we generate random masks M i to simulate missing channels, and apply them to covariance matrices Σ i computed from EEG trials X i . Then, we compare different means: a masked Riemannian mean RM trained on all matrices of the training set, a Riemannian mean Rmd trained only on the complete matrices (i.e. after applying a matrix deletion strategy on the raw training set), and a masked Euclidean mean ĒM . We rely on minimum-distance-to-mean (MDM) [START_REF] Barachant | Multiclass brain-computer interface classification by Riemannian geometry[END_REF] to embed these different means on a same classifier. This is a simple, robust, and online classifier that is part of the state-of-the-art for motor imagery BCI (Jayaram and Barachant, 2018). In this experiment, the gold standard is defined as the MDM trained on non-masked training set (with access to all information, with no masked or deleted variables). This gold standard thus indicates the results that could be obtained in optimal condition, without missing data. We compare the classification accuracies obtained with a 5-fold cross-validation for each subject, while varying the proportion of masked matrices, from 0 to 60% of matrices with missing variables. The classifiers are evaluated in the two cases: when p = 24 out of the c = 60 EEG channels are masked (thus 40% of the channels) and when p = 36 channels are masked (60% of channels). These numerical experiments encompass the different situations where there is a moderate corruption in few trials to the case where the majority of the dataset is corrupted. The MDM classifier that relies on masked Riemannian mean RM is called " RM +MDM". It is compared to MDM based on matrix deletion strategy, called " Rmd +MDM", and to MDM based on masked Euclidean mean, called " ĒM +MDM".

Accuracies are displayed on Fig. 5, with mean and standard deviation computed across subjects. The results clearly show that RM +MDM outperforms other classifiers. ĒM +MDM is not performing well, even while using a nan-mean is a common strategy during preprocessing. Euclidean MDM is known to perform poorly when compared to Riemannian MDM (Kalunga et al., 2015), as it is subject to the swelling effect. The matrix deletion strategy Rmd +MDM, as it uses only a subpart of the available information, achieves a better accuracy than ĒM +MDM, but is lower than RM +MDM. RM +MDM yields results that are almost at the level of the gold standard, and seems not sensitive to the proportion of matrices with missing channels.

Conclusion

This paper describes a method for the robust estimation of an average matrix from a set of SPD matrices with missing variables. This masked mean method could cope with any number and distribution of missing variables on the considered set. The proposed method solves an optimization problem formulated as a matrix factorization, where the missing variables are discarded to estimate the solution on a collection of embedded Riemannian submanifold. While the convergence of this method is attested in practice, we have proved the geodesic convexity of the problem on a restricted class of SPD matrices. The experiments on synthetic and real datasets exhibit the practical interest of this method. The applications target covariance matrices, but without loss of generality this method could be applied on any set of SPD matrices.

The applicative contributions of the proposed method are multiple: (i) it is a true alternative to data imputation, avoiding to exclude data from analysis or to complete data with mean or median values (Experiments 4.1 and 4.2), (ii) it gives a more robust calibration of Riemannian classifiers, to obtain better classification scores (Experiment 4.3), and (iii) it provides a more robust covariance average for spatial filters estimation (Yger et al., 2015), allowing a supervised dimension reduction for high-dimensional signals enhancing the separation of classes [START_REF] Barachant | MEG decoding using riemannian geometry and unsupervised classification[END_REF].

This method is flexible and paves the way to the extension of several other Riemannian approaches. For example, it could be used as a basis for extending SPD networks (Huang and Van Gool, 2017;Brooks et al., 2019) or dictionary of SPD matrices (Cherian and Sra, 2017) with missing data. Moreover, using a similar matrix decomposition, the proposed approach could go beyond missing variables as we are currently investigating the case of missing subspaces. A following work could address the data imputation problem with a compressive sensing flavor for SPD matrices.

Figure 1 :

 1 Figure 1: Illustration of the different types of missing data: (a) missing samples / observations in matrix X, (b) missing variables in matrices X and Σ, (c) missing elements in the matrix Σ.

Figure 2 :

 2 Figure 2: Illustration of the masked Riemannian mean Σ ∈ P c of two matrices, Σ 1 (resp.Σ 2 ) containing p 1 (resp. p 2 ) missing variables encoded by mask M 1 (resp. M 2 ).

Figure 3 :

 3 Figure 3: Convergence error for the masked Riemannian mean RM (in dark blue) and the masked Euclidean mean ĒM (light green). Distances to groundtruth Ḡ are evaluated with 10 repetitions, for p = {1, 2, 3} missing variables and for σ ∈ [0.05, 0.75].

Figure 4 :

 4 Figure 4: Riemannian mean estimation with different strategies: masked Riemannian mean RM , Riemannian mean computed after matrix deletion Rmd , and Riemannian mean after variable deletion Rvd . Distances to groundtruth Ḡ are evaluated with 10 repetitions, for p = {1, 2, 3} missing variables and for different proportions of matrices with missing variables (10 to 60 %).

Figure 5 :

 5 Figure 5: Comparison of MDM computed on non-masked matrices (gold standard, or R+MDM), MDM based on masked Riemannian mean ( RM +MDM), MDM based on Riemannian mean after matrix deletion ( Rmd +MDM), and MDM based on masked Euclidean mean ( ĒM +MDM). Left: Accuracy averaged for the 10 subjects, with 0 to 60% of matrices having 40% of masked channels (missing variables). Right: Same for 60% of masked channels.

  NH Bingham. Modelling and prediction of financial time series. Lu, S M Pizer, and S Joshi. Principal geodesic analysis for the study of nonlinear statistics of shape. IEEE Trans Med Imaging, 23(8):995-1005, 2004. W Förstner and B Moonen. A metric for covariance matrices. Technical report, Bonn University, 1999. M Fukuda, M Kojima, K Murota, and K Nakata. Exploiting sparsity in semidefinite programming via matrix completion I: General framework. SIAM J Optim, 11(3):647-674, 2001. M Harandi, M Salzmann, and R Hartley. Dimensionality reduction on SPD manifolds: The emergence of geometry-aware methods. IEEE Trans Pattern Anal Mach Intell, 40(1): 48-62, 2017. I Horev, F Yger, and M Sugiyama. Geometry-aware principal component analysis for symmetric positive definite matrices. In ACML, pages 1-16, 2015. Z Huang and L Van Gool. A Riemannian network for SPD matrix learning. In AAAI CAI, pages 2036-2042, 2017. V Jayaram and A Barachant. MOABB: trustworthy algorithm benchmarking for BCIs. J Neural Eng, 15(6):066011, 2018. C R Johnson. Matrix completion problems: a survey. In Matrix Theory and Applications, volume 40, pages 171-198, 1990. C R Johnson and P Tarazaga. Connections between the real positive semidefinite and distance matrix completion problems. Linear Algebra Appl, 223:375-391, 1995. E K Kalunga, S Chevallier, Q Barthélemy, K Djouani, Y Hamam, and E Monacelli. From Euclidean to Riemannian means: Information geometry for SSVEP classification. In GSI, pages 595-604. Springer, 2015. L Korczowski, M Congedo, and C Jutten. Single-trial classification of multi-user p300-based brain-computer interface using Riemannian geometry. In IEEE EMBC, pages 1769-1772, 2015. M Laurent. Polynomial instances of the positive semidefinite and Euclidean distance matrix completion problems. SIAM J Matrix Anal Appl, 22(3):874-894, 2001. R J A Little. Robust estimation of the mean and covariance matrix from data with missing values. J R Stat Soc Ser C Appl Stat, 37:23-38, 1998. K Lounici. High-dimensional covariance matrix estimation with missing observations. Bernoulli, 20(3):1029-1058, 2014. M Moakher. A differential geometric approach to the geometric mean of symmetric positivedefinite matrices. SIAM J Matrix Anal Appl, 26(3):735-747, 2005. S K Narang, A Gadde, and A Ortega. Signal processing techniques for interpolation in graph structured data. In IEEE ICASSP, pages 5445-5449, 2013. G Obozinski, B Taskar, and M I Jordan. Joint covariate selection and joint subspace selection for multiple classification problems. Stat Comput, 20(2):231-252, 2010. A Olinsky, S Chen, and L Harlow. The comparative efficacy of imputation methods for missing data in structural equation modeling. Eur J Oper Res, 151(1):53-79, 2003.

	Commun Stat Theory
	Methods, 43(7):1351-1361, 2014.
	W E Bishop and M Y Byron. Deterministic symmetric positive semidefinite matrix com-
	pletion. In Adv Neural Inf Process Syst, pages 2762-2770, 2014.
	N Boumal. An introduction to optimization on smooth manifolds, 2020. URL http:
	//www.nicolasboumal.net/book.
	D Brooks, O Schwander, F Barbaresco, J-Y Schneider, and M Cord. Riemannian batch
	normalization for SPD neural networks. In Adv Neural Inf Process Syst, pages 15463-
	15474, 2019.
	E J Candès and B Recht. Exact matrix completion via convex optimization. Found Comut
	Math, 9(6):717, 2009.
	Y Chen and Y Chi. Spectral compressed sensing via structured matrix completion. arXiv
	preprint arXiv:1304.4610, 2013.
	S Chevallier, EK Kalunga, Q Barthélemy, and E Monacelli. Review of Riemannian distances
	and divergences, applied to SSVEP-based BCI. Neuroinformatics, pages 1-14, 2020.
	H Choi, S Kim, and Y Shi. Geometric mean of partial positive definite matrices with missing
	entries. Linear Multilinear A, pages 1-26, 2019.
	M Congedo, A Barachant, and R Bhatia. Riemannian geometry for EEG-based brain-
	computer interfaces; a primer and a review. Brain-Computer Interfaces, 4(3):155-174,
	2017.
	M J Daniels and R E Kass. Shrinkage estimators for covariance matrices. Biometrics, 57
	(4):1173-1184, 2001.
	P T Fletcher, C

A Cherian and S Sra. Riemannian dictionary learning and sparse coding for positive definite matrices. IEEE Trans Neural Netw Learn Syst, 28(12):2859-2871, 2017. S Chevallier, EK Kalunga, Q Barthélemy, and F Yger. Brain Computer Interfaces Handbook: Technological and Theoretical Advances, chapter 19 -Riemannian classification for SSVEP based BCI : offline versus online implementations, pages 371-396. 2018. X Pennec, P Fillard, and N Ayache. A Riemannian framework for tensor computing. Int J Comput Vis, 66(1):41-66, 2006.

In theory, covariance matrices belong to the set of symmetric semi-definite matrices but making the realistic assumption that X is corrupted by a Gaussian noise leads to Σ ∈ Pc.

In the Euclidean geometry of definite positive matrices, the determinant of Σ 1 +Σ 2 2 can be greater than the determinants of Σ1 or Σ2. Hence, the Euclidean interpolation can bring some spurious information, contrary to the Riemannian geometry.

Using this transformation as a whitening, this invariance could be used in a transfer learning setup as inYger and Sugiyama (2015).

Hence, mask M is a semi-orthogonal matrix, i.e. M M = I (c-p) , and our setup could be generalized to missing subspaces in a compressive sensing manner using orthogonal matrices.

Acknowledgments

FY acknowledges the support of the ANR as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute). SC acknowledges that this work is not supported by any ANR or IDEX project but by the university recurrent funding.

B Recht. A simpler approach to matrix completion. J Mach Learn Res, 12:3413-3430, 2011