There will be a webinar on pyGIMLi hosted by SEG on March 19, 2024 at 4 pm CET. Register for free here.

# Matrices#

There is a large number of Matrix types that are all derived from the base class MatrixBase. They do not have to store elements but can be logical or wrappers. They just have to provide the functions A.cols(), A.rows() (column and row numbers), A.mult(x)$$=A\cdot x$$ and A.transMult(y)$$=A^T\cdot y$$.

# We start off with the typical imports
import numpy as np
import pygimli as pg


## Dense matrix Matrix#

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.

A = pg.Matrix(3, 4)
A[0, 0] = 1
A[2, 2] = -1
A[1] = np.arange(4)
print(A)
x = np.arange(4) + 5
A*x

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 pg.matrix.CMatrix.

## Index-based sparse matrix SparseMapMatrix#

Sparse matrix (most elements are zero) with single access and. Typical for traveltime Jacobian (only certain cells covered by individual rays) or constraint matrices.

A = pg.SparseMapMatrix(2, 3)
A.setVal(0, 0, -1)
A.setVal(0, 1, +1)  # first-order derivative
A.setVal(1, 0, -1)
A.setVal(1, 1, +2)
A.setVal(1, 2, -1)  # second-order derivative
x = [10, 12, 13]
print(A * x)
ax, _ = pg.show(A)

2 [2.0, 1.0]


Exists also as complex-valued variant pg.matrix.CSparseMapMatrix.

## Column-compressed matrix SparseMatrix#

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 pg.matrix.CSparseMatrix.

## Diagonal matrices#

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.

A = pg.matrix.IdentityMatrix(3)  # , 2.0)
x = [10, 11, 12]
A*x

3 [10.0, 11.0, 12.0]


## Weighted matrices#

MultLeftMatrix/MultRightMatrix/MultLeftRightMatrix

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.

A = pg.Matrix(3, 4)
A += 1
w = [2, 3, -1]
B = pg.matrix.MultLeftMatrix(A, w)
x = np.arange(4)
print(A*x)
print(B*x)

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 (Mult2Matrix) or addition (Add2Matrix).

A = pg.Matrix(2, 3)
A += 1
B = pg.SparseMapMatrix(3, 4)
C = pg.matrix.Mult2Matrix(A, B)
x = np.arange(4)
C * x

2 [0.0, 0.0]


## Block matrices#

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 appear multiply.

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.


## Matrix combinations#

There are also simpler types of matrix combinations: * H2Matrix/V2Matrix: two matrices below/next to each other * HNMatrix/VNMatrix: one matrix repeated N times horizontally/vertically * NDMatrix: block diagonal matrix

## Matrix wrappers#

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

## Matrix generators#

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

$C^{-0.5}\cdot x=Q\cdot D^{-0.5}\cdot Q^T \cdot x$

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.

Gallery generated by Sphinx-Gallery