Classifying with Uncertain Data Envelopment Analysis

Classifications are ubiquitous and essential constructs across the human experience, a reality motivated by the fact that we learn through association and disassociation. Indeed, our thoughts are replete with comparisons and contrasts among known and novel aspects of our understanding, and we routinely toil to organize and group information into like collections for myriad purposes. This overriding narrative to categorize information manifests itself in the joys of life, say by classifying literature, art, food, and sport; it promotes scientific progress by delineating the fields of study and by interpreting the results of experimentation; and it aids innovation in fields like engineering, healthcare, and politics as we study how the practical benefits of basic research can promote the general welfare.

The classification concept has not escaped mathematical and practical study, and the process of dividing entities into categories has a long and fruitful history. The literature on how to classify is consequently substantial, well studied, broad in application, and methodologically diverse as illustrated by the following compendious list.

decision trees [34]	integer programming [3]
neural networks [35]	robust statistics [11]
Bayesian networks [25]	discriminant analysis [15]
graph clustering [30]	data mining [23]

[]

Additional classification techniques, along with numerous citations, are reviewed in [6, 17, 19, 29]. This disciplinary range challenges having a categorical language that would promote shared efforts across disciplines, and while all of mathematics, computer science, statistics, probability, data science, operations research, and machine learning dabble with similar classification problems, they each scope their problems in terms apropos for their specific tasks, see [17, 20] for related comments. We use the term classification diffusely to simply mean that entities are partitioned into categories, and we refrain from related terms like clustering and grouping that have more refined definitions in some settings.

We propose a classification technique motivated by uncertain data envelopment analysis (uDEA) that scores a classification by how intra-equitable its categories are. Our fundamental assumptions are that:

1. classifications of the same entities are to be compared on a
   common set of characteristics, and

2. the characteristics are imperfect.

The assumption of imperfect characteristics is surely an anticipated general rule and not a rare oddity; after all, a lack of preciseness taints everything from scientific experimentation, to expert opinion, and to rule based delineation. We do not make any stochastic assumption, and hence, we bypass all hassle of verifying distributional suppositions. That said, we can interpret our model stochastically if such an interpretation appropriately aids analysis. We instead use the modern concept of uncertainty associated with robust optimization to account for imperfect information, and our model builds from the work in [8] to leverage uncertainty toward the goal of establishing equity within a classification’s categories. We comment that data envelopment analysis (DEA) is not new to the processes of clustering and classification, see [4, 26, 32] as illustrative examples.

Sections 2 and 3 define the proximity to equitable efficiency and give an algorithm to approximate its value. We present our classification scheme in Section 4, and we review examples in finance and medicine in Section 5. Note that all norms are $2$-norms unless otherwise noted, although we explicitly denote all norms in the proof of Theorem 1 to avoid confusion.

We define an assessment metric that models the quality of a classification. Assume $O=\{o^t:t=1,2,...T\}$ is a finite set of objects that are to be partitioned into the nonempty disjoint categories $C^s$, for $s=1,2,...,S$. So each $o^t$ is in some $C^s$, and $C=\{C^s:s=1,2,...,S\}$ partitions $O$. We assess a classification $C$ by defining a measure of the intra-similarity of the elements within the individual categories $C^s$, and the classification’s overall quality then equates to the sum of its categorical values.

Characteristics describe individual objects, and we assume that a classification’s quality coincides with some of these characteristics being small and some of them being large. Our language here agrees with the DEA literature, see [5, 10, 16, 22] as reviews, and characteristics whose decrease signifies improvement are called input-characteristics, or more simply inputs, and those whose increase signifies improvement are called output-characteristics, or more simply outputs. We let

We assume there are $N$ inputs and $M$ outputs, making $X^s$ an $N\times |C^s|$ matrix and $Y^s$ an $M\times |C^s|$ matrix.

We employ the concept of efficiency from DEA, and the objects within a category are our decision making units (DMUs). An object’s efficiency within its category is measured as a ratio of weighted sums, with the numerator being an aggregate of the output-characteristics and the denominator being an aggregate of the input-characteristics. The primary goal of DEA is to decide nonnegative weighting parameters to maximize an object’s efficiency while maintaining that each object’s efficiency is no greater than one with the same weights. There are several standard variants of this central theme, and we select the canonical input oriented model with variable returns to scale [5]. We specifically chose to calculate efficiency with the dual of the LP that maximizes the efficiency score for object $t$, which is

provided that

with $o^t\in C^s$ and $x^t$ and $y^t$ being the respective $t$-columns of $X^s$ and $Y^s$. The all ones column vector is $e$, with its length being decided by the context of its use. Note that

is always feasible, and hence, $E^t\le 1$. The value of $E^t$ is the efficiency score of object $t$ relative to its category, and if this score is less than one, then object $t$ is inefficient in its category. Inefficiency means that another object’s efficiency score within the category bests the efficiency score of object $t$ no matter how the inputs and outputs are weighted.

Our assumption of uncertain information indicates that $E^t$ is likely imprecise and that an inefficiency could be the result of inaccurate data. The uDEA problem in [8] addresses this concern by acknowledging and modeling uncertainty and by recasting the efficiency score as a robust linear program. The solution of the robust linear program maximizes an object’s efficiency score by selecting a best collection of data from among the uncertain possibilities. Thus an inefficient object based on the original data might become efficient with a prudent and reasonable selection of new data from a set of uncertain options. We note that several modes of robust DEA are now proposed and that applications are increasing, see [24] and its bibliography as a review, as well as the more recent publications [28, 33].

We use an ellipsoidal model of uncertainty for each row of $A^t$, and denoting the $i$-th row of $A^t$ by $A_i^t$, we set

The scalar ${\sigma }_i$ is assumed to be nonnegative, and the matrix $R_i$ has $|C^s|+2$ columns along with a compatible number of rows to align with the dimension of the vector $u$. The matrix $R_i$ defines the structure of uncertainty, whereas ${\sigma }_i$ scales this structure to shrink or enlarge the amount of uncertainty.

The format of $A^t$ imposes a concomitant format on $R_i$. The last two columns of $A^t$ duplicate the input- and output-characteristics of object $t$, and $R_i$ must maintain this format. Assume for illustrative purposes that $t=2\in C^s$ corresponds with the second column of data in $X^s$ and $Y^s$. Then, if

we have

Suppose our model of uncertainty for the inputs assumes a $2\times 2$ identity with regard to $X^s$. Such diagonal models are commonly motivated by probabilistic assumptions such as the characteristics being uncorrelated random variables. We have in this case that

from which

The last column of $R_3$ must be the negative of the second column for the structures of $A_3^2$ and $A_3^2+{\sigma }_3{u}^T{R}_3$ to agree as $u$ varies over the unit disc. Violations to this structural agreement would essentially, and erroneously, equate to the second input-characteristic of the second object having two values as the data ranged over the uncertainty set. Similar requirements for the output-characteristics are necessary, with illustrative (non-identity) options for $A^2$ being

We express $R_i$ as $\left[R_i^{\prime }|R_i^{\prime \prime }\right]$ in the forthcoming discussion, with $R_i^{\prime \prime }$ having two columns and being decided by $R_i^{\prime }$ as illustrated.

We let $\sigma$ be the $M+N$ vector whose $i$-th component is ${\sigma }_i$, and we define the robust efficiency score for object $t$ with respect to $\sigma$ as,

where the second math program is the standard second-order cone expression of the robust problem. The structure of each $R_i$ mandates that $E^t\left(\sigma \right)\le 1$ for the same reason that $E^t\le 1$, i.e. because $\hat{\eta }$ in (2) remains feasible.

We require that $E^t\left(\sigma \right)=1$ for sufficiently large $\sigma$, a property that ensures object $t$ is capable of efficiency, or more succinctly capable, if enough uncertainty is assumed. A mathematically contrived example in [8] without ellipsoidal uncertainty shows that $sup\left\{E^t\right(\sigma ):\sigma \ge 0\}$ can be strictly less than one, and hence, capability is not necessarily ensured and requires applicable justification. Theorem 1 gives serviceable assurance.

An immediate consequence of Theorem 1 is that invertible $R_i^{\prime }$ matrices ensure capability, and since common examples set $R_i^{\prime }$ to be positive diagonal matrices or covariance matrices, our assumption of capability is routine. Note also that the rank condition is sufficient but not necessary, and that examples not satisfying the rank condition can still result in each object being capable. Indeed, all but cleverly constructed examples lead to the capability of each object in the authors’ experience, making the assumption of capability prevalent and not inimical to practical application.

We classify objects into categories so that the objects within a category require a minimum amount of uncertainty for each of the objects to have a claim to efficiency. So our motivating question is, how much uncertainty does a category need for all of its objects to have a robust efficiency score of one? We define the proximity to equitable efficiency for category $s$ to answer this question, which is

The constraints of this problem necessitate that each object reaches its maximum efficiency score by selecting characteristics from the uncertainty sets $U_i^t\left({\sigma }_i\right)$, a requirement that imparts equity among the objects within the category relative to efficiency. Note that our proximity model does not require each object to be assessed in the same fashion because each object can combine their characteristic values differently to reach an efficiency score of one. A straightforward interpretation of $P^s$ is that it defines a minimum vicinity over which the characteristic data can vary so that each object is equally and maximally efficient if it were allowed to select values from the uncertain possibilities. We comment that a different sense of equity appears in the multiple objective literature [1, 18].

We say that a nonnegative $\sigma$ is feasible-for-classification for category $s$ if $E^t\left(\sigma \right)=1$ for all $t\in C^s$, that is, if $\sigma$ is feasible for the optimization problem defining $P^s$. Note that we have two properties from (2), those being:

Property 1:	if ${\sigma }^{\prime \prime }\ge {\sigma }^{\prime }$ and $E^t\left({\sigma }^{\prime }\right)=1$, then $E^t\left({\sigma }^{\prime \prime }\right)=1$, and
Property 2:	if ${\sigma }^{\prime }$ is feasible-for-classification and ${\sigma }^{\prime \prime }\ge {\sigma }^{\prime }$,
	then ${\sigma }^{\prime \prime }$ is also feasible-for-classification.

Property 1 shows that object $t$ maintains an efficiency score of one if $\sigma$ surpasses what is required to satisfy $E^t\left(\sigma \right)=1$, which suggests that we can construct a feasible-for-classification $\sigma$ by computing the componentwise maximums of solutions to

Theorem 2 confirms this suggestion and establishes both upper and lower bounds on $P^s$ based on this construction.

Let

and let $\hat{\sigma }\left(\Sigma \right)$ be the nonnegative vector whose components are

Theorem 2 establishes upper and lower bounds on $P^s$ in terms of $\hat{\sigma }\left(\Sigma \right)$.

Two corollaries follow, the first of which is an immediate consequence of Theorem 2.

Theorem 2 can also be strengthened to show that $P^s$ achieves its upper bound should $\hat{\sigma }\left(\Sigma \right)$ be decided by a single object.

We address how these mathematical results guide our computational effort to compute $P^s$ in the next section. However, we first note that we assess a classification by summing the proximity values of the individual categories, giving an overall proximity evaluation for the entire classification. So we appraise a classification $C$ with

and its value indicates how much uncertainty is needed so that all objects have a claim to efficiency within their categories. Classifications with low total proximity values are preferred because they partition objects into sets that are intra-category efficient relative to low amounts of uncertainty.

Solving uDEA problems is generally difficult due to a loss of convexity, see Example 3 in [8] as an illustration, and calculating the proximity to equitable efficiency inherits this difficulty because Corollary 2 equates the value of $P^s$ to solving a uDEA problem under appropriate conditions. A first order algorithm in [8] successfully solves uDEA problems, and we use this algorithm to compute a $\Sigma$ by iteratively solving for each $t\in C^s$,

The resulting $\hat{\sigma }\left(\Sigma \right)$ then bounds $P^s$ according to Theorem 2, and it even computes $P^s$ if $\hat{\sigma }\left(\Sigma \right)$ is defined by a single object, in which case $P^s=\parallel \hat{\sigma }\left(\Sigma \right)\parallel$.

We work to calculate $P^s$ by reducing $\hat{\sigma }\left(\Sigma \right)$ if $P^s$ is not decided by a single object. Let $\Gamma \left(\sigma \right)$ be the total sum of robust efficiency scores within the category,

Then,

where the last equality follows from the fact that $\Gamma \left(\sigma \right)$ is equal to the number of objects in $C^s$ if and only if $\sigma$ is feasible-for-classification. A unit-length improving direction $d$ from $\sigma$ is calculated by solving,

The first constraint ensures that $\sigma +\alpha d$ is marginally feasible-for-classification, and the objective minimizes the directional derivative of $\parallel \sigma {\parallel }^2$ along $d$ – note that we have removed the superfluous scalar of two.

We search along the solution $d$ if the optimal value of (5) is negative, and a straightforward calculation in this case shows that the step size $\alpha$ must satisfy

to ensure that $\sigma +\alpha d$ does not violate the lower bound in Theorem 2. We search this interval with the method of bisection to find the largest value of $\alpha$ guaranteeing $\Gamma (\sigma +\alpha d)\ge |C^s|-\epsilon$, where $\epsilon$ is a suitably small tolerance required by our computational study. We then update $\sigma$, solve (5), and repeat our search for $\alpha$. The process starts with $\sigma =\hat{\sigma }\left(\Sigma \right)$ and terminates once the solution to (5) is nonnegative.

We assumed the existence of $\nabla \Gamma$ in (5), but this vector is approximated componentwise with the finite differences

where $e_i$ is the vector of zeros except for a one in the $i$-th position. We note that this approximation exists even if the gradient doesn’t. The use of a backward difference is important since the more common forward difference would instead use $\Gamma (\sigma +\delta e_i)$ and $\Gamma \left(\sigma \right)$. However, $\sigma$ being feasible-for-classification implies that $\sigma +\delta e_i$ would also be feasible-for-classification, and hence, both $\Gamma (\sigma +\delta e_i)$ and $\Gamma \left(\sigma \right)$ would have the same value and our approximate gradient would be zero. The search direction from (5) would then be $d=-\sigma /\parallel \sigma \parallel$, making

The misguided computational interpretation would be that we could approximate $P^s$ by calculating any $\hat{\sigma }\left(\Sigma \right)$ and then scaling this vector down while satisfying

So the use of a forward difference would limit the search to compute $P^s$ from an initial $\hat{\sigma }\left(\Sigma \right)$ to a single search direction due to the erroneous assumption that $\Gamma \left(\sigma \right)$ is constant if $\sigma$ is in some small neighborhood of $\hat{\sigma }\left(\Sigma \right)$. The use of the backward differences removes this concern and provides a more accurate search strategy.

Our algorithmic approach to compute $P^s$ is:

Step 0:	Initialize	Use the algorithm in [8] to compute a $\Sigma$ and a resulting $\hat{\sigma }\left(\Sigma \right)$. If $\hat{\sigma }\left(\Sigma \right)$ is defined by a single object, then set $P^s=\parallel \hat{\sigma }\left(\Sigma \right)\parallel$ and stop. Otherwise, initialize $\sigma$ to $\hat{\sigma }\left(\Sigma \right)$.
Step 1:	Search Direction	Compute an approximate $\nabla \Gamma \left(\sigma \right)$ with (7) and then use this approximation to calculate a search direction $d$ by solving (5). Stop if $d^T\sigma \ge 0$.
Step 2:	Line Search	Use the method of bisection to calculate the largest $\alpha$ satisfying $\Gamma (\sigma +\alpha d)\ge S-\epsilon$, where $\alpha$ satisfies (6). Update $\sigma$ to $\sigma +\alpha d$ and return to Step 1.

[8pt]

We could add a termination criterion to stop if $\Gamma (\sigma +\alpha d)$ reached its lower bound in Step 2; however, theory indicates that the directional derivative from (5) would be nonnegative in this case, and we use the termination criterion in Step 1 to validate our computational scheme’s agreement with theory. This approach has been computationally successful.

We comment that each evaluation of $\Gamma \left(\sigma \right)$ requires $|C^s|$ robust optimization problems to be solved, and while the algorithmic structure above is standard, assembling the computational effort is less so. Robust solvers are increasingly trustworthy, but still, they lack some of the dependability associated with other problem classes like linear programming [2, 14]. We use Gurobi with many of the numeric options set to their most stringent possibilities, and we validate the optimal status of each solve. Success is routine with the increased numeric focus but is more questionable with default settings. If a problem fails our stringent expectations, then we iteratively relax convergent tolerances while guaranteeing optimal values, which are efficiency scores, to at least six decimal places – so an efficiency score of at least 0.999999 equates to 1.0. The vast majority of solves more stringently reaches 16 decimal places of accuracy.

Our classification problem partitions the $T$ objects in $O$ into the $S$ categories in $C$ in a way that (approximately) minimizes the total proximity $P\left(C\right)$, i.e. we seek to solve

The search space of this problem is combinatorially problematic for even modest problems because its size is the Sterling number of the second kind that counts the number of ways to disjointly partition $T$ objects into $S$ nonempty sets. To illustrate, the search spaces of the examples in the next section are

The immensity of these spaces fundamentally challenges a classification scheme independent of how we assess individual classifications, which prompts the need for heuristics to navigate a search. We divide our search for a (near) optimal classification into an initialization procedure, which efficiently solves several adapted $p$-median problems, and an iterative improvement strategy, which is a neighborhood search that reassigns one object per iteration.

The $p$-median problem partitions objects into sets so that the aggregate distance from the objects within the sets to their representative medians is as small as possible. We begin by considering the entire collection of objects as a single classification, i.e. we do not assume a partition and originally set $C=\left\{O\right\}$. The $p$-median problem depends on a distance between the objects, and we let the distance between objects $t_i$ and $t_j$ be

where ${\sigma }^i$ and ${\sigma }^j$ respectively solve

Again, these $\sigma$ vectors are relative to the entire collection of objects and not to a proposed partition of $S$ categories; so there is no intra-category consideration in the initialization process. The idea is that calculating the minimum amount of uncertainty for each object against the whole set of objects should be a reasonable proxy to seed our search to minimize the sum of the intra-category uncertainties comprising $P\left(C\right)$.

We adapt the $p$-median problem to include constraints that mandate the sizes of the sets. For instance, if we seek to classify $30$ objects into three categories, then we solve an adapted $p$-median problem for sets whose sizes are $(2,2,26)$, $(2,3,25)$, $(2,4,24)$, etc. If the set sizes are $(p_1,p_2,\dots ,p_S)$, then we solve

where $i$ and $j$ index over $\{1,2,\dots ,S\}$, $(p_1^{\prime },p_2^{\prime },\dots ,p_{S^{\prime }}^{\prime })$ lists the unique elements of $(p_1,p_2,\dots ,p_S)$, and $k$ indexes over $\{1,2,\dots ,S^{\prime }\}$. The first three constraints model the traditional $p$-median problem, whereas the next three ensure that the partitioning sets have the specified cardinalities. We note that we do not need to consider singletons because any set with only one or two elements trivially has a proximity value of zero. So we solve the adapted $p$-median for all other summative partitions of $S$, and we calculate $P\left(C\right)$ for each. Our initial classification is the one with the smallest value of $P\left(C\right)$.

Suppose $C=\{C^1,C^2,\dots ,C^S\}$ is our initial incumbent classification. We define a neighbor of $C$ to be any classification that moves a single object from one category to another. So $C^{\prime }$ is a neighbor of $C$ if for some $s^{\prime }$ and $s^{\prime \prime }$, and for some object $o^t\in C^{s^{\prime }}$,

From the current incumbent classification $C$, we solve

by iteratively swapping each object from its current category to the other categories. We accept a solution to (9) as the next incumbent if its value is less than $P\left(C\right)$, but we otherwise terminate with a final classification. This process repeats with the updated incumbent if improvement is found.

We conclude with a few comments about our computational process. Our initialization and subsequent neighborhood search does not guarantee a solution to (8). However, our adapted $p$-median problem efficiently and efficaciously seeds our neighborhood search algorithm so that it converges in a few iterations. Moreover, and beyond combinatorial concerns, shrewd and judicious navigation through the collection of all possible classifications is prudent due to the computational burden of calculating $P\left(C\right)$, which requires thousands of SOCP solves even on modestly sized problems. Our two phase approach demonstrates persuasively on the examples of the next section, and we advocate it as a reasonable and practically justified approximation algorithm. Lastly, we have experimented with numerous minor variations, e.g. by adjusting the distances in our adapted $p$-median problem to $d(t_1,t_2)=\parallel {\sigma }^1-{\sigma }^2{\parallel }_p$ or by tweaking the definition of a classification’s neighborhood. We considered such adjustments in an attempt to expediting our computational task without degrading fidelity to the goal of identifying an optimal classification, but all adjustments either diminished outcomes or lacked computational improvement.

We apply our classification scheme to two problems to demonstrate its employ. The first example partitions the 30 stocks in the Down Jones Industrial Average into three performance tiers, and the second categorizes 42 prostrate treatments that were originally analyzed with DEA in [22] and subsequently with uDEA in [8, 31]. The objects, i.e. stocks and prostate treatments, have one input and one output characteristic in each case. Single input and output examples lend themselves to graphical inspection to help verify efficacy, and they guarantee that the upper and lower bounds on $P\left(C\right)$ are within $\sqrt{2}$ of each other from Theorem 2. Both examples set $R_i^{\prime }=I$, so a feasible-for-classification $\sigma$ exists from Theorem 1.

A few comments about the computational burden are worthwhile. Our code is written in python and is available as a supplement to this article. We use Pyomo [12] to model all optimization problems and Gurobi to solve all models. Our initial computational testing suggested that both examples would take about a week of serial calculation on a modern laptop, which seemed extreme. A parallel version of the code that distributes the calculation over twelve nodes has significantly reduced the calculation time to several hours, with the stock example solving in about six hours and the prostate example solving in about eight. The computational burden of classifying with regard to uncertain data is nonetheless significant and requires careful consideration.

The input and output characteristics of each stock are respectively the semi-deviation and the average annual return, both of which are uncertain estimates of future values. Semi-deviation is a common measure of risk and is the standard deviation of downside returns, whereas the average annual return is a routine measure of reward. The risk reward trade-off is steeped in investing, with risk adverse individuals accepting lower average returns to help safeguard principal and risk accepting individuals jeopardizing principal as they pursue higher returns. The input and output characteristics are calculated from the twelve annual returns ending at the beginning of each month in 2021.

Classifying stocks relative to efficiency with regard to risk and reward does not preference risk or reward but instead considers them as uncertain co-equal trade-offs, and hence, the categories tier stocks into performant classes independent of risk tolerance. A reasonable interpretation is that Tier 1 stocks should be considered over Tier 2 or 3 stocks regardless of one’s palatability toward risk. Tier 2 stocks should likewise be considered over Tier 3 stocks.

Figures 4 through 4 show the progression of our classification process. Tier 1 stocks are green squares, Tier 2 stocks are blue squares, and Tier 3 stocks are red diamonds. The outcome from the initialization process is in Figure 4, and we note that Tier 3 can obviously add stocks from Tier 2 while maintaining a proximity value of zero. The first two improvements identify this fact and reduce the uncertainty of Tier 2 by moving two of its stocks to Tier 3. The third and final improvement moves a Tier 1 stock to Tier 2, and the algorithm identifies no subsequent improvement. Table 1 lists the categories throughout the classification process along with the proximity values. Stocks are listed by their ticker symbol.