Derandomization with Pseudorandomness
Emin Karayel

Legal Aspects for Software Developers Interested in Generative AI Applications
    Joel Mittel
    April 29, 2024
==============================================================================


Derandomization techniques are often used within advanced randomized algorithms.
In particular, pseudorandom objects, such as hash families and expander graphs, are key components of such algorithms, but their verification presents a challenge.
This work shows how such algorithms can be expressed and verified in Isabelle and presents a pseudorandom objects library that abstracts away the involved deep algebraic/analytic results.
Moreover, it presents examples that show how the library eases and enables the verification of advanced randomized algorithms.
Highlighting the value of this framework is that it was recently used to verify the optimal-space distinct elements algorithm by Błasiok from 2018,
which relies on the combination of many derandomization techniques to achieve its optimality.



§ INTRODUCTION

Computer scientists are familiar with pseudorandom number generators (PRNGs), functions that consume a small sequence of seed bits and output a much larger sequence of pseudorandom bits.
Well-known examples are the Lehmer random number generator <cit.>, Blum Blum Shub <cit.>, and the Mersenne-Twister <cit.>.
The randomness of these functions can only be verified empirically <cit.>.
But there are no mathematical, in particular, probability-theoretic theorems about them that would make them useful in the design of provably correct randomized algorithms.[
There are conditional results, , with the assumption of the difficulty of factoring.
Moreover, there are rigorous results about the period of some PRNGs, but they are usually not helpful with respect to the correctness of algorithms using such PRNGs.]

On the other hand, there are pseudorandom objects[The term was coined by Vadhan <cit.> but it is applied here more generally than he does.], such as k-independent hash families and expander graphs, fulfilling such theorems, for example, limited independence or Chernoff-type tail bounds.
They are often critical components of advanced randomized algorithms.
Similar to PRNGs, they are functions taking a small sequence of seed bits and returning a much larger sequence of output bits. Often, the result is represented as a function or a vector.

Pseudorandom objects are used as a derandomization technique and thus improve the randomness usage of simpler algorithms.
This is especially interesting in cases where the saving in coin-flips translates into the reduction of space usage, because the algorithm needs to preserve the coin-flips during its operation.
But, as we will see in this work, the contruction of pseudorandom objects relies on advanced mathematical foundations, such as (higher-order) finite fields and spectral graph theory.
Making their formalization in proof assistants a challenge, as well as, the verification of randomized algorithms relying on them.
This work shows how advanced randomized algorithms can be expressed and verified in Isabelle <cit.> (and in general on every proof assistant), through the use of a composable pseudorandom objects library.



  
A motivating example: To illustrate the utility of pseudorandom objects, let us look into the following simple randomized streaming algorithm. For readers less familiar with Isabelle syntax: It is possible to skip these. There is always a precise description of these algorithms in the surrounding text.

fun example0pt1  nat  nat list  real pmf  where example0pt1 n xs      do 0pt      h  prod0ptpmf 0n  pmf0ptof0ptset 11      return0ptpmf x  xs h x2    
This is a randomized approximation algorithm for the second-frequency moment[The second frequency moment is the square sum of the occurrence counts of each value in the stream.
It is useful in assessing the skew of the stream, which can, for example, be needed in monitoring applications or query optimization <cit.>.
The algorithm <cit.> is very concise, which makes it useful as an illustrative example.] F_2 of a sequence x_1,…,x_m where the sequence values are assumed to be in the set U := {0,…,n-1}.

In the first step, the algorithm randomly chooses a function h from U (the set containin the sequence values) into the two-element set {-1,1}, , it assigns each element of U a random signed unit.
The algorithm then returns the square of the sum of the signed units assigned to each sequence element, : (∑_i=1^m h(x_i))^2.
It is effortless to verify that the algorithm returns an approximation of F_2. Indeed, the expected result is:

    (∑_i=1^m h(x_i))^2 = (∑_u∈ U c(u) h(u))^2 = ∑_u ∈ U∑_v ∈ U c(u) c(v)  h(u) h(v) = ∑_u ∈ U c(u)^2 =: F_2 .

Here, c(u) denotes the occurrence count of u in the stream x_1,…,x_m.
The key step in the above derivation is the third equality, where the expectation h(u) h(v) is 0 if u ≠ v and 1 otherwise.
The last equality is the definition of the second frequency moment.
A similar argument can be used to derive:

    (∑_i=1^m h(x_i))^2 ≤ 2 F_2^2 .

The algorithm above is, however, not very efficient: It requires (n + ln m) bits of memory.
At least n bits are required to store the function h and another ln m bits for the accumulation register, while evaluating the sum.
As promised, it can be improved using the “magic of pseudorandomness”.
Instead of choosing the function h completely randomly, which necessarily requires (n) bits of storage space, it can be selected using a pseudorandom object with a small seed space.
To achieve that the function has to be stored implicitly using the seed, instead of, as a dictionary.
Indeed, a careful examination of Eq. <ref>, shows that the values of the function h only have to be pair-wise independent.
(A similar observation can be made for Eq. <ref>, where only 4-wise independence is sufficient.)

This means the pseudorandom object, we choose is required to form a 4-wise independent family of functions.
The following example is the formalization of the updated algorithm:

fun example0pt2  nat  nat list  real pmf  where example0pt2 n xs      do 0pt      h  sample0ptpro  4 n  11      return0ptpmf x  xs h x2    
The new algorithm has the same properties as example1 (, fulfilling Eq. <ref> and <ref>) but requires only (ln n + ln m) bits of memory.
The only difference is in Line 4, where h is now obtained using the pseudorandom object specified by the term  4 n  11, which requires (ln n) bits of input as seed and returns a function from {0,…,n-1} to {-1,1} (which is stored implicitly using the (ln n) bits of space, instead of as a dictionary).
If the seed is chosen uniformly at random, then the evaluation of the function at any point in {0,…,n-1} forms a random variable.
Crucially, the family of all n random variables is 4-wise independent, and each is distributed uniformly on {-1,1}.
The operation samplepro interprets a pseudorandom object as a probability space.
It uniformly chooses a random seed from the seed space of the pseudorandom object and evaluates the object at that seed. 

The expression  4 n  11 from above represents a combined pseudorandom object. 
It is defined using two constructs: The expression  11 describes a pseudorandom object with a single-bit input that outputs -1 or 1 depending on the bit.
The construct  k n P is a k-independent hash family, with domain {0,…,n-1}, and the range P, where P is another pseudorandom object.
This means it is easy to construct application-specific pseudorandom objects from a small set of primitives.
The presented ones are part of the libraries <cit.> (available in the AFP <cit.> for Isabelle.)

The step from example1 to example2 demonstrates derandomization with pseudorandom objects.
The above is just a teaser example and will be further expanded on in Section <ref>, where combinations of multiple types of pseudorandom objects, in particular, including expander graphs, will take the stage.



  
Contributions: The main novel results presented in this work are:


  * A novel approach to constructing pseudorandom objects based on simple primitives and the introduction of notation for them that allows a succinct representation within algorithms using them.

  * Formal verification of an executable library for computations in non-primitive finite fields, , GF(p^n) for p prime and any n ≥ 1.
(Prior work in proof assistants <cit.> focused on the much simpler case of n=1.)
This is the basis for the first pseudorandom object — k-independent hash families — detailed out in Section <ref>.

  * Formal verification of expander graphs and a construction of a strongly explicit family for every degree and spectral bound and the verification of tail bounds for random walks in expander graphs.
This is the basis for the second pseudorandom object — random walks in expander graphs — detailed out in Section <ref>.

It should be noted that the examples are verified in Isabelle <cit.>.
The supporting libraries require more than 16400 lines of code[Appendix <ref> provides a comprehensive of the formalizations mentioned here.] because the construction of the pseudorandom objects builds on mathematical foundations such as higher-order finite fields and spectral graph theory.
However, the libraries can be used easily based on high-level results that abstract away implementation details.

The following sections will discuss related work and additional important background.
The formalization of pseudorandom objects in an abstract way is detailed out in Section <ref>, followed by the concrete pseudorandom objects k-independent hash families (in Section <ref>) and expander graphs (in Section <ref>).
Section <ref> presents an application of the framework, where multiple types of pseudorandom objects are going to combined to improve the previous introductory example.
Section <ref> concludes with a summary and future research opportunities.



§ RELATED WORK



 §.§ Derandomization with Pseudorandomness

As the abstract was alluding to, an exciting use case for this framework is the distinct elements problem — estimating the number of distinct elements in a stream in a single pass with limited mutable state.
The research on the problem started in 1985 by Flajolet <cit.>, followed by successive improvements <cit.> until an optimal solution was found by Błasiok in 2018 <cit.>.
In 2023 <cit.>, I wrote about the formalization of it using the framework for pseudorandom objects presented here.
It is an improved version of Błasiok’s solution.
The concrete pseudorandom object used there is an expander random walk chained with a second expander random walk chained with a product of pseudorandom objects for k-independent hash families.
Indeed, the formalization of this algorithm was a big part of the motivation for the design of the library presented in this work.
At the current state of research, pseudorandom objects are essential for the distinct elements problem.
A notable exception is the solution by Chakraborty  <cit.> that does not require pseudorandom objects, but its space complexity is asymptotically worse than the optimal solution.

There are, of course, many further randomized algorithms, that rely on pseudorandom objects, such as heavy hitters <cit.>, frequency moments <cit.> and other stream statistics <cit.>, graph matching <cit.>, graph sparsification <cit.> and low-distortion-embeddings <cit.>.

Another motivating influence for this work is Goldreich’s chapter <cit.> on samplers: Algorithms that probe a function at select points and obtain an estimate for the average of the function.
The performance of samplers is measured with respect to the number of probing points, the number of random bits consumed, and the quality of the estimate.
Goldreich also relies mainly on the two fundamental pseudorandom objects mentioned here: expander graphs and hash families.
In contrast to this work, he does not treat pseudorandom objects as first-class objects but the samplers, which use pseudorandom objects as a component.
However, that viewpoint prevents some of the applications I describe here, including the example from the introduction.
(Conversely, samplers are easily represented with the framework presented here.)

An ancestor of the notation developed here for pseudorandom objects was the work of Kane  <cit.>, who introduced a notation for hash families and forming products.
However, they do not consider chaining pseudorandom objects for which there is no prior art.

Vadhan <cit.> has written a monograph and survey on pseudorandomness.
As far as I can tell, he coined the term pseudorandom object, and his work stands out by covering a large number of derandomization techniques and highlighting connections and differences.
Important to note is that in this work, the term pseudorandom object is used more broadly than in Vadhan’s monograph.



 §.§ Randomized Algorithms in Proof Assistants

The best way for the verification of randomized algorithms in Isabelle/HOL is based on the Giry monad <cit.>, which is the method applied in this work.
Eberl  <cit.> have done foundational research and verified several randomized algorithms demonstrating the feasibility of the approach.
It should however be noted that this is still an active research area <cit.>.
Similarly there is work towards transfering the approach to other proof assistants <cit.>.

There is however not yet much research on derandomization techniques (even k-independent hash families).
The only research I know of is the formalization of an approximate model counting algorithm that relies on the construction of a 3-universal XOR-based hash-family by Tan  <cit.> in Isabelle/HOL.
Their work uses some libraries I developed during this work.
Conversely, Yong Kiam Tan made several, much appreciated, improvements to my AFP entry on the median method <cit.>.

In 2021 <cit.>, I presented the formalization of randomized algorithms for frequency moments.
It only uses hash families (and only based on finite fields with simple prime order) and does not use other pseudorandom objects or combinations of these.
Moreover, the focus of this work is on derandomization techniques, as opposed to the verification of concrete algorithms.



 §.§ Other Derandomization Methods

Besides the use of pseudorandom objects, there are also other derandomization methods. The two most prominent are:


  * Non-uniform (non-constructive) derandomization

  * Method of conditional expectations (or probabilities)

The first — non-uniform derandomization — is based on Adleman's theorem <cit.>, who showed that any randomized algorithm with a finite set S of possible inputs can be derandomized entirely.
This works by first amplifying it until the failure probability is smaller than 1/|S|.
For such an algorithm it is possible to show that there must exist a fixed choice of coin flips, for which the algorithm must succeed for every input. This follows from a union-bound.
The result is sometimes abbreviated as: RP⊆P/poly.
The theorem does not provide a way to find such a sequence of coin-flips.
The topic is less relevant for formal verification: Constructing such an algorithm would require a brute-force search, and it is trivial to verify the correctness of an algorithm with a finite input set.

The second — method of conditional expectations <cit.> — can be applied in some specific cases.
Given a randomized algorithm, it is sometimes possible to derive a deterministic algorithm, that performs at least as well as the randomized algorithm on average.
This requires that the probability space can be decomposed recursively into a search tree, where conditional expectations can be computed for each node.
I have formalized two algorithms based on this technique. It turned out that it is straightforward to verify such algorithms.
This contrasts pseudorandom objects, which require large foundational libraries to work with them.



§ BACKGROUND ON RANDOMIZED ALGORITHMS

This section provides a very small introduction, into the representation of (randomized) algorithms in Isabelle.

Isabelle is a meta-logic for which several instantiations exist <cit.>.
The most prominent, and the framework used in this work, is HOL. Other instantiations are, for example, First-order logic (FOL) or Zermelo-Fraenkel Set Theory (ZF).
It can be described as a typed version of classical mathematical logic.[In my opinion HOL is actually closer to informal mathematics used today than, , Isabelle/ZF. For example a textbook proof would never consider the equality of a group and a neighborhood basis of a topology, they are, informally, considered different types of objects, even though formally everything is a set.]

The type system is similar to that of many functional programming languages, such as Haskell or ML.
The primitives nat, real represent natural (resp. real) numbers.
It is also possible to construct more complex types from simple ones, for example nat  real represents a function from natural to real numbers.
A bit unexpected is that type constructors are written post-fix, , nat list represents a list of natural numbers.
Type variables are denoted using primed letters, : a, b.
The double colon is used signify membership of an expression or name in a particular type, : n  nat.

Isabelle's expressions resemble standard mathematical notation closely, but there are sometimes slight deviations. The following table summarizes such cases:




Description     Isabelle Term     Classical Math 

Cardinality of a set S     card S     |S| 

Summation over a list xs  x1xm     x  xs f x     ∑_i=1^m f(x_i) 

Summation over indices     i  n f i     ∑_i=0^n-1 f(i)

Interval with strict bounds     ab     {x | a < x < b}

Interval with inclusive bounds     ab     {x | a ≤ x ≤ b} 

Strictly upper bounded set     b     {x | x < b}





The interval notation is flexible. The above are just examples, but any combination of unbounded, strict or inclusive lower and upper bounds are possible.

An important type constructor for this work is the one describing probability mass functions (PMFs): pmf.
They are used to represent finite/discrete probability spaces and form a special class of general probability spaces, with the condition that the σ-algebra is the powerset of the universe of the type. (Eberl <cit.> provides a more thorough introduction.)
This greatly reduces the difficulty of working with them because measurability becomes trivial.
Because distributions of randomized algorithms will always be supported on a countable set, the condition on the σ-algebra is not a restriction.
For example real pmf represents a discrete probability space, whose elements are real numbers.
The following table provides a brief overview of the terms that can be used to construct PMFs.




Term     Description 

pmfofset S     For a finite set S, the uniform probability space on S. (Every element of S is equiprobable.) 

pmfofmultiset S     For a multiset S, the probability space on the elements of S, with the probability of elements being proportional to their multiplicity. 

pairpmf A B     Given probability spaces A and B, the probability space supported on the product of the supports of A and B. (With each pair (x,y) having the probability of the product of the probabilities of x and y.) 

prodpmf I M     Finite product of multiple probability spaces. The elements are maps with a domain I, the index set, and values in M(i) (depending on the index). 

mappmf f A     The probability space representing the distribution of the random variable f over the probability space A. 

returnpmf x     The probability space of the singleton {x}. 

bindpmf p f     This represents algorithmically the composition of two randomized algorithms p and f, where f is executed using the result of p. Mathematically, bindpmf p f is the probability space for which an event E has the probability _bindpmf p f E := ∫_f x (E)   d p(x). 

Note also that: mappmf f p  bindpmf p returnpmf  f 





The functions bindpmf and returnpmf form the Giry-monad for probability mass functions,
because of that, randomized algorithms can be represented more intuitively using the do-notation.
As an example  bindpmf p x returnpmf q x can be written as:

  do     x  p    returnpmf (q x)  
which represents the distribution of a randomized algorithm that obtains x from the probablity space p and returns q x.



§ PSEUDORANDOM OBJECTS 

This section is less exciting, but it is nevertheless necessary to tackle the basic definitions before discussing the advanced objects.
Let us take a brief look at the representation of pseudorandom objects in Isabelle <cit.>.
This is done using a record:

record a pseudorandomobject   prolast  nat  proselect  nat  a

Note that the record type pseudorandomobject is parameterized over the type of the elements of the pseudorandom object a.
A record in Isabelle can be thought of as a tuple type, but with the advantage that the components can be accessed using named selectors. In this case, these are: prolast, proselect.

The selector prolast represents the object's size minus 1.
The selector proselect is a function that returns an element from the object, given an index between 0 and prolast (inclusive).
Note that this means that pseudorandom objects always have at least one element.
Effectively, a pseudorandom object is just a non-empty multiset from which it is easy to sample.
The fact that empty pseudorandom objects are disallowed comes in handy because this means that we can associate with any pseudorandom object: a probability space.
The function samplepro does exactly that.
It is defined as:

samplepro P = mappmf proselect P pmfofset 0prolast P
While we will keep explicitly using the term samplepro, it is actually defined as a coercion. Whenever a pseudorandom object is used in a place that requires a probability space, Isabelle would automatically convert it to a probability space by introducing an application of samplepro. (Similar to the automatic coercion of integers to real numbers.)

Note that the function proselect is not required to be injective.
(Such a condition would incur a performance cost for the more complex pseudorandom objects, such as expander graphs.).
For an element x of a pseudorandom object, we call the count of indices in 0prosize P mapped to x the multiplicity of x.

The following are trivial pseudorandom objects and combinators defined in the library <cit.>:


  * If xs is a non-empty list, then  xs is a pseudorandom object. The multiplicity of an element is its occurrence count in the list. In particular: samplepro  xs  pmfofmultiset mset xs.

  * If n > 0, then  n is the pseudorandom object for the natural numbers {0,…,n-1} (each with multiplicity 1). In particular: samplepro  n  pmfofset 0n.

  * If P and Q are pseudorandom objects, then P P Q denotes the product pseudorandom object. Its size is the product of the sizes of P and Q and consists of pairs whose first component is from P and whose second component is from Q. The multiplicity of a pair (x,y) is the product of the multiplicities of x in P and y in Q. 
 In particular: 
samplepro P P Q  pairpmf samplepro P samplepro P.

In the following two sections, we will dive into the more exciting pseudorandom objects. The following section is dedicated to k-independent hash familes and Section <ref> to random walks in expander graphs.



§ HASH FAMILIES 

The textbook <cit.> definition of k-independent hash families is as a family of functions from a finite domain D to a finite range R.
For any domain element d ∈ D, we can associate a random variable X_d — the evaluation of the function at the point d, where the probability space is the hash family.
The key condition k-independent hash families fulfill is that the random variables X_d are k-wise independent and uniformly distributed in the range R.

The classic construction for such families is due to Wegman and Carter <cit.>, where the hash functions are just polynomials (with degrees strictly smaller than k) over finite fields. It is straightforward to verify the k-independence, which follows from the algebraic fact that there is exactly one polynomial (of degree smaller than k) interpolating a given set of k points. (See also Figure <ref>.)



In particular, using the Galois field GF(q), where q is a prime power, we obtain a family in {0,…,q-1}→{0,…,q-1}, requiring a seed argument from {0,…,q^k-1} (the coefficients of a degree k polynomial).

It is straightforward to construct k-independent hash families with different domains and ranges; for example, if q' divides q, then composition with x → x  q' does not affect the k-wise independence or uniform distribution property.
Similarly, both properties are preserved if the domain is restricted arbitrarily.
This means we can construct a k-independent family for functions between any domain of size n and range of size m, where n is arbitrary, and m must be a prime power.[This works by choosing q as p^k where k = max(log_p n, log_p m) and p is the unique prime dividing m.]

The approach taken in the formalization <cit.> goes beyond that because hash families may be used to select the seed of another pseudorandom object.
This is key to the compositionality of the pseudorandom objects library.
Instead of selecting a finite set for the range, in the formalization, the range space can be any other pseudorandom object (if its seed space's size is a prime power).
(Internally, this just represents composition with the selection function of the new pseudorandom object.)
Notationally, this is represented by the expression  k n P where P must be a pseudorandom object with a seed space whose size is a prime power.
The above is not a restriction because we can, of course, just write:  k n  m to obtain a family from {0,…,n-1} to {0,…,m-1} (if m is a prime power.)

Note that the seed space of the Carter-Wegman hash family is always a prime power. This means it is possible to chain hash families.  

The defining property of k-independent hash families is expressed in the following way:

  mappmf restrict J samplepro  k n P  prodpmf J  samplepro P
under the condition that J  0n and card J  k. , if the domains are restricted to any k element set, the functions from the pseudorandom object behave as if they are random functions. More precisely, they behave as if they are sampled independetly using the inner pseudorandom object P.

The fact that the right-hand side is expressed with respect to the probability space of the chained pseudorandom object samplepro P permits the derivation of results in the case of more complex combined pseudorandom objects. For interested readers, Section <ref> contains a step-by-step proof using combinations of pseudorandom objects.  



 §.§ Construction of Finite Fields

Critical to the construction of the Carter-Wegman-Hash families is the construction of finite fields.
For any prime p, the factor ring GF(p) := ℤ/pℤ is a finite field, but not all finite fields can be constructed in this way.
The finite fields of order p^k, where k > 1, can be constructed using factors of polynomial rings.
In particular, if f is an irreducible polynomial of degree k with coefficients in a finite field F, then the factor ring F[X]/ f F[X] becomes a finite field of order |F|^k.
This means the construction of the finite field GF(p^k) reduces to finding an irreducible polynomial of degree k with coefficients in GF(p).

The latter is helped by the fact that the probability of a random polynomial of degree k ≥ 1 being irreducible is at least 1/2k — a consequence of Gauss’s formula counting the number of irreducible polynomials of a given degree over a finite field GF(p):

    1/k∑_d | kμ(d) p^k/d

This follows from a deep result by Gauss, showing that in a finite field F an irreducible polynomial divides X^|F|^n-X if its degree divides n, as well as the fact that X^|F|^n-X is square-free.
Another consequence of the above is Rabin’s irreducibility test <cit.>; a polynomial f is irreducible if it divides X^|F|^deg(f)-X, but is coprime w.r.t. X^|F|^deg(f)/k-X for any prime k dividing deg(f).
The latter suggests a practical algorithm to obtain irreducible polynomials by sampling polynomials randomly and checking for irreducibility. And indeed, this is a practical strategy. (The algorithm's expected iteration count is (k)). By the way, as far as I know, there are no deterministic algorithms that have a similar performance <cit.>.

The AFP entry <cit.> contains a formalization of these results, in particular, a complete characterization of the finite fields, , that every finite field is of prime power order and finite fields with the same order are isomorphic <cit.>, Gauss formula, the Frobenius homomorphism and Rabin’s irreducibility test.
It also introduces efficient algorithms for field operations, , computing the inverse using the extended Euclidean algorithm, as well as long division and fast binary modular exponentiation. 
For interested readers, the following code fragment is the sampling algorithm for irreducible polynomials <cit.>:

partialfunction smpf sampleirreduciblepoly  nat  nat  nat list  nat  where sampleirreduciblepoly p n     do       k  spmfofset pk      let poly  enummonicpoly modring p n k      if rabintest modring p poly        then returnspmf poly1        else tickspmf sampleirreduciblepoly p n    
Note that the algorithm does not terminate unconditionally, but only almost surely, , the probability that it does not terminate is 0.
This is why the sampling algorithm is represented as a partial function in the SPMF monad. Such functions have a different semantics.
In particular ordinary induction proofs can not be applied but must be replaced with a probabilistic convergence proof.
This relies on the framework by Lochbihler <cit.>. The entry also contains several illustrative examples.

The algorithm example2 from the introduction used a deterministic construction of the finite fields.
(Internally, it chooses the lexicographically smallest irreducible polynomial of a given degree.)
However, as we just saw, a randomized selection of the irreducible polynomial is more efficient.
To use that version, the example needs to be modified slightly:

50pt

fun example0pt3  nat  nat list  real pmf  where example0pt3 n xs      do 0pt      h  sample0ptpro  P 4 n  11      return0ptpmf x  xs h x2    
, we replace function applications with monadic applications  and use the constructor P instead of the deterministic version .
Note that  is just the flipped version of the more classical bind operator .

This actually increases the number of coin flips, but the expected number of coin flips remains the same asymptotically.



§ EXPANDER GRAPHS 

The next pseudorandom object is based on random walks in expander graphs <cit.>.
They are complementary to hash families.
While k-independent hash families approximate independent sampling, in the sense that any k element set forms an independent product, random walks in expander graphs form a more analytic approximation of independence. In the sense that an l-step random walk in an expander graph fulfills Chernoff-bounds similar to that of an l-wise indepedendent product. What this exactly means will be explored at the end of this section.

In the example in Section <ref>, we will see how the combination of these two types of pseudorandom objects can lead to impressive results.
Let us first briefly review the concept of expander graphs before presenting the notion of strongly explicit expander graphs, concluding with the construction of the pseudorandom object as random walks in them.

A family of expander graphs is an infinite sequence G_1, G_2,… of undirected d-regular graphs[A graph is d-regular if all of its nodes have degree d.] fulfilling a common expansion property (detailed below).
The size of the graphs increases monotonically but not too fast.
Hoory  require |G_i| < |G_i+1| ≤ |G_i|^2 <cit.>, which we follow here, but note that the growth condition is by no means standard.
Moreover, it is essential to note that parallel edges and self-loops are allowed. (A possible mathematical model for such graphs is as a pair of finite sets (V, E) where E is a subset of V × V × L, here L can be an arbitrary space of labels for the edges.
An alternative would be to use the more classical E ⊆ V × V  but allow E to be a multiset.)



For each graph, we can associate an adjacency matrix.
Given an ordering of the vertices, each cell of the adjacency matrix A = (a_ij) denotes the number of edges between the i-th vertex and the j-th vertex.
The fact that the graphs must be undirected and regular, , the degree of all nodes being the same, means that the matrix is symmetric, and each row and column has the same sum corresponding to the degree of the nodes.
The normalized matrix, , scaling the adjacency matrix by the inverse degree of the graph, is called the stochastic matrix.
It corresponds to the transition matrix of a random walk in the graph.
It can be shown that the eigenvalues of both the stochastic and adjacency matrices are all on the real line in the above conditions.
Additionally, the eigenvalues of the stochastic matrix will be in the interval [-1,1], and 1 is always an eigenvalue, which corresponds to the stationary distribution.
The corresponding eigenvector represents the stationary distribution of the Markov chain.
It is a vector with all components being identical non-zero.[This follows from the regularity of the graph.]

Regarding the expansion condition, different authors use different definitions; commonly, one of the following expansion conditions is employed:


  *  The second largest absolute eigenvalue of the stochastic matrices of all the graphs in the family is bounded by a common value c < 1, ,  λ_a(G_i) < c for all i.

  *  The second largest eigenvalues of the stochastic matrices are bounded by a value c < 1, , λ_2(G_i) < c for all i.

  *  There is a common strictly positive lower bound on the edge expansions h(G_i), where h(G) is the minimum ratio between the number of edges crossing the cut and the (smaller) size of the cut, , 
h(G) = min_0<|S|<|G|E(S, G ∖ S)/min(|S|, |G∖ S|).  

See also Figure <ref> for an example graph with its expansion properties.
It is possible to infer between these three different definitions.
For example, a family fulfilling property <ref> also fulfills property <ref>.
Using the Cheeger inequality, it is possible to derive a connection between spectral expansion and edge expansion:

    d/2 (1 -  λ_2(G)) ≤ h(G) ≤ d √(2(1 - λ_2(G))).

There are the following alternative variational definition for λ_a(G),  λ_2(G), which are useful:


    λ_a(G)     :=    max_x=1, x ⊥ u |x^T A x| = max_x=1, x ⊥ uA x
    λ_2(G)     :=    max_x=1, x ⊥ u x^T A x

where u is the eigenvector (u_i := |G|^-1/2) corresponding to the stationary distribution.

For interested readers, it is straightforward to conclude that random walks in these graphs are rapidly mixing Markov chains <cit.>.
Any of the three expansion conditions automatically imply that the graph must be connected.

The reason we are interested in expander graphs is that random walks in them approximate independent sampling, fulfilling Chernoff bounds similar to ones for independent random variables.

For example, let us assume a μ-fraction of the vertices of an expander graph with n vertices are marked by 1, and the others are marked by 0, , f : G →{0,1} where μ =|{x∈ G. h(x) = 1}|/n. 
If we choose l vertices independently, then the number of vertices hitting the μ-fraction will be concentrated around μ n, with a sharp tail bound fulfilling:


    _v ∈ G^l( | 1/l∑_i=1^l f(v_i) - μ| ≥ c ) ≤ 2 exp(- l c^2 )  for  c ≥ 0


Note that this follows from the classic Hoeffding inequality.
Instead of choosing the vertices independently, we can start at a random vertex and perform a random walk, where each edge of the walk is chosen uniformly from those that are incident to the current vertex.
As alluded to above, this will fulfill a similar Chernoff bound:


    _v ∈Walk(G,l)( | 1/l∑_i=1^l f(v_i) - μ| ≥ (c +  λ_a(G)) ) ≤ 2 exp(- l c^2 )	 for  c ≥ 0


The key advantage of random walks in expander graphs, compared to choosing vertices randomly, is that the number of coin flips necessary to select such a walk is only ln |G| + l ln d — as opposed to l ln |G|, , we can reduce the number of coin flips, paying with the additional additive term λ_a(G) in Eq. <ref>. (See also Figure <ref>.)


 §.§ Formalization



Because of the many possible definitions, the formalization <cit.> is as flexible as possible.
It introduces regular graphs as a specialization of general graphs, introduced in the AFP-library for Isabelle by Noschinski <cit.>.
For such graphs all three types of expansion conditions are introduced.
As well as theorems such as the Cheeger inequality, Eq. <ref> or <ref>, and the previously mentioned Chernoff bound for random walks (Eq. <ref>.)

The proofs rely heavily on matrix algebra.
Note that Isabelle has two major representations of matrices and vectors with different approaches.

Thiemann  <cit.> developed a tuple-based representation, where vectors and matrices are functions paired with their dimension.
This approach is especially useful, when performing induction over the dimension of the vector or matrix, , verifying an algorithm that computes the Jordan normal form.
The main drawback is that they require explicit proofs of dimensional consistency.
For example to apply associativiy A (B C) = (A B) C, we need to check that the row count of A equals the column count of B and that the row count of B matches the column count of C.
This has a detrimental influnce on proof automation.
 
The alternative representation — the one I chose — is part of the HOL-Library and relies on numeral types.
Vectors are represented as functions from a finite type with a size corresponding to the dimension of the vector space to the real/complex numbers.
This representation has the advantage of enabling type-level verification of dimensional consistency, which works seamlessly.
Indeed, it often helps avoid type annotations as they can be inferred automatically.
The approach has been developed by Harrison <cit.>.

To be able to use the type-based representation, it was necessary to rely on the types-to-set mechanism <cit.>.
This works by introducing an assumption that a type exists, whose size matches the vertex set of the graph.
With the assumption it is possible to define the stochastic matrix, and perform algebraic arguments.
The assumption for the existence of the type can then be discharged for theorems that do not refer to it.
This means that the non-algebraic results about expander graphs are available even without such an assumption.
A trick I employed is to introduce λ_a(G), λ_2(G) using the variational definitions (Eq. <ref> and <ref>), , without referring to the stochastic matrix, and show they are equivalent to the spectral definition, , with respect to the multiset of eigenvalues. This allows the expression of all the relevant theorems, without explicit reference to the stochastic matrix, and hence the type pre-condition could be discharged for them.



 §.§ Diagonalization

We will not go into all of the proofs. Nevertheless, the following is an interesting proof technique <cit.> that highlights some of the fine issues when formalizing textbook proofs. Let us recall Eq. <ref>:

    λ_a(G) = max_x=1, x ⊥ u |x^T A x|

where A is the stochastic matrix of a regular graph G, and u is the vector whose components are all |G|^-1/2.

The easiest way to see that this is true is to consider the diagonalization of A, , there exists a unitary matrix U s.t. A = U^T D U where D is a diagonal matrix consisting of the eigenvalues (with their respective multiplicities) of A.
The rows of U consist of eigenvectors of A.
More precisely, the i-th row is the eigenvector, corresponding to the i-th eigenvalue in the diagonal matrix D.
This means the right-hand side of Eq. <ref> equals max_y=1, y ⊥ Uu |y^T D y|.
One of the diagonal elements of D (, d_ii) must be a 1 for which the corresponding row in U should be the stationary eigenvector u.
(Alternatively, Uu will be the i-th basis vector.) From this, we can derive that the right-hand side is indeed λ_a(G), the second largest eigenvalue of A.

Unfortunately, that argument fails, in general, if the multiplicity of the eigenvalue 1 is larger than 1.
(Something that textbook proofs usually ignore.)
However, there always exists a decomposition of the above form, where one of the rows of U is the eigenvector u, but that is not something we can immediately derive from the classical decomposition theorem.

The best strategy to solve this problem is to use simultaneous diagonalization.
For that, observe that A commutes with the matrix J, whose entries are all |G|^-1.
This follows from the regularity property.
On the other hand, for J, the multiplicity of the eigenvalue 1 is 1.
(The only other eigenvalue of J is 0 with a multiplicity of |G|-1.)
Now let us obtain a simultaneous diagonalization of A and J (this is possible because both are symmetric, real, and they commute with each other <cit.>)[The result has been formalized by Echenim <cit.>.]:

    A = U^T D U        J = U^T E U

Now, there is an index i s.t. e_ii = 1, corresponding to the u in U, because of the decomposition of J.
Thus, we obtain a sequence of eigenvectors of A, which necessarily contains the eigenvector u.

It should be noted that the same result could be obtained by introducing a correction to the diagonalization of A since the span of the eigenspace for the eigenvalue 1 will contain the eigenvector u.
However, the approach above is much more elegant. 



 §.§ Strongly Explicit Expander Graphs

The class of expander graphs that are needed for derandomization are strongly explicit expander graphs.
For such a graph, it is possible to sample a random vertex uniformly and compute the edges incident to a given vertex algorithmically, , it is possible to sample a random walk without having to represent the graph in memory.
This makes them very useful for the application in space-constrained settings.
This is in contrast to randomized constructions of expander graphs <cit.>, which is much easier but would require the representation of the graph in memory.[Edmonds and Paulson <cit.> have done foundational work on the formalization of probability theoretic existence proofs for (hyper-)graphs.]

The first strongly explicit construction was by Margulis <cit.>, who could show for a family of graphs that they have a spectral bound strictly smaller than 1.
But the bound itself was unknown.
A few years later, Gabber and Galil <cit.> managed to show, for a slightly different family of graphs, a concrete spectral bound, to be precise: λ_a(G_n) ≤5/8√(2).

The description of their family of graphs — now usually called the Margulis-Gaber-Galil family — is pretty elementary. For every n > 0, the graph G_n has n^2 vertices, which we label with the elements of ℤ/nℤ×ℤ/nℤ. Then each vertex (x,y) ∈ G_n is connected to the vertices:

    (x+1,y), (x-1,y) (x,y+1), (x,y-1), (x,y+x), (x,y-x), (x+y,y), (x-y,y)

The elementary description is, however, quite deceiving because the derivation of the spectral bound is far from it.
There were multiple efforts towards reducing the complexity of the proof.
Indeed, the initial proof by Gabber and Galil was improved by Jimbo and Marouka <cit.> in 1985.
Ravi Boppana further improved it in his lectures.
Hoory  <cit.> present a version of his proof.
But even for that last version, the authors admit that it is “subtle and mysterious”
It heavily relies on the discrete Fourier transform on the torus ℤ/nℤ×ℤ/nℤ, followed by a number of careful decoupling steps, a partial order, many case splits, and symmetry arguments. 

After that, there was more research.
For example, Lee <cit.> shows a weaker spectral bound for the same family but with a less complex proof.

The formalized proof <cit.> obtains the best-known spectral bound 5/8√(2) and follows the previously mentioned proof by Hoory 
(Although I could improve on their proof slightly by using one more symmetry argument to reduce the number of case splits.)

After Margulis's discovery, many more constructions have been found <cit.>.
An interesting solution was discovered by Reingold  <cit.> called the zig-zag construction, which works by starting with small fixed expander graphs and performing various graph operations to construct a family of larger expander graphs.
This would be easy to formalize, although the Margulis construction is much more elegant.
Moreover, the fact that the vertex counts of successive Margulis graphs do not grow by more than a factor of two is helpful, as we will see later in this section.

Other constructions rely on the Cayley graph of non-abelian groups <cit.>.
The proofs heavily rely on representation theory, which would require much more groundwork in Isabelle/HOL.
In particular, constructions of Ramanujan graphs <cit.>, which are expander graphs with the best-possible spectral bound (w.r.t. to their degree.), belong to this group.
I believe the formalization of such a construction is a significant open problem.



 §.§ Graph Operations

The Margulis-Gabber-Galil family consists of graphs for all vertex counts that are square. Using a method described by Murtagh  <cit.>, it is possible to obtain expander graphs for every vertex count.
This works by vertex identification, specifically given an expander graph G with degree d, a spectral bound λ_a and n vertices v_0,…,v_n-1, we identify distinct vertices (, vertex v_0 with v_m, v_1 with v_m+1 etc. and v_n-m-1 with v_n-1) for some m s.t. n/2≤ m ≤ n.
Note that edges are preserved, even if this introduces parallel edges.
If there is an edge between identified vertices, the edge becomes a self-loop.
To preserve regularity, we introduce d self-loops for vertices that have not been paired (for v_n-m,…,v_m-1).
This means the resulting contracted graph will have degree 2d.

Then, the resulting graph will have a spectral bound of at most λ_a(G)+1/2.

While the idea is due to Murtagh , they only show that contracting a graph with a spectral bound of at most 1/4 results in a graph with a spectral bound of at most 1/2.
Note that the result stated above, which is part of the formalization <cit.>, is a generalization.

Let us go through the proof, which can be verified using a generalization of Eq. <ref>.

    |x^T A x| ≤λ_a(G)  x ^2 + (1 - λ_a(G)) (x^T u)^2  for allx.

where u is defined as in Eq. <ref>.

Let x ∈ℝ^m, x  ≠ 0 and x ⊥ u' and let A denote the stochastic matrix of G and A' be the stochastic matrix of the contracted graph G'. 
Similarly u ∈ℝ^n, u_i = n^-1/2 and u' ∈ℝ^m, u'_i = m^-1/2. 
Moreover let y ∈ℝ^n s.t. y_i := x_i  m and β = |∑_i=n-m^m-1 x_i^2|. Then:


    2 |x^T A' x|     =    |∑_i,j=0^n-1 x_i  m a_ij x_j  m + ∑_i=n-m^m-1 x_i^2| ≤ |y^T A y| + β
       ≤ Eq. <ref>   λ_a(G) y^2 + (1-λ_a(G)) (y^T u)^2 + β
        =    λ_a(G) (y^2 + β) + (1-λ_a(G)) ((y^T u)^2 + β) 
        =     2 λ_a(G) x^2 + (1-λ_a(G)) ((∑_i=0^n-1 x_i  m n^-1/2)^2 + β) 
       =x ⊥ u'    2 λ_a(G) x^2 + (1-λ_a(G)) ((∑_i=m^n-1 x_i  m n^-1/2)^2 + β) 
       ≤ Cauchy-Schwarz     2 λ_a(G) x^2 + (1-λ_a(G)) ((∑_i=m^n-1 x_i  m^2)(∑_i=m^n-11/n) + β) 
        =     2 λ_a(G) x^2 + (1-λ_a(G)) (n-m/n(∑_i=0^n-m-1 x_i^2) + β) 
       ≤(n-m)/n ≤ 1    2 λ_a(G) x^2 + (1-λ_a(G)) x^2 .



Using the Margulis construction and the fact that there is a square between any n and 2n (inclusive), it follows that there is a strongly explicit expander graph with spectral bound (5 √(()2) + 8)/16 and degree 16 for every vertex count. 

Another graph operation is the power operation <cit.>: Let G be a graph, then the k-th power of G is a graph — denoted as G^k — consisting of the same set of vertices, where the number of edges between two vertices v and w in G^k is the number of paths of length k between v and w in G. Note that algebraically the stochastic matrix of G^k is A^k if A is the stochastic matrix of G.
It is easy to verify that λ_a(G^k) ≤λ_a(G)^k.

The previous two observations, in combination, imply that there exists a strongly explicit expander graph for every vertex count and any spectral bound λ_a>0 with a degree fulfilling: ln d ∈( ln (λ_a^-1)) <cit.>.
At this point, we are finally at a stage where we can introduce the second pseudorandom object, denoted by  l  P.
It first selects the strongly explicit expander graph constructed above with a number of vertices n, matching the size of the seed space of P and a spectral bound of at most .
Then, it samples random walks in that graph of length l.
Note that because of the regularity of the graph, there are exactly n d^l-1 random walks, and given a number between {0,…,nd^l-1-1}, we can output the vertices representing that path algorithmically.
Here, it is important that multiplicities are taken into account, which is also the reason some vertex tuples may have a multiplicity higher than 1.
   
Note that, for the users of the library, it is possible to ignore the complex internal structure of the pseudorandom object. Let us consider the classic Hoeffding bound for independent {0,1}-valued random variables:

theorem classicchernoffbound  assumes AE x in measurepmf p f x  01  defines  = expectation p f  shows w in prodpmf0l  p jl fw jllc2expnc2
If we instead wanted to choose expander random walks, we just have to substitute the sampling space and use the Chernoff bound for these <cit.>:


theorem expanderchernoffbound  assumes AE x in samplepro P f x  01  defines  = expectation samplepro P f  shows w in samplepro  l  P jl fw jllc2expnc2

Like before, the right-hand-side is expressed with respect to the probability space of the chained pseudorandom object: samplepro P, which allows the derivation of results with respect to more complex combined pseudorandom objects seamlessly. As mentioned before Section <ref> will demonstrate this on a concrete example.  



 §.§ Constructiveness

The pseudorandom objects presented are implemented in an entirely constructive way.
This means, for example, that the algorithms can be exported into executable efficient code <cit.> in various languages such as Haskell, Scala or SML. 

Here, it is essential to note that, for randomized algorithms, this doesn’t work directly because they are represented in Isabelle using the Giry-Monad. 
However, there is a framework introduced by Eberl and myself <cit.>, which enables exporting code for randomized algorithms.
This works by externalizing the source of randomness using a lazy stream of coin flips and establishing a Scott-continuous monad-morphism between the Giry-monad representation and the representation as a randomness consumer.
The resulting code can then be executed in the target language using any source of randomness (for example, a Hardware random number generator.)



§ COMBINING PSEUDORANDOM OBJECTS 

This section presents a concrete example where a combination of multiple pseudorandom objects are used to improve an algorithm's space complexity.
To achieve that, we are going to start with example2 from the introduction and improve its accuracy (relative error with respect to the correct value) and its probability of success, ,
the probability with which the algorithm returns a result with the desired accuracy. To make the presentation accessible, we will do this step-by-step and also briefly point to the 
main points of the correctness proofs. The examples are all also available in the AFP <cit.>. 
Let us recall:

fun example0pt2  nat  nat list  real pmf  where example0pt2 n xs      do 0pt      h  sample0ptpro  4 n  11      return0ptpmf x  xs h x2    
We had already established that the result has the expectation F_2 and a variance of at most 2F_2^2.
A standard technique to improve the accuracy of such an algorithm is to run it in parallel multiple times and obtain the average.

fun example0pt4  nat  nat list  real pmf  where example0pt4 n xs      do 0pt      let s  82      h  prodpmf 0s  sample0ptpro  4 n  11      return0ptpmf js x  xs h j x2s    
It is possible to show that the above algorithm returns an estimate of F_2 with relative error ε with probability at least 3/4, or equivalently the probability that the result does not have the desired accuracy is at most 1/4.

lemma example4correct  assumes set xs  0n  shows  in example4   F2 xs    F2 xs  14

Usually the first step of the verification is to express the result in terms of random variables over the probability space of the pseudorandom object(s):

  define s where s  82  define R where R h  js x  xs h x2s  let q  prodpmf 0s  sample0ptpro  4 n  11
With these definition it is now straightforward to check that the left-hand-side of the lemma equals:

  measure q h R h  F2 xs    F2 xs
which can be shown to be at most 1/4 using Chebyshev's inequality and
the facts that:
 
  
  * The expectation of R is F2 xs.
  
  * The variance is at most 2  F2 xs24.

Both results follow using the previously mentioned result about the expectation and variance of example2.
The first one follows directly using linearity of expectations. The second condition requires the application of Bienaymé's identity, , the variance of the sum of pairwise independent random variables equals the sum of the variances of the random variables:

lemma bienaymesidentitypairwiseindep2[<cit.>]  fixes f b  a  real  assumes finite I  assumes i i  I  f i  borelmeasurable M  assumes i i  I  integrable M  f i 2  assumes J J  I  card J  2  indepvars  borel f J  shows variance  i  I f i   i  I variance f i
The above solution requires ( ε^-2 (ln n + ln m) ) bits of space.
However, it is possible to improve on that considerably because Bienayme’s identity only requires pairwise independence.
Something we already know how to handle. Instead of choosing the seeds of the hash functions independently, it is enough to choose them pairwise independently, which can again be done using hash families.
The following is the improved version, where we replaced the product construction with the pseudorandom object for 2-independent hash families:

fun example0pt5  nat  nat list  real pmf  where example0pt5 n xs      do 0pt      let s  82      h  sample0ptpro  2 s  4 n  11      return0ptpmf js x  xs h j x2s    
The correctness proof for example5 and example4 are, of course, very similar.
The additional lemmas <cit.> needed in the proof include the fact that the outer hash family is pairwise independent:

lemma hashpropmfkindep  assumes isprimepower prosize S  shows probspacekwiseindepvars samplepro  k n S k  discrete i   i n
and the fact that mapping to a specific component is like sampling directly from the inner pseudorandom object.

lemma hashpropmfcomponent  assumes isprimepower prosize S  assumes i  n  assumes k  0  shows mappmf f f i samplepro  k n S  samplepro S

The latter theorem illustrates a common approach — when working with composed pseudorandom objects.
First a result is established for the inner pseudorandom object, then a theorem like the above can be used to lift it to a statement for the combined object. 
Indeed, the theorems about pseudorandom objects are specifically designed in a way that makes this possible. 
 
Note that the resulting algorithm now requires only ( ln n + ε^-2ln m ) bits.
This algorithm has a failure probability of 1/4. In the next step, we would like to improve the failure probability.

A straightforward method would be to increase the parameter s (the number of repetitions), , to obtain a mean of a larger sample.
Indeed, by setting s = ⌈ 2 ε^-2δ^-1⌉, the algorithm would return an ε-approximation with a probability of at least (1-δ). However, there is a much more efficient strategy.

Let us first observe that the median of a sequence will be an ε-approximation if more than half of the sequence elements are.
(Note that under the assumption, when we sort the sequence, the sequence elements that are ε-approximations will form a consecutive subsequence, whose length would be longer than half the length of the sequence. 
Thus, the subsequence would necessarily include the median.)
This is reflected in the following theorem <cit.>:

lemma medianest  assumes interval I  assumes 2card k k  n  f k  I  n  shows median n f I

In proofs the converse is more applicable, because of the upper-bound-approach mentioned before:

lemma medianestrev interval I  median n f I  2card k knf k  In

Now, what we can do is to repeat the algorithm example5 ⌈ 8 ln (δ^-1)⌉ times independently and obtain the median:

fun example0pt6  nat  nat list  real pmf  where example0pt6 n xs      do 0pt      let s  82 let t  8ln 1      h  prodpmf 0t  sample0ptpro  2 s  4 n  11      return0ptpmf median t i js x  xs h i j x2s    
Because the failure probability of example5 was 1/4, the probability that more than half of the t independent repetitions will fail (and hence the median) is at most exp(-t/8) ≤δ. This follows from the one-sided version of Hoeffdings’ inequality for {0,1}-valued indepedendent random variables.

The above leads to an algorithm with the space usage (ln(δ^-1) (ln n + ε^-2ln m)).
However, even the above can be improved using the second pseudorandom object discussed in Section <ref>: Random walks in expander graphs.
In Section <ref>, we pointed out that they fulfill a similar Chernoff bound. Hence, we can again replace the independent sampling with a pseudorandom object. The result is the following algorithm:


fun example0pt7  nat  nat list  real pmf  where example0pt7 n xs      do 0pt      let s  82 let t  8ln 1      h  sample0ptpro  t 18  2 s  4 n  11      return0ptpmf median t i js x  xs h i j x2s    
for which we can show:

lemma example7correct  assumes set xs  0n  shows  in example7   F2 xs    F2 xs  
Let us briefly go through the main parts of the formalized proof:

proof 
We set s, t as it is in the algorithm and introduce I the interval for the result:


  define s t where sdef s = 82 and tdef t = 32  ln 1/  define I where I =    F2 xs    F2 xs  define R where R h  js x  xs h x2s


The random variable R represents the estimates from example5 over the pseudorandom object:
 2 s  4 n  11.
Based on our previous observation, the random variable R will be outside the desired interval with probability at most 1/4.

  have measure  2 s  4 n  11 h R h  I  14

This can be stated as the expectation of an indicator variable:

  hence h ofboolR h  I  2 s  4 n  11  14

Let us introduce an abbreviation for the pseudorandom object of example7. It represents sampling using random walks from the pseudorandom object of example5.

  let q   t 18  2 s  4 n  11
Now we start an equation chain. First, we can express the left-hand-side as follows: 

  have L  measure q h median t i R h i  I
As mentioned earlier, the median can only be outside the interval if at least half of the estimates are outside:

  also have   measure q h j  t ofboolR h j  It  12
The following step is critical. Using the Chernoff bound for expander walks, we can conclude that the probability that half of the estimates will be outside the interval is less than exp(-2 l (1/8)^2). Here, we used the previous observation .

  also have   measure q h j  t ofboolR h j  It  14  18  18  also have   exp2t182
The right-hand side is less than δ, which concludes the proof.

  also have     finally show thesis by simpqed

Note that the resulting algorithm has a space complexity of (ln n + ε^-2ln(δ^-1) ln m).
In 1999 Alon  <cit.> had presented a solution for the estimation of the second frequency moment using (ε^-2ln(δ^-1) (ln n + ln m)).
That solution was relying only on a single hash-family.
This section demonstrates an asymptotically more space-efficient algorithm, and its construction was possible simply by using additional pseudorandom objects.

The above was a minimal example to demonstrate, how the library works and how combining psuedorandom objects can be used to improve a randomized algorithm.
In 2023, I had presented the formalization of an optimal space algorithm for the distinct elements problem <cit.>.
Since we are familiar with the notation, it makes sense to replicate the pseudorandom object I used there:


   m   l   2 n  nexpP  2 n  C7b2 P  k C7b2  b

A detailed explanation of all the parameters can be found in the source material <cit.>. However, this highlights that depending on the algorithm, very different
constructions may be required and points out the importance of the flexibility of the library. A pseudorandom object that is used in the above term, that is new is the primitive  n
which has a geometric distribution, i.e., the probability of the value 0 is 1/2, the probability of the value 1 is 1/4 etc.



§ CONCLUSION

This work presents pseudorandom objects as versatile components in randomized algorithms that can improve the asymptotic space and randomness complexity of randomized algorithms.
We discussed the mathematical foundations and the formal verification of two such objects in depth: expander graphs and k-independent hash families. 

Novel in this work is the viewpoint that they can be represented as a composable library that allows the construction of application-specific pseudorandom objects from primitive ones.
The notation introduced by this work allows the concise presentation of these algorithms informally and formally in Isabelle/HOL.
The high-level theorems provided by the library allow small and easy proofs (even if complex pseudorandom objects are needed).

Similarly, the formalization of higher-order finite fields and expander graphs has — as far as I know — not been done before.

The formalization of the optimal space distinct elements problem <cit.> highlights how this framework helps tackle advanced and complex randomized algorithms (that required many decades of research to develop.)

While this work has been done in Isabelle/HOL <cit.>, the approach is easily transferrable to other theorem provers/checkers/proof assistants such as Agda, Coq, Idris, Lean, etc.    

There is still much work that could be done in this field. In the context of expander graphs, there are still important constructions that could be formalized.
Another task is formalizing Chernoff bounds for random walks in expander graphs for non-discrete functions <cit.>.
Beyond that, further pseudorandom objects are open formalization tasks, such as randomness extractors and list-decodable codes. 



§ ACKNOWLEDGEMENTS

I would like to thank Tobias Nipkow and Yong Kiam Tan for feedback and many helpful suggesions and remarks on earlier versions of this manuscript.

plainurl




§ AFP ENTRIES

The following list summarizes the AFP entries that contain the results and definitions for pseudorandom objects, with a total of 16466 LoC.


  Finite Fields <cit.>
Formalization of finite fields and algorithms for efficient computations in and constructions of them. (7184 LoC)

  Interpolation Polynomials (in HOL-Algebra) <cit.>
Existence and count of interpolation polynomials for a given set of points and degree bound. This is used in the construction of hash families. (712 LoC)

  Universal Hash Families <cit.>
Formalization and construction of k-independent hash families as pseudrandom objects. Section <ref> summarizes this AFP entry. (1444 LoC)

  Expander Graphs <cit.>
Formalization and construction of (strongly explicit) expander graphs. Section <ref> summarizes this AFP entry. (7065 LoC)

The following list summarizes some AFP entries relevant for the verification of randomized algorithms:


  Median Method <cit.>
Concept of the median (and order statistics in general) and their properties as a random variable. It verifies that order statistics are measurable.
In particular, it contains the results about the median, we used in Section <ref>. (484 LoC)

  Concentration Inequalities <cit.>
Joint work with Yong Kiam Tan. Collection of advanced concentration inequalities, such as the Efron-Stein or Paley-Zygmund inequality. Includes also Bienaymé's identity, which we used in Section <ref>. (2019 LoC)

The following entries are verifications of randomized algorithms that rely on pseudorandom objects:


  Distributed Distinct Elements <cit.>
Verification of an optimal space cardinality estimation algorithm <cit.>. (4700 LoC)

  Formalization of Randomized Approximation Alg. for Frequency Moments <cit.>
Verification of several approximation algorithms for frequency moments on streams <cit.>. (3789 LoC)

  Approximate Model Counting <cit.>
By Yong Kiam Tan and Jiong Yang. Verification of a randomized approximation algorithm for the count of satisfying assignments to a boolean formula <cit.>.
It relies on the construction of a 3-universal XOR-hash-family. (See also <cit.>.)
 (5890 LoC)


