# Enhanced Monte Carlo Estimation of the Fisher Information Matrix with
Independent Perturbations for Complex Problems

Enhanced Monte Carlo Estimation of the Fisher Information Matrix with
Independent Perturbations for Complex Problems

Xuan Wu

Department of Applied Math and Statistics
The Johns Hopkins University
Baltimore, Maryland 21218, USA

Abstract

The Fisher information matrix provides a way to measure the amount of
information given observed data based on parameters of interest. Many
applications of the FIM exist in statistical modeling, system identification,
and parameter estimation. We sometimes use the Monte Carlo-based method to
estimate the FIM because its analytical form is often impossible or difficult
to be computed in real-world models. In this paper, we review the basic method
based on simultaneous perturbations and present an enhanced resampling-based
method with independent simultaneous perturbations to estimate the Fisher
information matrix. We conduct theoretical and numerical analysis to show its
accuracy via variance reduction from $O(1/N)$ to $O(1/(nN))$, where $n$ is the
sample size of the data and $N$ is a measure of the Monte Carlo averaging. We
also consider the trade-off between accuracy and computational cost.
Key Words: Monte Carlo simulation; Simultaneous perturbation; Fisher
information matrix; Variance reduction

##  1 Introduction

The Fisher Information Matrix (FIM) summarizes the amount of information about
data related to parameters of interest, bearing in mind that the central role
of data analysis is to extract information from data. We can see that the FIM
plays a critical role in the theory and practice of statistical modeling,
including system identification [1], parameter estimation [2], information
theory [3], and other areas. Let $\boldsymbol{Z}$ be the collection of $n$
random vectors $[\boldsymbol{z}_{1},\ldots,\boldsymbol{z}_{n}]$. In the
following discussion, assume that $p(\boldsymbol{Z}|\boldsymbol{\theta})$ is
the density function for the model of interest, where $\boldsymbol{\theta}$
represents the unknown $p\times 1$ parameter vector, and
$L(\boldsymbol{\theta}|\boldsymbol{Z})$ is the log-likelihood function of
$p(\boldsymbol{Z}|\boldsymbol{\theta})$. Then the $p\times p$ FIM is defined
as:

$\boldsymbol{F}(\boldsymbol{\theta})=\mathbb{E}\bigg{[}\frac{\partial
L}{\partial\boldsymbol{\theta}}\frac{\partial
L}{\partial\boldsymbol{\theta}^{T}}\bigg{]}.$

In general, it is difficult to compute the expectation of a product of
nonlinear score functions with multiple parameters. On the other hand, if the
second partial derivatives of $L$ exists, i.e., the Hessian matrix
$\boldsymbol{H}(\boldsymbol{Z}|\boldsymbol{\theta})$, then under specific
regularity conditions [4], the FIM can be represented as:

$\boldsymbol{F}(\boldsymbol{\theta})=-\mathbb{E}\bigg{[}\frac{\partial^{2}L}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{T}}\bigg{]}=-\mathbb{E}[\boldsymbol{H}(\boldsymbol{\theta})],$

which is generally easier to calculate compared with the basic definitional
form.

In that there are no closed forms for $\boldsymbol{F}(\boldsymbol{\theta})$ in
many practical problems and given the importance of the FIM, one way is to use
averages of the estimated Hessian matrix based on pseudo data as the Monte
Carlo estimation of the FIM. Moreover, some research studies present methods
to estimate the FIM in general or specific cases (e.g., [5], [6]). In
particular, we can estimate the FIM via measurements of a Hessian matrix based
on the simultaneous perturbation approach. Spall in [5] discusses a basic
resampling-based approach. Later several papers present improved approaches,
including a feedback-based approach in [7] and a method where information on
some of the FIM is available in [8]. In this paper, we introduce another
enhanced Monte Carlo method: the independent perturbation approach. Ref. [7]
also mentions the independent perturbation approach, but we do a more thorough
analysis here.

In Section 2, we demonstrate theories of the basic resampling-based method and
the independent perturbation method. Simultaneously, we show that the
independent perturbation method reduces the variance of entries of the
estimated FIM by from $O(1/N)$ to $O(1/(nN))$, where $n$ is the sample size of
the data and $N$ is a measure of the Monte Carlo averaging. Section 3
illustrates the efficiency of the enhanced method through two numerical
examples, respectively. The signal-plus-noise example verifies the theoretical
ratio about variance reduction and compares the performance at each diagonal
entries of the estimated FIM for the two methods. The mixture Gaussian example
utilizes the relative error based on the spectral norm to show the obvious
improvement of the independent perturbation method. Section 4 gives a
conclusion and ideas about future work.

##  2 Theoretical Analysis

###  2.1 Standard Method

We summarize aspects of the basic method in [5] below. Assume that
$\boldsymbol{Z}_{\mathrm{pseudo}}(i)$ is a set of $n$ independent random
vectors generated by the Monte Carlo simulation that follow the distribution
with parameters $\boldsymbol{\theta}$, and $\hat{\boldsymbol{H}}_{k|i}$ is the
$k$th estimate of $\boldsymbol{H}(\boldsymbol{\theta})$ at
$\boldsymbol{Z}_{\mathrm{pseudo}}(i)$. Additionally, suppose that
$\boldsymbol{\Delta}_{k|i}=\left[\Delta_{k1|i},\Delta_{k2|i},\ldots,\Delta_{kp|i}\right]^{T}$
is a random vector with a zero mean where the scalar elements are independent
and identical symmetrically distributed random variables that are uniformly
bounded and satisfy
$\mathbb{E}\left(\left|1/\Delta_{kj|i}\right|\right)<\infty$, where
$j=1,\ldots,p$, and
$\boldsymbol{\Delta}_{k|i}^{-1}\equiv[\Delta_{k1|i}^{-1},\Delta_{k2|i}^{-1},\ldots,\Delta_{kp|i}^{-1}]^{T}$.

Here we generate an estimate of $\boldsymbol{H}(\boldsymbol{\theta})$ by
simultaneous perturbation (SPSA) methods:

$\hat{\boldsymbol{H}}_{k|i}=\frac{1}{2}\left[\frac{\delta\boldsymbol{g}_{k|i}}{2c}(\boldsymbol{\Delta}_{k|i}^{-1})^{T}+\left(\frac{\delta\boldsymbol{g}_{k|i}}{2c}(\boldsymbol{\Delta}_{k|i}^{-1})^{T}\right)^{T}\right],$
(1)

where
$\delta\boldsymbol{g}_{k|i}=\boldsymbol{g}(\boldsymbol{\theta}+c\boldsymbol{\Delta}_{k|i}|\boldsymbol{Z}_{\mathrm{pseudo}}(i))-\boldsymbol{g}(\boldsymbol{\theta}-c\boldsymbol{\Delta}_{k|i}|\boldsymbol{Z}_{\mathrm{pseudo}}(i))$,
and $\boldsymbol{g}(\cdot)$ is the (estimated if there is not enough
information) gradient function of $L$, and $c>0$ is a small constant.

The Monte Carlo-based estimate of $\boldsymbol{F}(\boldsymbol{\theta})$ in
[5], denoted $\bar{\boldsymbol{F}}_{M,N}(\boldsymbol{\theta})$, is

$\bar{\boldsymbol{F}}_{M,N}(\boldsymbol{\theta})\equiv-\frac{1}{N}\sum_{i=1}^{N}\frac{1}{M}\sum_{k=1}^{M}\hat{\boldsymbol{H}}_{k|i}.$
(2)

The first “inner” average calculates Hessian estimates at a given
$\boldsymbol{Z}_{\mathrm{pseudo}}(i)$ ($i=1,2,\ldots,N$) based on $M$ values
of $\hat{\boldsymbol{H}}_{k|i}$, and the second “outer” average sums these
sample means of Hessian estimates across $N$ values of pseudo data.

To get the most intuitive comparison, we set $M=N=1$ without loss of
generality since (2) indicates that Hessian estimates do not depend on $N$.
Additionally, while $\hat{\boldsymbol{H}}_{k|i}(k=1,\ldots,M)$ are estimated
based on the same pseudo data with fixed $i$, [10] shows $M=1$ is the optimal
choice if the $n$ vectors entering each $\boldsymbol{Z}_{\mathrm{pseudo}}(i)$
are mutually independent. Then the variance of the $jj$th entry of the
estimate $\bar{\boldsymbol{F}}_{M,N}(\boldsymbol{\theta})$ is

$\mathrm{var}\left\\{\left[\bar{\boldsymbol{F}}_{M,N}\right]_{jj}\right\\}=\mathrm{var}\left(\hat{H}_{jj}\right),$
(3)

where $\hat{H}_{jj}$ denotes the $jj$th entry of
$\hat{\boldsymbol{H}}=\hat{\boldsymbol{H}}_{1|1}$. We can use the calculation
at $M=N=1$ to readily extend to arbitrary matrix. Let $O_{(\cdot)}(c^{2})$
denote a random “big-$O$” term, where the subscript represents the relevant
randomness. For example, $O_{\boldsymbol{Z},\boldsymbol{\Delta}_{1}}(c^{2})$
denotes a random “big-$O$” term depends on
$\boldsymbol{Z}_{\mathrm{pseudo}}(1)$ and $\boldsymbol{\Delta}_{1}$ such that
$O_{\boldsymbol{Z},\boldsymbol{\Delta}_{1}}(c^{2})/c^{2}$ is bounded almost
surely as $c\rightarrow 0$. Then, by [8], the $jj$th entry of
$\hat{\boldsymbol{H}}$ is

$\hat{{H}}_{jj}={H}_{jj}+\sum_{l\neq
j}{H}_{jl}\frac{\Delta_{1l}}{\Delta_{1j}}+O_{\boldsymbol{Z},\boldsymbol{\Delta}_{1}}(c^{2}).$
(4)

Assume that

$\mathbb{E}(\Delta_{1l}/\Delta_{1j})=0,\;\mathrm{var}(\Delta_{1l}/\Delta_{1j})=v;$
(5)

where $v$ is a constant, then, given the independence of the
$\left\\{\Delta_{1j}\right\\}$, and the fact that
$\mathbb{E}({H}_{jl})=-{F}_{jl}$ ($j,l=1,2,\ldots,p$), we have:

$\displaystyle\mathrm{var}(\hat{{H}}_{jj})_{\text{basic}}=$
$\displaystyle\mathbb{E}\big{(}[{H}_{jj}+\sum_{l\neq
j}{H}_{jl}\frac{\Delta_{1l}}{\Delta_{1j}}+O_{\boldsymbol{Z},\boldsymbol{\Delta}_{1}}(c^{2})]^{2}\big{)}-\big{(}\mathbb{E}[{H}_{jj}+\sum_{l\neq
j}{H}_{jl}\frac{\Delta_{1l}}{\Delta_{1j}}+O_{\boldsymbol{Z},\boldsymbol{\Delta}_{1}}(c^{2})]\big{)}^{2}$
$\displaystyle=$
$\displaystyle\mathbb{E}\big{(}{H}_{jj}^{2}+2{H}_{jj}\sum_{l\neq
j}{H}_{jl}\frac{\Delta_{1l}}{\Delta_{1j}}+[\sum_{l\neq
j}{H}_{jl}\frac{\Delta_{1l}}{\Delta_{1j}}]^{2}+O_{\boldsymbol{Z},\boldsymbol{\Delta}_{1}}(c^{2})\big{)}-\big{(}\mathbb{E}\big{[}{H}_{jj}\big{]}+O(c^{2})\big{)}^{2}$
$\displaystyle=$ $\displaystyle\mathbb{E}({H}_{jj}^{2})+\sum_{l\neq
j}\mathbb{E}\bigg{[}{H}_{jl}\frac{\Delta_{1l}}{\Delta_{1j}}\bigg{]}^{2}-\big{(}\mathbb{E}\big{[}{H}_{jj}\big{]}\big{)}^{2}+O(c^{2})$
$\displaystyle=$ $\displaystyle\mathbb{E}({H}_{jj}^{2})+v\sum_{l\neq
j}\mathbb{E}({H}_{jl}^{2})-{F}_{jj}^{2}+O(c^{2})$ $\displaystyle=$
$\displaystyle\mathrm{var}({H}_{jj})+v\sum_{l\neq
j}\left[\mathrm{var}({H}_{jl})+{F}_{jl}^{2}\right]+O(c^{2})$ $\displaystyle=$
$\displaystyle\mathrm{var}({H}_{jj})+v\sum_{l\neq
j}\mathrm{var}({H}_{jl})+v\sum_{l\neq j}{F}_{jl}^{2}+O(c^{2}).$

###  2.2 Implementation with Independent Perturbations

The estimation of $\boldsymbol{F}(\boldsymbol{\theta})$ can be enhanced if the
$n$ vectors entering each $\boldsymbol{Z}_{\mathrm{pseudo}}(i)$ are mutually
independent. Note that the $n$ vectors do not have to be identically
distributed — just independent. By the independence, we have

$\boldsymbol{F}(\boldsymbol{\theta})=\sum_{t=1}^{n}\boldsymbol{F}^{(t)},$

where $\boldsymbol{F}^{(t)}$ denotes the FIM for each $t$.

In particular, the variance of the entries of each Hessian estimate
$\hat{\boldsymbol{H}}_{k|i}$ in the standard method can be reduced by
decomposing $\hat{\boldsymbol{H}}_{k|i}$ into a sum of $n$ independent
estimates $\hat{\boldsymbol{H}}_{k|i}^{(t)}$, each corresponding to one vector
in the pseudo data. To distinguish the way that $\boldsymbol{\Delta}_{k|i}$ is
generated in the basic method and the independent perturbation method,
$\boldsymbol{\Delta}_{I}$ is used here to represent the perturbation term
$\boldsymbol{\Delta}_{1|1}$ in the independent perturbation method, where
$\boldsymbol{\Delta}_{I}=\sum_{t=1}^{n}\boldsymbol{\Delta}_{I}^{(t)}$. Then a
separate perturbation vector $\boldsymbol{\Delta}_{I}^{(t)}$ can be applied to
each independent estimate, which results in variance reduction in the estimate
$\bar{\boldsymbol{F}}_{M,N}$. Hence, based on assumption (5), similar to the
above calculation process, we have:

$\displaystyle\mathrm{var}(\hat{{H}}_{jj})_{\text{indep}}=\mathrm{var}(\sum_{t=1}^{n}\hat{{H}}_{jj}^{(t)})=\sum_{t=1}^{n}\mathrm{var}(\hat{{H}}_{jj}^{(t)})$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{n}\\{\mathbb{E}([{H}_{jj}^{(t)}+\sum_{l\neq
j}{H}_{jl}^{(t)}\frac{\Delta_{Il}^{(t)}}{\Delta_{Ij}^{(t)}}+O_{\boldsymbol{Z},\boldsymbol{\Delta}_{I}}(c^{2})]^{2})-(\mathbb{E}[{H}_{jj}^{(t)}+\sum_{l\neq
j}{H}_{jl}^{(t)}\frac{\Delta_{Il}^{(t)}}{\Delta_{Ij}^{(t)}}+O_{\boldsymbol{Z},\boldsymbol{\Delta}_{I}}(c^{2})])^{2}\\}$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{n}\\{\mathbb{E}[({H}_{jj}^{(t)})^{2}]+v\sum_{l\neq
j}\mathbb{E}[({H}_{jl}^{(t)})^{2}]-({F}_{jj}^{(t)})^{2}+O(c^{2})\\}$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{n}\\{\mathrm{var}\left({H}_{jj}^{(t)}\right)+v\sum_{l\neq
j}[\mathrm{var}({H}_{jl}^{(t)})+({F}_{jl}^{(t)})^{2}]+O(c^{2})\\}$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{n}\\{\mathrm{var}({H}_{jj}^{(t)})+v\sum_{l\neq
j}\mathrm{var}({H}_{jl}^{(t)})\\}+v\sum_{t=1}^{n}[\sum_{l\neq
j}({F}_{jl}^{(t)})^{2}]+O(c^{2}).$

The variance of the $j$th diagonal entry of $\hat{\boldsymbol{H}}$ generated
from the standard method can be written as:

$\displaystyle\mathrm{var}(\hat{{H}}_{jj})_{\text{basic}}=$
$\displaystyle\mathrm{var}(\sum_{t=1}^{n}{H}_{jj}^{(t)})+v\sum_{l\neq
j}\mathrm{var}(\sum_{t=1}^{n}{H}_{jl}^{(t)})+v\sum_{l\neq
j}[\sum_{t=1}^{n}{F}_{jl}^{(t)}]^{2}+O(c^{2})$ $\displaystyle=$
$\displaystyle\sum_{t=1}^{n}\\{\mathrm{var}\left({H}_{jj}^{(t)}\right)+v\sum_{l\neq
j}\mathrm{var}({H}_{jl}^{(t)})\\}+v\sum_{l\neq
j}[\sum_{t=1}^{n}{F}_{jl}^{(t)}]^{2}+O(c^{2}).$

Therefore, the difference of the variance of the diagonal entries in the
estimate of $\boldsymbol{F}(\boldsymbol{\theta})$ between the standard method
and the independent perturbation method is

$\mathrm{var}(\hat{{H}}_{jj})_{\text{basic}}-\mathrm{var}(\hat{{H}}_{jj})_{\text{indep}}=2v\sum_{l\neq
j}\sum_{t_{1}<t_{2}}{F}_{jl}^{(t_{1})}{F}_{jl}^{(t_{2})}+O(c^{2}).$

###  2.3 Theoretical Comparison via Variance Reduction

Denoting $A=\sum_{t=1}^{n}\\{\mathrm{var}({H}_{jj}^{(t)})+v\sum_{l\neq
j}\mathrm{var}({H}_{jl}^{(t)})\\}$, the result above shows that the variance
difference between these two method is independent of $A$. That is,
$\mathrm{var}(\hat{{H}}_{jj})_{\text{indep}}$ and
$\mathrm{var}(\hat{{H}}_{jj})_{\text{basic}}$ contain the same components that
are made up of the variance of entries of the true Hessian matrix. Hence, we
should mainly consider the difference between $v\sum_{t=1}^{n}[\sum_{l\neq
j}({F}_{jl}^{(t)})^{2}]$ and $v\sum_{l\neq
j}[\sum_{t=1}^{n}{F}_{jl}^{(t)}]^{2}$. However, the result above cannot
directly reflect the degree of variance reduction. Let us consider the ratio
$\mathrm{var}(\hat{{H}}_{jj})_{\text{indep}}$ to
$\mathrm{var}(\hat{{H}}_{jj})_{\text{basic}}$, which can be written as:

$\frac{\text{var}(\hat{{H}}_{jj})_{\text{indep}}}{\text{var}(\hat{{H}}_{jj})_{\text{basic}}}=\frac{A+v\sum_{t=1}^{n}\left[\sum_{l\neq
j}({F}_{jl}^{(t)})^{2}\right]+O(c^{2})}{A+v\sum_{l\neq
j}[\sum_{t=1}^{n}{F}_{jl}^{(t)}]^{2}+O(c^{2})}.$ (6)

Usually $c$ is small (for example, 0.0001), so let us focus on the second
part: off-diagonal entries of the true Fisher matrix. With a fixed $j$, let us
for convenience denote entries ${F}_{jl}^{(t)}=f_{l}^{(t)}$, where $l\neq j$.
Under some assumptions, the ratio of the second part in (6),
$v\sum_{t=1}^{n}[\sum_{l\neq j}({F}_{jl}^{(t)})^{2}]/v\sum_{l\neq
j}[\sum_{t=1}^{n}{F}_{jl}^{(t)}]^{2}$, is $O(1/n)$. That is, assume that:

  1. 1.

There are at least 2 non-zero off-diagonal entries in each true Fisher matrix
$\boldsymbol{F}^{(t)}$.

  2. 2.

For every $j=1,2,\ldots,p$, the non-zero sequence
$\\{f_{l}^{(t)}\\}_{t=1}^{n}$ is limited within
$[f_{l}^{(\text{min})},f_{l}^{(\text{max})}]$ while fixing $l$, where
$f_{l}^{(\text{min})}$ is the lower bound of non-zero entries in
$\\{|f_{l}^{(t)}|\\}_{t=1}^{n}$, and $f_{l}^{(\text{max})}$ is the upper bound
of $\\{|f_{l}^{(t)}|\\}_{t=1}^{n}$. For convenience, denote
$r_{1_{j}}=\sum_{l\neq
j}^{p}(f_{l}^{(\text{max})}/{f_{l}^{(\text{min})}})^{2}$.

  3. 3.

Denote $f_{l}^{(+)}$ as the sum of all positive elements in the sequence, and
$f_{l}^{(-)}$ as the sum of all negative elements in the sequence, the ratio
is
$r_{2_{j}}=(f_{l}^{(+)}+|f_{l}^{(-)}|)^{2}/{(f_{l}^{(+)}-|f_{l}^{(-)}|)^{2}}$
such that $r_{1_{j}}r_{2_{j}}=O(1)$ as $n$ goes to infinity.

Based on the above conditions, we have:

Proposition 1: Suppose that conditions (1), (2), and (3) hold for all $j$.
Then the independent perturbations above reduce the variance of the diagonal
entries in the estimate of $\boldsymbol{F}(\boldsymbol{\theta})$ by
$O(1/n)+O(c^{2})$.

The details about the proof have been shown in [10], so we get:

$\displaystyle\frac{\text{var}(\hat{{H}}_{jj})_{\text{indep}}}{\text{var}(\hat{{H}}_{jj})_{\text{basic}}}\leq$
$\displaystyle\frac{A+\frac{\text{C}}{n}B_{\text{basic}}+O(c^{2})}{A+B_{\text{basic}}+O(c^{2})},$

where $C=r_{1_{j}}r_{2_{j}}$ is a constant,
$A=\sum_{t=1}^{n}\\{\mathrm{var}({H}_{jj}^{(t)})+v\sum_{l\neq
j}\mathrm{var}({H}_{jl}^{(t)})\\}=O(n)$, and $B_{\text{basic}}=v\sum_{l\neq
j}^{p}[\sum_{t=1}^{n}f_{l}^{(t)}]^{2}\leq n^{2}v\sum_{l\neq
j}^{p}[f_{l}^{(\text{max})}]^{2}=O(n^{2})$. Therefore,

$\frac{\text{var}\left(\hat{{H}}_{jj}\right)_{\text{indep}}}{\text{var}\left(\hat{{H}}_{jj}\right)_{\text{basic}}}=O\bigg{(}\frac{1}{n}\bigg{)}+O(c^{2}),$
(7)

where $O(c^{2})$ is often so small that it can be ignored (see Example 1 in
Section 3).

In addition, similar to Section 4.2 in [5], [10] shows that in the independent
case, $M>1$ is considered only when it is costly to generate new pseudodata
vectors. This is despite the fact that we mentioned earlier that the average
of $M$ estimated Hessian matrices, which is based on the same pseudodata
vector $\boldsymbol{Z}_{\text{pseudo}}$, is usually used to estimate a Fisher
information matrix.

That is, in general, if it is allowed to compute $C=MN$ Hessian matrix
estimates, we can maximize the accuracy of $\bar{\boldsymbol{F}}_{M,N}$ when
each estimate $\hat{\boldsymbol{H}}_{k|i}$ is generated based on a new
$\boldsymbol{Z}_{\text{pseudo}}$ (i.e., $M=1$).

Proposition 2: Assume that elements $\boldsymbol{\Delta}_{k|i}^{(t)}$ and
$\boldsymbol{Z}^{(t)}_{\text{pseudo}}(i)$ are mutually independent, where
$k=1,\ldots,M$, $i=1,\ldots,N$, and $t=1,\ldots,n$. For a fixed $C=MN$, the
variance of each entry in $\bar{\boldsymbol{F}}_{M,N}(\boldsymbol{\theta})$ is
minimized when $M=1$.

##  3 Numerical Examples

Here we introduce two examples, which are also used in [11]. The first one is
under the signal-plus noise setting. We calculate the variance ratio to verify
the theoretical result and analyze the distinction among the variance ratio
for different diagonal entries. The second one compares the performance for
the two methods when the pseudo data are mixture Gaussian distributed. At the
same time, we show the improvement of the enhanced method from time analysis
in both examples.

Both of the mixture distribution and the signal-plus-noise settings arise
regularly in the literature. Mixture problems are thoroughly reviewed in [12]
and [13], with applications in information theory in [14] and bivariate
quantile estimation in [15]. The signal-plus-noise problem with non-identical
noise distributions arises in practical problems where measurements are
collected with varying quality of information across the sample. Some
practical implications are discussed in [16] relative to the initial
conditions in a Kalman filter model and [17] in the context of outlier
analysis.

###  3.1 Multivariate Normal with Signal-Plus-Noise

Assume that the $\boldsymbol{z}_{i}$ are independently and normally
distributed with mean $\boldsymbol{\mu}$ and covariance
$\boldsymbol{\Sigma}+\boldsymbol{P}_{i}$ for all $i$, where $\boldsymbol{\mu}$
and $\boldsymbol{\Sigma}$ are unknown parameters to be estimated. It can be
interpreted that $\boldsymbol{z}_{i}$ is the observed values of the
$N(\boldsymbol{\mu},\boldsymbol{\Sigma})$ distributed signal with independent
$N(\boldsymbol{0},\boldsymbol{P}_{i})$ distributed noise. This setting has
been widely discussed; for example, [16] considers estimating the initial mean
vector and covariance matrix of a state-space model and [18] demonstrates
parameter estimation of random-coefficient linear models.

Let us consider a 3-dimensional case: $\text{dim}(\boldsymbol{z}_{i})=3$,
$n=30$, $\boldsymbol{\mu}=\boldsymbol{0}$, $\boldsymbol{\Sigma}$ is a matrix
with all 0.5’s except 2’s on the diagonal (which is a little different from
that in [10]), and
$\boldsymbol{P}_{i}=\sqrt{i}\boldsymbol{U}^{T}\boldsymbol{U}$, where
$\boldsymbol{U}$ is a $3\times 3$ matrix with entries generated from uniform
(0,1) distribution. $\boldsymbol{\theta}$ is the collection of unique elements
in $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$, so the FIM is $p\times p$
where $p=9$. Reference [16] gives the gradient function of log-likelihood and
the analytical form of the FIM in this setting.

In particular, assume that components of the perturbation vector
$\boldsymbol{\Delta}_{k|i}$ follow Bernoulli $\pm 1$ distribution. In both
methods, we estimate the Hessian matrix by using the gradient of the log-
likelihood function under $M=N=1$. The ratio values for estimated variance of
the diagonal entries of the estimated FIM from the standard method to that
from the independent perturbation method based on 2 million replicates are
shown in Table 1. Furthermore, we change $n$ from 30 to 200 to observe the
performance.

Table 1: Variance of $\bar{\boldsymbol{F}}_{M,N}(\boldsymbol{\theta})_{jj}$ Diagonal | Independent | Standard | Ratio
---|---|---|---
$j=1$ | 8.24 | 60.12 | 0.14
$j=2$ | 13.32 | 56.58 | 0.23
$j=3$ | 4.88 | 24.48 | 0.20
$j=4$ | 5.00 | 11.18 | 0.45
$j=5$ | 23.49 | 39.13 | 0.60
$j=6$ | 14.85 | 26.69 | 0.56
$j=7$ | 8.69 | 13.84 | 0.63
$j=8$ | 12.91 | 19.70 | 0.65
$j=9$ | 3.04 | 4.98 | 0.61

Table 1 indicates that there is an obvious reduction of the variance of
entries in the estimated FIM over the enhanced method compared with that of
the standard method for every $j$, which matches our theoretical conclusion in
Section 2.

Table 2: Time cost/$seconds$ | $n=30$ | $n=100$ | $n=200$
---|---|---|---
Independent Method | 8726 | 28914 | 57066
Standard Method | 5259 | 16978 | 33541
Ratio | 0.6027 | 0.5872 | 0.5878

On the other hand, researchers usually need to trade off accuracy and time
cost in experiments. Results with high accuracy generally take more time to
obtain, but the payoff may not be worthwhile if the computational cost is too
high. Therefore, under the same setting as above, Figures 1 and 2, and Table 2
show the performance of the standard method and the enhanced method as $n$ is
equal to 30, 100, and 200 based on the variance ratio
$\text{var}(\hat{H}_{jj})_{\text{indep}}/\text{var}(\hat{H}_{jj})_{\text{basic}}$
and time cost. To better compare the variance ratio with $O(1/n)$ as $n$
grows, we add the reference lines $3.5/n$ and $13/n$ in Figures 1 and 2,
separately. Figure 1 and 2 reveal that the change of variance ratio as $n$
grows is similar to the trend of $O(1/n)$. That is, the curves follow the
overall slope given by the reference line. But the time cost ratio in Table 2
is stable over the range of $n$ considered.

Figure 1: Variance ratio for $j=1,2,3$ Figure 2: Variance ratio for
$j=4,\ldots,9$

In addition, under specific regularity conditions [4], the $jj$th entry of
$\boldsymbol{F}(\boldsymbol{\theta})$ is the expectation of
$\partial^{2}L/\partial{\theta}_{j}^{2}$, where ${\theta}_{j}(j=1,\ldots,9)$
corresponds to three elements in $\boldsymbol{\mu}$ and six unique entries in
$\boldsymbol{\Sigma}$ respectively. We can see that the estimates of diagonal
entries of $\boldsymbol{F}(\boldsymbol{\theta})$ corresponding to
$\boldsymbol{\mu}$ part are more stable than those corresponding to the
$\boldsymbol{\Sigma}$ part. Furthermore, when $n=30$, the variance ratio
values in Figure 1 are significantly lower than 0.6027, which is the time cost
ratio in Table 2. Part of the variance ratio values in Figure 2, however, are
higher than 0.6027.

Note that compared with [10], we only change the diagonal values of the true
covariance matrix $\boldsymbol{\Sigma}$. But in [10], all variance ratio
values are much lower than 0.9586, which is the time cost ratio as $n=30$.
Here, the performance of variance ratios that corresponds to the
$\boldsymbol{\mu}$ part is much better than that corresponding to the
$\boldsymbol{\Sigma}$ part. That is, the variance ratio values corresponding
to the $\boldsymbol{\mu}$ part are lower and closer to $1/n$. This might be
because the true variance value for each dimension becomes large, which tend
to expand the range of generated values and affects the estimation for the
covariance part. Then the fluctuation of corresponding diagonal entries (for
the covariance part) in the estimated FIM are relatively obvious compared with
that for the mean part.

Although Table 2 indicates that it may not be worth the time (compare the time
cost ratio with part of the variance ratio values) to improve accuracy when
$n=30$, the results for $n=100,200$ from the independent method achieve a
significant reduction on variance with acceptable and reasonable time cost,
which means that our enhanced method is applicable in practice.

###  3.2 Mixture Gaussian Distribution

Denote the generated vector $\boldsymbol{Z}=[z_{1},z_{2},\ldots,z_{n}]^{T}$,
and assume that $z_{i}\;(i=1,\ldots,n)$ are independently and identically
distributed with probability density function:

$f(z,\boldsymbol{\theta})=\lambda\exp(-(z-\mu_{1})^{2}/(2\sigma_{1}^{2}))/\sqrt{2\pi\sigma_{1}^{2}}+(1-\lambda)\exp(-(z-\mu_{2})^{2}/(2\sigma_{2}^{2}))/\sqrt{2\pi\sigma_{2}^{2}},$

where
$\boldsymbol{\theta}=[\lambda,\mu_{1},\sigma_{1}^{2},\mu_{2},\sigma_{2}^{2}]^{T}$.
There are $5(5+1)/2=15$ unique entries in the FIM. Furthermore, we use the
average of the negative Hessian matrix based on the sample vector over
$10^{6}$ independent replicates to approximate (since [19] showed that the
Hessian matrix is attainable) the true FIM, which is not computable in this
example.

Here consider the same setting as that in [11], where
$\boldsymbol{\theta}=[0.2,0,4,1,9]^{T}$, $M=2$, $N=40000$, $c=0.0001$, and
elements of the perturbation $\boldsymbol{\Delta}_{k|i}$ follow Bernoulli $\pm
1$ distribution for all $k\;(k=1,\ldots,N)$ and $i\;(i=1,\ldots,M)$. We
estimate the true FIM $\boldsymbol{F}(\boldsymbol{\theta})$ given the gradient
of the log-likelihood function, and measure the performance of methods by
computing the sample mean of
$\|\boldsymbol{F}_{\text{est}}(\boldsymbol{\theta})-\boldsymbol{F}(\boldsymbol{\theta})\|/\|\boldsymbol{F}(\boldsymbol{\theta})\|$
based on 50 independent replicates, where $\|\bullet\|$ represents the
spectral norm. Table 3 gives the typical estimated FIM based on the standard
method and the independent method, respectively. By typical, we mean it is the
estimated FIM corresponding to the 25th value of
$\|\boldsymbol{F}_{\text{est}}(\boldsymbol{\theta})-\boldsymbol{F}(\boldsymbol{\theta})\|/\|\boldsymbol{F}(\boldsymbol{\theta})\|$
in descending order given 50 independent replicates. Note that the true FIM is
approximated by sample averages over $10^{6}$ independent replicates. Table 4
summarizes the performances (accuracy and time cost) between the standard
method and the independent approach.

Table 3: Simulation Results for Example 3.2 True FIM | $\left[\begin{array}[]{ccccc}796.8750&0&0&0&0\\\ 0&7.5000&1.5000&0&0\\\ 0&1.5000&3.1876&0&0\\\ 0&0&0&3.3333&-0.0741\\\ 0&0&0&-0.0741&0.1523\end{array}\right]$
---|---
Typical $\bar{\boldsymbol{F}}(\boldsymbol{\theta})_{\text{basic}}$ | $\left[\begin{array}[]{ccccc}796.8752&1.5414&1.8820&0.3499&-0.7470\\\ 1.5414&7.5004&1.4958&0.01550&-0.0075\\\ 1.8820&1.4958&3.1868&0.0157&-0.02251\\\ 0.3499&0.0155&0.0157&3.3332&-0.0863\\\ -0.7470&-0.0075&-0.02251&-0.0863&0.1521\end{array}\right]$
Typical $\bar{\boldsymbol{F}}(\boldsymbol{\theta})_{\text{indep}}$ | $\left[\begin{array}[]{ccccc}796.8752&0.2536&0.0278&0.4273&-0.0930\\\ 0.2536&7.4949&1.4922&-0.0043&-0.0009\\\ 0.0278&1.4922&3.1768&-0.0019&-0.0026\\\ 0.4273&-0.0043&-0.0019&3.3323&-0.0754\\\ -0.0930&-0.0009&-0.0026&-0.0754&0.1509\end{array}\right]$
Table 4: Comparison between Two Methods | Independent | Standard | Ratio
---|---|---|---
Accuracy ($\|\boldsymbol{F}_{\text{est}}(\boldsymbol{\theta})-\boldsymbol{F}(\boldsymbol{\theta})\|/\|\boldsymbol{F}(\boldsymbol{\theta})\|$) | 0.00063 | 0.00330 | 0.19
Time Cost/$seconds$ | 1760 | 1212 | 1.45

Table 4 shows the obvious advantage of the enhanced method with independent
perturbations while taking time cost into account. It is worth noting that our
accuracy of the basic approach is higher than that shown in [11] because for
each time, i.e., fixing $k$ and $i$, [11] only generates one sample following
the mixture Gaussian distribution, but we generate $n$ samples, which might
result in a more accurate estimation.

##  4 Conclusion and Future Work

It is often hard to estimate the FIM in real-world models because of high
dimension or nonlinear functions. This paper presents an enhanced Monte Carlo
method with independent perturbations to estimate the Fisher information
matrix in problems where it is not analytically available. The theoretical
analysis and numerical examples show that this approach improves the
estimation performance by reducing variance of elements compared to the basic
Monte Carlo approach.

Moreover, the above results indicate that the performance might change even if
we only change certain parameters in the model; for example, consider the time
cost ratios in Example 1 and [10]. However, as $n$ increases, the advantage of
the enhanced method is apparent with no doubt.

In future work, in addition to considering the performance of estimated FIM in
other practical problems (for example, Bayesian statistics [20] and the EM
algorithm [8]), it is important to account for the relative benefit and cost
of reducing variance, especially when $n$ is very large, where we cannot
ignore the cost of generating the additional perturbation vectors in the
independent perturbation method.

## References

  * [1] Ljung, L., System Identification—Theory for the User, Prentice-Hall, Upper Saddle River, N.J., 2nd edition, 1999, pp. 212–221.
  * [2] Spall, J. C., “Identification for Systems with Binary Subsystems,” IEEE Transactions on Automatic Control, 59(1), 2014, pp. 3–17.
  * [3] Park, S.H., Simeone, O., Sahin, O. and Shitz, S.S., “Fronthaul Compression for Cloud Radio Access Networks: Signal Processing Advances Inspired by Network Information Theory,” IEEE Signal Processing Magazine, 31(6), 2014, pp. 69–79.
  * [4] Wilks, S. S., Mathematical Statistics, Wiley, New York, 1962.
  * [5] Spall, J. C., “Monte Carlo Computation of the Fisher Information Matrix in Nonstandard Settings,” Journal of Computational and Graphical Statistics, 14, 2005, pp. 889–909.
  * [6] Louis, T. A., “Finding the Observed Information Matrix When Using the EM Algorithm,” Journal of the Royal Statistical Society. Series B (Methodological), 44, 1982, pp. 226–233.
  * [7] Spall, J. C., “Improved Methods for Monte Carlo Estimation of the Fisher Information Matrix,” Proceedings of the American Control Conference, Seattle, WA, 2008, pp. 2395–2400.
  * [8] Meng, L. and Spall, J. C., “Efficient Computation of the Fisher Information Matrix in the EM Algorithm,” Proceedings of the 51st Annual Conference on Information Sciences and Systems (CISS), 2017, pp. 1–6.
  * [9] Spall, J. C., “Adaptive Stochastic Approximation by the Simultaneous Perturbation Method,” IEEE Transactions on Automatic Control, 45(10), 2000, pp. 1839–1853.
  * [10] Wu, X. and Spall, J. C., “Improved Monte Carlo Estimation of the Fisher Information Matrix with Independent Perturbations,” Proceedings of the 55th Annual Conference on Information Sciences and Systems, online conference due to covid-19, 24–26 March 2021.
  * [11] Cao, X., “Demonstration of Enhanced Monte Carlo Computation of the Fisher Information for Complex Problems,” Proceedings of the American Control Conference, 2013, pp. 4003–4008.
  * [12] Redner, R. A. and Walker, H. F., “Mixture Densities, Maximum Likelihood, and the EM Algorithm,” SIAM Review, vol. 26(2), 1984, pp. 195–239.
  * [13] Titterington, D. M., Smith, A. F. M., and Makov, U. E., Statistical Analysis of Finite Mixture Distributions, Wiley, New York, 1985.
  * [14] Spall, J. C. and Hill, S. D., “Least-Informative Bayesian Prior Distributions for Finite Samples Based on Information Theory,” IEEE Transactions on Automatic Control, vol. 35(5), 1990, pp. 580–583.
  * [15] Spall, J. C. and Maryak, J. L., “A Feasible Bayesian Estimator of Quantiles for Projectile Accuracy from Non-i.i.d. Data,” Journal of the American Statistical Association, vol. 87, 1992, pp. 676–681.
  * [16] Shumway, R. H., Olsen, D. E., and Levy, L. J., “Estimation and Tests of Hypotheses for the Initial Mean and Covariance in the Kalman Filter Model,” Communications in Statistics—Theory and Methods, 1981, pp. 1625–1641.
  * [17] Spall, J. C. and Chin, D. C., “First-Order Data Sensitivity Measures with Applications to a Multivariate Signal-Plus-Noise Problem,” Computational Statistics and Data Analysis, vol. 9, 1990, pp. 297–307.
  * [18] Sun, F. K., “A Maximum Likelihood Algorithm for the Mean and Covariance of Nonidentically Distributed Observations,” IEEE Transactions on Automatic Control, 1982, pp. 245–247.
  * [19] Boldea, O. and Magnus, J. R., “Maximum Likelihood Estimation of the Multivariate Normal Mixture Model,” Journal of the American Statistical Association 104(488), 2009, pp. 1539–1549.
  * [20] Das, S., Spall, J.C. and Ghanem, R., “Efficient Monte Carlo Computation of Fisher Information Matrix Using Prior Information,” Computational Statistics & Data Analysis, 54(2), 2010, pp. 272–289.