# Numerical Linear Algebra Magma has some limited capabilities with numerical linear algebra, that is, linear algebra over a real or complex field. The basis of these is the `RQDecomposition` (appropriate as Magma is row-based), and its companion `QLDecomposition`. Internally, Magma uses these functions to compute (e.g.) matrix inverses and determinants over **R** and **C**. All of the intrinsics starting with `Numeric` can also be accessed with this word omitted (by extended types on a real/complex matrix), though for some purposes it is useful to remind oneself that numerical linear algebra is not always the same. ## `RQDecomposition(M): Mtrx[RngReCom] -> Mtrx, AlgMatElt` Given a matrix $M$ over a real or complex field, write $M=RQ$ where $R$ is of triangular form and $Q$ orthogonal (or unitary in the complex case). In particular, $R$ has the same dimensions as $M$; when $M$ has more columns $c$ than rows $r$, then $R$ is of the form $[Z|U]$ where $Z$ is a zero matrix that has $(c-r)$ columns and $U$ is upper triangular. When $r>c$ then the bottom $c$ rows of $R$ form an upper triangular matrix. No normalisations are made except that $Q$ has determinant 1. ## `QLDecomposition(M): Mtrx[RngReCom] -> AlgMatElt, Mtrx` Given a matrix $M$ over a real or complex field, write $M=QL$ where $L$ is of triangular form and $Q$ orthogonal (or unitary in the complex case). This function works simply by applying `RQDecomposition` to the transpose. ## `NumericalInverse(M): Mtrx[RngReCom] -> AlgMatElt` ```magma Epsilon: FldReElt Default: ``` Given a square positive-dimension matrix over a real or complex field, compute its inverse. The `Epsilon` parameter allows one to specify what size of values should be considered to be zero (resulting in an inversion failure). By default, the value of `Epsilon` is computed in a way that is relative to the size of the entries and their precision. ## `Example: Numerlinalg Examples (ex-5a8543)` We give some examples of these decompositions. ```magma > RR := RealField(10); > r := 4; > c := 5; > A := [RR!Random([-2^10..2^10])/2^9 : i in [1..r*c]]; > M := Matrix(r,c,A); > M; [-0.1796875000 -1.394531250 -1.345703125 -0.4414062500 0.0000000000] [1.396484375 -1.496093750 -0.2617187500 -1.882812500 -1.679687500] [-1.169921875 -0.9062500000 1.707031250 -0.6582031250 1.480468750] [0.8867187500 -0.3359375000 0.5839843750 0.9707031250 1.896484375] > R,Q := RQDecomposition(M); R, Q; [0.0000000000 1.767478576 0.8447303787 -0.05853296288 0.3765439489] [0.0000000000 0.0000000000 2.924272555 -0.2048205303 1.424771353] [0.0000000000 0.0000000000 0.0000000000 2.589389831 -1.011949267] [0.0000000000 0.0000000000 0.0000000000 0.0000000000 -2.403971597] [-0.1278237197 -0.4893124575 0.3070157953 0.6619247514 -0.4602513884] [-0.3369945746 -0.5416048144 -0.7236111740 0.05012569151 0.2587917417] [0.6155218516 -0.6003864917 0.06838441636 -0.4759773097 -0.1715752984] [-0.5959647961 -0.2953736282 0.5643042371 -0.4119965068 0.2634387749] [-0.3688557516 0.1397427076 -0.2429248231 -0.4037914284 -0.7888963320] > Max([Abs(x) : x in Eltseq(R*Q-M)]); 3.492459655E-10 10 > Max([Abs(x) : x in Eltseq(Transpose(Q)*Q-1)]); 1.527951099E-10 14 > // > CC := ComplexField(5); > r := 5; > c := 3; > A := [CC!Random([-2^10..2^10])/2^9 : j in [1..r*c]]; > B := [CC!Random([-2^10..2^10])/2^9 : j in [1..r*c]]; > M := Matrix(r,c,A) + i*Matrix(r,c,B); > R,Q := RQDecomposition(M); R, Q; [-1.1471 + 1.0971*i -2.6950 - 0.38148*i 0.47294 - 0.95428*i] [-0.25971 - 0.0084722*i 0.61157 + 0.49147*i 1.5717 - 0.61096*i] [0.30902 + 0.59626*i 0.79478 + 1.0406*i -1.1157 - 1.6800*i] [0.00000 1.2132 + 1.4143*i -0.95447 + 0.14362*i] [0.00000 0.00000 3.0078 - 0.24966*i] [0.42713 - 0.0038340*i 0.51414 - 0.56570*i -0.39550 + 0.27712*i] [0.54343 + 0.60175*i -0.087580 + 0.0038354*i 0.56113 + 0.14178*i] [0.33511 + 0.21874*i -0.62188 + 0.14580*i -0.65721] > Max([Abs(x) : x in Eltseq(R*Q-M)]); 0.00020472 3 > Max([Abs(x) : x in Eltseq(ConjugateTranspose(Q)*Q-1)]); 0.00012207 9 > // > r := 5; > c := 5; > A := [RR!Random([-2^10..2^10])/2^9 : j in [1..r*c]]; > M := Matrix(r,c,A); > MI := NumericalInverse(M); MI; [0.5488038571 -0.3734925483 0.005386589048 0.4034292318 -0.6410457874] [0.3082898339 -0.5126616665 0.1612965318 -0.2971166681 0.07521900340] [-0.4446389269 0.4962564367 0.1565165202 -0.4722494193 0.2177560581] [-0.05815013416 -0.3372341073 0.03559373548 0.5402886617 -0.7996494168] [-0.9780538200 1.225635935 -0.4735623930 -1.060353816 1.356842221] > Max([Abs(x) : x in Eltseq(M*MI-1)]); 4.656612873E-10 4 ``` ## Rank, Kernel, Solution, and Pseudoinverse To compute the kernel (thus also rank) of a matrix $M$ or a solution to a system $vM=w$, the implementation first performs a QL decomposition internally. When $L$ is of a strict triangular form, namely $L=\bigl({Z\over A|T}\bigr)$ where $Z$ is zero and $T$ is not only lower triangular but also invertible, then the information can be read off immediately. In particular, the rank $r$ is the dimension of $T$, the kernel is the first $k$ rows of $Q$ where $k$ is number of rows in $Z$, while $vM=w$ requires $vQL=w$ so $v=w_rT^{-1}Q_r^T$ where $w_r$ and $Q_r$ are the rightmost $r$ columns of $w$ and $Q$; one checks $vQ_rA$ is (numerically) equal to the other columns of $w$. When $L$ is not of this form, then an arbitrary orthogonal/unitary transformation is performed on the input (on the right), and the algorithm is re-tried. This case occurs when the $r$ rightmost columns of $M$ are not linearly independent (Magma is internally working on the transpose of $M$ in the $RQ$ decomposition, so that the rightmost columns become the lowermost rows). Similarly, if a condition number for $T$ is too large, such a randomising transform can be profitably applied. Note that the kernel is returned as a matrix, and indeed has orthonormal rows. The same is true for the image. The computation of the `NumericalPseudoinverse` is similar, being based on the `QLDecomposition` and possibly also transforming the columns on the right, to ensure that $L$ is of the nice triangular form (this is internally called `QLCDecomposition`). The pseudoinverse of the row-independent (bottom) part $B$ of $L$ is then computed as the solution $\tilde P$ of $\tilde P BB^\star=B^\star$, and one row-extends $\tilde P$ to the left (by zero columns) to get the full (numerical) pseudoinverse $P$. ### `NumericalRank(M): Mtrx[RngReCom] -> RngIntElt` ```magma Epsilon: FldReElt Default: ``` Given a matrix over a real or complex field, compute its (numerical) rank. The `Epsilon` parameter allows one to specify what size of values should be considered to be zero. By default, the value of `Epsilon` is computed in a way that is relative to the size of the entries and their precision. ### `NumericalKernel(M): Mtrx[RngReCom] -> Mtrx` ```magma Epsilon: FldReElt Default: ``` Given a matrix over a real or complex field, compute its (numerical) kernel. The `Epsilon` parameter allows one to specify what size of values should be considered to be zero. By default, the value of `Epsilon` is computed in a way that is relative to the size of the entries and their precision. Exceptionally, `Kernel` returns a space, while `NumericalKernel` returns a basis matrix of it. ### `NumericalImage(M): Mtrx[RngReCom] -> Mtrx` ```magma Epsilon: FldReElt Default: ``` Given a matrix over a real or complex field, compute its (numerical) image, that is, its row space. The `Epsilon` parameter allows one to specify what size of values should be considered to be zero. By default, the value of `Epsilon` is computed in a way that is relative to the size of the entries and their precision. Exceptionally, `Image` returns a space, while `NumericalImage` returns a basis matrix of it. ### `NumericalSolution(M, w): Mtrx[RngReCom], Mtrx[RngReCom] -> Mtrx, Mtrx` ```magma Epsilon: FldReElt Default: ``` Given a matrix $M$ and a compatible matrix (or vector) $w$ over a real or complex field, solve numerically for $vM=w$. This gives an error if there is no solution. The second return value is the kernel of $M$. The `Epsilon` parameter allows one to specify what size of values should be considered to be zero. By default, the value of `Epsilon` is computed in a way that is relative to the size of the entries and their precision. Exceptionally, the second return value is a matrix for `NumericalSolution` and a space for `Solution`. ### `NumericalIsConsistent(M, w): Mtrx[RngReCom], Mtrx[RngReCom] -> BoolElt, Mtrx, Mtrx` ```magma Epsilon: FldReElt Default: ``` Given a matrix $M$ and a compatible matrix (or vector) $w$ over a real or complex field, determine whether $vM=w$ is numerically solvable. The second return value is a solution (when one exists)(, and the third return value is the (numerical) kernel of $M$. The `Epsilon` parameter allows one to specify what size of values should be considered to be zero. By default, the value of `Epsilon` is computed in a way that is relative to the size of the entries and their precision. Exceptionally, the third return value is a matrix for `NumericalIsConsistent` and a space for `IsConsistent`. ### `NumericalPseudoinverse(M): Mtrx[RngReCom] -> Mtrx` ```magma Epsilon: FldReElt Default: ``` Given a matrix $M$ over a real or complex field, compute its pseudoinverse $P$ such that $PMP=P$, $MPM=M$, $(MP)^\star=MP$ and $(PM)^\star=PM$. The `Epsilon` parameter allows one to specify what size of values should be considered to be zero. By default, the value of `Epsilon` is computed in a way that is relative to the size of the entries and their precision. Note that the Magma intrinsic `PseudoInverse` is quite different (largely applicable for matrices over Euclidean domains). ### `Example: Numerlinalg Examples2 (ex-b3a16a)` We give some examples for ranks, kernels, solutions, and pseudoinverses. ```magma > RR := RealField(10); > r := 3; > c := 5; > S := [-2^10..2^10]; > A := [Vector([RR!Random(S)/2^9 : i in [1..c]]) : j in [1..r]]; > M := Matrix(A cat [A[1]+A[2]-A[3]]); > M; [-1.347656250 -1.980468750 0.2070312500 0.2734375000 -1.273437500] [0.6464843750 1.726562500 -0.3769531250 0.2011718750 -1.646484375] [-0.1660156250 1.236328125 -1.500000000 1.119140625 -1.195312500] [-0.5351562500 -1.490234375 1.330078125 -0.6445312500 -1.724609375] > NumericalRank(M); 3 > NumericalKernel(M); // scaled to norm 1 [-0.5000000001 -0.4999999998 0.4999999999 0.5000000002] > v := Vector([RR!Random(S)/2^9 : i in [1..c]]); > NumericalIsConsistent(M,v); false > w := &+[M[i] * Random(S)/2^9 : i in [1..r]]; > NumericalIsConsistent(M,w); true (1.174316406 0.4165039062 0.7807617187 0.8100585938) > PI := NumericalPseudoinverse(M); PI; [-0.1991067080 0.1335911770 -0.1240987264 0.05858319545] [-0.1885126384 0.1881125824 0.02268111094 -0.02308116681] [-0.1331315985 0.1012317424 -0.2422989451 0.2103990890] [0.1423299946 -0.08451596130 0.1982404871 -0.1404264538] [-0.07615077666 -0.2393032832 -0.09708672579 -0.2183673341] > Max([Abs(x) : x in Eltseq(PI*M*PI-PI)]); 7.275957614E-11 9 > // > CC := ComplexField(30); > r := 5; > c := 4; > A := [CC!Random([-2^30..2^30])/2^28 : j in [1..r*c]]; > B := [CC!Random([-2^30..2^30])/2^28 : j in [1..r*c]]; > M := Matrix(r,c,A) + i*Matrix(r,c,B); > NumericalRank(M); 4 > NumericalKernel(M); [0.208965559769595845191678654168 + 0.0334373308308175861091994295635*i -0.0227045638063900267279802631268 - 0.354081499548224953341852663885*i -0.220008165304096296179634856382 - 0.815082824683314746916890536889*i 0.0146505331723745379341299824850 - 0.155641405430086914216739561265*i 0.209604845631814749141329584762 - 0.219520964338492435702429703559*i] > PI := NumericalPseudoinverse(M); > Max([Abs(x) : x in Eltseq(PI*M*PI-PI)]); 1.60521564392120817260990102530E-30 1 > Max([Abs(x) : x in Eltseq(M*PI*M-M)]); 2.54399997817916409183980874314E-29 18 > Max([Abs(x) : x in Eltseq(M*PI-ConjugateTranspose(M*PI))]); 7.88860905221011805411728565283E-30 25 > Max([Abs(x) : x in Eltseq(PI*M-ConjugateTranspose(PI*M))]); 5.03961120551682881427899748034E-30 8 ``` ## Eigenvalues and the Singular Value Decomposition Computing the eigenvalues of a real or complex square matrix is done by first computing the Hessenberg form (almost upper-triangular, but with possible nonzero entries on the immediate subdiagonal), and then applying iterative methods to achieve the Schur form (upper-triangular except for 2-by-2 blocks with complex eigenvalues in the real case). Special methods are used in the case where the matrix is symmetric. The computation of the singular value decomposition works by first computing the bidiagonal form, and then similarly iterating. ### `NumericalHessenbergForm(M): Mtrx[RngReCom] -> Mtrx, Mtrx` Given a square matrix $M$ over a real/complex field, compute a Hessenberg matrix $H$ that is upper-triangular except for possible nonzero entries on the immediate subdiagonal, and an orthogonal/unitary transformation matrix $Q$ such that $H=QMQ^{\star}$ and $QQ^{\star}=I$ (both equalities being numerically approximate). A procedure similar to the RQ decomposition is used, working on both rows and columns. ### `NumericalSchurForm(M): Mtrx[RngReCom] -> Mtrx, Mtrx` ```magma Transform: BoolElt Default: true ``` Given a square matrix $M$ over a real/complex field, compute the Schur form $S$ that is upper-triangular except for possible nonzero entries on the immediate subdiagonal corresponding to 2-by-2 real blocks with complex eigenvalues, and an orthogonal/unitary transformation matrix $Q$ such that $S=QMQ^{\star}$ and $QQ^{\star}=I$ (both equalities being numerically approximate). First the Hessenberg form is calculated, and then the unwanted subdiagonal terms are removed by a iterative process. ### `NumericalEigenvalues(M): Mtrx[RngReCom] -> SeqEnum` ```magma Balance: BoolElt Default: true ``` Given a square matrix $M$ over a real/complex field, compute numerical approximations to its eigenvalues. This operates by calculating the Schur form, and then returning the relevant eigenvalues. A form of Parlett-Reinsch balancing is used by default, to increase the precision of the answer. Exceptionally, `NumericalEigenvalues` returns a sequence of reals or complexes, while `Eigenvalues` on real/complex matrix returns a tuple of eigenvalues and multiplicities, for reasons of backwards compatibility (the multiplicities really cannot be guaranteed in the numerical case). Note that `NumericalEigenvalues` always returns a full sequence of eigenvalues for a real input, while `Eigenvalues` will return only those which are (approximately) real. ### `NumericalEigenvectors(M, e): Mtrx, FldComElt -> SeqEnum` Given a square matrix $M$ that is coercible into the complexes and an approximation $e$ to an eigenvalue of it, attempt to find eigenvectors. This function uses `NumericalKernel`, though returns the result as a sequence of (row) vectors. ### `NumericalBidiagonalForm(M): Mtrx[RngReCom] -> Mtrx, Mtrx, Mtrx` Given a matrix $M$ over a real/complex field, compute a bidiagonal form $B$ and orthogonal/unitary transformations $U,V$ such that $B=UMV^\star$. Here $B$ is upper bidiagonal if the number of columns is at least as large as the number of rows, and else is lower bidiagonal. ### `NumericalSingularValueDecomposition(M): Mtrx[RngReCom] -> Mtrx, Mtrx, Mtrx` Given a matrix $M$ over a real/complex field, compute the singular value decomposition $S$ and orthogonal/unitary transformations $U,V$ such that $S=UMV^\star$. Here $S$ is diagonal in its upper-left corner, with the diagonal entries being nonnegative and sorted in increasing size. ### `Example: Numerlinalg Examples3 (ex-604853)` We give some examples for eigenvalues. ```magma > RR := RealField(10); > r := 5; > c := 5; > A := [RR!Random([-2^10..2^10])/2^9 : i in [1..r*c]]; > M := Matrix(r,c,A); > H := NumericalHessenbergForm(M); H; [-1.387758248 0.2244522681 0.1050196798 -0.6292371533 1.500649611] [-1.737800256 -0.6560018402 -0.1292093616 -0.4045870847 -0.1501588456] [-0.0000000000 3.503611888 -0.4200132582 -1.521474700 -0.4771033714] [-0.0000000000 0.0000000000 1.241929396 -1.887789155 -1.398580785] [-0.0000000000 0.0000000000 0.0000000000 2.063240612 -1.882812500] > S := NumericalSchurForm(M); S; // need not be diag in real case [-0.5276818045 2.970316193 -0.8183434354 -1.227178569 1.207265137] [-0.6560286404 0.1374330601 0.7474073924 -0.7405697671 0.4496312767] [0.0000000000 0.0000000000 -2.733003988 0.5326130259 -0.2389927500] [0.0000000000 0.0000000000 0.0000000000 -1.878704188 2.041990652] [0.0000000000 0.0000000000 0.0000000000 -2.574033133 -1.232418078] > NumericalEigenvalues(M); [ -2.733003988, -1.555561133 - 2.269742311*$.1, -1.555561133 + 2.269742311*$.1, -0.1951243722 - 1.355735243*$.1, -0.1951243722 + 1.355735243*$.1 ] > M := M + Transpose(M); // make symmetric > NumericalEigenvalues(M); // all real [ -5.927297871, -4.230690093, -3.248785626, -1.661515789, 2.599539376 ] > // > // example with companion matrix of a polynomial > // > f:=Polynomial([Random(10^(12-j)) : j in [0..8]] cat [1]); f; x^9 + 8569*x^8 + 98753*x^7 + 533140*x^6 + 3583980*x^5 + 52066941*x^4 + 183889636*x^3 + 5743452355*x^2 + 22923098485*x + 725907276071 > NumericalEigenvalues(ChangeRing(CompanionMatrix(f),RR)); [ -8557.467295, -10.76902934 - 4.819656295*$.1, -10.76902934 + 4.819656295*$.1, -4.602845057 - 8.104590519*$.1, -4.602845057 + 8.104590519*$.1, 2.007741523 - 8.923992632*$.1, 2.007741523 + 8.923992632*$.1, 7.597780247 - 5.110282521*$.1, 7.597780247 + 5.110282521*$.1 ] > Sort([r[1] : r in Roots(f,ComplexField(10))],func); [ -8557.467295, -10.76902933 - 4.819656325*$.1, -10.76902933 + 4.819656325*$.1, -4.602845059 - 8.104590521*$.1, -4.602845059 + 8.104590521*$.1, 2.007741527 - 8.923992640*$.1, 2.007741527 + 8.923992640*$.1, 7.597780251 - 5.110282532*$.1, 7.597780251 + 5.110282532*$.1 ] ``` ### `Example: Numerlinalg Examples4 (ex-4e7c97)` We give some examples for bidiagonal form and singular value decomposition. ```magma > CC := ComplexField(5); > r := 4; > c := 5; > A := [CC!Random([-2^10..2^10])/2^9 : i in [1..r*c]]; > B := [CC!Random([-2^10..2^10])/2^9 : i in [1..r*c]]; > M := Matrix(r,c,A) + i*Matrix(r,c,B); > X := NumericalBidiagonalForm(M); X; > NumericalBidiagonalForm(M); // matrix over CC [2.1983 3.7789 0.00000 0.00000 0.00000] [0.00000 2.0496 1.8677 0.00000 0.00000] [0.00000 0.00000 2.6879 2.4315 0.00000] [0.00000 0.00000 0.00000 3.5035 0.00000] > NumericalSingularValueDecomposition(M); [4.9265 0.00000 0.00000 0.00000 0.00000] [0.00000 4.5508 0.00000 0.00000 0.00000] [0.00000 0.00000 2.5876 0.00000 0.00000] [0.00000 0.00000 0.00000 0.73134 0.00000] > [Sqrt(Real(x)) : x in NumericalEigenvalues(M*ConjugateTranspose(M))]; [ 0.73122, 2.5876, 4.5506, 4.9265 ] ```