Matrix inversion

Many signal processing algorithms are described using inverses and square rootsof matrices. However, their actual computation is very rarely required and oneshould resort to alternatives in practical implementations instead. Below, we willdescribe 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 examplex_=tf.random.normal([10,1])h=tf.random.normal([10,10])y=tf.linalg.matmul(h,x_)# Solve via matrix inversionh_inv=tf.linalg.inv(h)x=tf.linalg.matmul(h_inv,y)

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

\[\mathbf{H}\mathbf{x} = \mathbf{y}\]

which looks in code like this:

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

When\(\mathbf{H}\) is a Hermitian positive-definite matrix, we can leveragetheCholeskydecomposition\(\mathbf{H}=\mathbf{L}\mathbf{L}^{\mathsf{H}}\), where\(\mathbf{L}\) isa lower-triangular matrix, for aneven faster implementation:

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

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

Correlated random vectors

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

\[\mathbf{x} = \mathbf{R}^{\frac12} \mathbf{w}\]

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

# Create covariance matrixr=tf.constant([[1.0,0.5,0.25],[0.5,1.0,0.5],[0.25,0.5,1.0]])# Cholesky decompositionl=tf.linalg.cholesky(r)# Create batch of correlated random vectorsw=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}^{-\frac12}\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}^{-\frac12}\), it is sufficient to compute\(\mathbf{L}^{-1}\), which can beachieved by solving the linear system\(\mathbf{L} \mathbf{X} = \mathbf{I}\),exploiting the diagonal structure of the Cholesky decomposition\(\mathbf{L}\):

# Create covariance matrixr=tf.constant([[1.0,0.5,0.25],[0.5,1.0,0.5],[0.25,0.5,1.0]])# Cholesky decompositionl=tf.linalg.cholesky(r)# Inverse of Ll_inv=tf.linalg.triangular_solve(l,tf.eye(3),lower=True)