\typearea

Double shrinkage priors for a normal mean matrix

Takeru Matsuda The University of Tokyo RIKEN Center for Brain Science Fumiyasu Komaki The University of Tokyo RIKEN Center for Brain Science William E. Strawderman Rutgers University

Abstract

We consider estimation of a normal mean matrix under the Frobenius loss. Motivated by the Efron–Morris estimator, a generalization of Stein’s prior has been recently developed, which is superharmonic and shrinks the singular values towards zero. The generalized Bayes estimator with respect to this prior is minimax and dominates the maximum likelihood estimator. However, here we show that it is inadmissible by using Brown’s condition. Then, we develop two types of priors that provide improved generalized Bayes estimators and examine their performance numerically. The proposed priors attain risk reduction by adding scalar shrinkage or column-wise shrinkage to singular value shrinkage. Parallel results for Bayesian predictive densities are also given.

1 Introduction

Suppose that we have independent matrix observations $Y^{(1)},\dots,Y^{(N)}\in\mathbb{R}^{n\times p}$ whose entries are independent normal random variables $Y^{(t)}_{ij}\sim{\rm N}(M_{ij},1)$ , where $M\in\mathbb{R}^{n\times p}$ is an unknown mean matrix. In the notation of Gupta and Nagar (2000), this is expressed as $Y^{(t)}\sim{\rm N}_{n,p}(M,I_{n},I_{p})$ for $t=1,\dots,N$ , where $I_{k}$ denotes the $k$ -dimensional identity matrix. We consider estimation of $M$ under the Frobenius loss

\displaystyle l(M,\hat{M})=\|\hat{M}-M\|_{{\rm F}}^{2}=\sum_{a=1}^{n}\sum_{i=1}^{p}(\hat{M}_{ai}-M_{ai})^{2}.

By sufficiency reduction, it suffices to consider the average ${Y}=(Y^{(1)}+\dots+Y^{(N)})/N\sim{\rm N}(M,I_{n},N^{-1}I_{p})$ in estimation of $M$ . We assume $n-p-1>0$ in the following. Note that vectorization reduces this problem to estimation of a normal mean vector ${\rm vec}(M)$ from ${\rm vec}({Y})\sim{\rm N}_{np}({\rm vec}(M),N^{-1}I_{np})$ under the quadratic loss, which has been well studied in shrinkage estimation theory (Fourdrinier et al., 2018).

Efron and Morris (1972) proposed an empirical Bayes estimator:

\displaystyle\hat{M}_{{\rm EM}}=Y\left(I_{p}-\frac{n-p-1}{N}(Y^{\top}Y)^{-1}\right).

(1)

This estimator can be viewed as a generalization of the James–Stein estimator ( $p=1$ ) for a normal mean vector. Efron and Morris (1972) showed that $\hat{M}_{{\rm EM}}$ is minimax and dominates the maximum likelihood estimator $\hat{M}=Y$ under the Frobenius loss. Let $Y=U\Lambda V^{\top}$ , $U\in\mathbb{R}^{n\times p}$ , $V\in\mathbb{R}^{p\times p}$ , $\Lambda={\rm diag}(\sigma_{1}(Y),\ldots,\sigma_{p}(Y))$ be the singular value decomposition of $Y$ , where $U^{\top}U=V^{\top}V=I_{p}$ and $\sigma_{1}(Y)\geq\cdots\geq\sigma_{p}(Y)\geq 0$ are the singular values of $Y$ . Stein (1974) pointed out that $\hat{M}_{{\rm EM}}$ does not change the singular vectors but shrinks the singular values of $Y$ towards zero:

\displaystyle\hat{M}_{{\rm EM}}=U\hat{\Lambda}_{{\rm EM}}V^{\top},\quad\hat{\Lambda}_{{\rm EM}}={\rm diag}(\sigma_{1}(\hat{M}_{{\rm EM}}),\ldots,\sigma_{p}(\hat{M}_{{\rm EM}})),

where

\displaystyle\sigma_{i}(\hat{M}_{{\rm EM}})=\left(1-\frac{n-p-1}{N\sigma_{i}(Y)^{2}}\right)\sigma_{i}(Y),\quad i=1,\ldots,p.

See Tsukuma and Kubokawa (2020); Yuasa and Kubokawa (2023a, b) for recent developments around the Efron–Morris estimator.

As a Bayesian counterpart of $\hat{M}_{{\rm EM}}$ , Matsuda and Komaki (2015) proposed a singular value shrinkage prior

\displaystyle\pi_{{\rm SVS}}(M)=\det(M^{\top}M)^{-(n-p-1)/2},

(2)

and showed that the generalized Bayes estimator $\hat{M}_{{\rm SVS}}$ with respect to $\pi_{{\rm SVS}}$ dominates the maximum likelihood estimator $\hat{M}=Y$ under the Frobenius loss. This prior can be viewed as a generalization of Stein’s prior $\pi(\mu)=\|\mu\|^{2-n}$ for a normal mean vector $\mu$ ( $p=1$ ) by Stein (1974). Similarly to $\hat{M}_{{\rm EM}}$ in (1) , $\hat{M}_{{\rm SVS}}$ shrinks the singular values towards zero. Thus, it works well when the true matrix is close to low-rank. See Matsuda and Strawderman (2022) and Matsuda (2023) for details on the risk behavior of $\hat{M}_{{\rm EM}}$ and $\hat{M}_{{\rm SVS}}$ .

In this paper, we show that the generalized Bayes estimator with respect to the singular value shrinkage prior $\pi_{{\rm SVS}}$ in (2) is inadmissible under the Frobenius loss. Then, we develop two types of priors that provide improved generalized Bayes estimators asymptotically. The first type adds scalar shrinkage while the second type adds column-wise shrinkage. We conduct numerical experiments and confirm the effectiveness of the proposed priors in finite samples. We also provide parallel results for Bayesian prediction as well as a similar improvement of the blockwise Stein prior, which was conjectured by Brown and Zhao (2009).

This paper is organized as follows. In Section 2, we prove the inadmissibility of the generalized Bayes estimator with respect to the singular value shrinkage prior $\pi_{{\rm SVS}}$ in (2). In Sections 3 and 4, we provide two types of priors that asymptotically dominate the singular value shrinkage prior $\pi_{{\rm SVS}}$ in (2) by adding scalar or column-wise shrinkage, respectively. Numerical results are also given. In Section 5, we provide parallel results for Bayesian prediction. Technical details and similar results for the blockwise Stein prior are given in the Appendix.

2 Inadmissibility of the singular value shrinkage prior

Here, we show that the generalized Bayes estimator with respect to the singular value shrinkage prior $\pi_{{\rm SVS}}$ in (2) is inadmissible under the Frobenius loss. Since $N$ does not affect admissibility results, we fix $N=1$ for convenience in this section.

For estimation of a normal mean vector under the quadratic loss, Brown (1971) derived the following sufficient condition for inadmissibility of generalized Bayes estimators.

Lemma 2.1.

(Brown, 1971) In estimation of $\theta$ from $Y\sim\mathrm{N}_{d}(\theta,I_{d})$ under the quadratic loss, the generalized Bayes estimator of $\theta$ with respect to a prior $\pi(\theta)$ is inadmissible if

\int_{c}^{\infty}{r^{1-d}\underline{m}(r)}{\rm d}r<\infty

for some $c>0$ , where

\underline{m}(r)=\int\frac{1}{m_{\pi}(y)}{\rm d}U_{r}(y),

m_{\pi}(y)=\int p(y\mid\theta)\pi(\theta){\rm d}\theta,

p(y\mid\theta)=\frac{1}{(2\pi)^{d/2}}\exp\left(-\frac{\|y-\theta\|^{2}}{2}\right),

and $U_{r}$ is the uniform measure on the sphere of radius $r$ in $\mathbb{R}^{d}$ .

After vectorization, estimation of a normal mean matrix $M$ from $Y\sim{\rm N}_{n,p}(M,I_{n},I_{p})$ under the Frobenius loss reduces to estimation of a normal mean vector ${\rm vec}(M)$ from ${\rm vec}(Y)\sim{\rm N}_{np}({\rm vec}(M),I_{np})$ under the quadratic loss. Then, by using Brown’s condition in Lemma 2.1, we obtain the following.

Theorem 2.1.

When $p\geq 2$ , the generalized Bayes estimator with respect to $\pi_{{\rm SVS}}$ in (2) is inadmissible under the Frobenius loss.

Proof.

From $n-p-1>0$ and the AM-GM inequality

\displaystyle\left(\prod_{i=1}^{p}\sigma_{i}(M)^{2}\right)^{1/p}\leq\frac{1}{p}\sum_{i=1}^{p}\sigma_{i}(M)^{2},

we have

	$\displaystyle\pi_{{\rm SVS}}(M)=\left(\prod_{i=1}^{p}\sigma_{i}(M)^{2}\right)^{-(n-p-1)/2}$	$\displaystyle\geq\left(\frac{1}{p}\sum_{i=1}^{p}\sigma_{i}(M)^{2}\right)^{-p(n-p-1)/2}$
		$\displaystyle=A_{n,p}\\|M\\|_{{\rm F}}^{-p(n-p-1)},$

where $A_{n,p}=p^{p(n-p-1)/2}$ . Therefore,

$\displaystyle m_{{\rm SVS}}(Y)$	$\displaystyle=\int\pi_{{\rm SVS}}(M)p(Y\mid M){\rm d}M$
	$\displaystyle\geq A_{n,p}\int\\|M\\|_{{\rm F}}^{-p(n-p-1)}p(Y\mid M){\rm d}M$
	$\displaystyle=A_{n,p}{\rm E}[\\|Y+Z\\|_{{\rm F}}^{-p(n-p-1)}]$
	$\displaystyle\geq A_{n,p}{\rm E}[(\\|Y\\|_{{\rm F}}+\\|Z\\|_{{\rm F}})^{-p(n-p-1)}],$	(3)

where $Z=M-Y\sim{\rm N}_{n,p}(O,I_{n},I_{p})$ and we used the triangle inequality. As $\|Y\|_{{\rm F}}\to\infty$ ,

\displaystyle\|Y\|_{{\rm F}}^{p(n-p-1)}{\rm E}[(\|Y\|_{{\rm F}}+\|Z\|_{{\rm F}})^{-p(n-p-1)}]={\rm E}\left[\left(1+\frac{\|Z\|_{{\rm F}}}{\|Y\|_{{\rm F}}}\right)^{-p(n-p-1)}\right]\to 1,

which yields

\displaystyle{\rm E}[(\|Y\|_{{\rm F}}+\|Z\|_{{\rm F}})^{-p(n-p-1)}]=O(\|Y\|_{{\rm F}}^{-p(n-p-1)}).

(4)

Now, we apply Lemma 2.1 by noting that estimation of a normal mean matrix $M$ from $Y\sim{\rm N}_{n,p}(M,I_{n},I_{p})$ under the Frobenius loss is equivalent to estimation of a normal mean vector ${\rm vec}(M)$ from ${\rm vec}(Y)\sim{\rm N}_{np}({\rm vec}(M),I_{np})$ under the quadratic loss. Let $U_{r}$ be the uniform measure on the sphere of radius $r$ in $\mathbb{R}^{n\times p}$ , where the Frobenius norm is adopted for radius. Then, from (3) and (4),

\displaystyle\underline{m}_{{\rm SVS}}(r)

\displaystyle=\int\frac{1}{m_{{\rm SVS}}(Y)}{\rm d}U_{r}(Y)\leq Cr^{p(n-p-1)}

for some constant $C$ . Therefore, since $-p^{2}-p+1<-1$ when $p\geq 2$ ,

\int_{1}^{\infty}r^{1-np}\underline{m}_{{\rm SVS}}(r){\rm d}r\leq C\int_{1}^{\infty}r^{-p^{2}-p+1}{\rm d}r<\infty.

From Lemma 2.1, it implies the inadmissibility of the generelized Bayes estimator with respect to $\pi_{{\rm SVS}}$ under the Frobenius loss. ∎

3 Improvement by additional scalar shrinkage

Here, motivated by the result of Efron and Morris (1976), we develop a class of priors for which the generalized Bayes estimators asymptotically dominate that with respect to the singular value shrinkage prior $\pi_{{\rm SVS}}$ in (2). Efron and Morris (1976) proved that the estimator

\displaystyle\hat{M}_{{\rm MEM}}=Y\left(I_{p}-\frac{n-p-1}{N}(Y^{\top}Y)^{-1}-\frac{p^{2}+p-2}{N\|Y\|_{\mathrm{F}}^{2}}I_{p}\right)

(5)

dominates $\hat{M}_{{\rm EM}}$ in (1) under the Frobenius loss. This estimator shrinks the singular values of $Y$ more strongly than $\hat{M}_{{\rm EM}}$ :

\displaystyle\hat{M}_{{\rm MEM}}=U\hat{\Lambda}_{{\rm MEM}}V^{\top},\quad\hat{\Lambda}_{{\rm MEM}}={\rm diag}(\sigma_{1}(\hat{M}_{{\rm MEM}}),\ldots,\sigma_{p}(\hat{M}_{{\rm MEM}})),

where

\displaystyle\sigma_{i}(\hat{M}_{{\rm MEM}})=\left(1-\frac{n-p-1}{N\sigma_{i}(Y)^{2}}-\frac{p^{2}+p-2}{N\|Y\|_{\mathrm{F}}^{2}}\right)\sigma_{i}(Y),\quad i=1,\ldots,p.

(6)

In other words, $\hat{M}_{{\rm MEM}}$ adds scalar shrinkage to $\hat{M}_{{\rm EM}}$ . Konno (1990, 1991) showed corresponding results in the unknown covariance setting. By extending these results, Tsukuma and Kubokawa (2007) derived a general method for improving matrix mean estimators by adding scalar shrinkage.

Motivated by $\hat{M}_{{\rm MEM}}$ in (5), we construct priors by adding scalar shrinkage to $\pi_{{\rm SVS}}$ in (2):

\displaystyle\pi_{{\rm MSVS1}}(M)=\pi_{{\rm SVS}}(M)\|M\|_{{\rm F}}^{-\gamma},

(7)

where $\gamma\geq 0$ . Note that Tsukuma and Kubokawa (2017) studied this type of prior in the context of Bayesian prediction. Let

\displaystyle m_{{\rm MSVS1}}(Y)=\int p(Y\mid M)\pi_{{\rm MSVS1}}(M){\rm d}M

be the marginal density of $Y$ under the prior $\pi_{{\rm MSVS1}}(M)$ .

Lemma 3.1.

If $0\leq\gamma<p^{2}+p$ , then $m_{{\rm MSVS1}}(Y)<\infty$ for every $Y$ .

Proof.

Since $m_{{\rm MSVS1}}(Y)$ is interpreted as the expectation of $\pi_{{\rm MSVS1}}(M)$ under $M\sim{\rm N}_{n,p}(Y,I_{n},I_{p})$ , it suffices to show that $\pi_{{\rm MSVS1}}(M)$ is locally integrable at every $M$ .

First, consider $M\neq O$ . Since

\displaystyle m_{{\rm SVS}}(Y)

\displaystyle=\int\pi_{{\rm SVS}}(M)p(Y\mid M){\rm d}M<\infty

for every $Y$ from Lemma 1 of Matsuda and Komaki (2015), $\pi_{{\rm SVS}}(M)$ is locally integrable at $M$ . Also, $\|M\|_{{\rm F}}>c$ for some $c>0$ in a neighborhood of $M$ . Thus, $\pi_{{\rm MSVS1}}(M)=\pi_{{\rm SVS}}(M)\|M\|_{{\rm F}}^{-\gamma}$ is locally integrable at $M$ if $\gamma\geq 0$ .

Next, consider $M=O$ and take its neighborhood $A=\{M\mid\|M\|_{\mathrm{F}}\leq\varepsilon\}$ for $\varepsilon>0$ . To evaluate the integral on $A$ , we use the variable transformation from $M$ to $(r,U)$ , where $r=\|M\|_{\mathrm{F}}$ and $U=M/r$ so that $M=rU$ . We have ${\rm d}M=r^{np-1}{\rm d}r{\rm d}U$ . Also, from $\det(M^{\top}M)=r^{2p}\det(U^{\top}U)$ ,

\displaystyle\pi_{{\rm MSVS1}}(M)=r^{-p(n-p-1)-\gamma}\det(U^{\top}U)^{-(n-p-1)/2}.

Thus,

		$\displaystyle\int_{A}\pi_{{\rm MSVS1}}(M){\rm d}M$
	$\displaystyle=$	$\displaystyle\int_{0}^{\varepsilon}r^{-p(n-p-1)-\gamma+np-1}{\rm d}r\int\det(U^{\top}U)^{-(n-p-1)/2}{\rm d}U$
	$\displaystyle=$	$\displaystyle\int_{0}^{\varepsilon}r^{p^{2}+p-\gamma-1}{\rm d}r\int\det(U^{\top}U)^{-(n-p-1)/2}{\rm d}U.$

The integral with respect to $r$ is finite if $p^{2}+p-\gamma-1>-1$ , which is equivalent to $\gamma<p^{2}+p$ . The integral with respect to $U$ is finite due to the local integrability of $\pi_{\mathrm{SVS}}$ , which corresponds to $\gamma=0$ , at $M=O$ . Therefore, $\pi_{{\rm MSVS1}}(M)$ is locally integrable at $M=O$ if $0\leq\gamma<p^{2}+p$ .

Hence, $\pi_{{\rm MSVS1}}(M)$ is locally integrable at every $M$ if $0\leq\gamma<p^{2}+p$ . ∎

From Lemma 3.1, the generalized Bayes estimator with respect to $\pi_{{\rm MSVS1}}$ is well-defined when $0\leq\gamma<p^{2}+p$ . We denote it by $\hat{M}_{{\rm MSVS1}}$ .

Theorem 3.1.

For every $M$ ,

\displaystyle N^{2}({\rm E}_{M}[\|\hat{M}_{{{\rm MSVS1}}}-M\|_{{\rm F}}^{2}]-{\rm E}_{M}[\|\hat{M}_{{{\rm SVS}}}-M\|_{{\rm F}}^{2}])\to\frac{\gamma(\gamma-2p^{2}-2p+4)}{{\rm tr}(M^{\top}M)}

(8)

as $N\to\infty$ . Therefore, if $p\geq 2$ and $0<\gamma<p^{2}+p$ , then the generalized Bayes estimator with respect to $\pi_{{\rm MSVS1}}$ in (7) asymptotically dominates that with respect to $\pi_{{\rm SVS}}$ in (2) under the Frobenius loss.

Proof.

Let $K=M^{\top}M$ and $K^{ij}$ be the $(i,j)$ th entry of $K^{-1}$ . From

\displaystyle\frac{\partial K_{jk}}{\partial M_{ai}}=\delta_{ik}M_{aj}+\delta_{ij}M_{ak},\quad\frac{\partial}{\partial K_{ij}}\det K=K^{ij}\det K,

(9)

we have

\displaystyle\frac{\partial}{\partial M_{ai}}\det K

\displaystyle=\sum_{j,k}\frac{\partial K_{jk}}{\partial M_{ai}}\frac{\partial}{\partial K_{jk}}\det K=2\sum_{j}M_{aj}K^{ij}\det K.

Therefore,

\displaystyle\frac{\partial}{\partial M_{ai}}\log\pi_{{\rm SVS}}(M)=-(n-p-1)\sum_{j}M_{aj}K^{ij}.

(10)

Let $\pi_{{\rm S}}(M)=\|M\|_{{\rm F}}^{-\gamma}=({\rm tr}K)^{-\gamma/2}$ . Since

\displaystyle\frac{\partial}{\partial M_{ai}}{\rm tr}K=2M_{ai}

from (9), we have

	$\displaystyle\frac{\partial}{\partial M_{ai}}\log\pi_{{\rm S}}(M)$	$\displaystyle=-\gamma M_{ai}({\rm tr}K)^{-1},$		(11)
	$\displaystyle\frac{\partial^{2}}{\partial M_{ai}^{2}}\log\pi_{{\rm S}}(M)$	$\displaystyle=-\gamma({\rm tr}K-2M_{ai}^{2})({\rm tr}K)^{-2}.$		(12)

By using (10), (11), and (12), we obtain

	$\displaystyle{\rm tr}(\widetilde{\nabla}\log\pi_{{\rm SVS}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm S}}(M))$	$\displaystyle=\gamma p(n-p-1)({\rm tr}K)^{-1},$
	$\displaystyle{\rm tr}(\widetilde{\nabla}\log\pi_{{\rm S}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm S}}(M))$	$\displaystyle=\gamma^{2}({\rm tr}K)^{-1},$
	$\displaystyle{\rm tr}(\widetilde{\Delta}\log\pi_{{\rm S}}(M))$	$\displaystyle=-\gamma(np-2)({\rm tr}K)^{-1},$

where we used the matrix derivative notations (28) and (29). Therefore, from Lemma A.2,

		$\displaystyle{\rm E}_{M}[\\|\hat{M}_{{\rm MSVS1}}-M\\|_{{\rm F}}^{2}]-{\rm E}_{M}[\\|\hat{M}_{{\rm SVS}}-M\\|_{{\rm F}}^{2}]$
	$\displaystyle=$	$\displaystyle\frac{1}{N^{2}}{\rm tr}(2\widetilde{\nabla}\log\pi_{{\rm SVS}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm S}}(M)+\widetilde{\nabla}\log\pi_{{\rm S}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm S}}(M)+2\widetilde{\Delta}\log\pi_{{\rm S}}(M))$
		$\displaystyle\quad+o(N^{-2})$
	$\displaystyle=$	$\displaystyle\frac{1}{N^{2}}\gamma(\gamma-2p^{2}-2p+4)({\rm tr}K)^{-1}+o(N^{-2}).$

Hence, we obtain (8). ∎

From (8), the choice $\gamma=p^{2}+p-2$ attains the minimum risk among $0<\gamma<p^{2}+p$ . Note that $p^{2}+p-2$ also appears in the singular value decomposition form of the modified Efron–Morris estimator $\hat{M}_{{\rm MEM}}$ in (6).

Now, we examine the performance of $\pi_{{\rm MSVS1}}$ in (7) by Monte Carlo simulation. Figure 1 plots the Frobenius risk of generalized Bayes estimators with respect to $\pi_{{\rm MSVS1}}$ in (7) with $\gamma=p^{2}+p-2$ , $\pi_{{\rm SVS}}$ in (2) and $\pi_{{\rm S}}(M)=\|M\|_{{\rm F}}^{2-np}$ , which is Stein’s prior on the vectorization of $M$ , for $n=10$ , $p=3$ and $N=1$ . We computed the generalized Bayes estimators by using the random-walk Metropolis–Hastings algorithm with Gaussian proposal of variance $0.1$ . Note that the Frobenius risk of these estimators depends only on the singular values of $M$ due to the orthogonal invariance. Similarly to the Efron–Morris estimator and $\pi_{{\rm SVS}}$ , $\pi_{{\rm MSVS1}}$ works well when $M$ is close to low-rank. Also, $\pi_{{\rm MSVS1}}$ attains large risk reduction when $M$ is close to the zero matrix like $\pi_{{\rm S}}$ . Thus, $\pi_{{\rm MSVS1}}$ has the best of both worlds. Figure 2 plots the Frobenius risk for $n=10$ , $p=3$ and $N=10$ , computed by the random walk Metropolis–Hastings algorithm with proposal variance $0.005$ . The risk behavior is similar to Figure 1. Figure 3 plots the Frobenius risk for $n=20$ , $p=3$ and $N=2$ , computed by the random walk Metropolis–Hastings algorithm with proposal variance $0.01$ . Again, the risk behavior is similar to Figure 1. Note that the value of $np/N=30$ is the same with Figure 1.

Figure 1: Frobenius risk of generalized Bayes estimators for

n=10

p=3

and

N=1

. Left:

\sigma_{2}=\sigma_{3}=0

. Right:

\sigma_{1}=10

\sigma_{3}=0

. solid:

\pi_{{\rm MSVS1}}

with

\gamma=p^{2}+p-2

, dashed:

\pi_{{\rm SVS}}

, dotted: Stein’s prior

\pi_{{\rm S}}

, Note that the minimax risk is

np/N=30

Figure 2: Frobenius risk of generalized Bayes estimators for

n=10

p=3

and

N=10

. Left:

\sigma_{2}=\sigma_{3}=0

. Right:

\sigma_{1}=10

\sigma_{3}=0

. solid:

\pi_{{\rm MSVS1}}

with

\gamma=p^{2}+p-2

, dashed:

\pi_{{\rm SVS}}

, dotted: Stein’s prior

\pi_{{\rm S}}

, Note that the minimax risk is

np/N=3

Figure 3: Frobenius risk of generalized Bayes estimators for

n=20

p=3

and

N=2

. Left:

\sigma_{2}=\sigma_{3}=0

. Right:

\sigma_{1}=10

\sigma_{3}=0

. solid:

\pi_{{\rm MSVS1}}

with

\gamma=p^{2}+p-2

, dashed:

\pi_{{\rm SVS}}

, dotted: Stein’s prior

\pi_{{\rm S}}

, Note that the minimax risk is

np/N=30

Improvement by additional scalar shrinkage holds even under the matrix quadratic loss (Matsuda and Strawderman, 2022; Matsuda, 2024).

Theorem 3.2.

For every $M$ ,

	$\displaystyle N^{2}({\rm E}_{M}[(\hat{M}_{{{\rm MSVS1}}}-M)^{\top}$	$\displaystyle(\hat{M}_{{{\rm MSVS1}}}-M)]-{\rm E}_{M}[(\hat{M}_{{{\rm SVS}}}-M)^{\top}(\hat{M}_{{{\rm SVS}}}-M)])$
		$\displaystyle\to\gamma({\rm tr}K)^{-2}(-2(p+1)({\rm tr}K)I_{p}+(\gamma+4)K)$		(13)

as $N\to\infty$ . Therefore, if $p\geq 2$ and $0<\gamma<2p-2$ , then the generalized Bayes estimator with respect to $\pi_{{\rm MSVS1}}$ in (7) asymptotically dominates that with respect to $\pi_{{\rm SVS}}$ in (2) under the matrix quadratic loss.

Proof.

We use the same notation with the proof of Theorem 3.1. By using (10), (11), and (12), we obtain

	$\displaystyle(\widetilde{\nabla}\log\pi_{{\rm SVS}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm S}}(M))$	$\displaystyle=\gamma(n-p-1)({\rm tr}K)^{-1}I_{p},$
	$\displaystyle(\widetilde{\nabla}\log\pi_{{\rm S}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm S}}(M))$	$\displaystyle=\gamma^{2}({\rm tr}K)^{-2}K,$
	$\displaystyle(\widetilde{\Delta}\log\pi_{{\rm S}}(M))$	$\displaystyle=-n\gamma({\rm tr}K)^{-1}I_{p}+2\gamma({\rm tr}K)^{-2}K.$

Therefore, from Lemma A.3,

		$\displaystyle{\rm E}_{M}[(\hat{M}_{{{\rm MSVS1}}}-M)^{\top}(\hat{M}_{{{\rm MSVS1}}}-M)]-{\rm E}_{M}[(\hat{M}_{{{\rm SVS}}}-M)^{\top}(\hat{M}_{{{\rm SVS}}}-M)]$
	$\displaystyle=$	$\displaystyle\frac{1}{N^{2}}(2\widetilde{\nabla}\log\pi_{{\rm SVS}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm S}}(M)+\widetilde{\nabla}\log\pi_{{\rm S}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm S}}(M)+2\widetilde{\Delta}\log\pi_{{\rm S}}(M))$
		$\displaystyle\quad+o(N^{-2})$
	$\displaystyle=$	$\displaystyle\frac{1}{N^{2}}\gamma({\rm tr}K)^{-2}(-2(p+1)({\rm tr}K)I_{p}+(\gamma+4)K)+o(N^{-2}).$

Hence, we obtain (13). Since $K\preceq({\rm tr}K)I_{p}$ from $K\succeq O$ ,

\displaystyle-2(p+1)({\rm tr}K)I_{p}+(\gamma+4)K\preceq(\gamma-2p+2)({\rm tr}K)I_{p}\prec O

if $0<\gamma<2p-2$ . ∎

The generalized Bayes estimator with respect to $\pi_{{\rm MSVS1}}$ in (7) attains minimaxity in some cases as follows.

Theorem 3.3.

If $p\geq 2$ , $p+2\leq n<2p+2-2/p$ and $0<\gamma\leq-np+2p^{2}+2p-2$ , then the generalized Bayes estimator with respect to $\pi_{{\rm MSVS1}}$ in (7) is minimax under the Frobenius loss.

Proof.

From Proposition B.1,

\displaystyle\Delta\pi_{{\rm MSVS1}}(M)

\displaystyle=\gamma(\gamma+np-2p^{2}-2p+2)\|M\|_{{\rm F}}^{-2}\pi_{{\rm MSVS1}}(M)\leq 0.

Thus, $\pi_{{\rm MSVS1}}(M)$ is superharmonic, which indicates the minimaxity of the generalized Bayes estimator with respect to $\pi_{{\rm MSVS1}}$ in (7) under the Frobenius loss from Stein’s classical result (Stein, 1974; Matsuda and Komaki, 2015). ∎

It is an interesting problem whether the generalized Bayes estimator with respect to $\pi_{{\rm MSVS1}}$ in (7) attains admissibility or not. In addition to Lemma 2.1, Brown (1971) derived the following sufficient condition for admissibility of generalized Bayes estimators, which may be useful here. While the condition (14) can be verified by using a similar argument to Theorem 2.1, the verification of the uniform boundedness of $\|\nabla\log m_{\pi}(y)\|$ seems difficult. We leave further investigation for future work.

Lemma 3.2.

\displaystyle\int_{c}^{\infty}\frac{r^{1-d}}{\accentset{{\cc@style\underline{\mskip 10.0mu}}}{m}(r)}{\rm d}r=\infty,

(14)

where

\accentset{{\cc@style\underline{\mskip 10.0mu}}}{m}(r)=\int{m_{\pi}(y)}{\rm d}U_{r}(y)

and $U_{r}$ is the uniform measure on the sphere of radius $r$ in $\mathbb{R}^{d}$ .

4 Improvement by additional column-wise shrinkage

Here, instead of scalar shrinkage, we consider priors with additional column-wise shrinkage:

\displaystyle\pi_{{\rm MSVS2}}(M)=\pi_{{\rm SVS}}(M)\prod_{i=1}^{p}\|M_{\cdot i}\|^{-\gamma_{i}},

(15)

where $\gamma_{i}\geq 0$ for every $i$ and $\|M_{\cdot i}\|$ denotes the norm of the $i$ -th column vector of $M$ . Let

\displaystyle m_{{\rm MSVS2}}(Y)=\int p(Y\mid M)\pi_{{\rm MSVS2}}(M){\rm d}M

be the marginal density of $Y$ under the prior $\pi_{{\rm MSVS2}}(M)$ .

Lemma 4.1.

If $0\leq\gamma_{i}\leq p$ for every $i$ , then $m_{{\rm MSVS2}}(Y)<\infty$ for every $Y$ .

Proof.

Similarly to Lemma 3.1, it suffices to show that $\pi_{{\rm MSVS2}}(M)$ is locally integrable at $M=O$ . Consider the neighborhood of $M=O$ defined by $A=\{M\mid\|M\|_{\mathrm{F}}\leq\varepsilon\}$ for $\varepsilon>0$ . To evaluate the integral on $A$ , we use the variable transformation from $M$ to $(r_{1},\dots,r_{p},u_{1},\dots,u_{p})$ , where each $r_{i}\in[0,\infty)$ and $u_{i}\in\mathbb{R}^{n}$ with $\|u_{i}\|=1$ are defined by $r_{i}=\|M_{\cdot i}\|$ and $u_{i}=M_{\cdot i}/r_{i}$ for the $i$ -th column vector $M_{\cdot i}$ of $M$ so that $M_{\cdot i}=r_{i}u_{i}$ (polar coordinate). Since ${\rm d}M_{\cdot i}=r_{i}^{n-1}{\rm d}r_{i}{\rm d}u_{i}$ ,

\displaystyle{\rm d}M=r_{1}^{n-1}\dots r_{p}^{n-1}{\rm d}r_{1}\dots\rm dr_{p}{\rm d}u_{1}\dots\rm du_{p}.

Also, from $M^{\top}M=D(U^{\top}U)D$ with $D={\rm diag}(r_{1},\dots,r_{p})$ and $U=(u_{1}\dots u_{p})$ ,

	$\displaystyle\pi_{{\rm MSVS2}}(M)$	$\displaystyle=\det(D)^{-(n-p-1)}\det(U^{\top}U)^{-(n-p-1)/2}\prod_{i=1}^{p}r_{i}^{-\gamma_{i}}$
		$\displaystyle=\prod_{i=1}^{p}r_{i}^{-(n-p-1)-\gamma_{i}}\cdot\det(U^{\top}U)^{-(n-p-1)/2}.$

Thus,

		$\displaystyle\int_{A}\pi_{{\rm MSVS2}}(M){\rm d}M$
	$\displaystyle=$	$\displaystyle\int_{\\|r\\|\leq\varepsilon}\left(\prod_{i=1}^{p}r_{i}^{p-\gamma_{i}}\right){\rm d}r_{1}\dots\rm dr_{p}\int\det(U^{\top}U)^{-(n-p-1)/2}{\rm d}u_{1}\dots\rm du_{p}.$		(16)

By variable transformation from $r=(r_{1},\dots,r_{p})$ to $s=\|r\|$ and $v=r/s$ , the first integral in (16) is reduced to

\displaystyle\int_{0}^{\varepsilon}s^{p^{2}-\sum_{i=1}^{p}\gamma_{i}+p-1}{\rm d}s\int_{\|v\|=1}\left(\prod_{i=1}^{p}v_{i}^{p-\gamma_{i}}\right){\rm d}v.

The integral with respect to $s$ is finite if $p^{2}-\sum_{i=1}^{p}\gamma_{i}+p-1>-1$ , which is equivalent to $\sum_{i=1}^{p}\gamma_{i}<p^{2}+p$ . The integral with respect to $v$ is finite if $p-\gamma_{i}\geq 0$ for every $i$ . On the other hand, the second integral in (16) is finite due to the local integrability of $\pi_{\mathrm{SVS}}$ , which corresponds to $\gamma=0$ , at $M=O$ . Therefore, $\pi_{{\rm MSVS2}}(M)$ is locally integrable at $M=O$ if $0\leq\gamma_{i}\leq p$ for every $i$ . ∎

From Lemma 4.1, the generalized Bayes estimator with respect to $\pi_{{\rm MSVS2}}$ is well-defined when $0\leq\gamma\leq p$ . We denote it by $\hat{M}_{{\rm MSVS2}}$ .

Theorem 4.1.

For every $M$ ,

\displaystyle N^{2}({\rm E}_{M}[\|\hat{M}_{{{\rm MSVS2}}}-M\|_{{\rm F}}^{2}]-{\rm E}_{M}[\|\hat{M}_{{{\rm SVS}}}-M\|_{{\rm F}}^{2}])\to\sum_{i=1}^{p}\gamma_{i}(\gamma_{i}-2p+2)\|M_{\cdot i}\|^{-2}

(17)

as $N\to\infty$ . Therefore, if $p\geq 2$ and $0<\gamma_{i}\leq p$ for every $i$ , then the generalized Bayes estimator with respect to $\pi_{{\rm MSVS2}}$ in (15) asymptotically dominates that with respect to $\pi_{{\rm SVS}}$ in (2) under the Frobenius loss.

Proof.

Let

\displaystyle\pi_{{\rm CS}}(M)=\prod_{i=1}^{p}\|M_{\cdot i}\|^{-\gamma_{i}}.

Then,

\displaystyle\frac{\partial}{\partial M_{ai}}\log\pi_{{\rm CS}}(M)=-\gamma_{i}M_{ai}\|M_{\cdot i}\|^{-2},

(18)

\displaystyle\frac{\partial^{2}}{\partial M_{ai}^{2}}\log\pi_{{\rm CS}}(M)=-\gamma_{i}\left(\|M_{\cdot i}\|^{2}-2M_{ai}^{2}\right)\|M_{\cdot i}\|^{-4}.

(19)

From (10), (18), and (19),

	$\displaystyle{\rm tr}(\widetilde{\nabla}\log\pi_{{\rm SVS}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm CS}}(M))=(n-p-1)\sum_{i=1}^{p}\gamma_{i}\\|M_{\cdot i}\\|^{-2},$		(20)
	$\displaystyle{\rm tr}(\widetilde{\nabla}\log\pi_{{\rm CS}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm CS}}(M))=\sum_{i=1}^{p}\gamma_{i}^{2}\\|M_{\cdot i}\\|^{-2},$		(21)
	$\displaystyle{\rm tr}(\widetilde{\Delta}\log\pi_{{\rm CS}}(M))=-(n-2)\sum_{i=1}^{p}\gamma_{i}\\|M_{\cdot i}\\|^{-2},$		(22)

where we used the matrix derivative notations (28) and (29). Therefore, from Lemma A.2,

		$\displaystyle{\rm E}_{M}[\\|\hat{M}_{{\rm MSVS2}}-M\\|_{{\rm F}}^{2}]-{\rm E}_{M}[\\|\hat{M}_{{\rm SVS}}-M\\|_{{\rm F}}^{2}]$
	$\displaystyle=$	$\displaystyle\frac{1}{N^{2}}{\rm tr}(2\widetilde{\nabla}\log\pi_{{\rm SVS}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm CS}}(M)+\widetilde{\nabla}\log\pi_{{\rm CS}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm CS}}(M)+2\widetilde{\Delta}\log\pi_{{\rm CS}}(M))$
		$\displaystyle\quad+o(N^{-2})$
	$\displaystyle=$	$\displaystyle\frac{1}{N^{2}}\sum_{i=1}^{p}\gamma_{i}(\gamma_{i}-2p+2)\\|M_{\cdot i}\\|^{-2}+o(N^{-2}).$

Hence, we obtain (17). ∎

From (17), the choice $\gamma_{1}=\dots=\gamma_{p}=p-1$ attains the minimum risk among $0<\gamma_{i}\leq p$ .

Now, we examine the performance of $\pi_{{\rm MSVS2}}$ in (15) by Monte Carlo simulation. Figure 4 plots the Frobenius risk of generalized Bayes estimators with respect to $\pi_{{\rm MSVS2}}$ in (15) with $\gamma_{1}=\dots=\gamma_{p}=p-1$ , $\pi_{{\rm MSVS1}}$ in (7) with $\gamma=p^{2}+p-2$ , $\pi_{{\rm SVS}}$ in (2) and $\pi_{{\rm S}}(M)=\|M\|_{{\rm F}}^{2-np}$ , which is Stein’s prior on the vectorization of $M$ . We computed the generalized Bayes estimators by using the random walk Metropolis–Hastings algorithm with proposal variance $0.1$ . We set $M=U\Sigma$ , where $U^{\top}U=I_{p}$ and $\Sigma={\rm diag}(\sigma_{1},\dots,\sigma_{p})$ . For this $M$ , the Frobenius risk of the estimators compared here depends only on the singular values $\sigma_{1},\dots,\sigma_{p}$ of $M$ . Overall, $\pi_{{\rm MSVS2}}$ performs better than $\pi_{{\rm SVS}}$ . Also, $\pi_{{\rm MSVS2}}$ even dominates $\pi_{{\rm MSVS1}}$ and $\pi_{{\rm S}}$ except when $\sigma_{1}$ is sufficiently small. This is understood from the column-wise shrinkage effect of $\pi_{{\rm MSVS2}}$ . Figure 5 plots the Frobenius risk for $n=10$ , $p=3$ and $N=10$ , computed by the random walk Metropolis–Hastings algorithm with proposal variance $0.01$ . The risk behavior is similar to Figure 4. Figure 6 plots the Frobenius risk for $n=20$ , $p=3$ and $N=2$ , computed by the random walk Metropolis–Hastings algorithm with proposal variance $0.01$ . Again, the risk behavior is similar to Figure 4. Note that the value of $np/N=30$ is the same with Figure 4.

Figure 4: Frobenius risk of generalized Bayes estimators for

n=10

p=3

and

N=1

where

M=U\Sigma

with

U^{\top}U=I_{p}

and

\Sigma={\rm diag}(\sigma_{1},\dots,\sigma_{p})

. Left:

\sigma_{2}=\sigma_{3}=0

. Right:

\sigma_{1}=10

\sigma_{3}=0

. solid:

\pi_{{\rm MSVS2}}

with

\gamma_{1}=\dots=\gamma_{p}=p-1

, dash-dotted:

\pi_{{\rm MSVS1}}

with

\gamma=p^{2}+p-2

, dashed:

\pi_{{\rm SVS}}

, dotted: Stein’s prior

\pi_{{\rm S}}

. Note that the minimax risk is

np/N=30

Figure 5: Frobenius risk of generalized Bayes estimators for

n=10

p=3

and

N=10

where

M=U\Sigma

with

U^{\top}U=I_{p}

and

\Sigma={\rm diag}(\sigma_{1},\dots,\sigma_{p})

. Left:

\sigma_{2}=\sigma_{3}=0

. Right:

\sigma_{1}=10

\sigma_{3}=0

. solid:

\pi_{{\rm MSVS2}}

with

\gamma_{1}=\dots=\gamma_{p}=p-1

, dash-dotted:

\pi_{{\rm MSVS1}}

with

\gamma=p^{2}+p-2

, dashed:

\pi_{{\rm SVS}}

, dotted: Stein’s prior

\pi_{{\rm S}}

. Note that the minimax risk is

np/N=3

Figure 6: Frobenius risk of generalized Bayes estimators for

n=20

p=3

and

N=2

where

M=U\Sigma

with

U^{\top}U=I_{p}

and

\Sigma={\rm diag}(\sigma_{1},\dots,\sigma_{p})

. Left:

\sigma_{2}=\sigma_{3}=0

. Right:

\sigma_{1}=10

\sigma_{3}=0

. solid:

\pi_{{\rm MSVS2}}

with

\gamma_{1}=\dots=\gamma_{p}=p-1

, dash-dotted:

\pi_{{\rm MSVS1}}

with

\gamma=p^{2}+p-2

, dashed:

\pi_{{\rm SVS}}

, dotted: Stein’s prior

\pi_{{\rm S}}

. Note that the minimax risk is

np/N=30

Improvement by additional column-wise shrinkage holds even under the matrix quadratic loss (Matsuda and Strawderman, 2022; Matsuda, 2024).

Theorem 4.2.

For every $M$ ,

	$\displaystyle N^{2}({\rm E}_{M}[(\hat{M}_{{{\rm MSVS2}}}-M)^{\top}$	$\displaystyle(\hat{M}_{{{\rm MSVS2}}}-M)]-{\rm E}_{M}[(\hat{M}_{{{\rm SVS}}}-M)^{\top}(\hat{M}_{{{\rm SVS}}}-M)])$
		$\displaystyle\to-2(p-1)D+DM^{\top}MD$		(23)

as $N\to\infty$ , where $D={\rm diag}(\gamma_{1}\|M_{\cdot 1}\|^{-2},\dots,\gamma_{p}\|M_{\cdot p}\|^{-2})$ . Therefore, if $p\geq 2$ and $0<\gamma_{1}=\dots=\gamma_{p}<2-2/p$ , then the generalized Bayes estimator with respect to $\pi_{{\rm MSVS2}}$ in (15) asymptotically dominates that with respect to $\pi_{{\rm SVS}}$ in (2) under the matrix quadratic loss.

Proof.

We use the same notation with the proof of Theorem 4.1. By using (10), (18), and (19), we obtain

	$\displaystyle\widetilde{\nabla}\log\pi_{{\rm SVS}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm CS}}(M)=(n-p-1)D,$
	$\displaystyle\widetilde{\nabla}\log\pi_{{\rm CS}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm CS}}(M)=DM^{\top}MD,$
	$\displaystyle\widetilde{\Delta}\log\pi_{{\rm CS}}(M)=-(n-2)D.$

Therefore, from Lemma A.3,

		$\displaystyle{\rm E}_{M}[(\hat{M}_{{{\rm MSVS2}}}-M)^{\top}(\hat{M}_{{{\rm MSVS2}}}-M)]-{\rm E}_{M}[(\hat{M}_{{{\rm SVS}}}-M)^{\top}(\hat{M}_{{{\rm SVS}}}-M)]$
	$\displaystyle=$	$\displaystyle\frac{1}{N^{2}}(2\widetilde{\nabla}\log\pi_{{\rm SVS}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm CS}}(M)+\widetilde{\nabla}\log\pi_{{\rm CS}}(M)^{\top}\widetilde{\nabla}\log\pi_{{\rm CS}}(M)+2\widetilde{\Delta}\log\pi_{{\rm CS}}(M))$
		$\displaystyle\quad+o(N^{-2})$
	$\displaystyle=$	$\displaystyle-2(p-1)D+DM^{\top}MD.$

Hence, we obtain (23).

Suppose that $p\geq 2$ and $0<\gamma_{1}=\dots=\gamma_{p}<2-2/p$ . Let $\|\cdot\|$ be the operater norm. Since $D\succeq O$ is diagonal and $\|M\|^{2}\leq\|M\|_{\mathrm{F}}^{2}=\sum_{i=1}^{p}\|M_{\cdot i}\|^{2}$ ,

	$\displaystyle\\|D^{1/2}M^{\top}MD^{1/2}\\|$	$\displaystyle\leq\\|D^{1/2}\\|\\|M^{\top}M\\|\\|D^{1/2}\\|$
		$\displaystyle=(\max_{i}D_{ii})\\|M\\|^{2}$
		$\displaystyle=\gamma_{1}\frac{\\|M\\|^{2}}{\min_{i}\\|M_{\cdot i}\\|^{2}}$
		$\displaystyle\leq\gamma_{1}\frac{\sum_{i=1}^{p}\\|M_{\cdot i}\\|^{2}}{\min_{i}\\|M_{\cdot i}\\|^{2}}$
		$\displaystyle\leq\gamma_{1}p$
		$\displaystyle<2(p-1),$

which yields $D^{1/2}M^{\top}MD^{1/2}\prec 2(p-1)I_{p}$ . Therefore,

\displaystyle-2(p-1)D+DM^{\top}MD=D^{1/2}(-2(p-1)I_{p}+D^{1/2}M^{\top}MD^{1/2})D^{1/2}\prec O.

∎

The generalized Bayes estimator with respect to $\pi_{{\rm MSVS2}}$ in (15) attains minimaxity in some cases as follows. It is an interesting future work to investigate its admissibility.

Theorem 4.3.

If $p\geq 3$ , $p+2\leq n<2p$ and $0<\gamma\leq-n+2p$ , then the generalized Bayes estimator with respect to $\pi_{{\rm MSVS2}}$ in (15) is minimax under the Frobenius loss.

Proof.

From Proposition B.2,

\displaystyle\Delta\pi_{{\rm MSVS2}}(M)

\displaystyle=\gamma(\gamma+n-2p)\left(\sum_{i=1}^{p}\|M_{\cdot i}\|^{-2}\right)\pi_{{\rm MSVS2}}(M)\leq 0.

Thus, $\pi_{{\rm MSVS2}}(M)$ is superharmonic, which indicates the minimaxity of the generalized Bayes estimator with respect to $\pi_{{\rm MSVS2}}$ in (15) under the Frobenius loss from Stein’s classical result (Stein, 1974; Matsuda and Komaki, 2015). ∎

5 Bayesian prediction

Here, we consider Bayesian prediction and provide parallel results to those in Sections 3 and 4. Suppose that we observe $Y\sim{\rm N}_{n,p}(M,I_{n},N^{-1}I_{p})$ and predict $\widetilde{Y}\sim{\rm N}_{n,p}(M,I_{n},I_{p})$ by a predictive density $\hat{p}(\widetilde{Y}\mid Y)$ . We evaluate predictive densities by the Kullback–Leibler loss:

\displaystyle D({p}(\cdot\mid M),\hat{p}(\cdot\mid Y))=\int{p}(\widetilde{Y}\mid M)\log\frac{{p}(\widetilde{Y}\mid M)}{\hat{p}(\widetilde{Y}\mid Y)}{\rm d}\widetilde{Y}.

The Bayesian predictive density based on a prior $\pi(M)$ is defined as

\displaystyle\hat{p}_{\pi}(\widetilde{Y}\mid Y)

\displaystyle=\int{p}(\widetilde{Y}\mid M)\pi(M\mid Y){\rm d}M,

where $\pi(M\mid Y)$ is the posterior distribution of $M$ given $Y$ , and it minimizes the Bayes risk (Aitchison, 1975):

\displaystyle\hat{p}_{\pi}(\widetilde{Y}\mid Y)

\displaystyle=\operatorname*{arg\,min}_{\hat{p}}\int D({p}(\cdot\mid M),\hat{p}(\cdot\mid Y))p(Y\mid M)\pi(M){\rm d}Y{\rm d}M.

The Bayesian predictive density with respect to the uniform prior is minimax. However, it is inadmissible and dominated by Bayesian predictive densities based on superharmonic priors (Komaki, 2001; George, Liang and Xu, 2006). In particular, the Bayesian predictive density based on the singular value shrinkage prior $\pi_{{\rm SVS}}$ in (2) is minimax and dominates that based on the uniform prior (Matsuda and Komaki, 2015).

The asymptotic expansion of the difference between the Kullback–Leibler risk of two Bayesian predictive densities is obtained as follows.

Lemma 5.1.

As $N\to\infty$ , the difference between the Kullback–Leibler risk of $p_{\pi_{1}}(\widetilde{Y}\mid Y)$ and $p_{\pi_{1}\pi_{2}}(\widetilde{Y}\mid Y)$ is expanded as

	$\displaystyle{\rm E}_{M}[D(p(\widetilde{Y}\mid M),p_{\pi_{1}\pi_{2}}(\widetilde{Y}\mid Y))]-{\rm E}_{M}[D(p(\widetilde{Y}\mid M),p_{\pi_{1}}(\widetilde{Y}\mid Y))]$
$\displaystyle=$	$\displaystyle\frac{1}{2N^{2}}{\rm tr}(2(\widetilde{\nabla}\log\pi_{1}(M))^{\top}(\widetilde{\nabla}\log\pi_{2}(M))+(\widetilde{\nabla}\log\pi_{2}(M))^{\top}(\widetilde{\nabla}\log\pi_{2}(M))+2\widetilde{\Delta}\log\pi_{2}(M))$
	$\displaystyle\quad+o(N^{-2}).$	(24)

Proof.

For the normal model with known covariance, the information geometrical quantities (Amari, 1985) are given by

\displaystyle g_{ij}=g^{ij}=\delta_{ij},\quad\Gamma_{ij}^{k}=0,\quad T_{ijk}=0.

Also, the Jeffreys prior coincides with the uniform prior $\pi(M)\equiv 1$ . Therefore, from equation (3) of Komaki (2006), the Kullback–Leibler risk of the Bayesian predictive density $p_{\pi}(\widetilde{Y}\mid Y)$ based on a prior $\pi(M)$ is expanded as

		$\displaystyle{\rm E}_{M}[D(p(\widetilde{Y}\mid M),p_{\pi}(\widetilde{Y}\mid Y))]$
	$\displaystyle=$	$\displaystyle\frac{np}{2N}+\frac{1}{2N^{2}}{\rm tr}((\widetilde{\nabla}\log\pi(M))^{\top}(\widetilde{\nabla}\log\pi(M))+2\widetilde{\Delta}\log\pi(M))+g(M)+o(N^{-2}),$		(25)

where $g(M)$ is a function independent of $\pi(M)$ . Substituting $\pi=\pi_{1}\pi_{2}$ and $\pi=\pi_{1}$ into (25) and taking difference, we obtain (24). ∎

By comparing Lemma 5.1 to Lemma A.2, we obtain the following connection between estimation and prediction.

Proposition 5.1.

For every $M$ ,

		$\displaystyle\lim_{N\to\infty}N^{2}({\rm E}_{M}[D(p(\widetilde{Y}\mid M),\hat{p}_{\pi_{1}\pi_{2}}(\widetilde{Y}\mid Y))]-{\rm E}_{M}[D(p(\widetilde{Y}\mid M),\hat{p}_{\pi_{1}}(\widetilde{Y}\mid Y))])$
	$\displaystyle=$	$\displaystyle\frac{1}{2}\lim_{N\to\infty}N^{2}({\rm E}_{M}[\\|\hat{M}^{\pi_{1}\pi_{2}}-M\\|_{\mathrm{F}}^{2}]-{\rm E}_{M}[\\|\hat{M}^{\pi_{1}}-M\\|_{\mathrm{F}}^{2}]).$

Therefore, if $\hat{\theta}^{\pi_{1}\pi_{2}}$ asymptotically dominates $\hat{\theta}^{\pi_{1}}$ under the quadratic loss, then $\hat{p}_{\pi_{1}\pi_{2}}(\widetilde{Y}\mid Y)$ asymptotically dominates $\hat{p}_{\pi_{1}}(\widetilde{Y}\mid Y)$ under the Kullback–Leibler loss.

Therefore, Theorems 3.1 and 4.1 are extended to Bayesian prediction as follows. Other results in the previous sections can be extended to Bayesian prediction similarly.

Theorem 5.1.

For every $M$ ,

		$\displaystyle N^{2}({\rm E}_{M}[D(p(\cdot\mid M),\hat{p}_{{\rm MSVS1}}(\cdot\mid Y))]-{\rm E}_{M}[D(p(\cdot\mid M),\hat{p}_{{\rm SVS}}(\cdot\mid Y))])$
	$\displaystyle\to$	$\displaystyle\frac{\gamma(\gamma-2p^{2}-2p+4)}{2{\rm tr}(M^{\top}M)}$

as $N\to\infty$ . Therefore, if $p\geq 2$ and $0<\gamma<p^{2}+p$ , then the Bayesian predictive density with respect to $\pi_{{\rm MSVS1}}$ in (7) asymptotically dominates that with respect to $\pi_{{\rm SVS}}$ in (2) under the Kullback–Leibler loss.

Theorem 5.2.

For every $M$ ,

		$\displaystyle N^{2}({\rm E}_{M}[D(p(\cdot\mid M),\hat{p}_{{\rm MSVS2}}(\cdot\mid Y))]-{\rm E}_{M}[D(p(\cdot\mid M),\hat{p}_{{\rm SVS}}(\cdot\mid Y))])$
	$\displaystyle\to$	$\displaystyle\frac{1}{2}\sum_{i=1}^{p}\gamma_{i}(\gamma_{i}-2p+2)\\|M_{\cdot i}\\|^{-2}$

as $N\to\infty$ . Therefore, if $p\geq 2$ and $0<\gamma\leq p$ , then the Bayesian predictive density with respect to $\pi_{{\rm MSVS2}}$ in (15) asymptotically dominates that with respect to $\pi_{{\rm SVS}}$ in (2) under the Kullback–Leibler loss.

Figures 7 and 8 plot the Kullback–Leibler risk of Bayesian predictive densities in similar settings to Figures 1 and 4, respectively. They show that the risk behavior in prediction is qualitatively the same with that in estimation, which is compatible with Theorems 5.1 and 5.2.

Figure 7: Kullback–Leibler risk of Bayesian predictive densities for

n=10

p=3

and

N=1

. Left:

\sigma_{2}=\sigma_{3}=0

. Right:

\sigma_{1}=10

\sigma_{3}=0

. solid:

\pi_{{\rm MSVS1}}

with

\gamma=p^{2}+p-2

, dashed:

\pi_{{\rm SVS}}

, dotted: Stein’s prior

\pi_{{\rm S}}

Figure 8: Kullback–Leibler risk of Bayesian predictive densities for

n=10

p=3

and

N=1

where

M=U\Sigma

with

U^{\top}U=I_{p}

and

\Sigma={\rm diag}(\sigma_{1},\dots,\sigma_{p})

. Left:

\sigma_{2}=\sigma_{3}=0

. Right:

\sigma_{1}=10

\sigma_{3}=0

. solid:

\pi_{{\rm MSVS2}}

with

\gamma_{1}=\dots=\gamma_{p}=p-1

, dash-dotted:

\pi_{{\rm MSVS1}}

with

\gamma=p^{2}+p-2

, dashed:

\pi_{{\rm SVS}}

, dotted: Stein’s prior

\pi_{{\rm S}}

Acknowledgments

We thank the referees for helpful comments. This work was supported by JSPS KAKENHI Grant Numbers 19K20220, 21H05205, 22K17865 and JST Moonshot Grant Number JPMJMS2024.

References

Aitchison (1975) Aitchison, J. (1975). Goodness of prediction fit. Biometrika 62, 547–554.
Amari (1985) Amari, S. (1985). Differential-Geometrical Methods in Statistics. New York: Springer.
Brown (1971) Brown, L. D. (1971). Admissible estimators, recurrent diffusions, and insoluble boundary value problems. Annals of Mathematical Statistics 36, 855–903.
Brown and Zhao (2009) Brown, L. D. & Zhao, L. H. (2009). Estimators for Gaussian models having a block-wise structure. Statistica Sinica 19, 885–903.
Efron and Morris (1972) Efron, B. & Morris, C. (1972). Empirical Bayes on vector observations: an extension of Stein’s method. Biometrika 59, 335–347.
Efron and Morris (1976) Efron, B. & Morris, C. (1976). Multivariate empirical Bayes and estimation of covariance matrices. Annals of Statistics 4, 22–32.
Fourdrinier et al. (2018) Fourdrinier, D., Strawderman, W. E. & Wells, M. (2018). Shrinkage Estimation. Springer, New York.
George, Liang and Xu (2006) George, E. I., Liang, F. & Xu, X. (2006). Improved minimax predictive densities under Kullback–Leibler loss. Annals of Statistics 34, 78–91.
Gupta and Nagar (2000) Gupta, A. K. & Nagar, D. K. (2000). Matrix Variate Distributions. New York: Chapman & Hall.
Komaki (2001) Komaki, F. (2001). A shrinkage predictive distribution for multivariate normal observables. Biometrika 88, 859–864.
Komaki (2006) Komaki, F. (2006). Shrinkage priors for Bayesian prediction. Annals of Statistics 34, 808–819.
Konno (1990) Konno, Y. (1990). Families of minimax estimators of matrix of normal means with unknown covariance matrix. J. Japan Statist. Soc. 20, 191–201.
Konno (1991) Konno, Y. (1991). On estimation of a matrix of normal means with unknown covariance matrix. Journal of Multivariate Analysis 36, 44–55.
Matsuda (2023) Matsuda, T. (2023). Adapting to arbitrary quadratic loss via singular value shrinkage. IEEE Transactions on Information Theory, accepted.
Matsuda (2024) Matsuda, T. (2024). Matrix quadratic risk of orthogonally invariant estimators for a normal mean matrix. Japanese Journal of Statistics and Data Science, accepted.
Matsuda and Komaki (2015) Matsuda, T. & Komaki, F. (2015). Singular value shrinkage priors for Bayesian prediction. Biometrika 102, 843–854.
Matsuda and Strawderman (2019) Matsuda, T. & Strawderman, W. E. (2019). Improved loss estimation for a normal mean matrix. Journal of Multivariate Analysis 169, 300–311.
Matsuda and Strawderman (2022) Matsuda, T. & Strawderman, W. E. (2022). Estimation under matrix quadratic loss and matrix superharmonicity. Biometrika 109, 503–519.
Stein (1974) Stein, C. (1974). Estimation of the mean of a multivariate normal distribution. Proc. Prague Symp. Asymptotic Statistics 2, 345–381.
Tsukuma and Kubokawa (2007) Tsukuma, H. & Kubokawa, T. (2007). Methods for improvement in estimation of a normal mean matrix. Journal of Multivariate Analysis 98, 1592–1610.
Tsukuma and Kubokawa (2017) Tsukuma, H. & Kubokawa, T. (2017). Proper Bayes and minimax predictive densities related to estimation of a normal mean matrix. Journal of Multivariate Analysis 159, 138–150.
Tsukuma and Kubokawa (2020) Tsukuma, H. & Kubokawa, T. (2020). Shrinkage estimation for mean and covariance matrices. Springer.
Yuasa and Kubokawa (2023a) Yuasa, R. & Kubokawa, T. (2023). Generalized Bayes estimators with closed forms for the normal mean and covariance matrices. Journal of Statistical Planning and Inference 222, 182–194.
Yuasa and Kubokawa (2023b) Yuasa, R. & Kubokawa, T. (2023). Weighted shrinkage estimators of normal mean matrices and dominance properties. Journal of Multivariate Analysis 194, 1–17.

Appendix A Asymptotic expansion of risk

Here, we provide asymptotic expansion formulas for estimators of a normal mean vector. Consider the problem of estimating $\theta$ from the observation $Y\sim{\rm N}_{d}(\theta,N^{-1}I_{d})$ under the quadratic loss $l(\theta,\hat{\theta})=\|\hat{\theta}-\theta\|^{2}$ . As shown in Stein (1974), the generalized Bayes estimator $\hat{\theta}^{\pi}$ with respect to a prior $\pi(\theta)$ is expressed as

\displaystyle\hat{\theta}^{\pi}(y)=y+\frac{1}{N}\nabla_{y}\log m_{\pi}(y),

where

\displaystyle m_{\pi}(y)=\int p(y\mid\theta)\pi(\theta){\rm d}\theta.

The asymptotic difference between the quadratic risk of two generalized Bayes estimators as $N\to\infty$ is given as follows.

Lemma A.1.

As $N\to\infty$ , the difference between the quadratic risk of $\hat{\theta}^{\pi_{1}}$ and $\hat{\theta}^{\pi_{1}\pi_{2}}$ is expanded as

		$\displaystyle{\rm E}_{\theta}[\\|\hat{\theta}^{\pi_{1}\pi_{2}}-\theta\\|^{2}]-{\rm E}_{\theta}[\\|\hat{\theta}^{\pi_{1}}-\theta\\|^{2}]$
	$\displaystyle=$	$\displaystyle\frac{1}{N^{2}}(2\nabla\log\pi_{1}(\theta)^{\top}\nabla\log\pi_{2}(\theta)+\\|\nabla\log\pi_{2}(\theta)\\|^{2}+2\Delta\log\pi_{2}(\theta))+o(N^{-2}).$		(26)

Proof.

By using Stein’s lemma (Fourdrinier et al., 2018) and $m_{\pi}(y)=\pi(y)+o(1)$ as $N\to\infty$ , the quadratic risk of the generalized Bayes estimator $\hat{\theta}^{\pi}$ is calculated as

	$\displaystyle{\rm E}_{\theta}[\\|\hat{\theta}^{\pi}(y)-\theta\\|^{2}]$
$\displaystyle=$	$\displaystyle{\rm E}_{\theta}[\\|y-\theta\\|^{2}]+\frac{2}{N}{\rm E}_{\theta}[(y-\theta)^{\top}\nabla\log m_{\pi}(y)]+\frac{1}{N^{2}}{\rm E}_{\theta}[\\|\nabla\log m_{\pi}(y)\\|^{2}]$
$\displaystyle=$	$\displaystyle\frac{d}{N}+\frac{1}{N^{2}}{\rm E}_{\theta}[\\|\nabla\log m_{\pi}(y)\\|^{2}+2\Delta\log m_{\pi}(y)]$
$\displaystyle=$	$\displaystyle\frac{d}{N}+\frac{1}{N^{2}}\left(\\|\nabla\log\pi(\theta)\\|^{2}+2\Delta\log\pi(\theta)\right)+o(N^{-2}).$	(27)

Substituting $\pi=\pi_{1}\pi_{2}$ and $\pi=\pi_{1}$ into (27) and taking difference, we obtain (26). ∎

We extend the above formula to matrices by using the matrix derivative notations from Matsuda and Strawderman (2022). For a function $f:\mathbb{R}^{n\times p}\to\mathbb{R}$ , its matrix gradient $\widetilde{\nabla}f:\mathbb{R}^{n\times p}\to\mathbb{R}^{n\times p}$ is defined as

\displaystyle(\widetilde{\nabla}f(X))_{ai}=\frac{\partial}{\partial X_{ai}}f(X).

(28)

For a $C^{2}$ function $f:\mathbb{R}^{n\times p}\to\mathbb{R}$ , its matrix Laplacian $\widetilde{\Delta}f:\mathbb{R}^{n\times p}\to\mathbb{R}^{p\times p}$ is defined as

\displaystyle(\widetilde{\Delta}f(X))_{ij}=\sum_{a=1}^{n}\frac{\partial^{2}}{\partial X_{ai}\partial X_{aj}}f(X).

(29)

Then, the above formulas can be straightforwardly extended to matrix-variate normal distributions as follows.

Lemma A.2.

As $N\to\infty$ , the difference between the Frobenius risk of $\hat{M}^{\pi_{1}}$ and $\hat{M}^{\pi_{1}\pi_{2}}$ is expanded as

		$\displaystyle{\rm E}_{M}[\\|\hat{M}^{\pi_{1}\pi_{2}}-M\\|_{\mathrm{F}}^{2}]-{\rm E}_{M}[\\|\hat{M}^{\pi_{1}}-M\\|_{\mathrm{F}}^{2}]$
	$\displaystyle=$	$\displaystyle\frac{1}{N^{2}}{\rm tr}(2\widetilde{\nabla}\log\pi_{1}(M)^{\top}\widetilde{\nabla}\log\pi_{2}(M)+\widetilde{\nabla}\log\pi_{2}(M)^{\top}\widetilde{\nabla}\log\pi_{2}(M)+2\widetilde{\Delta}\log\pi_{2}(M))$
		$\displaystyle\quad+o(N^{-2}).$

Lemma A.3.

As $N\to\infty$ , the difference between the matrix quadratic risk of $\hat{M}^{\pi_{1}}$ and $\hat{M}^{\pi_{1}\pi_{2}}$ is expanded as

		$\displaystyle{\rm E}_{M}[(\hat{M}^{\pi_{1}\pi_{2}}-M)^{\top}(\hat{M}^{\pi_{1}\pi_{2}}-M)-(\hat{M}^{\pi_{1}}-M)^{\top}(\hat{M}^{\pi_{1}}-M)]$
	$\displaystyle=$	$\displaystyle\frac{1}{N^{2}}(2\widetilde{\nabla}\log\pi_{1}(M)^{\top}\widetilde{\nabla}\log\pi_{2}(M)+\widetilde{\nabla}\log\pi_{2}(M)^{\top}\widetilde{\nabla}\log\pi_{2}(M)+2\widetilde{\Delta}\log\pi_{2}(M))$
		$\displaystyle\quad+o(N^{-2}).$

Komaki (2006) derived the asymptotic expansion of the Kullback–Leibler risk of Bayesian predictive densities. For the normal model as discussed in Section 5, the result shows that Stein’s prior dominates the Jeffreys prior in $O(N^{-1})$ term at the origin and $O(N^{-2})$ term at other points, which is reminiscent of superefficiency theory. A similar phenomenon should exist in estimation as well. Unlike Stein’s prior, the priors for a normal mean matrix such as $\pi_{\mathrm{SVS}}$ diverge at many points such as low-rank matrices. It is an interesting future problem to investigate the asymptotic risk of such priors in detail.

Appendix B Laplacian of $\pi_{{\rm MSVS1}}$ and $\pi_{{\rm MSVS2}}$

Lemma B.1.

(Stein, 1974; Matsuda and Strawderman, 2019) Suppose that $f:\mathbb{R}^{n\times p}\to\mathbb{R}$ is represented as $f(X)=\widetilde{f}(\sigma)$ , where $n\geq p$ and $\sigma=(\sigma_{1}(X),\ldots,\sigma_{p}(X))$ denotes the singular values of $X$ . If $f$ is twice weakly differentiable, then its Laplacian is

\displaystyle\Delta f

\displaystyle=\sum_{a=1}^{n}\sum_{i=1}^{p}\frac{\partial^{2}f}{\partial X_{ai}^{2}}=2\sum_{i<j}\frac{\sigma_{i}{\partial\widetilde{f}}/{\partial\sigma_{i}}-\sigma_{j}{\partial\widetilde{f}}/{\partial\sigma_{j}}}{\sigma_{i}^{2}-\sigma_{j}^{2}}+(n-p)\sum_{i=1}^{p}\frac{1}{\sigma_{i}}\frac{\partial\widetilde{f}}{\partial\sigma_{i}}+\sum_{i=1}^{p}\frac{\partial^{2}\widetilde{f}}{\partial\sigma_{i}^{2}}.

Proposition B.1.

The Laplacian of $\pi_{{\rm MSVS1}}$ in (7) is given by

\displaystyle\Delta\pi_{{\rm MSVS1}}(M)

\displaystyle=\gamma(\gamma+np-2p^{2}-2p+2)\|M\|_{{\rm F}}^{-2}\pi_{{\rm MSVS1}}(M).

Proof.

Let

\displaystyle\widetilde{f}(\sigma)=\left(\prod_{i=1}^{p}\sigma_{i}^{-(n-p-1)}\right)\left(\sum_{i=1}^{p}\sigma_{i}^{2}\right)^{-\gamma/2}

so that $\pi_{{\rm MSVS1}}(M)=\widetilde{f}(\sigma)$ with $\sigma=(\sigma_{1}(M),\dots,\sigma_{p}(M))$ . From Lemma B.1 and

\displaystyle\frac{\partial\widetilde{f}}{\partial\sigma_{i}}=-(n-p-1)\sigma_{i}^{-1}\widetilde{f}-\gamma\sigma_{i}\left(\sum_{j=1}^{p}\sigma_{j}^{2}\right)^{-1}\widetilde{f},

\displaystyle\frac{\partial^{2}\widetilde{f}}{\partial\sigma_{i}^{2}}=(n-p)(n-p-1)\sigma_{i}^{-2}\widetilde{f}+(2n-2p-3)\gamma\left(\sum_{j=1}^{p}\sigma_{j}^{2}\right)^{-1}\widetilde{f}+\gamma(\gamma+2)\sigma_{i}^{2}\left(\sum_{j=1}^{p}\sigma_{j}^{2}\right)^{-2}\widetilde{f},

we have

	$\displaystyle\Delta\pi_{{\rm MSVS1}}(M)$	$\displaystyle=2\sum_{i<j}\frac{\sigma_{i}{\partial\widetilde{f}}/{\partial\sigma_{i}}-\sigma_{j}{\partial\widetilde{f}}/{\partial\sigma_{j}}}{\sigma_{i}^{2}-\sigma_{j}^{2}}+(n-p)\sum_{i=1}^{p}\frac{1}{\sigma_{i}}\frac{\partial\widetilde{f}}{\partial\sigma_{i}}+\sum_{i=1}^{p}\frac{\partial^{2}\widetilde{f}}{\partial\sigma_{i}^{2}}$
		$\displaystyle=-2\cdot\frac{p(p-1)}{2}\gamma\left(\sum_{i=1}^{p}\sigma_{i}^{2}\right)^{-1}\widetilde{f}-(n-p)(n-p-1)\left(\sum_{i=1}^{p}\sigma_{i}^{-2}\right)\widetilde{f}$
		$\displaystyle\quad-p(n-p)\gamma\left(\sum_{i=1}^{p}\sigma_{i}^{2}\right)^{-1}\widetilde{f}+(n-p)(n-p-1)\left(\sum_{i=1}^{p}\sigma_{i}^{-2}\right)\widetilde{f}$
		$\displaystyle\quad+\gamma(p(2n-2p-3)+\gamma+2)\left(\sum_{i=1}^{p}\sigma_{i}^{2}\right)^{-1}\widetilde{f}$
		$\displaystyle=\gamma(\gamma+np-2p^{2}-2p+2)\left(\sum_{i=1}^{p}\sigma_{i}^{2}\right)^{-1}\widetilde{f}.$

∎

Proposition B.2.

The Laplacian of $\pi_{{\rm MSVS2}}$ in (15) is given by

\displaystyle\Delta\pi_{{\rm MSVS2}}(M)

\displaystyle=\gamma(\gamma+n-2p)\left(\sum_{i=1}^{p}\|M_{\cdot i}\|^{-2}\right)\pi_{{\rm MSVS2}}(M).

Proof.

From

\displaystyle\Delta\log f=\frac{\Delta f}{f}-\|\nabla\log f\|^{2}

and (20), (21) and (22),

	$\displaystyle\frac{\Delta\pi_{{\rm MSVS2}}(M)}{\pi_{{\rm MSVS2}}(M)}$	$\displaystyle=\Delta\log\pi_{{\rm MSVS2}}(M)+\\|\nabla\log\pi_{{\rm MSVS2}}(M)\\|^{2}$
		$\displaystyle=\Delta\log\pi_{{\rm SVS}}(M)+\Delta\log\pi_{{\rm CS}}(M)+\\|\nabla\log\pi_{{\rm SVS}}(M)+\nabla\log\pi_{{\rm CS}}(M)\\|^{2}$
		$\displaystyle=\frac{\Delta\pi_{{\rm SVS}}(M)}{\pi_{{\rm SVS}}(M)}+(-(n-2)\gamma+2(n-p-1)\gamma+\gamma^{2})\sum_{i=1}^{p}\\|M_{\cdot i}\\|^{-2}$
		$\displaystyle=\gamma(\gamma+n-2p)\sum_{i=1}^{p}\\|M_{\cdot i}\\|^{-2},$

where we used $\Delta\pi_{{\rm SVS}}(M)=0$ (Theorem 2 of Matsuda and Komaki, 2015). ∎

Appendix C Improving on the block-wise Stein prior

Here, we develop priors that asymptotically dominate the block-wise Stein prior in estimation and prediction. Suppose that we observe $Y\sim{\rm N}_{d}(\theta,N^{-1}I_{d})$ and estimate $\theta$ or predict $\widetilde{Y}\sim{\rm N}_{d}(\theta,I_{d})$ . We assume that the $d$ -dimensional mean vector $\theta$ is split into $B$ disjoint blocks $\theta^{(1)},\dots,\theta^{(B)}$ with size $d_{1},\dots,d_{B}$ , where $d_{1}+\dots+d_{B}=d$ . For example. such a situation appears in balanced ANOVA and wavelet regression (Brown and Zhao, 2009). Then, the block-wise Stein prior is defined as

\displaystyle\pi_{{\rm BS}}(\theta)=\prod_{b=1}^{B}\|\theta^{(b)}\|^{R_{b}},\quad R_{b}=-(d_{b}-2)_{+},

(30)

which puts Stein’s prior on each block. Since it is superharmonic, the generalized Bayes estimator $\hat{\theta}^{\pi_{{\rm BS}}}$ with respect to $\pi_{{\rm BS}}$ is minimax. However, Brown and Zhao (2009) showed that $\hat{\theta}^{\pi_{{\rm BS}}}$ is inadmissible and dominated by an estimator with additional James–Stein type shrinkage defined by

\displaystyle\hat{\theta}(y)=\hat{\theta}^{\pi_{{\rm BS}}}(y)-\frac{R_{\#}+d-2}{\|y\|^{2}}y,

where $R_{\#}=\sum_{b}R_{b}>2-d$ . From this result, Brown and Zhao (2009) conjectured that the block-wise Stein prior can be improved by multiplying a Stein-type shrinkage prior in Remark 3.2. Following their conjecture, we construct priors by adding scalar shrinkage to the block-wise Stein priors:

\displaystyle\pi_{{\rm MBS}}(\theta)=\pi_{{\rm BS}}(\theta)\|\theta\|^{-\gamma},

(31)

where $\gamma\geq 0$ . Let

\displaystyle m_{{\rm MBS}}(y)=\int p(y\mid\theta)\pi_{{\rm MBS}}(\theta){\rm d}\theta.

Lemma C.1.

If $0\leq\gamma<B(R_{b}+d_{b})$ for every $b$ , then $m_{{\rm MBS}}(y)<\infty$ for every $y$ .

Proof.

Since $m_{{\rm MBS}}(y)$ is interpreted as the expectation of $\pi_{{\rm MBS}}(\theta)$ under $\theta\sim{\rm N}_{d}(y,I_{d})$ , it suffices to show that $\pi_{{\rm MBS}}(\theta)$ is locally integrable at every $\theta$ .

First, consider $\theta\neq 0$ . Since $m_{{\rm BS}}(y)<\infty$ for every $y$ (Brown and Zhao, 2009), $\pi_{{\rm BS}}(\theta)$ is locally integrable at $\theta$ . Also, $\|\theta\|>c$ for some $c>0$ in a neighborhood of $\theta$ . Thus, $\pi_{{\rm MBS}}(\theta)=\pi_{{\rm BS}}(\theta)\|\theta\|^{-\gamma}$ is locally integrable at $\theta$ .

Next, consider $\theta=0$ and take the neighborhood $A=\{\theta\mid\|\theta^{(1)}\|\leq s,\dots,\|\theta^{(B)}\|\leq s\}$ for $s>0$ . From the AM-GM inequality,

\displaystyle\|\theta\|^{2}=\sum_{b}\|\theta^{(b)}\|^{2}\geq B\left(\prod_{b}\|\theta^{(b)}\|^{2}\right)^{1/B}=B\prod_{b}\|\theta^{(b)}\|^{2/B}.

Thus,

\displaystyle\int_{A}\pi_{{\rm MBS}}(\theta){\rm d}\theta\leq\

\displaystyle C\int_{0}^{s}\cdots\int_{0}^{s}\prod_{b}r_{b}^{R_{b}+d_{b}-1-\gamma/B}{\rm d}r_{1}\cdots{\rm d}r_{B},

where $r_{b}=\|\theta^{(b)}\|$ and $C$ is a constant. Therefore, $\pi_{{\rm MBS}}(\theta)$ is locally integrable at $\theta=0$ if $R_{b}+d_{b}-1-\gamma/B>-1$ for every $b$ , which is equivalent to $\gamma<B(R_{b}+d_{b})$ for every $b$ . ∎

From Lemma C.1, the generalized Bayes estimator with respect to $\pi_{{\rm MBS}}$ is well-defined when $0\leq\gamma<B(R_{b}+d_{b})$ for every $b$ . We denote it by $\hat{\theta}_{{\rm MBS}}$ .

Theorem C.1.

For every $M$ ,

\displaystyle N^{2}({\rm E}_{\theta}[\|\hat{\theta}_{{{\rm MBS}}}-\theta\|^{2}]-{\rm E}_{\theta}[\|\hat{\theta}_{{{\rm BS}}}-\theta\|^{2}])\to\gamma(\gamma-2(R_{\#}+d-2))\|\theta\|^{-2}

(32)

as $N\to\infty$ . Therefore, if $0<\gamma<2(R_{\#}+d-2)$ , then the generalized Bayes estimator with respect to $\pi_{{\rm MBS}}$ in (31) asymptotically dominates that with respect to $\pi_{{\rm BS}}$ in (30) under the Frobenius loss.

Proof.

Let $\pi_{{\rm S}}(\theta)=\|\theta\|^{-\gamma}$ . By straightforward calculation, we obtain

	$\displaystyle\nabla\log\pi_{{\rm BS}}(\theta)^{\top}\nabla\log\pi_{{\rm S}}(\theta)=-\gamma R_{\#}\\|\theta\\|^{-2},$
	$\displaystyle\nabla\log\pi_{{\rm S}}(\theta)^{\top}\nabla\log\pi_{{\rm S}}(\theta)=\gamma^{2}\\|\theta\\|^{-2},$
	$\displaystyle\Delta\log\pi_{{\rm S}}(\theta)=-\gamma(d-2)\\|\theta\\|^{-2}.$

Therefore, from Lemma A.1,

		$\displaystyle{\rm E}_{\theta}[\\|\hat{\theta}_{{{\rm MBS}}}-\theta\\|^{2}]-{\rm E}_{\theta}[\\|\hat{\theta}_{{{\rm BS}}}-\theta\\|^{2}]$
	$\displaystyle=$	$\displaystyle\frac{1}{N^{2}}\left(2\nabla\log\pi_{{\rm BS}}(\theta)^{\top}\nabla\log\pi_{{\rm S}}(\theta)+\\|\nabla\log\pi_{{\rm S}}(\theta)\\|^{2}+2\Delta\log\pi_{{\rm S}}(\theta)\right)+o(N^{-2})$
	$\displaystyle=$	$\displaystyle\frac{1}{N^{2}}\gamma(\gamma-2(R_{\#}+d-2))\\|\theta\\|^{-2}+o(N^{-2}).$

Hence, we obtain (32). ∎

From (32), the choice $\gamma=R_{\#}+d-2$ is optimal. As discussed in Section 5, Theorem C.1 is extended to Bayesian prediction as follows.

Theorem C.2.

For every $M$ ,

\displaystyle N^{2}({\rm E}_{\theta}[D(p(\cdot\mid\theta),\hat{p}_{{\rm MBS}}(\cdot\mid y))]-{\rm E}_{\theta}[D(p(\cdot\mid\theta),\hat{p}_{{\rm BS}}(\cdot\mid y))])\to\frac{\gamma(\gamma-2(R_{\#}+d-2))}{2}\|\theta\|^{-2}

as $N\to\infty$ . Therefore, if $p\geq 2$ and $0<\gamma<p^{2}+p$ , then the Bayesian predictive density with respect to $\pi_{{\rm MBS}}$ in (31) asymptotically dominates that with respect to $\pi_{{\rm BS}}$ in (30) under the Kullback–Leibler loss.

	$\displaystyle\\|D^{1/2}M^{\top}MD^{1/2}\\|$	$\displaystyle\leq\\|D^{1/2}\\|\\|M^{\top}M\\|\\|D^{1/2}\\|$
		$\displaystyle=(\max_{i}D_{ii})\\|M\\|^{2}$
		$\displaystyle=\gamma_{1}\frac{\\|M\\|^{2}}{\min_{i}\\|M_{\cdot i}\\|^{2}}$
		$\displaystyle\leq\gamma_{1}\frac{\sum_{i=1}^{p}\\|M_{\cdot i}\\|^{2}}{\min_{i}\\|M_{\cdot i}\\|^{2}}$
		$\displaystyle\leq\gamma_{1}p$
		$\displaystyle<2(p-1),$

	$\displaystyle\frac{\Delta\pi_{{\rm MSVS2}}(M)}{\pi_{{\rm MSVS2}}(M)}$	$\displaystyle=\Delta\log\pi_{{\rm MSVS2}}(M)+\\|\nabla\log\pi_{{\rm MSVS2}}(M)\\|^{2}$
		$\displaystyle=\Delta\log\pi_{{\rm SVS}}(M)+\Delta\log\pi_{{\rm CS}}(M)+\\|\nabla\log\pi_{{\rm SVS}}(M)+\nabla\log\pi_{{\rm CS}}(M)\\|^{2}$
		$\displaystyle=\frac{\Delta\pi_{{\rm SVS}}(M)}{\pi_{{\rm SVS}}(M)}+(-(n-2)\gamma+2(n-p-1)\gamma+\gamma^{2})\sum_{i=1}^{p}\\|M_{\cdot i}\\|^{-2}$
		$\displaystyle=\gamma(\gamma+n-2p)\sum_{i=1}^{p}\\|M_{\cdot i}\\|^{-2},$

		$\displaystyle{\rm E}_{\theta}[\\|\hat{\theta}_{{{\rm MBS}}}-\theta\\|^{2}]-{\rm E}_{\theta}[\\|\hat{\theta}_{{{\rm BS}}}-\theta\\|^{2}]$
	$\displaystyle=$	$\displaystyle\frac{1}{N^{2}}\left(2\nabla\log\pi_{{\rm BS}}(\theta)^{\top}\nabla\log\pi_{{\rm S}}(\theta)+\\|\nabla\log\pi_{{\rm S}}(\theta)\\|^{2}+2\Delta\log\pi_{{\rm S}}(\theta)\right)+o(N^{-2})$
	$\displaystyle=$	$\displaystyle\frac{1}{N^{2}}\gamma(\gamma-2(R_{\#}+d-2))\\|\theta\\|^{-2}+o(N^{-2}).$

Double shrinkage priors for a normal mean matrix

Abstract

1 Introduction

2 Inadmissibility of the singular value shrinkage prior

Lemma 2.1.

Theorem 2.1.

Proof.

3 Improvement by additional scalar shrinkage

Lemma 3.1.

Proof.

Theorem 3.1.

Proof.

Theorem 3.2.

Proof.

Theorem 3.3.

Proof.

Lemma 3.2.

4 Improvement by additional column-wise shrinkage

Lemma 4.1.

Proof.

Theorem 4.1.

Proof.

Theorem 4.2.

Proof.

Theorem 4.3.

Proof.

5 Bayesian prediction

Lemma 5.1.

Proof.

Proposition 5.1.

Theorem 5.1.

Theorem 5.2.

Acknowledgments

References

Appendix A Asymptotic expansion of risk

Lemma A.1.

Proof.

Lemma A.2.

Lemma A.3.

Appendix B Laplacian of πMSVS1\pi_{{\rm MSVS1}} and πMSVS2\pi_{{\rm MSVS2}}

Lemma B.1.

Proposition B.1.

Proof.

Proposition B.2.

Proof.

Appendix C Improving on the block-wise Stein prior

Lemma C.1.

Proof.

Theorem C.1.

Proof.

Theorem C.2.

Appendix B Laplacian of $\pi_{{\rm MSVS1}}$ and $\pi_{{\rm MSVS2}}$