Zero-Knowledge Proofs: Verifying Computation and Preserving Privacy

Expanding Verifiability Beyond Merkle Trees in the Age of AI

From my previous post, you're already familiar with Merkle Trees as a powerful data structure for efficient and secure validation of contents. You know they achieve this by hashing data and then hashing those hashes up a tree, allowing for Merkle proofs that can verify data inclusion or consistency without revealing the entire dataset. This capability is becoming even more crucial as AI-generated content makes verifiable content more imperative. 

But what if you need to prove something more complex than just data inclusion? What if you need to prove a computation was performed correctly, or that you know a secret, without revealing that secret? This is where Zero-Knowledge Proofs (ZKPs) come into play, offering new dimensions of verifiability and privacy.

What are Zero-Knowledge Proofs?

A Zero-Knowledge Proof is a cryptographic protocol where a prover can convince a verifier that a statement is true, without revealing any information beyond the truth of the statement itself. Think of it like proving you're over 18 without showing your ID or revealing your name and address. ZKPs bring two main "primitives" or building blocks:

• Computational Integrity (Succinctness): They allow you to create proofs of computations that are significantly easier and faster to verify than to perform the original computation. This means the proof itself remains small, regardless of how complex the computation being proven is. Just as a Merkle proof is small compared to the original data, a ZKP is small compared to the computation it verifies.
• Zero-Knowledge (Privacy): They provide the option to hide parts of the computation (like sensitive inputs or even parts of the model) while still proving its correctness.

While generating ZK proofs can be very computationally intensive, advancements in cryptography, hardware, and distributed systems are making them feasible for increasingly complex computations. This expansion of capabilities opens up a vast "design space for new applications".

Programming ZKPs: A Shift in Mindset

Unlike traditional programming, which focuses on how to compute, programming ZKPs (often called circuits) focuses on defining a set of constraints. These constraints are mathematical rules that the computation must satisfy. For example, you might constrain that two secret numbers multiplied together equal a public number, without ever revealing the secret numbers.

The typical workflow for building a ZKP involves:

1. Writing the circuit: Defining the constraints of your computation.
2. Building the circuit: Compiling it into a binary form and WebAssembly.
3. Trusted Setup: A crucial pre-processing step that generates a proving key (for the prover) and a verification key (for the verifier). 
4. Generating the proof: Using your private inputs (the "witness"), the compiled circuit, and the proving key.
5. Verifying the proof: Using the verification key, the public output, and the generated proof.

Concepts like hash functions are fundamental in ZKPs, just as they are in Merkle Trees. However, ZKPs often use "ZK-Friendly hash functions" like Poseidon, which are optimized for use within ZKP circuits, offering significant performance gains compared to traditional hashes like SHA-256 due to their arithmetic-based implementation. 

Commitments, a cryptographic primitive allowing you to "commit" to a secret value without revealing it, are also crucial, often built using these hash functions. These are key building blocks for applications like digital signatures and more advanced concepts like group signatures, where you can prove you are part of a group without revealing your specific identity.

Programming ZKPs: An Example

Lets walk through a basic example of proving that we know two numbers whose product is 36 without revealing what those numbers are.

Write the circuit using Circom. c <== a * b; is the constraint, two numbers multiplied together equals a third number.

Circom compiles it into a Wasm file which we'll use to generate a witness for specifying our private inputs when creating the proof and a Rank 1 Constraint System binary file mathematically defining our single constraint.

Then we perform a trusted setup. The generated Common Reference String (CRS) consists of a proving key and a verification key. These keys can then be used every time we want to generate and verify proofs, respectively. They can be shared publicly and I provide mine here as part of the example:

• example1_verification_key.json
• example1_0001.zkey

Finally, we generate the proof using snarkjs with the Wasm file, proving key and private input that might be something like 9 and 4.

We get proof and public output JSON files.

We've proven that we know two secret values, a and b, whose product is 36. You can verify the proof (assuming you trust my verification key) with snarkjs.

$ snarkjs groth16 verify example1_verification_key.json public.json proof.json
[INFO]  snarkJS: OK!

If you change public.json to contain a different number the proof will no longer be valid. I've no longer proved I know the factors of this new number.

ZKPs and Blockchains, and Machine Learning (ZKML)

The convergence of ZKPs, blockchains (Web3), and machine learning is a rapidly advancing area with significant potential.

Blockchain use cases include:

• Scaling Blockchains: Public blockchains have limited computational power. ZKPs enable computations to be executed off-chain, with only a small ZK proof verified on-chain. This scales blockchains without sacrificing decentralization or security. Examples include ZK rollups like Polygon zkEVM and zkSync.
• Privacy-Preserving Applications: The zero-knowledge property is ideal for creating applications that protect users' privacy and personal data when making cryptographic attestations. Aztec Network, for instance, uses a ZK rollup for Ethereum where users' balances and transactions are completely hidden.
• Identity Primitives and Data Provenance: Projects like WorldID use ZKPs for privacy-preserving proof-of-personhood protocols, allowing a person to prove they are a unique human without revealing their identity.

ZKML is about applying ZK proofs to machine learning models, specifically focusing on the inference step. The core motivations for ZKML include:

• Verifying AI-Generated Content: With AI content becoming indistinguishable from human-created content, ZKPs can help determine that a particular piece of content was produced by applying a specific model to a given input.
• Privacy-Preserving Inference: ZKPs allow you to apply an ML model to sensitive data, where a user can get the result of the model's inference without revealing their input to any third party.

While proving something as large as current LLMs with ZKPs is not currently feasible, there's significant progress on creating proofs for smaller models. Teams are actively working on improving ZK technology, including specialized hardware and proof system architectures, to allow proving bigger models on less powerful machines in less time.

Summary

While Merkle Trees excel at verifying data inclusion and consistency, ZKPs extend this idea to verifying computations and knowledge with the added benefit of privacy. This makes them incredibly powerful for building the next generation of scalable and private applications on blockchains, especially as AI-generated content and privacy concerns continue to grow. The future of verifiable content, whether data or computation, is increasingly intertwined with these advanced cryptographic proofs.

Update

A few days after I published this I saw Opening up ‘Zero-Knowledge Proof’ technology to promote privacy in age assurance from Google showing that some well-known players are active in this space as well.

Sources

https://zkintro.com/articles/programming-zkps-from-zero-to-hero

https://world.org/blog/engineering/intro-to-zkml

Merkle Trees

Merkle Trees are a data structure that allows for efficient and secure validation of contents. They are typically implemented as binary trees. Given N pieces of data, they would have 2N nodes and log(N) height. Each leaf node is the hash of a piece of data and every non-leaf node is a hash of its children's hashes. With only a hash stored at each node Merkle trees have a small, predictable size compared to the data of a large N.

Examples of Merkle tree usage include detecting data inconsistencies between replicas in NoSQL databases, ensuring file integrity in distributed storage systems, and verifying blockchain transactions. We can easily create one with pymerkle.

└─9f42c047...

    ├──50fcd75a...

    │   ├──4c4b77fe...

    │   │   ├──e8bcd97e...

    │   │   │   ├──2215e8ac...

    │   │   │   └──fa61e3de...

    │   │   └──9c769ac2...

    │   │       ├──906c5d24...

    │   │       └──11e1f558...

    │   └──fed7af7d...

    │       ├──2b15ae18...

    │       │   ├──53304f5e...

    │       │   └──3bf9c81c...

    │       └──8007dd69...

    │           ├──797427cf...

    │           └──195f58bc...

    └──555da077...

        ├──85224a5c...

        └──5c889ef4...

A "Merkle proof" is a list of intermediate hashes from the path between a leaf node representing the data you want to prove and the root of the tree. Generating the proof is like a modified DFS. The beauty of this is that anyone can verify the data is included in the tree without the whole tree being revealed. 

{

    "metadata": {

        "algorithm": "sha256",

        "security": true,

        "size": 10

    },

    "rule": [

        0,

        0,

        1,

        0,

        0

    ],

    "subset": [],

    "path": [

        "53304f5e3fd4bcd20b39abdef2fe118031cc5ae8217bcea008dea7e27869348a",

        "3bf9c81c231cae70b678d3f3038f9f4f6d6b9d7adcf9b378f25919ae53d17686",

        "8007dd69b92a67ea6410098635fa8ba53c44a5994c7e5d92b99e27f0711c626f",

        "4c4b77fe3fc6cfb92e4d3c90b5ade42f059a1f112a49827f07edbb7bd4540e7b",

        "555da077fcadba1f23e0f2bfac8793e6a3c79a0d605902df34ab43d3e0fb487c"

    ]

}

Verifying the proof requires only calculating the root hash from the provided proof and the data being verified. If the calculated root matches the known root of the tree then the data is present in the tree. We are using less space and less compute than if we were iterating a list to check if data is present.

That was an inclusion or audit proof. We can also do a consistency proof. In an append-only tree we can verify earlier versions of the tree against later versions to make sure no tampering has occurred. The later version must include everything in the earlier version, in the same order, and all new entries come after old entries.

{

    "metadata": {

        "algorithm": "sha256",

        "security": true,

        "size": 10

    },

    "rule": [

        1,

        0,

        0,

        0,

        0

    ],

    "subset": [

        0,

        1,

        0,

        0,

        0

    ],

    "path": [

        "fa61e3dec3439589f4784c893bf321d0084f04c572c7af2b68e3f3360a35b486",

        "2215e8ac4e2b871c2a48189e79738c956c081e23ac2f2415bf77da199dfd920c",

        "9c769ac26f8d61ff40859e5201537845555136f0fd7ab604f7033180fbe76af9",

        "fed7af7d64bf0a73fcad018df1219928dbafa4d96b5d78f8a5e9be66ff0ada38",

        "555da077fcadba1f23e0f2bfac8793e6a3c79a0d605902df34ab43d3e0fb487c"

    ]

}

Given how powerful and pervasive Merkle trees are, I'm surprised they aren't discussed more along with other common data structures. It seems that with AI-generated content making verifiable content more imperative, their usage will only increase going forwards.

CUDA Kernels: Speeding up Matrix Multiplication

The Python ecosystem benefits greatly from be able to use libraries written in C/C++ and Rust (see my previous post) to increase performance. Increasingly, though, I've been running code on GPUs instead. Libraries like PyTorch simplify this, but how do they do it? Previously I could only answer something like "it's CUDA" without knowing what that really means. In this post we'll dive deeper and see what it takes to create our own CUDA kernel. We'll (re-)implement matrix multiplication, run it on a GPU, and compare to numpy and PyTorch performance.

First, what is CUDA? CUDA is a general purpose parallel computing platform and API. It gives direct access to NVIDIA GPU's virtual instruction set and parallel computation. It works with C, C++, Fortran, has Java and Python wrappers, and is supposed to be easier to use than earlier APIs like OpenCL. CUDA allows us to execute kernels which are functions compiled for GPUs, separately from the main program. This is how the majority of deep learning taking place today functions.

Second, what makes GPUs so special? GPUs are generally slower than CPUs at one operation but excel at operations that can be parallelized. They can have thousands of cores and higher memory bandwidth. If we want to tackle more than embarrassingly parallel problems, like matrix multiplication, this means we need a parallel algorithm to make use of the GPU.

Third, how do we access it from Python? It turns out there are several options, or several Python wrappers for this. I chose to go with Numba, a just-in-time (JIT) compiler that can target both CPUs and GPUs. It lets you write kernels directly using a subset of Python. It does not implement the complete CUDA API, but supports enough to tackle many problems. There are other (uninvestigated) options like CuPy and PyCUDA as well.

Lastly, before we get to the code, is the topic of writing efficient kernels. While this is beyond the scope of my post, I can say that I've at least learned it requires understanding several additional concepts beyond general concurrent programming. CUDA kernels are not for the faint of heart. It's not something you can just pick up in a couple of hours or from following one tutorial. One of the things you'll first notice, as a very simple example, is the need to explicitly move data back and forth between CPUs and GPUs. And kernels can't return values, you can only pass inputs and outputs.

For a better intro, I recommend reading the four-part CUDA by Numba Examples series (part 2, part 3, part 4).

Our baseline for this experiment will be numpy's highly-optimized matmul function. It supports vectorized array operations and some multi-threading by releasing the GIL per Parallel Programming with numpy and scipy. The underlying BLAS routines will be optimized for your hardware. This method of matrix multiplication has been tuned over the past several decades.

Here we'll try our own kernel. It's moving the computation from the CPU to the GPU, but it's likely lacking many optimizations that a battle-tested library has. Still, for medium-sized, two-dimensional matrices, we get some performance improvement.

Now we'll use PyTorch's matmul function. It's highly optimized like numpy and has the GPU benefit like the kernel. It's amazingly fast and works with larger, higher-dimensional matrices.

NVIDIA has a profiling tool called Nsight Systems that we can use to see GPU utilization. GPUs are expensive, we would want them to be fully utilized. From the reports, I see that the PyTorch implementation used more threads so that's consistent with higher parallelism. It also seems to have a higher ratio of memory operations vs. kernel operations. I'm not sure I understand what that means, but the kernel operations are sgemm which looks to be a matrix multiplication algorithm akin to numpy using BLAS.

Creating a CUDA kernel has become accessible enough that I could do it in a couple of hours on a laptop, yet my implementation remains far from the top implementations. Matrix multiplication is obviously common and great libraries exist. For less common operations, even if there's a known parallel algorithm, I would hesitate going the custom kernel route again. It's not a one-off thing you would just try to improve performance, it requires a way of thinking and deep optimization knowledge that most software developers don't have. The underlying libraries used by PyTorch are, for example, optimizing use of memory caches, shared vs. global memory accesses, thread utilization, and probably tuning kernel parameters for my specific GPU.

UPDATE:

NVIDIA Finally Adds Native Python Support to CUDA.

Rust to Python and Fibonacci Numbers

While Python continues to be the language of choice for ML projects, I'm increasingly seeing mentions of Rust in packages I use. Hugging Face fast tokenizers are written in Rust to speed up training and tokenization. They say it "Takes less than 20 seconds to tokenize a GB of text on a server's CPU." I recently used Polars, which is also written in Rust, for a quick task where I wanted to run a SQL query on 20 GB of CSV files. It only took about 60 seconds on my laptop. It's a way of bypassing Python's GIL to write code that is parallelizable. This is in contrast to packages like numpy that traditionally have been written in C/C++ for performance reasons. 

How does one turn Rust code into something you can call from Python? It turns out I've already done something similar, compiling Rust into WebAssembly and using it in a JavaScript project. For Python, the process is similar. I followed Calling Rust from Python using PyO3 and re-created the Fibonacci numbers experiment (apparently I'm not the only one with this idea) from my previous post. It's a toy example, just intended to show the ease with which Rust can be leveraged. It ignores more complex data types and any actual multi-threaded code. Perhaps that will be the topic of a future post.

The function in Rust looks like this:

And the Python code that calls it and times it:

Using the same n=35 as in my JavaScript experiment, I'm seeing about .05 seconds for Rust vs. 7.31 seconds for pure Python. The likelihood of wanting/needing this in a Python project seems greater than with JavaScript. My guess is that the trend continues.

Queueing Simulations

Waiting in lines is something we've all experienced. But how long will the wait for a table really be? How many checkouts should the grocery store have open? Queueing theory gives us the tools to be able to answer questions like these and more. It's also relevant to understanding software under heavy load. I've seen it come up recently in The Amazon Builder's Library Avoiding Insurmountable Queue Backlogs and Site Reliability Engineering's Addressing Cascading Failures chapter. 

Since you are reading this blog it is assumed that you are already familiar with queues and the idea of using them in software systems. To summarize, queues can improve availability and throughput at the expense of latency and resource consumption. Here we will see what the field of queueing theory can teach us about how and when to best use them. 

The Builder's Library article warns us of the bi-modal behavior of queue-based systems: the behavior can be very different based on whether there is a backlog or not. Recovery time after an outage can be dramatically increased and time can be wasted doing work that is no longer useful. We then get the first hints of how to manage this. Using more than one queue helps shape traffic. LIFO-ish behavior can be more desirable than FIFO. The SRE book discusses throttling via small queue sizes so requests are rejected when the incoming request rate is too high to be sustained. Traffic patterns, queue sizes, and the number of threads removing work from the queue are all intertwined. Can we better quantify these recommendations though?

I first learned of queueing theory in The Essential Guide to Queueing Theory e-book from VividCortex and wanted to revisit it in the context of the SRE reading. Right from the beginning we are warned that while queues are linear data structures, queueing systems behave non-linearly. Queueing theory is probabilistic. We won't understand the behavior exactly, but will know about wait times on average or a distribution of wait times.

The first thing they suggest we must understand is that queueing happens even when there is enough capacity to do the work. This is because work arrives to the queue in irregular sizes and at irregular intervals. Queueing gets worse at high utilization, when there is high variability, and when fewer workers are servicing the queues.

Any system can be decomposed into networks of queues and workers. These are the parameters and metrics commonly used to understand these systems:

• arrival rate to the queue (λ or A) - in a stable system this is the same as throughput
• average wait time or residence time in the queue (W or W_q)
• average post-queue service time (S)
• service rate (μ) - inverse of service time
• latency or residence time (R) - sum of wait time and service time
• worker utilization (ρ or U) - 0 to 1
• average number of requests waiting or in service concurrently (L or L_q)
• number of workers (M)

Little's Law states that concurrency is arrival rate times residence time: L = λR and L_{q =} λW_q

The Utilization Law states that utilization is throughput times service time: ρ = λS or ρ = λ/μ or with M workers ρ = λS/M

Kendall's Notation is a shorthand for describing queue systems. Normally we'll be working with M/M/m systems. The first M says events arrive randomly and independently (memoryless) meaning they are generated by a Poisson process. The second M says the service times are exponentially distributed. The last m is the number of workers or M from above. These are safe assumptions unless you know otherwise.

For a M/M/1 system R = S / (1 - ρ) which is a fast-growing curve (the non-linear mentioned above). Residence time is inversely proportional to idle capacity. By rearranging the equations above we can solve for other parameters as well. For systems beyond M/M/1 it gets more complicated. For M/M/m, for example, we need to bring in Erlang's C-formula.

To do that, and to see these laws in action, it's time to write some code. There appear to be several Python packages for simulating queueing systems. I will use Ciw for this simulation. We'll assume we have 1,000 requests per minute, several servers pulling work from a queue, and need enough capacity such that 2 servers could be offline without the system falling behind.

servers	utilization %	latency (sec)
3	83	.35
4	63	.18
5	50	.15

We can use this simulation to figure out what "levers to pull" to get the behavior we want. It could even become a Linear Programming optimization problem. Based on the arrival rate and service rate it's easy to see we need at least 3 servers, but not as easy to see how the latency would change going from 3 to 5 or what the minimum queue size is.

While not always intuitive, especially with multiple queues, there are ways to understand how queue-based systems will behave. If you are interested in the topic, I suggest reading the linked resources as they go into much more detail. If you have other examples of queueing theory in action, please share in the comments.

Java Updates Versions 9-19

The past few years I haven't written much Java code and when I did, it was Java 8. Many projects, it seems, have stuck with Java 8 which was released back in 2014. Per the roadmap, Java 8 is designated as LTS, but so are Java 11 and Java 17. In fact, Java 19 is available as of last month and many interesting features have been introduced in the past 8 years. This post is an overview of what's changed. The highlights, in my opinion, so we're up-to-date. I think it's enough to be interesting but not so much that it can't be picked up quickly if you have experience with older Java versions.

First, some name conventions. Java EE is now Jakarta EE. Definitely don't call it J2EE anymore. And since Java 11 Oracle JDK and OpenJDK are basically the same.

Tooling Updates

I won't go into too much detail here. If we are interested in using any of these they are explained and documented well elsewhere. For awareness:

• Java got a REPL called JShell for learning and prototyping interactively (9)
• Java Module System as a new package abstraction (9) - my impression is that this more for the JDK itself and some libraries while OSGi continues to be applicable for regular, modular apps
• Single-file programs can be executed directly and don't need the intermediate javac step (11)
• Multiple new garbage collection options for improved and more consistent performance (11)
• ZGC for low latency
• Epsilon GC no-op
• Linux container awareness for GC, thread pool sizes, etc. (11)
• Application class data sharing (CDS) for improved startup times and lower memory footprints (12) - this also seems to have implications for apps running in containers

API Updates

We want to start incorporating these into our code where applicable, so I created examples to help get used to some of these updates.

Private interface methods (9). Helps to encapsulate code in default methods and create more reusable code.

Variables can have implicit types, including in lambdas, to reduce the verbosity of code (10 and 11).

Switch expressions to simplify code and prepare for pattern matching in the future (12).

Text blocks as way to simplify code with multi-line strings (13).

Simple pattern matching (14). I was introduced to pattern matching when programming in Scala and this seems to continue a trend of Scala features making their way, in some form, to Java. It looks like more is coming in terms of pattern matching options.

Record keyword for immutable data classes (14). Getters, a public constructor, plus equals, hashCode, and toString methods are generated automatically. Lombok is still more flexible, but this is nice for simple cases.

Sealed classes for fine-grained inheritance control (15). Super-classes that are widely accessible but not widely extensible.

Paradigm Updates

For lack of a better name I'll call these paradigm updates as they relate more to programming models. 

Flow API as an implementation of the Reactive Streams Specification (9). This seems to be a way to get the specification interfaces into the JDK but not necessarily replace libraries with better implementations.

Vector API introduces vectorization to Java (16). It looks like it doesn't happen automatically and requires special code, so to be more useful I think it needs to be included in a common library like NumPy in Python.

Virtual threads and structured concurrency (19). The one-to-one kernel to user thread mapping is broken enabling easier asynchronous programming. Read/watch Project Loom: Revolution in Java Concurrency or Obscure Implementation Detail? The tl;dr is we'll still need a higher level of abstraction like reactive programming unless you want to relearn all the low-level concurrency structures.

Probabilistic Graphical Models

Deep learning and neural networks get a lot of (deserved) attention, but there is another class of ML models called Probabilistic Graphical Models (PGMs) that can also be used for inference and prediction. They have applications in fields such as medical diagnosis, image understanding, and speech recognition. Think decision making based on incomplete or insufficient knowledge.

More formally, PGMs use graphs to encode joint probability distributions as opposed to the more traditional ML approach of learning a function that directly maps input to a target variable. This post isn't a technical introduction though. Rather, it is more of an introduction-by-example and a summary of pgmpy's excellent notebooks.

Given a simple graph for flower type:

Our two approaches would look something like this:

Bayesian networks

In this section I'll use a more complex graph for student grades:

For problems with many features and/or high cardinality features, inference will be difficult because the size of the joint probability distribution increases exponentially. PGMs can compactly represent it by exploiting conditional independence. They provide us efficient methods for doing inference over these joint distributions.

In this graph we have cardinalities of 2 for each node except Letter which is 3. The joint distribution would require storing 48 values (2*2*2*2*3) while the PGM only requires 26 (see notebook 1 for details).

This is what's known as a Bayesian network, which is always represented as a directed acyclic graph. Each node is parameterized by a conditional probability distribution (CPD) like P(node|parents). For example, the Grade node has the CPD P(G|D,I). Bayesian networks are used when you want to represent causal relationships between random variables. Naive Bayes is a special case where all random variables are assumed to be independent of each other, each only directly affecting the target variable.

Given tabular data and a graph structure, CPDs can be estimated using Maximum Likelihood Estimation (MLE). It's similar to what was done with the Iris data in the first code block above. It's also fragile because it is so dependent on the amount and quality of the observed data (see notebook 10 for details). This explains why that code breaks with some random seeds. 

A better solution is Bayesian Parameter Estimation. There you start with CPDs based on your prior beliefs (or uniform priors) and update them based on the observed data.

One method of exact inference in PGMs is variable elimination. It efficiently avoids computing the entire joint probability distribution (see notebooks 2 and 5 for details). For larger graphs there are other, approximate algorithms because an exact solution would be intractable.

Making predictions is similar. Instead of getting a distribution we get the most probably state.

Markov networks

Markov networks are represented by undirected graphs. They represent non-causal relationships. They can, however, represent dependencies that a Bayesian model can't, like cycles and bi-directional dependencies. Factors describe connected variable affinity, or how much two nodes agree with each other. The joint probability distribution is the product of all factors.

A quick note because the names sound similar. Markov chains are not PGMs because the nodes are not random variables. They can be represented as as Bayesian networks and PGM algorithms would be available.

Sampling

Sampling algorithms approximate exact inference by generating a large number of samples that will converge to the original distribution. One of these is Hamiltonian Monte Carlo. It is a Markov Chain Monte Carlo (MCMC) that proposes future states in the Markov Chain using Hamilton dynamics from physics (see notebook 8 for details). Other MCMC algorithms you may encounter are Metropolis-Hastings and Gibb's Sampling. See Monte Carlo Approximation Methods: Which one should you choose and when? for a comparison of these methods.

Another interesting find that fits in at this point is the PyMC3 library and the Probabilistic Programming and Bayesian Methods for Hackers open source book. 

I also think this is a nice writeup on Bayesian Logistic Regression using Pyro, another probabilistic programming library, and MCMC.

Learning networks

Learning a Bayesian network can be done as an optimization problem by scoring networks on how well they fit a data set, and searching through the space of all possible models. For non-trivial graphs where an exhaustive search is not possible, hill climbing can be used (see notebook 11 for details).

Wrap-up

I've only scratched the surface here, but I think it's a more intuitive introduction to the topic than most of the material in this space. We could build up to more complex graphs and problems from here.

Hidden Markov Models Explained with a Real Life Example and Python code is a nice extension, for example.

And to bring things full circle on where PGMs fit in to the ML landscape, here is an opinion from well-known ML researcher Ian Goodfellow:

The two aren’t mutually exclusive. Most applications of neural nets can be considered graphical models that use neural nets to provide some of the conditional probability distributions. You could argue that the graphical model perspective is growing less useful because so many recent neural models have such simple graph structure … These graphs are not very structured compared to neural models that were popular a few years ago like … But there are some recent models that make a little bit of use of graph structure, like VAEs with auxiliary variables.

Plus a tweet from the Standford NLP group:

There is increasing convergence between this decade’s neural models and last decade’s probabilistic graphical models… https://t.co/MKlIP4Sayx

— Stanford NLP Group (@stanfordnlp) May 15, 2018

Thus is would seem that knowing these concepts will continue to be useful even if we don't directly use PGMs or focus solely on PGMs.

ML at Scale Part 3: distributed compute

In part 2 I focused on ML when your data won't fit in memory. This post will move on to slow, or compute bound ML instead. I'll continue to use Dask and explore how it can help us.

Dask leverages multiple CPU cores to enable efficient parallel computation on a single machine. It can also run on a thousand-machine cluster, breaking up computations and routing them efficiently. The Comparison to Spark documentation is a great reference for understanding Dask in the context of an older tool and that older post. Perhaps the most interesting difference is Spark just being an extension of the MapReduce paradigm and Dask being able to implement more sophisticated algorithms by being generic task scheduling-based.

Sticking with the scikit-learn examples, here is replacing a parallel algorithm's Joblib backend with Dask to (potentially) spread work out across a cluster.

The Dask documentation is quick to point out in their best practices that not everyone needs distributed ML as it has some overhead. Compiling code with Numba or Cython could help, as could intelligently sampling some of your data. In this post I got huge speedups by vectorizing some code that was looping through large matrices.

Across this three part series we've now seen how to speed up reading large datasets, work with datasets that don't fit entirely in memory, and distribute processing across multiple machines. There's obviously a lot more to this, but I wanted to develop a better intuition for how to approach these types of issues and at least know where to start. Hopefully you learned something too.

UPDATE:

A specific tool isn't supposed to be the focus of this post. Dask was used here to illustrate the idea and show how simple it can be, but there are other options in this space. Here are a couple of other examples that I've come across:

• tune-sklearn, which is part of the Ray framework for building distributed applications, can do distributed hyperparameter tuning
• Distributed training with TensorFlow shows how you can distribute model training within the TensorFlow ecosystem

UPDATE 2:

Both Pandas 2.0 and Polars now use Apache Arrow as a memory model. Polars, a relatively new entrant to this space, is implemented in Rust and exposes a Python API. It is created specifically for fast data processing, not ML, but overlaps enough with this series that it is worth checking out.

UPDATE 3:

Since I mentioned Numba earlier, I'll also mention Jax another, newer library in that space. Jax has a NumPy-like interface, works with GPUs, offers JIT compilation, and supports automatic differentiation, vectorization, and parallelization. Check out their simple NN example to see it in action. It's not replacing distributed compute or even competing with Dask, the point is more that there are now a lot of amazing tools that make working with data easier and faster.

ML at Scale Part 2: memory

When data we want to train a machine learning model on becomes too big to fit in memory, we need to find a way to work on subsets of the data. I hinted at this in part 1. We can use Pandas to read chunks of a file, but that is fairly primitive and slow.

Libraries like Vaex and Dask attempt to abstract this away. 

Vaex provides lazy, out-of-core (not all in memory at once) DataFrames via memory-mapping. Pre-processing and feature engineering are more efficient, and memory is freed up for model training. It also has a vaex.ml package which provides a scikit-learn wrapper.

Dask provides large, parallel DataFrames composed of smaller Pandas DataFrames. This helps with data too big to fit in memory because the individual Pandas DataFrames can be stored on disk. DaskML provides estimators designed to work with Dask DataFrames.

The accuracies both came out to 96%. Similar ideas, different implementation. In fact, these DataFrames remind me a little bit of the persistent data structures covered in my Exploring Immutability post.

Even with these fancy DataFrames, many machine learning algorithms are designed train on all the data at once. If our data is too big to fit in memory then that's going to be a problem. As in the code above, we need to use online learning or incremental algorithms to solve this problem. The Incremental and IncrementalPredictor classes handle "streaming" the data (in batches) and we specify upfront the possible classes.

DaskML adds several generalized linear model implementations. scikit-multiflow is designed for actual streaming data and adds several other online learning algorithms. Neural networks are also trained in this manner, often being fed "mini-batches", so they are good candidates for datasets that don't fit in memory as well.

Stay tuned for part 3 where I'll get into being compute-constrained instead of memory-constrained.

UPDATE:

There's a lot of I/O and memory-related work coming out of the TensorFlow community as well. Checkout Better performance with the tf.data API as it overlaps nicely with Parts 1 and 2 of this series.

ML at Scale Part 1: I/O

I recently came across a somewhat large dataset from a Kaggle competition where the data was provided as an approximately 6 GB CSV file. A frequent comment in the discussion forum was how long it took just to read this file. A few GBs is large enough where this starts to become noticeable, but it's not really that big. If a laptop can have a 1+ TB drive and 32 GB memory then this isn't even in the realm of "big data". That's good, though, because it means there are some simple tricks we can use to cut down on that read time.

The pandas read_csv() method takes about 63 seconds for me. That's our baseline.

First we try reducing the precision. This gets us to 59 seconds. Not great.

Next we try reading the files in chunks. This is actually slower, but the technique could help us if the file didn't fit entirely in memory. More on that in future posts.

Then we try Dask which will spread the work across multiple processers. 30 seconds. Better, but I still don't want to wait that long. More on Dask in future posts as well.

Finally we convert the CSV file to a different format. I tried Apache Parquet but there are others. It's a binary columnar format (remember Column-oriented Database Basics?). Stored in this manner the data is 2.5 GB. And this gets us to just 3 seconds for reading the whole file!

Converting our data to the binary file format and possibly reducing the precision or using Dask as well would really shorten our feedback loop while training a ML model. It would seem that any cleaning of the data or preprocessing that we can do ahead of time would make sense to do once, before converting the file format, when the data is this size.

UPDATE:

Apache Arrow is another project to checkout in this space along with memory-mapped files. Reading this data in the Feather file format is even faster than Parquet.

Rust to WebAssembly and Fibonacci Numbers

WebAssembly (wasm) is a virtual assembly language. It's a binary instruction format that can be executed at near native speeds by JavaScript engines like V8, which means it can run in a browser or Node.js app. Not intended as a JavaScript replacement, it instead works with it for performance critical pieces of code. You can make calls from JavaScript to WebAssembly and vice versa. It can be useful for games, VR, and AR, for example. 

Rust is a fast and memory-efficient programming language with interesting type and thread-safety characteristics. I've never used it, but have wanted to check it out for a while. Along with Go and C, it can easily be compiled into WebAssembly. This post just touches the surface of Rust.

First, to learn how to use Rust in a JavaScript app, I followed Compiling from Rust to WebAssembly. They guide you through compiling Rust code, generating a JavaScript wrapper package, then using npm and webpack to run it. The end user only needs a modern browser and they can execute the code originally written in Rust none the wiser.

Building on that, as a simple experiment to see how much faster Rust code compiled to WebAssembly can be versus vanilla JavaScript, I compared calculating the nth Fibonacci number in each. I used the inefficient recursive algorithm intentionally, wanting it to be slow.

The function in Rust looks like this:

And the JavaScript code that calls it, times it, and compares it:

Surprisingly, although maybe it shouldn't be, the WebAssembly version of this simple function is orders of magnitude faster when executed in the browser. At n=35, where they really start to diverge, I'm seeing about .06 seconds versus 3.5 seconds. 

I don't often have to write CPU-intensive code like this in JavaScript, but with such a stark performance difference possible it's a good tool to have in the toolbox. 

GraphQL API Gateway Prototype

I've been meaning to check out GraphQL for a while. At work I'm seeing calls to REST APIs being made to get one small piece of data, most of the response discarded. Or, multiple HTTP requests, usually sequential, whose responses must be combined to do something useful with the retrieved data. Their granularity is wrong for the use case, but maybe not for someone else's. Can GraphQL help with this? It seems like it can. Model the business domain as a graph, like our mental models and object-oriented programming. Engines for many languages. A type system to enable good developer tools. It sounds great.

At the risk of adding another network hop, another moving, part, another layer, it would be nice to be able to call something to get exactly the data I need without changing the existing APIs. Mobile devices or slow internet connections could benefit from the reduced number of round trips. This extra layer could abstract away the different (read poor) uses of HTTP status codes and different response body styles. What I think I really want is a GraphQL API gateway.

The API gateway is a single entry point for all clients (the backends for frontends variation of the pattern has a different gateway optimized for each type of client). Requests can be proxied straight through to a single microservice or fanned out to several microservices. Responses can be aggregated and/or modified. This is the overlap point with GraphQL and why I think they would go well together. The API gateway can also centralize several cross-cutting concerns like throttling, routing, circuit-breaking, input validation, authentication (authorization stays in the business logic), etc.

Apollo is the GraphQL implementation I went with to prototype this. From their tutorial I started with the rocket launch API and extended it with a made up weather API to see how the two could be chained together. To the client, it's seamless. Both "services" are part of the same graph.

This is a lot of code to show in one shot, but I'll explain it below and then show some example GraphQL queries.

First, typeDefs defines the schema. You can query a list of rocket launches or a single launch by ID. dataSources specifies, of course, where the data is coming from, like a database or REST API. resolvers stitches these together. Notice how the Weather type takes a site from its parent, a Launch. This is how they are linked, or chained together.

With the Apollo server running you can try it out at http://localhost:4000/ in a browser. The GraphQL queries are on the left and the responses are on the right.

It's cool to see the different responses without having had to write any code to specifically handle them. A natural extension to this prototype would be to add mutation resolvers so clients could also update the graph.

Finally, the ThoughtWorks Tech Radar cautions against trying to create a universal, canonical, centralized data model. I think the bounded context ideas from DDD would apply. In their zero trust architecture blurb they also mention that a network perimeter isn't a security boundary anymore. That makes me question thinking of the API gateway as a place to shift all those cross-cutting concerns to. Do users have to go through the gateway or can they hit microservices directly? They mention service mesh as a solution and that would seem to be an API gateway alternative, but they aren't mutually exclusive either.

Column-oriented Database Basics

This is a short post on column-oriented databases. I'll barely scratch the surface, but among the types of NoSQL databases—document, key-value, column-oriented and graph—I've always thought column-oriented was the most difficult to wrap my head around. Hopefully we can get past that initial hurdle here, run a few queries, and they will seem like less of a mystery.

A relational database is optimized for retrieving rows of data. This works well for transactional applications. A column-oriented database is optimized for retrieving columns of data. This works well for analytical applications and some queries, like aggregations, become really fast because much less data needs to be read from disk to retrieve the whole column. There seems to be a lot of overlap between column-oriented databases and data warehouses.

A relational database would store 3 rows of data like this:

1:a,b,c;2:d,e,f;3:g,h,i

While a column-oriented database would store the same data like this:

a:1,d:2,g:3;b:1,e:2,h:3;c:1,f:2,i:3

For a low-cardinality column, compression algorithms work very well. Something like a:2,a:3 becomes a:2,3. In some ways this is like normalizing a relational database to reduce data duplication, but I don't think that enables the same level of compression or gives you the same data locality benefits.

Of course column-oriented databases aren't good for all workloads. They aren't optimized for queries that touch many fields. Writes can also be slow since they aren't just appending to the end of a file.

So where do they fit into a big data architecture? When thinking about this I remembered the article How to Move Beyond a Monolithic Data Lake to a Distributed Data Mesh. It avoids mentioning specific products, but they way I interpret it is that their first generation was dedicated data warehouse products where you did ETL. The second generation was more ELT, maybe using Hadoop and Spark. The third generation, the eponymous data mesh, unifies batch and stream processing, perhaps adding Kafka to the above mix. Using the CQRS pattern, for example, a column-oriented database could fulfill some of the Q (query) operations as a sort of data warehouse.

For my first column-oriented database adventure, after that background research, I chose Amazon Redshift and I followed their Getting Started with Amazon Redshift guide. It only took about an hour to spin up a Redshift cluster, upload their sample data to an S3 bucket, copy that into Redshift, then run a few SQL queries.

The big surprise and takeaway was that this felt just like using a relational database. The data I uploaded was in text files (basically CSV files) and was used to create tables. The queries I ran were SQL queries that joined tables, selected different fields and aggregated the results. The details of how the database works under the covers is interesting, but doesn't inform its usage. There isn't much mystery after all.

SELECT * FROM event;

SELECT firstname, lastname, total_quantity FROM (SELECT buyerid, sum(qtysold) total_quantity FROM sales GROUP BY buyerid ORDER BY total_quantity desc limit 10) Q, users WHERE Q.buyerid = userid ORDER BY Q.total_quantity desc;

A couple of other column-oriented databases you may have heard of are Cassandra and HBase. NoSQL databases are built to scale horizontally so you also need to consider how the CAP Theorem applies to your situation when choosing among them as they each make different trade-offs in addition to their unique features. The best choice will be highly data and workload specific.

Dynamic Programming

Optimal substructure means that the optimal solution can be constructed from optimal solutions of sub-problems. Think recursion. When we combine optimal solutions to non-overlapping sub-problems the algorithmic technique is called divide and conquer. Think of a merge sort where each sub-sort can happen independently and they are combined at the end. When we combine optimal solutions to overlapping sub-problems the algorithmic technique we can use is dynamic programming. Because the sub-problems overlap, we can improve the running time by remembering the results of already computed sub-problems and not computing them again.

An easy way to visualize this is in calculating Fibonacci numbers. The standard recursive solution can be memoized to significantly improve time complexity at the cost of using more space to store sub-problem results. This is top-down dynamic programming—recursion and memoization. Top-down dynamic programming can be easier to program and anecdotally it is more commonly used.

Bottom-up dynamic programming is instead iterative and we solve sub-problems first, building them into bigger sub-problems. Without the overhead of recursive calls we don't necessarily increase space complexity, but still have the decrease in time complexity. For calculating the nth Fibonacci number this looks like:

This is time complexity O(n) and space complexity O(1), a huge improvement on the naive recursive solution that's O(2ⁿ) time and O(n) space. For details on big-O of recursive algorithms read this.

Some other interesting applications of dynamic programming are:

• Dijkstra's algorithm (shortest paths)
• Tower of Hanoi (teaching game)
• Knapsack problem (combinatorial optimization)

The knapsack problem, for example, has a top-down solution but I think the bottom-up solution is especially appealing. Using a matrix to store sub-problem solutions we can make the O(2ⁿ) time recursive algorithm O(nW) time and space:

But wait, there's more. We only need to remember a part of the matrix for the next iteration, which means we don't even need the matrix. It can be further optimized to only use O(W) space:

Not, in my opinion, an obvious solution by any means, but a very elegant one.

Quantum Teleportation

I recently checked out the IBM Quantum Experience and was really impressed with the content available for getting started with quantum computing. There are simulators for executing quantum circuits, runnable both locally and in the cloud, but the highlight is being able to submit jobs to run on a real quantum computer. To create these circuits we use the Python DSL Qiskit.

For gaining a better understanding of Qiskit I recommend starting with the Coding with Qiskit video series. Following along with the videos and coding my first circuit then led down a rabbit hole of trying to understand what I was really building and how it could be useful.

Quantum computing uses phenomena from quantum mechanics like superposition and entanglement to perform computation. In classical computing we have bits which can either be on or off, 1 or 0. In quantum computing we have qubits which can be on or off at the same time, or somewhere in between. This is superposition. I like the coin analogy from Quantum computers will change the world (if they work) of this being a coin flip (classical) versus a spinning coin (quantum).

A qubit is represented as a vector of two complex numbers with length 1 where these represent the probability of a 0 or 1 state. Measurements destroy quantum state, collapsing the qubit to a classical bit. A measurement cannot be reversed. A quantum state cannot be copied like a classical bit's state can. Instead, quantum teleportation is the transfer of state from one qubit to another via a linkage known as entanglement. Even if moved far apart, entangled particles transmit information to one another. One coin flip mirroring another, linked flip.

The probabilistic nature of quantum computing enables many calculations to be processed simultaneously, thus making some intractable problems tractable. Quantum algorithms use superposition or entanglement to solve problems like factorization or search faster than classical algorithms. Qiskit has some amazing Jupyter notebooks showing how these can be used. Three I thought looked the most interesting, for example:

• quantum-enhanced support vector machines
• traveling salesman problem
• portfolio diversification

My first circuit is much simpler, of course, but I feel good about getting it running. Based on the videos linked above, I teleport a 1 state from one qubit to another and then measure it to confirm. I simulate this locally, in the cloud, and then run it on one of IBM's quantum computers. The actual quantum computer run includes some noise which is from the imperfect nature of today's quantum computers.

UPDATE:

TensorFlow Quantum is another library that could be used if you want to leverage Google's quantum computing resources.