Rémi Monasson

Inference of principal components of noisy correlation matrices with prior information

Keywords: Random Matrix Theory, Spiked Covariance Model, Prior Information, Replica Method, Phase Transitions

The problem of infering the top component of a noisy sample covariance matrix with prior information about the distribution of its entries is considered, in the framework of the spiked covariance model. Using the replica method of statistical physics the computation of the overlap between the top components of the sample and population covariance matrices is formulated as an explicit optimization problem for any kind of entry-wise prior information. The approach is illustrated on the case of top components including large entries, and the corresponding phase diagram is shown. The calculation predicts that the maximal sampling noise level at which the recovery of the top population component remains possible is higher than its counterpart in the spiked covariance model with no prior information.

I. INTRODUCTION

In the era of big data inferring features of complex systems, characterized by many degrees of freedom, is of crucial importance. The high-dimensional setting, where the number of features to extract is not small compared to the number of available data, makes this task statistically or computationally hard. One case of practical interest is the inference of the largest component (eigenvector) of correlation matrices. Consider T independently drawn observations of N interacting Gaussian variables, i.e. such that the population covariance matrix C is not the identity matrix. If T is much larger than N the empirical covariance matrix Ĉ computed from the T observations converges to C, and recovering the top eigenvector is easy. The case where both N, T are large (sent to infinity at fixed ratio r = N/T ) has received a lot of attention, both theoretically and practically [START_REF] Johnstone | High dimensional statistical inference and random matrices[END_REF]. From a theoretical point of view, it has been shown, in the case of a covariance matrix C with one (or few compared to N ) eigenvalues larger than unity, say, γ, that recovery is possible if r is smaller than the critical value r c = (γ -1) 2 [START_REF] Hoyle | Principal-component-analysis eigenvalue spectra from data with symmetry breaking structure[END_REF], [START_REF] Baik | Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices[END_REF]. For larger sampling noise (r > r c ), the top eigenvector of Ĉ is essentially orthogonal to the top component of C, and is therefore not informative. It is reasonable to expect that the situation will improve in the presence of additional, prior information about the structure of the top component to be recovered, and that recovery will be possible even when r is (not too much) larger than r c . That this is indeed the case has been rigorously shown when all entries are nonnegative [START_REF] Montanari | Non-negative Principal Component Analysis: Message Passing Algorithms and Sharp Asymptotics[END_REF], and is supported by strong numerical evidence when the top component is known to have large entries (finite as N → ∞) [START_REF] Villamaina | Estimating the principal components of correlation matrices from all their empirical eigenvectors[END_REF]. In the present work, using techniques from statistical physics we propose explicit conjectures about the critical noise level and its dependence on the signal eigenvalue (γ) and on prior knowledge. The framework is general and can be applied to any kind of entry-wise prior probability, i.e. factorized over the entries ξ i of the top component ξ. We show how rigorous results in the nonnegative case of [START_REF] Montanari | Non-negative Principal Component Analysis: Message Passing Algorithms and Sharp Asymptotics[END_REF] are recovered, and present new results for the large entry prior.

The motivation to consider the latter prior stems from computational biology, more precisely, from the study of coevolution between amino acids in protein families. Sequences of proteins diverged from a common ancestor widely differ across many organisms, while the protein structure and function are often very well conserved. The constraints induced by structural and functional conservation manifest themselves as correlations between amino acids (the N variables, where N is the protein length) across the different organisms (the T observations). Recently, it was shown that the eigenmodes ξ of the amino-acid covariance matrix corresponding to low eigenvalues were informative about three-dimensional contacts on the protein structure [START_REF] Cocco | From principal component to direct coupling analysis of coevolution in proteins: Low-eigenvalue modes are needed for structure prediction[END_REF]. These modes show large entries on the protein sites and amino-acid types in contact; as the other entries contain diffuse, non-structural signal [START_REF] Jacquin | Benchmarking inverse statistical approaches for protein structure and design with exactly solvable models[END_REF], the components ξ cannot be thought of as being sparse. The presence of large entries in structurally-informative components was empirically assessed through the so-called inverse participation ratio, i ξ 4 i (for normalized ξ), a quantity that remains finite for components with (few) large entries and otherwise vanishes for N → ∞. We hereafter use this quantity as a prior over the components to facilitate their recovery.

II. PROBABILISTIC FRAMEWORK

A. Spiked covariance model

We consider the popular Spiked Covariance Model, in which the entries of N -dimensional vectors, x = (x 1 , x 2 , ..., x N ), are Gaussian random variables with zero means and population covariance matrix C. All eigenvalues of C but one are equal to unity, while the remaining eigenvalue is γ = 1, with associated eigenvector u. As usual we choose γ > 1 but our results could be transposed to the case γ < 1 with minor modifications. We draw T independent samples x t , t = 1, 2, ..., T , and define the sample covariance matrix Ĉ, with entries Ĉij = 1 T t x t i x t j . The top eigenvector of Ĉ is denoted by ξ. In the double limit N, T → ∞ at fixed ratio r = N/T , there exists a phase transition at a critical value of the sampling noise r c = (γ-1) 2 separating the high-noise regime, r > r c , in which ξ and u have asymptotically zero squared dot product, and the lownoise regime, r < r c , where the squared dot product between ξ and u is strictly positive with high probability [START_REF] Hoyle | Principal-component-analysis eigenvalue spectra from data with symmetry breaking structure[END_REF], [START_REF] Baik | Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices[END_REF].

B. Likelihood of principal component ξ

The sample covariance matrix Ĉ obeys a Wishart distribution, determined by C, N and T . Using Bayes formula we may write the likelihood (density of probability) for the normalized top component ξ of Ĉ as follows

ρ(ξ) ∝ exp   r β 2 i,j ξ i Ĉij ξ j   δ(ξ 2 -1) , (1)
up to a normalization coefficient. Parameter β above is equal to β Bayes = 1 -1 γ . However, it is convenient to consider β as a free parameter. The β → ∞ limit corresponds to Maximum Likelihood inference, while working at low values of β may be useful to ensure rapid mixing of Monte Carlo Markov Chain sampling of distribution ρ, especially in the presence of prior information, see below.

C. Prior information over ξ

We now assume that prior information over the population eigenvector u is available under the form of a potential V acting on the entries of ξ. The posterior distribution over the top component now reads

ρ(ξ) ∝ exp   r β 2 i,j ξ i Ĉij ξ j - i V (ξ i )   δ(ξ 2 -1) , (2)
up to a normalization coefficient. Three choices for the potential V are shown in Fig. 1. Motivated by previous works on protein sequence analysis, see Introduction, we will hereafter mostly concentrate on V (ξ) = -V 0 ξ 4 , with V 0 ≥ 0 (Fig. 1(a)). This potential favors the presence of large entries in the top component, but does not rule out the existence of many entries with small magnitude (typically, of the order of N -1/2 ). It is therefore different from sparsity-enforcing potentials, such as V (ξ) ∝ |ξ| in Fig. 1(b). Exact results for the location of the transition in the nonnegative case (Fig. 1(c)) were recently derived [START_REF] Montanari | Non-negative Principal Component Analysis: Message Passing Algorithms and Sharp Asymptotics[END_REF]. Our formalism finds back those results, and can be applied to any potential V as shown below.

III. CALCULATION OF PHASE DIAGRAM

A. General replica calculation

We assume that the logarithm, divided by N , of the normalization coefficient of ρ in Eq. ( 2), the replica method, a non rigorous technique commonly used in the statistical physics of disordered systems [START_REF] Engel | Statistical Mechanics of Learning[END_REF], see [START_REF] Advani | Statistical mechanics of optimal convex inference in high dimensions[END_REF], [START_REF] Nakanishi-Ohno | Sparse approximation based on a random overcomplete basis[END_REF], [START_REF] Krzakala | Probabilistic reconstruction in compressed sensing: algorithms, phase diagrams, and threshold achieving matrices[END_REF] for recent applications to related high-dimensional inference problem. We obtain

Z( Ĉ) = ξ 2 =1 dξ exp   r β 2 i,j ξ i Ĉij ξ j - i V (ξ i )   (3) is concentrated around its expectation value L ≡ 1 N E Ĉ[log Z] in the N, T → ∞ limit,

and compute the latter with the help of

L(r, γ, β) = Extr - 1 2r log 1 β -q 0 + q 1 -q 0 q0 + q 1 2r( 1 β -q 0 + q 1 ) + q0 -q1 1 + (γ -1)p 2 + μ + q 1 q1 + pp - 1 2 log μ + q1 2μ + L(μ, p, u) , (4)
where the Extremum is computed over the order parameters

p 2 = 1 N E Ĉ[ ( i u i ξ i ) 2 Ĉ], q 0 = 1 N i,j C ij E Ĉ[ ξ i ξ j Ĉ], q 1 = 1 N i,j C ij E Ĉ[ ξ i Ĉ ξ j Ĉ],
and the conjugated Lagrange multipliers q0 , q1 , μ. Here, • Ĉ denotes the expectation over the ρ distribution over ξ in Eq. ( 2). The term L in Eq. ( 4) depends on the prior potential V and on the structure of the normalized population top component u, more precisely, on how its entries scale with N . We now consider two cases of interest.

B. Case of nonnegative entries

We assume first that the components of u scales as ũi √ N , with ũi finite, and denote by ϕ(ũ) the distribution of the ũi . We focus on the nonnegative entry prior, for which V (ξ) = +∞ for ξ < 0 and 0 for ξ ≥ 0. We obtain

L = p2 2μ + 1 2μ dũ ϕ(ũ) log erfc -p ũ -z √ 2q 1 2 √ μ z , ( 5
)
where erfc is the complementary error function, and [F (z)] z denotes the average of F (z) over the Gaussian measure, e -z 2 /2 / √ 2π. After some elementary algebra, we obtain the expression for the overlap

p = [ũ(x ũ + z)] + z [(x ũ + z) 2 ] + z with x = p √ 2q 1 = (γ -1) p r(1 + (γ -1)p 2 ) (6) where [F (z))] + z = [max(F (z), 0)] z .
These equations correspond to Eqs. ( 7), ( 8), ( 9), ( 21) and (23) of [START_REF] Montanari | Non-negative Principal Component Analysis: Message Passing Algorithms and Sharp Asymptotics[END_REF].

C. Case of large entries

We now assume that u has only K 'large' entries, u 1 , ..., u K , different from zero in the N → ∞ limit (with finite K), and that the other entries decay fast enough with N , e.g. are of the order of 1

√ N . Then L = - K i=1 min -1≤ξ≤1 μ ξ 2 + p u i ξ + V (ξ) , (7)
While the above formula is valid for generic V we restrict ourselves, in the remaining part of this article, to the potential V (ξ) = -V 0 ξ 4 . Furthermore, we assume that the K finite entries of u are all equal to u = 1 √ K ; The calculation can be easily extended to u < 1 √ K , or to the case of nonhomogeneous entries u i . In addition we assume that (A1) all ξ i s take identical values in Eq. ( 4) (homogeneous regime); (A2) ξ has no large (finite when N → ∞) entry ξ i on sites such that u i = 0. The validity of these assumptions will be discussed in the next Section.

After some elementary algebra the extremization conditions reduce to the following set of K + 1 coupled equations over μ and the first K entries of ξ:

ξ i - 2 V 0 μ ξ 3 i = (γ -1)p r √ K 2μ β -1 (i = 1, ..., K) , K i=1 ξ 2 i = 1 - 1 2μ - 1 + (γ -1)p 2 -1 2μ r 2μ β -1 2 (8)
where

p = 1 √ K K i=1 ξ i . (9)
Note that the K variables ξ i s obey the same cubic equation and, hence, can take at most three different values as i varies.

The equations above admit the solution p = ξ i = 0, μ = 1 2 , corresponding to the 'unaligned' phase. In addition, in some well-defined regions of the four-dimensional control parameter space (r, γ, β, V 0 ) solutions with p = 0 may be found. We give in Section IV below results for three cases: (A) no prior (V 0 = 0); (B) weak exploitation of many data with prior information (small β, r for finite V 0 ); (C) maximum a posteriori inference (finite r, and infinite β and V 0 at fixed ratio V 0 /β).

IV. RESULTS FOR 'LARGE ENTRY' PRIOR

A. Warm-up: no prior

We start with the case V 0 = 0. Extremization conditions over the parameters in Eq. ( 4) give the value of the squared overlap p 2 between u and ξ for any β. We find: p 2 = 0 for r > r c = (1 -γ) 2 whatever the value of β, and, for r < r c ,

p 2 = 1 - β(r) β 1 -r r 2 c 1 + r rc if β > β(r) ≡ r r + γ -1 .
(10) Those expressions perfectly agree, in the β → ∞ limit, with known results for the spiked covariance model [START_REF] Hoyle | Principal-component-analysis eigenvalue spectra from data with symmetry breaking structure[END_REF], [START_REF] Baik | Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices[END_REF]. In addition our formalism gives access to the value of p 2 for finite β. Note that, for r < r c , inference of the direction of u is possible, i.e. p 2 > 0, even for β ≤ β Bayes (but larger than β(r)). At the critical noise, β(r c ) = β Bayes .

B. Inference at low β with prior information

The above results imply that, in the absence of prior information (V 0 = 0), inference of the top component direction is possible at low β → 0 provided the sampling noise r is smaller than β/(1 -γ). In other words, when both β and r tend to zero with a fixed ratio β = β/r, the aligned and notaligned phases correspond, respectively, to β > βc (0) = 1 γ-1 , and β < βc (0).

We expect the critical ratio βc (V 0 ) to be a drecreasing function of the prior strength V 0 , as stronger prior should facilitate the recovery of the large-entry top component u. Resolution of Eq. ( 8) gives the phase diagram shown in Fig. 2. Several regimes can be identified, depending on V 0 :

• For 0 < V 0 ≤ K
4 , the critical ratio βc (V 0 ) remains unchanged, see Fig. 2, and equal to 1 γ-1 . As β crosses this critical value the overlap p continuously increases from 0 to a positive value, see Fig. 3.

• For K 4 < V 0 ≤ K, the aligned phase (p = 0) exist for β > βc,1 = 4 √ V0/K 1- √ V0/K γ-1
, see dashed line in Fig. 2. As β crosses βc,1 the squared overlap p 2 discontinuously jumps from 0 to 1-

1 2 √ V0 > 0, see Fig. 3. • For K 4 < V 0 ≤ V + 0 × K, the aligned phase (p = 0
) is thermodynamically stable, meaning that the value L + of L in Eq. ( 4) attached to this phase is larger than the one of the nonligned phase (p = 0), L 0 = 1 2 (1 + β), for β > βc,2 , see full line in Fig. 2. The value of βc,2 and of the overlap p 2 at the phase transition are the roots of the two coupled equations

0 = p 2 2 1 -p 2 2 + 1 2 log 1 -p 2 2 - V 0 K p 4 2 , βc,2 = 1 1 -p 2 2 -4 V 0 K p 2 2 , (11)
where the first equation actually implements the condition L + = L 0 . The phase transition is illustrated in the middle panel of Fig. 3 for a specific value of V 0 . The value of V + 0 1.227703... is defined from

V + 0 = min V 0 s.t. max ξ =0 log(1 -ξ 2 ) + 2V 0 ξ 4 > 0 (12) • For V 0 ≥ V + 0 ×
K the prior strength is so strong that the inferred component ξ has few large entries whatever the value of β. For β > 0 it is aligned (p = 0) with u (Fig. 3), while for β < 0, it is not, see Fig. 2. Assumption (A1), see Section III-C, is trivially valid for K = 1, but is not necessarily correct for K ≥ 2 and strong prior strength, for which we expect that ξ will condensate and one component, say, ξ 1 , will be larger than the other components, say, ξ i , with i = 2, ..., K (nonhomogeneous regime). The transition line between these two regimes is (1, 1, ..., 1, 0, 0, ..., 0). Dots locate the points (1, 1 4 ), ( 89 , 4 9 ) and (0, 1.228...), see text. Transitions between phases are shown by full lines. The dashed lines show the limit of existence of the Aligned phase, while the dot-dashed lines separate the regions with (above) and without (below line) homogeneity breaking among the K large entries of ξ. identified upon imposing that the cubic equation over ξ i in Eq. ( 4) admits a two-fold degenerate solution ξ, that is, 6V 0 ξ 2 = μ. The transition line is plotted in the phase diagram of Fig. 2 (dot-dashed line), and ends up in the point of coordinate ( 4 9 , 8 9 ) lying on the existence line (dashed line).

C. Inference at high β with strong prior information

We now focus on MAP inference at finite sampling noise r, whereas β and V 0 are both sent to infinity, with a fixed ratio S = r V 0 /(β). Parameter S, hereafter referred to as the slope, controls the relative magnitude of the Ĉ-dependent and prior terms in ρ, see Eq. ( 2), while the multiplicative factor r is introduced in the definition of S to compensate for the explicit dependence upon r in ρ. For simplicity we present results for K = 1 only, the extension to larger K being rather straightforward. Equations ( 8) admit the solution p = 0, μ = 1 2 , and another solution, with μ → ∞ as β, V 0 → ∞, with ratios β/μ, V 0 /μ having finite limits. After some simple algebra we obtain the following expresson for the slope as a function of the squared overlap for the latter solution:

S(p 2 ) = r -(γ -1)y 1 + y 4 p 2 y with y = r(1 -p 2 ) 1 + (γ -1)p 2 .
(13) We show in Fig. 4 the representative curve of p 2 vs. the slope S, for r below and above the critical noise level, r c = (γ -1) 2 , in the absence of prior. For r < r c the squared overlap is an increasing function of S, starting from a non zero value for S = 0: the population eigenvector direction can be estimated without prior at low sampling noise [START_REF] Baik | Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices[END_REF], but the overlap is increased in the presence of prior. For r > r c a discontinuous jump is observed from p = 0 to p > 0 at a critical value of the slope, S -= min

p 2 >0 S(p 2 ) (r > r c ) , (14)
while the overlap further increases as S exceeds S -. Remarkably, even for large sampling noise values, the presence of a sufficiently strong prior allows us to infer u. The value of S - as a function of the noise level r is shown in Fig. 5; for large noise levels we have S -= r 4 + 3 4 4/3 γ 1/3 r 2/3 + O(r 1/3 ). This aligned phase competes with a nonaligned, but condensed phase, in which assumption (A2), see Section III-C, is violated. In other words, for S and r sufficiently large, ξ has few large entries (ξ 2 j > 0 in the N → ∞ limit), but not along the directions i corresponding to the large components of u; hence, p 2 = 0. To describe this new phase we set u i to 0 in the expression for L in Eq. [START_REF] Jacquin | Benchmarking inverse statistical approaches for protein structure and design with exactly solvable models[END_REF]. The corresponding optimization equations can be solved in the β, V 0 → ∞ limit, with the result that the nonaligned, condensed regime exists for S larger than The value of S + as a function of the noise level r is shown in Fig. 5. For small r (slightly above r c ) we observe that S + is larger than S -, as intuitively expected: it is favorable to condense ξ along the direction of u rather than any other direction. It can be checked that the value of L in Eq. ( 4) is higher for this phase than for the aligned condensed phase. Hence, as soon as S exceeds S + the overlap p vanishes.

S + = min 0<b< √ r 1+ √ r r 2 (1 -b) 2 4b r(1 -b) 2 -b 2 . ( 15
) 0 0.5 1 1.5 2
For large noise levels, however, we have S + = r 4 + 3 4 4/3 r 2/3 + O(r 1/3 ), which is asymptotically smaller than S -, see insert of Fig. 5. Indeed, the threshold slopes S -and S + cross at a well-defined value of the noise, r d , which depends on the top eigenvalue γ. We show in Fig. 6 the behaviour of r d vs. γ, and compare it to the critical noise r c in the absence of prior. We observe the presence of a region in the (γ, r) plane, above the critical line r c , where the direction of u can be inferred thanks to the large-entry prior. Our replica symmetric theory predicts that the benefit of the prior does not extend to very large values of the signal eigenvalue γ and of the noise r, see Fig. 6.

V. CONCLUSION

The non rigorous calculations done in this paper suggest that inference of the top component of the population covariance matrix is possible in the presence of prior information, even above the critical noise level of the spiked covariance model, in agreement with rigorous results for the nonnegative case [START_REF] Montanari | Non-negative Principal Component Analysis: Message Passing Algorithms and Sharp Asymptotics[END_REF] and numerical investigations for the large-entry case [START_REF] Villamaina | Estimating the principal components of correlation matrices from all their empirical eigenvectors[END_REF]. Many interesting questions have not been investigated in the present paper: how hard is the recovery problem from a computational point of view? Are there 'dynamical' phase transitions separating subregions in the aligned phase, where the top component can be recovered in polynomial time or not? If so how do these line compare to the 'static' critical lines derived in this paper? In addition it would be interesting to investigate the validity of the replica-symmetric hypothesis used to derive the results above [START_REF] Engel | Statistical Mechanics of Learning[END_REF]. Though replica symmetry is generally expected to be correct for convex optimization problems what happens in nonconvex situations is not clear. For instance, inference of the top component with the nonnegative prior gives rise to a nonconvex optimization problem [START_REF] Montanari | Non-negative Principal Component Analysis: Message Passing Algorithms and Sharp Asymptotics[END_REF], but all rigorous results are exactly found back within our replica symmetric approach, see Section III.B. From this perspective it would be useful to investigate whether the present results are robust against replica symmetry breaking.

Fig. 1 .

 1 Fig. 1. Three potentials V (ξ) corresponding to three prior information about the entries ξ of the principal component: (a) large entries are present, (b) L 1 penalty favoring zero entries, (c) all entries are nonnegative.

2 .

 2 Fig. 2. Phase diagram of top component recovery in the ( β, V 0 ) plane; axis are rescaled by (γ -1) and 1 K factors, where K is the number of nonzero components in u = 1 √ K

2 Fig. 3 .

 23 Fig.3. Average squared overlap between top components of the population and sample covariance matrices, p 2 ≡ E Ĉ[(u • ξ) 2 ], vs. control parameter (γ -1) β for three prior strengths and K = 1.

5 Fig. 4 .

 54 Fig. 4. Average squared overlap between top components of the population and sample covariance matrices, p 2 , vs. slope S for sampling noises r = 0.2 (top) and 0.4 (bottom). Note the presence of the discontinuous transition at S - 2.19 in the latter case. The randomly condensed solution appears for S > S + 3.68. Here, rc = (γ -1) 2 = 0.25.

Fig. 5 .

 5 Fig. 5. Behaviours of S -(full line) and S + (dashed line) as functions of the sampling noise r. Insert: S + -S -vs. r. The difference vanishes in r d 6.54. Here, rc = 0.25, i.e. γ = 1.5.

Fig. 6 .

 6 Fig.6. Replica-symmetric phase diagram of top component inference with large-entry prior in the (γ, r) plane. The dashed line rc = (γ -1) 2 divides the plane into the weak noise region (below line), where inference is possible withour prior, and the strong noise region (above line). The full line shows the value of r d , at which S -= S + , as a function of γ. In between the dashed and full lines, inference of the top component is possible in the presence of a prior with appropriate strength. The two lines merge for γ 5.3, r 18.5.Here, K = 1.

ACKNOWLEDGMENT I am grateful to M.R. McKay for the invitation to the Asilomar 2016 conference, and to S. Cocco and D. Villamaina for useful discussions. This work has benefited from the financial support of the CNRS Inphyniti Inferneuro project.


C Evripidou

A Burns
email: alan.burns@york.ac.uk

Scheduling for Mixed-criticality Hypervisor Systems in the Automotive Domain

Virtualisation has been proposed for use in the automotive domain as it has the potential to reduce the number of ECUs (Electronic Control Units) that are required in a modern vehicle. In this paper we first introduce a visualisation architecture that makes use of two different types of execution-time servers to provide separation, low run-time overheads but short response-times for event-triggered computation. This model is then extended to mixed-criticality systems and utilises a run-time switch between the task-server mapping to enhance schedulability (at the cost of extended response-times). An industrial case study is used to evaluate the approach.

I. INTRODUCTION

The automotive industry has been using software in cars for over 30 years, but this is now increasing at a very fast pace [START_REF] Broy | Challenges in automotive software engineering[END_REF], [START_REF] Broy | Engineering automotive software[END_REF]. The software in a vehicle typically runs within Electronic Control Units (ECUs), which are embedded hardware platforms responsible for controlling different subsystems within a vehicle. Examples of ECUs are door control units, transmission control units and engine control units.

As stated by Broy [START_REF] Broy | Challenges in automotive software engineering[END_REF] in 2006, about 40% of the production cost of a vehicle was spent on electronics and software. After the 30-year growth of the volume of software in vehicles, a modern premium car may contain in excess of 100 processors spread across 70 ECUs [START_REF] Abelein | Complexity, quality and robustness-the challenges of tomorrow's automotive electronics[END_REF], [START_REF] André | Operating systems technology for converged ECUs[END_REF], [START_REF] Broy | Challenges in automotive software engineering[END_REF], [START_REF] Broy | Engineering automotive software[END_REF]. Nolte [START_REF] Nolte | Hierarchical scheduling of complex embedded real-time systems[END_REF] studied this trend; there is a clear indication of the increased complexity and the added hardware costs that are prominent in modern vehicles.

When they were first introduced, ECUs were functionally independent and were connected solely to sensors and actuators [START_REF] Broy | Engineering automotive software[END_REF]. As ECUs were required to provide additional functionality, there was a need to establish communication channels between them. This change has resulted in multiple ECUs cooperating to provide a certain piece of functionality. In addition to numbers, ECUs can be heterogeneous and be responsible for different types of tasks; hard real-time (vehicle control) and soft real-time (infotainment). Given the large number of ECUs per vehicle, the dependencies between them and their lack of homogeneity, it can be inferred that they form a complex system that is hard to reason about.

Apart from the structural complexity of the electronic parts of modern vehicles, another challenge in software engineering for the automotive domain is the difference in lifetime of car models and ECU hardware [START_REF] Broy | Engineering automotive software[END_REF], [START_REF] Pretschner | Software engineering for automotive systems: A roadmap[END_REF]. Specifically, a car model typically has a production lifetime of about 7 years, whereas a microprocessor 5 years. In addition to the 7 years of production lifetime, a car manufacturer (OEM) needs to provide service and spare parts for an additional 15 years. The result of this difference in lifetimes is that during a car's lifetime it is very likely that some ECU hardware components stop being available in the market. Given that the ECU software is usually highly optimised for the underlying hardware, porting to a newer platform can be a hard and expensive task.

In order to deal with the increased complexity that characterises modern vehicles, OEMs, suppliers and other relevant companies formed a worldwide development partnership. The result of this partnership is the Automotive Open System Architecture (AUTOSAR) [START_REF]AUTomotive Open System ARchitecture[END_REF]. AUTOSAR aims to provide a common architecture as well as a methodology that will help with the understanding of the interaction of ECUs, allow software reuse and enable the combination of multiple functions on a single ECU [START_REF] Kirschke-Billern | AUTOSAR -A worldwide standard: Current developments, roll-out and outlook[END_REF].

In this paper we consider the use of virtualisation (see Section II) for the integration of multiple ECUs into a single hardware platform. Section III outlines the general scheduling approach followed for virtualisation. The key characteristic of the scheduling approach is the use of a deferrable and a periodic server per partition for the execution of eventdriven and time-driven tasks; event-driven tasks execute on deferrable servers for low response times and time-driven tasks execute on periodic servers for the improved utilisation and lower overheads. The main contribution of this paper is the extension of the proposed scheduling approach for virtualisation to support two modes of degradation (D1 and D2). The mixed criticality model is presented in Section IV All the work presented was done in close cooperation with the automotive company ETAS. They provided the requirements and application code that was used to form the presented model, and derive the case study of Section V, which was used to evaluate the proposed approach.

II. VIRTUALISATION

Virtualisation is a technique, initially developed in the early 60's [START_REF] Buzen | The evolution of virtual machine architecture[END_REF], where logical resources are created in order to allow one or more applications to execute on the same hardware platform. The logical resources are created and managed by the hypervisor (HV), also referred to as virtual machine manager (VMM). From our experience with ETAS, there is increasing interest in the automotive industry for the use of virtualisation to alleviate some of their current problems.

The main use case for HV technology in the automotive domain is the reduction of ECU count by combining multiple ECUs on a single hardware platform. The key properties that must hold in a HV system is spatial and temporal isolation of the VMs. Spatial isolation is achieved by prohibiting the VMs from accessing memory areas outside of their memory space. Temporal isolation, which is the focus of the research reported in this paper, is the property under which a VM's behaviour cannot cause another VM to violate its real-time properties.

Early work on virtualisation by Popek and Goldberg [START_REF] Popek | Formal requirements for virtualizable third generation architectures[END_REF] in 1974 identify three characteristics of an HV. First, an HV needs to be able to provide its hosted virtual machines (VMs) with an execution environment which is indistinguishable from real hardware. Second, execution is efficient by mapping a large subset of the virtual processor instruction set to a physical processor. Third, the HV has complete control over all hardware resources and is able to allow or prohibit VMs access, according to their configuration.

Popek and Goldberg [START_REF] Popek | Formal requirements for virtualizable third generation architectures[END_REF] also identify a set of properties for HVs:

• Efficiency: all non-privileged instructions are executed directly on hardware. • Resource control: it is impossible for a VM to interfere with any system resources that are not allocated to it. • Equivalence: a VM produces the same results when executing as if it was executing without a HV. The HV characteristics by Popek and Goldberg [START_REF] Popek | Formal requirements for virtualizable third generation architectures[END_REF] refer to full virtualisation. Full virtualisation implies that the applications in VMs can be executed without requiring any modifications [START_REF] Kaiser | Combining partitioning and virtualization for safetycritical systems[END_REF]. In order to maintain the properties identified above in a fully virtualised environment it is necessary to have adequate hardware support. Specifically, allowing a VM to execute directly most of the time on the underlying hardware for efficiency requires that the HV will be able to identify attempts to execute privileged instructions. The HV is then responsible for checking whether the VM is allowed to perform the operation it attempted to and act accordingly.

Paravirtualisation was introduced as means of alleviating the lack of hardware support and to simplify the development of the HV. Specifically, in a paravirtualised environment, VMs execute directly on hardware using modified versions of their application code [START_REF] Kaiser | Combining partitioning and virtualization for safetycritical systems[END_REF], using HV calls to replace the functionality of privileged instructions. Examples of virtualisation platforms are: OKL4 [START_REF]OKL4 microvisor[END_REF], XtratuM [START_REF] Xtratum | [END_REF], Xen [START_REF] Barham | Xen and the art of virtualization[END_REF] and PikeOS [START_REF] Sysgo | PikeOS RTOS and virtualization concept[END_REF].

III. REQUIREMENTS AND APPROACH TO VIRTUALISATION

The driving requirements for the architecture described in this paper is: the need to partition the code of the applications, the need to have low run-time overheads (the automotive industry is very dependent of the efficient use of hardware resources), and the need to provide very fast response-times for a certain class of application code (namely, event-triggered computation). The latter two requirements are contradictory, fast response times require small containers, efficient implementation dictates the use of larger containers. This dilemma is resolved by the use of two different types of container.

A. Task Model

Traditional hypervisor scheduling approaches were developed assuming no visibility at the task level in the partitions. In this section we define a flexible task model, which allows exploiting visibility, where that is available, taking into account implementation overheads. First, we classify tasks as synchronous or asynchronous. Synchronous tasks are strictly periodic, whereas asynchronous are event-triggered tasks with a known minimum inter-arrival time. The main motivation behind the proposed scheduling method is that asynchronous tasks require quick response times, while synchronous ones can be serviced in a more efficient, lower overheads approach. All operations performed by the hypervisor are executed in a non-preemptive manner and are described as highly predictable pieces of code. The execution of all tasks is defined to be structured as shown in Figure 1. A task τ i is defined by the tuple

(C cs1 i , C i , C cs2 i , T i , P i ): • C cs1
i : the scheduling and context switching overheads required before the execution of the main task body. In asynchronous tasks, this section is performed by the hypervisor and is therefore non-preemptible. Synchronous tasks can be preempted by the hypervisor during this section but not by other synchronous tasks of the same partition, since this operation is performed by the partition.

• C i : the time required for the task's main body to execute.

Partition tasks (synchronous and asynchronous) can be preempted while in this section, however hypervisor tasks run non-preemptively. • C cs2 i : the overheads of terminating the execution of the task. The preemption rules that apply for C cs1 i also apply during the execution of this section.

• T i : the period or minimum inter-arrival time of the task.

• P i : the priority level of the task.

B. Execution Servers

Our approach aims to minimise the response time for eventtriggered tasks, while at the same time maximise schedulability and enforce temporal protection between the different partitions 1 . The CPU time is shared between the partitions using execution servers. Temporal protection is achieved by prohibiting partitions to execute for more than their servers' capacity. Each server is associated with a hypervisor task, which is responsible for replenishing its capacity.

Asynchronous tasks are released in response to events and therefore need to be serviced with the lowest possible latency. To facilitate this requirement, asynchronous tasks execute using deferrable servers that are assumed to always have enough capacity to service all event-driven tasks, given their WCET, hypervisor overheads and period. With the use of a deferrable server, no server capacity is expended when the system is idle, and events are serviced as they arrive, provided they have the highest priority in the system.

Deferrable servers (DS) offer low response time for the asynchronous tasks, however they are inferior in terms of schedulability in comparison to periodic servers. The timedriven or synchronous tasks in the system are execute using a periodic server (PS) in order to alleviate this trade off, therefore improving schedulability, without compromising on the low latency required by event-driven tasks.

The association between servers and tasks is defined using the association matrix M . The rows of the matrix represent the tasks in the system, whereas the columns are the servers. All elements can take the values 0 or 1. If a task τ i is serviced by s j , then M τi,sj = 1, otherwise M τi,sj = 0. Moreover, a task can be serviced by exactly one server, which implies that the sum of each row results is at most 1.

M =             HV DS 0 DS 1 P S 0 P S 1 τ 0 1 0 0 0 0 τ 1 1 0 0 0 0 τ 2 1 0 0 0 0 τ 3 1 0 0 0 0 τ 4 0 0 1 0 0 τ 5 0 1 0 0 0 τ 6 0 0 0 1 0 τ 7 0 0 0 0 1             (1)
Equation ( 1) is an example configuration of an association matrix of a simple system with two partitions (applications), p 0 and p 1 . Each of the partitions is associated with a deferrable and a periodic server. Specifically, DS 0 and P S 0 are associated with partition p 0 and DS 1 and P S 1 are associated with partition p 1 . Each of the servers requires a hypervisor task so that its capacity is replenished periodically. The hypervisor tasks that are responsible for replenishing the server's capacity are τ 0 , τ 1 , τ 2 and τ 3 . Partition p 0 has two tasks, τ 5 and τ 6 , and p 0 also has two tasks, τ 4 and τ 7 . τ 4 and τ 5 are asynchronous tasks and are therefore associated with DS 0 and DS 1 respectively. Similarly, τ 6 and τ 7 are synchronous tasks and are associated with P S 0 and P S 1 respectively.

C. Priority Space

Figure 2 shows the relationship between the execution modes in terms of their corresponding priority levels. The hypervisor executes in hypervisor mode at the highest system priority region. Since the hypervisor's code is trusted and is typically consisted of short, highly predictable non-preemptive tasks. The motivation behind this approach is to allow eventtriggered tasks, such as the service routines of exceptions or interrupts, to execute at a priority level that is strictly higher than any time-triggered tasks.

Synchronous tasks have the lowest priority range in the system. Specifically, a partition executes in synchronous mode if no event is pending and eligible to be handled. The eligibility of handling an event is directly associated with the asynchronous budget. A partition executing in synchronous mode is responsible for doing its own task scheduling and may therefore maintain its own internal priority space.

The analysis employed with the defined model uses standard response-time analysis to ensure:

• each server (DS and PS) is schedulable, and • each task within each server is schedulable. For space reasons the analysis is not given in this paper (see [START_REF] Evripidou | Scheduling for Mixed-criticality Hypervisor Systems in the Automotive Domain[END_REF]).

IV. APPROACH TO MIXED CRITICALITY

Having now introduced the approach to virtualisation we can describe how this model is extended to support mixed criticality. Separation is already supported. We now wish to provide a means by which criticality mode changes can be accommodated.

The introduction of safety in the automotive industry with the ISO26262 standard on functional safety of road vehicles [START_REF]road vehicles -functional safety[END_REF] gave rise to concerns regarding the integration of components with different ASIL (Automotive Safety Integrity Level) in the same ECU. AUTOSAR enables the integration of SWC (Software Components) from different vendors, requiring modifications to the configuration of the RTE (Runtime Environment) and BSW (Basic Software) [START_REF] Morgan | Safety and security with hypervisor technology[END_REF]. From a safetycritical perspective AUTOSAR lacks the required separation mechanisms. Specifically a failure in an AUTOSAR SWC typically results in an ECU reset. In the case where all SWCs are of the same ASIL this is acceptable, however in the case of mixed-criticality this could potentially allow a low criticality component to interfere with a higher criticality one. This directly violates freedom from interference (FFI), as dictated by ISO 26262-6:2011 Annex D [START_REF]road vehicles -functional safety[END_REF]. Moreover, Esper et al. [START_REF] Esper | How realistic is the mixed-criticality real-time system model[END_REF] state that suspending the execution of lower criticality tasks may hinder the certification progress due to lack of isolation unless the suspended tasks are non-critical.

A possible way of achieving FFI is to ensure that components of different criticality levels are located on separate physical ECUs. Failures in lower criticality components would therefore be isolated and not propagate to other ECUs of higher criticality. This approach could potentially increase the number of physical ECUs, which, as identified earlier, is one of the main drivers to the high development costs for software and hardware in vehicles. The use of a hypervisor can provide the necessary isolation between its partitions to allow the integration of multiple ECU images on a single physical ECU, while enforcing FFI.

The proposed mixed-criticality model supports three levels of criticality, LO, M I, HI, with HI > M I > LO. The task model is defined in a similar manner as the model proposed by Vestal [START_REF] Vestal | Preemptive scheduling of multi-criticality systems with varying degrees of execution time assurance[END_REF]. A task τ i is defined by the tuple ( C cs1 i , C i , C cs2 i , T i , P i , L i ), where: • C cs1 i : A vector containing the WCET of the implementation overheads before the execution of the task's main body for each criticality level, C cs1 i (LO), C cs1 i (M I) and C cs1 i (HI). • C i : A vector containing the WCET of the task's main body at each criticality level.

C i (LO) ≤ C i (M I) ≤ C i (HI).
• C cs2 i : The WCET of the overheads after the execution of the task's main body for each criticality level, C cs2 i (LO), C cs2 i (M I) and C cs2 i (HI). • T i : The period of the task at each level of criticality, T (LO) ≥ T (M I) ≥ T (HI). • P i : The priority level of the task at each criticality level, P i (LO), P i (M I), P i (HI). • L i : The criticality level of the task (ie. LO, M I, HI).

A. Execution Modes

Figure 3 summarises the mode transitions performed by the system to support multiple levels of criticality. A requirement of implementing the mixed-criticality model is to monitor the execution time of all partition tasks. As the figure illustrates as a task executes for more than its C(LO) value the first degraded mode is entered; if necessary a second level of degradation is entered. A distinctive aspect of the model is that D1 does not involve the dropping of any tasks. Rather, it invokes an architectural change in which the event-triggered tasks executing in deferrable servers are migrated to their corresponding partitions' periodic servers. Hence all tasks continue to execute but the event-triggered work has increased latency; the periodic work continues to meet all deadlines. If a further degradation is required then D2 does involve the dropping of low criticality work (following the model of Vestal [START_REF] Vestal | Preemptive scheduling of multi-criticality systems with varying degrees of execution time assurance[END_REF]).

V. EVALUATION

The proposed approach was heavily motivated by the requirements of the automotive industry and takes into account implementation overheads. In this section a case study, derived from application code provided by ETAS Ltd., is used to evaluate the proposed approach.

The case study evaluation was performed using ECU application code that was provided by ETAS Ltd. The application code consisted of a set of AUTOSAR TASKs of a Mercedes-Benz M160 engine controller. The functionality of the application code includes controlling the air flow, idle speed, and fuel injection. Each TASK has a unique period and a set of sub-tasks that are to be executed at that period. The provided taskset was split into two partitions with different criticality levels by an expert with respect to the individual task functionality. Each task was then classified as asynchronous or synchronous, with respect to their real-time requirements.

A. Task Measurement 1) Application Tasks: The next step for the case study is to analyse the provided ECU code, in order to realise the proposed model. The source contained only application code, without the underlying OS. The minimum information that is required for the proposed model is the period of each task, which was provided and the execution time. The execution times of the tasks were obtained using a measurement-based approach. The first part of the timing analysis was to study the provided code. Upon inspection, the code was primarily linear with minimal branching. Numeric calculations and variable conversions were the primary operations performed by the code. The input/output of the tasks was made by reading/writing to external variable, each with a clearly defined range of valid values.

Measuring the task execution times required modifications to the provided code, since the OS code was not available. Specifically the code was modified with definitions of all the external variables and structures, that were initialised using valid values. The granularity of the timer was insufficient for the measurement of the short task execution times, therefore each task execution time was estimated by measuring the amount of time required to execute it one million times in sequence. Each such measurement was then taken 500 times, therefore obtaining a sample of 500 estimated execution times per task.

Table I summarises the sampled execution times in ns. From the table, there is a small difference between max and mean as well as a small standard deviation of the execution times. This confirms that there is little variation between each measured execution time, as expected from the minimal branching that was observed during the inspection of the code. Analysis using Pearson's r showed no correlation between the WCET and the period of tasks, r(30) = -0.0622, p = 0.735 (two-tailed test). To calculate the hypervisor overheads, a partial hypervisor implementation was built. The hypervisor overheads in the system are responsible for the replenishment of the server capacity and handling forwarding and returning from interrupts. The WCETs of the developed hypervisor components were calculated using the maximum cycles required for the ARM1176JZF-S instructions to execute [START_REF]ARM1176JZF-S TM technical reference manual[END_REF]. Table II summarises the hypervisor overheads. These values are deemed applicable to all criticality modes. Therefore, the WCET for all operations performed by the hypervisor are not scaled up for the M I and HI criticality levels.

B. Hypervisor Overheads

C. Experiments

The case study was used to undertake sensitivity analysis using the response-time analysis developed for the two server type model. The case study taskset specification in conjunction with the system overhead parameters are used to produce a system configuration that is marginally schedulable at all criticality levels.

Sensitivity analysis is then repeated for the two degraded modes, D1 and D2, to see what improvement in schedulability results. The processor used for these experiments was considerable more powerful that the one from which the code WCET times were derived. High scaling factors were therefore observed. A scaling factor of 100 implies that all code could be increased 100 times and the system would still (just) be schedulable.

Three methods of supporting D1 were investigated:

• 2-partition -here each of the two applications uses just one DS and one PS; in D1 all tasks in the DS are moved to the PS. • 3-partition -here each application uses just one DS and one PS; in D1 the DS is converted to a new, distinct, PS (i.e. another partition is added). • 8-partition -here each application uses a number of DSs and PSs; in D1 all DSs are converted to PSs.

D. Results

In this section we discuss the results of the sensitivity analysis on the three identified system configurations.

Figure 4 shows the scaling that was achieved by each system configuration for all criticality levels. The 2-partition configuration while the system executes in N mode achieved a scaling factor of 103. At the degraded modes there is significant capacity loss, which makes the proposed model ineffective for the 2-partition configuration. The reason for the poor performance in D1 and D2 is the large variation of the temporal requirements of the application tasks served by individual periodic servers. This variation in temporal requirements results in wasted server capacity.

The 3-partition configuration achieves the same scaling as the 2-partition configuration while the system executes in N mode. After the mode switch to D1 there is a significant increase of the WCET scaling, by a factor of 3. Assigning the asynchronous tasks in a separate partition eliminated the server parameter shortcoming of the 2-partition configuration. Specifically, in the degraded modes the asynchronous tasks execute using the original periodic server. All asynchronous tasks in the system share the same minimum inter-arrival time of 1ms. This allows for no wasted server capacity, therefore increasing the achieved scaling in the degraded modes. Switching to the second degraded mode, D2, there is a small increase in the scaling factor. The small increase is as expected, since the 95% of application task utilisation is used by HI criticality tasks.

The scaling achieved with the 8-partition configuration was 80, which is significantly lower than the 2 and 3-partition configurations. The lower scaling achieved for the 8-partition configuration was attributed to the large number of hypervisor tasks, which reside at the highest priority band in the system. Specifically, server capacity replenishment tasks are assigned strictly higher priorities than application tasks, therefore having a greater impact on the optimality of the priority assignment algorithm. Another contributing factor to this is the context switching overheads of asynchronous tasks. After the switch to the first degraded mode, D1, the achieved scaling is 310, which is the highest achieved in all three configurations. The use of solely periodic servers servicing tasks of the same period made the priority assignment algorithm very effective. Switching to D2 results in more spare capacity than the 3partition configuration, since additional capacity is freed due to the number of server replenishment tasks that are no longer in use.

VI. CONCLUSIONS

In this paper we have shown how two forms of executiontime servers (namely deferrable server and periodic server) can be used (with a hypervisor) to support multiple applications (partitions) running on the same processor/core. The approach derived has been in response to the needs of the automotive domain. The servers provide the necessary level of separation, and the different characteristics of the two servers allow low latencies for event-triggered tasks and low overheads for time-triggered (periodic) tasks to be manifest. The proposed approach was extended to support three levels of criticality, featuring two degraded modes (D1 and D2). The mixedcriticality model aims to provide the partitions with additional CPU capacity by providing an alternative task-server mapping (D1) before the suspension of lower criticality tasks (D2). An industrial case study was used to evaluate the approach. This demonstrated the need to keep the overheads of the hypervisor as low as possible.

Future work will move the scheme to multi-core platforms, this is a strong requirement of the automotive domain, although it presents a number of significant challenges.

Fig. 1 .

 1 Fig. 1. Task structure.

Fig. 2 .

 2 Fig. 2. Example of a k-partition system priority space.

Fig. 3 .

 3 Fig. 3. State transitions for the mixed-criticality model.

Fig. 4 .

 4 Fig. 4. Application task WCET scaling with 2, 3 and 8-partition configurations.

A similar approach by Missimer et al.[START_REF] Missimer | Mixed-criticality scheduling with I/O[END_REF] using sporadic and priority inheritance bandwidth preserving servers was published during the development of the approach discussed in this paper. Although they also use servers, they do not have the same approach to mixed-criticality.