Extrapolating Expected Accuracies for Large Multi-Class Problems
Charles Zheng Rakesh Achanta Yuval Benjamini
stat.ML cs.CV cs.LG








Extrapolating Expected Accuracies for Large Multi-Class Problems
     Charles Zheng
 charles.y.zheng@gmail.com 

       Rakesh Achanta 
 rakesha@stanford.edu 

       Yuval Benjamini 
 yuval.benjamini@mail.huji.ac.il
    December 30, 2023
=======================================================================================================================================================




The difficulty of multi-class classification generally increases with
the number of classes.  Using data from a subset of the classes, can
we predict how well a classifier will scale with an increased number
of classes?  Under the assumptions that the classes are sampled
identically and independently from a population, and that the classifier is based on independently learned
scoring functions, we show that the expected accuracy when the
classifier is trained on k classes is the k-1st moment of a
certain distribution that can be estimated from data.  We present an
unbiased estimation method based on the theory, and demonstrate its
application on a facial recognition example.







§ INTRODUCTION


Many machine learning tasks are interested in recognizing or
identifying an individual instance within a large set of possible
candidates. These problems are usually modeled as multi-class
classification problems, with a large and possibly complex label
set. Leading examples include detecting the speaker from his voice
patterns <cit.>, identifying the author from her
written text <cit.>, or labeling the object
category from its image
<cit.>.  In all
these examples, the algorithm observes an input x, and uses the
classifier function h to guess the label y from a large label set
𝒮.

There are multiple practical challenges in developing classifiers for large label sets. Collecting high quality training data is
perhaps the main obstacle, as the costs scale with the number of
classes.  It can be affordable to first collect data for a small set
of classes, even if the long-term goal is to generalize to a larger
set.  Furthermore, classifier development can be accelerated by
training first on fewer classes, as each training cycle may require
substantially less resources.  Indeed, due to interest in how
small-set performance generalizes to larger sets, such comparisons can
found in the literature <cit.>.
A natural question is: how does changing the size of the label set
affect the classification accuracy?

We consider a pair of classification problems on finite
label sets: a source task with label set 𝒮_k_1 of size
k_1, and a target task with a larger label set 𝒮_k_2
of size k_2 > k_1.  For each label set 𝒮_k, one
constructs the classification rule h^(k):𝒳→𝒮_k.  Supposing that in each task, the test example
(X^*, Y^*) has a joint distribution, define the generalization
accuracy for label set 𝒮_k as

    GA_k = [h^(k)(X^*) = Y^*].

The problem of performance extrapolation is the following:
using data from only the source task 𝒮_k_1, 
predict the accuracy for a target task with a larger unobserved 
label set 𝒮_k_2.

A natural use case for performance extrapolation would be in the
deployment of a facial recognition system.  Suppose a system was
developed in the lab on a database of k_1 individuals. Clients would
like to deploy this system on a new larger set of k_2
individuals. Performance extrapolation could allow the lab to predict
how well the algorithm will perform on the client's problem,
accounting for the difference in label set size.

Extrapolation should be possible when the source and target
classifications belong to the same problem domain.  In
many cases, the set of categories 𝒮 is to some degree a
random or arbitrary selection out of a larger, perhaps infinite, set
of potential categories 𝒴. Yet any specific experiment
uses a fixed finite set.  For example, categories in the classical
Caltech-256 image recognition data set <cit.> were
assembled by aggregating keywords proposed by students and then
collecting matching images from the web.  The arbitrary nature of the
label set is even more apparent in biometric applications (face
recognition, authorship, fingerprint identification) where the labels
correspond to human individuals <cit.>.  In all these cases, the number of the
labels used to define a concrete data set is therefore an experimental
choice rather than a property of the domain.  Despite the arbitrary
nature of these choices, such data sets are viewed as representing the
larger problem of recognition within the given domain, in the sense
that success on such a data set should inform performance on similar
problems.

In this paper, we assume that both 𝒮_k_1 and
𝒮_k_2 are independent identically distributed (i.i.d.)
samples from a population (or prior distribution) of labels π, which is
defined on the label space 𝒴.  These assumptions help
concretely analyze the generalization accuracy, although both are only
approximate characterizations of the label selection process, which is
often at least partially manual. Since we assume the label set is
random, the generalization accuracy of a given classifier becomes a
random variable.  Performance extrapolation then becomes the problem
of estimating the average generalization accuracy AGA_k of an
i.i.d. label set 𝒮_k of size k.  The condition of
i.i.d. sampling of labels ensures that the separation of labels in a
random set 𝒮_k_2 can be inferred by looking at the
empirical separation in 𝒮_k_1, and therefore that some
estimate of the average accuracy on 𝒮_k_2 can be
obtained.  We also make the assumption that the classifiers train a
separate model for each class.  This convenient property allows us to
characterize the accuracy of the classifier by selectively
conditioning on one class at a time.

Our paper presents two main contributions related to extrapolation.
First, we present a theoretical formula describing how average accuracy
for smaller k is linked to average accuracy for label set of size
K>k.  We show that accuracy at any size depends on a
discriminability function D, which is determined by properties
of the data distribution and the classifier.  Second, we propose an
estimation procedure that allows extrapolation of the observed average
accuracy curve from k_1-class data to a larger number of classes,
based on the theoretical formula. Under certain conditions, the
estimation method has the property of being an unbiased estimator of
the average accuracy.

The paper is organized as follows.  In the rest of this section, we
discuss related work.  The framework of randomized classification is
introduced in Section <ref>, and there we also
introduce a toy example which is revisited throughout the
paper. Section <ref> develops our theory of
extrapolation, and Section <ref> we
suggest an estimation method. We evaluate our method using simulations in Section 4.
In Section
<ref>, we demonstrate our method on a facial
recognition problem, as well as an optical character recognition problem. In Section <ref> we
discuss modeling choices and limitations of our theory, as well as
potential extensions.




 §.§ Related Work


Linking performance between two different but related classification
tasks can be considered an instance of transfer learning
<cit.>. Under 's terminology,
our setup is an example of multi-task learning, because the source
task has labeled data, which is used to predict performance on a
target task that also has labeled data. Applied examples of transfer learning from one label set to another include <cit.>,
<cit.>, <cit.>. 
However, there is little theory for predicting the behavior of the learned classifier on a new label set. Instead, most research classification for large label sets deal with the computational challenges of jointly optimizing the many parameters
required for these models for specific classification algorithms <cit.>. <cit.>
presents a method for estimating the accuracy of a classifier which can
be used to improve performance for general classifiers, but doesn't apply for different set sizes.  

The theoretical framework we adopt is one where there exists a family
of classification problems with increasing number of classes. This
framework can be traced back to <cit.>, who considered the error rate of a random codebook, which is a special case of randomized classification. More recently, a number of authors have considered the problem of high-dimensional feature selection for multiclass
classification with a large number of classes <cit.>.  All of these works
assume specific distributional models for classification compared to
our more general setup. However, we do not deal with the problem of
feature selection.

Perhaps the most similar method that deals with extrapolation of classification error to a larger number of classes can be found in <cit.>. They trained a classifier for identifying the observed stimulus from a functional MRI scan of brain activity, and were interested in its performance on larger stimuli sets. They proposed an extrapolation algorithm as a heuristic with little theoretical discussion. In Section <ref> we interpret their method within our theory, and discuss cases where it performs well compared to our algorithm.





§ RANDOMIZED CLASSIFICATION


The randomized classification model we study has the following
features.  We assume that there exists an infinite, perhaps
continuous, label space 𝒴 and a example space 𝒳∈ℝ^p.  We assume there exists a prior distribution π
on the label space 𝒴.  And for each label y ∈𝒴, there exists a distribution of examples F_y. In other
words, for an example-label pair (X, Y), the conditional
distribution of X given Y = y is given by F_y.



A random classification task can be generated as follows.  The label
set 𝒮 = {Y^(1),, Y^(k)} is generated by
drawing labels Y^(1),, Y^(k) i.i.d. from π.  For each
label, we sample a training set and a test set.  The training set is
obtained by sampling r_train observations X_j, train^(i)
i.i.d. from F_Y^(i) for j = 1,, r_train and i =
1,, k.  The test set is likewise obtained by sampling r
observations X_j^(i) i.i.d. from F_Y^(i) for j = 1,,
r.

We assume that the classifier h(x) works by assigning a score to
each label y^(i)∈𝒮, then choosing the label with the
highest score.  That is, there exist real-valued score
  functions m_y^(i)(x) for each label y^(i)∈𝒮.  Since the classifier is allowed to depend on the
training data, it is convenient to view it (and its associated score
functions) as random.  We write H(x) when we wish to work with the
classifier as a random function, and likewise M_y(x) to denote the
score functions whenever they are considered as random.

For a fixed instance of the classification task with labels
𝒮 = {y^(i)}_i=1^k and associated score functions
{m_y^(i)}_i=1^k, recall the definition of the k-class
generalization error (<ref>).  Assuming that there are no
ties, it can be written in terms of score functions as

    GA_k(h) = 1/k∑_i=1^k  [m_y^(i)(X^(i)) = max_j
    m_y^(j)(X^(i))],

where X^(i)∼ F_y^(i) for i =1,, k.  However, when
we consider the labels {Y^(i)}_i=1^k and associated score
functions to be random, the generalization accuracy also becomes a
random variable.

Suppose we specify k but do not fix any of the random quantities in
the classification task.  Then the k-class average
  generalization accuracy of a classifier is the expected value of
the generalization accuracy GA_k(H) resulting from a random
set of k labels, Y^(1), , Y^(k)iid∼π, and their associated score functions,

    AGA_k    = 1/k∑_i=1^k [M_Y^(i)(X^(i)) = max_j
    M_Y^(j)(X^(i))]
    
       = [M_Y^(1)(X^(1)) = max_j M_Y^(j)(X^(1))].

The last line follows from noting that all k summands in the
previous line are identical.  The definition of average
  generalization accuracy is illustrated in Figure
  <ref>.





 §.§ Marginal Classifier


In our analysis, we do not want the classifier to rely too strongly on
complicated interactions between the labels in the set. We therefore
propose the following property of marginal separability for
classification models:


The classifier H(x) is called a marginal classifier if the
score function M_y^(i)(x) only depends on the label y^(i)
and the class training set X_j, train^(i); that is, for some function g,

    M_y^(i)(x) = g(x; y^(i),X_1, train^(i),...,X_r_train, train^(i)).


This means that the score function for y^(i) does not depend on
other labels y^(j) or their training samples.  Therefore, each
M_y can be considered to have been drawn from a distribution
ν_y.  Classes “compete” only through selecting the highest
score, but not in constructing the score functions.  The operation of
a marginal classifier is illustrated in Figure
<ref>.




The marginal property allows us to prove
strong results about the accuracy of the classifier under
i.i.d. sampling assumptions.


Comments:


  * If H is a marginal classifier then 
M_Y^(i) is independent of Y^(j) and M_Y^(j) for i ≠ j.

  * Estimated Bayes classifiers are primary examples of marginal
  classifiers. Let f̂_̂ŷ be a density estimate of the example
  distribution under label y obtained from the empirical
  distribution F̂_̂ŷ. Then, we can use the estimated density to
  produce the score functions:

    M^EB_y(x) = log(f̂_̂ŷ(x)).

The resulting empirical approximation for the Bayes classifier would
be

    H^EB(x) = argmax_Y ∈𝒮(M^EB_Y(x)).


  * Both Quadratic Discriminant Analysis (QDA) and na ive Bayes
  classifiers can be seen as specific instances of an estimated Bayes
  classifier.
[QDA is the special case of the estimated Bayes classifier
  when f̂_̂ŷ is obtained as the multivariate Gaussian density
  with mean and covariance parameters estimated from the data.  Naive
  Bayes is the estimated Bayes classifier when f̂_̂ŷ is obtained
  as the product of estimated componentwise marginal distributions of
  p(x_i|y).]  For QDA, the score function is given by

    m_y^QDA(x) = -(x - μ(F̂_y))^T Σ(F̂_y)^-1 (x-μ(F̂_y)) - log(Σ(F̂_y)),

where μ(F) = ∫ y dF(y) and Σ(F) = ∫ (y-μ(F))(y-μ(F))^T dF(y).
In Naive Bayes, the score function is

    m^NB_y(x) = ∑_j=1^p logf̂_y, j(x),

where f̂_y, j is a density estimate for the j-th component of
F̂_y.

  * For some classifiers, M_y is a deterministic function of y
  (and therefore ν_y is degenerate). A prime example is when
  there exist fixed or pre-trained embeddings g, g̃ that map
  labels y and examples x into R^p. Then

    M_y^embed = -g(y) - g̃(x)_2.


  * There are many classifiers which do not satisfy the marginal
  property, such as multinomial logistic regression, multilayer neural
  networks, decision trees, and k-nearest neighbors.


Notational remark.  Henceforth, we shall relax the assumption
that the classifier H(x) is based on a training set.  Instead, we
assume that there exist score functions {M_Y^(i)}_i=1^k
associated with the random label set {Y^(i)}_i=1^k, and that
the score functions M_Y^(i) are independent of the test set.
The classifier H(x) is marginal if and only if M_Y^(i) are
independent of both Y^(j) and M_Y^(j) for j ≠ i.



 §.§ Estimation of Average Accuracy

Before tackling extrapolation, it is useful to discuss a simpler task of generalizing accuracy results when the target set is not larger than the source set. Suppose we have test data for a classification task with k_1
classes.  That is, we have a label set 𝒮_k_1 =
{y^(i)}_i=1^k_1 and its associated set of score functions
M_y^(i), as well as test observations (x_1^(i),,
x_r^(i)) for i = 1,, k_1.  What would be the predicted
accuracy for a new randomly sampled set of k_2 ≤ k_1 labels?

Note that AGA_k_2 is the expected value of the accuracy on
the new set of k_2 labels.  Therefore, any unbiased estimator of
AGA_k_2 will be an unbiased predictor for the accuracy on
the new set.

Let us start with the case k_2 = k_1 = k.  For each test observation
x_j^(i), define the ranks of the candidate classes ℓ =
1,, k by

    R_j^i, ℓ = ∑_s = 1^k I{m_y^(ℓ)(x_j^(i)) ≥ m_y^(s)(x_j^(i))}.

The test accuracy is the fraction of observations for which the
correct class also has the highest rank

    TA_k = 1/r k∑_i=1^k∑_j=1^r I{R_j^i,i = k}.

Taking expectations over both the test set and the random labels, the
expected value of the test accuracy is AGA_k.  Therefore, in this special case, TA_k provides an unbiased estimator for AGA_k_2.

Next, let us consider the case where k_2 < k_1.  Consider label set
𝒮_k_2 obtained by sampling k_2 labels uniformly
without replacement from 𝒮_k_1. Since
𝒮_k_2 is unconditionally an i.i.d. sample from the
population of labels π, the test accuracy of 𝒮_k_2
is an unbiased estimator of AGA_k_2.  However, we can get a
better unbiased estimate of AGA_k_2 by averaging over all
the possible subsamples 𝒮_k_2⊂𝒮_k_1.
This defines the average test accuracy over subsampled tasks,
ATA_k_2.

Remark.  Naïvely, computing ATA_k_2 requires us
to train and evaluate k_1k_2 classification rules.
However, for marginal classifiers, retraining the classifier is not
necessary.  Looking at the rank R_j^i,i of the correct label i
for x_j^(i), allows us to determine how many subsets
𝒮_2 will result in a correct classification. Specifically,
there are R_j^i,i - 1 labels with a lower score than the correct
label i.  Therefore, as long as one of the classes in
𝒮_2 is i, and the other k_2-1 labels are from the set
of R_j^i,i-1 labels with lower score than i, the
classification of x_j^(i) will be correct.  This implies that
there are R_j^i,i-1k_2-1 such subsets 𝒮_2
where x_j^(i) is classified correctly, and therefore the average
test accuracy for all k_1k_2 subsets 𝒮_2 is

    ATA_k_2 = 1/k_1k_21/r k_2∑_i=1^k_1∑_j=1^rR_j^i,i-1k_2-1.






 §.§ Toy Example: Bivariate Normal



Let us illustrate these ideas using a toy example.  Let (Y, X) have
a bivariate normal joint distribution,

    (Y, X) ∼ N([ 0; 0 ], [ 1 ρ; ρ 1 ]),

as illustrated in Figure <ref>(a).  Therefore, for a given
randomly drawn label Y, the conditional
distribution of X for that label is univariate normal with mean ρ Y and variance 1-ρ^2,

    X|Y = y ∼ N(ρ Y, 1-ρ^2).

Supposing we draw k = 3 labels {y_1,y_2, y_3}, the classification
problem will be to assign a test instance X^* to the correct label.
The test instance X^* would be drawn with equal probability from one
of three conditional distributions X | Y=y^(i), as illustrated in
Figure <ref>(b, top).  The Bayes rule assigns X^* to the
class with the highest density p(x|y_i), as illustrated by Figure
<ref>(b, bottom): it is therefore a marginal classifier, with
score function

    M_y^(i)(x) = log(p(x|y^(i))) = -(x - ρ y)^2/2(1-ρ^2)  + const.







§ EXTRAPOLATION




For this model, the generalization accuracy of the Bayes rule for any
label set {y^(1),, y^(k)} is given by

    GA_k(y_1,, y_k)    = 1/k∑_i=1^k _X ∼ p(x|y_i)[p(X|y_i) = max_j=1^k p(X|y_j)]
    
       = 1/k∑_i=1^k Φ(y^[i+1] - y^[i]/2√(1-ρ^2)) - Φ(y^[i-1] - y^[i]/2√(1-ρ^2)),

where Φ is the standard normal cdf, y^[1] < ⋯ < y^[k]
are the sorted labels, y^[0] = -∞ and y^[k+1] =
∞.  We numerically computed GA_k(Y_1,, Y_k) for
randomly drawn labels Y_1,, Y_k iid∼ N(0, 1), and 
the distributions of GA_k for k = 2,, 10 are
illustrated in Figure <ref>.  The mean of the distribution of
GA_k is the k-class average accuracy, AGA_k. The
theory presented in the next section deals with how to analyze the
average accuracy AGA_k as a function of k.




The section is organized as follows.  We begin by introducing an
explicit formula for the average accuracy AGA_k.  The
formula reveals that AGA_k is determined by moments of a
one-dimensional function D(u).  Using this formula, we can estimate 
D(u) using subsampled accuracies.  These estimates allow us to extrapolate the average
generalization accuracy to an arbitrary number of labels.

The result of our analysis is to expose the average accuracy
AGA_k as the weighted average of a function D(u),
where D(u) is independent of k, and where k only changes
the weighting.  The result is stated as follows.


Suppose π, {F_y}_y ∈𝒴, and score functions M_y
satisfy the tie-breaking condition.  Then, there exists a cumulative
distribution function D(u) defined on the interval [0,1]
such that

    AGA_k = 1 - (k-1) ∫D(u) u^k-2 du.



The tie-breaking allows us
to neglect specifying the case when
margins are tied.

Tie-breaking condition: for all x ∈𝒳,
M_Y(x) ≠ M_Y'(x)
with probability one for Y, Y' independently drawn from π.

In practice, one can simply break ties randomly,
which is mathematically equivalent to adding a small amount of random
noise ϵ to the function ℳ.



 §.§ Analysis of Average Accuracy


For the following discussion, we often consider a random label with
its associated score function and example vector. Explicitly, this
sampling can be written:

    Y ∼π,  M_Y|Y ∼ν_Y,  X|Y ∼ F_Y.



Similarly we use (Y',M_Y',X') and (Y^*,M_Y^*,X^*) for two more
triplets with independent and identical distributions. Specifically,
X^* will typically note the test example, and therefore Y^* the
true label and M_Y^* its score function.

The function D is related to a favorability
function. Favorability measures the probability that the score for the
example x^* is going to be maximized by a particular score function m_y,
compared to a random competitor M_Y'.  Formally, we write

    U_x^*(m_y) = [m_y(x^*) > M_Y'(x^*)].


Note that for fixed example x^*, favorability is monotonically
increasing in m_y(x^*).  If m_y(x^*) > m_y^†(x^*), then
U_x^*(y) > U_x^*(y^†), because the event {m_y(x^*) >
M_Y'(x^*)} contains the event {m_y^†(x^*) >
M_Y'(x^*)}.

Therefore, given labels y^(1),,y^(k) and test instance
x^*, we can think of the classifier as choosing the label with the
greatest favorability:

    ŷ = _y^(i)∈𝒮 m_y^(i)(x^*) = _y^(i)∈𝒮 U_x^*(m_y^(i)).

Furthermore, via a conditioning argument, we see that this is still
the case even when the test instance and labels are random:

    Ŷ = _Y^(i)∈𝒮 M_Y^(i)(X^*) = _Y^(i)∈𝒮 U_X^*(M_Y^(i)).



The favorability takes values between 0 and 1, and when any of its
arguments are random, it becomes a random variable with a distribution
supported on [0,1].  In particular, we consider the following two
random variables:


  a. the incorrect-label favorability U_x^*(M_Y)
  between a given fixed test instance x^*, and the score function of
  a random incorrect label M_Y, and

  b. the correct-label favorability U_X^*(M_Y^*)
  between a random test instance X^*, and the score function of the
  correct label, M_Y^*.



  §.§.§ Incorrect-Label Favorability

The incorrect-label favorability can be written explicitly as

    U_x^*(M_Y) = [M_Y(x^*) > M_Y'(x^*)|M_Y].

Note that M_Y and M_Y' are identically distributed, and are both
are unrelated to x^* that is fixed. This leads to the following
result:

Under the tie-breaking condition, the incorrect-label favorability
U_x^*(M_Y) is uniformly distributed for any x^* ∈𝒳,
meaning 
    [U_x^*(M_Y) ≤ u] = u

for all u ∈ [0,1].


 Write U_x^*(M_Y) = [Z > Z'|Z], where Z = M_Y(x) and Z' =
M_Y'(x) for Y, Y' i.i.d.∼π.
The tie-breaking condition implies that [Z=Z']=0.  Now observe that for independent random variables Z, Z' with Z D= Z' and [Z=Z']=0, the conditional probability
[Z > Z'|Z] is uniformly distributed.






  §.§.§ Correct-Label Favorability


The correct-label favorability is 

    U^* = U_X^*(M_Y^*) = [M_Y^*(X^*) > M_Y'(X^*)|Y^*,M_Y^*,X^*].

The distribution of U^* will depend on π, {F_y}_y ∈𝒮 and {ν_y}_y ∈𝒮, and
generally cannot be written in a closed form.  However, this
distribution is central to our analysis–indeed, we will see that the
function D appearing in theorem <ref>
is defined as the cumulative distribution function of U^*.

The special case of k=2 shows the relation between the distribution
of U^* and the average generalization accuracy, AGA_2. In
the two-class case, the average generalization accuracy is the
probability that a random correct label score function gives a larger
value than a random distractor:

    AGA_2 = [M_Y^*(X^*) > M_Y'(X^*)].

where Y^* is the correct label, and Y' is a random incorrect
label.  If we condition on Y^*, M_Y^* and X^*, we get

    AGA_2 = [[M_Y^*(X^*) > M_Y'(X^*)|Y^*, M_Y^*, X^*]].

Here, the conditional probability inside the expectation is the
correct-label favorability.  Therefore,

    AGA_2 = [U^*] = ∫D(u) du,

where D(u) is the cumulative distribution function of U^*,
D(u) = [U^* ≤ u].  Theorem <ref>
extends this to general k; we now give the proof.




 Without loss of generality, suppose that the true label is Y^* and
the incorrect labels are Y^(1),, Y^(k-1).  We have

    AGA_k = [M_Y^*(X^*) > max_i=1^k-1 M_Y^(i)(X^*)]
    = [U^* > max_i=1^k-1 U_X^*(M_Y^(i))],

recalling that U^* = U_X^*(M_Y^*).  Now, if we condition on X^*
= x^*, Y^* = y^* and M_Y^* = m_y^*, then the random variable
U^* becomes fixed, with value

    u^* = U_x^*(m_y^*).

Therefore,

    AGA_k    =[[U^* > max_i=1^k-1 U_X^*(M_Y^(i))|X^* = x^*, Y^* = y^*, M_Y^* = m_y^*]]
    
       = [[U^* > max_i=1^k-1 U_X^*(M_Y^(i))|X^* = x^*, U^* = u^*]].

Now define U_max, k-1 = max_i=1^k-1 U_X^*(M_Y^(i)). 
Since by Lemma <ref>,
U_X^*(M_Y^(i)) are i.i.d. uniform conditional on X^* = x^*, we know that

    U_max, k-1|X^* = x^* ∼Beta(k-1, 1).

Furthermore, U_max, k-1 is independent of U^* conditional on
X^*.  Therefore, the conditional probability can be computed as

    [U^* > U_max, k-1|X^* = x^*, U^* = u^*] = ∫_u^*^1 (k-1) u^k-2 du.

Consequently,

    AGA_k    = [[U^* > max_i=1^k-1 U_X^*(M_Y^(i))|X^* = x^*, U^* = u^*]]
    
       = [∫_0^U^* (k-1) u^k-2 du|U^* = u^*]
    
       = [∫_0^1 I{u ≤ U^*} (k-1) u^k-2 du ]
    
       = (k-1) ∫_0^1 [U^* ≥ u] u^k-2 du
    
       = 1 - (k-1) ∫_0^1 [U^* ≤ u] u^k-2 du.

By defining D(u) as the cumulative distribution function of
U^* on [0,1],

    D(u) = [U_X^*(M_Y^*) ≤ u],

and substituting this definition into (<ref>), we obtain the identity (<ref>).



Theorem <ref> expresses the average accuracy
as a weighted integral of the function D(u).  Essentially, this theoretical result allows us
to reduce the problem of estimating AGA_k to one of estimating D(u).
But how shall we estimate D(u) from data?
We propose using non-parametric regression for this purpose in Section <ref>.




 §.§ Favorability and Average Accuracy for the Toy Example


Recall that for the toy example from Section <ref>, the
score function M_y was a non-random function of y that measures
the distance between x and ρ y

    M_y(x^*) = log(p(x^*|y)) = -(x^* - ρ y)^2/2(1-ρ^2) .


For this model, the favorability function U_x^*(m_y) compares the
distance between x^* and ρ y to the distance between x^* and
ρ Y' for a randomly chosen distractor Y'∼ N(0,1):

    U_x^*(m_y)    = [|ρ y - x^*|> |ρ Y' - x^*|] 
    
       = Φ(x^* + |ρ y - x^*|/ρ) - Φ(x^* - |ρ y - x^*|/ρ),

where Φ is the standard normal cumulative distribution function.
Figure <ref>(a) illustrates the level sets of the function
U_x^*(m_y).  The highest values of U_x^*(m_y) are near the
line x^* = ρ y corresponding to the conditional mean of X|Y, and as
one moves farther from the line, U_x^*(m_y) decays.  Note, however,
that large values of x^* and y (with the same sign) result in
larger values of U_x^*(m_y) since it becomes unlikely for Y' ∼
N(0,1) to exceed Y = y.

Using the formula above, we can calculate the correct-label
favorability U^* = U_X^*(M_Y^*) and its cumulative distribution
function D(u).  The function D is illustrated in Figure
<ref>(b) for the current example with ρ = 0.7.  The red
curve in Figure <ref> was computed using the formula

    AGA_k = 1-(k-1) ∫D(u) u^k-2 du.


It is illuminating to consider how the average accuracy curves and the
D(u) functions vary as we change the parameter ρ.  Higher
correlations ρ lead to higher accuracy, as seen in Figure
<ref>(a), where the accuracy curves are shifted upward as
ρ increases from 0.3 to 0.9.  The favorability U_x^*(m_y)
tends to be higher on average as well, which leads to lower values of
the cumulative distribution function–as we see in Figure
<ref>(b), where the function D(u) becomes smaller as
ρ increases.







 §.§ Estimation


Next, we discuss how to use data from smaller classification tasks to
extrapolate average accuracy.  Assume that we have data from a
k_1-class random classification task, and would like to estimate the
average accuracy AGA_k_2 for k_2>k_1 classes.  
Our estimation method will use the k-class average test accuracies,
ATA_2,...,ATA_k_1 (see Eq <ref>), for
its inputs.

The key to understanding the behavior of the average accuracy
AGA_k is the function D.  We adopt a linear model

    D(u) = ∑_ℓ = 1^m β_ℓ h_ℓ(u),

where h_ℓ(u) are known basis functions, and β_ℓ are the
linear coefficients to be estimated.  Since our proposed method is
based on the linearity assumption (<ref>), we refer to it as ClassExReg,
meaning Classification Extrapolation using Regression.

Conveniently, AGA_k can also be expressed in terms of the β_ℓ coefficients.
If we plug in the assumed linear model (<ref>) into the
identity (<ref>), then we get

    1 - AGA_k    = (k-1)∫D(u) u^k-2 du
    
       = (k-1)∫_0^1 ∑_ℓ = 1^m β_ℓ h_ℓ(u) u^k-2 du
    
       = ∑_ℓ = 1^m β_ℓ H_ℓ,k,

where

    H_ℓ,k = (k-1) ∫_0^1 h_ℓ(u) u^k-2 du.

The constants H_ℓ, k are moments of the basis function
h_ℓ.  Note that H_ℓ, k can be precomputed numerically for any k ≥ 2.

Now, since the test accuracies ATA_k are unbiased estimates
of AGA_k, this implies that the regression estimate

    β̂ = _β∑_k=2^k_1( (1 - ATA_k) - ∑_ℓ=1^m β_ℓ H_ℓ, k)^2,

is unbiased for β. The estimate of AGA_k_2 is similarly obtained
from (<ref>), via

    AGA_k_2 = 1 - ∑_ℓ=1^m β̂_ℓ H_ℓ, k_2.




 §.§ Model Selection


Accurate extrapolation using ClassExReg depends on a good
fit between the linear model (<ref>) and the true
discriminability function D(u).  However, since the function D(u)
depends on the unknown joint distribution of the data, it makes sense
to let the data help us choose a good basis {h_u} from a set of
candidate bases.

Let B_1,, B_s be a set of candidate bases, with B_i = {h_u^(i)}_u=1^m_i.  Ideally, we would like our model selection procedure to choose the B_i that obtains the best root-mean-squared error (RMSE) on the extrapolation from k_1 to k_2 classes.  As an approximation, we estimate the RMSE of extrapolation from k_1/2 source classes to k_1 target classes, by means of the “bootstrap principle.”  This amounts to a resampling-based model selection approach, where we perform extrapolations from k_0 = ⌊k_1/2⌋ classes to k_1 classes, and evaluate methods based on how closely the predicted ABA_k_1 matches the test accuracy ATA_k_1.  To elaborate, our model selection procedure is as follows.



  * For ℓ=1,,L resampling steps:


  * Subsample 𝒮_k_0^(ℓ) from 𝒮_k_1 uniformly with replacement.

  * Compute average test accuracies ATA_2^(ℓ),,ATA_k_0^(ℓ) from the subsample 𝒮_k_0^(ℓ).

  * For each candidate basis B_i, with i = 1,, s:


  * Compute β̂^(i,ℓ) by solving the least-squares problem

    β̂^(i,ℓ) = _β∑_k=2^k_0( (1 - ATA_k^(ℓ)) - ∑_j=1^m_iβ_j H_j, k^(i))^2.


  * Estimate AGA_k_1^(i,ℓ) by

    AGA_k_1^(i,ℓ) = ∑_ℓ=1^m_iβ̂_j^(i,ℓ) H_j, k_1^(i).




  * Select the basis B_i^* by

    i^* = argmin_i=1^s ∑_ℓ=1^L (AGA_k_1^(i,ℓ) - ATA_k_1)^2.


  * Use the basis B_i^* to extrapolate from k_1 classes (the full data) to k_2 classes.




§ SIMULATION STUDY


We ran simulations to check how the proposed extrapolation method, ClassExReg,
performs in different settings.  The results are displayed in Figure
<ref>. 
We varied the number of classes k_1 in the
source data set, the difficulty of classification, and the basis
functions. We generated data according to a mixture of isotropic
multivariate Gaussian distributions: labels Y were sampled from Y
∼ N(0, I_10), and the examples for each label sampled from X|Y
∼ N(Y, σ^2 I_10). The noise-level parameter σ
determines the difficulty of classification. Similarly to the
real-data example, we consider a 1-nearest neighbor classifier, which
is given a single training instance per class.




For the estimation, we use the model selection procedure described
in section <ref> to select the parameter h of the
“radial basis”

    h_ℓ(u) = Φ(Φ^-1(u) - t_ℓ/h).

where t_ℓ are a set of regularly spaced knots which are
determined by h and the problem parameters. Additionally, we add a
constant element to the basis, equivalent to adding an intercept to
the linear model (<ref>). 

The rationale behind the radial basis is to model the density of
Φ^-1(U^*) as a mixture of gaussian kernels with variance h^2.
To control overfitting, the knots are separated by at least a distance
of h/2, and the largest knots have absolute value Φ^-1(1 -
1/rk_1^2). The size of the maximum knot is set this way since
rk_1^2 is the number of ranks that are calculated and used by our
method.  Therefore, we do not expect the training data to contain
enough information to allow our method to distinguish between more
than rk_1^2 possible accuracies, and hence we set the maximum knot
to prevent the inclusion of a basis element that has on average a
higher mean value than u = 1-1/rk_1^2.  However, in
simulations we find that the performance of the basis depends only
weakly on the exact positioning and maximum size of the knots, as long
as sufficiently large knots are included. As is the case throughout non-parametric statistics, the bandwidth h is the
most crucial parameter.  In the simulation, we use a grid h = {0.1,
0.2, , 1} for bandwidth selection.





 §.§ Comparison to Kay


In their paper,[The KDE extrapolation method is described in page 29 of supplement to <cit.>.  While the method is only described for
a one-nearest neighbor classifier and for the setting where there is at most
one test observation per class, we have taken the liberty of extending it to a generic multi-class classification problem.]
<cit.> proposed a method for extrapolating classification
accuracy to a larger number of classes. The method depends on repeated
kernel-density estimation (KDE) steps. Because the method is only
briefly motivated in the original text, we present it in our
notation. 

For k_1 observed classes, let m_y^(j)(x^(i)_ℓ)  1≤ i,j≤ k_1 be
the observed score comparing feature vector of the ℓ'th test example of the i'th class x^(i)_ℓ
to the model trained for the j'th class m_y^(j). 
For each feature-vector x^(i)_ℓ, the density of wrong-class scores is estimated by smoothing the observed scores
with a kernel function K(·,·) with bandwidth h,

    f̂_ℓ^(i)(m) = 1/k_1-1∑_j≠ i K_h(m_y^(j)(x^(i)_ℓ), m).

An estimate of favorability for x^(i)_ℓ against a
single competitor is obtained by integrating the density below the observed true score m_y^(i)(x^(i)_ℓ), 

    acc_2(x_ℓ^(i)) = ∫_0^m_y^(i)(x^(i)_ℓ)f̂_ℓ^(i)(m)dm,

and the probability of accurate classification against K-1 random
competitors is acc_K(x^(i)_ℓ) = (acc_2(x^(i)_ℓ))^K-1. The average of these probabilities is the estimate for generalization accuracy

    ÂĜÂ_K^(KDE) = 1/k_1 r∑_i=1^k_1∑_ℓ=1^r acc_K(x^(i)_ℓ).


Note that the KDE method depends non-trivially on the smoothing bandwidth used in the density estimation step: when the kernel bandwidth is too small compared to the number of classes, the extrapolated accuracy will be upwardly biased.  To see
this, consider a feature vector x^(i)_ℓ that is correctly classified in the original class set. That is, m_y^(i)(x^(i)_ℓ)>m_y^(j)(x^(i)_ℓ) for every
j≠ i. If we use for f̂_i the unsmoothed empirical density,
acc_2(x^(i)_ℓ)=1 and therefore acc_K(x^(i)_ℓ) = 1 as well. The method relies on smoothing of each class to generate a the tail density that exceeds m_y^(i)(x^(i)_ℓ), and therefore it is highly dependent on the choice of kernel bandwidth. 

<cit.> recommended using a Gaussian
kernel, with the bandwidth chosen via pseudolikelihood cross-validation <cit.>.
In our simulation, we tested our own implementation of the KDE method,
using the two methods for cross-validated KDE estimation provided in
the stats package in the R statistical computing
environment: biased cross-validation and unbiased cross-validation <cit.>.




 §.§ Simulation Results


We see in Figure <ref> that ClassExReg and the KDE methods with unbiased and biased
cross-validation (KDE-UCV, KDE-BCV) perform comparably in the Gaussian
simulations.  We studied how the difficulty of extrapolation relates to both the
absolute size of the number of classes and the extrapolation factor k_2/k_1.  Our simulation has two settings for k_1 = {500,5000}, and
within each setting we have extrapolations to 2 times, 4 times, 10
times, and 20 times the number of classes.



Within each problem setting defined by the number of
source and target classes (k_1,k_2), we use the maximum RMSE across
all signal-to-noise settings to quantify the overall performance of
the method, as displayed in Table <ref>.

The results also indicate that more accurate extrapolation
appears to be possible for smaller extrapolation ratios
k_2/k_1 and larger k_1.
ClassExReg improves in worst-case RMSE when moving from k_1 =500 to k_1 = 5000 while keeping the extrapolation factor fixed, most dramatically in the case k_2/k_1 = 2 when it improves from a
maximum RMSE of 0.032±0.001 (k_1 = 500) to 0.009 ± 0.000 (k_1 = 5000), which is 3.5-fold reduction in worst-case RMSE, but
also benefiting from at least a 1.8-fold reduction in RMSE when going from
the smaller problem to the larger problem in the other three cases.

The kernel-density method produces comparable results, but is seen to depend strongly on the choice of bandwidth selection: KDE-UCV and KDE-BCV show very different performance profiles, although they differ only in the method used to choose the bandwidth.
Also, the KDE methods show significant estimation bias, as can be seen from Figure <ref>.  
This can be explained by the fact that
the KDE method ignores the bias introduced by exponentiation.
That is, even if p̂ is an unbiased estimator of p, p̂^k may not be a very good estimate of p^k, since

    p^k = [p̂]^k ≠[p̂^k].

unless p̂ is a degenerate random variable
(constant). Otherwise, for large k, p̂^k may be an extremely
biased estimate of p.  ClassExReg avoids this source of
bias by estimating the (k-1)st moment of D(u) directly.  As we see
in Figure <ref>, correcting for the bias of
exponentiation helps greatly to reduce the overall bias.   Indeed, while
ClassExReg shows comparable bias for the 500 to 10000 extrapolation,
the bias is very well-controlled in all of the k_1 = 5000 extrapolations.
































§ EXPERIMENTAL EVALUATION


We demonstrate the extrapolation of average accuracy in two data examples:
(i) predicting the
accuracy of a face recognition on a large set of labels from the
system's accuracy on a smaller subset, and (ii) extrapolating the performance of various classifiers on an optical character recognition (OCR) problem in the Telugu script, which has over 400 glyphs.

The face-recognition example takes data from the “Labeled Faces in the Wild” data set (<cit.>), where we selected the 1672 individuals with at least 2 face photos.  We form a
data set consisting of photo-label pairs (x_j^(i), y^(i))
for i = 1,, 1672 and j = 1,2 by randomly selecting 2 face
photos for each individual. 
We used the OpenFace (<cit.>) embedding for feature
extraction.[For each photo x, a 128-dimensional feature vector
g(x) is obtained as follows.  The computer vision library DLLib is
used to detect landmarks in x, and to apply a nonlinear
transformation to align x to a template.  The aligned photograph is
then downsampled to a 96 × 96 image. The downsampled image is
fed into a pre-trained deep convolutional neural network to obtain the
128-dimensional feature vector g(x). More details are found in
<cit.>.]
In order to identify a new photo x^*, we obtain the feature
vector g(x^*) from the OpenFace network, and guess the label ŷ
with the minimal Euclidean distance between g(y^(i)) and g(x^*),
which implies a score function

    M_y^(i)(x^*) = -||g(x_1^(i)) - g(x^*)||^2.

In this way, we can compute the test accuracy on all 1672 classes, TA_1672, but we also subsample k_1 = {100,200,400} classes in order to extrapolate from k_1 to 1672 classes.







In the Telugu optical character recognition example
(<cit.>), we consider
the use of three different classifiers: logistic
regression, linear support-vector machine (SVM), and a
deep convolutional neural network.[The network architecture is
  as follows: 48x48-4C3-MP2-6C3-8C3-MP2-32C3-50C3-MP2-200C3-SM.] 
The full data consists of 400 classes with 50 training and 50 test observations for each class. We create a nested hierarchy of subsampled data sets consisting of (i) a subset of 100 classes uniformly sampled without replacement from the 400 classes, and (ii) a subset consisting of 20 classes uniformly sampled without replacement from the size-100 subsample.   We therefore study three different prediction extrapolation problems:


  * Predicting the accuracy on k_2 = 100 classes from k_1 = 20 classes, comparing the predicted accuracy to the test accuracy of the classifier on the 100-class subsample as ground truth.

  * Same as (1), but setting k_2 = 400 and k_1 = 20, and using the full data set for the ground truth.

  * Same as (2), but setting k_2 = 400 and k_1 = 100.

Note that unlike in the case of the face recognition example,
here the assumption of marginal classification is satisfied for none of
the classifiers.  
We compare the result of our model to the ground
truth obtained by using the full data set.



 §.§ Results


The extrapolation results for the face recognition problem can be seen in Figure
<ref>, which plots the extrapolated accuracy
curves for each method for 100 different subsamples of size k_1.  As
can be seen, for all three methods, the variances decrease rapidly as
k_1 increases. 

The root-mean-square errors between at k_2=1672 can be seen in Table
<ref>.  KDE-BCV achieves the best extrapolation for
all three cases k_1= {100,200,400} with KDE-UCV consistently
achieving second place.  These results differ from the ranking of the
RMSEs for the analagous simulation when predicting k_2 = 2000 from
k_1 = 500 for accuracies around 0.45: in the first row and second
column of Figure <ref>, where the true accuracy is 0.43
(from setting σ^2=0.2), the lowest RMSE belongs to KDE-UCV
(RMSE=0.0361 ± 0.001), followed closely by ClassExReg (RMSE=0.0372
± 0.002), and KDE-BCV (RMSE=0.0635±0.001) having the highest
RMSE.  These discrepancies could be explained by differences between
the data distributions between the simulation and the face recognition
example, and also by the fact that we only have access to the k_2 =
1672-class ground truth for the real data example.





The results for Telugu OCR classification are displayed in Table <ref>.  If we rank the three extrapolation methods in terms of distance to the ground truth accuracy, we see a consistent pattern of rankings between the 20-to-100 extrapolation and the 100-to-400 extrapolation.  As we remarked in the simulation, the difficulty of extrapolation appears to be primarily sensitive to the extrapolation ratio k_2/k_1, which are similar (5 versus 4) in the 20-to-100 and 100-to-400 problems.  In both settings, ClassExReg comes closest to the ground truth for the Deep CNN and the SVM, but KDE-BCV comes closest to ground truth for the Logistic regression.  However, even for logistic regression, ClassExReg does better or comparably to KDE-UCV.


In the 20-to-400 extrapolation, which has the highest extrapolation ratio (k_2/k_1 = 20), none of the three extrapolation methods performs consistently well for all three classifiers.  It could be the case that the variability is a dominating effect given the small training set, making it difficult to compare extrapolation methods using the 20-to-400 extrapolation task.

Unlike in the face recognition example, we did not resample training classes here, because that would require retraining all of the classifiers–which would be prohibitively time-consuming for the Deep CNN.  Thus, we cannot comment on the robustness of the comparisons from this example, though it is likely that we would obtain different rankings under a new resampling of the training classes.





§ DISCUSSION


In this work, we suggest treating the class set in a classification
task as random, in order to extrapolate classification performance on
a small task to the expected performance on a larger unobserved task.
We show that average generalized accuracy decreases with increased
label set size like the (k-1)th moment of a distribution function.
Furthermore, we introduce an algorithm for estimating this underlying
distribution, that allows efficient computation of higher order
moments. Code for the methods and the simulations can be found in <https://github.com/snarles/ClassEx>.

There are many choices and simplifying assumptions used in
the description of the method.  In this discussion, we discuss these
decisions and map some alternative models or strategies for future
work.




Since our analysis is currently restricted to i.i.d. sampling of classes, 
one direction for future work is to generalize the sampling mechanism,
such as to cluster sampling.  More broadly, the assumption that the labels in 𝒮_k are a
random sample from a homogeneous distribution π may be inappropriate.  Many
natural classification problems arise from hierarchically partitioning
a space of instances into a set of labels.  Therefore, rather than
modeling 𝒮_k as a random sample, it may be more suitable
to model it as a random hierarchical partition of 𝒴, such as one arising from an optional Pólya tree process
<cit.>.
Finally, note that we assume no knowledge about the new class-set
except for its size. Better accuracy might be achieved if some partial information 
is known.



Also, we only discussed extrapolating the
classification accuracy–or equivalently, the risk for the zero-one
cost function.  However, it is possible to extend our analysis to risk
functions with arbitrary cost functions, which is the subject of
forthcoming work.



A third direction of exploration is impose additional modeling assumptions for specific problems.  ClassExReg adopts a non-parametric model of the discriminability function D(u), in the sense that D(u) was defined via a spline expansion.   However, an alternative approach is to assume a parametric family for D(u) defined by a small number of parameters.  In forthcoming work, we show that under certain limiting conditions, D(u) is well-described by a two-parameter family.  This substantially increases the efficiency of estimation in cases where the limiting conditions are well-approximated.




§ ACKNOWLEDGMENTS


We thank Jonathan Taylor, Trevor Hastie, John Duchi, Steve
  Mussmann, Qingyun Sun, Robert Tibshirani, Patrick McClure, and Gal Elidan for useful discussion.  CZ
  is supported by an NSF graduate research fellowship, and would also
  like to thank the European Research Council under the ERC grant agreement n^∘[PSARPS-294519]  for travel support.




