Optimal High-Dimensional Shrinkage Covariance Estimation for Elliptical Distributions
\IEEEauthorblockN{Esa Ollila
\IEEEauthorblockN{Michael Shell\IEEEauthorrefmark{1
We derive an optimal shrinkage sample covariance matrix (SCM) estimator which is suitable for high dimensional problems %(large dimension problems) and when sampling from an unspecified elliptically symmetric distribution. Specifically, we derive the optimal (oracle) shrinkage parameters that obtain the minimum mean-squared error (MMSE) between the shrinkage SCM and the true covariance matrix when sampling from an elliptical distribution. Subsequently, we show how the oracle shrinkage parameters can be consistently estimated under the random matrix theory regime. Simulations show the advantage of the proposed estimator over the conventional shrinkage SCM estimator due to Ledoit and Wolf (2004). The proposed shrinkage SCM estimator often provides significantly better performance than the Ledoit-Wolf estimator and has the advantage that consistency is guaranteed over the whole class of elliptical distributions with finite 4th order moments.
definitions_esa.tex
%%%%% ESA's DEFINITIONS start here  %%%%%%%%%%

\newcommand{\pdim}{p}
\newcommand{\ndim}{n}
\newcommand{\paino}{\sl} 
\newcommand{\hop}{\top}              % Hermitian transpose 
\newcommand{\covm}{\bom \Sigma} %\bom{\mathcal C}} % Covariance matrix 
\newcommand{\M}{\covm} %\bom{\mathcal C}} % Covariance matrix 

\newcommand{\E}{\mathbb{E}}                  
\newcommand{\bo}[1]{\mathbf{#1}}              % boldface math 
\newcommand{\bom}[1]{\bm{#1}}%    {\boldsymbol{#1}}    % boldface math (for greek letters)
\newcommand{\RE}{\mathrm{RE}}
\newcommand{\CE}{\mathrm{CE}}
\newcommand{\CN}{\C \mathcal N}
\renewcommand{\bar}{\overline}
\newcommand{\be}{\beta}
\newcommand{\pr}{\partial}
\newcommand{\iidsim}{\overset{iid}{\sim}}

\newcommand{\al}{\alpha}  
\newcommand{\Id}{\bo I}  
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\bmat}{\begin{pmatrix}}
\newcommand{\emat}{\end{pmatrix}}
\newcommand{\RR}{{\mathbb{R}}}

%%%%% FRED's DEFINITIONS start here%%%

%\newcommand{\sq}[1]{\ensuremath{\mathsf{#1}}}
\newcommand{\pw}[1]{\psframebox[linewidth=0.4pt,framearc=.3]{#1}}
\def\Gr#1{\mathbf #1}
\def\tred#1{{\color{red}#1}}

% Shorthand notations for sets N, Z, R and C
\newcommand{\ZZ}{{\mathbb{Z}}}
\newcommand{\CC}{{\mathbb{C}}}
\newcommand{\cD}{{\mathcal{D}}}
\newcommand{\herm}{{\sf H}}
\newcommand{\trans}{{\sf T}}
\newcommand{\uR}{\underline\R}
\renewcommand{\P}{{\bf P}}
\newcommand{\N}{{\bf N}}
\newcommand{\V}{{\bf V}}
\newcommand{\J}{{\bf J}}
\renewcommand{\S}{{\bf S}}
\newcommand{\ASV}{\mathrm{ASV}}
\newcommand{\sca}{\sigma}

\newcommand{\I}{\bo I} 
\newcommand{\mD}{{\bf D}}
\newcommand{\im}{j} %\jmath}  % imaginary unit      

\renewcommand{\dim}{\pdim} 

\newcommand{\wei}{u} 
\newcommand{\R}{\mathbb{R}}
\newcommand{\F}{\mathbb{F}}
\newcommand{\X}{{\bf X}}
\newcommand{\mC}{{\bf C}}

\newcommand{\Z}{{\bf Z}}
\newcommand{\G}{{\bf G}}
\newcommand{\Q}{{\bf Q}}
\newcommand{\Y}{{\bf Y}}
\newcommand{\W}{{\bf W}}
\newcommand{\A}{{\bf A}}
\newcommand{\B}{{\bf B}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\T}{{\bf T}}
\renewcommand{\L}{{\bf \Lambda}}
\newcommand{\U}{{\bf U}}
\newcommand{\K}{{\bf K}}
\newcommand{\Mn}{\M}  
\newcommand{\commat}{\bo K_\pdim}  % COMMUTATION matrix 
\newcommand{\ve}{\mathrm{vec}} % vectorize a matrix 

\renewcommand{\H}{{\bf H}}
%\newcommand{\I}{{\bf I}}
\newcommand{\x}{{\bf x}}
\renewcommand{\u}{{\bf u}}

\newcommand{\h}{{\bf h}}
\newcommand{\s}{{\bf s}}
\renewcommand{\a}{{\bf a}}
\newcommand{\z}{{\bf x}}
\newcommand{\p}{{\bf p}}
\newcommand{\Mal}{\bo C_{\be}}  
\newcommand{\blambda}{{\bm \lambda}}
\newcommand{\bnu}{{\bm \nu}}
\renewcommand{\b}{{\bf b}}
\newcommand{\y}{{\bf y}}
\renewcommand{\v}{{\bf v}}
\newcommand{\uu}{{\bf u}}
\newcommand{\e}{{\bf e}}
\newcommand{\w}{{\bf w}}
\newcommand{\f}{{\bf f}}
\newcommand{\Fr}{\mathrm{F}}
\newcommand{\Mor}{\bo S}
\newcommand{\Mcl}{\bo S_{\al,\be}}
\newcommand{\MSE}{\mathrm{MSE}}
\newcommand{\PDH}{\mathcal{S}}   % Class of PDH matrices
\newcommand{\qva}{t}

\newcommand{\ka}{\kappa}

\newcommand{\n}{{\bf n}}
\newcommand{\oh}{{\frac{1}{2}}}
\newcommand{\qdim}{Q}
\newcommand{\bGamma}{\bar{\Gamma}}
\newcommand{\um}{\underline{m}}
\newcommand{\uY}{\underline{\Y}}
\newcommand{\uF}{\underline{F}}
\newcommand{\hum}{\underline{\hat{m}}}
\newcommand{\tm}{\tilde{m}}
\newcommand{\htm}{\hat{\tilde{m}}}
\newcommand{\asto}{\overset{\rm a.s.}{\longrightarrow}}
\newcommand{\eqdef}{\overset{\Delta}{=}}
\newcommand{\EE}{{\rm E}}
\newcommand{\kdim}{K}

\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator{\Tr}{tr}
\DeclareMathOperator{\diag}{diag}
\DeclareMathOperator{\cov}{cov} 



\newcommand{\dI}{d_{{\rm KL}}}
\newcommand{\dE}{d_{{\rm E}}}
\newcommand{\dS}{d_{{\rm S}}}
\newcommand{\dR}{d_{{\rm R}}}
\newcommand{\dX}{d_{{\rm X}}}
\newcommand{\dF}{d_{{\rm F}}}


rscmEUSIPCO17final.tex

%% bare_conf.tex
%% V1.4b
%% 2015/08/26
%% by Michael Shell
%% See:
%% http://www.michaelshell.org/
%% for current contact information.
%%
%% This is a skeleton file demonstrating the use of IEEEtran.cls
%% (requires IEEEtran.cls version 1.8b or later) with an IEEE
%% conference paper.
%%
%% Support sites:
%% http://www.michaelshell.org/tex/ieeetran/
%% http://www.ctan.org/pkg/ieeetran
%% and
%% http://www.ieee.org/

%%*************************************************************************
%% Legal Notice:
%% This code is offered as-is without any warranty either expressed or
%% implied; without even the implied warranty of MERCHANTABILITY or
%% FITNESS FOR A PARTICULAR PURPOSE! 
%% User assumes all risk.
%% In no event shall the IEEE or any contributor to this code be liable for
%% any damages or losses, including, but not limited to, incidental,
%% consequential, or any other damages, resulting from the use or misuse
%% of any information contained here.
%%
%% All comments are the opinions of their respective authors and are not
%% necessarily endorsed by the IEEE.
%%
%% This work is distributed under the LaTeX Project Public License (LPPL)
%% ( http://www.latex-project.org/ ) version 1.3, and may be freely used,
%% distributed and modified. A copy of the LPPL, version 1.3, is included
%% in the base LaTeX documentation of all distributions of LaTeX released
%% 2003/12/01 or later.
%% Retain all contribution notices and credits.
%% ** Modified files should be clearly indicated as such, including  **
%% ** renaming them and changing author support contact information. **
%%*************************************************************************


% *** Authors should verify (and, if needed, correct) their LaTeX system  ***
% *** with the testflow diagnostic prior to trusting their LaTeX platform ***
% *** with production work. The IEEE's font choices and paper sizes can   ***
% *** trigger bugs that do not appear when using other class files.       ***                          ***
% The testflow support page is at:
% http://www.michaelshell.org/tex/testflow/



\documentclass[conference]{IEEEtran}
% Some Computer Society conferences also require the compsoc mode option,
% but others use the standard conference format.
%
% If IEEEtran.cls has not been installed into the LaTeX system files,
% manually specify the path to it like:
% \documentclass[conference]{../sty/IEEEtran}


\usepackage{amsmath,amsfonts,graphicx,bm,amssymb,amsthm} 
%enumerate,epsfig,psfrag}

\usepackage{amsopn}
\usepackage{subfigure}

\include{definitions_esa}
\newcommand{\bmu}{\bom \mu}

\newcounter{ctheorem}
\newtheorem{theorem}[ctheorem]{Theorem}
\newcounter{cremark}
\newtheorem{remark}[cremark]{Remark}
\newcounter{clemma}
\newtheorem{lemma}[clemma]{Lemma}
\newcounter{ccorollary}
\newtheorem{corollary}[ccorollary]{Corollary}


% Some very useful LaTeX packages include:
% (uncomment the ones you want to load)


% *** MISC UTILITY PACKAGES ***
%
%\usepackage{ifpdf}
% Heiko Oberdiek's ifpdf.sty is very useful if you need conditional
% compilation based on whether the output is pdf or dvi.
% usage:
% \ifpdf
%   % pdf code
% \else
%   % dvi code
% \fi
% The latest version of ifpdf.sty can be obtained from:
% http://www.ctan.org/pkg/ifpdf
% Also, note that IEEEtran.cls V1.7 and later provides a builtin
% \ifCLASSINFOpdf conditional that works the same way.
% When switching from latex to pdflatex and vice-versa, the compiler may
% have to be run twice to clear warning/error messages.






% *** CITATION PACKAGES ***
%
%\usepackage{cite}
% cite.sty was written by Donald Arseneau
% V1.6 and later of IEEEtran pre-defines the format of the cite.sty package
% \cite{} output to follow that of the IEEE. Loading the cite package will
% result in citation numbers being automatically sorted and properly
% "compressed/ranged". e.g., [1], [9], [2], [7], [5], [6] without using
% cite.sty will become [1], [2], [5]--[7], [9] using cite.sty. cite.sty's
% \cite will automatically add leading space, if needed. Use cite.sty's
% noadjust option (cite.sty V3.8 and later) if you want to turn this off
% such as if a citation ever needs to be enclosed in parenthesis.
% cite.sty is already installed on most LaTeX systems. Be sure and use
% version 5.0 (2009-03-20) and later if using hyperref.sty.
% The latest version can be obtained at:
% http://www.ctan.org/pkg/cite
% The documentation is contained in the cite.sty file itself.






% *** GRAPHICS RELATED PACKAGES ***
%
\ifCLASSINFOpdf
  % \usepackage[pdftex]{graphicx}
  % declare the path(s) where your graphic files are
%   \graphicspath{{./pics/}}
  % and their extensions so you won't have to specify these with
  % every instance of \includegraphics
  % \DeclareGraphicsExtensions{.pdf,.jpeg,.png}
\else
  % or other class option (dvipsone, dvipdf, if not using dvips). graphicx
  % will default to the driver specified in the system graphics.cfg if no
  % driver is specified.
  % \usepackage[dvips]{graphicx}
  % declare the path(s) where your graphic files are
  % \graphicspath{{../eps/}}
  % and their extensions so you won't have to specify these with
  % every instance of \includegraphics
  % \DeclareGraphicsExtensions{.eps}
\fi
% graphicx was written by David Carlisle and Sebastian Rahtz. It is
% required if you want graphics, photos, etc. graphicx.sty is already
% installed on most LaTeX systems. The latest version and documentation
% can be obtained at: 
% http://www.ctan.org/pkg/graphicx
% Another good source of documentation is "Using Imported Graphics in
% LaTeX2e" by Keith Reckdahl which can be found at:
% http://www.ctan.org/pkg/epslatex
%
% latex, and pdflatex in dvi mode, support graphics in encapsulated
% postscript (.eps) format. pdflatex in pdf mode supports graphics
% in .pdf, .jpeg, .png and .mps (metapost) formats. Users should ensure
% that all non-photo figures use a vector format (.eps, .pdf, .mps) and
% not a bitmapped formats (.jpeg, .png). The IEEE frowns on bitmapped formats
% which can result in "jaggedy"/blurry rendering of lines and letters as
% well as large increases in file sizes.
%
% You can find documentation about the pdfTeX application at:
% http://www.tug.org/applications/pdftex





% *** MATH PACKAGES ***
%
%\usepackage{amsmath}
% A popular package from the American Mathematical Society that provides
% many useful and powerful commands for dealing with mathematics.
%
% Note that the amsmath package sets \interdisplaylinepenalty to 10000
% thus preventing page breaks from occurring within multiline equations. Use:
%\interdisplaylinepenalty=2500
% after loading amsmath to restore such page breaks as IEEEtran.cls normally
% does. amsmath.sty is already installed on most LaTeX systems. The latest
% version and documentation can be obtained at:
% http://www.ctan.org/pkg/amsmath





% *** SPECIALIZED LIST PACKAGES ***
%
%\usepackage{algorithmic}
% algorithmic.sty was written by Peter Williams and Rogerio Brito.
% This package provides an algorithmic environment fo describing algorithms.
% You can use the algorithmic environment in-text or within a figure
% environment to provide for a floating algorithm. Do NOT use the algorithm
% floating environment provided by algorithm.sty (by the same authors) or
% algorithm2e.sty (by Christophe Fiorio) as the IEEE does not use dedicated
% algorithm float types and packages that provide these will not provide
% correct IEEE style captions. The latest version and documentation of
% algorithmic.sty can be obtained at:
% http://www.ctan.org/pkg/algorithms
% Also of interest may be the (relatively newer and more customizable)
% algorithmicx.sty package by Szasz Janos:
% http://www.ctan.org/pkg/algorithmicx




% *** ALIGNMENT PACKAGES ***
%
%\usepackage{array}
% Frank Mittelbach's and David Carlisle's array.sty patches and improves
% the standard LaTeX2e array and tabular environments to provide better
% appearance and additional user controls. As the default LaTeX2e table
% generation code is lacking to the point of almost being broken with
% respect to the quality of the end results, all users are strongly
% advised to use an enhanced (at the very least that provided by array.sty)
% set of table tools. array.sty is already installed on most systems. The
% latest version and documentation can be obtained at:
% http://www.ctan.org/pkg/array


% IEEEtran contains the IEEEeqnarray family of commands that can be used to
% generate multiline equations as well as matrices, tables, etc., of high
% quality.




% *** SUBFIGURE PACKAGES ***
%\ifCLASSOPTIONcompsoc
%  \usepackage[caption=false,font=normalsize,labelfont=sf,textfont=sf]{subfig}
%\else
%  \usepackage[caption=false,font=footnotesize]{subfig}
%\fi
% subfig.sty, written by Steven Douglas Cochran, is the modern replacement
% for subfigure.sty, the latter of which is no longer maintained and is
% incompatible with some LaTeX packages including fixltx2e. However,
% subfig.sty requires and automatically loads Axel Sommerfeldt's caption.sty
% which will override IEEEtran.cls' handling of captions and this will result
% in non-IEEE style figure/table captions. To prevent this problem, be sure
% and invoke subfig.sty's "caption=false" package option (available since
% subfig.sty version 1.3, 2005/06/28) as this is will preserve IEEEtran.cls
% handling of captions.
% Note that the Computer Society format requires a larger sans serif font
% than the serif footnote size font used in traditional IEEE formatting
% and thus the need to invoke different subfig.sty package options depending
% on whether compsoc mode has been enabled.
%
% The latest version and documentation of subfig.sty can be obtained at:
% http://www.ctan.org/pkg/subfig




% *** FLOAT PACKAGES ***
%
%\usepackage{fixltx2e}
% fixltx2e, the successor to the earlier fix2col.sty, was written by
% Frank Mittelbach and David Carlisle. This package corrects a few problems
% in the LaTeX2e kernel, the most notable of which is that in current
% LaTeX2e releases, the ordering of single and double column floats is not
% guaranteed to be preserved. Thus, an unpatched LaTeX2e can allow a
% single column figure to be placed prior to an earlier double column
% figure.
% Be aware that LaTeX2e kernels dated 2015 and later have fixltx2e.sty's
% corrections already built into the system in which case a warning will
% be issued if an attempt is made to load fixltx2e.sty as it is no longer
% needed.
% The latest version and documentation can be found at:
% http://www.ctan.org/pkg/fixltx2e


%\usepackage{stfloats}
% stfloats.sty was written by Sigitas Tolusis. This package gives LaTeX2e
% the ability to do double column floats at the bottom of the page as well
% as the top. (e.g., "\begin{figure*}[!b]" is not normally possible in
% LaTeX2e). It also provides a command:
%\fnbelowfloat
% to enable the placement of footnotes below bottom floats (the standard
% LaTeX2e kernel puts them above bottom floats). This is an invasive package
% which rewrites many portions of the LaTeX2e float routines. It may not work
% with other packages that modify the LaTeX2e float routines. The latest
% version and documentation can be obtained at:
% http://www.ctan.org/pkg/stfloats
% Do not use the stfloats baselinefloat ability as the IEEE does not allow
% \baselineskip to stretch. Authors submitting work to the IEEE should note
% that the IEEE rarely uses double column equations and that authors should try
% to avoid such use. Do not be tempted to use the cuted.sty or midfloat.sty
% packages (also by Sigitas Tolusis) as the IEEE does not format its papers in
% such ways.
% Do not attempt to use stfloats with fixltx2e as they are incompatible.
% Instead, use Morten Hogholm'a dblfloatfix which combines the features
% of both fixltx2e and stfloats:
%
% \usepackage{dblfloatfix}
% The latest version can be found at:
% http://www.ctan.org/pkg/dblfloatfix




% *** PDF, URL AND HYPERLINK PACKAGES ***
%
%\usepackage{url}
% url.sty was written by Donald Arseneau. It provides better support for
% handling and breaking URLs. url.sty is already installed on most LaTeX
% systems. The latest version and documentation can be obtained at:
% http://www.ctan.org/pkg/url
% Basically, \url{my_url_here}.




% *** Do not adjust lengths that control margins, column widths, etc. ***
% *** Do not use packages that alter fonts (such as pslatex).         ***
% There should be no need to do such things with IEEEtran.cls V1.6 and later.
% (Unless specifically asked to do so by the journal or conference you plan
% to submit to, of course. )


% correct bad hyphenation here
\hyphenation{op-tical net-works semi-conduc-tor}


\begin{document}
%
% paper title
% Titles are generally capitalized except for words such as a, an, and, as,
% at, but, by, for, in, nor, of, on, or, the, to and up, which are usually
% not capitalized unless they are the first or last word of the title.
% Linebreaks \\ can be used within to get better formatting as desired.
% Do not put math or special symbols in the title.
\title{Optimal High-Dimensional Shrinkage Covariance Estimation for Elliptical Distributions}


% author names and affiliations
% use a multiple column layout for up to three different
% affiliations
\author{\IEEEauthorblockN{Esa Ollila}
\IEEEauthorblockA{Department of Signal Processing and Acoustics\\
Aalto University, Finland
%Email: http://www.michaelshell.org/contact.html}
%\and
%\IEEEauthorblockN{Homer Simpson}
%\IEEEauthorblockA{Twentieth Century Fox\\
%Springfield, USA\\
%Email: homer@thesimpsons.com}
%\and
%\IEEEauthorblockN{James Kirk\\ and Montgomery Scott}
%\IEEEauthorblockA{Starfleet Academy\\
%San Francisco, California 96678--2391\\
%Telephone: (800) 555--1212\\
%Fax: (888) 555--1212
}}

% conference papers do not typically use \thanks and this command
% is locked out in conference mode. If really needed, such as for
% the acknowledgment of grants, issue a \IEEEoverridecommandlockouts
% after \documentclass

% for over three affiliations, or if they all won't fit within the width
% of the page, use this alternative format:
% 
%\author{\IEEEauthorblockN{Michael Shell\IEEEauthorrefmark{1},
%Homer Simpson\IEEEauthorrefmark{2},
%James Kirk\IEEEauthorrefmark{3}, 
%Montgomery Scott\IEEEauthorrefmark{3} and
%Eldon Tyrell\IEEEauthorrefmark{4}}
%\IEEEauthorblockA{\IEEEauthorrefmark{1}School of Electrical and Computer Engineering\\
%Georgia Institute of Technology,
%Atlanta, Georgia 30332--0250\\ Email: see http://www.michaelshell.org/contact.html}
%\IEEEauthorblockA{\IEEEauthorrefmark{2}Twentieth Century Fox, Springfield, USA\\
%Email: homer@thesimpsons.com}
%\IEEEauthorblockA{\IEEEauthorrefmark{3}Starfleet Academy, San Francisco, California 96678-2391\\
%Telephone: (800) 555--1212, Fax: (888) 555--1212}
%\IEEEauthorblockA{\IEEEauthorrefmark{4}Tyrell Inc., 123 Replicant Street, Los Angeles, California 90210--4321}}




% use for special paper notices
%\IEEEspecialpapernotice{(Invited Paper)}


% make the title area
\maketitle

% As a general rule, do not put math, special symbols or citations
% in the abstract
\begin{abstract}
We derive an optimal shrinkage sample covariance matrix (SCM) estimator 
which is suitable for high dimensional problems %(large dimension problems)  
and when sampling from an unspecified elliptically symmetric  distribution. 
Specifically, we derive the optimal (oracle) shrinkage parameters that obtain the minimum mean-squared error (MMSE) between the shrinkage SCM and the true covariance matrix when sampling from an elliptical distribution.  Subsequently, we show how the oracle shrinkage parameters can be 
consistently estimated under the random matrix theory regime. 
Simulations show the advantage of the proposed estimator over the conventional shrinkage SCM estimator due to Ledoit and Wolf (2004). The proposed shrinkage SCM estimator often provides significantly better performance than the Ledoit-Wolf estimator and has the advantage that consistency is guaranteed over the whole class of elliptical distributions with finite 4th order moments. 
\end{abstract}

% no keywords




% For peer review papers, you can put extra information on the cover
% page as needed:
% \ifCLASSOPTIONpeerreview
% \begin{center} \bfseries EDICS Category: 3-BBND \end{center}
% \fi
%
% For peerreview papers, this IEEEtran command inserts a page break and
% creates the second title. It will be ignored for other modes.
\IEEEpeerreviewmaketitle



\section{Introduction}\label{sec:intro}

We consider the problem of estimating the covariance matrix based on a sample  $\x_1,\ldots, \x_n$ of independent and identically distributed (i.i.d.)  random vectors from an unspecified $\pdim$-variate distribution $\x \sim F$ with mean vector  $\E[\x] = \bo 0$  and $\pdim \times \pdim$ positive definite covariance matrix $\M= \E[\x \x^\top]$.  
%The problem is fundamental in signal processing and data mining as many classical  techniques require an estimate of the covariance matrix or some nonlinear function of it, e.g., the inverse covariance marix or its eigenvectors/eigenvalues. 
The  {\paino sample covariance matrix} (SCM) $ \S=\frac{1}{\ndim}\sum_{i=1}^\ndim \z_i \z_i^\top $ 
is the most commonly used estimator of the covariance matrix,    and when random sampling from a multivariate Gaussian 
$\mathcal N_\pdim(\bo 0,\covm)$ distribution, it is also the optimal maximum likelihood estimator (MLE).  
Estimation of high-dimensional (HD) covariance matrix when the sample size $\ndim$ is smaller, or not much larger than the dimension $\pdim$, has attracted a significant research interest in recent years. Indeed since such data  problems are becoming increasingly common in finance \cite{ledoit_wolf:2004}, genomics  or classification, %\cite{friedman1989regularized}, 
for example.  
 % The density generator is $g_\nu(t)  =  (1+t/(\nu-2))^{-(\pdim+\nu)/2}$ and we denote this case by $\x \sim t_\nu(\bo 0 , \M)$. 
Insufficient number of samples causes significant estimation errors in 
%many conventional estimators of the covariance matrix including 
 the SCM. Moreover, if  $p > n$, the SCM  $\S$ is always singular, i.e., not invertible even if the true covariance matrix $\M$ is known to be positive definite and hence non-singular. 
The commonly used approach %for improving the estimation accuracy of high-dimension covaraince matrix based on insufficient sample size 
is then to  use shrinkage regularization as in   \cite{ledoit_wolf:2004,ollila2014regularized,zhang2016automatic,chen2010shrinkage,couillet_mckay:2014,pascal2014generalized}, for example.  

One of the most ommonly used estimator in "large $\pdim$ compared to sample size $\ndim$ problems" is the {\paino regularized SCM (RSCM)}, 
\beq \label{eq:regSCM} 
\S_{\al, \be} = \be \S + \al    \I,   
\eeq 
where %$\bo S=\frac{ 1}{\ndim } \sum_{i=1}^\ndim \x_i \x_i ^\top$ denotes  the SCM and 
$\al, \be >0$ denotes the {\paino shrinkage (regularization) parameters}.   Optimal RSCM estimator is often defined as one that is based on  {\paino oracle 
shrinkage parameters} that minimize the mean squared error (MSE),   
\beq \label{eq:RegSCM_oracle} 
(\al_o,\be_o)=   \underset{\al,\be >0} {\arg \min} \Big \{ \mathrm{MSE}(\S_{\al,\be}) =  \E \Big[ \big\| \M - \S_{\al,\be}  \big\|^2_{\rm F}  \Big] \Big\} ,   
\eeq 
where $\| \cdot \|_\Fr$ denotes the Frobenius matrix norm ($\| \A \|_\Fr^2 = \tr(\A^\top \A)=\tr(\A \A^\top)$ for any matrix $\A$).  
The solution  $(\al_o,\be_o)$ are called "oracle" shrinkage parameters as they will obviously  depend on the true unknown covariance matrix $\Mn$ and hence can not be used in practise. The widely popular Ledoit-Wolf (LW-)RSCM  \cite{ledoit_wolf:2004}   is based on  consistent estimators  $(\hat \al_o^{{\tiny \mbox{LW}}}, \hat \be_o^{{\tiny \mbox{LW}}})$ of $(\al_o,\be_o)$ under the random matrix theory (RMT) regime. However, more accurate finite sample estimation performance %One can impose also more stringent  assumption such as 
can be obtained by assuming that the observations are from a specific $\pdim$-variate distribution, e.g., the multivariate normal distribution,  as has been shown in \cite{chen2010shrinkage}.  In this paper, we derive  consistent estimators of the oracle shrinkage parameters $( \al_o,\be_o)$ under the RMT regime when sampling from an {\it unspecified} elliptically symmetric distribution.  Elliptical distributions (see  \cite{muirhead:1982,fang_etal:1990,ollila2012complex}) constitute a large class of distributions that include e.g., the multivariate normal distribution, generalized Gaussian and all compound Gaussian distributions as special cases. 

%We are interested in obtaining  optimal shrinkage estimators of high-dimensional covariance matrix under 
The  RMT regime refers to the case that 
\begin{itemize} 
\item[(R1)] %\beq \label{eq:RMTregime}
$\ndim, \pdim \to \infty$ and $\pdim/\ndim \to c$, where $0 < c < \infty$.  
\end{itemize} 
% \quad \mbox{and} \quad  \pdim/\ndim \to c, \quad 0 < c < \infty.   %\tag{R1}
%\eeq 
Furthermore, we assume that the set of eigenvalues of $\Mn$ converge to a fixed spectrum, and that 
%\beq \label{eq:RMTregime2} 
\begin{itemize} 
\item[(R2)] As $\pdim \to \infty$,   $\eta_i  = \tr(\Mn^i)/\pdim \to  \eta_i^o$,   $ 0 < \eta_i^0< \infty$ for $i=1, \ldots,4$ %\tag{R2}
\end{itemize} 
%\eeq 
 Our numerical examples illustrate that  the RSCM estimator that is based on  the proposed consistent estimators $(\hat \al^{Ell}_o,\hat \be^{Ell}_o)$   
%perform equally well as the optimal Gaussian RSCM estimators when the data follows the Gaussian distributions but 
 outperform its competitors, e.g., the LW-RSCM estimator, when sampling from an elliptical population. 
  %them otherwise when the data is from a non-Gaussian (but elliptically symmetric) distribution. Furthermore, our proposed estimator outperform the popular  Ledoit-Wolf estimator in all considered examples. 

The paper is organized as follows.  In Section~\ref{sec:optim1} and Section~\ref{sec:optim2} we derive the optimal shrinkage parameters $(\al_o,\be_o)$ under the general assumption of sampling from any general   $\pdim$-variate distribution and  an  elliptical distribution with finite 4th order moments, respectively.  In Section~\ref{sec:estim}, consistent estimators  of $(\al_o,\be_o)$  are proposed under assumptions (R1) and (R2) when sampling from an unspecified elliptical distribution. Simulation studies of Section~\ref{sec:simul} illustrate that the proposed shrinkage estimator always outperforms the LW estimator when the samples are drawn from an elliptical population. 

%\subsection{Notation}
{\it Notation:} Let $\PDH_p$ be the open cone of positive definite $p \times p$ symmetric matrices, and let $\I$ be the identity matrix of proper dimension,  $\ve(\cdot)$ denotes an operator that transforms a matrix into a vector by stacking the columns of the matrix, $\tr(\cdot)$ denotes the matrix trace operator,   and  $\otimes$ denote the {\paino Kronecker product}: 
for any matrix $\A$ and $\B$, $\A \otimes \B$ is a block matrix with $(i,j)$-block being equal to $a_{ij}\B$. A {\paino commutation matrix} $\commat$ is a $\dim^2 \times \dim^2$ block matrix with $(i,j)$-block 
equal to a $\dim \times \dim$ matrix that has a $1$ at entry $(j,i)$ and $0$'s elsewhere. 
It has the following important property \cite{magnus_neudecker:1999}: $\commat \ve(\A)= \ve(\A^\top) $ for any $\pdim \times \pdim$ matrix $\A$.  

%We use notation $\E_{\bo 0,\M} [ F(\x) ]$ to denote that expected value of some function $F(\x)$ of a random vector  $\x \sim \mathcal E_\pdim(\bo 0 ,\M,g)$.  



% An example of a floating figure using the graphicx package.
% Note that \label must occur AFTER (or within) \caption.
% For figures, \caption should occur after the \includegraphics.
% Note that IEEEtran v1.7 and later has special internal code that
% is designed to preserve the operation of \label within \caption
% even when the captionsoff option is in effect. However, because
% of issues like this, it may be the safest practice to put all your
% \label just after \caption rather than within \caption{}.
%
% Reminder: the "draftcls" or "draftclsnofoot", not "draft", class
% option should be used if it is desired that the figures are to be
% displayed while in draft mode.
%
%\begin{figure}[!t]
%\centering
%\includegraphics[width=2.5in]{myfigure}
% where an .eps filename suffix will be assumed under latex, 
% and a .pdf suffix will be assumed for pdflatex; or what has been declared
% via \DeclareGraphicsExtensions.
%\caption{Simulation results for the network.}
%\label{fig_sim}
%\end{figure}

% Note that the IEEE typically puts floats only at the top, even when this
% results in a large percentage of a column being occupied by floats.


% An example of a double column floating figure using two subfigures.
% (The subfig.sty package must be loaded for this to work.)
% The subfigure \label commands are set within each subfloat command,
% and the \label for the overall figure must come after \caption.
% \hfil is used as a separator to get equal spacing.
% Watch out that the combined width of all the subfigures on a 
% line do not exceed the text width or a line break will occur.
%
%\begin{figure*}[!t]
%\centering
%\subfloat[Case I]{\includegraphics[width=2.5in]{box}%
%\label{fig_first_case}}
%\hfil
%\subfloat[Case II]{\includegraphics[width=2.5in]{box}%
%\label{fig_second_case}}
%\caption{Simulation results for the network.}
%\label{fig_sim}
%\end{figure*}
%
% Note that often IEEE papers with subfigures do not employ subfigure
% captions (using the optional argument to \subfloat[]), but instead will
% reference/describe all of them (a), (b), etc., within the main caption.
% Be aware that for subfig.sty to generate the (a), (b), etc., subfigure
% labels, the optional argument to \subfloat must be present. If a
% subcaption is not desired, just leave its contents blank,
% e.g., \subfloat[].


% An example of a floating table. Note that, for IEEE style tables, the
% \caption command should come BEFORE the table and, given that table
% captions serve much like titles, are usually capitalized except for words
% such as a, an, and, as, at, but, by, for, in, nor, of, on, or, the, to
% and up, which are usually not capitalized unless they are the first or
% last word of the caption. Table text will default to \footnotesize as
% the IEEE normally uses this smaller font for tables.
% The \label must come after \caption as always.
%
%\begin{table}[!t]
%% increase table row spacing, adjust to taste
%\renewcommand{\arraystretch}{1.3}
% if using array.sty, it might be a good idea to tweak the value of
% \extrarowheight as needed to properly center the text within the cells
%\caption{An Example of a Table}
%\label{table_example}
%\centering
%% Some packages, such as MDW tools, offer better commands for making tables
%% than the plain LaTeX2e tabular which is used here.
%\begin{tabular}{|c||c|}
%\hline
%One & Two\\
%\hline
%Three & Four\\
%\hline
%\end{tabular}
%\end{table}


% Note that the IEEE does not put floats in the very first column
% - or typically anywhere on the first page for that matter. Also,
% in-text middle ("here") positioning is typically not used, but it
% is allowed and encouraged for Computer Society conferences (but
% not Computer Society journals). Most IEEE journals/conferences use
% top floats exclusively. 
% Note that, LaTeX2e, unlike IEEE journals/conferences, places
% footnotes above bottom floats. This can be corrected via the
% \fnbelowfloat command of the stfloats package.

\section{Optimal oracle shrinkage parameters}  \label{sec:optim1}

Define {\paino scale} measures of $\M \in \PDH_\pdim$ as 
\beq \label{eq:scale} 
\eta =  \tr( \M )/\pdim \quad \mbox{and} \quad \eta_2 =  \tr( \M^2 )/\pdim . 
\eeq 
%Note that $\M$ has a decomposition $\M =\eta \V$, where $\V$ is {\paino shape matrix} (or normalized covariance matrix), 
%\beq \label{eq:shape} 
%\V =  \M/\eta =   \pdim \M/ \tr(\Mn). 
%\eeq 
%verifying $\tr(\V)=\pdim$. 
An important measure in our future developments is the following measure of sphericity \cite{srivastava2005some},  
\beq \label{eq:gamma} 
\gamma = \frac{\eta_2}{\eta^2} = \frac{\pdim \tr(\M^2)}{\tr(\M)^2}   . 
\eeq 
Statistic $\gamma$  measures how close the covariance matrix is to a scaled identity matrix.  
%From Cauchy-Schwartz inequality it follows that 
%\[
%\tr(\M) ^2= | \langle \bom \lambda,  \bo 1  \rangle |^2 =\Big( \sum_{i=1}^\pdim \lambda_i \cdot 1 \Big)^2  \leq \pdim  \cdot \sum_{i=1}^\pdim \lambda_i^2  = \pdim \tr(\M^2), 
%\] 
%where $\bom \lambda = (\lambda_1,\ldots,\lambda_\pdim)^\top$ is a $\pdim$-vector that contains the eigenvectors of $\Mn$. 
%Consequently, 
 It verifies $\gamma \geq 1$ and $\gamma=1$ if and only if $\Mn= c \I$ %(i.e., $\bom \lambda = c \bo 1$) 
 for some $c>0$. 
 

%\begin{remark}  
 The parameters $\eta$ and $\gamma$ are elemental in our developments. As is shown in Theorem~\ref{th:beta0_ELL}, the optimal shrinkage parameter pair $(\al_o,\be_o)$ for elliptical distributions depends on the true covariance matrix $\M$ only through $\eta$ and $\gamma$.  Simple "plug-in" estimates of $(\al_o,\be_o)$ can then be  obtained by simply replacing $(\eta, \gamma)$ with  their estimates.  Finding accurate and consistent estimators of the shrinkage parameters is then a considerably simpler task than in the general case of Theorem~\ref{th:oracle_LW}. %     since we only need to find estimates of $\eta$ and $\gamma$ which are  simple scalar functions of $\M$. 
% \end{remark} 

Next theorem provides the expression for the oracle shrinkage parameters in the  case of sampling from an unspecified $\pdim$-variate distribution with finite $4$th order moments. 


\begin{theorem}  \label{th:oracle_LW}    
%Let $\{\x_i\}_{i=1}^\ndim \iidsim  \mathcal E_\pdim(\bo 0, \M, g)$ 
%and assume that elliptical population possesses finite    4th-order moments.   
Let $\{\x_i\}_{i=1}^\ndim$ denote a random sample from any $\pdim$-variate distribution (not necessarily elliptical distribution) with finite 4th order moments.  Then the oracle parameters in \eqref{eq:RegSCM_oracle} are  
%Then the oracle parameters in \eqref{eq:RegSCM_oracle} are  
\begin{align} \label{eq:beta0ver2} 
\be_o  %&=  \frac{ \| \M -  \eta  \I \|^2_{\rm F}}{  \| \M -  \eta  \I \|^2_{\rm F}  +  \E\big[ \| \S -  \M \|^2_{\rm F}\big] } =  \frac{ \| \M -  \eta  \I \|^2_{\rm F}}{\E\big[ \| \S -  \eta \I \|^2_{\rm F}\big] }    \notag \\�
&=  \frac{   \pdim (\gamma-1) \eta^2  }{   \E \big[\tr \big(\S^2\big) \big]  -  \pdim \eta^2 }    
%\end{align} 
\quad \mbox{and} \quad  
%\beq \label{eq:alpha0}
\al_o =  (1 - \be_o) \eta % (1 - \be_o) \tr(\Mn)/\pdim. 
\end{align} 
where $\eta$  and $\gamma$ are defined in \eqref{eq:scale} and \eqref{eq:gamma}, respectively.  
The value of MSE at the optimum is 
\beq \label{eq:MSEopt}
\MSE(\bo S_{\al_o,\be_o})  = \| \M - \eta \I \|^2_{\Fr}  ( 1- \be_o). 
\eeq 
The optimal   $\be_o$ is  always in the range $[0,1)$. %The above results are valid for both the real-valued and complex-valued observations. 
\end{theorem} 

\begin{proof}  
It was shown in \cite[Theorem~2.1]{ledoit_wolf:2004}  that 
\begin{align} \label{eq:beta0ver2_apu} 
\be_o  &=  \frac{ \| \M -  \eta  \I \|^2_{\rm F}}{  \| \M -  \eta  \I \|^2_{\rm F}  +  \E\big[ \| \S -  \M \|^2_{\rm F}\big] } 
%=  \frac{ \| \M -  \eta  \I \|^2_{\rm F}}{\E\big[ \| \S -  \eta \I \|^2_{\rm F}\big] }    
\end{align} 
and $\al_o =  (1 - \be_o) \eta$. The form of $\be_o$ in \eqref{eq:beta0ver2_apu} implies that $\be_o \in [0,1)$.  
We now show that \eqref{eq:beta0ver2_apu} can be expressed in the form \eqref{eq:beta0ver2}. 
First, we observe that 
\begin{align} 
 a_1 &= \E \big[ \| \Mor - \Mn \|^2_{\Fr} \big] =   \E \big[\tr \! \big(\Mor^2\big) \big] - 2 \E \big[ \tr\big(\Mor \Mn\big)\big] + \tr\big(\Mn^2\big) \notag 
 \\ 
 &=  \E \big[\tr \big(\Mor^2\big) \big]  -  \tr(\Mn^2) \label{eq:varrho}
\end{align}
where we used that $\E[\tr(\Mor \Mn)] = \Tr(\E[\Mor]\Mn) =  \tr(\Mn^2)$. 
The numerator of $\be_o$ in \eqref{eq:beta0ver2_apu} is 
\begin{align} 
a_2 &=\| \Mn -  \eta  \I \|^2_{\rm F} = \tr(\Mn^2) - (1/\pdim) \big\{\tr(\Mn) \big\}^2   \notag \\ 
&= \pdim( \eta_2 - \eta^2)  =  \pdim(\gamma-1) \eta^2  \label{eq:th_oracle_LW_apu}
\end{align} 
which shows that denominator of $\be_o$ is $a_1  + a_2 = \E \big[\tr\! \big(\Mor^2\big) \big]  - (1/\pdim) \big\{\tr(\Mn) \big\}^2  =  \E \big[\tr\! \big(\Mor^2\big) \big]  - \pdim \eta^2$.   These expressions for numerator and denominator of $\beta_o$ yield the assertion \eqref{eq:beta0ver2}  for $\be_o$. 
Write $L(\al,\be)=\E \big[ \| \Mcl  - \Mn \|_{\rm F}^2 \big]$ for the MSE. Note that 
\begin{align} 
L(\al,\be) %&= \E \Big[ \big \| \be \Mor  + \al \I  - \Mn \big\|_{\rm F}^2 \Big]  \notag \\ 
&=  \E \Big[ \big \| \al \I  + \be( \Mor - \Mn) - (1-\be) \Mn  \big\|_{\rm F}^2 \Big]  \notag
 \\ &= \al^2 \pdim + \be^2  a_1  + (1-\be)^2 \eta_2 \pdim   - 2\al(1-\be) \eta \pdim.  \label{eq:oracle_LW_apu} 
\end{align} 
The MSE at the optimum is 
\begin{align*} 
\MSE(\bo S_{\al_o,\be_o})   &= L((1-\be_o)\eta,\be_o)\\ &= \be_o^2 a_1 +  (1-\be_o)^2 \eta_2 \pdim  - (1-\be_o)^2 \eta^2 \pdim  \\ 
&= (1-\be_o^2)^2 \underbrace{ \{ \pdim(\eta_2 - \eta^2) \}} _{ \mbox{$=a_2$ by \eqref{eq:th_oracle_LW_apu}}}  + \be_o^2 a_1 \\ &= (1-\be_o)^2a_2 + (1-\be_o) \be_o a_2  \\
&= (1-\be_o) a_2,
\end{align*} 
where the 3rd identity follows as $\be_o a_1 = (1-\be_o) a_2$. This completes the proof.  
\end{proof}

%\begin{remark} \label{rem:CV} First,
Theorem has important implications.  First, since $\al_o = (1-\be_o)\eta$  is determined by the value of $\be_o \in [0,1)$, the optimal RSCM can be expressed simply 
as
\[ 
\bo S_{\al_o,\be_o}  =   \be_o \frac{1}{n} \sum_{i=1}^\ndim  \z_i \z_i^\top      +  (1-\be_o)\eta   \I.   
\] 
Since $\hat \eta=\tr(\S)/\pdim$ is a consistent estimator of $\eta=\tr(\M)/\pdim$ both in the conventional (fixed $\pdim$) and RMT asymptotic regime, we need  to simply  focus on finding a consistent estimator $\hat \be_o$ of $\be_o$. Consistent estimator of $\al_o$ is determined simply as 
$\hat \al_o = (1-\hat \be_o) \tr(\S)/\pdim$.   
%Thus only a single  penalty parameter $\be \in [0,1)$ needs to be tuned. 
%Finding an consistent estimator of  $\be_o$ then generally requires finding estimate of  
%$ \gamma$  and  $ \E \big[\tr \big(\Mor^2\big) \big]$ as well.   
%\end{remark} 


Ledoit and Wolf \cite{ledoit_wolf:2004} showed  that 
the following estimate 
% \beq  \label{eq:LW_betastar} 
\[
\hat \be_{{\tiny \mbox{LW}}}^* =  1 -  \frac{ \sum_{i=1}^\ndim \| \x_i \x_i^\top - \S \|_\Fr^2 }{\pdim \ndim^2 (\hat \gamma-1)}    
=  1 -  \dfrac{   \frac{1}{\ndim\pdim} \sum_{i=1}^\ndim \| \x_i \|_2^4  -  \hat \eta_2}{\ndim (\hat \gamma-1)}    
\]
%\eeq
where $ \hat \gamma=  \hat \eta_2/\hat \eta=p\tr(\S^2)/\tr(\S)$ and $\hat \eta_2 =  \tr( \S^2 )/\pdim$, converges to $\be_o$ in \eqref{eq:beta0ver2} in   probability   under RMT regime (R1) and  (R2) 
when sampling from a distribution $\x \sim F$ with finite $4$th-order moments. 
The authors of \cite{ledoit_wolf:2004}  then proposed to estimate the shrinkage parameters using  
\begin{align*} %\label{eq:LW_beta} 
\hat \be_0^{{\tiny \mbox{LW}}} &= \max(0, \hat \be_{{\tiny \mbox{LW}}}^*)  \quad \mbox{and} \quad \hat \al_0^{{\tiny \mbox{LW}}}   = (1-  \hat \be_0^{{\tiny \mbox{LW}}}) \tr(\S)/\pdim , 
\end{align*} 
where the $\max$ constraint ensures that the final estimate  remains on the interval $[0,1]$. The RSCM based on the above penalty parameters %is denoted by $\S_{\mathrm{LW}}$ and 
is referred hereafter as {\paino LW-RSCM estimator}.   

\section{Optimal oracle shrinkage parameters: the elliptical case} \label{sec:optim2}


% no \IEEEPARstart
Assume now that  $\x_1,\ldots, \x_n$ are independent and identically distributed (i.i.d.)  random vectors from 
a centered elliptical distribution with  mean vector $\E[\x]=\bo 0$ and positive definite covariance matrix $\M = \E[ \x \x^\top] $, denoted $ \mathcal E_\pdim(\bo 0,\M,g)$.  For a review of elliptical distributions, see  \cite{muirhead:1982,fang_etal:1990,ollila2012complex}.  %fang_zhang:1990}
The probability density function (p.d.f.) of   $\x \sim \mathcal E_\pdim(\bo 0,\M,g)$ is 
% denote a random vector from a centered elliptically symmetric distribution with  mean vector $\bom \mu$ and positive definite covariance matrix $\M \in \PD$. The p.d.f. of $\x$ is then of the form 
\[
f(\z) = C_{\pdim,g} 
  |\M|^{-1/2} g\big( \z^\top \M^{-1} \z\big)
 \] 
where  $g: [0,\infty) \to [0,\infty)$ is  a fixed function, called the  {\paino density generator}, that is independent of $\x$ and $\M$, and $C_{\pdim,g}$ is a normalizing constant ensuring that $f(\z)$ integrates to 1. Let $g$ be defined so that $\M$ represents the covariance matrix  of $\x$.  For example, the $\pdim$-variate Gaussian distribution, denoted $\x \sim \mathcal N_\pdim(\bo 0,\M)$, is a member in this class with density generator $g(t)=\exp(-t/2)$.  
% which means that 
%$ \int_0^\infty t^{\dim/2} \,  g(t)   \mathrm{d} t  = \pdim$. 
%If the first moments exists, then $\E[\x ] = \bo 0$
%We consider the case that $\M$ 
%The density generator $g$ is assumed to verify $\int_0^\infty  \int_0^\infty t^{\dim/2+1} \,  g(t)   \mathrm{d} t  = \pdim$ which is equivalent to saying that 2nd-order moments exists and  that $\M= \E[\x \x^\top]$.  
%In this paper, we %estrict to a  subclass of 
As earlier in Theorem~\ref{th:oracle_LW},  we   
assume that the elliptical population possesses finite 4th-order moments. 
%Technically, this implies  that %the following  restriction on $g$: 
%This is equivalent to stating that the density generator needs to verify  
%\beq \label{eq:finite_4mom}
 %\int_0^\infty t^{\dim/2+1} \,  g(t)   \mathrm{d} t  < \infty ,
%\eeq 
%holds. 
%and $\int_0^\inty  \int_0^\infty t^{\dim/2+1} \,  g(t)   \mathrm{d} t  = k$. 
% Also  $\pdim$-variate $t$-distribution, denoted  $\x \sim t_\nu(\bo 0 , \M)$,  when the degrees of freedom 
 %$\nu>4$, is a valid  distribution in this (sub)class of elliptical distiributions.  



Recall that the {\paino kurtosis} of a zero mean random variable $x$  is defined as 
\[
 \mathrm{kurt}(x) = \frac{\E[ x^4]}{ (\E[x^2])^2} - 3 .
\]
The {\paino elliptical kurtosis parameter} \cite{muirhead:1982}  $\ka$ of a random vector $\x=(x_1,\ldots,x_\pdim)^\top \sim \mathcal E_\pdim(\bo 0,\M,g)$ is defined as 
 \beq  \label{eq:kappa} 
 \kappa =  \frac{ \E[r^4] }{\pdim(\pdim+2)}  -1  =  \frac 1  3 \cdot  \mathrm{kurt}(x_i),  
 \eeq 
 where  $r$ denotes the (2nd order) modular variate of the elliptical distribution,  defined as $r= \sqrt{ \x^\top \M^{-1} \x}$.  
 The elliptical kurtosis %generalization of kurtosis for elliptical random vectors and 
 shares properties similar to kurtosis of a real random variable.  Especially, if    
$\x \sim \mathcal N_\pdim(\bo 0,\M)$, then $\ka=0$.  This is  obvious since the marginal distributions are Gaussian and  hence $\kappa = (1/3) \, \mathrm{kurt}(x_i)=0$.  
Another way to derive this is by noting that $ r^2=\x^\top \M^{-1} \x  \sim \chi^2_\pdim$ 
 and hence %the 2nd-moment of $\chi^2_\pdim$-distribution  is $\pdim(\pdim+2) $, so 
 $\E[r^4] = \pdim(\pdim+2)$. 
   %  and  $x_i$ denotes  (any) element of $\x \sim F_{\bo 0, \M}= \mathcal E(\bo 0, \M,g)$ 
    % where $\M=\cov(\x)$.   
    The importance of elliptical kurtosis parameter $\kappa$ is due to the fact that the $\pdim^2 \times \pdim^2$ covariance matrix of $\ve(\S)$ can be expressed as 
      \cite{muirhead:1982}:
 \begin{gather}
 \cov( \ve(\S))  =  \notag \\  
  \frac{ (1+\ka)}{\ndim}  (\I + \commat)  (\M \otimes \M) +  \frac{\ka}{\ndim} \,  \ve(\M) \ve(\M)^\top, \label{eq:cov_vecS} 
 \end{gather} 
 where $\commat$�  denotes the commutation matrix defined in the Introduction. Thus the elliptical kurtosis parameter $\ka$ along with the true covariance matrix $\M$ provide a complete description of the covariances between elements $S_{ij}$ and $S_{kl}$   of the SCM $\S$.  

In the next Lemma we derive the MSE of the SCM. 

\begin{lemma} \label{lem:NMSE}   Let $ \{ \x_i \}_{i=1}^\ndim \iidsim \mathcal E_\pdim(\bo 0,\M,g)$, where $\M = \cov(\x_i)$
and 4th-order moments exists. 
Then the MSE of $\S$ is 
\[
\MSE(\S)=  \frac \pdim \ndim \cdot  \eta^2  \Big\{  \kappa( 2 \gamma +  \pdim) +   \gamma+\pdim \Big\},
\]
and the {\paino normalized mean squared error (NMSE)} is 
\[
\mathrm{NMSE}(\S)  =   \frac{\E \big[ \| \S  - \M \|_\Fr^2 \big] }{ \| \M \|^2_\Fr} =  \frac 1 \gamma \cdot  \frac 1 \ndim  \Big\{  \kappa( 2 \gamma +  \pdim) +  \gamma+\pdim \Big\}.
\]
Furthermore, 
\begin{align*} 
\E [ \tr(\S^2)]  &= \MSE(\S) + \pdim \eta_2. 
%&= \frac{1 + \ka}{\ndim} \big( 2 \tr(\M^2)  + [\tr(\M)]^2 \big) + \Big ( 1 -  \frac 1 \ndim \Big) \tr(\M^2)    \\  %\notag 
%   &=   \frac \pdim  \ndim \cdot  \eta^2    \{  \ka ( 2 \gamma + \pdim) + (\ndim+1)\gamma + \pdim \}  % \label{eq:exp_trS2} �%\label{eq:exp_trSsq} 
%\E [ \tr(\S)^2] & =  \frac{1 + \ka}{\ndim} \big( 2 \tr(\M^2)  + [\tr(\M)]^2 \big) + \Big ( 1 -  \frac 1 \ndim \Big) \tr(\M)^2 % \label{eq:exp_trSsq} 
 \end{align*} 
Above  $\eta, \gamma$ and  $\ka$ are defined in \eqref{eq:scale}, \eqref{eq:gamma} and \eqref{eq:kappa}, respectively. 
%where $\gamma$ and $\kappa$ are defined in \eqref{eq:gamma} and \eqref{eq:kappa}, respectively.  
\end{lemma} 

\begin{proof} 
Since $\S$ is unbiased, so $\E[\S] = \M$, it holds that 
\beq \label{eq:lem:NMSE:apu1} 
\mathrm{MSE}(\S) =  \tr \{ \cov( \ve(\S)) \} ,
\eeq 
where $\cov(\ve(\S))$ has the expression stated in \eqref{eq:cov_vecS}.  
Then recall the following results: 
$\tr(\A \otimes \B) = \tr(\A) \tr(\B) $,  $\tr \{  \ve(\A) \ve(\B)^\top \} = \tr(\A \B)$ for any square matrices $\A$ and $\B$ of same order; see e.g., \cite{magnus_neudecker:1999}. 
These imply that 
\beq \label{eq:lem:NMSE:apu2} 
\tr( \M \otimes \M)= \tr(\M)^2, \, \,  \tr \{\ve(\M) \ve(\M)^\top \} = \tr(\M^2).
\eeq 
It is also easy to  show that  
\beq \label{eq:lem:NMSE:apu3} 
\tr\big\{ \commat (\M \otimes \M) \big\} = \tr(\M^2)
\eeq
by recalling the definition of the commutation matrix and  the property $\tr(\A \otimes \B) = \tr(\A) \tr(\B)$. % We omit the proof of this straightforward result.
 Using  \eqref{eq:lem:NMSE:apu1} -   \eqref{eq:lem:NMSE:apu3}, then yield the stated expression for $\MSE(\S)$. The expression for NMSE is obtained by dividing $\MSE(\S)$  by $\tr(\M^2) = \pdim \eta_2$. The last results follows as
\begin{align*} 
\MSE(\S) &=\E[ \| \S - \M \|_\Fr^2 ] = \E[ \tr\{ (\S - \M) (\S - \M)\}] \\ &= \E[ \tr(\S^2) - 2 \tr(\S \M) + \tr(\M^2)]  \\ &= \E[\tr(\S^2)] - \pdim \eta_2
\end{align*} 
by using that $\E[\tr(\S \M)] = \tr(\E[\S] \M) = \tr(\M^2) = \pdim \eta_2$.
% Hence we get 
%\begin{align*} 
%\E[\tr(\S^2)]  &= \MSE(\S) +  \pdim \eta_2 \\ 
%&= \frac \pdim \ndim \cdot  \eta^2  \Big\{  \kappa( 2 \gamma +  \pdim) +   (\gamma+\pdim) \Big\} + \pdim \eta_2 \\ 
%&=  \frac \pdim \ndim \cdot  \eta^2  \Big\{  \kappa( 2 \gamma +  \pdim) +   (\gamma+\pdim) + \ndim (\eta_2/\eta^2 ) \Big\} 
%\end{align*} 
%which is the stated expression since $\eta_2/\eta^2 = \gamma$. 
\end{proof}  

Next theorem states that the oracle parameters derived in Theorem~\ref{th:oracle_LW} can be written in a much simpler form when sampling from an 
elliptically symmetric distribution.  
 
 \begin{theorem}  \label{th:beta0_ELL} Let $\{\x_i\}_{i=1}^\ndim \iidsim  \mathcal E_\pdim(\bo 0, \M, g)$ 
and assume that elliptical population possesses finite    4th-order moments.   
Then the oracle  parameters $(\al_o,\be_o)$  that minimize the MSE 
are  
\begin{align*} 
\be_o^{Ell} &= \dfrac{ \gamma-1}{   (\gamma-1) \, + \, \gamma \cdot  \mathrm{NMSE}(\S) }  \\
&= \dfrac{ \gamma-1 }{   \gamma-1 +  (1/\ndim)\{  \kappa( 2 \gamma +  \pdim) +  \gamma+\pdim \} }  
\end{align*}
and $\al_o^{Ell} = (1-\be_o^{Ell}) \eta$.   
%where $\eta=\tr(\M)/\pdim$ and $\gamma=\pdim \tr(\M^2)/\tr(\M)^2$ and $\ka$ is defined in \eqref{eq:kappa}.  
%The value of MSE at the optimum is 
%$\MSE(\S_{\al_o,\be_o})  = \| \M - \eta \I \|^2_{\Fr}  ( 1- \be_o^{Ell})$. 
\end{theorem} 

\begin{proof} Using Lemma~\ref{lem:NMSE}, the denominator of $\be_o$ is 
\begin{align*} 
\E\big[&\tr\! \big(\S^2\big) \big]  -  \pdim \eta^2  \\ 
&= \MSE(\S) +\pdim  \eta_2 - \pdim \eta^2 \\ &= \pdim \eta^2 \{ \MSE(\S)/(\pdim \eta^2) + \gamma-1    \}  \\�
&=  \pdim \eta^2 \{  \gamma \cdot \mathrm{NMSE}(\S)  +  \gamma-1 \}, 
%&=  (\pdim/\ndim)  \eta^2    \{  (1 + \ka)  \big( 2 \gamma  +  \pdim \big) + \big ( n -  1 \big) \gamma \}  -   \pdim \eta^2 \\ 
%&=(2/\ndim)  \pdim \eta^2    \{  (\ndim/2) (\gamma-1) +  (1/\ndim)\{ \kappa (2 \gamma + \pdim) + (\gamma + \pdim) \}
\end{align*} 
where the last idenitity follows as $\mathrm{NMSE}(\S)= \MSE(\S)/\| \M\|_F^2 = \MSE(\S)/(\pdim\eta_2)$  and recalling that $\gamma=\eta_2/\eta^2$. 
Substituting this expression into \eqref{eq:beta0ver2} yields the first assertion for  $\be_o$. The second assertion follows by 
recalling the expression for $\mathrm{NMSE}(\S)$�  from Lemma~\ref{lem:NMSE}.  
%The complex-valued case follows similarly. 
\end{proof} 

It is not surprising that $\beta_o$ and hence also $\al_o$ depend on the functional form of the elliptical distribution (i.e., on density generator $g$) only via elliptical kurtosis parameter $\ka$.  Specifying the elliptical distribution (e.g., Gaussian, $t$-distribution, etc), also specifies the value of  $\kappa$.   
For example, when sampling from the Gaussian distribution, the elliptical kurtosis parameter is $\kappa=0$, 
%and consequently $\be_o^{Ell}$ in Theorem~\ref{th:beta0_ELL} reduces to  
%\beq \label{eq:beta0Gau}
%\be_o^{Gau} = \frac{ (\gamma-1)}{ (\gamma-1)  +  (\gamma +  \pdim)/\ndim } .
%\eeq 
%An estimator of  $\be_o^{Gau}$  is obtained by substituting an estimator $\hat \gamma$ %and $\hat \eta$ 
%in place of $\gamma$ in \eqref{eq:beta0Gau}. Recall that an estimator of $\al_o$ is then simply obtained as  $\hat \al_o = (1-\hat \be_o^{Gau}) \tr(\S)/\pdim$.  
%In the general  case, when one does not specify the underlying elliptical density, 
but since we do not assume any particular elliptical distribution, we need to %substitute  
find a  consistent estimator of the elliptical  kurtosis parameter  $\hat \kappa$ as well.  
% in place of $\kappa$ to obtain an estimate 
%$\hat \be_o^{Ell}$. 


\section{Consistent estimation of the oracle parameters} \label{sec:estim} 

Let $\{\x_i\}_{i=1}^\ndim \iidsim \mathcal E_\pdim(\bo 0, \M , g)$, where $\cov(\x)=\M$ and assume that the $4$th-order moments exists. 
In this section, we address the important topic of how to obtain consistent estimators of the unknown parameters $\eta, \gamma$ and $\kappa$. 

%Note that $\hat \gamma^{Gau}$ in   \eqref{eq:hatgamma} is consistent to $\gamma$  under Gaussian assumption.  
%We now provide two consistent alternative estimators for $\gamma$. 
First we recall that the sample sign covariance matrix,   
%\cite{visuri_etal:2000} is  
defined 
as 
\[
\S_{sgn}  = \frac{1}{\ndim} \sum_{i=1}^\ndim \frac{\x_i \x_i^\top}{\| \x_i \|^2} , 
\]
%The sample sign covariance matrix 
is well-known to be highly robust although it is not a consistent estimator of the covariance matrix \cite{croux2002sign}.  However,  the following result from \cite[Lemma~4.1]{zhang2016automatic} shows that it can be used to estimate the parameter $\gamma$. 

\begin{lemma}  \label{lem:eta2_ell1}  Let $\{\x_i\}_{i=1}^\ndim \iidsim \mathcal  E_\pdim(\bo 0,\M,g)$. Then  
\beq \label{eq:gammahat}
\hat \gamma = \pdim \tr \big(\S_{sgn}^2 \big) - (\pdim/\ndim) 
\eeq 
is a consistent estimator of $\gamma = \pdim \tr(\M^2)/\tr(\M)^2$ under assumption (R1) and (R2). 
\end{lemma}

Note that $\hat \gamma$ is a robust and distribution-free estimator of $\gamma$. % and does not, in general, require assumption (R2). 
The optimum parameter $\beta_o^{Ell}= \beta_o^{Ell}(\gamma, \ka)$ depends on $\gamma$ and $\ka$. Hence a {\paino plug-in estimator}, 
\[ 
\hat \beta_o^{Ell}= \beta_o^{Ell}(\hat \gamma, \hat \ka),
\] 
where $\hat \gamma$ and $\hat \ka$ are consistent estimators of $\gamma$ and $\kappa$, is a consistent estimator of  $\beta_o$ as well. % of Theorem~\ref{th:beta0_ELL}.  
%Thus we still need to 
 %$\ka$ and is not directly useful unless we assume a known elliptical distribution. 
  %in which case $\ka$ can be calculated analytically. 
  A natural estimate of $\ka$  is the  conventional sample average, 
\beq \label{eq:kappahat}
\hat \kappa =  \max \Big( - \frac{2}{p+2} \, ,  \,  \frac{1}{3 \pdim} \sum_{j=1}^\pdim  \hat k_j \Big), % \quad  \bigg(\mbox{resp. } \hat \kappa =  \frac{1}{2 \pdim} \sum_{j=1}^\pdim  \hat k_j, \bigg)
\eeq 
where   $\hat k_j =  m^{(4)}_j/ \big(m_j^{(2)} \big)^2 - 3$   is the sample kurtosis of the $j$th variable  and 
%f the  kurtosis of the $j$th variable, i.e.,  
%\[ 
%\hat k_j  = \left\{ \frac{\ndim^2 (n+1)}{(n-1)(n-2)(n-3)} \right\} \frac{ m^{(4)}_j}{m_j^{(2)}}   - 3 \frac{(n-1)^2}{(n-2)(n-3)} 
%\]
$m^{(q)}_j = \frac 1 \ndim \sum_{i=1}^\ndim  (x_{ij})^q$ 
denotes the $q$th order sample moment, %(and note that $s^2_j = m^{(2)}_j$ is the sample variance of $j^{th}$ variable),
 $j=1,\ldots, \pdim$.   Above  the $\max$ constraint ensures that the final estimate  $\hat \ka$ does not exceed the theoretical lower bound \cite{bentler1986greatest},  $-2/(p+2)$ of elliptical kurtosis parameter $\ka$.  
 The estimate $\hat \kappa$ is a consistent estimator of the elliptical kurtosis $\kappa$ both in the conventional and RMT regime.
 
 
   %We then substitute the estimate  $\hat \ka$ in place of $\ka$ in  equation for $\beta_o^{Ell}$ of Theorem~\ref{th:beta0_ELL}. 
 
%Thus we expect that $\hat \gamma^{Ell1}$  provides better accuracy for heavy-tailed distributions.  
 
We can now define the {\paino Ell-RSCM estimator} as the regularized SCM  based on the following estimated optimal shrinkage parameters 
\begin{align} 
 \hat \beta_o^{Ell} &=  \max \left(0, \dfrac{T}{  T +  (1/\ndim)\{ \hat \kappa (  2 \hat \gamma +  \pdim) +   \hat \gamma + \pdim \} } \right)   \label{eq:Ell-RSCM}  \\ 
 \hat \al_o^{Ell} &= (1-  \hat \be_0^{Ell}) \tr(\S)/\pdim  \notag
\end{align} 
where $T = \hat \gamma -1$ and $\hat \gamma$  and $\hat \kappa$ are defined in \eqref{eq:gammahat} and  \eqref{eq:kappahat}, respectively,

 %(resp. $T^{Ell1} = n( \hat \gamma^{Ell1} -1)$) in the real-valued (resp. complex-valued) case. 
%We define  $\hat \beta_o^{Ell2}$ similarly as in \eqref{eq:Ell-RSCM} but 
%$T^{Ell1}$ and $\hat \gamma^{Ell1}$ replaced by  $T^{Ell2}$ and $\hat \gamma^{Ell2}$, respectively. 


\section{Simulation study} \label{sec:simul} 

We conduct a small simulation study to investigate the performance of RSCM estimators in terms of their finite sample NMSE. %normalized mean squared error (NMSE).  We include the LW-RSCM and Ell-RSCM estimators. % defined in  \eqref{eq:LW_beta} %\eqref{eq:Gau-RSCM} 
%and \eqref{eq:Ell-RSCM}, respectively.  
Each simulation is repeated 10000 times and the NMSE is computed (averaged of Monte-Carlo runs) for each  RSCM estimator.  Theoretical oracle MSE value derived in  \eqref{eq:MSEopt} and normalized by $\| \M \|^2_{\Fr}$ is used as a benchmark lower bound for empirical NMSE values.  This is shown in the figures as solid black line. 






%We also compute the Oracle-RSCM estimator 
%defined as RSCM estimator $\S_{\al,\be}$ that uses the  true oracle value  $(\al_o^{Ell},\be_o^{Ell})$  given in \eqref{eq:Ell-RSCM}.  
%Naturally, we will expect that the Oracle-RSCM estimator will  have similar performance with the theoretical NMSE value, but it gives better benchmark (for any $\ndim$ and $\pdim$) for the  best attainable performance. 

\subsection{AR(1) covariance matrix} 

%  $\S_{\hat \be}$, where $\hat \be_o$ 
%is Ell1, Ell2, Gau or LW estimator of $\be_0$ as defined earlier.  

 In the first experiment, an autoregressive covariance structured  is used. We let $\M$ be the covariance
matrix of a Gaussian AR$(1)$ process, 
%Toeplitz form
\[
[\M]_{ij}  = \varrho^{|i-j|},  \quad r \in (0,1) . 
\]
Note that $\M$ %is in fact a shape matrix as it 
verifies $\eta = \tr(\M)/\pdim = 1$. When $\varrho$ is close to $0$, then $\M$ is 
close to an identity matrix and when $\varrho$  tends to $1$, $\M$ tends to  a singular matrix of rank 1. % $\bo 1 \bo 1^\top$  
Thus the theoretical value $\be_o$ is close to $0$ for small values of $\varrho$, i.e., when the true covariance matrix is close to the target $\I$, 
and $\be_0 \approx 1$ for $\varrho$  close to $1$. Dimension is  fixed at $\pdim = 100$  and  $\ndim$ is allowed to vary from $0.2 \cdot \pdim $ to $1.2 \cdot \pdim$. 

Figure~\ref{fig:AR1_p100}  depicts the NMSE performance  %averaged over MC trials 
when the samples are drawn from a Gaussian distribution (upper panel) and 
 a multivariate $t_\nu$-distribution with $\nu=8$ degrees of freedom (lower panel).  
%as well as the  the graph of the average estimate $\hat \be_o$ (lower panel) as a function $c^{-1}=\pdim/\ndim$ for $r=0.1, 0.4$ and $r=0.8$.    
Several conclusions can be drawn from these figures. First, when $\varrho=0.1$ and thus $\M$ is close to the shrinkage target matrix $\I$,  Ell-RSCM 
%and Ell2-RSCM 
estimators outperform the LW-RSCM estimator.  Especially, when the ratio $\ndim/\pdim$ is small (i.e., $\pdim$ larger than $\ndim$), we observe the largest performance differences in favor of Ell-RSCM. 
%Gau-RSCM estimator is performing only slightly better compared to Ell-RSCM estimators. This can also be seen on the $\hat \be_o$ curves, where Gau-RSCM estimator is consistently closer to the oracle value $\be_o^{Ell}$. 
%Second, for all considered values of $r$, Gau-RSCM and both of the Ell-RSCM estimators are observed to have essentially equivalent performance in terms of NMSE. 
Second, when the true $\M$ starts to deviate significantly from the identity target matrix $\I$ (i.e., $\varrho=0.4$), LW-RSCM and Ell-RSCM estimator have similar performance especially for large values of $\ndim/\pdim$.  
%This is not suprising, but as it should be, i.e., for $\ndim > \pdim$, shrinkage is not as necessary and $\be_o$ is in general large.  
 Third, when the samples are drawn from $t_8$-distribution, the performance of LW-RSCM estimator is seen to deterioritate in comparison to the proposed Ell-RSCM estimator.  Indeed very large differences are witnessed in NMSE between the estimators especially when $\ndim/\pdim < 0.5$.  

%This is  somewhat expected, since now the shrinkage target matrix deviates  significantly from the true covariance matrix $\M$, so very similar large values of $\be_o$ are chosen for all estimators. As a consequence, all of the RSCM estimators are getting closer to the SCM. Indeed this can be seen in the graphs of $\be_o$  which are overlapping for $r= 0.8$, but differ considerably for $r=0.1$. 
%In conclusion, one can state that the tuned estimators for elliptical samples, so Gau-RSCM and Ell-RSCM estimators have better performance of LW estimator. Tuned %estimators have signicantly better performance when the shrinkage target matrix is close to the true covariance matrix (the case $r=0.1$). 

%../Springer_book/Part3/pics
\begin{figure}
\subfigure[$\varrho=0.1$]{\includegraphics[width=0.24\textwidth]{gausim_p100_rho0dot1_aveNMSE.pdf}}
\subfigure[$\varrho=0.4$]{\includegraphics[width=0.24\textwidth]{gausim_p100_rho0dot4_aveNMSE.pdf}}
\subfigure[$\varrho=0.1$]{\includegraphics[width=0.24\textwidth]{v8_p100_rho0dot1_aveNMSE_new.pdf}}
\subfigure[$\varrho=0.4$]{\includegraphics[width=0.24\textwidth]{v8_p100_rho0dot4_aveNMSE_new.pdf}}
% \includegraphics[width=0.25\textwidth]{./pics/gausim_p100_rho0dot4_beta_mean_new.pdf}}
%\includegraphics[width=0.25\textwidth]{./pics/gausim_p100_rho0dot1_beta_mean_new.pdf}
%\centerline{ (a) $r=0.1$}
%\centerline{\includegraphics[width=0.25\textwidth]{./pics/gausim_p100_rho0dot4_aveNMSE.pdf} \includegraphics[width=0.25\textwidth]{./pics/gausim_p100_rho0dot4_beta_mean_new.pdf}}
%\centerline{ (b) $r=0.4$}
\caption{AR(1) process: Comparison of covariance estimators when $\pdim = 100$ and  $\varrho \in \{ 0.1, 0.4\}$ and the samples are from Gaussian distribution (upper panel) 
and $t_\nu$-distribution with $\nu=8$ degrees of freedom (lower panel).}   \label{fig:AR1_p100}
\end{figure}


\subsection{Largely varying spectrum} 

Our next study follows the set-up in \cite{zhang2016automatic} in which $\M$ has one
(or a few) large eigenvalues.
In the first set-up,  $\M$ is a diagonal matrix of size 
$50 \times 50$, where $m$ eigenvalues are equal to $1$ and the remaining $50-m$ eigenvalues are 0.01.
For the case $\ndim = \pdim=50$, Figure~\ref{fig:tengsim1} depicts the NMSE as a function of $m$  
%averaged over 10,000 Monte Carlo runs   
when sampling from a $t_\nu$ distribution with $\nu=8$ degrees of freedom. 
Ell-RSCM has excellent performance as its NMSE curve is essentially overlapping with  the theoretical NMSE curve. 
LW-RSCM estimator is performing poorly for all values of $m$ except at the extremes, i.e, when $m$ is either small or  large, in which case the covariance matrix $\M$  is close to an (scaled) identity matrix. 
  
%Figure~\ref{fig:tengsim1}(b) depicts the average $\hat \beta_0$ values along with the theoretical oracle value $\beta_o$. As can be seen Ell1-RSCM is consistenly attaining the %oracle $\beta_o$ value for all values of $m$. LW-RSCM is far from the oracle value for intermediate values of $m$.  

\begin{figure}
\centerline{\includegraphics[width=0.45\textwidth]{v8_p50_tengsim1_new_v2.pdf}}
%\subfigure{\includegraphics[width=0.48\textwidth]{../Springer_book/Part3/pics/v8_p50_tengsim1_beta_mean.pdf}}
\vspace{-0.2cm}
\caption{The covariance matrix $\M$ has $m$ eigenvalues equal to $1$ and $50-m$ eigenvalues equal to $0.01$. The samples are from $t_\nu$-distribution with $\nu=8$ degrees of freedom and $\ndim = \pdim = 50$. }   \label{fig:tengsim1}
\end{figure}


Next simulation set-up considers a very challenging scenario in which the spectrum of $\M$ 
consists of several different eigenvalues. %This set-up also  is similar to \cite{zhang2016automatic}. 
We consider the case that $\pdim=100$ and the covariance matrix  $\M$ has 30 eigenvalues equal to $100$, $40$ eigenvalues equal to $1$
and 30 eigenvalues of $0.01$. Samples are drawn from $t_\nu$ distribution with $\nu=8$ degrees of freedom. 
% and results are averaged over 10,000 Monte Carlo runs. 
The NMSE curves shown in Figure~\ref{fig:tengsim1} illustrate the huge advantage of the proposed Ell-RSCM  over the LW-RSCM estimator. In fact, in this scenario the LW estimator fails and it assigns  $\hat \beta_o= 0$ for all values of $\ndim$. 
 %Again the Ell1-RSCM estimator has slight but consistent advantage over Ell2-RSCM estimator and 
 Again the Ell-RSCM estimator  reaches  near oracle performance and thus there is not much space for improvements.
 % Figure~\ref{fig:tengsim1}(b) which depicts the curve of the average estimate $\hat \beta_o$ shows the solid performance of our proposed estimators.  
  
\begin{figure}
\centerline{\includegraphics[width=0.45\textwidth]{v8_p100_tengsim2_v2.pdf}}
\vspace{-0.2cm}
%\subfigure{\includegraphics[width=0.48\textwidth]{../Springer_book/Part3/pics/v8_p100_tengsim2_beta_mean.pdf}}
\caption{ The covariance matrix  $\M$ has 30 eigenvalues equal to $100$, $40$ eigenvalues equal to $1$
and 30 eigenvalues equal to $0.01$.  The samples are from $t_\nu$-distribution with $\nu=8$ degrees of freedom and �$\pdim=100$.}   \label{fig:tengsim2}
\end{figure}



\section{Conclusion}

We proposed an optimal regularized sample covariance matrix estimator, called Ell-RSCM estimator,   
which is suitable for high-dimensional problems and when sampling from an unspecified elliptically symmetric  distribution. The estimator is based on consistent estimators (under RMT regime) of the optimal shrinkage parameters that minimize the MSE. It smartly exploits  elliptical theory such  as the knowledge of the form of  MSE of the SCM when sampling from an elliptical population.  Our simulation studies illustrated the advantage of the proposed Ell-RSCM over the Ledoit-Wolf (LW-)RSCM estimator.  The performance differences were often significant.  



%\appendix

% conference papers do not normally have an appendix


% use section* for acknowledgment
%\section*{Acknowledgment}
%The authors would like to thank...





% trigger a \newpage just before the given reference
% number - used to balance the columns on the last page
% adjust value as needed - may need to be readjusted if
% the document is modified later
%\IEEEtriggeratref{8}
% The "triggered" command can be changed if desired:
%\IEEEtriggercmd{\enlargethispage{-5in}}

% references section

% can use a bibliography generated by BibTeX as a .bbl file
% BibTeX documentation can be easily obtained at:
% http://mirror.ctan.org/biblio/bibtex/contrib/doc/
% The IEEEtran BibTeX style support page is at:
% http://www.michaelshell.org/tex/ieeetran/bibtex/
%\bibliographystyle{IEEEtran}
% argument is your BibTeX string definitions and bibliography database(s)
%\bibliography{IEEEabrv,../bib/paper}
%
% <OR> manually copy in the resultant .bbl file
% set second argument of \begin to the number of references
% (used to reserve space for the reference number labels box)

%\bibliography{../IEEE/bibtex/IEEEabrv,../mybib/STATabrv,../mybib/STATbib,../mybib/ENGbib}

\input{rscmEUSIPCO17final.bbl}
% that's all folks
\end{document}




\begin{figure}
\centerline{\includegraphics[width=0.25\textwidth]{pics/v8_p100_rho0dot1_aveNMSE_new.pdf} \includegraphics[width=0.25\textwidth]{pics/v8_p100_rho0dot1_beta_mean_newer.pdf}}
\centerline{ (a) $r=0.1$}
\centerline{\includegraphics[width=0.25\textwidth]{pics/v8_p100_rho0dot4_aveNMSE_new.pdf}\includegraphics[width=0.25\textwidth]{pics/v8_p100_rho0dot4_beta_mean_newer.pdf}}
\centerline{ (b) $r=0.4$}
%\centerline{\includegraphics[width=0.4\textwidth]{../Springer_book/Part3/pics/v8_p100_rho0dot8_aveNMSE_new.pdf}\includegraphics[width=0.4\textwidth]{../Springer_book/Part3/pics/v8_p100_rho0dot8_beta_mean_newer.pdf}}
%\centerline{ (c) $r=0.8$}
%\centerline{\includegraphics[width=0.32\textwidth]{../Springer_book/Part3/pics/v8_p100_rho0dot1_aveNMSE_new.pdf} \includegraphics[width=0.32\textwidth]{../Springer_book/Part3/pics/v8_p100_rho0dot4_aveNMSE_new.pdf} \includegraphics[width=0.32\textwidth]{../Springer_book/Part3/pics/v8_p100_rho0dot8_aveNMSE_new.pdf}}   
%\subfigure[$r=0.1$]{\includegraphics[width=0.32\textwidth]{../Springer_book/Part3/pics/v8_p100_rho0dot1_beta_mean_newer.pdf}} 
%\subfigure[$r=0.4$]{\includegraphics[width=0.32\textwidth]{../Springer_book/Part3/pics/v8_p100_rho0dot4_beta_mean_newer.pdf}}
%\subfigure[$r=0.8$]{\includegraphics[width=0.32\textwidth]{../Springer_book/Part3/pics/v8_p100_rho0dot8_beta_mean_newer.pdf}}
%}
\caption{AR(1) process: Comparison of covariance estimators when $\pdim = 100$ and  $r = 0.1, 0.4$ and the samples are from $t_\nu$-distribution with $\nu=6$ degrees of freedom.}   \label{fig:AR1_p100_v6}
\end{figure}

\section{Proof of Lemma~\ref{lem:trace_S} } 
\subsection{Auxiliary results} %Proof of Lemma~\ref{lem:trace_S} } 


%An important identity involving the Kronecker product and  vec-operator is 
%\beq \label{eq:vec}
%\ve(\A \B \mC)=(\mC^\top \otimes \A) \ve(\B), 
%\eeq 
%where $\A,\B$ and $\mC$ are matrices such that the product $\A\B\mC$ is properly defined. 
%Also, some useful, well-known properties of Kronecker product are: 
%$(\A \otimes \B)(\mC \otimes \mD) = \A\mC \otimes \B\mD$ (where matrices $\A,\B,\mC$ and $\mD$ are such that the products $\A\mC$ and $\B\mD$ 
%are properly defined),  $(\A \otimes \B)^\hop = \A^\hop \otimes \B^\hop$, 
%$a \otimes \A= a \A$ and $\bo a^\hop \otimes \bo b= \bo b \bo a^\hop$. 
%The matrix derives its name from the property: 
%\beq
%\commat (\A \otimes \B) = (\B \otimes \A) \commat. \label{eq:com} 
%\eeq 
Let $S^{\pdim } =\{\z \in \R^\pdim \, : \, \| \z \| = 1 \}$ 
%and by $\C S^{\pdim -1} =\{\z \in \C^\pdim \, : \, \| \z \| = 1 \}$ 
denote the unit real %and complex 
 $\pdim$-sphere.   We now derive moments and expectations of quadratic forms of a random vector $\u$ uniformly distributed on $\S^\pdim$, denoted $\u \sim U(S^\pdim)$. 

 
  %, respectively.  

%In order to prove we need to derive moments of quadratic forms of random vector $\u$ possessing  with a uniform distribution on  $S^{\pdim-1}$.
 
\begin{lemma}[See e.g. { \cite[Lemma A.1]{ollila_oja_croux:2003}}]  \label{lem:u_4th_order}
Let $\bo u$ be a random vector with a uniform distribution on  $S^{\pdim-1}$.  Then 
$\E[\u ] = \bo 0$ and $\E[ \u \u^\top ] = (1/\pdim) \I$ 
and the only non-zero 4th order moments of $\u$ are 
\[
 \E[u_i^4] = \frac{3}{\pdim(\pdim+2)}   \quad  \mbox{and} \quad  \E[u_i^2 u_j^2] = \frac{1}{\pdim(\pdim+2)} %\quad \mbox{where $i\neq j$},  
\]
for all $i,j \in \{1,\ldots,\pdim\}$, $i \neq j$.    
\end{lemma}

\begin{lemma} \label{LemU} 
Let $\u_1$ and $\u_2$ be independent random vectors with a uniform distribution on  $S^{\pdim-1}$. 
\begin{itemize} 
\item[a)] $\E[ \u_1 \u_1^\top \otimes \u_1 \u_1^\top  ]  =   \frac{1}{\pdim(\pdim +2)}  \big\{  \Id +  \ve(\I)\ve(\I)^\top +  \commat \big\}$ 
\item[b)] $\E[ \u_1 \u_1^\top \otimes  \u_2 \u_2^\top  ]   =   (1/\pdim)^2 \Id $.  
\item[c)] For any $\pdim \times \pdim$ matrix $\M \in \PDH_\pdim$, it holds that  
\begin{align*}
\E \big[ (\u^\top_1 \M \u_1)^2\big] %= \frac{1}{\pdim(\pdim+2)}  \bigg\{ 2 \tr(\M^2) + [ \Tr(\M)]^2 \bigg\}  
&= \frac{1}{\pdim(\pdim+2)}\{ 2 \pdim \eta_2 + \pdim^2 \eta^2\}  \\
   \E \big[ (\u^\top_1 \M \u_2)^2\big]  &= \frac{1}{\pdim^2} \Tr(\M^2)   =  \frac{\eta_2}{\pdim} 
\end{align*}
where $\eta_2 = \tr(\M^2)/\pdim$ and $\eta = \tr(\M)/\pdim$. 
\end{itemize} 
\end{lemma}

\begin{proof}
a) and b) part follow directly  by applying Lemma~\ref{lem:u_4th_order} to the components of the the matrices $\E[ \u_1 \u_1^\top \otimes \u_1 \u_1^\top  ]$ and 
$\E[ \u_1 \u_1^\top \otimes  \u_2 \u_2^\top  ] $. 
 %The latter expection also follows from  [b]-part of the Lemma and the fact that $\u_1$ and $\u_2$ are independent, so $\E[u_i^2 \tilde u_j^2] =   \E[u_i^2] \cdot \E[ \tilde u_j^2]  = (1/\pdim)^2$,  where the last identity follows using [a]-part of the Lemma. 
 Next  note that 
\begin{align*} 
&\E\big[ (\u^\top_1 \M \u_1)^2\big]   \cdot \pdim(\pdim +2)  
= \E [ \tr(  \u_1 \u_1^\top \M  \u_1 \u_1^\top \M) ] \cdot \pdim(\pdim +2)   \\
 &\overset{1)}= \ve(\M)^\top \E[ \u_1 \u_1^\top \otimes  \u_1 \u_1^\top  ] \ve(\M)     \cdot \pdim(\pdim +2)   \\ 
%&=  \ve(\M)^\top \big(   \E[u_i^2 u_j^2]  \Id + (\E[u_i^4]  -  2 \E[u_i^2 u_j^2] ) \ve(\I)\ve(\I)^\top + \E[ u_i^2 u_j^2]  \commat \big) \ve(\M)  
&\overset{2)}= \ve(\M)^\top \big\{  \Id +  \ve(\I)\ve(\I)^\top +  \commat \big\} \ve(\M)   \\
&=    \ve(\M)^\top \ve(\M) +  [\ve(\M)^\top\ve(\I) ]^2 +  \ve(\M)^\top  \commat \ve(\M)  \\ 
&\overset{3)}=  2 \tr(\M^2) +  [\tr(\M)]^2 
\end{align*}  
where in 1) we used the property $\tr(\A \B \A \B)=\ve(\B)^\top ( \A \otimes \A) \ve(\B)$, in 2) we used a)-part of the Lemma, and in 3) we utilized the 
following identities $\commat \ve(\A)= \ve(\A^\top) $ and $\tr(\A \B)=\ve(\A)^\top \ve(\B)$. Similarly, we get 
\begin{align*} 
\E \big[ (\u^\top_1 \M \u_2)^2\big]   &= \ve(\M)^\top \E[ \u_1 \u^\top_1 \otimes  \u_2 \u^\top_2  ] \ve(\M)   \\ &\overset{1)}= (1/\pdim)^2 \ve(\M)^\top \ve(\M) 
= (1/\pdim)^2 \tr(\M^2)
\end{align*}
where in 1) we used  b)-part of the Lemma.   The last identity follows using the property $\tr(\A \B)=\ve(\A)^\top \ve(\B)$.  
%\hspace{\stretch{1}} $\blacksquare$
 \end{proof} 
 
\subsection{Proof}
 
Let  $\{\x_i \}_{i=1}^\ndim  \iidsim \mathcal E_\pdim(\bo 0, \M, g)$,  By stochastic representation theorem,
 $\x_i =_d r_i \M^{1/2} \u_i$,  % \M^{-1/2} \x_i/\|\M^{-1/2} \x_i\| $ 
 where $\u_i$ has a uniform distribution on $S^{\pdim-1}$,  $\bo u_i$ is independent of modular variate $r_i$ and $r_i =_d \x_i^\top \M^{-1} \x_i$.  Hence we have that 
 %are   independent of $\u_i$ and $\u_i$-s form a random sample  
\begin{align} \label{eq:EtraceC^2}
 &\E \big[\tr \! \big(\S^2\big) \big] \cdot \ndim^2  \notag \\ &=  \E \bigg[\tr \!  \bigg\{  \M^{1/2}  \sum_{i=1}^\ndim r^2_i \u_i \u_i^\hop  \M^{1/2}  \cdot \M^{1/2}  \sum_{j=1}^\ndim r^2_j \u_j \u_j^\hop  \M^{1/2}  \bigg\} \bigg]  \notag \\ 
 &=    \sum_{i=1}^\ndim \sum_{j=1}^\ndim    \E \big[   r^2_i r^2_j \big]   \E \bigg[ \underbrace{ \tr \!   \bigg\{   \u_i \u_i^\hop  \M  \u_j \u_j^\hop   \M  \bigg\}  }_{  {\displaystyle =  ( \u_i^\hop \M \u_j )^2} } \bigg]   \notag  \\ 
 &=    \sum_{i=1}^\ndim  \E \big[   r^4 _i \big]   \E \big[   ( \u_i^\top \M \u_i)^2 \big] + 
 \sum_{i \neq j }\E \big[  r^2 _i \big]     \E \big[  r^2_j \big] \E \big[   |\u_i^\top \M \u_j|^2 \big] \   \notag 
  \\ 
 &=  \ndim  \E \big[  r^4_1 \big]   \E \big[   ( \u^\hop_1 \M \u_1)^2 \big] + 
 \ndim(\ndim- 1) \pdim^2   \E \big[  ( \u^\hop_1 \M \u_2)^2 \big] .
 \end{align} 
%In obtaining the second identity, we used the fact that $r_i$ is independent of $\u_i$ and $\tr(\A \B) = \tr(\B \A)$, in 3rd identity we used  that $r_i$ is independent of $r_j$, $i \neq %j$, and in obtaining the 4th identity, 
Above we only used that $\{\u_i\}$ and $\{r_i\}$ are i.i.d. samples and that $\u_i$ is also independent of $r_i$. In the last identity we used the fact that 
$\E[r_i^2] = \E[ \x_i^\top \M^{-1} \x_i] = \E[\tr(\M^{-1} \x_i \x_i^\top)] =  \tr(\M^{-1} \E[\x_i \x_i^\top] ) = \tr(\I)=\pdim$. 
 % as r.v. $r = \| \x \|$ where $\x \sim F_{0,\I}= \mathcal E_\pdim(\bo 0, \I, g)$ and 
%$\E[r_i^2]=\E_{\bo 0, \I} [\| \x \|^2] = \pdim$.  
 Now the stated equation for $\E [ \tr(\S^2)]$ follows simply by using Lemma~\ref{LemU}c) 
for  $ \E \big[   ( \u^\top_1 \M \u_2)^2 \big]$ and $\E \big[  (\u^\top_1 \M \u_1)^2 \big] $ 
and  recalling the expression for  $\sigma_1=1+\ka$ from \eqref{eq:psi1_SCM}. 
%\end{proof}




rscmEUSIPCO17final.bbl
% Generated by IEEEtran.bst, version: 1.14 (2015/08/26)
\begin{thebibliography}{10}
\providecommand{\url}[1]{#1}
\csname url@samestyle\endcsname
\providecommand{\newblock}{\relax}
\providecommand{\bibinfo}[2]{#2}
\providecommand{\BIBentrySTDinterwordspacing}{\spaceskip=0pt\relax}
\providecommand{\BIBentryALTinterwordstretchfactor}{4}
\providecommand{\BIBentryALTinterwordspacing}{\spaceskip=\fontdimen2\font plus
\BIBentryALTinterwordstretchfactor\fontdimen3\font minus
  \fontdimen4\font\relax}
\providecommand{\BIBforeignlanguage}[2]{{%
\expandafter\ifx\csname l@#1\endcsname\relax
\typeout{** WARNING: IEEEtran.bst: No hyphenation pattern has been}%
\typeout{** loaded for the language `#1'. Using the pattern for}%
\typeout{** the default language instead.}%
\else
\language=\csname l@#1\endcsname
\fi
#2}}
\providecommand{\BIBdecl}{\relax}
\BIBdecl

\bibitem{ledoit_wolf:2004}
O.~Ledoit and M.~Wolf, ``A well-conditioned estimator for large-dimensional
  covariance matrices,'' \emph{J. Mult. Anal.}, vol.~88, pp. 365--411, 2004.

\bibitem{ollila2014regularized}
E.~Ollila and D.~E. Tyler, ``Regularized {$M$}-estimators of scatter matrix,''
  \emph{{IEEE} Trans. Signal Process.}, vol.~62, no.~22, pp. 6059--6070, 2014.

\bibitem{zhang2016automatic}
T.~Zhang and A.~Wiesel, ``Automatic diagonal loading for tyler's robust
  covariance estimator,'' in \emph{IEEE Statistical Signal Processing Workshop
  (SSP'16)}, 2016, pp. 1--5.

\bibitem{chen2010shrinkage}
Y.~Chen, A.~Wiesel, Y.~C. Eldar, and A.~O. Hero, ``Shrinkage algorithms for
  mmse covariance estimation,'' \emph{{IEEE} Trans. Signal Process.}, vol.~58,
  no.~10, pp. 5016--5029, 2010.

\bibitem{couillet_mckay:2014}
R.~Couillet and M.~McKay, ``Large dimensional analysis and optimization of
  robust shrinkage covariance matrix estimators,'' \emph{J. Mult. Anal.}, vol.
  131, pp. 99--120, 2014.

\bibitem{pascal2014generalized}
F.~Pascal, Y.~Chitour, and Y.~Quek, ``Generalized robust shrinkage estimator
  and its application to stap detection problem,'' \emph{{IEEE} Trans. Signal
  Process.}, vol.~62, no.~21, pp. 5640--5651, 2014.

\bibitem{muirhead:1982}
R.~J. Muirhead, \emph{Aspects of Multivariate Statistical Theory}.\hskip 1em
  plus 0.5em minus 0.4em\relax New York: Wiley, 1982, 704 pages.

\bibitem{fang_etal:1990}
K.-T. Fang, S.~Kotz, and K.~W. Ng, \emph{Symmetric Multivariate and Related
  Distributions}.\hskip 1em plus 0.5em minus 0.4em\relax London: Chapman and
  hall, 1990.

\bibitem{ollila2012complex}
E.~Ollila, D.~E. Tyler, V.~Koivunen, and H.~V. Poor, ``Complex elliptically
  symmetric distributions: survey, new results and applications,'' \emph{{IEEE}
  Trans. Signal Process.}, vol.~60, no.~11, pp. 5597--5625, 2012.

\bibitem{magnus_neudecker:1999}
J.~R. Magnus and H.~Neudecker, \emph{Matrix Differential Calculus with
  Applications in Statistics and Econometrics}, 2nd~ed.\hskip 1em plus 0.5em
  minus 0.4em\relax Chichester: Wiley, 1999, 422 pages.

\bibitem{srivastava2005some}
M.~S. Srivastava, ``Some tests concerning the covariance matrix in high
  dimensional data,'' \emph{Journal of the Japan Statistical Society}, vol.~35,
  no.~2, pp. 251--272, 2005.

\bibitem{croux2002sign}
C.~Croux, E.~Ollila, and H.~Oja, ``Sign and rank covariance matrices:
  statistical properties and application to principal components analysis,'' in
  \emph{Statistical data analysis based on the L1-norm and related
  methods}.\hskip 1em plus 0.5em minus 0.4em\relax Birkh{\"a}user Basel, 2002,
  pp. 257--269.

\bibitem{bentler1986greatest}
P.~M. Bentler and M.~Berkane, ``Greatest lower bound to the elliptical theory
  kurtosis parameter,'' \emph{Biometrika}, vol.~73, no.~1, pp. 240--241, 1986.

\end{thebibliography}
