text:
# Untitled Document

On Reducing the Complexity of Matrix Clocks

Lúcia M. A. Drummond

Universidade Federal Fluminense

Instituto de Computação

Rua Passo da Pátria, 156

24210-240 Niterói - RJ, Brazil

lucia@dcc.ic.uff.br

Valmir C. Barbosa

Universidade Federal do Rio de Janeiro

Programa de Engenharia de Sistemas e Computação, COPPE

Caixa Postal 68511

21941-972 Rio de Janeiro - RJ, Brazil

valmir@cos.ufrj.br

Abstract

Matrix clocks are a generalization of the notion of vector clocks that allows the local representation of causal precedence to reach into an asynchronous distributed computation’s past with depth $`x`$, where $`x1`$ is an integer. Maintaining matrix clocks correctly in a system of $`n`$ nodes requires that every message be accompanied by $`O(n^x)`$ numbers, which reflects an exponential dependency of the complexity of matrix clocks upon the desired depth $`x`$. We introduce a novel type of matrix clock, one that requires only $`nx`$ numbers to be attached to each message while maintaining what for many applications may be the most significant portion of the information that the original matrix clock carries. In order to illustrate the new clock’s applicability, we demonstrate its use in the monitoring of certain resource-sharing computations.

Keywords: Vector clocks, matrix clocks.

Address for correspondence and proofs:

Valmir C. Barbosa, valmir@cos.ufrj.br

Programa de Sistemas, COPPE/UFRJ

Caixa Postal 68511

21941-972 Rio de Janeiro - RJ

Brazil

1. Introduction and background

We consider an undirected graph $`G`$ on $`n`$ nodes. Each node in $`G`$ stands for a computational process and undirected edges in $`G`$ represent the possibilities for bidirectional point-to-point communication between pairs of processes. A fully asynchronous distributed computation carried out by the distributed system represented by $`G`$ can be viewed as a set of events occurring at the various nodes. An event is the sending of a message by a node to another node to which it is directly connected in $`G`$ (a neighbor in $`G`$), or the reception of a message from such a neighbor, or yet the occurrence at a node of any internal state transition of relevance (“state” and “relevance” here are highly dependent upon the particular computation at hand, and are left unspecified).

The standard framework for analyzing such a system is the partial order, often denoted by $``$, that formalizes the usual “happened-before” notion of distributed computing . This partial order is the transitive closure of the more elementary relation to which the ordered pair $`(v,v^{})`$ of events belongs if either $`v`$ and $`v^{}`$ are consecutive events at a same node in $`G`$ or $`v`$ and $`v^{}`$ are, respectively, the sending and receiving of a message between neighbors in $`G`$.

At node $`i`$, we let local time be assessed by the number $`t_i`$ of events that already occurred. We interchangeably adopt the notations $`t_i=\mathrm{𝑡𝑖𝑚𝑒}_i(v)`$ and $`v=\mathrm{𝑒𝑣𝑒𝑛𝑡}_i(t_i)`$ to indicate that $`v`$ is the $`t_i`$th event occurring at node $`i`$, provided $`t_i1`$. Another important partial order on the set of events is the relation that gives the predecessor at node $`ji`$ of an event $`v^{}`$ occurring at node $`i`$. We say that an event $`v`$ is such a predecessor, denoted by $`v=\mathrm{𝑝𝑟𝑒𝑑}_j(v^{})`$, if $`v`$ is the event occurring at node $`j`$ such that $`vv^{}`$ for which $`\mathrm{𝑡𝑖𝑚𝑒}_j(v)`$ is greatest. If no such $`v`$ exists, then $`\mathrm{𝑝𝑟𝑒𝑑}_j(v^{})`$ is undefined and $`\mathrm{𝑡𝑖𝑚𝑒}_j\left(\mathrm{𝑝𝑟𝑒𝑑}_j(v^{})\right)`$ is assumed to be zero.

The relation $`\mathrm{𝑝𝑟𝑒𝑑}_j`$ allows the definition of vector clocks \[11–14, 17\], as follows. The vector clock of node $`i`$ at time $`t_i`$ (that is, following the occurrence of $`t_i`$ events at node $`i`$), denoted by $`V^i(t_i)`$, is a vector whose $`j`$th component, for $`1jn`$, is either equal to $`0`$, if $`t_i=0`$, or given by

$$V_j^i(t_i)=\{\begin{array}{cc}t_i,\hfill & \text{if }j=i\text{;}\hfill \\ \mathrm{𝑡𝑖𝑚𝑒}_j\left(\mathrm{𝑝𝑟𝑒𝑑}_j\left(\mathrm{𝑒𝑣𝑒𝑛𝑡}_i(t_i)\right)\right),\hfill & \text{if }ji\text{,}\hfill \end{array}$$
$`(1)`$
if $`t_i1`$. In other words, $`V_j^i(t_i)`$ is either the current time at node $`i`$, if $`j=i`$, or is the time at node $`j`$ that results from the occurrence at that node of the predecessor of the $`t_i`$th event of node $`i`$, otherwise. If no such event exists (i.e., $`t_i=0`$), then $`V_j^i(t_i)=0`$.

Vector clocks evolve following two simple rules:

$``$Upon sending a message to one of its neighbors, node $`i`$ attaches $`V^i(t_i)`$ to the message, where $`t_i`$ is assumed to already account for the message that is being sent.

$``$Upon receiving a message from node $`k`$ with attached vector clock $`V^k`$, node $`i`$ sets $`V_i^i(t_i)`$ to $`t_i`$ (which is assumed to already reflect the reception of the message) and $`V_j^i(t_i)`$ to $`\mathrm{max}\{V_j^i(t_i1),V_j^k\}`$, for $`ji`$.

It is a simple matter to verify that these rules do indeed maintain vector clocks consistently with their definition \[11–14, 17\]. Under these rules or variations thereof, vector clocks have proven useful in a variety of distributed algorithms to detect some types of global predicates .

For large $`n`$, attaching a vector clock to every message that is sent is likely to become burdensome, so the question arises whether less costly implementations are possible. Under the very general assumptions we have made concerning the nature of $`G`$ as a distributed system, the answer is negative: a result similar to Dilworth’s theorem on partially ordered sets  establishes that the size-$`n`$ attachments are necessary . However, it is possible to use more economical attachments if the edges of $`G`$ provide FIFO communication , or if certain aspects of the structure of $`G`$ can be taken into account , or yet if the full capabilities of vector clocks are not needed .

One natural generalization of the notion of vector clocks is the notion of matrix clocks . For an integer $`x1`$, the $`x`$-dimensional matrix clock of node $`i`$ at time $`t_i`$, denoted by $`M^i(t_i)`$, has $`O(n^x)`$ components. For $`1j_1,\mathrm{},j_xn`$, component $`M_{j_1,\mathrm{},j_x}^i(t_i)`$ of $`M^i(t_i)`$ is only defined for $`i=j_1=\mathrm{}=j_x`$ and for $`ij_1\mathrm{}j_x`$. As in the definition of $`V^i(t_i)`$, $`M_{j_1,\mathrm{},j_x}^i(t_i)=0`$ if $`t_i=0`$. For $`t_i1`$, on the other hand, we have

$$M_{j_1,\mathrm{},j_x}^i(t_i)=\{\begin{array}{cc}t_i,\hfill & \text{if }i=j_1=\mathrm{}=j_x\text{;}\hfill \\ \mathrm{𝑡𝑖𝑚𝑒}_{j_x}\left(\mathrm{𝑝𝑟𝑒𝑑}_{j_x}\mathrm{}\mathrm{𝑝𝑟𝑒𝑑}_{j_1}\left(\mathrm{𝑒𝑣𝑒𝑛𝑡}_i(t_i)\right)\right),\hfill & \text{if }ij_1\mathrm{}j_x\text{,}\hfill \end{array}$$
$`(2)`$
which, for $`ij_1\mathrm{}j_x`$, first takes the predecessor at node $`j_1`$ of the $`t_i`$th event occurring at node $`i`$, then the predecessor at node $`j_2`$ of that predecessor, and so on through node $`j_x`$, whose local time after the occurrence of the last predecessor in the chain is assigned to $`M_{j_1,\mathrm{},j_x}^i(t_i)`$. Should any of these predecessors be undefined, the ones that follow it in the remaining nodes are undefined as well, and the local time that results at the end is zero. It is straightforward to see that, for $`x=1`$, this definition is equivalent to the definition of a vector clock in (1). Similarly, the maintenance of matrix clocks follows rules entirely analogous to those used to maintain vector clocks .

While the $`j`$th component of the vector clock following the occurrence of event $`v^{}`$ at node $`ij`$ gives the time resulting at node $`j`$ from the occurrence of $`\mathrm{𝑝𝑟𝑒𝑑}_j(v^{})`$, the analogous interpretation that exists for matrix clocks requires the introduction of additional notation. Specifically, the definition of a set of events encompassing the possible $`x`$-fold compositions of the relation $`\mathrm{𝑝𝑟𝑒𝑑}_j`$ with itself, denoted by $`\mathrm{𝑃𝑟𝑒𝑑}_j^{(x)}`$, is necessary. If an event $`v`$ occurs at node $`j`$, then we say that $`v\mathrm{𝑃𝑟𝑒𝑑}_j^{(x)}(v^{})`$ for an event $`v^{}`$ occurring at node $`i`$ if one of the following holds:

$``$$`x=1`$, $`ji`$, and $`v=\mathrm{𝑝𝑟𝑒𝑑}_j(v^{})`$.

$``$$`x>1`$ and there exists $`ki`$ such that an event $`\overline{v}`$ occurs at node $`k`$ for which $`\overline{v}=\mathrm{𝑝𝑟𝑒𝑑}_k(v^{})`$ and $`v\mathrm{𝑃𝑟𝑒𝑑}_j^{(x1)}(\overline{v})`$.

Note that this definition requires $`jk`$, in addition to $`ki`$, for $`v\mathrm{𝑃𝑟𝑒𝑑}_j^{(x)}(v^{})`$ to hold when $`x=2`$.

For $`x>1`$ and $`ij_1\mathrm{}j_x=j`$, this definition allows the following interpretation of entry $`j_1,\mathrm{},j_x`$ of the matrix clock that follows the occurrence of event $`v^{}`$ at node $`i`$, that is, of $`M_{j_1,\mathrm{},j_{x1},j}^i\left(\mathrm{𝑡𝑖𝑚𝑒}_i(v^{})\right)`$. It gives the time resulting at node $`j`$ from the occurrence of an event that is in the set $`\mathrm{𝑃𝑟𝑒𝑑}_j^{(x)}(v^{})`$, so long as it is nonempty. Of course, the number of possibilities for such an event is $`O(n^{x1})`$ in the worst case, owing to the several possible combinations of $`j_1,\mathrm{},j_{x1}`$.

Interestingly, applications that require the full capabilities of matrix clocks have yet to be identified. In fact, it seems a simple matter to argue that a slightly more sophisticated use of vector clocks or simply the use of two-dimensional matrix clocks (the $`x=2`$ case) suffices to tackle some of the problems that have been offered as possible applications of higher-dimensional matrix clocks . What we do in this paper is to demonstrate how the distributed monitoring of certain resource-sharing computations can benefit from the use of matrix clocks and that, at least for such computations, it is possible to employ much less complex matrix clocks (that is, matrix clocks with many fewer components), which nonetheless retain the ability to reach into the computation’s past with arbitrary depth.

The key to this reduction in complexity is the use of one single event in place of each set $`\mathrm{𝑃𝑟𝑒𝑑}_j^{(y)}(v^{})`$, for $`1yx`$. In other words, for each node $`j`$ the matrix clock we introduce retains only one of the $`O(n^{y1})`$ components of each of the $`x`$ $`y`$-dimensional original matrix clocks—one component from the vector clock, one from the two-dimensional matrix clock, and so on through one component of the $`x`$-dimensional matrix clock. As we argue in Section 3, following a brief discussion of the resource-sharing context in Section 2, this simplification leads to matrix clocks of size $`nx`$, therefore considerably less complex than the $`O(n^x)`$-sized original matrix clocks.

The single other attempt at reducing the complexity of a matrix clock that we know of was given for the $`x=2`$ case specifically and culminated with the introduction of two techniques . The first one requires attachments of size $`O(n)`$ (a considerable reduction from the original $`O(n^2)`$), but is only applicable if the full asynchronism we have been assuming does not hold; it is therefore of little interest in our context. The other technique is somewhat closer to our own in spirit, since it aims at approximating the two-dimensional matrix clock by retaining only $`k`$ of the $`O(n)`$ components that correspond to each node. However, the criterion to select the components to be retained is to choose the $`k`$ greatest components, which seems unrelated to our own criterion, based as it is on the $`\mathrm{𝑃𝑟𝑒𝑑}_j`$ sets.

We give concluding remarks in Section 4, after a discussion of how the technique of Section 3 can successfully solve the problem posed in Section 2.

2. Monitoring resource-sharing computations

The resource-sharing computation we consider is one of the classical solutions to the paradigmatic Dining Philosophers Problem (DPP)  in generalized form . In this case, $`G`$’s edge set is constructed from a given set of resources and from subsets of that set, one for each node, indicating which resources can be ever needed by that node. This construction places an edge between nodes $`i`$ and $`j`$ if the sets of resources ever to be needed by $`i`$ and $`j`$ have a nonempty intersection. Notice that this construction is consonant with the interpretation of edges as bidirectional communication channels, because it deploys edges between every pair of nodes that may ever compete for a same resource and must therefore be able to communicate with each other to resolve conflicts.

In DPP, the computation carried out by a node makes it cycle endlessly through three states, which are identified with the conditions of being idle, being in the process of acquiring exclusive access to the resources it needs, and using those resources for a finite period of time. While in the idle state, the node starts acquiring exclusive access to resources when the need arises to compute on shared resources. It is a requirement of DPP that the node must acquire exclusive access to all the resources it shares with all its neighbors, so it suffices for the node to acquire a token object it shares with each of its neighbors (the “fork,” as it is called), each object representing all the resources it shares with that particular neighbor. When in possession of all forks, the node may start using the shared resources .

The process of collecting forks from neighbors follows a protocol based on the sending of request messages by the node that needs the forks and the sending of the forks themselves by the nodes that have them. More than one protocol exists, each implementing a different rule to ensure the absence of deadlocks and lockouts during the computation. The solution we consider in this section is based on relative priorities assigned to nodes. Another prominent solution is also based on the assignment of relative priorities, but to the resources instead of to the nodes .

The priority scheme of interest to us is based on the graph-theoretic concept of an acyclic orientation of $`G`$, which is an assignment of directions to the edges of $`G`$ in such a way that directed cycles are not formed. Such an acyclic orientation is then a partial order of the nodes of $`G`$, and is as such suitable to indicate, given a pair of neighbors, which of the two has priority over the other. Most of the details of this priority scheme are not relevant to our present discussion, but before stating its essentials we do mention that the lockout-freedom requirement leads the acyclic orientation to be changed dynamically (so that relative priorities are never fixed), which in turn leads to a rich dynamics in the set of all the acyclic orientations of $`G`$ and to important concurrency-related issues \[3–5\].

Once a priority scheme is available over the set of nodes, what matters to us is how it is used in the fork-collecting protocol. When a request for fork arrives at node $`j`$ from node $`i`$, $`j`$ sends $`i`$ the fork they share if $`j`$ is either idle or is also collecting forks but does not have priority over $`i`$. If $`j`$ is also collecting forks and has priority over $`i`$, or if $`j`$ is using shared resources, then the sending of the fork to $`i`$ is postponed to until $`j`$ has finished using the shared resources. Note that two types of wait may happen here. If $`j`$ is using shared resources when the request arrives, then the wait is independent of $`n`$. If $`j`$ is also collecting forks, then the wait for $`j`$ to start using the shared resources and ultimately to send $`i`$ the fork is in the worst case $`n1`$ . The reason for this is simple: $`j`$ is waiting for a fork from another node, which may in turn be waiting for a fork from yet another node, and so on. Because the priority scheme is built on acyclic orientations of $`G`$, such a chain of waits does necessarily end and is $`n1`$ nodes long in the worst case.

Whether such long waits occur or not is of course dependent upon the details of each execution of the resource-sharing computation. But if they do occur, one possibility for reducing the average wait is to increase the availability of certain critical resources so that $`G`$ becomes less dense . Perhaps another possibility would be to fine-tune some of the characteristics of each node’s participation in the overall computation, such as the duration of its idle period, which could be subject to a mandatory lower bound, for example. In any event, the ability to locally detect long waits (a global property, since it relates to the notion of time in fully asynchronous distributed systems ) and identify the nodes at which the wait chains end is crucial.

To see how this relates to the formalism of Section 1, suppose our resource-sharing computation consists of the exchange of fork-bearing messages only (that is, request messages and all other messages involved in the computation, such as those used to update acyclic orientations, are ignored). For distinct nodes $`i`$ and $`j`$, and for an event $`v^{}`$ occurring at node $`i`$ at the reception of a fork, the set $`\mathrm{𝑃𝑟𝑒𝑑}_j^{(x)}(v^{})`$ for $`x1`$ is either empty or only contains events that correspond to the sending of forks by node $`j`$. Now suppose we systematically investigate the sets $`\mathrm{𝑃𝑟𝑒𝑑}_j^{(x)}(v^{})`$ for every appropriate $`j`$ by letting $`x`$ increase from $`1`$ through $`n1`$. Suppose also that, for each $`j`$, we record the first $`x`$ that is found such that $`\mathrm{𝑃𝑟𝑒𝑑}_j^{(x)}(v^{})`$ contains the sending of a fork as response to the reception of a request message without any wait for fork collection on the part of $`j`$. The greatest such value to be recorded, say $`x^{}`$, has special significance: it means that the eventual reception of a fork by node $`i`$ through $`v^{}`$ is the result of an $`x^{}`$-long chain of waits.

If this were the only information of interest, then Lamport’s clocks  could be used trivially to discover it. However, taking corrective measures may require a wider array of long wait chains to be found (not simply the longest), as well as the nodes at which those chains end. The matrix clock that we introduce in Section 3 is capable of conveying this information to node $`i`$ succinctly, so long as there exists enough redundancy in the sets $`\mathrm{𝑃𝑟𝑒𝑑}_j^{(x)}(v^{})`$ that attaching only $`nX`$ integers to forks suffices, where $`X`$ is a threshold in the interval $`[1,n1]`$ beyond which wait chains are known not to occur, given the structure of $`G`$ and the initial arrangement of priorities.1 Discovering the value of $`X`$ is not necessarily a simple task, but some empirical knowledge has already been accumulated for modestly-sized systems ; in any case, in the likely event that $`X`$ cannot be determined with certainty, there is always the possibility of adaptation as the resource-sharing computation is observed on the system at hand. It so happens that such redundancy clearly exists: the only events in $`\mathrm{𝑃𝑟𝑒𝑑}_j^{(x)}(v^{})`$ that matter are those corresponding to the sending of forks without any wait for fork collection. Detecting any one of them suffices, so we may as well settle for the latest, that is, one single event from the whole set. We return to this in Section 4.

3. A simpler matrix clock

For $`x1`$, the new matrix clock we introduce is an $`x\times n`$ matrix. For node $`i`$ at time $`t_i`$ (i.e., following the occurrence of $`t_i`$ events at node $`i`$), it is denoted by $`C^i(t_i)`$. For $`1yx`$ and $`1jn`$, component $`C_{y,j}^i(t_i)`$ of $`C^i(t_i)`$ is defined as follows. If $`t_i=0`$, then $`C_{y,j}^i(t_i)=0`$, as for vector clocks and the matrix clocks of Section 1. If $`t_i1`$, then we have

$$C_{y,j}^i(t_i)=\{\begin{array}{cc}t_i,\hfill & \text{if }y=1\text{ and }j=i\text{;}\hfill \\ \mathrm{max}_{ij_1\mathrm{}j_y=j}\left\{\mathrm{𝑡𝑖𝑚𝑒}_{j_y}\left(\mathrm{𝑝𝑟𝑒𝑑}_{j_y}\mathrm{}\mathrm{𝑝𝑟𝑒𝑑}_{j_1}\left(\mathrm{𝑒𝑣𝑒𝑛𝑡}_i(t_i)\right)\right)\right\},\hfill & \text{if }y>1\text{ or }ji\text{.}\hfill \end{array}$$
$`(3)`$
Note, first of all, that for $`y=1`$ this definition is equivalent to the definition of a vector clock in (1). Thus, the first row of $`C^i(t_i)`$ is the vector clock $`V^i(t_i)`$. For $`y>1`$, the definition in (3) implies, according to the interpretation of matrix clocks that follows our definition in (2), that

$$C_{y,j}^i(t_i)=\{\begin{array}{cc}\mathrm{max}\left\{\mathrm{𝑡𝑖𝑚𝑒}_j(v)\right|v\mathrm{𝑃𝑟𝑒𝑑}_j^{(y)}(v^{})\},\hfill & \text{if }\mathrm{𝑃𝑟𝑒𝑑}_j^{(y)}(v^{})\mathrm{}\text{;}\hfill \\ 0,\hfill & \text{if }\mathrm{𝑃𝑟𝑒𝑑}_j^{(y)}(v^{})=\mathrm{}\text{,}\hfill \end{array}$$
$`(4)`$
where $`v^{}`$ is the $`t_i`$th event occurring at node $`i`$. What this means is that, of all the $`O(n^{y1})`$ events that may exist in $`\mathrm{𝑃𝑟𝑒𝑑}_j^{(y)}(v^{})`$, only one (the one to have occurred latest at node $`j`$) makes it to the matrix clock $`C^i(t_i)`$. Note also that (3) implies the equality in (4) for $`y=1`$ as well, so long as $`ji`$. In this case, $`\mathrm{𝑃𝑟𝑒𝑑}_j^{(y)}(v^{})`$, if nonempty, is the singleton $`\left\{\mathrm{𝑝𝑟𝑒𝑑}_j(v^{})\right\}`$.

Before we derive the update rules for our new matrix clock, let us pause and examine an example. Consider Figure 1, where a computation on six nodes is illustrated by means of the usual message diagram that forms the basis of the relation $``$. Nodes are numbered $`1`$ through $`6`$, and in the figure local time elapses from left to right independently for each node. Filled circles represent events and the arrows connecting two events indicate messages.

Figure 1. A computation fragment on six nodes.

Three events are singled out in Figure 1, namely $`v^{}`$, $`v_1`$, and $`v_2`$. They are related to one another in such a way that $`\mathrm{𝑃𝑟𝑒𝑑}_6^{(4)}(v^{})=\{v_1,v_2\}`$, that is, $`v_1`$ and $`v_2`$ are the “depth-$`4`$” predecessors of $`v^{}`$ at node $`6`$. Recalling that $`\mathrm{𝑡𝑖𝑚𝑒}_1(v^{})=3`$, there are two components of the four-dimensional matrix clock $`M^1(3)`$ that reflect this relationship among the three events, namely

$$M_{3,4,5,6}^1(3)=\mathrm{𝑡𝑖𝑚𝑒}_6(v_1)=2$$

and

$$M_{2,3,4,6}^1(3)=\mathrm{𝑡𝑖𝑚𝑒}_6(v_2)=3.$$

These follow directly from (2) with $`x=4`$. By (3) with $`y=4`$, the same diagram of Figure 1 yields

$$\begin{array}{cc}\hfill C_{4,6}^1(3)& =\mathrm{max}\{\mathrm{𝑡𝑖𝑚𝑒}_6\left(\mathrm{𝑝𝑟𝑒𝑑}_6\mathrm{𝑝𝑟𝑒𝑑}_5\mathrm{𝑝𝑟𝑒𝑑}_4\mathrm{𝑝𝑟𝑒𝑑}_3(v^{})\right),\mathrm{𝑡𝑖𝑚𝑒}_6\left(\mathrm{𝑝𝑟𝑒𝑑}_6\mathrm{𝑝𝑟𝑒𝑑}_4\mathrm{𝑝𝑟𝑒𝑑}_3\mathrm{𝑝𝑟𝑒𝑑}_2(v^{})\right)\}\hfill \\ & =\mathrm{max}\{\mathrm{𝑡𝑖𝑚𝑒}_6(v_1),\mathrm{𝑡𝑖𝑚𝑒}_6(v_2)\}\hfill \\ & =3,\hfill \end{array}$$

which is clearly also in consonance with (4).

Let us now look at the update rules. For $`y=1`$, these rules are the same as those given for vector clocks in Section 1. For $`y>1`$, we know from the definition of $`\mathrm{𝑃𝑟𝑒𝑑}_j^{(y)}`$ that $`v\mathrm{𝑃𝑟𝑒𝑑}_j^{(y)}(v^{})`$ if and only if there exists $`ki`$ such that an event $`\overline{v}`$ occurs at node $`k`$ for which $`\overline{v}=\mathrm{𝑝𝑟𝑒𝑑}_k(v^{})`$ and $`v\mathrm{𝑃𝑟𝑒𝑑}_j^{(y1)}(\overline{v})`$. Thus, if we let $`t_k=\mathrm{𝑡𝑖𝑚𝑒}_k(\overline{v})`$ and consider the matrix clocks resulting from the occurrence of $`\overline{v}`$ and of $`v^{}`$ ($`C^k(t_k)`$ and $`C^i(t_i)`$, respectively), then it follows from (4) that

$$C_{y,j}^i(t_i)=\{\begin{array}{cc}\mathrm{max}_{kK^i(t_i,y,j)}C_{y1,j}^k(t_k),\hfill & \text{if }K^i(t_i,y,j)\mathrm{}\text{;}\hfill \\ 0,\hfill & \text{if }K^i(t_i,y,j)=\mathrm{}\text{,}\hfill \end{array}$$
$`(5)`$
where $`K^i(t_i,y,j)`$ is the set comprising every appropriate $`k`$. Notice, in (5), that $`k=j`$ can never occur as the maximum is taken if $`y=2`$: aside from $`i`$, node $`j`$ is the only node that cannot possibly be a member of $`K^i(t_i,2,j)`$, by (3).

According to (5), and recalling once again the special case of $`y=1`$, we are then left with the following two rules for the evolution of our new matrix clocks:

$``$Upon sending a message to one of its neighbors, node $`i`$ attaches $`C^i(t_i)`$ to the message, where $`t_i`$ is assumed to already account for the message that is being sent.

$``$Upon receiving a message from node $`k`$ with attached matrix clock $`C^k`$, and assuming that $`t_i`$ already reflects the reception of the message, node $`i`$ sets $`C_{y,j}^i(t_i)`$ to

$$\{\begin{array}{cc}t_i,\hfill & \text{if }y=1\text{ and }j=i\text{;}\hfill \\ \mathrm{max}\{C_{y,j}^i(t_i1),C_{y,j}^k\},\hfill & \text{if }y=1\text{ and }ji\text{;}\hfill \\ C_{y,j}^i(t_i1),\hfill & \text{if }y=2\text{ and }j=k\text{;}\hfill \\ \mathrm{max}\{C_{y,j}^i(t_i1),C_{y1,j}^k\},\hfill & \text{if }1<yx\text{, provided }y>2\text{ or }jk\text{.}\hfill \end{array}$$

According to these rules, every message carries an attachment of $`nx`$ integers.

4. Discussion and concluding remarks

We are now in a position to return to the problem posed in Section 2, namely the problem of monitoring executions of the solution to DPP that employs a partial order on $`G`$’s set of nodes to establish priorities. As we discussed in that section, the overall goal is to allow nodes to detect locally, upon receiving a fork, whether the delivery of that fork is the result of a chain of fork deliveries that started too far back in the past. In the affirmative case, the wait since the fork was requested will have been too long, in terms of the usual notions of time complexity in asynchronous distributed computations.

More specifically, if $`v^{}`$ is the event corresponding to the reception of a fork at node $`i`$, then the goal is for $`i`$ to be able to detect the existence of an event in $`\mathrm{𝑃𝑟𝑒𝑑}_j^{(y)}(v^{})`$ that corresponds to the sending of a fork by node $`j`$ either immediately upon the reception of the request for that fork or, if node $`j`$ was using shared resources when the request arrived, immediately upon finishing. Here $`1yX`$ and $`j`$ is any node, provided $`y=1`$ and $`j=i`$ do not occur in conjunction. The value of $`X`$ is such that $`1Xn1`$, and is chosen as a bound to reflect the maximum possible chain of waits. As a consequence, the sets $`\mathrm{𝑃𝑟𝑒𝑑}_j^{(y)}(v^{})`$ must include fork-related events only.

The new matrix clocks introduced in Section 3 can be used for this detection with only minimal adaptation. The key ingredients are:

$``$The only messages sent by node $`i`$ to be tagged with the matrix clock $`C^i`$ are forks. If the sending of a fork by node $`i`$ does not depend on further fork collection by $`i`$, then every component of $`C^i`$ other than $`C_{1,i}^i`$ is reset to zero before it is attached to the fork. Matrix clocks are $`X\times n`$, and all further handling of them follows the rules given in Section 3.

$``$Upon receiving a fork with attached matrix clock, and having updated $`C^i`$ accordingly, node $`i`$ looks for components of $`C^i`$ that contain nonzero values. If $`C_{y,j}^i`$ is one such component for $`y>1`$ or $`ji`$, then a wait chain of length $`y`$ that ends at node $`j`$ has been discovered and can be checked against a certain threshold $`X^{}<X`$ representing the maximum allowable chain length. If $`C^i`$ contains zero components at all positions but $`(1,i)`$, then it is certain that the value of $`X`$ has to be revised, as clearly a wait chain exists whose length surpasses $`X`$. A greater value for $`X`$ is then needed.

This strategy reflects the general method of Section 3 when applied to events that relate to the flow of forks only. Whenever the request for a fork is received and the fork can be sent without the need for any forks to be received by the node in question, say node $`j`$, zeroes get inserted into the matrix clock at all positions but $`(1,j)`$ and are sent along with the fork. The sending of this fork may unleash the sending of forks by other nodes in a chain of events, and along the way the original value of $`C_{1,j}^j`$ may never be reset to zero, reflecting the increasing length of the wait chain rooted at node $`j`$. The reception of a fork whose matrix clock has such a nonzero component beyond row $`X^{}`$ is then an indication that such reception is part of chains whose lengths are considered too long.

It is now worth returning to Figure 1 in order to view its diagram as representing the fork-bearing messages in a fragment of a DPP computation of the type we have been considering. For such, consider Figure 2, where two six-node graphs are depicted. Figure 2(a) shows the graph $`G`$ corresponding to a certain resource-sharing pattern for the six nodes; it shows an acyclic orientation of the edges of $`G`$ as well, indicating the initial arrangement of priorities. In a situation of heavy demand for resources by all six nodes, only node $`6`$ can acquire all the (three) forks it needs and proceed. All others must wait along the length-$`5`$ wait chain shown in Figure 2(b): node $`1`$ for node $`2`$, $`2`$ for $`3`$, and so on through node $`6`$.

Figure 2. $`G`$ oriented acyclically (a) and a possible wait chain (b).

Eventually, it must happen that node $`6`$ sends out its three forks. If we assume that the corresponding request messages arrive when node $`6`$ already holds all three forks, then upon sending them out all components of $`C^6`$ are sent as zeroes, except for $`C_{1,6}^6`$, sent as $`1`$, $`2`$, or $`3`$, depending on the destination. For $`X=5`$, at the occurrence of $`v^{}`$ node $`1`$ updates its matrix clock in such a way that $`C_{5,6}^1=2`$, which characterizes the length-$`5`$ wait chain ending at node $`6`$.

In summary, we have in this paper introduced a novel notion of matrix clocks. Similarly to the original matrix clocks , whose definition appears in (2), our matrix clock has the potential of reflecting causal dependencies in the flow of messages that stretch as far as depth $`x`$ into the past. Unlike those original matrix clocks, however, ours increases with $`x`$ as $`nx`$, while in the original case the growth is according to an $`O(n^x)`$ function, that is, exponentially.

We have illustrated the applicability of the new matrix clocks with an example from the area of resource-sharing problems. What we have demonstrated is a means of collecting information locally during the resource-sharing computation so that exceedingly long global waits can be detected, possibly indicating the need for overall system re-structuring, as described in Section 2.

Acknowledgments

The authors acknowledge partial support from CNPq, CAPES, the PRONEX initiative of Brazil’s MCT under contract 41.96.0857.00, and a FAPERJ BBP grant.

References

1. M. Ahuja, T. Carlson, and A. Gahlot, “Passive-space and time view: vector clocks for achieving higher performance, program correction, and distributed computing,” IEEE Trans. on Software Engineering 19 (1993), 845–855.

2. V. C. Barbosa, An Introduction to Distributed Algorithms, The MIT Press, Cambridge, MA, 1996.

3. V. C. Barbosa, An Atlas of Edge-Reversal Dynamics, Chapman & Hall/CRC, London, UK, 2000.

4. V. C. Barbosa, “The combinatorics of resource sharing,” in R. Corrêa et alii (Eds.), Models for Parallel and Distributed Computation: Theory, Algorithmic Techniques and Applications, 27–52, Kluwer Academic Publishers, Dordrecht, The Netherlands, 2002.

5. V. C. Barbosa and E. Gafni, “Concurrency in heavily loaded neighborhood-constrained systems,” ACM Trans. on Programming Languages and Systems 11 (1989), 562–584.

6. K. M. Chandy and J. Misra, “The drinking philosophers problem,” ACM Trans. on Programming Languages and Systems 6 (1984), 632–646.

7. B. Charron-Bost, “Concerning the size of logical clocks in distributed systems,” Information Processing Letters 39 (1991), 11–16.

8. E. W. Dijkstra, “Hierarchical ordering of sequential processes,” Acta Informatica 1 (1971), 115–138.

9. R. P. Dilworth, “A decomposition theorem for partially ordered sets,” Annals of Mathematics 51 (1950), 161–165.

10. L. M. A. Drummond and V. C. Barbosa, “Distributed breakpoint detection in message-passing programs,” J. of Parallel and Distributed Computing 39 (1996), 153–167.

11. C. J. Fidge, “Timestamps in message-passing systems that preserve partial ordering,” Proc. of the 11th Australian Computer Science Conference, 56–66, 1988.

12. C. J. Fidge, “Logical time in distributed computing systems,” IEEE Computer 24 (1991), 28–33.

13. C. J. Fidge, “Fundamentals of distributed system observation,” IEEE Software 13 (1996), 77–83.

14. V. K. Garg, Principles of Distributed Systems, Kluwer Academic Publishers, Boston, MA, 1996.

15. L. Lamport, “Time, clocks, and the ordering of events in a distributed system,” Comm. of the ACM 21 (1978), 558–565.

16. N. A. Lynch, “Upper bounds for static resource allocation in a distributed system,” J. of Computer and System Sciences 23 (1981), 254–278.

17. F. Mattern, “Virtual time and global states in distributed systems,” in M. Cosnard et alii (Eds.), Parallel and Distributed Algorithms: Proc. of the Int. Workshop on Parallel and Distributed Algorithms, 215–226, North-Holland, Amsterdam, The Netherlands, 1989.

18. M. Raynal, “Illustrating the use of vector clocks in property detection: an example and a counter-example,” in P. Amestoy et alii (Eds.), Euro-Par’99—Parallel Processing, 806–814, Lecture Notes in Computer Science 1685, Springer-Verlag, Berlin, Germany, 1999.

19. L. E. T. Rodrigues and P. Veríssimo, “Causal separators for large-scale multicast communication,” Proc. of the 15th Int. Conf. on Distributed Computing Systems, 83–91, 1995.

20. F. Ruget, “Cheaper matrix clocks,” in G. Tel and P. Vitányi (Eds.), Distributed Algorithms: Proc. of the 8th Int. Workshop on Distributed Algorithms, 355–369, Lecture Notes in Computer Science 857, Springer-Verlag, Berlin, Germany, 1994.

21. S. K. Sarin and L. Lynch, “Discarding obsolete information in a replicated database system,” IEEE Trans. on Software Engineering SE-13 (1987), 39–46.

22. M. Singhal and A. Kshemkalyani. “An efficient implementation of vector clocks,” Information Processing Letters 43 (1992), 47–52.

23. J. L. Welch and N. A. Lynch, “A modular drinking philosophers algorithm,” Distributed Computing 6 (1993), 233–244.

24. G. T. J. Wuu and A. J. Bernstein, “Efficient solutions to the replicated log and dictionary problems,” Proc. of the 3rd Annual ACM Symposium on Principles of Distributed Computing, 233–242, 1984.


format: text

text:
# Quantum lower bounds for the set equality problems

## 1 Introduction, Motivation and Results

The Shor’s integer factoring quantum algorithm provides exponential speed-up over the best known classical algorithm. This motivates to search other quantum algorithms with great speed-up. However, proving quantum lower bounds for such problems is not trivial, for example, proving the exponential quantum lower bound for NP-Complete problems will imply $`PNP`$.

One of the problems quantum computer could have an exponential speed-up over classical computer is the Graph Isomorphism problem. One way to attack this problem could be by the reduction to the set equality problem. Notice the sets of all permutations over vertexes for given graphs. If these sets are equal, then there is an isomorphism between the graphs, but if there is not isomorphism between graphs, then these sets are strictly disjoint.

Denote the set $`\{1,2,\mathrm{},n\}`$ by $`[n]`$.

###### Definition 1

Let $`a:[n][m]`$ and $`b:[n][m]`$ be a functions. Let $`A`$ be a set of all $`a^{}s`$ images $`A=\{a(1),a(2),\mathrm{},a(n)\}`$ and $`B=\{b(1),b(2),\mathrm{},b(n)\}`$. There is a promise that either $`A=B`$ or $`AB=Ø`$.

Call the general set equality problem to distinguish these two cases.

Finding quantum query lower bound for general set equality problem was posed an open problem by Shi.

We will show that Ambainis’  adversary method imply $`\mathrm{\Omega }(\sqrt{n})`$ lower bound for the general set equality problem. The proof uses the possibility to have many preimages for some image. However, graph theorists think that the Graph Isomorphism problem, when graphs are promised not to be equal with themselves by any nonidentical permutation, still is very complex task. Now reduction lead us to the set equality where $`a`$ and $`b`$ are one-to-one functions.

###### Definition 2

Call the general set equality problem to be a one-to-one set equality problem if $`a(i)a(j)`$ and $`b(i)b(j)`$ for all $`ij`$.

The proof that worked for the general set equality problem does not work for one-to-one set equality problem, because it uses that fact that there can be very many preimages for any element of the sets. However, we will prove lower bound for a problem between these problems.

###### Definition 3

Call the general set equality problem to be a f(n) set equality problem if $`|a^1(x)|=O(f(n))`$ and $`|b^1(x)|=O(f(n))`$ for all images $`xa([n])b([n])`$ and for some function $`f`$.

We will prove $`\mathrm{\Omega }(\frac{n^{1/3}}{\mathrm{log}^{1/3}n})`$ lower bound for the $`log(n)`$ set equality problem.

The first result for lower bounds of the set equality like problem was done by Aaronson . He showed $`\mathrm{\Omega }(n^{1/6})`$ lower bound for so called set comparison problem: to decide whether two sets are equal or disjoint on a constant fraction of elements. He also assumed that both $`a`$ and $`b`$ are one-to-one functions. In this paper, we will study lower bound of problem when these sets $`A`$ and $`B`$ are strictly disjoint or equal, however $`a`$ and $`b`$ is not a one-to-one.

## 2 Preliminaries

### 2.1 Quantum Query algorithms

The most popular model of quantum computing is a query (oracle) model where the input is given by a black box. For more details, see a survey by Ambainis  or textbook by Gruska . In this paper we are able to skip them because our proof will be built on reduction to solved problems.

One of the most amazing quantum algorithms is a Grover’s search algorithm. It shows how a given $`x_1\{0,1\},x_2\{0,1\},\mathrm{},x_n\{0,1\}`$ to find the $`i`$ that $`x_i=1`$ with $`O(\sqrt{n})`$ queries.

This algorithm can be generalized to so called amplitude amplification . Using amplitude amplification one can make good quantum algorithms for many problems till the quadratic speed-up over classical algorithms.

By straightforward use of amplitude amplification we get a quantum algorithm with $`O(\sqrt{n})`$ queries for the general set equality problem and a quantum algorithm with $`O(n^{1/3})`$ queries for the one-to-one set equality problem.

### 2.2 Quantum query lower bounds

There are two main approaches to get good quantum lower bounds. The first is Ambainis’  quantum adversary method. The other is lower bound by polynomials introduced by Beals et al.  and substantially generalized by Aaronson  and Shi . Although explicitly we will use only Ambainis’ method, main result we will get by a reduction to problem, solved by polynomials’ method.

The basic idea of adversary method is that, if we can construct relation $`RA\times B`$, where $`A`$ and $`B`$ consisting of 0-instances and 1-instances and there is a lot of ways how to get from an instance in $`A`$ to an instance in $`B`$ that is in the relation and back by flipping various variables, then query complexity must be high.

###### Theorem 1

 Let $`f(x_1,\mathrm{},x_N)`$, be a function of n {0, 1}-valued variables and $`X,Y`$ be two sets of inputs such that $`f(x)f(y)`$ if $`xX`$ and $`yY`$. Let $`RX\times Y`$ be such that

* For every $`xX`$, there exist at least $`m`$ different $`yY`$ such that $`(x,y)R`$.
* For every $`yY`$, there exist at least $`m^{}`$ different $`xX`$ such that $`(x,y)R`$.
* For every $`xX`$ and $`i\{1,\mathrm{},n\}`$, there are at most $`l`$ different $`yY`$ such that $`(x,y)R`$ and $`x_iy_i`$.
* For every $`yY`$ and $`i\{1,\mathrm{},n\}`$, there are at most $`l^{}`$ different $`xX`$ such that $`(x,y)R`$ and $`x_iy_i`$.

Then, any quantum algorithm computing $`f`$ uses $`\mathrm{\Omega }(\sqrt{\frac{mm^{}}{ll^{}}})`$ queries.

### 2.3 The collision problem

Finding $`w(1)`$ quantum lower bound for the collision problem was an open problem since 1997. In 2001 Scott Aaronson  solved it showing polynomial lower bound. Later his result was improved by Yaoyun Shi . Newly Shi’s result was extended by Samuel Kutin  and by Andris Ambainis  in another directions.

Below is exact formulation of collision problem due to Shi.

###### Definition 4

Let $`n>0`$ and $`r2`$ be integers with $`r|n`$, and let a function of domain size n be given as an oracle with the promise that it is either one-to-one or r-to-one. Call the r-to-one collision problem the problem to distinguishing these two cases.

Shi  showed quantum lower bound for r-to-one collision problem.

###### Theorem 2

 Any error-bounded quantum algorithm to solve r-to-one collision must evaluate the function $`\mathrm{\Omega }((n/r)^{1/3})`$ times.

## 3 Results

### 3.1 Lower bound for the general set equality problem

###### Theorem 3

Any quantum algorithm which solves the general set equality problem makes $`\mathrm{\Omega }(\sqrt{n})`$ queries.

*Proof.* Simple use of Ambainis’ Theorem 1. Since Ambainis’ Theorem 1 deals with boolean functions, we will modify any quantum algorithm that solves the general set equality problem to an algorithm, that computes boolean function.

We will prove this theorem even in a restricted case, when functions returns only two values, let say $`0`$ and $`1`$. So we have a problem, given two functions $`a:n\{0,1\}`$ and $`b:n\{0,1\}`$ answer either the sets $`A=\{a(1),\mathrm{},a(n)\}`$ and $`B=\{b(1),\mathrm{},b(n)\}`$ are equal or disjoin under the promise that one of this is the case.

Let $`f:\{0,1\}^{2n}\{0,1\}`$ be partially defined function, such that

$`f(a_1,a_2,\mathrm{},a_n,b_1,b_2,\mathrm{},b_n)=\{\begin{array}{cc}1,\hfill & \text{if }\{a_1,\mathrm{},a_n\}=\{b_1,\mathrm{},b_n\}\text{;}\hfill \\ 0,\hfill & \text{if }\{a_1,\mathrm{},a_n\}\{b_1,\mathrm{},b_n\}=Ø\text{.}\hfill \end{array}`$

It is easy to see, that if we can solve a general set equality problem, we can compute this function with constant slowdown, too.

Let construct the relation $`R`$ from Ambainis’ Theorem 1 with $`X=\{0^n1^n\}`$ and $`Y=\{0^i10^{ni1}1^i01^{ni1}:0i<n\}`$ as follows:

$$R=X\times Y=\{(0^n1^n,0^i10^{ni1}1^i01^{ni1}):0i<n\}.$$

One can check that $`R`$ is well defined and $`m=n`$, $`m^{}=1`$, $`l=1`$ and $`l^{}=1`$. Thus any quantum algorithm computing $`f`$ makes $`\mathrm{\Omega }(\sqrt{n})`$ queries. $`\mathrm{}`$

### 3.2 Lower bound for the log(n) set equality problem

Now we will prove the main result in this paper.

###### Theorem 4

Any error-bounded quantum algorithm which solves the log(n) set equality problem makes $`\mathrm{\Omega }(\frac{n^{1/3}}{\mathrm{log}^{1/3}n})`$ queries.

Proof:

To prove Theorem 4 we will reduce r-to-one collision problem to the log(n) set equality problem. We are given function $`f:[n][m]`$, under promise to be either r-to-one or one-to-one and $`r=\mathrm{log}n`$ and $`r|n`$. We randomly choose two sets $`A`$ and $`B`$ such that $`|A|=|B|=n/2`$ and $`AB=[n]`$ and $`AB=Ø`$. Denote $`A^{}=f(A)`$ and $`B^{}=f(B)`$. It is obviously that if $`f`$ is one-to-one then $`A^{}B^{}=Ø`$.

If $`f`$ is r-to-one then the situation is more complicate. In the next subsection we will prove that with big probability holds that $`A^{}`$ and $`B^{}s`$ includes all images of f, thus $`A^{}=B^{}=f([n])`$.

Let the functions $`a`$ and $`b`$ from Theorem 4 be the same as $`f`$ but domain for $`a`$ is $`A`$ and domain for $`b`$ is $`B`$.

Denote the set of all preimages of x in the set $`A`$ by $`f_{}^{1}{}_{A}{}^{}(x)=f^1(x)A`$. Since $`|f^1(x)|=r`$ for every $`xf([n])`$, it is clear that also $`|a^1(x)|=O(\mathrm{log}n)`$ and $`|b^1(x)|=O(\mathrm{log}n)`$ for every $`xf([n])`$.

So with constant probability we get the log(n) set equality problem with domain size $`n/2`$. Now Theorem 2 implies Theorem 4. $`\mathrm{}`$

### 3.3 Reduction

###### Lemma 1

From all possible divisions of a set $`[n]`$ into two equal sized parts $`A`$ and $`B`$ such that $`AB=Ø`$, only few of them are such that for some $`xf([n])`$ there is no preimage either in $`A`$ or $`B`$.

Proof:

Total count of all (possibly uniform) divisions are

$$𝒞_n^{n/2}=\frac{n!}{(n/2)!(n/2)!}.$$

Total count of such divisions is at most the count of images ($`n/r`$) multiplied by count of divisions where one fixed $`xf([n])`$ has no preimage either in $`A`$ or $`B`$.

Assume that all preimages of $`x`$ is in $`A`$, thus $`B`$ has not any of them. Number of ways how we can choose residual elements is

$$𝒞_{nr}^{n/2r}=\frac{(nr)!}{(n/2r)!(n/2)!}.$$

Analysis of an opposite assumption is similar, so probability to choose division which is bad on $`x`$ is

$$\frac{2𝒞_{nr}^{n/2r}}{𝒞_n^{n/2}}=\frac{2(nr)!(n/2)!(n/2)!}{(n/2r)!(n/2)!n!}=\frac{2(n/2)(n/21)(n/22)\mathrm{}(n/2r+1)}{n(n1)(n2)\mathrm{}(nr+1)}=$$

$$=\frac{n/21}{n1}\frac{n/22}{n2}\mathrm{}\frac{n/2r+1}{nr+1}(\frac{1}{2})^{r1}.$$

So probability to choose bad division for any $`xf([n])`$ is at most

$$(\frac{1}{2})^{r1}\frac{n}{r}=\frac{2n}{2^rr}.$$

Since $`r=\mathrm{log}n`$

$$\frac{2n}{2^rr}\frac{2}{\mathrm{log}n}$$

which is small for large $`n`$.

$`\mathrm{}`$

## 4 Conclusion

Finding lower bound for the set equality problem is one of the most challenging today’s task in theory of quantum query lower bounds. We have solved this problem partially. One can argue that to solve the set equality problem can be easier when functions are promised to be with a small range of preimages for all images. Our paper shows that the difference between the general set equality problem and the $`\mathrm{log}n`$ set equality problem is very small, respectively $`\mathrm{\Omega }(\sqrt{n})`$ and $`\mathrm{\Omega }(\frac{n^{1/3}}{\mathrm{log}^{1/3}n})`$. This enforce opinion that quantum computer probably cannot solve the one-to-one set equality problem with only polylogarithmic number of questions.

## 5 Acknowledgments

I am very grateful to Andris Ambainis for introducing me with quantum algorithms and also with this problem as well as comments during writing this paper.

Research supported by Grant No.01.0354 from the Latvian Council of Science, and Contract IST-1999-11234 (QAIP) from the European Commission.


text:
# A Study of the Branching Ratio of 𝐻→𝑐𝑐̄ at a Future 𝑒⁺𝑒⁻ Linear Collider

## I Introduction

In the theory of elementary particles, Higgs boson is introduced to elucidate the origin of particle masses. The Standard Model (SM) proposed the single Higgs doublet that gives rise to a scalar particle, while extended supersymmetric models hypothesized that two Higgs doublets give separate vacuum expectation values to the up-type and down-type quarks. Among the extended models beyond SM, the Minimal Supersymmetric Standard Model (MSSM) sets the upper bound of the lightest neutral Higgs mass around 130 GeV/$`c^2`$susy . Since the mass of SM Higgs boson is expected to be between 114.4 GeV/$`c^2`$ and 211 GeV/$`c^2`$ by LEP experimentlep , the Higgs boson should be identified if it is found. Therefore a precise measurement of Higgs couplings to fermions at a future $`e^+e^{}`$ linear collider is indispensable to unveil physics of the Higgs sector. In particular, a measurement of the branching ratio of $`Hc\overline{c}`$ would provide a unique opportunity to study Higgs to up-type quark coupling. However, the $`Hb\overline{b}`$ is the dominant process for Higgs mass below 140 GeV/$`c^2`$ .

A high efficient charm quark identification is necessary for the precise branching ratio study and it is possible only with a high performance vertex detector. Development of this vertex detector has been a major issue for the linear collider lcrnd . To this end, understanding detector performance is extremely important to physics.

In this paper, we introduce our method of measuring the branching ratio of $`Hc\overline{c}`$ against SM $`e^+e^{}`$ processes and other decay branches from Higgs. In addition, we would like to see the significance of this precision study with four layers of Charge Coupled Device (CCD) and an additional inner layer of vertex detector.

## II HIGGS SIMULATION

A mass of SM Higgs boson is assumed to be 120 GeV/$`c^2`$ and the center-of-mass energy of 250 GeV is selected by the highest $`s`$-channel production cross-section at this energy, 226 fb. It is concerned with an integrated luminosity of 500 $`\mathrm{fb}^1`$ for the first few years running of the future $`e^+e^{}`$ linear collider experiment.

Events were generated using PYTHIA 5.7 pythia  and only the Higgsstrahlung process ($`e^+e^{}Z^{}Z^{}H`$) was taken into account. The branching ratios of the SM Higgs boson was estimated by the HDecay program hdecay . Depending on the decay modes of $`Z^{}`$, the events are categorized into 4-jet mode ($`Z^{}q\overline{q}`$ and $`Hq\overline{q}`$), 2-jet mode ($`Z^{}\nu \overline{\nu }`$ and $`Hq\overline{q}`$), and charged lepton pair mode ($`Z^{}\mathrm{}^+\mathrm{}^{}`$ and $`Hq\overline{q}`$). Since the branching ratio of Higgs to $`c\overline{c}`$ is estimated to be very small ($``$3%), we concentrated on the 2-jet and 4-jet modes only. As background processes from $`e^+e^{}`$ collision, $`W^+W^{},Z^{}Z^{}`$, and $`q\overline{q}`$ events were considered with corresponding cross-sections of 15460 fb, 1250 fb, and 47300 fb, respectively.

The detector simulation was performed using a fast parameterized simulator jlc , which is implemented in the Joint Linear Collider (currently called as Global Linear Collider) detector report . In this simulator five parameters of the helical track and their error matrices including non-diagonal elements were generated; thus the quality of the vertices is similar to that of a full simulation. The JLC vertex detector is equipped with four layers of CCD at the radius from 2.4 cm to 6 cm of which intrinsic spatial resolutions are 4 $`\mu \mathrm{m}`$ in $`\mathrm{r}\varphi `$ and z directions. The solenoidal magnet field of the detector is 3 tesla. With a vertex detector constraint, the impact parameter resolution in xy plane ($`\sigma _{\mathrm{r}\varphi }`$) is $`\sqrt{(25/\mathrm{psin}^{2/3}\theta )^2+4^2}\mu \mathrm{m}`$, and the momentum resolution for charged track ($`\mathrm{\Delta }\mathrm{p}_\mathrm{t}/\mathrm{p}_\mathrm{t}`$) is $`\sqrt{(1\times 10^4\mathrm{p}_\mathrm{t})^2+10^3}`$. This study was done with four/five CCD layers of vertex detector parameter, repectively, to see the influence of the vertex detector options on the physics result. The vertex detector parameter sets are shown in the Table 1.

In the event reconstruction, jets were reconstructed using JADE clustering algorithm jade  and the vertices were reconstructed using the topological vertex reconstruction algorithm (ZVTOP program) zvtop . Based on kinematic variables and reconstructed vertex information, an artificial neural network (NN) nn  is formed to identify $`Hc\overline{c}`$ events more efficiently.

We considered the non-$`c\overline{c}`$ Higgs decays ($`Hb\overline{b},gg,WW^{}`$; Higgs background) and other processes from $`e^+e^{}`$ collision ($`e^+e^{}Z^{}Z^{},W^+W^{},q\overline{q}`$; non-Higgs background) as backgrounds of $`Hc\overline{c}`$ measurement.

## III Selection of $`Hc\overline{c}`$

### III.1 2-jet mode

In 2-jet mode all clustered particles are forced to make 2 jets by adjusting the maximum value of ycut 333$`\mathrm{y}_{\mathrm{ij}}=\frac{2\mathrm{E}_\mathrm{i}\mathrm{E}_\mathrm{j}(1\mathrm{cos}\theta _{\mathrm{ij}})}{\mathrm{E}_{\mathrm{visible}}^2}`$, where $`\mathrm{E}_\mathrm{i}\mathrm{and}\mathrm{E}_\mathrm{j}`$ are energy of i-th and j-th cluster, respectively, and $`\theta _{ij}`$ is angle between two clusters., Ymax. The following selections are applied to reduce non-Higgs background in the sample: (1)visible energy between 110 GeV and 143 GeV, (2)missing transverse momentum between 25 GeV/$`c`$ and 70 GeV/$`c`$, (3)Higgs mass between 105 GeV/$`c^2`$ and 125 GeV/$`c^2`$, (4)recoiled $`Z^{}`$ mass between 82 GeV/$`c^2`$ and 120 GeV/$`c^2`$, (5)thrust between 0.75 and 0.99, (6)Ymax between 0.72 and 0.84, and (7)mass of each jet between 2 GeV/$`c^2`$ and 40 GeV/$`c^2`$. Here the Higgs mass is an invariant mass of all the observed particles and the recoiled $`Z^{}`$ mass is calculated disregarding the initial state radiation. We also require a successful secondary vertex reconstruction and at least one $`\mathrm{P}_\mathrm{t}`$ corrected invariant mass (MSPTM, $`\sqrt{\mathrm{M}_{}^{2}{}_{\mathrm{VTX}}{}^{}+\mathrm{P}_{\mathrm{t}}^{}{}_{}{}^{2}}+|\mathrm{P}_\mathrm{t}|`$) of the secondary vertex to be between 0.1 GeV/$`c^2`$ and 7.0 GeV/$`c^2`$, where $`\mathrm{P}_\mathrm{t}`$ is the total transverse momentum of the secondary tracks with respect to the flight direction of the vertex. The separation between signal and non-Higgs background is clearly seen in Fig. 1, while Higgs background is not distinguishable in this aspect.

After the sample selection, NN is trained to maximize the selection quality. Two sets of training are used; one against non-Higgs backgrounds (Background NN-training) and the other against Higgs backgrounds ( Higgs NN-training). For each NN-training $`HWW^{}`$ and $`e^+e^{}q\overline{q}`$ processes are not used due to the featureless patterns and small statistics, respectively.

The normalized input patterns for NN-trainings were given from the ZVTOP program; the number of vertices, MSPTM, transverse momentum of the secondary vertex, decay length, invariant mass of the secondary vertex, number of tracks in the secondary vertex, and momentum of the secondary vertex divided by total momentum (corrected secondary momentum). Distinctive input patterns for Higgs NN-training and Background NN-training are shown in Fig. 2 and Fig. 3, respectively. In the Background NN-training Ymax and jet decay angle were included in the input patterns especially for this 2-jet mode.

Using both results of Higgs NN-training and Background NN-training, the signal region is evaluated in two dimensional plane for best significance. The signal selection area is shown in the top right corner of the Fig. 4.

### III.2 4-jet mode

Here we used forced 4-jet clustering for Higgs jet pair and $`Z^{}`$ jet pair. The jets are identified by minimum $`\chi ^2`$, which is defined as the squared sum of the differences between reconstructed invariant mass of jet pair and the expected mass divided by each mass resolution:

$`\chi ^2=({\displaystyle \frac{M_Z^{}91.2}{\mathrm{width}}})^2+({\displaystyle \frac{M_{\mathrm{recoil}Z^{}}91.2}{\mathrm{width}}})^2`$
$`+({\displaystyle \frac{M_H120}{\mathrm{width}}})^2+({\displaystyle \frac{M_{\mathrm{recoil}H}120}{\mathrm{width}}})^2.`$
The selection cuts were made: (1) visible energy greater than 210 GeV, (2) $`\mathrm{cos}\theta `$ of thrust axis(thrust angle) between $`0.85`$ and 0.85, (3) Higgs mass between 103 GeV/$`c^2`$ and 130 GeV/$`c^2`$, (4) $`Z^{}`$ mass between 78 GeV/$`c^2`$ and 102 GeV/$`c^2`$, (5) mass recoiled to the Higgs jet pair greater than 80 GeV/$`c^2`$, (6) number of particles in each jet of Higgs pair greater than 5, and (7) Ymax value greater than 0.01.

Only 2 jets which are assigned to Higgs are used for the vertex tagging. We also required at least one successful secondary vertex reconstruction with MSPTM between 0.1 GeV/$`c^2`$ and 7.0 GeV/$`c^2`$ out of 2 jets from Higgs. The distributions of Higgs mass and thrust angle are shown in the Fig. 5. With the same procedure as the 2-jet study, we trained the NN against non-Higgs background and Higgs background separately using normalized vertex information from the ZVTOP program. The MSPTM and the corrected secondary momentum distributions of the $`c\overline{c}`$ and the non-$`c\overline{c}`$ Higgs decays are seen in the Fig. 6. We did not use $`HWW^{}`$ events in Higgs NN-training due to their featureless patterns. In the Background NN-training, thrust angle is added in the input pattern. As seen in the Fig. 7, the thrust angle shows a clear separation between the signal and non-Higgs backgrounds. We use both results of Higgs NN-training and Background NN-training same as 2-jet mode. The 2D NN-results are shown in the Fig. 8 with signal box at the top right.

## IV Discussion

The results of $`c\overline{c}`$ event selection are summarized in Table 2. As seen in the table, roughly 120 and 310 $`Hc\overline{c}`$ events are selected in 2-jet and 4-jet modes, respectively, with a reasonable significance despite its small branching ratio of $`Hc\overline{c}`$ and large backgrounds from non-Higgs processes and non-$`c\overline{c}`$ Higgs decays.

In the $`Hc\overline{c}`$ measurement in 2-jet mode, major backgrounds are those from $`Hb\overline{b}`$ where $`b`$ is misidentified as $`c`$ and $`W^+W^{}`$ events which are reconstructed as 2-jet events due to imperfect detector acceptance. $`e^+e^{}e\nu W`$ is not considered here due to relatively small cross-section than other background processes.

In the case of 4-jet mode, $`W^+W^{}`$ and $`Z^{}Z^{}`$ events are increased after event selection because two $`c`$ jets can be produced in the final state of these processes. In addition to the c-jet contamination from backgrounds, there is additional ambiguity in selecting two Higgs jets out of four jets. Combining these effects, the significance in the 4-jet mode is less than 3$`\sigma `$ and worse than the 2-jet mode. The statistical uncertainty combining 2-jet and 4-jet analyses is estimated as 25.7% in the case of four CCD layers of the vertex detector. Similar analysis with different detector configuration was reportJB  for assumption of the same Higgs mass and luminosity at $`\sqrt{s}`$=350 GeV and 500 GeV that the relative braching ratio errors be 19% and 39%, respectively.

When an additional CCD layer of the vertex detector is included near the interaction point, the impact parameter resolution for low momentum track is improved since the lever arm for track extrapolation to the interaction point is reduced. Thus we can reconstruct secondary vertices much closer to the interaction point. The better separation of $`b`$\- and $`c`$\- jets reduces backgrounds in the event selection, especially in 2-jet mode. Improvements in 4-jet selection is not decisive, suffering ambiguities in Higgs jet selections. Combining 2-jet and 4-jet, relatively 20% improvement in background reduction is achieved compared to the study with four CCD layers of the vertex detector.

## V Summary

In this study we focused on the $`Hc\overline{c}`$ measurements in 2-jet and 4-jet modes in the case of 120 GeV/$`c^2`$ Higgs mass at the center-of-mass energy of 250 GeV in the future $`e^+e^{}`$ linear collider. In the study, the topological vertex finding algorithm was used for tagging the $`c`$-jet and the neural network was used to optimize the $`Hc\overline{c}`$ selection.

With consideration of 500 $`\mathrm{fb}^1`$ data, we obtained the statistical uncertainty of 25.7% for the measurement of $`Hc\overline{c}`$ with four CCD layers vertex detector. The statistical uncertainty would be improved to 20.1% with an additional CCD layer of vertex detector near the interaction point.

## VI Acknowledgment

We would like to thank the members of ACFA Joint Linear Collider Physics and Detector Working group acfa  for valuable discussions during the course of the analysis. We also thank Dr. David Jackson for allowing us to use his ZVTOP program. We acknowledge the support by Korea Research Foundation Grant funded by Korea Government (MOEHRD) (KRF-2005-070-C000321) and Korea Science and Engineering Foundation. Akiya Miyamoto is partially supported by Japan-Europe Research Cooperative Program.


text:
# An Algorithm for Optimal Partitioning of Data on an Interval

## I Introduction: The Problem

A variety of signal processing and related problems can be viewed as the search for an optimal partition of data given on a time interval $`I`$. For example, one may estimate a segmented model by maximizing some measure of model fitness111This concept goes by many names, including *fitness*, *cost*, *goodness of fit*, *loss*, *penalty*, *objective function*, *risk* etc. Here we assume that the data analysis problem can be phrased in terms of some fitness that is to be maximized. defined on partitions of $`I`$. Since the space of all partitions of a continuum is infinite, it is advantageous to discretize the interval. Often the data points themselves, say

$$X_n,n=1,2,\mathrm{},N,$$
(1)
naturally subdivide $`I`$ into subintervals – which we call data cells. We avoid a precise definition because many types of data cells are possible. Common examples are counts in bins, measurements at a set of sample times (evenly spaced or not), and event or point data. The underlying idea is that restricting consideration to the finite space of partitions whose elements are sets of data cells will result in no significant loss of information or of resolution in the independent variable.

A partition $`𝐏`$ of an interval $`I`$ is a set of $`M`$ blocks,

$$𝐏(I)=\{B_m,m\text{M}\},\text{M}\{1,2,\mathrm{}M\},$$
(2)
where the blocks are sets of data cells defined by index sets $`\text{N}_m`$:

$$B_m=\{X_n,n\text{N}_m\}$$
(3)
satisfying the usual conditions, $`\underset{m}{}B_m=I`$ and $`B_mB_m^{}=\mathrm{}`$ if $`mm^{}`$. Here $`\mathrm{}`$ means a set of zero area, since the 1D boundaries are not relevant here. We impose connectedness too—i.e. that there be no gaps between the cells comprising a given block. $`M`$, the number of blocks, must satisfy $`0MN`$. Partitions will be denoted in boldface, and refer to the interval $`I`$ unless otherwise stated. Define $`P^{}`$ to be the (finite) set of all partitions of $`I`$ into blocks.

Take as given an additive fitness function that assigns a value to any partition $`𝐏P^{}`$ in the form

$$V(𝐏)=\underset{m=1}{\overset{M}{}}g(B_m),$$
(4)
where $`g(B_m)`$ is the fitness of block $`B_m`$. Computationally, the data cells must be represented by a data structure that contains sufficient statistics for the model – i.e. all information necessary to determine $`g`$ for any block.

We exhibit an efficient $`O(N^2)`$ dynamic programming algorithm that finds an optimal partition $`𝐏^{\text{max}}P^{}`$: $`V(𝐏^{\text{max}})V(𝐏)`$ for all partitions $`𝐏P^{}`$.

Scargle  proposed two greedy iterative algorithms for finding near-optimal partitions: one top-down (optimally divide $`I`$ into two parts; recursively do the same to each such part) the other bottom-up (merge adjacent data cells). In both cases Bayesian model comparison provides effective fitness functions and halting criteria, implementing an $`O(N^2)`$ procedure for data spaces of 1, 2 and higher dimensions—hence the term Bayesian Blocks . But in practice these greedy algorithms often find significantly suboptimal partitions, motivating the development reported here.

The dynamic programming idea for this kind of problem seems to have originated with Richard Bellman , a paper which influenced Kay’s work . An extensive discussion of precisely the same problem addressed here, but with a different approach to its solution, is in . Work by Hubert , with applications to meteorology, influenced Kehagias and co-workers , who developed a dynamic programming algorithm much like ours, for applications such as text segmentation (see also ), where the raw data are provided in the form of a similarity matrix.  gives an $`O(kN^2)`$ dynamic programming algorithm for finding the optimal partition of an interval into $`k`$ blocks, for a given $`k`$. See also  for related work. Thus while similar algorithms have been previously developed, ours finds the exact global optimum for any block-additive fitness function, automatically determines the number of segments, and can be used in either real-time or retrospective analysis.

## II Dynamic Programming: Finding Optimal Partitions

We describe an $`O(N^2)`$ algorithm that is guaranteed to solve the above problem by finding an exact global optimum, for any fitness function $`V`$ that is additive in the sense of Eq. (4). There is a large ($`2^{N1}`$) but finite number of partitions in $`P^{}`$. Dynamic programming  is an intelligent method of searching this space of all possible solutions. Our algorithm can be applied whenever any subpartition of an optimal partition is optimal.

###### Theorem 1 (Principle of Optimality)

Let $`𝐏^{\text{max}}`$ be an optimal partition of $`I`$ and $`𝐏_1=\{B_m,ma\}`$ be any subset of the blocks of $`𝐏^{\text{max}}`$. Then $`𝐏_1`$ is an optimal partition of the part of $`I`$ it covers, namely $`I_1=\underset{ma}{}B_m`$.

Intuitively this result follows from the contradiction that a better subpartition of $`I_1`$ could be used to construct a partition of $`I`$ better than $`𝐏^{\text{max}}`$. The proof relies on the fact that the block-additivity of the fitness function implies that it is also additive on subpartitions. To see this, divide partition $`𝐏`$ into any two disjoint parts, $`𝐏_1=\{B_m,ma\}`$ and $`𝐏_2=\{B_m,mb\}`$, with $`𝐏_1𝐏_2=𝐏`$ and $`ab=\text{M}`$. Then the additivity of $`V`$ yields

$`V(𝐏)`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{M}{}}}g(B_m)`$ (5)
$`=`$ $`{\displaystyle \underset{ma}{}}g(B_m)+{\displaystyle \underset{mb}{}}g(B_m)`$
$`=`$ $`V(𝐏_1)+V(𝐏_2).`$
Proof 1: As above, denote by $`𝐏_2`$ the subpartition of $`𝐏^{\text{max}}`$, consisting of the blocks $`\{B_m,m\text{M}a\}`$ in $`𝐏^{\text{max}}`$ that are not in $`𝐏_1`$. Let $`𝐏_3`$ be any other partition of $`I_1`$. Since $`𝐏^{\text{max}}`$ is an optimal partition of $`I`$ and $`𝐏_3𝐏_2`$ is also a partition of $`I`$ it follows that $`V(𝐏^{\text{max}})=V(𝐏_1)+V(𝐏_2)V(𝐏_3𝐏_2)=V(𝐏_3)+V(𝐏_2)`$ so $`V(𝐏_1)V(𝐏_3).`$ Thus $`𝐏_1`$ is an optimal partition of $`I_1`$.

Dynamic programming is a recursive procedure that can be used to efficiently find the solution to many kinds of combinatorial optimization problems. Our algorithm derives the optimal partition of the first $`n+1`$ data points using previously obtained optimal partitions, *i.e.* those of the first $`1,2,\mathrm{}n`$ data points. At each iteration we must consider all possible starting locations $`j`$, $`1jn`$ of the last block of the optimal partition. For each putative $`j`$ the fitness function is – by the principle of optimality – the fitness of the optimal subpartition prior to $`j`$ plus the fitness of the last block itself. The former was stored at previous iterations, and the latter is a simple evaluation of $`V`$. The desired new optimal partition corresponds to the maximum over all $`j`$.

More precisely, define opt$`(n)`$ to be the value of the fitness function of the optimal partition $`𝐏_n^{\text{max}}`$ of the first $`n`$ cells of $`I`$, for $`1nN`$. The following dynamic programming algorithm finds the optimal partition $`𝐏_N^{\text{max}}`$:

1. Define opt$`(0)=0`$
2. Given that opt$`(j)`$ has been determined for $`j=0,1,\mathrm{},n`$:

* Define end$`(j,n+1)=g(B_{j,n+1})`$; $`B_{j,n+1}`$ is the union of cells $`j,j+1,\mathrm{},n+1`$
    * Then compute

    $$\mathrm{𝑜𝑝𝑡}(n+1)=\underset{j}{\text{Max}}\{opt(j1)+\mathrm{𝑒𝑛𝑑}(j,n+1)\},$$

    for $`j=1,2,\mathrm{},n+1`$.
    * The value of $`j`$ where this maximum occurs is stored as lastchange$`(n+1)`$.
3. Repeat 2 until $`n+1=N`$, when opt$`(N)`$, the optimal partition fitness for all $`N`$ cells, has been obtained.
4. Backtrack using the lastchange vector to identify the start points of individual blocks of the optimal partition $`𝐏^{\text{max}}`$ in the following way. Let $`n_1=\mathrm{𝑙𝑎𝑠𝑡𝑐ℎ𝑎𝑛𝑔𝑒}(N),n_2=\mathrm{𝑙𝑎𝑠𝑡𝑐ℎ𝑎𝑛𝑔𝑒}(n_11)`$, etc. Then the last block in $`𝐏^{\text{max}}`$ contains cells $`n_1,n_1+1,\mathrm{},N`$, the next-to-last block in $`𝐏^{\text{max}}`$ contains cells $`n_2,n_2+1,\mathrm{},n_11,`$ and so on.

###### Theorem 2

This deterministic $`O(N^2)`$ dynamic programming algorithm finds the partition of $`I`$ that maximizes the (additive) fitness function.

Proof 2: The proof is by mathematical induction. Clearly $`\mathrm{𝑜𝑝𝑡}(1)=\text{Max}\{\mathit{0}+\mathrm{𝑒𝑛𝑑}(1,1)\}=g(B_{1,1})`$ is the fitness of the only possible (and therefore optimal) partition of the set comprising the first cell. At iteration $`n+1`$, assume not only that we have found the optimum partition of $`𝐏_n^{\text{max}}`$, but also that for $`i=1,2,\mathrm{},n,`$ we have stored the corresponding fitness for this and all previous iterations in the array $`\mathrm{𝑜𝑝𝑡}(i)`$, and the index of the cell beginning this partition’s last block in array $`\mathrm{𝑙𝑎𝑠𝑡𝑐ℎ𝑎𝑛𝑔𝑒}(i)`$. Let $`F(j)=opt(j1)+\mathrm{𝑒𝑛𝑑}(j,n+1)`$; then the principle of optimality shows that when $`j`$ indexes the first cell of the last block of the desired partition $`𝐏_{n+1}^{\text{max}},F(j)`$ is the corresponding maximum fitness. Further, for any $`j,F(j)`$ is the fitness of a legitimate partition of $`𝐏_{n+1}^{\text{max}}`$, namely that consisting of the optimal partition of the cells prior to $`j`$ followed by the single block $`B_{j,n+1}`$. These two facts combine to prove that the maximum of $`F(j)`$ specified in Eq. (2) gives the desired optimum partition at iteration $`n+1`$. Identification of the corresponding optimal blocks – starting with the last one and working backwards, as in part (4) of the algorithm – can be validated with straightforward recursive application of the principle of optimality. Finally, since end$`(j,n+1)=g(B_{j,n+1})`$ the algorithm requires $`1+2+\mathrm{}+N=O(N^2)`$ evaluations of the function $`g`$. It also requires $`O(N^2)`$ additions and $`O(N^2)`$ comparisons in determining the maximums.

## III Applications

These results apply to any segmented modeling of 1-dimensional data defined by a fitness function that satisfies Eq. (4). Piecewise constant, or step functions form the most natural model class. However the nature of the model depends on the application, and many other forms are possible, including piecewise linear and piecewise exponential. The key is that the fitness function must not depend on any model parameters other than the changepoint locations. We have found excellent results with the posterior probability of the model for each segment, given the data in that segment, marginalized over all parameters but these locations . Using logs turns the product posterior (resulting from statistical independence of the blocks) into the required additive form Eq.(4).

A comment about smoothness constraints is in order. With some fitness functions the algorithm produces the degenerate solution assigning a segment to each data point. In the Bayesian setting described above, a natural way to address this problem is to adopt a prior distribution for the number of segments, for example giving higher weight to smaller numbers. Indeed, the geometric prior  corresponds to a constant term in the fitness function for each block, so there is no problem with maintaining additivity. This artifice controls model complexity, but without imposing an explicit smoothness condition with concomitant loss of time resolution. Time series features on any time scale, no matter how short, can be found if the data support them in a statistically significant way.

Finally, we mention a few sample applications. Implementing density estimation with piecewise constant Poisson models yields histograms in which the bins are not constrained to be equal. The number of bins and their sizes and locations are determined by the data. The same model provides denoising and structure estimation for time series of events or counts of events in bins .

Further, almost all of the results described here can be easily extended – almost without change – to data of higher dimensionality, as will be described in future papers. In this setting cluster analysis can be effected as a post-processing of segmented models – piecing the blocks together into clusters – and similarly with unsupervised classification and other data mining procedures.

## IV Conclusion

As we have seen dynamic programming gives a good (polynomial) algorithm for finding an optimal partition of data on an interval for any fitness function $`V`$ satisfying the additive property \[see Eq.(4)\]. Ironically it has the same $`O(N^2)`$ complexity as the greedy algorithm.

In comparing the use of our algorithm to detect and characterize clusters (collections of blocks) with some of the standard clustering techniques , we note that our method inherently compares partitions that have different numbers of blocks, so the number of blocks is automatically determined by the data. This is to be contrasted with most standard clustering techniques, in which $`k`$, the fixed number of clusters must be specified ahead of time. One often seeks to minimize the maximum diameter (defined as the maximum distance between any pair of points in the cluster) of the clusters, or to maximize the minimum separation between the clusters. In dimension 1, there are well-known $`O(kN^2)`$ dynamic programming algorithms for finding the best partitions into k clusters. For dimension 2 and higher it is known that these standard problems are NP-complete. We don’t yet know if our problem is NP-complete in dimension 2 and higher.

In addition, considered as a density estimation or signal detection technique, our approach does not introduce any explicit smoothing of the data. Structure on any time scale, no matter how short, will be detected if it is supported by the data. While the parameter in the geometric prior discussed above controls to some extent the number of blocks – and thus affects the roughness of the optimized model – it is not explicitly a smoothing parameter. Another feature is that the incremental way the algorithm operates on the data makes a real-time mode trivial to implement. This mode has found to be very useful in the rapid detection of change points in a data stream. And since opt$`(i+1)`$ is calculated from opt$`(j),j=1,2,\mathrm{},i,`$ some of the necessary calculations can be performed as the data are still being collected. Also it is easy to modify the dynamic programming to yield the optimal partition with blocks of a minimum size (each block contains at least d data points, for a given positive integer d). These and other features are described in more detail at an algorithm repository at:
`http://astrophysics.arc.nasa.gov/~pgazis/CodeArchiveServer.html`

We are grateful to Steve Kay, Thanasis Kehagias, and the referees for helpful comments and pointers to relevant earlier work.


text:
# Double-transverse spin asymmetries at NLO11footnote 1Invited talk presented by W. Vogelsang at the “Conference on the Intersections of Particle and Nuclear Physics (CIPANP 2003)”, New York City, May 19-24, 2003.

BNL-NT-03/18

DO-TH 03/14

RBRC-330

## Abstract

We report on a next-to-leading order QCD calculation of the cross section and the spin asymmetry for isolated large-$`p_T`$ prompt photon production in collisions of transversely polarized protons. Corresponding measurements may be used at RHIC to determine the transversity parton distributions of the proton.

The partonic structure of spin-1/2 targets at the leading-twist level is characterized by the unpolarized, longitudinally polarized, and transversely polarized distribution functions $`f`$, $`\mathrm{\Delta }f`$, and $`\delta f`$, respectively ref:jaffeji ; ref:ralston . These non-perturbative parton densities can be probed universally in a multitude of inelastic scattering processes, for which it is possible to separate (“factorize”) the long-distance physics relating to nucleon structure from a partonic short-distance scattering that is amenable to QCD perturbation theory.

In contrast to the $`f`$ and $`\mathrm{\Delta }f`$, the “transversity” distributions $`\delta f`$ are unmeasured thus far. They are presently the focus of much experimental activity. For example, information should soon be gathered from transversely polarized proton-proton collisions at the BNL Relativistic Heavy Ion Collider (RHIC) ref:rhic . The potential of RHIC in accessing transversity in measurements of transverse double-spin asymmetries $`A_{\mathrm{TT}}`$ was examined in ref:attlo  for high transverse momentum $`p_T`$ prompt photon and jet production (for earlier studies, see ref:attold ; ref:artru ; ref:jaffesaito ). All of these calculations were performed only at the lowest order (LO) approximation for the underlying partonic hard-scattering. As is well known, next-to-leading order (NLO) QCD corrections are generally indispensable in order to arrive at a firmer theoretical prediction for hadronic cross sections and spin asymmetries. The NLO calculation for $`A_{\mathrm{TT}}^\gamma `$ for isolated high-$`p_T`$ prompt photon production, $`pp\gamma X`$, was recently completed ref:msv ; here we give a brief report on those results.

Interesting new technical questions arise beyond the LO in case of transverse polarization. Unlike for longitudinally polarized cross sections where the spin vectors are aligned with momentum, transverse spin vectors specify extra spatial directions, giving rise to non-trivial dependence of the cross section on the azimuthal angle of the observed photon. As is well-known ref:ralston , for $`A_{\mathrm{TT}}`$ this dependence is always of the form $`\mathrm{cos}(2\mathrm{\Phi })`$, if the $`z`$ axis is defined by the direction of the initial protons in their center-of-mass system (c.m.s.), and the spin vectors are taken to point in the $`\pm x`$ direction. Integration over the photon’s azimuthal angle is therefore not appropriate. On the other hand, standard techniques developed in the literature for performing NLO phase-space integrations usually rely on integration over the full azimuthal phase space, and also on the choice of particular reference frames that are related in complicated ways to the one just specified. In ref:msv , a new general technique was introduced which facilitates NLO calculations with transverse polarization by conveniently projecting on the azimuthal dependence of the matrix elements in a covariant way. The key point here is to recognize that the the factor $`\mathrm{cos}(2\mathrm{\Phi })`$ in the cross section actually results from the covariant expression

$$(p_\gamma ,s_a,s_b)=\frac{s}{tu}\left[2(p_\gamma s_a)(p_\gamma s_b)+\frac{tu}{s}(s_as_b)\right],$$
(1)
with $`s_a,s_b`$ the initial spin vectors and $`p_\gamma `$ the photon momentum. $``$ reduces to $`\mathrm{cos}(2\mathrm{\Phi })`$ in the hadronic c.m.s. frame. One may thus integrate over all phase space without obtaining a vanishing result if one simply multiplies the squared matrix element by the factor $`(p_\gamma ,s_a,s_b)`$. Integration over terms involving the $`s_a,s_b`$ can be carried out in a covariant way by using standard tensor decompositions. After this step, there are no scalar products involving the $`s_i`$ left in the squared matrix element. For the ensuing integration over all azimuthal phase space we can now employ techniques familiar from the corresponding calculations in the unpolarized and longitudinally polarized cases.

At NLO, there are two subprocesses that contribute for transverse polarization, $`q\overline{q}\gamma X`$ and $`qq\gamma X`$. The first one is already present at LO, where $`X=g`$. At NLO, one has virtual corrections to the Born cross section ($`X=g`$), but also $`23`$ real emission diagrams, with $`X=gg+q\overline{q}+q^{}\overline{q}^{}`$. For the second subprocess, $`X=qq`$.

Owing to the presence of ultraviolet, infrared, and collinear singularities at intermediate stages of the calculation, it is necessary to introduce a regularization. Our choice is dimensional regularization, that is, the calculation is performed in $`d=42\epsilon `$ space-time dimensions. Ultraviolet poles in the virtual diagrams are removed by the renormalization of the strong coupling constant. Infrared singularities cancel in the sum between virtual and real-emission diagrams. After this cancelation, only collinear poles are left. These result for example from a parton in the initial state splitting collinearly into a pair of partons, corresponding to a long-distance contribution in the partonic cross section. From the factorization theorem it follows that such contributions need to be factored into the parton distribution functions. In our calculations ref:msv , we have imposed on the photon the isolation cut proposed in ref:frixione . All final-state collinear singularities then cancel. The isolation constraint was implemented analytically by assuming a narrow isolation cone.

For our numerical predictions we model the $`\delta f`$ by assuming that the Soffer inequality ref:soffer  is saturated at some low input scale $`\mu _00.6\mathrm{GeV}`$. For $`\mu >\mu _0`$ the transversity densities $`\delta f(x,\mu )`$ are then obtained by solving the appropriate QCD evolution equations. Our numerical predictions apply for prompt photon measurements with the Phenix detector at RHIC. Figure 1 shows our results for the transversely polarized prompt photon production cross sections at NLO and LO for two different c.m.s. energies. The lower part of the figure displays the so called “$`K`$-factor”, $`K=d\delta \sigma ^{\mathrm{NLO}}/d\delta \sigma ^{\mathrm{LO}}`$. One can see that the NLO corrections are somewhat smaller for $`\sqrt{S}=500\mathrm{GeV}`$ and increase with $`p_T`$. The shaded bands in the upper panel of Fig. 1 indicate the uncertainties from varying the factorization and renormalization scales in the range $`p_T/2\mu _R=\mu _F2p_T`$. The solid and dashed lines are always for the choice where all scales are set to $`p_T`$, and so is the $`K`$ factor underneath. One can see that the scale dependence becomes much weaker at NLO, as expected. The corresponding spin asymmetries $`A_{\mathrm{TT}}^\gamma =d\delta \sigma /d\sigma `$ may be found in Fig. 2 of Ref. ref:msv ; they are generally smaller at NLO than at LO.

W.V. is grateful to RIKEN, Brookhaven National Laboratory and the U.S. Department of Energy (contract number DE-AC02-98CH10886) for providing the facilities essential for the completion of this work. This work is supported in part by the “Bundesministerium für Bildung und Forschung (BMBF)” and the “Deutsche Forschungsgemeinschaft (DFG)”.


text:
# Systematic effects induced on IBIS detectors by background and inhomogeneity of the spatial response Based on observations with INTEGRAL, an ESA project with instruments and science data centre funded by ESA member states (especially the PI countries: Denmark, France, Germany, Italy, Switzerland, Spain), Czech Republic and Poland, and with the participation of Russia and the USA.

## 1 Introduction

The IBIS instrument on board INTEGRAL (Winkler et al. (2003)) is the first large area, space born $`\gamma `$-ray telescope carrying a pixellated, multi-layer detector system (Ubertini et al. (2003)). Due to its large aperture and large collecting area, the scientific performances are strongly dependent on background reduction. Optimal levels of background intensity and stability in time have been reached after tuning and calibration of instrumental parameters. This has been achieved during in-flight commissioning (Ubertini et al. (2003)). In IBIS, different data acquisition modes working simultaneously ensure a good capability of background rejection. Important features are the selection of single and multiple events in the high energy detector (PICsIT, Labanti et al. (2003)) and of multi-layer coincidence events (Compton mode). Furthermore, the two detector layers ISGRI (the hard X-ray detector, Lebrun et al. (2003)) and PICsIT are actively shielded by a multi-module VETO system, based on one lateral and one bottom arrays of BGO blocks (Quadrini et al. (2003)). ISGRI is also shielded from the diffuse hard X-ray background outside the field-of-view by a composite passive shield protecting the whole telescope from mask to detector base (see e.g. Natalucci & Caroli (1996)).

The IBIS data processing software (Goldwurm et al. (2003)) is capable of recovering many time varying systematic effects as telemetry losses, noisy pixels, temporary switch-off of ISGRI modules, dead time etc., which are well quantified and properly taken into account in the analysis. Once the data have been corrected for these effects we still see important spatial structures, which affect both detector layers. These are easily detected as different average count-rates in detector modules/semi-modules, enhancements at the edge of the modules and local effects induced by the readout electronics. There are several factors known to produce this inhomogeneity of the spatial response (see also Stephen et al. (2003); Terrier et al. (2003)): a) a spatially dependent hadronic background component, depending on the external payload and spacecraft mass distribution; b) the efficiency of detection of single and multiple events, which is expected to be dependent on the detector position; c) the individual response of each detector module; d) the intrinsic response of each pixel as a function of its position, i.e. its proximity to a module, an ASIC (Application Specific Integrated Circuit) readout element and/or a dead pixel; e) variations of the performance of the VETO modules, induced especially by temperature variations. These effects are, in turn, energy dependent. Moreover a given pixel may undergo periods of unusual behaviour. For ISGRI, defective pixels including noisy elements which are switched off permanently or for a relatively short time period may induce local variations in the detector itself and even in PICsIT. PICsIT pixels have also intrinsic noise levels that could have local effects in ISGRI. PICsIT data are normally transmitted as spectral histograms (SI) already accumulated on board. Varying count rates of pixels can influence the local response of detectors for less than the typical SI integration time ($`2150s`$), and hence cannot be properly followed in PICsIT. For this reason, local corrections based on analysis of single pixel light curves cannot be always applied efficiently. Therefore, methods of flattening the images independently of the time behaviour of pixels are being studied (see Sect. 2).

The difficulty of correcting data by using maps obtained by long background exposures is mainly related to the background being variable with time, and to different systematics in the counting behaviour of pixels in ”observation” and ”background” exposures. In the following sections, an overview of the main problems related to background induced effects is given and some results are reported in more detail, especially regarding the PICsIT detector. For PICsIT, an important source of background is the excess counts caused by the decay of long-lived phosphorescence states of the CsI crystals (Hurley (1978)). These counts cannot be vetoed because they give rise to single events, separated by intervals of the order of a few $`10^5`$s. However, they can be easily identified and removed in photon by photon data, and their spatial distribution can be separately studied (see Sect. 4). Other spatial inhomogeneity effects which are background independent are described in Sect. 3.

## 2 Flat-fielding methods

For IBIS, flat-fielding prior to image deconvolution (Goldwurm et al. (2001)) is essential especially at high energies, as long as the size of the fluctuations in the background level, propagated through the imaging process exceed the observed celestial source intensity. The approach of applying a flat-field based on Fourier filtering is difficult, because a substantial power of the background distribution is seen at high spatial frequencies, where most of the signal from multiplexing the source counts by the coded mask is also found. An alternative approach is to characterize each pixel by averaging on pixel categories, defined on the basis of the pixel position in the detector, i.e. proximity to a module edge and/or to a noisy/dead pixel, position within an ASIC etc. This has the form of a model map applicable to a given energy band. The model is obtained from real data by measuring average count rates for each pixel category (see Fig. 1). Correction maps defined for eight standard energy bands, covering the range 0.2-6.5 MeV, have been used to subtract a background shape from the single event shadowgrams. The resulting detector image is defined by a subtraction of the model from the original image: $`D_{ff}`$=($`D_{raw}`$-$`Mf`$)+$`D_{raw}`$, where $`M`$ is the model map and $`f`$=$`D_{raw}`$/$`M`$ is the ratio between average pixel values in the raw image and model map.

The standard deviation in the pixel counts of the background subtracted shadowgram is compared to that of the original image, and a flat-field efficiency is evaluated as:

$`ϵ_{ff}`$ = ($`\sigma _{raw}\sigma _{ff}`$) / ($`\sigma _{raw}\sqrt{D_{raw}}`$)

We have used a model map based on the sum of 3 science windows, for a total exposure of 6450 s of an empty field observation, to flatten a series of spectral histograms obtained from a Crab exposure. In Fig. 2, the average flat-field efficiency obtained by this processing is plotted as a function of energy. We remark that in the first two energy bands, the low efficiency obtained is not intrinsic to the flat-field process, but is rather due to the effect of phosphorescence counts (see Sect. 4), causing the average deviates to be well above the Poisson statistical value. Apart from this effect, the efficiency is energy dependent.

Adding more input data in building the model maps is not found to improve the efficiency. This is probably due to the time variations in the background. Improvements to the model are currently under study, as long as more systematic effects are identified.

## 3 Response to $`\gamma `$-ray source

Besides the known effect of background non-uniformity, scattering of celestial source $`\gamma `$-rays associated with a non-ideal response of the event selection logic and/or non-ideal geometry of the detector can produce a non-homogeneity of the source count distribution. This component of the non-uniformity (which is intrinsically different from the background induced one) is being studied by Monte Carlo simulations using a numerical code based on a detailed geometric and mass model of the telescope, including electronics and event selection logic (Laurent et al. (2001)). For this purpose, we have simulated a parallel flux of $`\gamma `$-rays incident on the IBIS detector with zero $`\mathrm{deg}`$ inclination, at different monochromatic photon energies and for a power law source. The efficiency and response of single pixels are considered uniform in the simulation, so inhomogeneities cannot be ascribed to their local behaviour. For a monochromatic source, it was verified that the image is perfectly uniform when selecting events under the full-energy peak. This is expected as far as there is no re-distribution of the event energy in pixels other than the incidence pixel. In contrast, the count distribution shows differences, especially in the counting rate at the edge of modules, when the energy range selected covers part of the Compton continuum. This is true in particular for ISGRI (see Fig 3) and PICsIT multiple events. This non-uniformity should be ideally corrected by dividing the residual detector image after background subtraction or flat-field (see Sect. 2) by an efficiency map which takes into account the scattering induced effects. This is, of course, dependent on the input source spectrum. A set of reference efficiency maps are being computed for a limited set of input spectra models (based on power law or cutoff power laws) to be tested and compared. We consider, however, that correcting real data using a non-uniformity model of the source illumination is expected to be of secondary importance compared to an efficient background subtraction.

## 4 Effects of phosphorescence induced events

The PICsIT detector is found to show significant count-rate increase caused by fake events, induced by decay of phosphorescent states in the CsI crystals. These states are excited by passage of particles in the active body of the PICsIT detector, which occur at a rate $`1000`$ times lower than the normal background rate. One burst of fake events can be produced by either the passage of a primary particle, or by a secondary shower (cascade) producing a track or other figure of well-defined spatial pattern (Segreto et al. (2003)). All these counts are mostly detected at energies below $`300`$ keV and can be easily recognized in photon by photon data, in which we have full timing information. Since each individual passage of particles produces a relatively large number of these events, the count rate statistics is strongly modified. The two panels of Fig. 4 show the fluctuations of the summed count rate for 16 pixels of an ASIC, in the energy range 170-285 keV. On the upper plot, representing the total counts, the fluctuations are clearly non-Poissonian. Once the phosphorescence triggers are eliminated, the fluctuations are fairly compatible with Poissonian noise (lower panel). It is quite easy to recognize, by comparing the two plots, that some time bins are not affected at all by phosphorescence triggers. This effect can be quantified by computing the ratio, $`R`$, between the observed standard deviation of the counts and the Poisson statistical error, as a function of energy and for different integration times. In Fig. 5 (upper panel) are shown the values of $`R`$ obtained from the same set of data (16 pixels, $`8.2\times 10^4`$ s). From the figure, it is clearly seen that the loss of sensitivity per channel does not depend on integration time for $`\mathrm{\Delta }`$$`10`$s.

We have then investigated the effect of these fake events in spectral imaging data, which is the PICsIT standard mode of operation, using an empty field exposure (integration time: 2150s). As the background intensity always shows a long-term trend, expected values of rates were computed by a linear fit to the data, after selecting an interval of constant linear variation of 25 hours. For this test, we have used total count rates and verified that a large spreading of the values is also present in the energy channels below $`250`$ keV (see Fig. 5, lower panel).

The ratio $`R`$ changes significantly when rebinning in energy. For example, we estimate $`R4.5`$ and $`R1.7`$ for the energy bands 170-220 keV, and 220-280 keV respectively. These values are sensibly higher than the ones plotted in Fig. 5 for the single energy channels. This is actually expected, as each track produced by a particle leaves a number of phosphorescent counts which are widely distributed in energy, thus giving rise to highly correlated count rate fluctuations in different energy channels.

Finally, we investigated the spatial distribution of the fake events by using the photon by photon data. The images obtained in the energy band 170-285 keV clearly show that the non-uniformity in this range cannot be ascribed to the phosphorescence triggers (see Fig. 6).

## 5 Conclusions

Several effects which reduce the sensitivity of the imaging methods used by IBIS, especially in the PICsIT high energy layer have been described. Some of them are ascribed to image non-uniformity, and can be corrected by a large extent by a suitable flat-field. Some background components (namely, the delayed phosphorescence emission of the CsI crystals) are found to compromise the sensitivity in the energy range 180-250 keV, when IBIS is operated in its Standard Mode. Image non-uniformity problems are described in the framework of finding good methods for background subtraction. The presence of local systematic effects due to variable or unexpected behaviour of single pixels and/or ASICs is hampering the application of standard flat-field techniques. Present efforts are devoted to understanding/modelling this type of systematics, at the same time refining the already available methods of background correction.

###### Acknowledgements.

The IBIS project is partially granted by Italian Space Agency (ASI). AJB is funded by PPARC grant GR/2002/00446.