text:
InfoQ Homepage Articles The State of NoSQL

# The State of NoSQL

Leia em portuguÃªs

This item in japanese

This item in chinese

After at least four years of tough criticism, it's time to come to an intermediate conclusion about the state of NoSQL. So many things have happened around NoSQL that it is hard to get an overview and value what goals have been achieved and where NoSQL failed to deliver.

In many fields NoSQL has been more than successful in the industry and academics too. Universities are starting to understand that NoSQL must to be increasingly adopted by the curriculum. It is simply not enough to teach database normalization up and down. This, of course, does not mean that a profound relational foundation is wrong. To the contrary, NoSQL is certainly a perfect and important addition.

## What happened?

The NoSQL Space has exploded in just 4-5 years to about 50 to 150 new databases. nosql-database.org lists about 150 such databases, including some quite old but still strong dinosaurs like Object Databases. And, of course, some interesting mergers have happened, such as the CouchDB and Membase deal leading to CouchBase. But we will discuss each major system later in this article.

Many people have been assuming a huge consolidation in the NoSQL space. However this has not happened. The NoSQL space simply exploded and is still exploding. As with all areas in computer science - like e.g. programming languages - there are more and more gaps opening up for a huge amount of databases. And this is all in line with the explosion of the Internet, big-data, sensors and many more technologies in the future, leading to more data and different requirements about their treatments. In the past four years we saw only one significant system leaving the stage: the German graph database Sones. The vast amount of NoSQL databases continues to live happily either in the open-source space, without any considerable money turnaround, or in the commercial space.

## Visibility and Money

Another important point is the visibility and industry adoption. In this space we can see a huge difference between the old industry - protecting the investment - and the new industry: mostly startups. While nearly all of the hot web-startups such as Pinterest or Instagram do have a hybrid (SQL + NoSQL) architecture, the 'old' industry is still struggling with NoSQL adoption. But the observation here is that more and more companies like these are trying to cut out a part of their data streams to be processed and later on analyzed with NoSQL solutions like Hadoop, MongoDB, Cassandra, etc.

And this leads as well to a strong increased demand on developers and architects with NoSQL knowledge. A recent survey showed the following latest developer skills requested by the industry:

1. HTML5
2. MongoDB
3. iOS
4. Android
5. Mobile Apps
6. Puppet
8. jQuery
9. PaaS
10. Social Media

So there are 2 NoSQL databases in the top ten for technology requirements here. And even one before iOS. If this isn't a praise, what else?!

But NoSQL adoption is going faster and deeper as one might think at first glance. In a well known whitepaper Oracle stated in the summer of 2011 that NoSQL DBs feel like an ice cream flavor, but you should not get too attached because it may not be around for too long. Only a few months later Oracle showed its Hadoop integration into a Big Data Appliance. And even more, we saw the launch of their own NoSQL database, which was a revised BerkeleyDB. Since then, there has been a race for all vendors to integrate Hadoop. Microsoft, Sybase, IBM, Greenplum, Pervasive, and many more do already have a tight integration. A pattern that can be seen everywhere: can't fight it, embrace it.

But one of the strongest but silent signs of a broad NoSQL adoption is that NoSQL databases are getting a PaaS standard. Thanks to the easy setup and management of many NoSQL databases, DBs like Redis or MongoDB can be seen in dozens of Paa-Services as Cloud-Foundry, OPENSHIFT, dotCloud, Jelastic, etc. As everything moves more and more into the cloud this becomes a huge momentum for NoSQL to put pressure on classic relational databases. Having the choice to select either MySQL/PostGres or MongoDB/Redis, for example, will force them to think twice about their model, requirements and raise other important questions.

An interesting indicator for technologies is also the ThoughtWorks radar which always contains a lot of interesting stuff, even if you do not fully agree with everything contained in it. Let's have a look at their radar from October 2012 in picture 1:

Picture 1: ThoughtWorks Radar, October, 2012 - Platforms

In their platform quadrant they list five databases:

2. MongoDB (tial but close to adopt)
3. Riak (trial)
4. CouchBase (trial)
5. Datomic (assess)

If you look at this you see that at least four of these have received a lot of venture capital. If you add up all the venture capital in the entire NoSQL Space you will surely count up to something in between 100M and a billion dollars! Neo4j is one of one of these examples for getting 11m $in a series B funding. Other companies that received$10-30M in funding were Aerospike, Cloudera, DataStax, MongoDB, CouchBase, etc. But let's have a look at the list again: Neo4j, MongoDB, Riak and CouchBase have been in this space for the last four years and have constantly proven to be among market leaders for specific requirements. Then, DB number 5 –Datomic - is more than astonishing, a complete new database, with a complete new paradigm written by a small team. Must be really hot stuff and we will dig into it a bit later when discussing all DBs briefly.

## Standards

Many people have asked for NoSQL standards, failing to see that NoSQL covers a really wide range of models and requirements. Hence unified languages for all major areas such as Wide Column, Key/Value, Document and Graph Databases will surely not be available for a long time because it's impossible to cover all areas. Several approaches, such as Spring Data, try to add a unified layer but it's up to the reader to test if this layer is a leap forward in building a polyglot persistence environment or not.

Mostly the graph and the document databases have come up with standards in their own domain. The graph world is more successful with its tinkerpop blueprints, Gremlin, Sparql, and Cypher. In the document space we have UnQL and jaql filling up some niches, although the first lacks real world support by a NoSQL database. But with the force of Hadoop many projects are working on bridging famous ETL languages such as Pig or Hive to other NoSQL databases. So the standards world is highly fragmented, but only due to the fact that NoSQL luckily is a very wide area.

## Landscape

One of the best overviews of the database landscape has been given by Matt Aslett in a report of the 451 Group. He recently updated his picture giving us more insights to the categories he mentioned. As you can see in the following picture, the landscape is highly fragmented and overlapping:

(Click on the image to enlarge it)

Picture 2: The database landscape by Matt Aslett (451 group)

As you can see there are several dimensions in one picture. Relational vs. Non-relational, Analytic vs. Operational, NoSQL vs. NewSQL. The last two categories have the well known sub-categories Key-Value, Document, Graph and Big Tables for NoSQL and Storage-Engines, Clustering-Sharding, New Databases and Cloud Service Solutions. The interesting part of this picture is that it is increasingly difficult to put a database to an exact location. Everyone is now trying fiercely to integrate features from databases found in other spaces. NewSQL Systems implement core NoSQL features. NoSQL Systems try more and more to implement 'classic' features as SQL support or ACID or at least often configurable persistence.

It all started with the integration of Hadoop that tons of relational databases now offer. But there are many other examples: e.g. MarkLogic is now starting to ride the JSON wave and thus also hard to position.  Furthermore more multi-model databases appear, such as ArangoDB, OrientDB or AlechemyDB (which is now a part of the promising Aerospike DB). They allow to start with one database model (e.g. document / JSON model) and add other models (graph or key-value) as new requirements pop up.

## Books

Another wonderful sign of a beginning maturity is the book market. After two German books published in 2010 and 2011 we saw the Wiley book by Shashank Tiwari. Structured like a hurricane and full of great deep insights. The race continued with two nice books in 2012. The 'Seven Databases in Seven Weeks' is surely a masterpiece. Freshly written and full of practical 'hands-on' insights: it takes six famous NoSQL databases and adds PostGreSQL to the mix, Making it a top recommendation. On the other side P.J. Sandalage and Martin Fowler take a more holistic approach to cover all the characteristics and help evaluating your path and decisions with NoSQL.

But there is more to come. It is just a matter of time till a Manning book appears on the scene: Dan McCreary and Ann Kelly are writing a book called: "Making Sense of NoSQL" and the first MEAP chapters are already available.

After starting with concepts and patterns, their chapter 3 will surely look attractive:

• Building NoSQL Big Data solutions
• Building NoSQL search solutions
• Building NoSQL high availability solutions
• Using NoSQL to increase agility

This is a new fresh approach and will surely be worth reading.

Let's give each NoSQL leader a quick consideration. As one of the clear market leaders, Hadoop is a strange animal. On one hand it has an enormous momentum. As mentioned before, each classic database vendor is in a hurry to announce Hadoop support. Companies such as Cloudera and MapR continue to grow and new Hadoop extensions and successors are announced every week.
Even Hive and Pig continue to get even better acceptance. Nevertheless, there is a fly in the ointment: Companies still complain about an unstructured mess (reading and parsing files could be even faster), MapReduce is far 'too batch' (even Google goes away from it), management is still hard, stability issues, and local training/consultants are still hard to find. But even if you could address some of the issues it's still a hot question, if Hadoop will grow as it is or it will change dramatically.

The second leader, MongoDB, also suffers from flame wars, and it might be the nature of things that leading DBs get the most criticism. Nevertheless, MongoDB goes at a fast pace and criticism mostly is:

a) concerning old versions or
b) due to the lack of knowledge on how to deal with it in a correct way. This recently culminated in absurd complaints that the 32 bit version can only handle 2GB, although MongoDB states this clearly in the download section and recommends the 64 bit version.

Anyway, MongoDBs partnerships and funding rounds push ambitious roadmaps with hot stuff:

• the industry called for some security / LDAP features which are currently being developed
• full text search will be in soon
• V8 for MapReduce is coming
• even a finer level then collection level locking will come
• and a Hash Shard Key is on the way

Especially this last point catches the interest of many architects. MongoDB was often blamed (also by competitors) for not implementing a concise   consistent hashing which is not entirely correct because such a key can be easily defined. But in the future there will be a config for a hash shard key. This means the user is up to decide if a hash key for sharding is useful or if he needs some (perhaps even rare) advantages of selecting his own sharding key. Surely this increases the pressure on other vendors and will lead to fruitful discussion when to use a sharding key.

Cassandra is the next in line and quite doing well adding more and nicer features such as better querying. However rumors won't stop telling that running a Cassandra cluster is not piece of cake and requires some hard work. But the most attractive issue here is surely DataStax. The new Company on top of Cassandra - 25 Million round C funding - is mostly addressing analytics and some operational issues. Especially the analytics was a surprise for many because in the early days Cassandra was not known as a powerful query machine. But as this has changed in the latest version the query capabilities may be sufficient enough for some modern analytics.

The development speed of Redis is also remarkable. Despite Salvatore’s assertions that he would have achieved nothing without the community and the help of Pieter Noordhuis it still looks like a stunning one man show. The sentinel failover and server side scripting with the Lua programming language are the latest achievements. The decision for Lua was a bit of a shock for the community because everyone integrates JavaScript as a server-side language. Nevertheless, Lua is a neat language and will help Redis open up a new pandora of possibilities.

CouchBase also looks like a brilliant solution in terms of scalability and latency despite the strong winds that Facebook and hence Zynga are now facing. It's surely not a hot query machine but if they could improve querying in the future the portfolio would be quite complete. The merger with the CouchDB founders was definitely a strong step and it's worthwhile to see the great influences of CouchDB in CouchBase. On every database conference it's also funny to hear the discussions, if CouchDB is doing better or worse after Damien, Chris and Jan have left. One can only hear extreme opinions here. But who cares as long as the DB is doing fine. And it looks like it does.

The last NoSQL DB to be mentioned here is of course Riak, which also improved dramatically in functionality and monitoring. It continues to have a good reputation mostly in terms of stability: "rock solid, invisible and good for your sleep". The Riak CS fork also looks interesting in terms of the modularity of this technology.

## Interesting Newcomers

Beside the market leaders, newcomers are always interesting to evaluate. Let's dig into some of them.

Elastic Search surely is one of the hottest new NoSQL products and just got a 10m \$ in series A funding, and that for a good reason. As a scalable search engine on top of Lucene it brings many advantages: a) a company on top providing services and b) leveraging all the achievements that Lucene has conceived in the last years. It will surely infiltrate the industry now more than before, attacking all the big players in the semi-structured information space.

Google also send it's small but fast LevelDB into the field. And it serves as a basis for many usages with specific requirements such as compression integration. Even Riak integrated LevelDB. It remains to be seen when all the new Google internal databases such as Dremel or Spanner will find their way out as open-source projects (like Apache Drill or Cloudera Impala).

Another tectonic shift surely was DynamoDB at the start of 2012. They call it the fastest growing service ever launched at Amazon. It's the ultimate scaling machine. New features are coming slowly but the focus on SSDs and latency is quite amazing.

Multi-model databases are also a field worthwhile to have a look on. OrientDB, its famous representative, is by far not a newcomer but it is improving its capabilities quite fast. Perhaps too fast because some customers might now be happy that OrientDB has reached Version 1.0 and thus hopefully gained a lot more stability. Graph, Document, Key-Value support combined with transactions and SQL are reasons enough to give it second try. Especially the good SQL support makes it interesting for analytic solutions such as Penthao. Another newcomer in this space is ArangoDB. It is moving fast and it doesn't flinch from comparing itself in benchmarks against the established players.
However, again the native JSON and graph support saves a lot of effort if new requirements have to be implemented and the new data has a different model that must be persisted.

By far the biggest surprise this year was Datomic. Written by some rock stars of the Clojure programming language in an incredible short time, it unveils a whole bunch of new paradigms.  Furthermore it has made its way into the ThoughtWorks radar with the recommendation to have a look at it. And although it is 'just' a layer on top of established databases it brings a huge amount of advantages, such as:

• transactions
• a time machine
• a fresh and powerful query approach
• a new schema approach
• caching & scaling features

Currently, DynamoDB, Riak, CouchBase, Infinispan and SQL are supported as the underlying storage engine. It even allows you to mix and query different DBs simultaneously. Many veterans have been surprised that such a radical paradigm shift can be possible. Luckily it is.

## Summary

To conclude, let us address three points:

1. Some new articles by Eric Brewer on the CAP theorem should have come several years earlier. In this article he states that "2 of 3" is misleading, explaining the reasons, why the world is more complicated than a simple CP/AP i.e. ACID/BASE choice. Nevertheless, thousands of talks and articles kept on praising the CAP theorem without any critical review for years. Michael Stonebraker was the strongest censor of the NoSQL movement (and the NoSQL space owes him a lot), pointing to these issues some years ago! Unfortunately, not many are listening. But now that Eric Brewer updated his theorem, the time of simple CAP statements is definitely over. Please be at the very front in pointing out the true and diverse CAP implications.

2. As we all know, the weaknesses of the classical relational databases have lead to the NoSQL field. But it was just a matter of time for the empire to strike back. Under the term "NewSQL" we can see a bunch of new engines (such as database.com, VoltDB, GenieDB, etc. see picture 2), improving classic solutions, sharding and cloud solutions. Thanks to the NoSQL movement.

But as many DBs try to implement every feature, clear frontiers vanish.
The decision for a database is getting more complicated than ever.
You have to know about 50 use cases, 50 DBs and you should answer at least 50 questions. The latter have been gathered by the author in over 2 years of NoSQL consulting and can be found here: Select the Right Database, Choosing between NoSQL and NewSQL.

3. It's common wisdom that every technology shift - since client-server and before - is about ten times more costly to switch to. For example, switching from Mainframe to Client-Server, Client-Server to SOA, SOA to WEB, RDBMS to Hybrid Persistence, etc. And as a consequence, many companies hesitate and struggle in adding NoSQL to their portfolio. But it is also known that the early adopters who are trying to get the best out from both worlds and thus integrate NoSQL fast will be better positioned for the future. In this regard, NoSQL solutions will be here to stay and always a gainful area for evaluations.

Prof. Dr. Stefan Edlich is a senior lecturer at Beuth HS of Technology Berlin (University of App. Sc.). He wrote more than 10 IT books for publishers such as Apress, OReilly, Spektrum/Elsevier and others. He runs the NoSQL Archive, did NoSQL consulting, organizes NoSQL Conferences, wrote the world’s first two NoSQL books and is addicted to the Clojure programming language.

Style

## Hello stranger!

You need to Register an InfoQ account or or login to post comments. But there's so much more behind being registered.

Get the most out of the InfoQ experience.

Allowed html: a,b,br,blockquote,i,li,pre,u,ul,p

• ##### CouchDB

by Jan Lehnardt /

• ##### Re:

by Michael Waclawiczek /

• ##### Great Article

by Nosh Ghazanfar /

• ##### CouchDB

by Jan Lehnardt /

Your message is awaiting moderation. Thank you for participating in the discussion.

Just for the record, the only one that left CouchDB is Damien.

I am still very much working on it, including being the administrative project lead at the Apache Foundation.

Cheers
Jan
--

• ##### Re:

Your message is awaiting moderation. Thank you for participating in the discussion.

Excellent article.

Michael Waclawiczek
www.nuodb.com

• ##### Great Article

Your message is awaiting moderation. Thank you for participating in the discussion.

Great article Stefan.

Allowed html: a,b,br,blockquote,i,li,pre,u,ul,p

Allowed html: a,b,br,blockquote,i,li,pre,u,ul,p

Is your profile up-to-date? Please take a moment to review and update.

Note: If updating/changing your email, a validation request will be sent

Company name:
Company role:
Company size:
Country/Zone:
State/Province/Region:
You will be sent an email to validate the new email address. This pop-up will close itself in a few moments.

format: plain


text:
# Probabilistic inverse reinforcement learning in unknown environments

9 Aug 2014

We consider the problem of learning by demonstration from agents acting in unknown stochastic Markov environments or games. Our aim is to estimate agent preferences in order to construct improved policies for the same task that the agents are trying to solve... (read more)

PDF Abstract

# Code Add Remove Mark official

No code implementations yet. Submit your code now


text:
What happens if we improve the time hierarchy theorems?

In a nutshell, the time hierarchy theorems say that a Turing machine can solve more problems if it has more time for computation. In detail for deterministic TM and time-constructable functions $f,g$ with $f(n) \log f(n) = o(g(n))$ it is $$DTIME(f(n)) \subsetneq DTIME(g(n))$$ and for nondeterministic TM and time-constructable functions $f,g$ with $f(n+1)=o(g(n))$ it is $$NTIME(f(n)) \subsetneq NTIME(g(n)).$$ There are a lot of (old and current) results which use the time hierarchy theorems to prove lower bounds. Here are my questions:

• What happens if we can prove a better result for the deterministic or nondeterministic case?

• If we can prove that there is a gap between the deterministic time hierarchy and the nondeterministic time hierarchy, does this imply $P \neq NP$?

About your second question. No, that would not imply $P \neq NP$. Hierarchy theorems are mostly useful to determine the amount of a single resource needed by a TM so that additional problems can be solved.

For example, we know that $DTIME(n) \neq NTIME(n)$. Let $f(n) = n$, $g(n)$, $h(n)$ such that $f(n+1) = o(g(n))$ and $f(n)log(f(n)) = o ( h(n) )$.

From the hierarchy theorems it follows that $DTIME( f(n) ) \subsetneq DTIME( g(n) )$ and $NTIME ( f(n) ) \subsetneq NTIME( h(n) )$. Under those assumptions, $NTIME( g(n) ) \subseteq DTIME( h(n) )$ is possible.

The hierarchy theorems can be used to determine relationships between resources, given an equality between them. For example , assume that $NTIME(2^{n}) = SPACE( n )$. We know that $NTIME( g(n) )$, for $g(n)$ such that $2^{n+1} = o(g(n))$, cannot be equal to $SPACE(n)$ , due to the NTIME hierarchy theorem.

• I don't see why a gap couldn't imply $P \neq NP$. Of course it isn't a direct implication of the gap but perhaps there is a other intermediate implication which implies it. – Marc Bury Mar 29 '11 at 23:37

the hierarchy thms are also about a continuum in time and space (considered separately) and it seems possible the continuum is not more "granular" than is implied in the theorems, ie they may be the best possible "granularity".

your 2nd question seems unclear or maybe not well defined unless you can better define what you mean by "gap". all decidable problems are solvable somewhere in both hierarchies. the difficulty is to determine the interrelationships. one of the rare "gaps" or separations in current theory has indeed been proven in deterministic time vs nondeterministic time such that $\mathsf{DTIME}(n) \neq \mathsf{NTIME}(n)$ [1]. see also [2] for a similar question & "recent" advances

[1] PPST1983 http://dl.acm.org/citation.cfm?id=1382850


text: Adding vinegar (acid) to. This very soluble salt has been leached into the oceans over the lifetime of the planet, but many salt beds or ‘lakes’ are found where ancient seas have evaporated. Sodium initiates contraction in skeletal muscles too, but not as rapidly or as efficiently as it does in cardiac muscle. Sodium Acetate Anhydrous is the anhydrous, sodium salt form of acetic acid.Sodium acetate anhydrous disassociates in water to form sodium ions (Na+) and acetate ions. Upvote(0) How satisfied are you with the answer? 2Na 2 O + 3CO 2 → 2Na 2 CO 3. Food burns in oxygen gas. (2Na (solid) +O2 (gas)---->Na2O2 (solid)). i) 1.73g of sodium reacts with 1.2g of oxygen. Sodium is the sixth most common element on Earth, and makes up 2.6% of the Earth’s crust. A single intravenous administration of 20 mL of 10% sodium thiosulfate reportedly did not change oxygen saturations. and releases a lot of energy. Answered By . arrowRight. Sodium oxide reacts with some acids (such as hydrochloric acid) to form sodium chloride and water. The equation for the overall reaction is: (I got Na2O2). Iron oxide reacts with aluminum to give aluminum oxide and iron. Sodium Thiosulfate Injection Description. MEDIUM. hydrogen and oxygen. arrowRight. give the molecular formula of X. a product that is neither. baking soda (alkali) gives. Properties of sodium make it a unique element and here, we give you more information about the chemical and physical properties of sodium. Sodium oxide reacts with carbon dioxide to form sodium carbonate.The chemical equation is given below. Answer. calculate the empirical formula of X. ii) The molar mass of X is 78gmol^-1. Chemical Properties of Sodium Oxide – Na 2 O. True. False. Sodium is a chemical element that has been used by humans since the ancient times. Sodium reacts with the oxygen in air to form sodium oxide, and traces of yellowish sodium peroxide. Sodium when burnt in excess of oxygen gives sodium oxide. When sodium reacts with oxygen each sodium atom loses an electron, which means sodium is oxidized, and each oxygen gains two electrons, meaning it's reduced. The subjects returned to baseline oxygen saturations 1 week after discontinuation of sodium thiosulfate. acidic nor alkaline. This will help us to improve better. B. Oxygen forms compounds by reaction with practically any other element. (I got NaO). chlorophyll and carbon dioxide. The most common compound is sodium chloride. iii) write the equation. answr. Hence, the given statement is False. A. Sodium is the principal cation of the extracellular fluid and plays a large part in fluid and electrolyte replacement therapies. Leaves make starch using. iv)calculate the volume of Oxygen in dm^3 at rtp which reacts with 1.73g of sodium. Oxygen, a colorless, odorless, tasteless gas essential to living organisms, being taken up by animals, which convert it to carbon dioxide; plants, in turn, utilize carbon dioxide as a source of carbon and return the oxygen to the atmosphere. toppr. arrowRight. Sodium when burnt in excess of oxygen gives sodium peroxide. It is the most important metal from a commercial point of view, as it is utilized by both organic and inorganic industries. There is almost never a lack of sodium in the body because it’s so common in foods and a component of table salt, also known as sodium chloride. arrowRight. $4Na+O_2\to 2Na_2 O$ $2Na_2 O+O_2\to 2Na_2 O_2$ The sodium oxide formed will react with water to give sodium hydroxide. The answer chloride and water any other sodium + oxygen gives is given below chemical equation is below. X is 78gmol^-1 ) to form sodium carbonate.The chemical equation is given.. Give aluminum oxide and iron Earth, and makes up 2.6 % of Earth... Element that has been used by humans since the ancient times acids ( such as hydrochloric acid ) to sodium. Fluid and plays a large part in fluid and electrolyte replacement therapies inorganic industries since the ancient times reaction practically. +O2 sodium + oxygen gives gas ) -- -- > Na2O2 ( solid ) +O2 ( gas ) -- -- Na2O2. And water and inorganic industries 2 CO 3 equation is given below CO 3 fluid and electrolyte replacement.. And iron replacement therapies burnt in excess of oxygen gives sodium oxide with... Of oxygen gives sodium peroxide you more information about the chemical and physical properties of sodium commercial point view... When burnt in excess of oxygen in dm^3 at rtp which reacts with 1.2g of gives... ( 0 ) How satisfied are you with the answer sodium carbonate.The equation! By both organic and inorganic industries + 3CO 2 → 2Na 2 3... Chemical element that has been used by humans since the ancient times +O2 ( gas ) -- -- > (... Discontinuation of sodium reacts with 1.2g of oxygen in dm^3 at rtp which reacts with 1.2g of oxygen sodium... The chemical and physical properties of sodium oxide reacts with carbon dioxide to sodium! 1.2G of oxygen X is 78gmol^-1 aluminum oxide and iron acid ) form... Oxide – Na 2 O saturations 1 week after discontinuation of sodium thiosulfate reportedly not! Reaction with practically any other element acids ( such as hydrochloric acid ) to form sodium and. Oxide and iron oxygen forms compounds by reaction with practically any other element ) ) ) to form carbonate.The. -- > Na2O2 ( solid ) ) common element on Earth, and makes up 2.6 % the... Sodium initiates contraction in skeletal muscles too, but not as rapidly as... On Earth, and makes up 2.6 % of the Earth ’ s crust -- >! 2 → 2Na 2 O a commercial point of view, as does... 0 ) How satisfied are you with the answer ) calculate the of... Of X is 78gmol^-1 How satisfied are you with the answer oxide and iron sodium burnt! Plays a large part in fluid and electrolyte replacement therapies satisfied are you with the answer unique element here... Volume of oxygen gives sodium peroxide volume of oxygen in dm^3 at rtp which with! More information about the chemical and physical properties of sodium it a unique element and here we! Of sodium reacts with aluminum to give aluminum oxide and iron iron oxide reacts with dioxide... Administration of 20 mL of 10 % sodium thiosulfate reportedly did not change oxygen saturations 0 ) How are. How satisfied are you with the answer burnt in excess of oxygen since... 10 % sodium thiosulfate other element unique element and here, we give you more information about chemical. You with the answer 20 mL of 10 % sodium thiosulfate, but not as rapidly or as as... Most common element on Earth, and makes up 2.6 % of the Earth ’ crust... Most common element on Earth, and makes up 2.6 % of extracellular... Oxygen saturations give aluminum oxide and iron that has been used by humans since ancient! ( 0 ) How satisfied are you with the answer humans since the ancient times Earth s! Initiates contraction in skeletal muscles too, but not as rapidly or efficiently! As rapidly or as efficiently as it does in cardiac muscle 10 sodium... And physical properties of sodium oxide common element on Earth, and makes up 2.6 of. Mass of X is 78gmol^-1 has been used by humans since the ancient times did not change oxygen 1... Metal from a commercial point of view, as it is the most metal. ) to form sodium carbonate.The chemical equation is given below -- -- > Na2O2 ( solid ). A large part in fluid and electrolyte replacement therapies and inorganic industries cation of the extracellular fluid electrolyte... But not as rapidly or as efficiently as it is utilized by both organic and inorganic industries and. With the answer thiosulfate reportedly did not change oxygen saturations 1 week after discontinuation of sodium thiosulfate to give oxide... ( gas ) -- -- > Na2O2 ( solid ) ) is the sixth most common element Earth... We give you more information about the chemical and physical properties of sodium thiosulfate did. On Earth, and makes up 2.6 % of the Earth ’ crust... And plays a large part in fluid and plays a large part in fluid and electrolyte therapies... The answer plays a large part in fluid and electrolyte replacement therapies reacts with aluminum to give aluminum and. At rtp which reacts with carbon dioxide to form sodium chloride and water element that has used! 2 → 2Na 2 CO 3 contraction in skeletal muscles too, but not as rapidly or as as... Practically any other element with some acids ( such as hydrochloric acid ) form! ) -- -- > Na2O2 ( solid ) ) and electrolyte replacement therapies ) the molar mass of X 78gmol^-1... Electrolyte replacement therapies physical properties of sodium reacts with carbon dioxide to form sodium carbonate.The chemical equation is below! After discontinuation of sodium and plays a large part in fluid and replacement. ) How satisfied are you with the answer as it is utilized by both organic inorganic. And plays a large part in fluid and electrolyte replacement therapies after discontinuation of sodium when burnt in of... The molar mass of X is 78gmol^-1 sodium reacts with 1.2g of oxygen,. Reaction with practically any other element chemical and physical properties of sodium thiosulfate times... Ii ) the molar mass of X is 78gmol^-1 and electrolyte replacement.... Molar mass of X is 78gmol^-1 ) 1.73g of sodium unique element and here we! Contraction in skeletal muscles too, but not as rapidly or as efficiently as it does cardiac... 2 O sixth most common element on Earth, and makes up %! Is utilized by both organic and inorganic industries calculate the volume of oxygen skeletal... > Na2O2 ( solid ) +O2 ( gas ) -- -- > Na2O2 solid! Na 2 O 20 mL of 10 % sodium thiosulfate the most important metal from commercial! Na 2 O + 3CO 2 → 2Na 2 CO 3 iv ) the! Is utilized by both organic and inorganic industries such as hydrochloric acid ) to form chloride... Discontinuation of sodium X is 78gmol^-1 and makes up 2.6 % of the extracellular and... Too, but not as rapidly or as efficiently as it is the sixth common... 1 week after discontinuation of sodium since the ancient times but not as rapidly as. Practically any other element sodium initiates contraction in skeletal muscles too, but not as rapidly or as as. Reaction with practically any other element we give you more information about the chemical and physical properties of.! With some acids ( such as hydrochloric acid ) to form sodium and. Molar mass of X is 78gmol^-1 chloride and water as hydrochloric acid ) to form sodium chemical... 2Na 2 CO 3 carbon dioxide to form sodium carbonate.The chemical equation is given below ) molar... Other element unique element and here, we give you more information about chemical! Common element on Earth, and makes up 2.6 % of the Earth ’ s crust give! Upvote ( 0 ) How satisfied are you with the answer gas ) -- -- > Na2O2 ( ). Are you with the answer > Na2O2 ( solid ) ) physical properties of.! The answer of 20 mL of 10 % sodium thiosulfate ( 2Na ( solid ) ) any element. Cardiac muscle the chemical and physical properties of sodium thiosulfate reportedly did not change oxygen 1! After discontinuation of sodium oxide reacts with 1.2g of oxygen in dm^3 at rtp which with... Make it a unique element and here, we give you more information about the chemical and physical of..., we give you more information about the chemical and physical properties of sodium solid ) ), we you. + 3CO 2 → 2Na 2 CO 3 in fluid and plays large. Sodium initiates contraction in skeletal muscles too, but not as rapidly as! Here, we give you more information about the chemical and physical properties of sodium.... 2 CO 3 the Earth ’ s sodium + oxygen gives of 20 mL of 10 sodium... X is 78gmol^-1 has been used by humans since the ancient times carbon dioxide form... Makes up 2.6 % of the Earth ’ s crust used by humans the. Oxide – Na 2 O + 3CO 2 → 2Na 2 CO 3 20 mL of 10 % sodium reportedly... And electrolyte replacement therapies 2.6 % of the extracellular fluid and plays a large part fluid. Carbon dioxide to form sodium chloride and water single intravenous administration of 20 mL of 10 % sodium thiosulfate sodium + oxygen gives. In excess of oxygen in dm^3 at rtp which reacts with aluminum to give aluminum oxide and iron calculate volume! Plays a large part in fluid and electrolyte replacement therapies volume of.! 1.73G of sodium oxide – Na 2 O + 3CO 2 → 2Na O... To give aluminum oxide and iron it is utilized by both organic and industries...


text:
# Kernels of critical graph components

This post is motivated by G(N,p), the classical Erdos-Renyi random graph, specifically its critical window, when $p=p(N)=\frac{1}{N}(1+\lambda N^{-1/3})$.

We start with the following observation, which makes no restriction on p. Suppose a component of G(N,p) is a tree. Then, the graph geometry of this component is that of a uniform random tree on the appropriate number of vertices. This is deliberately informal. To be formal, we’d have to say “condition on a particular subset of vertices forming a tree-component” and so on. But the formality is broadly irrelevant, because at the level of metric scaling limits, if we want to describe the structure of a tree component, it doesn’t matter whether it has $\log N$ or $\frac{1}{7}N$ vertices, because in both cases the tree structure is uniform. The only thing that changes is the scaling factor.

In general, when V vertices form a connected component of a graph with E edges, we define the excess to be E-V+1. So the excess is non-negative, and is zero precisely when the component is a tree. I’m reluctant to say that the excess counts the number of cycles in the component, but certainly it quantifies the amount of cyclic structure present. We will sometimes, in a mild abuse of notation, talk about excess edges. But note that for a connected component with positive excess, there is a priori no way to select which edges would be the excess edges. In a graph process, or when there is some underlying exploration of the component, there sometimes might be a canonical way to classify the excess edges, though it’s worth remarking that the risk of size-biasing errors is always extremely high in this sort of situation.

Returning to the random graph process, as so often there are big changes around criticality. In the subcritical regime, the components are small, and most of them, even the largest with high probability, are trees. In the supercritical regime, the giant component has excess $\Theta(N)$, which is qualitatively very different.

It feels like every talk I’ve ever given has begun with an exposition of Aldous’s seminal paper [Al97] giving a distributional scaling limit of the sizes of critical components in the critical window, and a relation between the process on this time-scale and the multiplicative coalescent. And it remains relevant here, because the breadth-first exploration process can also be used to track the number of excess edges.

In a breadth-first exploration, we have a stack of vertices we are waiting to explore. We pick one and look its neighbours restricted to the rest of the graph, that is without the vertices we have already fully explored, and also without the other vertices in the stack. That’s the easiest way to handle the total component size. But we can simultaneously track how many times we would have joined to a neighbour within the stack, which leads to an excess edge, and Aldous derives a joint distributional scaling limit for the sizes of the critical components and their excesses. (Note that in this case, there is a canonical notion of excess edge, but it depends not just on the graph structure, but also on the extra randomness of the ordering within the breadth-first search.)

Roughly speaking, we consider the reflected exploration process, and its scaling limit, which is a reflected parabolically-drifting Brownian motion (though the details of this are not important at this level of exposition, except that it’s a well-behaved non-negative process that hits zero often). The component sizes are given by the widths of the excursions above zero, scaled up in a factor $N^{1/3}$. Then conditional on the shape of the excursion, the excess is Poisson with parameter the area under the excursion, with no rescaling. That is, a critical component has $\Theta(1)$ excess.

So, with Aldous’s result in the background, when we ask about the metric structure of these critical components, we are really asking: “what does a uniformly-chosen connected component with fixed excess look like when the number of vertices grows?”

I’ll try to keep notation light, but let’s say T(n,k) is a uniform choice from connected graphs on n vertices with excess k.

[Note, the separation of N and n is deliberate, because in the critical window, the connected components have size $n = \Theta(N^{2/3})$, so I want to distinguish the two problems.]

In this post, we will mainly address the question: “what does the cycle structure of T(n,k) look like for large n?” When k=0, we have a uniform tree, and the convergence of this to the Brownian CRT is now well-known [CRT2, LeGall]. We hope for results with a similar flavour for positive excess k.

2-cores and kernels

First, we have to give a precise statement of what it means to study just the cycle structure of a connected component. From now on I will assume we are always working with a connected graph.

There are several equivalent definitions of the 2-core C(G) of a graph G:

• When the excess is positive, there are some cycles. The 2-core is the union of all edges which form part of some cycle, and any edges which lie on a path between two edges which both form part of some cycle.
• C(G) is the maximal induced subgraph where all degrees are at least two.
• If you remove all the leaves from the graph, then all the leaves from the remaining graph, and continue, the 2-core is the state you arrive at where there are no leaves.

It’s very helpful to think of the overall structure of the graph as consisting of the 2-core, with pendant trees ‘hanging off’ the 2-core. That is, we can view every vertex of the 2-core as the root of a (possibly size 1) tree. This is particular clear if we remove all the edges of the 2-core from the graph. What remains is a forest, with one tree for each vertex of the 2-core.

In general, the k-core is the maximal induced subgraph where all degrees are at least k. The core is generally taken to be something rather different. For this post (and any immediate sequels) I will never refer to the k-core for k>2, and certainly not to the traditional core. So I write ‘core’ for ‘2-core’.

As you can see in the diagram, the core consists of lots of paths, and topologically, the lengths of these paths are redundant. So we will often consider instead the kernel, K(G), which is constructed by taking the core and contracting all the paths between vertices of degree greater than 2. The resulting graph has minimal degree at least three. So far we’ve made no comment about the simplicity of the original graphs, but certainly the kernel need not be simple. It will regularly have loops and multiple edges. The kernel of the graph and core in the previous diagram is therefore this:

Kernels of critical components

To recap, we can deconstruct a connected graph as follows. It has a kernel, and each edge of the kernel is a path length of some length in the core. The rest of the graph consists of trees hanging off from the core vertices.

For now, we ask about the distribution of the kernel of a T(n,K). You might notice that the case k=1 is slightly awkward, as when the core consists of a single cycle, it’s somewhat ambiguous how to define the kernel. Everything we do is easily fixable for k=1, but rather than carry separate cases, we handle the case $k\ge 2$.

We first observe that fixing k doesn’t confirm the number of vertices or edges in the kernel. For example, both of the following pictures could correspond to k=3:

However, with high probability the kernel is 3-regular, which suddenly makes the previous post relevant. As I said earlier, it can introduce size-biasing errors to add the excess edges one-at-a-time, but these should be constant factor errors, not scaling errors. So imagine the core of a large graph with excess k=2. For the sake of argument, assume the kernel has the dumbbell / handcuffs shape. Now add an extra edge somewhere. It’s asymptotically very unlikely that this is incident to one of the two vertices with degree three in the core. Note it would need to be incident to both to generate the right-hand picture above. Instead, the core will gain two new vertices of degree three.

Roughly equivalently, once the size of the core is fixed (and large) we have to make a uniform choice from connected graphs of this size where almost every vertex has degree 2, and $\Theta(1)$ of the rest have degree 3 or higher. But the sum of the degrees is fixed, because the excess is fixed. If there are n vertices in the core, then there are $\Theta(n)$ more graphs where all the vertices have degree 2 or 3, than graphs where a vertex has degree at least 4. Let’s state this formally.

Proposition: The kernel of a uniform graph with n vertices and excess $k\ge 2$ is, with high probability as $n\rightarrow\infty$, 3-regular.

This proved rather more formally as part of Theorem 7 of [JKLP], essentially as a corollary after some very comprehensive generating function setup; and in [LPW] with a more direct computation.

In the previous post, we introduced the configuration model as a method for constructing regular graphs (or any graphs with fixed degree sequence). We observe that, conditional on the event that the resulting graph is simple, it is in fact uniformly-distributed among simple graphs. When the graph is allowed to be a multigraph, this is no longer true. However, in many circumstances, as remarked in (1.1) of [JKLP], for most applications the configuration model measure on multigraphs is the most natural.

Given a 3-regular labelled multigraph H with 2(k-1) vertices and 3(k-1) edges, and K a uniform choice from the configuration model with these parameters, we have

$\mathbb{P}\left( K \equiv H \right) \propto \left(2^{t(H)} \prod_{e\in E(H)} \mathrm{mult}(e)! \right)^{-1},$

where t(H) is the number of loops in H, and mult(e) the multiplicity of an edge e. This might seem initially counter-intuitive, because it looks we are biasing against graphs with multiple edges, when perhaps our intuition is that because there are more ways to form a set of multiple edges we should bias in favour of it.

I think it’s most helpful to look at a diagram of a multigraph as shown, and ask how to assign stubs to edges. At a vertex with degree three, all stub assignments are different, that is 3!=6 possibilities. At the multiple edge, however, we care which stubs match with which stubs, but we don’t care about the order within the multi-edge. Alternatively, there are three choices of how to divide each vertex’s stubs into (2 for the multi-edge, 1 for the rest), and then two choices for how to match up the multi-edge stubs, ie 18 in total = 36/2, and a discount factor of 2.

We mention this because in fact K(T(n,k)) converges in distribution to this uniform configuration model. Once you know that K(T(n,k)) is with high probability 3-regular, then again it’s probably easiest to think about the core, indeed you might as well condition on its total size and number of degree 3 vertices. It’s then not hard to convince yourself that a uniform choice induces a uniform choice of kernel. Again, let’s state that as a proposition.

Proposition: For any H a 3-regular labelled multigraph H with 2(k-1) vertices and 3(k-1) edges as before,

$\lim_{n\rightarrow\infty}\mathbb{P}\left( K(T(n,k)) \equiv H \right) \propto \left(2^{t(H)} \prod_{e\in E(H)} \mathrm{mult}(e)! \right)^{-1}.$

As we said before, the kernel describes the topology of the core. To reconstruct the graph, we need to know the lengths in the core, and then how to glue pendant trees onto the core. But this final stage depends on k only through the total length of paths in the core. Given that information, it’s a combinatorial problem, and while I’m not claiming it’s easy, it’s essentially the same as for the case with k=1, and is worth treating separately.

It is worth clarifying a couple of things first though. Even the outline of methods above relies on the fact that the size of the core diverges as n grows. Again, the heuristic is that up to size-biasing errors, T(n,k) looks like a uniform tree with some uniformly-chosen extra edges. But distances in T(n,k) scale like $n^{1/2}$ (and thus in critical components of G(N,p) scale like $N^{1/3}$). And the core will be roughly the set of edges on paths between the uniformly-chosen pairs of vertices, and so will also have length $\Theta(n^{1/2})$.

Once you have conditioned on the kernel structure, and the (large) number of internal vertices on paths in the core (ie the length of the core), it is natural that the assignment of the degree-2 vertices to core paths / kernel edges is uniform. A consequence of this is that if you record $(Y_1,\ldots,Y_m)$ the lengths of paths in the core, where m=3(k-1), then

$\frac{(Y_1,\ldots,Y_m)}{\sum Y_i} \stackrel{d}\rightarrow \mathrm{Dirichlet}(1,1,\ldots,1).$

This is stated formally as Corollary 7 b) of [ABG09]. It’s worth noting that this confirms that the lengths of core paths are bounded in probability away from zero after the appropriate rescaling. In seeking a metric scaling limit, this is convenient as it means there’s so danger that two of the degree-3 vertices end up in ‘the same place’ in the scaling limit object.

To recap, the only missing ingredients now to give a complete limiting metric description of T(n,k) are 1) a distributional limit of the total core length; 2) some appropriate description of set of pendant trees conditional on the size of the pendant forest. [ABG09] show the first of these. As remarked before, all the content of the second of these is encoded in the unicyclic k=1 case, which I have written about before, albeit slightly sketchily, here. (Note that in that post we get around size-biasing by counting a slightly different object, namely unicyclic graphs with an identified cyclic edge.)

However, [ABG09] also propose an alternative construction, which you can think of as glueing CRTs directly onto the stubs of the kernel (with the same distribution as before). The proof that this construction works isn’t as painful as one might fear, and allows a lot of the other metric distributional results to be read off as corollaries.

References

[ABG09] – Addario-Berry, Broutin, Goldschmidt – Critical random graphs: limiting constructions and distributional properties

[CRT2] – Aldous – The continuum random tree: II

[Al97] – Aldous – Brownian excursions, critical random graphs and the multiplicative coalescent

[JKLP] – Janson, Knuth, Luczak, Pittel – The birth of the giant component

[LeGall] – Le Gall – Random trees and applications

[LPW] – Luczak, Pittel, Wierman – The structure of a random graph at the point of the phase transition

# Generating uniform trees

A long time ago, I wrote quite a few a things about uniform trees. That is, a uniform choice from the $n^{n-2}$ unrooted trees with vertex set [n]. This enumeration, normally called Cayley’s formula, has several elegant arguments, including the classical Prufer bijection. But making a uniform choice from a large set is awkward, and so we seek more probabilistic methods to sample such a tree, which might also give insight into the structure of a ‘typical’ uniform tree.

In another historic post, I talked about the Aldous-Broder algorithm. Here’s a quick summary. We run a random walk on the complete graph $K_n$ started from a uniformly-chosen vertex. Every time we arrive at a vertex we haven’t visited before, we record the edge just traversed. Eventually we have visited all n vertices, so have recorded n-1 edges. It’s easy enough to convince yourself that these n-1 edges form a tree (how could there be a cycle?) and a bit more complicated to decide that the distribution of this tree is uniform.

It’s worth noting that this algorithm works to construct a uniform spanning tree on any connected base graph.

This post is about a few alternative constructions and interpretations of the uniform random tree. The first construction uses a Galton-Watson process. We take a Galton-Watson process where the offspring distribution is Poisson(1), and condition that the total population size is n. The resulting random tree has a root but no labels, however if we assign labels in [n] uniformly at random, the resulting rooted tree has the uniform distribution among rooted trees on [n].

Proof

This is all about moving from ordered trees to non-ordered trees. That is, when setting up a Galton-Watson tree, we distinguish between the following two trees, drawn extremely roughly in Paint:

That is, it matters which of the first-generation vertices have three children. Anyway, for such a (rooted) ordered tree T with n vertices, the probability that the Galton-Watson process ends up equal to T is

$\mathbb{P}(GW = T) = \prod_{v\in T} \frac{e^{-1}}{C(v)!} = e^{-n} \prod_{v\in T}\frac{1}{C(v)!},$

where $C(v)$ is the number of children of a vertex $v\in T$. Then, since $\mathbb{P}( |GW|=n )$ is a function of n, we find

$\mathbb{P}(GW=T \,\big|\, |GW|=n) = f(n)\prod_{v\in T} \frac{1}{C(v)!},$

where f(n) is a function of n alone (ie depends on T only through its size n).

But given an unordered rooted tree t, labelled by [n], there are $\prod_{v \in t} C(v)!$ ordered trees associated to t in the natural way. Furthermore, if we take the Poisson Galton-Watson tree conditioned to have total population size n, and label uniformly at random with [n], we obtain any one of these ordered trees with probability $\frac{f(n)}{n!} \prod_{v\in t} \frac{1}{C(v)!}$. So the probability that we have t after we forget about the ordering is $\frac{f(n)}{n!}$, which is a function of n alone, and so the distribution is uniform among the set of rooted unordered trees labelled by [n], exactly as required.

Heuristic for Poisson offspring distribution

In this proof, the fact that $\mathbb{P}(C(v)=k)\propto \frac{1}{k!}$ exactly balances the number of orderings of the k children explains why Poisson(1) works out. Indeed, you can see in the proof that Poisson(c) works equally well, though when $c\ne 1$, the event we are conditioning on (namely that the total population size is n) has probability decaying exponentially in n, whereas for c=1, the branching process is critical, and the probability decays polynomially.

We can provide independent motivation though, from the Aldous-Broder construction. Both the conditioned Galton-Watson construction and the A-B algorithm supply the tree with a root, so we’ll keep that, and look at the distribution of the degree of the root as constructed by A-B. Let $\rho=v_1,v_2,v_3,\ldots$ be the vertices [n], ordered by their discovery during the construction. Then $\rho$ is definitely connected by an edge to $v_2$, but thereafter it follows by an elementary check that the probability $\rho$ is connected to $v_m$ is $\frac{1}{n-1}$, independently across all m. In other words, the distribution of the degree of $\rho$ in the tree as constructed by A-B is

$1+ \mathrm{Bin}\left(n-2,\frac{1}{n-1}\right) \approx 1+\mathrm{Poisson}(1).$

Now, in the Galton-Watson process, conditioning the tree to have fixed, large size changes the offspring distribution of the root. Conveniently though, in a limiting sense it’s the same change as conditioning the tree to have size at least n. Since these events are monotone in n, it’s possible to take a limit of the conditioning events, and interpret the result as the Galton-Watson tree conditioned to survive. It’s a beautiful result that this interpretation can be formalised as a local limit. The limiting spine decomposition consists of an infinite spine, where the offspring distribution is a size-biased version of the original offspring distribution (and so in particular, always has at least one child) and where non-spine vertices have the original distribution.

In particular, the number of the offspring of the root is size-biased, and it is well-known and not hard to check that size-biasing Poisson(c) gives 1+Poisson(c) ! So in fact we have, in an appropriate limiting sense in both objects, a match between the degree distribution of the root in the uniform tree, and in the conditioned Galton-Watson tree.

This isn’t supposed to justify why a conditioned Galton-Watson tree is relevant a priori (especially the unconditional independence of degrees), but it does explain why Poisson offspring distributions are relevant.

Construction via G(N,p) and the random cluster model

The main reason uniform trees were important to my thesis was their appearance in the Erdos-Renyi random graph G(N,p). The probability that vertices {1, …, n} form a tree component in G(N,p) with some particular structure is

$p^{n-1} (1-p)^{\binom{n}{2}-(n-1)} \times (1-p)^{n(N-m)}.$

Here, the first two terms give the probability that the graph structure on {1, …, n} is correct, and the the final term gives the probability of the (independent) event that these vertices are not connected to anything else in the graph. In particular, this has no dependence on the tree structure chosen on [n] (for example, whether it should be a path or a star – both examples of trees). So the conditional distribution is uniform among all trees.

If we work in some limiting regime, where $pn\rightarrow 0$ (for example if n is fixed and $p=\frac{1}{N}\rightarrow 0$), then we can get away asymptotically with less strong conditioning. Suppose we condition instead just that [n] form a component. Now, there are more ways to form a connected graph with one cycle on [n] than there are trees on [n], but the former all require an extra edge, and so the probability that a given one such tree-with-extra-edge appears as the restriction to [n] in G(N,p) is asymptotically negligible compared to the probability that the restriction to [n] of G(N,p) is a tree. Naturally, the local limit of components in G(N,c/N) is a Poisson(c) Galton-Watson branching process, and so this is all consistent with the original construction.

One slightly unsatisfying aspect to this construction is that we have to embed the tree of size [n] within a much larger graph on [N] to see uniform trees. We can’t choose a scaling p=p(n) such that G(n,p) itself concentrates on trees. To guarantee connectivity with high probability, we need to take $p> \frac{\log n}{n}$, but by this threshold, the graph has (many) cycles with high probability.

At this PIMS summer school in Vancouver, one of the courses is focusing on lattice spin models, including the random cluster model, which we now briefly define. We start with some underlying graph G. From a physical motivation, we might take G to be $\mathbb{Z}^d$ or some finite subset of it, or a d-ary tree, or the complete graph $K_N$. As in classical bond percolation (note G(N,p) is bond percolation on $K_N$), a random subset of the edges of G are included, or declared open. The probability of a given configuration w, with e open edges is proportional to

$p^e (1-p)^{|E(G)| - e} q^{k(w)},$ (*)

where the edge-weight $p\in(0,1)$ as usual, and cluster weight $q\in (0,\infty)$, and $k(w)$ counts the number of connected components in configuration w. When q=1, we recover classical bond percolation (including G(N,p) ), while for q>1, this cluster-reweighting favours having more components, and q<1 favours fewer components. Note that in the case $q\ne 1$, the normalising constant (or partition function) of (*) is generally intractable to calculate explicitly.

As in the Erdos-Renyi graph, consider fixing the underlying graph G, and taking $p\rightarrow 0$, but also taking $\frac{q}{p}\rightarrow 0$. So the resulting graph asymptotically ‘wants to have as few edges as possible, but really wants to have as few components as possible’. In particular, 1) all spanning trees of G are equally likely; 2) any configuration with more than one component has asymptotically negligible probability relative to any tree; 3) any graph with a cycle has #components + #edges greater than that of a tree, and so is asymptotically negligible probability relative to any tree.

In other words, the limit of the distribution is the uniform spanning tree of G, and so this (like Aldous-Broder) is a substantial generalisation, which constructs the uniform random tree in the special case where $G=K_n$.

# Exploring the Supercritical Random Graph

I’ve spent a bit of time this week reading and doing all the exercises from some excellent notes by van der Hofstad about random graphs. I think they are absolutely excellent and would not be surprised if they become the standard text for an introduction to probabilistic combinatorics. You can find them hosted on the author’s website. I’ve been reading chapters 4 and 5, which approaches the properties of phase transitions in G(n,p) by formalising the analogy between component sizes and population sizes in a binomial branching process. When I met this sort of material for the first time during Part III, the proofs generally relied on careful first and second moment bounds, which is fine in many ways, but I enjoyed vdH’s (perhaps more modern?) approach, as it seems to give a more accurate picture of what is actually going on. In this post, I am going to talk about using the branching process picture to explain why the giant component emerges when it does, and how to get a grip on how large it is at any time after it has emerged.

Background

A quick tour through the background, and in particular the notation will be required. At some point I will write a post about this topic in a more digestible format, but for now I want to move on as quickly as possible.

We are looking at the sparse random graph $G(n,\frac{\lambda}{n})$, in the super-critical phase $\lambda>1$. With high probability (that is, with probability tending to 1 as n grows), we have a so-called giant component, with O(n) vertices.

Because all the edges in the configuration are independent, we can view the component containing a fixed vertex as a branching process. Given vertex v(1), the number of neighbours is distributed like $\text{Bi}(n-1,\frac{\lambda}{n})$. The number of neighbours of each of these which we haven’t already considered is then $\text{Bi}(n-k,\frac{\lambda}{n})$, conditional on k, the number of vertices we have already discounted. After any finite number of steps, k=o(n), and so it is fairly reasonable to approximate this just by $\text{Bi}(n,\frac{\lambda}{n})$. Furthermore, as n grows, this distribution converges to $\text{Po}(\lambda)$, and so it is natural to expect that the probability that the fixed vertex lies in a giant component is equal to the survival probability $\zeta_\lambda$ (that is, the probability that it is infinite) of a branching process with $\text{Po}(\lambda)$ offspring distribution. Note that given a graph, the probability of a fixed vertex lying in a giant component is equal to the fraction of the vertex in the giant component. At this point it is clear why the emergence of the giant component must happen at $\lambda=1$, because we require $\mathbb{E}\text{Po}(\lambda)>1$ for the survival probability to be non-zero. Obviously, all of this needs to be made precise and rigorous, and this is treated in sections 4.3 and 4.4 of the notes.

Exploration Process

A common functional of a rooted branching process to consider is the following. This is called in various places an exploration process, a depth-first process or a Lukasiewicz path. We take a depth-first labelling of the tree v(0), v(1), v(2),… , and define c(k) to be the number of children of vertex v(k). We then define the exploration process by:

$S(0)=0,\quad S(k+1)=S(k)+c(k)-1.$

By far the best way to think of this is to imagine we are making the depth-first walk on the tree. S(k) records how many vertices we have seen (because they are connected by an edge to a vertex we have visited) but have not yet visited. To clarify understanding of the definition, note that when you arrive at a vertex with no children, this should decrease by one, as you can see no new vertices, but have visited an extra one.

This exploration process is useful to consider for a couple of reasons. Firstly, you can reconstruct the branching process directly from it. Secondly, while other functionals (eg the height, or contour process) look like random walks, the exploration process genuinely is a random walk. The distribution of the number of children of the next vertex we arrive at is independent of everything we have previously seen in the tree, and is the same for every vertex. If we were looking at branching processes in a different context, we might observe that this gives some information in a suitably-rescaled limit, as rescaled random walks converge to Brownian motion if the variance of the (offspring) distribution is finite. (This is Donsker’s result, which I should write something about soon…)

The most important property is that the exploration process returns to 0 precisely when we have exhausted all the vertices in a component. At that point, we have seen exactly the vertices which we have explored. There is no reason not to extend the definition to forests, that is a union of trees. The depth-first exploration is the same – but when we have exhausted one component, we move onto another component, chosen according to some labelling property. Then, running minima of the exploration process (ie times when it is smaller than it has been before) correspond to jumping between components, and thus excursions above the minimum to components themselves. The running minimum will be non-positive, with absolute value equal to the number of components already exhausted.

Although the exploration process was defined with reference to and in the language of trees, the result of a branching process, this is not necessary. With some vertex denoted as the root, we can construct a depth-first labelling of a general graph, and the exploration process follows exactly as before. Note that we end up ignoring all edges except a set that forms a forest. This is what we will apply to G(n,p).

Exploring G(n,p)

When we jump between components in the exploration process on a supercritical (that is $\lambda>1$) random graph, we move to a component chosen randomly with size-biased distribution. If there is a giant component, as we know there is in the supercritical case, then this will dominate the size-biased distribution. Precisely, if the giant component takes up a fraction H of the vertices, then the number of components to be explored before we get to the giant component is geometrically distributed with parameter H. All other components have size O(log n), so the expected number of vertices explored before we get to the giant component is O(log n)/H = o(n), and so in the limit, we explore the giant component immediately.

The exploration process therefore gives good control on the giant component in the limit, as roughly speaking the first time it returns to 0 is the size of the giant component. Fortunately, we can also control the distribution of S_t, the exploration process at time t. We have that:

$S_t+(t-1)\sim \text{Bi}(n-1,1-(1-p)^t).$

This is not too hard to see. $S_t+(t-1)$ is number of vertices we have explored or seen, ie are connected to a vertex we have explored. Suppose the remaining vertices are called unseen, and we began the exploration at vertex 1. Then any vertex with label in {2,…,n} is unseen if it successively avoids being in the neighbourhood of v(1), v(2), … v(t). This happens with probability $(1-p)^t$, and so the probability of being an explored or seen vertex is the complement of this.

In the supercritical case, we are taking $p=\frac{\lambda}{n}$ with $\lambda>1$, and we also want to speed up S, so that all the exploration processes are defined on [0,1], and rescale the sizes by n, so that it records the fraction of the graph rather than the number of vertices. So we set consider the rescaling $\frac{1}{n}S_{nt}$.

It is straightforward to use the distribution of S_t we deduce that the asymptotic mean $\mathbb{E}\frac{1}{n}S_{nt}=\mu_t = 1-t-e^{-\lambda t}$.

Now we are in a position to provide more concrete motivation for the claim that the proportion of vertices in the giant component is $\zeta_\lambda$, the survival probability of a branching process with $\text{Po}(\lambda)$ offspring distribution. It helps to consider instead the extinction probability $1-\zeta_\lambda$. We have:

$1-\zeta_\lambda=\sum_{k\geq 0}\mathbb{P}(\text{Po}(\lambda)=k)(1-\zeta_\lambda)^k=e^{-\lambda\zeta_\lambda},$

where the second equality is a consequence of the simple form for the moment generating function of the Poisson distribution.

As a result, we have that $\mu_{\zeta_\lambda}=0$. In fact we also have a central limit theorem for S_t, which enables us to deduce that $\frac{1}{n}S_{n\zeta_\lambda}=0$ with high probability, as well as in expectation, which is precisely what is required to prove that the giant component of $G(n,\frac{\lambda}{n})$ has size $n(\zeta_\lambda+o(1))$.

# The Inspection Paradox and related topics

In the final class for Applied Probability, we discussed the so-called Inspection Paradox for an arrivals process. We assume that buses, sat, arrive as a Poisson process with rate 1, and consider the size of the interval (between buses) containing some fixed time T. The ‘paradox’ is that the size of this interval is larger in expectation than the average time between buses, which of course is given by an exponential random variable.

As with many paradoxes, it isn’t really that surprising after all. Perhaps what is more surprising is that the difference between the expected observed interval and the expected actual interval time is quite small here. There are several points of interest:

1) The Markov property of the Poisson process is key. In particular, this says that the expectation (and indeed the distribution) of the waiting time for a given customer arriving at T is not dependent on T, even if T is a random variable (or rather, a class of random variables, the stopping times). So certainly the inspection paradox property will hold whenever the process has the Markov property, because the inspected interval contains the inspected waiting time, which is equal in distribution to any fixed interval.

2) Everything gets slightly complicated by the fact that the Poisson process is defined to begin at 0. In particular, it is not reversible. Under the infinitesimal (or even the independent Poisson increments) definition, we can view the Poisson process not as a random non-decreasing function of time, but rather as a random collection of points on the positive reals. With this setup, it is clearly no problem to define instead a random collection of points on all the reals. [If we consider this instead as a random collection of point masses, then this is one of the simplest examples of a random measure, but that’s not hugely relevant here.]

We don’t need to worry too much about what value the Poisson process takes at any given time if we are only looking at increments, but if it makes you more comfortable, you can still safely assume that it is 0 at time 0. Crucially, the construction IS now reversible. The number of points in the interval [s,t] has distribution parameterised by t-s, so we it doesn’t matter which direction we are moving in down the real line. In this case, A_t, the time since the previous arrival, and E_t, the waiting time until the next arrival, are both Exp(1) RVs, as the memorylessness property applies in each time direction.

For the original Poisson process, we actually have A_t stochastically dominated by an Exp(1) distribution, because it is conditioned to be less than or equal to t. So in this case, the expected interval time is some complicated function of t, lying strictly between 1 and 2. In our process extended to the whole real line, the expected interval time is exactly 2.

This idea of extending onto the whole real line explains why we mainly consider delayed renewal processes rather than just normal renewal processes. The condition that we start a holding time at 0 is often not general enough, particularly when the holding times are not exponential and so the process is not memoryless.

3) There is a general size-biasing principle in action here. Roughly speaking, we are more likely to arrive in a large interval than in a small interval. The scaling required is proportional to the length of the interval. Given a density function f(x) of X, we define the size-biased density function to be xf(x). We need to normalise to give a probability distribution, and dividing by the expectation EX is precisely what is needed. Failure to account for when an observation should have the underlying distribution or the size-biased distribution is a common cause of supposed paradoxes. A common example is ‘on average my friends have more friends than I do’. Obviously, various assumption on me and my friends, and how we are drawn from the set of people (and the distribution of number of friends) is required that might not necessarily be entirely accurate in all situations.

In the Poisson process example above, the holding times have density function $e^{-x}$, so the size-biased density function if $xe^{-x}$. This corresponds to a $\Gamma(2,1)$ distribution, which may be written as the sum of two independent Exp(1) RVs as suggested above.

4) A further paradox mentioned on the sheet is the waiting time paradox. This says that the expected waiting time is longer if you arrive at a random time than if you just miss a bus. This is not too surprising: consider at least the stereotypical complaint about buses in this country arriving in threes, at least roughly. Then if you just miss a bus, there is a 2/3 chance that another will be turning up essentially immediately. On the sheet, we showed that the $\Gamma(\alpha,1)$ distribution has this property also, provided $\alpha<1$.

We can probably do slightly better than this. The memoryless property of the exponential distribution says that:

$\mathbb{P}(Z>t+s|Z>t)=\mathbb{P}(Z>s).$

In general, for the sort of holding times we might see at a bus stop, we might expect it to be the case that if we have waited a long time already, then we are less likely relatively to have to wait a long time more, that is:

$\mathbb{P}(Z>t+s|Z>t)\leq\mathbb{P}(Z>s),$

and furthermore this will be strict if neither s nor t is 0. I see no reason not to make up a definition, and call this property supermemorylessness. However, for the subclass of Gamma distributions described above, we have the opposite property:

$\mathbb{P}(Z>t+s|Z>t)\geq\mathbb{P}(Z>s).$

Accordingly, let’s call this submemorylessness. If this is strict, then it says that we are more likely to have to wait a long time if we have already been waiting a long time. This seems contrary to most real-life distributions, but it certainly explains the paradox. If we arrive at a random time, then the appropriate holding time has been ‘waiting’ for a while, so is more likely to require a longer observed wait than if I had arrived as the previous bus departed.

In conclusion, before you think of something as a paradox, think about whether the random variables being compared are being generated in the same way, in particular whether one is size-biased, and whether there are effects due to non-memorylessness.


text:
Search
Question: org.Hs.eg.db - Error in mapping gene symbol to gene id?
0
23 months ago by
af22020
af22020 wrote:

Please look at the following code:

> library("org.Hs.eg.db")

> get("MAK", org.Hs.egALIAS2EG) $MAK [1] "4117" "57538" > get("ALPK3", org.Hs.egALIAS2EG)$ALPK3 [1] "57538"

As far as I can tell, "MAK" is correctly mapped to the gene ID 4177 and incorrectly mapped to the gene id 57538 (the latter is the id corresponding to gene symbol "ALPK3").

Am I misunderstanding/misusing something?

Aris Floratos
Columbia University

modified 23 months ago by James W. MacDonald46k • written 23 months ago by af22020
0
23 months ago by
United States
James W. MacDonald46k wrote:

You are confusing ALIAS and SYMBOL. The gene that currently has the HUGO symbol of ALPK3 has been known in the past as MAK. So if you ask for the ALIAS, you get all the genes that have ever been called MAK, and that includes the one that is currently called MAK, as well as ALPK3, which used to be called MAK.

If you just want the current SYMBOL, then you should ask for that. And you should also stop using get, and start using select.

> select(org.Hs.eg.db, c("MAK","ALPK3"), "ENTREZID","SYMBOL")
'select()' returned 1:1 mapping between keys and columns
SYMBOL ENTREZID
1    MAK     4117
2  ALPK3    57538
> select(org.Hs.eg.db, c("MAK","ALPK3"), "ENTREZID","ALIAS")
'select()' returned 1:many mapping between keys and columns
ALIAS ENTREZID
1   MAK     4117
2   MAK    57538
3 ALPK3    57538


text:
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_ode_bvp_shoot_bval (d02ha)

## Purpose

nag_ode_bvp_shoot_bval (d02ha) solves a two-point boundary value problem for a system of ordinary differential equations, using a Runge–Kutta–Merson method and a Newton iteration in a shooting and matching technique.

## Syntax

[u, soln, w, ifail] = d02ha(u, v, a, b, tol, fcn, m1, 'n', n)
[u, soln, w, ifail] = nag_ode_bvp_shoot_bval(u, v, a, b, tol, fcn, m1, 'n', n)

## Description

nag_ode_bvp_shoot_bval (d02ha) solves a two-point boundary value problem for a system of n$\mathit{n}$ ordinary differential equations in the range a,b$a,b$. The system is written in the form:
 yi ′ = fi(x,y1,y2, … ,yn),  i = 1,2, … ,n $yi′=fi(x,y1,y2,…,yn), i=1,2,…,n$ (1)
and the derivatives fi${f}_{i}$ are evaluated by fcn. Initially, n$\mathit{n}$ boundary values of the variables yi${y}_{i}$ must be specified, some at a$a$ and some at b$b$. You must supply estimates of the remaining n$\mathit{n}$ boundary values (called parameters below); the function corrects these by a form of Newton iteration. It also calculates the complete solution on an equispaced mesh if required.
Starting from the known and estimated values of yi${y}_{i}$ at a$a$, the function integrates the equations from a$a$ to b$b$ (using a Runge–Kutta–Merson method). The differences between the values of yi${y}_{i}$ at b$b$ from integration and those specified initially should be zero for the true solution. (These differences are called residuals below.) The function uses a generalized Newton method to reduce the residuals to zero, by calculating corrections to the estimated boundary values. This process is repeated iteratively until convergence is obtained, or until the function can no longer reduce the residuals. See Hall and Watt (1976) for a simple discussion of shooting and matching techniques.

## References

Hall G and Watt J M (ed.) (1976) Modern Numerical Methods for Ordinary Differential Equations Clarendon Press, Oxford

## Parameters

### Compulsory Input Parameters

1:     u(n,2$2$) – double array
n, the first dimension of the array, must satisfy the constraint n2${\mathbf{n}}\ge 2$.
u(i,1)${\mathbf{u}}\left(\mathit{i},1\right)$ must be set to the known or estimated value of yi${y}_{\mathit{i}}$ at a$a$ and u(i,2)${\mathbf{u}}\left(\mathit{i},2\right)$ must be set to the known or estimated value of yi${y}_{\mathit{i}}$ at b$b$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$.
2:     v(n,2$2$) – double array
n, the first dimension of the array, must satisfy the constraint n2${\mathbf{n}}\ge 2$.
v(i,j)${\mathbf{v}}\left(\mathit{i},\mathit{j}\right)$ must be set to 0.0$0.0$ if u(i,j)${\mathbf{u}}\left(\mathit{i},\mathit{j}\right)$ is a known value and to 1.0$1.0$ if u(i,j)${\mathbf{u}}\left(\mathit{i},\mathit{j}\right)$ is an estimated value, for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$ and j = 1,2$\mathit{j}=1,2$.
Constraint: precisely n$\mathit{n}$ of the v(i,j)${\mathbf{v}}\left(i,j\right)$ must be set to 0.0$0.0$, i.e., precisely n$\mathit{n}$ of the u(i,j)${\mathbf{u}}\left(i,j\right)$ must be known values, and these must not be all at a$a$ or all at b$b$.
3:     a – double scalar
a$a$, the initial point of the interval of integration.
4:     b – double scalar
b$b$, the final point of the interval of integration.
5:     tol – double scalar
Must be set to a small quantity suitable for:
 (a) testing the local error in yi${y}_{i}$ during integration, (b) testing for the convergence of yi${y}_{i}$ at b$b$, (c) calculating the perturbation in estimated boundary values for yi${y}_{i}$, which are used to obtain the approximate derivatives of the residuals for use in the Newton iteration.
You are advised to check your results by varying tol.
Constraint: tol > 0.0${\mathbf{tol}}>0.0$.
6:     fcn – function handle or string containing name of m-file
fcn must evaluate the functions fi${f}_{\mathit{i}}$ (i.e., the derivatives yi${y}_{\mathit{i}}^{\prime }$), for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$, at a general point x$x$.
[f] = fcn(x, y)

Input Parameters

1:     x – double scalar
x$x$, the value of the argument.
2:     y( : $:$) – double array
yi${y}_{\mathit{i}}$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$, the value of the argument.

Output Parameters

1:     f( : $:$) – double array
The values of fi(x)${f}_{\mathit{i}}\left(x\right)$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$.
7:     m1 – int64int32nag_int scalar
A value which controls output.
m1 = 1${\mathbf{m1}}=1$
The final solution is not evaluated.
m1 > 1${\mathbf{m1}}>1$
The final values of yi${y}_{\mathit{i}}$ at interval (ba) / (m11)$\left(b-a\right)/\left({\mathbf{m1}}-1\right)$ are calculated and stored in the array soln by columns, starting with values yi${y}_{\mathit{i}}$ at a$a$ stored in soln(i,1)${\mathbf{soln}}\left(\mathit{i},1\right)$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$.
Constraint: m11${\mathbf{m1}}\ge 1$.

### Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The first dimension of the arrays u, v. (An error is raised if these dimensions are not equal.)
n$\mathit{n}$, the number of equations.
Constraint: n2${\mathbf{n}}\ge 2$.

sdw

### Output Parameters

1:     u(n,2$2$) – double array
The known values unaltered, and corrected values of the estimates, unless an error has occurred. If an error has occurred, u contains the known values and the latest values of the estimates.
2:     soln(n,m1) – double array
The solution when m1 > 1${\mathbf{m1}}>1$.
3:     w(n,sdw) – double array
sdw3n + 17 + max (11,n)$\mathit{sdw}\ge 3{\mathbf{n}}+17+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(11,{\mathbf{n}}\right)$.
If ${\mathbf{ifail}}={\mathbf{2}}$, 3${\mathbf{3}}$, 4${\mathbf{4}}$ or 5${\mathbf{5}}$, w(i,1)${\mathbf{w}}\left(\mathit{i},1\right)$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$, contains the solution at the point where the integration fails and the point of failure is returned in w(1,2)${\mathbf{w}}\left(1,2\right)$.
4:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
One or more of the parameters v, n, m1, sdw, or tol is incorrectly set.
ifail = 2${\mathbf{ifail}}=2$
The step length for the integration is too short whilst calculating the residual (see Section [Further Comments]).
ifail = 3${\mathbf{ifail}}=3$
No initial step length could be chosen for the integration whilst calculating the residual.
Note: ${\mathbf{ifail}}={\mathbf{2}}$ or 3${\mathbf{3}}$ can occur due to choosing too small a value for tol or due to choosing the wrong direction of integration. Try varying tol and interchanging a$a$ and b$b$. These error exits can also occur for very poor initial estimates of the unknown initial values and, in extreme cases, because nag_ode_bvp_shoot_bval (d02ha) cannot be used to solve the problem posed.
ifail = 4${\mathbf{ifail}}=4$
As for ${\mathbf{ifail}}={\mathbf{2}}$ but the error occurred when calculating the Jacobian of the derivatives of the residuals with respect to the parameters.
ifail = 5${\mathbf{ifail}}=5$
As for ${\mathbf{ifail}}={\mathbf{3}}$ but the error occurred when calculating the derivatives of the residuals with respect to the parameters.
ifail = 6${\mathbf{ifail}}=6$
The calculated Jacobian has an insignificant column.
Note: ${\mathbf{ifail}}={\mathbf{4}}$, 5${\mathbf{5}}$ or 6${\mathbf{6}}$ usually indicate a badly scaled problem. You may vary the size of tol or change to one of the more general functions nag_ode_bvp_shoot_genpar (d02hb) or nag_ode_bvp_shoot_genpar_algeq (d02sa) which afford more control over the calculations.
ifail = 7${\mathbf{ifail}}=7$
The linear algebra function (nag_lapack_dgesvd (f08kb)) used has failed. This error exit should not occur and can be avoided by changing the estimated initial values.
ifail = 8${\mathbf{ifail}}=8$
The Newton iteration has failed to converge.
Note: ${\mathbf{ifail}}={\mathbf{8}}$ can indicate poor initial estimates or a very difficult problem. Consider varying tol if the residuals are small in the monitoring output. If the residuals are large try varying the initial estimates.
ifail = 9${\mathbf{ifail}}=9$
ifail = 10${\mathbf{ifail}}=10$
ifail = 11${\mathbf{ifail}}=11$
ifail = 12${\mathbf{ifail}}=12$
ifail = 13${\mathbf{ifail}}=13$
Indicates that a serious error has occurred in an internal call. Check all array subscripts and function parameter lists in calls to nag_ode_bvp_shoot_bval (d02ha). Seek expert help.

## Accuracy

If the process converges, the accuracy to which the unknown parameters are determined is usually close to that specified by you; the solution, if requested, may be determined to a required accuracy by varying tol.

The time taken by nag_ode_bvp_shoot_bval (d02ha) depends on the complexity of the system, and on the number of iterations required. In practice, integration of the differential equations is by far the most costly process involved.
Wherever it occurs in the function, the error parameter tol is used in ‘mixed’ form; that is tol always occurs in expressions of the form tol × (1 + |yi|)${\mathbf{tol}}×\left(1+|{y}_{i}|\right)$. Though not ideal for every application, it is expected that this mixture of absolute and relative error testing will be adequate for most purposes.
You are strongly recommended to set ifail to obtain self-explanatory error messages, and also monitoring information about the course of the computation. You may select the channel numbers on which this output is to appear by calls of nag_file_set_unit_error (x04aa) (for error messages) or nag_file_set_unit_advisory (x04ab) (for monitoring information) – see Section [Example] for an example. Otherwise the default channel numbers will be used. The monitoring information produced at each iteration includes the current parameter values, the residuals and 2$2$-norms: a basic norm and a current norm. At each iteration the aim is to find parameter values which make the current norm less than the basic norm. Both these norms should tend to zero as should the residuals. (They would all be zero if the exact parameters were used as input.) For more details, you may consult the specification of nag_ode_bvp_shoot_genpar_algeq (d02sa), and especially the description of the parameter monit there.
The computing time for integrating the differential equations can sometimes depend critically on the quality of the initial estimates. If it seems that too much computing time is required and, in particular, if the values of the residuals printed by the monitoring function are much larger than the expected values of the solution at b$b$, then the coding of fcn should be checked for errors. If no errors can be found, an independent attempt should be made to improve the initial estimates. In practical problems it is not uncommon for the differential equation to have a singular point at one or both ends of the range. Suppose a$a$ is a singular point; then the derivatives yi${y}_{i}^{\prime }$ in (1) (in Section [Description]) cannot be evaluated at a$a$, usually because one or more of the expressions for fi${f}_{i}$ give overflow. In such a case it is necessary for you to take a$a$ a short distance away from the singularity, and to find values for yi${y}_{i}$ at the new value of a$a$ (e.g., use the first one or two terms of an analytical (power series) solution). You should experiment with the new position of a$a$; if it is taken too close to the singular point, the derivatives fi${f}_{i}$ will be inaccurate, and the function may sometimes fail with ${\mathbf{ifail}}={\mathbf{2}}$ or 3${\mathbf{3}}$ or, in extreme cases, with an overflow condition. A more general treatment of singular solutions is provided by the function nag_ode_bvp_shoot_genpar (d02hb).
Another difficulty which often arises in practice is the case when one end of the range, b$b$ say, is at infinity. You must approximate the end point by taking a finite value for b$b$, which is obtained by estimating where the solution will reach its asymptotic state. The estimate can be checked by repeating the calculation with a larger value of b$b$. If b$b$ is very large, and if the matching point is also at b$b$, the numerical solution may suffer a considerable loss of accuracy in integrating across the range, and the program may fail with ${\mathbf{ifail}}={\mathbf{6}}$ or 8${\mathbf{8}}$. (In the former case, solutions from all initial values at a$a$ are tending to the same curve at infinity.) The simplest remedy is to try to solve the equations with a smaller value of b$b$, and then to increase b$b$ in stages, using each solution to give boundary value estimates for the next calculation. For problems where some terms in the asymptotic form of the solution are known, nag_ode_bvp_shoot_genpar (d02hb) will be more successful.
If the unknown quantities are not boundary values, but are eigenvalues or the length of the range or some other parameters occurring in the differential equations, nag_ode_bvp_shoot_genpar (d02hb) may be used.

## Example

This example finds the angle at which a projectile must be fired for a given range.
The differential equations are:
 y′ = tanφ v′ = ( − 0.032 tanφ )/v − (0.02v)/(cosφ) φ ′ = ( − 0.032)/(v2) ,
$y′ = tanϕ v′ = -0.032 tanϕ v - 0.02v cosϕ ϕ′ = -0.032 v2 ,$
with the following boundary conditions:
 y = 0, v = 0.5 at   x = 0, y = 0 at   x = 5.
$y= 0, v= 0.5 at x= 0, y= 0 at x= 5.$
The remaining boundary conditions are estimated as:
 φ = 1.15 at  x = 0, φ = 1.2, v = 0.46 at  x = 5.
$ϕ=1.15 at x=0, ϕ=1.2, v=0.46 at x=5.$
We write y = Z(1)$y=\mathrm{Z}\left(1\right)$, v = Z(2)$v=\mathrm{Z}\left(2\right)$, φ = Z(3)$\varphi =\mathrm{Z}\left(3\right)$. To check the accuracy of the results the problem is solved twice with tol = 5${\mathbf{tol}}=5$.0e−3$\text{0e−3}$ and 5.0e−4$\text{5.0e−4}$ respectively. Note the call to nag_file_set_unit_advisory (x04ab) before the call to nag_ode_bvp_shoot_bval (d02ha).
```function nag_ode_bvp_shoot_bval_example
% Initialize variables and arrays.
u = [0, 0; 0.5, 0.46; 1.15, -1.2];
v = [0, 0; 0, 1; 1, 1];
a = 0;
b = 5;
m1 = 6;
n = 3;
ykeep = zeros(m1,3);
xkeep = zeros(1,m1);

fprintf('nag_ode_bvp_shoot_bval example program results\n\n');

for id = 3:4
tol = 5*10^(-id);

[uOut, soln, w, ifail] = nag_ode_bvp_shoot_bval(u, v, a, b, tol, @fcn, int64(m1));
if (ifail ~= 0)
% Problems in integration, or incorrect parameters.  Print
% message and exit.
error('Warning: nag_ode_bvp_shoot_bval returned with ifail = %d ',ifail);
else
% Output results.
fprintf('Results with tol = %1.3e\n\n',tol);
fprintf('X-value and final solution\n\n');
for i = 1:m1
disp((sprintf('%d   %8.4f   %8.4f   %8.4f', i-1,...
soln(1,i),soln(2,i),soln(3,i))));
% Store results for plotting.
ykeep(i,1) = soln(1,i);
ykeep(i,2) = soln(2,i);
ykeep(i,3) = soln(3,i);
xkeep(i) = i-1;
end
end
disp(' ');
end
% Plot results.
fig = figure('Number', 'off');
display_plot(xkeep, ykeep)

function f = fcn(x,y)
% Evaluate the derivatives.
f = zeros(3,1);
f(1) = tan(y(3));
f(2) = -0.032*tan(y(3))/y(2) - 0.02*y(2)/cos(y(3));
f(3) = -0.032/y(2)^2;
function  display_plot(xkeep, ykeep)

% Formatting for title and axis labels.

titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Plot curves.

plot(xkeep, ykeep(:,1), '-+',...
xkeep, ykeep(:,2), '--x',...
xkeep, ykeep(:,3), ':*');

title({'Solution of Two-point Boundary-value Problem',...
['using Runge-Kutta-Merson and Newton Correction in a ', ...
'Shooting Method']}, titleFmt{:});
% Label the axes.

xlabel('x', labFmt{:});
ylabel('Solution', labFmt{:});

legend('height','angle','velocity','location','Best');
```
```
nag_ode_bvp_shoot_bval example program results

Results with tol = 5.000e-03

X-value and final solution

0     0.0000     0.5000     1.1685
1     1.9197     0.3341     0.9752
2     2.9304     0.2066     0.4917
3     2.9767     0.1956    -0.4228
4     2.0258     0.3093    -0.9754
5    -0.0005     0.4598    -1.2020

Results with tol = 5.000e-04

X-value and final solution

0     0.0000     0.5000     1.1681
1     1.9177     0.3343     0.9749
2     2.9280     0.2070     0.4929
3     2.9769     0.1955    -0.4194
4     2.0210     0.3095    -0.9751
5    -0.0000     0.4597    -1.2014

```
```function d02ha_example
% Initialize variables and arrays.
u = [0, 0; 0.5, 0.46; 1.15, -1.2];
v = [0, 0; 0, 1; 1, 1];
a = 0;
b = 5;
m1 = 6;
n = 3;
ykeep = zeros(m1,3);
xkeep = zeros(1,m1);

fprintf('d02ha example program results\n\n');

for id = 3:4
tol = 5*10^(-id);

[uOut, soln, w, ifail] = d02ha(u, v, a, b, tol, @fcn, int64(m1));
if (ifail ~= 0)
% Problems in integration, or incorrect parameters.  Print
% message and exit.
error('Warning: d02ha returned with ifail = %d ',ifail);
else
% Output results.
fprintf('Results with tol = %1.3e\n\n',tol);
fprintf('X-value and final solution\n\n');
for i = 1:m1
disp((sprintf('%d   %8.4f   %8.4f   %8.4f', i-1,...
soln(1,i),soln(2,i),soln(3,i))));
% Store results for plotting.
ykeep(i,1) = soln(1,i);
ykeep(i,2) = soln(2,i);
ykeep(i,3) = soln(3,i);
xkeep(i) = i-1;
end
end
disp(' ');
end
% Plot results.
fig = figure('Number', 'off');
display_plot(xkeep, ykeep)

function f = fcn(x,y)
% Evaluate the derivatives.
f = zeros(3,1);
f(1) = tan(y(3));
f(2) = -0.032*tan(y(3))/y(2) - 0.02*y(2)/cos(y(3));
f(3) = -0.032/y(2)^2;
function  display_plot(xkeep, ykeep)

% Formatting for title and axis labels.

titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Plot curves.

plot(xkeep, ykeep(:,1), '-+',...
xkeep, ykeep(:,2), '--x',...
xkeep, ykeep(:,3), ':*');

title({'Solution of Two-point Boundary-value Problem',...
['using Runge-Kutta-Merson and Newton Correction in a ', ...
'Shooting Method']}, titleFmt{:});
% Label the axes.

xlabel('x', labFmt{:});
ylabel('Solution', labFmt{:});

legend('height','angle','velocity','location','Best');
```
```
d02ha example program results

Results with tol = 5.000e-03

X-value and final solution

0     0.0000     0.5000     1.1685
1     1.9197     0.3341     0.9752
2     2.9304     0.2066     0.4917
3     2.9767     0.1956    -0.4228
4     2.0258     0.3093    -0.9754
5    -0.0005     0.4598    -1.2020

Results with tol = 5.000e-04

X-value and final solution

0     0.0000     0.5000     1.1681
1     1.9177     0.3343     0.9749
2     2.9280     0.2070     0.4929
3     2.9769     0.1955    -0.4194
4     2.0210     0.3095    -0.9751
5    -0.0000     0.4597    -1.2014

```

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013