{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Machine Learning and the Physical World\n", "=======================================\n", "\n", "### [Neil D. Lawrence](http://inverseprobability.com), Amazon Cambridge\n", "\n", "and University of Sheffield\n", "\n", "### 2018-12-10" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Abstract**: Machine learning is a data driven endeavour, but real\n", "world systems are physical and mechanistic. In this talk we will review\n", "approaches to integrating machine learning with real world systems. Our\n", "focus will be on emulation (otherwise known as surrogate modeling)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Centrifugal Governor\n", "------------------------\n", "\n", "\\[edit\\]\n", "\n", "\n", "\n", "Figure: Centrifugal governor as held by “Science” on Holborn\n", "Viaduct" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Boulton and Watt’s Steam Engine\n", "-------------------------------\n", "\n", "\\[edit\\]\n", "\n", "\n", "\n", "Figure: Watt’s Steam Engine which made Steam Power Efficient and\n", "Practical.\n", "\n", "James Watt’s steam engine contained an early machine learning device. In\n", "the same way that modern systems are component based, his engine was\n", "composed of components. One of which is a speed regulator sometimes\n", "known as *Watt’s governor*. The two balls in the center of the image,\n", "when spun fast, rise, and through a linkage mechanism.\n", "\n", "The centrifugal governor was made famous by Boulton and Watt when it was\n", "deployed in the steam engine. Studying stability in the governor is the\n", "main subject of James Clerk Maxwell’s paper on the theoretical analysis\n", "of governors (Maxwell, 1867). This paper is a founding paper of control\n", "theory. In an acknowledgment of its influence, Wiener used the name\n", "[*cybernetics*](https://en.wikipedia.org/wiki/Cybernetics) to describe\n", "the field of control and communication in animals and the machine\n", "(Wiener, 1948). Cybernetics is the Greek word for governor, which comes\n", "from the latin for helmsman.\n", "\n", "A governor is one of the simplest artificial intelligence systems. It\n", "senses the speed of an engine, and acts to change the position of the\n", "valve on the engine to slow it down.\n", "\n", "Although it’s a mechanical system a governor can be seen as automating a\n", "role that a human would have traditionally played. It is an early\n", "example of artificial intelligence.\n", "\n", "The centrifugal governor has several parameters, the weight of the balls\n", "used, the length of the linkages and the limits on the balls movement.\n", "\n", "Two principle differences exist between the centrifugal governor and\n", "artificial intelligence systems of today.\n", "\n", "1. The centrifugal governor is a physical system and it is an integral\n", " part of a wider physical system that it regulates (the engine).\n", "2. The parameters of the governor were set by hand, our modern\n", " artificial intelligence systems have their parameters set by *data*.\n", "\n", "\n", "\n", "Figure: The centrifugal governor, an early example of a decision\n", "making system. The parameters of the governor include the lengths of the\n", "linkages (which effect how far the throttle opens in response to\n", "movement in the balls), the weight of the balls (which effects inertia)\n", "and the limits of to which the balls can rise.\n", "\n", "This has the basic components of sense and act that we expect in an\n", "intelligent system, and this system saved the need for a human operator\n", "to manually adjust the system in the case of overspeed. Overspeed has\n", "the potential to destroy an engine, so the governor operates as a safety\n", "device.\n", "\n", "The first wave of automation did bring about sabotoage as a worker’s\n", "response. But if machinery was sabotaged, for example, if the linkage\n", "between sensor (the spinning balls) and action (the valve closure) was\n", "broken, this would be obvious to the engine operator at start up time.\n", "The machine could be repaired before operation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What is Machine Learning?\n", "-------------------------\n", "\n", "\\[edit\\]\n", "\n", "Machine learning allows us to extract knowledge from data to form a\n", "prediction.\n", "\n", "$$\\text{data} + \\text{model} \\stackrel{\\text{compute}}{\\rightarrow} \\text{prediction}$$\n", "\n", "A machine learning prediction is made by combining a model with data to\n", "form the prediction. The manner in which this is done gives us the\n", "machine learning *algorithm*.\n", "\n", "Machine learning models are *mathematical models* which make weak\n", "assumptions about data, e.g. smoothness assumptions. By combining these\n", "assumptions with the data, we observe we can interpolate between data\n", "points or, occasionally, extrapolate into the future.\n", "\n", "Machine learning is a technology which strongly overlaps with the\n", "methodology of statistics. From a historical/philosophical view point,\n", "machine learning differs from statistics in that the focus in the\n", "machine learning community has been primarily on accuracy of prediction,\n", "whereas the focus in statistics is typically on the interpretability of\n", "a model and/or validating a hypothesis through data collection.\n", "\n", "The rapid increase in the availability of compute and data has led to\n", "the increased prominence of machine learning. This prominence is\n", "surfacing in two different but overlapping domains: data science and\n", "artificial intelligence." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From Model to Decision\n", "----------------------\n", "\n", "\\[edit\\]\n", "\n", "The real challenge, however, is end-to-end decision making. Taking\n", "information from the environment and using it to drive decision making\n", "to achieve goals." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Process Automation\n", "------------------\n", "\n", "\\[edit\\]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from IPython.lib.display import YouTubeVideo\n", "YouTubeVideo('3HJtmx5f1Fc')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Figure: An actual Santa’s sleigh. Amazon’s new delivery drone.\n", "Machine learning algorithms are used across various systems including\n", "sensing (computer vision for detection of wires, people, dogs etc) and\n", "piloting. The technology is necessarily a combination of old and new\n", "ideas. The transition from vertical to horizontal flight is vital for\n", "efficiency and requires sophisticated machine learning to achieve." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Artificial Intelligence and Data Science\n", "----------------------------------------\n", "\n", "\\[edit\\]\n", "\n", "Artificial intelligence has the objective of endowing computers with\n", "human-like intelligent capabilities. For example, understanding an image\n", "(computer vision) or the contents of some speech (speech recognition),\n", "the meaning of a sentence (natural language processing) or the\n", "translation of a sentence (machine translation)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Supervised Learning for AI\n", "\n", "The machine learning approach to artificial intelligence is to collect\n", "and annotate a large data set from humans. The problem is characterized\n", "by input data (e.g. a particular image) and a label (e.g. is there a car\n", "in the image yes/no). The machine learning algorithm fits a mathematical\n", "function (I call this the *prediction function*) to map from the input\n", "image to the label. The parameters of the prediction function are set by\n", "minimizing an error between the function’s predictions and the true\n", "data. This mathematical function that encapsulates this error is known\n", "as the *objective function*.\n", "\n", "This approach to machine learning is known as *supervised learning*.\n", "Various approaches to supervised learning use different prediction\n", "functions, objective functions or different optimization algorithms to\n", "fit them.\n", "\n", "For example, *deep learning* makes use of *neural networks* to form the\n", "predictions. A neural network is a particular type of mathematical\n", "function that allows the algorithm designer to introduce invariances\n", "into the function.\n", "\n", "An invariance is an important way of including prior understanding in a\n", "machine learning model. For example, in an image, a car is still a car\n", "regardless of whether it’s in the upper left or lower right corner of\n", "the image. This is known as translation invariance. A neural network\n", "encodes translation invariance in *convolutional layers*. Convolutional\n", "neural networks are widely used in image recognition tasks.\n", "\n", "An alternative structure is known as a recurrent neural network (RNN).\n", "RNNs neural networks encode temporal structure. They use auto regressive\n", "connections in their hidden layers, they can be seen as time series\n", "models which have non-linear auto-regressive basis functions. They are\n", "widely used in speech recognition and machine translation.\n", "\n", "Machine learning has been deployed in Speech Recognition (e.g. Alexa,\n", "deep neural networks, convolutional neural networks for speech\n", "recognition), in computer vision (e.g. Amazon Go, convolutional neural\n", "networks for person recognition and pose detection).\n", "\n", "The field of data science is related to AI, but philosophically\n", "different. It arises because we are increasingly creating large amounts\n", "of data through *happenstance* rather than active collection. In the\n", "modern era data is laid down by almost all our activities. The objective\n", "of data science is to extract insights from this data.\n", "\n", "Classically, in the field of statistics, data analysis proceeds by\n", "assuming that the question (or scientific hypothesis) comes before the\n", "data is created. E.g., if I want to determine the effectiveness of a\n", "particular drug, I perform a *design* for my data collection. I use\n", "foundational approaches such as randomization to account for\n", "confounders. This made a lot of sense in an era where data had to be\n", "actively collected. The reduction in cost of data collection and storage\n", "now means that many data sets are available which weren’t collected with\n", "a particular question in mind. This is a challenge because bias in the\n", "way data was acquired can corrupt the insights we derive. We can perform\n", "randomized control trials (or A/B tests) to verify our conclusions, but\n", "the opportunity is to use data science techniques to better guide our\n", "question selection or even answer a question without the expense of a\n", "full randomized control trial (referred to as A/B testing in modern\n", "internet parlance).\n", "\n", "- There is a gap between the world of data science and AI.\n", "- The mapping of the virtual onto the physical world.\n", "- E.g. Causal understanding." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Supply Chain\n", "------------\n", "\n", "\\[edit\\]\n", "\n", "\n", "\n", "Figure: Packhorse Bridge under Burbage Edge. This packhorse route\n", "climbs steeply out of Hathersage and heads towards Sheffield. Packhorses\n", "were the main route for transporting goods across the Peak District. The\n", "high cost of transport is one driver of the ‘smith’ model, where there\n", "is a local skilled person responsible for assembling or creating goods\n", "(e.g. a blacksmith). \n", "\n", "On Sunday mornings in Sheffield, I often used to run across Packhorse\n", "Bridge in Burbage valley. The bridge is part of an ancient network of\n", "trails crossing the Pennines that, before Turnpike roads arrived in the\n", "18th century, was the main way in which goods were moved. Given that the\n", "moors around Sheffield were home to sand quarries, tin mines, lead mines\n", "and the villages in the Derwent valley were known for nail and pin\n", "manufacture, this wasn’t simply movement of agricultural goods, but it\n", "was the infrastructure for industrial transport.\n", "\n", "The profession of leading the horses was known as a Jagger and leading\n", "out of the village of Hathersage is Jagger’s Lane, a trail that headed\n", "underneath Stanage Edge and into Sheffield.\n", "\n", "The movement of goods from regions of supply to areas of demand is\n", "fundamental to our society. The physical infrastructure of supply chain\n", "has evolved a great deal over the last 300 years." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Cromford\n", "--------\n", "\n", "\\[edit\\]\n", "\n", "\n", "\n", "Figure: Richard Arkwright is regarded of the founder of the modern\n", "factory system. Factories exploit distribution networks to centralize\n", "production of goods. Arkwright located his factory in Cromford due to\n", "proximity to Nottingham Weavers (his market) and availability of water\n", "power from the tributaries of the Derwent river. When he first arrived\n", "there was almost no transportation network. Over the following 200 years\n", "The Cromford Canal (1790s), a Turnpike (now the A6, 1816-18) and the\n", "High Peak Railway (now closed, 1820s) were all constructed to improve\n", "transportation access as the factory blossomed.\n", "\n", "Richard Arkwright is known as the father of the modern factory system.\n", "In 1771 he set up a [Mill](https://en.wikipedia.org/wiki/Cromford_Mill)\n", "for spinning cotton yarn in the village of Cromford, in the Derwent\n", "Valley. The Derwent valley is relatively inaccessible. Raw cotton\n", "arrived in Liverpool from the US and India. It needed to be transported\n", "on packhorse across the bridleways of the Pennines. But Cromford was a\n", "good location due to proximity to Nottingham, where weavers where\n", "consuming the finished thread, and the availability of water power from\n", "small tributaries of the Derwent river for Arkwright’s [water\n", "frames](https://en.wikipedia.org/wiki/Spinning_jenny) which automated\n", "the production of yarn from raw cotton.\n", "\n", "By 1794 the [Cromford\n", "Canal](https://en.wikipedia.org/wiki/Cromford_Canal) was opened to bring\n", "coal in to Cromford and give better transport to Nottingham. The\n", "construction of the canals was driven by the need to improve the\n", "transport infrastructure, facilitating the movement of goods across the\n", "UK. Canals, roads and railways were initially constructed by the\n", "economic need for moving goods. To improve supply chain.\n", "\n", "The A6 now does pass through Cromford, but at the time he moved there\n", "there was merely a track. The High Peak Railway was opened in 1832, it\n", "is now converted to the High Peak Trail, but it remains the highest\n", "railway built in Britain.\n", "\n", "Cooper (1991)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Containerization\n", "----------------\n", "\n", "\\[edit\\]\n", "\n", "\n", "\n", "Figure: The container is one of the major drivers of globalization,\n", "and arguably the largest agent of social change in the last 100 years.\n", "It reduces the cost of transportation, significantly changing the\n", "appropriate topology of distribution networks. The container makes it\n", "possible to ship goods halfway around the world for cheaper than it\n", "costs to process those goods, leading to an extended distribution\n", "topology.\n", "\n", "Containerization has had a dramatic effect on global economics, placing\n", "many people in the developing world at the end of the supply chain.\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "Figure: Wild Alaskan Cod, being solid in the Pacific Northwest, that\n", "is a product of China. It is cheaper to ship the deep frozen fish\n", "thousands of kilometers for processing than to process locally.\n", "\n", "For example, you can buy Wild Alaskan Cod fished from Alaska, processed\n", "in China, sold in North America. This is driven by the low cost of\n", "transport for frozen cod vs the higher relative cost of cod processing\n", "in the US versus China. Similarly,\n", "Scottish\n", "prawns are also processed in China for sale in the UK.\n", "\n", "This effect on cost of transport vs cost of processing is the main\n", "driver of the topology of the modern supply chain and the associated\n", "effect of globalization. If transport is much cheaper than processing,\n", "then processing will tend to agglomerate in places where processing\n", "costs can be minimized.\n", "\n", "Large scale global economic change has principally been driven by\n", "changes in the technology that drives supply chain.\n", "\n", "Supply chain is a large-scale automated decision making network. Our aim\n", "is to make decisions not only based on our models of customer behavior\n", "(as observed through data), but also by accounting for the structure of\n", "our fulfilment center, and delivery network.\n", "\n", "Many of the most important questions in supply chain take the form of\n", "counterfactuals. E.g. “What would happen if we opened a manufacturing\n", "facility in Cambridge?” A counter factual is a question that implies a\n", "mechanistic understanding of a system. It goes beyond simple smoothness\n", "assumptions or translation invariants. It requires a physical, or\n", "*mechanistic* understanding of the supply chain network. For this\n", "reason, the type of models we deploy in supply chain often involve\n", "simulations or more mechanistic understanding of the network.\n", "\n", "In supply chain Machine Learning alone is not enough, we need to bridge\n", "between models that contain real mechanisms and models that are entirely\n", "data driven.\n", "\n", "This is challenging, because as we introduce more mechanism to the\n", "models we use, it becomes harder to develop efficient algorithms to\n", "match those models to data.\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "UNCERTAINTY QUANTIFICATION\n", "==========================\n", "\n", "\\[edit\\]\n", "\n", "Machine learning aims to replicate processes through the direct use of\n", "data. When deployed in the domain of ‘artificial intelligence,’ the\n", "processes that it is replicating, or *emulating*, are cognitive\n", "processes.\n", "\n", "The first trick in machine learning is to convert the process itself\n", "into a *mathematical function*. That function has a set of parameters\n", "which control its behaviour. What we call learning is the adaption of\n", "these parameters to change the behavior of the function. The choice of\n", "mathematical function we use is a vital component of the model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Emukit Playground\n", "-----------------\n", "\n", "\\[edit\\]\n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", "Adam Hirst\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", "Cliff McCollum\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Emukit playground is a software toolkit for exploring the use of\n", "statistical emulation as a tool. It was built by [Adam\n", "Hirst](https://twitter.com/_AdamHirst), during his software engineering\n", "internship at Amazon and supervised by [Cliff\n", "McCollum](https://www.linkedin.com/in/cliffmccollum/).\n", "\n", "\n", "\n", "Figure: Emukit playground is a tutorial for understanding the\n", "simulation/emulation relationship.\n", "https://amzn.github.io/emukit-playground/\n", "\n", "\n", "\n", "Figure: Tutorial on Bayesian optimization of the number of taxis\n", "deployed from Emukit playground.\n", "https://amzn.github.io/emukit-playground/#!/learn/bayesian_optimization\n", "\n", "You can explore Bayesian optimization of a taxi simulation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Uncertainty Quantification\n", "--------------------------\n", "\n", "\\[edit\\]\n", "\n", "> Uncertainty quantification (UQ) is the science of quantitative\n", "> characterization and reduction of uncertainties in both computational\n", "> and real world applications. It tries to determine how likely certain\n", "> outcomes are if some aspects of the system are not exactly known.\n", "\n", "We will to illustrate different concepts of [Uncertainty\n", "Quantification](https://en.wikipedia.org/wiki/Uncertainty_quantification)\n", "(UQ) and the role that Gaussian processes play in this field. Based on a\n", "simple simulator of a car moving between a valley and a mountain, we are\n", "going to illustrate the following concepts:\n", "\n", "- **Systems emulation**. Many real world decisions are based on\n", " simulations that can be computationally very demanding. We will show\n", " how simulators can be replaced by *emulators*: Gaussian process\n", " models fitted on a few simulations that can be used to replace the\n", " *simulator*. Emulators are cheap to compute, fast to run, and always\n", " provide ways to quantify the uncertainty of how precise they are\n", " compared the original simulator.\n", "\n", "- **Emulators in optimization problems**. We will show how emulators\n", " can be used to optimize black-box functions that are expensive to\n", " evaluate. This field is also called Bayesian Optimization and has\n", " gained an increasing relevance in machine learning as emulators can\n", " be used to optimize computer simulations (and machine learning\n", " algorithms) quite efficiently.\n", "\n", "- **Multi-fidelity emulation methods**. In many scenarios we have\n", " simulators of different quality about the same measure of interest.\n", " In these cases the goal is to merge all sources of information under\n", " the same model so the final emulator is cheaper and more accurate\n", " than an emulator fitted only using data from the most accurate and\n", " expensive simulator." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Mountain Car Simulator\n", "----------------------\n", "\n", "\\[edit\\]\n", "\n", "To illustrate the above mentioned concepts we we use the [mountain car\n", "simulator](https://github.com/openai/gym/wiki/MountainCarContinuous-v0).\n", "This simulator is widely used in machine learning to test reinforcement\n", "learning algorithms. The goal is to define a control policy on a car\n", "whose objective is to climb a mountain. Graphically, the problem looks\n", "as follows:\n", "\n", "\n", "\n", "Figure: The mountain car simulation from the Open AI gym.\n", "\n", "The goal is to define a sequence of actions (push the car right or left\n", "with certain intensity) to make the car reach the flag after a number\n", "$T$ of time steps.\n", "\n", "At each time step $t$, the car is characterized by a vector\n", "$\\mathbf{ x}_{t} = (p_t,v_t)$ of states which are respectively the the\n", "position and velocity of the car at time $t$. For a sequence of states\n", "(an episode), the dynamics of the car is given by\n", "\n", "$$\n", "\\mathbf{ x}_{t+1} = f(\\mathbf{ x}_{t},\\textbf{u}_{t})\n", "$$\n", "\n", "where $\\textbf{u}_{t}$ is the value of an action force, which in this\n", "example corresponds to push car to the left (negative value) or to the\n", "right (positive value). The actions across a full episode are\n", "represented in a policy $\\textbf{u}_{t} = \\pi(\\mathbf{ x}_{t},\\theta)$\n", "that acts according to the current state of the car and some parameters\n", "$\\theta$. In the following examples we will assume that the policy is\n", "linear which allows us to write $\\pi(\\mathbf{ x}_{t},\\theta)$ as\n", "\n", "$$\n", "\\pi(\\mathbf{ x},\\theta)= \\theta_0 + \\theta_p p + \\theta_vv.\n", "$$ For $t=1,\\dots,T$ now given some initial state $\\mathbf{ x}_{0}$ and\n", "some some values of each $\\textbf{u}_{t}$, we can **simulate** the full\n", "dynamics of the car for a full episode using\n", "[Gym](https://gym.openai.com/envs/). The values of $\\textbf{u}_{t}$ are\n", "fully determined by the parameters of the linear controller.\n", "\n", "After each episode of length $T$ is complete, a reward function\n", "$R_{T}(\\theta)$ is computed. In the mountain car example the reward is\n", "computed as 100 for reaching the target of the hill on the right hand\n", "side, minus the squared sum of actions (a real negative to push to the\n", "left and a real positive to push to the right) from start to goal. Note\n", "that our reward depend on $\\theta$ as we make it dependent on the\n", "parameters of the linear controller." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Emulate the Mountain Car\n", "------------------------" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install gym" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import gym" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "env = gym.make('MountainCarContinuous-v0')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our goal in this section is to find the parameters $\\theta$ of the\n", "linear controller such that\n", "\n", "$$\n", "\\theta^* = arg \\max_{\\theta} R_T(\\theta).\n", "$$\n", "\n", "In this section, we directly use Bayesian optimization to solve this\n", "problem. We will use [EmuKit](https://emukit.github.io) so we first\n", "define the objective function." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import urllib.request" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "urllib.request.urlretrieve('https://raw.githubusercontent.com/lawrennd/talks/gh-pages/mountain_car.py','mountain_car.py')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import mountain_car as mc\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def target_function(X):\n", " \"\"\"Run the Mountain Car simulaton for each set of controller parameters in the matrix.\"\"\"\n", " simulation_function = lambda x: mc.run_simulation(env, x)[0]\n", " return np.asarray([simulation_function(np.atleast_2d(x)) for x in X])[:, np.newaxis]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For each set of parameter values of the linear controller we can run an\n", "episode of the simulator (that we fix to have a horizon of $T=500$) to\n", "generate the reward. Using as input the parameters of the controller and\n", "as outputs the rewards we can build a Gaussian process emulator of the\n", "reward.\n", "\n", "We start defining the input space, which is three-dimensional:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.core import ContinuousParameter, ParameterSpace" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "position_domain = [-1.2, +1]\n", "velocity_domain = [-1/0.07, +1/0.07]\n", "constant_domain = [-1, +1]\n", "\n", "space = ParameterSpace(\n", " [ContinuousParameter('position_parameter', *position_domain), \n", " ContinuousParameter('velocity_parameter', *velocity_domain),\n", " ContinuousParameter('constant', *constant_domain)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To initalize the model we start sampling some initial points for the\n", "linear controller randomly." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.core.initial_designs import RandomDesign" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "design = RandomDesign(space)\n", "n_initial_points = 25\n", "initial_design = design.get_samples(n_initial_points)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now run the simulation 25 times across our initial design." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "y = target_function(initial_design)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we start any optimization, lets have a look to the behaviour of\n", "the car with the first of these initial points that we have selected\n", "randomly." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "random_controller = initial_design[0,:]\n", "_, _, _, frames = mc.run_simulation(env, np.atleast_2d(random_controller), render=True)\n", "anim=mc.animate_frames(frames, 'Random linear controller')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from IPython.core.display import HTML" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "HTML(anim.to_jshtml())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mc.save_frames(frames, \n", " diagrams='./uq', \n", " filename='mountain-car-random.html')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Figure: Random linear controller for the Mountain car. It fails to\n", "move the car to the top of the mountain.\n", "\n", "As we can see the random linear controller does not manage to push the\n", "car to the top of the mountain. Now, let’s optimize the regret using\n", "Bayesian optimization and the emulator for the reward. We try 50 new\n", "parameters chosen by the expected improvement acquisition function.\n", "\n", "First we initizialize a Gaussian process emulator." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import GPy" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "kern = GPy.kern.RBF(3)\n", "model_gpy = GPy.models.GPRegression(initial_design, y, kern, noise_var=1e-10)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.model_wrappers.gpy_model_wrappers import GPyModelWrapper" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model_emukit = GPyModelWrapper(model_gpy, n_restarts=5)\n", "model_emukit.optimize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In Bayesian optimization an acquisition function is used to balance\n", "exploration and exploitation to evaluate new locations close to the\n", "optimum of the objective. In this notebook we select the expected\n", "improvement (EI). For further details have a look at the review paper of\n", "Shahriari et al. (2016)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.bayesian_optimization.acquisitions import ExpectedImprovement" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "acquisition = ExpectedImprovement(model_emukit)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.bayesian_optimization.loops.bayesian_optimization_loop import BayesianOptimizationLoop" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bo = BayesianOptimizationLoop(space, model_emukit, acquisition=acquisition)\n", "bo.run_loop(target_function, 50)\n", "results= bo.get_results()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we visualize the result for the best controller that we have found\n", "with Bayesian optimization." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "_, _, _, frames = mc.run_simulation(env, np.atleast_2d(results.minimum_location), render=True)\n", "anim=mc.animate_frames(frames, 'Best controller after 50 iterations of Bayesian optimization')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "HTML(anim.to_jshtml())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mc.save_frames(frames, \n", " diagrams='./uq', \n", " filename='mountain-car-simulated.html')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Figure: Mountain car simulator trained using Bayesian optimization\n", "and the simulator of the dynamics. Fifty iterations of Bayesian\n", "optimization are used to optimize the controler.\n", "\n", "The car can now make it to the top of the mountain! Emulating the reward\n", "function and using expected improvement acquisition helped us to find a\n", "linear controller that solves the problem." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Data Efficient Emulation\n", "------------------------\n", "\n", "\\[edit\\]\n", "\n", "In the previous section we solved the mountain car problem by directly\n", "emulating the reward but no considerations about the dynamics $$\n", "\\mathbf{ x}_{t+1} =g(\\mathbf{ x}_{t},\\textbf{u}_{t})\n", "$$ of the system were made.\n", "\n", "We ran the simulator 25 times in the initial design, and 50 times in our\n", "Bayesian optimization loop. That required us to call the dynamics\n", "simulation $500\\times 75 =37,500$ times, because each simulation of the\n", "car used 500 steps. In this section we will show how it is possible to\n", "reduce this number by building an emulator for $g(\\cdot)$ that can later\n", "be used to directly optimize the control.\n", "\n", "The inputs of the model for the dynamics are the velocity, the position\n", "and the value of the control so create this space accordingly." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import gym" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "env = gym.make('MountainCarContinuous-v0')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.core import ContinuousParameter, ParameterSpace" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "position_dynamics_domain = [-1.2, +0.6]\n", "velocity_dynamics_domain = [-0.07, +0.07]\n", "action_dynamics_domain = [-1, +1]\n", "\n", "space_dynamics = ParameterSpace(\n", " [ContinuousParameter('position_dynamics_parameter', *position_dynamics_domain), \n", " ContinuousParameter('velocity_dynamics_parameter', *velocity_dynamics_domain),\n", " ContinuousParameter('action_dynamics_parameter', *action_dynamics_domain)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we sample some input parameters and use the simulator to compute\n", "the outputs. Note that in this case we are not running the full\n", "episodes, we are just using the simulator to compute $\\mathbf{ x}_{t+1}$\n", "given $\\mathbf{ x}_{t}$ and $\\textbf{u}_{t}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.core.initial_designs import RandomDesign" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "design_dynamics = RandomDesign(space_dynamics)\n", "n_initial_points = 500\n", "initial_design_dynamics = design_dynamics.get_samples(n_initial_points)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import mountain_car as mc" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "### --- Simulation of the (normalized) outputs\n", "y_dynamics = np.zeros((initial_design_dynamics.shape[0], 2))\n", "for i in range(initial_design_dynamics.shape[0]):\n", " y_dynamics[i, :] = mc.simulation(initial_design_dynamics[i, :])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Normalize the data from the simulation\n", "y_dynamics_normalisation = np.std(y_dynamics, axis=0)\n", "y_dynamics_normalised = y_dynamics/y_dynamics_normalisation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The outputs are the velocity and the position. Our model will capture\n", "the change in position and velocity on time. That is, we will model\n", "\n", "$$\n", "\\Delta v_{t+1} = v_{t+1} - v_{t}\n", "$$\n", "\n", "$$\n", "\\Delta x_{t+1} = p_{t+1} - p_{t}\n", "$$\n", "\n", "with Gaussian processes with prior mean $v_{t}$ and $p_{t}$\n", "respectively. As a covariance function, we use `Matern52`. We need\n", "therefore two models to capture the full dynamics of the system." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import GPy" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "kern_position = GPy.kern.Matern52(3)\n", "position_model_gpy = GPy.models.GPRegression(initial_design_dynamics, y_dynamics[:, 0:1], kern_position, noise_var=1e-10)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "kern_velocity = GPy.kern.Matern52(3)\n", "velocity_model_gpy = GPy.models.GPRegression(initial_design_dynamics, y_dynamics[:, 1:2], kern_velocity, noise_var=1e-10)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.model_wrappers.gpy_model_wrappers import GPyModelWrapper" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "position_model_emukit = GPyModelWrapper(position_model_gpy, n_restarts=5)\n", "velocity_model_emukit = GPyModelWrapper(velocity_model_gpy, n_restarts=5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In general we might use much smarter strategies to design our emulation\n", "of the simulator. For example, we could use the variance of the\n", "predictive distributions of the models to collect points using\n", "uncertainty sampling, which will give us a better coverage of the space.\n", "For simplicity, we move ahead with the 500 randomly selected points.\n", "\n", "Now that we have a data set, we can update the emulators for the\n", "location and the velocity." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "position_model_emukit.optimize()\n", "velocity_model_emukit.optimize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now have a look to how the emulator and the simulator match.\n", "First, we show a contour plot of the car aceleration for each pair of\n", "can position and velocity. You can use the bar bellow to play with the\n", "values of the controler to compare the emulator and the simulator." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from IPython.html.widgets import interact" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "control = mc.plot_control(velocity_model_emukit)\n", "interact(control.plot_slices, control=(-1, 1, 0.05))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see how the emulator is doing a fairly good job approximating the\n", "simulator. On the edges, however, it struggles to captures the dynamics\n", "of the simulator.\n", "\n", "Given some input parameters of the linear controlling, how do the\n", "dynamics of the emulator and simulator match? In the following figure we\n", "show the position and velocity of the car for the 500 time steps of an\n", "episode in which the parameters of the linear controller have been fixed\n", "beforehand. The value of the input control is also shown." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# change the values of the linear controller to observe the trajectories.\n", "controller_gains = np.atleast_2d([0, .6, 1]) " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mc.emu_sim_comparison(env, controller_gains, \n", " [position_model_emukit, velocity_model_emukit], \n", " max_steps=500, diagrams='./uq')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Figure: Comparison between the mountain car simulator and the\n", "emulator.\n", "\n", "We now make explicit use of the emulator, using it to replace the\n", "simulator and optimize the linear controller. Note that in this\n", "optimization, we don’t need to query the simulator anymore as we can\n", "reproduce the full dynamics of an episode using the emulator. For\n", "illustrative purposes, in this example we fix the initial location of\n", "the car.\n", "\n", "We define the objective reward function in terms of the simulator." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "### --- Optimize control parameters with emulator\n", "car_initial_location = np.asarray([-0.58912799, 0])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def target_function_emulator(X):\n", " \"\"\"Run the Mountain Car simulation for each set of controller parameters in the matrix using the emulation.\"\"\"\n", " emulation_function = lambda x: mc.run_emulation([position_model_emukit, velocity_model_emukit], x, car_initial_location)[0]\n", " return np.asarray([emulation_function(np.atleast_2d(x)) for x in X])[:, np.newaxis]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "And as before, we use Bayesian optimization to find the best possible\n", "linear controller.\n", "\n", "\n", "\n", "And we optimize using Bayesian optimzation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.core import ContinuousParameter, ParameterSpace" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "position_domain = [-1.2, +1]\n", "velocity_domain = [-1/0.07, +1/0.07]\n", "constant_domain = [-1, +1]\n", "\n", "space = ParameterSpace(\n", " [ContinuousParameter('position_parameter', *position_domain), \n", " ContinuousParameter('velocity_parameter', *velocity_domain),\n", " ContinuousParameter('constant', *constant_domain)])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from emukit.core.initial_designs import RandomDesign" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "design = RandomDesign(space)\n", "n_initial_points = 25\n", "initial_design = design.get_samples(n_initial_points)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "n_initial_points = 25 random_design = RandomDesign(design_space)\n", "initial_design = random_design.get_samples(n_initial_points) acquisition\n", "= GPyOpt.acquisitions.AcquisitionEI(model, design_space,\n", "optimizer=aquisition_optimizer) evaluator =\n", "GPyOpt.core.evaluators.Sequential(acquisition)}" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bo_multifidelity = GPyOpt.methods.ModularBayesianOptimization(model, design_space, objective_multifidelity, acquisition, evaluator, initial_design)\n", "bo_multifidelity.run_optimization(max_iter=50)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "_, _, _, frames = mc.run_simulation(env, np.atleast_2d(bo_multifidelity.x_opt), render=True)\n", "anim=mc.animate_frames(frames, 'Best controller with multi-fidelity emulator')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from IPython.core.display import HTML" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "HTML(anim.to_jshtml())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mc.save_frames(frames, \n", " diagrams='./uq', \n", " filename='mountain-car-multi-fidelity.html')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Best Controller with Multi-Fidelity Emulator\n", "\n", "\n", "\n", "Figure: Mountain car learnt with multi-fidelity model. Here 250\n", "observations of the high fidelity simulator and 250 observations of the\n", "low fidelity simulator are used to learn the controller.\n", "\n", "And problem solved! We see how the problem is also solved with 250\n", "observations of the high fidelity simulator and 250 of the low fidelity\n", "simulator." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Emukit\n", "======\n", "\n", "\\[edit\\]\n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", "Javier Gonzalez\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "The Emukit software we will be using across the next part of this module\n", "is a python software library that facilitates emulation of systems. The\n", "software’s origins go back to work done by [Javier\n", "Gonzalez](https://javiergonzalezh.github.io/) as part of his\n", "post-doctoral project at the University of Sheffield. Javier led the\n", "design and build of a Bayesian optimization software. The package\n", "`GPyOpt` worked with the SheffieldML software GPy for performing\n", "Bayesian optimization.\n", "\n", "`GPyOpt` has a modular design that allows the user to provide their own\n", "surrogate models, the package is build with `GPy` as a surrogate model\n", "in mind, but other surrogate models can also be wrapped and integrated.\n", "\n", "However, `GPyOpt` doesn’t allow the full flexibility of surrogate\n", "modelling for domains like experimental design, sensitivity analysis\n", "etc.\n", "\n", "Emukit was designed and built for a more general approach. The software\n", "is MIT licensed and its design and implementation was led by Javier\n", "Gonzalez and [Andrei Paleyes](https://www.linkedin.com/in/andreipaleyes)\n", "at Amazon. Building on the experience of `GPyOpt`, the aim with Emukit\n", "was to use the modularisation ideas embedded in `GPyOpt`, but to extend\n", "them beyond the modularisation of the surrogate models to modularisation\n", "of the acquisition function.\n", "\n", "\n", "\n", "Figure: The Emukit software is a set of software tools for emulation\n", "and surrogate modeling.\n", "https://emukit.github.io/emukit/" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install gpy" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install pyDOE" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install emukit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", " \n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", "Javier Gonzalez\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", "Andrei Paleyes\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", "Mark Pullin\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", "Maren Mahsereci\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "The software was initially built by the team in Amazon. As well as\n", "Javier Gonzalez (ML side) and Andrei Paleyes (Software Engineering)\n", "included Mark Pullin, Maren Mahsereci, Alex Gessner, Aaron Klein, Henry\n", "Moss, David-Elias Künstle as well as management input from Cliff\n", "McCollum and myself." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "MXFusion: Modular Probabilistic Programming on MXNet\n", "----------------------------------------------------\n", "\n", "\\[edit\\]\n", "\n", "One challenge for practitioners in Gaussian processes, is flexible\n", "software that allows the construction of the relevant GP modle. With\n", "this in mind, the Amazon Cambridge team has developed MXFusion. It is a\n", "modular probabilistic programming language focussed on efficient\n", "implementation of hybrid GP-neural network models, but with additional\n", "probabilistic programming capabilities.\n", "\n", "\n", "\n", "Figure: MXFusion is a probabilistic programming language targeted\n", "specifically at Gaussian process models and combining them with\n", "probaiblistic neural network. It is available through the MIT license\n", "and we welcome contributions throguh the Github repository\n", "https://github.com/amzn/MXFusion.\n", "\n", "We developed the framework for greater ease of transitioning models from\n", "‘science’ to ‘production,’ our aim was to have code that could be\n", "created by scientists, but deployed in our systems through solutions\n", "such as AWS SageMaker.\n", "\n", "\\\\ericMeissner{15%}\\\\zhenwenDai{15%}\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "- Work by Eric Meissner and Zhenwen Dai.\n", "- Probabilistic programming.\n", "- Available on [Github](https://github.com/amzn/mxfusion)\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "Figure: The MXFusion software." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "MxFusion\n", "--------" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Why another framework?\n", "----------------------" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Key Requirements\n", "----------------\n", "\n", "Specialized inference methods + models, without requiring users to\n", "reimplement nor understand them every time. Leverage expert knowledge.\n", "Efficient inference, flexible framework. Existing frameworks either did\n", "one or the other: flexible, or efficient." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What does it look like?\n", "-----------------------\n", "\n", "**Modelling**\n", "\n", "**Inference**\n", "\n", "``` python\n", "m = Model()\n", "m.mu = Variable()\n", "m.s = Variable(transformation=PositiveTransformation())\n", "m.Y = Normal.define_variable(mean=m.mu, variance=m.s)\n", "```\n", "\n", "- Variable\n", "\n", "- Distribution\n", "\n", "- Function\n", "\n", "- `log_pdf`\n", "\n", "- `draw_samples`\n", "\n", "- Variational Inference\n", "\n", "- MCMC Sampling (*soon*) Built on MXNet Gluon (imperative code, not\n", " static graph)\n", "\n", "``` python\n", "infr = GradBasedInference(inference_algorithm=MAP(model=m, observed=[m.Y]))\n", "infr.run(Y=data)\n", "```\n", "\n", "- Model + Inference together form building blocks.\n", " - Just doing modular modeling with universal inference doesn’t\n", " really scale, need specialized inference methods for specialized\n", " modelling objects like non-parametrics." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "PILCO: A Model-based Policy Search\n", "----------------------------------\n", "\n", "\\[edit\\]\n", "\n", "Common reinforcement learning methods suffer from data inefficiency,\n", "which can be a issue in real world applications where gathering\n", "sufficiently large amounts of data pose economic issues and may be\n", "impossible. propose a model-based policy search method known as PILCO in\n", "part to address this issue. PILCO uses a Gaussian process (GP) for\n", "learning the dynamics of the environment and optimizes a parametric\n", "policy function using the learned dynamics model.\n", "\n", "We construct an implementation of PILCO using MXFusion. This\n", "implementation follows the main idea of PILCO and has a few enhancements\n", "in addition to the published implementation. The enhancements are as\n", "follows: \\* **Use Monte Carlo integration instead of moment\n", "estimation.** We approximate the expected reward using Monte Carlo\n", "integration instead of the proposed moment estimation approach. This\n", "removes the bias in the expected reward computation and enables a wide\n", "range of choices of kernels and policy functions. In the original work,\n", "only RBF and linear kernel and only linear and RBF network policy can be\n", "used. \\* **Use automatic differentiation.** Thanks to automatic\n", "differentiation, no gradient derivation is needed. \\* **A unified\n", "interface of Gaussian process.** MXFusion provides an unified inferface\n", "of GP modules. We allows us to easily switch among plan GP, variational\n", "sparse GP and stocastic variational GP implementations.\n", "\n", "This notebook depends on MXNet, MXFusion and Open AI Gym. These packages\n", "can be installed into your Python environment by running the following\n", "commands.\n", "\n", "Set the global configuration." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "import gym\n", "import mxnet\n", "import mxfusion\n", "os.environ['MXNET_ENGINE_TYPE'] = 'NaiveEngine'\n", "from mxfusion.common import config\n", "config.DEFAULT_DTYPE = 'float64'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We use the inverted pendulum swingup problem as an example. We use the\n", "[Pendulum-v0](https://gym.openai.com/envs/Pendulum-v0/) environment in\n", "Open AI Gym. The task is to swing the pendulum up and balance it at the\n", "inverted position. This is a classical control problem and is known to\n", "be unsolvable with a linear controller.\n", "\n", "To solve this problem with PILCO, we need three components:\n", "\n", "- Execute a policy in an real environment (an Open AI Gym simulator in\n", " this example) and collect data.\n", "- Fit a GP model as the model for the dynamics of the environment.\n", "- Optimize the policy given the dynamics model learned from all the\n", " data that have been collected so far.\n", "\n", "The overall PILCO algorithm is to iterate the above three steps until a\n", "policy that can solve the problem is found." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Execute the Environment\n", "-----------------------" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import gym\n", "env = gym.make('Pendulum-v0')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The state of the pendulum environment is a 3D vector. The first two\n", "dimensions are the 2D location of the end point of the pendulum. The\n", "third dimension encodes the angular speed of the pendulum. The action\n", "space is a 1D vector in \\[-2, 2\\].\n", "\n", "We write a helper function for executing the environment with a given\n", "policy." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "from matplotlib import animation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def run_one_episode(env, policy, initial_state=None, max_steps=200, verbose=False, render=False):\n", " \"\"\"\n", " Drives an episode of the OpenAI gym environment using the policy to decide next actions.\n", " \"\"\"\n", " observation = env.reset()\n", " if initial_state is not None:\n", " env.env.state = initial_state\n", " observation = env.env._get_obs()\n", " env._max_episode_steps = max_steps\n", " step_idx = 0\n", " done = False\n", " total_reward = 0\n", " frames = []\n", " all_actions = []\n", " all_observations = [observation]\n", " while not done:\n", " if render:\n", " frames.append(env.render(mode = 'rgb_array'))\n", " if verbose:\n", " print(observation)\n", " action = policy(observation)\n", " observation, reward, done, info = env.step(action)\n", " all_observations.append(observation)\n", " all_actions.append(action)\n", " total_reward += reward\n", " step_idx += 1\n", " if done or step_idx>=max_steps-1:\n", " print(\"Episode finished after {} timesteps because {}\".format(step_idx+1, \"'done' reached\" if done else \"Max timesteps reached\"))\n", " break\n", " if render:\n", " fig = plt.figure()\n", " ax = fig.gca()\n", " fig.tight_layout()\n", " patch = ax.imshow(frames[0])\n", " ax.axis('off')\n", " def animate(i):\n", " patch.set_data(frames[i])\n", " anim = animation.FuncAnimation(plt.gcf(), animate, frames = len(frames), interval=20)\n", " return total_reward, np.array(all_observations, dtype=np.float64,), np.array(all_actions, dtype=np.float64), anim\n", " else:\n", " return total_reward, np.array(all_observations, dtype=np.float64,), np.array(all_actions, dtype=np.float64)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We first apply a random policy and see how the environment reacts. The\n", "random policy uniformly samples in the space of action." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def random_policy(state):\n", " return env.action_space.sample()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The animation is generated with the following commands:\n", "\n", "``` python\n", "anim = run_one_episode(env, random_policy, max_steps=500, render=True, verbose=False)[-1]\n", "\n", "with open('animation_random_policy.html', 'w') as f:\n", " f.write(anim.to_jshtml())\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pendulum\n", "--------\n", "\n", "\n", "\n", "Figure: Random policy for the control of the pendulum.\n", "\n", "The dynamics model of pendulum can be written as $$\n", "p(y_{t+1}|y_t, a_t)\n", "$$ where $y_t$ is the state vector at the time $t$ and $a_t$ is the\n", "action taken at the time $t$.\n", "\n", "PILCO uses a Gaussian process to model the above dynamics." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def prepare_data(state_list, action_list, win_in):\n", " \"\"\"\n", " Prepares a list of states and a list of actions as inputs to the Gaussian Process for training.\n", " \"\"\"\n", " \n", " X_list = []\n", " Y_list = []\n", " \n", " for state_array, action_array in zip(state_list, action_list):\n", " # the state and action array shape should be aligned.\n", " assert state_array.shape[0]-1 == action_array.shape[0]\n", " \n", " for i in range(state_array.shape[0]-win_in):\n", " Y_list.append(state_array[i+win_in:i+win_in+1])\n", " X_list.append(np.hstack([state_array[i:i+win_in].flatten(), action_array[i:i+win_in].flatten()]))\n", " X = np.vstack(X_list)\n", " Y = np.vstack(Y_list)\n", " return X, Y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example, we do a maximum likelihood estimate for the model\n", "hyper- parameters. In `MXFusion`, Gaussian process regression model is\n", "available as a module, which includes a dediated inference algorithm." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import mxnet as mx\n", "from mxfusion import Model, Variable\n", "from mxfusion.components.variables import PositiveTransformation\n", "from mxfusion.components.distributions.gp.kernels import RBF\n", "from mxfusion.modules.gp_modules import GPRegression\n", "from mxfusion.inference import GradBasedInference, MAP" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define and fit the model\n", "------------------------" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def fit_model(state_list, action_list, win_in, verbose=True):\n", " \"\"\"\n", " Fits a Gaussian Process model to the state / action pairs passed in. \n", " This creates a model of the environment which is used during\n", " policy optimization instead of querying the environment directly.\n", " \n", " See mxfusion.gp_modules for additional types of GP models to fit,\n", " including Sparse GP and Stochastic Varitional Inference Sparse GP.\n", " \"\"\"\n", " X, Y = prepare_data(state_list, action_list, win_in)\n", "\n", " m = Model()\n", " m.N = Variable()\n", " m.X = Variable(shape=(m.N, X.shape[-1]))\n", " m.noise_var = Variable(shape=(1,), transformation=PositiveTransformation(),\n", " initial_value=0.01)\n", " m.kernel = RBF(input_dim=X.shape[-1], variance=1, lengthscale=1, ARD=True)\n", " m.Y = GPRegression.define_variable(\n", " X=m.X, kernel=m.kernel, noise_var=m.noise_var,\n", " shape=(m.N, Y.shape[-1]))\n", " m.Y.factor.gp_log_pdf.jitter = 1e-6\n", "\n", " infr = GradBasedInference(\n", " inference_algorithm=MAP(model=m, observed=[m.X, m.Y]))\n", " infr.run(X=mx.nd.array(X, dtype='float64'),\n", " Y=mx.nd.array(Y, dtype='float64'),\n", " max_iter=1000, learning_rate=0.1, verbose=verbose)\n", " return m, infr, X, Y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Policy\n", "------\n", "\n", "PILCO computes the expected reward of a policy given the dynamics model.\n", "First, we need to define the parametric form of the policy. In this\n", "example, we use a neural network with one hidden layer. As the action\n", "space is \\[-2, 2\\], we apply a `tanh` transformation and multiply the\n", "come with two. This enforces the returned actions stay within the range.\n", "\n", "We define a neural network with one hidden layer and and output\n", "constrained between \\[-2,2\\] for the policy." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from mxnet.gluon import HybridBlock\n", "from mxnet.gluon.nn import Dense\n", "\n", "class NNController(HybridBlock):\n", " \"\"\"Define a neural network policy network.\"\"\"\n", " def __init__(self, prefix=None, params=None):\n", " super(NNController, self).__init__(prefix=prefix, params=params)\n", " self.dense1 = Dense(100, in_units=len(env.observation_space.high), dtype='float64', activation='relu')\n", " self.dense2 = Dense(1, in_units=100, dtype='float64', activation='tanh')\n", "\n", " def hybrid_forward(self, F, x):\n", " out = self.dense2(self.dense1(x))*2 # Scale up the output\n", " return out \n", " \n", "policy = NNController()\n", "policy.collect_params().initialize(mx.initializer.Xavier(magnitude=1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To compute the expected reward, we also need to define a reward\n", "function. This reward function is defined by us according to the task.\n", "The main component is the height of the pendulum. We also penalize the\n", "force and the angular momentum." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "class CostFunction(mx.gluon.HybridBlock):\n", " \"\"\"\n", " The goal is to get the pendulum upright and stable as quickly as possible.\n", "\n", " Taken from the code for Pendulum.\n", " \"\"\"\n", " def hybrid_forward(self, F, state, action):\n", " \"\"\"\n", " :param state: [np.cos(theta), np.sin(theta), ~ momentum(theta)]\n", " a -> 0 when pendulum is upright, largest when pendulum is hanging down completely.\n", " b -> penalty for taking action\n", " c -> penalty for pendulum momentum\n", " \"\"\"\n", " a_scale = 2.\n", " b_scale = .001\n", " c_scale = .1\n", " a = F.sum(a_scale * (state[:,:,0:1] -1) ** 2, axis=-1)\n", " b = F.sum(b_scale * action ** 2, axis=-1)\n", " c = F.sum(c_scale * state[:,:,2:3] ** 2, axis=-1)\n", " return (a + c + b)\n", " \n", "cost = CostFunction()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The expected reward function can be written as $$\n", "R = \\mathbb{E}_{p(y_T, \\ldots,\n", "y_0)}\\left(\\sum_{t=0}^\\top r(y_t)\\right)\n", "$$ where $r(\\cdot)$ is the reward function, $p(y_T, \\ldots, y_0)$ is the\n", "joint distribution when applying the policy to the dynamics model: $$\n", "p(y_T, \\ldots, y_0) = p(y_0) \\prod_{t=1}^\\top p(y_t|y_{t-1}, a_{t-1}),\n", "$$ where $a_{t-1} = \\pi(y_{t-1})$ is the action taken at the time $t-1$,\n", "which is the outcome of the policy $\\pi(\\cdot)$.\n", "\n", "The expected reward function is implemented as follows." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obtaining the policy gradients\n", "------------------------------" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from mxfusion.inference.inference_alg import SamplingAlgorithm" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "class PolicyUpdateGPParametricApprox(SamplingAlgorithm):\n", " \"\"\"Class for the policy update for PILCO.\"\"\"\n", " def compute(self, F, variables):\n", " \n", " s_0 = self.initial_state_generator(self.num_samples)\n", " a_0 = self.policy(s_0)\n", " a_t_plus_1 = a_0\n", " x_t = F.expand_dims(F.concat(s_0, a_0, dim=1), axis=1)\n", "\n", " gp = self.model.Y.factor\n", " sample_func = gp.draw_parametric_samples(F, variables, self.num_samples, self.approx_samples)\n", " cost = 0\n", " for t in range(self.n_time_steps):\n", " s_t_plus_1 = sample_func(F, x_t)\n", " cost = cost + self.cost_function(s_t_plus_1, a_t_plus_1)\n", " a_t_plus_1 = mx.nd.expand_dims(self.policy(s_t_plus_1), axis=2)\n", " x_t = mx.nd.concat(s_t_plus_1, a_t_plus_1, dim=2)\n", " total_cost = F.mean(cost)\n", " return total_cost, total_cost" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We optimize the policy with respect to the expected reward by using a\n", "gradient optimizer." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from mxfusion.inference import GradTransferInference" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def optimize_policy(policy, cost_func, model, infr, model_data_X, model_data_Y,\n", " initial_state_generator, num_grad_steps,\n", " learning_rate=1e-2, num_time_steps=100, \n", " num_samples=10, verbose=True):\n", " \"\"\"\n", " Takes as primary inputs a policy, cost function, and trained model.\n", " Optimizes the policy for num_grad_steps number of iterations.\n", " \"\"\"\n", " mb_alg = PolicyUpdateGPParametricApprox(\n", " model=model, observed=[model.X, model.Y], cost_function=cost_func,\n", " policy=policy, n_time_steps=num_time_steps,\n", " initial_state_generator=initial_state_generator,\n", " num_samples=num_samples)\n", "\n", " infr_pred = GradTransferInference(\n", " mb_alg, infr_params=infr.params, train_params=policy.collect_params())\n", " infr_pred.run(\n", " max_iter=num_grad_steps,\n", " X=mx.nd.array(model_data_X, dtype='float64'),\n", " Y=mx.nd.array(model_data_Y, dtype='float64'),\n", " verbose=verbose, learning_rate=learning_rate)\n", " return policy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Loop\n", "--------\n", "\n", "We need to define a function that provides random initial states." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def initial_state_generator(num_initial_states):\n", " \"\"\"\n", " Starts from valid states by drawing theta and momentum\n", " then computing np.cos(theta) and np.sin(theta) for state[0:2].s\n", " \"\"\"\n", " return mx.nd.array(\n", " [env.observation_space.sample() for i in range(num_initial_states)],\n", " dtype='float64')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "num_episode = 20 # how many model fit + policy optimization episodes to run\n", "num_samples = 100 # how many sample trajectories the policy optimization loop uses\n", "num_grad_steps = 1000 # how many gradient steps the optimizer takes per episode\n", "num_time_steps = 100 # how far to roll out each sample trajectory\n", "learning_rate = 1e-3 # learning rate for the policy optimization\n", "\n", "all_states = []\n", "all_actions = []" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for i_ep in range(num_episode):\n", " # Run an episode and collect data.\n", " if i_ep == 0:\n", " policy_func = lambda x: env.action_space.sample()\n", " else:\n", " policy_func = lambda x: policy(mx.nd.expand_dims(mx.nd.array(x, dtype='float64'), axis=0)).asnumpy()[0]\n", " total_reward, states, actions = run_one_episode(\n", " env, policy_func, max_steps=num_time_steps)\n", " all_states.append(states)\n", " all_actions.append(actions)\n", "\n", " # Fit a model.\n", " model, infr, model_data_X, model_data_Y = fit_model(\n", " all_states, all_actions, win_in=1, verbose=True)\n", "\n", " # Optimize the policy.\n", " policy = optimize_policy(\n", " policy, cost, model, infr, model_data_X, model_data_Y,\n", " initial_state_generator, num_grad_steps=num_grad_steps,\n", " num_samples=num_samples, learning_rate=learning_rate,\n", " num_time_steps=num_time_steps)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Policy after the first episode (random exploration):\n", "\n", "\n", "\n", "Figure: PILCO policy for control of the animation after first episode\n", "(using random exploration).\n", "\n", "Policy after the 5th episode:\n", "\n", "\n", "\n", "Figure: PILCO policy for control of the animation after the fifth\n", "episode.\n", "\n", "https://github.com/amzn/mxfusion" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Long term Aim\n", "-------------\n", "\n", "- Simulate/Emulate the components of the system.\n", " - Validate with real world using multifidelity.\n", " - Interpret system using e.g. sensitivity analysis.\n", "- Perform end to end learning to optimize.\n", " - Maintain interpretability." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Thanks!\n", "-------\n", "\n", "For more information on these subjects and more you might want to check\n", "the following resources.\n", "\n", "- twitter: [@lawrennd](https://twitter.com/lawrennd)\n", "- podcast: [The Talking Machines](http://thetalkingmachines.com)\n", "- newspaper: [Guardian Profile\n", " Page](http://www.theguardian.com/profile/neil-lawrence)\n", "- blog:\n", " [http://inverseprobability.com](http://inverseprobability.com/blog.html)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "References\n", "----------" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Cooper, B., 1991. Transformation of a valley: Derbyshire derwent.\n", "Scarthin Books.\n", "\n", "Maxwell, J.C., 1867. On governors. Proceedings of the Royal Society of\n", "London 16, 270–283.\n", "\n", "Perdikaris, P., Raissi, M., Damianou, A., Lawrence, N.D., Karniadakis,\n", "G.E., 2017. Nonlinear information fusion algorithms for data-efficient\n", "multi-fidelity modelling. Proceedings of the Royal Society of London A:\n", "Mathematical, Physical and Engineering Sciences 473.\n", "\n", "\n", "Shahriari, B., Swersky, K., Wang, Z., Adams, R.P., de Freitas, N., 2016.\n", "Taking the human out of the loop: A review of Bayesian optimization.\n", "Proceedings of the IEEE 104, 148–175.\n", "\n", "\n", "Wiener, N., 1948. Cybernetics: Control and communication in the animal\n", "and the machine. MIT Press, Cambridge, MA." ] } ], "nbformat": 4, "nbformat_minor": 5, "metadata": {} }