    Seismic wave extrapolation using lowrank symbol approximation  Next: EXAMPLES Up: Fomel, Ying, & Song: Previous: WAVE EXTRAPOLATION

# LOWRANK APPROXIMATION

The key idea of the lowrank decomposition is decomposing the wave extrapolation matrix (11)

for a fixed into a separated representation (12)

Representation (12) speeds up the computation of since     (13)

The evaluation of the last formula is effectively equivalent to applying inverse Fast Fourier Transforms. Physically, a separable lowrank approximation amounts to selecting a set of representative spatial locations and representative wavenumbers.

In order to discuss the construction of approximation (12), let us view it as a matrix decomposition problem (14)

where is the matrix with entries , is the submatrix of that consists of the columns associated with , is the submatrix that consists of the rows associated with , and . In practice, we find that the matrix has a low rank separated representation provided that is sufficiently small, which, in the case of smooth models, can be partially explained by the separation of terms in the Taylor series 5. Let be a prescribed accuracy of this separated representation, and be the numerical rank of . The construction of the separated representation in equation (14) follows the method of Engquist and Ying (2007, 2009) and is detailed in the appendix. The main observation is that the columns of and the rows of should span the column space and row space of , respectively, as well as possible. The algorithm for computing (14) takes the following steps:
1. Pick a uniformly random set of columns of where is chosen to be 3 or 4 in practice. Perform the pivoted QR factorization to (Golub and Van Loan, 1996). The first pivoted columns correspond to rows of the matrix . Define to be the submatrix of that consists of these rows and set with to be the corresponding values of these rows.
2. Pick a uniformly random set of rows of and perform the pivoted QR factorization to . Define to be the submatrix of that consists of these columns and set with to be the corresponding values of these columns.
3. Set the middle matrix where stands for the pseudoinverse.
4. Combining the result of the previous three steps gives the required separated representation .
The algorithm does not require, at any step, access to the full matrix , only to its selected rows and columns. Once the decomposition is complete, it can be used at every time step during the wave extrapolation process. In multiple-core implementations, the matrix operations in equation (12) are easy to parallelize. The algorithm details are outlined in the appendix.

The cost of the algorithm is operations per time step, where refers to the cost of the Fourier transform. In comparison, the cost of finite-difference wave extrapolation is , where is the size of the finite-difference stencil. Song et al. (2011) present an application of the proposed lowrank approximation algorithm for devising accurate finite-different schemes. There is a natural trade-off in the selection of : larger values lead to a more accurate wave representation but require a longer computational time. In the examples of the next section, we select these parameters based on an estimate of the approximation accuracy and generally aiming for the relative accuracy of . The resulting is typically smaller than the number of Fourier transforms required for pseudo-spectral algorithms such as pseudo-spectral implementations of the rapid expansion method (Pestana and Stoffa, 2011).    Seismic wave extrapolation using lowrank symbol approximation  Next: EXAMPLES Up: Fomel, Ying, & Song: Previous: WAVE EXTRAPOLATION

2013-08-31