MINIMAL PERCOLATING SETS IN BOOTSTRAP PERCOLATION

合集下载

Discovering episodes with compact minimal windows-DMKD

The common wisdom is that ﬁnding frequent patterns is not enough. Discovering frequent patterns with high threshold will result to trivial patterns, omitting many interesting patterns, while using a low threshold will result in a pattern explosion. This phenomenon has led to many ranking methods for itemsets, the most well-studied pattern type. Unlike for itemsets, ranking episodes is heavily under-developed. Existing statistical approaches for ranking episodes are mostly based on the number of ﬁxedsize windows (see more detailed discussion in Sect. 6). However, a natural way of measuring the goodness of an episode is the average length of its instances—a good episode should have compact minimal windows. Hence, our goal and contribution is a measure based directly on the average length of minimal windows.

最近鲁棒优化进展Recent Advances in Robust Optimization and Robustness An Overview

Recent Advances in Robust Optimization and Robustness:An OverviewVirginie Gabrel∗and C´e cile Murat†and Aur´e lie Thiele‡July2012AbstractThis paper provides an overview of developments in robust optimization and robustness published in the aca-demic literature over the pastﬁve years.1IntroductionThis review focuses on papers identiﬁed by Web of Science as having been published since2007(included),be-longing to the area of Operations Research and Management Science,and having‘robust’and‘optimization’in their title.There were exactly100such papers as of June20,2012.We have completed this list by considering 726works indexed by Web of Science that had either robustness(for80of them)or robust(for646)in their title and belonged to the Operations Research and Management Science topic area.We also identiﬁed34PhD disserta-tions dated from the lastﬁve years with‘robust’in their title and belonging to the areas of operations research or management.Among those we have chosen to focus on the works with a primary focus on management science rather than system design or optimal control,which are broadﬁelds that would deserve a review paper of their own, and papers that could be of interest to a large segment of the robust optimization research community.We feel it is important to include PhD dissertations to identify these recent graduates as the new generation trained in robust optimization and robustness analysis,whether they have remained in academia or joined industry.We have also added a few not-yet-published preprints to capture ongoing research efforts.While many additional works would have deserved inclusion,we feel that the works selected give an informative and comprehensive view of the state of robustness and robust optimization to date in the context of operations research and management science.∗Universit´e Paris-Dauphine,LAMSADE,Place du Mar´e chal de Lattre de Tassigny,F-75775Paris Cedex16,France gabrel@lamsade.dauphine.fr Corresponding author†Universit´e Paris-Dauphine,LAMSADE,Place du Mar´e chal de Lattre de Tassigny,F-75775Paris Cedex16,France mu-rat@lamsade.dauphine.fr‡Lehigh University,Industrial and Systems Engineering Department,200W Packer Ave Bethlehem PA18015,USA aure-lie.thiele@2Theory of Robust Optimization and Robustness2.1Deﬁnitions and BasicsThe term“robust optimization”has come to encompass several approaches to protecting the decision-maker against parameter ambiguity and stochastic uncertainty.At a high level,the manager must determine what it means for him to have a robust solution:is it a solution whose feasibility must be guaranteed for any realization of the uncertain parameters?or whose objective value must be guaranteed?or whose distance to optimality must be guaranteed? The main paradigm relies on worst-case analysis:a solution is evaluated using the realization of the uncertainty that is most unfavorable.The way to compute the worst case is also open to debate:should it use aﬁnite number of scenarios,such as historical data,or continuous,convex uncertainty sets,such as polyhedra or ellipsoids?The answers to these questions will determine the formulation and the type of the robust counterpart.Issues of over-conservatism are paramount in robust optimization,where the uncertain parameter set over which the worst case is computed should be chosen to achieve a trade-off between system performance and protection against uncertainty,i.e.,neither too small nor too large.2.2Static Robust OptimizationIn this framework,the manager must take a decision in the presence of uncertainty and no recourse action will be possible once uncertainty has been realized.It is then necessary to distinguish between two types of uncertainty: uncertainty on the feasibility of the solution and uncertainty on its objective value.Indeed,the decision maker generally has different attitudes with respect to infeasibility and sub-optimality,which justiﬁes analyzing these two settings separately.2.2.1Uncertainty on feasibilityWhen uncertainty affects the feasibility of a solution,robust optimization seeks to obtain a solution that will be feasible for any realization taken by the unknown coefﬁcients;however,complete protection from adverse realiza-tions often comes at the expense of a severe deterioration in the objective.This extreme approach can be justiﬁed in some engineering applications of robustness,such as robust control theory,but is less advisable in operations research,where adverse events such as low customer demand do not produce the high-proﬁle repercussions that engineering failures–such as a doomed satellite launch or a destroyed unmanned robot–can have.To make the robust methodology appealing to business practitioners,robust optimization thus focuses on obtaining a solution that will be feasible for any realization taken by the unknown coefﬁcients within a smaller,“realistic”set,called the uncertainty set,which is centered around the nominal values of the uncertain parameters.The goal becomes to optimize the objective,over the set of solutions that are feasible for all coefﬁcient values in the uncertainty set.The speciﬁc choice of the set plays an important role in ensuring computational tractability of the robust problem and limiting deterioration of the objective at optimality,and must be thought through carefully by the decision maker.A large branch of robust optimization focuses on worst-case optimization over a convex uncertainty set.The reader is referred to Bertsimas et al.(2011a)and Ben-Tal and Nemirovski(2008)for comprehensive surveys of robust optimization and to Ben-Tal et al.(2009)for a book treatment of the topic.2.2.2Uncertainty on objective valueWhen uncertainty affects the optimality of a solution,robust optimization seeks to obtain a solution that performs well for any realization taken by the unknown coefﬁcients.While a common criterion is to optimize the worst-case objective,some studies have investigated other robustness measures.Roy(2010)proposes a new robustness criterion that holds great appeal for the manager due to its simplicity of use and practical relevance.This framework,called bw-robustness,allows the decision-maker to identify a solution which guarantees an objective value,in a maximization problem,of at least w in all scenarios,and maximizes the probability of reaching a target value of b(b>w).Gabrel et al.(2011)extend this criterion from aﬁnite set of scenarios to the case of an uncertainty set modeled using intervals.Kalai et al.(2012)suggest another criterion called lexicographicα-robustness,also deﬁned over aﬁnite set of scenarios for the uncertain parameters,which mitigates the primary role of the worst-case scenario in deﬁning the solution.Thiele(2010)discusses over-conservatism in robust linear optimization with cost uncertainty.Gancarova and Todd(2012)studies the loss in objective value when an inaccurate objective is optimized instead of the true one, and shows that on average this loss is very small,for an arbitrary compact feasible region.In combinatorial optimization,Morrison(2010)develops a framework of robustness based on persistence(of decisions)using the Dempster-Shafer theory as an evidence of robustness and applies it to portfolio tracking and sensor placement.2.2.3DualitySince duality has been shown to play a key role in the tractability of robust optimization(see for instance Bertsimas et al.(2011a)),it is natural to ask how duality and robust optimization are connected.Beck and Ben-Tal(2009) shows that primal worst is equal to dual best.The relationship between robustness and duality is also explored in Gabrel and Murat(2010)when the right-hand sides of the constraints are uncertain and the uncertainty sets are represented using intervals,with a focus on establishing the relationships between linear programs with uncertain right hand sides and linear programs with uncertain objective coefﬁcients using duality theory.This avenue of research is further explored in Gabrel et al.(2010)and Remli(2011).2.3Multi-Stage Decision-MakingMost early work on robust optimization focused on static decision-making:the manager decided at once of the values taken by all decision variables and,if the problem allowed for multiple decision stages as uncertainty was realized,the stages were incorporated by re-solving the multi-stage problem as time went by and implementing only the decisions related to the current stage.As theﬁeld of static robust optimization matured,incorporating–ina tractable manner–the information revealed over time directly into the modeling framework became a major area of research.2.3.1Optimal and Approximate PoliciesA work going in that direction is Bertsimas et al.(2010a),which establishes the optimality of policies afﬁne in the uncertainty for one-dimensional robust optimization problems with convex state costs and linear control costs.Chen et al.(2007)also suggests a tractable approximation for a class of multistage chance-constrained linear program-ming problems,which converts the original formulation into a second-order cone programming problem.Chen and Zhang(2009)propose an extension of the Afﬁnely Adjustable Robust Counterpart framework described in Ben-Tal et al.(2009)and argue that its potential is well beyond what has been in the literature so far.2.3.2Two stagesBecause of the difﬁculty in incorporating multiple stages in robust optimization,many theoretical works have focused on two stages.Regarding two-stage problems,Thiele et al.(2009)presents a cutting-plane method based on Kelley’s algorithm for solving convex adjustable robust optimization problems,while Terry(2009)provides in addition preliminary results on the conditioning of a robust linear program and of an equivalent second-order cone program.Assavapokee et al.(2008a)and Assavapokee et al.(2008b)develop tractable algorithms in the case of robust two-stage problems where the worst-case regret is minimized,in the case of interval-based uncertainty and scenario-based uncertainty,respectively,while Minoux(2011)provides complexity results for the two-stage robust linear problem with right-hand-side uncertainty.2.4Connection with Stochastic OptimizationAn early stream in robust optimization modeled stochastic variables as uncertain parameters belonging to a known uncertainty set,to which robust optimization techniques were then applied.An advantage of this method was to yield approaches to decision-making under uncertainty that were of a level of complexity similar to that of their deterministic counterparts,and did not suffer from the curse of dimensionality that afﬂicts stochastic and dynamic programming.Researchers are now making renewed efforts to connect the robust optimization and stochastic opti-mization paradigms,for instance quantifying the performance of the robust optimization solution in the stochastic world.The topic of robust optimization in the context of uncertain probability distributions,i.e.,in the stochastic framework itself,is also being revisited.2.4.1Bridging the Robust and Stochastic WorldsBertsimas and Goyal(2010)investigates the performance of static robust solutions in two-stage stochastic and adaptive optimization problems.The authors show that static robust solutions are good-quality solutions to the adaptive problem under a broad set of assumptions.They provide bounds on the ratio of the cost of the optimal static robust solution to the optimal expected cost in the stochastic problem,called the stochasticity gap,and onthe ratio of the cost of the optimal static robust solution to the optimal cost in the two-stage adaptable problem, called the adaptability gap.Chen et al.(2007),mentioned earlier,also provides a robust optimization perspective to stochastic programming.Bertsimas et al.(2011a)investigates the role of geometric properties of uncertainty sets, such as symmetry,in the power ofﬁnite adaptability in multistage stochastic and adaptive optimization.Duzgun(2012)bridges descriptions of uncertainty based on stochastic and robust optimization by considering multiple ranges for each uncertain parameter and setting the maximum number of parameters that can fall within each range.The corresponding optimization problem can be reformulated in a tractable manner using the total unimodularity of the feasible set and allows for aﬁner description of uncertainty while preserving tractability.It also studies the formulations that arise in robust binary optimization with uncertain objective coefﬁcients using the Bernstein approximation to chance constraints described in Ben-Tal et al.(2009),and shows that the robust optimization problems are deterministic problems for modiﬁed values of the coefﬁcients.While many results bridging the robust and stochastic worlds focus on giving probabilistic guarantees for the solutions generated by the robust optimization models,Manuja(2008)proposes a formulation for robust linear programming problems that allows the decision-maker to control both the probability and the expected value of constraint violation.Bandi and Bertsimas(2012)propose a new approach to analyze stochastic systems based on robust optimiza-tion.The key idea is to replace the Kolmogorov axioms and the concept of random variables as primitives of probability theory,with uncertainty sets that are derived from some of the asymptotic implications of probability theory like the central limit theorem.The authors show that the performance analysis questions become highly structured optimization problems for which there exist efﬁcient algorithms that are capable of solving problems in high dimensions.They also demonstrate that the proposed approach achieves computationally tractable methods for(a)analyzing queueing networks,(b)designing multi-item,multi-bidder auctions with budget constraints,and (c)pricing multi-dimensional options.2.4.2Distributionally Robust OptimizationBen-Tal et al.(2010)considers the optimization of a worst-case expected-value criterion,where the worst case is computed over all probability distributions within a set.The contribution of the work is to deﬁne a notion of robustness that allows for different guarantees for different subsets of probability measures.The concept of distributional robustness is also explored in Goh and Sim(2010),with an emphasis on linear and piecewise-linear decision rules to reformulate the original problem in aﬂexible manner using expected-value terms.Xu et al.(2012) also investigates probabilistic interpretations of robust optimization.A related area of study is worst-case optimization with partial information on the moments of distributions.In particular,Popescu(2007)analyzes robust solutions to a certain class of stochastic optimization problems,using mean-covariance information about the distributions underlying the uncertain parameters.The author connects the problem for a broad class of objective functions to a univariate mean-variance robust objective and,subsequently, to a(deterministic)parametric quadratic programming problem.The reader is referred to Doan(2010)for a moment-based uncertainty model for stochastic optimization prob-lems,which addresses the ambiguity of probability distributions of random parameters with a minimax decision rule,and a comparison with data-driven approaches.Distributionally robust optimization in the context of data-driven problems is the focus of Delage(2009),which uses observed data to deﬁne a”well structured”set of dis-tributions that is guaranteed with high probability to contain the distribution from which the samples were drawn. Zymler et al.(2012a)develop tractable semideﬁnite programming(SDP)based approximations for distributionally robust individual and joint chance constraints,assuming that only theﬁrst-and second-order moments as well as the support of the uncertain parameters are given.Becker(2011)studies the distributionally robust optimization problem with known mean,covariance and support and develops a decomposition method for this family of prob-lems which recursively derives sub-policies along projected dimensions of uncertainty while providing a sequence of bounds on the value of the derived policy.Robust linear optimization using distributional information is further studied in Kang(2008).Further,Delage and Ye(2010)investigates distributional robustness with moment uncertainty.Speciﬁcally,uncertainty affects the problem both in terms of the distribution and of its moments.The authors show that the resulting problems can be solved efﬁciently and prove that the solutions exhibit,with high probability,best worst-case performance over a set of distributions.Bertsimas et al.(2010)proposes a semideﬁnite optimization model to address minimax two-stage stochastic linear problems with risk aversion,when the distribution of the second-stage random variables belongs to a set of multivariate distributions with knownﬁrst and second moments.The minimax solutions provide a natural distribu-tion to stress-test stochastic optimization problems under distributional ambiguity.Cromvik and Patriksson(2010a) show that,under certain assumptions,global optima and stationary solutions of stochastic mathematical programs with equilibrium constraints are robust with respect to changes in the underlying probability distribution.Works such as Zhu and Fukushima(2009)and Zymler(2010)also study distributional robustness in the context of speciﬁc applications,such as portfolio management.2.5Connection with Risk TheoryBertsimas and Brown(2009)describe how to connect uncertainty sets in robust linear optimization to coherent risk measures,an example of which is Conditional Value-at-Risk.In particular,the authors show the link between polyhedral uncertainty sets of a special structure and a subclass of coherent risk measures called distortion risk measures.Independently,Chen et al.(2007)present an approach for constructing uncertainty sets for robust opti-mization using new deviation measures that capture the asymmetry of the distributions.These deviation measures lead to improved approximations of chance constraints.Dentcheva and Ruszczynski(2010)proposes the concept of robust stochastic dominance and shows its applica-tion to risk-averse optimization.They consider stochastic optimization problems where risk-aversion is expressed by a robust stochastic dominance constraint and develop necessary and sufﬁcient conditions of optimality for such optimization problems in the convex case.In the nonconvex case,they derive necessary conditions of optimality under additional smoothness assumptions of some mappings involved in the problem.2.6Nonlinear OptimizationRobust nonlinear optimization remains much less widely studied to date than its linear counterpart.Bertsimas et al.(2010c)presents a robust optimization approach for unconstrained non-convex problems and problems based on simulations.Such problems arise for instance in the partial differential equations literature and in engineering applications such as nanophotonic design.An appealing feature of the approach is that it does not assume any speciﬁc structure for the problem.The case of robust nonlinear optimization with constraints is investigated in Bertsimas et al.(2010b)with an application to radiation therapy for cancer treatment.Bertsimas and Nohadani (2010)further explore robust nonconvex optimization in contexts where solutions are not known explicitly,e.g., have to be found using simulation.They present a robust simulated annealing algorithm that improves performance and robustness of the solution.Further,Boni et al.(2008)analyzes problems with uncertain conic quadratic constraints,formulating an approx-imate robust counterpart,and Zhang(2007)provide formulations to nonlinear programming problems that are valid in the neighborhood of the nominal parameters and robust to theﬁrst order.Hsiung et al.(2008)present tractable approximations to robust geometric programming,by using piecewise-linear convex approximations of each non-linear constraint.Geometric programming is also investigated in Shen et al.(2008),where the robustness is injected at the level of the algorithm and seeks to avoid obtaining infeasible solutions because of the approximations used in the traditional approach.Interval uncertainty-based robust optimization for convex and non-convex quadratic programs are considered in Li et al.(2011).Takeda et al.(2010)studies robustness for uncertain convex quadratic programming problems with ellipsoidal uncertainties and proposes a relaxation technique based on random sampling for robust deviation optimization sserre(2011)considers minimax and robust models of polynomial optimization.A special case of nonlinear problems that are linear in the decision variables but convex in the uncertainty when the worst-case objective is to be maximized is investigated in Kawas and Thiele(2011a).In that setting,exact and tractable robust counterparts can be derived.A special class of nonconvex robust optimization is examined in Kawas and Thiele(2011b).Robust nonconvex optimization is examined in detail in Teo(2007),which presents a method that is applicable to arbitrary objective functions by iteratively moving along descent directions and terminates at a robust local minimum.3Applications of Robust OptimizationWe describe below examples to which robust optimization has been applied.While an appealing feature of robust optimization is that it leads to models that can be solved using off-the-shelf software,it is worth pointing the existence of algebraic modeling tools that facilitate the formulation and subsequent analysis of robust optimization problems on the computer(Goh and Sim,2011).3.1Production,Inventory and Logistics3.1.1Classical logistics problemsThe capacitated vehicle routing problem with demand uncertainty is studied in Sungur et al.(2008),with a more extensive treatment in Sungur(2007),and the robust traveling salesman problem with interval data in Montemanni et al.(2007).Remli and Rekik(2012)considers the problem of combinatorial auctions in transportation services when shipment volumes are uncertain and proposes a two-stage robust formulation solved using a constraint gener-ation algorithm.Zhang(2011)investigates two-stage minimax regret robust uncapacitated lot-sizing problems with demand uncertainty,in particular showing that it is polynomially solvable under the interval uncertain demand set.3.1.2SchedulingGoren and Sabuncuoglu(2008)analyzes robustness and stability measures for scheduling in a single-machine environment subject to machine breakdowns and embeds them in a tabu-search-based scheduling algorithm.Mittal (2011)investigates efﬁcient algorithms that give optimal or near-optimal solutions for problems with non-linear objective functions,with a focus on robust scheduling and service operations.Examples considered include parallel machine scheduling problems with the makespan objective,appointment scheduling and assortment optimization problems with logit choice models.Hazir et al.(2010)considers robust scheduling and robustness measures for the discrete time/cost trade-off problem.3.1.3Facility locationAn important question in logistics is not only how to operate a system most efﬁciently but also how to design it. Baron et al.(2011)applies robust optimization to the problem of locating facilities in a network facing uncertain demand over multiple periods.They consider a multi-periodﬁxed-charge network location problem for which they ﬁnd the number of facilities,their location and capacities,the production in each period,and allocation of demand to facilities.The authors show that different models of uncertainty lead to very different solution network topologies, with the model with box uncertainty set opening fewer,larger facilities.?investigate a robust version of the location transportation problem with an uncertain demand using a2-stage formulation.The resulting robust formulation is a convex(nonlinear)program,and the authors apply a cutting plane algorithm to solve the problem exactly.Atamt¨u rk and Zhang(2007)study the networkﬂow and design problem under uncertainty from a complexity standpoint,with applications to lot-sizing and location-transportation problems,while Bardossy(2011)presents a dual-based local search approach for deterministic,stochastic,and robust variants of the connected facility location problem.The robust capacity expansion problem of networkﬂows is investigated in Ordonez and Zhao(2007),which provides tractable reformulations under a broad set of assumptions.Mudchanatongsuk et al.(2008)analyze the network design problem under transportation cost and demand uncertainty.They present a tractable approximation when each commodity only has a single origin and destination,and an efﬁcient column generation for networks with path constraints.Atamt¨u rk and Zhang(2007)provides complexity results for the two-stage networkﬂow anddesign plexity results for the robust networkﬂow and network design problem are also provided in Minoux(2009)and Minoux(2010).The problem of designing an uncapacitated network in the presence of link failures and a competing mode is investigated in Laporte et al.(2010)in a railway application using a game theoretic perspective.Torres Soto(2009)also takes a comprehensive view of the facility location problem by determining not only the optimal location but also the optimal time for establishing capacitated facilities when demand and cost parameters are time varying.The models are solved using Benders’decomposition or heuristics such as local search and simulated annealing.In addition,the robust networkﬂow problem is also analyzed in Boyko(2010),which proposes a stochastic formulation of minimum costﬂow problem aimed atﬁnding network design andﬂow assignments subject to uncertain factors,such as network component disruptions/failures when the risk measure is Conditional Value at Risk.Nagurney and Qiang(2009)suggests a relative total cost index for the evaluation of transportation network robustness in the presence of degradable links and alternative travel behavior.Further,the problem of locating a competitive facility in the plane is studied in Blanquero et al.(2011)with a robustness criterion.Supply chain design problems are also studied in Pan and Nagi(2010)and Poojari et al.(2008).3.1.4Inventory managementThe topic of robust multi-stage inventory management has been investigated in detail in Bienstock and Ozbay (2008)through the computation of robust basestock levels and Ben-Tal et al.(2009)through an extension of the Afﬁnely Adjustable Robust Counterpart framework to control inventories under demand uncertainty.See and Sim (2010)studies a multi-period inventory control problem under ambiguous demand for which only mean,support and some measures of deviations are known,using a factor-based model.The parameters of the replenishment policies are obtained using a second-order conic programming problem.Song(2010)considers stochastic inventory control in robust supply chain systems.The work proposes an inte-grated approach that combines in a single step dataﬁtting and inventory optimization–using histograms directly as the inputs for the optimization model–for the single-item multi-period periodic-review stochastic lot-sizing problem.Operation and planning issues for dynamic supply chain and transportation networks in uncertain envi-ronments are considered in Chung(2010),with examples drawn from emergency logistics planning,network design and congestion pricing problems.3.1.5Industry-speciﬁc applicationsAng et al.(2012)proposes a robust storage assignment approach in unit-load warehouses facing variable supply and uncertain demand in a multi-period setting.The authors assume a factor-based demand model and minimize the worst-case expected total travel in the warehouse with distributional ambiguity of demand.A related problem is considered in Werners and Wuelﬁng(2010),which optimizes internal transports at a parcel sorting center.Galli(2011)describes the models and algorithms that arise from implementing recoverable robust optimization to train platforming and rolling stock planning,where the concept of recoverable robustness has been deﬁned in。

Multicamera People Tracking with a Probabilistic Occupancy Map

Multicamera People Tracking witha Probabilistic Occupancy MapFranc¸ois Fleuret,Je´roˆme Berclaz,Richard Lengagne,and Pascal Fua,Senior Member,IEEE Abstract—Given two to four synchronized video streams taken at eye level and from different angles,we show that we can effectively combine a generative model with dynamic programming to accurately follow up to six individuals across thousands of frames in spite of significant occlusions and lighting changes.In addition,we also derive metrically accurate trajectories for each of them.Our contribution is twofold.First,we demonstrate that our generative model can effectively handle occlusions in each time frame independently,even when the only data available comes from the output of a simple background subtraction algorithm and when the number of individuals is unknown a priori.Second,we show that multiperson tracking can be reliably achieved by processing individual trajectories separately over long sequences,provided that a reasonable heuristic is used to rank these individuals and that we avoid confusing them with one another.Index Terms—Multipeople tracking,multicamera,visual surveillance,probabilistic occupancy map,dynamic programming,Hidden Markov Model.Ç1I NTRODUCTIONI N this paper,we address the problem of keeping track of people who occlude each other using a small number of synchronized videos such as those depicted in Fig.1,which were taken at head level and from very different angles. This is important because this kind of setup is very common for applications such as video surveillance in public places.To this end,we have developed a mathematical framework that allows us to combine a robust approach to estimating the probabilities of occupancy of the ground plane at individual time steps with dynamic programming to track people over time.This results in a fully automated system that can track up to six people in a room for several minutes by using only four cameras,without producing any false positives or false negatives in spite of severe occlusions and lighting variations. As shown in Fig.2,our system also provides location estimates that are accurate to within a few tens of centimeters, and there is no measurable performance decrease if as many as20percent of the images are lost and only a small one if 30percent are.This involves two algorithmic steps:1.We estimate the probabilities of occupancy of theground plane,given the binary images obtained fromthe input images via background subtraction[7].Atthis stage,the algorithm only takes into accountimages acquired at the same time.Its basic ingredientis a generative model that represents humans assimple rectangles that it uses to create synthetic idealimages that we would observe if people were at givenlocations.Under this model of the images,given thetrue occupancy,we approximate the probabilities ofoccupancy at every location as the marginals of aproduct law minimizing the Kullback-Leibler diver-gence from the“true”conditional posterior distribu-tion.This allows us to evaluate the probabilities ofoccupancy at every location as the fixed point of alarge system of equations.2.We then combine these probabilities with a color and amotion model and use the Viterbi algorithm toaccurately follow individuals across thousands offrames[3].To avoid the combinatorial explosion thatwould result from explicitly dealing with the jointposterior distribution of the locations of individuals ineach frame over a fine discretization,we use a greedyapproach:we process trajectories individually oversequences that are long enough so that using areasonable heuristic to choose the order in which theyare processed is sufficient to avoid confusing peoplewith each other.In contrast to most state-of-the-art algorithms that recursively update estimates from frame to frame and may therefore fail catastrophically if difficult conditions persist over several consecutive frames,our algorithm can handle such situations since it computes the global optima of scores summed over many frames.This is what gives it the robustness that Fig.2demonstrates.In short,we combine a mathematically well-founded generative model that works in each frame individually with a simple approach to global optimization.This yields excellent performance by using basic color and motion models that could be further improved.Our contribution is therefore twofold.First,we demonstrate that a generative model can effectively handle occlusions at each time frame independently,even when the input data is of very poor quality,and is therefore easy to obtain.Second,we show that multiperson tracking can be reliably achieved by processing individual trajectories separately over long sequences.. F.Fleuret,J.Berclaz,and P.Fua are with the Ecole Polytechnique Fe´de´ralede Lausanne,Station14,CH-1015Lausanne,Switzerland.E-mail:{francois.fleuret,jerome.berclaz,pascal.fua}@epfl.ch..R.Lengagne is with GE Security-VisioWave,Route de la Pierre22,1024Ecublens,Switzerland.E-mail:richard.lengagne@.Manuscript received14July2006;revised19Jan.2007;accepted28Mar.2007;published online15May2007.Recommended for acceptance by S.Sclaroff.For information on obtaining reprints of this article,please send e-mail to:tpami@,and reference IEEECS Log Number TPAMI-0521-0706.Digital Object Identifier no.10.1109/TPAMI.2007.1174.0162-8828/08/$25.00ß2008IEEE Published by the IEEE Computer SocietyIn the remainder of the paper,we first briefly review related works.We then formulate our problem as estimat-ing the most probable state of a hidden Markov process and propose a model of the visible signal based on an estimate of an occupancy map in every time frame.Finally,we present our results on several long sequences.2R ELATED W ORKState-of-the-art methods can be divided into monocular and multiview approaches that we briefly review in this section.2.1Monocular ApproachesMonocular approaches rely on the input of a single camera to perform tracking.These methods provide a simple and easy-to-deploy setup but must compensate for the lack of 3D information in a single camera view.2.1.1Blob-Based MethodsMany algorithms rely on binary blobs extracted from single video[10],[5],[11].They combine shape analysis and tracking to locate people and maintain appearance models in order to track them,even in the presence of occlusions.The Bayesian Multiple-BLob tracker(BraMBLe)system[12],for example,is a multiblob tracker that generates a blob-likelihood based on a known background model and appearance models of the tracked people.It then uses a particle filter to implement the tracking for an unknown number of people.Approaches that track in a single view prior to computing correspondences across views extend this approach to multi camera setups.However,we view them as falling into the same category because they do not simultaneously exploit the information from multiple views.In[15],the limits of the field of view of each camera are computed in every other camera from motion information.When a person becomes visible in one camera,the system automatically searches for him in other views where he should be visible.In[4],a background/foreground segmentation is performed on calibrated images,followed by human shape extraction from foreground objects and feature point selection extraction. Feature points are tracked in a single view,and the system switches to another view when the current camera no longer has a good view of the person.2.1.2Color-Based MethodsTracking performance can be significantly increased by taking color into account.As shown in[6],the mean-shift pursuit technique based on a dissimilarity measure of color distributions can accurately track deformable objects in real time and in a monocular context.In[16],the images are segmented pixelwise into different classes,thus modeling people by continuously updated Gaussian mixtures.A standard tracking process is then performed using a Bayesian framework,which helps keep track of people,even when there are occlusions.In such a case,models of persons in front keep being updated, whereas the system stops updating occluded ones,which may cause trouble if their appearances have changed noticeably when they re-emerge.More recently,multiple humans have been simulta-neously detected and tracked in crowded scenes[20]by using Monte-Carlo-based methods to estimate their number and positions.In[23],multiple people are also detected and tracked in front of complex backgrounds by using mixture particle filters guided by people models learned by boosting.In[9],multicue3D object tracking is addressed by combining particle-filter-based Bayesian tracking and detection using learned spatiotemporal shapes.This ap-proach leads to impressive results but requires shape, texture,and image depth information as input.Finally, Smith et al.[25]propose a particle-filtering scheme that relies on Markov chain Monte Carlo(MCMC)optimization to handle entrances and departures.It also introduces a finer modeling of interactions between individuals as a product of pairwise potentials.2.2Multiview ApproachesDespite the effectiveness of such methods,the use of multiple cameras soon becomes necessary when one wishes to accurately detect and track multiple people and compute their precise3D locations in a complex environment. Occlusion handling is facilitated by using two sets of stereo color cameras[14].However,in most approaches that only take a set of2D views as input,occlusion is mainly handled by imposing temporal consistency in terms of a motion model,be it Kalman filtering or more general Markov models.As a result,these approaches may not always be able to recover if the process starts diverging.2.2.1Blob-Based MethodsIn[19],Kalman filtering is applied on3D points obtained by fusing in a least squares sense the image-to-world projections of points belonging to binary blobs.Similarly,in[1],a Kalman filter is used to simultaneously track in2D and3D,and objectFig.1.Images from two indoor and two outdoor multicamera video sequences that we use for our experiments.At each time step,we draw a box around people that we detect and assign to them an ID number that follows them throughout thesequence.Fig.2.Cumulative distributions of the position estimate error on a3,800-frame sequence(see Section6.4.1for details).locations are estimated through trajectory prediction during occlusion.In[8],a best hypothesis and a multiple-hypotheses approaches are compared to find people tracks from 3D locations obtained from foreground binary blobs ex-tracted from multiple calibrated views.In[21],a recursive Bayesian estimation approach is used to deal with occlusions while tracking multiple people in multiview.The algorithm tracks objects located in the intersections of2D visual angles,which are extracted from silhouettes obtained from different fixed views.When occlusion ambiguities occur,multiple occlusion hypotheses are generated,given predicted object states and previous hypotheses,and tested using a branch-and-merge strategy. The proposed framework is implemented using a customized particle filter to represent the distribution of object states.Recently,Morariu and Camps[17]proposed a method based on dimensionality reduction to learn a correspondence between the appearance of pedestrians across several views. This approach is able to cope with the severe occlusion in one view by exploiting the appearance of the same pedestrian on another view and the consistence across views.2.2.2Color-Based MethodsMittal and Davis[18]propose a system that segments,detects, and tracks multiple people in a scene by using a wide-baseline setup of up to16synchronized cameras.Intensity informa-tion is directly used to perform single-view pixel classifica-tion and match similarly labeled regions across views to derive3D people locations.Occlusion analysis is performed in two ways:First,during pixel classification,the computa-tion of prior probabilities takes occlusion into account. Second,evidence is gathered across cameras to compute a presence likelihood map on the ground plane that accounts for the visibility of each ground plane point in each view. Ground plane locations are then tracked over time by using a Kalman filter.In[13],individuals are tracked both in image planes and top view.The2D and3D positions of each individual are computed so as to maximize a joint probability defined as the product of a color-based appearance model and2D and 3D motion models derived from a Kalman filter.2.2.3Occupancy Map MethodsRecent techniques explicitly use a discretized occupancy map into which the objects detected in the camera images are back-projected.In[2],the authors rely on a standard detection of stereo disparities,which increase counters associated to square areas on the ground.A mixture of Gaussians is fitted to the resulting score map to estimate the likely location of individuals.This estimate is combined with a Kallman filter to model the motion.In[26],the occupancy map is computed with a standard visual hull procedure.One originality of the approach is to keep for each resulting connex component an upper and lower bound on the number of objects that it can contain. Based on motion consistency,the bounds on the various components are estimated at a certain time frame based on the bounds of the components at the previous time frame that spatially intersect with it.Although our own method shares many features with these techniques,it differs in two important respects that we will highlight:First,we combine the usual color and motion models with a sophisticated approach based on a generative model to estimating the probabilities of occu-pancy,which explicitly handles complex occlusion interac-tions between detected individuals,as will be discussed in Section5.Second,we rely on dynamic programming to ensure greater stability in challenging situations by simul-taneously handling multiple frames.3P ROBLEM F ORMULATIONOur goal is to track an a priori unknown number of people from a few synchronized video streams taken at head level. In this section,we formulate this problem as one of finding the most probable state of a hidden Markov process,given the set of images acquired at each time step,which we will refer to as a temporal frame.We then briefly outline the computation of the relevant probabilities by using the notations summarized in Tables1and2,which we also use in the following two sections to discuss in more details the actual computation of those probabilities.3.1Computing the Optimal TrajectoriesWe process the video sequences by batches of T¼100frames, each of which includes C images,and we compute the most likely trajectory for each individual.To achieve consistency over successive batches,we only keep the result on the first 10frames and slide our temporal window.This is illustrated in Fig.3.We discretize the visible part of the ground plane into a finite number G of regularly spaced2D locations and we introduce a virtual hidden location H that will be used to model entrances and departures from and into the visible area.For a given batch,let L t¼ðL1t;...;L NÃtÞbe the hidden stochastic processes standing for the locations of individuals, whether visible or not.The number NÃstands for the maximum allowable number of individuals in our world.It is large enough so that conditioning on the number of visible ones does not change the probability of a new individual entering the scene.The L n t variables therefore take values in f1;...;G;Hg.Given I t¼ðI1t;...;I C tÞ,the images acquired at time t for 1t T,our task is to find the values of L1;...;L T that maximizePðL1;...;L T j I1;...;I TÞ:ð1ÞAs will be discussed in Section 4.1,we compute this maximum a posteriori in a greedy way,processing one individual at a time,including the hidden ones who can move into the visible scene or not.For each one,the algorithm performs the computation,under the constraint that no individual can be at a visible location occupied by an individual already processed.In theory,this approach could lead to undesirable local minima,for example,by connecting the trajectories of two separate people.However,this does not happen often because our batches are sufficiently long.To further reduce the chances of this,we process individual trajectories in an order that depends on a reliability score so that the most reliable ones are computed first,thereby reducing the potential for confusion when processing the remaining ones. This order also ensures that if an individual remains in the hidden location,then all the other people present in the hidden location will also stay there and,therefore,do not need to be processed.FLEURET ET AL.:MULTICAMERA PEOPLE TRACKING WITH A PROBABILISTIC OCCUPANCY MAP269Our experimental results show that our method does not suffer from the usual weaknesses of greedy algorithms such as a tendency to get caught in bad local minima.We thereforebelieve that it compares very favorably to stochastic optimization techniques in general and more specifically particle filtering,which usually requires careful tuning of metaparameters.3.2Stochastic ModelingWe will show in Section 4.2that since we process individual trajectories,the whole approach only requires us to define avalid motion model P ðL n t þ1j L nt ¼k Þand a sound appearance model P ðI t j L n t ¼k Þ.The motion model P ðL n t þ1j L nt ¼k Þ,which will be intro-duced in Section 4.3,is a distribution into a disc of limited radiusandcenter k ,whichcorresponds toalooseboundonthe maximum speed of a walking human.Entrance into the scene and departure from it are naturally modeled,thanks to the270IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,VOL.30,NO.2,FEBRUARY 2008TABLE 2Notations (RandomQuantities)Fig.3.Video sequences are processed by batch of 100frames.Only the first 10percent of the optimization result is kept and the rest is discarded.The temporal window is then slid forward and the optimiza-tion is repeated on the new window.TABLE 1Notations (DeterministicQuantities)hiddenlocation H,forwhichweextendthemotionmodel.The probabilities to enter and to leave are similar to the transition probabilities between different ground plane locations.In Section4.4,we will show that the appearance model PðI t j L n t¼kÞcan be decomposed into two terms.The first, described in Section4.5,is a very generic color-histogram-based model for each individual.The second,described in Section5,approximates the marginal conditional probabil-ities of occupancy of the ground plane,given the results of a background subtractionalgorithm,in allviewsacquired atthe same time.This approximation is obtained by minimizing the Kullback-Leibler divergence between a product law and the true posterior.We show that this is equivalent to computing the marginal probabilities of occupancy so that under the product law,the images obtained by putting rectangles of human sizes at occupied locations are likely to be similar to the images actually produced by the background subtraction.This represents a departure from more classical ap-proaches to estimating probabilities of occupancy that rely on computing a visual hull[26].Such approaches tend to be pessimistic and do not exploit trade-offs between the presence of people at different locations.For instance,if due to noise in one camera,a person is not seen in a particular view,then he would be discarded,even if he were seen in all others.By contrast,in our probabilistic framework,sufficient evidence might be present to detect him.Similarly,the presence of someone at a specific location creates an occlusion that hides the presence behind,which is not accounted for by the hull techniques but is by our approach.Since these marginal probabilities are computed indepen-dently at each time step,they say nothing about identity or correspondence with past frames.The appearance similarity is entirely conveyed by the color histograms,which has experimentally proved sufficient for our purposes.4C OMPUTATION OF THE T RAJECTORIESIn Section4.1,we break the global optimization of several people’s trajectories into the estimation of optimal individual trajectories.In Section 4.2,we show how this can be performed using the classical Viterbi’s algorithm based on dynamic programming.This requires a motion model given in Section 4.3and an appearance model described in Section4.4,which combines a color model given in Section4.5 and a sophisticated estimation of the ground plane occu-pancy detailed in Section5.We partition the visible area into a regular grid of G locations,as shown in Figs.5c and6,and from the camera calibration,we define for each camera c a family of rectangular shapes A c1;...;A c G,which correspond to crude human silhouettes of height175cm and width50cm located at every position on the grid.4.1Multiple TrajectoriesRecall that we denote by L n¼ðL n1;...;L n TÞthe trajectory of individual n.Given a batch of T temporal frames I¼ðI1;...;I TÞ,we want to maximize the posterior conditional probability:PðL1¼l1;...;L NÃ¼l NÃj IÞ¼PðL1¼l1j IÞY NÃn¼2P L n¼l n j I;L1¼l1;...;L nÀ1¼l nÀ1ÀÁ:ð2ÞSimultaneous optimization of all the L i s would beintractable.Instead,we optimize one trajectory after theother,which amounts to looking for^l1¼arg maxlPðL1¼l j IÞ;ð3Þ^l2¼arg maxlPðL2¼l j I;L1¼^l1Þ;ð4Þ...^l NÃ¼arg maxlPðL NÃ¼l j I;L1¼^l1;L2¼^l2;...Þ:ð5ÞNote that under our model,conditioning one trajectory,given other ones,simply means that it will go through noalready occupied location.In other words,PðL n¼l j I;L1¼^l1;...;L nÀ1¼^l nÀ1Þ¼PðL n¼l j I;8k<n;8t;L n t¼^l k tÞ;ð6Þwhich is PðL n¼l j IÞwith a reduced set of the admissiblegrid locations.Such a procedure is recursively correct:If all trajectoriesestimated up to step n are correct,then the conditioning onlyimproves the estimate of the optimal remaining trajectories.This would suffice if the image data were informative enoughso that locations could be unambiguously associated toindividuals.In practice,this is obviously rarely the case.Therefore,this greedy approach to optimization has un-desired side effects.For example,due to partly missinglocalization information for a given trajectory,the algorithmmight mistakenly start following another person’s trajectory.This is especially likely to happen if the tracked individualsare located close to each other.To avoid this kind of failure,we process the images bybatches of T¼100and first extend the trajectories that havebeen found with high confidence,as defined below,in theprevious batches.We then process the lower confidenceones.As a result,a trajectory that was problematic in thepast and is likely to be problematic in the current batch willbe optimized last and,thus,prevented from“stealing”somebody else’s location.Furthermore,this approachincreases the spatial constraints on such a trajectory whenwe finally get around to estimating it.We use as a confidence score the concordance of theestimated trajectories in the previous batches and thelocalization cue provided by the estimation of the probabil-istic occupancy map(POM)described in Section5.Moreprecisely,the score is the number of time frames where theestimated trajectory passes through a local maximum of theestimated probability of occupancy.When the POM does notdetect a person on a few frames,the score will naturallydecrease,indicating a deterioration of the localizationinformation.Since there is a high degree of overlappingbetween successive batches,the challenging segment of atrajectory,which is due to the failure of the backgroundsubtraction or change in illumination,for instance,is met inseveral batches before it actually happens during the10keptframes.Thus,the heuristic would have ranked the corre-sponding individual in the last ones to be processed whensuch problem occurs.FLEURET ET AL.:MULTICAMERA PEOPLE TRACKING WITH A PROBABILISTIC OCCUPANCY MAP2714.2Single TrajectoryLet us now consider only the trajectory L n ¼ðL n 1;...;L nT Þof individual n over T temporal frames.We are looking for thevalues ðl n 1;...;l nT Þin the subset of free locations of f 1;...;G;Hg .The initial location l n 1is either a known visible location if the individual is visible in the first frame of the batch or H if he is not.We therefore seek to maximizeP ðL n 1¼l n 1;...;L n T ¼l nt j I 1;...;I T Þ¼P ðI 1;L n 1¼l n 1;...;I T ;L n T ¼l nT ÞP ðI 1;...;I T Þ:ð7ÞSince the denominator is constant with respect to l n ,we simply maximize the numerator,that is,the probability of both the trajectories and the images.Let us introduce the maximum of the probability of both the observations and the trajectory ending up at location k at time t :Èt ðk Þ¼max l n 1;...;l nt À1P ðI 1;L n 1¼l n 1;...;I t ;L nt ¼k Þ:ð8ÞWe model jointly the processes L n t and I t with a hidden Markov model,that isP ðL n t þ1j L n t ;L n t À1;...Þ¼P ðL n t þ1j L nt Þð9ÞandP ðI t ;I t À1;...j L n t ;L nt À1;...Þ¼YtP ðI t j L n t Þ:ð10ÞUnder such a model,we have the classical recursive expressionÈt ðk Þ¼P ðI t j L n t ¼k Þ|ﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄ}Appearance modelmax P ðL n t ¼k j L nt À1¼ Þ|ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ}Motion modelÈt À1ð Þð11Þto perform a global search with dynamic programming,which yields the classic Viterbi algorithm.This is straight-forward,since the L n t s are in a finite set of cardinality G þ1.4.3Motion ModelWe chose a very simple and unconstrained motion model:P ðL n t ¼k j L nt À1¼ Þ¼1=Z Áe À k k À k if k k À k c 0otherwise ;&ð12Þwhere the constant tunes the average human walkingspeed,and c limits the maximum allowable speed.This probability is isotropic,decreases with the distance from location k ,and is zero for k k À k greater than a constantmaximum distance.We use a very loose maximum distance cof one square of the grid per frame,which corresponds to a speed of almost 12mph.We also define explicitly the probabilities of transitions to the parts of the scene that are connected to the hidden location H .This is a single door in the indoor sequences and all the contours of the visible area in the outdoor sequences in Fig.1.Thus,entrance and departure of individuals are taken care of naturally by the estimation of the maximum a posteriori trajectories.If there are enough evidence from the images that somebody enters or leaves the room,then this procedure will estimate that the optimal trajectory does so,and a person will be added to or removed from the visible area.4.4Appearance ModelFrom the input images I t ,we use background subtraction to produce binary masks B t such as those in Fig.4.We denote as T t the colors of the pixels inside the blobs and treat the rest of the images as background,which is ignored.Let X tk be a Boolean random variable standing for the presence of an individual at location k of the grid at time t .In Appendix B,we show thatP ðI t j L n t ¼k Þzﬄﬄﬄﬄﬄﬄﬄﬄﬄ}|ﬄﬄﬄﬄﬄﬄﬄﬄﬄ{Appearance model/P ðL n t ¼k j X kt ¼1;T t Þ|ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ}Color modelP ðX kt ¼1j B t Þ|ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ}Ground plane occupancy:ð13ÞThe ground plane occupancy term will be discussed in Section 5,and the color model term is computed as follows.4.5Color ModelWe assume that if someone is present at a certain location k ,then his presence influences the color of the pixels located at the intersection of the moving blobs and the rectangle A c k corresponding to the location k .We model that dependency as if the pixels were independent and identically distributed and followed a density in the red,green,and blue (RGB)space associated to the individual.This is far simpler than the color models used in either [18]or [13],which split the body area in several subparts with dedicated color distributions,but has proved sufficient in practice.If an individual n was present in the frames preceding the current batch,then we have an estimation for any camera c of his color distribution c n ,since we have previously collected the pixels in all frames at the locations272IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,VOL.30,NO.2,FEBRUARY2008Fig.4.The color model relies on a stochastic modeling of the color of the pixels T c t ðk Þsampled in the intersection of the binary image B c t produced bythe background subtraction and the rectangle A ck corresponding to the location k .。

德国工业4.0原版

z
Intense research activities in universities and other research institutions Drastically increasing number of publications in recent years Large amount of funding by the German government
Model predictive control (MPC)
Modern, optimization-based control technique Successful applications in many industrial ﬁelds Can handle hard constraints on states and inputs Optimization of some performance criterion Applicable to nonlinear, MIMO systems
A system is strictly dissipative on a set W ⊆ Z with respect to the supply rate s if there exists a storage function λ such that for all (x , u ) ∈ W it holds that λ(f (x , u )) − λ(x ) ≤ s (x , u ) − ρ(x ) with ρ > 0.
k =0 x (k |t + 1) x (t + 1) state x input u t+1 u (k |t + 1) k =N
Basic MPC scheme

Dirichlet boundary control of semilinear parabolic equation

126
N. Arada and J.-P. Raymond
initial condition belongs to Wad ⊂ L ∞ ( ), and A is a second-order elliptic operator of the form Ay(x ) = − iN , j =1 Di (aij ( x ) D j y ( x )) (where Di denotes the partial derivate with respect to xi ). We are interested in optimality conditions in the form of Pontryagin’s principles for the control problem (P) where J ( y , u , v, w) =
Linear-quadratic problems of the form ( P ) (when (x , t , y , u ) is of the form a (x , t ) y + u and J is quadratic) have been studied for 30 years [18], [14]. More recently the Dynamic Programming approach has been developed for Dirichlet boundary controls of nonlinear parabolic equations [5], [4], [12]. In this approach the state variable is a time-dependent function with values in a Hilbert space. Very recently, growing interest has been taken in the corresponding problems with pointwise state constraints [6], [9], [10], [17], [20], [21]. For such problems, it is natural to look for solutions in L ∞ ( Q ) or in Cb ( Q ) (the space of bounded continuous functions in Q ). The main purpose of Part 1 is to develop new tools to analyze the state equation when solutions are bounded and continuous in Q , and next to obtain Pontryagin’s principles. The corresponding problem with pointwise state constraints is studied in Part 2. Taking advantage of the fact that v belong to L ∞ ( ), we prove existence and regularity results for (1) by a method different from the ones previously mentioned. We deﬁne the solution y of (1) by the transposition method and we prove that it is the limit of a sequence of solutions for equations with a Robin boundary condition (Theorem 3.9). This kind of property is well known for solutions of linear hyperbolic equations [15], [19], or linear elliptic equations [7]. Here the novelty is that we use the same method for nonlinear equations to prove that solutions belong to Cb ( Q ∪ T ). Moreover, following [16], we think that this method may also be interesting for numerical approximations. We prove that optimal solutions satisfy three decoupled Pontryagin’s principles, one for the distributed control, one for the boundary control, and the last one for the control in the initial condition (Theorem 2.1). These optimality conditions are obtained with a method of diffuse perturbations. Expansions for diffuse perturbations of a boundary | d (x , ) > τ }× ]τ, T [ (see control are obtained only in subcylinders Q τ = {x ∈ Theorem 5.2). Such tools will be very useful to study the corresponding problems with pointwise state constraints in Part 2 [2]. The plan of the paper is as follows. Assumptions and the main result (Theorem 2.1) are stated in Section 2. The state equation, the linearized state equation, and estimates needed for the proof of Theorem 2.1 are studied in Section 3. We recall some results for the adjoint equation in Section 4. Taylor expansions for diffuse perturbations are obtained in Section 5 and the proof of Theorem 2.1 is given in Secttions and Main Result

药物经济学嵌入决策树的马尔可夫模型

药物经济学嵌入决策树的马尔可夫模型嵌入决策树的马尔可夫模型是一种常用的经济模型，它结合了决策树和马尔可夫链两个工具，可以更全面地评估药物治疗方案的效果和成本。

决策树被用来描述治疗决策的不同选项和可能的结果，而马尔可夫链则用来描述状态的转移和概率。

在嵌入决策树的马尔可夫模型中，决策树的根节点代表治疗方案的选择，每个分支节点代表其中一种治疗方案的不同决策路径，叶节点代表治疗结束后的结果。

根据治疗路径的选择，马尔可夫链描述了不同状态之间的转移概率。

每个状态代表疾病的不同阶段或治疗效果的不同等级。

在模型中，除了决策路径和转移概率，还需要考虑治疗方案的成本和效果。

成本包括直接医疗费用和间接费用，效果通常通过一些评估指标来衡量，例如生命质量调整年（quality-adjusted life years, QALYs）或无生命质量调整年（life years gained, LYs）等。

通过概率和成本效果的计算，可以评估每个治疗方案的成本效益比。

在实际应用中，嵌入决策树的马尔可夫模型可以用于评估不同药物治疗方案的成本效益，并帮助决策者做出更明智的决策。

例如，对于其中一种疾病的治疗方案，模型可以比较不同治疗路径的成本效益，帮助决策者选择最经济和最有效的方案。

此外，模型还可以用于评估长期治疗策略的成本效益。

由于疾病状态和治疗效果可能发生改变，模型可以预测不同治疗方案在长期使用下的效果和成本。

这有助于决策者考虑治疗决策的长期影响，并制定适当的管理策略。

总之，嵌入决策树的马尔可夫模型是一种有力的工具，可用于评估药物治疗方案的成本效益。

它综合了决策树和马尔可夫链的优势，可帮助决策者制定更明智和经济的医疗决策。

随着药物经济学的不断发展，这种模型将在药物研发和临床实践中发挥越来越重要的作用。

Empirical processes of dependent random variables

2
Preliminaries
n i=1
from R to R. The centered G -indexed empirical process is given by (P n − P )g = 1 n
n
the marginal and empirical distribution functions. Let G be a class of measurabrocesses that have been discussed include linear processes and Gaussian processes; see Dehling and Taqqu (1989) and Cs¨ org˝ o and Mielniczuk (1996) for long and short-range dependent subordinated Gaussian processes and Ho and Hsing (1996) and Wu (2003a) for long-range dependent linear processes. A collection of recent results is presented in Dehling, Mikosch and Sorensen (2002). In that collection Dedecker and Louhichi (2002) made an important generalization of Ossiander’s (1987) result. Here we investigate the empirical central limit problem for dependent random variables from another angle that avoids strong mixing conditions. In particular, we apply a martingale method and establish a weak convergence theory for stationary, causal processes. Our results are comparable with the theory for independent random variables in that the imposed moment conditions are optimal or almost optimal. We show that, if the process is short-range dependent in a certain sense, then the limiting behavior is similar to that of iid random variables in that the limiting distribution is a Gaussian process and the norming √ sequence is n. For long-range dependent linear processes, one needs to apply asymptotic √ expansions to obtain n-norming limit theorems (Section 6.2.2). The paper is structured as follows. In Section 2 we introduce some mathematical preliminaries necessary for the weak convergence theory and illustrate the essence of our approach. Two types of empirical central limit theorems are established. Empirical processes indexed by indicators of left half lines, absolutely continuous functions, and piecewise diﬀerentiable functions are discussed in Sections 3, 4 and 5 respectively. Applications to linear processes and iterated random functions are made in Section 6. Section 7 presents some integral and maximal inequalities that may be of independent interest. Some proofs are given in Sections 8 and 9.

A New Approach for Filtering Nonlinear Systems

computational overhead as the number of calculations demanded for the generation of the Jacobian and the predictions of state estimate and covariance are large. In this paper we describe a new approach to generalising the Kalman ﬁlter to systems with nonlinear state transition and observation models. In Section 2 we describe the basic ﬁltering problem and the notation used in this paper. In Section 3 we describe the new ﬁlter. The fourth section presents a summary of the theoretical analysis of the performance of the new ﬁlter against that of the EKF. In Section 5 we demonstrate the new ﬁlter in a highly nonlinear application and we conclude with a discussion of the implications of this new ﬁlter1
Tቤተ መጻሕፍቲ ባይዱ
= = =
δij Q(i), δij R(i), 0, ∀i, j.
(3) (4) (5)

序贯最小化方法

这里 ui = ∑ α
j >0
(2)
y jα j k ( x j , xi ) b,
(3)
b由式(1)计算
罗林开
模式识别与智能系统研究所
2010-102010-10-21
7
大规模凸二次规划问题的求解
但是对于大规模问题（即决策变量个数或样本维数较大），传统的方法（如有效集方法)基本上难以有效求解. 大规模凸二次规划问题的求解思路：通过求解一系列小规模（即原决策变量的一个合适子集）的凸二次规划问题，获得原问题的解，此类方法统称为分块或分解的方法.
与第一个乘子的结合，应该使第二个乘子的迭代步长较大。根据(4)式，第二个乘子的迭代步长大致正比于|E1-E2|，因此应选择因此应选择第二个乘子最大化|E1-E2|，，即当E1为正时第二个乘子最大化，，即当，，即当最小化E 当为负E 时最大化E 最小化 2，当为负 1时最大化 2.
罗林开
Osuna方法（Osuna，1997）
– 固定Chunking工作集的大小，比如每次迭代时加入1个违反KKT条件的乘子，则需在原工作集中去除一个乘子。
SMO方法（Platt，1998）
– 固定Chunking工作集的大小为2(最小).
罗林开
模式识别与智能系统研究所
2010-102010-10-21
11
两个Langrange乘子变量问题的求解
1 1 2 2 min M (α1 , α 2 ) = k11α1 + k22α 2 + y1 y2 k12α1α 2 α1 ,α 2 2 2 (α1 + α 2 ) + y1α1v1 + y2α 2 v2 + 常数 S .T . y1α1 + y2α 2 = ∑ yiα i = 常数

DB33∕T 1136-2017 建筑地基基础设计规范

5
地基计算 ....................................................................................................................... 14 5.1 承载力计算......................................................................................................... 14 5.2 变形计算 ............................................................................................................ 17 5.3 稳定性计算......................................................................................................... 21
主要起草人：施祖元刘兴旺潘秋元陈云敏王立忠李冰河（以下按姓氏拼音排列）蔡袁强陈青佳陈仁朋陈威文陈舟樊良本胡凌华胡敏云蒋建良李建宏王华俊刘世明楼元仓陆伟国倪士坎单玉川申屠团兵陶琨叶军徐和财许国平杨桦杨学林袁静主要审查人：益德清龚晓南顾国荣钱力航黄茂松朱炳寅朱兆晴赵竹占姜天鹤赵宇宏童建国浙江大学参编单位：（排名不分先后）浙江工业大学温州大学华东勘测设计研究院有限公司浙江大学建筑设计研究院有限公司杭州市建筑设计研究院有限公司浙江省建筑科学设计研究院汉嘉设计集团股份有限公司杭州市勘测设计研究院宁波市建筑设计研究院有限公司温州市建筑设计研究院温州市勘察测绘院中国联合工程公司浙江省电力设计院浙江省省直建筑设计院浙江省水利水电勘测设计院浙江省工程勘察院大象建筑设计有限公司浙江东南建筑设计有限公司湖州市城市规划设计研究院浙江省工业设计研究院浙江工业大学工程设计集团有限公司中国美术学院风景建筑设计研究院华汇工程设计集团股份有限公司

晶体结构测定程序简介

晶体结构测定程序简介一. 程序包SHELXTL PC 应包含如下主要程序:1. 数据准备程序 XPREP.EXE另外还需要有衍射数据:*.HKL 文件, 晶胞参数和分子式.2. 起始套产生程序 XS.EXE3. 原子挑选,命名 XP.EXE4. 数据还原及最小二乘修正程序 XLS.EXE 等二. 数据准备程序 XPREP.EXE该程序用于对原始衍射数据进行预处理: 当输入晶包参数后,程序会根据数据的消光规律测定空间群; 输入分子式后会产生元素表等,最后产生 *.PRP *.INS 文件,以待往下解结构用. 具体操作如下:在DOS状态下执行C:>XPREP NAME 回车, 会出现如下一系列提问, 请按例子回答:+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XPREP - RECIPROCAL SPACE EXPLORATION - Version 4.20.1 for MSDOS ++ Copyright 1990 Siemens Analytical Xray Inst. Inc., All Rights Reserved ++ Licensed to: UNIV OF CAL SAN DIEGO ++ CUPY started at 22:27:41 on 9 Aug 1995 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++Original cell in Angstroms and degrees:*注: 输入晶胞参数 a b c αβγ:15.623 23.756 11.036 90.00 90.00 90.004069 reflections read from file CUPY.HKL Mean(I/sigma) = 7.81. Lattice exceptions: P A B C I F ObvRev AllN (total) = 0 2026 2028 2030 2040 3042 2710 2710 4069N (int>3sigma) = 0 1131 961 1022 1026 1557 1480 1479 2203Mean intensity = 0.0 50.6 42.7 31.1 53.3 41.4 68.4 72.7 69.0Mean int/sigma = 0.0 6.9 6.0 5.2 6.9 6.0 7.9 7.9 7.8Lattice type P chosen Volume: 4095.90-------------------------------------------------------------------------------Determination of reduced (Niggli) cellTransformation from original cell (HKLF-matrix):0.0000 0.0000 1.0000 1.0000 0.0000 0.0000 0.0000 1.00000.0000Unitcell: 11.036 15.623 23.756 90.00 90.00 90.00Niggli form: a.a = 121.79 b.b = 244.08 c.c = 564.35b.c = 0.00 a.c = 0.00 a.b = 0.00-------------------------------------------------------------------------------Search for higher METRIC symmetry-------------------------------------------------------------------------------Option A: FOM = 0.000 deg. ORTHORHOMBIC P-lattice R(int) = 0.000 [ 0]Cell: 11.036 15.623 23.756 90.00 90.00 90.00 Volume: 4095.90Matrix: 0.0000 0.0000 1.0000 1.0000 0.0000 0.0000 0.0000 1.0000 0.0000Option A selected-------------------------------------------------------------------------------Space group determinationLattice exceptions: P A B C I F ObvRev AllN (total) = 0 2030 2026 2028 2040 3042 2713 2710 4069N (int>3sigma) = 0 1022 1131 961 1026 1557 1477 1480 2203Mean intensity = 0.0 31.1 50.6 42.7 53.3 41.4 70.0 68.4 69.0Mean int/sigma = 0.0 5.2 6.9 6.0 6.9 6.0 7.7 7.9 7.8Crystal system O and Lattice type P selectedMean |E*E-1| = 1.003 [expected .968 centrosym and .736 non-centrosym]*注: |E*E-1|值接近0.968为中心对称, 接近0.736为非中心对称.Chiral flag NOT setSystematic absence exceptions:b-- c-- n-- 21-- -c- -a- -n- -21- --a --b --n--21N 210 211 217 7 152 149 155 9 100 100 104 13N I>3s 6 58 56 1 51 52 3 1 74 73 69 0<I> 1.5 13.2 12.7 1.4 28.1 28.5 1.5 1.1 128.1 77.5 99.21.1<I/s> 0.6 2.9 2.8 1.0 4.7 4.9 0.5 0.8 14.7 10.8 12.3 0.3Option Space Group No. Type Axes CSD R(int) N(eq) Syst. Abs. CFOM[A] Pna2(1) # 33 non-cen 2 903 0.000 0 1.0 / 2.811.69[B] Pnma # 62 centro 3 894 0.000 0 1.0 / 2.84.68Option [B] chosen2.-------------------------------------------------------------------------------*注: 自动选择CFOM值小的空间群.No. 号为该空间群在国际表中的序号.CSD 为该空间群在实际晶体中出现的几率.Determination of unit-cell contentsFormula: [Cu2(C5H5NCH2CO2)4(H2O)Cl](ClO4)3.2H2O*注: 分子式的原子是双字母的,第二字母要小写.Formula weight = 1065.35Tentative Z (number of formula units/cell) = 4.0 giving rho = 1.728,non-H atomic volume = 16.8 and following cell contents and analysis:C 112.00 31.57 % H 128.00 3.03 % N 16.00 5.26 % O 91.20 34.24 % Cl 16.80 13.98 % Cu 8.00 11.93 % F(000) = 2159.2 Mo-K(alpha) radiation Mu (mm-1) = 1.40-------------------------------------------------------------------------------File CUPY.INS set up as follows:TITL CUPY in PnmaCELL 0.71073 15.623 23.756 11.036 90.00 90.00 90.00ZERR 4.00 0.003 0.005 0.002 0.00 0.00 0.00LATT 1SYMM .5-X, -Y, .5+ZSYMM -X, .5+Y, -ZSYMM .5+X, .5-Y, .5-ZSFAC C H N O CL CUUNIT 112 128 16 91.2 16.8 8TREFHKLF 4END-------------------------------------------------------------------------------*注:1. TREF 为直接法, 可在其后加 200～1000 的数, 表示直接法计算200～1000次.对有重金属的分子, 也可换用帕特逊法PATT.2. HKL 4 后若有转向矩阵数据要抄下, 第一次XLS修正完后要加回INS文件中.3. 若有提示:" *** polar space group -- origin needs to be fixed in Y *** ".则做完XS后, 把一个原子(通常是重原子)的Y座标加10作固定.三. 原子挑选,命名 XP.EXE该程序主要用于对XS.EXE 或XLS.EXE 产生的分子结构粗结果数据 NAME.RES 进行挑选,对其中产生的新原子命名, 删除不合理原子; 看分子立体图形, 输出图形文件, 看键长键角, 及面间夹角等, 具体操作及常用命令如下:在DOS状态下执行 C:>XP NAME 回车, 表示已调入NAME.RES 数据,3.(若在执行 C:>XP NAME.INS 回车, 则表示已调入NAME.INS 数据,) 在提示符:" # " 状态下可执行如下常用命令: (也可键入 HELP XP;HELPFACILITY 等看求助)# REAP NAME2.INS 把NAME2.INS数据文件读入;# READ NAME2 把NAME2.RES数据文件非Q值读入;# FMOL 把所有原子座标(包括未命名的原子Q??)读入内存, 以便处理.# FMOL LESS $Q 除了未命名的各原子Q??外, 把所有原子座标读入内存读入内存. # MPLN/N 找最小二乘平面;# PROJ 看分子的立体图. 其中右边菜单命令:看所有原子的名称;看立体图, 不打名称;看向左旋转的立体图;看向右旋转的立体图;看从下向上旋转的立体图;看从上向下旋转的立体图;看平面右旋的立体图;看平面左旋的立体图;退出, 回 # 状态. 也可按ESC键退出.# PICK 逐个挑选、命名原子: 当某个原子名在闪时, 键入新名,按回车. 若是: 按空格键: 留下原名不变; 按ENTER键: 删除;按←BACKSPACE键: 退回前一个原子; 按ESC键: 以前所做的命令、挑选全作废;按 /键: 提前结束挑选工作, 回 # 状态;# INFO C?? 显示C1～C99 的座标参数;# PERS 看球棒状分子图.# ENVI C1 看C1原子与周围原子相连的键长、键角.# ENVI CU 2 看CU原子半径+2 以内之周边原子;# UNIQ C1 以后操作只处理与C1相连的那一片分子. ( 按FMOL 后又可处理所有分子.)# NAME Q1 C1 把Q1命名为C1.# NAME Q?? C?? 把所有Q原子(Q1至Q99)对应地命名为各C原子(C1至C99).# HADD 给整个分子加氢. ( 可能有个别不合理.)# HADD 1 给季胺加 H;# HADD 2 给次甲基加 H;# HADD 3 给甲基加 H;# HADD 4 C1 C2 C3 C4 C5 C6 给苯环(C1 C2 C3 C4 C5 C6)加 H;# HADD 5 给季胺加 H;# HADD 8 给 O 加 H : O-H;# HADD 9 给双键端基 =C ,=N 等加 H;# KILL C1 C2 …删除C1 C2 …等各原子.# KILL Q10 TO Q20 只删除Q10至Q20的原子.(要看FMOL后的表来删.)# KILL $Q 删所有的Q原子.# KILL Q1? 删Q1和Q10至Q19的原子.4.看环1(C1 C2 C3 C4 C5 C6)与环2(C7 C8 C9 C10 C11 C12)之间的夹角:# MPLN C1 C2 C3 C4 C5 C6 回车,# MPLN C7 C8 C9 C10 C11 C12 回车则显示环1与环2的夹角,及环2(C7 C8 C9 C10 C11 C12)的共面(好坏)因子; # SORT C1 C2 C3 …, 原子排序;# SYEN 5655 Ci,C2 把Ci C2 原子按对数号5655找出;# VIEW 看立体彩图;# UNDO C1 C2 不连C1与C2;# LINK C1 C2 连C1与C2;# GROW CU1 只操作CU1所连的一群;若O1原子边(0.5)以内有另一原子C4A, 可把O1与C4A取平均:# ENVI O1 找O1的对数卡(如4656);# SGEN 4656 O1; 把O1对称性相联的另一O1找出; 得:# O1A# CENT/X O1 O1A 取O1 O1A 平均值,得:# X1A …# KILL O1 O1A# NAME X1A O1 。

2005-A global Malmquist productivity index

A global Malmquist productivity indexJesu ´s T.Pastor a ,C.A.Knox Lovell b ,TaCentro de Investigacio ´n Operativa,Universidad Miguel Herna ´ndez,03206Elche (Alicante),SpainbDepartment of Economics,University of Georgia,Athens,GA 30602,USA Received 2June 2004;received in revised form 24January 2005;accepted 16February 2005Available online 23May 2005AbstractThe geometric mean Malmquist productivity index is not circular,and its adjacent period components can provide different measures of productivity change.We propose a global Malmquist productivity index that is circular,and that gives a single measure of productivity change.D 2005Elsevier B.V .All rights reserved.Keywords:Malmquist productivity index;Circularity JEL classification:C43;D24;O471.IntroductionThe geometric mean form of the contemporaneous Malmquist productivity index,introduced by Caves et al.(1982),is not circular.Whether this is a serious problem depends on the powers of persuasion of Fisher (1922),who dismissed the test,and Frisch (1936),who endorsed it.The index averages two possibly disparate measures of productivity change.Fa ¨re and Grosskopf (1996)state sufficient conditions on the adjacent period technologies for the index to satisfy circularity,and to average the same measures of productivity change.When linear programming techniques are used to compute and decompose the index,infeasibility can occur.Whether this is a serious problem depends on0165-1765/$-see front matter D 2005Elsevier B.V .All rights reserved.doi:10.1016/j.econlet.2005.02.013T Corresponding author.Tel.:+17065423689;fax:+17065423376.E-mail address:knox@ (C.A.K.Lovell).Economics Letters 88(2005)266–271/locate/econbasethe structure of the data.Xue and Harker(2002)provide necessary and sufficient conditions on the datafor LP infeasibility not to occur.We demonstrate that the source of all three problems is the specification of adjacent periodtechnologies in the construction of the index.We show that it is possible to specify a base periodtechnology in a way that solves all three problems,without having to impose restrictive conditions oneither the technologies or the data.Berg et al.(1992)proposed an index that compares adjacent period data using technology from a baseperiod.This index satisfies circularity and generates a single measure of productivity change,but it paysfor circularity with base period dependence,and it remains susceptible to LP infeasibility.Shestalova(2003)proposed an index having as its base a sequential technology formed from data ofall producers in all periods up to and including the two periods being compared.This index is immune toLP infeasibility,and it generates a single measure of productivity change,but it fails circularity and itprecludes technical regress.Thus no currently available Malmquist productivity index solves all three problems.We propose anew global index with technology formed from data of all producers in all periods.This index satisfiescircularity,it generates a single measure of productivity change,it allows technical regress,and it isimmune to LP infeasibility.In Section2we introduce and decompose the circular global index.Its efficiency change componentis the same as that of the contemporaneous index,but its technical change component is new.In Section3we relate it to the contemporaneous index.In Section4we provide an empirical illustration.Section5concludes.2.The global Malmquist productivity indexConsider a panel of i=1,...,I producers and t=1,...,T time periods.Producers use inputs x a R N+toproduce outputs y a R P+.We define two technologies.A contemporaneous benchmark technology isdefined as T c t={(x t,y t)|x t can produce y t}with k T c t=T c t,t=1,...,T,k N0.A global benchmarktechnology is defined as T c G=conv{T c1v...v T c T}.The subscript b c Q indicates that both benchmark technologies satisfy constant returns to scale.A contemporaneous Malmquist productivity index is defined on T c s asM scx t;y t;x tþ1;y tþ1ÀÁ¼D scx tþ1;y tþ1ðÞD scx t;y tðÞ;ð1Þwhere the output distance functions D c s(x,y)=min{/N0|(x,y//)a T c s},s=t,t+1.Since M c t(x t,y t,x t+1, y t+1)p M c t+1(x t,y t,x t+1,y t+1)without restrictions on the two technologies,the contemporaneous index is typically defined in geometric mean form as M c(x t,y t,x t+1,y t+1)=[M c t(x t,y t,x t+1,y t+1)ÂM c t+1(x t,y t,x t+1, y t+1)]1/2.A global Malmquist productivity index is defined on T c G asM Gcx t;y t;x tþ1;y tþ1ÀÁ¼D Gcx tþ1;y tþ1ðÞD Gcx t;y tðÞ;ð2Þwhere the output distance functions D c G(x,y)=min{/N0|(x,y//)a T c G}.J.T.Pastor,C.A.K.Lovell/Economics Letters88(2005)266–271267Both indexes compare (x t +1,y t +1)to (x t ,y t ),but they use different benchmarks.Since there is only one global benchmark technology,there is no need to resort to the geometric mean convention when defining the global index.M cGdecomposes as M G c x t ;y t ;x t þ1;f y t þ1ÀÁ¼D t þ1c x t þ1;y t þ1ðÞD t c x t ;y t ðÞÂD G c x t þ1;y t þ1ðÞD t þ1c x t þ1;y t þ1ðÞÂD t cx t ;y t ðÞD Gc x t ;y t ðÞ&'¼TE t þ1c x t þ1;y t þ1ðÞTE t c x t ;y t ðÞÂD G c Àx t þ1;y t þ1=D t þ1c x t þ1;y t þ1ðÞÁD G c x t ;y t =D t cx t ;y t ðÞÀÁ()¼EC c ÂBPG G ;t þ1cx t þ1;y t þ1ðÞBPG cx t ;y tðÞ()¼EC c ÂBPC c ;ð3Þwhere EC c is the usual efficiency change indicator and BPG c G,s V 1is a best practice gap between T c Gand T c s measured along rays (x s ,y s),s =t ,t +1.BPC c is the change in BPG c ,and provides a new measure of technical change.BPC c f 1indicates whether the benchmark technology in period t +1in the region[(x t +1,y t +1/D ct +1(x t +1,y t +1))]is closer to or farther away from the global benchmark technology than is the benchmark technology in period t in the region [(x t ,y t /D ct (x t ,y t ))].M c G has four virtues.First,like any fixed base index,M cGis circular,and since EC c is circular,so is BPC c .Second,each provides a single measure,with no need to take the geometric mean of disparate adjacent period measures.Third,but not shown here,the decomposition in (3)can be extended to generate a three-way decomposition that is structurally identical to the Ray and Desli (1997)decomposition of the contemporaneous index.M cGand M c share a common efficiency change component,but they have different technical change and scale components,and so M c Gp M c without restrictions on the technologies.Finally,the technical change and scale components of M c Gare immune to the LP infeasibility problem that plagues these components of M c .paring the global and contemporaneous indexes The ratioM G c =M c¼M G c =M t þ1cÀÁÂM G c =M t cÀÁÂÃ1=2¼D G cx t þ1;y t þ1=D t þ1c x t þ1;y t þ1ðÞÀÁD G c x t ;y t =D t þ1c x t ;y t ðÞÀÁ"#ÂD G c x t þ1;y t þ1=D t c x t þ1;y t þ1ðÞÀÁD G c x t ;y t =D t c x t ;y t ðÞÀÁ"#()1=2¼BPG G ;t þ1cx t þ1;y t þ1ðÞBPG G ;t þ1cx t ;y tðÞ"#ÂBPG G ;t c xt þ1;y t þ1ðÞBPG G ;t c x t ;y tðÞ"#()1=2ð4Þis the geometric mean of two terms,each being a ratio of benchmark technology gaps along differentrays.M c G /M c f 1as projections onto T c t and T c t +1of period t +1data are closer to,equidistant from,orfarther away from T c G than projections onto T c t and T ct +1of period t data are.J.T.Pastor,C.A.K.Lovell /Economics Letters 88(2005)266–271268J.T.Pastor,C.A.K.Lovell/Economics Letters88(2005)266–271269 Table1Electricity generation data,annual means1977198219871992 Output(000MW h)13,70013,86016,18017,270 Labor(#FTE)1373179719952021 Fuel(billion BTU)1288144116671824 Capital(To¨rnqvist)44,756211,622371,041396,386 M c G=M c if BPG c G,s(x t+1,y t+1)=BPG c G,s(x t,y t),s=t,t+1.From the first equality in(4),this condition is equivalent to the condition M c G=M c s,s=t,t+1.If this condition holds for all s,it is equivalent to the condition M c t=M c1for all t.Althin(2001)has shown that a sufficient condition for base period independence is that technical change be Hicks output-neutral(HON).Hence HON is also sufficient for M c G=M c.4.An empirical illustrationWe summarize an application intended to illustrate the behavior of M c G,and to compare its performance with that of M c.We analyze a panel of93US electricity generating firms in four years (1977,1982,1997,1992).The firms use labor(FTE employees),fuel(BTUs of energy)and capital(a multilateral To¨rnqvist index)to generate electricity(net generation in MW h).The data are summarized in Table1.Electricity generation increased by proportionately less than each input did.The main cause of the rapid increase in the capital input was the enactment of environmental regulations mandating the installation of pollution abatement equipment.We are unable to disaggregate the capital input into its productive and abatement components.Empirical findings are summarized in Table2.The first three rows report decomposition(3)of M c G, and the final three rows report M c and its two adjacent period components.Columns correspond to time periods.M c G shows a large productivity decline from1977to1982,followed by weak productivity growth. Cumulative productivity in1992was25%lower than in1977.M c G calculated using1992and1977data generates the same value,verifying that it is circular.The efficiency change component EC c of M c G(and M c)is also circular,and cumulates to an18% improvement.Best practice change,BPC c,is also circular,and declined by35%.Capital investment in Table2Global and contemporaneous Malmquist productivity indexes1977–19821982–19871987–1992Cumulative productivity1977–1992 M c G0.685 1.064 1.0390.7570.757EC c 1.163 1.0890.929 1.176 1.176 BPC c0.5890.977 1.1180.6440.644M c0.4310.895 1.0390.4000.592M c t0.7130.902 1.0530.678 1.333M c t+10.2600.887 1.0240.2360.263pollution abatement equipment generated cleaner air but not more electricity.Consequently catching up with deteriorating best practice was relatively easy.Turning to the contemporaneous index M c reported in the final three rows,the story is not so clear.Cumulative productivity in 1992was 60%lower than in 1977.However calculating M c using 1992and 1977data generates a smaller 40%decline,verifying that M c is not circular.Neither figure is close to the25%decline reported by M cG,verifying that technical change was not HON,but (pollution abatement)capital-using.The lack of circularity is reflected in the frequently large differences between M ct and M c t +1,which give conflicting signals when computed using 1992and 1977data,with M c tsignaling productivitygrowth and M ct +1signaling productivity decline.Although not reported in Table 2,we have calculated three-way decompositions of M cG and M c .All three components of M c G are circular,and LP infeasibility does not occur.In contrast,the technical change and scale components of M c are not circular,and infeasibility occurs for 13observations.The circular global index M cGtells a single story about productivity change,and its decomposition is intuitively appealing in light of what we know about the industry during the cking circularity,M c and its two adjacent period components tell different stories that are often contradictory.Thedifferences between M cGand M c are a consequence of the capital-using bias of technical change,which was regressive due to the mandated installation of pollution abatement equipment,augmented perhaps by the rate base padding that was prevalent during the period.5.ConclusionsThe contemporaneous Malmquist productivity index is not circular,its adjacent period components can give conflicting signals,and it is susceptible to LP infeasibility.The global Malmquist productivity index and each of its components is circular,it provides single measures of productivity change and its components,and it is immune to LP infeasibility.The global index decomposes into the same sources of productivity change as the contemporaneous index does.A sufficient condition for equality of the two indexes,and their respective components,is Hicks output neutrality of technical change.The global index must be recomputed when a new time period is incorporated.Diewert’s (1987)assertion that b ...economic history has to be rewritten ...Q when new data are incorporated is the base period dependency problem revisited.The problem can be serious when using base periods t =1and t =T ,but it is likely to be benign when using global base periods {1,...,T }and {1,...,T +1}.While new data may change the global frontier,the rewriting of history is likely to be quantitative rather than qualitative.ReferencesAlthin,R.,2001.Measurement of productivity changes:two Malmquist index approaches.Journal of Productivity Analysis 16,107–128.Berg,S.A.,Førsund,F.R.,Jansen,E.S.,1992.Malmquist indices of productivity growth during the deregulation of Norwegian banking,1980–89.Scandinavian Journal of Economics 94,211–228(Supplement).Caves,D.W.,Christensen,L.R.,Diewert,W.E.,1982.The economic theory of index numbers and the measurement of input output,and productivity.Econometrica 50,1393–1414.J.T.Pastor,C.A.K.Lovell /Economics Letters 88(2005)266–271270J.T.Pastor,C.A.K.Lovell/Economics Letters88(2005)266–271271 Diewert,W.E.,1987.Index numbers.In:Eatwell,J.,Milgate,M.,Newman,P.(Eds.),The New Palgrave:A Dictionary of Economics,vol.2.The Macmillan Press,New York.Fa¨re,R.,Grosskopf,S.,1996.Intertemporal Production Frontiers:With Dynamic DEA.Kluwer Academic Publishers,Boston. Fisher,I.,1922.The Making of Index Numbers.Houghton Mifflin,Boston.Frisch,R.,1936.Annual survey of general economic theory:the problem of index numbers.Econometrica4,1–38.Ray,S.C.,Desli,E.,1997.Productivity growth,technical progress,and efficiency change in industrialized countries:comment.American Economic Review87,1033–1039.Shestalova,V.,2003.Sequential Malmquist indices of productivity growth:an application to OECD industrial activities.Journal of Productivity Analysis19,211–226.Xue,M.,Harker,P.T.,2002.Note:ranking DMUs with infeasible super-efficiency in DEA models.Management Science48, 705–710.。

Advances in

Advances in Geosciences,4,17–22,2005 SRef-ID:1680-7359/adgeo/2005-4-17 European Geosciences Union©2005Author(s).This work is licensed under a Creative CommonsLicense.Advances in GeosciencesIncorporating level set methods in Geographical Information Systems(GIS)for land-surface process modelingD.PullarGeography Planning and Architecture,The University of Queensland,Brisbane QLD4072,Australia Received:1August2004–Revised:1November2004–Accepted:15November2004–Published:9August2005nd-surface processes include a broad class of models that operate at a landscape scale.Current modelling approaches tend to be specialised towards one type of pro-cess,yet it is the interaction of processes that is increasing seen as important to obtain a more integrated approach to land management.This paper presents a technique and a tool that may be applied generically to landscape processes. The technique tracks moving interfaces across landscapes for processes such as waterﬂow,biochemical diffusion,and plant dispersal.Its theoretical development applies a La-grangian approach to motion over a Eulerian grid space by tracking quantities across a landscape as an evolving front. An algorithm for this technique,called level set method,is implemented in a geographical information system(GIS).It ﬁts with aﬁeld data model in GIS and is implemented as operators in map algebra.The paper describes an implemen-tation of the level set methods in a map algebra program-ming language,called MapScript,and gives example pro-gram scripts for applications in ecology and hydrology.1IntroductionOver the past decade there has been an explosion in the ap-plication of models to solve environmental issues.Many of these models are speciﬁc to one physical process and of-ten require expert knowledge to use.Increasingly generic modeling frameworks are being sought to provide analyti-cal tools to examine and resolve complex environmental and natural resource problems.These systems consider a vari-ety of land condition characteristics,interactions and driv-ing physical processes.Variables accounted for include cli-mate,topography,soils,geology,land cover,vegetation and hydro-geography(Moore et al.,1993).Physical interactions include processes for climatology,hydrology,topographic landsurface/sub-surfaceﬂuxes and biological/ecological sys-Correspondence to:D.Pullar(d.pullar@.au)tems(Sklar and Costanza,1991).Progress has been made in linking model-speciﬁc systems with tools used by environ-mental managers,for instance geographical information sys-tems(GIS).While this approach,commonly referred to as loose coupling,provides a practical solution it still does not improve the scientiﬁc foundation of these models nor their integration with other models and related systems,such as decision support systems(Argent,2003).The alternative ap-proach is for tightly coupled systems which build functional-ity into a system or interface to domain libraries from which a user may build custom solutions using a macro language or program scripts.The approach supports integrated models through interface speciﬁcations which articulate the funda-mental assumptions and simpliﬁcations within these models. The problem is that there are no environmental modelling systems which are widely used by engineers and scientists that offer this level of interoperability,and the more com-monly used GIS systems do not currently support space and time representations and operations suitable for modelling environmental processes(Burrough,1998)(Sui and Magio, 1999).Providing a generic environmental modeling framework for practical environmental issues is challenging.It does not exist now despite an overwhelming demand because there are deep technical challenges to build integrated modeling frameworks in a scientiﬁcally rigorous manner.It is this chal-lenge this research addresses.1.1Background for ApproachThe paper describes a generic environmental modeling lan-guage integrated with a Geographical Information System (GIS)which supports spatial-temporal operators to model physical interactions occurring in two ways.The trivial case where interactions are isolated to a location,and the more common and complex case where interactions propa-gate spatially across landscape surfaces.The programming language has a strong theoretical and algorithmic basis.The-oretically,it assumes a Eulerian representation of state space,Fig.1.Shows a)a propagating interface parameterised by differ-ential equations,b)interface fronts have variable intensity and may expand or contract based onﬁeld gradients and driving process. but propagates quantities across landscapes using Lagrangian equations of motion.In physics,a Lagrangian view focuses on how a quantity(water volume or particle)moves through space,whereas an Eulerian view focuses on a localﬁxed area of space and accounts for quantities moving through it.The beneﬁt of this approach is that an Eulerian perspective is em-inently suited to representing the variation of environmen-tal phenomena across space,but it is difﬁcult to conceptu-alise solutions for the equations of motion and has compu-tational drawbacks(Press et al.,1992).On the other hand, the Lagrangian view is often not favoured because it requires a global solution that makes it difﬁcult to account for local variations,but has the advantage of solving equations of mo-tion in an intuitive and numerically direct way.The research will address this dilemma by adopting a novel approach from the image processing discipline that uses a Lagrangian ap-proach over an Eulerian grid.The approach,called level set methods,provides an efﬁcient algorithm for modeling a natural advancing front in a host of settings(Sethian,1999). The reason the method works well over other approaches is that the advancing front is described by equations of motion (Lagrangian view),but computationally the front propagates over a vectorﬁeld(Eulerian view).Hence,we have a very generic way to describe the motion of quantities,but can ex-plicitly solve their advancing properties locally as propagat-ing zones.The research work will adapt this technique for modeling the motion of environmental variables across time and space.Speciﬁcally,it will add new data models and op-erators to a geographical information system(GIS)for envi-ronmental modeling.This is considered to be a signiﬁcant research imperative in spatial information science and tech-nology(Goodchild,2001).The main focus of this paper is to evaluate if the level set method(Sethian,1999)can:–provide a theoretically and empirically supportable methodology for modeling a range of integral landscape processes,–provide an algorithmic solution that is not sensitive to process timing,is computationally stable and efﬁcient as compared to conventional explicit solutions to diffu-sive processes models,–be developed as part of a generic modelling language in GIS to express integrated models for natural resource and environmental problems?The outline for the paper is as follow.The next section will describe the theory for spatial-temporal processing us-ing level sets.Section3describes how this is implemented in a map algebra programming language.Two application examples are given–an ecological and a hydrological ex-ample–to demonstrate the use of operators for computing reactive-diffusive interactions in landscapes.Section4sum-marises the contribution of this research.2Theory2.1IntroductionLevel set methods(Sethian,1999)have been applied in a large collection of applications including,physics,chemistry,ﬂuid dynamics,combustion,material science,fabrication of microelectronics,and computer vision.Level set methods compute an advancing interface using an Eulerian grid and the Lagrangian equations of motion.They are similar to cost distance modeling used in GIS(Burroughs and McDonnell, 1998)in that they compute the spread of a variable across space,but the motion is based upon partial differential equa-tions related to the physical process.The advancement of the interface is computed through time along a spatial gradient, and it may expand or contract in its extent.See Fig.1.2.2TheoryThe advantage of the level set method is that it models mo-tion along a state-space gradient.Level set methods start with the equation of motion,i.e.an advancing front with velocity F is characterised by an arrival surface T(x,y).Note that F is a velocityﬁeld in a spatial sense.If F was constant this would result in an expanding series of circular fronts,but for different values in a velocityﬁeld the front will have a more contorted appearance as shown in Fig.1b.The motion of thisinterface is always normal to the interface boundary,and its progress is regulated by several factors:F=f(L,G,I)(1)where L=local properties that determine the shape of advanc-ing front,G=global properties related to governing forces for its motion,I=independent properties that regulate and inﬂu-ence the motion.If the advancing front is modeled strictly in terms of the movement of entity particles,then a straightfor-ward velocity equation describes its motion:|∇T|F=1given T0=0(2) where the arrival function T(x,y)is a travel cost surface,and T0is the initial position of the interface.Instead we use level sets to describe the interface as a complex function.The level set functionφis an evolving front consistent with the under-lying viscosity solution deﬁned by partial differential equa-tions.This is expressed by the equation:φt+F|∇φ|=0givenφ(x,y,t=0)(3)whereφt is a complex interface function over time period 0..n,i.e.φ(x,y,t)=t0..tn,∇φis the spatial and temporal derivatives for viscosity equations.The Eulerian view over a spatial domain imposes a discretisation of space,i.e.the raster grid,which records changes in value z.Hence,the level set function becomesφ(x,y,z,t)to describe an evolv-ing surface over time.Further details are given in Sethian (1999)along with efﬁcient algorithms.The next section de-scribes the integration of the level set methods with GIS.3Map algebra modelling3.1Map algebraSpatial models are written in a map algebra programming language.Map algebra is a function-oriented language that operates on four implicit spatial data types:point,neighbour-hood,zonal and whole landscape surfaces.Surfaces are typ-ically represented as a discrete raster where a point is a cell, a neighbourhood is a kernel centred on a cell,and zones are groups of mon examples of raster data include ter-rain models,categorical land cover maps,and scalar temper-ature surfaces.Map algebra is used to program many types of landscape models ranging from land suitability models to mineral exploration in the geosciences(Burrough and Mc-Donnell,1998;Bonham-Carter,1994).The syntax for map algebra follows a mathematical style with statements expressed as equations.These equations use operators to manipulate spatial data types for point and neighbourhoods.Expressions that manipulate a raster sur-face may use a global operation or alternatively iterate over the cells in a raster.For instance the GRID map algebra (Gao et al.,1993)deﬁnes an iteration construct,called do-cell,to apply equations on a cell-by-cell basis.This is triv-ially performed on columns and rows in a clockwork manner. However,for environmental phenomena there aresituations Fig.2.Spatial processing orders for raster.where the order of computations has a special signiﬁcance. For instance,processes that involve spreading or transport acting along environmental gradients within the landscape. Therefore special control needs to be exercised on the order of execution.Burrough(1998)describes two extra control mechanisms for diffusion and directed topology.Figure2 shows the three principle types of processing orders,and they are:–row scan order governed by the clockwork lattice struc-ture,–spread order governed by the spreading or scattering ofa material from a more concentrated region,–ﬂow order governed by advection which is the transport of a material due to velocity.Our implementation of map algebra,called MapScript (Pullar,2001),includes a special iteration construct that sup-ports these processing orders.MapScript is a lightweight lan-guage for processing raster-based GIS data using map alge-bra.The language parser and engine are built as a software component to interoperate with the IDRISI GIS(Eastman, 1997).MapScript is built in C++with a class hierarchy based upon a value type.Variants for value types include numeri-cal,boolean,template,cells,or a grid.MapScript supports combinations of these data types within equations with basic arithmetic and relational comparison operators.Algebra op-erations on templates typically result in an aggregate value assigned to a cell(Pullar,2001);this is similar to the con-volution integral in image algebras(Ritter et al.,1990).The language supports iteration to execute a block of statements in three ways:a)docell construct to process raster in a row scan order,b)dospread construct to process raster in a spreadwhile(time<100)dospreadpop=pop+(diffuse(kernel*pop))pop=pop+(r*pop*dt*(1-(pop/K)) enddoendwhere the diffusive constant is stored in thekernel:Fig.3.Map algebra script and convolution kernel for population dispersion.The variable pop is a raster,r,K and D are constants, dt is the model time step,and the kernel is a3×3template.It is assumed a time step is deﬁned and the script is run in a simulation. Theﬁrst line contained in the nested cell processing construct(i.e. dospread)is the diffusive term and the second line is the population growth term.order,c)doﬂow to process raster byﬂow order.Examples are given in subsequent sections.Process models will also involve a timing loop which may be handled as a general while(<condition>)..end construct in MapScript where the condition expression includes a system time variable.This time variable is used in a speciﬁc fashion along with a system time step by certain operators,namely diffuse()andﬂuxﬂow() described in the next section,to model diffusion and advec-tion as a time evolving front.The evolving front represents quantities such as vegetation growth or surface runoff.3.2Ecological exampleThis section presents an ecological example based upon plant dispersal in a landscape.The population of a species follows a controlled growth rate and at the same time spreads across landscapes.The theory of the rate of spread of an organism is given in Tilman and Kareiva(1997).The area occupied by a species grows log-linear with time.This may be modelled by coupling a spatial diffusion term with an exponential pop-ulation growth term;the combination produces the familiar reaction-diffusion model.A simple growth population model is used where the reac-tion term considers one population controlled by births and mortalities is:dN dt =r·N1−NK(4)where N is the size of the population,r is the rate of change of population given in terms of the difference between birth and mortality rates,and K is the carrying capacity.Further dis-cussion of population models can be found in Jrgensen and Bendoricchio(2001).The diffusive term spreads a quantity through space at a speciﬁed rate:dudt=Dd2udx2(5) where u is the quantity which in our case is population size, and D is the diffusive coefﬁcient.The model is operated as a coupled computation.Over a discretized space,or raster,the diffusive term is estimated using a numerical scheme(Press et al.,1992).The distance over which diffusion takes place in time step dt is minimally constrained by the raster resolution.For a stable computa-tional process the following condition must be satisﬁed:2Ddtdx2≤1(6) This basically states that to account for the diffusive pro-cess,the term2D·dx is less than the velocity of the advancing front.This would not be difﬁcult to compute if D is constant, but is problematic if D is variable with respect to landscape conditions.This problem may be overcome by progressing along a diffusive front over the discrete raster based upon distance rather than being constrained by the cell resolution.The pro-cessing and diffusive operator is implemented in a map al-gebra programming language.The code fragment in Fig.3 shows a map algebra script for a single time step for the cou-pled reactive-diffusion model for population growth.The operator of interest in the script shown in Fig.3is the diffuse operator.It is assumed that the script is run with a given time step.The operator uses a system time step which is computed to balance the effect of process errors with efﬁ-cient computation.With knowledge of the time step the it-erative construct applies an appropriate distance propagation such that the condition in Eq.(3)is not violated.The level set algorithm(Sethian,1999)is used to do this in a stable and accurate way.As a diffusive front propagates through the raster,a cost distance kernel assigns the proper time to each raster cell.The time assigned to the cell corresponds to the minimal cost it takes to reach that cell.Hence cell pro-cessing is controlled by propagating the kernel outward at a speed adaptive to the local context rather than meeting an arbitrary global constraint.3.3Hydrological exampleThis section presents a hydrological example based upon sur-face dispersal of excess rainfall across the terrain.The move-ment of water is described by the continuity equation:∂h∂t=e t−∇·q t(7) where h is the water depth(m),e t is the rainfall excess(m/s), q is the discharge(m/hr)at time t.Discharge is assumed to have steady uniformﬂow conditions,and is determined by Manning’s equation:q t=v t h t=1nh5/3ts1/2(8)putation of current cell(x+ x,t,t+ ).where q t is theﬂow velocity(m/s),h t is water depth,and s is the surface slope(m/m).An explicit method of calcula-tion is used to compute velocity and depth over raster cells, and equations are solved at each time step.A conservative form of aﬁnite difference method solves for q t in Eq.(5). To simplify discussions we describe quasi-one-dimensional equations for theﬂow problem.The actual numerical com-putations are normally performed on an Eulerian grid(Julien et al.,1995).Finite-element approximations are made to solve the above partial differential equations for the one-dimensional case ofﬂow along a strip of unit width.This leads to a cou-pled model with one term to maintain the continuity ofﬂow and another term to compute theﬂow.In addition,all calcu-lations must progress from an uphill cell to the down slope cell.This is implemented in map algebra by a iteration con-struct,called doﬂow,which processes a raster byﬂow order. Flow distance is measured in cell size x per unit length. One strip is processed during a time interval t(Fig.4).The conservative solution for the continuity term using aﬁrst or-der approximation for Eq.(5)is derived as:h x+ x,t+ t=h x+ x,t−q x+ x,t−q x,txt(9)where the inﬂow q x,t and outﬂow q x+x,t are calculated in the second term using Equation6as:q x,t=v x,t·h t(10) The calculations approximate discharge from previous time interval.Discharge is dynamically determined within the continuity equation by water depth.The rate of change in state variables for Equation6needs to satisfy a stability condition where v· t/ x≤1to maintain numerical stabil-ity.The physical interpretation of this is that aﬁnite volume of water wouldﬂow across and out of a cell within the time step t.Typically the cell resolution isﬁxed for the raster, and adjusting the time step requires restarting the simulation while(time<120)doflow(dem)fvel=1/n*pow(depth,m)*sqrt(grade)depth=depth+(depth*fluxflow(fvel)) enddoendFig.5.Map algebra script for excess rainfallﬂow computed over a 120minute event.The variables depth and grade are rasters,fvel is theﬂow velocity,n and m are constants in Manning’s equation.It is assumed a time step is deﬁned and the script is run in a simulation. Theﬁrst line in the nested cell processing(i.e.doﬂow)computes theﬂow velocity and the second line computes the change in depth from the previous value plus any net change(inﬂow–outﬂow)due to velocityﬂux across the cell.cycle.Flow velocities change dramatically over the course of a storm event,and it is problematic to set an appropriate time step which is efﬁcient and yields a stable result.The hydrological model has been implemented in a map algebra programming language Pullar(2003).To overcome the problem mentioned above we have added high level oper-ators to compute theﬂow as an advancing front over a land-scape.The time step advances this front adaptively across the landscape based upon theﬂow velocity.The level set algorithm(Sethian,1999)is used to do this in a stable and accurate way.The map algebra script is given in Fig.5.The important operator is theﬂuxﬂow operator.It computes the advancing front for waterﬂow across a DEM by hydrologi-cal principles,and computes the local drainageﬂux rate for each cell.Theﬂux rate is used to compute the net change in a cell in terms ofﬂow depth over an adaptive time step.4ConclusionsThe paper has described an approach to extend the function-ality of tightly coupled environmental models in GIS(Ar-gent,2004).A long standing criticism of GIS has been its in-ability to handle dynamic spatial models.Other researchers have also addressed this issue(Burrough,1998).The con-tribution of this paper is to describe how level set methods are:i)an appropriate scientiﬁc basis,and ii)able to perform stable time-space computations for modelling landscape pro-cesses.The level set method provides the following beneﬁts:–it more directly models motion of spatial phenomena and may handle both expanding and contracting inter-faces,–is based upon differential equations related to the spatial dynamics of physical processes.Despite the potential for using level set methods in GIS and land-surface process modeling,there are no commercial or research systems that use this mercial sys-tems such as GRID(Gao et al.,1993),and research systems such as PCRaster(Wesseling et al.,1996)offerﬂexible andpowerful map algebra programming languages.But opera-tions that involve reaction-diffusive processing are speciﬁc to one context,such as groundwaterﬂow.We believe the level set method offers a more generic approach that allows a user to programﬂow and diffusive landscape processes for a variety of application contexts.We have shown that it pro-vides an appropriate theoretical underpinning and may be ef-ﬁciently implemented in a GIS.We have demonstrated its application for two landscape processes–albeit relatively simple examples–but these may be extended to deal with more complex and dynamic circumstances.The validation for improved environmental modeling tools ultimately rests in their uptake and usage by scientists and engineers.The tool may be accessed from the web site .au/projects/mapscript/(version with enhancements available April2005)for use with IDRSIS GIS(Eastman,1997)and in the future with ArcGIS. It is hoped that a larger community of users will make use of the methodology and implementation for a variety of environmental modeling applications.Edited by:P.Krause,S.Kralisch,and W.Fl¨u gelReviewed by:anonymous refereesReferencesArgent,R.:An Overview of Model Integration for Environmental Applications,Environmental Modelling and Software,19,219–234,2004.Bonham-Carter,G.F.:Geographic Information Systems for Geo-scientists,Elsevier Science Inc.,New York,1994. Burrough,P.A.:Dynamic Modelling and Geocomputation,in: Geocomputation:A Primer,edited by:Longley,P.A.,et al., Wiley,England,165-191,1998.Burrough,P.A.and McDonnell,R.:Principles of Geographic In-formation Systems,Oxford University Press,New York,1998. Gao,P.,Zhan,C.,and Menon,S.:An Overview of Cell-Based Mod-eling with GIS,in:Environmental Modeling with GIS,edited by: Goodchild,M.F.,et al.,Oxford University Press,325–331,1993.Goodchild,M.:A Geographer Looks at Spatial Information Theory, in:COSIT–Spatial Information Theory,edited by:Goos,G., Hertmanis,J.,and van Leeuwen,J.,LNCS2205,1–13,2001.Jørgensen,S.and Bendoricchio,G.:Fundamentals of Ecological Modelling,Elsevier,New York,2001.Julien,P.Y.,Saghaﬁan,B.,and Ogden,F.:Raster-Based Hydro-logic Modelling of Spatially-Varied Surface Runoff,Water Re-sources Bulletin,31(3),523–536,1995.Moore,I.D.,Turner,A.,Wilson,J.,Jenson,S.,and Band,L.:GIS and Land-Surface-Subsurface Process Modeling,in:Environ-mental Modeling with GIS,edited by:Goodchild,M.F.,et al., Oxford University Press,New York,1993.Press,W.,Flannery,B.,Teukolsky,S.,and Vetterling,W.:Numeri-cal Recipes in C:The Art of Scientiﬁc Computing,2nd Ed.Cam-bridge University Press,Cambridge,1992.Pullar,D.:MapScript:A Map Algebra Programming Language Incorporating Neighborhood Analysis,GeoInformatica,5(2), 145–163,2001.Pullar,D.:Simulation Modelling Applied To Runoff Modelling Us-ing MapScript,Transactions in GIS,7(2),267–283,2003. Ritter,G.,Wilson,J.,and Davidson,J.:Image Algebra:An Overview,Computer Vision,Graphics,and Image Processing, 4,297–331,1990.Sethian,J.A.:Level Set Methods and Fast Marching Methods, Cambridge University Press,Cambridge,1999.Sklar,F.H.and Costanza,R.:The Development of Dynamic Spa-tial Models for Landscape Ecology:A Review and Progress,in: Quantitative Methods in Ecology,Springer-Verlag,New York, 239–288,1991.Sui,D.and R.Maggio:Integrating GIS with Hydrological Mod-eling:Practices,Problems,and Prospects,Computers,Environ-ment and Urban Systems,23(1),33–51,1999.Tilman,D.and P.Kareiva:Spatial Ecology:The Role of Space in Population Dynamics and Interspeciﬁc Interactions.Princeton University Press,Princeton,New Jersey,USA,1997. Wesseling C.G.,Karssenberg, D.,Burrough,P. A.,and van Deursen,W.P.:Integrating Dynamic Environmental Models in GIS:The Development of a Dynamic Modelling Language, Transactions in GIS,1(1),40–48,1996.。

minimalmarkerset...

Minimal Marker Sets to Discriminate Among Seedlines Thomas C.Hudson,Ann E.Stapleton,Amy M.CurleyUniversity of North Carolina at WilmingtonDepartments of Computer Science and Biological Sciences{hudsont,stapletona}@AbstractRaising seeds for biological experiments is prone to er-ror;a careful experimenter will test in the lab to verify that plants are of the intended strain.Choosing a minimal set of tests that will discriminate between all known seedlines is an instance of Minimal Test Set,a NP-complete problem. Similar biological problems,such as minimizing the num-ber of haplotype tag SNPs,require complex nondeterminis-tic heuristics to solve in reasonable timeframes over modest datasets.However,selecting the minimal marker set to dis-criminate among seedlines is less complicated than other problems considered in the literature;we show that a sim-ple heuristic approach works well in practice.Finding all minimal sets of tests to identify91Zea mays recombinant inbred lines would require months of CPU time;our heuris-tic gives a result less than twice the minimal possible size in underﬁve seconds,with similar performance on Arabidop-sis thaliana recombinant inbred lines.1.IntroductionWhen a plant geneticist wants to conduct an experiment, she needs samples of a plant.Frequently,she will grow the plant herself from seeds kept in her laboratory.However, raising these plants is a labor-intensive,error-prone proce-dure:seeds can be wrongly sown,ﬁelds wrongly marked, natural pollination occur unintentionally,collected seeds mislabelled in theﬁeld or stored incorrectly in the lab.A cautious scientist will perform tests on the plants she takes her experimental samples from to conﬁrm that they are from the intended seedline.To verify the genotype of the sample,the scientist selects markers,extracts DNA from the sample plants,and ampli-ﬁes each test region;these regions have known detectable differences in length.In the case of recombinant inbred lines,there are only two possibilities for each marker,con-ventionally referred to as size“A”and“B”.Our poster reports on heuristic algorithms developed to help minimize the expense of testing.Finding the optimum set of markers to use is a problem that can take months or years of CPU time;this software produces near-optimum answers in under a minute.The algorithms discussed in our poster have been im-plemented in Java and are available under an open-source license at /csc/bioinformatics/.2.Heuristic SolutionA randomized greedy algorithm gives a reasonableﬁrst answer for the problem ofﬁnding minimal marker sets to distinguish among the seedlines:1.Shufﬂe the markers into a random order2.Examine each marker in order(a)Remove it from the set of markers if the resul-tant set is still able to discriminate among all theseedlinesIn our experiments on Zea mays(134markers,91seed-lines)and Arabidopsis thaliana(99markers,162seedlines), this random greedy approach produces answers no more than twice the size of the theoretical optimum;repeated tri-als show that the results are roughly normally distributed (see our poster).If there are N seedlines and M markers, the theoretical complexity is O(M2N2);the algorithm runs in seconds on those datasets.These distributions imply that random sampling of the search space could yield reasonable results.The quality of the result of random sampling is very sensitive to the input: some subsets of the full data have many minimal-length an-swers,making random discovery likely,while others have only one.However,in practice they seem to have a large number of solutions requiring one marker more than min-imal,which are reasonably likely to be found by random search.As problems grow larger–more seedlines are de-veloped and more markers are identiﬁed–larger and largersamples of the search space will be necessary to have a rea-sonable likelihood ofﬁnding a good solution.Sorting according to simple metrics does not yield anyimprovement on random ordering,but provides consistency.Assigning a large negative value to a marker for every seed-line about which the marker returns an inconclusive resultgives a coarse ordering.If A and B appear with dissimilarfrequency,adding a small positive value to the marker’s rat-ing for every seedline on which it returns the less-commonresult gives aﬁner ordering.Neither of these metrics out-performs random ordering;both typically give a result com-parable to the median result returned in one thousand tri-als of random ordering.However,they do so in a singletrial(underﬁve seconds for both Zea mays and Arabidopsisthaliana),which gives us good input for the second stage of our algorithm.We thenﬁlter the data.If the initial greedy heuristic re-turns a solution S containing K markers,we run the greedyalgorithm K additional times.Let S i be the i th marker in S;on the i th additional execution in thisﬁltering pass,we remove S i from the set of possible markers.Whether we start with a random or sorted list of mark-ers,running the basic greedy algorithm and then one passofﬁltering gives us an answer of the same size as the bestanswers ever returned by the randomized algorithm.Addi-tional passes of theﬁltering algorithm do not yield furtherimprovement.For both Zea mays and Arabidopsis thaliana,this is roughly one pointﬁve times the length of the smallestpossible answer.In essence,this algorithm performs a heuristic searchof the M-dimensional space of possible answers toﬁnda candidate answer,and then exhaustively explores its K-dimensional immediate neighborhood looking for a localminimum.Weﬁnd that the solution initially reported bythe greedy heuristic is rarely a local minimum,but that itconsistently has an adjacent local minimum.Over the datacurrently available,the two-stage approach gives reliablygood results about a minute.3.Exact SolutionA heuristic solution to the problem is not strictly nec-essary.The minimal discriminating set of markers canbe found by examining all potentially discriminating sets.However,this requires an exhaustive search over a largesearch space.For N seedlines and M markers,there areMJsub-sets of markers of size J.For each subset that we exam-ine,a straightforward determination of whether the subset distinguishes between each pair of seedlines takes O(JN2) time.The total predicted time is O(M K N2K),where K is the size of the minimal discriminating marker set;K>=log2(N).To verify this O()characterization,we implemented an exact solver for the minimal discriminating marker set prob-lem and ran it over subsets of the Zea mays data.Graphs of the time performance of the exact solver can be found on the poster.A trial run of the exact solver on a dedicated 2.4GHz Xeon CPU examined only1.33%of the possible size-7solutions for Zea mays in17.5CPU hours;if there isa size7answer,it would take us54days toﬁnd.4.Theory and ContextFinding the minimal discriminating set of markers is an instance of a well-known NP-complete problem,Minimal Test Set[1].In Garey and Johnson’s formulation,the asso-ciated decision problem is:INSTANCE:A collection C of subsets of aﬁniteset S,a positive integer K≤||C||.QUESTION:Is there a subcollection C ⊆Cwith||C ||<K such that for each pair of distinctelements u,v∈S,there exists some set c∈Cthat contains exactly one of u and v?[3]is a comprehensive survey of approaches to the Minimal Test Set problem.This problem looks similar to another question intensely studied in bioinformatics,Haplotype Tag Selection.Al-though the decision problems are only subtly different,this difference signiﬁcantly increases the complexity of algo-rithms that solve Haplotype Tag Selection.Approaches like ours to Minimal Test Set are not sufﬁcient to solve Haplo-type Tag Selection.[2]is a survey of current work on Haplotype Tag Selec-tion.The authorsﬁt21published Haplotype Tag Selection algorithms into a three-stage framework:evaluating each SNP for how well it describes other nearby SNPs,evaluat-ing a candidate set of SNPs for how well they classify the entire set of data,and constructing aﬁnal minimal set of SNPs.Our algorithm performs three analogous activities, albeit in a different order:ﬁltering to minimize the set of results,sorting metrics,and a greedy minimization phase.References[1]M.R.Garey and puters and Intractability:A Guide to the Theory of NP-Completeness.W.H.Freemanand Company,San Francisco,1979.[2] B.V.Halldorsson,S.Istrail,and F.M.De La Vega.Opti-mal selection of snp markers for disease association studies.Human Heredity,58:190–202,2004.[3] B.Moret and H.Shapiro.On minimizing a set of tests.SIAMJournal of Scientiﬁc Computing,6(4):983–1003,1985.Below is given annual work summary, do not need friends can download after editor deleted Welcome to visit againXXXX annual work summaryDear every leader, colleagues:Look back end of XXXX, XXXX years of work, have the joy of success in your work, have a collaboration with colleagues, working hard, also have disappointed when encountered difficulties and setbacks. Imperceptible in tense and orderly to be over a year, a year, under the loving care and guidance of the leadership of the company, under the support and help of colleagues, through their own efforts, various aspects have made certain progress, better to complete the job. For better work, sum up experience and lessons, will now work a brief summary.To continuously strengthen learning, improve their comprehensive quality. With good comprehensive quality is the precondition of completes the labor of duty and conditions. A year always put learning in the important position, trying to improve their comprehensive quality. Continuous learning professional skills, learn from surrounding colleagues with rich work experience, equip themselves with knowledge, the expanded aspect of knowledge, efforts to improve their comprehensive quality.The second Do best, strictly perform their responsibilities. Set up the company, to maximize the customer to the satisfaction of the company's products, do a good job in technical services and product promotion to the company. And collected on the properties of the products of the company, in order to make improvement in time, make the products better meet the using demand of the scene.Three to learn to be good at communication, coordinating assistance. On‐site technical service personnel should not only have strong professional technology, should also have good communication ability, a lot of a product due to improper operation to appear problem, but often not customers reflect the quality of no, so this time we need to find out the crux, and customer communication, standardized operation, to avoid customer's mistrust of the products and even the damage of the company's image. Some experiences in the past work, mentality is very important in the work, work to have passion, keep the smile of sunshine, can close the distance between people, easy to communicate with the customer. Do better in the daily work to communicate with customers and achieve customer satisfaction, excellent technical service every time, on behalf of the customer on our products much a understanding and trust.Fourth, we need to continue to learn professional knowledge, do practical grasp skilled operation. Over the past year, through continuous learning and fumble, studied the gas generation, collection and methods, gradually familiar with and master the company introduced the working principle, operation method of gas machine. With the help of the department leaders and colleagues, familiar with and master the launch of the division principle, debugging method of the control system, and to wuhan Chen Guchong garbage power plant of gas machine control system transformation, learn to debug, accumulated some experience. All in all, over the past year, did some work, have also made some achievements, but the results can only represent the past, there are some problems to work, can't meet the higher requirements. In the future work, I must develop the oneself advantage, lack of correct, foster strengths and circumvent weaknesses, for greater achievements. Looking forward to XXXX years of work, I'll be more efforts, constant progress in their jobs, make greater achievements. Every year I have progress, the growth of believe will get greater returns, I will my biggest contribution to the development of the company, believe inyourself do better next year!I wish you all work study progress in the year to come.。

经验风险最小化在特征选择中的应用

经验风险最小化在特征选择中的应用在机器学习和数据挖掘领域，特征选择是一个重要的任务，它用于从原始数据中选择出最具有代表性的特征，以提高模型的性能和减少计算复杂度。

经验风险最小化（ERM）是一种常用的优化方法，它通过最小化经验风险来选择最佳的模型参数。

本文将探讨经验风险最小化在特征选择中的应用，并介绍一些常用的特征选择算法。

首先，让我们了解一下经验风险最小化的基本概念。

经验风险是指模型在训练集上的平均误差，它可以用来衡量模型的拟合能力。

最小化经验风险意味着选择能够在训练集上表现最好的模型参数。

然而，仅仅依靠经验风险来选择模型参数可能会导致过拟合的问题，即模型在训练集上表现良好，但在测试集上表现较差。

为了解决这个问题，我们需要引入正则化项来惩罚模型的复杂度。

在特征选择中，经验风险最小化可以用来选择最具有代表性的特征子集。

特征选择的目标是去除冗余和噪声特征，以提高模型的泛化能力。

一种常用的特征选择方法是基于过滤的方法，它通过计算每个特征与目标变量之间的相关性来选择特征。

具体而言，我们可以使用皮尔逊相关系数、互信息等指标来衡量特征与目标变量之间的关联程度。

然后，我们可以根据这些指标来选择具有最高相关性的特征。

除了基于过滤的方法，还有一些基于包装的特征选择方法，它们通过训练模型来评估特征的重要性。

这些方法通常使用交叉验证技术来评估模型的性能，并根据模型的性能来选择特征。

其中，递归特征消除（RFE）是一种常用的包装方法，它通过反复训练模型并剔除最不重要的特征来选择特征子集。

RFE的核心思想是，如果一个特征对模型的性能影响较小，那么剔除这个特征不会对模型的性能造成太大的影响。

此外，基于嵌入的特征选择方法也是一种常见的选择方法。

嵌入方法将特征选择与模型训练过程融合在一起，通过优化模型的目标函数来选择特征。

经典的嵌入方法包括L1正则化和决策树算法。

L1正则化可以实现稀疏性，即将某些特征的权重设为0，从而实现特征选择的目的。

Apprenticeship learning via inverse reinforcement learning

Apprenticeship Learning via Inverse Reinforcement Learning Pieter Abbeel pabbeel@Andrew Y.Ng ang@ Computer Science Department,Stanford University,Stanford,CA94305,USAAbstractWe consider learning in a Markov decisionprocess where we are not explicitly given a re-ward function,but where instead we can ob-serve an expert demonstrating the task thatwe want to learn to perform.This settingis useful in applications(such as the task ofdriving)where it may be diﬃcult to writedown an explicit reward function specifyingexactly how diﬀerent desiderata should betraded oﬀ.We think of the expert as try-ing to maximize a reward function that is ex-pressible as a linear combination of knownfeatures,and give an algorithm for learningthe task demonstrated by the expert.Our al-gorithm is based on using“inverse reinforce-ment learning”to try to recover the unknownreward function.We show that our algorithmterminates in a small number of iterations,and that even though we may never recoverthe expert’s reward function,the policy out-put by the algorithm will attain performanceclose to that of the expert,where here per-formance is measured with respect to the ex-pert’s unknown reward function.1.IntroductionGiven a sequential decision making problem posed in the Markov decision process(MDP)formalism,a num-ber of standard algorithms exist forﬁnding an optimal or near-optimal policy.In the MDP setting,we typi-cally assume that a reward function is given.Given a reward function and the MDPs state transition prob-abilities,the value function and optimal policy are ex-actly determined.The MDP formalism is useful for many problems be-cause it is often easier to specify the reward function than to directly specify the value function(and/or op-timal policy).However,we believe that even the re-ward function is frequently diﬃcult to specify manu-ally.Consider,for example,the task of highway driv-ing.When driving,we typically trade oﬀmany dif-Appearing in Proceedings of the21st International Confer-ence on Machine Learning,Banﬀ,Canada,2004.Copyright2004by the authors.ferent desiderata,such as maintaining safe following distance,keeping away from the curb,staying far from any pedestrians,maintaining a reasonable speed,per-haps a slight preference for driving in the middle lane, not changing lanes too often,and so on....To specify a reward function for the driving task,we would have to assign a set of weights stating exactly how we would like to trade oﬀthese diﬀerent factors.Despite being able to drive competently,the authors do not believe they can conﬁdently specify a speciﬁc reward function for the task of“driving well.”1In practice,this means that the reward function is of-ten manually tweaked(cf.reward shaping,Ng et al., 1999)until the desired behavior is obtained.From con-versations with engineers in industry and our own ex-perience in applying reinforcement learning algorithms to several robots,we believe that,for many problems, the diﬃculty of manually specifying a reward function represents a signiﬁcant barrier to the broader appli-cability of reinforcement learning and optimal control algorithms.When teaching a young adult to drive,rather than telling them what the reward function is,it is much easier and more natural to demonstrate driving to them,and have them learn from the demonstration. The task of learning from an expert is called appren-ticeship learning(also learning by watching,imitation learning,or learning from demonstration).A number of approaches have been proposed for ap-prenticeship learning in various applications.Most of these methods try to directly mimic the demonstrator by applying a supervised learning algorithm to learn a direct mapping from the states to the actions.This literature is too wide to survey here,but some ex-amples include Sammut et al.(1992);Kuniyoshi et al.(1994);Demiris&Hayes(1994);Amit&Mataric (2002);Pomerleau(1989).One notable exception is given in Atkeson&Schaal(1997).They considered the 1We note that this is true even though the reward func-tion may often be easy to state in English.For instance, the“true”reward function that we are trying to maximize when driving is,perhaps,our“personal happiness.”Thepractical problem however is how to model this(i.e.,our happiness)explicitly as a function of the problems’states, so that a reinforcement learning algorithm can be applied.problem of having a robot arm follow a demonstrated trajectory,and used a reward function that quadrat-ically penalizes deviation from the desired trajectory. Note however,that this method is applicable only to problems where the task is to mimic the expert’s tra-jectory.For highway driving,blindly following the ex-pert’s trajectory would not work,because the pattern of traﬃc encountered is diﬀerent each time.Given that the entireﬁeld of reinforcement learning is founded on the presupposition that the reward func-tion,rather than the policy or the value function,is the most succinct,robust,and transferable deﬁnition of the task,it seems natural to consider an approach to apprenticeship learning whereby the reward function is learned.2The problem of deriving a reward function from ob-served behavior is referred to as inverse reinforcement learning(Ng&Russell,2000).In this paper,we assume that the expert is trying(without necessar-ily succeeding)to optimize an unknown reward func-tion that can be expressed as a linear combination of known“features.”Even though we cannot guarantee that our algorithms will correctly recover the expert’s true reward function,we show that our algorithm will nonethelessﬁnd a policy that performs as well as the expert,where performance is measured with respect to the expert’s unknown reward function.2.PreliminariesA(ﬁnite-state)Markov decision process(MDP)is a tu-ple(S,A,T,γ,D,R),where S is aﬁnite set of states;A is a set of actions;T={P sa}is a set of state transition probabilities(here,P sa is the state transition distribu-tion upon taking action a in state s);γ∈[0,1)is a discount factor;D is the initial-state distribution,from which the start state s0is drawn;and R:S→A is the reward function,which we assume to be bounded in absolute value by1.We let MDP\R denote an MDP without a reward function,i.e.,a tuple of the form(S,A,T,γ,D).We assume that there is some vector of featuresφ: S→[0,1]k over states,and that there is some“true”reward function R∗(s)=w∗·φ(s),where w∗∈R k.3 2A related idea is also seen in the biomechanics and cog-nitive science,where researchers have pointed out that sim-ple reward functions(usually ones constructed by hand)of-ten suﬃce to explain complicated behavior(policies).Ex-amples include the minimum jerk principle to explain limb movement in primates(Hogan,1984),and the minimum torque-change model to explain trajectories in human mul-tijoint arm movement.(Uno et al.,1989)Related examples are also found in economics and some other literatures. (See the discussion in Ng&Russell,2000.)3The case of state-action rewards R(s,a)oﬀers no ad-ditional diﬃculties;using features of the formφ:S×A→[0,1]k,and our algorithms still apply straightforwardly.In order to ensure that the rewards are bounded by1,we also assume w∗ 1≤1.In the driving domain,φmight be a vector of features indicating the diﬀerent desiderata in driving that we would like to trade oﬀ, such as whether we have just collided with another car,whether we’re driving in the middle lane,and so on.The(unknown)vector w∗speciﬁes the relative weighting between these desiderata.A policyπis a mapping from states to probability distributions over actions.The value of a policyπisE s0∼D[Vπ(s0)]=E[ ∞t=0γt R(s t)|π](1)=E[ ∞t=0γt w·φ(s t)|π](2)=w·E[ ∞t=0γtφ(s t)|π](3) Here,the expectation is taken with respect to the ran-dom state sequence s0,s1,...drawn by starting from a state s0∼D,and picking actions according toπ. We deﬁne the expected discounted accumulated fea-ture value vectorµ(π),or more succinctly the feature expectations,to beµ(π)=E[ ∞t=0γtφ(s t)|π]∈R k.(4) Using this notation,the value of a policy may be writ-ten E s0∼D[Vπ(s0)]=w·µ(π).Given that the reward R is expressible as a linear combination of the fea-turesφ,the feature expectations for a given policyπcompletely determine the expected sum of discounted rewards for acting according to that policy.LetΠdenote the set of stationary policies for an MDP. Given two policiesπ1,π2∈Π,we can construct a new policyπ3by mixing them together.Speciﬁcally,imag-ine thatπ3operates byﬂipping a coin with biasλ,and with probabilityλpicks and always acts according to π1,and with probability1−λalways acts according toπ2.From linearity of expectation,clearly we have thatµ(π3)=λµ(π1)+(1−λ)µ(π2).Note that the randomization step selecting betweenπ1andπ2occurs only once at the start of a trajectory,and not on ev-ery step taken in the MDP.More generally,if we have found some set of policiesπ1,...,πd,and want toﬁnd a new policy whose feature expectations vector is a convex combination n i=1λiµ(πi)(λi≥0, iλi=1) of these policies’,then we can do so by mixing together the policiesπ1,...,πd,where the probability of picking πi is given byλi.We assume access to demonstrations by some expert πE.Speciﬁcally,we assume the ability to observe trajectories(state sequences)generated by the expert starting from s0∼D and taking actions according to πE.It may be helpful to think of theπE as the optimal policy under the reward function R∗=w∗Tφ,though we do not require this to hold.For our algorithm,we will require an estimate of the expert’s feature expectationsµE=µ(πE).Speciﬁ-cally,given a set of m trajectories{s(i)0,s(i)1,...}m i=1 generated by the expert,we denote the empirical esti-mate forµE by4ˆµE=1m m i=1 ∞t=0γtφ(s(i)t).(5)In the sequel,we also assume access to a reinforcement learning(RL)algorithm that can be used to solve an MDP\R augmented with a reward function R=w Tφ. For simplicity of exposition,we will assume that the RL algorithm returns the optimal policy.The general-ization to approximate RL algorithms oﬀers no special diﬃculties;see the full paper.(Abbeel&Ng,2004) 3.AlgorithmThe problem is the following:Given an MDP\R,a feature mappingφand the expert’s feature expecta-tionsµE,ﬁnd a policy whose performance is close to that of the expert’s,on the unknown reward function R∗=w∗Tφ.To accomplish this,we willﬁnd a policy ˜πsuch that µ(˜π)−µE 2≤ .For such a˜π,we would have that for any w∈R k( w 1≤1),|E[ ∞t=0γt R(s t)|πE]−E[ ∞t=0γt R(s t)|˜π]|(6) =|w Tµ(˜π)−w TµE|(7)≤ w 2 µ(˜π)−µE 2(8)≤1· = (9) Theﬁrst inequality follows from the fact that|x T y|≤ x 2 y 2,and the second from w 2≤ w 1≤1. So the problem is reduced toﬁnding a policy˜πthat induces feature expectationsµ(˜π)close toµE.Our apprenticeship learning algorithm forﬁnding such a policy˜πis as follows:1.Randomly pick some policyπ(0),compute(or approx-imate via Monte Carlo)µ(0)=µ(π(0)),and set i=1.pute t(i)=max w: w2≤1min j∈{0..(i−1)}w T(µE−µ(j)),and let w(i)be the value of w that attains this maximum.3.If t(i)≤ ,then terminate.ing the RL algorithm,compute the optimal policyπ(i)for the MDP using rewards R=(w(i))Tφ.pute(or estimate)µ(i)=µ(π(i)).6.Set i=i+1,and go back to step2.Upon termination,the algorithm returns{π(i):i= 0...n}.Let us examine the algorithm in detail.On iteration i,we have already found some policiesπ(0),...,π(i−1). The optimization in step2can be viewed as an inverse reinforcement learning step in which we try to guess 4In practice we truncate the trajectories after aﬁnite number H of steps.If H=H =logγ( (1−γ))is the -horizon time,then this introduces at most error into the approximation.µ(π(0))µ)(2)Figure1.Three iterations for max-margin algorithm. the reward function being optimized by the expert. The maximization in that step is equivalently written max t,w t(10) s.t.w TµE≥w Tµ(j)+t,j=0,...,i−1(11) ||w||2≤1(12) From Eq.(11),we see the algorithm is trying to ﬁnd a reward function R=w(i)·φsuch thatE s0∼D[VπE(s0)]≥E s0∼D[Vπ(i)(s0)]+t.I.e.,a re-ward on which the expert does better,by a“margin”of t,than any of the i policies we had found previously. This step is similar to one used in(Ng&Russell,2000), but unlike the algorithms given there,because of the 2-norm constraint on w it cannot be posed as a linear program(LP),but only as a quadratic program.5 Readers familiar with support vector machines (SVMs)will also recognize this optimization as be-ing equivalent toﬁnding the maximum margin hyper-plane separating two sets of points.(Vapnik,1998)The equivalence is obtained by associating a label1with the expert’s feature expectationsµE,and a label−1 with the feature expectations{µ(π(j)):j=0..(i−1)}. The vector w(i)we want is the unit vector orthogonal to the maximum margin separating hyperplane.So,an SVM solver can also be used toﬁnd w(i).(The SVM problem is a quadratic programming problem(QP),so we can also use any generic QP solver.)In Figure1we show an example of what theﬁrst three iterations of the algorithm could look like geo-metrically.Shown are several exampleµ(π(i)),and the w(i)’s given by diﬀerent iterations of the algorithm. Now,suppose the algorithm terminates,with t(n+1)≤ .(Whether the algorithm terminates is discussed in Section4.)Then directly from Eq.(10-12)we have:∀w with w 2≤1∃i s.t.w Tµ(i)≥w TµE− .(13) Since||w∗||2≤||w∗||1≤1,this means that there is at least one policy from the set returned by the al-gorithm,whose performance under R∗is at least as good as the expert’s performance minus .Thus,at this stage,we can ask the agent designer to manually test/examine the policies found by the algorithm,and 5Although we previously assumed that the w∗specify-ing the“true”rewards statisfy w∗ 1≤1(and our theo-retical results will use this assumption),we still implement the algorithm using w 2≤1,as in Eq.(12).pick one with acceptable performance.A slight exten-sion of this method ensures that the agent designer has to examine at most k+1,rather than all n+1, diﬀerent policies(see footnote6).If we do not wish to ask for human help to select a policy,alternatively we canﬁnd the point closest to µE in the convex closure ofµ(0),...,µ(n)by solving the following QP:min||µE−µ||2,s.t.µ= iλiµ(i),λi≥0, iλi=1. BecauseµE is“separated”from the pointsµ(i)by a margin of at most ,we know that for the solutionµwe have||µE−µ||2≤ .Further,by“mixing”together the policiesπ(i)according to the mixture weightsλi as discussed previously,we obtain a policy whose feature expectations are given byµ.Following our previous discussion(Eq.6-9),this policy attains performance near that of the expert’s on the unknown reward func-tion.6Note that although we called one step of our algorithm an inverse RL step,our algorithm does not necessarily recover the underlying reward function correctly.The performance guarantees of our algorithm only depend on(approximately)matching the feature expectations, not on recovering the true underlying reward function.3.1.A simpler algorithmThe algorithm described above requires access to a QP(or SVM)solver.It is also possible to change the algorithm so that no QP solver is needed.We will call the previous,QP-based,algorithm the max-margin method,and the new algorithm the projec-tion method.Brieﬂy,the projection method replaces step2of the algorithm with the following:-Set¯µ(i−1)=¯µ(i−2)+(µ(i−1)−¯µ(i−2))T(µE−¯µ(i−2))(µ(i−1)−¯µ(i−2))T(µ(i−1)−¯µ(i−2))(µ(i−1)−¯µ(i−2))(This computes the orthogonal projection ofµE ontothe line through¯µ(i−2)andµ(i−1).)-Set w(i)=µE−¯µ(i−1)-Set t(i)= µE−¯µ(i−1) 2In theﬁrst iteration,we also set w(1)=µE−µ(0)and ¯µ(0)=µ(0).The full justiﬁcation for this method is deferred to the full paper(Abbeel and Ng,2004),but in Sections4and5we will also give convergence results for it,and empirically compare it to the max-margin 6In k-dimensional space,any point that is a convex combination of a set of N points,with N>k+1,can be written as a convex combination of a subset of only k+1points of the original N points(Caratheodory’s Theorem,Rockafeller,1970).Applying this toµ=arg minµ∈Co{µ(π(i))}n i=0||µE−µ||2and{µ(π(i))}n i=0we ob-tain a set of k+1policies which is equally close to the expert’s feature expectations and thus have same perfor-mance guarantees.(Co denotes convex hull.)π(0))=µ(0)µ((2))µFigure2.Three iterations for projection algorithm. algorithm.An example showing three iterations of the projection method is shown in Figure2.4.Theoretical resultsMost of the results in the previous section were predi-cated on the assumption that the algorithm terminates with t≤ .If the algorithm sometimes does not ter-minate,or if it sometimes takes a very(perhaps ex-ponentially)large number of iterations to terminate, then it would not be useful.The following shows that this is not the case.Theorem1.Let an MDP\R,featuresφ:S→[0,1]k, and any >0be given.Then the apprenticeship learn-ing algorithm(both max-margin and projection ver-sions)will terminate with t(i)≤ after at mostn=O k22log k (14) iterations.The previous result(and all of Section3)had assumed thatµE was exactly known or calculated.In prac-tice,it has to be estimated from Monte Carlo samples (Eq.5).We can thus ask about the sample complex-ity of this algorithm;i.e.,how many trajectories m we must observe of the expert before we can guarantee we will approach its performance.Theorem2.Let an MDP\R,featuresφ:S→[0,1]k, and any >0,δ>0be given.Suppose the appren-ticeship learning algorithm(either max-margin or pro-jection version)is run using an estimateˆµE forµE obtained by m Monte Carlo samples.In order to en-sure that with probability at least1−δthe algorithm terminates after at most a number of iterations n given by Eq.(14),and outputs a policy˜πso that for any true reward R∗(s)=w∗Tφ(s)( w∗ 1≤1)we haveE[ ∞t=0γt R∗(s t)|˜π]≥E[ ∞t=0γt R∗(s t)|πE]− ,(15) it suﬃces thatm≥2k2log2k.The proofs of these theorems are in Appendix A.In the case where the true reward function R∗does not lie exactly in the span of the basis functionsφ,the algorithm still enjoys a graceful degradation of perfor-mance.Speciﬁcally,if R∗(s)=w∗·φ(s)+ε(s)forsome residual(error)termε(s),then our algorithmwill have performance that is worse than the expert’s by no more than O( ε ∞).5.Experiments5.1.GridworldIn ourﬁrst set of experiments,we used128by128 gridworlds with multiple/sparse rewards.The reward is not known to the algorithm,but we can sample tra-jectories from an expert’s(optimal)policy.The agent has four actions to try to move in each of the four com-pass directions,but with30%chance an action fails and results in a random move.The grid is divided into non-overlapping regions of16by16cells;we call these 16x16regions“macrocells.”A small number of the re-sulting64macrocells have positive rewards.For each value of i=1,...,64,there is one featureφi(s)indi-cating whether that state s is in macrocell i.Thus,the rewards may be written R∗=(w∗)Tφ.The weights w∗are generated randomly so as to give sparse rewards, which leads to fairly interesting/rich optimal policies.7 In the basic version,the algorithm is run using the 64-dimensional featuresφ.We also tried a version in which the algorithm knows exactly which macrocells have non-zero rewards(but not their values),so that the dimension ofφis reduced to contain only features corresponding to non-zero rewards.In Figure3,we compare the max-margin and projec-tion versions of the algorithm,whenµE is known ex-actly.We plot the margin t(i)(distance to expert’s pol-icy)vs.the number of iterations,using all64macro-cells as features.The expert’s policy is the optimal policy with respect to the given MDP.The two al-gorithms exhibited fairly similar rates of convergence, with the projection version doing slightly better. The second set of experiments illustrates the perfor-mance of the algorithm as we vary the number m of sampled expert trajectories used to estimateµE.The performance measure is the value of the best policy in the set output by the algorithm.We ran the al-gorithm once using all64features,and once using only the features that truly correspond to non-zero rewards.8We also report on the performance of three 7Details:We usedγ=0.99,so the expected horizon is of the order of the gridsize.The true reward function was generated as follows:For each macrocell i(i=1,...,64), with probability0.9the reward there is zero(w∗i=0), and with probability0.1a weight w∗i is sampled uniformly from[0,1].Finally,w∗is renormalized so that w∗ 1= 1.Instances with fewer than two non-zero entries in w∗are non-interesting and were discarded.The initial state distribution is uniform over all states.8Note that,as in the text,our apprenticeship learning algorithm assumes the ability to call a reinforcement learn-ing subroutine(in this case,an exact MDP solver using value iteration).In these experiments,we are interestedFigure the max-margin and projection versions of the algorithm on a 128x128grid.Euclidean distance to the expert’s feature expectations is plotted as a function of the number of it-erations.We rescaled the feature expectations by(1−γ) such that they are in[0,1]k.The plot shows averages over 40runs,with1s.e.errorbars.Figure tra-jectories from the expert.(Shown in color,where avail-able.)Averages over20instances are plotted,with1s.e. errorbars.Note the base-10logarithm scale on the x-axis. other simple algorithms.The“mimic the expert”al-gorithm picks whatever action the expert had taken if itﬁnds itself in a state in which it had previously observed the expert,and picks an action randomly oth-erwise.The“parameterized policy stochastic”uses a stochastic policy,with the probability of each action constant over each macrocell and set to the empiri-cal frequency observed for the expert in the macrocell. The“parameterized policy majority vote”algorithm takes deterministically the most frequently observed action in the macrocell.Results are shown in ing our algorithm,only a few sampled expert trajectories—far fewer than for the other methods—are needed to attain performance approaching that of the expert.(Note log scale on x-axis.)9Thus,by mainly in the question of how many times an expert must demonstrate a task before we learn to perform the same task.In particular,we do not rely on the expert’s demon-strations to learn the state transition probabilities.9The parameterized policies never reach the expert’s performance,because their policy class is not rich enough. Their restricted policy class is what makes them do betterFigure5.Screenshot of driving simulator. learning a compact representation of the reward func-tion,our algorithm signiﬁcantly outperforms the other methods.We also observe that when the algorithm is told in advance which features have non-zero weight in the true reward function,it is able to learn using fewer expert trajectories.5.2.Car driving simulationFor our second experiment,we implemented a car-driving simulation,and applied apprenticeship learn-ing to try to learn diﬀerent“driving styles.”A screen-shot of our simulator is shown in Figure5.We are driving on a highway at25m/s(56mph),which is faster than all the other cars.The MDP hasﬁve diﬀerent ac-tions,three of which cause the car to steer smoothly to one of the lanes,and two of which cause us to drive oﬀ(but parallel to)the road,on either the left or the right side.Because our speed isﬁxed,if we want to avoid hitting other cars it is sometimes necessary to drive oﬀthe road.The simulation runs at10Hz,and in the experiments that follow,the expert’s features were estimated from a single trajectory of1200samples(corresponding to 2minutes of driving time).There were features in-dicating what lane the car is currently in(including oﬀroad-left and oﬀroad-right,for a total ofﬁve fea-tures),and the distance of the closest car in the current lane.10Note that a distance of0from the nearest car implies a collision.When running the apprenticeship learning algorithm,the step in which reinforcement learning was required was implemented by solving a discretized version of the problem.In all of our exper-iments,the algorithm was run for30iterations,and a policy was selected by inspection(per the discussion in Section3).We wanted to demonstrate a variety of diﬀerent driv-ing styles(some corresponding to highly unsafe driv-ing)to see if the algorithm can mimic the same“style”in every instance.We consideredﬁve styles:1.Nice:The highest priority is to avoid collisions than the“mimic the expert”algorithm initially.10More precisely,we used the distance to the single clos-est car in its current lane,discretized to the nearest car length between-7to+2,for a total of10features.with other cars;we also prefer the right lane over the middle lane over the left lane,over driving oﬀ-road.2.Nasty:Hit as many other cars as possible.3.Right lane nice:Drive in the right lane,but gooﬀ-road to avoid hitting cars in the right lane. 4.Right lane nasty:Drive oﬀ-road on the right,butget back onto the road to hit cars in the right lane.5.Middle lane:Drive in the middle lane,ignoringall other cars(thus crashing into all other cars in the middle lane).After each style was demonstrated to the algorithm (by one of the authors driving in the simulator for2 minutes),apprenticeship learning was used to try to ﬁnd a policy that mimics demonstrated style.Videos of the demonstrations and of the resulting learned poli-cies are available at/~pabbeel/irl/In every instance,the algorithm was qualitatively able to mimic the demonstrated driving style.Since no “true”reward was ever speciﬁed or used in the experi-ments,we cannot report on the results of the algorithm according to R∗.However,Table1shows,for each of theﬁve driving styles,the feature expectations of the expert(as estimated from the2minute demonstra-tion),and the feature expectations of the learned con-troller for the more interesting features.Also shown are the weights w used to generate the policy shown. While our theory makes no guarantee about any set of weights w found,we note that the values there gener-ally make intuitive sense.For instance,in theﬁrst driving style,we see negative rewards for collisions and for driving oﬀroad,and larger positive rewards for driving in the right lane than for the other lanes.6.Discussion and ConclusionsWe assumed access to demonstrations by an expert that is trying to maximize a reward function express-ible as a linear combination of known features,and pre-sented an algorithm for apprenticeship learning.Our method is based on inverse reinforcement learning,ter-minates in a small number of iterations,and guaran-tees that the policy found will have performance com-parable to or better than that of the expert,on the expert’s unknown reward function.Our algorithm assumed the reward function is express-ible as a linear function of known features.If the set of features is suﬃciently rich,this assumption is fairly unrestrictive.(In the extreme case where there is a separate feature for each state-action pair,fully gen-eral reward functions can be learned.)However,it remains an important problem to develop methods for learning reward functions that may be non-linear func-tions of the features,and to incorporate automatic fea-。

LinguisticCyclesIntroduction

a. no/ne
early Old English
b. ne (na wiht/not) after 900, esp S
c. (ne) not
after 1350
d. not > -not/-n’t after 1400
20
Old English: (1) Men ne cunnon secgan to soðe ... hwa
so. (4) You maybe you've done it but have forgotten. (5) Me, I was flying economy, but the plane, … was
guzzling gas
8
Doubling and cliticization
(1) Me, I've tucking had it with the small place. (2) %Him, he .... (3) %Her, she shouldn’t do that (not
Then why me also I not-have not the right to fill-with-smoke the others some minutes in a bar
(4) tu
vas où
Colloquial French
2S go where
(5) nta tu vas travailler Arabic-French
Cognitive Economy (or UG) principles
help the learner, e.g:
Phrase > head (minimize structure)
Avoid too much movement

the bootstrap samples method

the bootstrap samples method
Bootstrap是一种统计学中的一种重要方法，也被称为自助法。

它的核心思想是通过重抽样技术模拟原始样本的分布，从而估计统计量的变异性并进行统计量区间估计。

具体步骤如下：
1. 从原始样本中有放回地抽取一定数量的样本，这个数量可以自己设定，通常与原始样本数量相同。

2. 根据抽出的样本计算待估计的统计量。

3. 重复上述过程多次（通常大于1000次），得到多个统计量。

4. 计算这些统计量的样本方差，以此估计统计量的方差。

Bootstrap的优势在于它不需要对总体分布做出任何假设，因此在样本量较小或总体分布未知的情况下也能够进行可靠的统计推断。

在实际应用中，Bootstrap方法被广泛用于参数估计和置信区间的计算、假设检验、模型选择、预测误差的估计等方面。

然而，Bootstrap方法也存在一些缺点。

例如在计算上可能会比较耗时，特别是在样本量较大时。

此外，Bootstrap方法对于样本分布的偏斜性和尾重性比较敏感，这可能会影响到Bootstrap估计的准确性。

并且，Bootstrap方法并不适用于所有的统计问题，特别是在样本量较小或总体分布明显偏离正态分布的情况下，其效果可能会受到限制。

Bundle adjustment

Bundle adjustmentA sparse matrix obtained when solving a modestly sized bundle adjustment problem. This is the sparsity pattern of a 992×992 normal equations (i.e. approximate Hessian) matrix.Black regions correspond to nonzero blocks.Given a set of images depicting a number of 3D points from differentviewpoints, bundle adjustment can be defined as the problem ofsimultaneously refining the 3D coordinates describing the scenegeometry as well as the parameters of the relative motion and theoptical characteristics of the camera(s) employed to acquire theimages, according to an optimality criterion involving thecorresponding image projections of all points.Bundle adjustment is almost always used as the last step of everyfeature-based 3D reconstruction algorithm. It amounts to anoptimization problem on the 3D structure and viewing parameters (i.e.,camera pose and possibly intrinsic calibration and radial distortion), toobtain a reconstruction which is optimal under certain assumptionsregarding the noise pertaining to the observed image features: If the image error is zero-mean Gaussian, then bundle adjustment is the Maximum Likelihood Estimator. Its name refers to the "bundles" of light rays originating from each 3D feature and converging on eachcamera's optical center, which are adjusted optimally with respect to both the structure and viewing parameters.Bundle adjustment was originally conceived in the field of photogrammetry during 1950s and has increasingly been used by computer vision researchers during recent years.Bundle adjustment boils down to minimizing the reprojection error between the image locations of observed and predicted image points, which is expressed as the sum of squares of a large number of nonlinear, real-valued functions. Thus, the minimization is achieved using nonlinear least-squares algorithms. Of these,Levenberg –Marquardt has proven to be one of the most successful due to its ease of implementation and its use of an effective damping strategy that lends it the ability to converge quickly from a wide range of initial guesses. By iteratively linearizing the function to be minimized in the neighborhood of the current estimate, the Levenberg –Marquardt algorithm involves the solution of linear systems known as the normal equations. When solving the minimization problems arising in the framework of bundle adjustment, the normal equations have a sparse block structure owing to the lack of interaction among parameters for different 3D points and cameras. This can be exploited to gain tremendous computational benefits by employing a sparse variant of the Levenberg –Marquardt algorithm which explicitly takes advantage of the normal equations zeros pattern, avoiding storing and operating on zero elements.Mathematical definitionBundle adjustment amounts to jointly refining a set of initial camera and structure parameter estimates for finding the set of parameters that most accurately predict the locations of the observed points in the set of available images.More formally, assume that 3D points are seen in m views and let be the projection of the i th point on image . Let denote the binary variables that equal 1 if point is visible in image and 0 otherwise. Assume alsothat each camera is parameterized by a vectorand each 3D point by a vector . Bundle adjustmentminimizes the total reprojection error with respect to all 3D point and camera parameters, specificallywhere is the predicted projection of point i on image j and denotes the Euclidean distance between the image points represented by vectors and . Clearly, bundle adjustment is by definition tolerant to missing image projections and minimizes a physically meaningful criterion.Software•sba [1]: A Generic Sparse Bundle Adjustment C/C++ Package Based on the Levenberg–Marquardt Algorithm (C, Matlab)•ssba [2]: Simple Sparse Bundle Adjustment package based on the Levenberg–Marquardt Algorithm (C) with LGPL license.References• B. Triggs; P. McLauchlan and R. Hartley and A. Fitzgibbon (1999). "Bundle Adjustment — A Modern Synthesis". ICCV '99: Proceedings of the International Workshop on Vision Algorithms. Springer-Verlag.pp. 298–372. doi:10.1007/3-540-44480-7_21. ISBN 3-540-67973-1.•M.I.A. Lourakis and A.A. Argyros (2009). "SBA: A Software Package for Generic Sparse Bundle Adjustment".ACM Transactions on Mathematical Software (ACM) 36 (1): 1–30. doi:10.1145/1486525.1486527.•R.I. Hartley and A. Zisserman (2004). Multiple View Geometry in computer vision (2nd ed.). Cambridge University Press. ISBN 0521540518.External links• B. Triggs, P. McLauchlan, R. Hartley and A. Fitzgibbon, Bundle Adjustment — A Modern Synthesis, Vision Algorithms: Theory and Practice, 1999 [3].• A. Zisserman. Bundle adjustment [4]. CV Online.References[1]http://www.ics.forth.gr/~lourakis/sba/[2]http://www.inf.ethz.ch/personal/chzach/opensource.html[3]http://lear.inrialpes.fr/pubs/2000/TMHF00/Triggs-va99.pdf[4]/rbf/CVonline/LOCAL_COPIES/ZISSERMAN/bundle/bundle.htmlArticle Sources and Contributors3 Article Sources and ContributorsBundle adjustment Source: /w/index.php?oldid=400293908 Contributors: Andy M. Wang, BenFrantzDale, Daviddoria, Lourakis, Michael Hardy, Physchim62, Thrapper,12 anonymous editsImage Sources, Licenses and ContributorsImage:Bundle adjustment sparse matrix.png Source: /w/index.php?title=File:Bundle_adjustment_sparse_matrix.png License: Creative Commons Attribution 3.0Contributors: Lourakis, 1 anonymous editsLicenseCreative Commons Attribution-Share Alike 3.0 Unported/licenses/by-sa/3.0/。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

a rX iv:mat h /7237v1[mat h.CO]13Fe b27MINIMAL PERCOLATING SETS IN BOOTSTRAP PERCOLATION ROBERT MORRIS Abstract.In standard bootstrap percolation,a subset A of the grid [n ]2is initially infected .A new site is then infected if at least two of its neighbours are infected,and an infected site stays infected forever.The set A is said to percolate if eventually the entire grid is infected.A percolating set is said to be minimal if none of its subsets percolate.Answering a question of Bollob´a s,we show that there exists a minimal percolating set of size 4n 2/33+o (n 2),but there does not exist one larger than (n +2)2/6.1.Introduction Consider the following deterministic process on a (ﬁnite,connected)graph G .Given an initial set of ‘infected’sites,A ⊂V (G ),a vertex becomes infected if at least r of its neighbours are already infected,and infected sites remain infected forever.This process is known as r -neighbour bootstrap percolation on G .If eventually the entire vertex set becomes infected,we say that the set A percolates on G .For a given graph G ,we would like to know which sets percolate.Bootstrap percolation is an example of a cellular automaton (see [13]and [16]),since its update rule is homogeneous (i.e.,it is the same for every vertex)and local (it depends only on the states of the neigh-bours of a vertex).As such,it is related to interacting systems of particles;indeed,the bootstrap process was introduced (on the inﬁnite grid)by Chalupa,Leith and Riech [9]in 1979in order to study ferro-magnetism.The ﬁrst rigorous results for the inﬁnite case were proved by van Enter [10]and Schonmann [14],[15],and for the ﬁnite case by Aizenman and Lebowitz [2].For more on the various physical applica-tions of bootstrap percolation,we refer the reader to the survey article of Adler and Lev [1],and the references therein.In this paper we shall study 2-neighbour bootstrap percolation on the m ×n grid,i.e.,we shall set r =2,and our base graph,G (m,n ),will have vertex set V =[m ]×[n ],and edge set E ={uv :|u 1−v 1|+|u 2−2ROBERT MORRISv2|=1}.The question of which sets percolate on G(n,n)when the elements of A are chosen independently at random,with someﬁxed probability p=p(n),has been extensively studied[2],[5],[11].In particular,Holroyd[11]proved that if|A|=λn2log nthen as n→∞,P(A percolates)→ 0ifλ<π2/181ifλ>π2/18.For results on other graphs,and for other values of r,see[3],[4],[?], [7],[8]and[12],for example.Another interesting question about percolating sets is,how small can they be?For the graph G(n,n),this is just an exercise(which we leave to the reader!),but for many other graphs(and values of r)it is open, and seems to be a fairly challenging problem.Bollob´a s asked the following,somewhat diﬀerent question about per-colating sets.Call a percolating set A minimal if A\{v}does not percolate for each vertex v∈A.What can we say about the minimal percolating sets(MinPS)of G(n,n);in particular,how large can such a set be?Let us deﬁneE(m,n)=max{|A|:A⊂[m]×[n]is a MinPS of G(m,n)}, and write E(n)=E(n,n).Thus our problem is to determine E(m,n) for every m,n∈N.It is easy to construct a MinPS with about2(m+n)/3elements. For example(assuming for simplicity that m,n≡0(mod3)),take A={(k,1):k≡0,2(mod3)}∪{(1,ℓ):ℓ≡0,2(mod3)}.It is easy to see that A percolates,and that if x∈A,then A\{x}does not percolate.For example,if x=(3,1)then the3rd and4th columns of V=[m]×[n]are empty.However,it is non-trivial toﬁnd a MinPS with more than2(m+n)/3elements,and one is easily tempted to suspect that in fact E(m,n)=⌊2(m+n)/3⌋.(The interested reader is encouraged to stop at this point and try to construct a minimal percolating set with more than this many elements!)As it turns out, however,the correct answer is rather a long way from this.In fact,there exist MinPS which are(essentially)as large as they possibly could be. The following theorem is the main result of this paper.MINIMAL PERCOLATING SETS IN BOOTSTRAP PERCOLATION3Theorem1.For every2 m,n∈N,we have4mnm−m3/2 E(m,n)(m+2)(n+2)33−9n3/2 E(n) (n+2)2m 66and33√n2exists.The rest of the paper is organised as follows.In Section2we shall deﬁne corner-avoiding minimal percolating sets,which will be instru-mental in the proofs of Theorems1and2,and use them to prove various inequalities satisﬁed by E(m,n).In Section3we shall describe our construction which gives the lower bound;in Section4we shall prove the upper bound;in Section5we shall prove Theorem2;and in Section6we shall discuss possible further work,including the equiva-lent problem on[n]d.2.Corner-avoiding setsWe being with some simple deﬁnitions.Let m,n∈N,and let V= [m]×[n].Given any two vertices x,y∈V,deﬁned(x,y)=|x1−y1|+|x2−y2|,and for sets X,Y⊂V,deﬁned(X,Y)=min{d(x,y):x∈X,y∈Y}.A rectangle is a set[(a,b),(c,d)]={(x,y):a x c,b y d},where a,b,c,d∈N.For any rectangle R=[(a,b),(c,d)],deﬁne dim(R)=(width(R),height(R))=(c−a+1,d−b+1). Now,given a set X⊂V,write[X]for the set of points which are eventually infected if the initial set is X.Observe that in G(m,n),[X] is always a union of rectangles.If Y⊂[X]then we shall say that X spans Y,and if moreover Y⊂[X∩Y],then we say that X internally4ROBERT MORRISspans Y.Note that A percolates on G(m,n)if and only if A spans [m]×[n].The following easy lemma is from[3].Lemma3.A set A⊂[m]×[n]does not percolate on G(m,n)if,and only if,there exists a set R1,...,R t of rectangles such that A⊂ i R i, and d(R i,R j) 3for every i=j.Moreover,[A]⊂ i R i.We now make an important deﬁnition.Let m,n∈N and V=[m]×[n].The top-left corner of V is the rectangle J L=[(1,n−1),(2,n)]andthe bottom-right corner of V is the rectangle J R=[(m−1,1),(m,2)]. Call a minimal percolating set A⊂V corner-avoiding if wheneverv∈A,we have[A\{v}]∩(J L∪J R)=∅,i.e.,if the initially infected sites are a(proper)subset of A,then thetop-left and bottom-right corners remain uninfected.DeﬁneE c(m,n)=max{|A|:A⊂[m]×[n]is a corner-avoiding MinPS of G(m,n)} if such sets exist,and let E c(m,n)=0otherwise.As before,writeE c(n)=E c(n,n).E(m,n)and E c(m,n)are related in the following way.Lemma4.If m,n∈N,with n 4,thenE(m,n) E c(m+16,n+8) E(m+16,n+8).The proof of Lemma4makes use of the following simple structures. Given a set A⊂[m]×[n],and integers k,ℓ∈N,deﬁneA+(k,ℓ)={(i,j)∈N2:(i−k,j−ℓ)∈A}.Now,let P be the pair of points{(1,1),(1,3)},and for each k∈N letL(k)=k−1i=0 P+(0,3i) .Furthermore,for each a,b∈N letL(k;a,b)=L(k)+(a−1,b−1).Observe that[L(k;a,b)]=[(a,b),(a,b+3k−1)],and that L(k;a,b)is a minimal spanning set for[L(k;a,b)].We are now ready to prove Lemma4.Proof of Lemma4.Let m,n∈N,with n 4,and let A be a minimal percolating set of G(m,n)with|A|=E(m,n).Let m′=m+16and n′=n+8,let V=[m′]×[n′],let N=⌊n+2MINIMAL PERCOLATING SETS IN BOOTSTRAP PERCOLATION 5Now,let B ′be the set obtained by rotating B through 180◦,and placing the top-right corner at the point (m ′,n ′),i.e.,B ′={(x,y ):(m ′−x +1,n ′−y +1)∈B }.Finally,letC =B ∪(A +(8,4))∪B ′,so C ⊂V (see Figure 1).We claim that C is a corner-avoiding minimal percolating set of G (m ′,n ′)....................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................[A +(8,4)].............Figure 1:The set CClearly C percolates (note that we needed n 4here,so that 3N >3),and if v ∈C \(A +(8,4))then C \{v }does not percolate,and avoids the (top-left and bottom-right)corners.So let v ∈A ,and let R 1,...,R t be rectangles satisfying (A \{v })+(8,4)⊂ k R k ,and d (R i ,R j ) 3for each 1 i <j t .Such rectangles exist by Lemma 3,since A is minimal.Let u =(7,5)and w =(m +10,n +4).If (9,5)∈R i for some i ∈[t ]then let R ′i =[R i ∪{u }],and note that d (R ′i ,R k ) 3for each i =k ∈[t ].Similarly,if (m +8,n +4)∈R j for some j ∈[t ]then let R ′j =[R j ∪{w }],and note that d (R ′j ,R k ) 3for each j =k ∈[t ].Note also that i =j ,since A \{v }does not percolate.Now,if (7,n +4)∈R ′i ,then let R ′′i =[[(1,1),(5,n +4)]∪R ′i ],and note that [(9,5),(9,n +4)]⊂R i .Thus d (R ′′i ,R k ) 3for every i =k ∈[t ],and d (R ′′i ,R ′j ) 3.Similarly,if (m +10,5)∈R ′j ,then6ROBERT MORRISdeﬁne R′′j=[[(m+10,5),(m′,n′)]∪R′j],and note that d(R′′j,R k) 3 for every j=k∈[t],and also d(R′i,R′′j) 3and d(R′′i,R′′j) 3. Thus,in each case we canﬁnd rectangles T1,...,T t′in[m′]×[n′] such that C⊂ T k,d(T i,T j) 3for every1 i<j t′,and each T i avoids the corners of V.Thus C is a corner-avoiding minimal percolating set,as claimed.By Lemma4,corner-avoiding minimal percolating sets exist in G(m,n) if m 17and n 12.The next lemma explains our interest in corner-avoiding minimal percolating sets.Lemma5.Let17 m,m′,n∈N.ThenE c(m+m′+3,n+2) E c(m,n)+E c(m′,n)+2.Proof.Let B⊂[m]×[n]and C⊂[m′]×[n]be corner-avoiding MinPS, with|B|=E c(m,n)and|C|=E c(m′,n).Note that B and C exist by Lemma4,since m,m′,n 17.Now,letC′=C+(m+3,2)⊂[m+m′+3]×[n+2],and letA=B∪{(m+1,1),(m+3,n+2)}∪C′.Then A is a corner-avoiding MinPS in[m+m′+3]×[n+2],and |A|=E(m,n)+E(m′,n)+2. In fact we can generalize this construction in the following way. Lemma6.Let k,m,n∈N.ThenE c(km+3(k−1),n+2(k−1)) kE c(m,n).Proof.If E c(m,n)=0then the result is trivial,so let us assume that E c(m,n) 1,and let B be a corner-avoiding minimal percolating set of G(m,n)with|B|=E c(m,n).Now,let m′=km+3(k−1), n′=n+2(k−1),V=[m′]×[n′],and for each i∈[k],letB i=B+((i−1)(m+3),2(i−1))⊂V.Also,let C={(m+1,1),(m+3,n+2)},and for each i∈[k−1],letC i=C+((i−1)(m+3),2(i−1)).We claim thatk i=1B i∪k−1 i=1C iA=is a corner-avoiding minimal percolating set of G(m′,n′)(see Figure2).MINIMAL PERCOLATING SETS IN BOOTSTRAP PERCOLATION7...................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................[B 1][B 2][B k ]Figure 2:The set AIt is clear that A percolates,and that if v ∈ C i then A \{v }does not percolate,since there are two adjacent empty columns.So let v ∈B i for some i ∈[k ].By symmetry,it is enough to consider the cases i =1and i =2.Suppose v ∈B 1,let u =(m +1,1),and let R 1,...,R t be the rect-angles given by Lemma 3.Since B is corner-avoiding,d (R i ,u ) 3for each i ∈[t ].Let R t +1={u },and let R t +2=[(m +3,3),(m ′,n ′)].Then R 1,...,R t +2satisfy A \{v }⊂ R i ,and d (R i ,R j ) 3for each 1 i <j t +2,so by Lemma 3,A \{v }does not percolate,and avoids the corners.The proof for v ∈B 2is very similar.Let u =(m +3,n +2),w =(2m +4,3),and let R 1,...,R t be the rectangles given by Lemma 3.LetR t +1=[(1,1),(m +1,n )],R t +2={u },R t +3={w },and R t +4=[(2m +6,5),(m ′,n ′)].Now R 1,...,R t +4satisfy A \{v }⊂ R i ,and d (R i ,R j ) 3for each 1 i <j t +4,so we are again done by Lemma 3. We shall now deduce from Lemmas 4and 5the other recurrences that we shall use to prove Theorems 1and 2.First however we need to make a couple of observations.For each a ∈R ,let g a (x )=2x +a .Observation 7.g t a (x )=2t (x +a )−a for every t ∈N .8ROBERT MORRISThe next lemma,which we shall use several times,says that E(m,n) is increasing in both m and n.Lemma8.If k m andℓ n,then E(k,ℓ) E(m,n).Proof.It is clearly enough to prove the lemma in the case that n=ℓand m=k+1.So let V=[m−1]×[n],and let A⊂V be a MinPS in G(m−1,n).Since A percolates,(m−1,a)∈A for some a∈[n];we claim that either B=A∪{(m,a)}or C=A∪{(m,a)}\{(m−1,a)} is a MinPS for G(m,n).Supposeﬁrst that C\{u}percolates for some u∈C.It is easy to see that A\{u}also percolates(in G(m−1,n)).Since(m,a)is the only element in column m,u=(m,a),so this is a contradiction.So suppose that C does not percolate,but B\{v}does percolate. It is again easy to see that A\{v}also percolates,and that v/∈{(m−1,a),(m,a)},since C doesn’t percolate.This contradiction completes the proof. We can now easily deduce the desired recurrence relations for E(m,n).Lemma9.Let m,n,t∈N,with m,n 17.Then(a)E c(2t(m+19)−3),n+2t+8) 2t E(m,n),and(b)E(2t(n+37)−37) 22t E(n).Proof.We apply Lemma5t times to E(m,n).To be precise,Lemma5 with m=m′gives E c(2m+3,n+2) 2E c(m,n),and it follows that E c(g t3(m),n+2t) 2t E c(m,n).But nowE c(2t(m+19)−3,n+2t+8)=E c(g t3(m+16),n+2t+8)2t E c(m+16,n+8) 2t E(m,n)by Lemma4and Observation7.This proves part(a).To prove part(b),note that since E c(2m+3,n+2) 2E c(m,n)and E c(m+2,2n+3) 2E c(m,n),we have E c(2m+5,2n+7) 4E c(m,n), for every m,n∈N.ThusE(2n+37) E c(2n+37,2n+23) 4E c(n+16,n+8) 4E(n).So,by Observation7,we haveE(2t(n+37)−37)=E(g t37(n)) 4t E(n),as required.MINIMAL PERCOLATING SETS IN BOOTSTRAP PERCOLATION 93.A large minimal setWe shall now use the results of the previous section to construct a minimal percolating set in G (m,n )of size (4/33+o (1))mn .The construction will have three stages.First,we use the ideas from the proof of Lemma 4toﬁndacorner-avoiding MinPS in G (8,3k +2)for all k ∈N .Next,using Lemma 6,we put about√11 n −2 m +33 411−7n .Proof.Let M = m +33 .The result is trivial if M (N +1) 0,so assume that M 1and N 0.Applying Lemma 6with m =8,n =3N +2and k =M ,we obtainE c (11M −3,3N +2M ) ME c (8,3N +2) 4M (N +1),by Lemma 10.But m 11M −3and n 3N +2N ,so by Lemma 8we haveE (m,n ) E c (11M −3,3N +2M ).The ﬁnal inequality is trivial.We are now ready to prove the lower bound in Theorem 1.Proof of the lower bound in Theorem 1.Let m,n ∈N .We may as-sume that n 8√2log 2m ⌉,x =⌊(m +3)/2t ⌋,M =x −19and N =n −2t −8.By10ROBERT MORRISLemma 11,we know thatE (M,N ) S (M,N )411−7N>411−26n .Now,by Lemma 9,and noting that 2t x >m −2t ,x <√m ,we obtainE (m,n ) E (2t (M +19)−3,N +2t +8)2t E (M,N )>2t411−26n >4n (m −2t )363−7n√33−8n √3(b )E (m,2)= 2(m +2)3Proof.The lower bounds are easy,so we shall only prove the upper bounds.In each case,let A be a minimal percolating set.To prove part (a ),simply note that A ⊂[m ]×[1]can contain at most two out of three consecutive points.For part (b ),observe that if A ⊂[m ]×[2]percolates,there must exist s,t ∈[m ]such that (s,1),(t,2)∈A and |s −t | 1.Indeed,it is easy to see that if no such s and t exist,then [{(k,1)∈A }]and[{(k,2)∈A }]are at distance at least 3.There are thus two cases.IfMINIMAL PERCOLATING SETS IN BOOTSTRAP PERCOLATION11 s=t then(i,j)/∈A for i∈{s−1,s+1},j∈{1,2},and A can contain at most two points from any three consecutive columns.Thus|A| 2(s−1)3 2m+43 +4+ 2(m−s−1)3 .The reader can easily check that when s 3or s m−1,the calcu-lation is exactly the same.Part(c)requires a little more work,and will be proved by induction on m.Observe that the result follows by parts(a)and(b)if m 2, and that E(3,3)=E(4,3)=4.So let m 5,and assume that the result holds for all smaller m.Supposeﬁrst that there exists an internally spanned rectangle R, with dim(R)=(k,3),which does not contain either the(m−1)st or the m th column of V=[m]×[3].Then either[m−3]×[3]or[m−2]×[3] must be internally spanned.In the former case,we have|A| E(m−3,3)+2 2m3 ,while in the latter case we have|A| E(m−2,3)+1 2(m+1)3 .So assume that no such rectangle R exists(and similarly for the1st and2nd columns of V),and observe that there must therefore exist some internally spanned rectangle T with dim(T)=(2,2),since otherwise A would not percolate,as before.Without loss of generality,we may assume(since m 5)that T does not intersect either the(m−1)st or the m th column of V.Now,by allowing T to grow one block at a time,weﬁnd that either [m−3]×[2]is internally spanned,or[m−2]×[2]is internally spanned, or there exists an internally spanned rectangle T′,with dim(T′)=(ℓ,2) for someℓ∈[m−4],such that d(A\T′,T′) 3.If[m−2]×[2]is internally spanned,then|A| E(m−2,2)+2 2m3 .12ROBERT MORRISAlso,if[m−3]×[2]is internally spanned but[m−2]×[2]is not,then |A| E(m−3,2)+2 2(m−1)3 , since if|A∩[(m−1,1),(m,3)]| 3,then[(m−1,1),(m,3)]is internally spanned,which contradicts our earlier assumption.So,without loss of generality,T′=[ℓ]×[2]is internally spanned, and d(A\T′,T′) 3,for someℓ∈[m−4].But then the rectangle [(ℓ+2,2),(m,3)]must be internally spanned,since A percolates and there is no internally spanned k×3rectangle R in V.Thus |A| E(ℓ,2)+E(m−ℓ−1,2)2(ℓ+2)3 2(m+3).6Let<R be the following partial order on rectangles in[m]×[n].First, given a,c∈[m]and b,d∈[n],let(a,b)<R(c,d)if min{m−a,n−b}> min{m−c,n−d},or min{m−a,n−b}=min{m−c,n−d}and max{m−a,n−b}>max{m−c,n−d}.Now,given rectangles S and T,let S<R T if and only if dim(S)<R dim(T).Observation14.If(p,q) R(k,ℓ),then kℓ+p(n−ℓ)+q(m−k) mn. Proof.Note thatkℓ+p(n−ℓ)+q(m−k)=mn+(m−k)(n−ℓ)−(m−p)(n−ℓ)−(m−k)(n−q). Now,if p k then(m−k)(n−ℓ) (m−p)(n−ℓ),while if p>k, then q<ℓ,and so(m−k)(n−ℓ) (m−k)(n−q).In either case,the result follows. Proof of the upper bound in Theorem1.If2 min{m,n} 3then the result follows by Corollary13,and note that the result also holds if m=n=1.So let m,n∈N,with m,n 4,let A be a minimal percolating set in V=[m]×[n],and assume that if[p]×[q] V and p,q 2,then E(p,q) (p+2)(q+2)/6.We shall show that |A| (m+2)(n+2)/6.In order to aid the reader’s understanding, we shall let a=1/6,b=1/3and c=2/3,so that(m+2)(n+2)/6= amn+b(m+n)+c.MINIMAL PERCOLATING SETS IN BOOTSTRAP PERCOLATION13 Let S be a maximal(in the order<R)internally spanned rectangle in V,other than V itself,and let dim(S)=(k,ℓ).We shall distinguish several cases.Case1:Either k=m orℓ=n.Suppose that k=m.Since A percolates,there cannot be two con-secutive empty rows,so since S is maximal,we must haveℓ=n−2or n−1,which means that|A\S|=1(since A is minimal).Hence,by the induction hypothesis,|A| E(m,n−1)+1am(n−1)+b(m+n−1)+c+1=amn+b(m+n)+c−(am+b−1)amn+b(m+n)+c,since m 4,so am+b 1.The proof ifℓ=n is identical.Assume from now on that m−k,n−ℓ 1,and let B=A\S. Since S is maximal and A is minimal,it follows that either|B|=1,or d(S,B) 3.Now let T=[B],and note that T is a rectangle,since some rectangle T1⊂T must have d(S,T1) 2,and thus[S∪T1]=V by maximality of S.Since A is minimal,we must have B\T1=∅,andso T=T1.Since S is maximal,it contains some corner of V,so without loss of generality we can assume that S=[k]×[ℓ].Say that S and T overlap rows if they both contain an element of some row of V,and say that they overlap columns if they both contain an element of some column. Case2:S and T neither overlap rows,nor overlap columns.It is clear that T⊂[k+1,m]×[ℓ+1,n],so|A| E(k,ℓ)+E(m−k,n−ℓ)a(kℓ+(m−k)(n−ℓ))+b(k+ℓ+(m−k)+(n−ℓ))+2c=amn+b(m+n)+c−(a(k(n−ℓ)+ℓ(m−k))−c)amn+b(m+n)+csince k,m−k 1,so k(n−ℓ)+ℓ(m−k) n 4,and4a−c=0.14ROBERT MORRISSo assume,without loss of generality,that S and T overlap columns. Thus|B| 2,so d(S,B) 3and n−ℓ 3.Furthermore,both of the dimensions of T must be at least two,and B is a MinPS for T.Thus, (for exactly the same reasons that S and T must exist),there exist disjointly internally spanned rectangles P and Q,such that P,Q=T, but[P∪Q]=T.Choose P and Q so that min{|P|,|Q|}is minimal subject to these conditions.Moreover,given min{|P|,|Q|},choose P and Q to each have both dimensions at least two if possible. Suppose that d(S,P) 2.Then[S∪P]=V(since S is maximal), and so B⊂P(since A is minimal),so P=T.But we chose P=T, so in fact we must have d(S,P) 3,and similarly d(S,Q) 3.Let dim(P)=(p,s)and dim(Q)=(t,q).We claim that either min{p,q,s,t} 2,or|Q|=1and both di-mensions of P are at least two(or vice-versa),or|B| 3.Indeed, suppose that q=1but t>1,say.We have d(P,Q) 2;supposeﬁrst that some internally spanned proper sub-rectangle Q′ Q also satisﬁes d(P,Q′) 2.Since q=1,we may choose Q′so that|B∩(Q\Q′)|=1. But P′=[P∪Q′]=T,since A is minimal,so we could have chosen P and Q with|Q|=1and both dimensions of P at least two.This contradicts our choice of P and Q.So suppose next that no internally spanned sub-rectangle Q′ Q satisﬁes d(P,Q′) 2.The only possibility is that q=1and t=3. Now either|P|=1,in which case|B|=3,or p=1and there exists a internally spanned proper sub-rectangle P′ P with d(P′,Q) 2, in which case we may repeat the argument of the previous paragraph. Thus,we conclude that either min{p,q,s,t} 2,or(without loss of generality)|Q|=1and both dimensions of P are at least two,or |B| 3,as claimed.Case3:|B| 3.Recall that m−k 1and n−ℓ 3,so|A| E(m−1,n−3)+3a(m−1)(n−3)+b(m+n−4)+c+3=amn+b(m+n)+c−(a(3m+n−3)+4b−3)amn+b(m+n)+c,since m,n 4,so a(3m+n−3)+4b−3 13a+4b−3>0.MINIMAL PERCOLATING SETS IN BOOTSTRAP PERCOLATION15So assume from now on that|B| 4,and so by the comments above, both dimensions of P are at least two,and either|Q|=1or both di-mensions of Q are least two.Case4:|Q|=1,and S and P overlap columns.Since d(S,P) 3,S and P cannot overlap both rows and columns,so they don’t overlap rows.It follows that|A| E(k,ℓ)+E(p,n−ℓ−2)+1a(kℓ+pn−pℓ−2p))+b(k+p+n−2)+2c+1=a(pn+ℓ(k−p))+b(k+n)−p(2a−b)−(2b−c)+c+1amn+b(m+n)+cif an(m−p)+b(m−k) aℓ(k−p)+1,since2a=b and2b=c.Note that m−p 1(since S was maximal in the order<R).But now ifp k then we can replace aℓ(k−p)by0,and if p<k then we can replace it by an(k−p),and in either case the inequality follows from the fact that an+b 1.Case5:|Q|=1,and S and P do not overlap columns.Since S and T overlap columns,S and Q must overlap columns. But|Q|=1and d(P,Q) 2,so S and P cannot overlap rows.Alsod(S,P) 3,so(k+1,ℓ+1)∈P.Thus either|A| E(k,ℓ)+E(m−k−1,n−ℓ)+1a(kℓ+(m−k−1)(n−ℓ))+b(m+n−1)+2c+1=amn+b(m+n)+c−(a((k+1)(n−ℓ)+ℓ(m−k))+b−c−1), or|A| E(k,ℓ)+E(m−k,n−ℓ−1)+1amn+b(m+n)+c−(a(k(n−ℓ)+(ℓ+1)(m−k))+b−c−1). Now,n−ℓ 3,so(k+1)(n−ℓ)+ℓ(m−k) m+2k+3 9,andk(n−ℓ)+(ℓ+1)(m−k) 2m+k 9.But9a+b−c−1 0,so in either case we have|A| amn+b(m+n)+c.So assume from now on that min{p,q,s,t} 2.Note that n−ℓ 4, and assume that S and P overlap columns.Note that height(T)q+1,since otherwise we could have chosen P and Q with|P|=1.16ROBERT MORRISSimilarly,width(T) p+1.Case6:S and T overlap both rows and columns.Since d(S,P) 3and we assumed that S and P overlap columns,it follows that S and P do not overlap rows,and that S and Q overlap rows but not columns.Therefore,|A| E(k,ℓ)+E(p,n−ℓ−2)+E(q,m−k−2)a(kℓ+p(n−ℓ−2)+q(m−k−2))+b(m+n+p+q−4)+3c amn+b(m+n)−(p+q)(2a−b)−2(2b−c)+c=amn+b(m+n)+csince2a=b and2b=c,and by Observation14.There is one remaining case to consider.Case7:S and T overlap columns but not rows.Note that q n−ℓ−1,since height(T) q+1,and p m−2, since width(T) p+1.We have|A| E(k,ℓ)+E(p,n−ℓ−2)+E(q,m−k)a(kℓ+p(n−ℓ−2)+q(m−k))+b(m+n+p+q−2)+3camn+b(m+n)−(p+q)(2a−b)−2(2b−c)+c=amn+b(m+n)+cifM=(m−p)(n−ℓ)+(m−k)(n−q)−(m−k)(n−ℓ) 2(q+b/a), since kℓ+p(n−ℓ)+q(m−k)=mn−M.Now,if k p+2,thenM (k−p)(n−ℓ)+(ℓ+1)(m−k) 2(n−ℓ)+2 2(q+b/a).But if k p+1 width(T),thenℓ n−ℓ−1(else S<R T),soM (m−p)(n−ℓ)+(m−k)(2ℓ+1−n) 3(n−ℓ) 2(q+b/a)if p m−3.However,if p=m−2,then k=m−1(sinceℓ n−4, so if k m−2then S<R T),so T contains no element of rowℓ+1, and so q height(T)−1 n−ℓ−2.ThusM (m−p)(n−ℓ)+(m−k)(2ℓ+2−n) 2(n−ℓ) 2(q+b/a), and we are done.MINIMAL PERCOLATING SETS IN BOOTSTRAP PERCOLATION 175.Proof of Theorem 2In this section we shall prove that the sequence E (n )/n 2converges.The proof uses only the basic techniques from Section 2.Lemma 15.Let n,N ∈N and ε>0,with n (20/ε)2and N n 3/2 213.ThenE (N )n 2.Proof.With foresight,choose t ∈N so thatNn 3/2,and letx = N +37n −3,and recall that 2t N n ,sokn +2k 2+21k x −3n +2n +21k x,since √18ROBERT MORRISand22t k2=22t x2t n−38n2−10·2t Nn2,since2t 2Nn20N2>(1−ε/2)E(n)n2?We can,of course,ask exactly the same questions for graphs other than G(m,n),and for diﬀerent values of r.Thus,for each graph G, and each r∈N,deﬁneE(G,r)=max{|A|:A⊂V(G)is a minimal percolating set of Gin r-neighbour bootstrap percolation}. Having obtained results for G(n,n),it is natural to ask whether similar bounds hold for the n×n torus,T(n,n).The following theorem says that they do.MINIMAL PERCOLATING SETS IN BOOTSTRAP PERCOLATION19 Theorem16.E(T(n,n),2)=E(G(n,n),2)+o(n2).Sketch of proof.We begin by showing that for every m,n∈N,E(T(m+4,n+4),2) E c(m,n)+4.Let A⊂[m]×[n]be a corner-avoiding MinPS of G(m,n)with|A|= E c(m,n).Let m′=m+4,n′=n+4,and letB=(A+(2,2))∪{(2,n+2),(4,n+4),(m+2,2),(m+4,4)}⊂[m′]×[n′]. It is easy to see that B is a minimal percolating set for T(m′,n′),and so E(T(m+4,n+4),2) E c(m,n)+4,as claimed.Now,let A be a MinPS for T(n,n),and build up the infected sites as follows.Let A0=A,and form A i+1from A i by choosing two entirely infected rectangles S and T in A i,with d(S,T) 2,and infecting all sites of[S∪T](see[3]for a more detailed description of this process). Let t∈N be deﬁned as theﬁrst time at which A t contains a rectangle which loops all the way around the torus,and consider A t−1.It consists of two entirely spanned rectangles,X and Y,say,each of which is a minimally internally spanned subset of G(n,n),and at most n extra infected sites.Thus,by choosing n suﬃciently large,we canﬁnd an arbitrarily large rectangle inside T(n,n)which is minimally internally spanned(as a subset of G(n,n)),and has about the same(or greater) density as A.Now,by an argument similar to the proof of Theorem2, it follows that T(n,n) E(n)+o(n2). For other graphs however,the problem is likely to be more diﬃcult. Problem1.Give bounds on E(G,r)for other graphs G and integers r.In particular,determine(a)E([n]d,2),where[n]d is the d-dimensional grid,(b)E(G n,p,r),where G n,p is the random graph,and(c)E(G n,d,r),where G n,d is the random d-regular graph.The method of this paper allows us to make a start on part(a). Theorem17.For every n,d∈N,we haveE([n]d,2) C(d)n d,where C(d)>0is independent of n.Sketch of proof.Divide[n]d into n hyperplanes,each of co-dimension 1.Take a(corner-avoiding)construction for d−1in every fourth hyperplane,and put a single point in opposite corners of everyﬁrst and third(mod4)hyperplane,along the lines of Lemma6.This set is now a corner-avoiding MinPS of[n]d.。