Matrices

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:

[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:

[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:

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:

These are the properties available for matrices:

And these are the supported methods:

These are the operators available for matrices:

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

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:

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.

Optimisations

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:

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:

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:

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:

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:

mat[, 0]

Omitting the column number yields a whole row:

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:

mat[^1, ^1]

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

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:

-- 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:

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:

> 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:

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:

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:

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

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:

-- 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:

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

Reference

Other Resources