Numerical Linear Algebra
Mathwrist White Paper Series

One fundamental built-in block of NPL is the Matrix class. We don’t introduce a separate data type for vectors. Instead, they are just 1-row or 1-column instances of the matrix class. Matrix row and column index starts from 0.

A matrix object is always associated with some underlying memory storage. The actual memory allocation scheme is handled by the library as implementation details. From a user perspective, it is sufficient to know that dynamic memory is allocated to store matrix elements. When a matrix object goes out of scope, its memory is freed automatically.

The standard use case of creating a $m\times n$ matrix is by calling the constructor Matrix::Matrix(m,n), which allocates a $(m\times n)$-element buffer in heap. In this case, matrix elements are contiguous double values placed in row-wise major. Internally, the library could have different memory representation scheme for special matrices. But only the standard use case is exposed to API level. Most of linear algebra functionality is implemented as member functions of the Matrix class in matrix.h.

Very often, we will need work on a part or a transformed representation of an original matrix. For efficiency purpose, we introduce the concept of matrix views that can be conveniently used just like regular matrix objects. A matrix view in this case shares the same data storage as the original matrix. Hence it is important to make sure the original matrix object is alive while using its views.

There are two types of matrix views, Matrix::View and Matrix::ViewConst both as nested classes of the Matrix class. From a constant matrix object, one can only obtain a view of Matrix::ViewConst type, which ensures the content of an original matrix cannot be modified by the view. On the other hand, a non-constant matrix could be modified by its views. The following code example illuatrates this.

Attention needs be paid when we copy matrix and view objects using the assignment operator. The behavior of this assignment operation depends on the actual types of its left hand side(LHS) and right hand side(RHS) parameters.

1. 1.

   LHS is a matrix: This matrix object is first resized to the same size as the RHS object and then an elemen-wise copy (deep copy) is performed.

2. 2.

   LHS is a non-constant view:

   1. (a)

      RHS is a matrix: Because a view has a predefined size, an element-wise copy happens if both sides have the same size. Otherwise, an exception is thrown.

   2. (b)

      RHS is also a non-constant view: No deep copy is involved. LHS becomes the same view as the RHS view object.

   3. (c)

      RHS is a constant view: RHS is implicitly converted to constant matrix object, then 2.a applies.

3. 3.

   LHS is a constant view: RHS can only be a constant view as well. In this case, both LHS and RHS become the same view on the original matrix. All other RHS types will cause a comilation error.

The following code example reflects the behavior of view assignment operations described above.

Matrix views can be passed as arguments anywhere when a matrix object is expected. We can also call a matrix member function through a view object using the operator -$>$ notation.

A special type of matrix views are matrix expressions, introduced as the following matrix nested classes that are all derived from Matrix::View.

• •

  Matrix::Negation, unary expression;

• •

  Matrix::Addition, binary expression;

• •

  Matrix::Substraction, binary expression;

• •

  Matrix::Multiplication, binary expression;

These expressions are designed with efficiency consideration so as to provide some syntax sugar for programmers. The C++ language ensures arithmetic precedence when multiple operators appear in one single expression. We can write cleaner code using matrix expressions illustrated by the code example below:

Matrix expressions are special in the sense that each expression object has its own memory allocated to store the elements of expression result. A regular view on the contrary just shares the same data storage as its original matrix source.

Matrix permutation reflects a sequence of row interchanges or column interchanges on a target matrix. Conceptually, we can think of permutation itself as a square matrix, where there is exactly one element equal to 1 in each row and column. All other elements are 0. When this permutation matrix left or right multiplies a regular matrix $M$, effectively row or column interchange are applied to $M$.

We represent matrix permutation as a nested class Matrix::Permutation in matrix.h. It can left multiply or right multiply a regular matrix object to generate a permutated view of the original matrix. The permutated view then can be used conveniently as a regular matrix.

Multiplication can also happen between two permutation objects. The end result is a new permutation object the reflects the accumulated row/column interchanges. Let $P_{i,j}$ and $P_{k,l}$ denote $(i,j)$ and $(k,l)$ interchange respectively. $P_{k,l}\times P_{i,j}\times M$ is to first perform $(i,j)$ row interchange on matrix $M$ and then followed by $(k,l)$ row interchange on the result of the first permutation. To perform the same sequence of interchanges on columns, we will use right multiplication $M\times P_{i,j}\times P_{k,l}$.

Finally, when a permutation is applied to a constant matrix, it results to a Matrix::ViewConst object, a constant view of the matrix as the code example below.

The elements of a matrix can be accessed via the Matrix::operator() member function. Each of such access is a separate function call, which could cause a lot of overhead if an application needs frequently access to matrix elements. In most cases though, we will be working with matrix objects that have elements contiguous in memory, i.e. directly on an original matrix object, or on a non-permutated view of the matrix. In these typical use cases, we can use matrix iterators for more efficient element access.

Matrix iterators are designed by the following mechanism inside the definition of Matrix class.

Essentially, an iterator support is a shared pointer of dense matrix traits that implements iterators by pointer arithmetic. When we obtain an iterator support from a non-permutated view or from a matrix directly, this shared pointer is not null. From that point on, we can simply use iterators to access matrix elements. Matrix iterator implementation consist of very light inline functions which are more efficient than Matrix::operator() member access function.

From a non-constant matrix we can obtain an iterator_support that allows us to both read and write the value of matrix elements. From a constant matrix, we obtain const_iterator_support that supports read access only. The code example below demonstrates how to use iterators.

Let $𝐀$ be the $n\times n$ matrix object that we want to compute its inverse ${𝐀}^{-1}$. The Matrix::inv() function takes an optional Matrix::Form parameter that helps to choose what algorithm to use, based on the form of $𝐀$. There are three special forms that are relevant in this context.

Given a $n\times n$ matrix $𝐀$, its eigen decomposition is computed in two phases. The first phase involves the order of $𝒪(n^3)$ number of floating point operations. The second phase generally requires $𝒪(n^3)$ floating point operations in order to converge to machine precision. However, when $𝐀$ is symmetric, the second phase can be computed much faster, even in the order of $𝒪(n^2)$ if only eigen values are required. Details of eigen algorithms can be found in chapter V in the book [1] by Lloyd N. Trefethen and David Bau. Our Matrix::eigen() function uses the Matrix::Form parameter to distinguish the symmetric case from the general case.

Given a generate $m\times n$ matrix $𝐀$, the svd() function computes a compact form of singular value decomposition $𝐀={\mathrm{𝐔𝐒𝐕}}^T$. Let $r=\min (m,n)$. $𝐔$ is resized to $m\times r$. $𝐒$ is resized to $r\times 1$ and $𝐕$ is resized to $n\times r$. The singular values in $𝐒$ are non-negative and arranged in descending order.

Given a $m\times n$ matrix $𝐀$, let $r=\min (m,n)$. The QR decomposition computes a reduced size $m\times r$ orthonormal basis matrix $𝐐$ and a $r\times n$ upper triangular matrix $𝐑$ such that $\mathrm{𝐐𝐑}=𝐀$. In this reduced version, when $m>n$, $𝐐$ consists of only the first $n$ orthonormal basis in columns. The full QR decomposition on the other hand computes a full $m\times m$ basis expansion matrix $𝐐$ using Householder reflector algorithm.

The reduced size QR and full QR have the same behavior on matrix $𝐑$. When $m>n$, $𝐑$ is a $n\times n$ upper triangular matrix. When , $𝐑$ is $m\times n$ an upper trapezoid matrix.

Our Cholesky factorization method is based on Algorithm 23.1 outlined in the book [1] by Lloyd N. Trefethen and David Bau. In addition, we also implemented a few other Cholesky decompositin methods of different flavors. They are used internally in our numerical optimization algorithms and could be relevant to other practical problems.

When , the linear system is underdetermined and has infinite number of solutions. The QR factorization method effectively produces a minimum norm solution of $𝐱$. We first compute , where $𝐘$ is the $n\times m$ range space and $𝐙$ is the $n\times (n-m)$ null space.

Since $𝐐$ is an orthonormal basis, we can write solution vector $𝐱={\mathrm{𝐘𝐱}}_Y+{\mathrm{𝐙𝐱}}_Z$ in terms of a range space component ${𝐱}_Y$ and null space component ${𝐱}_Z$. Clearly,

, where ${𝐱}_Y$ can be found uniquely by forward substitution on equation ${𝐑}^T{𝐱}_Y=𝐛$.

Different choice of ${𝐱}_Z$ produces different solution of $𝐱$. In practice, we let ${𝐱}_Z=0$ to obtain a minimum norm solution. If the QR factorization is known, we can directly use the range-space solver LinearSystem::range_space_solve().

When $m>n$, the linear system is overdetermined and doesn’t have an exact solution. When $𝐀$ is singular, the linear system either doesn’t have an exact solution or has infinite many solutions. The solver effectively computes the least square solution of $𝐱$ using SVD.

For any symmetric matrix $𝐀$, there exists a ${\mathrm{𝐋𝐃𝐋}}^T$ factorization where $𝐋$ is a lower triangular unit matrix and $𝐃$ is a diagonal matrix. Further, if $𝐀$ is positive definite, all elements in $𝐃$ are positive, $𝐀$ then has a Cholesky factor.

When $𝐀$ is a very large $n\times n$ positive definite matrix, computing the exact solution of $𝐱$ could be very time consuming. Instead, one might be satisfied with a “near” solution $\overline{𝐱}$ such that the residual $𝐫=𝐀\overline{𝐱}-𝐛$ is small enough.

In this situation, we can use conjugate gradient (CG) update method to generate a sequence of $({𝐱}_k,{𝐩}_k,{𝐫}_k)$. For some initial guess ${𝐱}_0$, we compute the residual vector ${𝐫}_0={\mathrm{𝐀𝐱}}_0-𝐛$. The initial search direction ${𝐩}_0$ is taken as ${𝐩}_0=-{𝐫}_0$. We always pass $({𝐱}_k,{𝐩}_k,{𝐫}_k)$ as input at the $k$-th iteration. Subsequently, the CG method updates ${𝐱}_{k+1}$ as improved “guess” with smaller residual ${𝐫}_{k+1}$. It also updates the search direction ${𝐩}_{k+1}$ that should be taken at the next iteration. Users can stop the iterations when ${𝐫}_k$ is considered small enough. In the worst case, the CG method ensures the exact soluton can be found after $n$ number of iterations. Our implementation is based on [2], chapter 5.

Suppose initially we know a LU decomposition of matrix ${𝐀}_0={𝐋}_0{𝐔}_0$. Subsequently, matrix ${𝐀}_i$ is obtained by pivoting the $j$-th column of ${𝐀}_{i-1}$ with a vector ${𝐩}_i$. In this situation, the LU decomposition of ${𝐀}_i$ can be efficiently computed by updating ${𝐋}_{i-1}$ and ${𝐔}_{i-1}$. LinearSystem::LU provides a linear system solver to solve ${𝐀}_i𝐱=𝐛$. Our implementation is based on [5].