\documentstyle[jair,twoside,11pt,theapa,psfig]{article}

\input epsf %% at top of file

\jairheading{18}{2003}{445-490}{04/02}{05/03}
\ShortHeadings{Searching Bayesian network structures in the space of RPDAGs}
{Acid \& de Campos}

\firstpageno{445}

\def\qed{\hbox{${\vcenter{\vbox{   % this is the cube
    \hrule height 0.4pt\hbox{\vrule width 0.4pt height 5pt
    \kern5pt\vrule width 0.4pt}\hrule height 0.4pt}}}$}}
\def\leaderfill{\leaders\hbox to 0.5em{\hss.\hss}\hfill}


\newtheorem{theo}{Theorem}
\newtheorem{defin}{Definition}
\newtheorem{prop}{Proposition}
%\graphicspath{./figuras}

\begin{document}

\title{Searching for Bayesian Network Structures in the Space of Restricted Acyclic Partially Directed Graphs}

\author{\name Silvia Acid \email acid@decsai.ugr.es \\
        \name Luis M. de Campos \email lci@decsai.ugr.es \\ 
\addr Departamento de Ciencias de la Computaci\'on e I.A. \\
E.T.S.I. Inform\'atica, Universidad de Granada \\
18071 - Granada, SPAIN
}


\maketitle

\begin{abstract}

Although many algorithms have been designed to construct Bayesian network structures using different approaches and principles, they all employ only two methods: those based on independence criteria, and those based on a scoring function and a search procedure (although some methods combine the two). Within the score+search paradigm, the dominant approach uses local search methods in the space of directed acyclic graphs (DAGs), where the usual choices for defining the elementary modifications (local changes) that can be applied are arc addition, arc deletion, and arc reversal. In this paper, we propose a new local search method that uses a different search space, and which takes account of the concept of equivalence between network structures: restricted acyclic partially directed graphs (RPDAGs). In this way, the number of different configurations of the search space is reduced, thus improving efficiency. Moreover, although the final result must necessarily be a local optimum given the nature of the search method, the topology of the new search space, which avoids making early decisions about the directions of the arcs, may help to find better local optima than those obtained by searching in the DAG space. Detailed results of the evaluation of the proposed search method on several test problems, including the well-known Alarm Monitoring System, are also presented.
\end{abstract}

\section{Introduction} \label{intro}

Nowadays, the usefulness of {\it Bayesian networks}~\cite{pe88} in representing knowledge with uncertainty and efficient reasoning is widely accepted. A Bayesian network consists of a qualitative part, a {\it directed acyclic graph} (DAG), and a quantitative one, a collection of numerical parameters, usually conditional probability tables.
The knowledge represented in the graphical component is expressed in terms of dependence and independence relationships between variables. These relationships are encoded using the presence or absence of links between nodes in the graph. The knowledge represented in the numerical part quantifies the dependences encoded in the graph, and allows us to introduce uncertainty into the model. All in all, Bayesian networks provide a very intuitive graphical tool for representing available knowledge.

Another attraction of Bayesian networks is their ability to efficiently perform reasoning tasks \cite{Jensen96,pe88}. The independences represented in the DAG are the key to this ability, reducing changes in the knowledge state to local computations. In addition, important savings in storage requirements are possible since independences allow a factorization of the global numerical representation (the joint probability distribution).

There has been a lot of work in recent years on the automatic learning of Bayesian networks from data. Consequently, there are a great many learning algorithms which may be subdivided into two general approaches: methods based on {\it conditional independence tests} (also called {\it constraint-based}), and methods based on a {\it scoring function}\footnote{Some authors also use the term {\it scoring metric}.} and a {\it search} procedure.

The algorithms based on independence tests perform a qualitative study of the dependence and independence relationships among the variables in the domain, and attempt to find a network that represents these relationships as far as possible. They therefore take a list of conditional independence relationships (obtained from the data by means of conditional independence tests) as the input, and generate a network that represents most of these relationships. The computational cost of these algorithms is mainly due to the number and complexity of such tests, which can also cause unreliable results. Some of the algorithms based on this approach obtain simplified models  \cite{campos98,camposhuete97,Geiger90,Geiger93,Huete93}, whereas other are designed for general DAGs \cite{camposhuete00,Cheng97,Meek95,Spirtes93,Verma90,Wermuth83}.

The algorithms based on a scoring function attempt to find a graph that maximizes the selected score; the scoring function is usually defined as a measure of fit between the graph and the data. All of them use a scoring function in combination with a search method in order to measure the goodness of each explored structure from the space of feasible solutions. During the exploration process, the scoring function is applied in order to evaluate the fitness of each candidate structure to the data. Each algorithm is characterized by the specific scoring function and search procedure used. The scoring functions are based on different principles, such as entropy  \cite{kuta,chow,campos98,Rebane87}, Bayesian approaches \cite{buntine94,buntine96,CooHersk:k2,friedmankolleruai00,friedmannpuai99,geigerheck,heckerman96,heckgeigerchick,madigan94,ramoni,steck00}, or the Minimum Description Length \cite{Bouckaert93,friedman96,Lam94,suzuki93,suzuki96,tian00}.

There are also hybrid algorithms that use a combination of constraint-based and scoring-based methods: In several works \cite{Singh93,SinVal95,spirtesmeek95,dashdruzdzel99,campospuertaIJIS02} the independence-based and scoring-based algorithms are maintained as separate processes, which are combined in some way, whereas the hybridization proposed by  \citeA{pkdd,ijar01} is based on the development of a scoring function that quantifies the discrepancies between the independences displayed by the candidate network and the database, and the search process is limited by the results of some independence tests. 

In this paper, we focus on the scoring+search approach. Although algorithms in this category have commonly used local search methods \cite{Buntine91,CooHersk:k2,Chigeihec95,campospuertaIJIS02,heckgeigerchick}, due to the exponentially large size of the search space, there is a growing interest in other heuristic search methods, i.e. simulated annealing \cite{Chigeihec95}, tabu search \cite{Bouckaert95,muntenau00}, branch and bound \cite{tian00}, genetic algorithms and evolutionary programming \cite{larra:pami,Myers99,Wong99}, Markov chain Monte Carlo \cite{kockacastelo01,Myers99}, variable neighborhood search \cite{campospuertaISAS01,Puertatesis}, ant colony optimization \cite{ijar02,Puertatesis}, greedy randomized adaptive search procedures \cite{campospuertaIberamia02}, and estimation of distribution algorithms \cite{Blanco02}.

All of these employ different search methods but the same search space: the space of DAGs. A possible alternative is the space of the orderings of the variables \cite{merida,ipmu00,campospuertaESQ01,friedmankolleruai00,larra_orden}. In this paper, however, we are more interested in the space of equivalence classes of DAGs \cite{jverma90}, i.e. classes of DAGs with each representing a different set of probability distributions. There is also a number of learning algorithms that carry out the search in this space \cite{andersson97,chick:le,dashdruzdzel99,madigan96,spirtesmeek95}. This feature reduces the size of the search space, although recent results \cite{Gillispie01} confirm that this reduction is not as important in terms of the DAG space as previously hoped (the ratio of the number of DAGs to the number of equivalence classes is lower than four). The price we have to pay for this reduction is that the evaluation of the candidate structures does not take advantage of an important property of many scoring functions, namely decomposability, and therefore the corresponding algorithms are less efficient.

In this paper we propose a new search space which is closely related to the space of equivalence classes of DAGs, and which we have called the space of {\it restricted acyclic partially directed graphs} (RPDAGs). We define a local search algorithm in this space, and show that by using a decomposable scoring function, we can evaluate locally the score of the structures in the neighborhood of the current RPDAG, thus obtaining an efficient algorithm while retaining many of the advantages of using equivalence classes of DAGs. After the original submission of this paper, \citeA{Chickering02} proposed another learning algorithm that searches in the space of equivalence classes of DAGs and which can also score the candidate structures locally, using a canonical representation scheme for equivalence classes, called {\it completed acyclic partially directed graphs} (CPDAGs).


The rest of the paper is organized as follows: section \ref{classes} discusses some preliminaries and the advantages and disadvantages of carrying out the search process in the spaces of DAGs and equivalence classes of DAGs. Section \ref{RPDAGs} describes the graphical objects, RPDAGs, that will be included in the proposed search space. In Section \ref{search}, a detailed description of the local search method used to explore this space is provided. Section \ref{decompo} shows how we can evaluate RPDAGs efficiently using a decomposable scoring function. Section \ref{resulta} contains the experimental results of the evaluation of the proposed algorithm on the Alarm \cite{alarm}, Insurance \cite{insurance97} and Hailfinder \cite{hailfinder} networks, as well as on databases from the UCI Machine Learning Repository. We also include an empirical comparison with other state-of-the-art learning algorithms. Finally, Section \ref{futuro} contains the concluding remarks and some proposals for future work.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{DAGs and Equivalence Classes of DAGs}
\label{classes}

The search procedures used within Bayesian network learning algorithms usually operate on the space of DAGs. In this context, the problem can be formally expressed as: Given a complete training set (i.e. we do not consider missing values or latent variables) $D=\{{\bf u}^1,\ldots , {\bf u}^m\}$ of instances of a finite set of $n$ variables, ${\mathcal{U}}$, find the DAG $H^*$ such that 
%
\begin{equation}
H^*=\arg\max_{H \in {\mathcal H}_n} g(H:D)
\end{equation}
%
where $g(H:D)$ is a scoring function measuring the fitness of any candidate DAG $H$ to the dataset $D$, and ${\mathcal H}_n$ is the family of all the DAGs with $n$ nodes\footnote{For reasons of simplicity, the set of nodes which are in one-to-one correspondence with the variables in ${\mathcal{U}}$ will also be denoted by ${\mathcal{U}}$.}.

Many of the search procedures, including the commonly used local search methods, rely on a neighborhood structure that defines the local rules (operators) used to move within the search space. The standard neighborhood in the space of DAGs uses the operators of arc addition, arc deletion and arc reversal, thereby avoiding (in the first and the third case) the inclusion of directed cycles in the graph. 

The algorithms that search in the space of DAGs using local methods are efficient mainly because of the decomposability property that many scoring functions exhibit. A scoring function $g$ is said to be {\it decomposable} if the score of any Bayesian network structure may be expressed as a product (or a sum in the log-space) of local scores involving only one node and its parents:
%
\begin{eqnarray}
g(H:D) = \sum_{y\in {\mathcal{U}}} g_D(y,Pa_H(y)) \\ \nonumber \\ 
g_D(y,Pa_H(y)) = g(y,Pa_H(y):N_{y,Pa_H(y)})
\label{decomp}
\end{eqnarray}
%
\noindent where $N_{y,Pa_H(y)}$ are the statistics of the variables $y$ and $Pa_H(y)$ in $D$, i.e. the number of instances in $D$ that match each possible instantiation of $y$ and $Pa_H(y)$. $Pa_H(y)$ will denote the parent set of $y$ in the DAG $H$, i.e. $Pa_H(y)=\{t\in {\cal U} \,|\, t \rightarrow y \in H\}$.

A procedure that changes one arc at each move can efficiently evaluate the improvement obtained by this change. Such a procedure can reuse the computations carried out at previous stages, and only the statistics corresponding to the variables whose parent sets have been modified need to be recomputed. The addition or deletion of an arc $x\rightarrow y$ in a DAG $H$ can therefore be evaluated by computing only one new local score, $g_D(y,Pa_{H\cup\{x\}}(y))$ or $g_D(y,Pa_{H\setminus\{x\}}(y))$, respectively. The evaluation of the reversal of an arc $x\rightarrow y$ requires the computation of two new local scores, $g_D(y,Pa_{H\setminus\{x\}}(y))$ and $g_D(x,Pa_{H\cup\{y\}}(x))$.


It should be noted that each structure in the DAG space is not always different from the others in terms of its representation capability: if we interpret the arcs in a DAG as causal interactions between variables, then each DAG represents a different model; however, if we see a DAG as a set of dependence/independence relationships between variables (that permits us to factorize a joint probability distribution), then different DAGs may represent the same model. Even in the case of using a causal interpretation, if we use observation-only data (as opposed to experimental data where some variables may be manipulated), it is quite common for two Bayesian networks to be empirically indistinguishable. When two DAGs $H$ and $H'$ can represent the same set of conditional independence assertions, we say that these DAGs are {\it equivalent}\footnote{Several authors also use the term {\it independence equivalent}, and reserve the term {\it distribution equivalent} (wrt some family of distributions ${\cal F}$) for the more restricted case where the two DAGs can represent the same probability distributions. In the common situation where all the variables are discrete and ${\cal F}$ is the family of unrestricted multinomial distributions, these two concepts of equivalence coincide.} \cite{jverma90}, and we denote this as $H\simeq H'$.

When learning Bayesian networks from data using scoring functions, two different (but equivalent) DAGs may be indistinguishable, due to the existence of invariant properties on equivalent DAGs, yielding equal scores. We could take advantage of this in order to get a more reduced space of structures to be explored.

The following theorem provides a graphical criterion for determining the equivalence of two DAGs:
%
\begin{theo}{\bf\cite{jverma90}}
\label{theo1}
Two DAGs are equivalent if and only if they have the same skeleton and the same v-structures.
\end{theo}
%
The {\it skeleton} of a DAG is the undirected graph that results from ignoring the directionality of every edge. A {\it v-structure} in a DAG $H$ is an ordered triplet of nodes, $(x,y,z)$, such that (1) $H$ contains the arcs $x\rightarrow y$ and $y \leftarrow z$, and (2) the nodes $x$ and $z$ are not adjacent in $H$. A {\it head-to-head pattern} (shortened h-h) in a DAG $H$ is an ordered triplet of nodes, $(x,y,z)$, such that $H$ contains the arcs $x\rightarrow y$ and $y \leftarrow z$. Note that in an h-h pattern $(x,y,z)$ the nodes $x$ and $z$ can be adjacent.

Another characterization of equivalent DAGs was presented by \citeA{chick:eq}, together with proof that several scoring functions used for learning Bayesian networks from data give the same score to equivalent structures (such functions are called {\it score equivalent functions}).

The concept of equivalence of DAGs partitions the space of DAGs into a set of equivalence classes. Whenever a score equivalent function is used, it seems natural to search for the best configuration in this new space of equivalence classes of DAGs. This change in the search space may bring several advantages:
\begin{itemize}
\item The space of equivalent classes is more reduced than the space of DAGs (although it is still enormous). We could therefore expect to obtain better results (with the same search effort).
%
\item As we do not spend time generating (using the operators defined to move between neighboring configurations in the search space) and evaluating (using the scoring function) equivalent DAGs, we could obtain more efficient algorithms. However, as the ratio of the number of equivalence classes to the number of DAGs seems (empirically) to asymptote to 0.267 \cite{Gillispie01}, the efficiency improvement may be small.
%
\item The search in the space of DAGs may be easily trapped in a local optimum, and the situation worsens as the operators defined for this space can move between configurations corresponding to equivalent DAGs (which will be evaluated with the same score). This difficulty can be partially avoided if we search in the space of equivalence classes.
\end{itemize}

The disadvantages are that, in this space of equivalence classes, it is more expensive to generate neighboring configurations, because we may be forced to perform some kind of consistency check, in order to ensure that these configurations represent equivalence classes\footnote{Note that the operators of addition and reversal of an arc in the DAG space also need a consistency check, but in this case we simply test the absence of directed cycles.}; in addition, the evaluation of the neighboring configurations may also be more expensive if we are not able to take advantage of the decomposability property of the scoring function. Finally, the new search space might introduce new local maxima that are not present in DAG space.


In order to design an exploring process for the space of equivalence classes we could use two distinct approaches: the first consists in considering that an equivalence class is represented by any of its components (in this case, it is necessary to avoid evaluating more than one component per class); and the second consists in using a canonical representation scheme for the classes.

The graphical objects commonly used to represent equivalence classes of DAGs are {\it acyclic partially directed graphs} \cite{jverma90} (known as {\it PDAGs}). These graphs contain both directed (arcs) and undirected (links) edges, but no directed cycles. Given a PDAG $G$ defined on a finite set of nodes ${\cal U}$ and a node $y\in {\cal U}$, the following subsets of nodes are defined:
%
\begin{itemize}
\item $Pa_G(y)=\{t\in {\cal U} \,|\, t \rightarrow y \in G\}$, the set of {\it parents} of $y$.
\item $Ch_G(y)=\{t\in {\cal U} \,|\, y \rightarrow t \in G\}$, the set of {\it children} of $y$. 
\item $Ne_G(y)=\{t\in {\cal U} \,|\, y\mbox{---} t \in G\}$, the set of {\it neighbors} of $y$.
\item $Ad_G(y)=\{t\in {\cal U} \,|\, t \rightarrow y \in G, \mbox{ or } y\rightarrow t \in G \mbox{ or } y\mbox{---} t \in G\}$, the set of {\it adjacents} to $y$. Obviously $Ad_G(y)=Pa_G(y)\cup Ch_G(y)\cup Ne_G(y)$.
\end{itemize}
%
An arc $x\rightarrow y$ in a DAG $H$ is {\it compelled} if it appears in every DAG belonging to the same equivalence class as $H$. An arc $x\rightarrow y$ in $H$ is said to be {\it reversible} if it is not compelled, i.e. there is a DAG $H'$ equivalent to $H$ that contains the arc $x\leftarrow y$. As every DAG in a particular equivalence class has the same set of compelled and reversible arcs, a canonical representation of an equivalence class is the PDAG consisting of an arc for every compelled arc in the equivalence class, and a link for every reversible arc. This kind of representation has been given several names: {\it pattern} \cite{spirtesmeek95}, {\it completed PDAG} ({\it CPDAG}) \cite{chick:le}, {\it essential graph} \cite{andersson97,dashdruzdzel99}. As a consequence of theorem \ref{theo1}, a completed PDAG possesses an arc $x\rightarrow y$ if and only if a triplet of nodes $(x,y,z)$ forms a v-structure or the arc $x\rightarrow y$ is required to be directed due to other v-structures (to avoid forming a new v-structure or creating a directed cycle) (see Figure \ref{ejemplo}). 

\begin{figure}[hbt]
\centerline{\psfig{figure=./figuras/midagypdag.eps,height=.3\textwidth}}
\caption{(a) Dag and (b) completed PDAG. The arcs $z\rightarrow x$, $y\rightarrow x$ and $x\rightarrow u$ are compelled; the arc $y\rightarrow w$ is reversible}
\label{ejemplo}
\end{figure}

Note that an arbitrary PDAG does not necessarily represent some equivalence class of DAGs, although there is a one-to-one correspondence between completed PDAGs and equivalence classes of DAGs. Nevertheless, completed PDAGs are considerably more complicated than general PDAGs. A characterization of the specific properties that a PDAG must verify in order to be a completed PDAG was obtained by \citeA{andersson97}:
%
\begin{theo}{\bf \cite{andersson97}}
\label{theo2}
A PDAG $G$ is a completed PDAG if and only if it satisfies the following conditions:
\begin{itemize}
\item[1] $G$ is a chain graph, i.e. it contains no (partially) directed cycles.
%
\item[2] The subgraph induced by every chain component\footnote{A chain component of $G$ is any connected component of the undirected graph obtained from $G$ by removing all the arcs.} of $G$ is chordal (i.e. on every undirected cycle of length greater than or equal to 4 there are two non-consecutive nodes connected by a link).
%
\item[3] The configuration $x\rightarrow y\mbox{---}z$ does not occur as an induced subgraph of $G$.
%
\item[4] Every arc $x\rightarrow y\in G$ occurs in at least one of the four configurations displayed in Figure~\ref{4conf} as an induced subgraph of $G$.
\end{itemize}
\end{theo}


\begin{figure}[hbt]
\centerline{\psfig{figure=./figuras/4conf.eps,width=.67\textheight}} %height=.3\textwidth}}
\caption{The four different configurations containing an arc $x\rightarrow y$ in a completed PDAG}
\label{4conf}
\end{figure}

Let us illustrate the advantages of searching in the space of equivalence classes of DAGs rather than the space of DAGs with a simple example. Figure \ref{pattern2} displays the set of possible DAGs involving three nodes $\{x,y,z\}$, with arcs between $z$ and $x$, and between $y$ and $x$. The first three DAGs are equivalent. In terms of independence information, they lead to the same independence statement $I(y,z|x)$ ($y$ and $z$ are conditionally independent given $x$), whereas the statement $I(y,z|\emptyset)$ ($y$ and $z$ are marginally independent) corresponds to the fourth one.  The four DAGs may be summarized by only two different completed PDAGs, shown in Figure~\ref{equival2}.

\begin{figure}[hbt]
\centerline{\psfig{figure=./figuras/mipattern2.eps,height=.3\textwidth}}
\caption{Four different DAGs with three nodes and two arcs}
\label{pattern2}
\end{figure}


\begin{figure}[hbt]
\centerline{\psfig{figure=./figuras/miequival2.eps,height=.15\textwidth}}
\caption{Two different equivalence classes of DAGs}
\label{equival2}
\end{figure}

As we can see, the search space may be reduced by using PDAGs to represent the classes: in our example, to two classes instead of four configurations; it can be seen \cite{andersson97} that the ratio of the number of DAGs to the number of classes is $25\,/\,11$ for three nodes, $543\,/\,185$ for four nodes and $29281\,/\,8792$ for five nodes; in more general terms, the results obtained by \citeA{Gillispie01} indicate that this ratio approaches a value of less than four as the number of nodes increases. The use of equivalence classes therefore entails convenient savings in exploration and evaluation effort, although the gain is not spectacular.

On the other hand, the use of a canonical representation scheme allows us to explore the space progressively and systematically, without losing any unexplored configuration unnecessarily. Returning to our example, let us suppose that the true model is the DAG displayed in Figure~\ref{equival2}.b and we start the search with an empty graph (with no arcs). Let us also assume that the search algorithm identifies that an edge between $x$ and $y$ produces the greatest improvement in the score. At this moment, the two alternatives, $x\rightarrow y$ and $x\leftarrow y$ (case 1 and case 2 in Figure \ref{twooptions}, respectively), are equivalent. Let us now suppose that we decide to connect the nodes $x$ and $z$; again we have two options: $x\rightarrow z$ or $x\leftarrow z$. Nevertheless, depending on the previous selected configuration, we obtain different outcomes that are no longer equivalent (see Figure \ref{twooptions}). 

If we had chosen case 1 (thus obtaining either case 1.1 or case 1.2), we would have eliminated the possibility of exploring the DAG $z\rightarrow x\leftarrow y$, and therefore the exploring process would have been trimmed. As the true model is precisely this DAG (case 2.1 in Figure \ref{twooptions}), then the search process would have to include another arc connecting $y$ and $z$ (cases 1.1.1, 1.2.1 or 1.2.2), because $y$ and $z$ are conditionally dependent given $x$. At this moment, any local search process would stop (in a local optimum), because every local change (arc reversal or arc removal) would make the score worse.

\begin{figure}[htb]

\centerline{\psfig{figure=./figuras/trapped.eps,height=.32\textwidth}}
\caption{Local search in the space of DAGs is trapped at a local optimum}
\label{twooptions}
\end{figure}



Consequently, our purpose consists in adding or removing edges (either links or arcs) to the structure without pruning the search space unnecessarily. We could therefore introduce links instead of arcs (when there is not enough information to distinguish between different patterns of arcs), which would serve as templates or dynamic linkers to equivalence patterns. They represent any valid combination of arcs which results in a DAG belonging to the same equivalence class.

Looking again at the previous example, we would proceed as follows: assuming that in our search space the operators of link addition and creation of h-h patterns are available, we would first include the link $x\mbox{---}y$; secondly, when considering the inclusion of a connection between $x$ and $z$, we would have two options, shown in Figure~\ref{equival2}: the h-h pattern $z \rightarrow x\leftarrow y$ and the pattern $z \mbox{---}x \mbox{---}y$. In this case the scoring function would assign the greatest score to the h-h pattern $z\rightarrow x\leftarrow y$, thus obtaining the correct DAG.


\section{Restricted Acyclic Partially Directed Graphs}
\label{RPDAGs}

The scheme of representation that we will use is slightly different from the formalism of completed PDAGs. It is not necessary for each configuration of our search space (which we call {\it restricted PDAG} or {\it RPDAG}) to correspond to a different equivalence class; two different RPDAGs may correspond to the same equivalence class. The main reason for this is efficiency: by allowing an equivalence class to be represented (only in some cases) by different RPDAGs, we will gain in efficiency to explore the space. Before explaining this in greater detail, let us define the concept of RPDAG: 

\begin{defin}[restricted PDAG] \label{defiRPDAG}
A PDAG $G$ is a restricted PDAG (RPDAG) if and only if it satisfies the following conditions:
\begin{itemize}
\item[1] $\forall y\in {\cal U}$, $Pa_G(y)\neq\emptyset \, \Rightarrow \, Ne_G(y)=\emptyset$.
\item[2] $G$ does not contain any directed cycle.
\item[3] $G$ does not contain any completely undirected cycle, i.e. a cycle containing only links.
\item[4] If $x\rightarrow y$ exists in $G$ then either $|Pa_G(y)|\geq 2$ or $Pa_G(x)\neq\emptyset$. 

This condition states that an arc $x\rightarrow y$ exists in $G$ only if it is either part of an h-h pattern or there is another arc (originated by an h-h pattern) going to $x$.
\end{itemize}
\end{defin}

As an RPDAG is a PDAG, it could be considered to be a representation of a set of (equivalent) DAGs. We therefore must define which the set of DAGs is represented by a given RPDAG $G$, i.e. how direction may be given to the links in $G$ in order to extend it to a DAG. The following definition formalizes this idea.

\begin{defin}[Extension of a PDAG] \label{defiextension}
Given a PDAG $G$, we say that a DAG $H$ is an extension of $G$ if and only if:
\begin{itemize}
\item[1] $G$ and $H$ have the same skeleton.
\item[2] If $x\rightarrow y$ is an arc in $G$ then $x\rightarrow y$ is also an arc in $H$ (no arc is redirected).
\item[3] $G$ and $H$ have the same h-h patterns (i.e. the process of directing the links in $G$ in order to produce $H$ does not generate new h-h patterns).
\end{itemize}
\end{defin}

We will use $Ext(G)$ to denote the set of DAGs that are extensions of a given PDAG $G$.

\begin{prop} \label{prop1}
Let $G$ be an RPDAG. Then:
\begin{itemize}
\item[(a)] $Ext(G)\neq\emptyset$ ($G$ can be extended to obtain a DAG, i.e. the extension of an RPDAG is well-defined).
\item[(b)] $\forall H,H'\in Ext(G) \;\; H\simeq H'$ (i.e. all the different DAGs that can be obtained from $G$ by extending it are equivalent).
\end{itemize}
\end{prop}

\noindent {\bf Proof:} \\
(a) As $G$ has no directed cycle (condition 2 in Definition~\ref{defiRPDAG}), then either $G$ is already a DAG or it has some links. Let us consider an arbitrary link $x\mbox{---}y$. Using condition 1 in Definition~\ref{defiRPDAG}, neither $x$ nor $y$ can have a parent. We can then direct the link $x\mbox{---}y$ in either direction without creating an h-h pattern. If we direct the link $x\mbox{---}y$ as $x\rightarrow y$ and $y$ is part of another link $y\mbox{---}z$, then we direct it as $y\rightarrow z$ (in order to avoid a new h-h pattern). We can continue directing the links in a chain in this way, and this process cannot generate a directed cycle because, according to condition 3 in Definition~\ref{defiRPDAG}, $G$ has no completely undirected cycle.\\
(b) The extension process of $G$ does not modify the skeleton and does not create new h-h patterns. Therefore, all the extensions of $G$ have the same skeleton and the same v-structures (a v-structure is a particular case of h-h pattern), hence they are equivalent. \qed

\vspace{0.5cm}
It should be noted that condition 4 in Definition \ref{defiRPDAG} is not necessary to prove the results in Proposition \ref{prop1}. This condition is included to ensure that the type of PDAG used to represent subsets of equivalent DAGs is as general as possible. In other words, condition 4 guarantees that an RPDAG is a representation of the greatest number of equivalent DAGs, subject to the restrictions imposed by conditions 1-3 in Definition \ref{defiRPDAG}. As we will see in the next proposition, this is achieved by directing the minimum number of edges. For example, $z\mbox{---}x\rightarrow y\rightarrow u$ would not be a valid RPDAG. The RPDAG that we would use in this case is $z\mbox{---}x\mbox{---}y\mbox{---}u$.

\begin{prop} \label{prop2}
Let $G$ be a PDAG verifying the conditions 1-3 in Definition \ref{defiRPDAG}. There is then a single RPDAG $R$ such that $Ext(G)\subseteq Ext(R)$.
\end{prop}

\noindent {\bf Proof:} \\
The proof is constructive. We shall build the RPDAG $R$ as follows: the skeleton and the h-h patterns of $R$ are the same as those in $G$. An arc $x\rightarrow y$ in $G$ shall now be considered such that $Pa_G(x)=\emptyset$ and $Pa_G(y)=\{x\}$ (if such an arc does not exist, then $G$ itself would be an RPDAG): we convert the arc $x\rightarrow y$ into the link $x\mbox{---}y$. This process is then repeated. Obviously, the PDAG $R$ obtained in this way has no directed cycle and verifies condition 4 in Definition \ref{defiRPDAG}. Moreover, we cannot obtain a configuration $z\rightarrow x\mbox{---}y$ as a subgraph of $R$ because $Pa_G(x)=\emptyset$ (we only remove the direction of arcs whose initial nodes have no parent). In addition, $R$ cannot have any completely undirected cycle because either the arc $x\rightarrow y$ is not part of any cycle in $G$ or it is part of a cycle in $G$ that must contain at least one h-h pattern (and the directions of the arcs in this pattern will never be removed). $R$ is therefore an RPDAG. 

Let us now prove that $Ext(G)\subseteq Ext(R)$: if $H\in Ext(G)$ then $H$ and $G$ have the same skeleton and h-h patterns, hence $H$ and $R$ also have the same skeleton and h-h patterns. Moreover, as all the arcs in $R$ are also arcs in $G$, if $x\rightarrow y\in R$ then $x\rightarrow y\in G$, which in turn implies that $x\rightarrow y\in H$. Therefore, according to Definition~\ref{defiextension}, $H\in Ext(G)$.

Finally, let us prove the uniqueness of $R$: we already know that any other RPDAG $R^{\prime}$ verifying that $Ext(G)\subseteq Ext(R^{\prime})$ has the same skeleton and h-h patterns as $R$. According to condition 1 in Definition \ref{defiRPDAG}, the edges that are not part of any of these h-h patterns but are incident to the middle node $y$ in any h-h pattern $x\rightarrow y\leftarrow z$ must be directed away from $y$ (in order to avoid new h-h patterns). The remaining edges that are not part of any h-h pattern must be undirected, in order to satisfy condition 4 in Definition \ref{defiRPDAG}. There is therefore only one RPDAG that matches a given skeleton and a set of h-h patterns, so $R$ is the only RPDAG verifying that $Ext(G)\subseteq Ext(R)$. Figure~\ref{ejprop2} shows an example of the construction process. \qed


\begin{figure}[hbt]
\centerline{\psfig{figure=./figuras/miejprop2bis.eps,height=.32\textwidth}}
\caption{Illustrating the construction process in Proposition~\ref{prop2}: (a) PDAG $G$; (b) undirecting the arc $a\rightarrow y$; (c) undirecting the arc $y\rightarrow x$; (d) undirecting the arc $x\rightarrow e$, thus obtaining the RPDAG $R$}
\label{ejprop2}
\end{figure}

The following proposition ensures that the concept of RPDAG allows us to define a partition in the space of DAGs.

\begin{prop} \label{propnu}
Let $G$ and $G^{\prime}$ be two different RPDAGs. Then $Ext(G)\cap Ext(G^{\prime})=\emptyset$.
\end{prop}

\noindent {\bf Proof:} \\
Let $H$ be any DAG. Then $H$ itself is a PDAG and obviously $H=Ext(H)$. By applying the result in Proposition~\ref{prop2}, we can assert that there is a single RPDAG $G$ such that $H\subseteq Ext(G)$. \qed

\vspace{0.5cm}
In the proposition below, we show the properties which are common to all the DAGs belonging to the same extension of an RPDAG.

\begin{prop} \label{propc}
Two DAGs belong to the extension of the same RPDAG if and only if they have the same skeleton and the same h-h patterns.
\end{prop}

\noindent {\bf Proof:} \\
The necessary condition is obvious. Let us prove the sufficient condition: let $H$ and $H^{\prime}$ be two DAGs with common skeleton and h-h patterns. We shall construct a PDAG $G$ as follows: the skeleton and the h-h patterns of $G$ are the same as those in $H$ and $H^{\prime}$; the edges that have the same orientation in $H$ and $H^{\prime}$ are directed in $G$ in the same way; the other edges in $G$ remain undirected. From Definition~\ref{defiextension}, it is clear that $H,H^{\prime}\in Ext(G)$.

$G$ has no directed cycles because $H$ and $H^{\prime}$ are DAGs. $G$ has no completely undirected cycles, since all the cycles in $H$ and $H^{\prime}$ share at least the h-h patterns. In addition, $x\rightarrow y\mbox{---}z$ cannot be a subgraph of $G$ because this would imply the existence of the subgraphs $x\rightarrow y\leftarrow z$ and $x\rightarrow y\rightarrow z$ in $H$ and $H^{\prime}$, respectively, and therefore these two DAGs would not have the same h-h patterns. 

Therefore, the PDAG $G$ satisfies conditions 1-3 in Definition~\ref{defiRPDAG}. By applying Proposition~\ref{prop2}, we can then build a single RPDAG $R$ such that $Ext(G)\subseteq Ext(R)$, hence $H,H^{\prime}\in Ext(R)$. \qed

\vspace{0.5cm}
A characterization of the extension of an RPDAG that will be useful later is:
%
\begin{prop} \label{prop3}
Given an RPDAG $G$ and a DAG $H$, then $H$ is an extension of $G$ if and only if the following conditions hold:
\begin{itemize}
\item[1] $G$ and $H$ have the same skeleton.
\item[2] $\forall y\in {\mathcal{U}}$, if $Pa_G(y)\neq\emptyset$ then $Pa_H(y)=Pa_G(y)$.
\item[3] $\forall y\in {\mathcal{U}}$, if $Pa_G(y)=\emptyset$ and $Pa_H(y)\neq\emptyset$ then $|Pa_H(y)|=1$.
\end{itemize}
\end{prop}

\noindent {\bf Proof:} \\
\noindent $\bullet$ {\it Necessary condition}:

\noindent -- $Pa_G(y)\neq\emptyset$: Let $x\in Pa_G(y)$, i.e., $x\rightarrow y\in G$. Then, from condition 2 in Definition~\ref{defiextension}, $x\rightarrow y\in H$, i.e., $x\in Pa_H(y)$. Moreover, $z\in Pa_G(y) \Leftrightarrow$ $x\rightarrow y\leftarrow z$ is an h-h pattern in $G$. From condition 3 in Definition~\ref{defiextension}, this occurs if and only if $x\rightarrow y\leftarrow z$ is an h-h pattern in $H$, which is equivalent to $z\in Pa_H(y)$. Therefore, $Pa_H(y)=Pa_G(y)$.

\noindent -- $Pa_G(y)=\emptyset$ and $Pa_H(y)\neq\emptyset$: Let $x\in Pa_H(y)$. If there is another node $z\in Pa_H(y)$, then $x\rightarrow y\leftarrow z$ is an h-h pattern in $H$ and therefore it is also an h-h pattern in $G$, which contradicts the fact that $Pa_G(y)=\emptyset$. So, $y$ cannot have more than one parent in $H$, hence $|Pa_H(y)|=1$.

\vspace{0.25cm}
\noindent $\bullet$ {\it Sufficient condition}:

\noindent -- If $x\rightarrow y\in G$ then $Pa_G(y)\neq\emptyset$. From condition 2 we have $Pa_H(y)=Pa_G(y)$, hence $x\rightarrow y\in H$.

\noindent -- If $x\rightarrow y\leftarrow z$ is an h-h pattern in $G$, once again from condition 2, $Pa_H(y)=Pa_G(y)$ and therefore $x\rightarrow y\leftarrow z$ is an h-h pattern in $H$.

\noindent -- If $x\rightarrow y\leftarrow z$ is an h-h pattern in $H$, then $|Pa_H(y)|\neq1$ and $Pa_H(y)\neq\emptyset$. So, from condition 3, we obtain $Pa_G(y)\neq\emptyset$ and, from condition 2, $Pa_G(y)=Pa_H(y)$. Therefore $x\rightarrow y\leftarrow z$ is an h-h pattern in $G$. \ \ \qed

\subsection{Restricted PDAGs and Completed PDAGs}
%\vspace{0.5cm}
Let us now examine the main differences between the different representations: a representation based on PDAGs ensures that every equivalence class has a unique representation, but there are PDAGs that do not correspond to any equivalence class (in other words, the mapping from equivalence classes to PDAGs is injective). On the other hand, our representation based on RPDAGs guarantees that every RPDAG corresponds to an equivalence class (proposition \ref{prop1}) but does not ensure that every equivalence class has a single representation (the mapping from equivalence classes to RPDAGs is onto). However, the mapping from equivalence classes to CPDAGs is bijective. Figures \ref{figdiferencias}.a, \ref{figdiferencias}.b and \ref{figdiferencias}.c show the three RPDAGs corresponding to the same equivalence class; the associated completed PDAG is shown in Figure \ref{figdiferencias}.d. In this example, the number of DAGs in the equivalence class is 12.

\begin{figure}[htb]
\centerline{\psfig{figure=./figuras/midiferenciasbis.eps,height=.26\textwidth}}
\caption{(a), (b) and (c) Three RPDAGs corresponding to the same equivalence class; (d) the associated completed PDAG}
\label{figdiferencias}
\end{figure}

As we can see, the difference appears when there are triangular structures. If we compare the definition of RPDAG (Definition \ref{defiRPDAG}) with the characterization of CPDAGs (Theorem~\ref{theo2}), we may observe that the essential difference is that a CPDAG may have completely undirected cycles, but these cycles must be {\em chordal}. In RPDAGs, undirected cycles are therefore forbidden, whereas in CPDAGs undirected non-chordal cycles are forbidden.

It should be noted that we could also define RPDAGs by replacing condition 3 in Definition \ref{defiRPDAG} for its equivalent: {\it The subgraph induced by every chain component of $G$ is a tree} (which is a specific type of chordal graph). In this way, the similarities and differences between CPDAGs and RPDAGs are even clearer. Any of the RPDAGs in the same equivalence class is obtained from the corresponding CPDAG by removing some of the links (converting them into arcs) in order to obtain a tree structure.

Examining the problem from another perspective, from Theorem \ref{theo1} and Proposition \ref{propc} we can see that the role played by the v-structures in CPDAGs is the same as that played by the h-h patterns in RPDAGs.

It is also interesting to note that the number of DAGs which are extensions of a given RPDAG, $G$, can be calculated very easily: the subgraph induced by each chain component of $G$ is a tree, and this tree can be directed in different ways by selecting any of the nodes as the root node. Moreover, we can proceed independently within each chain component. The number of DAGs in $Ext(G)$ is therefore equal to the product of the number of nodes within each chain component of $G$. Regarding the number of RPDAGs that represent the same equivalence class, this number grows exponentially with the size of the undirected cliques in the CPDAG. For example, if the subgraph induced by a chain component in a CPDAG consists of a complete subgraph of $m$ nodes, then the number of RPDAG representations is $\frac{m!}{2}$. This obviously does not mean that a search method based on RPDAGs must explore all these equivalent representations. 

Our reason for using RPDAGs is almost exclusively practical. In fact, RPDAGs do not have a clear semantics (they are a somewhat hybrid creature, between DAGs and completed PDAGs). We can only say that RPDAGs would correspond to sets of equivalent DAGs which share all the causal patterns where an effect node has at least two causes (and only the causal patterns where a single cause provokes an effect are not determined). This is not problematic when we are performing model {\it selection} but it becomes critical if we are doing model {\it averaging}: without a semantic understanding of the class of RPDAGs, it will be quite difficult to assign a prior to them.

Our intention is to trade the uniqueness of the representation of equivalence classes of DAGs (CPDAGs) for a more manageable one (RPDAGs): testing whether a given PDAG $G$ is an RPDAG is easier than testing whether $G$ is a completed PDAG. In the first case, the consistency check involves testing for the absence of directed and completely undirected cycles (the complexity of these tests and those necessary to verify whether a directed graph is a DAG is exactly the same), whereas in the second case, in addition to testing for the absence of directed and partially directed cycles, we also need to perform chordality tests and check that each arc is part of one of the induced subgraphs displayed in Figure~\ref{4conf}. The price we have to pay for using RPDAGs is that we may occasionally need to evaluate an equivalence class more than once. In the next section, we will examine how a local search method which uses RPDAGs can also take advantage of the decomposability property of a scoring function in order to efficiently evaluate neighboring structures.



\section{The Search Method}
\label{search}

We will use a local method to explore the search space of RPDAGs. The starting point of the search process will be an empty RPDAG (corresponding to an empty DAG). Nevertheless, we could start from another configuration if we have some prior knowledge about the presence or absence of some edges or v-structures. We must then define the operators to move from one configuration to another neighboring configuration.   

\subsection{Overview}

Our basic operators are the inclusion of an edge between a pair of non-adjacent nodes and the removal of an existing edge between a pair of adjacent nodes in the current configuration. These edges may be either directed or undirected. 

The inclusion of an isolated link $x\mbox{---}y$ will serve as a template for the arcs $x\rightarrow y$ and $x\leftarrow y$; however, the link $x\mbox{---}y$ together with another link $x\mbox{---}z$ represent any combination of arcs except those that create new h-h patterns (the DAGs (a), (b) and (c) in Figure~\ref{pattern2}). 
In the case of adding an arc, we may obtain several different neighboring configurations, depending on the topology of the current RPDAG and the direction of the arc being included. As we will see, if we are testing the inclusion of an edge between two nodes $x$ and $y$, this may involve testing some of the different valid configurations obtained by the inclusion of the link $x\mbox{---}y$, the arc $x\rightarrow y$, the arc $x\leftarrow y$, the h-h pattern $x\rightarrow y\leftarrow z$ or the h-h pattern $z\rightarrow x\leftarrow y$ (where $z$ in the last two cases would be any node such that either the link $y\mbox{---}z$ or the link $z\mbox{---}x$ exists in the current configuration). However, the removal of an edge will always result in only one neighboring configuration.
Other operators, such as arc reversal \cite{chick:le}, will not be used by our search method.
The set of neighboring configurations of a given RPDAG $G$ will therefore be the set of all the different RPDAGs obtained from $G$ by adding or deleting a single edge (either directed or undirected).


Before explaining the details of the search method, let us illustrate the main ideas by means of the following example: consider the RPDAG in Figure~\ref{config1}, which represents the current configuration of the search process (this figure only displays the part of the RPDAG corresponding to the neighborhood of the nodes $x$ and $y$), and assume that we shall include an edge between the nodes $x$ and $y$. 

\begin{figure}[hbt]
\centerline{\psfig{figure=./figuras/miconfig1f.eps,height=.18\textwidth}}
\caption{Node $x$ has one parent and one child, $y$ does not have any parents and has two neighbors}
\label{config1}
\end{figure}

In this situation, we cannot introduce the link $x\mbox{---}y$ because we would violate one of the conditions defining RPDAGs (condition 1). We may introduce the arc $x\rightarrow y$ and in this case, again in order to preserve condition 1, the two neighbors of $y$ must be converted into children. We can also include the arc $x\leftarrow y$. Finally, we may include two different h-h patterns $x\rightarrow y\leftarrow z$, where $z$ is a neighbor of $y$ (the other neighbor must be converted into a child, once again in order to preserve condition 1). These four different configurations are displayed in Figure~\ref{config2}. 

\begin{figure}[hbt]
\centerline{\psfig{figure=./figuras/miconfig2f.eps,height=.4\textwidth}}
\caption{Neighboring configurations obtained by including an edge between $x$ and $y$ in the configuration of Fig.~\ref{config1}}
\label{config2}
\end{figure}

\subsection{Neighboring Configurations}

In order to design a systematic way to determine which neighboring RPDAGs arise from the inclusion or the removal of an edge in an RPDAG, it is sufficient to consider some local parameters of the two nodes to be connected. First, some additional notation is introduced. If $|\cdot|$ represents the cardinality of a set, given a node $x$ in a PDAG $G$, we define:



\begin{itemize}
\item $p_G(x) = |Pa_G(x)|$, $ \ \ \bullet \ \ $ $c_G(x) = |Ch_G(x)|$
\item $n_G(x) = |Ne_G(x)|$, $ \ \ \bullet \ \ $ $a_G(x) = |Ad_G(x)|$
\end{itemize}
%
Observe that for any RPDAG $G$, the following two properties hold:
\begin{itemize}
\item $p_G(x)+c_G(x)+n_G(x) = a_G(x)$
\item if $p_G(x) \neq 0 \Rightarrow n_G(x)=0$ (hence $p_G(x)+c_G(x) = a_G(x)$)
\end{itemize}

\subsubsection{Adding Edges}

The number and type of neighboring configurations that can be obtained from the inclusion in the current RPDAG of an edge between $x$ and $y$ can be determined from the parameters above. The resultant casuistry may be reduced to seven states, which we have labeled from $A$ to $G$. In order to facilitate its description, we shall use the decision tree shown in Figure~\ref{arb_stats}\footnote{This tree could be organized differently in order to improve the efficiency in the location of the current state. However, this particular tree was selected so as to clarify the presentation.}.

\begin{figure}[htb]
\centerline{\psfig{figure=./figuras/miarbolestados.eps,height=1.1\textwidth}}
%height=1.2\textwidth}}
\caption{The tree of possible states that may result by adding an edge between nodes $x$ and $y$}
\label{arb_stats}
\end{figure}

In this tree, the lower box of each non-terminal vertex contains a test (about the number of parents or the number of neighbors of nodes $x$ and $y$). The lower box of each terminal vertex contains the label of the state resulting from following the path from the root to that terminal vertex. The description of each state (i.e. the different neighboring configurations that can be obtained in this case) can be found in Table~\ref{table_stats}. The upper boxes of all the vertices in the tree show the restrictions imposed on each intermediate or terminal state. For example, state B corresponds to a situation where both nodes $x$ and $y$ do not have neighbors and at least one of them has some parent. Although the tree has seven different states, there are only five truly different states, since states D and E, and states F and G are symmetrical.

\begin{table}[htb]
\begin{center}
\begin{tabular}{|c|c|l|c|c|c|}
\hline\hline
State  & Number of      & Added edges & Directed & Undirected & Completing ? \\
       & configurations &  & Cycles ? & Cycles ?   &              \\
\hline\hline
A & 1 & \begin{tabular}{c} 
               $x$---$y$ \\ 
        \end{tabular}  & No & No & No \\ \hline
B & 2 & \begin{tabular}{l}
               $x\rightarrow y$ \\
	       $x\leftarrow y$ \\
        \end{tabular}  & 
                   \begin{tabular}{c} Yes$^{(1)}$ \\ Yes$^{(2)}$ \\ \end{tabular} & 
                   \begin{tabular}{c} No \\ No \\ \end{tabular} &
                   \begin{tabular}{c} No \\ No \\ \end{tabular} \\ \hline

C & $n_G(x)+n_G(y)+1$ & \begin{tabular}{l} 
                               $x$---$y$ \\ 
                               $x\rightarrow y\leftarrow z$ \\
                               $z\rightarrow x\leftarrow y$ \\
                         \end{tabular} &
         \begin{tabular}{c} No \\ Yes$^{(3)}$ \\ Yes$^{(4)}$ \\ \end{tabular} & 
         \begin{tabular}{c} Yes \\ No \\ No \\ \end{tabular} &
         \begin{tabular}{c} No \\ Yes$^{(3)}$ \\ Yes$^{(4)}$ \\ \end{tabular} \\ \hline
D & $n_G(y)+1$ & \begin{tabular}{l} 
                        $x$---$y$ \\ 
                        $x\rightarrow y\leftarrow z$ \\
                 \end{tabular} &
               \begin{tabular}{c} No \\ No \\ \end{tabular} & 
               \begin{tabular}{c} No \\ No \\ \end{tabular} &
               \begin{tabular}{c} No \\ Yes$^{(3)}$ \\ \end{tabular} \\ \hline
E & $n_G(x)+1$ & \begin{tabular}{l} 
                        $x$---$y$ \\ 
                        $z\rightarrow x\leftarrow y$ \\
                 \end{tabular} &
                   \begin{tabular}{c} No \\ No \\ \end{tabular} & 
                   \begin{tabular}{c} No \\ No \\ \end{tabular} &
                   \begin{tabular}{c} No \\ Yes$^{(4)}$ \\ \end{tabular} \\ \hline
F & $n_G(y)+2$ & \begin{tabular}{l}
                        $x\leftarrow y$ \\
                        $x\rightarrow y$ \\
	                $x\rightarrow y\leftarrow z$ \\
                 \end{tabular} &
               \begin{tabular}{c} No \\ Yes \\ Yes$^{(3)}$ \\ \end{tabular} & 
               \begin{tabular}{c} No \\ No \\ No \\ \end{tabular} &
               \begin{tabular}{c} No \\ Yes \\ Yes$^{(3)}$ \\ \end{tabular} \\ \hline
G & $n_G(x)+2$ & \begin{tabular}{l}
                        $x\rightarrow y$ \\
                        $x\leftarrow y$ \\
	                $z\rightarrow x\leftarrow y$ \\
                 \end{tabular}  &
                 \begin{tabular}{c} No \\ Yes \\ Yes$^{(4)}$ \\ \end{tabular} & 
                 \begin{tabular}{c} No \\ No \\ No \\ \end{tabular} &
                 \begin{tabular}{c} No \\ Yes \\ Yes$^{(4)}$ \\ \end{tabular} \\
\hline\hline
\multicolumn{3}{|l}{(1) only if $p_G(x)\neq 0$ and $c_G(y)\neq 0$} & 
\multicolumn{3}{l|}{(3) only if $n_G(y)\geq 2$} \\
\multicolumn{3}{|l}{(2) only if $p_G(y)\neq 0$ and $c_G(x)\neq 0$} &
\multicolumn{3}{l|}{(4) only if $n_G(x)\geq 2$} \\
\hline\hline
\end{tabular}
\end{center}
\caption{Table of states that may result by adding an edge between nodes $x$ and $y$}
\label{table_stats}
\end{table}

In Table~\ref{table_stats}, each row corresponds to a state: the first column contains the labels of the states; the second column displays the total number of neighboring configurations that can be obtained for each state; the third column shows the different types of edges that, for each state, can be added to the current configuration; columns four, five and six will be discussed later.

Using the example in Figure~\ref{config1}, we shall explain the use of the decision tree as well as the instantiation of the information in Table~\ref{table_stats}. Following the decision tree, at level 1 (the root vertex), the test is false since $y$ has two neighbors. At level 2, the test is also false as $x$ has one parent. At level 3 the test is true, since $x$ has no neighbor. At level 4 we reach a terminal vertex. Our current configuration therefore corresponds to state F. Then, by examining state F in Table~\ref{table_stats}, we can confirm that we reach four different configurations ($n_G(y)=2$): $G\cup \{x \leftarrow y\}$, $G\cup \{x \rightarrow y\}$ without new h-h patterns, and two $G\cup \{x \rightarrow y\leftarrow z\}$ which produce new h-h patterns. So, these are the only structures that our algorithm must evaluate when considering the inclusion of the candidate edge $x\mbox{---}y$. In Figure~\ref{configS}, we show an example for each of the five non-symmetrical states (once again, these examples only display the part of the RPDAGs corresponding to the neighborhood of the nodes $x$ and $y$).


We therefore have a systematic way to explore all the neighboring configurations which result from adding an edge. However, it will sometimes be necessary to perform some additional steps since the configurations obtained must be RPDAGs:
%
\begin{itemize}
\item First, we must maintain the condition~1 ($p_G(y)\neq 0\Rightarrow n_G(y)=0$). It is therefore necessary to {\it complete} the configuration for some of the described states, i.e. some of the links must be converted into arcs. The completing process consists in firing an orientation in cascade, starting from the links $y\mbox{---}t$ such that the arc just introduced is $x\rightarrow y$. Let us consider the situation in Figure~\ref{config3}, where we want to connect the nodes $x$ and $y$, a case corresponding to state D. Among the three possible neighboring configurations, let us suppose that we are testing the one which introduces the h-h pattern $x\rightarrow y\leftarrow z$. The RPDAG obtained from the completing process is also displayed in Figure~\ref{config3}. The sixth column in Table~\ref{table_stats} shows which states and neighboring configurations may require the completing process.
\end{itemize}

\begin{figure}[hbt]
\centerline{\psfig{figure=./figuras/miconfig3.eps,height=.25\textwidth}}
\caption{Transformation of a configuration after including the pattern $x\rightarrow y\leftarrow z$ and completing}\label{config3}
\end{figure}

\begin{itemize}
\item Secondly, it is possible that some of the neighboring configurations must be rejected, as they give rise to directed or completely undirected cycles (conditions 2 and 3 defining RPDAGs). For example, let us consider the situation displayed in Figure~\ref{config4}, which corresponds to state F. In this case, the configuration obtained after including the arc $x\rightarrow y$ and completing would generate a directed cycle. This configuration must therefore be rejected. The columns four and five in Table~\ref{table_stats} show which states and configurations may require a detection of directed or completely undirected cycles, respectively.
\end{itemize}

\begin{figure}[hbt]
\centerline{\psfig{figure=./figuras/miconfig4.eps,height=.3\textwidth}}
\caption{Neighboring configuration that gives rise to a directed cycle}
\label{config4}
\end{figure}

\subsubsection{Deleting Edges}

The other basic operator, the removal of an edge (either link or arc) is much simpler than the addition of an edge, since only one neighboring configuration is obtained when we delete an edge. Moreover, it is not necessary to perform any test for detecting directed or undirected cycles. However, in this case we need to preserve condition 4 in the definition of RPDAG (if $x\rightarrow y$ exists in $G \Rightarrow |Pa_G(y)|\geq 2$ or $Pa_G(x)\neq\emptyset$), that could be violated after an arc is deleted. This situation may appear only when we are going to remove an arc $x\rightarrow y$ and either $Pa_G(y)=\{x\}$ or $Pa_G(y)=\{x,u\}$: in the first case, all the children of $y$ that do not have other parents than $y$ must be converted into neighbors of $y$, and this process is repeated starting from each of these children; in the second case, if $Pa_G(u)=\emptyset$ then in addition to the previous process, the arc $u\rightarrow y$ must be converted into the link $u\mbox{---}y$. Figure~\ref{desorientacion} illustrates these situations. This procedure of transforming arcs into links is exactly the same as the one described in the proof of Proposition~\ref{prop2}.

\begin{figure}[hbt]
\centerline{\psfig{figure=./figuras/desorientacion.eps,height=.25\textwidth}}
\caption{Transforming arcs into links after removing the arc $x\rightarrow y$}
\label{desorientacion}
\end{figure}

\begin{figure}[htbp]
\centerline{\psfig{figure=./figuras/miejemplo.eps,width=.8\textwidth}}
\caption{For each state of the decision tree in Fig.~\ref{arb_stats}, an example of the neighboring configurations of an RPDAG that can be obtained after adding an edge between nodes $x$ and $y$}
\label{configS}
\end{figure}

\subsection{The Operators}
\label{operators}

Although the previous description of the operators that define the neighborhood of an RPDAG is quite convenient from a practical (implementation) point of view, for the sake of clarity, we shall describe them in another way. In fact, we shall use five operators: 
%
\begin{itemize}
\item $A\_arc(x,y)$, addition of an arc $x\rightarrow y$. 
\item $A\_link(x,y)$, addition of a link $x\mbox{---}y$. 
\item $D\_arc(x,y)$, deletion of an arc $x\rightarrow y$. 
\item $D\_link(x,y)$, deletion of a link $x\mbox{---}y$. 
\item $A\_hh(x,y,z)$, addition of an arc $x\rightarrow y$ and creation of the h-h pattern $x\rightarrow y\leftarrow z$ by transforming the link $y\mbox{---}z$ into the arc $y\leftarrow z$. 
\end{itemize}

The conditions that the current RPDAG $G$ must verify so that each of these operators can be applied in order to obtain a valid neighboring RPDAG are shown in Table~\ref{condiciones}. These conditions can be easily derived from the information in Figure~\ref{arb_stats} and Table~\ref{table_stats}. In Table~\ref{condiciones}, $UC(x\mbox{---}y)$ represents a test for detecting completely undirected cycles after inserting the link $x\mbox{---}y$ in the current RPDAG. Note that we can perform this test very easily without actually inserting the link: it is only necessary to check the existence of a path between $x$ and $y$ exclusively formed by the links. Similarly, $DC(x\rightarrow y)$ and $DC(x\rightarrow y\leftarrow z)$ represent tests for detecting directed cycles after inserting the arc $x\rightarrow y$ and the h-h pattern $x\rightarrow y\leftarrow z$ in the current RPDAG, respectively (and perhaps completing). It should also be noted that we can perform these tests without inserting the arc or the h-h pattern: in this case we only need to check the existence of a path from $y$ to $x$ containing only either links or arcs directed away from $y$ (a partially directed path from $y$ to $x$). Table~\ref{condiciones} also shows which operators may require a post-processing step in order to ensure that the corresponding neighboring configuration of $G$ is an RPDAG. In Table~\ref{condiciones}, $Complete(y)$ and $Undo(y)$ refer to the procedures that preserve conditions 1 and 4 in Definition~\ref{defiRPDAG}, respectively. Note that both $UC(x\mbox{---}y)$ and $Complete(y)$ take time $O(l_y)$ in the worst case, where $l_y$ is the number of links in the subgraph induced by the chain component of $G$ that contains $y$; $Undo(y)$ takes time $O(d_y)$ in the worst case, where $d_y$ is the number of arcs in the subgraph induced by the set of descendants of $y$ in $G$ that only have one parent; $DC(x\rightarrow y)$ and $DC(x\rightarrow y\leftarrow z)$ both take time $O(dl_y)$ in the worst case, where $dl_y$ is the number of edges (either arcs or links) in the subgraph induced by the nodes in the chain component of $G$ that contains $y$ together with their descendants.

\begin{table}[htb]
\begin{tabular}{|l|l|l|} \hline
Operator & Conditions & Post-processing\\ \hline\hline

 & $x\not\in Ad_G(y)$; $p_G(x)\neq 0$ {\small or} $p_G(y)\neq 0$; & if $p_G(x)\neq 0 \mbox{ {\small and} } n_G(y)\neq 0$ \\ 
%
$A\_arc(x,y)$ & if $\left( p_G(x)\neq 0 \mbox{ {\small and} } (c_G(y)\neq 0 \mbox{ {\small or} } n_G(y)\neq 0)\right)$ & \ \ \ then $Complete(y)$ \\
%
 & \ \ \ then $DC(x\rightarrow y)=False$ & \\ \hline
%
 & $x\not\in Ad_G(y)$; $p_G(x)=0$ {\small and} $p_G(y)=0$; & \\
%
$A\_link(x,y)$ & if $\left(n_G(x)\neq 0 \mbox{ {\small and} } n_G(y)\neq 0\right)$ & None \\ 
 & \ \ \ then $UC(x\mbox{---}y)=False$ & \\ \hline
%
$D\_arc(x,y)$ & $x\in Pa_G(y)$ & if $p_G(y)\leq 2$ \\
%
 & & \ \ \ then $Undo(y)$ \\ \hline
%
$D\_link(x,y)$  & $x\in Ne_G(y)$ & None \\ \hline
%
 & $x\not\in Ad_G(y)$; $z\in Ne_G(y)$; & \\
%
$A\_hh(x,y,z)$ & $p_G(y)=0$ {\small and} $n_G(y)\neq 0$; & if $n_G(y)\geq 2$ \\
%
 & if $\left(n_G(y)\geq 2 \mbox{ {\small and} } (p_G(x)\neq 0 \mbox{ {\small or} } n_G(x)\neq 0\right)$) & \ \ \ then $Complete(y)$ \\
%
 & \ \ \ then $DC(x\rightarrow y\leftarrow z)=False$ & \\ \hline
\end{tabular}
\caption{The operators, their conditions of applicability, and post-processing requirements}
\label{condiciones}
\end{table}



\section{The Exploring Process and the Evaluation of Candidate Structures}
\label{decompo}

The search method we have described may be applied in combination with any score equivalent function $g$ (for example the AIC, BIC, MDL and BDe scoring functions are score equivalent). An easy (but inefficient) way to integrate our search method with a score equivalent function would be as follows: given an RPDAG $G$ to be evaluated, select any extension $H$ of $G$ and compute $g(H:D)$. We could also use other (non-equivalent) scoring functions, although the score of $G$ would depend on the selected extension. 

However, let us consider the case of a decomposable scoring function $g$: the DAG obtained by adding or removing an arc from the current DAG $H$ can be evaluated by modifying only one local score: 
\begin{eqnarray} \nonumber
g(H\cup\{x\rightarrow y\}:D)=g(H:D)-g_D(y, Pa_H(y))+g_D(y, Pa_H(y)\cup\{x\}) \\ \nonumber
g(H\setminus\{x\rightarrow y\}:D)=g(H:D)-g_D(y, Pa_H(y))+g_D(y, Pa_H(y)\setminus\{x\}) 
\end{eqnarray}

\citepunct{}{\&}{and}{, }{; }{, }{}{}{.}

Using decomposable scoring functions, the process of selecting, given an RPDAG, a representative DAG and then evaluating it may be quite inefficient, since we would have to recompute the local scores for all the nodes instead of only one local score. This fact can make a learning algorithm that searches in the space of equivalence classes of DAGs considerably slower than an algorithm that searches in the space of DAGs (this is the case of the algorithm proposed by \cite{chick:le}).

\citepunct{(}{\&}{and}{, }{; }{, }{)}{}{.}

Our search method can be used for decomposable scoring functions so that: (1) it is not necessary to transform the RPDAG into a DAG, the RPDAG can be evaluated directly, and (2) the score of any neighboring RPDAG can be obtained by computing at most two local scores. All the advantages of the search methods on the space of DAGs are therefore retained, but a more reduced and robust search space is used. 

Before these assertions are proved, let us examine an example. Consider the RPDAG $G$ in Figure~\ref{mievaluating} and the three neighboring configurations produced by the inclusion of an edge between $x$ and $y$, $G_1$, $G_2$ and $G_3$ (also displayed in Figure~\ref{mievaluating}). 

\begin{figure}[hbtp]
\centerline{\psfig{figure=./figuras/mievaluating.eps,height=.3\textwidth}}
\caption{An RPDAG $G$ and three neighboring configurations $G_1$, $G_2$ and $G_3$}
\label{mievaluating}
\end{figure}

The score of each of these RPDAGs is equal to the score of any of their extensions. Figure~\ref{extensions} displays one extension for each neighboring configuration. 

\begin{figure}[htb]
\centerline{\psfig{figure=./figuras/miextensions.eps,height=.3\textwidth}}
\caption{Extensions $H_1$, $H_2$ and $H_3$ of the RPDAGs $G_1$, $G_2$ and $G_3$ in Fig.~\ref{mievaluating}}
\label{extensions}
\end{figure}

We can therefore write:
%
{\small
\begin{eqnarray} \nonumber
g(G_1:D)=g(H_1:D) &=& g_D(x,\emptyset)+g_D(a,\{xb\})+g_D(b,c)+g_D(c,y)+g_D(y,x)+g_D(d,y) \\ \nonumber
g(G_2:D)= g(H_2:D) &=& g_D(x,\emptyset)+g_D(a,\{xb\})+g_D(b,c)+g_D(c,y)+g_D(y,\{xd\})+g_D(d,\emptyset) \\ \nonumber
g(G_3:D)= g(H_3:D) &=& g_D(x,\emptyset)+g_D(a,\{xb\})+g_D(b,c)+g_D(c,\emptyset)+g_D(y,\{xc\})+g_D(d,y)
\end{eqnarray} }

For each extension $H_i$ of any neighboring configuration $G_i$, it is always possible to find an extension $H_{Gi}$ of the current RPDAG $G$ such that the scores of $H_i$ and $H_{Gi}$ only differ in one local score (Figure~\ref{extensions2} displays these extensions). We can then write:

{\small 
\begin{eqnarray} \nonumber
g(G:D)=g(H_{G1}:D)&=&g_D(x,\emptyset)+g_D(a,\{xb\})+g_D(b,c)+g_D(c,y)+g_D(y,\emptyset)+g_D(d,y) \\ \nonumber
g(G:D)=g(H_{G2}:D)&=&g_D(x,\emptyset)+g_D(a,\{xb\})+g_D(b,c)+g_D(c,y)+g_D(y,d)+g_D(d,\emptyset) \\ \nonumber
g(G:D)=g(H_{G3}:D)&=&g_D(x,\emptyset)+g_D(a,\{xb\})+g_D(b,c)+g_D(c,\emptyset)+g_D(y,c)+g_D(d,y)
\end{eqnarray} }
%

\begin{figure}[htb]
\centerline{\psfig{figure=./figuras/miextensions2.eps,height=.3\textwidth}}
\caption{Three different extensions $H_{G1}$, $H_{G2}$ and $H_{G3}$ of the RPDAG $G$ in Fig.~\ref{mievaluating}}
\label{extensions2}
\end{figure}

Taking into account the previous expressions, we obtain:
%
\begin{eqnarray} \nonumber
g(G_1:D)&=&g(G:D)-g_D(y,\emptyset)+g_D(y,x) \\ \nonumber
g(G_2:D)&=&g(G:D)-g_D(y,d)+g_D(y,\{xd\}) \\ \nonumber
g(G_3:D)&=&g(G:D)-g_D(y,c)+g_D(y,\{xc\})
\end{eqnarray}

Therefore, the score of any neighboring configuration may be obtained from the score of $G$ by computing only two local scores. Note that some of these local scores may have already been computed at previous iterations of the search process: for example, $g_D(y,\emptyset)$ had to be used to score the initial empty RPDAG, and either $g_D(y,d)$ or $g_D(y,c)$ could have been computed when the link $y\mbox{---}d$ or $y\mbox{---}c$ was inserted into the structure.

\begin{prop} \label{prop4}
Let $G$ be an RPDAG and $G'$ be any RPDAG obtained by applying one of the operators described in Table~\ref{condiciones} to $G$. Let $g$ be a score equivalent and decomposable function. 
\begin{itemize}
\item[(a)] If the operator is $A\_link(x,y)$ then
\[
g(G':D)=g(G:D)-g_D(y,\emptyset)+g_D(y,\{x\})
\]
\item[(b)] If the operator is $A\_arc(x,y)$ then 
\[
g(G':D)=g(G:D)-g_D(y,Pa_G(y))+g_D(y,Pa_G(y)\cup\{x\})
\]
\item[(c)] If the operator is $A\_hh(x,y,z)$ then 
\[
g(G':D)=g(G:D)-g_D(y,\{z\})+g_D(y,\{x,z\})
\]
\item[(d)] If the operator is $D\_link(x,y)$ then
\[
g(G':D)=g(G:D)-g_D(y,\{x\})+g_D(y,\emptyset)
\]
\item[(e)] If the operator is $D\_arc(x,y)$ then 
\[
g(G':D)=g(G:D)-g_D(y,Pa_G(y))+g_D(y,Pa_G(y)\setminus\{x\})
\]
\end{itemize}
\end{prop}
%
\noindent {\bf Proof:} \\

\noindent (1) First, we shall prove that we can construct an extension $H'$ of $G'$ and another extension $H$ of $G$, such that $H$ and $H'$ differ in only one arc (this arc being $x\rightarrow y$).

\vspace{0.25cm}
\noindent $\bullet$ Consider the cases (a), (b), and (c), which correspond to the addition of an edge between $x$ and $y$: in case (a), $G'=G\cup\{x\mbox{---}y\}$ and let $H'$ be an extension of $G'$ that contains the arc $x\rightarrow y$; in case (b), where $G'=G\cup\{x\rightarrow y\}$, and in case (c), where $G'=(G\setminus\{y\mbox{---}z\})\cup\{x\rightarrow y\leftarrow z\}$, let $H'$ be any extension of $G'$ (which will contain the arc $x\rightarrow y$). In all three cases, let $H=H'\setminus\{x\rightarrow y\}$. We shall prove that $H$ is an extension of $G$:

-- First, it is obvious that $G$ and $H$ have the same skeleton. 

-- Secondly, if $u\rightarrow v\in G$ (in either case $u\rightarrow v\neq x\rightarrow y$), then $u\rightarrow v\in G'$. As $H'$ is an extension of $G'$, then $u\rightarrow v\in H'$, and this implies that $u\rightarrow v\in H$. Therefore, all the arcs in $G$ are also arcs in $H$. This result also ensures that every h-h pattern in $G$ is also an h-h pattern in $H$.

-- Thirdly, if $u\rightarrow v\leftarrow w$ is an h-h pattern in $H$ (in either case $u\rightarrow v\leftarrow w\neq x\rightarrow y\leftarrow w$), then $u\rightarrow v\leftarrow w\in H'$. Once again, as $H'$ is an extension of $G'$, we can see that $u\rightarrow v\leftarrow w\in G'$, and then $u\rightarrow v\leftarrow w\in G$. So, $G$ and $H$ have the same h-h patterns.

$H$ is therefore an extension of $G$, according to Definition~\ref{defiextension}. Note that $\forall u\neq y \;\, Pa_H(u)=Pa_{H'}(u)$ and $Pa_H(y)=Pa_{H'}(y)\setminus\{x\}$.

\vspace{0.25cm}
\noindent $\bullet$ Let us now consider cases (d) and (e), which correspond to the deletion of an edge between $x$ and $y$ (either a link or an arc, respectively): in case (d), let $H$ be an extension of $G$ containing the arc $x\rightarrow y$; in case (e), let $H$ be any extension of $G$. In both cases, let $H'=H\setminus\{x\rightarrow y\}$. We will prove that $H'$ is an extension of $G'$:

-- First, it is clear that $G'$ and $H'$ have the same skeleton. 

-- Secondly, if $u\rightarrow v\in G'$ (note that $u\rightarrow v\neq x\rightarrow y$), then $u\rightarrow v\in G$. As $H$ is an extension of $G$, then $u\rightarrow v\in H$, and therefore $u\rightarrow v\in H'$. So, all the arcs in $G'$ are also arcs in $H'$. Moreover, every h-h pattern in $G'$ is also an h-h pattern in $H'$.

-- Thirdly, if $u\rightarrow v\leftarrow w$ is an h-h pattern in $H'$ (and we know that $u\rightarrow v\leftarrow w\neq x\rightarrow y\leftarrow w$), then $u\rightarrow v\leftarrow w\in H$. As $H$ is an extension of $G$, then $u\rightarrow v\leftarrow w\in G$. Therefore, $u\rightarrow v\leftarrow w\in G'$ (the removal of the arc $x\rightarrow y$ cannot destroy any h-h pattern where $x\rightarrow y$ is not involved). So, $G'$ and $H'$ have the same h-h patterns.

In this way, $H'$ is an extension of $G'$. Moreover, we can see that $\forall u\neq y \;\, Pa_{H'}(u)=Pa_H(u)$ and $Pa_{H'}(y)=Pa_H(y)\setminus\{x\}$.

\vspace{0.25cm}
\noindent (2) The scores of $G$ and $G'$ are the same as the scores of $H$ and $H'$ respectively, since $g$ is score equivalent. Moreover, as $g$ is decomposable, we can write 
%
\begin{equation}
\begin{array}{r}
g(G':D)= g(H':D) = \sum_u g_D(u,Pa_{H'}(u)) = \sum_{u\neq y} g_D(u,Pa_{H'}(u))+ g_D(y,Pa_{H'}(y)) = \\
\sum_{u\neq y} g_D(u,Pa_{H}(u))+ g_D(y,Pa_{H}(y))- g_D(y,Pa_{H}(y))+ g_D(y,Pa_{H'}(y)) = \\
\sum_u g_D(u,Pa_{H}(u))- g_D(y,Pa_{H}(y))+ g_D(y,Pa_{H'}(y))= \\ 
 g(H:D)- g_D(y,Pa_{H}(y))+ g_D(y,Pa_{H'}(y))= \\
g(G:D)- g_D(y,Pa_H(y))+ g_D(y,Pa_{H'}(y))
\end{array}
\label{qq}
\end{equation}

Let us now consider the five different cases:

\noindent (a) In this case, we know from Table~\ref{condiciones} that $Pa_G(y)=\emptyset$. Moreover, $Pa_{G'}(y)=\emptyset$ (because we are inserting a link) and $Pa_{H'}(y)\neq\emptyset$ (because $H'$ is an extension of $G'$ that contains the arc $x\rightarrow y$). Then, from Proposition~\ref{prop3} we obtain $|Pa_{H'}(y)|=1$, i.e. $Pa_{H'}(y)=\{x\}$. Moreover, $Pa_H(y)=Pa_{H'}(y)\setminus\{x\}=\emptyset$. So, Eq.~(\ref{qq}) becomes
\[
g(G':D)= g(G:D)-g_D(y,\emptyset)+g_D(y,\{x\})
\]

\noindent (b) From Table~\ref{condiciones} we get $Pa_G(y)\neq\emptyset$ or $Pa_G(x)\neq\emptyset$. 

If $Pa_G(y)\neq\emptyset$, from Proposition~\ref{prop3} we obtain $Pa_H(y)=Pa_G(y)$. Moreover, $Pa_{H'}(y)=Pa_H(y)\cup\{x\}=Pa_G(y)\cup\{x\}$.

If $Pa_G(y)=\emptyset$ then $Pa_{G'}(y)=\{x\}$ (because we are adding the arc $x\rightarrow y$). From Proposition~\ref{prop3} we obtain $Pa_{H'}(y)=Pa_{G'}(y)=\{x\}=Pa_G(y)\cup\{x\}$. Moreover, $Pa_H(y)=Pa_{H'}(y)\setminus\{x\}=\emptyset=Pa_G(y)$.

In either case, Eq.~(\ref{qq}) becomes
\[
g(G':D)=g(G:D)-g_D(y,Pa_G(y))+g_D(y,Pa_G(y)\cup\{x\})
\]

\noindent (c) In this case, $Pa_G(y)=\emptyset$ and $Pa_{G'}(y)=\{x,z\}$. From Proposition~\ref{prop3} we obtain $Pa_{H'}(y)=\{x,z\}$. Moreover, $Pa_H(y)=Pa_{H'}(y)\setminus\{x\}=\{z\}$. Then, Eq.~(\ref{qq}) becomes
\[
g(G':D)=g(G:D)-g_D(y,\{z\})+g_D(y,\{x,z\})
\]

\noindent (d) As $Pa_G(y)=\emptyset$ and $H$ is an extension of $G$ containing the arc $x\rightarrow y$, from Proposition~\ref{prop3} we get $Pa_H(y)=\{x\}$. Moreover, $Pa_{H'}(y)=Pa_H(y)\setminus\{x\}=\emptyset$. In this case Eq.~(\ref{qq}) becomes
\[
g(G':D)=g(G:D)-g_D(y,\{x\})+g_D(y,\emptyset)
\]

\noindent (e) In this case, as $Pa_G(y)\neq\emptyset$, Proposition~\ref{prop3} asserts that $Pa_H(y)=Pa_G(y)$. Moreover, $Pa_{H'}(y)=Pa_H(y)\setminus\{x\}=Pa_G(y)\setminus\{x\}$. Therefore, Eq.~(\ref{qq}) becomes
\[
g(G':D)=g(G:D)-g_D(y,Pa_G(y))+g_D(y,Pa_G(y)\setminus\{x\})  \ \ \ \qed
\]

\subsection{Comparison with Other Approaches}

As we have already mentioned, there are several works devoted to learning Bayesian networks, within the score+search approach, which use the space of completed PDAGs to carry out the search process. There is a slight difference between the operators considered in the different works:
the addition and deletion of edges is considered by \citeA{madigan96}, within a Markov Chain Monte Carlo process, which also performs Monte Carlo sampling from the space of the orderings of the variables compatible with the current CPDAG. Edge addition and deletion is also used by \citeA{spirtesmeek95}, but within a greedy process that first grows the structure by adding edges and then thins it by deleting edges. Additional operators are considered by \citeA{chick:le}, including arc reversal and creation of v-structures. 

All these methods move through the space of completed PDAGs in the following way: given the current CPDAG $G$, after selecting an operator, applying it to $G$ and obtaining a neighboring PDAG $G'$, they generate a {\em consistent extension} $H'$ of $G'$ (a DAG belonging to the equivalence class represented by the PDAG), if one exists. If this is the case (otherwise $G'$ is not a valid configuration), then $G'$ is evaluated by computing the score of $H'$, $g(H':D)$. The completed PDAG representation of $G'$ is then recovered from its consistent extension $H'$.

The process of checking the existence of a consistent extension and generating it is carried out with a procedure called {\it PDAG-to-DAG} \cite{Dor92}, which runs in time $O(n\cdot e)$ in the worst case, where $e$ denotes the number of edges in the PDAG. Another procedure, called {\it DAG-to-PDAG}, is invoked in order to obtain the completed PDAG representation of the new valid configuration. There are different implementations of DAG-to-PDAG \cite{andersson97,chick:eq,Meek95,jverma90}. For example, the time complexity of the algorithm proposed by \citeA{chick:eq} is $O(e)$ on the average and $O(n\cdot e)$ in the worst case.

Our search method does not need to use any of these two procedures: in order to check the validity of a neighboring configuration of an RPDAG $G$, it is only necessary, in some cases, to perform a test to detect either an undirected path or a partially directed path between two nodes in $G$ (implemented by procedures $UC()$ and $DC()$ in Section~\ref{operators}). On the other hand, once the search process has explored the neighborhood of $G$ and determined the best neighboring configuration $G'$, $G'$ is not always an RPDAG, and we must generate its RPDAG representation. This generation procedure is also very simple: it consists in firing, starting from a single node $y$, a cascaded process that either directs links away from $y$ or undirects arcs (implemented by procedures $Complete()$ and $Undo()$ in Section~\ref{operators}). Note that all these procedures used by our search method are less time-consuming than PDAG-to-DAG and DAG-to-PDAG.

More importantly, our search method can take advantage of the decomposability of many scoring functions, and each RPDAG (except the initial one) can be evaluated by computing only two local scores. However, the methods based on completed PDAGs need to recompute all the local scores, although the algorithm proposed by \citeA{muntenau00}, which operates on completed PDAGs and uses three insertion operators (for arcs, links and v-structures) is also able to score any neighboring configuration using two local scores; however, the validity conditions of some of these operators are not correct. 

Finally, \citeA{Chickering02}\footnote{This work appeared after the original submission of this paper.} describes an algorithm that searches in the space of completed PDAGs and is also able to evaluate configurations by computing only (up to four) local scores. It uses six operators, link and arc addition, link and arc deletion, creation of v-structures by directing two already existing links, and reversal of arcs. All the operators can be evaluated using two local scores, except reversal and creation of v-structures, that require four local scores. The validity conditions of the operators are established essentially in terms of two conditions: (1) the absence of semi-directed or undirected paths between two nodes that do not pass through certain set of nodes, $S$, and (2) the fact that a certain set of nodes forms a clique. Link insertion and creation of v-structures need the first type of condition, link and arc deletion need the second one, whereas arc insertion and arc reversal require both conditions. The ``path'' validity conditions take time $O(|S|+e)$ in the worst case, and the ``clique'' conditions take time $O(|S|^2)$, also in the worst case. This algorithm also requires the {\it PDAG-to-DAG} and {\it DAG-to-PDAG} procedures to be used.

So, although the validity conditions of the operators in Chickering's algorithm and their postprocessing are somewhat more complex than ours, the advantage is that this algorithm does not have any duplicate representations of the equivalence classes. Whether the computational cost of moves in the CPDAG space can compensate for the larger number of RPDAGs (and the larger number of local scores to be computed) is a matter of empirical evaluation, that will possibly depend on the ``sparseness'' of the specific domain problem considered.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Experimental Results} \label{resulta}

In this section we shall describe the experiments carried out with our algorithm, the obtained results, and a comparative study with other algorithms for learning Bayesian networks. We have selected nine different problems to test our algorithm, all of which only contain discrete variables: Alarm (Figure \ref{figalarm}), Insurance (Figure~\ref{figinsurance}), Hailfinder (Figure~\ref{fighailfinder}), Breast-Cancer, crx, Flare2, House-Votes, Mushroom, and Nursery.

The Alarm network displays the relevant variables and relationships for the Alarm Monitoring System \cite{alarm}, a diagnostic application for patient monitoring. This network, which contains 37 variables and 46 arcs, has been considered as a benchmark for evaluating Bayesian network learning algorithms. The input data commonly used are subsets of the Alarm database built by \citeA{hersktesis}, which contains 20000 cases that were stochastically generated using the Alarm network. In our experiments, we have used three databases of different sizes (the first $k$ cases in the Alarm database, for $k=3000, 5000$ and $10000$).

\begin{figure}[ht]
\centerline{\psfig{figure=./figuras/alarm.eps,height=.40\textwidth}}
\caption{The Alarm network}
\label{figalarm}
\end{figure}

Insurance \cite{insurance97} is a network for evaluating car insurance risks. The Insurance network contains 27 variables and 52 arcs. In our experiments, we have used five databases containing 10000 cases, generated from the Insurance Bayesian network.

\begin{figure}[ht]
\centerline{\psfig{figure=./figuras/insurance.epsi,height=.40\textheight,width=.80\textwidth}} 
\caption{The Insurance network}
\label{figinsurance}
\end{figure}

Hailfinder~\cite{hailfinder} is a normative system that forecasts severe summer hail in northeastern Colorado. The Hailfinder network contains 56 variables and 66 arcs. In this case, we have also used five databases with 10000 cases generated from the Hailfinder network.

\begin{figure}[ht]
\centerline{\psfig{figure=./figuras/hailfinder.eps,height=.87\textheight,width=.99\textwidth}} %width}}
\caption{The Hailfinder network}
\label{fighailfinder}
\end{figure}

Breast-Cancer, crx, Flare2, House-Votes, Mushroom, and Nursery are databases available from the UCI Machine Learning Repository. Breast-Cancer contains 10 variables (9 attributes, two of which have missing values, and a binary class variable) and 286 instances. The crx database concerns credit card applications. It has 490 cases and 16 variables (15 attributes and a class variable), and seven variables have missing values. Moreover, six of the variables in the crx database are continuous and were discretized using the MLC++ system \cite{mlc}. Flare2 uses 13 variables (10 attributes and 3 class variables, one for the number of times a certain type of solar flare occured in a 24-hour period) and contains 1066 instances, without missing values. House-Votes stores the votes for each of the U.S. House of Representatives Congressmen on 16 key votes; it has 17 variables and 435 records and all the variables except two have missing values. Mushroom contains 8124 cases corresponding to species of gilled mushrooms in the Agaricus and Lepiota Family; there are 23 variables (a class variable, stating whether the mushroom is edible or poisonous, and 22 attribute variables) and only one variable has missing values. Nursery contains data relative to the evaluation of applications for nursery schools, and has 9 variables and 12960 cases, without missing values. In all of the cases, missing values are not discarded but treated as a distinct state.

In the first series of experiments, we aim to compare the behavior of our RPDAG-based local search method ({\sc rpdag}) with the classical local search in the space of DAGs ({\sc dag}). The scoring function selected is BDeu~\cite{heckgeigerchick} (which is score equivalent and decomposable), with the parameter representing the equivalent sample size set to 1 and a uniform structure prior. The starting point of the search is the empty graph in both cases.


We have collected the following information about the experiments:
\begin{description}
\item [BDeu.-] The BDeu score (log version) of the learned network.  
%
\item [Edg.-] The number of edges included in the learned network.
%
\item [H.-] The Hamming distance, H=A+D+I, i.e. the number of different edges, added (A), deleted (D), or wrongly oriented (without taking into account the differences between equivalent structures) (I), in the learned network with respect to the gold standard network (the original model). This measure is only computed for the three test domains where a gold standard exists.
%
\item [Iter.-] The number of iterations carried out by the algorithm to reach the best network, i.e. the number of operators used to transform the initial graph into a local optimum.
%
\item [Ind.-] The number of individuals (either DAGs or RPDAGs) evaluated by the algorithm.
% 
\item [EstEv.-] The number of {\it different} statistics evaluated during the execution of the algorithm. This is a useful value to measure the efficiency of the algorithms, because most of the running time of a scoring-based learning algorithm is spent in the evaluation of statistics from the database.
%
\item [TEst.-] The total number of statistics used by the algorithm. Note that this number can be considerably greater than EstEv. By using hashing techniques we can store and efficiently retrieve any previously calculated statistics. It is not therefore necessary to recompute them by accessing the database, thus gaining in efficiency.
%
\item [NVars.-] The average number of variables that intervene in the different statistics (i.e. the values $N_{y,Pa_H(y)}$ in Eq.~(\ref{decomp})) computed. This value is also important because the time required to compute a statistic increases exponentially with the number of variables involved.
%
\item [Time.-] The time, measured in {\it seconds}, employed by the algorithm to learn the network. Our implementation is written in the JAVA programming language and runs under Linux. This value is only a rough measure of the efficiency of the algorithms, because there are many circumstances that may influence the running time (external loading in a networked computer, caching or any other aspect of the computer architecture, memory paging, use of virtual memory, threading, different code, etc.). Nevertheless, we have tried to ensure that the two algorithms run under the same conditions as far as possible, and the two implementations share most of the code. In fact, the two algorithms have been integrated into the Elvira package (available at {\tt http://www.leo.ugr.es/\verb=~=elvira}).

\end{description}

For the Insurance and Hailfinder domains, the reported results are the average values across the five databases considered. The results of our experiments for synthetic data, i.e. Alarm, Insurance and Hailfinder, are displayed in Tables \ref{tab:alarm}, \ref{tab:insurance} and \ref{tab:hailfinder}, respectively, where we also show the BDeu values for the true ($T_{\mbox{{\scriptsize D}}}$) and the empty ($\emptyset_{\mbox{{\scriptsize D}}}$) networks (with parameters re-trained from the corresponding database $D$), which may serve as a kind of scale. The results obtained for real data are displayed in Table~\ref{tab:uci}.


{\small 
\begin{table}[htb] % ALARM 
\begin{center}
\begin{tabular}{|l|r|r|r|r|r|r|r|r|r|r|r|r|} \hline 
 & \multicolumn{1}{c|}{BDeu} & \multicolumn{1}{c|}{Edg} & \multicolumn{1}{c|}{H} & \multicolumn{1}{c|}{A} & \multicolumn{1}{c|}{D} & \multicolumn{1}{c|}{I} & \multicolumn{1}{c|}{Iter} & \multicolumn{1}{c|}{Ind} & \multicolumn{1}{c|}{EstEv} & \multicolumn{1}{c|}{TEst} & \multicolumn{1}{c|}{NVar} & \multicolumn{1}{c|}{Time} \\
\hline\hline
\multicolumn{13}{|c|}{Alarm-3000} \\ \hline
{\scriptsize RPDAG} & -33101 & 46 & 2 & 1 & 1 & 0 & 49 & 63124 & 3304 & 123679 & 2.98 & 111 \\
{\scriptsize DAG}   & -33109 & 47 & 7 & 3 & 2 & 2 & 58 & 72600 & 3300 & 145441 & 2.88 &  117 \\ \cline{4-13}
$T_{\mbox{{\scriptsize Alarm3}}}$ & -33114 & 46 & \multicolumn{10}{c|}{}\\ 
$\emptyset_{\mbox{{\scriptsize Alarm3}}}$ & -59890 & 0 & \multicolumn{10}{c|}{} \\ \hline
\multicolumn{13}{|c|}{Alarm-5000} \\ \hline
{\scriptsize RPDAG} & -54761 & 46 &  2 & 1 & 1 & 0 & 49 & 62869 & 3326 & 123187 & 2.97 & 179 \\
{\scriptsize DAG}   & -54956 & 54 & 16 & 9 & 1 & 6 & 60 & 76212 & 3391 & 152663 & 2.93 & 194 \\ \cline{4-13}
$T_{\mbox{{\scriptsize Alarm5}}}$ & -54774 & 46 & \multicolumn{10}{c|}{}\\ 
$\emptyset_{\mbox{{\scriptsize Alarm5}}}$ & -99983 & 0 & \multicolumn{10}{c|}{} \\ \hline
\multicolumn{13}{|c|}{Alarm-10000} \\ \hline
{\scriptsize RPDAG} & -108432 & 45 &  1 & 0 & 1 & 0 & 48 & 61190 & 3264 & 120049 & 2.97 & 346 \\
{\scriptsize DAG}   & -108868 & 52 & 13 & 7 & 1 & 5 & 60 & 75504 & 3449 & 151251 & 2.94 & 380 \\ \cline{4-13}
$T_{\mbox{{\scriptsize Alarm10}}}$ & -108452 & 46 & \multicolumn{10}{c|}{}\\ 
$\emptyset_{\mbox{{\scriptsize Alarm10}}}$ & -199920 & 0 & \multicolumn{10}{c|}{} \\ \hline\hline
\end{tabular}
\end{center}
\caption{Results for the Alarm databases}
\label{tab:alarm}
\end{table} }

{\small 
\begin{table}[htb] % INSURANCE
\begin{center}
\begin{tabular}{|l|r|r|r|r|r|r|r|r|r|r|r|r|} \hline 
 & \multicolumn{1}{c|}{BDeu} & \multicolumn{1}{c|}{Edg} & \multicolumn{1}{c|}{H} & \multicolumn{1}{c|}{A} & \multicolumn{1}{c|}{D} & \multicolumn{1}{c|}{I} & \multicolumn{1}{c|}{Iter} & \multicolumn{1}{c|}{Ind} & \multicolumn{1}{c|}{EstEv} & \multicolumn{1}{c|}{TEst} & \multicolumn{1}{c|}{NVar} & \multicolumn{1}{c|}{Time} \\
\hline\hline
{\scriptsize RPDAG} & -133071 & 45 & 18 & 4 & 10 & 4 & 48 & 36790 & 1990 & 69965 & 2.95 & 202 \\ 
{\scriptsize DAG}   & -133205 & 49 & 25 & 7 & 10 & 8 & 58 & 36178 & 2042 & 72566 & 3.03 & 214 \\ \cline{4-13}
$T_{\mbox{{\scriptsize Insur}}}$ & -133040 & 52 & \multicolumn{10}{c|}{}\\ 
$\emptyset_{\mbox{{\scriptsize Insur}}}$ & -215278 & 0 & \multicolumn{10}{c|}{} \\ \hline\hline
\end{tabular}
\end{center}
\caption{Average results for the Insurance domain across 5 databases of size 10000}
\label{tab:insurance}
\end{table} }


{\small 
\begin{table}[htb] % HAILFINDER
\begin{center}
\begin{tabular}{|l|r|r|r|r|r|r|r|r|r|r|r|r|} \hline 
 & \multicolumn{1}{c|}{BDeu} & \multicolumn{1}{c|}{Edg} & \multicolumn{1}{c|}{H} & \multicolumn{1}{c|}{A} & \multicolumn{1}{c|}{D} & \multicolumn{1}{c|}{I} & \multicolumn{1}{c|}{Iter} & \multicolumn{1}{c|}{Ind} & \multicolumn{1}{c|}{EstEv} & \multicolumn{1}{c|}{TEst} & \multicolumn{1}{c|}{NVar} & \multicolumn{1}{c|}{Time} \\
\hline\hline
{\scriptsize RPDAG} & -497872 & 67 & 24 & 12 & 10 &  2 & 68 & 235374 & 7490 & 459436 & 2.92 & 847 \\
{\scriptsize DAG}   & -498395 & 75 & 45 & 21 & 13 & 11 & 81 & 240839 & 7313 & 482016 & 2.81 & 828 \\ \cline{4-13}
$T_{\mbox{{\scriptsize Hail}}}$ & -503230 & 66 & \multicolumn{10}{c|}{}\\ 
$\emptyset_{\mbox{{\scriptsize Hail}}}$ & -697826 & 0 & \multicolumn{10}{c|}{} \\ \hline\hline
\end{tabular}
\end{center}
\caption{Average results for the Hailfinder domain across 5 databases of size 10000}
\label{tab:hailfinder}
\end{table} }

{\small 
\begin{table}[htb] % UCI
\begin{center}
\begin{tabular}{|l|r|r|r|r|r|r|r|r|} \hline 
 & \multicolumn{1}{c|}{BDeu} & \multicolumn{1}{c|}{Edg} & \multicolumn{1}{c|}{Iter} & \multicolumn{1}{c|}{Ind} & \multicolumn{1}{c|}{EstEv} & \multicolumn{1}{c|}{TEst} & \multicolumn{1}{c|}{NVar} & \multicolumn{1}{c|}{Time} \\
\hline\hline
\multicolumn{9}{|c|}{Breast-Cancer} \\ \hline
RPDAG & -2848 & 6 & 7 & 619 & 151 & 1232 & 2.26 & 0.673 \\
DAG   & -2848 & 6 & 7 & 619 & 148 & 1284 & 2.26 & 0.686  \\ \hline
\multicolumn{9}{|c|}{crx} \\ \hline
RPDAG & -5361 & 19 & 20 & 5318 & 545 & 10020 & 2.75 & 3.03 \\
DAG   & -5372 & 19 & 20 & 4559 & 510 &  9208 & 2.61 & 2.85 \\ \hline
\multicolumn{9}{|c|}{Flare2} \\ \hline
RPDAG & -6728 & 15 & 16 & 2637 & 329 & 4887 & 2.71 & 3.66 \\
DAG   & -6733 & 13 & 14 & 2012 & 310 & 4093 & 2.45 & 3.23 \\ \hline
\multicolumn{9}{|c|}{House-Votes} \\ \hline
RPDAG & -4629 & 22 & 23 & 6370 & 591 & 11883 & 2.80 & 3.10  \\
DAG   & -4643 & 23 & 24 & 6094 & 621 & 12289 & 2.66 & 3.39 \\ \hline
\multicolumn{9}{|c|}{Mushroom} \\ \hline
RPDAG & -77239 & 92 &  97 & 43121 & 2131 & 78459 & 3.99 & 432 \\
DAG    & -77208 & 87 & 103 & 39944 & 2173 & 80175 & 3.96 & 449 \\ \hline
\multicolumn{9}{|c|}{Nursery} \\ \hline
RPDAG & -125717 & 8 & 9 & 415 & 115 &  803 & 2.75 & 11.91 \\
DAG   & -125717 & 8 & 9 & 611 & 133 & 1269 & 2.59 & 13.50 \\  
\hline\hline
\end{tabular}
\end{center}
\caption{Results for the UCI databases}
\label{tab:uci}
\end{table} }

As we consider there to be a clear difference between the results obtained for the synthetic and the UCI domains, we shall discuss them separately.
%
\begin{itemize}
\item For the synthetic domains (Tables \ref{tab:alarm}, \ref{tab:insurance} and \ref{tab:hailfinder}):
\begin{itemize}
\item Our RPDAG-based algorithm outperforms the DAG-based one with respect to the value of the scoring function used to guide the search: we always obtain better results on the five databases considered. Note that we are using a logarithmic version of the scoring function, so that the differences are much greater in a non-logarithmic scale. These results support the idea that RPDAGs are able to find new and better local maxima within the score+search approach for learning Bayesian networks in this type of highly structured domains.
%
\item Our search method is also preferable from the point of view of the Hamming distances, which are always considerably lower than the ones obtained by using the DAG space.
%
\item Moreover, our search method is generally more efficient: it carries out fewer iterations (on the five cases), evaluates fewer individuals (on four cases), computes fewer different statistics from the databases (on three cases), uses fewer statistics (on the five cases), and runs faster (on four cases). On the contrary, the average number of variables involved in the statistics is slightly greater (on four cases).
%
\end{itemize}
%
\item For the UCI domains (Table~\ref{tab:uci}):
\begin{itemize}
\item The results in this case are not as conclusive about the advantages of the RPDAG-based method with respect to the DAG-based one in terms of effectiveness: both algorithms reach the same solution on two cases from six, {\sc rpdag} is better than {\sc dag} on three cases, and {\sc dag} is better on one case.
%
\item With respect to the efficiency of the two algorithms, the situation is similar: neither algorithm clearly outperforms the other with respect to any of the five efficiency measures considered.
\end{itemize}
%
\end{itemize}


In a second series of experiments, we aim to test the behavior of the search in the RPDAG space when used in combination with a search heuristic which is more powerful than a simple greedy search. The heuristic selected is {\it Tabu Search} \cite{Glover,Bouckaert95}, which tries to escape from a local maximum by selecting a solution that minimally decreases the value of the scoring function; immediate re-selection of the local maximum just visited is prevented by maintaining a list of solutions that are forbidden, the so-called {\it tabu list} (although for practical reasons the tabu list stores forbidden operators not solutions, and consequently, solutions which have not been visited previously may also become forbidden).

We have implemented two simple versions of tabu search: {\sc ts-rpdag} and {\sc ts-dag}, which explore the RPDAG and DAG spaces, respectively, using the same operators as their respective greedy versions. The parameters used by these algorithms are the length $tll$ of the tabu list and the number $tsit$ of iterations required to stop the search process. In our experiments, these values have been fixed as follows: $tll=n$ and $tsit=n(n-1)$, $n$ being the number of variables in the domain. The scoring function and the initial graph are the same as in previous experiments, as well as the collected performance measures, with one exception: as the number of iterations is now fixed (Iter=$tsit$), we compute the iteration where the best graph was found (BIter) instead. The results of these experiments are displayed in Tables \ref{tab:tsalarm}, \ref{tab:tsinsurance}, \ref{tab:tshailfinder} and \ref{tab:tsuci}.

{\small 
\begin{table}[htb] % ALARM-TS
\begin{center}
\begin{tabular}{|l|r|r|r|r|r|r|r|r|r|r|r|r|} \hline 
 & \multicolumn{1}{c|}{BDeu} & \multicolumn{1}{c|}{Edg} & \multicolumn{1}{c|}{H} & \multicolumn{1}{c|}{A} & \multicolumn{1}{c|}{D} & \multicolumn{1}{c|}{I} & \multicolumn{1}{c|}{BIter} & \multicolumn{1}{c|}{Ind} & \multicolumn{1}{c|}{EstEv} & \multicolumn{1}{c|}{TEst} & \multicolumn{1}{c|}{NVar} & \multicolumn{1}{c|}{Time} \\
\hline\hline
\multicolumn{13}{|c|}{Alarm-3000} \\ \hline
{\scriptsize TS-RPDAG} & -33101 & 46 &  2 & 1 & 1 & 0 &  48 & 1779747 & 9596 & 3044392 & 3.63 & 286 \\
{\scriptsize TS-DAG}   & -33115 & 51 & 11 & 7 & 2 & 2 & 129 & 1320696 & 8510 & 2645091 & 3.59 & 391 \\  \hline
\multicolumn{13}{|c|}{Alarm-5000} \\ \hline
{\scriptsize TS-RPDAG} & -54761 & 46 & 2 & 1 & 1 & 0 &  48 & 1779579 & 10471 & 3031966 & 3.61 & 421 \\
{\scriptsize TS-DAG}   & -54762 & 47 & 3 & 2 & 1 & 0 & 720 & 1384990 & 11113 & 2773643 & 3.58 & 541 \\ \hline
\multicolumn{13}{|c|}{Alarm-10000} \\ \hline
{\scriptsize TS-RPDAG} & -108432 & 45 & 1 & 0 & 1 & 0 &  47 & 1764670 & 10671 & 3020165 & 3.66 & 735 \\
{\scriptsize TS-DAG}   & -108442 & 50 & 6 & 5 & 1 & 0 & 284 & 1385065 & 11014 & 2773795 & 3.60 & 862 \\ \hline\hline
\end{tabular}
\end{center}
\caption{Results for the Alarm databases using Tabu Search}
\label{tab:tsalarm}
\end{table} }

{\small 
\begin{table}[htb] % INSURANCE-TS
\begin{center}
\begin{tabular}{|l|r|r|r|r|r|r|r|r|r|r|r|r|} \hline 
 & \multicolumn{1}{c|}{BDeu} & \multicolumn{1}{c|}{Edg} & \multicolumn{1}{c|}{H} & \multicolumn{1}{c|}{A} & \multicolumn{1}{c|}{D} & \multicolumn{1}{c|}{I} & \multicolumn{1}{c|}{BIter} & \multicolumn{1}{c|}{Ind} & \multicolumn{1}{c|}{EstEv} & \multicolumn{1}{c|}{TEst} & \multicolumn{1}{c|}{NVar} & \multicolumn{1}{c|}{Time} \\
\hline\hline
{\scriptsize TS-RPDAG} & -133070 & 45 & 18 & 4 & 10 & 4 &  58 & 458973 & 2823 & 751551 & 3.44 & 182 \\ 
{\scriptsize TS-DAG}   & -132788 & 47 & 18 & 5 & 10 & 3 & 415 & 352125 & 5345 & 706225 & 4.16 & 428 \\  \hline\hline
\end{tabular}
\end{center}
\caption{Average results for the Insurance domain using Tabu Search }
\label{tab:tsinsurance}
\end{table} }

{\small 
\begin{table}[htb] % HAILFINDER-TS
\begin{center}
\begin{tabular}{|l|r|r|r|r|r|r|r|r|r|r|r|r|} \hline 
 & \multicolumn{1}{c|}{BDeu} & \multicolumn{1}{c|}{Edg} & \multicolumn{1}{c|}{H} & \multicolumn{1}{c|}{A} & \multicolumn{1}{c|}{D} & \multicolumn{1}{c|}{I} & \multicolumn{1}{c|}{BIter} & \multicolumn{1}{c|}{Ind} & \multicolumn{1}{c|}{EstEv} & \multicolumn{1}{c|}{TEst} & \multicolumn{1}{c|}{NVar} & \multicolumn{1}{c|}{Time} \\
\hline\hline
{\scriptsize TS-RPDAG} & -497872 & 67 & 24 & 12 & 10 & 2 &   67 & 9189526 & 19918 & 1.5223387E7 & 4.07 & 3650 \\ 
{\scriptsize TS-DAG}   & -498073 & 70 & 35 & 17 & 13 & 5 & 1631 & 7512114 & 22184 & 1.5031642E7 & 4.07 & 4513  \\  \hline\hline
\end{tabular}
\end{center}
\caption{Average results for the Hailfinder domain using Tabu Search }
\label{tab:tshailfinder}
\end{table} }


{\small 
\begin{table}[htb] % UCI-TS
\begin{center}
\begin{tabular}{|l|r|r|r|r|r|r|r|r|} \hline 
 & \multicolumn{1}{c|}{BDeu} & \multicolumn{1}{c|}{Edg} & \multicolumn{1}{c|}{BIter} & \multicolumn{1}{c|}{Ind} & \multicolumn{1}{c|}{EstEv} & \multicolumn{1}{c|}{TEst} & \multicolumn{1}{c|}{NVar} & \multicolumn{1}{c|}{Time} \\
\hline\hline
\multicolumn{9}{|c|}{Breast-Cancer} \\ \hline
TS-RPDAG & -2848 & 6 & 6 & 8698 & 345 & 14209 & 3.03 & 1.96 \\
TS-DAG   & -2848 & 6 & 6 & 6806 & 316 & 13892 & 2.87 & 1.98 \\ \hline
\multicolumn{9}{|c|}{crx} \\ \hline
TS-RPDAG & -5361 & 19 & 19 & 61175 &  908 & 98574 & 3.20 &  9.26 \\
TS-DAG   & -5362 & 20 & 29  & 44507 & 1176 & 89714 & 3.17 & 12.23 \\ \hline
\multicolumn{9}{|c|}{Flare2} \\ \hline
TS-RPDAG & -6728 & 15 &  15 & 23363 & 616 & 37190 & 3.47 & 9.70 \\
TS-DAG   & -6726 & 15 & 129 & 18098 & 681 & 36665 & 3.27 & 10.60\\ \hline
\multicolumn{9}{|c|}{House-Votes} \\ \hline
TS-RPDAG & -4622 & 24 & 180 & 73561 & 1144 & 121206 & 3.32 & 11.14  \\
TS-DAG   & -4619 & 23 & 252 & 56570 & 1364 & 113905 & 3.29 & 15.30 \\ \hline
\multicolumn{9}{|c|}{Mushroom} \\ \hline
TS-RPDAG  & -77002 & 99 & 495 & 209556 & 4021 & 350625 & 4.83 &  883 \\
TS-DAG    & -77073 & 90 & 450 & 157280 & 3455 & 315975 & 4.57 & 1725 \\ \hline
\multicolumn{9}{|c|}{Nursery} \\ \hline
TS-RPDAG & -125717 & 8 & 8 & 4352 & 251 & 6525 & 3.09 & 28.05\\
TS-DAG   & -125717 & 8 & 8 & 3898 & 237 & 7991 & 2.95 & 26.20\\  
\hline\hline
\end{tabular}
\end{center}
\caption{Results for the UCI databases using Tabu Search}
\label{tab:tsuci}
\end{table} }


For the synthetic domains, in all the cases, except in one of the insurance databases, the results obtained by {\sc ts-rpdag} and {\sc rpdag} are the same. This phenomenon also appears for the UCI databases, where only in two databases does {\sc ts-rpdag} improve the results of {\sc rpdag}. Therefore, the Tabu Search does not contribute significantly to improving the greedy search in the RPDAG space (at least using the selected values for the parameters $tll$ and $tsit$). This is in contrast with the situation in the DAG space, where {\sc ts-dag} improves the results obtained by {\sc dag}, with the exception of two UCI databases (equal results) and Alarm-3000 (where {\sc dag} performs better than {\sc ts-dag}).

With respect to the comparison between {\sc ts-rpdag} and {\sc ts-dag}, we still consider {\sc ts-rpdag} to be preferable to {\sc ts-dag} on the synthetic domains, although in this case {\sc ts-dag} performs better on the insurance domain. For the UCI databases, the two algorithms perform similarly: each algorithm is better than the other on two domains, and both algorithms perform equally on the remaining two domains. {\sc ts-rpdag} is somewhat more efficient than {\sc ts-dag} with respect to running time.


We have carried out a third series of experiments to compare our learning algorithm based on RPDAGs with other algorithms for learning Bayesian networks. In this case, the comparison is only intended to measure the quality of the learned network. In addition to the DAG-based local and tabu search previously considered, we have also used the following algorithms:
%
\begin{itemize}
\item {\sc pc}~\cite{Spirtes93}, an algorithm based on independence tests. We used an independence test based on the measure of conditional mutual information \cite{Kullback68}, with a fixed confidence level equal to 0.99. 
%
\item The K2 search method~\cite{CooHersk:k2}, in combination with the BDeu scoring function ({\sc k2}). Note that {\sc k2} needs an ordering of the variables as the input. We used an ordering consistent with the topology of the corresponding networks.
%
\item Another algorithm, BN Power Constructor ({\sc bnpc}), that uses independence tests~\cite{Cheng97,bnpc98}. 
%
\end{itemize}
%
The two independence-based algorithms, {\sc pc} and {\sc bnpc}, operate on the space of equivalence classes, whereas {\sc k2} explores the space of DAGs which are compatible with a given ordering. We have included the algorithm {\sc k2} in the comparison, using a correct ordering and the same scoring function as the RPDAG and DAG-based search methods, in order to test whether our method can outperform the results obtained with a more informed algorithm. The results for the algorithms {\sc pc} and {\sc k2} have been obtained using our own implementations (which are also included in the Elvira software). For {\sc bnpc}, we used the software package available at {\tt http://www.cs.ualberta.ca/\verb=~=jcheng/bnsoft.htm}. 

The test domains included in these experiments are Alarm, Insurance, and Hailfinder. In addition to the BDeu values, the number of edges in the learned networks and the Hamming distances, we have collected two additional performance measures:
%
\begin{description}
\item [BIC.-] The value of the BIC (Bayesian Information Criterion) scoring function \cite{Schwarz78} for the learned network. This value measures the quality of the network using maximum likelihood and a penalty term. Note that BIC is also score-equivalent and decomposable.
%
\item [KL.-] The Kullback-Leibler distance (cross-entropy) \cite{Kullback68} between the probability distribution, $P$, associated to the database (the empirical frequency distribution) and the probability distribution associated to the learned network, $P_G$. Notice that this measure is actually the same as the log probability of the data. We have in fact calculated a decreasing monotonic linear transformation of the Kullback-Leibler distance, because this one has exponential complexity and the transformation can be computed very efficiently:  
If $P_{G}$ is the joint probability distribution associated to a network $G$, then the KL distance can be written in the following way~\cite{campos98,Lam94}:
\begin{equation}
KL(P,P_G)= -H_P({\cal U})+\sum_{x\in {\cal U}} H_P(x)-\!\!\!\!\!\!\!\!\sum_{x\in {\cal U},Pa_G(x)\neq\emptyset} 
\!\!\!\!\!\!\!\! MI_P(x,Pa_G(x))
\label{kull}
\end{equation}
where $H_P(Z)$ denotes Shannon entropy with respect to the distribution $P$ for the subset of variables $Z$ and $MI_P(x,Pa_G(x))$ is the measure of mutual information between the two sets of variables $\{x\}$ and $Pa_G(x)$. As the first two terms of the expression above do not depend on the graph $G$, our transformation consists in calculating only the third term in equation~(\ref{kull}). So, the interpretation of our transformation of the Kullback-Leibler distance is: the higher this value is, the better the network fits the data. However, this measure should be handled with caution, since a high KL value may also indicate overfitting (a network with many edges will probably have a high KL value).
%
\end{description}

Although for those algorithms whose goal is to optimize the Bayesian score, BDeu is really the metric that should be used to evaluate them, we have also computed BIC and KL because two of the algorithms considered use independence tests instead of a scoring function.

The results of these experiments are displayed in Table~\ref{tab:otros}. The best value for each performance measure and each database is written in bold, and the second best value in italics.
These results indicate that our search method in the RPDAG space, in combination with the BDeu scoring function, is competitive with respect to other algorithms: only the {\sc ts-dag} algorithm, which uses a more powerful (and more computationally intensive) search heuristic in the DAG space, and, to a lesser extent, the more informed {\sc k2} algorithm, perform better than {\sc rpdag} in some cases. Observe that both {\sc ts-dag} and {\sc k2} perform better than {\sc rpdag} in terms of KL on four cases from five.

%\vspace{-1.50cm}
{\small 
\begin{table}[htb] % comparativa 
\begin{center}
\begin{tabular}{|l|r|r|r|r|r|r|r|r|} \hline 
 & \multicolumn{1}{c|}{BDeu} & \multicolumn{1}{c|}{BIC} &  \multicolumn{1}{c|}{KL} & \multicolumn{1}{c|}{Edg} & \multicolumn{1}{c|}{H} & \multicolumn{1}{c|}{A} & \multicolumn{1}{c|}{D} & \multicolumn{1}{c|}{I}  \\
\hline\hline
\multicolumn{9}{|c|}{Alarm-3000} \\ \hline
RPDAG & {\bf-33101} & {\bf -33930} & {\it 9.23055} & 46 & {\bf 2} & 1 & 1 & 0 \\
DAG   & {\it -33109} & {\it -33939} & 9.23026 & 47 & {\it 7} & 3 & 2 & 2 \\
TS-DAG& -33115 & -33963 & 9.23047 & 51 & 11 & 7 & 2 & 2 \\
PC    & -36346 & -36691 & 8.06475 & 37 & 10 & 0 & 9 & 1 \\ 
BNPC  & -33422 & -35197 & 9.11910 & 43 &  {\it 7} & 2 & 5 & 0 \\
K2    & -33127 & -34351 & {\bf 9.23184} & 46 & {\bf 2} & 1 & 1 & 0\\ \hline
\multicolumn{9}{|c|}{Alarm-5000} \\ \hline
RPDAG & {\bf -54761} & {\bf -55537} & 9.25703 & 46 & {\bf 2} & 1 & 1 & 0 \\
DAG   & -54956 & -55831 & 9.25632 & 54 & 16 & 9 & 1 & 6 \\
TS-DAG& {\it -54762} & {\it -55540} & {\it 9.25736} & 47 & {\it 3} & 2 & 1 & 0 \\
PC    & -61496 & -61822 & 7.85435 & 38 & 16 & 2 & 10 & 4 \\
BNPC  & -55111 & -55804 & 9.16787 & 42 & 4 & 0 & 4 & 0 \\
K2    & -54807 & -55985 & {\bf 9.25940} & 47 & {\it 3} & 2 & 1 & 0 \\ \hline
\multicolumn{9}{|c|}{Alarm-10000} \\ \hline
RPDAG & {\bf -108432} & {\bf -109165} & 9.27392 & 45 &  {\bf 1} & 0 & 1 & 0 \\
DAG   & -108868 & -110537 & {\bf 9.27809} & 52 & 13 & 7 & 1 & 5 \\
TS-DAG& {\it -108442} & {\it -109188} & 9.27439 & 50 & 6 & 5 & 1 & 0 \\
PC    & -117661 & -117914 & 8.31704 & 38 & 11 & 1 & 9 & 1 \\
BNPC  & -109164 & -109827 & 9.18884 & 42 &  4 & 0 & 4 & 0 \\
K2    & -108513 & -109647 & {\it 9.27549} & 46 &  {\it 2} & 1 & 1 & 0 \\ \hline
\multicolumn{9}{|c|}{Insurance} \\ \hline
RPDAG & -133071 & -134495 & 8.38502 & 45 & {\it 18} & 4 & 10 & 4 \\ 
DAG   & -133205 & -135037 & 8.39790 & 49 & 25 & 7 & 10 & 8 \\
TS-DAG& {\it -132788} & {\it -134414} & {\it 8.41467} & 47 & {\it 18} & 5 & 10 & 3    \\  
PC    & -139101 & -141214 & 7.75574 & 33 & 23 & 0 & 20 & 3 \\
BNPC  & -134726 & -135832 & 8.21606 & 37 & 26 & 3 & 18 & 5 \\
K2    & {\bf -132615} & {\bf -134095} & {\bf 8.42471} & 44 & {\bf 10} & 1 & 9 & 0 \\ \hline
\multicolumn{9}{|c|}{Hailfinder} \\ \hline
RPDAG & {\bf -497872} & -531138 & 20.53164 & 67 & {\it 24} & 12 & 10 &  2 \\
DAG   & -498395 & -531608 & 20.48503 & 75 & 45 & 21 & 13 & 11 \\ 
TS-DAG& {\it -498073} & {\it -516100} & {\it 20.59372} & 70 & 35 & 17 & 13 & 5 \\
PC    & -591507 & -588638 & 12.65981 & 49 & 50 & 16 & 33 & 1 \\
BNPC  & -503440 & {\bf -505160} & {\bf 20.61581} & 64 & 28 & 12 & 15 & 1 \\
K2    & -498149 & -531373 & 20.51822 & 67 & {\bf 23} & 12 & 11 & 0 \\ 
\hline\hline
\end{tabular}
\end{center}
\caption{Performance measures for different learning algorithms}
\label{tab:otros}
\end{table} }

The fourth series of experiments attempts to evaluate the behavior of the same algorithms on a dataset which is different from the training set used to learn the network. In order to do so, we have computed the BDeu, BIC, and KL values of the network structure learned using a database, with respect to a different database: for the Alarm domain, the training set is the Alarm-3000 database used previously, and the test set is formed by the 3000 next cases in the Alarm database; for both Insurance and Hailfinder, we selected one of the five databases that we have been using as the training set and another of these databases as the test set. The results are shown in Table~\ref{tab:testset}. 
We can observe that they are analogous to the results obtained in Table~\ref{tab:otros}, where the same databases were used for training and testing. However, in this case {\sc rpdag} also performs better than {\sc ts-dag} and {\sc k2} in terms of KL. Therefore, the good behavior of our algorithm cannot be attributed to overfitting.

{\small 
\begin{table}[htb] % comparativa 
\begin{center}
\begin{tabular}{|l|r|r|r|} \hline 
 & \multicolumn{1}{c|}{BDeu} & \multicolumn{1}{c|}{BIC} &  \multicolumn{1}{c|}{KL} \\
\hline\hline
\multicolumn{4}{|c|}{Alarm-3000} \\ \hline
RPDAG & {\bf -32920} & {\bf -33750} & {\it 9.35839} \\
DAG   & {\it -32938} & {\it -33769} & 9.35488 \\
TS-DAG& -32947 & -33793 & 9.35465 \\
PC    & -36286 & -36632 & 8.15227\\ 
BNPC  & -33309 & -35051 & 9.23576 \\
K2    & -32951 & -34165 & {\bf 9.36168} \\ \hline
\multicolumn{4}{|c|}{Insurance} \\ \hline
RPDAG & -132975 & -134471 & {\bf 8.44891} \\ 
DAG   & -133004 & -134952 & {\it 8.44745} \\
TS-DAG& {\bf -132810} & {\bf -134281} & 8.44538 \\
PC    & -140186 & -143011 & 7.70182 \\
BNPC  & -135029 & -136125 & 8.23063 \\
K2    & {\it -132826} & {\it -134309} & 8.44671 \\ \hline
\multicolumn{4}{|c|}{Hailfinder} \\ \hline
RPDAG & {\bf -497869} & -531165 & {\it 20.58267} \\
DAG   & -498585 & -531913 & 20.49730 \\ 
TS-DAG& {\it -497983} & {\bf -505592} & 20.55420 \\
PC    & -587302 & -584939 & 12.68256 \\
BNPC  & -506680 & {\it -508462} & {\bf 20.71069} \\
K2    & -498118 & -531356 & 20.57876 \\ 
\hline\hline
\end{tabular}
\end{center}
\caption{Performance measures for the learning algorithms using a different test set}
\label{tab:testset}
\end{table} }

Finally, we have carried out another series of experiments, which aim to compare our {\sc rpdag} algorithm with the algorithm proposed by \citeA{Chickering02}, that searches in the CPDAG space. In this case, we have selected the House-Votes and Mushroom domains (which were two of the datasets used by \citeauthor{Chickering02}). In order to approximate our experimental conditions to those described in Chickering's work, we used the BDeu scoring function with a prior equivalent sample size of ten, and a structure prior of $0.001^f$, where $f$ is the number of free parameters in the DAG; moreover, we used five random subsets of the original databases, each containing approximately 70\% of the total data (304 cases for House-Votes and 5686 for Mushroom). Table \ref{compara} displays the average values across the five datasets of the relative improvement of the per-case score obtained by {\sc rpdag} to the per-case score of {\sc dag}, as well as the ratio of the time spent by {\sc dag} to the time spent by {\sc rpdag}. We also show in Table \ref{compara} the corresponding values obtained by Chickering (using only one dataset) for the comparison between his {\sc cpdag} algorithm and {\sc dag}.

{\small 
\begin{table}[htb] % comparacion con Chickering 
\begin{center}
\begin{tabular}{l|c|c||c|c|} \cline{2-5}
 & \multicolumn{2}{|c||}{RPDAG versus DAG} & \multicolumn{2}{|c|}{CPDAG versus DAG} \\ \hline %\cline{2-5}
\multicolumn{1}{|c|}{} & Relative & Time  & Relative & Time \\
\multicolumn{1}{|c|}{Dataset} & Improv.  & Ratio & Improv.  & Ratio \\ \hline
\multicolumn{1}{|c|}{House-Votes} & 0.0000 & 1.041 &  0.0000 & 1.27 \\
\multicolumn{1}{|c|}{Mushroom}    & 0.0158 & 1.005 & -0.0382 & 2.81 \\
\hline\hline
\end{tabular}
\end{center}
\caption{Comparison with Chickering's work on completed PDAGs}
\label{compara}
\end{table} }

We may observe that the behavior of {\sc rpdag} and {\sc cpdag} is somewhat different: although both algorithms are more efficient than {\sc dag}, it seems to us that {\sc cpdag} runs faster than {\sc rpdag}. With respect to effectiveness, both {\sc rpdag} and {\sc cpdag} obtain exactly the same solution as {\sc dag} in the House-Votes domain (no relative improvement); however, in the other domain, {\sc rpdag} outperforms {\sc dag} (on the five datasets considered) whereas {\sc cpdag} performs worse than {\sc dag}. In any case, the differences are small (they could be a result of differences in the experimental setup) and a much more systematic experimentation with these algorithms would be necessary in order to establish general conclusions about their comparative behavior.

\section{Concluding Remarks} \label{futuro}

We have developed a new local search algorithm, within the score+search approach for learning Bayesian network structures from databases. The main feature of our algorithm is that it does not search in the space of DAGs, but uses a new form of representation, restricted PDAGs, that allows us to search efficiently in a space similar to the space of equivalence classes of DAGs. For the common situation in which a decomposable scoring function is used, the set of operators that define the neighborhood structure of our search space can be scored locally (as it happens in the space of DAGs), i.e. we can evaluate any neighboring restricted PDAG by computing at most two local scores. In this way, we maintain the computational efficiency that the space of DAGs offers and, at the same time, we explore a more reduced search space, with a smoother landscape, which avoids some early decisions on edge directions. These characteristics may help to direct the search process towards better network structures.

The experimental results show that our search method based on restricted PDAGs can efficiently and accurately recover complex Bayesian network structures from data, and can compete with several state of the art Bayesian network learning algorithms, although it does not significantly improve them. Our experiments in Section \ref{resulta}, as well as those conducted by \citeA{Chickering02}, seem to point out that search algorithms based on PDAGs can obtain slightly better results, with respect to both effectiveness and efficiency, than search methods based on DAGs, especially for highly structured models (i.e., models that can be (almost) perfectly represented by a DAG). We believe that PDAGs can also be useful in domains which are complex (contain many variables and complicated dependence patterns) and sparse (represent many independence relationships).

For future research, we are planning to integrate the techniques developed in this paper within more powerful search methods, such as the ones considered by \citeA{Blanco02,ijar02} or \citeA{campospuertaISAS01}. Additionally, in the light of the results obtained by our method in combination with Tabu Search, it may be interesting to incorporate another operator, which could either be a classical arc reversal or some kind of specific operator to destroy h-h patterns. We also intend to work on the adaptation and application of our algorithm to real problems in the field of classification. 


%\vspace{0.4cm}
%\noindent {\large\bf Acknowledgments} 

\acks{This work has been supported by the Spanish Ministerio de Ciencia y Tecnolog\'{\i}a (MCYT) and the Junta de Comunidades de Castilla-La Mancha under Projects TIC2001-2973-CO5-01 and PBC-02-002, respectively. We are grateful to our colleague Jos\'e M. Puerta for his invaluable help with the implementation of several algorithms. We are also grateful to the three anonymous reviewers for useful comments and suggestions.}

 
\bibliographystyle{theapa}
%\bibliography{referencesfin.tex}
\input{referencesfin.tex}
\end{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
