Numerical Linear Algebra
Mathwrist C++ API Series

(January 20, 2023)

Abstract

This document details the C++ programming interface of the numerical linear algebra (LA) features implemented in Mathwrist’s Numerical Programming Library (NPL). For a technical overview of algorithms, please refer to our white paper series.

1 Matrix Basics

The Matrix class is one of the most important built-in blocks of the library. The majority part of linear algebra implementation is encapsulated as member functions of the Matrix class defined in C++ header file matrix.h.

1.1 Enumeration of Matrix Forms

Matrix class defines the following enumeration of matrix forms. These forms help to choose appropriate algorithms in certain context.

1 enum Form

2 {

3 SYMMETRIC, /**< symmetric matrix */

4 POS_DEFINITE, /**< symmetric positive definite */

5 POS_SEMI, /**< symmetric positive semi-definite */

6 INDEFINITE, /**< symmetric indefinite */

7 SINGULAR, /**< matrix is possibly singular */

8 GENERAL /**< general matrix non-singular */

9 };

1.2 Object Creation

A $m\times n$ matrix object created by the following constructors has contiguous memory allocated in heap to store elements. All matrix elements are initialized to 0 by the constructor. We shall call such a matrix object as “physically allocated” matrix. Later on in the view section, we will introduce “view represented” matrix objects. No matter what is the actual memory allocation scheme, when a matrix object goes out of scope, its internal resouces are freed by the destructor.

1 Matrix();

2 Matrix(size_t nRows, size_t nCols);

3 Matrix(size_t n);

4 Matrix(const Matrix& other);

5 ~Matrix();

Description:

1.

Default contructor to create an empty matrix object.
2.

Create a matrix object of nRows number of rows and nCols number of columns;
3.

Create a $n\times n$ square matrix;
4.

Copy constructor, deep copy all elements from the other matrix.
5.

Destructor, free all resouces used by the matrix object.

1.3 Static Initialization Methods

1 static void identity (size_t n, Matrix& I);

2 static void diagonal(size_t n, double lambda, Matrix& D);

3 static void diagonal(const Matrix& dv, Matrix& D,

4 bool inverse = false);

Description:

1.

Initialize matrix I to an identity matrix of size n.
2.

Diagonize a matrix D with a given constant factor lambda.
3.

Diagonize a matrix D with a given vector dv. If parameter inverse is passed true, the inverse of vector element of dv are assigned to the diagonal of D.

1.4 Element Access

All functions listed below are Matrix class memeber functions. Matrix row and column index starts from 0.

1 double& operator() (size_t i, size_t j);

2 const double& operator() (size_t i, size_t j) const;

Description:

1.

Return a reference to the element at the $i$ -th row and $j$ -th column in the matrix. Element can be modified.
2.

For read only, return a constant reference to the element at the $i$ -th row and $j$ -th column in the matrix.

1.5 Dimension and Size Query

All functions listed below are Matrix class memeber functions.

1 bool has_dimension(size_t n_rows, size_t n_cols) const;

2 bool has_same_dim (const Matrix& other) const;

3 bool is_defined() const;

4 bool is_square() const;

5 bool is_symmetric() const;

6 bool is_row_vector() const;

7 bool is_col_vector() const;

8 bool is_vector() const;

9 size_t row_size () const;

10 size_t col_size () const;

11 size_t size () const;

Description:

In the same order of the above code block, the function returns,

1.

if the matrix has n_rows number of rows and n_cols number of columns.
2.

if the matrix has the same dimension size as the other matrix.
3.

if the matrix is empty.
4.

if the matrix is a square matrix.
5.

if the matrix is a symmetric matrix.
6.

if the matrix is a row vector.
7.

if the matrix is a column vector.
8.

if the matrix is a row or column vector.
9.

the total number of rows in the matrix.
10.

the total number of columns in the matrix.
11.

the total number of elements in the matrix.

1.6 Comparison and Content Query

All functions listed below are Matrix class memeber functions.

1 bool operator==(const Matrix& other) const;

2 bool equals_to(const Matrix& other, double eps=EPSILON()) const;

3 bool equals_to(const Matrix& other, const Matrix& eps) const;

4 bool is_zero(double eps=EPSILON()) const;

5 bool is_zero(const Matrix& eps) const;

6 bool is_diag(double eps=EPSILON()) const;

7 bool is_lower_triangular(double eps=EPSILON()) const;

8 bool is_upper_triangular(double eps=EPSILON()) const;

9 bool is_identity(double eps=EPSILON()) const;

Description:

In the same order of the above code block, based on numerical comparison with tolerance $\epsilon$ (default to $\epsilon=1.e^{-12}$ ), the function returns if the current matrix object numerically is,

1.

equal to the other matrix.
2.

equal to the other matrix with a single tolerance eps
3.

equal to the other matrix with an element-wise tolerance matrix eps.
4.

equal to 0, with a single tolerance eps.
5.

equal to 0, with an element-wise tolerance matrix eps.
6.

a diagonal matrix, with a single tolerance eps.
7.

a lower triangular matrix, with a single tolerance eps.
8.

a upper triangular matrix, with a single tolerance eps.
9.

an identity matrix, with a single tolerance eps.

1.7 Assignment Operators

All functions listed below are Matrix class memeber functions.

1 Matrix& operator= (const Matrix& other);

2 Matrix& operator= (double v);

Description:

1.

Matrix element-wise copy $\mathbf{A}=\mathbf{B}$ , where $\mathbf{A}$ is the current matrix object, $\mathbf{B}$ is passed by reference from parameter other. If $\mathbf{A}$ is actually a view representation, it must have the same dimension size as $\mathbf{B}$ . Otherwise, an exception will be thrown. If $\mathbf{A}$ is a physically allocated matrix object and has difference size from $\mathbf{B}$ , it will be resized before copying.
2.

Set all elements in the current matrix to value parameter v.

1.8 Manipulation Functions

All functions listed below are Matrix class memeber functions.

1 void resize (size_t nRows, size_t nCols, bool to_zeros = true);

2 void element_copy(const Matrix& other, int sign = 1);

3 void zeros();

4 void ones();

5 void e(size_t i);

6 void negate();

7 void diagonalize (double lambda);

8 void diagonalize (const Matrix& diag);

9 void diagonalize (const ftl::Buffer<double>& diag);

10 void triangular_clean (bool clean_lower);

11 void symmetric_copy(bool upper_to_lower = true);

12 void diagonal_copy(const Matrix& other);

13 void transpose();

14 void mean_centered();

15 void swap_col(size_t i, size_t j);

16 void swap_row(size_t i, size_t j);

17 void shift_col_to_right();

18 void shift_col_to_left();

19 void shift_row_to_top();

20 void shift_row_to_bottom();

Description:

Let $\mathbf{A}$ denote the current matrix object. In the same order in the above code block,

1.

Change the size of $\mathbf{A}$ to nRows number of rows and nCols number of columns. If the parameter to_zeros flag is passed to true, all elements of $\mathbf{A}$ are set to 0 after resize. This function involves memory allocation if $\mathbf{A}$ doesn’t have the desired size. Hence, if $\mathbf{A}$ is actually a view representation, an exception is thrown when $\mathbf{A}$ has different size.
2.

REQUIRED: Matrix sizes are identical. It element-wise copies all values from the other matrix (passed by reference) to $\mathbf{A}$ . This function is mostly useful for views because the assignment operator for views doesn’t copy elements, instead view assignment creates a view copy.
3.

Sets all elements in $\mathbf{A}$ equal to 0.
4.

Sets all elements in $\mathbf{A}$ equal to 1.
5.

REQUIRED: $\mathbf{A}$ is a row vector or a column vector. It sets the $i$ -th element of $\mathbf{A}$ to 1 and all other elements to 0, i.e. $\mathbf{A}=e_{i}$ .
6.

Flips the sign of all elements in $\mathbf{A}$ .
7.

REQUIRED: $\mathbf{A}$ is a square matrix. It sets $\mathbf{A}$ to be a diagonal matrix where all diagonals are equal to the given parameter value lambda.
8.

REQUIRED: $\mathbf{A}$ is a square matrix. It sets $\mathbf{A}$ to be a diagonal matrix where the diagonal elements are passed as a vector from parameter diag.
9.

REQUIRED: $\mathbf{A}$ is a square matrix. It sets $\mathbf{A}$ to be a diagonal matrix where the diagonal elements are passed as a buffer from parameter diag.
10.

REQUIRED: $\mathbf{A}$ is a square matrix. It sets all elements below diagonal to 0 when flag clean_lower is passed as true. When it is false, all elements above diagonal are set to 0.
11.

REQUIRED: $\mathbf{A}$ is a square matrix. When flag upper_to_lower is passed as true, it copies all above diagonal elements to below diagonal elements so that the matrix is symmetric. When flase is passed, it copies all below diagonal elements to above diagonal.
12.

REQUIRED: $\mathbf{A}$ is a square matrix. The other matrix is required to be either a square matrix of the same size, or a vector of same row or column size. When other is a matrix, its diagonal elements are copied to the diagonal of $\mathbf{A}$ . When it is a vector, its element are copied to the diagonal of $\mathbf{A}$ .
13.

Change the current matrix to its transpose.
14.

Let $\mathbf{a}_{j}$ denote the $j$ -th column of $\mathbf{A}$ . This function computes the mean of $\mathbf{a}_{j},\forall j$ and then subtract all elements of $\mathbf{a}_{j}$ by this mean.
15.

Swaps the $i$ -th column of $\mathbf{A}$ with the $j$ -th column.
16.

Swaps the $i$ -th row of $\mathbf{A}$ with the $j$ -th row.
17.

Shift all columns of $\mathbf{A}$ to the right by 1 column. The content of the original first column is not erased.
18.

Shift all columns of $\mathbf{A}$ to the left by 1 column. The content of the original last column is not erased.
19.

Shift all rows to the top by 1 row. The content of the original last (bottom) row is not erased.
20.

Shift all rows to the bottom by 1 row. The content of the original first (top) row is not erased.

1.9 Norm, Min and Max Functions

All functions listed below are Matrix class memeber functions.

1 double inner_product(const Matrix& b) const;

2 double vector_norm_1() const;

3 double vector_norm_2() const;

4 double vector_norm_max(size_t* _i_max = 0) const;

5 double vector_min(size_t* _i = 0) const;

6 double vector_max(size_t* _i = 0) const;

7 double vector_min_abs(size_t* _i = 0) const;

8 double vector_sum() const;

9 double vector_max_diff(const Matrix& b) const;

Description:

Let $\mathbf{a}$ denote the current matrix object. The above function requires $\mathbf{a}$ to be a row or column vector of size $n$ . In the same order of the above code block,

1.

REQUIRED: Parameter $\mathbf{b}$ is a vector object of size $n$ . The function computes the inner product of $\mathbf{a}$ and $\mathbf{b}$ .
2.

The function computes the vector 1-norm $\|\mathbf{a}\|_{1}$ .
3.

The function computes the vector 2-norm $\|\mathbf{a}\|_{2}$ .
4.

The function computes the vector $\infty$ -norm $\|\mathbf{a}\|_{\infty}$ . If the pointer parameter i_max is not null, *i_max holds the index of the max norm element on return.
5.

The function returns the min element in the vector. If the pointer parameter _i is not null, *_i holds the index of the min element on return.
6.

The function returns the max element in the vector. If the pointer parameter _i is not null, *_i holds the index of the max element on return.
7.

The function returns the minimum of $|\mathbf{a}_{i}|,i=0,\cdots,n-1$ . If the pointer parameter _i is not null, *_i holds the index of this element on return.
8.

The functions returns the sum of all elements.
9.

The function returns the max of $|\mathbf{a}_{i}-\mathbf{b}_{i}|,i=0,\cdots,n-1$ .

2 Matrix Arithmetics

All functions listed in this section are Matrix class memeber functions.

2.1 Arithmetic Operators

1 Matrix& operator+= (const Matrix& other);

2 Matrix& operator+= (double v);

3 Matrix& operator-= (const Matrix& other);

4 Matrix& operator-= (double v);

5 Matrix& operator*= (const Matrix& other);

6 Matrix& operator*= (double v);

Let $(\cdot)$ denote the relevant arithmetic operation. The functions listed above perform $\mathbf{A}=\mathbf{A}(\cdot)\mathbf{B}$ operation on the current matrix object $\mathbf{A}$ with operand $\mathbf{B}$ and returns a reference to $\mathbf{A}$ . An exception will be thrown if $\mathbf{A}$ and $\mathbf{B}$ have incompatible size.

Description:

1.

$(\cdot)$ is addition, $\mathbf{B}$ is passed by reference from parameter other.
2.

$(\cdot)$ is addition, $\mathbf{B}$ has all element value being parameter v.
3.

$(\cdot)$ is subtraction, $\mathbf{B}$ is passed by reference from parameter other.
4.

$(\cdot)$ is subtraction, $\mathbf{B}$ has all element value being parameter v.
5.

$(\cdot)$ is multiplication, $\mathbf{B}$ is passed by reference from parameter other.
6.

$(\cdot)$ is multiplication, $\mathbf{B}$ has all element value being parameter v.

2.2 Element-wise Multiplication and Division

1 void dot_times(const Matrix& other);

2 void dot_divide(const Matrix& other);

Description:

Mostly used for vectorized arithmetics, these two function perform element-wise operation $\mathbf{A}(i,j)=\mathbf{A}(i,j)(\cdot)\mathbf{B}(i,j),\forall i,j$ , where $\mathbf{A}$ is the current matrix and $\mathbf{B}$ is passed by reference from parameter other. The operation $(\cdot)$ here is multiplication and division respectively. If $\mathbf{A}$ and $\mathbf{B}$ have different size, an exception is thrown.

2.3 Matrix Multiplication and Inverse

1 void multiply(const Matrix& B, Matrix& C) const;

2 void AtA(Matrix& C) const;

3 void inv(Matrix& Inv, Matrix::Form f = GENERAL) const;

Description:

Let $\mathbf{A}$ denote the current $m\times n$ matrix object.

1.

computes $\mathbf{C}=\mathbf{A}\times\mathbf{B}$ , where the operand matrix $\mathbf{B}$ is required to have size $n\times l$ . The result matrix $\mathbf{C}$ is resized to be a $m\times l$ matrix.
2.

computes $\mathbf{C}=\mathbf{A}^{T}\times\mathbf{A}$ , where the result matrix $\mathbf{C}$ is resized to be a $n\times n$ matrix.
3.

computes the inverse of the current matrix, where $m=n$ .

Parameters
1. (a)
  
  Inv will be resized to hold the inverse matrix on return.
2. (b)
  
  f indicates what is the form of the current matrix object, which helps to choose appropriate algorithms to compute the inverse. Please refer to our white paper for the algorithm details based on different matrix forms.

3 Transformation Matrices

Linear transformations mathematically can be written as matrix multiplications. Programmatically, transformation matrices often deserve special treatment for both storage and performance efficiency. We gives details on two special transfromations, permutation and Givens rotaion, in this section.

3.1 Permutation

Matrix::Permuation is a nested class of the Matrix class. Conceptually, it represents a $n\times n$ square matrix $\mathbf{P}$ where each row and column has exactly one element equal to 1. All other element are 0s. This allows us to only store the position of those 1-element in order to represent the matrix.

Given a matrix $\mathbf{M}$ , the left multiplication $\mathbf{P}\times\mathbf{M}$ generates a row interchanged view of $\mathbf{M}$ . The right multiplication $\mathbf{M}\times\mathbf{P}$ generates a column interchanged view of $\mathbf{M}$ . We will discuss matrix views in a separate section.

3.1.1 Constructors

1 Permutation();

2 Permutation(size_t n, size_t i = 0, size_t j = 0);

3 Permutation(const ftl::Buffer<size_t>& encoding);

4 Permutation(const size_t* _encoding, size_t n);

Description:

In the same order of the above code block,

1.

Creates an empty permutataion object.
2.

Creates a permutation object that represents $i$ and $j$ row or column interchange. Parameter $n$ defines the size of the permutation matrix.
3.

Creates a permutation object $\mathbf{P}$ where element $\mathbf{P}(i,\mbox{encoding}[i])=1$ . The encoding buffer has size $n$ , which also defines the size of the permutation matrix.
4.

Same as the constructor that takes an encoding buffer as argument, this function takes a pointer and size $n$ parameter instead.

3.1.2 Two Useful Factory Methods:

The functions listed below are Matrix::Permutation class static member functions.

1 static Permutation i_shift_row(size_t n, size_t i, size_t j);

2 static Permutation i_shift_col(size_t n, size_t i, size_t j);

Description:

These two functions can be best explained by examples. Let $\mathbf{M}$ be a $6\times 6$ matrix. Let $\mathbf{m}_{k}$ denote the $k$ -th row of $\mathbf{M}$ , $k=0,\cdots,5$ . In other words,

\mathbf{M}=\begin{pmatrix}\mathbf{m}_{0}\cr\mathbf{m}_{1}\cr\mathbf{m}_{2}\cr% \mathbf{m}_{3}\cr\mathbf{m}_{4}\cr\mathbf{m}_{5}\end{pmatrix}

If we want to shift rows $\mathbf{m}_{k=2,3,4}$ one row up and move $\mathbf{m}_{1}$ to $\mathbf{m}_{4}$ ’s original position, it involves a sequence of row interchanges. Let $\mathbf{P}$ be the permutation object that performs this sequence of interchanges.

\mathbf{P}\mathbf{M}=\begin{pmatrix}\mathbf{m}_{0}\cr\mathbf{m}_{2}\cr\mathbf{% m}_{3}\cr\mathbf{m}_{4}\cr\mathbf{m}_{1}\cr\mathbf{m}_{5}\end{pmatrix}

This permutation object $\mathbf{P}$ can be obtained by calling the i_shift_row() factory method with arguments i=1 and j=4 as the code example below,

1 // m_1 moves to m_4’s position. m_2, m_3, m_4 shift up.

2 Permutation P = Permutation::i_shift_row(6, 1, 4);

If we want to shift rows $\mathbf{m}_{k=1,2,3}$ one row down and move $\mathbf{m}_{4}$ to $\mathbf{m}_{1}$ ’s position such that

\mathbf{P}\mathbf{M}=\begin{pmatrix}\mathbf{m}_{0}\cr\mathbf{m}_{4}\cr\mathbf{% m}_{1}\cr\mathbf{m}_{2}\cr\mathbf{m}_{3}\cr\mathbf{m}_{5}\end{pmatrix}

, this $\mathbf{P}$ can be obtained by

1 // m_4 moves to m_1’s position. m_1, m_2, m_3 shift down.

2 Permutation P = Permutation::i_shift_row(6, 4, 1);

Effectively, when index $i<j$ , i_shift_row shifts all rows with index $h,i<h\leq j$ one row up. When $i>j$ , the shift is one row down. In both cases, the original $i$ -th row is simultaneously moved to the $j$ -th row.

The i_shift_col() function behaves in a similar way for column interchanges. Now let $\mathbf{m}_{k}$ be the columns of $\mathbf{M}$ . When parameter $i<j$ , columns $\mathbf{m}_{k=i+1,\cdots,j}$ shift to left by one column and column $\mathbf{m}_{i}$ moves to the original position of $\mathbf{m}_{j}$ . When $i>j$ , columns $\mathbf{m}_{k=j,\cdots,i-1}$ shift to right by one column and $\mathbf{m}_{i}$ moves to the original position of $\mathbf{m}_{j}$ .

When parameter $i=j$ , these two function both return a permutation corresponding to an identity matrix.

3.1.3 Multiplications

All functions listed below, except the 6th one, are Matrix::Permutation class member functions. The 6th function is a friend binary operator where the permutation object appears as the right hand side operand of the multiplication.

1 Permutation Permutation::operator*(const Permutation& Q) const;

2 Permutation& Permutation::operator*=(const Permutation& Q);

3 View Permutation::operator*(Matrix& M) const;

4 ViewConst Permutation::operator*(const Matrix& M) const;

5 size_t Permutation::operator*(size_t i) const;

6 friend size_t operator*(size_t j, const Permutation& P);

7 View Permutation::PMPt(Matrix& M) const;

8 ViewConst Permutation::PMPt(const Matrix& M) const;

9 View Permutation::PtMP(Matrix& M) const;

10 ViewConst Permutation::PtMP(const Matrix& M) const;

Description:

Let $\mathbf{P}$ and $\mathbf{Q}$ be two permutation matrices of $n$ rows. Let $\mathbf{M}$ be a matrix under the permutation operation. In the same order of the above code block,

1.

Returns a new permutation $\hat{\mathbf{P}}=\mathbf{P}\times\mathbf{Q}$ .

Generalize this to $\hat{\mathbf{P}}=\mathbf{P}_{n}\cdots\mathbf{P}_{1}\mathbf{P}_{0}$ , where each $\mathbf{P}_{i},i=0,\cdots,n$ is a single interchange operation. $\hat{\mathbf{P}}\times\mathbf{M}$ is equivalent to applying the row interchange sequence $i=0,\cdots,n$ . $\mathbf{M}\times\hat{\mathbf{P}}$ is to apply the column interchanges in reverse order.
2.

Performs $\mathbf{P}=\mathbf{P}\times\mathbf{Q}$ .
3.

REQUIRED: $\mathbf{M}$ has $n$ number of rows. This function returns a row-interchanged view of $\mathbf{M}$ , where the row interchanges are represented by $\mathbf{P}$ . Please refer to the Matrix class member multiplication operators, i.e. Matrix::operator*(Permutation&) for column-interchanged permutation.
4.

Same as the previous function except that now $\mathbf{P}$ is applied to a constant matrix and returns a constant view.
5.

Return the final position of the $i$ -th row after the permutation. This is when $\mathbf{P}$ appears as the left hand side multiplication operand.
6.

Return the final position of the $j$ -th column after the permutation. This is when $\mathbf{P}$ appears as the right hand side multiplication operand.
7.

REQUIRED: $\mathbf{M}$ is $n\times n$ square matrix. The function performs symmetric row and column interchanges on $\mathbf{M}$ and returns a vew that represents the multiplication $\mathbf{P}\mathbf{M}\mathbf{P}^{T}$ .
8.

Same as the previous function exception now $\mathbf{M}$ is a constant matrix and the function returns a constant view.
9.

REQUIRED: $\mathbf{M}$ is $n\times n$ square matrix. The function performs symmetric row and column interchanges on $\mathbf{M}$ and returns a vew that represents the multiplication $\mathbf{P}^{T}\mathbf{M}\mathbf{P}$ .
10.

Same as the previous function exception now $\mathbf{M}$ is a constant matrix and the function returns a constant view.

3.1.4 Query Functions:

The functions listed below are Matrix::Permutation class member functions.

1 bool operator== (const Permutation& other) const;

2 size_t size() const;

Description:

Let $\mathbf{P}$ be the current $n\times n$ permutation object. In the same order of the above code block,

1.

Compares if $\mathbf{P}$ effectively performs the same operation as the other permutation.
2.

Return the size $n$ .

3.1.5 Manipulation Functions:

The functions listed below are Matrix::Permutation class member functions.

1 void resize(size_t n);

2 void reset();

3 void transpose();

4 Permutation get_transpose() const;

Description:

Let $\mathbf{P}$ be the current $m\times m$ permutation object. In the same order of the above code block,

1.

Change $\mathbf{P}$ to size $n\times n$ and initialize it to be effectively an identity matrix.
2.

Reset $\mathbf{P}$ to be effectively an identity matrix of the same size.
3.

Change $\mathbf{P}$ to be $\mathbf{P}^{T}$ . For a general $\mathbf{P}=\mathbf{P}_{n}\cdots\mathbf{P}_{1}\mathbf{P}_{0}$ where each $\mathbf{P}_{i},i=0,\cdots,n$ is a single interchange operation, $\mathbf{P}^{T}=\mathbf{P}_{0}\mathbf{P}1\cdots\mathbf{P}_{n}$ . The transpose effectively reverses the order of the sequence.
4.

Return a new permutation object that is $\mathbf{P}^{T}$ .

3.2 Givens Rotation

A $n\times n$ Givens rotation matrix $\mathbf{G}$ is a unitary matrix with two special index $i$ and $j$ . It has the following elements,

•

All diagonal elements $\mathbf{G}_{k,k}=1,\forall k\neq i,j$ ;
•

Diagonal elements $\mathbf{G}_{i,i}=\mathbf{G}_{j,j}=c$
•

All off-diagonal elements are 0 except $\mathbf{G}_{i,j}=-s$ and $\mathbf{G}_{j,i}=s$ ;

, where $c=\cos\theta$ and $s=\sin\theta$ for some $\theta$ . When $\mathbf{G}$ left multiplies a column vector $\mathbf{x}$ . Only the $i$ -th and $j$ -th element of the result vector $\bar{\mathbf{x}}=\mathbf{G}\mathbf{x}$ are different from $\mathbf{x}$ .

We represent the Givens rotation matrix by nested class Matrix::Givens defined in header file “matrix.h”.

3.2.1 Constructors:

1 Givens();

2 Givens(size_t i, size_t j, double a, double b);

Description:

1.

Default constructor, create an empty object.
2.

Initialize the rotation object with index $i$ , $j$ , $\sin\theta=-\frac{b}{r}$ and $\cos\theta=\frac{a}{r}$ , where $r=\sqrt{a^{2}+b^{2}}$ . For a column vector $\mathbf{x}$ where $\mathbf{x}_{i}=a$ and $\mathbf{x}_{j}=b$ . This construction of $\mathbf{G}$ sets $\bar{\mathbf{x}}_{i}=r$ and $\bar{\mathbf{x}}_{j}=0$ in the result vector $\bar{\mathbf{x}}=\mathbf{G}\mathbf{x}$ .

3.2.2 Rotation Operation:

The functions listed below are Givens class member functions.

1 Matrix& operator* (Matrix& M) const;

2 void apply(double& a, double& b) const;

Description:

Let $\mathbf{G}$ be the current rotation object.

1.

The left multiplication operator* takes operand mathrix $\mathbf{M}$ and applies the operation $\mathbf{M}=\mathbf{G}\mathbf{M}$ on all columns of $\mathbf{M}$ . After the rotation, a reference of $\mathbf{M}$ is returned. The right multiplication operation $\mathbf{M}=\mathbf{M}\mathbf{G}$ is implemented as a class member function of the Matrix,

Matrix& Matrix::operator*= (const Givens& G);

It modifies the $i$ -th and $j$ -th column elements on all rows of $\mathbf{M}$ and returns a reference of $\mathbf{M}$ after rotation.
2.

The apply function conceptually performs the operation $\mathbf{x}=\mathbf{G}\mathbf{x}$ on vector $\mathbf{x}$ . Since only the $i$ -th and $j$ -th element of $\mathbf{x}$ are changed, $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ are passed by reference at parameter a and b. $\mathbf{G}$ then directly modifies the values of $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ .

3.2.3 Query Functions:

The functions listed below are Givens class member functions.

1 size_t index_a() const;

2 size_t index_b() const;

3 double r() const;

4 double c() const;

5 double s() const;

Description:

These query functions correspond to the constructor that takes $(i,j,a,b)$ as arguments to construct the rotation $\mathbf{G}$ .

1.

Return the index $i$ .
2.

Return the index $j$ .
3.

Return $r=\sqrt{a^{2}+b^{2}}$ .
4.

Return $c=\frac{a}{r}$ .
5.

Return $s=-\frac{b}{r}$ .

3.2.4 Manipulation Functions:

The functions listed below are Givens class member functions.

1 void increment_shift(size_t m);

2 void decrement_shift(size_t m);

Description:

These two funtions are useful when applying the same rotation to vector elements $\mathbf{x}_{k}$ and $\mathbf{x}_{l}$ where $k$ and $l$ are offsetted from $i$ and $j$ .

1.

Increase the index $i$ and $j$ by $m$ , i.e. $k=i+m$ and $l=j+m$ .
2.

Decrease the index $i$ and $j$ by $m$ , i.e. $k=i-m$ and $l=j-m$ .

4 Matrix Views

Matrix views are alternative representation of matrices in Mathwrist NPL. They are defined as nested classes of the Matrix class. A view object can be used anywhere in a context where a regular matrix object is expected, i.e. passed as function parameters or expression operands. In addition, a view object provides pointer notation so that it can be used as the pointer of a matrix object to access the Matrix class memeber functions.

Views can be copied by the copy constructor and assignment operator. The copy operation on views effectively creates indentical views instead of copying data elements. There are two view types, Matrix::ViewConst and Matrix::View respectively. A constant view provides read-only access to matrix elements. A non-constant view supports both read and write access to matrix elements. When a modification is applied to the non-constant view, it is the original matrix being modified through the view.

4.1 Constructors:

1 ViewConst();

2 ViewConst(size_t nRows, size_t nCols, const double* buffer);

3 View();

4 View(size_t nRows, size_t nCols, double* buffer);

Description:

1.

Creates an empty constant view.
2.

Creates a constant view of nRows by nCols matrix where the data elements are provided in a constant buffer in row-wise major.
3.

Creates an empty non-constant view.
4.

Creates a non-constant view of nRows by nCols matrix where the data elements are read from or written to a given buffer in row-wise major.

Example:

1 // Create a 3 x 3 matrix view from the given data buffer.

2 const double _DATA[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9 };

3 Matrix::ViewConst M(3, 3, _DATA);

5 // Create a 3 x 3 identity (physically allocated) matrix.

6 Matrix I;

7 Matrix::identity(3, I);

9 // M passed as a regular matrix argument.

10 Matrix N;

11 I.multiply(M, N);

12 assert(N.equals_to(M));

14 // Pointer notation to access Matrix member funciton.

15 M->multiply(I, N);

16 assert(M->equals_to(N));

18 // Access data elements same way as a matrix object

19 std::cout << M(0,1) << std::endl; // output is 2

4.2 Conversion Operators:

1 ViewConst::operator const Matrix& () const;

2 View::operator const Matrix& () const;

3 View::operator Matrix& ();

Description:

1.

In the context where a constant matrix is expected, convert a constant view to a constant matrix, i.e. line 11 and 12 in the code example above.
2.

In the context where a constant matrix is expected, convert a non-constant view to a constant matrix.
3.

In the context where a non-constant matrix is expected, convert a non-constant view to a non-constant matrix.

4.3 Pointer Operators:

1 const Matrix* ViewConst::operator->() const;

2 const Matrix* View::operator->() const;

3 Matrix* View::operator->();

Description:

1.

Allows a constant view to call Matrix class constant member functions, i.e. line 15 and 16 in the code example above.
2.

Allows a non-constant view to call Matrix class constant member functions.
3.

Allows a non-constant view to call Matrix class non-constant member functions.

4.4 Assignment Operators:

1 void View::operator= (const Matrix& m);

2 void View::operator= (double v);

Description:

These two functions are the member funcitons of only Matrix::View class. They copy elements from a source to the current view. If the current view mapps to an original matrix, the content of the original matrix is modified.

1.

REQUIRED: The current view has same size as the source matrix ${\color[rgb]{0,0,1}\small{m}}$ . This function copies all elements from m to the current view.
2.

This function sets all elements in the current view equal to the parameter value v.

Example:

1 Matrix m(3, 3);

2 Matrix n(1, 2);

4 n(0, 0) = 0.2;

5 n(0, 1) = 0.3;

7 // v is a view of the leading 1 x 2 sub matrix of m.

8 Matrix::View v = m.sub(0, 0, 1, 2);

10 // Copy elements from source matrix n to the

11 // sub matrix view of m.

12 v = n;

14 // It is the original matrix m being modified.

15 std::cout << m(0,0) << std::endl; // output 0.2

16 std::cout << m(0,1) << std::endl; // output 0.3

The above code example illustrates the assignment operation to a view object v. Note that if view v has different size from matrix n, the assignment at line# 12 will trigger an exception.

4.5 Copy Operation:

The copy constructor and assignment operator from views to views are not implemented in the library. It is the compiler’s default version (shallow copy) that is used when we copy two objects of the same view type. In which case we obtain two identical views. The code example below illustrates this.

Example:

1 Matrix m(3, 3);

2 Matrix n(2, 2);

4 // v is the leading 2 x 2 sub matrix view of m.

5 Matrix::View v = m.sub(0, 0, 2, 2);

7 // u is a view of the 1st row of n.

8 Matrix::View u = n.row(0);

10 // Now, v is a view of the 1st row of n as well.

11 // u and v are identical views.

12 v = u;

Seconly, because View and ViewConst are two different types, copy from View to ViewConst will cause compile time error. On the other hand, copy from ViewCont to View can be resolved by C++ standard implicit type conversion rules. First, the ViewConst class has a conversion operator that can take a ViewConst object to anywhere const Matrix& is expected. Secondly, the View class has an assignment operator to perform element-wise copy from a const Matrix& matrix source. Therefore, copy from ViewCont to View involves an implicit conversion and then an element wise copy operation. The code example below illustrates this.

Example:

1 // Make m a 2 x 2 identity matrix.

2 Matrix m(2,2);

3 m.diagonalize(1);

5 const Matrix& mr = m;

7 // r0 is a constant view of the 1st row of m.

8 Matrix::ViewConst r0 = mr.row(0);

10 // r1 is a non-constant view of the 2nd row of m.

11 Matrix::View r1 = m.row(1);

13 // Copying from a constant view to a non-constant view

14 // involves:

15 // 1) r0 implicitly converted to a matrix object.

16 // 2) element wise copy from r0 to r1.

17 r1 = r0;

19 // Now m’s 2nd row should be same as the 1st row.

20 ftl::Assert(m(1, 0) == 1);

21 ftl::Assert(m(1, 1) == 0);

4.6 Element Access:

1 double ViewConst::operator() (size_t i, size_t j) const;

2 double View::operator() (size_t i, size_t j) const;

3 double& View::operator() (size_t i, size_t j);

Description:

Access to the matrix element at the $i$ -th row and $j$ -th column (the same way as calling from a regular matrix object).

1.

Read only access from a constant view.
2.

Read only access from a non-constant view.
3.

Read and write access from a non-constant view.

4.7 Transposed View:

1 ViewConst ViewConst::operator*() const;

2 View View::operator*();

3 ViewConst View::operator*() const;

4 View Matrix::operator* ();

5 ViewConst Matrix::operator* () const;

Description:

Note that these functions are member functions of class ViewConst, View and Matrix respectively.

Let $\mathbf{A}$ be a matrix and $\mathbf{V}$ be a view of $\mathbf{A}$ . The function calls $*\mathbf{V}$ and $*\mathbf{A}$ both return a new view that represents the transpose matrix $\mathbf{A}^{T}$ .

1.

$\mathbf{V}$ is a constant view. The transpose of $\mathbf{V}$ is always a constant view.
2.

$\mathbf{V}$ is a non-constant view. The transpose of $\mathbf{V}$ is a non-constant view.
3.

$\mathbf{V}$ is a non-constant view. The transpose of $\mathbf{V}$ can be a constant view when $\mathbf{V}$ has constant modifier.
4.

$\mathbf{A}$ is a non-constant matrix. The transpose of $\mathbf{A}$ is a non-constant view.
5.

The transpose of $\mathbf{A}$ is a constant view when $\mathbf{A}$ has constant modifier.

Example 1:

1 // Create a view of 2 x 1 column vector.

2 const double _DATA[] = { 1, 2 };

3 Matrix::ViewConst V(2, 1, _DATA);

4 assert(V->is_col_vector());

6 // *V creates a transposed view of V.

7 Matrix::ViewConst U = *V;

8 assert(U->is_row_vector());

9 std::cout << U(0,0) << std::endl; // output 1

10 std::cout << U(0,1) << std::endl; // output 2

Example 2:

1 Matrix m(2,2);

3 // v is a non-constant view of the 1st row or m.

4 Matrix::View v = m.row(0);

6 // u is also non-constant, transpose v to a column.

7 Matrix::View u = *v;

9 // Apply constant modifier to v.

10 const Matrix::View& vr = v;

12 // Now the transpose returns a constant view.

13 Matrix::ViewConst w = *vr;

4.8 Permutated View:

1 ViewConst ViewConst::operator*(const Permutation&) const;

2 View View::operator*(const Permutation&);

3 ViewConst View::operator*(const Permutation&) const;

4 View Matrix::operator* (const Permutation& P);

5 ViewConst Matrix::operator* (const Permutation& P) const;

Description:

Note that these functions are member functions of class ViewConst, View and Matrix respectively.

Let $\mathbf{A}$ be a matrix and $\mathbf{V}$ be a view of $\mathbf{A}$ . Let $\mathbf{P}$ be a permutation. These multiplication operators return a new view that is equivalent to $\mathbf{A}\times\mathbf{P}$ . For the permutated view equivalent to $\mathbf{P}\times\mathbf{A}$ , please refer to the multiplication operators defined in class Matrix::Permutation.

REQUIRED: $\mathbf{P}$ , $\mathbf{A}$ and $\mathbf{V}$ have compatible size for multiplication.

1.

Return a permutated view that is still a constant view.
2.

Return a permutated view that is a non-constant view. Elements in this permutated view can be modified.
3.

Return a permutated view that could be a constant view when $\mathbf{V}$ has constant modifier.
4.

Return a non-constant view (elements can be modified).
5.

Same as above, except that the return is a constant view.

4.9 Arithmetic Operators:

1 void View::operator+= (const Matrix& M);

2 void View::operator+= (double v);

3 void View::operator-= (const Matrix& M);

4 void View::operator-= (double v);

5 void View::operator*= (const Matrix& M);

6 void View::operator*= (double v);

7 void View::operator*= (const Givens& G);

Description:

These functions are the member functions of class View, since they modify the content of the underlying matrix. Let $\mathbf{V}$ be the view representation of a matrix.

1.

Performs $\mathbf{V}=\mathbf{V}+\mathbf{M}$ . REQUIRED: The operand matrix $\mathbf{M}$ has same size as $\mathbf{V}$ .
2.

Add all element in $\mathbf{V}$ by the same number $v$ .
3.

Performs $\mathbf{V}=\mathbf{V}-\mathbf{M}$ . REQUIRED: The operand matrix $\mathbf{M}$ has same size as $\mathbf{V}$ .
4.

Subtract all element in $\mathbf{V}$ by the same number $v$ .
5.

Performs $\mathbf{V}=\mathbf{V}\times\mathbf{M}$ . REQUIRED: The operand matrix $\mathbf{M}$ has compatible size.
6.

Scale all element in $\mathbf{V}$ by the same number $v$ .
7.

Performs $\mathbf{V}=\mathbf{V}\times\mathbf{G}$ . REQUIRED: The operand Givens rotation matrix $\mathbf{G}$ has compatible size.

4.10 Sub Matrix Views:

The functions listed below are Matrix class member functions.

1 View row(size_t i);

2 ViewConst row(size_t i) const;

3 View col(size_t j);

4 ViewConst col(size_t j) const;

5 View sub(size_t i1, size_t j1, size_t i2, size_t j2);

6 ViewConst sub(size_t i1, size_t j1, size_t i2, size_t j2) const;

7 View self();

8 ViewConst self() const;

Description:

All these functions are member functions of class Matrix. Let $\mathbf{M}$ be the current matrix object. In the same order of the above code block,

1.

Return a non-constant view of the $i$ -th row of $\mathbf{M}$ .
2.

Same as above, except that the return is a constant view.
3.

Return a non-constant view of the $j$ -th column of $\mathbf{M}$ .
4.

Same as above, except that the return is a constant view.
5.

Return a non-constant view of the sub matrix of $\mathbf{M}$ . The sub matrix starts from row $i_{1}$ and column $j_{1}$ (inclusive). It stops at row $i_{2}$ and column $j_{2}$ (exclusive).
6.

Same as above, except that the return is a constant view.
7.

Return $\mathbf{M}$ itself as a non-constant view.
8.

Same as above, except that the return is a constant view.

4.11 Matrix Expressions

Matrix expressions are nested classes of class Matrix, all designed as derived types from the View base class. Because they are views, an expression object can appear anywhere a regular matrix is expected. In extra, they have their own memory to store data elements. Expression objects are created by the following operators.

1 Negation operator-(const Matrix& A);

2 Addition operator+(const Matrix& A, const Matrix& B);

3 Subtraction operator-(const Matrix& A, const Matrix& B);

4 Multiplication operator*(const Matrix& A, const Matrix& B);

5 Multiplication operator*(const Matrix& A, double f);

6 Multiplication operator*(double f, const Matrix& A);

Description:

Note that these operators are library scope operators. They can be exposed to the global scope with the using statement at a client application.

Let $\mathbf{A}$ and $\mathbf{B}$ be two matrix objects. Let $f$ be a scalar.

REQUIRED: $\mathbf{A}$ and $\mathbf{B}$ have compatible size relevent to the expression.

In the same order of the above code block,

1.

Negation: unary operator - on operand $\mathbf{A}$ creates a negation expression $-\mathbf{A}$ .
2.

Addition: binary operator + on operands $\mathbf{A}$ and $\mathbf{B}$ creates an addition expression $\mathbf{A}+\mathbf{B}$ .
3.

Subtraction: binary operator - on operands $\mathbf{A}$ and $\mathbf{B}$ creates an subtraction expression $\mathbf{A}-\mathbf{B}$ .
4.

Multiplication: binary operator * on operands $\mathbf{A}$ and $\mathbf{B}$ creates an multiplication expression $\mathbf{A}*\mathbf{B}$ .
5.

Multiplication: binary operator * on operands $\mathbf{A}$ and $f$ creates an multiplication expression $\mathbf{A}*f$ .
6.

Multiplication: binary operator * on operands $f$ and $\mathbf{A}$ creates an multiplication expression $f*\mathbf{A}$ .

5 Matrix Iterators

Matrix element iterators allow more efficient element access than the Matrix class member function operator()(size_t i, size_t j). In order to use iterators, we first need obtain an iterator support from a matrix or a view object. An iterator support is effectively a unique pointer of some iterator traits. If we are calling the iterator support funtion from a regular matrix object or from a view object other than permutated views, this traits pointer will be not null. We can then ask for iterators from the traits.

5.1 Iterator Supports:

Iterator supports are defined as the following inside the Matrix class:

1 typedef std::unique_ptr<iterator_traits> iterator_support;

2 typedef std::unique_ptr<const_iterator_traits>

3 const_iterator_support;

Iterator traits is a template nested struct of the Matrix class. The template type argument _Element is double for iterator_traits and const double for const_iterator_traits. The implementation details of iterator traits is not a focus here. From programming interface perspective, it is sufficient to know that:

1.

iterator_traits has row and column iterator types that allow both read and write access;

1 Matrix::iterator_traits::row_iterator;

2 Matrix::iterator_traits::col_iterator;
2.

const_iterator_traits has row and column iterator types with read access only.

1 Matrix::const_iterator_traits::row_iterator;

2 Matrix::const_iterator_traits::col_iterator;

Each of these two iterator traits types has the following member functions that return row or column traverse iterators.

1 row_iterator row_begin(size_t i) const;

2 row_iterator row_end(size_t i) const;

3 col_iterator col_begin(size_t j) const;

4 col_iterator col_end(size_t j) const;

Class Matrix has two member functions to return iterator supports.

1 Matrix::iterator_support Matrix::it_support();

2 Matrix::const_iterator_support Matrix::it_support() const;

1.

From a non-constant matrix object, we can obtain non-const iterator support, then non-const iterators to read or write elements.
2.

From a constant matrix, we can only obtain const iterator support, then const iterators for read only.

5.2 Iterator Operations:

Let it denote an iterator object obtained from an iterator traits. Let _Element be the template type argument that instantiates the relevant iterator traits. Let T denote the traits as abbreviation. Table 1 summarizes iterator operations. When iterator it reaches to a row end or column end, those non-comparison operators have undefined behavior.

Table 1: Iterator Operators

Operator	Syntax	Inside class definition	Remark
Deference	$*$ it	_Element& T::operator $*$ ();	Return a reference to the matrix
			element it points to.
Pre-increment	$++$ it	_Element& T::operator $++$ ();	Move it to the next element
			in the same row if it is a row
			iterator, or the next element in
			the same column if it is a column
			iterator.
			After the move, return a reference
			to the element that it now points
			to.
Pre-decrement	$--$ it	_Element& T::operator $--$ ();	Move it to the previous element
			in the same row if it is a row
			iterator, or the previous element in
			the same column if it is a column
			iterator.
			After the move, return a reference
			to the element that it now points
			to.
Post-increment	it $++$	_Element T::operator $++$ (int);	Same as pre-increment, except that
			the operation first saves the element
			value and then move. The operation
			returns the saved previous element
			value.
Post-decrement	it $--$	_Element T::operator $--$ (int);	Same as pre-decrement, except that
			the operation first saves the element
			value and then move. The operation
			returns the saved previous element
			value.
Jump-increment	it += n	void T::operator+=(size_t);	Move it to the element in the same
			row and n column greater if it is a
			row iterator,
			or the element in the same column
			and n row greater if it is a
			column iterator.
Jump-decrement	it -= n	void T::operator-=(size_t);	Move it to the element in the same
			row and n column smaller if it is a
			row iterator,
			or the element in the same column
			and n row smaller if it is a
			column iterator.
Comparision	it1 $==$ it2	bool T::operator $==$ (	Compare if two iterators pointing
Equal		const self_type&) const;	to the same element.
Comparision	it1 $!=$ it2	bool T::operator $!=$ (	Compare if two iterators pointing
Not Equal		const self_type&) const;	to different elements.
Comparision	it1 $<$ it2	bool T::operator $<$ (	Compare if iterator it1 pointing
Less		const self_type&) const;	to an element that has smaller
			memory address than the element
			pointed by it2.

6 Matrix Factorization

All functions listed in this section are Matrix class member functions. Let $\mathbf{A}$ be the matrix object to factorize.

6.1 Eigen

1 void eigen(Matrix& eigen_vect, Matrix& eigen_value,

2 Matrix::Form f = GENERAL,

3 Matrix* eigen_img = 0) const;

Description:

This function computes the eigen vectors and eigen values of $\mathbf{A}$ .

REQUIRED: $\mathbf{A}$ is a $n\times n$ square matrix.

Parameters:

1.

eigen_vect will be resized to $n\times n$ matrix that holds the right eigen vectors in columns.
2.

eigen_value will be resized to $n\times 1$ to hold the eigen values on return.
3.

f indicates if $\mathbf{A}$ is symmetric.
4.

eigen_img is relevant to non-symmetric matrix. If supplied, it will be resized to hold the imaginary part of eigen values.

6.2 SVD

1 void svd (Matrix& U, Matrix& S, Matrix& V) const;

Description:

Given a $m\times n$ matrix $\mathbf{A}$ , this function computes a compact form of singular value decomposition $\mathbf{A}=\mathbf{U}\mathbf{S}\mathbf{V}^{T}$ . Let $r=\min(m,n)$ . $\mathbf{U}$ has size $m\times r$ . $\mathbf{S}$ has size $r\times r$ . $\mathbf{V}$ has size $n\times r$ .

Parameters:

1.

U is resized to $m\times r$ to hold matrix $\mathbf{U}$ on return.
2.

S is resized to $r\times 1$ to hold $\mathbf{S}$ on return. All elements in $\mathbf{S}$ are non-negative and are arranged in descending order.
3.

V is resized to $n\times r$ to hold the matrix $\mathbf{V}$ , not $\mathbf{V}^{T}$ on return.

6.3 Cholesky

1 bool cholesky(Matrix& L) const;

Description:

For a positive definite matrix $\mathbf{A}$ , this function computes the Cholesky decomposition $\mathbf{A}=\mathbf{L}\mathbf{L}^{T}$ . If $\mathbf{A}$ is not positive definite, the function returns false.

Parameters:

L will be resized and hold the lower triangular matrix $\mathbf{L}$ .

6.4 Reduced Cholesky

1 bool cholesky_reduced(const Matrix& W, Matrix& L) const;

Description:

Given a $n\times n$ positive definite matrix $\mathbf{A}$ and a $n\times m$ matrix $\mathbf{W}$ , $m<n$ , this function computes the Cholesky factor,

\mathbf{W}^{T}\mathbf{A}\mathbf{W}=\mathbf{L}\mathbf{L}^{T}

It returns false if $\mathbf{A}$ is actually not positive definite. Note that when $\mathbf{A}$ is indefinite, it is possible that $\mathbf{W}^{T}\mathbf{A}\mathbf{W}$ could still be positive definite. This reduced Cholesky function only works when $\mathbf{A}$ is positive definite.

Parameters:

L will be resized to $m\times m$ and hold the lower triangular matrix $\mathbf{L}$ .

6.5 Modified Cholesky

1 void cholesky_modified(Matrix& L, Matrix* _E = 0) const;

Description:

Given a $n\times n$ matrix $\mathbf{A}$ that could be indefinite, this function finds a miniminzed norm diagonal remedy matrix $\mathbf{E}$ such that

\mathbf{A}+\mathbf{E}=\mathbf{L}\mathbf{L}^{T}

Parameters:

1.

L will be resized to $n\times n$ and hold the lower triangular matrix $\mathbf{L}$ .
2.

_E is an optional argument. When provided, the deferenced matrix object (*_E) will be resized to $n\times 1$ to holds the diagonal values of $\mathbf{E}$ .

6.6 Partial Cholesky

1 size_t cholesky_partial(Matrix& L, Matrix& B,

2 Matrix::Permutation& P, double nu = 0.5) const;

Description:

Given a $n\times n$ matrix $\mathbf{A}$ that could be indefinite, this function finds a permutation $\mathbf{P}$ and computes the partial Cholesky factoring,

\mathbf{P}^{T}\mathbf{A}\mathbf{P}=\mathbf{L}\mathbf{B}\mathbf{L}^{T}

, where $\mathbf{L}$ and $\mathbf{B}$ have the following block structure.

	$\displaystyle\mathbf{L}$	$\displaystyle=$	$\displaystyle\left(\begin{array}[]{cc}\hidden@noalign{}\hfil\mathbf{L}_{11}&% \cr\hidden@noalign{}\hfil\mathbf{L}_{21}&\mathbf{I}\end{array}\right)$
	$\displaystyle\mathbf{B}$	$\displaystyle=$	$\displaystyle\left(\begin{array}[]{cc}\hidden@noalign{}\hfil\mathbf{D}&\cr&% \mathbf{K}\end{array}\right)$

$\mathbf{L}_{11}$ is a lower triangular unitary matrix. $\mathbf{D}$ is a diagonal matrix with positive elements. The function returns the size of $\mathbf{L}_{11}$ .

Parameters:

1.

L will be resized to $n\times n$ to hold matrix $\mathbf{L}$ .
2.

B will be resized to $n\times n$ to holds matrix $\mathbf{B}$ .
3.

P will be resized to hold the permutation $\mathbf{P}$ .
4.

nu is a number between 0 and 1 to control the termination condition. Please refer to our white paper on this subject for details.

6.7 LU

1 void lu(Matrix& L, Matrix& U, Permutation& p) const;

Description:

Given a general $m\times n$ matrix $\mathbf{A}$ , this function computes the LU decomposition $\mathbf{A}=\mathbf{P}\mathbf{L}\mathbf{U}$ , where $\mathbf{P}$ is the permutation introduced by pivoting, $\mathbf{L}$ is $m\times r$ , $\mathbf{U}$ is $r\times n$ and $r=\min(m,n)$ . $\mathbf{L}$ is lower triangular if $m\leq n$ or lower trapezoidal if $m>n$ . $\mathbf{U}$ is upper triangular if $m\geq n$ or upper trapezoidal if $m<n$ .

Parameters:

1.

L will be resized to $m\times r$ to hold matrix $\mathbf{L}$ .
2.

U will be resized to $r\times n$ to hold matrix $\mathbf{U}$ .
3.

P will hold the permutation on return.

6.8 QR

1 void qr(Matrix& R, Matrix* _Q = 0) const;

Description:

Given a general $m\times n$ matrix $\mathbf{A}$ , this function computes a compact form of QR decomposition $\mathbf{A}=\mathbf{Q}\mathbf{R}$ , where $\mathbf{Q}$ is $m\times r$ , $\mathbf{R}$ is $r\times n$ and $r=\min(m,n)$ . $\mathbf{Q}$ has orthonomal vectors in columns. $\mathbf{R}$ is upper triangular if $m\geq n$ or upper trapezoidal if $m<n$ .

Parameters:

1.

R will be resized to $r\times n$ to hold matrix $\mathbf{R}$ .
2.

_Q is optional. If provided, the deferenced matrix object (*_Q) will be resized to $m\times r$ to hold matrix $\mathbf{Q}$ ,

6.9 Full QR

1 void full_qr (Matrix& R, Matrix* _Q = 0) const;

Description:

Given a general $m\times n$ matrix $\mathbf{A}$ , this function computes a full QR decomposition $\mathbf{A}=\mathbf{Q}\mathbf{R}$ , where $\mathbf{Q}$ is $m\times m$ and $\mathbf{R}$ is $m\times n$ . $\mathbf{Q}$ has orthonomal vectors in columns. If $m\geq n$ , $\mathbf{R}$ has $n\times n$ upper triangular at the top and the rest rows be 0. If $m<n$ , $\mathbf{R}$ is upper trapezoidal.

Parameters:

1.

R will be resized to $m\times n$ to hold matrix $\mathbf{R}$ .
2.

_Q is optional. If provided, the deferenced matrix object (*_Q) will be resized to $m\times m$ to hold matrix $\mathbf{Q}$ ,

7 Linear System

LinearSystem class includes a group of static functions and two nested classes LinearSystem::Tridiagonal and LinearSystem::LU to solve a linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$ . The C++ header file is linear_system.h.

7.1 General Solver

1 static void LinearSystem::solve(

2 const Matrix& A,

3 const Matrix& b,

4 Matrix& x,

5 Matrix::Form formA = Matrix::GENERAL);

Description:

This function is the most general solver and allows $\mathbf{A}$ to be a $m\times n,m\neq n$ matrix. Depending on the nature of $\mathbf{A}$ , $\mathbf{x}$ could be an exact solution or a least square solution. Parameter formA describes the form of $\mathbf{A}$ , which is used to choose appropriate algorithms. Please refer to our white paper for technical details. On function return, parameter x is resized to $n\times 1$ to hold the solution $\mathbf{x}$ .

7.2 Special Cases

7.2.1 Forward Substitution

1 static void LinearSystem::fwd_subst(

2 const Matrix& L,

3 const Matrix& b,

4 Matrix& x);

5 static void LinearSystem::fwd_subst_e_j(

6 const Matrix& L,

7 size_t j,

8 Matrix& x);

When $\mathbf{A}$ is a $n\times n$ lower triangular matrix $\mathbf{L}$ , these two functions solves $\mathbf{L}\mathbf{x}=\mathbf{b}$ and $\mathbf{L}\mathbf{x}=\mathbf{e}_{j}$ respectively. $\mathbf{e}_{j}$ is a unit vector with the $j$ -th element equal to 1 and all other elements are 0. Parameter x is resized to hold solution $\mathbf{x}$ on return.

7.2.2 Backward Substitution

1 static void LinearSystem::bwd_subst (

2 const Matrix& U,

3 Matrix& y);

4 static void LinearSystem::bwd_subst_e_j(

5 const Matrix& U,

6 size_t j,

7 Matrix& y);

When $\mathbf{A}$ is a $n\times n$ upper triangular matrix $\mathbf{U}$ . The first function solves $\mathbf{U}\mathbf{y}=\mathbf{b}$ , $\mathbf{b}$ is passed as parameter y on call. On return, parameter y holds the solution vector $\mathbf{y}$ . REQUIRED: Parameter y is a $n$ -element vector.

The second function solves $\mathbf{U}\mathbf{y}=\mathbf{e}_{j}$ . $\mathbf{e}_{j}$ is a unit vector with the $j$ -th element equal to 1 and all other elements are 0. Parameter y is resized to hold the solution $\mathbf{y}$ on return.

7.2.3 Conjugate Gradient Update

1 static void LinearSystem::cg_step_update(

2 const Matrix& A,

3 Matrix& x_k,

4 Matrix& p_k,

5 Matrix& r_k);

Description:

When $\mathbf{A}$ is a large $n\times n$ positive definite matrix, the conjugate gradient (CG) method can solve the linear system by generating a sequence of improved guess $\mathbf{x}_{k}$ , $k=0,\cdots,n$ . At each iteration, it also outputs the by-product of a residual vector $\mathbf{r}_{k}$ and a search direction $\mathbf{p}_{k}$ to be taken at the $(k+1)$ -th iteration.

Initially, users need choose an initial guess $\mathbf{x}_{0}$ . The residual vector should be computed as $\mathbf{r}_{0}=\mathbf{A}\mathbf{x}_{0}-\mathbf{b}$ . The initial search direction should be taken as $\mathbf{p}_{0}=-\mathbf{r}_{0}$ . Subsequently, $\left(\mathbf{x}_{k},\mathbf{p}_{k},\mathbf{r}_{k}\right)$ are updated by the function. User can stop the update when $\mathbf{r}_{k}$ is considered as small enough. In the worst case, the CG method finds the exact solution after $n$ iterations.

REQUIRED: Parameter x_k, p_k and r_k for $\mathbf{x}_{k}$ , $\mathbf{p}_{k}$ and $\mathbf{r}_{k}$ are $n$ -element column vectors.

7.2.4 Range Space Solver

1 static void LinearSystem::range_space_solve(

2 const Matrix& Q,

3 const Matrix& R,

4 const Matrix& b,

5 Matrix& x,

6 bool transpose = false);

Description:

Given a $m\times n$ matrix $\mathbf{A}$ , $m\leq n$ , if the QR decomposition of $\mathbf{A}^{T}=\mathbf{Q}\mathbf{R}$ is known, we can use this range space solver to solve $\mathbf{A}\mathbf{x}=\mathbf{b}$ when parameter transpose is false, or $\mathbf{A}^{T}\mathbf{x}=\mathbf{b}$ when parameter transpose is passed true.

When $m<n$ , $\mathbf{A}\mathbf{x}=\mathbf{b}$ is underdetermined. This solution vector $\mathbf{x}$ being found is in the span of the range space of $\mathbf{A}^{T}$ . On the other hand, $\mathbf{A}^{T}\mathbf{x}=\mathbf{b}$ is overdetermined, the $\mathbf{x}$ computed is a least square solution.

Parameter x is resized to hold the solution $\mathbf{x}$ on return.

7.2.5 Tridiagonal

When $\mathbf{A}$ is a $n\times n$ tridiagonal matrix, LinearSystem::Tridiagonal solver can be used to compute the solution $\mathbf{x}$ in $2n$ floating point operations.

LinearSystem::Tridiagonal a nested class of LinearSystem class. It has the following member functions.

Constructor and Member Functions:

1 Tridiagonal(size_t n = 0);

2 size_t size() const { return m_D.size(); }

3 void resize(size_t n = 0);

4 double sub_diag(size_t i) const;

5 double& sub_diag(size_t i);

6 double diag(size_t i) const;

7 double& diag(size_t i);

8 double supper_diag(size_t i) const;

9 double& supper_diag(size_t i);

10 void multiply(const Matrix& X, Matrix& Y) const;

11 void solve(const Matrix& B, Matrix& X) const;

Description:

1.

Constructor, create a solver for $n\times n$ tridiagonal system.
2.

Return the size $n$ of the system.
3.

Resize the solver to a new size.
4.

Return the value of the sub-diagonal element at the $i$ -th row.
5.

Return a reference to the sub-diagonal element at the $i$ -th row.
6.

Return the value of the diagonal element at the $i$ -th row.
7.

Return a referece to the diagonal element at the $i$ -th row.
8.

Return the value of the super-diagonal element at the $i$ -th row.
9.

Return a reference to the super-diagonal element at the $i$ -th row.
10.

Compute the multiplication $\mathbf{Y}=\mathbf{A}\mathbf{X}$ . Parameter Y is resized to hold the multiplication result $\mathbf{Y}$ on return. REQUIRED: : Parameter X for matrix $\mathbf{X}$ has $n$ rows.
11.

Solve the linear system $\mathbf{A}\mathbf{X}=\mathbf{B}$ , where $\mathbf{B}$ could be a matrix. Parameter X is resized to hold the solution matrix $\mathbf{X}$ on return. Column $\mathbf{X}_{i}$ is the solution of right hand side column $\mathbf{B}_{i}$ . REQUIRED: Parameter B for matrix $\mathbf{B}$ has $n$ rows.

7.2.6 LU Update

LinearSystem::LU class is a nested class of the LinearSystem class. It is useful in the situation when

1.

$\mathbf{A}_{0}$ is a $n\times n$ matrix and;
2.

Subsequently $\mathbf{A}_{i+1}$ is obtained by pivoting a column of $\mathbf{A}_{i}$ with another vector $\mathbf{p}$ and
3.

We want to repeatedly solve $\mathbf{A}_{i+1}\mathbf{x}=\mathbf{b}$ .

We initially compute the LU decomposition of $\mathbf{A}_{0}=\mathbf{L}_{0}\mathbf{U}_{0}$ . The LU decomposition of $\mathbf{A}_{i+1}$ can economically updated from $\mathbf{L}_{i}$ and $\mathbf{U}_{i}$ , hence obtaining the solution of updated linear system.

Constructor and Member Functions:

1 LU();

2 LU(const LU& other);

3 LU(size_t n, size_t t);

4 ~LU();

5 void resize(size_t n, size_t t);

6 void reset(const Matrix& A);

7 void reset_diagonal(const Matrix& D);

8 void pivot (const Matrix& p, size_t j);

9 const Matrix& current_basis() const;

10 void operator()(const Matrix& b, Matrix& x, bool trans=false);

Description:

1.

Default constructor, create an empty object.
2.

Copy constructor, deep copy all LU update solver internal details.
3.

Initialize the LU update solver for a $n\times n$ system with at most $t$ number of pivoting operations to perform.
4.

Destructor, free all internal resouces.
5.

Resize the solver to a new $n\times n$ linear system and $t$ number of pivoting.
6.

Parameter A is takens as $\mathbf{A}_{0}$ , compute the initial LU decomposition. REQUIRED: A has correct matrix size.
7.

Parameter D is takens as $\mathbf{A}_{0}$ , which is a diagonal matrix. Initialize $\mathbf{L}_{0}$ to an identity matrix, and $\mathbf{U}_{0}=\mathbf{A}_{0}$ . REQUIRED: D is a $n$ -element column vector.
8.

Perform the pivoting operation to swap the $j$ -th column in the current $\mathbf{A}_{i}$ with a vector $\mathbf{p}$ to obtain $\mathbf{A}_{i+1}$ . $\mathbf{L}_{i+1}$ and $\mathbf{U}_{i+1}$ are updated. REQUIRED: p is $n$ -element column vector.
9.

Return a constant reference to the current basis matrix $\mathbf{A}_{i}$ .
10.

If parameter trans is passed false, solve the linear system $\mathbf{A}_{i}\mathbf{x}=\mathbf{b}$ ; otherwise, solve $\mathbf{A}_{i}^{T}\mathbf{x}=\mathbf{b}$ . Parameter x is resized and holds the solution on return. REQUIRED: b is a $n$ -element vector.

Numerical Linear Algebra Mathwrist C++ API Series

Abstract

1 Matrix Basics

1.1 Enumeration of Matrix Forms

1.2 Object Creation

1.3 Static Initialization Methods

1.4 Element Access

1.5 Dimension and Size Query

1.6 Comparison and Content Query

1.7 Assignment Operators

1.8 Manipulation Functions

1.9 Norm, Min and Max Functions

2 Matrix Arithmetics

2.1 Arithmetic Operators

2.2 Element-wise Multiplication and Division

2.3 Matrix Multiplication and Inverse

3 Transformation Matrices

3.1 Permutation

3.1.1 Constructors

3.1.2 Two Useful Factory Methods:

3.1.3 Multiplications

3.1.4 Query Functions:

3.1.5 Manipulation Functions:

3.2 Givens Rotation

3.2.1 Constructors:

3.2.2 Rotation Operation:

3.2.3 Query Functions:

3.2.4 Manipulation Functions:

4 Matrix Views

4.1 Constructors:

4.2 Conversion Operators:

4.3 Pointer Operators:

4.4 Assignment Operators:

4.5 Copy Operation:

4.6 Element Access:

4.7 Transposed View:

4.8 Permutated View:

4.9 Arithmetic Operators:

4.10 Sub Matrix Views:

4.11 Matrix Expressions

5 Matrix Iterators

5.1 Iterator Supports:

5.2 Iterator Operations:

6 Matrix Factorization

6.1 Eigen

6.2 SVD

6.3 Cholesky

6.4 Reduced Cholesky

6.5 Modified Cholesky

6.6 Partial Cholesky

6.7 LU

6.8 QR

6.9 Full QR

7 Linear System

7.1 General Solver

7.2 Special Cases

7.2.1 Forward Substitution

7.2.2 Backward Substitution

7.2.3 Conjugate Gradient Update

7.2.4 Range Space Solver

7.2.5 Tridiagonal

7.2.6 LU Update

Numerical Linear Algebra
Mathwrist C++ API Series