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
template
QSMatrix& 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
template
T& 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.
