#  Estimation of Constrained Mean-Covariance of Normal Distributions

Anupam Kundu  Department of Statistics, Texas A&M University Mohsen Pourahmadi
Department of Statistics, Texas A&M University

111Author address: Department of Statistics, Texas A&M University, 3143 TAMU,
College Station, TX 77843-3143
Email<EMAIL_ADDRESS>(Corresponding Author<EMAIL_ADDRESS>
Abstract: Estimation of the mean vector and covariance matrix is of central
importance in the analysis of multivariate data. In the framework of
generalized linear models, usually the variances are certain functions of the
means with the normal distribution being an exception. We study some
implications of functional relationships between covariance and the mean by
focusing on the maximum likelihood and Bayesian estimation of the mean-
covariance under the joint constraint $\bm{\Sigma}\bm{\mu}=\bm{\mu}$ for a
multivariate normal distribution. A novel structured covariance is proposed
through reparameterization of the spectral decomposition of $\bm{\Sigma}$
involving its eigenvalues and $\bm{\mu}$. This is designed to address the
challenging issue of positive-definiteness and to reduce the number of
covariance parameters from quadratic to linear function of the dimension. We
propose a fast (noniterative) method for approximating the maximum likelihood
estimator by maximizing a lower bound for the profile likelihood function,
which is concave. We use normal and inverse gamma priors on the mean and
eigenvalues, and approximate the maximum aposteriori estimators by both MH
within Gibbs sampling and a faster iterative method. A simulation study shows
good performance of our estimators.

Abstract Keywords: Joint Mean-Covariance Estimation, Structured Covariance,
Shrinkage Inverse Wishart, Normal-Inverse Gamma, Metropolis-Hastings within
Gibbs Sampling.

##  1 Introduction

Mean and covariance estimation of multivariate normal distribution are
essential in almost every area of classical multivariate statistics (Bibby et
al., , 1979). The range of modern applications includes astrophysics
(Hamimeche and Lewis, , 2009), economics (Ledoit and Wolf, 2004a, ),
environmental sciences (Eguchi et al., , 2010), climatology (Guillot et al., ,
2015) and genetics (Schäfer and Strimmer, , 2005).

Estimation under joint constraints on the mean vector and covariance matrix of
data from a $N_{p}(\bm{\mu},\bm{\Sigma})$ distribution is relatively uncommon
in multivariate statistics (Bibby et al., , 1979). It is well-known that
estimation of a covariance matrix alone is a daunting task because of some
standard constraints e.g. 1) positive definiteness and 2) the number of
unknown covariance parameters growing quadratically with the dimension. Many
strategies are developed to bypass the notorious positive definiteness e.g.
spectral, Cholesky decomposition and factor models see Chiu et al.,  (1996);
Pourahmadi,  (1999); Fan et al.,  (2008). Here due to the nature of the
constraint on the mean and covariance, we rely on the spectral decomposition
and introduce a new class of structured covariance matrices with several
desirable properties. To the best of our knowledge this class of covariance
matrices seems to be new and has not been studied before.

We consider the following joint constraint on the mean vector and covariance
matrix of a multivariate normal distribution:

$\displaystyle\bm{\Sigma}\bm{\mu}$
$\displaystyle=\bm{\mu},\qquad|\bm{\Sigma}|=1,$ (1.1)

which appeared first in (Paine et al., , 2018) in the context of spherical
data. Interestingly, the first constraint forces the mean vector to be an
eigenvector of the covariance matrix corresponding to the eigenvalue one and
is more consequential. The second which constraints the product of the
remaining eigenvalues is less stringent and can be realized by rescaling, so
that without loss of generality it is ignored from here on.

Nevertheless, the joint constraint on the mean and covariance matrix will
definitely impact their estimators and the shape of the contours of the
multivariate normal density function as gleaned from the spectral
decomposition of the covariance matrix:

$\displaystyle\bm{\Sigma}=\bm{PDP}^{\top}$
$\displaystyle=\sum_{i=1}^{p}\lambda_{i}\bm{P}_{i}\bm{P}_{i}^{\top}=\sum_{i=1}^{p-1}\lambda_{i}\bm{P}_{i}\bm{P}_{i}^{\top}+\bm{\mu}\bm{\mu}^{\top}$
(1.2)

where $\bm{D}=\text{diag}(1,\lambda_{1},\lambda_{2},\dots,\lambda_{p-1})$ is
the diagonal matrix of ordered eigenvalues other than 1 and
$\bm{P}=[\bm{\mu},\bm{P}_{1},\bm{P}_{2},\dots,\bm{P}_{p-1}]$ is an orthogonal
matrix of eigenvectors.

Though these constraints arise in the context of directional data analysis
(Paine et al., , 2018), they seem to resonate with some deep classical issues
in statistical estimation theory. Note that the presence of the quadratic term
$\bm{\mu}\bm{\mu}^{\top}$ in (1.2) has resemblance with classical estimation
in the $N(\theta,\theta^{2})$ distribution. It also encourages modeling the
covariance matrix as a parsimonious quadratic function of the mean vector
similar to Hoff and Niu,  (2012).

We propose a simple and novel structured covariance model based on the
spectral decomposition (1.3) which overcomes several challenges in covariance
estimation. Let $\bm{u}=\frac{\bm{\mu}}{\|\bm{\mu}\|}$ be the direction of the
mean vector so that $\bm{\mu}=c_{0}\bm{u}$ where $c_{0}\in\mathbb{R}^{+}$ and
$\bm{u}$ lies on a unit sphere; $c_{0}$ can be interpreted as the radius of
the sphere on which the mean vector lie. Our structured covariance model is

$\displaystyle\bm{\Sigma}$
$\displaystyle=\bm{\Sigma}(\bm{u},\lambda_{1},\dots,\lambda_{p-1})=\bm{P}(\bm{u})\bm{D}\bm{P}^{\top}(\bm{u})$
(1.3)

where $\bm{D}=\text{ diag }(1,\lambda_{1},\dots,\lambda_{p-1})=\text{ diag
}(1,\bm{\lambda})$ is the matrix of eigenvalues, and for a given value of the
mean direction $\bm{u}$, the orthogonal matrix of eigenvectors is
$\bm{P}(\bm{u})=[\bm{u},\bm{V}(\bm{u})]$ where
$\bm{V}=\bm{V}(\bm{u})\in{R}^{p\times(p-1)}$ is a known matrix function so
that orthogonality of $\bm{P}$ is ensured. A simple and prominent example of
such a $\bm{V}$ is obtained by an application of Gram-Schmidt procedure
(Trefethen and Bau III, , 1997) to the set
$\\{\bm{u},\bm{e}_{1},\bm{e}_{2},\dots,\bm{e}_{p-1}\\}$ where
$\bm{e}_{i},i=1,2,\dots,p-1$ denote the canonical basis of $\mathcal{R}^{p}$
and $u_{p}\neq 0$. With $\bm{V}$ assumed known in (1.3), the unknown
parameters then are $(\bm{\mu},\bm{\lambda})$, the vectors of mean and the
$(p-1)$ eigenvalues, so that the number of parameters drastically reduces from
quadratic to linear i.e. $(2p-1)$ in the dimension. Moreover, the first
constraint is automatically satisfied, since using
$\bm{\mu}\perp\bm{V}_{i}(\bm{u})=\bm{V}_{i}$ (columns of the $\bm{V}$ matrix)
for $i=1,2,\dots,(p-1)$, it follows that

$\displaystyle\bm{\Sigma}\bm{\mu}$
$\displaystyle=\bm{P}(\bm{u})\bm{\Lambda}\bm{P}^{\top}(\bm{u})\bm{\mu}=\sum_{i=1}^{p-1}\lambda_{i}\bm{V}_{i}\bm{V}^{\top}_{i}\bm{\mu}+\bm{u}\bm{u}^{\top}\bm{\mu}=c_{0}\bm{u}\bm{u}^{\top}\bm{u}=\bm{\mu}.$
(1.4)

Our approximate MLE of the parameters in Section 2 provides positive estimates
for the eigenvalues and hence guarantees the positive-definiteness of the
estimated structured covariance matrix.

Next, we will provide two concrete examples to elucidate the idea behind the
model and its components, particularly the eigenvectors as functions of the
mean vector.

Example 1: Consider the case of equal means $\bm{\mu}=c\mathbb{1}$, then using
Gram-Schmidt procedure (Trefethen and Bau III, , 1997), it follows that
$\bm{P}(\bm{u})$ is of the form:

$\displaystyle\bm{P}(\bm{u})$
$\displaystyle=\begin{bmatrix}\mathbb{1}_{p}&\bm{z}^{(p)}_{p}&\bm{z}^{(p)}_{p-1}&\bm{z}^{(p)}_{p-2}&\dots&\bm{z}^{(p)}_{2}\end{bmatrix}$
(1.5)

where

$\displaystyle\bm{z}^{(p)}_{s}$
$\displaystyle=\left(0,\dots,0,\bm{z}^{(s)}_{s}\right)^{\top}\quad\quad\quad\forall
s\in\left\\{2,3\dots,p\right\\}$ $\displaystyle\bm{z}^{(s)}_{s}$
$\displaystyle=\left(\frac{s-1}{\sqrt{s(s-1)}},\frac{-1}{\sqrt{s(s-1)}},\dots,\frac{-1}{\sqrt{s(s-1)}}\right)^{\top}$
(1.6)

and
$\mathbb{1}_{p}=\left(\frac{1}{\sqrt{p}},\frac{1}{\sqrt{p}},\dots,\frac{1}{\sqrt{p}}\right)$.
The $c$ drops due to normalization.

Example 2: As a generalization of Example 1, take

$\bm{\mu}=(\mu_{1},\mu_{2},\mu_{3},\dots,\mu_{3})^{\top}$

where the first three entries of the mean vector are different and rest of
them are the same as the third entry of the vector. Then, using the
construction method described in the paragraph following (1.3), the $\bm{P}$
matrix takes the form:

$\displaystyle\bm{P}(\bm{u})=\begin{bmatrix}\bm{u}&\bm{w}_{1}&\bm{w}_{2}&z_{p-2}^{(p)}&\dots
z_{2}^{(p)}\end{bmatrix}$ (1.7)

where $z_{s}^{(p)}$’s are same as in Example 1, and

$\displaystyle\bm{w}_{1}$
$\displaystyle=\left(\mu_{0}^{2},-\mu_{1}\mu_{2},-\mu_{1}\mu_{3},\dots,-\mu_{1}\mu_{3}\right)^{\top}/\left(\mu_{0}\|\mu\|\right)$
$\displaystyle\bm{w}_{2}$
$\displaystyle=\left(0,(p-2)\mu_{3},-\mu_{2},\dots,-\mu_{2}\right)^{\top}/\left((p-2)\mu_{0}\right),$
(1.8)

where $\mu_{0}^{2}=\mu_{2}^{2}+(p-2)\mu_{3}^{2}$. When $p=3$, $\bm{w}_{1}$ and
$\bm{w}_{2}$ reduce to the $\tilde{\xi}_{1}$ and $\tilde{\xi}_{2}$ in Paine et
al.,  (2018, §2.3), so that this example is a generalization of their
construction to the case of $p>3$ or more general spherical data.

Spherical data arise in many scientific disciplines like shape analysis,
geology, meteorology (e.g. Mardia and Jupp,  (2000)), text analysis (e.g.
Hamsici and Martinez,  (2007)), etc. Paine et al.,  (2020) study spherical
regression using their elliptically symmetric angular Gaussian distribution.
Since our constraint is similar without the spherical nature of the
distribution, we can always project them in the corresponding $\mathbb{R}^{p}$
(Mardia and Jupp, , 2000, §3.5.6) and fit a normal model with the constraint
(1.1). This is helpful in developing multivariate techniques like discriminant
analysis, clustering, etc. for spherical data. We have applied it on earth’s
historic magnetic pole data by Schmidt,  (1976).

Our goal is to compute the maximum likelihood estimator and provide a Bayesian
framework for estimation of the parameters in (1.3). It turns out that
computing the maximum likelihood estimator is challenging due to the
nonlinearity and intractability of $\bm{P}(\bm{u})$. However, the MLE of the
eigenvalues and $c_{0}$ have closed forms as a function of the mean direction,
thus making it possible to compute its profile likelihood. At this stage, we
approximate the MLE of the mean direction by first finding a lower bound for
the concave profile log-likelihood function of the mean direction (for given
$\bm{\lambda}$) and then maximizing it. In the Bayesian context, the maximum
aposteriori (MAP) does not have a closed form, even though the posterior
distribution obtained by using a Gaussian prior for the mean vector and
inverse gamma on the eigenvalues does. These priors are suggested by extending
the shrinkage inverse Wishart prior of Berger et al.,  (2020) to the case of
nonzero means. We show in Section 3.4 that one can generate from the posterior
distribution quite easily using Metropolis-Hastings within Gibbs sampling. But
if we are only interested in point estimation (e.g. MAP), then we propose a
simpler way of computing it. The idea is to follow the calculation of the
approximate MLE by providing a lower bound for the posterior and maximize it
using a modified version of Newton’s method.

The paper is organized as follows: Section 2 delves into the detailed
computation of an approximate maximum likelihood estimator. Section 3
discusses the standard priors such as normal-inverse Wishart and normal-
shrinkage inverse Wishart, motivates the normal-inverse gamma prior selection,
provide details of the MH within Gibbs sampling process for generation from
the posterior and imitates the approximate MLE calculation to propose a fast
method to compute the Bayesian point estimate. Section 4 computes the risks of
these methods of approximation through simulation and exhibits the performance
of our method. Finally in section 5, we discuss the key take away for this
paper.

##  2 MLE of $\mu$ and $\lambda$ in Model (1.3)

Let $\bm{x}_{1},\bm{x}_{2},\dots,\bm{x}_{n}$ be a sample of size $n$ from
$N_{p}(\bm{\mu},\bm{\Sigma})$, where $\bm{\Sigma}$ is parameterized as in
(1.3) with $\bm{\mu}=c_{0}\bm{u}$ and $\bm{X}$ is the $n\times p$ data matrix.
Then, the log-likelihood can be written as

$\displaystyle l(\bm{u},c_{0},\bm{\lambda}\mid\bm{X})$
$\displaystyle\propto-\frac{n}{2}\sum_{i=1}^{p-1}\log(\lambda_{i})-\frac{1}{2}\left[\sum_{j=1}^{n}(c_{0}\bm{u}-\bm{x}_{i})^{\top}\bm{\Sigma}^{-1}(c_{0}\bm{u}-\bm{x}_{i})\right]$
(2.1)
$\displaystyle\propto-\frac{n}{2}\sum_{i=1}^{p-1}\log(\lambda_{i})-\frac{1}{2}\text{Tr}\left[\bm{D}^{-1}\bm{B}(\bm{X},c_{0},\bm{u})\right]$
$\displaystyle\propto-\frac{n}{2}\sum_{i=1}^{p-1}\log(\lambda_{i})-\frac{1}{2}\left[\bm{B}(\bm{X},c_{0},\bm{u})_{11}+\sum_{i=1}^{p-1}\frac{\bm{B}(\bm{X},c_{0},\bm{u})_{(i+1)(i+1)}}{\lambda_{i}}\right],$
(2.2)

where
$\bm{B}(\bm{X},c_{0},\bm{u})=\bm{P}^{\top}(\bm{u})\bm{A}(c_{0}\bm{u})\bm{P}(\bm{u})$
and
$\bm{A}(\bm{\mu})=\sum_{i=1}^{n}(\bm{x}_{i}-\bm{\mu})(\bm{x}_{i}-\bm{\mu})^{\top}$.

Evidently, maximization of the log-likelihood function with respect to the
mean direction is more challenging than those with respect to $c_{0}$ and
$\bm{\lambda}$. This suggests and we propose approximating the MLE of the
parameters using the following three (non-iterative) steps:

  1. 1.

Differentiate the log-likelihood function to obtain the MLE of $c_{0}$ and
$\lambda_{i}$’s as a function of the mean direction vector $\bm{u}$.

  2. 2.

Compute the profile log-likelihood.

  3. 3.

Find a lower bound for the profile log-likelihood and maximize it to obtain
the approximate MLE of the mean direction $\bm{u}$. Use the latter to compute
the approximate MLE of $c_{0}$ and the eigenvalues.

Next, we provide further details about implementing these three steps.

Step 1: For a given $\bm{u}$, differentiating the log-likelihood in (2.1) with
respect to $c_{0}$ (Bibby et al., , 1979, §4.2.9) and in (2.2) with respect to
$\lambda_{i}$, we obtain

$\displaystyle\widehat{c}_{0}=\frac{\bm{u}^{\top}\bm{P}(\bm{u})\bm{D}^{-1}\bm{P}^{\top}(\bm{u})\bar{\bm{x}}}{\bm{u}^{\top}\bm{P}(\bm{u})\bm{D}^{-1}\bm{P}^{\top}(\bm{u})\bm{u}},\quad$
$\displaystyle\text{ and
}\quad\widehat{\lambda}_{i}=\frac{\bm{B}(\bm{X},c_{0},\bm{u})_{(i+1)(i+1)}}{n}.$
(2.3)

The expression for the $\widehat{c}_{0}$ in (2.3) reduces to

$\displaystyle\widehat{c}_{0}$
$\displaystyle=\frac{\bm{e}^{\top}_{1}\bm{D}^{-1}\bm{P}^{\top}(\bm{u})\bar{\bm{x}}}{\bm{e}_{1}^{\top}\bm{D}^{-1}\bm{e}_{1}}=\frac{\frac{1}{\lambda_{1}}\bm{e}^{\top}_{1}\bm{P}^{\top}(\bm{u})\bar{\bm{x}}}{\frac{1}{\lambda_{1}}}=\bm{u}^{\top}\bar{\bm{x}},$
(2.4)

and note that the estimates of the eigenvalues are always positive.

Step 2: Substituting for $\widehat{c}_{0}$ and $\widehat{\lambda}_{i}$ in the
log-likelihood (2.2), the profile likelihood of the mean direction $\bm{u}$
turns out to be

$\displaystyle l(\bm{u},\widehat{c}_{0},\widehat{\bm{\lambda}}\mid\bm{X})$
$\displaystyle\propto-\frac{n}{2}\sum_{i=1}^{p-1}\log\left(\frac{\bm{B}(\bm{X},\widehat{c}_{0},\bm{u})_{(i+1)(i+1)}}{n}\right)-\frac{1}{2}\left[\bm{B}(\bm{X},\widehat{c}_{0},\bm{u})_{11}+n(p-1)\right]$
(2.5)

Denoting the columns of $\bm{V}$ by $\bm{V}_{i}$ for $i=1,2,\dots,(p-1)$, then
the diagonal entries of the $\bm{B}(\bm{X},\widehat{c}_{0},\bm{u})$ are

$\displaystyle\left\\{\bm{B}(\bm{X},\widehat{c}_{0},\bm{u})_{11},\bm{B}(\bm{X},\widehat{c}_{0},\bm{u})_{22},\dots,\bm{B}(\bm{X},\widehat{c}_{0},\bm{u})_{p,p}\right\\}$
$\displaystyle=\left\\{\bm{u}^{\top}\bm{A}(\widehat{c}_{0}\bm{u})\bm{u},\bm{V}^{\top}_{1}\bm{A}(\widehat{c}_{0}\bm{u})\bm{V}_{1},\dots,\bm{V}^{\top}_{p-1}\bm{A}(\widehat{c}_{0}\bm{u})\bm{V}_{p-1}\right\\}.$

Since $\bm{u}\perp\bm{V}_{i}$ for $i=1,2,\dots,{p-1}$, one may further
simplify the diagonal entries of $\bm{B}(\bm{X},\widehat{c}_{0},\bm{u})$ to
the following:

$\displaystyle\bm{V}_{i}^{\top}\bm{A}(\widehat{c}_{0}\bm{u})\bm{V}_{i}$
$\displaystyle=\bm{V}_{i}^{\top}\left[\bm{A}(\bm{0})-n\widehat{c}_{0}\bar{\bm{x}}\bm{u}^{\top}-n\widehat{c}_{0}\bm{u}\bar{\bm{x}}^{\top}+n\widehat{c}_{0}^{2}\bm{u}\bm{u}^{\top}\right]\bm{V}_{i}=\bm{V}_{i}^{\top}\bm{A}(\bm{0})\bm{V}_{i}$
$\displaystyle\bm{u}^{\top}\bm{A}(\widehat{c}_{0}\bm{u})\bm{u}$
$\displaystyle=\bm{u}^{\top}\bm{A}(\bar{\bm{x}})\bm{u}.$ (2.6)

In spite of this simplification, the profile log-likelihood is hard to
differentiate as a function of $\bm{u}$, in general, when we are not assuming
any specific form for the matrix $\bm{V}$.

Step 3: We find a workable lower bound for (2.5) and maximize it with respect
to the mean direction to obtain an approximate MLE for $\bm{u}$.
Alternatively, this amounts to assuming that the first sum in the profile
likelihood (2.7) is constant. The expression in equation (2) is a quadratic
form of the matrix $\bm{A(0)}=\sum_{j=1}^{n}\bm{x}_{j}\bm{x}_{j}^{\top}$ and
can be bounded by its largest eigenvalue (Rao, , 1973, §1f.2.1). Consequently,
using (2) and that columns of $\bm{V}$ are orthonormal simplify the profile
log-likelihood (2.5) and we arrive at the following lower bound:

$\displaystyle
l(\bm{u},\widehat{c}_{0},\bm{\lambda}\mid\bm{X})\propto-\frac{n}{2}\sum_{i=1}^{p-1}\log\left[\frac{\bm{V}_{i}^{\top}\bm{A}(\bm{0})\bm{V}_{i}}{n}\right]-\frac{1}{2}\left[\bm{u}^{\top}\bm{A}(\bar{\bm{x}})\bm{u}+n(p-1)\right]$
$\displaystyle\geq-\frac{n(p-1)}{2}\log\left[\frac{\lambda_{1}\left\\{\bm{A}(\bm{0})\right\\}}{n}\right]-\frac{1}{2}\left[\bm{u}^{\top}\bm{A}(\bar{\bm{x}})\bm{u}+n(p-1)\right]=h(\bm{u})$
(2.7)

Evidently, the lower bound denoted by $h(\bm{u})$ is maximized by minimizing
the quadratic expression inside the second bracket. Since $\bm{u}\neq 0$ and
$h$ is a quadratic function of $\bm{u}$, the maximum occurs when $\bm{u}$ is
the eigenvector corresponding to the smallest eigenvalue of the matrix
$\bm{A}(\bar{\bm{x}})$. As such the approximate MLE of $\bm{u}$ is unique only
up to a sign.

Even though the computed mean direction $\widehat{\bm{u}}$ is not the exact
MLE, it is still a good approximation as confirmed by the simulation results
in Table LABEL:tab:mlemapNR of Section 4. The approximate maximum likelihood
estimate of the eigenvalues and the constant $c_{0}$ are obtained by plugging
in the approximate MLE of $\bm{u}$ in (2.3) and (2.4). Fortunately, this idea
of approximatiing MLE can also be replicated in our posterior MAP estimator
approximation developed in the next section, but the maximization of the lower
bound is not as straightforward and requires employing a version of Newton’s
iterative method.

##  3 Bayesian Estimation of $(\mu,\lambda)$ in (1.3)

In this section we develop a Bayesian methodology for estimating the
parameters using two sets of priors reviewed in the next two subsections.

###  3.1 Normal-Inverse Wishart Priors:

For the parameters of a multivariate normal distribution the most popular
prior is normal-inverse Wishart described in Gelman et al.,  (2013, Section
3.6) and Barnard et al.,  (2000) with the hyperparameters
$(\bm{\mu}_{0},\kappa_{0};\nu_{0},\bm{\Lambda}_{0})$:

$\displaystyle\bm{\mu}\mid\bm{\Sigma}\sim
N_{p}(\bm{\mu}_{0},\bm{\Sigma}/\kappa_{0}),\qquad\bm{\Sigma}\sim\pi_{IW}(\nu_{0},\bm{\Lambda}^{-1}_{0}),$

where $\nu_{0}$ and $\bm{\Lambda}_{0}$ are the degrees of freedom and the
scale matrix. The suggested values for the hyperparameters in Gelman and Hill,
(2006) are $\bm{\Lambda}_{0}=\bm{I}_{p}$ and $\nu_{0}=p+1$. The remaining
hyperparameters are the prior mean, $\bm{\mu}_{0}$ and $\kappa_{0}$ on the
$\bm{\Sigma}$ scale. Due to conjugacy, the posterior density has the following
parameters:

$\displaystyle\bm{\mu}_{n}$
$\displaystyle=\frac{\kappa_{0}}{\kappa_{0}+n}\bm{\mu}_{0}+\frac{n}{\kappa_{0}+n}\bar{\bm{x}},$
$\displaystyle\bm{\Lambda}_{n}$
$\displaystyle=\bm{\Lambda}_{0}+\bm{A}+\frac{n\kappa_{0}}{\kappa_{0}+n}(\bar{\bm{x}}-\bm{\mu}_{0})(\bar{\bm{x}}-\bm{\mu}_{0})^{\top}.$
(3.1)

As a point estimator posterior mode or the maximum aposteriori probability
(MAP) estimator for $\bm{\Sigma}$ is generally the most popular choice
(Murphy, , 2012, §5.2.1). The posterior marginal distribution of the mean
parameter $\bm{\mu}$ is multivariate
$t_{\nu_{n}-p+1}\left(\bm{\mu}_{n},\frac{\bm{\Lambda}_{n}}{\kappa_{n}(\nu_{n}-p+1)}\right)$
with the MAP estimator $\bm{\mu}_{n}$, and the MAP estimator for $\bm{\Sigma}$
is given by

$\displaystyle\widehat{\bm{\Sigma}}_{map}$
$\displaystyle=\frac{\bm{\Lambda}_{n}}{\nu_{n}+p+2},$

where $\nu_{n}=\nu_{0}+n$ (proof is shown in appendix 1).

###  3.2 Shrinkage Inverse Wishart Priors

The most natural prior for mean and covariance i.e. normal-inverse Wishart
reviewed above is known to overdisperse the eigenvalues of the posterior
covariance estimator. In fact, Yang and Berger,  (1994) describes that the
normal-inverse Wishart prior has the term
$\prod_{i<j}(\lambda_{i}-\lambda_{j})$ in the density where $\lambda_{i}$’s
are the ordered eigenvalues of $\bm{\Sigma}$, thus forcing the eigenvalues
apart and increasing the variability (Berger et al., , 2020). This is one
major motivation for using the normal-shrinkage inverse Wishart prior.

The shrinkage-inverse Wishart (SIW) prior has the same density as inverse
Wishart prior save the extra term $\prod_{i<j}(\lambda_{i}-\lambda_{j})$ in
the denominator:

$\displaystyle\pi_{SIW}(\bm{\Sigma}\mid\nu_{0},b,\bm{\Lambda}^{-1}_{0})$
$\displaystyle\propto\mid\bm{\Sigma}\mid^{-\frac{\nu_{0}+p+1}{2}}\frac{\exp{\left[-\frac{1}{2}Tr(\bm{\Lambda}_{0}\bm{\Sigma}^{-1})\right]}}{\prod_{i<j}(\lambda_{i}-\lambda_{j})^{b}},$
(3.2)

where $\nu_{0}$ is a real constant, $b\in[0,1]$ and $\bm{\Lambda}_{0}$ is a
positive semi-definite matrix. It is interesting to notice that $b=0$
corresponds to some common priors like inverse Wishart, reference, Jefferey’s
etc. which also contain the term $\prod_{i<j}(\lambda_{i}-\lambda_{j})$ , see
Berger et al.,  (2020). The posterior density for $(\bm{D},\bm{P})$,

using the one-to-one transformation from $\bm{\Sigma}$ to
$\bm{D}=\text{diag}(\lambda_{1},\dots,\lambda_{p})$ and the orthogonal
eigenvector matrix $\bm{P}$, turns out to be

$\displaystyle\pi_{SIW}(\bm{D},\bm{P}\mid\nu_{0},b,\bm{\Lambda}^{-1}_{0})$
$\displaystyle\propto\mid\bm{\Sigma}\mid^{-\frac{\nu_{0}+p+1}{2}}\frac{\exp{\left[-\frac{1}{2}Tr(\bm{\Lambda}_{0}\bm{P}\bm{D}^{-1}\bm{P}^{\top})\right]}}{\prod_{i<j}(\lambda_{i}-\lambda_{j})^{(b-1)}}\mathbb{1}_{\lambda_{1}\geq\lambda_{2}\geq\dots\lambda_{p}}.$
(3.3)

Note that the posterior is zero whenever the eigenvalues are close together so
that effectively it forces the eigenvalues apart. Moreover, when $b=1$, the
term in question goes away retaining the conjugacy property, even with the
additional normal prior on the nonzero mean (Proof is shown in the Appendix 2)
and the same posterior parameters as before. The next important consequence of
choosing $b=1$, from (3.3), is that the conditional distribution of the
eigenvalues given the orthogonal matrix $\bm{P}$ is an ordered inverse gamma
distribution (Berger et al., , 2020, §3.2). This plays a central role in
developing a computational scheme for computing the MAP.

Recall that the structured covariance model in (1.3) has its orthogonal matrix
$\bm{P}$ determined by a non-zero mean vector. This requires working with
nonzero mean vectors in contrast to the mean zero framework in Berger et al.,
(2020), while maintaining the conditional distribution of the eigenvalues
given the mean vector as an inverse gamma distribution, see also Berger et
al.,  (2020). This motivates us to select an inverse gamma prior for the
eigenvalues as in (3.4) which, fortunately, in turn implies that the
conditional distribution of eigenvalues is the inverse gamma (see section
3.4).

###  3.3 The Mean-Eigenvalue Priors

A convenient prior distribution on the mean vector is the multivariate normal
and the inverse gamma on the ordered eigenvalues. The latter prior comes
naturally from the shrinkage inverse Wishart prior of (Berger et al., , 2020,
§3.2) as the conditional distribution of eigenvalues given the matrix of
eigenvectors (see section 3.2 for detailed discussion), see also Hoff,  (2009)
(section 3.3). More specifically, our proposed prior is:

$\displaystyle\bm{\mu}\mid\bm{D}\sim N_{p}(\bm{\mu}_{0},\bm{D}/\kappa_{0})$
$\displaystyle\text{ and }\bm{\lambda}_{i}\sim\text{Inverse-
gamma}(a-1,c_{i}/2),$ (3.4)

where $a,c_{i}$’s are the hyperparameters. Setting
$\bm{H}_{0}=\text{diag}(1,c_{1},c_{2},\dots,c_{p-1})$, then the posterior
distribution has the form

$\displaystyle p(\bm{\mu},\bm{\lambda}\mid\bm{X})$
$\displaystyle\propto\left(\prod_{i=1}^{p-1}\lambda_{i}\right)^{-\frac{n+1+2a}{2}}\exp\left[-\frac{1}{2}Tr\left\\{\bm{D}^{-1}\bm{H}_{N}\right\\}\right]$
(3.5)

with
$\bm{H}_{N}=\bm{B}(\bm{X},\bm{\mu})+\kappa_{0}(\bm{\mu}-\bm{\mu}_{0})(\bm{\mu}-\bm{\mu}_{0})^{\top}+\bm{H}_{0}$.
Here $\bm{B}(\bm{X},\bm{\mu})$ has the same expression as
$\bm{B}(\bm{X},c_{0},\bm{u})$ which appeared in Section 2.

One needs to generate from this posterior distribution when computing the
Bayesian point estimates, credible intervals and various other quantities. We
show that Metropolis-Hastings within Gibbs sampling is possible here and
suitable for our goals (Brooks et al., , 2011, §1.12.10).

###  3.4 Gibbs Sampling

The Gibbs sampler (Gelfand and Smith, , 1990) method is a numerical technique
for sampling from the joint posterior distribution. Given an initial vector,
the Gibbs sampling proceeds by sampling from the conditional posterior
distribution. More generally, if the full conditional posterior distribution
in any Gibbs step is of a non-standard form, using MH (Metropolis-Hastings)
step is convenient (Brooks et al., , 2011). The technique is known as
Metropolis-Hastings within Gibbs Sampling or alternatively as single-component
Metropolis-Hastings sampling as follows:

  1. 1.

Since $p(\bm{D}\mid\bm{\mu},\bm{X})\propto p(\bm{\mu},\bm{\lambda}\mid\bm{X})$
has an Inverse-Gamma distribution, the orderd eigenvalues are generated from
independent Inverse Gamma $\left(\frac{n+2a-1}{2},c^{*}_{i}/2\right)$ where
$c^{*}_{i}$ is the i-th diagonal entry of $\bm{H}_{N}$ matrix.

  2. 2.

Generate from $p(\bm{\mu}\mid\bm{D},\bm{X})$ using the MH algorithm with a
proposal distribution $q$ selected as

$\displaystyle
q\left(.\mid\bm{\mu}^{(i-1)}\right)=N_{p}\left(\bm{\mu}^{(i-1)},\frac{1}{n}\bm{P}\left(\bm{\mu}^{(i-1}\right)\bm{D}\bm{P}^{\top}\left(\bm{\mu}^{(i-1)}\right)\right).$

More concretely, the steps of MH within Gibbs algorithm in our context are as
follows:

Algorithm 1 MH within Gibbs for generating samples from the posterior
distribution

1:Start with $\bm{\mu}^{(0)}=\bar{\bm{x}}$.

2:Repeat $s$ times: $j-th$ step

3: Generate $\bm{\lambda}^{(j)}_{i}\sim
IG\left(\frac{n+2a-1}{2},c^{*}_{i}/2\right)$ to form $\bm{D}^{(j)}$.

4: Start with $\bm{\mu}^{(j-1)}$, Repeat MH Step $l$ times: $k-th$ step

5: Generate $\bm{\mu}^{*}\sim q\left(.\mid\bm{\mu}^{(k-1)}\right)$

6: Calculate
$r\left(\bm{\mu}^{*},\bm{\mu}^{(k-1)}\right)=\min\left\\{1,\frac{p\left(\bm{\mu}^{*},\bm{\lambda}^{(j)}\mid\bm{X}\right)q\left(\bm{\mu}^{(k-1)}\mid\bm{\mu}^{*}\right)}{q\left(\bm{\mu}^{*}\mid\bm{\mu}^{(k-1)}\right)p\left(\bm{\mu}^{(k-1)},\bm{\lambda}^{(j)}\mid\bm{X}\right)}\right\\}$

7: Generate $u\sim U(0,1)$

8: Set $\bm{\mu}^{(k)}=\bm{\mu}^{*}$ if
$u<r\left(\bm{\mu}^{*},\bm{\mu}^{(k-1)}\right)$ else
$\bm{\mu}^{(k)}=\bm{\mu}^{(k-1)}$

9: Set $\bm{\mu}^{(j)}=\bm{\mu}^{(l)}$

10:Collect $(\bm{\mu}^{(j)},\bm{\lambda}^{(j)})$ for $j=1,2,\dots,s$

11:End

An interesting observation is that the initial value of $\bm{\lambda}$ is
irrelevant. In the MCMC, the chain for the eigenvalues depend on the previous
iteration only through the generation of mean vector. This is happening
because the parameters for the conditional distribution of the eigenvalues are
$(n+2a-1)/2$ and $c^{*}_{i}/2$ where $c^{*}_{i}$’s are the diagonal entries of
the $\bm{H}_{N}$ matrix and hence a function of the mean vector.

To illustrate the algorithm, we have generated 100 pairs of
$(\bm{\mu},\bm{\lambda})$ from the posterior distribution by MH within Gibbs
sampling. We find out the posterior maximizing $(\bm{\mu},\bm{\lambda})$ pair
by exhaustive computation of posterior density at each pair. This is possible
because the posterior density has a closed form barring the normalizing
constants. We take the value of the mean vector from the posterior maximizing
pair, set it as our MAP estimate for the mean vector discarding the
corresponding $\bm{\lambda}$ and make a final update on $\bm{\lambda}$. Given
the estimate of the mean vector, the posterior conditional distribution of
$\bm{\lambda}$ is inverse gamma and has a known mode. Using this, we can
calculate the best eigenvalue vector $\widehat{\bm{\lambda}}$ to be
$\widehat{\lambda}_{i}=\frac{c_{i}^{*}}{n+1+2a}$ for $i=1,2,\dots,(p-1)$. The
newly calculated value $\widehat{\lambda}$, accompanied by the MAP estimate of
$\bm{\mu}$ constitute our MAP estimator
$\left(\widehat{\mu},\widehat{\lambda}\right)$. The performance of such an
estimator is assessed through a simulation study (see Section 4.2.1) and Table
LABEL:tab:Gibbstable of the simulation Section 4.

###  3.5 Approximation of MAP through a Lower Bound

The Gibbs sampling is computationally challenging and the time complexity
increases exponentially with the dimensions. We follow the steps Section 2 to
approximate the MAP estimator using a lower bound for the log posterior
density. The main difference is in the maximization step of the lower bound as
a function of $\bm{u}$ where it is not as straightforward as the MLE case and
does not have a closed-form. So, we resort to an iterative modified Newton-
Raphson algorithm. Using $\bm{\mu}=c_{0}\bm{u}$ where $c_{0}\in\mathbb{R}$ as
before we arrive at

$\displaystyle\log p(\bm{u},c_{0},\bm{\lambda}\mid\bm{X})$
$\displaystyle\propto-t\sum_{i=1}^{p-1}\log\lambda_{i}-\frac{1}{2}\left[f(c_{0})+\text{
Tr }(\bm{D}^{-1}\bm{H}_{0})\right]$ (3.6)
$\displaystyle\propto-t\sum_{i=1}^{p-1}\log\lambda_{i}-\frac{1}{2}\left[\left(\bm{H}_{N}\right)_{11}+\sum_{i=1}^{p-1}\frac{\left(\bm{H}_{N}\right)_{i+1,i+1}}{\lambda_{i}}\right]$
(3.7)

where $t=\frac{n+1+2a}{2}$ and
$f(c_{0})=\sum_{i=1}^{n}(\bm{x}_{i}-c_{0}\bm{u})^{\top}\bm{\Sigma}^{-1}(\bm{x}_{i}-c_{0}\bm{u})+\kappa_{0}(c_{0}\bm{u}-\bm{\mu}_{0})^{\top}\bm{D}^{-1}(c_{0}\bm{u}-\bm{\mu}_{0}).$
Following similar calculations as in (2.3), the estimates of $c_{0}$ and the
eigenvalues in terms of the mean direction are as follows:

$\displaystyle\widehat{c}_{0}=\frac{n\bm{u}^{\top}\bar{\bm{x}}+\kappa_{0}\bm{u}^{\top}\bm{D}^{-1}\bm{\mu}_{0}}{n+\kappa_{0}\bm{u}^{\top}\bm{D}^{-1}\bm{u}},\quad$
$\displaystyle\quad\widehat{\lambda}_{i}=\frac{\left(\bm{H}_{N}\right)_{i+1,i+1}}{2t}$
(3.8)

where
$\bm{H}_{N}=\bm{B}(\bm{X},c_{0},\bm{u})+\kappa_{0}(c_{0}\bm{u}-\bm{\mu}_{0})^{\top}\bm{D}^{-1}(c_{0}\bm{u}-\bm{\mu}_{0})+\bm{H}_{0}$.
Here the expression of $\widehat{c}_{0}$ and $\widehat{\lambda}_{i}$ are
intertwined. So we are doing the following calculation with the aim to propose
an iterative algorithm for the MAP estimators. We will use equation (3.8) and
(3.15) as our updating equation.

Using the expression of $\widehat{\lambda}_{i}$ from equation (3.8) in the log
posterior likelihood function we obtain the following analogue of the profile
likelihood in (2.5):

$\displaystyle\log p(\bm{u},c_{0},\widehat{\bm{\lambda}}\mid\bm{X})$
$\displaystyle\propto-t\sum_{i=1}^{p-1}\log\left(\frac{\left(\bm{H}_{N}\right)_{i+1,i+1}}{2t}\right)-\frac{1}{2}\left[\left(\bm{H}_{N}\right)_{11}+(p-1)t\right]$
(3.9)

where the diagonal entries of $\bm{H}_{N}$ are

$\displaystyle\left(\bm{H}_{N}\right)_{i+1,i+1}$
$\displaystyle=\begin{cases}\bm{u}^{\top}\bm{A}(c_{0}\bm{u})\bm{u}+\kappa_{0}\left(c_{0}\bm{u}-\bm{\mu}_{0}\right)^{2}_{1}+\left(\bm{H}_{0}\right)_{1,1}\quad\quad\quad\text{
for }i=0\\\
\bm{V}_{i}^{\top}\bm{A}(0)\bm{V}_{i}+\kappa_{0}\left(c_{0}\bm{u}-\bm{\mu}_{0}\right)^{2}_{i+1}+\left(\bm{H}_{0}\right)_{i+1,i+1}\quad\quad\quad\text{
for }i\neq 0\text{ and }\forall c_{0}\end{cases}.$ (3.10)

This follows from the expression of the matrix $\bm{B}(\bm{X},c_{0},\bm{u})$
in equation in (2) and the definition of $\bm{H}_{N}$. We observe the
following

$\displaystyle\left(\bm{H}_{N}\right)_{i+1,i+1}$
$\displaystyle\leq\lambda_{1}\left\\{\bm{A}(\bm{0})\right\\}+\kappa_{0}\|c_{0}\bm{u}-\bm{\mu}_{0}\|^{2}+\left(\bm{H}_{0}\right)_{i+1,i+1}$
$\displaystyle=m_{i}+\kappa_{0}\|c_{0}\bm{u}-\bm{\mu}_{0}\|^{2}\quad\text{ for
}i\neq 0\text{ and }\forall c_{0}\in\mathbb{R}.$ (3.11)

For a given value of $c_{0}$, equation (3.9), provides a lower bound
$h(\bm{u})$, which we will maximize to approximate $\bm{u}$. This is
equivalent to minimizing $-h(\bm{u})$. We apply Newton - Raphson algorithm
(Lange, , 2013) to obtain an update for $\bm{u}$. The calculation of
derivatives are shown below.

$\displaystyle\log p(\bm{u},c_{0},\widehat{\bm{\lambda}}\mid\bm{X})\geq
h(\bm{u}),$
$\displaystyle-h(\bm{u})=t\sum_{i=1}^{p-1}\log\left[\frac{m_{i}+\kappa_{0}\|c_{0}\bm{u}-\bm{\mu}_{0}\|^{2}}{2t}\right]+\frac{1}{2}\left[\bm{u}^{\top}\bm{A}(c_{0}\bm{u})\bm{u}+\kappa_{0}\|c_{0}\bm{u}-\bm{\mu}_{0}\|^{2}_{1}+m_{1}\right],$
(3.12) $\displaystyle-\nabla
h(\bm{u})=t\sum_{i=1}^{p-1}\frac{4t\kappa_{0}c_{0}(c_{0}\bm{u}-\bm{\mu}_{0})}{m_{i}+\kappa_{0}\|c_{0}\bm{u}-\bm{\mu}_{0}\|^{2}}+\left[\bm{A}(0)\bm{u}-c_{0}n\bm{\bar{x}}+\kappa_{0}c_{0}(c_{0}\bm{u}-\bm{\mu}_{0})\right],$
(3.13)
$\displaystyle-\nabla^{2}h(\bm{u})=t\sum_{i=1}^{p-1}\frac{4t\kappa_{0}c_{0}\left[(m_{i}+\kappa_{0}\|c_{0}\bm{u}-\bm{\mu}_{0}\|^{2})\bm{I}-2(c_{0}\bm{u}-\bm{mu}_{0})(c_{0}\bm{u}-\bm{\mu}_{0})^{\top}\right]}{\left\\{m_{i}+\kappa_{0}\|c_{0}\bm{u}-\bm{\mu}_{0}\|^{2}\right\\}^{2}},$
$\displaystyle\quad\quad\quad\quad\quad\quad+\left[\bm{A}(0)+\kappa_{0}c^{2}_{0}(c_{0}\bm{u}-\bm{\mu}_{0})(c_{0}\bm{u}-\bm{\mu}_{0})^{\top}\right]$
(3.14)

Thus, the updating equation for the mean direction is

$\displaystyle\bm{u}^{(k+1)}=\bm{u}^{(k)}-\left[\nabla^{2}h(\bm{u})\right]^{-1}\nabla
h(\bm{u})$ (3.15)

Due to the concavity of the lower bound, Newton’s method is an appealing
choice in our case. However there are two potential problems with Newton’s
method (Lange, , 2013, §10.3). First, it may be computationally expensive to
invert the second derivative matrix in each step. Second, the Newton’s method
is not really a ascent algorithm in the sense that
$h\left(\bm{u}^{(k+1)}\right)>h\left(\bm{u}^{(k)}\right)$ for a concave
function. The second problem can be remedied by modifying the increment such
that it is a partial step in the ascent direction. Let us denote:

$\displaystyle\bm{v}=\left[\nabla^{2}h(\bm{u})\right]^{-1}\nabla h(\bm{u})$
(3.16)

The idea is to take a sufficiently short increment in the direction of
$\bm{v}$. If $\left[\bm{u}^{(k)}-\alpha\bm{v}\right]$ shows increment in
$h(\bm{u})$ value, we update our mean vector in the iteration, otherwise we
look at $\left[\bm{u}^{(k)}-\alpha^{j}\bm{v}\right]$ for $j=1,2,\dots$ until
we observe an increment. Due to the good performance of the MLE, we use the
MLE of $c_{0}$ and $\bm{u}$ as our initial values which make the convergence
fast. Based on (3.8) and (3.15), an iterative algorithm for computing the
approximate MAP is summarized in the following algorithm:

Algorithm 2 MAP Approximation

1:Initialize: Start with $c^{(0)}_{0}=$ MLE of $c_{0}$, $\bm{u}^{(0)}=$ MLE of
$\bm{u}$ and
$\lambda^{(0)}_{i}=\frac{\bm{H}_{N}\left(c^{(0)},\bm{u}^{(0)}\right)_{i+1,i+1}}{2t}$

2:For  $k\rightarrow(k+1)$

3: Update $c^{(k+1)}$ and $\lambda^{(k+1)}_{i}$ from equation (3.8)

4: For $j\rightarrow(j+1)$

5: $\bm{u}^{(j)}=\bm{u}^{(k)}$ and
$\bm{v}^{(j)}=\left[\nabla^{2}h(\bm{u}^{(j)})\right]^{-1}\nabla
h(\bm{u}^{(j)})$

6: For  $l\rightarrow(l+1)$

7: $\bm{u}^{(j+1)}=\bm{u}^{(j)}-\alpha^{l}\bm{v}^{(j)}$

8: If  $h(\bm{u}^{(j+1)})>h(\bm{u}^{(j)})$

9: Accept the update of $\bm{u}$

10: Break

11: Else  Change $\alpha^{l}\rightarrow\alpha^{l+1}$

12: If  $\|\bm{u}^{(j+1)}-\bm{u}^{(j)}\|<\epsilon$

13: $\bm{u}^{(k+1)}=\bm{u}^{(j+1)}$

14: Break

15: If  $h(\bm{u}^{k+1})<h(\bm{u}^{(k)})$

16: Break

17:Get the final value of $c_{0}$, $\bm{\lambda}$ and $\bm{u}$.

18:End

This algorithm produces the MAP approximation of the eigenvalues and the mean
vector. The continuity of the estimate of the covariance matrix, shown in
Lemma 1, ensures that the covariance estimate will converge with the
convergence of the mean direction vector $\bm{u}$.

######  Lemma 1.

If $\bm{V}(\bm{u})$, described after (1.3) is continuous, then the approximate
MAP of the covariance matrix $\bm{\Sigma}$ is also a continuous function of
$\bm{u}$.

######  Proof.

If $\bm{V}(\bm{u})$ is a continuous function of $\bm{u}$, then each column of
$\bm{V}$, $\bm{V}_{j}$ is also continuous. This implies diagonal entries of
$\bm{H}_{N}$ and in turn the MAP estimate of the eigenvalues are also a
continuous function of $\bm{u}$. The claim follows since the MAP of
$\bm{\Sigma}$ given by

$\displaystyle\widehat{\bm{\Sigma}}$
$\displaystyle=\bm{u}\bm{u}^{\top}+\sum_{i=1}^{p-1}\widehat{\lambda}_{i}\bm{V}_{i}\bm{V}^{\top}_{i}=\bm{u}\bm{u}^{\top}+\sum_{i=1}^{p-1}\frac{\left(\bm{H}_{N}\right)_{i+1,i+1}}{n}\bm{V}_{i}\bm{V}^{\top}_{i}$
(3.17)

is sum of continuous functions of $\bm{u}$.

##  4 Simulations

In this section, we perform several simulations to assess the performance of:

  1. 1.

MLE approximation through lower bound maximization (see section 2) of the
profile likelihood function,

  2. 2.

MAP estimator with normal-inverse gamma prior (see section 3.3) approximated
through Gibbs sampling in section 3.4, and

  3. 3.

MAP approximation (see section 3.5) through posterior lower bound

by computing their risks. Then the three estimators are compared with the MAP
estimator with normal-inverse Wishart prior (see section 3.1). We have used
RStudio 1.3.1093 and R 4.0.3 on a 64 bit 4 Core Windows 10 laptop for all our
simulations.

###  4.1 The Simulation Set up:

We have taken the sample sizes and dimensions to be $n=50,100,300$ and
$p=3,5,10$, respectively. In all cases the data generation mechanism and the
risk function are kept the same. We have used Frobenius loss as our default
loss function and risks are approximated by averaging the losses for 100
independent replications in each of the nine combinations of $(n,p)$.

The parameters of the Gaussian distributions used for data generation are
selected in the following way:

  * •

the entries of the mean vector $\bm{\mu}$ are taken to be independent standard
Gaussian variables,

  * •

the covariance matrix is generated from $\bm{\Psi}=\bm{LL}^{\top}$ where

$\displaystyle L$ $\displaystyle=\begin{cases}L_{ij}\sim N(0,1)\quad\text{ for
}i\neq j,\\\ L_{ij}\sim N(5,1)\quad\text{ for }i=j.\end{cases}$

The larger diagonal entries of $\bm{L}$ ensure positive-definiteness of the
modified covariance matrix $\bm{\Sigma}$.

  * •

Since such $(\bm{\mu},\bm{\Psi})$ does not necessarily satisfy conditions
(1.1), the covariance matrix $\bm{\Sigma}$ is constructed from applying Lemma
4 of Kundu and Pourahmadi,  (2020) on $(\bm{\mu},\bm{\Psi})$.

The performance of the estimators is assessed using the Frobenius (scaled
$L_{2}$) risk as in Ledoit and Wolf, 2004b (, §3.1):

$R(\bm{\mu},\widehat{\bm{\mu}}^{*})=\mathbb{E}\left[\frac{1}{p}\left\|\widehat{\bm{\mu}}^{*}-\bm{\mu}\right\|^{2}_{\mathcal{F}}\right]\qquad,\qquad
R(\bm{\Sigma},\widehat{\bm{\Sigma}}^{*})=\mathbb{E}\left[\frac{1}{p}\left\|\widehat{\bm{\Sigma}}^{*}-\bm{\Sigma}\right\|^{2}_{\mathcal{F}}\right],$

where $\widehat{\bm{\mu}}^{*}$ and $\widehat{\bm{\Sigma}}^{*}$ are the final
estimators in each of the cases.

###  4.2 Simulation Results:

Here we provide the details of the simulation studies for the MLE
approximation, MAP estimator computed through Gibbs sampling and using a lower
bound of the posterior density. We have used normal-inverse Wishart as a
yardstick to compare the risk of our estimators (MLE approximation and the two
MAP’s).

####  4.2.1 Gibbs Sampling

The risk of the MAP estimator is computed through Gibbs sampling as described
in Section 3.4. Let us denote the initial values to be
$\left(\bm{\mu}^{(0)},\bm{\lambda}^{(0)}\right)$. We have chosen
$\bm{\mu}^{(0)}=\bm{\mu}_{0}=\bar{\bm{x}},\bm{H}_{0}=\bm{I},\kappa_{0}=1.5$
and $a=(p+1)$. The first iteration takes in $\bm{\mu}^{(0)}$ and calculates
the value of $\bm{H}_{N}$ as described after equation (3.5). Using this, we
can generate the value of $\bm{\lambda}^{(1)}$ from the distribution
$p(\bm{D}\mid\bm{\mu},\bm{X})$ (see step 1 of Gibbs sampling in section 3.4)
which by construction makes sure that each eigenvalue is positive. We can
easily see that the initial value of $\bm{\lambda}$ i.e. $\bm{\lambda}^{(0)}$
does not affect the Gibbs sampling.

We have generated 100 samples from the posterior distribution using the MH
within Gibbs algorithm described in section 3.4. We repeat this experiment 100
times and approximate the risk. The performance of the Gibbs sampling is
displayed in the following table LABEL:tab:Gibbstable. We can observe that the
performance is improving as the number of observations $(n)$ increase. For low
dimensional cases the risk of the MAP estimator of the covariance matrix
calculated from Gibbs sampling is equivalent and better as compared to its
normal-inverse Wishart counterpart.

A well-known problem with Gibbs sampling is that as the dimension increases
the time complexity of the MAP approximation increases significantly. The
acceptance rate of the MH step within Gibbs sampling is reasonable and
decreases with in increment of dimension.

Table 1: Risk ratio of MAP approximation (from MH within Gibbs sampling) relative to the MAP of normal-inverse Wishart |  |  | Normal-Inverse
---
Gamma Risk
Acc Rate (in MH) |  Time (Sec.)
n | p | Mean | Sigma
50 | 3 | 1.1173 | 0.8928 | 0.4135 | 980.52
50 | 5 | 1.0824 | 1.0963 | 0.2822 | 2029.15
50 | 10 | 1.2423 | 1.2798 | 0.1218 | 13574.59
100 | 3 | 1.0542 | 0.9872 | 0.4352 | 910.92
100 | 5 | 1.1677 | 1.3662 | 0.2908 | 1901.41
100 | 10 | 1.1963 | 1.6229 | 0.1245 | 12183.56
300 | 3 | 1.1661 | 1.2546 | 0.4287 | 1256.81
300 | 5 | 1.5729 | 2.0163 | 0.2895 | 3215.97
300 | 10 | 1.5533 | 2.4053 | 0.1294 | 9385

####  4.2.2 MLE and MAP Approximation using a Lower Bound

The structured covariance matrix defined in equation (1.3) does not allow us
to calculate the maximum likelihood estimate and the MAP estimator directly
due to the intractability of the likelihood function. Hence, we approximate
them through a lower bound described in Sections 2 and 3.5, respectively. The
approximation of MLE performs well.

Table 2: Risk ratio of MLE and MAP approximation (through a lower bound) relative to the MAP of normal-inverse Wishart |  |  | MLE Approx.
---
Risk
| MAP Approx.
---
Risk
Time (Sec)
n | p | Mean | Sigma | Mean | Sigma
50 | 3 | 0.4253 | 1.1331 | 0.4253 | 1.1331 | 3.25
50 | 5 | 0.6625 | 1.1528 | 0.6625 | 1.1528 | 2.59
50 | 10 | 1.5009 | 1.195 | 1.5009 | 1.195 | 5.22
100 | 3 | 0.3481 | 1.4383 | 0.3481 | 1.4383 | 2.52
100 | 5 | 0.5342 | 1.5065 | 0.5342 | 1.5065 | 2.69
100 | 10 | 1.401 | 1.63 | 1.401 | 1.63 | 2.84
300 | 3 | 0.3159 | 2.3753 | 0.3159 | 2.3753 | 2.7
300 | 5 | 0.5797 | 2.6429 | 0.5797 | 2.6429 | 2.95
300 | 10 | 1.303 | 2.6493 | 1.303 | 2.6493 | 3.41

A modified version of the Newton-Raphson method (described in the paragraph
following equation (3.15)) converges relatively fast i.e. two to five
iterations for MAP when initialized with the approximation of MLE of $c_{0}$
and $\bm{u}$.

The performance of these approximations (both MLE and MAP), assessed by the
risk, are equivalent in dimension three and five to the Gibbs sampling
technique discussed in the last section with the added benefit of being
significantly faster (see the last column of table LABEL:tab:mlemapNR). For
$p=10$, the performance is worse than Gibbs sampling (see Table
LABEL:tab:mlemapNR). The usual observations like the decline of risk as $n$
increases and $p$ decreases are generally true. One interesting observation is
that the performance of MLE approximation is not significantly different from
the MAP approximation.

####  4.2.3 An Example: Estimates of the Historic Position of Earth’s Magnetic
Pole

The dataset collected by Schmidt,  (1976) contains the site mean direction
estimates of the Earth’s historic magnetic pole collected from 33 different
sites in Tasmania. The longitude and latitudes from the data set is
transformed to $X_{1},X_{2},\dots,X_{33}$ on a three dimensional unit sphere
(Preston and Paine, , 2017). The angualr gaussian distribution family is the
marginal directional component of a multivariate normal distribution with ESAG
distribution as a subfamily. Paine et al.,  (2018) provided strong evidence in
favor of ESAG distribution which satisfy the constraint over isotropic angular
Gaussian distribution while analyzing this dataset. This inspire us to make
normality assumption under the constraint similar to ESAG distribution
disregarding the spherical nature of the transformed dataset. The constrained
maximum likelihood estimate under the structured covariance model is:

$\widehat{\bm{u}}^{\top}=\begin{bmatrix}-0.44&0.32&0.75\end{bmatrix},\quad\widehat{\bm{\Sigma}}=\begin{bmatrix}1.63&-0.91&-2.13\\\
-0.91&1.09&1.53\\\ -2.13&1.53&4.02\end{bmatrix}$

which are comparable to the maximum likelihood estimate calculated using
elliptically symmetric angular Gaussian distribution with the parametrization
in three dimension (Paine et al., , 2018) i.e.
$\widehat{\bm{u}}^{\top}=(-0.56,0.24,0.79)$.

##  5 Discussion

This article demonstrates how simple ideas can be used to solve a challenging
joint estimation problem under constraint (1.1). The structured covariance
model takes into account of the constraint, reduces dimensions significantly
and do not suffer from the usual issue of positive definiteness. The
intractability of the model is resolved through a lower bound of the profile
likelihood to obtain MLE, MAP and sampling from posterior density is performed
through Metropolis within Gibbs. We have proposed a simple prior for our
structured model i.e. normal on mean vector and inverse gamma on eigenvalues
inspired by Berger et al.,  (2020). Although this model works nicely here for
comparatively smaller dimensions, it is possible to replicate the idea in
higher dimensions with some modifications. Other priors can also be explored
under such a model in future.

## References

  * Barnard et al., (2000) Barnard, J., McCulloch, R., and Meng, X.-L. (2000).  Modeling covariance matrices in terms of standard deviations and correlations, with application to shrinkage.  Statistica Sinica, pages 1281–1311.
  * Berger et al., (2020) Berger, J. O., Sun, D., Song, C., et al. (2020).  Bayesian analysis of the covariance matrix of a multivariate normal distribution with a new class of priors.  Annals of Statistics, 48(4):2381–2403.
  * Bibby et al., (1979) Bibby, J., Kent, J., and Mardia, K. (1979).  Multivariate analysis.  Academic Press, London.
  * Brooks et al., (2011) Brooks, S., Gelman, A., Jones, G., and Meng, X.-L. (2011).  Handbook of markov chain monte carlo.  CRC press.
  * Chiu et al., (1996) Chiu, T. Y., Leonard, T., and Tsui, K.-W. (1996).  The matrix-logarithmic covariance model.  Journal of the American Statistical Association, 91(433):198–210.
  * Eguchi et al., (2010) Eguchi, N., Saito, R., Saeki, T., Nakatsuka, Y., Belikov, D., and Maksyutov, S. (2010).  A priori covariance estimation for co2 and ch4 retrievals.  Journal of Geophysical Research: Atmospheres, 115(D10).
  * Fan et al., (2008) Fan, J., Fan, Y., and Lv, J. (2008).  High dimensional covariance matrix estimation using a factor model.  Journal of Econometrics, 147(1):186–197.
  * Gelfand and Smith, (1990) Gelfand, A. E. and Smith, A. F. (1990).  Sampling-based approaches to calculating marginal densities.  Journal of the American statistical association, 85(410):398–409.
  * Gelman and Hill, (2006) Gelman, A. and Hill, J. (2006).  Data analysis using regression and multilevel/hierarchical models.  Cambridge university press.
  * Gelman et al., (2013) Gelman, A., Stern, H. S., Carlin, J. B., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013).  Bayesian data analysis.  Chapman and Hall/CRC.
  * Guillot et al., (2015) Guillot, D., Rajaratnam, B., Emile-Geay, J., et al. (2015).  Statistical paleoclimate reconstructions via markov random fields.  The Annals of Applied Statistics, 9(1):324–352.
  * Hamimeche and Lewis, (2009) Hamimeche, S. and Lewis, A. (2009).  Properties and use of cmb power spectrum likelihoods.  Physical Review D, 79(8):083012.
  * Hamsici and Martinez, (2007) Hamsici, O. C. and Martinez, A. M. (2007).  Spherical-homoscedastic distributions: The equivalency of spherical and normal distributions in classification.  Journal of Machine Learning Research, 8(7).
  * Hoff, (2009) Hoff, P. D. (2009).  A hierarchical eigenmodel for pooled covariance estimation.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(5):971–992.
  * Hoff and Niu, (2012) Hoff, P. D. and Niu, X. (2012).  A covariance regression model.  Statistica Sinica, pages 729–753.
  * Kundu and Pourahmadi, (2020) Kundu, A. and Pourahmadi, M. (2020).  Mle of jointly constrained mean-covariance of multivariate normal distributions.  arXiv preprint arXiv:2012.11826.
  * Lange, (2013) Lange, K. (2013).  Elementary optimization.  In Optimization, pages 1–17. Springer.
  * (18) Ledoit, O. and Wolf, M. (2004a).  Honey, i shrunk the sample covariance matrix.  The Journal of Portfolio Management, 30(4):110–119.
  * (19) Ledoit, O. and Wolf, M. (2004b).  A well-conditioned estimator for large-dimensional covariance matrices.  Journal of multivariate analysis, 88(2):365–411.
  * Mardia and Jupp, (2000) Mardia, K. and Jupp, P. (2000).  Directional statistics.
  * Murphy, (2012) Murphy, K. P. (2012).  Machine learning: a probabilistic perspective.  MIT press.
  * Paine et al., (2018) Paine, P., Preston, S. P., Tsagris, M., and Wood, A. T. (2018).  An elliptically symmetric angular gaussian distribution.  Statistics and Computing, 28(3):689–697.
  * Paine et al., (2020) Paine, P. J., Preston, S., Tsagris, M., and Wood, A. T. (2020).  Spherical regression models with general covariates and anisotropic errors.  Statistics and Computing, 30(1):153–165.
  * Pourahmadi, (1999) Pourahmadi, M. (1999).  Joint mean-covariance models with applications to longitudinal data: Unconstrained parameterisation.  Biometrika, 86(3):677–690.
  * Preston and Paine, (2017) Preston, S. and Paine, P. (2017).  Analysis of spherical data with ESAG.
  * Rao, (1973) Rao, C. R. (1973).  Linear statistical inference and its applications.  Wiley New York, second edition.
  * Schäfer and Strimmer, (2005) Schäfer, J. and Strimmer, K. (2005).  A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics.  Statistical applications in genetics and molecular biology, 4(1).
  * Schmidt, (1976) Schmidt, P. (1976).  The non-uniqueness of the australian mesozoic palaeomagnetic pole position.  Geophysical Journal International, 47(2):285–300.
  * Trefethen and Bau III, (1997) Trefethen, L. N. and Bau III, D. (1997).  Numerical linear algebra, volume 50.  Siam.
  * Yang and Berger, (1994) Yang, R. and Berger, J. O. (1994).  Estimation of a covariance matrix using the reference prior.  The Annals of Statistics, pages 1195–1211.

##  Appendix A Proofs of Results:

  1. 1.

In case of normal-inverse Wishart prior the posterior density with parameter
$(\bm{\mu}_{n},\kappa_{n},\nu_{n},\bm{\Lambda}_{n})$ is the following:

$\displaystyle
p(\bm{\mu},\bm{\Sigma}\mid\bm{\mu}_{n},\kappa_{n},\nu_{n},\bm{\Lambda}_{n})$
$\displaystyle\propto\kappa_{n}^{-\frac{p}{2}}\mid\bm{\Sigma}\mid^{-\left(\frac{\nu_{n}+p}{2}+1\right)}\exp\left[-\frac{1}{2}\left\\{\text{Tr}\left(\bm{\Sigma}^{-1}\bm{\Lambda}_{n}\right)+\kappa_{n}\left(\bm{\mu}-\bm{\mu}_{n}\right)^{\top}\bm{\Sigma}^{-1}\left(\bm{\mu}-\bm{\mu}_{n}\right)\right\\}\right]$

After taking logarithm and differentiating we get :

$\displaystyle\frac{\partial\log
p(\bm{\mu},\bm{\Sigma}\mid\bm{\mu}_{n},\kappa_{n},\nu_{n},\bm{\Lambda}_{n})}{\partial\bm{\Sigma}^{-1}}$
$\displaystyle=\left(\frac{\nu_{n}+p}{2}+1\right)\bm{\Sigma}-\frac{1}{2}\bm{\Lambda}_{n}-\frac{\kappa_{n}}{2}\left(\bm{\mu}-\bm{\mu}_{n}\right)\left(\bm{\mu}-\bm{\mu}_{n}\right)^{\top}=0$
$\displaystyle\widehat{\bm{\Sigma}}_{map}$
$\displaystyle=\frac{\bm{\Lambda}_{n}+\kappa_{n}\left(\bm{\mu}-\bm{\mu}_{n}\right)\left(\bm{\mu}-\bm{\mu}_{n}\right)^{\top}}{\nu_{n}+p+2}$

By applying Theorem 4.2.1 of Bibby et al.,  (1979) we can say that this is
indeed the MAP estimator.

  2. 2.

Let $\bm{x}_{1},\bm{x}_{2},\dots,\bm{x}_{n}\sim N_{p}(\bm{\mu},\bm{\Sigma})$
and the prior on $(\bm{\mu},\bm{\Sigma})$ is normal-shrinkage inverse Wishart
prior i.e. $\bm{\mu}\mid\bm{\Sigma}\sim
N_{p}(\bm{\mu}_{0},\bm{\Sigma}/\kappa_{0})$ and
$\bm{\Sigma}\sim\pi_{SIW}(\nu_{0},b,\bm{\Lambda}_{0}^{-1})$. Here w are
interested with $b=1$ as stated earlier. The posterior distribution is
calculated as follows:

$\displaystyle\pi(\bm{\mu},\bm{\Sigma}\mid\bm{x}_{1},\dots,\bm{x}_{n},\bm{\mu}_{0},\kappa_{0},\nu_{0},\bm{\Lambda}_{0})$
$\displaystyle\propto\prod_{i=1}^{n}N_{p}(\bm{x}_{i};\bm{\mu},\bm{\Sigma})N_{p}(\bm{\mu};\bm{\mu}_{0},\bm{\Sigma}/\kappa_{0})\pi_{SIW}(\bm{\Sigma};\nu_{0},1,\bm{\Lambda}_{0}^{-1})$
$\displaystyle\propto\prod_{i=1}^{n}N_{p}(\bm{x}_{i};\bm{\mu},\bm{\Sigma})N_{p}(\bm{\mu};\bm{\mu}_{0},\bm{\Sigma}/\kappa_{0})\frac{\pi_{IW}(\bm{\Sigma};\nu_{0},\bm{\Lambda}_{0}^{-1})}{\prod_{i<j}(\lambda_{i}-\lambda_{j})}$

This is because shrinkage-inverse-Wishart can be obtained from inverse-Wishart
density by dividing it with $\prod_{i<j}(\lambda_{i}-\lambda_{j})$. We know
that the normal-inverse Wishart prior is conjugate with posterior parameters
described in (3.1), $\kappa_{n}=\kappa_{0}+n$ and $\nu_{n}=\nu_{0}+n$
respectively. The numerator is the posterior distribution corresponding to the
normal-inverse Wishart prior and gives us the following:

$\displaystyle\pi(\bm{\mu},\bm{\Sigma}\mid\bm{x}_{1},\dots,\bm{x}_{n},\bm{\mu}_{0},\kappa_{0},\nu_{0},\bm{\Lambda}_{0})$
$\displaystyle\propto\frac{\pi_{NIW}(\bm{\mu},\bm{\Sigma}\mid\bm{\mu}_{n},\kappa_{n},\nu_{n},\bm{\Lambda}_{n}^{-1})}{\prod_{i<j}(\lambda_{i}-\lambda_{j})}$
$\displaystyle\propto
N_{p}(\bm{\mu}\mid\bm{\Sigma},\bm{\mu}_{n},\kappa_{n})\frac{\pi_{IW}(\bm{\Sigma}\mid\nu_{n},\bm{\Lambda}_{n}^{-1})}{\prod_{i<j}(\lambda_{i}-\lambda_{j})}$
$\displaystyle\propto
N_{p}(\bm{\mu}\mid\bm{\Sigma},\bm{\mu}_{n},\kappa_{n})\pi_{SIW}(\bm{\Sigma}\mid\nu_{n},1,\bm{\Lambda}_{n}^{-1})$

From the above calculation we can see that the posterior is normal-shrinkage
inverse Wishart with the same parameters as in normal-inverse Wishart proving
our hypothesis.

  3. 3.

Calculation of the Posterior Distribution:

$\displaystyle p(\bm{\mu},\bm{\lambda}\mid\bm{X})$
$\displaystyle\propto\left(\prod_{i=1}^{p-1}\lambda_{i}\right)^{-\frac{n}{2}}\exp\left[-\frac{1}{2}Tr\left\\{\bm{\Lambda}^{-1}B(\bm{X},\bm{\mu})\right\\}\right]\left(\prod_{i=1}^{p-1}\lambda_{i}\right)^{-\frac{1}{2}}$
$\displaystyle\quad\quad\quad\exp\left[-\frac{1}{2}Tr\left\\{\bm{\Lambda}^{-1}(\bm{\mu}-\bm{\mu}_{0})(\bm{\mu}-\bm{\mu}_{0})^{\top}\right\\}\right]\left(\prod_{i=1}^{p-1}\lambda_{i}\right)^{-a}\exp\left[-\frac{1}{2}Tr\left\\{\bm{\Lambda}^{-1}\bm{H}_{0}\right\\}\right]$
$\displaystyle\propto\left(\prod_{i=1}^{p-1}\lambda_{i}\right)^{-\frac{n+2(a+1)}{2}}\exp\left[-\frac{1}{2}Tr\left\\{\bm{\Lambda}^{-1}\left(B(\bm{X},\bm{\mu})+(\bm{\mu}-\bm{\mu}_{0})(\bm{\mu}-\bm{\mu}_{0})^{\top}+\bm{H}_{0}\right)\right\\}\right]$
$\displaystyle\propto\left(\prod_{i=1}^{p-1}\lambda_{i}\right)^{-\frac{n+2(a+1)}{2}}\exp\left[-\frac{1}{2}Tr\left\\{\bm{\Lambda}^{-1}\bm{H}_{N}\right\\}\right]$

where
$B(\bm{X},\bm{\mu})=\bm{P}(\bm{\mu})^{\top}\bm{A}(\bm{\mu})\bm{P}(\bm{\mu})$.


# Spinoza, Leibniz, Kant, and Weyl

Michael H. Freedman  Michael H. Freedman
Microsoft Research, Station Q, and Department of Mathematics
University of California, Santa Barbara
Santa Barbara, CA 93106

I walk into a House of Worship and the guy up front is in the middle of an
oration. I don’t expect to stay but catch a few words, “… so my friends, reach
down and uplift thy less fortunate neighbor—to the $\frac{1}{2}$-power, for so
it is written, and so may it be derived.” Startled, I take a seat; this sounds
like theology on a wavelength I can receive. This was my experience meeting
(virtually) Glen Weyl and reading his paper with Buterin and Hitzig [bhw]. The
paper describes a funding mechanism called _Quadratic Finance_ (QF) and
deploys a bit of calculus to show that within a very clean and simple linear
model QF maximizes social utility. They differentiate the social utility
function. The mathematical content of this note is that by taking one further
derivative, one may also deduce that QF is the _unique_ solution. But what
drew me into this house of worship, is a line of reasoning in [BHW], made
explicit in a phone call with Glen, that Kant’s _categorical imperative_ (CI),
properly interpreted, becomes a differential equation uniquely solved by QF.
It was this second derivation from Kant’s “axiom” that caused me to take a
seat, and for a very personal reason.

My father, Benedict Freedman [bf1, bf2], was a novelist, screenwriter, and
educator, trained in philosophy and mathematical logic. He lived in the realm
spanning Spinoza to Gödel; he sought, all his life, to expand upon Spinoza’s
Ethics, an astonishingly ambitious project to logically deduce how to act in
our world from a minimal axiomatic foundation of core beliefs. As
psychologists know, beliefs are generally context dependent and fuzzy. What
kind of logical machinery can they survive? Is it realistic to deduce ethics?
One thought (I hope to develop later), which I believe my father would
endorse, is to use a logic with a web-like rather than a linear structure, a
pluralist’s logic, a logic which is able to weigh preponderance of evidence
from various pathways much as Feynmann diagrams in perturbative quantum field
theory must be summed over to reach a reliable conclusion.

As I learned from Jerry Muller [muller], Western political and social
structure, both in theory and practice, can be viewed as a secularization of
Judeo-Christian religious tradition. Our quest for understanding has this
broader context. This view may also be applied to the sciences, and
comfortably so, in the case of economics. Marxism is a prime exemplar, with
its emphasis on a goal-directed, teleological, historiography replete with
catastrophism and a utopian eschatology. My father once said, “All the great
Jews of history were atheists: Spinoza, Marx, Einstein.” He added that these
three each demonstrated deep religious sensibilities. (My father, a modest
man, merely considered himself a Jewish agnostic. When I asked him as a
teenager how a Jew could be an atheist, he answered that atheism did not come
from nowhere—it was invented by Jews, it is a sect within Judaism.) Spinoza’s
religious perspective encouraged his belief that there are anchor points
strong enough to sustain deduction. Alas, I was never able to follow Spinoza’s
deductions, a shortcoming, I suspect, of a too narrowly mathematical mind,
always demanding clarity of definition. For fifty years I gave all this little
thought, until seeing last month, in the deduction of QF from CI, a first
glimmer of success for Spinoza’s (and my father’s) project. That is, a success
that I could understand. Imagining the future, Leibniz wrote, “When there are
disputes among persons, we can simply say, ‘Let us calculate’, and without
further ado, see who is right.” In Leibniz’s future, mathematics may usefully
manipulate ethical propositions. I had not encountered a convincing example of
this prior to CI $\rightarrow$ QF. Is this future closer than we thought? New
tools open new domains: Can an economic perspective, certainly foreign to
Spinoza and Leibniz, lead to an ethical calculus?

Kant’s CI is more subtle than the much remembered, “A deed which must be done
though the heavens fall.” Kant wrote in _Groundwork of the Metaphysics of
Morals_ , “Act only according to that maxim whereby you can at the same time,
will that it should become universal law.” This rejection of exceptionalism
clearly has roots in the Golden Rule but is now couched unconditionally in a
more abstract form. All those interested in distributed architectures and
shared governance can feel the Kantian vibe. To accept CI requires a degree of
a liberal, rationalist mind set, but superficially, it seems consistent with a
wide spectrum of social and political views. But once CI is accepted, I will
argue, following [bhw], that it implies QF, an apparently more generous
doctrine which in the orator’s words, uplifts the weaker players according the
$\frac{1}{2}$-power. The surprise is that a calculation is required—Leibniz
take note—to learn the moral scope of one’s own maxim. Curiously, the
conclusion is fine-tuned, one is not lead to extremes of egalitarianism, all
social differences are not leveled, but rather a certain formula for adjusting
then, involving square roots and squares emerges. There is a context to the
argument, standard assumptions of linearity and smoothness and concavity of
individual utility functions are made, as well as the insensitivity of
individual preferences to some level of taxation (to cover the shortfall
between contributions and dispersed funds). Nevertheless, it is remarkable how
little beyond CI is required to deduce QF. The model is simple and thus far
from complete, yet seeing this deduction within the model was a moving
experience, if there are other examples of “ethical calculus,” I would like to
learn about them.

Let’s review QF. QF is an answer to a question: How should public goods be
financed. A “public good,” such as a public park, will typically be a good
that is inefficient to price and restrict, and one from which individuals
receive value greater than their personal contributions. Below is a brief
summary of the model presented in [bhw].

_Society_ consists of $n$ well-defined citizens, $i=1,\dots,n$. $p\in P$ are
the public goods in question. (In our treatment we simply consider a single
good, and eventually drop the $p$-index. It has no real effect on the
calculus.)

Let $V_{i}^{p}(F^{p})$ be the currency-equivalent utility citizen $i$ receives
if the public funding of $p$ is $F^{p}$. These functions are assumed to be
smooth, increasing and concave.111Even this assumption can often be relaxed to
monotonicity. Remarkably, nothing else about the functional form of the
citizen’s values needs to be assumed. We further assume all these values
(utilities) are independent, i.e. just add, and ignore any issues of partial
information and timing of decision making within the model.

$\vec{c}=\\{c_{i}^{p}\\}$ is the vector of individual contributions of the ith
citizen toward the $p$th public good, which are all assumed non-negative (a
non-trivial assumption); and $\vec{F}$ is the vector of funding with
components $F^{p}$ for each public good. Let $C$ and $F$ denote the vector
spaces holding the vectors $\vec{c}$ and $\vec{F}$.

The problem considered: Find a (and eventually all) funding mechanism
$\Phi:C\rightarrow F$ which maximizes the total social welfare
$W=\sum_{i,p}V_{i}^{p}(F^{p})-\sum_{p}F^{p}$. We do not consider within the
model the equity of the mechanism, which must be treated externally to the
model. And very importantly we do not analyze here the back reaction on
citizen contributions c which taxation inevitably entails. Such an analysis is
initiated in section 4.5 of [bhw], and is shown to perturb the QF extremum,
but we have no further insight to add here. Taxation $\\{t_{i}\\}$ (according
to some other mechanism) is, of course, required to balance the budget:

$\sum_{i}t_{i}=\sum_{p}\left(F^{p}-\sum_{i}c_{i}^{p}\right)$

and the individual’s tax-corrected utility will be

$U_{i}^{t}=\sum_{p}V_{i}^{p}(F^{p})-c_{i}^{p}-t_{i}$

but I ignore taxes entirely in the calculations (as well as eventually
dropping the index $p$ for visual clarity).

## Optimality condition

Fixing $p$ and differentiating $W^{p}=\sum_{i}V_{i}^{p}(F^{p})-F^{p}$ w.r.t.
$F^{p}$ and setting this equal to zero, we find the marginal value derived
from good $p$ should equal 1, provided $F^{p}$ is positive at 0,
$\sum_{i}(V_{i}^{p})^{\prime}=1$. Note this sum is written ${V^{p}}^{\prime}$
in [bhw].

## Individual utility

$U_{i}=\sum_{p}V_{i}^{p}(F^{p})-c_{i}=\sum_{p}V_{i}^{p}(g(\sum_{j}h(c_{j}^{p}))-c_{i}^{p})$,
where the model presumes that the functions $F^{p}$ are built as indicated
from the internal functions $h$ and $g$, explained in the following paragraph.
Below, the superscript $t$ is dropped since we are now neglecting taxes, and
$p$ added to indicate the utility of a fixed good. Upon differentiating and
setting $\frac{\partial U_{i}^{p}}{\partial c_{i}}=0$, we obtain:

(1) $\frac{\partial V_{i}^{p}}{\partial g}\cdot\frac{\partial g}{\partial
h(c_{i}^{p})}\cdot\frac{dh(c_{i}^{p})}{dc_{i}^{p}}=1$

In what follows differential calculations will not be interrupted to consider
special case where $c_{i}^{p}=0$, nor to reiterated the obvious boundary
condition that $\Phi(0)=0$, no contributions implies no funding.

The possible mechanisms $\Phi$ are posited to be _democratic_ in that the
funding function $F$ is assumed to be symmetric in its $n$ variables, and
beyond that I adopted the _simplifying linearity assumption_ of [bhw] that
$\Phi(\vec{c})=\\{F^{P}(\vec{c})\\}=\\{g(\sum_{i}h(c_{i}^{p}))\\}$ for some
analytic function $g$ and $h$ on the positive Reals. Our job is to find $g$
and $h$ maximizing social utility. The internal function $h$ is the _weight_
of contribution $c$, and $g$ converts total weight into funding, it is the
_funding lever_. Notice in the definition of $\Phi^{\mathrm{QF}}$ below that
$h$ and $g$ scale reciprocally so funding choices are independent of the units
of the currency. However, this is an outcome not an assumption.

######  Definition 1.

$\Phi^{\mathrm{QF}}(c^{p})=\left(\sum_{i}(c_{i}^{p})^{\frac{1}{2}}\right)^{2}$,
that is, for every good $p$ its level of funding is the square of the sum of
the square roots of the individual contributions.

At this point, going forward, I drop the superscript $p$ and the set notation
$\\{\\}$. I will shortly set up the natural extremal problem for $\Phi$ and
show that $\Phi^{\mathrm{QF}}$ is its unique critical point, but first let’s
show how $\Phi^{\mathrm{QF}}$ may be deduced from CI.

The insight I learned from Glen, reflected the differential equation (10) of
[bhw], reproduced below, is that CI implies that if citizen $j$ deems it
proper to perturb her weighted contribution $h(c_{j})$, say by increasing it
1%, she should be following, not her limited self-interest, but be justified
in expecting all her peers to also see the virtue of such a similar
proportional increase in their weighted contribution—“act … whereby … it
should become universal law.” So, mathematically we may write:

(2) $\frac{\partial g(\sum_{i}h(c_{i}))}{\partial
c_{j}}=\sum_{i}\frac{h(c_{i})}{h(c_{j})}$

That is, the funding should respond to the imputed community wide judgement
that additional matched resources are required for this good. Importantly the
individual’s judgement and its imputed extension must be expressed in terms of
_weighted_ contributions, for it is the weighting function h that codifies how
any given funding mechanism (obeying our assumptions) hears the preferences of
its citizens.

From (2) QF follows directly. Differentiating we obtain

(3)
$g^{\prime}\left(\sum_{i}h(c_{i})\right)h^{\prime}(c_{j})=\sum_{i}\frac{h^{\prime}(c_{i})}{h(c_{j})}$

Equation 3 is actually two independent equations as the factors involving only
$c_{j}$ must separately be proportional, as must the factors involving all the
$\\{c_{i}\\}$. Separating and solving one finds that for some positive
constant $k$:

(4) $\displaystyle g^{\prime}(x)=kx\ \text{ and }\
h^{\prime}(y)=\frac{1}{k}h^{-1}(y),\text{ so}$ (5) $\displaystyle
g(x)=\frac{k}{2}x^{2}+m\ \text{ and }\ h(y)=\frac{2}{k}y^{\frac{1}{2}}+n$

The boundary conditions tells us that $m=n=0$, and $k$ must be set so that for
a society with a single citizen, we obtain $F(c)=c$. Thus we find:

(6) $g(x)=x^{2}\ \text{ and }\ h(y)=y^{\frac{1}{2}}$

Kant’s IC has been interpreted as equation 2; which, as we have seen, has QF
as its unique solution.

Let’s now return now to the calculus problem of extremizing social utility
$U=\sum_{i}(U_{i})$ by varying our funding mechanism $\Phi$ though our choice
of the internal functions $h$ and $g$. To do this, first on line (7) we
rewrite (1), and then on line (8) sum (7) over $i$, apply the optimality
condition, and then differentiate w.r.t. any one of the citizen contributions,
say $c_{1}$ (by symmetry there is no real choice here). The result is:

(7) $\displaystyle\frac{\partial V_{i}}{\partial g}=\frac{1}{\frac{\partial
g(\Sigma)}{\partial
h(c_{i})}\cdot\frac{dh(c_{i})}{dc_{i}}}=\frac{1/\frac{dh(c_{i})}{dc_{i}}}{\frac{\partial
g(\Sigma)}{\partial h(c_{i})}},\ \text{where
}\Sigma\coloneqq\sum_{i=1}^{n}h(c_{i})=\text{ total weight}$ (8)
$\displaystyle\frac{\partial}{\partial
c_{1}}\left(\sum_{i=1}^{n}\frac{\partial V_{i}}{\partial
g}\right)=\frac{\partial}{\partial
c_{1}}\left(\frac{1/\frac{dh(c_{i})}{dc_{i}}}{\frac{\partial
g(\Sigma)}{\partial h(c_{i})}}\right)=0$

Now expand the outer partial using
$\left(\frac{u}{v}\right)^{\prime}=\frac{u^{\prime}v-uv^{\prime}}{v^{2}}$,
keeping only the numerator, and breaking out the first variable $c_{1}$ from
the rest of the sum, we obtain:

(9) $0=\left(\frac{1}{h^{\prime}(c_{1})}\right)^{\prime}\cdot
g^{\prime}(\Sigma)-\frac{g^{\prime\prime}(\Sigma)h^{\prime}(c_{1})}{h^{\prime}(\Sigma)}+0-\sum_{i=2}^{n}\frac{1}{h^{\prime}(c_{i})}g^{\prime\prime}(\Sigma)h^{\prime}(c_{1})$

Collecting terms:

(10) $\displaystyle\left(\frac{1}{h^{\prime}(c_{1})}\right)^{\prime}\cdot
g^{\prime}(\Sigma)=g^{\prime\prime}(\Sigma)\left(\sum_{i=1}^{n}\frac{h^{\prime}(c_{1})}{h^{\prime}(c_{i})}\right),\text{
or}$ (11)
$\displaystyle\frac{g^{\prime}(\Sigma)}{g^{\prime\prime}(\Sigma)}=\frac{\sum_{i=1}^{n}\frac{h^{\prime}(c_{1})}{h^{\prime}(c_{i})}}{\left(\frac{1}{h^{\prime}}\right)^{\prime}(c_{1})},\text{
or}$ (12) $\displaystyle(\log
g^{\prime})^{-1}=\frac{h^{\prime}(c_{1})}{-\frac{h^{\prime\prime}(c_{1})}{(h^{\prime}(c_{1}))^{2}}}\left(\sum_{i=1}^{n}\frac{1}{h^{\prime}(c_{i})}\right),\text{
or}$ (13) $\displaystyle(\log
g^{\prime})^{-1}=\left[-\frac{(h^{\prime})^{3}}{h^{\prime\prime}}(c_{1})\right]\left(\sum_{i=1}^{n}\frac{1}{h^{\prime}(c_{i})}\right)$

The factor in brackets depends only on the chosen variable, $c_{1}$, whereas
$\log(g^{\prime})(\sum_{i=1}^{n}\frac{1}{h^{\prime}(c_{i})})$ is symmetric in
all $n$ variables. So assuming $n>1$, there are more than a single citizen in
the society, the factor in brackets must be constant. This leads to the
equation:

(14) $(h^{\prime})^{3}=-kh^{\prime\prime}$

for some constant $k$.

Equation 14 may be solved recursively in a Taylor series expansion around any
positive value of the variable, yielding:

(15) $h(y)=ay^{\frac{1}{2}}+b,\text{ and }b=0\text{ by the boundary
condition.}$

Thus, $\frac{1}{h^{\prime}(c_{i})}=\frac{2}{a}y^{\frac{1}{2}}$, so
$\sum_{i=1}^{n}\frac{1}{h^{\prime}(c_{i})}=\frac{a^{2}}{2}\Sigma$, and (13)
becomes:

(16) $\frac{g^{\prime}(\Sigma)}{g^{\prime\prime}(\Sigma)}=\mathrm{const.}\
\Sigma$

which again by solving the Taylor series must be

(17) $g(x)=\mathrm{const.}\ x^{2}+\mathrm{const.}^{\prime}$

with $\mathrm{const.}^{\prime}=0$ by the boundary condition, and const. $=1$
to match self-funding in the limit of one positive contribution. So we have
recovered the $g$ and $h$ of QF:

(18) $g(x)=x^{2}\text{ and }h(y)=y^{\frac{1}{2}}$

To complete the general picture of deriving QF by maximizing social utility,
we can rederive the main conclusion (18) more quickly by _assuming_
homogeneity. If one presumes in advance that $g$ and $h$ are homogeneous,
$g(\Sigma)=\Sigma^{q}$ and $h(c)=c^{p}$, scale-invariance of currency implies
$p=\frac{1}{q}$. Now consider the individual utility at the maximum for each
$i$,

(19) $\displaystyle
U_{i}=V_{i}\left(\sum_{j}c_{j}^{\frac{1}{q}}\right)^{q}-c_{i}$ (20)
$\displaystyle\frac{qV_{i}^{\prime}(\sum_{j}c_{j}^{\frac{1}{q}})^{q-1}}{qc_{i}\frac{q-1}{q}}=1,\text{
or
}V_{i}^{\prime}=\frac{c_{i}^{\frac{q-1}{q}}}{\left(\sum_{j=1}^{n}c_{j}^{\frac{1}{q}}\right)^{1-q}}\text{
at criticality}$

The optimality condition requires:

(21)
$1=\sum_{i=1}^{n}V_{i}^{\prime}=\frac{\sum_{i=1}^{n}c_{i}^{\frac{q-1}{q}}}{\left(\sum_{j=1}^{n}c_{j}^{\frac{1}{q}}\right)^{q-1}}$

By Hölder’s inequality, except when all $c_{j}$ are equal, (20) cannot hold
except at $q=2$, again yielding

(22) $g(x)=x^{2}\text{ and }h(y)=y^{\frac{1}{2}}$

Indeed, there are many roads to Quadratic Finance, but Kant’s categorical
imperative is the most exciting—to your expositor, to Spinoza, to Leibnitz, to
Kant, and to my father.

Our secular rabbi, priest, minister, and iman has done his work well and found
a convert. And I doubt that the modifier _secular_ is even germane. As Jerry
Muller has shown us, we all lie within a single broad and evolving tradition
of thought. We can use both reason and love, faith for those who possess it,
and the calculus of Newton and Leibnitz, for those who possess that, to
navigate to a fairer, less contentious world.

## Acknowledgment

I thank the Aspen Center for Physics for their hospitality.

## References