text:
\section{Introduction}
\IEEEPARstart{F}{isher} information \cite{Fisher22} \cite{Fisher25} is a classical measure in information, estimation and decision theory with various applications. It can be used in order to predict the performance of efficient unbiased estimation algorithms in a compact way \cite{Rao45,Cram46,Kay93,Trees07} or for the construction of strong statistical tests \cite{Wald43,Engle83,Barbe11}. Further, it describes the covariance of maximum-likelihood estimates around the true parameter in the asymptotic regime where the amount of observations $N$ is large \cite{VanDerVaart} and finds application in Bayesian inference through the Bernstein-von Mises theorem \cite{LeCam}. While access to the Fisher information is useful in different situations, calculation of the information measure itself can turn out to be quite complex or even impossible. In particular this is the case when the evaluation of the log-likelihood function is intractable, when the integration of the score function is difficult or when the underlying statistical model is unknown. For example for a nonlinear squaring device with Gaussian input, the probability density function is in general of the noncentral chi-squared type and contains a modified Bessel function making the evaluation of Fisher information nontrivial. If one imagines a multivariate Gaussian distribution with $N$ dimensions which is element-wise hard-limited, the output likelihood function requires the orthant probabilities of the multivariate Gaussian model for which up to the present day no general closed-form exists \cite{Sinn11}. The required integration over the score function of the hard-limited Gaussian variable additionally requires to calculate a sum with $2^N$ components, resulting in prohibitively high computational complexity when $N$ is large. Considering measurement systems in practical scenarios, the exact analytical representation of the system model can rarely be deduced due to a high number of nonlinear and random effects like amplification, quantization, filtering, internal noise sources and rounding. The only possibility to access the Fisher information measure then is to approximate the parametric output distribution by means of empirical methods like histograms \cite{Barbe14}. For multivariate problems this is demanding and requires a high amount of memory.
\subsection{Related Work}
Early works dealing with the analysis of nonlinear systems are \cite{Wiener58} \cite{Brilliant58} \cite{Volterra59} and concentrate on the problem of representing the output directly in terms of the values of the input by using polynomials or series expansions with orthogonal functions. \cite{Oehlert92} considers the approximation of moments by polynomial representations of nonlinear functions. Another group of researchers, mainly concerned with communication issues, focuses on other important output characteristics like the power spectrum \cite{Rice44} \cite{Bennett48} \cite{Vleck66} or correlation functions \cite{Bussgang52} \cite{Price58} \cite{Baum69} of nonlinear devices while  \cite{Davenport53} discusses signal-to-noise ratio and \cite{Rowe82} distortion-to-signal power ratio at the output of nonlinearities. Classical attempts to derive the output distributions of nonlinear systems are found in \cite{Kac47} \cite{Meyer54}, while a recent result with applications in the field of Bayesian inference is \cite{Julier04}. A discussion on the properties of Fisher information is found in \cite{Stam59}, while \cite{Boek77} considers similar problems for a generalized version of Fisher information. The problem of constraint distributions providing minimum Fisher information is discussed in \cite{UhrKli95} \cite{Zivo97} \cite{Berch09}, whereas \cite{Stoica11} \cite{Park13} focus on the Gaussian assumption for additive systems under an estimation theoretic perspective.
\subsection{Motivation}
Anticipating that for future signal processing and communication systems analog complexity will be shifted into the digital domain in order to exploit the consequences of Moore's law and to obtain an energy- and hardware-efficient design, we have identified that a framework is required which allows to analyze the possible signal processing performance with challenging nonlinear stochastic models. Trying to circumvent exact calculation of the required Fisher information, we have recently shown \cite{Stein14} \cite{Stein15} that compact generic lower bounds for the information measure can be obtained. These approximations benefit from the fact that they exclusively rely on the first two \cite{Stein14} or four moments \cite{Stein15} of the output. Statistical moments are usually better tractable than the probability density or mass function itself and can be deduced by simple measurements \cite{Gallant96}. Through lower bounds, the derived expressions are guaranteed to be pessimistic. This makes them suitable as figure of merit for the problem of system design and optimization \cite{SteinWCL14}. While our last approach \cite{Stein15} on bounding the Fisher information measure turns out to provide tight results for various cases (Bernoulli, exponential, Gaussian or Poisson distribution), we have also identified an example where the information bound is loose (Laplace distribution). Looking at the properties of the example distributions, it becomes obvious that for the tight cases the sufficient statistics are the first two moments $z$ or $z^2$. In contrast, the zero-mean Laplace distribution has the sufficient statistic $|z|$. Such a statistic is not well captured by the first two moments, used in the construction of our bounding approach \cite{Stein15}. Might this be the reason why the approximation \cite{Stein15} fails to generate a tight bound under the Laplacian distribution?
\subsection{Outline}
Addressing this question, we pick the log-normal distribution with known scale parameter $\sigma$ and the Weibull distribution with known shape parameter $k$ as additional examples. Like the Laplace distribution, the log-normal has sufficient statistics which can not be represented by a finite number of raw moments, i.e., $\ln z$ and $\ln^2 z$. The Weibull distribution has the property that it's sufficient statistic is $z^k$. Therefore, the Fisher information bound \cite{Stein15} should be tight for the cases $k=1$ and $k=2$ and loose in any other configuration in order to confirm an existing connection between the construction of the information bound and the sufficient statistics.

A second interesting observation is the fact that all example distributions discussed in \cite{Stein15} belong to the class of the exponential family. Therefore, in order to obtain a better understanding on Fisher information and possible lower bounds, we adapt to this level of abstraction and study the properties of the Fisher information measure evaluated for distributions of the exponential family. Restricting to such a class of distributions, allows to provide an identity connecting the Fisher information measure and a weighted sum of the derivatives of the expected sufficient statistics. The weight of each derivative is the derivative of the associated natural parameter. The derived identity provides a guideline for the construction of strong Fisher information lower bounds for any kind of stochastic system. To this end, the original system is replaced by a counterpart in the exponential family which is equivalent with respect to a set of auxiliary statistics. We show that the Fisher information of this replacement is always dominated by the information measure of the original system such that an optimization of the replacement model will lead to a strong lower bound. The presented information bound has the advantage that instead of full characterization of the model likelihood, exclusively the expected values, the derivatives of the expected values and the covariance of the used auxiliary statistics are required in order to evaluate the estimation theoretic quality of the system.

We utilize this generic result in order to formulate a specific information bound involving the derivatives of the first $L$ raw moments. For the initial example of a log-normal and a Weibull distribution, we test the quality of this approximation of the Fisher information. By constructing a second bound which takes into consideration raw moments, the expected absolute value and the expected log-value, we show how to use the information bound in order to learn informative statistics of a parametric system with unknown output model and how to determine the minimum guaranteed interference capability of the model by calibrated measurements. In order to emphasize the practical impact of our discussion, we demonstrate this aspect with the Rapp model, which is popular for modeling the saturation effect of solid-state power amplifiers. Finally, based on the learned statistics and their optimized weights, we demonstrate that consistent estimators can be obtained from compressed observations of the system output which achieve a performance equivalent to the inverse of our pessimistic approximation for the Fisher information measure.  By reformulation of the proposed estimation algorithm, we reveal that the exponential replacement forms a conservative framework which allows to derive Hansen's famous estimation approach known as {generalized method of moments} \cite{Hansen82} by a maximum-likelihood argument. Special cases of the main result and further examples are found in our conference publication \cite{SteinI2MTC}. In \cite{SteinWSA16} we use the presented method in order to discuss the problem of parametric covariance estimation with $1$-bit hard-limiting.
\section{Approximations for the Fisher Information}
Consider a parametrized family of probability measures, characterized by a probability density or mass function $p(z;\theta)$, with random variable $z\in\mathcal{Z}$ and a deterministic parameter $\theta\in\Theta$. $\mathcal{Z}\subseteq\fieldR$ is the support of the random variable $z$ and $\Theta\subseteq\fieldR$ the parameter space of $\theta$. Throughout our discussion we assume that all integrands are absolutely integrable on $\mathcal{Z}$. All density or mass functions $p(z;\theta)$ exhibit regularity and are differentiable with respect to the parameter $\theta$.
\subsection{Fisher Information Bounds}
If the first
\begin{align}
\mu_1(\theta)&=\int_{\mathcal{Z}} z p_{z}(z;\theta) {\rm d}z
\end{align}
and the second central moment
\begin{align}
\mu_2(\theta)&=\int_{\mathcal{Z}} \big( z-\mu_1(\theta) \big)^2 p_{z}(z;\theta) {\rm d}z,
\end{align}
as well as the third
\begin{align}
\bar{\mu}_3(\theta)&=\int_{\mathcal{Z}} \bigg( \frac{ z-\mu_1(\theta) } {\sqrt{\mu_2(\theta)}} \bigg)^3 p_{z}(z;\theta) {\rm d}z
\end{align}
and the fourth central moment in normalized form
\begin{align}
\bar{\mu}_4(\theta)&=\int_{\mathcal{Z}} \bigg( \frac{ z-\mu_1(\theta) } {\sqrt{\mu_2(\theta)}} \bigg)^4 p_{z}(z;\theta) {\rm d}z
\end{align}
are at hand, it can be shown that the Fisher information
\begin{align}
F(\theta)&=\int_{\mathcal{Z}} \bigg(\frac{\partial \ln{p(z;\theta)}}{\partial\theta}\bigg)^2 p(z;\theta)  {\rm d}z
\label{measure:fisher}
\end{align}
is in general bounded from below by the expression \cite{Stein15}
\begin{align}
F(\theta) &\geq S(\theta)\notag\\
&=\frac{1}{{\mu_2(\theta)}} \frac{\Big(\frac{\partial \mu_1(\theta)}{\partial\theta} + \frac{b(\theta)}{\sqrt{\mu_2(\theta)}} \frac{\partial \mu_2(\theta)}{\partial\theta} \Big)^2}{1+ 2 b(\theta) \bar{\mu}_3(\theta) + {b}^2(\theta) (\bar{\mu}_4(\theta)-1)}.
\label{bound:fisher:two}
\end{align}
The weighting factor $b(\theta)$ which provides the maximum value on the right-hand side of \eqref{bound:fisher:two} is
\begin{align}
b^{\star}(\theta)&=\frac{\frac{\partial \mu_1(\theta)}{\partial\theta} \sqrt{\mu_2(\theta)}\bar{\mu}_3 (\theta) - \frac{\partial \mu_2(\theta)}{\partial\theta}}
{ \frac{\partial \mu_2(\theta)}{\partial\theta} \bar{\mu}_3 (\theta) -\frac{\partial \mu_1(\theta)}{\partial\theta} \sqrt{\mu_2(\theta)} (\bar{\mu}_4 (\theta)-1)}.
\end{align}
Setting the weighting factor $b(\theta)$ to zero, gives a simple special case of the information bound
\begin{align}\label{bound:fisher:one}
F(\theta) &\geq \frac{1}{{\mu_2(\theta)}} {\Big(\frac{\partial \mu_1(\theta)}{\partial\theta} \Big)^2}.
\end{align}
The inequality \eqref{bound:fisher:one} was the starting point of our discussion on lower bounds for the Fisher information in \cite{Stein14}. In our subsequent work \cite{Stein15} it was shown that \eqref{bound:fisher:two} produces tight results for Bernoulli, exponential, Gaussian and Poisson distributions and is loose for the Laplacian distribution.
\subsection{Example - Log-normal Distribution}
In order to provide further cases where the bound \eqref{bound:fisher:two} attains loose results, we start with the example of a log-normal distribution. The log-normal distribution with known scale parameter $\sigma>0$ is characterized by the probability density
\begin{align}
p(z;\theta)=\frac{1}{\sqrt{2\pi \sigma^2} z} \expe{-\frac{1}{2\sigma^2}{(\log{z}-\theta)^2}}
\label{model:lognormal}
\end{align}
with $z>0$. The log-likelihood function is given by
\begin{align}
\ln p(z;\theta)=-\frac{1}{2}\ln{2\pi \sigma^2} -\ln{z} -\frac{1}{2\sigma^2}{(\log{z}-\theta)^2}.
\end{align}
The score function is
\begin{align}
\frac{\partial \ln p(z;\theta)}{\partial \theta}= \frac{1}{\sigma^2}(\log{z}-\theta),
\end{align}
with it's derivative
\begin{align}
\frac{\partial^2 \ln p(z;\theta)}{\partial \theta^2}= -\frac{1}{\sigma^2}.
\end{align}
Therefore, the Fisher information measure with respect to the parameter $\theta$ is given by
\begin{align}\label{fisher:information:lognormal}
F(\theta)&=\exdi{z;\theta}{ \bigg(\frac{ \partial \ln p(z;\theta)}{\partial \theta}\bigg)^2 }\notag\\
&=-\exdi{z;\theta}{\frac{ \partial^2 \ln p(z;\theta)}{\partial \theta^2}}=\frac{1}{\sigma^2}.
\end{align}
For the evaluation of the moments we use the fact that the $l$-th raw moment of the log-normal distribution follows from
\begin{align}
\exdi{z;\theta}{z^l}&=\tilde{\mu}_l(\theta)\notag\\
&=\expe{l\theta+\frac{1}{2} l^2 \sigma^2}.
\end{align}
Accordingly, the first moment and its derivative are given by
\begin{align}
\mu_1(\theta)&=\exdi{z;\theta}{z}\notag\\
&=\expe{\theta+\frac{1}{2}\sigma^2},\\
\frac{\partial \mu_1(\theta)}{\partial \theta}&=\expe{\theta+\frac{1}{2}\sigma^2}.
\end{align}
The second central moment and its derivative are
\begin{align}
\mu_2(\theta)&=\exdi{z;\theta}{\big(z-\mu_1(\theta)\big)^2} \notag\\
&=\tilde{\mu}_2(\theta)-\mu_1^2(\theta),\notag\\
&= \expe{2\theta+\sigma^2}(\expe{\sigma^2}-1)\\
\frac{\partial \mu_2(\theta)}{\partial \theta}&= 2\expe{2\theta+\sigma^2}(\expe{\sigma^2}-1).
\end{align}
The third normalized central moment is obtained by
\begin{align}
\bar{\mu}_3(\theta)&=\exdi{z;\theta}{\bigg(\frac{z-\mu_1(\theta)}{\sqrt{\mu_2(\theta)}}\bigg)^3}\notag\\
&=\frac{ \tilde{\mu}_3 -3\mu_1(\theta)\tilde{\mu}_2(\theta)+2\mu_1^3(\theta) }{\mu_2^{\frac{3}{2}}(\theta)}\notag\\
&={  \sqrt{(\expe{\sigma^2}-1)} (\expe{\sigma^2}+2)  }
\end{align}
and the fourth normalized central moment is found to be
\begin{align}
\bar{\mu}_4(\theta)&=\exdi{z;\theta}{\bigg(\frac{z-\mu_1(\theta)}{\sqrt{\mu_2(\theta)}}\bigg)^4}\notag\\
&=\frac{\tilde{\mu}_4 - 4 \mu_1(\theta)\tilde{\mu}_3(\theta)+6\mu_1^2(\theta)\tilde{\mu}_2(\theta)-3\mu_1^4(\theta) }{\mu_2^2(\theta)}\notag\\
&=3\expe{2\sigma^2}+2\expe{3\sigma^2}+\expe{4\sigma^2}-3.
\end{align}
For the assessment of the quality we define the ratio between the approximation \eqref{bound:fisher:two} and the exact information measure \eqref{fisher:information:lognormal}
\begin{align}\label{info:loss:lbs}
\chi(\theta)=\frac{S(\theta)}{F(\theta)}.
\end{align}
As the gap \eqref{info:loss:lbs} is independent of $\theta$, in Fig. \ref{LogNormal_Loss}, we depict $\chi(\theta)$ for different values of the known scale parameter $\sigma$. It can be observed that for $\sigma>1$ the difference between $S(\theta)$ and the exact Fisher information $F(\theta)$ becomes large.
\pgfplotsset{legend style={rounded corners=2pt,nodes=right}}
\begin{figure}[!htbp]
\begin{tikzpicture}[scale=1.0]

  	\begin{axis}[ylabel=$10 \log{\chi(\theta)}$,
  			xlabel=$\sigma$,
			grid,
			ymin=-12,
			ymax=0,
			xmin=0.2,
			xmax=2,
			legend pos=south east]

    			\addplot[black, style=solid, line width=0.75pt,smooth, every mark/.append style={solid}, mark=diamond*, mark repeat=50] table[x index=0, y index=1]{LogNormal_Loss.txt};
	\end{axis}
\end{tikzpicture}
\caption{Log-Normal Distribution - Information Loss}
\label{LogNormal_Loss}
\end{figure}
\subsection{Example - Weibull Distribution}
As a second example, we study the quality of the Fisher information bound $S(\theta)$ for the case of a Weibull distribution. The Weibull distribution with known shape parameter $k$ is given by the probability density function
\begin{align}
p(z;\theta)= \frac{k}{\theta} \Big( \frac{z}{\theta} \Big)^{k-1} \expe{- (\frac{z}{\theta})^k }
\label{model:weilbull}
\end{align}
with $z,\theta>0$. The log-likelihood of the distribution is
\begin{align}
\ln p(z;\theta)&= \ln{k}-\ln{\theta} + (k-1)\ln{z}\notag\\
&- (k-1)\ln \theta - \Big( \frac{z}{\theta} \Big)^{k}
\end{align}
and the score function is given by
\begin{align}
\frac{ \partial \ln p(z;\theta)}{\partial \theta}
&= \frac{k}{\theta} \bigg( \Big( \frac{z}{\theta} \Big)^{k} -1\bigg).
\end{align}
The derivative of the score function has the form
\begin{align}
\frac{ \partial^2 \ln p(z;\theta)}{\partial \theta^2}
&= -\frac{k (k+1)}{\theta^2} \Big( \frac{z}{\theta} \Big)^{k} +\frac{k}{\theta^2}.
\end{align}
Consequently, the Fisher information measure is given by
\begin{align}
F(\theta)
&=\frac{k}{\theta^2}\bigg(\frac{(k+1)}{\theta^{k}} \exdi{z;\theta}{z^k}-1\bigg) \notag\\
&=\frac{k}{\theta^2}\Big((k+1) \Gamma(2)-1\Big)=\bigg(\frac{k}{\theta}\bigg)^2,
\end{align}
where we used the property that the $l$-th raw moment of the Weibull distribution is
\begin{align}
\tilde{\mu}_l&=\exdi{z;\theta}{z^l}\notag\\
&=\theta^l\Gamma_l
\end{align}
with the shorthand notational convention
\begin{align}
\Gamma_l=\Gamma\Big(1+\frac{l}{k}\Big)
\end{align}
for the Gamma function
\begin{align}
\Gamma(x)=\int_{0}^{\infty} w^{x-1} \exp{-w} {\rm d}w.
\end{align}
For the information bound $S(\theta)$ we require the first moment
\begin{align}
\mu_1(\theta)&=\exdi{z;\theta}{z}\notag\\
&=\theta\Gamma_1,
\end{align}
it's derivative
\begin{align}
\frac{\partial \mu_1(\theta)}{\partial \theta}&=\Gamma_1,
\end{align}
the second central moment
\begin{align}
\mu_2(\theta)&=\exdi{z;\theta}{\big(z-\mu_1(\theta)\big)^2} \notag\\
&=\tilde{\mu}_2(\theta)-\mu_1^2(\theta)\notag\\
&= \theta^2 \big( \Gamma_2 - \Gamma_1^2 \big),
\end{align}
it's derivative
\begin{align}
\frac{\partial \mu_2(\theta)}{\partial \theta}&= 2\theta \big( \Gamma_2 - \Gamma_1^2 \big),
\end{align}
the third normalized central moment
\begin{align}
\bar{\mu}_3(\theta)&=\exdi{z;\theta}{\bigg(\frac{z-\mu_1(\theta)}{\sqrt{\mu_2(\theta)}}\bigg)^3}\notag\\
&=\frac{ \tilde{\mu}_3 -3\mu_1(\theta)\tilde{\mu}_2(\theta)+2\mu_1^3(\theta) }{\mu_2^{\frac{3}{2}}(\theta)}\notag\\
&=\frac{  \Gamma_3 -3 \Gamma_1\Gamma_2+2\Gamma_1^3 }{  \big( \Gamma_2 - \Gamma_1^2 \big)^{\frac{3}{2}} }
\end{align}
and the fourth normalized central moment
\begin{align}
\bar{\mu}_4(\theta)&=\exdi{z;\theta}{\bigg(\frac{z-\mu_1(\theta)}{\sqrt{\mu_2(\theta)}}\bigg)^4}\notag\\
&=\frac{\tilde{\mu}_4 - 4 \mu_1(\theta)\tilde{\mu}_3(\theta)+6\mu_1^2(\theta)\tilde{\mu}_2(\theta)-3\mu_1^4(\theta) }{\mu_2^2(\theta)}\notag\\
&=\frac{\Gamma_4 -    4\Gamma_1\Gamma_3 +6 \Gamma_1^2\Gamma_2-3\Gamma_1^4   }{ \big( \Gamma_2 - \Gamma_1^2 \big)^2}.
\end{align}
In Fig. \ref{Weibull_Loss} we plot the information loss $\chi(\theta)$ for different shape values $k=1,\ldots,5$. Note, that also in this example the loss $\chi(\theta)$ is independent of the parameter $\theta$.
\begin{figure}[!htbp]
\begin{tikzpicture}[scale=1.0]

  	\begin{axis}[ylabel=$10 \log{\chi(\theta)}$,
  			xlabel=$k$,
			grid,
			ymin=-0.1,
			ymax=0,
			ytick={-0.10,-0.08,-0.06,-0.04,-0.02,0.00},
			xmin=1,
			xmax=5,
			tick label style={/pgf/number format/fixed},
			legend pos=south east]

    			\addplot[black, style=solid, line width=0.75pt,smooth, every mark/.append style={solid}, mark=diamond*, mark repeat=50] table[x index=0, y index=1]{Weibull_Loss.txt};

\end{axis}
\end{tikzpicture}
\caption{Weibull Example - Information Loss}
\label{Weibull_Loss}
\end{figure}
It can be observed that the information bound $S(\theta)$ is tight for the cases where the shape of the distribution is $k=1$ and $k=2$ while the quality degrades significantly for the cases $k>2$.
\subsection{Connection to the Sufficient Statistics}
The Weibull example is of special interest for the analysis of the quality of the approximation $S(\theta)$ as it contains cases where the bound is tight and cases where it is not. In \cite{Stein15} it has been verified that $S(\theta)$ is tight for the Bernoulli, exponential, Gaussian and Poisson distributions. These distributions have in common that their sufficient statistics are $z$ and $z^2$. Interestingly, for the Weibull distribution the sufficient statistic is given by $z^k$ and the information bound $S(\theta)$ is only tight for $k=1$ and $k=2$.
Additionally, the approximation of the Fisher information $S(\theta)$ is loose for the log-normal distribution and the Laplacian distribution \cite{Stein15}, both cases where the sufficient statistics are distinct from $z$ or $z^2$. This provides an indication for the existence of a connection between the sufficient statistics and the quality of Fisher information bounds. Note, that in \cite{Stein15} the information bound $S(\theta)$ was constructed through the Cauchy-Schwarz inequality
\begin{align}\label{ineq:cs}
\int_{\mathcal{Z}} {f^2(z;\theta)} p(z;\theta){\rm d}z&\geq \frac{\Big(\int_{\mathcal{Z}} f(z;\theta)g(z;\theta)p(z;\theta) {\rm d}z\Big)^2}{\int_{\mathcal{Z}}g^2(z;\theta)p(z;\theta) {\rm d}z}
\end{align}
by using the functions
\begin{align}
f(z;\theta)&=\frac{\partial \ln{p(z;\theta)}}{\partial\theta}
\end{align}
and
\begin{align}
g(z;\theta) = \bigg(\frac{z-\mu_1(\theta)}{\sqrt{\mu_2(\theta)}}\bigg) +b(\theta) \bigg(\frac{z-\mu_1(\theta)}{\sqrt{\mu_2(\theta)}}\bigg)^2 -b(\theta).
\label{func:bound:g}
\end{align}
The latter contains the statistics $z$ and $z^2$ in normalized central form. This corresponds to the sufficient statistics of the distributions where (\ref{bound:fisher:two}) provides tight results.
\section{Fisher Information and the\\Exponential Family}
In order to strengthen the insights on a possible connection of the sufficient statistics and strong Fisher information bounds, we require an approach that allows to investigate and cover a brought class of distributions. All example distributions discussed so far belong to the univariate exponential family $z\in\fieldR$ with a single parameter $\theta\in\fieldR$. Therefore, in the following we study the properties of the Fisher information measure for distributions from this particular class. In order to provide the results in a generalized form we focus in the following on the multivariate case $\ve{z}\in\fieldR^{N}$ and a multi-dimensional parameter where $\ve{\theta}\in\fieldR^{M}$.
\subsection{Multivariate Exponential Family}
The multivariate exponential family with a parameter $\ve{\theta}\in\fieldR^{M}$ is the set of probability density or mass functions, which can be factorized
\begin{align}\label{exp:family:vec}
p(\ve{z};\ve{\theta})=\exp{\sum_{l=1}^{L} w_l(\ve{\theta}) t_l(\ve{z}) - \lambda(\ve{\theta})+\kappa(\ve{z})}.
\end{align}
$w_l(\ve{\theta})\in\fieldR$ is the $l$-th natural parameter, $t_l(\ve{z})\in\fieldR$ is the associated sufficient statistic, $\lambda(\ve{\theta})\in\fieldR$ is the log-normalizer and $\kappa(\ve{z})\in\fieldR$ is the so-called carrier measure. The log-likelihood function of the exponential family is given by
\begin{align}
\ln p(\ve{z};\ve{\theta}) &= \sum_{l=1}^{L} w_l(\ve{\theta}) t_l(\ve{z}) - \lambda(\ve{\theta})+\kappa(\ve{z}),
\end{align}
while the score function attains the structure
\begin{align}\label{exp:score:vec}
\frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} &= \sum_{l=1}^{L} \frac{\partial w_l(\ve{\theta})}{\partial \ve{\theta}} t_l(\ve{z}) - \frac{\partial\lambda(\ve{\theta})}{\partial \ve{\theta}}.
\end{align}
Note, that we strictly follow the notational convention
\begin{align}
\left[ \frac{\partial \ve{x}(\ve{y})}{\partial \ve{y}} \right]_{ij} = \frac{\partial {x}_i(\ve{y})}{\partial y_j}.
\end{align}
An essential property of the score function \eqref{exp:score:vec} is that it's expected value vanishes, i.e.,
\begin{align}\label{exp:score:vec:zero}
\exdi{\ve{z};\ve{\theta}}{\frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}}}=\ve{0}^{\T}.
\end{align}
Therefore, with \eqref{exp:score:vec}
\begin{align}\label{identity: log:normalizer}
\sum_{l=1}^{L} \frac{\partial w_l(\ve{\theta})}{\partial \ve{\theta}} \exdi{\ve{z};\ve{\theta}}{t_l(\ve{z})} = \frac{\partial\lambda(\ve{\theta})}{\partial \ve{\theta}}.
\end{align}
\subsection{Fisher Information for the Exponential Family}
The Fisher information measure for the case of a multi-dimensional parameter $\ve{\theta}\in\fieldR^{M}$ has matrix form $\ve{F}(\ve{\theta})\in\fieldR^{M \times M}$ and is defined by
\begin{align}\label{def:fisher:matrix}
\ve{F}(\ve{\theta})=\int_{\mathcal{Z}} \bigg( \frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} \bigg)^{\T} \bigg( \frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} \bigg) p(\ve{z};\ve{\theta}) {\rm d} \ve{z}.
\end{align}
Note that the Fisher information measure \eqref{def:fisher:matrix} exists for all probability laws of the form \eqref{exp:family:vec}.
\begin{identity}[Fisher Identity]
For any parametric probability density or mass function $p(\ve{z};\ve{\theta})$ belonging to the exponential family \eqref{exp:family:vec}, the Fisher information matrix is given by
\begin{align}\label{theorem:fish:identity:vec}
\ve{F}(\ve{\theta})&=\sum_{l=1}^{L}   \bigg( \frac{\partial \exdi{\ve{z};\ve{\theta}}{ t_l(\ve{z}) } }{\partial \ve{\theta}}  \bigg)^{\T}   \frac{\partial w_l(\ve{\theta})}{\partial \ve{\theta}}.
\end{align}
\end{identity}
\begin{IEEEproof}
With the definition of the Fisher information \eqref{def:fisher:matrix} we obtain
\begin{align}
\ve{F}(\ve{\theta})&=\int_{\mathcal{Z}} \bigg( \frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} \bigg)^{\T} \bigg( \frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} \bigg) p(\ve{z};\ve{\theta}) {\rm d} \ve{z}\notag\\
&=\int_{\mathcal{Z}} \bigg( \frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} \bigg)^{\T} \bigg( \sum_{l=1}^{L} \frac{\partial w_l(\ve{\theta})}{\partial \ve{\theta}} t_l(\ve{z})  \bigg) p(\ve{z};\ve{\theta}) {\rm d} \ve{z}\notag\\
&-\int_{\mathcal{Z}} \bigg( \frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} \bigg)^{\T} \frac{\partial\lambda(\ve{\theta})}{\partial \ve{\theta}} p(\ve{z};\ve{\theta}) {\rm d} \ve{z},
\end{align}
where we have used \eqref{exp:score:vec} in the second step to substitute one of the involved score functions. Using \eqref{exp:score:vec:zero}, we obtain
\begin{align}
\ve{F}(\ve{\theta})&=  \sum_{l=1}^{L}  \bigg( \int_{\mathcal{Z}} \frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} t_l(\ve{z}) p(\ve{z};\ve{\theta}) {\rm d} \ve{z} \bigg)^{\T} \frac{\partial w_l(\ve{\theta})}{\partial \ve{\theta}} \notag\\
&-\Bigg( \int_{\mathcal{Z}}  \bigg( \frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} \bigg)^{\T}  p(\ve{z};\ve{\theta}) {\rm d} \ve{z} \Bigg) \frac{\partial\lambda(\ve{\theta})}{\partial \ve{\theta}}\notag\\
&=\sum_{l=1}^{L}   \bigg(\int_{\mathcal{Z}} \frac{\partial  p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} t_l(\ve{z})  {\rm d} \ve{z} \bigg)^{\T}   \frac{\partial w_l(\ve{\theta})}{\partial \ve{\theta}}\notag\\
&=\sum_{l=1}^{L}   \bigg( \frac{\partial \exdi{\ve{z};\ve{\theta}}{ t_l(\ve{z}) } }{\partial \ve{\theta}}  \bigg)^{\T}   \frac{\partial w_l(\ve{\theta})}{\partial \ve{\theta}}.
\end{align}
\end{IEEEproof}
Defining the vector of sufficient statistics
\begin{align}
\ve{t}(\ve{z})=\begin{bmatrix} t_1(\ve{z}) &t_2(\ve{z}) &\ldots &t_L(\ve{z})\end{bmatrix}^{\T},
\end{align}
it's expected value
\begin{align}
\ve{\mu}_{\ve{t}}(\ve{\theta}) = \exdi{\ve{z};\ve{\theta}}{\ve{t}(\ve{z})},
\end{align}
and the vector of natural parameters
\begin{align}
\ve{w}(\ve{\theta})=\begin{bmatrix} w_1(\ve{\theta}) &w_2(\ve{\theta}) &\ldots &w_L(\ve{\theta})\end{bmatrix}^{\T},
\end{align}
we can reformulate the identity \eqref{theorem:fish:identity:vec} in a compact form
\begin{align}\label{fish:identity:vec:compact}
\ve{F}(\ve{\theta})&= \bigg( \frac{\partial \ve{\mu}_{\ve{t}}(\ve{\theta})}{\partial \ve{\theta}} \bigg)^{\T} \frac{\partial \ve{w}(\ve{\theta})}{\partial \ve{\theta}}.
\end{align}

\iffalse

\begin{align}\label{def:exp:family}
p(z;\theta)=\exp{\sum_{l=1}^{L} w_l(\theta) t_l(z) - \lambda(\theta)+\kappa(z)}
\end{align}

\begin{align}
\ln p(z;\theta) &= \sum_{l=1}^{L} w_l(\theta) t_l(z) - \lambda(\theta)+\kappa(z),
\end{align}
while the score function attains the structure
\begin{align}\label{exp:score}
\frac{\partial \ln p(z;\theta)}{\partial \theta} &= \sum_{l=1}^{L} \frac{\partial w_l(\theta)}{\partial \theta} t_l(z) - \frac{\partial\lambda(\theta)}{\partial \theta}.
\end{align}

\begin{identity}[Fisher Identity]
For any probability density or mass function $p(z;\theta)$ belonging to an exponential family \eqref{def:exp:family} the Fisher information is given by
\begin{align}\label{exp:fisher:ident}
F(\theta)=\sum_{l=1}^{L} \frac{\partial w_l(\theta)}{\partial \theta} \frac{\partial \exdi{z;\theta}{t_l(z)}}{\partial\theta}.
\end{align}
\end{identity}
\begin{IEEEproof}
By the definition of the Fisher information and the score function for exponential families \eqref{exp:score} we have
\begin{align}
&F(\theta)=\exdi{z;\theta}{ \bigg(\frac{\partial \ln{p(z;\theta)}}{\partial\theta}\bigg)^2 }\notag\\
&=\int_{\mathcal{Z}} \bigg(\frac{\partial \ln{p(z;\theta)}}{\partial\theta}\bigg) \bigg(\frac{\partial \ln{p(z;\theta)}}{\partial\theta}\bigg) p(z;\theta)  {\rm d}z\notag\\
&=\int_{\mathcal{Z}} \bigg(\frac{\partial \ln{p(z;\theta)}}{\partial\theta}\bigg) \bigg( \sum_{l=1}^{L} \frac{\partial w_l(\theta)}{\partial \theta} t_l(z)- \frac{\partial\lambda(\theta)}{\partial \theta} \bigg) p(z;\theta)  {\rm d}z.
\end{align}
By a straightforward manipulation, we obtain
\begin{align}
F(\theta)&=\sum_{l=1}^{L} \frac{\partial w_l(\theta)}{\partial \theta}  \int_{\mathcal{Z}} \frac{\partial \ln{p(z;\theta)}}{\partial\theta} t_l(z) p(z;\theta)  {\rm d}z\notag\\
&-\frac{\partial\lambda(\theta)}{\partial \theta} \int_{\mathcal{Z}} \frac{\partial \ln{p(z;\theta)}}{\partial\theta} p(z;\theta)  {\rm d}z.
\end{align}
Using the fact that with the regularity of $p(z;\theta)$ we can interchange differentiation and integration
\begin{align}\label{exp:log:zero}
\int_{\mathcal{Z}} \frac{\partial {\ln p(z;\theta)}}{\partial\theta} p(z;\theta) {\rm d}z&=\int_{\mathcal{Z}} \frac{\partial {p(z;\theta)}}{\partial\theta}  {\rm d}z\notag\\
&= \frac{\partial }{\partial\theta}  \int_{\mathcal{Z}} {p(z;\theta)} {\rm d}z\notag\\
&=0,
\end{align}
we finally obtain
\begin{align}
F(\theta)&=\sum_{l=1}^{L} \frac{\partial w_l(\theta)}{\partial \theta}  \int_{\mathcal{Z}} \frac{\partial {p(z;\theta)}}{\partial\theta} t_l(z) {\rm d}z\notag\\
&=\sum_{l=1}^{L} \frac{\partial w_l(\theta)}{\partial \theta} \frac{\partial }{\partial\theta} \int_{\mathcal{Z}}  t_l(z) p(z;\theta) {\rm d}z\notag\\
&=\sum_{l=1}^{L} \frac{\partial w_l(\theta)}{\partial \theta} \frac{\partial \exdi{z;\theta}{t_l(z)}}{\partial\theta}.
\end{align}
\end{IEEEproof}
The identity \eqref{exp:fisher:ident} shows that the Fisher information measure for an exponential family can be derived from a weighted sum of the derivatives of the $L$ expected sufficient statistics $\frac{\partial \exdi{z;\theta}{t_l(z)}}{\partial\theta}$. The weighting factors are the derivatives of the associated $L$ natural parameters $\frac{\partial w_l(\theta)}{\partial \theta}$ with respect to the system parameter $\theta$.
\fi
\subsection{Example - Parametric Gaussian Model}
Consider as an example the univariate Gaussian distribution
\begin{align}\label{parametric:gauss}
p(z;\theta)=\frac{1}{\sqrt{2\pi{\sigma}^2(\theta)}} {\rm e}^{- \frac{(z-{\mu}(\theta))^2}{2{\sigma}^2(\theta)} }
\end{align}
with $z\in\fieldR$ and parameter $\theta\in\fieldR$. The natural parameters are
\begin{align}
w_1(\theta)&=\frac{{\mu}(\theta)}{{\sigma}^2(\theta)},\\
w_2(\theta)&=-\frac{1}{2{\sigma}^2(\theta)}
\end{align}
and the two corresponding sufficient statistics are
\begin{align}
t_1(z)&=z,\\
t_2(z)&=z^2.
\end{align}
Therefore, the expectations of the sufficient statistics are
\begin{align}
\exdi{z;\theta}{t_1(z)}&={\mu}(\theta),\\
\exdi{z;\theta}{t_2(z)}&={\sigma}^2(\theta)+{\mu}^2(\theta).
\end{align}
With the derivatives of the natural parameters
\begin{align}
\frac{\partial w_1(\theta)}{\partial \theta}&=\frac{1}{{\sigma}^2(\theta)}\frac{\partial {\mu}(\theta)}{\partial \theta}   -\frac{{\mu}(\theta)}{{\sigma}^4(\theta)}\frac{\partial  {\sigma}^2(\theta) }{\partial \theta} ,\\
\frac{\partial w_2(\theta)}{\partial \theta}&=\frac{1}{2{\sigma}^4(\theta)}\frac{\partial  {\sigma}^2(\theta) }{\partial \theta}
\end{align}
and the derivatives of the expected sufficient statistics
\begin{align}
\frac{\partial \exdi{z;\theta}{t_1(z)}}{\partial \theta}&=\frac{\partial \mu(\theta)}{\partial \theta},\notag\\
\frac{\partial\exdi{z;\theta}{t_2(z)}}{\partial \theta}&=\frac{\partial{\sigma}^2(\theta)}{\partial \theta}+2{\mu}(\theta)\frac{\partial \mu(\theta)}{\partial \theta}
\end{align}
using the identity for the exponential family \eqref{theorem:fish:identity:vec} we obtain
\begin{align}\label{fish:identity:gauss}
F(\theta)&=\frac{\partial w_1(\theta)}{\partial \theta}\frac{\partial \exdi{z;\theta}{t_1(z)}}{\partial \theta}+\frac{\partial w_2(\theta)}{\partial \theta}\frac{\partial\exdi{z;\theta}{t_2(z)}}{\partial \theta}\notag\\
&=\frac{1}{{\sigma}^2(\theta)} \Big(\frac{\partial {\mu}(\theta)}{\partial \theta}\Big)^2 +\frac{1}{2{\sigma}^4(\theta)} \Big(\frac{\partial  {\sigma}^2(\theta) }{\partial \theta}\Big)^2.
\end{align}
It can be verified \cite[pp. 47]{Kay93} that \eqref{fish:identity:gauss} is the actual Fisher information measure for the parametric Gaussian model \eqref{parametric:gauss}.
\section{Strong Fisher Information Bounds}
\subsection{The Exponential Replacement}
If the parametric model $p(\ve{z};\ve{\theta})$ belongs to the exponential family, the $L$ natural parameters $w_l(\ve{\theta})$ and the associated expected values of the sufficient statistics $\exdi{\ve{z};\ve{\theta}}{  t_l(\ve{z}) }$ are known, the identity \eqref{theorem:fish:identity:vec} shows that the Fisher information measure can be computed by evaluating a simple sum. In the inconvenient situation where it is unclear if the model $p(\ve{z};\ve{\theta})$ belongs to the exponential family and the sufficient statistics $\ve{t}(\ve{z})$ or the natural parameters $\ve{w}(\ve{\theta})$ are unknown, the identity \eqref{fish:identity:vec:compact} can't be applied. The conceptual idea behind our approach for such a scenario is to replace the original system by an equivalent distribution $\tilde{p}(\ve{z};\ve{\theta})$ which is member of the exponential family. To this end, we select an arbitrary set of $L$ auxiliary statistics $\phi_l(\ve{z}): \fieldR^N \to \fieldR$, determine the expected values $\exdi{\ve{z};\ve{\theta}}{  \phi_l(\ve{z}) }$ with the original system $p(\ve{z};\ve{\theta})$ and choose the exponential family distribution $\tilde{p}(\ve{z};\ve{\theta})$ with sufficient statistics $t_l(\ve{z})=\phi_l(\ve{z})$ and equivalent expected values $\exdi{\ve{\tilde{z}};\ve{\theta}}{  t_l(\ve{z}) }=\exdi{\ve{z};\ve{\theta}}{  \phi_l(\ve{z}) }$.
\subsection{Fisher Information Bound}
Through the covariance inequality \cite[pp. 113]{Lehman03} we show that the Fisher information of the equivalent exponential family replacement $\tilde{p}(\ve{z};\ve{\theta})$ is always dominated by the information measure of the original system  $p(\ve{z};\ve{\theta})$.
\begin{theorem}[Fisher Information Bound]
For any probability density or mass function $p(\ve{z};\ve{\theta})$, any set of $L$ auxiliary statistics $\phi_l(\ve{z}): \fieldR^N \to \fieldR$ and with arbitrary weightings $\ve{b}_l(\ve{\theta})\in\fieldR^{M}$, the Fisher information matrix dominates
\begin{align}\label{ineq:covar:fisher}
&\ve{F}(\ve{\theta}) \succeq
\Bigg(\sum_{l=1}^{L} \bigg( \frac{\partial  \exdi{\ve{z};\ve{\theta}}{  \phi_l(\ve{z}) }}{\partial \ve{\theta}} \bigg)^{\T}  \ve{b}^{\T}_l(\ve{\theta})  \Bigg)  \notag\\
&\cdot\Bigg(\exdi{\ve{z};\ve{\theta}}{ \bigg( \sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \phi_l(\ve{z})\bigg) \bigg( \sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \phi_l(\ve{z}) \bigg)^{\T} }\notag\\
&-\bigg(\sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \exdi{\ve{z};\ve{\theta}}{  \phi_l(\ve{z}) } \bigg) \bigg( \sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \exdi{\ve{z};\ve{\theta}}{ \phi_l(\ve{z})}\bigg)^{\T}\Bigg)^{-1}\notag\\
&\cdot \bigg(\sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \frac{\partial  \exdi{\ve{z};\ve{\theta}}{  \phi_l(\ve{z}) }}{\partial \ve{\theta}} \bigg).
\end{align}
\end{theorem}
\begin{IEEEproof}
see Appendix \ref{append:fish:bound:vector}.
\end{IEEEproof}
Defining the matrix
\begin{align}
\ve{B}(\ve{\theta})=\begin{bmatrix} \ve{b}_1(\ve{\theta}) &\ve{b}_2(\ve{\theta}) &\ldots &\ve{b}_L(\ve{\theta})\end{bmatrix}^{\T},
\end{align}
the vector of auxiliary statistics
\begin{align}
\ve{\phi}(\ve{z})=\begin{bmatrix} \phi_1(\ve{z}) &\phi_2(\ve{z}) &\ldots &\phi_L(\ve{z})\end{bmatrix}^{\T},
\end{align}
and it's expected value
\begin{align}\label{mean:aux:stat}
\ve{\mu}_{\ve{\phi}}(\ve{\theta})=\exdi{\ve{z};\ve{\theta}}{\ve{\phi}(\ve{z})},
\end{align}
it is possible to write
\begin{align}
\sum_{l=1}^{L} \bigg( \frac{\partial  \exdi{\ve{z};\ve{\theta}}{  \phi_l(\ve{z}) }}{\partial \ve{\theta}} \bigg)^{\T}  \ve{b}^{\T}_l(\ve{\theta}) &= \bigg(\frac{\partial \ve{\mu}_{\ve{\phi}}(\ve{\theta})}{ \partial \ve{\theta}} \bigg)^{\T} \ve{B}(\ve{\theta}),\\
\sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \phi_l(\ve{z}) &= \ve{B}^{\T}(\ve{\theta})\ve{\phi}(\ve{z}),\\
\sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \exdi{\ve{z};\ve{\theta}}{  \phi_l(\ve{z}) } &= \ve{B}^{\T}(\ve{\theta})\ve{\mu}_{\ve{\phi}}(\ve{\theta}).
\end{align}
Further, defining
\begin{align}\label{covariance:aux:stat}
\ve{R}_{\ve{\phi}}(\ve{\theta}) &=\exdi{\ve{z};\ve{\theta}}{ \ve{\phi}(\ve{z}) \ve{\phi}^{\T}(\ve{z}) } - \ve{\mu}_{\ve{\phi}}(\ve{\theta}) \ve{\mu}^{\T}_{\ve{\phi}}(\ve{\theta})
\end{align}
the bound \eqref{ineq:covar:fisher} can be reformulated
\begin{align}\label{ineq:covar:fisher:compact}
&\ve{F}(\ve{\theta}) \succeq \notag\\
&\bigg(\frac{\partial \ve{\mu}_{\ve{\phi}}(\ve{\theta})}{ \partial \ve{\theta}} \bigg)^{\T} \ve{B}(\ve{\theta})\big(\ve{B}^{\T}(\ve{\theta}) \ve{R}_{\ve{\phi}}(\ve{\theta})\ve{B}(\ve{\theta})\big)^{-1}\ve{B}^{\T}(\ve{\theta}) \frac{\partial \ve{\mu}_{\ve{\phi}}(\ve{\theta})}{ \partial \ve{\theta}}.
\end{align}
\subsection{Tight Information Bound}
Note, that it is left to choose the auxiliary statistics $\ve{\phi}(\ve{z})$ and optimize the associated weighting factors $\ve{B}(\ve{\theta})$ such that the right hand side of the information bound \eqref{ineq:covar:fisher:compact} is maximized in the matrix sense. If the underlying statistical model $p(\ve{z};\ve{\theta})$ is from the exponential family type  \eqref{exp:family:vec} and the sufficient statistics $\ve{\phi}(\ve{z})=\ve{t}(\ve{z})$ are used for the approximation \eqref{ineq:covar:fisher:compact}, it is possible to obtain a tight information bound.
\begin{theorem}[Tightness of the Fisher Information Bound]
If the probability density or mass function $p(\ve{z};\ve{\theta})$ belongs to the exponential family \eqref{exp:family:vec} and the $L$ auxiliary functions $\ve{\phi}(\ve{z})$ are the sufficient statistics $\ve{t}(\ve{z})$ of the statistical model $p(\ve{z};\ve{\theta})$, the optimization of the weighting matrix $\ve{B}^{\star}(\ve{\theta})$ in \eqref{ineq:covar:fisher:compact} will lead to a tight Fisher information bound.
\end{theorem}
\begin{IEEEproof}
see Appendix \ref{append:fish:bound:vector:tight}.
\end{IEEEproof}
\section{Optimization of the Information Bound}
In the following we discuss how to perform the optimization of the right hand side of \eqref{ineq:covar:fisher} when the sufficient statistics $\ve{t}(\ve{z})$ are unknown and we have to resort to the auxiliary statistics $\ve{\phi}(\ve{z})$. Substituting
\begin{align}\label{matrix:substitution}
\ve{B}(\ve{\theta})=\ve{R}_{\ve{\phi}}^{-\frac{1}{2}}(\ve{\theta}) \ve{B}'(\ve{\theta})
\end{align}
in \eqref{ineq:covar:fisher:compact} under the constraint that
\begin{align}\label{matrix:constraint}
\ve{B}'^{\T}(\ve{\theta})\ve{B}'(\ve{\theta})=\ve{I},
\end{align}
we obtain a modified bound
\begin{align}
&\ve{F}(\ve{\theta}) \succeq \notag\\
&\bigg(\frac{\partial \ve{\mu}_{\ve{\phi}}(\ve{\theta})}{ \partial \ve{\theta}} \bigg)^{\T} \ve{R}_{\ve{\phi}}^{-\frac{1}{2}}(\ve{\theta}) \ve{B}'(\ve{\theta})\ve{B}'^{\T}(\ve{\theta}) \ve{R}_{\ve{\phi}}^{-\frac{1}{2}}(\ve{\theta}) \frac{\partial \ve{\mu}_{\ve{\phi}}(\ve{\theta})}{ \partial \ve{\theta}}.
\end{align}
The right hand side is maximized in the matrix sense under the constraint \eqref{matrix:constraint} with
\begin{align}\label{matrix:maximizer}
&\ve{B}'(\ve{\theta})=\notag\\
&\ve{R}_{\ve{\phi}}^{-\frac{1}{2}}(\ve{\theta}) \frac{\partial \ve{\mu}_{\ve{\phi}}(\ve{\theta})}{ \partial \ve{\theta}} \Bigg(\bigg(\frac{\partial \ve{\mu}_{\ve{\phi}}(\ve{\theta})}{ \partial \ve{\theta}} \bigg)^{\T}  \ve{R}_{\ve{\phi}}^{-1}(\ve{\theta}) \frac{\partial \ve{\mu}_{\ve{\phi}}(\ve{\theta})}{ \partial \ve{\theta}} \Bigg)^{-\frac{1}{2}}.
\end{align}
\begin{theorem}[Strong Fisher Information Bound]
For any probability density or mass function $p(\ve{z};\ve{\theta})$ and any set of $L$ auxiliary statistics $\ve{\phi}(\ve{z})$, with the definitions \eqref{mean:aux:stat} and \eqref{covariance:aux:stat}, the Fisher information matrix $F(\theta)$ dominates
\begin{align}\label{bound:eigenvalues}
\ve{F}(\ve{\theta}) \succeq \bigg(\frac{\partial \ve{\mu}_{\ve{\phi}}(\ve{\theta})}{ \partial \ve{\theta}} \bigg)^{\T} \ve{R}_{\ve{\phi}}^{-1}(\ve{\theta}) \frac{\partial \ve{\mu}_{\ve{\phi}}(\ve{\theta})}{ \partial \ve{\theta}}.
\end{align}
\end{theorem}
\begin{IEEEproof}
Follows from using \eqref{matrix:maximizer} and \eqref{matrix:substitution} in \eqref{ineq:covar:fisher:compact}.
\end{IEEEproof}
Note, that due to the tightness of the bound \eqref{bound:eigenvalues} for exponential family models \eqref{exp:family:vec}, besides \eqref{fish:identity:vec:compact} we have an additional identity for such kind of system models
\begin{align}
\ve{F}(\ve{\theta}) = \bigg(\frac{\partial \ve{\mu}_{\ve{t}}(\ve{\theta})}{ \partial \ve{\theta}} \bigg)^{\T} \ve{R}_{\ve{t}}^{-1}(\ve{\theta}) \frac{\partial \ve{\mu}_{\ve{t}}(\ve{\theta})}{ \partial \ve{\theta}}.
\end{align}
In the following we will demonstrate applications of the presented main result \eqref{bound:eigenvalues}. For simplicity, we focus on univariate problems $z\in\fieldR$ with a single parameter $\theta\in\fieldR$.
\section{Information Bound with $L$ Moments}
In order to test a generalization of the approach \eqref{bound:fisher:two} with the derivatives of $L$ raw moments, we use the obtained result \eqref{bound:eigenvalues} under the convention
\begin{align}
\phi_l(z)=z^l
\end{align}
and apply the resulting information bound to the log-normal \eqref{model:lognormal} and the Weilbull model \eqref{model:weilbull}. The required expectations \eqref{mean:aux:stat} and \eqref{covariance:aux:stat} of the auxiliary statistics $\phi_l(z)$ are directly available by the fact that for the log-normal distribution the $l$-th raw moment is given by
\begin{align}
\exdi{z;\theta}{z^l}&=\expe{l\theta+\frac{1}{2} l^2 \sigma^2}
\end{align}
and therefore it's derivative is defined by
\begin{align}
\frac{\partial \exdi{z;\theta}{z^l}}{\partial \theta}&=l\expe{l\theta+\frac{1}{2} l^2 \sigma^2}.
\end{align}
Accordingly, for the Weibull distribution we have
\begin{align}
\exdi{z;\theta}{z^l}=\theta^l\Gamma_l,
\end{align}
and
\begin{align}
\frac{\partial \exdi{z;\theta}{z^l}}{\partial \theta}=l\theta^{l-1}\Gamma_l.
\end{align}
\begin{figure}[!htbp]
\begin{tikzpicture}[scale=1.0]

  	\begin{axis}[ylabel=$10 \log{\chi(\theta)}$,
  			xlabel=$k$,
			grid,
			ymin=-0.2,
			ymax=0.1,
			xmin=1,
			xmax=10,
			legend pos=south east]

    			\addplot[black, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=otimes*, mark repeat=5] table[x index=0, y index=1]{Weibull_Loss_Lmo_K1.txt};
			\addlegendentry{$L=1$};

			\addplot[green, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=diamond*, mark repeat=5] table[x index=0, y index=1]{Weibull_Loss_Lmo_K2.txt};
			\addlegendentry{$L=2$};

    			\addplot[blue, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=square*, mark repeat=5] table[x index=0, y index=1]{Weibull_Loss_Lmo_K3.txt};
			\addlegendentry{$L=3$};

    			\addplot[red, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=star, mark repeat=5] table[x index=0, y index=1]{Weibull_Loss_Lmo_K4.txt};
			\addlegendentry{$L=4$};

\end{axis}
\end{tikzpicture}
\caption{Information Bound with $L$ Moments  - Weibull Distribution}
\label{LMo_Weibull_Loss}
\end{figure}
\begin{figure}[!htbp]
\begin{tikzpicture}[scale=1.0]

  	\begin{axis}[ylabel=$10 \log{\chi(\theta)}$,
  			xlabel=$\sigma$,
			grid,
			ymin=-3.1,
			ymax=0.1,
			xmin=0.2,
			xmax=1.2,
			legend pos=south west]

    			\addplot[black, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=otimes*, mark repeat=5] table[x index=0, y index=1]{LogNormal_Loss_Lmo_L1.txt};
			\addlegendentry{$L=1$};

			\addplot[green, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=diamond*, mark repeat=5] table[x index=0, y index=1]{LogNormal_Loss_Lmo_L2.txt};
			\addlegendentry{$L=2$};

    			\addplot[blue, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=square*, mark repeat=5] table[x index=0, y index=1]{LogNormal_Loss_Lmo_L3.txt};
			\addlegendentry{$L=3$};

    			\addplot[red, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=star, mark repeat=5] table[x index=0, y index=1]{LogNormal_Loss_Lmo_L4.txt};
			\addlegendentry{$L=4$};

\end{axis}
\end{tikzpicture}
\caption{Information Bound with $L$ Moments - Log-normal Distribution}
\label{LMo_LogNormal_Loss}
\end{figure}
In Fig. \ref{LMo_Weibull_Loss}  and Fig. \ref{LMo_LogNormal_Loss} we depict the approximation loss
\begin{align}
\chi(\theta)=\frac{\big(\frac{\partial \ve{\mu}_{\ve{\phi}}(\theta)}{ \partial {\theta}} \big)^{\T} \ve{R}_{\ve{\phi}}^{-1}(\theta) \frac{\partial \ve{\mu}_{\ve{\phi}}(\theta)}{ \partial {\theta}}}{F(\theta)}
\end{align}
for different values $L$.
For the Weibull distribution, for which the result is depicted in Fig. \ref{LMo_Weibull_Loss}, we observe, that the bound with $L=1$ \cite{Stein14} and $L=2$ \cite{Stein15} can be significantly improved by incorporating the derivatives of the third and the fourth moment. In contrast for the log-normal distribution (see Fig. \ref{LMo_LogNormal_Loss}) taking into account more than the first two raw moments results only in a slight performance improvement.
In order to visualize the result of the optimization of the bound \eqref{ineq:covar:fisher:compact}, in Fig. \ref{4Mo_Weibull_Weights} the normalized absolute weights
\begin{align}
\bar{b}^{\star}_{l}(\theta)= \frac{|{b}^{\star}_{l}(\theta)|}{\sum_{l=1}^{L}|{b}^{\star}_{l}(\theta)|}
\end{align}
are plotted for the Weibull example, where with \eqref{matrix:substitution} and \eqref{matrix:maximizer}
\begin{align}
{b}_{l}^{\star}(\theta)&=\frac{ \Big[ \ve{R}_{\ve{\phi}}^{-1}(\theta)\frac{\partial \ve{\mu}_{\ve{\phi}}(\theta)}{\partial \theta} \Big]_l}{\sqrt{\big(\frac{\partial \ve{\mu}_{\ve{\phi}}(\theta)}{\partial \theta}\big)^{\T}  \ve{R}_{\ve{\phi}}^{-1}(\theta)\frac{\partial \ve{\mu}_{\ve{\phi}}(\theta)}{\partial \theta}}}.
\end{align}
The individual normalized weights $\bar{b}^{\star}_{l}(\theta)$ indicate the importance of the corresponding auxiliary statistic $\phi_l(\theta)=z^l$ in the approximation of the Fisher information. It can be seen, that for the Weibull distribution the sufficient statistic $z^k$ attains the full weight in the cases $k=1,2,3,4$.
\begin{figure}[!htbp]
\begin{tikzpicture}[scale=1.0]

  	\begin{axis}[ylabel=$\beta^{\star}_l(\theta)$,
			xlabel=$k$,
			grid,
			ymin=0,
			ymax=1,
			xmin=1,
			xmax=10,
			legend pos=north east]

    			\addplot[black, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=otimes*, mark repeat=5] table[x index=0, y index=1]{Weibull_4mo_Weight.txt};
			\addlegendentry{$\bar{b}^{\star}_1(\theta)$};

			\addplot[green, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=diamond*, mark repeat=5] table[x index=0, y index=2]{Weibull_4mo_Weight.txt};
			\addlegendentry{$\bar{b}^{\star}_2(\theta)$};

    			\addplot[blue, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=square*, mark repeat=5] table[x index=0, y index=3]{Weibull_4mo_Weight.txt};
			\addlegendentry{$\bar{b}^{\star}_3(\theta)$};

    			\addplot[red, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=star, mark repeat=5] table[x index=0, y index=4]{Weibull_4mo_Weight.txt};
			\addlegendentry{$\bar{b}^{\star}_4(\theta)$};

\end{axis}
\end{tikzpicture}
\caption{Information Bound with $L$ Moments  - Weibull Weights ($L=4$)}
\label{4Mo_Weibull_Weights}
\end{figure}
In contrast for the log-normal distribution it is observed in Fig. \ref{4Mo_LogNormal_Weights} that none of the moments obtains full weight. However, the first moment plays a dominant role, in particular when $\sigma>0.6$. Note, that for the log-normal distribution with known scale parameter $\sigma$, $\ln z$ is a sufficient statistic. Incorporating this statistics into the approximation by using $\phi_1(z)=\ln z$ would change the situation and provide a tight approximation for $F(\theta)$.
\begin{figure}[!htbp]
\begin{tikzpicture}[scale=1.0]

  	\begin{axis}[ylabel=$\beta^{\star}_l(\theta)$,
			xlabel=$\sigma$,
			grid,
			ymin=0,
			ymax=1,
			xmin=0.2,
			xmax=1.2,
			legend pos=north east]

    			\addplot[black, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=otimes*, mark repeat=10] table[x index=0, y index=1]{LogNormal_4mo_Weights.txt};
			\addlegendentry{$\bar{b}^{\star}_1(\theta)$};

			\addplot[green, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=diamond*, mark repeat=10] table[x index=0, y index=2]{LogNormal_4mo_Weights.txt};
			\addlegendentry{$\bar{b}^{\star}_2(\theta)$};

    			\addplot[blue, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=square*, mark repeat=10] table[x index=0, y index=3]{LogNormal_4mo_Weights.txt};
			\addlegendentry{$\bar{b}^{\star}_3(\theta)$};

    			\addplot[red, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=star, mark repeat=10] table[x index=0, y index=4]{LogNormal_4mo_Weights.txt};
			\addlegendentry{$\bar{b}^{\star}_4(\theta)$};

\end{axis}
\end{tikzpicture}
\caption{Information Bound with $L$ Moments - Log-normal Weights ($L=4$)}
\label{4Mo_LogNormal_Weights}
\end{figure}
\section{Learning Nonlinear Statistical Models}
The previous section indicates an interesting application of the presented result \eqref{bound:eigenvalues}. Being able to describe the expectations \eqref{mean:aux:stat} and \eqref{covariance:aux:stat} for an arbitrary statistical model $p(z;\theta)$ with an arbitrary set of auxiliary functions $\ve{\phi}(z)$, the optimization result \eqref{matrix:maximizer} can be used in order to identify, among the auxiliary statistics $\ve{\phi}(z)$, candidates for the sufficient statistics of the system model. Further, together with \eqref{mean:aux:stat} and \eqref{covariance:aux:stat}, the information bound \eqref{bound:eigenvalues} allows to specify the minimum inference capability that can be guaranteed to be achievable for the model of interest.
\subsection{Practical Impact}
This has high practical relevance as in real-world applications technical systems are subject to various random and non-linear effects. Under such circumstances an accurate analytical description of the probability density or mass function $p(z;\theta)$ is usually hard to obtain. Nevertheless, to be able to identify transformations of the data exhibiting high information content is attractive in such a situation. Such functions can be used for applications like efficient data compression and for the formulation of high-resolution estimation algorithms. Further, a conservative approximation of the Fisher information measure like \eqref{bound:eigenvalues} allows to benchmark the performance of estimation algorithms on the system under investigation or to identify system layouts of high estimation theoretic quality. If the system parameter $\theta$ can be controlled in a calibrated setup, the entities \eqref{mean:aux:stat} and \eqref{covariance:aux:stat} can be determined for any system $p(z;\theta)$ by measurements at the system output $z$.
\subsection{Example - Nonlinear Amplification Device}
We demonstrate such a measurement-driven learning approach in a calibrated setup with the example of a solid-state power amplifier. The system parameter $\theta$ of interest is assumed to be the direct-current offset (the mean) at the input of the nonlinear device. For the mapping from the input $x$ to the output $y$ of the amplifier, we apply the Rapp model \cite{Rapp91}
\begin{align}
y=\frac{x}{(1+|x|^{2\rho})^{\frac{1}{2\rho}}},
\label{def:Rapp}
\end{align}
where $\rho$ is a smoothness factor. Fig. \ref{Rapp:io} depicts the input-to-output relation of this nonlinear system model with different values of $\rho$.
\begin{figure}[!htbp]
\begin{tikzpicture}[scale=1.0]

  	\begin{axis}[ylabel=$z$,
  			xlabel=$x$,
			grid,
			ymin=-1.0,
			ymax=1.0,
			xmin=-4,
			xmax=4,
			legend pos=south east]

			\addplot[red, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=otimes*, mark repeat=3] table[x index=0, y index=5]{Rapp_IO.txt};
			\addlegendentry{$\rho=8.00$}

			\addplot[magenta, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=diamond*, mark repeat=3] table[x index=0, y index=4]{Rapp_IO.txt};
			\addlegendentry{$\rho=2.00$}

			\addplot[green, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=square*, mark repeat=3] table[x index=0, y index=3]{Rapp_IO.txt};
			\addlegendentry{$\rho=1.00$}

			\addplot[blue, style=solid, line width=0.75pt,smooth,smooth,every mark/.append style={solid}, mark=star, mark repeat=3] table[x index=0, y index=2]{Rapp_IO.txt};
			\addlegendentry{$\rho=0.50$}

			\addplot[black, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=triangle*, mark repeat=3] table[x index=0, y index=1]{Rapp_IO.txt};
			\addlegendentry{$\rho=0.25$}

\end{axis}
\end{tikzpicture}
\caption{Rapp Model - Input-to-Output}
\label{Rapp:io}
\end{figure}
We apply a Gaussian input
\begin{align}
x=\theta+\eta
\end{align}
with $\eta\sim\mathcal{N}(0,1)$ to the nonlinear system \eqref{def:Rapp}, we set
\begin{align}
\ve{\phi}(z)=
\begin{bmatrix}
z &z^2 &z^3 &z^4 &|z| &\ln |z|
\end{bmatrix}^{\T}
\end{align}
and for each value of $\theta$ we approximate the expectations \eqref{mean:aux:stat} and \eqref{covariance:aux:stat}, with $10^8$ realizations of the system output, by their sample mean.
\begin{figure}[!htbp]
\begin{tikzpicture}[scale=1.0]

  	\begin{axis}[ylabel=$10 \log{\chi(\theta)}$,
  			xlabel=$\theta$,
			grid,
			ymin=-6.0,
			ymax=0,
			xmin=0,
			xmax=4.0,
			legend pos=south west]

    			\addplot[black, style=solid, line width=0.75pt, every mark/.append style={solid}, mark=diamond*, mark repeat=5] table[x index=0, y index=1]{LearnRapp_2.txt};

\end{axis}
\end{tikzpicture}
\caption{Performance Loss - Rapp model with $\rho=2.0$}
\label{Rapp_Loss}
\end{figure}
Fig. \ref{Rapp_Loss} shows the obtained performance loss
\begin{align}
\chi(\theta)=\frac{\big(\frac{\partial \ve{\mu}_{\ve{\phi}}(\theta)}{ \partial {\theta}} \big)^{\T} \ve{R}_{\ve{\phi}}^{-1}(\theta) \frac{\partial \ve{\mu}_{\ve{\phi}}(\theta)}{ \partial {\theta}}}{F_{x}(\theta)}
\end{align}
which is introduced by the nonlinear Rapp model with smoothness factor $\rho=2.0$. Note, that $F_x(\theta)$ is the Fisher information with respect to $\theta$ at the input $x$ of the nonlinear Rapp model. It is observed that for an input mean $\theta>2.0$, the saturation of the non-linear Rapp model introduces a significant information loss. Fig. \ref{Rapp_Weights} shows the normalized absolute weights $\bar{b}^{\star}_l$ associated with each statistic, which have been attained by the optimization \eqref{matrix:maximizer}.
\begin{figure}[!htbp]
\begin{tikzpicture}[scale=1.0]

  	\begin{axis}[ylabel=$\beta^{\star}_l(\theta)$,
			xlabel=$\theta$,
			grid,
			ymin=0,
			ymax=0.5,
			xmin=0,
			xmax=4,
			legend pos=south east]

    			\addplot[black, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=otimes*, mark repeat=3] table[x index=0, y index=2]{LearnRapp_2.txt};
			\addlegendentry{$\bar{b}^{\star}_1(\theta)$};

			\addplot[green, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=diamond*, mark repeat=3] table[x index=0, y index=3]{LearnRapp_2.txt};
			\addlegendentry{$\bar{b}^{\star}_2(\theta)$};

    			\addplot[blue, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=square*, mark repeat=3] table[x index=0, y index=4]{LearnRapp_2.txt};
			\addlegendentry{$\bar{b}^{\star}_3(\theta)$};

    			\addplot[magenta, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=star, mark repeat=3] table[x index=0, y index=5]{LearnRapp_2.txt};
			\addlegendentry{$\bar{b}^{\star}_4(\theta)$};

			\addplot[red, style=solid, line width=0.75pt,smooth,every mark/.append style={solid}, mark=triangle*, mark repeat=3] table[x index=0, y index=6]{LearnRapp_2.txt};
			\addlegendentry{$\bar{b}^{\star}_5(\theta)$};

			\addplot[cyan, style=solid, line width=0.75pt,smooth,every mark/.append style={solid, fill=white}, mark=*, mark repeat=3] table[x index=0, y index=7]{LearnRapp_2.txt};
			\addlegendentry{$\bar{b}^{\star}_6(\theta)$};

\end{axis}
\end{tikzpicture}
\caption{Optimization Weights - Rapp model with $\rho=2.0$}
\label{Rapp_Weights}
\end{figure}
It can be seen that the second moment and the expected absolute value play a dominant role in the approximation of the Fisher information.
\section{Inference from Compressed Observations}
Finally, we address the question how to perform estimation of the system parameter $\theta$ under an unknown statistical system model $p(z;\theta)$ after having learned the properties $\ve{\mu}_{\ve{\phi}}(\theta)$ and $\ve{R}_{\ve{\phi}}(\theta)$ with an arbitrary set of $L$ auxiliary statistics $\ve{\phi}(z)$ by calibrated measurements.
\subsection{Nonlinear Compression}
Observing $N$ samples at the system output of $p(z;\theta)$, the data vector $\ve{z}\in\fieldR^N$ is available. First we apply compression by using the auxiliary statistics $\ve{\phi}(z)$ to form the sample mean
\begin{align}
\ve{\tilde{\phi}}=\frac{1}{N}\sum_{n=1}^{N} \ve{\phi}(z_n)
\end{align}
and subsequently discarding the original data $\ve{z}$. Note that this reduces the size of the data by a factor of $\frac{N}{L}$.
\subsection{Conservative Maximum-Likelihood Estimation}
If the analytic characterization of $p(z;\theta)$ is available, given the data $\ve{z}$ one usually resorts to the asymptotically optimum estimator based on the maximization of the log-likelihood
\begin{align}
\hat{\theta}_{\text{ML}}(\ve{z})&=\arg\max_{\theta \in \Theta} \ln p(\ve{z};\theta)\notag\\
&=\arg\max_{\theta \in \Theta} \sum_{n=1}^{N}\ln p(z_n;\theta).
\end{align}
This is not possible if no description of the model $p(z;\theta)$ is available. We propose to replace the original system $p(z;\theta)$ by a distribution of the exponential family \eqref{exp:family:vec} with log-likelihood function
\begin{align}\label{replacement:distribution}
\ln \tilde{p}(z;\theta)= \ve{\beta}^{\T}(\theta) \ve{\phi}(z) - \lambda({\theta})+\kappa({z}).
\end{align}
Note, that the replacement \eqref{replacement:distribution} has as sufficient statistics the auxiliary functions $\ve{\phi}(z)$. Assuming arbitrary weighting factors $\ve{b}(\theta)=\frac{\partial \ve{\beta}(\theta)}{\partial \theta}$, the score function takes the form
\begin{align}
\frac{\partial \ln \tilde{p}(z;\theta)}{\partial \theta} &= \ve{b}^{\T}(\theta) \big(\ve{\phi}(z) - \ve{\mu}_{\ve{\phi}}(\theta)\big).
\end{align}
A conservative maximum-likelihood estimate (CMLE) is then found by setting the score of the data to zero, i.e.
\begin{align}\label{approx:score:zero}
\left.\frac{\partial \ln \tilde{p}(\ve{z};{\theta})}{\partial \theta} \right|_{\theta=\hat{\theta}} = 0.
\end{align}
For a memoryless system, the receive score can be written
\begin{align}
\frac{\partial \ln \tilde{p}(\ve{z};\theta)}{\partial \theta} &= \sum_{n=1}^{N} \frac{\partial \ln \tilde{p}({z}_{n};\theta)}{\partial \theta}\notag\\
&=\ve{b}^{\T}(\theta) \sum_{n=1}^{N}  \big(\ve{\phi}(z_n) - \ve{\mu}_{\ve{\phi}}(\theta)\big)\notag\\
&=  N\ve{b}^{\T}(\theta) \big(\ve{\tilde{\phi}} - \ve{\mu}_{\ve{\phi}}(\theta)\big),
\end{align}
such that the CMLE $\hat{\theta}(\ve{\tilde{\phi}})$ is found by solving
\begin{align}\label{cmle:calculation}
 \ve{b}^{\T}({\theta}) \big(\ve{\tilde{\phi}} - \ve{\mu}_{\ve{\phi}}({\theta})\big)=0
\end{align}
with respect to $\theta$. Note, that for the calculation of the CMLE \eqref{cmle:calculation}, access to the original data $\ve{z}$ is not required.
\begin{theorem}[Asymptotic Performance and Consistency]
If $p(z;\theta)$ is the data-generating model, the CMLE $\hat{\theta}(\ve{\tilde{\phi}})$, asymptotically in $N$, produces consistent estimates which are Gaussian distributed
\begin{align}\label{amle:performance:unopt}
\hat{\theta} \sim \mathcal{N}\left({\theta},\frac{1}{N}\frac{\ve{b}^{\T}(\theta) \ve{R}_{\ve{\phi}}(\theta) \ve{b}(\theta)}{ \big(\frac{\partial \ve{\mu}_{\ve{\phi}}(\theta)}{\partial \theta}\big)^{\T} \ve{b}(\theta) \ve{b}^{\T}(\theta) \frac{\partial \ve{\mu}_{\ve{\phi}}(\theta)}{\partial \theta}}\right).
\end{align}
\end{theorem}
\begin{IEEEproof}
see Appendix \ref{append:amle:perf:unbias}.
\end{IEEEproof}
Using the best weighting
\begin{align}
\ve{b}^{\star}({\theta})=\ve{R}_{\ve{\phi}}^{-1}({\theta}) \frac{\partial \ve{\mu}_{\ve{\phi}}({\theta})}{ \partial {\theta}} \Bigg(\bigg(\frac{\partial \ve{\mu}_{\ve{\phi}}({\theta})}{ \partial {\theta}} \bigg)^{\T}  \ve{R}_{\ve{\phi}}^{-1}({\theta}) \frac{\partial \ve{\mu}_{\ve{\phi}}({\theta})}{ \partial {\theta}} \Bigg)^{-\frac{1}{2}}
\end{align}
for the auxiliary statistics, the CMLE is found by solving
\begin{align}
\frac{\big(\frac{\partial \ve{\mu}_{\ve{\phi}}({\theta})}{ \partial {\theta}} \big)^{\T}\ve{R}_{\ve{\phi}}^{-1}({\theta})}{\sqrt{\Big(\big(\frac{\partial \ve{\mu}_{\ve{\phi}}({\theta})}{ \partial {\theta}} \big)^{\T}  \ve{R}_{\ve{\phi}}^{-1}({\theta}) \frac{\partial \ve{\mu}_{\ve{\phi}}({\theta})}{ \partial {\theta}} \Big)}} \big(\ve{\tilde{\phi}} - \ve{\mu}_{\ve{\phi}}({\theta})\big)
=0.
\end{align}
The estimator then achieves a performance equivalent to the inverse of the approximation for the Fisher information measure \eqref{bound:eigenvalues}.
\begin{corollary}[CMLE - Optimized Asymptotic Performance]
If $p(z;\theta)$ is the data-generating model, the CMLE $\hat{\theta}(\ve{\tilde{\phi}})$ calculated with optimized weights $ \ve{b}^{\star}({\theta})$, asymptotically in $N$, produces consistent estimates which are Gaussian distributed
\begin{align}\label{amle:performance:opt}
\hat{\theta} \sim \mathcal{N}\Bigg({\theta},\frac{1}{N} \frac{1}{\big(\frac{\partial \ve{\mu}_{\ve{\phi}}({\theta})}{ \partial {\theta}} \big)^{\T} \ve{R}_{\ve{\phi}}^{-1}({\theta}) \frac{\partial \ve{\mu}_{\ve{\phi}}({\theta})}{ \partial {\theta}}}\Bigg).
\end{align}
\end{corollary}
\begin{IEEEproof}
Follows from the fact that with \eqref{matrix:maximizer}
\begin{align}
\frac{\ve{b}^{\star\T}(\theta) \ve{R}_{\ve{\phi}}(\theta) \ve{b}^{\star}(\theta)}{ \big(\frac{\partial \ve{\mu}_{\ve{\phi}}(\theta)}{\partial \theta}\big)^{\T} \ve{b}^{\star}(\theta) \ve{b}^{\star\T}(\theta) \frac{\partial \ve{\mu}_{\ve{\phi}}(\theta)}{\partial \theta}}
=\frac{1}{\big(\frac{\partial \ve{\mu}_{\ve{\phi}}(\ve{\theta})}{ \partial \ve{\theta}} \big)^{\T} \ve{R}_{\ve{\phi}}^{-1}(\ve{\theta}) \frac{\partial \ve{\mu}_{\ve{\phi}}(\ve{\theta})}{ \partial \ve{\theta}}}.
\end{align}
\end{IEEEproof}
\subsection{Connection to the Generalized Method of Moments}
By squaring the CMLE \eqref{cmle:calculation} can be reformulated
\begin{align}
\hat{\theta}(\ve{\tilde{\phi}})&=\arg \min_{\theta\in\Theta}\big(\ve{\tilde{\phi}} - \ve{\mu}_{\ve{\phi}}({\theta})\big)^{\T}\ve{b}({\theta})\ve{b}^{\T}({\theta}) \big(\ve{\tilde{\phi}} - \ve{\mu}_{\ve{\phi}}({\theta})\big)
\end{align}
which is identified as special case of Hansen's estimator \cite{Hansen82}
\begin{align}\label{estimator:hansen}
\hat{\theta}(\ve{z})&=\arg \min_{\theta\in\Theta} \Big(\frac{1}{N} \sum_{n=1}^{N} \ve{f}(z_n;\theta)\Big)^{\T}\ve{D}({\theta}) \Big(\frac{1}{N} \sum_{n=1}^{N} \ve{f}(z_n;\theta)\Big)
\end{align}
with the moment condition
\begin{align}
\ve{f}(z;\theta)=\ve{\phi}(z)-\ve{\mu}_{\ve{\phi}}({\theta})
\end{align}
and an optimized weighting matrix
\begin{align}
\ve{D}^{\star}({\theta})=\frac{ \ve{R}_{\ve{\phi}}^{-1}({\theta}) \frac{\partial \ve{\mu}_{\ve{\phi}}({\theta})}{ \partial {\theta}}\big(\frac{\partial \ve{\mu}_{\ve{\phi}}({\theta})}{ \partial {\theta}} \big)^{\T}\ve{R}_{\ve{\phi}}^{-1}({\theta})}{{\Big(\big(\frac{\partial \ve{\mu}_{\ve{\phi}}({\theta})}{ \partial {\theta}} \big)^{\T}  \ve{R}_{\ve{\phi}}^{-1}({\theta}) \frac{\partial \ve{\mu}_{\ve{\phi}}({\theta})}{ \partial {\theta}} \Big)}}.
\end{align}
The generalized method of moments \eqref{estimator:hansen} is an extension of the classical method of moments \cite{Pearson94}, derived by considering orthogonality conditions \cite{Hansen82}
\begin{align}
\exdi{z;\theta_t}{\ve{f}(z;\theta_t)}=\ve{0}
\end{align}
under the true parameter $\theta_t$. It is interesting to observe, that we obtain the method as a straightforward maximum-likelihood estimator after approximating the original system model $p({z};{\theta})$ through a set of auxiliary statistics $\ve{\phi}(z)$ by the closest (in the sense of the Fisher information measure) equivalent distribution $\tilde{p}({z};{\theta})$ within the exponential family. Therefore the equivalent exponential replacement provides a potential unifying link, like subtly requested by \cite{Bera02}, between Pearson's method of moments \cite{Pearson94} and Fisher's competing concept of likelihood \cite{Fisher22}.
\section{Conclusion}
We have established a generic approach for the construction of strong lower bounds for the Fisher information measure. For an arbitrary statistical model $p(\ve{z};\ve{\theta})$ and a set of auxiliary statistics $\ve{\phi}(\ve{z})$, we approximate the original likelihood function by using the likelihood function $\tilde{p}(\ve{z};\ve{\theta})$ of a distribution which belongs to the exponential family and exhibits the chosen set of auxiliary statistics $\ve{\phi}(\ve{z})$ as it's sufficient statistics $\ve{t}(\ve{z})$. Such a replacement model exhibits lower Fisher information than the original data-generating probability law. Through optimization of the derivative of the natural parameters of the replacement model we find the closest (in the Fisher information sense) replacement in the exponential family and therefore obtain an accurate pessimistic approximation for the Fisher information measure. The presented method has the advantage that the statistical output model $p(z;\theta)$ does not have to be known. The required expected values of the statistics can be learned from calibrated measurements at the system output. We have demonstrated that based on the learned auxiliary statistics and their optimized weights, consistent estimators can be formulated which achieve a performance that is equivalent to the inverse of the presented pessimistic approximation of the Fisher information measure. This forms a framework for the problem of parameter estimation with nonlinear systems. For any system model the results can be used to conservatively predict the theoretically possible estimation performance by analytical derivation or calibrated measurement of the mean and covariance of the auxiliary statistics. In addition our discussion provides a straightforward derivation of a data-efficient estimation method linked to the method of moments which allows to achieve the pessimistically approximated asymptotic system performance in practice.
\appendices

\section{Fisher Information Bound}\label{append:fish:bound:vector}
\begin{IEEEproof}
With $L$ functions $\phi_l(\ve{z})$, $L$ weighting vectors $\ve{b}_l(\ve{\theta})$ and a normalizer $\ve{\alpha}(\ve{\theta})$, we have
\begin{align}
&\int_{\mathcal{Z}} \bigg( \frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} \bigg)^{\T} \bigg( \sum_{l=1}^{L} \ve{b}^{\T}_l(\ve{\theta}) \phi_l(\ve{z}) - \ve{\alpha}^{\T}(\ve{\theta}) \bigg) p(\ve{z};\ve{\theta}) {\rm d}\ve{z}\notag\\
&=\quad \sum_{l=1}^{L} \bigg( \int_{\mathcal{Z}}  \frac{\partial p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} \phi_l(\ve{z}) {\rm d}\ve{z} \bigg)^{\T} \ve{b}^{\T}_l(\ve{\theta})\notag\\
&=\quad \sum_{l=1}^{L}   \bigg(\frac{\partial \exdi{\ve{z};\ve{\theta}}{\phi_l(\ve{z})} )}{\partial \ve{\theta}} \bigg)^{\T} \ve{b}^{\T}_l(\ve{\theta})
\end{align}
and
\begin{align}
&\int_{\mathcal{Z}}  \bigg( \sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \phi_l(\ve{z}) - \ve{\alpha}(\ve{\theta}) \bigg) \bigg( \frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} \bigg) p(\ve{z};\ve{\theta}) {\rm d}\ve{z}\notag\\
&=\quad \sum_{l=1}^{L}  \ve{b}_l(\ve{\theta}) \frac{\partial \exdi{\ve{z};\ve{\theta}}{\phi_l(\ve{z})} )}{\partial \ve{\theta}}.
\end{align}
Supposing that the best choice for $\ve{\alpha}(\ve{\theta})$ is
\begin{align}
\ve{\alpha}(\ve{\theta})=\sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \exdi{\ve{z};\ve{\theta}}{  \phi_l(\ve{z}) },
\end{align}
it holds that
\begin{align}
&\exdi{} {\bigg( \sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \phi_l(\ve{z}) - \ve{\alpha}(\ve{\theta}) \bigg) \bigg( \sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \phi_l(\ve{z}) - \ve{\alpha}(\ve{\theta}) \bigg)^{\T} }\notag\\
&=\exdi{\ve{z};\ve{\theta}}{ \bigg( \sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \phi_l(\ve{z})\bigg) \bigg( \sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \phi_l(\ve{z}) \bigg)^{\T} }\notag\\&-\bigg(\sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \exdi{\ve{z};\ve{\theta}}{  \phi_l(\ve{z}) } \bigg) \bigg( \sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \exdi{\ve{z};\ve{\theta}}{ \phi_l(\ve{z})}\bigg)^{\T}.
\end{align}
Through the covariance inequality \cite{Lehman03}, we obtain
\begin{align}\label{proof:covariance:ineq:fish}
&\ve{F}(\ve{\theta}) \succeq
\Bigg(\sum_{l=1}^{L} \bigg( \frac{\partial  \exdi{\ve{z};\ve{\theta}}{  \phi_l(\ve{z}) }}{\partial \ve{\theta}} \bigg)^{\T}  \ve{b}^{\T}_l(\ve{\theta})  \Bigg)  \notag\\
&\cdot\Bigg(\exdi{\ve{z};\ve{\theta}}{ \bigg( \sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \phi_l(\ve{z})\bigg) \bigg( \sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \phi_l(\ve{z}) \bigg)^{\T} }\notag\\
&-\bigg(\sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \exdi{\ve{z};\ve{\theta}}{  \phi_l(\ve{z}) } \bigg) \bigg( \sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \exdi{\ve{z};\ve{\theta}}{ \phi_l(\ve{z})}\bigg)^{\T}\Bigg)^{-1}\notag\\
&\cdot \bigg(\sum_{l=1}^{L} \ve{b}_l(\ve{\theta}) \frac{\partial  \exdi{\ve{z};\ve{\theta}}{  \phi_l(\ve{z}) }}{\partial \ve{\theta}} \bigg).
\end{align}
\end{IEEEproof}
\section{Tightness of the Information Bound}\label{append:fish:bound:vector:tight}
\begin{IEEEproof}
Note, that with $\ve{\phi}(\ve{\theta})=\ve{t}(\ve{\theta})$ and optimized weightings $\ve{B}^{\star}(\ve{\theta})$, due to the definition \eqref{covariance:aux:stat}, we have
\begin{align}
\ve{B}^{\star\T}(\ve{\theta}) \ve{R}_{\ve{\phi}}(\ve{\theta}) \ve{B}^{\star}(\ve{\theta}) &=\exdi{\ve{z};\ve{\theta}}{ \ve{B}^{\star\T}(\ve{\theta})\ve{t}(\ve{\theta})\ve{t}^{\T}(\ve{\theta})\ve{B}^{\star}(\ve{\theta})}\notag\\ &-\ve{B}^{\star\T}(\ve{\theta})\ve{\mu}_{\ve{t}}(\ve{\theta}) \ve{\mu}^{\T}_{\ve{t}}(\ve{\theta})  \ve{B}^{\star}(\ve{\theta}).
\end{align}
Now let us assume that one possible optimizer is $\ve{B}^{\star}(\ve{\theta})=\frac{\partial \ve{w}(\ve{\theta}) }{\partial  \ve{\theta}} $. Then
\begin{align}
\ve{B}^{\star\T}(\ve{\theta}) \ve{R}_{\ve{\phi}}(\ve{\theta}) \ve{B}^{\star}(\ve{\theta}) &=\exdi{\ve{z};\ve{\theta}}{ \Big( \frac{\partial \ve{w}(\ve{\theta}) }{\partial  \ve{\theta}} \Big)^{\T} \ve{t}(\ve{\theta}) \ve{t}^{\T}(\ve{\theta}) \frac{\partial \ve{w}(\ve{\theta}) }{\partial  \ve{\theta}} }\notag\\
&-\Big(\frac{\partial \ve{w}(\ve{\theta}) }{\partial  \ve{\theta}} \Big)^{\T} \ve{\mu}_{\ve{t}}(\ve{\theta}) \ve{\mu}^{\T}_{\ve{t}}(\ve{\theta}) \frac{\partial \ve{w}(\ve{\theta}) }{\partial  \ve{\theta}}.
\end{align}
With
\begin{align}
\frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} &= \sum_{l=1}^{L} \frac{\partial w_l(\ve{\theta})}{\partial \ve{\theta}} t_l(\ve{z}) - \frac{\partial\lambda(\ve{\theta})}{\partial \ve{\theta}}\notag\\
&= \ve{t}^{\T}(\ve{z}) \frac{\partial \ve{w}(\ve{\theta}) }{\partial  \ve{\theta}} - \frac{\partial\lambda(\ve{\theta})}{\partial \ve{\theta}}
\end{align}
and
\begin{align}
\exdi{\ve{z};\ve{\theta}}{\frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}}} &= \ve{\mu}_{\ve{t}}^{\T}(\ve{\theta}) \frac{\partial \ve{w}(\ve{\theta}) }{\partial  \ve{\theta}}  - \frac{\partial\lambda(\ve{\theta})}{\partial \ve{\theta}}\notag\\
&=\ve{0}^{\T}
\end{align}
we obtain
\begin{align}\label{proof:tight:bound:vec}
&\ve{B}^{\star\T}(\ve{\theta}) \ve{R}_{\ve{\phi}}(\ve{\theta}) \ve{B}^{\star}(\ve{\theta}) =\notag\\
&=\ex{ \bigg( \frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} + \frac{\partial\lambda(\ve{\theta})}{\partial \ve{\theta}} \bigg)^{\T}  \bigg(\frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} + \frac{\partial\lambda(\ve{\theta})}{\partial \ve{\theta}}  \bigg)}\notag\\
&-\bigg( \frac{\partial\lambda(\ve{\theta})}{\partial \ve{\theta}} \bigg)^{\T} \bigg(  \frac{\partial\lambda(\ve{\theta})}{\partial \ve{\theta}} \bigg)\notag\\
&=\exdi{\ve{z};\ve{\theta}}{ \bigg( \frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}} \bigg)^{\T}  \bigg(\frac{\partial \ln p(\ve{z};\ve{\theta})}{\partial \ve{\theta}}  \bigg) }
\notag\\
&=\bigg( \frac{\partial \ve{\mu}_{\ve{t}}(\ve{\theta})}{\partial \ve{\theta}} \bigg)^{\T} \frac{\partial \ve{w}(\ve{\theta})}{\partial \ve{\theta}}.
\end{align}
Using \eqref{proof:tight:bound:vec} in \eqref{ineq:covar:fisher:compact} we obtain
\begin{align}
\ve{F}(\ve{\theta}) &\succeq \bigg( \frac{\partial \ve{\mu}_{\ve{t}}(\ve{\theta})}{\partial \ve{\theta}} \bigg)^{\T} \frac{\partial \ve{w}(\ve{\theta})}{\partial \ve{\theta}},
\end{align}
which holds with equality due to the identity \eqref{theorem:fish:identity:vec}.
\end{IEEEproof}
\section{CMLE Performance and Consistency}\label{append:amle:perf:unbias}
\begin{IEEEproof}
In order to analyze the performance of the CMLE, we proceed according to \cite[p. 212]{Kay93} and use a Taylor expansion of the used score function around the true parameter $\theta_t$
\begin{align}
\left.\frac{\partial \ln \tilde{p}(\ve{z};{\theta})}{\partial \theta} \right|_{\theta=\hat{\theta}} &= \left.\frac{\partial \ln \tilde{p}(\ve{z};{\theta})}{\partial \theta} \right|_{\theta={\theta}_t}\notag\\
&+\left.\frac{\partial^2 \ln \tilde{p}(\ve{z};{\theta})}{\partial \theta^2} \right|_{\theta={\theta}_t} (\hat{\theta} - {\theta}_t).
\end{align}
Due to the property \eqref{approx:score:zero}
\begin{align}
\left.\frac{\partial \ln \tilde{p}(\ve{z};{\theta})}{\partial \theta} \right|_{\theta={\theta}_t}=-\left.\frac{\partial^2 \ln \tilde{p}(\ve{z};{\theta})}{\partial \theta^2} \right|_{\theta={\theta}_t} (\hat{\theta} - {\theta}_t)
\end{align}
such that
\begin{align}\label{amle:performance}
\sqrt{N}(\hat{\theta} - {\theta}_t)=
\frac{ \frac{1}{\sqrt{N}} \left.\frac{\partial \ln \tilde{p}(\ve{z};{\theta})}{\partial \theta} \right|_{\theta={\theta}_t}}{-\frac{1}{N}\left.\frac{\partial^2 \ln \tilde{p}(\ve{z};{\theta})}{\partial \theta^2} \right|_{\theta={\theta}_t}}.
\end{align}
The denominator of \eqref{amle:performance}
\begin{align}
&-\frac{1}{N}\left.\frac{\partial^2 \ln \tilde{p}(\ve{z};{\theta})}{\partial \theta^2} \right|_{\theta={\theta}_t} =-\frac{1}{N} \sum_{n=1}^{N}\left.\frac{\partial^2 \ln \tilde{p}(z_n;{\theta})}{\partial \theta^2} \right|_{\theta={\theta}_t}
\end{align}
converges towards the constant value
\begin{align}
- \exdi{z;\theta}{\left.\frac{\partial^2 \ln \tilde{p}({z};{\theta})}{\partial \theta^2} \right|_{\theta={\theta}_t}}&= \sum_{l=1}^{L} b_l(\theta_t) \frac{ \partial \exdi{z;\theta_t}{\phi_l(z)} }{\partial \theta }\notag\\
&=\ve{b}^{\T}(\theta_t)\frac{\partial \ve{\mu}_{\ve{\phi}}(\theta_t)}{\partial \theta}
\end{align}
where we used the derivative of the replacement score
\begin{align}
\frac{\partial^2 \ln \tilde{p}(z;\theta)}{\partial \theta^2} &= \sum_{l=1}^{L} \frac{\partial b_l(\theta)}{\partial \theta} \big(\phi_l(z) - \exdi{z;\theta}{\phi_l(z)}\big)\notag\\
&- b_l(\theta) \frac{ \partial \exdi{z;\theta}{\phi_l(z)} }{\partial \theta }
\end{align}
and the property
\begin{align}\label{approx:score:derivative:zero}
\exdi{z;\theta}{\frac{\partial^2 \ln \tilde{p}(z;\theta)}{\partial \theta^2}} &=  -  \sum_{l=1}^{L} b_l(\theta) \frac{ \partial \exdi{z;\theta}{\phi_l(z)} }{\partial \theta }.
\end{align}
Due to the central limit theorem and the property
\begin{align}
\exdi{z;\theta}{\frac{\partial \ln \tilde{p}(z;\theta)}{\partial \theta}} = 0,
\end{align}
the nominator of \eqref{amle:performance}
\begin{align}
 \frac{1}{\sqrt{N}} \left.\frac{\partial \ln \tilde{p}(\ve{z};{\theta})}{\partial \theta} \right|_{\theta={\theta}_t} =  \frac{1}{\sqrt{N}}  \sum_{n=1}^{N} \left.\frac{\partial \ln \tilde{p}(z_n;{\theta})}{\partial \theta} \right|_{\theta={\theta}_t}
\end{align}
converges to a Gaussian random variable with zero mean
\begin{align}
&\exdi{\ve{z};\theta}{\frac{1}{\sqrt{N}}  \sum_{n=1}^{N} \left.\frac{\partial \ln \tilde{p}(z_n;{\theta})}{\partial \theta} \right|_{\theta={\theta}_t}}\notag\\
&=\sqrt{N} \exdi{z;\theta}{ \left.\frac{\partial \ln \tilde{p}(z;{\theta})}{\partial \theta} \right|_{\theta={\theta}_t}}=0
\end{align}
and variance
\begin{align}
&\exdi{\ve{z};\theta}{ \bigg( \frac{1}{\sqrt{N}}  \sum_{n=1}^{N} \left.\frac{\partial \ln \tilde{p}(z_n;{\theta})}{\partial \theta} \right|_{\theta={\theta}_t} \bigg)^2}\notag\\
&=\exdi{z;\theta}{ \bigg(\left.\frac{\partial \ln \tilde{p}(z;{\theta})}{\partial \theta}  \bigg)^2 \right|_{\theta={\theta}_t}}\notag\\
&=\exdi{z;\theta}{\bigg(\sum_{l=1}^{L} b_l(\theta) \big(\phi_l(z) - \exdi{z;\theta}{\phi_l(z)}\big)\bigg)^{2}}\notag\\
&=\exdi{z;\theta}{\bigg(\sum_{l=1}^{L} b_l(\theta) \phi_l(z) \bigg)^2 } - \bigg(\sum_{l=1}^{L} b_l(\theta) \exdi{z;\theta}{\phi_l(z)} \bigg)^2\notag\\
&=\ve{b}^{\T}(\theta) \ve{R}_{\ve{\phi}}(\theta) \ve{b}(\theta).
\end{align}
With Slutsky's theorem \cite[pp. 255]{Cram46}, it follows that asymptotically
\begin{align}
\sqrt{N}(\hat{\theta} - {\theta}_t)\sim\mathcal{N}\Bigg(0,\frac{\ve{b}^{\T}(\theta) \ve{R}_{\ve{\phi}}(\theta) \ve{b}(\theta)}{ \big(\frac{\partial \ve{\mu}_{\ve{\phi}}(\theta)}{\partial \theta}\big)^{\T} \ve{b}(\theta)  \ve{b}^{\T}(\theta) \frac{\partial \ve{\mu}_{\ve{\phi}}(\theta)}{\partial \theta}}\Bigg).
\end{align}
\end{IEEEproof}
\bibliographystyle{IEEEbib}