Data recovery: a toy example

In Ethereum 2.0's data availability model, we have a need to be able to reconstruct chunks of data when up to half of that data has been lost or withheld. I recently implemented some code to do this, but I always like to understand what's really going on, so I put together this worked example.

This example is going to be super unambitious. We are going to start with two items of data. We will encode these into four items of data, and then throw two of them away. Finally we will recover the original two items.

At this scale, you don't need to be very sophisticated. If our data are

$[a,b]$, we can extend it to

$[a,b,a+b,a-b]$: now we can lose any two of these values and still recover the original numbers (try it!). But for much larger scales, we need a general technique that can be efficiently implemented.

Such a technique is described by Vitalik Buterin in a post Reed-Solomon erasure code recovery in n*log^2(n) time with FFTs. My example here essentially follows that description and notation. It will seem massively over-complex when applied to our meager two data, but remember that the method is both general and efficient.

If you would like to play with it yourself, I've coded up the toy example below in this repo.

The basic equipment

The basic items of equipment we shall be using are polynomials, roots of unity, and Fourier transforms. These are all interrelated in very cool ways.

Polynomials

If we have four items of data,

$[d_0,d_1,d_2,d_3]$, we can make a polynomial out of them like this:

$P(x)=d_3{x}^3+d_2{x}^2+d_{1}x+d_0$. The polynomial can then be evaluated for whatever value of

$x$ we like.

Roots of unity

Powers of primitive roots of unity are the key to making this whole thing work efficiently. They allow us to take lots of shortcuts, and are the points at which we will evaluate our polynomials. If a number is an

$n$th root of unity, then raising it to the power of

$n$ results in

$1$. If it is a primitive

$n$th root of unity, then no lower power of

$n$ equals

$1$.

For our example, we will be using four powers of a primitive fourth root of unity. These are

$[1,i,-1,-i]$, where

$i$ is a square-root of minus one, an "imaginary" number. The quantity

$i$ is a primitive fourth root of unity, and the four powers are generated from from it thus:

$[i^0,i^1,i^2,i^3]$. We could have chosen

$-i$ as our generator–-it is also a primitive fourth root–-and ended up with

$[1,-i,-1,i]$ which would work just fine, but we will stick with convention. Collectively, to save words, I'll be referring to this list below as "the roots of unity" (which they are, albeit not all primitive), but bear in mind that they are generated from a single primitive root.

Fourier transforms

We're going to need to evaluate our polynomial at each of these four roots. Evaluating by hand gives the following relations,

$$\begin{array}{rl}{e}_0& =P(i^0)=P(1)=d_0+d_2+d_1+d_3\\ e_1& =P(i^1)=P(i)=d_0-d_2+i(d_1-d_3)\\ e_2& =P(i^2)=P(-1)=d_0+d_2-d_1-d_3\\ e_3& =P(i^3)=P(-i)=d_0-d_2-i(d_1-d_3)\end{array}$$

where I have grouped terms suggestively. The upshot of this is that, if we have the coefficients of a polynomial (the

$d_j$) and want to know the polynomial's values at certain roots of unity (the

$e_j$), we can just apply a transformation to the coefficients, which is a shortcut when compared to evaluating the polynomial every time.

Magically, it turns out that this is identical to taking the Fourier transform of the coefficients, and a very fast way of doing this is known: the fast Fourier transform (FFT).

Taking the Fourier transform of a polynomial's coefficients gives you its evaluations at a set of roots of unity.

Interestingly, we can do some arithmetic to invert the simultaneous equations above, with the result:

$$\begin{array}{rl}{d}_0& =\frac{1}{4}(e_0+e_2+e_1+e_3)\\ d_1& =\frac{1}{4}(e_0-e_2-i(e_1-e_3))\\ d_2& =\frac{1}{4}(e_0+e_2-e_1-e_3)\\ d_3& =\frac{1}{4}(e_0-e_2+i(e_1-e_3))\end{array}$$

These equations give us a way to reconstruct (interpolate) a polynomial, by calculating its coefficients given only its values at the roots of unity. This is identical to an inverse Fourier transform, which can also be done very efficiently (you can see that it's almost the same as a forward Fourier transform).

Taking the inverse Fourier transform of a polynomial's evaluations at a set of roots of unity gives you its coefficients.

Multiplication and division

Other useful tools are polynomial multiplication and division. Polynomial multiplication is taking, say,

$x^2+3x+2$ and multiplying by

$x-1$ to get

$x^3+2x^2-x-2$. For small polynomials, this is easily done via long multiplication; for large polynomials, we can use fast Fourier transforms for a big speed up.

The convolution theorem tells us that, if we have polynomials

$P(x)$ and

$Q(x)$, then the product

$(P\ast Q)(x)$ is given by

$\mathtt{i}\mathtt{n}\mathtt{t}\mathtt{e}\mathtt{r}\mathtt{p}\mathtt{o}\mathtt{l}\mathtt{a}\mathtt{t}\mathtt{e}\mathtt{(}\mathtt{e}\mathtt{v}\mathtt{a}\mathtt{l}\mathtt{u}\mathtt{a}\mathtt{t}\mathtt{e}\mathtt{(}\mathit{P}\mathtt{)}\times \mathtt{e}\mathtt{v}\mathtt{a}\mathtt{l}\mathtt{u}\mathtt{a}\mathtt{t}\mathtt{e}\mathtt{(}\mathit{Q}\mathtt{)}\mathtt{)}$ where

$\times$ is point-wise multiplication. As we saw above, the interpolate and evaluate steps can be done quickly with fast Fourier transforms. If

$P$ and

$Q$ are of order

$n$ in degree, then the long multiplication method takes

$\mathcal{O}(n^2)$ operations, whereas the convolution method takes

$\mathcal{O}(n\log n)$ operations.

Polynomial division can be done in the same way by replacing pointwise multiplication with pointwise division. Once we have the Fourier transforms, it's easier than implementing polynomial long division.

Finite fields

As a final note about the toolkit: if we were doing this "for real", we'd be using a finite field rather than the complex numbers for our coefficients and evaluations. Appropriately constructed finite fields have roots of unity that we can use to construct Fourier transforms (sometimes called number theoretic transforms in that context), and all the maths carries over. The full implementation uses a 256-bit finite field. I've stuck to complex numbers for this example as it seems a little more intuitive to me.

Update: I've added an example using the integers modulo 17 to the repo. In that field, 4 is a primitive fourth root of unity (we could also use 13), and our roots become [1, 4, 16, 13]. Everything else is unchanged.

The example

With our toolkit complete, let's attack our toy example.

These are my starting data, the data I wish to make recoverable:

$[5,7]$

Step 1: encode the data

To encode these data we first extend with zeros from a length of two to a length of four. In the general case, we double the length with zeros: it is using this redundancy to extend our domain that is the essence of being able to later recover the data after losing up to half of it.

The padded data:

$$[5,7,0,0]$$

We treat this padded data as a polynomial,

$D(x)=7x+5$, and encode the data by evaluating the polynomial at the roots of unity. We can do this via our forward FFT.

The encoded data, which is the evaluation of

$D(x)$ at the roots of unity:

$$[12,5+7i,-2,5-7i]$$

By substituting

$d_0=5$,

$d_1=7$, and

$d_2=d_3=0$, you can see how this relates to the explanation of evaluation above:

$$\begin{array}{rl}{e}_0& =D(i^0)=D(1)=5+7\\ e_1& =D(i^1)=D(i)=5+7i\\ e_2& =D(i^2)=D(-1)=5-7\\ e_3& =D(i^3)=D(-i)=5-7i\end{array}$$

Step 2: throw away some data!

For this demo, we shall discard two elements of the encoded data. It could be any two, but we will lose indices 1 and 2, and replace them with zeros.

Encoded data with mising elements:

$$[12,0,0,5-7i]$$

In Vitalik's article This is the evaluation of a polynomial

$E(x)$ that we never actually reconstruct.

Step 3: make the "zero polynomial"

For reasons that will hopefully be clear in a moment, we need to create a polynomial that evaluates to zero at the same places as the missing indices. We'll call this the zero polynomial,

$Z(x)$, and in our example it must evaluate to zero at the first and second powers of the root of unity,

$x=i^1=i$ and

$x=i^2=-1$, corresponding to indices 1 and 2.

Making zero polynomials is conceptually simple: we just multiply up the

$(x-x_j)$ for all the

$x_j$ values that we want to evaluate to zero. In our case:

$Z(x)=(x-i)(x+1)=x^2+(1-i)x-i$. In practice, for large sizes, calculating the zero polynomial ends up being the slowest part of the whole procedure.

$$Z(x)\text{: }[-i,1-i,1,0]$$

Step 4: create the polynomial

$(E\ast Z)(x)$

Recalling polynomial multiplication, if we have the evaluation of

$E$ (we do), and the evaluation of

$Z$ (we already have the two zeros, and we can work out the remainder), then we can interpolate their pointwise multiplications to get the polynomial

$(E\ast Z)(x)$.

Evaluation of

$Z$:

$$[2-2i,0,0,-2-2i]$$

Pointwise multiplication of the evaluation of

$Z$ with the evaluation of

$E$:

$$[24-24i,0,0,-24+4i]$$

Finally, we interpolate this using the inverse Fourier transform to obtain the polynomial

$(E\ast Z)(x)=7x^3+(12-7i)x^2+(5-12i)x-5i$. This is represented as,

$$(E\ast Z)(x)\text{: }[-5i,5-12i,12-7i,7]$$

Step 5: ponder the magic

Now the magic happens. We are trying to reconstruct our original polynomial,

$D(x)$, so that we can get our data back. Consider

$(E\ast Z)(x)$ and

$(D\ast Z)(x)$. By construction, both of these polymonials evaluate to zero at the missing indices, 1 (

$x=i$) and 2 (

$x=-1$). And both evaluate to the same values at the non-missing indices, 0 (

$x=1$) and 3 (

$x=-i$). Third degree polynomials are fully determined by four values, which means that these two are the same polynomial. Take a moment to convince yourself of this.

As a result of all this, if we can divide the polynomial

$(E\ast Z)(x)=(D\ast Z)(x)$ by the polyomial

$Z(x)$, then we recover

$D(x)$ and our work is complete.

Step 6: divide

$(D\ast Z)(x)$ by

$Z(x)$

If only it were that simple! We still have a little hurdle to jump. We would normally do the polynomial division by using the values of the polynomials at the roots of unity. Unfortunately, we know that for indices 1 and 2, both polynomials evaluate to zero because of how we made them, and we have no way to assign values to the

$0/0$ result.

Step 6a: scaling the polynomials

Thankfully, we can use a little trick to avoid making evaluations at these zero points. The trick is to scale the polynomial by transforming

$x\mapsto x/k$ for some constant

$k$. (This is called shifting the polynomial in the original article, but it isn't really a shift.)

This is reliable because we know that

$Z(x)$ is only zero at two points. As long as our scaled version is not evaluated at those points then we'll not be trying to divide by zero. This works out fine as long as the scale factor,

$k$, is not

$0$,

$1$, or any of our roots.

The shift is easy to implement: we just multiply the polynomial's coefficients by the appropriate power of

$k$. Shifting by

$k=2$ gives us the following:

$$\begin{array}{rl}(D\ast Z)(x/2)\text{: }& [-5i,10-24i,48-28i,56]\\ Z(x/2)\text{: }& [-i,2-2i,4,0]\end{array}$$

Step 6b: division via convolution

To do the division, we first form the evaluations of the two polynomials at the roots of unity, using the forward Fourier transform method we are already familiar with.

Evaluation of scaled

$D\ast Z$:

$$[114-57i,-24-23i,-18-9i,-72+69i]$$

Evaluation of scaled

$Z$:

$$[6-3i,-2+i,2+i,-6-3i]$$

Now we need to divide the first by the second, entry by entry. Division of complex numbers is a little non-obvious, but well-defined:

$$[19,5+14i,-9,5-14i]$$

As the last step in the division, we interpolate this back to a polynomial via the inverse Fourier transform:

$$D(x/2)\text{: }[5,14,0,0]$$

Step 6c: unscaling

We're almost there! We just need to unscale

$D(x/2)$ by dividing each coefficient by the appropriate power of

$k=2$.

$$D(x)\text{: }[5,7,0,0]$$

Step 7: profit!

We have recovered the polynomial

$D(x)$. And voila, our original data is intact there, in the first two places:

$[5,7]$.

Conclusion

As noted at the top, all this palaver is a bit ridiculous for our tiny example with only two data items. However, the real value is that the method described above can scale efficiently to large sizes, using essentially the same toolkit. My full code takes about a second to recover 2MB of data when fully half of the encoded data has been lost. About half of that time is spent in generating the zero polynomial - possibly the simplest part conceptually, but apparently hard to do quickly.

As mentioned above, demo code implementing the above toy example is available on my GitHub.

Related work

• This toy example coded up.
• Vitalik: Reed-Solomon erasure code recovery in n*log^2(n) time with FFTs
• Ethereum Data Availability Sampling Phase 1 Proposal
• Dankrad Feist: polynomial_reconstruction.py. An implementation in Python.
• Protolambda: go-kzg. Includes an implementation in Go, see recover_from_samples.go.
• Me: c-kzg. Largely follows go-kzg. See recover.c. Work in progress.
• Vitalik's artice on fast Fourier transforms has lots of helpful insights.