## Matrices Module

This module implements algorithms for matrices on rational numbers. To this end, it first implements a type for rational and complex numbers, and based on these, algorithms for inverting a matrix, for solving linear (in)equation systems, and for computing eigenvalues (e.g. to solve linear differential equation systems) are provided.

### Types

 Type Description A general simplex problem is given by a set of linear ineqalities, i.e. lw[i] <= sum{j=0..maxCol} m[i,j] * x[j] <= up[i] for i=0..maxRows. The initial linear equations are used to define "basic variables" y[i] in terms of the nonbasic variables x[j]. Since the role of basic and nonbasic variables is changed during the simplex algorithm, we use maps basic and nonbasic to map indices to the current variable, hence, e.g., basic[i]=j means that variable j is currently the i-th basic variable. Variables are generally identified with integers that are mapped via idx to their names. The arrays lower and upper denote lower and upper bounds, if these exists for the variables, and eval denotes their current value in the solution process. This data structure represents a linear inequality system used for Fourier-Motzkin variable elimination. The linear inequality system is then simply an array of this record type that holds an optional lower and upper bound and a vector of coefficients for the linear expression sum_{i=0..L-1} vector[i] * x[i] for variables x[i]. Complex numbers cplx(re,im) are constructed by two floating point numbers which are the real and imaginary part of the complex number. Rational numbers rat(n,d) are constructed by two integers, the nominator n and the denominator d. We interprete numbers rat(0,d) as zero (where d!=0), and rat(n,0) as infinity.

### Functions and values

 Function or value Description  CompanionMatrixOfPolynomial p  Full Usage: CompanionMatrixOfPolynomial p Parameters: p : double[] Returns: float[,] Given a polynomial p(x) = sum{i=0..n} a[i]*x^i, compute its companion matrix, i.e., a matrix whose characteristic polynomial is p. One such matrix is the one that is zero except for the entries m[i+1,i]:=1 and m[i,n-1]:=-a[i]/a[n]. This way, computing the roots of the polynomial can be done by computing the eigenvalues of its matrix. p : double[] Returns: float[,]  DecomposeQR printSteps ostr dimRow dimCol m  Full Usage: DecomposeQR printSteps ostr dimRow dimCol m Parameters: printSteps : bool ostr : TextWriter dimRow : int dimCol : int m : double[,] Returns: float[,] * double[,] This function computes a QR-decomposition of a given matrix by means of Givens rotations, i.e., for a given matrix m, it computes an orthogonal matrix q (the column vectors of q are orthogonal unit length vectors) and an upper triangular matrix r such that m=q*r holds. Despite the real size of the matrix m, it will only consider the upper left submatrix of size dimRow x dimCol. Note that the inverse of an orthogonal matrix is its transpose. printSteps : bool ostr : TextWriter dimRow : int dimCol : int m : double[,] Returns: float[,] * double[,]  DecomposeQRC printSteps ostr dimRow dimCol m  Full Usage: DecomposeQRC printSteps ostr dimRow dimCol m Parameters: printSteps : bool ostr : TextWriter dimRow : int dimCol : int m : cplx[,] Returns: cplx[,] * cplx[,] This function computes a QR-decomposition of a given matrix by means of Givens rotations, i.e., for a given matrix m, it computes an orthogonal matrix q (the column vectors of q are orthogonal unit length vectors) and an upper triangular matrix r such that m=q*r holds. Despite the real size of the matrix m, it will only consider the upper left submatrix of size dimRow x dimCol. Note that the inverse of an orthogonal matrix is its transpose. See the comments on DecomposeQR to understand the implementation. printSteps : bool ostr : TextWriter dimRow : int dimCol : int m : cplx[,] Returns: cplx[,] * cplx[,]  EigenvaluesDoubleShift printSteps ostr m  Full Usage: EigenvaluesDoubleShift printSteps ostr m Parameters: printSteps : bool ostr : TextWriter m : double[,] Returns: cplx[] option This function computes an iteration for approximating the eigenvalues of a given matrix m over real numbers. To this end, we define the QR iteration with double shifts as follows: m[0] := m m[i+1] := trans(q[i])*m[i]*q[i] for the QR-decomposition q[i]*r[i] := m[i]*m[i] - (s1+s2)*m[i] + s1*s2*I using shifts s1,s2 Note that this is indeed a double shift with the sift values s1 and s2: To see this, define the following values as would be done with two linear shifts with values s1 and s2: b := r1*q1+s1*I with QR decomposition q1*r1=m[i]-s1*I m[i+1] := r2*q2+s2*I with QR decomposition q2*r2=b-s2*I It will compute the same m[i+1] as our definition above: We have m[i+1] := trans(q[i])*m[i]*q[i] with q[i]:=q1*q2: m[i+1] = r2*q2+s2*I = (trans(q2)*(b-s2*I))*q2+s2*I = trans(q2)*b*q2 - s2*trans(q2)*q2+s2*I = trans(q2)*b*q2 = trans(q2)*(r1*q1+s1*I)*q2 = trans(q2)*(trans(q1)*(m[i]-s1*I)*q1+s1*I)*q2 = trans(q2)*trans(q1)*m[i]*q1*q2 - s1*trans(q2)*(trans(q1)*q1*q2 + s1*trans(q2)*q2 = trans(q2)*trans(q1)*m[i]*q1*q2 = trans(q)*m[i]*q Second, we have q1*q2*r2*r1 = (m[i]-s2*I)*(m[i]-s1*I): q1*q2*r2*r1 = q1*(q2*r2)*r1 = q1*(b-s2*I)*r1 = (m[i]-s2*I)*q1*r1 = (m[i]-s2*I)*(m[i]-s1*I), where q1*(b-s2*I) = (m[i]-s2*I)*q1 can be seen as follows: q1*(b-s2*I) = (m[i]-s2*I)*q1 iff q1*b-s2*q1 = m[i]*q1-s2*q1 iff q1*b = m[i]*q1 iff q1*(r1*q1+s1*I) = m[i]*q1 iff (q1*r1)*q1+s1*q1 = m[i]*q1 iff (m[i]-s1*I)*q1+s1*q1 = m[i]*q1 iff m[i]*q1-s1*q1+s1*q1 = m[i]*q1 iff m[i]*q1 = m[i]*q1 Third, note that (m[i]-s2*I)*(m[i]-s1*I) is nothing else but the polynomial m[i]*m[i] - (s1+s2)*m[i] + s1*s2*I used in the first iteration. For complex conjugate numbers s1,s2, both s1+s2 and s1*s2 are real numbers, so that only real number computations are required. The above steps define therefore a double step QR iteration where for a real valued matrix only computations with real numbers are required even though some of the eigenvalues may be complex numbers. For the shift values, one considers the eigenvalues of the lower right 2x2 submatrix. printSteps : bool ostr : TextWriter m : double[,] Returns: cplx[] option  EigenvaluesLinearShiftC printSteps ostr m  Full Usage: EigenvaluesLinearShiftC printSteps ostr m Parameters: printSteps : bool ostr : TextWriter m : cplx[,] Returns: cplx[] option This function computes an iteration for approximating the eigenvalues of a given matrix m over complex numbers. To this end, we define the QR iteration with linear shifts as follows: m[0] := m m[i+1] := r[i]*q[i]+s[i]*I for the QR-decomposition q[i]*r[i] := m[i]-s[i]*I using some shift value s[i] (a complex number). Since trans(q[i]) is the inverse of q[i], the QR decomposition implies that r[i]=trans(q[i])*(m[i]-s[i]*I) holds, and thus m[i+1] = r[i]*q[i]+s[i]*I = trans(q[i])*(m[i]-s[i]*I)*q[i]+s[i]*I = trans(q[i])*m[i]*q[i]-s[i]*trans(q[i])*q[i]+s[i]*I = trans(q[i])*m[i]*q[i]. Since m[i+1] = trans(q[i])*m[i]*q[i] holds, it follows that m[i] and m[i+1] have the same eigenvalues: Note that m[i+1]*x = lambda*x implies trans(q[i])*m[i]*q[i]*x = lambda*x, thus, m[i]*q[i]*x = lambda*(q[i]*x), and m[i]*x = lambda*x implies trans(q[i])*m[i]*x = lambda*trans(q[i])*x, which in turn implies trans(q[i])*m[i]*q[i]*x' = lambda*x' for x' := trans(q[i])*x. Hence, if the sequence of matrices m[i] will converge, it will converge to a triangular matrix, where we can read its eigenvalues from its diagonal. However, convergence is poor unless we will make use of shifts. Good values for a shift are approximations of eigenvalues, in particular, we make use of the Wilkinson shift. printSteps : bool ostr : TextWriter m : cplx[,] Returns: cplx[] option  FourierMotzkin printSteps ostr vars liq  Full Usage: FourierMotzkin printSteps ostr vars liq Parameters: printSteps : bool ostr : TextWriter vars : string[] liq : LinearInequality[] Returns: rat[] option This function implements the Fourier-Motzkin variable elimination to determine a solution of a linear inequality system. printSteps : bool ostr : TextWriter vars : string[] liq : LinearInequality[] Returns: rat[] option  GaussJordan printSteps ostr m  Full Usage: GaussJordan printSteps ostr m Parameters: printSteps : bool ostr : TextWriter m : rat[,] GaussJordan transforms the given matrix into a diagonal matrix of using in-place operations that add multiples of another row to a row. It is the basis for solving linear equation systems or for inverting matrices. printSteps : bool ostr : TextWriter m : rat[,]  InverseMatrixRat printSteps ostr m  Full Usage: InverseMatrixRat printSteps ostr m Parameters: printSteps : bool ostr : TextWriter m : rat[,] Returns: rat[,] InverseMatrixRat inverts a given matrix in that it concatenates the unit matrix and applies then the Gauss-Jordan algorithm. The given matrix can be inverted if the left half is finally the unit matrix, and in that case the right half contains the inverse of the given matrix. printSteps : bool ostr : TextWriter m : rat[,] Returns: rat[,]  LinInequality2Simplex liq  Full Usage: LinInequality2Simplex liq Parameters: liq : LinearInequality array Returns: GeneralSimplexProblem Convert a linear inequality system to a GeneralSimplexProblem liq : LinearInequality array Returns: GeneralSimplexProblem  MatrixFloatToCmplx m  Full Usage: MatrixFloatToCmplx m Parameters: m : double[,] Returns: cplx[,] convert a matrix of floats to a matrix of complex numbers m : double[,] Returns: cplx[,]  MatrixTool qscoll  Full Usage: MatrixTool qscoll Parameters: qscoll : NameValueCollection This function is to be called by the html page. It expects in the qscoll for the key "MatrixTask" one of the following DecomposeQR, EigenvaluesShift EigenvaluesDoubleShift, GaussJordan, Inverse, LinearEquationSystem and will then call one of the functions of this module with the appropriate html output. qscoll : NameValueCollection  ParseGeneralSimplexProblem s  Full Usage: ParseGeneralSimplexProblem s Parameters: s : string Returns: GeneralSimplexProblem Parse a GeneralSimplexProblem from a string as the following let s = "-oo <= 5 1 <= 5 -oo <= 2 1 <= 4 -oo <= 1/2 1 <= 2 -oo <= 1/5 1 <= 1 1 <= 1/2 1 <= +oo 3/2 <= 1 1 <= +oo" s : string Returns: GeneralSimplexProblem  ParseLinearInequalitySystem s  Full Usage: ParseLinearInequalitySystem s Parameters: s : string Returns: LinearInequality[] Parse a string to a linear inequality system. s : string Returns: LinearInequality[]  ParseMatrixFloat s  Full Usage: ParseMatrixFloat s Parameters: s : string Returns: float[,] Parse a matrix from a string where rows are separated by newline symbols or semicolons, and entries within a row are separated by blanks or tabs. The entries of the string are expected as floating point numbers. s : string Returns: float[,]  ParseMatrixRat s  Full Usage: ParseMatrixRat s Parameters: s : string Returns: rat[,] Parse a matrix from a string where rows are separated by newline symbols or semicolons, and entries within a row are separated by blanks or tabs. The entries of the string are expected as rational numbers. s : string Returns: rat[,]  PolynomialRootsTool qscoll  Full Usage: PolynomialRootsTool qscoll Parameters: qscoll : NameValueCollection This function is to be called by the html page. It will parse a polynomial and will compute its roots via the eigenvalues of its companion matrix qscoll : NameValueCollection  PrintGeneralSimplexProblemASCII ostr gs  Full Usage: PrintGeneralSimplexProblemASCII ostr gs Parameters: ostr : TextWriter gs : GeneralSimplexProblem print a GeneralSimplex in ASCII ostr : TextWriter gs : GeneralSimplexProblem  PrintGeneralSimplexProblemLatex ostr gs  Full Usage: PrintGeneralSimplexProblemLatex ostr gs Parameters: ostr : TextWriter gs : GeneralSimplexProblem print a GeneralSimplex in Latex ostr : TextWriter gs : GeneralSimplexProblem  PrintGeneralSimplexProblemMathML ostr gs  Full Usage: PrintGeneralSimplexProblemMathML ostr gs Parameters: ostr : TextWriter gs : GeneralSimplexProblem print a GeneralSimplex in MathML ostr : TextWriter gs : GeneralSimplexProblem  PrintLinearInequalitySystem ostr vars liq  Full Usage: PrintLinearInequalitySystem ostr vars liq Parameters: ostr : TextWriter vars : string[] liq : LinearInequality[] Print a linear inequality system to output stream ostr. The additional array vars holds the names of the variables in the linear inequality system. ostr : TextWriter vars : string[] liq : LinearInequality[]  PrintLinearInequalitySystemLatex ostr vars liq  Full Usage: PrintLinearInequalitySystemLatex ostr vars liq Parameters: ostr : TextWriter vars : string[] liq : LinearInequality[] Print a linear inequality system to output stream ostr in Latex. The additional array vars holds the names of the variables in the linear inequality system. ostr : TextWriter vars : string[] liq : LinearInequality[]  PrintLinearInequalitySystemMathML ostr vars liq  Full Usage: PrintLinearInequalitySystemMathML ostr vars liq Parameters: ostr : TextWriter vars : string[] liq : LinearInequality[] Print a linear inequality system to output stream ostr in MathML. The additional array vars holds the names of the variables in the linear inequality system. ostr : TextWriter vars : string[] liq : LinearInequality[]  PrintMatrixCplx ostr m  Full Usage: PrintMatrixCplx ostr m Parameters: ostr : TextWriter m : cplx[,] print a matrix to output stream ostr ostr : TextWriter m : cplx[,]  PrintMatrixCplxToHtml ostr m  Full Usage: PrintMatrixCplxToHtml ostr m Parameters: ostr : TextWriter m : cplx[,] ostr : TextWriter m : cplx[,]  PrintMatrixCplxToLatex ostr indent m  Full Usage: PrintMatrixCplxToLatex ostr indent m Parameters: ostr : TextWriter indent : string m : cplx[,] ostr : TextWriter indent : string m : cplx[,]  PrintMatrixCplxToMathML ostr m  Full Usage: PrintMatrixCplxToMathML ostr m Parameters: ostr : TextWriter m : cplx[,] ostr : TextWriter m : cplx[,]  PrintMatrixFloat ostr m  Full Usage: PrintMatrixFloat ostr m Parameters: ostr : TextWriter m : double[,] print a matrix to output stream ostr ostr : TextWriter m : double[,]  PrintMatrixFloatToHtml ostr m  Full Usage: PrintMatrixFloatToHtml ostr m Parameters: ostr : TextWriter m : double[,] ostr : TextWriter m : double[,]  PrintMatrixFloatToLatex ostr indent m  Full Usage: PrintMatrixFloatToLatex ostr indent m Parameters: ostr : TextWriter indent : string m : double[,] ostr : TextWriter indent : string m : double[,]  PrintMatrixFloatToMathML ostr m  Full Usage: PrintMatrixFloatToMathML ostr m Parameters: ostr : TextWriter m : double[,] ostr : TextWriter m : double[,]  PrintMatrixRat ostr m  Full Usage: PrintMatrixRat ostr m Parameters: ostr : TextWriter m : rat[,] print a matrix to output stream ostr ostr : TextWriter m : rat[,]  PrintMatrixRatToHtml ostr m  Full Usage: PrintMatrixRatToHtml ostr m Parameters: ostr : TextWriter m : rat[,] ostr : TextWriter m : rat[,]  PrintMatrixRatToLatex ostr indent m  Full Usage: PrintMatrixRatToLatex ostr indent m Parameters: ostr : TextWriter indent : string m : rat[,] ostr : TextWriter indent : string m : rat[,]  PrintMatrixRatToMathML ostr m  Full Usage: PrintMatrixRatToMathML ostr m Parameters: ostr : TextWriter m : rat[,] ostr : TextWriter m : rat[,]  SimplexIntSolve optMaxCalls printSteps ostr gs  Full Usage: SimplexIntSolve optMaxCalls printSteps ostr gs Parameters: optMaxCalls : int option printSteps : bool ostr : TextWriter gs : GeneralSimplexProblem Returns: bool option SimplexIntSolve solves the integer linear programming problem given by gs. Essentially, it first calls SimplexRatSolve to determine a rational solution if there is one, and narrows then the constraints to integer bounds until an integer solution is found. It therefore solves an integer linear programming problem which is known to be NP-complete. Note however, that it is possible that this method will not terminate since there are inequation systems with unbounded rational solutions that have no integer solutions. For this reason, one can specify an optional maximal number of recursion steps, if that is reached, the function returns None, otherwise, it will find a solution and optMaxCalls : int option printSteps : bool ostr : TextWriter gs : GeneralSimplexProblem Returns: bool option  SimplexRatSolve printSteps ostr gs  Full Usage: SimplexRatSolve printSteps ostr gs Parameters: printSteps : bool ostr : TextWriter gs : GeneralSimplexProblem Returns: bool Simplex algorithm as decision procedure: Given a general Simplex problem (see the data type of this module, this function performs the required pivoting steps to determine a variable assignment in the contained eval component that satisfies the linear equation system of the tableau and also the constraints. It returns true if a solution has been found, and false if there is no solution on the rational numbers. The computations are done in-place in the given gs. Note that the Simplex algorithm may have an exponential runtime (which shows up rarely in practise), but it is known that the problem can be decided on the rational numbers in polynomial time (Kharmarkar's algorithm, and Kachian's ellipsoid method). printSteps : bool ostr : TextWriter gs : GeneralSimplexProblem Returns: bool  SimplexTool qscoll  Full Usage: SimplexTool qscoll Parameters: qscoll : NameValueCollection This function is to be called by a teaching tool's html page. It will parse a general simplex problem, will solve it and print the results to the html page. qscoll : NameValueCollection  SolveLinearEquationSystem printSteps ostr m  Full Usage: SolveLinearEquationSystem printSteps ostr m Parameters: printSteps : bool ostr : TextWriter m : rat[,] Returns: Set * Set * rat[,] Solves a linear equation system A*x=0 or A*x=b, where in the first case the given matrix is just A, and in the second case, the given matrix should be matrix (A|b). The result is a tuple (defVarSet,parVarSet,kern) where defVarSet is the set of variable indices that are defined in terms of the other ones contained in the set of parameters parVarSet, and kern is a matrix whose image is the kernel of the given matrix m. For the homogeneous equation system A*x=0, this means that all linear combinations of the column vectors of kern are solutions, and there is just the trivial solution x=0 iff the kernel is a dimCol x 0 matrix. For the inhomogeneous equation system A*x=b, there are no solutions if the variable associated with the rightmost column is a defined variable. Otherwise, the solutions of the inhomogeneous system are obtained from the solutions of the homogeneous system by setting the last parameter -1 and removing the last rows in the kernel's column vectors. printSteps : bool ostr : TextWriter m : rat[,] Returns: Set * Set * rat[,]  TensorProduct a b  Full Usage: TensorProduct a b Parameters: a : float[][] b : float[][] Returns: float[][] tensor product a : float[][] b : float[][] Returns: float[][]  bwdSubst uMat y  Full Usage: bwdSubst uMat y Parameters: uMat : rat[][] y : rat[] Returns: rat[] solving a linear equation system with an upper triangular matrix uMat and right hand side vector y uMat : rat[][] y : rat[] Returns: rat[]  fwdSubst lMat b  Full Usage: fwdSubst lMat b Parameters: lMat : rat[][] b : rat[] Returns: rat[] solving a linear equation system with a lower triangular matrix lMat and right hand side vector b lMat : rat[][] b : rat[] Returns: rat[]  luDecompose m  Full Usage: luDecompose m Parameters: m : rat[][] Returns: rat[][] * rat[][] Perform a LU-decomposition of the given matrix which is not modified by this function unlike in function luDecomposeInline. m : rat[][] Returns: rat[][] * rat[][]  luDecomposeInline m  Full Usage: luDecomposeInline m Parameters: m : rat[][] Perform a LU-decomposition of the given matrix where the matrices L and U are stored in the given matrix (L has an additional diagonal of 1s). m : rat[][]  lupDecompose m  Full Usage: lupDecompose m Parameters: m : rat[][] Returns: rat[][] * rat[][] * rat[][] Perform a LUP-decomposition of the given matrix which is not modified by this function unlike in function luDecomposeInline. The function returns the matrices L,U,P. m : rat[][] Returns: rat[][] * rat[][] * rat[][]  lupDecomposeInline m  Full Usage: lupDecomposeInline m Parameters: m : rat[][] Returns: int[] Perform a LUP-decomposition of the given matrix where the matrices L and U are stored in the given matrix (L has an additional diagonal of 1s) and the permutation matrix is returned as a permutation. m : rat[][] Returns: int[]  lupSolve P L U b  Full Usage: lupSolve P L U b Parameters: P : rat[][] L : rat[][] U : rat[][] b : rat[] Returns: rat[] solving a linear equation system A*x=b by a LUP decomposition P*A=L*U where the LUP-decomposition is given P : rat[][] L : rat[][] U : rat[][] b : rat[] Returns: rat[]  matMatProd a b  Full Usage: matMatProd a b Parameters: a : rat[][] b : rat[][] Returns: rat[][] matrix * matrix product a : rat[][] b : rat[][] Returns: rat[][]  matVecProd m x  Full Usage: matVecProd m x Parameters: m : rat[][] x : rat[] Returns: rat[] matrix * vector product m : rat[][] x : rat[] Returns: rat[]