# Count-min sketch with variable number of hash functions: an experimental
study

Éric Fusy  LIGM, CNRS, Univ. Gustave Eiffel, Marne-la-Vallée, France
<EMAIL_ADDRESS>Gregory Kucherov  LIGM, CNRS,
Univ. Gustave Eiffel, Marne-la-Vallée, France
<EMAIL_ADDRESS>
###### Abstract

Conservative Count-Min, a stronger version of the popular Count-Min sketch
[Cormode, Muthukrishnan 2005], is an online-maintained hashing-based sketch
summarizing element frequency information of a stream. Although several works
attempted to analyze the error of conservative Count-Min, its behavior remains
poorly understood. In [Fusy, Kucherov 2022], we demonstrated that under the
uniform distribution of input elements, the error of conservative Count-Min
follows two distinct regimes depending on its load factor.

In this work, we present a series of results providing new insights into the
behavior of conservative Count-Min. Our contribution is twofold. On one hand,
we provide a detailed experimental analysis of Count-Min sketch in different
regimes and under several representative probability distributions of input
elements. On the other hand, we demonstrate improvements that can be made by
assigning a variable number of hash functions to different elements. This
includes, in particular, reduced space of the data structure while still
supporting a small error.

##  1 Introduction

In most general terms, Count-Min sketch is a data structure for representing
an associative array of numbers indexed by elements (keys) drawn from a large
universe, where the array is provided through a stream of (key,value) updates
so that the current value associated to a key is the sum of all previous
updates of this key. Perhaps the most common setting for applying Count-Min,
that we focus on in this paper, is the counting setting where all update
values are $+1$. In this case, the value of a key is its count telling how
many times this key has appeared in the stream. In other words, Count-Min can
be seen as representing a multiset, that is a mapping of a subset of keys to
non-negative integers. With this latter interpretation in mind, each update
will be called insertion. The main supported query of Count-Min is retrieving
the count of a given key, and the returned estimate may not be exact, but can
only overestimate the true count.

The counting version of Count-Min is applied to different practical problems
related to data stream mining and data summarization. One example is tracking
frequent items (heavy hitters) in streams [26, 9, 13]. It occurs in network
traffic monitoring [20], optimization of cache usage [19]. It also occurs in
non-streaming big data applications, e.g. in bioinformatics [28, 3, 32].

Count-Min relies on hash functions but, unlike classic hash tables, does not
store elements but only count information (hence the term sketch). It was
proposed in [14], however a very similar data structure was proposed earlier
in [11] under the name Spectral Bloom filter. The latter, in turn, is closely
related to Counting Bloom filters [22]. In this work, we adopt the definition
of [11] but still call it Count-Min to be consistent with the name commonly
adopted in the literature. A survey on Count-Min can be found e.g. in [12].

In this paper, we study a stronger version of Count-Min called conservative.
This modification of Count-Min was introduced in [20] under the name
conservative update, see [12]. It was also discussed in [11] under the name
minimal increase. Conservative Count-Min provides strictly tighter count
estimates using the same memory and thus strictly outperforms the original
version. The price to pay is the impossibility to deal with deletions
(negative updates), whereas the original Count-Min can handle deletions as
well, provided that the cumulative counts remain non-negative (condition known
as strict turnstile model [26]).

Analysis of error of conservative Count-Min is a difficult problem having
direct consequences on practical applications. Below in Section 2.2 we survey
known related results in more details. In our previous work [23], we
approached this problem through the relationship with random hypergraphs. We
proved, in particular, that if the elements represented in the data structure
are uniformly distributed in the input, the error follows two different
regimes depending on the peelability property of the underlying hash
hypergraph. While properties of random hypergraphs have been known to be
crucially related to some data structures (see Section 2.3), this had not been
known for Count-Min.

Starting out from these results, in this paper we extend and strengthen this
analysis in several ways, providing experimental demonstrations in support of
our claims. Our first goal is to provide a fine analysis of the “anatomy” of
conservative Count-Min, describing its behavior in different regimes. Our main
novel contribution is the demonstration that assigning different number of
hash functions to different elements can significantly improve the error, and,
as a consequence, lead to memory saving. Another major extension concerns the
probability distribution of input elements: here we study non-uniform
distributions as well, in particular step distribution and Zipf’s
distribution, and analyze the behavior of Count-Min for these distributions.
This analysis is important not only because non-uniform distributions commonly
occur in practice, but also because this provides important insights for the
major application of Count-Min: detection of most frequent elements (sometimes
called heavy hitters problem [26, 9, 13]). In particular, we consider the
“small memory regime” (supercritical, in our terminology) when the number of
distinct represented elements is considerably larger than the size of the data
structure, and analyse conditions under which most frequent elements are
evaluated with negligible error. This has direct applications to the frequent
elements problem.

To conclude the introduction, we note that the experimental character of our
analysis does not restrict the generality of our results that hold for a wide
range of parameters. This follows from the general nature of tested
hypotheses, as well as from theoretical justifications based on previous
works.

##  2 Background and related work

###  2.1 Conservative Count-Min: definitions

A Count-Min sketch is a counter array $A$ of size $n$ together with a set of
hash functions mapping elements (keys) of a given universe $U$ to $[1..n]$. In
this work, each element $e\in U$ can in general be assigned a different number
$k_{e}$ of hash functions. Hash functions are assumed fully random, therefore
we assume w.l.o.g. that an element $e$ is assigned hash functions
$h_{1},\ldots,h_{k_{e}}$.

At initialization, counters $A[i]$ are set to $0$. When processing an
insertion of an input element $e$, basic Count-Min increments by $1$ each
counter $A[h_{i}(e)]$, $1\leq i\leq k_{e}$. The conservative version of Count-
Min increments by $1$ only the smallest of all $A[h_{i}(e)]$. That is,
$A[h_{i}(e)]$ is incremented by $1$ if and only if $A[h_{i}(e)]=\min_{1\leq
j\leq k_{e}}\\{A[h_{j}(e)]\\}$ and is left unchanged otherwise.

In both versions, the estimate of the number of occurrences of a queried
element $e$ is computed by $c(e)=\min_{1\leq i\leq k_{e}}\\{A[h_{i}(e)]\\}$.
It is easily seen that for any input sequence of elements, the estimate
computed by original Count-Min is greater than or equal to the one computed by
the conservative version. This follows from the observation that on the same
input, an entry of counter array $A$ under conservative update can never get
larger than the same entry under Count-Min.

In this work, we study the conservative version of Count-Min. Let $H$ denote a
selection of hash functions $H=\\{h_{1},h_{2},\ldots\\}$. Consider an input
sequence $I$ of $N$ insertions and let $E$ be the set of distinct elements
appearing in $I$. The relative error of an element $e$ is defined by
$\mathrm{err}_{H,I}(e)=(c(e)-occ(e))/occ(e)$, where $occ(e)$ is the number of
occurrences of $e$ in the input. The combined error is an average error over
all elements in $I$ weighted by the number of occurrences, i.e.

$\mathrm{err}_{H,I}=\frac{1}{N}\sum_{e\in
E}occ(e)\cdot\mathrm{err}(e)=\frac{1}{N}\sum_{e\in E}(c(e)-occ(e)).$

We assume that $I$ is an i.i.d. random sequence drawn from a probability
distribution on a set of elements $E\subseteq U$. A key parameter is the size
of $E$ relative to the size $n$ of $A$. By analogy to hash tables,
$\lambda=|E|/n$ is called the load factor, or simply the load.

###  2.2 Analysis of conservative Count-Min: prior works

Motivated by applications to traffic monitoring, [7] was probably the first
work devoted to the analysis of conservative Count-Min in the counting
setting. Their model assumed that all $\binom{n}{k}$ counter combinations are
equally likely, where $k$ hash functions are applied to each element. This
implies the regime when $|E|\gg n$. The focus of [7] was on the analysis of
the growth rate of counters, i.e. the average number of counter increments per
insertion, using a technique based on Markov chains and differential
equations. Another approach proposed in [19] simulates a conservative Count-
Min sketch by a hierarchy of ordinary Bloom filters. Obtained error bounds are
not explicit but are expressed via a recursive relation based on false
positive rates of corresponding Bloom filters.

Recent works [5, 4] propose an analytical approach for computing error bounds
depending on element probabilities assumed independent but not necessarily
uniform, in particular leading to improved precision bounds for detecting
heavy hitters. However the efficiency of this technique is more limited when
all element probabilities are small. In particular, if the input distribution
is uniform, their approach does not bring out any improvement over the general
bounds known for original Count-Min.

In our recent work [23], we proposed an analysis of conservative Count-Min
based on its relationship with random hypergraphs. We summarize the main
results of this work below in Section 2.4.

###  2.3 Hash hypergraph

Many hashing-based data structures are naturally associated with hash
hypergraphs so that hypergraph properties are directly related to the proper
functioning of the data structure. This is the case with Cuckoo hashing [30]
and Cuckoo filters [21], Minimal Perfect Hash Functions and Static Functions
[27], Invertible Bloom Lookup Tables [24], and some others. [33] provides an
extended study of relationships between hash hypergraphs and some of those
data structures.

A Count-Min sketch is associated with a hash hypergraph $H=(V,E)$ where
$V=\\{1..n\\}$ and $E=\\{\\{h_{1}(e),...h_{k_{e}}(e)\\}\\}$ over all distinct
input elements $e$. We use notation $\mathcal{H}_{n,m}$ for hypergraphs with
$n$ vertices and $m$ edges, and $\mathcal{H}_{n,m}^{k}$ for $k$-uniform such
hypergraphs, where all edges have cardinality $k$. In the latter case, since
our hash functions are assumed fully random, a hash hypergraph is a
$k$-uniform Erdős-Rényi random hypergraph.

As inserted elements are assumed to be drawn from a random distribution, it is
convenient to look at the functioning of a Count-Min sketch as a stochastic
process on the associated hash hypergraph [23]. Each vertex holds a counter
initially set to zero, and therefore each edge is associated with a set of
counters held by the corresponding vertices. Inserting an element consists in
incrementing the minimal counters of the corresponding edge, and retrieving
the estimate of an element returns the minimum value among the counters of the
corresponding edge. From now on in our presentation, we will interchangeably
speak of distinct elements and edges of the associated hash hypergraph, as
well as of counters and vertices. Thus, we will call the vertex value the
value of the corresponding counter, and the edge value the estimate of the
corresponding element. Also, we will speak about the load of a hypergraph
understood as the density $|E|/|V|$.

###  2.4 Hypergraph peelability and phase transition of error

A hypergraph $H=(V,E)$ is called peelable if iterating the following step
starting from $H$ results in the empty graph: if the graph has a vertex of
degree $1$ or $0$, delete this vertex together with the incident edge (if
any). Like many other properties of random hypergraphs, peelability undergoes
a phase transition. Consider the Erdős-Rényi $k$-uniform hypergraph model
where graphs are drawn from $\mathcal{H}_{n,m}^{k}$ uniformly at random. It is
shown in [29] that a phase transition occurs at a (computable) peelability
threshold $\lambda_{k}$: a random graph from $\mathcal{H}_{n,\lambda n}^{k}$
is with high probability (w.h.p.) peelable if $\lambda<\lambda_{k}$, and
w.h.p. non-peelable if $\lambda>\lambda_{k}$. The first values are
$\lambda_{2}=0.5$, $\lambda_{3}\approx 0.818$, $\lambda_{4}\approx 0.772$,
etc., $\lambda_{3}$ being the largest. Note that the case $k=2$ makes an
exception to peelability: for $\lambda<\lambda_{2}$, a negligible fraction of
vertices remain after peeling.

Peelability is known to be directly relevant to certain constructions of
Minimal Perfect Hash Functions [27] as well as to the proper functioning of
Invertible Bloom filters [24]. In [23], we proved that it is relevant to
Count-Min as well.

######  Theorem 1 ([23]).

Consider a conservative Count-Min where each element is hashed using $k$
random hash functions. Assume that the input $I$ of length $N$ is drawn from a
uniform distribution on a set $E\subseteq U$ of elements and let
$\lambda=|E|/n$, where $n$ is the number of counters. If
$\lambda<\lambda_{k}$, then for a randomly chosen element $e$, the error
$(c(e)-occ(e))/occ(e)$ is $o(1)$ w.h.p. when both $n$ and $N/n$ grow.

In the complementary regime $\lambda>\lambda_{k}$, we showed in [23], under
some additional assumptions, that $\mathrm{err}_{H,I}$ is $\Theta(1)$. Thus,
the peelability threshold for random hash hypergraphs corresponds to phase
transition in the error produced by conservative Count-Min for uniform
distribution of input. We call regimes $\lambda<\lambda_{k}$ and
$\lambda>\lambda_{k}$ subcritical and supercritical, respectively.

###  2.5 Variable number of hash functions: mixed hypergraphs

The best peelability threshold $\lambda_{3}\approx 0.818$ can be improved in
at least two different ways. One way is to use a carefully defined class of
hash functions which replace uniform sampling of $k$-edges by a specific non-
uniform sampling. Thus, [17] showed that the peelability threshold can be
increased to $\approx 0.918$ for $k=3$ and up to $\approx 0.999$ for larger
$k$’s if a special class of hypergraphs is used.

Another somewhat surprising idea, that we apply in this paper, is to apply a
different number of hash functions to differents elements, that is to consider
non-uniform hypergraphs. Following [16], [31] showed that non-uniform
hypergraphs may have a larger peelability threshold than uniform ones. More
precisely, [31] showed that mixed hypergraphs with two types of edges of
different cardinalities, each constituting a constant fraction of all edges,
may have a larger peelability threshold: for example, hypergraphs with a
fraction of $\approx 0.887$ of edges of cardinality $3$ and the remaining
edges of cardinality $21$ have the peelability threshold $\approx 0.920$,
larger than the best threshold $0.818$ achieved by uniform hypergraphs. We
adopt the notation of [31] for mixed hypergraphs: by writing $k=(k_{1},k_{2})$
we express that the hypergraph contains edges of cardinality $k_{1}$ and
$k_{2}$, and $k=(k_{1},k_{2};\alpha)$ specifies in addition that the fraction
of $k_{1}$-edges is $\alpha$.

The idea of using different number of hash functions for different elements
has also appeared in data structures design. [8] proposed weighted Bloom
filters which apply a different number of hash functions depending on the
frequency with which elements are queried and on probabilities for elements to
belong to the set. The idea is to assign more hash functions to elements with
higher query frequency and to those with smaller occurrence probability. It is
shown that this leads to a reduced false positive probability, where the
latter is defined to be weighted by query frequencies. This idea was further
refined in [34], and then further in [6], under the name Daisy Bloom filter.

##  3 Results

###  3.1 Uniform distribution

We start with the case where input elements are uniformly distributed, i.e.
edges of the associated hash hypergraph have equal probabilities to be
processed for updates.

####  3.1.1 Subcritical regime

Theorem 1 in conjunction with the results of Section 2.5 leads to the
assumption that using a different number of hash functions for different
elements one could “extend” the regime of $o(1)$ error of Count-Min sketch,
which can be made into a rigorous statement (for simplicity we only give it
with two different edge cardinalities).

######  Theorem 2 (Extends Theorem 1 to mixed cardinalities).

Consider a conservative Count-Min with $n$ counters. Assume that the input of
length $N$ is drawn from a uniform distribution on $E\subseteq U$ and let
$\lambda<\lambda_{k}$. Assume further that elements of $E$ are hashed
according to a mixed hypergraph model $k=(k_{1},k_{2};\alpha)$. Let $c_{k}$ be
the peelability ratio associated to $k$. Then, when $\lambda<c_{k}$, the error
rate of a randomly chosen key is $o(1)$ w.h.p., as both $n$ and $N/n$ grow.

###### Proof.

The crucial point is that the analysis of [23] first establishes a
deterministic result (Lemma 5 in the extended version) which can be stated as
follows: the error rate of a peelable edge $e$ in an arbitrary fixed
hypergraph is $o(1)$, and the rate of convergence to $0$ is bounded in terms
of the number $M_{e}$ of edges whose peeling has to precede the peeling of
$e$. This ensures that, for any model of random hypergraphs, if a random edge
$e$ in the hypergraph is peelable w.h.p. and if $M_{e}$ is $O(1)$ w.h.p., then
the error rate of $e$ is $o(1)$. Then, the analysis performed in [16] ensures
that the model of random hypergraphs with mixed cardinalities has this
property below the peelability threshold. ∎

Figure 1 shows the average relative error as a function of the load factor for
three types of hypergraphs: 2-uniform, 3-uniform and mixed hypergraph where a
$0.885$ fraction of edges are of cardinality 3 and the remaining ones are of
cardinality 14. 2-uniform and 3-uniform hypergraphs illustrate phase
transitions at load factors approaching respectively $0.5$ and $\approx
0.818$, peelability thresholds for $2$-uniform and $3$-uniform hypergraphs
respectively. It is clearly seen that the phase transition for the mixed
hypergraphs occurs at a larger value approaching $\approx 0.898$ which is the
peelability threshold for this class of hypergraphs [31].

Figure 1: $\mathrm{err}_{H,I}$ for small $\lambda=m/n$, for uniform
distribution and different types of hypergraphs: 2-uniform, 3-uniform and
(3,14)-mixed with a fraction of $0.885$ of 3-edges (parameters borrowed from
[31]). Data obtained for $n=1000$. The input size in each experiment is
$5,000$ times the number of edges. Each average is taken over 10 random
hypergraphs.

While this result follows by combining results of [31] and [23], it has not
been observed earlier and has an important practical consequence: using a
varying number of hash functions in Count-Min sketch allows one to increase
the load factor while keeping negligibly small error. In particular, for the
same input, this leads to space saving compared to the uniform case.

Note that parameters $k=(3,14;0.885)$ are borrowed from [31] in order to make
sure that the phase transition corresponds to the peelability threshold
obtained in [31]. In practice, “simpler” parameters can be chosen, for example
we found that $k=(2,5;0.5)$ produces essentially the same curve as
$k=(3,14;0.885)$.

####  3.1.2 Supercritical regime

When the load factor becomes large (supercritical regime), the situation
changes drastically. When the load factor just surpasses the threshold, some
edges are still evaluated with small or zero error, whereas for the other
edges, the error becomes large. This “intermediate regime” has been
illustrated in [23]. Interestingly, edge values are distributed in this regime
in a very peculiar way, concentrating around several values (see Figure 3 in
[23] for illustration). These values must be explained by some graph
structural patterns which remain to be elucidated.

When the load factor goes even larger, the multi-level pattern of edge values
disappears and all edge values become concentrated around the same value. We
call this phenomenon saturation. For example, for $k=3$ saturation occurs at
around $\lambda=6$ (data not shown). Under this regime, the hash hypergraph is
dense enough so that its specific topology is likely to be irrelevant and the
largest counter level “percolates” into all vertex counters. In other words,
all counters grow at the same rate, without any of them “lagging behind”
because of particular graph structural patterns (such as edges containing leaf
vertices).

Bianchi et al. [7] did their analysis under the assumption that each of
$\binom{n}{k}$ edges is equally likely to be processed at each step. This
emulates the situation where the load factor is very large and the hypergraph
is saturated. The focus of [7] is on the growth rate which is the expected
number of counter increments when processing an edge. Obviously, this number
varies between $1$ and $k$, and [7] establishes that larger values of $k$
imply larger growth rates. Note that the growth rate determines the slope of
the linear dependence of $\mathrm{err}_{H,I}$ on $\lambda$.

The case $k=1$ has not been considered by [7]. In this case, the growth rate
is trivially $1$, as each insertion increments exactly one counter.
Furthermore, $\mathrm{err}_{H,I}$ can be easily inferred in this case, as the
error of a given key is defined by the number of keys hashed to the same
counter. The number of such keys is approximated by the Poisson distribution
with parameter $\frac{m}{n}=\lambda$, with expected value $\lambda$. Since
keys are uniformly distributed, $\mathrm{err}_{H,I}=\lambda$ is a good
estimator of the average relative error.

However, it is non-trivial to see how this compares to the error for $k=2$ in
the case of moderate load factors. Figure 2 clarifies the situation. It shows
that $k=1$ produces larger estimates for load factors until a certain point
before going below the estimate for $k=2$. This “intermediate regime” roughly
corresponds to the $n\ln n$ coupon collector bound, i.e. to the regime where a
significant number of counters remain zero. In this case, even if the total
sum of counters is smaller for $k=1$ than for $k=2$, it is more evenly
distributed for the latter case, similar to the well-known power of choice
phenomenon in resource allocation [2]. Since empty counters are irrelevant for
the average relative error, this results in a smaller error for $k=2$ compared
to $k=1$. Note that for $k=1$ we observe that the error curve fits with the
diagonal line. This can be explained by the fact that, for any given key $x$,
the number of other keys falling in the same counter as $x$ (which determines
the error in the estimate for $x$) is classically approximated by a Poisson
random variable of parameter $\lambda=m/n$, of expected value $\lambda$.

Figure 2: $\mathrm{err}_{H,I}$ as a function of $\lambda=m/n$, for uniform
distribution and supercritical regime.

Note that the configuration $k=1$ does not benefit from the advantage of the
conservative update strategy over the regular Count-Min and has a limited
interest from the practical viewpoint. It is however interesting to observe
that $k=2$ still outperforms $k=1$ for limited $\lambda$.

####  3.1.3 Mixed hypergraphs

We have just seen that for sufficiently large loads, the configuration $k=1$
results in the smallest error. Our next question is whether using a varying
number of hash functions (mixed hypergraphs) can make the error even smaller,
by analogy to the subcritical regime where it extends the range of load values
supporting $o(1)$ error. At first glance, this question is not relevant, as
the configuration $k=1$ has obviously the smallest possible sum of counters,
since every insertion increments at least one counter and therefore $k=1$
seems to yield smallest possible estimates.

This argument, however, applies to the regime when the hypergraph is saturated
(all counters are hit), and $\lambda$ is large enough so that counters are
concentrated. Note that as mentioned earlier, counter values approximately
follow a $\mathrm{Poisson}(k\lambda)$ distribution, which is concentrated
around $k\lambda$ when $\lambda$ gets large. Perhaps surprisingly, it turns
out that the error produced by $k=1$ can be made smaller for a large interval
of $\lambda$ by using a varying number of hash functions, before $\lambda$
reaches the saturation point. Figure 3 illustrates this phenomenon. It
compares the error produced by the uniform case $k=1$ and a mixed
configuration with $80\%$ of edges of cardinality $1$ and $20\%$ of edges of
cardinality $3$. The latter produces a smaller error for values of $\lambda$
up to about $50$. Similar to the previous section, this is explained by the
power of choice effect: 3-edges “smooth out” counter values making the
combined error smaller, which in the not-yet-asymptotic regime $\lambda<50$
has a slightly more significant effect than the fact that the presence 3-edges
makes the sum of counter values larger (as some insertions increment more than
one counter).

Figure 3: $\mathrm{err}_{H,I}$ as a function of $\lambda$, for uniform
distribution and supercritical regime. Uniform configuration $k=1$ (blue
curve) vs. mixed configuration $k=(1,3;0.8)$ (orange curve). For clarity, the
difference between the former and the latter is shown in the right plot.

###  3.2 Step distribution

The analysis of the behavior of conservative Count-Min under uniform
distribution of input elements shows that in the supercritical regime, the
error made by the sketch grows linearly with the load factor. This implies a
limited practical utility of the sketch in this regime. On the other hand, in
many practical situations, input elements are not occurring with the same
frequency. This motivates the application of Count-Min to non-uniform
distributions and, in particular, to detection and analysis of frequent
elements in the input stream. One popular problem here is computation of heavy
hitters, where Count-Min sketch have been previously applied [14].

In this section, we focus on the simplest non-uniform distribution – step
distribution – in order to examine the behavior of Count-Min sketch in
presence of elements with different frequencies. Our model is as follows. We
assume that input elements are classified into two groups that we call hot and
cold, where a hot element has a larger appearance probability than a cold one.
Note that we assume that we have a prior knowledge on whether a given element
belongs to hot or cold ones. This setting is similar to the one studied for
Bloom filters augmented with prior membership and query probabilities [6].
Note that our definition of $\mathrm{err}_{H,I}$ assumes that the query
probability of an element and its appearance probability in the input are
equal.

We assume that the load factors of hot and cold elements are $\lambda_{h}$ and
$\lambda_{c}$ respectively. That is, there are $\lambda_{h}n$ hot and
$\lambda_{c}n$ cold edges in the hash hypergraph. $G>1$, called gap factor,
denotes the ratio between probabilities of a hot and a cold element
respectively. Let $p_{h}$ (resp. $p_{c}$) denote the probability for an input
element to be hot (resp. cold). Then $p_{h}/p_{c}=G\lambda_{h}/\lambda_{c}$,
and since $p_{h}+p_{c}=1$, we have

$p_{h}=\frac{G\lambda_{h}}{\lambda_{c}+G\lambda_{h}},~{}~{}p_{c}=\frac{\lambda_{c}}{\lambda_{c}+G\lambda_{h}}.$

For example, if there are 10 times more distinct cold elements than hot ones
($\lambda_{h}/\lambda_{c}=0.1$) but each hot element is 10 times more frequent
than a cold one ($G=10$), than we have about the same fraction of hot and cold
elements in the input ($p_{h}=p_{c}=0.5$).

In the rest of this section, we will be interested in the combined error of
hot elements alone, denoted $\mathrm{errhot}_{H,I}$. If $E_{h}\subseteq E$ is
the subset of hot elements, and $N_{h}$ is the total number of occurrences of
hot elements in the input, then $\mathrm{errhot}_{H,I}$ is defined by

$\mathrm{errhot}_{H,I}=\frac{1}{N_{h}}\sum_{e\in
E_{h}}occ(e)\cdot\mathrm{err}(e)=\frac{1}{N_{h}}\sum_{e\in
E_{h}}(c(e)-occ(e)).$

####  3.2.1 “Interaction” of hot and cold elements

A partition of elements into hot and cold induces the partition of the
underlying hash hypergraph into two subgraphs that we call hot and cold
subgraphs respectively. Since hot elements have larger counts, one might
speculate that counters associated with hot edges are larger than counts of
cold elements and therefore are not incremented by those. Then,
$\mathrm{errhot}_{H,I}$ is entirely defined by the hot subgraph, considered
under the uniform distribution of elements. In particular,
$\mathrm{errhot}_{H,I}$ as a function of $\lambda_{h}$ should behave the same
way as $\mathrm{err}_{H,I}$ for the uniform distribution (see Section 3.1).

This conjecture, however, is not true in general. One reason is that there is
a positive probability that all nodes of a cold edge are incident to hot edges
as well. As a consequence, “hot counters” (i.e. those incident to hot edges)
gain an additional increment due to cold edges, and the latter contribute to
the overestimate of hot edge counts. Figure 4(a) illustrates this point. It
shows, for $k=3$, $\mathrm{errhot}_{H,I}$ as a function of $\lambda_{h}$ in
presence of cold elements with $\lambda_{c}=5$, for the gap value $G=20$. For
the purpose of comparison, the orange curve shows the error for the uniform
distribution (as in Figure 1), that is the error that hot elements would have
if cold elements were not there. We clearly observe the contribution of cold
elements to the error, even in the load interval below the peelability
threshold.

(a) $G=20$

(b) $G=50$

(c) $G=5$

Figure 4: $\mathrm{errhot}_{H,I}$ for $k=3$ depending on $\lambda_{h}$, in
presence of cold elements with $\lambda_{c}=5$ (blue curves) and without any
cold elements (orange curve).

Figure 4(b) illustrates that when the gap becomes larger (here, $G=50$), the
contribution of cold elements diminishes and the curve approaches the one of
the uniform distribution. A larger gap leads to larger values of hot elements
and, as a consequence, to a smaller relative impact of cold ones.

Another reason for which the above conjecture may not hold is the following:
even if the number of hot elements is very small but the gap factor is not
large enough, the cold edges may cause the counters to become large if
$\lambda_{c}$ is large enough, in particular in the saturation regime
described in Section 3.1.2. As a consequence, the “background level” of
counters created by cold edges may be larger than true counts of hot edges,
causing their overestimates. As an example, consider again the configuration
of Figure 4 with $k=3$ and $\lambda_{c}=5$. Figure 2 shows that the cold
elements taken alone would have an error of about 6 on average ($\approx
6.25$, to be precise, data not shown) which means an about $7\times$
overestimate. Since the graph is saturated in this regime (see Section 3.1.2),
this means that most of the counters will be about 7 times larger than counts
of cold edges. Now, if a hot element is only $5$ times more frequent than a
cold one, those will be about $1.4\times$ overestimated, i.e. will have an
error of about 0.4. Referring to Figure 4, this means that the blue curve
“starts off” at a positive error value in this case, instead of zero. This
situation is illustrated in Figure 4(c).

####  3.2.2 Mixed hypergraphs

The analysis above shows that in presence of a “background” formed by large
number of cold elements, the error of hot elements starts growing for much
smaller load factors than without cold elements, even if the latter are much
less frequent than the former. Inspired by results of Section 3.1.1, one may
ask if the interval of negligible error can be extended by employing the idea
of variable number of hash functions. Note that here this idea applies more
naturally by assigning a different number of hash functions to hot and cold
elements.

Figure 5: $\mathrm{errhot}_{H,I}$ as a function of $\lambda_{h}$ for $k=3$,
$\lambda_{c}=5$ and $G=20$ (blue curve, same as in Figure 4(a)) vs. $k=(2,5)$
for hot and cold elements respectively (red curve)

Figure 5 illustrates that this is indeed possible by assigning a smaller
number of hash functions to hot elements and a larger number to cold ones. It
is clearly seen that the interval supporting close-to-zero errors is extended.
This happens because when the hot subgraph is not too dense, increasing the
cardinality of cold edges leads to a higher probability that at least one of
the vertices of such an edge is not incident to a hot edge. As a consequence,
this element does not affect the error of hot edges. For the same reason,
decreasing the cardinality of hot edges (here, from 3 to 2) improves the
error, as this increases the fraction of vertices non-incident to hot edges.

####  3.2.3 Saturation in supercritical regime

In Section 3.1.2 we discussed the saturation regime occurring for large load
values: when the load grows sufficiently large, i.e. the hash hypergraph
becomes sufficiently dense, all counters reach the same level, erasing
distinctions between edges. In this regime, assuming a fixed load (graph
density) and the uniform distribution of input, the edge value depends only on
input size and not on the graph structure (with high probability). It is in
this context that Bianchi et al. [7] studied the growth rate of edge values
depending on input size.

It is an interesting, natural and practically important question whether this
saturation phenomenon holds for non-uniform distributions as well, as it is
directly related to the capacity of distinguishing elements of different
frequency. A full and precise answer to this question is not within the scope
of this work. We believe that the answer is positive at least when the
distribution is piecewise uniform, when edges are partitioned into several
classes and are equiprobable within each class, provided that each class takes
a linear fraction of all elements. Here we illustrate this thesis with the
step distribution.

(a) $2$-uniform

(b) $3$-uniform

(c) $(2,5)$-mixed

Figure 6: Convergence of average estimates of hot and cold elements for
2-uniform (6(a)), 2-uniform (6(b)) and (2,5)-mixed (6(c)) hypergraphs.
$x$-axis shows the total load $\lambda=\lambda_{h}+\lambda_{c}$ with
$\lambda_{h}=0.1\cdot\lambda$ and $\lambda_{c}=0.9\cdot\lambda$ and $G=10$ in
all cases. All experiences were run for $n=1000$ and the input stream of
length $5,000$ times the number of edges (that is, $5\cdot
10^{6}\cdot\lambda$). Each average is taken over 10 random hypergraphs.

Figure 6 illustrates the saturation phenomenon by showing average values of
hot and cold edges ($G=10$) with three different configurations: 2-uniform,
3-uniform, and (2,5)-mixed. Note that the $x$-axis shows here the total load
$\lambda=\lambda_{h}+\lambda_{c}$, where $\lambda_{h}=0.1\cdot\lambda$ and
$\lambda_{c}=0.9\cdot\lambda$. That is, the number of both hot and cold edges
grows linearly when the total number of edges grows.

One can observe that in all configurations, values of hot and cold edges
converge, which is a demonstration of the saturation phenomenon.
Interestingly, the “convergence speed” heavily depends on the configuration:
the convergence is “slower” for uniform configurations, whereas in the mixed
configuration, it occurs right after the small error regime for hot edges.

###  3.3 Zipf’s distribution

Power law distributions are omnipresent in practical applications. The
simplest of those is Zipf’s distribution which is often used as a test case
for different algorithms including Count-Min sketches [15, 7, 19, 10, 4]. We
also mention a recent learning-based variant of CountMin [25] (learning heavy
hitters) and its study under a Zipfian distribution [1, 18]. Under Zipf’s
distribution, element probabilities in descending order are proportional to
$1/i^{\beta}$, where $i$ is the rank of the key and $\beta\geq 0$ is the
skewness parameter. Note that for $\beta=0$, Zipf’s distribution reduces to
the uniform one.

Zipf’s distribution is an important test case for our study as well, as it
forces several (few) most frequent elements to have very large counts and a
large number of elements (heavy tail) to have small counts whose values
decrease only polynomially on the element rank and are therefore of the same
order of magnitude. Bianchi et al. [7, Fig. 1] observed that for Zipf’s
distribution in the supercritical regime, the estimates follow the “waterfall-
type behavior”: the most frequent elements have essentially exact estimates
whereas the other elements have all about the same estimate regardless of
their frequency. Figure 7 illustrates this phenomenon for different skewness
values.

(a) $\beta=0.7$

(b) $\beta=0.5$

(c) $\beta=0.3$

Figure 7: Exact (blue) and estimated (orange) edge values for Zipf’s
distribution as a function on the element frequency rank, plotted in double
log scale. All plots obtained for $n=1000$, $\lambda=5$, $k=2$, and the input
size $50\cdot 10^{6}$. Estimates are averaged over 10 hash function draws.

The waterfall-type behavior for Zipf’s distribution is well explained by the
analysis we developed in the previous sections. The “waterfall pool level” of
values (called error floor in [7]) is the effect of saturation formed by heavy
tail elements. The few “exceptionally frequent” elements are too few to affect
the saturation level (their number is $\ll n$), they turn out to constitute
“peaks” above the level and are thus estimated without error. Naturally,
smaller skewness values make the distribution less steep and reduce the number
of “exceptionally frequent” elements. For example, according to Figure 7, for
$\lambda=5$ and $k=2$, about 50 most frequent elements are evaluated without
error for $\beta=0.7$, about 40 for $\beta=0.5$ and only 5 for $\beta=0.3$.

(a) $k=2$

(b) $k=(2,5;0.2)$

Figure 8: Exact (blue) and estimated (orange) edge values for Zipf’s
distribution as a function on the element frequency rank, for $n=1000$,
$\lambda=2$, and input size $20\cdot 10^{6}$. Values for ranks 20 to 120 only
are shown. Estimates are averaged over 50 hash function draws.

Following our results from previous sections, we studied whether using a
varying number of hash functions can extend the range of frequent elements
estimated with small error. We found that for moderate loads $\lambda$, this
is possible indeed. Figure 8 illustrates this for $\lambda=2$. It shows a
“zoom” around the “break point” (see Figure 7) for $k=2$ vs. $k=(2,5;0.2)$.
For clarity, plots are shown in regular scale and for elements of rank 20 to
120 only. The Figure demonstrates that the case $k=(2,5;0.2)$ provides a
sharper break (see also Figure 6). As a result, even if the saturation level
is higher for $k=(2,5;0.2)$ than for $k=2$, about 70 most frequent elements
are evaluated with small error with $k=(2,5;0.2)$, vs. about 50 with $k=2$.

##  4 Conclusions

In this paper, we presented a series of experimental results providing new
insights into the behavior of conservative Count-Min sketch. Some of them have
direct applications to practical usage of this data structure. Main results
can be summarized as follows.

  * •

For the uniform distribution of input elements, assigning a different number
of hash functions to different elements extends the subcritical regime (range
of load factors $\lambda$) that supports asymptotically vanishing relative
error. This immediately implies space saving for Count-Min configurations
verifying this regime. For non-uniform distributions, varying number of hash
functions allows extending the regime of negligible error for most frequent
elements,

  * •

Under “sufficiently uniform distributions”, including uniform and step
distributions, a Count-Min sketch reaches a saturation regime when $\lambda$
becomes sufficiently large. In this regime, counters become concentrated
around the same value and elements with different frequency become
indistinguishable,

  * •

Frequent elements that can be estimated with small error can be seen as those
which surpass the saturation level formed by the majority of other elements.
For example, in case of Zipf’s distribution, those elements are a few
“exceptionally frequent elements”, whereas the saturation is insured by the
heavy-tail elements. Applying a varying number of hash functions can increase
the number of those elements for moderate loads $\lambda$.

Many of those results lack a precise mathematical analysis. Perhaps the most
relevant to practical usage of Count-Min is the question of estimating the
“waterfall pool level” for different distributions, which is a fundamental
information for heavy-hitter type of applications. Indeed, in a “very
supercritical” regime, estimates tend to converge to the same level, as it was
illustrated in Section 3.2.3 for the step distribution. On the other hand, a
few elements with “abnormally high” frequencies that exceed this level are
those which are evaluated with small error. This is illustrated with Zipf’s
distribution (Section 3.3). Bianchi et al. [7] observed that in the case of
non-uniform distributions of input elements, the “waterfall pool level” is
upper-bounded by the saturation level of counter values when the input
distribution is modeled by a uniform choice among all $\binom{n}{k}$ possible
edges. This latter is computed in [7] using a method based on Markov chains
and differential equations. We believe that this method can be extended to the
case of mixed graphs as well and leave it for future work. However, computing
the “waterfall pool level” for non-uniform distributions including Zipf’s
distribution is an open problem.

## References

  * [1] Anders Aamand, Piotr Indyk, and Ali Vakilian.  (Learned) frequency estimation algorithms under Zipfian distribution.  CoRR, abs/1908.05198, 2019.  URL: http://arxiv.org/abs/1908.05198, arXiv:1908.05198.
  * [2] Yossi Azar, Andrei Z. Broder, Anna R. Karlin, and Eli Upfal.  Balanced allocations.  SIAM J. Comput., 29(1):180–200, 1999.  doi:10.1137/S0097539795288490.
  * [3] Sairam Behera, Sutanu Gayen, Jitender S Deogun, and NV Vinodchandran.  KmerEstimate: A streaming algorithm for estimating k-mer counts with optimal space usage.  In Proceedings of the 2018 ACM International Conference on Bioinformatics, Computational Biology, and Health Informatics, pages 438–447, 2018.
  * [4] Younes Ben Mazziane, Sara Alouf, and Giovanni Neglia.  A Formal Analysis of the Count-Min Sketch with Conservative Updates.  In IEEE INFOCOM WNA 2022 - The second Workshop on Networking Algorithms (WNA), New York, United States, May 2022.  also in arXiv.2203.14549.  doi:10.1109/INFOCOMWKSHPS54753.2022.9798146.
  * [5] Younes Ben Mazziane, Sara Alouf, and Giovanni Neglia.  Analyzing count min sketch with conservative updates.  Computer Networks, 217:109315, 2022.  URL: https://www.sciencedirect.com/science/article/pii/S1389128622003607.
  * [6] Ioana O. Bercea, Jakob Bæk Tejs Houen, and Rasmus Pagh.  Daisy Bloom filters.  CoRR, abs/2205.14894, 2022.  arXiv:2205.14894.
  * [7] Giuseppe Bianchi, Ken Duffy, Douglas J. Leith, and Vsevolod Shneer.  Modeling conservative updates in multi-hash approximate count sketches.  In 24th International Teletraffic Congress, ITC 2012, Kraków, Poland, September 4-7, 2012, pages 1–8. IEEE, 2012.  URL: https://ieeexplore.ieee.org/document/6331813/.
  * [8] Jehoshua Bruck, Jie Gao, and Anxiao Jiang.  Weighted Bloom filter.  In Proceedings 2006 IEEE International Symposium on Information Theory, ISIT 2006, The Westin Seattle, Seattle, Washington, USA, July 9-14, 2006, pages 2304–2308. IEEE, 2006.  doi:10.1109/ISIT.2006.261978.
  * [9] Moses Charikar, Kevin Chen, and Martin Farach-Colton.  Finding frequent items in data streams.  Theoretical Computer Science, 312(1):3–15, 2004.
  * [10] Peiqing Chen, Yuhan Wu, Tong Yang, Junchen Jiang, and Zaoxing Liu.  Precise error estimation for sketch-based flow measurement.  In Proceedings of the 21st ACM Internet Measurement Conference, IMC ’21, page 113–121, New York, NY, USA, 2021. Association for Computing Machinery.  doi:10.1145/3487552.3487856.
  * [11] Saar Cohen and Yossi Matias.  Spectral Bloom filters.  In Alon Y. Halevy, Zachary G. Ives, and AnHai Doan, editors, Proceedings of the 2003 ACM SIGMOD International Conference on Management of Data, San Diego, California, USA, June 9-12, 2003, pages 241–252. ACM, 2003\.  doi:10.1145/872757.872787.
  * [12] Graham Cormode.  Count-Min sketch.  In Ling Liu and M. Tamer Özsu, editors, Encyclopedia of Database Systems, Second Edition. Springer, 2018.  doi:10.1007/978-1-4614-8265-9_87.
  * [13] Graham Cormode and Marios Hadjieleftheriou.  Finding frequent items in data streams.  Proceedings of the VLDB Endowment, 1(2):1530–1541, 2008.
  * [14] Graham Cormode and S. Muthukrishnan.  An improved data stream summary: the count-min sketch and its applications.  Journal of Algorithms, 55(1):58–75, April 2005.  doi:10.1016/j.jalgor.2003.12.001.
  * [15] Graham Cormode and S. Muthukrishnan.  Summarizing and mining skewed data streams.  In Hillol Kargupta, Jaideep Srivastava, Chandrika Kamath, and Arnold Goodman, editors, Proceedings of the 2005 SIAM International Conference on Data Mining, SDM 2005, Newport Beach, CA, USA, April 21-23, 2005, pages 44–55. SIAM, 2005.  doi:10.1137/1.9781611972757.5.
  * [16] Martin Dietzfelbinger, Andreas Goerdt, Michael Mitzenmacher, Andrea Montanari, Rasmus Pagh, and Michael Rink.  Tight thresholds for cuckoo hashing via XORSAT.  In Samson Abramsky, Cyril Gavoille, Claude Kirchner, Friedhelm Meyer auf der Heide, and Paul G. Spirakis, editors, Automata, Languages and Programming, 37th International Colloquium, ICALP 2010, Bordeaux, France, July 6-10, 2010, Proceedings, Part I, volume 6198 of Lecture Notes in Computer Science, pages 213–225. Springer, 2010.  doi:10.1007/978-3-642-14165-2_19.
  * [17] Martin Dietzfelbinger and Stefan Walzer.  Dense peelable random uniform hypergraphs.  In Michael A. Bender, Ola Svensson, and Grzegorz Herman, editors, 27th Annual European Symposium on Algorithms, ESA 2019, September 9-11, 2019, Munich/Garching, Germany, volume 144 of LIPIcs, pages 38:1–38:16. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2019.  doi:10.4230/LIPIcs.ESA.2019.38.
  * [18] Elbert Du, Franklyn Wang, and Michael Mitzenmacher.  Putting the "learning" into learning-augmented algorithms for frequency estimation.  In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, ICML 2021, 18-24 July 2021, Virtual Event, volume 139 of Proceedings of Machine Learning Research, pages 2860–2869. PMLR, 2021.  URL: http://proceedings.mlr.press/v139/du21d.html.
  * [19] Gil Einziger and Roy Friedman.  A formal analysis of conservative update based approximate counting.  In International Conference on Computing, Networking and Communications, ICNC 2015, Garden Grove, CA, USA, February 16-19, 2015, pages 255–259. IEEE Computer Society, 2015.  doi:10.1109/ICCNC.2015.7069350.
  * [20] Cristian Estan and George Varghese.  New directions in traffic measurement and accounting.  In Matthew Mathis, Peter Steenkiste, Hari Balakrishnan, and Vern Paxson, editors, Proceedings of the ACM SIGCOMM 2002 Conference on Applications, Technologies, Architectures, and Protocols for Computer Communication, August 19-23, 2002, Pittsburgh, PA, USA, pages 323–336. ACM, 2002.  doi:10.1145/633025.633056.
  * [21] Bin Fan, Dave G. Andersen, Michael Kaminsky, and Michael D. Mitzenmacher.  Cuckoo filter: Practically better than Bloom.  In Proceedings of the 10th ACM International on Conference on Emerging Networking Experiments and Technologies, CoNEXT ’14, page 75–88, New York, NY, USA, 2014. Association for Computing Machinery.  doi:10.1145/2674005.2674994.
  * [22] Li Fan, Pei Cao, J. Almeida, and A.Z. Broder.  Summary cache: a scalable wide-area web cache sharing protocol.  IEEE/ACM Transactions on Networking, 8(3):281–293, 2000.  doi:10.1109/90.851975.
  * [23] Éric Fusy and Gregory Kucherov.  Phase transition in count approximation by Count-Min sketch with conservative updates.  In Marios Mavronicolas, editor, Proc. of the 13th International Conference on Algorithms and Complexity (CIAC’23), Larnaca (Cyprus), June 13-16, 2023, volume 13898 of Lecture Notes in Computer Science, pages 232 – 246. Springer, 2023.  full version in arxiv:2203.15496, 2022.
  * [24] Michael T Goodrich and Michael Mitzenmacher.  Invertible Bloom lookup tables.  In 2011 49th Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages 792–799. IEEE, 2011.
  * [25] Chen-Yu Hsu, Piotr Indyk, Dina Katabi, and Ali Vakilian.  Learning-based frequency estimation algorithms.  In 7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019. OpenReview.net, 2019.  URL: https://openreview.net/forum?id=r1lohoCqY7.
  * [26] Hongyan Liu, Yuan Lin, and Jiawei Han.  Methods for mining frequent items in data streams: an overview.  Knowledge and information systems, 26(1):1–30, 2011.
  * [27] Bohdan S Majewski, Nicholas C Wormald, George Havas, and Zbigniew J Czech.  A family of perfect hashing methods.  The Computer Journal, 39(6):547–554, 1996.
  * [28] Hamid Mohamadi, Hamza Khan, and Inanc Birol.  ntCard: a streaming algorithm for cardinality estimation in genomics data.  Bioinformatics, 33(9):1324–1330, 2017.
  * [29] Michael Molloy.  Cores in random hypergraphs and Boolean formulas.  Random Structures & Algorithms, 27(1):124–135, 2005.
  * [30] Rasmus Pagh and Flemming Friche Rodler.  Cuckoo hashing.  Journal of Algorithms, 51(2):122–144, 2004.
  * [31] Michael Rink.  Mixed hypergraphs for linear-time construction of denser hashing-based data structures.  In Peter van Emde Boas, Frans C. A. Groen, Giuseppe F. Italiano, Jerzy R. Nawrocki, and Harald Sack, editors, SOFSEM 2013: Theory and Practice of Computer Science, 39th International Conference on Current Trends in Theory and Practice of Computer Science, Špindlerův Mlýn, Czech Republic, January 26-31, 2013. Proceedings, volume 7741 of Lecture Notes in Computer Science, pages 356–368. Springer, 2013.  doi:10.1007/978-3-642-35843-2_31.
  * [32] Yoshihiro Shibuya and Gregory Kucherov.  Set-Min sketch: a probabilistic map for power-law distributions with application to $k$-mer annotation.  bioRxiv, page 2020.11.14.382713, November 2020.  doi:10.1101/2020.11.14.382713.
  * [33] Stefan Walzer.  Random hypergraphs for hashing-based data structures.  PhD thesis, Technische Universität Ilmenau, Germany, 2020.  URL: https://www.db-thueringen.de/receive/dbt_mods_00047127.
  * [34] Xiujun Wang, Yusheng Ji, Zhe Dang, Xiao Zheng, and Baohua Zhao.  Improved weighted Bloom filter and space lower bound analysis of algorithms for approximated membership querying.  In Matthias Renz, Cyrus Shahabi, Xiaofang Zhou, and Muhammad Aamir Cheema, editors, Database Systems for Advanced Applications - 20th International Conference, DASFAA 2015, Hanoi, Vietnam, April 20-23, 2015, Proceedings, Part II, volume 9050 of Lecture Notes in Computer Science, pages 346–362. Springer, 2015.  doi:10.1007/978-3-319-18123-3_21.


Buyun Liang1<EMAIL_ADDRESS>
Tim Mitchell2<EMAIL_ADDRESS>
Ju Sun1<EMAIL_ADDRESS>
1Computer Science & Engineering, University of Minnesota, Minneapolis, USA
2Queens College, City University of New York, New York City, USA

#  NCVX: A General-Purpose Optimization Solver for
Constrained Machine and Deep Learning

###### Abstract

Imposing explicit constraints is relatively new but increasingly pressing in
deep learning, stimulated by, e.g., trustworthy AI that performs robust
optimization over complicated perturbation sets and scientific applications
that need to respect physical laws and constraints. However, it can be hard to
reliably solve constrained deep learning problems without optimization
expertise. The existing deep learning frameworks do not admit constraints.
General-purpose optimization packages can handle constraints but do not
perform auto-differentiation and have trouble dealing with nonsmoothness. In
this paper, we introduce a new software package called NCVX, whose initial
release contains the solver PyGRANSO, a PyTorch-enabled general-purpose
optimization package for constrained machine/deep learning problems, _the
first of its kind_. NCVX inherits auto-differentiation, GPU acceleration, and
tensor variables from PyTorch, and is built on freely available and widely
used open-source frameworks. NCVX is available at https://ncvx.org, with
detailed documentation and numerous examples from machine/deep learning and
other fields.

##  1 Introduction

Mathematical optimization is an indispensable modeling and computational tool
for all science and engineering fields, especially for machine/deep learning.
To date, researchers have developed numerous foolproof techniques, user-
friendly solvers, and modeling languages for convex (CVX) problems, such as
SDPT3 (Toh et al., 1999), Gurobi (Gurobi Optimization, LLC, 2021), Cplex
(Cplex, 2009), TFOCS (Becker et al., 2011), CVX(PY) (Grant et al., 2008;
Diamond and Boyd, 2016), AMPL (Gay, 2015), YALMIP (Lofberg, 2004). These
developments have substantially lowered the barrier of CVX optimization for
non-experts. However, practical problems, especially from machine/deep
learning, are often nonconvex (NCVX), and possibly also constrained (CSTR) and
nonsmooth (NSMT).

There are methods and packages handling NCVX problems in restricted settings:
PyTorch (Paszke et al., 2019) and TensorFlow (Abadi et al., 2015) can solve
large-scale NCVX, NSMT problems without constraints. CSTR problems can be
heuristically turned into penalty forms and solved as unconstrained, but this
may not produce feasible solutions for the original problems. When the
constraints are simple, structured methods such as projected (sub)gradient and
Frank-Wolfe (Sra et al., 2012) can be used. When the constraints are
differentiable manifolds, one can consider manifold optimization methods and
packages, e.g., (Py)manopt (Boumal et al., 2014; Townsend et al., 2016),
Geomstats (Miolane et al., 2020), McTorch (Meghwanshi et al., 2018), and
Geoopt (Kochurov et al., 2020). For general CSTR problems, KNITRO (Pillo and
Roma, 2006) and IPOPT (Wächter and Biegler, 2006) implement interior-point
methods, while ensmallen (Curtin et al., 2021) and GENO (Laue et al., 2019)
rely on augmented Lagrangian methods. However, moving beyond smooth (SMT)
constraints, both of these families of methods, at best, handle only special
types of NSMT constraints. Finally, packages specialized for machine learning,
e.g., scikit-learn (Pedregosa et al., 2011), MLib (Meng et al., 2016) and Weka
(Witten et al., 2011), often use problem-specific solvers that cannot be
easily extended to new formulations.

##  2 The GRANSO and NCVX packages

GRANSO111http://www.timmitchell.com/software/GRANSO/ is among the first
optimization packages that can handle general NCVX, NSMT, CSTR problems
(Curtis et al., 2017):

$\displaystyle\begin{split}&\min_{{\bm{x}}\in\mathbb{R}^{n}}f({\bm{x}}),\;\;\text{s.t.
}c_{i}({\bm{x}})\leq
0,\;\forall\;i\in\mathcal{I};\;\;c_{i}({\bm{x}})=0,\;\forall\;i\in\mathcal{E}.\end{split}$
(1)

Here, the objective $f$ and constraint functions $c_{i}$’s are only required
to be almost everywhere continuously differentiable. GRANSO is based on quasi-
Newton updating with sequential quadratic programming (BFGS-SQP) and has the
following advantages: (1) unified treatment of NCVX problems: no need to
distinguish CVX vs NCVX and SMT vs NSMT problems, similar to typical nonlinear
programming packages; (2) reliable step-size rule: specialized methods for
NSMT problems, such as subgradient and proximal methods, often entail tricky
step-size tuning and require the expertise to recognize the structures (Sra et
al., 2012), while GRANSO chooses step sizes adaptively via a gold-standard
line search; (3) principled stopping criterion: GRANSO stops its iteration by
checking a theory-grounded stationarity condition for NMST problems, whereas
specialized methods are usually stopped when reaching ad-hoc iteration caps.

However, GRANSO users must derive gradients analytically222GRANSO is
implemented in MATLAB and does not support auto-differentiation, although
recent versions of MATLAB have included primitive auto-differentiation
functionalities.  and then provide code for these computations, a process
which is often error-prone in machine learning and impractical for deep
learning. Furthermore, as part of the MATLAB software ecosystem, GRANSO is
generally not compatible with popular machine/deep learning frameworks—mostly
in Python and R—and users’ own existing toolchains. To overcome these issues
and facilitate both high performance and ease of use in machine/deep learning,
we introduce a new software package called NCVX, whose initial release
contains the solver PyGRANSO, a PyTorch-port of GRANSO with several new and
key features: (1) auto-differentiation of all gradients, a critical feature to
make PyGRANSO user-friendly; (2) support for both CPU and GPU computations for
improved hardware acceleration and massive parallelism; (3) support for
general tensor variables including vectors and matrices, as opposed to the
single vector of concatenated optimization variables that GRANSO uses; (4)
integrated support for OSQP (Stellato et al., 2020) and other QP solvers for
respectively computing search directions and the stationarity measure on each
iteration. OSQP generally outperforms commercial QP solvers in terms of
scalability and speed. All of these enhancements are crucial for solving
large-scale machine/deep learning problems. NCVX, licensed under the AGPL V3,
is built entirely on freely available and widely used open-source frameworks;
see https://ncvx.org for documentation and examples.

##  3 Usage examples: dictionary learning and neural perceptual attack

In order to make NCVX friendly to non-experts, we strive to keep the user
input minimal. The user is only required to specify the optimization variables
(names and dimensions of variables) and define the objective and constraint
functions. Here, we briefly demonstrate the usage of PyGRANSO solver on a
couple of machine/deep learning problems.

#### Orthogonal dictionary learning (ODL, Bai et al. (2018))

One hopes to find a “transformation” ${\bm{q}}\in\mathbb{R}^{n}$ to sparsify a
data matrix ${\bm{Y}}\in\mathbb{R}^{n\times m}$:

$\displaystyle\min_{{\bm{q}}\in\mathbb{R}^{n}}\;f({\bm{q}})\doteq
1/m\cdot\norm{{\bm{q}}^{\intercal}{\bm{Y}}}_{1},\quad\text{s.t.}\;\norm{{\bm{q}}}_{2}=1,$
(2)

where the constraint $\norm{{\bm{q}}}_{2}=1$ is to avoid the trivial solution
${\bm{q}}=\bm{0}$. Eq. 2 is NCVX, NSMT, and CSTR: nonsmoothness comes from the
objective, and nonconvexity comes from the constraint. Demos 1 & 2 show the
implementations of ODL in GRANSO and PyGRANSO, respectively. Note that the
analytical gradients of the objective and constraint functions are not
required in PyGRANSO.

⬇

function[f,fg,ci,cig,ce,ceg]=fn(q)

f = 1/m*norm(q’*Y, 1);%obj

fg = 1/m*Y*sign(Y’*q);%obj grad

ci = [];cig = [];%no ineq constr

ce = q’*q - 1; % eq constr

ceg = 2*q; % eq constr grad

end

soln = granso(n,fn);

Demo 1: GRANSO for ODL

⬇

def fn(X_struct):

q = X_struct.q

f = 1/m*norm(q.T@Y, p=1) # obj

ce = pygransoStruct()

ce.c1 = q.T@q - 1 # eq constr

return [f,None,ce]

var_in = {"q": [n,1]}# def variable

soln = pygranso(var_in, fn)

Demo 2: PyGRANSO for ODL

#### Neural perceptual attack (NPA, Laidlaw et al. (2020))

The CSTR deep learning problem, NPA, is shown below:

$\displaystyle\max_{\widetilde{{\bm{x}}}}\;{\mathcal{L}}\left(f\left(\widetilde{{\bm{x}}}\right),y\right),\;\;\text{s.t.}\;\;d\left({\bm{x}},\widetilde{{\bm{x}}}\right)=\norm{\phi\left({\bm{x}}\right)-\phi\left(\widetilde{{\bm{x}}}\right)}_{2}\leq\epsilon.$
(3)

Here, ${\bm{x}}$ is an input image, and the goal is to find its perturbed
version $\widetilde{{\bm{x}}}$ that is perceptually similar to ${\bm{x}}$
(encoded by the constraint) but can fool the classifier $f$ (encoded by the
objective). The loss ${\mathcal{L}}\left(\cdot,\cdot\right)$ is the margin
loss used in Laidlaw et al. (2020). Both $f$ in the objective and $\phi$ in
the constraint are deep neural networks with ReLU activations, making both the
objective and constraint functions NSMT and NCVX. The
$d\left({\bm{x}},\widetilde{{\bm{x}}}\right)$ distance is called the Learned
Perceptual Image Patch Similarity (LPIPS) (Laidlaw et al., 2020; Zhang et al.,
2018b). Demo 3 is the PyGRANSO example for solving Eq. 3. Note that the codes
for data loading, model specification, loss function, and LPIPS distance are
not included here. It is almost impossible to derive analytical subgradients
for Eq. 3, and thus the auto-differentiation feature in PyGRANSO is necessary
for solving it.

⬇

def comb_fn(X_struct):

adv_inputs = X_struct.x_tilde

f = MarginLoss(model(adv_inputs),labels) # obj

ci = pygransoStruct()

ci.c1 = lpips_dists(adv_inputs) - 0.5 # ineq constr. bound eps=0.5

return [f,ci,None] # No eq constr

var_in = {"x_tilde": list(inputs.shape)} # define variable

soln = pygranso(var_in,comb_fn)

Demo 3: PyGRANSO for NPA

##  4 Constrained deep learning applications

In this section, we highlight $4$ families of constrained deep learning
problems with highly nontrivial, often nonsmooth, constraints. These
constraints cannot be easily built into the underlying neural networks.
PyGRANSO can directly solve these problems, and detailed codes and tutorials
are available on https://ncvx.org/examples.

###  4.1 Robustness of deep learning models

In visual recognition, deep neural networks (DNNs) are not robust against
perturbations—either adversarially constructed or naturally occurring—that are
easily discounted by human perception Goodfellow et al. (2015); Hendrycks and
Dietterich (2019). To formalize robustness, one popular way is the
_adversarial loss_ defined as Madry et al. (2017)

$\displaystyle\begin{split}\max_{\mathbf{x}^{\prime}}\ell\left(\mathbf{y},f_{\mathbf{\theta}}(\mathbf{x}^{\prime})\right)\quad\operatorname{s.t.}\;\mathbf{x}^{\prime}\in\Delta(\mathbf{x})=\\{\mathbf{x}^{\prime}\in[0,1]^{n}:d\left(\mathbf{x},\mathbf{x}^{\prime}\right)\leq{\epsilon}\\},\end{split}$
(4)

where $f_{\theta}$ is the DNN model, and $\Delta(\mathbf{x})$ is the set of
allowable perturbations with a radius ${\epsilon}$ measured with respect to
the metric $d$. Early works assume that $\Delta(\mathbf{x})$ is the $\ell_{p}$
norm ball intersected with the natural image box, i.e.,
$\\{\mathbf{x}^{\prime}\in[0,1]^{n}:\norm{\mathbf{x}-\mathbf{x}^{\prime}}_{p}\leq{\epsilon}\\}$,
where $p=1,2,\infty$ are popular choices Madry et al. (2017); Goodfellow et
al. (2015). To capture visually realistic perturbations, recent work has also
modeled nontrivial transformations Hendrycks and Dietterich (2019) that use
non-$\ell_{p}$ metrics. For empirical robustness evaluation (RE), solving Eq.
4 leads to the worst perturbations that could fool $f_{\mathbf{\theta}}$.

An alternative formalism of robustness is the _robustness radius_ (or minimum
distortion radius), defined as the minimal level of perturbation that can
cause $f_{\mathbf{\theta}}$ to change its predicted class:

$\displaystyle\min_{\mathbf{x}^{\prime}\in[0,1]^{n}}\;d\left(\mathbf{x},\mathbf{x}^{\prime}\right)\quad\operatorname{s.t.}\;\max_{i\neq
y}f_{\mathbf{\theta}}^{i}(\mathbf{x}^{\prime})\geq
f_{\mathbf{\theta}}^{y}(\mathbf{x}^{\prime}),$ (5)

where the superscript for $f_{\mathbf{\theta}}$ indexes its elements, i.e.,
the output logits. Solving Eq. 5 produces not only a minimally distorted
perturbation $\mathbf{x}^{\prime}$ but also a robustness radius, making it
another popular choice for RE Croce and Hein (2020a, b). For small and
restricted $f_{\mathbf{\theta}}$ and selected $d$, Eq. 5 can be solved exactly
by mixed integer programming Tjeng et al. (2017); Katz et al. (2017); Bunel et
al. (2020). For general $f_{\mathbf{\theta}}$ and selected $d$, lower bounds
of the robustness radius can be computed Weng et al. (2018a); Zhang et al.
(2018a); Weng et al. (2018b); Lyu et al. (2020). But in general, Eq. 5 is
heuristically solved via gradient-based methods or iterative linearization
Szegedy et al. (2013); Moosavi-Dezfooli et al. (2015); Hein and Andriushchenko
(2017); Carlini and Wagner (2016); Rony et al. (2019); Croce and Hein (2020a);
Pintor et al. (2021).

Eqs. 4 and 5 are difficult constrained problems, especially when $d$ is
sophisticated—necessary for modeling realistic perturbations Hendrycks and
Dietterich (2019); Engstrom et al. (2019); Wong et al. (2019); Laidlaw and
Feizi (2019); Hosseini and Poovendran (2018). Popular exisiting solvers for
them rely on explicit projections onto simple sets, and they will not work
when $d$ is a non-$\ell_{p}$ metric. Also, previous work has shown that the
quality of the solution using these handcrafted methods is sensitive to key
hyperparameters: e.g., step-size schedule and iteration budget Carlini et al.
(2019); Croce and Hein (2020b). With PyGRANSO, we can reliably solve Eqs. 4
and 5 with general $d$’s with minimal hyperparameter tuning Liang et al.
(2022); see https://arxiv.org/pdf/2210.00621 for more details.

###  4.2 Neural structural optimization

Designing physical structures such as bridges, optical devices and airplanes
often boils down to structural optimization, a fundamental family of problems
in physical science and engineering Christensen and Klarbring (2008); Hoyer et
al. (2019); Chandrasekhar et al. (2021); Chandrasekhar and Suresh (2021). A
primitive version of structural optimization takes the form333This form is
also called _topology optimization_.:

$\displaystyle\begin{split}\min_{{\bm{x}},\bm{u}}\bm{u}^{\intercal}{\bm{K}}({\bm{x}})\bm{u}\quad\operatorname{s.t.}{\bm{K}}({\bm{x}})\bm{u}=\bm{f},{\bm{V}}({\bm{x}})\leq
v_{0},{\bm{x}}\in\left\\{0,1\right\\}^{d},\end{split}$ (6)

where ${\bm{x}}\in\mathbb{R}^{d}$ is the vectorized binary design variable
that indicates where the material should be put to form the structure,
$\bm{u}\in\mathbb{R}^{n}$ is the displacement vector, $\mathbf{K}$ is the
global stiffness matrix, and $\bm{f}$ is the vector of external force. For the
constraints, ${\bm{K}}({\bm{x}})\bm{u}=\bm{f}$ encodes the physical laws
(e.g., Hooke’s law) and is often called the _equilibrium constraint_ , and
${\bm{V}}({\bm{x}})=v_{0}$ limits the amount of material to be used.

Recent work Hoyer et al. (2019) has introduced the deep image prior (DIP)
Ulyanov et al. (2018) idea into Eq. 6: $\mathbf{x}$ is reparametrized as
${\bm{x}}=g_{\mathbf{\theta}}(\mathbf{\beta})$, where $g_{\mathbf{\theta}}$ is
a deep neural network parametrized by $\mathbf{\theta}$, and $\mathbf{\beta}$
is a fixed random input. DIP has shown great promise for solving difficult
visual inverse problems (see, e.g., discussions in Li et al. (2021); Wang et
al. (2021); Zhuang et al. (2022)), and can implicitly promote spatial
continuity of the design—which is not explicitly modeled in Eq. 6444The state-
of-the-art methods for structural optimization use smoothing filters during
iteration heuristically to encourage spatial continuity; see, e.g.,
Chandrasekhar et al. (2021).. Also, the overparametrization in DIP tends to
ease the global optimization. Relaxing the integer constraint
${\bm{x}}\in\left\\{0,1\right\\}^{d}$ into a box constraint
${\bm{x}}\in[0,1]^{d}$ as is typically done in structural optimization, we
arrive at

$\displaystyle\begin{split}\min_{\mathbf{\theta},\bm{u}}\bm{u}^{\intercal}{\bm{K}}(g_{\mathbf{\theta}}(\mathbf{\beta}))\bm{u}\quad\operatorname{s.t.}{\bm{K}}(g_{\mathbf{\theta}}(\mathbf{\beta}))\bm{u}=\bm{f},{\bm{V}}(g_{\mathbf{\theta}}(\mathbf{\beta}))\leq
v_{0},g_{\mathbf{\theta}}(\mathbf{\beta})\in[0,1]^{d}\end{split}.$ (7)

In order to solve Eq. 7, recent work Hoyer et al. (2019) transforms Eq. 7 into
an unconstrained problem, which includes 1) eliminating $\mathbf{u}$ from the
physical constraint by solving the linear system
${\bm{K}}(g_{\mathbf{\theta}}(\mathbf{\beta}))\bm{u}=\bm{f}$; 2) enforcing the
couple of constraints ${\bm{V}}(g_{\mathbf{\theta}}(\mathbf{\beta}))\leq
v_{0}$ and $g_{\mathbf{\theta}}(\mathbf{\beta})\in[0,1]^{d}$ via
reparametrization, binary search, and implicit differentiation. By contrast,
using PyGRANSO we can directly solve Eq. 7 without any of the problem-specific
tricks. For realistic design problems, $\mathbf{x}$ needs to be discrete-
valued (it can have more than $2$ discrete values in multi-material design),
and the constraint
${\bm{K}}(g_{\mathbf{\theta}}(\mathbf{\beta}))\bm{u}=\bm{f}$ becomes
nonlinear—$\mathbf{K}$ becomes a nonlinear operator acting on $\mathbf{u}$ due
to the governing PDEs Chandrasekhar et al. (2021); Chandrasekhar and Suresh
(2021). PyGRANSO can easily handle these general cases also, whereas the
tricks used in Hoyer et al. (2019) cannot be generalized.

###  4.3 Orthogonal recurrent neural networks (RNNs)

The exploding and vanishing gradient issues are common in RNNs, and they could
occur whenever the recurrent kernel does not have unit eigenvalues Lezcano-
Casado and Martınez-Rubio (2019); Arjovsky et al. (2016); Harandi and Fernando
(2016); Bansal et al. (2018); Henaff et al. (2016), i.e., the associated
weight matrix is not orthogonal. Recent work proposes imposing orthogonality
constraint on the weight matrix directly Lezcano-Casado and Martınez-Rubio
(2019); Helfrich et al. (2018):

$\displaystyle\begin{split}\min_{\mathbf{\theta}}{\mathcal{L}}\left(f_{\mathbf{\theta}}({\bm{x}}),{\bm{y}}\right)\quad\operatorname{s.t.}{\bm{W}}_{hh}^{\intercal}(\mathbf{\theta}){\bm{W}}_{hh}(\mathbf{\theta})={\bm{I}},\det{\bm{W}}_{hh}(\mathbf{\theta})=1,\end{split}$
(8)

where $f_{\mathbf{\theta}}$ is the RNN parameterized by $\mathbf{\theta}$, and
${\bm{W}}_{hh}$ is the recurrent kernel of the RNN (i.e., a subvector of
$\mathbf{\theta}$). Since the constraints define the famous special orthogonal
group which is a smooth manifold, manifold optimization methods can be
developed to solve Eq. 8 Lezcano-Casado and Martınez-Rubio (2019); Helfrich et
al. (2018). But these methods entail heavy mathematics from differential
geometry, and hence is impractical for non-experts to implement and build on
(there could be additional constraints in real applications). Again, PyGRANSO
is able to directly handle the nonlinear constraints, whether the user
recognizes or not that the constraint forms the special orthogonal group.

###  4.4 Knowledge-aware machine learning (KAML)

KAML concerns learning that incorporates prior knowledge, e.g., physical laws
as constraints. As an example, it can take the general form of minimizing a
loss subject to PDE constraints McClenny and Braga-Neto (2020); Dener et al.
(2020):

$\displaystyle\begin{split}\min_{u(\mathbf{x})}\;{\mathcal{L}(u(\mathbf{x}))}\quad\operatorname{s.t.}\begin{cases}f\left({\bm{x}};\frac{\partial
u}{\partial x_{1}},\ldots,\frac{\partial u}{\partial
x_{d}};\frac{\partial^{2}u}{\partial x_{1}\partial
x_{1}},\ldots,\frac{\partial^{2}u}{\partial x_{1}\partial
x_{d}};\ldots\right)=0,\quad\forall\,{\bm{x}}\in\Omega\\\
{\mathcal{B}}\left(u,{\bm{x}}\right)=0,\quad\forall\,{\bm{x}}\in\partial\Omega\end{cases}\end{split}.$
(9)

Here, $\Omega$ and $\partial\Omega$ denote the domain and its boundary
respectively, and the loss $\mathcal{L}$ depends on the functional variable
$u(\mathbf{x})$ over the domain $\Omega$. The constraints are PDEs, with
boundary and/or initial conditions ${\mathcal{B}}\left(\cdot\right)=0$. If the
loss is a constant, this reduces to solving the PDE problem about
$u(\mathbf{x})$. In supervised learning scenarios, the loss could take the
form of $1/N\cdot\sum_{i=1}^{N}\ell(\mathbf{y}_{i},u(\mathbf{x}_{i}))$, where
$\\{(\mathbf{x}_{i},\mathbf{y}_{i})\\}$ is the training set and $u$ is the
predictor to be learned from the training set.

To solve Eq. 9, one can parametrize $u(\mathbf{x})$ as a neural network, i.e.,
$u(\mathbf{x};\mathbf{\theta})$. This is natural in modern supervised
learning, and is called physics-informed neural networks (PINNs) in the
numerical PDE community Dener et al. (2020); Cuomo et al. (2022); Nandwani et
al. (2019); Chen et al. (2019); Huang et al. (2021); Lu et al. (2021). In
numerical PDEs, classical methods use finite-difference approximations for all
the partial derivatives, whereas the PINN idea directly works with continuous
functions that allow partial derivatives to be computed via auto-
differentiation. This “mesh-free” nature of PINNs holds the promise for high-
precision solutions even for high-dimensional PDEs.

The state-of-the-art methods for solving the DNN-parametrized version of Eq. 9
use penalty methods, Lagrangian methods, and augmented Lagrangian methods Lu
et al. (2021); Cuomo et al. (2022); Dener et al. (2020); McClenny and Braga-
Neto (2020), which often involve delicate tuning of multiple hyperparameters
and could lead to infeasible solutions. By contrast, PyGRANSO stops iterations
by rigorous check of constraint violation and stationarity.

##  5 Roadmap

Although NCVX, with the PyGRANSO solver, already has many powerful features,
we plan to further improve it by adding several major components: (1)
symmetric rank one (SR1): SR1, another major type of quasi-Newton methods,
allows less stringent step-size search and tends to help escape from saddle
points faster by taking advantage of negative curvature directions (Dauphin et
al., 2014); (2) stochastic algorithms: in machine learning, computing with
large-scale datasets often involves finite sums with huge number of terms,
calling for stochastic algorithms for reduced per-iteration cost and better
scalability (Sun, 2019); (3) conic programming (CP): semidefinite programming
and second-order cone programming, special cases of CP, are abundant in
machine learning, e.g., kernel machines (Zhang et al., 2019); (4) minimax
optimization (MMO): MMO is an emerging modeling technique in machine learning,
e.g., generative adversarial networks (GANs) (Goodfellow et al., 2020) and
multi-agent reinforcement learning (Jin et al., 2020).

## References

  * Abadi et al. (2015) Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  URL https://www.tensorflow.org/.  Software available from tensorflow.org.
  * Arjovsky et al. (2016) Martin Arjovsky, Amar Shah, and Yoshua Bengio.  Unitary evolution recurrent neural networks.  In _International conference on machine learning_ , pages 1120–1128. PMLR, 2016.
  * Bai et al. (2018) Yu Bai, Qijia Jiang, and Ju Sun.  Subgradient descent learns orthogonal dictionaries.  _arXiv preprint arXiv:1810.10702_ , October 2018.
  * Bansal et al. (2018) Nitin Bansal, Xiaohan Chen, and Zhangyang Wang.  Can we gain more from orthogonality regularizations in training deep networks?  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Becker et al. (2011) Stephen R. Becker, Emmanuel J. Candès, and Michael C. Grant.  Templates for convex cone problems with applications to sparse signal recovery.  _Mathematical programming computation_ , 3(3):165, 2011.
  * Boumal et al. (2014) Nicolas Boumal, Bamdev Mishra, Pierre-Antoine Absil, and Rodolphe Sepulchre.  Manopt, a Matlab toolbox for optimization on manifolds.  _Journal of Machine Learning Research_ , 15(42):1455–1459, 2014.
  * Bunel et al. (2020) Rudy Bunel, P Mudigonda, Ilker Turkaslan, P Torr, Jingyue Lu, and Pushmeet Kohli.  Branch and bound for piecewise linear neural network verification.  _Journal of Machine Learning Research_ , 21(2020), 2020\.
  * Carlini and Wagner (2016) Nicholas Carlini and David Wagner.  Towards evaluating the robustness of neural networks.  _arXiv:1608.04644_ , August 2016.
  * Carlini et al. (2019) Nicholas Carlini, Anish Athalye, Nicolas Papernot, Wieland Brendel, Jonas Rauber, Dimitris Tsipras, Ian Goodfellow, Aleksander Madry, and Alexey Kurakin.  On evaluating adversarial robustness.  _arXiv:1902.06705_ , February 2019.
  * Chandrasekhar and Suresh (2021) Aaditya Chandrasekhar and Krishnan Suresh.  Tounn: topology optimization using neural networks.  _Structural and Multidisciplinary Optimization_ , 63(3):1135–1149, 2021.
  * Chandrasekhar et al. (2021) Aaditya Chandrasekhar, Saketh Sridhara, and Krishnan Suresh.  Auto: a framework for automatic differentiation in topology optimization.  _Structural and Multidisciplinary Optimization_ , 64(6):4355–4365, 2021.
  * Chen et al. (2019) Di Chen, Yiwei Bai, Wenting Zhao, Sebastian Ament, John M Gregoire, and Carla P Gomes.  Deep reasoning networks: Thinking fast and slow.  _arXiv preprint arXiv:1906.00855_ , 2019.
  * Christensen and Klarbring (2008) Peter W. Christensen and Anders Klarbring.  _An Introduction to Structural Optimization_.  Springer Netherlands, 2008.  10.1007/978-1-4020-8666-3.
  * Cplex (2009) IBM ILOG Cplex.  V12. 1: User’s manual for CPLEX.  _International Business Machines Corporation_ , 46(53):157, 2009.
  * Croce and Hein (2020a) Francesco Croce and Matthias Hein.  Minimally distorted adversarial examples with a fast adaptive boundary attack.  In _International Conference on Machine Learning_ , pages 2196–2205. PMLR, 2020a.
  * Croce and Hein (2020b) Francesco Croce and Matthias Hein.  Reliable evaluation of adversarial robustness with an ensemble of diverse parameter-free attacks.  In _International conference on machine learning_ , pages 2206–2216. PMLR, 2020b.
  * Cuomo et al. (2022) Salvatore Cuomo, Vincenzo Schiano Di Cola, Fabio Giampaolo, Gianluigi Rozza, Maizar Raissi, and Francesco Piccialli.  Scientific machine learning through physics-informed neural networks: Where we are and what’s next.  _arXiv preprint arXiv:2201.05624_ , 2022.
  * Curtin et al. (2021) Ryan R. Curtin, Marcus Edel, Rahul Ganesh Prabhu, Suryoday Basak, Zhihao Lou, and Conrad Sanderson.  The ensmallen library for flexible numerical optimization.  _Journal of Machine Learning Research_ , 22(166):1–6, 2021.
  * Curtis et al. (2017) Frank E Curtis, Tim Mitchell, and Michael L Overton.  A BFGS-SQP method for nonsmooth, nonconvex, constrained optimization and its evaluation using relative minimization profiles.  _Optimization Methods and Software_ , 32(1):148–181, 2017.
  * Dauphin et al. (2014) Yann N. Dauphin, Razvan Pascanu, Caglar Gulcehre, Kyunghyun Cho, Surya Ganguli, and Yoshua Bengio.  Identifying and attacking the saddle point problem in high-dimensional non-convex optimization.  In _Proceedings of the 27th International Conference on Neural Information Processing Systems - Volume 2_ , NIPS’14, page 2933–2941, Cambridge, MA, USA, 2014. MIT Press.
  * Dener et al. (2020) Alp Dener, Marco Andres Miller, Randy Michael Churchill, Todd Munson, and Choong-Seock Chang.  Training neural networks under physical constraints using a stochastic augmented lagrangian approach.  _arXiv preprint arXiv:2009.07330_ , 2020.
  * Diamond and Boyd (2016) Steven Diamond and Stephen Boyd.  CVXPY: A Python-embedded modeling language for convex optimization.  _Journal of Machine Learning Research_ , 17(83):1–5, 2016.
  * Engstrom et al. (2019) Logan Engstrom, Brandon Tran, Dimitris Tsipras, Ludwig Schmidt, and Aleksander Madry.  Exploring the landscape of spatial robustness.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pages 1802–1811. PMLR, 09–15 Jun 2019.  URL https://proceedings.mlr.press/v97/engstrom19a.html.
  * Gay (2015) David M. Gay.  The AMPL modeling language: An aid to formulating and solving optimization problems.  In Mehiddin Al-Baali, Lucio Grandinetti, and Anton Purnama, editors, _Numerical analysis and optimization_ , pages 95–116. Springer International Publishing, Cham, 2015.  ISBN 978-3-319-17689-5.
  * Goodfellow et al. (2020) Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial networks.  _Communications of the ACM_ , 63(11):139–144, 2020.
  * Goodfellow et al. (2015) Ian J. Goodfellow, Jonathon Shlens, and Christian Szegedy.  Explaining and harnessing adversarial examples.  _CoRR_ , abs/1412.6572, 2015.
  * Grant et al. (2008) Michael Grant, Stephen Boyd, and Yinyu Ye.  CVX: Matlab software for disciplined convex programming.  http://cvxr.com/cvx, 2008.
  * Gurobi Optimization, LLC (2021) Gurobi Optimization, LLC.  Gurobi Optimizer Reference Manual, 2021.  URL https://www.gurobi.com.
  * Harandi and Fernando (2016) Mehrtash Harandi and Basura Fernando.  Generalized backpropagation,$\backslash$’$\\{$E$\\}$ tude de cas: Orthogonality.  _arXiv preprint arXiv:1611.05927_ , 2016.
  * Hein and Andriushchenko (2017) Matthias Hein and Maksym Andriushchenko.  Formal guarantees on the robustness of a classifier against adversarial manipulation.  _arXiv:1705.08475_ , May 2017.
  * Helfrich et al. (2018) Kyle Helfrich, Devin Willmott, and Qiang Ye.  Orthogonal recurrent neural networks with scaled cayley transform.  In _International Conference on Machine Learning_ , pages 1969–1978. PMLR, 2018.
  * Henaff et al. (2016) Mikael Henaff, Arthur Szlam, and Yann LeCun.  Recurrent orthogonal networks and long-memory tasks.  In _International Conference on Machine Learning_ , pages 2034–2042. PMLR, 2016.
  * Hendrycks and Dietterich (2019) Dan Hendrycks and Thomas Dietterich.  Benchmarking neural network robustness to common corruptions and perturbations.  _arXiv preprint arXiv:1903.12261_ , 2019.
  * Hosseini and Poovendran (2018) Hossein Hosseini and Radha Poovendran.  Semantic adversarial examples.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , pages 1614–1619, 2018.
  * Hoyer et al. (2019) Stephan Hoyer, Jascha Sohl-Dickstein, and Sam Greydanus.  Neural reparameterization improves structural optimization.  _arXiv preprint arXiv:1909.04240_ , 2019.
  * Huang et al. (2021) Dou Huang, Haoran Zhang, Xuan Song, and Ryosuke Shibasaki.  Differentiable projection for constrained deep learning.  _arXiv preprint arXiv:2111.10785_ , 2021.
  * Jin et al. (2020) Chi Jin, Praneeth Netrapalli, and Michael Jordan.  What is local optimality in nonconvex-nonconcave minimax optimization?  In Hal Daumé III and Aarti Singh, editors, _Proceedings of the 37th International Conference on Machine Learning_ , volume 119 of _Proceedings of Machine Learning Research_ , pages 4880–4889. PMLR, 13–18 Jul 2020.
  * Katz et al. (2017) Guy Katz, Clark Barrett, David Dill, Kyle Julian, and Mykel Kochenderfer.  Reluplex: An efficient smt solver for verifying deep neural networks.  _arXiv:1702.01135_ , February 2017.
  * Kochurov et al. (2020) Max Kochurov, Rasul Karimov, and Serge Kozlukov.  Geoopt: Riemannian optimization in PyTorch.  _arXiv preprint arXiv:2005.02819_ , May 2020.
  * Laidlaw and Feizi (2019) Cassidy Laidlaw and Soheil Feizi.  Functional adversarial attacks.  _Advances in neural information processing systems_ , 32, 2019.
  * Laidlaw et al. (2020) Cassidy Laidlaw, Sahil Singla, and Soheil Feizi.  Perceptual adversarial robustness: Defense against unseen threat models.  _arXiv preprint arXiv:2006.12655_ , June 2020.
  * Laue et al. (2019) Sören Laue, Matthias Mitterreiter, and Joachim Giesen.  GENO – GENeric Optimization for classical machine learning.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.
  * Lezcano-Casado and Martınez-Rubio (2019) Mario Lezcano-Casado and David Martınez-Rubio.  Cheap orthogonal constraints in neural networks: A simple parametrization of the orthogonal and unitary group.  In _International Conference on Machine Learning_ , pages 3794–3803. PMLR, 2019.
  * Li et al. (2021) Taihui Li, Zhong Zhuang, Hengyue Liang, Le Peng, Hengkang Wang, and Ju Sun.  Self-validation: Early stopping for single-instance deep generative priors.  _arXiv:2110.12271_ , October 2021.
  * Liang et al. (2022) Hengyue Liang, Buyun Liang, Ying Cui, Tim Mitchell, and Ju Sun.  Optimization for robustness evaluation beyond $\ell_{p}$ metrics.  _arXiv preprint arXiv:2210.00621_ , 2022.
  * Lofberg (2004) Johan Lofberg.  YALMIP : a toolbox for modeling and optimization in MATLAB.  In _2004 IEEE International Conference on Robotics and Automation (IEEE Cat. No.04CH37508)_ , pages 284–289. IEEE, 2004.  10.1109/cacsd.2004.1393890.
  * Lu et al. (2021) Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Karniadakis.  Deepxde: A deep learning library for solving differential equations.  _SIAM Review_ , 63(1):208–228, 2021.
  * Lyu et al. (2020) Zhaoyang Lyu, Ching-Yun Ko, Zhifeng Kong, Ngai Wong, Dahua Lin, and Luca Daniel.  Fastened crown: Tightened neural network robustness certificates.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pages 5037–5044, 2020.
  * Madry et al. (2017) Aleksander Madry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and Adrian Vladu.  Towards deep learning models resistant to adversarial attacks.  _arXiv preprint arXiv:1706.06083_ , 2017.
  * McClenny and Braga-Neto (2020) Levi McClenny and Ulisses Braga-Neto.  Self-adaptive physics-informed neural networks using a soft attention mechanism.  _arXiv preprint arXiv:2009.04544_ , 2020.
  * Meghwanshi et al. (2018) Mayank Meghwanshi, Pratik Jawanpuria, Anoop Kunchukuttan, Hiroyuki Kasai, and Bamdev Mishra.  McTorch, a manifold optimization library for deep learning.  _arXiv preprint arXiv:1810.01811_ , October 2018.
  * Meng et al. (2016) Xiangrui Meng, Joseph Bradley, Burak Yavuz, Evan Sparks, Shivaram Venkataraman, Davies Liu, Jeremy Freeman, DB Tsai, Manish Amde, Sean Owen, Doris Xin, Reynold Xin, Michael J. Franklin, Reza Zadeh, Matei Zaharia, and Ameet Talwalkar.  MLlib: Machine learning in Apache Spark.  _Journal of Machine Learning Research_ , 17(34):1–7, 2016.
  * Miolane et al. (2020) Nina Miolane, Nicolas Guigui, Alice Le Brigant, Johan Mathe, Benjamin Hou, Yann Thanwerdas, Stefan Heyder, Olivier Peltre, Niklas Koep, Hadi Zaatiti, Hatem Hajri, Yann Cabanes, Thomas Gerald, Paul Chauchat, Christian Shewmake, Daniel Brooks, Bernhard Kainz, Claire Donnat, Susan Holmes, and Xavier Pennec.  Geomstats: A Python package for Riemannian geometry in machine learning.  _Journal of Machine Learning Research_ , 21(223):1–9, 2020.
  * Moosavi-Dezfooli et al. (2015) Seyed-Mohsen Moosavi-Dezfooli, Alhussein Fawzi, and Pascal Frossard.  Deepfool: a simple and accurate method to fool deep neural networks.  _arXiv:1511.04599._ , November 2015.
  * Nandwani et al. (2019) Yatin Nandwani, Abhishek Pathak, and Parag Singla.  A primal dual formulation for deep learning with constraints.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  PyTorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32:8026–8037, 2019.
  * Pedregosa et al. (2011) Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake Vanderplas, Alexandre Passos, David Cournapeau, Matthieu Brucher, Matthieu Perrot, and Édouard Duchesnay.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Pillo and Roma (2006) Gianni Pillo and Massimo Roma.  _Large-Scale Nonlinear Optimization_ , volume 83 of _Nonconvex Optimization and Its Applications_.  Springer Science & Business Media, Boston, MA, 2006.
  * Pintor et al. (2021) Maura Pintor, Fabio Roli, Wieland Brendel, and Battista Biggio.  Fast minimum-norm adversarial attacks through adaptive norm constraints.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Rony et al. (2019) Jérôme Rony, Luiz G Hafemann, Luiz S Oliveira, Ismail Ben Ayed, Robert Sabourin, and Eric Granger.  Decoupling direction and norm for efficient gradient-based l2 adversarial attacks and defenses.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 4322–4330, 2019.
  * Sra et al. (2012) Suvrit Sra, Sebastian Nowozin, and Stephen J. Wright.  _Optimization for Machine Learning_.  The MIT Press, Cambridge, MA, 2012.
  * Stellato et al. (2020) Bartolomeo Stellato, Goran Banjac, Paul Goulart, Alberto Bemporad, and Stephen Boyd.  OSQP: an operator splitting solver for quadratic programs.  _Mathematical Programming Computation_ , 12(4):637–672, 2020.  10.1007/s12532-020-00179-2.
  * Sun (2019) Ruoyu Sun.  Optimization for deep learning: theory and algorithms.  _arXiv preprint arXiv:1912.08957_ , December 2019.
  * Szegedy et al. (2013) Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian Goodfellow, and Rob Fergus.  Intriguing properties of neural networks.  _arXiv preprint arXiv:1312.6199_ , 2013.
  * Tjeng et al. (2017) Vincent Tjeng, Kai Xiao, and Russ Tedrake.  Evaluating robustness of neural networks with mixed integer programming.  _arXiv:1711.07356_ , November 2017.
  * Toh et al. (1999) Kim-Chuan Toh, Michael J. Todd, and Reha H. Tütüncü.  SDPT3—a Matlab software package for semidefinite programming, Version 1.3.  _Optimization Methods and Software_ , 11(1-4):545–581, 1999.
  * Townsend et al. (2016) James Townsend, Niklas Koep, and Sebastian Weichwald.  Pymanopt: A Python toolbox for optimization on manifolds using automatic differentiation.  _Journal of Machine Learning Research_ , 17(137):1–5, 2016.
  * Ulyanov et al. (2018) Dmitry Ulyanov, Andrea Vedaldi, and Victor Lempitsky.  Deep image prior.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 9446–9454, 2018.
  * Wang et al. (2021) Hengkang Wang, Taihui Li, Zhong Zhuang, Tiancong Chen, Hengyue Liang, and Ju Sun.  Early stopping for deep image prior.  _arXiv:2112.06074_ , December 2021.
  * Weng et al. (2018a) Lily Weng, Huan Zhang, Hongge Chen, Zhao Song, Cho-Jui Hsieh, Luca Daniel, Duane Boning, and Inderjit Dhillon.  Towards fast computation of certified robustness for relu networks.  In _International Conference on Machine Learning_ , pages 5276–5285. PMLR, 2018a.
  * Weng et al. (2018b) Tsui-Wei Weng, Huan Zhang, Pin-Yu Chen, Jinfeng Yi, Dong Su, Yupeng Gao, Cho-Jui Hsieh, and Luca Daniel.  Evaluating the robustness of neural networks: An extreme value theory approach.  _arXiv preprint arXiv:1801.10578_ , 2018b.
  * Witten et al. (2011) Ian H. Witten, Eibe Frank, and Mark A. Hall.  _Data Mining: Practical Machine Learning Tools and Techniques_.  Morgan Kaufmann, Burlington, MA, third edition, 2011.
  * Wong et al. (2019) Eric Wong, Frank R. Schmidt, and J. Zico Kolter.  Wasserstein adversarial examples via projected sinkhorn iterations.  _arXiv:1902.07906_ , February 2019.
  * Wächter and Biegler (2006) Andreas Wächter and Lorenz T. Biegler.  On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming.  _Mathematical Programming_ , 106(1):25–57, 2006\.  10.1007/s10107-004-0559-y.
  * Zhang et al. (2018a) Huan Zhang, Tsui-Wei Weng, Pin-Yu Chen, Cho-Jui Hsieh, and Luca Daniel.  Efficient neural network robustness certification with general activation functions.  _Advances in neural information processing systems_ , 31, 2018a.
  * Zhang et al. (2018b) Richard Zhang, Phillip Isola, Alexei A. Efros, Eli Shechtman, and Oliver Wang.  The unreasonable effectiveness of deep features as a perceptual metric.  In _2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 586–595, 2018b.
  * Zhang et al. (2019) Richard Y. Zhang, Cédric Josz, and Somayeh Sojoudi.  Conic optimization for control, energy systems, and machine learning: Applications and algorithms.  _Annual Reviews in Control_ , 47:323–340, 2019.
  * Zhuang et al. (2022) Zhong Zhuang, Taihui Li, Hengkang Wang, and Ju Sun.  Blind image deblurring with unknown kernel size and substantial noise.  _arXiv:2208.09483_ , August 2022.