% !TEX TS-program = LuaLaTeX
% !TEX encoding = UTF-8 Unicode 
% !TEX root = ../../mal-core.tex

\gdef\mujnazevCS{Odhad varianční matice ve vysoké dimenzi}
\gdef\mujnazevEN{Covariance matrix estimation in high-dimensional problems}
\gdef\mujnazevPR{\mujnazevCS}
\gdef\mujnazevDR{\mujnazevEN}
\gdef\mujauthor{Marie Turčičová, Kryštof Eben}

\bonushyper

\nazev{Odhad varianční matice\\[2pt] ve vysoké dimenzi}

\nazev{Covariance matrix estimation in\\[2pt] high-dimensional problems}

\author{Marie Turčičová$^{1}$, Kryštof Eben$^2$}

\Adresa{$^1$KPMS MFF UK, Sokolovská 49/83, 186\,00, Praha~3 -- Karlín \\
$^{1,2}$Ústav informatiky AV ČR, Pod Vodárenskou věží 271/2, 182\,07, Praha~8}

\Email{$^1$turcicova@cs.cas.cz, $^2$eben@cs.cas.cz}

%\enlargethispage{\baselineskip}

\Abstrakt{%
V~řadě statistických aplikací, kde je dimenze náhodného vektoru vysoká v~porovnání s~počtem dostupných měření, je velkým problémem odhad varianční matice. Klasická výběrová varianční matice má v~takovém případě řadu nežádoucích vlastností, zejména nízkou hodnost a~malou spolehlivost odhadu jednotlivých prvků. Tento článek obsahuje přehled metod, které se v~tomto případě k~odhadu varianční matice používají. 
Pozornost je nejdříve věnována výpočetně jednoduchým metodám pracujícím po prvcích, mezi které patří například metoda smrštění (shrinkage), posílení diagonály (tapering) a~další. Dále je uveden přehled složitějších přístupů, které používají parametrické modely založené na různých dodatečných předpokladech o~vlastnostech náhodného vektoru, zejména normality%gaussovskosti
, kovarianční stacionarity nebo markovské vlastnosti. Parametrické modely se používají jak k~popisu poklesu vlastních čísel, tak k~přímému modelování varianční matice či její inverze. Parametry příslušných modelů lze odhadovat standardními statistickými postupy.}
\smallskip
\KlicovaSlova{varianční matice, odhad, vysoká dimenze, regularizace}

\Abstract{In many statistical applications, where the dimension of a~random vector highly exceeds the number of available measurements, the estimation of covariance matrix poses a~challenge. The sample covariance matrix has several undesirable properties in this case, specifically low rank and poor accuracy of estimation of its individual elements. This paper provides an overview of methods that are used for covariance matrix estimation in high-dimensional problems. First, we pay attention to computationally simple methods which usually work element-wise, such as shrinkage, tapering, etc. Further, more complex approaches are presented, which employ parametric models based on additional assumptions about the properties of the random vector, especially normality, covariance stationarity and Markov property. Parametric models are used to describe the decay of eigenvalues or to model the covariance matrix or its inverse. Parameters of the corresponding models can be estimated by standard statistical techniques.}
\smallskip
\KeyWords{covariance matrix, estimator, high-dimension, regularization}


\section{Úvod}\label{sec01}
V~tomto článku se budeme zabývat metodami odhadu varianční matice náhodného vektoru $\mathbf{X}=(X_1,\ldots,X_p)^\top$ délky $p$ založenými na náhodném výběru, jehož rozsah $N$ je výrazně menší než $p$. Nejdůležitější aplikace těchto metod se týkají vícedimenzionálních náhodných polí, zejména v~prostorové statistice, kdy například v~každém z~500\,000 ($=p$) bodů třídimenzionální mřížky o~rozměrech $100 \times 100 \times 50$ zjišťujeme hodnotu určité proměnné, zatímco dostupný náhodný výběr má rozsah v~řádu desítek. Taková situace je typická pro atmosférické vědy, kde 
stav atmosféry popisovaný numerickým modelem má vysokou dimenzi, ale výpočetní náročnost modelu nedovoluje pracovat s~velkým počtem \uv{scénářů} a~tvořit náhodné výběry o~velkém rozsahu, které by dobře popisovaly stochastické vlastnosti modelových polí.   

Kromě prostorové statistiky se problém odhadování varianční matice v~situaci malého výběru objevuje v~mnoha dalších oborech, např. finančnictví, genetice,  lékařské statistice, epidemiologii nebo zpracování obrazu. Jednou z~typických úloh, kde odhad varianční matice hraje klíčovou roli, je začleňování aktuálních pozorování do dynamického systému (tzv. datová asimilace). 

Metody odhadu uvedené v~tomto článku jsou určeny zejména pro tuto aplikaci, tak jak se používá ve vědách o~atmosféře. Příslušné náhodné pole pak obsahuje např. hodnoty různých meteorologických veličin. 
Při modelování variančních matic vícerozměrných polí se obvykle body mřížky uspořádají do vektoru, například u~dvourozměrného pole se zkonstruuje vektor obsahující postupně všechny sloupce matice mřížky. Vektor $\mathbf{X}$ v~dalším textu obvykle značí takto uspořádané náhodné pole. Pak je možné použít  různé metody odhadu včetně standardních,  je však třeba vzít v~úvahu speciální strukturu varianční matice vyplývající z~povahy původního problému. 

V~celém článku budeme předpokládat, že rozdělení vektoru $\mathbf{X}$
má konečné druhé momenty a~označíme $\operatorname{var}(\mathbf{X})\equiv\Sigma$. Standardním odhadem varianční matice $\Sigma$ založeným na náhodném výběru $\mathbf{X}_1,\ldots,\mathbf{X}_N$ je výběrová varianční matice 
\begin{equation}
S~= \frac{1}{N-1} \sum_{k=1}^N (\mathbf{X}_k - \overline{\mathbf{X}})(\mathbf{X}_k - \overline{\mathbf{X}})^\top,
\label{sample_cov}
\end{equation}
kde $\overline{\mathbf{X}} = \frac{1}{N} \sum_{k=1}^N \mathbf{X}_k$. V~naší situaci malého počtu pozorování je to odhad nekvalitní a~trpí celou řadou nedostatků. Je-li rozdělení $\mathbf{X}$ spojité a~$p/N<1$, ale ne dostatečně blízké nule, je sice matice $S$ skoro jistě invertibilní, ale je obvykle špatně podmíněná, což mimo jiné znamená, že její inverze $S^{-1}$ je nekvalitním odhadem matice $\Sigma^{-1}$. V~případě $p>>N$ má již matice $S$  nízkou hodnost a~výpočet její inverze není vůbec možný. Odhad matice přesnosti~$\Sigma^{-1}$ je přitom v~řadě situací klíčový. Ke špatné podmíněnosti matice~$S$ se přidává jako další nežádoucí jev výskyt tzv. rušivých kovariancí (angl. spurious covariances), což jsou nereálně vysoké odhady kovariancí vzniklé pouze v~důsledku malého rozsahu výběru.

V~tomto článku se zaměříme na metody regularizace odhadu varianční matice, kde regularizací rozumíme libovolnou techniku vedoucí k~regulární matici, která bude dobrým odhadem skutečné varianční matice.  
Metody regularizace můžeme dělit na neparametrické a~parametrické, dále rozlišovat podle toho, zda probíhají v~původním (fyzickém) prostoru nebo v~transformovaném prostoru (spektrálním, waveletovém), a~konečně také podle toho, zda je regularizována varianční matice nebo její inverze.  Níže popsané neparametrické metody se často dají analogicky použít i~po transformaci, proto je   v~sekci \ref{Neparametrické metody} uvedeme jen ve fyzickém prostoru. V~sekci  \ref{Parametrické metody} popíšeme různé parametrické modely založené na dodatečných předpokladech o~struktuře odhadované varianční matice, a~to jak ve fyzickém prostoru, tak v~transformovaném prostoru a~pro inverzní varianční matici. 


\section{Neparametrické metody} \label{Neparametrické metody}
\subsection{Odhad metodou smrštění (tzv. shrinkage odhad)}\label{shrinkage}
Je známo, že deficit v~rozsahu výběru má za následek zkreslení struktury vlastních čísel výběrové varianční matice v~tom smyslu, že největší (nej\-men\-ší) vlastní čísla jsou při odhadu nadhodnocena (podhodnocena). Celá řada metod se zabývá korekcí tohoto fenoménu. Příkladem mohou být odhady metodou smrštění, tzv. \textit{shrinkage} odhady. 
Asi nejznámější z~nich je lineární \textit{shrinkage} odhad definovaný v~článku \cite{Ledoit-2004-WEL}, který je dobře podmíněný a~vede ke zpřesnění odhadu varianční matice, aniž by se předpokládala jakákoliv speciální struktura. Tento odhad je tvaru 
\begin{equation}
S_1 = \rho \nu I~+ (1-\rho) S,
\label{eq:shrinkage}
\end{equation}
kde  $0 < \rho < 1$ a~$ \nu > 0$. 
Výsledný odhad tedy vzniká smrštěním výběrové varianční matice $S$ k~diagonální matici. V~článku \cite{Ledoit-2004-WEL} se optimální odhad parametrů $\rho, \nu$ hledá prostřednictvím minimalizace střední kvadratické ztráty 
\begin{align*}
\min_{\nu,\rho} \operatorname{E} \left\Vert \rho \nu I~+ (1-\rho) S~- \Sigma \right\Vert_F^2,
\end{align*}
kde $\left\Vert A \right\Vert_F = \sqrt{ \operatorname{tr}(AA^\top)}$ značí Frobeniovu normu matice $A$. Při výpočtu minimalizace se využije, že $\operatorname{E}S = \Sigma$.
Autoři zjistili, že optimální koeficienty lineární kombinace (\ref{eq:shrinkage}) jsou
\begin{align}
	\nu&= \frac{1}{p} \operatorname{tr}(\Sigma), & \rho &= \frac{\beta^2}{\alpha^2 + \beta^2}= \frac{\beta^2}{\delta^2},
	\label{nu_and_rho}
\end{align}
kde
\begin{align*}
\beta^2& = \frac{1}{p} \operatorname{E}\left\Vert S-\Sigma \right\Vert_F^2, &
\alpha^2& = \frac{1}{p} \left\Vert \Sigma-\nu I~\right\Vert_F^2, & \delta^2& = \frac{1}{p} \operatorname{E}\left\Vert S-\nu I~\right\Vert_F^2,
\end{align*}
které bohužel závisejí na neznámé varianční matici $\Sigma$. V~článku \cite{Ledoit-2004-WEL} 
lze ale nalézt pro parametry $\nu,\ \beta,\ \alpha$ a~$\delta$ konzistentní odhady. Jejich dosazením do (\ref{eq:shrinkage}), s~využitím vztahu (\ref{nu_and_rho}) pro $\rho$, pak získáme pozitivně definitní odhad matice~$\Sigma$. 

Důvod, proč odhad (\ref{eq:shrinkage}) funguje dobře, lze nahlédnout z~toho, že výběrová vlastní čísla (tj. vlastní čísla výběrové varianční matice $S$) jsou více rozptýlena kolem svého průměru než skutečná vlastní čísla $\Sigma$.  V~\cite{Ledoit-2004-WEL} je dokázáno, že použitím konvexní kombinace (\ref{eq:shrinkage}) dojde ke smrštění výběrových vlastních čísel k~jejich průměru, a~tím k~celkovému zlepšení odhadu. Myšlenka vylepšení odhadu varianční matice prostřednictvím smrštění výběrových vlastních čísel k~jejich průměru je podrobněji popsána v~přehledovém článku \cite{Muirhead-DEE-1987}. 

Vzorec (\ref{eq:shrinkage}) může mít obecnější tvar
\begin{equation}
S^*_1 = \rho T + (1-\rho) S,
\label{eq:shrinkage_gen}
\end{equation}
kde $T$ je cílová matice (tzv. \textit{target matrix}), která disponuje kýženými vlastnostmi varianční matice, jako je plná hodnost a~pozitivní definitnost. Cílová matice $T$ se často volí diagonální. 

Místo přímého odhadu koeficientu $\rho$ ve vzorci (\ref{eq:shrinkage_gen}), popř. (\ref{eq:shrinkage}), lze smrštění dosáhnout také tak, že se dostupný náhodný výběr rozšíří o~další simulovaná pozorování, která se nagenerují z~rozdělení s~varianční maticí $T$. Z~takto vzniklého souboru se pak vypočte výběrová varianční matice (\ref{sample_cov}). Tato procedura má podobný efekt jako smrštění původní výběrové varianční matice~$S$ k~cílové matici $T$. Konstrukce nových členů se odvíjí od předpokládaných statistických vlastností náhodného vektoru $\mathbf{X}$.

\subsection{Posílení diagonály matice (tzv. tapering)}
Efektivní a~zároveň jednoduchý způsob jak potlačit rušivé kovariance ve výběrové varianční matici je vynásobit $S$ po prvcích řídkou pozitivně definitní maticí $M$. Označíme-li výsledný odhad jako $S_2$, pak
\begin{equation}
S_2 = S~\circ M,
\label{tapering}
\end{equation}
kde $\circ$ značí Schurův součin. Tato metoda bývá v~anglické literatuře označována jako \textit{tapering} (v~\cite{Pourahmadi} jako \textit{banding}) a~matice $S_2$ jako \textit{tapered matrix}. Díky řídkosti matice $M$ je řídká i~matice $S_2$, což přináší zejména v~datové asimilaci velké výpočetní výhody. 

Je-li $M$ pozitivně definitní matice, pak díky větě o~Schurově součinu (např. \cite{Horn_Johnson}, str. 479, Věta 7.5.3) je zajištěno, že $S_2$ bude také pozitivně definitní. 

V~kontextu datové asimilace byla regularizace kovariance pomocí Schurova součinu navržena v~\cite{Houtekamer-2001-SEK} a~\cite{Hamill-2001-DDF} a~vychází z~principů modelování kovariance náhodného pole pomocí kovarianční funkce. Nechť tedy $\mathcal{X}(\boldsymbol{s}) $ je spojité náhodné pole definované na omezené doméně $\mathfrak{D} \subset \mathbb{R}^3$ pokrývající část atmosféry Země
a~$\mathbf{X}=(X_1,\ldots,X_p)^\top$ představuje nějakou diskretizaci  $\mathcal{X}$, tj.  $X_i = \mathcal{X}(\boldsymbol{s}_i), \; \boldsymbol{s}_i \in \mathfrak{D}$.  Funkce  $C(\boldsymbol{s},  \boldsymbol{t}) = \mathtt{cov} (\mathcal{X}(\boldsymbol{s}),  \mathcal{X}(\boldsymbol{t} ))$ se nazývá kovarianční funkce pole $ \mathcal{X} $. Platí, že pro každou diskretizaci  $\mathbf{X} $ je matice  s~prvky   $C(\boldsymbol{s_i},  \boldsymbol{s_j}) $
pozitivně semidefinitní a~naopak, je-li pro libovolný výběr bodů $\boldsymbol{s_i}$ kovarianční matice s~prvky   $C(\boldsymbol{s_i},  \boldsymbol{s_j}) $  pozitivně semidefinitní, definuje $C$ kovarianční funkci.

Matice $M$ je v~\cite{Houtekamer-2001-SEK} konstruována tak, že je pozitivně definitní a~její $(i,j)$-tý prvek $M_{ij}$ je roven $\varrho(\left\Vert \boldsymbol{s}_i - \boldsymbol{s}_j \right\Vert_{\mathbb{R}^3})$, kde $\varrho \colon [0, \infty) \to [0,1]$ je funkce s~kompaktním nosičem a~$\left\Vert \cdot \right\Vert_{\mathbb{R}^3}$ je norma v~prostoru $\mathbb{R}^3$. %dán korelační funkcí $\varrho$ s kompaktním nosičem, která monotónně klesá se vzdáleností mezi $X_i$ a $X_j$. 
Jinými slovy,  $\varrho$ je možné ztotožnit s~dobře definovanou korelační funkcí a~(\ref{tapering}) tedy modeluje kovarianci pole jako Schurův součin výběrové kovariance s~korelační maticí danou funkcí   $\varrho$.
Zkonstruovat takovou funkci (zejména s~kompaktním nosičem) je netriviální problém, kterým se zabývá článek \cite{Gaspari-Cohn}, z~něhož je v~\cite{Hamill-2001-DDF} vzata funkce $\varrho$ hladká a~monotónně klesající k~nule, konkrétně jde o~polynom pátého řádu \cite[vztah~(4.10)]{Gaspari-Cohn}. Graf tvarem připomíná pravou polovinu Gaussovy křivky, která však nabývá nulové hodnoty v~konečné vzdálenosti, což zajišťuje řídkost~$M$.

Již pod názvem \textit{tapering} se tato regularizační metoda objevuje v~článku~\cite{Furrer-2007-EHP}, kde se matice $M$ určuje z~minimalizace vzdálenosti měřené střední čtvercovou chybou
\begin{equation} 
\mathrm{MSE}(S_2)= \operatorname{E} \left\Vert \Sigma - S_2 \right\Vert_F^2=\operatorname{E} \left( \operatorname{tr} \left( (\Sigma - S~\circ M )^2 \right) \right).
\end{equation}
Minimalizace by se měla provádět přes třídu pozitivně definitních matic $M$, což je výpočetně náročné. Toto omezení se proto ignoruje a~výsledná matice~$M$, pro jejíž prvky lze při daném $\Sigma$ získat explicitní formuli, se dodatečně převede na pozitivně definitní matici dalšími heuristickými postupy (např. se v~jejím spektrálním rozkladu ponechají pouze kladná vlastní čísla a~ostatní se položí rovna nějakému $\varepsilon >0$).

Druhou možností je pak v~souladu s~články \cite{Houtekamer-2001-SEK} a~\cite{Hamill-2001-DDF} matici $M$ parametrizovat pomocí validní kovarianční funkce, která bude popisovat dosah korelací, a~prostřednictvím minimalizace MSE odhadovat její paramet\-ry. Tam je již pozitivní definitnost zaručena, řídkost ale zajištěna  není, což dělá tuto metodu poněkud méně výpočetně atraktivní. 

\subsection{Vynulování zanedbatelných prvků (tzv. tresholding)}
Pro úplnost na závěr zmiňme metodu, pro niž se v~anglické odborné literatuře používá termín \textit{tresholding}. Myšlenka vychází z~toho, že vysokodimenzionální varianční matice, které se vyskytují v~praktických aplikacích, jsou obvykle řídké, a~obsahují tudíž mnoho nulových prvků. Vylepšený odhad, navržený původně v~\cite{Bickel-2008}, pak je
\begin{equation}
T_t(S) \equiv \left( s_{ij} \mathbf{1}_{[|s_{ij}| \geq t]}: i,j=1,\ldots,p \right),
\end{equation}
kde $T_t(S)$ značí operátor tresholdingu aplikovaný na výběrovou varianční matici a~$t > 0$ je zvolený práh, pod nímž již prvky pokládáme za zanedbatelné. Autoři \cite{Bickel-2008} doporučují $t$ zvolit podle následujícího postupu: Výběr se náhodně rozdělí na dvě části o~velikostech $N_1 = N \left( 1- \frac{1}{\log N} \right)$ a~$N_2 = \frac{N}{\log N}$ a~vypočtou se výběrové varianční matice $S_{N_1}$ a~$S_{N_2}$. Tento postup se opakuje $K$-krát  a~$t$ se zvolí tak, aby minimalizovalo
\begin{equation}
R(t)=\frac{1}{K} \sum_{k=1}^K \left\Vert T_t(S_{N_1,k})-S_{N_2,k}\right\Vert_F^2.
\end{equation}

Velkou předností metod jako tapering a~tresholding je jednoduchá implementace. Nevýhodou však může být to, že pozitivní definitnost výsledku často není zajištěna. Pro tresholding je v~\cite{Bickel-2008} alespoň ukázána stejnoměrná konzistence v~operátorové normě na množině řídkých matic, a~to za podmínky, že $\frac{\log p}{N} \to 0.$

\section{Parametrické metody} \label{Parametrické metody}
Pokud náhodný vektor $\mathbf{X}$ má specifické vlastnosti, které se promítají do tvaru jeho varianční matice, lze jejich zohledněním ve formě předpokladů odhad výrazně vylepšit. 
Například u~meteorologických veličin je často možné předpokládat normalitu, kovarianční stacionaritu  nebo prostorovou markovskou vlastnost. Každá z~těchto vlastností poskytuje možnost vylepšení odhadu varianční matice pomocí specifického parametrického modelu, jehož parametry lze odhadnout standardními statistickými metodami. Přesnost výsledného odhadu a~jeho kvalita, která se projeví v~dalším použití (např. v~datové asimilaci), se pak odvíjí od toho, jak realistické tyto dodatečné předpoklady byly. 

U~parametrických metod je nezbytné vyhnout se příliš komplexním modelům s~velkým počtem parametrů. Proto řada přístupů využívá různé transformace, které vedou k~přibližné diagonalitě varianční matice v~transformovaném prostoru. 
Tím dochází k~velké redukci počtu parametrů. 


\subsection{Regularizace ve spektrálním a~waveletovém prostoru}
\label{reg_spektr_wavelet}

Tento přístup je založen na reprezentaci vektoru $\mathbf{X}$ rozvojem 
\begin{equation}
\mathbf{X} = \operatorname{E}\mathbf{X} +   \sum_{i=1}^p d_i^{1/2} \xi_i \mathbf{v}_i ,
\label{eq:expansion}
\end{equation}
kde $d_i$ jsou koeficienty,  $\xi_i$ jsou nezávislé náhodné veličiny s~nulovou střední hodnotou a~jednotkovým rozptylem a~$ \mathbf{v_i} $ jsou ortonormální vektory v~$\mathbb{R}^p$. Pro varianč\-ní matici pak máme 
\begin{equation}
    \label{eq:FDF}
    \Sigma  = FDF^\top ,  \quad \text{ kde }D = \text{diag} (d_1, \ldots , d_p)
\end{equation}
a~matice $F$ má ve sloupcích vektory  $\mathbf{v}_i$. Je-li (\ref{eq:FDF}) spektrální rozklad varianční matice, dostáváme známý rozklad do hlavních komponent (Karhunen-Lo\`{e}ve). Ve velké dimenzi je ale většina výběrových vlastních čísel nulová a~odhad teoretického rozvoje (\ref{eq:expansion}) je obtížný. Pro regularizaci se proto využívají vhodné deterministicky dané báze (nebo obecněji i~framy jako v~případě waveletové transformace).

Modelujeme-li pole $\mathbf{X}$ pomocí (\ref{eq:expansion}),  můžeme  matici $\Sigma$ odhadnout prostřednictvím odhadu diagonální matice $D = F^\top \Sigma F$. Pokud nepřijmeme žádné další předpoklady, je možné podobně jako v~předchozí sekci použít některou neparametrickou metodu na výběrovou varianční matici $S$ pole $\mathbf{X}$ převedenou do spektrálního prostoru, tj. na $F^\top S F$. Tato matice přirozeně není přesně diagonální,  ale obsahuje (v~praxi často malé) nenulové mimodiagonální členy.  Nejjednodušším způsobem regularizace odhadu je pak jejich zanedbání \cite{Kasanicky-2015}. Jedná se vlastně o~\textit{tapering} spektrální varianční matice podle vzorce (\ref{tapering}) pro $M= I$. 

Za předpokladu konkrétního rozdělení vektoru $\mathbf{X}$ lze diagonální prvky $\{d_{ii} \}_{i=1}^p$ spektrální varianční matice $D$ odhadnout, například metodou maximální věrohodnosti založenou na náhodném výběru $\mathbf{X}_1,\ldots,\mathbf{X}_N$. V~případě normálního rozdělení dostáváme \cite{Turcicova-2019}
\begin{equation}
\widehat{d}_{ii} = \frac{1}{N} \sum_{k=1}^N X_{ki}^2, \quad i=1,\ldots,p,
\end{equation}
kde $X_{ki}$ značí $i$-tý prvek vektoru $\mathbf{X}_k$, $k=1,\ldots,N $.  

\subsubsection{Regularizace ve spektrálním prostoru}
\label{reg_spektr}

Doposud jsme nepředpokládali nic o~výběru ortogonální matice $F$. Podle očekávání častou volbou pro $F$ bude Fourierova báze. V~meteorologii podobně jako ve fyzice obecně se často vyskytují děje difúzního typu a~je přirozené, že se modely popisující difúzi úspěšně používají i~pro popis stochastických vlastností některých fyzikálních veličin. Tak se ukazuje souvislost mezi stacionaritou pole, Fourierovou transformací a~Laplaceovou rovnicí, a~je možné založit na ní regularizační techniky.  Pokud např.  mů\-že\-me považovat vektor  $\mathbf{X}$ za pozorování omezeného časového úseku kovariančně stacionárního procesu $\{ X_t ,\, t \in \mathbb{Z}\}$ pro $t=1,\ldots, p $,  vede diskrétní Fourierova transformace vektoru  $\mathbf{X}$  k~přibližné diagonalitě\footnote{Matice $\Sigma$ je v~transformovaném prostoru přesně diagonální pouze v~případě, kdy je proces $\{ X_t ,\, t \in \mathbb{Z}\}$ periodický s~periodou $p$, tj. $X_t = X_{t+p}, \; \forall t \in \mathbb{Z}$.} varianční matice v~transformovaném prostoru \cite{Dwivedi-2011}, tj. můžeme předpokládat model
\begin{equation}
    \label{eq:Fourier}
    \Sigma  = FDF^\top,
\end{equation}
kde matice $F$ provádí Fourierovu transformaci. 

V~meteorologických aplikacích vektor $\mathbf{X}$ často představuje diskretizaci spojité omezené domény, tedy oblasti pokrývající část atmosféry Země, a~index má význam souřadnice v~prostoru. Předpokládejme pro jednoduchost, že spojitou doménou je úsečka $[0,1]$ a~vektor $\mathbf{X}$ délky $p$ ji rozděluje na $p+1$ dílků délky $h = \frac{1}{p+1}$. Nyní můžeme využít faktu, že matice $F$ je tvořena vlastními vektory diskrétního Laplaceova operátoru $L \colon \mathbb{R}^p \rightarrow \mathbb{R}^p$, který  je v~jedné dimenzi identický s~operátorem druhé derivace  a~lze ho reprezentovat maticí  \begin{equation}
L = \frac{1}{h^2} \begin{pmatrix}
-2 & 1 & 0 & \ldots & 0 \\
1 & -2 & 1 & \ldots & 0 \\
& & \ddots & & \\
0 & \ldots & 1 & -2 & 1 \\
0 & \ldots & 0 & 1 & -2
\end{pmatrix}.
\label{def:diskr_Laplace}
\end{equation}
Pak pro $\mathbf{u}=(u_1,\ldots,u_p)^\top \in \mathbb{R}^p$ je $j$-tá složka vektoru $L(\mathbf{u}) \equiv L \mathbf{u} $ rovna
\begin{equation*}
L(\mathbf{u})_j = \frac{u_{j-1} - 2 u_j + u_{j+1}}{h^2},
\end{equation*}
kde $j=1,\ldots,p$. Pro body $u_0$ a~$u_{p+1}$ (ležící na krajích úsečky $[0,1]$) se uvažují různé okrajové podmínky. Definice (\ref{def:diskr_Laplace}) odpovídá tzv. Dirichletově okrajové podmínce, při níž je $u_0 = u_{p+1}=0$. Vlastní vektory  $\mathbf{v}_i$ mají v~tomto případě tvar 

\begin{equation}
\mathbf{v}_i =  \sqrt{\frac{2}{p+1}} \begin{pmatrix}
\sin\big( \frac{1}{p+1} i\pi \big) \\[5pt]
\sin\big( \frac{2}{p+1} i\pi \big) \\
\vdots \\
\sin\big( \frac{p}{p+1} i\pi \big)
\end{pmatrix}, \quad i=1,\ldots,p
\end{equation}
a~dostáváme (jak je možno se přesvědčit) spektrální rozklad  $ L = F \Lambda F^\top$, kde  $\Lambda$ je diagonální matice s~vlastními čísly
\begin{equation}
\lambda_i = - 4(p+1)^2 \sin^2 \left( \frac{\pi}{2 (p+1)}i \right),  \quad i=1,\ldots,p
\label{eigval_Laplace}
\end{equation}
na diagonále. 

Máme tedy $ L = F \Lambda F^\top$ a~$ \Sigma = F D F^\top$ a~je přirozené modelovat $ \Sigma $ 
jako funkci Laplaceova operátoru. 
Je-li totiž $ L = F \Lambda F^\top$,  pak pro spojitou funkci~$f$ lze definovat matici $f(L) $ pomocí spektrálního rozkladu $Ff(\Lambda)F^\top$, kde $f(\Lambda) = \operatorname{diag}(f(\lambda_{1}),f(\lambda_{2}),\ldots,f(\lambda_{p}))$. 

Vlastní čísla $d_{ii}$ můžeme modelovat jako vhodnou funkci vlastních čísel (\ref{eigval_Laplace}) Laplaceova operátoru, tj. $d_{ii} = f(\lambda_i)$,  $i=1,\ldots,p$. 
Jelikož $\mathbf{X}$ představuje diskretizaci spojitého náhodného pole $\mathcal{X}$ na kompaktní oblasti $\mathfrak{D}$, lze varian\-ční matici $\Sigma$ pokládat za diskretizaci variančního operátoru $\mathcal{C}$ pole $\mathcal{X}$. 
Za předpokladu spojitosti kovarianční funkce $\operatorname{cov}(\mathcal{X}(\boldsymbol{s}),\mathcal{X}(\boldsymbol{t}))$, $\boldsymbol{s},\boldsymbol{t} \in \mathfrak{D}$ pole $\mathcal{X}$ lze pomocí metod známých z~numerické matematiky ukázat, že vlastní čísla $\{d_{11},\ldots,d_{pp} \}$ matice  $\Sigma$ konvergují pro $p \to \infty$ k~vlastním číslům variančního operátoru $\mathcal{C}$ \cite{Atkinson-1975,Osborn-1975}.  
Jelikož každý varianční operátor musí mít konečnou stopu (\cite{Kuo}, Věta 2.1), 
musejí $\{d_{ii} \}_{i=1}^p$ s~rostoucím $i$ rychle klesat k~nule i~pro $p$ konečné.  

Předpokládejme tedy $\Sigma = Ff(\Lambda)F^\top$ pro nějakou vhodnou $f$, a~tudíž $d_{ii} = f(\lambda_i)$.
Kvůli požadavku na konečnost stopy varianční matice musí funkce $f$ pro $ \lambda \rightarrow  -\infty$ rychle klesat k~nule, neboť $\{\lambda_i\}_{i=1}^p$ je klesající posloupnost záporných čísel. 

V~praxi se užívá zejména polynomiální a~exponenciální model:
\begin{align*}
d_{ii} &=c (-\lambda_i)^{-\alpha}, \quad \alpha >1 \quad \text{(polynomiální model)} \\
d_{ii} &= c \ee^{\alpha \lambda_i} \quad \text{(exponenciální model)},
\end{align*}
kde $i=1,\ldots,p$. Koeficienty $c,\,\alpha$ jsou parametry daného modelu. To představuje v~daném případě snížení počtu parametrů z~$p$ na 2, přesto i~takto jednoduchý model může být přijatelný. Za předpokladu konkrétního rozdělení~$\mathbf{X}$ lze pak parametry odhadnout například metodou maximální věrohodnosti~\cite{Turcicova-2019}. Použití takového modelu koriguje strukturu vlastních čísel podobně jako shrinkage (viz začátek sekce \ref{shrinkage}), vyhlazuje také průběh odhadů $\widehat{d}_{ii}$ a~tím přispívá k~redukci šumu. 

Náš příklad byl pouze ilustrativní. Zajímavé a~hluboké vztahy mezi určitým typem stochastické rovnice difúze, 
náhodnými poli s~(prostorovou) markovskou vlastností a~kovariancemi z~tzv. Matérnovy třídy popisuje článek  \cite{Lindgren-2011-ELG}
a~též přehledný článek \cite{Simpson-2012}. Použití těchto metod v~praxi pak ilustruje  \cite{Mirouze-RCF-2010}. 

\subsubsection{Regularizace ve waveletovém prostoru}
Reprezentuje-li $\mathbf{X}$ náhodný vektor či náhodné pole v~prostoru jako je tomu u~meteorologických veličin, je zřejmé, že takové pole může mít jak charakteristiky vlnové povahy tak i~lokální vlastnosti vázané na umístění v~prostoru. Fourierova transfor\-ma\-ce postihne frekvenční komponenty vektoru $\mathbf{X}$, ale nezachytí lokální jevy. Nacházíme zde podobně jako v~řadě dalších aplikací určitou komplementaritu spektrální a~prostorové reprezentace. 

Nástrojem ke \uv{kompromisní} reprezentaci takového náhodného pole je tzv. waveletová transformace \cite{Burrus}, která je realizována podobně jako u~Fourierovy transformace rozvojem typu  (\ref{eq:expansion}) a~rozkladem varianční matice (\ref{eq:FDF}). Matice transformace se skládá z~bázových funkcí, tvořených různě umístěnými a~škálovanými vlnkami (odtud název \textit{wavelets}). 
Obecněji je možno v~rozvoji~(\ref{eq:expansion}) použít místo ortonormální báze tzv. frame, který  má některé vlastnosti báze, ale může mít i~více než $p$~členů. Různé volby framů pak dávají waveletové reprezentace různých typů.  

Po transformaci varianční matice do waveletového prostoru lze pak zanedbáním příslušných koeficientů vyřadit bázové funkce s~nízkou vypovídací schopností. Často je možné podobně jako u~Fourierovy transformace zanedbat mimodiagonální prvky transformované matice, což odpovídá lokálnímu průměrování prostorových korelací \cite{Pannekoucke-2007}. Můžeme také předpokládat vhodnou řídkou strukturu kovarianční matice ve waveletovém prostoru. 
 Po transformaci zpět do původního prostoru pak dostáváme varianční matici očištěnou o~šum. 
 
Reprezentace ve waveletovém prostoru často přináší dramatické snížení dimenze problému. 
Koeficienty bázových funkcí lze pak odhadovat pomocí vhodného parametrického nebo i~semiparametrického modelu. 
 
\subsubsection{Využití dalších typů rozkladu varianční matice k~jejímu odhadu}
Kromě spektrální a~waveletové transformace lze ke konstrukci pozitivně definitního odhadu varianční matice využít i~jiné typy rozkladu. 
Jednoduchý rozklad
\begin{equation}
\Sigma = BRB
\end{equation}
na matici směrodatných odchylek $B$ a~korelační matici $R$ má velký praktický význam, neboť oba faktory rozkladu jsou snadno interpretovatelné. To umožňuje odhadovat $B$ i~$R$ zvlášť, což může mít význam zejména v~aplikacích, kde je jeden z~těchto prvků důležitější. 

Druhým typem rozkladu, který se při odhadování s~úspěchem používá, je Choleského rozklad
\begin{equation}
\Sigma = C C^\top,
\end{equation}
kde $C$ je (pro pozitivně definitní $\Sigma$) jednoznačně daná dolní trojúhelníková matice s~pozitivními prvky na diagonále. V~určitých případech může být výhodnější modifikovaný Choleského rozklad
\begin{equation}
\Sigma = L D^2 L^\top,
\label{cholesky_2}
\end{equation}
kde $L=CD^{-1}$ je dolní trojúhelníková matice a~$D^2$ je diagonální matice, kterou lze opět modelovat zvlášť. 

Při těchto typech odhadu lze s~výhodou využít znalost (prostorových nebo časových) vazeb mezi složkami vektoru  $\mathbf{X}$. Zafixováním veličin ve vhodném pořadí lze totiž docílit pásové struktury $\Sigma$ a~následně řídkosti některých členů ve zvoleném rozkladu. 


\subsection{Lineární model kovariance}
Výše popsané metody nepředpokládaly žádnou speciální strukturu varianční matice nebo její inverze a~zabývaly se převážně regularizací výběrové kovarianční matice v~původním nebo transformovaném  prostoru. Alternativou k~těmto metodám je do jisté míry komplementární přístup, který vychází ze znalosti speciální struktury varianční matice nebo její inverze. Taková struktura může být dána prostorovými nebo i~frekvenčními charakteristikami výchozího náhodného pole. Často lze strukturu varianční matice dobře postihnout lineárním modelem \cite{Anderson-1973}  
\begin{equation}
\Sigma = \alpha_1 U_1 + \ldots + \alpha_q U_q, 
\label{lin_cov_mod}
\end{equation}
kde $U_i$ jsou známé, lineárně nezávislé (v~prostoru matic), symetrické matice a~$\alpha_i$ jsou neznámé parametry takové, že výsledná matice je pozitivně definitní. Tento model je zcela obecný, neboť pro  $q=\frac{p(p+1)}{2}$ lze každou varianční matici zapsat ve tvaru (\ref{lin_cov_mod}) s~$\alpha_{ij}= \sigma_{ij}=\sigma_{ji}$ a~maticemi $U_{ij}$, které jsou nulové až na jednotky na pozicích $(i,j)$ a~$(j,i)$. 
Teoreticky tak pokrývá libovolnou odhadovací metodu pracující po prvcích, jako je např. \textit{tapering} nebo \textit{banding}. Podrobnější přehled metod vycházejících z~lineárního modelu lze  najít v~\cite{Pourahmadi}.

Parametry  $\boldsymbol{\alpha}=(\alpha_1,\ldots,\alpha_q)^\top$ je možné odhadovat standardními statistickými metodami, jejich použití ale značně komplikuje podmínka, že výsledná lineární kombinace musí být pozitivně definitní matice.
Nutné a~postačující podmínky pro existenci explicitního maximálně věrohodného odhadu jsou k~dispozici v~článku \cite{Szatrowski-1980}. K~výpočtu řešení věrohodnostní rovnice je přitom použita iterativní metoda navržená v~\cite{Anderson-1973}. 

Jednou z~možností jak zajistit pozitivní definitnost odhadu je použít lineární model na logaritmus varianční matice
\begin{equation}
\log \Sigma = \alpha_1 U_1 + \ldots + \alpha_q U_q ,
\label{log_lin_cov_mod}
\end{equation}
kde $U_i$ jsou opět známé matice, ale na $\alpha_i$ nejsou kladeny žádné další požadavky. Nevýhodou tohoto přístupu je to, že $\log \Sigma$ nemá žádnou statistickou interpretaci.
Autoři \cite{Chiu-1996} uvádějí maximálně věrohodný odhad pro parametry modelu  (\ref{log_lin_cov_mod}), a~to včetně jeho asymptotických vlastností. 


\subsection{Regularizace v~inverzním prostoru}
Motivací pro regularizaci varianční matice skrze její inverzi, tzv. matici přesnosti, je fakt, že u~normálně rozdělených vektorů $\mathbf{X}$ má inverze korelační matice na $(i,j)$-tém místě prvek úměrný parciálnímu korelačnímu koeficientu mezi $X_i$ a~$X_j$ při vyloučení všech ostatních prvků. 

Má-li tedy vektor $\mathbf{X}$ normální rozdělení, pak  $(i,j)$-tý prvek matice přesnosti je nulový právě tehdy, když $X_i$ a~$X_j$ jsou nezávislé podmíněně na všech zbývajících prvcích \cite{Rue_Held}. Tento výsledek lze s~výhodou využít u~polí, která jsou navíc tzv. prostorově markovská. Vektor $\mathbf{X}$ má prostorovou markovskou vlastnost, pokud pro kaž\-dé $i=1,\ldots,p$ je podmíněné rozdělení $X_i$ závislé pouze na bodech z~jeho okolí. 
Na ostatních bodech (které jsou daným okolím od bodu $X_i$ odděleny) je $X_i$ podmíněně nezávislé. Odtud vyplývá, že matice přesnosti gaussovského markovského pole je řídká pásová matice. Toho je možno využít pro regularizaci.
Obrázek \ref{fig:structure} ilustruje nenulové prvky matice přesnosti dvourozměrného markovského pole $\mathbf{X}$, jehož každý bod závisí pouze na 8 nejbližších sousedech. 
    \begin{figure}[!h]
    \centering 
    \begin{tikzpicture}[remember picture]
    %\begin{minipage}[c]{9cm}
    \node (hlavni) {\includegraphics[width=0.7\textwidth]{fig1_inv_cov_structure9-crop.pdf}};
     %\end{minipage}
    %\begin{minipage}[c]{2cm}
    \node[anchor=west, xshift=1mm, yshift=0mm] at (hlavni.east) {\includegraphics[width=3cm]{fig1_8neighbours.pdf}};
    %\end{minipage}
    \end{tikzpicture}
    \caption{Matice přesnosti dvourozměrného markovského pole o~rozměru $10 \times 10$ bodů (body jsou uspořádány do sloupce), v~němž kaž\-dý bod závisí na 8 nejbližších sousedech. Vpravo 
    %V %horním 
    %rohu 
    je 
    schéma bodu a~jeho 8 nejbližších sousedů.}
    \label{fig:structure}
    \end{figure}

Metody použité v~\cite{Lindgren-2011-ELG}, \cite{Mirouze-RCF-2010} a~\cite{Simpson-2012}, o~kterých jsme se zmiňovali v~odstavci \ref{reg_spektr} v~souvislosti se stochastickou rovnicí difúze, se týkají gaussovských polí s~markovskou vlastností a~spadají tedy svojí povahou do modelování v~inverzním prostoru. Níže uvedeme ještě další použitelné techniky. 



\subsubsection{Lineární model pro matici přesnosti}
%Ueno
Lineární model lze použít jak pro varianční matici tak pro její inverzi, což je pro gaussovské markovské pole velmi výhodné. Uvažujme model 
\begin{equation}
\Sigma^{-1} = \beta_1 A_1 + \ldots + \beta_r A_r,
\end{equation} 
kde matice $A_1,\ldots,A_r$ jsou známé řídké matice a~$\boldsymbol{\beta}=(\beta_1,\ldots,\beta_r)^\top$ jsou neznámé parametry.  
Vhodnou volbou pro matice $A_1,\ldots,A_r$ mohou být opět matice $A_{ij}$, které mají na místě $(i,j)$ a~$(j,i)$ jedničku a~jinde nuly.  Obdobně je možno použít symetrické matice obsahující nenulové hodnoty pouze na vybraných subdiagonálách, nebo jejich úsecích. 

Parametry $\boldsymbol{\beta}$ lze odhadnout například metodou maximální věrohodnosti \cite{Ueno}.  Výsledný odhad ale není dán explicitně a~je výsledkem numerické maximalizace věrohodnosti. Explicitní vzorec  pro odhad parametrů (viz \cite{smef7}) lze získat metodou tzv. score-matchingu \cite{Hyvarinen}, která, stejně jako metoda maximální věrohodnosti,  poskytuje za určitých předpokladů konzistentní odhad.


\section{Závěr}
Pokrok současné měřicí a~výpočetní techniky vede často k~situaci, kdy rozsah naměřených či modelových dat je velký, ale možnost opakování experimentu je velmi omezená. Problematika odhadu varianční matice na základě náhodného výběru, jehož rozsah je malý v~porovnání s~dimenzí jednotlivých vektorů, se tak objevuje v~mnoha statistických aplikacích. 
V~tomto článku jsme se zaměřili zejména na metody odhadu používané při datové asimilaci, pro niž je kvalitní odhad varianční matice klíčovým prvkem ovlivňujícím kvalitu následné predikce, například při numerické předpovědi počasí.  
Metody jsou však vyvíjeny  v~souvislosti s~mnoha dalšími obory,  často souběžně v~různých vědeckých komunitách,  a~je zde řada otevřených problémů.  

Výčet metod uvedených v~tomto článku není jistě vyčerpávající, ale doufáme, že představuje přehled základních přístupů a~směrů výzkumu v~dané oblasti. 


%\bigskip

\Podekovani{%
Vznik článku byl částečně podporován z~grantu TA~ČR č. 
TL01000238.
Děkujeme Pavlu Krčovi, který zachránil text článku před jeho kompletním ztracením. 
%
}

%\bigskip
%\newpage

\selectlanguage{english}
\renewcommand{\refname}{Literatura}
\begin{thebibliography}{99}
\setlength\itemsep{-1pt}

\bibitem{Anderson-1973} 
Anderson, T.\,W. (1973): Asymptotically Efficient Estimation of
Covariance Matrices with Linear Structure. {\it Ann. Stat.} {\bf 1}(1), 135--141,\\
doi: \href{https://doi.org/10.1214/aos/1193342389}{10.1214/aos/1193342389}.

\bibitem{Atkinson-1975} 
Atkinson, K. (1975): Convergence Rates for Approximate Eigenvalues of Compact Integral Operators. {\it SIAM J. Numer. Anal.} {\bf 12}(2), 213--222.

\bibitem{Bickel-2008}
Bickel, P.\,J., Levina, E. (2008): Covariance Regularization by Thresholding. {\it Ann. Stat.} {\bf 36}(6), 2577--2604,
doi: \href{https://doi.org/10.1214/08-AOS600}{10.1214/08-AOS600}.

\bibitem{Burrus} 
Burrus, C.\,S., Gopinath, R.\,A.,  Guo, H. (1998): {\it Introduction to Wavelets and Wavelet Transforms: A~Primer}. Prentice Hall, New Jersey.

\bibitem{Chiu-1996} 
Chiu, T.\,Y.\,M., Leonard, T., Tsui, K.\,W. (1996): The Matrix-Logarithmic Covariance Model. {\it J. Amer. Statist. Assoc.} 
{\bf 91}, 198--210,\\
doi: \href{https://doi.org/10.2307/2291396}{10.2307/2291396}.

\bibitem{Dwivedi-2011} 
Dwivedi, Y., Rao, S.\,S. (2011): A~Test for Second-Order Stationarity of a~Time Series Based on the Discrete Fourier Transform. {\it J. Time Ser. Anal.} {\bf 32}, 68--91,
doi: \href{https://doi.org/10.1111/j.1467-9892.2010.00685.x}{10.1111/j.1467-9892.2010.00685.x}.

\bibitem{Furrer-2007-EHP} 
Furrer, R., Bengtsson, T. (2007): Estimation of High-Dimensional Prior and Posterior Covariance Matrices in Kalman Filter Variants. {\it J. Multivar. Anal.} {\bf 98}(2), 227--255,
doi: \href{https://doi.org/10.5167/uzh-21542}{10.5167/uzh-21542}.

\bibitem{Gaspari-Cohn}
Gaspari, G., Cohn, S.\,E. (1999): Construction of Correlation Functions in Two and Three Dimensions. {\it Q. J. R. Meteorol. Soc.} {\bf 125}(554), 723--757, 
doi: \href{https://doi.org/10.1002/qj.49712555417}{10.1002/qj.49712555417}.

\bibitem{Hamill-2001-DDF} 
Hamill, T.\,M., Whitaker, J.S., Snyder, Ch. (2001): Distance-Dependent Filtering of Background Error Covariance Estimates in an Ensemble Kalman Filter. {\it Mon. Weather Rev.} {\bf 129}(11), 2776--2790, \\
doi: \href{https://doi.org/10.1175/1520-0493(2001)129<2776:DDFOBE>2.0.CO;2}{10.1175/1520-0493(2001)129{\tt<}2776:DDFOBE{\tt>}2.0.CO;2}.

\bibitem{Horn_Johnson} 
Horn, R.\,A., Johnson, Ch.\,R. (2013): {\it Matrix Analysis}. Cambridge University Press, Second Edition.

\bibitem{Houtekamer-2001-SEK} 
Houtekamer, P.\,L., Mitchell, H.\,L. (2001): A~Sequential Ensemble Kal\-man Filter for Atmospheric Data Assimilation. {\it Mon. Weather Rev.} {\bf 129}(1), 123--137, \\
doi: \href{https://doi.org/10.1175/1520-0493(2001)129<0123:ASEKFF>2.0.CO;2}{10.1175/1520-0493(2001)129{\tt<}0123:ASEKFF{\tt>}2.0.CO;2}

\bibitem{Hyvarinen}
Hyvärinen, A. (2005): Estimation of Non-Normalized Statistical Models by Score Matching. {\it J. Mach. Learn. Res.} {\bf 6}(24), 695--709,\\
doi: \href{https://doi.org/10.5194/npg-22-485-2015}{10.5194/npg-22-485-2015}.

\bibitem{Kasanicky-2015} 
Kasanický I., Mandel J., Vejmelka M. (2015): Spectral Diagonal Ensemble
Kalman Filters. {\it Nonlinear Process. Geophys.} {\bf 2}, 115--143,\\
doi: \href{https://doi.org/10.5194/npg-22-485-2015}{10.5194/npg-22-485-2015}.

\bibitem{Kuo} 
Kuo, H.\,H. (1975):. {\it Gaussian Measures in Banach Spaces}. Lecture Notes in Mathematics. Vol.~463, Springer-Verlag, Berlin,\\
doi: \href{https://doi.org/10.1007/BFb0082007}{10.1007/BFb0082007}.

\bibitem{Lindgren-2011-ELG} 
Lindgren, F., Rue, H., Lindström, J. (2011): An Explicit Link Between Gaussian Fields and Gaussian Markov Random Fields: The Stochastic Partial Differential Equation Approach. {\it J. R. Stat. Soc.
Ser. B Stat. Methodol.} {\bf 73}(4),  423--498,\\
doi: \href{https://doi.org/10.1111/j.1467-9868.2011.00777.x}{10.1111/j.1467-9868.2011.00777.x}.

\bibitem{Ledoit-2004-WEL} 
Ledoit, O., Wolf, M. (2004): A~Well-Conditioned Estimator for Large-Dimensional Covariance Matrices. {\it J. Multivar. Anal.} {\bf 88}(2), 365--411,\\
doi: \href{https://doi.org/10.1016/S0047-259X(03)00096-4}{10.1016/S0047-259X(03)00096-4}.

\bibitem{Mirouze-RCF-2010} 
Mirouze, I., Weaver, A.\,T. (2010): Representation of Correlation Functions in Variational Assimilation Using
an Implicit Diffusion Operator. {\it  Q. J. R. Meteorol. Soc.} {\bf 136}, 1421--1443,
doi: \href{https://doi.org/10.1002/qj.643}{10.1002/qj.643}.

\bibitem{Muirhead-DEE-1987} 
Muirhead, R.\,J. (1987): Developments in Eigenvalue Estimation. {\it Adv. Multivariate Statist. Anal.} 277--288,\\
doi: \href{https://doi.org/10.1007/978-94-017-0653-7_14}{10.1007/978-94-017-0653-7\_14}.

\bibitem{Osborn-1975} 
Osborn, J.\,E. (1975): Spectral Approximation for Compact Operators. {\it Math. Comput.} {\bf 29}(131), 712--725,
doi: \href{https://doi.org/10.2307/2005282}{10.2307/2005282}.

\bibitem{Pannekoucke-2007} 
Pannekoucke O., Berre, L., Desroziers, G. (2007): Filtering Properties of Wavelets for Local Background-Error Correlations. {\it Q. J. R. Meteorol. Soc.} {\bf 133}, 363--379, doi: \href{https://doi.org/10.1002/qj.33}{10.1002/qj.33}.

\bibitem{Pourahmadi} 
Pourahmadi, M. (2011): Covariance Estimation: The GLM and Regularization Perspectives. {\it Stat. Sci.} {\bf 26}(3), 369--387,\\
doi: \href{https://doi.org/10.1214/11-STS358}{10.1214/11-STS358}.

\bibitem{Rue_Held} 
Rue, H., Held, L. (2005): Gaussian Markov Random Fields: Theory And Applications, {\it Monographs on Statistics and Applied Probability}, Chapman \& Hall/CRC.

\bibitem{Simpson-2012} 
Simpson, D., Lindgren, F., Rue, H. (2012): Think Continuous: Markovian Gaussian Models in Spatial Statistics. {\it Spat. Stat.} {\bf 1}, 16--29,\\
doi: \href{https://doi.org/10.1016/j.spasta.2012.02.003}{10.1016/j.spasta.2012.02.003}.

\bibitem{Szatrowski-1980} 
Szatrowski, T.\,H. (1980): Necessary and Sufficient Conditions
for Explicit Solutions in the Multivariate Normal Estimation Problem for Patterned Means and Covariances. {\it Ann. Stat.} {\bf 8}(4), 802--810,\\
doi: \href{https://doi.org/10.1214/aos/1176345072}{10.1214/aos/1176345072}.

\bibitem{Turcicova-2019} 
Turčičová, M., Mandel J., Eben, K. (2019): Multilevel Maximum Likelihood Estimation with Application to Covariance Matrices. {\it Commun. Stat. -- Theory and Methods} {\bf 48}(4),  909--925,\\
doi: \href{https://doi.org/10.1080/03610926.2017.1422755}{10.1080/03610926.2017.1422755}.

\bibitem{smef7}
Turčičová, M., Mandel, J., Eben, K. (2020): Score Matching Filters for Gaussian Markov Fields with Linear Model of Precision Matrix. {\it Článek se připravuje.}

\bibitem{Ueno}
Ueno, G., Tsuchiya, T. (2009): Covariance Regularization in Inverse Space. {\it Q. J. R. Meteorol. Soc.} {\bf 135}(642), 1133--1156,\\
doi: \href{https://doi.org/10.1002/qj.445}{10.1002/qj.445}.

\end{thebibliography}



