the Multivariate linear regression model is given by

$\underset{n \times d}{\mathbf{Y}} = \underset{n \times k}{\mathbf{X}} \hspace{2mm}\underset{k \times d}{\mathbf{B}} + \underset{n \times d}{\mathbf{E}}$

With $\mathbf{Y}$ is data, $\mathbf{X}$ is a non stochastic design matrix, $\mathbf{E}$ are normally distributed model errors, and $\mathbf{B}$ are the regression coefficients to be estimated. Also $n$ is the number of observation (e.g.,time samples in my specific context), $d$ the number of variables (e.g., time series at different pixels of an image) and $k$ the number of explanatory variables.

All the books i read until now treat this problems using the assumption that the obervations (the different rows of $\mathbf{E}$) are independants whereas the different variables (the different columns of $\mathbf{E}$) are dependant with the $d \times d$ covariance structure $\boldsymbol\Sigma_d$.

In this context, the model can be reparametrized as a multiple linear regression model: $\underset{nd \times 1}{\mathbf{y}} = \underset{nd \times nk}{\mathbf{D}} \hspace{2mm} \underset{nk \times 1}{\boldsymbol{\beta}} + \underset{nd \times 1}{\mathbf{e}}$ with $\mathbf{e}\sim \mathcal{N}(\mathbf{0},\mathbf{I}_n \otimes\boldsymbol\Sigma_d$) where $\otimes$ is the kronecker product and solved using Generalized Least Square or Feasible GLS if $\boldsymbol\Sigma_d$ is unknown.


1-What happens when the observation are not independants ?.

2-Can we reparametrize the model using with $\mathbf{e}\sim \mathcal{N}(\mathbf{0},\mathbf{\Omega}_n \otimes\boldsymbol\Sigma_d$) and solve it using the GLS when $\mathbf{\Omega}_n$ and $\boldsymbol\Sigma_d$ are known .

3- Why this general problem (non independance of observations) is not treated in textbooks . Is it because the FGLS is too complex (i.e., find an estimate of both $\mathbf{\Omega}_n$ and $\boldsymbol\Sigma_d$)


  • $\begingroup$ Why do you think that GLS would not work? $\endgroup$
    – user83346
    Commented Aug 31, 2016 at 14:42
  • $\begingroup$ I think it should work but because I am new to the domain i am asking the community its opinion. I was very surprised that textbooks just describes the multivariate problem using the assumption of independant observations but never dare to explain the problem with both dependant variables and observations . $\endgroup$
    – Alexis
    Commented Aug 31, 2016 at 14:56
  • $\begingroup$ It's probably not a question of daring. The model you suggest in point 2. is common and is a version of the matrix normal distribution. Most work I've seen on that model is concerned with maximum likelihood. Note that it is just a special case (Kronecker-decomposable cov. matrix) of a GLS model where $\mathrm{vec}Y$ is normal with an unrestricted covariance matrix. $\endgroup$
    – KOE
    Commented Aug 31, 2016 at 14:59
  • 1
    $\begingroup$ Sure, here is a pre-print of a paper that deals with the existence of MLEs of the covariance matrix in the matrix normal. Their references includes some applied work. Note that most of these are focused on estimation of the covariance matrix in that model though, and what they call space and time indices in your case is indices for observations and pixels, if I understand your question correctly. arxiv.org/pdf/1410.2118.pdf $\endgroup$
    – KOE
    Commented Aug 31, 2016 at 15:11
  • 2
    $\begingroup$ Thx you very much, i will read it carefully. PS: By the way, this paper treats of time/space correlation structure of EEG data and my research interest is close to it ( functional near infrared spectroscopy), so I deal with time series datasets among different channels too. $\endgroup$
    – Alexis
    Commented Aug 31, 2016 at 15:19

1 Answer 1


1. If the rows of $\mathbf{Y}$ are not independent, there is in general nothing one can say about the structure of $\mathbf{y}:=\mathrm{vec}(\mathbf{Y})$ other than $\mathbf{y}\sim N\left((\mathbf{I} \otimes \mathbf{X})\mathrm{vec}(\mathbf{B}), \mathbf{\Psi}\right)$ for some covariance matrix $\mathbf{\Psi}$, which is the classical GLS model.

2. Yes you can, given that you are willing to assume this more restrictive structure and that you also know the matrices. Regarding restrictions, you may notice that this decomposition implies that the autocovariance function for any one pixel is proportional to that of any other.

If you are indeed willing to make this assumption and know the matrices, then just set $\mathbf{\Psi} = \mathbf{\Omega}_n \otimes \mathbf{\Sigma}_d$ and your GLS estimator of $\mathbf{b} = \mathrm{vec}(\mathbf{B})$ is given by $\mathbf{Z}(\mathbf{Z}'\mathbf{\Psi}^{-1}\mathbf{Z})^{-1}\mathbf{Z}'\mathbf{\Psi}^{-1}\mathbf{y}$, where $\mathbf{Z} = \mathbf{I} \otimes \mathbf{X}$.

3. I believe your guess at this point is pretty much correct. The model you suggest for $\mathbf{e}$ is equivalent to assuming a matrix normal distribution for $\mathbf{E}$. This paper considers estimation of the component covariance matrices $\mathbf{\Omega}_n$ and $\mathbf{\Sigma}_d$ using maximum likelihood. Paywall link to published version in Journal of Multivariate Analysis, 2016.


There is in fact a newer paper (here is an arxiv version) which settles the uniqueness and existence question for maximum likelihood in this model.

Their Theorem 1 in section 4 gives the following result. Suppose we observe $N$ independent copies of the matrix normal variable $\mathbf{Y}$ and $\mathbf{X} = \mathbf{1}_n$. Then:

  1. If $N < \max(n/d, d/n) + 1$ the MLE of $\mathbf{\Psi}$ does not exist.
  2. If $N > d/n + n/d + 1$ the MLE exists and is unique almost surely.
  3. For $N \in [\max(n/d, d/n) + 1, d/n + n/d + 1]$ there are examples of when the MLE exists and when it does not, so no general statement can be made.

In particular, note that your case corresponds to having just $N=1$ observation.

  • $\begingroup$ Thx for the update, very usefull $\endgroup$
    – Alexis
    Commented Oct 24, 2016 at 15:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.