### Introduction

Application of various branches of mathematics is essential in writing complex automated trading systems. One of such branches is linear algebra.

There are currently no extensive publicly available libraries in MQL4 that implement various methods of linear algebra (particularly, the theory of matrices and determinants).

This article describes the LibMatrix library in MQL4 that includes the implementation of the most common matrix operations.

Four months after being written in C++ the entire code was rewritten in MQL4 with some modifications.

### 1. Library Structure

Let's first consider some peculiarities of working with matrices in the proposed library.

First of all, matrices are stored in one-dimensional arrays. It has to do with the fact that even though MQL4 offers features for creating multi-dimensional arrays, it is only possible to change the size of the first dimension. To store a matrix of N rows and M columns, we need a one-dimensional N*M array. The matrix is packed into an array row by row, i.e. the first row elements in a data array are followed by the second row elements and so on.

Second, when creating the library it was decided not to store additional information on matrix size in a data array as sizes are integer values, while elements are real numbers (storing sizes as real numbers could dramatically affect the operation of the library).

So a matrix in the format processed by the proposed library consists of three variables: a one-dimensional array of the **double** type and two integer variables (e.g. of the **int** type) for storing information on the number of matrix rows and columns.

The library consists of two files: LibMatrix.mqh and LibMatrix.mq4. The first file is to be included as required. It contains function prototypes that can be imported from the library. The second file contains function implementations (program code). That is what we are going to talk about further below.

The functions implemented in the library can be divided into several groups:

- general mathematical functions (
**MathLerp**,**MathInRangeRandom**,**MathDoublesEqual**) - auxiliary functions for working with matrices (
**MatrIndiciesToOffset**,**MatrCopy**,**MatrSetSize**,**MatrResize**) - functions for working with rows and columns (
**MatrSwapRows**,**MatrSwapCols**,**MatrCopyRow**,**MatrCopyCol**,**MatrRowIsZero**,**MatrColIsZero**) - functions for general check of conditions (
**MatrIsSquare**,**MatrIsElemOnMainDiagonal**,**MatrCompatiblityCheckAdd**,**MatrCompatiblityCheckMul**) - content initialization functions (
**MatrLoadZero**,**MatrLoadIdentity**,**MatrLoadInRangeRandom**) - functions for checking conditions by content type (
**MatrIsZero**,**MatrIsDiagonal**,**MatrIsIdentity**,**MatrIsSymmetric**,**MatrIsAntisymmetric**,**MatrEqual**) - functions for element-wise scalar operations (
**MatrAddScalar**,**MatrSubScalar**,**MatrMulByScalar**,**MatrDivByScalar**) - functions for basic matrix operations (
**MatrAddMatr**,**MatrSubMatr**,**MatrMulMatr**,**MatrTrace**,**MatrTranspose**) - other functions (
**MatrGaussianElimination**,**MatrGJBatchSolve**,**MatrMinor**,**MatrAlgebraicComplement**,**MatrInvertUsingMinors**,**MatrInvertUsingGJ**,**MatrDet**,**MatrDetTriang**,**MatrRank**,**MatrRankTriang**,**MatrComputeConnectedMatr**,**MatrLerpMatr**) - input/output functions (
**MatrPrint**,**FileWriteMatr**,**FileReadMatr**)

Let's have a closer look at each of the groups.

### 2. Function Descriptions

**2.1. General Mathematical Functions**

We will first consider the group of general mathematical functions. It includes functions that are not directly associated with matrices but are used in the library.

double MathLerp(double rangeLowLimit, double rangeHighLimit, double balance);

This function performs linear interpolation between two values using the following formula: **rangeLowLimit + balance * (rangeHighLimit - rangeLowLimit)**. The **balance** parameter value should be within the range [0;1].

double MathInRangeRandom(double rangeLowLimit, double rangeHighLimit);

This function returns uniformly distributed random numbers from the range [**rangeLowLimit**;**rangeHighLimit**]. The standard **MathRand** function is used for generation.

bool MathDoublesEqual(double value1, double value2, double tolerance);

This function is used to compare two values, **value1** and **value2**, of the **double** type with the specified accuracy (**tolerance**). Use of approximate comparisons is required due to possible loss of accuracy. The **tolerance **parameter value should not be negative (as the default value, you can use the **DEFAULT_TOLERANCE** constant defined in the library).

**2.2. Auxiliary Functions for Working with Matrices**

Let's proceed to the group of auxiliary functions for working with matrices. This group covers functions intended for simplification of operations related to packing matrices in one-dimensional arrays.

int MatrIndiciesToOffset(int row, int col, int numRows, int numCols)

This function calculates the offset of an element with respect to the beginning of the array in which the matrix is packed. Among the parameters it receives are the number of the row (**row**) and the column (**col**) at the intersection of which the element lies, as well as the matrix size (**numRows** - number of matrix rows, **numCols** - number of matrix columns). The **row** and **col** values should lie within the ranges [0;**numRows**-1] and [0;**numCols**-1], respectively.

void MatrCopy(double& src[], double& dst[]);

This function copies all elements of the **src** matrix to the **dst** matrix. Matrix data arrays are the only parameters passed. Upon completion of copying the **dst** matrix size becomes equal to the **src** matrix size.

void MatrSetSize(double& matr[], int numRows, int numCols);

This function changes the size of the **matr **data array so that it can store a matrix of **numRows** rows and **numCols** columns. Matrix data integrity is not guaranteed. The **numRows** and **numCols **parameter values should be strictly positive.

void MatrResize(double& matr[], int numRowsOld, int numColsOld, int numRowsNew, int numColsNew);

This function changes the size of the **matr** data array of the existing matrix. The **numRowsOld** and **numColsOld** parameters must be equal to the initial matrix size. New size is set by the **numRowsNew** and **numColsNew** parameters. Matrix data integrity is guaranteed subject to new size, except for the cases where one or both matrix dimensions are decreased. New elements are initialized to zeros.

**2.3. Functions for Working with Rows and Columns**

Let's have a look at the functions for working with matrix rows and columns.

void MatrSwapRows(double& matr[], int numRows, int numCols, int row1, int row2);

This function swaps the content of the **row1** and **row2** rows for the **matr** matrix consisting of **numRows** rows and **numCols** columns. The **row1** and **row2** parameter values should lie within the range [0;**numRows**-1].

void MatrSwapCols(double& matr[], int numRows, int numCols, int col1, int col2);

This function swaps the content of the **col1** and **col2** columns. The **col1** and **col2** parameter values should lie within the range [0;**numCols**-1].

void MatrCopyRow(double& matr[], int numRows, int numCols, int rowToCopy, int rowToReplace);

This function copies the content of the **rowToCopy** row to the **rowToReplace** row. The **rowToCopy** and **rowToReplace** parameter values should be within the range [0;**numRows**-1].

void MatrCopyCol(double& matr[], int numRows, int numCols, int colToCopy, int colToReplace);

This function copies the content of the **colToCopy** column to the **colToReplace** column. The **colToCopy** and **colToReplace** parameter values should lie within the range [0;**numCols**-1].

bool MatrRowIsZero(double& matr[], int numRows, int numCols, int row, double tolerance);

This function checks whether the **row** row is null with the **tolerance** degree of accuracy. If the row is null, the returned value will be **true**.

bool MatrColIsZero(double& matr[], int numRows, int numCols, int col, double tolerance);

This function checks whether the **col** column is null with the **tolerance** degree of accuracy. If the column is null, the returned value will be **true**.

**2.4. Functions for General Check of Conditions**

The group of functions for general check of conditions is intended for improvement of the code readability.

bool MatrIsSquare(int numRows, int numCols);

This function checks whether the matrix of **numRows** rows and **numCols** columns is square (i.e. it checks whether **numRows**=**numCols**). It returns **true** for square matrices.

bool MatrIsElemOnMainDiagonal(int row, int col);

This function checks whether the element with the [**row**][**col**] indices is on the main diagonal of the matrix. It returns **true** for the elements on the main diagonal.

bool MatrCompatiblityCheckAdd(int numRows1, int numCols1, int numRows2, int numCols2);

This function checks compatibility of two matrices for purposes of addition-like operations. Dimensions are checked pairwise for equality (**numRows1**=**numRows2** and **numCols1**=**numCols2**); **true** is returned in case of equality.

bool MatrCompatiblityCheckMul(int numRows1, int numCols1, int numRows2, int numCols2);

This function checks compatibility of two matrices for purposes of multiplication. It checks whether the number of columns in the first matrix (**numRows1** by **numCols1** matrix) is equal to the number of rows in the second matrix (**numRows2** by **numCols2** matrix). Multiplication is possible if **true** is returned.

**2.5. Matrix Content Initialization Functions**

The group of matrix content initialization functions is intended for loading frequently used matrices (e.g. null and identity matrices). Since memory is not allocated, the **MatrSetSize** or **MatrResize** functions should be called before calling the group functions.

void MatrLoadZero(double& matr[], int numRows, int numCols);

This function initializes the **matr** data array to zeros thus creating a null matrix of **numRows **rows and **numCols** columns.

void MatrLoadIdentity(double& matr[], int numRows, int numCols);

This function initializes all elements of the **matr** data array to zeros, except for those that are on the main diagonal of the matrix. The initialized matrix of **numRows** rows and **numCols** columns should be square.

void MatrLoadInRangeRandom(double& matr[], int numRows, int numCols, double rangeLowLimit, double rangeHighLimit);

This function initializes the data array of the **matr** matrix to random numbers from the range [**rangeLowLimit**;**rangeHighLimit**] (random numbers are generated using the **MathInRangeRandom** function).

**2.6. Functions for Checking Conditions by Content Type**

Let's now consider the group of functions for checking conditions by content type. The functions of this group perform a comparison with the **tolerance** degree of accuracy (this parameter value must not be negative) and return **true** if the condition under consideration is true.

bool MatrIsZero(double& matr[], int numRows, int numCols, double tolerance);

This function checks if the matrix is null.

bool MatrIsDiagonal(double& matr[], int numRows, int numCols, double tolerance);

This function checks if the matrix is diagonal. The matrix will be considered diagonal if it is square and all its elements outside the main diagonal are zero.

bool MatrIsIdentity(double& matr[], int numRows, int numCols, double tolerance);

This function checks if the matrix is identity.

bool MatrIsSymmetric(double& matr[], int numRows, int numCols, double tolerance);

This function checks if the matrix is symmetric (i.e. a[i][j]=a[j][i] is true for any i and j).

bool MatrIsAntisymmetric(double& matr[], int numRows, int numCols, double tolerance);

This function checks if the matrix is skew-symmetric/antisymmetric (i.e. a[i][j]=-a[j][i] is true for any i and j and a[k][k]=0 is true for any k).

bool MatrEqual(double& matr1[], double& matr2[], int numRows, int numCols, double tolerance);

This function checks element-wise equality between matrices with the same dimensions.

**2.7. Functions for Element-wise Scalar Operations**

Let's proceed to the group of functions for element-wise scalar operations. This group is peculiar in that function results are not placed in the initial array of the **matr** matrix but rather into the new **result** array. However, if the same array is passed as both the **result** and **matr** parameters, the results will be placed to the initial matrix.

void MatrAddScalar(double& matr[], int numRows, int numCols, double scalar, double& result[]);

This function adds the **scalar** value to each matrix element.

void MatrSubScalar(double& matr[], int numRows, int numCols, double scalar, double& result[]);

This function subtracts the **scalar** value from each matrix element.

void MatrMulByScalar(double& matr[], int numRows, int numCols, double scalar, double& result[]);

This function multiplies each matrix element by the **scalar** value.

void MatrDivByScalar(double& matr[], int numRows, int numCols, double scalar, double& result[]);

This function divides each matrix element by the **scalar** value.

**2.8. Functions for Basic Matrix Operations**

The group of functions for basic matrix operations consists of functions for performing elementary actions - calculations of sum, product and trace, as well as transposition operations. The results of function operation (except for **MatrTrace**) are returned to the **result** array.

void MatrAddMatr(double& matr1[], double& matr2[], int numRows, int numCols, double& result[]);

This function adds up two matrices with the same dimensions.

void MatrSubMatr(double& matr1[], double& matr2[], int numRows, int numCols, double& result[]);

This function subtracts the **matr2** matrix from the **matr1** matrix. The matrices should be of equal size.

void MatrMulMatr( double& matr1[], int numRows1, int numCols1, double& matr2[], int numRows2, int numCols2, double& result[], int& numRowsRes, int& numColsRes);

This function multiplies the **matr1 **matrix by the **matr2** matrix. The matrices should be compatible for multiplication, i.e. the number of columns (**numCols1**) in the first matrix should be equal to the number of rows (**numRows2**) in the second matrix. The size of the resulting matrix is returned in variables references to which are passed in the **numRowsRes **and **numColsRes** parameters.

double MatrTrace(double& matr[], int numRows, int numCols);

This function calculates the matrix trace (sum of the diagonal elements). The matrix must be square.

void MatrTranspose(double& matr[], int numRows, int numCols, double& result[], int& numRowsRes, int& numColsRes);

This function transposes the matrix (swaps columns with rows). The new matrix size is returned in variables references to which are passed in the **numRowsRes** and **numColsRes** parameters.

**2.9. Other Functions**

The group of other function offers possibilities for matrix inversion, calculation of ranks and determinants, etc.

int MatrGaussianElimination(double& matr[], int numRows, int numCols, double& result[], double tolerance);

This function represents the implementation of Gaussian elimination with pivoting. The initial **matr** matrix of **numRows **rows and **numCols** columns is reduced to a trapezoidal (triangular) form, with the results being placed in the **result** array.

The returned value is the number of linearly independent rows of the matrix (i.e. its rank). Since calculations can lead to function rounding errors, it is required to pass a non-negative **tolerance** parameter value that defines the degree of accuracy of the comparison.

bool MatrGJBatchSolve( double& matr[], int numRows, int numCols, double& rhs[], int numRowsRHS, int numColsRHS, double& roots[], int& numRowsRoots, int& numColsRoots, double tolerance );

This function represents the implementation of Gauss-Jordan elimination intended for solving sets of several equation systems with one and the same system matrix.

The common system matrix is passed in the **matr** (data array), **numRows **and **numCols** (dimensions) parameters. A set of right hand side vectors is also passed as a matrix (the **rhs** parameter), where each column represents a right hand side of the system under consideration.

In other words, the **numRowsRHS **parameter contains the number of equations in the system (number of rows in the system matrix), while the **numColsRHS** parameter contains the number of right hand side vectors (number of systems to be solved).

The results are returned in the form of the **roots** matrix of **numRowsRoots** rows and **numColsRoots** columns. The solution for each system of equations is placed in the matrix column and corresponds to the solution for the system of the given system matrix and the column with the relevant number from the **rhs** matrix. The function returns **true** if the solution for the set of equation systems has been found.

void MatrMinor( double& matr[], int numRows, int numCols, int rowToExclude, int colToExclude, double& result[], int& numRowsRes, int& numColsRes);

This function finds the minor of the matrix lying at the intersection of the **rowToExclude** row and the **colToExclude** column. The resulting matrix elements are returned in the **result** array, while the dimensions are returned in variables references to which are passed in the **numRowsRes** and **numColsRes** parameters.

double MatrAlgebraicComplement(double& matr[], int numRows, int numCols, int elemRow, int elemCol, double tolerance);

This function calculates the algebraic cofactor (minor with a place sign) for the element lying at the intersection of the **elemRow** row and **elemCol** column. To calculate the determinant, we use Gaussian elimination (in the form of the **MatrGaussianElimination** function) with the degree of accuracy set by the **tolerance** parameter.

bool MatrInvertUsingMinors(double& matr[], int numRows, int numCols, double& result[], double tolerance);

This function calculates the inverse of the square **matr** matrix using a union matrix (generated from algebraic cofactors of the corresponding matrix elements). Asymptotic complexity of the employed algorithm: O(N^5). The union matrix is calculated using the **MatrAlgebraicComplement** function (which in turn employs Gaussian elimination) and for this reason the function needs to be passed a non-negative **tolerance** parameter value that defines the degree of accuracy. The elements of the resulting matrix are placed in the **result** array. If the inverse of the matrix is successfully calculated (the determinant is other than zero), the function returns **true**.

bool MatrInvertUsingGJ(double& matr[], int numRows, int numCols, double& result[], double tolerance);

This function calculates the inverse of the given square **matr** matrix by adjoining an identity matrix, i.e. by solving N (N=**numRows**=**numCols**) systems of linear equations using the **MatrGJBatchSolve** function (the matrix of right hand sides is identity; the degree of accuracy is set by a non-negative **tolerance** parameter value). Asymptotic complexity of the employed algorithm: O(N^3). The elements of the resulting matrix are placed in the **result** array. If the inverse of the matrix is successfully calculated (the determinant is other than zero), the function returns **true**.

double MatrDet(double& matr[], int numRows, int numCols, double tolerance);

This function calculates the determinant of the square **matr** matrix by reducing it to a trapezoidal (triangular) form. For this purpose, we use Gaussian elimination with **tolerance ** degree of accuracy (this parameter must be non-negative). Asymptotic complexity of the algorithm: O(N^3).

double MatrDetTriang(double& matr[], int numRows, int numCols);

This function calculates the determinant of the triangular **matr** matrix.

double MatrRank(double& matr[], int numRows, int numCols, double tolerance);

This function calculates the rank of the **matr** matrix by reducing it to a trapezoidal (triangular) form. For this purpose, we use Gaussian elimination with **tolerance ** degree of accuracy (this parameter must be non-negative). Asymptotic complexity of the algorithm: O(N^3).

double MatrRankTriang(double& matr[], int numRows, int numCols, double tolerance);

This function calculates the rank of the **matr** matrix. The matrix should be reduced to a trapezoidal (triangular) form row by row (e.g. by calling the **MatrGaussianElimination** function). The non-negative **tolerance** value sets the degree of accuracy of the first null row calculation.

void MatrComputeConnectedMatr(double& matr[], int numRows, int numCols, double& result[], double tolerance);

This function calculates a union matrix for the square **matr** matrix. The elements of the resulting matrix are placed in the **result** array. The non-negative **tolerance** value defines the degree of accuracy in calculating algebraic cofactors.

void MatrLerpMatr(double& matr1[], double& matr2[], int numRows, int numCols, double balance, double& result[]);

This function performs element-wise linear interpolation between the **matr1** and **matr2** matrices of the same size. The **balance** parameter represents the linear interpolation coefficient and should lie within the range [0;1]. The elements of the resulting matrix are returned in the **result** array.

**2.10. Input/Output Functions**

The group of input-output functions is intended for saving/loading matrices, as well as writing the debugging output of a matrix to the log.

void MatrPrint(double& matr[], int numRows, int numCols);

This function is intended for writing the debugging output to the log (Expert Advisors tab in the Terminal panel). It outputs the matrix row by row taking into consideration that the log is displayed in the terminal in the inverse order (i.e. new entries appear at the top; the function prints matrix rows starting from the end for the convenience of visual analysis of matrices in the log).

void FileWriteMatr(double& matr[], int numRows, int numCols, int handle);

This function saves a matrix (including its dimensions) to the file whose handle is passed in the **handle** parameter. The file should be open in the **FILE_BIN|FILE_WRITE** mode.

void FileReadMatr(double& matr[], int& numRows, int& numCols, int handle);

This function loads a matrix from a file. First, matrix dimensions are read to the variables references to which are passed in the **numRows** and **numCols** parameters. The **matr** array is then resized according to the matrix size which is followed by loading the content of the matr array from the file. The handle of the file from which the data is read is passed in the **handle** parameter; the file should be open in the **FILE_BIN|FILE_READ** mode.

Now that the reader has become familiar with the function descriptions, we can proceed to solving practical issues using the library.

### 3. Example of Use

Let's have a look at an example of creating a polynomial regression on a series of price values using the proposed library.

The process of creating a polynomial regression consists in finding polynomial coefficients f(x)=a[0]+a[1]*x+...+a[degree]*x^degree of the degree degree. This is carried out by solving a system of linear algebraic equations where elements of the system matrix A[degree+1][degree+1] are defined as follows: A[i][j]=(x[0]^(i+j)+x[1]^(i+j)+...+x[numPoints]^(i+j))/numPoints, while elements of the right-hand side vector B[degree+1][1] are defined using the following formula: B[i]=(y[0]*x[0]^i+y[1]*x[1]^i+...+y[numPoints]*x[numPoints]^i)/numPoints.

To solve the task at hand, we have a script (the **LibMatrixEx.mq4** file in the attached archive) that creates a polynomial and displays it on the initial interval and to its right (i.e. extrapolation). Polynomial values on the extrapolation interval can be used to predict the direction of the price movement.

The script is controlled using three vertical lines: two of them are intended for selecting the interval to be analyzed and the third one sets the rightmost point up to which the polynomial is displayed.

To enable the operation of the script, you need to drag it onto the chart and set the required parameters: **delay** - chart refresh rate (in ms), **degree** - degree of polynomial, **linesMargin** - initial distance between the control lines, **linesWidth** - line width of the polynomial chart. You can also select colors for the vertical control lines (**colVLineInt **and **colVLineExt** parameters) and chart lines (**colInt** and **colExt** parameters).

Script operation example

Let's review the key function of the script that deals with matrices.

// creating polynomial regression bool Regression(double& x[], double& y[], int numPoints, int polyDegree, double& poly[]) { // create system matrix double A[]; int numRowsA = polyDegree + 1; int numColsA = polyDegree + 1; MatrSetSize(A, numRowsA, numColsA); // fill the matrix for (int row = 0; row < numRowsA; row++) { for (int col = 0; col < numColsA; col++) { int offset = MatrIndiciesToOffset(row, col, numRowsA, numColsA); A[offset] = SXY(x, y, numPoints, row + col, 0); } } // create a right hand side vector double B[]; int numRowsB = polyDegree + 1; int numColsB = 1; MatrSetSize(B, numRowsB, numColsB); // fill the right hand side vector for (row = 0; row < numRowsB; row++) { offset = MatrIndiciesToOffset(row, 0, numRowsB, numColsB); B[offset] = SXY(x, y, numPoints, row, 1); } // solve a system of linear algebraic equations int numRowsX, numColsX; bool status = MatrGJBatchSolve( A, numRowsA, numColsA, B, numRowsB, numColsB, poly, numRowsX, numColsX, DEFAULT_TOLERANCE ); if (!status) { Print("Error solving the system"); } return (status); }

The function takes on 5 arguments. The first two are **x** and **y** arrays that are intended for storing the coordinates of the points through which the polynomial is plotted. The number of the points is passed in the **numPoints** parameter. The fourth argument sets the degree of the polynomial under consideration (it must be less than the number of points by at least 1). The fifth argument is the reference to the array that receives polynomial coefficients.

At the beginning of the function, the** A** matrix is created, resized and filled using the above formulas. To reference the matrix elements through two-dimensional indices, we use the **MatrIndiciesToOffset** function that calculates one-dimensional offset with respect to the beginning of the array.

The vector-column **B** is then filled in a similar way. This is followed by a call of the **MatrGJBatchSolve** function that solves the created system, thus finding the polynomial coefficients.

The rest of the code that checks line positions and outputs the polynomial should not be of any difficulty to the reader.

### Conclusion

The article has introduced the library for working with matrices and discussed the implemented functions and their peculiarities.

The use of the library has been demonstrated by the example of creating a polynomial regression on a series of candlestick closing price values.

Translated from Russian by MetaQuotes Software Corp.

Original article: https://www.mql5.com/ru/articles/1365

**Attached files**|