THE INVERSE PROBLEMS OF SOME MATHEMATICAL

PROGRAMMING PROBLEMS

HUANG Siming

Institute of Policy and Management,

Academia Sinica, Beijing 100080, P.R. China.

Abstract

The non-convex quadratic programming problem and non-monotone linear
complementarity problem are NP-complete problems. In this paper we first show
that the inverse problem of determining a KKT point of the non-convex
quadratic programming is polynomial. We then show that the inverse problems of
non-monotone linear complementarity problem are polynomial solvable in some
cases, and in another case is NP-hard. Therefore we solve an open question
raised by Heuberger on inverse NP-hard problems and prove the CoNP=NP.

Key words: Inverse problem, quadratic programming, Linear complementarity
problem, Polynomial algorithm.

Abbreviated title: Inverse problems of mathematical programming problems.

1\. Introduction

Since Burton and Toint [3] first introduced inverse shortest path problem,
many researchers have contributed to the growing literature on inverse
optimization problems. For example, Burton and Toint [3],[4], and Burton,
Pulleyblank and Toint [5] have discussed inverse shortest path problem; Huang
and Liu [9],[10] have considered inverse linear programming problem and
applied it to inverse matching problem and inverse minimum cost flow problem
respectively; Zhang, Liu and Ma [19], Sokkalingam, Ahuja and Orlin [18], and
Ahuja and Orlin [2] studied inverse minimum spanning tree problem, etc. For a
more complete survey on inverse combinatorial optimization problems we refer
the reader to Heuberger [7]. Some researchers also did some work on the
inverse mathematical programming problems. For example, Huang and Liu [9], and
Zhang and Liu [20] studied inverse linear programming problem, Diao and Ding
[6] studied inverse convex programming problem. Most of the inverse problems
studied so far are polynomial solvable problems, and most of their inverse
problems are polynomial problems too. A natural open question arises: is there
a NP-complete problem such that it’s inverse problem is polynomial solvable?
This question was proposed in Huang and Liu [9] as an open question for future
research. In a survey paper written by Heuberger [7] on inverse combinatorial
optimization problems, he proposed three open questions for future research.
The first of them is: is there a NP-hard problem such that its inverse problem
is polynomial solvable? If so one proves CoNP=NP, a long standing open problem
in the theory of computational complexity. Recently, Huang [11] have shown
that the inverse Knapsack problem and the inverse problem of integer
programming with fixed number of constraints are pseudo polynomial. In this
paper, we will further show that the inverse problem of determining the KKT
points of a non-convex quadratic programming problem is polynomial solvable;
we then show that the inverse problem of determining the optimal solutions of
non-monotone linear complementarity problem in some cases are polynomial
solvable, in another case is NP-hard.

The theory of inverse optimization can also be viewed as the generalization or
complement to the theory of sensitivity analysis and the theory of stability
of optimal solutions of optimization. For example, the sensitivity analysis of
linear programming (LP) study the range of changes of each parameter such that
the optimal solution will not change if the changes of each parameter is
within the range. If the change of one of the parameter is outside the range,
then the optimal solution will change. The inverse linear programming [9]
study the problem of given any feasible solution $x^{0}$, how to change the
parameters such that $x^{0}$ becomes an optimal solution of the LP with the
change of the parameters as small as possible under the sense of $l_{p}$ norm.
Therefore it can be viewed as the generalization of the sensitivity analysis
of linear programming. It has been shown in [9] that the inverse problem of
linear programming can be solved efficiently.

The sensitivity analysis of NP-complete problems are generally viewed as
difficult since there is no polynomial algorithm to solve them. The results in
this paper show that for some NP-complete problems, for example the KKT
problem of a non-convex quadratic programming and non-monotone LCP problem,
the inverse problems of them can be efficiently solved. Therefore it is
relatively easy to do the sensitivity analysis for these problems.

The paper is organized as follows: in section 2 we introduce the quadratic
programming and show that the inverse problem of finding a KKT point of the
non-convex quadratic programming is polynomial solvable. In section 3, we show
that the inverse problem of finding an optimal solution of a class of non-
monotone linear complementarity problem in some cases are polynomial solvable,
in other case is NP-hard. In section 4 we will show that CoNP=NP using the
results in section 2. We will give some concluding remarks in section 5.

2\. The Inverse Problem of a KKT Point of Non-convex QP

The concept of inverse problem of an optimization problem can be found in [8].
For completion, we state it here again.Given an optimization problem:

$\hbox{min}\\{f(c,x)|x\in D\\},$ $None$

where $c\in R^{n}$ is a parameter vector, $D$ is the feasible region of $x$,
$f(c,x)$ is the objective function. Given a feasible solution $x^{0}$ of
$(1)$, is there $\bar{c}\in R^{n}$ such that $x^{0}$ is the optimal solution
of $(1)$ with $\bar{c}$ as parameter vector? Formally, let

$F(x^{0})=\\{\bar{c}\in R^{n}|\hbox{min}\\{f(\bar{c},x)|x\in
D\\}=f(\bar{c},x^{0})\\},$

if $F(x^{0})\not=\emptyset$, define

$\hbox{min}\\{\|c-\bar{c}\||\bar{c}\in F(x^{0})\\},$ $None$

where $\|.\|$ denotes the norm of the vector, the popular choices for the
norms are $l_{1}$, $l_{2}$ and $l_{\infty}$. we call $(2)$ the inverse problem
of $(1)$.

Now we consider the following general non-convex quadratic programming
problem:

7 (QP)&minimize$Q(x)={1\over 2}x^{T}Qx+c^{T}x$ subject to $Ax=b,\ x\geq 0,$

where $Q\in R^{n\times n}$ is an indefinite matrix, $A\in R^{m\times n}$ and
$c\in R^{n}$. Its dual is:

7 (QD)maximize$d(x,y)=b^{T}y-{1\over 2}x^{T}Qx-c^{T}x$ subject
to$A^{T}y+s-Qx=c,\ x,s\geq 0,$

Let $F(P)=\\{x|Ax=b,x\geq 0\\}$ and $F(D)=\\{(x,y,s)|A^{T}y+s-Qx=c,x,s\geq
0\\}$ be the feasible regions of the $(QP)$ and $(QD)$ respectively, then the
KKT conditions of the $(QP)-(QD)$ are:

$\eqalignno{Ax&=b,\cr A^{T}y+s-Qx&=c,\cr x^{T}s&=0.\cr}$

Where $(x,y,s)\in(R^{n}_{+},R^{m},R^{n}_{+})$ is a feasible solution of
$(QP)-(QD)$. As proved in [15], finding a KKT point of the nonconvex quadratic
programming is NP-complete.

Therefore the inverse problem of determining the KKT point of the $(QP)$ can
be stated as follows: given any $x^{0}\in F(P)$, find a matrix $\bar{Q}$ and a
vector $\bar{c}$ that are closest to $Q$ and $c$ such that $x^{0}$ is a KKT
point of $(QP)$ with the parameters $\bar{Q}$ and $\bar{c}$. Let $\bar{Q}$
denotes the matrix after changing the parameters in $Q$, $\bar{c}$ denotes the
vector after changing the parameters in $c$, then the inverse problem of a KKT
point of $(QP)$ with $l_{2}$ norm is:

7

(IQP)Min$\|Q-\bar{Q}\|^{2}_{f}+\|c-\bar{c}\|^{2}$ s.t.$Ax^{0}=b$
$A^{T}y+s-{\bar{Q}}x^{0}=\bar{c}$ ${x^{0}}^{T}s=0$ $s\geq 0$.

Where $\|A\|_{f}=\sqrt{\sum_{i,j}A^{2}_{ij}}$ denote the Frobenius norm of the
matrix $A$. Let $X=Q-\bar{Q}$, $z=c-\bar{c}$, $c^{\prime}=Qx^{0}+c$. We define
the $n\times n^{2}$ block diagonal matrix $X^{0}$ as follows:

$X^{0}=\pmatrix{{x^{0}}^{T}&&\cr&{x^{0}}^{T}&\cr&\dots&\cr&&{x^{0}}^{T}\cr}.$
$None$

That is, the first n entries of the first row of $X^{0}$ are the row vector of
$x^{0}$ and rest entries of the first row are $0$; the $n$ to $2n$ entries of
the second row of $X^{0}$ are the row vector of $x^{0}$ and rest entries are
$0$; and so on. Let $A=(a_{1},...,a_{n})$ be any $n\times n$ matrix, where
$a_{j}$ $(j=1,...,n)$ are columns of $A$, we define the vector of $A$ as
follows:

$vec(A)=(a^{T}_{1},...,a^{T}_{n})^{T}.$

Then (IQP) becomes:

7

(IQP1)Min$\|X\|^{2}_{f}+\|z\|^{2}$

s.t.$A^{T}y+s+{X^{0}}vec(X^{T})+z=c^{\prime}$

${x^{0}}^{T}s=0$

$s\geq 0$.

It is easy to see that the above problem is a convex quadratic optimization
problem with $n^{2}+2n+m$ variables. In order to apply the interior point
algorithms for solving $(IQP1)$, we do some transformations for $(IQP1)$
first.

We define $w=(vec(X^{T})^{T},z^{T},y^{T})^{T}$ to be a $\\{n^{2}+n+m\\}$
dimensional vector, $A^{\prime}=(X^{0},I,A^{T})$ to be a
$n\times\\{n^{2}+n+m\\}$ matrix, $P=diag(1,...1,0,...0)$ to be a
$\\{n^{2}+n+m\\}\times\\{n^{2}+n+m\\}$ diagonal matrix with first $n^{2}+n$
diagonal entries being $1$, rest $m$ diagonal entries being $0$. Note that
$\|X\|_{f}^{2}=vet(X)^{T}vec(X)=vec(X^{T})^{T}vec(X^{T})$. Then $(IQP1)$ can
be further transformed into:

7

(IQP2)Min$w^{T}Pw$

s.t.$s+A^{\prime}w=c^{\prime}$,

${x^{0}}^{T}s=0$,

$s\geq 0$.

Define non-negative vectors $w^{+}$ and $w^{-}$ as follows:

$w^{+}_{i}=\cases{w_{i},&if $w_{i}\geq 0$;\cr 0,&otherwise;\cr}\,\qquad
w^{-}_{i}=\cases{-w_{i},&if $w_{i}\leq 0$;\cr 0,&otherwise.\cr}$ $None$

Let $\bar{w}=({w^{+}}^{T},{w^{-}}^{T})^{T}$,
$\bar{A}=(A^{\prime},-A^{\prime})$, and

$\bar{P}=\pmatrix{2P&-2P\cr-2P&2P\cr}.$

Then $(IQP2)$ becomes:

7

(IQP3)Min${1\over 2}{\bar{w}}^{T}\bar{P}\bar{w}$

s.t$s+{\bar{A}}{\bar{w}}=c^{\prime}$,

${x^{0}}^{T}s=0,$

$s\geq 0$, ${\bar{w}}\geq 0$.

It is easy to show that $\bar{P}$ is positive semidefinite since $P$ is
positive semidefinite. Therefore $(IQP3)$ is a standard convex quadratic
programming problem. Hence one can use any polynomial interior point algorithm
for convex quadratic programming, for example [15], to solve this convex
quadratic programming. We conclude this section by following theorem.

Theorem 2.1 The inverse problem of determining a KKT point of the non-convex
quadratic programming $(IQP)$ can be solved in polynomial time.

Note that if the $l_{1}$ norm is used in (IQP), then it can be transformed
into a linear programming problem. We leave it to the reader to try it.

3\. The Inverse Problem of Linear Complementarity Problem

In this section, we discuss the inverse problem of non-monotone linear
complementarity problem. Given $M\in R^{n\times n}$ and $q\in R^{n}$, the
linear complementarity problem (LCP) is to find a pair $(x,s)\in R^{n}\times
R^{n}$ such that

$s=Mx+q,\,\ (x,s)\geq 0,\,\ \hbox{and}\,\ x^{T}s=0.$

We denote the feasible set of (LCP) as:

$F=\\{(x,s)|s=Mx+q,\,\ x,s\geq 0\\}.$

The LCP problem can also be formulated as an optimization problem:

7(LCP)Min$x^{T}s$ s.t.$s=Mx+q,\,\ x,s\geq 0$.

If $M$ (may not be symmetric) is positive semi-definite, then the LCP is
called a monotone LCP. Otherwise it is called a non-monotone LCP. Monotone
LCPs are “easy” problems since there are polynomial algorithms for solving the
monotone LCPs. See, for example, [13] and [14]. If $M$ is not positive semi-
definite, then the LCP becomes a “hard” problem, i.e, NP-complete problem,
since it can be transformed into the following non-convex quadratic
programming problem:

7Min$x^{T}Mx+q^{T}x$ s.t.$s=Mx+q,\,\ x,s\geq 0$.

In fact, any feasibility problem of a mixed integer programming problem can be
transformed into a LCP problem, see [8]. Therefore non-monotone LCP problems
are NP-hard.

To study the inverse problem of non-monotone LCP, we consider three different
cases, depends on what information were available to us:

(a): Given $x^{0}\geq 0$, find parameters $M^{\prime}$ and $q^{\prime}$ and
$s^{0}\geq 0$ such that $(x^{0},s^{0})$ is an optimal solution of the LCP
problem with parameters of $M^{\prime}$, $q^{\prime}$ and that $M^{\prime}$,
$q^{\prime}$ are closest to $M$, $q$.

(b): Given a pair of non-feasible complementary point $(x^{0},s^{0})\geq 0$
(i.e. ${x^{0}}^{T}s^{0}=0$, $s^{0}\not=Mx^{0}+q$), find parameters
$M^{\prime}$ and $q^{\prime}$ such that $(x^{0},s^{0})$ is an optimal solution
of the LCP problem with parameters of $M^{\prime}$, $q^{\prime}$ and that
$M^{\prime}$, $q^{\prime}$ are closest to $M$, $q$.

(c): Given $s^{0}\geq 0$, find parameters $M^{\prime}$ and $q^{\prime}$ and
$x^{0}\geq 0$ such that $(x^{0},s^{0})$ is an optimal solution of the LCP
problem with parameters of $M^{\prime}$, $q^{\prime}$ and that $M^{\prime}$,
$q^{\prime}$ are closest to $M$, $q$.

We now discuss the case (a) first. The inverse problem of LCP in this case
(with $l_{2}$ norm) is:

7

(ILCP(a))Min$\|M-M^{\prime}\|^{2}_{f}+\|q-q^{\prime}\|^{2}$

s.t.$s=M^{\prime}x^{0}+q^{\prime},$

${x^{0}}^{T}s=0,\,\ s\geq 0.$

Similar to the discussions in last section, we can transform the $(ILCP(a))$
into a convex programming problem. Let $X=M-M^{\prime}$, $z=q-q^{\prime}$,
matrix $X^{0}$ be defined as in $(3)$, $c=Mx^{o}+q$, then the $(ILCP(a))$
becomes:

7

(ILCP(a)1)Min$\|X\|^{2}_{f}+\|z\|^{2}$

s.t.$s+X^{0}vec(X^{T})+z=c$,

${x^{0}}^{T}s=0$,

$s\geq 0.$

Let $w^{T}=(vec(X^{T})^{T},z^{T})$ be a $\\{n^{2}+n\\}$ dimensional vector,
$A^{\prime}=(X^{0},I)$ be a $n\times\\{n^{2}+n\\}$ matrix, then (ILCP(a)1)
becomes:

7

(ILCP(a)2)Min$w^{T}w$

s.t.$s+A^{\prime}w=c$,

${x^{0}}^{T}s=0$,

$s\geq 0$.

Define non-negative vectors $w^{+}$ and $w^{-}$ as in $(4)$, and
$\bar{w}=({w^{+}}^{T},{w^{-}}^{T})^{T}$, $\bar{A}=(A^{\prime},-A^{\prime})$,
then $(ILCP(a)2)$ becomes:

7

(ILCP(a)3)Min${\bar{w}}^{T}\bar{w}$

s.t$s+{\bar{A}}{\bar{w}}=c$,

${x^{0}}^{T}s=0,$

$s\geq 0$, ${\bar{w}}\geq 0$.

It is easy to see that $(ILCP(a)3)$ is a convex quadratic programming problem,
therefore it can be solved by any polynomial interior point algorithms for
quadratic programming. See, for example, [14]. Hence we have the following
conclusion.

Theorem 3.1. Given any $x^{0}\geq 0$, the inverse problems of $(LCP)$ in case
(a), i.e., $ILCP(a)$ can be solved in polynomial time.

For the case (b), the inverse problem (with $l_{2}$ norm) becomes:

7

(ILCP(b))Min$\|M-M^{\prime}\|^{2}_{f}+\|q-q^{\prime}\|^{2}$

s.t.$s^{0}=M^{\prime}x^{0}+q^{\prime},$

${x^{0}}^{T}s^{0}=0.$

Similarly, let $X=M-M^{\prime}$, $z=q-q^{\prime}$, matrix $X^{0}$ be defined
as in $()$, vector $c=Mx^{0}+q-s^{0}$, $w^{T}=(vec(X^{T})^{T},z^{T})$, non-
negative vectors $w^{+}$ and $w^{-}$ be defined as in $(4)$, and
$\bar{w}=({w^{+}}^{T},{w^{-}}^{T})^{T}$, $A^{\prime}=(X^{0},I)$ and
$\bar{A}=(A^{\prime},-A^{\prime})$, then the $(ILCP(b))$ becomes:

7

(ILCP(b)1)Min${\bar{w}}^{T}\bar{w}$

s.t.$\bar{A}\bar{w}=c$,

$\bar{w}\geq 0$.

It is also easy to see that $(ILCP(b)1)$ is a convex quadratic programming
problem, therefore it can also be solved by any polynomial interior point
algorithms for quadratic programming. Hence we have the following theorem.

Theorem 3.2. Given any complementary pair $(x^{0},s^{0})\geq 0$, the inverse
problem of $(LCP)$ in case (b), i.e., $ILCP(b)$, can be solved in polynomial
time.

For the case (c), the inverse problem with $l_{1}$ norm becomes:

7

(ILCP(c))Min$\|M-M^{\prime}\|+\|q-q^{\prime}\|$

s.t.$s^{0}=M^{\prime}x+q^{\prime},$

$x^{T}s^{0}=0,$

$x\geq 0$.

Let $Y=M-M^{\prime}$, $z=q-q^{\prime}$, $c=q-s^{0}$, then (ILCP(c)) becomes:

7

(ILCP(c)1)Min$\sum^{n}_{i,j=1}|Y_{ij}|+\sum^{n}_{i=1}|z_{i}|$

s.t.$Yx-Mx+z=c$,

$x^{T}s^{0}=0,$

$x\geq 0$.

Let $Y_{j}$ denotes the jth column of the matrix $Y$ $(j=1,...,n)$, then
$Y^{T}_{j}$ is the jth row of the matrix $Y$ $(j=1,...,n)$. Therefore the
constraints

$Yx-Mx+z=c$

in (ILCP(c)1) can be written in $n$ quadratic constraints:

$x^{T}(Y^{T})_{j}-x^{T}(M^{T})_{j}+z_{j}=c_{j},\,\ j=1,...,n.$

Hence (ILCP(c)1) is NP-hard since an optimization problem with linear
objective function and one non-convex quadratic constraint is a NP-complete
problem. See for example, [11]. Now we can make following conclusion:

Theorem 3.3. Given any $s^{0}\geq 0$, the inverse problem of (LCP) in case (c)
is NP-hard.

4\. The Proof of CoNP=NP

As pointed out by Heuberger in [7] as one of the open problems: if one can
find a NP-hard optimization problem such that its inverse problem can be
solved in polynomial time, then one proves CoNP=NP, a long standing open
problem in computational complexity. But there he only gave a brief
discussion, no formal proof. In this section we give a formal proof of CoNP=NP
using the results obtained in above sections. Let us first review some of the
useful concepts and results related to CoNP. We will follow the notations in
[17].

Definition 4.1:(See P. 384 of [17]) Problem $\bar{A}$ is the complement of
Problem $A$ if the set of strings with symbols in $\sum$ that are encodings of
yes instances of $\bar{A}$ are exactly those that are not encodings of yes
instances of $A$.

Definition 4.2:(P. 385 of [17]) The class CoNP is the class of all problems
that are complements of problems in NP.

Theorem 4.3:(P. 385 of [17]) If the complement of an NP-complete problem is in
NP, then NP=CoNP.

Now let us consider the following decision problem:

P1 Given an instance of $(Q,c,A,b)$, is there a feasible solution $x^{0}$ such
that it is a KKT point of $(QP)$?

It is well known that $P1$ is NP-complete [16]. By Definition 4.1 the
complement of $P1$ is:

P2: All the feasible solutions of $(QP)$ are not the KKT points of $(QP)$, or
equivalently $(QP)$ has no KKT point.

Regarding the $P2$ we have the following theorem.

Theorem 4.4: $P2$ is in NP. Proof: Given any feasible solution $x^{0}$ of
$(QP)$, we can check whether it is a KKT point of $(QP)$ in polynomial time by
solving the inverse problem of $(QP)$, i.e., $(IQP)$. $x^{0}$ is not a KKT
point of $(QP)$ if and only if $\bar{Q}\not=Q$ or $\bar{c}\not=c$. Since
$(QP)$ have many feasible solutions (more than exponential number), so $P2$ is
in NP.

Using theorem 4.3 and theorem 4.4 we have the following main result of the
section.

Theorem 4.5 CoNP=NP.

5\. Conclusions

In this paper we have shown that the inverse problem of determining the KKT
point of the non-convex quadratic programming is polynomial solvable. We also
showed that the inverse problems of determining the optimal solutions of the
non-monotone linear complementarity problem in cases (a) and (a) are
polynomial, in case (c) is NP-hard. We showed that the inverse problem of
determining a KKT point of the non-convex $(QP)$ with $l_{2}$ norm is
equivalent to solve a convex quadratic programming, therefore can be solved by
polynomial interior point algorithms. We then showed that the inverse problems
of determining an optimal solution of the non-monotone linear complementarity
problem with $l_{2}$ norm in cases (a) and (b) are also equivalent to solve a
convex quadratic programming problem, hence can also be solved by polynomial
interior point algorithms. In case (c), we showed that the inverse problem of
non-monotone LCP with $l_{1}$ norm is equivalent to solve an optimization
problem with linear objective function and $n$ non-convex quadratic
constraints, therefore is NP-hard. The results in this paper solved an open
question raised by Heuberger in [7] on whether there exists an NP-hard problem
such that its inverse problem is polynomial solvable. If so one proves
CoNP=NP, a long standing open problem in the theory of computational
complexity.

It is also interesting to observe that finding the KKT points of a general
quadratic programming problem and solving the non-convex linear complementary
problem are the only NP-complete optimization problems we have found so far
whose inverse problems are polynomial solvable. Are there any other NP-
complete optimization problems whose inverse problems are polynomial solvable?
The main difficulty is that there are no optimality conditions for most of the
NP-complete optimization problems.

Acknowledgement: This research was partially supported by National Science
Fund of China Under No. 19731010, No. 70171023, No. 79970052, and a fund from
the Institute of Policy and Management, Chinese Academy of Sciences.

REFERENCES

[1] Ahuja, R.K. and J.B. Orlin, Combinatorial algorithms for inverse network
flow problems, Working Paper, Sloan School of Management, MIT, Cambridge, MA,
1998.

[2] Ahuja, R.K. and J.B. Orlin, A fast algorithm for the inverse spanning tree
problem, J. of Algorithms, 34, 177-193(2000).

[3] Burton, D. and Ph.L. Toint, On an instance of the inverse shortest paths
problem, Mathematical Programming, 53, 45-61(1992).

[4] Burton, D. and Ph.L. Toint, On the use of an inverse shortest paths
algorithm for recovering linearly correlated costs, Mathematical Programming,
63, 1-22(1994).

[5] Burton, D., B. Pulleyblank and Ph. L. Toint, The inverse shortest paths
problem with upper bounds on shortest paths costs, In Network Optimization,
edited by P. Pardalos, D.W. Hearn and W.H. Hager, Lecture notes in Economics
and Mathematical Systems, Volumn 450, pp. 156-171(1997).

[6] Diao, Z. and M. Ding, Models and Algorithms for inverse convex quadratic
programming problem, Chinese J. of OR, Vol. 4(4)(2000), 88-94.

[7] Heuberger, C., Inverse combinatorial optimization: A survey on problems,
methods, and results, J. of Combinatorial Optimization, 8(3), 329-361, 2003.

[8] Horst, R., P.M. Pardalos and N.V. Thoai, Introduction to Global
Optimization, Kluwer Academic Publisher, 1995.

[9] Huang, S. and Z. Liu, On the inverse problem of linear programming and its
application to minimum weight perfect k-matching, European Journal of
Operational Research, 112, 421-426(1999).

[10] Huang, S. and Z. Liu, On the inverse minimum cost flow problem, Advances
in Operations Research and Systems Engineering, World Publishing, Co., 30-37,
1998.

[11] Huang, S., Inverse Problems of some NP-complete Problems, Lecture Notes
in Computer Sciences, LNCS 3521, N. Megiddo, Y. Xu and B. Zhu (eds).,
Springer-Verlag, Berlin, Heidelberg, 422-429, 2005.

[12] Garry, M. R. and D. S. Johnson, COMPUTERS AND INTRACTABILITY, A Guide to
the Theory of NP-Complete, W. H. Freeman and Company, 1979.

[13] Ji, J., F. A. Potra and S. Huang, Predictor-corrector method for linear
complementarity problems with polynomial complexity and superlinear
convergence, J. of Optimization Theory and Applications, 84, 187-199, 1995.

[14] Kojima, M., S. Mizuno and A. Yoshise, A polynomial-time algorithm for a
class of linear complementarity problems, Math. Programming, 44, 1-26, 1989.

[15] Monteiro, R. D. C. and I. Adler, Interior path following primal-dual
algorithms: Part II: Convex quadratic programming, Math. Programming, 44,
43-66, 1989.

[16] Murty, K. G. and S. N. Kabadi, Some NP-complete problems in quadratic and
nonlinear programming, Math. Programming, 39: 117-129(1987).

[17] Papadimitriou, C.H. and K. Steiglitz, Combinatorial Optimization:
Algorithms and Complexity, Prentice-Hall, Inc., Englewood Cliffs, New Jersey,
1982.

[18] Sokkalingam, P.T., R.K. Ahuja and J.B Orlin, Solving inverse spanning
tree problems through network flow techniques, Operations Research, 47,
291-298(1999).

[19] Zhang, J., Z., Liu and Z., Ma, On the inverse problem of minimum spanning
tree with partition constraints, Mathematical Methods of Operations Research,
44, 171-188(1996).

[20] Zhang, J., Z., Liu, Calculating some inverse linear programming problems,
J. Comput. Appl. Math., 72(1996), 261-273.


# Synthetic Data:
Can We Trust Statistical Estimators?

###### Abstract

The increasing interest in data sharing makes synthetic data appealing.
However, the analysis of synthetic data raises a unique set of methodological
challenges. In this work, we highlight the importance of inferential utility
and provide empirical evidence against naive inference from synthetic data
(that handles these as if they were really observed). We argue that the rate
of false-positive findings (type 1 error) will be unacceptably high, even when
the estimates are unbiased. One of the reasons is the underestimation of the
true standard error, which may even progressively increase with larger sample
sizes due to slower convergence. This is especially problematic for deep
generative models. Before publishing synthetic data, it is essential to
develop statistical inference tools for such data.

Alexander Decruyenaere°*†, Heidelinde Dehaene°*†, Paloma
Rabaey$\mathsection$†, Christiaan Polet*†,
Johan Decruyenaere*†, Stijn Vansteelandt†, Thomas Demeester$\mathsection$†

*Ghent University Hospital – SYNDARA $\mathsection$imec †Ghent University

<EMAIL_ADDRESS>
##  1 Introduction

In recent years, the use of synthetic data as a privacy-preserving substitute
for real sensitive data has been gaining attention (Drechsler, 2011; Raab et
al., 2016; Raghunathan, 2021). Synthetic data are artificial data that mimic
the original data in terms of statistical properties, without revealing any
individual records (Chen et al., 2021). As such, synthetic data might be able
to replace the original data in statistical analysis, while preserving the
privacy of the individuals of the original dataset and thereby fulfilling the
regulatory privacy constraints. This is especially appealing for medical
research (Jordon et al., 2022; Zhang et al., 2020). Synthetic data can be
generated using statistical techniques or deep learning (DL) approaches such
as GANs and VAEs (Hernandez et al., 2022).

Statistical inference is typically based on $\sqrt{N}$-consistency and
asymptotic normality (Boos and Stefanski, 2013). In this work, we empirically
provide evidence against naive inference from synthetic data, whereby
synthetic data are treated as if they were actually observed. First, the extra
uncertainty coming from the fact that the data were generated by a predictive
model themselves should be acknowledged. Otherwise, this will yield estimators
with standard errors (SEs) that are too optimistically small (Raab et al.,
2016). Second, the regularisation bias inherent to data-adaptive (DL)
techniques (i.e., their bias-variance trade-off is optimised with respect to
the prediction error instead of the error in the estimator) may diminish too
slowly as the sample size grows larger, causing excess variability that is
difficult to systematically account for (Brain and Webb, 1999; Hines et al.,
2022). This may deliver estimators that are biased and/or whose SE converges
slower than $1/\sqrt{N}$.

The remainder of this paper is structured as follows. In Section 2, we
introduce a framework that allows us to empirically assess the behaviour of
statistical estimators when it is (incorrectly) assumed that synthetic data
can be treated like original data. We consider a setting with low-dimensional
tabular data that is commonly featured in applied research. Section 3
discusses the statistical properties of the estimators and addresses how to
correct for some of the added variability inherent to synthetic data in the
model-based SE, as originally proposed by Raab et al. (2016). Section 4
presents the results of our simulation study, draws attention to the
difference in inferential utility between synthetic data generated through
statistical vs. DL techniques, and elaborates on why the corrected SE fails
with DL techniques. Finally, Section 5 provides some concluding remarks and
hints at future work.

##  2 Experimental Setup

Figure 1: Visual representation of our experimental setup.

Figure 1 shows the steps that are taken in order to assess the finite sample
performance of various statistical estimators in the context of synthetic
data. The code to reproduce our simulation study is available on
https://github.com/syndara-lab/inferential-utility-workshop.

#### Step 1: Ground truth data distribution

We simulate low-dimensional tabular toy data sampled from an arbitrary ground
truth population. A generic setting with five variables is considered: _age_
(continuous), _disease stage_ (ordinal with four categories), _biomarker_
(continuous), _therapy_ (binary), and _death_ (binary). The dependency
structure of these variables is encoded by a directed acyclic graph (DAG),
shown in Figure 1. For details on the distributions of the variables, we refer
to Appendix A.1.

#### Step 2: Sample original data $\bm{D_{orig}}$

In a Monte Carlo simulation, we randomly draw $N$ samples from the population,
which form an original dataset $D_{orig}^{i}$. To study the behaviour of
statistical estimators with increasing sample sizes, we vary $N$ log-uniformly
between 50 and 5000 ($N$ $\in\\{50,160,500,1600,5000\\}$). For each value of
$N$, 200 Monte Carlo simulation runs are used.

#### Step 3: Train generative model $\bm{G}$

Each dataset $D^{i}_{orig}$ of size $N$ is used to train several types of
generative models. The considered state-of-the-art generative models fall into
two categories: statistical and DL approaches (Hernandez et al., 2022). The
statistical approaches we investigate include Synthpop (Nowok et al., 2016) (R
package which facilitates fitting a series of conditional predictive models,
the type of which is inferred from the data, the order provided by the user)
and Bayesian network with user-specified DAG (BN) (Ankan and Panda, 2015). For
the DL approaches, we consider two popular deep generative models adapted for
tabular synthetic data, namely CTGAN and TVAE (Xu et al., 2019). Additional
information on these methods can be found in Appendix A.2.

#### Step 4: Generate synthetic data $\bm{D_{synth}}$

Using the trained generative model $G$, we generate a single synthetic dataset
$D_{synth}^{i}$ of size $M$. In our study, we set $M$ = $N$ to retain the
sample size of the original dataset, to test whether the user would get
similar results on both. For a particular choice of $N$, this leaves us with
200 synthetic datasets for each of the four generator types.

#### Step 5: Assess statistical properties of estimators based on
$\bm{D_{synth}}$

The last step concerns the main contribution of this work and is explained in-
depth in the next section.

##  3 Evaluating Statistical Properties based on Synthetic Data

We assess the behaviour of some standard statistical estimators when it is
(incorrectly) assumed that $D_{synth}$ can be treated like $D_{orig}$. These
estimators are estimated on each synthetic dataset $D_{synth}^{i}$, resulting
in an estimate of the population parameter in question. For each generator
type, we then have 200 estimates per sample size $N$, for a variety of sample
sizes $N$, allowing us to evaluate their empirical bias (i.e., the average
difference between the estimates and the population parameter) and empirical
SE (i.e., the standard deviation of the estimates). Standard statistical
analysis assumes that the bias converges faster than the SE with the latter
diminishing at a rate of $1/\sqrt{N}$. The corrected SE discussed later in
this section, though taking into account some of the added variability of the
synthetic data generating process, still relies on the same assumption, thus
remaining invalid in the context of deviating convergence rates.

As a result of the excess variability induced by the generation process that
may, in addition, decrease slower than $\sqrt{N}$, we foresee atypical
behaviour of the estimators (which will be confirmed by our simulation results
in Section 4). This impacts the type 1 error rate (the probability to find a
significant effect when in truth there is none) of statistical analyses
performed on synthetic data.

Motivated by commonly used estimators in applied research, our analysis in
Section 4 initially focuses on two statistical estimators: sample mean of
_age_ (“ _mean age_ ”) and a coefficient (“ _effect age on death_ ”) of a
multiple logistic regression model with _death_ as outcome and _age_ , _stage_
and _therapy_ as covariates. We report results for the other coefficients of
this regression model (effect of _therapy_ and _stage_ on _death_) in Appendix
A.3.

As will be shown empirically in Section 4, the model-based SE of the
considered estimators systematically underestimates the empirical SE by
ignoring the extra variability caused by the synthetic data generation
process. To partially account for this added variability, we will use the
following SE estimator proposed by Raab et al. (2016) and adapted to our
setting:

$\sigma_{\hat{\theta}\text{, corrected}}=\sigma_{\hat{\theta}\text{,
naive}}\sqrt{1+\frac{M}{N}}$

where $\sigma_{\hat{\theta}\text{, naive}}$ is the model-based SE in the
synthetic data, $M$ the sample size of the synthetic data $D_{synth}$, and $N$
the sample size of the original data $D_{orig}$. Note that this is a minimal
correction since it only applies to $\sqrt{N}$-consistent estimators, hence
the added variability resulting from the regularisation bias in data-adaptive
techniques (i.e., DL approaches), is not accounted for.

##  4 Simulation study

Figure 2: Simulation study results for the two estimators. Each dot is an
estimate per Monte Carlo run (200 in total per $N$). The population parameter
is represented by the horizontal dashed line. The funnel indicates the
behaviour of an unbiased and $\sqrt{N}$-consistent estimator based on observed
data.

#### Bias and SE

The estimates are depicted in Figure 2 for the two estimators and for all four
generator types. This figure allows a qualitative and empirical assessment of
two key properties of statistical estimators: bias and SE. The funnel
indicates the behaviour of an unbiased and $\sqrt{N}$-consistent estimator
based on original data. Table 1 lists the same information numerically,
summarising the relative bias ($\text{RE}_{\hat{\theta}}$) and the relative
underestimation of the empirical SE by the naive model-based SE
($\text{RE}_{\hat{\sigma}_{\hat{\theta}}}$). Figure A1 and Table A1 in the
appendix show these metrics for all estimators.

First, generative model misspecification will introduce bias. This is, for
example, noticeable in the _effect age on death_ coefficient generated by BN,
which is incorrectly estimated due to discretisation of the continuous
variable $age$ during the generation of $death$. Interestingly, the DL
approaches (CTGAN and TVAE), despite being more flexible and also tuned to
prevent overfitting, fail to adequately fit the joint distribution in our
simulation study, leading to bias for several estimators (e.g., biased null
effect for _effect age on death_ by CTGAN).

Second, the empirical SEs are larger for $D_{synth}$ than for $D_{orig}$ and
may vary over generative models. Larger empirical SEs reflect the additional
(predictive) uncertainty in the generation of synthetic data, which seems most
pronounced with DL approaches. This uncertainty is discarded in the naive use
of model-based SEs, leading to underestimation of the empirical SE (as is
evident from Table 1).

Table 1:  Relative error (RE) for the two estimators, averaged over 200 Monte
Carlo runs. $\text{RE}_{\hat{\theta}}$ is the relative bias of the estimates
$\hat{\theta}$. $\text{RE}_{\hat{\sigma}_{\hat{\theta}}}$ is the RE between
the naive model-based SE ($\hat{\sigma}_{\hat{\theta},naive})$ and the
empirical SE. Positive and negative values indicate a relative over- and
underestimation.

| Synthpop | BN
---|---|---
| $\text{RE}_{\hat{\theta}}$ | $\text{RE}_{\hat{\sigma}_{\hat{\theta}}}$ | $\text{RE}_{\hat{\theta}}$ | $\text{RE}_{\hat{\sigma}_{\hat{\theta}}}$
Estimator | $N=50$ | $N=5000$ | $N=50$ | $N=5000$ | $N=50$ | $N=5000$ | $N=50$ | $N=5000$
Mean age | 0.20 | 0.01 | -36.21 | -31.23 | 2.59 | -0.10 | -6.17 | -2.32
Effect age on death | 30.55 | 0.83 | -38.11 | -28.42 | -41.76 | -57.81 | -9.98 | -3.13
| CTGAN | TVAE
| $\text{RE}_{\hat{\theta}}$ | $\text{RE}_{\hat{\sigma}_{\hat{\theta}}}$ | $\text{RE}_{\hat{\theta}}$ | $\text{RE}_{\hat{\sigma}_{\hat{\theta}}}$
Estimator | $N=50$ | $N=5000$ | $N=50$ | $N=5000$ | $N=50$ | $N=5000$ | $N=50$ | $N=5000$
Mean age | 0.21 | 0.13 | -51.19 | -79.15 | 1.81 | 0.13 | -48.92 | -75.99
Effect age on death | -100.83 | -98.41 | 24.45 | -29.73 | -29.65 | -37.95 | -18.03 | -53.96

Table 2:  Estimated exponent $a$ [$95\%$ CI] for the power law convergence
rate $N^{-a}$ for empirical SE. Deviations from $0.50$ indicate that the
estimator is not $\sqrt{N}$-consistent.

| Mean age | Effect age on death
---|---|---
Original | 0.50 [0.49; 0.50] | 0.49 [0.47; 0.52]
Synthpop | 0.45 [0.43; 0.47] | 0.43 [0.36; 0.49]
BN | 0.50 [0.49; 0.50] | 0.50 [0.47; 0.52]
CTGAN | 0.29 [0.25; 0.32] | 0.53 [0.47; 0.58]
TVAE | 0.39 [0.33; 0.46] | 0.40 [0.36; 0.44]

Table 3:  Type 1 error rates (%) of the one-sample t-test at $\alpha=5\%$ for
the population mean of $age$ with naive model-based SEs and corrected SEs. The
null hypothesis states that the population mean of $age$ is equal to the
ground truth.

| Naive SE | Corrected SE
---|---|---
| $N=50$ | $N=5000$ | $N=50$ | $N=5000$
Synthpop | 20.0 | 15.0 | 9.0 | 6.5
BN | 16.5 | 6.5 | 5.5 | 0.5
CTGAN | 39.5 | 78.0 | 15.0 | 68.0
TVAE | 34.5 | 62.5 | 19.0 | 49.0

#### Convergence rate

Assuming a power law $N^{-a}$ in convergence rate for the empirical SE, we
estimated the exponent $a$ from five logarithmically spaced sample sizes $N$
between 50 and 5000, shown for the two estimators in Table 3 and all
estimators in Table A2 and Figure A2 in the appendix. As expected, the
empirical SE of estimators based on original data indeed converges at a rate
of $1/\sqrt{N}$ (i.e., $a_{SE}=0.5$). Fully parametric generative models are
also expected to yield estimators of which the SE decreases at a rate of
$1/\sqrt{N}$. This is largely confirmed for the statistical approaches we
studied, especially BN, though Synthpop converges a bit slower (i.e.,
$a_{SE}<0.5$), as non-parametric components are built into its generation
process. By contrast, the SEs produced by the more data-adaptive DL approaches
converge much slower (i.e. $a_{SE}<<0.5$), except for the logistic regression
coefficients that are estimated as (biased) null effects on synthetic data
generated by CTGAN. The slower-than-$\sqrt{N}$-convergence leads to a
progressively increasing underestimation of the empirical SE by the model-
based SE (which assumes $\sqrt{N}$-convergence) as the sample size grows
larger, as seen in Table 1.

The convergence rate of the bias (or equivalently, convergence rate of the
estimator) was estimated analogously and is presented in Table A2 in the
appendix. As opposed to default behaviour, the bias may converge slower than
the SE ($a_{bias}\leq a_{SE}$) for some estimators and generators.

#### Inferential utility

The null hypothesis significance testing framework typically uses an estimate
and its associated (un)certainty (reflected by the SE) as a test statistic to
reject a null hypothesis. Due to (1) slower convergence of (the bias of) the
estimator compared to the SE and/or (2) underestimation of the empirical SE by
the naive model-based SE, we indeed observe an inflation of the type 1 error
rate beyond the nominal level of $\alpha=5\%$, as shown in Table 3 and Figure
A3 in the appendix: the more the empirical SE is underestimated by the naive
model-based SE (which is especially the case for the DL approaches with
increasing sample size), the larger the inflation. Use of corrected SEs (as
introduced in Section 3) will control the type 1 error rate at approximately
$5\%$, but only for the (parametric) statistical approaches. By contrast, the
corrected SE does not account sufficiently for the predictive uncertainty of
the DL approaches due to (slower-than-$\sqrt{N}$-convergence of) their
regularisation bias, resulting in an actual type 1 error rate still higher
than the nominal level.

##  5 Conclusion

Naive inference from synthetic data leads to a flurry of false-positive
findings, even when the estimates are unbiased. One of the reasons is the
underestimation of the true standard error, which may even progressively
increase with larger sample sizes due to slower convergence. This is
especially problematic for deep generative models. Before publishing synthetic
data to enable data sharing and accelerate research, it is essential to
develop statistical inference tools for such data.

## Acknowledgements

This research was funded by a grant received from the Fund for Innovation and
Clinical Research of Ghent University Hospital. Paloma Rabaey’s research is
funded by the Research Foundation Flanders (FWO-Vlaanderen) with grant number
1170122N. This research also received funding from the Flemish government
under the “Onderzoeksprogramma Artificiële Intelligentie (AI) Vlaanderen”
programme.

## References

  * Ankan and Panda (2015) Ankur Ankan and Abinash Panda.  pgmpy: Probabilistic graphical models using python.  In _Proceedings of the 14th Python in Science Conference (SCIPY 2015)_ , 2015.
  * Boos and Stefanski (2013) Dennis D Boos and Leonard A Stefanski.  _Essential statistical inference: theory and methods_ , volume 591.  Springer, 2013.
  * Brain and Webb (1999) Damien Brain and Geoffrey I Webb.  On the effect of data set size on bias and variance in classification learning.  In _Proceedings of the Fourth Australian Knowledge Acquisition Workshop, University of New South Wales_ , pages 117–128, 1999.
  * Chen et al. (2021) Richard Chen, Ming Lu, Tiffany Chen, Drew Williamson, and Faisal Mahmood.  Synthetic data in machine learning for medicine and healthcare.  _Nature Biomedical Engineering_ , 5:1–5, 06 2021.
  * Drechsler (2011) Jörg Drechsler.  _Synthetic datasets for statistical disclosure control: theory and implementation_ , volume 201.  Springer Science & Business Media, 2011.
  * Hernandez et al. (2022) Mikel Hernandez, Gorka Epelde, Ane Alberdi, Rodrigo Cilla, and Debbie Rankin.  Synthetic data generation for tabular health records: A systematic review.  _Neurocomputing_ , 493:28–45, 2022.
  * Hines et al. (2022) Oliver Hines, Oliver Dukes, Karla Diaz-Ordaz, and Stijn Vansteelandt.  Demystifying statistical learning based on efficient influence functions.  _The American Statistician_ , 76:1–48, 01 2022.
  * Jordon et al. (2022) James Jordon, Lukasz Szpruch, Florimond Houssiau, Mirko Bottarelli, Giovanni Cherubin, Carsten Maple, Samuel N. Cohen, and Adrian Weller.  Synthetic data – what, why and how?, 2022.
  * Nowok et al. (2016) Beata Nowok, Gillian M Raab, and Chris Dibben.  synthpop: Bespoke creation of synthetic data in r.  _Journal of statistical software_ , 74:1–26, 2016.
  * Patki et al. (2016) Neha Patki, Roy Wedge, and Kalyan Veeramachaneni.  The synthetic data vault.  In _IEEE International Conference on Data Science and Advanced Analytics (DSAA)_ , pages 399–410, 2016.
  * Raab et al. (2016) Gillian M Raab, Beata Nowok, and Chris Dibben.  Practical data synthesis for large samples.  _Journal of Privacy and Confidentiality_ , 7(3):67–97, 2016.
  * Raghunathan (2021) Trivellore E Raghunathan.  Synthetic data.  _Annual review of statistics and its application_ , 8:129–140, 2021.
  * Xu et al. (2019) Lei Xu, Maria Skoularidou, Alfredo Cuesta-Infante, and Kalyan Veeramachaneni.  Modeling tabular data using conditional GAN.  _Advances in neural information processing systems_ , 32, 2019.
  * Zhang et al. (2020) Ziqi Zhang, Chao Yan, Diego A Mesa, Jimeng Sun, and Bradley A Malin.  Ensuring electronic medical record simulation through better training, modeling, and evaluation.  _Journal of the American Medical Informatics Association_ , 27(1):99–108, 2020.

##  Appendix A Details experimental setup

###  A.1 Ground truth data distribution

It was opted to work with low-dimensional tabular data given its frequent use
in applied medical research. Different traditional and commonly used
regression models were taken into account when choosing the nature of the
variables. We wish to have a mix of variables that are continuous, binary,
ordinal, normally distributed or skewed. To obtain these requirements and
reflect a generic clinical setting, the data generating process consists of
the following five variables: _age_ (continuous), _disease stage_ (ordinal
with four categories or stages), _biomarker_ (continuous), _therapy_ (binary),
and _death_ (binary). It is assumed that a patient is observed at a given
point in time, which is further referred to as the baseline time point. At
this time, patient data about _age_ , _disease stage_ , _biomarker_ , and the
random assignment of _therapy_ is gathered. The binary outcome variable
_death_ is evaluated at a later time point. This design is a simplification of
reality since we assume that there are no missing data and we do not consider
the data as longitudinal.

_Age_ follows a normal distribution with mean $50$ and standard deviation
$10$. _Disease stage_ was generated according to a proportional odds
cumulative logit model where an increase in _age_ causes an increase in the
odds of having a _disease stage_ higher than a given stage $k$,
$\nu_{age}=-0.05$. The variable _biomarker_ is a quantification of the
_disease stage_ and was also based on a Generalized Linear Model (GLM), where
_biomarker_ follows a gamma distribution and its mean changes in function of
_disease stage_. It was constructed in such a way that a higher _disease
stage_ results in higher values for the _biomarker_ ($\gamma_{0}=4$,
$\gamma_{stage}=\\{0,-1,-2,-3\\}$ for stage I-IV, respectively). _Therapy_ is
considered to be 1:1 randomly assigned and is therefore sampled from a
Bernouilli distribution with a probability of $0.50$. The last variable,
_death_ , is generated by using a binomial logistic regression model in which
the odds of _death_ increase with an increasing _age_ ($\beta_{age}=0.05$), a
higher _disease stage_ ($\beta_{stage}=\\{0,0.50,1.00,1.50\\}$ for stage I-IV,
respectively), and absence of _therapy_ ($\beta_{therapy}=-0.50$).

###  A.2 Generative models

The hyperparameters of CTGAN and TVAE were tuned using the Tree-structured
Parzen Estimator algorithm with the 5-fold cross-validated average inverse of
the Kullback–Leibler divergence (IKLD) between the original and the synthetic
dataset as objective score, averaged over five seed initializations. We use
the Python package SDV (Patki et al., 2016) to train our CTGAN and TVAE
models. Numerical features are preprocessed using a cluster-based normalizer.

Our _synthpop_ statistical approach relies on the synthpop package for R,
which provides a routine to generate synthetic data (Nowok et al., 2016). This
framework encompasses both parametric and non-parametric methods to
sequentially fit a series of conditional joint distributions, based on the
observed data. We restrict ourselves to the default parametric methods and
provide information of the dependency structure of our data via specification
of a directed acyclic graph (DAG). As a result, this defines the order of the
sequence and which variables need to be included as predictors in the
conditional models.

Our Bayesian Networks were implemented using the Python package pgmpy (Ankan
and Panda, 2015). The ground-truth DAG is provided upfront, while the
conditional probability distributions (CPDs) are estimated by Maximum
Likelihood Estimation (MLE) and synthetic data are generated via forward
sampling.

###  A.3 Full results simulation study

Figure A1 shows the simulation results for all estimators (both mean of _age_
and all logistic regression coefficients for the variable _death_). Table A1
summarises the relative bias ($\text{RE}_{\hat{\theta}}$) and the relative
underestimation of the empirical SE by the naive model-based SE
($\text{RE}_{\hat{\sigma}_{\hat{\theta}}}$) for all estimators. Table A2 show
the rate at which the empirical SE and bias decrease in function of $N$ on the
log-scale for the same set of estimators. Figure A2 visually compares the rate
of convergence of the SE across all generators.

Figure A1: Simulation study results for all estimators. Each dot is an
estimate per Monte Carlo run (200 dots in total per value of $N$). The
population parameter is represented by the horizontal dashed line. The funnel
indicates the behaviour of an unbiased and $\sqrt{N}$-consistent estimator
based on observed data. Table A1:  Relative error (RE) for all estimators,
averaged over 200 Monte Carlo runs. $\text{RE}_{\hat{\theta}}$ is the relative
bias of the estimates $\hat{\theta}$.
$\text{RE}_{\hat{\sigma}_{\hat{\theta}}}$ is the RE between the naive model-
based SE ($\hat{\sigma}_{\hat{\theta},naive})$ and the empirical SE. Positive
and negative values indicate a relative over- and underestimation.

| Synthpop | BN
---|---|---
| $\text{RE}_{\hat{\theta}}$ | $\text{RE}_{\hat{\sigma}_{\theta}}$ | $\text{RE}_{\hat{\theta}}$ | $\text{RE}_{\hat{\sigma}_{\theta}}$
Estimator | $N=50$ | $N=5000$ | $N=50$ | $N=5000$ | $N=50$ | $N=5000$ | $N=50$ | $N=5000$
Mean |  |  |  |  |  |  |  |
Mean age | 0.20 | 0.01 | -36.21 | -31.23 | 2.59 | -0.10 | -6.17 | -2.32
Logistic regression |  |  |  |  |  |  |  |
Effect age on death | 30.55 | 0.83 | -38.11 | -28.42 | -41.76 | -57.81 | -9.98 | -3.13
Effect therapy on death | 23.94 | -0.26 | -34.81 | -30.13 | -50.16 | -55.24 | -11.63 | -15.83
Effect stage II on death | 31.98 | -0.04 | -30.26 | -30.46 | -49.14 | -54.56 | -12.22 | -6.64
Effect stage III on death | 23.83 | 2.31 | -31.76 | -32.47 | -42.65 | -51.54 | -6.69 | 0.39
Effect stage IV on death | 7.29 | 1.15 | -25.42 | -36.90 | -45.20 | -53.25 | -9.68 | -16.36
| CTGAN | TVAE
| $\text{RE}_{\hat{\theta}}$ | $\text{RE}_{\hat{\sigma}_{\theta}}$ | $\text{RE}_{\hat{\theta}}$ | $\text{RE}_{\hat{\sigma}_{\theta}}$
Estimator | $N=50$ | $N=5000$ | $N=50$ | $N=5000$ | $N=50$ | $N=5000$ | $N=50$ | $N=5000$
Mean |  |  |  |  |  |  |  |
Mean age | 0.21 | 0.13 | -51.19 | -79.15 | 1.81 | 0.13 | -48.92 | -75.99
Logistic regression |  |  |  |  |  |  |  |
Effect age on death | -100.83 | -98.41 | 24.45 | -29.73 | -29.65 | -37.95 | -18.03 | -53.96
Effect therapy on death | -99.92 | -96.48 | 10.41 | -17.28 | -17.43 | -28.77 | -21.86 | -57.59
Effect stage II on death | -104.41 | -99.37 | 38.07 | 6.12 | -37.74 | -5.12 | -19.34 | -51.34
Effect stage III on death | -104.58 | -98.90 | 30.65 | 12.75 | -34.81 | -19.28 | -10.76 | -53.25
Effect stage IV on death | -104.90 | -99.79 | 25.92 | 10.90 | -41.73 | -21.86 | -9.56 | -63.30

Table A2: Estimated exponent $a$ [$95\%$ CI] for the power law convergence
rate $N^{-a}$ for empirical SE and bias of all estimators. Standard
statistical analysis assumes that the bias converges faster than the SE with
the latter diminishing at a rate of $1/\sqrt{N}$ ($a_{SE}=0.50$). Bold values
indicate instances where the bias converges slower than the SE for the same
generator and estimator.

|  | Generator
---|---|---
| Original | Synthpop | BN | CTGAN | TVAE
Estimator, SE |  |  |  |  |
Mean |  |  |  |  |
Mean age | 0.50 [0.49; 0.50] | 0.45 [0.43; 0.47] | 0.50 [0.49; 0.50] | 0.29 [0.25; 0.32] | 0.39 [0.33; 0.46]
Logistic regression |  |  |  |  |
Effect age on death | 0.49 [0.47; 0.52] | 0.43 [0.36; 0.49] |  0.50 [0.47; 0.52] |  0.53 [0.47; 0.58] | 0.40 [0.36; 0.44]
Effect therapy on death | 0.50 [0.48; 0.51] | 0.44 [0.39; 0.49] |  0.49 [0.47; 0.51] |  0.51 [0.47; 0.56] | 0.41 [0.38; 0.43]
Effect stage II on death | 0.49 [0.47; 0.50] | 0.44 [0.40; 0.48] |  0.49 [0.46; 0.51] |  0.52 [0.50; 0.55] | 0.41 [0.40; 0.43]
Effect stage III on death | 0.49 [0.48; 0.51] | 0.44 [0.39; 0.49] |  0.49 [0.46; 0.51] |  0.53 [0.51; 0.55] | 0.42 [0.40; 0.43]
Effect stage IV on death | 0.49 [0.48; 0.51] | 0.44 [0.41; 0.47] |  0.49 [0.47; 0.50] |  0.54 [0.51; 0.56] | 0.41 [0.38; 0.43]
Estimator, bias |  |  |  |  |
Mean |  |  |  |  |
Mean age | 0.93 [0.81; 1.04] | 0.93 [0.79; 1.07] | 0.68 [0.44; 0.93] | 0.59 [0.33; 0.85] | 0.65 [0.53; 0.78]
Logistic regression |  |  |  |  |
Effect age on death | 0.88 [0.66; 1.10] | 0.81 [0.65; 0.97] |  0.32 [0.18; 0.47] |  0.24 [0.15; 0.32] | 0.51 [0.20; 0.82]
Effect therapy on death | 0.91 [0.82; 0.99] | 0.91 [0.87; 0.96] |  0.41 [0.26; 0.57] |  0.33 [0.21; 0.44] | 0.51 [0.29; 0.73]
Effect stage II on death | 0.99 [0.74; 1.24] | 1.03 [0.85; 1.22] |  0.46 [0.29; 0.63] |  0.36 [0.23; 0.49] | 0.66 [0.54; 0.79]
Effect stage III on death | 0.90 [0.68; 1.13] | 0.84 [0.64; 1.04] |  0.37 [0.23; 0.52] |  0.27 [0.17; 0.37] | 0.47 [0.34; 0.60]
Effect stage IV on death | 0.84 [0.75; 0.92] | 0.80 [0.67; 0.93] |  0.32 [0.19; 0.45] |  0.22 [0.14; 0.29] | 0.42 [0.31; 0.52]

Figure A2: Convergence rate of the empirical standard error (SE). If the SE is
of the form SE $=cN^{-a}$, where $c$ is a constant, then
$\log{(SE)}=c+(-a)\log{N}$. Therefore slope $a$ represents the convergence
rate and the vertical offset $c$ indicates the asymptotic variance. The dashed
line indicates the behaviour of the SE of an unbiased and
$\sqrt{N}$-consistent estimator based on observed data, whereas the dotted
line indicates the assumed behaviour of the SE of the same estimator based on
synthetic data, following the correction proposed by Raab et al. (2016). Note
that the asymptotic variance of the logistic regression coefficients (“
_effect [predictor] on death_ ”) estimated on synthetic data generated by
CTGAN is smaller than on the observed data, as CTGAN delivers a biased null
effect due to generative model misspecification. Figure A3: Type 1 error rate
of the one-sample t-test at $\alpha=5\%$ for the population mean of $age$ with
naive model-based and corrected SEs. The null hypothesis states that the
population mean of $age$ is equal to the ground truth.