Linear Algebra Implementations

DUBFI provides a linear algebra interface, see dubfi.linalg.

Note

If you intent to run an inversion with many data points, you will usually use the “mpipad” implementation. This implementation relies on some approximations that is only valid when configuring some localization and cutoff scales correctly.

Dense matrices (“dense”)

This implementation only uses numpy arrays. It is focussed on simplicity and can have acceptable performance for very small datasets. It quickly becomes inefficient and requires huge amounts of memory when working with larger datasets.

This implementation serves as a reference that can be compared to other implementations. It aims for a simple, readable and reliable implementation and not on performance.

CSC sparse matrices (“csc”)

This implementation uses the compressed-column representation of sparse matrices. This reduces memory requirements compared to dense matrices if the matrices are sparse. This implementation has limited performance but remains numerically exact.

MPI with dense matrices (“mpipad”)

This default implementation is optimized for solving the inversion problem using many parallel worker processes. It is optimized for good performance but requires a careful configuration (entry segment_buffer) to make sure that the results are correct.

The MPI parallelization is based on a central approximation: Matrix elements that are approximately zero in the localization matrix (used for constructing the R matrix) are negligible in all square matrices in observation space. This is physically motivated by the assumption that there are not long-range correlations in space or time. Based on this approximation, the observation space is split into chunks that are processed independent of each other. The linear algebra implementation using MPI parallelization is thus not numerically exact, but it is scalable to large observation datasets while preserving a good performance.

This implementation consists of one main process and a fix number of worker processes, each operating on one chunk of the observation space. The main process provides the functions required by the abstract linear interface. But the computation happens on the worker processes. These worker processes only follow instructions sent from the main process.

The main process first instructs the workers to create their own data providers of input data based on the configuration file. The workers then create the required matrices and vectors in observation space (e.g. the vector of model-observation mismatch). They can then do all required computations as requested by the main process.

Example

For example, the inversion may contain an operation like \(y^T R^{-1} y\) where \(y\) is a vector in observation space (observations minus model prediction) and \(R\) is a matrix. Each worker should already have its chunk of \(y\) and \(R\). The main process will send an instruction x = R.solve(y) to the workers, followed by y.apply(x). Then it waits for the results. The worker solves \(R x = y\) for \(x\) and assigns the result to a new vector x. Then it computes the inner product \(y^T x\) and sends the result back to the main process. The main process adds up the results to obtain the approximate solution of \(y^T R^{-1} y\) for approximately diagonal \(R\). (Note that the names will be more compicated in a real application.)

The linear algebra framework can usually be used agnostic of the implementation. However, care must be taken when using the MPI parallelization with separate worker processes. The user shall define how many worker processes shall be created (environment variable MPI_WORKERS) and ensure that the worker processes end when the main process ends (call dubfi.quit_clean() in the main process).