4.13 Matrices

A matrix is a two-dimensional array of expressions. The elements of a matrix with m rows and n columns are accessed with two unsigned indices, the first one in the range 0…m-1, the second one in the range 0…n-1.

There are a couple of ways to construct matrices, with or without preset elements. The constructor

matrix::matrix(unsigned r, unsigned c);

creates a matrix with ‘r’ rows and ‘c’ columns with all elements set to zero.

The fastest way to create a matrix with preinitialized elements is to assign a list of comma-separated expressions to an empty matrix (see below for an example). But you can also specify the elements as a (flat) list with

matrix::matrix(unsigned r, unsigned c, const lst & l);

The function

ex lst_to_matrix(const lst & l);

constructs a matrix from a list of lists, each list representing a matrix row.

There is also a set of functions for creating some special types of matrices:

ex diag_matrix(const lst & l);
ex unit_matrix(unsigned x);
ex unit_matrix(unsigned r, unsigned c);
ex symbolic_matrix(unsigned r, unsigned c, const string & base_name);
ex symbolic_matrix(unsigned r, unsigned c, const string & base_name,
                   const string & tex_base_name);

diag_matrix() constructs a diagonal matrix given the list of diagonal elements. unit_matrix() creates an ‘x’ by ‘x’ (or ‘r’ by ‘c’) unit matrix. And finally, symbolic_matrix constructs a matrix filled with newly generated symbols made of the specified base name and the position of each element in the matrix.

Matrices often arise by omitting elements of another matrix. For instance, the submatrix S of a matrix M takes a rectangular block from M. The reduced matrix R is defined by removing one row and one column from a matrix M. (The determinant of a reduced matrix is called a Minor of M and can be used for computing the inverse using Cramer's rule.)

ex sub_matrix(const matrix&m, unsigned r, unsigned nr, unsigned c, unsigned nc);
ex reduced_matrix(const matrix& m, unsigned r, unsigned c);

The function sub_matrix() takes a row offset r and a column offset c and takes a block of nr rows and nc columns. The function reduced_matrix() has two integer arguments that specify which row and column to remove:

{
    matrix m(3,3);
    m = 11, 12, 13,
        21, 22, 23,
        31, 32, 33;
    cout << reduced_matrix(m, 1, 1) << endl;
    // -> [[11,13],[31,33]]
    cout << sub_matrix(m, 1, 2, 1, 2) << endl;
    // -> [[22,23],[32,33]]
}

Matrix elements can be accessed and set using the parenthesis (function call) operator:

const ex & matrix::operator()(unsigned r, unsigned c) const;
ex & matrix::operator()(unsigned r, unsigned c);

It is also possible to access the matrix elements in a linear fashion with the op() method. But C++-style subscripting with square brackets ‘[]’ is not available.

Here are a couple of examples for constructing matrices:

{
    symbol a("a"), b("b");

    matrix M(2, 2);
    M = a, 0,
        0, b;
    cout << M << endl;
     // -> [[a,0],[0,b]]

    matrix M2(2, 2);
    M2(0, 0) = a;
    M2(1, 1) = b;
    cout << M2 << endl;
     // -> [[a,0],[0,b]]

    cout << matrix(2, 2, lst(a, 0, 0, b)) << endl;
     // -> [[a,0],[0,b]]

    cout << lst_to_matrix(lst(lst(a, 0), lst(0, b))) << endl;
     // -> [[a,0],[0,b]]

    cout << diag_matrix(lst(a, b)) << endl;
     // -> [[a,0],[0,b]]

    cout << unit_matrix(3) << endl;
     // -> [[1,0,0],[0,1,0],[0,0,1]]

    cout << symbolic_matrix(2, 3, "x") << endl;
     // -> [[x00,x01,x02],[x10,x11,x12]]
}

The method matrix::is_zero_matrix() returns true only if all entries of the matrix are zeros. There is also method ex::is_zero_matrix() which returns true only if the expression is zero or a zero matrix.

There are three ways to do arithmetic with matrices. The first (and most direct one) is to use the methods provided by the matrix class:

matrix matrix::add(const matrix & other) const;
matrix matrix::sub(const matrix & other) const;
matrix matrix::mul(const matrix & other) const;
matrix matrix::mul_scalar(const ex & other) const;
matrix matrix::pow(const ex & expn) const;
matrix matrix::transpose() const;

All of these methods return the result as a new matrix object. Here is an example that calculates A*B-2*C for three matrices A, B and C:

{
    matrix A(2, 2), B(2, 2), C(2, 2);
    A =  1, 2,
         3, 4;
    B = -1, 0,
         2, 1;
    C =  8, 4,
         2, 1;

    matrix result = A.mul(B).sub(C.mul_scalar(2));
    cout << result << endl;
     // -> [[-13,-6],[1,2]]
    ...
}

The second (and probably the most natural) way is to construct an expression containing matrices with the usual arithmetic operators and pow(). For efficiency reasons, expressions with sums, products and powers of matrices are not automatically evaluated in GiNaC. You have to call the method

ex ex::evalm() const;

to obtain the result:

{
    ...
    ex e = A*B - 2*C;
    cout << e << endl;
     // -> [[1,2],[3,4]]*[[-1,0],[2,1]]-2*[[8,4],[2,1]]
    cout << e.evalm() << endl;
     // -> [[-13,-6],[1,2]]
    ...
}

The non-commutativity of the product A*B in this example is automatically recognized by GiNaC. There is no need to use a special operator here. See section Non-commutative objects, for more information about dealing with non-commutative expressions.

Finally, you can work with indexed matrices and call simplify_indexed() to perform the arithmetic:

{
    ...
    idx i(symbol("i"), 2), j(symbol("j"), 2), k(symbol("k"), 2);
    e = indexed(A, i, k) * indexed(B, k, j) - 2 * indexed(C, i, j);
    cout << e << endl;
     // -> -2*[[8,4],[2,1]].i.j+[[-1,0],[2,1]].k.j*[[1,2],[3,4]].i.k
    cout << e.simplify_indexed() << endl;
     // -> [[-13,-6],[1,2]].i.j
}

Using indices is most useful when working with rectangular matrices and one-dimensional vectors because you don't have to worry about having to transpose matrices before multiplying them. See section Indexed objects, for more information about using matrices with indices, and about indices in general.

The matrix class provides a couple of additional methods for computing determinants, traces, characteristic polynomials and ranks:

ex matrix::determinant(unsigned algo=determinant_algo::automatic) const;
ex matrix::trace() const;
ex matrix::charpoly(const ex & lambda) const;
unsigned matrix::rank() const;

The ‘algo’ argument of determinant() allows to select between different algorithms for calculating the determinant. The asymptotic speed (as parametrized by the matrix size) can greatly differ between those algorithms, depending on the nature of the matrix' entries. The possible values are defined in the ‘flags.h’ header file. By default, GiNaC uses a heuristic to automatically select an algorithm that is likely (but not guaranteed) to give the result most quickly.

Matrices may also be inverted using the ex matrix::inverse() method and linear systems may be solved with:

matrix matrix::solve(const matrix & vars, const matrix & rhs,
                     unsigned algo=solve_algo::automatic) const;

Assuming the matrix object this method is applied on is an m times n matrix, then vars must be a n times p matrix of symbolic indeterminates and rhs a m times p matrix. The returned matrix then has dimension n times p and in the case of an underdetermined system will still contain some of the indeterminates from vars. If the system is overdetermined, an exception is thrown.