s2fieldsofstudy:
Computer Science
Business

text:
Digital Ecosystems: Evolving Service-Orientated Architectures

We view Digital Ecosystems to be the digital counterparts of biological ecosystems, exploiting the self-organising properties of biological ecosystems, which are considered to be robust, self-organising and scalable architectures that can automatically solve complex, dynamic problems. Digital Ecosystems are a novel optimisation technique where the optimisation works at two levels: a first optimisation, migration of agents (representing services) which are distributed in a decentralised peer-to-peer network, operating continuously in time; this process feeds a second optimisation based on evolutionary computing that operates locally on single peers and is aimed at finding solutions to satisfy locally relevant constraints. We created an Ecosystem-Oriented Architecture of Digital Ecosystems by extending Service-Oriented Architectures with distributed evolutionary computing, allowing services to recombine and evolve over time, constantly seeking to improve their effectiveness for the user base. Individuals within our Digital Ecosystem will be applications (groups of services), created in response to user requests by using evolutionary optimisation to aggregate the services. These individuals will migrate through the Digital Ecosystem and adapt to find niches where they are useful in fulfilling other user requests for applications. Simulation results imply that the Digital Ecosystem performs better at large scales than a comparable Service-Oriented Architecture, suggesting that incorporating ideas from theoretical ecology can contribute to useful self-organising properties in digital ecosystems.

I. INTRODUCTION
Is mimicking ecosystems the future of information systems? A key challenge in modern computing is to develop systems that address complex, dynamic problems in a scalable and efficient way, because the increasing complexity of software makes designing and maintaining efficient and flexible systems a growing challenge [2], [3], [4]. What with the ever expanding number of services being offered online from Application Programming Interfaces (APIs) being made public, there is an ever growing number of computational units available to be combined in the creation of applications. However, this is currently a task done manually by programmers, and it has been argued that current software development techniques have hit a complexity wall [5], which can only be overcome by automating the search for new algorithms. There are several existing efforts aimed at achieving this automated service composition [6], [7], [8], [9], the most prevalent of which is Service-Oriented Architectures and its associated standards and technologies [10], [11].
Alternatively, nature has been in the research business for 3.8 billion years and in that time has accumulated close to This work is supported by the European Commission under the EU project Digital Business Ecosystems (contract number 507953 [1]). 30 million well-adjusted solutions to a plethora of design challenges that humankind struggles to address with mixed results [12]. Biomimicry is a discipline that seeks solutions by emulating nature's designs and processes, and there is considerable opportunity to learn elegant solutions for human-made problems [12]. Biological ecosystems are thought to be robust, scalable architectures that can automatically solve complex, dynamic problems, possessing several properties that may be useful in automated systems. These properties include self-organisation, self-management, scalability, the ability to provide complex solutions, and automated composition of these complex solutions [13].
Therefore, an approach to the aforementioned challenge would be to develop Digital Ecosystems, artificial systems that aim to harness the dynamics that underlie the complex and diverse adaptations of living organisms in biological ecosystems. While evolution may be well understood in computer science under the auspices of evolutionary computing [14], ecological models are not. The possible connections between Digital Ecosystems and their biological counterparts are yet to be closely examined, so potential exists to create an Ecosystem-Oriented Architecture with the essential elements of biological ecosystems, where the word ecosystem is more than just a metaphor. We propose that an ecosystem inspired approach, would be more effective at greater scales than traditionally inspired approaches, because it would be built upon the scalable and self-organising properties of biological ecosystems [13].

II. SERVICE-ORIENTED ARCHITECTURES
Our approach to evolving high-level software applications requires a modular reusable paradigm to software development. Service-oriented architectures (SOAs) are the current state-of-the-art approach, being the current iteration of interface/component-based design from the 1990s, which was itself an iteration of event-oriented design from the 1980s, and before then modular programming from the 1970s [15], [16]. Service-oriented computing promotes assembling application components into a loosely coupled network of services, to create flexible, dynamic business processes and agile applications that span organisations and computing platforms [17]. This is achieved through a SOA, an architectural style that guides all aspects of creating and using business processes throughout their life-cycle, packaged as services. This includes defining and provisioning the infrastructure that allows different applications to exchange data and participate in business processes, arXiv:0712.4102v6 [cs.NE] 5 Oct 2009 loosely coupled from the operating systems and programming languages underlying the applications [18]. Hence, a SOA represents a model in which functionality is decomposed into distinct units (services), which can be distributed over a network, and can be combined and reused to create business applications [17].
A SOA depends upon service-orientation as its fundamental design principle. In a SOA environment, independent services can be accessed without knowledge of their underlying platform implementation [18]. Services reflect a service-oriented approach to programming that is based on composing applications by discovering and invoking networkavailable services to accomplish some task. This approach is independent of specific programming languages or operating systems, because the services communicate with each other by passing data from one service to another, or by coordinating an activity between two or more services [17]. So, the concepts of SOAs are often seen as built upon, and the development of, the concepts of modular programming and distributed computing [16].
SOAs allow for an information system architecture that enables the creation of applications that are built by combining loosely coupled and interoperable services [18]. They typically implement functionality most people would recognise as a service, such as filling out an online application for an account, or viewing an online bank statement [16]. Services are intrinsically un-associated units of functionality, without calls to each other embedded in them. Instead of services embedding calls to each other in their source code, protocols are defined which describe how services can talk to each other, in a process known as orchestration, to meet new or existing business system requirements [19]. This is allowing an increasing number of third-party software companies to offer software services, such that SOA systems will come to consist of such third-party services combined with others created in-house, which has the potential to spread costs over many users and uses, and promote standardisation both in and across industries [20]. For example, the travel industry now has a well-defined and documented, set of both services and data, sufficient to allow any competent software engineer to create travel agency software using entirely off-the-shelf software services [21], [22]. Other industries, such as the finance industry, are also making significant progress in this direction [23].
The vision of SOAs assembling application components from a loosely coupled network of services that can create dynamic business processes and agile applications that span organisations and computing platforms, is visualised in Figure  1. It will be made possible by creating compound solutions that use internal organisational software assets, including enterprise information and legacy systems, and combining these solutions with external components residing in remote networks [24]. The great promise of SOAs is that the marginal cost of creating the n-th application is virtually zero, as all the software required already exists to satisfy the requirements of other applications. Only their combination and orchestration are required to produce a new application [25], [26]. The key is that the interactions between the chunks, are not specified within the chunks themselves. Instead, the interaction of services (all of whom are hosted by un-associated peers) is specified by users in an ad-hoc way, with the intent driven by newly emergent business requirements [27].
The pinnacle of SOA interoperability, is the exposing of services on the internet as web services [18]. A web service is a specific type of service that is identified by a Uniform Resource Identifier (URI), whose service description and transport utilise open Internet standards. Interactions between web services typically occur as Simple Object Access Protocol (SOAP) calls carrying eXtensible Markup Language (XML) data content. Interface descriptions of the web services are expressed using the Web Services Definition Language (WSDL) [28]. The Universal Description Discovery and Integration (UDDI) standard defines a protocol for directory services that contain web service descriptions. UDDI enables web service clients to locate candidate services and discover their details. Service clients and service providers utilise these standards to perform the basic operations of SOAs [28]. Service aggregators can then use the Business Process Execution Language (BPEL) to create new web services by defining corresponding compositions of the interfaces and internal processes of existing services [28].
SOA services inter-operate based on a formal definition (or contract, e.g. WSDL) that is independent of the underlying platform and programming language. Service descriptions are used to advertise the service capabilities, interface, behaviour, and quality [28]. The publication of such information about available services provides the necessary means for discovery, selection, binding, and composition of services [28]. The (expected) behaviour of a service during its execution is described by its behavioural description (for example, as a workflow process). Also, included is a quality of service (QoS) description, which publishes important functional and non-functional service quality attributes, such as service metering and cost, performance metrics (response time, for instance), security attributes, integrity (transactional), reliability, scalability, and availability [28]. Service clients (end-user organisations that use some service) and service aggregators (organisations that consolidate multiple services into a new, single service offering) utilise service descriptions to achieve their objectives [28]. One of the most important and continuing developments in SOAs is the use of semantic descriptions for service discovery, so that a client can discover the services semantically, and then apply transformations to adapt the interface of the services to the interface expected, using already available client software [29].
There are multiple standards available and still being developed for SOAs [11], most notably of recent being REpresentational State Transfer (REST) [19]. The software industry now widely implements a thin SOAP/WSDL/UDDI veneer atop existing applications or components that implement the web services paradigm [24], but the choice of technologies could change with time. Therefore, SOAs and its services are best defined generically, because SOAs are technology agnostic and need not be tied to a specific technology [17]. Within the current and future scope of SOAs, there is clearly potential to evolve complex high-level software applications from the modular services of SOAs, instead of the instruction level evolution currently prevalent in genetic programming [30].

III. DISTRIBUTED EVOLUTIONARY COMPUTING
The fact that evolutionary computing manipulates a population of independent solutions actually makes it well suited for parallel computation architectures [31]. The motivation for using parallel or distributed evolutionary algorithms is twofold. First, improving the speed of evolutionary processes by conducting concurrent evaluations of individuals in a population. Second, improving the problem-solving process by overcoming difficulties that face traditional evolutionary algorithms, such as maintaining diversity to avoid premature convergence [32], [33]. There are several variants of distributed evolutionary computing, leading some to propose a taxonomy for their classification [34], with there being two main forms of models [31], [33]: • multiple-population/coarse-grained migration/island • single-population/fine-grained diffusion/neighbourhood In the coarse-grained island models [35], [31], evolution occurs in multiple parallel sub-populations (islands), each running a local evolutionary algorithm, evolving independently with occasional migrations of highly fit individuals among sub-populations. The core parameters for the evolutionary algorithm of the island-models are as follows [35]: • number of the sub-populations: 2, 3, 4, more • sub-population homogeneity size, crossover rate, mutation rate, migration interval • topology of connectivity: ring, star, fully-connect, random • static or dynamic connectivity • migration mechanisms: isloated/synchronous/asynchronous how often migrations occur which individuals migrate Fine-grained diffusion models [36], [33] assign one individual per processor. A local neighbourhood topology is assumed, and individuals are allowed to mate only within their neighbourhood, called a deme. The demes overlap by an amount that depends on their shape and size, and in this way create an implicit migration mechanism. Each processor runs an identical evolutionary algorithm which selects parents from the local neighbourhood, produces an offspring, and decides whether to replace the current individual with an offspring. However, even with the advent of multi-processor  [35], [31]: There are different probabilities of going from island 1 to island 2 , as there is of going from island 2 to island 1 . This mirrors the naturally inspired quality that although two populations have the same physical separation, it may be easier to migrate in one direction than the other, i.e. fish migration is easier downstream than upstream.
computers, and more recently multi-core processors, which provide the ability to execute multiple threads simultaneously [4], this approach would still prove impractical in supporting the number of agents necessary to create a Digital Ecosystem. Therefore, we shall further consider the island models.
An example island-model [35], [31] is visualised in Figure  2, in which there are different probabilities of going from island 1 to island 2 , as there is of going from island 2 to island 1 . This allows maximum flexibility for the migration process, and mirrors the naturally inspired quality that although two populations have the same physical separation, it may be easier to migrate in one direction than the other, i.e. fish migration is easier downstream than upstream. The migration of the island models is like the notion of migration in nature, being similar to the metapopulation models of theoretical ecology [37]. This model has also been used successfully in the determination of investment strategies in the commercial sector, in a product known as the Galapagos toolkit [38], [39]. However, all the islands in this approach work on exactly the same problem, which makes it less analogous to biological ecosystems in which different locations can be environmentally different [40]. We will take advantage of this property later when defining the Ecosystem-Oriented Architecture of Digital Ecosystems.

IV. THE DIGITAL ECOSYSTEM
We are concerned with the digital counterpart of biological ecosystems. However, the term digital ecosystem has been used to describe a variety of concepts, which it now makes sense to review. Some of these refer to the existing networking infrastructure of the internet [41], [42], [43], while several companies offer a digital ecosystem service or solution, which involves enabling customers to use existing e-business solutions [44], [45], [46]. The term is also being increasingly linked, yet undefined, to the future developments of Information and Communications Technology (ICT) adoption for e-business and e-commerce, to create so called business ecosystems [47], [48], [49]. However, perhaps the most frequent references to digital ecosystems arise in Artificial Life research, where they are created primarily to  Fig. 3.
Digital Ecosystem: Optimisation architecture in which agents (representing services) travel along the peer-to-peer connections; in every node (habitat) local optimisation is performed through an evolutionary algorithm, where the search space is determined by the agents present at the node.
investigate aspects of biological and other complex systems [50], [51], [52]. The extent to which these disparate systems resemble biological ecosystems varies, and frequently the word ecosystem is merely used for branding purposes without any inherent ecological properties.
We consider Digital Ecosystems [53], [54], [55] to be software systems that exploit the properties of biological ecosystems, which are robust, scalable, and self-organising [13]. So, Digital Ecosystems provide a two-level optimisation scheme inspired by natural ecosystems, in which a decentralised peer-to-peer network forms an underlying tier of distributed agents. These agents then feed a second optimisation level based on an evolutionary algorithm that operates locally on single habitats (peers), aiming to find solutions that satisfy locally relevant constraints. The local search is sped up through this twofold process, providing better local optima as the distributed optimisation provides prior sampling of the search space by making use of computations already performed in other peers with similar constraints [53]. The agents consist of an executable component and an ontological description [56]. So, the Digital Ecosystem can be considered a Multi-Agent System (MAS) [56] which uses distributed evolutionary computing [31], [33] to combine suitable agents in order to meet user requests for applications.
The landscape, in energy-centric biological ecosystems, defines the connectivity between habitats [40]. Connectivity of nodes in the digital world is generally not defined by geography or spatial proximity, but by information or semantic proximity. For example, connectivity in a peer-to-peer network is based primarily on bandwidth and information content, and not geography. The island-models of distributed evolutionary computing use an information-centric model for the connectivity of nodes (islands) [35]. However, because it is generally defined for one-time use (to evolve a solution to one problem and then stop) it usually has a fixed connectivity between the nodes, and therefore a fixed topology [31]. So, supporting evolution in the Digital Ecosystem, with a multiobjective selection pressure (fitness landscape [57] with many peaks), requires a re-configurable network topology, such that habitat connectivity can be dynamically adapted based on the observed migration paths of the agents between the users within the habitat network. Based on the island-models of distributed evolutionary computing [35], each connection between the habitats is bi-directional and there is a probability associated with moving in either direction across the connection, with the connection probabilities affecting the rate of migration of the agents. However, additionally, the connection probabilities will be updated by the success or failure of agent migration using the concept of Hebbian learning [58]: the habitats which do not successfully exchange agents will become less strongly connected, and the habitats which do successfully exchange agents will achieve stronger connections. This leads to a topology that adapts over time, resulting in a network that supports and resembles the connectivity of the user base. If we consider a business ecosystem, network of Small and Medium sized Enterprises, as an example user base; such business networks are typically small-world networks [59], [60]. They many strongly connected clusters (communities), called sub-networks (quasi-complete graphs), with a few connections between these clusters (communities) [61]. Graphs with this topology have a very high clustering coefficient and small characteristic path lengths [61]. So, the Digital Ecosystem will take on a topology similar to that of the user base.
The novelty of our approach comes from the evolving populations being created in response to similar requests. So whereas in the island-models of distributed evolutionary computing there are multiple evolving populations in response to one request [35], here there are multiple evolving populations in response to similar requests. In our Digital Ecosystems different requests are evaluated on separate islands (populations), and so adaptation is accelerated by the sharing of solutions between evolving populations (islands), because they are working to solve similar requests (problems).
The users will formulate queries to the Digital Ecosystem by creating a request as a semantic description, like those being used and developed in SOAs [29], specifying an application they desire and submitting it to their local peer (habitat). This description defines a metric for evaluating the fitness of a composition of agents, as a distance function between the semantic description of the request and the agents' ontological descriptions. A population is then instantiated in the user's habitat in response to the user's request, seeded from the agents available at their habitat. This allows the evolutionary optimisation to be accelerated in the following three ways: first, the habitat network provides a subset of the agents available globally, which is localised to the specific user it represents; second, making use of agent-sequences previously evolved in response to the user's earlier requests; and third, taking advantage of relevant agentsequences evolved elsewhere in response to similar requests by other users. The population then proceeds to evolve the optimal agent-sequence(s) that fulfils the user request, and as the agents are the base unit for evolution, it searches the available agent-sequence combination space. For an evolved agent-sequence that is executed (instantiated) by the user, it then migrates to other peers (habitats) becoming hosted where it is useful, to combine with other agents in other populations to assist in responding to other user requests for applications.

V. SIMULATION AND RESULTS
An important measure for determining the success of the Digital Ecosystem is its relative performance to a traditional SOA based system. So, we simulated a simple SOA with a distributed UDDI service registry [17], with redirects (links to other nodes) at each node for service descriptions not stored locally, and without caching (analogous to the Digital Ecosystem). We then compared it against a typical simulation run of the Digital Ecosystem. The time available to the SOA based system for fulfilling a request from the distributed registry was limited to the time the Digital Ecosystem required to respond to the same user request, and worked by searching for the optimal services based on the optimal segmentation of the request, because an exhaustive combinatorial search would have been impractical. We simulated the Digital Ecosystem, following the Ecosystem-Oriented Architecture from the previous section. Throughout the simulations we assumed a hundred users, which meant that at any time the number of users joining the network equalled those leaving. The habitats of the users were randomly connected at the start, to simulate the users going online for the first time. The users then produced agents (services) and requests for business applications. Initially, the users each deployed five agents to their habitats, for migration (distribution) to any habitats connected to theirs (i.e. their community). Users were simulated to deploy a new agent after the submission of three requests for business applications, and were chosen at random to submit their requests.
A simulated user request consisted of an abstract semantic description, as a list of sets of numeric tuples to represent the properties of a desired business application. The use of the numeric tuples made it comparable to the semantic descriptions of the services represented by the agents; while the list of sets (two level hierarchy) and a much longer length provided sufficient complexity to support the sophistication of business applications.
In the Digital Ecosystem user requests were handled by the habitats instantiating evolving populations, which used evolutionary computing to find the optimal solution(s), agentsequence(s). It was assumed that the users made their requests for business applications accurately, and always used the response (agent-sequence) provided. Populations of agents, [A 1 , A 1 , A 2 , ...], were evolved to solve user requests, seeded with agents and agent-sequences from the agent-pool of the habitats in which they were instantiated. A dynamic population size was used to ensure exploration of the available combinatorial search space, which increased with the average length of the population's agent-sequences. The optimal combination of agents (agent-sequence) was evolved to the user request, R, by an artificial selection pressure created by a fitness function generated from the user request, R. An individual (agent-sequence) of the population consisted of a set of attributes, a 1 , a 2 , ..., and a user request essentially consisted of a set of required attributes varied according to a Gaussian distribution, r 1 , r 2 , .... So, the fitness function for evaluating an individual agent-sequence, A, relative to a user request, R, was, where a is the member of A such that the difference to the required attribute r was minimised. Equation 1 was used to assign fitness values between 0.0 and 1.0 to each individual of the current generation of the population, directly affecting their ability to replicate into the next generation. The evolutionary computing process was encoded with a Graph of the performance of the Digital Ecosystem against a traditional SOA based system: Both the Digital Ecosystem and the SOA based system performed as expected, providing better responses to user requests as more services became available. The SOA reference system initially performed better than the Digital Ecosystem, but with the increasing number of services the Digital Ecosystem outperformed the reference system. low mutation rate, a fixed selection pressure and a nontrapping fitness function (i.e. did not get trapped at local optima). The type of selection used fitness-proportional and non-elitist. Fitness-proportional meaning that the fitter the individual the higher its probability of surviving to the next generation [62]. Non-elitist meaning that the best individual from one generation was not guaranteed to survive to the next generation; it had a high probability of surviving into the next generation, but it was not guaranteed as it might have been mutated, [14]. Crossover (recombination) was then applied to a randomly chosen 10% of the surviving population, a one-point crossover, by aligning two parent individuals and picking a random point along their length, and at that point exchanging their tails to create two offspring [14]. Mutations were then applied to a randomly chosen 10% of the surviving population; point mutations were randomly located, consisting of insertions (an agent was inserted into an agent-sequence), replacements (an agent was replaced in an agent-sequence), and deletions (an agent was deleted from an agent-sequence) [63]. The issue of bloat was controlled by augmenting the fitness function with a parsimony pressure [64] which biased the search to shorter agent-sequences, evaluating longer than average length agent-sequences with a reduced fitness, and thereby providing a dynamic control limit which adapted to the average length of the ever-changing evolving agent populations.
In Figure 4 we graphed the percentage match to the user requests for typical runs, as determined by a distance function between the request and the service descriptions. Both the Digital Ecosystem and the SOA based system performed as expected, providing better responses to user requests as more services became available. The SOA reference system initially performed better than the Digital Ecosystem, but with the increasing number of services the Digital Ecosystem outperformed the reference system. This was anticipated as the the Digital Ecosystem was expected to be more effective at larger scales [53].

VI. CONCLUSION
The experimental results indicate that under simulation conditions the Digital Ecosystem outperforms the comparison system based on a traditional Service-Oriented Architecture. Service-oriented architectures promise to provide potentially huge numbers of services that programmers can combine via standardised interfaces, to create increasingly sophisticated and distributed applications [28]. The Digital Ecosystem extends this concept with the automatic combining of available and applicable services in a scalable architecture to meet user requests for applications. This is made possible by a fundamental paradigm shift, from a pull-oriented approach to a push-oriented approach. So, instead of the pull-oriented approach of generating applications only upon request in Service-Oriented Architectures [19], the Digital Ecosystem follows a push-oriented approach of distributing and composing applications pre-emptively, as well as upon request. Although the use of Service-Oriented Architectures in the definition of Digital Ecosystems provides a predisposition to business [16], it does not preclude other more general uses. The Ecosystem-Oriented Architecture definition of Digital Ecosystems is intended to be inclusive and interoperable with other technologies, in the same way that the definition of Service-Oriented Architectures is with grid computing and other technologies [19]. For example, habitats could be executed using a distributed processing arrangement, such as cloud computing [65], which would be possible because the habitat network topology is information-centric (instead of location-centric).


Genetic Programming for Kernel-based Learning with Co-evolving Subsets Selection

Support Vector Machines (SVMs) are well-established Machine Learning (ML) algorithms. They rely on the fact that i) linear learning can be formalized as a well-posed optimization problem; ii) non-linear learning can be brought into linear learning thanks to the kernel trick and the mapping of the initial search space onto a high dimensional feature space. The kernel is designed by the ML expert and it governs the efficiency of the SVM approach. In this paper, a new approach for the automatic design of kernels by Genetic Programming, called the Evolutionary Kernel Machine (EKM), is presented. EKM combines a well-founded fitness function inspired from the margin criterion, and a co-evolution framework ensuring the computational scalability of the approach. Empirical validation on standard ML benchmark demonstrates that EKM is competitive using state-of-the-art SVMs with tuned hyper-parameters.

Introduction
Kernel methods, including the so-called Support Vector Machines (SVMs), are well-established learning approaches with both strong theoretical foundations and successful practical applications [1]. SVMs rely on two main advances in statistical learning. First, the linear supervised machine learning task is set as a well-posed (quadratic) optimization problem. Second, the above setting is extended to non-linear learning via the kernel trick : given a (manually designed) change of representation Φ mapping the initial space onto the so-called feature space, linear hypotheses are characterized in terms of the scalar product in the feature space, or kernel. These hypotheses correspond to non-linear hypotheses in the initial space. Although many specific kernels have been proposed in the literature, designing a kernel well suited for an application domain or a dataset so far remains an art more than a science. This paper proposes a system, the Evolutionary Kernel Machine (EKM), for the automatic design of data-specific kernels. EKM applies Genetic Programming (GP) [2] to construct symmetric functions (kernels), and optimizes a fitness function inspired from the margin criterion [3]. Kernels are assessed within a Nearest Neighbor classification process [4,5]. In order to cope with computational complexity, a cooperative co-evolution governs the prototype subset selection and the GP kernel design, while the fitness case subset selection undergoes a competitive co-evolution.
The paper is organized as follows. Section 2 introduces the formal background and notations on kernel methods. Sections 3 and 4 respectively describe the GP representation and the fitness function proposed for the EKM. Scalability issues are addressed in the co-evolutionary framework introduced in Section 5.
Results on benchmark problems are given in Section 6. Finally, related works are discussed in Section 7 before concluding the paper in Section 8.

Formal Background and Notations
Supervised machine learning takes as input a dataset E = {(x i , y i ), i = 1 . . . n, x i ∈ X, y i ∈ Y }, made of n examples; x i and y i respectively stand for the description and the label of the i-th example. The goal is to construct a hypothesis h(x) mapping X onto Y with minimal generalization error. Only vectorial domains (X = IR d ) are considered throughout this paper; further, only binary classification problems (Y = {1, −1}) are considered in the rest of this section.
Due to space limitations, the reader is referred to [6] for a comprehensive presentation of SVMs. In the simplest (linear separable) case, the hyper-plane h(x) maximizing the geometrical margin (distance to the closest examples) is constructed. The label associated to example x is the sign of h(x), with: where < x, x i > denotes the scalar product of x and x i . Let Φ denotes a mapping from the instance space X onto the feature space and let the kernel K(x, x ′ ) be defined as: Under some conditions (the kernel trick), non-linear classifiers on X are constructed as in the linear case, and characterized as h( Besides SVMs, the kernel trick can be used to revisit all learning methods involving a distance measure. In the paper, the kernel nearest neighbor (Kernel-NN) algorithm [5], which revisits the k-nearest neighbors (k-NN) [4], is considered. Given a distance (or dissimilarity) function d(x, x ′ ) defined on the instance space X, given a set of labelled examples E = {(x 1 , y 1 ), . . . , (x n , y n )} and an instance x to be classified, the k-NN algorithm: i) determines the k examples closest to x according to d(x, x ′ ); ii) outputs the majority class of these k examples. Kernel-NN proceeds as k-NN, where distance d K (x, x ′ ) is defined after the kernel K(x, x ′ ) (more on this in Section 4).
Standard kernels on X = IR d include Gaussian and polynomial kernels 3 . It must be noted that the addition, multiplication and compositions of kernels are kernels, and therefore the standard SVM machinery can find the optimal value of hyper-parameters (e.g. σ, c or k) among a finite set. Quite the opposite, the functional (symbolic) optimization of K(x, x ′ ) cannot be tackled to our best knowledge except by Genetic Programming.

Genetic Programming of Kernels
The Evolutionary Kernel Machine applies GP to determine symmetric functions K(x, x ′ ) on IR d × IR d best suited to the dataset at hand. As shown in Table 1, the main difference compared to standard symbolic regression is that terminals are symmetric expressions of x and x ′ (e.g. , enforcing the symmetry of the kernels (K(x, x ′ ) = K(x ′ , x)).
The initialization of GP individuals is done using a ramped half and half procedure [2]. The selection probability of terminals A i , M i , I i and S i (respectively C i,j ) is divided by 1/d (resp. 2/d(d + 1)), where d is the dimension of the initial instance space (X = IR d ).
Indeed the kernel functions built after Table 1 might not satisfy Mercer's condition (K(x, x) ≤ 0 ⇒ x = 0) required for SVM optimization [6]. However these kernels will be assessed along a Kernel-NN classification rule [5]; therefore the fact that they are not necessarily positive is not a limitation. Quite the contrary, EKM kernels can achieve feature selection; typically, terminals associated to non-informative features should disappear along evolution. The use of EKM for feature selection will be examined in a future work.

Fitness Measure
Every kernel K(x, x ′ ) is assessed after the Kernel-NN classification rule, using the dissimilarity d K defined as Given a prototype set E p = {(x 1 , y 1 ), . . . , (x ℓ , y ℓ )} and a training example e = (x, y), let us assuming that E p is ordered by increasing dissimilarity to ). Let p(e) denotes the minimum rank over all prototype examples in the same class as e (p(e) = min{i, y i = y, i = 1 . . . ℓ}); let n(e) denotes the minimum rank over all other prototype examples (not belonging to the same class as e, n(e) = min{i, y i = y, i = 1 . . . ℓ}).
As noted by [3], the quality of the Kernel-NN classification of e can be assessed from δ K (e) = n(e) − p(e). The higher δ K (e), the more confident the classification of e is, e.g. with respect to perturbations of E p or d K ; δ K (e) measures the margin of e with respect to Kernel-NN.
Crossed multiplication-addition between the i th and j th com- Accordingly, given a prototype set E p = {(x 1 , y 1 ), . . . , (x ℓ , y ℓ )} and a fit- The computation of F has linear complexity in the number ℓ of prototypes and in the number m of fitness cases. In a standard setting, E p and E s both coincide with the whole training set E (ℓ = m = n). However the quadratic complexity of the fitness computation with respect to the number n of training examples is incompatible with the scalability of the approach.

Tractability Through Co-evolution
EKM scalability is obtained along two directions, by i) reducing the number ℓ of prototypes used for classification, and ii) reducing the size m of the fitness case subset considered during each generation. ; set E t p to the best offspring; (d) Generate λs offsprings of E t−1 s , by replacing a fraction ρs of the fitness case (uniform stratified sampling); assess these offsprings after h t, * and E t p ; set E t s to the best offspring. 5. Output h * * , selected among h t, * , t = 0 . . . T as the one minimizing the 1-NN error rate on the whole training set E using the associated E t p prototype subset. More precisely, a co-evolutionary framework involving three species is considered, as detailed in Figure 1. The first species includes the GP kernels. The second species includes the prototype subset (fixed-size subsets of the training set), subject to a cooperative co-evolution [7] with the GP kernels. The third species includes the fitness case subset (fixed-size subsets of the training set), subject to a competitive host-parasite co-evolution [8] with the GP kernels.
The prototype species is evolved to find good prototypes such that they maximize the fitness of the GP kernels. The fitness case species is evolved to find hard and challenging examples, such that they minimize the kernel fitness. Of course there is a danger that the fitness case subset ultimately capture the noisy examples, as observed in the boosting framework [9] (see Section 6.2).
Both prototype and selection species are initialized using a stratified uniform sampling with no replacement (the class distribution in the sample is the same as in the whole dataset and all examples are distinct). Both species are evolved using a (1, λ) evolution strategy; in each generation, λ offsprings are generated using a uniform stratified replacement of a given fraction of the parent subset, and assessed after the best kernel in the current kernel population. The parent subset is replaced by the best offspring. In each generation, the kernels are assessed after the current prototype and fitness case individuals.

Experimental Validation
This section reports on the experimental validation of EKM, on a standard set of benchmark problems [10], detailed in Table 2. The system is implemented using the Open BEAGLE framework 4 for evolutionary computation [11].

Experimental Setting
The parameters used in EKM are reported in Table 3. The average evolution time for one run is less than one hour (AMD Athlon 2800+). On each problem, EKM has been evaluated along the standard 10-fold cross validation methodology. The whole data set is partitioned into 10 (stratified) subsets; the training set is made of all subsets but one; the best hypothesis learned from this training set is evaluated on the remaining subset, or test set. The accuracy is averaged over the 10 folds (as the test set ranges over the 10 subsets of the whole dataset); for each fold, EKM is launched 10 times; the 5 best hypotheses (after their accuracy on the training set) are assessed on the test set; the reported accuracy is the average over the 10 folds of these 5 best hypotheses on the test set. In total, EKM is launched 100 times on each problem.
EKM is compared to state of the art algorithms, including k-nearest neighbor and SVMs with Gaussian kernels, similarly assessed using 10-fold cross validation. For k-NN, the underlying distance is the Euclidean one, and scaling Table 3. Tableau of the evolutions parameters.

Parameter
Description and parameter values GP kernel functions evolution parameters Primitives See Table 1 normalization option has been considered; the k parameter has been varied in {1, 3, 5}; the best setting has been kept. For Gaussian SVMs, the Torch3 implementation has been used [12]; the error cost (parameter C) has been varied in {10 i , i = −3 . . . 4}, the σ parameter is set to 10, and the best setting has been similarly retained. Table 4 shows the results obtained by EKM compared with k-NN and Gaussian SVM, together with the optimal parameters for the latter algorithms. The size of the best GP kernel (last column) shows that no bloat occurred, thanks to the lexicographic parsimony pressure. Each algorithm is shown to be the best performing on the half or more of the tested datasets, with frequent ties according to a paired Student's t-test.

Results
Typically, the problems where Gaussian SVMs perform well are those where the optimal C value for cost error is high, suggesting that the noise level in these datasets is high too. Indeed, the fitness case subset selection embedded in EKM might favor the selection of noisy examples, as those are more challenging to GP Table 4. Comparative 10-fold results of k-NN, Gaussian SVM and EKM on the UCI data sets, with optimal settings (k and scaling for k-NN, C for SVM). The reported test error is averaged over the 10 folds. For each fold tested with the EKM, the 5 solutions out of 10 runs with best training error are assessed on the test set, and their error is averaged. Test error rates in bold denotes the statistically best results according to a 95% two-tails paired Student's t-test. "Average rank" column gives the test error ranking obtained for EKM compared to k-NN and SVM averaged over the 10 folds. kernels. A more progressive selection mechanism, taking into account all kernels in the GP population to better filter out noisy examples and outliers, will be considered in further research.
The k-NN outperforms SVM and EKM on the bos problem, where the noise level appears to be very low. Indeed, the optimal value for the number k of nearest neighbors is k = 1, while the optimal cost error is 10 −3 , suggesting that the error rate is also low. Still, the fact that the error rate is close to 23% might be explained as the target concept is complex and/or many examples lie close to its frontier. On bcw, the differences between the three algorithms are not statistically different and the test error rate is about 2%, suggesting that the problem is rather easy.
EKM is found to outperform the other algorithms on bld, demonstrating that Kernel-based dissimilarity can improve on Euclidean distance with and without rescaling. Last, EKM behaves like k-NN on the pid problems. Further, it must be noted that EKM classifies the test examples using a 50-examples prototype set, whereas k-NN uses the whole training set (above 300 examples in the bld problem and 690 in the pid problem).
As the well-known No Free Lunch theorem applies to Machine Learning too, no learning method is expected to be universally competent. Rather, the above experimental validation demonstrates that the GP-evolved kernels can improve on standard kernels in some cases.
The most relevant work to EKM is the Genetic Kernel Support Vector Machine (GK-SVM) [13]. GK-SVM similarly uses GP within an SVM-based approach, with two main differences compared to EKM. On one hand, GK-SVM focuses on feature construction, using GP to optimize mapping Φ (instead of the kernel). On the other hand, the fitness function used in GK-SVM suffers from a quadratic complexity in the number of training examples. Accordingly, all datasets but one considered in the experimentations are small (less than 200 examples). On a larger dataset, the authors acknowledge that their approach does not improve on a standard SVM with well chosen parameters. Another related work similarly uses GP for feature construction, in order to classify time series [14]. The set of features (GP trees) is further evolved using a GA, where the fitness function is based on the accuracy of an SVM classifier. Most other works related to evolutionary optimization within SVMs (see [15]) actually focus on parametric optimization, e.g. achieving features selection or tuning some parameters.
Another related work is proposed by Weinberger et al. [16], optimizing a Mahalanobis distance based on the k-NN margin criterion inspired from [3] and also used in EKM. However, restricted to linear changes of representation, the optimization problem is tackled by semi-definite programming in [16]. Lastly, EKM is also inspired by the Dynamic Subset Selection first proposed by Gathercole and Ross [17] and further developed by [18] to address scalability issues in EC-based Machine Learning.

Conclusion
The Evolutionary Kernel Machine proposed in this paper aims to improve kernelbased nearest neighbor classification [5], combining two original aspects. First, EKM implicitly addresses the feature construction problem by designing a new representation of the application domain better suited to the dataset at hand. However, in contrast with [13,14], EKM takes advantage of the kernel trick, using GP to optimize the kernel function. Secondly, EKM proposes a co-evolution framework to ensure the scalability of the approach and control the computational complexity of the fitness computation. The empirical validation demonstrates that this new approach is competitive with well-founded learning algorithms such as SVM and k-NN using tuned hyper-parameters.
A limitation of the approach, also observed in the well-known boosting algorithm [9], is that the competitive co-evolution of kernels and examples tends to favor noisy validation examples. A perspective for further research is to exploit the evolution archive, to estimate the probability for an example to be noisy and achieve a sensitivity analysis. Another perspective is to incorporate ensemble learning, typically bagging and boosting, within EKM. Indeed the diversity of the solutions constructed along population-based optimization enables ensemble learning almost for free.


extfieldsofstudy:
Computer Science
Mathematics
Physics

text:
Boolean Functions, Projection Operators and Quantum Error Correcting Codes

This paper describes a fundamental correspondence between Boolean functions and projection operators in Hilbert space. The correspondence is widely applicable, and it is used in this paper to provide a common mathematical framework for the design of both additive and non-additive quantum error correcting codes. The new framework leads to the construction of a variety of codes including an infinite class of codes that extend the original ((5,6,2)) code found by Rains [21]. It also extends to operator quantum error correcting codes.

I. INTRODUCTION
The additive or stabilizer construction of quantum error correcting codes (QECC) takes a classical binary code that is self-orthogonal with respect to a certain symplectic inner product, and produces a quantum code, with minimum distance determined by the classical code (for more details see [7], [8] and [14]). The first non-additive quantum error-correcting code was constructed by Rains et al. [21]. This code was constructed numerically by building a projection operator with a given weight distribution. Grassl and Beth [13] generalized this construction by introducing union quantum codes, where the codes are formed by taking the sum of subspaces generated by two quantum codes. Roychowdhury and Vatan [23] gave some sufficient conditions for the existence of nonadditive codes, and Arvind et al. [5] developed a theory of non-additive codes based on the Weyl commutation relations. Most recently, Kribs et al. [16] introduced operator quantum error correction (OQEC) which unifies the standard error correction model, the method of decoherence-free subspaces, and that of noiseless subsystems.
We will describe, what we believe to be the first mathematical framework for code design that encompasses both additive and non-additive quantum error correcting codes. It is based on a correspondence between Boolean functions and projection operators in Hilbert space that is described in Sections II and III. We have used an initial version of this correspondence to construct Grassmannian packings [1] and space-time codes for wireless communication [3]. However, the correspondence in Section III applies to a larger class of projection operators and includes the correspondence described in [3] as a special This work was supported in part by AFOSR under contract 00852833. The material in this paper was presented in part at the IEEE International Symposium on Information Theory, Nice, France, June 2007. case (see Section IV). We note that prior work by Danielson [11] interpreted Boolean functions as quantum states and developed a correspondence between Boolean functions and zero-dimensional quantum codes.
After introducing the fundamentals of quantum error correcting codes in Section V, we will derive in Section VI sufficient conditions for existence of QECC in terms of existence of certain Boolean function. This paper goes beyond deriving sufficient conditions, and constructs the quantum code if these properties are satisfied. Hence, we convert the problem of finding a quantum code into a problem of finding Boolean function satisfying certain properties. We also see how certain well-known codes fit into this scheme. We focus on nondegenerate codes which is defensible given that we know of no parameters k, M and d for which there exists a ((k, M, d)) degenerate QECC but not a ((k, M, d)) non-degenerate QECC (see [2]). Further, in Section VII, we describe how this scheme fits into a general framework of operator quantum error correcting codes. More precisely, we give sufficient conditions for the existence of ((k, M, N, d)) stabilizer OQEC and also construct the code if these conditions are satisfied.

II. BOOLEAN FUNCTION
A Boolean function is defined as a mapping f : where y j is the value of the Boolean function at the decimal index j and c 0 (j), c 1 (j), .... , c m−1 (j) ∈ {0, 1} are the coordinates in the binary representation for j (with c m−1 as the most significant bit and c 0 as the least significant bit) with v 1 j = v j and v 0 j =v j (Theorem 7.7, [18]).

Example 1:
The truth vector of the three-variable Boolean Definition 2 ( [20]): Let ⊕ denote modulo two addition. The (periodic) autocorrelation function of a Boolean function f (v) at a is the inner product of f with a shift of f by a. More precisely, r(a) = This means that for any element a in the The complementary set links distinguishability in the quantum world (orthogonality of subspaces) with properties of Boolean functions. The quantity f (v ⊕ a) is the counterpart in the quantum world of the quantum subspace after the error has occurred, which is to be orthogonal to the original subspace corresponding to f (v) as will be described in later sections.
If the supports of f (v), f (v⊕a) intersect in N decimal indices then Hence, N = 0 and a ∈ Cset f .

III. BOOLEAN FUNCTIONS AND A LOGIC OF PROJECTION OPERATORS
The authors of [3] connected Boolean logic to projection operators derived from the Heisenberg-Weyl group. In this section, we generalize these results to a larger class of projection operators.
Let B(H) be the set of bounded linear operators on a Hilbert space H. An operator P ∈ B(H) is called a projection operator (sometimes we will use the terms orthogonal projection operator and self-adjoint projection operator) on H iff P = P P † . We denote the set of projection operators on H by P(H) and the set of all subspaces of H by L(H).
It is easy to see that ∨S is the smallest subspace in H containing S. 2) If S ⊆ H, the orthogonal complement of S is defined as S ⊥ = {x ∈ H|x⊥s for all s ∈ S}. 3) If S is a collection of subsets of H, we write ∨ S∈S S = ∨(∪ S∈S S).
Definition 5: Let P ∈ P(H) and let K = image(P ) = {P x|x ∈ H}. We call P the projection of H onto K. Two projections P and Q onto K and L are orthogonal (denoted P ⊥Q) if P Q = 0. It is easy to verify that P Q = 0 ⇔ K⊥L ⇔ QP = 0. (Theorem 5B.9, [10]) Definition 6: Let P, Q ∈ P(H) with K = image(P ) and L = image(Q). Then The structure (P(H), , ⊥) is a logic with unit I H (identity map on H) and zero Z H (zero map on H) (Theorem 5B.18, [10]). This logic is called Projection Logic. Lemma 2 (Theorem 5B.18, [10]): The map P → image(P ) from P(H) to L(H) is a bijection that preserves order, orthogonality, meet(∧) and join(∨).

Lemma 4 ([10]):
If P, Q ∈ P(H), then Lemma 5: If P and Q are commutative operators, then the distributive law holds (and this law fails to hold for noncommutative operators). Also, in this case, Hence, image(P Q) = image(P ∧ Q) and by Lemma 2, P ∧ Q = P Q.

2) We have
where (a) follows from Lemma 3, (b) follows from Lemma 4 and (c) follows directly from definition of P ⊕ Q. 3)P = I − P follows directly from Definition 6.

4) We have
where (d) follows from Lemma 3 since P ⊕ Q and P ∧ Q are orthogonal ((P + Q − 2P Q)P Q = 0), (e) follows from Lemma 4, and (f ), (g) follows from the distributive laws. Hence, P ∨ Q = P + Q − P Q.
Next we define projection functions following [3].

Definition 7:
Given an arbitrary Boolean function f (v 1 , ...., v m ), we define the projection function f (P 1 , ..., P m ) in which v i in the Boolean function is replaced by P i , multiplication in the Boolean logic is replaced by the meet operation in the projection logic, summation in the Boolean logic (or the or function) is replaced by the join operation in the projection logic and the not operation in Boolean logic is replaced by the tilde (P ) operation in the projection logic.
As is standard when writing Boolean functions, we use xor (modulo 2 addition, represented by ⊕) in place of or, hence by above definition, we will replace the xor in the Boolean logic by the xor operation in the projection logic.
Theorem 1: If (P 1 , ..., P m ) are pairwise commutative projection operators of dimension 2 m−1 such that P 1 P 2 ..P m , P 1 P 2 ..P m , ...P 1P2 ..P m are all one-dimensional projection operators and H is of dimension 2 m , then P f = f (P 1 , ....P m ) is an orthogonal projection on a subspace of dimension T r(P f ) = wt(f ), where wt(f ) is the Hamming weight of the Boolean function f .
Proof: By definition of f (P 1 , ....P m ), we have a representation of P f in terms of meet, join and tilde operations in the corresponding projection logic. By Lemma 2, every function of projection operators in terms of meet, join and tilde will be present in the projection logic. Hence, P f is an orthogonal projection operator and this proves the first part of the theorem. Now, we will find the dimension of this projection operator. f as described in Section II.
If wt(f ) = M , then M terms of y i are 1 and the remaining terms are 0. Also, in this case, . We know by the statement of the theorem that the dimension of P c0 , and all these subspaces are orthogonal. Also, the minimum subspace containing all these operators is the whole Hilbert space. So, the dimension of P f will be the sum of dimensions of for all i (which is 1 when y i = 1, and 0 otherwise). Hence, the dimension of P f is M .
Theorem 1 is a generalization of the Theorem 1 of [3] because we consider any pairwise commutative projection operators, while in [3], a special case of commutative projection operators using Heisenberg-Weyl group was used. This special case is described in Section IV. Hence, to prove Theorem 1, we use abstract properties of projection logic [10] rather than the properties of a particular commutative subgroup.

IV. THE CONSTRUCTION OF COMMUTATIVE PROJECTION OPERATORS FROM THE HEISENBERG-WEYL GROUP
Let X, Y , and Z be the Pauli matrices, given by We form the Heisenberg-Weyl group (sometimes in the literature this group is referred to as an extraspecial 2-group or as the Pauli group) E m of order 4 m+1 , which is realized as the group of linear operators αE, α = ±1, ±i. (For a detailed description of the Heisenberg-Weyl group and its use to construct quantum codes see [7], [8].) Next we define the symplectic product of two vectors and the symplectic weight of a vector.

Definition 8:
The symplectic inner product of vectors (a, b), (a ′ , b ′ ) ∈ F 2m q is given by (1) symplectic weight of a vector (a, b) is the number of indices i at which either a i or b i is nonzero.
The center of the group E m is {±I 2 m , ±iI 2 m } and the quotient group E m is isomorphic to the binary vector space F 2m 2 . We associate with binary vectors (a, b) ∈ F 2m 2 operators E (a,b) defined by where Lemma 6: Lemma 7: Thus E (a,b) and E (a ′ ,b ′ ) commute iff (a, b) and (a ′ , b ′ ) are orthogonal with respect to the symplectic inner product (1).
We will now describe how to construct commutative projection operators. Take m linearly independent vectors y 1 , y 2 , ..., y m of length 2m bits with the property that the symplectic product between any pair is equal to zero. If we take P i = 1 2 (I + E yi ), then P 1 , ... P m satisfy all the properties of Theorem 1 and hence, f (P 1 , ...P m ) is an orthogonal projection operator [3].
Take y 1 , y 2 and y 3 as (1, 0, 0, 0, 1, 0) , (0, 1, 1, 1, 1, 0) and (0, 0, 1, 0, 1, 1) respectively which are linearly independent with all pairwise symplectic products equal to zero. Then P f = P 1 ⊕ P 1 P 2 ⊕ P 3 = P 1 + P 3 − 2P 1 P 3 − P 1 P 2 + 2P 1 P 2 P 3 where P i = 1 2 (I + E yi ), that is M )) quantum error correcting code is an Mdimensional subspace of C 2 k . The parameter k is the codelength and the parameter M is the dimension or the size of the code. Let Q be the quantum code, and P be the corresponding orthogonal projection operator on Q. (For a detailed description, see [4].) Definition 10: An error operator E is called detectable iff P EP = c E P , where c E is a constant that depends only on E.
Following [12], we restrict attention to the errors in the Heisenberg-Weyl group. Next, we define the minimum distance of the code.

Definition 11:
The minimum distance of Q is the maximum integer d such that any error E, with symplectic weight at most d − 1, is detectable.
The parameters of the quantum error correcting code are written ((k, M, d)) where the third parameter d is the minimum distance of Q. We say that a ((k, M, d)) quantum error correcting code exists if there exists a ((k, M )) quantum error correcting code with minimum distance ≥ d. We assume d ≥ 2 throughout the paper. We also focus on non-degenerate ((k, M, d)) codes, for which P EP = 0 for all errors E of symplectic weight ≤ d − 1, which is a sufficient condition for existence of the quantum code.
For any quantum code Q, we define the stabilizer H Q as where E k is the Heisenberg-Weyl group defined in Section IV. Then H Q is an abelian group and is isomorphic to GF(2) m , for some m. A quantum code is called additive or a stabilizer code if it is defined by its stabilizer H Q , i.e.
A quantum code is non-additive if it is not equivalent to an additive code [22].

VI. QUANTUM ERROR CORRECTING CODES WITH MINIMUM DISTANCE d
We use * to denote the standard binary inner product. The projection operator corresponding to the QECC is obtained as follows: (i) Construct the matrix A f as above.
(ii) Define k projection operators each of the form 1 2 (I +E y ) where y is a row of the matrix A f , with P k corresponding to the 1 st row, P k−1 corresponding to the 2 nd row and so on, so that P 1 corresponds to the last row. (iii) Transform the Boolean function f into the projection operator P f using Definition 7 where the commutative projection operators P 1 .... P k are determined by the matrix A f . Proof: Consider a Boolean function f (v) satisfying conditions 1) and 2). It follows easily from Section III and IV that P f constructed as above is an M -dimensional projection operator. It remains to prove that the minimum distance is at least d, so we need to show that P f ηP f η = 0 for any error η in E k with symplectic weight at most d − 1.
An error η in E k transforms the projection operator P f to P ′ f = ηP f η, and the condition P f ηP f η = 0 means that P ′ f is orthogonal to P f . Denote by η i the error represented by the binary 2k-tuple with entry 1 in position i and zeros elsewhere. We emphasize that the subscripts i in x i , η i and A j,i ((j, i) th entry in the matrix A f ) are read modulo 2k, so that x 2k+1 is just x 1 .
If A 1,k+1 = 0 then η 1 commutes with P k and η 1 P k η 1 = P k , and if A 1,k+1 = 1 then η 1 P k η 1 =P k . In general, if A k+1−j,k+i = 0 then η i P j η i = P j , and if A k+1−j,k+i = 1 then η i P j η i =P j . Let η i P j η i = Q i,j where Q i,j = P j orP j and observe that Q i,j = P j if and only if entry (k + 1 − j) of x k+i is zero. Then η i P f η i = f (Q i,1 , Q i,2 , ..., Q i,k ) and the entries of x k+i determine η i P f η i . In fact, this correspondence can easily be understood in terms of the fundamental correspondence between between Boolean functions and projection operators, since the operator η i P f η i corresponds to the Boolean function f (v ⊕ x k+i ).
When d = 2, we need to take care of all errors of symplectic weight 1 by showing P f η i P f η i = 0 and P f η i η i+k P f η i η i+k = 0. Applying the fundamental correspondence between Boolean functions and projection operators, this is equivalent to show- This follows from the assumption that x k+i and x k+i ⊕ x i are in the complementary set Cset f ..
In general we need to show that P f ηP f η = 0 for all errors η of symplectic weight at most d − 1. We write η = i∈A η i , apply the fundamental correspondence, and find that P f ηP f η corresponds to the Boolean function f (v ⊕ ( ⊕ i∈A x i+k )). By assumption, ⊕ i∈A x i+k is in the complementary set Cset f , so x i+k )) = 0 for all v, and hence P f ηP f η = 0.
Note that for M ≥ 1 this construction only gives ((k, M, d)) quantum error correcting codes for which the minimum distance d is at most k+3 2 . This is because any k+1 columns of the matrix A f are linearly dependent, which means that there is a 2k bit vector w of symplectic weight at most k+1 2 such that [x 1 , x 2 , ....x 2k ] * w T = 0, and the zero vector is never in Cset f . A ((k, M, d)) additive QECC exists when 1) M = 2 m for some m 2) There are 2k binary k-tuples x 1 , x 2 , ..., x 2k such that

Lemma 8:
Cset have pairwise symplectic product zero and are linearly independent.

Remark 1:
The projection operator corresponding to the QECC is k i=m+1 1 2 (I + E yi ) where y i is k + 1 − i th row of A f . The quantum code obtained in this way is that formed in the stabilizer framework using E y k , E y k−1 , ..., E ym+1 as the stabilizers of the code.
Proof: By Theorem 2 there exists a ((k, M, d))-QECC. The construction method of Theorem 2 gives the corresponding projection operator as P f = k i=m+1 P i = k i=m+1 1 2 (I + E yi ). Any vector in the code subspace is given by |x >= P f |u > for some |u >∈ H. Since E yi and E yj are commutative, we have E yi |x >= |x > for m < i ≤ k. Hence, E y k , E y k−1 , ..., E ym+1 are the stabilizers of the quantum code and the quantum code is additive.
Remark 2: If the boolean function can be represented as a single monomial, it gives an additive code. The converse is not true in general; see for example, [22], where it is shown that every ((4, 4, 2)) code is equivalent to an additive code.
We can also see that the second property is satisfied, so we have constructed a ((2m, 4 m−1 , 2)) QECC that is not additive.
Lemma 9: 1) If there exists a ((k, M, 2)) QECC, then there exists a ((k + 2, 4M, 2)) QECC determined by  (1, 0, x 2k )). All the columns and the sum of columns i and i + k are in Cset f ′ . The symplectic product of any two rows is zero and all the rows are linearly independent, since this was true for A f = (x 1 , x 2 , . . . , x 2k ) 2) Given this choice of f ′ (v), we have Cset f ′ ⊇ Cset f , and this means that the same matrix A f ′ = A f will satisfy all the earlier properties.
Consider the Boolean function f It is a function of 2m+1 variables with weight 3 × 2 2m−3 .
Example 9: The perfect ((5, 2, 3)) additive code of R. Laflamme et al. [17] can be obtained by the above approach. it is easy to see that all rows are linearly independent, and that the symplectic inner product of any two rows is zero. Note that the stabilizers corresponding to the code are ZXXZI, IZXXZ, ZIZXX, and XZIZX.

VII. OPERATOR QUANTUM ERROR CORRECTION (OQEC)
The theory of operator quantum error correction [16] uses the framework of noiseless subsystems to improve the performance of decoding algorithms which might help improve the threshold for fault-tolerant quantum computation. It requires a fixed partition of the systems Hilbert space H = A ⊗ B ⊕ C ⊥ . Information is encoded on the A subsystem; the logical quantum state ρ A ∈ B A is encoded as ρ A ⊗ ρ B ⊕ 0 C ⊥ with an arbitrary ρ B ∈ B B (where B A and B B are the sets of all endomorphisms on subsystems A and B respectively). We say that the error E is correctable on subsystem A (called the logical subsystem) when there exists a physical map R that reverses its action, up to a transformation on the B subsystem (called the Gauge subsystem). In other words, this error correcting procedure may induce some nontrivial action on the B subsystem in the process of restoring information encoded in the A subsystem. This leads to recovery routines which explicitly make use of the subsystem structure [6] [24]. In the case of standard quantum error correcting codes, the dimension of B is 1. A ((k, M, N, d))-OQEC is defined as a OQEC in C 2 k with M and N as the dimension of the logical and gauge subsystems.


extfieldsofstudy: Mathematics Biology Physics

text:
Canalization and Symmetry in Boolean Models for Genetic Regulatory Networks

Canalization of genetic regulatory networks has been argued to be favored by evolutionary processes due to the stability that it can confer to phenotype expression. We explore whether a significant amount of canalization and partial canalization can arise in purely random networks in the absence of evolutionary pressures. We use a mapping of the Boolean functions in the Kauffman N-K model for genetic regulatory networks onto a k-dimensional Ising hypercube to show that the functions can be divided into different classes strictly due to geometrical constraints. The classes can be counted and their properties determined using results from group theory and isomer chemistry. We demonstrate that partially canalized functions completely dominate all possible Boolean functions, particularly for higher k. This indicates that partial canalization is extremely common, even in randomly chosen networks, and has implications for how much information can be obtained in experiments on native state genetic regulatory networks.

Introduction
To preserve the identity of a species, biological organisms must be capable of maintaining relatively stable phenotype expression in the face of a variety of environmental factors and a certain level of genetic randomness. Experimental observations have shown that certain developmental traits appear to control the expression of other traits. Waddington [1] termed the control of one trait by another "canalization," a name derived from the analogy that the developmental pathway of the organism is like one particular canal in a floodplain, and the further development of the organism is completely constrained by that canal. Canalization produces robustness because it suppresses those changes in phenotype expression that would require development to deviate from the canalized pathway. For this reason it has been suggested that organisms evolve to be canalized.
The significance of canalization and how it might evolve remains a subject of debate [2]. Since canalization suppresses the expression of genetic variability, experimental detection of the existence of a canalized trait generally involves perturbing the organism out of the canalized state [3]. There is good evidence for the existence of canalization in a variety of organisms [4,5,6]; however, the microscopic mechanism for canalization is not well established. Presumably canalization is produced genetically by the complex interactions between genes known as the genetic regulatory network (GRN). As we shall see, however, a certain amount of canalization is expected to appear in GRNs even in the absence of an evolutionary preference for canalization. An open question is whether or not real GRNs contain more canalization than would be expected from a random graph, which could indicate that evolution favors canalization.
Genetic regulatory networks have been proposed as the mechanisms through which identical genetic information is expressed as different cell types within the same organism, and they can also control distinct stages in the life cycle of an individual cell. Depending on the conditions experienced by a given cell and the regulatory interactions between genes, at any moment a distinct subset of all possible genes are activated. The state or temporal pattern of expression produces particular cell types. Organisms with larger numbers of genes have a larger number of potentially realizable cell types. There has been a recent dramatic increase in the amount of experimental information available on the structure of genetic regulatory networks in a range of organisms, including E. coli [7], budding yeast S. cervisiae [8,9], Drosophila species [10], Xenopus [11], and the embryo of the sea urchin S. purpuratus [12]. In the simplest representation, the nodes of the network are genes and the links between genes describe their interactions. Generally, the interactions are directional, so that the expression of gene A may depend on, that is "listen to," the expression of gene B, but the expression of gene B doesn't necessarily depend directly on the expression of gene A. The connectivity of a gene indicates how many other genes it "listens to" when determining whether to be in an active or inactive state. Analysis of the connectivity of E. coli [13,14] and other GRNs shows a broad distribution of connectivity among the genes, with a significant amount Table 1. An example of a k = 3 Boolean function.
of negative autoregulation. In the context of canalization, several questions arise. What types of structures in a GRN produce canalization of a trait? Do these structures arise randomly, or do they only appear because of a special evolutionary preference? How significant is canalization on the scale of the entire regulatory network?
The easiest way to approach such questions is through a simplified model for a genetic regulatory network such as the Kauffman N-K model [15], which represents the GRN as a random Boolean network. The N-K model has been studied extensively [16,17,18,19,20,21,22,23]. Certain features of real GRNs, including the ability of a single network to produce multiple cell types (which appear as multiple attractors for the network), are captured by the N-K model. In this model, each gene is represented by a single binary element which can be either on or off in the state 1 or 0. Every gene receives input from a fixed set of k other genes that are randomly chosen when the network is constructed, where k = K. Depending on the states of its input genes, a given gene determines whether to express the state 1 or 0 according to a randomly chosen Boolean function of k variables. An example of a Boolean function for k = 3 is given in table 1. The value of k may vary from gene to gene. The system evolves in discrete time steps and all genes update their states simultaneously. The entire network eventually settles into an attractor cycle which produces a specific series of network states as a function of time. The initial conditions of the states of the genes in the network determine which of the available attractors the network will reach. The different attractors are interpreted as representing different cell types expressed by a given set of genes.
A gene with connectivity k employs one of the 2 2 k possible Boolean functions to determine its response to its k inputs. Canalization occurs in a Boolean function if the output of the gene is fixed by a particular value of at least one of its input genes, regardless of the values of any other inputs to that gene. In this case the input that fixes the output of the regulated gene is a canalizing input. Note that one value of an input gene, say value 0, can be a canalizing input even if the other value 1 from the same input gene is not canalizing. Canalization also occurs in a Boolean function if particular values of two or more inputs together suffice to guarantee the next state of the regulated gene, regardless of the values of any other inputs to the gene. In this case, the inputs that together fix the output of the regulated gene are said to be collectively canalizing inputs. How canalizing a particular Boolean function is can be quantified by the set of numbers P n , n = 0, 1, . . . , k − 1, which are the fraction of sets of n individual input values that are canalizing or collectively canalizing. Note that Boolean functions with P 0 = 1 have a fixed output state regardless of their inputs. Boolean functions can also be characterized by their internal homogeneity p which is defined as the fraction of 1s or 0s, whichever is larger, output by the function due to all of the possible sets of input [24].
A consequence of canalization is that some of the interactions between genes may become irrelevant. As an extreme example, if the Boolean function of a particular gene has P 0 = 1, this gene will be insensitive to the state of the rest of the network and its interactions with its input genes are irrelevant. The number of canalizing functions as a function of k was derived recently in Ref. [25]. Although the behaviour of canalizing functions would certainly contribute to the stability of a network that is subjected to random perturbations, such an extreme behaviour has a detrimental effect on the ability of the network to respond to changing conditions. In contrast, other Boolean functions successfully maintain a degree of stability while retaining the ability to change. For these Boolean functions, which we will refer to as "partially canalizing," the gene may ignore one or more of its inputs under certain conditions. In some cases, the gene completely ignores n inputs at all times, so that its effective connectivity is k ef f = k − n. In other cases, if a particular input has the value 1, for example, the gene ignores the remaining inputs, but if that same input has the value 0, the gene listens to its other inputs. Here, the effective connectivity of the gene depends on the current state of the network. More complex categories are also possible, such as the nested canalizing functions proposed by Kauffman [19,26]. Since the fraction of canalizing functions drops rapidly with k, as shown in Ref. [27], it has been assumed that canalization plays a less important role at high k compared to small k [28]. The class of partially canalizing functions is considerably larger than the class of canalizing functions; however, is it large enough to dominate all classes of functions? As we will show below on mathematical grounds, the partially canalizing functions completely dominate the class of all possible Boolean functions as k increases, so that the emergence of canalization is essentially unavoidable in a complex network.

Results
Our approach is to examine the properties of individual Boolean functions and to determine the amount of canalization expected from a random sample of functions. Since the number of possible Boolean functions explodes combinatorially with k, we employ powerful techniques from group theory and isomer chemistry to classify the various functions and help obtain their properties. We provide results through k = 5 with these methods. The techniques can be applied readily to higher k, but become increasingly complicated. Therefore, to find results for larger k through k = 8 we employ For small enough k, the canalization properties of the functions can be obtained directly from inspection. When there are two inputs for each gene, k = 2, as shown in table 2, there are only 16 possible functions which fall into four classes: fixed (or completely canalizing) with P 0 = 1; sensitive to both inputs with P 0 = 0 and P 1 = 0; and the partially canalizing cases with P 0 = 0 and P 1 = 1/2: sensitive to only one input; sensitive to one or two inputs depending on the value of one input.
Inspection becomes a less viable option as k increases. In a simulation study of the evolution properties of the different Boolean functions, Bassler et al. [29] observed that functions with k = 3 inputs fell into 14 distinct classes. In their study all of the functions that were members of the same class evolved, on average, with equal probability. Upon examining representative functions from each class, they categorized the functions according to their canalization properties P n . The triple of numbers P 0 , P 1 , and P 2 possible for k = 3 was nearly enough, but not quite enough, to distinctly identify each class of function. Whether the function was symmetric about its midpoint also needed to be considered in defining the classes. These observations about the structures of the functions belonging to each class were essentially empirical. Class membership could be important for determining the properties of real networks since we expect that on average all functions in the same class will evolve with equal probability.
Here, we demonstrate that there is a fundamental geometric reason for the existence of distinct function classes. In the N-K model, a given function is normally represented by a Boolean string of numbers, such as 1001, of length 2 2 k . Comparing different functions by inspection amounts to comparing strings of numbers with each other. Instead of using this representation, we consider an alternative, equivalent representation of each function in the form of a unit k-dimensional Ising hypercube. Each axis of the kdimensional hypercube (or simply a k-hypercube) represents one of the k input variables. The coordinates on a given axis indicate the state of the corresponding input variable. Each vertex of the k-hypercube represents an output state of the gene. In figure 1 we illustrate the mapping of the input states onto a square and cube for k = 2 and k = 3, respectively. The output state of the gene corresponding to an input represented by a particular vertex can be indicated by colouring the vertex white or black to represent the values 0 or 1. It is important to note that this system obeys parity symmetry:  replacing all 0's with 1's and vice versa results in the same canalization properties for the function.
With this hypercube mapping, it becomes clear that functions which belong to the same class have the geometric property that they can be rotated into each other by symmetry operations on the k-hypercube plus parity. In mathematical terms, the classes that were identified empirically in Ref. [29] are group orbits. We illustrate the mapping for the sixteen k = 2 functions in figure 2, where the rotational plus parity symmetries of the functions belonging to each of the four classes are obvious. In figure  3 we illustrate one representative cube for each of the 14 function classes in k = 3. The remaining functions in each class are obtained by applying all possible rotations plus parity to the cube. In the hypercube representation, the canalization properties of a Boolean function correspond to the fraction of homogeneous hypersurfaces. That is, for a Boolean function with k inputs P n is proportional to the fraction of the k − n dimensional hypersurfaces that have all vertices the same. For the two classes with P 1 = 1/2 in figure 2, 2 of the 4 one-dimensional sides are uniformly coloured.
We can now employ results from group theory to obtain information about the class structure of functions at all values of k, not merely those values of k which are small enough to permit direct inspection of all functions. The total number of classes for a given k can be obtained by an application of the orbit-counting theorem, Here, the symmetry group G of the set X contains |G| symmetry operators g, which together include all transformations of the hypercube onto itself. The set of elements in X that are left invariant by g is denoted as X g . In order to find X g , note that the mapping of all k-hypercube vertices onto themselves by a given symmetry operator g can be written as a permutation of the vertex numbers. As a result, each operator g can be expressed in terms of its cycle structure x b 1 1 x b 2 2 ...x bn n , where n = k. This notation indicates that g contains b 1 cycles of length 1, b 2 cycles of length 2, ... b n cycles of length n [31]. For example, the k = 2 permutation (14)(2)(3) has the cycle representation x 2 1 x 2 since it has 2 cycles of length 1 and a single cycle of length 2. To apply the orbit-counting theorem, we must first construct all of the operators of our group, sum the number of functions left invariant by each operator (the fixed points of that operator), and divide by the total number of operators.
The number of functions in, or size of, a class, is given by the number of elements in the group divided by the number of elements in the isometry group of the functions in the class. The isometry group of a class is the subgroup of the full group that describes the symmetry of a function in the class. Note that the particular isometry group will vary from function to function in the class, but the size of the isometry group will remain invariant.
First, consider the number of symmetry operations |G| in our group, which is the khypercube crossed with parity. The symmetry group for the k-hypercube is isomorphic to the hyperoctahedral group O n with n = k, which has n!2 n symmetry transformations [30]. As an example, there are eight operators for the k = 2 square. There is one operator with cycle structure  (3) and (23)(1)(4). When these operators are combined with parity, which doubles the number of symmetry operators, we obtain Table 3. Cycle polynomials for k = 1 through 5 and the number of classes P G for each k.
k Cycle polynomial a total of |G| = 16 operators. For each operator without parity, the number of functions left invariant is equal to 2 Nc , where N c = k i=1 b i is the total number of cycles in the operator. Parity must be treated separately; no functions are left invariant by the parity operator with any k-hypercube operator containing at least one cycle of length 1. Thus there are 2 Np functions left invariant for the eight operators which include parity, where and Θ is the Heaviside step function. Applying the orbitcounting theorem produces the correct number of classes, but only if parity is included. In the case of k = 2 without parity, P G = (1/8)(2 4 + 2(2) + 3(2 2 ) + 2(2 3 )) = 6 classes. Including parity gives P G = (1/16)(2 4 + 2(2) + 3(2 2 ) + 2(2 3 ) + 2(2) + 3(2 2 )) = 4 classes, which is the correct result. We now face the task of identifying all operators g for the k-hypercube. This can be performed by simple inspection in k = 2 and k = 3, but it becomes more complicated to keep track of higher-dimensional rotation symmetries. Fortunately, this problem was solved in the middle of the last century, when Harrison [32] derived a formula that produces the complete cycle representation for all k in the group of interest to us, called the "Zyklenzeiger" in the notation of Ref. [32]. We have used this formula to obtain the cycle representations through k = 5, shown in table 3. The number of classes for each k, P G , is also listed in table 3. Clearly, obtaining the properties of the classes by inspection is not feasible by the time k = 5.
What we are interested in is not simply how many different classes are present, but also the size of each class, and the structure, particularly the canalization properties, of the functions belonging to them. For example, how many classes are there which have the same internal inhomogeneity p? To find this, we use an application of Pólya's theorem [33] which is frequently used in isomer chemistry [34]. In isomer chemistry, for a molecule composed of exactly two different types of atoms, the terms A and B can be used to represent the different types of atoms. In our case, A and B represent 0 and 1, such that either A = 0 and B = 1, or B = 0 and A = 1. Using the generating polynomial, substitute in a term of the form A a + B a for each x a . Divide the result by the total number of operators, including parity. Then, drop all terms in the result where the exponent of B exceeds that of A, as these terms are already accounted for by parity.  Table 5. Class structure for k = 3.
The multiplicity of each term indicates how many of the classes are of that form. For example, representatives of classes of the form A 2 B 2 are 1010 and 1100. This gives us the desired result of how many classes N h there are for each value of the internal homogeneity p. Since we also know that the total number N f (m, n) of functions of the form A m B n is simply N f (m, n) = (2 − δ m,n )(2 k )!/(m!n!), where m + n = 2 k , we can estimate the number of functions in each class by the average size of a class, S c = N f /N h . The class structure and average class size for k = 1 through 5 are listed in tables 4 through 7. We note that the actual size S c of each class is given by the number of operators that preserve the symmetry of that particular function class. Thus, as discussed earlier, the maximum class size S max c is equal to the total number of operators S max c = 16 for k = 2, 96 for k = 3, 768 for k = 4, and 7680 for k = 5. This is consistent with the average class sizes which we obtain.
We also note that isomer chemistry provides a simple means for determining whether two randomly selected functions belong to the same class. Construct the adjacency matrix for the k-hypercube. Along the diagonal, place the values A or B corresponding to the colours of the vertices of one of the functions under consideration, and then find the determinant of the resulting matrix. Each function class has a unique determinant, so performing this procedure on both functions provides an immediate test of whether the two functions fall into the same class.
To show that canalization remains important even as k increases, we measured the average number of homogeneous d-dimensional sides present in a series of randomly generated functions for different k. We denote the number of d-dimensional homogeneous sides (which produce canalization) that a k-dimensional Boolean function has as C(d, k). The total number of d-dimensional sides is Note that the canalization properties P n discussed in the beginning of this paper are related to C(d, k) as P n = C(k − n, k)/N k−n (k). In figure 4  The mapping of Boolean functions onto k-hypercubes we have described here provides a means of constructing k + 1 functions recursively from pairs of k functions. A k + 1 function can be composed by stacking together two k functions. Depending on the symmetry properties of the two k functions chosen, there may be only one possible class of k + 1 functions that can be constructed from those k functions, or there may be several classes that depend on the relative orientation of the k functions when they are stacked together. This allows us to bound the amount of canalization present. When we assemble a k + 1 function out of two k functions, we must have The lower bound is obtained from the fact that the value of C(d, k + 1) must be at least as large as the sum of the values C i (d, k) of the two functions that have been combined. This is simply a consequence of the fact that homogeneous d-dimensional sides cannot be destroyed as a result of combining two functions. The new sides that are added when the functions are joined may or may not be homogeneous, depending on which two k functions are combined and how they are oriented with respect to each other. It is possible that none of the new sides would be homogeneous, in which case the lower limit of Eq. (4) would apply. The internal homogeneity p(k + 1) of the composite function is given simply by

Discussion
Now that we have used the mapping of the functions onto k-hypercubes in order to obtain information about the class structures of the functions for several values of k, we can make some observations regarding how prevalent partial canalization is among all possible functions. Previous estimates of the fraction of canalizing functions indicated that canalization was of less and less importance as k increased. These estimates used a very narrow definition of canalization, however. Rather than counting the number of partially canalizing functions, consider the number of completely uncanalizing functions. These functions have the property that they are sensitive to all values of all inputs. There are exactly 2 such functions for each k, regardless of k. All of the remaining functions are at least partially canalizing. This means that partial canalization completely dominates the classes of functions, especially as k increases. The rampant occurrence of partial canalization has important implications for recent work on mapping of genetic regulatory networks. The experiments typically map only those connections between genes which are active in the native state of the organism. Here, "active" means that a change in one gene directly affects the second gene. This technique will not detect many of the partially canalizing interactions that could exist between genes. In the case where the partial canalization is of the form that a gene completely ignores one or more of its inputs, the actual value of k for that gene is larger than the apparent value of k. This could potentially impact the distributions of k that have been extracted from experimental measurements. A far more dangerous case is a partially canalizing interaction between genes in which a gene ignores one or more inputs when a canalizing input has a value of 1 (for example), but responds to the other inputs when that same canalizing input has a value of 0. If the gene ignores its inputs in the native state, the connection between that gene and its ignored inputs will not be detected experimentally. Suppose that the canalizing gene is identified as causing a disease state. Consultation of the experimentally determined genetic network map indicates that this gene does not appear to control anything else of importance. If, however, the canalizing gene is treated and switches to the state opposite from its canalizing value, the gene that received the canalizing input will suddenly start to respond to the values of its other inputs. This could result in unexpected side effects or worse effects. Thus, from a purely combinatorial point of view, it is important to consider all possible interactions between genes, and not merely those which are expressed in the native state.
The natural predominance of canalization as k increases suggests that the canalization observed experimentally could be due simply to the high fraction of the available Boolean functions which are canalizing, rather than evolutionary pressure to develop canalizing functions. It is, however, unclear how much canalization is present in real genetic regulatory networks, as discussed in Ref. [35]. It is possible that there is in fact a special evolutionary preference for canalization, which could result in real networks having even higher levels of canalization than would be expected from random selection at increasing k. In order to answer this question it would be necessary to measure the excess canalization, which is the difference between the P n s observed in real networks and that in random networks [36]. The existing experimental data on genetic regulatory networks is not extensive enough to determine whether an interaction between genes is canalizing or partially canalizing. As noted above, the difference between the two types of interactions can become important when the network is perturbed away from its native state. More experimental work is needed in order to determine the prevalence of canalization and/or partial canalization in actual genetic regulatory networks. The Boolean models can offer guidance in determining how likely it would be to observe any type of canalization in a random network.

Conclusion
In conclusion, we have used a mapping of the Boolean functions in the Kauffman model for genetic regulatory networks onto a k−hypercube to obtain information about the classes into which the functions can be divided. These classes arise due to geometrical constraints, and can be constructed by applying all possible rotations of the k−hypercube plus parity to each function. The classes can be counted and their properties determined using results from group theory and isomer chemistry. We emphasize that partially canalizing functions completely dominate all possible functions, particularly for higher k. This indicates that partial canalization should be extremely common, even in a randomly chosen network, and has implications for how much information can be obtained in experiments on native state genetic regulatory networks.