Basic operators and adjoints

Causal and leaky integration

Causal integration is defined as:

$$y(t) \eq \int_{-\infty}^t x(\tau ) d\tau$$ (17)

Leaky integration is defined as:

$$y(t) \eq \int_0^\infty x(t-\tau) e^{-\alpha\tau} d\tau$$ (18)

As $\alpha \rightarrow 0$, leaky integration becomes causal integration. The word ``leaky'' comes from electrical circuit theory in which the voltage on a capacitor would be the integral of the current if the capacitor did not leak electrons.

Sampling the time axis gives a matrix equation that we should call ``causal summation,'' but we often call it ``causal integration.'' Equation (19) represents causal integration for $\rho=1$ and leaky integration for $0<\rho <1$.

$$\bold y \eq \left[ \begin{array}{c} y_0 \\ y_1 \\ y_2... ..._4 \\ x_5 \\ x_6 \end{array} \right] \eq \bold C \bold x$$ (19)

(The discrete world is related to the continuous by $\rho = e^{-\alpha\Delta\tau}$ and in some applications, the diagonal is 1/2 instead of 1.) Causal integration is the simplest prototype of a recursive operator. The coding is trickier than that for the operators we considered earlier. Notice when you compute $y_5$ that it is the sum of 6 terms, but that this sum is more quickly computed as $y_5 = \rho y_4 + x_5$. Thus, equation (19) is more efficiently thought of as the recursion

$$y_t \eq \rho\; y_{t-1} + x_t \quad \quad \quad t {\rm increasing}$$ (20)

(which may also be regarded as a numerical representation of the differential equation $dy/dt+y (1-\rho)/\Delta t=x(t)$.)

When it comes time to think about the adjoint, however, it is easier to think of equation (19) than of equation (20). Let the matrix of equation (19) be called $\bold C$. Transposing to get $\bold{C}\T$ and applying it to $\bold y$ gives us something back in the space of $\bold x$, namely $\tilde{\bold x} = \bold C\T \bold y$. From it we see that the adjoint calculation, if done recursively, needs to be done backward, as in:

$$\tilde x_{t-1} \eq \rho \tilde x_{t} + y_{t-1} \quad \quad \quad t {\rm decreasing}$$ (21)

Thus, the adjoint of causal integration is anticausal integration.

A module to do these jobs is leakint. The code for anticausal integration is not obvious from the code for integration and the adjoint coding tricks we learned earlier. To understand the adjoint, you need to inspect the detailed form of the expression $\tilde{\bold x} = \bold C\T \bold y$ and take care to get the ends correct. Figure 5 illustrates the program for $\rho=1$.

api/c/causint.c

    t = 0.;
    if (adj) {
	for (i=nx-1; i >= 0; i-) {
	    t += yy[i];
	    xx[i] += t;
	}
    } else {
	for (i=0; i <= nx-1; i++) {
	    t += xx[i];
	    yy[i] += t;
	}
    }

causint Figure 5. in1 is an input pulse. C in1 is its causal integral. C' in1 is the anticausal integral of the pulse. A separated doublet is in2. Its causal integration is a box, and its anticausal integration is a negative box. CC in2 is the double causal integral of in2. How can an equilateral triangle be built?

Later, we consider equations to march wavefields up toward the Earth surface, a layer at a time, an operator for each layer. Then, the adjoint starts from the Earth surface and marchs down, a layer at a time, into the Earth.

Basic operators and adjoints