Armadillo: C++ linear algebra library

Preamble

If something needs to be clarified or fixed in this documentation, the developers are interested in hearing about it.

To aid the conversion of Matlab/Octave programs, there is a syntax conversion table.

First time users may want to have a look at a short example program.

There is a small but important API change from earlier versions of Armadillo (< 0.9.x). Previously, some multiplication operations directly converted result matrices with a size of 1x1 into scalars. This is no longer the case. If you know the result of an expression will be a 1x1 matrix and wish to treat it as a scalar, use the as_scalar() wrapping function.

Matrix, Vector, Cube and Field Classes

Mat<type>, mat and cx_mat

Col<type>, colvec and vec

Member Functions

save/load (matrices & cubes)

Other Classes

running_stat_vec<type>

wall_clock

Generated Vectors/Matrices/Cubes

Functions Individually Applied to Each Element of a Matrix/Cube

abs

eps

misc functions (exp, log, log10, pow, sqrt, square)

trigonometric functions (cos, sin, ...)

Scalar Valued Functions of Vectors/Matrices/Cubes

Scalar/Vector Valued Functions of Vectors/Matrices

statistics (mean, stddev, ...)

Vector/Matrix/Cube Valued Functions of Vectors/Matrices/Cubes

Decompositions, Inverse and Related Functions

Miscellaneous

math constants (pi, ...)

physical constants (speed of light, ...)

log_add

s32/u32

cx_float/cx_double

Matlab/Armadillo syntax differences

example program

Matrix, Vector, Cube and Field Classes

Mat<type>
mat
cx_mat

The root template matrix class is Mat<type>, where type can be one of: char, int, float, double, std::complex<double>, etc.

For convenience the following typedefs have been defined:

imat	=	Mat<s32>
umat	=	Mat<u32>
fmat	=	Mat<float>
mat	=	Mat<double>
cx_fmat	=	Mat<cx_float>
cx_mat	=	Mat<cx_double>

Data is stored with column-major ordering (i.e. column by column).

Functions which are wrappers for LAPACK or ATLAS functions (generally matrix decompositions) are only valid for the following types: fmat, mat, cx_fmat, cx_mat.

In this documentation the mat type is used for convenience. It is possible to use other types instead, e.g. fmat.

Constructors:
- mat()
- mat(n_rows, n_cols)
- mat(mat)
- mat(rowvec)
- mat(colvec)
- cx_mat(mat,mat) (for constructing a complex matrix out of two real matrices)

Advanced constructors:
- mat(aux_mem*, n_rows, n_cols, copy_aux_mem = true)
  
  Create a matrix using data from writeable auxiliary memory. By default the matrix allocates its own memory and copies data from the auxiliary memory (for safety). However, if copy_aux_mem is set to false, the matrix will instead directly use the auxiliary memory (i.e. no copying). This is faster, but can be dangerous unless you know what you're doing!
- mat(const aux_mem*, n_rows, n_cols)
  
  Create a matrix by copying data from read-only auxiliary memory.

Examples:

mat A = rand<mat>(5,5);
double x = A(1,2);
mat B = A + A;
mat C = A * B;
mat D = A % B;
B.zeros();
B.set_size(10,10);
B.zeros(10,10);

cx_mat X(A,B);

Caveat: For mathematical correctness, scalars are treated as 1x1 matrices during initialisation. As such, the code below will not generate a 5x5 matrix with every element equal to 123.0:
Use the following code instead:
Or:

Col<type>
colvec
vec

Classes for column vectors (matrices with one column).

The Col<type> class is derived from the Mat<type> class and inherits most of the member functions.

For convenience the following typedefs have been defined:

ivec, icolvec	=	Col<s32>
uvec, ucolvec	=	Col<u32>
fvec, fcolvec	=	Col<float>
vec, colvec	=	Col<double>
cx_fvec, cx_fcolvec	=	Col<cx_float>
cx_vec, cx_colvec	=	Col<cx_double>

In this documentation the colvec type is used for convenience and to accentuate that it is a column vector. It is possible to use the shorter version, i.e. vec, as well as other types, e.g. fvec, fcolvec.

Functions which take Mat as input can generally also take Col as input. Main exceptions are functions which require square matrices.

Constructors
- colvec(n_elem=0)
- colvec(colvec)
- colvec(mat)

An exception is thrown if the given matrix has more than one column.

Examples:

colvec x(10);
colvec y = zeros<colvec>(10,1);

mat    A = rand<mat>(10,10);
colvec z = A.col(5); // extract a column vector

Caveat: For mathematical correctness, scalars are treated as 1x1 matrices during initialisation. As such, the code below will not generate a column vector with every element equal to 123.0:
Use the following code instead:
Or:

Row<type>
rowvec

Classes for row vectors (matrices with one row)

The template Row<type> class is derived from the Mat<type> class and inherits most of the member functions.

For convenience the following typedefs have been defined:

irowvec	=	Row<s32>
urowvec	=	Row<u32>
frowvec	=	Row<float>
rowvec	=	Row<double>
cx_frowvec	=	Row<cx_float>
cx_rowvec	=	Row<cx_double>

In this documentation the rowvec type is used for convenience. It is possible to use other types instead, e.g. frowvec.

Functions which take Mat as input can generally also take Row as input. Main exceptions are functions which require square matrices.

Constructors

An exception is thrown if the given matrix has more than one row.

Examples:

rowvec x(10);
rowvec y = zeros<mat>(1,10);

mat A = rand<mat>(10,10);
rowvec z = A.row(5); // extract a row vector

Caveat: For mathematical correctness, scalars are treated as 1x1 matrices during initialisation. As such, the code below will not generate a row vector with every element equal to 123.0:
Use the following code instead:
Or:

Cube<type>
cube
cx_cube

Classes for cubes (aka "3D matrices" or sets of matrices with contiguous memory)

The root template cube class is Cube<type>, where type can be one of: char, int, float, double, std::complex<double>, etc.

For convenience the following typedefs have been defined:

icube	=	Cube<s32>
ucube	=	Cube<u32>
fcube	=	Cube<float>
cube	=	Cube<double>
cx_fcube	=	Cube<cx_float>
cx_cube	=	Cube<cx_double>

In this documentation the cube type is used for convenience. It is possible to use other types instead, e.g. fcube.

Functions which take Mat as input can generally also take cube slices as input.

Constructors:

Advanced constructors:
- cube(aux_mem*, n_rows, n_cols, n_slices, copy_aux_mem = true)
  
  Create a cube using data from writeable auxiliary memory. By default the cube allocates its own memory and copies data from the auxiliary memory (for safety). However, if copy_aux_mem is set to false, the cube will instead directly use the auxiliary memory (i.e. no copying). This is faster, but can be dangerous unless you know what you're doing!
- cube(const aux_mem*, n_rows, n_cols, n_slices)
  
  Create a cube by copying data from read-only auxiliary memory.

Examples:

cube x(1,2,3);
cube y = rand<cube>(4,5,6);

mat A = y.slice(1);  // extract a slice from the cube
                     // (each slice is a matrix)

mat B = rand<mat>(4,5);
y.slice(2) = B;     // set a slice in the cube

cube q = y + y;     // cube addition
cube r = y % y;     // element-wise cube multiplication

Caveats
- The size of individual slices can't be changed. For example, the following will not work:
- For mathematical correctness, scalars are treated as 1x1x1 cubes during initialisation. As such, the code below will not generate a cube with every element equal to 123.0:
  Use the following code instead:
  Or:

field<object type>

Class for one and two dimensional fields of arbitrary objects

Constructors (where object type is another class, e.g. std::string, mat, colvec, rowvec, etc):

Examples:

// create a field of strings
field<std::string> S(3,2);

S(0,0) = "hello";
S(1,0) = "there";

// string fields can be saved as plain text files
S.save("string_field");

// create a colvec field with 3 rows and 2 columns
field<colvec> F(3,2);

// access components of the field
F(0,0) = colvec(5);
F(1,1) = rand<colvec>(6);
F(2,0).set_size(7);

// access element 1 of the colvec stored at 2,0
double x = F(2,0)(1);

// copy rows
F.row(0) = F.row(2);

// extract a row of colvecs from F
field<colvec> G = F.row(1);

// print the field to the standard output
G.print("G =");

// save the field to a binary file
G.save("colvec_field");

See also:
- subfields
- saving/loading fields
- string class in the standard C++ library (cplusplus.com)

Member Functions

.copy_size(A)

Member function of Mat, Col, Row, Cube and field.

Set the size to be the same as object A.

Object A must be of the same root type as the object being modified (e.g. you can't set the size of a matrix by providing a cube).

Examples:

mat A = rand<mat>(5,6);
mat B;
B.copy_size(A);

cout << B.n_rows << endl;
cout << B.n_cols << endl;

.diag(k=0)

Member function of Mat.

Read/write access to the k-th diagonal in a matrix.

The argument k is optional -- by default the main diagonal is accessed (k=0).

For k > 0, the k-th super-diagonal is accessed (top-right corner). For k < 0, the k-th sub-diagonal is accessed (bottom-left corner).

An extracted a diagonal is interpreted as a column vector.

Examples:

mat    X = rand<mat>(5,5);

colvec a = X.diag();
colvec b = X.diag(1);
colvec c = X.diag(-2);

X.diag() = rand<colvec>(5);

See also: diagvec()

element/object access via (), [] and .at()

Provide access to individual elements or objects stored in a container object (i.e., Mat, Col, Row, Cube, field).

(n)		For colvec and rowvec, access the n-th element. For mat, cube and field, access the n-th element/object under the assumption of a flat layout, with column-major ordering of data (i.e. column by column). An exception is thrown if the requested element is out of bounds. The bounds check can be optionally disabled if the ARMA_NO_DEBUG or NDEBUG macro is defined.

[n]		As for (n), but without a bounds check. Not recommended for use unless your code has been thoroughly debugged.

(i,j)		For mat and field classes, access the element/object stored at the i-th row and j-th column. An exception is thrown if the requested element is out of bounds. The bounds check can be optionally disabled if the ARMA_NO_DEBUG or NDEBUG macro is defined.

.at(i,j)		As for (i,j), but without a bounds check. Not recommended for use unless your code has been thoroughly debugged.

(i,j,k)		Cube only: access the element stored at the i-th row, j-th column and k-th slice. An exception is thrown if the requested element is out of bounds. The bounds check can be optionally disabled if the ARMA_NO_DEBUG or NDEBUG macro is defined.

.at(i,j,k)		As for (i,j,k), but without a bounds check. Not recommended for use unless your code has been thoroughly debugged.

Examples:

mat A = rand<mat>(10,10);
A(9,9) = 123.0;
double x = A.at(9,9);
double y = A[99];

colvec p = rand<colvec>(10,1);
p(9) = 123.0;
double z = p[9];

.fill(value)

Member function of Mat, Col, Row and Cube classes.

Sets the elements to a specified value. The type of value must match the type of elements used by the container object (e.g. for mat the type is double).

Examples:

See also:
- .ones()
- .zeros()

.is_finite()

Member function of Mat, Col, Row and Cube classes.

Returns true if all elements of the object are finite.

Returns false if one or more elements of the object is non-finite.

Examples:

mat A = rand<mat>(5,5);
mat B = rand<mat>(5,5);

B(1,1) = 0.0 / 0.0;  // create a NaN

cout << A.is_finite() << endl;
cout << B.is_finite() << endl;

iterators (matrices & vectors)

STL-style iterators and associated member functions of the Mat, Col and Row classes.

iterator types:

mat::iterator colvec::iterator rowvec::iterator		random access iterators, for read/write access to elements; the elements are ordered column by column

mat::const_iterator colvec::const_iterator rowvec::const_iterator		random access iterators, for read-only access to elements; the elements are ordered column by column

mat::col_iterator colvec::col_iterator rowvec::col_iterator		random access iterators, for read/write access to the elements of a specific column

mat::const_col_iterator colvec::const_col_iterator rowvec::const_col_iterator		random access iterators, for read-only access to the elements of a specific column

mat::row_iterator		rudimentary forward iterator, for read/write access to the elements of a specific row

mat::const_row_iterator		rudimentary forward iterator, for read-only access to the elements of a specific row

colvec::row_iterator rowvec::row_iterator		random access iterators, for read/write access to the elements of a specific row

colvec::const_row_iterator rowvec::const_row_iterator		random access iterators, for read-only access to the elements of a specific row

Member functions:

.begin()		provide an iterator referring to the first element
.end()		provide an iterator referring to the past-the-end element
.begin_row(row_number)		provide an iterator referring to the first element of the specified row
.end_row(row_number)		provide an iterator referring to the past-the-end element of the specified row
.begin_col(column_number)		provide an iterator referring to the first element of the specified column
.end_col(column_number)		provide an iterator referring to the past-the-end element of the specified column

Examples:

mat X = rand<mat>(5,5);


mat::iterator a = X.begin();
mat::iterator b = X.end();

for(mat::iterator i=a; i!=b; ++i)
  {
  cout << *i << endl;
  }


mat::col_iterator c = X.begin_col(1);  // start of column 1
mat::col_iterator d = X.end_col(3);    // end of column 3

for(mat::col_iterator i=c; i!=d; ++i)
  {
  cout << *i << endl;
  (*i) = 123.0;
  }

See also: iterator at cplusplus.com

iterators (cubes)

STL-style iterators and associated member functions of the Cube class.

iterator types:

cube::iterator		random access iterator, for read/write access to elements; the elements are ordered slice by slice; the elements within each slice are ordered column by column

cube::const_iterator		random access iterators, for read-only access to elements

cube::slice_iterator		random access iterator, for read/write access to the elements of a particular slice; the elements are ordered column by column

cube::const_slice_iterator		random access iterators, for read-only access to the elements of a particular slice

Member functions:

.begin()		provide an iterator referring to the first element
.end()		provide an iterator referring to the past-the-end element
.begin_slice(slice_number)		provide an iterator referring to the first element of the specified slice
.end_slice(slice_number)		provide an iterator referring to the past-the-end element of the specified slice

Examples:

cube X = rand<cube>(2,3,4);


cube::iterator a = X.begin();
cube::iterator b = X.end();

for(cube::iterator i=a; i!=b; ++i)
  {
  cout << *i << endl;
  }


cube::slice_iterator c = X.begin_slice(1);  // start of slice 1
cube::slice_iterator d = X.end_slice(2);    // end of slice 2

for(cube::slice_iterator i=c; i!=d; ++i)
  {
  cout << *i << endl;
  (*i) = 123.0;
  }

See also: iterator at cplusplus.com

.n_rows		(number of rows)
.n_cols		(number of columns)
.n_elem		(number of elements)
.n_slices		(number of slices)
.n_elem_slice		(number of elements per slice)

n_rows, n_cols and n_elem are read-only member variables of Mat, Col, Row, Cube and field classes.

n_slices and n_elem_slice are read-only member variables of the Cube class.

The variables are of type u32 (unsigned integer).

Examples:

colvec q(5);
cout << "q has " << q.n_cols << " columns" << endl;

.ones()
.ones(n_elem)
.ones(n_rows, n_cols)
.ones(n_rows, n_cols, n_slices)

Set the elements of an object to one, optionally first resizing to specified dimensions

.ones() and .ones(n_elem) are member function of Col and Row

.ones() and .ones(n_rows, n_cols) are member functions of Mat

.ones() and .ones(n_rows, n_cols, n_slices) are member functions of Cube

Examples:

mat A = rand<mat>(5,10);
A.ones();      // sets all elements to one
A.ones(10,20); // sets the size to 10 cols and 20 rows
               // followed by setting all elements to one

See also:
- ones() (stand alone function)
- .zeros()
- .fill()

operators + - * / % == != <= >= < >

Overloaded operators for mat, colvec, rowvec and cube classes.

Meanings:

+		Addition of two objects.
-		Subtraction of one object from another or negation of an object.
/		Element-wise division of an object by another object or a scalar.
*		Matrix multiplication of two objects. Not applicable to the cube class unless multiplying a cube by a scalar.
%		Schur product: element-wise multiplication of two objects.
==		Element-wise equality evaluation of two objects. Generates a matrix of type umat with entries that indicate whether at a given position the two elements from the two objects are equal (1) or not equal (0).
!=		Element-wise non-equality evaluation of two objects.
>=		As for ==, but the check is for "greater than or equal to".
<=		As for ==, but the check is for "less than or equal to".
>		As for ==, but the check is for "greater than".
<		As for ==, but the check is for "less than".

An exception is thrown if incompatible object sizes are used.

If the +, - and % operators are chained, Armadillo will try to avoid the generation of temporaries. No temporaries are generated if all given objects are of the same type and size.

If the * operator is chained, Armadillo will try to find an efficient ordering of the matrix multiplications.

Caveat: operators involving an equality comparison (i.e., ==, !=, >=, <=) may not work as expected for floating point element types (i.e., float, double) due to the necessarily limited precision of these types. In other words, these operators are not recommended for matrices of type mat or fmat.

Examples:

mat A = rand<mat>(5,10);
mat B = rand<mat>(5,10);
mat C = rand<mat>(10,5);

mat P = A + B;
mat Q = A - B;
mat R = -B;
mat S = A / 123.0;
mat T = A % B;
mat U = A * C;

// V is constructed without temporaries
mat V = A + B + A + B;

imat AA = "1 2 3; 4 5 6; 7 8 9;";
imat BB = "3 2 1; 6 5 4; 9 8 7;";

// compare elements
umat ZZ = (AA >= BB);

.print(header="")
.print(stream, header="")

Member function of Mat, Col, Row, Cube and field.

The first form prints the contents of an object to the cout stream, with an optional header line. The second form prints to a user specified stream.

It's also possible print objects using the << stream operator.

Elements of a field can only be printed if there is an associated operator<< function defined.

Examples:

mat A = rand<mat>(5,5);
mat B = rand<mat>(6,6);

A.print();

// "B =" is the optional header line
B.print("B =");

cout << A << endl;
cout << "B =" << endl << B << endl;

.print_trans(header="")
.print_trans(stream, header="")

Member function of Mat, Col and Row.

Similar to .print(), with the difference being that a transposed version of the object is printed.

.raw_print(header="")
.raw_print(stream, header="")

Member function of Mat, Col, Row and Cube.

Similar to the .print() member function. The difference is that no formatting of the output is done -- i.e. the user can set the stream's parameters such as precision, cell width, etc.

If the cell width is set to zero, a space is printed between the elements.

Examples:

mat A = rand<mat>(5,5);

cout.precision(11);
cout.setf(ios::fixed);

A.raw_print(cout, "A =");

.raw_print_trans(header="")
.raw_print_trans(stream, header="")

Member function of Mat, Col and Row.

Similar to .raw_print(), with the difference being that a transposed version of the object is printed.

.reset()

Member function of Mat, Col, Row, Cube and field.

Causes an object to have no elements.

Examples:

See also: .set_size()

.save(name, file_type = arma_binary)
.save(stream, file_type = arma_binary)

.load(name, file_type = auto_detect)
.load(stream, file_type = auto_detect)

Store/retrieve data in files or streams.

Member functions of Mat, Col, Row and Cube classes.

load() will reset the object to have no elements if the loading process failed (e.g. wrong filename, unsupported file type, etc).

The following file formats are supported:

auto_detect		load(): try to automatically detect the file type as one of the formats described below. This is the default operation.
raw_ascii		Numerical data stored in raw ASCII format (no header), with the numbers separated by whitespace. The number of columns must be the same in each row. load() can read data which was saved in Matlab/Octave using the -ascii option, except for complex numbers. Complex numbers are stored in standard C++ notation (a tuple surrounded by brackets: e.g. (1.23,4.56) indicates 1.24 + 4.56i). Cubes are loaded as one slice.
arma_ascii		Numerical data stored in human readable text format, with a simple header to speed up loading. The header indicates the type of matrix as well as the number of rows and columns. For cubes, the header additionally specifies the number of slices.
arma_binary		Numerical data stored in machine dependent binary format, with a simple header to speed up loading. The header indicates the type of matrix as well as the number of rows and columns. For cubes, the header additionally specifies the number of slices.
pgm_binary		Image data stored in Portable Gray Map (PGM) format. Applicable to Mat only. Saving int, float or double matrices is a lossy operation, as each element is copied and converted to an 8 bit representation. As such the matrix should have values in the [0,255] interval, otherwise the resulting image may not display correctly.
ppm_binary		Image data stored in Portable Pixel Map (PPM) format. Applicable to Cube only. Saving int, float or double matrices is a lossy operation, as each element is copied and converted to an 8 bit representation. As such the cube/field should have values in the [0,255] interval, otherwise the resulting image may not display correctly.

Examples:

mat A = rand<mat>(5,5);

A.save("A1.mat");  // default save format is arma_binary
A.save("A2.mat", arma_ascii);

mat B;
// automatically detect format type
B.load("A1.mat");

mat C;
// force loading in the arma_ascii format
C.load("A2.mat", arma_ascii);


// example of saving/loading using a stream

std::stringstream s;
A.save(s);

mat D;
D.load(s);

See also: saving/loading fields

.save(name, file_type = arma_binary)
.save(stream, file_type = arma_binary)

.load(name, file_type = auto_detect)
.load(stream, file_type = auto_detect)

Store/retrieve fields in files or stream.

Member functions of the field class.

load() will reset the object to have no elements if the loading process failed (e.g. wrong filename, unsupported file type, etc).

Fields with objects of type std::string are saved and loaded as raw text files. The text files do not have a header. Each string is separated by a whitespace. The .load() function will only accept text files that have the same number of strings on each line. The strings can have variable lengths.

Other than storing string fields as text files, the following file formats are supported:

auto_detect		load(): try to automatically detect the field format type as one of the formats described below. This is the default operation.
arma_binary		Objects are stored in machine dependent binary format. Default type for fields of type Mat, Col or Row. Only applicable to fields of type Mat, Col or Row.
ppm_binary		Image data stored in Portable Pixmap Map (PPM) format. Only applicable to fields of type Mat, Col or Row. .load(): Loads the specified image and stores the red, green and blue components as three separate matrices. The resulting field is comprised of the three matrices, with the red, green and blue components in the first, second and third matrix, respectively. .save(): Saves a field with exactly three matrices of equal size as an image. It is assumed that the red, green and blue components are stored in the first, second and third matrix, respectively. Saving int, float or double matrices is a lossy operation, as each matrix element is copied and converted to an 8 bit representation.

See also: saving/loading matrices and cubes

.set_size(n_elem)		(member function of Col, Row, and field)
.set_size(n_rows, n_cols)		(member function of Mat and field)
.set_size(n_rows, n_cols, n_slices)		(member function of Cube)

Changes the size of an object.

If the old number of elements is equal to the requested number of elements, old memory is reused.

Examples:

mat A;
A.set_size(5,10);

colvec q;
q.set_size(100);

See also: .reset()

submatrix views

A collection of member functions of Mat, Col and Row classes that provide submatrix views.

Examples:

mat A = zeros<mat>(5,10);
mat B = A.submat(0,1,2,3);
A.submat(0,1,2,3) = rand<mat>(3,3);
A.row(1) = rand<mat>(1,10);

See also: join_rows() & join_cols()

subcube views and slices

A collection of member functions of the Cube class that provide subcube views.

An individual slice can be interpreted as a matrix (i.e. an instance of the Mat class).

The size of individual slices can't be changed.

Examples:

cube A = rand<cube>(2,3,4);
mat B = A.slice(1);
A.slice(0) = rand<mat>(2,3);
A.subcube(0,0,1,  1,1,2) = rand<cube>(2,2,2);
A.slice(0)(1,2) = 99.0;

subfield views

A collection of member functions of the field class that provide subfield views.

.swap_rows(row1, row2)
.swap_cols(col1, col2)

Member functions of Mat, Col and Row classes.

Swap the contents of specified rows or columns.

Examples:

mat X = rand<mat>(5,5);
X.swap_rows(0,4);

See also: fliplr() & flipud()

.zeros()
.zeros(n_elem)
.zeros(n_rows, n_cols)
.zeros(n_rows, n_cols, n_slices)

Set the elements of an object to zero, optionally first resizing to specified dimensions

.zeros() and .zeros(n_elem) are member function of Col and Row

.zeros() and .zeros(n_rows, n_cols) are member functions of Mat

.zeros() and .zeros(n_rows, n_cols, n_slices) are member functions of Cube

Examples:

mat A = rand<mat>(5,10);
A.zeros();      // sets all elements to zero
A.zeros(10,20); // sets the size to 10 columns and 20 rows
                // followed by setting all elements to zero

See also:
- zeros() (stand alone function)
- .ones()
- .fill()
- .reset()

Other Classes

running_stat<type>

Class for keeping statistics of a continuously sampled one dimensional process/signal. Useful if the storage of individual samples (scalars) is not necessary or desired. Also useful if the number of samples is not known beforehand or exceeds available memory.

type should be one of: float, double, cx_float, cx_double

Member functions:

.operator()(scalar)		update the statistics so far using the given scalar
.mean()		get the mean or average value so far
.var(norm_type=0)		get the variance so far
.stddev(norm_type=0)		get the standard deviation so far
.min()		get the minimum value so far
.max()		get the maximum value so far
.reset()		reset all statistics and set the number of samples to zero

For the .var() and .stddev() functions, the default norm_type=0 performs normalisation using N-1 (where N is the number of samples so far), providing the best unbiased estimator. Using norm_type=1 causes normalisation to be done using N, which provides the second moment around the mean.

Examples:

running_stat<double> stats;

for(u32 i=0; i<10000; ++i)
  {
  double sample = double(rand())/RAND_MAX;
  stats(sample);
  }

cout << "mean = " << stats.mean() << endl;
cout << "var  = " << stats.var()  << endl;
cout << "min  = " << stats.min()  << endl;
cout << "max  = " << stats.max()  << endl;

See also:
- statistics functions
- running_stat_vec

running_stat_vec<type>(calc_cov = false)

Class for keeping statistics of a continuously sampled multi-dimensional process/signal. Useful if the storage of individual samples is not necessary or desired. Also useful if the number of samples is not known beforehand or exceeds available memory.

This class is similar to running_stat, with the difference being that vectors are processed instead of single values.

type must match the element type used by the sample vectors (i.e. it must be one of float, double, cx_float, cx_double)

type can be inferred via the use of the elem_type member typedef of a vector class. For example, fcolvec::elem_type will be interpreted as float, while cx_colvec::elem_type will be interpreted as cx_double.

Member functions:

.operator()(vector)		update the statistics so far using the given vector
.mean()		get the mean vector so far
.var(norm_type=0)		get the vector of variances so far
.stddev(norm_type=0)		get the vector of standard deviations so far
.min()		get the vector of minimum values so far
.max()		get the vector of maximum values so far
.cov(norm_type=0)		get the covariance matrix so far NOTE: this only works if calc_cov is set to true during the construction of the class
.reset()		reset all statistics and set the number of samples to zero

For the .var() and .stddev() functions, the default norm_type=0 performs normalisation using N-1 (where N is the number of samples so far), providing the best unbiased estimator. Using norm_type=1 causes normalisation to be done using N, which provides the second moment around the mean.

Examples:

running_stat_vec<rowvec::elem_type> stats;

rowvec sample;

for(u32 i=0; i<10000; ++i)
  {
  sample = rand<rowvec>(5);
  stats(sample);
  }

cout << "mean = " << stats.mean() << endl;
cout << "var  = " << stats.var()  << endl;
cout << "min  = " << stats.min()  << endl;
cout << "max  = " << stats.max()  << endl;

//
//

running_stat_vec<rowvec::elem_type> more_stats(true);

for(u32 i=0; i<20; ++i)
  {
  sample = rand<rowvec>(3);
  
  sample(1) -= sample(0);
  sample(2) += sample(1);
  
  more_stats(sample);
  }

cout << "covariance matrix = " << endl;
cout << more_stats.cov() << endl;

rowvec sd = more_stats.stddev();

cout << "correlations = " << endl;
cout << more_stats.cov() / (trans(sd) * sd);

See also:

wall_clock

Simple wall clock timer class, for measuring the number of elapsed seconds between two intervals.

Examples:

wall_clock timer;

mat A = rand<mat>(4,4);
mat B = rand<mat>(4,4);
mat C;

timer.tic();
for(u32 i=0; i<100000; ++i)
  C = A + B + A + B;

double n_secs = timer.toc();
cout << "took " << n_secs << " seconds" << endl;

Generated Vectors/Matrices

eye(n_rows, n_cols)

Generate the identity matrix (diagonal elements set to one and off-diagonal elements set to zero).

Usage:
- matrix_type X = eye<matrix_type>(n_rows, n_cols)
- n_rows must be equal to n_cols
- Two arguments are used to make it explicitly clear that a square matrix is generated

Examples:

mat A = eye<mat>(5,5);
mat B = 123.0 * eye<mat>(5,5);

See also:

linspace(start, end, N)

Generate a vector with N elements; the values of the elements linearly increase from start upto (and including) end.

Usage:
- vector_type v = linspace<vector_type>(start, end, N)
- matrix_type X = linspace<matrix_type>(start, end, N)

If a matrix_type is specified, the resultant matrix will have one column.

Examples:

colvec q = linspace<colvec>(10, 20, 5);
rowvec r = linspace<rowvec>(10, 20, 5);

mat X = linspace<mat>(10, 20, 5);

ones(n_elem)
ones(n_rows, n_cols)
ones(n_rows, n_cols, n_slices)

Generate a vector, matrix or cube with all elements set to one.

Usage:
- vector_type v = ones<vector_type>(n_elem)
- matrix_type X = ones<matrix_type>(n_rows, n_cols)
- cube_type X = ones<cube_type>(n_rows, n_cols, n_slices)

Examples:

mat A = ones<mat>(5,6);
mat B = 123.0 * ones<mat>(5,6);
cube C = 123.0 * ones<cube>(5,6,7);

See also:
- .ones() (member function of Mat, Col, Row and Cube)
- eye()

rand(n_elem)
rand(n_rows, n_cols)
rand(n_rows, n_cols, n_slices)

randn(n_elem)
randn(n_rows, n_cols)
randn(n_rows, n_cols, n_slices)

Generate a vector, matrix or cube with the elements set to random values.

rand() uses a uniform distribution in the [0,1] interval.

randn() uses a normal/Gaussian distribution with zero mean and unit variance.

Usage:
- vector_type v = rand<vector_type>(n_elem)
- matrix_type X = rand<matrix_type>(n_rows, n_cols)
- cube_type C = rand<cube_type>(n_rows, n_cols, n_slices)

Examples:

colvec q = rand<colvec>(5);
rowvec r = rand<rowvec>(5);
mat A = rand<mat>(5,6);
cube C = rand<cube>(5,6,7);

repmat(mat A, n_copies_per_row, n_copies_per_col)

Generate a matrix by replicating matrix A in a block-like fashion.
The generated matrix has a size of (n_copies_per_row * A.n_rows, n_copies_per_col * A.n_cols)

Examples:

mat A = rand<mat>(2, 3);
mat B = repmat(A, 4, 5);

zeros(n_elem)
zeros(n_rows, n_cols)
zeros(n_rows, n_cols, n_slices)

Generate a vector, matrix or cube with the elements set to zero.

Usage:
- vector_type v = zeros<vector_type>(n_elem)
- matrix_type X = zeros<matrix_type>(n_rows, n_cols)
- cube_type X = zeros<cube_type>(n_rows, n_cols, n_slices)

Examples:

colvec q = zeros<colvec>(5);
rowvec r = zeros<rowvec>(5);
mat A = zeros<mat>(5,6);
cube C = zeros<cube>(5,6,7);

See also:
- .zeros() (member function of Mat, Col, Row and Cube)
- .ones() (member function of Mat, Col, Row and Cube)
- ones()

Functions Individually Applied to Each Element of a Matrix/Cube

abs(mat)
abs(cube)
abs(cx_mat)
abs(cx_cube)

Obtain the magnitude of each element.

Usage for non-complex matrices:
- matrix_type Y = abs(X)
- X and Y must have the same matrix_type

Usage for non-complex cubes:
- cube_type Y = abs(X)
- X and Y must have the same cube_type

Usage for complex matrices:
- non_complex_matrix_type Y = abs(X)
- X must be a have complex matrix type, e.g., cx_mat or cx_fmat
- The type of Y must be related to the type of X, e.g., if X has the type cx_mat, then the type of Y must be mat

Usage for complex cubes:
- non_complex_cube_type Y = abs(X)
- X must be a have complex cube type, e.g., cx_cube or cx_fcube
- The type of Y must be related to the type of X, e.g., if X has the type cx_cube, then the type of Y must be cube

Examples:

mat A = rand<mat>(5,5);
mat B = abs(A); 

cx_mat X = rand<cx_mat>(5,5);
mat    Y = abs(X);

eps(mat)
eps(cx_mat)
eps(float)
eps(double)
eps(cx_float)
eps(cx_double)

Obtain the positive distance of the absolute value of each element to the next largest representable floating point number.

Examples:

mat A = rand<mat>(4,5);

mat X = eps(A);

See also:

misc functions (exp, log, log10, pow, sqrt, square)

Apply a function to each element.

Usage:

matrix_type Y = misc_fn(X)
cube_type Y = misc_fn(X)
X and Y must have the same matrix_type/cube_type

misc_fn(X) is one of:

exp(X)		base-e exponential of each element
log(X)		natural log of each element
log10(X)		base-10 log of each element
pow(X, p)		raise each element to the power of p
sqrt(X)		square root of each element
square(X)		square of each element

Examples:

mat X = rand<mat>(5,5);
mat Y = exp(X);

trigonometric functions (cos, sin, tan, ...)

Apply a trigonometric function to each element.

Usage:
- matrix_type Y = trig_fn(X)
- cube_type Y = trig_fn(X)
- X and Y must have the same matrix_type/cube_type
- trig_fn is one of:
  - cos family: cos, acos, cosh, acosh
  - sin family: sin, asin, sinh, asinh
  - tan family: tan, atan, tanh, atanh

Examples:

mat X = rand<mat>(5,5);
mat Y = cos(X);

Scalar Valued Functions of Vectors/Matrices/Cubes

accu(mat)
accu(cube)

Accumulate (sum) all elements

Examples:

mat A = rand<mat>(5,5);
double x = accu(A);

mat B = rand<mat>(5,5);
double y = accu(A % B);

// operator % performs element-wise multiplication,
// hence accu(A % B) is a "multiply-and-accumulate"
// operation

See also: sum()

as_scalar(expression)

Evaluate an expression that results in a 1x1 matrix, followed by converting the 1x1 matrix to a scalar.

If a binary or trinary expression is given (i.e. 2 or 3 terms), the function will try to exploit the fact that the result is a 1x1 matrix by using optimised expression evaluations.

Examples:

rowvec r = rand<rowvec>(5);
colvec q = rand<colvec>(5);
mat    X = rand<mat>(5,5);

// examples of some expressions
// for which optimised implementations exist

double a = as_scalar(r*q);
double b = as_scalar(r*X*q);
double c = as_scalar(r*diagmat(X)*q);
double d = as_scalar(r*inv(diagmat(X))*q);

See also:
- conv_to()
- dot()

det(mat)

Determinant of a square matrix.

An exception is thrown if the given matrix is not square.

Caveat: for large matrices you may want to use log_det() instead.

Examples:

mat A = rand<mat>(5,5);
double x = det(A);

See also:
- determinant in MathWorld
- log_det()

dot(A, B)
norm_dot(A, B)

dot(A,B): dot product of A and B, under the assumption that A and B are vectors with the same number of elements

norm_dot(A,B): normalised version of dot(A,B)

Examples:

colvec a = rand<colvec>(10);
colvec b = rand<colvec>(10);

double x = dot(a,b);

See also:
- as_scalar()
- cross()

log_det(val, sign, mat)

Log determinant of a square matrix, such that the determinant is equal to exp(val)*sign.

An exception is thrown if the given matrix is not square.

Examples:

mat A = rand<mat>(5,5);

double val;
double sign;

log_det(val, sign, A);

See also:
- determinant in MathWorld
- det()

norm(X, p)

Compute the p-norm of X, under the assumption that X is a vector.

p is an integer >= 1, or one of: "inf", "-inf".

For p = "inf" the maximum norm is computed, while for p = "-inf", the minimum norm is calculated.

Examples:

colvec q = rand<colvec>(5);
double x = norm(q, 2);
double y = norm(q, "inf");

See also:
- Norm in Wikipedia
- Vector Norm in MathWorld

rank(mat)
rank(cx_mat)

Returns the rank of a matrix.

The computation is based on singular value decomposition.

Any singular values less than default tolerance are treated as zero.

For a matrix A with m rows and n columns the default tolerance is max(m,n)*norm(A)*math::eps().

Examples:

mat A = rand<mat>(4,5);

u32 r = rank(A);

See also:
- math::eps()
- Matrix Rank in MathWorld

trace(mat)

Sum of the diagonal elements of a square matrix.

An exception is thrown if the given matrix is not square.

Examples:

mat A = rand<mat>(5,5);
double x = trace(A);

Scalar/Vector Valued Functions of Vectors/Matrices

diagvec(mat, k=0)

Extract the k-th diagonal from a matrix.

The argument k is optional -- by default the main diagonal is extracted (k=0).

For k > 0, the k-th super-diagonal is extracted (top-right corner). For k < 0, the k-th sub-diagonal is extracted (bottom-left corner).

An extracted a diagonal is interpreted as a column vector.

Examples:

mat A = rand<mat>(5,5);

colvec d = diagvec(A);

See also:
- .diag()
- diagmat()

min(mat, dim=0)
min(rowvec)
min(colvec)

max(mat, dim=0)
max(rowvec)
max(colvec)

For a matrix argument, return the minimum/maximum value for each column (dim=0), or each row (dim=1).

For a vector argument, return the minimum/maximum value.

Examples:

colvec q = rand<colvec>(10,1);
double x = max(q);

mat A = rand<mat>(10,10);
rowvec b = max(A);

// same result as max(A)
// the 0 explicitly indicates
// "traverse across rows"
rowvec c = max(A,0); 

// the 1 explicitly indicates
// "traverse across columns"
colvec d = max(A,1);

// find the overall maximum value
double y = max(max(A));

prod(mat, dim=0)
prod(rowvec)
prod(colvec)

For a matrix argument, return the product of elements in each column (dim=0), or each row (dim=1).

For a vector argument, return the product of all elements.

Examples:

colvec q = rand<colvec>(10,1);
double x = prod(q);

mat A = rand<mat>(10,10);
rowvec b = prod(A);

// same result as prod(A)
// the 0 explicitly indicates
// "traverse across rows"
rowvec c = prod(A,0);

// the 1 explicitly indicates
// "traverse across columns"
colvec d = prod(A,1);

// find the overall product
double y = prod(prod(A));

See also: Schur product

sum(mat, dim=0)
sum(rowvec)
sum(colvec)

For a matrix argument, return the sum of elements in each column (dim=0), or each row (dim=1).

For a vector argument, return the sum of all elements.

Examples:

colvec q = rand<colvec>(10,1);
double x = sum(q);

mat A = rand<mat>(10,10);
rowvec b = sum(A);

// same result as sum(A)
// the 0 explicitly indicates
// "traverse across rows"
rowvec c = sum(A,0);

// the 1 explicitly indicates
// "traverse across columns"
colvec d = sum(A,1);

// find the overall sum
double y = sum(sum(A));

See also: accu()

statistics: mean, median, stddev, var

mean(mat, dim=0) mean(colvec) mean(rowvec)		mean (average value)
median(mat, dim=0) median(colvec) median(rowvec)		median
stddev(mat, norm_type=0, dim=0) stddev(colvec, norm_type=0) stddev(rowvec, norm_type=0)		standard deviation
var(mat, norm_type=0, dim=0) var(colvec, norm_type=0) var(rowvec, norm_type=0)		variance

For a matrix argument, find a particular statistic for each column (dim=0), or each row (dim=1)

For a vector argument, return a particular statistic calculated using all the elements of the vector.

For the var() and stddev() functions, the default norm_type=0 performs normalisation using N-1 (where N is the number of samples), providing the best unbiased estimator. Using norm_type=1 causes normalisation to be done using N, which provides the second moment around the mean.

Examples:

mat A    = rand<mat>(5,5);
mat B    = mean(A);
mat C    = var(A);
double m = mean(mean(A));

colvec q = rand<colvec>(5);
double v = var(q);

See also:

Vector/Matrix/Cube Valued Functions of Vectors/Matrices/Cubes

conv_to<type>::from(X)

A form of casting

Convert between matrix types (e.g. mat to fmat) as well as cube types (e.g. cube to fcube).

Conversion of mat into colvec or rowvec is only possible if the object has one column or one row, respectively.

Examples:

mat  A = rand<mat>(5,5);

fmat B = conv_to<fmat>::from(A);

See also: as_scalar()

conj(cx_mat)
conj(cx_cube)

Obtain the complex conjugate of each element in a complex matrix/cube.

Examples:

cx_mat X = rand<cx_mat>(5,5);
cx_mat Y = conj(X);

See also: htrans()

cor(mat, mat, norm_type=0)
cor(rowvec, rowvec, norm_type=0)
cor(colvec, colvec, norm_type=0)
cor(colvec, rowvec, norm_type=0)
cor(rowvec, colvec, norm_type=0)

cor(mat, norm_type=0)
cor(rowvec, norm_type=0)
cor(colvec, norm_type=0)

For two matrix arguments X and Y, if each row of X and Y is an observation and each column is a variable, the (i,j)-th entry of cor(X,Y) is the correlation coefficient between the i-th variable in X and the j-th variable in Y.

For vector arguments, the type of vector is ignored and each element in the vector is treated as an observation.

For matrices, X and Y must have the same dimensions.

For vectors, X and Y must have the same number of elements.

cor(X) is equivalent to cor(X, X), also called autocorrelation.

The default norm_type=0 performs normalisation of the correlation matrix using N-1 (where N is the number of observations). Using norm_type=1 causes normalisation to be done using N.

Examples:

mat X = rand<mat>(4,5);
mat Y = rand<mat>(4,5);

mat R = cor(X,Y);

See also:
- Correlation in MathWorld
- Autocorrelation in MathWorld

cov(mat, mat, norm_type=0)
cov(rowvec, rowvec, norm_type=0)
cov(colvec, colvec, norm_type=0)
cov(colvec, rowvec, norm_type=0)
cov(rowvec, colvec, norm_type=0)

cov(mat, norm_type=0)
cov(rowvec, norm_type=0)
cov(colvec, norm_type=0)

For two matrix arguments X and Y, if each row of X and Y is an observation and each column is a variable, the (i,j)-th entry of cov(X,Y) is the covariance between the i-th variable in X and the j-th variable in Y.

For vector arguments, the type of vector is ignored and each element in the vector is treated as an observation.

For matrices, X and Y must have the same dimensions.

For vectors, X and Y must have the same number of elements.

cov(X) is equivalent to cov(X, X).

The default norm_type=0 performs normalisation using N-1 (where N is the number of observations), providing the best unbiased estimation of the covariance matrix (if the observations are from a normal distribution). Using norm_type=1 causes normalisation to be done using N, which provides the second moment matrix of the observations about their mean.

Examples:

mat X = rand<mat>(4,5);
mat Y = rand<mat>(4,5);

mat C = cov(X,Y);

See also:

cross(A, B)

Calculate the cross product between A and B, under the assumption that A and B are 3 dimensional vectors.

Examples:

colvec a = rand<colvec>(3);
colvec b = rand<colvec>(3);

colvec c = cross(a,b);

See also:

diagmat(mat)
diagmat(rowvec)
diagmat(colvec)

Interpret a matrix or vector as a diagonal matrix.

For mat, given matrix must be square. The main diagonal is copied and all other elements in the generated matrix are set to zero.

For colvec and rowvec, elements of the vector are placed on the main diagonal in the generated matrix and all other elements are set to zero.

Examples:

mat A = rand<mat>(5,5);
mat B = diagmat(A);
mat C = A*diagmat(A);

rowvec q = rand<rowvec>(5);
colvec r = rand<colvec>(5);
mat X = diagmat(q)*diagmat(r);

See also: diagvec()

find(X, k=0, s="first")

Return a column vector (with elements of type u32) of the indices of nonzero elements of X.

If X is a matrix, it is intepreted as a vector, with column-by-column ordering of the elements of X.

If k=0 (default), return the indices of all nonzero elements, otherwise return at most k of their indices.

If s="first" (default), return at most the first k indices of the nonzero elements.

If s="last", return at most the last k indices of the nonzero elements.

Relational operators can be used instead of X, though the current implementation is not yet optimal in terms of speed.

Examples:

mat     A  = rand<mat>(10,10);
ucolvec q1 = find(A > 0.5);
ucolvec q2 = find(A > 0.5, 10, "last");

See also: sort_index()

fliplr(mat)
flipud(mat)

fliplr: generate a copy of the input matrix, with the order of the columns reversed

flipud: generate a copy of the input matrix, with the order of the rows reversed

Examples:

mat A = rand<mat>(5,5);

mat B = fliplr(A);
mat C = flipud(A);

See also: .swap_rows() & .swap_cols() (member functions of Mat, Col and Row classes)

htrans(cx_mat)
htrans(cx_colvec)
htrans(cx_rowvec)

Hermitian transpose.

This operation is equivalent to taking the transpose and then taking the complex conjugate of each element.

Examples:

cx_mat A = rand<cx_mat>(5,10);
cx_mat B = htrans(A);

See also:
- trans()
- conj()

imag(cx_mat)
imag(cx_cube)

real(cx_mat)
real(cx_cube)

Extract the imaginary/real part of a complex matrix/cube.

Examples:

cx_mat A = rand<cx_mat>(5,5);
mat B = imag(A);
mat C = real(A);

join_rows(A,B)
join_cols(A,B)

join_rows: for two input matrices A and B, append each row of B to its respective row of A, provided that A and B have the same number of rows.

join_cols: for two input matrices A and B, append each column of B to its respective column of A, provided that A and B have the same number of columns.

Examples:

mat A = rand<mat>(4,5);
mat B = rand<mat>(4,6);
mat C = rand<mat>(6,5);

mat X = join_rows(A,B);
mat Y = join_cols(A,C);

See also: submatrix views

kron(mat,mat)
kron(cx_mat,cx_mat)
kron(mat,cx_mat)
kron(cx_mat,mat)

Kronecker tensor product.

Using matrix A (with n rows and p columns) and matrix B (with m rows and q columns), kron(A,B) returns a matrix (with nm rows and pq columns) which denotes the tensor product of A and B.

Examples:

mat A = rand<mat>(4,5);
mat B = rand<mat>(5,4);

mat K = kron(A,B);

See also:
- Kronecker Product in MathWorld

reshape(mat, n_rows, n_cols, dim=0)
reshape(cube, n_rows, n_cols, n_slices, dim=0)

Return a matrix/cube sized according to given size specifications, whose elements are taken from the given matrix/cube, either column-wise (dim=0) or row-wise (dim=1).

For matrices, the number of elements in the given matrix must be equal to n_rows * n_cols.

For cubes, the number of elements in the given cube must be equal to n_rows * n_cols * n_slices.

Examples:

mat A = rand<mat>(10, 5);
mat B = reshape(A, 5, 10);

shuffle(mat, dim=0)
shuffle(rowvec, dim=0)
shuffle(colvec, dim=0)

Shuffle the rows (dim=0) or columns (dim=1) of a matrix or vector.

Examples:

mat A = rand<mat>(4,5);
mat B = shuffle(A);

sort(mat, sort_type=0, dim=0)
sort(rowvec, sort_type=0)
sort(colvec, sort_type=0)

For a matrix argument, return a matrix with the elements of the input matrix sorted in each column (dim=0), or each row (dim=1)

sort_type=0 (default) indicates an ascending sort

sort_type=1 indicates a descending sort

For a vector argument, return a vector which is a sorted version of the input vector

Examples:

mat A = rand<mat>(10,10);
mat B = sort(A);

sort_index(colvec, sort_type=0)
sort_index(rowvec, sort_type=0)

Sort the given vector and return a vector (with elements of type u32) which describes the sorted order (i.e. it contains the indices of the given vector's elements).

sort_type=0 (default) indicates an ascending sort

sort_type=1 indicates a descending sort

Examples:

colvec  q       = rand<colvec>(10);
ucolvec indices = sort_index(q);

See also: find()

trans(mat)
trans(colvec)
trans(rowvec)

Matrix transpose.

Caveat: for complex matrices this function does not take the conjugate of the elements -- you may wish to use htrans() instead.

Examples:

mat A = rand<mat>(5,10);
mat B = trans(A);

See also: htrans()

Decompositions and Related Functions

chol(R, X)
R = chol(X)

Cholesky decomposition of X, such that trans(R)*R = X.

X must be a symmetric, positive-definite matrix.

Examples:

mat X = rand<mat>(5,5);
mat Y = trans(X)*X;

mat R = chol(Y);

See also: Cholesky decomposition in MathWorld

eig_sym(colvec eigval, mat X)
eig_sym(colvec eigval, cx_mat X)

colvec eigval = eig_sym(mat X)
colvec eigval = eig_sym(cx_mat X)

eig_sym(colvec eigval, mat eigvec, mat X)
eig_sym(colvec eigval, cx_mat eigvec, cx_mat X)

Eigen decomposition of symmetric/hermitian matrix X.

The eigenvalues and corresponding eigenvectors are stored in eigval and eigvec, respectively.

The eigenvalues are not guaranteed to be ordered.

An exception is thrown if X is not square.

There is currently no check for symmetry of X.

Examples:

mat A = rand<mat>(10,10);
mat B = trans(A)*A;  // generate a symmetric matrix

colvec eigval;
mat    eigvec;

eig_sym(eigval, eigvec, B);

See also:
- princomp_cov()
- eigen decomposition in MathWorld

eig_gen(cx_colvec eigval, cx_mat eigvec, mat X, side='r')
eig_gen(cx_colvec eigval, cx_mat eigvec, cx_mat X, side='r')

eig_gen(cx_colvec eigval, mat l_eigvec, mat r_eigvec, mat X)
eig_gen(cx_colvec eigval, cx_mat l_eigvec, cx_mat r_eigvec, cx_mat X)

Eigen decomposition of general (non-symmetric/non-hermitian) square matrix X.

The eigenvalues and corresponding eigenvectors are stored in eigval and eigvec, respectively.

For the first two forms, side='r' (default) specifies that right eigenvectors are computed, while side='l' specifies that left eigenvectors are computed.

For the last two forms, both left and right eigenvectors are computed.

The eigenvalues are not guaranteed to be ordered.

An exception is thrown if X is not square.

Examples:

mat A = rand<mat>(10,10);

cx_colvec eigval;
cx_mat    eigvec;

eig_gen(eigval, eigvec, A);

See also: eigen decomposition in MathWorld

inv(mat)

Inverse of a square matrix.

An exception is thrown if the given matrix is not square.

If the matrix appears to be singular, a warning message is printed and the output matrix has no elements.

Examples:

mat A = rand<mat>(5,5);
mat B = inv(A);

// Diagonal elements in C are set to the
// reciprocal of the corresponding elements in A.
// Off-diagonal elements in C are set to zero.
mat C = inv(diagmat(A));

See also:

lu(mat L, mat U, mat X)
lu(mat L, mat U, mat P, mat X)

Lower-upper decomposition (with partial pivoting) of matrix X.

L is a lower-triangular matrix, U is an upper-triangular matrix and P is an optional permutation matrix.

The second form provides the permutation matrix P, such that P*X = L*U and X = trans(P)*L*U.

The second form is currently considerably slower than the first form.

Examples:

mat A = rand<mat>(10,10);

mat L;
mat U;
mat P;

lu(L, U, P, A);

mat B = trans(P)*L*U;

See also: LU decomposition in Wikipedia

pinv(mat A, tol = 0)
pinv(cx_mat A, tol = 0)

Moore-Penrose pseudo-inverse of matrix A.

The computation is based on singular value decomposition (SVD). If the SVD fails, an error message is printed and the output matrix has no elements.

Any singular values less than tol are treated as zero.

For matrix A with m rows and n columns, the default tolerance is max(m,n)*norm(A)*math::eps(), where math::eps() denotes the difference between 1 and the least value greater than 1 that is representable.

Examples:

mat A = rand<mat>(4,5);

mat B = pinv(A);

See also:

mat coeff = princomp(mat X)
cx_mat coeff = princomp(cx_mat X)

princomp(mat coeff, mat score, mat X)
princomp(cx_mat coeff, cx_mat score, cx_mat X)

princomp(mat coeff, mat score, colvec latent, mat X)
princomp(cx_mat coeff, cx_mat score, colvec latent, cx_mat X)

princomp(mat coeff, mat score, colvec latent, colvec tsquared, mat X)
princomp(cx_mat coeff, cx_mat score, colvec latent, cx_colvec tsquared, cx_mat X)

Principal component analysis of matrix X.

Each row of X is an observation and each column is a variable.

coeff: principal component coefficients.

score: projected data.

latent: eigenvalues of the covariance matrix of X.

tsquared: Hotteling's statistic for each sample.

The computation is based on singular value decomposition (SVD). If the SVD fails, an error message is printed and the output matrices & vectors have no elements.

Examples:

mat A = rand<mat>(5,4);

mat coeff;
mat score;
colvec latent;
colvec tsquared;

princomp(coeff, score, latent, tsquared, A);

See also:

mat coeff = princomp_cov(mat X)
cx_mat coeff = princomp_cov(cx_mat X)

princomp_cov(mat coeff, colvec latent, mat X)
princomp_cov(cx_mat coeff, colvec latent, cx_mat X)

princomp_cov(mat coeff, colvec latent, colvec explained, mat X)
princomp_cov(cx_mat coeff, colvec latent, colvec explained, cx_mat X)

Principal component analysis of the covariance matrix of the input matrix X.

coeff: principal component coefficients.

latent: principal component variances.

explained: percentage of the total variance explained by each principal component.

The computation is based on singular value decomposition (SVD). If the SVD fails, an error message is printed and the output matrices & vectors have no elements.

Examples:

mat A = rand<mat>(5,4);

mat coeff;
colvec latent;
colvec explained; 
princomp_cov(coeff, latent, explained, A);

See also:

qr(Q,R,X)

Decomposition of matrix X into an orthogonal (Q) and a right triangular matrix (R), such that Q*R = X.

Examples:

mat X = rand<mat>(5,5);
mat Q, R;

qr(Q,R,X);

See also: QR decomposition in MathWorld

solve(colvec x, mat A, colvec b)
solve(mat X, mat A, mat B)

colvec x = solve(mat A, colvec b)
mat X = solve(mat A, mat B)

Solve a system of linear equations, i.e., A*X = B, where X is unknown.

For a square matrix A, this function is conceptually the same as X = inv(A)*B, but is done more efficiently.

Similar functionality to the "\" (left division operator) operator in Matlab/Octave, i.e., X = A \ B.

The number of rows in A and B must be the same.

This function will also try to provide approximate solutions to under-determined as well as over-determined systems (where matrix A is non-square).

If no solution is found:
- the two argument version will output a zero length vector/matrix
- the three argument versions output a bool set to false

NOTE: Old versions of the ATLAS library (e.g. 3.6) have a bug which can corrupt memory and crash your program -- later versions of ATLAS (e.g. 3.8) should work. Using the standard LAPACK library works without problems.

Examples:

mat A    = rand<mat>(5,5);
colvec b = rand<colvec>(5);
mat B    = rand<mat>(5,5);

colvec x = solve(A,b);
mat X    = solve(A,B);

colvec x2;
bool status = solve(x2, A, b);

See also:
- inv()
- linear system of equations in MathWorld

svd(colvec s, mat X),
svd(colvec s, cx_mat X)

colvec s = svd(mat X)
colvec s = svd(cx_mat X)

svd(mat U, colvec s, mat V, mat X)
svd(cx_mat U, colvec s, cx_mat V, cx_mat X)

The single and two argument versions compute the singular values of X

The four argument version computes the singular value decomposition of X, such that:
- X = U*diagmat(s)*trans(V), if X has real numbers
- X = U*diagmat(s)*htrans(V), if X has complex numbers

If the decomposition fails:
- the single argument version will output a zero length vector
- the two and four argument versions output a bool set to false

NOTE: Old versions of the ATLAS library (e.g. 3.6) have a bug which can corrupt memory and crash your program -- later versions of ATLAS (e.g. 3.8) should work. Using the standard LAPACK library works without problems.

Examples:

mat X = rand<mat>(5,5);

mat U;
colvec s;
mat V;
bool status = svd(U,s,V,X);

See also: singular value decomposition in MathWorld

Miscellaneous

math constants (pi, e, euler, gratio, sqrt2, eps, log_min, log_max)

Collection of constants, with their precision and/or value dependant on the numerical type and/or machine used.

The constants are stored as static functions in the Math<type> class, where type is either float or double.

For convenience, Math<float> has been typedefed as fmath, while Math<double> has been typedefed as math.

Meaning of the constants:

math::pi()		π, the ratio of any circle's circumference to its diameter
math::e()		base of the natural logarithm
math::euler()		Euler's constant, aka Euler-Mascheroni constant
math::gratio()		golden ratio
math::sqrt2()		square root of 2
math::eps()		the difference between 1 and the least value greater than 1 that is representable
math::log_min()		log of minimum non-zero value
math::log_max()		log of maximum value

Examples:

cout << "2.0 * pi = " << 2.0 * math::pi() << endl;

cout << "log_max for floats = ";
cout << fmath::log_max() << endl;

cout << "log_max for doubles = ";
cout << math::log_max() << endl;

physical constants (speed of light, etc)

Collection of fundamental physical constants, mainly taken from NIST and some from WolframAlpha on 2009-06-23.

Constants from NIST are in turn sourced from the 2006 CODATA values.

The constants are stored as static functions in the Phy<type> class, where type is either float or double.

For convenience, Phy<float> has been typedefed as fphy, while Phy<double> has been typedefed as phy.

Meaning of the constants:

phy::m_u()		atomic mass constant (in kg)
phy::N_A()		Avogadro constant
phy::k()		Boltzmann constant (in joules per kelvin)
phy::k_evk()		Boltzmann constant (in eV/K)
phy::a_0()		Bohr radius (in meters)
phy::mu_B()		Bohr magneton
phy::Z_0()		characteristic impedance of vacuum (in ohms)
phy::G_0()		conductance quantum (in siemens)
phy::k_e()		Coulomb's constant (in meters per farad)
phy::eps_0()		electric constant (in farads per meter)
phy::m_e()		electron mass (in kg)
phy::eV()		electron volt (in joules)
phy::e()		elementary charge (in coulombs)
phy::F()		Faraday constant (in coulombs)
phy::alpha()		fine-structure constant
phy::alpha_inv()		inverse fine-structure constant
phy::K_J()		Josephson constant
phy::mu_0()		magnetic constant (in henries per meter)
phy::phi_0()		magnetic flux quantum (in webers)
phy::R()		molar gas constant (in joules per mole kelvin)
phy::G()		Newtonian constant of gravitation (in newton square meters per kilogram squared)
phy::h()		Planck constant (in joule seconds)
phy::h_bar()		Planck constant over 2 pi, aka reduced Planck constant (in joule seconds)
phy::m_p()		proton mass (in kg)
phy::R_inf()		Rydberg constant (in reciprocal meters)
phy::c_0()		speed of light in vacuum (in meters per second)
phy::sigma()		Stefan-Boltzmann constant
phy::R_k()		von Klitzing constant (in ohms)
phy::b()		Wien wavelength displacement law constant

Examples:

cout << "speed of light = " << phy::c_0() << endl;

See also: physical constant entry in Wikipedia.

log_add(log_a, log_b)

Safe replacement for log(exp(log_a) + exp(log_b))

Usage:
- scalar_type log_c = log_add(log_a, log_b)
- scalar_type is either float or double
- log_a, log_b and log_c must have the same type

s32, u32

s32 is a typedef for a 32 bit wide signed int

u32 is a typedef for a 32 bit wide unsigned int

See also: C++ variable types

cx_float, cx_double

cx_float is a typedef for std::complex<float>

cx_double is a typedef for std::complex<double>

See also: complex numbers in the standard C++ library

conversion table between Matlab/Octave and Armadillo syntax

Matlab/Octave	Armadillo	Notes

A(1, 1)	A(0, 0)	indexing in Armadillo starts at 0
A(k, k)	A(k-1, k-1)

size(A,1)	A.n_rows	read only
size(A,2)	A.n_cols
size(Q,3)	Q.n_slices	Q is a cube (3D array)
numel(A)	A.n_elem

A(:,k)	A.col(k-1)
A(k, :)	A.row(k-1)
A(:, p:q)	A.cols(p-1, q-1)
A(p:q, :)	A.rows(p-1, q-1)
A(p:q, r:s)	A.submat(p-1, r-1, q-1, s-1)	A.submat(first_row, first_col, last_row, last_col)

Q(:, :, k)	Q.slice(k-1)	Q is a cube (3D array)
Q(:, :, t:u)	Q.slices(t-1, u-1)
Q(p:q, r:s, t:u)	Q.subcube(p-1, r-1, t-1, q-1, s-1, u-1)	.subcube(first_row, first_column, first_slice, last_row, last_column, last_slice)

A.'	trans(A)	simple transpose (the conjugate of each element is not taken)
A'	htrans(A)	Hermitian transpose (for complex matrices, the conjugate of each element is taken)

A = zeros(size(A))	A.zeros()
A = ones(size(A))	A.ones()
A = zeros(k)	A = zeros<mat>(k,k)
A = ones(k)	A = ones<mat>(k,k)

C = complex(A,B)	cx_mat C = cx_mat(A,B)

A .* B	A % B	element-wise multiplication
A ./ B	A / B	element-wise division
A \ B	solve(A,B)
A = A + 1;	A++
A = A - 1;	A--

A = [ 1 2; 3 4; ]	A = "1 2; 3 4;"
X = [ A B ]	X = join_rows(A,B)
X = [ A; B ]	X = join_cols(A,B)

A	cout << A << endl; or A.print("A =");

save -ascii 'A.dat' A	A.save("A.dat", raw_ascii);	Matlab/Octave matrices saved as ascii are readable by Armadillo (and vice-versa)
load -ascii 'A.dat'	A.load("A.dat", raw_ascii);

S = { 'abc'; 'def' }	field<std::string> S(2); S(0) = "abc"; S(1) = "def";	fields can store arbitrary objects

example program

If you save the program below as arma_prog.cpp, under Linux you can compile it using:
g++ arma_prog.cpp -o arma_prog -O1 -larmadillo

You may also want to have a look at the example programs that come with the Armadillo archive.

#include <iostream>
#include <armadillo>

using namespace std;
using namespace arma;

int main(int argc, char** argv)
  {
  mat A = rand<mat>(4,5);
  mat B = rand<mat>(4,5);
  
  cout << "A*trans(B) =" << endl;
  cout << A*trans(B) << endl;
  
  return 0;
  }

Reference for Armadillo 0.9.8