Build up you own Matrix Computation Lib in C++

Last week professor Namini posted a new C++ programming project task: to implement the Mean-Variance and Black-Litterman methods to optimize an equity portfolio. Of couse, this project would involve quite amount of matrix operation, which is a little bit complex to achieve in C++. Surely, there are some existing libraries could help us to do that (like Quantlib).

But an alternative way is to to develop our own matrix computation library, like what what I decided on a....weird day, yes it is weird.

And here we go~

The first thing when designing a new things is to figure out what functions or features we need our product to have.

#1. here are what we need!

(1) Generate a simple (column or row) vector;

(2) Matrix transpose;

(3) Scale a matrix (by a constant)

(4) Display a matrix

(5) Matrix addition/subtraction;

(6) Matrix multiplication;

(7) Generate variance-covariance matrix

(8) Matrix inversion;

#2. Decide which data structure to use

Of course I am using vector. Its iterator is big convenience for me and when calculating small-size data, vector is almost as fast as array :)

Here we are considering two dimension matrix (it has rows, columns and that's all). So I used two dimension vectors to save a matrix.

i.e

********

vector<vector<double>> matrix;

********

Every element (vector) in the outer vector represents a row while the elements of the inner vector (double numbers) are the matrix elements contained in the specific row. Below is an example.

0 1 2 3

0 O O O O

1 O O O O

2 O O O O

this is a 3 x 4 metric, the outer vector is a vector that contains three vectors (0,1,2); there are 3 inner vectors contained respectively in the position 0 , 1, 2 . Of the outer vector, each of them contains 4 metric element O,O,O,O

#3 Implementation of matrix operations!

(1) generate a simple vector

******************************************************************************************

vector<vector<double>> Matrix::generator(int num, string type) { vector<vector<double>> single; if (type == "r") { vector<double> temp; for (int i = 0; i < num; i++) { temp.push_back(1); } single.push_back(temp);

} else { for (int i = 0; i < num; i++) { vector<double> temp; temp.push_back(1); single.push_back(temp); } } return single;

}

There are two inputs of this function: type decides if we are going to generate a single row vector or single column vector; the num decides the length of the vector;

if type="r" that means we are going to generate a single row vector with num elements (1), to achieve this we put num 1s in a vector since there are "num" elements in this row, and then put this vector in a vector; The last step seems trivial since there is only one element in the outer vector, however, we have to do this to guarantee all the matrix we generate and compute has the same structure.

---

(2) matrix transpose

********************************************************************************************

vector<vector<double>> Matrix::trans(vector<vector<double>>x){ vector<vector<double>>y; int oldrows = x.size(); int oldcolumns = x[0].size(); for (int j = 0; j < oldcolumns; j++) { vector<double> temp; for (int i = 0; i < oldrows; i++) { temp.push_back(x[i][j]); } y.push_back(temp); } return y;

}

this function takes a metric as input and generate the its transpose metric:

Remember, first, a metric is essentially a vector, so when we implement .size() operation on a matrix, it will return the size of the outer vector, i.e the number of rows in this matrix; second, every row must have the same number (number of columns) of elements, so it will return the number of columns if we implement .size() operation on any sub-`vector of the outer vector, here we take outvector[0].

after obtaining the size of original matrix, we can generate a new matrix as presented above and in (1) but in a different order of row and column to achieve the transpose operation (it is easy) .

---

(3) scale a matrix

*************************************************************************************************

vector<vector<double>>Matrix::scale(double x, vector<vector<double>>xx){ vector<vector<double>>y = xx; vector<vector<double>>::iterator iter = y.begin();

for (; iter != y.end(); iter++) { vector<double>::iterator iter2 = (*iter).begin(); for (; iter2 != ((*iter).end()); iter2++) (*iter2) = (*iter2)*x; } return y; }

This function takes the scale factor and a matrix as inputs and outputs the scaled matrix. The idea behind this function is very easy: we just iterate every element in the matrix and multiply them by the scale factor. The outer for loop goes though every row and inner for loop visits every element in the current row.

---

(4) display a matrix

**********************************************************************************************

void Matrix::show_matrix(vector<vector<double>> x){ vector<vector<double>>::iterator iter = x.begin();

for (; iter != x.end(); iter++) { vector<double>::iterator iter2 = (*iter).begin(); for (; iter2 != ((*iter).end()); iter2++) cout << (*iter2) << " "; cout << endl; }

}

This is nothing more than iterate matrix. each time we move to an element, we output it to screen by standard output "cout" and when we finish the iteration of one row (finish the iteration of every element in a inner vector), we move to next element of outer vector and change line (and refresh the

buffer) by cout<<endl;

---

(5) matrix addition/subtraction

*********************************************************************************************

// addition method

vector<vector<double>> Matrix::subtraction(vector<vector<double>> x, vector<vector<double>> y){ vector<vector<double>> add1; vector<vector<double>> add2; vector<vector<double>> add_result; add1 = x; add2 = y; int rows = add1.size(); int columns = add1[0].size(); for (int i = 0; i < rows; i++) { vector<double>temp; for (int j = 0; j < columns; j++) { temp.push_back(add1[i][j] - add2[i][j]); } add_result.push_back(temp);

}

return add_result; }

To perform +/- on two matrixes, we simply iterate them at the same time (since they must have same size of rows and columns) and +/- their correspond element and put the result in the same position in a new matrix.

---

(6) matrix multiplication

***********************************************************************************************

vector<vector<double>> Matrix::multiplication(vector<vector<double>> x, vector<vector<double>> y){ vector<vector<double>> multi1 = x; vector<vector<double>> multi2 = y; vector<vector<double>> multi_result; int rowsof1 = multi1.size(); int columnsof2 = multi2[0].size(); int columnsof1 = multi1[0].size(); for (int i = 0; i < rowsof1; i++) { vector<double> temprow; for (int j = 0; j < columnsof2; j++) { double sum = 0; for (int k = 0; k < columnsof1; k++) { sum += multi1[i][k] * multi2[k][j]; //cout << sum<<"~"; } temprow.push_back(sum); } multi_result.push_back(temprow); }

return multi_result; }

A (Nxk) multipy by (kxM) will yield a NxM matrix. So in order to generate the result matrix, we have to first gauge its size by obtaining the row size of the first matrix and the column of the second size by method introduced in (2).

To achieve our multiplication purpose, we have to go through 3 layer for loop: the our two loops will iterate (rows and columns) and put the result element in each position of the matrix; the most inner loop will compute the result in that position----recall the knowledge of linear algebra, each element in the new matrix is the dot product of a single line vector:

i.e

O

[O, O, O] .[ O ]= X

O

So, we definite need a loop to achieve this :)

---

(7) generate variance-covariance matrix

**********************************************************************************************************

vector<vector<double>>Matrix::covariance(vector<vector<double>>x){ vector<vector<double>>ratevec = x; for (int m = 0; m < ratevec.size(); m++) { double sum = 0; for (int i = 0; i < ratevec[0].size(); i++) { sum += ratevec[m][i]; } sum /= ratevec[0].size(); vector<double>temp; temp.push_back(sum); average.push_back(temp); }

vector<vector<double>> covmatrix; for (int i = 0; i < ratevec.size(); i++) { vector<double>temp; for (int j = 0; j < ratevec.size(); j++) { double sum = 0; for (int k = 0; k < ratevec[0].size(); k++) { sum += (ratevec[i][k] - average[i][0])*(ratevec[j][k] - average[j][0]); } temp.push_back(sum / ratevec[0].size()); } covmatrix.push_back(temp); } return covmatrix;

}

This is a big guy....

First we have to know the formula to compute covariance cov(x,y)=E[(x-E[x])(y-E[y])], when the sample size is large enough according to large number theorem, it is roughly ok to substitute Expectation with "average" and this is what we are going to do.

Our input is a matrix, each row of which represents an asset, while the elements in that row is the market data of that asset. For example, if we have 40 asset in portfolio with 250 daily data of 1 year then we have a 40X250 input matrix.

As the formula implied, we have to compute the average of each asset, this is what I did in the frist two nested for loop. Note the result is saved in a matrix called "average", which is actually a public variable of the matrix class so I didn't declare it in the function body.

Now we have a matrix of average, what to do next is to apply the formula and to generate the variance-covariance matrix. To do this, I performed a 3-layer for loop, the outer two loops generate the matrix body (i.e 40 x 40), the inner loop compute the variance or covriance using the the value of asset datas and information contained in the "average matrix"

---

(8) matrix inversion

************************************************************************************************

vector<vector<double>> Matrix::inverse(vector<vector<double>>vec) { vector<vector<double>> raw; vector<vector<double>> result; raw = vec; size = vec.size(); int i, j; data[30][60] = { 0 }; for (i = 0; i<size; i++) { vector<double> temp = raw[i]; for (j = 0; j<2 * size; j++) { if (j<size) { data[i][j] = temp[j]; } else if (j == i + size) data[i][j] = 1.0; else data[i][j] = 0.0; } } int ii, jj, kk; int maxI = 0; for (ii = 1; ii<size; ii++) { if (fabs(data[maxI][0])<fabs(data[ii][0])) maxI = ii; } if (maxI != 0) { double temp; for (jj = 0; jj<2 * size; jj++) { temp = data[0][jj]; data[0][jj] = data[maxI][jj]; data[maxI][jj] = temp; } } double temp2; for (ii = 0; ii<size; ii++) { if (data[ii][ii] != 0) temp2 = 1.0 / data[ii][ii]; else { cout << "this matrix has no inverse !" << endl;

} for (jj = 0; jj<2 * size; jj++) data[ii][jj] *= temp2; for (jj = 0; jj<size; jj++) { if (jj != ii) { double temp3 = data[jj][ii]; for (kk = 0; kk<2 * size; kk++) data[jj][kk] -= temp3*data[ii][kk]; } } } int m, n; for (m = 0; m < size; m++) { vector<double> temp; for (n = size; n<2 * size; n++) { temp.push_back(data[m][n]); } result.push_back(temp); } return result; }

This may well be the most complex part of this matrix (and actually it is not difficult at all).Here is a reference for matrix inversion

http://www.mathwords.com/i/inverse_of_a_matrix.htm

Basically, I was using the Augmented matrix method, which is explained in that web page.

actually this is the first function I wrote in the matrix computation class, I was intened to increase the speed by using build-in type array...this cause me sometime..and I am feeling lazy to narrative it again :(

Explore it! please!

---