Click here to download the full example code
There is a large number of Matrix types that are all derived from the
MatrixBase. They do not have to store elements but can be
logical or wrappers. They just have to provide the functions
A.rows() (column and row numbers),
A.mult(x)$=Acdot x$ and
# We start off with the typical imports import numpy as np import pygimli as pg
All elements are stored column-wise, i.e. all rows
A[i] are of type
pg.Vector. This matrix is used for storing dense data (like ERT
Jacobians) and doing simple algebra.
RMatrix: 3 x 4 4 [1.0, 0.0, 0.0, 0.0] 4 [0.0, 1.0, 2.0, 3.0] 4 [0.0, 0.0, -1.0, 0.0] 3 [5.0, 44.0, -7.0]
Exists also as complex matrix under
Index-based sparse matrix
Sparse matrix (most elements are zero) with single access and. Typical for traveltime Jacobian (only certain cells covered by individual rays) or constraint matrices.
2 [2.0, 1.0]
Exists also as complex-valued variant
Used for numerical approximation of partial differential equations like
finite-element or finite volume. Not typically used unless efficiency is
of importance. It exists also complex-valued as
First, there is an identity matrix
IdentityMatrix. No elements
stored at all. Important for constraint matrices when combined into
BlockMatrix (see below). More generally, there is a diagonal matrix
DiagonalMatrix where only its diagonal is stored as vector.
3 [10.0, 11.0, 12.0]
Often, matrices are weighted from either side by a vector, e.g. data error weighting of the Jacobian matrix, data and model transformations, or weighting individual smoothness parts according to the roughness so that only the weighting is changed and not the matrix.
3 [6.0, 6.0, 6.0] 3 [12.0, 18.0, -6.0]
Combinations of matrices#
Logical matrices can combine different other matrices (of arbitrary
type) avoiding double memory storage by multiplication (
or addition (
2 [0.0, 0.0]
The most important type is the
BlockMatrix, where arbitrary matrices
are combined into a logical matrix. This is of importance for inversion
frameworks: * joint inversion: Jacobian matrices are concatenated *
combination of different constrains: combining different regularization
* laterally, spatially or temporally constrained inversion:
regularization between model cells of each frame but also between the
frames Note that the matrices only have to be defined once and can
A = pg.BlockMatrix() A1 = pg.Matrix(2, 2) A.addMatrix(A1, 0, 0) A.addMatrix(A1, 4, 0, scale=2.0) A2 = pg.matrix.IdentityMatrix(5) A.addMatrix(A2, 1, 2) A3 = pg.SparseMapMatrix(2, 3) A.addMatrix(A3, 2, 3) print(A) ax, _ = pg.show(A)
pg.matrix.BlockMatrix of size 6 x 7 consisting of 4 submatrices.
There are also simpler types of matrix combinations: *
V2Matrix: two matrices below/next to each other *
VNMatrix: one matrix repeated N times
NDMatrix: block diagonal matrix
TransposedMatrix avoids transposing any matrix by exchanging left/right mult. SquaredMatrix keeps only the matrix A but works as A^T @ A SquaredTransposeMatrix keeps only the matrix A but works as A @ A.T RealNumpyMatrix holds a real-valued numpy array ComplexNumpyMatrix holds a complex-valued numpy array in a real pg matrix NumpyMatrix
Often, several matrices or even the same one have to be combined. RepeatHMatrix, RepeatVMatrix, RepeatDMatrix hold a single matrix that is repeated horizontally, vertically or diagonally. NDMatrix, FrameConstraintMatrix is a special generator for constraining cells of every (e.g. timelapse) frame and moreover the frames with each other.
F = pg.matrix.FrameConstraintMatrix(A3, 3) ax, _ = pg.show(F) print(F)
pg.matrix.BlockMatrix of size 12 x 9 consisting of 7 submatrices.
Geostatistical constraint matrix#
For geostatistical constraints, a correlation matrix is computed using correlation lengths and angles to define their directions. To access its inverse root in a way that avoids matrix inversion, an eigenvalue decomposition is done and the eigenvalues $D$ and -vectors $Q$ are stored so that the operator
Jordi, C., Doetsch, J., Günther, T., Schmelzbach, C. & Robertsson, J.O.A. (2018): Geostatistical regularisation operators for geophysical inverse problems on irregular meshes. Geophysical Journal International 213, 1374- 1386, doi:10.1093/gji/ggy055.