AUSTRA

AUSTRA matrices are represented by the matrix class. They are implemented as row-first, double precision dense matrices.

The AUSTRA matrix class is based on three different C# structures: Matrix, LMatrix, and RMatrix. The compiler takes automatically care of any conversions when needed.

A matrix can be constructed by enclosing its components inside brackets:

Austra

`[1, 2, 3; 2, 3, 4; 3, 4, 5]`

Rows must be separated by semicolons (;), and items in a row must be separated by commas. This syntax does not allows writing a matrix with only one row, since the compiler would not be able to tell it from a vector. A workaround is writing a matrix with only one column, and transposing it:

Austra

`[1; 2; 3; 4]'`

You can also create a new matrix by concatenating two existing matrices, or a matrix and a vector. You can use either vertical or horizontal concatenation:

Austra

**let** m = [1,2;3,4], v=[1, 1];
-- Horizontal concatenation (2x4 matrix).
[m, m];
-- Horizontal concatenation (2x3 matrix).
[m, v];
[v, m];
-- Vertical concatenation (4x2 matrix).
[m; m];
-- Vertical concatenation (3x2 matrix).
The vector is handled as a row vector.
[m; v];
[v; m];

These class methods are available for creating matrices:

matrix::new | Overloaded constructor (see below). |

matrix::rows | Creates a matrix given its rows as vectors. |

matrix::cols | Creates a matrix given its cols as vectors. |

matrix::eye | Creates an identity matrix given its size. |

matrix::diag | Creates a diagonal matrix given the diagonal as a vector. |

matrix::random | Creates a matrix with random values from a uniform distribution. |

matrix::nrandom | Creates a matrix with random values from a normal standard distribution. |

matrix::lrandom | Creates a lower-triangular matrix with random values from a uniform distribution. |

matrix::lnrandom | Creates a lower-triangular matrix with random values from a standard normal distribution. |

matrix::cov | Creates a covariance matrix given a list of series. |

matrix::corr | Creates a correlation matrix given a list of series. |

These are the properties available for matrices:

amax | Gets the absolute maximum. See AMax. |

amin | Gets the absolute minimum. See AMin. |

chol | Calculates the Cholesky factorization. See Cholesky. |

cols | Gets the number of columns. See Cols. |

det | Calculates the determinant. See Determinant. |

diag | Gets the main diagonal as a vector. See Diagonal. |

evd | Calculates the EigenValues Decomposition. See EVD. |

inverse | Gets the inverse of this matrix. See Inverse. |

isSymmetric | Verifies if the matrix is a symmetric one. See IsSymmetric. |

max | Gets the maximum value from the cells. See Maximum. |

min | Gets the minimum value from the cells. See Minimum. |

rows | Gets the number of rows. See Rows. |

trace | Gets the sum of the main diagonal. See Trace. |

stats | Returns statistics on cells. See Stats. |

And these are the supported methods:

all | Checks if all cells satisfy a lambda predicate. See All. |

any | Checks if exists a cell satisfying a lambda predicate. See Any. |

getCol | Gets a column by its index. See GetColumn. |

getRow | Gets a row by its index. See GetRow. |

map | Creates a new matrix with transformed cells. See Map. |

redim | Creates a new matrix a different size. See Redim. |

These are the operators available for matrices:

+ | Adds two matrices, or a matrix and a scalar. |

- | Subtracts two matrices, or a matrix and a scalar. It is also used as a unary operator. |

* | Matrix * matrix = matrix multiplication. Matrix * number = matrix scale. Matrix * vector = vector transformation. Vector * matrix = vector is transposed and then transformed. |

.* | Pointwise multiplication of two matrices. |

./ | Pointwise quotient of two matrices. |

/ | Divides a matrix by a scalar, but also divides either a vector or a matrix by a matrix, for solving linear equations. |

' | Unary sufix operator for matrix transpose. |

This examples shows how to solve linear equations for a vector, using division by a matrix:

Austra

**let** m = **matrix::** random(5) + 0.01,
v = **vec::** random(5),
answer = v / m **in**
m * answer - v

Solving equations for a matrix is also possible:

Austra

**let** m = **matrix::** random(5) + 1;
**matrix::** eye(5) / m - m.inverse

Internally, the LU factorization of the matrix is used for equation solving, for the general case. When the matrix at the left is a triangular matrix, a most efficient algorithm is used.

The compiler performs some optimisations for matrix operations. For instance, these two expressions yield the same result, but the second one avoids one matrix transpose:

Austra

**let** m = **matrix::** random(10), v = **vec::** random(10);
-- Transpose a matrix and then transform a vector:
m' * v;
-- Changing the order of operands avoids a transpose:
v * m;

The compiler also detects when the second matrix in a matrix multiplication is being transposed:

Austra

**let** m1 = **matrix::** random(10), m2 = **matrix::** random(10) **in**
m1 * m2'

This pattern is handled by the MultiplyTranspose method, which not only saves the time spent in the transpose but also avoids a temporal memory allocation.

These two operation patterns are also detected and implemented with a single method call:

Austra

**let** m = **matrix::** random(10);
**let** v1 = **vec::** random(10);
**let** v2 = **vec::** random(10);
**let** scaleFactor = 0.1;
m * v1 ± v2;
m * v1 ± scaleFactor * v2;

These special operations are implemented by the MultiplyAdd and MultiplySubtract group of overloaded methods.

Individual cells are accessed using the row and column inside brackets:

Austra

```
mat[0, 0];
mat[mat.rows - 1, mat.cols - 1]
```

All indexes starts from zero. If the row index is omitted, a whole column is returned:

Austra

`mat[, 0]`

Omitting the column number yields a whole row:

Austra

```
mat[0,];
mat[0]
```

Carets can also be used in any of the two indexes, to count positions from the end. For instance, this expression returns the rightmost lower cell of the matrix:

Austra

`mat[^1, ^1]`

Columns and rows can also be extracted as vectors using relative indexes:

Austra

`mat[, ^2]; `-- Next to last column.
mat[^2,] -- Next to last row.

Ranges are accepted for both dimensions, and can be combined with indexes too:

Austra

-- Remove last row and last column.
mat[0..^1, 0..^1];
-- Last row without first and last items.
mat[^1, 1..^1]

An eigenvalue λ and its associated eigenvector ν are any pair of values that satisfy this equation for a square matrix M:

*Mν=λv*

It means that, when the matrix transforms an eigenvector, the result is the same vector, except for a scale factor.

The Eigenvalue Decomposition of a matrix is an algorithm that identifies all the pairs of eigenvalues and eigenvector for a given square matrix. You can efficiently find eigenvalues and eigenvectors in AUSTRA applying the evd method on a matrix:

Austra

**let** m = **matrix::** random(10);
**let** e = mat.evd;
-- This is a matrix with all eigenvectors as columns.
e.vectors;
-- This is a complex vectors with all eigenvalues.
e.values;
-- An alternative representation of eigenvalues, using a real matrix.
e.d

A real matrix can have both real and complex eigenvalues, and that is why the values property of the decomposition is a complex vector. When there exists a complex eigenvalue, its complex conjugate is also an eigenvalue of the matrix. Aside from values, the eigenvalues are also returned in a d property, which is a real block diagonal matrix. Each real eigenvalue is placed in the main diagonal, and complex eigenvalues are represented as 2x2 blocks in the diagonal.

The best way to visualize how evd returns eigenvalues is to show an example. These are the eigenvalues of a 4x4 random matrix:

Austra

```
> e.values
ans ∊ ℂ(4)
<2.02631; 0>
<-0.229546; 0.093745>
<-0.229546; -0.093745>
<0.364586; 0>
> e.d
ans ∊ ℝ(4⨯4)
2.02631 0 0 0
0 -0.229546 0.093745 0
0 -0.093745 -0.229546 0
0 0 0 0.364586
```

The second and third eigenvalues are a conjugated pair of complex numbers. See how they are represented in the block diagonal matrix, using four cells.

When *m* is a square matrix, the following mathematical equivalence must hold:

Austra

`m * m.evd.vectors = m.evd.vectors * m.evd.d`

In practice, however, we must take the lost of precision into account. For verifying an EVD, you can use this formula in AUSTRA:

Austra

**let** m = matrix::random(32), e = m.evd **in**
(m * e.vectors - e.vectors * e.d).amax <= 1e-12;
**let** lm = matrix::lrandom(32), m = lm * lm', e = m.evd **in**
(m * e.vectors - e.vectors * e.d).amax <= 1e-12;

Symmetric matrices, as the one generated for the second example above, are decomposed using a more efficient algorithm, so Austra checks symmetry first, before applying any of these algorithms.

The LU (lower-upper) factorization algorithm takes a square matrix and generates a lower-triangular matrix L and an upper-triangular U that, when multiplied, regenerates the original matrix. The algorithm may also reorder rows for the sake of numerical stability.

The lu property, when applied to a square matrix, returns an LU structure from the Austra library that provides these properties:

det | Gets the determinant of the original matrix. |

lower | Gets the lower-triangular matrix from the factorization. |

perm | Gets an integer vector with the permutations. |

size | Gets the number of rows and columns from the original matrix. |

upper | Gets the upper-triangular matrix from the factorization. |

More relevant for us are the two overloads of the solve method of this structure:

solve(vec) | Solves the equations Mx=v or Mx=m, where M is the factorised matrix, v is a vector and m is another matrix. |

The solve method from the original matrix does the same work, but it delegates the solution to a LU factorization created on the fly. If we need to solve more than one equation involving the same matrix at the left side, it is more efficient to perform the LU factorisation once, and reuse the result for each linear system, as this code show:

Austra

-- This a matrix whose LU factorisation will be reused.

**let** m = **matrix::** random(10);
-- The LU factorisation is computed here.

**let** lu = m.lu;
-- Now we generate a vector and a matrix for the right sides.

**let** n = **matrix::** random(m.rows), v = **vec::** random(m.rows);
-- Solve m*x=n and m*y=v, and check the accuracy of the results.
(m * lu.solve(n) - n).amax;
(m * lu.solve(v) - v).amax;

The accuracy, for the above example, is near 1e-16 or 1e-15.

Another important factorisation is the Cholesky decomposition. It requires a square matrix, but this time, the matrix must be a symmetric one. The Cholesky algorithm finds, for a given *M* matrix, a lower-triangular matrix *C* such that:

*M=C∙C'*

Of course, matrices generated by multiplying a lower-triangular matrix by its transposed are symmetric, as is easy to demonstrate. Not only that: the resulting matrix must be a positive-definite matrix, with its determinant greater than zero since the determinant of a matrix product is the product of the determinants. *C*, when exists, can be considered as a sort of square root of the original matrix *M*.

A matrix provides two properties related to the Cholesky decomposition. The chol property returns the lower-triangular matrix when it exists or throws an exception otherwise. The cholesky property, on the other hand, returns an object that encapsulates the Cholesky matrix. The reason for having this apparent detour is that the returned object implements these two overloads of a solve method:

solve(vec) | Solves the equations Mx=v or Mx=m, where M is the factorised matrix, v is a vector and m is another matrix. |

The lower-triangular matrix computed by the decomposition can also be retrieved from the object return by cholesky using its lower property.