17. Real/Complex Inversion

In this tutorial we will discuss two equivalent approaches to the solution of inverse problems with real-valued model vector and complex-valued data vector. In other words, we consider a modelling operator \(\mathbf{A}:\mathbb{F}^m \to \mathbb{C}^n\) (which could be the case for example for the real FFT).

Mathematically speaking, this problem can be solved equivalently by inverting the complex-valued problem:

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

or the real-valued augmented system

\[\begin{split}\DeclareMathOperator{\Real}{Re} \DeclareMathOperator{\Imag}{Im} \begin{bmatrix} \Real(\mathbf{y}) \\ \Imag(\mathbf{y}) \end{bmatrix} = \begin{bmatrix} \Real(\mathbf{A}) \\ \Imag(\mathbf{A}) \end{bmatrix} \mathbf{x}\end{split}\]

Whilst we already know how to solve the first problem, let’s see how we can solve the second one by taking advantage of the real method of the pylops.LinearOperator object. We will also wrap our linear operator into a pylops.MemoizeOperator which remembers the last N model and data vectors and by-passes the computation of the forward and/or adjoint pass whenever the same pair reappears. This is very useful in our case when we want to compute the real and the imag components of

import matplotlib.pyplot as plt
import numpy as np

import pylops

plt.close("all")
np.random.seed(0)

To start we create the forward problem

n = 5
x = np.arange(n) + 1.0

# make A
Ar = np.random.normal(0, 1, (n, n))
Ai = np.random.normal(0, 1, (n, n))
A = Ar + 1j * Ai
Aop = pylops.MatrixMult(A, dtype=np.complex128)
y = Aop @ x

Let’s check we can solve this problem using the first formulation

A1op = Aop.toreal(forw=False, adj=True)
xinv = A1op.div(y)

print(f"xinv={xinv}\n")

xinv=[1.+9.99200722e-16j 2.-1.55431223e-15j 3.+1.33226763e-15j
 4.-8.88178420e-16j 5.+6.66133815e-16j]

Let’s now see how we formulate the second problem

Amop = pylops.MemoizeOperator(Aop, max_neval=10)
Arop = Amop.toreal()
Aiop = Amop.toimag()

A1op = pylops.VStack([Arop, Aiop])
y1 = np.concatenate([np.real(y), np.imag(y)])
xinv1 = np.real(A1op.div(y1))

print(f"xinv1={xinv1}\n")

xinv1=[1. 2. 3. 4. 5.]