wehub-resource-sync
570 lines
14 KiB
ReStructuredText
570 lines
14 KiB
ReStructuredText
Matrices
|
|
========
|
|
|
|
Creating matrices
|
|
-----------------
|
|
|
|
Basic methods
|
|
.............
|
|
|
|
Matrices in mpmath are implemented using dictionaries. Only non-zero values are
|
|
stored, so it is cheap to represent sparse matrices.
|
|
|
|
The most basic way to create one is to use the ``matrix`` class directly. You
|
|
can create an empty matrix specifying the dimensions::
|
|
|
|
>>> from mpmath import (matrix, ones, zeros, randmatrix, nprint, chop, iv,
|
|
... lu_solve, residual, fp, lu, diag, eye, eps, qr)
|
|
>>> matrix(2)
|
|
matrix(
|
|
[['0.0', '0.0'],
|
|
['0.0', '0.0']])
|
|
>>> matrix(2, 3)
|
|
matrix(
|
|
[['0.0', '0.0', '0.0'],
|
|
['0.0', '0.0', '0.0']])
|
|
|
|
Calling ``matrix`` with one dimension will create a square matrix.
|
|
|
|
To access the dimensions of a matrix, use the ``rows`` or ``cols`` keyword::
|
|
|
|
>>> A = matrix(3, 2)
|
|
>>> A
|
|
matrix(
|
|
[['0.0', '0.0'],
|
|
['0.0', '0.0'],
|
|
['0.0', '0.0']])
|
|
>>> A.rows
|
|
3
|
|
>>> A.cols
|
|
2
|
|
|
|
You can also change the dimension of an existing matrix. This will set the
|
|
new elements to 0. If the new dimension is smaller than before, the
|
|
concerning elements are discarded::
|
|
|
|
>>> A.rows = 2
|
|
>>> A
|
|
matrix(
|
|
[['0.0', '0.0'],
|
|
['0.0', '0.0']])
|
|
|
|
Internally ``convert`` is applied every time an element is set. This is
|
|
done using the syntax A[row,column], counting from 0::
|
|
|
|
>>> A = matrix(2)
|
|
>>> A[1,1] = 1 + 1j
|
|
>>> print(A)
|
|
[0.0 0.0]
|
|
[0.0 (1.0 + 1.0j)]
|
|
|
|
A more comfortable way to create a matrix lets you use nested lists::
|
|
|
|
>>> matrix([[1, 2], [3, 4]])
|
|
matrix(
|
|
[['1.0', '2.0'],
|
|
['3.0', '4.0']])
|
|
|
|
Advanced methods
|
|
................
|
|
|
|
Convenient functions are available for creating various standard matrices::
|
|
|
|
>>> zeros(2)
|
|
matrix(
|
|
[['0.0', '0.0'],
|
|
['0.0', '0.0']])
|
|
>>> ones(2)
|
|
matrix(
|
|
[['1.0', '1.0'],
|
|
['1.0', '1.0']])
|
|
>>> diag([1, 2, 3]) # diagonal matrix
|
|
matrix(
|
|
[['1.0', '0.0', '0.0'],
|
|
['0.0', '2.0', '0.0'],
|
|
['0.0', '0.0', '3.0']])
|
|
>>> eye(2) # identity matrix
|
|
matrix(
|
|
[['1.0', '0.0'],
|
|
['0.0', '1.0']])
|
|
|
|
You can even create random matrices::
|
|
|
|
>>> randmatrix(2) # doctest:+SKIP
|
|
matrix(
|
|
[['0.53491598236191806', '0.57195669543302752'],
|
|
['0.85589992269513615', '0.82444367501382143']])
|
|
|
|
Vectors
|
|
.......
|
|
|
|
Vectors may also be represented by the ``matrix`` class (with rows = 1 or cols = 1).
|
|
For vectors there are some things which make life easier. A column vector can
|
|
be created using a flat list, a row vectors using an almost flat nested list::
|
|
|
|
>>> matrix([1, 2, 3])
|
|
matrix(
|
|
[['1.0'],
|
|
['2.0'],
|
|
['3.0']])
|
|
>>> matrix([[1, 2, 3]])
|
|
matrix(
|
|
[['1.0', '2.0', '3.0']])
|
|
|
|
Optionally vectors can be accessed like lists, using only a single index::
|
|
|
|
>>> x = matrix([1, 2, 3])
|
|
>>> x[1]
|
|
mpf('2.0')
|
|
>>> x[1,0]
|
|
mpf('2.0')
|
|
|
|
Other
|
|
.....
|
|
|
|
Like you probably expected, matrices can be printed::
|
|
|
|
>>> print(randmatrix(3)) # doctest:+SKIP
|
|
[ 0.782963853573023 0.802057689719883 0.427895717335467]
|
|
[0.0541876859348597 0.708243266653103 0.615134039977379]
|
|
[ 0.856151514955773 0.544759264818486 0.686210904770947]
|
|
|
|
Use ``nstr`` or ``nprint`` to specify the number of digits to print::
|
|
|
|
>>> nprint(randmatrix(5), 3) # doctest:+SKIP
|
|
[2.07e-1 1.66e-1 5.06e-1 1.89e-1 8.29e-1]
|
|
[6.62e-1 6.55e-1 4.47e-1 4.82e-1 2.06e-2]
|
|
[4.33e-1 7.75e-1 6.93e-2 2.86e-1 5.71e-1]
|
|
[1.01e-1 2.53e-1 6.13e-1 3.32e-1 2.59e-1]
|
|
[1.56e-1 7.27e-2 6.05e-1 6.67e-2 2.79e-1]
|
|
|
|
As matrices are mutable, you will need to copy them sometimes::
|
|
|
|
>>> A = matrix(2)
|
|
>>> A
|
|
matrix(
|
|
[['0.0', '0.0'],
|
|
['0.0', '0.0']])
|
|
>>> B = A.copy()
|
|
>>> B[0,0] = 1
|
|
>>> B
|
|
matrix(
|
|
[['1.0', '0.0'],
|
|
['0.0', '0.0']])
|
|
>>> A
|
|
matrix(
|
|
[['0.0', '0.0'],
|
|
['0.0', '0.0']])
|
|
|
|
Finally, it is possible to convert a matrix to a nested list. This is very useful,
|
|
as most Python libraries involving matrices or arrays (namely NumPy or SymPy)
|
|
support this format::
|
|
|
|
>>> B.tolist()
|
|
[[mpf('1.0'), mpf('0.0')], [mpf('0.0'), mpf('0.0')]]
|
|
|
|
|
|
Matrix operations
|
|
-----------------
|
|
|
|
You can add and subtract matrices of compatible dimensions::
|
|
|
|
>>> A = matrix([[1, 2], [3, 4]])
|
|
>>> B = matrix([[-2, 4], [5, 9]])
|
|
>>> A + B
|
|
matrix(
|
|
[['-1.0', '6.0'],
|
|
['8.0', '13.0']])
|
|
>>> A - B
|
|
matrix(
|
|
[['3.0', '-2.0'],
|
|
['-2.0', '-5.0']])
|
|
>>> A + ones(3)
|
|
Traceback (most recent call last):
|
|
File "<stdin>", line 1, in <module>
|
|
File "...", line 238, in __add__
|
|
raise ValueError('incompatible dimensions for addition')
|
|
ValueError: incompatible dimensions for addition
|
|
|
|
It is possible to multiply or add matrices and scalars. In the latter case the
|
|
operation will be done element-wise::
|
|
|
|
>>> A * 2
|
|
matrix(
|
|
[['2.0', '4.0'],
|
|
['6.0', '8.0']])
|
|
>>> A / 4
|
|
matrix(
|
|
[['0.25', '0.5'],
|
|
['0.75', '1.0']])
|
|
>>> A - 1
|
|
matrix(
|
|
[['0.0', '1.0'],
|
|
['2.0', '3.0']])
|
|
|
|
Of course you can perform matrix multiplication, if the dimensions are
|
|
compatible::
|
|
|
|
>>> A * B
|
|
matrix(
|
|
[['8.0', '22.0'],
|
|
['14.0', '48.0']])
|
|
>>> matrix([[1, 2, 3]]) * matrix([[-6], [7], [-2]])
|
|
matrix(
|
|
[['2.0']])
|
|
|
|
You can raise powers of square matrices::
|
|
|
|
>>> A**2
|
|
matrix(
|
|
[['7.0', '10.0'],
|
|
['15.0', '22.0']])
|
|
|
|
Negative powers will calculate the inverse::
|
|
|
|
>>> A**-1
|
|
matrix(
|
|
[['-2.0', '1.0'],
|
|
['1.5', '-0.5']])
|
|
>>> nprint(A * A**-1, 3)
|
|
[ 1.0 1.08e-19]
|
|
[-2.17e-19 1.0]
|
|
|
|
Matrix transposition is straightforward::
|
|
|
|
>>> A = ones(2, 3)
|
|
>>> A
|
|
matrix(
|
|
[['1.0', '1.0', '1.0'],
|
|
['1.0', '1.0', '1.0']])
|
|
>>> A.T
|
|
matrix(
|
|
[['1.0', '1.0'],
|
|
['1.0', '1.0'],
|
|
['1.0', '1.0']])
|
|
|
|
|
|
Norms
|
|
.....
|
|
|
|
Sometimes you need to know how "large" a matrix or vector is. Due to their
|
|
multidimensional nature it's not possible to compare them, but there are
|
|
several functions to map a matrix or a vector to a positive real number, the
|
|
so called norms.
|
|
|
|
.. autofunction :: mpmath.norm
|
|
|
|
.. autofunction :: mpmath.mnorm
|
|
|
|
|
|
Linear algebra
|
|
--------------
|
|
|
|
Determinant and Rank
|
|
....................
|
|
|
|
.. autofunction :: mpmath.det
|
|
|
|
.. autofunction :: mpmath.rank
|
|
|
|
|
|
Decompositions
|
|
..............
|
|
|
|
.. autofunction :: mpmath.cholesky
|
|
|
|
|
|
Linear equations
|
|
................
|
|
|
|
Basic linear algebra is implemented; you can for example solve the linear
|
|
equation system::
|
|
|
|
x + 2*y = -10
|
|
3*x + 4*y = 10
|
|
|
|
using ``lu_solve``::
|
|
|
|
>>> A = matrix([[1, 2], [3, 4]])
|
|
>>> b = matrix([-10, 10])
|
|
>>> x = lu_solve(A, b)
|
|
>>> x
|
|
matrix(
|
|
[['30.0'],
|
|
['-20.0']])
|
|
|
|
If you don't trust the result, use ``residual`` to calculate
|
|
the residual `||A x-b||`::
|
|
|
|
>>> residual(A, x, b)
|
|
matrix(
|
|
[['3.46944695195361e-18'],
|
|
['3.46944695195361e-18']])
|
|
>>> str(eps)
|
|
'2.22044604925031e-16'
|
|
|
|
As you can see, the solution is quite accurate. The error is caused by the
|
|
inaccuracy of the internal floating-point arithmetic. Though, it's even smaller
|
|
than the current machine epsilon, which basically means you can trust the
|
|
result.
|
|
|
|
If you need more speed, use NumPy, or use ``fp`` instead ``mp`` matrices
|
|
and methods::
|
|
|
|
>>> A = fp.matrix([[1, 2], [3, 4]])
|
|
>>> b = fp.matrix([-10, 10])
|
|
>>> fp.lu_solve(A, b)
|
|
matrix(
|
|
[['29.999999999999996'],
|
|
['-19.999999999999996']])
|
|
|
|
``lu_solve`` accepts overdetermined systems. It is usually not possible to solve
|
|
such systems, so the residual is minimized instead. Internally this is done
|
|
using Cholesky decomposition to compute a least squares approximation. This means
|
|
that that ``lu_solve`` will square the errors. If you can't afford this, use
|
|
``qr_solve`` instead. It is twice as slow but more accurate, and it calculates
|
|
the residual automatically.
|
|
|
|
.. autofunction:: mpmath.lu_solve
|
|
|
|
Matrix factorization
|
|
....................
|
|
|
|
The function ``lu`` computes an explicit LU factorization of a matrix::
|
|
|
|
>>> P, L, U = lu(matrix([[0,2,3],[4,5,6],[7,8,9]]))
|
|
>>> print(P)
|
|
[0.0 0.0 1.0]
|
|
[1.0 0.0 0.0]
|
|
[0.0 1.0 0.0]
|
|
>>> print(L)
|
|
[ 1.0 0.0 0.0]
|
|
[ 0.0 1.0 0.0]
|
|
[0.571428571428571 0.214285714285714 1.0]
|
|
>>> print(U)
|
|
[7.0 8.0 9.0]
|
|
[0.0 2.0 3.0]
|
|
[0.0 0.0 0.214285714285714]
|
|
>>> print(P.T*L*U)
|
|
[0.0 2.0 3.0]
|
|
[4.0 5.0 6.0]
|
|
[7.0 8.0 9.0]
|
|
|
|
The function ``qr`` computes a QR factorization of a matrix::
|
|
|
|
>>> A = matrix([[1, 2], [3, 4], [1, 1]])
|
|
>>> Q, R = qr(A)
|
|
>>> print(Q)
|
|
[-0.301511344577764 0.861640436855329 0.408248290463863]
|
|
[-0.904534033733291 -0.123091490979333 -0.408248290463863]
|
|
[-0.301511344577764 -0.492365963917331 0.816496580927726]
|
|
>>> print(R)
|
|
[-3.3166247903554 -4.52267016866645]
|
|
[ 0.0 0.738548945875996]
|
|
[ 0.0 0.0]
|
|
>>> print(Q * R)
|
|
[1.0 2.0]
|
|
[3.0 4.0]
|
|
[1.0 1.0]
|
|
>>> print(chop(Q.T * Q))
|
|
[1.0 0.0 0.0]
|
|
[0.0 1.0 0.0]
|
|
[0.0 0.0 1.0]
|
|
|
|
|
|
The singular value decomposition
|
|
................................
|
|
|
|
The routines ``svd_r`` and ``svd_c`` compute the singular value decomposition
|
|
of a real or complex matrix A. ``svd`` is an unified interface calling
|
|
either ``svd_r`` or ``svd_c`` depending on whether *A* is real or complex.
|
|
|
|
Given *A*, two orthogonal (*A* real) or unitary (*A* complex) matrices *U* and *V*
|
|
are calculated such that
|
|
|
|
.. math ::
|
|
|
|
A = U S V, \quad U' U = 1, \quad V V' = 1
|
|
|
|
where *S* is a suitable shaped matrix whose off-diagonal elements are zero.
|
|
Here ' denotes the hermitian transpose (i.e. transposition and complex
|
|
conjugation). The diagonal elements of *S* are the singular values of *A*,
|
|
i.e. the square roots of the eigenvalues of `A' A` or `A A'`.
|
|
|
|
Examples::
|
|
|
|
>>> from mpmath import mp
|
|
>>> A = mp.matrix([[2, -2, -1], [3, 4, -2], [-2, -2, 0]])
|
|
>>> S = mp.svd_r(A, compute_uv = False)
|
|
>>> print(S)
|
|
[6.0]
|
|
[3.0]
|
|
[1.0]
|
|
>>> U, S, V = mp.svd_r(A)
|
|
>>> print(mp.chop(A - U * mp.diag(S) * V))
|
|
[0.0 0.0 0.0]
|
|
[0.0 0.0 0.0]
|
|
[0.0 0.0 0.0]
|
|
|
|
|
|
The Schur decomposition
|
|
.......................
|
|
|
|
This routine computes the Schur decomposition of a square matrix *A*.
|
|
Given *A*, a unitary matrix *Q* is determined such that
|
|
|
|
.. math ::
|
|
|
|
Q' A Q = R, \quad Q' Q = Q Q' = 1
|
|
|
|
where *R* is an upper right triangular matrix. Here ' denotes the
|
|
hermitian transpose (i.e. transposition and conjugation).
|
|
|
|
Examples::
|
|
|
|
>>> from mpmath import mp
|
|
>>> A = mp.matrix([[3, -1, 2], [2, 5, -5], [-2, -3, 7]])
|
|
>>> Q, R = mp.schur(A)
|
|
>>> mp.nprint(R, 3)
|
|
[2.0 0.417 2.53]
|
|
[0.0 4.0 4.74]
|
|
[0.0 0.0 9.0]
|
|
>>> print(mp.chop(A - Q * R * Q.transpose_conj()))
|
|
[0.0 0.0 0.0]
|
|
[0.0 0.0 0.0]
|
|
[0.0 0.0 0.0]
|
|
|
|
|
|
The eigenvalue problem
|
|
......................
|
|
|
|
The routine ``eig`` solves the (ordinary) eigenvalue problem for a real or complex
|
|
square matrix *A*. Given *A*, a vector *E* and matrices *ER* and *EL* are calculated such that
|
|
|
|
.. code ::
|
|
|
|
A ER[:,i] = E[i] ER[:,i]
|
|
EL[i,:] A = EL[i,:] E[i]
|
|
|
|
*E* contains the eigenvalues of *A*. The columns of *ER* contain the right eigenvectors
|
|
of *A* whereas the rows of *EL* contain the left eigenvectors.
|
|
|
|
|
|
Examples::
|
|
|
|
>>> from mpmath import mp
|
|
>>> A = mp.matrix([[3, -1, 2], [2, 5, -5], [-2, -3, 7]])
|
|
>>> E, ER = mp.eig(A)
|
|
>>> print(mp.chop(A * ER[:,0] - E[0] * ER[:,0]))
|
|
[0.0]
|
|
[0.0]
|
|
[0.0]
|
|
>>> E, EL, ER = mp.eig(A,left = True, right = True)
|
|
>>> E, EL, ER = mp.eig_sort(E, EL, ER)
|
|
>>> mp.nprint(E)
|
|
[2.0, 4.0, 9.0]
|
|
>>> print(mp.chop(A * ER[:,0] - E[0] * ER[:,0]))
|
|
[0.0]
|
|
[0.0]
|
|
[0.0]
|
|
>>> print(mp.chop( EL[0,:] * A - EL[0,:] * E[0]))
|
|
[0.0 0.0 0.0]
|
|
|
|
|
|
See also [Stoer]_ and [Kresser]_.
|
|
|
|
|
|
The symmetric eigenvalue problem
|
|
................................
|
|
|
|
The routines ``eigsy`` and ``eighe`` solve the (ordinary) eigenvalue problem
|
|
for a real symmetric or complex hermitian square matrix *A*.
|
|
``eigh`` is an unified interface for this two functions calling either
|
|
``eigsy`` or ``eighe`` depending on whether *A* is real or complex.
|
|
|
|
Given *A*, an orthogonal (*A* real) or unitary matrix *Q* (*A* complex) is
|
|
calculated which diagonalizes A:
|
|
|
|
.. math ::
|
|
|
|
Q' A Q = \operatorname{diag}(E), \quad Q Q' = Q' Q = 1
|
|
|
|
Here diag(*E*) a is diagonal matrix whose diagonal is *E*.
|
|
' denotes the hermitian transpose (i.e. ordinary transposition and
|
|
complex conjugation).
|
|
|
|
The columns of *Q* are the eigenvectors of *A* and *E* contains the eigenvalues:
|
|
|
|
.. code ::
|
|
|
|
A Q[:,i] = E[i] Q[:,i]
|
|
|
|
Examples::
|
|
|
|
>>> from mpmath import mp
|
|
>>> A = mp.matrix([[3, 2], [2, 0]])
|
|
>>> E = mp.eigsy(A, eigvals_only = True)
|
|
>>> print(E)
|
|
[-1.0]
|
|
[ 4.0]
|
|
>>> A = mp.matrix([[1, 2], [2, 3]])
|
|
>>> E, Q = mp.eigsy(A) # alternative: E, Q = mp.eigh(A)
|
|
>>> print(mp.chop(A * Q[:,0] - E[0] * Q[:,0]))
|
|
[0.0]
|
|
[0.0]
|
|
>>> A = mp.matrix([[1, 2 + 5j], [2 - 5j, 3]])
|
|
>>> E, Q = mp.eighe(A) # alternative: E, Q = mp.eigh(A)
|
|
>>> print(mp.chop(A * Q[:,0] - E[0] * Q[:,0]))
|
|
[0.0]
|
|
[0.0]
|
|
|
|
|
|
See also [Golub]_, [GolubWelsch]_, [Stoer]_ and [Stroud]_.
|
|
|
|
|
|
Determinant
|
|
...........
|
|
|
|
The determinant of a square matrix is computed by the
|
|
function ``det``::
|
|
|
|
>>> from mpmath import mp
|
|
>>> A = mp.matrix([[7, 2], [1.5, 3]])
|
|
>>> print(mp.det(A))
|
|
18.0
|
|
|
|
|
|
Interval and double-precision matrices
|
|
--------------------------------------
|
|
|
|
The ``iv.matrix`` and ``fp.matrix`` classes convert inputs
|
|
to intervals and Python floating-point numbers respectively.
|
|
|
|
Interval matrices can be used to perform linear algebra operations
|
|
with rigorous error tracking::
|
|
|
|
>>> a = iv.matrix([['0.1','0.3','1.0'],
|
|
... ['7.1','5.5','4.8'],
|
|
... ['3.2','4.4','5.6']])
|
|
>>>
|
|
>>> b = iv.matrix(['4','0.6','0.5'])
|
|
>>> c = iv.lu_solve(a, b)
|
|
>>> print(c)
|
|
[ [5.2582327113062393041, 5.2582327113062749951]]
|
|
[[-13.155049396267856583, -13.155049396267821167]]
|
|
[ [7.4206915477497212555, 7.4206915477497310922]]
|
|
>>> print(a*c)
|
|
[ [3.9999999999999866773, 4.0000000000000133227]]
|
|
[[0.59999999999972430942, 0.60000000000027142733]]
|
|
[[0.49999999999982236432, 0.50000000000018474111]]
|
|
|
|
Matrix functions
|
|
----------------
|
|
|
|
.. autofunction :: mpmath.expm
|
|
.. autofunction :: mpmath.cosm
|
|
.. autofunction :: mpmath.sinm
|
|
.. autofunction :: mpmath.sqrtm
|
|
.. autofunction :: mpmath.logm
|
|
.. autofunction :: mpmath.powm
|