# Supplementary File for “Addressing patient heterogeneity in disease
predictive model development”

Xu Gao  Weining Shen  Jing Ning  Ziding Feng  and Jianhua Hu

##  1 Additional numerical results

We present additional simulation results from Section 4 of the main paper
here.

###  1.1 Additional simulation results for normal mixture models

We provide the tuning parameter $C$ selection results for normal mixture model
with $m_{0}=1,2,3$ in Tables S1–S3. In each table, the selected value for $C$
is highlighted in bold text, e.g., $C=3.0$ in Table S1; and this value is
chosen based on the type I error (closest to the nominal levels).

Table S1: Tuning parameter $C$ selection: rejection proportion for mixture-of-normal example with $m_{0}=1$. C | $\alpha=.01$ | $\alpha=.05$ | $\alpha=.1$
---|---|---|---
0.2 | .013 | .078 | .145
0.4 | .012 | .071 | .133
0.6 | .010 | .064 | .127
0.8 | .010 | .061 | .122
1.0 | .010 | .060 | .115
1.5 | .009 | .058 | .110
2.0 | .008 | .058 | .106
3.0 | .008 | .056 | .105
8.0 | .007 | .054 | .101
12.0 | .007 | .053 | .100
Table S2: Tuning parameter $C$ selection: rejection proportion for mixture-of-normal example with $m_{0}=2$. C | $\alpha=.01$ | $\alpha=.05$ | $\alpha=.1$
---|---|---|---
0.2 | .017 | .065 | .137
0.4 | .016 | .060 | .117
0.6 | .014 | .055 | .107
0.8 | .013 | .051 | .099
1.0 | .012 | .045 | .088
1.5 | .009 | .039 | .081
2.0 | .008 | .033 | .072
Table S3: Tuning parameter $C$ selection: rejection proportion for mixture-of-normal example with $m_{0}=3$. C | $\alpha=.01$ | $\alpha=.05$ | $\alpha=.1$
---|---|---|---
0.2 | .021 | .094 | .175
0.4 | .019 | .084 | .157
0.6 | .014 | .082 | .140
0.8 | .013 | .079 | .132
1.0 | .012 | .072 | .123
1.5 | .010 | .066 | .106
2.0 | .008 | .056 | .096
3.0 | .007 | .046 | .080
5.0 | .006 | .038 | .066

###  1.2 Results for logistic mixture models

Here we present the simulation results for mixture of logistic regression
models with $m_{0}=1,2,3$ here. We consider three simulation scenarios.

Scenario 1. We let the true number of subgroups $m_{0}=1$ and $\alpha=1$,
i.e., there is only one subgroup. The response variable $Y$ is generated from
a logistic regression model with two covariates $X_{1},X_{2}$ sampled from an
uniform distribution on $(5,10)$, and the coefficient is
$\bm{\theta}=(.4,.6)$.

Scenario 2. The response variable $Y$ is generated from a mixture of two
logistic regression models, with two subgroup-specific covariates
$X_{1},X_{2}$ sampled from an uniform distribution on $(5,10)$. Here
$\bm{\alpha}=(.6,.4)$, the subgroup-specific regression coefficients are
$\bm{\theta_{1}}=(-.1,-.2)$ and $\bm{\theta_{2}}=(.2,.1)$.

Scenario 3. We consider $m_{0}=3$ with weights $\bm{\alpha}=(.4,.3.,3)^{T}$.
The response variable $Y$ is generated from a mixture of three logistic
regression models, with $X_{1},X_{2}$ sampled from an uniform distribution on
$(5,10)$. We set the subgroup-specific regression coefficients as
$\bm{\theta_{1}}=(-.1,-.2)$, $\bm{\theta_{2}}=(0,.3)$,
$\bm{\theta_{3}}=(.3,0)$.

We let the sample size $n$ take the values of $500$, $1000$, and $1500$ and
examine the nominal levels $\alpha$ of $0.01$ and $0.05$, respectively. We
repeat the EM test procedure for three iterations ($K=3$) for each case. To
decide the tuning parameter $C$ in the penalty function $p(\cdot)$, we
generate a separate dataset under each simulation scenario, consider a range
of different values of $C$ and choose the best one that provides the type I
errors that are closest to the nominal values. Based on the results summarized
in Tables S4-S6, we choose $C=1.8$ for $m=1$, $C=1.0$ for $m=2$, and $C=2.0$
for $m=3$. From those tables we also find that the type I errors are quite
robust to different values of $C$.

We summarize the type-I error based on 5000 Monte-Carlo replications in the
first part of Table S7. It can be seen that type-I error rates are generally
close to the nominal levels across different scenarios, regardless of
different sample sizes and nominal levels.

Table S4: Tuning parameter $C$ selection: rejection proportion for mixture-of-logistic example with $m_{0}=1$. C | $\alpha=.01$ | $\alpha=.05$ | $\alpha=.1$
---|---|---|---
0.4 | .019 | .079 | .151
0.8 | .014 | .066 | .122
1.4 | .013 | .055 | .106
1.8 | .013 | .050 | .101
2.2 | .013 | .048 | .098
3.0 | .012 | .045 | .089
Table S5: Tuning parameter $C$ selection: rejection proportion for mixture-of-logistic example with $m_{0}=2$. C | $\alpha=.01$ | $\alpha=.05$ | $\alpha=.1$
---|---|---|---
0.1 | .013 | .067 | .131
0.5 | .013 | .062 | .112
.9 | .012 | .051 | .096
1.0 | .012 | .049 | .095
1.5 | .010 | .045 | .083
2.0 | .010 | .036 | .078
Table S6: Tuning parameter $C$ selection: rejection proportion for mixture-of-logistic example with $m_{0}=3$. C | $\alpha=.01$ | $\alpha=.05$ | $\alpha=.1$
---|---|---|---
0.2 | .015 | .097 | .202
0.6 | .013 | .081 | .166
1.0 | .008 | .054 | .120
1.5 | .006 | .045 | .108
1.9 | .006 | .043 | .101
2.0 | .005 | .042 | .099
Table S7: Type-I errors and power for various sample sizes and nominal levels based on 5000 simulations for mixture-of-logistic example. | Type I error
---|---
| $m_{0}=1$ | $m_{0}=2$ | $m_{0}=3$
$n$ | $\alpha=.01$ | $\alpha=.05$ | $\alpha=.01$ | $\alpha=.05$ | $\alpha=.01$ | $\alpha=.05$
$500$ | .007 | .033 | .012 | .049 | .003 | .034
$1000$ | .009 | .0375 | .012 | .067 | .009 | .048
$1500$ | .013 | .050 | .013 | .075 | .010 | .062
| Power (nominal level $=.1$)
| $m_{0}=1$ | $m_{0}=2$ | $m_{0}=3$
$n$ | weak | strong | weak | strong | weak | strong
$500$ | .442 | .472 | .172 | .367 | .121 | .131
$1000$ | .608 | .686 | .273 | .496 | .224 | .251
$1500$ | .681 | .834 | .455 | .620 | .277 | .354

To evaluate the power, for each simulation scenario, we consider two cases
that represent different levels of heterogeneity in terms of the parameter
value differences between the subgroups.

  * •

For Scenario 1, we consider two alternatives, a “weak” case, with weights
$(.5,.5)$ and coefficients $\bm{\theta_{1}}=(-.4,-.2),\bm{\theta_{2}}=(0,0)$;
and a “strong” case, with weights $(.9,.1)$ and coefficients
$\bm{\theta_{1}}=(-.4,-.2),\bm{\theta_{2}}=(.2,.4)$.

  * •

For Scenario 2, we consider two alternatives, a “weak” case, with weights
$(.1,.1,.1,.7)$ and coefficients
$\bm{\theta_{1}}=(-.2,-.4),\bm{\theta_{2}}=(-.2,0),\bm{\theta_{3}}=(-.2,.2),\bm{\theta_{4}}=(.4,.2)$;
and a “strong” case, with weights $(.25,.25,.25,.25)$, and coefficients
$\bm{\theta_{1}}=(-.1,-.2),\bm{\theta_{2}}=(-.1,0),\bm{\theta_{3}}=(1,0),\bm{\theta_{4}}=(.2,.2)$.

  * •

For Scenario 3, we consider two alternatives, a “weak” case, with weights
$(.2,.2,.2,.2,.1,.1)$ and coefficients
$\bm{\theta_{1}}=(-.1,-.2),\bm{\theta_{2}}=(0,.3),\bm{\theta_{3}}=(.3,0),\bm{\theta_{4}}=(.1,-.2),\bm{\theta_{5}}=(-.2,.2),\bm{\theta_{6}}=(-.3,0)$;
and a “strong” case, with weights $(.1,.2,.2,.2,.2,.1)$, and coefficients
$\bm{\theta_{1}}=(-.1,-.2),\bm{\theta_{2}}=(.1,-.2),\bm{\theta_{3}}=(-.3,.1),\bm{\theta_{4}}=(-.1,.1),\bm{\theta_{5}}=(0,.3),\bm{\theta_{6}}=(.2,.3)$.

The remaining data generation, including response and covariate distributions,
is the same as the setting of type-I error examination. We summarize the
proportion of rejecting the null hypothesis based on $1000$ simulations in the
second part of Table S7. It can be seen that the power of the test increases
up rapidly for $m_{0}=1$ and $2$, but does not perform as well for $m_{0}=3$.
There are several possible explanations. In general, the mixture of logistic
regression model is more difficult to fit compared to mixture models with a
continuous outcome, and hence requires a much larger sample size to achieve a
comparable performance with the normal mixture model. When evaluating the
power for $m_{0}=3$, the data is simulated from six subgroups, which makes the
sample size for each subgroup very small. Therefore it becomes more
challenging to detect the difference across different subgroups. Another
possible explanation is due to the violation of Assumption (C8) as discussed
in Section 2 of the Web Supplementary File.

###  1.3 Comparison with Shen2015’s method

We compare the performance of our EM test with the method proposed in Shen2015
in Table S8 and S9 here.

Table S8: Type-I errors for our method and Shen2015’s method. | Our method | Shen2015
---|---|---
$n$ | $\alpha=.01$ | $\alpha=.05$ | $\alpha=.1$ | $\alpha=.01$ | $\alpha=.05$ | $\alpha=.1$
$60$ | .011 | .040 | .076 | .012 | .045 | .094
$100$ | .010 | .045 | .088 | .011 | .050 | .100
Table S9: Power for for our method and Shen2015’s method with nominal level $=.05$. $n$ | $a$ | $b$ | $c$ | Our method | Shen2015
---|---|---|---|---|---
60 | 0.5 | 1 | 1 | .346 | .730
| 0.5 | 0 | 1 | .386 | .186
| 0.5 | 1 | 0 | .336 | .302
| 1.0 | 1 | 1 | .490 | .820
| 1.0 | 0 | 1 | .691 | .504
| 1.0 | 1 | 0 | .448 | .454
100 | 0.5 | 1 | 1 | .525 | .970
| 0.5 | 0 | 1 | .584 | .376
| 0.5 | 1 | 0 | .506 | .420
| 1.0 | 1 | 1 | .704 | .980
| 1.0 | 0 | 1 | .887 | .700
| 1.0 | 1 | 0 | .659 | .654

###  1.4 Comparison with tree method

|  |
---|---|---
Figure S1: Classification tree results for Scenario (i)–(iii). Table S10: Type I error for our method under the simulation scenarios 1–3 in Section 4.2 of the main paper based on 1000 replications. Scenario | $\alpha=.01$ | $\alpha=.05$ | $\alpha=.1$
---|---|---|---
1 | .013 | .051 | .111
2 | .024 | .060 | .086
3 | .012 | .049 | .093

##  2 Assumptions

We state the assumptions that are used for the asymptotic results in Section 3
of the main paper. For simplicity, let $D=(Y,\bm{X},\bm{Z})$,
$D_{i}=(Y_{i},\bm{X}_{i},\bm{Z}_{i})$. We assume the following regularity
conditions,

  * (C1)

We assume a set of conditions on $f$ for consistency of the MLE.

    * (a)

$\text{E}|\log f(D;\Psi_{0},\bm{\gamma}_{0})|<\infty$ for any values of $D$.

    * (b)

There exist compact sets $\Theta$ and $\Gamma$ such that
$\text{E}\log\\{1+\sup_{\|\bm{\theta}-\bm{\theta}_{0}\|\leq\epsilon,\|\bm{\gamma}-\bm{\gamma}_{0}\|\leq\epsilon}f(D;\bm{\theta},\bm{\gamma})\\}<\infty$
for any $\bm{\theta}_{0}\in\Theta,\bm{\gamma}_{0}\in\Gamma$ and a sufficiently
small $\epsilon>0$.

    * (c)

$\text{E}\log\\{1+\sup_{\|\bm{\theta}\|_{2}\geq M,\|\bm{\gamma}\|_{2}\geq
M}f(D;\bm{\theta},\bm{\gamma})\\}<\infty$ for a sufficiently large $M$.

    * (d)

$\lim_{\|\bm{\theta}\|_{2},\|\bm{\gamma}\|_{2}\rightarrow\infty}f(D;\bm{\theta},\bm{\gamma})=0$
almost surely.

  * (C2)

Assume that $f$ is four times continuously differentiable with respect to
$\bm{\theta}$ and $\bm{\gamma}$.

  * (C3)

For any two mixting distributions $\Psi_{1}$ and $\Psi_{2}$, if

$\displaystyle\int f(D;\bm{\theta},\bm{\gamma})d\Psi_{1}(\bm{\theta})=\int
f(D;\bm{\theta},\bm{\gamma})d\Psi_{2}(\bm{\theta})$

for any $D$, then $\Psi_{1}=\Psi_{2}$.

  * (C4)

There exists a small $\epsilon>0$ and an integrable function $g$ such that for
$h=1,\ldots,m_{0}$, and every $\bm{\theta}$ satisfying
$\|\bm{\theta}-\bm{\theta}_{0h}\|\leq\epsilon$,

$\displaystyle|\Delta_{ih}|^{3}\leq g(\bm{X}_{i}),~{}|u|^{3}\leq
g(\bm{X}_{i}),$ $\displaystyle|v|^{3}\leq
g(\bm{X}_{i}),~{}\left|\frac{\partial{v}}{\partial\theta_{k}}\right|^{3}\leq
g(\bm{X}_{i}),~{}\left|\frac{\partial{v}^{2}}{\partial\theta_{k}\partial\theta_{k}^{\prime}}\right|^{3}\leq
g(\bm{X}_{i}),$

where $u$ is any element in $\bm{Y}_{i}(\bm{\theta})$, $v$ is any element in
$\bm{Z}_{i}(\bm{\theta})$, and $k,k^{\prime}$ are indexes from $1,\ldots,p$.

  * (C5)

The covariance matrix of $\bm{b}_{i}$, denoted by $\bm{B}$, is finite and
positive definite.

  * (C6)

The penalty function $p(\beta)$ is continuous and satisfies $p(\beta)\leq 0$
for $\beta\in[0,1]$, $p(.5)=0$, $p(0)=p(1)=-\infty$.

  * (C7)

For every $k=1,\ldots,p$, there exists a constant $C>0$ and an index
$j\in\\{1,\ldots,m_{0}\\}$ such that
$|\bm{\theta}_{0jk}-\bm{\theta}_{01k}|>C$.

  * (C8)

We assume
$\frac{\partial^{2}f}{\partial\theta_{i}\theta_{j}}=o_{p}(\frac{\partial^{2}f}{\partial\theta_{i}^{2}})$
for every $i,j=1,\ldots,p$ and $i\neq j$.

To simplify the notation, we drop the subscript $n$ from $\lambda_{n}$
throughout the proof. Conditions (C1)–(C6) are similar with those in Li2010
with modifications to accommodate the regression coefficients $\bm{\theta}$
and $\bm{\gamma}$. Condition (C7) essentially rules out the situation where a
covariate $x_{k}$ has the same effect on $y$ in every subgroup.  Condition
(C8) essentially assumes the cross-derivatives of the likelihood function
$\frac{\partial^{2}f}{\partial\theta_{i}\theta_{j}}$ are $0$ for $i\neq j$.
This is satisfied for a normal linear mixture model where the covariates are
uncorrelated with each other, which can hold in our case given that one can
first apply a PCA for the covariate space and then perform the usual linear
regression on the PC components. For the logistic regression model, this
condition is satisfied given covariates are uncorrelated and the predicted
individual probabilities $\hat{p}_{i}(1-\hat{p}_{i})$ are the same for all the
patients $i$ in the same subgroup. This condition seems reasonable given that
the variability is assumed to be homogeneous for each subgroup. It is unclear
how to relax/generalize those assumptions under the current framework, but the
simulation results seem to suggest a reasonable level of robustness in our
method under the violation of those assumptions. For example, for the logistic
mixture model (Section 1.2 of the Web Supplementary File), we find that the
type I error is in general very close to the nominal level for most simulation
scenarios, although there seems a loss of power especially when $m_{0}=2$ or
$3$ due to the violation of this assumption.

##  3 Proof of Theorems 1–3

###  3.1 Proof of Theorem 1

The proof extends Theorem 1 of Li2010 to allow multivariate $\bm{X}$. The main
idea is to use mathematical induction. In Lemmas 1 and 2, we show that Theorem
1 holds for $k=1$. In Lemma 3, we show that if Theorem holds for $k=k_{0}$,
then it also holds for $k=k_{0}+1$. The proof of these lemmas goes to the next
section.

######  Lemma 1.

Suppose the Conditions (C1)–(C3), (C6) and (C7) hold. Under the null
distribution $f(D;\Psi_{0},\bm{\gamma}_{0})$, the following holds for every
$\bm{\beta}_{0}\in\mathcal{B}^{m_{0}}$,

$\displaystyle\bm{\alpha}^{(1)}-\bm{\alpha}_{0}=o_{p}(1),~{}\bm{\beta}^{(1)}-\bm{\beta}_{0}=o_{p}(1),~{}\bm{\theta}_{1}^{(1)}-\bm{\theta}_{0}=o_{p}(1),~{}\bm{\theta}_{2}^{(1)}-\bm{\theta}_{0}=o_{p}(1),~{}\bm{\gamma}^{(1)}-\bm{\gamma}_{0}=o_{p}(1).$

######  Lemma 2.

Suppose the Conditions (C1)–(C8) hold. Under the null distribution
$f(D;\Psi_{0},\bm{\gamma}_{0})$, the following holds for every
$\bm{\beta}_{0}\in\mathcal{B}^{m_{0}}$ ,

$\displaystyle\bm{\alpha}^{(1)}-\bm{\alpha}_{0}=O_{p}(n^{-1/2}),~{}\bm{\beta}^{(1)}-\bm{\beta}_{0}=O_{p}(n^{-1/6}),~{}\bm{m}_{1}^{(1)}=O_{p}(n^{-1/2}),$
$\displaystyle\bm{\theta}_{1}^{(1)}-\bm{\theta}_{0}=O_{p}(n^{-1/4}),~{}\bm{\theta}_{2}^{(1)}-\bm{\theta}_{0}=O_{p}(n^{-1/4}),~{}\bm{\gamma}^{(1)}-\bm{\gamma}_{0}=O_{p}(n^{-1/2}).$

######  Lemma 3.

Suppose that Conditions (C1)–(C8) hold. Under the null distribution
$f(D;\Psi_{0},\bm{\gamma}_{0})$, the following holds for every
$\bm{\beta}_{0}\in\mathcal{B}^{m_{0}}$ and some $k\geq 1$,

$\displaystyle\bm{\alpha}^{(k)}-\bm{\alpha}_{0}=O_{p}(n^{-1/2}),~{}\bm{\beta}^{(k)}-\bm{\beta}_{0}=O_{p}(n^{-1/6}),~{}\bm{m}_{1}^{(k)}=O_{p}(n^{-1/2}),$
$\displaystyle\bm{\theta}_{1}^{(k)}-\bm{\theta}_{0}=O_{p}(n^{-1/4}),~{}\bm{\theta}_{2}^{(k)}-\bm{\theta}_{0}=O_{p}(n^{-1/4}),~{}\bm{\gamma}^{(k)}-\bm{\gamma}_{0}=O_{p}(n^{-1/2}).$

then

$\displaystyle\bm{\alpha}^{(k+1)}-\bm{\alpha}_{0}=O_{p}(n^{-1/2}),~{}\bm{\beta}^{(k+1)}-\bm{\beta}_{0}=O_{p}(n^{-1/6}),~{}\bm{m}_{1}^{(k+1)}=O_{p}(n^{-1/2}),$
$\displaystyle\bm{\theta}_{1}^{(k+1)}-\bm{\theta}_{0}=O_{p}(n^{-1/4}),~{}\bm{\theta}_{2}^{(k+1)}-\bm{\theta}_{0}=O_{p}(n^{-1/4}),~{}\bm{\gamma}^{(k+1)}-\bm{\gamma}_{0}=O_{p}(n^{-1/2}).$

###  3.2 Proof of Theorem 2

We drop $\bm{\gamma}_{0}$ from $f$ for convenience. The proof contains two
parts. We first show that $EM_{n}^{(K)}\leq\sup_{\bm{v}\geq
0}\\{2\bm{v}^{T}\sum_{i=1}^{n}\tilde{\bm{b}}_{2i}-n\bm{v}^{T}\tilde{\bm{B}}_{22}\bm{v}\\}+o_{p}(1)$.
Then in the second part, we show the upper bound can be attained.

(1) Consider the expansion of $l_{n}(\hat{\Psi}_{0})$ at $\Psi_{0}$:

$\displaystyle
R_{0n}=2\\{l_{n}(\hat{\Psi}_{0})-l_{n}(\Psi_{0})\\}=(\sum_{i=1}^{n}\bm{b}_{1i})^{T}(n\bm{B}_{11})^{-1}(\sum_{i=1}^{n}\bm{b}_{1i})+o_{p}(1).$

Define

$\displaystyle
R_{1n}(\Psi^{(k)}(\bm{\beta}_{0}))=2\\{pl_{n}(\Psi^{(k)}(\bm{\beta}_{0}))-pl_{n}^{0}(\Psi_{0})\\}$

By going through the arguments in the proof of Lemma 2, we find that (14)
holds for $R_{1n}(\Psi^{(k)}(\bm{\beta}_{0}))$ as well, that is

$\displaystyle R_{1n}(\Psi^{(k)}(\bm{\beta}_{0}))\leq
2\bar{\bm{t}}^{T}\sum_{i=1}^{n}\bm{b}_{i}-n\bar{\bm{t}}^{T}\bm{B}\bar{\bm{t}}\\{1+o_{p}(1)\\}+o_{p}(1),$
(1)

where

$\displaystyle\bar{\bm{t}}^{T}=\left(\alpha_{1}^{(k)}-\alpha_{01},\ldots,\alpha_{m_{0}-1}^{(k)}-\alpha_{0m_{0}-1};\alpha_{1}^{(k)}\\{\bm{m}_{11}^{(k)}\\}^{T},\ldots,\alpha_{m_{0}}^{(k)}\\{\bm{m}_{1m_{0}}^{(k)}\\}^{T};\alpha_{1}^{(k)}\\{\bm{m}_{21}^{(k)}\\}^{T},\ldots,\alpha_{m_{0}}^{(k)}\\{\bm{m}_{2m_{0}}^{(k)}\\}^{T}\right),$
$\displaystyle\bm{b}_{i}=(\bm{b}_{1i}^{T},\bm{b}_{2i}^{T})^{T},~{}\bm{b}_{1i}=(\Delta_{i1},\ldots,\Delta_{im_{0}-1},\bm{Y}_{i}(\bm{\theta}_{01})^{T},\ldots,\bm{Y}_{i}(\bm{\theta}_{0m_{0}}))^{T},~{}\bm{b}_{2i}=(\bm{Z}_{i}(\bm{\theta}_{01})^{T},\ldots,\bm{Z}_{i}(\bm{\theta}_{0m_{0}})^{T})^{T}.$
(2)

where $\bm{m}_{1h}^{(k)}$ is a $p$-dimensional column vector

$\displaystyle\bm{m}_{1h}^{(k)}=\beta_{h}^{(k)}(\bm{\theta}_{1h}^{(k)}-\bm{\theta}_{0h})+(1-\beta_{h}^{(k)})(\bm{\theta}_{2h}^{(k)}-\bm{\theta}_{0h}),$

$\bm{m}_{2h}^{(k)}$ is a $p$-dimensional vector with elements

$\displaystyle\beta_{h}^{(k)}\left\\{\bm{\theta}_{1h}^{(k)}(u)-\bm{\theta}_{0h}(u)\right\\}^{2}+(1-\beta_{h}^{(k)})\left\\{\bm{\theta}_{2h}^{(k)}(u)-\bm{\theta}_{0h}(u)\right\\}^{2}$

for $u=1,\ldots,p$. Recall

$\displaystyle\Delta_{ih}=\frac{f(D_{i};\bm{\theta}_{0h})-f(D_{i};\bm{\theta}_{0m_{0}})}{f(D_{i};\Psi_{0})},$
$\displaystyle\bm{Y}_{i}(\bm{\theta})=\frac{f^{\prime}(D_{i};\bm{\theta})}{f(D_{i};\Psi_{0})}\in\mathbb{R}^{p\times
1}~{}\text{for any}~{}\bm{\theta}\in\mathbb{R}^{p\times 1}$
$\displaystyle\bm{Z}_{i}(\bm{\theta})=\frac{1}{f(D_{i};\Psi_{0})}\times\left(\frac{\partial^{2}f(D_{i};\bm{\theta})}{\partial\theta_{1}\partial\theta_{1}},\ldots,\frac{\partial^{2}f(D_{i};\bm{\theta})}{\partial\theta_{p}\partial\theta_{p}}\right)^{T}\in\mathbb{R}^{p\times
1}.$

Then due to (C8),

$\displaystyle\delta_{i}=\sum_{h=1}^{m_{0}-1}(\alpha_{h}^{(1)}-\alpha_{0h})\Delta_{ih}+\sum_{h=1}^{m_{0}}\left\\{\alpha_{h}^{(1)}\\{\bm{m}_{1h}^{(1)}\\}^{T}\bm{Y}_{i}(\bm{\theta}_{0h})+\alpha_{h}^{(1)}\\{\bm{m}_{2h}^{(1)}\\}^{T}\bm{Z}_{i}(\bm{\theta}_{0h})\right\\}+\sum_{i=1}^{m_{0}}\epsilon_{ih}.$

Here $\bar{\bm{t}}^{T}$ is a $(2m_{0}p+m_{0}-1)$-dimensional vector,
$\bm{b}_{1i}$ is $(m_{0}p+m_{0}-1)$-dimensional, $\bm{b}_{2i}$ is
$m_{0}p$-dimensional. Since $\sum_{i=1}^{n}\bm{b}_{i}=O_{p}(n^{1/2})$, we have

$\displaystyle R_{1n}(\Psi^{(k)}(\bm{\beta}_{0}))\leq
2\bar{\bm{t}}^{T}\sum_{i=1}^{n}\bm{b}_{i}-n\bar{\bm{t}}^{T}\bm{B}\bar{\bm{t}}\\{1+o_{p}(1)\\}+o_{p}(1),$
(3)

Let $\bm{t}=(\bm{t}_{1}^{T},\bm{t}_{2}^{T})^{T}$ with

$\displaystyle\bm{t}_{1}=\left(\alpha_{1}-\alpha_{01},\ldots,\alpha_{m_{0}-1}-\alpha_{0m_{0}-1};\alpha_{1}\\{\bm{m}_{11}\\}^{T},\ldots,\alpha_{m_{0}}\\{\bm{m}_{1m_{0}}\\}^{T}\right)^{T},$
$\displaystyle\bm{t}_{2}=\left(\alpha_{1}\bm{m}_{21}^{T},\ldots,\alpha_{m_{0}}\bm{m}_{2m_{0}}^{T}\right)^{T}.$

Then we replace $\bar{\bm{t}}$ in (3) by an arbitrary $\bm{t}$ and take the
supremum, the following holds,

$\displaystyle
R_{1n}(\Psi^{(k)}(\bm{\beta}_{0}))-R_{0n}\leq\sup_{\bm{t}}\left\\{2\bm{t}^{T}\sum_{i=1}^{n}\bm{b}_{i}-n\bm{t}^{T}\bm{B}\bm{t}\right\\}-(\sum_{i=1}^{n}\bm{b}_{1i})^{T}(n\bm{B}_{11})^{-1}(\sum_{i=1}^{n}\bm{b}_{1i})+o_{p}(1).$

Let $\tilde{\bm{t}}_{1}=\bm{t}_{1}+\bm{B}_{11}^{-1}\bm{B}_{12}\bm{t}_{2}$,
then because of the orthogonality,

$\displaystyle
2\bm{t}^{T}\sum_{i=1}^{n}\bm{b}_{i}-n\bm{t}^{T}\bm{B}\bm{t}=2\tilde{\bm{t}}_{1}^{T}\sum_{i=1}^{n}\bm{b}_{1i}-n\tilde{\bm{t}}_{1}^{T}\bm{B}_{11}\tilde{\bm{t}}_{1}+2\tilde{\bm{t}}_{2}^{T}\sum_{i=1}^{n}\bm{b}_{2i}-n\tilde{\bm{t}}_{2}^{T}\bm{B}_{22}\tilde{\bm{t}}_{2}.$

Then

$\displaystyle EM_{n}^{(k)}\leq$
$\displaystyle\sup_{\tilde{\bm{t}}_{1}}\left\\{2\tilde{\bm{t}}_{1}^{T}\sum_{i=1}^{n}\bm{b}_{1i}-n\tilde{\bm{t}}_{1}^{T}\bm{B}_{11}\tilde{\bm{t}}_{1}\right\\}+\sup_{\tilde{\bm{t}}_{2}}\left\\{2\tilde{\bm{t}}_{2}^{T}\sum_{i=1}^{n}\bm{b}_{2i}-n\tilde{\bm{t}}_{2}^{T}\bm{B}_{22}\tilde{\bm{t}}_{2}\right\\}$
$\displaystyle~{}~{}~{}-\left(\sum_{i=1}^{n}\bm{b}_{1i}\right)^{T}(n\bm{B}_{11})^{-1}(\sum_{i=1}^{n}\bm{b}_{1i})+o_{p}(1)$
$\displaystyle=\left(\sum_{i=1}^{n}\bm{b}_{1i}\right)^{T}(n\bm{B}_{11})^{-1}(\sum_{i=1}^{n}\bm{b}_{1i})+\sup_{\tilde{\bm{t}}_{2}}\left\\{2\tilde{\bm{t}}_{2}^{T}\sum_{i=1}^{n}\bm{b}_{2i}-n\tilde{\bm{t}}_{2}^{T}\bm{B}_{22}\tilde{\bm{t}}_{2}\right\\}$
$\displaystyle~{}~{}~{}-\left(\sum_{i=1}^{n}\bm{b}_{1i}\right)^{T}(n\bm{B}_{11})^{-1}(\sum_{i=1}^{n}\bm{b}_{1i})+o_{p}(1)$
$\displaystyle=\sup_{\tilde{\bm{t}}_{2}}\left\\{2\tilde{\bm{t}}_{2}^{T}\sum_{i=1}^{n}\bm{b}_{2i}-n\tilde{\bm{t}}_{2}^{T}\bm{B}_{22}\tilde{\bm{t}}_{2}\right\\}+o_{p}(1).$

Note that all the elements of $\bm{t}_{2}$ are non-negative, this proves the
first part.

(2) Next we show that the upper bound can be attained. Since EM algorithm
always increases the $pl_{n}$, it is good enough to consider $k=1$. In the
definition of $\bar{\bm{t}}$, there are a total of free
$\left\\{2m_{p}p+m_{0}-1\right\\}$ parameters. Hence we can find $\alpha_{h}$,
$\bm{\theta}_{1h},\bm{\theta}_{2h}$ such that the corresponding
$\bar{\bm{t}}_{2}$ attained the supremum in
$\sup_{\tilde{\bm{t}}_{2}}\left\\{2\tilde{\bm{t}}_{2}^{T}\sum_{i=1}^{n}\bm{b}_{2i}-n\tilde{\bm{t}}_{2}^{T}\bm{B}_{22}\tilde{\bm{t}}_{2}\right\\}$.
Let
$\bar{\bm{t}}_{1}=(n\bm{B}_{11})^{-1}\sum_{i=1}^{n}\bm{b}_{1i}-\bm{B}_{11}^{-1}\bm{B}_{12}$.
Then the resulting $EM_{n}^{(1)}$ attains the upper bound, hence concludes the
proof.

###  3.3 Proof of Theorem 3

Let $\hat{\bm{v}}=\text{argmax}_{\bm{v}\geq
0}\\{2\bm{v}^{T}\sum_{i=1}^{n}\tilde{\bm{b}}_{2i}-n\bm{v}^{T}\tilde{\bm{B}}_{22}\bm{v}\\}$.
By Lemma 6 of Li2010, which is first given by Kudo1963, we have
$(\bm{w}^{T}-\hat{\bm{v}}^{T}\tilde{\bm{B}}_{22})\hat{\bm{v}}=0$. Therefore

$\displaystyle\sup_{\bm{v}\geq
0}(2\bm{v}^{T}\bm{w}-\bm{v}^{T}\tilde{\bm{B}}_{22}\bm{v})=\hat{\bm{v}}^{T}\tilde{\bm{B}}_{22}\hat{\bm{v}}.$
(4)

Now we only need to find out the distribution of
$\hat{\bm{v}}^{T}\tilde{\bm{B}}_{22}\hat{\bm{v}}$. Note $\hat{\bm{v}}$ has
$m_{0}p$ elements. Depending on the positiveness of these elements, we could
form $2^{m_{0}p}$ regions for $\bm{w}$. More precisely, let
$\hat{\bm{v}}=(\hat{\bm{v}}_{1},\hat{\bm{v}}_{2})^{T}$ be a partition of
$\hat{\bm{v}}$, in which $\hat{\bm{v}}_{1}$ is a vector of zeros of dimension
$\\{m_{0}p-s\\}$, and $\hat{\bm{v}}_{2}$ is a vector of rest positive elements
of dimension $s$. We can partition $\tilde{\bm{B}}_{22}$ correspondingly as

$\displaystyle\tilde{\bm{B}}_{22}=\left(\begin{array}[]{cc}\bm{W}_{11}&\bm{W}_{12}\\\
\bm{W}_{21}&\bm{W}_{22}\\\ \end{array}\right).$

By Lemma 6 of Li2010, $\hat{\bm{v}}_{1}=0$ and $\hat{\bm{v}}_{2}>0$ is
equivalent with

$\displaystyle\bm{W}_{22}^{-1}\bm{w}_{2}>0,~{}~{}~{}\bm{W}_{12}\bm{W}_{22}^{-1}\bm{w}_{2}-\bm{w}_{1}>0.$

Therefore the partition $\hat{\bm{v}}=(\hat{\bm{v}}_{1},\hat{\bm{v}}_{2})^{T}$
is uniquely related with the region of $\bm{w}$, which is defined by

$\displaystyle\mathcal{U}=\\{\bm{w}:\bm{W}_{22}^{-1}\bm{w}_{2}>0,~{}\bm{W}_{12}\bm{W}_{22}^{-1}\bm{w}_{2}-\bm{w}_{1}>0\\}.$

Denote the collection of all possible $\mathcal{U}$ by $\mathcal{S}$, then for
any $x$,

$\displaystyle\text{P}(\hat{\bm{v}}^{T}\tilde{\bm{B}}_{22}\hat{\bm{v}}\leq
x)=\sum_{\mathcal{U}\in\mathcal{S}}\text{P}(\hat{\bm{v}}^{T}\tilde{\bm{B}}_{22}\hat{\bm{v}}\leq
x|\bm{w}\in\mathcal{U})\text{P}(\bm{w}\in\mathcal{U}).$

Because
$\hat{\bm{v}}^{T}\tilde{\bm{B}}_{22}\hat{\bm{v}}=\bm{w}_{2}^{T}\bm{W}_{22}^{-1}\bm{w}_{2}$

$\displaystyle\text{P}(\hat{\bm{v}}^{T}\tilde{\bm{B}}_{22}\hat{\bm{v}}\leq
x|\bm{w}\in\mathcal{U})$
$\displaystyle=\text{P}(\bm{w}_{2}^{T}\bm{W}_{22}^{-1}\bm{w}_{2}\leq
x|\bm{W}_{22}^{-1}\bm{w}_{2}>0，~{}\bm{W}_{12}\bm{W}_{22}^{-1}\bm{w}_{2}-\bm{w}_{1}>0)$
$\displaystyle=\text{P}(\bm{w}_{2}^{T}\bm{W}_{22}^{-1}\bm{w}_{2}\leq
x|\bm{W}_{22}^{-1}\bm{w}_{2}>0)$

because of the normality of $\bm{w}$, and the fact that $\bm{w}_{2}$ and
$\bm{W}_{12}\bm{W}_{22}^{-1}\bm{w}_{2}-\bm{w}_{1}$ are independent. Moreover,
$\bm{w}_{2}^{T}\bm{W}_{22}^{-1}\bm{w}_{2}$ and
$\bm{W}_{22}^{-1/2}\bm{w}_{2}/\sqrt{\bm{w}_{2}^{T}\bm{W}_{22}^{-1}\bm{w}_{2}}$
are independent. Hence

$\displaystyle\text{P}(\hat{\bm{v}}^{T}\tilde{\bm{B}}_{22}\hat{\bm{v}}\leq
x|\bm{w}\in\mathcal{U})$
$\displaystyle=\text{P}(\bm{w}_{2}^{T}\bm{W}_{22}^{-1}\bm{w}_{2}\leq
x|\bm{W}_{22}^{-1/2}\bm{w}_{2}/\sqrt{\bm{w}_{2}^{T}\bm{W}_{22}^{-1}\bm{w}_{2}}>0)$
$\displaystyle=\text{P}(\bm{w}_{2}^{T}\bm{W}_{22}^{-1}\bm{w}_{2}\leq x).$

In other words, $\hat{\bm{v}}^{T}\tilde{\bm{B}}_{22}\hat{\bm{v}}$ follows
$\chi_{s}^{2}$ distribution where $s$ is the number of positive elements in
$\hat{\bm{v}}$. Let $\mathcal{U}_{s}$ be the collection of $\bm{w}$ whose
corresponding $\hat{\bm{v}}$ has $s$ positive components. Define
$a_{s}=P(\mathcal{U}_{s})$ and $\bar{s}=m_{0}p$, then
$\hat{\bm{v}}^{T}\tilde{\bm{B}}_{22}\hat{\bm{v}}\sim\sum_{s=0}^{\bar{s}}a_{s}\chi_{s}^{2}$,
the conclusion follows.

##  4 Proof of Lemmas 1–3

###  4.1 Proof of Lemma 1

Note that $\hat{\Psi}_{0}\in\Omega_{2m_{0}}(\bm{\beta}_{0})$. Hence
$\Psi^{(1)}\in\Omega_{2m_{0}}(\bm{\beta}_{0})$ and it maximizes $pl_{n}(\Psi)$
within the class. Therefore $pl_{n}(\Psi^{(1)})\geq pl_{n}(\hat{\Psi}_{0})$.
Hence

$l_{n}(\Psi^{(1)})+p(\bm{\beta}_{0})=pl_{n}(\Psi^{(1)})\geq
pl_{n}(\hat{\Psi}_{0})\geq l_{n}(\Psi_{0})+p(\bm{\beta}_{0}),$

which gives $l_{n}(\Psi^{(1)})\geq l_{n}(\Psi_{0})$. Following the same
arguments of the consistency of MLE in the proof of Lemma 1 in Li2010, we
conclude that $\|\Psi^{(1)}(\bm{\beta}_{0})-\Psi_{0}\|\rightarrow 0$ almost
surely with $\|F(\bm{\theta})\|=\sup_{\bm{\theta}\in\Theta}|F(\bm{\theta})|$.
Then we have the consistency of $\alpha^{(1)}$, $\bm{\beta}^{(1)}$. Given the
consistency of $\Psi$, We also have the consistency of
$f(D;\bm{\theta}^{(1)},\bm{\gamma}^{(1)})$, hence the consistency of
$\bm{\theta}_{1}^{(1)}$, $\bm{\theta}_{2}^{(1)}$ and $\bm{\gamma}^{(1)}$ holds
as well.

###  4.2 Proof of Lemma 2

Let
$\\{l_{n}(\Psi^{(1)}(\bm{\beta}_{0}))-l_{n}(\Psi_{0})\\}=\sum_{i=1}^{n}\log(1+\delta_{i})$,
where

$\displaystyle\delta_{i}$
$\displaystyle=\frac{f(D_{i};\Psi^{(1)}(\bm{\beta}_{0}))-f(D_{i};\Psi_{0})}{f(D_{i};\Psi_{0})}$
$\displaystyle=\sum_{h=1}^{m_{0}}(\alpha_{h}^{(1)}-\alpha_{0h})\frac{f(Di;\bm{\theta}_{0h})}{f(D_{i};\Psi_{0})}$
$\displaystyle~{}~{}~{}~{}~{}~{}+\sum_{h=1}^{m_{0}}\alpha_{h}^{(1)}\left\\{\beta_{h}^{(1)}\frac{f(D_{i};\bm{\theta}_{1h}^{(1)})-f(D_{i};\bm{\theta}_{0h})}{f(D_{i};\Psi_{0})}+(1-\beta_{h}^{(1)})\frac{f(D_{i};\bm{\theta}_{2h}^{(1)})-f(D_{i};\bm{\theta}_{0h})}{f(D_{i};\Psi_{0})}\right\\}$
$\displaystyle=\sum_{h=1}^{m_{0}-1}(\alpha_{h}^{(1)}-\alpha_{0h})\frac{f(D_{i};\bm{\theta}_{0h})-f(D_{i};\bm{\theta}_{0m_{0}})}{f(D_{i};\Psi_{0})}$
$\displaystyle~{}~{}~{}~{}~{}~{}+\sum_{h=1}^{m_{0}}\alpha_{h}^{(1)}\left\\{\beta_{h}^{(1)}\frac{f(D_{i};\bm{\theta}_{1h}^{(1)})-f(D_{i};\bm{\theta}_{0h})}{f(D_{i};\Psi_{0})}+(1-\beta_{h}^{(1)})\frac{f(D_{i};\bm{\theta}_{2h}^{(1)})-f(D_{i};\bm{\theta}_{0h})}{f(D_{i};\Psi_{0})}\right\\}.$
(5)

Note that

$\displaystyle f(D_{i};\bm{\theta}_{1h}^{(1)})-f(D_{i};\bm{\theta}_{0h})$
$\displaystyle=(\bm{\theta}_{1h}^{(1)}-\bm{\theta}_{0h})^{T}f^{\prime}(D_{i};\bm{\theta}_{0h})+\frac{1}{2}(\bm{\theta}_{1h}^{(1)}-\bm{\theta}_{0h})^{T}f^{\prime\prime}(D_{i};\bm{\theta}_{0h})(\bm{\theta}_{1h}^{(1)}-\bm{\theta}_{0h})+\epsilon_{ih}.$

Define

$\displaystyle\bm{m}_{1h}^{(k)}=\beta_{h}^{(k)}(\bm{\theta}_{1h}^{(k)}-\bm{\theta}_{0h})+(1-\beta_{h}^{(k)})(\bm{\theta}_{2h}^{(k)}-\bm{\theta}_{0h}).$
(6)

For index $u=1,\ldots,p$, let $\bm{m}_{2h}^{(k)}$ be a $p$-dimensional vector
with elements

$\displaystyle{\color[rgb]{0,0,1}\beta_{h}^{(k)}\left\\{\bm{\theta}_{1h}^{(k)}(u)-\bm{\theta}_{0h}(u)\right\\}^{2}+(1-\beta_{h}^{(k)})\left\\{\bm{\theta}_{2h}^{(k)}(u)-\bm{\theta}_{0h}(u)\right\\}^{2}.}$

Let $\bm{m}_{1}^{(k)}=(\bm{m}_{11}^{(k)},\ldots,\bm{m}_{1m_{0}}^{(k)})$ and
$\bm{m}_{2}^{(k)}=(\bm{m}_{21}^{(k)},\ldots,\bm{m}_{2m_{0}}^{(k)})$ be two
matrices. Define

$\displaystyle\Delta_{ih}=\frac{f(D_{i};\bm{\theta}_{0h})-f(D_{i};\bm{\theta}_{0m_{0}})}{f(D_{i};\Psi_{0})},$
$\displaystyle\bm{Y}_{i}(\bm{\theta})=\frac{f^{\prime}(D_{i};\bm{\theta})}{f(D_{i};\Psi_{0})}\in\mathbb{R}^{p\times
1}~{}\text{for any}~{}\bm{\theta}\in\mathbb{R}^{p\times 1}$
$\displaystyle\bm{Z}_{i}(\bm{\theta})=\frac{1}{f(D_{i};\Psi_{0})}\times{\color[rgb]{0,0,1}\left(\frac{\partial^{2}f(D_{i};\bm{\theta})}{\partial\theta_{1}\partial\theta_{1}},\ldots,\frac{\partial^{2}f(D_{i};\bm{\theta})}{\partial\theta_{p}\partial\theta_{p}}\right)^{T}\in\mathbb{R}^{p\times
1},}$
$\displaystyle\bm{b}_{1i}=(\Delta_{i1},\ldots,\Delta_{im_{0}-1},\bm{Y}_{i}(\bm{\theta})^{T},\ldots,\bm{Y}_{i}(\bm{\theta})^{T})^{T}\in\mathbb{R}^{(m_{0}-1+m_{0}p)\times
1},$
$\displaystyle\bm{b}_{2i}=(\bm{Z}_{i}(\bm{\theta}_{01})^{T},\ldots,\bm{Z}_{i}(\bm{\theta}_{0m_{0}})^{T})^{T}\in\mathbb{R}^{m_{0}p\times
1},$ (7)

and $\bm{b}_{i}=(\bm{b}_{1i}^{T},\bm{b}_{2i}^{T})^{T}$. Then

$\displaystyle\delta_{i}=\sum_{h=1}^{m_{0}-1}(\alpha_{h}^{(1)}-\alpha_{0h})\Delta_{ih}+\sum_{h=1}^{m_{0}}\left\\{\alpha_{h}^{(1)}\\{\bm{m}_{1h}^{(1)}\\}^{T}\bm{Y}_{i}(\bm{\theta}_{0h})+\alpha_{h}^{(1)}\\{\bm{m}_{2h}^{(1)}\\}^{T}\bm{Z}_{i}(\bm{\theta}_{0h})\right\\}+\sum_{i=1}^{m_{0}}\epsilon_{ih}.$

For $i=1,\ldots,n$ and index $u$, we define a $p$-dimensional vector
$\bm{U}_{ih}(\bm{\theta})$ as the collection of elements

$\displaystyle\frac{f(D_{i};\bm{\theta})-f(D_{i};\bm{\theta}_{0h})-f^{\prime}(D_{i};\bm{\theta}_{0h})(\bm{\theta}-\bm{\theta}_{0h})-(\bm{\theta}-\bm{\theta}_{0h})^{T}f^{\prime\prime}(D_{i};\bm{\theta}_{0h})(\bm{\theta}-\bm{\theta}_{0h})/2}{f(D_{i};\Psi_{0})\left\\{\bm{\theta}(u)-\bm{\theta}_{0h}(u)\right\\}^{3}},$

where $\bm{\theta}(u)$ is the $u$-th element of $\bm{\theta}$. Then Condition
(C4) implies that for every $h=1,\ldots,m_{0}$, the process
$n^{-1/2}\sum_{i=1}^{n}\bm{U}_{ih}(\bm{\theta})$ is tight for any
$\bm{\theta}\in B(\bm{\theta}_{0h},\epsilon_{0})$ (Billing1968). Then

$\displaystyle
n^{-1/2}\sum_{i=1}^{n}\bm{U}_{ih}(\tilde{\bm{\theta}})=O_{p}(1)$ (8)

when $\tilde{\bm{\theta}}$ is a consistent estimator of $\bm{\theta}$.

Again, for index $u=1,\ldots,p$, define $\bm{\theta}^{3}$ as a $p$-dimensional
vector with elements of $\bm{\theta}(u)^{3}$. Then

$\displaystyle\epsilon_{ih}=\alpha_{h}^{(1)}\left\\{\beta_{h}^{(1)}\\{(\bm{\theta}_{1h}^{(1)}-\bm{\theta}_{0h})^{3}\\}^{T}\bm{U}_{ih}(\bm{\theta}_{1h}^{(1)})+(1-\beta_{h}^{(1)})\\{(\bm{\theta}_{2h}^{(1)}-\bm{\theta}_{0h})^{3}\\}^{T}\bm{U}_{ih}(\bm{\theta}_{2h}^{(1)})\right\\}.$
(9)

Then

$\displaystyle\sum_{i=1}^{n}\delta_{i}=\sum_{i=1}^{n}\left\\{\sum_{h=1}^{m_{0}-1}(\alpha_{h}^{(1)}-\alpha_{0h})\Delta_{ih}+\sum_{h=1}^{m_{0}}\left\\{\alpha_{h}^{(1)}\\{\bm{m}_{1h}^{(1)}\\}^{T}\bm{Y}_{i}(\bm{\theta}_{0h})+\alpha_{h}^{(1)}\\{\bm{m}_{2h}^{(1)}\\}^{T}\bm{Z}_{i}(\bm{\theta}_{0h})\right\\}\right\\}+\epsilon_{n},$

where $|\epsilon_{n}|$ equals to

$\displaystyle\left|\sum_{i=1}^{n}\sum_{h=1}^{m_{0}}\epsilon_{ih}\right|$
$\displaystyle\leq
n^{1/2}\sum_{h=1}^{m_{0}}\alpha_{h}^{(1)}\left\\{\beta_{h}^{(1)}\\{|\bm{\theta}_{1h}^{(1)}-\bm{\theta}_{0h}|^{3}\\}^{T}n^{-1/2}\sum_{i=1}^{n}\bm{U}_{ih}(\bm{\theta}_{1h}^{(1)})+(1-\beta_{h}^{(1)})\\{|\bm{\theta}_{2h}^{(1)}-\bm{\theta}_{0h}|^{3}\\}^{T}n^{-1/2}\sum_{i=1}^{n}\bm{U}_{ih}(\bm{\theta}_{2h}^{(1)})\right\\}$
$\displaystyle=O_{p}(n^{1/2})\sum_{h=1}^{m_{0}}\alpha_{h}^{(1)}\left\\{\beta_{h}^{(1)}\\{|\bm{\theta}_{1h}^{(1)}-\bm{\theta}_{0h}|^{3}\\}^{T}+(1-\beta_{h}^{(1)})\\{|\bm{\theta}_{2h}^{(1)}-\bm{\theta}_{0h}|^{3}\\}^{T}\right\\}\bm{1}$
$\displaystyle=o_{p}(n^{1/2})\sum_{h=1}^{m_{0}}\\{\bm{m}_{2h}^{(1)}\\}^{T}\bm{1}$
$\displaystyle\leq
o_{p}(1)+o_{p}(n)\sum_{h=1}^{m_{0}}\|\bm{m}_{2h}^{(1)}\|_{\infty}^{2}.$ (10)

Let

$\displaystyle\bar{\bm{t}}=\left(\alpha_{1}^{(1)}-\alpha_{01},\ldots,\alpha_{m_{0}-1}^{(1)}-\alpha_{0m_{0}-1};\alpha_{1}^{(1)}\\{\bm{m}_{11}^{(1)}\\}^{T},\ldots,\alpha_{m_{0}}^{(1)}\\{\bm{m}_{1m_{0}}^{(1)}\\}^{T};\alpha_{1}^{(1)}\\{\bm{m}_{21}^{(1)}\\}^{T},\ldots,\alpha_{m_{0}}^{(1)}\\{\bm{m}_{2m_{0}}^{(1)}\\}^{T}\right)^{T}$

and $\bm{b}_{i}=(\bm{b}_{1i}^{T},\bm{b}_{2i}^{T})^{T}$. Then

$\displaystyle\sum_{i=1}^{n}\delta_{i}=\sum_{i=1}^{n}\bar{\bm{t}}^{T}\bm{b}_{i}+\epsilon_{n}.$
(11)

Using the fact that $\log(1+x)\leq x-x^{2}/2+x^{3}/3$, we have

$\displaystyle 2\sum_{i=1}^{n}\log(1+\delta_{i})$ $\displaystyle\leq
2\sum_{i=1}^{n}\delta_{i}-\sum_{i=1}^{n}\delta_{i}^{2}+\frac{2}{3}\sum_{i=1}^{n}\delta_{i}^{3}$
$\displaystyle=2\sum_{i=1}^{n}\bar{\bm{t}}^{T}\bm{b}_{i}-\sum_{i=1}^{n}(\bar{\bm{t}}^{T}\bm{b}_{i})^{2}+\frac{2}{3}\sum_{i=1}^{n}(\bar{\bm{t}}^{T}\bm{b}_{i})^{3}+O_{p}(\epsilon_{n}).$
(12)

By Condition (C5),

$\displaystyle\sum_{i=1}^{n}(\bar{\bm{t}}^{T}\bm{b}_{i})^{2}=n\bar{\bm{t}}^{T}\bm{B}\bar{\bm{t}}\\{1+o_{p}(1)\\},$
$\displaystyle\sum_{i=1}^{n}(\bar{\bm{t}}^{T}\bm{b}_{i})^{3}=o_{p}(n)\bar{\bm{t}}^{T}\bar{\bm{t}}.$
(13)

Note that $\epsilon_{n}=o_{p}(1)+o_{p}(n)\bar{\bm{t}}^{T}\bar{\bm{t}}$, hence

$\displaystyle 2\sum_{i=1}^{n}\log(1+\delta_{i})\leq
2\bar{\bm{t}}^{T}\sum_{i=1}^{n}\bm{b}_{i}-n\bar{\bm{t}}^{T}\bm{B}\bar{\bm{t}}\\{1+o_{p}(1)\\}+o_{p}(1).$
(14)

Because $pl_{n}(\Psi^{(1)})\geq pl_{n}(\Psi_{0})$, we have

$\displaystyle
2\\{l_{n}(\Psi^{(1)})-l_{n}(\Psi_{0})\\}+2n\left\\{q_{\lambda}(\Psi_{0})-q_{\lambda}(\Psi^{(1)})\right\\}\geq
2p(\bm{\beta_{0}}).$ (15)

Observe that

$\displaystyle|q_{\lambda}(\Psi_{0})-q_{\lambda}(\Psi^{(1)}(\bm{\beta}_{0}))|\leq\lambda\times$
$\displaystyle\sum_{k=1}^{p}\frac{|\sum_{j=2}^{m_{0}}(\theta_{0jk}-\theta_{01k})^{2}+\sum_{j^{\prime}=1}^{m_{0}}(\theta_{0j^{\prime}k}-\theta_{01k})^{2}-\sum_{j=2}^{m_{0}}(\theta_{1jk}-\theta_{11k})^{2}-\sum_{j^{\prime}=1}^{m_{0}}(\theta_{2j^{\prime}k}-\theta_{11k})^{2}|}{\\{\sum_{j=2}^{m_{0}}(\theta_{0jk}-\theta_{01k})^{2}+\sum_{j^{\prime}=1}^{m_{0}}(\theta_{0j^{\prime}k}-\theta_{01k})^{2}\\}^{1/2}+\\{\sum_{j=2}^{m_{0}}(\theta_{1jk}-\theta_{11k})^{2}+\sum_{j^{\prime}=1}^{m_{0}}(\theta_{2j^{\prime}k}-\theta_{11k})^{2}\\}^{1/2}}$
$\displaystyle\leq
C_{1}\lambda\left|\sum_{j=2}^{m_{0}}(\theta_{0jk}-\theta_{01k})^{2}+\sum_{j^{\prime}=1}^{m_{0}}(\theta_{0j^{\prime}k}-\theta_{01k})^{2}-\sum_{j=2}^{m_{0}}(\theta_{1jk}-\theta_{11k})^{2}-\sum_{j^{\prime}=1}^{m_{0}}(\theta_{2j^{\prime}k}-\theta_{11k})^{2}\right|$
$\displaystyle\leq C_{2}|\bar{\bm{t}}^{T}\bm{1}|$ (16)

for some positive constants $C_{1},C_{2}$ because of the Condition (C8). Note
that $n\lambda=o(n^{1/2})$ due to Condition (C7), (15) and (14) together
implies

$\displaystyle 2p(\bm{\beta}_{0})\leq
2\bar{\bm{t}}^{T}\left\\{\sum_{i=1}^{n}\bm{b}_{i}+o_{p}(n^{1/2})\bm{1}\right\\}-n\bar{\bm{t}}^{T}\bm{B}\bar{\bm{t}}\\{1+o_{p}(1)\\}+o_{p}(1)$
(17)

Because $n^{-1/2}\sum_{i=1}^{n}\bm{b}_{i}=O_{p}(1)$ and $p(\bm{\beta}_{0})$ is
fixed, we have $\bar{\bm{t}}=O_{p}(n^{-1/2})$. This implies
$\alpha_{1}^{(s)}-\alpha_{0s}=O_{p}(n^{-1/2})$ for $s=1,\ldots,m_{0}-1$,
$\bm{m}_{1t}^{(1)}=O_{p}(n^{-1/2})$ and $\bm{m}_{2t}^{(1)}=O_{p}(n^{-1/2})$
for every $t=1,\ldots,m_{0}$. The $n^{-1/2}$ convergence rate of
$\bm{m}_{2t}^{(1)}$ implies the $n^{-1/2}$ convergence rate of
$\bm{\theta}_{1h}^{(k)}$ and $\bm{\theta}_{2h}^{(k)}$. For the rate of
$\bm{\beta}^{(1)}$, it holds trivially because
$\bm{\beta}^{(1)}=\bm{\beta}_{0}$ by definition. Once we obtain the
convergence rate of the parameters in $\Psi^{(1)}$, the convergence rate of
$\bm{\gamma}^{(1)}$ is obvious. This completes the proof.

###  4.3 Proof of Lemma 3

We divide the proof into three parts.

(1) First, we show $\bm{\beta}^{(k+1)}-\bm{\beta}_{0}=O_{p}(n^{-1/6})$. Here
we first proof a useful result,

$\displaystyle\sum_{i=1}^{n}\frac{f(D_{i};\bm{\theta}_{1h}^{(k)})}{f(D_{i};\Psi^{(k)}(\bm{\beta}_{0}))}=n\\{1+O_{p}(n^{-1/6})\\}.$
(18)

Note that for every $j=1,2$ and $h=1,\ldots,p$,

$\displaystyle\frac{f(D_{i};\bm{\theta}_{jh}^{(k)})-f(D_{i};\bm{\theta}_{0h})}{f(D_{i};\Psi_{0})}=\bm{Y}_{i}(\bm{\theta}_{0h})^{T}(\bm{\theta}_{jh}^{(k)}-\bm{\theta}_{0h})+(\bm{\theta}_{jh}^{(k)}-\bm{\theta}_{0h})^{T}\tilde{\bm{Z}}_{i}(\tilde{\bm{\theta}}_{jh}^{(k)})(\bm{\theta}_{jh}^{(k)}-\bm{\theta}_{0h}),$
(19)

where $\tilde{\bm{Z}}_{i}(\tilde{\bm{\theta}}_{jh}^{(k)})$ is a $p$ by $p$
diagnoal matrix defined by

$\displaystyle\tilde{\bm{Z}}_{i}(\bm{\theta})=\frac{1}{f(D_{i};\Psi_{0})}\text{diag}\left\\{\frac{\partial^{2}f(D_{i};\bm{\theta})}{\partial\theta_{1}^{2}},\ldots,\frac{\partial^{2}f(D_{i};\bm{\theta})}{\partial\theta_{p}^{2}}\right\\},$
(20)

and $\tilde{\bm{\theta}}_{jh}^{(k)}$ is in the neighborhood of
$\bm{\theta}_{0h}$. Then

$\displaystyle\delta_{i}$
$\displaystyle=\frac{f(D_{i};\Psi^{(k)}(\bm{\beta}_{0}))-f(D_{i};\Psi_{0})}{f(D_{i};\Psi_{0})}$
$\displaystyle=\sum_{h=1}^{m_{0}-1}(\alpha_{h}^{(k)}-\alpha_{0})\Delta_{ih}+\sum_{h=1}^{m_{0}}\alpha_{h}^{(k)}\\{\bm{m}_{1h}^{(k)}\\}^{T}\bm{Y}_{i}(\bm{\theta}_{0h})+\sum_{h=1}^{m_{0}}\alpha_{h}^{(k)}\beta_{h}^{(k)}(\bm{\theta}_{1h}^{(k)}-\bm{\theta}_{0h})^{T}\tilde{\bm{Z}}_{i}(\tilde{\bm{\theta}}_{1h}^{(k)})(\bm{\theta}_{1h}^{(k)}-\bm{\theta}_{0h})$
$\displaystyle~{}~{}+\sum_{h=1}^{m_{0}}\alpha_{h}^{(k)}\beta_{h}^{(k)}(\bm{\theta}_{2h}^{(k)}-\bm{\theta}_{0h})^{T}\tilde{\bm{Z}}_{i}(\tilde{\bm{\theta}}_{2h}^{(k)})(\bm{\theta}_{2h}^{(k)}-\bm{\theta}_{0h}).$

Because $\Delta_{ih}$ has a constant upper bound, $\max_{1\leq i\leq
n}\sum_{h=1}^{m_{0}-1}|\alpha_{h}^{(k)}-\alpha_{0})\Delta_{ih}|=O_{p}(n^{-1/4})$.
For the second term of $\delta_{i}$, because of Condition (C4), $\max_{1\leq
i\leq n}\|\bm{Y}_{i}(\bm{\theta}_{0h})\|_{\infty}=O_{p}(n^{1/3})$. From the
proof of Lemma 2, $\bm{m}_{1h}^{(k)}=O_{p}(n^{-1/2})$. Therefore the second
term of $\delta_{i}$ is $O_{p}(n^{-1/6})$ uniformly. For the third and fourth
terms, Condition (C4) implies $\max_{1\leq i\leq
n}\|\tilde{\bm{Z}}_{i}(\tilde{\bm{\theta}}_{2h}^{(k)})\|_{2}=O_{p}(n^{1/3})$.
Combine with the assumption that
$\bm{\theta}_{j}^{(k)}-\bm{\theta}_{0}=O_{p}(n^{-1/4})$ for $j=1,2$, both
third and fourth terms of $\delta_{i}$ are $O_{p}(n^{-1/6})$. Therefore
$\max_{1\leq i\leq n}|\delta_{i}|=O_{p}(n^{-1/6})$. In other words, uniformly
in $i$,

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\frac{f(D_{i};\bm{\theta}_{1h}^{(k)})}{f(D_{i};\Psi^{(k)}(\bm{\beta}_{0})}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{f(D_{i};\bm{\theta}_{1h}^{(k)})}{f(D_{i};\Psi_{0})}\frac{f(D_{i};\Psi_{0})}{f(D_{i};\Psi^{(k)}(\bm{\beta}_{0})}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{f(D_{i};\bm{\theta}_{1h}^{(k)})}{f(D_{i};\Psi_{0})}\\{1+O_{p}(n^{-1/6})\\}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left\\{\frac{f(D_{i};\bm{\theta}_{0h})}{f(D_{i};\Psi_{0})}+\bm{Y}_{i}(\bm{\theta}_{0h})^{T}(\bm{\theta}_{jh}^{(k)}-\bm{\theta}_{0h})+(\bm{\theta}_{jh}^{(k)}-\bm{\theta}_{0h})^{T}\tilde{\bm{Z}}_{i}(\tilde{\bm{\theta}}_{jh}^{(k)})(\bm{\theta}_{jh}^{(k)}-\bm{\theta}_{0h})\right\\}\\{1+O_{p}(n^{-1/6})\\}$
$\displaystyle=1+O_{p}(n^{-1/6}).$

This shows (18). By the definition, $\beta_{h}^{(k+1)}$ maximizes

$\displaystyle
Q_{nh}(\beta)=\sum_{i=1}^{n}w_{i1h}^{(k)}\log(\beta)+\sum_{i=1}^{n}w_{i2h}^{(k)}\log(1-\beta)+p(\beta)$

with

$\displaystyle\sum_{i=1}^{n}w_{i1h}^{(k)}=\sum_{i=1}^{n}\alpha_{h}^{(k)}\beta_{h}^{(k)}\frac{f(D_{i};\bm{\theta}_{1h}^{(k)})}{f(D_{i};\Psi^{(k)}(\bm{\beta}_{0})}=n\alpha_{h}^{(k)}\beta_{h}^{(k)}\\{1+O_{p}(n^{-1/6})\\},$
$\displaystyle\sum_{i=1}^{n}w_{i2h}^{(k)}=n\alpha_{h}^{(k)}-\sum_{i=1}^{n}w_{i1h}^{(k)}=n\alpha_{h}^{(k)}(1-\beta_{h}^{(k)})\\{1+O_{p}(n^{-1/6})\\}$
(21)

Let

$\displaystyle
H_{nh}(\beta)=\sum_{i=1}^{n}w_{i1h}^{(k)}\log(\beta)+\sum_{i=1}^{n}w_{i2h}^{(k)}\log(1-\beta).$
(22)

Its maximizer is

$\displaystyle\hat{\beta}_{h}=\beta_{h}^{(k)}\\{1+O_{p}(n^{-1/6})\\}.$ (23)

Then consider any $\beta^{*}$ such that
$n^{-1/6}\leq|\beta^{*}-\hat{\beta}_{h}|\leq n^{-1/3}$, by Taylor expansion,

$\displaystyle
H_{nh}(\hat{\beta}_{h})-H_{nh}(\beta^{*})\geq\epsilon\alpha_{h}^{(k)}n(\beta^{*}-\hat{\beta}_{h})^{2}\geq\epsilon\alpha_{h}^{(k)}n^{2/3}$

for some $\epsilon>0$. Then

$\displaystyle Q_{nh}(\beta^{*})-Q_{nh}(\hat{\beta}_{h})\leq
p(\beta^{*})-p(\hat{\beta}_{h})-\epsilon\alpha_{h}^{(k)}n^{2/3}<0$

when $n$ is large enough. Therefore the maximizer of $Q_{nh}$ must be within
an $O_{p}(n^{-1/6})$-neighborhood of $\hat{\beta}_{h}$.

(2) Next, we show the consistency of $\bm{\alpha}^{(k+1)}$,
$\bm{\theta}_{1}^{(k+1)}$,$\bm{\theta}_{2}^{(k+1)}$ and $\bm{\gamma}^{(k+1)}$,
i.e.,

$\displaystyle\bm{\alpha}^{(k+1)}-\bm{\alpha}_{0}=o_{p}(1),~{}\bm{\theta}_{1}^{(k+1)}-\bm{\theta}_{0}=o_{p}(1),~{}\bm{\theta}_{2}^{(k+1)}-\bm{\theta}_{0}=o_{p}(1),~{}\bm{\gamma}^{(k+1)}-\bm{\gamma}_{0}=o_{p}(1).$

For $l=1,\ldots,m_{0}$, define a new mixture distribution such that its $l$-th
support points have been updated by the EM step while the rest points keep the
same with $\Psi^{(k)})$, i.e.,

$\displaystyle\Psi_{l}^{(k+1)}(\theta)$ $\displaystyle=\sum_{1\leq h\leq
m_{0};h\neq
l}\alpha_{h}^{(k)}\left\\{\beta_{h}^{(k)}I(\text{P}\bm{\theta}_{1h}^{(k)}\leq\theta)+(1-\beta_{h}^{(k)})I(\text{P}\bm{\theta}_{2h}^{(k)}\leq\theta)\right\\}$
$\displaystyle~{}~{}+\alpha_{l}^{(k)}\left\\{\beta_{l}^{(k)}I(\text{P}\bm{\theta}_{1l}^{(k+1)}\leq\theta)+(1-\beta_{l}^{(k)})I(\text{P}\bm{\theta}_{2l}^{(k+1)}\leq\theta)\right\\}.$

Because EM increases the likelihood (Wu1981), we have

$\displaystyle l_{n}(\Psi_{l}^{(k+1)})\geq pl_{n}(\Psi_{l}^{(k+1)})\geq
pl_{n}(\Psi_{l}^{(k)})\geq
pl_{n}(\hat{\Psi}_{0})=l_{n}(\hat{\Psi}_{0})+p(\bm{\beta}_{0})\geq
l_{n}(\Psi_{0})+p(\bm{\beta}_{0}).$

This implies $l_{n}(\Psi_{l}^{(k+1)})\geq l_{n}(\Psi_{0})$. Similarly with the
proof of Lemma 1, we have the consistency of $\Psi_{l}^{(k+1)}$ for
$\Psi_{0}$. This implies the consistency of $\bm{\theta}_{1h}^{(k+1)}$ and
$\bm{\theta}_{2h}^{(k+1)}$ are consistent for $\bm{\theta}_{0l}$. This holds
for every $l=1,\ldots,m_{0}$. Hence $\bm{\theta}_{1}$ and $\bm{\theta}_{2}$
are consistent. Similar idea can be applied on
$\alpha^{(k+1)}$,$\bm{\gamma}^{(k+1)}$ as well to obtain their consistency.

(3) For the last part, we derive the convergence rate of
$\bm{\alpha}^{(k+1)}$, $\bm{\theta}_{1}^{(k+1)}$, $\bm{\theta}_{2}^{(k+1)}$
and $\bm{m}_{1}^{(k+1)}$ based on consistency result in part (2). This
proceeds in the same way with the proof of Lemma 2. Let
$\\{l_{n}(\Psi^{(k+1)}(\bm{\beta}_{0}))-l_{n}(\Psi_{0})\\}=\sum_{i=1}^{n}\log(1+\delta_{i})$,
where

$\displaystyle\delta_{i}$
$\displaystyle=\frac{f(D_{i};\Psi^{(k+1)}(\bm{\beta}_{0}))-f(D_{i};\Psi_{0})}{f(D_{i};\Psi_{0})}$
$\displaystyle=\sum_{h=1}^{m_{0}-1}(\alpha_{h}^{(k+1)}-\alpha_{0h})\Delta_{ih}+\sum_{h=1}^{m_{0}}\left\\{\alpha_{h}^{(k+1)}\\{\bm{m}_{1h}^{(k+1)}\\}^{T}\bm{Y}_{i}(\bm{\theta}_{0h})+\alpha_{h}^{(k+1)}\\{\bm{m}_{2h}^{(k+1)}\\}^{T}\bm{Z}_{i}(\bm{\theta}_{0h})\right\\}+\sum_{i=1}^{m_{0}}\epsilon_{ih}$

with

$\displaystyle\epsilon_{ih}=\alpha_{h}^{(k+1)}\left\\{\beta_{h}^{(k+1)}\\{(\bm{\theta}_{1h}^{(k+1)}-\bm{\theta}_{0h})^{3}\\}^{T}\bm{U}_{ih}(\bm{\theta}_{1h}^{(k+1)})+(1-\beta_{h}^{(k+1)})\\{(\bm{\theta}_{2h}^{(k+1)}-\bm{\theta}_{0h})^{3}\\}^{T}\bm{U}_{ih}(\bm{\theta}_{2h}^{(k+1)})\right\\}$

and the rest quantities defined in the same way before. Let $\bar{\bm{t}}$ be

$\displaystyle\left(\alpha_{1}^{(k+1)}-\alpha_{01},\ldots,\alpha_{m_{0}-1}\\{\bm{m}_{11}^{(k+1)}\\}^{T},\ldots,\alpha_{m_{0}}^{(k+1)}\\{\bm{m}_{1m_{0}}^{(k+1)}\\}^{T};\alpha_{1}^{(k+1)}\\{\bm{m}_{21}^{(k+1)}\\}^{T},\ldots,\alpha_{m_{0}}^{(k+1)}\\{\bm{m}_{2m_{0}}^{(k+1)}\\}^{T}\right)^{T}$

Then
$\sum_{i=1}^{n}\delta_{i}=\sum_{i=1}^{n}\bar{\bm{t}}^{T}\bm{b}_{i}+\epsilon_{n}$
with $\epsilon_{n}=o_{p}(1)+o_{p}(n)\bar{\bm{t}}^{T}\bar{\bm{t}}$. By the same
arguments, we have $0\leq
2\bar{\bm{t}}^{T}\sum_{i=1}^{n}\bm{b}_{i}-n\bar{\bm{t}}^{T}\bm{B}\bar{\bm{t}}\\{1+o_{p}(1)\\}+o_{p}(1).$
Therefore $\alpha_{k+1}^{(s)}-\alpha_{0s}=O_{p}(n^{-1/2})$ for
$s=1,\ldots,m_{0}-1$, $\bm{m}_{1t}^{(k+1)}=O_{p}(n^{-1/2})$ and
$\bm{m}_{2t}^{(k+1)}=O_{p}(n^{-1/2})$ for every $t=1,\ldots,m_{0}$. This
completes the proof of lemma.


††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Semi-device-independent quantum random number generator with a broadband
squeezed state of light

Jialin Cheng1 Shaocong Liang1 Jiliang Qin1,2 Jiatong Li1 Zhihui Yan1,2 Xiaojun
Jia1,2<EMAIL_ADDRESS>Changde Xie1,2 Kunchi Peng1,2 1State Key Laboratory of
Quantum Optics and Quantum Optics Devices, Institute of Opto-Electronics,
Shanxi University, Taiyuan 030006, China
2Collaborative Innovation Center of Extreme Optics, Shanxi University,
Taiyuan, Shanxi 030006, China
Email<EMAIL_ADDRESS>
###### Abstract

Random numbers are a basic ingredient of simulation algorithms and
cryptography, and play a significant part in computer simulation and
information processing. One prominent feature of a squeezed light is its lower
fluctuation and more randomness in a pair of orthogonal oriented quadratures,
thus it prompts a significant application in not only quantum information and
quantum precision measurement but also an excellent entropy source for true
random number generation. Here we report a generation of a high-efficiency
semi-device-independent quantum random number based on a broadband squeezed
light, where a reliable randomness source is unnecessary and a noisy local
oscillator is allowed for homodyne detection. The equivalent generation of
private random bits is at a rate of 580.7 Mbps. In addition, the use of
squeezed light at 1.3 $\upmu$m enables the transmission of entropy sources and
local oscillators at the metropolitan scale, thus expanding the potential
applications of quantum random number generators based on non-classical state
of light.

Quantum information, Nonlinear optics, Quantum random number generator

## INTRODUCTION

As an important resource, random number plays an important role in many
fields, such as statistical sampling, computer simulation, cryptography and
lotteries. The random numbers based on classical computer algorithms are
usually known as pseudo-random numbers. In order to generate true random
numbers, the inherent randomness of quantum effects is exploited, which is
usually called quantum random number generation 7Herrero ; 8Ma . Compared with
pseudo-random number generators based on computer algorithms or a lack of
knowledge of the original state, quantum random number generator (QRNG)
exploits a quantum stochastic process as the sources of randomness, which
guarantees the generation of truly random numbers. The usual QRNGs include two
parts, entropy source and corresponding measure-and-process system. The most
simple and effective entropy source is the polarization of arbitrary
polarization state. Series of work about QRNGs have been realized by measuring
the polarization of single photon with a single-photon detector 7Herrero ; 8Ma
. However, the performances of the usual direct-detection system including
dead time and detection efficiency seriously limit the rate of random number
generation. Then several ways have been proposed and realized to enhance the
random number generation rate, such as measuring the temporal or spatial mode
of a photon Yan10 , the distribution of photon number of a laser pulse Apple
and so on Wahl ; LouQ . The quadrature of the quantum state or the phase noise
of the laser could be a good source for random number generation.

The usage of real quantum state as an entropy source can obviously enhance the
security of random number generation obviously and generate true random
numbers 8Ma . Thus, there are multifarious QRNGs based on various quantum
resources, in which the most high-profile is device-independent (DI) QRNG
Pironio ; Bierhorst ; Liu2 ; LiuW ; LiMH ; Acin . A fully DI QRNG can be
achieved solely by observing a violation of the Bell’s inequality Bell ; CHSH
, which does not rely on any physical implementations and the
characterizations of the realistic instruments 8Ma . In the DI QRNG protocol,
the loophole-free Bell violation implies that the predetermined quantum
measurements are performed while the outcomes are fundamentally unpredictable
Pivoluska ; WuD , for example, the typical Bell-type game Liu2 . This protocol
does not assume anything inside the devices used apart from space-like
separation and free choice TBBTC , such approach can be seen as DI. However,
the low generation rate and demanding physical conditions seriously limit the
practical application of DI QRNG. Other QRNGs, based on the conventional
protocols, are usually achieved in prepare-and-measure scenarios. In these
protocols, some devices for preparation or measurement are assumed to be
trusted in some condition. For instance, the source-independent (SI) QRNG
13Marangon ; 14Avesani ; ZhengZ ; LinX ; AvesaniM  and measurement-device-
independent (MDI) QRNG 17Cao ; 258Piran ; YuY  need complete characterizations
of the measurement device and source respectively, which are commonly called
semi-device-independent (SDI) QRNG, while there are also protocols based on
other weak assumptions, such as the dimension of the produced quantum states
LiHWY ; LunghiT , their overlap BraskJB ; TebyanianH , their energy HimbeeckTV
; RuscaD , and so on. By contrast, SDI QRNG generates faster random bits with
less stringent experimental requirements. In this condition, the rate of
generated random number is significantly improved by appropriately assuming
and modeling the source or the measurement device, which may meet the
requirements for practical applications.

With the emergence of SDI quantum random number protocols, high-speed QRNGs
can now be designed based on quadratures measurement of continuous-variable
(CV) quantum states. Because the high-performance macroscopic phtodetector is
used for the measurement, the generation rate of QRNG can be greatly improved
in this frame 253Weed . The common and economical way is measuring the
quadrature of the vacuum or a squeezed state with a high-efficiency balanced
homodyne detector (BHD) as shown in Fig. 1. The vacuum state can be easily
prepared and measured and the measurement performance can also be simply
compensated by optimizing the local oscillator (LO) power. By allocating the
analog result of a single measurement from the homodyne detector to multiple
bits during digitization, the system can generate multiple random bits per
measurement, enabling high-speed quantum random number generation with
relatively simple equipment. The quadratures can be measured via BHD with a
strong LO, and the resulting output current is filtered, digitized and
recorded with an oscilloscope (OSC). The finite resolution of the OSC limits
the discretization of the continuous variable measurements, thus it is
necessary to divide the measurement range of the OSC into $2^{n}$ bins.
Considering a trial with a total of $N$ measurements, the raw random numbers
are obtained when $N$ outcomes fall into different bins that are coded with
$n$ bits ($n$ is determined by the resolution of the OSC). The discretization
of quantum measurement corresponds to the partitions ${\\{(k\delta
o,(k+1)\delta o]\\}}_{k=-\infty}^{k=+\infty}$ with integer $k$, where $\delta
o$ is the precision of the quadrature $\hat{O}$
($\hat{O}=\\{\hat{P},\hat{Q}\\}$) measurement. In fact, this is a process that
Alice applies the coarse-grained POVMs ${\\{\hat{O}^{k}_{\delta o}}\\}$ with
elements $\hat{O}^{k}_{\delta o}={{\int_{k\delta o}^{(k+1)\delta
o}}}do|o\rangle\langle o|$ and stores the outcomes $o_{k}$. Compared with the
vacuum state, a squeezed state with an orthogonal oriented quadrature noise
below the shot noise limit (SNL) can also be used as an entropy source for
QRNG 16Michel ; Zhu3 . Furthermore, the entropy of the anti-squeezed
quadrature of the squeezed state is obviously larger than that of the vacuum
state, which means more randomness can be extracted when the noise of the
anti-squeezed quadrature is utilized for the QRNG as shown in Fig. 1c- 1d.

Figure 1:  The generation of raw random numbers with different entropy sources
in CVs. a-b, the vacuum state as the entropy source. c-d, A squeezed state as
the entropy source. The signal field (the vacuum or squeezed state) and a
strong LO interfere on a 50/50 beam splitter. Then the difference signal from
two photodetectors is digitized and recorded with an OSC. A flatter and wider
noise distribution of the anti-squeezed quadrature of the squeezed state
contains more randomness that can be extracted.

In this article, we experimentally demonstrate a SDI QRNG exploiting a
broadband squeezed state of light. Based on a short optical parametric
amplifier (OPA), a squeezed state ranging from 3.0 MHz to about 200.0 MHz is
obtained and used as the entropy source of the QRNG. Moreover, we introduce a
security analysis of the LO (as a part of the measurement device) by parameter
measurement of detection system, real-time monitoring of the LO fluctuations
to eliminate some eavesdropper’s interference and attacks on the measurement
device by means of a noisy LO. The usage of the untrusted source and noisy LO
at 1.3 $\upmu$m greatly improves the portability of the SDI QRNG. Therefore,
its security is higher than usual SI QRNGs and the amount of extractable
private randomness per unit bandwidth is significantly improved.

## RESULTS

Figure 2:  The entropies of different states as functions of the noise level
and resolution. a, The dashed line ($H_{\rm squ}$) represents the lower bound
on the conditional min-entropies of the anti-squeezed quadrature of the
squeezed state, the light blue ($H_{\rm min}$) and orange curves ($H_{\rm
low}$) represent the min-entropy of the thermal state and the lower bound on
the conditional min-entropy of the thermal state respectively. All the curves
are estimated with a precision $\delta=0.01536$, and the noise levels are set
in units of shot noise. The noise level represents the additional noise for
thermal state and pure quantum noise for squeezed state, respectively. b, The
lower bound on the conditional min-entropy versus the resolution of the
oscilloscope for the vacuum state, a 3.8 dB squeezed state and 8.0 dB squeezed
state. The precision is set as $\delta=0.01536$. The marks show the amount of
extractable secure bits in our experiment.

The security of QRNG and entropic uncertainty principle

To ensure the security of a QRNG, it is paramount to safeguard the generated
quantum random numbers against eavesdropping attempts. For instance, when
Alice employs a QRNG, she must consider potential eavesdropper Eve, who may
employ both classical and quantum methods to eavesdrop the generated random
numbers.

The overall security of a QRNG can be defined as follows. Let $\rho_{OE}$
denote the classical-quantum state of Alice’s random bits $O$ and a possible
eavesdropper $E$. Such a state can always be written as

$\rho_{OE}={\textstyle\sum_{o\in O}}p(o)|o\rangle\langle
o|\otimes\rho_{E}^{o},$ (1)

where $p(o)$ is the probability distribution of the random bits,
$\rho_{E}^{o}$ is the state of the adversary’s system given that $O$ takes the
value $o$. For $\epsilon>0$, we then call the protocol $\epsilon$-secure if
for any eavesdropper $E$

$\frac{1}{2}\left\|\rho_{OE}-\rho_{U}\otimes\rho_{E}\right\|_{1}\leq\epsilon$
(2)

holds. Here, $\left\|\cdot\right\|_{1}$ is the trace norm, $\rho_{U}$ is the
uniform distribution over $O$ and $\rho_{E}$ is the reduced state of
$\rho_{OE}$.

The entropic uncertainty principle (EUP), which originates from the well-known
Heisenberg uncertainty principle, establishes a lower bound on the (quantum)
conditional min-entropy of measuring quadrature for a CV quantum state 9Furrer
; 10Berta ; 11Vallone ; 23Furrer . In fact, randomness determined by general
min-entropy does not eliminate the possibility of eavesdropping, which may
lead to partial insecure randomness offered by eavesdropper in the bits after
post-processing. However, the conditional min-entropy gives the amount of
private and secure randomness in the presence of eavesdroppers 13Marangon ;
16Michel . In accordance with the convention of security analysis in QRNG,
Alice denotes the legitimate user generating random numbers, while Eve, the
malicious eavesdropper, tries to guess correctly the outcomes of measurements
by adopting an optimal strategy. The EUP plays a central role in SI QRNG
protocol for CV quantum systems. The early EUP ColesPJ

$\displaystyle H(P)+H(Q)\geq\log_{2}{\frac{1}{c_{0}}}$ (3)

originates from the famous Heisenberg uncertainty principle $\bigtriangleup
P\cdot\bigtriangleup
Q\geq\frac{1}{2}\left|\left\langle\left[P,Q\right]\right\rangle\right|$, where
the term $\frac{1}{c_{0}}$ quantifies the complementarity of the observables.
In Ref. 9Furrer , the form of EUP for QRNG described by quantum conditional
entropy is given by

$\displaystyle H_{\min}(P_{\delta p}|E)+H_{\max}(Q_{\delta
q})\geq-\log_{2}c(\delta q,\delta p)$ (4)

with

$c(\delta q,\delta p)=\frac{1}{2\pi}\delta q\delta
pS_{0}^{\left(1\right)}(1,\frac{\delta q\delta p}{4})^{2}$ (5)

where $H_{\min}(P_{\delta p}|E)$ is the conditional min-entropy of $P_{\delta
p}$. $H_{\max}(Q_{\delta q})$ is the max-entropy expressing Alice’s lack of
knowledge about outcomes of measuring quadrature $\hat{Q}$. Here $\hat{P}$,
$\hat{Q}$ are defined as data quadrature and check quadrature, respectively.
The term $c(\delta q,\delta p)$ denotes “incompatibility” of the two
quadratures $\hat{P}$ and $\hat{Q}$. $\delta p$ and $\delta q$ represent the
measurement precision of quadratures $\hat{P}$ and $\hat{Q}$ in phase space
with an actual digital device, and normally $\delta p=\delta q=\delta$ is set
for the continuity and convenience of quadratures measurement. $S_{0}^{(1)}$
is the 0th radial prolate spheroidal wave function of the first kind 9Furrer .

The smooth min- and max-entropies are usually exploited in operational
settings, since they result from the non-smoothed versions by an optimization
procedure over states close to the original state 23Furrer . The generalized
residual hash lemma indicates the existence of a family of hash functions such
that for any state $\rho_{AE}$, secure random numbers that are almost uniform
and independent of $E$ can be obtained. Therefore, exploiting the variation of
the conditional min-entropy, smooth min-entropy, the ultimate generation rate
of the secure random bits is further bounded as 13Marangon ; 24Eberle :

$\displaystyle H_{\min}^{\epsilon}(P_{\delta p}|E)$ $\displaystyle\geq$
$\displaystyle-\log_{2}c(\delta q,\delta p)-H_{\max}(Q_{\delta
q}))-\frac{1}{\sqrt{n_{p}}}\Delta$ (6) $\displaystyle\equiv$ $\displaystyle
H_{\rm low}^{\epsilon}(P_{\delta p}|E)$

with

$\Delta=4\sqrt{\log_{2}(\frac{2}{\epsilon^{2}})}\log_{2}(2^{1+\frac{H_{\max}(Q_{\delta
q})}{2}}+1),$ (7)

where $H_{\min}^{\epsilon}(P_{\delta p}|E)$ is the smooth min-entropy, $n_{p}$
is the number of measurements for quadrature $\hat{P}$ and $\epsilon$ is the
security parameter. The protocol is called $\epsilon$-secure, which means that
it is $\epsilon$-close to true random number with uniform distribution and
independent of $E$.

Advantages of a broadband squeezed light as entropy source

Following the SI QRNG protocol 13Marangon ; 16Michel  and ignoring hardware
and post-processing, the key factors affecting the generation rate of QRNG
mainly include the precision of coarse-grained measurement $\delta$, the
bandwidth of measured signal and the purity of quantum states (more details
can be found in Supplementary Note 3). Compared with the vacuum, exploiting a
broadband squeezed light as the entropy source can enhance the generation rate
due to more (secure) randomness from the outcomes of the anti-squeezed
quadrature measurements. As shown in Fig. 1, the outcome of anti-squeezed
quadrature is a flatter and wider distribution. Moreover, an extended
measurement bandwidth can significantly improve the generation rate of random
numbers (see Supplementary Note 3). In short, it is an appreciable improvement
by exploiting a wider-bandwidth and higher-level squeezed state as the entropy
source. The min-entropy of the thermal state, (quantum) conditional min-
entropy of the thermal state and the conditional min-entropy of the anti-
squeezed quadrature of squeezed state as functions of the noise level have
been theoretically analyzed respectively, as shown in Fig. 2a. The noise level
represents the additional noise for thermal state and pure quantum noise for
squeezed state, respectively. More details about the additional noise can be
found in Supplementary Note 3. With the increase of the additional noise in
the measured quadrature, the min-entropy of the thermal state $H_{\rm min}$
increases, in which the quadrature seems to contain more randomness. However,
the increased randomness cannot be used to generate secure random numbers
because the additional noise decreases the purity of the state and may be
insecure. The additional noise may be correlated with Eve and give side
information to Eve. The content of secure randomness can be estimated by the
conditional min-entropy, therefore it is clear that the extractable secure
randomness of the thermal state $H_{\rm low}$ decreases with increasing noise
level. In contrast, the increasing noise of the anti-squeezed quadrature of
the squeezed state, which comes from the property of the quantum state and is
pure quantum noise so that it can generate secure random bits. The lower bound
on the conditional min-entropy of the anti-squeezed quadrature of the squeezed
state is the same as the min-entropy of the thermal state with identical noise
level.

Without regard to the sample size, the lower bound on the conditional min-
entropy as a function of the resolution of the OSC is shown in Fig. 2b, where
the vacuum, a 3.8 dB-squeezed state (in our case) and an 8.0 dB-squeezed state
are chosen as entropy sources, respectively. The precision is set as
$\delta=0.01536$ in our experiment. It can be seen that no secure random bits
can be extracted when the resolution is very low for all these quantum states.
When the resolution increases to a certain value, the lower bound on the
conditional min-entropy increases linearly with the resolution.

Figure 3:  Schematic of experimental configuration for SDI QRNG. The whole
space can be divided into an untrusted space and a private space. A broadband
squeezed state is generated form the OPA, and detected with a BHD and the LO
transmitted from the untrusted space. The check quadrature, electronic noise
and the LO fluctuation are used to calibrate the conditional min-entropy. The
data quadrature is used to generate raw random bits. Very few terminal random
bits are injected into the locking box and chopper to determine the quadrature
type and measure electronic noise respectively. EC: electric cable; DC: data
cable; HR: mirror with high reflectivity; BS: 50/50 beam splitter; DBS:
dichroic beam splitter; OPA: optical parametric amplifier; BHD: balanced
homodyne detector; BA: broadband amplifier; LB: locking box; OSC:
oscilloscope; EEB: entropy extraction box; EB: entropy estimation box; NSU:
noise source unit; EOM: electro-optic modulator.

The security analysis of LO

There is not a delicate characterization of the noisy LO in CV SI QRNG
framework. Simply setting a power bound on the LO and incessantly monitoring
the LO power may not eliminate most security threats to the measurement
device. A notable threat is that Eve modulates the LO and adds extra
controlled classical noise as side information to the measured quadrature
furtively, which leads to an overestimation of the content of private
randomness. Although there is no such threat for a ideal BHD 27Huang . The
security of ultimate private bits is also compromised since Alice cannot
notice the existence of Eve by monitoring LO power in this case.

In order to close this backdoor in the measurement device of QRNG, it is
necessary to record real-time LO fluctuations and perform a further
calibration of the lower bound on the conditional min-entropy in CV SI QRNG
scenario, which ensures that this protocol is capable of eliminating some
side-channel attacks on the measurement device. On the basis of the theory
derived in Supplementary Note 4, we quantify the untrusted noise caused by LO
fluctuation in units of shot noise. The total noise $\sigma_{\rm M-vac}^{2}$,
$\sigma_{\rm M-squ}^{2}$ and $\sigma_{\rm M-ant}^{2}$ measured on quadratures
of the vacuum, the squeezed quadrature $\hat{Q}$ and anti-squeezed quadrature
$\hat{P}$ of the squeezed state are respectively given by:

$\displaystyle\sigma_{\rm M-vac}^{2}$ $\displaystyle=$
$\displaystyle\sigma_{\rm SNL}^{2}+\sigma_{\rm LO}^{2}+\sigma_{\rm E}^{2},$
$\displaystyle\sigma_{\rm M-squ}^{2}$ $\displaystyle=$
$\displaystyle\sigma_{\rm squ}^{2}+\sigma_{\rm LO}^{2}+\sigma_{\rm E}^{2},$
$\displaystyle\sigma_{\rm M-ant}^{2}$ $\displaystyle=$
$\displaystyle\sigma_{\rm ant}^{2}+\sigma_{\rm LO}^{2}+\sigma_{\rm E}^{2}.$
(8)

The terms $\sigma_{\rm SNL}^{2}$, $\sigma_{\rm LO}^{2}$, and $\sigma_{\rm
E}^{2}$ denote the measured pure vacuum noise, the untrusted noise introduced
by LO fluctuation and the intrinsic electronic noise of the measurement
device, respectively. $\sigma_{\rm squ}^{2}$, $\sigma_{\rm ant}^{2}$ denote
the noise of squeezed and anti-squeezed quadratures, respectively.

The electronic noise can be measured by blocking the signal field and LO.
Through carefully measuring the relevant parameters of the homodyne detector
and recording the LO fluctuations in real time, the untrusted noise introduced
by LO fluctuation can be obtained. Then the lower bound on the conditional
min-entropy can be further calibrated.

Experimental setup.

The schematic of experimental setup of SDI QRNG with a broadband squeezed
state is shown in Fig. 3, which consists of two parts, an untrusted space and
a private space. The two spaces can be connected by standard telecommunication
fiber for the portability of the QRNG. The untrusted space contains an
uncertified broadband squeezed state and an incompletely characterized LO. The
LO can be noisy since Eve may modulate the LO using an electro-optic modulator
(EOM) with a controlled noise source unit (NSU) for side information. The
laser at 1342 nm and 671 nm from a dual-wavelength laser (not shown in Fig. 3)
are injected into an OPA to generate a broadband squeezed state of light.
Eventually, the measured result for the noise power of broadband squeezed
state of light is shown as Fig. 4a. A $\hat{Q}$-squeezed state with $6.5\pm
0.2$ dB squeezing at 3.0 MHz and $3.8\pm 0.2$ dB squeezing at 200.0 MHz is
prepared. More details about the OPA can be found in the Methods section.

The second part is a trusted data acquisition and processing system in the
private space. The LO fluctuations are monitored and analyzed in real time via
a detector in order to resist the attacks on the LO. The relative phase
between the LO and signal field is locked at $\uppi/2$ and randomly switched
to 0 via a lock box (LB) according to the random seed, which determines the
obtained data type, the data quadrature $\hat{P}$ for raw random numbers and
the check quadrature $\hat{Q}$ for conditional min-entropy estimation
respectively. The chopper is used to randomly block the signal field and LO to
measure the electronic noise. The signals are filtered, amplified and acquired
by filters, a broadband amplifier (BA) and an OSC. The untrusted noise
introduced by LO fluctuation, the measured check quadrature and electronic
noise are fed into the conditional min-entropy estimation box (EB) to
calibrate the entropy extraction box (EEB). The measured data quadrature is
fed into the EEB to extract secure random bits. A random seed from the
terminal random bits in the memory is injected into the locking box and
chopper to determine the relative phase of LO and signal field and measure
electronic noise respectively.

Figure 4:  Experimental results. a, The measured noise power of the broadband
squeezed state of light. Black curve: electronic noise level. Blue curve:
shot-noise-limit (SNL) with the total LO power of 7.8 mW. Red curve: noise
power of squeezed quadrature amplitude $\hat{Q}$. Green curve: noise power of
anti-squeezed quadrature phase $\hat{P}$. The measured frequency is from 3.0
MHz to 200.0 MHz, and a $6.5\pm 0.2$ dB of squeezing at 3.0 MHz and $3.8\pm
0.2$ dB of squeezing at 200.0 MHz below the corresponding SNL is observed. b,
The results of a typical run of the NIST test. In the case of multiple tests
in a category, the smallest have been reported.

Compared with the previous SI QRNG protocol 16Michel , a security analysis of
the noisy LO is used for eliminating some Eve’s interference and attacks on
the detectors in our protocol, which acts as an essential process for
improving overall security. In the beginning of the experiment, the various
noise levels of the LO with different powers over the whole analysis frequency
range are measured and the relevant parameters of the devices inside private
space are calibrated, and then the precision $\delta$ is estimated. The max-
entropy $H_{\max}\left(Q_{\delta q}\right)$ is calculated by the recorded data
of check quadrature using frequentist estimator, while the data of the other
quadrature are used to generate random numbers. This protocol adopts an
extremely conservative approach to security, treating all classical noise and
unnecessary quantum noise as impurities within the input state see
Supplementary Note 3. Any observed “mixing” is attributed to a specific
quantum eavesdropping strategy to safeguard the security of QRNG. According to
real-time feedback of the $H_{\max}\left(Q_{\delta q}\right)$, the electronic
noise of the detector and the untrusted noise caused by the LO fluctuation,
the lower bound on the quantum conditional min-entropy can be estimated. The
Toeplitz-matrix hashing algorithm is applied to extract the private random
bits from the raw random numbers. Abort protocol whenever the data are out of
the range of the devices or the homodyne detector is saturated. Finally, the
generation rate of secure random bits of our protocol is 580.7 Mbps.

It is noteworthy that we need to emphasize a trivial assumption that these
internal devices and hardware in the private space of Fig. 3 are absolutely
secure, and there is no possibility that Eve changes the reliable parts of the
experimental setup (the reflectivity of the beam splitter, the quantum
efficiency of photodiodes and the circuitry part of the detector, etc.) and
even surreptitiously installs backdoor malware. However, Eve knows all of
these parameters and the running process of the setup and protocol, and is
capable of regulating the temperature, the electric and magnetic fields around
the private space to control and manipulate the electronic noise for side
information during the generation of random numbers.

## DISCUSSION

Unlike most CV quantum key distribution (QKD) schemes FhX , we use continuous
LO instead of optical pulses, and there is no need for synchronization of RF
signals for source preparation and quadrature measurement, which greatly
simplify the experimental system and avoid many attacks exist in CV-QKD. In
addition, the LO fluctuation is inappreciable in our SDI QRNG, whose effect on
the measured noise is small to distinguish with our parameters of the homodyne
detector as calculated in the Methods section. However, if the ratio of the
reflectivity to transmission is far from 1:1 for a substandard BHD, very
little classical fluctuation of the LO will give much side information to Eve.
Alice has no way to combat such an eavesdropping strategy by monitoring LO
power and aborting the protocol with noticed detector saturation. Moreover, we
have taken note of the antenna attack problem which has been recently proposed
as an area of concern SmithPR . It may be beneficial to eliminate this antenna
attack by monitoring the detector for saturation and random measurements of
the detector’s electronics noise.

SI QRNG shows a good prospect in both classical and quantum communication
networks. The securities and generation rates of SI QRNGs need more detailed
description and further improvement. Our SDI QRNG protocol aims to improve not
only security but generation rate of QRNG. A broadband squeezed light is
exploited to improve generation rate with a maintained measured bandwidth.
Based on the traditional SI protocol, a security analysis of the LO is
introduced for increasing the overall security. Finally, the generation rate
of secure random bits is estimated to 580.7 Mbps, which can meet the
requirements of an ordinary short-range QKD system and a quantum secure direct
communication system Hujy ; ZhangHET . In addition, the use of squeezed light
at 1.3 $\upmu$m enables the transmission of entropy sources and LOs at the
metropolitan scale, and makes it possible for the generation of the random
numbers certified by Bell’s theorem. Recently, microring resonators are
exploited as optical parametric oscillator LuJ ; BWA ; OkY  and some quantum
light sources have been produced on a photonic chip Zhaoy ; Cas ; Singh ;
luoWET ; LiXET , which mean a smaller SDI QRNG and wider squeezing bandwidth
are achievable. Meanwhile, there have been some reports of random number
generations with a photonic integrated chip BaiBB ; AbellanC  and an optical
chip for self-testing quantum random number generation with a silicon
photomultiplier LeoneN . On-chip quantum random number generation is a major
leap towards miniaturization and integration of QRNG. Moreover, a broadband
squeezed light has inestimable advantages in information processing and pure
state preparation, and shows extraordinary potential in quantum information
30Kashiwazaki ; 31Politi ; 32Masada  and quantum measurement 33Yu  due to its
wide bandwidth as well as higher purity at higher frequency.

## METHODS

Details about the experimental setup

The OPA contains a $1\times 2\times 6$ $\rm mm^{3}$ periodically poled KTiOPO4
(PPKTP) crystal and a piezo-actuated concave mirror with a radius of curvature
of 150 mm which works as the output coupler. The length of the cavity is about
10 mm. The front surface of the PPKTP is coated with a film of reflectivity
$R>99.9\%$ at 1342 nm and transmission $T=85.3\%$ at 671 nm. The rear surface
of the PPKTP is coated with anti-reflection (AR) at both 1342 nm and 671 nm.
The output coupler is coated with reflectivity $R=88.2\%$ at 1342 nm and
$R=99.9\%$ at 671 nm. The finesse of the OPA is about 25. Compared with our
previous squeezed light sources 28Hou ; YanZJia ; Zuoxj , the extended
bandwidth of the OPA is achieved mainly by redesigning the structure of the
OPA and shortening the cavity length. With a film-coated PPKTP crystal, the
usual input mirror can be omitted in this semi-monolithic configuration
5ZhouYY . Then the half width at half maximum of the OPA can be estimated to
about 202 MHz because the FSR of the cavity is approximately 10.1 GHz. When
the length of the cavity and working condition are well controlled with
corresponding locking systems, the intracavity losses are $0.3\%$ and the
parametric gain is about 20 at the pump power of 96 mW. Eventually a bright
squeezed state with the bandwidth of over 200.0 MHz is obtained.

The homodyne detection is achieved with a homemade detector included two
photodiodes FD150 (Fermionics Opto-Technology company) to meet the
requirements of wide measured bandwidth and high quantum efficiency for
broadband squeezed states. The interference signal of the signal field and the
LO is exploited as the error signal to achieve the locking of the quadrature
phase (i.e., the data quadrature), and provides a high-speed switching signal
to the fiber-coupled electro-optic modulator (not shown in Fig. 3) for
achieving quadrature amplitude switching according to the random seed. The
direct current of the detector is the interference signal, while the
alternating current is used to acquire data for further steps. Furthermore,
all frequencies of the modulating signals for frequency stabilization and
phase locking are set below 1 MHz to offer an appropriate frequency range for
the SDI QRNG protocol. For instance, the frequency stabilization of the mode
cleaners is achieved by exploiting a homemade locking box with a modulating
signal of 10 kHz. The propagation efficiency, fringe visibility, and
photodiodes quantum efficiency are 0.99, 0.99 and 0.88, respectively. The
electronic noise level is about -13 dB compared to the SNL with a LO intensity
of 7.8 mW at 200.0 MHz. The results of measured noise power of the generated
broadband squeezed state of light in Fig. 4a are obtained after filtering out
low-frequency noise via a 3.0 MHz high-pass filter. There is substantial
classical noise including relaxation oscillation noise and modulating signals
existing in the range of 0 to 3 MHz, thus the measured results in this range
have to be discarded for the generation of secure random bits. On the other
hand, only squeezing of 3.8 dB in the full frequency range (from 3.0 MHz to
200.0 MHz) is used to ensure the security of QRNG in the most conservative
situation. Otherwise, the amount of secure randomness will be overestimated.

Experimental results and data processing

The nominal resolution of the OSC is over 11 bits. The signals are acquired at
a rate of 1.5 GSamples per second and then downsampled to 200 MSamples per
second. The relative phase between the LO and signal beam is randomly switched
to quadrature $\hat{Q}$ and collect the check data in every 20 $\mu$s, and the
electronic noise is randomly measured in every 0.4 s. We set a ratio of 1/20
between the total measurements and the check measurements, and so does the
electronic noise.

In a round of random number generation, the measured vacuum noise contains
(5.75 $\pm$ 0.02)% of electronic noise and less than $10^{-4}$ % of the
untrusted noise introduced by the LO fluctuation. The precision in phase space
is estimated to be $\delta=\left(1.536\pm 0.002\right)\times 10^{-2}$, a
detailed process can be found in Supplementary Note 5. The duration of
randomly measuring check quadrature and electronic noise are set as 1$\upmu$s
and $0.02$ s respectively. $2\times 10^{5}$ samples are obtained and the
security parameter is set as $\epsilon=10^{-6}$. The lower smooth conditional
min-entropy is estimated to be $H_{\min}^{\epsilon}(P_{\delta p}|E)=7.21\pm
0.04$ bits per sample. Taking into account the effective number of bits of the
OSC, we use a conservative bit depth of 6 bits, and the entropy is reduced to
$3.32\pm 0.04$ bits. A $2$ Mb $\times$ $656$ kb Toeplitz matrix hashing is
applied to extract random bits. In addition, we need to consume very few
private bits for the data-type switching. Finally, the equivalent generation
of private random bits is at a rate greater than 580.7 Mbps.

For an undisturbed LO, it can be regarded as a coherent state. However, if Eve
modulates the LO, the LO fluctuations will be superimposed on the measured
noise with a ratio of $K_{1}/K_{2}$ (see Supplementary Note 4 for more
details). For instance, if the reflectivity and transmission of the
substandard beam splitter are about $21.2\%$ and $78.9\%$ respectively, and
the LO noise modulated by Eve increases by 3 dB for an eavesdropping strategy,
then the measured vacuum noise is not pure and contains $50\%$ classical
noise, i.e., $\sigma_{\rm LO}^{2}=\sigma_{\rm SNL}^{2}$. Eventually, one can
estimate that $31.02\%$ of the extracted random bits are not secure and
private.

In order to check for the randomness in terminal private bits, we test them
with the NIST suite 44Bass . The results of a typical run are reported in Fig.
4b. In the case of multiple tests in a category, the smallest has been
reported. The private bits pass all the NIST statistical tests.

In addition, the squeezing of CV quantum state rapidly decays with losses, and
the generation rate decreases with transmission distance. However, as long as
the squeezed state can be maintained in transmitting process, the bit rate of
a squeezed state is better than that of the vacuum. With the single mode loss
of 0.35 dB per kilometer at 1.3 $\upmu$m in fiber, the squeezing level reduces
to 1.4 dB at 200 MHz when the transmission distance is about 9.27 km.
Therefore, there is still obvious advantage over the vacuum in generation
rates.

DATA AVAILABILITY

These data for experimental results can be found at
https://github.com/chengjialin007/bssQRNG. The code is available from the
corresponding author on request.

ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China
(Grants No. 61925503, No. 62122044, No. 12147215, and No. 11834010), the Key
Project of the National Key R&D program of China (Grant No. 2016YFA0301402),
the Program for the Innovative Talents of Higher Education Institutions of
Shanxi, the Program for the Outstanding Innovative Teams of Higher Learning
Institutions of Shanxi and the fund for Shanxi “1331 Project” Key Subjects
Construction.

COMPETING INTERESTS

The authors declare no competing interests.

AUTHOR CONTRIBUTIONS

X.J. and C.X. conceived the original idea. J.C., S.L., Z.Y., X.J., and K.P.
designed the experiment. J.C., S.L., J.Q., and J.L. constructed and performed
the experiment. J.C., S.L., and X.J. accomplished theoretical calculation and
the data analysis. J.C., Z.Y., X.J., C.X., and K.P. wrote the paper. All the
authors reviewed the manuscript. J.C. and S.L. contributed equally to this
work.

## References

  * (1) Herrero-Collantes, M. & Garcia-Escartin, J. C. Quantum random number generators. Rev. Mod. Phys. 89, 015004 (2017).
  * (2) Ma, X., Yuan, X., Cao, Z., Qi, B. & Zhang, Z. Quantum random number generation. npj Quantum Inf. 2, 16021 (2016).
  * (3) Yan, Q., Zhao, B., Liao, Q. & Zhou, N. Multi-bit quantum random number generation by measuring positions of arrival photons. Rev. Sci. Instrum. 85, 103116 (2014).
  * (4) Applegate, M. et al. Efficient and robust quantum random number generation by photon number detection. Appl. Phys. Lett. 107, 071106 (2015).
  * (5) Wahl, M. et al. An ultrafast quantum random number generator with provably bounded output bias based on photon arrival time measurements. Appl. Phys. Lett. 98, 171105 (2011).
  * (6) Lou, Q. et al. Quantum random number generator based on single-photon emitter in gallium nitride. Opt. Lett.  45, 4224-4227 (2020).
  * (7) Pironio, S. et al. Random numbers certified by Bell theorem. Nature 464, 1021-1024 (2010).
  * (8) Bierhorst, P. et al. Experimentally generated randomness certified by the impossibility of superluminal signals. Nature 556, 223-226 (2018).
  * (9) Liu, Y. et al. Device-independent quantum random-number generation. Nature 562, 548-551 (2018).
  * (10) Liu, W. et al. Device-independent randomness expansion against quantum side information. Nat. Phys. 17, 448-451 (2021).
  * (11) Li, M.-H. et al. Experimental realization of device-independent quantum randomness expansion. Phys. Rev. Lett. 126, 050503 (2021).
  * (12) Acin, A. & Masanes L. Certified randomness in quantum physics. Nature 540, 213–219 (2016).
  * (13) Bell, J. S. On the Einstein-Podolsky-Rosen paradox, Physics 1 195 (1964).
  * (14) Clauser, J. F., Horne, M. A., Shimony, A. & Holt, R. A. Proposed experiment to test local hidden-variable theories. Phys. Rev. Lett. 23, 880 (1969).
  * (15) Pivoluska, M. et al. Semi-device-independent random number generation with flexible assumptions. npj Quantum Inf.  7, 50 (2021).
  * (16) Wu, D. et al. Closing the locality and detection loopholes in multiparticle entanglement self-testing. Phys. Rev. Lett. 128, 250401 (2022).
  * (17) The BIG Bell Test Collaboration. Challenging local realism with human choices. Nature 557, 212-216 (2018).
  * (18) Marangon, D. G., Vallone, G. & Villoresi, P. Source-device-independent Ultra-fast Quantum Random Number Generation. Phys. Rev. Lett. 118, 060503 (2017).
  * (19) Avesani, M., Marangon, D. G., Vallone, G. & Villoresi, P. Source-device-independent heterodyne-based quantum random number generator at 17 Gbps, Nat. Commun. 9, 5365 (2018).
  * (20) Zheng, Z. et al. Bias-free source-independent quantum random number generator. Opt. Express  28, 22388-22398 (2020).
  * (21) Lin, X. et al. Imperfection-insensitivity quantum random number generator with untrusted daily illumination. Opt. Express  30, 25474–25485 (2022).
  * (22) Avesani, M., Tebyanian, H., Villoresi, P. & Vallone, G. Unbounded randomness from uncharacterized sources. Commun. Phys. 5, 273 (2022).
  * (23) Cao, Z., Zhou, H. & Ma, X. Loss-tolerant measurement-device-independentquantum random number generation. New J. Phys. 17, 125011 (2015).
  * (24) Pirandola, S. et al. High-rate measurement-device-independentquantum cryptography. Nat. Photon. 9, 397-402 (2015).
  * (25) Yu, Y. et al. Measurement-device-independent verification of a quantum memory. Phys. Rev. Lett. 127, 160502 (2021).
  * (26) Li, H.-W. et al. Semi-device-independent random-number expansion without entanglement. Phys. Rev. A 84, 034301 (2011).
  * (27) Lunghi, T. et al. Self-testing quantum random number generator. Phys. Rev. Lett. 114, 150501 (2015).
  * (28) Brask, J. B. et al. Megahertz-rate semi-device-independent quantum random number generators based on unambiguous state discrimination. Phys. Rev. Appl. 7, 054018 (2017).
  * (29) Tebyanian, H. et al. Semi-device independent randomness generation based on quantum state’s indistinguishability. Quantum Sci. Technol. 6, 045026 (2021).
  * (30) Himbeeck, van T., Woodhead, E., Cerf, N. J., García-Patrón, R. & Pironio, S. Semi-device-independent framework based on natural physical assumptions. Quantum 1, 33 (2017).
  * (31) Rusca, D. et al. Self-testing quantum random-number generator based on an energy bound. Phys. Rev. Appl. 100, 062338 (2019).
  * (32) Weedbrook, C. et al. Gaussian quantum information. Rev. Mod. Phys. 84, 621 (2012).
  * (33) Zhu, Y., He, G. & Zeng, G. Unbiased quantum random number generation based on squeezed vacuum state. Int. J. Quantum Inf. 10, 1250012 (2012).
  * (34) Michel, T. et al. Real-time source independent quantum random number generator with squeezed states. Phys. Rev. Appl. 12, 034017 (2019).
  * (35) Furrer, F., Berta, M., Tomamichel, M., Scholz, V. B. & Christandl, M. Position-momentum uncertainty relations in the presence of quantum memory, J. Math.Phys. 55, 122205 (2014).
  * (36) Berta, M., Christandl, M., Colbeck, R., Renes, J. M. & Renner, R. The uncertainty principle in the presence of quantum memory. Nat. Phys. 6, 659-662 (2010).
  * (37) Vallone, G., Marangon, D. G., Tomasin, M. & Villoresi, P. Quantum randomness certified by the uncertainty principle. Phys. Rev. A 90, 052327 (2014).
  * (38) Furrer, F., Aberg, J. & Renner, R. C. Min- and max-entropy in infinite dimensions. Commun. Math. Phys. 306, 165-186 (2011).
  * (39) Coles, P. J., Berta, M., Tomamichel, M. & Wehner, S. Entropic uncertainty relations and their applications. Rev. Mod. Phys. 89, 015002 (2017).
  * (40) Eberle, T. et al. Gaussian entanglement for quantum key distribution from a single-mode squeezing source. New J. Phys. 15, 053049 (2013).
  * (41) Huang, W. et al. Practical security analysis of a continuous-variable quantum random-number generator with a noisy local oscillator. Phys. Rev. A 102, 012422 (2020).
  * (42) Xu, F., Ma, X., Zhang, Q., Lo, H.-W. & Pan, J.-W. Secure quantum key distribution with realistic devices. Rev. Mod. Phys. 92, 025002 (2020).
  * (43) Smith, P. R., Marangon, D. G., Lucamarini, M., Yuan, Z. L. & Shields, A. J. Out-of-band electromagnetic injection attack on a quantum random number generator. Phys. Rev. Applied 15, 044044 (2021).
  * (44) Hu, J. Y. et al. Experimental quantum secure direct communication with single photons. Light Sci. Appl. 5, e16144 (2016).
  * (45) Zhang, H. et al. Realization of quantum secure direct communication over 100 km fiber with time-bin and phase quantum states. Light Sci. Appl. 11, 83 (2022).
  * (46) Lu, J. et al. Ultralow-threshold thin-film lithium niobate optical parametric oscillator. Optica 8, 539–544 (2021).
  * (47) Bruch, A. W., Liu, X., Surya, J. B., Zou, C.-L. & Tang, H.-X. Onchip $\chi^{(2)}$ microring optical parametric oscillator. Optica 6, 1361–1366 (2019).
  * (48) Okawachi, Y. et al. Dynamic control of photon lifetime for quantum random number generation. Optica 8, 1458–1461 (2021).
  * (49) Zhao, Y. et al. Near-Degenerate Quadrature-Squeezed Vacuum Generation on a Silicon-Nitride Chip. Phys. Rev. Lett. 124, 193601 (2020).
  * (50) Caspani, L. et al. Integrated sources of photon quantum states based on nonlinear optics. Light Sci. Appl. 6, e17100 (2017).
  * (51) Singh, N., Raval, M., Ruocco, A. & Watts, M. R. Broadband 200-nm second-harmonic generation in silicon in the telecom band. Light Sci. Appl. 9, 17 (2020).
  * (52) Luo, W. et al. Recent progress in quantum photonic chips for quantum communication and internet. Light Sci. Appl. 12, 175 (2023).
  * (53) Li, X. et al. Bright semiconductor single-photon sources pumped by heterogeneously integrated micropillar lasers with electrical injections. Light Sci. Appl. 12, 65 (2023).
  * (54) Bai, B. et al. 18.8 Gbps real-time quantum random number generator with a photonic integrated chip. Appl. Phys. Lett. 118, 264001 (2021).
  * (55) Abellan, C. et al. Quantum entropy source on an InP photonic integrated circuit for random number generation. Optica 3, 989-994 (2016).
  * (56) Leone, N. et al. An optical chip for self-testing quantum random number generation. APL Photon. 5, 101301 (2020).
  * (57) Kashiwazaki, T. et al. Continuous-wave 6-dB-squeezed light with 2.5-THz-bandwidth from single-mode PPLN waveguide. APL Photon. 5, 036104 (2020).
  * (58) Politi, A., Cryan, M. J., Rarity, J. G., Yu, S. & O’Brien, J. L. Silica-on-silicon wave guide quantum circuits. Science 320, 646-649 (2008).
  * (59) Masada, G. et al. Continuous-variable entanglement on a chip. Nat. Photonics 9, 316-319 (2015).
  * (60) Yu, J. et al. Quantum enhanced optical phase estimation with a squeezed thermal state. Phys. Rev. Appl. 13, 024037 (2020).
  * (61) Huo, M. et al. Deterministic quantum teleportation through fiber channels. Sci. Adv. 4, eaas9401 (2018).
  * (62) Yan, Z. & Jia, X. Teleportation goes to Hertz rate. Light Sci. Appl. 12, 167 (2023).
  * (63) Zuo, X. et al. Quantum Interferometer Combining Squeezing and Parametric Amplification. Phys. Rev. Lett. 124, 173602 (2020).
  * (64) Zhou, Y. et al. Quantum secret sharing among four players using multipartite bound entanglement of an optical field. Phys. Rev. Lett. 121, 150502 (2018).
  * (65) Bassham III, L. E. et al. Sp 800-22 rev. 1a. A statistical test suite for random and pseudorandom number generators for cryptographic applications.(2010). URL https://www.nist.gov/publications/statistical-test-suite-random-and-pseudorandom-number-generators-cryptographic.