Bayesian sparse multiple regression for simultaneous rank reduction and variable selection
Antik Chakraborty
Anirban Bhattacharya
Bani K. Mallick
Antik Chakraborty
Anirban Bhattacharya
Bani K. Mallick
%\baselineskip=16pt \vskip 10mm \begin{center {\Large{Abstract \end{center \baselineskip=14pt We develop a Bayesian methodology aimed at simultaneously estimating low-rank and row-sparse matrices in a high-dimensional multiple-response linear regression model. We consider a carefully devised shrinkage prior on the matrix of regression coefficients which obviates the need to specify a prior on the rank, and shrinks the regression matrix towards low-rank and row-sparse structures. We provide theoretical support to the proposed methodology by proving minimax optimality of the posterior mean under the prediction risk in ultra-high dimensional settings where the number of predictors can grow sub-exponentially relative to the sample size. A one-step post-processing scheme induced by group lasso penalties on the rows of the estimated coefficient matrix is proposed for variable selection, with default choices of tuning parameters. We additionally provide an estimate of the rank using a novel optimization function achieving dimension reduction in the covariate space. We exhibit the performance of the proposed methodology in an extensive simulation study and a real data example. %
sparse_reduced_rank.tex
\documentclass[11pt, oneside]{article}   
\usepackage[utf8]{inputenc}
\usepackage{amsmath,amssymb,amsthm,algorithmicx,algorithm,setspace,url}
\usepackage{graphicx,bbm}
\usepackage[top=0.6in,margin=1in]{geometry}
\usepackage[authoryear, sort]{natbib}
%\usepackage{natbib}
\usepackage{authblk,booktabs,epstopdf,color,lscape,float}
\usepackage[hidelinks]{hyperref}
\usepackage{multirow}
\usepackage[space]{grffile}
%opening
\date{}
\def\Gauss{{\mathrm{N}}}
\def\Dmat{\mathrm{D}}
\def\prime{\mathrm{\scriptscriptstyle{T}}}
\def\eyen{\mathrm{I_n}}
\def\mr{\mathrm}
\def\ind{\mathbbm{1}}
\def\LS{{\textbf{LASSO}}}
\def\HS{{\textbf{HS}}}
\def\T{\mathrm{\scriptscriptstyle{T}}}
\def\tcr{\color{red}}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{assumption}[theorem]{Assumption}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{rem}{Remark}[section]
\newcommand{\norm}[1]{\left\lVert #1 \right\rVert}
\DeclareMathOperator*{\argmin}{arg\,min}
\begin{document}
\thispagestyle{empty}
\baselineskip=28pt
\begin{center}
{\LARGE{\bf  Bayesian sparse multiple regression for simultaneous rank reduction and variable selection}}

\end{center}
%\title{Bayesian sparse multiple regression for simultaneous rank reduction and variable selection}
\baselineskip=12pt

\vskip 10mm
\begin{center}
Antik Chakraborty\\
Department of Statistics, Texas A\&M University, College Station\\
3143 TAMU, TX 77843-3143, USA\\
antik@stat.tamu.edu\\
\hskip 5mm\\
Anirban Bhattacharya\\
Department of Statistics, Texas A\&M University, College Station\\
3143 TAMU, TX 77843-3143, USA\\
anirbanb@stat.tamu.edu\\
\hskip 5mm\\
Bani K. Mallick\\
Department of Statistics, Texas A\&M University, College Station\\
3143 TAMU, TX 77843-3143, USA\\
bmallick@stat.tamu.edu\\

\end{center}

%\author{Antik Chakraborty, Anirban Bhattacharya, Bani K. Mallick}
%\affil{Department of Statistics\\Texas A\&M University\\College Station, Texas}
\baselineskip=22pt

%\maketitle
%\begin{abstract}
%\baselineskip=16pt
\vskip 10mm
\begin{center}
{\Large{\bf Abstract}}
\end{center}
\baselineskip=14pt
We develop a Bayesian methodology aimed at simultaneously estimating low-rank and row-sparse matrices in a high-dimensional multiple-response linear regression model. We consider a carefully devised shrinkage prior on the matrix of regression coefficients which obviates the need to specify a prior on the rank, and shrinks the regression matrix towards low-rank and row-sparse structures. We provide theoretical support to the proposed methodology by proving minimax optimality of the posterior mean under the prediction risk in ultra-high dimensional settings where the number of predictors can grow sub-exponentially relative to the sample size. A one-step post-processing scheme induced by group lasso penalties on the rows of the estimated coefficient matrix is proposed for variable selection, with default choices of tuning parameters. We additionally provide an estimate of the rank using a novel optimization function achieving dimension reduction in the covariate space. We exhibit the performance of the proposed methodology in an extensive simulation study and a real data example.
%\end{abstract}
\vspace{1.5in}


\baselineskip=12pt
\par\vfill\noindent
\underline{\bf Key Words}: Bayesian; High dimension; Shrinkage prior; Posterior concentration; Dimension reduction; Variable selection.  

\par\medskip\noindent
\underline{\bf Short title}: Bayesian sparse multi-task learner

\clearpage\pagebreak\newpage
\pagenumbering{arabic}
\newlength{\gnat}
\setlength{\gnat}{22pt}
\baselineskip=\gnat

%\pagebreak

\section{Introduction}
Studying the relationship between multiple response variables and a set of predictors has broad applications ranging from bioinformatics, econometrics, time series analysis to growth curve models. The least squares solution in a linear multiple response regression problem is equivalent to performing separate least squares on each of the responses \citep{anderson1984multivariate} and ignores any potential dependence among the responses. 
%Moreover, it lacks interpretability due to the large number of parameters involved in the model. 
%Dimension reduction techniques such as principal component analysis, canonical correlation, factor models have gained wide popularity in multivariate exploratory data analysis due mainly to their ease of interpretation. 
In the context of multiple response regression, a popular technique to achieve parsimony and interpretability is to consider a reduced-rank decomposition of the coefficient matrix, commonly known as reduced rank regression \citep{anderson1951estimating,izenman1975reduced,velu2013multivariate}. 
%, wherein a low rank structure of the coefficient matrix is considered. 
%The study of such procedures date back to \cite{anderson1951estimating}, \cite{izenman1975reduced}; a thorough review of the developments in the last century is documented in the monograph by \cite{velu2013multivariate}. 
Although many results exist about the asymptotic properties of reduced rank estimators \citep{anderson2002specification}, formal statistical determination of the rank remains difficult even with fixed number of covariates and large sample size. 
%The problem is substantially harder when the number of covariates is large and thus limits the application of classical reduced rank procedures to routine modern statistical applications such as genomics, financial time series etc.\citep{bunea2011optimal}. 
%This motivated a series of developments aimed at simultaneously estimating the rank and handle large number of covariates. One unifying theme of this line of work relies on penalizing matrices with high ranks which is achieved either by penalizing the singular values of the coefficient matrix \citep{yuan2007dimension,chen2013reduced} or by penalizing the rank itself \citep{bunea2011optimal}. 
The problem becomes substantially harder when a large number of covariates are present, and has motivated a series of recent work on penalized estimation of low rank matrices, where either the singular values of the coefficient matrix \citep{yuan2007dimension,chen2013reduced}, or the rank itself \citep{bunea2011optimal} is penalized. 
Theoretical evaluations of these estimators focusing on adaptation to the oracle convergence rate when the true coefficient matrix is of low rank has been conducted \citep{bunea2011optimal}. It has also been noted \citep{bunea2012joint} that the convergence rate can be improved when the true coefficient matrix has zero rows and variable selection is incorporated within the estimation procedure. Methods that simultaneously handle rank reduction and variable selection include \cite{yuan2007dimension,bunea2012joint,chen2012sparse}. To best of our knowledge, uncertainty characterization for the parameter estimates from these procedures is currently not available. 
%Several authors \citep{tibshirani1996regression,yuan2006model} have additionally noted that the convergence rate can be improved further if one accounts for selecting variables  when the coefficient matrix has many zero rows and low rank. Methods that simultaneously handle low rank and row sparse matrices were then proposed \citep{yuan2007dimension,bunea2012joint,chen2012sparse} penalizing matrices with high rank and dense rows. 
%Uncertainty characterization in the estimation of parameters obtained from these procedures is currently not available to the best of our knowledge.

The first fully systematic Bayesian treatment of reduced rank regression was carried out in \cite{geweke1996bayesian}, where  conditioned on the rank, independent Gaussian priors were placed on the elements of the coefficient matrix. While formal Bayesian model selection can be performed to determine the rank \citep{geweke1996bayesian}, calculation of marginal likelihoods for various candidate ranks gets computationally burdensome with increasing dimensions.
%\cite{geweke1996bayesian} advocated using Bayes factors for determining the rank which get computationally burdensome with increasiang dimensions. 
%Simplified analytic approximations of the Bayes factor, such as the Savage-Dickey ratio \citep{dickey1971weighted}, exist in the literature; \cite{marin2010resolving} pointed out a measure theoretic issue with such procedures and proposed a modification. or,
%Simplified analytic approximations to the Bayes factor, such as the Savage-Dickey ratio \citep{dickey1971weighted}, suffer from well-kn
The problem of choosing the rank is not unique to reduced rank regression and is ubiquitous in situations involving low rank decompositions, with factor models being a prominent example. \cite{lopes2004bayesian} placed a prior on the number of factors and proposed a computationally intensive reversible jump algorithm \citep{green1995reversible} for model fitting. 
%see also \cite{godsill2001relationship} for related methods. However, reversible jump MCMC is usually expensive in terms of computation time and mixing.  
As an alternative, \cite{bhattacharya2011sparse} proposed to increasingly shrink the factors starting with a conservative upper bound and adaptively collapsing redundant columns inside their MCMC algorithm.
Recent advancements in Bayesian matrix factorization have taken a similar approach; see for example, \cite{lim2007variational,salakhutdinov2008bayesian,babacan2011variational,alquier2013bayesian}. 
%However, it is not clear how choices of hyper parameters involved in these priors may influence the estimation of the rank.

From a Bayesian point of view, a natural way to select variables in a single-response regression framework is to use point mass mixture priors \citep{george1993variable,scott2010bayes} which allow a subset of the regression coefficients to be exactly zero. These priors were also adapted to multiple response regression by several authors \citep{brown1998multivariate,lucas2006sparse,wang2010sparse,bhadra2013joint}.
%where independent point mass mixture priors were placed on the rows of the coefficient matrix. Various Gibbs sampling strategies to implement these priors are discussed in \cite{dellaportas2000bayesian}. 
Posterior inference with such priors involves a stochastic search over an exponentially growing model space and is computationally expensive even in moderate dimensions. To alleviate the computational burden, a number of continuous shrinkage priors have been proposed in the literature which mimic the operating characteristics of the discrete mixture priors. Such priors can be expressed as Gaussian scale mixtures \citep{polson2010shrink}, leading to block updates of model parameters; see \cite{bhattacharya2016fast} for a review of such priors and efficient implementations in high-dimensional settings. To perform variable selection with these continuous priors, several methods for post-processing the posterior distribution have been proposed \citep{bondell2012consistent,kundu2013bayes,hahn2015decoupling}. 



%On the other hand, by applying continuous shrinkage priors \citep{polson2010shrink} one can mimic the operating characteristics of point mass mixtures while avoiding the computational bottleneck. The scale mixture of Gaussian representation of these priors renders straightforward Gibbs update of the regression coefficients in a block, with priors in \cite{bhattacharya2015dirichlet} and \cite{carvalho2010horseshoe} having proven rapid mixing and convergence \citep{khare2013geometric,pal2014geometric}. However, such block updates from a high dimensional Gaussian distribution can be very expensive \citep{bhattacharya2015fast}.
%Also, these priors do not achieve variable selection. \cite{hahn2015decoupling} proposed posterior summary aimed at selecting variables with high predictive power to remedy this deficiency. Such posterior summaries have also been proposed in \cite{bondell2012consistent} and \cite{kundu2013bayes}.

In this article we simultaneously address the problems of dimension reduction and variable selection in high-dimensional reduced rank models from a Bayesian perspective. We develop a novel shrinkage prior on the coefficient matrix which encourages shrinkage towards low-rank and row-sparse matrices. The shrinkage prior is induced from appropriate shrinkage priors on the components of a full-rank decomposition of the coefficient matrix, and hence bypasses the need to specify a prior on the rank. We provide theoretical understanding into the operating characteristics of the proposed prior in terms of a novel prior concentration result around rank-reduced and low-sparse matrices. The prior concentration result is utilized to prove minimax concentration rates of the posterior under the fractional posterior framework of \cite{bhattacharya2016bayesian} in a ultrahigh-dimensional setting where the number of predictor variables can grow sub-exponentially in the sample size. 

The continuous nature of the prior enables efficient block updates of parameters inside a Gibbs sampler. In particular, we adapt an algorithm for sampling structured multivariate Gaussians from \cite{bhattacharya2016fast} to efficiently sample a high-dimensional matrix in a block. We propose two independent post-processing schemes to achieve row sparsity and rank reduction with encouraging performance. A key feature of our post-processing schemes is to exploit the posterior summaries to offer careful default choices of tuning parameters, resulting in a procedure which is completely free of tuning parameters. The resulting row-sparse and rank-reduced coefficient estimate is called a Bayesian sparse multi-task learner (BSML). We illustrate the superiority of BSML over its competitors through a detailed simulation study and the methodology is applied to a Yeast cell cycle data set. 
%To the best of our knowledge, no existing Bayesian methodology handles these two aspects at the same time. We develop our methodology with three main building blocks. First, instead of specifying a prior distribution on the rank, we shrink the entire coefficient matrix towards the space of low-rank and row-sparse matrices using continuous shrinkage priors, resulting in sufficient computational gains relative to priors that aim to select variables and/or the rank. We develop an exact sampling scheme to update the entire coefficient matrix in a block leading to cheaper computational cost and better mixing. Posterior summaries of the coefficient matrix, for example the pointwise mean or median, are then processed via two independent screening procedures: one which achieves row-sparsity, another which reduces dimensionality via row sparsity. One of the key features of our post-processing schemes is to come up with careful choices of tuning parameters based on the posterior, resulting in a procedure which is completely free of tuning parameters. We call the resulting procedure Bayesian sparse multi-task learner (BSML).

\section{Bayesian sparse multitask learner}
\subsection{Model and Prior Specification}\label{sec:mod}
Suppose, for each observational unit $i = 1, \ldots, n$, we have a multivariate response $y_i \in \Re^q$ on $q$ variables of interest, along with information on $p$ possible predictors $x_i \in \Re^p$, a subset of which are assumed to be important in predicting the $q$ responses. Let $X \in \Re^{n \times p}$ denote the design matrix whose $i$th row is $x_i^{\T}$, and $Y \in \Re^{n\times q}$ the matrix of responses with the $i$th row as $y_i^{\T}$. The multivariate linear regression model is,
\begin{align}\label{eq:lin_reg}
Y=XC+E, \hspace{0.2in}E=(e_1^\T,\ldots,e_n^\T)^\T,
\end{align}
where we follow standard practice to center the response and exclude the intercept term. The rows of the error matrix are independent, with $e_i \sim \Gauss(0,\Sigma)$. Our main motivation is the high-dimensional case where $p \ge \max\{n, q\}$, although the method trivially applies to $p < n$ settings as well. We shall also assume the dimension of the response $q$ to be modest relative to the sample size. 

The basic assumption in reduced rank regression is that $\text{rank}(C)=r\leq \mathrm{min}(p,q)$, whence $C$ admits a decomposition $C = B_*A_*^\T$ with $B_* \in \Re^{p \times r}$ and $A_* \in \Re^{q \times r}$. While it is possible to treat $r$ as a parameter and assign it a prior distribution inside a hierarchical formulation, posterior inference on $r$ requires calculation of intractable marginal likelihoods or resorting to complicated reversible jump Markov chain Monte Carlo algorithms. To avoid specifying a prior on $r$, we work within a parameter-expanded framework \citep{liu1999parameter} to consider a potentially full-rank decomposition $C = B A^{\T}$ with $B\in \Re^{p \times q}$ and $A \in \Re^{q \times q}$, and assign shrinkage priors to $A$ and $B$ to shrink out the redundant columns when $C$ is indeed low rank. This formulation embeds all reduced-rank models inside the full model; if a conservative upper bound $q^\ast \le q$ on the rank is known, the method can be modified accordingly. The role of the priors on $B$ and $A$ is important to encourage appropriate shrinkage towards reduced-rank models, which is discussed below. 

We consider independent standard normal priors on the entries of $A$. As an alternative, a uniform prior on the Stiefel manifold \citep{hoff2007gibbs} of orthogonal matrices can be used. However, our numerical results suggested significant gains in computation time using the Gaussian prior over the uniform prior with no discernible difference in statistical performance. The Gaussian prior allows an efficient block update of $\mbox{vec}(A)$, whereas the algorithm of \cite{hoff2007gibbs} involves conditional Gibbs update of each column of $A$. Our theoretical results also suggest that the shrinkage provided by the Gaussian prior is optimal when $q$ is modest relative to $n$, the regime we operate in. We shall henceforth denote $\Pi_A$ to denote the prior on $A$, i.e., $a_{hk} \sim \mbox{N}(0, 1)$ independently for $h, k = 1, \ldots, q$. 

Recalling that the matrix $B$ has dimension $p \times q$, with $p$ potentially larger than $n$, stronger shrinkage is warranted on the columns of $B$. We use independent horseshoe priors \citep{carvalho2010horseshoe} on the columns of $B$, which can be represented hierarchically as 
\begin{align}\label{eq:HS}
b_{jh} \mid \lambda_{jh},\tau_h \sim \Gauss(0,\lambda_{jh}^2 \tau_h^2), \quad \lambda_{jh} \sim \mbox{Ca}_+(0, 1), \quad \tau_h \sim \mbox{Ca}_+(0, 1),
\end{align}
independently for $j = 1, \ldots, p$ and $h = 1, \ldots, q$, where $\mbox{Ca}_+(0, 1)$ denotes the truncated standard half-Cauchy distribution with density proportional to $(1+t^2)^{-1} \mathbbm{1}_{(0, \infty)}(t)$. We shall denote the prior on the matrix $B$ induced by the hierarchy in \eqref{eq:HS} by $\Pi_B$. 


We shall primarily restrict attention to settings where $\Sigma$ is diagonal, $\Sigma = \mbox{diag}(\sigma_1^2, \ldots, \sigma_q^2)$, noting that extensions to non-diagonal $\Sigma$ can be incorporated in a straightforward fashion. For example, for moderate $q$, a conjugate inverse-Wishart prior can be used as a default. Furthermore, if $\Sigma$ has a factor model or Gaussian Markov random field structure, they can also be incorporated using standard techniques \citep{bhattacharya2011sparse,rue2001fast}. The cost-per-iteration of the Gibbs sampler retains the same complexity as in the diagonal $\Sigma$ case; see \S \ref{section:post_comp} for more details. In the diagonal case, we assign independent improper priors $\pi(\sigma_h^2) \propto \sigma_h^{-2}, \,h = 1, \ldots, q$ on the diagonal elements, and call the resulting prior $\Pi_{\Sigma}$. 

The model augmented with the above priors now takes the shape 
\begin{align}
&Y = XBA^\T+E, \quad e_i \sim \Gauss(0,\Sigma), \label{eq:c_factor} \\
& B \sim \Pi_B, \quad A \sim \Pi_A, \quad \Sigma \sim \Pi_{\Sigma}. \label{eq:prior}
\end{align}
We shall refer to the induced prior on $C = B A^{\T}$ by $\Pi_C$. 

\subsection{Concentration results}
In this section, we establish a minimax posterior concentration result under the prediction loss when the number of covariates are allowed to grow sub-exponentially in $n$. To best of our knowledge, this is the first such result in Bayesian reduced rank regression models. We are also not aware of a similar result in any model involving the horseshoe or another polynomial tailed shrinkage prior in ultrahigh-dimensional settings. \cite{armagan2013posterior} applied the general theory of posterior consistency \citep{ghosal2000} to linear models with growing number of covariates and established consistency for the horseshoe prior with a sample size dependent hyperparameter choice when $p = o(n)$. Results \citep{van2014horseshoe,ghosh2017asymptotic} that quantify rates of convergence focus {\em exclusively} on the normal means problem, with their proofs crucially exploiting an exact conjugate representation of the posterior mean. 

%Despite the empirical success of the horseshoe and other shrinkage priors, theoretical investigations into their large sample frequentist properties have been mostly limited to the normal means problem , where a specific conjugate nature of the posterior can be exploited to directly analyze the frequentist properties of the posterior mean. The general theory for posterior consistency and concentration \citep{ghosal2000} requires a delicate balance between complexity of the parameter space and the prior tails, which leads to sub-optimal rates of convergence when applied to heavy-tailed priors like the horseshoe. \cite{armagan2013posterior} applied the general theory to linear models with growing number of covariates and established posterior consistency for the horseshoe prior with a sample size dependent hyperparameter choice when $p = o(n)$. 
%%Since they used concentration properties of the ordinary least squares estimator, it is not clear whether their theory can be extended to $p > n$ settings. 
%We are not aware of a posterior concentration result for the horseshoe prior in $p > n$ settings for any statistical model. 
%Indeed, \cite{pati2014posterior} is the 
%Indeed, barring the Dirichlet--Laplace prior \citep{pati2014posterior}, we are not aware of any convergence results for any shrinkage prior in the $p > n$ setting. 

In this article, we study concentration properties of the posterior distribution in model \eqref{eq:c_factor} in $p \gg n$ settings. To aid the theoretical analysis, we adopt the fractional posterior framework of \cite{bhattacharya2016bayesian}, where a fractional power of the likelihood function is combined with a prior using the usual Bayes formula to arrive at a fractional posterior distribution. It is sufficient to verify a prior mass condition on Kullback--Leibler neighborhoods of the true parameter for the fractional posterior to concentrate at the true parameter, obviating the construction of sieves and verification of prior tail conditions to prove consistency of the actual posterior \citep{ghosal2000}. Our key contribution is the development of a novel non-asymptotic prior concentration result for the horseshoe prior around sparse vectors in Lemma S1 in the supplemental document. The only other instance of a similar prior concentration result in $p \gg n$ settings that we are aware of is for the Dirichlet--Laplace prior \citep{pati2014posterior}. We believe our prior concentration result will be of independent interest in various other models involving the horseshoe prior. We first list the assumptions needed for our concentration result. 
\begin{assumption}[Growth of number of covariates]
$\log p/n \le \gamma$ for some $\gamma \in (0,1)$. 
\end{assumption}
%\begin{assumption}[Entries of design matrix] The elements of the design matrix $X$ are generated from a distribution that is independent of $r,p,q$. 
%\end{assumption}
\begin{assumption}[Response covariance]
%The covariance matrix of $e_i$ is $\mathrm{I}_q$, that is $\Sigma=\mathrm{I}_q$.
The error covariance $\Sigma$ is known, and without loss of generality, assume $\Sigma = \mathrm{I}_q$. 
\end{assumption}
\begin{assumption}[True coefficient matrix] The true coefficient matrix $C_0$ admits the decomposition $C_0=B_0A_0^\T$ where $B_0 \in \Re^{p\times r_0}$ and $A_0 \in \Re^{q \times r_0}$ for some $r_0=\kappa q$, $\kappa\in (0,1]$. We additionally assume that $A_0$ is semi-orthogonal, i.e. $A_0^\T A_0= \mathrm{I}_{r_0}$, and all but $s$ rows of $B_0$ are identically zero for some $s=o(p)$. 
%The true rank and number of responses satisfy $r_0=O(q)$ and $q=O(s\log p)$.
Finally, $\underset{j,h}{\max}\mid C_{0jh}\mid <T$ for some $T > 0$. 
\end{assumption}
\begin{assumption} The number of response variables $q$ is fixed. 
%$q = O(\min\{s \log p, n\})$.  
\end{assumption}
\begin{assumption}[Design matrix]
For $X_j$ the $j$th column of $X$, $\max_{1 \le j \le p} \norm{X_j} \asymp n$. 
\end{assumption}


Assumption 1 allows the number of covariates $p$ to grow at a sub-exponential rate of $e^{n^{\gamma}}$ for some $\gamma \in (0, 1)$. Assumption 2 is made for technical convenience and is not critical as $\Sigma$ is only a $q \times q$ matrix. Assumption 3 posits that the true coefficient matrix $C_0$ admits a reduced-rank decomposition with the matrix $B_0$ row-sparse. Assumption 4 can be relaxed to let $q$ grow slowly with $n$. Assumption 5 is a standard minimal assumption on the design matrix, and is satisfied with large probability if the elements of the design matrix are independently drawn from a fixed probability distribution, such as $\mbox{N}(0, 1)$. 

Fix $\alpha \in (0, 1)$. Recall the prior $\Pi_C$ on $C$ defined after equation \eqref{eq:prior}. The $\alpha$-fractional posterior for $C$ under model \eqref{eq:c_factor} is given by 
\begin{align}\label{eq:frac_post}
\Pi_{n, \alpha}(C) \propto e^{- \alpha \, \mbox{tr}\{ (Y - X C) (Y - X C)^{\T} \}/2 } \, \Pi_C(C), 
\end{align}
where we used Assumption 2 to set $\Sigma  = \mathrm{I}_q$. Assuming the data is generated with a true coefficient matrix $C_0$, we wish to show that the $\alpha$-fractional posterior $\Pi_{n, \alpha}(\cdot)$ concentrates around $C_0$ under the averaged prediction loss $\norm{X C - X C_0}_F^2/(n q)$, where $\norm{\cdot}_F^2$ denotes the Frobenius norm of a matrix. Theorem \ref{th:theorem1} states a non-asymptotic concentration bound for $\Pi_{n, \alpha}$ below. 
\begin{theorem}\label{th:theorem1}
Fix $\alpha \in (0,1)$ and let $\Pi_{n, \alpha}$ be defined as in \eqref{eq:frac_post}. Suppose assumptions 1-5 are satisfied. Let $\epsilon_n=\{(q r_0+r_0s\log p)/n\}^{1/2}$. Then for any $D\geq 2$ and $t>0$,
$$\Pi_{n,\alpha}\left(C\in \Re^{p\times q}: \frac{1}{nq}\| XC-XC_0\|_F^2\geq \dfrac{2(D+3t)}{\alpha(1-\alpha)}\epsilon_n^2 \mid Y\right)\leq e^{-tn\epsilon_n^2}$$
holds with probability at least $1-2/\{(D-1+t)n\epsilon_n^2\}$ under the true distribution.
\end{theorem} 
A proof of Theorem \ref{th:theorem1} is provided in the Appendix. An immediate consequence is that the posterior mean $\overline{C} = \int C \, \Pi_{n, \alpha}(dC)$ satisfies $\|X \overline{C} - X C_0\|_F^2/(n q) \lesssim K \epsilon_n^2$ with large probability for some constant $K > 0$. It then follows from \cite{bunea2012joint} that the posterior mean achieves the minimax rate as a frequentist estimate. 

A main ingredient of the proof is a novel non-asymptotic prior concentration result for the prior $\Pi_C$. Building on a sequence of results, we prove in Lemma S4 of the supplementary document that if $C_0$ satisfies Assumption 3, then 
$$
\Pi_C(C: \|C - C_0\|_F < \epsilon_n) \ge e^{- K n \epsilon_n^2},
$$
where $K$ is a positive constant. The implication is that the prior places sufficient mass to minimax neighborhoods of the truth, which drives the concentration of the fractional posterior $\Pi_{n, \alpha}$. 

%En route, we derive a prior concentration result for the horseshoe prior in the supplementary document which can possibly be of independent interest. The only other instance of a similar prior concentration result in $p \gg n$ settings that we are aware of is for the Dirichlet--Laplace prior \citep{pati2014posterior}. 

%horseshoe prior, which is provided in a supplemental document due to its length. 

We expect the concentration results to hold as well for the actual posterior, i.e., when $\alpha = 1$, which, however, is well beyond the scope of the current paper. For model \eqref{eq:c_factor}, raising the likelihood to a fractional power only results in a change in the (co)variance term, and hence our Gibbs sampler can be easily adapted to sample from the fractional posterior. Our numerical investigations found the point estimates to be virtually indistinguishable from the actual posterior when the fractional power $\alpha$ is close to one, e.g., 0.95 or 0.99, justifying our choice of the theoretical environment. It was also theoretically shown in 
\cite{bhattacharya2016bayesian} that as $\alpha\rightarrow 1-$, the fractional posterior converges weakly to the usual posterior. 




%If additional structural information on $\Sigma$ is available in the form of $\Sigma$ 
%We assume that the rows of the error matrix are independent with 
% \Sigma=\text{diag}(\sigma_1^2,\ldots,\sigma_q^2)$, thus allowing the $q$ responses to be measured with different precision. 
%Depending on the application $\Sigma$ can assume different structures; we shall restrict attention to the case $\Sigma = \mbox{diag}(\sigma_1^2, \ldots, \sigma_q^2)$ in this article, noting that extensions to auto-regressive or low rank structures can be accommodated easily, especially when $q$ is large. For instance, we can write $\Sigma=\Lambda \Lambda^\T+D$ where $\Lambda \in \Re^{q\times k}$ for some integer $k<q$ and $D \in \Re^{q\times q}$ is a diagonal matrix. Sparse priors for the resulting factor model are developed in \cite{bhattacharya2011sparse}, \cite{lopes2004bayesian}. When $q$ is below 20, the Inverse-Wishart$(q,\mathrm{I}_q)$ prior on $\Sigma$ can be an adequate alternative. 

%In reduced rank regression problems with $\text{rank}(C)=r\leq \mathrm{min}(p,q)$, $C$ admits a rank decomposition $C = B_*A_*^\T$, where $B_* \in \Re^{p \times r}$ and $A_* \in \Re^{q \times r}$. While it is possible to treat $r$ as a parameter and assign it a prior distribution inside the hierarchical formulation, posterior inference on $r$ either requires calculation of marginal likelihoods, which is intractable for the present problem, or resorting to complicated reversible jump Markov chain Monte Carlo algorithms. As an alternative, we write $C = B A^{\T}$ where $B\in \Re^{p \times q}$ and $A \in \Re^{q \times q}$, and work in a parameter expanded framework \citep{liu1999parameter} with the aim of shrinking out redundant columns in $A$ and $B$. Restricting $A$ to be in the Stiefel manifold of orthogonal matrices forces $B$ to have the same row sparsity as the coefficient matrix \citep{chen2012sparse}. Uniform priors over the Stiefel manifold may then be used and computational strategies for posterior inference with such priors were developed in \cite{hoff2007gibbs}.  We found in our numerical experiments that posterior updates involving uniform priors resulted in relatively slow computation time as the algorithm suggested by \cite{hoff2007gibbs} involves conditional Gibbs update of each column of an orthogonal matrix. To this end, we assign independent Gaussian priors, $\Pi_A$ on the entries of $A$ noting that this prior specification allows block update of $\text{vec}(A)$, and assign shrinkage priors, $\Pi_B$ on the columns of $B$ to encourage column sparsity. Our hierarchical prior is,
%
%\begin{align}
%&Y = XBA^\T+E, \quad e_i \sim \Gauss(0,\Sigma), \label{eq:c_factor}\\
%& \Pi_B: b_{jh} \mid \lambda_{jh},\tau_h \sim \Gauss(0,\lambda_{jh}^2 \tau_h^2), \\
%&\lambda_{jh} \sim f, \quad \tau_{h} \sim g, \quad j=1,\ldots,p, \, h=1,\ldots,q, \\
%& \Pi_A: a_{hk} \sim N(0, 1), \quad h, k=1,\ldots q, \\
%& \pi(\sigma_h^2) \propto \sigma_h^{-2}, \quad h = 1, \ldots, q, 
%\end{align}
%where $f$ and $g$ are densities on the positive real line. Choices of $f$ and $g$ crucially impact the performance in sparse situations \citep{polson2010shrink,pati2014posterior}. We set $f$ and $g$ to be the Half-Cauchy distribution which leads to the Horseshoe prior \citep{carvalho2010horseshoe}. We write the induced prior on $C$ as $\Pi_C$, defined as $\Pi_C=\Pi_B \otimes \Pi_A$.
% 

%\begin{proof}
%See appendix.
%\end{proof}

%Theorem \ref{th:theorem1}, to the best of our knowledge, is the first result on posterior consistency of the Horseshoe prior in general high dimensional regression model. Earlier, \cite{armagan2013posterior} proved posterior consistency for a class of shrinkage priors. However, their results were obtained for$p<n$, limiting its application to high-dimensional problems. Moreover, the rate of contraction was not specified by the authors. Our result characterizes the rate as $qr_0+r_0s\log p$ which is the same minimax rate obtained in \cite{bunea2012joint}. We impose fairly minimal restrictions on the design matrix since the problem of prediction does not require identifiability constraints on the design matrix in high dimensional settings. See section 4 of \cite{castillo2015bayesian} for a more detailed discussion on this. The proof techniques reveal that the rate is optimal upto $\log q$, a price one pays for starting with $C \in \Re^{p\times q}$.  

\section{Posterior Computation}\label{section:post_comp}
Exploiting the conditional conjugacy of the proposed prior, we develop a straightforward and efficient Gibbs sampler to update the model parameters in \eqref{eq:c_factor} from their full conditional distributions. We use vectorization to update parameters in blocks. Specifically, in what follows, we will make multiple usage of the following identity. For matrices $\Phi_1, \Phi_2, \Phi_3$ with appropriate dimensions, and $\text{vec}(A)$ denoting column-wise vectorization, we have, 
\begin{equation}\label{eq:vectorization}
\text{vec}(\Phi_1\Phi_2\Phi_3)=(\Phi_3^\T \otimes \Phi_1)\text{vec}(\Phi_2)=(\Phi_3^\T \Phi_2^\T \otimes I_k)\text{vec}(\Phi_1),
\end{equation}
where the matrix $\Phi_1$ has $k$ rows and $\otimes$ denotes the Kronecker product. 

Letting $\theta\mid -$ denote the full conditional distribution of a parameter $\theta$ given other parameters and the data, the Gibbs sampler cycles through the following steps, sampling parameters from their full conditional distributions: 

\noindent \textbf{Step 1. } To sample $B \mid -$, use \eqref{eq:vectorization} to vectorize $Y=XBA^\T+E$ to obtain,
\begin{equation}\label{eq:vec_lin_reg1}
y =(X\otimes A) \beta + e, 
\end{equation} 
where $\beta = \text{ vec}(B^\T) \in \Re^{pq \times 1}$, $y = \text{vec}(Y^\T) \in \Re^{n q \times 1}$, and $e=\text{vec}(E^\T) \sim \Gauss_{nq}(0,\widetilde{\Sigma})$ with $\widetilde{\Sigma}=\text{diag}(\Sigma,\ldots,\Sigma)$. Multiplying both sides of \eqref{eq:vec_lin_reg1} by $\widetilde{\Sigma}^{-1/2}$ yields $\widetilde{y}=\widetilde{X}\beta+\widetilde{e}$ where $\widetilde{y}=\widetilde{\Sigma}^{-1/2} y$, $\widetilde{X}=\widetilde{\Sigma}^{-1/2}(X \otimes A)$ and $\widetilde{e}=\widetilde{\Sigma}^{-1/2}e \sim \Gauss_{nq}(0,\mathrm{I}_{nq})$. Thus, the full conditional distribution $\beta \mid - \sim \Gauss_{pq}(\Omega_B^{-1}\widetilde{X}^\T\widetilde{y},\Omega_B^{-1})$, where $\Omega_B=(\widetilde{X}^\T \widetilde{X}+\Lambda^{-1})$ with $\Lambda=\text{diag}(\lambda_{11}^2\tau_1^2,\ldots,\lambda_{1q}^2\tau_q^2,\ldots,\lambda_{p1}^2\tau_1^2,\ldots,\lambda_{pq}^2\tau_q^2)$. 

Naively sampling from the full conditional of $\beta$ has complexity $O(p^3q^3)$ which becomes highly expensive for moderate values of $p$ and $q$. \cite{bhattacharya2016fast} recently developed an algorithm to sample from a class of structured multivariate normal distributions whose complexity scales linearly in the ambient dimension. 
%because of the presence of $\Omega_B^{-1} \in \Re^{pq \times pq}$ and complexity of standard algorithms for generating random vectors from multivariate Gaussian distributions scale in the cubic order of the dimension. 
We adapt the algorithm in \cite{bhattacharya2016fast} as follows: 

\,

%\begin{algorithm}[h!]
%\caption{ } 
%\label{alg:al1}
%\vspace*{8pt}
%\begin{tabbing}
 \noindent (i) Sample $u \sim \Gauss(0,\Lambda)$ and $\delta \sim \Gauss(0,\mr I_{nq})$ independently. \\
   \noindent (ii) Set $v=\widetilde{X} u+\delta$.\\
   \noindent (iii) Solve $(\widetilde{X} \Lambda\widetilde{X}^{\T}+\mr I_{nq})w = (\tilde{y} - v)$ to obtain $w$.\\
   \noindent (iv) Set $\beta=u+\Lambda \widetilde{X}^{\T}w$.
%\end{tabbing}
%\vspace*{-2pt}
%\end{algorithm} 

It follows from \cite{bhattacharya2016fast} that $\beta$ obtained from steps (i) - (iv) above produce a sample from the desired full conditional distribution. One only requires matrix multiplications and linear system solvers to implement the above algorithm, and no matrix decomposition is required. It follows from standard results \citep{golub2012matrix} that the above steps have a combined complexity of $O(q^3 \max\{n^2, p\})$, a substantial improvement over $O(p^3 q^3)$ when $p \gg \max\{n, q\}$. 

\,
\noindent \textbf{Step 2.} To sample $A \mid -$, once again vectorize $Y=XBA^\T+E$, but this time use the equality of the first and the third terms in \eqref{eq:vectorization} to obtain,
\begin{equation}\label{eq:vec_lin_reg2}
y =(XB \otimes \mathrm{I}_q) a + e,
\end{equation}
where $e$ and $y$ are the same as in step 1, and $a = \text{vec}(A) \in \Re^{q^2 \times 1}$. The full conditional posterior distribution $a \mid - \sim \Gauss(\Omega_A^{-1}X_*\widetilde{y},\Omega_A^{-1})$,  where $\Omega_A = (X_*^\T X_*+\mathrm{I}_{q^2})$, $X_*=\widetilde{\Sigma}^{-1/2}(XB \otimes \mathrm{I}_{q^2})$ and $\widetilde{y} = \widetilde{\Sigma}^{-1/2} y$. To sample from the full conditional of $a$, we use the algorithm from \S 3.1.2 of  \cite{rue2001fast}. Compute the Cholesky decomposition $(X_*^{\T} X_* + I_{q^2}) = LL^\T$. Solve the system of equations: $L v = X_*^{\T}\tilde{y}$, $L^{\T} m = v$, and $L^{\T} w = z$, where $z \sim \Gauss(0,\mr{I_{q^2}})$. Finally obtain a sample as $a=m+w$.
\,

\noindent \textbf{Step 3.} To sample $\sigma_h^2 \mid -$, observe that $\sigma_h^2 \mid - \sim \text{inverse-Gamma}(n/2,S_h/2)$ independently across $h$, where $S_h=\{ Y_h-(XBA^\T)_h \}^\T \{ Y_h-(XBA^\T)_h \}$, with $\Phi_h$ denoting the $h$th column of a matrix $\Phi$.

\, 

\noindent \textbf{Step 4.} The global and local scale parameters $\lambda_{jh}$'s and $\tau_h$'s have independent conditional posteriors across $j$ and $h$, which can be sampled via a slice sampling scheme provided in the online supplement to \cite{polson2014bayesian}. We illustrate the sampling technique for a generic local shrinkage parameter $\lambda_{jh}$; a similar scheme works for $\tau_h$. 
%We have $\pi(\lambda_{jh} \mid -) \propto \lambda_{jh}^{-1}\exp\big \{- b_{jh}^2/ (2\lambda_{jh}^2\tau_h^2)\big\}\dfrac{1}{1+\lambda_{jh}^2} \propto \exp\big \{-\dfrac{b_{jh}^2}{2\tau_h^2}\eta_{jh}\big\}\dfrac{1}{1+\eta_{jh}},$
%where $\eta_{jh}=1/\lambda_{jh}^2$. 
Setting $\eta_{jh} = \lambda_{jh}^{-2}$, the slice sampler proceeds by sampling $u_{jh} \mid \eta_{jh} \sim \text{Unif}(0,1/(1+\eta_{jh}))$ and then sampling $\eta_{jh} \mid u_{jh} \sim \text{Exp}(2\tau_h^2/b_{jh}^2)\mathrm{I}\{ \eta_{jh}<(1-u_{jh})/u_{jh} \}$, a truncated exponential distribution. 

The Gibbs sampler above can be trivially modified to accommodate non-diagonal $\Sigma$. Steps 1-2 do not assume any structure for $\Sigma$. The matrix $\Sigma^{-1/2}$ can be computed in $O(q^3)$ steps using standard algorithms, which does not increase the overall complexity of steps 1 and 2 since since $q<n\ll p$ by assumption. Step 3 will change depending on the prior on $\Sigma$; for example, a conjugate update is retained under an inverse-Wishart prior on $\Sigma$. Modifications to situations where $\Sigma$ has a graphical/factor model structure is also straightforward. 
%Steps 1-4 in no way assumes any structure on the error covariance matrix $\Sigma$. As a result, for non-diagonal $\Sigma$ adding one more Gibbs step at the end of these four steps to update $\Sigma$ would suffice. In steps 1 and 2, $\Sigma^{-1/2}$ can be computed with $O(q^3)$ complexity by standard matrix factorization algorithms which does not increase the overall complexity of the algorithm since $q<n\ll p$ by assumption. 

Point estimates of $C$, such as the posterior mean, or element-wise posterior median, are readily obtained from the Gibbs sampler along with a natural uncertainty quantification, which can be used for point and interval predictions. However, the continuous nature of our prior implies that such point estimates will be non-sparse and full rank with probability one, and hence not directly amenable for variable selection and rank estimation. Motivated by our concentration result in Theorem \ref{th:theorem1} that the posterior mean $X \overline{C}$ increasingly concentrates around $X C_0$, we propose two simple post-processing schemes for variable selection and rank estimation below. The procedures are completely automated and do not involve any input of tuning parameters from the user's end. 

\subsection{Post processing for variable selection}
%Point estimates $\widehat{C}_B$, posterior mean or element-wise posterior median, of $C$ obtained from the Gibb's sampler above can be readily used for prediction purposes, with a natural quantification of predictive uncertainty. However, the continuous nature of our prior implies that such point estimates will be non-sparse and full rank with probability one. Theorem \ref{th:theorem1} established that the posterior distribution of $X\overline{C}$ will be centered around $XC_0$ with high probability. We expect a similar phenomenon for the posterior distribution of $\overline{C}$ to be centered around $C_0$ and develop a two step automated procedure to arrive at an estimate of $C$ which is row sparse and rank reduced. 
We first focus on variable selection. We define a row-sparse estimate $\widehat{C}_R$ for $C$ as the solution to the optimization problem  
\begin{equation}\label{eq:utility12}
\widehat{C}_R=\argmin_{\Gamma \in \Re^{p \times q}} \bigg\{ \| X\overline{C}-X\Gamma \|_F^2 + \sum_{j=1}^p \mu_j \| \Gamma^{(j)} \|_2 \bigg\},
\end{equation} 
where $\Phi^{(j)}$ represents the $j^{th}$ row of a matrix $\Phi$, and the $\mu_j$s are predictor specific regularization parameters. The objective function aims to find a row-sparse solution close to the posterior mean in terms of the prediction loss, with the sparsity driven by the group lasso penalty \citep{yuan2006model}. For a derivation of the objective function in \eqref{eq:utility12} from a utility function perspective as in \cite{hahn2015decoupling}, refer to the Appendix. 

%To solve \eqref{eq:utility12}, 
%
%
%
%If the true $C_0$ has zero rows, the corresponding rows of $\overline{C}$ are expected to be shrunk aggressively


%The first step is to zero out rows of $\widehat{C}_B$ aimed at removing irrelevant predictors. To achieve such row sparsity, we consider an utility function based approach which minimizes prediction error in future observations with an added penalty for model complexity. Specifically, using the posterior mean $\overline{C}$ as our choice for $\widehat{C}_B$, the row sparse estimate $\widehat{C}_R$ for $C$ is given by ,
%\begin{equation}\label{eq:utility12}
%\widehat{C}_R=\argmin_{\Gamma} \mid \mid X\overline{C}-X\Gamma \mid \mid_F^2 + \sum_{j=1}^p \mu_j \mid \mid \Gamma^{(j)} \mid \mid_2\\
%\end{equation} 
%where $\Phi^{(j)}$ represents the $j^{th}$ row of a matrix $\Phi$, and the $\mu_j$'s are predictor specific regularization parameters. The group lasso \citep{yuan2006model} penalty $\sum_{j=1}^p \mu_j \mid \mid \Gamma^{(j)} \mid \mid_2$ implies row sparsity in the resulting solution. The main difference from a usual group lasso estimator is achieved by replacing the noisy data matrix $Y$ by the fit $X \overline{C}$ from the shrinkage procedure. For a derivation of the objective function in \eqref{eq:utility12} from a utility function perspective as in \cite{hahn2015decoupling}, refer to the Appendix. 

To solve \eqref{eq:utility12}, we set the sub-gradient of \eqref{eq:utility12} with respect to $\Gamma^{(j)}$ to zero and replace $\|\Gamma^{(j)} \|$ by a data dependent quantity to obtain the soft thresholding estimate,
\begin{equation}\label{eq:subgrad2}
\widehat{C}_R^{(j)}=\dfrac{1}{X_j^\T X_j}\left(1-\dfrac{\mu_j}{2\| X_j^\T R_j \|}\right )_+ X_j^\T R_j, 
\end{equation}
where for $x\in \Re, x_+=\max(x,0)$, and $R_j$ is the residual matrix obtained after regressing $X\overline{C}$ on $X$ leaving out the $j^{th}$ predictor, $R_j= X \overline{C} - \sum_{k\neq j}X_k \widehat{C}_R^{(k)}$. For practical implementation, we use $\overline{C}$ as our initial estimate and make a single pass through each variable to update the initial estimate according to \eqref{eq:subgrad2}. With this initial choice, $R_j=X_j\overline{C}^{(j)}$ and $\| X_j^\T R_j\|=\|X_j\|^2 \| \overline{C}_j \|$. Derivation of \eqref{eq:subgrad2} is deferred to the Appendix. 

While the $p$ tuning parameters $\mu_j$ can be chosen by cross-validation, the computational cost explodes with $p$ to search over a grid in $p$ dimensions. Exploiting the presence of an optimal initial estimate in the form 
of $\overline{C}$, we recommend default choices for the hyperparameters as $\widehat{\mu}_j=1/\| \overline{C}_j\|^{-2}$. When predictor $j$ is not important, the minimax $\ell_2$-risk for estimating $C_0^{(j)}$ is $(\log q)/n$, so that $\| \overline{C}^{(j)}\| \asymp (\log q)/n$. Since $\|X_j\|^2 \asymp n$ by assumption, $\widehat{\mu}_j/\| X_j^\T R_j \| \asymp n^{1/2}/(\log q)^{3/2} \gg 1$, implying a strong penalty for all irrelevant predictors.

Following \cite{hahn2015decoupling}, posterior uncertainty in variable selection can be gauged if necessary by replacing $\overline{C}$ with the individual posterior samples for $C$ in \eqref{eq:utility12}. 


%Exploiting the We propose setting $\mu_j=1/\mid\mid \overline{C}_j\mid \mid^{-2}$ as default values of the tuning parameters. The choice is motivated by the initial estimate which is the posterior mean itself. With this initial choice, $R_j=X_j\overline{C}^{(j)}$ and $\mid \mid X_j^\T R_j\mid \mid=(X_j^\T X_j) \mid \mid \overline{C}_j \mid \mid$. When predictor $j$ is not important, the minimax $l_2$-risk for estimating $C_0^{(j)}$ is $(\log q)/n$, so that $\mid \mid \overline{C}^{(j)}\mid \mid \approx (\log q)/n$. By assumption 2, $X_j^\T X_j \approx n$ making $\mu_j/\mid \mid X_j^\T R_j \mid \mid= n^{1/2}/(\log q)^{3/2}>1$ for all irrelevant predictors.
%
%
%
%
%
%
% suggested inspection of a type of posterior summary plots to arrive at sparse estimates in the univariate regression context. 
%%It is not immediately clear how to extend the summary plots technique to matrix-valued coefficients.  
%Instead, we provide a soft thresholding scheme with default choices of thresholds. Specifically, we set the sub-gradient of \eqref{eq:utility12} with respect to $\Gamma^{(j)}$ to zero and replace $\mid \mid \Gamma^{(j)} \mid \mid$ by a data dependent quantity to obtain the soft thresholding estimate,
%\begin{equation}\label{eq:subgrad2}
%\widehat{C}_R^{(j)}=\dfrac{1}{X_j^\T X_j}\left(1-\dfrac{\mu_j}{2\mid \mid X_j^\T R_j \mid \mid}\right )_+ X_j^\T R_j, 
%\end{equation}
%where $R_j$ is the residual matrix obtained after regressing $X\overline{C}$ on $X$ leaving out the $j^{th}$ predictor, $R_j= X \overline{C} - \sum_{k\neq j}X_k \widehat{C}_R^{(k)}$ and for $x\in \Re, x_+=\max(x,0)$. We use $\overline{C}$ as our initial estimate and make a pass through each variable and update the initial estimate according to \eqref{eq:subgrad2}. Derivation of \eqref{eq:subgrad2} is deferred to the Appendix. 




\subsection{Post processing for rank estimation}\label{subsec:rank}
%As $\text{rank}(XC)\leq \min \{\text{rank}(X),\text{rank}(C)\}=\min (n,r)=r$, the singular values of $X\widehat{C}_R$ may show an elbow effect: few large singular values and the rest small but non-zero singular values. Hence to 
To estimate the rank, we threshold the singular values of $X \widehat{C}_R$, with $\widehat{C}_R$ obtained from \eqref{eq:subgrad2}. In situations where row sparsity is not warranted, $\overline{C}$ can be used instead of $\widehat{C}_R$. 
For $s_1, \ldots, s_q$ the singular values of $X \widehat{C}_R$, and a threshold $\omega > 0$, define the thresholded singular values as $\nu_h = s_h \, \mathrm{I}(s_h > \omega)$ for $h = 1, \ldots, q$. We estimate the rank as the number of nonzero thresholded singular values, that is, $\widehat{r}= \sum_{h=1}^q \mathrm{I}(\nu_h > 0) = \sum_{h=1}^q \mathrm{I}(s_h>\omega)$. We use the largest singular value of $Y-X\widehat{C}_R$ as the default choice of the threshold parameter $\omega$, a natural candidate for the maximum noise level in the model.


%
%
%
%consider
%\begin{equation}\label{eq:rank_utility}
%\argmin_{\nu_1,\ldots,\nu_q} \sum_{h=1}^q (s_h-\nu_h)^2+\omega \sum_{h=1}^q \mathrm{I}(\nu_h \neq 0).
%\end{equation}
%The solution of \eqref{eq:rank_utility} can be explicitly written as a hard thresholding scheme: $\nu_h=s_h \mathrm{I}(s_h > \omega)$. We estimate the rank as the number of nonzero thresholded singular values, that is, $\hat{r}=\sum_{h=1}^q \mathrm{I}(s_h>\omega)$. We use the largest singular value of $Y-X\widehat{C}_R$ as the default choice of the threshold parameter $\omega$, a natural candidate for the maximum noise level in the model.

\section{Simulation Results}

We performed a thorough simulation study to assess the performance of the proposed method across different settings. For all our simulation settings the sample size $n$ was fixed at $100$. We considered $3$ different $(p, q)$ combinations, $(p,q)=(500,10),(200,30),(1000,12)$. The data were generated from the model $Y = X C_0 + E$. Each row of the matrix $E$ was generated from a multivariate normal distribution with diagonal covariance matrix having diagonal entries uniformly chosen between $0.5$ and $1.75$. The columns of the design matrix $X$ were independently generated from $\Gauss(0,\Sigma_X)$. We considered two cases, $\Sigma_X=\mathrm{I}_p$, and $\Sigma_X=(\sigma_{ij}^X)$, $\sigma_{jj}^X=1$, $\sigma_{ij}^X=0.5$ for $i\neq j$. The true coefficient matrix $C_0 = B_*A_*^\T$, with $B_* \in \Re^{p \times r_0}$ and $A_* \in \Re^{r \times r_0}$, with the true rank $r_0 \in \{3, 5, 7\}$. The entries of $A_*$ were independently generated from a standard normal distribution. We generated the entries in the first $s = 10$ rows of $B_*$ independently from $\mbox{N}(0, 1)$, and the remaining $(p - s)$ rows were set equal to zero. 

%For each of these combinations the first 10 predictors were the active variables, so that the set of active variables $J=\{1,2,\ldots,10\}$, with $p_0=|J|=10$. The data were generated from the model $Y=X B_*A_*^\T+E$, with $B_* \in \Re^{p \times r}$ and $A_* \in \Re^{r \times r}$, with the rank $r \in \{3, 5, 7\}$. The elements of the first $p_0$ rows of $B_*$ and all entries of $A_*$ were independently generated from $N(0, 1)$, while the last $(p-p_0)$ rows of $B_*$ were set to $0$. The true coefficient matrix was set as $C_0=B_*A_*^\T$ . The columns of the design matrix $X$ were independently generated from $\Gauss(0,\Sigma_X)$. We considered two cases, $\Sigma_X=\mathrm{I}_p$, and $\Sigma_X=(\sigma_{ij}^X)$, $\sigma_{jj}^X=1$, $\sigma_{ij}^X=0.5$ for $i\neq j$. Each row of the matrix $E$ was generated from a multivariate normal distribution with diagonal covariance matrix having diagonal entries uniformly chosen between $0.5$ and $1.75$.

As a competitor, we considered the sparse partial least squares (SPLS) approach due to \cite{chun2010sparse}. Partial least squares minimizes the least square criterion between the response $Y$ and design matrix $X$ in a projected lower dimensional space where the projection direction is chosen to preserve the correlation between $Y$ and $X$ as well as the variation in $X$. \cite{chun2010sparse} suggested adding lasso type penalties while optimizing for the projection vectors for  sparse high dimensional problems. Since SPLS returns a coefficient matrix which is both row sparse and rank reduced, we create a rank reduced matrix $\widehat{C}_{RR}$ from $\widehat{C}_R$ for a fair comparison. Recalling that $\widehat{C}_R$ has zero rows, let $\widehat{S}_R$ denote the sub-matrix corresponding to the non-zero rows of $\widehat{C}_R$. Truncate the singular value decomposition of $\widehat{S}_R$ to the first $\widehat{r}$ terms as obtained in \S \ref{subsec:rank}. Insert back the zero rows corresponding to $\widehat{C}_R$ in the resulting matrix to obtain $\widehat{C}_{RR}$. Clearly, $\widehat{C}_{RR} \in \Re^{p \times q}$ so created is row sparse and has rank at most $\widehat{r}$; we shall refer to $\widehat{C}_{RR}$ as the Bayesian sparse multi-task learner (BSML). Matlab implementation of the proposed method can be found online at \href{http://www.stat.tamu.edu/~antik/software.html} {www.stat.tamu.edu/~antik/software.html}. 
 
%Letting $A^{(t)}$ and $B^{(t)}$ denote the value of $A$ and $B$ at the $t^{th}$ iteration of the MCMC chain we define $C^{(t)}=B^{(t)}A^{\T (t)}$. After leaving out the initial burnin samples the posterior mean $\overline{C}$ is set as $=1/T\sum_{t=1}^T C^{(t)}$ . The post processing steps for variable selection and rank estimation are then applied on $\overline{C}$ to obtain $\widehat{C}_{RR}$. We measure the performance of $\widehat{C}_{RR}$ in terms of the mean square error and mean square prediction error defined as $\mathrm{MSE}=\mid \mid \widehat{C}_{RR}-C_0 \mid \mid_F^2/pq \quad \text{and} \quad \mathrm{MSPE}=\mid \mid X\widehat{C}_{RR}-XC_0 \mid \mid_F^2/nq$. The results of the proposed method was compared with Sparse Partial Least Squares due to \cite{chun2010sparse}. Essentially, partial least squares (PLS) minimizes the least square criterion between the response $Y$ and design matrix $X$ in a projected lower dimensional space where the projection direction is chosen to preserve the correlation between $Y$ and $X$ as well as the variation in $X$. \cite{chun2010sparse} suggested adding lasso type penalties while optimizing for the projection vectors for  sparse high dimensional problems. Hereafter we refer their method as SPLS.

For an estimator $\widehat{C}$ of $C$, we consider the mean square error, $\mathrm{MSE}=\| \widehat{C} - C_0 \|_F^2/(pq)$, and the mean square prediction error, $\mathrm{MSPE}=\| X\widehat{C} -XC_0 \|_F^2/(nq)$ to measure its performance. The squared estimation and prediction errors of SPLS and $\widehat{C}_{RR}$ for different settings are reported in table \ref{tab:mse} along with the estimates of rank. In our simulations we used the default 10 fold cross validation in the \texttt{cv.spls} function from the \texttt{R} package \texttt{spls}. The SPLS estimator of the rank is the one for which the minimum cross validation error is achieved. We observed highly accurate estimates of the rank for the proposed method, whereas SPLS overestimated the rank in all the settings considered. The proposed method also achieved superior performance in terms of the two squared errors, improving upon SPLS by as much as 5 times in some cases. Additionally, we observed that the performance of SPLS deteriorated relative to BSML with increasing number of covariates. 

In terms of variable selection, both methods had specificity and sensitivity both close to one in all the simulation settings listed in table \ref{tab:mse}. Since SPLS consistently overestimated the rank, we further investigated the effect of the rank on variable selection. We focused on the simulation case $(p, q, r_0) = (1000, 12, 3)$, and fit both methods with different choices of the postulated rank between 3 and 9. For the proposed method, we set $q^\ast$ in \S \ref{sec:mod} to be the postulated rank, that is, considered $B \in \Re^{p \times q^*}$ and $A \in \Re^{q \times q^\ast}$ for $q^\ast \in \{3, \ldots, 9\}$. For SPLS, we simply input $q^\ast$ as the number of hidden components inside the function \texttt{spls}. Figure \ref{fig:sens_spec} plots the sensitivity and specificity of BSML and SPLS as a function of the postulated rank. While the specificity is robust for either method, the sensitivity of SPLS turned out to be highly dependent on the rank. The left panel of figure \ref{fig:sens_spec} reveals that at the true rank, SPLS only identifies $40 \%$ of the significant variables, and only achieves a similar sensitivity as BSML when the postulated rank is substantially overfitted. BSML, on the other hand, exhibits a decoupling effect wherein the overfitting of the rank does not impact the variable selection performance. 

%Variable selection performance of the proposed method and SPLS is measured in terms of the sensitivity and specificity. The sensitivity is defined as the ratio of the correctly selected variables and the true number of variables and the specificity is defined as the ratio between the number of variables removed and the true number of noise variables. Both methods had specificity and sensitivity close to 1 in all the simulation settings listed in table \ref{tab:mse}.  Since SPLS consistently overestimated the rank of the coefficient matrix we looked at the variable selection performance of both methods under rank misspecification. To that end, we carried out a simulation experiment with dimensions fixed at $(p,q)=(1000,12)$ and the true rank of coefficient matrix set to $r_0=3$. Matrices $C_0$, $X$ and $Y$ were generated as described previously. For BSML, we define the postulated rank as the number of columns of $B$ in the initial rank decomposition $C=BA^\T$. Prior (\ref{eq:HS}) is then modified accordingly by restricting the index $h$ in (\ref{eq:HS}) to vary from 1 to the postulated rank. For example, the postulated rank is $3$ if $C=BA^\T$ with $B \in \Re^{p\times 3}$, $A \in \Re^{q\times 3}$, $b_{jh} \sim (\ref{eq:HS})$ with $h=1,2,3$. The postulated rank for SPLS is the number of hidden components provided to the function \texttt{spls}. In the current simulation setting the postulated rank for both methods was varied from 3 to 9 with increments of 1. We computed the sensitivity and specificity across $50$ replicates for each postulated rank. The left panel of figure \ref{fig:sens_spec} shows that for SPLS variable selection is highly correlated to the postulated rank; to achieve good variable selection SPLS requires an overestimated rank. However, the proposed method exhibit a decoupling effect wherein the rank misspecification does not effect the variable selection performance.

We conclude this section with a simulation experiment carried out in a correlated response setting. Keeping the true rank $r_0$ fixed at 3, the data were generated similarly as before except that the individual rows $e_i$ of the matrix $E$ was generated from $\Gauss(0,\Sigma)$, with $\Sigma_{ii} = 1, \Sigma_{ij} = 0.5, 1\le i \ne j \le p$. 
To accommodate the non-diagonal error covariance, we placed a inverse-Wishart$(q,\mathrm{I}_q)$ prior on $\Sigma$. 
%We report the performance of the BSML estiamtor $\widehat{C}_{RR}$. 
An associate editor pointed out the recent article \citep{ruffieux2017efficient} which used spike-slab priors on the coefficients in a multiple response regression setting. They implemented a variational algorithm to posterior inclusion probabilities of each covariate, which is available from the \texttt{R} package \texttt{locus}. To select a model using the posterior inclusion probabilities, we used the median probability model \citep{barbieri2004optimal}; predictors with a posterior inclusion probability less than 0.5 were deemed irrelevant. We implemented their procedure with the prior average number of predictors to be included in the model conservatively set to 25 which worked favorably for their procedure. We observed a fair degree of sensitivity to this parameter, which when set to the true value $10$, resulted in comparatively poor performance. Table \ref{tab:var_sel} reports sensitivity and specificity of this procedure and ours, averaged over 50 replicates. While the two methods performed almost identically in the relatively low dimensional setting $(p,q)=(200,30)$, BSML consistently outperformed \cite{ruffieux2017efficient} when the dimension was higher. 


%which addressed the problem of variable selection in a multiple response regression setting. \cite{ruffieux2017efficient} uses spike slab prior on the coefficients and a variational Bayes algorithm to provide posterior inclusion probabilities of each covariate. To arrive at a model we used the median probabilty model \citep{barbieri2004optimal}; predictors with a posterior inclusion probabilty less than 0.5 were termed irrelevant. Their method is available to implement from the \texttt{R} package \texttt{locus}. In the current simulation setting we implemented it  with a prior average of predictors to be included in the model set to 25. Table \ref{tab:var_sel} reports sensitivity and specificity of the two methods averaged over 50 replicates. While the two methods perform almost equally well in the relatively low dimensional setting $(p,q)=(200,30)$, BSML consistently outperformed \cite{ruffieux2017efficient} in high dimensional settings $(p,q)=(500,10), (1000,12)$. We observed that  \cite{ruffieux2017efficient} is highly sensitive to the prior average number of predictors provided by the user with its performance deteriorating when this argument was set to the true model size $s=10$.

\begin{center}
\begin{table}[!h]
\caption{Estimation and predictive performance of the proposed method (BSML) versus SPLS across different simulation settings. We report the average estimated rank ($\hat{r}$), Mean Square Error, MSE ($\times 10^{-4}$) and Mean Square Predictive Error, MSPE, across 50 replications. For each setting the true number of signals were 10 and sample size was 100. For each combination of $(p,q,r_0)$ the columns of the design matrix were generated from $\Gauss(0,\Sigma_X)$. Two different choices of $\Sigma_X$ was considered. $\Sigma_X=\mathrm{I}_p$ (independent) and $\Sigma_X=(\sigma_{ij}^X)$,$\sigma_{jj}^X=1$,$\sigma_{ij}^X=0.5$ for $i\neq j$ (correlated). The method achieving superior performance for each setting is highlighted in bold.}
%\vspace{0.2in}
\huge
\centering
%\begin{flushleft}
\scalebox{0.43}{
\begin{tabular}{cccccccccccccc}\toprule
 &  & \multicolumn{12}{c}{{\bf (p,q)}}\\
 \cmidrule{3-14}
&  & \multicolumn{4}{c}{{\bf (200,30)}} & \multicolumn{4}{c}{{\bf (500,10)}} & \multicolumn{4}{c}{{\bf (1000,12)}}\\
 \cmidrule{3-14}
 &  & \multicolumn{2}{c}{\bf Independent} & \multicolumn{2}{c}{\bf Correlated} & \multicolumn{2}{c}{\bf Independent} & \multicolumn{2}{c}{\bf Correlated} & \multicolumn{2}{c}{\bf Independent} & \multicolumn{2}{c}{\bf Correlated}\\
 \cmidrule{3-14}
{{\bf Rank}} & \hspace{0.2in}{{\bf Measures}} & BSML & SPLS & BSML & SPLS & BSML & SPLS & BSML & SPLS & BSML & SPLS & BSML & SPLS \\
\cmidrule(lr){1-14}  %\cmidrule(lr){2-2}  \cmidrule{3-8}  \cmidrule{5-6}  \cmidrule{7-8}\\
& $\hat{r}$ & ${\bf 3.0}$ & $7.9$ & ${\bf 3.0}$  & $9.4$ & ${\bf 3.0}$ & $9.7$ & ${\bf 3.0}$ & $8.8$ & ${\bf 3.2}$ & $9.4$ & ${\bf 3.4}$ & $8.9$ \\
{\bf 3} & $\text{MSE}$ & ${\bf 3}$ & $14$ & ${\bf 5}$ & $15$ & ${\bf 3}$ & $7$ & ${\bf 5}$ & $30$ &  ${\bf 3}$ & $50$ & ${\bf 3}$ & $38$ \\
%${\bf 3}$ & $\text{MSE}_\text{me} $ & $3$ & & $5$ & & $2$ & & $4$ %& & $3$ & & $5$ & & $1$ & & $2$ & \\
& $\text{MSPE}$ & ${\bf 0.07}$ & $0.25$ & ${\bf 0.06}$ & ${\bf 0.17}$ & $0.22$ & ${\bf 0.15}$ & $0.34$ & $0.21$ & ${\bf 0.35}$ & $4.19$ & ${\bf 0.30}$ & $1.51$\\
% & $\text{MSPE}_\text{me}$ & $0.08$ & & $0.09$ & & $0.09$ & & $0.11$ & & $0.05$& & $0.05$ & & $0.12$ & & $0.16$ &\\
\cmidrule{1-14}\\
& $\hat{r}$ & ${\bf 5}$ & $9.7$ &  ${\bf 4.9}$ & $12.2$ & ${\bf 4.9}$ & $9.9$ & ${\bf 4.8}$ & $9.8$  & ${\bf 5.1}$ & $9.9$ & ${\bf 5.1}$ & $9.9$\\
{\bf 5}& $\text{MSE}$ & ${\bf 5}$ & $69$ & ${\bf 9}$ & $61$ & ${\bf 3}$ & $10$ & ${\bf 6}$ & $24$ & ${\bf 2}$ & $108$ & ${\bf 4}$ & $129$\\
%${\bf 5}$ & $\text{MSE}_\text{me} $ & $4$ & & $7$ & & $2$ & & $5$ & & $5$ & & $8$ & & $1$ & & $3$ &  \\
& $\text{MSPE}$ & ${\bf 0.11}$ & $3.8$ & ${\bf 0.09}$ & $4.6$ & ${\bf 0.17}$ & $0.41$ & ${\bf 0.20}$ & $0.38$ & ${\bf 0.32}$ & $9.54$ & ${\bf 0.32}$ & $4.63$\\
%& $\text{MSPE}_\text{me} $ & $0.10$ & & $0.12$ & & $0.11$ & & $0.15$ & & $0.09$ & & $0.08$ & & $0.14$ & & $0.19$ & \\
\cmidrule{1-14}\\
& $\hat{r}$ & ${\bf 6.9}$ & $10.3$ & ${\bf 6.9}$ & $15.8$ & ${\bf 6.8}$ & $10$ & ${\bf 6.7}$ & $9.7$ & ${\bf 6.8}$ & $10.2$ & ${\bf 6.6}$ & $11.5$\\
{\bf 7} & $\text{MSE}$ & ${\bf 6}$ & $116$ & ${\bf 10}$ & $112$ & ${\bf 3}$ & $20$ & ${\bf 5}$ & $49$ & ${\bf 2}$ & $195$ & ${\bf 4}$ & $261$ \\
%${\bf 7}$ & $\text{MSE}_\text{me}$ & $4$ & & $7$ & & $2$ & & $5$ & & $5$ & & $9$ & & $1$ & & $3$ &  \\
& $\text{MSPE} $ & ${\bf 0.12}$ & $10.81$ & ${\bf 0.11}$ & $9.01$ & ${\bf 0.16}$ & $0.72$ & ${\bf 0.16}$ & $0.92$ & ${\bf 0.32}$ & $16.70$ & ${\bf 0.31}$ & $7.44$ \\
% & $\text{MSPE}_\text{me} $ & $0.11$ & & $0.12$ & & $0.11$ & & $0.13$ & & $0.10$ & & $0.10$ & & $0.15$ & & $0.18$ & \\
 \hline
 \end{tabular}
}
%\end{flushleft}
\label{tab:mse}
\end{table}
\end{center}
\begin{center}
\begin{table}[!h]
\caption{Variable selection performance of the proposed method in a non-diagonal error structure setting; $e_i\sim \Sigma$, $\sigma_{ii}=1,\sigma_{ij}=0.5$. Sensitivity and specificity of BSML is compared with \cite{ruffieux2017efficient}.}
\scalebox{0.88}{
\begin{tabular}{ccccccc}\toprule
&  & & \multicolumn{2}{c}{BSML} & \multicolumn{2}{c}{\cite{ruffieux2017efficient}}\\
\cmidrule{1-7}\\
&  ${\bf (p,q)}$ & {\bf Measure}& {\bf Independent} & {\bf Correlated} & {\bf Independent} & {\bf Correlated}\\
 \cmidrule{1-7}\\
 \multirow{9}{*}{${\bf r_0=3}$} & {\bf \multirow{2}{*}{(200,30)}} & Sensitivity & 1 & 1 & 0.96 & 0.87  \\
  & & Specificity & 0.90 & 0.84 & 0.77 & 0.67\\
  \cmidrule{1-7}\\
 & {\bf \multirow{2}{*}{(500,10)}} &   Sensitivity & 1 & 0.99 & 0.9 & 0.8\\
   & & Specificity & 0.99 & 0.99 & 0.80 & 0.64\\
   \cmidrule{1-7}\\
 & {\bf \multirow{2}{*}{(1000,12)}} &   Sensitivity & 0.99 & 0.99 & 0.92 & 0.63\\
   & & Specificity & 0.99 & 0.99 & 0.80 & 0.64\\ 
   \hline  
\end{tabular}
}
\label{tab:var_sel}
\end{table}
\end{center}

%\begin{center}
%\begin{table}
%\huge
%\centering
%\begin{flushleft}
%\scalebox{0.43}{
%\begin{tabular}{ccccccc}\toprule
% &  & & \multicolumn{2}{c}{BSML} & \multicolumn{2}{c}{\cite{ruffieux2017efficient}}\\
% \cmidrule{1-7}\\
%& & & Independent & Correlated & Independent & Correlated\\
% \cmidrule{1-7}\\
%\multirow{3}{c}{r} & & & & & & \\
% \end{tabular}
%}
%\end{flushleft}
%\label{tab:mse}
%\end{table}
%\end{center}

\begin{figure}
\centering
%\begin{flushleft}
\caption{Average sensitivity and specificity across 50 replicates is plotted for different choices of the postulated rank. Here $(p,q,r_0)=(1000,12,3)$. Values for BSML (SPLS) are in bold (dashed).  }

\includegraphics[width=0.95\textwidth]{sens_spec_plot.eps}
%\end{flushleft}
\label{fig:sens_spec}
\end{figure}


\section{Yeast Cell Cycle Data}
%{\tcr Antik: Don't use unnecessary abbreviations for bka. So replace TF by transcription factors below. Also, even if you abbreviate, should be TFs, not TF's.}
Identifying transcription factors which are responsible for cell cycle regulation is an important scientific problem \citep{chun2010sparse}. The yeast cell cycle data from  \cite{spellman1998comprehensive} contains information from three different experiments on mRNA levels of 800 genes on an $\alpha$-factor based experiment. The response variable is the amount of transcription (mRNA) which was measured every 7 minutes in a period of 119 minutes, a total of 18 measurements $(Y)$ covering two cell cycle periods. The ChIP-chip data from \cite{lee2002transcriptional} on chromatin immunoprecipitation contains the binding information of the 800 genes for 106 transcription factors $(X)$. We analyze this data available publicly from the \texttt{R} package \texttt{spls} which has the above information completed for 542 genes. The yeast cell cycle data was also analyzed in \cite{chen2012sparse} via sparse reduced rank regression (SRRR). Scientifically 21 transcription factors of the 106 were verified by \cite{wang2007group} to be responsible for cell cycle regulation.


The proposed BSML procedure identified 33 transcription factors. Corresponding numbers for SPLS and SRRR were 48 and 69 respectively. Of the 21 verified transcription factors, the proposed method selected 14, whereas SPLS and SRRR selected 14 and 16 respectively.
%The effect of the 21 confirmed transcription factors as estimated from our procedure are drawn in figure \ref{fig:yeast1}. 
10 additional transcription factors that regulate cell cycle were identified by \cite{lee2002transcriptional}, out of which 3 transcription factors were selected by our proposed method. 
%We do not enumerate common transcription factors that are selected by SPLS, SRRR and the BSML for the sake of brevity. Some examples include, ACE2, SWI4 etc. 
Figure \ref{fig:yeast1} plots the posterior mean, BSML estimate $\widehat{C}_{RR}$, and 95 \% symmetric pointwise credible intervals for two common effects ACE2 and SW14 which are identified by all the methods. 
%Estimated effects as per the posterior mean, 95\% symmetric credible interval and BSML estimate of ACE2 and SWI4 are plotted in \ref{fig:yeast1}. 
Similar periodic pattern of the estimated effects are observed as well for all the other two methods in contention, perhaps unsurprisingly due to the two cell cycles during which the mRNA measurements were taken. Similar plots for the remaining 19 effects identified by our method are placed inside the supplemental document. 

The proposed automatic rank detection technique estimated a rank of 1 which is significantly different from SRRR $(4)$ and SPLS $(8)$. The singular values of $Y-X\widehat{C}_R$ showed a significant drop in magnitude after the first four values which agrees with the findings in \cite{chen2012sparse}. The 10-fold cross validation error with a postulated rank of 4 for BSML was 0.009 and that of SPLS was 0.19. 

We repeated the entire analysis with a non-diagonal $\Sigma$, which was assigned an inverse-Wishart prior. No changes in the identification of transcription factors or rank estimation were detected.  


%The rank estimator depends on the validity of the linear model assumed between the responses and the covariates and also the error covariance structure. We suspect that although the linear model assumption might be valid here, the error covariance matrix $\Sigma$ should not be assumed identity since the mRNA measurements were taken over a period of 120 minutes with equal intervals. However, the singular values of $Y-XC_R$ for the yeast data showed a significant drop in magnitude after the first four values which agrees with the findings in \cite{chen2012sparse}. The 10-fold cross validation error with a postulated rank at 4 for BSML was 0.009 and that of SPLS was 0.19.

\begin{figure}
\centering
\includegraphics[width=0.95\textwidth]{yeast_plot_new.eps}
\caption{Estimated effects of ACE2 and SWI4, two of 33 transcription factors with non-zero effects on cell cycle regulation. Both have been scientifically verified by \cite{wang2007group}. Dotted lines correspond to 95\% posterior symmetric credible intervals, bold lines represent the posterior mean and the dashed lines plot values of the BSML estimate $\widehat{C}_{RR}$.  }
\label{fig:yeast1}
\end{figure}

%\begin{appendix}
%\appendix
%\appendixone
\section*{Appendix}
\subsection*{Proof of theorem \ref{th:theorem1}}
\begin{proof}
For $C \in \Re^{p \times q}$, we write $p_C^{(n)}$ to denote the density of $Y|X$ which is proportional to $e^{- \alpha \, \mbox{tr}\{ (Y - X C) (Y - X C)^{\T} \}/2 }$. For any $C^*\in \Re^{p \times q}$ we define a $\epsilon$-neighborhood as,
\begin{equation}
B_n(C^*,\epsilon)=\left\lbrace C \in \Re^{p \times q}: \int p_{C^*}^{n} \log (p_{C^*}^{(n)}/p_{C}^{(n)})d\mu^{(n)}\leq n\epsilon^2,\int p_{C^*}^{n} \log^2 (p_{C^*}^{(n)}/p_{C}^{(n)})d\mu^{(n)}\leq n\epsilon^2\right\rbrace,
\end{equation}
and finally we let the R\'enyi divergence of order $\alpha$ as,
\begin{equation}
D_\alpha^{(n)}(C,C_0)=D_\alpha^{(n)}(p_C^{(n)},p_{C_0}^{(n)})=\frac{1}{\alpha-1} \log \int \{p_C^{(n)}\}^\alpha \{p_{C_0}^{(n)}\}^{1-\alpha} d\mu^{(n)}
\end{equation}
We start by characterizing $B_n(C_0,\epsilon)$ and $D_\alpha^{(n)}(C,C_0)$ in terms of $X$ and $C$. It is easy to observe that $B(C_0,\epsilon)\supset A(C_0,\epsilon)=\left\lbrace C\in \Re^{p\times q}:  \frac{1}{nq}\| XC-XC_0\|_F^2\leq \epsilon^2\right\rbrace$ for all $\epsilon>0$. We also have, 
$$D_\alpha^{(n)}(C,C_0)=\frac{\alpha}{2}\| XC-XC_0\|_F^2.$$

We now provide an upper bound to the mean prediction error $\|XC-XC_0\|_F^2$. We have,
\begin{align*}
\|X(C-C_0)\|_F^2\leq \|X\|_F^2 \|C-C_0\|_F^2\leq q \underset{1\leq j \leq p}{\max}\|X_j\|^2\|C-C_0\|_F^2=nq\|C-C_0\|_F^2, 
\end{align*}
where the first inequality follows from the Cauchy-Schwartz inequality and the last equality holds due to assumption 5, for a sufficiently large $n$.
%\begin{align*}
%\mid \mid XC-XC_0\mid \mid_F^2&=\mid \mid \widetilde{X}\gamma-\widetilde{X}\gamma_0\mid \mid_2^2= \mid \mid\sum_{j=1}^{pq}\widetilde{X}_j(\gamma_j-\gamma_{j0})\mid \mid_2^2\\
%&\leq \sum_{j=1}^{pq}\mid \mid \widetilde{X}_j(\gamma_j-\gamma_{j0})\mid \mid_2^2
% \leq \underset{1\leq j \leq pq}{\max}\mid \mid \widetilde{X}_j\mid \mid_2^2\mid \mid \gamma-\gamma_0\mid \mid_2^2=\mid \mid \widetilde{X}\mid \mid_2^2 \mid \mid C-C_0\mid \mid_F^2,
%\end{align*}
%where $\mid \mid \widetilde{X}\mid \mid=\underset{1\leq j \leq pq}{\max}\mid \mid \widetilde{X}_j\mid \mid_2^2$. Since by assumption the entries of the design matrix are independent of $(r,s,p,q)$, $\mid \mid \widetilde{X}\mid \mid $ will typically be of the order $nq$, for example if they are generated from $\Gauss(0,1)$ \citep{castillo2015bayesian}. As a result, we can alternatively write $\mid \mid XC-XC_0\mid \mid_F^2\leq M\mid \mid C-C_0\mid \mid_F^2$ for some positive constant $M$.

We proceed with a prior concentration result on $B_n(C_0,\epsilon_n)$ with $\epsilon_n=\{(qr_0+r_0s\log p)/n\}^{1/2}$ . We prove in Lemma S4 in the supplemental document that $\Pi_C(B_n)\geq \Pi_C(A_n)\geq  \Pi_C\left\lbrace C: \| C-C_0\|_F^2\leq \epsilon_n^2/M\right\rbrace\geq e^{-Kn\epsilon_n^2}$ for some positive constant $K$. The statement of the theorem then follows from Theorem 3.2 of Bhattacharya, Pati, Yang (2016).
\end{proof}
\subsection*{Derivation of \eqref{eq:utility12}}
Suppose $Y^* \in \Re^{n \times q}$ be $n$ future observations with design points $X$ so that given $C$, $Y^*$ can be decomposed into $Y^*=XC+E^*$ where $E^*$ has the same distribution as $E$ in \eqref{eq:lin_reg}. We define the utility function in terms of loss of predicting these $n$ new future observations. To encourage sparsity in rows of a coefficient matrix $\Gamma$ that balances the prediction we add a group lasso penalty \citep{yuan2006model} to this utility function. We define the utility function as,

\begin{equation}\label{eq:utility10}
\mathcal{L}(Y^*,\Gamma)=\| Y^*-X\Gamma \|_F^2 + \sum_{j=1}^p \mu_j \| \Gamma^{(j)} \|_2
\end{equation}
where the $p$ tuning parameters $\{\mu_j\}_{j=1}^p$ control the penalty for selecting each predictor variable and $\Phi^{(j)}$ represents the $j^{th}$ row of any matrix $\Phi$. Intuitively we want $\mu_j$ to be small if the $j^{th}$ predictor is important and vice versa.
The expected risk, $\mathbb{E}\{\mathcal{L}(Y^*,\Gamma)\}$, after integrating over the space of all such future observations given $C$ and $\Sigma$, is 
\begin{align}
\mathcal{L}(\Gamma,C,\Sigma)=q \text{ tr}(\Sigma)+ \| XC-X\Gamma \|_F^2 + \sum_{j=1}^p \mu_j \| \Gamma^{(j)} \|_2.
\end{align}

Finally we take expectation of this quantity with respect to $\pi(C,\Sigma \mid Y,X)$ and drop the constant terms to obtain \eqref{eq:utility12}.

\subsection*{Derivation of \eqref{eq:subgrad2}}
We let $\Phi_j$ and $\Phi^{(j)}$ denote the $j^{th}$ column and row of a generic matrix $\Phi$. Using the subgradient of \eqref{eq:utility12} with respect to $\Gamma^{(j)}$ \citep{friedman2007pathwise}, we have 
\begin{equation}\label{eq:subgrad}
2X_j^\T(X\Gamma-X\overline{C})+\mu_j \alpha_j=0, \quad j=1,\ldots,p,
\end{equation}
where $\alpha_j=\Gamma^{(j)}/\| \Gamma^{(j)} \|$ if $\| \Gamma^{(j)} \| \neq 0$ and $\| \alpha_j \| < 1$ when $\| \Gamma^{(j)} \|=0$. For $\Gamma^{(j)}=0$ we can rewrite \eqref{eq:subgrad} as, $2X_j^\T(\sum_{k\neq j}X_k\Gamma^{(k)}-X\overline{C})+\mu_j \alpha_j=0$ which imply that
$\alpha_j=-2X_j^\T R_j/\mu_j$, where $R_j$ is the residual matrix obtained after regressing $X\overline{C}$ on $X$ leaving out the $j^{th}$ predictor, $R_j=X\overline{C}-\sum_{k\neq j}X_k \Gamma^{(k)}$. We can use this to set $\Gamma^{(j)}$ to zero: if $\alpha_j<1$ set $\Gamma^{(j)}=0$. Otherwise we have $2X_j^\T(X_j\Gamma^{(j)}-R_j)+\mu_j \Gamma^{j}/\|\Gamma^{(j)} \|=0$. Solving for $\Gamma^{(j)}$ in the above equation we then get,
\begin{equation}\label{eq:subgrad1}
 \Gamma^{(j)}=\left(X_j^\T X_j+\dfrac{\mu_j}{2\| \Gamma^{(j)} \|}\right )^{-1} X_j^\T R_j.
\end{equation} 
This solution is dependent on the unknown quantity $\mid \mid \Gamma^{(j)} \mid \mid$. However, taking norm on both sides in \eqref{eq:subgrad1} we get a value of $\|\Gamma^{(j)} \|$ which does not involve any unknown quantities: $\| \Gamma^{(j)} \|=(\| X_j^\T R_j \|-\mu_j/2)/X_j^\T X_j$. Substituting this in \eqref{eq:subgrad1} we get, $\Gamma^{(j)}=(X_j^\T X_j)^{-1}\left(1-\mu_j/2\| X_j^\T R_j \|\right ) X_j^\T R_j$.  

Finally, combining the case when $\Gamma^{(j)}=0$, we have \eqref{eq:subgrad2}.





%\end{appendix}


 
\bibliography{sparse_reduced_rank_refs}
\bibliographystyle{biometrika}

\newpage
\setcounter{equation}{0}
\setcounter{page}{1}
\setcounter{table}{1}
\setcounter{section}{0}
\numberwithin{table}{section}
\renewcommand{\theequation}{S.\arabic{equation}}
\renewcommand{\thesubsection}{S.\arabic{section}.\arabic{subsection}}
\renewcommand{\thesection}{S.\arabic{section}}
%\renewcommand{\thepage}{S.\arabic{page}}
\renewcommand{\thetable}{S.\arabic{table}}
%\baselineskip=17pt

%\renewcommand{\theequation}{S\arabic{equation}}
\renewcommand{\thefigure}{S\arabic{figure}}
\renewcommand{\thelemma}{S\arabic{lemma}}
\renewcommand{\bibnumfmt}[1]{[S#1]}
\renewcommand{\citenumfont}[1]{S#1}

\begin{center}
\Large{\bf Supplementary material}
\end{center}

\vskip 10mm
\begin{center}
Antik Chakraborty\\
Department of Statistics, Texas A\&M University, College Station\\
3143 TAMU, TX 77843-3143, USA\\
antik@stat.tamu.edu\\
\hskip 5mm\\
Anirban Bhattacharya\\
Department of Statistics, Texas A\&M University, College Station\\
3143 TAMU, TX 77843-3143, USA\\
anirbanb@stat.tamu.edu\\
\hskip 5mm\\
Bani K. Mallick\\
Department of Statistics, Texas A\&M University, College Station\\
3143 TAMU, TX 77843-3143, USA\\
bmallick@stat.tamu.edu\\

\end{center}
%\author{Antik Chakraborty, Anirban Bhattacharya, Bani K. Mallick}
%\affil{Department of Statistics, Texas A\&M University, College Station, Texas, 77843, USA
%\email{antik@stat.tamu.edu}\email{anirbanb@stat.tamu.edu}\email{bmallick@stat.tamu.edu}}

\section*{Convention}
Equations defined in this document are numbered (S1), (S2) etc, while (1), (2) etc refer to those defined in the main document. Similar for lemma, theorems etc. 
\vspace{-0.2in}
\section{Proof of theorem \ref{th:theorem1}}
We establish few results in the following sequence of lemmas needed to prove Theorem 1. Recall $Ca^+(0, 1)$ denotes the half-Cauchy distribution with density proportional to $(1 + t^2)^{-1}\ind_{(0,\infty)}(t)$. For a vector $\beta \in \Re^p$, let $\Pi_\HS$ denote the prior on $\beta$ defined by the hierarchy,
\begin{align}\label{eq:horseshoe_prior}
\begin{aligned}
&\beta_j \mid \lambda_j, \tau \sim \Gauss(0,\lambda_j^2\tau^2), \quad (j=1,\ldots,p) \\ 
& \lambda_j \sim Ca^+(0,1), \quad (j=1,\ldots,p)  \\ 
& \tau \sim Ca^+(0,1). \\
\end{aligned}
\end{align} 
We first state a novel prior concentration result for the horseshoe prior which bounds from below the probability assigned to an Euclidean neighborhood of a sparse vector.
\begin{lemma}\label{lm:l1}
Suppose $\beta_0 \in l_0[s;p]$ where $l_0[s;p]=\{\beta_0 \in \Re^p: \#(1\leq j < p : \beta_{0j} \neq 0)\leq s\}$. Let $S=\{j:\beta_{0j}\neq 0\}$ so that $|S|\leq s$. Assume $s=o(p)$ and $\log p/n\leq \gamma$ for some $\gamma \in (0,1)$ and $\max\mid \beta_{0j}\mid \leq M$ for some $M>0$ and for $j \in S$. Suppose $\beta \sim \Pi_\HS$. Define $\delta=\{(s\log p)/n\}^{1/2}$. Then we have, for some positive constant $K$,
$$\Pi_\HS(\|\beta-\beta_0\|_2<\delta)\geq e^{-Ks\log p}.$$ 
\end{lemma}
\begin{proof}
Using the conditional formulation of prior (\ref{eq:horseshoe_prior}) we have,
\begin{equation}
\begin{split}
\Pi_\HS(\| \beta-\beta_0\|_2<\delta)&=\int_\tau pr(\| \beta-\beta_0\|_2<\delta\mid \tau) f(\tau) d\tau\\
& \geq \int_{I_{\tau_*}} pr(\| \beta-\beta_0\|_2<\delta\mid \tau) f(\tau) d\tau,
\end{split}
\end{equation}
where $I_{\tau_*}=[\tau_*/2,\tau_*]$ with $\tau_*=\left(s/p\right)^{3/2}(\log p/n)^{1/2}$.
We first provide lower bound of the conditional probability $pr(\| \beta-\beta_0\|_2<\delta\mid \tau \in I_{\tau_*})$ by dividing it into two terms and providing lower bounds for each of them. For $\tau \in I_{\tau^*}$ we have,
\begin{equation}\label{eq:split_prob}
\begin{split}
pr(\| \beta-\beta_0\| <\delta\mid \tau)&\geq pr\left(\|\beta_S-\beta_{0S}\|_2<\frac{\delta}{2}\mid \tau\right) pr\left(\| \beta_{S^c}\|_2<\frac{\delta}{2}\mid \tau\right)\\
&\geq \prod_{j\in S}pr\left(|\beta_j-\beta_{0j}|<\frac{\delta}{2\sqrt{s}}\mid \tau\right) \prod_{j \in S^c}pr\left(|\beta_j|<\frac{\delta}{2\sqrt{p}}\mid \tau\right)\\
& =\left\lbrace pr\left(|\beta_j-\beta_{0j}|<\frac{\delta}{2\sqrt{s}}\mid \tau \right)\right\rbrace^s \left\lbrace pr \left(|\beta_j|<\frac{\delta}{2\sqrt{p}}\mid \tau\right)\right\rbrace^{p-s}.
\end{split}
\end{equation}
%Now we will provide lower bounds for each of the terms $P(|\beta_j-\beta_{0j}|<\dfrac{\delta}{\sqrt{p}})$ and $P(|\beta_j|<\dfrac{\delta}{\sqrt{p}})$. 
We proceed with $pr(|\beta_j|<\delta /2\sqrt{p}\mid \tau)$ for $\tau \in I_{\tau^*}$. We have, with $R=1/\sqrt{2\pi^3}$ and $\delta_*=\delta/2\sqrt{p}$,
\begin{align*}
pr(|\beta_j|<\delta_*\mid \tau)& \geq \frac{R}{2\tau}\int_{-\delta_*}^{\delta_*} \log\left(1+\dfrac{4\tau^2}{\beta_j^2}\right)d\beta_j\\
&=\frac{R}{\tau}\int_{0}^{\delta_*} \log\left(1+\dfrac{4\tau^2}{\beta_j^2}\right)d\beta_j
\geq \frac{R}{\tau}\int_{0}^{\delta_*} \dfrac{4\tau^2/\beta_j^2}{1+4\tau^2/\beta_j^2} d\beta_j\\
& = 4R\tau \int_{0}^{\delta_*}\dfrac{1}{\beta_j^2+4\tau^2} d\beta_j
=2K\arctan{\dfrac{\delta_*}{2\tau}}\\
& \geq 2R\left(\dfrac{\pi}{2}-\dfrac{2\tau}{\delta_*}\right)
\geq R\left(\dfrac{\pi}{2}-\dfrac{2\tau_*}{\delta_*}\right) = R\pi\left(1-\dfrac{4}{\pi}\dfrac{s}{p}\right), 
\end{align*}
where we have used the inequalities $log(1+x)\geq x/(x+1)$ for $x>-1$ and $\arctan{x}>(\pi/2-1/x)$ in succession.

For $pr(\mid \beta_j-\beta_{0j} \mid <\delta_0|\tau)$ with $\tau \in I_{\tau_*}$, we have, letting $\delta_0=2^{-1}\{(\log p)/n\}^{1/2}$,
\begin{align*}
\begin{split}
pr(\mid \beta_j-\beta_{0j}\mid <\delta_0 \mid \tau)&=\sqrt{2}\pi^{-3/2}\int_{\lambda_j}\int_{\mid \beta_j-\beta_0 \mid <\delta_0} \exp\{{-\beta_j^2/(2\lambda_j^2\tau^2)}\} \dfrac{1}{\lambda_j\tau(1+\lambda_j^2)}d\lambda_j d\beta_j\\
& \geq \sqrt{2}\pi^{-3/2}\int_{|\beta_j-\beta_0|<\delta_0} \int_{1/\tau_*}^{4/\tau_*} \exp\{{-\beta_j^2/(2\lambda_j^2\tau^2)}\} \dfrac{1}{\lambda_j\tau(1+\lambda_j^2)}d\lambda_j d\beta_j\\
& \geq (2\pi)^{-3/2} \int_{|\beta_j-\beta_0|<\delta_0} \exp(-\beta_j^2/2) \int_{1/\tau_*}^{4/\tau_*} \dfrac{1}{1+\lambda_j^2} d\lambda_j d\beta_j,\\
\end{split}
\end{align*}
since for $\lambda_j \in [1/\tau_*, 4/\tau_*]$, $1/\lambda_j \tau \geq 1/4$ and $\exp\{{-\beta_j^2/(2\lambda_j^2\tau^2)}\} \geq \exp{(-\beta_j^2/2)}$. Continuing,
\begin{align*}
pr(\mid \beta_j-\beta_{0j} \mid <\delta_0\mid \tau)&\geq (2\pi)^{-3/2} \dfrac{3\tau_*}{16+\tau_*^2}\int_{|\beta_j-\beta_0|<\delta_0} \exp{(-\beta_j^2/2)}d\beta_j\\
& \geq (2\pi)^{-3/2} \dfrac{3\tau_*}{16+\tau_*^2} \exp\{{-(M+1)^2/2}\}\delta_0\\
& \approx  K_*\left(\dfrac{s}{p}\right)^{3/2}\dfrac{\log p}{n},\\
\end{align*}
for some positive constant $K_*$.
Substituting these bounds in \ref{eq:split_prob}, we have,
\begin{equation}
pr(\|\beta-\beta_0\|<\delta \mid \tau \in I_{\tau_*})\geq e^{-Ks\log p},
\end{equation}
where $K$ is a positive constant. The proof is completed upon observing that $pr(\tau\in I_{\tau_*})\approx \tau_*/2$, so that with a slight abuse of notation we get,
\begin{equation}
\Pi_\HS(\|\beta-\beta_0\|<\delta) \geq e^{-Ks\log p},
\end{equation} 
for some positive constant $K$.
\end{proof}

Recall the prior $\Pi_B$ from the main document. If a matrix $B \in \Re^{p \times q} \sim \Pi_B$ then each column of $B$ is a draw from $\Pi_{\HS}$ In the following Lemma we generalize Lemma \ref{lm:l1} to provide a lower bound on the probability the prior $\Pi_B$ assigns to Frobenius neighborhoods of $B_0 \in \Re^{p \times r_0}$. By assumption 3 the $h$-th column of $B_0$, $b_h \in l_0[s;p]$. In order to make the Frobenius neighborhood well defined we append $(q-r_0)$ zero columns to the right of $B_0$ and set $B_{0*}=(B_0\mid O^{p \times (q-r_0)})$.

\begin{lemma}\label{lm:l2}
Let the entries of $B_{0*} \in \Re^{p \times q}$ satisfy $\max \mid B_{0*} \mid \leq M$ for some positive constant $M$ . Suppose $B$ is a draw from $\Pi_B$. Define $\delta_B=\{(r_0 s \log p)/n\}^{1/2}$. Then for some positive constant $K$ we have,
$$\Pi_B(\| B-B_{0*} \|_F <\delta_B)\geq e^{-Kr_0s\log p}.$$
\end{lemma}
\begin{proof} 
Observe that, since by assumption 3, $r_0=\kappa q$ for some $\kappa\in (0,1]$,
\begin{align*}
 \Pi_B(& \|B- B_{0*}\|_F<\delta_B)  \geq \prod_{h=1}^q\Pi_\HS\left(\|b_h-b_{0h}\|_2<q^{-1/2}\delta_B\right)\\
 &=\prod_{h=1}^q\Pi_\HS\bigg[\| b_h-b_{0h}\|_2<\alpha\{(s\log p)/n\}^{1/2}\bigg]\\
 &= \prod_{h=1}^{r_0} \Pi_\HS\bigg[\| b_h-b_{0h}\|_2<\alpha\{(s\log p)/n\}^{1/2}\bigg]\prod_{h=r_0+1}^q \Pi_\HS\bigg[\| b_h-b_{0h}\|_2<\alpha\{(s\log p)/n\}^{1/2}\bigg]\\ 
 &= \prod_{h=1}^{r_0} \Pi_\HS\bigg[\| b_h-b_{0h}\|_2<\alpha\{(s\log p)/n\}^{1/2}\bigg]\prod_{h=r_0+1}^q \Pi_\HS\bigg[\| b_h\|_2<\alpha\{(s\log p)/n\}^{1/2}\bigg].\\
 \end{align*}
From Lemma \ref{lm:l1} we have, $\Pi_\HS[\| b_h-b_{0h}\|_2<\alpha\{(s\log p)/n\}^{1/2}]\geq e^{-K_1s\log p}$. Arguments along the same line of lemma \ref{lm:l1} can also be applied to obtain that, $\Pi_\HS\left(\| b_h \|_2<\alpha\{(\log p)/n\}^{1/2}\right)\geq e^{-K_2 \log p}$ for some positive $K_2$. Combining these two lower bounds in the above display we have,
$$\Pi_B(\| B-B_0\|_F<\delta_B)\geq e^{-K_1r_0s\log p} e^{-K_2(q-r_0)\log p}.$$
Since $r_0=\kappa q$, the result follows immediately with $K=K_1+(1/\kappa-1)K_2$.
\end{proof}

Similar to the previous lemma, the following result provides a lower bound on the probability assigned to Frobenius neighborhoods of $A_0$ by the prior $\Pi_A$. Again we append $(q-r_0)$ columns at the right of $A_0$ and set $A_{0*}=(A_0 \mid O^{q \times (q-r_0)})$.
\begin{lemma}\label{lm:l3}
Suppose the matrix $A \sim \Pi_A$. Let $\delta_A=(qr_0/n)^{1/2}$. Then for some positive constant $K$ we have,
$$\Pi_A(\| A-A_{0*}\|_F<\delta_A)\geq e^{-Kqr_0}.$$
\end{lemma}
\begin{proof}
First we use vectorization to obtain $\Pi_A(\|A-A_0\|_F<\delta_A)=\Pi_A(\|a-a_0\|_2<\delta_A)$, where $a,a_0\in \Re^{q^2}$.
Using Anderson's lemma \citep{bhattacharya2016suboptimality} for multivariate Gaussian distributions, we then have,
\begin{align*}
\Pi_A(\| a-a_0\|_2<\delta_A)&\geq e^{-\| a_0\|^2/2}pr(\| a \|_2<\delta_A/2)\\
& = e^{-r_0/2}pr (\| a \|_2<\delta_A/2).\\
\end{align*}
The quantity $pr(\| a \|<\delta_A/2)$ can be bounded from below as,
\begin{align*}
pr(\| a \|_2 <\delta_A/2)& \geq \{pr( \mid a_j \mid <\delta_A/q)\}^{q^2}
\geq m e^{-\delta_A^2} \left(\delta_A/q\right)^{q^2}\approx e^{-q^2},
\end{align*}
where $m$ is a positive constant. Since $r_0=\kappa q$, it follows that $\Pi_A(\|A-A_0 \|_F<\delta_A)\geq e^{-Kqr_0}$ for some positive $K$.
\end{proof}

Our final result will concern the prior mass assigned to Frobenius neighborhoods of $C\in \Re^{p \times q}$. As in the main document we write $\Pi_C$ for prior on $C=BA^\T$ induced from $\Pi_B$ and $\Pi_A$.

\begin{lemma}\label{lm:l4}
Suppose $C_0$ satisfies assumption 3. Let $C \sim \Pi_C$ with $\Pi_C$ as defined above. Define $\delta_C=\{(qr_0+r_0s\log p)/n\}^{1/2}$. Then for some positive constant $K$,
$$\Pi_C(\|C-C_0 \|_F<\delta_C)\geq e^{-K(qr_0+r_0s\log p)}.$$
\end{lemma}
\begin{proof}
Using the triangle inequality followed by the fact that for two matrices $P$ and $Q$, $\| PQ\|_F\leq s_{\max}(P)\|Q\|_F$ where $s_{\max}(P)$ is the largest singular value of $P$, we have,
\begin{align*}
\|C-C_0\|_F=\| BA^\T-B_0A_0^\T\|_F&=\| BA^\T -B_0A^\T +B_0A^\T -B_0A_0^\T\|_F\\
&=\| (B-B_{0*})A^\T+B_0(A-A_{0*})^\T\|_F\\
&\leq \|(B-B_0)A^\T\|_F+\| B_0(A-A_0)^\T\|_F\\
&\leq s_{\max}(A)\| (B-B_0)\|_F + s_{\max}(B_0) \| A-A_0\|_F,
\end{align*}
where $s_{\max}(A)$ and $s_{\max}(B_0)$ are the largest singular values of $A$ and $B_0$ respectively. From standard random matrix theory it is well known that for a random matrix of dimension $m \times n$ with independent Gaussian entries, the largest singular admits a high probability upper bound; for every $t\geq 0$, $s_{\max(A)}\leq \sqrt{m}+\sqrt{n}+t$ with probability at least $1-2\exp{(-t^2/2)}$ \citep{vershynin2010introduction}. Also since the elements of $B_0$ are bounded, so is $s_{\max}(B_0)$, say by $\xi$. For a sufficiently large positive number $L$ and for $A\in E=\{A: s_{\max}(A)\leq 2\sqrt{q}+L\}$ we then have,

$$\|C-C_0\|_F\leq (2\sqrt{q}+L)\| B-B_{0*}\| _F+\xi\| A-A_0\|_F.$$
Thus we have, $\Pi_C(\| C-C_0\|_F<\delta_C)\geq \Pi_C\{(2\sqrt{q}+L)\| B-B_0\|_F+\xi\| A-A_0\|_F<\delta_C\}$. Since $\delta_C=\{(qr_0+r_0s\log p)/n\}^{1/2}\geq 2^{-1/2}\{(r_0s\log p/n)^{1/2}+(qr_0)^{1/2}\}=2^{-1/2}(\delta_B+\delta_A)$, the probability $\Pi_C(\| C-C_0\|_F<\delta_C)\geq \Pi_B(\|B-B_0\|_F< m_1 \delta_B)\Pi_A(\|A-A_0\|_F<m_2\delta_A)$, where $m_1$ and $m_2$ are positive constants.
%\geq \Pi_B(\|B-B_0\|_F<\delta_C/2\sqrt{q})\Pi_A(\mid \mid A-%A_0\mid \mid_F<\delta_C/2\beta)$. Since $q=O(s\log p)$, it then %follows that, $\Pi_C(\mid \mid C-C_0\mid \mid_F<\delta_C)\geq \Pi_B(\mid \mid B-B_0\mid \mid_F<m_1\delta_A/\sqrt{q})\Pi_A(\mid \mid A-A_0\mid \mid_F<m_2\delta_A/\beta)$, where $\delta_A=(r_0s\log p/n)^{1/2}$ and $\delta_B=(qr_0/n)^{1/2}$ and $m_1$ and $m_2$ are positive constants.

From lemma \ref{lm:l2} it follows that, $\Pi_B(\|B-B_0\|_F<m_1\delta_A)\geq e^{-K_1r_0s\log p}$ and from lemma \ref{lm:l3} we have, $\Pi_A(\| A-A_0\|_F<m_2\delta_A/\beta)\geq e^{-K_2 qr_0}$. Hence $\Pi_C(E\cap \{C:\|C-C_0\|_F<\delta_C\})\geq e^{-K(qr_0+r_0s\log p)}$. Since for two sets $E_1$ and $E_2$, $pr(E_1\cup E_2)\geq pr(E_1)+pr(E_2)-1$, the Lemma is proved.
\end{proof}


\section*{Yeast cell cycle data}
The yeast cell cycle data consists of mRNA measurements $Y$, measured every 7 minutes in a period of 119 minutes. The covariates $X$ are binding information on 106 transcription factors. When applied to this data, the proposed method identified 33 transcription factors out of 106 that driving the variation in mRNA measurements. 14 of the identified transcription factors are among the 21 scientifically verified \citep{lee2002transcriptional}. In the main document we provided estimated effects of two of the 21 scientifically verified transcription factors. Here we plot the estimated effects of the remaining transcriptions factors that were scientifically verified.
\begin{figure}
\begin{flushleft}
%\centering
\caption{Estimated effects of the 19 of 21 scientifically verified transcription factors selected by the proposed method. Effects of other two, viz. ACE2 and SWI4 are included in the main manuscript. Red lines correspond to 95\% posterior symmetric credible intervals, black lines represent the posterior mean and the blue dashed line plots values of the BSML estimate $\widehat{C}_{RR}$.}
\label{fig:yeast2}
\includegraphics[height=6in,width=6in]{yeast_plot1.eps}
\end{flushleft}
\end{figure}
\end{document}
sparse_reduced_rank.bbl
\begin{thebibliography}{51}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi

\bibitem[{Alquier(2013)}]{alquier2013bayesian}
\textsc{Alquier, P.} (2013).
\newblock Bayesian methods for low-rank matrix estimation: Short survey and
  theoretical study.
\newblock In \textit{Algorithmic Learning Theory}. Springer.

\bibitem[{Anderson(1984)}]{anderson1984multivariate}
\textsc{Anderson, T.} (1984).
\newblock Multivariate statistical analysis.
\newblock \textit{Willey and Sons, New York, NY} .

\bibitem[{Anderson(2002)}]{anderson2002specification}
\textsc{Anderson, T.} (2002).
\newblock Specification and misspecification in reduced rank regression.
\newblock \textit{Sankhy{\=a}: The Indian Journal of Statistics, Series A} ,
  193--205.

\bibitem[{Anderson(1951)}]{anderson1951estimating}
\textsc{Anderson, T.~W.} (1951).
\newblock Estimating linear restrictions on regression coefficients for
  multivariate normal distributions.
\newblock \textit{The Annals of Mathematical Statistics} , 327--351.

\bibitem[{Armagan et~al.(2013)Armagan, Dunson, Lee, Bajwa \&
  Strawn}]{armagan2013posterior}
\textsc{Armagan, A.}, \textsc{Dunson, D.~B.}, \textsc{Lee, J.}, \textsc{Bajwa,
  W.~U.} \& \textsc{Strawn, N.} (2013).
\newblock Posterior consistency in linear models under shrinkage priors.
\newblock \textit{Biometrika} \textbf{100}, 1011--1018.

\bibitem[{Babacan et~al.(2011)Babacan, Molina \&
  Katsaggelos}]{babacan2011variational}
\textsc{Babacan, S.~D.}, \textsc{Molina, R.} \& \textsc{Katsaggelos, A.~K.}
  (2011).
\newblock Variational \uppercase{B}ayesian super resolution.
\newblock \textit{IEEE Transactions on Image Processing} \textbf{20}, 984--999.

\bibitem[{Barbieri \& Berger(2004)}]{barbieri2004optimal}
\textsc{Barbieri, M.~M.} \& \textsc{Berger, J.~O.} (2004).
\newblock Optimal predictive model selection.
\newblock \textit{The Annals of Statistics} \textbf{32}, 870--897.

\bibitem[{Bhadra \& Mallick(2013)}]{bhadra2013joint}
\textsc{Bhadra, A.} \& \textsc{Mallick, B.~K.} (2013).
\newblock Joint high-dimensional bayesian variable and covariance selection
  with an application to eqtl analysis.
\newblock \textit{Biometrics} \textbf{69}, 447--457.

\bibitem[{Bhattacharya et~al.(2016{\natexlab{a}})Bhattacharya, Chakraborty \&
  Mallick}]{bhattacharya2016fast}
\textsc{Bhattacharya, A.}, \textsc{Chakraborty, A.} \& \textsc{Mallick, B.~K.}
  (2016{\natexlab{a}}).
\newblock Fast sampling with gaussian scale mixture priors in high-dimensional
  regression.
\newblock \textit{Biometrika} \textbf{103}, 985--991.

\bibitem[{Bhattacharya \& Dunson(2011)}]{bhattacharya2011sparse}
\textsc{Bhattacharya, A.} \& \textsc{Dunson, D.} (2011).
\newblock {Sparse {B}ayesian infinite factor models}.
\newblock \textit{Biometrika} \textbf{98}, 291--306.

\bibitem[{Bhattacharya et~al.(2016{\natexlab{b}})Bhattacharya, Dunson, Pati \&
  Pillai}]{bhattacharya2016suboptimality}
\textsc{Bhattacharya, A.}, \textsc{Dunson, D.~B.}, \textsc{Pati, D.} \&
  \textsc{Pillai, N.~S.} (2016{\natexlab{b}}).
\newblock Sub-optimality of some continuous shrinkage priors.
\newblock \textit{Stochastic Processes and their Applications} \textbf{126},
  3828 -- 3842.
\newblock In Memoriam: Evarist Giné.

\bibitem[{Bhattacharya et~al.(2016{\natexlab{c}})Bhattacharya, Pati \&
  Yang}]{bhattacharya2016bayesian}
\textsc{Bhattacharya, A.}, \textsc{Pati, D.} \& \textsc{Yang, Y.}
  (2016{\natexlab{c}}).
\newblock Bayesian fractional posteriors.
\newblock \textit{arXiv preprint arXiv:1611.01125} .

\bibitem[{Bondell \& Reich(2012)}]{bondell2012consistent}
\textsc{Bondell, H.~D.} \& \textsc{Reich, B.~J.} (2012).
\newblock Consistent high-dimensional {B}ayesian variable selection via
  penalized credible regions.
\newblock \textit{Journal of the American Statistical Association}
  \textbf{107}, 1610--1624.

\bibitem[{Brown et~al.(1998)Brown, Vannucci \& Fearn}]{brown1998multivariate}
\textsc{Brown, P.~J.}, \textsc{Vannucci, M.} \& \textsc{Fearn, T.} (1998).
\newblock Multivariate {B}ayesian variable selection and prediction.
\newblock \textit{Journal of the Royal Statistical Society: Series B
  (Statistical Methodology)} \textbf{60}, 627--641.

\bibitem[{Bunea et~al.(2011)Bunea, She \& Wegkamp}]{bunea2011optimal}
\textsc{Bunea, F.}, \textsc{She, Y.} \& \textsc{Wegkamp, M.~H.} (2011).
\newblock Optimal selection of reduced rank estimators of high-dimensional
  matrices.
\newblock \textit{The Annals of Statistics} \textbf{39}, 1282--1309.

\bibitem[{Bunea et~al.(2012)Bunea, She, Wegkamp et~al.}]{bunea2012joint}
\textsc{Bunea, F.}, \textsc{She, Y.}, \textsc{Wegkamp, M.~H.} et~al. (2012).
\newblock Joint variable and rank selection for parsimonious estimation of
  high-dimensional matrices.
\newblock \textit{The Annals of Statistics} \textbf{40}, 2359--2388.

\bibitem[{Carvalho et~al.(2010)Carvalho, Polson \&
  Scott}]{carvalho2010horseshoe}
\textsc{Carvalho, C.}, \textsc{Polson, N.} \& \textsc{Scott, J.} (2010).
\newblock {The horseshoe estimator for sparse signals}.
\newblock \textit{Biometrika} \textbf{97}, 465--480.

\bibitem[{Chen et~al.(2013)Chen, Dong \& Chan}]{chen2013reduced}
\textsc{Chen, K.}, \textsc{Dong, H.} \& \textsc{Chan, K.-S.} (2013).
\newblock Reduced rank regression via adaptive nuclear norm penalization.
\newblock \textit{Biometrika} \textbf{100}, 901--920.

\bibitem[{Chen \& Huang(2012)}]{chen2012sparse}
\textsc{Chen, L.} \& \textsc{Huang, J.~Z.} (2012).
\newblock Sparse reduced-rank regression for simultaneous dimension reduction
  and variable selection.
\newblock \textit{Journal of the American Statistical Association}
  \textbf{107}, 1533--1545.

\bibitem[{Chun \& Kele{\c{s}}(2010)}]{chun2010sparse}
\textsc{Chun, H.} \& \textsc{Kele{\c{s}}, S.} (2010).
\newblock Sparse partial least squares regression for simultaneous dimension
  reduction and variable selection.
\newblock \textit{Journal of the Royal Statistical Society: Series B
  (Statistical Methodology)} \textbf{72}, 3--25.

\bibitem[{Friedman et~al.(2007)Friedman, Hastie, H{\"o}fling \&
  Tibshirani}]{friedman2007pathwise}
\textsc{Friedman, J.}, \textsc{Hastie, T.}, \textsc{H{\"o}fling, H.} \&
  \textsc{Tibshirani, R.} (2007).
\newblock Pathwise coordinate optimization.
\newblock \textit{The Annals of Applied Statistics} \textbf{1}, 302--332.

\bibitem[{George \& McCulloch(1993)}]{george1993variable}
\textsc{George, E.~I.} \& \textsc{McCulloch, R.~E.} (1993).
\newblock Variable selection via {G}ibbs sampling.
\newblock \textit{Journal of the American Statistical Association} \textbf{88},
  881--889.

\bibitem[{Geweke(1996)}]{geweke1996bayesian}
\textsc{Geweke, J.} (1996).
\newblock Bayesian reduced rank regression in econometrics.
\newblock \textit{Journal of econometrics} \textbf{75}, 121--146.

\bibitem[{Ghosal et~al.(2000)Ghosal, Ghosh \& van~der Vaart}]{ghosal2000}
\textsc{Ghosal, S.}, \textsc{Ghosh, J.~K.} \& \textsc{van~der Vaart, A.~W.}
  (2000).
\newblock Convergence rates of posterior distributions.
\newblock \textit{The Annals of Statistics} \textbf{28}, 500--531.

\bibitem[{Ghosh \& Chakrabarti(2017)}]{ghosh2017asymptotic}
\textsc{Ghosh, P.} \& \textsc{Chakrabarti, A.} (2017).
\newblock Asymptotic optimality of one-group shrinkage priors in sparse
  high-dimensional problems.
\newblock \textit{Bayesian Analysis} To appear.

\bibitem[{Golub \& van Loan(1996)}]{golub2012matrix}
\textsc{Golub, G.~H.} \& \textsc{van Loan, C.~F.} (1996).
\newblock \textit{Matrix Computations}.
\newblock John Hopkins University Press, 3rd ed.

\bibitem[{Green(1995)}]{green1995reversible}
\textsc{Green, P.~J.} (1995).
\newblock Reversible jump {M}arkov chain {M}onte {C}arlo computation and
  {B}ayesian model determination.
\newblock \textit{Biometrika} \textbf{82}, 711--732.

\bibitem[{Hahn \& Carvalho(2015)}]{hahn2015decoupling}
\textsc{Hahn, P.~R.} \& \textsc{Carvalho, C.~M.} (2015).
\newblock Decoupling shrinkage and selection in {B}ayesian linear models: a
  posterior summary perspective.
\newblock \textit{Journal of the American Statistical Association}
  \textbf{110}, 435--448.

\bibitem[{Hoff(2009)}]{hoff2007gibbs}
\textsc{Hoff, P.} (2009).
\newblock Simulation of the matrix {B}ingham-von {M}ises-{F}isher distribution,
  with applications to multivariate and relational data.
\newblock \textit{Journal of Computational and Graphical Statistics}
  \textbf{18}, 438--456.

\bibitem[{Izenman(1975)}]{izenman1975reduced}
\textsc{Izenman, A.~J.} (1975).
\newblock Reduced-rank regression for the multivariate linear model.
\newblock \textit{Journal of {M}ultivariate {A}nalysis} \textbf{5}, 248--264.

\bibitem[{Kundu et~al.(2013)Kundu, Baladandayuthapani \&
  Mallick}]{kundu2013bayes}
\textsc{Kundu, S.}, \textsc{Baladandayuthapani, V.} \& \textsc{Mallick, B.~K.}
  (2013).
\newblock Bayes regularized graphical model estimation in high dimensions.
\newblock \textit{arXiv preprint arXiv:1308.3915} .

\bibitem[{Lee et~al.(2002)Lee, Rinaldi, Robert, Odom, Bar-Joseph, Gerber,
  Hannett, Harbison, Thompson, Simon et~al.}]{lee2002transcriptional}
\textsc{Lee, T.~I.}, \textsc{Rinaldi, N.~J.}, \textsc{Robert, F.},
  \textsc{Odom, D.~T.}, \textsc{Bar-Joseph, Z.}, \textsc{Gerber, G.~K.},
  \textsc{Hannett, N.~M.}, \textsc{Harbison, C.~T.}, \textsc{Thompson, C.~M.},
  \textsc{Simon, I.} et~al. (2002).
\newblock Transcriptional regulatory networks in saccharomyces cerevisiae.
\newblock \textit{Science} \textbf{298}, 799--804.

\bibitem[{Lim \& Teh(2007)}]{lim2007variational}
\textsc{Lim, Y.~J.} \& \textsc{Teh, Y.~W.} (2007).
\newblock Variational {B}ayesian approach to movie rating prediction.
\newblock In \textit{Proceedings of KDD cup and workshop}, vol.~7. Citeseer.

\bibitem[{Liu \& Wu(1999)}]{liu1999parameter}
\textsc{Liu, J.~S.} \& \textsc{Wu, Y.~N.} (1999).
\newblock Parameter expansion for data augmentation.
\newblock \textit{Journal of the American Statistical Association} \textbf{94},
  1264--1274.

\bibitem[{Lopes \& West(2004)}]{lopes2004bayesian}
\textsc{Lopes, H.~F.} \& \textsc{West, M.} (2004).
\newblock Bayesian model assessment in factor analysis.
\newblock \textit{Statistica Sinica} , 41--67.

\bibitem[{Lucas et~al.(2006)Lucas, Carvalho, Wang, Bild, Nevins \&
  West}]{lucas2006sparse}
\textsc{Lucas, J.}, \textsc{Carvalho, C.}, \textsc{Wang, Q.}, \textsc{Bild,
  A.}, \textsc{Nevins, J.} \& \textsc{West, M.} (2006).
\newblock Sparse statistical modelling in gene expression genomics.
\newblock \textit{Bayesian Inference for Gene Expression and Proteomics}
  \textbf{1}, 0--1.

\bibitem[{Pati et~al.(2014)Pati, Bhattacharya, Pillai, Dunson
  et~al.}]{pati2014posterior}
\textsc{Pati, D.}, \textsc{Bhattacharya, A.}, \textsc{Pillai, N.~S.},
  \textsc{Dunson, D.} et~al. (2014).
\newblock Posterior contraction in sparse {B}ayesian factor models for massive
  covariance matrices.
\newblock \textit{The Annals of Statistics} \textbf{42}, 1102--1130.

\bibitem[{Polson \& Scott(2010)}]{polson2010shrink}
\textsc{Polson, N.~G.} \& \textsc{Scott, J.~G.} (2010).
\newblock Shrink globally, act locally: sparse {B}ayesian regularization and
  prediction.
\newblock \textit{Bayesian Statistics} \textbf{9}, 501--538.

\bibitem[{Polson et~al.(2014)Polson, Scott \& Windle}]{polson2014bayesian}
\textsc{Polson, N.~G.}, \textsc{Scott, J.~G.} \& \textsc{Windle, J.} (2014).
\newblock The {B}ayesian bridge.
\newblock \textit{Journal of the Royal Statistical Society: Series B
  (Statistical Methodology)} \textbf{76}, 713--733.

\bibitem[{Rue(2001)}]{rue2001fast}
\textsc{Rue, H.} (2001).
\newblock Fast sampling of gaussian markov random fields.
\newblock \textit{Journal of the Royal Statistical Society: Series B
  (Statistical Methodology)} \textbf{63}, 325--338.

\bibitem[{Ruffieux et~al.(2017)Ruffieux, Davison, Hager \&
  Irincheeva}]{ruffieux2017efficient}
\textsc{Ruffieux, H.}, \textsc{Davison, A.~C.}, \textsc{Hager, J.} \&
  \textsc{Irincheeva, I.} (2017).
\newblock Efficient inference for genetic association studies with multiple
  outcomes.
\newblock \textit{Biostatistics} , To appear.

\bibitem[{Salakhutdinov \& Mnih(2008)}]{salakhutdinov2008bayesian}
\textsc{Salakhutdinov, R.} \& \textsc{Mnih, A.} (2008).
\newblock Bayesian probabilistic matrix factorization using markov chain monte
  carlo.
\newblock In \textit{Proceedings of the 25th international conference on
  Machine learning}. ACM.

\bibitem[{Scott et~al.(2010)Scott, Berger et~al.}]{scott2010bayes}
\textsc{Scott, J.~G.}, \textsc{Berger, J.~O.} et~al. (2010).
\newblock Bayes and empirical-{B}ayes multiplicity adjustment in the
  variable-selection problem.
\newblock \textit{The Annals of Statistics} \textbf{38}, 2587--2619.

\bibitem[{Spellman et~al.(1998)Spellman, Sherlock, Zhang, Iyer, Anders, Eisen,
  Brown, Botstein \& Futcher}]{spellman1998comprehensive}
\textsc{Spellman, P.~T.}, \textsc{Sherlock, G.}, \textsc{Zhang, M.~Q.},
  \textsc{Iyer, V.~R.}, \textsc{Anders, K.}, \textsc{Eisen, M.~B.},
  \textsc{Brown, P.~O.}, \textsc{Botstein, D.} \& \textsc{Futcher, B.} (1998).
\newblock Comprehensive identification of cell cycle--regulated genes of the
  yeast saccharomyces cerevisiae by microarray hybridization.
\newblock \textit{Molecular biology of the cell} \textbf{9}, 3273--3297.

\bibitem[{van~der Pas et~al.(2014)van~der Pas, Kleijn \& van~der
  Vaart}]{van2014horseshoe}
\textsc{van~der Pas, S.}, \textsc{Kleijn, B.} \& \textsc{van~der Vaart, A.}
  (2014).
\newblock The horseshoe estimator: Posterior concentration around nearly black
  vectors.
\newblock \textit{Electronic Journal of Statistics} \textbf{8}, 2585--2618.

\bibitem[{Velu \& Reinsel(2013)}]{velu2013multivariate}
\textsc{Velu, R.} \& \textsc{Reinsel, G.~C.} (2013).
\newblock \textit{Multivariate reduced-rank regression: theory and
  applications}, vol. 136.
\newblock Springer Science \& Business Media.

\bibitem[{Vershynin(2010)}]{vershynin2010introduction}
\textsc{Vershynin, R.} (2010).
\newblock Introduction to the non-asymptotic analysis of random matrices.
\newblock \textit{arXiv preprint arXiv:1011.3027} .

\bibitem[{Wang(2010)}]{wang2010sparse}
\textsc{Wang, H.} (2010).
\newblock Sparse seemingly unrelated regression modelling: Applications in
  finance and econometrics.
\newblock \textit{Computational Statistics \& Data Analysis} \textbf{54},
  2866--2877.

\bibitem[{Wang et~al.(2007)Wang, Chen \& Li}]{wang2007group}
\textsc{Wang, L.}, \textsc{Chen, G.} \& \textsc{Li, H.} (2007).
\newblock Group scad regression analysis for microarray time course gene
  expression data.
\newblock \textit{Bioinformatics} \textbf{23}, 1486--1494.

\bibitem[{Yuan et~al.(2007)Yuan, Ekici, Lu \& Monteiro}]{yuan2007dimension}
\textsc{Yuan, M.}, \textsc{Ekici, A.}, \textsc{Lu, Z.} \& \textsc{Monteiro, R.}
  (2007).
\newblock Dimension reduction and coefficient estimation in multivariate linear
  regression.
\newblock \textit{Journal of the Royal Statistical Society: Series B
  (Statistical Methodology)} \textbf{69}, 329--346.

\bibitem[{Yuan \& Lin(2006)}]{yuan2006model}
\textsc{Yuan, M.} \& \textsc{Lin, Y.} (2006).
\newblock Model selection and estimation in regression with grouped variables.
\newblock \textit{Journal of the Royal Statistical Society: Series B
  (Statistical Methodology)} \textbf{68}, 49--67.

\end{thebibliography}
