Matrix Classes in C++ - The Source File

Matrix Classes in C++ - The Source File

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.