diff options
Diffstat (limited to 'circuitpython/extmod/ulab/docs/manual/source/numpy-linalg.rst')
| -rw-r--r-- | circuitpython/extmod/ulab/docs/manual/source/numpy-linalg.rst | 386 |
1 files changed, 386 insertions, 0 deletions
diff --git a/circuitpython/extmod/ulab/docs/manual/source/numpy-linalg.rst b/circuitpython/extmod/ulab/docs/manual/source/numpy-linalg.rst new file mode 100644 index 0000000..8439f33 --- /dev/null +++ b/circuitpython/extmod/ulab/docs/manual/source/numpy-linalg.rst @@ -0,0 +1,386 @@ + +numpy.linalg +============ + +Functions in the ``linalg`` module can be called by prepending them by +``numpy.linalg.``. The module defines the following seven functions: + +1. `numpy.linalg.cholesky <#cholesky>`__ +2. `numpy.linalg.det <#det>`__ +3. `numpy.linalg.eig <#eig>`__ +4. `numpy.linalg.inv <#inv>`__ +5. `numpy.linalg.norm <#norm>`__ +6. `numpy.linalg.qr <#qr>`__ + +cholesky +-------- + +``numpy``: +https://docs.scipy.org/doc/numpy-1.17.0/reference/generated/numpy.linalg.cholesky.html + +The function of the Cholesky decomposition takes a positive definite, +symmetric square matrix as its single argument, and returns the *square +root matrix* in the lower triangular form. If the input argument does +not fulfill the positivity or symmetry condition, a ``ValueError`` is +raised. + +.. code:: + + # code to be run in micropython + + from ulab import numpy as np + + a = np.array([[25, 15, -5], [15, 18, 0], [-5, 0, 11]]) + print('a: ', a) + print('\n' + '='*20 + '\nCholesky decomposition\n', np.linalg.cholesky(a)) + +.. parsed-literal:: + + a: array([[25.0, 15.0, -5.0], + [15.0, 18.0, 0.0], + [-5.0, 0.0, 11.0]], dtype=float) + + ==================== + Cholesky decomposition + array([[5.0, 0.0, 0.0], + [3.0, 3.0, 0.0], + [-1.0, 1.0, 3.0]], dtype=float) + + + + +det +--- + +``numpy``: +https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.det.html + +The ``det`` function takes a square matrix as its single argument, and +calculates the determinant. The calculation is based on successive +elimination of the matrix elements, and the return value is a float, +even if the input array was of integer type. + +.. code:: + + # code to be run in micropython + + from ulab import numpy as np + + a = np.array([[1, 2], [3, 4]], dtype=np.uint8) + print(np.linalg.det(a)) + +.. parsed-literal:: + + -2.0 + + + +Benchmark +~~~~~~~~~ + +Since the routine for calculating the determinant is pretty much the +same as for finding the `inverse of a matrix <#inv>`__, the execution +times are similar: + +.. code:: + + # code to be run in micropython + + from ulab import numpy as np + + @timeit + def matrix_det(m): + return np.linalg.inv(m) + + m = np.array([[1, 2, 3, 4, 5, 6, 7, 8], [0, 5, 6, 4, 5, 6, 4, 5], + [0, 0, 9, 7, 8, 9, 7, 8], [0, 0, 0, 10, 11, 12, 11, 12], + [0, 0, 0, 0, 4, 6, 7, 8], [0, 0, 0, 0, 0, 5, 6, 7], + [0, 0, 0, 0, 0, 0, 7, 6], [0, 0, 0, 0, 0, 0, 0, 2]]) + + matrix_det(m) + +.. parsed-literal:: + + execution time: 294 us + + + +eig +--- + +``numpy``: +https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html + +The ``eig`` function calculates the eigenvalues and the eigenvectors of +a real, symmetric square matrix. If the matrix is not symmetric, a +``ValueError`` will be raised. The function takes a single argument, and +returns a tuple with the eigenvalues, and eigenvectors. With the help of +the eigenvectors, amongst other things, you can implement sophisticated +stabilisation routines for robots. + +.. code:: + + # code to be run in micropython + + from ulab import numpy as np + + a = np.array([[1, 2, 1, 4], [2, 5, 3, 5], [1, 3, 6, 1], [4, 5, 1, 7]], dtype=np.uint8) + x, y = np.linalg.eig(a) + print('eigenvectors of a:\n', y) + print('\neigenvalues of a:\n', x) + +.. parsed-literal:: + + eigenvectors of a: + array([[0.8151560042509081, -0.4499411232970823, -0.1644660242574522, 0.3256141906686505], + [0.2211334179893007, 0.7846992598235538, 0.08372081379922657, 0.5730077734355189], + [-0.1340114162071679, -0.3100776411558949, 0.8742786816656, 0.3486109343758527], + [-0.5183258053659028, -0.292663481927148, -0.4489749870391468, 0.6664142156731531]], dtype=float) + + eigenvalues of a: + array([-1.165288365404889, 0.8029365530314914, 5.585625756072663, 13.77672605630074], dtype=float) + + + + +The same matrix diagonalised with ``numpy`` yields: + +.. code:: + + # code to be run in CPython + + a = array([[1, 2, 1, 4], [2, 5, 3, 5], [1, 3, 6, 1], [4, 5, 1, 7]], dtype=np.uint8) + x, y = eig(a) + print('eigenvectors of a:\n', y) + print('\neigenvalues of a:\n', x) + +.. parsed-literal:: + + eigenvectors of a: + [[ 0.32561419 0.815156 0.44994112 -0.16446602] + [ 0.57300777 0.22113342 -0.78469926 0.08372081] + [ 0.34861093 -0.13401142 0.31007764 0.87427868] + [ 0.66641421 -0.51832581 0.29266348 -0.44897499]] + + eigenvalues of a: + [13.77672606 -1.16528837 0.80293655 5.58562576] + + +When comparing results, we should keep two things in mind: + +1. the eigenvalues and eigenvectors are not necessarily sorted in the + same way +2. an eigenvector can be multiplied by an arbitrary non-zero scalar, and + it is still an eigenvector with the same eigenvalue. This is why all + signs of the eigenvector belonging to 5.58, and 0.80 are flipped in + ``ulab`` with respect to ``numpy``. This difference, however, is of + absolutely no consequence. + +Computation expenses +~~~~~~~~~~~~~~~~~~~~ + +Since the function is based on `Givens +rotations <https://en.wikipedia.org/wiki/Givens_rotation>`__ and runs +till convergence is achieved, or till the maximum number of allowed +rotations is exhausted, there is no universal estimate for the time +required to find the eigenvalues. However, an order of magnitude can, at +least, be guessed based on the measurement below: + +.. code:: + + # code to be run in micropython + + from ulab import numpy as np + + @timeit + def matrix_eig(a): + return np.linalg.eig(a) + + a = np.array([[1, 2, 1, 4], [2, 5, 3, 5], [1, 3, 6, 1], [4, 5, 1, 7]], dtype=np.uint8) + + matrix_eig(a) + +.. parsed-literal:: + + execution time: 111 us + + + +inv +--- + +``numpy``: +https://docs.scipy.org/doc/numpy-1.17.0/reference/generated/numpy.linalg.inv.html + +A square matrix, provided that it is not singular, can be inverted by +calling the ``inv`` function that takes a single argument. The inversion +is based on successive elimination of elements in the lower left +triangle, and raises a ``ValueError`` exception, if the matrix turns out +to be singular (i.e., one of the diagonal entries is zero). + +.. code:: + + # code to be run in micropython + + from ulab import numpy as np + + m = np.array([[1, 2, 3, 4], [4, 5, 6, 4], [7, 8.6, 9, 4], [3, 4, 5, 6]]) + + print(np.linalg.inv(m)) + +.. parsed-literal:: + + array([[-2.166666666666667, 1.500000000000001, -0.8333333333333337, 1.0], + [1.666666666666667, -3.333333333333335, 1.666666666666668, -0.0], + [0.1666666666666666, 2.166666666666668, -0.8333333333333337, -1.0], + [-0.1666666666666667, -0.3333333333333333, 0.0, 0.5]], dtype=float64) + + + + +Computation expenses +~~~~~~~~~~~~~~~~~~~~ + +Note that the cost of inverting a matrix is approximately twice as many +floats (RAM), as the number of entries in the original matrix, and +approximately as many operations, as the number of entries. Here are a +couple of numbers: + +.. code:: + + # code to be run in micropython + + from ulab import numpy as np + + @timeit + def invert_matrix(m): + return np.linalg.inv(m) + + m = np.array([[1, 2,], [4, 5]]) + print('2 by 2 matrix:') + invert_matrix(m) + + m = np.array([[1, 2, 3, 4], [4, 5, 6, 4], [7, 8.6, 9, 4], [3, 4, 5, 6]]) + print('\n4 by 4 matrix:') + invert_matrix(m) + + m = np.array([[1, 2, 3, 4, 5, 6, 7, 8], [0, 5, 6, 4, 5, 6, 4, 5], + [0, 0, 9, 7, 8, 9, 7, 8], [0, 0, 0, 10, 11, 12, 11, 12], + [0, 0, 0, 0, 4, 6, 7, 8], [0, 0, 0, 0, 0, 5, 6, 7], + [0, 0, 0, 0, 0, 0, 7, 6], [0, 0, 0, 0, 0, 0, 0, 2]]) + print('\n8 by 8 matrix:') + invert_matrix(m) + +.. parsed-literal:: + + 2 by 2 matrix: + execution time: 65 us + + 4 by 4 matrix: + execution time: 105 us + + 8 by 8 matrix: + execution time: 299 us + + + +The above-mentioned scaling is not obeyed strictly. The reason for the +discrepancy is that the function call is still the same for all three +cases: the input must be inspected, the output array must be created, +and so on. + +norm +---- + +``numpy``: +https://numpy.org/doc/stable/reference/generated/numpy.linalg.norm.html + +The function takes a vector or matrix without options, and returns its +2-norm, i.e., the square root of the sum of the square of the elements. + +.. code:: + + # code to be run in micropython + + from ulab import numpy as np + + a = np.array([1, 2, 3, 4, 5]) + b = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + + print('norm of a:', np.linalg.norm(a)) + print('norm of b:', np.linalg.norm(b)) + +.. parsed-literal:: + + norm of a: 7.416198487095663 + norm of b: 16.88194301613414 + + + + +qr +-- + +``numpy``: +https://numpy.org/doc/stable/reference/generated/numpy.linalg.qr.html + +The function computes the QR decomposition of a matrix ``m`` of +dimensions ``(M, N)``, i.e., it returns two such matrices, ``q``\ ’, and +``r``, that ``m = qr``, where ``q`` is orthonormal, and ``r`` is upper +triangular. In addition to the input matrix, which is the first +positional argument, the function accepts the ``mode`` keyword argument +with a default value of ``reduced``. If ``mode`` is ``reduced``, ``q``, +and ``r`` are returned in the reduced representation. Otherwise, the +outputs will have dimensions ``(M, M)``, and ``(M, N)``, respectively. + +.. code:: + + # code to be run in micropython + + from ulab import numpy as np + + A = np.arange(6).reshape((3, 2)) + print('A: \n', A) + + print('complete decomposition') + q, r = np.linalg.qr(A, mode='complete') + print('q: \n', q) + print() + print('r: \n', r) + + print('\n\nreduced decomposition') + q, r = np.linalg.qr(A, mode='reduced') + print('q: \n', q) + print() + print('r: \n', r) + +.. parsed-literal:: + + A: + array([[0, 1], + [2, 3], + [4, 5]], dtype=int16) + complete decomposition + q: + array([[0.0, -0.9128709291752768, 0.408248290463863], + [-0.447213595499958, -0.3651483716701107, -0.8164965809277261], + [-0.8944271909999159, 0.1825741858350553, 0.408248290463863]], dtype=float64) + + r: + array([[-4.47213595499958, -5.813776741499454], + [0.0, -1.095445115010332], + [0.0, 0.0]], dtype=float64) + + + reduced decomposition + q: + array([[0.0, -0.9128709291752768], + [-0.447213595499958, -0.3651483716701107], + [-0.8944271909999159, 0.1825741858350553]], dtype=float64) + + r: + array([[-4.47213595499958, -5.813776741499454], + [0.0, -1.095445115010332]], dtype=float64) + + + |