## Kubernetes services and load balancers

In my previous post, we have seen how we can use Kubernetes deployment objects to bring up a given number of pods running our Docker images in a cluster. However, most of the time, a pod by itself will not be able to operate – we need to connect it with other pods and the rest of the world, in other words we need to think about networking in Kubernetes.

To try this out, let us assume that you have a cluster up and running and that you have submitted a deployment of two httpd instances in the cluster. To easily get to this point, you can use the scripts in my GitHub repository as follows. These scripts will also open two ssh connections, one to each of the EC2 instances which are part of the cluster (this assumes that you are using a PEM file called eksNodeKey.pem as I have done it in my examples, if not you will have to adjust the script up.sh accordingly).

# Clone repository
$git clone https://github.com/christianb93/Kubernetes.git$ cd Kubernetes/cluster
# Bring up cluster and two EC2 nodes
$chmod 700 up.sh$ ./up.sh
# Bring down nginx controller
$kubectl delete svc ingress-nginx -n ingress-nginx # Deploy two instances of the httpd$ kubectl apply -f ../pods/deployment.yaml


Be patient, the creation of the cluster will take roughly 15 minutes, so time to get a cup of coffee. Note that we delete an object – the nginx ingress controller service – that my scripts generate and that we will use in a later post, but which blur the picture for today.

Now let us inspect our cluster and try out a few things. First, let us get the running pods.

$kubectl get pods -o wide  This should give you two pods, both running an instance of the httpd. Typically, Kubernetes will distribute these two pods over two different nodes in the cluster. Each pod has an IP address called the pod IP address which is displayed in the column IP of the output. This is the IP address under which the pod is reachable from other pods in the cluster. To verify this, let us attach to one of the pods and run curl to access the httpd running in the other pod. In my case, the first pod has IP address 192.168.232.210. We will attach to the second pod and verify that this address is reachable from there. To get a shell in the pod, we can use the kubectl exec command which executes code in a pod. The following commands extract the id of the second pod, opens a shell in this pod, installs curl and talks to the httpd in the first pod. Of course, you will have to replace the IP address of the first pod – 192.168.232.210 – with whatever your output gives you. $ name=$(kubectl get pods --output json | \ jq -r ".items[1].metadata.name")$ kubectl exec -it $name "/bin/bash" bash-4.4# apk add curl bash-4.4# curl 192.168.232.210 <h1>It works!</h1>  Nice. So we can reach a port on pod A from any other pod in the cluster. This is what is called the flat networking model in Kubernetes. In this model, each pod has a separate IP address. All containers in the pod share this IP address and one IP namespace. Every pod IP address is reachable from any other pod in the cluster without a need to set up a gateway or NAT. Coming back to the comparison of a pod with a logical host, this corresponds to a topology where all logical hosts are connected to the same IP network. In addition, Kubernetes assumes that every node can reach every pod as well. You can easily confirm this – if you log into the node (not the pod!) on which the second pod is running and use curl from there, directed to the IP address of the first pod, you will get the same result. Now Kubernetes is designed to run on a variety of platforms – locally, on a bare metal cluster, on GCP, AWS or Azure and so forth. Therefore, Kubernetes itself does not implement those rules, but leaves that to the underlying platform. To make this work, Kubernetes uses an interface called CNI (container networking interface) to talk to a plugin that is responsible for executing the platform specific part of the network configuration. On EKS, the AWS CNI plugin is used. We will get into the details in a later post, but for the time being simply assume that this works (and it does). So now we can reach every pod from every other pod. This is nice, but there is an issue – a pod is volatile. An application which is composed of several microservices needs to be prepared for the event that a pod goes down and is brought up again, be it due to a failure or simply due to the fact that an auto-scaler tries to empty a node. If, however, a pod is restarted, it will typically receive a different IP address. Suppose, for instance, you had a REST service that you want to expose within your cluster. You use a deployment to start three pods running the REST service, but which IP address should another service in the cluster use to access it? You cannot rely on the IP address of individual pods to be stable. What we need is a stable IP address which is reachable by all pods and which routes traffic to one instance of this REST service – a bit like a cluster-internal load balancer. Fortunately, Kubernetes services come to the rescue. A service is a Kubernetes object that has a stable IP address and is associated with a set of pods. When traffic is received which is directed to the service IP address, the service will select one of the pods (at random) and forward the traffic to it. Thus a service is a decoupling layer between the instable pod IP addresses and the rst of the cluster or the outer world. We will later see that behind the scenes, a service is not a process running somewhere, but essentially a set of smart NAT and routing rules. This sounds a bit abstract, so let us bring up a service. As usual, a service object is described in a manifest file. apiVersion: v1 kind: Service metadata: name: alpine-service spec: selector: app: alpine ports: - protocol: TCP port: 8080 targetPort: 80  Let us call this file service.yaml (if you have cloned my repository, you already have a copy of this in the directory network). As usual, this manifest file has a header part and a specification section. The specification section contains a selector and a list of ports. Let us look at each of those in turn. The selector plays a role similar to the selector in a deployment. It defines a set of pods that are assumed to be reachable via the service. In our case, we use the same selector as in the deployment, so that our service will send traffic to all pods brought up by this deployment. The ports section defines the ports on which the service is listening and the ports to which traffic is routed. In our case, the service will listen for TCP traffic on port 8080, and will forward this traffic to port 80 on the associated pods (as specified by the selector). We could omit the targetPort field, in this case the target port would be equal to the port. We could also specify more than one combination of port and target port and we could use names instead of numbers for the ports – refer to the documentation for a full description of all configuration options. Let us try this. Let us apply this manifest file and use kubectl get svc to list all known services. $ kubectl apply -f service.yaml
$kubectl get svc  You should now see a new service in the output, called alpine-service. Similar to a pod, this service has a cluster IP address assigned to it, and a port number (8080). In my case, this cluster IP address is 10.100.234.120. We can now again get a shell in one of the pods and try to curl that address $ name=$(kubectl get pods --output json | \ jq -r ".items[1].metadata.name")$ kubectl exec -it $name "/bin/bash" bash-4.4# apk add curl # might not be needed bash-4.4# curl 10.100.234.120:8080 <h1>It works!</h1>  If you are lucky and your container did not go down in the meantime, curl will already be installed and you can skip the apk add curl. So this works, we can actually reach the service from within our cluster. Note that we now have to use port 8080, as our service is listening on this port, not on port 80. You might ask how we can get the IP address of the service in real world? Well, fortunately Kubernetes does a bit more – it adds a DNS record for the service! So within the pod, the following will work bash-4.4# curl alpine-service:8080 <h1>It works!</h1>  So once you have the service name, you can reach the httpd from every pod within the cluster. ## Connecting to services using port forwarding Having a service which is reachable within a cluster is nice, but what options do we have to reach a cluster from the outside world? For a quick and dirty test, maybe the easiest way is using the kubectl port forwarding feature. This command allows you to forward traffic from a local port on your development machine to a port in the cluster which can be a service, but also a pod. In our case, let us forward traffic from the local port 5000 to port 8080 of our service (which is hooked up to port 80 on our pods). $ kubectl port-forward service/alpine-service 5000:8080


This will start an instance of kubectl which will bind to port 5000 on your local machine (127.0.0.1). You can now connect to this port, and behind the scenes, kubectl will tunnel the traffic through the Kubernetes master node into the cluster (you need to run this in a second terminal, as the kubectl process just started is still blocking your terminal).

$curl localhost:5000 <h1>It works!</h1>  A similar forwarding can be realized using kubectl proxy, which is designed to give you access to the Kubernetes API from your local machine, but can also be used to access services. ## Connect to a service using host ports Forwarding is easy and a quick solution, but most likely not what you want to do in a production setup. What other options do we have? One approach is to use host ports. Essentially, a host port is a port on a node that Kubernetes will wire up with the cluster IP address of a service. Assuming that you can reach the host from the outside world, you can then use the public IP address of the host to connect to a service. To create a host port, we have to modify our manifest file slightly by adding a host port specification. apiVersion: v1 kind: Service metadata: name: alpine-service spec: selector: app: alpine ports: - protocol: TCP port: 8080 targetPort: 80 type: NodePort  Note the additional line at the bottom of the file which instructs Kubernetes to open a node port. Assuming that this file is called nodePortService.yaml, we can again use kubectl to bring down our existing service and add the node port service. $ kubectl delete svc alpine-service
$kubectl apply -f nodePortService.yaml$ kubectl get svc


We see that Kubernetes has brought up our service, but this time, we see two ports in the line describing our service. The second port (32226 in my case) is the port that Kubernetes has opened on each node. Traffic to this port will be forwarded to the service IP address and port. To try this out, you can use the following commands to get the external IP address of the first node, adapt the AWS security group such that traffic to this node is allowed from your workstation, determine the node port and curl it. If your cluster is not called myCluster, replace every occurrence of myCluster with the name of your cluster.

$nodePort=$(kubectl get svc alpine-service --output json | jq ".spec.ports[0].nodePort")
$IP=$(aws ec2 describe-instances --filters Name=tag-key,Values=kubernetes.io/cluster/myCluster Name=instance-state-name,Values=running --output text --query Reservations[0].Instances[0].PublicIpAddress)
$SEC_GROUP_ID=$(aws ec2 describe-instances --filters Name=tag-key,Values=kubernetes.io/cluster/myCluster Name=instance-state-name,Values=running --output text --query Reservations[0].Instances[0].SecurityGroups[0].GroupId)
$myIP=$(wget -q -O- https://ipecho.net/plain)
$aws ec2 authorize-security-group-ingress --group-id$SEC_GROUP_ID --port $nodePort --protocol tcp --cidr "$myIP/32"
$curl$IP:$nodePort <h1>It works!</h1>  ## Connecting to a service using a load balancer A node port will allow you to connect to a service using the public IP of a node. However, if you do this, this node will be a single point of failure. For a HA setup, you would typically choose a different options – load balancers. Load balancers are not managed directly by Kubernetes. Instead, Kubernetes will ask the underlying cloud provider to create a load balancer for you which is then connected to the service – so there might be additional charges. Creating a service exposed via a load balancer is easy – just change the type field in the manifest file to LoadBalancer apiVersion: v1 kind: Service metadata: name: alpine-service spec: selector: app: alpine ports: - protocol: TCP port: 8080 targetPort: 80 type: LoadBalancer  After applying this manifest file, it takes a few seconds for the load balancer to be created. Once this has been done, you can find the external DNS name of the load balancer (which AWS will create for you) in the column EXTERNAL-IP of the output of kubectl get svc. Let us extract this name and curl it. This time, we use the jsonpath option of the kubectl command instead of jq. $ host=$(kubectl get svc alpine-service --output \ jsonpath="{.status.loadBalancer.ingress[0].hostname}")$ curl $host:8080 <h1>It works!</h1>  If you get a “couldn not resolve hostname” error, it might be that the DNS entry has not yet propagated through the infrastructure, this might take a few minutes. What has happened? Behind the scenes, AWS has created an elastic load balancer (ELB) for you. Let us describe this load balancer. $ aws elb describe-load-balancers --output json
{
{
"Policies": {
"OtherPolicies": [],
},
"AvailabilityZones": [
"eu-central-1a",
"eu-central-1c",
"eu-central-1b"
],
"Subnets": [
"subnet-06088e09ce07546b9",
"subnet-0d88f92baecced563",
],
"CreatedTime": "2019-03-17T11:07:41.580Z",
"SecurityGroups": [
"sg-055b253a63c7aba0a"
],
"Scheme": "internet-facing",
"VPCId": "vpc-060469b2a294de8bd",
"HealthCheck": {
"UnhealthyThreshold": 6,
"Interval": 10,
"Target": "TCP:30829",
"HealthyThreshold": 2,
"Timeout": 5
},
"BackendServerDescriptions": [],
"Instances": [
{
"InstanceId": "i-0cf7439fd8eb65858"
},
{
"InstanceId": "i-0fda48856428b9a24"
}
],
"SourceSecurityGroup": {
"OwnerAlias": "979256113747"
},
"ListenerDescriptions": [
{
"PolicyNames": [],
"Listener": {
"InstancePort": 30829,
"Protocol": "TCP",
"InstanceProtocol": "TCP"
}
}
],
"CanonicalHostedZoneNameID": "Z215JYRZR1TBD5"
}
]
}


This is a long output, let us see what this tells us. First, there is a list of instances, which are the instances of the nodes in your cluster. Then, there is the block ListenerDescriptions. This block specificies, among other things, the load balancer port (8080 in our case, this is the port that the load balancer exposes) and the instance port (30829). You will note that these are also the ports that kubectl get svc will give you. So the load balancer will send incoming traffic on port 8080 to port 30829 of one of the instances. This in turn is a host port as discussed before, and therefore will be connected to our service. Thus, even though technically not fully correct, the following picture emerges (technically, a service is not a process, but a collection of iptables rules on each node, which we will look at in more detail in a later post).

Using load balancers, however, has a couple of disadvantages, the most obvious one being that each load balancer comes with a cost. If you have an application that exposes tens or even hundreds of services, you clearly do not want to fire up a load balancer for each of them. This is where an ingress comes into play, which can distribute incoming HTTP(S) traffic across various services and which we will study in one of the next posts.

There is one important point when working with load balancer services – do not forget to delete the service when you are done! Otherwise, the load balancer will continue to run and create charges, even if it is not used. So delete all services before shutting down your cluster and if in doubt, use aws elb describe-load-balancers to check for orphaned load balancers.

## Creating services in Python

Let us close this post by looking into how services can be provisioned in Python. First, we need to create a service object and populate its metadata. This is done using the following code snippet.

service = client.V1Service()
service.api_version = "v1"
service.kind = "Service"


Now we assemble the service specification and attach it to the service object.

spec = client.V1ServiceSpec()
selector = {"app": "alpine"}
spec.selector = selector
port = client.V1ServicePort(
port = 8080,
protocol = "TCP",
target_port = 80 )
spec.ports = [port]
service.spec = spec


Finally, we authenticate, create an API endpoint and submit the creation request.

config.load_kube_config()
api = client.CoreV1Api()
api.create_namespaced_service(
namespace="default",
body=service)


If you have cloned my GitHub repository, you will find a script network/create_service.py that contains the full code for this and that you can run to try this out (do not forget to delete the existing service before running this).

## Python up an EKS cluster – part II

In the last post, we have seen how Python can be used to control the generation of an EKS cluster. However, an EKS cluster without any worker nodes – and hence without the ability to start pods and services – is of very limited use. Today, we therefore take a look at the process of adding worker nodes.

Again, we roughly follow the choreography of Amazons Getting started with EKS guide, but instead of going through the individual steps using the console first, we start using the API and the Python SDK right away.

This description assumes that you have started an EKS cluster as described in my last post. If not, you can do this by running the create_cluster Python script that executes all necessary steps.

Essentially, we will do two things. First, we will create an AWS Auto-Scaling group that we will use to manage our worker nodes. This can be done using the AWS API and the Python SDK. Part of this setup will be a new role, that, in a second step, we need to make known to Kubernetes so that our nodes can join the cluster. This will be done using the Kubernetes API. Let us now look at each of these steps in turn.

## Creating an auto-scaling group

An auto-scaling group is essentially a group of EC2 instances that are combined into one management unit. When you set up an auto-scaling group, you specify a scaling policy and AWS will apply that policy to make sure that a certain number of instances is automatically running in your group. If the number of instances drops below a certain value, or if the load increases (depending on the policy), then AWS will automatically spin up new instances for you.

We will use AWS CloudFormation to set up the auto-scaling group, using the template referenced in the “Getting started” guide. This is not an introduction into CloudFormation, but it is instructive to look at this template, stored at this URL.

If you look at this file, you will recognize a few key elements. First, there are parameters whose values we will have to provide when we use the template. These parameters are

• A unique name for the node group
• The name of an SSH key pair that we will later use to log into our nodes. I recommend to create a separate key pair for this, and I will assume that you have done this and that this key pair is called eksNodeKey
• An instance type that AWS will use to spin up nodes
• The ID of an AMI image that will be used to spin up the nodes. In this tutorial, we will use Kubernetes v1.11 and the specific AWS provided AMIs which are listed in the Getting started guide. I work in the region eu-frankfurt-1, and will therefore use the AMI ID ami-032ed5525d4df2de3
• The minimum, maximum and desired size of the auto-scaling group. We will use one for the minimum size, two for the desired size and three for the maximum size, so that our cluster will come up with two nodes initially
• The size of the volume attached to each node – we will not provide a value for this and stick to the default (20 GB at the time of writing)
• In addition, we will have to supply the VPC, the subnets and the security group that we used before

In addition to the parameters, the CloudFormation template contains a list of resources that will be created. We see that in addition to several security groups and the auto scaling group, a IAM role called the node instance role and a corresponding instance profile is created. This is the role that the nodes will assume to interact with the Kubernetes master, and behind the scenes, the AWS IAM authenticator will be used to map this IAM role to a Kubernetes user, so that the nodes can authenticate themselves against the Kubernetes API. This requires some configuration, more precisely a mapping between IAM roles and Kubernetes users, which we will set up in our second steps in the next section.

So how do we apply this cloud formation template? Again, the AWS API comes to the rescue. It provides an API endpoint (client) for CloudFormation which has a method create_stack. This method expects a Python data structure that contains the parameters that we want to use, plus the URL of the cloud formation template. It will then take the template, apply the parameters, and create the resources defined in the template.

Building the parameter structure is easy, assuming that you have already stored all configuration data in Python variables (we can use the same method as in the last post to retrieve cluster data like the VPC ID). The only tricky part is the list of subnets. In CloudFormation, a parameter value that is a list needs to be passed as a comma separated list, not as a JSON list (it took me some failed attempts to figure out this one). Hence we need to convert our list of subnets into CSV format. The code snippet to set up the parameter map then looks as follows.

params = [
{'ParameterKey' : 'KeyName' ,
'ParameterValue' : sshKeyPair },
{'ParameterKey' : 'NodeImageId' ,
'ParameterValue' : amiId },
{'ParameterKey' : 'NodeInstanceType' ,
'ParameterValue' : instanceType },
{'ParameterKey' : 'NodeAutoScalingGroupMinSize' ,
'ParameterValue' : '1' },
{'ParameterKey' : 'NodeAutoScalingGroupMaxSize' ,
'ParameterValue' : '3' },
{'ParameterKey' : 'NodeAutoScalingGroupDesiredCapacity' ,
'ParameterValue' : '2' },
{'ParameterKey' : 'ClusterName' ,
'ParameterValue' : clusterName },
{'ParameterKey' : 'NodeGroupName' ,
'ParameterValue' : nodeGroupName },
{'ParameterKey' : 'ClusterControlPlaneSecurityGroup' ,
'ParameterValue' : secGroupId },
{'ParameterKey' : 'VpcId' ,
'ParameterValue' : vpcId },
{'ParameterKey' : 'Subnets' ,
'ParameterValue' : ",".join(subnets) },
]


Some of the resources created by the template are IAM resources, and therefore we need to confirm explicitly that we allow for this. To do this, we need to specify a specific capability when calling the API. With that, our call now looks as follows.

cloudFormation = boto3.client('cloudformation')
cloudFormation.create_stack(
StackName = stackName,
TemplateURL = templateURL,
Parameters = params,
Capabilities = ['CAPABILITY_IAM'])


Here the template URL is the location of the CloudFormation template provided above, and stack name is some unique name that we use for our stack (it is advisable to include the cluster name, so that you can create and run more than one cluster in parallel).

When you run this command, several things will happen, and it is instructive to take a look at the EC2 console and the CloudFormation console while the setup is in progress. First, you will find that a new stack appears on the CloudFormation console, being in status “In Creation”.

Second, the EC2 console will display – after a few seconds – a newly created auto scaling group, with the parameters specified in the template, and a few additional security groups.

And when you list your instances, you will find that AWS will spin up two instances of the specified type. These instances will come up as usual and are ordinary EC2 instances.

Remember that we specified a key pair when we created the auto-scaling group? Hopefully, you did save the PEM file – if yes, you can now SSH into your nodes. For that purpose, you will have to modify the security group of the nodes and open the SSH port for inbound traffic from your own network. Once this has been done, you can simply use an ordinary SSH command to connect to your machines.

$ssh -i eksNodeKey.pem ec2-user@18.197.140.230  where of course you need to replace 18.197.140.230 with whatever public IP address the instance has. You can now inspect the node, for instance doing docker ps to see which docker containers are already running on the node, and you can use ps ax to see which processes are running. ## Attaching the auto-scaling group to Kubernetes At this point, the nodes are up and running, but in order to register with Kubernetes, they need to talk to the Kubernetes API. More specifically, the kubelet running on the node will try to connect to the cluster API and therefore act as a Kubernetes client. As every Kubernetes client, this requires some authorization. This field is a bit tricky, but let me try to give a short overview (good references for what follows are the documentation of the AWS IAM authenticator, this excellent blog post, and the Kubernetes documentation on authorization). Kubernetes comes with several possible authentication methods. The method used by EKS is called bearer token. With this authentication method, a client request contains a token in its HTTPS header which is used by Kubernetes to authorize the request. To create and verify this token, a special piece of software is used – the AWS IAM authenticator. This utility can operate in client mode – it then generates a token – and in server mode, where it verifies a token. The basic chain of events when a client like kubectl wants to talk to the API is as follows. • The client consults the configuration file ~/.kube/config to determine the authorization method to be used • In our case, this configuration refers to the AWS IAM authenticator, so the client calls this authenticator • The authenticator looks up the current AWS credentials, uses them to generate a token and returns that token to the client • The client sends the token to the API endpoint along with the request • The Kubernetes API server extracts the token and invokes the authenticator running on the Kubernetes management nodes (as a server) • The authenticator uses the token to determine the IAM role used by the client and to verify it • It then maps the IAM role to a Kubernetes user which is then used for authorization You can try this yourself. Log into one of the worker nodes and run the command $ aws sts get-caller-identity


This will return the IAM role ARN used by the node. If you compare this with the output of the CloudFormation stack that we have used to create our auto-scaling group, you will find that the node assumes the IAM role created when our auto-scaling group was generated. Let us now verify that this identity is also used when the AWS IAM authenticator creates a token. For that purpose, let us manually create a token by running – still on the node – the command

$aws-iam-authenticator token -i myCluster  This will return a dictionary in JSON format, containing a key called token. Extract the value of this key, this will be a base64 encoded string starting with “k8s-aws-v1”. Now, fortunately the AWS authenticator also offers an operation called verify that we can use to extract the IAM role again from the token. So run $ aws-iam-authenticator verify -i myCluster -d ""k8s-aws-v1..."


where the dots represent the full token extracted above. As a result, the authenticator will print out the IAM role, which we again recognize as the IAM role created along with the auto-scaling group.

Essentially, when the kubelet connects to the server and the AWS IAM authenticator running on the server receives the token, it will use the same operation to extract the IAM role. It then needs to map the IAM role to a Kubernetes role. But where does this mapping information come from?

The answer is that we have to provide it using the standard way to store configuration data in a Kubernetes cluster, i.e. by setting up a config map.

Again, Amazon provides a template for this config map in YAML format, which looks as follows.

apiVersion: v1
kind: ConfigMap
name: aws-auth
namespace: kube-system
data:
mapRoles: |
- rolearn:
groups:
- system:bootstrappers
- system:nodes


Thus our config map contains exactly one key-value pair, with key mapRoles. The value of this key is again a string in YAML format (note that the pipe is an escape character in YAML which instructs us to treat the next lines as a multi-line string, not as a list). Of course, we have to replace the placeholder by the ARN of the role that we want to map, i.e. the role created along with the auto-scaling group.

To create this config map, we will have to use the Python SDK for Kubernetes that will allow us to talk to the Kubernetes API. To install it, run

$pip3 install kubernetes  To be able to use the API from within our Python script, we of course have to import some components from this library and read our kubectl configuration file, which is done using the following code. from kubernetes import client, config config.load_kube_config() v1 = client.CoreV1Api()  The object that we call v1 is the client that we will use to create and submit API requests. To create our config map, we first need to instantiante an object of the type V1ConfigMap and populate it. Assuming that the role ARN we want to map is stored in the Python variable nodeInstanceRole, this can be done as follows. body = client.V1ConfigMap() body.api_version="v1" body.metadata = {} body.metadata['name'] = "aws-auth" body.metadata['namespace'] = "kube-system" body.data = {} body.data['mapRoles'] = "- rolearn: " + nodeInstanceRole + "\n username: system:node:{{EC2PrivateDNSName}}\n groups:\n - system:bootstrappers\n - system:nodes\n"  Finally, we call the respective method to actually create the config map. v1.create_namespaced_config_map("kube-system", body)  Once this method completes, our config map should be in place. You can verify this by running kubectl on your local machine to describe the config map. $ kubectl describe configMaps aws-auth -n kube-system
Name:         aws-auth
Namespace:    kube-system
Labels:
Annotations:

Data
====
mapRoles:
----
- rolearn: arn:aws:iam::979256113747:role/eks-auto-scaling-group-myCluster-NodeInstanceRole-HL95F9DBQAMD
groups:
- system:bootstrappers
- system:nodes

Events:


If you see this output, we are done. The EC2 instances should now be able to connect to the Kubernetes master nodes, and when running

$kubectl get nodes  you should get a list displaying two nodes. Congratulations, you have just set up your first Kubernetes cluster! The full Python script which automates these steps can as usual be retrieved from GitHub. Once you are done playing with your cluster, do not forget to delete everything again, i.e. delete the auto-scaling group (which will also bring down the nodes) and the EKS cluster to avoid high charges. The best approach to clean up is to delete the CloudFormation stack that we have used to bring up the auto-scaling group (in the AWS CloudFormation console), which will also terminate all nodes, and then use the AWS EKS console to delete the EKS cluster as well. This completes our post for today. In the next post in this series, we will learn how we can actually deploy containers into our cluster. ## Python up an EKS cluster – part I When you want to try out Kubernetes, you have several choices. You can install Kubernetes in a cluster, install it locally using Minikube, or use one of the Kubernets offerings of the major cloud providers like AWS, GCP or Azure. In this post, we set up a Kubernetes cluster on Amazons EKS platform. As I love Python, we will do this in two steps. First, we will go through the steps manually, as described in the AWS documentation, using the AWS console. Then, we will learn how to use the AWS Python SDK to automate the recurring steps using Python. A word of caution before we continue: when you follow these instructions, Amazon will charge you for the cluster as well as for the nodes that the cluster uses! So make sure that you stop all nodes and delete the cluster if you are done to avoid high charges! At the time of writing, Amazon charges 20 cents per hour per cluster plus the EC2 instances underlying the cluster, so if you play with this for one our two hours it will not cost you a fortune, but if you forget to shut everything down and let it run, it can easily add up! ## What is EKS? Before we get our hands dirty, let us quickly discuss what EKS actually is. Essentially, EKS offers you the option to set up a Kubernetes cluster in the Amazon cloud. You can connect to this cluster using the standard Kubernetes API and the standard Kubernetes tools. The cluster will manage worker nodes that are running on Amazons EC2 platform, using your EC2 account. The management nodes on which Kubernetes itself is running are not running on your EC2 instances, but are provided by Amazon. When you set up an EKS cluster, the cluster and your worker nodes will be running in a virtual private network (VPC). You can add load balancers to make your services reachable from the internet. In addition, EKS is integrated with Amazons IAM. This tutorial assumes that you own an AWS account and that you have followed the instructions in the IAM getting started guide to set up an IAM user with administrator privileges. Please make sure that you are logged into the AWS console with this IAM user, not your root account. To set up an EKS cluster, several preparational steps are required. First, you will need to set up an IAM role that EKS will use to control EC2 nodes on your behalf. Next, you will need a VPC in which your cluster will run. After these one time preparation steps have been completed, you can create a Kubernetes cluster, add worker nodes and start deployments. In the following, we will go through most of these steps one by one. ## One time preparational steps Before you can use EKS, you will need to add an IAM role. To do this, navigate to the EKS console (note that this link will automatically take you to your home region). As described on the Getting started with EKS page, open the IAM management console. In the menu on the left, chose “Roles” and then hit the button “Create role”. Choose “EKS” from the list of services and hit “Next: permissions” and then immediately “Next: tags” and “Next: Review”. As the role name, enter “EKSEC2UserRole” (you could of course choose any name you want, but to make the scripts that we will establish further below work, choose this name for this tutorial). At this point, your screen should look as follows. Now confirm the setup of the new role with “Create role” and observe how the new role appears in the list in your IAM console. The next ingredient that we need is the VPC in which the cluster will be running. To set this up, navigate to the CloudFormation page. Hit the button “Create new stack”. Select “Specify an Amazon S3 template URL” and enter https://amazon-eks.s3-us-west-2.amazonaws.com/cloudformation/2019-02-11/amazon-eks-vpc-sample.yaml . Then hit “Next”. On the next page, you will have to choose a stack name. Again, this can be any name, but let us follow the standard suggestion and call it “eks-vpc”. Then hit “Next” twice and then “Create”. AWS will now create the VPC for you, which might take a few minutes, you should see a screen as below while this is in progress. Once the stack has been created, open the “Outputs” tab and record the value of SecurityGroups, VpcID and SubnetIds. . In my case, the values are Key Value SecurityGroups sg-005cf103994878f27 VpcId vpc-060469b2a294de8bd SubnetIds subnet-06088e09ce07546b9, subnet-0e4a8fd662faadab6, subnet-0d88f92baecced563 Next, there is some software that you need on your PC. Please follows the instructions on the Getting started with EKS page to download and install kubectl, the aws-iam-authenticator and the aws CLI your your platform. On my Ubuntu Linux PC, this was done using the steps below. $ sudo snap install kubectl --classic
$curl -o aws-iam-authenticator https://amazon-eks.s3-us-west-2.amazonaws.com/1.11.5/2018-12-06/bin/linux/amd64/aws-iam-authenticator$ curl -o aws-iam-authenticator.sha256 https://amazon-eks.s3-us-west-2.amazonaws.com/1.11.5/2018-12-06/bin/linux/amd64/aws-iam-authenticator.sha256
$openssl sha1 -sha256 aws-iam-authenticator$ cat aws-iam-authenticator.sha256 # compare checksums!
$chmod +x ./aws-iam-authenticator$ mv aws-iam-authenticator $HOME/Local/bin # replace by your local bin directory$ pip3 install awscli --upgrade --user


Note that pip will install the AWS CLI in $HOME/.local/bin, so make sure to add this directory to your path. We can now test that all tools have been installed by entering $ kubectl help
$aws --version  At this point, the aws tool is not yet hooked up with your credentials. To do this, make sure that you have access to the IAM access key ID and access secret key (that you can generate at the IAM console) for your IAM admin user. Then run $ aws configure


Enter the keys and the required configuration information. When done, enter

$aws sts get-caller-identity  to confirm that you now access the AWS API with the admin user. ## Creating a cluster and adding worker nodes using the AWS console After all these preparations – which we of course have to go through only once – we are now in a position to create our first EKS cluster. For that purpose, go to the EKS console and look for the box called “Create EKS cluster”. Enter a cluster name – I have chosen “myCluster” (yes, I am very creative) – and hit the button “Next step”. This will take you to a configuration page. On this page, you will have to select an IAM role that the cluster will use to operate EC2, a VPC and a security group. Make sure to select those entities that we have set up above! Then hit “Create”. Your cluster will now be created, and you should be taken to the EKS cluster overview page that lists all your clusters along with their status. Creation of a cluster can take a few minutes, be patient. At some point, the status of your cluster should switch to active. Congratulations, you are now proud owner of a running Kubernetes cluster! However, this Kubernetes cluster is pretty useless – there are no worker nodes yet, and you are not yet able to communicate with your cluster using kubectl. We will create worker nodes and run services in the next post. To fix the second problem, we will have to hook up the IAM credential chain used by the AWS CLI tool with kubectl. To do this, enter $ aws eks --region eu-central-1 update-kubeconfig --name myCluster


Once this has completed, you should now be able to user kubectl to connect to your cluster. As an example, run

$kubectl get svc NAME TYPE CLUSTER-IP EXTERNAL-IP PORT(S) AGE kubernetes ClusterIP 10.100.0.1 443/TCP 6m  ## Creating a cluster automatically using the Python SDK Creating a cluster using the EKS console is a bit time consuming, and given the fact that Amazon will charge you for a running cluster, you will want to delete your cluster when you are done and recreate it easily when you need it again. For that purpose, it is very useful to be able to create a cluster using Python. To do this, let us first delete the cluster again (on the cluster overview page in the EKS console) that we have just created, and compose a Python script that recreates the cluster for us. First, of course, we need to install the Python SDK Boto3. This is as easy as $ pip3 install boto3


Note that the SDK uses the same credentials as the AWS CLI, so you still need to go through the steps above to install and configure the AWS CLI.

Our script will accept the name of a cluster (which was myCluster in the example above). It will then need to determine the IAM role, the VPC and subnets and the security group as we did it for the manual setup.

To determine the IAM role, we need an IAM client and use its method get_role to get role details and obtain the role ARN.

iam = boto3.client('iam')
r = iam.get_role(RoleName=eksRoleName)
roleArn = r['Role']['Arn']


Here we assume that the Python variable eksRoleName contains the name of the role that we have created as part of our one time setup.

Getting the VPC ID, the subnet IDs and the security group is a bit more tricky. For that purpose, we will refer back to the CloudFormation stack that we have created as part of our one-time setup. To get the VPC ID and the security group, we can use the following code snippet.

cloudFormation = boto3.client('cloudformation')
r = cloudFormation.describe_stack_resources(
StackName = "eks-vpc",
LogicalResourceId="VPC")
vpcId = r['StackResources'][0]['PhysicalResourceId']
r = cloudFormation.describe_stack_resources(
StackName = "eks-vpc",
LogicalResourceId="ControlPlaneSecurityGroup")
secGroupId = r['StackResources'][0]['PhysicalResourceId']


Once we have this, we can easily obtain the subnet IDs using the EC2 interface.

ec2 = boto3.resource('ec2')
vpc = ec2.Vpc(vpcId)
subnets = [subnet.id for subnet in vpc.subnets.all()]


We now have all the configuration data that we need to trigger the actual cluster setup. For that purpose, we will create an EKS client and call its method create_cluster, passing the information that we have collected so far as arguments. We then create a waiter object that is polling until the cluster status changes to “Active”.

eks = boto3.client("eks")
response = eks.create_cluster(
name="myCluster",
version="1.13",
roleArn = roleArn,
resourcesVpcConfig = {
'subnetIds' : subnets,
'securityGroupIds' : [secGroupId]})
waiter = eks.get_waiter("cluster_active")
waiter.wait(name="myCluster")


When this call completes, our cluster is created. However, we still have to update our Kubernetes configuration file. Of course, we could simply call the AWS CLI, but doing this in Python is more fun. The configuration file that the AWS CLI creates is stored in your home directory, in a subdirectory called kube, and formatted in the simplified markup language known as YAML. Inspecting this file, it is not difficult to translate this into a structure of Python dictionaries and lists, which contains four cluster-specific pieces of data.

• The cluster endpoint – this is the URL used by kubectl to submit an API call
• The Amazon resource name of your cluster (ARN)
• the cluster name
• a certificate used by kubectl to authorize requests

This information can be retrieved using the method describe_cluster of the EKS client and assembled into a Python data structure. We can then use the Python module yaml to turn this into a string in YAML format, which we can then write to disk. I will not go into detail, but you might want to take a look at the complete script to create a cluster on my GitHub page.

We are now able to automatically create a Kubernetes cluster and connect to it using kubectl. In the next post, we will add the meat – we will learn how to spin up worker nodes and deploy our first pods.

## Factoring integers on a quantum computer with Qiskit

After all the work done in the previous posts, we are now ready to actually implement Shor’s factoring algorithm on a real quantum computer, using once more IBMs Q Experience and the Qiskit framework.

First, recall that Shor’s algorithm is designed to factor an integer M, with the restriction that M is supposed to be odd and not a prime power. Thus the smallest meaningful value of M is M=15. Of course this is a toy example, but we will see that even this toy example can be challenging on the available hardware. Actually, 15 is also the number that was factored in the first implementation of Shor’s algorithm on a real quantum computer which used an NMR based device and was done in 2001 – see [1].

The quantum part of the algorithm is designed to find the period r of a chosen number a modulo M. There are several choices of a that are possible, the only condition we need is that a and M are co-prime. In this post, we follow the choice in [1] and use a = 11 and a = 7. As in [1] (and all other similar papers I am aware of, see also the discussion in this paper by Smolin, Smith which later appeared in Nature) we will also rely on prior knowledge of the result to build up our circuits, so we are actually cheating – this is not really an application of Shor’s algorithm, but merely a confirmation of a known factoring.

Of course, it is easy to compute the period directly in our case. The period of a = 11 is two, and the period of a = 7 is four. A smaller period means a simpler circuit, and therefore we start with the case a = 11. Let us see whether we can confirm the period.

## The easy case – a = 11

For our tests, we will use the description of Shor’s algorithm in terms of the quantum phase estimation procedure. Thus we have a primary register that we initialize in the state $|1\rangle$ and to which we apply a sequence of controlled multiplications by powers of a, and a working register, in which we do a quantum Fourier transform in the second step. The primary register needs to hold M=15, so we need four qubits there. For the working register, we use three qubits, which is the smallest possible choice, to keep the number of gates small.

The first step of our algorithm is controlled multiplication by a. However, this is greatly simplified by the fact that we apply this to exactly one state that we know upfront – the initial state $1 \rangle$. Thus we only need to implement a unitary transformation which conditionally maps the state $|1 \rangle$ to the state $|a = 11 \rangle$. In binary notation, eleven is 1011, and we therefore simply need to conditionally flip qubits 3 and 1. We also need a Pauli X gate in the primary register to build the initial state $|1 \rangle$ from the fiducial state $|0\rangle$ and Hadamard gates in the working register to create the initial superposition. Thus the first part of our circuit looks as follows.

What about higher powers of a? This is easy – we know that a2 is equal to one modulo M, and therefore the multiplication by a2 modulo M is simply the identity transformation. The same is true for all even powers and so we are already done with the first part of our circuit – this is why we call the case a = 11 the easy case. The only thing that is left is to add the quantum Fourier transformation and a measurement to the working register. This gives us the following circuit (which can obviously be simplified, more on this later).

Now we can run this on a simulator. Before we do this, let us try to understand what we expect. We know that after measuring the working register, the value that we obtain is a multiple of 23 / r = 8 / 2 = 4. Thus we expect to see peaks at 0 and four (binary 100). And this is actually what we get – the following histogram shows the output of an execution on the local QASM simulator integrated into Qiskit.

This is nice, our circuit seems to work correctly. So let us proceed and try this on real hardware. Before we do so, however, let us simplify our circuit a bit to reduce the number of gates needed and thus the noise level. First, we skip the final swap gates in the QFT circuit – this will change our expected output from 100 to 001, but we can keep track of this manually when setting up the measurement. Next, we have two Hadamard gates on w[2] that cancel and that we can therefore remove. And the Pauli X gate on p[0] is never really used and can be dropped. After these changes, we obtain the following circuit.

We can run this on the IBMQ hardware. As we require seven qubits in total, we need to do this on the 14 qubit model ibmq_16_melbourne. Here are the results of a test run, again as a histogram.

We see the expected peaks at 000 and 100, but also see some significant noise that is of course not present in the simulation. However, the probabilities to measure one of the theoretically expected values is still roughly twice as high as the probability of any of the other values.

## The hard case – a = 7

Having mastered the case a = 11, let us now turn to the case a = 7. The only change that we need to make is the circuit acting on the primary register. In the first part of the circuit, we can again make use of the fact that we only have to deal with one possible input, namely the initial state $|1 \rangle$. The conditional multiplication by a maps this to $|7\rangle$, and we can again implement this by two conditional bit flips, i.e. two CNOT gates.

The second part, the conditional multiplication by a2 = 4 modulo M, requires a bit more gates. On a bit level, multiplication by four is a bit shift by two bits to the left. We can realize this as a sequence of two conditional swap gates, swapping the bits zero and two and the bits one and three. A conditional swap gate can be implemented by three CNOT gates, or, in QASM syntax,

gate cswap a,b,c
{
cx c,b;
ccx a,b,c;
cx c,b;
}


We can also make a few simplifications by replacing CNOT gates with fixed values of the control qubit by either an inversion or the identity. We arrive at the following circuit.

This is already considerably more complex than for the case a = 11. When we add the QFT circuit and the measurements, we end up with the following circuit.

Now the period is four. Thus we expect peaks at all multiples of two, i.e. at all even values. And this is actually what we get when we run this on a simulator.

Now let us again see how we can optimize this circuit before trying it on real hardware. Again, there are a few additional simplifications that we can make to reduce the noise level. On w[2], we have a double Hadamard gate that we can remove. We can again dispose of the final swap gates and move the swap operation into the measurement. The last CNOT between p[1] and p[3] can be dropped as it does not affect the outcome. The same is true for the last CNOT between p[0] and p[2]. Thus we finally arrive at the following circuit.

And here is the result of running this on the 16 qubit IBM device.

Surprisingly, the noise level is not significantly worse than for the version with a = 11, even though we have a few more gates. We can clearly see that the probabilities to measure even values are significantly higher than the probabilities to measure odd values, corresponding to the period four.

So what did we learn from all this? First, our example was of course a toy model – we did code the circuit based on knowledge on the period and hard-coded the value of a. However, our example demonstrates the basic techniques to build circuits for the more general case. In real application, we could, given a number M, first determine a suitable value for a. The conditional multiplication by a is then acting as a permutation on the computational basis and can therefore be implemented by a sequence of (conditional) swap gates. The same is true for higher powers of a. So in theory, we have all the ingredients that we need to implement more general versions of the algorithm.

In practice, however, we will quickly reach a point where the number of gates exceeds the limit that we can reasonable implement given the current level of noise. So real applications of the algorithm in number ranges that are not tractable for classical digital computers will require advanced error correction mechanisms and significantly reduced noise levels. Still, it is nice (and fun) to see a toy version of the famous algorithm in action on a real device.

If you want to run the examples yourself, you can use my notebooks for the case a=11 and a=7.

## References

1. L. Vandersypen, M. Steffen, G. Breyta, C. Yannoni, M. Sherwood, I. Chuang, Experimental realization of Shor’s quantum factoring algorithm using nuclear magnetic resonance, Nature Vol. 414, (December 2001)
2. P. Shor, Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer, SIAM J.Sci.Statist.Comput. Vol. 26 Issue 5 (1997), pp 1484–1509, available as arXiv:quant-ph/9508027v2
3. R. Cleve, A. Ekert, C. Macchiavello, M. Mosca, Quantum Algorithms revisited, arXiv:9708016
4. A. Kitaev, Quantum measurements and the Abelian Stabilizer Problem, arXiv:quant-ph/9511026
5. J. Smolin, G. Smith, A. Vargo, Pretending to factor large numbers on a quantum computer, arXiv:1301.7007 [quant-ph]

## Implementing the quantum Fourier transform with Qiskit

The quantum Fourier transform is a key building block of many quantum algorithms, from Shor’s factoring algorithm over matrix inversion to quantum phase estimation and simulations. Time to see how this can be implemented with Qiskit.

Recall that the quantum Fourier transform (or, depending on conventions, its inverse) is given by

$|x \rangle \mapsto \frac{1}{\sqrt{2^n}} \sum_s \eta^{xs} |s \rangle$

where $\eta$ is an 2n-th root of unity with n being the number of qubits in the register, i.e.

$\eta = \exp \frac{2\pi i}{2^n}$

How can we effectively create this state with a quantum circuit? They key to this is the observation (see the references below) that the result of the quantum Fourier transform can be written as a product state, namely as

$\sum_s \eta^{xs} |s \rangle = (|0\rangle + \eta^{2^{n-1} x}|1\rangle)(|0 \rangle + \eta^{2^{n-2}x}|1 \rangle) \cdots (|0 \rangle + \eta^x|1 \rangle)$

which you can easily verify by multiplying out the product and collecting terms.

Here we use the tensor product order that is prescribed by OpenQASM, i.e. the most significant bit is q[n-1]. This bit is therefore given by

$|0 \rangle + \eta^{2^{n-1}x} |1 \rangle$

Let us analyse this expression further. For that purpose, we decompose x into its representation as a binary number, i.e. we write

$x = \sum_{i=0}^{n-1} 2^i x_i$

with xi being the binary digits of x. If we now multiply this by 2n-1, we will get $2^{n-1} x_0$ plus a multiple of 2n. As $\eta$ is a root of unity, this multiple cancels out and we obtain that

$\eta^{2^{n-1}x} = \eta^{2^{n-1}x_0} = e^{\pi i x_0} = (-1)^{x_0}$

Thus, we can write the most significant qubit of the Fourier transform as

$|0 \rangle + (-1)^{x_0} |1 \rangle$

which is nothing but

$H |x_0\rangle$

Thus we obtain the most significant qubit of the quantum Fourier transform by simply applying a Hadamard gate to the least significant qubit of the input.

This is nice and simple, but what about the next qubit, i.e. qubit n-2? From the decomposition above, we can see that this is simply

$|0 \rangle + \eta^{2^{n-2}x} |1 \rangle$

Using exactly the same arguments as for the most significant qubit, we easily find that this is

$|0 \rangle + \eta^{2^{n-2}x_0} H |x_1 \rangle$

Thus we obtain this qubit from qubit 1 by first applying a Hadamard gate and then a conditional phase gate, i.e. a conditional rotation around the z-axis, conditioned on the value of x0. In general, qubit n-j is

$|0 \rangle + \prod_{i < j - 1} \eta^{2^{n-j+i}x_i} H |x_{j-1}\rangle$

which is a Hadamard gate followed by a sequence of conditional rotations around the z-axis, conditioned on the qubits with lower significance.

So we find that each qubit of the Fourier transform is obtained by applying a Hadamard followed by a sequence of conditional rotations. However, the order of the qubits in the output is reversed, i.e. qubit n-j is obtained by letting gates act on qubit j. Therefore, at the end of the circuit, we need to revert the order of the qubits.

In OpenQASM and Qiskit, a conditional rotation around the z-axis is called CU1, and there are swap gates that we can use to implement the final reversing of the qubits. Thus, we can use the following code to build a quantum Fourier transformation circuit acting on n qubits.

def nBitQFT(q,c,n):
circuit = QuantumCircuit(q,c)
#
# We start with the most significant bit
#
for k in range(n):
j = n - k
# Add the Hadamard to qubit j-1
circuit.h(q[j-1])
#
# there is one conditional rotation for
# each qubit with lower significance
for i in reversed(range(j-1)):
circuit.cu1(2*np.pi/2**(j-i),q[i], q[j-1])
#
# Finally we need to swap qubits
#
for i in range(n//2):
circuit.swap(q[i], q[n-i-1])
return circuit


Here is the circuit that this code produces for n=4. We can clearly see the structure – on each qubit, we first act with a Hadamard gate, followed by a sequence of conditional rotations with decreasing angle, conditioned on the less significant qubits, and finally reorder the qubits.

This is already a fairly complex circuit, and we need to find a way to test it. Let us look at the options we have. First, a quantum circuit is a unitary transformation and can be described by a matrix. In our case, it is especially easy to figure out what this matrix should be. Looking at the formula for the quantum Fourier transform, we find that the matrix describing this transformation with respect to the computational basis has the elements

$U_{ij} = \frac{1}{\sqrt{2^n}} \eta^{ij}$

The Qiskit frameworks comes with a simulator called the unitary_simulator that accepts a quantum circuit as input and returns the matrix describing that circuit. Thus, one possible test approach could be to build the circuit, run the unitary simulator on it, and to compare the resulting unitary matrix with the expected result given by the formula above. In Python, the expected result is produced by the following code

def qftMatrix(n):
qft = np.zeros([2**n,2**n], dtype=complex)
for i in range(2**n):
for j in range(2**n):
qft[i,j] = np.exp(i*j*2*1j*np.pi/(2**n))
return 1/np.sqrt(2**n)*qft


and the test can be done using the following function

def testCircuit(n):
q = QuantumRegister(n,"x")
c = ClassicalRegister(n,"c")
circuit = nBitQFT(q,c,n)

backend = Aer.get_backend('unitary_simulator')
job = execute(circuit, backend)
actual = job.result().get_unitary()
expected = qftMatrix(n)
delta = actual - expected
print("Deviation: ", round(np.linalg.norm(delta),10))
return circuit


The outcome is reassuring – we find that the matrices are the same within the usual floating point rounding differences.

After passing this test, a next reasonable validation step could be to run the algorithm on a specific input. We know that the QFT will map the state $|0 \rangle$ into an equal superposition of all elements of the computational basis. Conversely, we therefore expect that if we start with such a superposition, the QFT will map the superposition onto $|0\rangle$, at least up to a phase.

Let us try this out. Our test circuit will consist of a layer of Hadamard gates to create the equal superposition, followed by the QFT circuit, followed by a measurement. The resulting circuit for n=4 is displayed below.

It we run this circuit on the QASM simulator embedded into Qiskit, the result is as expected – for 1024 shots, we get 1024 times the output ‘0000’. So our circuit works – at least theoretically. But what about real hardware?

Let us compile and run the circuit targeting the IBM Q Experience 14 qubit device. If we dump the QASM code after compilation, we see that the overall circuit will have roughly 140 gates. This is already significant, and we expect to see some noise. To see how bad it is, I have conducted several test runs and plotted the results as a histogramm (if you want to play with this yourself, you will find the notebook on Github). Here is the output of a run with n=4.

We still see a clear peak at the expected result, but also see that the noise level is close to making the result unusable – if we did not know the result upfront, we would probably not dare to postulate anything from this output. With only three qubits, the situation becomes slightly better but is still far from being satisfactory.

Of course we could now start to optimize the circuit – remove cancelling Hadamard gates, remove the final swap gates, reorder qubits to take the coupling map into account and so on – but it becomes clear that with the current noise level, we are quickly reaching a point where even comparatively simple circuit will inflict a level of noise that is at best difficult to handle. Hopefully, this is what you expected after reading my posts on quantum error correction, but I found it instructive to see noise in action in this example.

## References

1. M. Nielsen, I. Chuang, Quantum Computation and Quantum Information, Cambridge University Press 2010
2. R. Cleve, A. Ekert, C. Macchiavello, M. Mosca, Quantum Algorithms revisited, arXiv:9708016

## Running the Deutsch-Jozsa algorithm on IBMs Q experience

In one of the previous posts, we have looked at the basics of the Qiskit package that allows us to create and run quantum algorithms in Python. In this post, we will apply this to model and execute a real quantum algorithm – the Deutsch-Jozsa algorithm.

Recall that the Deutsch-Jozsa algorithm is designed to solve the following problem. We are given a boolean function f

$f \colon \{0,1\}^n \rightarrow \{0,1\}$

in n variables. We know that the function is either constant or it is balanced (i.e. the number of times it is zero is equal to the number of times it is equal to 1). The task is to determine which of the to properties – balanced or constant – the function f has.

The first choice that we need to make is of course a suitable function f. To keep things simple, we will use a function with two variables. Maybe the most straightforward example for a balanced function in two variables is the classical XOR function.

$f(x_0, x_1) = x_0 \oplus x_1$

Thus our first task is to develop a quantum equivalent of this function, i.e. a reversible version of f.

Recall that in general, given a boolean function f, we can define a unitary transformation Uf by the following rule

$U_f(|x \rangle |y \rangle) = |x \rangle |y \oplus f(x) \rangle$

Thus Uf flips the target qubit y if f(x) is one and leaves it unchanged otherwise. In our case, combining this with the definition of the XOR function implies that we are looking for a circuit acting on three qubits – two input qubits and one target qubit – that flips the target qubit if and only if the two input qubits are not equal. From this textual description, it is obvious that this can be constructed using a sequence of two CNOT gates.

Let us now go through the individual steps of the Deutsch-Jozsa algorithm, using the optimized version published 1997 in this paper by R. Cleve, A. Ekert, C. Macchiavello and M. Mosca (this version is also described in one of my earlier posts but is does not hurt to recap some of that). The first step of the algorithm is to prepare a superposition state

$|\psi_0 \rangle = \frac{1}{2} \sum_x |x \rangle$

This is accomplished by applying a Hadamard gate to each of the two qubits that make up our primary register in which the variable x lives. We then add an ancilla qubit that is in the state $H|1 \rangle$, so that our state is now

$\frac{1}{2} \sum_x |x \rangle \otimes H |1 \rangle = \frac{1}{2\sqrt{2}} \sum_x |x , 0 \rangle - |x, 1 \rangle$

So far our circuit looks as follows.

Next we apply the oracle Uf to our state. According to the definition of the oracle, we have

$U_f |x, 0 \rangle = |x\rangle |f(x) \rangle$

and

$U_f |x, 1 \rangle = |x\rangle |f(x) \oplus 1\rangle$

Therefore we find that

$U_f \colon: |x, 0 \rangle - |x, 1 \rangle \mapsto |x\rangle |f(x)\rangle - |x\rangle |f(x) \oplus 1 \rangle$

From this formula, we can read off how Uf acts depending on the value of f(x). If f(x) is 0, the right hand side is equal to the left hand side and Uf acts trivially. If, however, f(x) is one, the right hand side is simply minus the left hand side, and Uf acts as multiplication by -1. Combining this with our previous formula, we obtain

$U_f \colon: \frac{1}{2} \sum_x |x \rangle \otimes H |1 \rangle \mapsto \frac{1}{2} \sum_x (-1)^{f(x)}|x \rangle \otimes H |1 \rangle$

Next, we apply again the Hadamard operator followed by the Pauli X gate to the third qubit (the ancilla), i.e. we uncompute the ancilla qubit. From the expression above, we can read off directly that the result left in the first two qubits will be

$|\psi \rangle = \frac{1}{2} \sum_x (-1)^{f(x)}|x \rangle$

This is the vector that we will have in the first two qubits of our quantum register when we have processed the following circuit.

Let us now calculate the overlap, i.e. the scalar product, between this vector and the initial superposition $|\psi_0 \rangle$. Clearly,

$\langle \psi_0 |\psi \rangle = \frac{1}{4} \sum_x (-1)^{f(x)}$

We find that this overlap is zero if and only if the function f is balanced. But how can we measure this scalar product? To see this, recall that the initial state $|\psi_0 \rangle$ is the result of applying the tensor product H2 of the two Hadamard operators to the first two qubits in the fiducial state. Thus we can write our scalar product as

$\langle \psi_0 |\psi \rangle = \langle 0 | H^2 |\psi \rangle = \langle 0 | H^2 \psi \rangle$

In other words, we can determine the scalar product by again applying a Hadamard operator to each of the first two qubits and measuring the overlap of the resulting state with the basis state $|00 \rangle$ which is the same thing as the probability to measure $|00 \rangle$ when we perform a measurement in the computational basis. Thus we finally obtain the following circuit

and the function f is balanced if and only if the probability to measure $|00 \rangle$ is zero.

Let us now turn this into Python code using Qiskit. I found it useful to create subroutines that act on a given circuit and build parts of the circuit which can then be tested independently from each other on a simulator before combining them. Here is the code to create the oracle for our balanced function f.

def createOracleBalanced(q,c):
circuit = QuantumCircuit(q,c)
circuit.cx(q[0], q[2])
circuit.cx(q[1], q[2])
circuit.barrier(q)
return circuit


Similarly, we can write routines that create the initial state and join the ancilla qubit.

def createInitialState(circuit):
circuit.h(q[0])
circuit.h(q[1])
circuit.barrier(q)

circuit.x(q[2])
circuit.h(q[2])
circuit.barrier(q)


Finally, we need to be able to uncompute the ancilla and to add the final measurements.

def uncomputeAncilla(circuit):
circuit.h(q[2])
circuit.x(q[2])
circuit.barrier(q)

circuit.h(q[0])
circuit.h(q[1])
circuit.barrier(q)
circuit.measure(q[0], c[0])
circuit.measure(q[1], c[1])
circuit.barrier(q)


A word on barriers. In this example code, we have added barriers after each part of the circuit. Barriers are an element of the OpenQASM specification and instruct the compiler not to combine gates placed on different sides of a barrier during optimization. In Qiskit, barriers have the additional benefit of structuring the visualization of a circuit, and this is the main reason I have included them here. In an optimized version, it would probably be safe to remove them.

After all these preparations, it is now easy to compile the full circuit. This is done by the following code snippet – note that we can apply the + operator to two circuits to tell Qiskit to concatenate the circuits.

q = QuantumRegister(3,"q")
c = ClassicalRegister(3,"c")
circuit = QuantumCircuit(q,c)
createInitialState(circuit)
circuit = circuit + (createOracleBalanced(q,c))
uncomputeAncilla(circuit)


We can then compile this code for a simulator or a backend as explained in my last post on the IBM Q experience (you can also find the full source code here). The following histogramm shows the result of the final measurement after running this on the ibmqx4 5 qubit quantum computer.

We see, as in the experiments conducted before, that the theoretically expected result (confirmed by running the circuit on a simulator) is blurred by noise – we get the value 00 in a few instances, even though the function is balanced. Even though the outcome can still be derived with a reasonable likelihood, we start to get a feeling for the issues that we might have with noise for more complex functions and circuits.

It is instructive to extract the compiled QASM from the resulting qobj and load that code into the IBM Q Experience composer. After a few beautifications, the resulting code (which you can get here) is displayed as follows in the composer.

If we compare this to the original circuit, we see that the compiler has in fact rearranged the CNOT gates that make up our oracle. The reason is that the IBMQX4 device can only realize specific CNOT gates. It can, for instance, implement a CNOT gate with q[1] as control as q[0] as target, but not vice versa. Therefore the compiler has added Hadamard gates to swap control and target qubit. We also see that the compiler has respected the barriers and not cancelled some of the double Hadamard gates. If we remove the barriers and remove all Hadamard gates that clearly cancel each other, we finally obtain a greatly simplified version of the circuit which looks as follows.

Of course this simplification hinges on the special choice of the function f. Nevertheless, it is useful – fewer gates mean fewer errors. If we compare the error rates of the optimized circuit with the new circuit, we find that while the original version had an error (i.e. an unexpected amplitude of $|00 \rangle$) of roughly 10%, a run with the optimized circuit showed an error of only 5,5%.

This completes our post. The next algorithm that we will implement is Shor’s algorithm and the quantum Fourier transform.

## Using Python to access IBMs quantum computers

In a previous post, we have looked at IBMs Q experience and the graphical composer that you can use to build simple circuits and run them on the IBM hardware. Alternatively, the quantum hardware can be addressed using an API and a Python library called Qiskit which we investigate in this post.

## Installation and setup

To be able to use Qiskit, there are some setup steps that we need to complete. First, we obviously have to install Qiskit. This is easy – we can simply use pip.

pip install qiskit


This will download the latest version of Qiskit, in my case (at the time of writing this) this was version 0.6.1. The next thing that you need is an API token to be able to access the IBM API. Assuming that you are already a registered user, you can create your token on the advanced tab of your user profile page.

In order to easily access your token from a Python script, it is useful to store the token locally on your hard drive. Qiskit uses a file qiskitrc in the ~/.qiskit folder in your home directory to store your credentials. The easiest way to create this file is the following code snippet

python -c 'from qiskit import IBMQ ; IBMQ.save_account("your_token")'


where obviously you need to replace the text inside the quotes with your token. After running this, you should find a file ~/.qiskit/qiskitrc containing your saved token (in clear text).

Once that file is in place, you can now easily load the credentials from within a Python program using the following code snippet

from qiskit import IBMQ


The method IBMQ.active_accounts() will also return a list of currently available accounts which can be useful for debugging purposes. Loading an account is only needed if we want to connect to the IBM site to use their quantum hardware or the online simulator, not if we use a local backend – more on this later.

## Circuits, gates and measurements

Let us now take a look at the basic data structures of Qiskit. A QuantumCircuit is what you expect – a collection of gates and registers. In Qiskit, a circuit operates on a QuantumRegister and optionally contains a ClassicalRegister which holds the results of a measurement. The following code snippet will create a quantum register with two qubits, a classical register with two qubits and a quantum circuit based on those registers.

from qiskit import QuantumCircuit
from qiskit import ClassicalRegister
from qiskit import QuantumRegister
q = QuantumRegister(2,"q")
c = ClassicalRegister(2,"c")
circuit = QuantumCircuit(q,c)


Next, we need to add gates to our circuit. Adding gates is done by calling the appropriate methods of the circuit object. The gate model of Qiskit is based on the OpenQASM language described in [1].

First, there are one qubit gates. The basic one qubit gates are rotations, defined as usual, for instance

$R_X(\Theta) = \exp \left( -i \frac{\Theta}{2} \sigma_X \right) = \cos \frac{\Theta}{2} - i \sigma_X \sin \frac{\Theta}{2}$

and similarly for the other Pauli matrices. Now it is well known that any rotation of the Bloch sphere can be written as a product of three rotations around y- and z-axis, i.e. in the form

$R_Z(\Phi)R_Y(\Theta)R_Z(\lambda)$

which is denoted by

$U(\Theta,\Phi,\lambda)$

in OpenQASM and Qiskit. For instance, $U(0, 0, \lambda)$ is a rotation around the z-axis and so forth. A short calculation shows that

$U(\Theta,\Phi,\lambda) = \begin{pmatrix} \exp \left(-\frac{i}{2}(\Phi + \lambda)\right) \cos \frac{\Theta}{2} & - \exp \left(-\frac{i}{2}(\Phi - \lambda)\right) \sin \frac{\Theta}{2} \\ \exp \left(\frac{i}{2}(\Phi - \lambda)\right) \sin \frac{\Theta}{2} & \exp \left(\frac{i}{2}(\Phi + \lambda)\right) \cos \frac{\Theta}{2} \end{pmatrix}$

Other gates can then be built from this family of one qubit gates. When it comes to multi-qubit gates, the only multi-qubit gate specified by OpenQASM is the CNOT gate denoted by CX, which can then again be combined with other gates to obtain gates operating on three and more qubits.

For qiskit, the available gates are specified in QASM syntax in the file qelib1.inc (see https://github.com/Qiskit/qiskit-terra/blob/master/qiskit/qasm/libs/qelib1.inc). Note that global phases are suppressed, so if you carry out the calculations, you will sometimes find a difference in the (unimportant) global phase between the result of the calculation and the results in Qiskit.

There is a couple of gates that you will often use in your circuits and that are summarized in the following table.

Gate Description
X Pauli X gate
Y Pauli Y gate
Z Pauli Z gate
S Phase gate $\text{diag}(1,i)$
T T gate $\text{diag}(1,e^{i\frac{\pi}{4}})$
CX CNOT gate

Gates take arguments that specify the qubits on which the gates operate. Individual qubits in a register can be addressed using an array-like notation. For example, to implement a circuit that applies an X gate to the first (least significant) qubit and then a controlled-NOT gate with this qubit as control qubit and the second qubit as target qubit, the following code can be used.

q = QuantumRegister(2,"q")
c = ClassicalRegister(2,"c")
circuit = QuantumCircuit(q,c)
circuit.x(q[0])
circuit.cx(q[0], q[1])


On real hardware, CNOTs can only be applied to specific combinations of qubits as specified in the coupling map of the device. However, when a circuit is prepared for execution on a specific device – a process called compilation – the compiler will deal with that by either reordering the qubits or adding additional swap operations.

Now we have gates, but to be able to run the circuit and measure the outputs, we still need measurements. These can easily been added with the measure method of a circuit, which accepts two parameters – the quantum register to measure and the corresponding classical register.

circuit.measure(q,c)


When the measurement step is reached during the processing of the circuit, the measurement is done – resulting in the projection of the state vector to the corresponding subspace – and the results of the measurements are copied into the classical register from which they can then be retrieved.

A nice property of Qiskit is its ability to visualize a quantum circuit. For that purpose, several classes called drawers are available in qiskit.tools.visualization. The circuit above, for instance, can be plotted with only one command

from qiskit.tools.visualization import matplotlib_circuit_drawer as drawer
my_style = {'cregbundle': True}
drawer(circuit, style=my_style)


and gives the nice representation

## Compiling and running circuits

Let us now actually run the circuit. To do this, we need a Qiskit backend. Qiskit offers several backends. The Aer package contains a few simulators that are installed locally and can be executed directly from a Python script or a notebook. Some of these simulators calculate the actual state vectors (the unitary simulator and the state vector simulator), but cannot deal with measurements, others – the QASM simulator – only provide statistical results but can simulate the entire circuit including measurements.

The IBMQ package can be used to connect to the devices offered by the IBM Q experience program, including an online simulator with up to 32 qubits and the actual devices. As for the composer, accessing the IBM Q experience devices does obviously require an account and available units.

In order to run a circuit, we first compile the circuit, which will create a version of the circuit that is tailored for the specific hardware, and then submit the circuit as a job.

backend = IBMQ.get_backend('ibmq_16_melbourne')
from qiskit import compile
qobj = compile(circuit, backend=backend, shots=1024)
job = backend.run(qobj)


Once the job has been submitted, we can poll its status using job.status() every few seconds. When the job has completed, we can access the results using job.result(). Every job consists of a certain number of shots, i.e. individual executions, and the method result.get_counts() will return a hash map that lists the measured outcomes along with how often that outcome was obtained. The following gist shows a basic Python script that assembles a circuit, compiles it, submits a job to the Q experience and prints the results.

Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.

A few more features of the Qiskit package, including working with different simulators and visualization options as well as QASM output, are demonstrated in this script on my GitHub page. In one of the next posts, we will try to implement a real quantum algorithm, namely the Deutsch-Jozsa algorithm, and run it on IBMs device.

## References

1. Andrew W. Cross, Lev S. Bishop, John A. Smolin, Jay M. Gambetta, Open Quantum Assembly Language, arXiv:1707.03429v2 [quant-ph]
2. The Qiskit tutorial on basic quantum operations
3. The Qiskit documentation

## Shor’s quantum factoring algorithm

Until the nineties of the last century, quantum computing seemed to be an interesting theoretical possibility, but it was far from clear whether it could be useful to tackle computationally hard problems with high relevance for actual complications. This changed dramatically in 1994, when the mathematician P. Shor announced a quantum algorithm that could efficiently solve one of the most intriguing problems in applied mathematics – factoring large numbers into their constituent primes, which, for instance, can be used to break commonly used public key cryptography schemes like RSA.

Shor’s algorithm is significantly more complicated than the quantum algorithms that we have studied so far, so we start with a short overview and then look at the individual pieces in more detail.

Given a number M that we want to factorize, the first part of Shor’s algorithm is to find a number x which has no common divisor with M so that it is a unit modulo M. In practice, we can just guess some x and compute the greatest common divisor gcd(x,M) – if this is not one, we have found a factor of M and are done, if this is one, we have the number x that we need. This step can still be done efficiently on a classical computer and does not require a quantum computer.

The next part of the algorithm now uses a quantum algorithm to determine the period of x. The period is the smallest non-zero number r such that

$x^r \equiv 1 \mod M$

The core of this step is the quantum algorithm that we will study below. However, the quantum algorithm does not exactly return the number r, but it returns a number s which is close to a multiple of $\frac{N}{r}$, where N is a power of two. Getting r out of this information is again a classical step that uses the theory of continued fractions. The number N that appears here is N = 2n where n is the number of qubits that the quantum part requires, and needs to be chosen such that M2 can be represented with n bits.

Finally, the third part of the algorithm uses the period r to find a factor of M, which is again done classically using elementary number theory. Thus the overall layout of the algorithm is as follows.

• Find the smallest number n (the number of qubits that we will need) such that $M^2 \leq N = 2^n$ and a number $x < M$ such that gcd(x,M) = 1
• Use the quantum part of the algorithm to find a number s which is approximately an integer multiple of N / r
• Use the theory of continued fractions to extract the period r from this information
• Use the period to find a factor of M

We will know look at each of these steps in turn (yes, this is going to be a bit of a lengthy post). To make this more tangible, we use a real example and assume that we wanted to factor the number M = 21. This is of course a toy example, but it allows us to simulate and visualize the procedure efficiently on a classical computer.

## Determine n and x

The first step is easy. First, we need to determine the number n of qubits that we need. As mentioned above, this is simply the bit length of M2. In our case, M2 = 441, and the next power of two is 512 = 29, so we need n = 9 qubits.

The next step is to find the number x. This can easily be done by just randomly picking some x and checking that is has no common prime factor with M. In our example, let us choose x = 11 (which is a prime number, but this is just by accident and not needed). It is important to choose this rather randomly, as the algorithm might fail in some rare instances and we need to start over, but this only makes sense if we do not pick the same choice again for our second trial.

## The quantum part of the algorithm

Now the quantum part of the algorithm starts. We want to calculate the period of x = 11, i.e. the smallest number r such that xr – 1 is a multiple of M = 21.

Of course, as our numbers are small, we could easily calculate the period classically by taking successive powers of 11 and reducing modulo 21, and this would quickly tell us that the period is 6. This, however, is no longer feasible with larger numbers, and this is where our quantum algorithm comes into play.

The algorithm uses two quantum registers. The first register has n qubits, the second can have less, in fact any number of qubits will do as long as we can store all numbers up to M in it, i.e. the bit length of M will suffice. Initially, we bring the system into the superposition

$\frac{1}{\sqrt{N}} \sum_k |k \rangle |0 \rangle$

which we can for instance do by starting with the state with all qubits being zero and then applying the Hadamard-Walsh transformation to the first register.

Next, we consider the function f that maps a number k to xk modulo M. As for every classical function, we can again find a quantum circuit Uf that represents this function on the level of qubits and apply it to our state to obtain the state

$\frac{1}{\sqrt{N}} \sum_k |k \rangle |x^k \mod M \rangle$

In his original paper [2], Shor calls this part the modular exponentiation and shows that this is actually the part of the quantum algorithm where most gates are needed (not the quantum Fourier transform).

This state has already some periodicity built into it, as xk modulo M is periodic with period r. If we could measure all the amplitudes, we could easily determine r. However, every such measurement destroys the quantum state and we have to start again, so this algorithm will not be very efficient. So again, the measurement is an issue.

Now, Shor’s idea is to solve the measurement issue by first applying (the inverse of) a quantum Fourier transform to the first register and then measure the first register (we apply the inverse of the quantum Fourier transform while other sources will state that the algorithm uses the quantum Fourier transform itself, but this is just a matter of convention as to which transformation you call the Fourier transform). The outcome s of this measurement will then give us the period!

To get an idea why this is true, let us look at a simpler case. Assume that, before applying the quantum Fourier transform, we measure the value of the second register. Let us call this value y. Then, we can write y as a power of x modulo M. Let k0 be the smallest exponent such that

$x^{k_0} = y$

Then, due to the periodicity, all values of k such that xk = y modulo M are given by

$k = k_0 + jr$

Here the index j needs to be chosen such that k0 + jr is still smaller than M. Let A denote the number of possible choices for j. Then, after the measurement, our state will have collapsed to

$\frac{1}{\sqrt{A}} \sum_{j=0}^{A-1} |k_0 + jr \rangle |y \rangle$

Let us now apply the inverse of the quantum Fourier transform to this state. The result will be the state

$\frac{1}{\sqrt{AN}} \sum_{j=0}^{A-1} \sum_{s=0}^{N-1} \eta^{(x_0 + jr)s} |s \rangle$

Now let us measure the first register. From the expression above, we can read off the probability P(s) to measure a certain value of s – we just have to add up the squares of all amplitudes with this value of s. This gives us

$P(s) = \frac{1}{AN} \big| \sum_{j=0}^{A-1} \eta^{jrs} \big|^2$

This looks complicated, but in fact this is again a geometric series with coefficient $q = \eta^{rs}$. To see how the value of the series depends on s, let us assume for a moment that the period divides N (which is very unlikely in practice as N is a power of two, but let us assume this anyway just for the sake of argument), i.e. that N = r u with the frequency u being an integer. Thus, if s is a multiple of u, the coefficient q is equal to one (as $\eta^N = 1$) and the geometric series sums up to A, giving probability 1 / N to measure this value. If, however, s is not a multiple of u, the value of the geometric series is

$\frac{1 - q^A}{1 - q}$

But in our case, A is of course simply equal to u, and therefore qA is equal to one. Thus the amplitude is zero! We find – note the similarity to our analysis of the Fourier transform of a periodic sequence – that P(s) is sharply peaked at multiples of $u = \frac{N}{r}$!

We were able to derive this result using a few simplifications – an additional measurement and the assumption that the frequency is an integer. However, as carried out by Shor in [2], a careful analysis shows that these assumptions are not needed. In fact, one can show (if you want to see all the nitty-gritty details, you could look at Shor’s paper or at my notes on GitHub that are based on an argument that I have seen first in Preskill’s lecture notes) that with reasonably high probability, the result s of the measurement will be such that

$\big| \{sr\}_N \big| \leq \frac{r}{2}$

where $\{sr\}_N$ denotes the residual of sr modulo N. Intuitively, this means that with high probability, the residual is very small, i.e. rs is close to a multiple of N, i.e. s is close to a multiple of N / r. In other words, it shows that in fact, P(s) has peaks at multiples of N / r.

The diagram below plots the probability distribution P(s) for our example, i.e. N = 512 and r = 6 (this plot has been generated using the demo Shor.py available in my GitHub account which uses the numpy package to simulate a run of Shor’s algorithm on a classical computer)

As expected, we see sharp peaks, located roughly at multiples of 512 / 6 = 85.33. So when we measure the first register, the value s will be close to a multiple of 512 / 6 with a very high probability.

So the quantum algorithm can be summarized as follows.

• Prepare a superposition $\frac{1}{\sqrt{N}} \sum_k |k \rangle |x^k \mod M \rangle$
• Apply the (inverse of the) quantum Fourier transform to this state
• Measure the value of the first register and call the result s

When running the simulation during which the diagram above was created, I did in fact get a measurement at s = 427 which is very close to 5*512 / 6.

## Extracting the period

So having our measurement s = 427 in our hands, how can we use this to determine the period r? We know from the considerations above that s is close to a multiple of N / r, i.e. we know that there is an integer d such that

$\big| sr - dN \big| \leq \frac{r}{2}$

which we can rewrite as

$\big| \frac{s}{N} - \frac{d}{r} \big| \leq \frac{1}{2N} < \frac{1}{M^2}$

Thus we are given two rational numbers – s / N and d / r – which we know to be very close to each other. We have the first number s / N in our hands and want to determine the second number. We also know that the denominator r of the second number is smaller than M. Is this sufficient to determine d and r?

Surprisingly, the answer is “yes”. We will not go into details at this point and gloss over some of the number theory, but refer the reader to the classical reference [1] or to my notes for more details). The first good news is that two different fractions with denominators less than M need to be at least by 1 / M2 apart, so the number d / r is unique. The situation is indicated in the diagram below.

But how to find it? Luckily, the theory of continued fractions comes to the rescue. If you are not familiar with continued fractions, you can find out more in the appendix of my notes or on the very good Wikipedia page on this. Here, we will just go through the general procedure using our example at hand.

First, we write

$\cfrac{427}{512} = 0 + \cfrac{1}{\cfrac{512}{427}} = 0 + \cfrac{1}{1 + \cfrac{85}{427}}$

We can do the same with 85 / 427, i.e. we can write

$\cfrac{85}{427} = \cfrac{1}{\cfrac{427}{85}} = \cfrac{1}{5 + \cfrac{2}{85}}$

which will give us the decomposition

$\cfrac{427}{512} = 0 + \cfrac{1}{1 + \cfrac{1}{5 + \cfrac{2}{85}}}$

Driving this one step further, we finally obtain

$\cfrac{427}{512} = 0 + \cfrac{1}{1 + \cfrac{1}{5 + \cfrac{1}{42 + \cfrac{1}{2}}}}$

This is called the continued fraction expansion of the rational number 427 / 512. More generally, for every sequence $[a_0 ; a_1, a_2, \dots ]$, we can form the continued fraction

$a_0 + \cfrac{1}{a_1 + \cfrac{1}{a_2 + \dots}}$

given by that sequence of coefficients, which is obviously a rational number. One can show that in fact every rational number has a representation as a continued fraction, and our calculation has shown that

$\cfrac{427}{512} = [0; 1,5,42,2]$

This sequence has five coefficients. Now given a number m, we can of course look at the sequence that we obtain by looking at the cofficients up to index m only. For instance, for m = 3, this would give us the sequence

$[0; 1,5,42]$

The rational number represented by this sequence is

$0 + \cfrac{1}{1 + \cfrac{1}{5 + \cfrac{1}{42}}} = \cfrac{211}{253}$

and is called the m-th convergent of the original continued fraction. We have such a convergent for every m, and thus get a (finite) series of rational numbers with the last one being the original number.

One can now show that given any rational number x, the series of m-th convergents of its continued fraction expansion has the following properties.

• The convergents are in their lowest terms
• With increasing m, the difference between x and the m-th convergent gets smaller and smaller, i.e. the convergents form an approximation of x that gets better and better
• The denominators of the convergents are increasing

So the convergents can be used to approximate the rational number x by fractions with smaller denominator – and this is exactly what we need: we wish to approximate the rational number s / N by a ratio d / r with smaller denominator which we then know to be the period. Thus we need to look at the convergents of 427 / 512. These can be easily calculated and turn out to be

$0, 1, \cfrac{5}{6}, \cfrac{211}{253}, \cfrac{427}{512}$

The last convergent whose denominator is still smaller than M = 21 is 5 / 6, and thus we obtain r = 6. This is the period that we are looking for!

So in general, the recipe to get from the measured value s to r is to calculate the convergents of the rational number s / N and pick the denominator of the last convergent that has a denominator less than M. Again, if you want to see the exact algorithm, you can take a look at my script Shor.py.

## Find the factor

We are almost done. We have run the quantum algorithm to obtain an approximate multiple of N / r. We have then applied the theory of continued fractions to derive the period r of x from this measurement. The last step – which is again a purely classical step – is now to use this to find a factor of M. This is in fact comparatively easy.

Recall that – by definition of the period – we get one if we raise x to the power of r and than reduce module M. In other words, xr minus one is a multiple of M. Now assume that we are lucky and the period r is even. Then

$(x^{\frac{r}{2}} - 1)(x^{\frac{r}{2}} + 1) = (x^r - 1) \equiv 0 \mod M$

With a bit of elementary number theory, one can now show that this implies that the greatest common divisor gcd(xr/2-1, M) is a factor of M (unless, in fact, xr/2 is minus one modulo M, in which case the argument fails). So to get from the period to a potential factor of M, we simply calculate this greatest common divisor and check whether it divides M!

Let us do this for our case. Our period is r = 6. With x = 11, we have x3 = 1331, which is 8 module M. Thus

$\text{gcd}(x^{\frac{r}{2}} - 1 \mod M, M) = \text{gcd}(7,21) = 7$

which is the factor of M = 21 that we were looking for.

## Performance of the algorithm

In our derivation, we have ignored a few special cases which can make the algorithm fail. For instance, suppose we had not measured s = 427, but s = 341 after applying the Fourier transform. Then the corresponding approximation to 341 / 512 would have been 4 / 6. However, the continued fraction algorithm always produces a result that is in its lowest terms, i.e. it would give us not 4 / 6, but 2 / 3. Looking at this, we would infer that r = 3, which is not the correct result.

There are a few other things that can go wrong. For instance, we could find a period r which is odd, so that our step to derive a factor of M from r does not work, or we might measure an unlikely value of s.

In all these cases, we need to start over and repeat the algorithm. Fortunately, Shor could show that the probability that any of this happens is bounded from below. This bound is decreasing with larger values of M, but it is decreasing so slowly that the expected number of trials that we need grows at most logarithmically and does not destroy the overall performance of the algorithm.

Taking all these considerations into account and deriving bounds for the number of gates required to perform the quantum part of the algorithm, Shor was able to show that the number of steps to obtain a result grows at most polynomially with the number of bits that the number M has. This is obviously much better than the best classical algorithm that requires a bit less than exponential time to factor M. Thus, assuming that we are able to build a working quantum computer with the required number of gates and qubits, this algorithm would be able to factorize large numbers exponentially faster than any known classical algorithm.

Shor’s algorithm provides an example for a problem that is believed to be in the class NP (but not in P) on a classical computer, but in the class BQP on a quantum computer – this is the class of all problems that can be solved in polynomial time with a finite probability of success. However, even though factorization is generally believed not to be in P, i.e. not doable in polynomial time on classical hardware, there is not proof for that. And, even more important, it is not proved that factorization is NP-complete. Thus, Shor’s algorithm does not render every problem in NP solvable in polynomial time on a quantum computer. It does, however, still imply that all public key cryptography systems like RSA that rely on the assumption that large numbers are difficult to factor become inherently insecure once a large scale reliable quantum computer becomes available.

## References

1. G.H. Hardy, E.M. Wright, An introduction to the theory of numbers, Oxford University Press, Oxford, 1975
2. P. Shor, Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer, SIAM J.Sci.Statist.Comput. Vol. 26 Issue 5 (1997), pp 1484–1509, available as arXiv:quant-ph/9508027v2

## Grover’s algorithm – unstructured search with a quantum computer

In the last post, we have looked at the Deutsch-Jozsa algorithm that is considered to be the first example of a quantum algorithm that is structurally more efficient than any classical algorithm can probably be. However, the problem solved by the algorithm is rather special. This does, of course, raise the question whether a similar speed-up can be achieved for problems that are more relevant to practical applications.

In this post, we will discuss an algorithm of this type – Grover’s algorithm. Even though the speed-up provided by this algorithm is rather limited (which is of a certain theoretical interest in its own right), the algorithm is interesting due to its very general nature. Roughly speaking, the algorithm is concerned with an unstructured search. We are given a set of N = 2n elements, labeled by the numbers 0 to 2n-1, exactly one of which having a property denoted by P. We can model this property as a binary valued function P on the set $\{0,N-1\}$ that is zero on all but one elements. The task is to locate the element x0 for which P(x0) is true.

## Grover’s algorithm

Grover’s algorithm presented in [1] proceeds as follows to locate this element. First, we again apply the Hadamard-Walsh operator W to the state $|0 \rangle$ of an n-qubit system to obtain a superposition of all basis states. Then, we iteratively apply the following sequence of operations.

1. Apply a conditional phase shift S, i.e. apply the unique unitary transformation that maps $|x \rangle$ to $(-1)^{f(x)} |x \rangle$.
2. Apply the unitary transformation D called diffusion that we will describe below

Finally, after a defined number of outcomes, we perform a measurement which will collaps the system into one of the states $|x \rangle$. We claim – and will see why this is true below – that for the right number of iterations, this value x will, with a high likelihood, be the solution to our problem, i.e. equal to x0.

Before we can proceed, we need to define the matrix D. This matrix is $\frac{2}{N} - 1$ along the diagonal, with N = 2n, and $\frac{2}{N}$ away from the diagonal. In terms of basis vectors, the mapping is given by

$|i \rangle \mapsto \frac{2}{N} (\sum_j |j \rangle) - |i \rangle$

Consequently, we see that

$D \sum_i a_i |i \rangle = \sum_i (2 \bar{a} - a_i) |i \rangle$

where $\bar{a}$ is the average across the amplitudes $a_i$. Thus geometrically, the operation D performs an inversion around the average. Grover shows that this operation can be written as minus a Hadamard-Walsh operation followed by the operation that flips the sign for $|0 \rangle$, followed by a second Hadamard-Walsh transformation.

For the sake of completeness, let us also briefly discuss the first transformation employed by the algorithm, the conditional phase shift. We have already seen a similar transformation while studying the Deutsch-Jozsa algorithm. In fact, we have shown in the respective blog post that the circuit displayed below (with the notation slightly changed)

performs the required operation

$|\psi \rangle = \sum_x a_x |x \rangle \mapsto |\psi' \rangle = \sum_x a_x (-1)^{P(x)} |x \rangle$

Let us now see how why Grover’s algorithm works. Instead of going through the careful analysis in [1], we will use bar charts to visualize the quantum states (exploiting that all involved matrices are actually real valued).

It is not difficult to simulate the transformation in a simple Python notebook, at least for small values of N. This script performs several iterations of the algorithm and prints the result. The diagrams below show the outcome of this test.

Let us go through the diagrams one by one. The first diagram shows the initial state of the algorithm. I have used 3 qubits, i.e. n = 3 and N = 8. The initial state, after applying the Hadamard-Walsh transform to the zero state, is displayed in the first line. As expected, all amplitudes are equal to 1 over the square root of eight, which is approximately 0.35, i.e. we have a balanced superposition of all states.

We now apply one iteration of the algorithm. First, we apply the conditional phase flip. The element we are looking for is in this case located at x = 2. Thus, the phase flip will leave all basis vectors unchanged except for $|2 \rangle$ and it will change the amplitude of this vector to – 0.35. This will change the average amplitude to a value slightly below 0.35. If we now perform the inversion around the average, the amplitudes of all basis vectors different from $|2 \rangle$ will actually decrease, whereas the amplitude of $|2 \rangle$ will increase. The result is displayed in the second line of the diagram.

Thus, what really happens in this case is an amplitude amplification – we increase the amplitude of one component of the superposition while decreasing all the others.

The next few lines show the result of repeating these two steps. We see that after the second iteration, almost all of the amplitude is concentrated on the vector $|2 \rangle$, which represents the solution we are looking for. If we now perform a measurement, the result will be 2 with a very high probability!

It is interesting to see that when we perform one more iteration, the difference between the amplitude of the solution and the amplitudes of all other components decreases again. Thus the correct choice for the number of iterations is critical to make the algorithm work. In the last line, we have plotted the difference between the amplitude of $|2 \rangle$ and the other amplitudes (more precisely, the ratio between the amplitude of $|2 \rangle$ and the second largest amplitude) on the y-axis for the different number of iterations on the x-axis. We see that the optimal number of iterations is significantly below 10 (actually five iterations give the best result in this case), and more iterations decrease the likelihood of getting the correct result out of the measurement again. In fact, a careful analysis carried out in [2] shows that for large values of N, the best number of iterations is given by $\frac{\pi}{4} \sqrt{N}$, and that doubling the number of iterations does in general lead to a less optimal result.

## Generalizations and amplitude amplification

In a later paper ([3]), Grover describes a more general setup which is helpful to understand the basic reason why the algorithm works – the amplitude amplification. In this paper, Grover argues that given any unitary transformation U and a target state (in our case, the state representing the solution to the search problem), the probability to meet the target state by applying U to a given initial state can be amplified by a sequence of operations very much to the one considered above. We will not go into details, but present a graphical representation of the algorithm.

So suppose that we are given an n-qubit quantum system and two basis vectors – the vector t representing the target state and an initial state s. In addition, assume we are given a unitary transformation U. The goal is to reach t from s by subsequently applying U itself and a small number of additional gates.

Grover considers the two-dimensional subspace spanned by the vectors s and U-1t. If, within this subspace, we ever manage to reach U-1t, then of course one more application of U will move us into the desired target state.

Now let us consider the transformation

$Q = - I_s U^{-1} I_t U$

where Ix denotes a conditional phase shift that flips the phase on the vector $|x \rangle$. Grover shows that this transformation does in fact leave our two-dimensional subspace invariant, as indicated in the diagram below.

He then proceeds to show that for sufficiently small values of the matrix element Uts, the action of Q on this subspace can approximately be described as a rotation by the angle $\frac{2}{|U_{ts}|}$. Applying the operation Q n times will then approximately result in a superposition of the form

$\cos (\frac{2n}{|U_{ts}|})|s \rangle + a U^{-1}|t \rangle$

Thus if we can make the first coefficient very small, i.e. if $\frac{4n}{|U_{ts}|}$ is close to a multiple of $\pi$, then one application of U will take our state to a state very close to t.

Let us link this description to the version of the Grover algorithm discussed above. In this version, the initial state s is the state $|0 \rangle$. The transformation U is the Hadamard-Walsh transformation W. The target state is $|x_0 \rangle$ where x0 is the solution to the search problem. Thus the operation It is the conditional phase shift that we have denoted by S earlier. In addition, Grover shows in [1] already that the diffusion operator D can be expressed as -W I0 W. Now suppose we apply the transformation Q n times to the initial state and then apply U = W once more. Then our state will be

$W Q \dots Q |0 \rangle$

This can be written as

$W (-I_0 W S W) (-I_0 W S W) \dots (-I_0 W S W) |0 \rangle$

Regrouping this and using the relation D = -W I0 W, we see that this is the same as

$- (W I_0 W) S (- W I_0 W) S \dots \dots (- W I_0 W) S W |0 \rangle = D S \dots D S W |0 \rangle$

Thus the algorithm can equally well be described as applying W once to obtain a balanced superposition and then applying the sequence DS n times, which is the formulation of the algorithm used above. As $|U_{ts} | = \frac{1}{\sqrt{N}}$ in this case, we also recover the result that the optimal number of iterations is $\frac{\pi}{4} \sqrt{N}$ for large N.

## Applications

Grover’s algorithm is highly relevant for theoretical reasons – it applies to a very generic problem and (see the discussions in [1] and [2]) is optimal, in the sense that it provides a quadratic speedup compared to the best classical algorithm that requires O(N) operations, and that this cannot be improved further. Thus Grover’s algorithm provides an interesting example for a problem where a quantum algorithm delivers a significant speedup, but no exponential speedup as we will see it later for Shor’s algorithm.

However, as discussed in detail in section 9.6 of [4], the relevance of the algorithm for practical applications is limited. First, the algorithm applies to an unstructured search, i.e. a search over unstructured data. In most practical applications, we deal with databases that have some sort of structure, and then more efficient search algorithms are known. Second, we have seen that the algorithm requires $O(\sqrt{N})$ applications of the transformation UP. Whether this is better than the classical algorithm does of course depend on the efficiency with which we can implement this with quantum gates. If applying UP requires O(N) operations, the advantage of the algorithm is lost. Thus the algorithm only provides a significant speedup if the operation UP can be implemented efficiently and there is no additional structure that a classical algorithm could exploit.

Examples of such problems are brute forces searches as they appear in some crypto-systems. Suppose for instance we are trying to break a message that is encrypted with a symmetric key K, and suppose that we know the first few characters of the original text. We could then try to use an unstructured search over the space of all keys to find a key which matches at least the few characters that we know.

In [5], a more detailed analysis of the complexity in terms of qubits and gates that a quantum computer would have to attack AES-256 is made, arriving at a size of a few thousand logical quantum bits. Given the current ambition level, this does not appear to be completely out of reach. It does, however, not render AES completely unsecure. In fact, as Grover’s algorithm essentially results in a quadratic speedup, a code with a key length of n bits in a pre-quantum world is essentially as secure as the same code with a key length of 2n in a post-quantum world, i.e. roughly speaking, doubling the key length compensates the advantage of quantum computing in this case. This is the reason why the NIST report on post-quantum cryptography still classifies AES as inherently secure assuming increased key sizes.

In addition, the feasibility of a quantum algorithm is not only determined by the number of qubits required, but also by other factors like the depth, i.e. the number of operations required, and the number of quantum gates – and for AES, the estimates in [5] are significant, for instance a depth of more than 2145 for AES-256, which is roughly 1043. Even if we assume a switching time of only 10-12 seconds, we still would require astronomical 1031 seconds, i.e. in the order of 1023 years, to run the algorithm.

Even a much less sophisticated analysis nicely demonstrates the problem behind these numbers – the number of iterations required. Suppose we are dealing with a key length of n bits. Then we know that the algorithm requires

$\frac{\pi}{4} \sqrt{2^n} \approx 0.8 \sqrt{2}^n$

iterations. Taking the decimal logarithm, we see that this is in the order of 100.15*n. Thus, for n = 256, we need in the order of 1038 iterations – a number that makes it obvious that AES-256 can still be considered secure for all practical purposes.

So overall, there is no reason to be overly concerned about serious attacks to AES with sufficiently large keys in the near future. For asymmetric keys, however, we will soon see that the situation is completely different – algorithms like RSA or Elliptic curve cryptography are once and for all broken as soon as large-scale usable quantum computer become reality. This is a consequence of Shor’s algorithm that we will study soon. But first, we need some more preliminaries that we will discuss in the next post, namely quantum Fourier transforms.

## References

1. L.K. Grover, A fast quantum mechanical algorithm for database search, Proceedings, 28th Annual ACM Symposium on the Theory of Computing (STOC), May 1996, pages 212-219, available as arXiv:quant-ph/9605043v3
[2] M. Boyer, G. Brassard, P. Høyer, A. Tapp, Tight bounds on quantum searching, arXiv:quant-ph/9605034
3. L.K. Grover, A framework for fast quantum mechanical algorithms, arXiv:quant-ph/9711043
4. E. Rieffel, W. Polak, Quantum computing – a gentle introduction, MIT Press
5. M. Grassl, B. Langenberg, M. Roetteler, R. Steinwandt, Applying Grover’s algorithm to AES: quantum resource estimates, arXiv:1512.04965

## More on Paperspace Gradient

Its been a few days since I started to play with Paperspace, and I have come across a couple of interesting features that the platform has – enough for a second post on this topic.

First, GIT integration. Recall that the usual process is to zip the current working directory and submit the resulting file along with the job, the ZIP file is then unzipped in the container in which the job is running and the contents of the ZIP file constitute the working directory. However, if you want to run code that requires, for instance, custom libraries, it is much easier to instruct Paperspace to get the contents of the working directory from GitHub. You can do that by supplying a GIT URL using the --workspace switch. The example below, for instance, instructs Paperspace to pull my code for an RBM from GitHub and to run it as a job.

#!/bin/sh
#
# Run the RBM as a job on Paperspace. Assume that you have the paperspace NodeJS
# CLI and have done a paperspace login before to store your credentials
#
#
~/node_modules/.bin/paperspace jobs create  \
--workspace "git+https://github.com/christianb93/MachineLearning" \
--command "export MPLBACKEND=AGG ; python3 RBM.py \
--N=28 --data=MNIST \
--save=1 \
--tmpdir=/artifacts \
--hidden=128 \
--pattern=256 --batch_size=128 \
--epochs=40000 \
--run_samples=1 \
--sample_size=6,6 \
--beta=1.0 --sample=200000 \
--algorithm=PCDTF --precision=32" \
--machineType K80 \
--container "paperspace/tensorflow-python" \
--project "MachineLearning"


Be careful, the spelling of the URL must be exactly like this to be recognized as a GIT URL, i.e. “git+https” followed by the hostname without the “www”, if you use http instead of https or http://www.github.com instead of github.com, the job will fail (the documentation at this point could be better, and I have even had to look at the source code of the CLI to figure out the syntax). This is a nice feature, using that along with the job logs, I can easily reconstruct which version of the code has actually been executed, and it supports working in a team that is sharing GitHub repositories well.

Quite recently, Paperspace did apparently also add the option to use persistent storage in jobs to store data across job runs (see this announcement). Theoretically, the storage should be shared between notebooks and jobs in the same region, but as I have not yet found out how to start a notebook in a specific region, I could not try this out.

Another feature that I liked is that the container that you specify can actually be any container from the Docker Hub, for instance ubuntu. The only restriction is that Paperspace seems to overwrite the entrypoint in any case and will try to run bashinside the container to finally execute the command that you provide, so containers that do not have a bash in the standard execution path will not work. Still, you could use this to prepare your own containers, maybe with pre-installed data sets or libraries, and ask Paperspace to run them.

Finally, for those of us who are Python addicts, there is also a Python API for submitting and managing jobs in Paperspace. Actually, this API offers you two ways to run a Python script on Paperspace. First, you can import the paperspace package into your script and then, inside the script, do a paperspace.run(), as in the following example.

import paperspace
paperspace.run()
print('This will only be running on Paperspace')


What will happen behind the scenes is that the paperspace module takes your code, removes any occurrences of the paperspace package itself, puts the code into a temporary file and submits that as a job to Paperspace. You can then work with that job as with any other job, like monitoring it on the console or via the CLI.

That is nice and easy, but not everyone likes to hardcode the execution environment into the code. Fortunately, you can also simply import the paperspace package and use it to submit an arbitrary job, much like the NodeJs based CLI can do it. The code below demonstrates how to create a job using the Python API and download the output automatically (this script can also be found on GitHub).

import paperspace.jobs
from paperspace.login import apikey
import paperspace.config

import requests

#
# Define parameters
#
params = {}
#
# We want to use GIT, so we use the parameter workspaceFileName
# instead of workspace
#
params['workspaceFileName'] = "git+https://github.com/christianb93/MachineLearning"
params['machineType'] = "K80"
params['command'] = "export MPLBACKEND=AGG ; python3 RBM.py \
--N=28 --data=MNIST \
--save=1 \
--tmpdir=/artifacts \
--hidden=128 \
--pattern=256 --batch_size=128 \
--epochs=40000 \
--run_samples=1 \
--sample_size=6,6 \
--beta=1.0 --sample=200000 \
--algorithm=PCDTF --precision=32"
params['container'] = 'paperspace/tensorflow-python'
params['project'] = "MachineLearning"
params['dest'] = "/tmp"

#
# Get API key
#
apiKey = apikey()
print("Using API key ", apiKey)
#
# Create the job. We do NOT use the create method as it cannot
# handle the GIT feature, but assemble the request ourselves
#
http_method = 'POST'
path = '/' + 'jobs' + '/' + 'createJob'

r = requests.request(http_method, paperspace.config.CONFIG_HOST + path,
params=params, files={})
job = r.json()

if 'id' not in job:
print("Error, could not get jobId")
print(job)
exit(1)

jobId = job['id']
print("Started job with jobId ", jobId)
params['jobId']  = jobId

#
# Now poll until the job is complete

if job['state'] == 'Pending':
print('Waiting for job to run...')
job = paperspace.jobs.waitfor({'jobId': jobId, 'state': 'Running'})

print("Job is now running")
print("Use the following command to observe its logs: ~/node_modules/.bin/paperspace jobs logs --jobId ", jobId, "--tail")

job = paperspace.jobs.waitfor({'jobId': jobId, 'state': 'Stopped'})
print("Job is complete: ", job)

#
# Finally get artifacts
#