What is a matrix?

What is a matrix?

While we have CRTP base class DenseMatrix specifying the interface for matrices, this base is incomplete or wrong documented and even inconsistent. The specification is also not applied to all matrix-like classes neither in dune-common nor in dune-istl. So, I wanted to open a specific discussion on what is a matrix?

1. Matrix is a 2d container of values

On Wikipedia, a matrix is defined as "a rectangular array of numbers (or other mathematical objects)". Thus we have the following requirements:

• (a) A matrix stores elements/entries of a specific type. Seeing matrices as collections, this would mean we have a single element-type. For collections this type is typically denoted by value_type.
• "The size of a matrix is defined by the number of rows and columns that it contains". Actually, this definition specifies three size informations:
  1. (b) There is a type representing the sizes, typically denoted by size_type
  2. (c) There is number of rows and number of columns. In Dune this is (for historical reasons) denoted by N() and M() respectively.
  3. (d) There us the total size of a matrix. Since a matrix shape is denoted by n x m, I would suggest to define the size as the product n*m. This would allow to see single-row and single-column matrices as row-vector and column-vector, correspondingly. Compare, e.g., numpy.matrix.size, mtl::size(matrix), or Eigen::Matrix::size.
• (e) The elements of a matrix should be accessible by indices. Since the row and column indices are in the range [0,rows) and [0,cols), its index type should also be size_type.
• (f) Following the C-array syntax, elements are accessed by matrix[i][j].

2. Matrix is an element of a vector space

Matrices over a field build an algebraic set with vector-space properties:

• matrices for an additive abelian group (commutative group)
• matrices are scalable by scalars of the field type
• Multiplication with scalar is distributive with respect to addition.

These mathematical concepts result in the following requirements for the matrix type M:

• There is a field type associate to the matrix, i.e. Dune::FieldType<M>::field_type exists.
• There are arithmetic operations for matrices A and B:
  • (g) A + B, A - B
  • (h) A += B, A -= B
  • (i) -A, +A (optional)
• Multiplication with a scalar s of field type:
  • (j) A * s, s * A
  • (k) A *= s
  • (l) A / s, A /= s
• (m) Additionally, we typically require compound operations: A.axpy(s, B)

2.1. Normed space and inner-product space

By defining norms over matrices and inner products, the vector spcae could be turned into a normed vector space or an inner-product space

• Typical norm for matrices include:

  • (n) frobenius_norm()
  • (o) forbenius_norm2() (squared frobenius norm)
  • (p) infinity_norm()
  • (q) infinity_norm_real() (simplified infinity norm)

• Inner-products:

  • (r) frobenius inner-product A:B = tr(A^H A) (not implemented by any matrix)

3. Matrix is a linear map

A linear transformation or linear map between vector spaces can be represented by a matrix. The operations of a linear map are typically of the form of matrix-vector products:

• (s) A.mv(x, y): y = A*x
• (t) A.umv(x, y): y+= A*x
• (u) A.mmv(x, y): y-= A*x
• (v) A.usmv(alpha, x, y): y+= alpha * A*x

If the matrix additionally provide transposed operations:

• (w) A.mtv(x, y): y = A^T * x
• (x) A.umtv(x, y): y+= A^T * x
• (y) A.mmtv(x, y): y-= A^T * x
• (z) A.usmtv(alpha, x, y): y+= alpha * A^T * x

If the matrix additionally provide hermitian operations:

• (A) A.mhv(x, y): y = A^H * x (currently not implemented by any matrix)
• (B) A.umhv(x, y): y+= A^H * x
• (C) A.mmhv(x, y): y-= A^H * x
• (D) A.usmhv(alpha, x, y): y+= alpha * A^H * x

4. Other operations with matrices

• (E) Matrix multiplication: A in R^nxm and B in R^mxk then A * B in R^nxk
  • Matrices must be compatible
  • square matrices might form a ring under matrix addition and matrix multiplication this would allow A *= B
• (F) The transposed of a matrix: A in R^nxm -> A^T in R^mxn.
  • Either an operation A.transpose() or a free function transpose(A)?
  • The transposed matrix is itself a matrix
• (G) Trace of a matrix tr(A) = sum_i A_ii
  • Assumption: matrix is square
  • Either an operation A.trace() or a free function trace(A)?
  • Recursive definiton: tr(A) = sum_i tr(A_ii)
• (H) Determinant of a matrix
  • Assumption: matrix is square
  • Easy to define for small matrices
  • Generalization using, e.g., Laplace expansion
  • Special operation: gramian determinant
• (I) Inversion of a matrix or solution of a linear system
  • Assumption: matrix is square
  • easy to implement for small matrices
  • LU decomposition (or other forms of factorization)

5. Matrices can be traversed

Visiting all the (stored) elements of a matrix is a useful operation to generalize some algorithms

• Whats is a proper call for choosing the way/order of traversal?
• It exists begin() and end() functions.
  • Are not documented properly but typically mean traversal of rows.
• Additionally reverse traversal is given by beforeEnd() and beforeBegin()
  • Why not using rbegin() and rend() as for standard containers?
• Missing functions in the interface:
  • Traversal over rows first: e.g. begin_rows(), end_rows()
  • Traversal over columns first: e.g. begin_cols(), end_cols()

6. Special matrix types with more properties

• Sparse matrices:
  • Check whether elements exist in the store is required: A.exists(i,j)
  • Number of stored elements (nonzeros): A.nonzeroes()
  • Iterators only iterate over stored elements/nonzeros.
  • Often direct element-access not implemented (neither const nor mutable)
• Banded/Diagonal matrices:
  • Special type of sparse matrix
  • Additional access to diagonal(s): A.diagonal() (current interface does not allow return off-diagonals)
  • Direct access to diagonal entries: A.diagonal(i)
• Triangular matrices
• Tridiagonal matrices

Examples in dune-common

DenseMatrix

A DenseMatrix is the CRTB base class for some matrix implementation, but it lacks some of the properties from above. Also its documentation is partially wrong.

• (c) There are multiple functions returning number of rows/columns: A.N() and A.rows(), A.M() and A.cols(), where rows() -> mat_rows() and cols() -> mat_cols() call the implementation functions.
• (d) is defined differently. For some reason A.size() returns the number of rows.
• (g) arithmetic operations not completely available for matrices, +A not implemented (not required)
• (j) and (l) not implemented, except for A /= s
• (r) There is no inner product for matrices implemented
• (A) Missing lin-map operation A.mhv(x,y)
• (E) binary matrix-matrix product not implemented. For square matrices the mult-assign is implemented from left and right, i.e., A.leftmultiply(B) and A.rightmultiply(B). The general matrix-matrix multiplication is commented out.
• (G) Trace is not implemented.

FieldMatrix

Is an implementation of DenseMatrix thus inherits its properties. Some additional comments

• (c) Here rows and cols are enum constants hiding the interface functions rows() and cols()
• (g) FieldMatrix implements A + B and A - B
• (j) and (l) are completely implemented
• (E) Multiplication is implemented as A * B and leftmultiply[any]() and rightmultiply[any](). Additional operations are implemented in the namespace Dune::FMatrixHelp

DynamicMatrix

Is an implementation of DenseMatrix thus inherits its properties.

DiagonalMatrix

• (d) is defined differently. For some reason A.size() returns the number of rows.
• (f) The element access is only implemented for diagonal entries, i.e. matrix[i][j] throws if i != j
• (g), (j) Operations A + B, A - B, and -A not implemented.
• (i) compound assignment operators implemented for scalars: A += s and A -= s meaning elementwise addition and subtraction of that scalar.
• (m) Operation axpy not implemented.
• (r) There is no inner product for matrices implemented
• (A) Missing lin-map operation A.mhv(x,y)
• (E) No matrix-matrix products implemented
• (G) Trace is not implemented.

TransposedMatrixWrapper

The wrapper is created by transpose(matrix), but does not produce a matrix in any sense as described above. It has just one single operation: A * B^T. Note, this wrapper only works for transposing a FieldMatrix.

Test of the interface

There is a test in dune-common: checkmatrixinterface.hh but it does not test the complete matrix interface, but just the functions that were already implemented. Thus, this test is probably written after the implementation of the DenseMatrix class and not before.

An alternative test, that is implemented in terms of separating the various aspects of a matrix, is implemented in !953 (closed) in terms of c++20 concepts. These tests can run in the ci jobs for toolchains supporting this new c++ standard and can report interface failures.

Major problems with the current interface

The existing matrix implementations in dune-common (but also in dune-istl) only partially fulfil the interface definition from above and are inconsistent. Major problems are:

• Element access is not a defined property of a matrix. It is only guaranteed to succeed in some matrix implementations. This is mainly for performance reasons. This argument is a bit weak, because there are other ways to ensure performant data access, e.g. using iterators or specialized access functions.
• Traversal of matrices is not properly specified and cannot easily be extended to column-major traversal. begin() and end() have no knowledge about its traversal direction. There should be a better name for row-wise traversal and column-wise traversal. Proposal: begin_rows() and begin_cols(). Then, the user has to explicitly use a range-generator for the iteration, e.g., rows(matrix).begin() and columns(matrix).begin()
• The iterators itself are not properly specified. Is the dereferenced iterator again a range? What does it traverse?
• Some matrices implement a single index access matrix[r] as a row-view onto the matrix, i.e., this operation returns a vector-like container representing the rth row of the matrix. This is not clearly specified. Is this a required property of a matrix? What are the costs of this operation? Is the row a propery vector or just a proxy providing minimal access to the row entries?
• Especially in the implementation of geometries, one needs matrix-matrix multiplication. While this is implemented for FieldMatrix, all the other matrix types lack this operation. But, e.g., YaspGrid returns for its jacobian a DiagonalMatrix. Other grid implementations maybe something else. We need A * B, but also A^T * B and A * B^T for various combinations of matrices. This is partially implemented in some "Helper" namespaces in dune-geometry or for FieldMatrix in FMatrixHelper namespace. Multiplication seems not properly designed!
• Sparse matrices are a special type of matrices. The interface should mirror this mathematical property. Actually, sparse matrices are just matrices with a different storage type for its entries. Thus, sparse matrices should only extend the matrix interface with methods providing fast access to the stored elements. Either by iterators or something else. Currently, sparse matrices only implement a partial matrix interface and some methods even have different behaviour, e.g. matrix[i][j] is not guaranteed to succeed, not even for const-access, although the matrix value is clearly specified.
• Views over matrices, like transposed-view or hermitian-view are not properly implemented. They do not provide a matrix interface although this could easily be provided. One could discuss whether views are read-only or provide mutable access to the matrix elements as well.
• Currently, we can only test that a matrix supports all the different aspects as described above. Not the individual aspect. But, often an algorithm just needs some properties, e.g., in iterative krylov methods, one needs the linear-map properties for a matrix, for lu-decomposition maybe fast element access. In geometry parametrizations we need matrix-matrix multiplication and trace/Gramian determinant computation. It is required to have individual concept/interface checks.
• Single-row and Single-column matrices are not vectors. This often leads to problems, e.g., there is no two_norm() implemented and also dot() and axpy with another vector does not work, even if the shape would allow this. This problem is that vectors are not matrices and also not compatible in its interface to matrices.