The Shor Quantum Error Correcting Code (and a Monad for Heat)

Quantum computers promise great things, but that promise will only come to be if quantum computers can be made reliable. In order to achieve that reliability we need qubits that can maintain their state for extended periods of time. Unfortunately, just about any type of interaction between a qubit and its environment can render the qubit useless. But there is a solution - a quantum error correcting code. This is very similar to classical error correcting code in which a number of bits are stored, in a distributed fashion, over a larger set of bits, in such a way that if some of the bits are corrupted, the original state can be recovered. In the following will be some code to simulate the implementation of one of these codes: the nine bit Shor code.

Firstly, let me include the quantum monad that I've developed over earlier weeks:

> import Prelude hiding (flip)
> import Maybe
> import List
> import Data.Map (toList,fromListWith)
> import Complex
> import Data.Bits hiding (rotate)

> infixl 1 ->-, -><, =>=
> infixl 7 .*

Remember that W b a is the type of elements of a vector space over b with a basis whose elements are labelled by elements of a. I'm using the letter W as I'm thinking of elements of b as being weights.

> data W b a = W { runW :: [(a,b)] } deriving (Eq,Show,Ord)

It's useful to be able to apply transforms to all of the weights:

> mapW f (W l) = W $ map (\(a,b) -> (a,f b)) l

W is an instance of Num, Functor and Monad. Vector spaces form monads because you can specify a function between vectors spaces in terms of how it maps basis elements in the first vector space. The monad bind function lifts the function on basis elements to the entire vector space.

> instance Functor (W b) where
>   fmap f (W a) = W $ map (\(a,p) -> (f a,p)) a

> instance Num b => Monad (W b) where
>   return x = W [(x,1)]
>   l >>= f = W $ concatMap (\(W d,p) -> map (\(x,q)->(x,p*q)) d) (runW $ fmap f l)

(.*) gives us scalar multiplication:

> a .* b = mapW (a*) b

> instance (Eq a,Show a,Num b) => Num (W b a) where
>    W a + W b = W $ (a ++ b)
>    a - b = a + (-1) .* b
>    _ * _ = error "Num is annoying"
>    abs _ = error "Num is annoying"
>    signum _ = error "Num is annoying"
>    fromInteger a = if a==0 then W [] else error "fromInteger can only take zero argument"

It'd be cool if the bind function for W could simplfy our vectors by collecting together multiples of the same basis elements. Unfortunately this is tricky to do for reasons described in Eric Kidd's recent Set monad article. Instead I'll use an explicit function, collect, that collects our terms together. I've also used a crude method for throwing out terms that probably ought to be zero. A better solution might be to work only with rational numbers, or exact algebraic numbers.

> nonTrivial x = 10+abs x/=10
> collect :: (Ord a,Num b) => W b a -> W b a
> collect = W . filter (nonTrivial . snd) . toList . fromListWith (+) . runW

And here is the quantum monad:

> type Q a = W (Complex Double) a

There are a number of ways to apply functions to vector spaces. I'll be writing functions on the right so that I can apply operations to a quantum state reading from left to right. So here's a bunch of combinators to implement this. First there's ordinary function application:

> (->-) :: a -> (a -> b) -> b
> x ->- f = f x

Then there are functions that operate on the vector space by simply relabelling the basis:

> (=>=) :: Q a -> (a -> b) -> Q b
> x =>= f = fmap f x

Then there are functions that map from Q a to Q b but are defined by how they act on the basis elements of Q a:

> (-><) :: Q a -> (a -> Q b) -> Q b
> g ->< f = g >>= f

As I mentioned last post on entanglement, we can form states for joint independent systems using the tensor product:

> (<*>) :: Q a -> Q b -> Q (a,b)
> x <*> y = do
>   a <- x
>   b <- y
>   return (a,b)

The tensor product is a lot like a multiplication table. The coefficient of (a,b) in x <*> y is the product of the coefficients of a and b in x and y respectively. If two systems are unentangled then you can simulate the physics of each part separately. So it'd be nice to have a way to test whether or not a given state is unentangled, in other words, to test whether a state is a tensor product. As I'm not developing a serious application here, here's a quick and dirty algorithm that more or less does what anyone would do when faced with the task of recovering the row and column headings from a multiplcation table:

> norm x = sqrt $ sum $ map (\(_,a) -> magnitude a^2) $ runW $ collect x

> untensor :: (Ord a,Ord b) => Q (a,b) -> (Q a,Q b)
> untensor (W a) =
>     let u = nub $ sort $ map (fst . fst) a
>         v = nub $ sort $ map (snd . fst) a
>         p = W $ map (\x -> (x,fromJust (lookup (x,head v) a))) u
>         q = W $ map (\x -> (x,fromJust (lookup (head u,x) a))) v
>         in ((1/norm p :+ 0) .* p,(1/norm q :+ 0) .* q)

untensor recovers x and y from x <*> y (up to multiplication by a complex number with unit norm) but it gives garbage if the argument isn't a tensor product. So we need to test whether or not the original vector is the tensor product of the inputs. So here's a test of how close the original is to the product of the putative factors:

> independence x = collect (x-uncurry (<*>) (untensor x))

We'll be doing lots of stuff with tensor products so it's useful to have some tools to apply functions to one or other factor:

> mapFst f (a,b) = (f a,b)
> mapSnd f (a,b) = (a,f b)
> mapFstM f (a,b) = do
>     a' <- a ->- f
>     return (a',b)
> mapSndM f (a,b) = do
>     b' <- b ->- f
>     return (a,b')

If a bit is a Bool then by rights, a qubit should be a Q Bool. But I think 0 and 1 look nicer so I'm using Q Integer instead. So here's the fundamentally quantum rotation operator:

> rotate :: Double -> Integer -> Q Integer
> rotate theta 1 = let theta' = theta :+ 0
>   in cos (theta'/2) .* return 1 - sin (theta'/2) .* return 0
> rotate theta 0 = let theta' = theta :+ 0
>   in cos (theta'/2) .* return 0 + sin (theta'/2) .* return 1

Although rotate theta is a map Integer -> Q Integer, using ->< it becomes a map Q Integer to Q Integer. Note that the Double argument isn't float data inside a quantum device, it's a hard-coded feature of the hardware itself. For example, if we choose theta to be pi/2 and -pi/2 we get the standard square root of not operation and its inverse:

> hnot = rotate (pi/2)
> unhnot = rotate (-pi/2)

(You may be wondering about the use of theta/2 in my definition of rotate. There's an interesting story behind that involving what I think is one of the weirdest features of physics. For another day.)

I'm also goint to manipulate lots of lists of qubits so here are some tools for that:

> mapAt i f x = let (h,t:ts) = splitAt i x in h ++ (f t : ts)

> mapAtM :: Monad m => Int -> (a -> m a) -> [a] -> m [a]
> mapAtM i f x = let (h,t:ts) = splitAt i x in do
>     y <- f t
>     return $ h ++ (y : ts)

> flip a = 1-a
> flipAt i = mapAt i flip

So here's a statement of the problem. We have a qubit, ie. something of type Q Integer. But quantum systems are very delicate and our qubit may wander into a different state. In particular, after a time t it will be as if rotate theta has been applied to it, except that we can't control theta. We need to find a way to recover the original qubit. In fact, it's worse than that: our qubit is in some environment which we can consider to be a list of qubits. So our state is given by an instance of Q ([Integer],Integer) and the interaction might not simply be a rotation of the qubit but an operation that entangles the qubit with the environment. This might happen if the qubit isn't perfectly isolated.

Let's forget about the full generality here and consider just one type of data corruption - a qubit that might or might not have been flipped. How can we make a device robust against such occurences? Well this is practically identical to the problem of making classical data robust against noise. We can simply use a 3 bit repetition code to store our qubits as three qubits. So here's an implementation of a function to repair a codeword from the classical 3 bit repetition code:

> repairRepetition' [b,c,d] = (if a<3 then mapAt a flip else id) [b,c,d] where
>     a = case (b,c,d) of
>         (0,0,0) -> 3
>         (1,1,1) -> 3
>         (0,0,1) -> 2
>         (1,1,0) -> 2
>         (0,1,0) -> 1
>         (1,0,1) -> 1
>         (1,0,0) -> 0
>         (0,1,1) -> 0

But there are a couple of problems with this. Firstly this is a many-to-one map because [0,0,0] and [0,0,1], say, both get mapped to [0,0,0]. In quantum mechanics, all of our maps need to be unitary, and in particular, invertible. A many-to-one map isn't invertible. Secondly, we're creating a register a, and again this is not an invertible operation. These two problems will actually cancel each other out because the apparent destruction of information that comes about from our many-to-one map can be fixed by allowing a to carry that information away. And we can ensure reversibility by using a supplied register rather than creating a new one:

> repairRepetition'' (a,[b,c,d]) = (a,(if a<3 then mapAt a flip else id) [b,c,d]) where
>     a = case (b,c,d) of
>         (0,0,0) -> a `xor` 3
>         (1,1,1) -> a `xor` 3
>         (0,0,1) -> a `xor` 2
>         (1,1,0) -> a `xor` 2
>         (0,1,0) -> a `xor` 1
>         (1,0,1) -> a `xor` 1
>         (1,0,0) -> a `xor` 0
>         (0,1,1) -> a `xor` 0

It's easy to write an inverse for this function so it's a valid unitary transform. But it's actually a bit tedious to write out all those xors. So we can cheat a bit. Let's just assume there's an infinite supply of input zeroes and use them at our leisure. We no longer need to write the xors, we just note that every time we assign a value to a new register the hardware engineer building your device installs some means of import zero qubits. The a register in the above function will never be touched again. But we can't simply throw it away. It's a kind of side effect and we'll call these bits <em>ancillae</em>. We'll have lots of these side effects that we need to collect up, so here's the final quantum repetition codeword repair code:

> repairRepetition (a,[b,c,d]) = (i:a,(if i<3 then mapAt i flip else id) [b,c,d]) where
>     i = case (b,c,d) of
>         (0,0,0) -> 3
>         (1,1,1) -> 3
>         (0,0,1) -> 2
>         (1,1,0) -> 2
>         (0,1,0) -> 1
>         (1,0,1) -> 1
>         (1,0,0) -> 0
>         (0,1,1) -> 0

There's one last thing we could do to this code: the pairs being handed in and out of this function are essentially just instances of the Writer monad. So we could just squirrel away our ancillae inside a monad. I won't do that here because I think it's good to see everything explicitly when figuring this stuff out for the first time. But think about the bits that we'd be writing out. They correspond to information, or entropy, that we'd ultimately like to pump out of our device as exhaust. In other words they are heat. So this monad would, in effect, be a heat monad. I found this amusing given what I wrote in reason number eight in my list from last year.

So far so good. But that only corrects one type of error: a bit flip. But note that it also protects superpositions from bit flips. So let's take a quantum state built from repetition codewords placed in an empty environment:

> state1 = return [] <*> (return [0,0,0] + return [1,1,1])

Now corrupt the last bit:

> state1' = state1 =>= mapSnd (mapAt 2 flip)

And repair it:

> state1'' = state1' =>= repairRepetition

Although we have repaired the state of our system, the state of the environment has changed. But that's fine because this state is still an unentangles tensor product of the environment and our system. You can see this just by looking at the state, but for peace of mind use independence state1'' to see that the deviation from being a tensor product is the zero vector.

But that's still only a narrow class of error that we can fix. Another type of error maps return 0 to return 0 but maps return 1 to -return 1. This is known as a phase flip. One way to view this as as a flip in a rotated basis:

> phaseFlip x = x ->- hnot ->- fmap flip ->< unhnot
> pFlipAtM i = mapAtM i phaseFlip

There's a nice strategy to protect against phase flips: rotate our data into a basis where phase flips become ordinary flips, and use the repetition code. So we'll use these two states to represent 0 and 1:

> plus,minus :: Q [Integer]
> plus = return [1,1,1] ->< mapM unhnot 
> minus = return [0,0,0] ->< mapM unhnot 

So let's make a state that is a superposition of these states in some environment:

> state2 = return [] <*> (plus + minus) ->- collect

Corrupt the middle bit with a phase flip:

> state2' = state2 ->< mapSndM (mapAtM 1 phaseFlip)

And repair it:

> state2'' = state2 ->< mapSndM (mapM hnot) =>= repairRepetition ->< mapSndM (mapM unhnot) ->- collect

Comparing state2 and untensor state2'' we see that the repair is successful. Unfortunately, we can use this scheme to repair phase flips, but not plain flips. We need a scheme that can repair both of these. So now we turn to Shor's ingenious solution: use both of these schemes as the same time. So we'll use the scheme that protects against phase flips with the twist that instead of using three qubits we use three triplets of qubits. We can use the three triplets to fix phase flips and we can use the three bits within a triplet to fix flips. We'll represent 0 and 1 using shor 0 and shor 1 as defined here:

> tripletPlus = (1/sqrt 2) .* return [0,0,0] + (1/sqrt 2) .* return [1,1,1]
> tripletMinus = (1/sqrt 2) .* return [0,0,0] - (1/sqrt 2) .* return [1,1,1]

> shor 1 = do
>     a <- tripletPlus
>     b <- tripletPlus
>     c <- tripletPlus
>     return $ a ++ b ++ c

> shor 0 = do
>     a <- tripletMinus
>     b <- tripletMinus
>     c <- tripletMinus
>     return $ a ++ b ++ c

Repair is a little more complex before. First we use 6 ancilla bits to fix qubits within each triplet using the repetition code three times:

> ancilla (i,j,k) x = ((x!!i) `xor` (x!!j),(x!!j) `xor` (x!!k))

> fix (i,j,k) (e,x) = let (u,v) = ancilla (i,j,k) x
>                         b = case (u,v) of
>                             (0,0) -> x
>                             (0,1) -> flipAt k x
>                             (1,0) -> flipAt i x
>                             (1,1) -> flipAt j x
>                      in (u:v:e,b)

> flipRepair start = start
>             =>= fix (0,1,2)
>             =>= fix (3,4,5)
>             =>= fix (6,7,8)

We also use two ancilla bits to locate phase flips between the triplets. phaseFlipFix will be applied in a rotated basis inside phaseFlipRepair:

> phaseFlipFix :: ([Integer],[Integer]) -> Q ([Integer],[Integer])
> phaseFlipFix (e,y) = do
>     let u = (y!!0) `xor` (y!!1) `xor` (y!!2) `xor` (y!!3) `xor` (y!!4) `xor` (y!!5)
>     let v = (y!!3) `xor` (y!!4) `xor` (y!!5) `xor` (y!!6) `xor` (y!!7) `xor` (y!!8)
>     let c = case (u,v) of
>             (0,0) -> y
>             (0,1) -> y ->- flipAt 6
>             (1,0) -> y ->- flipAt 0
>             (1,1) -> y ->- flipAt 3
>     return (u:v:e,c)

> phaseFlipRepair x = x
>   ->< mapSndM (mapM hnot)
>   ->- collect
>   ->< phaseFlipFix
>   ->- collect
>   ->< mapSndM (mapM unhnot)
>   ->- collect

Now we can pull all of this together into a single process to repair both kinds of error:

> shorRepair x = x ->- flipRepair ->- collect ->- phaseFlipRepair ->- collect

So now all that's needed is to try corrupting some qubits to see if we can repair them.

> state3 = return [] <*> shor 1
> state3' = state3 ->< mapSndM (pFlipAtM 0)
> state3'' = state3' ->- shorRepair

Let's try a more complex superposition with a more complex corruption:

> state4 = return [] <*> (0.6.*shor 0 + 0.8.*shor 1) ->- collect
> state4' = state4 ->< mapSndM (pFlipAtM 5) =>= mapSnd (flipAt 5) ->- collect
> state4'' = state4' ->- shorRepair ->- collect

This is all well and good. But there are more ways to corrupt a qubit beyond simple flips. Let's mess up our state using two additional kinds of 'noise' that we haven't discussed. Firstly, instead of just flipping, we'll rotate one of the bits. Secondly, we'll additionally 'leak' that bit to the environment so as to cause decoherence. We've discussed solutions to neither of these problems, so let's see what happens.

> leak i x = do
>   (e,v) <- x
>   return ((v!!i):e,v)

> state5 = state4
> state5' = state5 ->- leak 4 ->< mapSndM (mapAtM 4 (rotate (pi/10))) ->- collect
> state5'' = state5' ->- shorRepair ->- collect

A minor miracle has taken place. The Shor code is robust against just about anything you can throw at a single qubit. In fact, any code that can survive flips and phase flips of a single qubit can survive any kind of corruption, including entanglement with the environment. (Despite having read the proof, I was still pretty surprised when the code actually worked!) And that gives researchers one hope that these methods might serve as ways to make quantum computing reliable.

I'll mention a couple more things in closing. Firstly, nowhere above have I mentioned wavefunction collapse, observations or superoperators. You don't need any of these things to define the Shor code and show that it can correct errors. (You probably do need some extra assumptions about quantum mechanics if you want to calculate probabilities of various outcomes in the event that the code fails - eg. due to multiple qubits being corrupted.) And secondly, despite it being my April 1st post last year, I still have the doubts I expressed here about quantum computing.

Oh...and expect those lines to take a few seconds to evaluate. When you simulate the Shor code you're performing computations in a 512-dimensional vector space!

Bibliography

1. Peter Shor's original paper has good discussion of the issues but I found the description of the correction algorithm a little laconic.
   
2. John Preskill's Lecture Notes were probably the most useful document.
   
3. A student term paper by Dheera Venkatraman was also pretty useful though I was confused for a while by what might be considered errors in it.
   

Update
Fixed a couple of minor typos. And remember this code should just run, without modifications, in GHCi, so you can play with it yourself.