What is the Cauchy combination test?

Set-up

Assume that we have p-values . Assume that they are computed from z-scores (test statistics following normal distributions). Let and let . Without loss of generality, assume that each test statistic has variance 1. With this, we can express the p-values as

where is the CDF function of the standard normal distribution.

We are interested in testing the global null hypothesis .

The Cauchy combination test

Assume that we have some random vector such that for all , and is independent of . The test statistic for the Cauchy combination test, proposed by Liu & Xie 2020 (Reference 1), is

Under the global null, for each , implying that has the standard Cauchy distribution (see this previous post). If the p-values are independent, then Proposition 1 in this other previous post implies that has the standard Cauchy distribution.

In turns out that even if the p-values are dependent, is still approximately Cauchy! This approximation is formalized in Theorem 1 of Reference 1:

Theorem. Suppose that for any , follows a bivariate normal distribution. Suppose also that . Let be a standard Cauchy random variable. Then for any fixed and any correlation matrix , we have

The theorem says that the test statistic has approximately a Cauchy tail even under dependency structures in . Knowing the (approximate) distribution of under the global null allows us to use it as a test statistic.

Other notes

• The “bivariate normal distribution” condition is a bit of a technical assumption, the authors claim it is a mild assumption.
• This result bears a lot of similarity with a result by Pillai & Meng 2016 (see this previous post for a description). Section 2.2 of Reference 1 discusses the similarities and the differences.
• Theorem 1 above holds for fixed (number of p-values). Section 2.3 of Reference 1 has a high-dimensional asymptotic result where with .
• The Cauchy combination test is especially powerful when only a small number of ‘s are large, or equivalently when a smaller number of ‘s are very small. We can see this intuitively: small ‘s become very large ‘s, so the test statistic will be dominated by a few very large p-values. See Section 4.2 of Reference 1 for a power comparison study.

References:

1. Liu, Y., and Xie, J. (2020). Cauchy Combination Test: A Powerful Test With Analytic p-Value Calculation Under Arbitrary Dependency Structures.

An interesting result involving the Cauchy distribution

If are independent variables and is a random vector independent of the ‘s with for all and , it is well-known that also has a distribution.

It is also well-known that if and are independent standard normal variables, then has a distribution.

Putting these two facts together, we have the following proposition:

Proposition 1. Let and be i.i.d. , where is a diagonal matrix with strictly positive entries on the diagonal. Let be a random vector independent of such that for all and . Then

Interestingly, the result above holds for any covariance matrix with positive entries on the diagonal! This was proved in Pillai & Meng 2016 (Reference 1). In other words, we can have any dependence structure between the elements of (and likewise for , but as long as the dependence structure for the two random vectors is the same, the linear combination of their ratios remains Cauchy-distributed! Here is the formal statement:

Proposition 2. Let and be i.i.d. , where is any covariance matrix with strictly positive entries on the diagonal. Let be a random vector independent of such that for all and . Then

References:

1. Pillai, N. S., and Meng, X.-L. (2016). An Unexpected Encounter with Cauchy and Lévy.

Converting a uniform distribution into a Cauchy distribution

The following proposition demonstrates how one can transform a uniform distribution on into a standard Cauchy distribution:

Proposition. If and , then has standard Cauchy distribution.

Proof: For any ,

where the second-last equality is due to the fact that , so the LHS inside the probability is always in . In other words, has the CDF of the standard Cauchy distribution.

We could also differentiate the CDF above w.r.t. to get the PDF of , which is more easily recognizable as the PDF of the standard Cauchy:

(Alternatively, one could use the formulas in this previous post with minor modifications since the transformation here is monotonically decreasing, not increasing.)

What is Adam and AdamW?

Adam (short for adaptive moment estimation) was proposed by Kingma & Ba (2014) (Reference 1) as a first-order gradient-based optimization method that combines the advantages of AdaGrad (Duchi et al. 2011, Reference 2) (ability to deal with sparse gradients) and RMSProp (Tieleman & Hinton 2012, Reference 3) (ability to deal with non-stationary objectives). While considered an old method by now, it’s still widely used as a baseline algorithm to compare more cutting-edge algorithms against.

According to the authors,

“Some of Adam’s advantages are that the magnitudes of parameter updates are invariant to rescaling of the gradient, its stepsizes are approximately bounded by the stepsize hyperparameter, it does not require a stationary objective, it works with sparse gradients, and it naturally performs a form of step size annealing.”

Notation

• is a noise objective function that is differentiable w.r.t. to its parameter . We want to minimize .
• are stochastic realizations of at timesteps .
• Let denote the gradient of w.r.t. .

Adam

At each timestep , Adam keeps track of , an estimate of the first moment of the gradient , and , an estimate of the second raw moment (uncentered variance) of the gradient ( represents element-wise multiplication). The Adam update is

where and are hyperparameters. The full description of Adam is in the Algorithm below:

Adam algorithm (Algorithm 1 of Kingma & Ba).

Comparison with other methods

In comparison with Adam, the vanilla SGD update is

and the classical momentum-based SGD update is

The basic version of AdaGrad has update rule

AdaGrad corresponds to Adam with , infinitesimal , and replacing with (see Section 5 of Reference 1 for a few more details).

AdamW

AdamW, proposed by Loshchilov & Hutter 2019 (Reference 4), is a simple modification of Adam which “improves regularization in Adam by decoupling the weight decay from the gradient-based update”.

The original weight decay update was of the form

where defines the rate of weight decay per step. Loshchilov & Hutter (2019) note that for standard SGD, this weight decay is equivalent to L2 regularization (Proposition 1 of Reference 4), but the two are not equivalent for the Adam update. They present Adam with L2 regularization and Adam with decoupled weight decay (AdamW) together in Algorithm 2 of the paper:

Overall, the authors found that AdamW performed better than Adam with L2 regularization.

PyTorch’s implementation of AdamW can be found here.

References:

1. Kingma, D. P., and Ba, J. (2014). Adam: A Method for Stochastic Optimization.
2. Duchi, J., et al. (2011). Adaptive Subgradient Methods for Online Learning and Stochastic Optimization.
3. Tieleman, T., and Hinton, G. (2012). Lecture 6.5 – RMSProp, Coursera: Neural Networks for Machine Learning. (If someone has a link/PDF for the slides, I’d be grateful!)
4. Loshchilov, I., and Hutter, F. (2019). Decoupled Weight Decay Normalization.

What is elastic weight consolidation?

Elastic weight consolidation (EWC) is a technique proposed by Kirkpatrick et al. (2016) (Reference 1) as a way for avoiding catastrophic forgetting in neural networks.

Set-up and notation

Imagine that we want to train a neural network to be able to perform well on several different tasks. In real-world settings, it’s not always possible to have all the data from all tasks available at the beginning of model training, so we want our model to be able to learn continually: “that is, [have] the ability to learn consecutive tasks without forgetting how to perform previously trained tasks.”

This turns out to be difficult for neural networks because of a phenomenon known as catastrophic forgetting. Assume that I have already trained a model to be good at one task, task A. Now, when I further train this model to be good at a second task, task B, performance on task A tends to degrade, with the performance drop often happening abruptly (hence “catastrophic”). This happens because “the weights in the network that are important for task A are changed to meet the objectives of task B”.

Let’s rewrite the previous paragraph with some mathematical notation. Assume that we have training data and for tasks A and B respectively. Let . Our neural network is parameterized by : model training corresponds to finding an “optimal” value of .

Let’s say that training on task A gives us the network parameters . Catastrophic forgetting says that when further training on task B (to reach new optimal parameters ), is too far away from and hence has poor performance on task A.

Elastic weight consolidation (EWC)

The intuition behind elastic weight consolidation (EWC) is the following: since large neural networks tend to be over-parameterized, it is likely that for any one task, there are many values of that will give similar performance. Hence, when further training on task B, we should constrain the model parameters to be “close” to in some sense. The figure below (Figure 1 from Reference 1) illustrates this intuition:How do we make this intuition concrete? Assume that we are in a Bayesian setting. What we would like to do is maximize the (log) posterior probability of given all the data :

Note that the equation above applies if we replace all instances of with or . Assuming the data for tasks A and B are independent, let’s do some rewriting of the equation above:

Thus, maximizing the log posterior probability of is equivalent to maximizing

Assume that we have already trained our model in some way to get to which does well on task A. We recognize the first term above as the log-likelihood function for task B’s training data: maximizing just the first term alone corresponds to the usual maximum likelihood estimation (MLE). For the second term, we can approximate the posterior distribution of based on task A’s data as a Gaussian distribution centered around with some covariance. (This is known as the Laplace approximation). Simplifying even further: instead of the distribution having precision matrix equal to the full Fisher information matrix , we assume that the precision matrix is diagonal having the same values as the diagonal of .

Putting it altogether, maximizing is equivalent to minimizing

where is the loss for task B only, indexes the weights and is an overall hyperparameter that balances the importance of the two tasks. EWC minimizes .

Why the name “elastic weight consolidation”? From the authors:

“This constraint is implemented as a quadratic penalty, and can therefore be imagined as a spring anchoring the parameters to the previous solution, hence the name elastic. Importantly, the stiffness of this spring should not be the same for all parameters; rather, it should be greater for those parameters that matter most to the performance during task A.”

References:

1. Kirkpartick, J., et al. (2016). Overcoming catastrophic forgetting in neural networks.

A simple probabilistic algorithm for estimating the number of distinct elements in a data stream

I just came across a really interesting and simple algorithm for estimating the number of distinct elements in a stream of data. The paper (Chakraborty et al. 2023) is available on arXiv; see this Quanta article (Reference 2) for a layman’s explanation.

Problem statement

Let’s state the problem formally. Let’s say we are given a stream where each , and let denote the number of distinct elements of . For any input parameters , we want an -approximation of , i.e. a number such that

A naive solution would be to maintain a buffer which keeps track of all distinct elements seen so far. As we inspect each new item, we check if the item is in : if it’s there we discard the new item, and if it’s not we add it to . After going through the items, the size of is exactly the number of distinct items. However, it has bad worst case memory and checking whether an element is in can be costly if is large (e.g. think of counting the number of distinct members visiting google.com today).

Algorithm

There is a long line of work seeking to do better than the naive solution (see Bibliographic Remarks in Reference 1); Chakraborty et al. claims to be the first state-of-the-art algorithm whose proof requires just basic probability theory. The algorithm is easy enough to be reproduced here; see Reference 2 for an excellent walkthrough of an example run of the algorithm.

Algorithm 1 from Chakraborty et al. (2023), annotations are mine.

At a high level, the algorithm maintains a buffer of seen items. The algorithm proceeds in rounds. In each round, we check whether each incoming item is in or not. If it is not in , add it with probability . If it is in , keep it in with probability . is made smaller and smaller as we go through more rounds. Whenever the buffer reaches its threshold number of items (threshold fixed throughout the algorithm), each element is removed from it with probability . When we run out of elements, our estimate of the number of distinct items is .

Here is the theorem for this algorithm:

Theorem. For any data stream and any , the algorithm above outputs an -approximation of . It uses space in the worst case.

The proof is not too long (about 3 pages) and while requiring some restatement of the problem via alternate algorithms, the most advanced tool it uses is Chernoff’s bound.

To obtain an -approximation for items, the algorithm requires a buffer of size . This grows logarithmically with . If we want our approximation to be twice as accurate, we need the buffer to be 4 times as large (). In practice the buffer can probably be much smaller: the authors are only interested in big-O type dependencies and so made the constants bigger wherever it made the proof simpler. (It’s worth noting that the algorithm gives the exact number of distinct items if the buffer happens to be larger than the true number of distinct items.)

Simulations

This algorithm is really easy to implement. Here is a quick implementation in R (definitely can be optimized further, just coded it up quickly for prototyping purposes):

# threshold: number of elements buffer can hold
# trajectory: if TRUE, return the estimate after every item in stream (default FALSE)
estimate_distinct_elements <- function(data_stream, threshold, trajectory = FALSE) {
  p <- 1
  m <- length(data_stream)
  buffer <- c()
  hit_fail_event <- FALSE
  estimates <- rep(NA, m)
  for (i in 1:m) {
    buffer <- setdiff(buffer, data_stream[i])
    if (runif(1) < p) buffer <- c(buffer, data_stream[i])
    if (length(buffer) == threshold) {
      buffer <- buffer[runif(threshold) < 0.5]
      p <- p / 2
      if (length(buffer) == threshold) {  # fail event
        hit_fail_event <- TRUE
      }
    }
    estimates[i] <- ifelse(hit_fail_event, -1, length(buffer) / p)
  }
  
  if (trajectory) {
    return(estimates)
  } else {
    return(estimates[m])
  }
}

Let’s take this for a run with some simulations. All simulation code can be found here. For a start, let’s generate a stream of length 100, where each item is an integer from 1 to 100 (inclusive). For the particular data stream we generated, there were 60 distinct elements. This means that for the naive solution, we would have buffer of length 60.

First, we run the algorithm 1000 times with buffer length 30, half of what the naive solution needed. The algorithm never hit the “fail” event, and here are some summary statistics of the distinct count estimate:

#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#  34.00   52.00   56.00   60.02   68.00  104.00 

The mean and median are pretty close to the true value! Here is a histogram of the estimates:

It’s worth noting that while the estimate is generally pretty good (50% chance of being within 8 of the true number), it can be wildly off as well!

Let’s rerun the simulation for the same data stream but with buffer length 15 (a quarter of the true number of distinct items). Again the algorithm never hit the “fail” event. The summary statistics and histogram are below and as one might expect, while the estimates are still centered at the right place, there is more variance in the estimate. (The histograms might look similar due to different scales on the x-axis.)

#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#   8.00   48.00   56.00   59.81   72.00  144.00 

Since this is a streaming algorithm, it provides an estimate of the number of distinct items after inspecting each item. Let’s look at how good the estimate is over time. Here are the plots for 10 estimate trajectories (true count in blue) for buffer length = 30 and buffer length = 15. Notice that the algorithm gives an exact estimate until the buffer is first filled. As expected, larger buffer length tends to give better estimates.

Let’s do one final simulation: Let’s look at a stream of length 10,000 with each element being a random integer from 1 to 4,000 (inclusive). For our data, the true number of distinct elements is 3,686.

Here is the plot of 10 trajectory estimates with buffer length 500 (about 14% of the true number of uniques) as well as the relative error of the estimate:

That’s not too bad! The estimates are pretty much within 10% of the true values. How about a much smaller buffer length of 50 (1.4% of the true number of uniques):

There’s much wider variation in the estimates now, probably too much to be useful for any real-world application.

References:

1. Chakraborty, S., et al. (2023). Distinct Elements in Streams: An Algorithm for the (Text) Book.
2. Nadis, S. (2024). Computer Scientists Invent an Efficient New Way to Count.

Vectorization and Kronecker products

If and , the Kronecker product is the block matrix defined by

For a matrix , the vectorization of is the column vector obtained by stacking ‘s columns on top of each other:

It turns out that vectorization and Kronecker products play nicely together as in the following proposition:

Proposition. If , and , then

Proof: This proof is taken from Reference 1. Let and denote the th row of and respectively. Let’s compute the th column of :

Stacking the columns of on top of each other:

which concludes the proof.

Here is one nice application of the proposition above. Consider the matrix equation . Vectorizing both sides (also known as the “vec trick”), we get

It can be shown that the Kronecker product is invertible if and only if and (or equivalently ) are invertible. Thus, it follows that has a unique solution if and only if and are invertible.

Here is a less trivial example (from commenter Jarle Tufto): consider the matrix equation . By applying the vec trick,

assuming is invertible. There are a few more examples of applying this trick to matrix equations in Reference 3.

References:

1. Kronecker Product and the vec Operator.
2. Wikipedia. Kronecker product.
3. Zhang, T. MATH 470 Independent Study in Matrix Theory: The Kronecker Product.

The mixed-product property of the Kronecker product

If and , the Kronecker product is the block matrix defined by

See the Wikipedia article (Reference 1) for a full list of properties of the Kronecker product. This post focuses on the mixed-product property of the Kronecker product. It’s a nice property that is often used to prove other properties of the Kronecker product. Here it is:

Proposition. For matrices , , , , we have

This is easy to prove if we use matrix multiplication for block matrices. We can write the RHS as

The LHS can be expanded as

Using block matrix multiplication, the th block of the LHS is

The mixed-product property can be used to prove many other properties. Here’s one corollary:

Corollary. For square matrices and and positive integer , we have

Proof: Proceed by induction on . The case of is trivial. Assume the statement is true for some . Replacing , , and in the mixed-product property with , , and , we get

which proves the induction step.

References:

1. Wikipedia. Kronecker product.

Rough approximations to move between logit and probability scales

The logit and expit functions are commonly used to move between the probability scale [0, 1] and the whole real number line, especially when working with supervised learning models with binary outcomes. These functions are very easy to define, here is an implementation in R:

logit <- function(p) {
  log(p / (1 - p))
}

expit <- function(x) {
  exp(x) / (1 + exp(x))
}

When working in this space, it’s helpful to have some approximations to move between the logit and probability scales. (As an analogy, it is helpful to know that for a normal distribution, the interval standard deviations around the mean covers about 95% of the distribution, while standard deviations covers about 99%.)

Here are some rough approximations to help you estimate probabilities from logit scores and vice versa:

Logit	Probability
-7	0.001
-5	0.01
-3	0.05
-2	0.1
-1	0.25
0	0.5
1	0.75
2	0.9
3	0.95
5	0.99
7	0.999

What is the Plackett-Luce model?

The Plackett-Luce model, attributed to Robert L. Plackett (Reference 1) and R. Duncan Luce (Reference 2), is a model for rankings data. In this post, we first describe the simpler Bradley-Terry model, then introduce the Plackett-Luce model as a generalization of the Bradley-Terry model.

The Bradley-Terry model

The Bradley-Terry model is a model that allows us to make pairwise comparisons between items. Given two items and , let denote “i is preferred to j”. The Bradley-Terry model assumes that pairwise preferences are generated according to the following probability model: there is some and associated with items and respectively, and

There are extensions to the Bradley-Terry model that allow for ties, home-field advantage and specific pairwise interaction terms; see for e.g. Davidson (1970) (Reference 3).

The Plackett-Luce model

Whereas the Bradley-Terry model works with pairwise preferences, the Plackett-Luce model works with listwise preferences. Let’s say we have some items , and we are given the ranked preference list . The Plackett-Luce model assumes that these preferences were generated according to a sequence of choices: the top-ranked item is chosen from the pool with some probability, then the second-ranked item is independently chosen from the remaining pool of items with some probability, and so on. Each item is associated with a real number . At each step, if denotes the remaining items, then the probability of choosing item from this pool of items is assumed to be

This is probably easier understood via an example. Assume that we have 3 items, and the preference we are given is . The probability associated with this ranked list is

If we only consider lists of length 2, the Plackett-Luce model becomes the Bradley-Terry model.

How to fit a Bradley-Terry model (and the Plackett-Luce model by extension)

For both the Bradley-Terry and Plackett-Luce models, the question is how to determine the ‘s for the items. This boils down to modeling the as a function of data/features associated with item . Assume that each item has a feature vector . We can write the model generically as

where for some parameter space . We then find the “optimal” based on some framework, most commonly maximum likelihood estimation (MLE).

Let’s write this out for the Bradley-Terry model for a specific model family (details with code in this previous blog post). Assume that is a linear combination of the features , and the only features we consider are the one-hot encoding features: “is this item ?” for all items . In other words, we assume

and our job is to estimate the . If our dataset consist of triplets of the form , then the likelihood function is

and the MLE is given by .

Practically, the PlackettLuce R package fits the Plackett-Luce model, and is able to handle generalizations such as ties and partial rankings (see the vignette Reference 4 for details).

References:

1. Plackett, R. L. (1975). “The Analysis of Permutations.”
2. Luce, R. D. (1959). Individual choice behavior.
3. Davidson, R. R. (1970). “On extending the bradley-terry model to accommodate ties in paired comparison experiments.”
4. Turner, H., et al. (2023). Introduction to PlackettLuce.