Stochastic Properties of Random Matrix Factors - 1

In this series of posts, we will be considering the stochastic properties of factors of random matrices. At first, we look at the simple case of Cholesky decomposition of Wishart matrices. More below.

Let $H$ be a matrix whose elements are random variables, each with a given distribution. Now assume that we factorize $H$ using one of the matrix factorization techniques such as QR, Eigenvalue, or Cholesky decomposition. The elements of the matrix factor are also random variables... but what are their distributions?

This is exactly what we hope to answer.

First and the obvious reason to find the distributions is to understand the characteristics of the matrix factor. For example, in a wireless communication MIMO channel, which is a zero-mean Gaussian random matrix, it helps us understand the stochastic properties of a "partial" or "inverse" channel.

The second reason is to speed up simulations. If we have a simulation that requires a random matrix factor which is computed numerically, the simulation may be sped up by generating the elements of the matrix factor directly with the given distributions.

Now we look at the stochastic properties of Cholesky decomposition of Wishart matrices.

Let $G \in \mathbb{C}^{n \times n}$ be a complex-valued matrix whose elements are Gaussian random variables i.e. $\mathcal{N}(0,1)$. Let $H = G^H G$. $H$ (Wishart matrix) is now a Hermitian symmetric complex-valued matrix.

It is well known that the Cholesky decomposition of $H$ is the upper triangular matrix from the QR decomposition of $G$, i.e. $H = R^H R$ where $G = QR$ for a Unitary matrix $Q$ and an upper triangular matrix with positive diagonal elements $R$.

Therefore, the distribution of $R$ is straightforward from [1] and is given as follows:
$$
\begin{bmatrix}
\chi_n & \mathcal{N} & \cdots & \mathcal{N} \\
           & \chi_{n-1} & \cdots & \mathcal{N} \\
           &                   & \ddots & \vdots \\
           &                   &            & \chi_1
\end{bmatrix},
$$

where $\mathcal{N}$ denotes Gaussian distribution and $\chi_k$ the Chi distribution with $k$ degrees of freedom.

Stochastic properties of the matrix factors may also be verified numerically. Let $G \in \mathbb{R}^{2 \times 2}$ be a real-valued matrix with elements distributed according to standard Gaussian distribution, i.e. $g_{ij}$'s are $\mathcal{N}(0,1)$. $G \leadsto R$ given by:
$$
\begin{bmatrix}
\mathcal{N}(0,1) & \mathcal{N}(0,1) \\
\mathcal{N}(0,1) & \mathcal{N}(0,1)
\end{bmatrix} \leadsto
\begin{bmatrix}
\chi_2 & \mathcal{N}(0,1) \\
 & \chi_1
\end{bmatrix}
$$

Figure below illustrates the distributions of elements of $R$ computed numerically using MATLAB. Normalized numerical quantities obtained over 100000 trials are plotted in the blue bar graph and the theoretical expectation is plotted in red. Numerical results closely match the theoretical expectation.

Numerical vs theoretical distribution of elements of R

In the next blog posts on this topic, I hope to be able to provide more examples for stochastic properties of matrix factors especially the ones we come across in wireless communication.

References:
[1] Edelman, Alan, and N. Raj Rao. "Random matrix theory." Acta Numerica 14 (2005): 233-297.

Inverse Cholesky Factor of a positive definite Hermitian symmetric Toeplitz matrix using Durbin recursions

Let $T \in \mathbb{C}^{n \times n}$ be a positive-definite Hermitian symmetric Toeplitz matrix defined with elements $t_{ij} = r_{i-j}$ for a set of scalars $r_{-n+1}, \cdots, r_0, \cdots, r_{n-1}$, where $r_{-t} = r_t^*$. Without loss of generality, assume $r_0 = 1$.

We have the Cholesky factor $R$ of the matrix $T$ such that $R^H R = T$. The inverse Cholesky factor is then given by $R^{-1}$. A closely related decomposition known as the LDL decomposition is given as follows: $\Lambda = \text{diag}(R)$, $L = R^H \Lambda^{-1}$, and therefore, $T = L \Lambda^2 L^H = L D L^H$. Matrix $L$ has ones along the principal diagonal. Let $\Sigma = LD$.

We have $T L^{-H} = \Sigma$, where $L^{-H}$ is the Hermitian transpose of the inverse with ones along the diagonal. As $\Sigma$ is lower triangular, we can solve for the strictly upper triangular parts of $L^{-H}$ as the R.H.S. is known to be zero.

For the $k$-th column of $L^{-H}$ we have the first $k-1$ elements $z_{k-1}$, upon simplification, as:

$$T_{k-1} z_{k-1} = - \begin{bmatrix} r_{k+1} \\ \vdots \\ r_1 \end{bmatrix}$$

For Hermitian symmetric Toeplitz matrix $T$ we have $T = E_n T^* E_n$, where $T^*$ is complex conjugate of the matrix and $E_n$ is the $n \times n$ exchange matrix. Using this property, and setting $y_{k-1} = E_k z^*_{k-1}$ we have:

$$T_{k-1} y_{k-1} = - \begin{bmatrix} r^*_1 \\ \vdots \\ r^*_{k+1} \end{bmatrix}$$,

which may be solved efficiently using the Durbin iterations [1]. The algorithm for finding $W = R^{-1} = [w_{i,j}]$ is given as follows:

$y_1 = -r^*_1; \enspace \alpha = -r^*_1; \enspace W = I_n$

$\beta = 1 - |\alpha|^2$

$w_{1,2} = y^*_1$

$w_{:,2} = w_{:,2}/\sqrt\beta$

$\text{for}~ k = 1:n-2$

$\quad \alpha = -(r^*_{k+1} + r^H_{k:-1:1} y_{1:k}) /\beta$

$\quad y_{1:k} = y_{1:k} + \alpha y^*_{k:-1:1}; \enspace y_{k+1} = \alpha$

$\quad \beta = (1 - |\alpha|^2) \beta$

$\quad w_{1:k+1,k+2} = y^*_{k+1:-1:1}$

$\quad w_{:,k+2} = w_{:,k+2}/\sqrt\beta$

$\text{end}$

The algorithm requires $2n^2$ flops [1] for the Durbin iterations, $n^2/2$ flops for scaling with $1/\sqrt\beta$ and $n$ inverse square root computations; and may be easily adapted to find $L^{-1}$ and $D$ of the $LDL$ decomposition.

A MATLAB reference code is available in [2].

References:

[1] Gene H. Golub, Charles F. Van Loan, Matrix Computations, Third Edition, Algorithm 4.7.1 (Durbin).

[2] http://de.mathworks.com/matlabcentral/fileexchange/52847-invchol-durbin-t-