Operations on Matrices

Extreme Numerics.NET provides methods for all common operations on matrices. Where applicable, overloaded operator methods are supplied for languages that support them, including element-wise operators in F#. In addition, a series of instance methods have been defined. These provide more efficient use of resources. All instance methods return a reference to the (modified) instance. This enables their use in compound expressions.

All operations can be performed on any combination of matrix types. The return value is always an instance of type Matrix<T>. The actual return value is of a derived type that is chosen for optimal performance.

Most operations come in 3 variants:

A static method in the Matrix class that names the operation, for example Matrix.Add<T>(Matrix<T>, Matrix<T>). These methods compute the result and return it as a new object. They take the element type of the matrices as a generic type argument. In nearly all cases, the element type can be inferred and the type argument can be omitted.
A static method in the Matrix class that names the operation and a Into suffix, for example Matrix.AddInto<T>(Matrix<T>, Matrix<T>, Matrix<T>). These methods again take the element type as a generic type argument that can usually be omitted. They are similar to the first kind, but take an additional parameter that specifies the matrix that is to hold the result. If this argument is not null, then the result of the operation is copied to this argument and the argument is returned. Otherwise, the method acts exactly like the first method.
An instance method in the Matrix<T> class that names the operation and an InPlace suffix, for example Matrix<T>.AddInPlace(Matrix<T>). These methods perform the operation in-place on the instance and return the instance.

Some operations that involve both vectors and matrices are covered in the section on Operations on vectors.

Many operations allow you to specify a transposition operation to apply to one or both operands before performing the actual operation. The transposition operation is specified using a TransposeOperation value. Its possible values are listed below:

Operators and methods for addition and subtraction

Value	Description
None	No operation is performed. The matrix is used 'as is.'
Transpose	The transpose of the matrix used.
Conjugate	The conjugate of the matrix is used. For real matrices, this is the same as None.
ConjugateTranspose	The conjugate transpose of the matrix is used. For real matrices, this is the same as Transpose.

Transposition and conjugation

The Transpose() returns the transpose of a matrix. The ConjugateTranspose() returns the conjugate transpose. If the elements are real, this is the same as Transpose(). The ConjugateTranspose() returns a matrix where each element is the complex conjugate of the original. If the elements are real, this method has no effect.

Addition and subtraction

Besides the usual binary operators, several common compound operations are supported. The following table lists the main methods for addition and subtraction:

Operators and methods for addition and subtraction

Method	Description
Add(a, b).	Adds the matrix b to the matrix a.
Add(a, op(A), b, op(B)).	Adds the matrix a, transformed by operation op(A) to the matrix b, transformed by operation op(B).
Add(a, f) Add(f, a).	Adds the scalar f to the matrix a.
AddScaled(a, f, b).	Adds the scalarf times the matrix b to the matrix a.
AddScaled(a, op(A), f, b, op(B)).	Adds the matrix a, transformed by operation op(A) to the matrix b, transformed by operation op(B).
AddProduct(a, b, c).	Adds the matrix b times the matrix c to the matrix a.
AddScaledProduct(a, f, b, c).	Adds the scalarf times the matrix b times the matrix c to the matrix a.
AddScaledProduct(a, f, b, op(b), c, op(c)).	Adds the scalarf times the matrix b transformed by the operation op(b) times the matrix c transformed by the operation op(c) to the matrix a.
Subtract(a, b). a.Subtract(b)	Subtracts the matrix b from the matrix a.
Subtract(a, f)	Subtracts the scalar f from the matrix a.
Subtract(f, a).	Subtracts the matrix a from the scalar f.
SubtractProduct(a, b, c).	Subtracts the matrix b times the matrix c from the matrix a.

Multiplication and Division

This section covers multiplication and division of matrices by scalars as well as element-wise operations, cross product and dot product.

Operators and methods for multiplication and division

Method	Description
Multiply(f, a).	Multiplies the matrix a by the scalar f.
Divide(a, f).	Divides the matrix a by the scalar f.
ElementwiseMultiply(a, b).	Multiplies each element of the matrix a by the corresponding element of the matrix b. In F#, the .* operator can be used instead.
ElementwiseConjugateMultiply(a, b).	Multiplies each element of the matrix a by the corresponding element of the matrix b. In F#, the .* operator can be used instead.
ElementwiseDivide(a, b).	Divides each element of the matrix a by the corresponding element of the matrix b. In F#, the ./ operator can be used instead.
ElementwiseDivide(f, a).	Divides the scalar f by each element of the matrix a. In F#, the .* operator can be used instead.
ElementwisePow(a, f).	Raises each element of the matrix a to the power specified the integer f. In F#, the .** operator can be used instead.
ElementwisePow(a, f).	Raises each element of the matrix a to the power specified the scalar f. In F#, the .** operator can be used instead.
ElementwisePow(a, b).	Raises each element of the matrix a to the power specified by the corresponding element in the matrix b. In F#, the .** operator can be used instead.

Broadcasting

Broadcasting refers to the process of creating a matrix from a vector by repeating the vector multiple times, and so broadcasting the vector's values along the rows or columns of the matrix. This operation is performed by the Broadcast(Dimension, Int32) method of the vector, or the BroadcastInto(Dimension, Int32, Matrix<T>) method if a result matrix is supplied. For example, broadcasting a three element vector 2 times along the rows, produces a 3x2 matrix. Broadcasting it along the columns produces a 2x3 matrix:

var v = Vector.Create(1.0, 2.0, 3.0);
var m1 = v.Broadcast(Dimension.Rows, 2);
// m1 = [[ 1, 1 ]
//       [ 2, 2 ]
//       [ 3, 3 ]]
var m2 = v.Broadcast(Dimension.Columns, 2);
// m1 = [[ 1, 2, 3 ]
//       [ 1, 2, 3 ]]

Visual Basic

Dim v = Vector.Create(1.0, 2.0, 3.0)
Dim m1 = v.Broadcast(Dimension.Rows, 2)
' m1 = [[ 1, 1 ]
'       [ 2, 2 ]
'       [ 3, 3 ]]
Dim m2 = v.Broadcast(Dimension.Columns, 2)
' m1 = [[ 1, 2, 3 ]
'       [ 1, 2, 3 ]]

Visual Basic

No code example is currently available or this language may not be supported.

Visual Basic

let v = Vector.Create(1.0, 2.0, 3.0)
let m1 = v.Broadcast(Dimension.Rows, 2)
// m1 = [[ 1, 1 ]
//       [ 2, 2 ]
//       [ 3, 3 ]]
let m2 = v.Broadcast(Dimension.Columns, 2)
// m1 = [[ 1, 2, 3 ]
//       [ 1, 2, 3 ]]
//

Broadcasting in this way is inefficient. Most operations involving two matrices support broadcasting directly. In these methods, one of the matrix arguments is replaced with two new arguments: a vector and a Dimension value that specifies whether the vector's elements should be broadcast along the rows or the columns. In-place variants have the same option. In the example below, we create a random matrix, and center each column by subtracting its mean. We use the generic Subtract<T>(Matrix<T>, Vector<T>, Dimension) method as well as the SubtractInPlace(Vector<T>, Dimension) instance method:

var data = Matrix.CreateRandom(20, 3);
var means = data.AggregateColumns(Aggregators.Mean);
var centered = Matrix.Subtract(data, means, Dimension.Columns);
data.SubtractInPlace(means, Dimension.Columns);

Visual Basic

Dim data = Matrix.CreateRandom(20, 3)
Dim means = data.AggregateColumns(Aggregators.Mean)
Dim centered = Matrix.Subtract(data, means, Dimension.Columns)
data.SubtractInPlace(means, Dimension.Columns)

Visual Basic

No code example is currently available or this language may not be supported.

Visual Basic

let data = Matrix.CreateRandom(20, 3)
let means = data.AggregateColumns(Aggregators.Mean)
let centered = Matrix.Subtract(data, means, Dimension.Columns)
data.SubtractInPlace(means, Dimension.Columns) |> ignore
//

Finally, it is possible to have both operands be broadcast vectors. In this case, the broadcast dimension for the second vector argument is omitted, as it must always be the complement of the first. The code below defines a function that returns a Vandermonde matrix. This is a matrix where each row contains a sequence of increasing powers of a value.

public static Matrix<double> Vandermonde(Vector<double> values, int degree)
{
    var exponents = Vector.CreateRange(degree + 1);
    return Matrix.ElementwisePow(values, Dimension.Rows, exponents);
}

Visual Basic

Public Shared Function Vandermonde(values As Vector(Of Double), degree As Integer) As Matrix(Of Double)
    Dim exponents = Vector.CreateRange(degree + 1)
    Return Matrix.ElementwisePow(values, Dimension.Rows, exponents)
End Function

Visual Basic

No code example is currently available or this language may not be supported.

Visual Basic

let Vandermonde (values : Vector<double>) degree = 
    let exponents = Vector.CreateRange(degree + 1)
    Matrix.ElementwisePow(values, Dimension.Rows, exponents)
//

A note on compound assignment operators

The C# language does not provide for explicit overloading of compound assignment operators. A compound operator is translated into the operator followed by an assignment. This means that an expression like

a += b;

is exactly equivalent to

a = a + b;

This is important when working with reference types such as matrices. To perform the above addition, a new matrix instance must be created to hold the sum of a and b. This new instance is then assigned to a, releasing the original value. This causes an unnecessary memory allocation. For large matrices, this may lead to excessive garbage collection, degrading overall performance.

The recommended way to use the InPlace instance methods discussed earlier.

Relational operators

The operators in this section are all element-wise. They return a Boolean matrix with the result of comparing the corresponding elements of two matrices, or the corresponding element of a matrix and a fixed scalar. Also listed here are the Max and Min methods, which follow the same pattern.

Each relational method comes in three flavors: one comparing a matrix to a matrix, one comparing a scalar to a matrix, and one comparing a matrix to a scalar. The table below lists the main functions.

Relational Operators

Method	Description
EqualTo(a, b).	Returns whether each element of a is equal to the corresponding element of the matrix b. In F#, the .= operator can be used instead.
EqualTo(a, b).	Returns whether each element of a is not equal to the corresponding element of the matrix b. In F#, the .<> operator can be used instead.
LessThan(a, b).	Returns whether each element of a is less than the corresponding element of the matrix b. In F#, the .< operator can be used instead.
LessThanOrEqualTo(a, b).	Returns whether each element of a is less than or equal to the corresponding element of the matrix b. In F#, the .<= operator can be used instead.
GreaterThan(a, b).	Returns whether each element of a is greater than the corresponding element of the matrix b. In F#, the .> operator can be used instead.
GreaterThanOrEqualTo(a, b).	Returns whether each element of a is greater than or equal to the corresponding element of the matrix b. In F#, the .>= operator can be used instead.
Max(a, b).	Returns the largest of two corresponding elements of the matrices a and b.
Min(a, b).	Returns the smallest of two corresponding elements of the matrices a and b.

Basic functions

The table below lists a number of basic functions that don't qualify as operators but also aren't elementary functions.

Relational Operators

Method	Description
Abs<T>.	Returns the absolute value of each element of the matrix.
Floor<T>.	Returns the largest integer less than or equal to each element of the matrix.
Ceiling<T>.	Returns the smallest integer less than or equal to each element of the matrix.
Conjugate<T>.	Returns the complex conjugate of each element of the matrix.
Map.	Applies a function to each element of the matrix.

Basic functions

The table below lists a number of basic functions that don't qualify as operators but also aren't elementary functions.

Relational Operators

Method	Description
Sqrt<T>(Matrix<T>).	Returns the square root of each element of the matrix.
Exp<T>(Matrix<T>).	Returns the exponential function of each element of the matrix.
Log<T>(Matrix<T>).	Returns the natural logarithm of each element of the matrix.
Log10<T>(Matrix<T>).	Returns the base 10 logarithm of each element of the matrix.
Hypot<T>(Matrix<T>, Matrix<T>).	Returns the square root of the sum of the squares of corresponding elements in two matrices.
Sin<T>(Matrix<T>).	Returns the sine of each element of the matrix.
Cos<T>(Matrix<T>).	Returns the cosine of each element of the matrix.
Tan<T>(Matrix<T>).	Returns the tangent of each element of the matrix.
Asin<T>(Matrix<T>).	Returns the inverse sine of each element of the matrix.
Acos<T>(Matrix<T>).	Returns the inverse cosine of each element of the matrix.
Atan<T>(Matrix<T>).	Returns the inverse tangent of each element of the matrix.
Atan2<T>(Matrix<T>, Matrix<T>).	Returns the 4 quadrant inverse tangent of the ratio of corresponding elements in two matrices.
Sinh<T>(Matrix<T>).	Returns the hyperbolic sine of each element of the matrix.
Cosh<T>(Matrix<T>).	Returns the hyperbolic cosine of each element of the matrix.
Tanh<T>(Matrix<T>).	Returns the hyperbolic tangent of each element of the matrix.
Asinh<T>(Matrix<T>).	Returns the inverse hyperbolic sine of each element of the matrix.
Acosh<T>(Matrix<T>).	Returns the inverse hyperbolic cosine of each element of the matrix.
Atanh<T>(Matrix<T>).	Returns the inverse hyperbolic tangent of each element of the matrix.

References

G. H. Golub, C. F. Van Loan, Matrix Computations (3rd Ed), Johns Hopkins University Press, 1996.