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
.
Matrix accessors
M[i,j]
read the(i,j)
-entry of matrixM
SetEntry(M,i,j, val)
set the(i,j)
-entry of matrixM
GetRow(M,i)
return thei
-th row ofM
as avector<RingElem>
GetCol(M,j)
return thej
-th col ofM
as avector<RingElem>
GetRows(M)
return the rows ofM
as avector<vector<RingElem>>
GetCols(M)
return the cols ofM
as avector<vector<RingElem>>
FlattenByRows(M)
return entries ofM
in avector<RingElem>
in order 1st row, 2nd row, etcFlattenByCols(M)
return entries ofM
in avector<RingElem>
in order 1st col, 2nd col, etc Note thatGetRow
,GetCol
,GetRows
,GetCols
,FlattenByRows
andFlattenByCols
make copies of the matrix entries.
Matrix Arithmetic
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 tolhs = M1*M2;
, note thatlhs
might be aSparseMatrix
(not yet implemented)operator*(M1, M2)
-- the productM1*M2
operator+(M1, M2)
-- the sumM1+M2
operator-(M1, M2)
-- the differenceM1-M2
power(M, n)
computen
-th power ofM
; ifn
is negative thenM
must be invertibleoperator*(n, M1)
-- scalar multiple ofM1
byn
(integer or RingElem)operator*(M1, n)
-- scalar multiple ofM1
byn
(integer or RingElem)operator/(M1, n)
-- scalar multiple ofM1
by1/n
(wheren
is integer or RingElem)operator-(M1)
-- scalar multiple ofM1
by -1
Matrix norms
Here are some matrix norms. The result is an element of the ring
containing the matrix elements. Note that FrobeniusNormSq
gives the
square of the Frobenius norm (so that the value surely lies in the
same ring).
FrobeniusNormSq(M)
-- the square of the Frobenius normOperatorNormInfinity(M)
-- the infinity norm, ring must be orderedOperatorNorm1(M)
-- the one norm, ring must be ordered
Sundry functions
Here are some fairly standard functions on matrices.
det(M)
-- determinant ofM
(M must be square)IsZeroDet(M)
-- equivalent toIsZero(det(M))
(but may be faster)HadamardBoundSq(M)
-- computes row and column bounds in astruct
(fieldsmyRowBound
andmyColBound
)rk(M)
-- rank ofM
(the base ring must be an integral domain)inverse(M)
-- inverse ofM
as aDenseMatrix
adj(M)
-- classical adjoint ofM
as aDenseMatrix
; sometimes called "adjugate".rref(M)
-- compute a reduced row echelon form ofM
(orig. matrix is unchanged); matrix must be over a fieldPseudoInverse(M)
-- PseudoInverse ofM
as aDenseMatrix
. I suspect that it requires that the matrix be of full rank.LinSolve(M,rhs)
-- solve forx
the linear systemM*x = rhs
; result is aDenseMatrix
; if no soln exists, result is the 0-by-0 matrixLinKer(M)
-- solve forx
the linear systemM*x = 0
; returns aDenseMatrix
whose columns are a base forker(M)
LinKerZZ(M)
-- solve forx
the linear systemM*x = 0
; returns aDenseMatrix
whose columns are a ZZ-base for integer points inker(M)
Further sundry functions
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 matricesdet3x3(M)
-- only for 3x3 matricesdet4x4(M)
-- only for 4x4 matricesdet5x5(M)
-- only for 5x5 matricesDetByGauss(M)
-- matrix must be over an integral domainRankByGauss(std::vector<long>& IndepRows, M)
InverseByGauss(M)
-- some restrictions (needs gcd)AdjointByDetOfMinors(M)
AdjointByInverse(M)
-- base ring must be integral domainLinSolveByGauss(M,rhs)
-- solve a linear system using gaussian elimination (base ring must be a field), result is aDenseMatrix
LinSolveByHNF(M,rhs)
-- solve a linear system using Hermite NormalForm (base ring must be a PID), result is aDenseMatrix
LinSolveByModuleRepr(M,rhs)
-- solve a linear system using module element representation, result is aDenseMatrix
void GrammSchmidtRows(MatrixView& M)
-- NYIvoid 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