# HDMM: Optimizing error of high-dimensional statistical queries under
differential privacy

Ryan McKenna  College of Information & Computer Sciences, The University of
Massachusetts, Amherst, MA 10002<EMAIL_ADDRESS>,  Gerome Miklau
College of Information & Computer Sciences, The University of Massachusetts,
Amherst, MA 10002<EMAIL_ADDRESS>,  Michael Hay  Department of Computer
Science, Colgate University, Hamilton, NY 13346<EMAIL_ADDRESS>and  Ashwin
Machanavajjhala  Department of Computer Science, Duke University, Durham, NC,
27708<EMAIL_ADDRESS>
###### Abstract.

In this work we describe the High-Dimensional Matrix Mechanism (HDMM), a
differentially private algorithm for answering a workload of predicate
counting queries. HDMM represents query workloads using a compact implicit
matrix representation and exploits this representation to efficiently optimize
over (a subset of) the space of differentially private algorithms for one that
is unbiased and answers the input query workload with low expected error. HDMM
can be deployed for both $\epsilon$-differential privacy (with Laplace noise)
and $(\epsilon,\delta)$-differential privacy (with Gaussian noise), although
the core techniques are slightly different for each. We demonstrate
empirically that HDMM can efficiently answer queries with lower expected error
than state-of-the-art techniques, and in some cases, it nearly matches
existing lower bounds for the particular class of mechanisms we consider.

##  1\. Introduction

Institutions like the U.S. Census Bureau and Medicare regularly release
summary statistics about individuals, including population statistics cross-
tabulated by demographic attributes [7, 31] and tables reporting on hospital
discharges organized by medical condition and patient characteristics [20].
These data have the potential to reveal sensitive information, especially
through joint analysis of multiple releases [26, 15, 35]. Differential privacy
[10, 11] offers a framework for releasing statistical summaries of sensitive
datasets, while providing formal and quantifiable privacy to the contributing
individuals.

We consider the problem of _batch query answering_ under differential privacy.
That is, our goal is to release answers to a given query _workload_ ,
consisting of a set of _predicate counting queries_ , while satisfying
differential privacy. A predicate counting query computes the number of
individuals in the dataset who satisfy an arbitrary predicate $\phi$ (e.g.,
how many individuals have Income $\geq\$50,\\!000$). Workloads of predicate
counting queries are quite versatile as they are capable of expressing
histograms, multi-dimensional range queries, group-by queries, data cubes,
marginals, and arbitrary combinations thereof. Answering a batch of predicate
counting queries has been widely studied by the research community. Past
results have established theoretical lower bounds [5, 17, 30, 24] as well as a
wealth of practical algorithms [46, 43, 25, 47, 39, 34, 21, 42, 41, 33, 32,
44, 40, 23, 8, 2, 38, 9, 22, 19, 4, 34, 45].

One of the simplest mechanisms for answering a workload of queries is to add
carefully calibrated Laplace or Gaussian noise directly to each of the
workload query answers. The noise magnitude is calibrated to a property of the
workload known as its _sensitivity_ , which can be large for some workloads,
resulting in a significant amount of noise. This method fails to adequately
exploit structure in the workload and correlation amongst queries, and thus it
often adds more noise than is strictly necessary to preserve differential
privacy, resulting in suboptimal utility.

A better approach generalizes the basic noise addition mechanism by first
selecting a new set of _strategy_ queries, then measuring the strategy queries
using a noise addition mechanism, and reconstructing answers to the workload
queries from the noisy measurements of the strategy queries. Choosing an
effective query answering strategy (different from the workload) can result in
orders-of-magnitude lower error than the Laplace mechanism, with no cost to
privacy. Many mechanisms for workload answering fall within the select-
measure-reconstruct paradigm [46, 43, 25, 47, 39, 34, 21, 42, 41, 33, 32, 44,
40, 23, 8, 2, 38, 9, 22, 19, 25], differing primarily in the strategy
selection step. We can characterize strategy selection as a search problem
over a space of strategies, distinguishing prior work in terms of key
algorithmic design choices: the search space, the cost function, and the type
of search algorithm (greedy, local, global, etc.). These design choices impact
the two key performance considerations: accuracy and scalability.

At one extreme are techniques that explore a narrow search space, making them
efficient and scalable but not particularly accurate (in particular, their
search space may include accurate strategies only for a limited class of
workloads). For example, HB [33] considers strategies consisting of
hierarchically structured interval queries. It performs a simple search to
find the branching factor of the hierarchical strategy that minimizes an error
measure that assumes the workload consists of all range queries (regardless of
the actual input workload). It is efficient and can scale to higher
dimensions, but it achieves competitive accuracy only when the workload
consists of range queries and the data is low dimensional.

At the other extreme are techniques that search a large space, and adapt to
the workload by finding a strategy within that space that offers low error on
the workload, thereby making them capable of producing a more accurate
strategy for the particular workload. However, this increased accuracy comes
at the cost of high runtime and poor scalability. This is exemplified by the
Matrix Mechanism [25]. The Matrix Mechanism represents the workload and
strategy as a matrix, and the data as a vector. With this representation, the
select, measure, and reconstruct steps can be completely defined in the
language of linear algebra. In addition, there is a simple formula for the
expected error of any selected strategy matrix in terms of elementary matrix
operations. This enables the Matrix Mechanism to select the optimal strategy
(i.e., the one that offers least expected error) by solving a numerical
optimization problem.

Using a matrix to represent a workload is appealing because the representation
is expressive enough to capture an arbitrary collection of predicate counting
queries, and it reveals any structure that may exist between the workload
queries. However, the size of the workload matrix is equal to the number of
queries times the size of the domain, and is infeasible to represent large
workloads defined over multi-dimensional domains as a matrix. Moreover,
solving the optimization problem underlying strategy selection is nontrivial
and expensive. In short, there is no prior work that is accurate for a wide
range of input workloads, and also capable of scaling to large multi-
dimensional domains.

### Overview of approach and contributions

This paper describes the High-Dimensional Matrix Mechanism (HDMM) which is a
practical instantiation of the Matrix Mechanism (MM), capable of scaling to
large multi-dimensional domains. HDMM offers the flexibility and workload-
adaptivity of the Matrix Mechanism, while offering the scalability of simpler
mechanisms. There are a number of innovations that distinguish HDMM from the
Matrix Mechanism:

  * •

The Matrix Mechanism represents query workloads _explicitly_ , as fully
materialized matrices, while HDMM uses a compact _implicit_ matrix
representation. This permits a lossless representation of queries that avoids
a representation exponential in the number of attributes. The implicit
representation consists of sub-workload matrices (usually one per attribute)
which are used as terms in a Kronecker product. Further, we allow the workload
to be expressed as unions of such Kronecker terms. This allows us to represent
large multi-dimensional workloads efficiently while maintaining the key
benefits that the explicit matrix representation offers.

  * •

The numerical optimization problem at the heart of the Matrix Mechanism is
practically infeasible, even for a single attribute with a domain of size 10.
HDMM introduces four optimization routines for strategy selection: $\mbox{\sf
OPT}_{0}$, $\mbox{\sf OPT}_{\otimes}$, $\mbox{\sf OPT}_{+}$, and $\mbox{\sf
OPT}_{\mathsf{M}}$. $\mbox{\sf OPT}_{0}$ is designed for _explicitly_
represented workloads, and can scale to domains as large as $8192$. $\mbox{\sf
OPT}_{\otimes}$, $\mbox{\sf OPT}_{+}$, and $\mbox{\sf OPT}_{\mathsf{M}}$ are
three different techniques for optimizing _implicitly_ represented workloads
(with implicitly represented strategies), and can scale to significantly
larger domains.111$\mbox{\sf OPT}_{\otimes}$ and $\mbox{\sf OPT}_{+}$ have
linear dependence on the number of attributes, while $\mbox{\sf OPT}_{0}$ and
any method that deals with explicitly represented workloads has an exponential
dependence on the number of attributes. These optimization routines differ in
the space of strategies they consider. In all cases, the strategy search space
is chosen so that is expressive enough to encode high-quality strategies,
while also enabling tractable optimization.

  * •

We also propose efficient algorithms for the measure and reconstruct steps of
HDMM. In the Matrix Mechanism, these steps are implemented by performing
matrix operations with the explicit workload and strategy matrices and the
data vector. HDMM exploits the implicit representation of the selected
strategies to significantly speed up these steps.

As a result of these innovations, HDMM achieves high accuracy on a wide
variety of realistic input workloads, in both low and high dimensions. In
fact, in our experiments, we find it has higher accuracy than all prior
select-measure-reconstruct techniques, even on input workloads for which the
prior techniques were specifically designed. It also achieves reasonable
runtime and scales more effectively than prior work that performs non-trivial
optimization (see Section 11 for a detailed scalability evaluation). The main
bottleneck of HDMM is _representing the data in vector form_ , which requires
space proportional to the domain size; HDMM can scale to domains as large as
$10^{9}$.222We show in this paper that in certain special cases, we can bypass
this fundamental limitation.

HDMM was first described by the authors in [27]. This paper provides a more
complete description of HDMM and adds several new technical contributions:

  1. (1)

We generalize and extend HDMM to support $(\epsilon,\delta)$-differential
privacy via Gaussian noise. This is a nontrivial extension: changing noise
distributions fundamentally changes the optimization problems underlying MM
and HDMM. We analyze this change and derive new optimization routines for
strategy selection in this regime. All four of our core optimization routines
$\mbox{\sf OPT}_{0},\mbox{\sf OPT}_{\otimes},\mbox{\sf OPT}_{+}$, and
$\mbox{\sf OPT}_{\mathsf{M}}$ require different changes to support Gaussian
noise.

  2. (2)

We provide new results on the SVD bound [24], a simple formula which provides
a lower bound on the achievable error of the Matrix Mechanism in terms of the
properties of the workload. Specifically, we show how the SVD bound can be
efficiently computed for implicitly represented workloads, and we use the SVD
bound to provide additional theoretical justification for our optimization
routines. We include the SVD bound in experiments to inform the optimality of
the strategies found by our optimization routines.

  3. (3)

We provide a complete description of $\mbox{\sf OPT}_{\mathsf{M}}$, the
optimization routine that searches over marginal query strategies. This is one
of most important optimization routines because it generally produces the best
strategies for marginal query workloads, one of the most common types of
workloads for multi-dimensional data. In addition, we show that for marginal
query workloads, it is sometimes possible to derive the optimal strategy in
closed form, a remarkable new result.

  4. (4)

We outline an approach to bypass the main bottleneck of HDMM, which is
representing the data in vector form and performing the RECONSTRUCT step. In
particular, we show that in certain special-but-common cases, representing the
data in vector form is not necessary. This is achieved by integrating HDMM
with Private-PGM [28], an alternate technique for RECONSTRUCT that produces a
compact factorized representation of the data. This modification allows HDMM
to scale to even higher dimensional domains, far beyond settings where the
data can be represented in vector form.

  5. (5)

We provide a more comprehensive set of experiments, in both low-dimensional
and high-dimensional settings, showing consistent utility improvements over
other mechanisms. We also provide a detailed analysis of the experimental
results.

### Organization

This paper is organized as follows. In Section 2 we provide background on the
data model and query representation. In Section 3 we provide background on
differential privacy, including the Matrix Mechanism. In Section 4, we
describe $\mbox{\sf OPT}_{0}$, an optimization routine that approximately
solves the Matrix Mechanism optimization problem for explicitly represented
workloads. In Section 5, we show how many common workloads over high-
dimensional domains can be implicitly represented in terms of Kronecker
products. In Section 6, we describe $\mbox{\sf OPT}_{\otimes}$ and $\mbox{\sf
OPT}_{+}$: two optimization routines that can effectively optimize implicitly
represented workloads. In Section 7, we show how marginal query workloads can
be represented implicitly, using an even more compact representation than the
one given in Section 5. In Section 8, we describe $\mbox{\sf
OPT}_{\mathsf{M}}$, an optimization routine that can efficiently optimize
implicitly represented workloads with marginal query strategies. In Section
10, we describe the remaining steps of HDMM, including efficient MEASURE and
RECONSTRUCT. We perform a thorough experimental study in Section 11. Related
work is discussed in Section 12 and proofs are provided in the appendix.

##  2\. Data model and query representation

In this section we introduce much of the notation and relevant background on
the data model and query representation required to understand this work. We
use as a running example, and motivating use case, the differentially private
release of a collection of 10 tabulations from the 2010 Summary File 1
(SF1)[7], an important data product based on the Census of Population and
Housing (CPH).

###  2.1. Notation

Symbol | Meaning
---|---
$A$ | Attribute
$t$ | Tuple (database item)
$\phi$ | Predicate
$\Phi$ | Set of predicates
$\mathbb{I}$ | Indicator function
$\mathcal{W}$ | Logical workload
$\bm{x}$ | Data vector
$\bm{q}$ | Query vector
${\bm{W}}$ | (Explicit) Workload matrix
$\bm{A}$ | (Explicit) Strategy matrix

Symbol | Meaning
---|---
$\epsilon,\delta$ | Privacy parameters
$\mathcal{K}$ | Noise addition mechanism
$\mathcal{K}=\mathcal{L}$ | Laplace mechanism
$\mathcal{K}=\mathcal{G}$ | Gaussian mechanism
$\left\|\bm{A}\right\|_{\mathcal{L}}$ |  $L_{1}$ sensitivity
$\left\|\bm{A}\right\|_{\mathcal{G}}$ |  $L_{2}$ sensitivity
$\left\lVert\bm{A}\right\rVert_{F}$ | Frobenius norm
$\otimes$ | Kronecker product
$\mathbb{W}$ | (Implicit) Workload matrix
$\mathbb{A}$ | (Implicit) Strategy matrix

Table 1. Table of notation.

A table of common notations is given in Table 1. In general, we adhere to the
following conventions. Scalars and tuples are lowercase, non-bold. Sets are
uppercase, non-bold. Vectors are lowercase, bold. Matrices are uppercase,
bold. Implicit matrices are uppercase, blackboard bold.

###  2.2. Data and schema

We assume a single-table relational schema $R(A_{1}\dots A_{d})$, where
$attr(R)$ denotes the set of attributes of $R$. Subsets of attributes are
denoted ${\mathcal{A}}\subseteq attr(R)$. Each attribute $A_{i}$ has a finite
domain $dom(A_{i})$ with size $|dom(A_{i})|=n_{i}$. The full domain of $R$ is
$dom(R)=dom(A_{1})\times\dots\times dom(A_{d})$, and has size
$n=\prod_{i}n_{i}$. An instance $I$ of relation $R$ is a multiset whose
elements are tuples in $dom(R)$.

######  Example 1.

The _Person_ relation has the following schema: six boolean attributes
describing Race, two boolean attributes for Hispanic Ethnicity and Sex, Age in
years between 0 and 114, and a Relationship-to-householder field that has 17
values. These queries are on a multidimensional domain of size $2^{6}\times
2\times 2\times 115\times 17=$500480$$. The data also includes a geographic
attribute encoding state (51 values including D.C.). The SF1+ queries are
defined on a domain of size $$500480$\times 51=$25524480$$.

###  2.3. Logical view of queries

Predicate counting queries are a versatile class, consisting of queries that
count the number of tuples satisfying any logical predicate. We define below a
natural logical representation of these queries, distinguished from a
subsequent vector representation.

######  Definition 1 (Predicate counting query).

A predicate on $R$ is a boolean function $\phi:dom(R)\rightarrow\\{0,1\\}$. A
predicate can be used as a counting query on instance $I$ of $R$ whose answer
is $\phi(I)=\sum_{t\in I}\phi(t)$.

A predicate corresponds to a condition in the WHERE clause of an SQL
statement, so in SQL a predicate counting query has the form: SELECT Count(*)
FROM R WHERE $\phi$.

When a predicate refers only to a subset of attributes ${\mathcal{A}}\subset
attr(R)$ we may say that it is defined with respect to ${\mathcal{A}}$ and
annotate it $\phi_{{\mathcal{A}}}:dom({\mathcal{A}})\rightarrow\\{0,1\\}$. If
$\phi_{\mathcal{A}}$ and $\phi_{\mathcal{B}}$ are predicates on attribute sets
${\mathcal{A}}$ and ${\mathcal{B}}$, then their conjunction is a predicate
$\phi_{\mathcal{A}}\wedge\phi_{\mathcal{B}}:dom({\mathcal{A}}\cup{\mathcal{B}})\rightarrow\\{0,1\\}$.

We assume that each query consists of _arbitrarily complex_ predicates on each
attribute, but require that they are combined across attributes with
conjunctions. In other words, each $\phi$ is of the form
$\phi=\phi_{A_{1}}\land\dots\land\phi_{A_{d}}$. This facilitates the compact
implicit representations described in Section 5. One approach to handling
disjunctions (and other more complex query features) is to transform the
schema by merging attributes. We illustrate this below in its application to
the SF1 workload and provide a more general approach to disjunctive queries in
Appendix A.

######  Example 2.

The SF1 workload consists of conjunctive conditions over its attributes, with
the exception of conditions on the six binary race attributes, which can be
complex disjunctions of conjunctions (such as “The number of Persons with two
or more races”). We simply merge the six binary race attributes and treat it
like a single $2^{6}=64$ size attribute (called simply Race). This schema
transformation does not change the overall domain size, but allows every SF1
query to be expressed as a conjunction.

###  2.4. Logical view of query workloads

A workload is a set of predicate counting queries
$\Phi=\\{\phi_{1},\dots,\phi_{m}\\}$. A workload may consist of queries
designed to support a variety of analyses or user needs, as is the case with
the SF1 workload described above. Workloads may also be built from the
sufficient statistics of models, or generated by tools that aid users in
exploring data, or a combination thereof. For the privacy mechanisms
considered here it is preferable for the workload to explicitly mention all
queries of interest, rather than a subset of the queries that could act like a
supporting view, from which the remaining queries of interest could be
computed. Enumerating all queries of interest allows error to be optimized
collectively. In addition, a workload query can be repeated, or equivalently,
weighted, to express the preference for greater accuracy on that query.

######  Example 3.

Our example workload is a subset of queries from SF1 that can be written as
predicate counting queries over a _Person_ relation. (We omit other queries in
SF1 that involve groups of persons organized into households; for brevity we
refer to our selected queries simply as SF1.) Our SF1 workload has 4151
predicate counting queries, each of the form SELECT Count(*) FROM _Person_
WHERE $\phi$, where $\phi$ specifies some combination of demographic
properties (e.g. number of Persons who are Male, over 18, and Hispanic) and
thus each query reports a count at the national level.

######  Example 4.

A workload we call SF1+ consists of the national level queries in SF1 as well
as the same queries at the state level for each of 51 states. We can
succinctly express the state level queries as an additional 4151 queries of
the form: SELECT state, Count(*) FROM _Person_ WHERE $\phi$ GROUP BY state.
Thus, SF1+ can be represented by a total of $4151+4151=8302$ SQL queries. (The
GROUP BY is a succinct way to represent a potentially large set of predicate
counting queries.) The SF1+ queries are defined on a domain of size
$$500480$\times 51=$25524480$$. In addition to their SQL representation, the
SF1 and SF1+ workloads can be naturally expressed in a logical form defined in
Section 5.1. We use $\mathcal{W}_{\mbox{\tiny SF1}}$and
$\mathcal{W}_{\mbox{\tiny SF1+}}$to denote the logical forms of SF1 and SF1+
respectively.

#### Structured multi-dimensional workloads

Multi-dimensional workloads are often defined in a structured form, as
products and unions of products, that we will exploit later in our implicit
representations. Following the notation above, we write $\Phi_{\mathcal{A}}$
to denote a set of predicates, each mentioning only attributes in
${\mathcal{A}}$. For example, the following are common predicate sets defined
over a single attribute $A$ of tuple $t$:
$I$ $\mbox{Identity}_{A}$ $=\\{\>\mathbb{I}[t_{A}=a]\mid a\in dom(A)\>\\}$ $P$
$\mbox{Prefix}_{A}$ $=\\{\>\mathbb{I}[t_{A}\leq a]\mid a\in dom(A)\>\\}$ $R$
$\mbox{AllRange}_{A}$ $=\\{\>\mathbb{I}[a\leq t_{A}\leq b]\mid a,b\in
dom(A),a\leq b\>\\}$ T $\mbox{Total}_{A}$ $=\\{\>\mathbb{I}[\text{True}]\>\\}$

Above, $\mathbb{I}$ is the indicator function, e.g.,
$\mathbb{I}[t_{A}=a]=\phi(t_{A})=\begin{cases}1&\text{if }t_{A}=a\\\
0&\text{otherwise}\end{cases}$.

$\mbox{Identity}_{A}$ contains one predicate for each element of the domain.
Both $\mbox{Prefix}_{A}$ and $\mbox{Range}_{A}$ rely on an ordered $dom(A)$;
they contain predicates defining a CDF (i.e. sufficient to compute the
empirical cumulative distribution function), and the set of all range queries,
respectively. The predicate set $\mbox{Total}_{A}$, consists of a single
predicate, returning True for any $a\in dom(A)$, and thus counting all
records.

We can construct multi-attribute workloads by taking the cross-product of
predicate sets defined for single attributes, and conjunctively combining
individual queries.

######  Definition 2 (Product).

The product of two predicate sets $\Phi_{\mathcal{A}}$ and
$\Phi_{\mathcal{B}}$ is another predicate set
$\Phi_{\mathcal{A}}\times\Phi_{\mathcal{B}}=\\{\phi_{\mathcal{A}}\wedge\phi_{\mathcal{B}}\mid\phi_{\mathcal{A}}\in\Phi_{\mathcal{A}},\phi_{\mathcal{B}}\in\Psi_{\mathcal{B}}\\}$
containing the conjunction of every pair of predicates.

We describe several examples of workloads constructed from products and unions
of products below.

######  Example 5 (Single query as product).

A predicate counting query in the SF1 workload is: SELECT Count(*) FROM
_Person_ WHERE sex=M AND age < 5. We can express this query as a product:
first, define predicate set $\Phi_{1}=\\{\mathbb{I}[t_{sex}=M]\\}$ and
predicate set $\Phi_{2}=\\{\mathbb{I}[t_{age}<5]\\}$. The query is expressed
as the product $\Phi_{1}\times\Phi_{2}$. (We omit Total on the other
attributes for brevity.)

######  Example 6 (GROUP BY query as product).

A GROUP BY query can be expressed as a product by including an Identity
predicate set for each grouping attribute and a singleton predicate set for
each attribute in the WHERE clause. The product would also include Total for
each attribute not mentioned in the query. For example, the query SELECT sex,
age, Count(*) FROM _Person_ WHERE hispanic = TRUE GROUP BY sex, age is
expressed as
$\mbox{I}_{\texttt{sex}}\times\mbox{I}_{\texttt{age}}\times\Phi_{3}$ where
$\Phi_{3}=\\{\mathbb{I}[t_{hispanic}=True]\\}$. This product contains $2\times
115$ counting queries, one for each possible setting of sex and age.

######  Example 7 (Marginal and Prefix-Marginal).

A marginal query workload is defined by the product of one or more Identity
predicates on selected attributes and Total on all other attributes. For
example, $\mbox{I}_{\tt Sex}\times\mbox{I}_{\tt Age}\times\mbox{I}_{\tt
Hispanic}$ contains predicates to compute one three-way marginal, and consists
of $2\times 115\times 2$ counting queries, one for each possible setting of
Sex, Age and Hispanic. This is equivalent to a GROUP BY query on Sex, Age, and
Hispanic with no WHERE clause.

A prefix-marginal query workload is a natural generalization of a marginal
query workload which can be obtained by replacing one or more of the Identity
predicates with Prefix. For example, $\mbox{I}_{\tt Sex}\times\mbox{P}_{\tt
Age}\times\mbox{I}_{\tt HIspanic}$ contains predicates to compute one three-
way prefix-marginal, and consists of $2\times 115\times 2$ counting queries of
the form $\mathbb{I}[t_{sex}=a,t_{age}\leq b,t_{hispanic}=c]$.

Marginals query workloads are very common workloads that make sense for
domains with categorical attributes, while prefix-marginals are more natural
for domains with both categorical and discretized numeric attributes.

######  Example 8 (SF1 Tabulation as Product).

Except for the population total, the queries in the P12 tabulation of the
Census SF1 workload [7] can be described by a single product: $\mbox{I}_{\tt
Sex}\times\mbox{R}_{\tt Age}$ where $\mbox{R}_{\tt Age}$ is a particular set
of range queries including $[0,114],[0,4]$, $[5,9]$, $[10,14],\dots[85,114]$.

#### Unions of products

Our workloads often combine multiple products as a union of the sets of
queries in each product. For example, the queries required to compute all
three-way marginals could be represented as a union of ${d\choose 3}$
workloads, each a product of the Identity predicate set applied to three
attributes. The input to the algorithms that follow is a logical workload
consisting of a union of products, each representing one or possibly many
queries.

######  Definition 3 (Logical workload).

A logical workload $\mathcal{W}=\\{Q_{1}\dots Q_{k}\\}$ consists of a set of
products $Q_{i}$ where each
$Q_{i}=\Phi^{(i)}_{A_{1}}\times\dots\times\Phi^{(i)}_{A_{d}}$.

######  Example 9 (SF1 as union of products).

The SF1 workload can be represented in a logical form, denoted
$\mathcal{W}_{\mbox{\tiny SF1}}$, that consists of a union of $k=4151$
products, each representing a single query. Because these queries are at the
national level, there is a Total predicate set on the State attribute. The
logical form of the SF1+ workload, denoted $\mathcal{W}_{\mbox{\tiny SF1+}}$,
includes those products, plus an additional 4151 products that are identical
except for replacing the Total on State with an Identity predicate set. There
are a total of $k=8302$ products, representing a total of $4151+51\times
4151=$215852$$ predicate counting queries. While this is a direct translation
from the SQL form, this representation can be reduced. First, we can reduce to
$k=4151$ products by simply adding $True$ to the Identity predicate set on
State to capture the national counts. Furthermore, through manual inspection,
we found that both $\mathcal{W}_{\mbox{\tiny SF1}}$ and
$\mathcal{W}_{\mbox{\tiny SF1+}}$ can be even more compactly represented as
the union of 32 products—we use $\mathcal{W}_{\mbox{\tiny SF1}}^{*}$ and
$\mathcal{W}_{\mbox{\tiny SF1+}}^{*}$ to denote more compact logical forms.
This results in significant space savings (as described shortly in Example 12)
and runtime improvements.

###  2.5. Explicit data and query vectorization

The vector representation of predicate counting queries (and the data they are
evaluated on) is central to the select-measure-reconstruct paradigm. The
vector representation of instance $I$ is denoted $\bm{x}_{I}$ (or simply
$\bm{x}$ if the context is clear) and called the data vector.

######  Definition 4 (Data vector).

The data vector representation of an instance $I$, denoted $\bm{x}_{I}$, is a
vector indexed by tuples $t\in dom(R)$, so that
$\bm{x}_{I}(t)=\sum_{t^{\prime}\in dom(R)}\mathbb{I}[t=t^{\prime}]$.

Informally, $\bm{x}(t)$ counts the number of occurences of $t$ in $I$. Note
that, throughout the paper, the representation of the data vector is always
explicit; it is the representation of queries that will be implicit. Every
predicate counting query $\phi$ also has a corresponding vector form.

######  Definition 5 (Vectorized query).

The vector representation of a predicate counting query $\phi$, denoted
$\bm{q}_{\phi}=vec(\phi)$ is a vector indexed by tuples $t\in dom(R)$, so that
$\bm{q}_{\phi}(t)=\phi(t)$.

The function $vec(\phi)$ which transforms a logical predicate into its
corresponding vector form has a simple implementation: compute $\phi(t)$ for
every $t\in dom(R)$ and store the results in a vector. Note that both the data
vector and the vectorized query have size $|dom(R)|=N$.

######  Proposition 1 (Query evaluation).

A predicate counting query can be answered by computing the dot product
between the query vector and data vector:
$\phi(I)=\bm{q}_{\phi}^{T}\bm{x}_{I}$.

To see why this works, observe that $\phi(I)=\sum_{t\in I}\phi(t)=\sum_{t\in
dom(R)}\phi(t)\bm{x}_{I}(t)=\bm{q}_{\phi}^{T}\bm{x}_{I}$. A single predicate
counting query is represented as a vector, so a workload of predicate counting
queries can be represented as a matrix in which queries are rows. For logical
workload $\mathcal{W}$, its (explicit) matrix form is written ${\bm{W}}$, and
the evaluation of the workload is equivalent to the matrix product
${\bm{W}}\bm{x}_{I}$. Note that the size of the workload matrix is $m\times n$
where $m$ is the number of queries, $\bm{x}_{I}$ is $n\times 1$, and the
vector of workload answers is $m\times 1$.

##  3\. Privacy background

Differential privacy is a property of a randomized algorithm that bounds the
difference in output distribution induced by changes to an individual’s data.
Let $nbrs(I)$ be the set of databases differing from $I$ in at most one
record.

######  Definition 6 (Differential Privacy [10]).

A randomized algorithm $\mathcal{K}$ is $(\epsilon,\delta)$-differentially
private if for any instance $I$, any $I^{\prime}\in nbrs(I)$, and any outputs
$O\subseteq Range(\mathcal{K})$,

$Pr[\mathcal{K}(I)\in O]\leq\exp(\epsilon)\times Pr[\mathcal{K}(I^{\prime})\in
O]+\delta$

When $\delta=0$, we say $\mathcal{K}$ is $\epsilon$-differentially private. In
this work we consider algorithms for answering a workload of predicate
counting queries, in which case we invoke the mechanism with the workload
matrix and data vector: $\mathcal{K}({\bm{W}},\bm{x})$. We focus on mechanisms
that give unbiased answers to the workload queries (sometimes called data-
independent mechanisms). Below we define two simple mechanisms of this form:
the Laplace mechanism ($\mathcal{K}=\mathcal{L}$) and the Gaussian mechanism
($\mathcal{K}=\mathcal{G}$).

######  Definition 7 (Laplace mechanism).

Given an $m\times n$ query matrix ${\bm{W}}$, and a noise magnitude $\sigma$,
the Laplace Mechanism $\mathcal{L}$ outputs the vector:
$\mathcal{L}({\bm{W}},\bm{x})={\bm{W}}\bm{x}+\left\|{\bm{W}}\right\|_{\mathcal{L}}\mbox{Lap}(b)^{m}$
where $\mbox{Lap}(b)^{m}$ is a vector of $m$ i.i.d. samples from a Laplace
distribution with scale $b$, and $\left\|{\bm{W}}\right\|_{\mathcal{L}}$
denotes the maximum $L_{1}$ norm ofthe columns of ${\bm{W}}$.

######  Definition 8 (Gaussian mechanism).

Given an $m\times n$ query matrix ${\bm{W}}$, and a noise magnitude $\sigma$,
the Gaussian mechanism $\mathcal{G}$ outputs the vector:
$\mathcal{G}(\bm{A},\bm{x})=\bm{A}\bm{x}+\left\|{\bm{W}}\right\|_{\mathcal{G}}\mbox{Gaus}(\sigma)^{m}$
where $\mbox{Gaus}(\sigma)^{m}$ is a vector of $m$ i.i.d. samples from a
Gaussian distribution with scale $\sigma$ and
$\left\|{\bm{W}}\right\|_{\mathcal{G}}$ denotes the maximum $L_{2}$ norm of
the columns of ${\bm{W}}$.

The quantities $\left\|{\bm{W}}\right\|_{\mathcal{L}}$ and
$\left\|{\bm{W}}\right\|_{\mathcal{G}}$ above are the $L_{1}$ sensitivity and
$L_{2}$ sensitivity of the query set defined by ${\bm{W}}$ respectively, since
they measure the maximum difference in the answers to the queries in
${\bm{W}}$ on any two databases that differ only by a single record [25]. In
the remainder of the paper, we will use
$\left\lVert{\bm{W}}\right\rVert_{\mathcal{K}}$ to denote this sensitivity
norm when referring to a general quantitity that applies for both Laplace
noise ($\mathcal{K}=\mathcal{L}$) and Gaussian noise
($\mathcal{K}=\mathcal{G}$). As long as $b$ or $\sigma$ is sufficiently large,
these two mechanisms can provide differential privacy. The precise conditions
are stated in Proposition 2.

######  Proposition 2 (Privacy of Laplace and Gaussian mechanisms [11, 3]).

The Laplace mechanism is $\epsilon$-differentially private as long as
$b\geq\frac{1}{\epsilon}$ and the Gaussian mechanism is $(\epsilon,\delta)$
differentially private as long as
$\delta\geq\Psi(\frac{1}{2\sigma}-\epsilon\sigma)-\exp{(\epsilon)}\Psi(-\frac{1}{2\sigma}-\epsilon\sigma)$,
where $\Psi$ is the cumulative distribution function of the standard gaussian
distribution.

Above, the condition on $\sigma$ for the Gaussian mechanism corresponds to the
so-called “analytic gaussian mechanism” [3], which _exactly_ calibrates the
minimum $\sigma$ needed to ensure $(\epsilon,\delta)$-differential privacy.
For this method, $\sigma$ can be obtained numerically with a root-finding
algorithm. Other simpler formulas exist for $\sigma$ that are not as tight,
using a classical analysis of the Gaussian mechanism [11] or through Rényi
differential privacy or concentrated differential privacy [29, 6]. It is
straightforward to analyze the error of these two mechanisms since they add
i.i.d. noise with the same (known) variance to all workload queries. In
particular, we use expected total squared error as our error metric:

######  Definition 9 (Expected Error [22]).

Given a $m\times n$ workload matrix ${\bm{W}}$ and a differentially-private
algorithm $\mathcal{K}$, the expected total squared error is:

$\text{TSE}({\bm{W}},\mathcal{K})=\mathbb{E}[\left\lVert{\bm{W}}\bm{x}-\mathcal{K}({\bm{W}},\bm{x})\right\rVert_{2}^{2}]$

where the expectation is taken over the randomness in the privacy mechanism
$\mathcal{K}$.

######  Proposition 3 (Error of Laplace and Gaussian mechanisms).

The Laplace and Gaussian mechanisms are unbiased and have the following
expected error:

$\displaystyle\text{TSE}({\bm{W}},\mathcal{L})$
$\displaystyle=2mb^{2}\left\|{\bm{W}}\right\|_{\mathcal{L}}^{2}$
$\displaystyle\text{TSE}({\bm{W}},\mathcal{G})$
$\displaystyle=m\sigma^{2}\left\|{\bm{W}}\right\|_{\mathcal{G}}^{2}$

As evident by Proposition 3, the expected error of these mechanisms depends
crucially on the sensitivity of the query matrix ${\bm{W}}$, and if this is
large then the total squared error will also be large. We now introduce a
generalization of these mechanisms that often has better expected error.

###  3.1. The matrix mechanism

The core idea of the matrix mechanism is to apply the mechanism $\mathcal{K}$
on a new query matrix $\bm{A}$, then use the noisy answers to the queries in
$\bm{A}$ to estimate answers to the queries in ${\bm{W}}$. The benefits of
this approach will become clear when we reason about the expected error.

######  Definition 10 (Matrix mechanism [22]).

Given a $m\times n$ workload matrix ${\bm{W}}$, a $p\times n$ strategy matrix
$\bm{A}$, and a differentially private algorithm $\mathcal{K}(\bm{A},\bm{x})$
that answers $\bm{A}$ on $\bm{x}$, the mechanism
$\mathcal{M}_{\bm{A},\mathcal{K}}$ outputs the following vector:
$\mathcal{M}_{\bm{A},\mathcal{K}}({\bm{W}},\bm{x})={\bm{W}}\bm{A}^{+}\mathcal{K}(\bm{A},\bm{x})$.

The privacy of the Matrix mechanism follows from the privacy of $\mathcal{K}$,
since that is the only part of the mechanism that has access to the true data.
In this paper, we assume $\mathcal{K}$ is either the Laplace mechanism or the
Gaussian mechanism, although in principle other noise-addition mechanisms are
also possible [13, 17, 25]. Under some mild conditions stated below, the
Matrix mechanism is unbiased and the error can be expressed analytically as
shown below:

######  Proposition 4 (Error of the Matrix mechanism [22]).

The matrix mechanism is unbiased, i.e.,
$\mathbb{E}[\mathcal{M}_{\bm{A},\mathcal{K}}({\bm{W}},\bm{x})]={\bm{W}}\bm{x}$,
and has the following expected error:

$\displaystyle\text{TSE}({\bm{W}},\mathcal{M}_{\bm{A},\mathcal{K}})$
$\displaystyle=\text{TSE}(\bm{A},\mathcal{K})\left\lVert{\bm{W}}\bm{A}^{+}\right\rVert_{F}^{2}$
$\displaystyle\propto\left\|\bm{A}\right\|_{\mathcal{K}}^{2}\left\lVert{\bm{W}}\bm{A}^{+}\right\rVert_{F}^{2}$

as long as ${\bm{W}}\bm{A}^{+}\bm{A}={\bm{W}}$ and $\mathcal{K}$ adds i.i.d.
noise with mean $0$.

When invoked with $\bm{A}={\bm{W}}$, the matrix mechanism is very similar to
base mechanism $\mathcal{K}$, but the error is typically lower. The term
$\left\lVert{\bm{W}}{\bm{W}}^{+}\right\rVert_{F}^{2}$ in the error formula is
equal to the rank of ${\bm{W}}$, which is bounded above by $m$. This implies
that the error of $\mathcal{M}_{{\bm{W}},\mathcal{K}}$ can never be higher
than $\mathcal{K}$. Further, there are often much better strategies to select
than $\bm{A}={\bm{W}}$.

Finding the best strategy $\bm{A}$ for a given workload ${\bm{W}}$ is the main
technical challenge of the matrix mechanism. This strategy selection problem
can be formulated as a constrained optimization problem. However, solving this
problem is computationally expensive, especially when
$\mathcal{K}=\mathcal{L}$, where it is generally infeasible to solve it for
any nontrivial input workload.

###  3.2. Lower bounds on error

An important theoretical question is to identify, or bound, how low the error
of the matrix mechanism can be for a given workload ${\bm{W}}$. This is useful
because finding the strategy with minimum error is a difficult (and often
intractable) problem, but computing a lower bound on error can be done
efficiently. Knowing how low the error can be allows one to compare the error
of a concrete strategy to the lower bound to see how close to optimal it is.
Additionally, the lower bound can be used to make important policy decisions,
such as setting the privacy budget, or whether it is worth investing the
resources to find a good strategy (as opposed to using other types of
mechanisms like data-dependent ones). Li and Miklau studied this problem and
derived the SVD bound [24].

######  Definition 11 (SVD Bound [24]).

Given a $m\times n$ workload ${\bm{W}}$, the singular value bound is:

$SVDB({\bm{W}})=\frac{1}{n}\big{(}\lambda_{1}+\dots+\lambda_{n}\big{)}^{2}$

where $\lambda_{1},\dots,\lambda_{n}$ are the singular values of ${\bm{W}}$.

The SVD bound gives a lower bound on the error achievable by the matrix
mechanism.

######  Proposition 5 (SVD Bound [24]).

Given a $m\times n$ workload ${\bm{W}}$ and a $p\times n$ strategy $\bm{A}$
that supports ${\bm{W}}$:

$SVDB({\bm{W}})\leq\left\|\bm{A}\right\|_{\mathcal{K}}^{2}\left\lVert{\bm{W}}\bm{A}^{+}\right\rVert_{F}^{2}$

The SVD Bound is known to be tight when $\mathcal{K}=\mathcal{G}$ and some
conditions on ${\bm{W}}$ are satisfied, meaning there is some strategy
$\bm{A}$ that achieves the equality. When $\mathcal{K}=\mathcal{L}$ the bound
may not be tight however. In this paper we use the SVD Bound to evaluate the
quality of the strategies found by our optimization routines. We also derive
expressions for efficiently computing the SVD Bound for implicitly represented
workloads, and use this analysis to theoretically justify our optimization
routines.

##  4\. Optimizing Explicit Workloads

In this section, we introduce the main optimization problem that underlies
strategy selection for the matrix mechanism, and describe $\mbox{\sf
OPT}_{0}$, our algorithm for solving it. We assume for now that the workload
is represented _explicitly_ as a dense matrix. The methods we describe are
useful by themselves for workloads defined over modest domains (namely, those
smaller than about $n=10^{4}$), and they are an essential building block for
the more scalable methods we describe in Section 6.

###  4.1. The optimization problem

Our goal is to find a query strategy that offers minimal expected error on the
workload. Using the analytic error formula from Proposition 4, this can be
defined as a constrained optimization problem. One formulation of this problem
is stated below:

######  Problem 1 (Matrix Mechanism Optimization [25]).

Given an $m\times n$ workload matrix ${\bm{W}}$:

$\displaystyle\underset{\bm{A}}{\text{minimize}}$
$\displaystyle\left\lVert\bm{A}\right\rVert_{\mathcal{K}}^{2}\left\lVert{\bm{W}}\bm{A}^{+}\right\rVert_{F}^{2}$
(1) subject to $\displaystyle{\bm{W}}\bm{A}^{+}\bm{A}={\bm{W}}$

For a number of reasons, this optimization problem is difficult to solve
exactly: it has many variables, it is not convex, and both the objective
function and constraints involve $\bm{A}^{+}$, which can be slow to compute.
In addition, $\left\lVert{\bm{W}}\bm{A}^{+}\right\rVert_{F}^{2}$ has points of
discontinuity near the boundary of the constraint
${\bm{W}}\\!\bm{A}^{+}\\!\bm{A}={\bm{W}}$. This problem was originally
formulated as a rank-constrained semi-definite program [22], and, while
algorithms exist to find the global optimum, they require
$O(m^{4}(m^{4}+n^{4}))$ time, making it infeasible in practice.

Gradient-based numerical optimization techniques can be used to find locally
optimal solutions to Problem 1. These techniques begin by guessing a solution
$\bm{A}_{0}$ and then iteratively improve it using the gradient of the
objective function to guide the search. The process ends after a number of
iterations are performed, controlled by a stopping condition based on
improvement of the objective function. The constraints complicate the problem
further, but even if we ignore them, gradient-based optimization is slow, as
the cost of computing the objective function for general $\bm{A}$ is
$O(n^{3})$, e.g. requiring more than $6$ minutes for $n=8192$.

In the next sections, we provide algorithms for efficiently solving Problem 1.
We separately consider the two cases of Gaussian noise and Laplace noise, as
the required techniques are quite different.

###  4.2. Strategy Optimization with Gaussian Noise

While Problem 1 with $\mathcal{K}=\mathcal{G}$ is not convex in its current
form, it can be reformulated into an equivalent problem that is convex [25,
43]. The key idea is that the objective function can be expressed in terms of
$\bm{X}=\bm{A}^{T}\bm{A}$, since
$\left\|\bm{A}\right\|_{\mathcal{G}}^{2}=max(diag(\bm{A}^{T}\bm{A}))$ and
$\left\lVert{\bm{W}}\bm{A}^{+}\right\rVert_{F}^{2}=tr[(\bm{A}^{T}\bm{A})^{+}({\bm{W}}^{T}{\bm{W}})]$.
This allows us to optimize $\bm{X}$ instead of $\bm{A}$, and then we can
recover $\bm{A}$ by performing Cholesky decomposition on $\bm{X}$. Remarkably,
the resulting problem is convex with respect to $\bm{X}$.

######  Definition 12 (Convex Reformulation [43]).

Given a workload matrix ${\bm{W}}$ of rank $n$, let $\mbox{\sf
OPT}_{0}({\bm{W}})=\bm{A}$ where $\bm{A}^{T}\bm{A}$ is a Cholesky
decomposition of $\bm{X}^{*}$ and:

$\displaystyle\bm{X}^{*}=$ $\displaystyle\underset{\bm{X}}{\text{minimize}}$
$\displaystyle tr[\bm{X}^{-1}({\bm{W}}^{T}{\bm{W}})]$ (2) subject to
$\displaystyle diag(\bm{X})=\mathbf{1}$ $\displaystyle\bm{X}\succ 0$

While the above problem is convex, it is still nontrivial to solve due to the
dependence on the matrix inverse and the constraint $\bm{X}\succ 0$ ($\bm{X}$
is positive definite). The equality constraint $diag(\bm{X})=1$ and
corresponding optimization variables $diag(\bm{X})$ can easily be eliminated
from the problem since they must always equal $1$. Additionally $\bm{X}$ must
be a symmetric matrix so we can optimize over the lower triangular entries,
essentially reducing the number of optimization variables by a factor of two.
The full details of a conjugate gradient algorithm for solving this problem
are available in [43]. The per-iteration runtime of their “COA” algorithm is
$O(n^{3})$, and it typically requires about 50 iterations to converge. In
practice it is able to scale up to about $n\approx 10^{4}$. The algorithm
works well in practice when $n$ is small, but is not particularly robust for
some workloads when $n$ is larger, which we observe empirically in Section 11
(e.g., COA on prefix workload in Table 5).

We thus design our own algorithm to solve the same optimization problem, which
is based on the same principles as the COA technique, but is more robust in
practice. There are two key differences in our implementation. First, we
initialize the optimization intelligently by setting
$\bm{X}=\bm{P}\mathbf{\sqrt{\Lambda}}\bm{P}^{T}$ where
$\bm{P}\mathbf{\Lambda}\bm{P}^{T}$ is the eigen-decomposition of
${\bm{W}}^{T}{\bm{W}}$. This an approximation to the optimal strategy based on
the SVD bound [24], and acts as a very good initialization. Second, instead of
using the custom-designed conjugate gradient algorithm proposed in [43], we
simply use scipy.optimize, an off-the-shelf optimizer. We heuristically ignore
the constraint $\bm{X}\succ 0$ during optimization, using a large loss value
when it is not satisfied. This makes the problem unconstrained, and readily
solvable by scipy.optimize. The constraint is verified to hold at the end of
the optimization. These changes lead to more robust optimization that produce
strategies nearly matching the SVD bound, as we show experimentally in Section
11.

###  4.3. Strategy Optimization with Laplace noise

Unfortunately, the techniques used in the previous section do not apply to the
Laplace noise setting, as the sensitivity norm
$\left\|\bm{A}\right\|_{\mathcal{L}}$ cannot be expressed in terms of
$\bm{A}^{T}\bm{A}$. In this section, we describe an alternate approach to
approximately solve Problem 1: parameterized strategies. Our key idea is to
judiciously restrict the search space of the optimization problem to simplify
the optimization while retaining expressivity of the search space. While our
approach does not necessarily produce a globally optimal solution to Problem
1, with good parameterizations it can still find state of the art strategies.
Below we describe the idea of parameterized strategies in its full generality.
Then we propose a specific parameterization that works well for a variety of
input workloads.

A parameterization is a function $\bm{A}(\bm{\theta})$ mapping a real-valued
parameter vector $\bm{\theta}$ to a strategy $\bm{A}$. Optimizing over a
parameterized strategy space can be performed by optimizing $\bm{\theta}$
rather than $\bm{A}$. In the extreme case, there may be one entry in
$\bm{\theta}$ for every entry in $\bm{A}$, but a more careful design of the
parameterization with fewer parameters and a smart mapping between the entries
of the parameter vector and the entries of the strategy matrix can lead to
more effective optimization. There are several design considerations for
setting up a good parameterization. First, the parameterization must be
expressive enough to encode high-quality strategies (this may depend on the
workload). Second, the parameterization should have structure that makes the
optimization problem more computationally tractable (such as eliminating
constraints). Third, the parameterization may encode domain expertise about
what a good strategy should look like, which could make it easier to find
high-quality local minima.

Several existing privacy mechanisms can be thought of as an instance of this
parameterization framework [33, 21, 9, 44, 23]. These mechanisms typically are
designed for a specific workload or workload class. For example, Qardaji et
al. consider the space of hierarchical strategies, which is parameterized by a
single parameter: the branching factor of the hierarchy, which is chosen to
minimize MSE on the workload of _all range queries_ [33]. Li et al. consider
the space of weighted hierarchical strategies (with fixed branching factor),
which is parameterized by vector of scaling factors for each query in the
hierarchy [21]. This method adapts to the input workload, but the hierarchical
parameterization only works well for workloads of range queries. Ding et al.
consider the space of marginal query strategies (and workloads of the same
form), the parameters can be thought of as a binary vector of length $2^{d}$
corresponding to which marginals should be included in the strategy [9]. These
parameters are optimized using a greedy heuristic. Li et al. consider the
space of strategies containing the eigenvectors of the workload gram matrix,
which is parameterized by a vector of scaling factors for each eigenvector
[23]. This parameterization works well when $\mathcal{K}=\mathcal{G}$, but not
when $\mathcal{K}=\mathcal{L}$, and it has recently been subsumed by the work
of Yuan et al. for $\mathcal{K}=\mathcal{G}$ [43].

We now present a new general-purpose parameterization, called $p$-Identity,
which handles these considerations without making strict assumptions about the
structure of the workload. It also out-performs all of the existing
parameterizations, even on the workloads for which they were designed. The
parameters of a p-Identity strategy are more naturally interpreted as a matrix
$\bm{\Theta}$ rather than a vector $\bm{\theta}$ so we instead use the
notation $\bm{A}(\bm{\Theta})$.

######  Definition 13 ($p$-Identity strategies).

Given a $p\times n$ matrix of non-negative values $\bm{\Theta}$, the
$p$-Identity strategy matrix $\bm{A}(\bm{\Theta})$ is defined as follows:

$A(\bm{\Theta})=\begin{bmatrix}\bm{I}\\\ \bm{\Theta}\end{bmatrix}\bm{D}$

where $\bm{I}$ is the identity matrix and
$\bm{D}=diag(1+\bm{1}^{T}\bm{\Theta})^{-1}$.

Above, the matrix $\bm{A}(\bm{\Theta})$ contains $n+p$ queries, including $n$
identity queries (i.e., queries that count the number of records in the
database for each domain element in $dom(R)$) and $p$ arbitrary queries that
depend on the parameters $\bm{\Theta}$. The diagonal matrix $\bm{D}$ is used
to re-weight the strategy so that $\left\|\bm{A}\right\|_{\mathcal{L}}=1$ and
each column of $\bm{A}$ has the same $L_{1}$ norm. This is an example of
domain knowledge incorporated into the parameterization, as it has been shown
that optimal strategies have uniform column norm [25]. 333If a strategy did
not, a query could be added to it without increasing sensitivity and addition
of that query would result in error less than or equal to that of the original
strategy.

######  Example 10.

For $p=2$ and $n=3$, we illustrate below how $\bm{A}(\bm{\Theta})$ is related
to its parameter matrix, $\bm{\Theta}$.

$\displaystyle\bm{\Theta}=\begin{bmatrix}1&2&3\\\ 1&1&1\end{bmatrix}$
$\displaystyle\bm{A}(\bm{\Theta})=\begin{bmatrix}0.33&0&0\\\ 0&0.25&0\\\
0&0&0.2\\\ 0.33&0.5&0.6\\\ 0.33&0.25&0.2\end{bmatrix}$

For this class of parameterized strategies, the resulting optimization problem
requires optimizing $\bm{\Theta}$ instead of $\bm{A}$ and is stated below; we
use $\mbox{\sf OPT}_{0}$ to denote the operator that solves this problem.

######  Definition 14 (parameterized optimization).

Given a workload matrix ${\bm{W}}$ and hyper-parameter $p$, let $\mbox{\sf
OPT}_{0}({\bm{W}})=\bm{A}(\bm{\Theta}^{*})$ where:

$\displaystyle\bm{\Theta}^{*}=\underset{\bm{\Theta}\in\mathbb{R}_{+}^{p\times
n}}{\text{argmin}}$
$\displaystyle\left\lVert{\bm{W}}\bm{A}(\bm{\Theta})^{+}\right\rVert_{F}^{2}$

This parameterization was carefully designed to simplify optimization. Because
$\bm{A}(\bm{\Theta})$ is full rank, the constraints are satisfied by
construction, and we can solve an unconstrained problem instead. Furthermore,
$\left\lVert\bm{A}(\bm{\Theta})\right\rVert_{1}=1$ for all $\bm{\Theta}$, so
that term can be removed from the objective. Additionally, due to the special
structure of $\bm{A}(\bm{\Theta})$ we can efficiently evaluate the objective
and its gradient in $O(pn^{2})$ time instead of $O(n^{3})$ time. For example,
when $n=8192$ it requires $>6$ minutes to evaluate the objective for a general
strategy $\bm{A}$, while it takes only $1.5$ seconds for a $p$-Identity
strategy (with $p=\frac{n}{16}$), which is a $240\times$ improvement. Despite
this imposed structure, $\bm{A}(\bm{\Theta})$ is still expressive enough to
encode high-quality strategies. Moreover, $p$ can always be tuned to balance
expressivity with efficiency.

###  4.4. Strategy Visualization

(a) Cholesky strategy (b) p-Identity strategy
(p=16) (c) Hierarchical strategy
(b=16)

Figure 1.  Visualization of three different strategies for answering the
workload of All Range queries on a domain of size $256$. Figures (a) and (b)
show strategies optimized by HDMM, and Figure (c) shows a hierarchical
strategy with branching factor $16$, a previously state-of-the-art strategy
for this workload [19, 33]. Each row is a query, and cells are color coded
according to their value in the strategy matrix. Figure (a) plots each query
as one very thin row, while Figure (b) and (c) only plot the non-trivial
queries as thicker rows. The diagonal queries are plotted separately as a
single row above the main plot, but it actually represents $256$ different
queries.

An interested reader may wonder what the optimized strategies actually look
like for some common workloads. In Figure 1, we plot three strategies designed
for the workload of All Range queries. Figures 1(a) and 1(b) show the
strategies produced by HDMM for Gaussian and Laplace noise, respectively, and
Figure 1(c) shows a previously state-of-the-art strategy for range queries.
These visualizations provide interesting insight into the nature of the
solution, and reveal why HDMM succeeds in reducing error.

In Figure 1(a), the strategy is an upper triangular square matrix; this is by
construction as it is obtained through a Cholesky decomposition.444There may
be equally good strategies that do not have this upper triangular structure,
but HDMM will always find one with this structure. In each row, weights are
largest near the main diagonal, and gradually decrease the further away it
gets. This can be observed from the smooth transition from yellow to green to
blue in Figure 1(a).

In Figure 1(b), the strategy was optimized with $p=16$, but only $13$ non-zero
queries were found. Each query has varying width, ranging between about $16$
and $64$, and has greatest weight towards the middle of the query, with
gradually decreasing weights away from the center. In most cases, the queries
overlap with half of the two neighboring queries. The weights on the identity
queries are approximately uniform throughout at around $0.5$, with higher
weights near the edges. There is a very natural reason why this structure
works well for range queries. Summing up adjacent queries leads to a bigger
query which looks approximately like a range query. It will have a long
uniform center, and decaying weights on the edges. These decaying weights can
be increased to match the uniform center by drawing on the answers from the
uniform queries. Thus, any range query can be answered by summing up a
relatively small number of strategy query answers.

This same intuition was used in the derivation of the hierarchical strategy in
Figure 1(c). However, this strategy answers some range queries more
effectively than others. It struggles for queries that require summing up many
identity queries. For example, it can answer the range $[0,16)$ using one
strategy query, but it requires summing up $16$ strategy queries to answer the
range $[8,24)$.

##  5\. Implicit representations for conjunctive query sets

The optimization methods we described in the previous section work well for
small and modest domain sizes; we were able to run them on domains as large as
$n=8192$ (see Section 11 for a scalability analysis). However, these methods
are fundamentally limited by the need to represent the workload and strategy
explicitly in matrix form. It requires $0.5$ gigabytes just to store a square
matrix of size $8192$ using $4$ byte floats, and it is time consuming to
perform nontrivial matrix operations on objects of this size. This limitation
is not unique to our mechanism, but is shared by _all possible_ methods that
optimize explicitly represented workload matrices.

To overcome this scalability limitation, we propose _implicit query matrices_
, which exploit structure in conjunctive query sets and offer a far more
concise representation than materialized explicit matrices, while still being
able to encode query sets containing an arbitrary collection of conjunctive
queries. These representations are lossless; they save space by avoiding
significant redundancy, rather than making approximations. As we will show
later in this section, many important matrix operations can be performed
efficiently using the implicit representation. This property of the
representation will be essential for solving the strategy optimization problem
efficiently on large domains.

###  5.1. Implicitly vectorized conjunctions

Consider a predicate defined on a single attribute, $A$, where
$|dom(A)|=n_{A}$. This predicate, $\phi_{A}$, can be vectorized with respect
to just the domain of $A$ (and not the full domain of all attributes)
similarly to Definition 5. When a predicate is formed from the conjunction of
such single-attribute predicates, its vectorized form has a concise implicit
representation in terms of the outer product between vectors.

######  Definition 15 (Outer product).

The outer product between two vectors $\bm{q}_{A}$ and $\bm{q}_{B}$, denoted
$\bm{q}_{A}\otimes\bm{q}_{B}$, is a vector indexed by pairs $t=(t_{A},t_{B})$
such that:
$(\bm{q}_{A}\otimes\bm{q}_{B})(t)=\bm{q}_{A}(t_{A})\bm{q}_{B}(t_{B})$.

The outer product is useful for representing conjunctions in vector form.

######  Theorem 1 (Implicit vectorization).

The vector representation of the conjunction $\phi=\phi_{A}\wedge\phi_{B}$ is
$vec(\phi)=vec(\phi_{A})\otimes vec(\phi_{B})$.

To see why this is true, observe that
$\phi(t)=\phi_{A}(t_{A})\phi_{B}(t_{B})$, since multiplication and logical AND
are equivalent for binary inputs, and this is the exact equation that defines
the outer product for vectors. While the explicit representation of
$vec(\phi)$ has size $n_{A}\cdot n_{B}$, the implicit representation,
$vec(\phi_{A})\otimes vec(\phi_{B})$, requires storing only $vec(\phi_{A})$
and $vec(\phi_{B})$, which has size $n_{A}+n_{B}$.

######  Example 11.

Recall that the workload $\mathcal{W}_{\mbox{\tiny SF1}}$ consists of 4151
queries, each defined on a data vector of size $500480$. Since explicitly
vectorized queries are the same size as the domain, the size of the explicit
workload matrix is $4151\times$500480$$, or $8.3$GB. Using the implicit
representation, each query can be encoded using $2+2+64+115+17=200$ values,
for a total of 3.3MB. For $\mathcal{W}_{\mbox{\tiny SF1+}}$, which consists of
$215852$ queries on a data vector of size $25524480$, the explicit workload
matrix would require 22TB of storage. In contrast, the implicit vector
representation would require 200MB.

###  5.2. Implicitly vectorized products

Product workloads (as in Definition 2 can be encoded even more efficiently
using the Kronecker product.

######  Definition 16 (Kronecker product).

For two matrices $\bm{A}\in\mathbb{R}^{m_{A}\times n_{A}}$ and
$\bm{B}\in\mathbb{R}^{m_{B}\times n_{B}}$, their Kronecker product is
$\bm{A}\otimes\bm{B}\in\mathbb{R}^{m_{A}m_{B}\times n_{A}n_{B}}$, where rows
are indexed by pairs $q=(q_{A},q_{B})$ and columns are indexed by pairs
$t=(t_{A},t_{B})$ such that:

$(\bm{A}\otimes\bm{B})(q,t)=\bm{A}(q_{A},t_{A})\bm{B}(q_{B},t_{B})$

The Kronecker product is a generalization of the outer product, where each row
is an outer product between a pair of rows from $\bm{A}$ and $\bm{B}$. Thus,
the same symbol $\otimes$ is used for both operations.

######  Theorem 2 (Implicit vectorization).

Let $\Phi_{A}$ and $\Phi_{B}$ be two predicate sets defined on attributes $A$
and $B$ respectively. Then $vec(\Phi_{A}\times\Phi_{B})=vec(\Phi_{A})\otimes
vec(\Phi_{B})$.

The proof of this claim follows immediately from Definition 2, since it
contains a cartesian product of conjunctions, and each conjunction is an outer
product. We implicitly represent a product workload in matrix form by storing
the factors of the Kronecker product ($\bm{A}=vec(\Phi_{A})$ and
$\bm{B}=vec(\Phi_{B})$). This requires $m_{A}n_{A}+m_{B}n_{B}$ space, rather
than $m_{A}m_{B}n_{A}n_{B}$ space, which is required by the explicit
representation, and $m_{A}m_{B}(n_{A}+n_{B})$ for the implicit representation
using a list of outer products. Thus, the savings can be quite substantial.
While Theorem 2 assumes the predicates are conjunctions, we show in Appendix A
that disjunctions can be handled in a similar manner.

###  5.3. Workload encoding algorithm

Given as input a logical workload $\mathcal{W}$ (as in Definition 3), the
ImpVec algorithm produces an implicitly represented workload matrix with the
following form:

$\mathbb{W}=\begin{bmatrix}w_{1}\mathbb{W}_{1}\\\ \vdots\\\
w_{k}\mathbb{W}_{k}\end{bmatrix}=\begin{bmatrix}w_{1}({\bm{W}}_{1}^{(1)}\otimes&\dots&\otimes{\bm{W}}_{d}^{(1)})\\\
\vdots&\ddots&\vdots\\\
w_{k}({\bm{W}}_{1}^{(k)}\otimes&\dots&\otimes{\bm{W}}_{d}^{(k)})\end{bmatrix}$
(3)

Here stacking sub-workloads is analogous to union and in formulas we will
sometimes abuse notation and write an implicit union-of-products workload as
$\mathbb{W}_{[k]}=w_{1}\mathbb{W}_{1}{+}\dots+w_{k}\mathbb{W}_{k}$. We use
blackboard bold font to distinguish an implicitly represented workload
$\mathbb{W}$ from an explicitly represented workload ${\bm{W}}$.

Algorithm 1: ImpVec
---
Input: Workload $\mathcal{W}=\\{q_{1}\dots q_{k}\\}$ and weights $w_{1}\dots
w_{k}$
Output: Implicit workload $\mathbb{W}$
1. | For each product $q_{i}\in\mathcal{W}$: $q_{i}=\Phi^{(i)}_{A_{1}}\times\dots\times\Phi^{(i)}_{A_{d}}$
2. |  For each $j\in[1..d]$
3. |  compute ${\bm{W}}^{(i)}_{j}=vec(\Phi^{(i)}_{A_{j}})$
4. |  Let $\mathbb{W}_{i}={\bm{W}}_{1}^{(i)}\otimes\dots\otimes{\bm{W}}_{d}^{(i)}$
5. |  Return: $w_{1}\mathbb{W}_{1}{+}\dots+w_{k}\mathbb{W}_{k}$

Note that line 3 of the ImpVec algorithm is _explicit_ vectorization, as in
Definition 5, of a set of predicates on a single attribute.

######  Example 12.

Recall from Example 9 that the $215852$ queries of $\mathcal{W}_{\mbox{\tiny
SF1+}}$can be represented as $k=8032$ products. If $\mathcal{W}_{\mbox{\tiny
SF1+}}$ is represented in this factored form, the ImpVec algorithm returns a
smaller implicit representation, reducing the 200MB (from Example 11) to 50
MB. If the workloads are presented in their manually factored format of $k=32$
products, the implicit representation of $\mathcal{W}_{\mbox{\tiny SF1}}^{*}$
requires only 335KB, and $\mathcal{W}_{\mbox{\tiny SF1+}}^{*}$ only 687KB.

###  5.4. Operations on vectorized objects

Reducing the size of the workload representation is only useful if critical
computations can be performed without expanding them to their explicit
representations. Standard properties of the Kronecker product [36] accelerate
strategy selection and reconstruction.

The following properties allow us to perform useful operations on Kronecker
products without materializing their full matrices.

######  Proposition 6 (Kronecker identities).

Kronecker products satisfy the following identities [36]:
Transpose: $(\bm{A}\otimes\bm{B})^{T}=\bm{A}^{T}\otimes\bm{B}^{T}$ Pseudo
Inverse: $(\bm{A}\otimes\bm{B})^{+}=\bm{A}^{+}\otimes\bm{B}^{+}$
Associativity:
$(\bm{A}\otimes\bm{B})\otimes\bm{C}=\bm{A}\otimes(\bm{B}\otimes\bm{C})$ Mixed
Product:
$(\bm{A}\otimes\bm{B})(\bm{C}\otimes\bm{D})=(\bm{A}\bm{C})\otimes(\bm{B}\bm{D})$

In addition to the standard properties of Kronecker product above, we can
prove additional properties about them which are useful for our privacy
mechanism.

######  Theorem 3 (Norm of a Kronecker product).

The following matrix norms decompose over the factors of the Kronecker product
$\bm{A}_{1}\otimes\dots\otimes\bm{A}_{d}$.

$\displaystyle\left\|\bm{A}_{1}\otimes\dots\otimes\bm{A}_{d}\right\|_{\mathcal{L}}$
$\displaystyle=\prod_{i=1}^{d}\left\|\bm{A}_{i}\right\|_{\mathcal{L}}$
$\displaystyle\left\|\bm{A}_{1}\otimes\dots\otimes\bm{A}_{d}\right\|_{\mathcal{G}}$
$\displaystyle=\prod_{i=1}^{d}\left\|\bm{A}_{i}\right\|_{\mathcal{G}}$
$\displaystyle\left\lVert\bm{A}_{1}\otimes\dots\otimes\bm{A}_{d}\right\rVert_{F}$
$\displaystyle=\prod_{i=1}^{d}\left\lVert\bm{A}_{i}\right\rVert_{F}$

The proof is largely routine algabraic manipulation and is given in the
appendix. We will later use these identities to efficiently evaluate TSE for
workloads and strategies built with Kronecker products.

##  6\. Optimizing conjunctive query workloads with conjunctive query
strategies

We now turn our attention to optimizing implicitly-represented conjunctive
query workloads. We assume the workload is a _union of Kronecker products_ ,
and that it takes the form shown in Equation 3 (restated below):

$\mathbb{W}=\begin{bmatrix}w_{1}\mathbb{W}_{1}\\\ \vdots\\\
w_{k}\mathbb{W}_{k}\end{bmatrix}=\begin{bmatrix}w_{1}({\bm{W}}_{1}^{(1)}\otimes&\dots&\otimes{\bm{W}}_{d}^{(1)})\\\
\vdots&\ddots&\vdots\\\
w_{k}({\bm{W}}_{1}^{(k)}\otimes&\dots&\otimes{\bm{W}}_{d}^{(k)})\end{bmatrix}$
(3)

For notational convenience, we will often write
$\mathbb{W}=w_{1}\mathbb{W}_{1}+\dots+w_{k}\mathbb{W}_{k}$, where
$\mathbb{W}_{i}={\bm{W}}_{1}^{(i)}\otimes\dots\otimes{\bm{W}}_{d}^{(i)}$ and
‘$+$’ serves the role of stacking matrices. The methods described in this
section will exploit the special structure of this class of workloads to scale
more effectively than techniques described in Section 4.

###  6.1. Optimizing product workloads.

We begin by considering a special case of the workload in Equation 3 that is a
_single Kronecker product_. We describe an elegant approach to optimize
workloads of this form by analyzing the optimization problem and utilizing
properties of Kronecker products.

######  Definition 17 ($\mbox{\sf OPT}_{\otimes}$).

Given a Kronecker product workload
$\mathbb{W}={\bm{W}}_{1}\otimes\dots\otimes{\bm{W}}_{d}$ and an optimization
oracle $\mbox{\sf OPT}_{0}$, let $\mbox{\sf
OPT}_{\otimes}(\mathbb{W})=\bm{A}_{1}\otimes\dots\otimes\bm{A}_{d}$ where
$\bm{A}_{i}=\mbox{\sf OPT}_{0}({\bm{W}}_{i})$.

Above, $\mbox{\sf OPT}_{0}$ may be any strategy optimization routine that
consumes an explicitly represented workload and returns a strategy matrix,
such as the techniques discussed in the previous section. $\mbox{\sf
OPT}_{\otimes}$ is very efficient: it simply requires solving $d$ small
subproblems rather than one large problem. This decomposition of the objective
function has a well-founded theoretical justifiication. Namely, if we restrict
the solution space to a (single) Kronecker product strategy of the form
$\mathbb{A}=\bm{A}_{1}\otimes\dots\otimes\bm{A}_{d}$, then the error of the
workload under $\mathbb{A}$ decomposes over the factors of the Kronecker
products as shown in the following theorem:

######  Theorem 4 (Error decomposition).

Given a workload $\mathbb{W}={\bm{W}}_{1}\otimes\dots\otimes{\bm{W}}_{d}$ and
a strategy $\mathbb{A}=\bm{A}_{1}\otimes\dots\otimes\bm{A}_{d}$, the error is:

$\left\|\mathbb{A}\right\|_{\mathcal{K}}^{2}\left\lVert\mathbb{W}\mathbb{A}^{+}\right\rVert_{F}^{2}=\prod_{i=1}^{d}\left\|\bm{A}_{i}\right\|_{\mathcal{K}}^{2}\left\lVert{\bm{W}}_{i}\bm{A}_{i}^{+}\right\rVert_{F}^{2}$

The overall error is minimized when $\bm{A}_{i}$ optimizes ${\bm{W}}_{i}$ for
each $i$, thus it makes sense to optimize each $\bm{A}_{i}$ separately. If we
expect the optimal strategy to be a single Kronecker product, then this
approach seems quite appealing. However it is possible that there exists a
strategy that is _not_ a single Kronecker product that offers lower error than
the best Kronecker product strategy. The following theorem shows that this is
not the case, and gives further justification for the above method, showing
that the SVD bound also decomposes over the factors of the Kronecker product.

######  Theorem 5 (SVD bound decomposition).

Given a workload $\mathbb{W}={\bm{W}}_{1}\otimes\dots\otimes{\bm{W}}_{d}$, the
SVD bound is:

$SVDB(\mathbb{W})=\prod_{i=1}^{d}SVDB({\bm{W}}_{i})$

If there exist strategies $\bm{A}_{1},\dots,\bm{A}_{d}$ that achieve the SVD
bound for ${\bm{W}}_{1},\dots,{\bm{W}}_{d}$ and we can find them, then we can
construct a Kronecker product strategy
$\mathbb{A}=\bm{A}_{1}\otimes\dots\otimes\bm{A}_{d}$ that achieves the SVD
bound for $\mathbb{W}={\bm{W}}_{1}\otimes\dots\otimes{\bm{W}}_{d}$. Since no
other strategy can have lower error than the SVD bound, this implies that _the
optimal strategy is a Kronecker product_. This gives excellent justification
for optimizing over the space of Kronecker products, especially when the SVD
bound is tight.

###  6.2. Optimizing union-of-product workloads.

The approach just described is principled and effective when the workload is a
single Kronecker product. We now turn our attention to the more general case
where the workload is a union of Kronecker products. Here, the right approach
is less clear. We define three approaches for optimizing implicit workloads in
the form of Equation 3. Each approach restricts the strategy to a different
region of the full strategy space for which optimization is tractable. The
first computes a strategy consisting of a single product; it generalizes
$\mbox{\sf OPT}_{\otimes}$. The second, $\mbox{\sf OPT}_{+}$, can generate
strategies consisting of unions of products. The third, $\mbox{\sf
OPT}_{\mathsf{M}}$, generates a strategy of weighted marginals. The best
approach to use will generally depend on the workload, and we will provide
some practical guidance and intuition to understand the situations in which
each method works best.

Single-product output strategy For weighted union of product workloads, if we
restrict the optimization problem to a single product strategy, then the
objective function decomposes as follows:

######  Theorem 6.

Given workload $\mathbb{W}=w_{1}\mathbb{W}_{1}{+}\dots+w_{k}\mathbb{W}_{k}$
and strategy $\mathbb{A}=\bm{A}_{1}\otimes\dots\otimes\bm{A}_{d}$, workload
error is:

$\begin{split}\left\|\mathbb{A}\right\|_{\mathcal{K}}^{2}\left\lVert\mathbb{W}\mathbb{A}^{+}\right\rVert_{F}^{2}&=\left\|\mathbb{A}\right\|_{\mathcal{K}}^{2}\sum_{j=1}^{k}w_{j}^{2}\left\lVert\mathbb{W}_{j}\mathbb{A}^{+}\right\rVert_{F}^{2}\\\
&=\sum_{j=1}^{k}w_{j}^{2}\prod_{i=1}^{d}\left\|\bm{A}_{i}\right\|_{\mathcal{K}}^{2}\left\lVert{\bm{W}}^{(j)}_{i}\bm{A}_{i}^{+}\right\rVert_{F}^{2}\end{split}$
(4)

This leads to the following optimization problem:

######  Definition 18 (Generalized $\mbox{\sf OPT}_{\otimes}$).

Given a workload $\mathbb{W}=w_{1}\mathbb{W}_{1}{+}\dots+w_{k}\mathbb{W}_{k}$,
let $\mbox{\sf
OPT}_{\otimes}(\mathbb{W})=\bm{A}_{1}\otimes\dots\otimes\bm{A}_{d}$ where:

$\displaystyle(\bm{A}_{1},\dots,\bm{A}_{d})$
$\displaystyle=\underset{\bm{A}_{1},\dots,\bm{A}_{d}}{\text{minimize}}$
$\displaystyle\sum_{j=1}^{k}w_{j}^{2}\prod_{i=1}^{d}\left\|\bm{A}_{i}\right\|_{\mathcal{K}}^{2}\left\lVert{\bm{W}}^{(j)}_{i}\bm{A}_{i}^{+}\right\rVert_{F}^{2}$

When $k=1$, the solution to Problem 18 is given in Definition 17, so we use
matching notation and allow the $\mbox{\sf OPT}_{\otimes}$ operator to accept
a single product or a union of products.

We can solve this problem efficiently by building on the optimization oracles
designed in the previous section. In particular, suppose we have a black box
optimization oracle $\mbox{\sf OPT}_{0}({\bm{W}})$ that accepts an explicitly
represented workload and gives back an explicitly represented strategy with
low (ideally minimal) error on that workload. Then we use a block method that
cyclically optimizes $\bm{A}_{1},\dots,\bm{A}_{d}$ until convergence. We begin
by initializing $\bm{A}_{i}=\bm{I}$ for all $i$. We then optimize one
$\bm{A}_{i}$ at a time, fixing the other $\bm{A}_{i^{\prime}}\neq\bm{A}_{i}$
using $\mbox{\sf OPT}_{0}$ on a carefully constructed surrogate workload
$\hat{{\bm{W}}}$ (Equation 5) that has the property that the error of any
strategy $\bm{A}_{i}$ on $\hat{{\bm{W}}}$ is the same as the error of
$\mathbb{A}$ on $\mathbb{W}$. Hence, the correct objective is being minimized.

$\displaystyle\hat{{\bm{W}}}_{i}=\begin{bmatrix}c_{1}{\bm{W}}_{i}^{(1)}\\\
\vdots\\\ c_{k}{\bm{W}}_{i}^{(k)}\end{bmatrix}$ $\displaystyle
c_{j}=w_{j}\prod_{i^{\prime}\neq
i}\left\|\bm{A}_{i^{\prime}}\right\|_{\mathcal{K}}\left\lVert{\bm{W}}_{i^{\prime}}^{(j)}\bm{A}_{i^{\prime}}^{+}\right\rVert_{F}$
(5)

The cost of running this optimization procedure is determined by the cost of
computing $\hat{{\bm{W}}}_{i}^{T}\hat{{\bm{W}}}_{i}$ and the cost of
optimizing it, which takes $O(n_{i}^{2}(p_{i}+k))$ and
$O(n_{i}^{2}p_{i}\cdot\text{\sc\\#iter})$ time respectively (assuming each
$({\bm{W}}^{T}{\bm{W}})_{i}^{(j)}$ has been precomputed). As before, this
method scales to arbitrarily large domains as long as the domain size of the
sub-problems allows $\mbox{\sf OPT}_{0}$ to be efficient.

Union-of-products output strategy For certain workloads, restricting to
solutions consisting of a single product, as $\mbox{\sf OPT}_{\otimes}$ does,
excludes good strategies, as shown in Example 13.

######  Example 13.

Consider the workload $\mathbb{W}=\mathbb{W}_{1}+\mathbb{W}_{2}$ where
$\mathbb{W}_{1}=\bm{P}\otimes\bm{T}$ and $\mathbb{W}_{2}=\bm{T}\otimes\bm{P}$
on a 2-dimensional domain of size $100\times 100$. Running $\mbox{\sf
OPT}_{\otimes}$ on this workload leads to an optimized strategy of the form
$\mathbb{A}=\bm{A}_{1}\otimes\bm{A}_{2}$. The expected error of this strategy
is $33385$, which is much higher than it should be for such a simple workload.
The poor expected error can be explained by the fact that to support the
workload, both $\bm{A}_{1}$ and $\bm{A}_{2}$ have to be full rank. This means
$\mathbb{A}$ has to include at least $100^{2}$ queries, even though
$\mathbb{W}$ only contains $200$ queries.

A better alternative would be to optimize $\mathbb{W}_{1}$ and
$\mathbb{W}_{2}$ separately using $\mbox{\sf OPT}_{\otimes}$. Doing this we
would end up with a strategy $\mathbb{A}=\mathbb{A}_{1}+\mathbb{A}_{2}$, where
$\mathbb{A}_{1}$ optimizes $\mathbb{W}_{1}$ and $\mathbb{A}_{2}$ optimizes
$\mathbb{W}_{2}$. The resulting strategy is much smaller because $\mbox{\sf
OPT}_{\otimes}(\mathbb{A}_{1})=\mbox{\sf OPT}_{0}(\bm{P})\otimes\mbox{\sf
OPT}_{0}(\bm{T})$, and $\mbox{\sf OPT}_{0}(\bm{T})=\bm{T}$. In fact, it only
contains $212$ queries and attains an expected error of $14252$, which is a
$2.34\times$ improvement.

Based on this example, we would like a principled approach to optimize over
the space of strategies that are a union of Kronecker products. Unfortunately,
computing the workload error exactly for a strategy of this form is
intractable, as the pseudo inverse may not be a union of Kronecker products.
This makes optimization over this space of strategies challenging. We thus
propose the following heuristic optimization routine inspired by Example 13.
This optimization routine individually optimizes each subworkload
$\mathbb{W}_{j}$ using $\mbox{\sf OPT}_{\otimes}$, and then combines the
strategies all together to form a single strategy. It simply requries calling
$\mbox{\sf OPT}_{\otimes}$ a number of times and computing appropriate weights
for each optimized strategy.

######  Definition 19 ($\mbox{\sf OPT}_{+}$).

Given a workload $\mathbb{W}=w_{1}\mathbb{W}_{1}+\dots+w_{k}\mathbb{W}_{k}$,
let $\mbox{\sf
OPT}_{+}(\mathbb{W})=c_{1}\mathbb{A}_{1}+\dots+c_{k}\mathbb{A}_{k}$ where
$\mathbb{A}_{j}=\mbox{\sf OPT}_{\otimes}(\mathbb{W}_{j})$ and

$c_{j}\propto\frac{1}{\left\|\mathbb{A}_{j}\right\|_{\mathcal{K}}}\begin{cases}\sqrt[3]{2E_{j}}&\text{
if }\mathcal{K}=\mathcal{L}\\\ \sqrt[4]{E_{j}}&\text{ if
}\mathcal{K}=\mathcal{G}\end{cases}$

for
$E_{j}=w_{j}^{2}\left\|\mathbb{A}_{j}\right\|_{\mathcal{K}}^{2}\left\lVert\mathbb{W}_{j}\mathbb{A}_{j}^{+}\right\rVert_{F}^{2}$.

Above, we assume that $\mathbb{A}_{j}$ will be used to answer
$\mathbb{W}_{j}$, and $c_{j}$ is the weight on $\mathbb{A}_{j}$: it
corresponds the portion of the privacy budget that will be spent to answer
those queries. It is chosen to minimize total workload error. Specifically, if
we allocate $c_{j}$ budget to answer $\mathbb{A}_{j}$, then the error will be
$E_{j}/c_{j}^{2}$. Thus, the choice of $c_{j}$ above is based on minimizing
$\sum E_{j}/c_{j}^{2}$ subject to the constraint $\sum|c_{j}|=1$ (for Laplace
noise) and $\sum c_{j}^{2}=1$ (for Gaussian noise). We solve this minimization
problem exactly, in closed from, using the method of Lagrange multipliers.

We remark that in the definition above, $\mathbb{W}$ is split up into $k$ sub-
workloads $\mathbb{W}_{1},\dots,\mathbb{W}_{k}$. Each subworkload
$\mathbb{W}_{j}$ is assumed to be a single Kronecker product, but the
optimization routine is still well defined even if $\mathbb{W}_{j}$ is a union
of Kronecker products. This opens up a nice opportunity: to group the
subworkloads into clusters which will be optimized together with $\mbox{\sf
OPT}_{\otimes}$. Intuitively, if two subworkloads are similar, it may make
sense to group them together to optimize collectively. We do not provide an
automated way to group subworkloads in this paper. This is a hard problem in
general, and is out of scope for this paper. A domain expert can work out good
clusterings on a case-by-case basis, or they can settle for the default
clustering (one Kronecker product per cluster).

##  7\. Implicit representations for marginal query sets

In the previous section we described $\mbox{\sf OPT}_{\otimes}$ and $\mbox{\sf
OPT}_{+}$, two methods for optimizing implicitly represented conjunctive query
workloads. These methods differ primarily in the space of strategies they
search over. Our final optimization method, $\mbox{\sf OPT}_{\mathsf{M}}$,
optimizes over the space of marginal query matrices, and offers a preferable
alternative to $\mbox{\sf OPT}_{\otimes}$ and $\mbox{\sf OPT}_{+}$ in some
settings.

In this section, we propose an implicit representation for marginal query sets
that is even more compact than our other representation for general
conjunctive query sets. We further show that these matrices can be operated on
efficiently, allowing us to solve the strategy optimization problem for large
multi-dimensional domains.

A marginal query matrix is a special case of a union of Kronecker products,
where each Kronecker product encodes the queries to compute a single marginal
(i.e., all the factors are either Identity or Total). First note that a
marginal on a $d$-dimensional domain can be specified by a subset of elements
of $\\{1,\dots,d\\}$. Hence, there are a total of $2^{d}$ possible marginals,
and each one can be specified by an element of the set
$[2^{d}]=\\{0,\dots,2^{d}-1\\}$. The most natural correspondence between these
integers and the associated marginals is based on the binary representation of
the integer. The query set required to compute the $a^{th}$ marginal would be
represented by $\bm{Q}_{1}\otimes\dots\otimes\bm{Q}_{d}$ where
$\bm{Q}_{i}=\bm{I}$ if the $i^{th}$ bit of the binary representation of $a$ is
$1$ and $\bm{Q}_{i}=\bm{T}$ otherwise. A collection of weighted marginals can
thus be represented as a vector $\bm{u}$ containing a weight for each
marginal. We refer to this marginal query matrix as $\mathbb{M}(\bm{u})$,
which is defined below.

######  Definition 20 (Marginal query matrix).

A marginal query matrix $\mathbb{M}(\bm{u})$ is defined by a vector of weights
$\bm{u}\in\mathbb{R}^{2^{d}}$ and is a special case of the query matrix shown
in Equation 3 where $k=2^{d}$, $w_{a+1}=\bm{u}(a)$, and

${\bm{W}}^{(a+1)}_{i}=\begin{cases}\bm{I}&a_{i}=1\\\
\bm{T}&a_{i}=0\end{cases}$

where $a\in\\{0,\dots,2^{d}-1\\}$,$i\in\\{1,\dots,d\\}$,and $a_{i}$ is the
$i^{th}$ bit of the binary representation of $a$.

For a marginal query matrix, the weight $\bm{u}(a)$ can be interpreted as the
relative _importance_ of the $a^{th}$ marginal. As shown in Proposition 7, it
is particularly simple to reason about the sensitivity of a marginal query
matrix.

######  Proposition 7.

The sensitivity of a marginal query matrix $\mathbb{M}(\bm{u})$ is:

$\displaystyle\left\|\mathbb{M}(\bm{u})\right\|_{\mathcal{L}}=\left\lVert\bm{u}\right\rVert_{1}$
$\displaystyle\left\|\mathbb{M}(\bm{u})\right\|_{\mathcal{G}}=\left\lVert\bm{u}\right\rVert_{2}$

Moving forward, it is convenient to work with the Gram matrix representation
of the marginal query matrix instead. As shown below, there is a simple
correspondence between the two:

######  Proposition 8 (Marginal Gram matrix).

Let $\mathbb{Q}=\mathbb{M}(\bm{u})$ be a marginal query matrix. Then the
correpsonding marginal Gram matrix is
$\mathbb{Q}^{T}\mathbb{Q}=\mathbb{G}(\bm{u}^{2})$, where:

$\mathbb{G}(\bm{v})=\sum_{a=0}^{2^{d}-1}\bm{v}(a)\mathbb{H}(a),\;\;\;\;\;\;\;\;\;\;\mathbb{H}(a)=\bigotimes_{i=1}^{d}[\bm{1}(a_{i}=0)+\bm{I}(a_{i}=1)]$

In the proposition above, the term $\bm{1}(a_{i}=0)+\bm{I}(a_{i}=1)$ is
shorthand notation for $\begin{cases}\bm{1}&a_{i}=0\\\
\bm{I}&a_{i}=1\end{cases}$. $\bm{1}$ and $\bm{I}$ are both $n_{i}\times n_{i}$
matrices, corresponding to the matrix of all ones, and the identity matrix
respectively. We will use this notation frequently in the section, so it is
important to understand the exact meaning. Another important object that will
appear repeatedly throughout this section is the so-called characteristic
vector 555this is not to be confused with the eigenvector, which is defined
below:

######  Definition 21 (Characteristic Vector).

The characteristic vector $\bm{c}\in\mathbb{R}^{2^{d}}$ is defined so that
each entry $\bm{c}(a)$ equals the number of entries in the $(\neg a)^{th}$
marginal

$\bm{c}(a)=\prod_{i=1}^{d}n_{i}(a_{i}=0)+1(a_{i}=1)$

The term $\neg a$ in the definition above is the bitwise negation of $a$, and
it is obtained by flipping each of the $d$ bits of the integer $a$. We will
rely heavily on this type of bitwise manipulation in this section to reason
about the behavior of marginal Gram matrices.

######  Example 14.

For a 3 dimensional domain, the query matrix for a 2 way marginal can be
expressed as $\mathbb{Q}=\bm{I}\otimes\bm{T}\otimes\bm{I}$, and
$\mathbb{Q}^{T}\mathbb{Q}=\bm{I}\otimes\bm{1}\otimes\bm{I}=\mathbb{H}(101_{2})=\mathbb{H}(5)$.

Now that we have introduced the necessary notation for marginal query and Gram
matrices, we are ready to show how to perform important matrix operations
while respecting the implicit representation. We begin with Theorem 7, which
shows that marginal Gram matrices interact nicely under matrix multiplication.

######  Theorem 7 (Multiplication of Marginal Gram Matrices).

For any $a,b\in[2^{d}]$,

$\mathbb{H}(a)\mathbb{H}(b)=\bm{c}(a|b)\mathbb{H}(a\&b)$

where $a|b$ denotes “bitwise or”, $a\&b$ denotes “bitwise and”, and $\bm{c}$
is the characteristic vector. Moreover, for any
$\bm{u},\bm{v}\in\mathbb{R}^{2^{d}}$

$\mathbb{G}(\bm{u})\mathbb{G}(\bm{v})=\mathbb{G}(\bm{X}(\bm{u})\bm{v})$

where $\bm{X}(\bm{u})$ is a $2^{d}\times 2^{d}$ triangular matrix with entries
$\bm{X}(\bm{u})(k,b)=\sum_{a:a\&b=k}\bm{u}(a)\bm{c}(a|b)$.

Theorem 7 allows us to efficiently multiply two matrices while maintaining the
compact implicit representation. Additionally, it follows immediately from the
proof of Theorem 7 that
$\mathbb{G}(\bm{u})\mathbb{G}(\bm{v})=\mathbb{G}(\bm{v})\mathbb{G}(\bm{u})$ —
i.e., matrix multiplication is commutative. We can apply Theorem 7 to
efficiently find the inverse or generalized inverse of $\mathbb{G}(\bm{u})$ as
well.

######  Theorem 8 (Inverse of Marginal Gram Matrices).

Let $\bm{X}(\bm{u})$ be the matrix defined in Theorem 7. If $\bm{X}(\bm{u})$
is invertible, then $\mathbb{G}(\bm{u})$ is invertible with inverse:

$\mathbb{G}^{-1}(\bm{u})=\mathbb{G}(\bm{X}^{-1}(\bm{u})\bm{z})$

where $\bm{z}(2^{d}-1)=1$ and $\bm{z}(a)=0$ for all other $a$. Moreover, if
$\bm{X}^{g}(\bm{u})$ is a generalized inverse of $\bm{X}(\bm{u})$, then a
generalized inverse of $\mathbb{G}(\bm{u})$ is given by:

$\mathbb{G}^{g}(\bm{u})=\mathbb{G}(\bm{X}^{g}(\bm{u})\bm{X}^{g}(\bm{u})\bm{u})$

.

Because $\bm{X}(\bm{u})$ is a triangular matrix, we can compute
$\bm{X}^{-1}(\bm{u})\bm{z}$ efficiently in $O(4^{d})$ time using back-
substitution (quadratic in the size of $\bm{z}$). Note that
$\mathbb{G}(\bm{u})$ and $\bm{X}(\bm{u})$ are invertible if and only if
$\bm{u}(2^{d}-1)>0$. The generalized inverse result holds even for non-
invertible matrices. This result is slightly more complicated, but is
important because we generally expect $\mathbb{G}$ to be singular (e.g., if it
is the Gram matrix of some workload of low-dimensional marginal query
matrices).

As we show in Theorem 9, we know the eigenvectors and eigenvalues of marginal
Gram matrices. Recall that $\bm{v}$ is an eigenvector with corresponding
eigenvalue $\lambda$ if $\mathbb{G}(\bm{w})\bm{v}=\lambda\bm{v}$ for some
real-valued $\lambda$. We use the term _eigenmatrix_ to refer to a matrix
where each column is an eigenvector that shares the same eigenvalue.

######  Theorem 9 (Eigenvectors and Eigenvalues of Marginal Gram Matrices).

Let $a\in[2^{d}]$ and let

$\mathbb{V}(a)=\bigotimes_{i=1}^{d}(a_{i}=0)\bm{T}+(a_{i}=1)(\bm{1}-n_{i}\bm{I})$

For any $b\in[2^{d}]$, $\mathbb{V}(a)$ is an eigenmatrix of $\mathbb{H}(b)$
with corresponding eigenvalue $\bm{\lambda}(a)=\bm{c}(b)$ if $a\&b=a$ and
$\bm{\lambda}(a)=0$ otherwise. Moreover, for any
$\bm{w}\in\mathbb{R}^{2^{d}}$, $\mathbb{V}(a)$ is an eigenmatrix of
$\mathbb{G}(\bm{w})$ with corresponding eigenvalue
$\bm{\kappa}(a)=\sum_{b:a\&b=a}\bm{w}(b)\bm{c}(b)$. That is,

$\displaystyle\mathbb{H}(b)\mathbb{V}(a)=\bm{\lambda}(a)\mathbb{V}(a)$
$\displaystyle\mathbb{G}(\bm{w})\mathbb{V}(a)=\bm{\kappa}(a)\mathbb{V}(a)$

Theorem 9 is remarkable because it states that the eigenmatrices (and hence
the eigenvectors) are the same for all marginal Gram matrices
$\mathbb{G}(\bm{w})$. Furthermore, the corresponding eigenvalues have a very
simple (linear) dependence on the weights $\bm{w}$. In fact, there is a
triangular matrix $\bm{Y}$ such that $\bm{\kappa}=\bm{Y}\bm{w}$.

##  8\. Optimizing conjunctive query workloads with marginal query strategies

In this section, we describe $\mbox{\sf OPT}_{\mathsf{M}}$, an optimization
operator that consumes a conjunctive query workload $\mathbb{W}$ and returns a
marginal query strategy $\mathbb{A}=\mathbb{M}(\bm{\theta})$.666We reserve the
symbol $\bm{\theta}$ for _strategies_ , and use $\bm{u},\bm{v}$ and $\bm{w}$
to refer to other marginal Gram matrices. Theorem 10 is the first key to our
approach for this problem. Intuitively, it states that for any conjunctive
query workload $\mathbb{W}$, there is a marginal Gram matrix
$\mathbb{G}(\bm{w})$ that is equivalent to $\mathbb{W}^{T}\mathbb{W}$ for the
purposes of optimization.

######  Theorem 10 (Marginal approximation of conjunctive query workload).

For any conjunctive query workload
$\mathbb{W}=w_{1}\mathbb{W}_{1}+\dots+w_{k}\mathbb{W}_{k}$, there is a
marginal Gram matrix $\mathbb{G}(\bm{w})$ such that
$tr[\mathbb{G}(\bm{u})\mathbb{W}^{T}\mathbb{W}]=tr[\mathbb{G}(\bm{u})\mathbb{G}(\bm{w})]$
for all $\bm{u}$.

$\mathbb{G}(\bm{u})$ in Theorem 10 represents $(\mathbb{A}^{T}\mathbb{A})^{+}$
in the expected error formula. We know this is a marginal Gram matrix by
Proposition 8 and Theorem 8. Theorem 10 allows us to reduce the problem of
optimizing an arbitrary conjunctive query workload to simply optimizing a
marginal query workload, which we can do efficiently. In fact, as we show in
Theorem 11, we can efficiently evaluate the matrix mechanism objective for
marginal query strategies, which is essential for efficient optimization.

######  Theorem 11 (Marginal parameterization objective function).

Let $\mathbb{W}=w_{1}\mathbb{W}_{1}+\dots+w_{k}\mathbb{W}_{k}$ be a
conjunctive query workload and let $\mathbb{G}(\bm{w})$ be the marginal
approximation of $\mathbb{W}^{T}\mathbb{W}$ (as in Theorem 10). For any
marginal query strategy $\mathbb{A}=\mathbb{M}(\bm{\theta})$, the matrix
mechanism objective function can be expressed as:

$\left\lVert\mathbb{A}\right\rVert_{\mathcal{K}}^{2}\left\lVert\mathbb{W}\mathbb{A}^{+}\right\rVert_{F}^{2}=\left\lVert\bm{\theta}\right\rVert_{\mathcal{K}}^{2}[\bm{1}^{T}\bm{X}^{+}(\bm{\theta}^{2})\bm{w}]$

where $\left\lVert\bm{\theta}\right\rVert_{\mathcal{K}}$ is the sensitivity
norm defined in Proposition 7, and $\bm{X}$ is the matrix defined in Theorem
7.

Theorem 11 shows that we can efficiently calculate the objective function in
terms of $\bm{w}$ and $\bm{\theta}$, without ever explicitly materializing
$\mathbb{G}(\bm{w})$ or $\mathbb{M}(\bm{\theta})$. This key idea will allow us
to solve the strategy selection problem efficiently. Problem 22 states the
main optimization problem that underlies $\mbox{\sf OPT}_{\mathsf{M}}$, which
immediately follows from Theorem 11.

######  Definition 22 (Marginals parameterization).

Given a conjunctive query workload
$\mathbb{W}=w_{1}\mathbb{W}_{1}+\dots+w_{k}\mathbb{W}_{k}$, let $\mbox{\sf
OPT}_{\mathsf{M}}(\mathbb{W})=\mathbb{M}(\bm{\theta}^{*})$ where

$\displaystyle\bm{\theta}^{*}=$
$\displaystyle\operatorname*{arg\,min}_{\bm{\theta}}$
$\displaystyle\left\lVert\bm{\theta}\right\rVert^{2}_{\mathcal{K}}[\bm{1}^{T}\bm{X}^{+}(\bm{\theta}^{2})\bm{w}]$
subject to
$\displaystyle\bm{X}^{+}(\bm{\theta}^{2})\bm{X}(\bm{\theta}^{2})\bm{w}=\bm{w}$

and $\mathbb{G}(\bm{w})$ is the marginal approximation of
$\mathbb{W}^{T}\mathbb{W}$ (as in Theorem 10).

Above, the constraint ensures that the strategy supports the workload. In
practice, this constraint can usually be ignored, and the resulting
unconstrained optimization problem can be solved instead. The constraint can
then be verified to hold at the end of the optimization. Intuitively, this is
because strategies that move closer to the boundary of the constraint will
have higher error, so the optimization will never approach it as long as
sufficiently small step sizes are taken. We use scipy.optimize to solve this
problem in practice.

We note that the number of parameters in the above optimization problem is
$2^{d}$ and that we can evaluate the objective in $O(4^{d})$ time (quadratic
in the number of parameters). Thus, it is feasible to solve this problem as
long as $d\leq 15$. Importantly, this means that the runtime complexity is
independent of the domain size of each attribute, so it will take the same
amount of time for $n_{i}=2$ (binary features), $n_{i}=10$, or any other
values of $n_{i}$.

In addition to being able to efficiently optimize over the space of marginal
query strategies, we can also efficiently compute the SVD bound for marginal
query workloads. Theorem 12 gives a remarkably simple formula for computing
the SVD bound for marginal query workloads.

######  Theorem 12 (SVD Bound for Marginal Query Workloads).

The SVD bound for a marginal query workload $\mathbb{W}$ with Gram matrix
$\mathbb{G}(\bm{w})$ is:

$SVDB(\mathbb{W})=\frac{1}{n}\Big{(}\sum_{a}\bm{c}(\neg
a)\sqrt{\sum_{b:a\&b=a}\bm{w}(b)\bm{c}(b)}\Big{)}^{2}$

Additionally, as a byproduct of this analysis, we give a similarly simple
formula to _find the optimal marginal query strategy in closed form_ in
Theorem 13, allowing us to bypass the need for numerical optimization in some
settings.

######  Theorem 13 (Closed form solution to Problem 22).

Let $\mathbb{W}$ be a workload with Gram matrix $\mathbb{G}(\bm{w})$ and let
$\bm{\theta}=\sqrt{\bm{Y}^{-1}\sqrt{\bm{Y}\bm{w}}}$, where $\bm{Y}$ is the
$2^{d}\times 2^{d}$ matrix with entries:

$\bm{Y}(a,b)=\begin{cases}\bm{c}(b)&a\&b=a\\\ 0&otherwise\\\ \end{cases}$

If $\bm{\theta}$ contains real-valued entries then the strategy
$\mathbb{A}=\mathbb{M}(\bm{\theta})$ attains the SVDB bound when
$\mathcal{K}=\mathcal{G}$, and is thus an optimal strategy. That is,
$\left\lVert\mathbb{A}\right\rVert^{2}_{\mathcal{G}}\left\lVert\mathbb{W}\mathbb{A}^{+}\right\rVert_{F}^{2}=SVDB(\mathbb{W})$.

While $\bm{\theta}$ may sometimes contain imaginary entries, we can always
fall back on numerical optimization to solve Problem 22. The formula in
Theorem 13 can still be used to initialize the optimization if the imaginary
entries of $\bm{\theta}$ are replaced with zeros. Li and Miklau derived
sufficient conditions for the SVD bound to be realizable [24], and marginal
query workloads satisfy those sufficient conditions. This implies that the SVD
bound should always be attainable for workloads of this form. If the
parameters in Theorem 13 contain imaginary entries, this suggests that the
optimal strategy that is _not_ a marginal query strategy. It is an interesting
open question to determine what the structure of the optimal strategy is when
Theorem 13 does not apply. In practice, even when the SVD bound is not
attained exactly by $\mbox{\sf OPT}_{\mathsf{M}}$, we get very close to it for
marginal query workloads, as we show empirically in Table 7 of the
experiments.

##  9\. The ${\mbox{\sf OPT}_{\sf HDMM}}$ strategy selection algorithm

Definition | Operator | Input workload | Output strategy | Complexity
---|---|---|---|---
§4 |  Definitions 12 and 14 | $\mbox{\sf OPT}_{0}$ | Explicit matrix ${\bm{W}}$ | Explicit matrix $\bm{A}$ | $O(n^{3})$
§6.1 | Definition 17 | $\mbox{\sf OPT}_{\otimes}$ | Kronecker Product | Kronecker product | $O(\sum_{i=1}^{d}n_{i}^{3})$
§6.2 | Definition 18 | $\mbox{\sf OPT}_{\otimes}$ | Union of Kronecker Products | Kronecker Product | $O(k\sum_{i=1}^{d}n_{i}^{3})$
§6.2 | Definition 19 | $\mbox{\sf OPT}_{+}$ | Union of Kronecker Products | Union of Kronecker Products | $O(k\sum_{i=1}^{d}n_{i}^{3})$
§8 | Definition 22 | $\mbox{\sf OPT}_{\mathsf{M}}$ | Union of Kronecker Products | Marginal Query Strategy | $O(4^{d})$

Table 2.  Summary of optimization operators: input and output types, and the
time complexity of objective/gradient calculations.

In this paper, we proposed four optimization routines: $\mbox{\sf OPT}_{0}$,
$\mbox{\sf OPT}_{\otimes}$, $\mbox{\sf OPT}_{+}$, and $\mbox{\sf
OPT}_{\mathsf{M}}$. In this section, we summarize these different approaches,
discuss the pros and cons of each one, and propose a meta-optimization
algorithm ${\mbox{\sf OPT}_{\sf HDMM}}$ that automatically chooses the best
one based on the workload. Table 2 summarizes the basic inputs and outputs of
each operator. $\mbox{\sf OPT}_{0}$ is designed to optimize an explicitly
represented workload, and returns an explicitly represented strategy. The
other optimization operators all operate in an implicit space however.

The time complexity of $\mbox{\sf OPT}_{0}$ is $O(n^{3})$ (where $n$ is the
domain size), and it generally feasible to run as long as $n\leq 10^{4}$. The
time complexity of $\mbox{\sf OPT}_{\otimes}$ and $\mbox{\sf OPT}_{+}$ is
$O(k\sum n_{i}^{3})$, where $k$ is the number of union terms in the workload,
and $n_{i}$ is the domain size of attribute $i$. It is generally feasible to
run as long as $\mbox{\sf OPT}_{0}$ is feasible on each of the individual
attributes (i.e., $n_{i}\leq 10^{4}$). In contrast to $\mbox{\sf OPT}_{0}$,
the total domain size for these operators can be arbitrarily large. The time
complexity of $\mbox{\sf OPT}_{\mathsf{M}}$ is $O(4^{d})$, which interestingly
does not depend on the domain size of individual attributes, only the number
of attributes. It is generally feasible to run as long as $d\leq 15$.

Each of the operators searches over a different space of strategies, and the
best one to use will ultimately depend on the workload. We illustrate the
behavior of each optimization operator on the simple workload of all 2-way
marginals in Example 15. This example highlights and summarizes the key
differences between $\mbox{\sf OPT}_{\otimes}$, $\mbox{\sf OPT}_{+}$, and
$\mbox{\sf OPT}_{\mathsf{M}}$. In this case, $\mbox{\sf OPT}_{\mathsf{M}}$ is
the best, which is not surprising because it is the most suitable for marginal
workloads. However, in general predicting which optimization operator will
yield the lowest error strategy requires domain expertise and may be
challenging for complex workloads. Since strategy selection is independent of
the input data and does not consume the privacy budget, we can just run each
optimization operator, keeping the output strategy that offers the smallest
expected error. Additionally, since the strategies found by each optimization
operator may depend on the initialization, we recommend running several random
restarts of each optimization operator, returning the best one.

By default, ${\mbox{\sf OPT}_{\sf HDMM}}$ invokes all three high-dimensional
optimization operators $\mbox{\sf OPT}_{\otimes}$, $\mbox{\sf OPT}_{+}$, and
$\mbox{\sf OPT}_{\mathsf{M}}$. ($\mbox{\sf OPT}_{0}$ may also be included for
lower-dimensional workloads). For $\mbox{\sf OPT}_{\otimes}$ and $\mbox{\sf
OPT}_{+}$ invoked with the p-Identity strategy we use the following convention
for setting the $p$ parameters: if an attribute’s subworkload is completely
defined in terms of $\bm{T}$ and $\bm{I}$, we set $p=1$ (this is a fairly
common case where more expressive strategies do not help), otherwise we set
$p=n_{i}/16$ for each attribute $A_{i}$ with size $n_{i}$.

######  Example 15 (Optimizing Marginal Query Workload).

Consider the workload containing queries to compute all 2-way marginals on a
domain of size $(2,5,50,100)$. This workload can be represented as a union of
$\binom{4}{2}=6$ Kronecker products. Table 3 gives the precise representation
of this workload, together with the optimized strategies found by $\mbox{\sf
OPT}_{\otimes}$, $\mbox{\sf OPT}_{+}$, and $\mbox{\sf OPT}_{\mathsf{M}}$. All
optimized strategies can be expressed in terms of the “Identity” ($\bm{I}$)
and “Total” ($\bm{T}$) building blocks. $\mbox{\sf OPT}_{\otimes}$ tries to
find the best single Kronecker product strategy, while $\mbox{\sf OPT}_{+}$
tries to find the optimal weight to assign to each of the six Kronecker
products that make up the workload. $\mbox{\sf OPT}_{\mathsf{M}}$ identifies a
different set of marginal queries from which all 2-way marginals can be
derived. Among the three optimization operators, $\mbox{\sf OPT}_{\mathsf{M}}$
is the best, followed by $\mbox{\sf OPT}_{+}$ and then $\mbox{\sf
OPT}_{\otimes}$. $\mbox{\sf OPT}_{\mathsf{M}}$ offers a $4.8\times$ reduction
in Expected TSE over the Identity baseline, and a $3.3\times$ reduction over
the Workload baseline.

| Query Matrix | Expected TSE
---|---|---
$\mathbb{W}$ |  |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ | $206,964$
|  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$
|  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$
|  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$
|  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$
|  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$
$\mathbb{I}$ |  |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ | $300,000$
$\mbox{\sf OPT}_{\otimes}(\mathbb{W})$ |  |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  $\begin{bmatrix}0.80\>{\color[rgb]{1,0,0}\bm{I}}\\\ 0.20\>{\color[rgb]{0.55,0.14,1}\bm{T}}\end{bmatrix}$ $\otimes$ |  $\begin{bmatrix}0.82\>{\color[rgb]{1,0,0}\bm{I}}\\\ 0.18\>{\color[rgb]{0.55,0.14,1}\bm{T}}\end{bmatrix}$ | $213,270$
$\mbox{\sf OPT}_{+}(\mathbb{W})$ | 0.39 |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ | $85,070$
0.18 |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$
0.14 |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$
0.13 |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$
0.11 |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$
0.05 |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$
$\mbox{\sf OPT}_{\mathsf{M}}(\mathbb{W})$ | 0.44 |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ | 62,886
0.31 |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$
0.25 |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{1,0,0}\bm{I}}$ $\otimes$ |  ${\color[rgb]{0.55,0.14,1}\bm{T}}$
Table 3.  A workload containing queries to compute _all 2-way marginals_ on a
four-dimensional domain of size $(2,5,50,100)$. The optimal strategy found by
each parameterization, and the respective error, is also shown. The Identity
and Total query matrices ${\color[rgb]{1,0,0}\bm{I}}$ and
${\color[rgb]{0.55,0.14,1}\bm{T}}$ are color coded for readability.

##  10\. Efficient MEASURE and RECONSTRUCT

Now that we have fully described how HDMM solves the strategy selection
problem, we are ready to discuss how to run the remainder of the mechanism.
Recall from Definition 10 that the matrix mechanism is defined as
$\mathcal{M}_{\bm{A},\mathcal{K}}({\bm{W}},\bm{x})={\bm{W}}\bm{A}^{+}\mathcal{K}(\bm{A},\bm{x})$.
With explicitly represented matrices, this computation can be done directly
without problem. However, HDMM replaces the explicitly represented matrices
${\bm{W}}$ and $\bm{A}$ with implicitly represented ones $\mathbb{W}$ and
$\mathbb{A}$, and it is no longer obvious how to run the mechanism.
Conceptually, these steps can be broken down as follows. In the MEASURE step,
we have to compute the noisy strategy query answers,
$\bm{y}=\mathbb{A}\bm{x}+\bm{\xi}$. In the RECONSTRUCT step, we have to
estimate the data vector and workload query ansewrs, i.e., compute
$\hat{\bm{x}}=\mathbb{A}^{+}\bm{y}$ and return $\mathbb{W}\hat{\bm{x}}$. A
necessary key subroutine to solve these problems is to compute matrix-vector
products where the matrix is a Kronecker product. Importantly, we must do this
without ever materializing $\mathbb{A}$ explicitly, as that is infeasible for
large domains.

######  Theorem 14 (Efficient matrix-vector multiplication).

Let $\mathbb{A}=\bm{A}_{1}\otimes\dots\otimes\bm{A}_{d}$ and let $\bm{x}$ be a
data vector of compatible shape. Then Algorithm 1 computes the matrix-vector
product $\mathbb{A}\bm{x}$. Furthermore, if
$\bm{A}_{i}\in\mathbb{R}^{n_{i}\times n_{i}}$ and $n=\prod n_{i}$ is the size
of $\bm{x}$ then Algorithm 1 runs in $O(n\sum n_{i})$ time.

Algorithm 1 is correct even if the factors of $\mathbb{A}$ are not square,
although the time complexity is not as clean when written down.

Algorithm 1 Kronecker Matrix-Vector Product

1:procedure kmatvec($\bm{A}_{1},\dots,\bm{A}_{d},\bm{x}$)

2: $m_{i},n_{i}=$ shape($\bm{A}_{i}$)

3: $r=\prod_{i=1}^{d}n_{i}$

4: $\bm{f}_{d+1}=\bm{x}$

5: for $i=d,\dots,1$ do

6: $\bm{Z}=$ reshape($\bm{f}_{i+1}$, $n_{i}$, $r/n_{i}$)

7: $r=r\cdot m_{i}/n_{i}$

8: $\bm{f}_{i}=$ reshape($\bm{A}_{i}\bm{Z}$, $r$, $1$)

9: end for

10: return $\bm{f}_{1}$

11:end procedure

Since all of the strategies found by our optimization routines are either
Kronecker products or unions of Kronecker products, we can directly apply
Algorithm 1 to efficiently implement the MEASURE step of HDMM. Note that
computing the matrix-vector product for a union of Kronecker products is a
trivial extension of Algorithm 1: it simply requires calling Algorithm 1 for
each Kronecker product and concatenating the results into a single vector.

We can also use Algorithm 1 to efficiently implement the RECONSTRUCT step of
HDMM. The main challenge is to compute $\mathbb{A}^{+}\bm{y}$, or a
pseudoinverse of $\mathbb{A}$ together with a matrix-vector product. This is
done slightly differently for each type of strategy:

  1. (1)

$\mathbb{A}=\mbox{\sf
OPT}_{\otimes}(\mathbb{W})=\bm{A}_{1}\otimes\dots\times\bm{A}_{d}$. From
Proposition 6 we know that
$\mathbb{A}^{+}=\bm{A}_{1}^{+}\otimes\dots\otimes\bm{A}_{d}^{+}$. That is, the
pseudoinverse of a Kronecker product is still a Kronecker product. Thus, we
can compute $\hat{\bm{x}}=\mathbb{A}^{+}\bm{y}$ efficiently using Algorithm 1.

  2. (2)

$\mathbb{A}=\mbox{\sf OPT}_{\mathsf{M}}(\mathbb{W})=\mathbb{M}(\bm{\theta})$.
From basic linear algebra, we know that
$\bm{A}^{+}=(\bm{A}^{T}\bm{A})^{+}\bm{A}^{T}$ for any matrix $\bm{A}$. Applied
to this setting, we have
$\mathbb{M}^{+}(\bm{\theta})=\mathbb{G}^{+}(\bm{\theta}^{2})\mathbb{M}^{T}(\bm{\theta})$,
since we know
$\mathbb{M}^{T}\mathbb{M}(\bm{\theta})=\mathbb{G}(\bm{\theta}^{2})$ by
Proposition 8. From Theorem 8 we know how to compute
$\mathbb{G}^{+}(\bm{\theta}^{2})$ efficiently, and we know that it equals
$\mathbb{G}(\bm{\eta})$ for some $\bm{\eta}$. We aim to compute
$\mathbb{M}^{+}(\bm{\theta})\bm{y}=\mathbb{G}^{+}(\bm{\theta}^{2})\mathbb{M}^{T}(\bm{\theta})\bm{y}$.
We can easily compute $\bm{v}=\mathbb{M}^{T}(\bm{\theta})\bm{y}$ using a
sequence of calls to Algorithm 1 by observing that
$\mathbb{M}^{T}(\bm{\theta})$ is a just a bunch of Kronecker products
_horizontally_ stacked together. In a similar fashion, we can compute
$\hat{\bm{x}}=\mathbb{G}^{+}(\bm{\theta}^{2})\bm{v}$ because
$\mathbb{G}^{+}(\bm{\theta}^{2})$ is just the sum of a bunch of Kronecker
products, which we can handle efficiently with repeated calls to Algorithm 1.

  3. (3)

$\mathbb{A}=\mbox{\sf
OPT}_{+}(\mathbb{W})=c_{1}\mathbb{A}_{1}+\dots+c_{k}\mathbb{A}_{k}$.
Unfortunately, for a strategy of this form, we do not have a way to
efficiently compute $\mathbb{A}^{+}\bm{y}$. While $\mathbb{A}$ is a union of
Kronecker products, the pseudoinverse is not necessarily, and we are not aware
of a simple formula for the pseudoinverse of $\mathbb{A}$ at all. However, we
can still produce an unbiased estimate of $\mathbb{W}\bm{x}$ by using _local
least squares_. To do this, we will compute
$\mathbb{W}_{j}\mathbb{A}^{+}_{j}\bm{y}_{j}$ for each $j=1,\dots,k$, where
$\bm{y}_{j}$ is the answers produced for sub-strategy $\mathbb{A}^{+}_{j}$.
Since $\mathbb{W}_{j}$ and $\mathbb{A}^{+}_{j}$ are both assumed to be
Kronecker products, this can be easily achieved using Algorithm 1. This
formula is an unbiased estimator for $\mathbb{W}_{j}\bm{x}$, but the main
drawback is that the workload query answers will not necessarily be consistent
between sub-workloads.

###  10.1. Improving scalability to large data vectors

Thus far, HDMM has addressed the fundamental limitation of the matrix matrix
mechanism — replacing explicit matrix representations with implicit ones, and
deriving efficient algorithms to solve the strategy optimization problem in
the implicit space. Our innovations allow HDMM to run in much higher-
dimensional settings than the matrix mechanism, but HDMM still has trouble
scaling to very high-dimensional settings, when the data vector no longer fits
in memory. Representing the data in vector form requires storing $n=\prod
n_{i}$ entries, which grows exponentially with the number of dimensions, and
quickly becomes infeasible for truly high-dimensional data. For example, a
$30$-dimensional dataset with binary attributes ($n_{i}=2$) would require
storing a data vector with $2^{30}$ entries, which is equivalent to
approximately $4\text{ GB}$ of space. Scaling beyond this point would be quite
challenging for HDMM.

It is important to note that the bottleneck of HDMM is MEASURE and
RECONSTRUCT, as these steps access and estimate the data vector. In the matrix
mechanism the main bottleneck is SELECT, as strategy optimization is the most
expensive step. HDMM can often still perform the SELECT step efficiently even
when MEASURE and RECONSTRUCT are intractable. In some special-but-common
cases, it may be possible for HDMM to bypass this bottleneck on MEASURE and
RECONSTRUCT, even scaling to settings where the data vector no longer fits in
memory, making it suitable for arbitrarily large domains.

The settings where HDMM can bypass this limitation depends crucially on the
strategy, and consequently the workload as well. If the workload is Identity
over the whole domain (or any other full rank workload), then very little can
be done because the vector of workload query answers (the output of HDMM) is
just as large as the data vector itself, and simply enumerating those answers
would require too much space. Thus, the number of queries in the workload
cannot be too large. A special-but-common case occurs when the workload
contains conjunctive queries over _small subsets of attributes_. The workload
may cover all attributes of the dataset, but it will generally consist of a
number of subworkloads, each which only cover a handful of attributes at a
time. With workloads of this form, strategies produced by HDMM ($\mbox{\sf
OPT}_{+}$ and $\mbox{\sf OPT}_{\mathsf{M}}$ in particular777strategies
produced by $\mbox{\sf OPT}_{\otimes}$ will generally be defined over the
whole domain, rather than over a small subset of attributes) will generally
contain queries that are also defined over small subsets of attributes. When
this is the case, MEASURE can be done by keeping the data in its natural
tabular format, and only vectorizing the data with respect to the relevant
attributes for each sub-workload or sub-strategy. Since these are assumed to
be defined over small subsets of attributes, these smaller data vectors can
easily be materialized explicitly and operated on accordingly. Thus, the main
remaining challenge is to RECONSTRUCT the workload query answers while
avoiding an explicit representation for $\hat{\bm{x}}$. This can be done using
a recently developed technique for efficient inference in differential privacy
called “Private-PGM” [28]. Private-PGM consumes as input a set of noisy
measurements defined over low-dimensional marginals, and produces a compact
implicit representation of $\hat{\bm{x}}$. It leverages _probabilistic
graphical models_ to compactly represent $\hat{\bm{x}}$ in terms of a product
of low-dimensional factors, and is able to scale to arbitrarily large domains
as long as the measurements allow it.

Using Private-PGM with HDMM does change the mechanism in some subtle but
important ways. The default HDMM method for RECONSTRUCT is based on standard
ordinary least squares, as it computes $\hat{\bm{x}}=\mathbb{A}^{+}\bm{y}$,
which is the solution to the minimization problem
$\hat{\bm{x}}=\operatorname*{arg\,min}_{\bm{x}}\left\lVert\mathbb{A}\bm{x}-\bm{y}\right\rVert_{2}^{2}$.
In contrast, Private-PGM is based on the related non-negative least squares
problem:
$\hat{\bm{x}}=\operatorname*{arg\,min}_{\bm{x}>0}\left\lVert\mathbb{A}\bm{x}-\bm{y}\right\rVert_{2}^{2}$.
We know that the true data vector is non-negative, so for this reason it seems
like the Private-PGM approach is more natural. However, non-negativity comes
at the cost of _bias_. An appealing property of the ordinary least squares
solution is that it produces an unbiased estimate of $\bm{x}$ under mild
conditions. Non-negative least squares does not share this same guarantee.
However, the introduction of bias often comes with reduced variance, and
overall error is usually better when enforcing non-negativity [25, 28]. Thus
Private-PGM can be used not only to improve scalability of HDMM, but also
utility. In practical settings where some bias can be tolerated for reduced
variance, we generally recommend incorporating Private-PGM post-processing
into HDMM to improve utility, even when it is not necessary for scalability
reasons.

##  11\. Experimental Evaluation

In this section we evaluate the accuracy and scalability of HDMM. We perform a
comprehensive comparison of HDMM with a variety of other mechanisms on low and
high-dimensional workloads, showing that it consistently offers lower error
than competitors and works in a broader range of settings than other
algorithms. We also evaluate the scalability of key components of HDMM,
showing that it is capable of scaling effectively to high-dimensional
settings.

In accuracy experiments, we report the Root Mean Squared Error (RMSE), which
is defined as $RMSE=\sqrt{\frac{1}{m}\text{TSE}({\bm{W}},\mathcal{K})}$ for an
algorithm $\mathcal{K}$. We compute this value analytically using the formulas
from Proposition 4 whenever possible. We separately report results for pure
differential privacy with Laplace noise and approximate differential privacy
with Gaussian noise. We use $\epsilon=1.0$ and $\delta=10^{-6}$ in all
experiments, but note that the ratio of errors between two data-independent
algorithms remains the same for all values of $\epsilon$ and $\delta$.

These experiments are meant to demonstrate that HDMM offers the best accuracy
in the _data-independent_ regime. It is possible that some data-dependent
mechanisms will outperform even the best data-independent mechanism, and this
will typically depend on the amount of data available and the privacy budget
[18, 37]. Data-independent mechanisms (like HDMM) are generally preferable
when there is an abundance of data and/or the privacy budget is not too small,
such as the U.S. Census decennial data release [1].

###  11.1. Evaluating $\mbox{\sf OPT}_{0}$ on Low Dimensional Workloads

$\bm{\epsilon}$-differential privacy (Laplace noise)
---
workload | domain | Identity | H2 | Wavelet | HB | GreedyH | LRM | OPT0 | SVDB
all-range | 64 | 6.63 | 11.28 | 10.11 | 6.63 | 6.34 | 7.02 | 5.55 | 3.22
256 | 13.11 | 16.27 | 14.87 | 8.90 | 9.72 | 15.73 | 8.07 | 4.07
1024 | 26.15 | 21.83 | 20.26 | 12.82 | 14.70 | - | 11.08 | 4.94
4096 | 52.27 | 27.90 | 26.18 | 16.19 | 22.21 | - | 14.38 | 5.82
prefix | 64 | 8.06 | 9.42 | 9.37 | 8.06 | 6.04 | 7.67 | 5.32 | 2.89
256 | 16.03 | 13.16 | 13.09 | 8.97 | 9.13 | 12.64 | 7.35 | 3.50
1024 | 32.02 | 17.29 | 17.20 | 12.87 | 14.32 | 15.43 | 9.58 | 4.11
4096 | 64.01 | 21.77 | 21.67 | 14.91 | 22.40 | - | 12.20 | 4.74
width32 | 64 | 8.00 | 12.02 | 11.09 | 8.00 | 7.32 | 9.44 | 5.88 | 2.75
256 | 8.00 | 15.50 | 13.57 | 7.41 | 8.00 | 25.81 | 6.34 | 3.26
1024 | 8.00 | 18.98 | 16.56 | 9.50 | 8.00 | 16.98 | 6.41 | 3.36
4096 | 8.00 | 22.45 | 19.58 | 10.96 | 8.00 | - | 6.46 | 3.38
permuted | 64 | 6.63 | 25.02 | 18.97 | 6.63 | 6.83 | 7.02 | 5.55 | 3.22
256 | 13.11 | 66.25 | 49.09 | 18.48 | 13.02 | 15.73 | 8.06 | 4.07
1024 | 26.15 | 157.50 | 117.06 | 37.07 | 23.94 | - | 11.08 | 4.94
4096 | 52.27 | 374.29 | 277.42 | 107.83 | 45.77 | - | 14.37 | 5.82

Table 4. Error of strategies for 1D workloads with $\epsilon=1.0$.

$(\bm{\epsilon},\bm{\delta})$-differential privacy (Gaussian noise)
---
workload | domain | Identity | H2 | Wavelet | HB | GreedyH | COA | OPT0 | SVDB
All Range | 64 | 19.82 | 12.74 | 11.42 | 19.82 | 14.64 | 9.73 | 9.73 | 9.62
256 | 39.18 | 16.20 | 14.81 | 18.80 | 23.34 | 12.26 | 12.26 | 12.15
1024 | 78.13 | 19.66 | 18.24 | 27.07 | 36.20 | 14.89 | 14.85 | 14.75
4096 | 156.14 | 23.12 | 21.69 | 27.92 | 56.21 | 17.92 | 17.46 | 17.38
prefix | 64 | 24.08 | 10.64 | 10.58 | 24.08 | 14.04 | 8.87 | 8.87 | 8.62
256 | 47.89 | 13.11 | 13.03 | 18.95 | 22.11 | 10.70 | 10.66 | 10.44
1024 | 95.64 | 15.57 | 15.49 | 27.18 | 35.59 | 16.29 | 12.49 | 12.29
4096 | 191.21 | 18.03 | 17.95 | 25.72 | 56.70 | 26.50 | 14.32 | 14.15
width32 | 64 | 23.90 | 13.57 | 12.52 | 23.90 | 17.30 | 8.79 | 8.74 | 8.23
256 | 23.90 | 15.44 | 13.52 | 15.65 | 23.90 | 12.24 | 9.93 | 9.73
1024 | 23.90 | 17.10 | 14.92 | 20.08 | 23.90 | 16.00 | 10.08 | 10.02
4096 | 23.90 | 18.60 | 16.22 | 18.90 | 23.90 | 18.38 | 10.11 | 10.09
permuted | 64 | 19.82 | 28.26 | 21.42 | 19.82 | 16.13 | 9.73 | 9.73 | 9.62
256 | 39.18 | 65.97 | 48.88 | 39.04 | 35.22 | 12.26 | 12.26 | 12.15
1024 | 78.13 | 141.86 | 105.44 | 78.30 | 60.60 | 14.89 | 14.85 | 14.75
4096 | 156.14 | 310.11 | 229.85 | 185.98 | 118.03 | 17.92 | 17.45 | 17.38

Table 5. Error of strategies for 1D workloads with $\epsilon=1.0$ and
$\delta=10^{-6}$.

We begin by studying the effectiveness of $\mbox{\sf OPT}_{0}$ in the one-
dimensional setting. Specifically, we evaluate the quality of the strategies
found by our optimization oracle compared with other _data-independent_
mechanisms designed for this setting. It is important to understand the
accuracy in the one-dimensional setting well, because $\mbox{\sf OPT}_{0}$ is
used as a sub-routine for the higher-dimensional optimization operators
$\mbox{\sf OPT}_{\otimes}$ and $\mbox{\sf OPT}_{+}$.

#### Workloads

We consider four different workloads: All Range, Prefix, Width 32 Range, and
Permuted Range, each defined over domain sizes ranging from $64$ to $4096$.
All Range contains every possible range query over the specified domain;
Prefix contains range queries defining an empirical CDF; Width 32 Range
contains all range queries of width 32. While the first three workloads are
subsets of range queries, the last workload, Permuted Range, is the result of
right-multiplying the workload of all range queries by a random permutation
matrix. Many proposed strategies have targeted workloads of range queries and
tend to work fairly well on subsets of range queries. Permuted Range poses a
challenge because the structure of the workload is hidden by the permutation,
requiring a truly adaptive method to find a good strategy.

Note the large size of some of these workloads: All Range and Permuted Range
have $\frac{n(n+1)}{2}$ queries. For large $n$ it is infeasible to write down
${\bm{W}}$ in matrix form, but we can still compute the expected error since
it only depends on the workload through its Gram matrix,
${\bm{W}}^{T}{\bm{W}}$, which is $n\times n$ and has special structure,
allowing it to be computed directly without materializing ${\bm{W}}$.

#### Mechanisms

We consider 8 competing mechanisms: Identity, Laplace, Gaussian, LRM [44], COA
[43], H2 [19], HB [33], Privelet [38], and GreedyH [21]. The first five
mechanisms are general purpose mechanisms, designed to support virtually any
workload. The last four mechanisms were specifically designed to offer low
error on range query workloads. We also report SVDB to understand the gap
between the error of the computed strategies and the best lower bound on error
we have (via the SVD bound).

#### Results and Findings

Table 4 and Table 5 report the error of various mechanisms in each setting,
for both Laplace and Gaussian noise respectively. We remind the reader that
these values do not depend on the true data $\bm{x}$, and thus they hold for
all $\bm{x}$. We report numbers for fixed $\epsilon=1.0$ and $\delta=10^{-6}$,
but we note that these privacy parameters only impact the error by a constant
factor, and hence the relationship between the errors of every pair of
mechanisms remains the same for all $(\epsilon,\delta)$. We have four main
findings from these results, enumerated below:

  1. (1)

$\mbox{\sf OPT}_{0}$ offers lower error than all competitors in all settings,
and the magnitude of the improvement offered by HDMM (over the next best
competitor) is as large as $3.18$ for Laplace noise (on Permuted Range) and
$1.61$ for Gaussian noise (on Width 32 Range). Interestingly, $\mbox{\sf
OPT}_{0}$ offers lower error than H2, HB, Privelet, and GreedyH on range query
workloads, even though these four mechanisms were designed specifically for
range queries. In additional, the second best method after $\mbox{\sf
OPT}_{0}$ differs in each setting, which shows that some competing algorithms
have specialized capabilities that allow them to perform well in some
settings, while HDMM performs well in a variety of settings as it does not
make strict assumptions about the workload.

  2. (2)

$\mbox{\sf OPT}_{0}$ gets within a factor of $2.57$ of the SVD bound for
Laplace noise and $1.01$ of the SVD bound for Gaussian noise on every tested
workload. The gap between $\mbox{\sf OPT}_{0}$ and SVDB is quite small for
Gaussian noise, suggesting that $\mbox{\sf OPT}_{0}$ is finding the best
possible strategy. Note that COA also finds a optimal strategy in many of the
settings, but it fails on the Prefix and Width 32 Range workloads when $n\geq
1024$. Thus, even though it is solving the same problem underlying $\mbox{\sf
OPT}_{0}$ in theory, the implementation is not as robust as ours. The gap
between $\mbox{\sf OPT}_{0}$ and SVDB is larger for Laplace noise, however,
and it is unclear if this gap is primarily due to looseness of the SVD bound
or suboptimality of the strategy. Nevertheless, even with Laplace noise the
ratio between $\mbox{\sf OPT}_{0}$ and SVDB is at most $2.57$.

  3. (3)

The error of $\mbox{\sf OPT}_{0}$ (and COA for $(\epsilon,\delta)$ privacy) is
the same on the All Range and Permuted Range workloads. Permuting the workload
doesn’t impact achievable error or our optimization algorithm in any
meaningful way. However, many of the methods we compared against perform well
on All Range but poorly on Permuted Range because they were specifically
designed for range queries. This shows that they exploit specific structure of
the input workload and have limited adaptivity.

  4. (4)

On these workloads, Laplace noise offers better error than Gaussian noise (for
appropriately conservative settings of $\delta$). This is because with
Gaussian noise there is an additional $\approx\sqrt{\log(1/\delta)}$ term in
the standard deviation of the noise, and this outweighs the benefit using the
$L_{2}$ sensitivity norm instead of the $L_{1}$ sensitivity norm, despite the
fact that we may be finding strategies that are further from optimal in the
$L_{1}$ case.

#### Scalability

Figure 2.  Time required to run $\mbox{\sf OPT}_{0}$ for $100$ iterations on
the AllRange workload for increasing domain sizes.

We now demonstrate the scalability of $\mbox{\sf OPT}_{0}$. Note that
optimization time dominates in the low-dimensional setting, and the time for
MEASURE and RECONSTRUCT is small in comparison to that. The per-iteration time
complexity only depends on the domain size, and not the contents of the
workload. While the number of iterations required for convergence may differ
slightly based on the queries in the workload, for simplicity we measure the
time required to run the optimization for 100 iterations on the All Range
workload.

Figure 2 shows the amount of time required to run $\mbox{\sf OPT}_{0}$ for
various domain sizes. It shows that $\mbox{\sf OPT}_{0}$ scales up to
$n=8192$, and runs for $n=1024$ in under 10 seconds for Laplace noise and 1
minute for Gaussian noise. This difference occurs because the per-iteration
time complexity is $O(pn^{2})$ under Laplace noise but $O(n^{3})$ under
Gaussian noise. For $n=8192$ it takes considerably longer, but is still
feasible to run. We remark that trading a few hours of computation time for a
meaningful reduction in error is typically a welcome trade-off in practice,
especially since workloads can be optimized once and the resulting strategies
reused many times. Additionally, we have a prototype implementation that uses
GPUs and PyTorch, and we found that it is possible (although very time
consuming) to scale up to $n=16384$. Beyond this point, it quickly becomes
infeasible to even represent the workload (or its Gram matrix) in matrix form,
let alone optimize it.

###  11.2. Evaluating $\mbox{\sf OPT}_{\otimes}$, $\mbox{\sf OPT}_{+}$, and
$\mbox{\sf OPT}_{\mathsf{M}}$ on Multi-Dimensional Workloads

$\bm{\epsilon}$-differential privacy (Laplace noise)
---
dataset | workload | Identity | Laplace | DataCube | OPT$\bm{\otimes}$ | OPT$\bm{+}$ | OPTM | HDMM | SVDB
Census (5D) | SF1 | 23.20 | 70.71 | - | 7.30 | 30.55 | 9.56 | 7.30 | -
SF1+ | 32.50 | 141.42 | - | 10.23 | 42.31 | 12.88 | 10.23 | -
CPS (5D) | All Marginals | 5.38 | 45.25 | 18.49 | 4.85 | 4.85 | 4.84 | 4.84 | 2.63
All Prefix-Marginals | 98.06 | 56568.54 | - | 40.59 | 40.59 | 69.38 | 40.59 | 9.32
Adult (14D) |  $\leq$ 3D Marginals | 5352117.26 | 664.68 | 494.06 | 872.58 | 306.33 | 225.35 | 225.35 | 15.08
2D Prefix-Marginals | 475602516.60 | 138602.83 | - | 1119.16 | 484.07 | 553.56 | 484.07 | -
Loans (12D) | Small Marginals | 3330650.46 | 265.87 | 113.98 | 654.35 | 204.17 | 100.92 | 100.92 | 11.61
Small Prefix-Marginals | 15340082.96 | 11013.90 | - | 1707.67 | 485.67 | 288.29 | 288.29 | -

Table 6. RMSE of HDMM strategies and baseline strategies on multi-dimensional
workloads (ranging from 5D to 14D) for $\epsilon=1.0$ with Laplace noise.

$(\bm{\epsilon},\bm{\delta})$-differential privacy (Gaussian noise)
---
dataset | workload | Identity | Laplace | DataCube | OPT$\bm{\otimes}$ | OPT$\bm{+}$ | OPTM | HDMM | SVDB
Census (5D) | SF1 | 69.33 | 29.87 | - | 9.80 | 75.66 | 14.31 | 9.80 | -
SF1+ | 97.08 | 42.25 | - | 10.90 | 84.16 | 15.88 | 10.90 | -
CPS (5D) | All Marginals | 16.08 | 23.90 | 19.53 | 7.85 | 7.85 | 7.86 | 7.85 | 7.85
All Prefix-Marginals | 292.93 | 844.94 | - | 29.48 | 29.49 | 104.09 | 29.48 | 27.85
Adult (14D) |  $\leq$ 3D Marginals | 15988375.02 | 91.59 | 77.36 | 82.42 | 899.04 | 46.44 | 46.44 | 45.06
2D Prefix-Marginals | 1420766966.19 | 1322.58 | - | 126.17 | 639.43 | 296.12 | 126.17 | -
Loans (12D) | Small Marginals | 9949649.11 | 57.92 | 37.37 | 81.51 | 631.40 | 34.91 | 34.91 | 34.67
Small Prefix-Marginals | 45825415.90 | 372.83 | - | 132.43 | 994.04 | 99.72 | 99.72 | -

Table 7. RMSE of HDMM strategies and baseline strategies on multi-dimensional
workloads (ranging from 5D to 14D) for $\epsilon=1.0$ and $\delta=10^{-6}$
with Gaussian noise.

We now shift our attention to the multi-dimensional setting.

#### Workloads

We consider four multi-dimensional schemas and two workloads for each schema.
The first schema, Census of Population and Housing (Census), has been used as
a running example throughout the paper. The second schema, Current Population