This is the second part in a two part article series on how to create a robust, extensible and reusable *matrix class* for the necessary numerical linear algebra work needed for algorithms in quantitative finance. In order to familiarise yourself with the *specification* and *declaration* of the class, please read the first article on the matrix header file before continuing.

Our task with the source file is to implement all of the methods outlined in the header file. In particular we need to implement methods for the following:

- Constructors (parameter and copy), destructor and assignment operator
- Matrix mathematical methods: Addition, subtraction, multiplication and the transpose
- Matrix/scalar element-wise mathematical methods: Addition, substraction, multiplication and division
- Matrix/vector multiplication methods
- Element-wise access (const and non-const)

We will begin with the construction, assignment and destruction of the class.

## Allocation and Deallocation

The first method to implement is the constructor, with paramaters. The constructor takes three arguments - the number of rows, the number of columns and an initial type value to populate the matrix with. Since the "vector of vectors" constructor has already been called at this stage, we need to call its `resize`

method in order to have enough elements to act as the row containers. Once the matrix `mat`

has been resized, we need to resize each individual vector within the rows to the length representing the number of columns. The `resize`

method can take an optional argument, which will initialise all elements to that particular value. Finally we adjust the private `rows`

and `cols`

unsigned integers to store the new row and column counts:

// Parameter Constructor template<typename T> QSMatrix<T>::QSMatrix(unsigned _rows, unsigned _cols, const T& _initial) { mat.resize(_rows); for (unsigned i=0; i<mat.size(); i++) { mat[i].resize(_cols, _initial); } rows = _rows; cols = _cols; }

The copy constructor has a straightforward implementation. Since we have not used any *dynamic memory allocation*, we simply need to copy each private member from the corresponding copy matrix `rhs`

:

// Copy Constructor template<typename T> QSMatrix<T>::QSMatrix(const QSMatrix<T>& rhs) { mat = rhs.mat; rows = rhs.get_rows(); cols = rhs.get_cols(); }

The destructor is even simpler. Since there is no dynamic memory allocation, we don't need to do anything. We can let the compiler handle the destruction of the individual type members (`mat`

, `rows`

and `cols`

):

// (Virtual) Destructor template<typename T> QSMatrix<T>::~QSMatrix() {}

The assignment operator is somewhat more complicated than the other construction/destruction methods. The first two lines of the method implementation check that the addresses of the two matrices aren't identical (i.e. we're not trying to assign a matrix to itself). If this is the case, then just return the dereferenced pointer to the current object (`*this`

). This is purely for performance reasons. Why go through the process of copying exactly the same data into itself if it is already identical?

However, if the matrix addresses differ, then we resize the old matrix to the be the same size as the `rhs`

matrix. Once that is complete we then populate the values element-wise and finally adjust the members holding the number of rows and columns. We then return the dereferenced pointer to `this`

. This is a common pattern for assignment operators and is considered good practice:

// Assignment Operator template<typename T> QSMatrix<T>& QSMatrix<T>::operator=(const QSMatrix<T>& rhs) { if (&rhs == this) return *this; unsigned new_rows = rhs.get_rows(); unsigned new_cols = rhs.get_cols(); mat.resize(new_rows); for (unsigned i=0; i<mat.size(); i++) { mat[i].resize(new_cols); } for (unsigned i=0; i<new_rows; i++) { for (unsigned j=0; j<new_cols; j++) { mat[i][j] = rhs(i, j); } } rows = new_rows; cols = new_cols; return *this; }

## Mathematical Operators Implementation

The next part of the implementation concerns the methods overloading the binary operators that allow matrix algebra such as addition, subtraction and multiplication. There are two types of operators to be overloaded here. The first is operation *without assignment*. The second is operation *with assignment*. The first type of operator method creates a new matrix to store the result of an operation (such as addition), while the second type applies the result of the operation into the left-hand argument. For instance the first type will produce a new matrix $C$, from the equation $C=A+B$. The second type will overwrite $A$ with the result of $A+B$.

The first operator to implementation is that for addition without assignment. A new matrix `result`

is created with initial value $0$. Then each element is iterated through to be the pairwise sum of the `this`

matrix and the new right hand side matrix `rhs`

. Notice that we use the pointer dereferencing syntax with `this`

when accessing the element values: `this->mat[i][j]`

. This is identical to writing `(*this).mat[i][j]`

. We must dereference the pointer before we can access the underlying object. Finally, we return the result:

*Note that this can be a particularly expensive operation. We are creating a new matrix for every call of this method. However, modern compilers are smart enough to make sure that this operation is not as performance heavy as it used to be, so for our current needs we are justified in creating the matrix here. Note again that if we were to return a matrix by reference and then create the matrix within the class via the* `new`

*operator, we would have an error as the matrix object would go out of scope as soon as the method returned.*

// Addition of two matrices template<typename T> QSMatrix<T> QSMatrix<T>::operator+(const QSMatrix<T>& rhs) { QSMatrix result(rows, cols, 0.0); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { result(i,j) = this->mat[i][j] + rhs(i,j); } } return result; }

The operation with assignment method for addition is carried out slightly differently. It DOES return a reference to an object, but this is fine since the object reference it returns is to `this`

, which exists outside of the scope of the method. The method itself makes use of the `operator+=`

that is bound to the *type* object. Thus when we carry out the line `this->mat[i][j] += rhs(i,j);`

we are making use of the types own operator overload. Finally, we return a dereferenced pointer to `this`

giving us back the modified matrix:

// Cumulative addition of this matrix and another template<typename T> QSMatrix<T>& QSMatrix<T>::operator+=(const QSMatrix<T>& rhs) { unsigned rows = rhs.get_rows(); unsigned cols = rhs.get_cols(); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { this->mat[i][j] += rhs(i,j); } } return *this; }

The two matrix subtraction operators `operator-`

and `operator-=`

are almost identical to the addition variants, so I won't explain them here. If you wish to see their implementation, have a look at the full listing below.

I will discuss the matrix multiplication methods though as their syntax is sufficiently different to warrant explanation. The first operator is that without assignment, `operator*`

. We can use this to carry out an equation of the form $C = A \times B$. The first part of the method creates a new `result`

matrix that has the same size as the right hand side matrix, `rhs`

. Then we perform the triple loop associated with matrix multiplication. We iterate over each element in the result matrix and assign it the value of `this->mat[i][k] * rhs(k,j)`

, i.e. the value of $A_{ik} \times B_{kj}$, for $k \in \{0,…,M-1\}$:

// Left multiplication of this matrix and another template<typename T> QSMatrix<T> QSMatrix<T>::operator*(const QSMatrix<T>& rhs) { unsigned rows = rhs.get_rows(); unsigned cols = rhs.get_cols(); QSMatrix result(rows, cols, 0.0); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { for (unsigned k=0; k<rows; k++) { result(i,j) += this->mat[i][k] * rhs(k,j); } } } return result; }

The implementation of the `operator*=`

is far simpler, but only because we are building on what already exists. The first line creates a new matrix called `result`

which stores the result of multiplying the dereferenced pointer to `this`

and the right hand side matrix, `rhs`

. The second line then sets `this`

to be equal to the result above. This is necessary as if it was carried out in one step, data would be overwritten before it could be used, creating an incorrect result. Finally the referenced pointer to `this`

is returned. Most of the work is carried out by the `operator*`

which is defined above. The listing is as follows:

// Cumulative left multiplication of this matrix and another templateQSMatrix & QSMatrix ::operator*=(const QSMatrix & rhs) { QSMatrix result = (*this) * rhs; (*this) = result; return *this; }

We also wish to apply scalar element-wise operations to the matrix, in particular element-wise scalar addition, subtraction, multiplication and division. Since they are all very similar, I will only provide explanation for the addition operator. The first point of note is that the parameter is now a `const T&M`

, i.e. a reference to a const type. This is the scalar value that will be added to all matrix elements. We then create a new `result`

matrix as before, of identical size to `this`

. Then we iterate over the elements of the result matrix and set their values equal to the sum of the individual elements of `this`

and our type value, `rhs`

. Finally, we return the result matrix:

// Matrix/scalar addition template<typename T> QSMatrix<T> QSMatrix<T>::operator+(const T& rhs) { QSMatrix result(rows, cols, 0.0); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { result(i,j) = this->mat[i][j] + rhs; } } return result; }

We also wish to allow (right) matrix vector multiplication. It is not too different from the implementation of matrix-matrix multiplication. In this instance we are returning a `std::vector`

and also providing a separate vector as a parameter. Upon invocation of the method we create a new `result`

vector that has the same size as the right hand side, `rhs`

. Then we perform a double loop over the elements of the `this`

matrix and assign the result to an element of the `result`

vector. Finally, we return the `result`

vector:

// Multiply a matrix with a vector template<typename T> std::vector<T> QSMatrix<T>::operator*(const std::vector<T>& rhs) { std::vector<T> result(rhs.size(), 0.0); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { result[i] = this->mat[i][j] * rhs[j]; } } return result; }

I've added a final matrix method, which is useful for certain numerical linear algebra techniques. Essentially it returns a vector of the diagonal elements of the matrix. Firstly we create the `result`

vector, then assign it the values of the diagonal elements and finally we return the `result`

vector:

// Obtain a vector of the diagonal elements template<typename T> std::vector<T> QSMatrix<T>::diag_vec() { std::vector<T> result(rows, 0.0); for (unsigned i=0; i<rows; i++) { result[i] = this->mat[i][i]; } return result; }

The final set of methods to implement are for accessing the individual elements as well as getting the number of rows and columns from the matrix. They're all quite simple in their implementation. They dereference `this`

and then obtain either an individual element or some private member data:

// Access the individual elements templateT& QSMatrix ::operator()(const unsigned& row, const unsigned& col) { return this->mat[row][col]; } // Access the individual elements (const) template const T& QSMatrix ::operator()(const unsigned& row, const unsigned& col) const { return this->mat[row][col]; } // Get the number of rows of the matrix template unsigned QSMatrix ::get_rows() const { return this->rows; } // Get the number of columns of the matrix template unsigned QSMatrix ::get_cols() const { return this->cols; }

## Full Source Implementation

Now that we have described all the methods in full, here is the full source listing for the `QSMatrix`

class:

#ifndef __QS_MATRIX_CPP #define __QS_MATRIX_CPP #include "matrix.h" // Parameter Constructor template<typename T> QSMatrix<T>::QSMatrix(unsigned _rows, unsigned _cols, const T& _initial) { mat.resize(_rows); for (unsigned i=0; i<mat.size(); i++) { mat[i].resize(_cols, _initial); } rows = _rows; cols = _cols; } // Copy Constructor template<typename T> QSMatrix<T>::QSMatrix(const QSMatrix<T>& rhs) { mat = rhs.mat; rows = rhs.get_rows(); cols = rhs.get_cols(); } // (Virtual) Destructor template<typename T> QSMatrix<T>::~QSMatrix() {} // Assignment Operator template<typename T> QSMatrix<T>& QSMatrix<T>::operator=(const QSMatrix<T>& rhs) { if (&rhs == this) return *this; unsigned new_rows = rhs.get_rows(); unsigned new_cols = rhs.get_cols(); mat.resize(new_rows); for (unsigned i=0; i<mat.size(); i++) { mat[i].resize(new_cols); } for (unsigned i=0; i<new_rows; i++) { for (unsigned j=0; j<new_cols; j++) { mat[i][j] = rhs(i, j); } } rows = new_rows; cols = new_cols; return *this; } // Addition of two matrices template<typename T> QSMatrix<T> QSMatrix<T>::operator+(const QSMatrix<T>& rhs) { QSMatrix result(rows, cols, 0.0); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { result(i,j) = this->mat[i][j] + rhs(i,j); } } return result; } // Cumulative addition of this matrix and another template<typename T> QSMatrix<T>& QSMatrix<T>::operator+=(const QSMatrix<T>& rhs) { unsigned rows = rhs.get_rows(); unsigned cols = rhs.get_cols(); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { this->mat[i][j] += rhs(i,j); } } return *this; } // Subtraction of this matrix and another template<typename T> QSMatrix<T> QSMatrix<T>::operator-(const QSMatrix<T>& rhs) { unsigned rows = rhs.get_rows(); unsigned cols = rhs.get_cols(); QSMatrix result(rows, cols, 0.0); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { result(i,j) = this->mat[i][j] - rhs(i,j); } } return result; } // Cumulative subtraction of this matrix and another template<typename T> QSMatrix<T>& QSMatrix<T>::operator-=(const QSMatrix<T>& rhs) { unsigned rows = rhs.get_rows(); unsigned cols = rhs.get_cols(); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { this->mat[i][j] -= rhs(i,j); } } return *this; } // Left multiplication of this matrix and another template<typename T> QSMatrix<T> QSMatrix<T>::operator*(const QSMatrix<T>& rhs) { unsigned rows = rhs.get_rows(); unsigned cols = rhs.get_cols(); QSMatrix result(rows, cols, 0.0); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { for (unsigned k=0; k<rows; k++) { result(i,j) += this->mat[i][k] * rhs(k,j); } } } return result; } // Cumulative left multiplication of this matrix and another template<typename T> QSMatrix<T>& QSMatrix<T>::operator*=(const QSMatrix<T>& rhs) { QSMatrix result = (*this) * rhs; (*this) = result; return *this; } // Calculate a transpose of this matrix template<typename T> QSMatrix<T> QSMatrix<T>::transpose() { QSMatrix result(rows, cols, 0.0); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { result(i,j) = this->mat[j][i]; } } return result; } // Matrix/scalar addition template<typename T> QSMatrix<T> QSMatrix<T>::operator+(const T& rhs) { QSMatrix result(rows, cols, 0.0); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { result(i,j) = this->mat[i][j] + rhs; } } return result; } // Matrix/scalar subtraction template<typename T> QSMatrix<T> QSMatrix<T>::operator-(const T& rhs) { QSMatrix result(rows, cols, 0.0); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { result(i,j) = this->mat[i][j] - rhs; } } return result; } // Matrix/scalar multiplication template<typename T> QSMatrix<T> QSMatrix<T>::operator*(const T& rhs) { QSMatrix result(rows, cols, 0.0); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { result(i,j) = this->mat[i][j] * rhs; } } return result; } // Matrix/scalar division template<typename T> QSMatrix<T> QSMatrix<T>::operator/(const T& rhs) { QSMatrix result(rows, cols, 0.0); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { result(i,j) = this->mat[i][j] / rhs; } } return result; } // Multiply a matrix with a vector template<typename T> std::vector<T> QSMatrix<T>::operator*(const std::vector<T>& rhs) { std::vector<T> result(rhs.size(), 0.0); for (unsigned i=0; i<rows; i++) { for (unsigned j=0; j<cols; j++) { result[i] = this->mat[i][j] * rhs[j]; } } return result; } // Obtain a vector of the diagonal elements template<typename T> std::vector<T> QSMatrix<T>::diag_vec() { std::vector<T> result(rows, 0.0); for (unsigned i=0; i<rows; i++) { result[i] = this->mat[i][i]; } return result; } // Access the individual elements template<typename T> T& QSMatrix<T>::operator()(const unsigned& row, const unsigned& col) { return this->mat[row][col]; } // Access the individual elements (const) template<typename T> const T& QSMatrix<T>::operator()(const unsigned& row, const unsigned& col) const { return this->mat[row][col]; } // Get the number of rows of the matrix template<typename T> unsigned QSMatrix<T>::get_rows() const { return this->rows; } // Get the number of columns of the matrix template<typename T> unsigned QSMatrix<T>::get_cols() const { return this->cols; } #endif

## Using the Matrix Class

We have the full listings for both the matrix header and source, so we can test the methods out with some examples. Here is the `main`

listing showing the matrix addition operator:

#include "matrix.h" #include <iostream> int main(int argc, char **argv) { QSMatrix<double> mat1(10, 10, 1.0); QSMatrix<double> mat2(10, 10, 2.0); QSMatrix<double> mat3 = mat1 + mat2; for (int i=0; i<mat3.get_rows(); i++) { for (int j=0; j<mat3.get_cols(); j++) { std::cout << mat3(i,j) << ", "; } std::cout << std::endl; } return 0; }

Here is the output of the code. We can see that the elements are all valued $3.0$, which is simply the element-wise addition of `mat1`

and `mat2`

:

3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,

We will make extensive use of this matrix class in our further numerical linear algebra routines, in particular with statistical analysis and finite difference methods.