Linear Algebra Library

Introduction

An overview of all the available functions can be obtained by using the MuPAD® function info. Here we see an extract of the functions available in the linear algebra library (we do not state the whole output generated by the call info(linalg), since the library contains more than 40 different functions):

info(linalg)
Library 'linalg': the linear algebra package

 -- Interface:
linalg::addCol,           linalg::addRow,
linalg::adjoint,          linalg::angle,
linalg::basis,            linalg::charmat,
linalg::charpoly,         linalg::col,
...

After being exported, library functions can also be used by their short notation. The function call use(linalg) exports all functions of linalg. After that one can use the function name gaussElim instead of linalg::gaussElim, for example.

Please note that user-defined procedures that use functions of the library linalg should always use the long notation linalg::functionname, in order to make sure that the unambiguity of the function name is guaranteed.

The easiest way to define a matrix A is using the command matrix. The following defines a 2×2 matrix:

A := matrix([[1, 2], [3, 2]])

Now, you can add or multiply matrices using the standard arithmetical operators of MuPAD:

A * A, 2 * A, 1/A

Further, you can use functions of the linalg library:

linalg::det(A)

The domain type returned by matrix is Dom::Matrix():

domtype(A)

which is introduced in the following section.

Data Types for Matrices and Vectors

The library linalg is based on the domains Dom::Matrix and Dom::SquareMatrix. These constructors enable the user to define matrices and they offer matrix arithmetic and several functions for matrix manipulation.

A domain created by Dom::Matrix represents matrices of arbitrary rows and columns over a specified ring. The domain constructor Dom::Matrix expects a coefficient ring of category Cat::Rng (a ring without unit) as argument. If no argument is given, the domain of matrices is created that represents matrices over the field of arithmetical expressions, i.e., the domain Dom::ExpressionField().

Be careful with calculations with matrices over this coefficient domain, because their entries usually do not have a unique representation (e.g., there is more than one representation of zero). You can normalize the components of such a matrix A with the command map(A, normal ).

The library Dom offers standard coefficient domains, such as the field of rational numbers (Dom::Rational), the ring of integers (Dom::Integer), the residue classes of integers (Dom::IntegerMod(n)) for an integer n, or even the rings of polynomials (such as Dom::DistributedPolynomial(ind,R) or Dom::Polynomial(R), where ind is the list of variables and R is the coefficient ring).

A domain created by the domain constructor Dom::SquareMatrix represents the ring of square matrices over a specified coefficient domain. Dom::SquareMatrix expects the number of rows of the square matrices and optionally a coefficient ring of category Cat::Rng.

There are several possibilities to define matrices of a domain created by Dom::Matrix or Dom::SquareMatrix. A matrix can be created by giving a two-dimensional array, a list of the matrix components, or a function that generates the matrix components:

delete a, b, c, d:
A := matrix([[a, b], [c, d]])

The command matrix actually is an abbreviation for the domain Dom::Matrix().

To create diagonal matrices one should use the option Diagonal (the third argument of matrix is either a function or a list):

B := matrix(2, 2, [2, -2], Diagonal)

The following two examples show the meaning of the third argument:

delete x: matrix(2, 2, () -> x), matrix(2, 2, x)

The MuPAD arithmetical operators are used to perform matrix arithmetic:

A * B - 2 * B

1/A

Next we create the 2×2 generalized Hilbert matrix (see also linalg::hilbert) as a matrix of the ring of two-dimensional square matrices:

MatQ2 := Dom::SquareMatrix(2, Dom::Rational)

H2 := MatQ2((i, j) -> 1/(i + j - 1))

A vector with n components is a 1×n matrix (a row vector) or a n×1 matrix (a column vector).

The components of a matrix or a vector are accessed using the index operator, i.e., A[i,j] returns the component of the row with index i and column with index j.

The input A[i, j]:= x sets the (i, j)-th component of the matrix A to the value of x.

The index operator can also be used to extract sub-matrices by giving ranges of integers as its arguments:

A := Dom::Matrix(Dom::Integer)(
  [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
)

A[1..3, 1..2], A[3..3, 1..3]

See also the function linalg::submatrix.

Remarks on Improving Runtime

The runtime of user-defined procedures that use functions of the linalg library and methods of the constructors Dom::Matrix and Dom::SquareMatrix can be considerably improved in certain cases.

  1. Some of the functions of the linalg library correspond to certain methods of the domain constructor Dom::Matrix in their name and functionality. These functions are implemented by calling relevant methods of the domain to that they belong, apart from additional argument checking. These functions enable an user-friendly usage on the interactive level after exporting.

    However, in user-defined procedures the methods of the corresponding domain should be used directly to avoid additionally calls of procedures.

    For example standard matrix manipulation functions such as deleting, extracting or swapping of rows and columns are defined as methods of the domain constructors Dom::Matrix and Dom::SquareMatrix.

    The method "gaussElim" offers a Gaussian elimination process for each domain created by these constructors.

  2. When creating a new matrix the method "new" is called. It converts each component of the matrix explicitly into a component the component ring, which may be time consuming.

    However, matrices and vectors are often the results of computations, whose components already are elements of the component ring. Thus, the conversion of the entries is not necessary. To take this into account, the domain constructors Dom::Matrix and Dom::SquareMatrix offer a method "create" to define matrices in the usual way but without the conversion of the components.

    Please note that this method does not test its arguments. Thus it should be used with caution.

  3. A further possibility of achieving better runtimes using functions of linalg or methods of the constructor Dom::Matrix is to store functions and methods that are called more than once in local variables. This enables a faster access of these functions and methods.

The following example shows how a user-defined procedure using functions of linalg and methods of the domain constructor Dom::Matrix may look like. It computes the adjoint of a square matrix defined over a commutative ring (see Cat::CommutativeRing):

adjoint := proc(A)
    local n, R, i, j, a, Ai, Mat,
          // local variables to store often used methods
          det, delRow, delCol, Rnegate;
begin
    if args(0) <> 1 then
        error("wrong number of arguments")
    end_if;

    Mat := A::dom;       // the domain of A
    R := Mat::coeffRing; // the component ring of A
    n := Mat::matdim(A); // the dimension of A; faster than calling
                         // 'linalg::matdim'!
    if testargs() then
        if Mat::hasProp(Cat::Matrix) <> TRUE then
            error("expecting a matrix")
        elif not R::hasProp( Cat::CommutativeRing ) then
            error("expecting matrix over a 'Cat::CommutativeRing'")
        elif n[1] <> n[2] then
            error("expecting a square matrix")
        end_if
    end_if;

    // store often used methods in local variables:
    det := linalg::det;
    delRow := Mat::delRow;  // faster than calling 'linalg::delRow'!
    delCol := Mat::delCol;  // faster than calling 'linalg::delCol'!
    Rnegate := R::_negate;  // faster than using the '-' operator!

    n := Mat::matdim(A)[1]; // faster than calling 'linalg::nrows'!
    a := array(1..n, 1..n);

    for i from 1 to n do
        Ai := delCol(A, i);
        for j from 1 to n do
            a[i, j] := det(delRow(Ai, j));
            if i + j mod 2 = 1 then
                a[i, j] := Rnegate(a[i, j])
            end_if
        end_for
    end_for;

    // create a new matrix: use Mat::create instead of Mat::new
    // because the entries of the array are already elements of R
    return(Mat::create(a))
end_proc:

We give an example:

MatZ6 := Dom::Matrix(Dom::IntegerMod(6)):
adjoint(MatZ6([[1, 5], [2, 4]]))

Was this topic helpful?