Matrix inversion

Many signal processing algorithms are described using inverses and square roots of matrices. However, their actual computation is very rarely required and one should resort to alternatives in practical implementations instead. Below, we will describe two representative examples.

Solving linear systems

One frequently needs to compute equations of the form

$$\mathbf{x}={\mathbf{H}}^{-1}\mathbf{y}$$

and would be tempted to implement this equation in the following way:

# Create random example
x_ = tf.random.normal([10, 1])
h = tf.random.normal([10, 10])
y = tf.linalg.matmul(h, x_)

# Solve via matrix inversion
h_inv = tf.linalg.inv(h)
x = tf.linalg.matmul(h_inv, y)

A much more stable and efficient implementation avoids the inverse computation and solves the following linear system instead

$$\mathbf{H}\mathbf{x}=\mathbf{y}$$

which looks in code like this:

# Solve as linar system
x = tf.linalg.solve(h, y)

When $\mathbf{H}$ is a Hermitian positive-definite matrix, we can leverage the Cholesky decomposition $\mathbf{H}=\mathbf{L}{\mathbf{L}}^{\mathsf{H}}$, where $\mathbf{L}$ is a lower-triangular matrix, for an even faster implementation:

# Solve via Cholesky decomposition
l = tf.linalg.cholesky(h)
x = tf.linalg.cholesky_solve(l, y)

This is the recommended approach for solving linear systems that we use throughout Sionna.

Correlated random vectors

Assume that we need to generate a correlated random vector with a given covariance matrix, i.e.,

$$\mathbf{x}={\mathbf{R}}^{\frac{1}{2}}\mathbf{w}$$

where $\mathbf{w}\sim \mathcal{CN}\mathcal{(}\mathbf{0}\mathcal{,}\mathbf{I}\mathcal{)}$ and $\mathbf{R}$ is known. One should avoid the explicit computation of the matrix square root here and rather leverage the Cholesky decomposition for a numerically stable and efficient implementation. We can compute $\mathbf{R}=\mathbf{L}{\mathbf{L}}^{\mathsf{H}}$ and then generate $\mathbf{x}=\mathbf{L}\mathbf{w}$, which can be implemented as follows:

# Create covariance matrix
r = tf.constant([[1.0, 0.5, 0.25],
                 [0.5, 1.0, 0.5],
                 [0.25, 0.5, 1.0]])

# Cholesky decomposition
l = tf.linalg.cholesky(r)

# Create batch of correlated random vectors
w = tf.random.normal([100, 3])
x = tf.linalg.matvec(l, w)

It also happens, that one needs to whiten a correlated noise vector

$$\mathbf{w}={\mathbf{R}}^{-\frac{1}{2}}\mathbf{x}$$

where $\mathbf{x}$ is random whith covariance matrix $\mathbf{R}=\mathbb{E}[\mathbf{x}{\mathbf{x}}^{\mathsf{H}}]$. Rather than computing ${\mathbf{R}}^{-\frac{1}{2}}$, it is sufficient to compute ${\mathbf{L}}^{-1}$, which can be achieved by solving the linear system $\mathbf{L}\mathbf{X}=\mathbf{I}$, exploiting the diagonal structure of the Cholesky decomposition $\mathbf{L}$:

# Create covariance matrix
r = tf.constant([[1.0, 0.5, 0.25],
                 [0.5, 1.0, 0.5],
                 [0.25, 0.5, 1.0]])

# Cholesky decomposition
l = tf.linalg.cholesky(r)

# Inverse of L
l_inv = tf.linalg.triangular_solve(l, tf.eye(3), lower=True)