|
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.
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.
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.
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.
We give some examples of these decompositions.
> 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<i> := 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
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=((Z/A|T))
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=wrT - 1QrT where wr and Qr
are the rightmost r columns of w and Q; one checks vQrA 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 BBstar=Bstar,
and one row-extends tilde P to the left (by zero columns)
to get the full (numerical) pseudoinverse P.
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.
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.
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.
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.
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.
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).
We give some examples for ranks, kernels, solutions, and pseudoinverses.
> 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<i> := 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
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.
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=QMQstar and QQstar=I (both equalities
being numerically approximate). A procedure similar to the RQ decomposition
is used, working on both rows and columns.
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=QMQstar and QQstar=I (both equalities being
numerically approximate). First the Hessenberg form is calculated, and
then the unwanted subdiagonal terms are removed by a iterative process.
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.
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.
Given a matrix M over a real/complex field, compute a bidiagonal form B
and orthogonal/unitary transformations U, V such that B=UMVstar.
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.
Given a matrix M over a real/complex field, compute the singular
value decomposition S and orthogonal/unitary transformations U, V
such that S=UMVstar. Here S is diagonal in its upper-left corner,
with the diagonal entries being nonnegative and sorted in increasing size.
We give some examples for eigenvalues.
> 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<x,y|Real(x)-Real(y)>);
[ -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 ]
We give some examples for bidiagonal form and singular value decomposition.
> CC<i> := 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 ]
[Next][Prev] [Right] [Left] [Up] [Index] [Root]
|