Harold Van Heerde
email: h.j.w.vanheerde@ewi.utwente.nl

Maarten Fokkinga
email: fokkinga@ewi.utwente.nl

Nicolas Anciaux
email: nicolas.anciaux@inria.fr

A framework to balance privacy and data usability using data degradation

Personal data is a valuable asset for service providers. To collect such data, free services are offered to users, for whom the risk of loosing privacy by subscribing to a service is often not clear. Although the services are free in terms of money, the user does not know how much he or she actually pays for a given service when allowing his or her data to be collected, unaware of taking a significant privacy risk by doing so. In practice, this risk is even not taken into account when deciding how long the data will be retained; the service provider simply wants to optimize the total worth of the stored data by retaining the data as long as possible. In this paper, we express the privacy risk for the user in terms of such a retention period; the user wants to optimize its privacy by allowing the data to be retained as short as possible. Now, in stead of only considering the interests of the service provider, we argue that we should optimize the common interest of both parties, and present a framework to reason about worth and privacy to find such optimum. Going one step further, we refine and generalize limited retention to data degradation, which prescribes to store data in progressively less accurate forms. Data degradation gives users and service providers a fine grained control over the price to be paid, in terms of privacy risks, and to optimize their common interest: balancing privacy and data usability.

I. INTRODUCTION

Privacy has become a popular topic, triggered by the vast amount of web services with an apparently unsatisfiable desire for their users' personal data. Acquiring personal data is big business, a new gold mine for Internet companies, boosting all kind of new web services, increasing the threat to privacy even further [START_REF] Conti | Googling Security: How Much Does Google Know About You[END_REF]. It works; Google can reach over half a billion unique individuals each year, collecting-among many different types of personal data-their search queries, which to an high extent encapsulate their daily lives' habits [START_REF] Cooper | A survey of query log privacy-enhancing techniques from a policy perspective[END_REF]. Google made in 2008 a revenue of $22.1 billion [START_REF] Finance | Key statistics for google inc[END_REF], indicating the worth of that personal data for the company. Google is not alone; in their footsteps many other companies followed, and many will follow.

What do the users get in return for their personal data? Indeed, they profit from all the services which ease their lives. The web has been made accessible thanks to search engines, and communicating with friends and relatives has never been easier. However, until the Internet era, transactions between producer and consumer have been much more transparent for both parties. The consumer pays the price which has been negotiated between producer and consumer, and the producer delivers the good or service. If the price is not satisfactory for both parties, the transaction will not take place. So, it pays off for the producer to be transparent. Today, business models are different. Services are offered for free-in terms of moneyto the user, such that, at first glance, there is no reason to negotiate anymore.

Here the privacy danger becomes apparent. Transparency is one of the key foundations of privacy [START_REF] Jiang | Approximate information flows: Socially-based modeling of privacy in ubiquitous computing[END_REF]; it must be clear for the user how his or her data is being handled, stored, and to whom it will be disclosed. In other words: the price a user has to pay for a service should be expressed in terms of privacy risks, where it was expressed in terms of money in the old days. If the service provider can argue that the data is needed to offer certain kind of services, the user may want to decide to allow the service provider to keep the data longer, paying a higher price, most probably benefiting from a better service.

So why is it a problem that companies store all these data about us? The fact is that, even if we put full trust in the service provider, this data can always be subject to data disclosure due to attacks, corrupt employees, governments demanding the data, et cetera. No access control mechanism has been proven to be both usable and fully secure; to give an example, even servers of the Pentagon [START_REF] Cnn | Hackers stole data on pentagon's newest fighter jet[END_REF] and FBI [START_REF] Post | Consultant breached fbi's computers[END_REF] have been hacked, credit card companies and mobile communication companies have lost personal data on several occasions [START_REF] Vnunet | T-mobile loses 17 million customer details[END_REF]. Moreover, human mistakes are hardly preventable: politicians and policemen lose usb sticks or other media with sensitive information [START_REF] Edri | Uk government loses personal data on 25 million citizens[END_REF], obsolete personal computers sold secondhand are subject to forensic analysis with sometimes shocking results [START_REF] Valli | Oops they did it again: The 2007 australian study of remnant data contained on 2nd hand hard disks[END_REF]. Finally, personal data is often weakly protected by obscure and loose privacy policies which are unjustly presumed to be good and acceptable for a given service.

The privacy violation will only increase with the growth of data which has been collected about us. All these data, even when "legally" obtained by the service providers themselves, foster ill-intended scrutiny and abusive usages justified by business interests, governmental pressures and inquisitiveness among people. Not only criminals and terrorists are threatened. Everyone may experience a particular event (e.g., accident, divorce, job or credit application) which suddenly makes her digital trail of interest for someone else. Moreover, identity fraud is nowadays becoming one of the most serious crimes, with huge consequences for the victims [START_REF] Gordon | Identity Fraud: a Critical National and Global Threat[END_REF]. The retention problem has become so important and the civil pressure so high that practices start changing. For instance, Google and other search engine companies announced to shorten the retention period of their query logs. Limiting the retention of personal data indeed reduces the privacy problems sketched above. Limited retention is a widely accepted privacy principle [START_REF] Agrawal | Hippocratic databases[END_REF], complementary to techniques such as access control, and is included in various privacy regulations [START_REF] Gemeenschap | Richtlijn 95/46/eg van het europees parlement en de raad[END_REF]. The principle prescribes that data should not be stored longer than necessary to fulfill the purpose for which the data has been collected [START_REF] Agrawal | Hippocratic databases[END_REF]. By limiting the time that data is stored, the impact of disclosure of a store is less severe [START_REF] Cooper | A survey of query log privacy-enhancing techniques from a policy perspective[END_REF].

Limited retention is however difficult to put in practice, because it is difficult to determine retention periods; we need to find a mechanism to balance privacy and usability. Otherwise we either end up with a lot of worth for the data collector and zero privacy for the user, or zero worth and full privacy [START_REF] Cooper | A survey of query log privacy-enhancing techniques from a policy perspective[END_REF]. In this paper we propose to put service provider and user at the negotiation table again, with equal rights for both parties. The goals are clear: the service provider wants to optimize the worth of the data it can collect, and the user wants to optimize his or her privacy. We will model those goals and argue to optimize the common interest of both parties. As a result, both parties will agree upon a retention period for storing the users' personal data in exchange for the service offered by the provider.

However, the all-or-nothing behavior of limited retention is too rigorous: after the retention period, the data will be destroyed completely. This makes it hard to balance data usability and privacy. In this paper we also propose a new technique named data degradation, which can bridge the gap between both ends. By using well-known generalization techniques [START_REF] Hilderman | Data mining in large databases using domain generalization graphs[END_REF] (not elaborated in this paper), data degradation degrades the accuracy of the data after predefined retention periods, such that although the usability of that data will decrease, the privacy sensitivity will also decrease. This technique is orthogonal to other privacy enhancing techniques such as k-anonymity [START_REF] Sweeney | k-anonymity: A model for protecting privacy[END_REF]; our aim is not to hide identities, but only to make the knowledge related to those identities less privacy sensitive. We will show that with data degradation, the common interest will be at least equal, but in many situations significantly higher, than which is possible with only limited retention.

In summary, the main contributions of our paper are these:

• a framework to reason about retention periods in order to optimize the common interest of both users and service providers, discussed in Section II; • the technique data degradation, which makes it possible to achieve a common interest that is higher than possible with limited retention, discussed in Section III and analyzed, with explanation of the benefits, in Section IV.

We also indicate the challenges data degradation raises with respect to traditional databases (Section V), mainly by referring to other work. Section VI compares data degradation to related work and indicates where data degradation is complementary to existing techniques. Finally, Section VII gives directions for further work in this new research area of limited retention techniques.

II. FINDING LIMITED RETENTION PERIODS

To let the reader get familiar with our notations and reasoning, we introduce our concepts for the limited retention principle here. In the next section we will reuse and extend those concepts when we refine limited retention to data degradation. Our goal is a qualitative framework which makes it possible to reason about retention periods. We model the interest of a service provider and, separately, the interest of the user in a way that is as simple as possible; a more elaborate model is left for future work. Then we combine these to a common interest, from which an optimal retention period can be derived. We make no quantitative statements about how exactly these interests will be expressed in practice; again, the practical implications will be left for future work.

A. Preliminaries

In the sequel, we consider only one user, which we refer to as "the user". Considering more users is no problem but would lead to the same result if we treat them on equal footing. Also, we shall not make a distinction between different kinds of data: we treat all data on equal footing.

The history of which datum is inserted into the store at which time, is called H; it is a set of pairs of a datum and its insertion time in the store. For example:

H = {(d, b), . . . , (d ′ , b ′ )}
We take H as a constant, we let d range over the set of data and t over the set of time points, and use letter b (birth) for 'the insertion time of a datum into the store'. We use N to measure an age, i.e., length of a time interval, and let a, δ range over ages. A simplifying assumption is that the insert rate of data is constant over time; that is, there exists a constant c such that during each interval of length a the amount of data inserted into the store is c × a:

∀t, a • #{(d, b) : H | t ≤ b < t + a} = c × a (1)
By its definition, limited retention bounds the interval during which a datum is stored by a fixed retention period δ. Hence, the store at time t depends on δ and is expressed as follows:

store(δ, t) = {(d, b) : H | b ≤ t < b + δ}

B. Service provider's interest

In practice, the worth of a datum for the service provider depends on multiple factors, such as the actual content of the datum, the context in which the datum has been acquired, the user from which the datum has been acquired, time of day, and possibly many others depending on the type of use. Our model does not limit the possibility to include those parameters; however, for the sake of simplicity we omit them. Thus, we assume that for the service provider the worth at time t of a datum with birth b depends only on the datum's age t-b:

worth((d, b), t) = wt(t-b)
where wt(a) is nonnegative, monotonic descending in a

The auxiliary function wt is monotonic descending since older data is assumed to be less valuable for the service provider. Figure 1(a) gives the typical shape of function wt(a).

The service provider wants, at each point t in time, to maximize the total worth of the data in the store:

totworth(δ, t) = definition (d,b):store(δ,t) worth((d, b), t) = (d,b):store(δ,t) wt(t -b) =     
A stored datum of age a contributes wt(a) to the sum. Thanks to the constant insert rate (1), the number of stored datums of age a is the same (namely c) for each a = 0, . . . , δ. c × (wt(0) + wt [START_REF] Aggarwal | Vision Paper: Enabling Privacy for the Paranoids[END_REF]

+ • • • + wt(δ)) = c × a=0...δ wt(a)
It turns out that totworth(δ) does not depend on t, hence we omit parameter t and simply write totworth(δ). gives the typical shape of function totworth(δ). Without other constraints, the service provider would achieve his goal by setting δ to infinity.

C. User's interest

For the user it is risky to have data stored at the service provider: in some way or another (hackers' attacks, for instance) the data might be disclosed. We assume that the harm for the user of a disclosure of his data is proportional to the amount of data. Below, we take the risk equal to the amount of data since this simplifies the formulas and nevertheless gives the same results. Again, there are multiple other factors which influence the risk of storing a datum. For example, the fact that a user searched for HIV is more risky than a search for flowers. Also, we assume that disclosure of old data is as harmful as disclosure of recent data. Our model does not restrict the possibility to take those factors into account, but we leave this for future work. Thus we define and simplify risk as follows:

risk(δ, t) = definition #{(d, b) : store(δ, t)} = #{(d, b) : H | b ≤ t < b + δ} = constant insert rate (1) c × δ
So, risk(δ, t) doesn't depend on t and we simply write risk(δ).

The goal of the user is to minimize risk(δ). It is equivalent to maximize the inverse: 1/risk(δ), which we call the privacy guarantee. It follows that the privacy guarantee is infinite when there is no data in the store, δ = 0 (and goes down to zero when retention is unlimited, δ = ∞). To escape mathematical problems (division by zero) and come to a slightly more realistic model of privacy, we apply a smoothing technique so that privacy cannot be infinite: add a constant to the denominator of 1/risk(δ). The smoothing constant s may have a reasonable interpretation; for example, the fact that there is no data in the store, might be interpreted as an indication that "the user has something to hide" and so his privacy guarantee is not infinite. Thus our definition reads:

priv(δ) = 1 s + risk(δ) = 1 s + c×δ Figure 1(c)
gives the typical shape of function priv(δ). The user wants to maximize priv(δ), which without further constraints is achieved by taking δ as small as possible.

D. Common interest

Above we defined the interests of both service provider and user. Those interests are conflicting; whereas the service provider benefits most when δ is large, the user aims for a δ as small as possible. We want a concept of common interest which both parties can agree upon. We expect that the common interest leads to a retention period which is both limited (δ < ∞) and non-zero (δ > 0).

To define the common interest CI(δ), we require two things. First, CI(δ) is proportional to the service provider's goal function totworth(δ) when the user's interest is viewed as constant. Second, CI(δ) is proportional to the user's goal function priv(δ) when the service provider's interest is viewed as constant. Since both goal functions are nonnegative, a suitable function CI(δ) is the product of these: Since worth is monotonic descending and risk is (almost) proportional to δ, it follows that CI( ) has a maximum, which it takes on argument δ opt , say. The existence of a maximum can be interpreted in the following way. By setting the retention period smaller than δ opt the user will gain more privacy, but the common interest will be lower because the decrease in stored data induces a greater loss of total worth for the service provider. Similarly, by setting the retention period larger than δ opt the service provider will gain more total worth of the stored data, but the common interest will be lower because of a larger decrease of the user's privacy.

CI(δ) = totworth(δ) × priv(δ)

III. THE CONCEPT OF DATA DEGRADATION

Recall that the principle of limited retention tries to satisfy the service provider by allowing to store data for at least the retention period δ, and it tries at the same time to satisfy the privacy concern of the user by ensuring that the data is stored for at most the retention period δ (and the previous section shows how to reason about the optimal retention period). The principle is a crude all-or-nothing approach: a datum either exists completely in the store or not at all. The principle of data degradation overcomes the all-or-nothing approach by storing data in progressively less accurate forms so as to make it less privacy sensitive over time while still providing some worth to the service provider. A well accepted form of data degradation is data generalization. This technique is often used in k-anonymity research, where it is used to generalize identities in order to make the data k-anonymous [START_REF] Sweeney | k-anonymity: A model for protecting privacy[END_REF], and applied in data mining and warehousing [START_REF] Hilderman | Data mining in large databases using domain generalization graphs[END_REF]. More about the techniques we use to generalize data-based on domain hierarchies and generalization trees-can be found elsewhere [START_REF] Anciaux | The life-cycle policy model[END_REF].

A. Life-cycle policies

To formalize the data degradation principle, we need some terminology and notation. First, we distinguish several levels of accuracy, say L 0 , L 1 , . . . , L n-1 in decreasing order of accuracy. Level L 0 denotes the most accurate level. Second, the degradation from L i-1 to L i is denoted τ i (τ is mnemonic for "transformation"). Third, the interval from the birth of a datum to its degradation to L i is denoted δ i ; it follows that δ 0 = 0 because a datum is supposed to enter the store in the most accurate level, and we let δ n be the interval from birth to removal from the store. The notation δ abbreviates the sequence δ 1 , . . . , δ n . Almost all this information is captured in a so-called life-cycle, as illustrated in Figure 2.

L 2 δ3 L 1 δ2 L 0 δ1 δ0 Fig. 2.
Graphical representation of a simple life-cycle. Edges denote transitions between Levels of accuracy after a retention period δ.

So, the store consists of data of age at most δ n , and degraded to the appropriate levels:

store( δ, t) = {(d, b):H | b+δ 0 ≤ t < b+δ 1 • (d, b, L 0 )} ∪ {(d, b):H | b+δ 1 ≤ t < b+δ 2 • (τ 1 (d), b, L 1 )} . . . ∪ {(d, b):H | b+δ n-1 ≤ t < b+δ n • (τ n-1 (d), b, L n-1 )}
The product of our framework is a life-cycle policy, which captures how and when data needs to be degraded, and to which the service provider has to comply with. As long as the service provider can be trusted, the life-cycle policy ensures that if the data store is attacked, the impact of disclosure will be less severe; only a small subset of the data will be stored in an accurate form, the rest will be either degraded or destroyed. Some of the technical challenges related to the implementation and enforcement of such policies on traditional database systems are discussed in Section V.

B. Interests revised

We assume that degraded data is less worthwhile for the service provider but also less risky for the user to store, and we will revise the definitions of worth and risk accordingly. The definitions of the total worth, privacy, and common interest in terms of worth and risk remain the same except for the replacement of δ by δ.

Worth and total worth: For the service provider, the worth of a datum depends not only on the age "t-b" but also on the level of accuracy l: worth((d, b, l), t) = wt l (t-b) wt l (a) is nonnegative, monotonic descending in a and l

The effect of degrading a datum to a level of lesser accuracy is that its worth for the service provider is decreased. Indeed, 'wt l (a) is monotonic descending in l' means that for all a: wt 0 (a) ≥ wt 1 (a) ≥ . . . ≥ wt n-1 (a)

As before, the total worth of the store at time t is the aggregation of the worth of all data in the store:

totworth( δ, t) = definition (d,b,l):store(t, δ) worth((d, b, l), t) = c × ( wt 0 (δ 0 ) + wt 0 (δ 0 +1) + . . . + wt 0 (δ 1 -1) + wt 1 (δ 1 ) + wt 1 (δ 1 +1) + . . . + wt 1 (δ 2 -1) + . . . wt n-1 (δ n-1 ) + wt n-1 (δ n-1 +1) + . . . + wt n-1 (δ n -1) ) = c × n-1 l=0 δ l+1 -1 a=δ l wt l (a)
It follows that totworth( δ, t) is independent of t, and we can simply write totworth( δ). Note that the contribution to totworth( δ) of the most degraded level may be negligible in comparison to the less degraded levels: both the age and the level are higher.

Risk and privacy:

The risk of having data in the store was proportional to the amount of data in the store. However, with data degradation, the assumption is that storing data in a more degraded level is less risky than for a more accurate level. To express this, we weight the risk of storing a datum in level L l with a factor r l such that 1 = r 0 ≥ r 1 ≥ . . . ≥ r n-1 . Hence, the definition reads:

risk( δ, t) = definition n-1 l=0 r l × #{(d, b, l ′ ) : store( δ, t) | l ′ = l} = c × n-1 l=0 r l × (δ l+1 -δ l ) Again, risk( δ, t
) is independent of t and we simply write risk( δ). The definition of privacy doesn't change:

priv( δ) = 1 / (s + risk( δ))
Common interest: The goals of both service provider and user remain the same: optimizing the common interest. Except for the replacement of δ by δ there is no change:

CI( δ) = totworth( δ) × priv( δ)
Note that if there is just one level of accuracy, n = 1, the newly defined notions coincide with the already existing notions (such as CI( δ) and CI(δ)) provided we take wt 0 (a) = wt(a) and δ 1 = δ n = δ.

IV. BENEFITS OF DATA DEGRADATION

The aim of data degradation is not to provide more privacy while ensuring the same amount of worth for the service provider, nor providing more worth while ensuring the same amount of privacy as what can be achieved with limited retention. Instead, we want to show that we can achieve an higher common interest with data degradation than with limited retention.

Since the contributions of δ 2 , δ 3 , . . . to CI are nonnegative, data degradation with n > 1 will outperform limited retention:

∀δ 1 , δ 2 , . . . • CI(δ 1 ) ≤ CI(δ 1 , . . . , δ n )
In the following we will show, using experimental examples, that for a set of examples which all comply with our assumptions, the common interest will be higher when we use data degradation. Hence, our target is to give insights in what we can gain by choosing data degradation in comparison to limited retention.

A. Analysis

We use Matlab as the platform for our analysis. We implemented a set of worth and privacy functions which simulate the functions which have a practical shape, complying with the assumptions posed in section II. Given those functions, we let Matlab optimize the common interest considering n = 1 . . . 4 possible degradation steps. Hence, choosing n = 1 means the same as applying limited retention, whereas n > 1 means we allow one or more degradation steps.

Although we can easily experiment with different types of worth functions to simulate a service provider's worth function-as long as they are monotonic descending-we choose the cosine function, as in Figure 3. In the worth functions we use weights w 1 ≥ w 2 ≥ . . . ≥ w n . Recall that r i are weight factors in the risk function. We choose w 1 = r 1 = 1 and w 2 . . . w n vary over 0 . . . 1, similarly for r i .

1) Increased common interest:

To show that using data degradation indeed can result in a higher common interest, we let our script find the δ for which common interest is maximal, with at most n = 4 degradation steps. We choose the following parameters: With those parameters, we obtain the following results (also shown in Figure 4). From this result we conclude for the chosen parameters: 1) Common interests is higher when progressively degrading (n > 1) the data than with limited retention of accurate data (n = 1) 2) With n > 1, δ 1 is smaller than the single δ with n = 1. It thus makes sense to degrade the data earlier to achieve an higher common interest. 3) When n = 4, it turns out that δ 3 = δ 4 , meaning that it is not possible to achieve an higher common interest with more than three degradation steps.

2) A closer look on weights and their effect on CI:

We already have seen that data degradation results in an higher common interest. An interesting question is how the ratios between weights w 1 . . . w n and r 1 . . . r n assigned to each level of accuracy L i influence the gain in common interest which can be achieved with data degradation. Our expectation is the For n = 2, we only vary over δ 2 . The first retention period δ 1 is chosen such that CI(δ 1 , δ 2 ) is optimal for the δ 2 where this plot reaches its highest point. The graph shows that δ 1 for n = 2 is shorter than the optimal δ 1 when n = 1, and that the common interest for n = 2 is higher than for n = 1.

following: with w i descending in i, and w i > r i , degradation results in a higher common interest compared to limited retention.

For this experiment we choose arbitrarily r 1 = 1, r 2 = 0.2, r 3 = 0.1, r 4 = 0.05, and vary over w i . We use a simple distance metric: dist( w) = n i=2 w 1w i to express the drop in worth when degrading the data. Note that the smaller w i is, the less worth is preserved on that accuracy level.

In Figure 5 we plotted for various weights w i the ratio between the common interest which can be achieved with only limited retention and the common interest which can be achieved by at most n = 4 degradation steps. Figure 5(b) shows only common interest points achieved with at least 3 degradation steps. The ratio indicates the fraction of common interest possible with limited retention compared to that of data degradation. Hence, this fraction cannot be higher than 1 since limited retention can never perform better than data degradation (see earlier this section).

We make the following observations:

1) When dist( w) ≃ dist( r) ≃ 2.65, the decrease in worth is similar to the increase in privacy, so that common interest hardly increases (CI ratio is close to 1). 2) When the decrease in worth becomes higher (dist(w) grows), data degradation hardly outperforms limited retention. 3) Most gain in common interest is achieved by applying one-step degradation: data is immediately degraded to a higher level, and fully degraded afterwards. This can be concluded from Figure 5(b), in which only multi-step degradation is allowed.

From the last observation, we conclude that for some (our current) parameter values it makes sense to generalize the data before storing it. Hence, a higher common interest can be achieved by handing over some accuracy to get much more privacy in return, especially when the privacy increase is much higher than the loss in worth. When the decrease in worth is more close to the increase in privacy, we showed again that multi-step degradation leads to an increase in common interest. 5. Ratio between the common interest which can be achieved with only limited retention and the common interest which can be achieved by allowing n = 4 degradation steps. A ratio equal to 1 means no increase in common interest, lower than 1 indicates that limited retention performs less than data degradation. This ratio can never be higher than 1; limited retention cannot perform better than data degradation.

V. IMPACT ON TRADITIONAL DATABASES

Since data degradation will take place within the data store, the first and legitimate question which comes in mind is how complex will the technology be to support it. Traditional databases are developed in order to efficiently and durably insert data, make updates to this data and to query them. The ACID properties ensure durability and correct execution of queries and updates to keep the database in a consistent state. In our context, we have to revise the ACID properties, such that the durability requirement ends when data should be degraded, even requiring that data can not be recovered. In order to apply the degradation model correctly, we have to ensure that even operations made outside control of the DBMS cannot reverse the degradation steps. A delete statement should therefor not only be made visible on the application level, but it should also be internally irreversible. However, in current traditional database implementations this is not the case [START_REF] Miklau | Securing history: Privacy and accountability in database systems[END_REF].

Hence, identifying the impact of making a database datadegradation aware leads to several important questions. For more details, and a comprehensive study on this topics, we refer to [START_REF] Anciaux | The life-cycle policy model[END_REF]. Here, we only give a brief overview. Updates to the data items due to degradation, as well as final removal from the database have to be enforced. As pointed out in [START_REF] Miklau | Securing history: Privacy and accountability in database systems[END_REF], traditional databases cannot even guarantee the non-recoverability of deleted data due to different forms of unintended retention in the data space, the indexes and the logs. With data degradation, the problem is particularly acute considering that each data item inserted in the database undergoes multiple degradation steps. The storage of degradable attributes, indexes and logs have to be revisited in this light.

Databases have been designed to speed up queries. Some workloads induce the need of few indexes on the most selective attributes to get the best trade-off between selection performance and insertion/update/deletion cost. For other application, insertions are done off-line, queries are complex and the data set is very large. This leads to multiple indexes to speed up even low selectivity queries. Data degradation can be useful in both contexts. However, data degradation changes the workload characteristics in the sense that queries become less selective when applied to degradable attributes. This introduces the need for indexing techniques supporting efficiently degradation.

User transactions inserting data items with degradable attributes generates effects all along the lifetime of the degradation process, that is from the transaction commit up to the time where all inserted data items have reached a final state for all their degradable attributes. This significantly impacts transaction atomicity and durability and even isolation considering potential conflicts between degradation steps and reader transactions.

VI. RELATED WORK

Proposals have been made to make the donor herself responsible for protecting her own data, not just relying on the philosophy of trusting organizations to protect the privacy. In their vision paper, Aggarwal et al proposed the P4P framework [START_REF] Aggarwal | Vision Paper: Enabling Privacy for the Paranoids[END_REF], in which the donor keeps control about which information to release to service providers. They consider the 'paranoid' user who doesn't trust the collecting organizations, in contrast to the users of policy based framework such as P3P [START_REF] W3c | Platform for privacy preferences (P3P) project[END_REF]. Although such solution is robust against server attacks, the accessibility for service providers is much lower, leading to high communication costs when data needs to be queried or updated and placing constraints on how applications are developed and deployed. In contrast, data degradation doesn't place these restrictions on applications, data can still be stored at the server side and by enforcing data degradation, donors are in control of the level of privacy risk they want to take in terms of retention periods.

Limited retention techniques have been proposed to ensure that data can no longer be subject to occasional disclosure. The limited retention principle is a key principle behind many privacy laws, and as such has been adopted in the work on Hippocratic databases [START_REF] Agrawal | Hippocratic databases[END_REF]. Implementation frameworks behind such a system, based on access control mechanisms, including handling of generalized forms of the data, already exist [START_REF] Agrawal | Hippocratic databases[END_REF], [START_REF] Bertino | Privacy-preserving database systems[END_REF], [START_REF] Byun | Micro-views, or on how to protect privacy while enhancing data usability: concepts and challenges[END_REF]. However, such systems are still based on trust [START_REF] Aggarwal | Vision Paper: Enabling Privacy for the Paranoids[END_REF]; trust which cannot be put forever on a system. Even when secure access control techniques are used, a database administrator can be or become malicious. Even if the chosen security regime can be proven successful now, it might be not in the future [START_REF] Gordon | CSI/FBI computer crime and security survey[END_REF].

In addition to access control, security measures for protecting a database server, such as data encryption, firewalls, and intrusion detection systems can be used. Those techniques make attacks more difficult without completely preventing them. Intrusion detection systems [START_REF] Cavusoglu | The value of intrusion detection systems in information technology security architecture[END_REF] are especially useful against repetitive attacks such as spying on a database, although it is still hard to find a good balance between false negative and false positive detections. However, used in addition to data degradation, intrusion detection systems would make it very hard for even a determined attacker to obtain a large consecutive history of accurate data.

Anonymization of data might be a solution to prevent disclosure of privacy sensitive-data. In fact, major companies such as Google already state they will adopt anonymization as a measure to improve privacy protection [START_REF] Fleischer | Taking steps to further improve our privacy practices[END_REF]. k-Anonymity [START_REF] Sweeney | k-anonymity: A model for protecting privacy[END_REF] is based on the idea of masking (parts of) the (quasi) identifier of a partly privacy-sensitive tuple, such that the sensitive part of the tuple will be hidden between k -1 potential identifier candidates within the same dataset. The work on l-diversity [START_REF] Machanavajjhala | l-diversity: Privacy beyond k-anonymity[END_REF] goes a step further by taking background knowledge into account, enforcing enough diversity between the privacy sensitive attributes. Although data degradation aims not on generalizing identities, the techniques used to generalize the data are comparable.

Usually, anonymization is applied to large datasets at once, making sure that for each tuple, the tuple shares the same identifier with k-1 others. In practice this could result in a strictly k-anonymous database at the cost of losing much usability. Although Byun et al provided a technique to update anonymized databases [START_REF] Byun | Secure anonymization for incremental datasets[END_REF], each time new data arrives, the database has to be sanitized into a k-anonymous state again, making it hard to obtain a clear view of the database from an application perspective, since old tuples might be sanitized at unpredictable times. Besides, correctly anonymizing the data is a hard problem [START_REF] Meyerson | On the complexity of optimal kanonymity[END_REF]. To illustrate, a good example of incorrect and insufficient use of anonymization has been given when American Online decided to put a large set of search queries online [START_REF] Hillyard | Issues around the protection or revelation of personal information[END_REF]. AOL anonymized the IP addresses of the computers from which the queries were issued, which was not enough to prevent attackers from inferring many privacy sensitive facts. Data degradation does not suffer from such kind of vulnerabilities, because the facts themselves are degraded, not the identities.

Data degradation can be complementary to all discussed techniques. Firstly, by limiting the impact of inevitable privacy breaches, data degradation is complementary to access control, since data which has been subject to degradation either has a lower level of accuracy and thus sensitivity, or has already been removed from the system. Moreover, although only on temporary basis, accurate data can still be protected against regular attacks with the use of access control techniques. Secondly, anonymization is good practice when datasets have to be made public without revealing the identity of the users; for example, when used for disclosing datasets for research purposes, and therefore it can be a complementary technique to data degradation. Data degradation is particularly useful when data needs to be accurate for some time, but where the details are less important when the data gets old. Moreover, within our degradation model, the identifier of the user can be kept intact; hence, user-oriented services can still exploit the information to the benefit of the user.

VII. FUTURE WORK AND CONCLUSION

In this paper, we presented a new framework to be able to reason about limited retention, and more specifically data degradation. We showed that using a function capturing the worth of storing data for the service provider, and a function for the risk involved for storing this data, we are able to find an optimal retention period such that the common interest for both service provider and user is respected.

The increase of common interest which can be achieved depends on multiple factors. When the decrease of worth is much higher than the increase in privacy, it is good to apply data degradation, even if this means that the data will never be stored accurately. In our analysis, we showed cases where it is useful to progressively degrade the data from accurate states to generalized states until final destruction of the data. We hope that our findings will lead to a break with a tradition where service providers collect everything they can, and store it as long as they can. This paper is a first attempt to achieve that goal. However, there is much future work to do:

• Firstly, we have to investigate the actual privacy guarantees of data degradation and its effect on the risk functions. How much 'privacy' can be provided by generalizing the data is an important question in order to correctly define the risk functions. Correctly defining privacy will always be subject of discussion, mainly because of its subjective nature. In that light, defining privacy as function based on risk, and to relate this risk to the amount of stored data is a promising first step.

• We introduced a model in which we only take retention periods of data into account. Indeed, both privacy and worth function can depend on various parameters other than retention periods. Users can choose different 'risk profiles' depending on the nature of the data items, and service providers can attach a lower worth to specific data based on the user, location, time of the day, et cetera. • An important question is if and how service providers will be able to express their worth functions. To put our framework into practice, it is necessary to provide tools which enable service providers to give transparency about their need to collect personal data. To conclude: our framework makes that it is indeed possible to reason about retention periods such that not only one but both parties will be satisfied. Together with data degradation, we presented a promising new approach to close the huge gap between the enormous amount of collection and storage of personal data, and the risk users have to take to be able to profit from all the new and exciting services offered to us today and in the future.

Figure 1

 1

Figure 1 (

 1 Figure 1(d) gives the typical shape of function CI(δ).Since worth is monotonic descending and risk is (almost) proportional to δ, it follows that CI( ) has a maximum, which it takes on argument δ opt , say. The existence of a maximum can be interpreted in the following way. By setting the retention period smaller than δ opt the user will gain more privacy, but the common interest will be lower because the decrease in stored data induces a greater loss of total worth for the service provider. Similarly, by setting the retention period larger than δ opt the service provider will gain more total worth of the stored data, but the common interest will be lower because of a larger decrease of the user's privacy.

Fig. 1 .

 1 Fig.1. The common interest function reaches its highest point at δopt meaning that a retention period of δopt gives the best balance between totworth for the service provider and priv for the user.

Fig. 3 .

 3 Fig. 3. wt l (a) = w l × 1 + cos a × π 180 , 0 ≤ a ≤ 180

Fig. 4 .

 4 Fig. 4. Common interest function for n = 1 (limited retention) and n = 2.For n = 2, we only vary over δ 2 . The first retention period δ 1 is chosen such that CI(δ 1 , δ 2 ) is optimal for the δ 2 where this plot reaches its highest point. The graph shows that δ 1 for n = 2 is shorter than the optimal δ 1 when n = 1, and that the common interest for n = 2 is higher than for n = 1.

  Fig.5. Ratio between the common interest which can be achieved with only limited retention and the common interest which can be achieved by allowing n = 4 degradation steps. A ratio equal to 1 means no increase in common interest, lower than 1 indicates that limited retention performs less than data degradation. This ratio can never be higher than 1; limited retention cannot perform better than data degradation.


Shouhei Hanaoka
email: hanaoka-tky@umin.ac.jp

Yoshitaka Masutani

Mitsutaka Nemoto

Yukihiro Nomura

Takeharu Yoshikawa

Naoto Hayashi

Naoki Yoshioka

Kuni Ohtomo

Probabilistic Modeling of Landmark Distances and Structure for Anomaly-proof Landmark Detection

Keywords: Landmark, Combinatorial optimization, MAP estimation, Anatomical anomaly, Computed tomography, Spine

A combinatorial optimization algorithm for detecting multiple anatomical landmarks is presented. It can determine the positions of over 100 landmarks concurrently, taking spatial correlations of all landmark pairs into account. Provided that a set of landmark candidate lists is given by sensitivity-optimized single-landmark detectors, the proposed algorithm can find the most probable combination of them through solving a MAP estimationbased combinatorial optimization problem. Additionally, it is designed to handle subjects with "segmentation anomaly of the spinal column," a common anatomical anomaly of the spine. The proposed system was evaluated with 156 landmarks in 50 datasets, using virtually created detector output sets. In the result, the algorithm achieved 97.6% of spinal anomaly estimation accuracy even with 50 points of candidates given per landmark, as well as 96.2% of accuracy in landmark candidate selection. From these results, usefulness of the proposed algorithm for subjects with spinal anomaly was suggested.

Introduction

Landmark point detection algorithms are extensively researched and widely used in various medical image processing applications. However, it is a difficult task to detect a large number of landmarks correctly, because the human body includes a lot of similar regions sharing their appearances. Even most of frequently-used landmarks do not have truly unique local shape or intensity. Consequently, detection results often include a certain number of false positive candidates.

Furthermore, some sort of important anatomical landmarks can intrinsically be non-existent, not only in patients due to pathological condition, but even in healthy subjects due to anatomical anomalies (fig. 1a). Confirming such a situation is very difficult by detecting each landmark independently and sequentially.

One solution is to determine the entire landmark positions simultaneously, taking their spatial relation into account. It can be done by dividing the whole problem into 2 sequential phases: the individual landmark detection phase and the combinatorial optimization phase (figs. 1b, 1c). In the former phase, each landmark is detected by a single detector. Each detector is optimized to maximize the sensitivity, not the specificity, so the detection result forms a candidate list which includes a lot of false positives. In the latter phase, the best combination of choice from all candidate lists is selected with the use of a priori knowledge on the inter-landmark relationship. Defect of any landmark can also be detected through this optimization phase. (b) A schema of the framework of the proposed method. Note that, for the non-existent landmark (□), all detected candidates are rejected by the following combinatorial optimization.

Among related studies, Seifert et al. reported a framework to automatically detect 19 landmarks. [START_REF] Seifert | Hierarchical parsing and semantic navigation of full body CT data[END_REF] In their approach, a belief propagation algorithm is used with prior knowledge about landmarks' geometric relationships such as "to the right of", "close to", etc. Though their result was excellent, the prior knowledge used and the way to build it were not described precisely.

In this paper, we introduce a framework to determine over 100 landmark positions reliably. Especially, we focus upon our method to solve the landmark-set combinatorial optimization problem with use of maximum a posteriori (MAP) estimation. The method was evaluated with 156 landmarks in 50 human body CT image datasets. Virtually-created detector outputs were used in this study, rather than those of real detectors [START_REF] Nemoto | A unified framework for concurrent detection of anatomical landmarks for medical image understanding[END_REF], in order to evaluate the ability of our framework to handle a large number of false candidates. Furthermore, handling of common anomaly (alteration in the number of thoracic/lumbar vertebrae) was also evaluated.

Methods

Definition

2.1.1. Detector output. Firstly, we defined a landmark detector generally. In our framework, each detector has to output not only a series of candidate positions, but also estimated probabilities of them. Based on this, outputs from one detector (for the m-th landmark) are defined as:

N m number of the candidates (0≤ N m ) c m i coordinates of the i-th candidate (1≤i≤N m ) p m i detector-estimated probability of the i-th candidate (0 < p m i ≤ 1) p m  probability that "none of candidates c m i is correct" (0 ≤ p m  ≤ 1)
The whole set of output is also defined as

I m =    m N m m N m m m p p p N m m , , , , , , , 1 1 
 c c . Any output I m must satisfy 1     i m m p p .
Letting the true landmark position be x m , the detector-estimated probability p m i can be interpreted as a conditional probability with a certain detector output I m . ." The proposed algorithm can consider such a situation, and handles it as one extra state which will be represented as "

 

  m x
" in this paper. The state  can be a "true answer" when the target landmark is out of the imaging range, or corresponding landmark does not exist anywhere (e.g., due to some anatomical anomaly or pathological condition). In this describing method, the probability p m  can also be interpreted as a conditional probability as

  m m m p p I |     x . (2)
2.1.2. The prior probability distribution of LM positions. Secondly, the prior probability function for all possible landmark position sets has to be defined in advance. Let x 1 , x 2 , …, x M be the positional vectors of total M landmarks, and X = (x 1 t x 2 t … x M t ) t be the concatenated form of them. Each element x m can be regarded as a stochastic variable whose domain is

  m N m m m c c x , , , 1   
. The aim is to approximate the prior probability distribution of X as a single function p(X), in order to use in the following MAP estimation.

In this study p(X) is defined as a function of the squared distances between all landmark pairs. When any x m (1≤m≤M) satisfies x m =, however, the corresponding term of p(X) is replaced by a conventionally-defined constant term. The calculation method of p(X) without considering  is discussed in the next chapter, which will be followed by the general definition of p(X).

(i) When x m ≠, m. Let d i,j = 2 j i x x  (1≤i<j≤M) (3)
be the squared distance between the i-th and j-th landmarks. Note that d i,j is defined only if both x i and x j are not . Then, the distance is normalized by its average E(d i,j ) and variance V(d i,j ) in the training datasets. The normalized distance g i,j follows the equation

      M j i d d d g j i j i j i j i j i          1 . , if V E , , , , x x (4) Let a vector   t ,M M j i , , g g g g 1 , 3 1 2 1     G
be the concatenated normalized squared distances between all landmark pairs. Note that the vector G has

  2 1 C 2   M M M of elements.
Then, the prior probability distribution p(X) was approximated by a multivariate normal distribution of the vector G. That is,

            G V G V X 1 C 2 1 exp ) 2 ( 1 2 t M p  (5)
where V represents the covariance matrix of G, which is calculated from training sets, and V is the determinant of V. The size of the matrix

V is M C 2 × M C 2 .
(ii) General definition. Eq. ( 5) can be written in an extended form

                                  1 1 1 1 1 , , , , 1 1 ,
C 2 1 exp ) 2 ( 1 2 M i M i j M k l k l k j i M k l j i g g p M V V X  (6)
where {V -1 } (i,j),(k,l) is the corresponding element of the matrix V -1 . Here, the summed term

      l k l k j i j i g g , , , , 1 ,    V
is not available when any of x i , x j , x k or x l is . In order to define p(X) generally, we replaced them as follows:

                                              ) ,
, , C 2 l k j i l k j i g g t t p l k j i l k l k j i j i l k j i M i M i j M k M k l l k j i M   x x x x V V X (7) 1 , ) , ( ), , ( 1 1 1 1 1 1 ,
When (i,j)≠(k,l), the summed term t (i,j),(k,l) evaluates how g i,j and g k,l are correlated as like that of training sets. When no prior information on g i,j (or g k,l ) is available, the expected value of this term should be zero. That is why it should be replaced by zero in case of .

On the other hand, when (i,j)=(k,l), the term t (i,j),(i,j) = { V -1 } (i,j),(i,j) • g i,j 2 always has a positive value which evaluates how the distance d i,j varies from that of the training sets. Therefore, we replaced it by a positive constant γ. Because increasing this parameter will reduce the probability p(X) for any x m =, γ can be regarded as a parameter which controls how the algorithm avoids the state , or how it prefers any detected candidate c m i , for each landmark. We empirically selected γ=2 based on our preliminary experiments (the data is not shown). 2.1.3. a posteriori probability distribution. Once the detector outputs are given, the probability a posteriori can be calculated with Bayes' theorem. From the series of detector outputs I 1 , I 2 , …, I M and the prior probability distribution p(X), the posterior probability can be calculated by the Bayes' as:

        M M M
3 2 1 3 2 1 X X X   (8)
The denominator is constant and independent of X. Therefore, the maximum a posteriori estimation of landmark position set X is as follows:

      X X X X X X p p p M M    | I ,..., I , I , I max arg I ,..., I , I , I | max arg ˆ3 2 1 3 2 1 (9)
We assumed that all detector outputs are independent of each other, as well as independent on the positions of the other landmarks. It means that each detector output is only dependent on the corresponding landmark position. Then the term

  X | I ,..., I , I , I 3 2 1 M p in (9) can be divided into the product of single-landmark conditional probabilities p(I m | x m ). That is,                  M m m m m m M m m m M p p p p p 1 1 3 2 1 I I | | I | I ,..., I , I , I x x x X . (10)
The term p(I m ) in this formula is independent of X, so it can be ignored through the MAP estimation. The denominator p(x m ) is also can be ignored, because the term is the probability distribution of one single landmark position x m without any prior information. Ignoring the term is equivalent to regarding it to be constant and homogeneous anywhere. The remaining term p(x m | I m ) is the detector-estimated probability as described in Eqs. ( 1) and (2). That is,

           m m i m m i m m m p i p p x c x x if , if I | . (11)
In conclusion, the MAP estimation will be performed by the following formula:

                       X x X X X X p p p M m m m M 1 3 2 1 I | max arg I ,..., I , I , I | max arg ˆ (12)
2.2. Implementation 2.2.1. Tikhonov's regularization. In practice, the dimension of G in Eq. ( 5) (or Eq. ( 7)), which equals

  2 1 C 2   M M M
, can be much greater than the number of training cases N. In such a case, the covariance matrix V will be rank-deficit and has no inverse matrix V -1 . To avoid them, Tikhonov's regularization [START_REF] Crimi | Prior knowledge regularization in statistical medical image tasks[END_REF] was performed by replacing V with its regularized matrix V reg as follows:

I V V    reg (13)
We empirically selected λ = 1.0 for this study.

Note that the calculation of Eq. ( 5) can be speeded up by using the following formula (derived from the Woodbury matrix identity),

      U I D U I I DU U I V V 1 1 1 2 1 1 1 1 '                       t t reg . ( 14
)
Here, D is a (N-1)×(N-1) diagonal matrix whose diagonal elements are nonzero eigenvalues of V, and U is a (N-1)× M C 2 matrix whose row vectors are corresponding eigenvectors. I and I' are identity matrices whose sizes are M C 2 and N-1, respectively.

Combinatorial optimization.

The maximization of Eq. ( 12) for all combinations of X={x m },

  m N m m m c c x , , , 1   
is performed by a Gibbs' samplerbased simulated annealing algorithm reported by Geman & Geman [START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF]. In the algorithm, a virtual temperature T is introduced to modify the probability distribution. The distribution in (12) was modified as follows:

      T M m m m p p Z T p 1 1 modified I | 1 ;                   X x X (15)
Z is a normalization factor in order to make the sum of probability 1.

When T=1, the modified distribution equals to the original in (12). With higher T, the modified distribution is almost homogenous for any X in its domain. However, it becomes sharper and more local with lower T. In the limiting case of T→+0, the modified distribution has nonzero probability only at its global maximum point. (Fig. 2) Fig. 2. A simple example how the probability distribution is modified by the temperature. The green curve illustrates an original distribution (a unary mixture Gaussian distribution for example). The modified distributions with T = 10 and 0.1 are illustrated with red and blue curves, respectively. Note that the distribution becomes flat with higher temperature, while it becomes sharper with lower temperature.

In the simulated annealing, each x m (m = 1, 2, …, M) is sequentially and repeatedly sampled from the modified distribution by a Gibbs' sampler. The sampling begins with very high T, which decreases gradually, and finally it gets so low that the system converges to the maximum point. In this study the simulated annealing were performed for 100 cycles with T=1000 (in order to cancel the effect of the initial condition; so-called burn-in), then it was gradually cooled down to T=0.01 through 1000 cycles.

Handling of vertebral anomalies.

A majority of human beings have 12 thoracic and 5 lumbar vertebrae. Segmentation anomaly of the spinal column is a common anatomical anomaly in which the subject has 11 or 13 thoracic, and/or 4 or 6 lumbar, vertebrae. The prevalence is, in a report, about 9 % [START_REF] Carrino | Effect of spinal segment variants on numbering vertebral levels at lumbar MR imaging[END_REF]. This anomaly is very problematic in both defining and detecting vertebral landmarks (fig. 1a). The "state " approach described in chapter 2.1. is not enough for them, because it is not only a local banishment of a single anatomical entity but causing a global morphological change in the spine.

To overcome this, a series of "anomaly landmark position set converters" are introduced. One converter can convert any landmark position set in a subject with a certain type of anomaly (e.g., 6 lumbar vertebrae, or "6L") into a virtually normalized landmark position set (i.e., a landmark position set as if she or he has only 5 lumbar vertebrae). It is simply performed by replacing each landmark coordinates by an appropriate internally dividing point between two of them (fig. 3). Because 7 types of anomalies (11T, 13T, 4L, 6L, 11T+6L, 13T+4L) were considered in this study, 7 different converters were designed. Through one of these converters, abnormal spines can be converted into a "normalized" one, with which the prior probability p(X) can be calculated. computing processor. The latter was utilized to speed up the calculation. The processing times were 9.7 minutes per case for comparison after all optimization strategy and 3.6 minutes for comparison in situ strategy with 100 candidates per landmark.

5 times of experiments were performed for each experimental setting, using different sets of virtual detector outputs.

The statistical model of squared distances (the mean E(d i,j ), variance V(d i,j ) and the covariance matrix V) was calculated by leave-one-case out method.

Results and discussion

The proposed algorithm was evaluated with 156 bony landmark points which were manually inputted in 50 thin-slice clinical human body CT datasets. The list of landmark used is available in [START_REF] Nemoto | A unified framework for concurrent detection of anatomical landmarks for medical image understanding[END_REF]. Among the 50 datasets, total 9 had a segmentation anomaly (three 13T, one 11T, three 6L, two 4L and none of 13T+4L or 11T+6L).

The summary of results is shown in Table 1. An example result is also shown in Figure 4. The accuracies of anomaly estimation were from 87.6 to 97.6%, varying among strategies and number of candidates. The accuracy of anomaly estimation was better in comparison after optimization strategy than in comparison in situ strategy. Though the former strategy takes approximately 3 times longer time, it seems to be useful especially for a larger number of candidates. The optimization result of each landmark in each case was also evaluated. Any landmark which existed in reality was classified as one of TP, FP candidate or FN. Within them, TP means that the algorithm correctly selected the true position candidate. FP candidate means that the algorithm selected any false position candidate. If the algorithm concluded that the target landmark did not exist, but in reality it existed, it was categorized as FN. In this study, the accuracy ratios (ratios of TP) for existent landmarks were varied from 89.64 to 97.65%.

Additionally, all of non-existing landmarks were classified as either TN or FP existence . TN means that the target landmark did not exist in reality and the algorithm correctly found it (either by adopting some anomaly hypothesis which does not include the target landmark, or by selecting the state  for the landmark.) If the algorithm chose any candidate other than , it was classified as FP existence . In this study all of non-existent landmarks were determined as TN, so the accuracy was 100%.

It is noticeable that, throughout all experimental condition, none of non-existing landmarks were mistaken as FP existence . In fact, without any exception, all mistakes in anomaly estimation were underestimation of thoracic / lumbar vertebral number. It may imply that our algorithm was not sufficiently optimized for the problem, favoring anomalies having less vertebral number, and yet to be investigated in the future work. Also, we are now planning to overcome this limitation by adding some additional term, which evaluates the regularity of vertebral bones' alignment, to the MAP estimation.

Conclusion

A novel combinatorial optimization algorithm for landmark detection and anomaly estimation was presented. The proposed method showed fair results even with a large number of landmark position candidates. Additionally, its feasibility to estimate segmentation anomaly of the vertebrae, which is one of the most common and problematic anomalies in detecting bone landmarks, has been shown. Therefore, we believe that our algorithm is useful in medical image analysis such as a pre-process for computer-assisted detection/diagnosis (CAD) applications. The future work will include improvement of accuracy in anomaly estimation and evaluation with real detector outputs with a large number of datasets.

Fig. 1 .

 1 Fig. 1. (a) An example of anatomical anomaly and landmark deficit in the 5th lumbar vertebra. (b) A schema of the framework of the proposed method. Note that, for the non-existent landmark (□), all detected candidates are rejected by the following combinatorial optimization.

   >0, it indicates that there are some possibilities of "no true landmark position is included in the candidate set  

Fig. 4 .

 4 Fig. 4. An example result of a case with 6 lumbar vertebra (6L) anomaly. (a) The false (green) and true (orange) candidates outputted by the virtual detector for the tip of right transverse process of the 4 th lumbar vertebra. (b) The optimization result for 156 landmarks.

Table 1 .

 1 The entire result of anomaly estimation and combinatorial optimization.

			No. of	anomaly estimation accuracy (%)	landmark determination: existent-in-reality LMs (%)	non-existent LMs (%)
			cands	anomaly cases	normal cases	overall	TP	FN	FPcandidate	TN	FPexi-stence
	comparison	after opt.	25 50 75	88.9 ±13.6 88.9 ±7.9 86.7 ±14.5	99.5 ±1.1 99.5 ±1.1 99.0 ±1.3	97.6 ±3.3 97.6 ±1.7 96.8 ±2.7	96.60 ±2.90 96.25 ±0.34 94.40 ±0.39	0.09 ±0.13 0.08 ±0.05 0.09 ±0.08	3.31 ±2.78 3.67 ±0.32 5.51 ±0.36
			100	77.8 ±11.1	94.6 ±4.0	91.6 ±3.0	91.33 ±1.03	0.26 ±0.09	8.42 ±1.00	100	0
	comparison	in situ	25 50 75	84.4 ±12.7 73.3 ±18.6 64.4 ±19.9	100 ±0 100 ±0 99.5 ±1.1	97.2 ±2.3 95.2 ±3.4 93.2 ±3.9	97.65 ±0.36 95.67 ±0.70 93.37 ±1.02	0.09 ±0.08 0.17 ±0.12 0.24 ±0.15	2.26 ±0.32 4.15 ±0.59 6.40 ±0.91	±0	±0
			100	48.9 ±14.9	96.1 ±2.8	87.6 ±2.6	89.64 ±0.61	0.40 ±0.05	9.96 ±0.57

Acknowledgement

This study is a part of the research project "Computational Anatomy for Computeraided Diagnosis and therapy: Frontiers of Medical Image Sciences", which is financially supported by the grant-in-aid for scientific research on innovative areas MEXT, Japan.

Anomaly conversion

Fig. 3. The anomaly conversion. In this example, (a) a landmark position set in a case of 4lumbar vertebra (4L) anomaly is virtually converted into (b) a "normal" landmark set with 5 vertebrae. With this conversion, the posterior probability p(X) of the given landmark position set can be calculated in the same manner as in the cases with (c) normal 5-lumbar spines.

Because the algorithm do not know which anomaly is correct, it must estimate it. It can be done by (i) hypothesizing one anomaly (or normal), (ii) calculating the posterior probability under the hypothesis and (iii) comparing the probabilities between all hypotheses. In detail, two different strategies were evaluated: (i) Comparison after all optimization. Firstly, a series of combinatorial optimization with all converters are performed. Then, the converter with the largest probability is chosen as the estimated anomaly (or normal). (ii) Comparison in situ. The comparison is performed at the end of every cycle in the simulated annealing. In other words, the state "which converter is currently selected" is also dealt with as one extra variable to be optimized.

Evaluation

Virtual detector output construction.

The detector outputs were virtually created for each landmark and for each CT dataset. If the target landmark existed in the subject's real body, it was included as a candidate, as well as 25, 50, 75 or 100 of false positive candidates. Each false positive candidate was determined randomly following a 3-D Gaussian probability distribution (the center corresponds to the true point, and the standard deviations were σ x = σ y = σ z = 3 1

•100 millimeters). Any false candidates within 20 millimeters from the true point were removed. If the target landmark did not exist in reality, only false positive points were added, using another adjacent landmark as the distribution center. The detector-estimated probability p m  was fixed to be 0.05, and all of p m i were set to be uniform.

Experimental settings.

The experiments were performed on a workstation with an Intel ® Core™ i7-2600 processor and one NVIDIA ® Tesla™ C2050 GPU


Shuyan Liu

Johan Segers

Multivariate MA(∞) processes with heavy tails and random coefficients

Keywords: extremes, heavy tails, regular variation, tail process

Many interesting processes share the property of multivariate regular variation. This property is equivalent to existence of the tail process introduced by B. Basrak and J. Segers [1] to describe the asymptotic behavior for the extreme values of a regularly varying time series. We apply this theory to multivariate MA(∞) processes with random coefficients.

Introduction

We interested in multivariate infinite-order moving average process with random coefficient matrices, defined for t ∈ Z by (1.1)

X t = ∞ i=0 C i (t)ξ t-i ,
where (ξ t ) t∈Z is a sequence of i.i.d. regularly varying random vector in R q and {C i (t), i ≥ 0, t ∈ Z} is an array of random d×q matrices. The tail behavior of such process is usually controlled by a moment condition on the matrix C i (t), this has been shown by Hult and Samorodnitsky [START_REF] Hult | Tail probabilities for infinite series of regularly varying random vectors[END_REF]. We adapt the conditions as follows.

INTRODUCTION

Condition 1.1 Suppose that P( i≥0 { C i (0) = 0}) = 0, and there is some ε ∈ (0, α)

such that i E C i (0) α-ε < ∞ and i E C i (0) α+ε < ∞, if α ∈ (0, 1) ∪ (1, 2); (1.2) E i C i (0) α-ε α+ε α-ε < ∞, if α ∈ {1, 2}; (1.3) E i C i (0) 2 α+ε 2 < ∞, if α ∈ (2, ∞); (1.4)
To apply the results in Hult and Samorodnitsky [START_REF] Hult | Tail probabilities for infinite series of regularly varying random vectors[END_REF], a 'predictability' assumption is required.

Condition 1.2 Suppose there is a filtration (F j , j ∈ Z) such that for all t ∈ Z and i ∈ N ∪ {0} (1.5) C i (t) ∈ F i-t , ξ t-i ∈ F i-t+1 , (or ξ j ∈ F 1-j ),

(1.6) F j is independent of σ(ξ -j , ξ -j-1 , . . .).

If the sequence (C i (t)) is independent of the sequence (ξ i ), we can set F j = σ((C k (t)) k≥0,t∈Z , (ξ k ) k≥1-j ).

Throughout the paper we will use the following notations, E q = [-∞, ∞] q \{0}, E q u = {x ∈ E q | x > u} and

S q-1 = {x ∈ R d | x = 1}. A function f is called regularly varying at infinity with index τ ∈ R, if lim x→∞ f (ux)/f (x) = u τ for all u > 0.
Definition 1.1. A random vector X in R d is regularly varying with index α > 0 if there exists a non-null Radon measure µ on E d and a regularly varying function V of index -α such that, as x → ∞,

(1.7) P x -1 X ∈ • V (x) v -→ µ(•),
where v -→ denote the vague convergence of measures. Relation (1.7) is equivalent to

P X > ux, x -1 X ∈ • P { X > x} v -→ u -α P {Θ ∈ •} , as x → ∞. The law of Θ is called spectral measure.
The random sequence (X t ) given by (1.1) is built from two components which satisfy the following basic assumptions.

A1

The q-dimensional random vector ξ t is regularly varying of index α ∈ (0, ∞) and with spectral measure L(Θ) on S q-1 , i.e. there exists a non-null Radon measure µ on E q such that, as x → ∞,

P x -1 ξ 0 ∈ • P { ξ 0 > x} v -→ µ(•),
where the measure µ has the following representation, for all u > 0,

µ x x > u, x x ∈ • = u -α P {Θ ∈ •} .
A2 The array of random matrices {C i (t), i ≥ 0, t ≥ Z} of dimension d×q is independent of the sequence (ξ t ) and stationary when index over t ∈ Z.

Joint regular variation

The following theorem is adapted version of Theorem 3.1 in [START_REF] Hult | Tail probabilities for infinite series of regularly varying random vectors[END_REF].

Theorem 2.1 [START_REF] Hult | Tail probabilities for infinite series of regularly varying random vectors[END_REF] Suppose that Condition 1.1 and 1.2 hold, then the series X t given by (1.1) converges a.s. and

(2.8) P x -1 X 0 ∈ • P { ξ 0 > x} v -→ ∞ i=0 E µ • C i (0) -1 (•) , as x → ∞ on E q .
Observe that by assumption A1, the function u → P { ξ 0 > u} is regularly varying with index -α < 0. Therefore if the limit on the right-hand side in (2.8) is a non-null Radon measure, then vector X 0 is regularly varying with index α. As we shall see slightly stronger condition is needed to ensure regular variation of vector X 0 . Consider for u > 0 Hence Condition 1.2 holds for the sequence Ci (s, t) and (ξ t ). In order to obtain the convergence of type (2.8) for the series X(s, t) defined by (2.10), it is sufficient to show that the array Ci (s, t) satisfy Condition 1.1. For this, the following lemma is needed.

E µ • C i (0) -1 (E d u ) = E [µ {x | C i (0)x > u}] = E 1I { C i (0)rθ >u} d(-r -α )P Θ (dθ) = u -α E C i (0)Θ α . Thus if (2.9) 0 < ∞ i=0 E C i (0)Θ α < ∞ the random
Lemma 2.2 Let A = (A 1 , . . . , A n ) be a block matrix of dimension nd × q with n blocks, A i be a d × q matrix with entries in R, i = 1, . . . , n. Then

(2.12) A ≤ n i=1 A i ,
where • denote the matrix norm.

Proof. Without loss of generality, we consider the vector norm (x 1 , . . . ,

x d ) d = max 1≤i≤d { x i } on R d . We may decompose the matrix A as (2.13) A = n i=1 B i
with B i is a nd × q matrix defined by

B i = (0, . . . , A i , . . . , 0)
where A i is in the ith position. We have, for all x ∈ R d

B i x nd = (0, . . . , A i x, . . . , 0) nd = A i x d , i = 1, . . . , n. Therefore (2.14) B i = sup x q ≤1 { B i x nd } = sup x q ≤1 { A i x d } = A i .
Inequality (2.12) follows from (2.13) and (2.14).

By Lemma 2.2 and the definition of the sequence Ci (s, t) (2.11), we obtain

(2.15) Ci (s, t) ≤ i j=(i-t+s)∨0 C j (t -i + j) .
Using (2.15) and the following inequality, for

a i > 0, i = 1, . . . , n, (2.16
)

n i=1 a i p ≤ c n i=1 a p i , where c = 1 if p ≤ 1, c = n p-1 if p > 1, we obtain, if α ∈ (0, 1) ∪ (1, 2), for some 0 < ε < α, ∞ i=0 E Ci (s, t) α-ε ≤ ∞ i=0 E   i j=(i-t+s)∨0 C j (t -i + j)   α-ε ≤ c 1 ∞ i=0 i j=(i-t+s)∨0 E C j (t -i + j) α-ε = c 1 t-s j=0 ∞ i=0 E C i (t -j) α-ε , (2.17)
where c 1 is a constant depending on α -ε and t -s. The latter inequality in combination with the condition (1.2) and the stationarity of the sequence

{C i (t)} implies ∞ i=0 E Ci (s, t) α-ε < ∞.
By the similar method for (2.17), we obtain, for some 0

< ε < α, ∞ i=0 E Ci (s, t) α+ε ≤ c 2 t-s j=0 ∞ i=0 E C i (t -j) α+ε , if α ∈ (0, 1) ∪ (1, 2), E ∞ i=0 Ci (s, t) α-ε α+ε α-ε ≤ c 3 t-s j=0 E ∞ i=0 C i (t -j) α-ε α+ε α-ε , if α ∈ {1, 2} and E ∞ i=0 Ci (s, t) 2 α+ε 2 ≤ c 4 t-s j=0 E ∞ i=0 C i (t -j) 2 α+ε 2 , if α ∈ (2, ∞),
where c 2 , c 3 and c 4 are the constants depending on α, ε and t -s. As a consequence, under the conditions of Theorem 2.1, it is not only the marginal distribution of (X t ) that is regularly varying, the same holds for all finite dimensional distributions.

Corollary 2.3 Under the conditions of Theorem 2.1, the stationary process (X t ) given by (1.1) satisfies, for s, t ∈ Z with s ≤ t

(2.18) P x -1 (X s , . . . , X t ) ∈ • P { ξ 0 > x} v -→ ∞ i=0 E µ • Ci (s, t) -1 (•) =: ν s,t (•), as x → ∞ on E (t-s+1)d
, where the sequence of random matrices Ci (s, t) is defined by (2.11). In particular, if X 0 is regularly varying, then the process (X t ) is jointly regularly varying, i.e. for all s, t ∈ Z withs ≤ t, random vector (X s , . . . , X t ) are regularly varying.

Tail process

The joint regular variation of the sequence (X t ) in particular means that there exists a process (Y t ) t∈Z in R d with P { Y 0 > y} = y -α for y ≥ 1 such that for all s, t ∈ Z with s ≤ t and as

x → ∞ L(x -1 (X s , . . . , X t ) | X 0 > x) L(Y s , . . . , Y t ),
see Theorem 2.1 in [START_REF] Basrak | Regularly varying multivariate time series[END_REF]. The process (Y t ) is called the tail process of (X t ). Moreover the so-called spectral process

Θ t = Yt Y 0 t∈Z is independent of Y 0 , see Theorem 3.1 in [1].
To understand the tail process in this case, it is helpful for instance to consider a set B in

E (t-s+1)d of the form B = E d bs × • • • × E d bt for arbitrary real constant b i > 0, i = s, .
. . , t. The limit of the probabilities, as x → ∞,

P x -1 (X s , . . . , X t ) ∈ B P { ξ 0 > x} is the sum ∞ i=0 E µ • Ci (s, t) -1 (B) = ∞ i=0 E µ x Ci (s, t)x ∈ B = ∞ i=t-s E [µ {x | C i-t+s (s)x > b s , . . . , C i (t)x > b t }] = ∞ i=t-s E P Θ (dθ) 1I { C i-t+s (s)rθ >bs} • • • 1I { C i (t)rθ >bt} d(-r -α ) = ∞ i=t-s E P Θ (dθ) max b s C i-t+s (s)θ , . . . , b t C i (t)θ -α = ∞ i=t-s E min C i-t+s (s)Θ α b α s , . . . , C i (t)Θ α b α t .
In particular, if s = t and b s = 1, one can apply this to obtain the following limit, as x → ∞,

(3.19) P { X 0 > x} P { ξ 0 > x} → ∞ i=0 E C i (0)Θ α .
The following result is the multivariate version of Breiman's lemma [START_REF] Breiman | On some limit theorems similar to the arc-sin law[END_REF] which appears as Proposition A. 1 in [START_REF] Basrak | Regular variation of GARCH processes[END_REF]. Lemma 3.1 Let Z be a q-dimensional random vector and let A be a d×q random matrix, independent of Z. Assume that Z is multivariate regularly varying of index α ∈ (0, ∞), i.e. there exists a non-null Radon measure µ on E q such that, as x → ∞,

P x -1 Z ∈ • P { Z > x} v -→ µ(•). If E A β < ∞ for some β > α, then in E d , as x → ∞, P x -1 AZ ∈ • P { Z > x} v -→ E µ • A -1 (•) .
Theorem 3.2 Let (X t ) t∈Z be a stationary process given by (1.1). Suppose that Condition 1.1 and 1.2 hold. If P (∩ i≥0 { C i (0)Θ = 0}) = 0, then for s, t ∈ Z with s ≤ 0 ≤ t and bounded and continuous function f

: (R d ) t-s+1 → R, E f X s X 0 , . . . , X t X 0 X 0 > x (3.20) → 1 ∞ i=0 E C i (0)Θ α ∞ i=0 E f C i+s (s)Θ C i (0)Θ , . . . , C i+t (t)Θ C i (0)Θ C i (0)Θ α , as x → ∞, where C i (•) = 0 if i < 0.
Proof. Let h : (R d ) t-s+1 → R be bounded and continuous. In view of convergence (3.19) and Corollary 2.3, as x → ∞,

E f X s x , . . . , X t x X 0 > x (3.21) → 1 ∞ i=0 E C i (0)Θ α h(y s , . . . , y t )1I { y 0 >1} ν s,t (dy). Note X(s, t) = (X s , . . . , X t ) = ∞ i=0 Ci (s, t)ξ t-i ,
where the sequence Ci (s, t) is defined by (2.11). By (2.15), (2.16), Condition 1.1 and the stationarity of sequence {C i (t)}, we have for some β > α,

E Ci (s, t) β ≤ c i j=(i-t+s)∨0 E C j (t -i + j) β = c i j=(i-t+s)∨0 E C j (0) β < ∞,
where c is a constant depending on β and t -s. Using Lemma 3.1, we obtain, for each i ≥ 0

(3.22) P x -1 Ci (s, t)ξ t-i ∈ • P { ξ 0 > x} v -→ E µ • Ci (s, t) -1 (•) =: ν (i) s,t (•), as x → ∞, in E (t-s+1)d
. By the definition of the measure ν s,t in (2.18),

(3.23) ν s,t (•) = ∞ i=0 ν (i) s,t (•).
We denote

H(x) = E h x -1 Ci (s, t)ξ t-i 1I { C i-t (0)ξ t-i >x} /P { ξ 0 > x} .
On the one hand, it follows from (3.22), as x → ∞,

H(x) = h(y s , . . . , y t )1I { y 0 >1} P x -1 Ci (s,t)ξ t-i (dy)/P { ξ 0 > x} → h(y s , . . . , y t )1I { y 0 >1} ν (i) s,t (dy). (3.24)
On the other hand, by independence between Ci (s, t) and ξ t-i and assumption A1, as x → ∞,

H(x) = E h Ci (s, t)y 1I { C i-t (0)y >1} P x -1 ξ t-i (dy)/P { ξ 0 > x} → ∞ 0 E h Ci (s, t)Θr 1I { C i-t (0)Θ r>1} d(-r -α ). (3.25) Since for i < t, ν (i) s,t {(y s , . . . , y t ) | y s = • • • = y 0 = 0 } = 1, (3.26) h(y s , . . . , y t )1I { y 0 >1} ν (i)
s,t (dy) = 0, i = 0, . . . , t -1.

In combination with (3.23), (3.24), (3.25) and (3.26), it follows that

h(y s , . . . , y t )1I { y 0 >1} ν s,t (dy) = ∞ i=t h(y s , . . . , y t )1I { y 0 >1} ν (i) s,t (dy) = ∞ i=t ∞ 0 E h Ci (s, t)Θr 1I { C i-t (0)Θ r>1} d(-r -α ). (3.27) Considering (3.21) and (3.27), we have, as x → ∞, E f X s x , . . . , X t x X 0 > x → I where I = 1 ∞ i=0 E C i (0)Θ α ∞ i=t ∞ 0 E h (C i-t+s (s)Θr, . . . , C i (t)Θr) 1I { C i-t (0)Θ r>1} d(-r -α ).
Applying this relation to the function h(x s , . . . , x t ) = f xs x 0 , . . . , xt x 0 to see that, as x → ∞, the left-hand side of (3.20) converges to

1 ∞ i=0 E C i (0)Θ α ∞ i=t ∞ 0 E f C i-t+s (s)Θ C i-t (0)Θ , . . . , C i (t)Θ C i-t (0)Θ 1I { C i-t (0)Θ r>1} d(-r -α ).
By Fubini's theorem, this is equal to the right-hand side of (3.20).

4 Applied models

X t = A t X t-1 + B t ξ t-1 + ξ t .
Iterating this equation backwards we arrive at the following MA(∞) representation of this process (4.29)

X t = ξ t + i≥1 (A t-i+1 + B t-i+1 ) i-2 j=0 A t-j ξ t-i ,
namely, the stationary solution can be represented as (1.1) with (4.30)

C i (t) =    Id, i = 0, (A t-i+1 + B t-i+1 ) i-2 j=0 A t-j , i ≥ 1.
From Theorem 3.1 in Hult and Samorodnitsky [START_REF] Hult | Tail probabilities for infinite series of regularly varying random vectors[END_REF] we obtain the following result.

Corollary 4.1 Suppose that ξ ∈ RV(α, µ) and there is some 0 < ε < α such that

(4.31) E A α+ε < 1, E B α+ε < ∞.
Then the series (4.29) converges a.s. and

(4.32) P{u -1 X 0 ∈ •} P{ ξ 0 > u} v -→ E µ • C j (0) -1 (•)
as u → ∞ on Rd \{0}.

Proof. Since C 0 (0) = Id, we have P( j≥0 { C j (0) = 0}) = 0. For α ∈ (0, 1) ∪ (1, 2), we have

E C i (0) α-ε = 1 + i≥1 E (A -i+1 + B -i+1 ) i-2 j=0 A -j α-ε ≤ 1 + i≥1 E A + B α-ε E i-2 j=0 A -j α-ε = 1 + E A + B α-ε i≥1 (E A α-ε ) i-1 = 1 + E A + B α-ε i≥0 (E A α-ε ) α+ε α-ε i α-ε
α+ε . (4.33) By Jensen's inequality, the last term is bounded by (4.34)

1 + E A + B α-ε i≥0 (E A α+ε ) i α-ε α+ε .
In the case of α ∈ (0, 1),

(4.35) E A + B α-ε < E A α-ε + E B α-ε .
In the case of α ∈ (1, 2)

(4.36) E A + B α-ε ≤ (E A + B α+ε ) α-ε α+ε ≤ E A α+ε 1 α+ε + E B α+ε 1 α+ε α-ε
.

Combining (4.33), (4.34), (4.35) and (4.36) proves

E C i (0) α-ε < ∞.
Similarly as in (4.33), we have (4.37)

E C i (0) α+ε ≤ 1 + E A + B α+ε i≥0 (E A α+ε ) i . If α + ε < 1, then (4.38) E A + B α+ε < E A α+ε + E B α+ε , else (4.39) E A + B α+ε ≤ (E A α+ε ) 1 α+ε + (E B α+ε ) 1 α+ε α+ε .
Combining (4.37), (4.38) and (4.39) proves

E C i (0) α+ε < ∞.
For α ∈ {1, 2}, using Lemma 3.2.1 in Kwapien and Woyczynski [START_REF] Kwapien | Random series and stochastic integrals: single and multiple[END_REF] it follows that

E C i (0) α-ε α+ε α-ε ≤ (E C i (0) α+ε ) α-ε α+ε α+ε α-ε =   1 + i≥1    E (A -i+1 + B -i+1 ) i-2 j=0 A -j α+ε    α-ε α+ε    α+ε α-ε ≤   1 + i≥1    E A + B α+ε E i-2 j=0 A -j α+ε    α-ε α+ε    α+ε α-ε =   1 + (E A + B α+ε ) α-ε α+ε i≥0 (E A α+ε ) i α-ε α+ε   α+ε α-ε . ( 4
E C i (0) α-ε α+ε α-ε < ∞.
For α ∈ (2, ∞), using Lemma 3.2.1 in Kwapien and Woyczynski [START_REF] Kwapien | Random series and stochastic integrals: single and multiple[END_REF] it follows that

E C i (0) 2 α+ε 2 ≤ (E C i (0) α+ε ) 2 α+ε α+ε 2 ≤   1 + (E A + B α+ε ) 2 α+ε i≥0 (E A α+ε ) i 2 α+ε   α+ε 2 . ( 4
E C i (0) 2 α+ε 2 < ∞.
Here is another proof of Corollary 4.1.

Proof. We denote (4.42) Xt = (X t , ξ t , ξ t+1 ) ,

Ãt =   A t B t Id 0 0 Id 0 0 0   and Z t = (0, 0, ξ t+1 ) .
Then the relation (4.28) can be wrote as

(4.43) Xt = Ãt Xt-1 + Z t .
Iterating this equation we get

Xt = ∞ i=0 C i (t)Z t-i
where C 0 (t) = Id, C 1 (t) = Ãt and

C i (t) =     i-1 j=0 A t-j B t-i+1 i-2 j=0 A t-j (A t-i+2 + B t-i+2 ) i-3 j=0 A t-j 0 0 0 0 0 0     , i ≥ 2.
It

≤ 1 + c i≥1   E i-1 j=0 A t-j α-ε + E   B t-i+1 i-2 j=0 A t-j   α-ε +E   A t-i+2 + B t-i+2 i-3 j=0 A t-j   α-ε   = 1 + c   i≥1 E A α-ε i + E B α-ε i≥0 E A α-ε i +E A + B α-ε i≥0 E A α-ε i + 2 α-ε   ,
where c is a constant depending on α -ε. By Jensen's inequality, Similarly we can prove i≥0 E C i (0) α+ε < ∞.

For α ∈ {1, 2}, using Lemma 3.2.1 in Kwapien and Woyczynski [START_REF] Kwapien | Random series and stochastic integrals: single and multiple[END_REF] it follows that

  .40) Considering (4.39) and (4.40) we get

( 4 .

 4 46) E A α-ε ≤ E A α+ε α-ε α+ε .Combining (4.45), (4.46), (4.38) and (4.39) proves i≥0 E C i (0) α-ε < ∞.

  16) and (4.44), (E C i (0) α+ε ) α-ε α+ε is bounded by the right-hand side of (4.45) which is bounded. Hence (4.48)E C i (0) α-ε α+ε α-ε < ∞.

  Assume that ((A t , B t , ξ t )) t is an i.i.d. sequence of random vectors in R d×d × R d×d × R d , for some d ≥ 1. Throughout we assume that ξ t 's are regularly varying random vectors. Suppose that a stationary sequence (X t ) t∈Z with value in R d satisfies a multivariate random coefficient ARMA(1,1) equation of the following form(4.28)

	4.1 Heavy tailed multivariate ARMA(1,1) process

  is equivalent to prove the convergence (4.32) of vector Xt defined by (4.42). By Lemma 2.2 and the properties of matrix norm, for i ≥ 1where A i = B i = Id if i > t. It follows from (2.16), if α ∈ (0, 1) ∪ (2, 1)

		i-1	i-2	i-3
	(4.44)	C i (t) ≤	A t-j + B t-i+1	A t-j + A t-i+2 + B t-i+2	A t-j
		j=0	j=0	j=0
	(4.45)	E C i (0) α-ε
		i≥0

 [START_REF] Breiman | On some limit theorems similar to the arc-sin law[END_REF]Research supported by IAP research network grant nr. P6/03 of the Belgian government (Belgian Science Policy).

By similar method of (4.48), we get