# MatrixOps

© 2009,2012 John Abbott
GNU Free Documentation License, Version 1.2

CoCoALib Documentation Index

## User documentation for MatrixOps

`MatrixOps` gathers together a number of operations on matrices; in most cases these operations are happy to accept a `MatrixView` (see `MatrixView`) as argument.

When not specified, a matrix argument is of type `ConstMatrixView`.

There are two ways of multiplying two matrices together. The infix operators return a `DenseMatrix`; the procedural version may be slightly faster than the infix operator.

• `mul(matrix& lhs, M1, M2)` -- a procedure equivalent to `lhs = M1*M2;`, note that `lhs` might be a `SparseMatrix` (not yet implemented)
• `operator*(M1, M2)` -- the product `M1*M2`
• `operator+(M1, M2)` -- the sum `M1+M2`
• `operator-(M1, M2)` -- the difference `M1-M2`
• `power(M, n)` compute `n`-th power of `M`; if `n` is negative then `M` must be invertible
• `operator*(n, M1)` -- scalar multiple of `M1` by `n` (integer or RingElem)
• `operator*(M1, n)` -- scalar multiple of `M1` by `n` (integer or RingElem)
• `operator/(M1, n)` -- scalar multiple of `M1` by `1/n` (where `n` is integer or RingElem)
• `operator-(M1)` -- scalar multiple of `M1` by -1

Here are some matrix norms. The result is an element of the ring containing the matrix elements. Note that `FrobeniusNorm2` gives the square of the Frobenius norm (so that the value surely lies in the same ring).

• `FrobeniusNorm2(M)` -- the square of the Frobenius norm
• `OperatorNormInfinity(M)` -- the infinity norm, ring must be ordered
• `OperatorNorm1(M)` -- the one norm, ring must be ordered

Here are some fairly standard functions on matrices.

• `det(M)` -- determinant of `M` (M must be square)
• `rk(M)` -- rank of `M` (the base ring must be an integral domain)
• `inverse(M)` -- inverse of `M` as a `DenseMatrix`
• `adj(M)` -- classical adjoint of `M` as a `DenseMatrix`; sometimes called "adjugate".
• `PseudoInverse(M)` -- PseudoInverse of `M` as a `DenseMatrix`. I suspect that it requires that the matrix be of full rank.
• `LinSolve(M,rhs)` -- solve for `x` the linear system `M*x = rhs`; result is a `DenseMatrix`; if no soln exists, result is the 0-by-0 matrix
• `LinKer(M)` -- solve for `x` the linear system `M*x = 0`; returns a `DenseMatrix` whose columns are a base for `ker(M)`

Here are some standard operations where the method used is specified explicitly. It would usually be better to use the generic operations above, as those should automatically select the most appropriate method for the given matrix.

• `det2x2(M)` -- only for 2x2 matrices
• `det3x3(M)` -- only for 3x3 matrices
• `det4x4(M)` -- only for 4x4 matrices
• `det5x5(M)` -- only for 5x5 matrices
• `DetByGauss(M)` -- matrix must be over an integral domain
• `RankByGauss(std::vector<long>& IndepRows, M)`
• `InverseByGauss(M)` -- some restrictions (needs gcd)
• `AdjointByDetOfMinors(M)`
• `AdjointByInverse(M)` -- base ring must be integral domain
• `LinSolveByGauss(M,rhs)` -- solve a linear system using gaussian elimination (base ring must be a field), result is a `DenseMatrix`
• `LinSolveByHNF(M,rhs)` -- solve a linear system using Hermite NormalForm (base ring must be a PID), result is a `DenseMatrix`
• `LinSolveByModuleRepr(M,rhs)` -- solve a linear system using module element representation, result is a `DenseMatrix`
• `void GrammSchmidtRows(MatrixView& M)` -- NYI
• `void GrammSchmidtRows(MatrixView& M, long row)` -- NYI

## Maintainer documentation for MatrixOps

Most impls are quite straightforward.

`power` is slightly clever with its iterative impl of binary powering.

`LinSolveByGauss` is a little complicated because it tries to handle all cases (e.g. full rank or not, square or more rows than cols or more cols than rows)

## Bugs, Shortcomings and other ideas

Can we make a common "gaussian elimination" impl which is called by the various algorithms needing it, rather than having several separate implementations?

Is the procedure `mul` really any faster than the infix operator?

## Main changes

2012

• June: Added negation, multiplication and division of a matrix by a scalar.
• April: Added LinSolve family (incl. LinSolveByGauss, LinSolveByHNF, LinSolveByModuleRepr)

2011

• May: Added power fn for matrices: cannot yet handle negative powers.
• March: added multiplication by RingElem