Vectorizing Regression With Kronecker Products

Why Vectorize a Multivariate Regression?

A multivariate regression writes several responses for each case in one matrix equation. That compact form makes the shared predictor structure visible, but it can hide the ordinary linear-system structure underneath. Vectorization stacks a matrix into a single column; Kronecker products then express the regression design and error covariance as matrices acting on that column.

The choice of stacking order matters. We use column-wise vectorization throughout, so the entries of the first response column come first, followed by those of the second response column, and so on. Keeping that order explicit prevents many otherwise plausible-looking but incorrect formulas.

Definition: For a matrix $A$ with columns $a_1,\ldots,a_k$, $\operatorname{vec}(A)$ is the column vector formed by stacking those columns in order: $$ \operatorname{vec}(A)= \begin{bmatrix} a_1\\ a_2\\ \vdots\\ a_k \end{bmatrix}. $$ For the multivariate regression model, let $Y$ be an $n\times q$ response matrix, $X$ an $n\times p$ design matrix, $B$ a $p\times q$ coefficient matrix, and $U$ an $n\times q$ error matrix. The model is $Y=XB+U$.

The dimensions determine the vectorized model: $\operatorname{vec}(Y)$ and $\operatorname{vec}(U)$ each have length $nq$, while $\operatorname{vec}(B)$ has length $pq$. For matrices $A$ and $C$, their Kronecker product $A\otimes C$ is the block matrix whose $(i,j)$ block is $A_{ij}C$. In particular, $I_q\otimes X$ is a block diagonal matrix with $q$ copies of $X$, each of size $n\times p$.

The Vectorized Regression Identity

For compatible matrices $A$, $M$, and $C$, column-wise vectorization satisfies $\operatorname{vec}(AMC)=(C^T\otimes A)\operatorname{vec}(M)$. Taking $A=X$, $M=B$, and $C=I_q$ yields the central identity for this model.

Theorem: If $Y=XB+U$, with dimensions as defined above, then $$ \operatorname{vec}(Y)=(I_q\otimes X)\operatorname{vec}(B)+\operatorname{vec}(U). $$ The vectorized design matrix is $I_q\otimes X$, with dimensions $nq\times pq$. Its action on $\operatorname{vec}(B)$ applies the same predictor matrix $X$ to the coefficient column for each response.

To see why the identity holds, write $B=[b_1,\ldots,b_q]$, $Y=[y_1,\ldots,y_q]$, and $U=[u_1,\ldots,u_q]$. The matrix equation gives $y_j=Xb_j+u_j$ for every $j=1,\ldots,q$. Stacking those $q$ equations in order gives the vectorized result: the first block of $I_q\otimes X$ multiplies $b_1$, the second multiplies $b_2$, and so forth.

Worked Example: Writing a Two-Response Model as One Linear System

Let $$ X=\begin{bmatrix}1&2\\0&1\end{bmatrix},\quad B=\begin{bmatrix}2&-1\\3&4\end{bmatrix},\quad U=\begin{bmatrix}0.5&-1\\1&2\end{bmatrix}. $$ Compute $Y=XB+U$, then verify the vectorized equation.

First calculate the fitted response matrix: $$ XB= \begin{bmatrix}1&2\\0&1\end{bmatrix} \begin{bmatrix}2&-1\\3&4\end{bmatrix} = \begin{bmatrix}8&7\\3&4\end{bmatrix}. $$ Therefore $$ Y=XB+U= \begin{bmatrix}8.5&6\\4&6\end{bmatrix}. $$

Column stacking gives $\operatorname{vec}(B)=[2,3,-1,4]^T$, $\operatorname{vec}(U)=[0.5,1,-1,2]^T$, and $\operatorname{vec}(Y)=[8.5,4,6,6]^T$. The vectorized design is $$ I_2\otimes X= \begin{bmatrix} 1&2&0&0\\ 0&1&0&0\\ 0&0&1&2\\ 0&0&0&1 \end{bmatrix}. $$ Multiplication gives $$ (I_2\otimes X)\operatorname{vec}(B) = \begin{bmatrix}8\\3\\7\\4\end{bmatrix}, \qquad (I_2\otimes X)\operatorname{vec}(B)+\operatorname{vec}(U) = \begin{bmatrix}8.5\\4\\6\\6\end{bmatrix} =\operatorname{vec}(Y). $$ The first two entries correspond to response 1, and the last two to response 2.

Error Covariance in Vectorized Form

Suppose the cases are independent, each error row has mean zero, and each row has the same response covariance matrix $\Sigma$, a $q\times q$ matrix. Thus, for rows $u_i^T$ and $u_j^T$ of $U$, the response errors within a case may be correlated, while errors from different cases are uncorrelated. These assumptions imply a particularly simple covariance for the column-stacked error vector.

Result: Under independent cases with $\operatorname{Cov}(u_i)=\Sigma$ for each case $i$, and zero cross-case covariances, $$ \operatorname{Cov}(\operatorname{vec}(U))=\Sigma\otimes I_n. $$ With this covariance, the vectorized model has design matrix $I_q\otimes X$ and error covariance matrix $\Sigma\otimes I_n$.

For a derivation, label an entry of the stacked vector by its response and case: $U_{ir}$ occurs at position $i+(r-1)n$. The entry of its covariance matrix at positions $i+(r-1)n$ and $j+(s-1)n$ is $\operatorname{Cov}(U_{ir},U_{js})$. Independence across cases makes this zero when $i\ne j$. When $i=j$, the within-case covariance assumption gives $\operatorname{Cov}(U_{ir},U_{is})=\Sigma_{rs}$. Thus, in all cases, $$ \operatorname{Cov}(U_{ir},U_{js})=\Sigma_{rs}\,\mathbf{1}\{i=j\}, $$ which is exactly the corresponding entry of $\Sigma\otimes I_n$. The order of the factors is tied to column stacking: response indices index the blocks, and case indices run within each block.

Worked Example: Building the Error Covariance Matrix

Suppose there are $n=2$ cases and $q=2$ responses, with $$ \Sigma= \begin{bmatrix}4&1\\1&9\end{bmatrix}. $$ Write the covariance matrix of $\operatorname{vec}(U)$ in the stated column-wise order.

The vector order is $[U_{11},U_{21},U_{12},U_{22}]^T$: the two observations of response 1 precede the two observations of response 2. By the result above, $$ \operatorname{Cov}(\operatorname{vec}(U)) =\Sigma\otimes I_2 = \begin{bmatrix} 4I_2&1I_2\\ 1I_2&9I_2 \end{bmatrix} = \begin{bmatrix} 4&0&1&0\\ 0&4&0&1\\ 1&0&9&0\\ 0&1&0&9 \end{bmatrix}. $$ For example, the covariance between $U_{11}$ and $U_{12}$, the two responses on case 1, is $1=\Sigma_{12}$. The covariance between $U_{11}$ and $U_{22}$, which are different responses on different cases, is zero. The diagonal entries $4,4,9,9$ are the variances of the corresponding response errors.

Translating Coefficient Operations

Many questions about a regression concern selected combinations of coefficients rather than the entire matrix $B$. The Kronecker identity translates these operations without changing their meaning. Let $L$ be an $a\times p$ matrix acting on predictor coefficients and $R$ a $q\times b$ matrix acting on responses. The transformed coefficient matrix $LBR$ has size $a\times b$.

Coefficient transformation: For compatible $L$, $B$, and $R$, $$ \operatorname{vec}(LBR)=(R^T\otimes L)\operatorname{vec}(B). $$ In particular, a scalar contrast $a^T B c$, for $a\in\mathbb{R}^p$ and $c\in\mathbb{R}^q$, is $$ a^T Bc=(c^T\otimes a^T)\operatorname{vec}(B). $$

The general formula follows by applying the vectorization identity with $A=L$, $M=B$, and $C=R$. For the scalar formula, set $L=a^T$ and $R=c$; then $R^T\otimes L=c^T\otimes a^T$, and vectorizing a $1\times1$ matrix leaves its single entry unchanged. This also shows the order of the contrast factors: the response contrast $c$ appears first in the Kronecker product, while the predictor contrast $a$ appears second.

Worked Example: Expressing a Coefficient Contrast as a Vector Product

Use $$ B=\begin{bmatrix}2&-1\\3&4\end{bmatrix},\quad a=\begin{bmatrix}1\\-1\end{bmatrix},\quad c=\begin{bmatrix}1\\2\end{bmatrix}. $$ Find the scalar $a^T Bc$ both directly and using $\operatorname{vec}(B)$.

First, $Bc=[0,11]^T$, so $$ a^TBc= \begin{bmatrix}1&-1\end{bmatrix} \begin{bmatrix}0\\11\end{bmatrix} =-11. $$

Column-wise vectorization gives $\operatorname{vec}(B)=[2,3,-1,4]^T$. The contrast row is $$ c^T\otimes a^T= \begin{bmatrix}1&2\end{bmatrix}\otimes\begin{bmatrix}1&-1\end{bmatrix} = \begin{bmatrix}1&-1&2&-2\end{bmatrix}. $$ Thus $$ (c^T\otimes a^T)\operatorname{vec}(B) =2-3-2-8=-11, $$ in agreement with the direct calculation.

This translation is useful whenever a coefficient hypothesis or summary is stated as a linear transformation of $B$. The matrix expression identifies which predictor and response combinations are involved; the vector expression gives the corresponding row or matrix acting on $\operatorname{vec}(B)$. The dimensions provide a quick check: $R^T\otimes L$ has size $ab\times pq$, matching $\operatorname{vec}(LBR)\in\mathbb{R}^{ab}$ and $\operatorname{vec}(B)\in\mathbb{R}^{pq}$.

Common Mistakes and Checks

The most frequent error is silently changing from column stacking to row stacking. With column stacking, $\operatorname{vec}(XB)=(I_q\otimes X)\operatorname{vec}(B)$, and the error covariance under independent cases is $\Sigma\otimes I_n$. Row stacking uses a different ordering and therefore requires different Kronecker factor orders. A formula from one convention cannot be used unchanged with the other.

Wrong reasoning: “There are $q$ responses, so the design must be $X\otimes I_q$.” This overlooks the specified stacking order. Correction: Write the response columns in order and check that the first block must map the first $p$ entries of $\operatorname{vec}(B)$ to the first $n$ entries of $\operatorname{vec}(Y)$. That gives $I_q\otimes X$.
Wrong reasoning: “Errors are correlated, so every entry in $\operatorname{Cov}(\operatorname{vec}(U))$ may be nonzero.” This confuses within-case response correlation with across-case dependence. Correction: Under the stated assumptions, entries for distinct cases have covariance zero; only same-case response pairs inherit the entries of $\Sigma$.
Wrong reasoning: “The contrast is $a^T\otimes c^T$, since it combines $a$ and $c$.” Correction: For column-wise stacking, $a^TBc=(c^T\otimes a^T)\operatorname{vec}(B)$. Verify the order by checking the scalar dimensions and, when in doubt, testing a small numerical matrix.

Key takeaway: Under column-wise stacking, the multivariate regression design is $I_q\otimes X$, the covariance of independent-case errors with response covariance $\Sigma$ is $\Sigma\otimes I_n$, and a coefficient transformation $B\mapsto LBR$ becomes $R^T\otimes L$ acting on $\operatorname{vec}(B)$.

Check Your Understanding

Use column-wise vectorization in each question.

If $n=5$, $p=3$, and $q=4$, what are the dimensions of $\operatorname{vec}(Y)$, $\operatorname{vec}(B)$, and $I_q\otimes X$?
For independent cases with response covariance $\Sigma$, explain why errors from two different cases have zero covariance even when the responses within a case are correlated.
Write the vectorized coefficient transformation corresponding to $LBR$, and state its dimension if $L$ is $2\times p$ and $R$ is $q\times3$.
For predictor contrast $a$ and response contrast $c$, which Kronecker-product row multiplies $\operatorname{vec}(B)$ to produce $a^TBc$?
What ordering does $\operatorname{vec}(U)$ use when $U$ has $n$ rows and $q$ columns?